red_tide.h

Go to the documentation of this file.

#undef  NEWMAN
#define LIMIT_INTERIOR
 
      MODULE biology_mod
!
!git $Id$
!======================================================== Ruoying He ===
!  Copyright (c) 2002-2025 The ROMS 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