#include "cppdefs.h"
      SUBROUTINE biology (ng,tile)
!
!========================================== Alexander F. Shchepetkin ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group                         !
!================================================== Hernan G. Arango ===
!                                                                      !
!  This routine computes the biological sources and sinks and adds     !
!  then the global biological fields.                                  !
!                                                                      !
!  Sarah Hinckley''s GOANPZ Code                                       !
!  Implemented by Craig Lewis (CVL)                                    !
!  Modified by Liz Dobbins and Sarah Hinckley                          !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_biology
      USE mod_forces
      USE mod_grid
      USE mod_ncparam
      USE mod_ocean
      USE mod_stepping
      USE mod_ice

#if defined CLIM_ICE_1D
      USE mod_clima
#endif

      integer, intent(in) :: ng, tile

#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)=__FILE__
      END IF
!
#ifdef PROFILE
      CALL wclock_on (ng, iNLM, 15, __LINE__, __FILE__)
#endif
      CALL biology_tile (ng, tile,                                      &
     &                   LBi, UBi, LBj, UBj,                            &
     &                   IminS, ImaxS, JminS, JmaxS,                    &
     &                   N(ng), NT(ng),                                 &
     &                   nnew(ng), nstp(ng),                            &
#ifdef MASKING
     &                   GRID(ng) % rmask,                              &
#endif
     &                   GRID(ng) % Hz,                                 &
     &                   GRID(ng) % z_r,                                &
     &                   GRID(ng) % z_w,                                &
     &                   FORCES(ng) % srflx,                            &
#if defined BENTHIC
     &                   OCEAN(ng) % bt,                                &

#endif
#if defined ICE_BIO
     &                   OCEAN(ng) % it,                                &
     &                   OCEAN(ng) % itL,                               &
# ifdef CLIM_ICE_1D
     &                   CLIMA(ng) % tclmG,                             &
     &                   CLIMA(ng) % tclm,                              &
# else
     &                   ICE(ng) % ti,                                  &
     &                   ICE(ng) % hi,                                  &
     &                   ICE(ng) % ai,                                  &
     &                   ICE(ng) % ageice,                              &
# endif
#endif
#ifdef STATIONARY
     &                   OCEAN(ng) % st,                                &
     &                   NTS(ng),                                       &
#endif
#ifdef STATIONARY2
     &                   OCEAN(ng) % st2,                               &
     &                   NTS2(ng),                                      &
#endif
#ifdef PROD3
     &                   OCEAN(ng) % pt3,                               &
     &                   NPT3(ng),                                      &
#endif
#ifdef PROD2
     &                   OCEAN(ng) % pt2,                               &
     &                   NPT2(ng),                                      &
#endif
#ifdef BIOFLUX
     &                   OCEAN(ng) % bflx,                              &
#endif
     &                   OCEAN(ng) % t)

#ifdef PROFILE
      CALL wclock_off (ng, iNLM, 15, __LINE__, __FILE__)
#endif
      RETURN
      END SUBROUTINE biology
!
!-----------------------------------------------------------------------
      SUBROUTINE biology_tile (ng, tile,                                &
     &                         LBi, UBi, LBj, UBj,                      &
     &                         IminS, ImaxS, JminS, JmaxS,              &
     &                         UBk, UBt,                                &
     &                         nnew, nstp,                              &
#ifdef MASKING
     &                         rmask,                                   &
#endif
     &                         Hz, z_r, z_w, srflx,                     &
#if defined BENTHIC
     &                         bt,                                      &
#endif
#if defined ICE_BIO
     &                         it,                                      &
     &                         itL,                                     &
# ifdef CLIM_ICE_1D
     &                         tclmG,                                   &
     &                         tclm,                                    &
# else
     &                         ti,                                      &
     &                         hi,                                      &
     &                         ai,                                      &
     &                         ageice,                                  &
# endif
#endif
#ifdef STATIONARY
     &                          st,                                     &
     &                          UBst,                                   &
#endif
#ifdef STATIONARY2
     &                          st2,                                    &
     &                          UBst2,                                  &
#endif
#ifdef PROD3
     &                          pt3,                                    &
     &                          UBpt3,                                  &
#endif
#ifdef PROD2
     &                          pt2,                                    &
     &                          UBpt2,                                  &
#endif
#ifdef BIOFLUX
     &                          bflx,                                   &
#endif
     &                          t)

!-----------------------------------------------------------------------
!
      USE mod_param
      USE mod_biology
      USE mod_scalars
      USE mod_ocean
      USE mod_grid
      USE mod_biology
      USE dateclock_mod,   ONLY : caldate
#if defined CLIM_ICE_1D
      USE mod_clima
#endif

      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: UBk, UBt
      integer, intent(in) :: nnew, nstp
!
#if defined STATIONARY
      integer, intent(in) :: UBst
#endif
#if defined STATIONARY2
      integer, intent(in) ::  UBst2
#  endif
#if defined PROD3
      integer, intent(in) :: UBpt3
#endif
#if defined PROD2
      integer, intent(in) ::  UBpt2
#endif

#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) :: srflx(LBi:,LBj:)
      real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)

# ifdef STATIONARY
      real(r8), intent(inout) :: st(LBi:,LBj:,:,:,:)
# endif
# ifdef STATIONARY2
      real(r8), intent(inout) :: st2(LBi:,LBj:,:,:)
# endif
# ifdef PROD3
      real(r8), intent(inout) :: pt3(LBi:,LBj:,:,:,:)
# endif
# ifdef PROD2
      real(r8), intent(inout) :: pt2(LBi:,LBj:,:,:)
# endif
# if defined BENTHIC
      real(r8), intent(inout) :: bt(LBi:,LBj:,:,:,:)
# endif
# if defined ICE_BIO
      real(r8), intent(inout) :: it(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: itL(LBi:,LBj:,:,:)
#  ifdef CLIM_ICE_1D
      real(r8), intent(inout) ::tclmG(LBi:,LBj:,:,:,:)
      real(r8), intent(inout) ::tclm(LBi:,LBj:,:,:)
#  else
      real(r8), intent(in) :: ti(LBi:,LBj:,:)
      real(r8), intent(in) :: hi(LBi:,LBj:,:)
      real(r8), intent(in) :: ai(LBi:,LBj:,:)
      real(r8), intent(in) :: ageice(LBi:,LBj:,:)
#  endif
#  ifdef BIOFLUX
      real(r8), intent(inout) :: bflx(:,:)
#  endif
# endif

#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) :: srflx(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
# ifdef STATIONARY
      real(r8), intent(inout) :: st(LBi:UBi,LBj:UBj,UBk,3,UBst)
# endif
# ifdef STATIONARY2
      real(r8), intent(inout) :: st2(LBi:UBi,LBj:UBj,3,UBst2)
# endif
# ifdef PROD3
      real(r8), intent(inout) :: pt3(LBi:UBi,LBj:UBj,UBk,3,UBpt3)
# endif
# ifdef PROD2
      real(r8), intent(inout) :: pt2(LBi:UBi,LBj:UBj,3,UBpt2)
# endif
# if defined BENTHIC
      real(r8), intent(inout) :: bt(LBi:UBi,LBj:UBj,UBk,3,1)
# endif
# if defined ICE_BIO
      real(r8), intent(inout) :: it(LBi:UBi,LBj:UBj,3,1)
      real(r8), intent(inout) :: itL(LBi:UBi,LBj:UBj,3,1)
#  ifdef CLIM_ICE_1D
      real(r8), intent(inout) ::tclmG(LBi:UBi,LBj:UBj,UBk,3,UBt+1)
      real(r8), intent(inout) ::tclm(LBi:UBi,LBj:UBj,UBk,UBt+1)
#  else
      real(r8), intent(in) :: ti(LBi:UBi,LBj:UBj,2)
      real(r8), intent(in) :: hi(LBi:UBi,LBj:UBj,2)
      real(r8), intent(in) :: ai(LBi:UBi,LBj:UBj,2)
      real(r8), intent(in) :: ageice(LBi:UBi,LBj:UBj,2)
#  endif
# endif
# ifdef BIOFLUX
      real(r8), intent(inout) :: bflx(UBt,UBt)
# endif
#endif
!
!  Local variable declarations.
!
      integer :: i, j, k, ibio, ibio2,itr, itrmx, itrc, itrc2
#ifdef BENTHIC
      integer :: ibioB
      real(r8) :: cff5,cff6,cff7,cff8,cff9,cff10
#endif
#ifdef ICE_BIO
      integer :: ibioBI
#endif
      integer :: Iter,is

      real(r8) :: cff1, cff2, cff3,cff4
      real(r8) :: Drate, Pmax, NOup, NHup
      real(r8) :: dtdays
      real(r8) :: LightLim, NOLim, NHLim, IronLim
      real(r8) :: yday, lat, k_phy, Dl, Par1
      real(r8) :: Sal1, Temp1
!      , TmaxPhS, KtBm_PhS, TmaxPhL, KtBm_PhL,TmaxMZS,KtBm_MZS,TmaxMZL,KtBm_MZL
      real(r8) :: ParMax,BasalMet
      real(r8) :: Iron1, kfePh,respPh
!      real(r8) :: PON,Pv0,PvT,Dep1,Nitrif,NH4R
      real(r8) :: PON,Dep1,Nitrif,NH4R
      real(r8) :: NitrifMax,DLNitrif

      real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
      real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: DBio
      real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_bak
      real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Prod
      real(r8), dimension(IminS:ImaxS,NT(ng)) :: Prod2
#if defined STATIONARY

       real(r8), dimension(IminS:ImaxS,N(ng),NTS(ng)) :: Stat3
#endif
#if defined STATIONARY2
       real(r8), dimension(IminS:ImaxS,NTS(ng)) :: Stat2
#endif

#if defined BENTHIC
      real(r8), dimension(IminS:ImaxS,NBL(ng),NBeT(ng)) :: BioB
      real(r8), dimension(IminS:ImaxS,NBL(ng),NBeT(ng)) :: DBioB
      real(r8), dimension(IminS:ImaxS,NBL(ng),NBeT(ng)) :: Bio_bakB
#endif
#if defined ICE_BIO
      real(r8), dimension(IminS:ImaxS,NIceT(ng)) :: BioBI
      real(r8), dimension(IminS:ImaxS,NIceT(ng)) :: DBioBI
      real(r8), dimension(IminS:ImaxS,NIceT(ng)) :: Bio_bakBI
#endif
#if defined BIOFLUX
      real(r8), dimension(NT(ng),NT(ng)) :: BioFlx
#endif
      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) :: PARsur
      real(r8), dimension(IminS:ImaxS,N(ng)) :: PAR
      real(r8), dimension(IminS:ImaxS,N(ng)) :: Dens
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TestVal
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TempFuncPhS
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TempFuncPhL
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TempFuncMZS
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TempFuncMZL
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TempFuncCop
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TempFuncNeo
      real(r8), dimension(IminS:ImaxS,N(ng)) :: TempFuncEup
#ifdef JELLY
      real(r8), dimension(PRIVATE_1D_SCRATCH_ARRAY,N(ng))::TempFuncJel
#endif
      real(r8), dimension(IminS:ImaxS,N(ng)) :: HzL
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: z_wL

      real(r8), dimension(IminS:ImaxS,N(ng)) :: sinkIN,sinkOUT,riseIN,riseOUT
#ifdef ICE_BIO
      real(r8) :: aiceIfrac,aiceNfrac,dhicedt,trs,cwi,twi
      real(r8) ::grow1, GROWAice,reN,fNO3,RAi0,RgAi
!      real(r8) :: alpha=0.8_r8, beta=0.018_r8,  inhib=1.46_r8
!      real(r8) :: ksnut1=1.0,ksnut2=4.0,mu0=1.44, R0i=0.05
!      real(r8) :: rg0=9.23e-4,rg=0.03,annit=6.2e-4,aidz=0.02
      real(r8) :: sb, gesi
      real(r8), dimension(PRIVATE_1D_SCRATCH_ARRAY,N(ng),3):: aib
!      logical :: IceBioInit=.FALSE.
      real(r8), dimension(PRIVATE_2D_SCRATCH_ARRAY) :: ice_thick
#endif

#ifdef DIAPAUSE
      logical :: downward = .false., upward = .false.
#endif
!
      real(r8), parameter :: eps  = 1.0E-20_r8
      real(r8), parameter :: minv = 0.0E-20_r8
!----------------------
#include "set_bounds.h"
!
      CALL caldate(tdays(ng), yd_dp=yday)
      dtdays = dt(ng)*sec2day/REAL(BioIter(ng),r8)
      k_phy = k_chl / ccr
!
      IF (yday.eq.34)THEN
        print*,'DAY =155'
      END IF

#ifdef DIAPAUSE
!  Based on date, determine if NCa are going upwards or downwards
      downward = .false.
      upward = .false.
      IF ( ( RiseStart.lt.RiseEnd .and.                                 &
     &       yday.ge.RiseStart .and. yday.le.RiseEnd ) .or.             &
     &     ( RiseStart.gt.RiseEnd .and.                                 &
     &      ( yday.ge.RiseStart .or. yday.le.RiseEnd ) ) )  THEN
        upward = .true.
      ELSE IF ( ( SinkStart.lt.SinkEnd .and.                            &
     &       yday.ge.SinkStart .and. yday.le.SinkEnd ) .or.             &
     &     ( SinkStart.gt.SinkEnd .and.                                 &
     &      ( yday.ge.SinkStart .or. yday.le.SinkEnd ) ) )  THEN
        downward = .true.
      END IF
#endif

!
! ----------------------------------------------------------------------
! Begin HORIZONTAL INDEX LOOPING
! ----------------------------------------------------------------------
!
     J_LOOP : DO j=Jstr,Jend
#if !defined ANA_BIOLOGY && defined BIOLOGY && defined BIO_GOANPZ
# ifdef BENTHIC
!if this is the first time step
        IF (iic(ng).eq.ntstart(ng)) THEN
          DO k=1,NBL(ng)
            DO i=Istr,Iend
              bt(i,j,k,1,1) = 8000_r8
              bt(i,j,k,2,1) = 8000_r8
              bt(i,j,k,3,1) = 8000_r8
              bt(i,j,k,1,2) = 000_r8
              bt(i,j,k,2,2) = 000_r8
              bt(i,j,k,3,2) = 000_r8

!              bt(i,j,k,2,itrc) = bt(i,j,k,1,itrc)
!              bt(i,j,k,3,itrc) = bt(i,j,k,1,itrc)
            END DO
          END DO
        END IF

# endif
#endif
        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

#ifdef ICE_BIO
        DO i=Istr,Iend

# if defined CLIM_ICE_1D
          ice_thick(i,j) = tclmG(i,j,42,1,15)

!         ice_thick(i,j) =   it(i,j,nnew,iIceZ)
# else

          ice_thick(i,j) = hi(i,j,nstp)/                                &
     &                    ai(i,j,nstp)
# endif

!g         IF (hi(i,j,nstp).eq.0.002.and.ai(i,j,nstp).eq.0.002)THEN

!g           ice_thick(i,j) =0_r8
!g         END IF
        END DO
#endif
!
!  Extract biological variables from tracer arrays, and
!  restrict their values to be positive definite.  Removed CVL''s
!  conservation of mass correction because conflicted with SPLINES.
!  For ROMS 2.2+, convert from "flux form" to concentrations by
!  dividing by grid cell thickness.
        DO itrc=1,NBT
          ibio=idbio(itrc)
          DO k=1,N(ng)
            DO i=Istr,Iend
              Bio_bak(i,k,ibio)=max(t(i,j,k,nstp,ibio),0.0_r8)
              Bio(i,k,ibio)=Bio_bak(i,k,ibio)
#ifdef PROD3
              Prod(i,k,ibio)=0.0_r8
#endif

              DBio(i,k,ibio)=0.0_r8
            END DO
          END DO
        END DO
#ifdef PROD2
        DO itrc=1,NBPT2
!          ibio=idbioP2(itrc)
          DO i=Istr,Iend
            Prod2(i,itrc)=0.0_r8
          END DO
        END DO
#endif

#ifdef STATIONARY
        DO itrc=1,UBst
          ibio=idbio3(itrc)
          DO i=Istr,Iend
!           Stat3(i,k,ibio)=st(i,j,k,nstp,ibio)
            Stat3(i,k,ibio)=0.0_r8
          END DO
        END DO
#endif
#ifdef STATIONARY2
        DO itrc=1,UBst2
          ibio=idbio2(itrc)
          DO i=Istr,Iend
            Stat2(i,ibio)=0.0_r8
          END DO
        END DO
#endif
#ifdef BENTHIC
        DO itrc=1,NBEN
          ibioB=idben(itrc)
          DO k=1,NBL(ng)
            DO i=Istr,Iend
              Bio_bakB(i,k,ibioB)=max(bt(i,j,k,nstp,ibioB),0.0_r8)
              BioB(i,k,ibioB)=bt(i,j,k,nstp,ibioB)
              DBioB(i,k,ibioB)=0.0_r8
            END DO
          END DO
        END DO
#endif

#if defined ICE_BIO
        DO i=Istr,Iend
          IF (ice_thick(i,j).ge.0.02)          THEN
            IF (itL(i,j,nstp,iIceLog).le.0 )  THEN
              !                print*,'Ice greater thatn 0.2'
!
! Initialize the ice biology
!
              it(i,j,1,iIcePhL) = 1.1638_r8  !0.1638_r8
              it(i,j,1,iIceNO3) = Bio(i,N(ng),iNO3)  !5.0_r8
              it(i,j,1,iIceNH4) = Bio(i,N(ng),iNH4)  !1.0_r8
              itL(i,j,1,iIceLog) =1.0_r8
!
              it(i,j,2,iIcePhL) = it(i,j,1,iIcePhL)
              it(i,j,2,iIceNO3) = it(i,j,1,iIceNO3)
              it(i,j,2,iIceNH4) = it(i,j,1,iIceNH4)
              itL(i,j,2,iIceLog) =itL(i,j,1,iIceLog)
            END IF
          END IF
        END DO
#endif
#ifdef ICE_BIO
        DO itrc=1,3
          ibioBI=idice(itrc)
          DO i=Istr,Iend
            Bio_bakBI(i,ibioBI)=max(it(i,j,nstp,ibioBI),0.0_r8)
            BioBI(i,ibioBI)=Bio_bakBI(i,ibioBI)
            DBioBI(i,ibioBI)=0.0_r8
          END DO
        END DO
#endif

#ifdef BIOFLUX
        IF (i.eq.3.and.j.eq.3) THEN
          DO itrc=1,UBst
          ibio=idbio3(itrc)
            DO itrc2=1,UBst
              ibio2=idbio3(itrc)
              BioFlx(ibio,ibio)=bflx(ibio,ibio2)
            END DO
          END DO
        ENDIF
#endif
        DO k=1,N(ng)
          DO i=Istr,Iend
            Bio(i,k,itemp)=t(i,j,k,nstp,itemp)
            Bio(i,k,isalt)=t(i,j,k,nstp,isalt)
            IF (Bio(i,k,itemp) .gt. 35._r8) THEN
              print *, 'Temperature: ',                                 &
     &             Bio(i,k,itemp),i, j, k,ng, yday
              print *,'Tracer: ',t(i,j,k,1,itemp),t(i,j,k,2,itemp),     &
     &              t(i,j,k,3,itemp),Hz(i,j,k),nnew
              print *,'Others: ', z_w(i,j,N(ng)),                       &
     &                  GRID(ng) % h(i,j)
            END IF
            IF ((grid(ng) % h(i,j) + z_w(i,j,N(ng))) .lt.0.0_r8) THEN
              print *, 'zeta & h: ',z_w(i,j,N(ng)),'  ',grid(ng) % h(i,j)
            END IF
          END DO
        END DO
!
! ----------------------------------------------------------------------
!  Calculate Day Length and Surface PAR
! ----------------------------------------------------------------------
!
!  Calculate Day Length
        DO i=Istr,Iend
#if defined DIURNAL_SRFLUX
!  Day Length is already accounted for in ANA_SWRAD so disable correction
          Dl = 24.0_r8
#else
!  Day Length calculation (orig from R. Davis) from latitude and declination.
!  cff2 is Solar declination from Oberhuber (1988) (COADS documentation)
!            lat = 58.00  orig from C code
!            lat = 45.00_r8  test for EPOC
          lat = GRID(ng) % latr(i,j)
          cff1 = 2.0_r8 * pi * ( yday-1.0_r8 ) / 365.0_r8
          cff2 = 0.006918_r8 - 0.399912_r8*cos(cff1)                    &
     &           + 0.070257_r8*sin(cff1) - 0.006758_r8*cos(2*cff1)      &
     &           + 0.000907_r8*sin(2*cff1) - 0.002697_r8*cos(3*cff1)    &
     &           + 0.00148_r8*sin(3*cff1)
          cff3 = lat * pi /180.0_r8
          IF ( abs( -tan(cff3)*tan(cff2) ) .le. 1.0_r8 ) THEN
            cff1 = acos( -tan(cff3)*tan(cff2) ) * 180.0_r8 / pi
            Dl = 2.0_r8 / 15.0_r8 * cff1
          ELSE
            IF ( yday.gt.90.0_r8 .and. yday.lt.270.0_r8 ) THEN
              Dl = 24.0_r8
            ELSE
              Dl = 0.0_r8
            END IF
          END IF
#endif
!  Calculate PAR at the surface
#ifdef KODIAK_IRAD
!  For PAR, Eyeball fit of data from Hinckley''s ezeroday.dat (E d-1 m-2)
          cff2 = 41.0_r8 - 35.0_r8                                   &
     &           * COS( ( 12.0_r8 + yday) * 2.0_r8 * pi / 365.0_r8 )
#else
!  For PAR, use Shortwave radiation ( = surface solar irradiance)
!  converted from deg C m/s to E/m2/day
          cff2 = srflx(i,j) * rho0 * Cp * 0.394848_r8
!         IF (i.eq.3)THEN
!           Stat2(i,1)=cff2
!         END IF
#endif
!-----------------------------------------------------------
!  George GIBSON''s version after Morel 1988 (in Loukos 1977)
!-----------------------------------------------------------
          PAR(i,N(ng)) = PARfrac(ng) * cff2                             &
     &      * exp( k_ext + k_chl*(Bio(i,N(ng),iPhS)/ccr +               &
     &                            Bio(i,N(ng),iPhL)/ccrPhL)**0.428      &
     &      * ( z_r(i,j,N(ng)) -  z_w(i,j,N(ng)) ) )
        END DO
!  Calculate light decay in the water column
#ifdef NEWSHADE
        DO k=N(ng)-1,1,-1
          DO i=Istr,Iend
            cff1 = k_ext * ( z_r(i,j,k) - z_r(i,j,k+1) )
            cff2 = (k_chl*(Bio(i,k+1,iPhS)/ccr +                        &
     &                       Bio(i,k+1,iPhL)/ccrPhL)**0.428_r8)         &
     &                 * ( z_w(i,j,k) - z_r(i,j,k+1) )
            cff3 = (k_chl*(Bio(i,k,iPhS)/ccr +                          &
     &                       Bio(i,k,iPhL)/ccrPhL)**0.428_r8)           &
                       * ( z_r(i,j,k) - z_w(i,j,k) )
            PAR(i,k) = PAR(i,k+1) * EXP(cff1+cff2+cff3)

          END DO
        END DO
#else
!  Version from Sarah''s old C code (probably wrong)
        DO k=N(ng),1,-1
          DO i=Istr,Iend
            cff3 = z_r(i,j,k)+2.5_r8
            IF ( cff3 .gt. -71.0_r8 ) THEN
              cff1 = k_ext + k_chl *                                  &
     &                  ( Bio(i,k,iPhS) + Bio(i,k,iPhL) ) / ccr
            ELSE
              cff1 = 0.077_r8
            END IF
            PAR(i,k) = PARfrac(ng) * cff2 * exp( cff1 * cff3 )
          END DO
        END DO
#endif

! These are read in now
!          TmaxPhS = 20_r8
!          KtBm_PhS = 0.069_r8
!          TmaxPhL = 20
!          KtBm_PhL = 0.069
!          TmaxMZS = 20
!          KtBm_MZS = 0.069
!          TmaxMZL = 20
!          KtBm_MZL = 0.069
        DO k=1,N(ng)
          DO i=Istr,Iend
            HzL(i,k) = Hz(i,j,k)
            Sal1 = Bio(i,k,isalt)
            Temp1 = Bio(i,k,itemp)

!------------------------------
!Compute sigma-t for each depth
!------------------------------
            Dens(i,k) = ComputeDensity(Temp1,Sal1)
          END DO
        END DO

        DO k=0,N(ng)
          DO i=Istr,Iend
            z_wL(i,k) = z_w(i,j,k)
          END DO
        END DO

!************************************************************************
!************************************************************************
! Begin BIOITER LOOP
!************************************************************************
!************************************************************************

        ITER_LOOP: DO Iter=1,BioIter(ng)
!         IF ( .not. downward) THEN
                  !------------------------------------
                  !Make Neocalanus go down if temp > 12
                  !------------------------------------
!           IF (Bio(i,k,itemp) .gt. 12._r8 .and. NCa(k) .gt. 0.2_r8) THEN
!             downward = .true.
!             goto 111
!           END IF
!         END IF
! 111           Continue
          LightLim = 1.0_r8
          NOLim = 1.0_r8
          NHLim = 1.0_r8
          IronLim = 1.0_r8

!=======================================================================
!  Nutrient uptake by Small Phytoplankton
!=======================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !------------------------
              !Growth rate computations
              !------------------------
              Drate = DiS * 10.0_r8 ** (DpS * Bio(i,k,itemp) )
!                Pmax = (2.0_r8 ** Drate - 1.0_r8 )
!Check with Ken - why remove day length fraction scalar
!
              Pmax = (2.0_r8 ** Drate - 1.0_r8 ) * Dl / 24.0_r8

!                 Pmax = 0.8_r8
!               old LightLim = TANH( alphaPhS * PAR(k) / Pmax / ccr )

              !------------------
              !Nitrate limitation
              !------------------
!g             NOLim = Bio(i,k,iNO3) * EXP( -psiPhS * Bio(i,k,iNH4) )  &
!g     &                    / ( k1PhS + Bio(i,k,iNO3) )

              NOLim = Bio(i,k,iNO3) *                                   &
     &                (1-(Bio(i,k,iNH4)/(k2PhS + Bio(i,k,iNH4))))       &
     &                    / ( k1PhS + Bio(i,k,iNO3) )
#ifdef IRON_LIMIT
              !--------------------------------------------------------
              ! Iron Limitation - disabled at concs of 2 micromol Fe m-3
              !---------------------------------------------------------
              IronLim = eps + Bio(i,k,iFe) / (kfePhS + Bio(i,k,iFe))    &
     &                            * (kfePhS + 2._r8) / 2._r8
#endif
              !--------------------------
              !Light limitation function
              !-------------------------
              Par1 = PAR(i,k)
!             IF (i.eq.3) THEN
!               Stat2(i,1)=Par1
!             END IF

              LightLim = GetLightLimIronSml(alphaPhS, Par1,             &
     &                  Pmax,ccr, IronLim)
!             IF (i.eq.3) THEN
!               Stat2(i,2)=LightLim
!             END IF
              !--------------
              !Nitrate uptake
              !--------------
              NOup = Bio(i,k,iPhS) * Pmax * LightLim * NOLim * IronLim
#ifdef STATIONARY

              Stat3(i,k,7)=alphaPhS
              Stat3(i,k,8)=ccr
              Stat3(i,k,9)=IronLim
              Stat3(i,k,10)=Par1


#endif
              !-------------------
              !Ammonium limitation
              !-------------------
              NHLim = Bio(i,k,iNH4) / ( k2PhS + Bio(i,k,iNH4) )

              !-----------------------------
              !Light limitation for ammonium
              !-----------------------------
              LightLim = GetLightLimSml(alphaPhS, Par1, Pmax, ccr)

              !---------------
              !Ammonium uptake
              !---------------
              NHup = Bio(i,k,iPhS) * Pmax * LightLim * NHLim
#ifdef STATIONARY
!                Stat3(i,k,13)=  Pmax * LightLim * NHLim
!                Stat3(i,k,14)= NHLim
#endif
              !-------------------------------
              !Change in nitrate concentration
              !-------------------------------
              DBio(i,k,iNO3) = DBio(i,k,iNO3) - xi * NOup * dtdays

              !--------------------------------
              !Change in ammonium concentration
              !--------------------------------
              DBio(i,k,iNH4) = DBio(i,k,iNH4) - xi * NHup * dtdays

              !----------------------------------------------
              !Change in concentration of small phytoplankton
              !----------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) + ( NOup + NHup ) * dtdays

#ifdef STATIONARY
              Stat3(i,k,11)=  LightLim
              Stat3(i,k,12)= NHLim
#endif
              !-----------------------------------------
              !Primary production of small phytoplankton
              !-----------------------------------------
#ifdef PROD3
              Prod(i,k,iPhS) = Prod(i,k,iPhS) + DBio(i,k,iPhS)
#endif
#ifdef IRON_LIMIT
              !----------------------------
              !Change in iron concentration
              !----------------------------
              DBio(i,k,iFe) = DBio(i,k,iFe) - FeC * NOup * dtdays
#endif
#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN

                bflx(iNO3,iPhS) = bflx(iNO3,iPhS) + NOup*xi
                bflx(iNH4,iPhS) = bflx(iNH4,iPhS) + NHup*xi
              ENDIF
#endif
            END DO
          END DO
!=========================================================================
!  Nutrient uptake by Large Phytoplankton
!=========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !------------------------
              !Growth rate computations
              !------------------------

              Drate = DiL * 10.0_r8 ** (DpL * Bio(i,k,itemp) )
!             Pmax = (2.0_r8 ** Drate - 1.0_r8 )


             !Ken - no day length fraction scalar

!             Pmax = (2.0_r8 ** Drate - 1.0_r8 )
              Pmax = (2.0_r8 ** Drate - 1.0_r8 ) * Dl / 24.0_r8

!             Pmax = 1.5_r8
!             Pmax = 2.2_r8
                !------------------
              !Nitrate limitation
              !------------------
!g            NOLim = Bio(i,k,iNO3) * EXP( -psiPhL * Bio(i,k,iNH4) )  &
!g     &                    / ( k1PhL + Bio(i,k,iNO3) )

              NOLim = Bio(i,k,iNO3) *                                   &
     &                (1-(Bio(i,k,iNH4)/(k2PhL + Bio(i,k,iNH4))))       &
     &                    / ( k1PhL + Bio(i,k,iNO3) )

#ifdef IRON_LIMIT
              !--------------------------------------------------------
              ! Iron Limitation - disabled at cons of 2 micromol Fe m-3
              !--------------------------------------------------------

              IronLim = eps + Bio(i,k,iFe) / (kfePhL + Bio(i,k,iFe))    &
     &                           * (kfePhL + 2._r8) / 2._r8
#endif
              !-------------------------
              !Light limitation function
              !-------------------------
              Par1 = Par(i,k)
              ParMax = Par(i,N(ng))

              !---------------------------------------------------
              !Ken uses composite light curve with iron limitation
              !---------------------------------------------------
!             LightLim = GetLightLimIron(alphaPhL,PAR1,Pmax,            &
!     &                  ccrPhL,IronLim,ParMax)
              LightLim = GetLightLimIron2(alphaPhL, PAR1, Pmax,         &
     &                  ccrPhL, IronLim)
!             LightLim=1.0_r8
!             IronLim=1.0_r8
!             NOLim=1.0_r8
              !--------------
              !Nitrate uptake
              !--------------
              NOup = Bio(i,k,iPhL) * Pmax * LightLim * NOLim * IronLim
#ifdef STATIONARY
              Stat3(i,k,1)=alphaPhL
              Stat3(i,k,2)=ccrPhL
              Stat3(i,k,3)=IronLim
              Stat3(i,k,4)=PAR1
#endif
              !-------------------
              !Ammonium limitation
              !-------------------
              NHLim = Bio(i,k,iNH4) / ( k2PhL + Bio(i,k,iNH4) )

!             NHLim=1.0_r8

              !-------------------------------------------------
              !Ken uses composite light curve without iron limitation
              !-------------------------------------------------
!             LightLim = GetLightLim(alphaPhL,PAR1,Pmax,                &
!     &                  ccrPhL,ParMax)

              !----------------------------
              !Use hyperbolic tangent curve
              !----------------------------
              LightLim = GetLightLim2(alphaPhL, PAR1, Pmax, ccrPhL)
!             LightLim=1.0_r8
!             NHLim=1.0_r8
              !---------------
              !Ammonium uptake
              !---------------
              NHup = Bio(i,k,iPhL) * Pmax * LightLim * NHLim
#ifdef STATIONARY
              Stat3(i,k,5)= LightLim
              Stat3(i,k,6)= NHLim
#endif
              !-------------------------------
              !Change in nitrate concentration
              !-------------------------------
              DBio(i,k,iNO3) = DBio(i,k,iNO3) - xi * NOup * dtdays

              !--------------------------------
              !Change in ammonium concentration
              !--------------------------------
              DBio(i,k,iNH4) = DBio(i,k,iNH4) - xi * NHup * dtdays

              !----------------------------------------------
              !Change in concentration of large phytoplankton
              !----------------------------------------------
              DBio(i,k,iPhL) = DBio(i,k,iPhL) + ( NOup + NHup ) * dtdays
!             DBio(i,k,iPhL) =dtdays
#ifdef STATIONARY
!             Stat3(i,k,7)= NOup + NHup
#endif
              !-----------------------------------------
              !Primary production of large phytoplankton
              !-----------------------------------------
#ifdef PROD3
              Prod(i,k,iPhL) = Prod(i,k,iPhL) + DBio(i,k,iPhL)
#endif
#ifdef IRON_LIMIT
              !----------------------------
              !Change in iron concentration
              !----------------------------
              DBio(i,k,iFe) = DBio(i,k,iFe) - FeC * NOup * dtdays
#endif
#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iNO3,iPhL) = bflx(iNO3,iPhL) + NOup*xi
                bflx(iNH4,iPhL) = bflx(iNH4,iPhL) + NHup*xi
              END IF
#endif
            END DO
          END DO

!=======================================================================
! Grazing by MZS
!=======================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpPhSMZS * Bio(i,k,iPhS)**2                       &
     &                 + fpPhLMZS * Bio(i,k,iPhL)**2

              !------------------
              !Food consumption
              !------------------
!             cff2 = eMZS * Bio(i,k,iMZS) / (fMZS**2 + cff1)
              cff2 = eMZS * Bio(i,k,iMZS) / (fMZS + cff1)
              !------------------------
              !Temperature correction
              !------------------------
!             cff3 = Q10MZS ** ( (Bio(i,k,itemp)-Q10MZST)/ 10.0_r8 )

              cff3 =1.0_r8

              !--------------------------------------------------------
              !Change in small and large phytoplankton due to predation
              !--------------------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) - fpPhSMZS *              &
     &                (Bio(i,k,iPhS)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) - fpPhLMZS *              &
     &                (Bio(i,k,iPhL)**2) * cff2 * cff3 * dtdays

              !----------------------------------------------------
              !Growth of small microzooplankton due to consumption
              !----------------------------------------------------
              DBio(i,k,iMZS) = DBio(i,k,iMZS) +                         &
     &                         gammaMZS * cff1 * cff2 * cff3 * dtdays

              !-------------------------------------
              !Production for small microzooplankton
              !-------------------------------------
#ifdef PROD3
              Prod(i,k,iMZS) = Prod(i,k,iMZS) + DBio(i,k,iMZS)
#endif
              !-------------------------------------------------
              ! Additions to detritus pool - unassimilated food
              !-------------------------------------------------
              DBio(i,k,iDet) = DBio(i,k,iDet) +                         &
     &                 (1.0_r8 - gammaMZS) * cff1 * cff2 * cff3 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iMZS) = bflx(iPhS,iMZS) +                     &
     &                 fpPhSMZS * (Bio(i,k,iPhS)**2) *                  &
     &                 cff2 * cff3 * dtdays*xi
                bflx(iPhL,iMZS) = bflx(iPhL,iMZS) +                     &
     &                 fpPhLMZS * (Bio(i,k,iPhL)**2) *                  &
     &                 cff2 * cff3 * dtdays*xi
                bflx(iMZS,iDet) = bflx(iMZS,iDet)  +                    &
     &                ( 1.0_r8-gammaMZS )*cff1*cff2* cff3 * dtdays*xi
              END IF
#endif
            END DO
          END DO

!========================================================================
! Grazing by MZL
!========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpPhSMZL * Bio(i,k,iPhS)**2                        &
     &               + fpPhLMZL * Bio(i,k,iPhL)**2                      &
     &               + fpMZSMZL * Bio(i,k,iMZS)**2

              !--------------------------------------
              !Food consumption
              !--------------------------------------
!             cff2 = eMZL * Bio(i,k,iMZL) / (fMZL**2 + cff1)
              cff2 = eMZL * Bio(i,k,iMZL) / (fMZL + cff1)
              !--------------------------------------
              !Temperature correction
              !--------------------------------------
!             cff3= Q10MZL ** ( (Bio(i,k,itemp)-Q10MZLT) / 10.0_r8 )

              cff3= 1.0_r8


              !--------------------------------------------------------
              !Change in small and large phytoplankton due to predation
              !--------------------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) - fpPhSMZL *              &
     &                 (Bio(i,k,iPhS)**2) * cff2 * cff3* dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) - fpPhLMZL *              &
     &                 (Bio(i,k,iPhL)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iMZS) = DBio(i,k,iMZS) - fpMZSMZL *              &
     &                 (Bio(i,k,iMZS)**2) * cff2 * cff3 * dtdays

              !--------------------------------
              !Growth of large microzooplankton
              !--------------------------------
              DBio(i,k,iMZL) = DBio(i,k,iMZL) +                         &
     &                        gammaMZL * cff1 * cff2 * cff3 * dtdays

              !------------------------------------
              !Production of large microzooplankton
              !------------------------------------
#ifdef PROD3
              Prod(i,k,iMZL) = Prod(i,k,iMZL) + DBio(i,k,iMZL)
#endif
              !------------------------------------------------
              ! Additions to detritus pool - unassimilated food
              !------------------------------------------------
              DBio(i,k,iDet) = DBio(i,k,iDet) +                         &
     &              (1.0_r8 - gammaMZL) * cff1 * cff2 * cff3 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iMZL) = bflx(iPhS,iMZL) +                     &
     &                      fpPhSMZL * (Bio(i,k,iPhS)**2) *             &
     &                      cff2 * cff3* dtdays*xi
                bflx(iPhL,iMZL) = bflx(iPhL,iMZL) +                     &
     &                     fpPhLMZL *  (Bio(i,k,iPhL)**2) *             &
     &                      cff2 * cff3 * dtdays*xi
                bflx(iMZS,iMZL) = bflx(iPhL,iMZL) +                     &
     &                     fpMZSMZL * (Bio(i,k,iMZS)**2) * cff2         &
     &                     * cff3 * dtdays*xi
                bflx(iMZL,iDet) = bflx(iMZL,iDet)  +                    &
     &                ( 1.0_r8-gammaMZL )*cff1*cff2* cff3 * dtdays*xi
              END IF
#endif
            END DO
          END DO

!==========================================================================
! Grazing and Predation by Copepods
!==========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpPhSCop * Bio(i,k,iPhS)**2                        &
     &               + fpPhLCop * Bio(i,k,iPhL)**2                      &
     &               + fpMZSCop * Bio(i,k,iMZS)**2                      &
     &               + fpMZLCop * Bio(i,k,iMZL)**2

              !--------------------------------------
              !Food consumption
              !--------------------------------------
!             cff2 = eCop * Bio(i,k,iCop) / (fCop**2 + cff1)
              cff2 = eCop * Bio(i,k,iCop) / (fCop + cff1)
              !--------------------------------------
              !Temperature correction
              !--------------------------------------
              cff3 = Q10Cop ** ( (Bio(i,k,itemp)-Q10CopT) / 10.0_r8 )

              !------------------------
              !Growth of small copepods
              !------------------------
              DBio(i,k,iCop) = DBio(i,k,iCop) +                         &
     &                        gammaCop * cff1 * cff2 * cff3 * dtdays
              !------------------
              !Copepod production
              !------------------
#ifdef PROD3
              Prod(i,k,iCop) = Prod(i,k,iCop) + DBio(i,k,iCop)
#endif
              !----------------------------------------------
              !Changes in prey concentration due to predation
              !----------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) -  fpPhSCop               &
     &                * (Bio(i,k,iPhS)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) -  fpPhLCop               &
     &                * (Bio(i,k,iPhL)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iMZS) = DBio(i,k,iMZS) -  fpMZSCop               &
     &                * (Bio(i,k,iMZS)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iMZL) = DBio(i,k,iMZL) -  fpMZLCop               &
     &                * (Bio(i,k,iMZL)**2) * cff2 * cff3 * dtdays

              !------------------------------------------------
              ! Additions to detritus pool - unassimilated food
              !---------------------------
              DBio(i,k,iDetF) = DBio(i,k,iDetF) +                       &
     &            (1.0_r8 - gammaCop) * cff1 * cff2 * cff3 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iCop)=bflx(iPhS,iCop)+                        &
     &                     fpPhSCop * (Bio(i,k,iPhS)**2) *              &
     &                     cff2 * cff3 * dtdays*xi
                bflx(iPhL,iCop)=bflx(iPhL,iCop)+                        &
     &                     fpPhLCop  * (Bio(i,k,iPhL)**2) * cff2 *      &
     &                     cff3 * dtdays*xi
                bflx(iMZS,iCop)=bflx(iMZS,iCop)+                        &
     &                     fpMZSCop  * (Bio(i,k,iMZS)**2) * cff2 *      &
     &                     cff3 * dtdays*xi
                bflx(iMZL,iCop)=bflx(iMZL,iCop)+                        &
     &                      fpMZLCop  * (Bio(i,k,iMZL)**2) * cff2 *     &
     &                     cff3 * dtdays*xi
                bflx(iCop,iDetF) = bflx(iCop,iDetF) +                   &
     &                 ( 1.0_r8-gammaCop )*cff1*cff2*cff3 * dtdays*xi
              END IF
#endif
            END DO
          END DO

!========================================================================
! Grazing and Predation by NCa initiated ON the shelf
!========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpPhSNCa * Bio(i,k,iPhS)**2                        &
     &               + fpPhLNCa * Bio(i,k,iPhL)**2                      &
     &               + fpMZSNCa * Bio(i,k,iMZS)**2                      &
     &               + fpMZLNCa * Bio(i,k,iMZL)**2

              !--------------------------------------
              !Food consumption
              !--------------------------------------

!             cff2 = eNCa * Bio(i,k,iNCaS) / (fNCa**2 + cff1)
              cff2 = eNCa * Bio(i,k,iNCaS) / (fNCa + cff1)
              !--------------------------------------
              !Temperature correction
              !--------------------------------------
              cff3 = Q10NCa ** ( (Bio(i,k,itemp)-Q10NCaT) / 10.0_r8 )

              !--------------------
              !Growth of Neocalanus
              !--------------------
              DBio(i,k,iNCaS) = DBio(i,k,iNCaS) +                       &
     &                        gammaNCa * cff1 * cff2 * cff3 * dtdays

              !---------------------
              !Neocalanus production
              !---------------------
#ifdef PROD3
              Prod(i,k,iNCaS) = Prod(i,k,iNCaS) + DBio(i,k,iNCaS)
#endif
              !----------------------------------------------
              !Changes in prey concentration due to predation
              !----------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) -                         &
     &             fpPhSNCa * (Bio(i,k,iPhS)**2) * cff2 * cff3 *dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) -                         &
     &             fpPhLNCa * (Bio(i,k,iPhL)**2) * cff2 * cff3 *dtdays
              DBio(i,k,iMZS) = DBio(i,k,iMZS) -                         &
     &             fpMZSNCa * (Bio(i,k,iMZS)**2) * cff2 * cff3 *dtdays
              DBio(i,k,iMZL) = DBio(i,k,iMZL) -                         &
     &             fpMZLNCa * (Bio(i,k,iMZL)**2) * cff2 * cff3 *dtdays


          !-------------------------------------------------
          ! Additions to Fast Sinking detritus pool - unassimilated food
          !-------------------------------------------------
              DBio(i,k,iDetF) = DBio(i,k,iDetF) +                       &
     &               (1.0_r8 - gammaNCa) * cff1 * cff2 * cff3 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iNCaS)=bflx(iPhS,iNCaS)+                      &
     &                fpPhSNCa * (Bio(i,k,iPhS)**2) * cff2 *            &
     &                cff3 *dtdays*xi
                bflx(iPhL,iNCaS)=bflx(iPhL,iNCaS)+                      &
     &                fpPhLNCa * (Bio(i,k,iPhL)**2) * cff2 *            &
     &                cff3 *dtdays*xi
                bflx(iMZS,iNCaS)=bflx(iMZS,iNCaS)+                      &
     &                fpMZSNCa * (Bio(i,k,iMZS)**2) * cff2 *            &
     &                cff3 *dtdays*xi
                bflx(iMZL,iNCaS)=bflx(iMZL,iNCaS)+                      &
     &                fpMZLNCa * (Bio(i,k,iMZL)**2) * cff2 *            &
     &                cff3 *dtdays*xi
                bflx(iNCaS,iDetF) = bflx(iNCaS,iDetF) +                 &
     &            ( 1.0_r8-gammaNCa )*cff1*cff2* cff3 * dtdays*xi
              END IF
#endif
            END DO
          END DO

!=========================================================================
! Grazing and Predation by Euphuasiids initiated ON the shelf
!=========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpPhSEup * Bio(i,k,iPhS)**2                        &
     &             + fpPhLEup * Bio(i,k,iPhL)**2                        &
     &             + fpMZSEup * Bio(i,k,iMZS)**2                        &
     &             + fpMZLEup * Bio(i,k,iMZL)**2                        &
     &             + fpCopEup * Bio(i,k,iCop)**2

              !--------------------------------------
              !Food consumption
              !--------------------------------------
!             cff2 = eEup * Bio(i,k,iEupS) / (fEup**2 + cff1)
              cff2 = eEup * Bio(i,k,iEupS) / (fEup + cff1)

              !--------------------------------------
              !Temperature correction
              !--------------------------------------
              cff3 = Q10Eup ** ( (Bio(i,k,itemp)-Q10EupT) / 10.0_r8 )

              !---------------------
              !Growth of Euphausiids
              !---------------------

              DBio(i,k,iEupS) = DBio(i,k,iEupS) +                     &
     &                       gammaEup * cff1 * cff2 * cff3 * dtdays

              !---------------------
              !Euphausiid production
              !---------------------
#ifdef PROD3
              Prod(i,k,iEupS) = Prod(i,k,iEupS) + DBio(i,k,iEupS)
#endif
              !----------------------------------------------
              !Changes in prey concentration due to predation
              !----------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) -                         &
     &             fpPhSEup * (Bio(i,k,iPhS)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) -                         &
     &             fpPhLEup * (Bio(i,k,iPhL)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iMZS) = DBio(i,k,iMZS) -                         &
     &             fpMZSEup * (Bio(i,k,iMZS)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iMZL) = DBio(i,k,iMZL) -                         &
     &             fpMZLEup * (Bio(i,k,iMZL)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iCop) = DBio(i,k,iCop) -                         &
     &             fpCopEup * (Bio(i,k,iCop)**2) * cff2 * cff3 * dtdays

              !-------------------------------------------------
              ! Additions to Fast Sinking detritus pool- unassimilated food
              !-------------------------------------------------
              DBio(i,k,iDetF) = DBio(i,k,iDetF) +                       &
     &            (1.0_r8 - gammaEup) * cff1 * cff2 * cff3 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iEupS)=bflx(iPhS,iEupS)+                      &
     &                   fpPhSEup * (Bio(i,k,iPhS)**2) * cff2 *         &
     &                   cff3 * dtdays*xi
                bflx(iPhL,iEupS)=bflx(iPhL,iEupS)+                      &
     &                    fpPhLEup * (Bio(i,k,iPhL)**2) * cff2 *        &
     &                   cff3 * dtdays*xi
                bflx(iMZS,iEupS)=bflx(iMZS,iEupS)+                      &
     &                    fpMZSEup * (Bio(i,k,iMZS)**2) * cff2 *        &
     &                   cff3 * dtdays*xi
                bflx(iMZL,iEupS)=bflx(iMZL,iEupS)+                      &
     &                    fpMZLEup * (Bio(i,k,iMZL)**2) * cff2 *        &
     &                   cff3 * dtdays*xi
                bflx(iCop,iEupS)=bflx(iCop,iEupS)+                      &
     &                    fpCopEup * (Bio(i,k,iCop)**2) * cff2 *        &
     &                   cff3 * dtdays*xi
                bflx(iEupS,iDetF) = bflx(iEupS,iDetF) +                 &
     &                ( 1.0_r8-gammaEup )*cff1*cff2* cff3 * dtdays*xi
              END IF
#endif
            END DO
         END DO

!========================================================================
! Grazing and Predation by NCa initiated OFF the shelf
!========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpPhSNCa * Bio(i,k,iPhS)**2                        &
     &             + fpPhLNCa * Bio(i,k,iPhL)**2                        &
     &             + fpMZSNCa * Bio(i,k,iMZS)**2                        &
     &             + fpMZLNCa * Bio(i,k,iMZL)**2

              !--------------------------------------
              !Food consumption
              !--------------------------------------

!             cff2 = eNCa * Bio(i,k,iNCaO) / (fNCa**2 + cff1)
              cff2 = eNCa * Bio(i,k,iNCaO) / (fNCa + cff1)
              !--------------------------------------
              !Temperature correction
              !--------------------------------------
              cff3 = Q10NCa ** ( (Bio(i,k,itemp)-Q10NCaT) / 10.0_r8 )

              !--------------------
              !Growth of Neocalanus
              !--------------------
              DBio(i,k,iNCaO) = DBio(i,k,iNCaO) +                       &
     &                        gammaNCa * cff1 * cff2 * cff3 * dtdays

              !---------------------
              !Neocalanus production
              !---------------------
#ifdef PROD3

              Prod(i,k,iNCaO) = Prod(i,k,iNCaO) + DBio(i,k,iNCaO)
#endif
              !----------------------------------------------
              !Changes in prey concentration due to predation
              !----------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) -                         &
     &             fpPhSNCa * (Bio(i,k,iPhS)**2) * cff2 * cff3 *dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) -                         &
     &             fpPhLNCa * (Bio(i,k,iPhL)**2) * cff2 * cff3 *dtdays
              DBio(i,k,iMZS) = DBio(i,k,iMZS) -                         &
     &             fpMZSNCa * (Bio(i,k,iMZS)**2) * cff2 * cff3 *dtdays
              DBio(i,k,iMZL) = DBio(i,k,iMZL) -                         &
     &             fpMZLNCa * (Bio(i,k,iMZL)**2) * cff2 * cff3 *dtdays

              !-------------------------------------------------
              ! Additions to detritus pool - unassimilated food
              !-------------------------------------------------
              DBio(i,k,iDetF) = DBio(i,k,iDetF) +                       &
     &                 (1.0_r8 - gammaNCa) * cff1 * cff2 * cff3 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iNCaO)=bflx(iPhS,iNCaO)+                      &
     &                   fpPhSNCa * (Bio(i,k,iPhS)**2) * cff2 *         &
     &                   cff3 *dtdays*xi
                bflx(iPhL,iNCaO)=bflx(iPhL,iNCaO)+                      &
     &                   fpPhLNCa * (Bio(i,k,iPhL)**2) * cff2 *         &
     &                   cff3 *dtdays *xi
                bflx(iMZS,iNCaO)=bflx(iMZS,iNCaO)+                      &
     &                   fpMZSNCa * (Bio(i,k,iMZS)**2) * cff2 *         &
     &                   cff3 *dtdays *xi
                bflx(iMZL,iNCaO)=bflx(iMZL,iNCaO)+                      &
     &                   fpMZLNCa * (Bio(i,k,iMZL)**2) * cff2 *         &
     &                   cff3 *dtdays*xi
                bflx(iNCaO,iDetF) = bflx(iNCaO,iDetF) +                 &
     &               ( 1.0_r8-gammaNCa )*cff1*cff2* cff3 * dtdays*xi
              END IF
#endif
            END DO
          END DO

!=========================================================================
! Grazing and Predation by Euphuasiids initiated OFF the shelf
!=========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpPhSEup * Bio(i,k,iPhS)**2                        &
     &             + fpPhLEup * Bio(i,k,iPhL)**2                        &
     &             + fpMZSEup * Bio(i,k,iMZS)**2                        &
     &             + fpMZLEup * Bio(i,k,iMZL)**2                        &
     &             + fpCopEup * Bio(i,k,iCop)**2

              !--------------------------------------
              !Food consumption
              !--------------------------------------
!             cff2 = eEup * Bio(i,k,iEupO) / (fEup**2 + cff1)
              cff2 = eEup * Bio(i,k,iEupO) / (fEup + cff1)

              !--------------------------------------
              !Temperature correction
              !--------------------------------------
              cff3 = Q10Eup ** ( (Bio(i,k,itemp)-Q10EupT) / 10.0_r8 )

              !---------------------
              !Growth of Euphausiids
              !---------------------

              DBio(i,k,iEupO) = DBio(i,k,iEupO) +                       &
     &                       gammaEup * cff1 * cff2 * cff3 * dtdays

              !---------------------
              !Euphausiid production
              !---------------------
#ifdef PROD3
              Prod(i,k,iEupO) = Prod(i,k,iEupO) + DBio(i,k,iEupO)
#endif
              !----------------------------------------------
              !Changes in prey concentration due to predation
              !----------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) -                         &
     &            fpPhSEup * (Bio(i,k,iPhS)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) -                         &
     &            fpPhLEup * (Bio(i,k,iPhL)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iMZS) = DBio(i,k,iMZS) -                         &
     &            fpMZSEup * (Bio(i,k,iMZS)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iMZL) = DBio(i,k,iMZL) -                         &
     &            fpMZLEup * (Bio(i,k,iMZL)**2) * cff2 * cff3 * dtdays
              DBio(i,k,iCop) = DBio(i,k,iCop) -                         &
     &            fpCopEup * (Bio(i,k,iCop)**2) * cff2 * cff3 * dtdays


              !-------------------------------------------------
              ! Additions to detritus pool- unassimilated food
              !-------------------------------------------------
              DBio(i,k,iDetF) = DBio(i,k,iDetF) +                       &
     &            (1.0_r8 - gammaEup) * cff1 * cff2 * cff3 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iEupO)=bflx(iPhS,iEupO)+                      &
     &                    fpPhSEup * (Bio(i,k,iPhS)**2) * cff2 *        &
     &                    cff3 * dtdays*xi
                bflx(iPhL,iEupO)=bflx(iPhL,iEupO)+                      &
     &                    fpPhLEup * (Bio(i,k,iPhL)**2) * cff2 *        &
     &                    cff3 * dtdays  *xi
                bflx(iMZS,iEupO)=bflx(iMZS,iEupO)+                      &
     &                    fpMZSEup * (Bio(i,k,iMZS)**2) * cff2 *        &
     &                    cff3 * dtdays  *xi
                bflx(iMZL,iEupO)=bflx(iMZL,iEupO)+                      &
     &                     fpMZLEup * (Bio(i,k,iMZL)**2) * cff2 *       &
     &                    cff3 * dtdays *xi
                bflx(iCop,iEupO)=bflx(iCop,iEupO)+                      &
     &                      fpCopEup * (Bio(i,k,iCop)**2) * cff2 *      &
     &                    cff3 * dtdays*xi
                bflx(iEupO,iDetF) = bflx(iEupO,iDetF) +                 &
     &                 ( 1.0_r8-gammaEup )*cff1*cff2* cff3 * dtdays*xi
              END IF
#endif
            END DO
          END DO
!=========================================================================
! Grazing and Predation by Jellyfish
!=========================================================================
#ifdef JELLY
          DO k=1,N(ng)
            DO i=Istr,Iend

              !----------------
              !Food preferences
              !----------------
              cff1 = fpCopJel * Bio(i,k,iCop) +                         &
     &               fpNCaJel * Bio(i,k,iNCaS) +                        &
     &               fpNCaJel * Bio(i,k,iNCaO) +                        &
     &               fpEupJel * Bio(i,k,iEupS) +                        &
     &               fpEupJel * Bio(i,k,iEupO)

              !--------------------------------------
              !Food consumption (Linear)
              !--------------------------------------
              cff2 = eJel

              !--------------------------------------
              !Temperature correction
              !--------------------------------------
              cff3= Q10Jele ** ((Bio(i,k,itemp)-Q10JelTe) / 10.0_r8)

              !---------------------
              !Growth of Jellies
              !---------------------
              DBio(i,k,iJel) =  DBio(i,k,iJel) +                        &
     &                 gammaJel * cff1 * cff2 * cff3 *                  &
     &                 Bio(i,k,iJel)*dtdays

              !---------------------
              !Jellyfish production
              !---------------------
#ifdef PROD3
              Prod(i,k,iJel) = Prod(i,k,iJel) + DBio(i,k,iJel)
#endif
#ifdef STATIONARY2
!             Stat2(i,1)=cff1 * cff2 * cff3 *  Bio(i,k,iJel)*dtdays
#endif
              !----------------------------------------------
              !Changes in prey concentration due to predation
              !----------------------------------------------

              DBio(i,k,iCop) = DBio(i,k,iCop) - fpCopJel *              &
     &                     Bio(i,k,iCop)*Bio(i,k,iJel)* cff2 *          &
     &                     cff3 * dtdays
              DBio(i,k,iEupS) = DBio(i,k,iEupS) - fpEupJel *            &
     &                     Bio(i,k,iEupS)*Bio(i,k,iJel)* cff2 *         &
     &                     cff3 * dtdays
              DBio(i,k,iNCaS) = DBio(i,k,iNCaS) - fpNCaJel *            &
     &                     Bio(i,k,iNCaS)*Bio(i,k,iJel)* cff2 *         &
     &                     cff3 * dtdays
              DBio(i,k,iEupO) = DBio(i,k,iEupO) - fpEupJel *            &
     &                     Bio(i,k,iEupO)*Bio(i,k,iJel)* cff2 *         &
     &                     cff3 * dtdays
              DBio(i,k,iNCaO) = DBio(i,k,iNCaO) - fpNCaJel *            &
     &                     Bio(i,k,iNCaO)*Bio(i,k,iJel)* cff2 *         &
     &                     cff3 * dtdays
              DBio(i,k,iDetF)=DBio(i,k,iDetF) +(1-gammaJel)             &
     &                     * cff1 * cff2 * cff3 * Bio(i,k,iJel)*dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iCop,iJel)=bflx(iCop,iJel)+                        &
     &                   fpCopJel * Bio(i,k,iCop)*Bio(i,k,iJel)* cff2 * &
     &                   cff3 * dtdays*xi
                bflx(iNCaS,iJel)=bflx(iNCaS,iJel)+                      &
     &                   fpNCaJel * Bio(i,k,iNCaS)*Bio(i,k,iJel)* cff2 *&
     &                   cff3 * dtdays*xi
                bflx(iEupS,iJel)=bflx(iEupS,iJel)+                      &
     &                   fpEupJel * Bio(i,k,iEupS)*Bio(i,k,iJel)* cff2 *&
     &                   cff3 * dtdays*xi
                bflx(iNCaO,iJel)=bflx(iNCaO,iJel)+                      &
     &                   fpNCaJel * Bio(i,k,iNCaO)*Bio(i,k,iJel)* cff2 *&
     &                   cff3 * dtdays*xi
                bflx(iEupO,iJel)=bflx(iEupO,iJel)+                      &
     &                   fpEupJel * Bio(i,k,iEupO)*Bio(i,k,iJel)* cff2 *&
     &                   cff3 * dtdays*xi
                bflx(iJel,iDetF) = bflx(iJel,iDetF) +                   &
     &                   ( 1.0_r8-gammaJel)*cff1*cff2* cff3             &
     &                   * Bio(i,k,iJel)*dtdays*xi
              END IF
#endif
            END DO
          END DO
#endif

!=======================================================================
! Phytoplankton Linear Mortality and Senescence Terms
!=======================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend
              cff1 = MAX( minmPhS , maxmPhS -                           &
     &                ( maxmPhS - minmPhS ) * Bio(i,k,iNO3) / NcritPhS )
              cff2 = MAX( minmPhL , maxmPhL -                           &
     &                ( maxmPhL - minmPhL ) * Bio(i,k,iNO3) / NcritPhL )
#ifdef STATIONARY
!             Stat3(i,k,8)=cff2
              Stat3(i,k,16)=cff1
#endif
              DBio(i,k,iPhS) = DBio(i,k,iPhS) -                         &
     &                       cff1 * Bio(i,k,iPhS) * dtdays
              DBio(i,k,iPhL) = DBio(i,k,iPhL) -                         &
     &                       cff2 * Bio(i,k,iPhL) * dtdays

              !-------------------------------------------------
              ! Additions to detritus pool - phytoplankton mort
              !-------------------------------------------------

              DBio(i,k,iDet) = DBio(i,k,iDet) +                         &
     &                     ( cff1 * Bio(i,k,iPhS)) * dtdays
              DBio(i,k,iDetF) = DBio(i,k,iDetF) +                       &
     &                     (cff2 * Bio(i,k,iPhL) ) * dtdays
#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iDet)= bflx(iPhS,iDet)                        &
     &                     + cff1*Bio(i,k,iPhS)* dtdays*xi
                bflx(iPhL,iDetF)= bflx(iPhL,iDetF)                      &
     &                     + cff2*Bio(i,k,iPhL)* dtdays*xi
              END IF
#endif
            END DO
          END DO
!=======================================================================
! Microzooplankton Mortality - use only linear OR QUADRATIC
!=======================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !---------
              ! Linear   (George)
              !---------
              DBio(i,k,iMZS) = DBio(i,k,iMZS) -                         &
     &                         mMZS * Bio(i,k,iMZS) * dtdays
              DBio(i,k,iMZL) = DBio(i,k,iMZL) -                         &
     &                         mMZL * Bio(i,k,iMZL) * dtdays
              !---------
              !Quadratic (Ken)
              !---------

!               DBio(i,k,iMZS) = DBio(i,k,iMZS) -                       &
!     &                         mpredMZS*dtdays*Bio(i,k,iMZS)**2
!               DBio(i,k,iMZL) = DBio(i,k,iMZL) -                       &
!     &                         mpredMZL*dtdays*Bio(i,k,iMZL)**2


              !-------------------------------------------------
              ! Additions to detritus pool - natural microzoo mortality
              !-------------------------------------------------
!if linear     (George)
              DBio(i,k,iDet) = DBio(i,k,iDet) +                         &
     &         (mMZS * Bio(i,k,iMZS) + mMZL * Bio(i,k,iMZL)) * dtdays
!if quadratic (Ken)
!               DBio(i,k,iDet) = DBio(i,k,iDet)                         &
!     &                   + (mpredMZS * Bio(i,k,iMZS)**2                &
!     &                   + mpredMZL * Bio(i,k,iMZL)**2 ) * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iMZS,iDet)= bflx(iMZS,iDet)                        &
     &                 + mMZS * Bio(i,k,iMZS) * dtdays*xi
                bflx(iMZL,iDet)= bflx(iMZL,iDet)                        &
     &                  + mMZL * Bio(i,k,iMZL) * dtdays*xi
              END IF
#endif
            END DO
          END DO

!==================================================================
! Mesozooplankton Mortality (Closure terms)
!==================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              DBio(i,k,iCop) = DBio(i,k,iCop) -                         &
     &                         mpredCop*dtdays*Bio(i,k,iCop)**2
              DBio(i,k,iNCaS) = DBio(i,k,iNCaS) -                       &
     &                         mpredNCa*dtdays*Bio(i,k,iNCaS)**2
              DBio(i,k,iEupS) = DBio(i,k,iEupS) -                       &
     &                         mpredEup*dtdays*Bio(i,k,iEupS)**2
              DBio(i,k,iNCaS) = DBio(i,k,iNCaS) -                       &
     &                         mpredNCa*dtdays*Bio(i,k,iNCaS)**2
              DBio(i,k,iEupS) = DBio(i,k,iEupS) -                       &
     &                         mpredEup*dtdays*Bio(i,k,iEupS)**2

              !---------------------------------
              !Detritus from nonlinear mortality
              !---------------------------------
              DBio(i,k,iDetF) = DBio(i,k,iDetF)                         &
     &                   +(mpredCop * Bio(i,k,iCop)**2                  &
     &                   + mpredNCa * Bio(i,k,iNCaS)**2                 &
     &                   + mpredEup * Bio(i,k,iEupS)**2                 &
     &                   + mpredNCa * Bio(i,k,iNCaO)**2                 &
     &                   + mpredEup * Bio(i,k,iEupO)**2                 &
     &                     ) * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iCop,iDetF)= bflx(iCop,iDetF)+                     &
     &                 mpredCop*dtdays*xi*Bio(i,k,iCop)**2
                bflx(iNcaS,iDetF)= bflx(iNCaS,iDetF)+                   &
     &                 mpredNCa*dtdays*xi*Bio(i,k,iNCaS)**2
                bflx(iEupS,iDetF)= bflx(iEupS,iDetF)+                   &
     &                 mpredEup*dtdays*xi*Bio(i,k,iEupS)**2
                bflx(iNcaO,iDetF)= bflx(iNCaO,iDetF)+                   &
     &                 mpredNCa*dtdays*xi*Bio(i,k,iNCaO)**2
                bflx(iEupO,iDetF)= bflx(iEupO,iDetF)+                   &
     &                 mpredEup*dtdays*xi*Bio(i,k,iEupO)**2
              END IF
#endif

#if defined JELLY
              DBio(i,k,iJel) = DBio(i,k,iJel) - mpredJel *              &
     &                              Bio(i,k,iJel) * dtdays
              DBio(i,k,iDetF) = DBio(i,k,iDetF)                         &
     &                   + mpredJel * Bio(i,k,iJel) * dtdays
#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iJel,iDet)= bflx(iJel,iDet)+                       &
     &                 mpredJel*dtdays*Bio(i,k,iJel)*xi
              END IF
#endif

#endif
            END DO
          END DO

!===============================================
!Phytoplankton respiration losses
!===============================================
          DO k=1,N(ng)
            DO i=Istr,Iend
              BasalMet = respPhS

#ifdef IRON_LIMIT
              !----------------------------------------------------------
              !Correct basal metabolism for iron limitation
              !----------------------------------------------------------

              Iron1 = Bio(i,k,iFe)
              respPh = respPhS
              kfePh = kfePhS
              BasalMet = GetBasalMetabolism(respPh,kfePh,Iron1)

#endif
              !---------------------------------
              !Arhonditsis temperature functions
              !---------------------------------
              TempFuncPhS(i,k) = GetPhytoResp2(Temp1,TmaxPhS,           &
     &              KtBm_PhS)

              !----------------------------------------------
              !Change in concentration of Small Phytoplankton
              !----------------------------------------------
              DBio(i,k,iPhS) = DBio(i,k,iPhS) -                         &
     &                  TempFuncPhS(i,k)*BasalMet*dtdays*Bio(i,k,iPhS)
#ifdef STATIONARY
              Stat3(i,k,16)=Stat3(i,k,16)+TempFuncPhS(i,k)*BasalMet
#endif

              !----------------------------------------------
              !I use for conservation of mass - Ken not using this
              !----------------------------------------------
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &                xi * TempFuncPhS(i,k)*BasalMet*dtdays*            &
     &                Bio(i,k,iPhS)

              !-----------------------------------------
              !Primary production of Small phytoplankton
              !-----------------------------------------
#ifdef PROD3
!g            Prod(i,k,iPHS) = Prod(i,k,iPHS) -                         &
!g     &               TempFuncPhS(i,k)*BasalMet*dtdays*Bio(i,k,iPhS)
#endif
              BasalMet = respPhL
#ifdef IRON_LIMIT

              !----------------------------------------------------------
              !Correct basal metabolism for iron limitation
              !----------------------------------------------------------
              respPh = respPhL
              kfePh = kfePhL
              BasalMet = GetBasalMetabolism(respPh,kfePh,Iron1)
#endif
              !---------------------------------
              !Arhonditsis temperature functions
              !---------------------------------
              TempFuncPhL(i,k) = GetPhytoResp2(Temp1,TmaxPhL,           &
     &              KtBm_PhL)

              !----------------------------------------------
              !Change in concentration of Large Phytoplankton
              !----------------------------------------------

              DBio(i,k,iPhL) = DBio(i,k,iPhL) -                         &
     &               TempFuncPhL(i,k)*BasalMet*dtdays*Bio(i,k,iPhL)
#ifdef STATIONARY
!             Stat3(i,k,8)=Stat3(i,k,8)+TempFuncPhL(i,k)*BasalMet
#endif
              !----------------------------------------------
              !I use for conservation of mass - Ken not using this
              !----------------------------------------------
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi * TempFuncPhL(i,k)*BasalMet*dtdays*Bio(i,k,iPhL)


              !-----------------------------------------
              !Primary production of Large phytoplankton
              !-----------------------------------------
#ifdef PROD3
!g            Prod(i,k,iPHL) = Prod(i,k,iPHL) -                         &
!g     &               TempFuncPhL(i,k)*BasalMet*dtdays*Bio(i,k,iPhL)
#endif


#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iPhS,iNH4)= bflx(iPhS,iNH4)+                       &
     &                TempFuncPhS(i,k)*BasalMet*dtdays*Bio(i,k,iPhS)*xi
                bflx(iPhL,iNH4)= bflx(iPhL,iNH4)+                       &
     &                TempFuncPhL(i,k)*BasalMet*dtdays*Bio(i,k,iPhL)*xi
              END IF
#endif

            END DO
          END DO

!======================================================
!Microzooplankton respiration losses
!======================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

          !----------------------
          !Small Microzooplankton
          !----------------------
              !---------------------------------
              !Arhonditsis temperature functions
              !---------------------------------
              TempFuncMZS(i,k) = GetPhytoResp2(Temp1,TmaxMZS,           &
     &              KtBm_MZS)

              BasalMet = respMZS

              !----------------------------------------------
              !Change in concentration of small microzooplankton
              !---------------------------------------------
              DBio(i,k,iMZS) = DBio(i,k,iMZS) -                         &
     &                TempFuncMZS(i,k)*BasalMet*dtdays*Bio(i,k,iMZS)

              !------------------------------
              !Small Microzooplankton production
              !------------------------------
#ifdef PROD3
              Prod(i,k,iMZS) = Prod(i,k,iMZS) -                         &
     &              TempFuncMZS(i,k)*BasalMet*dtdays*Bio(i,k,iMZS)
#endif
              !----------------------------------------------------------
              !Add ammonium to correct for excretion related to metabolism
              !----------------------------------------------------------
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &              xi*(TempFuncMZS(i,k)*BasalMet*dtdays*Bio(i,k,iMZS))

         !----------------------
         !Large Microzooplankton
         !----------------------
              !---------------------------------
              !Arhonditsis temperature functions
              !---------------------------------
              TempFuncMZL(i,k) = GetPhytoResp2(Temp1,TmaxMZL,           &
     &              KtBm_MZL)

              BasalMet = respMZL

              !----------------------------------------------
              !Change in concentration of large microzooplankton
              !---------------------------------------------
              DBio(i,k,iMZL) = DBio(i,k,iMZL) -                         &
     &               TempFuncMZL(i,k)*BasalMet*dtdays*Bio(i,k,iMZL)


              !----------------------
              !Large Microzooplankton net production
              !----------------------
#ifdef PROD3
              Prod(i,k,iMZL) = Prod(i,k,iMZL) -                         &
     &             TempFuncMZL(i,k)*BasalMet*dtdays*Bio(i,k,iMZL)
#endif
              !----------------------------------------------------------
              !Add ammonium to correct for excretion related to metabolism
              !----------------------------------------------------------
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi*(TempFuncMZL(i,k)*BasalMet*dtdays*Bio(i,k,iMZL))

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iMZS,iNH4)= bflx(iMZS,iNH4)+                       &
     &                TempFuncMZS(i,k)*BasalMet*dtdays*Bio(i,k,iMZS)*xi
                bflx(iMZL,iNH4)= bflx(iMZL,iNH4)+                       &
     &                TempFuncMZL(i,k)*BasalMet*dtdays*Bio(i,k,iMZL)*xi
              END IF
#endif
            END DO
          END DO

!======================================================
!Mesozooplankton respiration losses
!======================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !-------------------------------
              ! Copepod respiration correction
              !-------------------------------
              TempFuncCop(i,k) = GetCopepodResp(Temp1,respCop,          &
     &                             ktbmC,TrefC)

              !----------------------------------
              ! Neocalanus respiration correction
              !----------------------------------
              TempFuncNeo(i,k) = GetCopepodResp(Temp1,respNCa,          &
     &                             ktbmN,TrefN)

              !----------------------------------
              ! Euphausiid respiration correction
              !----------------------------------
              TempFuncEup(i,k) = GetCopepodResp(Temp1,respEup,          &
     &                             ktbmE,TrefE)

              !-----------------------------------------------------
              !Change in concentration from small copepod respiration
              !------------------------------------------------------
              DBio(i,k,iCop) = DBio(i,k,iCop) -                         &
     &                        TempFuncCop(i,k)*Bio(i,k,iCop)*dtdays

!             TestVal(i,k) = DBio(i,k,iNCa)
#ifdef PROD3
              Prod(i,k,iCop) = Prod(i,k,iCop) -                         &
     &            TempFuncCop(i,k)*Bio(i,k,iCop)*dtdays
#endif
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi*(TempFuncCop(i,k)*dtdays*Bio(i,k,iCop))

              !-------------------------------------------------------
              !Change in concentration from large copepods respiration
              !-------------------------------------------------------
              DBio(i,k,iNCaS) = DBio(i,k,iNCaS) -                       &
     &                        TempFuncNeo(i,k)*Bio(i,k,iNCaS)*dtdays
              DBio(i,k,iNCaO) = DBio(i,k,iNCaO) -                       &
     &                        TempFuncNeo(i,k)*Bio(i,k,iNCaO)*dtdays
#ifdef PROD3
              Prod(i,k,iNCaS) = Prod(i,k,iNCaS) -                       &
     &             TempFuncNeo(i,k)*Bio(i,k,iNCaS)*dtdays
              Prod(i,k,iNCaO) = Prod(i,k,iNCaO) -                       &
     &             TempFuncNeo(i,k)*Bio(i,k,iNCaO)*dtdays
#endif
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi*(TempFuncNeo(i,k)*dtdays*Bio(i,k,iNCaS))
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi*(TempFuncNeo(i,k)*dtdays*Bio(i,k,iNCaO))
              !---------------------------------------------------
              !Change in concentration from euphausiid respiration
              !---------------------------------------------------
              DBio(i,k,iEupS) = DBio(i,k,iEupS) -                       &
     &                        TempFuncEup(i,k)*Bio(i,k,iEupS)*dtdays
              DBio(i,k,iEupO) = DBio(i,k,iEupO) -                       &
     &                        TempFuncEup(i,k)*Bio(i,k,iEupO)*dtdays
#ifdef PROD3
              Prod(i,k,iEupS) = Prod(i,k,iEupS) -                       &
     &              TempFuncEup(i,k)*Bio(i,k,iEupS)*dtdays
              Prod(i,k,iEupO) = Prod(i,k,iEupO) -                       &
     &              TempFuncEup(i,k)*Bio(i,k,iEupO)*dtdays
#endif
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi*(TempFuncEup(i,k)*dtdays*Bio(i,k,iEupS))
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi*(TempFuncEup(i,k)*dtdays*Bio(i,k,iEupO))
#if defined JELLY
!             TempFuncJel(i,k) = =bmJ*( Q10Jelr **                      &
!     &                  ( (Temp(k)-Q10JelTr)/10.0_r8))

              TempFuncJel(i,k) = GetJelResp(Temp1,respJel,ktbmJ,TrefJ)
              DBio(i,k,iJel) = DBio(i,k,iJel) - TempFuncJel(i,k) *      &
     &                            Bio(i,k,iJel)*dtdays
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &               xi*(TempFuncJel(i,k)*dtdays*Bio(i,k,iJel))
# ifdef PROD3
              Prod(i,k,iJel) = Prod(i,k,iJel) -                         &
     &              TempFuncJel(i,k)*Bio(i,k,iJel)*dtdays
# endif
#endif

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iCop,iNH4)= bflx(iCop,iNH4)+                       &
     &                  TempFuncCop(i,k)*Bio(i,k,iCop)*dtdays*xi
                bflx(iNCaS,iNH4)= bflx(iNCaS,iNH4)+                     &
     &                  TempFuncNeo(i,k)*Bio(i,k,iNCaS)*dtdays*xi
                bflx(iNCaO,iNH4)= bflx(iNCaO,iNH4)+                     &
     &                  TempFuncNeo(i,k)*Bio(i,k,iNCaO)*dtdays*xi
                bflx(iEupS,iNH4)= bflx(iEupS,iNH4)+                     &
     &                  TempFuncEup(i,k)*Bio(i,k,iEupS)*dtdays  *xi
                bflx(iEupO,iNH4)= bflx(iEupO,iNH4)+                     &
     &                  TempFuncEup(i,k)*Bio(i,k,iEupO)*dtdays*xi
                bflx(iJel,iNH4)= bflx(iJel,iNH4)+                       &
     &                  TempFuncJel(i,k) * Bio(i,k,iJel)*dtdays*xi
              END IF
#endif
            END DO
          END DO
!=========================================================================
! Molting:
!=========================================================================

          !-----------------------------------------------------
          ! NOTE: It is unclear where molting equation came from.
          ! This is present only for euphausiids, not copepods
          !-----------------------------------------------------

          DO k=1,N(ng)
            DO i=Istr,Iend
              cff1 = 0.02_r8 / (10.0_r8 - 0.4_r8 * Bio(i,k,itemp))*     &
     &                       Bio(i,k,iEupS)
              DBio(i,k,iDet) = DBio(i,k,iDet) + cff1 * dtdays
              DBio(i,k,iEupS) = DBio(i,k,iEupS) - cff1 * dtdays

              cff1 = 0.02_r8 / (10.0_r8 - 0.4_r8 * Bio(i,k,itemp))*     &
     &                     Bio(i,k,iEupO)
              DBio(i,k,iDet) = DBio(i,k,iDet) + cff1 * dtdays
              DBio(i,k,iEupO) = DBio(i,k,iEupO) - cff1 * dtdays


#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iEupS,iDet)= bflx(iEupS,iDet) + cff1* dtdays*xi
                bflx(iEupO,iDet)= bflx(iEupO,iDet) + cff1* dtdays*xi
              END IF
#endif
            END DO
          END DO

!=========================================================================
! Detrital Remineralization   (Det -> NH4)
!=========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend
              Temp1 = Bio(i,k,itemp)
              !-------------------------
              !From Frost (1993).
              !-------------------------
              !cff1 = regen * dgrad * Bio(i,k,iDet)
              !DBio(i,k,iNH4) = DBio(i,k,iNH4) + xi * cff1 * dtdays
              !DBio(i,k,iDet) = DBio(i,k,iDet) - cff1 * dtdays

              !-------------------
              !From Kawamiya(2000)
              !-------------------
              PON = Bio(i,k,iDet)*xi  !Particulate organic nitrogen

              DBio(i,k,iDet) = DBio(i,k,iDet) -                       &
     &                     ((Pv0*exp(PvT*Temp1)*PON)/xi)*dtdays

              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                       &
     &                     ((Pv0*exp(PvT*Temp1)*PON))*dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iDet,iNH4)= bflx(iDet,iNH4) +                      &
     &                    ((Pv0*exp(PvT*Temp1)*PON))*dtdays
              END IF
#endif

              PON = Bio(i,k,iDetF)*xi  !Particulate organic nitrogen

              DBio(i,k,iDetF) = DBio(i,k,iDetF) -                       &
     &                     ((Pv0*exp(PvT*Temp1)*PON)/xi)*dtdays
              DBio(i,k,iNH4) = DBio(i,k,iNH4) +                         &
     &                     ((Pv0*exp(PvT*Temp1)*PON))*dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iDetF,iNH4)= bflx(iDetF,iNH4) +                    &
     &                    ((Pv0*exp(PvT*Temp1)*PON))*dtdays
              END IF
#endif
            END DO
          END DO
!=========================================================================
!Nitrification  (NH4 -> NO3)
!=========================================================================
          DO k=1,N(ng)
            DO i=Istr,Iend

              !-----------------------
              !Temperature dependance
              !-----------------------

              !Kawamiya 2000  NitrMax=Nitr0*exp(KnT*Temp1) - Ken

!             NitrifMax=GetNitrifMaxK(Temp1,KnT,Nitr0)

            ! Arhonditsis NitrMax=Nitr0*exp(-ktntr*(Temp1 - ToptNtr)**2)

              NitrifMax=GetNitrifMaxA(Nitr0, ktntr,Temp1,ToptNtr)

              !No temperaure effects - NitrMax is constant

!             NitrifMax=Nitr0

              !-----------------------
              !Light/Depth dependance
              !-----------------------

              !Fennel
              DLNitrif=GetNitrifLight(Par1,tI0,KI)

              !Denman
!             DLNitrif= (z_wL(k)**10_r8)/( (20_r8**10_r8) +             &
!     &                       z_wL(k)**10_r8)

              !No Depth/Ligt effects
!             DLNitrif= 1.0_r8

              !-----------------
              !Saturation
              !-----------------

              !Arhonditsis
              cff1 = Bio(i,k,iNH4)/(KNH4Nit +Bio(i,k,iNH4))

              !No saturation -ken
!             cff1 =1.0_r8


              DBio(i,k,iNH4) = DBio(i,k,iNH4)  - NitrifMax              &
     &                  * Bio(i,k,iNH4) * DLNitrif * cff1 * dtdays
              DBio(i,k,iNO3) = DBio(i,k,iNO3) + NitrifMax               &
      &                 * Bio(i,k,iNH4) * DLNitrif * cff1 * dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iNH4,iNO3)= bflx(iNH4,iNO3)                        &
     &                 + NitrifMax  * Bio(i,k,iNH4) * DLNitrif *        &
     &                 cff1 * dtdays
              END IF
#endif
            END DO
          END DO

!======================================================================
!Benthic Sub Model
!======================================================================
#   ifdef BENTHIC
          DO k=1,NBL(ng)
            DO i=Istr,Iend
              Temp1 = Bio(i,k,itemp)
!!Potential food available from bottom layer (k=1) of water column

              cff1=(prefD* (Bio(i,k,iDet)+Bio(i,k,iDetF))*Hz(i,j,k)/    &
     &          (prefD*(Bio(i,k,iDet)+Bio(i,k,iDetF)) *Hz(i,j,k)+LupP)) &
     &          *prefD*(Bio(i,k,iDet)+Bio(i,k,iDetF))*Hz(i,j,k)


!              cff1=(prefD* (Bio(i,k,iDet)+Bio(i,k,iDetF))*Hz(i,j,k)/   &
!     &          (prefD*(Bio(i,k,iDet)) *Hz(i,j,k)+LupP))               &
!     &          *prefD*(Bio(i,k,iDet))*Hz(i,j,k)


              cff2=(prefPL*Bio(i,k,iPhL)*Hz(i,j,k)/                     &
     &               (prefPL*Bio(i,k,iPhL)                              &
     &               *Hz(i,j,k)+LupP)) *prefPL*Bio(i,k,iPhL)*Hz(i,j,k)

              cff3=(prefPS*Bio(i,k,iPhS)*Hz(i,j,k)/                     &
     &               (prefPS*Bio(i,k,iPhS)                              &
     &               *Hz(i,j,k)+LupP)) *prefPS*Bio(i,k,iPhS)*Hz(i,j,k)

              cff4=cff1+cff2+cff3
              cff10=BioB(i,k,iBenDet)/(BioB(i,k,iBenDet)+LupD)
              cff5=q10**((Temp1-T0ben)/10.0_r8)
              cff9=(cff5*cff10*BioB(i,k,iBen)*Rup)/(cff10+KupD)

              !--------------------
              !Growth of Benthic Infauna
              !--------------------
              DBioB(i,k,iBenDet)=DBioB(i,k,iBenDet)-dtdays*cff9
              DBioB(i,k,iBen)=DBioB(i,k,iBen)                           &
     &                             + (cff9)*dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(NT(ng)+2,NT(ng)+1)=bflx(NT(ng)+2,NT(ng)+1)         &
     &            +dtdays*cff9*xi
              ENDIF
#endif

#ifdef STATIONARY2
!             Stat2(i,1)=(cff9)*dtdays
#endif
              !------------------
              ! Benthic Production
              !------------------
#ifdef PROD2
              Prod2(i,iBenPrd) = Prod2(i,iBenPrd) +  (cff9)*dtdays
#endif
              cff6=(cff5*cff1*BioB(i,k,iBen)*Rup)/(cff4+KupP)
              cff7=(cff5*cff2*BioB(i,k,iBen)*Rup)/(cff4+KupP)
              cff8=(cff5*cff3*BioB(i,k,iBen)*Rup)/(cff4+KupP)

              !----------------------------------
              !Changes due to benthic consumption
              !----------------------------------
              DBioB(i,k,iBen)=DBioB(i,k,iBen)                           &
     &                             + (cff6 +cff7+cff8)*dtdays

              DBio(i,k,iDet) =DBio(i,k,iDet)- dtdays*cff6/Hz(i,j,k)
              DBio(i,k,iPhL) =DBio(i,k,iPhL)- dtdays*cff7/Hz(i,j,k)
              DBio(i,k,iPhS) =DBio(i,k,iPhS)- dtdays*cff8/Hz(i,j,k)

#ifdef PROD2
              Prod2(i,iBenPrd) = Prod2(i,iBenPrd) +                     &
     &                     (cff6 +cff7+cff8)*dtdays
#endif
#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(iDet,NT(ng)+1)=  bflx(iDet,NT(ng)+1)               &
     &                +dtdays*cff6/Hz(i,j,k)*xi
                bflx(iPhL,NT(ng)+1)=  bflx(iPhL,NT(ng)+1)               &
     &                +dtdays*cff7/Hz(i,j,k)*xi
                bflx(iPhS,NT(ng)+1)=  bflx(iPhS,NT(ng)+1)               &
     &                +dtdays*cff8/Hz(i,j,k)*xi
              END IF
#endif
              !----------------------------
              !Excretion
              !----------------------------

              cff1=cff6*eexD
              cff2=cff7*eex
              cff3=cff8*eex
              cff4=cff9*eexD

              DBioB(i,k,iBen)=DBioB(i,k,iBen)-                          &
     &                   (cff1+cff2+cff3+cff4)*dtdays
              DBioB(i,k,iBenDet)=DBioB(i,k,iBenDet)                     &
     &                    +1.0_r8*dtdays*(cff1+cff2+cff3+cff4)
#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(NT(ng)+1,NT(ng)+2)=bflx(NT(ng)+1,NT(ng)+2)         &
     &                    + (cff1+cff2+cff3+cff4)*dtdays*xi
              ENDIF
#endif
              !---------------
              ! Respiration
              !---------------
              cff5= q10r**((Temp1-T0benr)/10.0_r8)
              cff3=cff5*BioB(i,k,iBen)*Rres
              cff4=Qres*(((1_r8-eexD)*cff6)+((1_r8-eex)*cff7)+          &
     &                 ((1_r8-eex)*cff8)+((1_r8-eexD)*cff9))
              cff6=cff3+cff4

              DBioB(i,k,iBen)=DBioB(i,k,iBen) -cff6*dtdays
              DBio(i,k,iNH4)= DBio(i,k,iNH4)+iremin*xi*dtdays*          &
     &                        cff6/Hz(i,j,k) !includes 20% loss to N2
#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(NT(ng)+1,iNH4)=bflx(NT(ng)+1,iNH4)                 &
     &                        +iremin*xi*dtdays*cff6/Hz(i,j,k)
              ENDIF
#endif
              !------------------
              ! Benthic Production
              !------------------
#ifdef PROD2
!g             Prod2(i,iBenPrd) = Prod2(i,iBenPrd) - cff6*dtdays
#endif
              !-----------
              ! Mortality
              !-----------
              cff1= rmort*BioB(i,k,iBen)
              DBioB(i,k,iBen) = DBioB(i,k,iBen)- cff1 * dtdays
              DBioB(i,k,iBenDet)= DBioB(i,k,iBenDet)+cff1*dtdays

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(NT(ng)+1,NT(ng)+2)=bflx(NT(ng)+1,NT(ng)+2)         &
     &                   + cff1*dtdays*xi
              ENDIF
#endif
              !------------
              ! Predation
              !------------
              DBioB(i,k,iBen) =DBioB(i,k,iBen)                          &
     &                  -BenPred *dtdays*BioB(i,k,iBen)**2

              !-----------------------------------
              !(Det -> NH4) temperature dependent
              !-----------------------------------

              PON = BioB(i,k,iBenDet)*xi/Hz(i,j,k)  !Benthic Particulate organic nitrogen

!              cff5= q10**((Temp1-T0ben)/10.0_r8)
!             cff1=Pv0*cff5*PON
              cff1=Pv0*exp(PvT*Temp1)*PON      !-Kawamiya 2000

              DBioB(i,k,iBenDet) =DBioB(i,k,iBenDet)                    &
     &                    - (cff1/xi)*dtdays*Hz(i,j,k)
              DBio(i,k,iNH4)= DBio(i,k,iNH4)                            &
     &                     +iremin*(cff1)*dtdays  !includes 20% loss to N2

#if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
                bflx(NT(ng)+2,iNH4)=bflx(NT(ng)+2,iNH4)                 &
     &                    + iremin*(cff1)*dtdays
              ENDIF
#endif

            END DO
          END DO
#   endif
!
!==================================================================
! Vertical Sinking of Particles(Phyto and Det)
!==================================================================

       !
       ! Use Sinking Code adapted from J. Warner ROMS sediment code
       !   - this is similar to the approach taken in Fasham
       !     biology code but adapted to be used in subroutine
       ! -Also addapted to return particle flux in and out of each
       !   level - used in BENTHIC sub model
       ! -Incorporated Liz Dobbins zlimit to determine sinking rate
       ! zlimit =1 for constant sinking - use for Phyto and Det
       ! zlimit =-1*NcmaxZ for Neocalanus- sink rate is then attenuated
       ! as max sinking depth is approached
       !
       !        G.Gibson  July 2008
       !
          CALL BIOSINK_1(ng,wPhS,Bio(IminS,1,iPhS),sinkIN,              &
     &            sinkOUT,HzL,dtdays,z_wL,1.0_r8,LBi,UBi,IminS, ImaxS)

          DO k=1,N(ng)
            DO i=Istr,Iend
              DBio(i,k,iPhS) =DBio(i,k,iPhS) +  (sinkIN(1,k)            &
     &                   -sinkOUT(1,k))/Hz(i,j,k)
#ifdef BENTHIC
              IF (k.eq.1) THEN
                DBioB(i,k,iBenDet)=DBioB(i,k,iBenDet)+sinkOUT(i,k)
                Stat2(i,1)=Stat2(i,1)+sinkOUT(i,k)
              END IF
#endif
            END DO
          END DO

          CALL BIOSINK_1(ng,wPhL,Bio(IminS,1,iPhL),sinkIN,              &
     &             sinkOUT,HzL,dtdays,z_wL,1.0_r8,LBi,UBi,IminS, ImaxS)

          DO k=1,N(ng)
            DO i=Istr,Iend
              DBio(i,k,iPhL) =DBio(i,k,iPhL) +  (sinkIN(1,k)            &
     &                -sinkOUT(1,k))/Hz(i,j,k)
#ifdef BENTHIC
              IF (k.eq.1) THEN
                DBioB(i,k,iBenDet)=DBioB(i,k,iBenDet)+sinkOUT(i,k)
                      Stat2(i,1)=Stat2(i,1)+sinkOUT(i,k)
              END IF
#endif
            END DO
          END DO

          CALL BIOSINK_1(ng,wDet,Bio(IminS,1,iDet),sinkIN,              &
     &             sinkOUT,HzL,dtdays,z_wL,1.0_r8,LBi,UBi,IminS, ImaxS)

          DO k=1,N(ng)
            DO i=Istr,Iend
              DBio(i,k,iDet) =DBio(i,k,iDet) +  (sinkIN(1,k)            &
     &                -sinkOUT(1,k))/Hz(i,j,k)
#ifdef BENTHIC
              IF (k.eq.1) THEN
                DBioB(i,k,iBenDet)=DBioB(i,k,iBenDet)+sinkOUT(i,k)
                      Stat2(i,1)=Stat2(i,1)+sinkOUT(i,k)
              END IF
#endif
            END DO
          END DO

          call BIOSINK_1(ng,wDetF,Bio(IminS,1,iDetF),sinkIN,            &
     &              sinkOUT,HzL,dtdays,z_wL,1.0_r8,LBi,UBi,IminS, ImaxS)


          DO k=1,N(ng)
            DO i=Istr,Iend
              DBio(i,k,iDetF) =DBio(i,k,iDetF) +  (sinkIN(1,k)          &
     &                  -sinkOUT(1,k))/Hz(i,j,k)
#ifdef BENTHIC
              IF (k.eq.1) THEN
                DBioB(i,k,iBenDet)=DBioB(i,k,iBenDet)+sinkOUT(i,k)
                Stat2(i,1)=Stat2(i,1)+sinkOUT(i,k)
              END IF
#endif
            END DO
          END DO

       ! Alternate sinking code -Ken
       ! Sinking code developed by C. Lewis and E. Dobbins
       ! Not sure where this came from
       !

!          call BIOSINK(ng,wPhS,Bio(LBi-1,1,iPhS),DBio(LBi-1,1,iPhS),    &
!     &                  HzL,dtdays,z_wL,1.0_r8,LBi,UBi)


!          call BIOSINK(ng,wPhL,Bio(LBi-1,1,iPhL),DBio(LBi-1,1,iPhL),    &
!     &                  HzL,dtdays,z_wL,1.0_r8,LBi,UBi)

             !------------------------------------------------------
             !ROMS seems to bomb with shallow water. Use of the
             !DetSINK option seems to have problems in shallow water
             !See if it will run OK for water deeper than 60
             ! (this check now happens inside DetSINK)
             !------------------------------------------------------
!          call DetSINK(ng,wDet,Bio(LBi-1,1,iDet),DBio(LBi-1,1,iDet),  &
!     &                       HzL,dtdays,z_wL,1.0_r8,Dens,LBi,UBi)


!=======================================================================
! Large Copepod Diapause
!=======================================================================

#ifdef DIAPAUSE
          IF ( downward ) THEN

          !new sinking code adapted from J. Warner+ Dobbins zlimit
          !G. Gibson July 2008

            call BIOSINK_1(ng,wNCsink,Bio(IminS,1,iNCa),sinkIN,         &
     &              sinkOUT,HzL,dtdays,z_wL,1.0_r8,LBi,UBi,IminS, ImaxS)

            DO  k=1,N(ng)
              DO i=Istr,Iend
                DBio(i,k,iNCa) =DBio(i,k,iNCa) +  (sinkIN(1,k)          &
     &                    -sinkOUT(1,k))/Hz(i,j,k)
              END DO
            END DO

!         original sinking code from Lewis and Dobbins
!
!             call BIOSINK_2(ng,wNCsink,Bio(LBi-1,1,iNCa),              &
!     &                DBio(LBi-1,1,iNCa),HzL,dtdays,z_wL,-1*NCmaxz,    &
!     &                LBi,UBi)

          ELSE IF (upward) THEN

        !New rising code adapted from J. Warner sink code + Dobbins zlimit
        !G. Gibson July 2008

            call BIORISE_1(ng,wNCsink,Bio(LBi-1,1,iNCa),riseOUT,        &
     &                riseIN,HzL,dtdays,z_wL,LBi,UBi,-1*NCmaxz,         &
     &                IminS, ImaxS)
            DO k=1,N(ng)
              DO i=Istr,Iend
                DBio(i,k,iNCa) =DBio(i,k,iNCa) +  (riseIN(1,k)          &
     &                   -riseOUT(1,k))/Hz(i,j,k)
              END DO
            END DO

!         original rising code from Lewis and Dobbins
!             call BIORISE(ng,wNCrise,Bio(LBi-1,1,iNCa),                &
!     &              DBio(LBi-1,1,iNCa),HzL,dtdays,z_wL,-1*NCmaxz,      &
!     &              LBi,UBi,6.0_r8)
          END IF
#endif

!==================================================
!Ice Sub Model
!==================================================
#ifdef ICE_BIO

          DO i=Istr,Iend
            IF (ice_thick(i,j).ge.0.02)THEN
!ICE CAN GROW
# if defined CLIM_ICE_1D
              Temp1 = Bio(i,N(ng),itemp)
# endif
!              Temp1 = ti(i,j,1)

!need to correct to use PAR in middleof light layer.maby not algae grow under the ice
              Par1 = Par(i,42)
              cff1=BioBI(i,iIceNO3)/(ksnut1+BioBI(i,iIceNO3))
              cff2=BioBI(i,iIceNH4)/(ksnut2+BioBI(i,iIceNH4))
              aiceNfrac=exp(-inhib*BioBI(i,iIceNH4))
              aiceIfrac=(1-exp(-alphaIb*Par1))*exp(-betaI*Par1)

!had to cap gesi to prevent from going negative
!when gesi is defined I find that Ice Algae does not grow

              sb=-3.9921-22.7* Temp1-1.0015*Temp1**2-0.02* Temp1**3
              gesi=max(0._r8,1.1e-2+3.012e-2*sb+1.0342e-3*sb*sb)

! growth of Ice Algae

!              grow1=mu0*exp(0.0633*Temp1)*gesi
              grow1=mu0*exp(0.0633*Temp1)
              NOup=cff1*min(aiceIfrac,aiceNfrac)
              NHup=cff2*aiceIfrac
              GROWAice=grow1*(NOup+NHup)

              DBioBI(i,iIcePhL)=DBioBI(i,iIcePhL)                       &
     &                   + GROWAice*BioBI(i,iIcePhL)* dtdays

              !-----------------------------------------
              !Primary production of ice algae
              !-----------------------------------------
# ifdef PROD2
              Prod2(i,iIAPrd) = Prod2(i,iIAPrd)                         &
     &                      + GROWAice*BioBI(i,iIcePhL)* dtdays
# endif

! respiration of Ice Algae
!              RAi0=R0i*grow1
              RAi0=R0i*GROWAice

              DBioBI(i,iIcePhL)=DBioBI(i,iIcePhL)                       &
     &                   -BioBI(i,iIcePhL)*RAi0*dtdays

! mortality of Ice Algae
              RgAi=rg0*exp(rg*Temp1)
              DBioBI(i,iIcePhL)=DBioBI(i,iIcePhL)                       &
     &                -BioBI(i,iIcePhL)*RgAi*dtdays

              reN=annit*BioBI(i,iIceNH4)

              cff2=BioBI(i,iIcePhL)*(GROWAice-RAi0)
              cff3=BioBI(i,iIcePhL)*RAi0-reN

              DBioBI(i,iIceNO3)=DBioBI(i,iIceNO3)                       &
     &                  -grow1*NOup*BioBI(i,iIcePhL)*xi*dtdays          &
     &                  +reN*dtdays

              DBioBI(i,iIceNH4)=DBioBI(i,iIceNH4)                       &
     &                     -grow1*NHup*BioBI(i,iIcePhL)*xi*dtdays       &
     &               +RAi0*BioBI(i,iIcePhL)*xi*dtdays                   &
     &                     +RgAi*BioBI(i,iIcePhL)*xi*dtdays             &
     &               -reN*dtdays

# if defined BIOFLUX && defined BIO_GOANPZ
              IF (i.eq.3.and.j.eq.3) THEN
!PLi->NH4
                bflx(NT(ng)+3,NT(ng)+5)=bflx(NT(ng)+3,NT(ng)+5)         &
     &           +RAi0*BioBI(i,iIcePhL)*xi*dtdays                       &
     &            +RgAi*BioBI(i,iIcePhL)*xi*dtdays

!NH4->NO3
                bflx(NT(ng)+5,NT(ng)+4)=bflx(NT(ng)+5,NT(ng)+4)         &
     &            +reN*dtdays

!NO3->iPL
                bflx(NT(ng)+4,NT(ng)+3)=bflx(NT(ng)+4,NT(ng)+3)         &
     &            +grow1*NOup*BioBI(i,iIcePhL)*xi*dtdays

!NH4->iPL
                bflx(NT(ng)+5,NT(ng)+3)=bflx(NT(ng)+5,NT(ng)+3)         &
     &            +grow1*NHup*BioBI(i,iIcePhL)*xi*dtdays
              ENDIF
# endif


!          st2(i,j,nnew,i2Stat3)=reN*dtdays
!          st2(i,j,nnew,i2Stat4)=annit*BioBI(i,iIceNH4)
!          st2(i,j,nnew,i2Stat5)=BioBI(i,iIceNH4)

# if defined CLIM_ICE_1D

!             dhicedt=tclmG(i,j,42,1,15) - tclmG(i,j,42,2,15)
              dhicedt=it(i,j,nnew,iIceZ)-it(i,j,nstp,iIceZ)
# else

              dhicedt=it(i,j,nnew,iIceZ)-it(i,j,nstp,iIceZ)
!              dhicedt=(hi(i,j,nstp)-hi(i,j,nnew))
# endif
              dhicedt=dhicedt*sec2day/dtdays !convert to m/s

              trs=9.667e-11+4.49e-6*dhicedt-1.39e-5*dhicedt**2
              trs=trs*86400   !convert to m/d
              twi=72*trs

              IF (dhicedt.gt.0) THEN
                trs=4.49e-6*dhicedt-1.39e-5*dhicedt**2
                trs=trs*86400
                twi=720*trs
              ENDIF

              !----------------------------------------------
              !Change in concentration of large phytoplankton/ Ice Algae
              !due to ice growth/melt
              !----------------------------------------------

              IF (twi.lt.0) THEN
                DBioBI(i,iIcePhL)=DBioBI(i,iIcePhL)                     &
     &                + BioBI(i,iIcePhL)*(twi/aidz)*                    &
     &                            86400*dtdays
                DBio(i,N(ng),iPhL) = DBio(i,N(ng),iPhL)                 &
     &            -twi*BioBI(i,iIcePhL)*dtdays*86400/                   &
     &                            Hz(i,j,N(ng))

# if defined BIOFLUX && defined BIO_GOANPZ
                IF (i.eq.3.and.j.eq.3) THEN
!PL->Phi
                  bflx(iPhL,NT(ng)+3)=bflx(iPhL,NT(ng)+3)               &
     &                +twi*BioBI(i,iIcePhL)*dtdays*86400/               &
     &                            Hz(i,j,N(ng)) *xi
                ENDIF
# endif

              ENDIF
! nutrient gradient between ice and water

              cff1=twi*(Bio(i,N(ng),iNO3)-BioBI(i,iIceNO3))*dtdays
              cff2=twi*(Bio(i,N(ng),iNH4)-BioBI(i,iIceNH4))*dtdays

              IF (twi.lt.0) THEN
                DBioBI(i,iIceNO3)=DBioBI(i,iIceNO3)                     &
     &                            + BioBI(i,iIceNO3)*twi*dtdays
                DBioBI(i,iIceNH4)=DBioBI(i,iIceNH4)                     &
     &                            + BioBI(i,iIceNH4)*twi*dtdays
                DBio(i,N(ng),iNO3) = DBio(i,N(ng),iNO3)+cff1/           &
     &                            Hz(i,j,N(ng))
                DBio(i,N(ng),iNH4) = DBio(i,N(ng),iNH4)+cff2/           &
     &                            Hz(i,j,N(ng))

# if defined BIOFLUX && defined BIO_GOANPZ
                IF (i.eq.3.and.j.eq.3) THEN
!NH4i->NH4
                  bflx(NT(ng)+4,iNO3)=bflx(NT(ng)+4,iNO3)               &
     &                  +  cff1/Hz(i,j,N(ng))
!NH4i->NH4
                  bflx(NT(ng)+5,iNH4)=bflx(NT(ng)+5,iNH4)               &
     &                  +  cff2/Hz(i,j,N(ng))
                ENDIF
# endif

              ELSE IF (twi.gt.0) THEN
                DBioBI(i,iIceNO3)=DBioBI(i,iIceNO3)+cff1/aidz
                DBioBI(i,iIceNH4)=DBioBI(i,iIceNO3)+cff2/aidz
                DBio(i,N(ng),iNO3) = DBio(i,N(ng),iNO3)-cff1/           &
     &                               Hz(i,j,N(ng))
                DBio(i,N(ng),iNH4) = DBio(i,N(ng),iNH4)-cff2/           &
     &                               Hz(i,j,N(ng))
# if defined BIOFLUX && defined BIO_GOANPZ
                IF (i.eq.3.and.j.eq.3) THEN
!NO3->NO3i
                  bflx(iNO3,NT(ng)+4)=bflx(iNO3,NT(ng)+4)               &
     &                  +cff1/aidz
!NH4->NH4i
                  bflx(iNH4,NT(ng)+5)=bflx(iNH4,NT(ng)+5)               &
     &                  +  cff2/aidz
                ENDIF
# endif

              ENDIF
# ifdef STATIONARY2
!        st2(i,j,nnew,i2Stat1)=trs
!        st2(i,j,nnew,i2Stat2)=twi


!        st2(i,j,nnew,i2Stat5)=it(i,j,nnew,iIceZ)
!        st2(i,j,nnew,i2Stat6)=it(i,j,nstp,iIceZ)
!        st2(i,j,nnew,i2Stat7)=tclmG(i,j,42,1,15)
!        st2(i,j,nnew,i2Stat8)=tclmG(i,j,42,2,15)
# endif

!              ENDIF

            ELSE IF (ice_thick(i,j).lt.0.02)      THEN
# if defined ICE_BIO
              IF (itL(i,j,nstp,iIceLog).eq.1_r8 ) THEN
                DBio(i,N(ng),iPhL) = DBio(i,N(ng),iPhL)                 &
     &                     + it(i,j,nnew,iIcePhL)*aidz/Hz(i,j,N(ng))
                DBio(i,N(ng),iNO3) = DBio(i,N(ng),iNO3)                 &
     &                     + it(i,j,nnew,iIceNO3)*aidz/Hz(i,j,N(ng))
                DBio(i,N(ng),iNH4) = DBio(i,N(ng),iNH4)                 &
     &                     + it(i,j,nnew,iIceNH4)*aidz/Hz(i,j,N(ng))
                itL(i,j,nnew,iIceLog) =-1_r8
#  if defined BIOFLUX && defined BIO_GOANPZ
                IF (i.eq.3.and.j.eq.3) THEN
!PLi->PhL
                  bflx(NT(ng)+3,iPhL)=bflx(NT(ng)+3,iPhL)               &
     &              + it(i,j,nnew,iIcePhL)*aidz/Hz(i,j,N(ng))*xi

!NO3i->NO3
                  bflx(NT(ng)+4,iNO3)=bflx(NT(ng)+4,iNO3)               &
     &              + it(i,j,nnew,iIceNO3)*aidz/Hz(i,j,N(ng))

!NH4i->NH4
                  bflx(NT(ng)+5,iNH4)=bflx(NT(ng)+5,iNH4)               &
     &              + it(i,j,nnew,iIceNH4)*aidz/Hz(i,j,N(ng))
                ENDIF
#  endif
                DO itrc=1,3 !NIB
                  ibioBI=idice(itrc)
                  DBioBI(i,ibioBI)=DBioBI(i,ibioBI)-BioBI(i,ibioBI)
                  it(i,j,nnew,ibioBI)=0_r8
                END DO
              ELSE
                DO itrc=1,3 !NIB
                  ibioBI=idice(itrc)
                  DBioBI(i,ibioBI)=0_r8
                END DO

!                DBio(i,N(ng),iPhL)=0_r8
!                DBio(i,N(ng),iNO3)=0_r8
!                DBio(i,N(ng),iNH4)=0_r8
              END IF
# endif
#endif
            END IF
          END DO


!=======================================================================
! Update Bio array
!=======================================================================

          DO itrc=1,NBT
            ibio=idbio(itrc)
            DO k=1,N(ng)
              DO i=Istr,Iend
!               DBio(i,k,6)=1
                Bio(i,k,ibio)=Bio(i,k,ibio)+DBio(i,k,ibio)
                IF (Bio(i,k,iNO3).gt.35) THEN
                  print*,'NO3=',Bio(i,k,iNO3)
                  print*,'i=',i,'j=',j,'k=',k
                END IF
              END DO
            END DO
          END DO
#ifdef BENTHIC
          DO itrc=1,NBEN
            ibioB=idben(itrc)
            DO k=1,NBL(ng)
              DO i=Istr,Iend
                BioB(i,k,ibioB)=BioB(i,k,ibioB)+DBioB(i,k,ibioB)
                BioB(i,k,ibioB)=BioB(i,k,ibioB)+dtdays
              END DO
            END DO
          END DO

#endif
#ifdef ICE_BIO
          DO itrc=1,3
            ibioBI=idice(itrc)
            DO i=Istr,Iend

              BioBI(i,ibioBI)=BioBI(i,ibioBI)+DBioBI(i,ibioBI)

            END DO
          END DO
#endif
        END DO ITER_LOOP

!=======================================================================
! Update global tracer variables (m Tunits).
!=======================================================================

        DO itrc=1,NBT
          ibio=idbio(itrc)
          DO k=1,N(ng)
            DO i=Istr,Iend
              t(i,j,k,nnew,ibio)=MAX(t(i,j,k,nnew,ibio)+                &
     &                               (Bio(i,k,ibio)-Bio_bak(i,k,ibio))  &
     &                               *Hz(i,j,k)                         &
     &                               ,0.0_r8)
            END DO
          END DO
        END DO
#   ifdef BENTHIC

        DO itrc=1,NBEN
          ibioB=idben(itrc)
          DO k=1,NBL(ng)
            DO i=Istr,Iend
              bt(i,j,k,nnew,ibioB)=MAX(bt(i,j,k,nnew,ibioB)+            &
     &                           (BioB(i,k,ibioB)-Bio_bakB(i,k,ibioB))  &
     &                               ,0.0_r8)
              bt(i,j,k,nstp,ibioB)= bt(i,j,k,nnew,ibioB)
            END DO
          END DO
        END DO

#   endif
#   ifdef ICE_BIO

        DO itrc=1,3 !NIB
          ibioBI=idice(itrc)
          DO i=Istr,Iend
!            IF (ice_thick(i,j).ge.0.02)THEN

            it(i,j,nnew,ibioBI)=MAX(it(i,j,nnew,ibioBI)+                &
     &                           BioBI(i,ibioBI)-Bio_bakBI(i,ibioBI)    &
     &                               ,0.0_r8)


!            ELSE
!                 it(i,j,nnew,ibioBI)=0_r8
!            END IF

          END DO
        END DO
#endif

#ifdef STATIONARY
        DO k=1,N(ng)
!         DO itrc=1,UBst
!          ibio=idbio3(itrc)
          DO i=Istr,Iend
!                  st(i,j,k,nstp, ibio)=0.0_r8
!               st(i,j,k,nstp, ibio)= Stat3(i,k,ibio)
!              END DO
            st(i,j,k,nstp,1) =    Stat3(i,k,1)
            st(i,j,k,nstp,2) =    Stat3(i,k,2)
            st(i,j,k,nstp,3) =    Stat3(i,k,3)
            st(i,j,k,nstp,4) =    Stat3(i,k,4)
            st(i,j,k,nstp,5) =    Stat3(i,k,5)
            st(i,j,k,nstp,6) =    Stat3(i,k,6)
            st(i,j,k,nstp,7) =    Stat3(i,k,7)
            st(i,j,k,nstp,8) =    Stat3(i,k,8)
            st(i,j,k,nstp,9) =    Stat3(i,k,9)
            st(i,j,k,nstp,10) =    Stat3(i,k,10)
            st(i,j,k,nstp,11) =    Stat3(i,k,11)
            st(i,j,k,nstp,12) =    Stat3(i,k,12)
            st(i,j,k,nstp,13) =    Stat3(i,k,13)
            st(i,j,k,nstp,14) =    Stat3(i,k,14)
            st(i,j,k,nstp,15) =    Stat3(i,k,15)
            st(i,j,k,nstp,16) =    Stat3(i,k,16)


!            st(i,j,k,nnew, i3Stat2) = st(i,j,k,nstp, i3Stat2)
!            st(i,j,k,nnew, i3Stat3) = st(i,j,k,nstp, i3Stat3)
!            st(i,j,k,nnew, i3Stat4) = st(i,j,k,nstp, i3Stat4)
!            st(i,j,k,nnew, i3Stat5) = st(i,j,k,nstp, i3Stat5)
!            st(i,j,k,nnew, i3Stat6) = st(i,j,k,nstp, i3Stat6)
!            st(i,j,k,nnew, i3Stat7) = st(i,j,k,nstp, i3Stat7)
!            st(i,j,k,nnew, i3Stat8) = st(i,j,k,nstp, i3Stat8)

          END DO
        END DO
#endif
#ifdef STATIONARY2

        DO i=Istr,Iend
          st2(i,j,nstp,1) =   Stat2(i,1)
!          print*,'Stat2=',Stat2(i,1)
          st2(i,j,nstp,2) =   Stat2(i,2)
          st2(i,j,nstp,3) =   Stat2(i,3)
          st2(i,j,nstp,4) =   Stat2(i,4)
          st2(i,j,nstp,5) =   Stat2(i,5)
          st2(i,j,nstp,6) =   Stat2(i,6)
          st2(i,j,nstp,7) =   Stat2(i,7)
          st2(i,j,nstp,8) =   Stat2(i,8)

        END DO
#endif
#ifdef PROD3
        DO k=1,N(ng)
          DO i=Istr,Iend
            pt3(i,j,k,nnew,iPhSprd) = pt3(i,j,k,nnew,iPhSprd) +           &
     &                             Prod(i,k,iPhS)
            pt3(i,j,k,nnew,iPhLprd) = pt3(i,j,k,nnew,iPhLprd) +           &
     &                             Prod(i,k,iPhL)
            pt3(i,j,k,nnew,iMZSprd) = pt3(i,j,k,nnew,iMZSprd) +           &
     &                             Prod(i,k,iMZS)
            pt3(i,j,k,nnew,iMZLprd) = pt3(i,j,k,nnew,iMZLprd) +           &
     &                             Prod(i,k,iMZL)
            pt3(i,j,k,nnew,iCopPrd) = pt3(i,j,k,nnew,iCopPrd) +           &
     &                             Prod(i,k,iCop)
            pt3(i,j,k,nnew,iNCaPrd) = pt3(i,j,k,nnew,iNCaPrd) +           &
     &                             Prod(i,k,iNCaS)+Prod(i,k,iNCaO)
            pt3(i,j,k,nnew,iEupPrd) = pt3(i,j,k,nnew,iEupPrd) +           &
     &                             Prod(i,k,iEupS)+Prod(i,k,iEupO)
#ifdef JELLY
            pt3(i,j,k,nnew,iJelPrd) = pt3(i,j,k,nnew,iJelPrd) +           &
     &                             Prod(i,k,iJel)
#endif
          END DO
        END DO
#endif
#ifdef PROD2

        DO i=Istr,Iend
#ifdef BENTHIC
          pt2(i,j,nnew,iBenPrd) = pt2(i,j,nnew,iBenPrd)+Prod2(i,iBenPrd)
#endif
          pt2(i,j,nnew,iIAPrd) = pt2(i,j,nnew,iIAPrd)+Prod2(i,iIAPrd)
        END DO
#endif

      END DO J_LOOP

      RETURN
      END SUBROUTINE biology_tile
 !=====================================================================
! BIOSINK_1  particle sinking subroutine After J. Warner sed sink code
! G. Gibson July 2008
!=====================================================================
      subroutine BIOSINK_1(ng,wBio,Bio,sinkIN,sinkOUT,HzL,dtdays,z_wL,  &
     &                       zlimit,LBi,UBi,IminS, ImaxS)

      USE mod_param
!
      implicit none
!

      integer, intent(in)     :: ng, LBi, UBi, IminS, ImaxS
      real(r8), intent(in)    :: wBio
      real(r8), intent(in)    :: zlimit

      real(r8), intent(in) :: z_wL(IminS:ImaxS,0:N(ng))
      real(r8), intent(inout) :: Bio(IminS:ImaxS,N(ng))
      real(r8), intent(in) :: HzL(IminS:ImaxS,N(ng))

      real(r8), intent(in)    :: dtdays


      real(r8), intent(out) :: sinkIN(UBi,N(ng)),sinkOUT(UBi,N(ng))

      integer :: i,k,ks
      real(r8) :: aL, aR, cff1, cff2
      real(r8) :: cffL, cffR, cu, dltL, dltR,cff


      real(r8):: dBio(0:N(ng)), wBiod(LBi:UBi,0:N(ng))
      real(r8) :: FC(IminS:ImaxS,0:N(ng))

      real(r8) :: Hz_inv(IminS:ImaxS,N(ng))
      real(r8) :: Hz_inv2(IminS:ImaxS,N(ng))
      real(r8) :: Hz_inv3(IminS:ImaxS,N(ng))


      integer :: ksource(IminS:ImaxS,N(ng))

      real(r8) :: qR(IminS:ImaxS,N(ng))
      real(r8) :: qL(IminS:ImaxS,N(ng))
      real(r8) :: WL(IminS:ImaxS,N(ng))
      real(r8) :: WR(IminS:ImaxS,N(ng))
      real(r8), dimension(IminS:ImaxS,N(ng)) :: qc

      IF ( zlimit .lt. 0 ) THEN
        DO k=0,N(ng)
          DO i=LBi, UBi
            IF ( z_wL(i,k) .ge. zlimit ) THEN
              wBiod(i,k) = wBio*exp( -1*(z_wL(i,k)-(zlimit/2))**2/   &
     &          (zlimit/2)**2 )
            ELSE
              wBiod(i,k) = 0.0_r8
            END IF
          END DO
        END DO
      ELSE
        DO k=0,N(ng)
          DO i=LBi, UBi
            wBiod(i,k) = wBio
          END DO
        END DO
      END IF
!
!
!  Compute inverse thickness to avoid repeated divisions.
!

      DO k=1,N(ng)
        DO i=LBi,UBi
          Hz_inv(i,k)=1.0_r8/HzL(i,k)
        END DO
      END DO
      DO k=1,N(ng)-1
        DO i=LBi,UBi
          Hz_inv2(i,k)=1.0_r8/(HzL(i,k)+HzL(i,k+1))
        END DO
      END DO
      DO k=2,N(ng)-1
        DO i=LBi,UBi
          Hz_inv3(i,k)=1.0_r8/(HzL(i,k-1)+HzL(i,k)+HzL(i,k+1))
        END DO
      END DO

      DO k=1,N(ng)
        DO i=LBi,UBi
          qc(i,k)=Bio(i,k)
        END DO
      END DO

!
!  Reconstruct vertical profile of suspended sediment "qc" in terms
!  of a set of parabolic segments within each grid box. Then, compute
!  semi-Lagrangian flux due to sinking.
!
      DO k=N(ng)-1,1,-1
        DO i=LBi,UBi
          FC(i,k)=(qc(i,k+1)-qc(i,k))*Hz_inv2(i,k)
!r         FC(i,k)=2_r8
        END DO
      END DO

      DO k=2,N(ng)-1
        DO i=LBi,UBi
          dltR=HzL(i,k)*FC(i,k)
          dltL=HzL(i,k)*FC(i,k-1)
          cff=HzL(i,k-1)+2.0_r8*HzL(i,k)+HzL(i,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 (qR,qL) 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 qL(k+1)-qR(k) may still have different sign than
!        qc(k+1)-qc(k).  This possibility is excluded, after qL and qR
!        are reconciled using WENO procedure.
!
          cff=(dltR-dltL)*Hz_inv3(i,k)
          dltR=dltR-cff*HzL(i,k+1)
          dltL=dltL+cff*HzL(i,k-1)
          qR(i,k)=qc(i,k)+dltR
          qL(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=LBi,UBi
          dltL=MAX(cff,WL(i,k  ))
          dltR=MAX(cff,WR(i,k+1))
          qR(i,k)=(dltR*qR(i,k)+dltL*qL(i,k+1))/(dltR+dltL)
          qL(i,k+1)=qR(i,k)
        END DO
      END DO

      DO i=LBi,UBi
        FC(i,N(ng))=0.0_r8              ! no-flux boundary condition
#  if defined LINEAR_CONTINUATION
        qL(i,N(ng))=qR(i,N(ng)-1)
        qR(i,N(ng))=2.0_r8*qc(i,N(ng))-qL(i,N(ng))
#  elif defined NEUMANN
        qL(i,N(ng))=qR(i,N(ng)-1)
        qR(i,N(ng))=1.5_r8*qc(i,N(ng))-0.5_r8*qL(i,N(ng))
#  else
        qR(i,N(ng))=qc(i,N(ng))         ! default strictly monotonic
        qL(i,N(ng))=qc(i,N(ng))         ! conditions
        qR(i,N(ng)-1)=qc(i,N(ng))
#  endif
#  if defined LINEAR_CONTINUATION
        qR(i,1)=qL(i,2)
        qL(i,1)=2.0_r8*qc(i,1)-qR(i,1)
#  elif defined NEUMANN
        qR(i,1)=qL(i,2)
        qL(i,1)=1.5_r8*qc(i,1)-0.5_r8*qR(i,1)
#  else
        qL(i,2)=qc(i,1)                 ! bottom grid boxes are
        qR(i,1)=qc(i,1)                 ! re-assumed to be
        qL(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=LBi,UBi
          dltR=qR(i,k)-qc(i,k)
          dltL=qc(i,k)-qL(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
          qR(i,k)=qc(i,k)+dltR
          qL(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.
!

      DO k=1,N(ng)
        DO i=LBi,UBi
          cff=dtdays*ABS(wBiod(i,k))
          FC(i,k-1)=0.0_r8
          WL(i,k)=z_wL(i,k-1)+cff
          WR(i,k)=HzL(i,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=LBi,UBi
            IF (WL(i,k).gt.z_wL(i,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=LBi,UBi
          ks=ksource(i,k)
          cu=MIN(1.0_r8,(WL(i,k)-z_wL(i,ks-1))*Hz_inv(i,ks))
          FC(i,k-1)=FC(i,k-1)+                                          &
     &                  HzL(i,ks)*cu*                                   &
     &                  (qL(i,ks)+                                      &
     &                   cu*(0.5_r8*(qR(i,ks)-qL(i,ks))-                &
     &                       (1.5_r8-cu)*                               &
     &                       (qR(i,ks)+qL(i,ks)-2.0_r8*qc(i,ks))))

            !G.Gibson      - FC is the flux into the level
            !              - should be 0 at the surface
          IF (k.eq.N(ng)) THEN
            FC(i,k)=0.0_r8
          END IF
        END DO
      END DO

      DO k=1,N(ng)
        DO i=LBi,UBi
!            The Bio variables are now updated in the main subroutine
!          Bio(i,k)=qc(i,k)+(FC(i,k)-FC(i,k-1))*Hz_inv(i,k)
          sinkIN(i,k)=FC(i,k)
          sinkOUT(i,k)=FC(i,k-1)
        END DO
      END DO

      RETURN
      END SUBROUTINE BIOSINK_1
!=====================================================================
! BIORISE_1  rising particle subroutine a reversal of J. Warner sed sink code
!      plus and attenuation of rise rate based on closeness to max sink depth
! G.Gibson July 2008
!=====================================================================
    subroutine BIORISE_1(ng,wBio,Bio,riseOUT,riseIN,HzL,dtdays,z_wL,  &
     &                       LBi,UBi,zlimit,IminS, ImaxS)

      USE mod_param
!
      implicit none
!
!      real(r8), intent(in)    :: z_w(0:N(ng))
      integer, intent(in)     :: ng, LBi, UBi, IminS, ImaxS
      real(r8), intent(in)    :: wBio
      real(r8), intent(in)    :: zlimit
!      real(r8), intent(inout) :: Bio(LBi:UBi,N(ng))
      real(r8), intent(in) :: z_wL(IminS:ImaxS,N(ng))
      real(r8), intent(inout) :: Bio(IminS:ImaxS,N(ng))
      real(r8), intent(in) :: HzL(IminS:ImaxS,N(ng))
!      real(r8), intent(in)    :: HzL(N(ng))
      real(r8), intent(in)    :: dtdays



      real(r8), intent(out) :: riseIN(UBi,N(ng)),riseOUT(UBi,N(ng))
!      real(r8) ::DensWk(1:N(ng))
      integer :: i,k,ks
      real(r8) :: aL, aR, cff1, cff2
      real(r8) :: cffL, cffR, cu, dltL, dltR,cff

!      real(r8):: FC(0:N(ng)), dBio(0:N(ng)),
      real(r8):: dBio(0:N(ng)), wBiod(LBi:UBi,N(ng))
      real(r8) :: FC(IminS:ImaxS,N(ng))

      real(r8) :: Hz_inv(IminS:ImaxS,N(ng))
      real(r8) :: Hz_inv2(IminS:ImaxS,N(ng))
      real(r8) :: Hz_inv3(IminS:ImaxS,N(ng))

!      real(r8):: Hz_inv(0:N(ng))
      integer :: ksource(IminS:ImaxS,N(ng))

!      real(r8):: Hz_inv3(N(ng))
!      real(r8):: Hz_inv2(N(ng))
!      real(r8):: WL(N(ng))
!      real(r8):: WR(N(ng))

      real(r8) :: qR(IminS:ImaxS,N(ng))
      real(r8) :: qL(IminS:ImaxS,N(ng))
      real(r8) :: WL(IminS:ImaxS,N(ng))
      real(r8) :: WR(IminS:ImaxS,N(ng))

!      real(r8):: qL(N(ng))
!      real(r8):: qR(N(ng))
!     real(r8):: qR, qL
      real(r8), dimension(IminS:ImaxS,N(ng)) :: qc



      IF ( zlimit .lt. 0 ) THEN
        DO k=1,N(ng)
          DO i=LBi,UBi
            IF ( z_wL(i,k-1).ge.zlimit/2) THEN
              wBiod(i,k) = wBio*exp( -1*(z_wL(i,k-1)-(zlimit/2))**2 /   &
     &          (zlimit/2)**2 )
            ELSE
              wBiod(i,k) = wBio
            END IF
! This check used to be outside the whole function call.
!           IF ( Bio(i,N(ng)).ge.dlimit) THEN
!             wBiod(i,k) = 0.0_r8
!           END IF
          END DO
        END DO
        DO i=LBi,UBi
          wBiod(i,N(ng)+1) = 0.0_r8
        END DO
      ELSE
        DO k=1,N(ng)+1
          DO i=LBi,UBi
            wBiod(i,k) = wBio
          END DO
        END DO
      END IF
!
!  Compute inverse thickness to avoid repeated divisions.
!
      DO k=1,N(ng)
        DO i=LBi,UBi
          Hz_inv(i,k)=1.0_r8/HzL(i,k)
        END DO
      END DO

      DO k=2,N(ng)+1
        DO i=LBi,UBi
          Hz_inv2(i,k)=1.0_r8/(HzL(i,k)+HzL(i,k-1))
        END DO
      END DO
      DO k=2,N(ng)+1
        DO i=LBi,UBi
          Hz_inv3(i,k)=1.0_r8/(HzL(i,k+1)+HzL(i,k)+HzL(i,k-1))
        END DO
      END DO

      DO k=1,N(ng)
        DO i=LBi,UBi
          qc(i,k)=Bio(i,k)
        END DO
      END DO

!
!  Reconstruct vertical profile of suspended sediment "qc" in terms
!  of a set of parabolic segments within each grid box. Then, compute
!  semi-Lagrangian flux due to sinking.
!
      DO k=2,N(ng)
        DO i=LBi,UBi
          FC(i,k)=(qc(i,k-1)-qc(i,k))*Hz_inv2(i,k)
        END DO
      END DO

      DO k=N(ng),2,-1
        DO i=LBi,UBi
          dltR=HzL(i,k)*FC(i,k)
          dltL=HzL(i,k)*FC(i,k+1)
          cff=HzL(i,k+1)+2.0_r8*HzL(i,k)+HzL(i,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 (qR,qL) 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 qL(k+1)-qR(k) may still have different sign than
!        qc(k+1)-qc(k).  This possibility is excluded, after qL and qR
!        are reconciled using WENO procedure.
!
          cff=(dltR-dltL)*Hz_inv3(i,k)
          dltR=dltR-cff*HzL(i,k-1)
          dltL=dltL+cff*HzL(i,k+1)

          qR(i,k)=qc(i,k)+dltR
          qL(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=N(ng),2,-1
        DO i=LBi,UBi
          dltL=MAX(cff,WL(i,k  ))
          dltR=MAX(cff,WR(i,k-1))
          qR(i,k)=(dltR*qR(i,k)+dltL*qL(i,k-1))/(dltR+dltL)
          qL(i,k-1)=qR(i,k)
        END DO
      END DO

      DO i=LBi,UBi
        FC(i,N(ng))=0.0_r8              ! no-flux boundary condition
#  if defined LINEAR_CONTINUATION
        qL(i,N(ng))=qR(i,N(ng)-1)
        qR(i,N(ng))=2.0_r8*qc(i,N(ng))-qL(i,N(ng))
#  elif defined NEUMANN
        qL(i,N(ng))=qR(i,N(ng)-1)
        qR(i,N(ng))=1.5_r8*qc(i,N(ng))-0.5_r8*qL(i,N(ng))
#  else
        qR(i,N(ng))=qc(i,N(ng))         ! default strictly monotonic
        qL(i,N(ng))=qc(i,N(ng))         ! conditions
        qR(i,N(ng)-1)=qc(i,N(ng))
#  endif
#  if defined LINEAR_CONTINUATION
        qR(i,1)=qL(i,2)
        qL(i,1)=2.0_r8*qc(i,1)-qR(i,1)
#  elif defined NEUMANN
        qR(i,1)=qL(i,2)
        qL(i,1)=1.5_r8*qc(i,1)-0.5_r8*qR(i,1)
#  else
        qL(i,2)=qc(i,1)                 ! bottom grid boxes are
        qR(i,1)=qc(i,1)                 ! re-assumed to be
        qL(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=N(ng),1,-1
        DO i=LBi,UBi
          dltR=qR(i,k)-qc(i,k)
          dltL=qc(i,k)-qL(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
          qR(i,k)=qc(i,k)+dltR
          qL(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.
!
      DO k=N(ng),1,-1
        DO i=LBi,UBi
          cff=dtdays*ABS(wBiod(i,k))
          FC(i,k+1)=0.0_r8
          WL(i,k)=z_wL(i,k+1)+cff
          WR(i,k)=HzL(i,k)*qc(i,k)
          ksource(i,k)=k
        END DO
      END DO
      DO k=N(ng),1,-1
        DO ks=N(ng),k,1
          DO i=LBi,UBi
            IF (WL(i,k).gt.z_wL(i,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=N(ng),1,-1
        DO i=LBi,UBi
          ks=ksource(i,k)

          cu=MIN(1.0_r8,(WL(i,k)-z_wL(i,ks+1))*Hz_inv(i,ks))
          FC(i,k+1)=FC(i,k+1)+                                          &
     &                  HzL(i,ks)*cu*                                   &
     &                  (qL(i,ks)+                                      &
     &                   cu*(0.5_r8*(qR(i,ks)-qL(i,ks))-                &
     &                       (1.5_r8-cu)*                               &
     &                       (qR(i,ks)+qL(i,ks)-2.0_r8*qc(i,ks))))
        END DO
      END DO


      DO k=N(ng),1,-1
        DO i=LBi,UBi
!           The Bio variables are now updated in the main biology subroutine
!             Bio(i,k)=qc(i,k)+(FC(i,k)-FC(i,k-1))*Hz_inv(i,k)
          riseIN(i,k)=FC(i,k)
          riseOUT(i,k)=FC(i,k+1)
          riseOUT(i,N)=0.0_r8
        END DO
      END DO

      RETURN
      END SUBROUTINE BIORISE_1


!=====================================================================
! BIOSINK_2      -original goanpz sink code (C.Lewis+ E. Dobbins)
!=====================================================================

!=====================================================================
!     BIORISE
!=====================================================================

!=======================================================================
!    DetSINK
!----------------
!====================================================================
        Function ComputeDensity(Temp1,Sal1)
!----------------------------------------------------------------
! Computes the water column density from salinity and temperature
! Returns sigma-t
!----------------------------------------------------------------
        USE mod_kinds

        Real(r8) ComputeDensity
        Real(r8) Temp1, Sal1
        Real(r8) Sig
        Sig = 999.842594 + 0.06793952 * Temp1
        Sig = Sig - 0.00909529 * Temp1 ** 2 +                      &
     &          0.0001001685 * Temp1 ** 3
        Sig = Sig - 0.000001120083 * Temp1 ** 4 +                  &
     &          0.000000006536332 * Temp1 ** 5
        Sig = Sig + 0.824493 * Sal1 - 0.0040899 * Temp1 * Sal1
        Sig = Sig + 0.000076438 * Temp1 ** 2 * Sal1 -              &
     &          0.00000082467 * Temp1 ** 3 * Sal1
        Sig = Sig + 0.0000000053875 * Temp1 ** 4 * Sal1 -          &
     &          0.00572466 * Sal1 ** (3 / 2)
        Sig = Sig + 0.00010227 * Temp1 * Sal1 ** (3 / 2) -         &
     &          0.0000016546 * Temp1 ** 2 * Sal1 ** (3 / 2)
        Sig = Sig + 0.00048314 * Sal1 ** 2
         ComputeDensity = Sig - 1000
        End Function ComputeDensity
!===============================================================
        Function GetPhytoResp2(Temp1, Tref, KbmPh)
!------------------------------------------------------
! Computes the temperature correction for phytoplankton
! respiration according to Arhonditsis 2005.
!------------------------------------------------------
        USE mod_kinds

        Real(r8) GetPhytoResp2
        Real(r8) Temp1      !Temperature, passed
        Real(r8) Tref       !Reference temperature
        Real(r8) KbmPh      !Half saturation, temperature

        Real(r8) Resp       !Returned variable

        Resp = exp(KbmPh * (Temp1 - Tref))
        GetPhytoResp2 = Resp
        Return
        End Function GetPhytoResp2
!=====================================================================
      FUNCTION GetLightLimIronSml(alphaPh, PAR1, Pmax1,          &
     &     CrChlRatio1,IronLim1)
!---------------------------------------------------------------
! Uses a normal hyperbolic tangent function for light limitation
!---------------------------------------------------------------
      USE mod_kinds
      USE mod_param
!
      implicit none

        Real(r8) :: alphaPh,Pmax1,IronLim1,PAR1,CrChlRatio1
        Real(r8) GetLightLimIronSml
        Real(r8) LightLim,OffSet
        OffSet = 0.0
          LightLim = TANH( alphaPh * MAX((PAR1 - OffSet),0.0_r8)   &
     &             / Pmax1 / CrChlRatio1 / IronLim1)
       GetLightLimIronSml = LightLim
      END FUNCTION GetLightLimIronSml
!=====================================================================
      FUNCTION GetLightLimSml(alphaPh, PAR1, Pmax1, CrChlRatio1)
!---------------------------------------------------------------
! Uses a normal hyperbolic tangent function for light limitation
!---------------------------------------------------------------
      USE mod_kinds
      USE mod_param
!
      implicit none

        Real(r8) :: alphaPh,Pmax1,IronLim1,PAR1,CrChlRatio1
        Real(r8) GetLightLimSml
        Real(r8) LightLim,OffSet
        OffSet = 0.0
          LightLim = TANH( alphaPh * MAX((PAR1 - OffSet),0.0_r8)   &
     &             / Pmax1 / CrChlRatio1 )
       GetLightLimSml = LightLim
      END FUNCTION GetLightLimSml
!=====================================================================
      FUNCTION GetLightLimIron(alphaPh, PAR1, Pmax1, CrChlRatio1,  &
     &     IronLim1, ParMax)
!------------------------------------------------------------------
! Light lim with varying alpha. Works with iron limitation. Alph is
! a function of the surface light intensity.
!------------------------------------------------------------------
      USE mod_kinds
      USE mod_param
!
      implicit none

       Real(r8) :: alphaPh,Pmax1,IronLim1,PAR1,CrChlRatio1
       Real(r8) GetLightLimIron
       Real(r8) Alpha,LightLim,OffSet,ParMax

       !Alpha = 1.e-8*EXP(0.48*ParMax) + 0.5
       !if (Alpha .gt. 10) Alpha = 10
       !--------------------------
       !Use a simple step function
       !--------------------------
       if (ParMax.lt.48_r8) then
          Alpha = 1._r8
       else
          Alpha = 4._r8
       end if
       LightLim = TANH( Alpha * Par1/Pmax1/CrChlRatio1/IronLim1)
       GetLightLimIron = LightLim

      END FUNCTION GetLightLimIron
!=======================================================================
      FUNCTION GetLightLim(alphaPh, PAR1, Pmax1, CrChlRatio1, ParMax)
!-----------------------------------------------------------------
! Generates a light lim with varying alphaPh without iron
!-----------------------------------------------------------------
      USE mod_kinds
      USE mod_param
!
      implicit none

       Real(r8) :: alphaPh,Pmax1,PAR1,CrChlRatio1
       Real(r8) GetLightLim
       Real(r8) Alpha,LightLim,OffSet,ParMax

!       Alpha = 1.e-8*EXP(0.48*ParMax) + 0.5
!       if (Alpha .gt. 10) Alpha = 10
       !--------------------------
       !Use a simple step function
       !--------------------------
       if (ParMax.lt.48_r8) then
          Alpha = 1._r8
       else
          Alpha = 4._r8
       end if
       LightLim = TANH(alphaPh * PAR1/Pmax1/CrChlRatio1)

       GetLightLim = LightLim
      END FUNCTION GetLightLim
!==============================================================
        Function GetCopepodResp(Temp1,respVal,ktbm,Tref)
!--------------------------------------------------------------
! Computes copepod respiration according to Arhonditsis (2005).
!--------------------------------------------------------------
        USE mod_kinds
        USE mod_param

        Real(r8) GetCopepodResp
        real(r8) :: respVal
        Real(r8) Temp1       !Passed variable
!        Real(r8) :: bm  = 0.04   !Basal metabolic rate day**-1
        Real(r8) :: ktbm  != 0.05 !Temperature response degrees C**-1
        Real(r8) :: Tref  != 20   !Reference temperature degrees C
        Real(r8) Resp        !Returned variable

        Resp = respVal * exp(ktbm * (Temp1 - Tref))
        GetCopepodResp = Resp
        Return
        End Function GetCopepodResp
!==============================================================
        Function GetJelResp(Temp1,bmJ,ktbmJ,TrefJ)

        USE mod_param
        implicit none
!
        real(r8) :: GetJelResp
        real(r8) :: Temp1          !Passed variable
        real(r8) :: bmJ   != 0.04   !Basal metabolic rate day**-1
        real(r8) :: ktbmJ != 0.05 !Temperature response degrees C**-1
        real(r8) :: TrefJ != 20   !Reference temperature degrees C
        real(r8) :: Resp           !Returned variable

        Resp = bmJ * exp(ktbmJ * (Temp1 - TrefJ))
        GetJelResp = Resp
        Return
        End Function GetJelResp
!=================================================================
        Function GetBasalMetabolism(respPh,kfePh,Iron1)
!---------------------------------------------------------
! Computes an iron correction for the basal metabolism for
! the phytoplankton respiration calculation
!---------------------------------------------------------
        USE mod_kinds

        Real(r8) GetBasalMetabolism
        Real(r8) Iron1     !Iron concentration
        Real(r8) kfePh     !Half saturation for iron
        Real(r8) respPh    !Phytoplankton uncorrected basal metabolism
        Real(r8) BaseMet   !Phytoplankton basal metabolism

        BaseMet = Iron1/(kfePh + Iron1)
        BaseMet = BaseMet * ((kfePh +2)/2)
        GetBasalMetabolism = BaseMet * respPh

        Return
        End Function GetBasalMetabolism
!=========================================================
        Function GetNitrif(Temp1,Dep1,NH4R)
!-------------------------------------------------------
! Computes the nitrification with respect to temperature
! according to Arhonditsis (2005).  Generates depth
! correction according to Denman (2003).
!-------------------------------------------------------
        USE mod_kinds
        Real(r8) GetNitrif
                               !--------------------------------
        Real(r8) Temp1, Dep1, NH4R !Passed variables
        Real(r8) NH4conv  /14/     !mg N/mol N
        Real(r8) KNH4Nit  /0.08/   !Half Sat Con mg N/m3/day
        Real(r8) KTNitr   /0.002/  !Temperature responce dec C^2
        Real(r8) ToptNtr  /28/     !Optimum nitrification temp
        Real(r8) Zox      /20/     !50% nitrification depth
        Real(r8) Nexp     /6/      !Exponent to adjust profile shape
        Real(r8) NitrMax  /0.011/  !Maximum nitrification (mM/m3/d
                               !--------------------------------

        Real(r8) Nitr, DepCor

        NH4R = NH4R * NH4conv
        Nitr = NH4R/(KNH4Nit + NH4R)
        Nitr = Nitr * exp(-KTNitr*(Temp1 - ToptNtr)**2)
        DepCor = (Dep1**Nexp)/( (Zox**Nexp) + Dep1**Nexp)
        Nitr = (Nitr * DepCor) * NitrMax
        GetNitrif = Nitr
        Return
        End Function GetNitrif
!========================================================================
        Function GetNitrif2(Temp1,PAR1,NH4R)
        !---------------------------------------------------------
        !Computes nitrificaton from Kawamiya with light correction
        !from Fennel; Kawamiya (2000), Fennel (2006)
        !---------------------------------------------------------
        USE mod_kinds
        USE mod_param
         Real(r8) GetNitrif2
                               !--------------------------------
        Real(r8) :: Temp1, NH4R       !Passed variables
        Real(r8) :: I0 = 0.0095     !Threshold,light inhibition, W m-2
        Real(r8) :: KI = 4.0        !Half Saturation light intensity, W m-2
        Real(r8) :: KN0 = 0.03       !Nitrification at 0 deg C, day-1
        Real(r8) :: KNT = 0.0693     !Temperature coefficient
        Real(r8) :: ParW              !Par in watts
        Real(r8) :: NitrMax           !Maximum nitrification
        Real(r8) :: Nitr              !Nitrification
                               !---------------------------------
        Real(r8) :: cff1, PAR1

        !-----------------------------------
        !Temperature dependent nitrification
        !-----------------------------------
        KI = 1.5
        KN0 = 0.15
        !KNT = 0.07
        NitrMax = (KN0*Exp(KNT*Temp1))*NH4R
        !-----------------------------------
        !Convert PAR in E m-2 d-1 to W day-1
        !-----------------------------------
        ParW = PAR1/0.394848_r8
        !---------------------------------
        !Light correction of nitrification
        !---------------------------------
        cff1 = (ParW-I0)/(KI+ParW-I0)
        Nitr = NitrMax*(1-MAX(0.0_r8,cff1))
        GetNitrif2 = Nitr
        Return
        End Function GetNitrif2
!-----------------------------------------------------------

        Function GetNitrifLight(Par1,tI0,KI)
         USE mod_kinds

        Real(r8) GetNitrifLight
        Real(r8) :: Par1               !--------------------------------
        Real(r8) :: tI0                 !Threshold,light inhibition, W m-2
        Real(r8) :: KI                 !Half Saturation light intensity, W m-2
        Real(r8) :: ParW               !Par in watts
        Real(r8) :: cff1

            ParW = Par1/0.394848_r8 !convert PAR back to watts
            cff1 = (ParW-tI0)/(KI+ParW-tI0)

            GetNitrifLight=(1-MAX(0.0_r8,cff1))

        Return
        End Function GetNitrifLight
!-----------------------------------------------
         Function GetNitrifMaxK(Temp1,KnT,Nitr0)
!-------------------------------------------------------
! Computes the nitrification with respect to temperature
! according to Arhonditsis (2005).  Generates depth
! correction according to Denman (2003).
!-------------------------------------------------------
        USE mod_kinds
        Real(r8) GetNitrifMaxK

        Real(r8) :: Temp1         !--------------------------------
        Real(r8) :: KnT           !Passed variables
        Real(r8) :: Nitr0         !--------------------------------


        GetNitrifMaxK=Nitr0*exp(KnT*Temp1)

        Return
        End Function GetNitrifMaxK

!-----------------------------------------------------------

        Function GetNitrifMaxA(Nitr0, ktntr,Temp1,ToptNtr)
         USE mod_kinds

        Real(r8) GetNitrifMaxA

        Real(r8) :: Temp1
        Real(r8) :: Ktntr           !--------------------------------
        Real(r8) :: Nitr0           !Passed variables
        Real(r8) :: ToptNtr         !--------------------------------


        GetNitrifMaxA=Nitr0*exp(-ktntr*(Temp1 - ToptNtr)**2)


        Return
        End Function GetNitrifMaxA

!-----------------------------------------------------------
!=====================================================================
      FUNCTION GetLightLimIron2(alphaPh, PAR1, Pmax1, CrChlRatio1, &
     &     IronLim1)
!---------------------------------------------------------------
! Uses a normal hyperbolic tangent function for light limitation
!---------------------------------------------------------------
      USE mod_kinds
      USE mod_param
!
      implicit none

        Real(r8) :: alphaPh,Pmax1,IronLim1,PAR1,CrChlRatio1
        Real(r8) GetLightLimIron2
        Real(r8) LightLim,OffSet
        OffSet = 0.0_r8
          LightLim = TANH( alphaPh * MAX((PAR1 - OffSet),0.0_r8)   &
     &             / Pmax1 / CrChlRatio1 / IronLim1)
       GetLightLimIron2 = LightLim
      END FUNCTION GetLightLimIron2
!=====================================================================
      FUNCTION GetLightLim2(alphaPh, PAR1, Pmax1, CrChlRatio1)
!---------------------------------------------------------------
! Uses a normal hyperbolic tangent function for light limitation
!---------------------------------------------------------------
      USE mod_kinds
      USE mod_param
!
      implicit none

        Real(r8) :: alphaPh,Pmax1,IronLim1,PAR1,CrChlRatio1
        Real(r8) GetLightLim2
        Real(r8) LightLim,OffSet
        OffSet = 0.0_r8
          LightLim = TANH( alphaPh * MAX((PAR1 - OffSet),0.0_r8)   &
     &             / Pmax1 / CrChlRatio1 )
       GetLightLim2 = LightLim
      END FUNCTION GetLightLim2
