pkg/darwin/darwin_forcing.F

8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002 #ifdef ALLOW_SEAICE
                0003 #include "SEAICE_OPTIONS.h"
                0004 #endif
                0005
                0006 CBOP
                0007 C !ROUTINE: DARWIN_FORCING
                0008 C !INTERFACE: ==========================================================
                0009       SUBROUTINE DARWIN_FORCING( Ptrdummy,
                0010      I                        bi, bj, iMin, iMax, jMin, jMax,
                0011      I                        myIter, myTime, myThid )
                0012
                0013 C !DESCRIPTION:
                0014
                0015 C !USES: ===============================================================
                0016       IMPLICIT NONE
                0017 #include "SIZE.h"
                0018 #include "GRID.h"
                0019 #include "EEPARAMS.h"
a8b28cdf62 Oliv*0020 #ifdef ALLOW_EXCH2
                0021 #include "W2_EXCH2_SIZE.h"
                0022 #include "W2_EXCH2_TOPOLOGY.h"
                0023 #include "W2_EXCH2_PARAMS.h"
                0024 #endif
8fbfd1f382 Oliv*0025 #include "PARAMS.h"
                0026 #include "DYNVARS.h"
                0027 #include "PTRACERS_SIZE.h"
                0028 #include "PTRACERS_PARAMS.h"
                0029 #include "PTRACERS_FIELDS.h"
                0030 #ifdef ALLOW_SEAICE
                0031 #include "SEAICE_SIZE.h"
                0032 #include "SEAICE.h"
                0033 #endif
                0034 #ifdef ALLOW_RADTRANS
                0035 #include "RADTRANS_SIZE.h"
                0036 #endif
                0037 #ifdef ALLOW_DARWIN
                0038 #include "GCHEM.h"
                0039 #include "DARWIN_SIZE.h"
                0040 #include "DARWIN_INDICES.h"
                0041 #include "DARWIN_DIAGS.h"
                0042 #include "DARWIN_PARAMS.h"
                0043 #include "DARWIN_TRAITS.h"
                0044 #include "DARWIN_FIELDS.h"
                0045 #include "DARWIN_EXF_FIELDS.h"
                0046 #endif
e79020cc6e Oliv*0047 #ifdef DARWIN_NUTRIENT_RUNOFF
                0048 #include "DARWIN_EXF_PARAMS.h"
                0049 #endif
8fbfd1f382 Oliv*0050
                0051 C !INPUT PARAMETERS: ===================================================
                0052 C  Ptrdummy             :: dummy for darwin2 compatibility
                0053 C  myThid               :: thread number
                0054       _RL Ptrdummy(*)
                0055       _RL myTime
                0056       INTEGER iMin, iMax, jMin, jMax, bi, bj, myIter, myThid
                0057 CEOP
                0058
                0059 #ifdef ALLOW_DARWIN
                0060
                0061 C!LOCAL VARIABLES: ====================================================
                0062 C  i,j                  :: loop indices
                0063 C  k                    :: vertical level
                0064       LOGICAL  DIAGNOSTICS_IS_ON
                0065       EXTERNAL DIAGNOSTICS_IS_ON
                0066       INTEGER i,j,k,kdn,iTr,l,isub
                0067       CHARACTER*8 diagname
                0068       _RL dTsub(Nr)
                0069       _RL midTime, subTime
63ffa46034 Oliv*0070 #if defined ALLOW_DIAGNOSTICS && ! defined DARWIN_ALLOW_CONS
                0071       _RL DARWIN_minFeLoss(sNx,sNy,Nr,nSx,nSy)
                0072 #endif
8fbfd1f382 Oliv*0073       _RL gPtr(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nDarwin)
                0074       _RL PAR(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0075       _RL diags(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr, darwin_nDiag)
                0076       _RL gDIC(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
                0077       _RL gALK(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
714cf988dd Oliv*0078 #ifdef DARWIN_DIAG_TENDENCIES
                0079       _RL gDICsurfForc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL gNO3surfForc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL gNO2surfForc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL gNH4surfForc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL gPO4surfForc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL gFeTsurfForc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL gSiO2surfForc(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
                0086 #ifdef DARWIN_ALLOW_CARBON
                0087       _RL gALKsurfForc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL gO2surfForc(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
                0089 #endif
                0090 #endif
8fbfd1f382 Oliv*0091       _RL gO2(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
3997029bbb Oliv*0092       _RL surfFe(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
8fbfd1f382 Oliv*0093       _RL freeFe(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr)
c8be5f4655 Oliv*0094       _RL sedFlxFe(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
c2ba01bd90 Oliv*0095       _RL tmp3d(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr)
82fd1e0266 Oliv*0096       _RL scavRate(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr)
fd32b5ab01 Oliv*0097       _RL scavLoss(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr)
66ade11a9e Oliv*0098       _RL scv,scav_pom
ad5342a8d6 Oliv*0099       _RL flx, POCl
6aa5674e10 Oliv*0100 #ifdef DARWIN_ALLOW_CARBON
                0101 C     Sediment fluxes
8a10163480 Oliv*0102       _RL OmegaC0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0103 #ifdef DARWIN_ALLOW_RADIv1
6aa5674e10 Oliv*0104       _RL DICFsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0105       _RL ALKFsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106       _RL O2Fsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL POCFsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0108       _RL CALFsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0109 C     Test variable for a diffusive temperature correction factor
                0110 C     (for now only implemented in the O2 flux)
                0111       _RL TcorrO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112       _RL TcorrDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL TcorrALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114 C     Dependent variables for sediment flux for validation
                0115       _RL fluxPIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116       _RL fluxPOC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117       _RL absT
                0118 #endif
8a10163480 Oliv*0119 #ifdef DARWIN_ALLOW_RADIv2
                0120       _RL DICFsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0121       _RL ALKFsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0122       _RL O2Fsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0123       _RL NO3Fsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0124       _RL PO4Fsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL NH4Fsediment(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126 C     Dependent variables for sediment flux for validation
                0127       _RL fluxPIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0128       _RL fluxPOC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0129       _RL fluxPOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0130       _RL fluxPON(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0131       _RL absT
                0132       _RL flxplnkC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0133       _RL flxplnkP(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0134       _RL flxplnkN(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0135       _RL sed_a1, sed_b1, sed_c1, sed_d1, sed_e1
                0136       _RL sed_a2, sed_b2, sed_c2, sed_d2, sed_e2
                0137       _RL sed_a3, sed_b3, sed_c3, sed_d3, sed_e3
5938ff335d Oliv*0138       _RL sed_a4, sed_b4, sed_c4, sed_d4, sed_e4, sed_f4, sed_g4
8a10163480 Oliv*0139       _RL sed_a5, sed_b5, sed_c5, sed_d5, sed_e5
                0140       _RL sed_a6, sed_b6, sed_c6, sed_d6, sed_e6
                0141 #endif
6aa5674e10 Oliv*0142 #endif
8fbfd1f382 Oliv*0143       _RL ptr(nDarwin), gtr(nDarwin), PARl(nlam)
                0144       _RL chlout(nPhoto)
                0145       _RL diagsl(darwin_nDiag)
5e7acb36b1 daat*0146       _RL sumpreyl(nplank)
                0147       _RL sumprey(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nplank)
                0148       _RL bioswimDVMup(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nplank)
                0149       _RL bioswimDVMdn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nplank)
                0150       _RL SwDr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nplank)
                0151       INTEGER max_loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nplank)
                0152       _RL sumPAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
8fbfd1f382 Oliv*0153       _RL photoTempFunc(nplank)
3f05298721 Oliv*0154       _RL hetTempFunc(nplank)
8fbfd1f382 Oliv*0155       _RL grazTempFunc(nplank)
                0156       _RL reminTempFunc
                0157       _RL mortTempFunc
                0158       _RL mort2TempFunc
                0159       _RL uptakeTempFunc
a092808e6b shlo*0160       _RL macromolTempFunc
ba0b6d5d33 Oliv*0161       _RL omegaCl
63ffa46034 Oliv*0162       _RL tmp, tmpFac
b61dca872c Oliv*0163       _RL sedFe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
2c93eb88ef Oliv*0164 #ifdef DARWIN_ALLOW_HYDROTHERMAL_VENTS
                0165       _RL ventFe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0166 #endif
6559acf5ff Oliv*0167       INTEGER k_debug
a8b28cdf62 Oliv*0168 #ifdef DARWIN_DEBUG
                0169       INTEGER iG, jG
                0170       INTEGER iBase, jBase
                0171       INTEGER iGjLoc, jGjLoc
                0172 #ifdef ALLOW_EXCH2
                0173       INTEGER tN
                0174 #endif /* ALLOW_EXCH2 */
                0175 #endif /* DARWIN_DEBUG */
8fbfd1f382 Oliv*0176
63ffa46034 Oliv*0177       DO k=1,Nr
                0178        DO j=1,sNy
0f72dca88c Oliv*0179         DO i=1,sNx
63ffa46034 Oliv*0180          DARWIN_minFeLoss(i,j,k,bi,bj) = 0 _d 0
                0181         ENDDO
                0182        ENDDO
                0183       ENDDO
                0184
bff492d006 Oliv*0185       DO l=1,darwin_nDiag
                0186        DO k=1,Nr
                0187         DO j=1,sNy
                0188          DO i=1,sNx
                0189           diags(i,j,k,l) = 0 _d 0
                0190          ENDDO
                0191         ENDDO
                0192        ENDDO
                0193       ENDDO
                0194
8fbfd1f382 Oliv*0195       DO k=1,Nr
                0196        dTsub(k) = PTRACERS_dTLev(k)/nsubtime
                0197       ENDDO
                0198
                0199 C     time at middle of sub-timestep
                0200       midTime = myTime - deltaTclock + .5*deltaTclock/nsubtime
                0201 C     time at end of sub-timestep
                0202       subTime = myTime - deltaTclock + deltaTclock/nsubtime
                0203       DO isub=1,nsubtime
                0204
                0205 C === reset tendencies =================================================
                0206       DO itr=1,nDarwin
                0207       DO k=1,Nr
                0208       DO j=jMin,jMax
                0209       DO i=iMin,iMax
                0210         gPtr(i,j,k,iTr) = 0.0 _d 0
                0211       ENDDO
                0212       ENDDO
                0213       ENDDO
                0214       ENDDO
                0215
                0216 C === light ============================================================
                0217 C     Initialize Chl from balanced-growth Chl:C if requested
                0218 C     and check Chl:C bounds.
                0219 C     Note: myIter has already been incremented
                0220 c      IF (myIter-1 .EQ. darwin_chlIter0) THEN
                0221 c        CALL DARWIN_INIT_CHL(bi, bj, subTime, myIter, myThid)
                0222 c      ENDIF
                0223
                0224 #ifdef ALLOW_SEAICE
                0225       IF (DARWIN_useSEAICE) THEN
                0226        DO j=jMin,jMax
                0227         DO i=iMin,iMax
                0228          iceFrac(i,j,bi,bj) = AREA(i,j,bi,bj)
                0229         ENDDO
                0230        ENDDO
                0231       ENDIF
                0232 #endif
                0233
                0234       CALL TIMER_START('DARWIN_LIGHT [DARWIN_FORCING]',myThid)
                0235 #ifdef ALLOW_RADTRANS
                0236       CALL DARWIN_LIGHT_RADTRANS(PAR,subTime,bi,bj,iMin,iMax,jMin,jMax,
                0237      &                        subTime,myIter,myThid)
                0238 #else
248275f1c4 Oliv*0239       CALL DARWIN_LIGHT(PAR, midTime, bi, bj, iMin, iMax, jMin, jMax,
8fbfd1f382 Oliv*0240      &               subTime, myIter, myThid)
                0241 #endif
                0242       CALL TIMER_STOP ('DARWIN_LIGHT [DARWIN_FORCING]',myThid)
                0243
                0244 C === dic ==============================================================
                0245 #ifdef DARWIN_ALLOW_CARBON
                0246 C carbon air-sea interaction
                0247       CALL TIMER_START('DARWIN_SURFFORCING [DARWIN_FORCING]',myThid)
                0248       CALL DARWIN_SURFFORCING(
                0249      O                    gDIC, gALK, gO2,
5e411acc9e Oliv*0250      I                    dTsub,bi,bj,iMin,iMax,jMin,jMax,
8fbfd1f382 Oliv*0251      I                    subTime,myIter,myThid)
                0252       CALL TIMER_STOP ('DARWIN_SURFFORCING [DARWIN_FORCING]',myThid)
                0253       DO j=jMin,jMax
                0254       DO i=iMin,iMax
e79020cc6e Oliv*0255         gPtr(i,j,1,iDIC) = gPtr(i,j,1,iDIC) + gDIC(i,j)
                0256         gPtr(i,j,1,iALK) = gPtr(i,j,1,iALK) + gALK(i,j)
                0257         gPtr(i,j,1,iO2)  = gPtr(i,j,1,iO2)  + gO2(i,j)
8fbfd1f382 Oliv*0258       ENDDO
                0259       ENDDO
c2ba01bd90 Oliv*0260
                0261 #ifdef NONLIN_FRSURF
                0262 C-    Account for change in level thickness
                0263 C     Other contributions to gPtr must come after these calls!
                0264       IF (nonlinFreeSurf.GT.0) THEN
                0265         CALL FREESURF_RESCALE_G(
                0266      I                           bi, bj, 1,
                0267      U                           gPtr(1-OLx,1-OLy,1,iDIC),
                0268      I                           myThid )
                0269         CALL FREESURF_RESCALE_G(
                0270      I                           bi, bj, 1,
                0271      U                           gPtr(1-OLx,1-OLy,1,iALK),
                0272      I                           myThid )
                0273         CALL FREESURF_RESCALE_G(
                0274      I                           bi, bj, 1,
                0275      U                           gPtr(1-OLx,1-OLy,1,iO2),
                0276      I                           myThid )
                0277       ENDIF
                0278 #endif /* NONLIN_FRSURF */
                0279 #endif /* DARWIN_ALLOW_CARBON */
8fbfd1f382 Oliv*0280
714cf988dd Oliv*0281 #ifdef DARWIN_DIAG_TENDENCIES
e79020cc6e Oliv*0282 C surface fluxes from E/P/runoff with constant/local concentrations
714cf988dd Oliv*0283       DO j=jMin,jMax
                0284       DO i=iMin,iMax
                0285         gDICsurfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iDIC)
                0286      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0287         gNO3surfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iNO3)
                0288      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0289         gNO2surfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iNO2)
                0290      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0291         gNH4surfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iNH4)
                0292      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0293         gPO4surfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iPO4)
                0294      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0295         gFeTsurfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iFeT)
                0296      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0297         gSiO2surfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iSiO2)
                0298      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0299 #ifdef DARWIN_ALLOW_CARBON
                0300         gALKsurfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iALK)
                0301      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0302         gO2surfForc(i,j) = surfaceForcingPTr(i,j,bi,bj,iO2)
                0303      &   * recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0304 #endif
                0305       ENDDO
                0306       ENDDO
                0307 #endif
                0308
e79020cc6e Oliv*0309 C === nutrient runoff ==================================================
                0310 #ifdef DARWIN_NUTRIENT_RUNOFF
16ebea87f8 Oliv*0311       CALL DARWIN_ADD_SURFFORC(
                0312      U                          gPtr(1-OLx,1-OLy,1,iDOC),
                0313      I                          DOCrunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0314      I                          bi, bj, iMin, iMax, jMin, jMax,
                0315      I                          myIter, myTime, myThid )
                0316
                0317       CALL DARWIN_ADD_SURFFORC(
                0318      U                          gPtr(1-OLx,1-OLy,1,iDON),
                0319      I                          DONrunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0320      I                          bi, bj, iMin, iMax, jMin, jMax,
                0321      I                          myIter, myTime, myThid )
                0322
                0323       CALL DARWIN_ADD_SURFFORC(
                0324      U                          gPtr(1-OLx,1-OLy,1,iDOP),
                0325      I                          DOPrunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0326      I                          bi, bj, iMin, iMax, jMin, jMax,
                0327      I                          myIter, myTime, myThid )
                0328
                0329       CALL DARWIN_ADD_SURFFORC(
                0330      U                          gPtr(1-OLx,1-OLy,1,iNO3),
                0331      I                          DINrunoff(1-OLx,1-OLy,bi,bj),
                0332      I                          R_NO3_DIN_runoff,
                0333      I                          bi, bj, iMin, iMax, jMin, jMax,
                0334      I                          myIter, myTime, myThid )
                0335
                0336       CALL DARWIN_ADD_SURFFORC(
                0337      U                          gPtr(1-OLx,1-OLy,1,iNO2),
                0338      I                          DINrunoff(1-OLx,1-OLy,bi,bj),
                0339      I                          R_NO2_DIN_runoff,
                0340      I                          bi, bj, iMin, iMax, jMin, jMax,
                0341      I                          myIter, myTime, myThid )
                0342
                0343       CALL DARWIN_ADD_SURFFORC(
                0344      U                          gPtr(1-OLx,1-OLy,1,iNH4),
                0345      I                          DINrunoff(1-OLx,1-OLy,bi,bj),
                0346      I                          R_NH4_DIN_runoff,
                0347      I                          bi, bj, iMin, iMax, jMin, jMax,
                0348      I                          myIter, myTime, myThid )
11d365acff Oliv*0349
                0350       CALL DARWIN_ADD_SURFFORC(
                0351      U                          gPtr(1-OLx,1-OLy,1,iNO3),
                0352      I                          NO3runoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0353      I                          bi, bj, iMin, iMax, jMin, jMax,
                0354      I                          myIter, myTime, myThid )
                0355
                0356       CALL DARWIN_ADD_SURFFORC(
                0357      U                          gPtr(1-OLx,1-OLy,1,iNO2),
                0358      I                          NO2runoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0359      I                          bi, bj, iMin, iMax, jMin, jMax,
                0360      I                          myIter, myTime, myThid )
                0361
                0362       CALL DARWIN_ADD_SURFFORC(
                0363      U                          gPtr(1-OLx,1-OLy,1,iNH4),
                0364      I                          NH4runoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0365      I                          bi, bj, iMin, iMax, jMin, jMax,
                0366      I                          myIter, myTime, myThid )
16ebea87f8 Oliv*0367
                0368       CALL DARWIN_ADD_SURFFORC(
                0369      U                          gPtr(1-OLx,1-OLy,1,iPO4),
                0370      I                          IPrunoff(1-OLx,1-OLy,bi,bj),
                0371      I                          R_DIP_IP_runoff,
                0372      I                          bi, bj, iMin, iMax, jMin, jMax,
                0373      I                          myIter, myTime, myThid )
                0374
                0375       CALL DARWIN_ADD_SURFFORC(
                0376      U                          gPtr(1-OLx,1-OLy,1,iSiO2),
                0377      I                          DSirunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0378      I                          bi, bj, iMin, iMax, jMin, jMax,
                0379      I                          myIter, myTime, myThid )
                0380
                0381       CALL DARWIN_ADD_SURFFORC(
                0382      U                          gPtr(1-OLx,1-OLy,1,iPOC),
                0383      I                          POCrunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0384      I                          bi, bj, iMin, iMax, jMin, jMax,
                0385      I                          myIter, myTime, myThid )
                0386
                0387       CALL DARWIN_ADD_SURFFORC(
                0388      U                          gPtr(1-OLx,1-OLy,1,iPON),
                0389      I                          PONrunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0390      I                          bi, bj, iMin, iMax, jMin, jMax,
                0391      I                          myIter, myTime, myThid )
                0392
                0393       CALL DARWIN_ADD_SURFFORC(
                0394      U                          gPtr(1-OLx,1-OLy,1,iPOP),
                0395      I                          POPrunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0396      I                          bi, bj, iMin, iMax, jMin, jMax,
                0397      I                          myIter, myTime, myThid )
                0398
                0399       CALL DARWIN_ADD_SURFFORC(
                0400      U                          gPtr(1-OLx,1-OLy,1,iDIC),
                0401      I                          DICrunoff(1-OLx,1-OLy,bi,bj), 1 _d 0,
                0402      I                          bi, bj, iMin, iMax, jMin, jMax,
                0403      I                          myIter, myTime, myThid )
                0404
                0405 #ifdef DARWIN_ALLOW_CARBON
                0406       CALL DARWIN_ADD_SURFFORC(
                0407      U                          gPtr(1-OLx,1-OLy,1,iALK),
                0408      I                          DICrunoff(1-OLx,1-OLy,bi,bj),
                0409      I                          R_ALK_DIC_runoff,
                0410      I                          bi, bj, iMin, iMax, jMin, jMax,
                0411      I                          myIter, myTime, myThid )
                0412 #endif
                0413
                0414       CALL DARWIN_ADD_SURFFORC(
                0415      U                          gPtr(1-OLx,1-OLy,1,iFeT),
                0416      I                          IPrunoff(1-OLx,1-OLy,bi,bj),
                0417      I                          R_DIP_IP_runoff*R_DFe_DIP_runoff,
                0418      I                          bi, bj, iMin, iMax, jMin, jMax,
                0419      I                          myIter, myTime, myThid )
                0420
                0421       CALL DARWIN_ADD_SURFFORC(
                0422      U                          gPtr(1-OLx,1-OLy,1,iDOFe),
                0423      I                          DOPrunoff(1-OLx,1-OLy,bi,bj),
                0424      I                          R_DOFe_DOP_runoff,
                0425      I                          bi, bj, iMin, iMax, jMin, jMax,
                0426      I                          myIter, myTime, myThid )
                0427
                0428       CALL DARWIN_ADD_SURFFORC(
                0429      U                          gPtr(1-OLx,1-OLy,1,iPOFe),
                0430      I                          POPrunoff(1-OLx,1-OLy,bi,bj),
                0431      I                          R_POFe_POP_runoff,
                0432      I                          bi, bj, iMin, iMax, jMin, jMax,
                0433      I                          myIter, myTime, myThid )
e79020cc6e Oliv*0434 #endif
                0435
8fbfd1f382 Oliv*0436 C === iron =============================================================
2c93eb88ef Oliv*0437
                0438 C     compute (and limit if DARWIN_MINFE) free iron
8fbfd1f382 Oliv*0439       CALL TIMER_START('DARWIN_FE_CHEM [DARWIN_FORCING]',myThid)
                0440       CALL DARWIN_FE_CHEM(
                0441      U                 Ptracer(1-OLx,1-OLy,1,bi,bj,iFeT),
63ffa46034 Oliv*0442      O                 freeFe(1-OLx,1-OLy,1),
                0443      U                 DARWIN_minFeLoss(1,1,1,bi,bj),
8fbfd1f382 Oliv*0444      I                 bi, bj, iMin, iMax, jMin, jMax, myThid)
                0445       CALL TIMER_STOP ('DARWIN_FE_CHEM [DARWIN_FORCING]',myThid)
                0446
2c93eb88ef Oliv*0447 C ----------------------------------------------------------------------
8fbfd1f382 Oliv*0448 C     scavenging
                0449       CALL TIMER_START('DARWIN_FE_SCAV [DARWIN_FORCING]',myThid)
                0450       DO k=1,Nr
                0451       DO j=jMin,jMax
                0452       DO i=iMin,iMax
66ade11a9e Oliv*0453 C       compute POM in mg/L
8fbfd1f382 Oliv*0454 #ifdef DARWIN_PART_SCAV_POP
66ade11a9e Oliv*0455         scav_pom = MAX(0 _d 0, Ptracer(i,j,k,bi,bj,iPOP))/scav_R_POPPOC
                0456         scv = scav_rat*scav_inter*(scav_pom**scav_exp)
8fbfd1f382 Oliv*0457 #elif defined(DARWIN_PART_SCAV)
66ade11a9e Oliv*0458         scav_pom = scav_POC_wgt*MAX(0 _d 0, Ptracer(i,j,k,bi,bj,iPOC))
28ac947a3e Oliv*0459      &           + scav_PSi_wgt*MAX(0 _d 0, Ptracer(i,j,k,bi,bj,iPOSi))
66ade11a9e Oliv*0460      &           + scav_PIC_wgt*MAX(0 _d 0, Ptracer(i,j,k,bi,bj,iPIC))
28ac947a3e Oliv*0461      &           + scav_degrPOM
66ade11a9e Oliv*0462         scv = scav_tau*scav_inter*(scav_pom**scav_exp)
8fbfd1f382 Oliv*0463 #else
                0464         scv = scav
                0465 #endif
82fd1e0266 Oliv*0466         scavRate(i,j,k) = scv
fd32b5ab01 Oliv*0467         scavLoss(i,j,k) = scv*MAX(0 _d 0, freefe(i,j,k))
                0468         gPtr(i,j,k,iFeT) = gPtr(i,j,k,iFeT) - scavLoss(i,j,k)
a59ebcf0a3 Oliv*0469 #ifdef DARWIN_ALLOW_CONS
5e411acc9e Oliv*0470         DARWIN_partScav(i,j,k,bi,bj) = DARWIN_partScav(i,j,k,bi,bj)
fd32b5ab01 Oliv*0471      &                               + dTsub(k)*scavLoss(i,j,k)
a59ebcf0a3 Oliv*0472 #endif
8fbfd1f382 Oliv*0473       ENDDO
                0474       ENDDO
                0475       ENDDO
                0476       CALL TIMER_STOP ('DARWIN_FE_SCAV [DARWIN_FORCING]',myThid)
                0477
2c93eb88ef Oliv*0478 C ----------------------------------------------------------------------
c59946c66d Oliv*0479 C     iron dust input
                0480       CALL TIMER_START('DARWIN_FE_DUST [DARWIN_FORCING]',myThid)
                0481       DO j=jMin,jMax
c2ba01bd90 Oliv*0482        DO i=iMin,iMax
                0483         tmp3d(i,j,1) = alpfe * inputFe(i,j,bi,bj)
                0484      &               * recip_drF(1) * recip_hFacC(i,j,1,bi,bj)
                0485        ENDDO
c59946c66d Oliv*0486       ENDDO
c2ba01bd90 Oliv*0487 #ifdef NONLIN_FRSURF
                0488 C-    Account for change in level thickness
                0489       IF (nonlinFreeSurf.GT.0) THEN
                0490         CALL FREESURF_RESCALE_G(
                0491      I                           bi, bj, 1,
                0492      U                           tmp3d,
                0493      I                           myThid )
                0494       ENDIF
                0495 #endif /* NONLIN_FRSURF */
                0496       DO j=jMin,jMax
                0497        DO i=iMin,iMax
                0498         gPtr(i,j,1,iFeT) = gPtr(i,j,1,iFeT) + tmp3d(i,j,1)
                0499        ENDDO
c59946c66d Oliv*0500       ENDDO
                0501       CALL TIMER_STOP ('DARWIN_FE_DUST [DARWIN_FORCING]',myThid)
                0502
2c93eb88ef Oliv*0503 C ----------------------------------------------------------------------
8fbfd1f382 Oliv*0504 C     iron sediment source (in bottom grid cell above kMaxFeSed)
                0505       CALL TIMER_START('DARWIN_FE_SED [DARWIN_FORCING]',myThid)
2c93eb88ef Oliv*0506       DO k=1,Nr
8fbfd1f382 Oliv*0507        DO j=jMin,jMax
                0508         DO i=iMin,iMax
2c93eb88ef Oliv*0509          sedFe(i,j,k) = 0 _d 0
                0510         ENDDO
                0511        ENDDO
                0512       ENDDO
                0513       DO j=jMin,jMax
                0514        DO i=iMin,iMax
                0515         k = kLowC(i,j,bi,bj)
                0516         IF (k .GE. kMinFeSed .AND. k .LE. kMaxFeSed) THEN
8fbfd1f382 Oliv*0517 #ifdef DARWIN_IRON_SED_SOURCE_VARIABLE
ad5342a8d6 Oliv*0518 # ifdef DARWIN_IRON_SED_SOURCE_POP
a092808e6b shlo*0519 C         mask is needed to make sure k-1 is not an ice shelf
8fbfd1f382 Oliv*0520           flx = fesedflux_pcm*wp_sink*R_CP_fesed*
                0521      &            MAX(0 _d 0, Ptracer(i,j,k-1,bi,bj,iPOP))
8e211684e0 Oliv*0522      &            *maskC(i,j,k-1,bi,bj)
ad5342a8d6 Oliv*0523 # else
                0524           POCl = MAX(0 _d 0, Ptracer(i,j,k,bi,bj,iPOC))
                0525           flx = fesedflux_pcm*wc_sink*POCl - fesedflux_min
                0526           flx = MAX(0. _d 0, flx)
                0527 # endif
8fbfd1f382 Oliv*0528 #else
                0529           flx = fesedflux
                0530 #endif
c8be5f4655 Oliv*0531           sedFlxFe(i,j) = flx
b61dca872c Oliv*0532           sedFe(i,j,k) = flx/(drF(k)*hFacC(i,j,k,bi,bj))
a59ebcf0a3 Oliv*0533 #ifdef DARWIN_ALLOW_CONS
5e411acc9e Oliv*0534           ironSedFlx(i,j,bi,bj) = ironSedFlx(i,j,bi,bj) + dTsub(k)*flx
a59ebcf0a3 Oliv*0535 #endif
c8be5f4655 Oliv*0536         ELSE
                0537          sedFlxFe(i,j) = 0 _d 0
2c93eb88ef Oliv*0538         ENDIF
8fbfd1f382 Oliv*0539        ENDDO
2c93eb88ef Oliv*0540       ENDDO
                0541       DO k=kMinFeSed, kMaxFeSed
c2ba01bd90 Oliv*0542 #ifdef NONLIN_FRSURF
                0543 C-    Account for change in level thickness
                0544        IF (nonlinFreeSurf.GT.0) THEN
                0545         CALL FREESURF_RESCALE_G(
                0546      I                           bi, bj, k,
                0547      U                           sedFe,
                0548      I                           myThid )
                0549        ENDIF
                0550 #endif /* NONLIN_FRSURF */
                0551        DO j=jMin,jMax
                0552         DO i=iMin,iMax
                0553           gPtr(i,j,k,iFeT) = gPtr(i,j,k,iFeT) + sedFe(i,j,k)
                0554         ENDDO
                0555        ENDDO
8fbfd1f382 Oliv*0556       ENDDO
                0557       CALL TIMER_STOP ('DARWIN_FE_SED [DARWIN_FORCING]',myThid)
                0558
2c93eb88ef Oliv*0559 C ----------------------------------------------------------------------
                0560 C     iron sediment source (in bottom grid cell above kMaxFeSed)
                0561 C     iron source due to hydrothermal input into bottom layer below
                0562 C     kMinFeVents (ventHe3 flux is mmol.m-2.s-1)
                0563 #ifdef DARWIN_ALLOW_HYDROTHERMAL_VENTS
                0564       IF (darwin_haveVentHe3) THEN
                0565        CALL TIMER_START('DARWIN_FE_VENT [DARWIN_FORCING]',myThid)
                0566        DO k=1,Nr
                0567         DO j=jMin,jMax
                0568          DO i=iMin,iMax
                0569           ventFe(i,j,k) = 0 _d 0
                0570          ENDDO
                0571         ENDDO
                0572        ENDDO
                0573        DO j=jMin,jMax
                0574         DO i=iMin,iMax
                0575          k = kLowC(i,j,bi,bj)
                0576          IF (k .GE. kMinFeVent) THEN
                0577           flx = solFeVent*R_FeHe3_vent*ventHe3(i,j,bi,bj)
                0578           ventFe(i,j,k) = flx*recip_drF(k)*recip_hFacC(i,j,k,bi,bj)
                0579 # ifdef DARWIN_ALLOW_CONS
                0580           ironVentFlx(i,j,bi,bj)=ironVentFlx(i,j,bi,bj) + dTsub(k)*flx
                0581 # endif
                0582          ENDIF
                0583         ENDDO
                0584        ENDDO
                0585        DO k=kMinFeVent,Nr
                0586 # ifdef NONLIN_FRSURF
                0587 C-    Account for change in level thickness
                0588         IF (nonlinFreeSurf.GT.0) THEN
                0589          CALL FREESURF_RESCALE_G(
                0590      I                            bi, bj, k,
                0591      U                            ventFe,
                0592      I                            myThid )
                0593         ENDIF
                0594 # endif /* NONLIN_FRSURF */
                0595         DO j=jMin,jMax
                0596          DO i=iMin,iMax
                0597            gPtr(i,j,k,iFeT) = gPtr(i,j,k,iFeT) + ventFe(i,j,k)
                0598          ENDDO
                0599         ENDDO
                0600        ENDDO
                0601        CALL TIMER_STOP ('DARWIN_FE_VENT [DARWIN_FORCING]',myThid)
                0602       ENDIF
                0603 #endif
                0604
8fbfd1f382 Oliv*0605 C === plankton =========================================================
                0606       CALL TIMER_START('DARWIN_PLANKTON [DARWIN_FORCING]',myThid)
a8b28cdf62 Oliv*0607
                0608 #ifdef DARWIN_DEBUG
                0609 C Set dimensions:
                0610 C-    default:
                0611       iGjLoc = 0
                0612       jGjLoc = 1
                0613       iBase = myXGlobalLo-1 + (bi-1)*sNx
                0614       jBase = myYGlobalLo-1 + (bj-1)*sNy
                0615 #ifdef ALLOW_EXCH2
                0616       IF ( W2_useE2ioLayOut ) THEN
                0617         tN = W2_myTileList(bi,bj)
                0618         iBase = exch2_txGlobalo(tN) - 1
                0619         jBase = exch2_tyGlobalo(tN) - 1
                0620         IF   ( exch2_mydNx(tN) .GT. exch2_global_Nx ) THEN
                0621 C-      face x-size larger than glob-size : fold it
                0622           iGjLoc = 0
                0623           jGjLoc = exch2_mydNx(tN) / exch2_global_Nx
                0624         ELSEIF ( exch2_tNy(tN) .GT. exch2_global_Ny ) THEN
                0625 C-      tile y-size larger than glob-size : make a long line
                0626           iGjLoc = exch2_mydNx(tN)
                0627           jGjLoc = 0
                0628         ELSE
                0629 C-      default (face fit into global-IO-array)
                0630           iGjLoc = 0
                0631           jGjLoc = 1
                0632         ENDIF
                0633       ENDIF
                0634 #endif /* ALLOW_EXCH2 */
                0635 #endif /* DARWIN_DEBUG */
                0636
8fbfd1f382 Oliv*0637       do k=1,Nr
                0638       do j=jMin,jMax
a8b28cdf62 Oliv*0639 #ifdef DARWIN_DEBUG
                0640        iG = iBase + (j-1)*iGjLoc
                0641        jG = jBase + (j-1)*jGjLoc + 1
                0642 #endif /* DARWIN_DEBUG */
                0643
                0644        do i=iMin,iMax
                0645        if (hFacC(i,j,k,bi,bj) .GT. 0.) then
6559acf5ff Oliv*0646         k_debug = 0
                0647 #ifdef DARWIN_DEBUG
                0648         IF (iG+i.eq.iDEBUG.and.jG.eq.jDEBUG) THEN
                0649          IF (kDEBUG.EQ.0 .OR. k.EQ.kDebug) THEN
                0650           k_debug = k
                0651          ENDIF
                0652         ENDIF
                0653 #endif
8fbfd1f382 Oliv*0654         CALL DARWIN_TEMPFUNC(Theta(i,j,k,bi,bj),
3f05298721 Oliv*0655      &         photoTempFunc, hetTempFunc, grazTempFunc,
                0656      &         reminTempFunc, mortTempFunc, mort2TempFunc,
a092808e6b shlo*0657      &         uptakeTempFunc, macromolTempFunc, myThid)
8fbfd1f382 Oliv*0658         DO iTr=1,nDarwin
                0659          ptr(iTr) = Ptracer(i, j, k, bi, bj, iTr)
                0660          gtr(iTr) = gPtr(i, j, k, iTr)
                0661         ENDDO
                0662         DO l=1,nlam
26960c5749 Oliv*0663          PARl(l) = MAX(0 _d 0, PAR(i, j, k, l))
8fbfd1f382 Oliv*0664         ENDDO
ba0b6d5d33 Oliv*0665 #ifdef DARWIN_ALLOW_CARBON
                0666         omegaCl = omegaC(i,j,k,bi,bj)
                0667 #else
                0668         omegaCl = 0 _d 0
                0669 #endif
8fbfd1f382 Oliv*0670         CALL DARWIN_PLANKTON(
                0671      I                 ptr,
                0672      U                 gtr,
5e7acb36b1 daat*0673      O                 chlout, diagsl, sumpreyl,
8fbfd1f382 Oliv*0674      I                 PARl, photoTempFunc,
3f05298721 Oliv*0675      I                 hetTempFunc,
8fbfd1f382 Oliv*0676      I                 grazTempFunc,
                0677      I                 reminTempFunc,
                0678      I                 mortTempFunc, mort2TempFunc,
a092808e6b shlo*0679      I                 uptakeTempFunc, macromolTempFunc,
ba0b6d5d33 Oliv*0680      I                 omegaCl,
6559acf5ff Oliv*0681      I                 k_debug,
8fbfd1f382 Oliv*0682      I                 subTime, myIter, myThid )
                0683         DO l=1,darwin_nDiag
                0684          diags(i, j, k, l) = diagsl(l)
                0685         ENDDO
5e7acb36b1 daat*0686
                0687         sumprey(i,j,k,:) = sumpreyl
                0688
8fbfd1f382 Oliv*0689 #ifdef DARWIN_ALLOW_CONS
5e411acc9e Oliv*0690         DARWIN_Nfix(i,j,k,bi,bj)   = DARWIN_Nfix(i,j,k,bi,bj)
                0691      &                             + dTsub(k)*diagsl(iNfix)
                0692         DARWIN_Ndenit(i,j,k,bi,bj) = DARWIN_Ndenit(i,j,k,bi,bj)
                0693      &                             + dTsub(k)*diagsl(iDenitN)
                0694         DARWIN_O2prod(i,j,k,bi,bj) = DARWIN_O2prod(i,j,k,bi,bj)
                0695      &                             + dTsub(k)*diagsl(iProdO2)
                0696         DARWIN_O2cons(i,j,k,bi,bj) = DARWIN_O2cons(i,j,k,bi,bj)
                0697      &                             + dTsub(k)*diagsl(iConsO2)
                0698         DARWIN_AlkSrc(i,j,k,bi,bj) = DARWIN_AlkSrc(i,j,k,bi,bj)
                0699      &                             + dTsub(k)*diagsl(iSrcAlk)
8fbfd1f382 Oliv*0700 #endif
                0701 #ifndef DARWIN_ALLOW_CHLQUOTA
                0702 #ifdef ALLOW_RADTRANS
                0703         DO l=1,nPhoto
                0704          chlPrev(i, j, k, bi, bj, l) = chlout(l)
                0705         ENDDO
                0706 #else
                0707         tmp = 0 _d 0
                0708         DO l=1,nPhoto
                0709          tmp = tmp + chlout(l)
                0710         ENDDO
                0711         chlPrev(i, j, k, bi, bj) = tmp
                0712 #endif
                0713 #endif
                0714         DO iTr=1,nDarwin
                0715          gPtr(i, j, k, iTr) = gtr(iTr)
                0716         ENDDO
a8b28cdf62 Oliv*0717        endif
                0718        enddo
8fbfd1f382 Oliv*0719       enddo
                0720       enddo
5e7acb36b1 daat*0721
8fbfd1f382 Oliv*0722       CALL TIMER_STOP('DARWIN_PLANKTON [DARWIN_FORCING]',myThid)
                0723
5e7acb36b1 daat*0724 C === DVM ==============================================================
                0725
                0726       CALL TIMER_START('DARWIN_DVM [DARWIN_FORCING]',myThid)
                0727
                0728 #ifdef DARWIN_ALLOW_DVM
                0729 C Find DVM swim direction to go to max food availablity
                0730       max_loc = MAXLOC(sumprey, DIM=3)
                0731       DO iTr=1,nplank
                0732        DO k=1,Nr
                0733         DO j=jMin,jMax
                0734          DO i=iMin,iMax
                0735           IF (k.LT.max_loc(i,j,iTr)) THEN
                0736            SwDr(i,j,k,iTr) = -1 _d 0
                0737           ELSEIF (k.GT.max_loc(i,j,iTr)) THEN
                0738            SwDr(i,j,k,iTr) = 1 _d 0
                0739           ELSE
                0740            SwDr(i,j,k,iTr) = 0 _d 0
                0741           ENDIF
                0742          ENDDO
                0743         ENDDO
                0744        ENDDO
                0745       ENDDO
                0746
                0747 C Find DVM swim direction based on light and max food availability
                0748 C Day
                0749 C If PAR at top >0 (or above some level in future),
                0750 C swim down (toward PARpref)
                0751 C If PAR at top >0 and PARtot < PARpref, swim up
                0752 C If PAR at top >0 and PARtot = PARpref, no swim
                0753
                0754 C Night--see logic in SwDr above
                0755 C If PAR at top 0 (or below some light level in future),
                0756 C and prey available at that depth < food max, swim toward food max
                0757 C If PAR at top 0 (or below some light level in future),
                0758 C and prey available at that depth = food max, don't swim; stay put
                0759
                0760       DO k=1,Nr
                0761        DO j=jMin,jMax
                0762         DO i=iMin,iMax
                0763          sumPAR(i,j,k) = 0 _d 0
                0764         ENDDO
                0765        ENDDO
                0766       ENDDO
                0767       DO l=1,nlam
                0768        DO k=1,Nr
                0769         DO j=jMin,jMax
                0770          DO i=iMin,iMax
                0771           sumPAR(i,j,k) = sumPAR(i,j,k) + PAR(i,j,k,l)
                0772          ENDDO
                0773         ENDDO
                0774        ENDDO
                0775       ENDDO
                0776
                0777       DO iTr=1,nplank
                0778        DO k=1,Nr
                0779         DO j=jMin,jMax
                0780          DO i=iMin,iMax
                0781           IF (sumPAR(i,j,1) .GT. 0) THEN
                0782 C          day
                0783            IF (sumPAR(i,j,k) .GT. PARpref(iTr)) THEN
                0784               bioswimDVMup(i,j,k,iTr) = 0 _d 0
                0785               bioswimDVMdn(i,j,k,iTr) = bioswimDVM(iTr)
                0786            ELSEIF (sumPAR(i,j,k). LT. PARpref(iTr)) THEN
                0787               bioswimDVMup(i,j,k,iTr) = bioswimDVM(iTr)
                0788               bioswimDVMdn(i,j,k,iTr) = 0 _d 0
                0789            ELSEIF (sumPAR(i,j,k) .EQ. PARpref(iTr)) THEN
                0790               bioswimDVMup(i,j,k,iTr) = 0 _d 0
                0791               bioswimDVMdn(i,j,k,iTr) = 0 _d 0
                0792            ENDIF
                0793           ELSE
                0794 C          night
                0795            IF (SwDr(i,j,k,iTr) .GT. 0.0) THEN
                0796                bioswimDVMup(i,j,k,iTr) = bioswimDVM(iTr)
                0797                bioswimDVMdn(i,j,k,iTr) = 0 _d 0
                0798            ELSEIF (SwDr(i,j,k,iTr) .EQ. 0.0) THEN
                0799                bioswimDVMup(i,j,k,iTr) = 0 _d 0
                0800                bioswimDVMdn(i,j,k,iTr) = 0 _d 0
                0801            ELSE
                0802                bioswimDVMup(i,j,k,iTr) = 0 _d 0
                0803                bioswimDVMdn(i,j,k,iTr) = bioswimDVM(iTr)
                0804            ENDIF
                0805           ENDIF
                0806          ENDDO
                0807         ENDDO
                0808        ENDDO
                0809       ENDDO
                0810 #else /* DARWIN_ALLOW_DVM */
                0811       DO iTr=1,nplank
                0812        DO k=1,Nr
                0813         DO j=jMin,jMax
                0814          DO i=iMin,iMax
                0815           bioswimDVMup(i,j,k,iTr) = 0 _d 0
                0816           bioswimDVMdn(i,j,k,iTr) = 0 _d 0
                0817          ENDDO
                0818         ENDDO
                0819        ENDDO
                0820       ENDDO
                0821 #endif /* DARWIN_ALLOW_DVM */
                0822
                0823       CALL TIMER_STOP('DARWIN_DVM [DARWIN_FORCING]',myThid)
                0824
5a97f77379 Oliv*0825 C === pre-sinking diags ================================================
                0826 #ifdef ALLOW_DIAGNOSTICS
                0827       IF (useDIAGNOSTICS) THEN
                0828        CALL TIMER_START('DIAGS_FILL [DARWIN_FORCING]',myThid)
                0829        DO iTr=1,nDarwin
                0830         diagname = 'gECO'//PTRACERS_ioLabel(iTr)
                0831         CALL DIAGNOSTICS_FILL(gPtr(1-OLx,1-OLy,1,iTr), diagname,
5e411acc9e Oliv*0832      &                        0,Nr,2,bi,bj,myThid)
5a97f77379 Oliv*0833        ENDDO
                0834        CALL TIMER_STOP ('DIAGS_FILL [DARWIN_FORCING]',myThid)
                0835       ENDIF
                0836 #endif
                0837
6aa5674e10 Oliv*0838 C === Sediment fluxes ==========================================================
                0839 #ifdef DARWIN_ALLOW_CARBON
8a10163480 Oliv*0840 C save omegaC above seafloor
                0841       DO j=jMin,jMax
                0842       DO i=iMin,iMax
                0843 C      grid cell just above sea floor
                0844        k = kLowC(i,j,bi,bj)
                0845        IF (k .GT. 0) THEN
                0846         OmegaC0(i,j) = omegaC(i,j,k,bi,bj)
                0847        ELSE
                0848         OmegaC0(i,j) = 0. _d 0
                0849        ENDIF
                0850       ENDDO
                0851       ENDDO
                0852
                0853 #ifdef DARWIN_ALLOW_RADIv1
6aa5674e10 Oliv*0854       DO j=jMin,jMax
                0855       DO i=iMin,iMax
                0856 C      grid cell just above sea floor
                0857        k = kLowC(i,j,bi,bj)
                0858        IF (k .GT. 0) THEN
                0859
                0860 C       Correction factor of diffusive flux with T, to potentially better
                0861 C       resolve sediment fluxes at coastal regions
                0862         absT = ABS(Theta(i,j,k,bi,bj))
                0863         TcorrO2(i,j) = (0.031558 _d 0+0.001428 _d 0*absT)/0.0335572 _d 0
                0864         TcorrDIC(i,j) = (0.015179 _d 0+0.000795 _d 0*absT)/0.016292 _d 0
                0865         TcorrALK(i,j) = (0.015179 _d 0+0.000795 _d 0*absT)/0.016292 _d 0
                0866
                0867 C       Dependent variables
086a45f245 Oliv*0868         fluxPIC(i,j) = wPIC_sink*MAX(0 _d 0,Ptracer(i,j,k,bi,bj,iPIC))
                0869         fluxPOC(i,j) = wC_sink*MAX(0 _d 0,Ptracer(i,j,k,bi,bj,iPOC))
6aa5674e10 Oliv*0870
                0871         DICFsediment(i,j) = (sed_a1*omegaC0(i,j)*sed_c +
                0872      &                       sed_b1*fluxPOC(i,j) +
                0873      &                       sed_c1*fluxPIC(i,j) +
                0874      &                       sed_d1*sed_c)*TcorrDIC(i,j)
                0875
                0876         ALKFsediment(i,j) = (sed_a2*omegaC0(i,j)*sed_c +
                0877      &                       sed_b2*fluxPOC(i,j) +
                0878      &                       sed_c2*fluxPIC(i,j) +
                0879      &                       sed_d2*sed_c)*TcorrALK(i,j)
                0880
                0881         O2Fsediment(i,j) = (sed_a3*omegaC0(i,j)*sed_c +
                0882      &                      sed_b3*fluxPOC(i,j) +
                0883      &                      sed_c3*fluxPIC(i,j) +
                0884      &                      sed_d3*sed_c)*TcorrO2(i,j)
                0885
                0886 C       Set the buried fluxes to 0 if positive
                0887         CALFsediment(i,j) = MIN(0 _d 0, -(sed_a5*omegaC0(i,j)*sed_c +
                0888      &    sed_b5*fluxPOC(i,j) + sed_c5*fluxPIC(i,j) + sed_d5*sed_c))
                0889
                0890         POCFsediment(i,j) = MIN(0 _d 0, -(sed_a4*omegaC0(i,j)*sed_c +
                0891      &    sed_b4*fluxPOC(i,j) + sed_c4*fluxPIC(i,j) + sed_d4*sed_c))
                0892
2a97638e3d Oliv*0893 #ifdef DARWIN_ALLOW_CONS
                0894 C       Accumulate for conservation check
                0895         radiFluxC(i,j,bi,bj) = radiFluxC(i,j,bi,bj) +
                0896      &       (DICFsediment(i,j) + POCFsediment(i,j) + CALFsediment(i,j))
                0897      &       *dTsub(k)*maskC(i,j,k,bi,bj)
                0898         radiFluxA(i,j,bi,bj) = radiFluxA(i,j,bi,bj) +
                0899      &        ALKFsediment(i,j)*dTsub(k)*maskC(i,j,k,bi,bj)
                0900         radiFluxO(i,j,bi,bj) = radiFluxO(i,j,bi,bj) +
                0901      &        O2Fsediment(i,j)*dTsub(k)*maskC(i,j,k,bi,bj)
                0902 #endif
                0903
6aa5674e10 Oliv*0904        ELSE
                0905
                0906         TcorrO2(i,j) = 0. _d 0
                0907         TcorrDIC(i,j) = 0. _d 0
                0908         TcorrALK(i,j) = 0. _d 0
                0909         DICFsediment(i,j) = 0. _d 0
                0910         ALKFsediment(i,j) = 0. _d 0
                0911         O2Fsediment(i,j) = 0. _d 0
                0912         CALFsediment(i,j) = 0. _d 0
                0913         POCFsediment(i,j) = 0. _d 0
3a1d328575 Oliv*0914         fluxPIC(i,j) = 0. _d 0
                0915         fluxPOC(i,j) = 0. _d 0
                0916         OmegaC0(i,j) = 0. _d 0
2a97638e3d Oliv*0917         radiFluxC(i,j,bi,bj) = 0. _d 0
                0918         radiFluxA(i,j,bi,bj) = 0. _d 0
                0919         radiFluxO(i,j,bi,bj) = 0. _d 0
6aa5674e10 Oliv*0920
                0921        ENDIF
                0922       ENDDO
                0923       ENDDO
c145d7be37 Oliv*0924
                0925 C     Calculate concentration changes at the deepest depth level
                0926 C     and rescale for nonlinear free surface if appropriate
                0927       CALL DARWIN_ADD_BOTFORC(
                0928      U                        gPtr(1-OLx,1-OLy,1,iDIC),
                0929      I                        DICFsediment, oneRL,
                0930      I                        bi, bj, iMin, iMax, jMin, jMax,
                0931      I                        myIter, myTime, myThid )
                0932
                0933       CALL DARWIN_ADD_BOTFORC(
                0934      U                        gPtr(1-OLx,1-OLy,1,iALK),
                0935      I                        ALKFsediment, oneRL,
                0936      I                        bi, bj, iMin, iMax, jMin, jMax,
                0937      I                        myIter, myTime, myThid )
                0938
                0939       CALL DARWIN_ADD_BOTFORC(
                0940      U                        gPtr(1-OLx,1-OLy,1,iO2),
                0941      I                        O2Fsediment, oneRL,
                0942      I                        bi, bj, iMin, iMax, jMin, jMax,
                0943      I                        myIter, myTime, myThid )
                0944
                0945       CALL DARWIN_ADD_BOTFORC(
                0946      U                        gPtr(1-OLx,1-OLy,1,iPIC),
                0947      I                        CALFsediment, oneRL,
                0948      I                        bi, bj, iMin, iMax, jMin, jMax,
                0949      I                        myIter, myTime, myThid )
                0950
                0951       CALL DARWIN_ADD_BOTFORC(
                0952      U                        gPtr(1-OLx,1-OLy,1,iPOC),
                0953      I                        POCFsediment, oneRL,
                0954      I                        bi, bj, iMin, iMax, jMin, jMax,
                0955      I                        myIter, myTime, myThid )
                0956
6aa5674e10 Oliv*0957 #endif
8a10163480 Oliv*0958
                0959 #ifdef DARWIN_ALLOW_RADIv2
                0960       DO j=jMin,jMax
                0961       DO i=iMin,iMax
                0962 C      grid cell just above sea floor
                0963        k = kLowC(i,j,bi,bj)
                0964        IF (k .GT. 0) THEN
                0965
                0966 C       Dependent variables
803871244c Oliv*0967         fluxPIC(i,j) = wPIC_sink*MAX(0 _d 0,Ptracer(i,j,k,bi,bj,iPIC))
                0968         fluxPOC(i,j) = wC_sink*MAX(0 _d 0,Ptracer(i,j,k,bi,bj,iPOC))
                0969         fluxPOP(i,j) = wP_sink*MAX(0 _d 0,Ptracer(i,j,k,bi,bj,iPOP))
                0970         fluxPON(i,j) = wN_sink*MAX(0 _d 0,Ptracer(i,j,k,bi,bj,iPON))
8a10163480 Oliv*0971         absT = ABS(Theta(i,j,k,bi,bj))
                0972
803871244c Oliv*0973 C       Reset flux to zero before accumulation
                0974         flxplnkC(i,j) = 0 _d 0
                0975         flxplnkP(i,j) = 0 _d 0
                0976         flxplnkN(i,j) = 0 _d 0
8a10163480 Oliv*0977         DO l = 1, nplank
803871244c Oliv*0978          flx = biosink(l)*MAX(0 _d 0, Ptracer(i,j,k,bi,bj,ic+l-1))
                0979          flxplnkC(i,j) = flxplnkC(i,j) + flx
                0980          flxplnkP(i,j) = flxplnkP(i,j) + flx*R_PC(l)
                0981          flxplnkN(i,j) = flxplnkN(i,j) + flx*R_NC(l)
                0982         ENDDO
8a10163480 Oliv*0983
                0984         sed_a1 = sed_globala1
                0985         sed_b1 = sed_globalb1
                0986         sed_c1 = sed_globalc1
                0987         sed_d1 = sed_globald1
                0988         sed_e1 = sed_globale1
                0989         sed_a2 = sed_globala2
                0990         sed_b2 = sed_globalb2
                0991         sed_c2 = sed_globalc2
                0992         sed_d2 = sed_globald2
                0993         sed_e2 = sed_globale2
                0994         sed_a3 = sed_globala3
                0995         sed_b3 = sed_globalb3
                0996         sed_c3 = sed_globalc3
                0997         sed_d3 = sed_globald3
                0998         sed_e3 = sed_globale3
                0999         sed_a4 = sed_globala4
                1000         sed_b4 = sed_globalb4
                1001         sed_c4 = sed_globalc4
                1002         sed_d4 = sed_globald4
                1003         sed_e4 = sed_globale4
5938ff335d Oliv*1004         sed_f4 = sed_globalf4
                1005         sed_g4 = sed_globalg4
8a10163480 Oliv*1006         sed_a5 = sed_globala5
                1007         sed_b5 = sed_globalb5
                1008         sed_c5 = sed_globalc5
                1009         sed_d5 = sed_globald5
                1010         sed_e5 = sed_globale5
                1011         sed_a6 = sed_globala6
                1012         sed_b6 = sed_globalb6
                1013         sed_c6 = sed_globalc6
                1014         sed_d6 = sed_globald6
                1015         sed_e6 = sed_globale6
                1016
                1017         DICFsediment(i,j) = sed_a1 +
                1018      &                      sed_b1*absT +
                1019      &                      sed_c1*omegaC0(i,j) +
                1020      &                      sed_d1*(fluxPOC(i,j)+flxplnkC(i,j)) +
                1021      &                      sed_e1*fluxPIC(i,j)
                1022
                1023         ALKFsediment(i,j) = sed_a2 +
                1024      &                      sed_b2*absT +
                1025      &                      sed_c2*omegaC0(i,j) +
                1026      &                      sed_d2*(fluxPOC(i,j)+flxplnkC(i,j)) +
                1027      &                      sed_e2*fluxPIC(i,j)
                1028
                1029         O2Fsediment(i,j)  = sed_a3 +
                1030      &                      sed_b3*absT +
                1031      &                      sed_c3*omegaC0(i,j) +
                1032      &                      sed_d3*(fluxPOC(i,j)+flxplnkC(i,j)) +
                1033      &                      sed_e3*fluxPIC(i,j)
                1034
                1035         NO3Fsediment(i,j) = sed_a4 +
5938ff335d Oliv*1036      &                      sed_b4*(fluxPON(i,j)+flxplnkN(i,j)) +
                1037      &                      sed_c4*DICFsediment(i,j) +
                1038      &                      sed_d4*O2Fsediment(i,j) +
                1039      &                      sed_e4*(absT*DICFsediment(i,j)) +
                1040      &                      sed_f4*(absT*O2Fsediment(i,j)) +
                1041      &                      sed_g4*(DICFsediment(i,j)*O2Fsediment(i,j))
8a10163480 Oliv*1042
                1043         PO4Fsediment(i,j) = sed_a5 +
                1044      &                      sed_b5*absT +
                1045      &                      sed_c5*omegaC0(i,j) +
                1046      &                      sed_d5*(fluxPOP(i,j)+flxplnkP(i,j)) +
                1047      &                      sed_e5*fluxPIC(i,j)
                1048
                1049         NH4Fsediment(i,j) = sed_a6 +
                1050      &                      sed_b6*absT +
                1051      &                      sed_c6*omegaC0(i,j) +
                1052      &                      sed_d6*(fluxPON(i,j)+flxplnkN(i,j)) +
                1053      &                      sed_e6*fluxPIC(i,j)
                1054
                1055 #ifdef DARWIN_ALLOW_CONS
                1056 C       Accumulate for conservation check
                1057         radiFluxC(i,j,bi,bj) = radiFluxC(i,j,bi,bj) +
                1058      &    DICFsediment(i,j)*dTsub(k)*maskC(i,j,k,bi,bj)
                1059         radiFluxA(i,j,bi,bj) = radiFluxA(i,j,bi,bj) +
                1060      &    ALKFsediment(i,j)*dTsub(k)*maskC(i,j,k,bi,bj)
                1061         radiFluxO(i,j,bi,bj) = radiFluxO(i,j,bi,bj) +
                1062      &    O2Fsediment(i,j)*dTsub(k)*maskC(i,j,k,bi,bj)
                1063         radiFluxP(i,j,bi,bj) = radiFluxP(i,j,bi,bj) +
                1064      &    PO4Fsediment(i,j)*dTsub(k)*maskC(i,j,k,bi,bj)
                1065         radiFluxN(i,j,bi,bj) = radiFluxN(i,j,bi,bj) +
                1066      &    (NO3Fsediment(i,j) + NH4Fsediment(i,j)) *
                1067      &    dTsub(k)*maskC(i,j,k,bi,bj)
                1068 #endif
                1069        ELSE
                1070
                1071         DICFsediment(i,j) = 0. _d 0
                1072         ALKFsediment(i,j) = 0. _d 0
                1073         O2Fsediment(i,j) = 0. _d 0
                1074         NO3Fsediment(i,j) = 0. _d 0
                1075         PO4Fsediment(i,j) = 0. _d 0
                1076         NH4Fsediment(i,j) = 0. _d 0
                1077         fluxPIC(i,j) = 0. _d 0
                1078         fluxPOC(i,j) = 0. _d 0
                1079         fluxPOP(i,j) = 0. _d 0
                1080         fluxPON(i,j) = 0. _d 0
                1081         flxplnkC(i,j) = 0. _d 0
                1082         flxplnkN(i,j) = 0. _d 0
                1083         flxplnkP(i,j) = 0. _d 0
                1084
                1085        ENDIF
                1086       ENDDO
                1087       ENDDO
                1088
                1089 C     Calculate concentration changes at the deepest depth level
                1090 C     and rescale for nonlinear free surface if appropriate
                1091       CALL DARWIN_ADD_BOTFORC(
                1092      U                        gPtr(1-OLx,1-OLy,1,iDIC),
                1093      I                        DICFsediment, oneRL,
                1094      I                        bi, bj, iMin, iMax, jMin, jMax,
                1095      I                        myIter, myTime, myThid )
                1096
                1097       CALL DARWIN_ADD_BOTFORC(
                1098      U                        gPtr(1-OLx,1-OLy,1,iALK),
                1099      I                        ALKFsediment, oneRL,
                1100      I                        bi, bj, iMin, iMax, jMin, jMax,
                1101      I                        myIter, myTime, myThid )
                1102
                1103       CALL DARWIN_ADD_BOTFORC(
                1104      U                        gPtr(1-OLx,1-OLy,1,iO2),
                1105      I                        O2Fsediment, oneRL,
                1106      I                        bi, bj, iMin, iMax, jMin, jMax,
                1107      I                        myIter, myTime, myThid )
                1108
                1109       CALL DARWIN_ADD_BOTFORC(
                1110      U                        gPtr(1-OLx,1-OLy,1,iNO3),
                1111      I                        NO3Fsediment, oneRL,
                1112      I                        bi, bj, iMin, iMax, jMin, jMax,
                1113      I                        myIter, myTime, myThid )
                1114
                1115       CALL DARWIN_ADD_BOTFORC(
                1116      U                        gPtr(1-OLx,1-OLy,1,iPO4),
                1117      I                        PO4Fsediment, oneRL,
                1118      I                        bi, bj, iMin, iMax, jMin, jMax,
                1119      I                        myIter, myTime, myThid )
                1120
                1121       CALL DARWIN_ADD_BOTFORC(
                1122      U                        gPtr(1-OLx,1-OLy,1,iNH4),
                1123      I                        NH4Fsediment, oneRL,
                1124      I                        bi, bj, iMin, iMax, jMin, jMax,
                1125      I                        myIter, myTime, myThid )
                1126 #endif
6aa5674e10 Oliv*1127 #endif
                1128
8fbfd1f382 Oliv*1129 C === sinking ==========================================================
                1130       CALL TIMER_START('DARWIN_SINKING [DARWIN_FORCING]',myThid)
5e7acb36b1 daat*1131       CALL DARWIN_SINKING( Ptracer,bioswimDVMup,bioswimDVMdn,
                1132      &                     gPtr,dTsub,bi,bj,
                1133      &                     iMin,iMax,jMin,jMax,subTime,myIter,myThid )
8fbfd1f382 Oliv*1134       CALL TIMER_STOP ('DARWIN_SINKING [DARWIN_FORCING]',myThid)
                1135
                1136 C === apply tendencies to tracers ======================================
                1137       CALL TIMER_START('DARWIN_STEP [DARWIN_FORCING]',myThid)
                1138       DO iTr=1,nDarwin
                1139       DO k=1,Nr
                1140       DO j=jMin,jMax
                1141       DO i=iMin,iMax
                1142         pTracer(i,j,k,bi,bj,iTr)=pTracer(i,j,k,bi,bj,iTr)
                1143      &                      +dTsub(k)*gPtr(i,j,k,iTr)*maskInC(i,j,bi,bj)
                1144       ENDDO
                1145       ENDDO
                1146       ENDDO
                1147       ENDDO
                1148       CALL TIMER_STOP ('DARWIN_STEP [DARWIN_FORCING]',myThid)
                1149
                1150 C === iron =============================================================
                1151 C     re-apply free iron limit to FeT
f48032ad03 Oliv*1152 #ifdef DARWIN_MINFE
8fbfd1f382 Oliv*1153       CALL TIMER_START('DARWIN_FE_CHEM [DARWIN_FORCING]',myThid)
                1154       CALL DARWIN_FE_CHEM(
                1155      U                 pTracer(1-OLx,1-OLy,1,bi,bj,iFeT),
                1156      O                 freeFe(1-OLx,1-OLy,1),
63ffa46034 Oliv*1157      U                 DARWIN_minFeLoss(1,1,1,bi,bj),
8fbfd1f382 Oliv*1158      I                 bi, bj, iMin, iMax, jMin, jMax, myThid)
                1159       CALL TIMER_STOP ('DARWIN_FE_CHEM [DARWIN_FORCING]',myThid)
f48032ad03 Oliv*1160 #endif
8fbfd1f382 Oliv*1161
                1162 C === diagnostics ======================================================
                1163 #ifdef ALLOW_DIAGNOSTICS
                1164       IF (useDIAGNOSTICS) THEN
                1165        CALL TIMER_START('DIAGS_FILL [DARWIN_FORCING]',myThid)
                1166        DO l = 1, nlam
                1167         WRITE(diagname, '(A,I3.3)') 'PAR', l
                1168         CALL DIAGNOSTICS_FILL(PAR(1-OLx,1-OLy,1,l),diagname,0,Nr,2,
5e411acc9e Oliv*1169      &                        bi,bj,myThid)
8fbfd1f382 Oliv*1170        ENDDO
                1171        IF (DIAGNOSTICS_IS_ON('PAR     ', myThid)) THEN
                1172         DO l=2,nlam
                1173          DO k=1,Nr
                1174           DO j=1,sNy
                1175            DO i=1,sNx
                1176             PAR(i,j,k,1) = PAR(i,j,k,1) + PAR(i,j,k,l)
                1177            ENDDO
                1178           ENDDO
                1179          ENDDO
                1180         ENDDO
                1181         WRITE(diagname, '(A)') 'PAR'
                1182         CALL DIAGNOSTICS_FILL(PAR,diagname,0,Nr,2,bi,bj,myThid)
                1183        ENDIF
865e47606a Oliv*1184 #ifndef DARWIN_ALLOW_CHLQUOTA
                1185        IF (DIAGNOSTICS_IS_ON('Chl     ', myThid)) THEN
                1186 #ifdef ALLOW_RADTRANS
                1187         DO k=1,Nr
                1188          DO j=1,sNy
                1189           DO i=1,sNx
                1190            tmp3d(i,j,k) = 0 _d 0
                1191           ENDDO
                1192          ENDDO
                1193         ENDDO
4f40d93b36 Oliv*1194         DO l=1,nPhoto
865e47606a Oliv*1195          DO k=1,Nr
                1196           DO j=1,sNy
                1197            DO i=1,sNx
                1198             tmp3d(i,j,k) = tmp3d(i,j,k) + chlPrev(i,j,k,bi,bj,l)
                1199            ENDDO
                1200           ENDDO
                1201          ENDDO
                1202         ENDDO
c2ba01bd90 Oliv*1203         CALL DIAGNOSTICS_FILL(tmp3d,'Chl     ',0,Nr,2,bi,bj,myThid)
865e47606a Oliv*1204 #else
                1205         CALL DIAGNOSTICS_FILL(chlPrev,'Chl     ',0,Nr,1,bi,bj,myThid)
                1206 #endif
                1207        ENDIF
                1208 #endif
8fbfd1f382 Oliv*1209        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPP),   'PP      ',
                1210      &       0,Nr,2,bi,bj,myThid)
                1211        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNfix), 'Nfix    ',
                1212      &       0,Nr,2,bi,bj,myThid)
                1213        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDenit),'Denit   ',
                1214      &       0,Nr,2,bi,bj,myThid)
3179881a2b Oliv*1215        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDenitN),'DenitN  ',
                1216      &       0,Nr,2,bi,bj,myThid)
c7b6c66d45 Oliv*1217 #ifdef DARWIN_ALLOW_CSTORE
                1218        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iEX),   'EXU     ',
                1219      &       0,Nr,2,bi,bj,myThid)
                1220        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iGW),   'BioSyn  ',
                1221      &       0,Nr,2,bi,bj,myThid)
                1222        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDN),   'DmdN    ',
                1223      &       0,Nr,2,bi,bj,myThid)
                1224        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDP),   'DmdP    ',
                1225      &       0,Nr,2,bi,bj,myThid)
                1226        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDFe),  'DmdFe   ',
                1227      &       0,Nr,2,bi,bj,myThid)
                1228        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDSi),  'DmdSi   ',
                1229      &       0,Nr,2,bi,bj,myThid)
                1230        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDmin), 'Dmdmin  ',
                1231      &       0,Nr,2,bi,bj,myThid)
                1232 #endif
8fbfd1f382 Oliv*1233
                1234 C      nutrient consumption diagnostics
                1235        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsDIN),'C_DIN   ',
b61dca872c Oliv*1236      &       0,Nr,2,bi,bj,myThid)
                1237        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsNO3),'C_NO3   ',
                1238      &       0,Nr,2,bi,bj,myThid)
                1239        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsNO2),'C_NO2   ',
                1240      &       0,Nr,2,bi,bj,myThid)
                1241        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsNH4),'C_NH4   ',
8fbfd1f382 Oliv*1242      &       0,Nr,2,bi,bj,myThid)
                1243        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsPO4),'C_PO4   ',
                1244      &       0,Nr,2,bi,bj,myThid)
b61dca872c Oliv*1245        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsSi), 'C_Si    ',
                1246      &       0,Nr,2,bi,bj,myThid)
                1247        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsFe), 'C_Fe    ',
                1248      &       0,Nr,2,bi,bj,myThid)
                1249        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsDIC),'C_DIC   ',
                1250      &       0,Nr,2,bi,bj,myThid)
                1251        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsDIC_PIC),
                1252      &           'C_DICPIC',0,Nr,2,bi,bj,myThid)
                1253        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iRespirDIC),'respDIC ',
8fbfd1f382 Oliv*1254      &       0,Nr,2,bi,bj,myThid)
b61dca872c Oliv*1255        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iReminDIC_DOC),
                1256      &           'rDIC_DOC',0,Nr,2,bi,bj,myThid)
                1257        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iReminDIC_POC),
                1258      &           'rDIC_POC',0,Nr,2,bi,bj,myThid)
                1259        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDisscDIC_PIC),
                1260      &           'dDIC_PIC',0,Nr,2,bi,bj,myThid)
b6fafcd98a Oliv*1261       CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ifIphavg),'fIph_avg ',
                1262      &       0,Nr,2,bi,bj,myThid)
                1263       CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ifTphavg),'fTph_avg ',
                1264      &       0,Nr,2,bi,bj,myThid)
                1265       CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ifnutavg),'fnut_avg ',
                1266      &       0,Nr,2,bi,bj,myThid)
b61dca872c Oliv*1267 #ifdef DARWIN_ALLOW_CARBON
                1268        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsALK),'C_ALK   ',
                1269      &       0,Nr,2,bi,bj,myThid)
                1270        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsO2), 'C_O2    ',
                1271      &       0,Nr,2,bi,bj,myThid)
                1272 #endif
8fbfd1f382 Oliv*1273
                1274 C      compute 'source' diagnostics: tendencies without consumption;
                1275 C      add full tendencies to consumption and store result back into
                1276 C      diags(iCons*) for convenience
                1277        DO k=1,Nr
                1278         DO j=1,sNy
                1279          DO i=1,sNx
                1280             diags(i,j,k,iConsDIN) = diags(i,j,k,iConsDIN)
                1281      &         + gPtr(i,j,k,iNH4) + gPtr(i,j,k,iNO2) + gPtr(i,j,k,iNO3)
b61dca872c Oliv*1282             diags(i,j,k,iConsNO3) = diags(i,j,k,iConsNO3)
                1283      &                            + gPtr(i,j,k,iNO3)
                1284             diags(i,j,k,iConsNO2) = diags(i,j,k,iConsNO2)
                1285      &                            + gPtr(i,j,k,iNO2)
                1286             diags(i,j,k,iConsNH4) = diags(i,j,k,iConsNH4)
                1287      &                            + gPtr(i,j,k,iNH4)
8fbfd1f382 Oliv*1288             diags(i,j,k,iConsPO4) = diags(i,j,k,iConsPO4)
                1289      &                            + gPtr(i,j,k,iPO4)
                1290             diags(i,j,k,iConsSi) = diags(i,j,k,iConsSi)
                1291      &                           + gPtr(i,j,k,iSiO2)
                1292             diags(i,j,k,iConsFe) = diags(i,j,k,iConsFe)
                1293      &                           + gPtr(i,j,k,iFeT)
b61dca872c Oliv*1294 #ifdef DARWIN_ALLOW_CARBON
                1295             diags(i,j,k,iConsALK) = diags(i,j,k,iConsALK)
                1296      &                            + gPtr(i,j,k,iALK)
                1297             diags(i,j,k,iConsO2)  = diags(i,j,k,iConsO2)
                1298      &                            + gPtr(i,j,k,iO2)
                1299 #endif
8fbfd1f382 Oliv*1300          ENDDO
                1301         ENDDO
                1302        ENDDO
                1303        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsDIN),'S_DIN   ',
b61dca872c Oliv*1304      &       0,Nr,2,bi,bj,myThid)
                1305        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsNO3),'S_NO3   ',
                1306      &       0,Nr,2,bi,bj,myThid)
                1307        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsNO2),'S_NO2   ',
                1308      &       0,Nr,2,bi,bj,myThid)
                1309        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsNH4),'S_NH4   ',
8fbfd1f382 Oliv*1310      &       0,Nr,2,bi,bj,myThid)
                1311        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsPO4),'S_PO4   ',
                1312      &       0,Nr,2,bi,bj,myThid)
b61dca872c Oliv*1313        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsSi), 'S_Si    ',
                1314      &       0,Nr,2,bi,bj,myThid)
                1315        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsFe), 'S_Fe    ',
8fbfd1f382 Oliv*1316      &       0,Nr,2,bi,bj,myThid)
b61dca872c Oliv*1317 #ifdef DARWIN_ALLOW_CARBON
                1318        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsALK),'S_ALK   ',
                1319      &       0,Nr,2,bi,bj,myThid)
                1320        CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iConsO2), 'S_O2    ',
8fbfd1f382 Oliv*1321      &       0,Nr,2,bi,bj,myThid)
b61dca872c Oliv*1322 #endif
8fbfd1f382 Oliv*1323
                1324        DO iTr=1,nDarwin
                1325         diagname = 'gDAR'//PTRACERS_ioLabel(iTr)
                1326         CALL DIAGNOSTICS_FILL(gPtr(1-OLx,1-OLy,1,iTr), diagname,
                1327      &          0,Nr,2,bi,bj,myThid)
                1328        ENDDO
b61dca872c Oliv*1329
ee4178dfbe Oliv*1330 #ifdef DARWIN_DIAG_PERTYPE
                1331        DO iTr=1,nplank
8fbfd1f382 Oliv*1332         WRITE(diagname, '(A,I4.4)') 'PP', iTr
                1333         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPPplank+iTr-1),
d5f8f1c735 Oliv*1334      &        diagname,0,Nr,2,bi,bj,myThid)
                1335         WRITE(diagname, '(A,I4.4)') 'PC', iTr
                1336         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPCplank+iTr-1),
a93256cee4 Oliv*1337      &        diagname,0,Nr,2,bi,bj,myThid)
                1338         WRITE(diagname, '(A,I4.4)') 'HP', iTr
                1339         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iHPplank+iTr-1),
                1340      &        diagname,0,Nr,2,bi,bj,myThid)
                1341         WRITE(diagname, '(A,I4.4)') 'HC', iTr
                1342         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iHCplank+iTr-1),
8fbfd1f382 Oliv*1343      &        diagname,0,Nr,2,bi,bj,myThid)
                1344         WRITE(diagname, '(A,I4.4)') 'GR', iTr
                1345         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iGRplank+iTr-1),
73aabcca42 Oliv*1346      &        diagname,0,Nr,2,bi,bj,myThid)
                1347         WRITE(diagname, '(A,I4.4)') 'GrGn', iTr
                1348         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iGrGn+iTr-1),
a93256cee4 Oliv*1349      &        diagname,0,Nr,2,bi,bj,myThid)
                1350         WRITE(diagname, '(A,I4.4)') 'GrGC', iTr
                1351         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iGrGC+iTr-1),
fd20fa810c Oliv*1352      &        diagname,0,Nr,2,bi,bj,myThid)
                1353         WRITE(diagname, '(A,I4.4)') 'Mort', iTr
                1354         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iMort+iTr-1),
                1355      &        diagname,0,Nr,2,bi,bj,myThid)
                1356         WRITE(diagname, '(A,I4.4)') 'Resp', iTr
                1357         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iResp+iTr-1),
                1358      &        diagname,0,Nr,2,bi,bj,myThid)
                1359         WRITE(diagname, '(A,I4.4)') 'fnut', iTr
                1360         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ifnut+iTr-1),
                1361      &        diagname,0,Nr,2,bi,bj,myThid)
                1362         WRITE(diagname, '(A,I4.4)') 'fIph', iTr
                1363         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ifIph+iTr-1),
                1364      &        diagname,0,Nr,2,bi,bj,myThid)
                1365         WRITE(diagname, '(A,I4.4)') 'fTph', iTr
                1366         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ifTph+iTr-1),
                1367      &        diagname,0,Nr,2,bi,bj,myThid)
                1368         WRITE(diagname, '(A,I4.4)') 'limN', iTr
                1369         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ilimN+iTr-1),
                1370      &        diagname,0,Nr,2,bi,bj,myThid)
                1371         WRITE(diagname, '(A,I4.4)') 'limP', iTr
                1372         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ilimP+iTr-1),
                1373      &        diagname,0,Nr,2,bi,bj,myThid)
                1374         WRITE(diagname, '(A,I4.4)') 'limF', iTr
                1375         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ilimF+iTr-1),
                1376      &        diagname,0,Nr,2,bi,bj,myThid)
                1377         WRITE(diagname, '(A,I4.4)') 'limS', iTr
                1378         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ilimS+iTr-1),
8fbfd1f382 Oliv*1379      &        diagname,0,Nr,2,bi,bj,myThid)
                1380        ENDDO
a092808e6b shlo*1381 # ifdef DARWIN_MACROMOLECULAR_GROWTH
                1382        DO iTr=1,nPhoto
                1383         WRITE(diagname, '(A,I4.4)') 'PChl', iTr
                1384         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPChl+iTr-1),diagname,
                1385      &          0,Nr,2,bi,bj,myThid)
                1386         WRITE(diagname, '(A,I4.4)') 'VN', iTr
                1387         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iVN+iTr-1),diagname,
                1388      &          0,Nr,2,bi,bj,myThid)
                1389         WRITE(diagname, '(A,I4.4)') 'VP', iTr
                1390         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iVP+iTr-1),diagname,
                1391      &          0,Nr,2,bi,bj,myThid)
                1392         WRITE(diagname, '(A,I4.4)') 'MODE', iTr
                1393         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iMODE+iTr-1),diagname,
                1394      &          0,Nr,2,bi,bj,myThid)
                1395         WRITE(diagname, '(A,I4.4)') 'Fe_C', iTr
                1396         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iFe_C+iTr-1),diagname,
                1397      &          0,Nr,2,bi,bj,myThid)
                1398         WRITE(diagname, '(A,I4.4)') 'exQc', iTr
                1399         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iexQc+iTr-1),diagname,
                1400      &          0,Nr,2,bi,bj,myThid)
                1401         WRITE(diagname, '(A,I4.4)') 'CChl', iTr
                1402         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iCChl+iTr-1),diagname,
                1403      &          0,Nr,2,bi,bj,myThid)
                1404         WRITE(diagname, '(A,I4.4)') 'NChl', iTr
                1405         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNChl+iTr-1),diagname,
                1406      &          0,Nr,2,bi,bj,myThid)
                1407         WRITE(diagname, '(A,I4.4)') 'NPho', iTr
                1408         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNPho+iTr-1),diagname,
                1409      &          0,Nr,2,bi,bj,myThid)
                1410         WRITE(diagname, '(A,I4.4)') 'NSyn', iTr
                1411         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNSyn+iTr-1),diagname,
                1412      &          0,Nr,2,bi,bj,myThid)
                1413         WRITE(diagname, '(A,I4.4)') 'NPrn', iTr
                1414         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNPrn+iTr-1),diagname,
                1415      &          0,Nr,2,bi,bj,myThid)
                1416         WRITE(diagname, '(A,I4.4)') 'NRNA', iTr
                1417         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNRNA+iTr-1),diagname,
                1418      &          0,Nr,2,bi,bj,myThid)
                1419         WRITE(diagname, '(A,I4.4)') 'NDNA', iTr
                1420         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNDNA+iTr-1),diagname,
                1421      &          0,Nr,2,bi,bj,myThid)
                1422         WRITE(diagname, '(A,I4.4)') 'NSTO', iTr
                1423         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNSTO+iTr-1),diagname,
                1424      &          0,Nr,2,bi,bj,myThid)
                1425         WRITE(diagname, '(A,I4.4)') 'NEXC', iTr
                1426         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iNEXC+iTr-1),diagname,
                1427      &          0,Nr,2,bi,bj,myThid)
                1428         WRITE(diagname, '(A,I4.4)') 'PRNA', iTr
                1429         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPRNA+iTr-1),diagname,
                1430      &          0,Nr,2,bi,bj,myThid)
                1431         WRITE(diagname, '(A,I4.4)') 'PDNA', iTr
                1432         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPDNA+iTr-1),diagname,
                1433      &          0,Nr,2,bi,bj,myThid)
                1434         WRITE(diagname, '(A,I4.4)') 'PTHY', iTr
                1435         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPTHY+iTr-1),diagname,
                1436      &          0,Nr,2,bi,bj,myThid)
                1437         WRITE(diagname, '(A,I4.4)') 'PCON', iTr
                1438         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPCON+iTr-1),diagname,
                1439      &          0,Nr,2,bi,bj,myThid)
                1440         WRITE(diagname, '(A,I4.4)') 'PSTO', iTr
                1441         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPSTO+iTr-1),diagname,
                1442      &          0,Nr,2,bi,bj,myThid)
                1443         WRITE(diagname, '(A,I4.4)') 'PEXC', iTr
                1444         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iPEXC+iTr-1),diagname,
                1445      &          0,Nr,2,bi,bj,myThid)
                1446         WRITE(diagname, '(A,I4.4)') 'FPHO', iTr
                1447         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iFPHO+iTr-1),diagname,
                1448      &          0,Nr,2,bi,bj,myThid)
                1449         WRITE(diagname, '(A,I4.4)') 'FSTO', iTr
                1450         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iFSTO+iTr-1),diagname,
                1451      &          0,Nr,2,bi,bj,myThid)
                1452         WRITE(diagname, '(A,I4.4)') 'FEXC', iTr
                1453         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iFEXC+iTr-1),diagname,
                1454      &          0,Nr,2,bi,bj,myThid)
                1455         WRITE(diagname, '(A,I4.4)') 'Y_RQ', iTr
                1456         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iY_RQ+iTr-1),diagname,
                1457      &          0,Nr,2,bi,bj,myThid)
                1458         WRITE(diagname, '(A,I4.4)') 'limC', iTr
                1459         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ilimC+iTr-1),diagname,
                1460      &          0,Nr,2,bi,bj,myThid)
                1461         WRITE(diagname, '(A,I4.4)') 'limL', iTr
                1462         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,ilimL+iTr-1),diagname,
                1463      &          0,Nr,2,bi,bj,myThid)
                1464        ENDDO
                1465 # endif
c7b6c66d45 Oliv*1466 #endif
                1467 #ifdef DARWIN_ALLOW_CSTORE
                1468 #ifdef DARWIN_ALLOW_CSTORE_DIAGS
                1469        DO iTr=1,nPhoto
                1470         WRITE(diagname, '(A,I4.4)') 'EXU', iTr
                1471         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iEXplank+iTr-1),
                1472      &        diagname,0,Nr,2,bi,bj,myThid)
                1473         WRITE(diagname, '(A,I4.4)') 'BS', iTr
                1474         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iGWplank+iTr-1),
                1475      &        diagname,0,Nr,2,bi,bj,myThid)
                1476         WRITE(diagname, '(A,I4.4)') 'DN', iTr
                1477         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDNplank+iTr-1),
                1478      &        diagname,0,Nr,2,bi,bj,myThid)
                1479         WRITE(diagname, '(A,I4.4)') 'DP', iTr
                1480         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDPplank+iTr-1),
                1481      &        diagname,0,Nr,2,bi,bj,myThid)
                1482         WRITE(diagname, '(A,I4.4)') 'DFe', iTr
                1483         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDFplank+iTr-1),
                1484      &        diagname,0,Nr,2,bi,bj,myThid)
                1485         WRITE(diagname, '(A,I4.4)') 'DSi', iTr
                1486         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDSplank+iTr-1),
                1487      &        diagname,0,Nr,2,bi,bj,myThid)
                1488         WRITE(diagname, '(A,I4.4)') 'Dmin', iTr
                1489         CALL DIAGNOSTICS_FILL(diags(1-OLx,1-OLy,1,iDminplank+iTr-1),
                1490      &        diagname,0,Nr,2,bi,bj,myThid)
                1491        ENDDO
                1492 #endif
ee4178dfbe Oliv*1493 #endif
82fd1e0266 Oliv*1494        CALL DIAGNOSTICS_FILL(scavRate,'scavRate',0,Nr,2,bi,bj,myThid)
fd32b5ab01 Oliv*1495        CALL DIAGNOSTICS_FILL(scavLoss,'scvLosFe',0,Nr,2,bi,bj,myThid)
342416534c Oliv*1496        CALL DIAGNOSTICS_FILL(sedFe,   'sedFe   ',0,Nr,2,bi,bj,myThid)
c8be5f4655 Oliv*1497        CALL DIAGNOSTICS_FILL(sedFlxFe,'sedFlxFe',0,1,2,bi,bj,myThid)
3997029bbb Oliv*1498        IF (DIAGNOSTICS_IS_ON('sfcSolFe', myThid)) THEN
                1499         DO j=jMin,jMax
                1500          DO i=iMin,iMax
                1501           IF (hFacC(i,j,1,bi,bj) .EQ. 0.) THEN
                1502            surfFe(i,j) = 0 _d 0
                1503           ELSE
                1504            surfFe(i,j) = alpfe * inputFe(i,j,bi,bj)
                1505           ENDIF
                1506          ENDDO
                1507         ENDDO
                1508         CALL DIAGNOSTICS_FILL(surfFe,'sfcSolFe',0,1,2,bi,bj,myThid)
                1509        ENDIF
b61dca872c Oliv*1510
8fbfd1f382 Oliv*1511 #ifdef DARWIN_ALLOW_CARBON
                1512        CALL DIAGNOSTICS_FILL(gDIC,   'gDICsurf',0,1,2,bi,bj,myThid)
                1513        CALL DIAGNOSTICS_FILL(gALK,   'gALKsurf',0,1,2,bi,bj,myThid)
                1514        CALL DIAGNOSTICS_FILL(gO2,    'gO2surf ',0,1,2,bi,bj,myThid)
8a10163480 Oliv*1515        CALL DIAGNOSTICS_FILL(OmegaC0,     'OmegCbot',0,1,2,bi,bj,myThid)
                1516 #ifdef DARWIN_ALLOW_RADIv1
6aa5674e10 Oliv*1517        CALL DIAGNOSTICS_FILL(DICFsediment,'DICFsed ',0,1,2,bi,bj,myThid)
                1518        CALL DIAGNOSTICS_FILL(ALKFsediment,'ALKFsed ',0,1,2,bi,bj,myThid)
                1519        CALL DIAGNOSTICS_FILL(O2Fsediment, 'O2Fsed  ',0,1,2,bi,bj,myThid)
                1520        CALL DIAGNOSTICS_FILL(POCFsediment,'POCFbur ',0,1,2,bi,bj,myThid)
                1521        CALL DIAGNOSTICS_FILL(CALFsediment,'CALFbur ',0,1,2,bi,bj,myThid)
                1522        CALL DIAGNOSTICS_FILL(fluxPIC,     'sedFlPIC',0,1,2,bi,bj,myThid)
                1523        CALL DIAGNOSTICS_FILL(fluxPOC,     'sedFlPOC',0,1,2,bi,bj,myThid)
                1524        CALL DIAGNOSTICS_FILL(TcorrO2,     'TcorrO2 ',0,1,2,bi,bj,myThid)
                1525        CALL DIAGNOSTICS_FILL(TcorrDIC,    'TcorrDIC',0,1,2,bi,bj,myThid)
                1526        CALL DIAGNOSTICS_FILL(TcorrALK,    'TcorrALK',0,1,2,bi,bj,myThid)
                1527        CALL DIAGNOSTICS_FILL(OmegaC0,     'OmegCbot',0,1,2,bi,bj,myThid)
                1528 #endif
8a10163480 Oliv*1529 #ifdef DARWIN_ALLOW_RADIv2
                1530        CALL DIAGNOSTICS_FILL(DICFsediment,'DICFsed ',0,1,2,bi,bj,myThid)
                1531        CALL DIAGNOSTICS_FILL(ALKFsediment,'ALKFsed ',0,1,2,bi,bj,myThid)
                1532        CALL DIAGNOSTICS_FILL(O2Fsediment, 'O2Fsed  ',0,1,2,bi,bj,myThid)
                1533        CALL DIAGNOSTICS_FILL(NO3Fsediment,'NO3Fsed ',0,1,2,bi,bj,myThid)
                1534        CALL DIAGNOSTICS_FILL(PO4Fsediment,'PO4Fsed ',0,1,2,bi,bj,myThid)
                1535        CALL DIAGNOSTICS_FILL(NH4Fsediment,'NH4Fsed ',0,1,2,bi,bj,myThid)
                1536        CALL DIAGNOSTICS_FILL(fluxPIC,     'sedFlPIC',0,1,2,bi,bj,myThid)
                1537        CALL DIAGNOSTICS_FILL(fluxPOC,     'sedFlPOC',0,1,2,bi,bj,myThid)
                1538        CALL DIAGNOSTICS_FILL(fluxPOP,     'sedFlPOP',0,1,2,bi,bj,myThid)
                1539        CALL DIAGNOSTICS_FILL(fluxPON,     'sedFlPON',0,1,2,bi,bj,myThid)
                1540        CALL DIAGNOSTICS_FILL(flxplnkC,    'sedFplkC',0,1,2,bi,bj,myThid)
                1541        CALL DIAGNOSTICS_FILL(flxplnkP,    'sedFplkP',0,1,2,bi,bj,myThid)
                1542        CALL DIAGNOSTICS_FILL(flxplnkN,    'sedFplkN',0,1,2,bi,bj,myThid)
                1543 #endif
8fbfd1f382 Oliv*1544 #endif
714cf988dd Oliv*1545
                1546 #ifdef DARWIN_DIAG_TENDENCIES
                1547        CALL DIAGNOSTICS_FILL(gDICsurfForc,'gDICEpr ',0,1,2,
                1548      &  bi,bj,myThid)
                1549        CALL DIAGNOSTICS_FILL(gNO3surfForc,'gNO3Epr ',0,1,2,
                1550      &  bi,bj,myThid)
                1551        CALL DIAGNOSTICS_FILL(gNO2surfForc,'gNO2Epr ',0,1,2,
                1552      &  bi,bj,myThid)
                1553        CALL DIAGNOSTICS_FILL(gNH4surfForc,'gNH4Epr ',0,1,2,
                1554      &  bi,bj,myThid)
                1555        CALL DIAGNOSTICS_FILL(gPO4surfForc,'gPO4Epr ',0,1,2,
                1556      &  bi,bj,myThid)
                1557        CALL DIAGNOSTICS_FILL(gFeTsurfForc,'gFeTEpr ',0,1,2,
                1558      &  bi,bj,myThid)
                1559        CALL DIAGNOSTICS_FILL(gSiO2surfForc,'gSiO2Epr',0,1,2,
                1560      &  bi,bj,myThid)
                1561 #ifdef DARWIN_ALLOW_CARBON
                1562        CALL DIAGNOSTICS_FILL(gALKsurfForc,'gALKEpr ',0,1,2,
                1563      &  bi,bj,myThid)
                1564        CALL DIAGNOSTICS_FILL(gO2surfForc, 'gO2Epr  ',0,1,2,
                1565      &  bi,bj,myThid)
                1566 #endif
                1567 #endif
                1568
8fbfd1f382 Oliv*1569        CALL TIMER_STOP ('DIAGS_FILL [DARWIN_FORCING]',myThid)
714cf988dd Oliv*1570 C     useDiagnostics
8fbfd1f382 Oliv*1571       ENDIF
714cf988dd Oliv*1572 #endif /* ALLOW_DIAGNOSTICS */
8fbfd1f382 Oliv*1573
                1574       midTime = midTime + deltaTclock/nsubtime
                1575       subTime = subTime + deltaTclock/nsubtime
                1576 C     isub
                1577       ENDDO
                1578
63ffa46034 Oliv*1579 #ifdef ALLOW_DIAGNOSTICS
                1580       IF (useDIAGNOSTICS) THEN
                1581        IF (DIAGNOSTICS_IS_ON('freeFeLs', myThid)) THEN
                1582         DO k=1,Nr
                1583          tmpFac = 1 _d 0 / PTRACERS_dTLev(k)
                1584          DO j=1,sNy
0f72dca88c Oliv*1585           DO i=1,sNx
63ffa46034 Oliv*1586            tmp3d(i,j,k) = tmpFac*DARWIN_minFeLoss(i,j,k,bi,bj)
                1587           ENDDO
                1588          ENDDO
                1589         ENDDO
0f72dca88c Oliv*1590         CALL DIAGNOSTICS_FILL(tmp3d,'freeFeLs',0,Nr,2,bi,bj,myThid)
63ffa46034 Oliv*1591        ENDIF
                1592       ENDIF
                1593 #endif
                1594
8fbfd1f382 Oliv*1595 #endif /* ALLOW_DARWIN */
                1596
                1597       RETURN
                1598       END