#undef NEWMAN #define LIMIT_INTERIOR MODULE biology_mod ! !git $Id$ !svn $Id: red_tide.h 1151 2023-02-09 03:08:53Z arango $ !======================================================== Ruoying He === ! Copyright (c) 2002-2023 The ROMS/TOMS Group ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md Hernan G. Arango ! !========================================== Alexander F. Shchepetkin === ! ! ! Red Tide Biological Model: Alexandrium fundyense ! ! ! ! This routine computes the biological sources and sinks and ! ! updates the global tracer array for dinoflagellates ! ! ! ! References: ! ! ! ! Stock, C.A., D.J. McGillicudy, A.R. Solow, and D.M. Anderson, ! ! 2005: Evaluating hypotheses for the initiation and development ! ! of Alexandrium fundyense blooms in the western Gulf of Maine ! ! using a coupled physical-biological model, Deep-Sea Research II,! ! 52, 2715-2744. ! ! ! ! He, R., D.J. McGillicuddy, B.A. Keafer, and D.M. Anderson, 2008: ! ! Historic 2005 toxic bloom of Alexandrium fundyense in the ! ! western Gulf of Maine: 2, Coupled biophysical modeling, J. ! ! Geophys. Res., 113, C07040, doi:10.1029/2007JC004602. ! ! ! !======================================================================= ! implicit none ! PRIVATE PUBLIC :: biology ! CONTAINS ! !*********************************************************************** SUBROUTINE biology (ng,tile) !*********************************************************************** ! USE mod_param USE mod_forces USE mod_grid USE mod_ncparam USE mod_ocean USE mod_stepping ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__ ! #include "tile.h" ! ! Set header file name. ! #ifdef DISTRIBUTE IF (Lbiofile(iNLM)) THEN #else IF (Lbiofile(iNLM).and.(tile.eq.0)) THEN #endif Lbiofile(iNLM)=.FALSE. BIONAME(iNLM)=MyFile END IF ! #ifdef PROFILE CALL wclock_on (ng, iNLM, 15, __LINE__, MyFile) #endif CALL red_tide_tile (ng, tile, & & LBi, UBi, LBj, UBj, N(ng), NT(ng), & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nnew(ng), & #ifdef MASKING & GRID(ng) % rmask, & #endif & GRID(ng) % Hz, & & GRID(ng) % z_r, & & GRID(ng) % z_w, & #ifdef DAILY_SHORTWAVE & FORCES(ng) % srflx_avg, & #endif & FORCES(ng) % srflx, & & OCEAN(ng) % CystIni, & & OCEAN(ng) % DIN_obs, & & OCEAN(ng) % t) #ifdef PROFILE CALL wclock_off (ng, iNLM, 15, __LINE__, MyFile) #endif ! RETURN END SUBROUTINE biology ! !----------------------------------------------------------------------- SUBROUTINE red_tide_tile (ng, tile, & & LBi, UBi, LBj, UBj, UBk, UBt, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & #ifdef MASKING & rmask, & #endif & Hz, z_r, z_w, & #ifdef DAILY_SHORTWAVE & srflx_avg, & #endif & srflx, CystIni, DIN_obs, & & t) !----------------------------------------------------------------------- ! USE mod_param USE mod_biology USE mod_ncparam USE mod_scalars ! USE dateclock_mod, ONLY : caldate ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nnew #ifdef ASSUMED_SHAPE # ifdef MASKING real(r8), intent(in) :: rmask(LBi:,LBj:) # endif real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(in) :: z_r(LBi:,LBj:,:) real(r8), intent(in) :: z_w(LBi:,LBj:,0:) real(r8), intent(in) :: CystIni(LBi:,LBj:) real(r8), intent(in) :: DIN_obs(LBi:,LBj:,:) # ifdef DAILY_SHORTWAVE real(r8), intent(in) :: srflx_avg(LBi:,LBj:) # endif real(r8), intent(in) :: srflx(LBi:,LBj:) real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:) #else # ifdef MASKING real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk) real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk) real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk) real(r8), intent(in) :: CystIni(LBi:UBi,LBj:UBj) # ifdef DAILY_SHORTWAVE real(r8), intent(in) :: srflx_avg(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj) real(r8), intent(in) :: DIN_ob(LBi:UBi,LBj:UBj,UBk) real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt) #endif ! ! Local variable declarations. ! integer, parameter :: Nswim = 1 integer, parameter :: NsedLayers = 10 integer :: Iter, i, ibio, iswim, itrc, j, k, ks, ksed integer, dimension(Nswim) :: idswim real(r8) :: Cell_Flux, C_depth, DIN, E_flux, EndoScale real(r8) :: Rad, RadScale real(r8) :: GermD, GermL, G_DIN, G_light, G_rate, M_rate real(r8) :: G_fac, S_fac, T_fac real(r8) :: dtdays, oNsedLayers, salt, temp, wmig real(dp) :: yday real(r8) :: alpha, cff, cffL, cffR, deltaL, deltaR, dz, wdt real(r8), parameter :: eps = 1.0E-8_r8 real(r8), dimension(Nswim) :: Wbio integer, dimension(IminS:ImaxS,N(ng)) :: ksource real(r8), dimension(IminS:ImaxS) :: Germ real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC real(r8), dimension(IminS:ImaxS,0:N(ng)) :: aL real(r8), dimension(IminS:ImaxS,0:N(ng)) :: aR real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dL real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dR real(r8), dimension(IminS:ImaxS,0:N(ng)) :: r real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv real(r8), dimension(IminS:ImaxS,N(ng)) :: Light ! ! Temperature dependent grow factor polynomial coefficients. ! !! real(r8), parameter :: TC0 = 0.382_r8 ! Stock el al, 2005 !! real(r8), parameter :: TC1 =-0.0867_r8 ! Eq. 5a !! real(r8), parameter :: TC2 = 0.0160_r8 !! real(r8), parameter :: TC3 =-0.000513_r8 !! real(r8), parameter :: TC0 = 0.379_r8 ! Revised real(r8), parameter :: TC1 =-0.0961_r8 ! Stock 8/15/2006 real(r8), parameter :: TC2 = 0.0169_r8 real(r8), parameter :: TC3 =-0.000536_r8 ! ! Salinity dependent grow factor polynomial coefficients. ! !! real(r8), parameter :: SC0 =-0.872_r8 ! Stock el al, 2005 !! real(r8), parameter :: SC1 = 0.220_r8 ! Eq. 6 !! real(r8), parameter :: SC2 =-0.00808_r8 !! real(r8), parameter :: SC3 = 0.0000882_r8 !! real(r8), parameter :: SC0 =-0.693_r8 ! Revised real(r8), parameter :: SC1 = 0.186_r8 ! Stock 8/15/2006 real(r8), parameter :: SC2 =-0.00622_r8 real(r8), parameter :: SC3 = 0.0000557_r8 #include "set_bounds.h" ! !----------------------------------------------------------------------- ! Add biological Source/Sink terms. !----------------------------------------------------------------------- ! ! Avoid computing source/sink terms if no biological iterations. ! IF (BioIter(ng).le.0) RETURN ! ! Set time-stepping size (days) according to the number of iterations. ! dtdays=dt(ng)*sec2day/REAL(BioIter(ng),r8) ! ! Set shortwave radiation scale. In ROMS all the fluxes are kinematic. ! RadScale=rho0*Cp ! Celsius m/s to Watts/m2 ! ! Set critical depth (m; negative) used in the growth function. ! C_depth=(LOG(G_r(ng)/(G_eff(ng)*srad_Cdepth(ng))))/AttW(ng) ! ! Set vertical swimming identification vector. ! idswim(1)=iDino ! Dinoflagellate ! ! Set vertical swimming velocity vector in the same order as the ! identification vector, IDSWIM. ! Wbio(1)=wDino(ng) ! Dinoflagellate ! ! Set scale for germination term. ! oNsedLayers=1.0_r8/REAL(NsedLayers,r8) ! ! Compute inverse thickness to avoid repeated divisions. ! J_LOOP : DO j=Jstr,Jend DO k=1,N(ng) DO i=Istr,Iend Hz_inv(i,k)=1.0_r8/Hz(i,j,k) END DO END DO ! ! Extract biological variables from tracer arrays, place them into ! scratch arrays, and restrict their values to be positive definite. ! At input, all tracers (index nnew) from predictor step have ! transport units (m Tunits) since we do not have yet the new ! values for zeta and Hz. These are known after the 2D barotropic ! time-stepping. ! DO itrc=1,NBT ibio=idbio(itrc) DO k=1,N(ng) DO i=Istr,Iend Bio_old(i,k,ibio)=MAX(0.0_r8,t(i,j,k,nstp,ibio)) Bio(i,k,ibio)=Bio_old(i,k,ibio) END DO END DO END DO ! ! Extract potential temperature and salinity. ! DO k=1,N(ng) DO i=Istr,Iend Bio(i,k,itemp)=MIN(t(i,j,k,nstp,itemp),36.0_r8) Bio(i,k,isalt)=MAX(0.0_r8,t(i,j,k,nstp,isalt)) END DO END DO ! ! Calculate endogenous clock scale used in the cysts germination ! term. The cysts germination rates are regulated by an endogenous ! circannual clock. ! CALL caldate (tdays(ng), yd_dp=yday) ! IF (yday.lt.Month_MidDay(1)) THEN cff=(365.0_r8-Month_MidDay(12)+yday)/ & & (365.0_r8-Month_MidDay(12)+Month_MidDay(1)) EndoScale=GPN(12)+cff*(GPN(1)-GPN(12)) ELSE IF (yday.ge.Month_MidDay(12)) THEN cff=(yday-Month_MidDay(12))/ & (365.0_r8-Month_MidDay(12)+Month_MidDay(1)) EndoScale=GPN(12)+cff*(GPN(1)-GPN(12)) ELSE DO i=1,11 IF ((yday.ge.Month_MidDay(i)).and. & & (yday.lt.Month_MidDay(i+1))) THEN cff=(yday-Month_MidDay(i))/ & & (Month_MidDay(i+1)-Month_MidDay(i)) EndoScale=GPN(i)+cff*(GPN(i+1)-GPN(i)) END IF END DO END IF ! !======================================================================= ! Start internal iterations to achieve convergence of the nonlinear ! backward-implicit solution. !======================================================================= ! ! The iterative loop below is to iterate toward an universal Backward- ! Euler treatment of all terms. So if there are oscillations in the ! system, they are only physical oscillations. These iterations, ! however, do not improve the accuaracy of the solution. ! ITER_LOOP: DO Iter=1,BioIter(ng) ! !----------------------------------------------------------------------- ! Add Cyst germination flux at the bottom layer !----------------------------------------------------------------------- ! ! Calculate Cyst germination rate at the top of the sediment layer ! as a function of bottom water temperature and non-spectral ! irradiance. ! DO i=Istr,Iend ! ! Calculate "light" and "dark" cyst germination rates as a function ! of bottom temperature. ! temp=Bio(i,1,itemp) ! bottom, k=1 GermL=(1.50_r8+ & & (8.72_r8-1.50_r8)*0.5_r8* & & (TANH(0.790_r8*temp-6.27_r8)+1.0_r8))*oNsedLayers GermD=(1.04_r8+ & & (4.26_r8-1.04_r8)*0.5_r8* & & (TANH(0.394_r8*temp-3.33_r8)+1.0_r8))*oNsedLayers ! ! Compute non-spectral irradiance flux at each sediment layer. Then, ! compute cyst germination rate according to the light regime. ! Germ(i)=0.0_r8 ! initialize DO ksed=1,NsedLayers # ifdef DAILY_SHORTWAVE E_flux=RadScale*srflx_avg(i,j)* & & EXP( AttW(ng)*z_w(i,j,0)- & & AttS(ng)*Dg(ng)*(REAL(ksed,r8)-0.5) ) # else E_flux=RadScale*srflx(i,j)* & & EXP( AttW(ng)*z_w(i,j,0)- & & AttS(ng)*Dg(ng)*(REAL(ksed,r8)-0.5) ) # endif IF (E_flux.gt.E_light(ng)) THEN Germ(i)=Germ(i)+GermL ELSE IF (E_flux.lt.E_dark(ng)) THEN Germ(i)=Germ(i)+GermD ELSE Germ(i)=Germ(i)+ & & (GermD+ & & (GermL-GermD)* & & ((E_flux-E_dark(ng))/ & & (E_light(ng)-E_dark(ng)))) END IF END DO ! ! Multiply by endogenous clock factor. The cyst germination are ! regulated by an endogenous circannual clock. The 100 factor here ! is because "Dg" is meters and we need centimeters. ! Germ(i)=Germ(i)*Dg(ng)*100.0_r8*EndoScale ! ! Convert percentage cysts/day into decimal fraction of cysts. ! Germ(i)=Germ(i)*0.01_r8 ! ! Calculate the flux of cells away from the bottom. It is referenced ! to the initial number of cysts to be consistent with laboratory ! experiments. ! Cell_Flux=CystIni(i,j)* & & Germ(i)*Hz_inv(i,1) ! cells/m3 ! ! Add cell flux at the bottom layer (k=1). ! Bio(i,1,iDino)=Bio(i,1,iDino)+Cell_Flux*dtdays END DO ! !----------------------------------------------------------------------- ! Compute growth term. !----------------------------------------------------------------------- ! ! The growth is dependent on temperature, salinity, non-spectral ! irradiance (light), and nutrient (Dissolved Inorganic Nutrient, DIN). ! DO k=1,N(ng) DO i=Istr,Iend temp=Bio(i,k,itemp) salt=Bio(i,k,isalt) ! ! Compute Alexandrium fundyense temperature dependent growth factor ! using a cubic polynomial fitted to available data. ! IF (temp.ge.Tmin_growth(ng)) THEN T_fac=TC0+temp*(TC1+temp*(TC2+temp*TC3)) ELSE ! linear extrapolation !! T_fac=TC0+temp*(TC1+temp*(TC22+temp*TC3))- & !! & 0.0343_r8*(5.0_r8-temp) ! Stock el al, 2005 !! ! Eq. 5b T_fac=0.254_r8-0.0327_r8*(5.0_r8-temp) ! Stock 8/15/2006 END IF ! ! Compute Alexandrium fundyense salinity dependent growth factor ! using a cubic polynomial to fit to available data. ! S_fac=SC0+salt*(SC1+salt*(SC2+salt*SC3)) ! ! Compute temperature and salinity growth factor. ! G_fac=T_fac*S_fac ! ! Compute light dependency factor (Platt and Jassby, 1976). ! # ifdef DAILY_SHORTWAVE Rad=srflx_avg(i,j)*RadScale*EXP(AttW(ng)*z_r(i,j,k)) # else Rad=srflx(i,j)*RadScale*EXP(AttW(ng)*z_r(i,j,k)) # endif IF (z_r(i,j,k).gt.C_depth) THEN cff=Gmax(ng)*G_fac+G_r(ng) G_light=MAX(0.0_r8,cff*TANH(G_eff(ng)*Rad/cff)-G_r(ng)) ELSE G_light=0.0_r8 END IF ! ! Compute dissolved inorganic nutrient (DIN) dependency. ! [JWilkin: This ELSE block below appears redundant because if ! z_r(i,j,k).le.C_depth then G_light=0.0_r8 (see above) and ! therefore G_rate will be set to zero (see below) regardless of ! the calculated value of G_DIN]. ! IF (z_r(i,j,k).gt.C_depth) THEN DIN=DIN_OBS(i,j,k) ELSE DIN=DIN_Cdepth(ng) END IF ! ! The nutrient dependence is modeled by the Monod formulation with ! half-saturation Kn. ! G_DIN=Gmax(ng)*G_fac*DIN/(MAX(Kn(ng),0.0_r8)+DIN) ! ! Compute growth term (implicit). The growth rate is either limited ! by the nutrient or light. The rate is capped to be positive. ! G_rate=MAX(MIN(G_light,G_DIN),0.0_r8) Bio(i,k,iDino)=Bio(i,k,iDino)/(1.0_r8-G_rate*dtdays) END DO END DO ! !----------------------------------------------------------------------- ! Compute mortality term. !----------------------------------------------------------------------- ! ! The mortality is modeled as function dependent on temperature ! (implicit). Use a Q10 mortality rate equation. ! DO k=1,N(ng) DO i=Istr,Iend temp=Bio(i,k,itemp) M_rate=Mor_a(ng)* & & Mor_Q10(ng)**((temp-Mor_T0(ng))*0.1_r8)+ & & Mor_b(ng) Bio(i,k,iDino)=Bio(i,k,iDino)/(1.0_r8+M_rate*dtdays) END DO END DO ! !----------------------------------------------------------------------- ! Vertical sinking/ascending term: dinoflagellate swimming !----------------------------------------------------------------------- ! ! Reconstruct vertical profile of selected biological constituents ! "Bio(:,:,iswim)" in terms of a set of parabolic segments within each ! grid box. Then, compute semi-Lagrangian flux due to vertical motion. ! Many thanks to Sasha Shchepetkin for the updated algorithm. ! SWIM_LOOP: DO iswim=1,Nswim ibio=idswim(iswim) DO k=N(ng)-1,1,-1 DO i=Istr,Iend FC(i,k)=(Bio(i,k+1,ibio)-Bio(i,k,ibio))/ & & (Hz(i,j,k+1)+Hz(i,j,k)) END DO END DO ! DO k=2,N(ng)-1 DO i=Istr,Iend deltaR=Hz(i,j,k)*FC(i,k ) deltaL=Hz(i,j,k)*FC(i,k-1) IF (deltaR*deltaL.lt.0.0_r8) THEN deltaR=0.0_r8 deltaL=0.0_r8 END IF cff=Hz(i,j,k-1)+2.0_r8*Hz(i,j,k)+Hz(i,j,k+1) cffR=cff*FC(i,k ) cffL=cff*FC(i,k-1) IF (ABS(deltaR).gt.ABS(cffL)) deltaR=cffL IF (ABS(deltaL).gt.ABS(cffR)) deltaL=cffR cff=(deltaR-deltaL)/(Hz(i,j,k-1)+Hz(i,j,k)+Hz(i,j,k+1)) deltaR=deltaR-cff*Hz(i,j,k+1) deltaL=deltaL+cff*Hz(i,j,k-1) ! aR(i,k)=Bio(i,k,ibio)+deltaR aL(i,k)=Bio(i,k,ibio)-deltaL ! dR(i,k)=(2.0_r8*deltaR-deltaL)**2 dL(i,k)=(2.0_r8*deltaL-deltaR)**2 END DO END DO #ifdef LIMIT_INTERIOR ! ! Apply boundary conditions for strictly monotonic option. The only way ! to avoid extrapolation toward the boundary is to assume that field is ! simply constant within topmost and bottommost grid boxes. ! DO i=Istr,Iend aR(i,N(ng))=Bio(i,N(ng),ibio) aL(i,N(ng))=Bio(i,N(ng),ibio) dR(i,N(ng))=0.0_r8 dL(i,N(ng))=0.0_r8 ! aR(i,1)=Bio(i,1,ibio) aL(i,1)=Bio(i,1,ibio) dR(i,1)=0.0_r8 dL(i,1)=0.0_r8 END DO #else ! ! Apply Neumann or linear continuation boundary conditions. Notice that ! for Neumann conditions, the extrapolate values aR(i,N(ng)) and aL(i,0) ! exceed corresponding box values. ! DO i=Istr,Iend aL(i,N(ng))=aR(i,N(ng)-1) # ifdef NEUMANN aR(i,N(ng))=1.5_r8*Bio(i,N(ng),ibio)-0.5_r8*aL(i,N(ng)) # else aR(i,N(ng))=2.0_r8*Bio(i,N(ng),ibio)-aL(i,N(ng)) # endif dR(i,N(ng))=(2.0_r8*aR(i,N(ng))+aL(i,N(ng))- & & 3.0_r8*Bio(i,N(ng),ibio))**2 dL(i,N(ng))=(3.0_r8*Bio(i,N(ng),ibio)- & & 2.0_r8*aL(i,N(ng))-aR(i,N(ng)))**2 ! aR(i,1)=aL(i,2) # ifdef NEUMANN aL(i,1)=1.5_r8*Bio(i,1,ibio)-0.5_r8*aR(i,1) # else aL(i,1)=2.0_r8*Bio(i,1,ibio)-aR(i,1) # endif dR(i,1)=(2.0_r8*aR(i,1)+aL(i,1)- & & 3.0_r8*Bio(i,1,ibio))**2 dL(i,1)=(3.0_r8*Bio(i,1,ibio)- & & 2.0_r8*aL(i,1)-aR(i,1))**2 END DO #endif ! ! Reconcile interfacial values aR and aL using Weighted Essentially ! Non-Oscillatory (WENO) procedure. ! DO k=1,N(ng)-1 DO i=Istr,Iend deltaL=MAX(dL(i,k ),eps) deltaR=MAX(dR(i,k+1),eps) r(i,k)=(deltaR*aR(i,k)+deltaL*aL(i,k+1))/ & & (deltaR+deltaL) END DO END DO DO i=Istr,Iend #ifdef NEUMANN r(i,N(ng))=1.5_r8*Bio(i,N(ng),ibio)-0.5_r8*r(i,N(ng)-1) r(i,0 )=1.5_r8*Bio(i,1 ,ibio)-0.5_r8*r(i,1 ) #else r(i,N(ng))=2.0_r8*Bio(i,N(ng),ibio)-r(i,N(ng)-1) r(i,0 )=2.0_r8*Bio(i,1 ,ibio)-r(i,1 ) #endif END DO ! ! Remapping step: This operation consists essentially of three stages: !---------------- (1) within each grid box compute averaged slope ! (stored as dR) and curvature (stored as dL); then (2) compute ! interfacial fluxes FC; and (3) apply these fluxes to complete ! remapping step. ! DO k=1,N(ng) DO i=Istr,Iend #ifdef LIMIT_INTERIOR deltaR=r(i,k)-Bio(i,k,ibio) ! Constrain parabolic deltaL=Bio(i,k,ibio)-r(i,k-1) ! segment monotonicity cffR=2.0_r8*deltaR ! like in PPM cffL=2.0_r8*deltaL IF (deltaR*deltaL.lt.0.0_r8) THEN deltaR=0.0_r8 deltaL=0.0_r8 ELSE IF (ABS(deltaR).gt.ABS(cffL)) THEN deltaR=cffL ELSE IF (ABS(deltaL).gt.ABS(cffR)) THEN deltaL=cffR END IF aR(i,k)=Bio(i,k,ibio)+deltaR aL(i,k)=Bio(i,k,ibio)-deltaL #else aR(i,k)=r(i,k ) aL(i,k)=r(i,k-1) #endif dL(i,k)=0.5_r8*(aR(i,k)-aL(i,k)) dR(i,k)=0.5_r8*(aR(i,k)+aL(i,k))-Bio(i,k,ibio) END DO END DO ! ! Compute interfacial fluxes. The convention is that Wbio is positive ! for upward motion (swimming, floating) and negative for downward motion ! (sinking). ! wdt=-Wbio(iswim)*dtdays DO k=1,N(ng)-1 DO i=Istr,Iend IF (wdt.gt.0.0_r8) THEN ! downward vertical alpha=Hz(i,j,k+1) ! motion (sinking) cff =aL(i,k+1) cffL=dL(i,k+1) cffR=dR(i,k+1) dz=wdt ELSE ! upward vertical alpha=-Hz(i,j,k) ! motion (swimming, cff =aR(i,k) ! migration) cffL=-dL(i,k) cffR=dR(i,k) dz=wdt !! IF (ABS(z_w(i,j,k)).lt.21.0_r8) THEN !! dz=wdt*(1.0_r8-TANH((21.0_r8+z_w(i,j,k))*0.1_r8)) !! ELSE !! dz=wdt !! END IF END IF alpha=dz/alpha ! Courant number FC(i,k)=dz*(cff+alpha*(cffL-cffR*(3.0_r8-2.0_r8*alpha))) END DO END DO DO i=Istr,Iend FC(i,0 )=0.0_r8 FC(i,N(ng))=0.0_r8 END DO ! ! Add semi-Lagrangian vertical flux. ! DO k=1,N(ng) DO i=Istr,Iend cff=(FC(i,k)-FC(i,k-1))*Hz_inv(i,k) Bio(i,k,ibio)=Bio(i,k,ibio)+cff END DO END DO END DO SWIM_LOOP END DO ITER_LOOP ! !----------------------------------------------------------------------- ! Update global tracer variables: Add increment due to BGC processes ! to tracer array in time index "nnew". Index "nnew" is solution after ! advection and mixing and has transport units (m Tunits) hence the ! increment is multiplied by Hz. Notice that we need to subtract ! original values "Bio_old" at the top of the routine to just account ! for the concentrations affected by BGC processes. This also takes ! into account any constraints (non-negative concentrations, carbon ! concentration range) specified before entering BGC kernel. If "Bio" ! were unchanged by BGC processes, the increment would be exactly ! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not ! bounded >=0 so that we can preserve total inventory of nutrients ! when advection causes tracer concentration to go negative. !----------------------------------------------------------------------- ! DO itrc=1,NBT ibio=idbio(itrc) DO k=1,N(ng) DO i=Istr,Iend cff=Bio(i,k,ibio)-Bio_old(i,k,ibio) t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*Hz(i,j,k) END DO END DO END DO END DO J_LOOP ! RETURN END SUBROUTINE red_tide_tile END MODULE biology_mod