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