C $Header: /u/gcmpack/MITgcm/pkg/bling/bling_production.F,v 1.18 2017/12/12 19:12:58 mmazloff Exp $
C $Name: $
#include "BLING_OPTIONS.h"
CBOP
subroutine BLING_PROD(
I PTR_NO3, PTR_PO4, PTR_FE,
I PTR_O2, PTR_DON, PTR_DOP,
#ifdef ADVECT_PHYTO
I PTR_PHY,
#endif
O G_NO3, G_PO4, G_FE,
O G_O2, G_DON, G_DOP,
O G_CACO3, NCP,
I bi, bj, imin, imax, jmin, jmax,
I myIter, myTime, myThid )
C =================================================================
C | subroutine bling_prod
C | o Nutrient uptake and partitioning between organic pools.
C | - Phytoplankton specific growth rate is calculated
C | as a function of light, nutrient limitation, and
C | temperature.
C | - Population growth is calculated as a function of the local
C | phytoplankton biomass.
C =================================================================
implicit none
C === Global variables ===
C phyto_sm :: Small phytoplankton biomass
C phyto_lg :: Large phytoplankton biomass
C phyto_diaz :: Diazotroph phytoplankton biomass
C *** if ADVECT_PHYTO, these are fraction of total biomass instead ***
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "BLING_VARS.h"
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#ifdef ALLOW_AUTODIFF
# include "tamc.h"
#endif
C === Routine arguments ===
C bi,bj :: tile indices
C iMin,iMax :: computation domain: 1rst index range
C jMin,jMax :: computation domain: 2nd index range
C myTime :: current time
C myIter :: current timestep
C myThid :: thread Id. number
INTEGER bi, bj, imin, imax, jmin, jmax
_RL myTime
INTEGER myIter
INTEGER myThid
C === Input ===
C PTR_NO3 :: nitrate concentration
C PTR_PO4 :: phosphate concentration
C PTR_FE :: iron concentration
C PTR_DON :: dissolved organic nitrogen concentration
C PTR_DOP :: dissolved organic phosphorus concentration
C PTR_O2 :: oxygen concentration
C PTR_PHYTO :: total phytoplankton biomass
_RL PTR_NO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_DON(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_DOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#ifdef ADVECT_PHYTO
_RL PTR_PHY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#endif
C === Output ===
C G_xxx :: Tendency of xxx
_RL G_NO3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL G_PO4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL G_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL G_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL G_DON (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL G_DOP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL G_CACO3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#ifdef ALLOW_BLING
C === Local variables ===
C i,j,k :: loop indices
C Phy_lg_local :: biomass in large phytoplankton
C Phy_sm_local :: biomass in small phytoplankton
C Phy_diaz_local:: biomass in diazotroph phytoplankton
C NO3_lim :: nitrate limitation
C PO4_lim :: phosphate limitation
C Fe_lim :: iron limitation for phytoplankton
C Fe_lim_diaz :: iron limitation for diazotrophs
C alpha_Fe :: initial slope of the P-I curve
C theta_Fe :: Chl:C ratio
C theta_Fe_max :: Fe-replete maximum Chl:C ratio
C irrk :: nut-limited efficiency of algal photosystems
C irr_inst :: instantaneous light
C irr_eff :: effective irradiance
C mld :: mixed layer depth
C Pc_m :: light-saturated max photosynthesis rate for phyt
C Pc_m_diaz :: light-saturated max photosynthesis rate for diaz
C Pc_tot :: carbon-specific photosynthesis rate
C expkT :: temperature function
C mu :: net carbon-specific growth rate for phyt
C mu_diaz :: net carbon-specific growth rate for diaz
C PtoN :: variable ratio of phosphorus to nitrogen
C FetoN :: variable ratio of iron to nitrogen
C N_uptake :: NO3 utilization by phytoplankton
C N_fix :: Nitrogen fixation by diazotrophs
C N_den_pelag :: Pelagic denitrification
C N_den_benthic :: Benthic denitrification
C P_uptake :: PO4 utilization by phytoplankton
C Fe_uptake :: dissolved Fe utilization by phytoplankton
C CaCO3_uptake :: Calcium carbonate uptake for shell formation
C CaCO3_diss :: Calcium carbonate dissolution
C DON_prod :: production of dissolved organic nitrogen
C DON_remin :: remineralization of dissolved organic nitrogen
C DOP_prod :: production of dissolved organic phosphorus
C DOP_remin :: remineralization of dissolved organic phosphorus
C O2_prod :: production of oxygen
C frac_exp :: fraction of sinking particulate organic matter
C N_spm :: particulate sinking of nitrogen
C P_spm :: particulate sinking of phosphorus
C Fe_spm :: particulate sinking of iron
C N_dvm :: vertical transport of nitrogen by DVM
C P_dvm :: vertical transport of phosphorus by DVM
C Fe_dvm :: vertical transport of iron by DVM
C N_recycle :: recycling of newly-produced organic nitrogen
C P_recycle :: recycling of newly-produced organic phosphorus
C Fe_recycle :: recycling of newly-produced organic iron
C N_reminp :: remineralization of particulate organic nitrogen
C P_reminp :: remineralization of particulate organic nitrogen
C Fe_reminsum :: iron remineralization and adsorption
C N_remindvm :: nitrogen remineralization due to diel vertical migration
C P_remindvm :: phosphorus remineralization due to diel vertical migration
C Fe_remindvm :: iron remineralization due to diel vertical migration
C POC_flux :: particulate organic carbon flux
C NPP :: net primary production
C NCP :: net community production
C
INTEGER i,j,k
_RL Phy_lg_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Phy_sm_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Phy_diaz_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL NO3_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PO4_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_lim_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL light_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL expkT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Pc_m
_RL Pc_m_diaz
_RL theta_Fe_max
_RL theta_Fe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL thFe_inv(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL irrk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL irr_inst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL irr_eff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL mld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL mu(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL mu_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PtoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL FetoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_fix(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_den_pelag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_den_benthic(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL P_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL CaCO3_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL CaCO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL DON_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL DOP_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL DON_remin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL DOP_remin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL O2_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL frac_exp
_RL N_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL P_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL P_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL P_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL P_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_reminsum(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL Fe_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL POC_flux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL NPP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL NCP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#ifdef ML_MEAN_PHYTO
_RL tmp_p_sm_ML
_RL tmp_p_lg_ML
_RL tmp_p_diaz_ML
_RL tmp_ML
#endif
CEOP
c Initialize output and diagnostics
DO j=jmin,jmax
DO i=imin,imax
mld(i,j) = 0. _d 0
ENDDO
ENDDO
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
G_NO3(i,j,k) = 0. _d 0
G_PO4(i,j,k) = 0. _d 0
G_Fe(i,j,k) = 0. _d 0
G_O2(i,j,k) = 0. _d 0
G_DON(i,j,k) = 0. _d 0
G_DOP(i,j,k) = 0. _d 0
G_CaCO3(i,j,k) = 0. _d 0
N_uptake(i,j,k) = 0. _d 0
N_fix(i,j,k) = 0. _d 0
N_den_pelag(i,j,k) = 0. _d 0
N_den_benthic(i,j,k)= 0. _d 0
P_uptake(i,j,k) = 0. _d 0
Fe_uptake(i,j,k) = 0. _d 0
N_dvm(i,j,k) = 0. _d 0
P_dvm(i,j,k) = 0. _d 0
Fe_dvm(i,j,k) = 0. _d 0
CaCO3_uptake(i,j,k) = 0. _d 0
DON_prod(i,j,k) = 0. _d 0
DOP_prod(i,j,k) = 0. _d 0
O2_prod(i,j,k) = 0. _d 0
mu_diaz(i,j,k) = 0. _d 0
irr_eff(i,j,k) = 0. _d 0
irr_inst(i,j,k) = 0. _d 0
irrk(i,j,k) = 0. _d 0
NO3_lim(i,j,k) = 0. _d 0
PO4_lim(i,j,k) = 0. _d 0
Fe_lim(i,j,k) = 0. _d 0
Fe_lim_diaz(i,j,k) = 0. _d 0
light_lim(i,j,k) = 0. _d 0
PtoN(i,j,k) = 0. _d 0
FetoN(i,j,k) = 0. _d 0
NPP(i,j,k) = 0. _d 0
N_reminp(i,j,k) = 0. _d 0
P_reminp(i,j,k) = 0. _d 0
Fe_reminsum(i,j,k) = 0. _d 0
N_remindvm(i,j,k) = 0. _d 0
P_remindvm(i,j,k) = 0. _d 0
ENDDO
ENDDO
ENDDO
c-----------------------------------------------------------
c avoid negative nutrient concentrations that can result from
c advection when low concentrations
#ifdef BLING_NO_NEG
CALL BLING_MIN_VAL( PTR_NO3, 1. _d -7 ,bi ,bj)
CALL BLING_MIN_VAL( PTR_PO4, 1. _d -8 ,bi ,bj)
CALL BLING_MIN_VAL( PTR_FE, 1. _d -11 ,bi ,bj)
CALL BLING_MIN_VAL( PTR_O2, 1. _d -11 ,bi ,bj)
CALL BLING_MIN_VAL( PTR_DON, 1. _d -11 ,bi ,bj)
CALL BLING_MIN_VAL( PTR_DOP, 1. _d -11 ,bi ,bj)
#endif
c-----------------------------------------------------------
c Phytoplankton size classes
c
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
#ifdef ADVECT_PHYTO
Phy_lg_local(i,j,k) = PTR_PHY(i,j,k)*phyto_lg(i,j,k,bi,bj)
Phy_sm_local(i,j,k) = PTR_PHY(i,j,k)*phyto_sm(i,j,k,bi,bj)
Phy_diaz_local(i,j,k) = PTR_PHY(i,j,k)*phyto_diaz(i,j,k,bi,bj)
#else
Phy_lg_local(i,j,k) = phyto_lg(i,j,k,bi,bj)
Phy_sm_local(i,j,k) = phyto_sm(i,j,k,bi,bj)
Phy_diaz_local(i,j,k) = phyto_diaz(i,j,k,bi,bj)
#endif
ENDDO
ENDDO
ENDDO
c-----------------------------------------------------------
c mixed layer depth calculation for light, phytoplankton and dvm
c do not need to calculate if not using ML_MEAN_LIGHT, ML_MEAN_PHYTO,
c and USE_BLING_DVM
c (with BLING_ADJOINT_SAFE flag, USE_BLING_DVM is undefined)
#if ( defined (ML_MEAN_LIGHT)
defined (ML_MEAN_PHYTO)
defined (USE_BLING_DVM) )
CALL BLING_MIXEDLAYER(
U mld,
I bi, bj, imin, imax, jmin, jmax,
I myIter, myTime, myThid)
#endif
c Phytoplankton mixing
c The mixed layer is assumed to homogenize vertical gradients of phytoplankton.
c This allows for basic Sverdrup dynamics in a qualitative sense.
c This has not been thoroughly tested, and care should be
c taken with its interpretation.
#ifdef ML_MEAN_PHYTO
DO j=jmin,jmax
DO i=imin,imax
tmp_p_sm_ML = 0. _d 0
tmp_p_lg_ML = 0. _d 0
tmp_p_diaz_ML = 0. _d 0
tmp_ML = 0. _d 0
DO k=1,Nr
IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
IF ((-rf(k+1) .le. mld(i,j)).and.
& (-rf(k+1).lt.200. _d 0)) THEN
tmp_p_sm_ML = tmp_p_sm_ML+Phy_sm_local(i,j,k)*drF(k)
& *hFacC(i,j,k,bi,bj)
tmp_p_lg_ML = tmp_p_lg_ML+Phy_lg_local(i,j,k)*drF(k)
& *hFacC(i,j,k,bi,bj)
tmp_p_diaz_ML = tmp_p_diaz_ML+Phy_diaz_local(i,j,k)*drF(k)
& *hFacC(i,j,k,bi,bj)
tmp_ML = tmp_ML + drF(k)
ENDIF
ENDIF
ENDDO
DO k=1,Nr
IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
IF ((-rf(k+1) .le. mld(i,j)).and.
& (-rf(k+1).lt.200. _d 0)) THEN
Phy_lg_local(i,j,k) = max(1. _d -8,tmp_p_lg_ML/tmp_ML)
Phy_sm_local(i,j,k) = max(1. _d -8,tmp_p_sm_ML/tmp_ML)
Phy_diaz_local(i,j,k) = max(1. _d -8,tmp_p_diaz_ML/tmp_ML)
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
#endif
c-----------------------------------------------------------
c light availability for biological production
CALL BLING_LIGHT(
I mld,
U irr_inst, irr_eff,
I bi, bj, imin, imax, jmin, jmax,
I myIter, myTime, myThid )
c phytoplankton photoadaptation to local light level
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
irr_mem(i,j,k,bi,bj) = irr_mem(i,j,k,bi,bj) +
& (irr_eff(i,j,k) - irr_mem(i,j,k,bi,bj))*
& min( 1. _d 0, gamma_irr_mem*PTRACERS_dTLev(k) )
ENDDO
ENDDO
ENDDO
c ---------------------------------------------------------------------
c Nutrient uptake and partitioning between organic pools
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
c ---------------------------------------------------------------------
c First, calculate the limitation terms for NUT and Fe, and the
c Fe-limited Chl:C maximum. The light-saturated maximal photosynthesis
c rate term (Pc_m) is simply the product of a prescribed maximal
c photosynthesis rate (Pc_0/Pc_2d), the Eppley temperature dependence, and a
c resource limitation term. The iron limitation term has a lower limit
c of Fe_lim_min and is scaled by (k_Fe2P + Fe2P_max) / Fe2P_max so that
c it approaches 1 as Fe approaches infinity. Thus, it is of comparable
c magnitude to the macro-nutrient limitation term.
c Macro-nutrient limitation
NO3_lim(i,j,k) = PTR_NO3(i,j,k)/(PTR_NO3(i,j,k)+k_NO3)
PO4_lim(i,j,k) = PTR_PO4(i,j,k)/(PTR_PO4(i,j,k)+k_PO4)
c Iron limitation
Fe_lim(i,j,k) = PTR_FE(i,j,k)
& / (PTR_FE(i,j,k)+k_Fe2d(i,j,bi,bj))
Fe_lim_diaz(i,j,k) = PTR_FE(i,j,k)
& / (PTR_FE(i,j,k)+k_Fe_diaz2d(i,j,bi,bj))
c ---------------------------------------------------------------------
c Diazotrophs are assumed to be more strongly temperature sensitive,
c given their observed restriction to relatively warm waters. Presumably
c this is because of the difficulty of achieving N2 fixation in an oxic
c environment. Thus, they have lower pc_0 and higher kappa_eppley.
c Taking the square root, to provide the geometric mean.
expkT(i,j,k) = exp(kappa_eppley * theta(i,j,k,bi,bj))
c Light-saturated maximal photosynthesis rate
Pc_m = Pc_2d(i,j,bi,bj) * expkT(i,j,k)
& * min(NO3_lim(i,j,k), PO4_lim(i,j,k), Fe_lim(i,j,k))
& * maskC(i,j,k,bi,bj)
Pc_m_diaz = Pc_2d_diaz(i,j,bi,bj)
& * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
& * min(PO4_lim(i,j,k), Fe_lim_diaz(i,j,k))
& * maskC(i,j,k,bi,bj)
c (Pc_m and Pc_m_diaz crash adjoint if get too small)
#ifdef BLING_ADJOINT_SAFE
Pc_m = max(Pc_m ,maskC(i,j,k,bi,bj)*1. _d -15)
Pc_m_diaz = max(Pc_m_diaz,maskC(i,j,k,bi,bj)*1. _d -15)
#endif
c ---------------------------------------------------------------------
c Fe limitation 1) reduces photosynthetic efficiency (alpha_Fe)
c and 2) reduces the maximum achievable Chl:C ratio (theta_Fe)
c below a prescribed, Fe-replete maximum value (theta_Fe_max),
c to approach a prescribed minimum Chl:C (theta_Fe_min) under extreme
c Fe-limitation.
theta_Fe_max = theta_Fe_max_lo+
& (theta_Fe_max_hi-theta_Fe_max_lo)*Fe_lim(i,j,k)
theta_Fe(i,j,k) = theta_Fe_max
& / (1. _d 0 + alpha_photo2d(i,j,bi,bj)
& *theta_Fe_max
& *irr_mem(i,j,k,bi,bj)/(epsln + 2. _d 0*Pc_m))
c for diagnostics: C:Chl ratio in g C / g Chl
IF ( theta_Fe(i,j,k) .EQ.0. ) THEN
thFe_inv(i,j,k) = 0.
ELSE
thFe_inv(i,j,k) = 1./theta_Fe(i,j,k)
ENDIF
c ---------------------------------------------------------------------
c Nutrient-limited efficiency of algal photosystems, irrk, is calculated
c with the iron limitation term included as a multiplier of the
c theta_Fe_max to represent the importance of Fe in forming chlorophyll
c accessory antennae, which do not affect the Chl:C but still affect the
c phytoplankton ability to use light (eg Stzrepek & Harrison, Nature 2004).
irrk(i,j,k) = Pc_m
& /(epsln + alpha_photo2d(i,j,bi,bj)*theta_Fe_max)
& + irr_mem(i,j,k,bi,bj)/2. _d 0
c Carbon-specific photosynthesis rate
mu(i,j,k) = Pc_m * ( 1. _d 0 - exp(-irr_eff(i,j,k)
& /(epsln + irrk(i,j,k))))
mu_diaz(i,j,k) = Pc_m_diaz * ( 1. _d 0 - exp(-irr_eff(i,j,k)
& /(epsln + irrk(i,j,k))))
light_lim(i,j,k) = ( 1. _d 0 - exp(-irr_eff(i,j,k)
& /(epsln + irrk(i,j,k))))
c Temperature threshold for diazotrophs
IF (theta(i,j,k,bi,bj) .lt. 14) THEN
mu_diaz(i,j,k) = 0. _d 0
ENDIF
c Stoichiometry
PtoN(i,j,k) = PtoN_min + (PtoN_max - PtoN_min) *
& PTR_PO4(i,j,k) / (k_PtoN + PTR_PO4(i,j,k))
FetoN(i,j,k) = FetoN_min + (FetoN_max - FetoN_min) *
& PTR_FE(i,j,k) / (k_FetoN + PTR_FE(i,j,k))
c Nutrient uptake
N_uptake(i,j,k) = mu(i,j,k)*(Phy_sm_local(i,j,k)
& + Phy_lg_local(i,j,k))
N_fix(i,j,k) = mu_diaz(i,j,k) * Phy_diaz_local(i,j,k)
P_uptake(i,j,k) = (N_uptake(i,j,k) +
& N_fix(i,j,k)) * PtoN(i,j,k)
Fe_uptake(i,j,k) = (N_uptake(i,j,k) +
& N_fix(i,j,k)) * FetoN(i,j,k)
c ---------------------------------------------------------------------
c Alkalinity is consumed through the production of CaCO3. Here, thi
c is simply a linear function of the implied growth rate of small
c phytoplankton, which gave a reasonably good fit to the global
c observational synthesis of Dunne (2009). This is consistent
c with the findings of Jin et al. (GBC,2006).
CaCO3_uptake(i,j,k) = Phy_sm_local(i,j,k)*phi_sm2d(i,j,bi,bj)
& * mu(i,j,k) * CatoN
ENDIF
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
Phy_lg_local(i,j,k) = Phy_lg_local(i,j,k) +
& Phy_lg_local(i,j,k)*(mu(i,j,k) - lambda_0
& *expkT(i,j,k) *
& (Phy_lg_local(i,j,k)/pivotal)**(1. / 3.))
& * PTRACERS_dTLev(k)
Phy_sm_local(i,j,k) = Phy_sm_local(i,j,k) +
& Phy_sm_local(i,j,k)*(mu(i,j,k) - lambda_0
& *expkT(i,j,k) * (Phy_sm_local(i,j,k)/pivotal) )
& * PTRACERS_dTLev(k)
Phy_diaz_local(i,j,k) = Phy_diaz_local(i,j,k) +
& Phy_diaz_local(i,j,k)*(mu_diaz(i,j,k) - lambda_0
& *expkT(i,j,k) * (Phy_diaz_local(i,j,k)/pivotal) )
& * PTRACERS_dTLev(k)
ENDIF
ENDDO
ENDDO
ENDDO
c Separate loop for adjoint stores
CADJ STORE Phy_sm_local = comlev1, key = ikey_dynamics, kind=isbyte
CADJ STORE Phy_lg_local = comlev1, key = ikey_dynamics, kind=isbyte
CADJ STORE Phy_diaz_local = comlev1, key = ikey_dynamics, kind=isbyte
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
#ifdef ADVECT_PHYTO
c update tracer here
PTR_PHY(i,j,k) = Phy_lg_local(i,j,k) + Phy_sm_local(i,j,k)
& + Phy_diaz_local(i,j,k)
c update fractional abundance
phyto_lg(i,j,k,bi,bj) = Phy_lg_local(i,j,k)/PTR_PHY(i,j,k)
phyto_sm(i,j,k,bi,bj) = Phy_sm_local(i,j,k)/PTR_PHY(i,j,k)
phyto_diaz(i,j,k,bi,bj) = Phy_diaz_local(i,j,k)/PTR_PHY(i,j,k)
#else
c update biomass
phyto_lg(i,j,k,bi,bj) = Phy_lg_local(i,j,k)
phyto_sm(i,j,k,bi,bj) = Phy_sm_local(i,j,k)
phyto_diaz(i,j,k,bi,bj) = Phy_diaz_local(i,j,k)
#endif
c use the diagnostic biomass to calculate the chl concentration
c in mg/m3 (carbon = 12.01 g/mol)
chl(i,j,k,bi,bj) = max(chl_min, CtoN * 12.01 * 1. _d 3 *
& theta_Fe(i,j,k) *
& (Phy_lg_local(i,j,k) + Phy_sm_local(i,j,k)
& + Phy_diaz_local(i,j,k)))
ENDIF
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
c ---------------------------------------------------------------------
c Partitioning between organic pools
c The uptake of nutrients is assumed to contribute to the growth of
c phytoplankton, which subsequently die and are consumed by heterotrophs.
c This can involve the transfer of nutrient elements between many
c organic pools, both particulate and dissolved, with complex histories.
c We take a simple approach here, partitioning the total uptake into two
c fractions - sinking and non-sinking - as a function of temperature,
c following Dunne et al. (2005).
c Then, the non-sinking fraction is further subdivided, such that the
c majority is recycled instantaneously to the inorganic nutrient pool,
c representing the fast turnover of labile dissolved organic matter via
c the microbial loop, and the remainder is converted to semi-labile
c dissolved organic matter. Iron and macro-nutrient are treated
c identically for the first step, but all iron is recycled
c instantaneously in the second step (i.e. there is no dissolved organic
c iron pool).
c sinking fraction: particulate organic matter
frac_exp = (phi_sm2d(i,j,bi,bj) + phi_lg2d(i,j,bi,bj) *
& (mu(i,j,k)/(epsln + lambda_0*expkT(i,j,k)))**2.)/
& (1. + (mu(i,j,k)/(epsln + lambda_0*expkT(i,j,k)))**2.)*
& exp(kappa_remin * theta(i,j,k,bi,bj))
#ifdef USE_BLING_DVM
N_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
& (N_uptake(i,j,k) + N_fix(i,j,k))
P_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
& P_uptake(i,j,k)
Fe_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
& Fe_uptake(i,j,k)
N_dvm(i,j,k) = frac_exp *
& (N_uptake(i,j,k) + N_fix(i,j,k)) - N_spm(i,j,k)
P_dvm(i,j,k) = frac_exp * P_uptake(i,j,k) -
& P_spm(i,j,k)
Fe_dvm(i,j,k) = frac_exp * Fe_uptake(i,j,k) -
& Fe_spm(i,j,k)
#else
N_spm(i,j,k) = frac_exp * (N_uptake(i,j,k) + N_fix(i,j,k))
P_spm(i,j,k) = frac_exp *P_uptake(i,j,k)
Fe_spm(i,j,k) = frac_exp * Fe_uptake(i,j,k)
N_dvm(i,j,k) = 0
P_dvm(i,j,k) = 0
Fe_dvm(i,j,k) = 0
#endif
c the remainder is divided between instantaneously recycled and
c long-lived dissolved organic matter.
DON_prod(i,j,k) = phi_DOM2d(i,j,bi,bj)*(N_uptake(i,j,k)
& + N_fix(i,j,k) - N_spm(i,j,k)
& - N_dvm(i,j,k))
DOP_prod(i,j,k) = phi_DOM2d(i,j,bi,bj)*(P_uptake(i,j,k)
& - P_spm(i,j,k) - P_dvm(i,j,k))
N_recycle(i,j,k) = N_uptake(i,j,k) + N_fix(i,j,k)
& - N_spm(i,j,k) - DON_prod(i,j,k)
& - N_dvm(i,j,k)
P_recycle(i,j,k) = P_uptake(i,j,k)
& - P_spm(i,j,k) - DOP_prod(i,j,k)
& - P_dvm(i,j,k)
Fe_recycle(i,j,k) = Fe_uptake(i,j,k)
& - Fe_spm(i,j,k) - Fe_dvm(i,j,k)
ENDIF
ENDDO
ENDDO
ENDDO
c-----------------------------------------------------------
c remineralization of sinking organic matter
CALL BLING_REMIN(
I PTR_NO3, PTR_FE, PTR_O2, irr_inst,
I N_spm, P_spm, Fe_spm, CaCO3_uptake,
U N_reminp, P_reminp, Fe_reminsum,
U N_den_benthic, CACO3_diss,
I bi, bj, imin, imax, jmin, jmax,
I myIter, myTime, myThid)
c-----------------------------------------------------------
c remineralization from diel vertical migration
#ifdef USE_BLING_DVM
CALL BLING_DVM(
I N_dvm,P_dvm,Fe_dvm,
I PTR_O2, mld,
O N_remindvm, P_remindvm, Fe_remindvm,
I bi, bj, imin, imax, jmin, jmax,
I myIter, myTime, myThid)
#endif
c-----------------------------------------------------------
c
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
c Dissolved organic matter slow remineralization
#ifdef BLING_NO_NEG
DON_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DON
& *PTR_DON(i,j,k),0. _d 0)
DOP_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DOP
& *PTR_DOP(i,j,k),0. _d 0)
#else
DON_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DON
& *PTR_DON(i,j,k)
DOP_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DOP
& *PTR_DOP(i,j,k)
#endif
c Pelagic denitrification
c If anoxic
IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
IF (PTR_NO3(i,j,k) .gt. oxic_min) THEN
N_den_pelag(i,j,k) = max(epsln, (NO3toN *
& ((1. _d 0 - phi_DOM2d(i,j,bi,bj))
& * (N_reminp(i,j,k)
& + N_remindvm(i,j,k)) + DON_remin(i,j,k) +
& N_recycle(i,j,k))) - N_den_benthic(i,j,k))
ENDIF
ENDIF
c oxygen production through photosynthesis
O2_prod(i,j,k) = O2toN * N_uptake(i,j,k)
& + (O2toN - 1.25 _d 0) * N_fix(i,j,k)
c-----------------------------------------------------------
C ADD TERMS
c Nutrients
c Sum of fast recycling, decay of sinking POM, and decay of DOM,
c less uptake, (less denitrification).
G_PO4(i,j,k) = -P_uptake(i,j,k) + P_recycle(i,j,k)
& + (1. _d 0 - phi_DOM2d(i,j,bi,bj))
& * (P_reminp(i,j,k)
& + P_remindvm(i,j,k)) + DOP_remin(i,j,k)
G_NO3(i,j,k) = -N_uptake(i,j,k)
IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
c Anoxic
G_NO3(i,j,k) = G_NO3(i,j,k)
& - N_den_pelag(i,j,k) - N_den_benthic(i,j,k)
ELSE
c Oxic
G_NO3(i,j,k) = G_NO3(i,j,k)
& + N_recycle(i,j,k)
& + (1. _d 0 - phi_DOM2d(i,j,bi,bj))
& * (N_reminp(i,j,k) + N_remindvm(i,j,k))
& + DON_remin(i,j,k)
ENDIF
c Iron
c remineralization, sediments and adsorption are all bundled into
c Fe_reminsum
G_FE(i,j,k) = -Fe_uptake(i,j,k) + Fe_reminsum(i,j,k)
& + Fe_remindvm(i,j,k) + Fe_recycle(i,j,k)
c Dissolved Organic Matter
c A fraction of POM remineralization goes into dissolved pools.
G_DON(i,j,k) = DON_prod(i,j,k) + phi_DOM2d(i,j,bi,bj)
& * (N_reminp(i,j,k) + N_remindvm(i,j,k))
& - DON_remin(i,j,k)
G_DOP(i,j,k) = DOP_prod(i,j,k) + phi_DOM2d(i,j,bi,bj)
& * (P_reminp(i,j,k) + P_remindvm(i,j,k))
& - DOP_remin(i,j,k)
c Oxygen:
c Assuming constant O2:N ratio in terms of oxidant required per mol of
c organic N. This implies a constant stoichiometry of C:N and H:N (where H is
c reduced, organic H). Because the N provided by N2 fixation is reduced from
c N2, rather than NO3-, the o2_2_n_fix is slightly less than the NO3- based
c ratio (by 1.25 mol O2/ mol N). Account for the organic matter respired
c through benthic denitrification by subtracting 5/4 times the benthic
c denitrification NO3 utilization rate from the overall oxygen consumption.
G_O2(i,j,k) = O2_prod(i,j,k)
c If oxic
IF (PTR_O2(i,j,k) .gt. oxic_min) THEN
G_O2(i,j,k) = G_O2(i,j,k)
& -O2toN * ((1. _d 0 - phi_DOM2d(i,j,bi,bj))
& * (N_reminp(i,j,k) + N_remindvm(i,j,k))
& + DON_remin(i,j,k) + N_recycle(i,j,k))
c If anoxic but NO3 concentration is very low
c (generate negative O2; proxy for HS-).
ELSEIF (PTR_NO3(i,j,k) .lt. oxic_min) THEN
G_O2(i,j,k) = G_O2(i,j,k)
& -O2toN * ((1. _d 0 - phi_DOM2d(i,j,bi,bj))
& * (N_reminp(i,j,k) + N_remindvm(i,j,k))
& + DON_remin(i,j,k) + N_recycle(i,j,k))
& + N_den_benthic(i,j,k) * 1.25 _d 0
ENDIF
G_CaCO3(i,j,k) = CaCO3_diss(i,j,k) - CaCO3_uptake(i,j,k)
c Carbon flux diagnostic
POC_flux(i,j,k) = CtoN * N_spm(i,j,k)
NPP(i,j,k) = (N_uptake(i,j,k) + N_fix(i,j,k)) * CtoN
NCP(i,j,k) = (N_uptake(i,j,k) + N_fix(i,j,k)
& - N_recycle(i,j,k)
& - (1. _d 0 - phi_DOM2d(i,j,bi,bj))
& * (N_reminp(i,j,k) + N_remindvm(i,j,k))
& - DON_remin(i,j,k) ) * CtoN
c for diagnostics: convert to mol C/m3
Phy_lg_local(i,j,k) = Phy_lg_local(i,j,k) * CtoN
Phy_sm_local(i,j,k) = Phy_sm_local(i,j,k) * CtoN
Phy_diaz_local(i,j,k) = Phy_diaz_local(i,j,k) * CtoN
c for constraints determine POC, assuming that phytoplankton carbon
c is 30% of POC
poc(i,j,k,bi,bj) = (Phy_lg_local(i,j,k) + Phy_sm_local(i,j,k) +
& Phy_diaz_local(i,j,k) ) * 3.33333 _d 0
ENDIF
ENDDO
ENDDO
ENDDO
c ---------------------------------------------------------------------
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
c 3d global variables
CALL DIAGNOSTICS_FILL(Phy_sm_local,'BLGPSM ',0,Nr,1,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(Phy_lg_local,'BLGPLG ',0,Nr,1,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(Phy_diaz_local,'BLGPDIA ',0,Nr,1,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(chl, 'BLGCHL ',0,Nr,1,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(irr_mem,'BLGIMEM ',0,Nr,1,bi,bj,myThid)
c 3d local variables
CALL DIAGNOSTICS_FILL(irrk, 'BLGIRRK ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(irr_eff, 'BLGIEFF ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(poc, 'BLGPOC ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(theta_Fe,'BLGCHL2C',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(thFe_inv,'BLGC2CHL',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(Fe_lim, 'BLGFELIM',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(NO3_lim, 'BLGNLIM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(PO4_lim, 'BLGPLIM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(light_lim,'BLGLLIM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(NPP, 'BLGNPP ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(NCP, 'BLGNCP ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(Fe_dvm, 'BLGFEDVM',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(Fe_spm, 'BLGFESPM',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(Fe_recycle, 'BLGFEREC',0,Nr,2,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(Fe_remindvm, 'BLGFERD ',0,Nr,2,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(Fe_uptake,'BLGFEUP ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(N_den_benthic,'BLGNDENB',0,Nr,2,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(N_den_pelag, 'BLGNDENP',0,Nr,2,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(N_dvm, 'BLGNDVM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(N_fix, 'BLGNFIX ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(DON_prod, 'BLGDONP ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(DON_remin,'BLGDONR ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(N_spm, 'BLGNSPM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(N_recycle,'BLGNREC ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(N_remindvm,'BLGNRD ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(N_reminp, 'BLGNREM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(N_uptake, 'BLGNUP ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(P_dvm, 'BLGPDVM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(DOP_prod, 'BLGDOPP ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(P_spm, 'BLGPSPM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(P_recycle,'BLGPREC ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(P_remindvm,'BLGPRD ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(P_reminp, 'BLGPREM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(P_uptake, 'BLGPUP ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(mu, 'BLGMU ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(mu_diaz, 'BLGMUDIA',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(CaCO3_diss, 'BLGCCdis',0,Nr,2,bi,bj,
& myThid)
CALL DIAGNOSTICS_FILL(CaCO3_uptake,'BLGCCpro',0,Nr,2,bi,bj,
& myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* ALLOW_BLING */
RETURN
END