MODULE biology_mod
!
!git $Id$
!svn $Id: npzd_Franks.h 1151 2023-02-09 03:08:53Z arango $
!================================================== Hernan G. Arango ===
! Copyright (c) 2002-2023 The ROMS/TOMS Group Craig V. Lewis !
! Licensed under a MIT/X style license !
! See License_ROMS.md !
!=======================================================================
! !
! Nutrient-Phytoplankton-Zooplankton-Detritus Model. !
! !
! This routine computes the biological sources and sinks and adds !
! then the global biological fields. !
! !
! Reference: !
! !
! Franks et al, 1986: Behavior of simple plankton model with !
! food-level acclimation by herbivores, Marine Biology, 91, !
! 121-129. !
! !
! Adapted from code written originally by Craig V. Lewis. !
! !
!=======================================================================
!
implicit none
!
PRIVATE
PUBLIC :: biology
!
CONTAINS
!
!***********************************************************************
SUBROUTINE biology (ng,tile)
!***********************************************************************
!
USE mod_param
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 npzd_franks_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, &
& OCEAN(ng) % t)
#ifdef PROFILE
CALL wclock_off (ng, iNLM, 15, __LINE__, MyFile)
#endif
!
RETURN
END SUBROUTINE biology
!
!-----------------------------------------------------------------------
SUBROUTINE npzd_franks_tile (ng, tile, &
& LBi, UBi, LBj, UBj, UBk, UBt, &
& IminS, ImaxS, JminS, JmaxS, &
& nstp, nnew, &
#ifdef MASKING
& rmask, &
#endif
& Hz, z_r, z_w, &
& t)
!-----------------------------------------------------------------------
!
USE mod_param
USE mod_biology
USE mod_ncparam
USE mod_scalars
!
! 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(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(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
#endif
!
! Local variable declarations.
!
integer, parameter :: Nsink = 1
integer :: Iter, i, ibio, isink, itrc, itrmx, j, k, ks
integer, dimension(Nsink) :: idsink
real(r8), parameter :: eps = 1.0e-16_r8
real(r8) :: cff, cff1, cff2, cff3, dtdays
real(r8) :: cffL, cffR, cu, dltL, dltR
real(r8), dimension(Nsink) :: Wbio
integer, dimension(IminS:ImaxS,N(ng)) :: ksource
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,N(ng)) :: Hz_inv
real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
#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 according to the number of iterations.
!
dtdays=dt(ng)*sec2day/REAL(BioIter(ng),r8)
!
! Set vertical sinking indentification vector.
!
idsink(1)=iSDet ! Small detritus
!
! Set vertical sinking velocity vector in the same order as the
! identification vector, IDSINK.
!
Wbio(1)=wDet(ng) ! Small detritus
!
! 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
DO k=1,N(ng)-1
DO i=Istr,Iend
Hz_inv2(i,k)=1.0_r8/(Hz(i,j,k)+Hz(i,j,k+1))
END DO
END DO
DO k=2,N(ng)-1
DO i=Istr,Iend
Hz_inv3(i,k)=1.0_r8/(Hz(i,j,k-1)+Hz(i,j,k)+Hz(i,j,k+1))
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)=t(i,j,k,nstp,ibio)
END DO
END DO
END DO
!
! Determine Correction for negativity.
!
DO k=1,N(ng)
DO i=Istr,Iend
cff1=MAX(0.0_r8,eps-Bio_old(i,k,iNO3_))+ &
& MAX(0.0_r8,eps-Bio_old(i,k,iPhyt))+ &
& MAX(0.0_r8,eps-Bio_old(i,k,iZoop))+ &
& MAX(0.0_r8,eps-Bio_old(i,k,iSDet))
!
! If correction needed, determine the largest pool to debit.
!
IF (cff1.gt.0.0) THEN
itrmx=idbio(1)
cff=t(i,j,k,nstp,itrmx)
DO ibio=idbio(2),idbio(NBT)
IF (t(i,j,k,nstp,ibio).gt.cff) THEN
itrmx=ibio
cff=t(i,j,k,nstp,ibio)
END IF
END DO
!
! Update new values.
!
DO itrc=1,NBT
ibio=idbio(itrc)
Bio(i,k,ibio)=MAX(eps,Bio_old(i,k,ibio))- &
& cff1*(SIGN(0.5_r8, &
& REAL(itrmx-ibio,r8)**2)+ &
& SIGN(0.5_r8, &
& -REAL(itrmx-ibio,r8)**2))
END DO
ELSE
DO itrc=1,NBT
ibio=idbio(itrc)
Bio(i,k,ibio)=Bio_old(i,k,ibio)
END DO
END IF
END DO
END DO
!
!=======================================================================
! Start internal iterations to achieve convergence of the nonlinear
! backward-implicit solution.
!=======================================================================
!
! During the iterative procedure a series of fractional time steps are
! performed in a chained mode (splitting by different biological
! conversion processes) in sequence of the main food chain. In all
! stages the concentration of the component being consumed is treated
! in fully implicit manner, so the algorithm guarantees non-negative
! values, no matter how strong s the concentration of active consuming
! component (Phytoplankton or Zooplankton). The overall algorithm,
! as well as any stage of it, is formulated in conservative form
! (except explicit sinking) in sense that the sum of concentration of
! all components is conserved.
!
! In the implicit algorithm, we have for example (N: nutrient,
! P: phytoplankton),
!
! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
! {Michaelis-Menten}
! below, we set
! The N in the numerator of
! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
! as N(new)
!
! so the time-stepping of the equations becomes:
!
! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
! consuming, divide by (1 + cff)
! and
!
! P(new) = P(old) + cff * N(new) (2) when adding a source term,
! growing, add (cff * source)
!
! Notice that if you substitute (1) in (2), you will get:
!
! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
!
! If you add (1) and (3), you get
!
! N(new) + P(new) = N(old) + P(old)
!
! implying conservation regardless how "cff" is computed. Therefore,
! this scheme is unconditionally stable regardless of the conversion
! rate. It does not generate negative values since the constituent
! to be consumed is always treated implicitly. It is also biased
! toward damping oscillations.
!
! 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)
!
! Nutrient uptake by phytoplankton.
!
cff1=dtdays*Vm_NO3(ng)
DO k=1,N(ng)
DO i=Istr,Iend
cff=Bio(i,k,iPhyt)* &
& cff1*EXP(K_ext(ng)*z_r(i,j,k))/ &
& (K_NO3(ng)+Bio(i,k,iNO3_))
Bio(i,k,iNO3_)=Bio(i,k,iNO3_)/ &
& (1.0_r8+cff)
Bio(i,k,iPhyt)=Bio(i,k,iPhyt)+ &
& Bio(i,k,iNO3_)*cff
END DO
END DO
!
! Phytoplankton grazing by Zooplankton and mortality to Detritus
! (rate: PhyMR).
!
cff1=dtdays*ZooGR(ng)
cff2=dtdays*PhyMR(ng)
cff3=K_phy(ng)*K_phy(ng)
DO k=1,N(ng)
DO i=Istr,Iend
cff=Bio(i,k,iZoop)*Bio(i,k,iPhyt)*cff1/ &
& (cff3+Bio(i,k,iPhyt)*Bio(i,k,iPhyt))
Bio(i,k,iPhyt)=Bio(i,k,iPhyt)/ &
& (1.0_r8+cff+cff2)
Bio(i,k,iZoop)=Bio(i,k,iZoop)+ &
& Bio(i,k,iPhyt)*cff*(1.0_r8-ZooGA(ng))
Bio(i,k,iSDet)=Bio(i,k,iSDet)+ &
& Bio(i,k,iPhyt)* &
& (cff2+cff*(ZooGA(ng)-ZooEC(ng)))
Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
& Bio(i,k,iPhyt)*cff*ZooEC(ng)
END DO
END DO
!
! Zooplankton excretion to nutrients and mortality to Detritus.
!
cff1=1.0_r8/(1.0_r8+dtdays*(ZooMR(ng)+ZooMD(ng)))
cff2=dtdays*ZooMR(ng)
cff3=dtdays*ZooMD(ng)
DO k=1,N(ng)
DO i=Istr,Iend
Bio(i,k,iZoop)=Bio(i,k,iZoop)*cff1
Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
& Bio(i,k,iZoop)*cff2
Bio(i,k,iSDet)=Bio(i,k,iSDet)+ &
& Bio(i,k,iZoop)*cff3
END DO
END DO
!
! Detritus breakdown to nutrients.
!
cff1=dtdays*DetRR(ng)
cff2=1.0_r8/(1.0_r8+cff1)
DO k=1,N(ng)
DO i=Istr,Iend
Bio(i,k,iSDet)=Bio(i,k,iSDet)*cff2
Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
& Bio(i,k,iSDet)*cff1
END DO
END DO
!
!-----------------------------------------------------------------------
! Vertical sinking terms.
!-----------------------------------------------------------------------
!
! Reconstruct vertical profile of selected biological constituents
! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
! grid box. Then, compute semi-Lagrangian flux due to sinking.
!
SINK_LOOP: DO isink=1,Nsink
ibio=idsink(isink)
!
! Copy concentration of biological particulates into scratch array
! "qc" (q-central, restrict it to be positive) which is hereafter
! interpreted as a set of grid-box averaged values for biogeochemical
! constituent concentration.
!
DO k=1,N(ng)
DO i=Istr,Iend
qc(i,k)=Bio(i,k,ibio)
END DO
END DO
!
DO k=N(ng)-1,1,-1
DO i=Istr,Iend
FC(i,k)=(qc(i,k+1)-qc(i,k))*Hz_inv2(i,k)
END DO
END DO
DO k=2,N(ng)-1
DO i=Istr,Iend
dltR=Hz(i,j,k)*FC(i,k)
dltL=Hz(i,j,k)*FC(i,k-1)
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)
!
! Apply PPM monotonicity constraint to prevent oscillations within the
! grid box.
!
IF ((dltR*dltL).le.0.0_r8) THEN
dltR=0.0_r8
dltL=0.0_r8
ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN
dltR=cffL
ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN
dltL=cffR
END IF
!
! Compute right and left side values (bR,bL) of parabolic segments
! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
!
! NOTE: Although each parabolic segment is monotonic within its grid
! box, monotonicity of the whole profile is not guaranteed,
! because bL(k+1)-bR(k) may still have different sign than
! qc(i,k+1)-qc(i,k). This possibility is excluded,
! after bL and bR are reconciled using WENO procedure.
!
cff=(dltR-dltL)*Hz_inv3(i,k)
dltR=dltR-cff*Hz(i,j,k+1)
dltL=dltL+cff*Hz(i,j,k-1)
bR(i,k)=qc(i,k)+dltR
bL(i,k)=qc(i,k)-dltL
WR(i,k)=(2.0_r8*dltR-dltL)**2
WL(i,k)=(dltR-2.0_r8*dltL)**2
END DO
END DO
cff=1.0E-14_r8
DO k=2,N(ng)-2
DO i=Istr,Iend
dltL=MAX(cff,WL(i,k ))
dltR=MAX(cff,WR(i,k+1))
bR(i,k)=(dltR*bR(i,k)+dltL*bL(i,k+1))/(dltR+dltL)
bL(i,k+1)=bR(i,k)
END DO
END DO
DO i=Istr,Iend
FC(i,N(ng))=0.0_r8 ! NO-flux boundary condition
#if defined LINEAR_CONTINUATION
bL(i,N(ng))=bR(i,N(ng)-1)
bR(i,N(ng))=2.0_r8*qc(i,N(ng))-bL(i,N(ng))
#elif defined NEUMANN
bL(i,N(ng))=bR(i,N(ng)-1)
bR(i,N(ng))=1.5_r8*qc(i,N(ng))-0.5_r8*bL(i,N(ng))
#else
bR(i,N(ng))=qc(i,N(ng)) ! default strictly monotonic
bL(i,N(ng))=qc(i,N(ng)) ! conditions
bR(i,N(ng)-1)=qc(i,N(ng))
#endif
#if defined LINEAR_CONTINUATION
bR(i,1)=bL(i,2)
bL(i,1)=2.0_r8*qc(i,1)-bR(i,1)
#elif defined NEUMANN
bR(i,1)=bL(i,2)
bL(i,1)=1.5_r8*qc(i,1)-0.5_r8*bR(i,1)
#else
bL(i,2)=qc(i,1) ! bottom grid boxes are
bR(i,1)=qc(i,1) ! re-assumed to be
bL(i,1)=qc(i,1) ! piecewise constant.
#endif
END DO
!
! Apply monotonicity constraint again, since the reconciled interfacial
! values may cause a non-monotonic behavior of the parabolic segments
! inside the grid box.
!
DO k=1,N(ng)
DO i=Istr,Iend
dltR=bR(i,k)-qc(i,k)
dltL=qc(i,k)-bL(i,k)
cffR=2.0_r8*dltR
cffL=2.0_r8*dltL
IF ((dltR*dltL).lt.0.0_r8) THEN
dltR=0.0_r8
dltL=0.0_r8
ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN
dltR=cffL
ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN
dltL=cffR
END IF
bR(i,k)=qc(i,k)+dltR
bL(i,k)=qc(i,k)-dltL
END DO
END DO
!
! After this moment reconstruction is considered complete. The next
! stage is to compute vertical advective fluxes, FC. It is expected
! that sinking may occurs relatively fast, the algorithm is designed
! to be free of CFL criterion, which is achieved by allowing
! integration bounds for semi-Lagrangian advective flux to use as
! many grid boxes in upstream direction as necessary.
!
! In the two code segments below, WL is the z-coordinate of the
! departure point for grid box interface z_w with the same indices;
! FC is the finite volume flux; ksource(:,k) is index of vertical
! grid box which contains the departure point (restricted by N(ng)).
! During the search: also add in content of whole grid boxes
! participating in FC.
!
cff=dtdays*ABS(Wbio(isink))
DO k=1,N(ng)
DO i=Istr,Iend
FC(i,k-1)=0.0_r8
WL(i,k)=z_w(i,j,k-1)+cff
WR(i,k)=Hz(i,j,k)*qc(i,k)
ksource(i,k)=k
END DO
END DO
DO k=1,N(ng)
DO ks=k,N(ng)-1
DO i=Istr,Iend
IF (WL(i,k).gt.z_w(i,j,ks)) THEN
ksource(i,k)=ks+1
FC(i,k-1)=FC(i,k-1)+WR(i,ks)
END IF
END DO
END DO
END DO
!
! Finalize computation of flux: add fractional part.
!
DO k=1,N(ng)
DO i=Istr,Iend
ks=ksource(i,k)
cu=MIN(1.0_r8,(WL(i,k)-z_w(i,j,ks-1))*Hz_inv(i,ks))
FC(i,k-1)=FC(i,k-1)+ &
& Hz(i,j,ks)*cu* &
& (bL(i,ks)+ &
& cu*(0.5_r8*(bR(i,ks)-bL(i,ks))- &
& (1.5_r8-cu)* &
& (bR(i,ks)+bL(i,ks)- &
& 2.0_r8*qc(i,ks))))
END DO
END DO
DO k=1,N(ng)
DO i=Istr,Iend
Bio(i,k,ibio)=qc(i,k)+(FC(i,k)-FC(i,k-1))*Hz_inv(i,k)
END DO
END DO
END DO SINK_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 concentractions 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 npzd_franks_tile
END MODULE biology_mod