Back to home page

darwin3

 
 

    


File indexing completed on 2025-09-13 12:07:41 UTC

view on githubraw file Latest commit 5e7acb36 on 2025-07-11 21:23:07 UTC
8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C !ROUTINE: DARWIN_SINKING
                0005 C !INTERFACE: ==========================================================
                0006       SUBROUTINE DARWIN_SINKING(
5e7acb36b1 daat*0007      I     Ptr,bioswimDVMup,bioswimDVMdn,
8fbfd1f382 Oliv*0008      U     gTr,
248275f1c4 Oliv*0009      I     dTsub,bi,bj,iMin,iMax,jMin,jMax,myTime,myIter,myThid)
8fbfd1f382 Oliv*0010 
                0011 C !DESCRIPTION:
                0012 C     compute tendencies from sinking of particulate organic matter
                0013 
                0014 C !USES: ===============================================================
                0015       IMPLICIT NONE
                0016 #include "SIZE.h"
dcd877c466 Oliv*0017 #include "EEPARAMS.h"
                0018 #include "PARAMS.h"
8fbfd1f382 Oliv*0019 #include "GRID.h"
dcd877c466 Oliv*0020 #include "SURFACE.h"
8fbfd1f382 Oliv*0021 #include "DARWIN_SIZE.h"
                0022 #include "DARWIN_INDICES.h"
                0023 #include "DARWIN_PARAMS.h"
                0024 #include "DARWIN_TRAITS.h"
248275f1c4 Oliv*0025 #include "DARWIN_FIELDS.h"
8fbfd1f382 Oliv*0026 
                0027 C !INPUT PARAMETERS: ===================================================
5e7acb36b1 daat*0028 C myThid       :: thread number
                0029 C Ptr          :: darwin model tracers
                0030 C bioswimDVMup :: extra upward velocity due to diel vertical migration
                0031 C bioswimDVMdn :: extra downward velocity due to diel vertical migration
8fbfd1f382 Oliv*0032       _RL Ptr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nDarwin)
5e7acb36b1 daat*0033       _RL bioswimDVMup(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nplank)
                0034       _RL bioswimDVMdn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nplank)
                0035 
248275f1c4 Oliv*0036       _RL dTsub(Nr)
8fbfd1f382 Oliv*0037       INTEGER bi,bj,iMin,iMax,jMin,jMax
                0038       INTEGER myThid, myIter
                0039       _RL myTime
                0040 
                0041 C !INPUT/OUTPUT PARAMETERS: ============================================
                0042 C  gTr    :: computed tendencies
                0043       _RL gTr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nDarwin)
                0044 CEOP
                0045 
                0046 #ifdef ALLOW_DARWIN
                0047 
                0048 c !LOCAL VARIABLES: ====================================================
                0049       INTEGER i,j,k,l
dcd877c466 Oliv*0050       _RL upfc, dnfc, flux, hFacCdn
                0051       _RL upf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0052       _RL dnf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
8fbfd1f382 Oliv*0053 
                0054       DO k=1,Nr-1
                0055        DO j=jMin,jMax
                0056         DO i=iMin,iMax
b409127548 Oliv*0057          IF (maskC(i,j,k,bi,bj) .NE. zeroRS .AND.
                0058      &       maskC(i,j,k+1,bi,bj) .NE. zeroRS) THEN
dcd877c466 Oliv*0059           upf(i,j,k)=recip_drF(k)*recip_hFacC(i,j,k,bi,bj)
                0060           dnf(i,j,k+1)=recip_drF(k+1)*recip_hFacC(i,j,k+1,bi,bj)
                0061          ELSE
                0062           upf(i,j,k) = 0 _d 0
                0063           dnf(i,j,k+1) = 0 _d 0
                0064          ENDIF
                0065         ENDDO
                0066        ENDDO
                0067 #ifdef NONLIN_FRSURF
                0068 C-    Account for change in level thickness
                0069        IF (nonlinFreeSurf.GT.0) THEN
                0070          CALL FREESURF_RESCALE_G(
                0071      I                            bi, bj, k,
                0072      U                            upf,
                0073      I                            myThid )
                0074          CALL FREESURF_RESCALE_G(
                0075      I                            bi, bj, k+1,
                0076      U                            dnf,
                0077      I                            myThid )
                0078        ENDIF
                0079 #endif /* NONLIN_FRSURF */
                0080        DO j=jMin,jMax
                0081         DO i=iMin,iMax
                0082          upfc = upf(i,j,k)
                0083          dnfc = dnf(i,j,k+1)
                0084          IF (dnfc .GT. 0 _d 0) THEN
086a45f245 Oliv*0085           flux = wPIC_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPIC))
dcd877c466 Oliv*0086           gTr(i,j,k  ,iPIC ) = gTr(i,j,k  ,iPIC ) - flux*upfc
                0087           gTr(i,j,k+1,iPIC ) = gTr(i,j,k+1,iPIC ) + flux*dnfc
086a45f245 Oliv*0088           flux = wC_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOC))
dcd877c466 Oliv*0089           gTr(i,j,k  ,iPOC ) = gTr(i,j,k  ,iPOC ) - flux*upfc
                0090           gTr(i,j,k+1,iPOC ) = gTr(i,j,k+1,iPOC ) + flux*dnfc
086a45f245 Oliv*0091           flux = wN_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPON))
dcd877c466 Oliv*0092           gTr(i,j,k  ,iPON ) = gTr(i,j,k  ,iPON ) - flux*upfc
                0093           gTr(i,j,k+1,iPON ) = gTr(i,j,k+1,iPON ) + flux*dnfc
086a45f245 Oliv*0094           flux = wP_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOP))
dcd877c466 Oliv*0095           gTr(i,j,k  ,iPOP ) = gTr(i,j,k  ,iPOP ) - flux*upfc
                0096           gTr(i,j,k+1,iPOP ) = gTr(i,j,k+1,iPOP ) + flux*dnfc
086a45f245 Oliv*0097           flux = wSi_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOSi))
dcd877c466 Oliv*0098           gTr(i,j,k  ,iPOSi) = gTr(i,j,k  ,iPOSi) - flux*upfc
                0099           gTr(i,j,k+1,iPOSi) = gTr(i,j,k+1,iPOSi) + flux*dnfc
086a45f245 Oliv*0100           flux = wFe_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOFe))
dcd877c466 Oliv*0101           gTr(i,j,k  ,iPOFe) = gTr(i,j,k  ,iPOFe) - flux*upfc
                0102           gTr(i,j,k+1,iPOFe) = gTr(i,j,k+1,iPOFe) + flux*dnfc
8fbfd1f382 Oliv*0103           DO l = 1, nplank
5e7acb36b1 daat*0104            flux = (biosink(l)+bioswimDVMdn(i,j,k,l))*MAX(0.0,
                0105      &           Ptr(i,j,k,bi,bj,ic+l-1))
dcd877c466 Oliv*0106            gTr(i,j,k  ,ic+l-1 )=gTr(i,j,k  ,ic+l-1 ) - flux*upfc
                0107            gTr(i,j,k+1,ic+l-1 )=gTr(i,j,k+1,ic+l-1 ) + flux*dnfc
5e7acb36b1 daat*0108            flux = (bioswim(l)+bioswimDVMup(i,j,k+1,l))*MAX(0.0,
                0109      &           Ptr(i,j,k+1,bi,bj,ic+l-1))
dcd877c466 Oliv*0110            gTr(i,j,k  ,ic+l-1 )=gTr(i,j,k  ,ic+l-1 ) + flux*upfc
                0111            gTr(i,j,k+1,ic+l-1 )=gTr(i,j,k+1,ic+l-1 ) - flux*dnfc
8fbfd1f382 Oliv*0112 #ifdef DARWIN_ALLOW_NQUOTA
5e7acb36b1 daat*0113            flux = (biosink(l)+bioswimDVMdn(i,j,k,l))*MAX(0.0,
                0114      &           Ptr(i,j,k,bi,bj,in+l-1))
dcd877c466 Oliv*0115            gTr(i,j,k  ,in+l-1 )=gTr(i,j,k  ,in+l-1 ) - flux*upfc
                0116            gTr(i,j,k+1,in+l-1 )=gTr(i,j,k+1,in+l-1 ) + flux*dnfc
5e7acb36b1 daat*0117            flux = (bioswim(l)+bioswimDVMup(i,j,k+1,l))*MAX(0.0,
                0118      &           Ptr(i,j,k+1,bi,bj,in+l-1))
dcd877c466 Oliv*0119            gTr(i,j,k  ,in+l-1 )=gTr(i,j,k  ,in+l-1 ) + flux*upfc
                0120            gTr(i,j,k+1,in+l-1 )=gTr(i,j,k+1,in+l-1 ) - flux*dnfc
5e7acb36b1 daat*0121 
8fbfd1f382 Oliv*0122 #endif
                0123 #ifdef DARWIN_ALLOW_PQUOTA
5e7acb36b1 daat*0124            flux = (biosink(l)+bioswimDVMdn(i,j,k,l))*MAX(0.0,
                0125      &           Ptr(i,j,k,bi,bj,ip+l-1))
dcd877c466 Oliv*0126            gTr(i,j,k  ,ip+l-1 )=gTr(i,j,k  ,ip+l-1 ) - flux*upfc
                0127            gTr(i,j,k+1,ip+l-1 )=gTr(i,j,k+1,ip+l-1 ) + flux*dnfc
5e7acb36b1 daat*0128            flux = (bioswim(l)+bioswimDVMup(i,j,k+1,l))*MAX(0.0,
                0129      &           Ptr(i,j,k+1,bi,bj,ip+l-1))
dcd877c466 Oliv*0130            gTr(i,j,k  ,ip+l-1 )=gTr(i,j,k  ,ip+l-1 ) + flux*upfc
                0131            gTr(i,j,k+1,ip+l-1 )=gTr(i,j,k+1,ip+l-1 ) - flux*dnfc
8fbfd1f382 Oliv*0132 #endif
                0133 #ifdef DARWIN_ALLOW_SIQUOTA
5e7acb36b1 daat*0134            flux = (biosink(l)+bioswimDVMdn(i,j,k,l))*MAX(0.0,
                0135      &           Ptr(i,j,k,bi,bj,isi+l-1))
dcd877c466 Oliv*0136            gTr(i,j,k  ,isi+l-1)=gTr(i,j,k  ,isi+l-1) - flux*upfc
                0137            gTr(i,j,k+1,isi+l-1)=gTr(i,j,k+1,isi+l-1) + flux*dnfc
5e7acb36b1 daat*0138            flux = (bioswim(l)+bioswimDVMup(i,j,k+1,l))*MAX(0.0,
                0139      &           Ptr(i,j,k+1,bi,bj,isi+l-1))
dcd877c466 Oliv*0140            gTr(i,j,k  ,isi+l-1)=gTr(i,j,k  ,isi+l-1) + flux*upfc
                0141            gTr(i,j,k+1,isi+l-1)=gTr(i,j,k+1,isi+l-1) - flux*dnfc
8fbfd1f382 Oliv*0142 #endif
                0143 #ifdef DARWIN_ALLOW_FEQUOTA
5e7acb36b1 daat*0144            flux = (biosink(l)+bioswimDVMdn(i,j,k,l))*MAX(0.0,
                0145      &           Ptr(i,j,k,bi,bj,ife+l-1))
dcd877c466 Oliv*0146            gTr(i,j,k  ,ife+l-1)=gTr(i,j,k  ,ife+l-1) - flux*upfc
                0147            gTr(i,j,k+1,ife+l-1)=gTr(i,j,k+1,ife+l-1) + flux*dnfc
5e7acb36b1 daat*0148            flux = (bioswim(l)+bioswimDVMup(i,j,k+1,l))*MAX(0.0,
                0149      &           Ptr(i,j,k+1,bi,bj,ife+l-1))
dcd877c466 Oliv*0150            gTr(i,j,k  ,ife+l-1)=gTr(i,j,k  ,ife+l-1) + flux*upfc
                0151            gTr(i,j,k+1,ife+l-1)=gTr(i,j,k+1,ife+l-1) - flux*dnfc
8fbfd1f382 Oliv*0152 #endif
                0153           ENDDO
                0154 #ifdef DARWIN_ALLOW_CHLQUOTA
                0155           DO l = 1, nPhoto
5e7acb36b1 daat*0156            flux = (biosink(l)+bioswimDVMdn(i,j,k,l))*MAX(0.0,
                0157      &           Ptr(i,j,k,bi,bj,iChl+l-1))
dcd877c466 Oliv*0158            gTr(i,j,k  ,iChl+l-1)=gTr(i,j,k  ,iChl+l-1)-flux*upfc
                0159            gTr(i,j,k+1,iChl+l-1)=gTr(i,j,k+1,iChl+l-1)+flux*dnfc
5e7acb36b1 daat*0160            flux = (bioswim(l)+bioswimDVMup(i,j,k+1,l))*MAX(0.0,
                0161      &           Ptr(i,j,k+1,bi,bj,iChl+l-1))
dcd877c466 Oliv*0162            gTr(i,j,k  ,iChl+l-1)=gTr(i,j,k  ,iChl+l-1)+flux*upfc
                0163            gTr(i,j,k+1,iChl+l-1)=gTr(i,j,k+1,iChl+l-1)-flux*dnfc
8fbfd1f382 Oliv*0164           ENDDO
                0165 #endif
                0166          ENDIF
                0167         ENDDO
                0168        ENDDO
                0169       ENDDO
                0170 
1fb46047a0 Oliv*0171 #ifdef DARWIN_BOTTOM_SINK
                0172       DO k=1,Nr
                0173        DO j=jMin,jMax
                0174         DO i=iMin,iMax
                0175          IF (k.LT.Nr) THEN
dcd877c466 Oliv*0176           hFacCdn = hFacC(i,j,k+1,bi,bj)
1fb46047a0 Oliv*0177          ELSE
dcd877c466 Oliv*0178           hFacCdn = 0 _d 0
1fb46047a0 Oliv*0179          ENDIF
dcd877c466 Oliv*0180          IF (hFacC(i,j,k,bi,bj).GT.0 _d 0 .AND. hFacCdn.EQ.0 _d 0) THEN
                0181           upf(i,j,k)=recip_drF(k)*recip_hFacC(i,j,k,bi,bj)
                0182          ELSE
                0183           upf(i,j,k)=0 _d 0
                0184          ENDIF
                0185         ENDDO
                0186        ENDDO
                0187 #ifdef NONLIN_FRSURF
                0188 C-    Account for change in level thickness
                0189        IF (nonlinFreeSurf.GT.0) THEN
                0190          CALL FREESURF_RESCALE_G(
                0191      I                            bi, bj, k,
                0192      U                            upf,
                0193      I                            myThid )
                0194        ENDIF
                0195 #endif /* NONLIN_FRSURF */
                0196        DO j=jMin,jMax
                0197         DO i=iMin,iMax
                0198          upfc = upf(i,j,k)
                0199          IF (upfc .GT. 0 _d 0) THEN
086a45f245 Oliv*0200           flux = wPIC_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPIC))
dcd877c466 Oliv*0201           gTr(i,j,k  ,iPIC ) = gTr(i,j,k  ,iPIC ) - flux*upfc
248275f1c4 Oliv*0202 #ifdef DARWIN_ALLOW_CONS
                0203           botSnkC(i,j,bi,bj) = botSnkC(i,j,bi,bj) + dTsub(k)*flux
                0204 #endif
                0205 C
086a45f245 Oliv*0206           flux = wC_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOC))
dcd877c466 Oliv*0207           gTr(i,j,k  ,iPOC ) = gTr(i,j,k  ,iPOC ) - flux*upfc
248275f1c4 Oliv*0208 #ifdef DARWIN_ALLOW_CONS
                0209           botSnkC(i,j,bi,bj) = botSnkC(i,j,bi,bj) + dTsub(k)*flux
                0210 #endif
                0211 C
086a45f245 Oliv*0212           flux = wN_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPON))
dcd877c466 Oliv*0213           gTr(i,j,k  ,iPON ) = gTr(i,j,k  ,iPON ) - flux*upfc
248275f1c4 Oliv*0214 #ifdef DARWIN_ALLOW_CONS
                0215           botSnkN(i,j,bi,bj) = botSnkN(i,j,bi,bj) + dTsub(k)*flux
                0216 #endif
                0217 C
086a45f245 Oliv*0218           flux = wP_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOP))
dcd877c466 Oliv*0219           gTr(i,j,k  ,iPOP ) = gTr(i,j,k  ,iPOP ) - flux*upfc
248275f1c4 Oliv*0220 #ifdef DARWIN_ALLOW_CONS
                0221           botSnkP(i,j,bi,bj) = botSnkP(i,j,bi,bj) + dTsub(k)*flux
                0222 #endif
                0223 C
086a45f245 Oliv*0224           flux = wFe_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOFe))
dcd877c466 Oliv*0225           gTr(i,j,k  ,iPOFe) = gTr(i,j,k  ,iPOFe) - flux*upfc
248275f1c4 Oliv*0226 #ifdef DARWIN_ALLOW_CONS
                0227           botSnkFe(i,j,bi,bj) = botSnkFe(i,j,bi,bj) + dTsub(k)*flux
                0228 #endif
                0229 C
086a45f245 Oliv*0230           flux = wSi_sink*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iPOSi))
248275f1c4 Oliv*0231           gTr(i,j,k  ,iPOSi) = gTr(i,j,k  ,iPOSi) - flux*upfc
                0232 #ifdef DARWIN_ALLOW_CONS
                0233           botSnkSi(i,j,bi,bj) = botSnkSi(i,j,bi,bj) + dTsub(k)*flux
                0234 #endif
1fb46047a0 Oliv*0235           DO l = 1, nplank
086a45f245 Oliv*0236            flux = biosink(l)*MAX(0 _d 0, Ptr(i,j,k,bi,bj,ic+l-1))
dcd877c466 Oliv*0237            gTr(i,j,k  ,ic+l-1 )=gTr(i,j,k  ,ic+l-1 ) - flux*upfc
248275f1c4 Oliv*0238 #ifdef DARWIN_ALLOW_CONS
                0239            botSnkC(i,j,bi,bj) = botSnkC(i,j,bi,bj) + dTsub(k)*flux
                0240      &                                         *(1 _d 0 + R_PICPOC(l))
                0241 # ifndef DARWIN_ALLOW_NQUOTA
                0242            botSnkN(i,j,bi,bj) = botSnkN(i,j,bi,bj)
                0243      &                        + dTsub(k)*flux*R_NC(l)
                0244 # endif
                0245 # ifndef DARWIN_ALLOW_PQUOTA
                0246            botSnkP(i,j,bi,bj) = botSnkP(i,j,bi,bj)
                0247      &                        + dTsub(k)*flux*R_PC(l)
                0248 # endif
                0249 # ifndef DARWIN_ALLOW_FEQUOTA
                0250            botSnkFe(i,j,bi,bj) = botSnkFe(i,j,bi,bj)
                0251      &                         + dTsub(k)*flux*R_FeC(l)
                0252 # endif
                0253 # ifndef DARWIN_ALLOW_SIQUOTA
                0254            botSnkSi(i,j,bi,bj) = botSnkSi(i,j,bi,bj)
                0255      &                         + dTsub(k)*flux*R_SiC(l)
                0256 # endif
                0257 #endif
                0258 C
1fb46047a0 Oliv*0259 #ifdef DARWIN_ALLOW_NQUOTA
086a45f245 Oliv*0260            flux = biosink(l)*MAX(0 _d 0, Ptr(i,j,k,bi,bj,in+l-1))
dcd877c466 Oliv*0261            gTr(i,j,k  ,in+l-1 )=gTr(i,j,k  ,in+l-1 ) - flux*upfc
248275f1c4 Oliv*0262 # ifdef DARWIN_ALLOW_CONS
                0263            botSnkN(i,j,bi,bj) = botSnkN(i,j,bi,bj) + dTsub(k)*flux
                0264 # endif
1fb46047a0 Oliv*0265 #endif
248275f1c4 Oliv*0266 C
1fb46047a0 Oliv*0267 #ifdef DARWIN_ALLOW_PQUOTA
086a45f245 Oliv*0268            flux = biosink(l)*MAX(0 _d 0, Ptr(i,j,k,bi,bj,ip+l-1))
dcd877c466 Oliv*0269            gTr(i,j,k  ,ip+l-1 )=gTr(i,j,k  ,ip+l-1 ) - flux*upfc
248275f1c4 Oliv*0270 # ifdef DARWIN_ALLOW_CONS
                0271            botSnkP(i,j,bi,bj) = botSnkP(i,j,bi,bj) + dTsub(k)*flux
                0272 # endif
1fb46047a0 Oliv*0273 #endif
248275f1c4 Oliv*0274 C
1fb46047a0 Oliv*0275 #ifdef DARWIN_ALLOW_FEQUOTA
086a45f245 Oliv*0276            flux = biosink(l)*MAX(0 _d 0, Ptr(i,j,k,bi,bj,ife+l-1))
dcd877c466 Oliv*0277            gTr(i,j,k  ,ife+l-1)=gTr(i,j,k  ,ife+l-1) - flux*upfc
248275f1c4 Oliv*0278 # ifdef DARWIN_ALLOW_CONS
                0279            botSnkFe(i,j,bi,bj) = botSnkFe(i,j,bi,bj) + dTsub(k)*flux
                0280 # endif
                0281 #endif
                0282 C
                0283 #ifdef DARWIN_ALLOW_SIQUOTA
086a45f245 Oliv*0284            flux = biosink(l)*MAX(0 _d 0, Ptr(i,j,k,bi,bj,isi+l-1))
248275f1c4 Oliv*0285            gTr(i,j,k  ,isi+l-1)=gTr(i,j,k  ,isi+l-1) - flux*upfc
                0286 # ifdef DARWIN_ALLOW_CONS
                0287            botSnkSi(i,j,bi,bj) = botSnkSi(i,j,bi,bj) + dTsub(k)*flux
                0288 # endif
1fb46047a0 Oliv*0289 #endif
                0290           ENDDO
                0291 #ifdef DARWIN_ALLOW_CHLQUOTA
                0292           DO l = 1, nPhoto
086a45f245 Oliv*0293            flux = biosink(l)*MAX(0 _d 0, Ptr(i,j,k,bi,bj,iChl+l-1))
dcd877c466 Oliv*0294            gTr(i,j,k  ,iChl+l-1)=gTr(i,j,k  ,iChl+l-1)-flux*upfc
1fb46047a0 Oliv*0295           ENDDO
                0296 #endif
                0297          ENDIF
                0298         ENDDO
                0299        ENDDO
                0300       ENDDO
                0301 #endif /* DARWIN_BOTTOM_SINK */
                0302 
8fbfd1f382 Oliv*0303 #endif /* ALLOW_DARWIN */
                0304 
                0305       RETURN
                0306       END
                0307