npzd_iron.h

Go to the documentation of this file.

      MODULE biology_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Jerome Fiechter   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  Nutrient-Phytoplankton-Zooplankton-Detritus Model,                  !
!  including Iron Limitation on Phytoplankton Growth.                  !
!                                                                      !
!  This routine computes the biological sources and sinks and adds     !
!  then the global biological fields.                                  !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Fiechter, J., A.M. Moore, C.A. Edwards, K.W. Bruland,             !
!      E. Di Lorenzo, C.V.W. Lewis, T.M. Powell, E. Curchitser         !
!      and K. Hedstrom, 2009: Modeling iron limitation of primary      !
!      production in the coastal Gulf of Alaska, Deep Sea Res. II      !
!                                                                      !
!=======================================================================
!
      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 npzd_iron_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
#ifdef IRON_LIMIT
     &                     grid(ng) % h,                                &
#endif
     &                     grid(ng) % Hz,                               &
     &                     grid(ng) % z_r,                              &
     &                     grid(ng) % z_w,                              &
     &                     forces(ng) % srflx,                          &
     &                     ocean(ng) % t)
 
#ifdef PROFILE
      CALL wclock_off (ng, inlm, 15, __line__, myfile)
#endif
!
      RETURN
      END SUBROUTINE biology
!
!-----------------------------------------------------------------------
      SUBROUTINE npzd_iron_tile (ng, tile,                              &
     &                           LBi, UBi, LBj, UBj, UBk, UBt,          &
     &                           IminS, ImaxS, JminS, JmaxS,            &
     &                           nstp, nnew,                            &
#ifdef MASKING
     &                           rmask,                                 &
#endif
#ifdef IRON_LIMIT
     &                           h,                                     &
#endif
     &                           Hz, z_r, z_w,                          &
     &                           srflx,                                 &
     &                           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
# ifdef IRON_LIMIT
      real(r8), intent(in) :: h(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) :: srflx(LBi:,LBj:)
      real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
#else
# ifdef MASKING
      real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
# ifdef IRON_LIMIT
      real(r8), intent(in) :: h(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) :: srflx(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
#endif
!
!  Local variable declarations.
!
      integer, parameter :: Nsink = 2
 
      integer :: Iter, i, ibio, isink, itime, itrc, iTrcMax, j, k, ks
 
      integer, dimension(Nsink) :: idsink
 
      real(r8), parameter :: MinVal = 1.0e-6_r8
 
      real(r8) :: Att, ExpAtt, Itop, PAR
      real(r8) :: cff, cff1, cff2, cff3, cff4, cff5, cff6, dtdays
      real(r8) :: cffL, cffR, cu, dltL, dltR
      real(r8) :: fac
#ifdef IRON_LIMIT
      real(r8) :: Nlimit, FNlim
      real(r8) :: FNratio, FCratio, FCratioE, Flimit
      real(r8) :: FeC2FeN, FeN2FeC
      real(r8) :: cffFe
# ifdef IRON_RELAX
      real(r8) :: FeNudgCoef
# endif
#endif
      real(r8), dimension(Nsink) :: Wbio
 
      integer, dimension(IminS:ImaxS,N(ng)) :: ksource
 
      real(r8), dimension(IminS:ImaxS) :: PARsur
 
      real(r8), dimension(NT(ng),2) :: BioTrc
 
      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)) :: Light
      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 size (days) according to the number of iterations.
!
      dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
 
#if defined IRON_LIMIT && defined IRON_RELAX
!
!  Set nudging coefficient for dissolved iron over the shelf.
!
      fenudgcoef=dt(ng)/(fenudgtime(ng)*86400.0_r8)
#endif
#ifdef IRON_LIMIT
!
!  Set Fe:N and Fe:C conversion ratio and its inverse.
!
          fen2fec=(16.0_r8/106.0_r8)*1.0e3_r8
          fec2fen=(106.0_r8/16.0_r8)*1.0e-3_r8
#endif
!
!  Set vertical sinking indentification vector.
!
      idsink(1)=iphyt                 ! Phytoplankton
      idsink(2)=isdet                 ! Small detritus
!
!  Set vertical sinking velocity vector in the same order as the
!  identification vector, IDSINK.
!
      wbio(1)=wphy(ng)                ! Phytoplankton
      wbio(2)=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
!
!  Restrict biological tracer to be positive definite. If a negative
!  concentration is detected, nitrogen is drawn from the most abundant
!  pool to supplement the negative pools to a lower limit of MinVal
!  which is set to 1E-6 above.
!
        DO k=1,n(ng)
          DO i=istr,iend
!
!  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)
              biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
              biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)*hz_inv(i,k)
            END DO
!
!  Impose positive definite concentrations.
!
            cff2=0.0_r8
            DO itime=1,2
              cff1=0.0_r8
              itrcmax=idbio(1)
#ifdef IRON_LIMIT
              DO itrc=1,nbt-2
#else
              DO itrc=1,nbt
#endif
                ibio=idbio(itrc)
                cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
                IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
                  itrcmax=ibio
                END IF
                biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
              END DO
              IF (biotrc(itrcmax,itime).gt.cff1) THEN
                biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
              END IF
#ifdef IRON_LIMIT
              DO itrc=nbt-1,nbt
                ibio=idbio(itrc)
                biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
              END DO
#endif
            END DO
!
!  Load biological tracers into local arrays.
!
            DO itrc=1,nbt
              ibio=idbio(itrc)
              bio_old(i,k,ibio)=biotrc(ibio,nstp)
              bio(i,k,ibio)=biotrc(ibio,nstp)
            END DO
 
#if defined IRON_LIMIT && defined IRON_RELAX
!
!  Relax dissolved iron at coast (h <= FeHim) to a constant value
!  (FeMax) over a time scale (FeNudgTime; days) to simulate sources
!  at the shelf.
!
            IF (h(i,j).le.fehmin(ng)) THEN
              bio(i,k,ifdis)=bio(i,k,ifdis)+                            &
     &                       fenudgcoef*(femax(ng)-bio(i,k,ifdis))
            END IF
#endif
          END DO
        END DO
!
!  Calculate surface Photosynthetically Available Radiation (PAR).  The
!  net shortwave radiation is scaled back to Watts/m2 and multiplied by
!  the fraction that is photosynthetically available, PARfrac.
!
        DO i=istr,iend
#ifdef CONST_PAR
!
!  Specify constant surface irradiance a la Powell and Spitz.
!
          parsur(i)=158.075_r8
#else
          parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
#endif
        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 a fully implicit manner, so the algorithm guarantees non-negative
!  values, no matter how strong 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)
!
!  Compute light attenuation as function of depth.
!
          DO i=istr,iend
            par=parsur(i)
            IF (parsur(i).gt.0.0_r8) THEN              ! day time
              DO k=n(ng),1,-1
!
!  Compute average light attenuation for each grid cell. Here, AttSW is
!  the light attenuation due to seawater and AttPhy is the attenuation
!  due to phytoplankton (self-shading coefficient).
!
                att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))*              &
     &              (z_w(i,j,k)-z_w(i,j,k-1))
                expatt=exp(-att)
                itop=par
                par=itop*(1.0_r8-expatt)/att    ! average at cell center
                light(i,k)=par
!
!  Light attenuation at the bottom of the grid cell. It is the starting
!  PAR value for the next (deeper) vertical grid cell.
!
                par=itop*expatt
              END DO
            ELSE                                       ! night time
              DO k=1,n(ng)
                light(i,k)=0.0_r8
              END DO
            END IF
          END DO
!
!  Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
!  The Michaelis-Menten curve is used to describe the change in uptake
!  rate as a function of nitrate concentration. Here, PhyIS is the
!  initial slope of the P-I curve and K_NO3 is the half saturation of
!  phytoplankton nitrate uptake.
#ifdef IRON_LIMIT
!
!  Growth reduction factors due to iron limitation:
!
!    FNratio     current Fe:N ratio [umol-Fe/mmol-N]
!    FCratio     current Fe:C ratio [umol-Fe/mol-C]
!                  (umol-Fe/mmol-N)*(16 M-N/106 M-C)*(1E3 mmol-C/mol-C)
!    FCratioE    empirical  Fe:C ratio
!    Flimit      Phytoplankton growth reduction factor due to Fe
!                  limitation based on Fe:C ratio
!
#endif
!
          cff1=dtdays*vm_no3(ng)*phyis(ng)
          cff2=vm_no3(ng)*vm_no3(ng)
          cff3=phyis(ng)*phyis(ng)
          DO k=1,n(ng)
            DO i=istr,iend
#ifdef IRON_LIMIT
!
!  Calculate growth reduction factor due to iron limitation.
!
              fnratio=bio(i,k,ifphy)/max(minval,bio(i,k,iphyt))
              fcratio=fnratio*fen2fec
              fcratioe=b_fe(ng)*bio(i,k,ifdis)**a_fe(ng)
              flimit=fcratio*fcratio/                                   &
     &               (fcratio*fcratio+k_fec(ng)*k_fec(ng))
!
              nlimit=1.0_r8/(k_no3(ng)+bio(i,k,ino3_))
              fnlim=min(1.0_r8,flimit/(bio(i,k,ino3_)*nlimit))
#endif
              cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
              cff=bio(i,k,iphyt)*                                       &
#ifdef IRON_LIMIT
     &            cff1*cff4*light(i,k)*fnlim*nlimit
#else
     &            cff1*cff4*light(i,k)/                                 &
     &            (k_no3(ng)+bio(i,k,ino3_))
#endif
              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
 
#ifdef IRON_LIMIT
!
!  Iron uptake proportional to growth.
!
              fac=cff*bio(i,k,ino3_)*fnratio/max(minval,bio(i,k,ifdis))
              bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+fac)
              bio(i,k,ifphy)=bio(i,k,ifphy)+                            &
     &                       bio(i,k,ifdis)*fac
!
!  Iron uptake to reach appropriate Fe:C ratio.
!
              cff5=dtdays*(fcratioe-fcratio)/t_fe(ng)
              cff6=bio(i,k,iphyt)*cff5*fec2fen
              IF (cff6.ge.0.0_r8) THEN
                cff=cff6/max(minval,bio(i,k,ifdis))
                bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+cff)
                bio(i,k,ifphy)=bio(i,k,ifphy)+                          &
     &                         bio(i,k,ifdis)*cff
              ELSE
                cff=-cff6/max(minval,bio(i,k,ifphy))
                bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
                bio(i,k,ifdis)=bio(i,k,ifdis)+                          &
     &                         bio(i,k,ifphy)*cff
              END IF
#endif
            END DO
          END DO
!
!  Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
!  formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
!  pool as function of "sloppy feeding" and metabolic processes
!  (ZooEEN and ZooEED fractions).
#ifdef IRON_LIMIT
!  The lost of phytoplankton to the dissolve iron pool is scale by the
!  remineralization rate (FeRR).
#endif
!
          cff1=dtdays*zoogr(ng)
          cff2=1.0_r8-zooeen(ng)-zooeed(ng)
          DO k=1,n(ng)
            DO i=istr,iend
              cff=bio(i,k,izoop)*                                       &
     &            cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/         &
     &            bio(i,k,iphyt)
              bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
              bio(i,k,izoop)=bio(i,k,izoop)+                            &
     &                       bio(i,k,iphyt)*cff2*cff
              bio(i,k,ino3_)=bio(i,k,ino3_)+                            &
     &                       bio(i,k,iphyt)*zooeen(ng)*cff
              bio(i,k,isdet)=bio(i,k,isdet)+                            &
     &                       bio(i,k,iphyt)*zooeed(ng)*cff
#ifdef IRON_LIMIT
              bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
              bio(i,k,ifdis)=bio(i,k,ifdis)+                            &
     &                       bio(i,k,ifphy)*cff*ferr(ng)
#endif
            END DO
          END DO
!
!  Phytoplankton mortality to nutrients (PhyMRNro rate), detritus
!  (PhyMRD rate), and if applicable dissolved iron (FeRR rate).
!
          cff3=dtdays*phymrd(ng)
          cff2=dtdays*phymrn(ng)
          cff1=1.0_r8/(1.0_r8+cff2+cff3)
          DO k=1,n(ng)
            DO i=istr,iend
              bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
              bio(i,k,ino3_)=bio(i,k,ino3_)+                            &
     &                       bio(i,k,iphyt)*cff2
              bio(i,k,isdet)=bio(i,k,isdet)+                            &
     &                       bio(i,k,iphyt)*cff3
#ifdef IRON_LIMIT
              bio(i,k,ifphy)=bio(i,k,ifphy)*cff1
              bio(i,k,ifdis)=bio(i,k,ifdis)+                            &
     &                       bio(i,k,ifphy)*(cff2+cff3)*ferr(ng)
#endif
            END DO
          END DO
!
!  Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
!  (ZooMRD rate).
!
          cff3=dtdays*zoomrd(ng)
          cff2=dtdays*zoomrn(ng)
          cff1=1.0_r8/(1.0_r8+cff2+cff3)
          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: remineralization (DetRR rate).
!
          cff2=dtdays*detrr(ng)
          cff1=1.0_r8/(1.0_r8+cff2)
          DO k=1,n(ng)
            DO i=istr,iend
              bio(i,k,isdet)=bio(i,k,isdet)*cff1
              bio(i,k,ino3_)=bio(i,k,ino3_)+                            &
     &                       bio(i,k,isdet)*cff2
            END DO
          END DO
!
!-----------------------------------------------------------------------
!  Vertical sinking terms: Phytoplankton and Detritus
!-----------------------------------------------------------------------
!
!  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_iron_tile
 
      END MODULE biology_mod