pkg/bling/bling_bio_nitrogen.F

e0f9a7ba0b Matt*0001 #include "BLING_OPTIONS.h"
a284455135 Matt*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
e0f9a7ba0b Matt*0005
                0006 CBOP
a284455135 Matt*0007       SUBROUTINE BLING_BIO_NITROGEN(
e0f9a7ba0b Matt*0008      I           PTR_O2, PTR_FE, PTR_PO4, PTR_DOP,
                0009      I           PTR_NO3, PTR_DON,
                0010 #ifdef USE_SIBLING
                0011      I           PTR_SI,
                0012 #endif
                0013 #ifdef ADVECT_PHYTO
                0014      I           PTR_PHY,
                0015 #endif
                0016      O           G_DIC, G_ALK, G_O2, G_FE,
                0017      O           G_PO4, G_DOP, G_NO3, G_DON,
                0018 #ifdef USE_SIBLING
                0019      O           G_SI,
                0020 #endif
                0021      I           bi, bj, imin, imax, jmin, jmax,
                0022      I           myTime, myIter, myThid)
                0023
                0024 C     =================================================================
                0025 C     | subroutine bling_bio_nitrogen
                0026 C     | o Nutrient uptake and partitioning between organic pools.
                0027 C     | - Phytoplankton specific growth rate is calculated
                0028 C     |   as a function of light, nutrient limitation, and
                0029 C     |   temperature.
                0030 C     | - Population growth is calculated as a function of the local
                0031 C     |   phytoplankton biomass.
                0032 C     | o Organic matter export and remineralization.
                0033 C     | - Sinking particulate flux and diel migration contribute to
                0034 C     |   export.
                0035 C     | - Benthic denitrification
                0036 C     | - Iron source from sediments
                0037 C     | - Iron scavenging
                0038 C     | o Diel Vertical Migration
                0039 C     =================================================================
                0040
                0041       IMPLICIT NONE
                0042
                0043 C     === Global variables ===
                0044 C     phyto_sm     :: small phytoplankton biomass
                0045 C     phyto_lg     :: large phytoplankton biomass
                0046 C     phyto_diaz   :: diazotroph phytoplankton biomass
                0047 C *** if ADVECT_PHYTO, these are fraction of total biomass instead ***
                0048
                0049 #include "SIZE.h"
                0050 #include "DYNVARS.h"
                0051 #include "EEPARAMS.h"
                0052 #include "PARAMS.h"
                0053 #include "GRID.h"
                0054 #include "BLING_VARS.h"
                0055 #include "PTRACERS_SIZE.h"
                0056 #include "PTRACERS_PARAMS.h"
a284455135 Matt*0057 #ifdef ALLOW_AUTODIFF_TAMC
e0f9a7ba0b Matt*0058 # include "tamc.h"
                0059 #endif
                0060
                0061 C     === Routine arguments ===
                0062 C     bi,bj         :: tile indices
                0063 C     iMin,iMax     :: computation domain: 1rst index range
                0064 C     jMin,jMax     :: computation domain: 2nd  index range
                0065 C     myTime        :: current time
                0066 C     myIter        :: current timestep
                0067 C     myThid        :: thread Id. number
                0068       INTEGER bi, bj, imin, imax, jmin, jmax
                0069       _RL     myTime
                0070       INTEGER myIter
                0071       INTEGER myThid
                0072 C     === Input ===
                0073 C     PTR_O2        :: oxygen concentration
                0074 C     PTR_FE        :: iron concentration
                0075 C     PTR_PO4       :: phosphate concentration
                0076 C     PTR_DOP       :: dissolved organic phosphorus concentration
                0077 C     PTR_NO3       :: nitrate concentration
                0078 C     PTR_DON       :: dissolved organic nitrogen concentration
                0079 C     PTR_SI        :: silica concentration
                0080 C     PTR_PHY       :: total phytoplankton biomass
                0081       _RL     PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0082       _RL     PTR_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0083       _RL     PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0084       _RL     PTR_DOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0085       _RL     PTR_NO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0086       _RL     PTR_DON(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0087 #ifdef USE_SIBLING
                0088       _RL     PTR_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0089 #endif
                0090 #ifdef ADVECT_PHYTO
                0091       _RL     PTR_PHY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0092 #endif
                0093
                0094 C     === Output ===
                0095 C     G_xxx         :: tendency term for tracer xxx
                0096       _RL     G_DIC     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0097       _RL     G_ALK     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0098       _RL     G_O2      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0099       _RL     G_FE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0100       _RL     G_PO4     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0101       _RL     G_DOP     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0102       _RL     G_NO3     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0103       _RL     G_DON     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0104 #ifdef USE_SIBLING
                0105       _RL     G_SI      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0106 #endif
                0107
                0108 #ifdef ALLOW_BLING
                0109 # ifndef USE_BLING_V1
                0110 C     === Local variables ===
                0111 C     i,j,k         :: loop indices
                0112 C     Phy_lg_local  :: biomass in large phytoplankton
                0113 C     Phy_sm_local  :: biomass in small phytoplankton
                0114 C     Phy_diaz_local:: biomass in diazotroph phytoplankton
                0115 C     NO3_lim       :: nitrate limitation
                0116 C     PO4_lim       :: phosphate limitation
                0117 C     Fe_lim        :: iron limitation for phytoplankton
                0118 C     Fe_lim_diaz   :: iron limitation for diazotrophs
                0119 C     light_lim     :: light limitation
                0120 C     expkT         :: temperature function
                0121 C     Pc_m          :: light-saturated max photosynthesis rate for phyto
                0122 C     Pc_m_diaz     :: light-saturated max photosynthesis rate for diaz
                0123 C     theta_Fe      :: Chl:C ratio
                0124 C     theta_Fe_inv  :: C:Chl ratio
                0125 C     theta_Fe_max  :: Fe-replete maximum Chl:C ratio
                0126 C     irrk          :: nut-limited efficiency of algal photosystems
                0127 C     irr_inst      :: instantaneous light
                0128 C     irr_eff       :: effective irradiance
                0129 C     mld           :: mixed layer depth
                0130 C     mu            :: net carbon-specific growth rate for phyt
                0131 C     mu_diaz       :: net carbon-specific growth rate for diaz
                0132 C     PtoN          :: variable ratio of phosphorus to nitrogen
                0133 C     FetoN         :: variable ratio of iron to nitrogen
                0134 C     N_uptake      :: NO3 utilization by phytoplankton
                0135 C     N_fix         :: Nitrogen fixation by diazotrophs
                0136 C     N_den_pelag   :: Pelagic denitrification
                0137 C     N_den_benthic :: Benthic denitrification
                0138 C     P_uptake      :: PO4 utilization by phytoplankton
                0139 C     Fe_uptake     :: dissolved Fe utilization by phytoplankton
                0140 C     CaCO3_uptake  :: Calcium carbonate uptake for shell formation
                0141 C     CaCO3_diss    :: Calcium carbonate dissolution
                0142 C     G_CaCO3       :: tendency of calcium carbonate
                0143 C     DON_prod      :: production of dissolved organic nitrogen
                0144 C     DOP_prod      :: production of dissolved organic phosphorus
                0145 C     DON_remin     :: remineralization of dissolved organic nitrogen
                0146 C     DOP_remin     :: remineralization of dissolved organic phosphorus
                0147 C     O2_prod       :: production of oxygen
                0148 C     frac_exp      :: fraction of sinking particulate organic matter
                0149 C     N_spm         :: particulate sinking of nitrogen
                0150 C     P_spm         :: particulate sinking of phosphorus
                0151 C     Fe_spm        :: particulate sinking of iron
                0152 C     N_dvm         :: vertical transport of nitrogen by DVM
                0153 C     P_dvm         :: vertical transport of phosphorus by DVM
                0154 C     Fe_dvm        :: vertical transport of iron by DVM
                0155 C     N_recycle     :: recycling of newly-produced organic nitrogen
                0156 C     P_recycle     :: recycling of newly-produced organic phosphorus
                0157 C     Fe_recycle    :: recycling of newly-produced organic iron
                0158 C     N_reminp      :: remineralization of particulate organic nitrogen
                0159 C     P_reminp      :: remineralization of particulate organic nitrogen
                0160 C     Fe_reminp     :: remineralization of particulate organic iron
                0161 C     Fe_reminsum   :: iron remineralization and adsorption
                0162 C     N_remindvm    :: nitrogen remineralization due to diel vertical migration
                0163 C     P_remindvm    :: phosphorus remineralization due to diel vertical migration
                0164 C     Fe_remindvm   :: iron remineralization due to diel vertical migration
                0165 C     POC_flux      :: particulate organic carbon flux
                0166 C     NPP           :: net primary production
                0167 C     NCP           :: net community production
                0168 C     depth_l       :: depth of lower interface
                0169 C     *flux_u       :: "*" flux through upper interface
                0170 C     *flux_l       :: "*" flux through lower interface
                0171 C     zremin        :: remineralization lengthscale for nutrients
                0172 C     zremin_caco3  :: remineralization lengthscale for CaCO3
                0173 C     wsink         :: speed of sinking particles
                0174 C     POC_sed       :: flux of particulate organic carbon to sediment
                0175 C     Fe_sed        :: sediment iron efflux
                0176 C     kFe_eq_lig    :: iron ligand stability constant
                0177 C     FreeFe        :: free (unbound) iron concentration
                0178 C     Fe_ads_inorg  :: iron adsorption
                0179 C     Fe_ads_org    :: iron adsorption  (2nd order)
                0180 C     NO3_btmflx    :: bottom flux of nitrate
                0181 C     PO4_btmflx    :: bottom flux of phosphate
                0182 C     O2_btmflx     :: bottom flux of oxygen
                0183 C     Fe_burial     :: flux of iron to sediment as particulate
                0184 C
                0185       INTEGER i,j,k
                0186       _RL Phy_lg_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0187       _RL Phy_sm_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0188       _RL Phy_diaz_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0189       _RL NO3_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0190       _RL PO4_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0191       _RL Fe_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0192       _RL Fe_lim_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0193       _RL light_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0194       _RL expkT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0195       _RL Pc_m
                0196       _RL Pc_m_diaz
                0197       _RL theta_Fe_max
                0198       _RL theta_Fe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0199       _RL theta_Fe_inv(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0200       _RL irrk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0201       _RL irr_eff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
a284455135 Matt*0202       _RL irr_inst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e0f9a7ba0b Matt*0203       _RL mld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0204       _RL mu(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0205       _RL mu_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0206       _RL PtoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0207       _RL FetoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0208       _RL N_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0209       _RL N_fix(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0210       _RL N_den_pelag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0211       _RL N_den_benthic(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0212       _RL P_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0213       _RL Fe_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0214       _RL CaCO3_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0215       _RL CaCO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0216       _RL G_CACO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0217       _RL DON_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0218       _RL DOP_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0219       _RL DON_remin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0220       _RL DOP_remin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0221       _RL O2_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0222       _RL frac_exp
                0223       _RL N_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0224       _RL P_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0225       _RL Fe_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0226       _RL N_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0227       _RL P_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0228       _RL Fe_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0229       _RL N_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0230       _RL P_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0231       _RL Fe_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0232       _RL N_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0233       _RL P_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0234       _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0235       _RL Fe_reminsum(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0236       _RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0237       _RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0238       _RL Fe_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0239       _RL POC_flux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0240       _RL NPP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0241       _RL NCP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0242       _RL depth_l
                0243       _RL PONflux_u
                0244       _RL PONflux_l
                0245       _RL POPflux_u
                0246       _RL POPflux_l
                0247       _RL PFEflux_u
                0248       _RL PFEflux_l
                0249       _RL CaCO3flux_u
                0250       _RL CaCO3flux_l
                0251       _RL zremin
                0252       _RL zremin_caco3
                0253       _RL wsink
                0254       _RL POC_sed
                0255       _RL Fe_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0256       _RL kFe_eq_lig
                0257       _RL FreeFe
                0258       _RL Fe_ads_inorg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0259       _RL Fe_ads_org(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0260       _RL log_btm_flx
                0261       _RL NO3_btmflx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0262       _RL PO4_btmflx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0263       _RL O2_btmflx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0264       _RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0265       _RL no3_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0266       _RL po4_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0267       _RL fe_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0268       _RL o2_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0269       _RL don_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0270       _RL dop_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0271 #ifdef ML_MEAN_PHYTO
                0272       _RL tmp_p_sm_ML
                0273       _RL tmp_p_lg_ML
                0274       _RL tmp_p_diaz_ML
                0275       _RL tmp_ML
                0276 #endif
                0277 #ifdef USE_BLING_DVM
                0278       INTEGER tmp_dvm
                0279       _RL o2_upper
                0280       _RL o2_lower
                0281       _RL dz_upper
                0282       _RL dz_lower
                0283       _RL temp_upper
                0284       _RL temp_lower
                0285       _RL z_dvm_regr
                0286       _RL frac_migr
                0287       _RL fdvm_migr
                0288       _RL fdvm_stat
                0289       _RL fdvmn_vint
                0290       _RL fdvmp_vint
                0291       _RL fdvmfe_vint
                0292       _RL z_dvm
                0293       _RL dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0294       _RL x_erfcc,z_erfcc,t_erfcc,erfcc
                0295 #endif
                0296 #ifdef USE_SIBLING
                0297       _RL Si_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0298       _RL Si_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
00fa2d4ddd mmaz*0299       _RL Si_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e0f9a7ba0b Matt*0300       _RL Si_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0301       _RL Si_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
00fa2d4ddd mmaz*0302       _RL SitoN_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e0f9a7ba0b Matt*0303       _RL Phy_diatom_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0304       _RL frac_diss_Si(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0305       _RL Pc_m_diatom
                0306       _RL mu_diatom(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0307       _RL zremin_Si
00fa2d4ddd mmaz*0308       _RL Siflux_u
                0309       _RL Siflux_l
e0f9a7ba0b Matt*0310       _RL si_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0311 #endif
                0312 #ifdef ADVECT_PHYTO
                0313       _RL phy_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0314 #endif
a284455135 Matt*0315 #ifdef SIZE_DEP_LIM
                0316       _RL light_lim_sm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0317       _RL light_lim_lg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0318       _RL NO3_lim_sm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0319       _RL NO3_lim_lg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0320       _RL PO4_lim_sm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0321       _RL PO4_lim_lg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0322       _RL Fe_lim_sm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0323       _RL Fe_lim_lg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0324       _RL Pc_m_sm
                0325       _RL Pc_m_lg
                0326       _RL mu_sm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0327       _RL mu_lg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0328 #endif
7c50f07931 Mart*0329 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0330 C     tkey   :: tape key (tile dependent)
                0331 C     ijkkey :: tape key (depends on i,j,k indices and tiles)
                0332       INTEGER tkey, ijkkey
7c50f07931 Mart*0333 #endif
e0f9a7ba0b Matt*0334 CEOP
                0335
a284455135 Matt*0336 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0337       tkey = bi + (bj - 1)*nSx + (ikey_dynamics - 1)*nSx*nSy
a284455135 Matt*0338 #endif
                0339
                0340 C  Initialize output and diagnostics
6df47d206d Jean*0341        DO j=1-OLy,sNy+OLy
                0342          DO i=1-OLx,sNx+OLx
e0f9a7ba0b Matt*0343             mld(i,j)             = 0. _d 0
                0344             NO3_btmflx(i,j)      = 0. _d 0
                0345             PO4_btmflx(i,j)      = 0. _d 0
                0346             O2_btmflx(i,j)       = 0. _d 0
                0347             Fe_burial(i,j)       = 0. _d 0
6df47d206d Jean*0348          ENDDO
e0f9a7ba0b Matt*0349        ENDDO
                0350        DO k=1,Nr
6df47d206d Jean*0351          DO j=1-OLy,sNy+OLy
                0352            DO i=1-OLx,sNx+OLx
e0f9a7ba0b Matt*0353               G_DIC(i,j,k)          = 0. _d 0
                0354               G_ALK(i,j,k)          = 0. _d 0
                0355               G_O2(i,j,k)           = 0. _d 0
                0356               G_Fe(i,j,k)           = 0. _d 0
                0357               G_PO4(i,j,k)          = 0. _d 0
                0358               G_DOP(i,j,k)          = 0. _d 0
                0359               G_NO3(i,j,k)          = 0. _d 0
                0360               G_DON(i,j,k)          = 0. _d 0
                0361               G_CaCO3(i,j,k)        = 0. _d 0
                0362               Phy_lg_local(i,j,k)   = 0. _d 0
                0363               Phy_sm_local(i,j,k)   = 0. _d 0
                0364               Phy_diaz_local(i,j,k) = 0. _d 0
                0365               NO3_lim(i,j,k)        = 0. _d 0
                0366               PO4_lim(i,j,k)        = 0. _d 0
                0367               Fe_lim(i,j,k)         = 0. _d 0
                0368               Fe_lim_diaz(i,j,k)    = 0. _d 0
                0369               light_lim(i,j,k)      = 0. _d 0
                0370               expkT(i,j,k)          = 0. _d 0
                0371               theta_Fe(i,j,k)       = 0. _d 0
                0372               irrk(i,j,k)           = 0. _d 0
                0373               irr_inst(i,j,k)       = 0. _d 0
                0374               irr_eff(i,j,k)        = 0. _d 0
                0375               mu(i,j,k)             = 0. _d 0
                0376               mu_diaz(i,j,k)        = 0. _d 0
                0377               PtoN(i,j,k)           = 0. _d 0
                0378               FetoN(i,j,k)          = 0. _d 0
                0379               N_uptake(i,j,k)       = 0. _d 0
                0380               N_fix(i,j,k)          = 0. _d 0
                0381               N_den_pelag(i,j,k)    = 0. _d 0
                0382               N_den_benthic(i,j,k)  = 0. _d 0
                0383               P_uptake(i,j,k)       = 0. _d 0
                0384               Fe_uptake(i,j,k)      = 0. _d 0
                0385               CaCO3_uptake(i,j,k)   = 0. _d 0
                0386               CaCO3_diss(i,j,k)     = 0. _d 0
                0387               DON_prod(i,j,k)       = 0. _d 0
                0388               DOP_prod(i,j,k)       = 0. _d 0
                0389               DON_remin(i,j,k)      = 0. _d 0
                0390               DOP_remin(i,j,k)      = 0. _d 0
                0391               O2_prod(i,j,k)        = 0. _d 0
                0392               N_spm(i,j,k)          = 0. _d 0
                0393               P_spm(i,j,k)          = 0. _d 0
                0394               Fe_spm(i,j,k)         = 0. _d 0
                0395               N_dvm(i,j,k)          = 0. _d 0
                0396               P_dvm(i,j,k)          = 0. _d 0
                0397               Fe_dvm(i,j,k)         = 0. _d 0
                0398               N_recycle(i,j,k)      = 0. _d 0
                0399               P_recycle(i,j,k)      = 0. _d 0
                0400               Fe_recycle(i,j,k)     = 0. _d 0
                0401               N_reminp(i,j,k)       = 0. _d 0
                0402               P_reminp(i,j,k)       = 0. _d 0
                0403               Fe_reminp(i,j,k)      = 0. _d 0
                0404               Fe_reminsum(i,j,k)    = 0. _d 0
                0405               N_remindvm(i,j,k)     = 0. _d 0
                0406               P_remindvm(i,j,k)     = 0. _d 0
                0407               Fe_remindvm(i,j,k)    = 0. _d 0
                0408               NCP(i,j,k)            = 0. _d 0
                0409               Fe_sed(i,j,k)         = 0. _d 0
                0410               Fe_ads_org(i,j,k)     = 0. _d 0
                0411               Fe_ads_inorg(i,j,k)   = 0. _d 0
                0412 #ifdef USE_BLING_DVM
                0413               dvm(i,j,k)            = 0. _d 0
                0414 #endif
                0415 #ifdef USE_SIBLING
00fa2d4ddd mmaz*0416               G_SI(i,j,k)             = 0. _d 0
e0f9a7ba0b Matt*0417               Si_lim(i,j,k)           = 0. _d 0
00fa2d4ddd mmaz*0418               Si_spm(i,j,k)           = 0. _d 0
e0f9a7ba0b Matt*0419               Si_uptake(i,j,k)        = 0. _d 0
                0420               Si_recycle(i,j,k)       = 0. _d 0
                0421               Si_reminp(i,j,k)        = 0. _d 0
00fa2d4ddd mmaz*0422               SitoN_uptake(i,j,k)     = 0. _d 0
e0f9a7ba0b Matt*0423               Phy_diatom_local(i,j,k) = 0. _d 0
                0424               frac_diss_Si(i,j,k)     = 0. _d 0
                0425               mu_diatom(i,j,k)        = 0. _d 0
                0426 #endif
6df47d206d Jean*0427            ENDDO
                0428          ENDDO
e0f9a7ba0b Matt*0429        ENDDO
                0430
a284455135 Matt*0431 C-----------------------------------------------------------
                0432 C  avoid negative nutrient concentrations that can result from
                0433 C  advection when low concentrations
e0f9a7ba0b Matt*0434
                0435 #ifdef BLING_NO_NEG
82e538d851 aver*0436 # ifdef ALLOW_AUTODIFF_TAMC
                0437 CADJ STORE PTR_O2, PTR_FE, PTR_PO4, PTR_DOP, PTR_NO3, PTR_DON
                0438 #  ifdef USE_SIBLING
                0439 CADJ &   , PTR_SI
                0440 #  endif
                0441 #  ifdef ADVECT_PHYTO
                0442 CADJ &   , PTR_PHY
                0443 #  endif
                0444 CADJ &     = comlev1_bibj, key=tkey, kind=isbyte
                0445 # endif /* ALLOW_AUTODIFF_TAMC */
e0f9a7ba0b Matt*0446       CALL BLING_MIN_VAL(PTR_O2,  1. _d -11, o2_adj, bi, bj)
                0447       CALL BLING_MIN_VAL(PTR_FE,  1. _d -11, fe_adj, bi, bj)
                0448       CALL BLING_MIN_VAL(PTR_PO4, 1. _d -8, po4_adj, bi, bj)
                0449       CALL BLING_MIN_VAL(PTR_DOP, 1. _d -11, dop_adj, bi, bj)
                0450       CALL BLING_MIN_VAL(PTR_NO3, 1. _d -7, no3_adj, bi, bj)
                0451       CALL BLING_MIN_VAL(PTR_DON, 1. _d -11, don_adj, bi, bj)
                0452 # ifdef USE_SIBLING
                0453       CALL BLING_MIN_VAL(PTR_SI , 1. _d -7, si_adj, bi, bj)
                0454 # endif
                0455 # ifdef ADVECT_PHYTO
                0456       CALL BLING_MIN_VAL(PTR_PHY, 1. _d -11, phy_adj, bi, bj)
                0457 # endif
                0458 #endif
                0459
a284455135 Matt*0460 C-----------------------------------------------------------
                0461 C  Phytoplankton size classes
e0f9a7ba0b Matt*0462
6df47d206d Jean*0463       DO k=1,Nr
                0464        DO j=jmin,jmax
                0465         DO i=imin,imax
e0f9a7ba0b Matt*0466 #ifdef ADVECT_PHYTO
                0467           Phy_lg_local(i,j,k)   = PTR_PHY(i,j,k)*phyto_lg(i,j,k,bi,bj)
                0468           Phy_sm_local(i,j,k)   = PTR_PHY(i,j,k)*phyto_sm(i,j,k,bi,bj)
                0469           Phy_diaz_local(i,j,k) = PTR_PHY(i,j,k)*phyto_diaz(i,j,k,bi,bj)
                0470 #else
                0471           Phy_lg_local(i,j,k)   = phyto_lg(i,j,k,bi,bj)
                0472           Phy_sm_local(i,j,k)   = phyto_sm(i,j,k,bi,bj)
                0473           Phy_diaz_local(i,j,k) = phyto_diaz(i,j,k,bi,bj)
                0474 #endif
                0475         ENDDO
                0476        ENDDO
6df47d206d Jean*0477       ENDDO
e0f9a7ba0b Matt*0478
a284455135 Matt*0479 C-----------------------------------------------------------
                0480 C  Mixed layer depth calculation for light, phytoplankton and dvm.
                0481 C  Do not need to calculate if not using ML_MEAN_LIGHT, ML_MEAN_PHYTO,
                0482 C  and USE_BLING_DVM
                0483 C  (with BLING_ADJOINT_SAFE flag, USE_BLING_DVM is undefined)
e0f9a7ba0b Matt*0484
                0485 #if ( defined (ML_MEAN_LIGHT) || \
                0486       defined (ML_MEAN_PHYTO) || \
                0487       defined (USE_BLING_DVM) )
6df47d206d Jean*0488       CALL BLING_MIXEDLAYER(
e0f9a7ba0b Matt*0489      U                         mld,
                0490      I                         bi, bj, imin, imax, jmin, jmax,
                0491      I                         myTime, myIter, myThid)
fe1862e69b Mart*0492 # ifdef ALLOW_AUTODIFF_TAMC
                0493 C     save result of s/r bling_mixedlayer to avoid calling it in AD-code
edb6656069 Mart*0494 CADJ STORE mld = comlev1_bibj, key=tkey, kind=isbyte
fe1862e69b Mart*0495 # endif
e0f9a7ba0b Matt*0496 #endif
                0497
a284455135 Matt*0498 C  Phytoplankton mixing
                0499 C  The mixed layer is assumed to homogenize vertical gradients of phytoplankton.
                0500 C  This allows for basic Sverdrup dynamics in a qualitative sense.
                0501 C  This has not been thoroughly tested, and care should be
                0502 C  taken with its interpretation.
e0f9a7ba0b Matt*0503
                0504 #ifdef ML_MEAN_PHYTO
82e538d851 aver*0505 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0506 CADJ STORE Phy_sm_local   = comlev1_bibj, key=tkey, kind=isbyte
                0507 CADJ STORE Phy_lg_local   = comlev1_bibj, key=tkey, kind=isbyte
                0508 CADJ STORE Phy_diaz_local = comlev1_bibj, key=tkey, kind=isbyte
82e538d851 aver*0509 # endif
e0f9a7ba0b Matt*0510       DO j=jmin,jmax
                0511        DO i=imin,imax
                0512
                0513         tmp_p_sm_ML = 0. _d 0
                0514         tmp_p_lg_ML = 0. _d 0
                0515         tmp_p_diaz_ML = 0. _d 0
                0516         tmp_ML  = 0. _d 0
                0517
                0518         DO k=1,Nr
                0519
a284455135 Matt*0520          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*0521           IF ((-rf(k+1) .le. mld(i,j)).and.
                0522      &               (-rf(k+1).lt.MLmix_max)) THEN
                0523            tmp_p_sm_ML = tmp_p_sm_ML+Phy_sm_local(i,j,k)*drF(k)
                0524      &                  *hFacC(i,j,k,bi,bj)
                0525            tmp_p_lg_ML = tmp_p_lg_ML+Phy_lg_local(i,j,k)*drF(k)
                0526      &                  *hFacC(i,j,k,bi,bj)
                0527            tmp_p_diaz_ML = tmp_p_diaz_ML+Phy_diaz_local(i,j,k)*drF(k)
                0528      &                  *hFacC(i,j,k,bi,bj)
                0529            tmp_ML = tmp_ML + drF(k)
                0530           ENDIF
                0531          ENDIF
                0532
                0533         ENDDO
                0534
                0535         DO k=1,Nr
                0536
a284455135 Matt*0537          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*0538           IF ((-rf(k+1) .le. mld(i,j)).and.
                0539      &               (-rf(k+1).lt.MLmix_max)) THEN
                0540
                0541            Phy_lg_local(i,j,k)   = max(1. _d -8,tmp_p_lg_ML/tmp_ML)
                0542            Phy_sm_local(i,j,k)   = max(1. _d -8,tmp_p_sm_ML/tmp_ML)
                0543            Phy_diaz_local(i,j,k) = max(1. _d -8,tmp_p_diaz_ML/tmp_ML)
                0544
                0545           ENDIF
                0546          ENDIF
                0547
                0548         ENDDO
                0549        ENDDO
                0550       ENDDO
82e538d851 aver*0551 # ifdef ALLOW_AUTODIFF_TAMC
                0552 CADJ STORE Phy_sm_local   = comlev1_bibj, key=tkey, kind=isbyte
                0553 CADJ STORE Phy_lg_local   = comlev1_bibj, key=tkey, kind=isbyte
                0554 CADJ STORE Phy_diaz_local = comlev1_bibj, key=tkey, kind=isbyte
                0555 # endif
e0f9a7ba0b Matt*0556 #endif
                0557
                0558 #ifdef USE_SIBLING
                0559       DO k=1,Nr
                0560        DO j=jmin,jmax
                0561         DO i=imin,imax
                0562          Phy_diatom_local(i,j,k) = Phy_lg_local(i,j,k)
                0563         ENDDO
                0564        ENDDO
                0565       ENDDO
                0566 #endif
                0567
a284455135 Matt*0568 C-----------------------------------------------------------
                0569 C  light availability for biological production
e0f9a7ba0b Matt*0570
6df47d206d Jean*0571       CALL BLING_LIGHT(
e0f9a7ba0b Matt*0572      I                    mld,
                0573      U                    irr_inst, irr_eff,
                0574      I                    bi, bj, imin, imax, jmin, jmax,
                0575      I                    myTime, myIter, myThid )
fe1862e69b Mart*0576 #ifdef ALLOW_AUTODIFF_TAMC
                0577 C     save results of s/r bling_light to avoid calling it in AD-code
edb6656069 Mart*0578 CADJ STORE irr_inst = comlev1_bibj, key=tkey, kind=isbyte
                0579 CADJ STORE irr_eff  = comlev1_bibj, key=tkey, kind=isbyte
fe1862e69b Mart*0580 #endif
e0f9a7ba0b Matt*0581
a284455135 Matt*0582 C  phytoplankton photoadaptation to local light level
                0583 C  This represents the fact that phytoplankton cells are adapted
                0584 C  to the averaged light field to which they've been exposed over
                0585 C  their lifetimes, rather than the instantaneous light. The
                0586 C  timescale is set by gamma_irr_mem.
e0f9a7ba0b Matt*0587
                0588       DO k=1,Nr
                0589        DO j=jmin,jmax
                0590         DO i=imin,imax
                0591
                0592           irr_mem(i,j,k,bi,bj) = irr_mem(i,j,k,bi,bj) +
                0593      &           (irr_eff(i,j,k) - irr_mem(i,j,k,bi,bj))*
                0594      &           min( 1. _d 0, gamma_irr_mem*PTRACERS_dTLev(k) )
                0595
                0596         ENDDO
                0597        ENDDO
                0598       ENDDO
                0599
a284455135 Matt*0600 C  Nutrient uptake and partitioning between organic pools
e0f9a7ba0b Matt*0601
                0602       DO k=1,Nr
                0603        DO j=jmin,jmax
                0604         DO i=imin,imax
a284455135 Matt*0605          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*0606
82e538d851 aver*0607 C--   Compute 3d fields that used more than once.
a284455135 Matt*0608 C ---------------------------------------------------------------------
                0609 C  First, calculate the limitation terms for NUT and Fe, and the
                0610 C  Fe-limited Chl:C maximum. The light-saturated maximal photosynthesis
                0611 C  rate term (Pc_m) is simply the product of a prescribed maximal
                0612 C  photosynthesis rate (Pc_0), the Eppley temperature dependence, and a
                0613 C  resource limitation term. The iron limitation term has a lower limit
                0614 C  of Fe_lim_min and is scaled by (k_Fe2P + Fe2P_max) / Fe2P_max so that
                0615 C  it approaches 1 as Fe approaches infinity. Thus, it is of comparable
                0616 C  magnitude to the macro-nutrient limitation term.
e0f9a7ba0b Matt*0617
a284455135 Matt*0618 C  Macro-nutrient limitation
e0f9a7ba0b Matt*0619
                0620           NO3_lim(i,j,k) = PTR_NO3(i,j,k)/(PTR_NO3(i,j,k)+k_NO3)
                0621
                0622           PO4_lim(i,j,k) = PTR_PO4(i,j,k)/(PTR_PO4(i,j,k)+k_PO4)
                0623
a284455135 Matt*0624 C  Iron limitation
e0f9a7ba0b Matt*0625
                0626           Fe_lim(i,j,k) = PTR_FE(i,j,k)
                0627      &                    / (PTR_FE(i,j,k)+k_Fe_2d(i,j,bi,bj))
                0628
                0629           Fe_lim_diaz(i,j,k) = PTR_FE(i,j,k)
                0630      &                         / (PTR_FE(i,j,k)+k_Fe_diaz_2d(i,j,bi,bj))
                0631
                0632 #ifdef USE_SIBLING
                0633           Si_lim(i,j,k) = PTR_Si(i,j,k)/(PTR_Si(i,j,k)+k_Si)
                0634 #endif
                0635
9edc0e3a85 aver*0636 #ifdef SIZE_DEP_LIM
82e538d851 aver*0637           NO3_lim_sm(i,j,k) = PTR_NO3(i,j,k)/(PTR_NO3(i,j,k)+k_NO3_sm)
                0638
                0639           PO4_lim_sm(i,j,k) = PTR_PO4(i,j,k)/(PTR_PO4(i,j,k)+k_PO4_sm)
                0640
                0641           Fe_lim_sm(i,j,k) = PTR_FE(i,j,k)/(PTR_FE(i,j,k)+k_Fe_sm)
                0642
                0643           NO3_lim_lg(i,j,k) = PTR_NO3(i,j,k)/(PTR_NO3(i,j,k)+k_NO3_lg)
                0644
                0645           PO4_lim_lg(i,j,k) = PTR_PO4(i,j,k)/(PTR_PO4(i,j,k)+k_PO4_lg)
                0646
                0647           Fe_lim_lg(i,j,k) = PTR_FE(i,j,k)/(PTR_FE(i,j,k)+k_Fe_lg)
9edc0e3a85 aver*0648 #endif
82e538d851 aver*0649 C  NB The temperature function of Eppley (1972) is used to represent the
                0650 C  effects of temperature on uptake, synthesis and ecosystem function.
                0651 C  The reference value is 0 C. Note that this differs slightly from the
                0652 C  more common Arrhenius form of the equation, but by a tiny amount (less
                0653 C  than 4% difference over the range -2 to 32 C for an Arrhenius constant
                0654 C  (Ae) of 5150).
                0655
a284455135 Matt*0656 C ---------------------------------------------------------------------
                0657 C  Temperature functionality of growth and grazing
e0f9a7ba0b Matt*0658
a284455135 Matt*0659 C  NB The temperature function of Eppley (1972) is used to represent the
                0660 C  effects of temperature on uptake, synthesis and ecosystem function.
                0661 C  The reference value is 0 C. Note that this differs slightly from the
                0662 C  more common Arrhenius form of the equation, but by a tiny amount (less
                0663 C  than 4% difference over the range -2 to 32 C for an Arrhenius constant
                0664 C  (Ae) of 5150).
e0f9a7ba0b Matt*0665
                0666           expkT(i,j,k) = exp(kappa_eppley * theta(i,j,k,bi,bj))
                0667
a284455135 Matt*0668 C  The light-saturated maximal photosynthesis rate term (pc_m) is
                0669 C  the product of a prescribed maximal photosynthesis rate (pc_0), the
                0670 C  Eppley temperature dependence, and a co-limitation term (the product
                0671 C  of Michaelis-Menten no3-limitation, po4-limitation, or iron-limitation).
                0672 C  Note this was altered from the Liebig limitation of BLINGv0, following
                0673 C  the arguments of Danger (Oikos, 2008) and Harpole (Ecol. Letters 2011).
                0674 C  Taking the cube root, to provide the geometric mean.
e0f9a7ba0b Matt*0675
a284455135 Matt*0676 C  Geometric mean is causing instability in b-SOSE. Returning to Liebig
e0f9a7ba0b Matt*0677
                0678 #ifdef MIN_NUT_LIM
                0679           Pc_m = Pc_0_2d(i,j,bi,bj) * expkT(i,j,k)
                0680      &           * min(NO3_lim(i,j,k), PO4_lim(i,j,k), Fe_lim(i,j,k))
                0681      &           * maskC(i,j,k,bi,bj)
                0682 #else
                0683           Pc_m = Pc_0_2d(i,j,bi,bj) * expkT(i,j,k)
                0684      &           * (NO3_lim(i,j,k) * PO4_lim(i,j,k) * Fe_lim(i,j,k))
                0685      &           **(1./3.) * maskC(i,j,k,bi,bj)
                0686 #endif
                0687
a284455135 Matt*0688 C  Diazotrophs are assumed to be more strongly temperature sensitive,
                0689 C  given their observed restriction to relatively warm waters. Presumably
                0690 C  this is because of the difficulty of achieving N2 fixation in an oxic
                0691 C  environment. Thus, they have lower pc_0 and higher kappa_eppley.
                0692 C  Taking the square root, to provide the geometric mean.
e0f9a7ba0b Matt*0693
                0694 #ifdef MIN_NUT_LIM
                0695           Pc_m_diaz = Pc_0_diaz_2d(i,j,bi,bj)
                0696      &           * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
                0697      &           * min(PO4_lim(i,j,k), Fe_lim_diaz(i,j,k))
                0698      &           * maskC(i,j,k,bi,bj)
                0699 #else
                0700           Pc_m_diaz = Pc_0_diaz_2d(i,j,bi,bj)
                0701      &           * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
                0702      &           * (PO4_lim(i,j,k) * Fe_lim_diaz(i,j,k))**(1./2.)
                0703      &           * maskC(i,j,k,bi,bj)
                0704 #endif
                0705
a284455135 Matt*0706 C  Diatoms are limited by the abundance of silica
                0707 C  Note: this is not in original SiBLING code
e0f9a7ba0b Matt*0708
                0709 #ifdef USE_SIBLING
                0710 # ifdef MIN_NUT_LIM
                0711           Pc_m_diatom = Pc_0_2d(i,j,bi,bj) * expkT(i,j,k)
                0712      &           * min(NO3_lim(i,j,k), PO4_lim(i,j,k), Fe_lim(i,j,k),
                0713      &           Si_lim(i,j,k)) * maskC(i,j,k,bi,bj)
                0714 # else
                0715           Pc_m_diatom = Pc_0_2d(i,j,bi,bj) * expkT(i,j,k)
                0716      &           * (NO3_lim(i,j,k) * PO4_lim(i,j,k) * Fe_lim(i,j,k)
                0717      &           * Si_lim(i,j,k)) **(1./4.) * maskC(i,j,k,bi,bj)
                0718 # endif
                0719 #endif
                0720
a284455135 Matt*0721 C  (Pc_m and Pc_m_diaz crash adjoint if get too small)
e0f9a7ba0b Matt*0722 #ifdef BLING_ADJOINT_SAFE
                0723           Pc_m        = max(Pc_m       ,maskC(i,j,k,bi,bj)*1. _d -15)
                0724           Pc_m_diaz   = max(Pc_m_diaz  ,maskC(i,j,k,bi,bj)*1. _d -15)
                0725 # ifdef USE_SIBLING
                0726           Pc_m_diatom = max(Pc_m_diatom,maskC(i,j,k,bi,bj)*1. _d -15)
                0727 # endif
                0728 #endif
                0729
a284455135 Matt*0730 C ---------------------------------------------------------------------
                0731 C  Fe limitation reduces the maximum achievable Chl:C ratio (theta_Fe)
                0732 C  below a prescribed, Fe-replete maximum value (theta_Fe_max),
                0733 C  to approach a prescribed minimum Chl:C (theta_Fe_min) under extreme
                0734 C  Fe-limitation.
e0f9a7ba0b Matt*0735
                0736           theta_Fe_max = theta_Fe_max_lo+
                0737      &                  (theta_Fe_max_hi-theta_Fe_max_lo)*Fe_lim(i,j,k)
                0738
                0739           theta_Fe(i,j,k) = theta_Fe_max
                0740      &                  / (1. _d 0 + alpha_photo_2d(i,j,bi,bj)
                0741      &                  *theta_Fe_max
                0742      &                  *irr_mem(i,j,k,bi,bj)/(epsln + 2. _d 0*Pc_m))
                0743
a284455135 Matt*0744 C ---------------------------------------------------------------------
                0745 C  Nutrient-limited efficiency of algal photosystems, irrk, is calculated
                0746 C  with the iron limitation term included as a multiplier of the
                0747 C  theta_Fe_max to represent the importance of Fe in forming chlorophyll
                0748 C  accessory antennae, which do not affect the Chl:C but still affect the
                0749 C  phytoplankton ability to use light (eg Stzrepek & Harrison, Nature 2004).
                0750 C  Use same thetamax for diazotrophs, since their extra Fe requirement
                0751 C  is for nitrogenase, not for the photosystem.
e0f9a7ba0b Matt*0752
                0753           irrk(i,j,k) = Pc_m
                0754      &             /(epsln + alpha_photo_2d(i,j,bi,bj)*theta_Fe_max)
                0755      &             + irr_mem(i,j,k,bi,bj)/2. _d 0
                0756
                0757           light_lim(i,j,k) = ( 1. _d 0 - exp(-irr_eff(i,j,k)
                0758      &               /(epsln + irrk(i,j,k))))
                0759
a284455135 Matt*0760 C  Carbon-specific photosynthesis rate
e0f9a7ba0b Matt*0761
                0762           mu(i,j,k) = Pc_m * light_lim(i,j,k)
                0763
                0764 #ifdef USE_SIBLING
                0765           mu_diatom(i,j,k) = Pc_m_diatom * light_lim(i,j,k)
                0766 #endif
                0767
a284455135 Matt*0768 C  Temperature threshold for diazotrophs
e0f9a7ba0b Matt*0769
                0770           IF (theta(i,j,k,bi,bj) .gt. 14) THEN
                0771            mu_diaz(i,j,k) = Pc_m_diaz * light_lim(i,j,k)
                0772           ELSE
                0773            mu_diaz(i,j,k) = 0. _d 0
                0774           ENDIF
                0775
a284455135 Matt*0776 C  Stoichiometry
e0f9a7ba0b Matt*0777
                0778           PtoN(i,j,k) = PtoN_min + (PtoN_max - PtoN_min) *
                0779      &                  PTR_PO4(i,j,k) / (k_PtoN + PTR_PO4(i,j,k))
                0780
                0781           FetoN(i,j,k) = FetoN_min + (FetoN_max - FetoN_min) *
                0782      &                   PTR_FE(i,j,k) / (k_FetoN + PTR_FE(i,j,k))
                0783
                0784 #ifdef USE_SIBLING
a284455135 Matt*0785 C  Si uptake by diatoms varies with both Si-limitation and growth-rate
                0786 C  limitation by other factors: increasing Si limitation decreases Si
                0787 C  uptake, but growth limitation by other factors under high Si conditions
                0788 C  results in increase Si:P of uptake.
e0f9a7ba0b Matt*0789
6ffd1aa797 Jean*0790           SitoN_uptake(i,j,k) = min(SitoN_uptake_max,
00fa2d4ddd mmaz*0791      &                   Si_lim(i,j,k) * max(SitoN_uptake_min,
                0792      &                   (1. / epsln + SitoN_uptake_scale + min(
e0f9a7ba0b Matt*0793      &                   PO4_lim(i,j,k), Fe_lim(i,j,k) ) * (1. -
a284455135 Matt*0794      &              exp(-irr_eff(i,j,k)/(epsln + irrk(i,j,k) )) ))
00fa2d4ddd mmaz*0795      &                   **SitoN_uptake_exp))
e0f9a7ba0b Matt*0796 #endif
                0797
a284455135 Matt*0798 C  Nutrient uptake
e0f9a7ba0b Matt*0799
                0800           N_uptake(i,j,k) = mu(i,j,k) * (Phy_sm_local(i,j,k)
                0801      &                      + Phy_lg_local(i,j,k))
                0802
a284455135 Matt*0803 #ifdef SIZE_DEP_LIM
                0804           Pc_m_sm = Pc_0_2d(i,j,bi,bj) * expkT(i,j,k)
                0805      &           * min(NO3_lim_sm(i,j,k), PO4_lim_sm(i,j,k),
                0806      &           Fe_lim_sm(i,j,k)) * maskC(i,j,k,bi,bj)
                0807
                0808           Pc_m_lg = Pc_0_2d(i,j,bi,bj) * expkT(i,j,k)
                0809      &           * min(NO3_lim_lg(i,j,k), PO4_lim_lg(i,j,k),
                0810      &           Fe_lim_lg(i,j,k)) * maskC(i,j,k,bi,bj)
                0811
                0812           light_lim_sm(i,j,k) = light_lim(i,j,k)
                0813
                0814           light_lim_lg(i,j,k) = light_lim(i,j,k)
                0815
                0816           mu_sm(i,j,k) = Pc_m_sm * light_lim_sm(i,j,k)
                0817           mu_lg(i,j,k) = Pc_m_lg * light_lim_lg(i,j,k)
                0818
                0819           N_uptake(i,j,k) = mu_sm(i,j,k) * Phy_sm_local(i,j,k)
                0820      &                      + mu_lg(i,j,k) * Phy_lg_local(i,j,k)
                0821 #endif /* SIZE_DEP_LIM */
                0822
e0f9a7ba0b Matt*0823           N_fix(i,j,k) = mu_diaz(i,j,k) * Phy_diaz_local(i,j,k)
                0824
                0825           P_uptake(i,j,k) = (N_uptake(i,j,k) +
                0826      &                      N_fix(i,j,k)) * PtoN(i,j,k)
                0827
                0828           Fe_uptake(i,j,k) = (N_uptake(i,j,k) +
                0829      &                       N_fix(i,j,k)) * FetoN(i,j,k)
                0830
                0831 #ifdef USE_SIBLING
                0832           Si_uptake(i,j,k) = mu(i,j,k) * Phy_diatom_local(i,j,k)
00fa2d4ddd mmaz*0833      &                       * SitoN_uptake(i,j,k)
e0f9a7ba0b Matt*0834 #endif
                0835
a284455135 Matt*0836 C  Calcium carbonate production
e0f9a7ba0b Matt*0837
a284455135 Matt*0838 C  Alkalinity is consumed through the production of CaCO3. Here, this is
                0839 C  simply a linear function of the implied growth rate of small
                0840 C  phytoplankton, which gave a reasonably good fit to the global
                0841 C  observational synthesis of Dunne (2009). This is consistent
                0842 C  with the findings of Jin et al. (GBC,2006).
e0f9a7ba0b Matt*0843
9edc0e3a85 aver*0844 #ifdef SIZE_DEP_LIM
a284455135 Matt*0845           CaCO3_uptake(i,j,k) = mu_sm(i,j,k) * Phy_sm_local(i,j,k)
                0846      &                          * phi_sm_2d(i,j,bi,bj) * CatoN
82e538d851 aver*0847 #else
                0848           CaCO3_uptake(i,j,k) = mu(i,j,k) * Phy_sm_local(i,j,k)
                0849      &                          * phi_sm_2d(i,j,bi,bj) * CatoN
a284455135 Matt*0850 #endif
                0851
e0f9a7ba0b Matt*0852          ENDIF
                0853         ENDDO
                0854        ENDDO
                0855       ENDDO
                0856
a284455135 Matt*0857 C ---------------------------------------------------------------------
                0858 C  Grazing is treated using the size-dependent grazing laws of
                0859 C  Dunne et al. (GBC, 2005), with alpha=1/3 (eq. 5b):
                0860 C   lambda_small = lambda0 * e^kT * (small/pivotal)
                0861 C   lambda_large = lambda0 * e^hT * (large/pivotal)^alpha
                0862 C  where small/large are the biomasses of small and large phytoplankton,
                0863 C  and lambda is mortality. The constants were determined empirically
                0864 C  from a global dataset.
                0865 C  Basically, small phytoplankton have a linear increase in mortality
                0866 C  with biomass, reflecting the rapid increase of small zooplankton as
                0867 C  growth rates increase, while large zooplankton respond weakly to
                0868 C  increases in the growth rates of large phytoplankton. This produces
                0869 C  a simple bloom dynamic.
                0870 C  Because diazotrophs have specialized grazers (O'Neil and Roman, 1992;
                0871 C  Mulholland, Biogeosciences 2007), they are treated independently from
                0872 C  other phytoplankton. This has an important impact, reducing direct
                0873 C  competition between diazotrophs and other phytoplankton.
                0874 C  Note that the phytoplankton concentrations here are not part of the
                0875 C  conserved nutrient pools: they are only used to keep track of the
                0876 C  biomass, which is required to calculate the uptake fluxes.
e0f9a7ba0b Matt*0877
                0878       DO k=1,Nr
                0879        DO j=jmin,jmax
                0880         DO i=imin,imax
a284455135 Matt*0881          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*0882
6df47d206d Jean*0883           Phy_lg_local(i,j,k) = Phy_lg_local(i,j,k) +
e0f9a7ba0b Matt*0884      &                Phy_lg_local(i,j,k)*(mu(i,j,k) - lambda_0
                0885      &                * expkT(i,j,k) *
                0886      &                (Phy_lg_local(i,j,k)/pivotal)**(1. / 3.))
                0887      &                * PTRACERS_dTLev(k)
                0888
                0889           Phy_sm_local(i,j,k) = Phy_sm_local(i,j,k) +
                0890      &                Phy_sm_local(i,j,k)*(mu(i,j,k) - lambda_0
                0891      &                * expkT(i,j,k) * (Phy_sm_local(i,j,k)/pivotal))
                0892      &                * PTRACERS_dTLev(k)
                0893
                0894           Phy_diaz_local(i,j,k) = Phy_diaz_local(i,j,k) +
82e538d851 aver*0895      &                Phy_diaz_local(i,j,k)*(mu_diaz(i,j,k)-20*lambda_0
e0f9a7ba0b Matt*0896      &                * expkT(i,j,k) * (Phy_diaz_local(i,j,k)/pivotal))
                0897      &                * PTRACERS_dTLev(k)
                0898
82e538d851 aver*0899 #ifdef ALLOW_AUTODIFF
                0900 C     avoid recomputation warnings
a284455135 Matt*0901          ENDIF
                0902         ENDDO
                0903        ENDDO
                0904       ENDDO
                0905 #ifdef ALLOW_AUTODIFF_TAMC
82e538d851 aver*0906 C     These stores avoid some recomputations, but may not be worth it.
                0907 cCADJ STORE Phy_sm_local   = comlev1_bibj, key=tkey, kind=isbyte
                0908 cCADJ STORE Phy_lg_local   = comlev1_bibj, key=tkey, kind=isbyte
                0909 cCADJ STORE Phy_diaz_local = comlev1_bibj, key=tkey, kind=isbyte
                0910 #endif
                0911       DO k=1,Nr
                0912        DO j=jmin,jmax
                0913         DO i=imin,imax
                0914          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
                0915 #endif
                0916            Phy_lg_local(i,j,k) =  max( epsln, Phy_lg_local(i,j,k) )
                0917            Phy_sm_local(i,j,k) =  max( epsln, Phy_sm_local(i,j,k) )
                0918            Phy_diaz_local(i,j,k) =  max( epsln, Phy_diaz_local(i,j,k) )
                0919          ENDIF
                0920         ENDDO
                0921        ENDDO
                0922       ENDDO
                0923
                0924 #if ( defined SIZE_DEP_LIM || defined USE_SIBLING )
                0925 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0926 CADJ STORE Phy_sm_local   = comlev1_bibj, key=tkey, kind=isbyte
                0927 CADJ STORE Phy_lg_local   = comlev1_bibj, key=tkey, kind=isbyte
82e538d851 aver*0928 CADJ STORE Phy_diaz_local = comlev1_bibj, key=tkey, kind=isbyte
                0929 # endif
a284455135 Matt*0930       DO k=1,Nr
                0931        DO j=jmin,jmax
                0932         DO i=imin,imax
                0933          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
82e538d851 aver*0934 # ifdef SIZE_DEP_LIM
a284455135 Matt*0935
                0936           Phy_lg_local(i,j,k) = Phy_lg_local(i,j,k) +
                0937      &                Phy_lg_local(i,j,k)*(mu_lg(i,j,k) - lambda_0
                0938      &                * expkT(i,j,k) *
                0939      &                (Phy_lg_local(i,j,k)/pivotal)**(1. / 3.))
                0940      &                * PTRACERS_dTLev(k)
                0941
                0942           Phy_sm_local(i,j,k) = Phy_sm_local(i,j,k) +
                0943      &                Phy_sm_local(i,j,k)*(mu_sm(i,j,k) - lambda_0
                0944      &                * expkT(i,j,k) * (Phy_sm_local(i,j,k)/pivotal))
                0945      &                * PTRACERS_dTLev(k)
                0946
82e538d851 aver*0947 # endif
a284455135 Matt*0948
82e538d851 aver*0949 # ifdef USE_SIBLING
a284455135 Matt*0950 C  For now: assuming that all large phyto are diatoms
                0951 C  Another option would be to have a separate diatom population and give
                0952 C  them a higher growth rate to be competitive. Would also need to add
                0953 C  a global variable "phyto_diatom"
e0f9a7ba0b Matt*0954           Phy_diatom_local(i,j,k) = Phy_lg_local(i,j,k)
                0955
                0956 c          Phy_diatom_local(i,j,k) = Phy_diatom_local(i,j,k) +
                0957 c     &                Phy_diatom_local(i,j,k)*(mu_diatom(i,j,k)-lambda_0
                0958 c     &                * expkT(i,j,k) * (Phy_diatom_local(i,j,k)/pivotal))
                0959 c     &                * PTRACERS_dTLev(k)
82e538d851 aver*0960 # endif
e0f9a7ba0b Matt*0961
                0962          ENDIF
                0963         ENDDO
                0964        ENDDO
                0965       ENDDO
a284455135 Matt*0966 #endif /* SIZE_DEP_LIM or USE_SIBLING */
e0f9a7ba0b Matt*0967
a284455135 Matt*0968 C  Separate loop for adjoint stores
                0969 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0970 CADJ STORE Phy_sm_local   = comlev1_bibj, key=tkey, kind=isbyte
                0971 CADJ STORE Phy_lg_local   = comlev1_bibj, key=tkey, kind=isbyte
                0972 CADJ STORE Phy_diaz_local = comlev1_bibj, key=tkey, kind=isbyte
a284455135 Matt*0973 #endif
e0f9a7ba0b Matt*0974
                0975       DO k=1,Nr
                0976        DO j=jmin,jmax
                0977         DO i=imin,imax
                0978
a284455135 Matt*0979          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*0980
                0981 #ifdef ADVECT_PHYTO
a284455135 Matt*0982 C  update tracer here
6df47d206d Jean*0983           PTR_PHY(i,j,k) = Phy_lg_local(i,j,k) + Phy_sm_local(i,j,k)
e0f9a7ba0b Matt*0984      &                     + Phy_diaz_local(i,j,k)
a284455135 Matt*0985 C  update fractional abundance
6df47d206d Jean*0986           phyto_lg(i,j,k,bi,bj) = Phy_lg_local(i,j,k)/PTR_PHY(i,j,k)
                0987           phyto_sm(i,j,k,bi,bj) = Phy_sm_local(i,j,k)/PTR_PHY(i,j,k)
                0988           phyto_diaz(i,j,k,bi,bj) = Phy_diaz_local(i,j,k)/PTR_PHY(i,j,k)
e0f9a7ba0b Matt*0989 #else
a284455135 Matt*0990 C  update biomass
6df47d206d Jean*0991           phyto_lg(i,j,k,bi,bj) = Phy_lg_local(i,j,k)
                0992           phyto_sm(i,j,k,bi,bj) = Phy_sm_local(i,j,k)
                0993           phyto_diaz(i,j,k,bi,bj) = Phy_diaz_local(i,j,k)
e0f9a7ba0b Matt*0994 #endif
                0995
a284455135 Matt*0996 C  use the diagnostic biomass to calculate the chl concentration
                0997 C  in mg/m3 (carbon = 12.01 g/mol)
e0f9a7ba0b Matt*0998
                0999           chl(i,j,k,bi,bj) = max(chl_min, CtoN * 12.01 * 1. _d 3 *
                1000      &           theta_Fe(i,j,k) *
                1001      &           (Phy_lg_local(i,j,k) + Phy_sm_local(i,j,k)
                1002      &           + Phy_diaz_local(i,j,k)))
                1003
                1004          ENDIF
                1005         ENDDO
                1006        ENDDO
                1007       ENDDO
82e538d851 aver*1008 #ifdef ALLOW_AUTODIFF_TAMC
                1009 CADJ STORE phyto_sm(:,:,:,bi,bj)   = comlev1_bibj, key=tkey, kind=isbyte
                1010 CADJ STORE phyto_lg(:,:,:,bi,bj)   = comlev1_bibj, key=tkey, kind=isbyte
                1011 CADJ STORE phyto_diaz(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                1012 #endif
e0f9a7ba0b Matt*1013
a284455135 Matt*1014 C ---------------------------------------------------------------------
                1015 C  Partitioning between organic pools
                1016
                1017 C  The uptake of nutrients is assumed to contribute to the growth of
                1018 C  phytoplankton, which subsequently die and are consumed by heterotrophs.
                1019 C  This can involve the transfer of nutrient elements between many
                1020 C  organic pools, both particulate and dissolved, with complex histories.
                1021 C  We take a simple approach here, partitioning the total uptake into two
                1022 C  fractions - sinking and non-sinking - as a function of temperature,
                1023 C  following Dunne et al. (2005).
                1024 C  Then, the non-sinking fraction is further subdivided, such that the
                1025 C  majority is recycled instantaneously to the inorganic nutrient pool,
                1026 C  representing the fast turnover of labile dissolved organic matter via
                1027 C  the microbial loop, and the remainder is converted to semi-labile
                1028 C  dissolved organic matter. Iron and macro-nutrient are treated
                1029 C  identically for the first step, but all iron is recycled
                1030 C  instantaneously in the second step (i.e. there is no dissolved organic
                1031 C  iron pool).
                1032 C  Assume that the fraction of uptake that produces sinking particulate
                1033 C  organic matter is the same for diazotrophic and non-diazotrophic
                1034 C  phytoplankton.
e0f9a7ba0b Matt*1035
                1036       DO k=1,Nr
                1037        DO j=jmin,jmax
                1038         DO i=imin,imax
                1039
a284455135 Matt*1040          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*1041
a284455135 Matt*1042 C  sinking fraction: particulate organic matter
                1043
                1044 #ifdef NEW_FRAC_EXP
9edc0e3a85 aver*1045 C  Fraction that does into POM is
                1046 C  (assume diaz has same export frac as small phyto)
                1047 C  e^kT * [phi_lg * B_lg/Btot + phi_sm * B_sm/Btot + phi_sm * Bdiaz/Bto]
a284455135 Matt*1048
                1049           frac_exp = (phi_sm_2d(i,j,bi,bj) * (phyto_sm(i,j,k,bi,bj) +
                1050      &      phyto_diaz(i,j,k,bi,bj)) + phi_lg_2d(i,j,bi,bj) *
                1051      &      phyto_lg(i,j,k,bi,bj)) /
                1052      &      (phyto_sm(i,j,k,bi,bj) + phyto_diaz(i,j,k,bi,bj) +
                1053      &      phyto_lg(i,j,k,bi,bj)) *
                1054      &      exp(kappa_remin * theta(i,j,k,bi,bj))
                1055
                1056 #else
                1057
9edc0e3a85 aver*1058 C  This is from blingv1, where B_sm and B_lg are diagnosed from
                1059 C  the value of mu...
e0f9a7ba0b Matt*1060           frac_exp = (phi_sm_2d(i,j,bi,bj) + phi_lg_2d(i,j,bi,bj) *
a284455135 Matt*1061      &           (mu(i,j,k)/(epsln + lambda_0*expkT(i,j,k)))**2)/
                1062      &           (1. + (mu(i,j,k)/
                1063      &           (epsln + lambda_0*expkT(i,j,k)))**2)*
e0f9a7ba0b Matt*1064      &           exp(kappa_remin * theta(i,j,k,bi,bj))
                1065
a284455135 Matt*1066 #endif
                1067
e0f9a7ba0b Matt*1068 #ifdef USE_BLING_DVM
                1069           N_dvm(i,j,k) = phi_dvm * frac_exp *
                1070      &             (N_uptake(i,j,k) + N_fix(i,j,k))
                1071
                1072           P_dvm(i,j,k) = phi_dvm * frac_exp * P_uptake(i,j,k)
                1073
                1074           Fe_dvm(i,j,k) = phi_dvm * frac_exp * Fe_uptake(i,j,k)
                1075
                1076           N_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
                1077      &                 (N_uptake(i,j,k) + N_fix(i,j,k))
                1078
                1079           P_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
                1080      &                 P_uptake(i,j,k)
                1081
                1082           Fe_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
                1083      &                  Fe_uptake(i,j,k)
                1084 #else
                1085           N_spm(i,j,k) = frac_exp * (N_uptake(i,j,k) + N_fix(i,j,k))
                1086           P_spm(i,j,k) = frac_exp * P_uptake(i,j,k)
                1087           Fe_spm(i,j,k) = frac_exp * Fe_uptake(i,j,k)
                1088
                1089           N_dvm(i,j,k) = 0
                1090           P_dvm(i,j,k) = 0
                1091           Fe_dvm(i,j,k) = 0
                1092 #endif
                1093
a284455135 Matt*1094 C  the remainder is divided between instantaneously recycled and
                1095 C  long-lived dissolved organic matter.
e0f9a7ba0b Matt*1096
                1097           DON_prod(i,j,k) = phi_DOM_2d(i,j,bi,bj)*(N_uptake(i,j,k)
                1098      &                      + N_fix(i,j,k) - N_spm(i,j,k)
                1099      &                      - N_dvm(i,j,k))
                1100
                1101           DOP_prod(i,j,k) = phi_DOM_2d(i,j,bi,bj)*(P_uptake(i,j,k)
                1102      &                      - P_spm(i,j,k) - P_dvm(i,j,k))
                1103
                1104           N_recycle(i,j,k) = N_uptake(i,j,k) + N_fix(i,j,k)
                1105      &                       - N_spm(i,j,k) - DON_prod(i,j,k)
                1106      &                       - N_dvm(i,j,k)
                1107
                1108           P_recycle(i,j,k) = P_uptake(i,j,k)
                1109      &                       - P_spm(i,j,k) - DOP_prod(i,j,k)
                1110      &                       - P_dvm(i,j,k)
                1111
                1112           Fe_recycle(i,j,k) = Fe_uptake(i,j,k)
                1113      &                        - Fe_spm(i,j,k) - Fe_dvm(i,j,k)
                1114
00fa2d4ddd mmaz*1115 #ifdef USE_SIBLING
6ffd1aa797 Jean*1116           frac_diss_Si(i,j,k) = exp(-SitoN_uptake(i,j,k) /
                1117      &                        (q_SitoN_diss *
00fa2d4ddd mmaz*1118      &                        exp(kappa_remin_Si * theta(i,j,k,bi,bj))))
                1119
                1120           Si_recycle(i,j,k) = Si_uptake(i,j,k) * frac_diss_Si(i,j,k)
                1121
                1122           Si_spm(i,j,k) = Si_uptake(i,j,k) - Si_recycle(i,j,k)
                1123 #endif
                1124
e0f9a7ba0b Matt*1125          ENDIF
                1126
                1127         ENDDO
                1128        ENDDO
                1129       ENDDO
                1130
82e538d851 aver*1131 CML#ifdef ALLOW_AUTODIFF_TAMC
                1132 CMLC     these stores avoid recomputing the previous loop, but avoiding
                1133 CMLC     that may not be necessary
                1134 CMLCADJ STORE N_spm     = comlev1_bibj, key = tkey, kind = isbyte
                1135 CMLCADJ STORE P_spm     = comlev1_bibj, key = tkey, kind = isbyte
                1136 CMLCADJ STORE Fe_spm    = comlev1_bibj, key = tkey, kind = isbyte
                1137 CMLCADJ STORE don_prod  = comlev1_bibj, key = tkey, kind = isbyte
                1138 CMLCADJ STORE n_recycle = comlev1_bibj, key = tkey, kind = isbyte
                1139 CML#endif
fe1862e69b Mart*1140
e0f9a7ba0b Matt*1141 C ---------------------------------------------------------------------
a284455135 Matt*1142 C  Remineralization
                1143 C  In general, the flux at the bottom of a grid cell should equal
                1144 C  Fb = (Ft + Prod*dz) / (1 + zremin*dz)
                1145 C  where Ft is the flux at the top, and prod*dz is the integrated
                1146 C  production of new sinking particles within the layer.
                1147 C  This is also multiplied by the fraction of the grid cell that is
                1148 C  not intercepted by the bottom at this level, according to the
                1149 C  subgridscale sediment interception fraction.
                1150 C  Since Ft=0 in the first layer,
e0f9a7ba0b Matt*1151
                1152 C$TAF LOOP = parallel
6df47d206d Jean*1153       DO j=jmin,jmax
e0f9a7ba0b Matt*1154 C$TAF LOOP = parallel
6df47d206d Jean*1155        DO i=imin,imax
e0f9a7ba0b Matt*1156
a284455135 Matt*1157 C  Initialize upper flux
e0f9a7ba0b Matt*1158
                1159         PONflux_u            = 0. _d 0
                1160         POPflux_u            = 0. _d 0
                1161         PFEflux_u            = 0. _d 0
                1162         CaCO3flux_u          = 0. _d 0
00fa2d4ddd mmaz*1163 #ifdef USE_SIBLING
                1164         Siflux_u             = 0. _d 0
                1165 #endif
e0f9a7ba0b Matt*1166
                1167         DO k=1,Nr
fe1862e69b Mart*1168 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1169          ijkkey = k + ( (i-1) + (j-1) * (2*OLx+sNx) ) * Nr
                1170      &          + (tkey - 1) * (2*OLx+sNx) * (2*OLy+sNy) * Nr
fe1862e69b Mart*1171 CADJ STORE PONflux_u   = comlev1_bibj_ijk, key = ijkkey, kind = isbyte
                1172 CADJ STORE POPflux_u   = comlev1_bibj_ijk, key = ijkkey, kind = isbyte
                1173 CADJ STORE PFEflux_u   = comlev1_bibj_ijk, key = ijkkey, kind = isbyte
                1174 CADJ STORE CaCO3flux_u = comlev1_bibj_ijk, key = ijkkey, kind = isbyte
                1175 # ifdef USE_SIBLING
82e538d851 aver*1176 CADJ STORE Siflux_u    = comlev1_bibj_ijk, key = ijkkey, kind = isbyte
fe1862e69b Mart*1177 # endif
                1178 #endif
e0f9a7ba0b Matt*1179
a284455135 Matt*1180 C  Initialization here helps taf
e0f9a7ba0b Matt*1181
                1182          Fe_ads_org(i,j,k)    = 0. _d 0
                1183
a284455135 Matt*1184          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*1185
a284455135 Matt*1186 C  Calculate the remineralization lengthscale matrix, zremin, a function
                1187 C  of z. Sinking rate (wsink) is constant over the upper wsink0_z metres,
                1188 C  then  increases linearly with depth.
                1189 C  The remineralization rate is a function of oxygen concentrations,
                1190 C  following a Holling type 2 dependence, decreasing to a minimum value
                1191 C  of remin_min. This is ad hoc, following work by Bianchi, Sarmiento,
                1192 C  Galbraith and Kwon (unpublished).
                1193 C  Sinking speed is evaluated at the bottom of the cell
e0f9a7ba0b Matt*1194
                1195           depth_l=-rF(k+1)
                1196           IF (depth_l .LE. wsink0z)  THEN
                1197            wsink = wsink0_2d(i,j,bi,bj)
                1198           ELSE
                1199            wsink = wsinkacc * (depth_l - wsink0z) + wsink0_2d(i,j,bi,bj)
                1200           ENDIF
                1201
                1202           zremin = gamma_POM_2d(i,j,bi,bj) * ( PTR_O2(i,j,k)**2 /
                1203      &               (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)
                1204      &               + remin_min )/(wsink + epsln)
                1205
                1206 C  Calcium remineralization relaxed toward the inverse of the
                1207 C  ca_remin_depth constant value as the calcite saturation approaches 0.
                1208
                1209           zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0 - min(1. _d 0,
                1210      &               omegaC(i,j,k,bi,bj) + epsln ))
                1211
00fa2d4ddd mmaz*1212 #ifdef USE_SIBLING
                1213           zremin_Si = gamma_Si_0 / wsink_Si
                1214 #endif
                1215
e0f9a7ba0b Matt*1216 C  POM flux leaving the cell
                1217
                1218           PONflux_l = (PONflux_u+N_spm(i,j,k)*drF(k)
                1219      &           *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
                1220      &           *hFacC(i,j,k,bi,bj))
                1221
                1222           POPflux_l = (POPflux_u+P_spm(i,j,k)*drF(k)
                1223      &           *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
                1224      &           *hFacC(i,j,k,bi,bj))
                1225
                1226 C  CaCO3 flux leaving the cell
                1227
                1228           CaCO3flux_l = (caco3flux_u+CaCO3_uptake(i,j,k)*drF(k)
                1229      &           *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k)
                1230      &           *hFacC(i,j,k,bi,bj))
                1231
00fa2d4ddd mmaz*1232 #ifdef USE_SIBLING
                1233           Siflux_l = (Siflux_u+Si_spm(i,j,k)*drF(k)
                1234      &           *hFacC(i,j,k,bi,bj))/(1+zremin_Si*drF(k)
                1235      &           *hFacC(i,j,k,bi,bj))
                1236 #endif
                1237
a284455135 Matt*1238 C  Check if we are not a bottom cell
e0f9a7ba0b Matt*1239
a284455135 Matt*1240           IF ( k.NE.kLowC(i,j,bi,bj) ) THEN
                1241 C  Start with cells that are not the deepest cells
e0f9a7ba0b Matt*1242
                1243 C  Nutrient accumulation in a cell is given by the biological production
                1244 C  (and instant remineralization) of particulate organic matter
                1245 C  plus flux thought upper interface minus flux through lower interface.
                1246 C  (Since not deepest cell: hFacC=1)
                1247
                1248            N_reminp(i,j,k) = (PONflux_u + N_spm(i,j,k)*drF(k)
                1249      &                    - PONflux_l)*recip_drF(k)
                1250
                1251            P_reminp(i,j,k) = (POPflux_u + P_spm(i,j,k)*drF(k)
                1252      &                    - POPflux_l)*recip_drF(k)
                1253
                1254            CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_uptake(i,j,k)
                1255      &                    *drF(k) - CaCO3flux_l)*recip_drF(k)
                1256
00fa2d4ddd mmaz*1257 #ifdef USE_SIBLING
                1258            Si_reminp(i,j,k) = (Siflux_u + Si_spm(i,j,k)*drF(k)
                1259      &                     - Siflux_l)*recip_drF(k)
                1260 #endif
                1261
e0f9a7ba0b Matt*1262            Fe_sed(i,j,k) = 0. _d 0
                1263
                1264           ELSE
a284455135 Matt*1265 C  Now do bottom layer
e0f9a7ba0b Matt*1266
                1267 C  If this layer is adjacent to bottom topography or it is the deepest
                1268 C  cell of the domain, then remineralize/dissolve in this grid cell
                1269 C  i.e. do not subtract off lower boundary fluxes when calculating remin
                1270
                1271            N_reminp(i,j,k) = PONflux_u*recip_drF(k)
                1272      &                    *recip_hFacC(i,j,k,bi,bj) + N_spm(i,j,k)
                1273
                1274            P_reminp(i,j,k) = POPflux_u*recip_drF(k)
                1275      &                    *recip_hFacC(i,j,k,bi,bj) + P_spm(i,j,k)
                1276
                1277            CaCO3_diss(i,j,k) = CaCO3flux_u*recip_drF(k)
                1278      &                  *recip_hFacC(i,j,k,bi,bj) + CaCO3_uptake(i,j,k)
                1279
00fa2d4ddd mmaz*1280 #ifdef USE_SIBLING
                1281            Si_reminp(i,j,k) = Siflux_u*recip_drF(k)
                1282      &                  *recip_hFacC(i,j,k,bi,bj) + Si_spm(i,j,k)
                1283 #endif
                1284
e0f9a7ba0b Matt*1285 C  Efflux Fed out of sediments
                1286 C  The phosphate flux hitting the bottom boundary
                1287 C  is used to scale the return of iron to the water column.
                1288 C  Maximum value added for numerical stability.
                1289
                1290            POC_sed = PONflux_l * CtoN
                1291
                1292            Fe_sed(i,j,k) = max(epsln, FetoC_sed * POC_sed * recip_drF(k)
                1293      &            *recip_hFacC(i,j,k,bi,bj))
                1294
                1295            log_btm_flx = 1. _d -20
                1296
                1297 CMM: this is causing instability in the adjoint. Needs debugging
                1298 #ifndef BLING_ADJOINT_SAFE
                1299            IF (POC_sed .gt. 0. _d 0) THEN
                1300
                1301 C  Convert from mol N m-2 s-1 to umol C cm-2 d-1 and take the log
                1302
                1303             log_btm_flx = log10(min(43.0 _d 0, POC_sed *
                1304      &           86400. _d 0 * 100.0 _d 0))
                1305
                1306 C  Metamodel gives units of umol C cm-2 d-1, convert to mol N m-2 s-1 and
                1307 C  multiply by no3_2_n to give NO3 consumption rate
                1308
6df47d206d Jean*1309             N_den_benthic(i,j,k) = min (POC_sed * NO3toN / CtoN,
e0f9a7ba0b Matt*1310      &         (10 _d 0)**(-0.9543 _d 0 + 0.7662 _d 0 *
                1311      &         log_btm_flx - 0.235 _d 0 * log_btm_flx * log_btm_flx)
                1312      &         / (CtoN * 86400. _d 0 * 100.0 _d 0) * NO3toN *
                1313      &         PTR_NO3(i,j,k) / (k_no3 + PTR_NO3(i,j,k)) ) *
                1314      &         recip_drF(k)
                1315
6df47d206d Jean*1316            ENDIF
e0f9a7ba0b Matt*1317 #endif
                1318
                1319 C  ---------------------------------------------------------------------
                1320 C  Calculate external bottom fluxes. Positive fluxes
                1321 C  are into the water column from the seafloor. For P, the bottom flux puts
                1322 C  the sinking flux reaching the bottom cell into the water column through
                1323 C  diffusion. For iron, the sinking flux disappears into the sediments if
                1324 C  bottom waters are oxic (assumed adsorbed as oxides). If bottom waters are
                1325 C  anoxic, the sinking flux of Fe is returned to the water column.
                1326 C
                1327 C  For oxygen, the consumption of oxidant required to respire the settling flux
                1328 C  of organic matter (in support of the no3 bottom flux) diffuses from the
                1329 C  bottom water into the sediment.
                1330
                1331 C  Assume all NO3 for benthic denitrification is supplied from the bottom water,
                1332 C  and that all organic N is also consumed under denitrification (Complete
                1333 C  Denitrification, sensu Paulmier, Biogeosciences 2009). Therefore, no NO3 is
                1334 C  regenerated from organic matter respired by benthic denitrification
                1335 C  (necessitating the second term in b_no3).
                1336
6df47d206d Jean*1337            NO3_btmflx(i,j) = PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
e0f9a7ba0b Matt*1338      &                   - N_den_benthic(i,j,k) / NO3toN
                1339
6df47d206d Jean*1340            PO4_btmflx(i,j) = POPflux_l*drF(k)*hFacC(i,j,k,bi,bj)
e0f9a7ba0b Matt*1341
                1342 C  Oxygen flux into sediments is that required to support non-denitrification
                1343 C  respiration, assuming a 4/5 oxidant ratio of O2 to NO3. Oxygen consumption
                1344 C  is allowed to continue at negative oxygen concentrations, representing
                1345 C  sulphate reduction.
                1346
b3fa08b734 Jean*1347            O2_btmflx(i,j) = -( O2toN*PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
e0f9a7ba0b Matt*1348      &                  - N_den_benthic(i,j,k)* 1.25)
                1349
                1350           ENDIF
                1351
a284455135 Matt*1352 C  Calculate free and inorganically associated iron concentration for
                1353 C  scavenging. We assume that there is a spectrum of iron ligands present
                1354 C  in seawater, with  varying binding strengths and whose composition
                1355 C  varies with light and iron concentrations. For example,
                1356 C  photodissocation of ligand complexes occurs under bright light,
                1357 C  weakening the binding strength (e.g. Barbeau et al., Nature 2001),
                1358 C  while at very low iron concentrations (order kfe_eq_lig_femin),
                1359 C  siderophores are thought to be produced as a response to extreme iron
                1360 C  stress. In anoxic waters, iron should be reduced, and therefore mostly
                1361 C  immune to scavenging; therefore, the feprime calculation is not made
                1362 C  if oxygen is less than 0.
e0f9a7ba0b Matt*1363
                1364           kFe_eq_lig = kFe_eq_lig_max-(kFe_eq_lig_max-kFe_eq_lig_min)
                1365      &             *(irr_inst(i,j,k)**2
                1366      &             /(kFe_eq_lig_irr**2+irr_inst(i,j,k)**2))
                1367      &             *max(epsln,min(1. _d 0,(PTR_FE(i,j,k)
                1368      &             -kFe_eq_lig_Femin)/
                1369      &             (PTR_FE(i,j,k)+epsln)*1.2 _d 0))
                1370
                1371           FreeFe = (-(1+kFe_eq_lig*(ligand-PTR_FE(i,j,k)))
                1372      &            +((1+kFe_eq_lig*(ligand-PTR_FE(i,j,k)))**2+4*
                1373      &            kFe_eq_lig*PTR_FE(i,j,k))**(0.5))/(2*
                1374      &            kFe_eq_lig)
                1375
                1376           IF (PTR_O2(i,j,k) .LT. oxic_min)  THEN
                1377            FreeFe = 0. _d 0
                1378           ENDIF
                1379
6df47d206d Jean*1380           Fe_ads_inorg(i,j,k) =
e0f9a7ba0b Matt*1381      &       kFe_inorg*(max(1. _d -8,FreeFe))**(1.5)
                1382
a284455135 Matt*1383 C  Calculate the Fe adsorption using this PONflux_l:
                1384 C  The absolute first order rate constant is calculated from the
                1385 C  concentration of organic particles, after Parekh et al. (2005).
                1386 C  (Galbraith's code limits the value to 1/2dt for numerical stability).
e0f9a7ba0b Matt*1387
6df47d206d Jean*1388           IF ( PONflux_l .GT. 0. _d 0 ) THEN
e0f9a7ba0b Matt*1389             Fe_ads_org(i,j,k) =
                1390      &           kFE_org*(PONflux_l/(epsln + wsink)
                1391      &             * MasstoN)**(0.58)*FreeFe
6df47d206d Jean*1392           ENDIF
e0f9a7ba0b Matt*1393
6df47d206d Jean*1394           PFEflux_l = (PFEflux_u+(Fe_spm(i,j,k)+Fe_ads_inorg(i,j,k)
e0f9a7ba0b Matt*1395      &            +Fe_ads_org(i,j,k))*drF(k)
                1396      &            *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
                1397      &            *hFacC(i,j,k,bi,bj))
                1398
                1399 C  Added the burial flux of sinking particulate iron here as a
                1400 C  diagnostic, needed to calculate mass balance of iron.
                1401 C  this is calculated last for the deepest cell
                1402
6df47d206d Jean*1403           Fe_burial(i,j) = PFEflux_l
e0f9a7ba0b Matt*1404
6df47d206d Jean*1405           IF ( PTR_O2(i,j,k) .LT. oxic_min ) THEN
e0f9a7ba0b Matt*1406             PFEflux_l = 0. _d 0
6df47d206d Jean*1407           ENDIF
e0f9a7ba0b Matt*1408
6df47d206d Jean*1409           Fe_reminp(i,j,k) = (PFEflux_u+(Fe_spm(i,j,k)
e0f9a7ba0b Matt*1410      &            +Fe_ads_inorg(i,j,k)
                1411      &            +Fe_ads_org(i,j,k))*drF(k)
                1412      &            *hFacC(i,j,k,bi,bj)-PFEflux_l)*recip_drF(k)
                1413      &            *recip_hFacC(i,j,k,bi,bj)
                1414
                1415           Fe_reminsum(i,j,k) = Fe_reminp(i,j,k) + Fe_sed(i,j,k)
                1416      &                 - Fe_ads_org(i,j,k) - Fe_ads_inorg(i,j,k)
                1417
                1418 C  Prepare the tracers for the next layer down
                1419
                1420           PONflux_u   = PONflux_l
                1421           POPflux_u   = POPflux_l
                1422           PFEflux_u   = PFEflux_l
                1423           CaCO3flux_u = CaCO3flux_l
00fa2d4ddd mmaz*1424 #ifdef USE_SIBLING
                1425           Siflux_u    = Siflux_l
                1426 #endif
e0f9a7ba0b Matt*1427
                1428          ENDIF
                1429
                1430         ENDDO
                1431        ENDDO
                1432       ENDDO
fe1862e69b Mart*1433 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1434 CADJ STORE N_reminp = comlev1_bibj, key=tkey, kind=isbyte
fe1862e69b Mart*1435 #endif
a284455135 Matt*1436 C-----------------------------------------------------------
                1437 C  remineralization from diel vertical migration
e0f9a7ba0b Matt*1438 #ifdef USE_BLING_DVM
                1439
a284455135 Matt*1440 C  DIEL VERTICAL MIGRATOR EXPORT
                1441 C  The effect of vertically-migrating animals on the export flux of organic
                1442 C  matter from the ocean surface is treated similarly to the scheme of
                1443 C  Bianchi et al., Nature Geoscience 2013.
                1444 C  This involves calculating the stationary depth of vertical migrators, using
                1445 C  an empirical multivariate regression, and ensuring that this remains
                1446 C  above the bottom as well as any suboxic waters.
                1447 C  The total DVM export flux is partitioned between a swimming migratory
                1448 C  component and the stationary component, and these are summed.
e0f9a7ba0b Matt*1449
                1450 C$TAF LOOP = parallel
                1451       DO j=jmin,jmax
                1452 C$TAF LOOP = parallel
                1453        DO i=imin,imax
                1454
a284455135 Matt*1455 C  Initialize
e0f9a7ba0b Matt*1456         o2_upper = 0.
                1457         o2_lower = 0.
                1458         dz_upper = 0.
                1459         dz_lower = 0.
                1460         temp_upper = 0.
                1461         temp_lower = 0.
                1462         z_dvm_regr = 0.
                1463         z_dvm     = 0.
                1464         frac_migr = 0.
                1465         fdvm_migr = 0.
                1466         fdvm_stat = 0.
                1467         fdvmn_vint = 0.
                1468         fdvmp_vint = 0.
                1469         fdvmfe_vint = 0.
                1470
                1471         DO k=1,Nr
                1472
a284455135 Matt*1473          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*1474
a284455135 Matt*1475 C  Calculate the depth of migration based on linear regression.
e0f9a7ba0b Matt*1476
6df47d206d Jean*1477           depth_l=-rF(k+1)
e0f9a7ba0b Matt*1478
a284455135 Matt*1479 C  Average temperature and oxygen over upper 35 m, and 140-515m.
                1480 C  Also convert O2 to mmol m-3.
e0f9a7ba0b Matt*1481
                1482           if ( abs(depth_l) .lt. 35.) then
                1483             dz_upper = dz_upper + drf(k)
                1484             temp_upper = temp_upper + theta(i,j,k,bi,bj)*drf(k)
                1485             o2_upper = o2_upper + PTR_O2(i,j,k) * drf(k)*1.0 _d 3
                1486           endif
                1487           if ( (abs(depth_l) .gt. 140.0 _d 0) .and.
                1488      &          (abs(depth_l) .lt. 515. _d 0)) then
                1489             dz_lower = dz_lower + drf(k)
                1490             temp_lower = temp_lower + theta(i,j,k,bi,bj)*drf(k)
                1491             o2_lower = o2_lower + PTR_O2(i,j,k) * drf(k)*1.0 _d 3
                1492           endif
                1493
                1494          ENDIF
                1495         ENDDO
                1496
                1497         o2_upper = o2_upper / (epsln + dz_upper)
                1498         temp_upper = temp_upper / (epsln + dz_upper)
                1499         o2_lower = o2_lower / (epsln + dz_lower)
                1500         temp_lower = temp_lower / (epsln + dz_lower)
                1501
a284455135 Matt*1502 C  Calculate the regression, using the constants given in Bianchi et al. (2013).
                1503 C  The variable values are bounded to lie within reasonable ranges:
                1504 C         O2 gradient   : [-10,300] mmol/m3
                1505 C         Log10 Chl     : [-1.8,0.85] log10(mg/m3)
                1506 C         mld           : [0,500] m
                1507 C         T gradient    : [-3,20] C
e0f9a7ba0b Matt*1508
                1509         z_dvm_regr = 398. _d 0
                1510      &   - 0.56 _d 0*min(300. _d 0,max(-10. _d 0,(o2_upper - o2_lower)))
                1511      &   - 115. _d 0*min(0.85 _d 0,max(-1.80 _d 0,
                1512      &   log10(max(chl(i,j,1,bi,bj),chl_min))))
                1513      &   + 0.36 _d 0*min(500. _d 0,max(epsln,mld(i,j)))
                1514      &   - 2.40 _d 0*min(20. _d 0,max(-3. _d 0,(temp_upper-temp_lower)))
                1515
a284455135 Matt*1516 C  Limit the depth of migration in polar winter.
                1517 C  Use irr_mem since this is averaged over multiple days, dampening the
                1518 C  diurnal cycle.
                1519 C  Tapers Z_DVM to the minimum when surface irradince is below a given
                1520 C  threshold (here 10 W/m2).
e0f9a7ba0b Matt*1521
                1522         if ( irr_mem(i,j,1,bi,bj) .lt. 10. ) then
                1523           z_dvm_regr = 150. _d 0 + (z_dvm_regr - 150. _d 0) *
                1524      &       irr_mem(i,j,1,bi,bj) / 10. _d 0
                1525         endif
                1526
a284455135 Matt*1527 C  Check for suboxic water within the column. If found, set dvm
                1528 C  stationary depth to 2 layers above it. This is not meant to
                1529 C  represent a cessation of downward migration, but rather the
                1530 C  requirement for aerobic DVM respiration to occur above the suboxic
                1531 C  water, where O2 is available.
e0f9a7ba0b Matt*1532
                1533         tmp_dvm = 0
                1534         DO k=1,Nr-2
                1535
a284455135 Matt*1536          IF ( maskC(i,j,k,bi,bj).EQ.oneRS .AND. tmp_dvm.EQ.0 ) THEN
e0f9a7ba0b Matt*1537
                1538           z_dvm = -rf(k+1)
                1539           if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp_dvm = 1
                1540
                1541          ENDIF
                1542
6df47d206d Jean*1543         ENDDO
e0f9a7ba0b Matt*1544
a284455135 Matt*1545 C  The stationary depth is constrained between 150 and 700, above any
                1546 C  anoxic waters found, and above the bottom.
e0f9a7ba0b Matt*1547
                1548         z_dvm = min(700. _d 0,max(150. _d 0,z_dvm_regr),z_dvm,-rf(k+1))
                1549
a284455135 Matt*1550 C  Calculate the fraction of migratory respiration that occurs
                1551 C  during upwards and downwards swimming. The remainder is
                1552 C  respired near the stationary depth.
                1553 C  Constants for swimming speed and resting time are hard-coded
                1554 C  after Bianchi et al, Nature Geoscience 2013.
e0f9a7ba0b Matt*1555
                1556         frac_migr = max( 0.0 _d 0, min( 1.0 _d 0, (2.0 _d 0 * z_dvm) /
                1557      &        (epsln + 0.05 _d 0 * 0.5 _d 0 * 86400. _d 0)))
                1558
a284455135 Matt*1559 C  Calculate the vertical profile shapes of DVM fluxes.
                1560 C  These are given as the downward organic flux due to migratory
                1561 C  DVM remineralization, defined at the bottom of each layer k.
e0f9a7ba0b Matt*1562
                1563         tmp_dvm = 0
                1564         DO k=1,Nr
                1565
a284455135 Matt*1566          IF ( maskC(i,j,k,bi,bj).EQ.oneRS .AND. tmp_dvm.EQ.0 ) THEN
e0f9a7ba0b Matt*1567
a284455135 Matt*1568 C  First, calculate the part due to active migration above
                1569 C  the stationary depth.
e0f9a7ba0b Matt*1570           if (-rf(k+1) .lt. z_dvm) then
                1571             fdvm_migr = frac_migr / (epsln + z_dvm - (-rf(2))) *
                1572      &            (z_dvm - (-rf(k+1)) )
                1573           else
                1574             fdvm_migr  = 0.0
                1575           endif
                1576
a284455135 Matt*1577 C  Then, calculate the part at the stationary depth.
e0f9a7ba0b Matt*1578
a284455135 Matt*1579 C  Approximation of the complementary error function
                1580 C  From Numerical Recipes (F90, Ch. 6, p. 216)
                1581 C  Returns the complementary error function erfc(x)
                1582 C  with fractional error everywhere less than 1.2e-7
e0f9a7ba0b Matt*1583
6df47d206d Jean*1584           x_erfcc = (-rf(k) - z_dvm) /
e0f9a7ba0b Matt*1585      &        ( (epsln + 2. _d 0 * sigma_dvm**2. _d 0)**0.5)
                1586
6df47d206d Jean*1587           z_erfcc = abs(x_erfcc)
e0f9a7ba0b Matt*1588
6df47d206d Jean*1589           t_erfcc = 1. _d 0/(1. _d 0+0.5 _d 0*z_erfcc)
e0f9a7ba0b Matt*1590
6df47d206d Jean*1591           erfcc = t_erfcc*exp(-z_erfcc*z_erfcc-1.26551223+t_erfcc*
e0f9a7ba0b Matt*1592      &       (1.00002368+t_erfcc*(0.37409196+t_erfcc*
                1593      &       (.09678418+t_erfcc*(-.18628806+t_erfcc*(.27886807+
                1594      &       t_erfcc*(-1.13520398+t_erfcc*(1.48851587+
                1595      &       t_erfcc*(-0.82215223+t_erfcc*0.17087277)))))))))
                1596
                1597           if (x_erfcc .lt. 0.0) then
                1598             erfcc = 2.0 - erfcc
                1599           endif
                1600
                1601           fdvm_stat = (1. _d 0 - frac_migr) / 2. _d 0 * erfcc
                1602
a284455135 Matt*1603 C  Add the shapes, resulting in the 3-d DVM flux operator. If the
                1604 C  current layer is the bottom layer, or the layer beneath the
                1605 C  underlying layer is suboxic, all fluxes at and below the current
                1606 C  layer remain at the initialized value of zero. This will cause all
                1607 C  remaining DVM remineralization to occur in this layer.
e0f9a7ba0b Matt*1608
a284455135 Matt*1609           IF (k.LT.Nr-1) THEN
e0f9a7ba0b Matt*1610             if (PTR_O2(i,j,k+2) .lt. (5. _d 0*oxic_min)) tmp_dvm = 1
                1611           ENDIF
                1612 c!!          if (k .eq. grid_kmt(i,j)) exit
                1613           dvm(i,j,k)  = fdvm_migr + fdvm_stat
                1614
                1615          ENDIF
                1616
6df47d206d Jean*1617         ENDDO
e0f9a7ba0b Matt*1618
a284455135 Matt*1619 C  Sum up the total organic flux to be transported by DVM
e0f9a7ba0b Matt*1620
fe1862e69b Mart*1621         do k = 1, Nr
e0f9a7ba0b Matt*1622           fdvmn_vint  = fdvmn_vint  + N_dvm(i,j,k)  * drf(k)
                1623           fdvmp_vint  = fdvmp_vint  + P_dvm(i,j,k)  * drf(k)
                1624           fdvmfe_vint = fdvmfe_vint + Fe_dvm(i,j,k) * drf(k)
                1625         enddo
                1626
a284455135 Matt*1627 C  Calculate the remineralization terms as the divergence of the flux
e0f9a7ba0b Matt*1628
                1629         N_remindvm(i,j,1)  = fdvmn_vint *  (1 - dvm(i,j,1)) /
                1630      &     (epsln + drf(1))
                1631         P_remindvm(i,j,1)  = fdvmp_vint *  (1 - dvm(i,j,1)) /
                1632      &     (epsln + drf(1))
                1633         Fe_remindvm(i,j,1) = fdvmfe_vint * (1 - dvm(i,j,1)) /
                1634      &     (epsln + drf(1))
                1635
fe1862e69b Mart*1636         do k = 2, Nr
e0f9a7ba0b Matt*1637           N_remindvm(i,j,k)  = fdvmn_vint  *
                1638      &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
                1639           P_remindvm(i,j,k)  = fdvmp_vint  *
                1640      &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
                1641           Fe_remindvm(i,j,k) = fdvmfe_vint *
                1642      &       (dvm(i,j,k-1) - dvm(i,j,k)) / (epsln + drf(k))
                1643         enddo
                1644
                1645        enddo
                1646       enddo
                1647
                1648 #endif
                1649
a284455135 Matt*1650 C-----------------------------------------------------------
e0f9a7ba0b Matt*1651
                1652       DO k=1,Nr
                1653        DO j=jmin,jmax
                1654         DO i=imin,imax
                1655
a284455135 Matt*1656          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
e0f9a7ba0b Matt*1657
a284455135 Matt*1658 C  Dissolved organic matter slow remineralization
e0f9a7ba0b Matt*1659
                1660 #ifdef BLING_NO_NEG
                1661           DON_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DON
                1662      &                       *PTR_DON(i,j,k),0. _d 0)
                1663
                1664           DOP_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DOP
                1665      &                       *PTR_DOP(i,j,k),0. _d 0)
                1666 #else
                1667           DON_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DON
                1668      &                       *PTR_DON(i,j,k)
                1669
                1670           DOP_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DOP
                1671      &                       *PTR_DOP(i,j,k)
                1672 #endif
                1673
a284455135 Matt*1674 C  Pelagic denitrification
                1675 C  If anoxic
e0f9a7ba0b Matt*1676
                1677           IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
                1678            IF (PTR_NO3(i,j,k) .gt. oxic_min) THEN
                1679             N_den_pelag(i,j,k) = max(epsln, (NO3toN *
                1680      &                    ((1. _d 0 - phi_DOM_2d(i,j,bi,bj))
                1681      &                    * (N_reminp(i,j,k)
                1682      &                    + N_remindvm(i,j,k)) + DON_remin(i,j,k) +
                1683      &                    N_recycle(i,j,k))) - N_den_benthic(i,j,k))
                1684            ENDIF
                1685           ENDIF
                1686
a284455135 Matt*1687 C  oxygen production through photosynthesis
e0f9a7ba0b Matt*1688
                1689           O2_prod(i,j,k) = O2toN * N_uptake(i,j,k)
                1690      &                     + (O2toN - 1.25 _d 0) * N_fix(i,j,k)
                1691
a284455135 Matt*1692 C-----------------------------------------------------------
e0f9a7ba0b Matt*1693 C  ADD TERMS
                1694
a284455135 Matt*1695 C  Nutrients
                1696 C  Sum of fast recycling, decay of sinking POM, and decay of DOM,
                1697 C  less uptake, (less denitrification).
e0f9a7ba0b Matt*1698
                1699           G_PO4(i,j,k) = -P_uptake(i,j,k) + P_recycle(i,j,k)
                1700      &                   + (1. _d 0 - phi_DOM_2d(i,j,bi,bj))
                1701      &                   * (P_reminp(i,j,k)
                1702      &                   + P_remindvm(i,j,k)) + DOP_remin(i,j,k)
                1703
                1704           G_NO3(i,j,k) = -N_uptake(i,j,k)
                1705           IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
a284455135 Matt*1706 C  Anoxic
e0f9a7ba0b Matt*1707            G_NO3(i,j,k) = G_NO3(i,j,k)
                1708      &                    - N_den_pelag(i,j,k) - N_den_benthic(i,j,k)
                1709           ELSE
a284455135 Matt*1710 C  Oxic
e0f9a7ba0b Matt*1711            G_NO3(i,j,k) = G_NO3(i,j,k)
                1712      &                   + N_recycle(i,j,k)
                1713      &                   + (1. _d 0 - phi_DOM_2d(i,j,bi,bj))
                1714      &                   * (N_reminp(i,j,k) + N_remindvm(i,j,k))
                1715      &                   + DON_remin(i,j,k)
a284455135 Matt*1716      &                   - N_den_benthic(i,j,k)
e0f9a7ba0b Matt*1717           ENDIF
                1718
a284455135 Matt*1719 C  Iron
                1720 C  remineralization, sediments and adsorption are all bundled into
                1721 C  Fe_reminsum
e0f9a7ba0b Matt*1722
                1723           G_FE(i,j,k) = -Fe_uptake(i,j,k) + Fe_reminsum(i,j,k)
                1724      &                  + Fe_remindvm(i,j,k) + Fe_recycle(i,j,k)
                1725
a284455135 Matt*1726 C  Dissolved Organic Matter
                1727 C  A fraction of POM remineralization goes into dissolved pools.
e0f9a7ba0b Matt*1728
                1729           G_DON(i,j,k) = DON_prod(i,j,k) + phi_DOM_2d(i,j,bi,bj)
                1730      &                   * (N_reminp(i,j,k) + N_remindvm(i,j,k))
                1731      &                   - DON_remin(i,j,k)
                1732
                1733           G_DOP(i,j,k) = DOP_prod(i,j,k) + phi_DOM_2d(i,j,bi,bj)
                1734      &                   * (P_reminp(i,j,k) + P_remindvm(i,j,k))
                1735      &                   - DOP_remin(i,j,k)
                1736
a284455135 Matt*1737 C  Oxygen:
                1738 C  Assuming constant O2:N ratio in terms of oxidant required per mol of
                1739 C  organic N. This implies a constant stoichiometry of C:N and H:N (where H is
                1740 C  reduced, organic H). Because the N provided by N2 fixation is reduced from
                1741 C  N2, rather than NO3-, the o2_2_n_fix is slightly less than the NO3- based
                1742 C  ratio (by 1.25 mol O2/ mol N). Account for the organic matter respired
                1743 C  through benthic denitrification by subtracting 5/4 times the benthic
                1744 C  denitrification NO3 utilization rate from the overall oxygen consumption.
e0f9a7ba0b Matt*1745
                1746           G_O2(i,j,k) = O2_prod(i,j,k)
a284455135 Matt*1747 C  If oxic
e0f9a7ba0b Matt*1748           IF (PTR_O2(i,j,k) .gt. oxic_min) THEN
6df47d206d Jean*1749            G_O2(i,j,k) = G_O2(i,j,k)
e0f9a7ba0b Matt*1750      &                  -O2toN * ((1. _d 0 - phi_DOM_2d(i,j,bi,bj))
                1751      &                  * (N_reminp(i,j,k) + N_remindvm(i,j,k))
                1752      &                  + DON_remin(i,j,k) +  N_recycle(i,j,k))
a284455135 Matt*1753      &                  + N_den_benthic(i,j,k) * 1.25 _d 0
                1754 C  If anoxic but NO3 concentration is very low
                1755 C  (generate negative O2; proxy for HS-).
e0f9a7ba0b Matt*1756           ELSEIF (PTR_NO3(i,j,k) .lt. oxic_min) THEN
6df47d206d Jean*1757            G_O2(i,j,k) = G_O2(i,j,k)
e0f9a7ba0b Matt*1758      &                  -O2toN * ((1. _d 0 - phi_DOM_2d(i,j,bi,bj))
                1759      &                  * (N_reminp(i,j,k) + N_remindvm(i,j,k))
                1760      &                  + DON_remin(i,j,k) +  N_recycle(i,j,k))
                1761           ENDIF
                1762
00fa2d4ddd mmaz*1763 #ifdef USE_SIBLING
                1764           G_Si(i,j,k) = -Si_uptake(i,j,k) + Si_recycle(i,j,k)
                1765      &                  + Si_reminp(i,j,k)
                1766 #endif
                1767
a284455135 Matt*1768 C  Carbon flux diagnostic
e0f9a7ba0b Matt*1769
                1770           NCP(i,j,k) = (N_uptake(i,j,k) + N_fix(i,j,k)
                1771      &                   - N_recycle(i,j,k)
                1772      &                   - (1. _d 0 - phi_DOM_2d(i,j,bi,bj))
                1773      &                   * (N_reminp(i,j,k) + N_remindvm(i,j,k))
                1774      &                   - DON_remin(i,j,k) ) * CtoN
                1775
                1776           G_CaCO3(i,j,k) = CaCO3_diss(i,j,k) - CaCO3_uptake(i,j,k)
                1777
                1778           G_ALK(i,j,k) = - G_NO3(i,j,k)
                1779      &              + 2. _d 0*G_CaCO3(i,j,k)
                1780
                1781           G_DIC(i,j,k) = -NCP(i,j,k) + G_CaCO3(i,j,k)
                1782
a284455135 Matt*1783 C  for diagnostics: convert to mol C/m3
e0f9a7ba0b Matt*1784
6df47d206d Jean*1785           Phy_lg_local(i,j,k) = Phy_lg_local(i,j,k) * CtoN
                1786           Phy_sm_local(i,j,k) = Phy_sm_local(i,j,k) * CtoN
                1787           Phy_diaz_local(i,j,k) = Phy_diaz_local(i,j,k) * CtoN
e0f9a7ba0b Matt*1788
a284455135 Matt*1789 C  for constraints determine POC, assuming that phytoplankton carbon
                1790 C  is 30% of POC
e0f9a7ba0b Matt*1791
6df47d206d Jean*1792           poc(i,j,k,bi,bj) = (Phy_lg_local(i,j,k) + Phy_sm_local(i,j,k)
82e538d851 aver*1793      &                 + Phy_diaz_local(i,j,k) ) * CtoN * 3.33333 _d 0
e0f9a7ba0b Matt*1794
                1795          ENDIF
                1796
                1797         ENDDO
                1798        ENDDO
                1799       ENDDO
                1800
a284455135 Matt*1801 C ---------------------------------------------------------------------
e0f9a7ba0b Matt*1802
                1803 #ifdef ALLOW_DIAGNOSTICS
                1804       IF ( useDiagnostics ) THEN
                1805
82e538d851 aver*1806       DO k=1,Nr
                1807        DO j=jmin,jmax
                1808         DO i=imin,imax
                1809
                1810          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
                1811
                1812           NPP(i,j,k) = (N_uptake(i,j,k) + N_fix(i,j,k)) * CtoN
                1813
                1814 c  for diagnostics: convert biomass to mol C/m3
                1815
                1816           Phy_lg_local(i,j,k) = Phy_lg_local(i,j,k) * CtoN
                1817           Phy_sm_local(i,j,k) = Phy_sm_local(i,j,k) * CtoN
                1818           Phy_diaz_local(i,j,k) = Phy_diaz_local(i,j,k) * CtoN
                1819
                1820 c  Carbon flux diagnostic
                1821
                1822           POC_flux(i,j,k) = CtoN * N_spm(i,j,k)
                1823
                1824 c  C:Chl ratio in g C / g Chl
                1825
                1826           IF ( theta_Fe(i,j,k).EQ.zeroRL ) THEN
                1827              theta_Fe_inv(i,j,k) = 0. _d 0
                1828           ELSE
                1829              theta_Fe_inv(i,j,k) = oneRL/theta_Fe(i,j,k)
                1830           ENDIF
                1831
                1832          ELSE
                1833              NPP(i,j,k) = 0. _d 0
                1834              Phy_lg_local(i,j,k) = 0. _d 0
                1835              Phy_sm_local(i,j,k) = 0. _d 0
                1836              Phy_diaz_local(i,j,k) = 0. _d 0
                1837              POC_flux(i,j,k) = 0. _d 0
                1838              theta_Fe_inv(i,j,k) = 0. _d 0
                1839          ENDIF
                1840
                1841         ENDDO
                1842        ENDDO
                1843       ENDDO
                1844
a284455135 Matt*1845 C 3d global variables
e0f9a7ba0b Matt*1846         CALL DIAGNOSTICS_FILL(chl,    'BLGCHL  ',0,Nr,1,bi,bj,myThid)
                1847         CALL DIAGNOSTICS_FILL(poc,    'BLGPOC  ',0,Nr,1,bi,bj,myThid)
                1848         CALL DIAGNOSTICS_FILL(irr_mem,'BLGIMEM ',0,Nr,1,bi,bj,myThid)
a284455135 Matt*1849 C 3d local variables
e0f9a7ba0b Matt*1850         CALL DIAGNOSTICS_FILL(G_DIC   ,'BLGBIOC ',0,Nr,2,bi,bj,myThid)
                1851         CALL DIAGNOSTICS_FILL(G_ALK   ,'BLGBIOAL',0,Nr,2,bi,bj,myThid)
                1852         CALL DIAGNOSTICS_FILL(G_O2    ,'BLGBIOO2',0,Nr,2,bi,bj,myThid)
                1853         CALL DIAGNOSTICS_FILL(G_Fe    ,'BLGBIOFE',0,Nr,2,bi,bj,myThid)
                1854         CALL DIAGNOSTICS_FILL(G_PO4   ,'BLGBIOP ',0,Nr,2,bi,bj,myThid)
                1855         CALL DIAGNOSTICS_FILL(G_NO3   ,'BLGBION ',0,Nr,2,bi,bj,myThid)
                1856         CALL DIAGNOSTICS_FILL(Phy_sm_local,'BLGPSM  ',0,Nr,2,bi,bj,
                1857      &       myThid)
                1858         CALL DIAGNOSTICS_FILL(Phy_lg_local,'BLGPLG  ',0,Nr,2,bi,bj,
                1859      &       myThid)
                1860         CALL DIAGNOSTICS_FILL(Phy_diaz_local,'BLGPDIA ',0,Nr,2,bi,bj,
                1861      &       myThid)
                1862         CALL DIAGNOSTICS_FILL(irrk,    'BLGIRRK ',0,Nr,2,bi,bj,myThid)
                1863         CALL DIAGNOSTICS_FILL(irr_eff, 'BLGIEFF ',0,Nr,2,bi,bj,myThid)
00fa2d4ddd mmaz*1864         CALL DIAGNOSTICS_FILL(irr_inst,'BLGIRRIS',0,Nr,2,bi,bj,myThid)
e0f9a7ba0b Matt*1865         CALL DIAGNOSTICS_FILL(theta_Fe,'BLGCHL2C',0,Nr,2,bi,bj,myThid)
                1866         CALL DIAGNOSTICS_FILL(theta_Fe_inv,'BLGC2CHL',0,Nr,2,bi,bj,
                1867      &       myThid)
                1868         CALL DIAGNOSTICS_FILL(Fe_lim,  'BLGFELIM',0,Nr,2,bi,bj,myThid)
                1869         CALL DIAGNOSTICS_FILL(NO3_lim, 'BLGNLIM ',0,Nr,2,bi,bj,myThid)
                1870         CALL DIAGNOSTICS_FILL(PO4_lim, 'BLGPLIM ',0,Nr,2,bi,bj,myThid)
                1871 #ifdef USE_SIBLING
00fa2d4ddd mmaz*1872         CALL DIAGNOSTICS_FILL(G_SI    ,'BLGBIOSI',0,Nr,2,bi,bj,myThid)
e0f9a7ba0b Matt*1873         CALL DIAGNOSTICS_FILL(Si_lim,  'BLGSILIM',0,Nr,2,bi,bj,myThid)
00fa2d4ddd mmaz*1874         CALL DIAGNOSTICS_FILL(Si_uptake,'BLGSIUP ',0,Nr,2,bi,bj,myThid)
                1875         CALL DIAGNOSTICS_FILL(Si_recycle,'BLGSIREC',0,Nr,2,bi,bj,myThid)
                1876         CALL DIAGNOSTICS_FILL(Si_reminp,'BLGSIREM',0,Nr,2,bi,bj,myThid)
                1877         CALL DIAGNOSTICS_FILL(SitoN_uptake,'BLGSI2N ',0,Nr,2,bi,bj,
                1878      &       myThid)
                1879         CALL DIAGNOSTICS_FILL(frac_diss_Si,'BLGSIDIS',0,Nr,2,bi,bj,
                1880      &       myThid)
e0f9a7ba0b Matt*1881 #endif
                1882         CALL DIAGNOSTICS_FILL(light_lim,'BLGLLIM ',0,Nr,2,bi,bj,myThid)
                1883         CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid)
                1884         CALL DIAGNOSTICS_FILL(PtoN,    'BLGP2N  ',0,Nr,2,bi,bj,myThid)
                1885         CALL DIAGNOSTICS_FILL(FetoN,   'BLGFE2N ',0,Nr,2,bi,bj,myThid)
                1886         CALL DIAGNOSTICS_FILL(NPP,     'BLGNPP  ',0,Nr,2,bi,bj,myThid)
                1887         CALL DIAGNOSTICS_FILL(NCP,     'BLGNCP  ',0,Nr,2,bi,bj,myThid)
                1888         CALL DIAGNOSTICS_FILL(Fe_dvm,  'BLGFEDVM',0,Nr,2,bi,bj,myThid)
                1889         CALL DIAGNOSTICS_FILL(Fe_spm,  'BLGFESPM',0,Nr,2,bi,bj,myThid)
                1890         CALL DIAGNOSTICS_FILL(Fe_recycle,  'BLGFEREC',0,Nr,2,bi,bj,
                1891      &       myThid)
                1892         CALL DIAGNOSTICS_FILL(Fe_remindvm, 'BLGFERD ',0,Nr,2,bi,bj,
                1893      &       myThid)
                1894         CALL DIAGNOSTICS_FILL(Fe_uptake,'BLGFEUP ',0,Nr,2,bi,bj,myThid)
                1895         CALL DIAGNOSTICS_FILL(N_den_benthic,'BLGNDENB',0,Nr,2,bi,bj,
                1896      &       myThid)
                1897         CALL DIAGNOSTICS_FILL(N_den_pelag, 'BLGNDENP',0,Nr,2,bi,bj,
                1898      &       myThid)
                1899         CALL DIAGNOSTICS_FILL(N_dvm,    'BLGNDVM ',0,Nr,2,bi,bj,myThid)
                1900         CALL DIAGNOSTICS_FILL(N_fix,    'BLGNFIX ',0,Nr,2,bi,bj,myThid)
                1901         CALL DIAGNOSTICS_FILL(DON_prod, 'BLGDONP ',0,Nr,2,bi,bj,myThid)
                1902         CALL DIAGNOSTICS_FILL(DON_remin,'BLGDONR ',0,Nr,2,bi,bj,myThid)
                1903         CALL DIAGNOSTICS_FILL(N_spm,    'BLGNSPM ',0,Nr,2,bi,bj,myThid)
                1904         CALL DIAGNOSTICS_FILL(N_recycle,'BLGNREC ',0,Nr,2,bi,bj,myThid)
                1905         CALL DIAGNOSTICS_FILL(N_remindvm,'BLGNRD  ',0,Nr,2,bi,bj,myThid)
                1906         CALL DIAGNOSTICS_FILL(N_reminp, 'BLGNREM ',0,Nr,2,bi,bj,myThid)
                1907         CALL DIAGNOSTICS_FILL(N_uptake, 'BLGNUP  ',0,Nr,2,bi,bj,myThid)
                1908         CALL DIAGNOSTICS_FILL(P_dvm,    'BLGPDVM ',0,Nr,2,bi,bj,myThid)
                1909         CALL DIAGNOSTICS_FILL(DOP_prod, 'BLGDOPP ',0,Nr,2,bi,bj,myThid)
                1910         CALL DIAGNOSTICS_FILL(DOP_remin,'BLGDOPR ',0,Nr,2,bi,bj,myThid)
                1911         CALL DIAGNOSTICS_FILL(P_spm,    'BLGPSPM ',0,Nr,2,bi,bj,myThid)
                1912         CALL DIAGNOSTICS_FILL(P_recycle,'BLGPREC ',0,Nr,2,bi,bj,myThid)
                1913         CALL DIAGNOSTICS_FILL(P_remindvm,'BLGPRD  ',0,Nr,2,bi,bj,myThid)
                1914         CALL DIAGNOSTICS_FILL(P_reminp, 'BLGPREM ',0,Nr,2,bi,bj,myThid)
                1915         CALL DIAGNOSTICS_FILL(P_uptake, 'BLGPUP  ',0,Nr,2,bi,bj,myThid)
                1916         CALL DIAGNOSTICS_FILL(mu,       'BLGMU   ',0,Nr,2,bi,bj,myThid)
                1917         CALL DIAGNOSTICS_FILL(mu_diaz,  'BLGMUDIA',0,Nr,2,bi,bj,myThid)
                1918         CALL DIAGNOSTICS_FILL(CaCO3_diss,  'BLGCCdis',0,Nr,2,bi,bj,
                1919      &       myThid)
                1920         CALL DIAGNOSTICS_FILL(CaCO3_uptake,'BLGCCpro',0,Nr,2,bi,bj,
                1921      &       myThid)
                1922         CALL DIAGNOSTICS_FILL(Fe_ads_org,  'BLGFEADO',0,Nr,2,bi,bj,
                1923      &       myThid)
                1924         CALL DIAGNOSTICS_FILL(Fe_ads_inorg,'BLGFEADI',0,Nr,2,bi,bj,
                1925      &       myThid)
                1926         CALL DIAGNOSTICS_FILL(Fe_sed,    'BLGFESED',0,Nr,2,bi,bj,myThid)
                1927         CALL DIAGNOSTICS_FILL(Fe_reminp, 'BLGFEREM',0,Nr,2,bi,bj,myThid)
                1928         CALL DIAGNOSTICS_FILL(N_reminp,  'BLGNREM ',0,Nr,2,bi,bj,myThid)
                1929         CALL DIAGNOSTICS_FILL(P_reminp,  'BLGPREM ',0,Nr,2,bi,bj,myThid)
                1930         CALL DIAGNOSTICS_FILL(no3_adj,   'BLGNONEN',0,Nr,2,bi,bj,myThid)
a284455135 Matt*1931 C 2d local variables
82e538d851 aver*1932         CALL DIAGNOSTICS_FILL(mld,       'BLGMLD  ',0,1,2,bi,bj,myThid)
e0f9a7ba0b Matt*1933         CALL DIAGNOSTICS_FILL(Fe_burial, 'BLGFEBUR',0,1,2,bi,bj,myThid)
                1934         CALL DIAGNOSTICS_FILL(NO3_btmflx,'BLGNSED ',0,1,2,bi,bj,myThid)
                1935         CALL DIAGNOSTICS_FILL(PO4_btmflx,'BLGPSED ',0,1,2,bi,bj,myThid)
                1936         CALL DIAGNOSTICS_FILL(O2_btmflx, 'BLGO2SED',0,1,2,bi,bj,myThid)
                1937       ENDIF
                1938 #endif /* ALLOW_DIAGNOSTICS */
                1939
                1940 #endif /* USE_BLING_V1 */
                1941
                1942 #endif /* ALLOW_BLING */
                1943
                1944       RETURN
                1945       END