pkg/darwin/darwin_plankton.F

8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: DARWIN_PLANKTON
                0005 C !INTERFACE: ==========================================================
                0006       SUBROUTINE DARWIN_PLANKTON(
                0007      I     Ptr,
                0008      U     gTr,
5e7acb36b1 daat*0009      O     chlout, diags, sumprey,
3f05298721 Oliv*0010      I     PAR, photoTempFunc, hetTempFunc, grazTempFunc, reminTempFunc,
a092808e6b shlo*0011      I     mortTempFunc, mort2TempFunc,
                0012      I     uptakeTempFunc, macromolTempFunc,
ba0b6d5d33 Oliv*0013      I     omegaCl,
6559acf5ff Oliv*0014      I     k_debug,myTime,myIter,myThid)
8fbfd1f382 Oliv*0015
                0016 C !DESCRIPTION:
                0017
                0018 C !USES: ===============================================================
                0019       IMPLICIT NONE
                0020 #ifdef ALLOW_RADTRANS
                0021 #include "RADTRANS_SIZE.h"
                0022 #endif
                0023 #include "DARWIN_SIZE.h"
                0024 #include "DARWIN_INDICES.h"
                0025 #include "DARWIN_DIAGS.h"
                0026 #include "DARWIN_RADTRANS.h"
                0027 #include "DARWIN_PARAMS.h"
                0028 #include "DARWIN_TRAITS.h"
                0029
                0030 C !INPUT PARAMETERS: ===================================================
                0031 C  Ptr    :: darwin model tracers
                0032 C  PAR    :: PAR in uEin/s/m2
                0033 C         :: either non-spectral (tlam=1) or waveband total
                0034 C  myTime :: current time
                0035 C  myIter :: current iteration number
                0036 C  myThid :: thread number
                0037       _RL Ptr(nDarwin)
                0038       _RL PAR(nlam)
                0039       _RL photoTempFunc(nplank)
3f05298721 Oliv*0040       _RL hetTempFunc(nplank)
8fbfd1f382 Oliv*0041       _RL reminTempFunc
                0042       _RL uptakeTempFunc
                0043       _RL grazTempFunc(nplank)
                0044       _RL mortTempFunc
                0045       _RL mort2TempFunc
a092808e6b shlo*0046       _RL macromolTempFunc
ba0b6d5d33 Oliv*0047       _RL omegaCl
8fbfd1f382 Oliv*0048       _RL myTime
6559acf5ff Oliv*0049       INTEGER k_debug
                0050       INTEGER myThid, myIter
8fbfd1f382 Oliv*0051
                0052 C !INPUT/OUTPUT PARAMETERS: ============================================
                0053 C  gTr    :: accumulates computed tendencies
                0054       _RL gTr(nDarwin)
                0055
                0056 C !OUTPUT PARAMETERS: ==================================================
5e7acb36b1 daat*0057 C  chlout  :: computed acclimated chlorophyll if not dynamic
                0058 C  diags   :: pass diagnostics back to darwin_forcing
                0059 C  sumprey :: pass total prey per zoo back for DVM
8fbfd1f382 Oliv*0060       _RL chlout(nPhoto)
                0061       _RL diags(darwin_nDiag)
5e7acb36b1 daat*0062       _RL sumprey(nplank)
8fbfd1f382 Oliv*0063 CEOP
                0064
                0065 #ifdef ALLOW_DARWIN
                0066
                0067 c !LOCAL VARIABLES: ====================================================
                0068       INTEGER j, l
                0069       INTEGER jz, jp
                0070
                0071       _RL DIC
                0072       _RL NH4
                0073       _RL NO2
                0074       _RL NO3
                0075       _RL PO4
                0076       _RL SiO2
                0077       _RL FeT
                0078       _RL DOC
                0079       _RL DON
                0080       _RL DOP
                0081       _RL DOFe
                0082       _RL POC
                0083       _RL PON
                0084       _RL POP
                0085       _RL POSi
                0086       _RL POFe
                0087       _RL pIC
                0088       _RL O2
                0089
                0090       _RL X(nplank)
                0091       _RL Qc(nplank)
c7b6c66d45 Oliv*0092       _RL Qctot(nplank)
8fbfd1f382 Oliv*0093       _RL Qn(nplank)
                0094       _RL Qp(nplank)
                0095       _RL Qsi(nplank)
                0096       _RL Qfe(nplank)
                0097 #ifdef DARWIN_ALLOW_CHLQUOTA
                0098       _RL QChl(nPhoto)
                0099 #endif
c7b6c66d45 Oliv*0100       _RL Qch(nplank)
                0101
6b11238516 Oliv*0102       _RL regQc
                0103       _RL regQn
                0104       _RL regQp
                0105       _RL regQfe
                0106       _RL regQsi
                0107
                0108       _RL MM_PO4
                0109       _RL MM_SiO2
                0110       _RL MM_FeT
                0111       _RL MM_NH4
                0112       _RL MM_NO2
                0113       _RL MM_NO3
8fbfd1f382 Oliv*0114
                0115       _RL limitpCO2
                0116       _RL limitNH4
                0117       _RL limitNO2
                0118       _RL limitNO3
                0119       _RL fracNH4
                0120       _RL fracNO2
                0121       _RL fracNO3
                0122       _RL limitn
                0123       _RL limitp
                0124       _RL limitsi
                0125       _RL limitfe
                0126       _RL limitnut
                0127       _RL limitI
                0128
                0129       _RL muPON
                0130       _RL muPOC
                0131       _RL muPOP
                0132       _RL muPOFe
                0133       _RL muDON
                0134       _RL muDOC
                0135       _RL muDOP
                0136       _RL muDOFe
                0137       _RL muO
                0138       _RL mu
                0139
                0140       _RL uptakeDIC
                0141       _RL uptakeNH4
                0142       _RL uptakeNO2
                0143       _RL uptakeNO3
                0144       _RL uptakeN
                0145       _RL uptakePO4
                0146       _RL uptakeSiO2
                0147       _RL uptakeFeT
                0148       _RL consumDIC
                0149       _RL consumDIC_PIC
                0150       _RL consumNH4
                0151       _RL consumNO2
                0152       _RL consumNO3
                0153       _RL consumPO4
                0154       _RL consumSiO2
                0155       _RL consumFeT
b6fafcd98a Oliv*0156       _RL fIph_sum
                0157       _RL fTph_sum
                0158       _RL fnut_sum
                0159       _RL X_photo
8fbfd1f382 Oliv*0160
                0161       _RL uptakePON
                0162       _RL uptakePOP
                0163       _RL uptakePOC
                0164       _RL uptakePOFe
                0165       _RL uptakeDON
                0166       _RL uptakeDOP
                0167       _RL uptakeDOC
                0168       _RL uptakeDOFe
                0169       _RL uptakeO2
                0170
                0171       _RL respPON
                0172       _RL respPOP
                0173       _RL respPOC
                0174       _RL respPOFe
                0175       _RL respPOSi
                0176       _RL respDON
                0177       _RL respDOP
                0178       _RL respDOC
                0179       _RL respDOFe
                0180
                0181       _RL hydrolPON
                0182       _RL hydrolPOP
                0183       _RL hydrolPOC
                0184       _RL hydrolPOFe
                0185       _RL solubilPON
                0186       _RL solubilPOP
                0187       _RL solubilPOC
                0188       _RL solubilPOFe
                0189
                0190       _RL consumPON
                0191       _RL consumPOP
                0192       _RL consumPOC
                0193       _RL consumPOFe
                0194       _RL consumPOSi
                0195       _RL consumDON
                0196       _RL consumDOP
                0197       _RL consumDOC
                0198       _RL consumDOFe
                0199       _RL consumO2
                0200
                0201       _RL inhibNH4
                0202
c7b6c66d45 Oliv*0203       _RL photoSyn
8fbfd1f382 Oliv*0204       _RL alpha_I
                0205       _RL alpha_I_growth
                0206       _RL PCm
                0207       _RL PC
                0208       _RL acclim
                0209       _RL chl2c
                0210       _RL growth
                0211       _RL rhochl
                0212       _RL Ek
                0213       _RL EkoverE
                0214
                0215       _RL synthChl
                0216
                0217       _RL reminDOC
                0218       _RL reminDON
                0219       _RL reminDOP
                0220       _RL reminDOFe
                0221       _RL reminPOC
                0222       _RL reminPON
                0223       _RL reminPOP
                0224       _RL reminPOSi
                0225       _RL reminPOFe
                0226       _RL disscPIC
                0227
                0228       _RL prodNO2
                0229       _RL prodNO3
                0230
                0231       _RL PARtot
                0232
                0233       _RL tmp
                0234
a092808e6b shlo*0235 #ifdef DARWIN_MACROMOLECULAR_GROWTH
                0236       _RL mmX, mmQn, mmQp, mmQfe
                0237       _RL PChl
                0238       _RL QC_chl_D0
                0239       _RL QC_chl_D1
                0240       _RL QN_chl_D0
                0241       _RL QN_chl_D1
                0242       _RL QN_photo_D0
                0243       _RL QN_photo_D1
                0244       _RL QN_bio_D1
                0245       _RL QN_pro_D0
                0246       _RL QN_pro_D1
                0247       _RL QN_rna_D1
                0248       _RL QN_rna_D2
                0249       _RL Qn_D0
                0250       _RL Qn_D1
                0251       _RL Qn_D2
                0252       _RL Qn_exc
                0253       _RL QC_chl
                0254       _RL QN_chl
                0255       _RL QN_pro_photo
                0256       _RL QN_pro_bio
                0257       _RL QN_pro_tot
                0258       _RL QN_RNA
                0259       _RL QnExc
                0260       _RL exQn
                0261       _RL QC_thy_D0
                0262       _RL QC_thy_D1
                0263       _RL QC_pro_D0
                0264       _RL QC_pro_D1
                0265       _RL QC_rna_D1
                0266       _RL QC_rna_D2
                0267       _RL Qc_D0
                0268       _RL Qc_D1
                0269       _RL Qc_D2
                0270       _RL Qc_exc
                0271 C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
                0272       _RL QP_thy_D0
                0273       _RL QP_thy_D1
                0274       _RL QP_rna_D1
                0275       _RL QP_rna_D2
                0276       _RL Qp_D0
                0277       _RL Qp_D1
                0278       _RL Qp_D2
                0279       _RL PC_n
                0280       _RL Qp_exc
                0281       _RL PC_p
                0282       _RL QP_thy
                0283       _RL QP_RNA
                0284       _RL QN_essential
                0285       _RL QP_essential
                0286       _RL QN_max
                0287       _RL QC_essential
                0288       _RL PC_c
                0289       _RL QP_store
                0290       _RL QP_excess
                0291       _RL QN_store
                0292       _RL QN_excess
                0293       _RL QpExc
                0294       _RL exQp
                0295       _RL VN,VP,VFe,VSi
                0296 C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
                0297       _RL exQc
                0298 C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
                0299       _RL Qfe_D1
                0300       _RL Qfe_D0
                0301       _RL Qfe_exc
                0302       _RL PC_fe
                0303       _RL Qfe_pro_photo
                0304       _RL Qfe_essential
                0305       _RL Qfe_Store
                0306       _RL Qfe_excess
                0307       _RL exQfe
                0308       _RL QfeExc
                0309       _RL MODE
                0310       _RL limN
                0311       _RL limP
                0312       _RL limF
                0313       _RL limC
                0314       _RL limL
                0315 C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
                0316       _RL QC_chlN
                0317       _RL QC_chlP
                0318       _RL QC_chlF
                0319       _RL Scal_coeff
                0320       _RL QC_store
                0321       _RL RQ
                0322       _RL QN_DNA_actl
                0323       _RL QP_other_actl
                0324       _RL QP_DNA_actl
                0325 C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
                0326 #endif /* DARWIN_MACROMOLECULAR_GROWTH */
                0327
8fbfd1f382 Oliv*0328 #ifdef DARWIN_ALLOW_CDOM
                0329       _RL CDOM
                0330       _RL reminPOC_CDOM
                0331       _RL reminPON_CDOM
                0332       _RL reminPOP_CDOM
                0333       _RL reminPOFe_CDOM
                0334       _RL degrCDOM_DOC
                0335       _RL degrCDOM_DON
                0336       _RL degrCDOM_DOP
                0337       _RL degrCDOM_DOFe
                0338 #endif
                0339
                0340 #ifdef DARWIN_ALLOW_DENIT
                0341       _RL denit, denitNH4
                0342 #endif
                0343
c7b6c66d45 Oliv*0344 #ifdef DARWIN_ALLOW_CSTORE
                0345       _RL Dmd_NC
                0346       _RL Dmd_PC
                0347       _RL Dmd_SiC
                0348       _RL Dmd_FeC
                0349       _RL Dmd_min
                0350       _RL excretC
                0351       _RL excC
                0352       _RL storeC
                0353 #endif
                0354
8fbfd1f382 Oliv*0355 C for grazing
                0356
5e7acb36b1 daat*0357       _RL sumpref, grazphy
8fbfd1f382 Oliv*0358
                0359       _RL preygraz   (nplank)
                0360       _RL preygrazexp(nplank)
a93256cee4 Oliv*0361       _RL grazrate   (nplank)
8fbfd1f382 Oliv*0362       _RL predgrazc  (nplank)
                0363 #ifdef DARWIN_ALLOW_NQUOTA
                0364       _RL predgrazn  (nplank)
                0365 #endif
                0366 #ifdef DARWIN_ALLOW_PQUOTA
                0367       _RL predgrazp  (nplank)
                0368 #endif
                0369 #ifdef DARWIN_ALLOW_FEQUOTA
                0370       _RL predgrazfe (nplank)
                0371 #endif
                0372
                0373       _RL predexpc, predexpn, predexpp, predexpfe
                0374       _RL graz2OC, graz2ON, graz2OP, graz2OFe
                0375       _RL graz2POC, graz2PON, graz2POP, graz2POSi, graz2POFe
                0376       _RL graz2PIC
                0377
                0378       _RL expfrac
                0379
                0380       _RL Xe
                0381       _RL mortX
                0382       _RL mortX2
                0383
                0384       _RL exude_DOC
                0385       _RL exude_DON
                0386       _RL exude_DOP
                0387       _RL exude_DOFe
                0388
                0389       _RL exude_PIC
                0390       _RL exude_POC
                0391       _RL exude_PON
                0392       _RL exude_POP
                0393       _RL exude_POSi
                0394       _RL exude_POFe
                0395
                0396       _RL mort_c(nplank)
                0397
                0398       _RL respir
                0399       _RL respir_c
8752612d6f Oliv*0400       _RL respirN
                0401       _RL respirP
                0402       _RL respirFe
                0403       _RL respirSi
8fbfd1f382 Oliv*0404
                0405 #ifdef DARWIN_ALLOW_CDOM
                0406       _RL graz2CDOM, exude_CDOM
                0407 #endif
                0408
5e7acb36b1 daat*0409       _RL mortDVM_c
                0410       _RL mortDVM_DOC
                0411       _RL mortDVM_DON
                0412       _RL mortDVM_DOP
                0413       _RL mortDVM_DOFe
                0414       _RL mortDVM_POC
                0415       _RL mortDVM_PON
                0416       _RL mortDVM_POP
                0417       _RL mortDVM_POFe
                0418       _RL mortDVM_POSi
                0419
b854db6bc4 Oliv*0420       _RL DARWIN_EPS
                0421       PARAMETER(DARWIN_EPS=1 _d -38)
                0422
c7b6c66d45 Oliv*0423 C=======================================================================
                0424 C==== Precompute a few things ==========================================
8fbfd1f382 Oliv*0425
                0426       PARtot = SUM(PAR)
                0427
c7b6c66d45 Oliv*0428 C==== Make all bio fields non-negative and compute quotas ==============
8fbfd1f382 Oliv*0429
086a45f245 Oliv*0430       DIC  = MAX(0 _d 0, Ptr(iDIC))
                0431       NH4  = MAX(0 _d 0, Ptr(iNH4))
                0432       NO2  = MAX(0 _d 0, Ptr(iNO2))
                0433       NO3  = MAX(0 _d 0, Ptr(iNO3))
                0434       PO4  = MAX(0 _d 0, Ptr(iPO4))
                0435       SiO2 = MAX(0 _d 0, Ptr(iSiO2))
                0436       FeT  = MAX(0 _d 0, Ptr(iFeT))
                0437
                0438       DOC  = MAX(0 _d 0, Ptr(iDOC))
                0439       DON  = MAX(0 _d 0, Ptr(iDON))
                0440       DOP  = MAX(0 _d 0, Ptr(iDOP))
                0441       DOFe = MAX(0 _d 0, Ptr(iDOFe))
                0442
                0443       PIC  = MAX(0 _d 0, Ptr(iPIC))
                0444       POC  = MAX(0 _d 0, Ptr(iPOC))
                0445       PON  = MAX(0 _d 0, Ptr(iPON))
                0446       POP  = MAX(0 _d 0, Ptr(iPOP))
                0447       POSi = MAX(0 _d 0, Ptr(iPOSi))
                0448       POFe = MAX(0 _d 0, Ptr(iPOFe))
8fbfd1f382 Oliv*0449 #ifdef DARWIN_ALLOW_CARBON
086a45f245 Oliv*0450       O2   = MAX(0 _d 0, Ptr(iO2))
8fbfd1f382 Oliv*0451 #endif
                0452 #ifdef DARWIN_ALLOW_CDOM
086a45f245 Oliv*0453       CDOM = MAX(0 _d 0, Ptr(iCDOM))
8fbfd1f382 Oliv*0454 #endif
                0455
                0456       DO j = 1, nplank
                0457 C fixed carbon quota, for now 1.0 (may change later)
                0458         Qc(j) = 1.0
c7b6c66d45 Oliv*0459         Qctot(j) = 1.0
086a45f245 Oliv*0460         X(j) = MAX(0 _d 0, Ptr(ic+j-1))/Qc(j)
c7b6c66d45 Oliv*0461 C other elements get quota from corresponding ptracer or set to fixed
8fbfd1f382 Oliv*0462 c ratio if not variable.
                0463 #ifdef DARWIN_ALLOW_NQUOTA
b854db6bc4 Oliv*0464         IF (X(j) .GT. DARWIN_EPS) THEN
c697368a52 Oliv*0465           Qn(j) = MAX(0 _d 0, Ptr(in+j-1))/X(j)
                0466         ELSE
                0467           Qn(j) = Qnmax(j)
                0468         ENDIF
8fbfd1f382 Oliv*0469 #else
                0470         Qn(j) = R_NC(j)
                0471 #endif
                0472 #ifdef DARWIN_ALLOW_PQUOTA
b854db6bc4 Oliv*0473         IF (X(j) .GT. DARWIN_EPS) THEN
c697368a52 Oliv*0474           Qp(j) = MAX(0 _d 0, Ptr(ip+j-1))/X(j)
                0475         ELSE
                0476           Qp(j) = Qpmax(j)
                0477         ENDIF
8fbfd1f382 Oliv*0478 #else
                0479         Qp(j) = R_PC(j)
                0480 #endif
                0481 #ifdef DARWIN_ALLOW_SIQUOTA
b854db6bc4 Oliv*0482         IF (X(j) .GT. DARWIN_EPS) THEN
c697368a52 Oliv*0483           Qsi(j) = MAX(0 _d 0, Ptr(isi+j-1))/X(j)
                0484         ELSE
                0485           Qsi(j) = Qsimax(j)
                0486         ENDIF
8fbfd1f382 Oliv*0487 #else
                0488         Qsi(j) = R_SiC(j)
                0489 #endif
                0490 #ifdef DARWIN_ALLOW_FEQUOTA
b854db6bc4 Oliv*0491         IF (X(j) .GT. DARWIN_EPS) THEN
c697368a52 Oliv*0492           Qfe(j) = MAX(0 _d 0, Ptr(ife+j-1))/X(j)
                0493         ELSE
                0494           Qfe(j) = Qfemax(j)
                0495         ENDIF
8fbfd1f382 Oliv*0496 #else
                0497         Qfe(j) = R_FeC(j)
                0498 #endif
                0499       ENDDO
                0500 #ifdef DARWIN_ALLOW_CHLQUOTA
                0501       DO j = 1, nPhoto
b854db6bc4 Oliv*0502         IF (X(j) .GT. DARWIN_EPS) THEN
c697368a52 Oliv*0503           QChl(j) = MAX(0 _d 0, Ptr(iChl+j-1))/X(j)
                0504         ELSE
                0505           QChl(j) = chl2cmax(j)
                0506         ENDIF
8fbfd1f382 Oliv*0507       ENDDO
                0508 #endif
c7b6c66d45 Oliv*0509 #ifdef DARWIN_ALLOW_CSTORE
                0510       DO j = 1, nPhoto
b854db6bc4 Oliv*0511         IF (X(j) .GT. DARWIN_EPS) THEN
c7b6c66d45 Oliv*0512           Qch(j) = MAX(0 _d 0, Ptr(ich+j-1))/X(j)
                0513         ELSE
                0514           Qch(j) = 0 _d 0
                0515         ENDIF
                0516         Qctot(j) = Qctot(j) + Qch(j)
                0517       ENDDO
                0518 #endif
8fbfd1f382 Oliv*0519
c7b6c66d45 Oliv*0520 C==== Initialize accumulators ==========================================
8fbfd1f382 Oliv*0521       consumDIC  = 0.0
                0522       consumDIC_PIC = 0.0
                0523       consumNH4  = 0.0
                0524       consumNO2  = 0.0
                0525       consumNO3  = 0.0
                0526       consumPO4  = 0.0
                0527       consumSiO2 = 0.0
                0528       consumFeT  = 0.0
                0529       consumPON  = 0.0
                0530       consumPOP  = 0.0
                0531       consumPOC = 0.0
                0532       consumPOFe  = 0.0
                0533       consumPOSi  = 0.0
                0534       consumDON  = 0.0
                0535       consumDOP  = 0.0
                0536       consumDOC = 0.0
                0537       consumDOFe  = 0.0
                0538       consumO2 = 0.0
                0539       reminPON  = 0.0
                0540       reminPOP  = 0.0
                0541       reminPOC = 0.0
                0542       reminPOFe  = 0.0
                0543       reminPOSi  = 0.0
                0544       reminDON  = 0.0
                0545       reminDOP  = 0.0
                0546       reminDOC = 0.0
                0547       reminDOFe  = 0.0
                0548       solubilPON  = 0.0
                0549       solubilPOP  = 0.0
                0550       solubilPOC = 0.0
                0551       solubilPOFe  = 0.0
                0552       prodNO2 = 0.0
                0553       prodNO3 = 0.0
8752612d6f Oliv*0554       respir = 0.0 _d 0
b6fafcd98a Oliv*0555       fIph_sum = 0.0 _d 0
                0556       fTph_sum = 0.0 _d 0
                0557       fnut_sum = 0.0 _d 0
                0558       X_photo = 0.0 _d 0
5e7acb36b1 daat*0559
c7b6c66d45 Oliv*0560 #ifdef DARWIN_ALLOW_CSTORE
                0561       excC = 0.0
                0562 #endif
8fbfd1f382 Oliv*0563
                0564       DO j = 1, nPhoto
                0565         chlout(j) = 0.0 _d 0
                0566       ENDDO
                0567
                0568       DO l=1,darwin_nDiag
                0569         diags(l) = 0.0
                0570       ENDDO
                0571
                0572 C=======================================================================
c7b6c66d45 Oliv*0573 C==== Phototrophy ======================================================
8fbfd1f382 Oliv*0574
                0575       DO j = 1, nPhoto
                0576        IF (isPhoto(j) .NE. 0) THEN
                0577
a092808e6b shlo*0578 #ifdef DARWIN_MACROMOLECULAR_GROWTH
                0579
                0580         IF (X(j) .GT. DARWIN_EPS) THEN
                0581          mmQn  = Qn(j)
                0582          mmQp  = Qp(j)
                0583          mmQfe = Qfe(j)
                0584          mmX   = X(j)
                0585         ELSE
                0586          mmQn = 0 _d 0
                0587          mmQp = 0 _d 0
                0588          mmQfe = 0 _d 0
                0589          mmX = 0 _d 0
                0590         ENDIF
                0591
                0592         limP = 0 _d 0
                0593         limN = 0 _d 0
                0594         limF = 0 _d 0
                0595         limC = 0 _d 0
                0596         limL = 0 _d 0
                0597
                0598 C       set up Chl and N requirements as a function of growth rate PC
                0599 C       then solve for PC
                0600
                0601 C       Carbohydrate fixation rate per chlorophyll
                0602 C       (molC (mol C in chl s)-1)
                0603         PChl = Sf(j)*VI_max(j)*(1 - EXP(-A_I(j)*PARtot))
                0604
                0605 C       Case 1: Enought light, growth rate is a function of macromol
                0606 C       allocation
                0607         IF ( PChl .GT. VI_min(j) ) THEN
                0608
                0609 C        chlorophyll quota (mol C cell-1)
                0610 C        Chl = ((1 _d 0 + E)*ls + m)/PChl
                0611 C        where ls = D*Qc               (molC s-1) Biomass synthesis rate
                0612          QC_chl_D0 = maintConsum(j)/PChl
                0613          QC_chl_D1 = (1 + ECo2Prod(j))/PChl
                0614
                0615 C        nitrogen pools in mol N (mol C)-1
                0616          QN_chl_D0 = QC_chl_D0*Y_NC_chl(j)
                0617          QN_chl_D1 = QC_chl_D1*Y_NC_chl(j)
                0618
                0619 C        photosynthesis related protein
                0620          QN_photo_D0 = QC_chl_D0*A_pho(j)/Y_CN_protein(j)
                0621          QN_photo_D1 = QC_chl_D1*A_pho(j)/Y_CN_protein(j)
                0622
                0623 C        biosynthesis related protein
                0624          QN_bio_D1 = A_bio(j)/Y_CN_protein(j)/macromolTempFunc
                0625
                0626 C        total protein
                0627          QN_pro_D0 = QN_photo_D0 + QN_pro_other(j)
                0628          QN_pro_D1 = QN_photo_D1 + QN_bio_D1
                0629
                0630 C        nitrogen in RNA
                0631          QN_rna_D1 = QN_pro_D0*AN_RNA(j)/macromolTempFunc
                0632          QN_rna_D2 = QN_pro_D1*AN_RNA(j)/macromolTempFunc
                0633
                0634 C        Total nitrogen
                0635          Qn_D0 = QN_chl_D0 + QN_pro_D0 + QN_RNA_min(j) + QN_DNA(j)
                0636          Qn_D1 = QN_chl_D1 + QN_pro_D1 + QN_rna_D1
                0637          Qn_D2 = QN_rna_D2
                0638
                0639
                0640 C        Phosphate restriction on growth rate
                0641          QP_rna_D1 = QN_rna_D1*Y_PN_nucacid(j)
                0642          QP_rna_D2 = QN_rna_D2*Y_PN_nucacid(j)
                0643          QP_thy_D0 = QC_chl_D0*Y_THY_P(j)
                0644          QP_thy_D1 = QC_chl_D1*Y_THY_P(j)
                0645
                0646          Qp_D0 = QP_thy_D0 + QP_RNA_min(j) + QP_DNA(j) + QP_other(j)
                0647          Qp_D1 = QP_thy_D1 + QP_rna_D1
                0648          Qp_D2 = QP_rna_D2
                0649
                0650
                0651 C        Iron restriction on growth
                0652          Qfe_D0 = QN_photo_D0*Y_FeN_photo(j)
                0653          Qfe_D1 = QN_photo_D1*Y_FeN_photo(j)
                0654
                0655
                0656 C        carbon restriction on growth rate
                0657          QC_thy_D0 = QP_thy_D0*Y_CP_Plip(j)
                0658          QC_thy_D1 = QP_thy_D1*Y_CP_Plip(j)
                0659          QC_pro_D0 = QN_pro_D0*Y_CN_protein(j)
                0660          QC_pro_D1 = QN_pro_D1*Y_CN_protein(j)
                0661          QC_rna_D1 = QN_rna_D1*Y_CN_RNA(j)
                0662          QC_rna_D2 = QN_rna_D2*Y_CN_RNA(j)
                0663
                0664          Qc_D0 = QC_chl_D0 + QC_pro_D0 + QC_RNA_min(j) + QC_thy_D0
                0665      &         + QC_DNA(j) + QC_other(j)
                0666          Qc_D1 = QC_chl_D1 + QC_pro_D1 + QC_thy_D1 + QC_rna_D1
                0667          Qc_D2 = QC_rna_D2
                0668
                0669
                0670 C        Determine what is limiting
                0671          Qn_exc = MAX(0 _d 0, mmQn - Qn_D0)
                0672          PC_n = 2*Qn_exc/(Qn_D1 + SQRT(Qn_D1*Qn_D1 + 4*Qn_D2*Qn_exc))
                0673          Qp_exc = MAX(0 _d 0, mmQp - Qp_D0)
                0674          PC_p = 2*Qp_exc/(Qp_D1 + SQRT(Qp_D1*Qp_D1 + 4*Qp_D2*Qp_exc))
                0675          Qfe_exc = MAX(0 _d 0, mmQfe - Qfe_D0)
                0676          PC_fe = Qfe_exc/Qfe_D1
                0677          Qc_exc = MAX(0 _d 0, Qc(j) - Qc_D0)
                0678          PC_c = 2*Qc_exc/(Qc_D1 + SQRT(Qc_D1*Qc_D1 + 4*Qc_D2*Qc_exc))
                0679
                0680          PC = MIN(PC_n,PC_p,PC_fe,PC_c)
                0681          IF (PC_p.EQ.PC) THEN
                0682           MODE = 1 _d 0
                0683           limP = 1 _d 0
                0684          ELSEIF (PC_n.EQ.PC) THEN
                0685           MODE = 2 _d 0
                0686           limN = 1 _d 0
                0687          ELSEIF (PC_fe.EQ.PC) THEN
                0688           MODE = 3 _d 0
                0689           limF = 1 _d 0
                0690          ELSEIF (PC_c.EQ.PC) THEN
                0691           MODE = 4 _d 0
                0692           limC = 1 _d 0
                0693          ELSE
                0694           MODE = 0 _d 0
                0695          ENDIF
                0696
                0697          QC_chl = QC_chl_D0 + QC_chl_D1*PC
                0698
                0699 C       Case 2: Light limitation, growth rate is zero, QC_chl is
                0700 C       calculated based on cellular mass balance
                0701         ELSE
                0702          MODE = 4 _d 0
                0703          limL = 1 _d 0
                0704          PC = 0. _d 0
                0705         ENDIF
                0706
                0707 C       If growth rate is zero, recompute maximum QC_chl possible with
                0708 C       given quotas
                0709         IF (PC .EQ. 0. _d 0) THEN
                0710          QC_chlN = (mmQn - QnNoChl(j))
                0711      &         /(Y_NC_chl(j) + A_pho(j)/Y_CN_protein(j))
                0712          QC_chlP = (mmQp - QpNoChl(j))/Y_THY_P(j)
                0713          QC_chlF = (mmQfe - QfeNoChl(j))
                0714      &         /(A_pho(j)/Y_CN_protein(j)*Y_FeN_photo(j))
                0715          QC_chl = MAX(0 _d 0, MIN(QC_chlN,QC_chlP,QC_chlF,QC_chlMax(j)))
                0716         ENDIF
                0717
                0718 C       convert to mg Chl m-3
                0719         chlout(j) = 893.49 _d 0*X(j)*QC_chl/55 _d 0
                0720
                0721 C       Calculate diagnostics and essentials
                0722         QN_chl       = QC_chl*Y_NC_chl(j)
                0723         QN_pro_photo = QC_chl*A_pho(j)/Y_CN_protein(j)
                0724         QN_pro_bio   = A_bio(j)/Y_CN_protein(j)*PC/macromolTempFunc
                0725         QN_pro_tot   = QN_pro_other(j) + QN_pro_photo + QN_pro_bio
                0726         QN_RNA       = QN_RNA_min(j)
                0727      &               + QN_pro_tot*AN_RNA(j)*PC/macromolTempFunc
                0728
                0729         QP_thy = QC_chl*Y_THY_P(j)
                0730         QP_RNA = QN_RNA*Y_PN_nucacid(j)
                0731
                0732         Qfe_pro_photo = QN_pro_photo*Y_FeN_photo(j)
                0733
                0734 C       Case 3: if Qn, Qp are lower then the minimum requirement for Chl
                0735 C       the remaining pools are scaled down to not exceed quotas
                0736         IF (QC_chl.EQ.0 _d 0) THEN
                0737          Scal_coeff  = MIN(mmQp/QpNoChl(j),mmQn/QnNoChl(j),1 _d 0)
                0738          QN_RNA      = QN_RNA_min(j)*Scal_coeff
                0739          QN_DNA_actl = QN_DNA(j)*Scal_coeff
                0740          QP_RNA      = QP_RNA_min(j)*Scal_coeff
                0741          QP_DNA_actl = QP_DNA(j)*Scal_coeff
                0742 C        excess N (if Qp.LT.QpNoChl) goes here
                0743          QN_pro_tot  = MIN(mmQn - QN_RNA - QN_DNA_actl, QN_pro_other(j))
                0744 C        excess P (if Qn.LT.QnNoChl) goes here
                0745          QP_other_actl = MIN(mmQp - QP_RNA - QP_DNA_actl, QP_other(j))
                0746         ELSE
                0747          QN_DNA_actl   = QN_DNA(j)
                0748          QP_DNA_actl   = QP_DNA(j)
                0749          QP_other_actl = QP_other(j)
                0750         ENDIF
                0751
                0752 C       total pools required for functioning (or scaled down)
                0753         QN_essential  = QN_chl + QN_pro_tot + QN_RNA + QN_DNA_actl
                0754         QN_max        = QN_essential + QN_sto_max(j)
                0755         QP_essential  = QP_thy + QP_RNA + QP_DNA_actl +  QP_other_actl
                0756         Qfe_essential = Qfe_pro_photo
                0757         QC_essential  = QC_chl + QN_pro_tot*Y_CN_protein(j)
                0758      &                + QP_thy*Y_CP_Plip(j) + QN_RNA*Y_CN_RNA(j)
                0759      &                + QN_DNA_actl*Y_CN_DNA(j) + QC_other(j)
                0760
                0761 C Storage calculation -----------------
                0762
                0763 C       N quota above QN_essential is store, but not more than
                0764 C       QN_sto_max(j) and make sure enough C for store
                0765         QN_store = MIN(mmQn - QN_essential,
                0766      &                 (Qc(j) - QC_essential)/Y_CN_cyano(j),
                0767      &                 QN_sto_max(j))
                0768         QN_store = MAX(0 _d 0, QN_store)
                0769         QN_excess = MAX(0 _d 0, mmQn - QN_essential - QN_store)
                0770
                0771 C       P quota above QP_essential is store, but keep total P under
                0772 C       Qp_max(j)
                0773         QP_store = MAX(0 _d 0, MIN(mmQp, Qp_max(j)) - QP_essential)
                0774         QP_excess = MAX(0 _d 0, mmQp - QP_essential - QP_store)
                0775
                0776 C       Fe quota above Qfe_essential is store, but keep total iron under
                0777 C       Qfe_max(j)
                0778         Qfe_store = MAX(0 _d 0, MIN(mmQfe, Qfe_max(j)) - Qfe_essential)
                0779         Qfe_excess = MAX(0 _d 0, mmQfe - Qfe_essential - Qfe_store)
                0780
                0781 C       deduce carbon store
                0782         QC_store = MAX(0 _d 0,
                0783      &                 1 _d 0 - QC_essential - QN_store*Y_CN_cyano(j))
                0784
                0785 C shlomit: new nutrient-uptake regulation terms:
                0786 C when store is full, phytoplankton reduce their nutrient-uptake:
                0787
                0788 C       Nitrogen:
                0789         IF ((QN_essential + QN_store) .GT. 0 _d 0) THEN
                0790          regQn = (1.1 _d 0*(QN_essential + QN_store) - mmQn)/
                0791      &             (0.1 _d 0*(QN_essential + QN_store))
                0792          regQn = MAX(0 _d 0, MIN(1 _d 0, regQn))
                0793         ELSE
                0794          regQn = 1 _d 0
                0795         ENDIF
                0796         MM_NO3 = NO3/(NO3 + ksatNO3(j))
                0797 C       Carbon-specific NO3 uptake rate
                0798         VN = vmaxNO3(j)*MM_NO3*regQn*uptakeTempFunc
                0799         uptakeN = VN*mmX
                0800         uptakeNO3 = VN*mmX
                0801         uptakeNO2 = 0 _d 0
                0802         uptakeNH4 = 0 _d 0
                0803
                0804 C       Phosphorous:
                0805         IF ((QP_essential+QP_store) .GT. 0 _d 0) THEN
                0806          regQp = (1.1 _d 0*(QP_essential+QP_store) - mmQp)/
                0807      &             (0.1 _d 0*(QP_essential+QP_store))
                0808          regQp = MAX(0 _d 0, MIN(1 _d 0, regQp))
                0809         ELSE
                0810          regQp = 1 _d 0
                0811         ENDIF
                0812         MM_PO4 = PO4/(PO4 + ksatPO4(j))
                0813 C       Carbon-specific PO4 uptake rate
                0814         VP = vmaxPO4(j)*MM_PO4*regQp*uptakeTempFunc
                0815         uptakePO4 = VP*mmX
                0816
                0817 C       Iron:
                0818         IF ((Qfe_essential + Qfe_Store) .GT. 0 _d 0) THEN
                0819          regQFe = (1.1 _d 0*(Qfe_essential + Qfe_Store) - mmQfe)/
                0820      &              (0.1 _d 0*(Qfe_essential + Qfe_Store))
                0821          regQFe = MAX(0 _d 0, MIN(1 _d 0, regQFe))
                0822         ELSE
                0823          regQFe = 1 _d 0
                0824         ENDIF
                0825         MM_FeT = FeT/(FeT + ksatFeT(j))
                0826 C       Carbon-specific FeT uptake rate
                0827         Vfe = vmaxFeT(j)*MM_FeT*regQFe*uptakeTempFunc
                0828         uptakeFeT = Vfe*mmX
                0829
                0830         exQc = 0 _d 0
                0831
                0832 C       shlomit: Calclate the O2:C ratio
                0833         RQ = ( QN_pro_tot*Y_CN_protein(j))* 1.56 _d 0
                0834      &     + ( QC_other(j)/2  + QC_store/2) * 1.34 _d 0
                0835      &     + ( Y_CP_Plip(j) * QP_thy)* 1.30 _d 0
                0836      &     + ( QC_other(j)/2 + QC_store/2) * 1.0 _d 0
                0837      &     + ( QC_chl) * 1.39 _d 0
                0838      &     + ( QN_DNA_actl*Y_CN_DNA(j)) * 1.62 _d 0
                0839      &     + ( QN_RNA*Y_CN_RNA(j)) *  1.57 _d 0
                0840      &     + ( QN_store*Y_CN_cyano(j)) * 1.85 _d 0
                0841
                0842 # ifdef DARWIN_DIAG_PERTYPE
                0843         diags(iPChl+j-1) = PChl
                0844         diags(iVN  +j-1) = VN
                0845         diags(iVP  +j-1) = VP
                0846         diags(iMODE+j-1) = MODE
                0847         diags(iFe_C+j-1) = mmQfe
                0848         diags(iexQc+j-1) = exQc
                0849         diags(iCChl+j-1) = QC_chl
                0850         diags(iNChl+j-1) = QN_chl
                0851         diags(iNPho+j-1) = QN_pro_photo
                0852         diags(iNSyn+j-1) = QN_pro_bio
                0853         diags(iNPrn+j-1) = QN_pro_tot
                0854         diags(iNRNA+j-1) = QN_RNA
                0855         diags(iNDNA+j-1) = QN_DNA_actl
                0856         diags(iNSTO+j-1) = QN_store
                0857         diags(iNEXC+j-1) = QN_excess
                0858         diags(iPRNA+j-1) = QP_RNA
                0859         diags(iPDNA+j-1) = QP_DNA_actl
                0860         diags(iPTHY+j-1) = QP_thy
                0861         diags(iPCON+j-1) = QP_other_actl
                0862         diags(iPSTO+j-1) = QP_store
                0863         diags(iPEXC+j-1) = QP_excess
                0864         diags(iFPHO+j-1) = Qfe_pro_photo
                0865         diags(iFSTO+j-1) = Qfe_Store
                0866         diags(iFEXC+j-1) = Qfe_excess
                0867         diags(iY_RQ+j-1) = RQ
                0868         diags(ilimN+j-1) = limN
                0869         diags(ilimP+j-1) = limP
                0870         diags(ilimF+j-1) = limF
                0871         diags(ilimC+j-1) = limC
                0872         diags(ilimL+j-1) = limL
                0873 # endif
                0874
                0875 #else /* DARWIN_MACROMOLECULAR_GROWTH */
                0876
c7b6c66d45 Oliv*0877 C==== Uptake and nutrient limitation ===================================
8fbfd1f382 Oliv*0878 C       for quota elements, growth is limiteed by available quota,
                0879 C       for non-quota elements, by available nutrients in medium
                0880
                0881 C       to not use PO4, ..., set ksatPO4=0 and vmaxPO4=0 (if DARWIN_ALLOW_PQUOTA)
                0882 C       or R_PC=0 (if not)
                0883 C       the result will be limitp = 1, uptakePO4 = 0
                0884
c7b6c66d45 Oliv*0885 C       CSTORE/exudation:
                0886 C       Dmd_XC means the amount of carbon needed for biosynthesis when X
                0887 C       is the limiting nutrient.  To find the lowest demand (maximum
                0888 C       achievable growth rate), start with a large value.  Quota
                0889 C       elements will be handled further down in a dedicatd CSTORE
                0890 C       section.
                0891 #ifdef DARWIN_ALLOW_CSTORE
                0892         Dmd_min = HUGE(Dmd_min)
                0893 #endif
                0894
8fbfd1f382 Oliv*0895 c phosphorus
c697368a52 Oliv*0896         IF (ksatPO4(j) .GT. 0 _d 0) THEN
6b11238516 Oliv*0897           MM_PO4 = PO4/(PO4 + ksatPO4(j))
c697368a52 Oliv*0898         ELSE
6b11238516 Oliv*0899           MM_PO4 = 1 _d 0
c697368a52 Oliv*0900         ENDIF
8fbfd1f382 Oliv*0901 #ifdef DARWIN_ALLOW_PQUOTA
086a45f245 Oliv*0902         regQp = MAX(0 _d 0, MIN(1 _d 0, (Qpmax(j)-Qp(j))/
                0903      &                                  (Qpmax(j)-Qpmin(j)) ))
a092808e6b shlo*0904         regQp = regQp**hillnumPO4(j)
6b11238516 Oliv*0905 C       carbon-specific PO4 uptake rate
                0906         uptakePO4 = vmaxPO4(j)*MM_PO4*regQp*uptakeTempFunc*X(j)
8fbfd1f382 Oliv*0907 c       normalized Droop limitation
086a45f245 Oliv*0908         limitp = MIN(1 _d 0, (1.0-Qpmin(j)/MAX(Qpmin(j), Qp(j)))/
                0909      &                       (1.0-Qpmin(j)/Qpmax(j)))
6b11238516 Oliv*0910 #else /* no PQUOTA */
c7b6c66d45 Oliv*0911 # ifdef DARWIN_ALLOW_CSTORE
                0912         IF (R_PC(j) .GT. 0 _d 0) THEN
                0913           Dmd_PC = vmaxPO4(j)*MM_PO4*uptakeTempFunc*X(j)/R_PC(j)
                0914           Dmd_min = MIN(Dmd_min, Dmd_PC)
                0915         ELSE
                0916 C       set to zero to avoid floating point error in diagonstic
                0917           Dmd_PC = 0 _d 0
                0918         ENDIF
                0919 # else
6b11238516 Oliv*0920         limitp = MM_PO4
c7b6c66d45 Oliv*0921 # endif
6b11238516 Oliv*0922 #endif /* pquota */
8fbfd1f382 Oliv*0923
                0924 c silica
c697368a52 Oliv*0925         IF (ksatSiO2(j) .GT. 0 _d 0) THEN
6b11238516 Oliv*0926           MM_SiO2 = SiO2/(SiO2 + ksatSiO2(j))
c697368a52 Oliv*0927         ELSE
6b11238516 Oliv*0928           MM_SiO2 = 1 _d 0
c697368a52 Oliv*0929         ENDIF
8fbfd1f382 Oliv*0930 #ifdef DARWIN_ALLOW_SIQUOTA
086a45f245 Oliv*0931         regQsi = MAX(0 _d 0, MIN(1 _d 0, (Qsimax(j) - Qsi(j))/
                0932      &                                   (Qsimax(j) - Qsimin(j)) ))
a092808e6b shlo*0933         regQsi = regQsi**hillnumSiO2(j)
6b11238516 Oliv*0934         uptakeSiO2 = vmaxSiO2(j)*MM_SiO2*regQsi*uptakeTempFunc*X(j)
8fbfd1f382 Oliv*0935 c       linear limitation
086a45f245 Oliv*0936         limitsi = MAX(0 _d 0, MIN(1 _d 0, (Qsi(j) - Qsimin(j))/
                0937      &                                    (Qsimax(j) - Qsimin(j)) ))
6b11238516 Oliv*0938 #else /* no SIQUOTA */
c7b6c66d45 Oliv*0939 # ifdef DARWIN_ALLOW_CSTORE
                0940         IF (R_SiC(j) .GT. 0 _d 0) THEN
                0941           Dmd_SiC = VmaxSiO2(j)*MM_SiO2*uptakeTempFunc*X(j)/R_SiC(j)
                0942           Dmd_min = MIN(Dmd_min, Dmd_SiC)
                0943         ELSE
                0944           Dmd_SiC = 0 _d 0
                0945         ENDIF
                0946 # else
6b11238516 Oliv*0947         limitsi = MM_SiO2
c7b6c66d45 Oliv*0948 # endif
8fbfd1f382 Oliv*0949 #endif
                0950
                0951 c iron
c697368a52 Oliv*0952         IF (ksatFeT(j) .GT. 0 _d 0) THEN
6b11238516 Oliv*0953           MM_FeT = FeT/(FeT + ksatFeT(j))
c697368a52 Oliv*0954         ELSE
6b11238516 Oliv*0955           MM_FeT = 1 _d 0
c697368a52 Oliv*0956         ENDIF
8fbfd1f382 Oliv*0957 #ifdef DARWIN_ALLOW_FEQUOTA
086a45f245 Oliv*0958         regQfe = MAX(0 _d 0, MIN(1 _d 0, (Qfemax(j)-Qfe(j))/
                0959      &                                   (Qfemax(j)-Qfemin(j)) ))
a092808e6b shlo*0960         regQfe = regQfe**hillnumFeT(j)
6b11238516 Oliv*0961         uptakeFeT = vmaxFeT(j)*MM_FeT*regQfe*uptakeTempFunc*X(j)
8fbfd1f382 Oliv*0962 c       normalized Droop limitation
086a45f245 Oliv*0963         limitfe = MIN(1 _d 0, (1.0-Qfemin(j)/MAX(Qfemin(j), Qfe(j)))/
                0964      &                        (1.0-Qfemin(j)/Qfemax(j)))
6b11238516 Oliv*0965 #else /* no FEQUOTA */
c7b6c66d45 Oliv*0966 C iron limit is applied in Geider
                0967 # ifdef DARWIN_ALLOW_CSTORE
                0968         IF (R_FeC(j) .GT. 0 _d 0) THEN
                0969           Dmd_FeC = vmaxFeT(j)*MM_FeT*uptakeTempFunc*X(j)/R_FeC(j)
                0970           Dmd_min = MIN(Dmd_min, Dmd_FeC)
                0971         ELSE
                0972           Dmd_FeC = 0 _d 0
                0973         ENDIF
                0974 # else
6b11238516 Oliv*0975         limitfe = MM_FeT
c7b6c66d45 Oliv*0976 # endif
8fbfd1f382 Oliv*0977 #endif
                0978
                0979 c nitrogen
                0980 #ifdef DARWIN_ALLOW_NQUOTA
                0981 c       have nitrogen quota
                0982         inhibNH4 = EXP(-amminhib(j)*NH4)
6b11238516 Oliv*0983         MM_NH4 = NH4/(NH4 + ksatNH4(j))
                0984         MM_NO2 = NO2/(NO2 + ksatNO2(j))*inhibNH4
                0985         MM_NO3 = NO3/(NO3 + ksatNO3(j))*inhibNH4
086a45f245 Oliv*0986         regQn = MAX(0 _d 0, MIN(1 _d 0, (Qnmax(j)-Qn(j))/
                0987      &                                  (Qnmax(j)-Qnmin(j)) ))
a092808e6b shlo*0988         regQn = regQn**hillnumDIN(j)
6b11238516 Oliv*0989         uptakeNH4 = vmaxNH4(j)*MM_NH4*regQn*uptakeTempFunc*X(j)
                0990         uptakeNO2 = vmaxNO2(j)*MM_NO2*regQn*uptakeTempFunc*X(j)
                0991         uptakeNO3 = vmaxNO3(j)*MM_NO3*regQn*uptakeTempFunc*X(j)
8fbfd1f382 Oliv*0992 #ifdef DARWIN_ALLOW_FEQUOTA
                0993 #ifdef DARWIN_NITRATE_FELIMIT
                0994         uptakeNO3 = uptakeNO3 * limitfe
                0995 #endif
                0996 #endif
                0997         uptakeN = MAX(uptakeNH4 + uptakeNO2 + uptakeNO3,
6b11238516 Oliv*0998      &                vmaxN(j)*regQn*uptakeTempFunc*X(j)*diazo(j))
8fbfd1f382 Oliv*0999
                1000 c       linear limitation
086a45f245 Oliv*1001         limitn = MAX(0 _d 0, MIN(1 _d 0, (Qn(j) - Qnmin(j))/
                1002      &                                   (Qnmax(j) - Qnmin(j)) ))
8fbfd1f382 Oliv*1003 #else /* not DARWIN_ALLOW_NQUOTA */
                1004         inhibNH4 = EXP(-amminhib(j)*NH4)
                1005         limitNH4 = useNH4(j)*NH4/(NH4 + ksatNH4(j))
                1006         limitNO2 = useNO2(j)*NO2/
                1007      &   (NO2 + combNO(j)*(NO3 + ksatNO3(j) - ksatNO2(j)) + ksatNO2(j))*
                1008      &   inhibNH4
                1009         limitNO3 = useNO3(j)*NO3/
                1010      &   (combNO(j)*NO2 + NO3 + ksatNO3(j))*inhibNH4
                1011         limitn = limitNH4 + limitNO2 + limitNO3
c697368a52 Oliv*1012 C       normalize to sum 1
b854db6bc4 Oliv*1013         IF (limitn .GT. DARWIN_EPS) THEN
c697368a52 Oliv*1014           fracNH4 = limitNH4/limitn
                1015           fracNO2 = limitNO2/limitn
                1016           fracNO3 = limitNO3/limitn
                1017         ELSE
                1018           fracNH4 = 0 _d 0
                1019           fracNO2 = 0 _d 0
                1020           fracNO3 = 0 _d 0
                1021         ENDIF
c7b6c66d45 Oliv*1022 # ifdef DARWIN_ALLOW_CSTORE
                1023 C       compute maximum N uptake and C demand by N
                1024         uptakeNH4 = vmaxNH4(j)*limitNH4*uptakeTempFunc*X(j)
                1025         uptakeNO2 = vmaxNO2(j)*limitNO2*uptakeTempFunc*X(j)
                1026         uptakeNO3 = vmaxNO3(j)*limitNO3*uptakeTempFunc*X(j)
                1027         uptakeN = MAX(uptakeNH4 + uptakeNO2 + uptakeNO3,
                1028      &                VmaxN(j)*uptakeTempFunc*diazo(j)*X(j))
                1029         IF (R_NC(j) .GT. 0 _d 0) THEN
                1030           Dmd_NC = uptakeN/R_NC(j)
                1031           Dmd_min = MIN(Dmd_min, Dmd_NC)
                1032         ELSE
                1033           Dmd_NC = 0 _d 0
                1034         ENDIF
                1035 # endif
8fbfd1f382 Oliv*1036 C if diazo, all fracN* == 0 but want no N limitation
6b11238516 Oliv*1037         limitn = MIN(1 _d 0, limitn)
                1038         IF (diazo(j) .GT. 0) THEN
                1039           limitn = 1 _d 0
                1040         ENDIF
8fbfd1f382 Oliv*1041 #endif /* DARWIN_ALLOW_NQUOTA */
                1042
c7b6c66d45 Oliv*1043 #ifdef DARWIN_ALLOW_CSTORE
                1044 C       will apply limits for quota elements to growth rate instead
                1045         limitnut = 1 _d 0
                1046 #else
8fbfd1f382 Oliv*1047         limitnut = MIN(limitn, limitp, limitsi)
9eed995475 Oliv*1048 #if !(defined DARWIN_ALLOW_GEIDER && defined DARWIN_ALLOW_FEQUOTA)
8fbfd1f382 Oliv*1049         limitnut = MIN(limitnut, limitfe)
c7b6c66d45 Oliv*1050 #endif
8fbfd1f382 Oliv*1051 #endif
                1052
                1053         limitpCO2 = 1.
                1054
c7b6c66d45 Oliv*1055 C==== Photosynthesis ===================================================
8fbfd1f382 Oliv*1056 #ifdef DARWIN_ALLOW_GEIDER
                1057
                1058         alpha_I = 0 _d 0
                1059         DO l = 1, nlam
                1060           alpha_I = alpha_I + alphachl(j,l)*PAR(l)
                1061         ENDDO
                1062 C       NB: for quota, PCmax(j) = Vmax_c(j)
                1063         PCm = PCmax(j)*limitnut*photoTempFunc(j)*limitpCO2
                1064
c7b6c66d45 Oliv*1065 C       chl:c for 'balanced growth' at given light
8fbfd1f382 Oliv*1066         IF (PCm .GT. 0.0) THEN
                1067           acclim = MAX(chl2cmin(j), MIN(chl2cmax(j),
                1068      &             chl2cmax(j)/(1+(chl2cmax(j)*alpha_I)/(2*PCm)) ))
                1069         ELSE
                1070           acclim = chl2cmin(j)
                1071         ENDIF
                1072
                1073 #ifdef DARWIN_ALLOW_CHLQUOTA
                1074 C       quotas are already relative to carbon
                1075         chl2c = QChl(j)
                1076 #else
                1077         chl2c = acclim
                1078 #endif
                1079
                1080         alpha_I_growth = alpha_I
                1081 C a la quota
                1082 #ifdef DARWIN_ALLOW_FEQUOTA
                1083         alpha_I_growth = alpha_I_growth*limitfe
                1084 #endif
                1085
                1086 C       carbon-specific growth rate
                1087         IF (PCm .GT. 0.0 .AND. PARtot .GT. PARmin) THEN
fd20fa810c Oliv*1088           limitI = 1 _d 0 - EXP(-alpha_I_growth*chl2c/PCm)
8fbfd1f382 Oliv*1089         ELSE
fd20fa810c Oliv*1090           limitI = 0 _d 0
8fbfd1f382 Oliv*1091         ENDIF
                1092
                1093         IF (inhibGeider(j) .GT. 0.0) THEN
                1094 C         "total" PAR:
                1095           tmp = alpha_I/alpha_mean(j)
                1096           Ek = PCm/(chl2c*alpha_mean(j))
2cf0d99c0f Oliv*1097           IF (tmp .GT. Ek) THEN
                1098             EkoverE = Ek / tmp
fd20fa810c Oliv*1099             limitI = limitI*EkoverE*inhibGeider(j)
8fbfd1f382 Oliv*1100           ENDIF
                1101         ENDIF
                1102
fd20fa810c Oliv*1103         PC = PCm*limitI
                1104
8fbfd1f382 Oliv*1105 #else /* not DARWIN_ALLOW_GEIDER */
                1106
                1107         IF (PARtot .GT. PARmin) THEN
                1108 C         only 1 waveband without DARWIN_ALLOW_GEIDER
                1109           limitI = (1.0 _d 0 - EXP(-PARtot*ksatPAR(j)))*
                1110      &             EXP(-PARtot*kinhPAR(j)) * normI(j)
                1111           PC = PCmax(j)*limitnut*limitI*photoTempFunc(j)*limitpCO2
                1112         ELSE
fd20fa810c Oliv*1113           limitI = 0 _d 0
8fbfd1f382 Oliv*1114           PC = 0.0 _d 0
                1115         ENDIF
                1116         synthChl = 0.0
                1117
                1118 #endif /* DARWIN_ALLOW_GEIDER */
                1119
a092808e6b shlo*1120 #ifdef DARWIN_DIAG_PERTYPE
                1121 C  which nutrient is limiting growth?
                1122         IF (limitnut==limitn) diags(ilimN+j-1) = 1 _d 0
                1123         IF (limitnut==limitp) diags(ilimP+j-1) = 1 _d 0
                1124         IF (limitnut==limitfe) diags(ilimF+j-1) = 1 _d 0
                1125         IF (limitnut==limitsi) diags(ilimS+j-1) = 1 _d 0
                1126 C factors that go into growth equation
                1127         diags(ifIph+j-1) = limitI
                1128         fIph_sum = fIph_sum + limitI * X(j)
                1129         diags(ifTph+j-1) = photoTempFunc(j)
                1130         fTph_sum = fTph_sum + photoTempFunc(j) * X(j)
                1131         diags(ifnut+j-1) = limitnut
                1132         fnut_sum = fnut_sum + limitnut * X(j)
                1133 #endif
                1134
                1135 #endif /* DARWIN_MACROMOLECULAR_GROWTH */
                1136
c7b6c66d45 Oliv*1137         photoSyn = PC*X(j)
8fbfd1f382 Oliv*1138
8752612d6f Oliv*1139 C==== Respiration ======================================================
                1140 C Respiration goes straight back to DIC, so we subtract from DIC uptake.
                1141 C Uptake of other non-quota elements is reduced according to plankton
                1142 C stoichiometry to ensure conservation of chemical elements.
                1143 C Zooplankton respiration is handled in the mortality/exudation loop.
                1144         Xe = MAX(0 _d 0, X(j) - Xmin(j))
                1145         respir_c = respRate(j)*Xe*reminTempFunc
                1146         respir = respir + respir_c
c7b6c66d45 Oliv*1147         uptakeDIC = photoSyn - respir_c
fd20fa810c Oliv*1148 #ifdef DARWIN_DIAG_PERTYPE
                1149         diags(iResp+j-1) = respir_c
                1150 #endif
c7b6c66d45 Oliv*1151
                1152 C==== Growth ===========================================================
                1153 C Growth is growth rate of functional internal carbon
                1154
                1155 #ifdef DARWIN_ALLOW_CSTORE
                1156 C To decide the growth rate according to the minimum of C demand and
                1157 C photosynthesis.  To decide the minimum of C demand by the four
                1158 C elements (N, P, Si, Fe), Dmd_XC is the maximum carbon-specific
                1159 C photosynthesis rate convertible into biomass when X is the limiting
                1160 C nutrient as computed above.
                1161         growth = MIN(uptakeDIC, Dmd_min)
086a45f245 Oliv*1162         excretC = MAX(0 _d 0, uptakeDIC - growth)*FracExudeC(j)
c7b6c66d45 Oliv*1163         excC = excC + excretC
                1164         storeC = uptakeDIC - growth - excretC
                1165 #else /* without CSTORE */
                1166         growth = uptakeDIC
                1167 #endif /* DARWIN_ALLOW_CSTORE */
                1168
                1169 C==== Geider 1997 Chlorophyll ==========================================
8fbfd1f382 Oliv*1170 #ifdef DARWIN_ALLOW_GEIDER
76f7482618 Oliv*1171 # ifdef DARWIN_ALLOW_CHLQUOTA
                1172 #  ifdef DARWIN_ALLOW_NQUOTA
8fbfd1f382 Oliv*1173 C       Geider 1998
                1174         IF (PARtot .GT. PARmin) THEN
                1175          IF (alpha_I*chl2c .GT. 0.0 _d 0) THEN
c9c2f3a0a4 Oliv*1176           rhochl = Chl2Nmax*PC/(alpha_I*chl2c)
8fbfd1f382 Oliv*1177          ELSE
                1178           rhochl = Chl2Nmax
                1179          ENDIF
                1180          synthChl = rhochl*uptakeN
                1181         ELSE
                1182          synthChl = 0 _d 0
                1183         ENDIF
9de618f89c Oliv*1184 C       synth cost goes straight back to DIC
8fbfd1f382 Oliv*1185         uptakeDIC = uptakeDIC - synthcost*uptakeN
c7b6c66d45 Oliv*1186         growth = growth - synthcost*uptakeN
76f7482618 Oliv*1187 #  else /* not DARWIN_ALLOW_NQUOTA */
                1188 #   ifdef DARWIN_GEIDER_RHO_SYNTH
8fbfd1f382 Oliv*1189         IF (alpha_I .GT. 0.0 _d 0 .AND. acclim .GT. 0.0 _d 0) THEN
c7b6c66d45 Oliv*1190 C         should be chl2c instead of acclim?
c9c2f3a0a4 Oliv*1191           rhochl = Chl2Cmax(j)*PC/(alpha_I*acclim)
8fbfd1f382 Oliv*1192         ELSE
                1193           rhochl = 0.0 _d 0    ! should be Chl2Cmax(j) ?????
                1194         ENDIF
c7b6c66d45 Oliv*1195         synthChl = rhochl*growth +
8fbfd1f382 Oliv*1196      &                      acclimtimescl(j)*(acclim-chl2c)*X(j)
76f7482618 Oliv*1197 #   else
c7b6c66d45 Oliv*1198         synthChl = acclim*growth +
8fbfd1f382 Oliv*1199      &                      acclimtimescl(j)*(acclim-chl2c)*X(j)
76f7482618 Oliv*1200 #   endif
                1201 #  endif /* not DARWIN_ALLOW_NQUOTA */
                1202 # else /* not DARWIN_ALLOW_CHLQUOTA */
8fbfd1f382 Oliv*1203         chlout(j) = X(j)*Qc(j)*chl2c
                1204         synthChl = 0.0
76f7482618 Oliv*1205 # endif /* DARWIN_ALLOW_CHLQUOTA */
8fbfd1f382 Oliv*1206 #endif /* DARWIN_ALLOW_GEIDER */
9de618f89c Oliv*1207
c7b6c66d45 Oliv*1208 C==== Uptake of non-quota elements =====================================
                1209 C       non-quota elements are taken up with fixed stoichiometry
9de618f89c Oliv*1210 #ifndef DARWIN_ALLOW_NQUOTA
c7b6c66d45 Oliv*1211         uptakeN = growth*R_NC(j)
8752612d6f Oliv*1212 C Maintain ratios of NH4, NO2 and NO3 uptake.
                1213         IF (uptakeN .GE. 0 _d 0) THEN
                1214          uptakeNH4 = uptakeN*fracNH4
                1215          uptakeNO2 = uptakeN*fracNO2
                1216          uptakeNO3 = uptakeN*fracNO3
                1217         ELSE
                1218 C Do not allow uptake of NH4 and NO2 to go negative.
                1219 C Excess respired N goes to NO3.
                1220          uptakeNH4 = 0 _d 0
                1221          uptakeNO2 = 0 _d 0
                1222          uptakeNO3 = uptakeN
                1223         ENDIF
9de618f89c Oliv*1224 #endif
                1225 #ifndef DARWIN_ALLOW_PQUOTA
c7b6c66d45 Oliv*1226         uptakePO4 = growth*R_PC(j)
9de618f89c Oliv*1227 #endif
                1228 #ifndef DARWIN_ALLOW_SIQUOTA
c7b6c66d45 Oliv*1229         uptakeSiO2 = growth*R_SiC(j)
9de618f89c Oliv*1230 #endif
                1231 #ifndef DARWIN_ALLOW_FEQUOTA
c7b6c66d45 Oliv*1232         uptakeFeT = growth*R_FeC(j)
9de618f89c Oliv*1233 #endif
                1234
c7b6c66d45 Oliv*1235 C==== Nutrient consumption =============================================
                1236         consumDIC_PIC = consumDIC_PIC + growth*R_PICPOC(j)
8fbfd1f382 Oliv*1237         consumDIC  = consumDIC  + uptakeDIC
                1238         consumNH4  = consumNH4  + uptakeNH4
                1239         consumNO2  = consumNO2  + uptakeNO2
                1240         consumNO3  = consumNO3  + uptakeNO3
                1241         consumPO4  = consumPO4  + uptakePO4
                1242         consumSiO2 = consumSiO2 + uptakeSiO2
                1243         consumFeT  = consumFeT  + uptakeFeT
                1244
9de618f89c Oliv*1245         diags(iPP) = diags(iPP) + uptakeDIC
ee4178dfbe Oliv*1246 #ifdef DARWIN_DIAG_PERTYPE
9de618f89c Oliv*1247         diags(iPPplank+j-1) = diags(iPPplank+j-1) + uptakeDIC
d5f8f1c735 Oliv*1248         diags(iPCplank+j-1) = diags(iPCplank+j-1) + PC
ee4178dfbe Oliv*1249 #endif
8fbfd1f382 Oliv*1250         IF (diazo(j) .GT. 0.0 _d 0) THEN
                1251          diags(iNfix)=diags(iNfix)+uptakeN-uptakeNH4-uptakeNO2-uptakeNO3
                1252         ENDIF
c7b6c66d45 Oliv*1253 #ifdef DARWIN_ALLOW_CSTORE
                1254         diags(iEX) = diags(iEX) + excretC
                1255         diags(iGW) = diags(iGW) + growth
                1256         diags(iDN) = diags(iDN) + Dmd_NC
                1257         diags(iDP) = diags(iDP) + Dmd_PC
                1258         diags(iDFe) = diags(iDFe) + Dmd_FeC
                1259         diags(iDSi) = diags(iDSi) + Dmd_SiC
                1260         diags(iDmin)= diags(iDmin)+ Dmd_min
                1261 #ifdef DARWIN_ALLOW_CSTORE_DIAGS
                1262         diags(iEXplank+j-1) = diags(iEXplank+j-1) + excretC
                1263         diags(iGWplank+j-1) = diags(iGWplank+j-1) + growth
                1264         diags(iDNplank+j-1) = diags(iDNplank+j-1) + Dmd_NC
                1265         diags(iDPplank+j-1) = diags(iDPplank+j-1) + Dmd_PC
                1266         diags(iDFplank+j-1) = diags(iDFplank+j-1) + Dmd_FeC
                1267         diags(iDSplank+j-1) = diags(iDSplank+j-1) + Dmd_SiC
                1268         diags(iDminplank+j-1) = diags(iDminplank+j-1) + Dmd_min
                1269 #endif
                1270 #endif
                1271
                1272 C==== Apply phototrophy tendencies =====================================
                1273
                1274         gTr(ic+j-1)=gTr(ic+j-1)  + growth
a092808e6b shlo*1275 #ifdef DARWIN_MACROMOLECULAR_GROWTH
                1276         gTr(ic+j-1)=gTr(ic+j-1)  - mmX*exQc
                1277         gTr(iDOC) = gTr(iDOC) + mmX*exQc
                1278 #endif
c7b6c66d45 Oliv*1279 #ifdef DARWIN_ALLOW_CSTORE
                1280         gTr(ich+j-1)=gTr(ich+j-1)  + storeC
                1281 #endif
8fbfd1f382 Oliv*1282 #ifdef DARWIN_ALLOW_NQUOTA
                1283         gTr(in+j-1)=gTr(in+j-1)  + uptakeN
                1284 #endif
                1285 #ifdef DARWIN_ALLOW_PQUOTA
                1286         gTr(ip+j-1)=gTr(ip+j-1)  + uptakePO4
                1287 #endif
                1288 #ifdef DARWIN_ALLOW_SIQUOTA
                1289         gTr(isi+j-1)=gTr(isi+j-1) + uptakeSiO2
                1290 #endif
                1291 #ifdef DARWIN_ALLOW_FEQUOTA
                1292         gTr(ife+j-1)=gTr(ife+j-1) + uptakeFeT
                1293 #endif
                1294 #ifdef DARWIN_ALLOW_CHLQUOTA
                1295         gTr(iChl+j-1)=gTr(iChl+j-1) + synthChl
                1296 #endif
                1297
                1298 #ifdef DARWIN_DEBUG
6559acf5ff Oliv*1299         IF (k_debug.GT.0) THEN
                1300          print*,'uptake',myiter,k_debug,j,
8fbfd1f382 Oliv*1301      &     uptakeDIC,
                1302      &     uptakeNH4,
                1303      &     uptakeNO2,
                1304      &     uptakeNO3,
                1305      &     uptakeN,
                1306      &     uptakePO4,
                1307      &     uptakeSiO2,
                1308      &     uptakeFeT
                1309         ENDIF
                1310 #endif
                1311
b6fafcd98a Oliv*1312         X_photo = X_photo + X(j)
                1313
8fbfd1f382 Oliv*1314 C      isPhoto(j)
                1315        ENDIF
                1316 C     j
                1317       ENDDO
                1318
b6fafcd98a Oliv*1319       IF (X_photo .GT. 0 _d 0) THEN
                1320        fIph_sum = fIph_sum / X_photo
                1321        fTph_sum = fTph_sum / X_photo
                1322        fnut_sum = fnut_sum / X_photo
                1323       ELSE
                1324        fIph_sum = 0 _d 0
                1325        fTph_sum = 0 _d 0
                1326        fnut_sum = 0 _d 0
                1327       ENDIF
                1328       diags(ifIphavg) = fIph_sum
                1329       diags(ifTphavg) = fTph_sum
                1330       diags(ifnutavg) = fnut_sum
a092808e6b shlo*1331
8fbfd1f382 Oliv*1332 C=======================================================================
c7b6c66d45 Oliv*1333 C==== Bacteria =========================================================
8fbfd1f382 Oliv*1334
                1335       DO j = 1, nplank
                1336        IF (bactType(j) .NE. 0) THEN
                1337
                1338         uptakeO2  = 0. _d 0
                1339         uptakeNO3 = 0. _d 0
                1340         uptakePOC = 0. _d 0
                1341         uptakePON = 0. _d 0
                1342         uptakePOP = 0. _d 0
                1343         uptakePOFe = 0. _d 0
                1344         uptakeDOC = 0. _d 0
                1345         uptakeDON = 0. _d 0
                1346         uptakeDOP = 0. _d 0
                1347         uptakeDOFe = 0. _d 0
                1348         hydrolPOC = 0. _d 0
                1349         hydrolPON = 0. _d 0
                1350         hydrolPOP = 0. _d 0
                1351         hydrolPOFe = 0. _d 0
                1352         respPOC = 0. _d 0
                1353         respPON = 0. _d 0
                1354         respPOP = 0. _d 0
                1355         respPOFe = 0. _d 0
                1356         respDOC = 0. _d 0
                1357         respDON = 0. _d 0
                1358         respDOP = 0. _d 0
                1359         respDOFe = 0. _d 0
                1360         growth = 0. _d 0
                1361
                1362         IF (isAerobic(j) .NE. 0) THEN
                1363           muO = yieldO2(j)*pcoefO2*O2
                1364         ELSEIF (isDenit(j) .NE. 0) THEN
3f05298721 Oliv*1365           muO = yieldNO3(j)*pmaxDIN*NO3/(NO3 + ksatDIN)*hetTempFunc(j)
8fbfd1f382 Oliv*1366         ENDIF
                1367
                1368 C       POM-consuming (particle-associated)
                1369         IF (bactType(j) .EQ. 1) THEN
                1370
3f05298721 Oliv*1371           PCm = yield(j)*PCmax(j)*hetTempFunc(j)
8fbfd1f382 Oliv*1372           muPON  = PCm*PON/(PON + ksatPON(j))
                1373           muPOC  = PCm*POC/(POC + ksatPOC(j))
                1374           muPOP  = PCm*POP/(POP + ksatPOP(j))
                1375           muPOFe = PCm*POFe/(POFe + ksatPOFe(j))
                1376           mu = MIN(muPON, muPOC, muPOP, muPOFe, muO)
                1377
                1378           growth = mu*X(j)
                1379
                1380           uptakePOC = alpha_hydrol*growth/yield(j)
                1381           uptakePON  = uptakePOC*R_NC(j)
                1382           uptakePOP  = uptakePOC*R_PC(j)
                1383           uptakePOFe = uptakePOC*R_FeC(j)
                1384 C         O2/NO3 is only used for the part of POC that is metabolized:
                1385           uptakeO2 = isAerobic(j)*growth/yieldO2(j)
                1386           uptakeNO3 = isDenit(j)*growth/yieldNO3(j)
                1387
                1388 C         This is the part of POM that is hydrolized into DOM:
                1389           hydrolPOC = (alpha_hydrol-1)*growth/yield(j)
                1390           hydrolPON  = hydrolPOC*R_NC(j)
                1391           hydrolPOP  = hydrolPOC*R_PC(j)
                1392           hydrolPOFe = hydrolPOC*R_FeC(j)
                1393
                1394 C         These are the bacteria products for remineralization of POM:
                1395           respPOC = growth*(1/yield(j)-1)
                1396           respPON  = respPOC*R_NC(j)
                1397           respPOP  = respPOC*R_PC(j)
                1398           respPOFe = respPOC*R_FeC(j)
                1399
                1400 C       DOM-consuming (free-living):
                1401         ELSEIF (bactType(j) .EQ. 2) THEN
                1402
3f05298721 Oliv*1403           PCm = yield(j)*PCmax(j)*hetTempFunc(j)
8fbfd1f382 Oliv*1404           muDON  = PCm*DON/(DON + ksatDON(j))
                1405           muDOC  = PCm*DOC/(DOC + ksatDOC(j))
                1406           muDOP  = PCm*DOP/(DOP + ksatDOP(j))
                1407           muDOFe = PCm*DOFe/(DOFe + ksatDOFe(j))
                1408           mu = MIN(muDON, muDOC, muDOP, muDOFe, muO)
                1409
                1410           growth = mu*X(j)
                1411
                1412           uptakeDOC = growth/yield(j)
                1413           uptakeDON  = uptakeDOC*R_NC(j)
                1414           uptakeDOP  = uptakeDOC*R_PC(j)
                1415           uptakeDOFe = uptakeDOC*R_FeC(j)
                1416           uptakeO2 = isAerobic(j)*growth/yieldO2(j)
                1417           uptakeNO3 = isDenit(j)*growth/yieldNO3(j)
                1418
                1419 C         DOC respired to DIC
                1420           respDOC = growth*(1/yield(j)-1)
                1421           respDON  = respDOC*R_NC(j)
                1422           respDOP  = respDOC*R_PC(j)
                1423           respDOFe = respDOC*R_FeC(j)
                1424
                1425         ENDIF
                1426
5e411acc9e Oliv*1427         diags(iDenitN) = diags(iDenitN) + uptakeNO3
ee4178dfbe Oliv*1428 #ifdef DARWIN_DIAG_PERTYPE
a93256cee4 Oliv*1429         diags(iHPplank+j-1) = diags(iHPplank+j-1) + growth
                1430         diags(iHCplank+j-1) = diags(iHCplank+j-1) + mu
ee4178dfbe Oliv*1431 #endif
8fbfd1f382 Oliv*1432
90de433d90 Oliv*1433         gTr(ic+j-1) = gTr(ic+j-1) + growth
                1434
                1435 C bacteria have fixed elemental ratios for now
                1436 #ifdef DARWIN_ALLOW_NQUOTA
                1437         gTr(in+j-1) = gTr(in+j-1) + growth*R_NC(j)
                1438 #endif
                1439 #ifdef DARWIN_ALLOW_PQUOTA
                1440         gTr(ip+j-1) = gTr(ip+j-1) + growth*R_PC(j)
                1441 #endif
                1442 #ifdef DARWIN_ALLOW_FEQUOTA
                1443         gTr(ife+j-1) = gTr(ife+j-1) + growth*R_FeC(j)
                1444 #endif
8fbfd1f382 Oliv*1445
                1446 C==== Cumulative consum, remin, and prod ===============================
                1447         consumNO3  = consumNO3  + uptakeNO3
                1448
                1449 C       add B consum and accumulating remin, and prod:
                1450         consumO2 = consumO2 + uptakeO2
                1451
                1452         consumDOC = consumDOC + uptakeDOC
                1453         consumDON = consumDON + uptakeDON
                1454         consumDOP = consumDOP + uptakeDOP
                1455         consumDOFe = consumDOFe + uptakeDOFe
                1456
                1457         consumPOC = consumPOC + uptakePOC
                1458         consumPON = consumPON + uptakePON
                1459         consumPOP = consumPOP + uptakePOP
                1460         consumPOFe = consumPOFe + uptakePOFe
                1461
                1462         reminPOC = reminPOC + respPOC
                1463         reminPON = reminPON + respPON
                1464         reminPOP = reminPOP + respPOP
                1465         reminPOFe = reminPOFe + respPOFe
                1466
                1467         solubilPOC = solubilPOC + hydrolPOC
                1468         solubilPON = solubilPON + hydrolPON
                1469         solubilPOP = solubilPOP + hydrolPOP
                1470         solubilPOFe = solubilPOFe + hydrolPOFe
                1471
                1472         reminDOC = reminDOC + respDOC
                1473         reminDON = reminDON + respDON
                1474         reminDOP = reminDOP + respDOP
                1475         reminDOFe = reminDOFe + respDOFe
                1476
                1477        ENDIF
c74953b924 Oliv*1478 C     bacteria j loop end
8fbfd1f382 Oliv*1479       ENDDO
                1480
                1481 C=======================================================================
c7b6c66d45 Oliv*1482 C==== Apply phototrophic and bacterial nutrient consumption ============
8fbfd1f382 Oliv*1483
                1484       gTr(iDIC )=gTr(iDIC ) - consumDIC - consumDIC_PIC
b61dca872c Oliv*1485
                1486       diags(iConsDIC) = consumDIC
                1487       diags(iConsDIC_PIC) = consumDIC_PIC
                1488
8fbfd1f382 Oliv*1489       gTr(iNH4 )=gTr(iNH4 ) - consumNH4
                1490       gTr(iNO2 )=gTr(iNO2 ) - consumNO2
                1491       gTr(iNO3 )=gTr(iNO3 ) - consumNO3
                1492       gTr(iPO4 )=gTr(iPO4 ) - consumPO4
                1493       gTr(iSiO2)=gTr(iSiO2) - consumSiO2
                1494       gTr(iFeT )=gTr(iFeT ) - consumFeT
                1495
c7b6c66d45 Oliv*1496 C==== Add parameterized remineralization ===============================
eff5247c2f Oliv*1497 #ifdef DARWIN_ALLOW_DENIT
                1498       IF (O2 .GE. O2crit .OR. NO3 .GE. NO3crit) THEN
                1499 #else
                1500       IF (.TRUE.) THEN
                1501 #endif
a092808e6b shlo*1502       IF (ksatO2remin .GT. 0 _d 0) THEN
                1503        tmp = O2/(O2 + ksatO2remin)
                1504       ELSE
                1505        tmp = 1 _d 0
                1506       ENDIF
8fbfd1f382 Oliv*1507 C parameterized remineralization; want to set all K except KPOSi to zero
                1508 C if running with bacteria
a092808e6b shlo*1509 C shlomit: add O2-dependece respiration
                1510       respDOC  = reminTempFunc*KDOC *DOC*tmp
                1511       respDON  = reminTempFunc*KDON *DON*tmp
                1512       respDOP  = reminTempFunc*KDOP *DOP*tmp
                1513       respDOFe = reminTempFunc*KDOFe*DOFe*tmp
                1514       respPOC  = reminTempFunc*KPOC *POC*tmp
                1515       respPON  = reminTempFunc*KPON *PON*tmp
                1516       respPOP  = reminTempFunc*KPOP *POP*tmp
                1517       respPOSi = reminTempFunc*KPOSi*POSi*tmp
                1518       respPOFe = reminTempFunc*KPOFe*POFe*tmp
8fbfd1f382 Oliv*1519
                1520       consumDOC  = consumDOC  + respDOC
                1521       consumDON  = consumDON  + respDON
                1522       consumDOP  = consumDOP  + respDOP
                1523       consumDOFe = consumDOFe + respDOFe
                1524       consumPOC  = consumPOC  + respPOC
                1525       consumPON  = consumPON  + respPON
                1526       consumPOP  = consumPOP  + respPOP
                1527       consumPOSi = consumPOSi + respPOSi
                1528       consumPOFe = consumPOFe + respPOFe
                1529
                1530       reminDOC  = reminDOC  + respDOC
                1531       reminDON  = reminDON  + respDON
                1532       reminDOP  = reminDOP  + respDOP
                1533       reminDOFe = reminDOFe + respDOFe
                1534       reminPOC  = reminPOC  + respPOC
                1535       reminPON  = reminPON  + respPON
                1536       reminPOP  = reminPOP  + respPOP
                1537       reminPOSi = reminPOSi + respPOSi
                1538       reminPOFe = reminPOFe + respPOFe
eff5247c2f Oliv*1539 C     endif O2, NO3 for DARWIN_ALLOW_DENIT
                1540       ENDIF
8fbfd1f382 Oliv*1541
                1542 #ifdef DARWIN_ALLOW_CARBON
6100320f2e Oliv*1543       IF (O2 .GE. O2crit) THEN
eff5247c2f Oliv*1544         consumO2  = consumO2  + respDOP*R_OP
8fbfd1f382 Oliv*1545 #ifndef DARWIN_ALLOW_CDOM
eff5247c2f Oliv*1546         consumO2  = consumO2  + respPOP*R_OP
8fbfd1f382 Oliv*1547 #endif
eff5247c2f Oliv*1548       ENDIF
8fbfd1f382 Oliv*1549 #endif
                1550
b5b44fea51 Oliv*1551       IF (darwin_disscSelect .EQ. 2) THEN
                1552 C      Naviaux et al. 2019, Marine Chemistry dissolution rate law
                1553        IF (omegaCl .LT. 1.0 _d 0) THEN
                1554         IF (omegaCl .GT. 0.8272 _d 0) THEN
                1555           disscPIC = PIC * 5.22 _d -9 * (1 _d 0 - omegaCl)**0.11 _d 0
                1556         ELSE
                1557           disscPIC = PIC * 1.65 _d -5 * (1 _d 0 - omegaCl)**4.7 _d 0
                1558         ENDIF
                1559        ELSE
                1560         disscPIC = 0.0 _d 0
                1561        ENDIF
                1562       ELSEIF (darwin_disscSelect .EQ. 1) THEN
ba0b6d5d33 Oliv*1563 C      Keir 1980
                1564        IF (omegaCl .LT. 1.0 _d 0) THEN
                1565         disscPIC = PIC*darwin_KeirCoeff
                1566      &                *(1 _d 0 - omegaCl)**darwin_KeirExp
                1567        ELSE
                1568         disscPIC = 0 _d 0
                1569        ENDIF
                1570       ELSE
                1571        disscPIC = Kdissc*PIC
                1572       ENDIF
8fbfd1f382 Oliv*1573
c7b6c66d45 Oliv*1574 C==== Nitrogen chemistry ===============================================
8fbfd1f382 Oliv*1575 c NH4 -> NO2 -> NO3 by bacterial action, parameterized
                1576       prodNO2 = knita*NH4
                1577       prodNO3 = knitb*NO2
                1578       IF (PAR_oxi .NE. 0.0 _d 0) THEN
086a45f245 Oliv*1579         prodNO2 = prodNO2*MAX(0 _d 0, 1.0 - PARtot/PAR_oxi)
                1580         prodNO3 = prodNO3*MAX(0 _d 0, 1.0 - PARtot/PAR_oxi)
8fbfd1f382 Oliv*1581       ENDIF
                1582
c7b6c66d45 Oliv*1583 C==== CDOM =============================================================
8fbfd1f382 Oliv*1584 #ifdef DARWIN_ALLOW_CDOM
1c72adde5c Oliv*1585 # ifdef DARWIN_CDOM_UNITS_CARBON
                1586       reminPOC_CDOM = fracCDOM*reminPOC
                1587       reminPOP_CDOM = R_PC_CDOM*reminPOC_CDOM
                1588       reminPON_CDOM = R_NC_CDOM*reminPOC_CDOM
                1589       reminPOFe_CDOM = R_FeC_CDOM*reminPOC_CDOM
                1590 # else
8fbfd1f382 Oliv*1591       reminPOP_CDOM = fracCDOM*reminPOP
                1592       reminPOC_CDOM = R_CP_CDOM*reminPOP_CDOM
                1593       reminPON_CDOM = R_NP_CDOM*reminPOP_CDOM
                1594       reminPOFe_CDOM = R_FeP_CDOM*reminPOP_CDOM
1c72adde5c Oliv*1595 # endif
eff5247c2f Oliv*1596 c degradation of  CDOM
                1597       tmp = CDOMdegrd
                1598 #ifdef DARWIN_ALLOW_DENIT
                1599       IF (O2 .LT. O2crit .AND. NO3 .LT. NO3crit) THEN
                1600         tmp = 0 _d 0
                1601       ENDIF
                1602 #endif
                1603 c increase when bleached by light
                1604       tmp = tmp + CDOMbleach*MIN(1.0 _d 0, PARtot/PARCDOM)
1c72adde5c Oliv*1605 # ifdef DARWIN_CDOM_UNITS_CARBON
                1606       degrCDOM_DOC = reminTempFunc*CDOM*tmp
                1607       degrCDOM_DOP  = R_PC_CDOM  * degrCDOM_DOC
                1608       degrCDOM_DON  = R_NC_CDOM  * degrCDOM_DOC
                1609       degrCDOM_DOFe = R_FeC_CDOM * degrCDOM_DOC
                1610 # else
eff5247c2f Oliv*1611       degrCDOM_DOP = reminTempFunc*CDOM*tmp
8fbfd1f382 Oliv*1612       degrCDOM_DOC  = R_CP_CDOM  * degrCDOM_DOP
                1613       degrCDOM_DON  = R_NP_CDOM  * degrCDOM_DOP
                1614       degrCDOM_DOFe = R_FeP_CDOM * degrCDOM_DOP
1c72adde5c Oliv*1615 # endif
8fbfd1f382 Oliv*1616 #endif
                1617
c7b6c66d45 Oliv*1618 C==== Apply remaining tendencies =======================================
8fbfd1f382 Oliv*1619
                1620 #ifdef DARWIN_ALLOW_CARBON
                1621 c production of O2 by photosynthesis
                1622       gTr(iO2  )=gTr(iO2  ) + R_OP*consumPO4
5e411acc9e Oliv*1623       diags(iProdO2) = diags(iProdO2) + R_OP*consumPO4
8fbfd1f382 Oliv*1624 c loss of O2 by remineralization
6100320f2e Oliv*1625       IF (O2 .GE. O2crit) THEN
8fbfd1f382 Oliv*1626         gTr(iO2)=gTr(iO2) - consumO2
5e411acc9e Oliv*1627         diags(iConsO2) = diags(iConsO2) + consumO2
8fbfd1f382 Oliv*1628       ENDIF
                1629
                1630       gTr(iALK)=gTr(iALK) - (prodNO3 - consumNO3)
                1631      &                    - 2.0 _d 0*(consumDIC_PIC - disscPIC)
5e411acc9e Oliv*1632       diags(iSrcAlk) = diags(iSrcAlk) - (prodNO3 - consumNO3)
                1633      &                    - 2.0 _d 0*(consumDIC_PIC - disscPIC)
8fbfd1f382 Oliv*1634 #endif /* DARWIN_ALLOW_CARBON */
                1635
                1636       gTr(iDIC )=gTr(iDIC ) + reminDOC + disscPIC
b61dca872c Oliv*1637
                1638       diags(iReminDIC_DOC) = reminDOC
                1639       diags(iDisscDIC_PIC) = disscPIC
                1640
8fbfd1f382 Oliv*1641       gTr(iNH4 )=gTr(iNH4 ) + reminDON - prodNO2
                1642       gTr(iNO2 )=gTr(iNO2 ) + prodNO2 - prodNO3
                1643       gTr(iNO3 )=gTr(iNO3 ) + prodNO3
                1644 #ifdef DARWIN_ALLOW_DENIT
                1645       IF (O2 .LT. O2crit) THEN
                1646         denitNH4 = reminDON
                1647         denit = denit_NP*reminDOP
                1648 #ifndef DARWIN_ALLOW_CDOM
                1649         denitNH4 = denitNH4 + reminPON
                1650         denit = denit + denit_NP*reminPOP
                1651 #endif
                1652         diags(iDenit) = denit
                1653         gTr(iNH4)=gTr(iNH4) - denitNH4
                1654         gTr(iNO3)=gTr(iNO3) - denit_NO3/denit_np*denit
                1655         gTr(iALK)=gTr(iALK) + denit_NO3/denit_np*denit
5e411acc9e Oliv*1656         diags(iDenitN)=diags(iDenitN)+denitNH4+denit_NO3/denit_np*denit
                1657         diags(iSrcAlk)=diags(iSrcAlk)+denit_NO3/denit_np*denit
8fbfd1f382 Oliv*1658       ENDIF
                1659 #endif /* DARWIN_ALLOW_DENIT */
                1660
                1661       gTr(iPO4 )=gTr(iPO4 ) + reminDOP
                1662       gTr(iFeT )=gTr(iFeT ) + reminDOFe
                1663       gTr(iSiO2)=gTr(iSiO2)             + reminPOSi
                1664
                1665 C     DOC is created by #4 PA-assoc solubilization and consumed by #5
                1666       gTr(iDOC )=gTr(iDOC ) + solubilPOC - consumDOC
c7b6c66d45 Oliv*1667 #ifdef DARWIN_ALLOW_CSTORE
                1668       gTr(iDOC )=gTr(iDOC ) + excC
                1669 #endif
8fbfd1f382 Oliv*1670       gTr(iDON )=gTr(iDON ) + solubilPON - consumDON
                1671       gTr(iDOP )=gTr(iDOP ) + solubilPOP - consumDOP
                1672       gTr(iDOFe)=gTr(iDOFe) + solubilPOFe - consumDOFe
                1673
                1674       gTr(iPIC )=gTr(iPIC ) - disscPIC
                1675       gTr(iPOC )=gTr(iPOC ) - consumPOC
                1676       gTr(iPON )=gTr(iPON ) - consumPON
                1677       gTr(iPOP )=gTr(iPOP ) - consumPOP
                1678       gTr(iPOFe)=gTr(iPOFe) - consumPOFe
                1679       gTr(iPOSi)=gTr(iPOSi) - consumPOSi
                1680
                1681 #ifdef DARWIN_ALLOW_CDOM
                1682       gTr(iDOC )=gTr(iDOC ) + reminPOC  - reminPOC_CDOM  + degrCDOM_DOC
                1683       gTr(iDON )=gTr(iDON ) + reminPON  - reminPON_CDOM  + degrCDOM_DON
                1684       gTr(iDOP )=gTr(iDOP ) + reminPOP  - reminPOP_CDOM  + degrCDOM_DOP
                1685       gTr(iDOFe)=gTr(iDOFe) + reminPOFe - reminPOFe_CDOM + degrCDOM_DOFe
                1686
1c72adde5c Oliv*1687 # ifdef DARWIN_CDOM_UNITS_CARBON
                1688       gTr(iCDOM)=gTr(iCDOM) + reminPOC_CDOM - degrCDOM_DOC
                1689 # else
8fbfd1f382 Oliv*1690       gTr(iCDOM)=gTr(iCDOM) + reminPOP_CDOM - degrCDOM_DOP
1c72adde5c Oliv*1691 # endif
8fbfd1f382 Oliv*1692 #else
                1693       gTr(iDIC )=gTr(iDIC ) + reminPOC
b61dca872c Oliv*1694
                1695       diags(iReminDIC_POC) = reminPOC
                1696
8fbfd1f382 Oliv*1697       gTr(iNH4 )=gTr(iNH4 ) + reminPON
                1698       gTr(iPO4 )=gTr(iPO4 ) + reminPOP
                1699       gTr(iFeT )=gTr(iFeT ) + reminPOFe
                1700 #endif /* DARWIN_ALLOW_CDOM */
                1701
b61dca872c Oliv*1702       diags(iConsALK) = (prodNO3 - consumNO3)
                1703      & +2.0 _d 0*(consumDIC_PIC - disscPIC)
8fbfd1f382 Oliv*1704       diags(iConsDIN) = consumNH4 + consumNO2 + consumNO3
b61dca872c Oliv*1705       diags(iConsNO3) = consumNO3
                1706       diags(iConsNO2) = consumNO2
                1707       diags(iConsNH4) = consumNH4
8fbfd1f382 Oliv*1708       diags(iConsPO4) = consumPO4
                1709       diags(iConsSi)  = consumSiO2
                1710       diags(iConsFe)  = consumFeT
                1711
c7b6c66d45 Oliv*1712 C=======================================================================
                1713 C==== Grazing ==========================================================
8fbfd1f382 Oliv*1714
                1715       DO j=1,nplank
                1716        preygraz(j)   = 0.0
                1717        preygrazexp(j) = 0.0
a93256cee4 Oliv*1718        grazrate(j)  = 0.0
8fbfd1f382 Oliv*1719        predgrazc(j)  = 0.0
5e7acb36b1 daat*1720        sumprey(j) = 0.0
8fbfd1f382 Oliv*1721 #ifdef DARWIN_ALLOW_NQUOTA
                1722        predgrazn(j)  = 0.0
                1723 #endif
                1724 #ifdef DARWIN_ALLOW_PQUOTA
                1725        predgrazp(j)  = 0.0
                1726 #endif
                1727 #ifdef DARWIN_ALLOW_FEQUOTA
                1728        predgrazfe(j) = 0.0
                1729 #endif
                1730       ENDDO
                1731       graz2POC  = 0.0
                1732       graz2PON  = 0.0
                1733       graz2POP  = 0.0
                1734       graz2POSI = 0.0
                1735       graz2POFE = 0.0
                1736       graz2OC   = 0.0
                1737       graz2ON   = 0.0
                1738       graz2OP   = 0.0
                1739       graz2OFE  = 0.0
                1740       graz2PIC  = 0.0
                1741
                1742       regQn  = 1.0
                1743       regQp  = 1.0
                1744       regQfe = 1.0
                1745       regQc  = 1.0
                1746
                1747 C=======================================================================
                1748       DO jz = 1, nplank
                1749        IF (isPred(jz).NE.0) THEN
                1750
                1751 C       regulate grazing near full quota
                1752         regQc = 1.0 _d 0
                1753 #ifdef DARWIN_ALLOW_NQUOTA
086a45f245 Oliv*1754         regQn = MAX(0 _d 0, MIN(1 _d 0, (Qnmax(jz)-Qn(jz))/
                1755      &                                  (Qnmax(jz)-Qnmin(jz)) ))
8fbfd1f382 Oliv*1756         regQc = MIN(regQc, 1.0 _d 0 - regQn)
                1757         regQn = regQn**hillnumGraz
                1758 #endif
                1759 #ifdef DARWIN_ALLOW_PQUOTA
086a45f245 Oliv*1760         regQp = MAX(0 _d 0, MIN(1 _d 0, (Qpmax(jz)-Qp(jz))/
                1761      &                                  (Qpmax(jz)-Qpmin(jz)) ))
8fbfd1f382 Oliv*1762         regQc = MIN(regQc, 1.0 _d 0 - regQp)
                1763         regQp = regQp**hillnumGraz
                1764 #endif
                1765 #ifdef DARWIN_ALLOW_FEQUOTA
086a45f245 Oliv*1766         regQfe= MAX(0 _d 0, MIN(1 _d 0, (Qfemax(jz)-Qfe(jz))/
                1767      &                                  (Qfemax(jz)-Qfemin(jz)) ))
8fbfd1f382 Oliv*1768         regQc = MIN(regQc, 1.0 _d 0 - regQfe)
                1769         regQfe=regQfe**hillnumGraz
                1770 #endif
                1771         regQc = regQc**hillnumGraz
                1772
                1773         sumpref = 0.0
                1774         DO jp = 1, nplank
                1775         IF (palat(jp,jz).NE.0 _d 0) THEN
5e7acb36b1 daat*1776          sumprey(jz) = sumprey(jz) + palat(jp,jz)*X(jp)
8fbfd1f382 Oliv*1777 #ifdef DARWIN_GRAZING_SWITCH
                1778          sumpref = sumpref + palat(jp,jz)*palat(jp,jz)*X(jp)*X(jp)
                1779 #else
                1780          sumpref = sumpref + palat(jp,jz)*X(jp)
                1781 #endif
                1782         ENDIF
                1783         ENDDO
5e7acb36b1 daat*1784         sumprey(jz) = MAX(0 _d 0, sumprey(jz) - phygrazmin)
8fbfd1f382 Oliv*1785         sumpref = MAX(phygrazmin, sumpref)
a93256cee4 Oliv*1786         tmp = grazemax(jz)*grazTempFunc(jz)**tempGraz(jz)*
5e7acb36b1 daat*1787      &    (sumprey(jz)**hollexp/
                1788      &     (sumprey(jz)**hollexp + kgrazesat(jz)**hollexp))*
                1789      &    (1.0 - EXP(-inhib_graz*sumprey(jz)))**inhib_graz_exp
8fbfd1f382 Oliv*1790
                1791         predexpc  = 0.0 _d 0
                1792         predexpn  = 0.0 _d 0
                1793         predexpp  = 0.0 _d 0
                1794         predexpfe = 0.0 _d 0
                1795         DO jp = 1, nplank
                1796          IF (palat(jp,jz).NE.0 _d 0) THEN
                1797 #ifdef DARWIN_GRAZING_SWITCH
                1798           grazphy = tmp*palat(jp,jz)*palat(jp,jz)*X(jp)*X(jp)/sumpref
                1799 #else
                1800           grazphy = tmp*palat(jp,jz)*X(jp)/sumpref
                1801 #endif
a93256cee4 Oliv*1802           grazrate(jz) = grazrate(jz) + grazphy*asseff(jp,jz)*regQc
                1803           grazphy = grazphy*X(jz)
8fbfd1f382 Oliv*1804
                1805           expFrac = ExportFracPreyPred(jp,jz)
                1806
                1807           preygraz(jp) = preygraz(jp) + grazphy
                1808           preygrazexp(jp) = preygrazexp(jp) + expFrac*grazphy
                1809
c7b6c66d45 Oliv*1810           predgrazc(jz) = predgrazc(jz) + grazphy*asseff(jp,jz)*
                1811      &                                    regQc*Qctot(jp)
                1812           predexpc = predexpc + expFrac*grazphy*asseff(jp,jz)*
                1813      &                                    regQc*Qctot(jp)
8fbfd1f382 Oliv*1814 #ifdef DARWIN_ALLOW_NQUOTA
                1815           predgrazn(jz) = predgrazn(jz) + grazphy*asseff(jp,jz)*
                1816      &                                    regQn*Qn(jp)
                1817           predexpn = predexpn + expFrac*grazphy*asseff(jp,jz)*
                1818      &                                    regQn*Qn(jp)
                1819 #endif
                1820 #ifdef DARWIN_ALLOW_PQUOTA
                1821           predgrazp(jz) = predgrazp(jz) + grazphy*asseff(jp,jz)*
                1822      &                                    regQp*Qp(jp)
                1823           predexpp = predexpp + expFrac*grazphy*asseff(jp,jz)*
                1824      &                                  regQp*Qp(jp)
                1825 #endif
                1826 #ifdef DARWIN_ALLOW_FEQUOTA
                1827           predgrazfe(jz) = predgrazfe(jz) + grazphy*asseff(jp,jz)*
                1828      &                                      regQfe*Qfe(jp)
                1829           predexpfe = predexpfe + expFrac*grazphy*asseff(jp,jz)*
                1830      &                                    regQfe*Qfe(jp)
                1831 #endif
                1832          ENDIF
                1833         ENDDO
                1834
                1835 C organic-matter gain will be total preygraz - predgraz
                1836         graz2OC   = graz2OC  - predgrazc(jz)
                1837         graz2POC  = graz2POC - predexpc
                1838
                1839 #ifdef DARWIN_ALLOW_NQUOTA
                1840         graz2ON   = graz2ON  - predgrazn(jz)
                1841         graz2PON  = graz2PON - predexpn
                1842 #else
                1843         graz2ON   = graz2ON  - predgrazc(jz)*Qn(jz)
                1844         graz2PON  = graz2PON - predexpc     *Qn(jz)
                1845 #endif
                1846
                1847 #ifdef DARWIN_ALLOW_PQUOTA
                1848         graz2OP   = graz2OP  - predgrazp(jz)
                1849         graz2POP  = graz2POP - predexpp
                1850 #else
                1851         graz2OP   = graz2OP  - predgrazc(jz)*Qp(jz)
                1852         graz2POP  = graz2POP - predexpc     *Qp(jz)
                1853 #endif
                1854
                1855 #ifdef DARWIN_ALLOW_FEQUOTA
                1856         graz2OFe  = graz2OFe  - predgrazfe(jz)
                1857         graz2POFe = graz2POFe - predexpfe
                1858 #else
                1859         graz2OFe  = graz2OFe  - predgrazc(jz)*Qfe(jz)
                1860         graz2POFe = graz2POFe - predexpc     *Qfe(jz)
                1861 #endif
                1862
                1863        ENDIF
                1864 C     end predator loop
                1865       ENDDO
                1866
                1867       DO jp = 1, nplank
                1868        IF (isPrey(jp).NE.0) THEN
c7b6c66d45 Oliv*1869           graz2OC  = graz2OC  + preygraz(jp)*Qctot(jp)
8fbfd1f382 Oliv*1870           graz2ON  = graz2ON  + preygraz(jp)*Qn (jp)
                1871           graz2OP  = graz2OP  + preygraz(jp)*Qp (jp)
                1872           graz2POSi = graz2POSi + preygraz(jp)*Qsi(jp)
                1873           graz2OFe = graz2OFe + preygraz(jp)*Qfe(jp)
                1874           graz2PIC = graz2PIC + preygraz(jp)*R_PICPOC(jp)
                1875
c7b6c66d45 Oliv*1876           graz2POC  = graz2POC  + preygrazexp(jp)*Qctot(jp)
8fbfd1f382 Oliv*1877           graz2PON  = graz2PON  + preygrazexp(jp)*Qn (jp)
                1878           graz2POP  = graz2POP  + preygrazexp(jp)*Qp (jp)
                1879           graz2POFe = graz2POFe + preygrazexp(jp)*Qfe(jp)
                1880        ENDIF
                1881       ENDDO
                1882
c7b6c66d45 Oliv*1883 C==== Apply grazing tendencies =========================================
8fbfd1f382 Oliv*1884
                1885       gTr(iDOC )=gTr(iDOC ) + graz2OC  - graz2POC
                1886       gTr(iDON )=gTr(iDON ) + graz2ON  - graz2PON
                1887       gTr(iDOP )=gTr(iDOP ) + graz2OP  - graz2POP
                1888       gTr(iDOFe)=gTr(iDOFe) + graz2OFe - graz2POFe
                1889       gTr(iPOC )=gTr(iPOC ) + graz2POC
                1890       gTr(iPON )=gTr(iPON ) + graz2PON
                1891       gTr(iPOP )=gTr(iPOP ) + graz2POP
                1892       gTr(iPOSi)=gTr(iPOSi) + graz2POSi
                1893       gTr(iPOFe)=gTr(iPOFe) + graz2POFe
                1894       gTr(iPIC )=gTr(iPIC ) + graz2PIC
                1895 #ifdef DARWIN_ALLOW_CDOM
1c72adde5c Oliv*1896 # ifdef DARWIN_CDOM_UNITS_CARBON
                1897       graz2CDOM = fracCDOM*(graz2OC - graz2POC)
                1898 # else
8fbfd1f382 Oliv*1899       graz2CDOM = fracCDOM*(graz2OP - graz2POP)
1c72adde5c Oliv*1900 # endif
8fbfd1f382 Oliv*1901       gTr(iCDOM)=gTr(iCDOM) + graz2CDOM
1c72adde5c Oliv*1902 # ifdef DARWIN_CDOM_UNITS_CARBON
                1903       gTr(iDOC )=gTr(iDOC ) -           graz2CDOM
                1904       gTr(iDON )=gTr(iDON ) - R_NC_CDOM*graz2CDOM
                1905       gTr(iDOP )=gTr(iDOP ) - R_PC_CDOM*graz2CDOM
                1906       gTr(iDOFe)=gTr(iDOFe) - R_FeC_CDOM*graz2CDOM
                1907 # else
8fbfd1f382 Oliv*1908       gTr(iDOC )=gTr(iDOC )             - R_CP_CDOM*graz2CDOM
                1909       gTr(iDON )=gTr(iDON )             - R_NP_CDOM*graz2CDOM
                1910       gTr(iDOP )=gTr(iDOP ) - graz2CDOM
                1911       gTr(iDOFe)=gTr(iDOFe)             - R_FeP_CDOM*graz2CDOM
1c72adde5c Oliv*1912 # endif
8fbfd1f382 Oliv*1913 #endif
                1914
                1915       DO jp = 1, nplank
                1916       IF (isPrey(jp).NE.0) THEN
                1917        gTr(ic+jp-1)= gTr(ic+jp-1) - preygraz(jp)
                1918 #ifdef DARWIN_ALLOW_NQUOTA
                1919        gTr(in+jp-1)=gTr(in+jp-1) - preygraz(jp)*Qn(jp)
                1920 #endif
                1921 #ifdef DARWIN_ALLOW_PQUOTA
                1922        gTr(ip+jp-1)=gTr(ip+jp-1) - preygraz(jp)*Qp(jp)
                1923 #endif
                1924 #ifdef DARWIN_ALLOW_SIQUOTA
                1925        gTr(isi+jp-1)=gTr(isi+jp-1) - preygraz(jp)*Qsi(jp)
                1926 #endif
                1927 #ifdef DARWIN_ALLOW_FEQUOTA
                1928        gTr(ife+jp-1)=gTr(ife+jp-1) - preygraz(jp)*Qfe(jp)
                1929 #endif
                1930       ENDIF
                1931       ENDDO
                1932 #ifdef DARWIN_ALLOW_CHLQUOTA
                1933       DO jp = 1, nPhoto
                1934       IF (isPrey(jp).NE.0) THEN
                1935        gTr(iChl+jp-1)=gTr(iChl+jp-1) - preygraz(jp)*QChl(jp)
                1936       ENDIF
                1937       ENDDO
                1938 #endif
c7b6c66d45 Oliv*1939 #ifdef DARWIN_ALLOW_CSTORE
                1940       DO jp = 1, nPhoto
                1941       IF (isPrey(jp).NE.0) THEN
                1942        gTr(ich+jp-1)=gTr(ich+jp-1) - preygraz(jp)*Qch(jp)
                1943       ENDIF
                1944       ENDDO
                1945 #endif
8fbfd1f382 Oliv*1946
                1947       DO jz = 1, nplank
                1948       IF (isPred(jz).NE.0) THEN
                1949        gTr(ic+jz-1)=gTr(ic+jz-1) + predgrazc(jz)
                1950 #ifdef DARWIN_ALLOW_NQUOTA
                1951        gTr(in+jz-1)=gTr(in+jz-1) + predgrazn(jz)
                1952 #endif
                1953 #ifdef DARWIN_ALLOW_PQUOTA
                1954        gTr(ip+jz-1)=gTr(ip+jz-1) + predgrazp(jz)
                1955 #endif
                1956 #ifdef DARWIN_ALLOW_FEQUOTA
                1957        gTr(ife+jz-1)=gTr(ife+jz-1) + predgrazfe(jz)
                1958 #endif
                1959       ENDIF
                1960       ENDDO
                1961
ee4178dfbe Oliv*1962 #ifdef DARWIN_DIAG_PERTYPE
                1963       DO jp = 1, nplank
8fbfd1f382 Oliv*1964         diags(iGRplank+jp-1) = preygraz(jp)
                1965       ENDDO
ee4178dfbe Oliv*1966       DO jz = 1, nplank
73aabcca42 Oliv*1967         diags(iGrGn+jz-1) = predgrazc(jz)
a93256cee4 Oliv*1968         diags(iGrGC+jz-1) = grazrate(jz)
73aabcca42 Oliv*1969       ENDDO
ee4178dfbe Oliv*1970 #endif
73aabcca42 Oliv*1971
c7b6c66d45 Oliv*1972 C=======================================================================
                1973 C==== Mortality and parameterized exudation ============================
8fbfd1f382 Oliv*1974       exude_DOC  = 0.0 _d 0
                1975       exude_POC  = 0.0 _d 0
                1976       exude_DON  = 0.0 _d 0
                1977       exude_PON  = 0.0 _d 0
                1978       exude_DOFe = 0.0 _d 0
                1979       exude_POFe = 0.0 _d 0
                1980       exude_DOP  = 0.0 _d 0
                1981       exude_POP  = 0.0 _d 0
                1982       exude_POSi = 0.0 _d 0
                1983       exude_PIC  = 0.0 _d 0
                1984       respir     = 0.0 _d 0
8752612d6f Oliv*1985       respirN    = 0.0 _d 0
                1986       respirP    = 0.0 _d 0
                1987       respirFe   = 0.0 _d 0
                1988       respirSi   = 0.0 _d 0
5e7acb36b1 daat*1989       mortDVM_DOC = 0.0 _d 0
                1990       mortDVM_DON = 0.0 _d 0
                1991       mortDVM_DOP = 0.0 _d 0
                1992       mortDVM_DOFe = 0.0 _d 0
                1993       mortDVM_POC = 0.0 _d 0
                1994       mortDVM_PON = 0.0 _d 0
                1995       mortDVM_POP = 0.0 _d 0
                1996       mortDVM_POFe = 0.0 _d 0
                1997       mortDVM_POSi = 0.0 _d 0
8fbfd1f382 Oliv*1998
                1999       DO jp = 1, nplank
                2000         Xe = MAX(0 _d 0, X(jp) - Xmin(jp))
7e97ec2956 Oliv*2001         mortX = mort(jp)*Xe*mortTempFunc**tempMort(jp)
                2002         mortX2= mort2(jp)*Xe*Xe*mort2TempFunc**tempMort2(jp)
8fbfd1f382 Oliv*2003
                2004         mort_c(jp) = mortX + mortX2
                2005
c7b6c66d45 Oliv*2006 #ifdef DARWIN_ALLOW_CSTORE
                2007         exude_DOC=exude_DOC+(1.-ExportFracMort(jp)) *mortX*Qctot(jp)
                2008      &                     +(1.-ExportFracMort2(jp))*mortX2*Qctot(jp)
                2009         exude_POC=exude_POC+    ExportFracMort(jp)  *mortX*Qctot(jp)
                2010      &                     +    ExportFracMort2(jp) *mortX2*Qctot(jp)
                2011 #else
8fbfd1f382 Oliv*2012         exude_DOC = exude_DOC + (1.-ExportFracMort(jp)) *mortX
                2013      &                        + (1.-ExportFracMort2(jp))*mortX2
                2014         exude_POC = exude_POC +     ExportFracMort(jp)  *mortX
                2015      &                        +     ExportFracMort2(jp) *mortX2
c7b6c66d45 Oliv*2016 #endif
8fbfd1f382 Oliv*2017
                2018         exude_DON = exude_DON + (1.-ExportFracMort(jp)) *mortX *Qn(jp)
                2019      &                        + (1.-ExportFracMort2(jp))*mortX2*Qn(jp)
                2020         exude_PON = exude_PON +     ExportFracMort(jp)  *mortX *Qn(jp)
                2021      &                        +     ExportFracMort2(jp) *mortX2*Qn(jp)
                2022
                2023         exude_DOP = exude_DOP + (1.-ExportFracMort(jp)) *mortX *Qp(jp)
                2024      &                        + (1.-ExportFracMort2(jp))*mortX2*Qp(jp)
                2025         exude_POP = exude_POP +     ExportFracMort(jp)  *mortX *Qp(jp)
                2026      &                        +     ExportFracMort2(jp) *mortX2*Qp(jp)
                2027
                2028         exude_DOFe= exude_DOFe+ (1.-ExportFracMort(jp)) *mortX *Qfe(jp)
                2029      &                        + (1.-ExportFracMort2(jp))*mortX2*Qfe(jp)
                2030         exude_POFe= exude_POFe+     ExportFracMort(jp)  *mortX *Qfe(jp)
                2031      &                        +     ExportFracMort2(jp) *mortX2*Qfe(jp)
                2032
                2033         exude_POSi = exude_POSi + mort_c(jp)*Qsi(jp)
                2034
8752612d6f Oliv*2035         exude_PIC = exude_PIC + mort_c(jp)*R_PICPOC(jp)
8fbfd1f382 Oliv*2036
8752612d6f Oliv*2037 C Respiration of non-phototrophs
                2038 C Respired carbon goes to DIC.  Other non-quota elements go to inorganic
                2039 C matter pools according to plankton stoichiometry.  PIC also goes back
                2040 C to DIC as if it had never been synthesized.
                2041         IF (isPhoto(jp) .EQ. 0) THEN
                2042          respir_c = respRate(jp)*Xe*reminTempFunc
                2043          gTr(ic+jp-1) = gTr(ic+jp-1) - respir_c
fd20fa810c Oliv*2044 #ifdef DARWIN_DIAG_PERTYPE
                2045          diags(iResp+jp-1) = respir_c
                2046 #endif
8752612d6f Oliv*2047          respir = respir + respir_c*(1 _d 0 + R_PICPOC(jp))
                2048 #ifndef DARWIN_ALLOW_NQUOTA
                2049          respirN = respirN + respir_c*R_NC(jp)
                2050 #endif
                2051 #ifndef DARWIN_ALLOW_PQUOTA
                2052          respirP = respirP + respir_c*R_PC(jp)
                2053 #endif
                2054 #ifndef DARWIN_ALLOW_FEQUOTA
                2055          respirFe = respirFe + respir_c*R_FeC(jp)
                2056 #endif
                2057 #ifndef DARWIN_ALLOW_SIQUOTA
                2058          respirSi = respirSi + respir_c*R_SiC(jp)
                2059 #endif
                2060         ENDIF
65ac4e45b6 Oliv*2061
8752612d6f Oliv*2062         gTr(ic+jp-1)=gTr(ic+jp-1)  - mort_c(jp)
8fbfd1f382 Oliv*2063 #ifdef DARWIN_ALLOW_NQUOTA
                2064         gTr(in+jp-1)=gTr(in+jp-1)  - mort_c(jp)*Qn(jp)
                2065 #endif
                2066 #ifdef DARWIN_ALLOW_PQUOTA
                2067         gTr(ip+jp-1)=gTr(ip+jp-1)  - mort_c(jp)*Qp(jp)
                2068 #endif
                2069 #ifdef DARWIN_ALLOW_SIQUOTA
                2070         gTr(isi+jp-1)=gTr(isi+jp-1) - mort_c(jp)*Qsi(jp)
                2071 #endif
                2072 #ifdef DARWIN_ALLOW_FEQUOTA
                2073         gTr(ife+jp-1)=gTr(ife+jp-1) - mort_c(jp)*Qfe(jp)
                2074 #endif
fd20fa810c Oliv*2075 #ifdef DARWIN_DIAG_PERTYPE
                2076         diags(iMort+jp-1) = mort_c(jp)
                2077 #endif
                2078
8fbfd1f382 Oliv*2079
                2080 #ifdef DARWIN_ALLOW_EXUDE
                2081         exude_DOC =
                2082      &  exude_DOC + (1.-ExportFracExude(jp))*kexcc(jp)*Xe
                2083         exude_POC =
                2084      &  exude_POC +     ExportFracExude(jp) *kexcc(jp)*Xe
                2085         exude_DON =
                2086      &  exude_DON + (1.-ExportFracExude(jp))*kexcn(jp)*Xe*Qn(jp)
                2087         exude_PON =
                2088      &  exude_PON +     ExportFracExude(jp) *kexcn(jp)*Xe*Qn(jp)
                2089         exude_DOP =
                2090      &  exude_DOP + (1.-ExportFracExude(jp))*kexcp(jp)*Xe*Qp(jp)
                2091         exude_POP =
                2092      &  exude_POP +     ExportFracExude(jp) *kexcp(jp)*Xe*Qp(jp)
                2093         exude_DOFe =
                2094      &  exude_DOFe + (1.-ExportFracExude(jp))*kexcfe(jp)*Xe*Qfe(jp)
                2095         exude_POFe =
                2096      &  exude_POFe +     ExportFracExude(jp) *kexcfe(jp)*Xe*Qfe(jp)
                2097         exude_POSi =
                2098      &  exude_POSi + kexcsi(jp)*Xe*Qsi(jp)
65ac4e45b6 Oliv*2099         exude_PIC = exude_PIC + kexcc(jp)*Xe*R_PICPOC(jp)
8fbfd1f382 Oliv*2100         gTr(ic+jp-1)=gTr(ic+jp-1)   - kexcc(jp)*Xe
                2101 #ifdef DARWIN_ALLOW_NQUOTA
                2102         gTr(in+jp-1)=gTr(in+jp-1)   - kexcn(jp)*Xe*Qn(jp)
                2103 #endif
                2104 #ifdef DARWIN_ALLOW_PQUOTA
                2105         gTr(ip+jp-1)=gTr(ip+jp-1)   - kexcp(jp)*Xe*Qp(jp)
                2106 #endif
                2107 #ifdef DARWIN_ALLOW_SIQUOTA
                2108         gTr(isi+jp-1)=gTr(isi+jp-1) - kexcsi(jp)*Xe*Qsi(jp)
                2109 #endif
                2110 #ifdef DARWIN_ALLOW_FEQUOTA
                2111         gTr(ife+jp-1)=gTr(ife+jp-1) - kexcfe(jp)*Xe*Qfe(jp)
                2112 #endif
                2113 #endif
5e7acb36b1 daat*2114
                2115 C====== Light-dependent mortality from DVM ====
                2116
                2117 #ifdef DARWIN_ALLOW_DVM
                2118        IF (isPred(jp) .EQ. 1) THEN
                2119         mortDVM_c = mortmaxDVM(jp)*Xe/(ksatDVM(jp)+Xe)*
                2120      &      ( 1 _d 0 - fracPARmort(jp) +
                2121      &        fracPARmort(jp)
                2122      &        *(PARtot/
                2123      &          ((ksatDVM(jp)/(ksatDVM(jp)+Xe))*ksatPARDVM(jp) + PARtot)
                2124      &         )
                2125      &      )*Xe
                2126
                2127         gTr(ic+jp-1)=gTr(ic+jp-1) - mortDVM_c
                2128 # ifdef DARWIN_ALLOW_NQUOTA
                2129         gTr(in+jp-1)=gTr(in+jp-1) - mortDVM_c*Qn(jp)
                2130 # endif
                2131 # ifdef DARWIN_ALLOW_PQUOTA
                2132         gTr(ip+jp-1)=gTr(ip+jp-1) - mortDVM_c*Qp(jp)
                2133 # endif
                2134 # ifdef DARWIN_ALLOW_SIQUOTA
                2135         gTr(isi+jp-1)=gTr(isi+jp-1) - mortDVM_c*Qsi(jp)
                2136 # endif
                2137 # ifdef DARWIN_ALLOW_FEQUOTA
                2138         gTr(ife+jp-1)=gTr(ife+jp-1) - mortDVM_c*Qfe(jp)
                2139 # endif
                2140
                2141         mortDVM_DOC = mortDVM_DOC + (1.-ExportFracDVM(jp))*mortDVM_c
                2142         mortDVM_POC = mortDVM_POC + ExportFracDVM(jp)*mortDVM_c
                2143
                2144         mortDVM_DON =
                2145      &  mortDVM_DON + (1.-ExportFracDVM(jp))*mortDVM_c*Qn(jp)
                2146         mortDVM_PON =
                2147      &  mortDVM_PON + ExportFracDVM(jp)*mortDVM_c*Qn(jp)
                2148         mortDVM_DOP =
                2149      &  mortDVM_DOP + (1.-ExportFracDVM(jp))*mortDVM_c*Qp(jp)
                2150         mortDVM_POP =
                2151      &  mortDVM_POP + ExportFracDVM(jp)*mortDVM_c*Qp(jp)
                2152         mortDVM_DOFe =
                2153      &  mortDVM_DOFe + (1.-ExportFracDVM(jp))*mortDVM_c*Qfe(jp)
                2154         mortDVM_POFe =
                2155      &  mortDVM_POFe + ExportFracDVM(jp)*mortDVM_c*Qfe(jp)
                2156         mortDVM_POSi =
                2157      &  mortDVM_POSi + mortDVM_c*Qsi(jp)
                2158        ENDIF
                2159 #endif /* DARWIN_ALLOW_DVM */
                2160
8fbfd1f382 Oliv*2161       ENDDO
                2162
                2163 #ifdef DARWIN_ALLOW_CHLQUOTA
                2164       DO jp = 1, nPhoto
                2165         gTr(iChl+jp-1)=gTr(iChl+jp-1) - mort_c(jp)*QChl(jp)
                2166       ENDDO
                2167 #endif
c7b6c66d45 Oliv*2168 #ifdef DARWIN_ALLOW_CSTORE
                2169       DO jp = 1, nPhoto
                2170         gTr(ich+jp-1)=gTr(ich+jp-1) - mort_c(jp)*Qch(jp)
                2171       ENDDO
                2172 #endif
8fbfd1f382 Oliv*2173
8752612d6f Oliv*2174       gTr(iDIC) = gTr(iDIC) + respir
                2175 #ifndef DARWIN_ALLOW_NQUOTA
                2176       gTr(iNO3) = gTr(iNO3) + respirN
                2177 #endif
                2178 #ifndef DARWIN_ALLOW_PQUOTA
                2179       gTr(iPO4) = gTr(iPO4) + respirP
                2180 #endif
                2181 #ifndef DARWIN_ALLOW_FEQUOTA
                2182       gTr(iFeT) = gTr(iFeT) + respirFe
                2183 #endif
                2184 #ifndef DARWIN_ALLOW_SIQUOTA
                2185       gTr(iSiO2) = gTr(iSiO2) + respirSi
                2186 #endif
8fbfd1f382 Oliv*2187
b61dca872c Oliv*2188       diags(iRespirDIC) = respir
                2189
5e7acb36b1 daat*2190       gTr(iDOC )=gTr(iDOC ) + exude_DOC + mortDVM_DOC
                2191       gTr(iDON )=gTr(iDON ) + exude_DON + mortDVM_DON
                2192       gTr(iDOP )=gTr(iDOP ) + exude_DOP + mortDVM_DOP
                2193       gTr(iDOFe)=gTr(iDOFe) + exude_DOFe + mortDVM_DOFe
8fbfd1f382 Oliv*2194
                2195       gTr(iPIC )=gTr(iPIC ) + exude_PIC
5e7acb36b1 daat*2196       gTr(iPOC )=gTr(iPOC ) + exude_POC + mortDVM_POC
                2197       gTr(iPON )=gTr(iPON ) + exude_PON + mortDVM_PON
                2198       gTr(iPOP )=gTr(iPOP ) + exude_POP + mortDVM_POP
                2199       gTr(iPOSi)=gTr(iPOSi) + exude_POSi + mortDVM_POSi
                2200       gTr(iPOFe)=gTr(iPOFe) + exude_POFe + mortDVM_POFe
8fbfd1f382 Oliv*2201 #ifdef DARWIN_ALLOW_CDOM
1c72adde5c Oliv*2202 # ifdef DARWIN_CDOM_UNITS_CARBON
b81bfeb644 Oliv*2203       exude_CDOM = fracCDOM*exude_DOC
1c72adde5c Oliv*2204       gTr(iDOC )=gTr(iDOC ) -           exude_CDOM
                2205       gTr(iDON )=gTr(iDON ) - R_NC_CDOM*exude_CDOM
                2206       gTr(iDOP )=gTr(iDOP ) - R_PC_CDOM*exude_CDOM
                2207       gTr(iDOFe)=gTr(iDOFe) - R_FeC_CDOM*exude_CDOM
                2208 # else
b81bfeb644 Oliv*2209       exude_CDOM = fracCDOM*exude_DOP
8fbfd1f382 Oliv*2210       gTr(iDOC )=gTr(iDOC )              - R_CP_CDOM*exude_CDOM
                2211       gTr(iDON )=gTr(iDON )              - R_NP_CDOM*exude_CDOM
                2212       gTr(iDOP )=gTr(iDOP ) - exude_CDOM
                2213       gTr(iDOFe)=gTr(iDOFe)              - R_FeP_CDOM*exude_CDOM
1c72adde5c Oliv*2214 # endif
b81bfeb644 Oliv*2215       gTr(iCDOM)=gTr(iCDOM) + exude_CDOM
8fbfd1f382 Oliv*2216 #endif
                2217
                2218 #endif /* ALLOW_DARWIN */
                2219
                2220       RETURN
                2221       END SUBROUTINE