MODULE biology_mod
!
!git $Id$
!svn $Id: ecosim.h 1151 2023-02-09 03:08:53Z arango $
!================================================== Hernan G. Arango ===
! Copyright (c) 2002-2023 The ROMS/TOMS Group !
! Licensed under a MIT/X style license !
! See License_ROMS.md !
!=================================================== W. Paul Bissett ===
! Copyright (c) 1997 W. Paul Bissett, FERI !
!=======================================================================
! !
! The EcoSim code has been developed for research purposes only. It !
! consists of unpublished, proprietary formulations protected under !
! U.S. copyright law. It is freely available on request from the !
! Florida Environmental Research Institute (FERI). Commercial usage !
! of these formulations is forbidden without express written !
! permission from FERI. All rights reserved. !
! !
!***********************************************************************
! !
! This routine computes the EcoSim sources and sinks and adds them !
! to the global biological fields. !
! !
! Reference: !
! !
! Bissett, W.P., J.J. Walsh, D.A. Dieterle, K.L. Carder, 1999: !
! Carbon cycling in the upper waters of the Sargasso Sea: I. !
! Numerical simulation of differential carbon and nitrogen !
! fluxes, Deep-Sea Res., 46, 205-269. !
! !
! Bissett, W.P., K.L. Carder, J.J. Walsh, D.A. Dieterle, 1999: !
! Carbon cycling in the upper waters of the Sargasso Sea: II. !
! Numerical simulation of apparent and inherent optical !
! properties, Deep-Sea Res., 46, 271-317 !
! !
! NOTES to EcoSim: !
! !
! * This version uses a descending index for depth that is different !
! than the original coding. !
! !
! * This version of the code has been modified by Bronwyn Cahill and !
! includes a semi-Lagrangian vertical sinking flux algorithm for !
! fecal material and bio_sediment subroutine which remineralizes !
! particulate nitrogen and returns it into dissolved nitrate pool. !
! !
!=======================================================================
!
implicit none
!
PRIVATE
PUBLIC :: biology
!
CONTAINS
!
!***********************************************************************
SUBROUTINE biology (ng, tile)
!***********************************************************************
!
USE mod_param
#ifdef DIAGNOSTICS_BIO
USE mod_diags
#endif
USE mod_forces
USE mod_grid
USE mod_ncparam
USE mod_ocean
USE mod_stepping
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
!
! Local variable declarations.
!
character (len=*), parameter :: MyFile = &
& __FILE__
!
#include "tile.h"
!
! Set header file name.
!
#ifdef DISTRIBUTE
IF (Lbiofile(iNLM)) THEN
#else
IF (Lbiofile(iNLM).and.(tile.eq.0)) THEN
#endif
Lbiofile(iNLM)=.FALSE.
BIONAME(iNLM)=MyFile
END IF
!
#ifdef PROFILE
CALL wclock_on (ng, iNLM, 15, __LINE__, MyFile)
#endif
CALL ecosim_tile (ng, tile, &
& LBi, UBi, LBj, UBj, N(ng), NT(ng), &
& IminS, ImaxS, JminS, JmaxS, &
& nstp(ng), nnew(ng), &
#ifdef MASKING
& GRID(ng) % rmask, &
# if defined WET_DRY && defined DIAGNOSTICS_BIO
& GRID(ng) % rmask_full, &
# endif
#endif
& GRID(ng) % Hz, &
& GRID(ng) % z_r, &
& GRID(ng) % z_w, &
& FORCES(ng) % SpecIr, &
& FORCES(ng) % avcos, &
#ifdef DIAGNOSTICS_BIO
& DIAGS(ng) % DiaBio3d, &
& DIAGS(ng) % DiaBio4d, &
#endif
& OCEAN(ng) % t)
#ifdef PROFILE
CALL wclock_off (ng, iNLM, 15, __LINE__, MyFile)
#endif
!
RETURN
END SUBROUTINE biology
!
!***********************************************************************
SUBROUTINE ecosim_tile (ng, tile, &
& LBi, UBi, LBj, UBj, UBk, UBt, &
& IminS, ImaxS, JminS, JmaxS, &
& nstp, nnew, &
#ifdef MASKING
& rmask, &
# if defined WET_DRY && defined DIAGNOSTICS_BIO
& rmask_full, &
# endif
#endif
& Hz, z_r, z_w, &
& SpecIr, avcos, &
#ifdef DIAGNOSTICS_BIO
& DiaBio3d, &
& DiaBio4d, &
#endif
& t)
!***********************************************************************
!
USE mod_param
USE mod_biology
USE mod_eclight
USE mod_scalars
USE mod_iounits
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: nstp, nnew
#ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:,LBj:)
# if defined WET_DRY && defined DIAGNOSTICS_BIO
real(r8), intent(in) :: rmask_full(LBi:,LBj:)
# endif
# 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) :: SpecIr(LBi:,LBj:,:)
real(r8), intent(in) :: avcos(LBi:,LBj:,:)
# ifdef DIAGNOSTICS_BIO
real(r8), intent(inout) :: DiaBio3d(LBi:,LBj:,:,:)
real(r8), intent(inout) :: DiaBio4d(LBi:,LBj:,:,:,:)
# endif
real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
#else
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# if defined WET_DRY && defined DIAGNOSTICS_BIO
real(r8), intent(in) :: rmask_full(LBi:UBi,LBj:UBj)
# endif
# 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) :: SpecIr(LBi:UBi,LBj:UBj,NBands)
real(r8), intent(in) :: avcos(LBi:UBi,LBj:UBj,NBands)
# ifdef DIAGNOSTICS_BIO
real(r8), intent(inout) :: DiaBio3d(LBi:UBi,LBj:UBj, &
& NDbands,NDbio3d)
real(r8), intent(inout) :: DiaBio4d(LBi:UBi,LBj:UBj,N(ng), &
& NDbands,NDbio4d)
# endif
real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
#endif
!
! Local variable declarations.
!
integer, parameter :: Msink = 30
integer :: i, j, k, ks
integer :: Iter, Tindex, ic, isink, ibio, id, itrc, ivar
integer :: ibac, iband, idom, ifec, iphy, ipig
integer :: Nsink
integer, dimension(Msink) :: idsink
integer, dimension(IminS:ImaxS,N(ng)) :: ksource
real(r8), parameter :: MinVal = 0.0_r8
real(r8) :: FV1, FV2, FV3, FV4, FV5, FV6, FV7, dtbio
real(r8) :: DOC_lab, Ed_tot, Nup_max, aph442, aPHYN_wa
real(r8) :: avgcos_min, par_b, par_s, photo_DIC, photo_DOC
real(r8) :: photo_decay, slope_AC, tChl, theta_m, total_photo
real(r8) :: WLE, factint
real(r8) :: Het_BAC
real(r8) :: N_quota, RelDOC1, RelDON1, RelDOP1, RelFe
real(r8) :: cff, cff1, cff2, cffL, cffR, cu, dltL, dltR
#ifdef DIAGNOSTICS_BIO
real(r8) :: fiter
#endif
real(r8), dimension(Msink) :: Wbio
real(r8), dimension(4) :: Bac_G
real(r8), dimension(NBands) :: dATT_sum
real(r8), dimension(N(ng),NBands) :: avgcos, dATT
real(r8), dimension(N(ng),NBands) :: specir_d
real(r8), dimension(N(ng),NBands) :: tot_ab, tot_b, tot_s
#ifdef BIO_OPTIC
real(r8), dimension(0:N(ng),NBands) :: specir_w
#endif
real(r8), dimension(N(ng),Nphy) :: C2CHL, C2CHL_w
real(r8), dimension(N(ng),Nphy) :: Gt_fl, Gt_ll, Gt_nl
real(r8), dimension(N(ng),Nphy) :: Gt_sl, Gt_pl
real(r8), dimension(N(ng),Nphy) :: alfa
real(r8), dimension(N(ng),Nphy) :: pac_eff
real(r8), dimension(N(ng),Nphy,Npig) :: Pigs_w
integer, dimension(IminS:ImaxS) :: Keuphotic
real(r8), dimension(IminS:ImaxS,N(ng)) :: E0_nz
real(r8), dimension(IminS:ImaxS,N(ng)) :: Ed_nz
real(r8), dimension(IminS:ImaxS,N(ng)) :: DOC_frac
real(r8), dimension(IminS:ImaxS,N(ng)) :: NitrBAC
real(r8), dimension(IminS:ImaxS,N(ng)) :: NH4toNO3
real(r8), dimension(IminS:ImaxS,N(ng)) :: NtoNBAC
real(r8), dimension(IminS:ImaxS,N(ng)) :: NtoPBAC
real(r8), dimension(IminS:ImaxS,N(ng)) :: NtoFeBAC
real(r8), dimension(IminS:ImaxS,N(ng)) :: totDOC_d
real(r8), dimension(IminS:ImaxS,N(ng)) :: totDON_d
real(r8), dimension(IminS:ImaxS,N(ng)) :: totDOP_d
real(r8), dimension(IminS:ImaxS,N(ng)) :: totFe_d
real(r8), dimension(IminS:ImaxS,N(ng)) :: totNH4_d
real(r8), dimension(IminS:ImaxS,N(ng)) :: totNO3_d
real(r8), dimension(IminS:ImaxS,N(ng)) :: totPO4_d
real(r8), dimension(IminS:ImaxS,N(ng)) :: totSiO_d
real(r8), dimension(IminS:ImaxS,N(ng),Nbac) :: GtBAC
real(r8), dimension(IminS:ImaxS,N(ng),Nbac) :: NupDOC_ba
real(r8), dimension(IminS:ImaxS,N(ng),Nbac) :: NupDON_ba
real(r8), dimension(IminS:ImaxS,N(ng),Nbac) :: NupDOP_ba
real(r8), dimension(IminS:ImaxS,N(ng),Nbac) :: NupFe_ba
real(r8), dimension(IminS:ImaxS,N(ng),Nbac) :: NupNH4_ba
real(r8), dimension(IminS:ImaxS,N(ng),Nbac) :: NupPO4_ba
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: C2fALG
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: C2nALG
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: C2pALG
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: C2sALG
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: GtALG
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: GtALG_r
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: NupDOP
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: NupDON
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: NupFe
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: NupNH4
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: NupNO3
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: NupPO4
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: NupSiO
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: graz_act
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_f
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_n
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_p
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_s
real(r8), dimension(IminS:ImaxS,N(ng),Nphy) :: refuge
real(r8), dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_C
real(r8), dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_F
real(r8), dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_N
real(r8), dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_P
real(r8), dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_S
real(r8), dimension(IminS:ImaxS,N(ng),NBands) :: specir_scal
real(r8), dimension(IminS:ImaxS,N(ng),Nphy,NBands) :: aPHYN_al
real(r8), dimension(IminS:ImaxS,N(ng),Nphy,NBands) :: aPHYN_at
real(r8), dimension(IminS:ImaxS,N(ng),NBands) :: aDET
real(r8), dimension(IminS:ImaxS,N(ng),NBands) :: aCDC
real(r8), dimension(IminS:ImaxS,N(ng),NBands) :: b_phy
real(r8), dimension(IminS:ImaxS,N(ng),NBands) :: s_phy
real(r8), dimension(IminS:ImaxS,N(ng),NBands) :: b_tot
real(r8), dimension(IminS:ImaxS,N(ng),NBands) :: s_tot
real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_new
real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
#include "set_bounds.h"
#ifdef DIAGNOSTICS_BIO
!
!-----------------------------------------------------------------------
! If appropriate, initialize time-averaged diagnostic arrays.
!-----------------------------------------------------------------------
!
IF (((iic(ng).gt.ntsDIA(ng)).and. &
& (MOD(iic(ng),nDIA(ng)).eq.1)).or. &
& ((iic(ng).ge.ntsDIA(ng)).and.(nDIA(ng).eq.1)).or. &
& ((nrrec(ng).gt.0).and.(iic(ng).eq.ntstart(ng)))) THEN
DO ivar=1,NDbio3d
DO k=1,NDbands
DO j=Jstr,Jend
DO i=Istr,Iend
DiaBio3d(i,j,k,ivar)=0.0_r8
END DO
END DO
END DO
END DO
DO ivar=1,NDbio4d
DO iband=1,NDbands
DO k=1,N(ng)
DO j=Jstr,Jend
DO i=Istr,Iend
DiaBio4d(i,j,k,iband,ivar)=0.0_r8
END DO
END DO
END DO
END DO
END DO
END IF
#endif
!
!=======================================================================
! Add EcoSim Source/Sink terms.
!=======================================================================
!
! Set internal time-stepping.
!
dtbio=dt(ng)/REAL(BioIter(ng),r8)
#ifdef DIAGNOSTICS_BIO
!
! A factor to account for the number of iterations in accumulating
! diagnostic rate variables.
!
fiter=1.0_r8/REAL(BioIter(ng),r8)
#endif
!
! Set vertical sinking identification and associated sinking velocity
! arrays.
!
ic=1
DO ifec=1,Nfec
idsink(ic)=iFecN(ifec)
Wbio(ic)=WF(ifec,ng)
ic=ic+1
idsink(ic)=iFecC(ifec)
Wbio(ic)=WF(ifec,ng)
ic=ic+1
idsink(ic)=iFecP(ifec)
Wbio(ic)=WF(ifec,ng)
ic=ic+1
idsink(ic)=iFecS(ifec)
Wbio(ic)=WF(ifec,ng)
ic=ic+1
idsink(ic)=iFecF(ifec)
Wbio(ic)=WF(ifec,ng)
ic=ic+1
END DO
DO iphy=1,Nphy
idsink(ic)=iPhyN(iphy)
Wbio(ic)=WS(iphy,ng)
ic=ic+1
idsink(ic)=iPhyC(iphy)
Wbio(ic)=WS(iphy,ng)
ic=ic+1
idsink(ic)=iPhyP(iphy)
Wbio(ic)=WS(iphy,ng)
IF (iPhyS(iphy).ne.0) THEN
ic=ic+1
idsink(ic)=iPhyS(iphy)
Wbio(ic)=WS(iphy,ng)
END IF
ic=ic+1
idsink(ic)=iPhyF(iphy)
Wbio(ic)=WS(iphy,ng)
ic=ic+1
END DO
Nsink=ic-1
!
!-----------------------------------------------------------------------
! Compute inverse thickness to avoid repeated divisions.
!-----------------------------------------------------------------------
!
J_LOOP : DO j=Jstr,Jend
DO k=1,N(ng)
DO i=Istr,Iend
Hz_inv(i,k)=1.0_r8/Hz(i,j,k)
END DO
END DO
DO k=1,N(ng)-1
DO i=Istr,Iend
Hz_inv2(i,k)=1.0_r8/(Hz(i,j,k)+Hz(i,j,k+1))
END DO
END DO
DO k=2,N(ng)-1
DO i=Istr,Iend
Hz_inv3(i,k)=1.0_r8/(Hz(i,j,k-1)+Hz(i,j,k)+Hz(i,j,k+1))
END DO
END DO
!
!-----------------------------------------------------------------------
! Extract biological variables from tracer arrays, place them into
! scratch arrays, and restrict their values to be positive definite.
!-----------------------------------------------------------------------
!
! At input, all tracers (index nnew) from predictor step have
! transport units (m Tunits) since we do not have yet the new
! values for zeta and Hz. These are known after the 2D barotropic
! time-stepping.
!
DO ibio=1,NBT
itrc=idbio(ibio)
DO k=1,N(ng)
DO i=Istr,Iend
Bio(i,k,itrc)=MAX(MinVal,t(i,j,k,nstp,itrc))
Bio_old(i,k,itrc)=Bio(i,k,itrc)
!!
!! HGA - The new tendency terms were not initialized. This gives
!! unexpected behavior on different computers since a variable
!! was used before it was assigned. This may explain earlier
!! problems with the algorithm. Perhaps, this time-stepping can
!! be modified latter to avoid unnecessary storage between
!! Bio, Bio_old, and Bio_new.
!!
Bio_new(i,k,itrc)=0.0_r8
END DO
END DO
END DO
!
! Extract potential temperature and salinity.
!
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)
END DO
END DO
!
!-----------------------------------------------------------------------
! Compute temperature and salinity dependent variables.
!-----------------------------------------------------------------------
!
! Refuge depth calculation.
!
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
refuge(i,k,iphy)=MinRefuge(iphy,ng)
END DO
END DO
END DO
!
! Initialize fecal regeneration arrays (N, P, and Fe from Moore et al.,
! DSRII 2001; silica is given by values from Bidle and Azam, Nature,
! 1999).
!
IF (Regen_flag(ng)) THEN
DO ifec=1,Nfec
DO k=1,N(ng)
DO i=Istr,Iend
FV1=EXP(RegTfac(ifec,ng)*(Bio(i,k,itemp)- &
& RegTbase(ifec,ng)))
Regen_C(i,k,ifec)=RegCR(ifec,ng)*FV1
Regen_N(i,k,ifec)=RegNR(ifec,ng)*FV1
Regen_P(i,k,ifec)=RegPR(ifec,ng)*FV1
Regen_F(i,k,ifec)=RegFR(ifec,ng)*FV1
Regen_S(i,k,ifec)=RegSR(ifec,ng)*FV1
END DO
END DO
END DO
END IF
!
! Calculate temperature dependent growth rate.
!
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
GtALG(i,k,iphy)=GtALG_max(iphy,ng)* &
& EXP(PhyTfac(iphy,ng)* &
& (Bio(i,k,itemp)-PhyTbase(iphy,ng)))
!
! Calculate mu_bar for droop equation.
!
FV1=maxC2nALG(iphy,ng)*(1.0_r8+GtALG(i,k,iphy))
mu_bar_n(i,k,iphy)=GtALG(i,k,iphy)* &
& FV1/(FV1-minC2nALG(iphy,ng))
IF (HsSiO(iphy,ng).lt.LARGER) THEN
FV1=maxC2SiALG(iphy,ng)*(1.0_r8+GtALG(i,k,iphy))
mu_bar_s(i,k,iphy)=GtALG(i,k,iphy)* &
& FV1/(FV1-minC2SiALG(iphy,ng))
ELSE
mu_bar_s(i,k,iphy)=LARGER
END IF
IF (HsPO4(iphy,ng).lt.LARGER) THEN
FV1=maxC2pALG(iphy,ng)*(1.0_r8+GtALG(i,k,iphy))
mu_bar_p(i,k,iphy)=GtALG(i,k,iphy)* &
& FV1/(FV1-minC2pALG(iphy,ng))
ELSE
mu_bar_p(i,k,iphy)=LARGER
END IF
IF (HsFe(iphy,ng).lt.LARGER) THEN
FV1=maxC2FeALG(iphy,ng)*(1.0_r8+GtALG(i,k,iphy))
mu_bar_f(i,k,iphy)=GtALG(i,k,iphy)* &
& FV1/(FV1-minC2FeALG(iphy,ng))
ELSE
mu_bar_f(i,k,iphy)=LARGER
END IF
END DO
END DO
END DO
!
! Bacterial growth rate from Fasham et al., 1990.
!
DO ibac=1,Nbac
DO k=1,N(ng)
DO i=Istr,Iend
GtBAC(i,k,ibac)=GtBAC_max(ibac,ng)* &
& EXP(BacTfac(ibac,ng)* &
& (Bio(i,k,itemp)-BacTbase(ibac,ng)))
END DO
END DO
END DO
!
! Grazing rate calculation.
! NOTE: ES1 included separation calculations for grazing beneath the
! zone of refuge (250 m). This has been removed and may
! result in differences in deeper waters.
!! Revisions, WPB 10/20/02. New grazing formulation that is better
!! representation of basal loss rates and biomass accumulations.
!
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
FV1=MAX(1.0_r8,(Bio(i,k,iPhyC(iphy))/refuge(i,k,iphy)))
graz_act(i,k,iphy)=HsGRZ(iphy,ng)*LOG(FV1)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
! Iterate biology source and sink terms.
!-----------------------------------------------------------------------
!
ITER_LOOP : DO Iter=1,BioIter(ng)
DO k=1,N(ng)
DO i=Istr,Iend
totNH4_d(i,k)=0.0_r8
totNO3_d(i,k)=0.0_r8
totPO4_d(i,k)=0.0_r8
totSiO_d(i,k)=0.0_r8
totFe_d (i,k)=0.0_r8
totDOC_d(i,k)=0.0_r8
totDON_d(i,k)=0.0_r8
totDOP_d(i,k)=0.0_r8
END DO
END DO
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
NupNH4(i,k,iphy)=0.0_r8
NupNO3(i,k,iphy)=0.0_r8
NupPO4(i,k,iphy)=0.0_r8
NupSiO(i,k,iphy)=0.0_r8
NupFe (i,k,iphy)=0.0_r8
NupDON(i,k,iphy)=0.0_r8
NupDOP(i,k,iphy)=0.0_r8
END DO
END DO
END DO
!
! Compute Ratio Arrays.
! (Calculating only those that are accessed more than once.)
!
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
C2nALG(i,k,iphy)=0.0_r8
IF (Bio(i,k,iPhyN(iphy)).gt.0.0_r8) THEN
C2nALG(i,k,iphy)=Bio(i,k,iPhyC(iphy))/ &
& Bio(i,k,iPhyN(iphy))
END IF
C2pALG(i,k,iphy)=0.0_r8
IF (Bio(i,k,iPhyP(iphy)).gt.0.0_r8) THEN
C2pALG(i,k,iphy)=Bio(i,k,iPhyC(iphy))/ &
& Bio(i,k,iPhyP(iphy))
END IF
C2sALG(i,k,iphy)=0.0_r8
IF (iPhyS(iphy).gt.0) THEN
IF (Bio(i,k,iPhyS(iphy)).gt.0.0_r8) THEN
C2sALG(i,k,iphy)=Bio(i,k,iPhyC(iphy))/ &
& Bio(i,k,iPhyS(iphy))
END IF
END IF
C2fALG(i,k,iphy)=0.0_r8
IF (Bio(i,k,iPhyF(iphy)).gt.0.0_r8) THEN
C2fALG(i,k,iphy)=Bio(i,k,iPhyC(iphy))/ &
& Bio(i,k,iPhyF(iphy))
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
! Daylight Computations.
!-----------------------------------------------------------------------
!
! Initialize.
!
DO i=Istr,Iend
Ed_nz(i,N(ng))=0.0_r8
E0_nz(i,N(ng))=0.0_r8
Keuphotic(i)=N(ng)+1
IF (SpecIr(i,j,21).gt.VSMALL) THEN
DO k=1,N(ng)-1
Ed_nz(i,k)=0.0_r8
E0_nz(i,k)=0.0_r8
END DO
DO iband=1,NBands
dATT_sum(iband)=0.0_r8
DO k=1,N(ng)
avgcos(k,iband)=0.0_r8
dATT(k,iband)=0.0_r8
#ifdef DIAGNOSTICS_BIO
aDET(i,k,iband)=0.0_r8
aCDC(i,k,iband)=0.0_r8
b_phy(i,k,iband)=0.0_r8
s_phy(i,k,iband)=0.0_r8
b_tot(i,k,iband)=0.0_r8
s_tot(i,k,iband)=0.0_r8
#endif
END DO
DO iphy=1,Nphy
DO k=1,N(ng)
aPHYN_at(i,k,iphy,iband)=0.0_r8
aPHYN_al(i,k,iphy,iband)=0.0_r8
END DO
END DO
END DO
!
! Calculate average cosine zenith angle at surface.
! (See equation 14 Morel, 1991 Prog. Ocean.)
!
Ed_tot=0.0_r8
DO iband=1,NBands
Ed_tot=Ed_tot+SpecIr(i,j,iband)*DLAM
avgcos(N(ng),iband)=avcos(i,j,iband)
END DO
!
! Total aph(442). adp(442) is set to 50% of aph(442).
! NOTE: choosing sbands=9 which is band 442 using v8r16
! sbands formulation. If spectral resolution changes, this
! value must change!
!
DO k=N(ng),1,-1
IF (Ed_tot.ge.1.0_r8) THEN
aph442=0.0_r8
tChl=0.0_r8
DO iphy=1,Nphy
IF (Bio(i,k,iPhyC(iphy)).gt.0.0_r8) THEN
tChl=tChl+Bio(i,k,iPigs(iphy,ichl))
pac_eff(k,iphy)=1.0_r8
IF (b_PacEff(iphy,ng).gt.SMALL) THEN
FV2=Bio(i,k,iPigs(iphy,ichl))/ &
& (Bio(i,k,iPhyC(iphy))*12.0_r8)
pac_eff(k,iphy)=MAX(0.5_r8, &
& (MIN(1.0_r8, &
& b_PacEff(iphy,ng)+ &
& mxPacEff(iphy,ng)* &
& (FV2- &
& b_C2Cl(iphy,ng)))))
END IF
iband=9
DO ipig=1,Npig
IF (iPigs(iphy,ipig).gt.0) THEN
aph442=aph442+ &
& Bio(i,k,iPigs(iphy,ipig))* &
& apigs(ipig,iband)*pac_eff(k,iphy)
END IF
END DO
END IF
END DO
!
! Calculate absorption.
! Calculating phytoplankton absorption for attenuation calculation.
! NOTE: 12 factor to convert to ugrams (mg m-3)
!
aph442=0.5_r8*aph442
DO iband=1,NBands
tot_ab=0.0_r8
DO iphy=1,Nphy
DO ipig=1,Npig
IF (iPigs(iphy,ipig).gt.0) THEN
aPHYN_at(i,k,iphy,iband)= &
& aPHYN_at(i,k,iphy,iband)+ &
& Bio(i,k,iPigs(iphy,ipig))* &
& apigs(ipig,iband)* &
& pac_eff(k,iphy)
END IF
END DO
tot_ab(k,iband)=tot_ab(k,iband)+ &
& aPHYN_at(i,k,iphy,iband)
#ifdef DIAGNOSTICS_BIO
DiaBio4d(i,j,k,iband,idaPHY)=DiaBio4d(i,j,k,iband,&
& idaPHY)+ &
& aPHYN_at(i,k,iphy, &
& iband)
#endif
!
! Removing absorption due to PPC for "alfa" calculation.
!
ipig=5
IF (iPigs(iphy,ipig).gt.0) THEN
aPHYN_al(i,k,iphy,iband)= &
& aPHYN_at(i,k,iphy,iband)- &
& Bio(i,k,iPigs(iphy,ipig))* &
& apigs(ipig,iband)* &
& pac_eff(k,iphy)
END IF
END DO
!
! Adding detrital absorption.
!
cff=aph442*EXP(0.011_r8* &
& (442.0_r8- &
& (397.0_r8+REAL(iband,r8)*DLAM)))
tot_ab(k,iband)=tot_ab(k,iband)+cff
#ifdef DIAGNOSTICS_BIO
aDET(i,k,iband)=aDET(i,k,iband)+cff
DiaBio4d(i,j,k,iband,idaDET)=DiaBio4d(i,j,k,iband, &
& idaDET)+ &
& aDET(i,k,iband)
#endif
!
! Calculate CDOC absorption.
! NOTE: 12 factor is to convert ugrams per liter, and 0.001 converts
! to mg/liter. Specific absorption
! coefficients were calculated as m-1 / (mg DOC/liters sw).
! net factor = (12*0.001) = 0.012
!
cff=0.012_r8*(Bio(i,k,iCDMC(ilab))* &
& aDOC(ilab,iband)+ &
& Bio(i,k,iCDMC(irct))* &
& aDOC(irct,iband))
tot_ab(k,iband)=tot_ab(k,iband)+cff+awater(iband)
#ifdef DIAGNOSTICS_BIO
aCDC(i,k,iband)=aCDC(i,k,iband)+cff
DiaBio4d(i,j,k,iband,idaCDC)=DiaBio4d(i,j,k,iband, &
& idaCDC)+ &
& aCDC(i,k,iband)
#endif
!
! Calculate scattering and backscattering (see equation 19 Morel, 1991,
! Prog. Ocean). Morel, 1988 puts spectral dependency in backscattering.
! Since Morel (1991) does not have a backscattering equation, use 1988
! paper. Morel 2001 has slight adjustment 0.01, rather than 0.02.
! This was altered, but never tested in ROMS 1.8 on 03/08/03.
!
par_s=0.3_r8*(tChl**0.62_r8) ! scattering
par_b=0.0_r8 ! backscattering
IF (tChl.gt.0.0_r8) THEN
par_b=par_s*(0.002_r8+ &
& 0.01_r8* &
& (0.5_r8-0.25_r8*LOG10(tChl))* &
& wavedp(iband))
END IF
par_b=MAX(par_b,0.0_r8)
#ifdef DIAGNOSTICS_BIO
s_phy(i,k,iband)=s_phy(i,k,iband)+par_s
b_phy(i,k,iband)=b_phy(i,k,iband)+par_b
DiaBio4d(i,j,k,iband,idsPHY)=DiaBio4d(i,j,k,iband, &
& idsPHY)+ &
& s_phy(i,k,iband)
DiaBio4d(i,j,k,iband,idbPHY)=DiaBio4d(i,j,k,iband, &
& idbPHY)+ &
& b_phy(i,k,iband)
#endif
!
! However, for omega0 calculation, "par_s" must be spectral, so use
! dependency from Sathy and Platt 1988.
!
tot_s(k,iband)=bwater(iband)+par_s*wavedp(iband)
#ifdef DIAGNOSTICS_BIO
s_tot(i,k,iband)=s_tot(i,k,iband)+tot_s(k,iband)
#endif
!
! Morel, 1988 instead of 1991. See methods.
!
tot_b(k,iband)=0.5_r8*bwater(iband)+par_b
#ifdef DIAGNOSTICS_BIO
b_tot(i,k,iband)=b_tot(i,k,iband)+tot_b(k,iband)
DiaBio4d(i,j,k,iband,idsTOT)=DiaBio4d(i,j,k,iband, &
& idsTOT)+ &
& s_tot(i,k,iband)
DiaBio4d(i,j,k,iband,idbTOT)=DiaBio4d(i,j,k,iband, &
& idbTOT)+ &
& b_tot(i,k,iband)
#endif
#ifdef BIO_OPTIC
!
! Next statement is Eq. 11 of Lee et al. (2005) from Gallegos
! equal to SPECKD in Gallegos. Notice that the in-water solar
! angle in "cff1" is in degrees.
!
cff1=1.0_r8+ &
& 0.005_r8*ACOS(avgcos(k,iband))*rad2deg
cff2=4.18_r8*(1.0_r8-0.52_r8* &
& EXP(-10.8_r8*tot_ab(k,iband)))
dATT(k,iband)=cff1*tot_ab(k,iband)+ &
& cff2*tot_b (k,iband)
#else
!
! Sathy and Platt JGR 1988. This is set with the average cosine of
! the box above, and used to calculate a new avgcos for this level.
! This new average cosine is then used to recalculate the attenuation
! coefficient.
!
dATT(k,iband)=(tot_ab(k,iband)+ &
& tot_b (k,iband))/avgcos(k,iband)
#endif
!
! See Mobley, 1995 for graphical depiction of this equation.
!
avgcos_min=avgcos(k,iband)+ &
& (0.5_r8-avgcos(k,iband))* &
& (tot_s(k,iband)/ &
& (tot_ab(k,iband)+tot_s(k,iband)))
!
! Calculate average cosine. Linear fit to average cosine versus optical
! depth relationship. The FV1 calculation keeps the denominator of the
! slope calculation from going negative and above 1.
!
FV1=MAX(1.0_r8, &
& 7.0_r8-dATT(k,iband)*ABS(z_r(i,j,k)))
slope_AC =MIN(0.0_r8, &
& (avgcos_min-avgcos(k,iband))/FV1)
avgcos(k,iband)=avgcos(k,iband)+ &
& slope_AC*dATT(k,iband)*Hz(i,j,k)
#ifdef BIO_OPTIC
!
! Next statement is Eq. 11 of Lee et al. (2005). Notice that "cff1" is
! recomputed because "avgcos" changed and "cff2" is the same as above.
!
cff1=1.0_r8+ &
& 0.005_r8*ACOS(avgcos(k,iband))*rad2deg
dATT(k,iband)=cff1*tot_ab(k,iband)+ &
& cff2*tot_b (k,iband)
#else
dATT(k,iband)=(tot_ab(k,iband)+ &
& tot_b (k,iband))/avgcos(k,iband)
#endif
!
! Set avgcos for next level.
!
IF (k.ne.1) THEN
avgcos(k-1,iband)=avgcos(k,iband)
END IF
!
! Calculate spectral irradiance with depth.
!
FV1=dATT(k,iband)*Hz(i,j,k)
FV2=dATT_sum(iband)+0.5_r8*FV1
dATT_sum(iband)=dATT_sum(iband)+FV1
specir_d(k,iband)=SpecIr(i,j,iband)* &
& EXP(-FV2)*DLAM
!
! Calculate spectral scalar irradiance. Morel, 1991 Prog. Ocean.
!
specir_scal(i,k,iband)=specir_d(k,iband)* &
& (dATT(k,iband)/ &
& tot_ab(k,iband))
E0_nz(i,k)=E0_nz(i,k)+specir_scal(i,k,iband)
!
! Calculate Ed_nz.
!
Ed_nz(i,k)=Ed_nz(i,k)+specir_d(k,iband)
#ifdef DIAGNOSTICS_BIO
DiaBio3d(i,j,iband,idSpIr)=DiaBio3d(i,j,iband, &
& idSpIr)+ &
& SpecIr(i,j,iband)
DiaBio4d(i,j,k,iband,iddIrr)=DiaBio4d(i,j,k,iband, &
& iddIrr)+ &
& specir_d(k,iband)
DiaBio4d(i,j,k,iband,idsIrr)=DiaBio4d(i,j,k,iband, &
& idsIrr)+ &
& specir_scal(i,k,iband)
DiaBio4d(i,j,k,iband,idLatt)=DiaBio4d(i,j,k,iband, &
& idLatt)+ &
& dATT(k,iband)
DiaBio4d(i,j,k,iband,idAcos)=DiaBio4d(i,j,k,iband, &
& idAcos)+ &
& avgcos(k,iband)
#endif
END DO
Ed_tot=E0_nz(i,k)
!
! Set bottom of the euphotic zone.
!
Keuphotic(i)=k
END IF
END DO
END IF
END DO
!
!-----------------------------------------------------------------------
! Bacterial nutrient uptake.
!-----------------------------------------------------------------------
!
DO ibac=1,Nbac
DO k=1,N(ng)
DO i=Istr,Iend
!
! DOM uptake.
!
IF ((Bio(i,k,iDOMC(ilab)).gt.0.0_r8).and. &
& (Bio(i,k,iDOMN(ilab)).gt.0.0_r8).and. &
& (Bio(i,k,iDOMP(ilab)).gt.0.0_r8)) THEN
NupDOC_ba(i,k,ibac)=GtBAC(i,k,ibac)* &
& Bio(i,k,iBacC(ibac))* &
& I_Bac_Ceff(ng)* &
& (Bio(i,k,iDOMC(ilab))/ &
& (HsDOC_ba(ibac,ng)+ &
& Bio(i,k,iDOMC(ilab))))
NupDON_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)* &
& Bio(i,k,iDOMN(ilab))/ &
& Bio(i,k,iDOMC(ilab))
NupDOP_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)* &
& Bio(i,k,iDOMP(ilab))/ &
& Bio(i,k,iDOMC(ilab))
ELSE
NupDOC_ba(i,k,ibac)=0.0_r8
NupDON_ba(i,k,ibac)=0.0_r8
NupDOP_ba(i,k,ibac)=0.0_r8
END IF
totDOC_d(i,k)=totDOC_d(i,k)+NupDOC_ba(i,k,ibac)
totDON_d(i,k)=totDON_d(i,k)+NupDON_ba(i,k,ibac)
totDOP_d(i,k)=totDOP_d(i,k)+NupDOP_ba(i,k,ibac)
!
! NH4 uptake.
!
NupNH4_ba(i,k,ibac)=GtBAC(i,k,ibac)* &
& Bio(i,k,iBacN(ibac))* &
& Bio(i,k,iNH4_)/ &
& (HsNH4_ba(ibac,ng)+Bio(i,k,iNH4_))
totNH4_d(i,k)=totNH4_d(i,k)+NupNH4_ba(i,k,ibac)
!
! PO4 uptake.
!
NupPO4_ba(i,k,ibac)=GtBAC(i,k,ibac)* &
& Bio(i,k,iBacP(ibac))* &
& Bio(i,k,iPO4_)/ &
& (HsPO4_ba(ibac,ng)+Bio(i,k,iPO4_))
totPO4_d(i,k)=totPO4_d(i,k)+NupPO4_ba(i,k,ibac)
!
! Fe uptake.
!
NupFe_ba(i,k,ibac)=GtBAC(i,k,ibac)* &
& Bio(i,k,iBacF(ibac))* &
& Bio(i,k,iFeO_)/ &
& (HsFe_ba(ibac,ng)+Bio(i,k,iFeO_))
totFe_d(i,k)=totFe_d(i,k)+NupFe_ba(i,k,ibac)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
! Phytoplankton dark nutrient uptake.
!-----------------------------------------------------------------------
!
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
IF (C2nALG(i,k,iphy).gt.C2nALGminABS(iphy,ng)) THEN
!
! NOTE: these are being saved to test for total nutrient uptake.
! If nutrient uptake is greater than maximum nutrient, then
! each of the uptakes are reduced by their fractional contri-
! bution to the total.
!
Nup_max=GtALG(i,k,iphy)
NupNO3(i,k,iphy)=(Bio(i,k,iNO3_)/ &
& (HsNO3(iphy,ng)+Bio(i,k,iNO3_))* &
& EXP(-BET_(iphy,ng)*Bio(i,k,iNH4_)))
NupNH4(i,k,iphy)=Bio(i,k,iNH4_)/ &
& (HsNH4(iphy,ng)+Bio(i,k,iNH4_))
!
! Test that Wroblewski equation does not exceed 1.0.
!
FV1=NupNO3(i,k,iphy)+NupNH4(i,k,iphy)
IF (FV1.gt.1.0_r8) THEN
FV1=1.0_r8/FV1
NupNO3(i,k,iphy)=NupNO3(i,k,iphy)*FV1
NupNH4(i,k,iphy)=NupNH4(i,k,iphy)*FV1
END IF
!
! Change from percentage of maximum to mass per second.
!
FV1=Nup_max*Bio(i,k,iPhyN(iphy))
NupNO3(i,k,iphy)=NupNO3(i,k,iphy)*FV1
NupNH4(i,k,iphy)=NupNH4(i,k,iphy)*FV1
!
! Test for DON uptake.
!
IF (C2nALG(i,k,iphy).gt.C2nNupDON(iphy,ng)) THEN
NupDON(i,k,iphy)=FV1* &
& Bio(i,k,iDOMN(ilab))/ &
& (HsDON(iphy,ng)+ &
& Bio(i,k,iDOMN(ilab)))
END IF
!
! Accumulate total demand for nutrients.
!
totNO3_d(i,k)=totNO3_d(i,k)+NupNO3(i,k,iphy)
totNH4_d(i,k)=totNH4_d(i,k)+NupNH4(i,k,iphy)
totDON_d(i,k)=totDON_d(i,k)+NupDON(i,k,iphy)
END IF
!
! Dark silica uptake, min C2Si test.
! The LARGER test can be removed after testing phase.
!
IF (HsSiO(iphy,ng).lt.LARGER) THEN
IF (C2sALG(i,k,iphy).gt.C2SiALGminABS(iphy,ng)) THEN
Nup_max=GtALG(i,k,iphy)
NupSiO(i,k,iphy)=Bio(i,k,iSiO_)/ &
& (HsSiO(iphy,ng)+Bio(i,k,iSiO_))
!
! Change from percentage of maximum to mass per second.
!
IF (iPhyS(iphy).gt.0) THEN
FV1=Nup_max*Bio(i,k,iPhyS(iphy))
NupSiO(i,k,iphy)=NupSiO(i,k,iphy)*FV1
ELSE
NupSiO(i,k,iphy)=0.0_r8
END IF
!
! Accumulate total demand for nutrients.
!
totSiO_d(i,k)=totSiO_d(i,k)+NupSiO(i,k,iphy)
END IF
END IF
!
! Dark phophorus uptake, min C2P test.
! The LARGER test can be removed after testing phase.
!
IF (HsPO4(iphy,ng).lt.LARGER) THEN
IF (C2pALG(i,k,iphy).gt.C2pALGminABS(iphy,ng)) THEN
Nup_max=GtALG(i,k,iphy)
NupPO4(i,k,iphy)=Bio(i,k,iPO4_)/ &
& (HsPO4(iphy,ng)+Bio(i,k,iPO4_))
!
! Change from percentage of maximum to mass per second.
!
FV1=Nup_max*Bio(i,k,iPhyP(iphy))
NupPO4(i,k,iphy)=NupPO4(i,k,iphy)*FV1
!
! Test for alk. phosphatase
!
IF (C2pALG(i,k,iphy).gt.C2pALKPHOS(iphy,ng)) THEN
NupDOP(i,k,iphy)=FV1* &
Bio(i,k,iDOMP(ilab))/ &
& (HsDOP(iphy,ng)+ &
& Bio(i,k,iDOMP(ilab)))
END IF
!
! Accumulate total demand for nutrients.
!
totPO4_d(i,k)=totPO4_d(i,k)+NupPO4(i,k,iphy)
totDOP_d(i,k)=totDOP_d(i,k)+NupDOP(i,k,iphy)
END IF
END IF
!
! Dark iron uptake, min C2Fe test.
! The LARGER test can be removed after testing phase.
!
IF (HsFe(iphy,ng).lt.LARGER) THEN
IF (C2fALG(i,k,iphy).gt.C2FeALGminABS(iphy,ng)) THEN
Nup_max=GtALG(i,k,iphy)
NupFe(i,k,iphy)=Bio(i,k,iFeO_)/ &
& (HsFe(iphy,ng)+Bio(i,k,iFeO_))
!
! Change from percentage of maximum to mass per second.
!
FV1=Nup_max*Bio(i,k,iPhyF(iphy))
NupFe(i,k,iphy)=NupFe(i,k,iphy)*FV1
!
! Accumulate total demand for nutrients.
!
totFe_d(i,k)=totFe_d(i,k)+NupFe(i,k,iphy)
END IF
END IF
END DO
END DO
END DO
!
! Calculate bacterial nitrification as a Michaelis-Menton function
! of ambient NH4 concentration, beneath the euphotic zone (light
! inhibits nitrification).
!
DO k=1,N(ng)
DO i=Istr,Iend
NitrBAC(i,k)=0.0_r8
NH4toNO3(i,k)=0.0_r8
NtoNBAC(i,k)=0.0_r8
NtoPBAC(i,k)=0.0_r8
NtoFeBAC(i,k)=0.0_r8
IF (k.lt.Keuphotic(i)) THEN
NH4toNO3(i,k)=RtNIT(ng)* &
& Bio(i,k,iNH4_)/(HsNIT(ng)+Bio(i,k,iNH4_))
!
! Nitrification fixes DIC into POC.
! Conversion factor of 7.0 from Kaplan 1983 "Nitrogen in the Sea"
! factor equals (1.0 / (7.0 * C2nBAC)). Adds NH4 uptake as biomass.
!
NitrBAC(i,k)=NH4toNO3(i,k)/7.0_r8
NtoNBAC(i,k)=NitrBAC(i,k)*N2cBAC(ng)
NtoPBAC(i,k)=NitrBAC(i,k)*P2cBAC(ng)
NtoFeBAC(i,k)=NitrBAC(i,k)*Fe2cBAC(ng)
totNH4_d(i,k)=totNH4_d(i,k)+NH4toNO3(i,k)+NtoNBAC(i,k)
totPO4_d(i,k)=totPO4_d(i,k)+NtoPBAC(i,k)
totFe_d (i,k)=totFe_d (i,k)+NtoFeBAC(i,k)
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
! Test that total nutrient demand does not exceed supply. If it does
! total demand is normalized to the total supply. Each species demand
! is reduced to its weighted average percentage of the supply.
!-----------------------------------------------------------------------
!
DO k=1,N(ng)
DO i=Istr,Iend
FV2=totNO3_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iNO3_)) THEN
FV1=(Bio(i,k,iNO3_)-VSMALL)/FV2
DO iphy=1,Nphy
NupNO3(i,k,iphy)=NupNO3(i,k,iphy)*FV1
END DO
END IF
!
FV2=totNH4_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iNH4_)) THEN
FV1=(Bio(i,k,iNH4_)-VSMALL)/FV2
DO iphy=1,Nphy
NupNH4(i,k,iphy)=NupNH4(i,k,iphy)*FV1
END DO
DO ibac=1,Nbac
NupNH4_ba(i,k,ibac)=NupNH4_ba(i,k,ibac)*FV1
END DO
NH4toNO3(i,k)=NH4toNO3(i,k)*FV1
NtoNBAC(i,k)=NtoNBAC(i,k)*FV1
END IF
!
FV2=totSiO_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iSiO_)) THEN
FV1=(Bio(i,k,iSiO_)-VSMALL)/FV2
DO iphy=1,Nphy
NupSiO(i,k,iphy)=NupSiO(i,k,iphy)*FV1
END DO
END IF
!
FV2=totPO4_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iPO4_)) THEN
FV1=(Bio(i,k,iPO4_)-VSMALL)/FV2
DO iphy=1,Nphy
NupPO4(i,k,iphy)=NupPO4(i,k,iphy)*FV1
END DO
DO ibac=1,Nbac
NupPO4_ba(i,k,ibac)=NupPO4_ba(i,k,ibac)*FV1
END DO
NtoPBAC(i,k)=NtoPBAC(i,k)*FV1
END IF
!
FV2=totFe_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iFeO_)) THEN
FV1=(Bio(i,k,iFeO_)-VSMALL)/FV2
DO iphy=1,Nphy
NupFe(i,k,iphy)=NupFe(i,k,iphy)*FV1
END DO
DO ibac=1,Nbac
NupFe_ba(i,k,ibac)=NupFe_ba(i,k,ibac)*FV1
END DO
NtoFeBAC(i,k)=NtoFeBAC(i,k)*FV1
END IF
!
! Bacteria are the only group to take up DOC. Remove BAC DON and
! BAC DOP uptake from total uptake; adjust uptake and add back.
!
FV2=totDOC_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iDOMC(ilab))) THEN
FV1=(Bio(i,k,iDOMC(ilab))-VSMALL)/FV2
totDOC_d(i,k)=totDOC_d(i,k)*FV1
DO ibac=1,Nbac
NupDOC_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)*FV1
totDON_d(i,k)=totDON_d(i,k)-NupDON_ba(i,k,ibac)
NupDON_ba(i,k,ibac)=NupDON_ba(i,k,ibac)*FV1
totDON_d(i,k)=totDON_d(i,k)+NupDON_ba(i,k,ibac)
totDOP_d(i,k)=totDOP_d(i,k)-NupDOP_ba(i,k,ibac)
NupDOP_ba(i,k,ibac)=NupDOP_ba(i,k,ibac)*FV1
totDOP_d(i,k)=totDOP_d(i,k)+NupDOP_ba(i,k,ibac)
END DO
END IF
!
! Remove BAC DON uptake from total uptake; adjust uptake and add back.
!
FV2=totDON_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iDOMN(ilab))) THEN
FV1=(Bio(i,k,iDOMN(ilab))-VSMALL)/FV2
totDON_d(i,k)=totDON_d(i,k)*FV1
totDOC_d(i,k)=totDOC_d(i,k)*FV1
DO iphy=1,Nphy
NupDON(i,k,iphy)=NupDON(i,k,iphy)*FV1
END DO
DO ibac=1,Nbac
NupDON_ba(i,k,ibac)=NupDON_ba(i,k,ibac)*FV1
NupDOC_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)*FV1
totDOP_d(i,k)=totDOP_d(i,k)-NupDOP_ba(i,k,ibac)
NupDOP_ba(i,k,ibac)=NupDOP_ba(i,k,ibac)*FV1
totDOP_d(i,k)=totDOP_d(i,k)+NupDOP_ba(i,k,ibac)
END DO
END IF
!
! Remove BAC DOP uptake from total uptake; adjust uptake and add back.
!
FV2=totDOP_d(i,k)*dtbio
IF (FV2.gt.Bio(i,k,iDOMP(ilab))) THEN
FV1=(Bio(i,k,iDOMP(ilab))-VSMALL)/FV2
totDOP_d(i,k)=totDOP_d(i,k)*FV1
totDOC_d(i,k)=totDOC_d(i,k)*FV1
DO iphy=1,Nphy
NupDOP(i,k,iphy)=NupDOP(i,k,iphy)*FV1
END DO
DO ibac=1,Nbac
NupDOP_ba(i,k,ibac)=NupDOP_ba(i,k,ibac)*FV1
totDON_d(i,k)=totDON_d(i,k)-NupDON_ba(i,k,ibac)
NupDON_ba(i,k,ibac)=NupDON_ba(i,k,ibac)*FV1
totDON_d(i,k)=totDON_d(i,k)+NupDON_ba(i,k,ibac)
NupDOC_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)*FV1
END DO
END IF
END DO
END DO
!
! Increase particulate nutrients by the amount of the uptake.
!
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
Bio_new(i,k,iPhyN(iphy))=Bio_new(i,k,iPhyN(iphy))+ &
& NupNO3(i,k,iphy)+ &
& NupNH4(i,k,iphy)+ &
& NupDON(i,k,iphy)
Bio_new(i,k,iPhyP(iphy))=Bio_new(i,k,iPhyP(iphy))+ &
& NupPO4(i,k,iphy)+ &
& NupDOP(i,k,iphy)
Bio_new(i,k,iPhyF(iphy))=Bio_new(i,k,iPhyF(iphy))+ &
& NupFe(i,k,iphy)
IF (iPhyS(iphy).gt.0) THEN
Bio_new(i,k,iPhyS(iphy))=Bio_new(i,k,iPhyS(iphy))+ &
& NupSiO(i,k,iphy)
END IF
!
! Update nutrient arrays for growth and budgets. Bacterial uptake
! included below.
!
Bio_new(i,k,iNO3_)=Bio_new(i,k,iNO3_)- &
& NupNO3(i,k,iphy)
Bio_new(i,k,iNH4_)=Bio_new(i,k,iNH4_)- &
& NupNH4(i,k,iphy)
Bio_new(i,k,iSiO_)=Bio_new(i,k,iSiO_)- &
& NupSiO(i,k,iphy)
Bio_new(i,k,iPO4_)=Bio_new(i,k,iPO4_)- &
& NupPO4(i,k,iphy)
Bio_new(i,k,iFeO_)=Bio_new(i,k,iFeO_)- &
& NupFe (i,k,iphy)
Bio_new(i,k,iDOMN(ilab))=Bio_new(i,k,iDOMN(ilab))- &
& NupDON(i,k,iphy)
Bio_new(i,k,iDOMP(ilab))=Bio_new(i,k,iDOMP(ilab))- &
& NupDOP(i,k,iphy)
END DO
END DO
END DO
!
! Nitrification fixes DIC into DOC.
!
DO k=1,N(ng)
DO i=Istr,Iend
Bio_new(i,k,iDIC_)=Bio_new(i,k,iDIC_)- &
& NitrBAC(i,k)
END DO
END DO
!
! Add nitrifying bacteria biomass to heterotrophic bacteria biomass.
! Adding PON, POP, POFe to BacC arrays at current C2_BAC ratios.
!
DO ibac=1,Nbac
DO k=1,N(ng)
DO i=Istr,Iend
Bio_new(i,k,iBacC(ibac))=Bio_new(i,k,iBacC(ibac))+ &
& NitrBAC(i,k)
Bio_new(i,k,iBacN(ibac))=Bio_new(i,k,iBacN(ibac))+ &
& NtoNBAC(i,k)
Bio_new(i,k,iBacP(ibac))=Bio_new(i,k,iBacP(ibac))+ &
& NtoPBAC(i,k)
Bio_new(i,k,iBacF(ibac))=Bio_new(i,k,iBacF(ibac))+ &
& NtoFeBAC(i,k)
END DO
END DO
END DO
!
! Update nutrient arrays for nitrification.
!
DO k=1,N(ng)
DO i=Istr,Iend
Bio_new(i,k,iNO3_)=Bio_new(i,k,iNO3_)+ &
& NH4toNO3(i,k)
Bio_new(i,k,iNH4_)=Bio_new(i,k,iNH4_)- &
& (NH4toNO3(i,k)+NtoNBAC(i,k))
Bio_new(i,k,iPO4_)=Bio_new(i,k,iPO4_)- &
& NtoPBAC(i,k)
Bio_new(i,k,iFeO_)=Bio_new(i,k,iFeO_)- &
& NtoFeBAC(i,k)
END DO
END DO
!
!-----------------------------------------------------------------------
! Light mediated carbon growth.
!-----------------------------------------------------------------------
!
DO i=Istr,Iend
DO k=N(ng),Keuphotic(i),-1
DO iphy=1,Nphy
IF (Bio(i,k,iPhyC(iphy)).gt.0.0_r8) THEN
!
! Calculate weighted average spectral absorption.
!
aPHYN_wa=0.0_r8
DO iband=1,NBands
aPHYN_wa=aPHYN_wa+(aPHYN_al(i,k,iphy,iband)* &
& specir_scal(i,k,iband))
END DO
!
! If Keuphotic(i) < N+1, and E0_nz(i,k)=0, this will cause pigments to
! blow up. This should never happen, unless Keuphotic is not calcuated
! properly. WPB
!
aPHYN_wa=aPHYN_wa/E0_nz(i,k)
!
! Calculate "alfa" for HTAN function of P vs. I.
! (conversion: Ein/microEin * 10e3)
!
alfa(k,iphy)=(aPHYN_wa/Bio(i,k,iPhyC(iphy)))* &
& qu_yld(iphy,ng)*0.001_r8
!
! Light limited growth rate.
!
FV1=MAX(0.0_r8,E0_nz(i,k)-E0_comp(iphy,ng))
FV2=E0_nz(i,k)-E0_inhib(iphy,ng)
IF (FV2.gt.0.0_r8) THEN
Gt_ll(k,iphy)=GtALG(i,k,iphy)* &
& TANH(alfa(k,iphy)*FV1/ &
& GtALG(i,k,iphy))* &
& EXP(-inhib_fac(iphy,ng)*FV2)
ELSE
Gt_ll(k,iphy)=GtALG(i,k,iphy)* &
& TANH(alfa(k,iphy)*FV1/ &
& GtALG(i,k,iphy))
END IF
!
! Nutrient limited growth rates.
!
! REMEMBER that sinking speed to be set by gradient of limiting
! nutrient, allowing for negative sinking. Try storing growth
! rate terms in an array and using MAXLOC for if test.
!
! Nitrogen limited growth rate.
!
IF (Bio(i,k,iPhyN(iphy)).gt.0.0_r8) THEN
FV1=Bio(i,k,iPhyC(iphy))/ &
& (Bio(i,k,iPhyN(iphy))+Bio_new(i,k,iPhyN(iphy)))
Gt_nl(k,iphy)=mu_bar_n(i,k,iphy)* &
& (1.0_r8-ImaxC2nALG(iphy,ng)*FV1)
Gt_nl(k,iphy)=MAX(0.0_r8, &
& MIN(Gt_nl(k,iphy), &
& GtALG(i,k,iphy)))
END IF
!
! Silica limited growth rate.
! Testing for silica incorporation.
!
IF (iPhyS(iphy).gt.0) THEN
IF ((HsSiO(iphy,ng).lt.LARGER).and. &
& (Bio(i,k,iPhyS(iphy)).gt.0.0_r8)) THEN
FV1=Bio(i,k,iPhyC(iphy))/ &
& (Bio(i,k,iPhyS(iphy))+ &
& Bio_new(i,k,iPhyS(iphy)))
Gt_sl(k,iphy)=mu_bar_s(i,k,iphy)* &
& (1.0_r8-ImaxC2SiALG(iphy,ng)*FV1)
Gt_sl(k,iphy)=MAX(0.0_r8, &
& MIN(Gt_sl(k,iphy), &
& GtALG(i,k,iphy)))
ELSE
Gt_sl(k,iphy)=LARGER
END IF
ELSE
Gt_sl(k,iphy)=LARGER
END IF
!
! Phosphorus limited growth rate.
!
IF ((HsPO4(iphy,ng).lt.LARGER).and. &
& (Bio(i,k,iPhyP(iphy)).gt.0.0_r8)) THEN
FV1=Bio(i,k,iPhyC(iphy))/ &
& (Bio(i,k,iPhyP(iphy))+Bio_new(i,k,iPhyP(iphy)))
Gt_pl(k,iphy)=mu_bar_p(i,k,iphy)* &
& (1.0_r8-ImaxC2pALG(iphy,ng)*FV1)
Gt_pl(k,iphy)=MAX(0.0_r8, &
& MIN(Gt_pl(k,iphy), &
& GtALG(i,k,iphy)))
ELSE
Gt_pl(k,iphy)=LARGER
END IF
!
! Iron limited growth rate
!
IF ((HsFe(iphy,ng).lt.LARGER).and. &
& (Bio(i,k,iPhyF(iphy)).gt.0.0_r8)) THEN
FV1=Bio(i,k,iPhyC(iphy))/ &
& (Bio(i,k,iPhyF(iphy))+Bio_new(i,k,iPhyF(iphy)))
Gt_fl(k,iphy)=mu_bar_f(i,k,iphy)* &
& (1.0_r8-ImaxC2FeALG(iphy,ng)*FV1)
Gt_fl(k,iphy)=MAX(0.0_r8, &
& MIN(Gt_fl(k,iphy), &
& GtALG(i,k,iphy)))
ELSE
Gt_fl(k,iphy)=LARGER
END IF
!
! Realized growth rate is minimum of light or nutrient limited rate.
!
GtALG_r(i,k,iphy)=MIN(Gt_ll(k,iphy),Gt_nl(k,iphy), &
& Gt_sl(k,iphy),Gt_pl(k,iphy), &
& Gt_fl(k,iphy))
IF (GtALG_r(i,k,iphy).ge.LARGER) THEN
GtALG_r(i,k,iphy)=0.0_r8
END IF
!
! Carbon growth calculations.
!
FV1=Bio(i,k,iPhyC(iphy))*GtALG_r(i,k,iphy)
Bio_new(i,k,iPhyC(iphy))=Bio_new(i,k,iPhyC(iphy))+ &
& FV1
Bio_new(i,k,iDIC_)=Bio_new(i,k,iDIC_)- &
& FV1
!
! Pigment growth calculations.
!
DO ipig=1,Npig
IF (iPigs(iphy,ipig).gt.0) THEN
itrc=iPigs(iphy,ipig)
IF (Bio(i,k,iPhyC(iphy)).gt.0.0_r8) THEN
FV1=Bio(i,k,itrc)*GtALG_r(i,k,iphy)
Bio_new(i,k,itrc)=Bio_new(i,k,itrc)+FV1
END IF
END IF
END DO
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
! Bacterioplankton carbon growth terms.
!-----------------------------------------------------------------------
!
DO k=1,N(ng)
DO i=Istr,Iend
Het_BAC=0.0_r8
RelDOC1=0.0_r8
RelDON1=0.0_r8
RelDOP1=0.0_r8
RelFe=0.0_r8
!
! NOTE: Only DOC2/DON2 formation is in this section.
! Take colored excretion off the top. 03/18/00
! also, not excreting any DOP or Fe
! REMEMBER, if excreting DOP and Fe, must address changes in growth if
! tests. (see DON equations). 03/21/00.
!
DO ibac=1,Nbac
FV1=NupDOC_ba(i,k,ibac)*ExBAC_c(ng)* &
& (1.0_r8-cDOCfrac_c(irct,ng))
FV2=NupDOC_ba(i,k,ibac)*ExBAC_c(ng)* &
& cDOCfrac_c(irct,ng)
FV3=NupDON_ba(i,k,ibac)*ExBAC_n(ng)
!
Bio_new(i,k,iDOMC(irct))=Bio_new(i,k,iDOMC(irct))+ &
& FV1
Bio_new(i,k,iCDMC(irct))=Bio_new(i,k,iCDMC(irct))+ &
& FV2
Bio_new(i,k,iDOMN(irct))=Bio_new(i,k,iDOMN(irct))+ &
& FV3
!
! As we are taking it off the top, must remove from DOMN1 now. No other
! organisms use DOMC1, so net term (totDOC_d) can be used in budgeting
! below. This saves cycles, but makes code difficult to read. WPB
!
Bio_new(i,k,iDOMN(ilab))=Bio_new(i,k,iDOMN(ilab))- &
& FV3
!
! Remove from uptake.
!
NupDOC_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)- &
& (FV1+FV2)
NupDON_ba(i,k,ibac)=NupDON_ba(i,k,ibac)- &
& FV3
!
! Determine growth limitation. Assuming 100% efficiency for N, P, Fe.
! If DOMC=0, or DOMN=0, or DOMP=0, then NupDOC_ba = NupDON_ba =
! NupDOP_ba = 0 and none of the divisions below are accessed. WPB
!
Bac_G(1)=NupDOC_ba(i,k,ibac)*Bac_Ceff(ng)
Bac_G(2)=(NupDON_ba(i,k,ibac)+ &
& NupNH4_ba(i,k,ibac))* &
& C2nBAC(ng)
Bac_G(3)=(NupDOP_ba(i,k,ibac)+ &
& NupPO4_ba(i,k,ibac))* &
& C2pBAC(ng)
Bac_G(4)=NupFe_ba(i,k,ibac)*C2FeBAC(ng)
!
! Energy limited case. All excess nutrients returned in inorganic form.
!
IF ((Bac_G(1).le.Bac_G(2)).and. &
& (Bac_G(1).le.Bac_G(3)).and. &
& (Bac_G(1).le.Bac_G(4))) THEN
Het_BAC=Bac_G(1)
FV1=Bac_G(1)*N2cBAC(ng)
FV2=Bac_G(1)*P2cBAC(ng)
FV3=Bac_G(1)*Fe2cBAC(ng)
Bio_new(i,k,iBacN(ibac))=Bio_new(i,k,iBacN(ibac))+ &
& FV1
Bio_new(i,k,iBacP(ibac))=Bio_new(i,k,iBacP(ibac))+ &
& FV2
Bio_new(i,k,iBacF(ibac))=Bio_new(i,k,iBacF(ibac))+ &
& FV3
!
! Uptake arrays should probably now be negative. If NH4 or PO4 is
! positive, then there is some uptake of inorganic forms, but this
! value will be less than the original Nup value because of IF test.
!
NupNH4_ba(i,k,ibac)=FV1-NupDON_ba(i,k,ibac)
NupPO4_ba(i,k,ibac)=FV2-NupDOP_ba(i,k,ibac)
!
! Because Fe is considered to be all inorganic, only net uptake of Fe
! is needed.
!
RelFe=NupFe_ba(i,k,ibac)-FV3
NupFe_ba(i,k,ibac)=FV3
!
! Nitrogen limited case. Excess nutrients returned in organic form
! first, inorganic second.
!
ELSE IF ((Bac_G(2).le.Bac_G(3)).and. &
& (Bac_G(2).le.Bac_G(4))) THEN
Het_BAC=Bac_G(2)
FV2=Bac_G(2)*P2cBAC(ng)
FV3=Bac_G(2)*Fe2cBAC(ng)
Bio_new(i,k,iBacN(ibac))=Bio_new(i,k,iBacN(ibac))+ &
& (NupDON_ba(i,k,ibac)+ &
& NupNH4_ba(i,k,ibac))
Bio_new(i,k,iBacP(ibac))=Bio_new(i,k,iBacP(ibac))+ &
& FV2
Bio_new(i,k,iBacF(ibac))=Bio_new(i,k,iBacF(ibac))+ &
& FV3
!
! Uptake arrays will now reflect release of inorganic and organic
! revision of uptake.
!
FV1=(Bac_G(1)-Bac_G(2))*I_Bac_Ceff(ng)
NupDOC_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)-FV1
RelDOC1=FV1
!
! To get accurate DOP from C2pDOC, must add back excreted DOC.
!
FV4=FV1*R_ExBAC_c(ng)* &
!! & DOC_frac(i,k)* &
& Bio(i,k,iDOMP(ilab))/ &
& Bio(i,k,iDOMC(ilab))
FV5=FV2-(NupDOP_ba(i,k,ibac)+ &
NupPO4_ba(i,k,ibac)-FV4)
!
! If FV5 is positive then released DOP is required for bacteria growth.
!
IF (FV5.lt.0.0_r8) THEN
RelDOP1=FV4
NupPO4_ba(i,k,ibac)=NupPO4_ba(i,k,ibac)+FV5
ELSE
RelDOP1=FV4-FV5
END IF
NupDOP_ba(i,k,ibac)=NupDOP_ba(i,k,ibac)-RelDOP1
!
! Release Fe.
!
RelFe=NupFe_ba(i,k,ibac)-FV3
NupFe_ba(i,k,ibac)=FV3
!
! Phosphorous limited case. Excess nutrients returned in organic form
! first, inorganic second.
!
ELSE IF (Bac_G(3).le.Bac_G(4)) THEN
Het_BAC=Bac_G(3)
FV2=Bac_G(3)*N2cBAC(ng)
FV3=Bac_G(3)*Fe2cBAC(ng)
Bio_new(i,k,iBacN(ibac))=Bio_new(i,k,iBacN(ibac))+ &
& FV2
Bio_new(i,k,iBacP(ibac))=Bio_new(i,k,iBacP(ibac))+ &
& (NupDOP_ba(i,k,ibac)+ &
& NupPO4_ba(i,k,ibac))
Bio_new(i,k,iBacF(ibac))=Bio_new(i,k,iBacF(ibac))+ &
& FV3
!
! Uptake arrays will now reflect release of inorganic and organic
! revision of uptake.
!
FV1=(Bac_G(1)-Bac_G(3))*I_Bac_Ceff(ng)
NupDOC_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)-FV1
RelDOC1=FV1
!
! To get accurate DON from C2nDOC, must add back excreted DOC.
!
FV4=FV1*R_ExBAC_c(ng)* &
!! & DOC_frac(i,k)* &
& (Bio(i,k,iDOMN(ilab))/ &
& Bio(i,k,iDOMC(ilab)))*Frac_ExBAC_n(ng)
FV5=FV2-(NupDON_ba(i,k,ibac)+ &
& NupNH4_ba(i,k,ibac)-FV4)
!
! If FV5 is positive then released DON is required for bacteria growth.
!
IF (FV5.lt.0.0_r8) THEN
RelDON1=FV4
NupNH4_ba(i,k,ibac)=NupNH4_ba(i,k,ibac)+FV5
ELSE
RelDON1=FV4-FV5
END IF
NupDON_ba(i,k,ibac)=NupDON_ba(i,k,ibac)-RelDON1
!
! Release Fe.
!
RelFe=NupFe_ba(i,k,ibac)-FV3
NupFe_ba(i,k,ibac)=FV3
!
! Fe limited case. Excess nutrients returned in organic form
! first, inorganic second.
!
ELSE
Het_BAC=Bac_G(4)
FV2=Bac_G(4)*N2cBAC(ng)
FV3=Bac_G(4)*P2cBAC(ng)
Bio_new(i,k,iBacN(ibac))=Bio_new(i,k,iBacN(ibac))+ &
& FV2
Bio_new(i,k,iBacP(ibac))=Bio_new(i,k,iBacP(ibac))+ &
& FV3
Bio_new(i,k,iBacF(ibac))=Bio_new(i,k,iBacF(ibac))+ &
& NupFe_ba(i,k,ibac)
!
! Uptake arrays will now reflect release of inorganic and organic
! revision of uptake.
!
FV1=(Bac_G(1)-Bac_G(4))*I_Bac_Ceff(ng)
NupDOC_ba(i,k,ibac)=NupDOC_ba(i,k,ibac)-FV1
RelDOC1=FV1
!
! To get accurate DON from C2nDOC, must add back excreted DOC.
!
FV4=FV1*R_ExBAC_c(ng)* &
!! & DOC_frac(i,k)* &
& Bio(i,k,iDOMN(ilab))/ &
& Bio(i,k,iDOMC(ilab))*Frac_ExBAC_n(ng)
FV5=FV2-(NupDON_ba(i,k,ibac)+ &
& NupNH4_ba(i,k,ibac)-FV4)
!
! If FV5 is positive then released DON is required for bacteria growth.
!
IF (FV5.lt.0.0_r8) THEN
RelDON1=FV4
NupNH4_ba(i,k,ibac)=NupNH4_ba(i,k,ibac)+FV5
ELSE
RelDON1=FV4-FV5
END IF
NupDON_ba(i,k,ibac)=NupDON_ba(i,k,ibac)-RelDON1
!
! To get accurate DOP from C2pDOC, must add back excreted DOC.
!
FV4=FV1*R_ExBAC_c(ng)* &
!! & DOC_frac(i,k)* &
& Bio(i,k,iDOMP(ilab))/ &
& Bio(i,k,iDOMC(ilab))
FV5=FV2-(NupDOP_ba(i,k,ibac)+ &
& NupPO4_ba(i,k,ibac)-FV4)
!
! If FV5 is positive then released DOP is required for bacteria growth.
!
IF (FV5.lt.0.0_r8) THEN
RelDOP1=FV4
NupPO4_ba(i,k,ibac)=NupPO4_ba(i,k,ibac)+FV5
ELSE
RelDOP1=FV4-FV5
END IF
NupDOP_ba(i,k,ibac)=NupDOP_ba(i,k,ibac)-RelDOP1
END IF
!
! Increment nutrient arrays.
!
Bio_new(i,k,iBacC(ibac))=Bio_new(i,k,iBacC(ibac))+ &
& Het_BAC
FV1=NupDOC_ba(i,k,ibac)-Het_BAC
Bio_new(i,k,iDIC_)=Bio_new(i,k,iDIC_)+ &
& FV1
!
! NOTE: to be strictly accurate we should remove RelDOC1 from DOCNP1,
! and then add it back, since NupDOC_ba is a net term. This should
! wash out in the budgeting.
!
Bio_new(i,k,iDOMC(ilab))=Bio_new(i,k,iDOMC(ilab))- &
& (totDOC_d(i,k)-RelDOC1)
!! & (totDOC_d(i,k)-RelDOC1)* &
!! & DOC_frac(i,k)
!! Bio_new(i,k,iCDMC(ilab))=Bio_new(i,k,iCDMC(ilab))- &
!! & (totDOC_d(i,k)-RelDOC1)* &
!! & (1.0_r8-DOC_frac(i,k))
!!
! This is inclusive of RelDOX1, excretion of DON1 removed above.
!
Bio_new(i,k,iDOMN(ilab))=Bio_new(i,k,iDOMN(ilab))- &
& NupDON_ba(i,k,ibac)
Bio_new(i,k,iDOMP(ilab))=Bio_new(i,k,iDOMP(ilab))- &
& NupDOP_ba(i,k,ibac)
Bio_new(i,k,iNH4_)=Bio_new(i,k,iNH4_)- &
& NupNH4_ba(i,k,ibac)
Bio_new(i,k,iPO4_)=Bio_new(i,k,iPO4_)- &
& NupPO4_ba(i,k,ibac)
Bio_new(i,k,iFeO_)=Bio_new(i,k,iFeO_)- &
& NupFe_ba(i,k,ibac)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
! Phytoplankton Losses.
!-----------------------------------------------------------------------
!
DO iphy=1,Nphy
DO k=1,N(ng)
DO i=Istr,Iend
!
! Excretion.
!
IF ((C2nALG(i,k,iphy).ge. &
& C2nALGminABS(iphy,ng)).and. &
& (C2pALG(i,k,iphy).ge. &
& C2pALGminABS(iphy,ng)).and. &
& (HsSiO(iphy,ng).gt.LARGER)) THEN
FV1=Bio(i,k,iPhyC(iphy))*ExALG(iphy,ng)
Bio_new(i,k,iPhyC(iphy))=Bio_new(i,k,iPhyC(iphy))- &
& FV1
!
! No excretion of CDOC.
!
Bio_new(i,k,iDOMC(ilab))=Bio_new(i,k,iDOMC(ilab))+ &
& FV1
ELSE IF ((C2nALG(i,k,iphy).ge. &
& C2nALGminABS(iphy,ng)).and. &
& (C2pALG(i,k,iphy).ge. &
& C2pALGminABS(iphy,ng)).and. &
& (C2sALG(i,k,iphy).ge. &
& C2SiALGminABS(iphy,ng))) THEN
FV1=Bio(i,k,iPhyC(iphy))*ExALG(iphy,ng)
Bio_new(i,k,iPhyC(iphy))=Bio_new(i,k,iPhyC(iphy))- &
& FV1
!
! No excretion of CDOC.
!
Bio_new(i,k,iDOMC(ilab))=Bio_new(i,k,iDOMC(ilab))+ &
& FV1
END IF
!
! Grazing.
!
IF (Bio(i,k,iPhyC(iphy)).gt.refuge(i,k,iphy)) THEN
!
! Carbon calculations.
!
FV1=graz_act(i,k,iphy)*Bio(i,k,iPhyC(iphy))
Bio_new(i,k,iPhyC(iphy))=Bio_new(i,k,iPhyC(iphy))- &
& FV1
Bio_new(i,k,iFecC(isfc))=Bio_new(i,k,iFecC(isfc))+ &
& FecPEL(iphy,isfc,ng)*FV1
Bio_new(i,k,iFecC(iffc))=Bio_new(i,k,iFecC(iffc))+ &
& FecPEL(iphy,iffc,ng)*FV1
FV3=FecDOC(iphy,ng)*FV1
Bio_new(i,k,iDOMC(ilab))=Bio_new(i,k,iDOMC(ilab))+ &
& (1.0_r8-cDOCfrac_c(ilab,ng))*&
& FV3
Bio_new(i,k,iCDMC(ilab))=Bio_new(i,k,iCDMC(ilab))+ &
& cDOCfrac_c(ilab,ng)*FV3
Bio_new(i,k,iDIC_)=Bio_new(i,k,iDIC_)+ &
& FecCYC(iphy,ng)*FV1
!
! Nitrogen calculations.
!
FV2=graz_act(i,k,iphy)*Bio(i,k,iPhyN(iphy))
Bio_new(i,k,iPhyN(iphy))=Bio_new(i,k,iPhyN(iphy))- &
& FV2
Bio_new(i,k,iFecN(isfc))=Bio_new(i,k,iFecN(isfc))+ &
& FecPEL(iphy,isfc,ng)*FV2
Bio_new(i,k,iFecN(iffc))=Bio_new(i,k,iFecN(iffc))+ &
& FecPEL(iphy,iffc,ng)*FV2
Bio_new(i,k,iDOMN(ilab))=Bio_new(i,k,iDOMN(ilab))+ &
& FecDOC(iphy,ng)*FV2
Bio_new(i,k,iNH4_)=Bio_new(i,k,iNH4_)+ &
& FecCYC(iphy,ng)*FV2
!
! Silica calculations.
!
IF (iPhyS(iphy).gt.0) THEN
FV2=graz_act(i,k,iphy)*Bio(i,k,iPhyS(iphy))
Bio_new(i,k,iPhyS(iphy))=Bio_new(i,k,iPhyS(iphy))- &
& FV2
!
! Assuming that the fraction of material lost via sloppy feeding/cell
! lysis also results in silica tests being put into FecS pool.
!
Bio_new(i,k,iFecS(isfc))=Bio_new(i,k,iFecS(isfc))+ &
& FecDOC(iphy,ng)*FV2
Bio_new(i,k,iFecS(iffc))=Bio_new(i,k,iFecS(iffc))+ &
& (1.0_r8-FecDOC(iphy,ng))* &
& FV2
END IF
!
! Phosphorus calculations.
!
FV2=graz_act(i,k,iphy)*Bio(i,k,iPhyP(iphy))
Bio_new(i,k,iPhyP(iphy))=Bio_new(i,k,iPhyP(iphy))- &
& FV2
Bio_new(i,k,iFecP(isfc))=Bio_new(i,k,iFecP(isfc))+ &
& FecPEL(iphy,isfc,ng)*FV2
Bio_new(i,k,iFecP(iffc))=Bio_new(i,k,iFecP(iffc))+ &
& FecPEL(iphy,iffc,ng)*FV2
Bio_new(i,k,iDOMP(ilab))=Bio_new(i,k,iDOMP(ilab))+ &
& FecDOC(iphy,ng)*FV2
Bio_new(i,k,iPO4_)=Bio_new(i,k,iPO4_)+ &
& FecCYC(iphy,ng)*FV2
!
! Iron calculations. Assuming no DOMF.
!
FV2=graz_act(i,k,iphy)*Bio(i,k,iPhyF(iphy))
Bio_new(i,k,iPhyF(iphy))=Bio_new(i,k,iPhyF(iphy))- &
& FV2
Bio_new(i,k,iFecF(isfc))=Bio_new(i,k,iFecF(isfc))+ &
& FecPEL(iphy,isfc,ng)*FV2
Bio_new(i,k,iFecF(iffc))=Bio_new(i,k,iFecF(iffc))+ &
& FecPEL(iphy,iffc,ng)*FV2
Bio_new(i,k,iFeO_)=Bio_new(i,k,iFeO_)+ &
& (FecCYC(iphy,ng)+ &
& FecDOC(iphy,ng))*FV2
END IF
END DO
END DO
END DO
!
! Pigment Grazing. No fecal or dissolved terms for pigments.
!
DO ipig=1,Npig
DO iphy=1,Nphy
IF (iPigs(iphy,ipig).gt.0) THEN
itrc=iPigs(iphy,ipig)
DO k=1,N(ng)
DO i=Istr,Iend
IF (Bio(i,k,iPhyC(iphy)).gt.refuge(i,k,iphy)) THEN
FV1=graz_act(i,k,iphy)*Bio(i,k,itrc)
Bio_new(i,k,itrc)=Bio_new(i,k,itrc) - FV1
END IF
END DO
END DO
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
! Bacterial losses.
!-----------------------------------------------------------------------
!
! NOTE: Bacterial growth is completely reminerialized.
!
DO ibac=1,Nbac
DO k=1,N(ng)
DO i=Istr,Iend
!
! Grazing calculation. (All fecal material to slow sinking pool.)
!
!! WPB - There appears to be some rounding errors that cause bacteria
!! populations to drop just below initialization values. Once
!! they do, they never recover and the new lower values propagate
!! through the model. Only evident in Bac_P1 at the moment.
!!
!! FV1=BacCYC(ng)*Bio_new(i,k,iBacC(ibac))
!! FV2=BacPEL(ng)*Bio_new(i,k,iBacC(ibac))
!! FV3=BacDOC(ng)*Bio_new(i,k,iBacC(ibac))
!! FV4=FV1+FV2+FV3
!
! Carbon calculations.
!
Bio_new(i,k,iBacC(ibac))=Bio_new(i,k,iBacC(ibac))- &
!! & FV4
& Bio_new(i,k,iBacC(ibac))
Bio_new(i,k,iFecC(isfc))=Bio_new(i,k,iFecC(isfc))+ &
!! & FV2
& Bio_new(i,k,iBacC(ibac))* &
& BacPEL(ng)
Bio_new(i,k,iDOMC(ilab))=Bio_new(i,k,iDOMC(ilab))+ &
& (1.0_r8-cDOCfrac_c(ilab,ng))* &
!! & FV3
& Bio_new(i,k,iBacC(ibac))* &
& BacDOC(ng)
Bio_new(i,k,iCDMC(ilab))=Bio_new(i,k,iCDMC(ilab))+ &
& cDOCfrac_c(ilab,ng)* &
!! & FV3
& Bio_new(i,k,iBacC(ibac))* &
& BacDOC(ng)
Bio_new(i,k,iDIC_)=Bio_new(i,k,iDIC_)+ &
!! & FV1
& Bio_new(i,k,iBacC(ibac))* &
& BacCYC(ng)
!
! Nitrogen calculations.
!
Bio_new(i,k,iBacN(ibac))=Bio_new(i,k,iBacN(ibac))- &
!! & N2cBAC(ng)*FV4
& Bio_new(i,k,iBacN(ibac))
Bio_new(i,k,iFecN(isfc))=Bio_new(i,k,iFecN(isfc))+ &
!! & N2cBAC(ng)*FV2
& Bio_new(i,k,iBacN(ibac))* &
& BacPEL(ng)
Bio_new(i,k,iDOMN(ilab))=Bio_new(i,k,iDOMN(ilab))+ &
!! & N2cBAC(ng)*FV3
& Bio_new(i,k,iBacN(ibac))* &
& BacDOC(ng)
Bio_new(i,k,iNH4_)=Bio_new(i,k,iNH4_)+ &
!! & N2cBAC(ng)*FV1
& Bio_new(i,k,iBacN(ibac))* &
& BacCYC(ng)
!
! Phosphorous calculations.
!
Bio_new(i,k,iBacP(ibac))=Bio_new(i,k,iBacP(ibac))- &
!! & P2cBAC(ng)*FV4
& Bio_new(i,k,iBacP(ibac))
Bio_new(i,k,iFecP(isfc))=Bio_new(i,k,iFecP(isfc))+ &
!! & P2cBAC(ng)*FV2
& Bio_new(i,k,iBacP(ibac))* &
& BacPEL(ng)
Bio_new(i,k,iDOMP(ilab))=Bio_new(i,k,iDOMP(ilab))+ &
!! & P2cBAC(ng)*FV3
& Bio_new(i,k,iBacP(ibac))* &
& BacDOC(ng)
Bio_new(i,k,iPO4_)=Bio_new(i,k,iPO4_)+ &
!! & P2cBAC(ng)*FV1
& Bio_new(i,k,iBacP(ibac))* &
& BacCYC(ng)
!
! Iron calculations.
!
Bio_new(i,k,iBacF(ibac))=Bio_new(i,k,iBacF(ibac))- &
!! & Fe2cBAC(ng)*FV4
& Bio_new(i,k,iBacF(ibac))
Bio_new(i,k,iFecF(isfc))=Bio_new(i,k,iFecF(isfc))+ &
!! & Fe2cBAC(ng)*FV2
& Bio_new(i,k,iBacF(ibac))* &
& BacPEL(ng)
Bio_new(i,k,iFeO_)=Bio_new(i,k,iFeO_)+ &
!! & Fe2cBAC(ng)*(FV1+FV3)
& Bio_new(i,k,iBacF(ibac))* &
& (BacDOC(ng)+BacCYC(ng))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
! Fecal pellet remineralization.
!-----------------------------------------------------------------------
!
DO ifec=1,Nfec
DO k=1,N(ng)
DO i=Istr,Iend
!
! Carbon calculations. All carbon goes to CO2.
!
FV3=Regen_C(i,k,ifec)*Bio(i,k,iFecC(ifec))
Bio_new(i,k,iFecC(ifec))=Bio_new(i,k,iFecC(ifec))- &
& FV3
Bio_new(i,k,iDIC_)=Bio_new(i,k,iDIC_)+ &
& FV3
!
! Nitrogen calculations. Nitrogen goes to NH4.
!
FV2=Regen_N(i,k,ifec)*Bio(i,k,iFecN(ifec))
Bio_new(i,k,iFecN(ifec))=Bio_new(i,k,iFecN(ifec))- &
& FV2
Bio_new(i,k,iNH4_)=Bio_new(i,k,iNH4_)+ &
& FV2
!
! Silica calculations.
!
FV2=Regen_S(i,k,ifec)*Bio(i,k,iFecS(ifec))
Bio_new(i,k,iFecS(ifec))=Bio_new(i,k,iFecS(ifec))- &
& FV2
Bio_new(i,k,iSiO_)=Bio_new(i,k,iSiO_)+ &
& FV2
!
! Phosphorous calculations.
!
FV2=Regen_P(i,k,ifec)*Bio(i,k,iFecP(ifec))
Bio_new(i,k,iFecP(ifec))=Bio_new(i,k,iFecP(ifec))- &
& FV2
Bio_new(i,k,iPO4_)=Bio_new(i,k,iPO4_)+ &
& FV2
!
! Iron calculations.
!
FV2=Regen_F(i,k,ifec)*Bio(i,k,iFecF(ifec))
Bio_new(i,k,iFecF(ifec))=Bio_new(i,k,iFecF(ifec))- &
& FV2
Bio_new(i,k,iFeO_)=Bio_new(i,k,iFeO_)+ &
& FV2
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
! CDMC photolysis calculations.
!-----------------------------------------------------------------------
!
IF (RtUVR_flag(ng)) THEN
DO i=Istr,Iend
!
! If Ed_nz(i,N(ng)) > zero, then there is sunlight. Standardizing rate
! to 1500 umol quanta m-2 s-1.
!
IF (Ed_nz(i,N(ng)).ge.0.01) THEN
!
FV1=RtUVR_DIC(ng)*Ed_nz(i,N(ng))/1500.0_r8
FV2=RtUVR_DOC(ng)*Ed_nz(i,N(ng))/1500.0_r8
!
! FV4 equals the CDMC1 absorption at 410 nm. 0.012 converts to g m-3.
! FV5 equals the CDMC2 absorption at 410 nm.
! Weighted average attenuation of UVB of water at 300 nm = 0.2 m-1.
!
FV4=Bio(i,N(ng),iCDMC(ilab))*0.012_r8*aDOC410(ilab)
FV5=Bio(i,N(ng),iCDMC(irct))*0.012_r8*aDOC410(irct)
photo_decay=0.5_r8*Hz(i,j,N(ng))* &
& (0.2_r8+(FV4+FV5)*aDOC300(ilab))
FV3=EXP(-photo_decay)
photo_decay=2.0_r8*photo_decay
!
! Do not photolyze below the euphotic zone.
!
DO k=N(ng),Keuphotic(i),-1
IF (FV3.gt.0.01_r8) THEN
FV6=FV5+FV4
IF (FV6.gt.0.0_r8) THEN
FV7=FV4/FV6
photo_DIC=FV3*FV1*FV6
photo_DOC=FV3*FV2*FV6
total_photo=photo_DIC+photo_DOC
!
! NOTE: not testing for excess photolysis (CDOC going negative).
!
FV4=(1.0_r8-FV7)*total_photo
Bio_new(i,k,iCDMC(irct))=Bio_new(i,k,iCDMC(irct))-&
& FV4
Bio_new(i,k,iDOMC(ilab))=Bio_new(i,k,iDOMC(ilab))+&
& photo_DOC
Bio_new(i,k,iCDMC(ilab))=Bio_new(i,k,iCDMC(ilab))-&
& FV7*total_photo
Bio_new(i,k,iDIC_)=Bio_new(i,k,iDIC_)+ &
& photo_DIC
END IF
!
! FV4 equals the CDMC1 absorption at 410 nm. 0.012 converts to g m-3.
! FV5 equals the CDMC2 absorption at 410 nm.
! Weighted average attenuation of UVB of water at 300 nm = 0.2 m-1.
!
FV4=Bio(i,k,iCDMC(ilab))*0.012_r8*aDOC410(ilab)
FV5=Bio(i,k,iCDMC(irct))*0.012_r8*aDOC410(irct)
FV7=photo_decay+ &
& 0.5_r8*Hz(i,j,k)*(0.2_r8+(FV4+FV5)*aDOC300(ilab))
!
! If k is greater than the bottom of the euphotic zone (and by
! by extension the bottom boundary) or the decay constant is
! greater than 4.61 (or < 1% photolysis zone) then exit do loop.
!
FV3=EXP(-FV7)
!
! Store value for passage through entire Hz(i,j,k).
!
photo_decay=photo_decay+2.0_r8*FV7
END IF
END DO
END IF
END DO
END IF
!
!-----------------------------------------------------------------------
! Create optimal pigment ratios.
!-----------------------------------------------------------------------
!
DO i=Istr,Iend
IF (Keuphotic(i).le.N(ng)) THEN
DO iphy=1,Nphy
!
! Carbon to chlorophyll a ratio
! This statement says that nutrient limitation of C2CHL ratio overides
! light adaptation. Minimum of two functions may be more ecologically
! accurate?
!
DO k=N(ng),Keuphotic(i),-1
IF (b_C2Cn(iphy,ng).lt.0.0_r8+SMALL) THEN
C2CHL_w(k,iphy)=MIN((b_C2Cl(iphy,ng)+ &
& mxC2Cl(iphy,ng)*E0_nz(i,k)), &
& C2CHL_max(iphy,ng))
ELSE IF (C2nALG(i,k,iphy).gt. &
& minC2nALG(iphy,ng)+SMALL) THEN
C2CHL_w(k,iphy)=b_C2Cn(iphy,ng)+ &
& mxC2Cn(iphy,ng)* &
& (C2nALG(i,k,iphy)- &
& minC2nALG(iphy,ng))
ELSE
C2CHL_w(k,iphy)=MIN((b_C2Cl(iphy,ng)+ &
& mxC2Cl(iphy,ng)*E0_nz(i,k)), &
& C2CHL_max(iphy,ng))
END IF
END DO
!
! Chlorophyll a concentation per species. form g CHL a / g C
!
DO k=N(ng),Keuphotic(i),-1
Pigs_w(k,iphy,ichl)=1.0_r8/C2CHL_w(k,iphy)
END DO
!
! Chlorophyll b concentration per species. form g CHL b / g C
!
IF (iPigs(iphy,2).gt.0) THEN
DO k=N(ng),Keuphotic(i),-1
Pigs_w(k,iphy,2)=b_ChlB(iphy,ng)+ &
& mxChlB(iphy,ng)* &
& (C2CHL_w(k,iphy)- &
& b_C2Cl(iphy,ng))
Pigs_w(k,iphy,2)=Pigs_w(k,iphy,2)* &
& Pigs_w(k,iphy,ichl)
END DO
END IF
!
! Chlorophyll c concentration per species. form g CHL c / g C
!
IF (iPigs(iphy,3).gt.0) THEN
DO k=N(ng),Keuphotic(i),-1
Pigs_w(k,iphy,3)=b_ChlC(iphy,ng)+ &
& mxChlC(iphy,ng)* &
& (C2CHL_w(k,iphy)- &
& b_C2Cl(iphy,ng))
Pigs_w(k,iphy,3)=Pigs_w(k,iphy,3)* &
& Pigs_w(k,iphy,ichl)
END DO
END IF
!
! Photosynthetic caroteniods per species. form g PSC / g C
!
IF (iPigs(iphy,4).gt.0) THEN
DO k=N(ng),Keuphotic(i),-1
Pigs_w(k,iphy,4)=b_PSC(iphy,ng)+ &
& mxPSC(iphy,ng)* &
& (C2CHL_w(k,iphy)- &
& b_C2Cl(iphy,ng))
Pigs_w(k,iphy,4)=Pigs_w(k,iphy,4)* &
& Pigs_w(k,iphy,ichl)
END DO
END IF
!
! Photoprotective caroteniods per species. form g PPC / g C
!
IF (iPigs(iphy,5).gt.0) THEN
DO k=N(ng),Keuphotic(i),-1
Pigs_w(k,iphy,5)=b_PPC(iphy,ng)+ &
& mxPPC(iphy,ng)* &
& (C2CHL_w(k,iphy)- &
& b_C2Cl(iphy,ng))
Pigs_w(k,iphy,5)=Pigs_w(k,iphy,5)* &
& Pigs_w(k,iphy,ichl)
END DO
END IF
!
! Low Urobilin Phycoeurythin concentration per species. g LPUB / g C
!
IF (iPigs(iphy,6).gt.0) THEN
DO k=N(ng),Keuphotic(i),-1
Pigs_w(k,iphy,6)=b_LPUb(iphy,ng)+ &
& mxLPUb(iphy,ng)* &
& (C2CHL_w(k,iphy)- &
& b_C2Cl(iphy,ng))
Pigs_w(k,iphy,6)=Pigs_w(k,iphy,6)* &
& Pigs_w(k,iphy,ichl)
END DO
END IF
!
! High Urobilin Phycoeurythin concentration per species (g HPUB / g C).
!
IF (iPigs(iphy,7).gt.0) THEN
DO k=N(ng),Keuphotic(i),-1
Pigs_w(k,iphy,7)=b_HPUb(iphy,ng)+ &
& mxHPUb(iphy,ng)* &
& (C2CHL_w(k,iphy)- &
& b_C2Cl(iphy,ng))
Pigs_w(k,iphy,7)=Pigs_w(k,iphy,7)* &
& Pigs_w(k,iphy,ichl)
END DO
END IF
END DO
!
! Calculate pigment ratio changes.
! NOTE: 12 factor to convert to ugrams (mg m-3)
!
DO ipig=1,Npig
DO iphy=1,Nphy
IF (iPigs(iphy,ipig).gt.0) THEN
itrc=iPigs(iphy,ipig)
DO k=N(ng),Keuphotic(i),-1
IF ((Bio(i,k,iPhyC(iphy)).gt.0.0_r8).and. &
& (Bio(i,k,itrc).gt.0.0_r8)) THEN
FV1=Bio(i,k,iPhyC(iphy))*12.0_r8
FV2=GtALG_r(i,k,iphy)
FV3=FV1/ &
& (FV2/Pigs_w(k,iphy,ipig)+ &
& FV1*(1.0_r8-FV2)/ &
& Bio(i,k,itrc))
Bio_new(i,k,itrc)=Bio_new(i,k,itrc)+ &
& (FV3-Bio(i,k,itrc))
END IF
END DO
END IF
END DO
END DO
END IF
END DO
!
!-----------------------------------------------------------------------
! Vertical sinking terms.
!-----------------------------------------------------------------------
!
! Reconstruct vertical profile of selected biological constituents
! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
! grid box. Then, compute semi-Lagrangian flux due to sinking.
!
SINK_LOOP: DO isink=1,Nsink
itrc=idsink(isink)
!
! Copy concentration of biological particulates into scratch array
! "qc" (q-central, restrict it to be positive) which is hereafter
! interpreted as a set of grid-box averaged values for biogeochemical
! constituent concentration.
!
DO k=1,N(ng)
DO i=Istr,Iend
qc(i,k)=Bio(i,k,itrc)
END DO
END DO
!
DO k=N(ng)-1,1,-1
DO i=Istr,Iend
FC(i,k)=(qc(i,k+1)-qc(i,k))*Hz_inv2(i,k)
END DO
END DO
DO k=2,N(ng)-1
DO i=Istr,Iend
dltR=Hz(i,j,k)*FC(i,k)
dltL=Hz(i,j,k)*FC(i,k-1)
cff=Hz(i,j,k-1)+2.0_r8*Hz(i,j,k)+Hz(i,j,k+1)
cffR=cff*FC(i,k)
cffL=cff*FC(i,k-1)
!
! Apply PPM monotonicity constraint to prevent oscillations within the
! grid box.
!
IF ((dltR*dltL).le.0.0_r8) THEN
dltR=0.0_r8
dltL=0.0_r8
ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN
dltR=cffL
ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN
dltL=cffR
END IF
!
! Compute right and left side values (bR,bL) of parabolic segments
! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
!
! NOTE: Although each parabolic segment is monotonic within its grid
! box, monotonicity of the whole profile is not guaranteed,
! because bL(k+1)-bR(k) may still have different sign than
! qc(i,k+1)-qc(i,k). This possibility is excluded,
! after bL and bR are reconciled using WENO procedure.
!
cff=(dltR-dltL)*Hz_inv3(i,k)
dltR=dltR-cff*Hz(i,j,k+1)
dltL=dltL+cff*Hz(i,j,k-1)
bR(i,k)=qc(i,k)+dltR
bL(i,k)=qc(i,k)-dltL
WR(i,k)=(2.0_r8*dltR-dltL)**2
WL(i,k)=(dltR-2.0_r8*dltL)**2
END DO
END DO
cff=1.0E-14_r8
DO k=2,N(ng)-2
DO i=Istr,Iend
dltL=MAX(cff,WL(i,k ))
dltR=MAX(cff,WR(i,k+1))
bR(i,k)=(dltR*bR(i,k)+dltL*bL(i,k+1))/(dltR+dltL)
bL(i,k+1)=bR(i,k)
END DO
END DO
DO i=Istr,Iend
FC(i,N(ng))=0.0_r8 ! NO-flux boundary condition
#if defined LINEAR_CONTINUATION
bL(i,N(ng))=bR(i,N(ng)-1)
bR(i,N(ng))=2.0_r8*qc(i,N(ng))-bL(i,N(ng))
#elif defined NEUMANN
bL(i,N(ng))=bR(i,N(ng)-1)
bR(i,N(ng))=1.5_r8*qc(i,N(ng))-0.5_r8*bL(i,N(ng))
#else
bR(i,N(ng))=qc(i,N(ng)) ! default strictly monotonic
bL(i,N(ng))=qc(i,N(ng)) ! conditions
bR(i,N(ng)-1)=qc(i,N(ng))
#endif
#if defined LINEAR_CONTINUATION
bR(i,1)=bL(i,2)
bL(i,1)=2.0_r8*qc(i,1)-bR(i,1)
#elif defined NEUMANN
bR(i,1)=bL(i,2)
bL(i,1)=1.5_r8*qc(i,1)-0.5_r8*bR(i,1)
#else
bL(i,2)=qc(i,1) ! bottom grid boxes are
bR(i,1)=qc(i,1) ! re-assumed to be
bL(i,1)=qc(i,1) ! piecewise constant.
#endif
END DO
!
! Apply monotonicity constraint again, since the reconciled interfacial
! values may cause a non-monotonic behavior of the parabolic segments
! inside the grid box.
!
DO k=1,N(ng)
DO i=Istr,Iend
dltR=bR(i,k)-qc(i,k)
dltL=qc(i,k)-bL(i,k)
cffR=2.0_r8*dltR
cffL=2.0_r8*dltL
IF ((dltR*dltL).lt.0.0_r8) THEN
dltR=0.0_r8
dltL=0.0_r8
ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN
dltR=cffL
ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN
dltL=cffR
END IF
bR(i,k)=qc(i,k)+dltR
bL(i,k)=qc(i,k)-dltL
END DO
END DO
!
! After this moment reconstruction is considered complete. The next
! stage is to compute vertical advective fluxes, FC. It is expected
! that sinking may occurs relatively fast, the algorithm is designed
! to be free of CFL criterion, which is achieved by allowing
! integration bounds for semi-Lagrangian advective flux to use as
! many grid boxes in upstream direction as necessary.
!
! In the two code segments below, WL is the z-coordinate of the
! departure point for grid box interface z_w with the same indices;
! FC is the finite volume flux; ksource(:,k) is index of vertical
! grid box which contains the departure point (restricted by N(ng)).
! During the search: also add in content of whole grid boxes
! participating in FC.
!
cff=dtbio*ABS(Wbio(isink))
DO k=1,N(ng)
DO i=Istr,Iend
FC(i,k-1)=0.0_r8
WL(i,k)=z_w(i,j,k-1)+cff
WR(i,k)=Hz(i,j,k)*qc(i,k)
ksource(i,k)=k
END DO
END DO
DO k=1,N(ng)
DO ks=k,N(ng)-1
DO i=Istr,Iend
IF (WL(i,k).gt.z_w(i,j,ks)) THEN
ksource(i,k)=ks+1
FC(i,k-1)=FC(i,k-1)+WR(i,ks)
END IF
END DO
END DO
END DO
!
! Finalize computation of flux: add fractional part.
!
DO k=1,N(ng)
DO i=Istr,Iend
ks=ksource(i,k)
cu=MIN(1.0_r8,(WL(i,k)-z_w(i,j,ks-1))*Hz_inv(i,ks))
FC(i,k-1)=FC(i,k-1)+ &
& Hz(i,j,ks)*cu* &
& (bL(i,ks)+ &
& cu*(0.5_r8*(bR(i,ks)-bL(i,ks))- &
& (1.5_r8-cu)* &
& (bR(i,ks)+bL(i,ks)- &
& 2.0_r8*qc(i,ks))))
END DO
END DO
DO k=1,N(ng)
DO i=Istr,Iend
Bio(i,k,itrc)=qc(i,k)+(FC(i,k)-FC(i,k-1))*Hz_inv(i,k)
END DO
END DO
#ifdef BIO_SEDIMENT
!
! Particulate flux reaching the seafloor is remineralized and returned
! to the dissolved nitrate pool. Without this conversion, particulate
! material falls out of the system. This is a temporary fix to restore
! total nitrogen conservation. It will be replaced later by a
! parameterization that includes the time delay of remineralization
! and dissolved oxygen.
!
DO ifec=1,Nfec
IF (itrc.eq.iFecN(ifec)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iNO3_)=Bio(i,1,iNO3_)+cff1
END DO
ELSE IF (itrc.eq.iFecC(ifec)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iDIC_)=Bio(i,1,iDIC_)+cff1
END DO
ELSE IF (itrc.eq.iFecP(ifec)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iPO4_)=Bio(i,1,iPO4_)+cff1
END DO
ELSE IF (itrc.eq.iFecS(ifec)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iSiO_)=Bio(i,1,iSiO_)+cff1
END DO
ELSE IF (itrc.eq.iFecF(ifec)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iFeO_)=Bio(i,1,iFeO_)+cff1
END DO
END IF
END DO
DO iphy=1,Nphy
IF (itrc.eq.iPhyN(iphy)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iNO3_)=Bio(i,1,iNO3_)+cff1
END DO
ELSE IF (itrc.eq.iPhyC(iphy)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iDIC_)=Bio(i,1,iDIC_)+cff1
END DO
ELSE IF (itrc.eq.iPhyP(iphy)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iPO4_)=Bio(i,1,iPO4_)+cff1
END DO
ELSE IF (itrc.eq.iPhyS(iphy)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iSiO_)=Bio(i,1,iSiO_)+cff1
END DO
ELSE IF (itrc.eq.iPhyF(iphy)) THEN
DO i=Istr,Iend
cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iFeO_)=Bio(i,1,iFeO_)+cff1
END DO
END IF
END DO
#endif
END DO SINK_LOOP
!
!-----------------------------------------------------------------------
! Update the tendency arrays
!-----------------------------------------------------------------------
!
DO ibio=1,NBT
itrc=idbio(ibio)
DO k=1,N(ng)
DO i=Istr,Iend
Bio(i,k,itrc)=Bio(i,k,itrc)+dtbio*Bio_new(i,k,itrc)
END DO
END DO
END DO
END DO ITER_LOOP
!
!-----------------------------------------------------------------------
! Update global tracer variables: Add increment due to BGC processes
! to tracer array in time index "nnew". Index "nnew" is solution after
! advection and mixing and has transport units (m Tunits) hence the
! increment is multiplied by Hz. Notice that we need to subtract
! original values "Bio_old" at the top of the routine to just account
! for the concentractions affected by BGC processes. This also takes
! into account any constraints (non-negative concentrations, carbon
! concentration range) specified before entering BGC kernel. If "Bio"
! were unchanged by BGC processes, the increment would be exactly
! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
! bounded >=0 so that we can preserve total inventory of nutrients
! when advection causes tracer concentration to go negative.
!-----------------------------------------------------------------------
!
DO ibio=1,NBT
itrc=idbio(ibio)
DO k=1,N(ng)
DO i=Istr,Iend
cff=Bio(i,k,itrc)-Bio_old(i,k,itrc)
t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+cff*Hz(i,j,k)
END DO
END DO
END DO
END DO J_LOOP
!
RETURN
END SUBROUTINE ecosim_tile
END MODULE biology_mod