Skip to content

GitLab

Commit 20015aff authored by Gwenaelle Gremion's avatar Gwenaelle Gremion
Browse files

Code commenté

parent c946efdb
Expand all Hide whitespace changes
Inline Side-by-side
nml/bio.nml
View file @ 20015aff
...@@ -10,20 +10,22 @@ ...@@ -10,20 +10,22 @@
! 4: Fasham et al. 1990 (7 variables) ! 4: Fasham et al. 1990 (7 variables)
! 5: IOW-ERGOM MaBenE version (9 variables) ! 5: IOW-ERGOM MaBenE version (9 variables)
! 6: ISMER model (9 variables) ! 6: ISMER model (9 variables)
! 7: GSJ model (11 variables) ! 7: The GSJ model, modified from ISMER
! 8: Ariadna Nocera's model (tbd) ! 8: The model of Ariadna Celina Nocera, modified from Fasham
! 9: The model with four detrital compartments
! 10: The model for the north Water polynya - NPZ4D (14 Variables)
! !
! bio_eulerian -> state variables are Eulerian (.true./.false.) ! bio_eulerian -> state variables are Eulerian (.true./.false.)
! !
! cnpar -> Cranck-Nicolson parameter for vertical diffusion ! cnpar -> Cranck-Nicolson parameter for vertical diffusion
! !
! w_adv_discr -> advection scheme for vertical motion ! w_adv_discr -> advection scheme for vertical motion - case as in adv_center.F90 (method)
! 1: first order upstream ! 1: first order upstream - case (UPSTREAM)
! 2: not coded yet ! 2: not coded yet -
! 3: third-order polynomial ! 3: third-order polynomial - case (P1)
! 4: TVD with Superbee limiter ! 4: TVD with Superbee limiter - case (Superbee)
! 5: TVD with MUSCL limiter ! 5: TVD with MUSCL limiter - case (MUSCL)
! 6: TVD with ULTIMATE QUICKEST ! 6: TVD with ULTIMATE QUICKEST - case ((P2),(P2_PDM))
! !
! ode_method -> ODE scheme for source and sink dynamics ! ode_method -> ODE scheme for source and sink dynamics
! 1: first-order explicit (not positive) ! 1: first-order explicit (not positive)
...@@ -45,16 +47,17 @@ ...@@ -45,16 +47,17 @@
! bio_lagrange_mean -> averaging Lagrangian conc. on output (.true./.false.) ! bio_lagrange_mean -> averaging Lagrangian conc. on output (.true./.false.)
! !
! bio_npar -> total number of Lagrangian particles ! bio_npar -> total number of Lagrangian particles
!
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
&bio_nml &bio_nml
bio_calc= .true. bio_calc= .true.
bio_model= 2 bio_model= 10
bio_eulerian= .true. bio_eulerian= .true.
cnpar= 1.0 cnpar= 1.0
w_adv_discr= 6 w_adv_discr= 6 ! 4
ode_method= 8 ode_method= 10 !11
split_factor= 1 split_factor= 1
bioshade_feedback= .true. bioshade_feedback= .true.
bio_lagrange_mean= .false. bio_lagrange_mean= .true.
bio_npar= 100000 bio_npar= 10000
/ /
nml/bio_polynow.nml
View file @ 20015aff
This diff is collapsed.
src/extras/bio/bio.F90
View file @ 20015aff
...@@ -48,8 +48,8 @@ ...@@ -48,8 +48,8 @@
use bio_mab, only : init_bio_mab,init_var_mab,var_info_mab use bio_mab, only : init_bio_mab,init_var_mab,var_info_mab
use bio_mab, only : light_mab,surface_fluxes_mab,do_bio_mab use bio_mab, only : light_mab,surface_fluxes_mab,do_bio_mab
use bio_polynow, only : init_bio_polynow,init_var_polynow,var_info_polynow !GG-A use bio_polynow, only : init_bio_polynow,init_var_polynow,var_info_polynow
use bio_polynow, only : light_polynow,do_bio_polynow !GG-A use bio_polynow, only : light_polynow,do_bio_polynow
use output, only : out_fmt,write_results,ts use output, only : out_fmt,write_results,ts
...@@ -360,7 +360,7 @@ ...@@ -360,7 +360,7 @@
call var_info_npzd4() call var_info_npzd4()
case (10) ! The model for the north Water polynya - NPZD- Sedimentation GG-A case (10) ! The model for the north Water polynya - NPZ4D- Sedimentation
call init_bio_polynow(namlst,'bio_polynow.nml',unit) call init_bio_polynow(namlst,'bio_polynow.nml',unit)
...@@ -514,9 +514,6 @@ ...@@ -514,9 +514,6 @@
if (present(w_)) w = w_ if (present(w_)) w = w_
if (present(w_adv_ctr_)) w_adv_ctr = w_adv_ctr_ if (present(w_adv_ctr_)) w_adv_ctr = w_adv_ctr_
return return
end subroutine set_env_bio end subroutine set_env_bio
!EOC !EOC
...@@ -594,6 +591,7 @@ ...@@ -594,6 +591,7 @@
! !LOCAL VARIABLES: ! !LOCAL VARIABLES:
integer, parameter :: adv_mode_0=0 integer, parameter :: adv_mode_0=0
integer, parameter :: adv_mode_1=1 integer, parameter :: adv_mode_1=1
integer, parameter :: adv_mode_3=3 ! Specific for Bio_model =10 , for Msn concentration and size
REALTYPE :: Qsour(0:nlev),Lsour(0:nlev) REALTYPE :: Qsour(0:nlev),Lsour(0:nlev)
REALTYPE :: RelaxTau(0:nlev) REALTYPE :: RelaxTau(0:nlev)
REALTYPE :: dt_eff REALTYPE :: dt_eff
...@@ -670,8 +668,29 @@ ...@@ -670,8 +668,29 @@
do i=1,nlev do i=1,nlev
cc(1,i) = nit(i) cc(1,i) = nit(i)
end do end do
else
! do advection step due to settling or rising ! Specific advection scheme for Bio_model =10 , for Msn concentration (10) and size (14)
else if (j==10 .and. bio_model ==10) then
call adv_center(nlev,dt,h,h,ws(j,:),flux, &
flux,_ZERO_,_ZERO_,w_adv_discr,adv_mode_3,cc(j,:))
else if (j==14 .and. bio_model ==10) then
call adv_center(nlev,dt,h,h,ws(j,:),flux, &
flux,_ZERO_,_ZERO_,w_adv_discr,adv_mode_3,cc(j,:))
! Specific advection scheme for Bio_model =10 , for nitrate (5) and ammonium (6)
else if (j ==5 .and. bio_model==10) then
do i=1,nlev
cc(j,i) = nit(i)
end do
else if (j ==6 .and. bio_model==10) then
do i=1,nlev
cc(j,i) = amm(i)
end do
else
! do advection step due to settling or rising for all the other variables
call adv_center(nlev,dt,h,h,ws(j,:),flux, & call adv_center(nlev,dt,h,h,ws(j,:),flux, &
flux,_ZERO_,_ZERO_,w_adv_discr,adv_mode_1,cc(j,:)) flux,_ZERO_,_ZERO_,w_adv_discr,adv_mode_1,cc(j,:))
...@@ -873,9 +892,6 @@ ...@@ -873,9 +892,6 @@
if (allocated(Flux_D2)) deallocate(Flux_D2) if (allocated(Flux_D2)) deallocate(Flux_D2)
if (allocated(Flux_D3)) deallocate(Flux_D3) if (allocated(Flux_D3)) deallocate(Flux_D3)
if (allocated(size_msnow)) deallocate(size_msnow)
if (allocated(w_msn_lev)) deallocate(w_msn_lev)
! The external provide arrays ! The external provide arrays
if (allocated(h)) deallocate(h) if (allocated(h)) deallocate(h)
if (allocated(nuh)) deallocate(nuh) if (allocated(nuh)) deallocate(nuh)
...@@ -1005,8 +1021,6 @@ ...@@ -1005,8 +1021,6 @@
allocate(ppnet(0:nlev),stat=rc) allocate(ppnet(0:nlev),stat=rc)
if (rc /= 0) stop 'init_bio(): Error allocating (ppnet)' if (rc /= 0) stop 'init_bio(): Error allocating (ppnet)'
allocate(flux_msn(0:nlev),stat=rc) allocate(flux_msn(0:nlev),stat=rc)
if (rc /= 0) stop 'init_bio(): Error allocating (flux_msn)' if (rc /= 0) stop 'init_bio(): Error allocating (flux_msn)'
...@@ -1019,14 +1033,6 @@ ...@@ -1019,14 +1033,6 @@
allocate(Flux_D3(0:nlev),stat=rc) allocate(Flux_D3(0:nlev),stat=rc)
if (rc /= 0) stop 'init_bio(): Error allocating (Flux_D3)' if (rc /= 0) stop 'init_bio(): Error allocating (Flux_D3)'
allocate(size_msnow(0:nlev),stat=rc)
if (rc /= 0) stop 'init_bio(): Error allocating (size_msnow)'
allocate(w_msn_lev(0:nlev),stat=rc)
if (rc /= 0) stop 'init_bio(): Error allocating (w_msn_lev)'
! The external provide arrays ! The external provide arrays
allocate(h(0:nlev),stat=rc) allocate(h(0:nlev),stat=rc)
if (rc /= 0) stop 'init_bio(): Error allocating (h)' if (rc /= 0) stop 'init_bio(): Error allocating (h)'
......
src/extras/bio/bio_polynow.F90
View file @ 20015aff
This diff is collapsed.
src/extras/bio/bio_save.F90
View file @ 20015aff
...@@ -93,7 +93,7 @@ ...@@ -93,7 +93,7 @@
!----------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------
!Prepare the new variable for the NETCDF output file as well as informations !Prepare the new variable for the NETCDF output file as well as informations
!----------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------
!Density of the fluid at each level !Density of the fluid at each level
iret = new_nc_variable(ncid,'rho',NF_REAL,4,dims,rho_id) iret = new_nc_variable(ncid,'rho',NF_REAL,4,dims,rho_id)
iret = set_attributes(ncid,rho_id,units='kg.m-3',long_name='Density') iret = set_attributes(ncid,rho_id,units='kg.m-3',long_name='Density')
...@@ -218,6 +218,7 @@ ...@@ -218,6 +218,7 @@
! Flux of FP going to msn ! Flux of FP going to msn
iret = new_nc_variable(ncid,'Flux_D3',NF_REAL,4,dims,Flux_D3_id) iret = new_nc_variable(ncid,'Flux_D3',NF_REAL,4,dims,Flux_D3_id)
iret = set_attributes(ncid,Flux_D3_id,units='1/day',long_name='Flux_D3') iret = set_attributes(ncid,Flux_D3_id,units='1/day',long_name='Flux_D3')
endif endif
!DD Diagnostic de npar (nb de particules lagrangiennes) pour bebogage !DD Diagnostic de npar (nb de particules lagrangiennes) pour bebogage
...@@ -240,11 +241,9 @@ ...@@ -240,11 +241,9 @@
iret = store_data(ncid,var_ids(n),XYZT_SHAPE,nlev,array=cc(n,:)) iret = store_data(ncid,var_ids(n),XYZT_SHAPE,nlev,array=cc(n,:))
end do end do
!Density of the fluid at each level !Density of the fluid at each level
iret = store_data(ncid,rho_id,XYZT_SHAPE,nlev,array=rho) iret = store_data(ncid,rho_id,XYZT_SHAPE,nlev,array=rho)
! Sedimentation rate of phytoplankton ! Sedimentation rate of phytoplankton
iret = store_data(ncid,wp_id,XYZT_SHAPE,nlev,array=secs_pr_day*ws(1,:)) iret = store_data(ncid,wp_id,XYZT_SHAPE,nlev,array=secs_pr_day*ws(1,:))
...@@ -259,7 +258,7 @@ ...@@ -259,7 +258,7 @@
if (bio_model .eq. 10) then if (bio_model .eq. 10) then
! Here it is needed to multiply by 86400(secs_pr_day) in order to have the model ouptu (in s) convert in per day ! Here it is needed to multiply by 86400(secs_pr_day) in order to have the model ouptut (in s) convert in per day
! Living zooplankton ! Living zooplankton
iret = store_data(ncid,wz_id,XYZT_SHAPE,nlev,array=secs_pr_day*ws(2,:)) iret = store_data(ncid,wz_id,XYZT_SHAPE,nlev,array=secs_pr_day*ws(2,:))
! Dead Phytoplankton ! Dead Phytoplankton
...@@ -274,36 +273,24 @@ ...@@ -274,36 +273,24 @@
! Stickiness ! Stickiness
! 2 Living phytoplankton ! 2 Living phytoplankton
iret = store_data(ncid,sti_2p_id,XYZT_SHAPE,nlev,array=sti_2p(:)) iret = store_data(ncid,sti_2p_id,XYZT_SHAPE,nlev,array=sti_2p(:))
! Living Phytoplankton with dead phytoplankton ! Living Phytoplankton with dead phytoplankton
iret = store_data(ncid,sti_pdph_id,XYZT_SHAPE,nlev,array=sti_pdph(:)) iret = store_data(ncid,sti_pdph_id,XYZT_SHAPE,nlev,array=sti_pdph(:))
! Living Phytoplankton with dead zooplankton ! Living Phytoplankton with dead zooplankton
iret = store_data(ncid,sti_pdzo_id,XYZT_SHAPE,nlev,array=sti_pdzo(:)) iret = store_data(ncid,sti_pdzo_id,XYZT_SHAPE,nlev,array=sti_pdzo(:))
! Living Phytoplankton with fecal pellets ! Living Phytoplankton with fecal pellets
iret = store_data(ncid,sti_pfp_id,XYZT_SHAPE,nlev,array=sti_pfp(:)) iret = store_data(ncid,sti_pfp_id,XYZT_SHAPE,nlev,array=sti_pfp(:))
! 2 Dead phytoplankton ! 2 Dead phytoplankton
iret = store_data(ncid,sti_2dph_id,XYZT_SHAPE,nlev,array=sti_2dph(:)) iret = store_data(ncid,sti_2dph_id,XYZT_SHAPE,nlev,array=sti_2dph(:))
! Dead Phytoplankton with dead zooplankton ! Dead Phytoplankton with dead zooplankton
iret = store_data(ncid,sti_dphdzo_id,XYZT_SHAPE,nlev,array=sti_dphdzo(:)) iret = store_data(ncid,sti_dphdzo_id,XYZT_SHAPE,nlev,array=sti_dphdzo(:))
! Dead Phytoplankton with fecal pellets ! Dead Phytoplankton with fecal pellets
iret = store_data(ncid,sti_dphfp_id,XYZT_SHAPE,nlev,array=sti_dphfp(:)) iret = store_data(ncid,sti_dphfp_id,XYZT_SHAPE,nlev,array=sti_dphfp(:))
! 2 Dead zooplankton ! 2 Dead zooplankton
iret = store_data(ncid,sti_2dzo_id,XYZT_SHAPE,nlev,array=sti_2dzo(:)) iret = store_data(ncid,sti_2dzo_id,XYZT_SHAPE,nlev,array=sti_2dzo(:))
! Dead zooplankton with fecal pellets ! Dead zooplankton with fecal pellets
iret = store_data(ncid,sti_dzofp_id,XYZT_SHAPE,nlev,array=sti_dzofp(:)) iret = store_data(ncid,sti_dzofp_id,XYZT_SHAPE,nlev,array=sti_dzofp(:))
! 2 fecal pellets ! 2 fecal pellets
iret = store_data(ncid,sti_2fp_id,XYZT_SHAPE,nlev,array=sti_2fp(:)) iret = store_data(ncid,sti_2fp_id,XYZT_SHAPE,nlev,array=sti_2fp(:))
! Rho_F
! iret = store_data(ncid,rho_F_id,XYZT_SHAPE,nlev,array=)
end if end if
...@@ -319,6 +306,7 @@ ...@@ -319,6 +306,7 @@
iret = store_data(ncid,ammlim2_id,XYZT_SHAPE,nlev,array=ammlim2(:)) iret = store_data(ncid,ammlim2_id,XYZT_SHAPE,nlev,array=ammlim2(:))
iret = store_data(ncid,ppnet_id,XYZT_SHAPE,nlev,array=ppnet(:)) iret = store_data(ncid,ppnet_id,XYZT_SHAPE,nlev,array=ppnet(:))
!-----Coagulation fluxes
if (bio_model .eq. 10) then if (bio_model .eq. 10) then
iret = store_data(ncid,flux_msn_id,XYZT_SHAPE,nlev,array=secs_pr_day*flux_msn(:)) iret = store_data(ncid,flux_msn_id,XYZT_SHAPE,nlev,array=secs_pr_day*flux_msn(:))
iret = store_data(ncid,Flux_P_id,XYZT_SHAPE,nlev,array=secs_pr_day*Flux_P(:)) iret = store_data(ncid,Flux_P_id,XYZT_SHAPE,nlev,array=secs_pr_day*Flux_P(:))
......
src/extras/bio/bio_var.F90
View file @ 20015aff
...@@ -17,34 +17,54 @@ ...@@ -17,34 +17,54 @@
public public
! !
! !PUBLIC DATA MEMBERS: ! !PUBLIC DATA MEMBERS:
!-----Model choice
integer :: bio_model integer :: bio_model
!-----Model variable number
integer :: numc,numcc integer :: numc,numcc
!-----Variables names
integer, dimension(:), allocatable :: var_ids
character(len=64), dimension(:), allocatable :: var_names
character(len=64), dimension(:), allocatable :: var_units
character(len=64), dimension(:), allocatable :: var_long
!-----?
REALTYPE, dimension(:), allocatable :: zlev REALTYPE, dimension(:), allocatable :: zlev
!-----PAR radiation
REALTYPE, dimension(:), allocatable :: par REALTYPE, dimension(:), allocatable :: par
REALTYPE, dimension(:), allocatable :: lumlim1,nitlim1,ammlim1,lumlim2,nitlim2,ammlim2 !DD integer :: par_id
!-----Nutrients limitations
REALTYPE, dimension(:), allocatable :: lumlim1,nitlim1,ammlim1,lumlim2,nitlim2,ammlim2
integer :: lumlim1_id,nitlim1_id,ammlim1_id,lumlim2_id,nitlim2_id,ammlim2_id
!-----Primary Production
REALTYPE, dimension(:), allocatable :: ppnet REALTYPE, dimension(:), allocatable :: ppnet
integer :: ppnet_id
!-----Coagulation fluxes
REALTYPE, dimension(:), allocatable :: flux_msn REALTYPE, dimension(:), allocatable :: flux_msn
integer :: flux_msn_id integer :: flux_msn_id
REALTYPE, dimension(:), allocatable :: Flux_P REALTYPE, dimension(:), allocatable :: Flux_P
integer :: Flux_P_id integer :: Flux_P_id
REALTYPE, dimension(:), allocatable :: Flux_D1 REALTYPE, dimension(:), allocatable :: Flux_D1
integer :: Flux_D1_id integer :: Flux_D1_id
REALTYPE, dimension(:), allocatable :: Flux_D2 REALTYPE, dimension(:), allocatable :: Flux_D2
integer :: Flux_D2_id integer :: Flux_D2_id
REALTYPE, dimension(:), allocatable :: Flux_D3 REALTYPE, dimension(:), allocatable :: Flux_D3
integer :: Flux_D3_id integer :: Flux_D3_id
!-----Stickiness
REALTYPE, dimension(:), allocatable :: sti_2p,sti_pdph,sti_pdzo,sti_pfp
integer :: sti_2p_id,sti_pdph_id,sti_pdzo_id,sti_pfp_id
REALTYPE, dimension(:), allocatable :: sti_2dph,sti_dphdzo,sti_dphfp
integer :: sti_2dph_id,sti_dphdzo_id,sti_dphfp_id
REALTYPE, dimension(:), allocatable :: sti_2dzo,sti_dzofp
integer :: sti_2dzo_id,sti_dzofp_id
REALTYPE, dimension(:), allocatable :: sti_2fp
integer :: sti_2fp_id
REALTYPE, dimension(:), allocatable :: size_msnow !-----?
integer :: size_msnow_id REALTYPE :: dt_eff
!-----Concentrations (cc) and Advection (ws)
REALTYPE, dimension(:), allocatable :: w_msn_lev
integer :: w_msn_lev_id
REALTYPE :: dt_eff
REALTYPE, dimension(:,:), allocatable :: cc,ws REALTYPE, dimension(:,:), allocatable :: cc,ws
integer :: wp_id,wz_id,wd1_id,wd2_id,wd3_id,wd4_id
!-----?
integer :: surface_flux_method=-1 integer :: surface_flux_method=-1
integer :: n_surface_fluxes=-1 integer :: n_surface_fluxes=-1
REALTYPE, dimension(:), allocatable :: sfl_read REALTYPE, dimension(:), allocatable :: sfl_read
...@@ -55,28 +75,7 @@ ...@@ -55,28 +75,7 @@
integer, dimension(:,:), allocatable :: particle_indx integer, dimension(:,:), allocatable :: particle_indx
REALTYPE, dimension(:,:), allocatable :: particle_pos REALTYPE, dimension(:,:), allocatable :: particle_pos
!-----Time unites
!Stickiness
REALTYPE, dimension(:), allocatable :: sti_2p,sti_pdph,sti_pdzo,sti_pfp
integer :: sti_2p_id,sti_pdph_id,sti_pdzo_id,sti_pfp_id
REALTYPE, dimension(:), allocatable :: sti_2dph,sti_dphdzo,sti_dphfp
integer :: sti_2dph_id,sti_dphdzo_id,sti_dphfp_id
REALTYPE, dimension(:), allocatable :: sti_2dzo,sti_dzofp
integer :: sti_2dzo_id,sti_dzofp_id
REALTYPE, dimension(:), allocatable :: sti_2fp
integer :: sti_2fp_id
integer, dimension(:), allocatable :: var_ids
integer :: wp_id,wz_id,wd1_id,wd2_id,wd3_id,wd4_id
integer :: par_id
integer :: lumlim1_id,nitlim1_id,ammlim1_id,lumlim2_id,nitlim2_id,ammlim2_id !DD
integer :: ppnet_id
character(len=64), dimension(:), allocatable :: var_names
character(len=64), dimension(:), allocatable :: var_units
character(len=64), dimension(:), allocatable :: var_long
REALTYPE, parameter :: secs_pr_day=86400. REALTYPE, parameter :: secs_pr_day=86400.
! external variables - i.e. provided by the calling program but ! external variables - i.e. provided by the calling program but
...@@ -85,10 +84,10 @@ ...@@ -85,10 +84,10 @@
REALTYPE, dimension(:), allocatable :: h REALTYPE, dimension(:), allocatable :: h
REALTYPE, dimension(:), allocatable :: t REALTYPE, dimension(:), allocatable :: t
REALTYPE, dimension(:), allocatable :: s REALTYPE, dimension(:), allocatable :: s
REALTYPE, dimension(:), allocatable :: nit !CHG3 REALTYPE, dimension(:), allocatable :: nit
REALTYPE, dimension(:), allocatable :: nprof !CHG3 REALTYPE, dimension(:), allocatable :: nprof
REALTYPE, dimension(:), allocatable :: amm !CHG5 REALTYPE, dimension(:), allocatable :: amm
REALTYPE, dimension(:), allocatable :: aprof !CHG5 REALTYPE, dimension(:), allocatable :: aprof
REALTYPE, dimension(:), allocatable :: hcb REALTYPE, dimension(:), allocatable :: hcb
REALTYPE, dimension(:), allocatable :: hcprof REALTYPE, dimension(:), allocatable :: hcprof
REALTYPE, dimension(:), allocatable :: rho REALTYPE, dimension(:), allocatable :: rho
...@@ -98,7 +97,7 @@ ...@@ -98,7 +97,7 @@
REALTYPE, dimension(:), allocatable :: rad REALTYPE, dimension(:), allocatable :: rad
REALTYPE :: wind REALTYPE :: wind
REALTYPE :: I_0 REALTYPE :: I_0
integer :: w_adv_ctr=0 integer :: w_adv_ctr= 0
! external variables updated by the biological models ! external variables updated by the biological models
! the variables are copied back to the calling program using ! the variables are copied back to the calling program using
......
src/util/adv_center.F90
View file @ 20015aff
...@@ -267,6 +267,34 @@ ...@@ -267,6 +267,34 @@
! vertical loop ! vertical loop
do k=1,N-1 do k=1,N-1
! --- For the model bio_polynow (Bio_model = 10 in bio.nml)
! a specific advection scheme is set for the marine snow concentration and size variables.
! This is refer in the following by mode =3 as set in bio.F90
if (mode.eq.3) then
! compute the slope ration
! for negative speed
! compute courant number
c=-ww(k)/float(it)*dt/(0.5*(h(k)+h(k+1)))
if (k .lt. N-1) then
Yu=Y(k+2) ! upstream value
else
Yu=Y(k+1)
endif
Yc=Y(k+1) ! central value
Yd=Y(k ) ! downstream value
! compute slope ratio
if (abs(Yc-Yd) .gt. 1e-10) then
r=(Yu-Yc)/(Yc-Yd)
else
r=(Yu-Yc)*1.e10
endif
else
! compute the slope ration ! compute the slope ration
if (ww(k) .gt. _ZERO_) then if (ww(k) .gt. _ZERO_) then
...@@ -312,6 +340,8 @@ ...@@ -312,6 +340,8 @@
end if end if
endif
! compute the flux-factor phi ! compute the flux-factor phi
x = one6th*(1.-2.0*c) x = one6th*(1.-2.0*c)
Phi = (0.5+x)+(0.5-x)*r Phi = (0.5+x)+(0.5-x)*r
...@@ -339,10 +369,20 @@ ...@@ -339,10 +369,20 @@
stop stop
end select end select
if (mode.eq.3) then
! --- For the model bio_polynow (Bio_model = 10 in bio.nml) for Msn size and concentration
! compute limited flux
cu(k) =ww(k)*(Yc+0.5*limit*(1-c)*(Yd-Yc))
cu(k-1)=ww(k)*(Yd+0.5*limit*(1-c)*(Yd-Y(k-1)))
else ! For the other variables of bio_polynow, or other biological models
! compute the limited flux ! compute the limited flux
cu(k)=ww(k)*(Yc+0.5*limit*(1-c)*(Yd-Yc)) cu(k)=ww(k)*(Yc+0.5*limit*(1-c)*(Yd-Yc))
endif
end do end do ! End of vertical loop
! do the upper boundary conditions ! do the upper boundary conditions
select case (Bcup) select case (Bcup)
...@@ -391,6 +431,14 @@ ...@@ -391,6 +431,14 @@
Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/ & Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/ &
h(k)-Y(k)*(ww(k)-ww(k-1))/h(k)) h(k)-Y(k)*(ww(k)-ww(k-1))/h(k))
enddo enddo
! --- For the model bio_polynow (Bio_model = 10 in bio.nml) for Msn size and concentration
elseif (mode.eq.3) then ! conservative
do k=1,N
! advection step
Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k))
enddo
! For the other variables of bio_polynow, or other biological models
else ! conservative else ! conservative
do k=1,N do k=1,N
Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k)) Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k))
......