pkg/bling/bling_bio.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(
e0f9a7ba0b Matt*0008      I           PTR_O2, PTR_FE, PTR_PO4, PTR_DOP,
                0009      O           G_DIC, G_ALK, G_O2, G_FE,
                0010      O           G_PO4, G_DOP,
                0011      I           bi, bj, imin, imax, jmin, jmax,
                0012      I           myTime, myIter, myThid)
                0013
                0014 C     =================================================================
                0015 C     | subroutine bling_bio_nitrogen
                0016 C     | o Nutrient uptake and partitioning between organic pools.
                0017 C     | - Phytoplankton biomass-specific growth rate is calculated
                0018 C     |   as a function of light, nutrient limitation, and
                0019 C     |   temperature.
                0020 C     | - A simple relationship between growth rate,
                0021 C     |   biomass, and uptake is derived by assuming that growth is
                0022 C     |   exactly balanced by losses.
                0023 C     | o Organic matter export and remineralization.
                0024 C     | - Calculate the nutrient flux to depth from bio activity
                0025 C     | - Iron source from sediments
                0026 C     | - Iron scavenging
                0027 C     =================================================================
                0028
                0029       IMPLICIT NONE
                0030
                0031 C     === Global variables ===
                0032 C     phyto_sm     :: Small phytoplankton biomass
                0033 C     phyto_lg     :: Large phytoplankton biomass
                0034 C     phyto_diaz   :: Diazotroph phytoplankton biomass
                0035 C *** if ADVECT_PHYTO, these are fraction of total biomass instead ***
                0036
                0037 #include "SIZE.h"
                0038 #include "DYNVARS.h"
                0039 #include "EEPARAMS.h"
                0040 #include "PARAMS.h"
                0041 #include "GRID.h"
                0042 #include "BLING_VARS.h"
                0043 #include "PTRACERS_SIZE.h"
                0044 #include "PTRACERS_PARAMS.h"
a284455135 Matt*0045 #ifdef ALLOW_AUTODIFF_TAMC
e0f9a7ba0b Matt*0046 # include "tamc.h"
                0047 #endif
                0048
                0049 C     === Routine arguments ===
                0050 C     bi,bj         :: tile indices
                0051 C     iMin,iMax     :: computation DOPain: 1rst index range
                0052 C     jMin,jMax     :: computation DOPain: 2nd  index range
                0053 C     myTime        :: current time
                0054 C     myIter        :: current timestep
                0055 C     myThid        :: thread Id. number
                0056       INTEGER bi, bj, imin, imax, jmin, jmax
                0057       _RL     myTime
                0058       INTEGER myIter
                0059       INTEGER myThid
                0060 C     === Input ===
                0061 C     PTR_O2        :: oxygen concentration
                0062 C     PTR_FE        :: iron concentration
                0063 C     PTR_PO4       :: phosphate concentration
                0064 C     PTR_DOP       :: dissolved organic phosphorus concentration
                0065       _RL     PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0066       _RL     PTR_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0067       _RL     PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0068       _RL     PTR_DOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0069 C     === Output ===
                0070 C     G_xxx         :: tendency term for tracer xxx
                0071       _RL     G_DIC     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0072       _RL     G_ALK     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0073       _RL     G_O2      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0074       _RL     G_FE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0075       _RL     G_PO4     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0076       _RL     G_DOP     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0077
                0078 #ifdef ALLOW_BLING
                0079 # ifdef USE_BLING_V1
                0080 C     === Local variables ===
                0081 C     i,j,k         :: loop indices
                0082 C     Phy_lg_local  :: biomass in large phytoplankton
                0083 C     Phy_sm_local  :: biomass in small phytoplankton
                0084 C     PO4_lim       :: phosphate limitation
                0085 C     Fe_lim        :: iron limitation for phytoplankton
                0086 C     FetoP_up      :: iron to phosphorus uptake ratio
                0087 C     light_lim     :: light limitation
                0088 C     expkT         :: temperature function
                0089 C     Pc_m          :: light-saturated max photosynthesis rate for phyto
                0090 C     Pc_tot        :: carbon-specific photosynthesis rate
                0091 C     theta_Fe      :: Chl:C ratio
                0092 C     theta_Fe_inv  :: C:Chl ratio
                0093 C     theta_Fe_max  :: Fe-replete maximum Chl:C ratio
                0094 C     alpha_Fe      :: initial slope of the P-I curve
                0095 C     irrk          :: nut-limited efficiency of algal photosystems
                0096 C     irr_inst      :: instantaneous light
                0097 C     irr_eff       :: effective irradiance
                0098 C     mld           :: mixed layer depth
                0099 C     mu            :: net carbon-specific growth rate for phyt
                0100 C     biomass_sm    :: nutrient concentration in small phyto biomass
                0101 C     biomass_lg    :: nutrient concentration in large phyto biomass
                0102 C     P_uptake      :: PO4 utilization by phytoplankton
                0103 C     Fe_uptake     :: dissolved Fe utilization by phytoplankton
                0104 C     CaCO3_uptake  :: Calcium carbonate uptake for shell formation
                0105 C     CaCO3_diss    :: Calcium carbonate dissolution
                0106 C     G_CaCO3       :: tendency of calcium carbonate
                0107 C     DOP_prod      :: production of dissolved organic phosphorus
                0108 C     DOP_remin     :: remineralization of dissolved organic phosphorus
                0109 C     frac_exp      :: fraction of sinking particulate organic matter
                0110 C     P_spm         :: particulate sinking of phosphorus
                0111 C     Fe_spm        :: particulate sinking of iron
                0112 C     P_recycle     :: recycling of newly-produced organic phosphorus
                0113 C     Fe_recycle    :: recycling of newly-produced organic iron
                0114 C     P_reminp      :: remineralization of particulate organic phosphorus
                0115 C     Fe_reminsum   :: iron remineralization and adsorption
                0116 C     POC_flux      :: particulate organic carbon flux
                0117 C     NPP           :: net primary production
                0118 C     NCP           :: net community production
                0119 C     depth_l       :: depth of lower interface
                0120 C     *flux_u       :: "*" flux through upper interface
                0121 C     *flux_l       :: "*" flux through lower interface
                0122 C     *_export      :: vertically-integrated export of "*"
                0123 C     zremin        :: remineralization lengthscale for nutrients
                0124 C     zremin_caco3  :: remineralization lengthscale for CaCO3
                0125 C     wsink         :: speed of sinking particles
                0126 C     POC_sed       :: flux of particulate organic carbon to sediment
                0127 C     Fe_sed        :: sediment iron efflux
                0128 C     kFe_eq_lig    :: iron ligand stability constant
                0129 C     FreeFe        :: free (unbound) iron concentration
                0130 C     Fe_ads_inorg  :: iron adsorption
                0131 C     Fe_ads_org    :: iron adsorption  (2nd order)
                0132 C
                0133       INTEGER i,j,k
                0134       INTEGER bottomlayer
                0135       _RL Phy_lg_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0136       _RL Phy_sm_local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0137       _RL PO4_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0138       _RL Fe_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0139       _RL FetoP_up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0140       _RL light_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0141       _RL expkT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0142       _RL Pc_m
                0143       _RL Pc_tot
                0144       _RL theta_Fe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0145       _RL theta_Fe_inv(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0146       _RL theta_Fe_max
                0147       _RL alpha_Fe
                0148       _RL irrk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0149       _RL irr_inst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0150       _RL irr_eff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0151       _RL mld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0152       _RL mu(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0153       _RL biomass_sm
                0154       _RL biomass_lg
                0155       _RL P_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0156       _RL Fe_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0157       _RL CaCO3_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0158       _RL CaCO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0159       _RL G_CaCO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0160       _RL DOP_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0161       _RL DOP_remin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0162       _RL frac_exp
                0163       _RL P_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0164       _RL Fe_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0165       _RL P_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0166       _RL Fe_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0167       _RL P_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0168       _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0169       _RL Fe_reminsum(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0170       _RL POC_flux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0171       _RL NPP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0172       _RL NCP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0173       _RL depth_l
                0174       _RL POPflux_u
                0175       _RL POPflux_l
                0176       _RL PFEflux_u
                0177       _RL PFEflux_l
                0178       _RL CaCO3flux_u
                0179       _RL CaCO3flux_l
                0180       _RL POP_export(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0181       _RL PFE_export(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0182       _RL CaCO3_export(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0183       _RL zremin
                0184       _RL zremin_caco3
                0185       _RL wsink
                0186       _RL POC_sed
                0187       _RL Fe_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0188       _RL kFe_eq_lig
                0189       _RL FreeFe
                0190       _RL Fe_ads_inorg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0191       _RL Fe_ads_org(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0192       _RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0193       _RL po4_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0194       _RL fe_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0195       _RL o2_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0196       _RL dop_adj(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0197 #ifdef ML_MEAN_PHYTO
                0198       _RL tmp_p_sm_ML
                0199       _RL tmp_p_lg_ML
                0200       _RL tmp_p_diaz_ML
                0201       _RL tmp_ML
                0202 #endif
4e4ad91a39 Jean*0203 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0204 C     tkey :: tape key (tile dependent)
                0205       INTEGER tkey
4e4ad91a39 Jean*0206 #endif
e0f9a7ba0b Matt*0207 CEOP
                0208
a284455135 Matt*0209 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0210       tkey = bi + (bj - 1)*nSx + (ikey_dynamics - 1)*nSx*nSy
a284455135 Matt*0211 #endif
                0212
e0f9a7ba0b Matt*0213 c ---------------------------------------------------------------------
                0214 c  Initialize output and diagnostics
                0215       DO j=jmin,jmax
                0216        DO i=imin,imax
                0217         mld(i,j)           = 0. _d 0
                0218         POP_export(i,j)    = 0. _d 0
                0219         PFE_export(i,j)    = 0. _d 0
                0220         CaCO3_export(i,j)  = 0. _d 0
                0221        ENDDO
                0222       ENDDO
                0223        DO k=1,Nr
                0224         DO j=jmin,jmax
                0225           DO i=imin,imax
                0226               G_DIC(i,j,k)          = 0. _d 0
                0227               G_ALK(i,j,k)          = 0. _d 0
                0228               G_O2(i,j,k)           = 0. _d 0
                0229               G_FE(i,j,k)           = 0. _d 0
                0230               G_PO4(i,j,k)          = 0. _d 0
                0231               G_DOP(i,j,k)          = 0. _d 0
                0232               G_CaCO3(i,j,k)        = 0. _d 0
                0233               Phy_lg_local(i,j,k)   = 0. _d 0
                0234               Phy_sm_local(i,j,k)   = 0. _d 0
                0235               PO4_lim(i,j,k)        = 0. _d 0
                0236               Fe_lim(i,j,k)         = 0. _d 0
                0237               FetoP_up(i,j,k)       = 0. _d 0
                0238               light_lim(i,j,k)      = 0. _d 0
                0239               expkT(i,j,k)          = 0. _d 0
                0240               theta_Fe(i,j,k)       = 0. _d 0
                0241               theta_Fe_inv(i,j,k)   = 0. _d 0
                0242               irrk(i,j,k)           = 0. _d 0
                0243               irr_inst(i,j,k)       = 0. _d 0
                0244               irr_eff(i,j,k)        = 0. _d 0
                0245               mu(i,j,k)             = 0. _d 0
                0246               P_uptake(i,j,k)       = 0. _d 0
                0247               Fe_uptake(i,j,k)      = 0. _d 0
                0248               CaCO3_uptake(i,j,k)   = 0. _d 0
                0249               CaCO3_diss(i,j,k)     = 0. _d 0
                0250               DOP_prod(i,j,k)       = 0. _d 0
                0251               DOP_remin(i,j,k)      = 0. _d 0
                0252               P_spm(i,j,k)          = 0. _d 0
                0253               Fe_spm(i,j,k)         = 0. _d 0
                0254               P_recycle(i,j,k)      = 0. _d 0
                0255               Fe_recycle(i,j,k)     = 0. _d 0
                0256               P_reminp(i,j,k)       = 0. _d 0
                0257               Fe_reminp(i,j,k)      = 0. _d 0
                0258               Fe_reminsum(i,j,k)    = 0. _d 0
                0259               POC_flux(i,j,k)       = 0. _d 0
                0260               NPP(i,j,k)            = 0. _d 0
                0261               NCP(i,j,k)            = 0. _d 0
                0262               Fe_sed(i,j,k)         = 0. _d 0
                0263               Fe_ads_org(i,j,k)     = 0. _d 0
                0264               Fe_ads_inorg(i,j,k)   = 0. _d 0
                0265           ENDDO
                0266        ENDDO
                0267       ENDDO
                0268
                0269 c-----------------------------------------------------------
                0270 c  avoid negative nutrient concentrations that can result from
                0271 c  advection when low concentrations
                0272
                0273 #ifdef BLING_NO_NEG
                0274       CALL BLING_MIN_VAL(PTR_O2,  1. _d -11, o2_adj,  bi, bj)
                0275       CALL BLING_MIN_VAL(PTR_FE,  1. _d -11, fe_adj,  bi, bj)
                0276       CALL BLING_MIN_VAL(PTR_PO4, 1. _d -8,  po4_adj, bi, bj)
                0277       CALL BLING_MIN_VAL(PTR_DOP, 1. _d -11, dop_adj, bi, bj)
                0278 #endif
                0279
                0280 c-----------------------------------------------------------
                0281 c  Phytoplankton size classes
                0282
                0283        DO k=1,Nr
                0284         DO j=jmin,jmax
                0285          DO i=imin,imax
                0286           Phy_lg_local(i,j,k) = phyto_lg(i,j,k,bi,bj)
                0287           Phy_sm_local(i,j,k) = phyto_sm(i,j,k,bi,bj)
                0288          ENDDO
                0289         ENDDO
                0290        ENDDO
                0291
                0292 c-----------------------------------------------------------
                0293 c  Mixed layer depth calculation for light, phytoplankton and dvm.
                0294 c  Do not need to calculate if not using ML_MEAN_LIGHT, ML_MEAN_PHYTO,
                0295 c  and USE_BLING_DVM
                0296 c  (with BLING_ADJOINT_SAFE flag, USE_BLING_DVM is undefined)
                0297
                0298 #if ( defined (ML_MEAN_LIGHT) || \
                0299       defined (ML_MEAN_PHYTO) || \
                0300       defined (USE_BLING_DVM) )
                0301        CALL BLING_MIXEDLAYER(
                0302      U                         mld,
                0303      I                         bi, bj, imin, imax, jmin, jmax,
                0304      I                         myTime, myIter, myThid)
                0305 #endif
                0306
                0307 c  Phytoplankton mixing
                0308 c  The mixed layer is assumed to homogenize vertical gradients of phytoplankton.
                0309 c  This allows for basic Sverdrup dynamics in a qualitative sense.
                0310 c  This has not been thoroughly tested, and care should be
                0311 c  taken with its interpretation.
                0312
                0313 #ifdef ML_MEAN_PHYTO
                0314       DO j=jmin,jmax
                0315        DO i=imin,imax
                0316
                0317         tmp_p_sm_ML = 0. _d 0
                0318         tmp_p_lg_ML = 0. _d 0
                0319         tmp_ML  = 0. _d 0
                0320
                0321         DO k=1,Nr
                0322
                0323          IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
                0324          IF ((-rf(k+1) .le. mld(i,j)).and.
                0325      &               (-rf(k+1).lt.MLmix_max)) THEN
                0326           tmp_p_sm_ML = tmp_p_sm_ML+Phy_sm_local(i,j,k)*drF(k)
                0327      &                  *hFacC(i,j,k,bi,bj)
                0328           tmp_p_lg_ML = tmp_p_lg_ML+Phy_lg_local(i,j,k)*drF(k)
                0329      &                  *hFacC(i,j,k,bi,bj)
                0330           tmp_ML = tmp_ML + drF(k)
                0331          ENDIF
                0332          ENDIF
                0333
                0334         ENDDO
                0335
                0336         DO k=1,Nr
                0337
                0338          IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
                0339          IF ((-rf(k+1) .le. mld(i,j)).and.
                0340      &               (-rf(k+1).lt.MLmix_max)) THEN
                0341
                0342           Phy_lg_local(i,j,k)   = max(1. _d -8,tmp_p_lg_ML/tmp_ML)
                0343           Phy_sm_local(i,j,k)   = max(1. _d -8,tmp_p_sm_ML/tmp_ML)
                0344
                0345          ENDIF
                0346          ENDIF
                0347
                0348         ENDDO
                0349        ENDDO
                0350       ENDDO
                0351
                0352 #endif
                0353
                0354 c-----------------------------------------------------------
                0355 c  light availability for biological production
                0356
                0357        CALL BLING_LIGHT(
                0358      I                    mld,
                0359      U                    irr_inst, irr_eff,
                0360      I                    bi, bj, imin, imax, jmin, jmax,
                0361      I                    myTime, myIter, myThid )
                0362
                0363 c  phytoplankton photoadaptation to local light level
                0364 c  This represents the fact that phytoplankton cells are adapted
                0365 c  to the averaged light field to which they've been exposed over
                0366 c  their lifetimes, rather than the instantaneous light. The
                0367 c  timescale is set by gamma_irr_mem.
                0368
                0369       DO k=1,Nr
                0370        DO j=jmin,jmax
                0371         DO i=imin,imax
                0372
                0373           irr_mem(i,j,k,bi,bj) = irr_mem(i,j,k,bi,bj) +
                0374      &           (irr_eff(i,j,k) - irr_mem(i,j,k,bi,bj))*
                0375      &           min( 1. _d 0, gamma_irr_mem*PTRACERS_dTLev(k) )
                0376
                0377         ENDDO
                0378        ENDDO
                0379       ENDDO
                0380
                0381 c ---------------------------------------------------------------------
                0382 c  Nutrient uptake and partitioning between organic pools
                0383
                0384 c  Phytoplankton are assumed to grow according to the general properties
                0385 c  described in Geider (1997). This formulation gives a biomass-specific
                0386 c  growthrate as a function of light, nutrient limitation, and
                0387 c  temperature. We modify this relationship slightly here, as described
                0388 c  below, and also use the assumption of steady state growth vs. loss to
                0389 c  derive a simple relationship between growth rate, biomass and uptake.
                0390
                0391       DO k=1,Nr
                0392        DO j=jmin,jmax
                0393         DO i=imin,imax
                0394
                0395          IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
                0396
                0397 c ---------------------------------------------------------------------
                0398 c  First, calculate the limitation terms for PO4 and Fe, and the
                0399 c  Fe-limited Chl:C maximum. The light-saturated maximal photosynthesis
                0400 c  rate term (Pc_m) is simply the product of a prescribed maximal
                0401 c  photosynthesis rate (Pc_0), the Eppley temperature dependence, and a
                0402 c  resource limitation term. The iron limitation term has a lower limit
                0403 c  of Fe_lim_min and is scaled by (k_Fe2P + Fe2P_max) / Fe2P_max so that
                0404 c  it approaches 1 as Fe approaches infinity. Thus, it is of comparable
                0405 c  magnitude to the macro-nutrient limitation term.
                0406
                0407 c  Macro-nutrient limitation
                0408
                0409           PO4_lim(i,j,k) = PTR_PO4(i,j,k)/(PTR_PO4(i,j,k)+k_PO4)
                0410
                0411 c  Iron to macro-nutrient uptake. More iron is utilized relative
                0412 c  to macro-nutrient under iron-replete conditions.
                0413
                0414           FetoP_up(i,j,k) = FetoP_max*PTR_FE(i,j,k)/(k_Fe+PTR_FE(i,j,k))
                0415
                0416 c  Iron limitation
                0417
                0418           Fe_lim(i,j,k) = Fe_lim_min + (1-Fe_lim_min)*
                0419      &                    (FetoP_up(i,j,k)/(k_FetoP + FetoP_up(i,j,k)))*
                0420      &                    (k_FetoP+FetoP_max)/FetoP_max
                0421
                0422 c ---------------------------------------------------------------------
                0423 c  Light-saturated maximal photosynthesis rate
                0424
                0425 c  NB: The temperature effect of Eppley (1972) is used instead of that in
                0426 c  Geider et al (1997) for both simplicity and to incorporate combined
                0427 c  effects on uptake, incorporation into organic matter and photorespiration.
                0428 c  Values of PCc_m are normalized to 0C rather than 20C in Geider et al. (1997)
                0429
                0430           expkT(i,j,k) = exp(kappa_eppley * theta(i,j,k,bi,bj))
                0431
                0432 c  For the effective resource limitation, there is an option to replace
                0433 c  the default Liebig limitation (the minimum of Michaelis-Menton
                0434 c  PO4-limitation, or iron-limitation) by the product (safer for adjoint)
                0435
                0436 #ifdef MULT_NUT_LIM
                0437           Pc_m = Pc_0*expkT(i,j,k)
                0438      &           *PO4_lim(i,j,k)*Fe_lim(i,j,k)*maskC(i,j,k,bi,bj)
                0439 #else
                0440           Pc_m = Pc_0*expkT(i,j,k)
                0441      &           *min(PO4_lim(i,j,k), Fe_lim(i,j,k))*maskC(i,j,k,bi,bj)
                0442 #endif
                0443
                0444 c ---------------------------------------------------------------------
                0445 c  Fe limitation 1) reduces photosynthetic efficiency (alpha_Fe)
                0446 c  and 2) reduces the maximum achievable Chl:C ratio (theta_Fe)
                0447 c  below a prescribed, Fe-replete maximum value (theta_Fe_max),
                0448 c  to approach a prescribed minimum Chl:C (theta_Fe_min) under extreme
                0449 c  Fe-limitation.
                0450
                0451           alpha_Fe = alpha_min + (alpha_max-alpha_min)*Fe_lim(i,j,k)
                0452
                0453           theta_Fe_max = theta_Fe_max_lo+
                0454      &                  (theta_Fe_max_hi-theta_Fe_max_lo)*Fe_lim(i,j,k)
                0455
                0456           theta_Fe(i,j,k) = theta_Fe_max/(1. _d 0 + alpha_Fe
                0457      &                  *theta_Fe_max*irr_mem(i,j,k,bi,bj)/
                0458      &                  (2. _d 0*Pc_m))
                0459
                0460 c  for diagnostics: C:Chl ratio in g C / g Chl
                0461
                0462           IF ( theta_Fe(i,j,k) .EQ.0. ) THEN
                0463            theta_Fe_inv(i,j,k) = 0.
                0464           ELSE
                0465            theta_Fe_inv(i,j,k) = 1./theta_Fe(i,j,k)
                0466           ENDIF
                0467
                0468 c ---------------------------------------------------------------------
                0469 c  Nutrient-limited efficiency of algal photosystems, irrk, is calculated
                0470 c  with the iron limitation term included as a multiplier of the
                0471 c  theta_Fe_max to represent the importance of Fe in forming chlorophyll
                0472 c  accessory antennae, which do not affect the Chl:C but still affect the
                0473 c  phytoplankton ability to use light (eg Stzrepek & Harrison, Nature 2004).
                0474
                0475           irrk(i,j,k) = Pc_m/(epsln + alpha_Fe*theta_Fe_max) +
                0476      &              irr_mem(i,j,k,bi,bj)/2. _d 0
                0477
                0478           light_lim(i,j,k) = ( 1. _d 0 - exp(-irr_eff(i,j,k)
                0479      &               /(epsln + irrk(i,j,k))))
                0480
                0481 c  Carbon-specific photosynthesis rate
                0482
                0483           Pc_tot = Pc_m * light_lim(i,j,k)
                0484
                0485 c ---------------------------------------------------------------------
                0486 c  Account for the maintenance effort that phytoplankton must exert in
                0487 c  order to combat decay. This is prescribed as a fraction of the
                0488 c  light-saturated photosynthesis rate, resp_frac. The result of this
                0489 c  is to set a level of energy availability below which net growth
                0490 c  (and therefore nutrient uptake) is zero, given by resp_frac * Pc_m.
                0491
                0492           mu(i,j,k) = max(0. _d 0, Pc_tot - resp_frac*Pc_m)
                0493
                0494 c ---------------------------------------------------------------------
                0495 c  In order to convert this net carbon-specific growth rate to nutrient
                0496 c  uptake rates, the quantities we are interested in, a biomass is required.
                0497 c  This is determined by the balance between growth and grazing.
                0498
                0499 c  Since there is no explicit biomass tracer, use the result of Dunne
                0500 c  et al. (GBC, 2005) to calculate an implicit biomass from the uptake
                0501 c  rate through the application of a simple idealized grazing law.
                0502
                0503 c  instantaneous nutrient concentration in phyto biomass
                0504
                0505           biomass_lg = pivotal*(mu(i,j,k)/(lambda_0
                0506      &                *expkT(i,j,k)))**3
                0507
                0508           biomass_sm = pivotal*(mu(i,j,k)/(lambda_0
                0509      &                *expkT(i,j,k)))
                0510
                0511 c  phytoplankton biomass diagnostic
                0512 c  for no lag: set gamma_biomass to 0
                0513
                0514           phy_sm_local(i,j,k) = phy_sm_local(i,j,k) +
                0515      &       (biomass_sm - phy_sm_local(i,j,k))
                0516      &       *min(1., gamma_biomass*PTRACERS_dTLev(k))
                0517
                0518           phy_lg_local(i,j,k) = phy_lg_local(i,j,k) +
                0519      &       (biomass_lg - phy_lg_local(i,j,k))
                0520      &       *min(1., gamma_biomass*PTRACERS_dTLev(k))
                0521
                0522 c  use the diagnostic biomass to calculate the chl concentration
                0523 c  in mg/m3 (carbon = 12.01 g/mol)
                0524
                0525           chl(i,j,k,bi,bj) = max(chl_min, CtoP * 12.01 * 1. _d 3 *
                0526      &           theta_Fe(i,j,k) *
                0527      &           (Phy_lg_local(i,j,k) + Phy_sm_local(i,j,k)))
                0528
                0529 c  Nutrient uptake
                0530
                0531           P_uptake(i,j,k) = mu(i,j,k)*(phy_sm_local(i,j,k)
                0532      &                        + phy_lg_local(i,j,k))
                0533
                0534 c  Iron is then taken up as a function of nutrient uptake and iron
                0535 c  limitation, with a maximum Fe:P uptake ratio of Fe2p_max
                0536
                0537           Fe_uptake(i,j,k) = P_uptake(i,j,k)*FetoP_up(i,j,k)
                0538
                0539          ENDIF
                0540         ENDDO
                0541        ENDDO
                0542       ENDDO
                0543
                0544 c  Separate loop for adjoint stores
a284455135 Matt*0545 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0546 CADJ STORE Phy_sm_local   = comlev1_bibj, key=tkey, kind=isbyte
                0547 CADJ STORE Phy_lg_local   = comlev1_bibj, key=tkey, kind=isbyte
a284455135 Matt*0548 #endif
e0f9a7ba0b Matt*0549
                0550       DO k=1,Nr
                0551        DO j=jmin,jmax
                0552         DO i=imin,imax
                0553
                0554          IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
                0555
                0556 c  update biomass
                0557           phyto_lg(i,j,k,bi,bj) = Phy_lg_local(i,j,k)
                0558           phyto_sm(i,j,k,bi,bj) = Phy_sm_local(i,j,k)
                0559
                0560          ENDIF
                0561         ENDDO
                0562        ENDDO
                0563       ENDDO
                0564
                0565 c ---------------------------------------------------------------------
                0566 c  Partitioning between organic pools
                0567
                0568 c  The uptake of nutrients is assumed to contribute to the growth of
                0569 c  phytoplankton, which subsequently die and are consumed by heterotrophs.
                0570 c  This can involve the transfer of nutrient elements between many
                0571 c  organic pools, both particulate and dissolved, with complex histories.
                0572 c  We take a simple approach here, partitioning the total uptake into two
                0573 c  fractions - sinking and non-sinking - as a function of temperature,
                0574 c  following Dunne et al. (2005).
                0575 c  Then, the non-sinking fraction is further subdivided, such that the
                0576 c  majority is recycled instantaneously to the inorganic nutrient pool,
                0577 c  representing the fast turnover of labile dissolved organic matter via
                0578 c  the microbial loop, and the remainder is converted to semi-labile
                0579 c  dissolved organic matter. Iron and macro-nutrient are treated
                0580 c  identically for the first step, but all iron is recycled
                0581 c  instantaneously in the second step (i.e. there is no dissolved organic
                0582 c  iron pool).
                0583
                0584       DO k=1,Nr
                0585        DO j=jmin,jmax
                0586         DO i=imin,imax
                0587
                0588          IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
                0589
                0590 c  sinking particulate organic matter
                0591
                0592           frac_exp = (phi_sm + phi_lg *
                0593      &               (mu(i,j,k)/(lambda_0*expkT(i,j,k)))**2.)/
                0594      &               (1. + (mu(i,j,k)/(lambda_0*expkT(i,j,k)))**2.)*
                0595      &               exp(kappa_remin * theta(i,j,k,bi,bj))
                0596
                0597           P_spm(i,j,k) = frac_exp * P_uptake(i,j,k)
                0598
                0599           Fe_spm(i,j,k) = P_spm(i,j,k)*FetoP_up(i,j,k)
                0600
                0601 c  the remainder is divided between instantaneously recycled and
                0602 c  long-lived dissolved organic matter.
                0603 c  (recycling = P_uptake - P_spm - DOP_prod)
                0604
                0605           DOP_prod(i,j,k) = phi_DOM*(P_uptake(i,j,k)
                0606      &                      - P_spm(i,j,k))
                0607
                0608           P_recycle(i,j,k) = P_uptake(i,j,k) - P_spm(i,j,k)
                0609      &                        - DOP_prod(i,j,k)
                0610
                0611           Fe_recycle(i,j,k) = Fe_uptake(i,j,k) - Fe_spm(i,j,k)
                0612
                0613 c  Carbon flux diagnostic
                0614
                0615           POC_flux(i,j,k) = CtoP*P_spm(i,j,k)
                0616
                0617 c ---------------------------------------------------------------------
                0618 c  Calcium carbonate production
                0619
                0620 c  Alkalinity is consumed through the production of CaCO3. Here, this is
                0621 c  simply a linear function of the implied growth rate of small
                0622 c  phytoplankton, which gave a reasonably good fit to the global
                0623 c  observational synthesis of Dunne (2009). This is consistent
                0624 c  with the findings of Jin et al. (GBC,2006).
                0625
                0626           CaCO3_uptake(i,j,k) = phy_sm_local(i,j,k)*phi_sm
                0627      &           *expkT(i,j,k)*mu(i,j,k)*CatoP
                0628
                0629          ENDIF
                0630         ENDDO
                0631        ENDDO
                0632       ENDDO
                0633
                0634 c ---------------------------------------------------------------------
                0635 c  Nutrients export/remineralization, CaCO3 export/dissolution
                0636 c
                0637 c  The flux at the bottom of a grid cell equals
                0638 C  Fb = (Ft + prod*dz) / (1 + zremin*dz)
                0639 C  where Ft is the flux at the top, and prod*dz is the integrated
                0640 C  production of new sinking particles within the layer.
                0641 C  Ft = 0 in the first layer.
                0642
                0643 C$TAF LOOP = parallel
                0644        DO j=jmin,jmax
                0645 C$TAF LOOP = parallel
                0646         DO i=imin,imax
                0647
                0648 C  Initialize upper flux
                0649
                0650         POPflux_u            = 0. _d 0
                0651         PFEflux_u            = 0. _d 0
                0652         CaCO3flux_u          = 0. _d 0
                0653
                0654         DO k=1,Nr
                0655
                0656 C Initialization here helps taf
                0657
                0658          Fe_ads_org(i,j,k)    = 0. _d 0
                0659
                0660 c  check if we are at a bottom cell
                0661
                0662          bottomlayer = 1
                0663           IF (k.LT.Nr) THEN
                0664            IF (hFacC(i,j,k+1,bi,bj).GT.0) THEN
                0665 c  not a bottom cell
                0666             bottomlayer = 0
                0667            ENDIF
                0668           ENDIF
                0669
                0670          IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN
                0671
                0672 C  Sinking speed is evaluated at the bottom of the cell
                0673
                0674           depth_l=-rF(k+1)
                0675           IF (depth_l .LE. wsink0z)  THEN
                0676            wsink = wsink0_2d(i,j,bi,bj)
                0677           ELSE
                0678            wsink = wsinkacc * (depth_l - wsink0z) + wsink0_2d(i,j,bi,bj)
                0679           ENDIF
                0680
                0681 C  Nutrient remineralization lengthscale
                0682 C  Not an e-folding scale: this term increases with remineralization.
                0683
                0684           zremin = gamma_POM_2d(i,j,bi,bj) * ( PTR_O2(i,j,k)**2 /
                0685      &               (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)
                0686      &               + remin_min )/(wsink + epsln)
                0687
                0688 C  Calcium remineralization relaxed toward the inverse of the
                0689 C  ca_remin_depth constant value as the calcite saturation approaches 0.
                0690
                0691           zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0 - min(1. _d 0,
                0692      &               omegaC(i,j,k,bi,bj) + epsln ))
                0693
                0694 C  POM flux leaving the cell
                0695
                0696           POPflux_l = (POPflux_u+P_spm(i,j,k)*drF(k)
                0697      &           *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
                0698      &           *hFacC(i,j,k,bi,bj))
                0699
                0700 C  CaCO3 flux leaving the cell
                0701
                0702           CaCO3flux_l = (caco3flux_u+CaCO3_uptake(i,j,k)*drF(k)
                0703      &           *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k)
                0704      &           *hFacC(i,j,k,bi,bj))
                0705
                0706 C  Begin iron uptake calculations by determining ligand bound and free iron.
                0707 C  Both forms are available for biology, but only free iron is scavenged
                0708 C  onto particles and forms colloids.
                0709
                0710           kFe_eq_lig = kFe_eq_lig_max-(kFe_eq_lig_max-kFe_eq_lig_min)
                0711      &             *(irr_inst(i,j,k)**2
                0712      &             /(kFe_eq_lig_irr**2+irr_inst(i,j,k)**2))
                0713      &             *max(epsln,min(1. _d 0,(PTR_FE(i,j,k)
                0714      &             -kFe_eq_lig_Femin)/
                0715      &             (PTR_FE(i,j,k)+epsln)*1.2 _d 0))
                0716
                0717 C  Use the quadratic equation to solve for binding between iron and ligands
                0718
                0719           FreeFe = (-(1+kFe_eq_lig*(ligand-PTR_FE(i,j,k)))
                0720      &            +((1+kFe_eq_lig*(ligand-PTR_FE(i,j,k)))**2+4*
                0721      &            kFe_eq_lig*PTR_FE(i,j,k))**(0.5))/(2*
                0722      &            kFe_eq_lig)
                0723
                0724 C  Iron scavenging does not occur in anoxic water (Fe2+ is soluble), so set
                0725 C  FreeFe = 0 when anoxic.  FreeFe should be interpreted the free iron that
                0726 C  participates in scavenging.
                0727
                0728           IF (PTR_O2(i,j,k) .LT. oxic_min)  THEN
                0729            FreeFe = 0. _d 0
                0730           ENDIF
                0731
                0732 C  Two mechanisms for iron uptake, in addition to biological production:
                0733 C  colloidal scavenging and scavenging by organic matter.
                0734
                0735            Fe_ads_inorg(i,j,k) =
                0736      &       kFe_inorg*(max(1. _d -8,FreeFe))**(1.5)
                0737
                0738 C  Scavenging of iron by organic matter:
                0739 C  The POM value used is the bottom boundary flux. This does not occur in
                0740 C  oxic waters, but FreeFe is set to 0 in such waters earlier.
                0741
                0742            IF ( POPflux_l .GT. 0. _d 0 ) THEN
                0743             Fe_ads_org(i,j,k) =
                0744      &           kFE_org*(POPflux_l/(epsln + wsink)
                0745      &             * MasstoN*NtoP)**(0.58)*FreeFe
                0746            ENDIF
                0747
                0748 C  If water is oxic then the iron is remineralized normally. Otherwise
                0749 C  it is completely remineralized (fe 2+ is soluble, but unstable
                0750 C  in oxidizing environments).
                0751
                0752            IF ( PTR_O2(i,j,k) .LT. oxic_min ) THEN
                0753             PFEflux_l = 0. _d 0
                0754            ELSE
                0755             PFEflux_l = (PFEflux_u+(Fe_spm(i,j,k)+Fe_ads_inorg(i,j,k)
                0756      &            +Fe_ads_org(i,j,k))*drF(k)
                0757      &            *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
                0758      &            *hFacC(i,j,k,bi,bj))
                0759            ENDIF
                0760
                0761 C  Nutrient accumulation in a cell is given by the biological production
                0762 C  (and instant remineralization) of particulate organic matter
                0763 C  plus flux thought upper interface minus flux through lower interface.
                0764 C  If this layer is adjacent to bottom topography or it is the deepest
                0765 C  cell of the domain, then remineralize/dissolve in this grid cell
                0766 C  i.e. do not subtract off lower boundary fluxes when calculating remin
                0767
                0768 C  For the deepest cells:
                0769
                0770           IF (bottomlayer.EQ.1) THEN
                0771
                0772            POPflux_l   = 0. _d 0
                0773            CACO3flux_l = 0. _d 0
                0774
                0775 C  Efflux Fed out of sediments
                0776 C  The phosphate flux hitting the bottom boundary
                0777 C  is used to scale the return of iron to the water column.
                0778 C  Maximum value added for numerical stability.
                0779
                0780            POC_sed = POPflux_l * CtoP
                0781
                0782            Fe_sed(i,j,k) = max(epsln, FetoC_sed * POC_sed * recip_drF(k)
                0783      &                   * recip_hFacC(i,j,k,bi,bj))
                0784
                0785           ELSE
                0786
                0787            Fe_sed(i,j,k) = 0. _d 0
                0788
                0789           ENDIF
                0790
                0791           P_reminp(i,j,k) = (POPflux_u + P_spm(i,j,k) * drF(k)
                0792      &                      * hFacC(i,j,k,bi,bj) - POPflux_l)
                0793      &                      * recip_drF(k) * recip_hFacC(i,j,k,bi,bj)
                0794
                0795           CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_uptake(i,j,k)
                0796      &                      * drF(k) * hFacC(i,j,k,bi,bj) - CaCO3flux_l)
                0797      &                      * recip_drF(k) * recip_hFacC(i,j,k,bi,bj)
                0798
                0799           Fe_sed(i,j,k) = 0. _d 0
                0800
                0801           Fe_reminp(i,j,k) = (PFEflux_u + (Fe_spm(i,j,k)
                0802      &                     + Fe_ads_inorg(i,j,k) + Fe_ads_org(i,j,k))
                0803      &                     * drF(k) * hFacC(i,j,k,bi,bj) - PFEflux_l)
                0804      &                     * recip_drF(k) * recip_hFacC(i,j,k,bi,bj)
                0805
                0806           Fe_reminsum(i,j,k) = Fe_reminp(i,j,k) + Fe_sed(i,j,k)
                0807      &                       - Fe_ads_org(i,j,k) - Fe_ads_inorg(i,j,k)
                0808
                0809 C  Added the burial flux of sinking particulate iron here as a
                0810 C  diagnostic, needed to calculate mass balance of iron.
                0811 C  this is calculated last for the deepest cell
                0812
                0813            Fe_burial(i,j) = PFEflux_l
                0814
                0815 C  Prepare the tracers for the next layer down
                0816
                0817            POPflux_u   = POPflux_l
                0818            PFEflux_u   = PFEflux_l
                0819            CaCO3flux_u = CaCO3flux_l
                0820
                0821          ENDIF
                0822
                0823         ENDDO
                0824        ENDDO
                0825       ENDDO
                0826
                0827 C-----------------------------------------------------------
                0828 C  add all tendencies
                0829
                0830        DO k=1,Nr
                0831          DO j=jmin,jmax
                0832           DO i=imin,imax
                0833
                0834 C  Dissolved organic matter slow remineralization
                0835
                0836 #ifdef BLING_NO_NEG
                0837            DOP_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DOP
                0838      &                    *PTR_DOP(i,j,k),0. _d 0)
                0839 #else
                0840            DOP_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DOP
                0841      &                    *PTR_DOP(i,j,k)
                0842 #endif
                0843
                0844 c  Tendencies
                0845
                0846            G_PO4(i,j,k) = -P_uptake(i,j,k) + P_recycle(i,j,k)
                0847      &                    + DOP_remin(i,j,k)
                0848      &                    + (1-phi_DOM) * P_reminp(i,j,k)
                0849
                0850            G_DOP(i,j,k) = DOP_prod(i,j,k) - DOP_remin(i,j,k)
                0851      &                   + phi_DOM * P_reminp(i,j,k)
                0852
                0853            if ( PTR_O2(i,j,k) .GT. oxic_min ) then
                0854              G_O2(i,j,k) = -O2toP*G_PO4(i,j,k)
                0855            else
                0856              G_O2(i,j,k) = 0. _d 0
                0857            endif
                0858
                0859            G_FE(i,j,k) = - Fe_uptake(i,j,k) + Fe_reminsum(i,j,k)
                0860      &                   + Fe_recycle(i,j,k)
                0861
                0862 C  Carbon system diagnostics
                0863 C  Change in DIC from primary production, from recycling and
                0864 C  remineralization, change in carbonate ions concentration
                0865 C  from biological activity:
                0866
                0867            G_CaCO3(i,j,k) = CaCO3_diss(i,j,k) - CaCO3_uptake(i,j,k)
                0868
                0869            NPP(i,j,k) = P_uptake(i,j,k) * CtoP
                0870
                0871            NCP(i,j,k) = -G_PO4(i,j,k)*CtoP
                0872
                0873            G_ALK(i,j,k) = 2. _d 0*G_CaCO3(i,j,k) - NtoP*G_PO4(i,j,k)
                0874
                0875            G_DIC(i,j,k) = -NCP(i,j,k) + G_CaCO3(i,j,k)
                0876
                0877 c  Carbon flux diagnostic
                0878
                0879            POC_flux(i,j,k) = CtoP * P_spm(i,j,k)
                0880
                0881 c  for diagnostics: convert to mol C/m3
                0882
                0883            Phy_lg_local(i,j,k) = Phy_lg_local(i,j,k) * CtoP
                0884            Phy_sm_local(i,j,k) = Phy_sm_local(i,j,k) * CtoP
                0885
                0886 c  for constraints determine POC, assuming that phytoplankton carbon
                0887 c  is 30% of POC
                0888
                0889            poc(i,j,k,bi,bj) = (Phy_lg_local(i,j,k) +
                0890      &                        Phy_sm_local(i,j,k)) * 3.33333 _d 0
                0891
                0892           ENDDO
                0893          ENDDO
                0894        ENDDO
                0895
                0896 c ---------------------------------------------------------------------
                0897
                0898 #ifdef ALLOW_DIAGNOSTICS
                0899       IF ( useDiagnostics ) THEN
                0900
                0901 c 3d global variables
                0902         CALL DIAGNOSTICS_FILL(chl,    'BLGCHL  ',0,Nr,1,bi,bj,myThid)
                0903         CALL DIAGNOSTICS_FILL(poc,    'BLGPOC  ',0,Nr,1,bi,bj,myThid)
                0904         CALL DIAGNOSTICS_FILL(irr_mem,'BLGIMEM ',0,Nr,1,bi,bj,myThid)
                0905 c 3d local variables
                0906         CALL DIAGNOSTICS_FILL(G_DIC   ,'BLGBIOC ',0,Nr,2,bi,bj,myThid)
                0907         CALL DIAGNOSTICS_FILL(G_ALK   ,'BLGBIOAL',0,Nr,2,bi,bj,myThid)
                0908         CALL DIAGNOSTICS_FILL(G_O2    ,'BLGBIOO2',0,Nr,2,bi,bj,myThid)
                0909         CALL DIAGNOSTICS_FILL(G_Fe    ,'BLGBIOFE',0,Nr,2,bi,bj,myThid)
                0910         CALL DIAGNOSTICS_FILL(G_PO4   ,'BLGBIOP ',0,Nr,2,bi,bj,myThid)
                0911         CALL DIAGNOSTICS_FILL(Phy_sm_local,'BLGPSM  ',0,Nr,2,bi,bj,
                0912      &       myThid)
                0913         CALL DIAGNOSTICS_FILL(Phy_lg_local,'BLGPLG  ',0,Nr,2,bi,bj,
                0914      &       myThid)
                0915         CALL DIAGNOSTICS_FILL(irrk,    'BLGIRRK ',0,Nr,2,bi,bj,myThid)
                0916         CALL DIAGNOSTICS_FILL(irr_eff, 'BLGIEFF ',0,Nr,2,bi,bj,myThid)
                0917         CALL DIAGNOSTICS_FILL(theta_Fe,'BLGCHL2C',0,Nr,2,bi,bj,myThid)
                0918         CALL DIAGNOSTICS_FILL(theta_Fe_inv,'BLGC2CHL',0,Nr,2,bi,bj,
                0919      &       myThid)
                0920         CALL DIAGNOSTICS_FILL(Fe_lim,  'BLGFELIM',0,Nr,2,bi,bj,myThid)
                0921         CALL DIAGNOSTICS_FILL(PO4_lim, 'BLGPLIM ',0,Nr,2,bi,bj,myThid)
                0922         CALL DIAGNOSTICS_FILL(light_lim,'BLGLLIM ',0,Nr,2,bi,bj,myThid)
                0923         CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid)
                0924         CALL DIAGNOSTICS_FILL(NPP,     'BLGNPP  ',0,Nr,2,bi,bj,myThid)
                0925         CALL DIAGNOSTICS_FILL(NCP,     'BLGNCP  ',0,Nr,2,bi,bj,myThid)
                0926         CALL DIAGNOSTICS_FILL(Fe_spm,  'BLGFESPM',0,Nr,2,bi,bj,myThid)
                0927         CALL DIAGNOSTICS_FILL(Fe_recycle,'BLGFEREC',0,Nr,2,bi,bj,
                0928      &       myThid)
                0929         CALL DIAGNOSTICS_FILL(Fe_uptake,'BLGFEUP ',0,Nr,2,bi,bj,myThid)
                0930         CALL DIAGNOSTICS_FILL(DOP_prod, 'BLGDOPP ',0,Nr,2,bi,bj,myThid)
                0931         CALL DIAGNOSTICS_FILL(DOP_remin,'BLGDOPR ',0,Nr,2,bi,bj,myThid)
                0932         CALL DIAGNOSTICS_FILL(P_spm,    'BLGPSPM ',0,Nr,2,bi,bj,myThid)
                0933         CALL DIAGNOSTICS_FILL(P_recycle,'BLGPREC ',0,Nr,2,bi,bj,myThid)
                0934         CALL DIAGNOSTICS_FILL(P_reminp, 'BLGPREM ',0,Nr,2,bi,bj,myThid)
                0935         CALL DIAGNOSTICS_FILL(P_uptake, 'BLGPUP  ',0,Nr,2,bi,bj,myThid)
                0936         CALL DIAGNOSTICS_FILL(mu,       'BLGMU   ',0,Nr,2,bi,bj,myThid)
                0937         CALL DIAGNOSTICS_FILL(CaCO3_diss,'BLGCCdis',0,Nr,2,bi,bj,
                0938      &       myThid)
                0939         CALL DIAGNOSTICS_FILL(CaCO3_uptake,'BLGCCpro',0,Nr,2,bi,bj,
                0940      &       myThid)
                0941 c 2d local variables
                0942         CALL DIAGNOSTICS_FILL(mld,       'BLGMLD  ',0,1,1,bi,bj,myThid)
                0943       ENDIF
                0944 #endif /* ALLOW_DIAGNOSTICS */
                0945
                0946 #endif /* USE_BLING_V1 */
                0947
                0948 #endif /* ALLOW_BLING */
                0949
                0950       RETURN
                0951       END