model/src/integr_continuity.F

6d54cf9ca1 Ed H*0001 #include "PACKAGES_CONFIG.h"
fb5eaa30cd Jean*0002 #include "CPP_OPTIONS.h"
52dd931a85 Timo*0003 #ifdef ALLOW_AUTODIFF
                0004 # include "AUTODIFF_OPTIONS.h"
                0005 #endif
aecc8b0f47 Mart*0006 #ifdef ALLOW_CTRL
                0007 # include "CTRL_OPTIONS.h"
                0008 #endif
fb5eaa30cd Jean*0009
                0010 CBOP
                0011 C     !ROUTINE: INTEGR_CONTINUITY
                0012 C     !INTERFACE:
                0013       SUBROUTINE INTEGR_CONTINUITY(
bfdbc242ac Jean*0014      I                             uFld, vFld,
fb5eaa30cd Jean*0015      I                             myTime, myIter, myThid )
                0016 C     !DESCRIPTION: \bv
                0017 C     *==========================================================*
21b8df1d84 Jean*0018 C     | SUBROUTINE INTEGR_CONTINUITY
                0019 C     | o Integrate the continuity Eq : compute vertical velocity
bfdbc242ac Jean*0020 C     |   and free surface "r-anomaly" (etaN,etaH) to satisfy
                0021 C     |   exactly the conservation of the Total Volume
52dd931a85 Timo*0022 C     | o Routine dummy_for_etan located here in order to print the
                0023 C     |   correct adjoint variable for etaN
fb5eaa30cd Jean*0024 C     *==========================================================*
                0025 C     \ev
                0026
                0027 C     !USES:
                0028       IMPLICIT NONE
                0029 C     == Global variables
                0030 #include "SIZE.h"
                0031 #include "EEPARAMS.h"
                0032 #include "PARAMS.h"
                0033 #include "DYNVARS.h"
                0034 #include "GRID.h"
                0035 #include "SURFACE.h"
                0036 #include "FFIELDS.h"
aecc8b0f47 Mart*0037 #ifdef ALLOW_AUTODIFF_TAMC
                0038 # include "tamc.h"
                0039 #endif
fb5eaa30cd Jean*0040
                0041 C     !INPUT/OUTPUT PARAMETERS:
                0042 C     == Routine arguments ==
bfdbc242ac Jean*0043 C     uFld     :: Zonal velocity ( m/s )
                0044 C     vFld     :: Meridional velocity ( m/s )
                0045 C     myTime   :: Current time in simulation
                0046 C     myIter   :: Current iteration number in simulation
                0047 C     myThid   :: my Thread Id. number
fb5eaa30cd Jean*0048       _RL myTime
                0049       INTEGER myIter
                0050       INTEGER myThid
e1fb02e8f0 Jean*0051       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0052       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
fb5eaa30cd Jean*0053
                0054 C     !LOCAL VARIABLES:
                0055 C     Local variables in common block
                0056
                0057 C     Local variables
bfdbc242ac Jean*0058 C     bi,bj    :: tile index
                0059 C     i,j,k    :: Loop counters
                0060 C     uTrans   :: Volume transports ( uVel.xA )
                0061 C     vTrans   :: Volume transports ( vVel.yA )
8dc89a3cca Jean*0062 C     hDivFlow :: Div. Barotropic Flow [transport unit m3/s]
e1fb02e8f0 Jean*0063       INTEGER k,bi,bj
55e9ea8a90 Jean*0064 #ifdef EXACT_CONSERV
e1fb02e8f0 Jean*0065       INTEGER i, j
b924705ae1 Matt*0066       INTEGER ks
e1fb02e8f0 Jean*0067       _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL hDivFlow(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d2a11ab670 Jean*0070       _RL facEmP, facMass
e1fb02e8f0 Jean*0071 #else /* EXACT_CONSERV */
                0072 # ifdef ALLOW_OBCS
                0073       INTEGER i, j
                0074 # endif
d2a11ab670 Jean*0075 #endif /* EXACT_CONSERV */
                0076 #ifndef ALLOW_ADDFLUID
                0077       _RL addMass(1)
                0078 #endif /* ndef ALLOW_ADDFLUID */
65de7f3d8d Jean*0079 #if (defined NONLIN_FRSURF) && !(defined DISABLE_RSTAR_CODE)
e1fb02e8f0 Jean*0080       _RL rStarDhDt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
aecc8b0f47 Mart*0081 # if (defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_DEPTH_CONTROL)
edb6656069 Mart*0082 C     tkey :: tape key (tile dependent)
                0083       INTEGER tkey
aecc8b0f47 Mart*0084 # endif
65de7f3d8d Jean*0085 #else
                0086       _RL rStarDhDt(1)
                0087 #endif
fb5eaa30cd Jean*0088 CEOP
                0089
bfdbc242ac Jean*0090 C--   Start bi,bj loops
                0091       DO bj=myByLo(myThid),myByHi(myThid)
                0092        DO bi=myBxLo(myThid),myBxHi(myThid)
                0093
fb5eaa30cd Jean*0094 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0095
                0096 #ifdef EXACT_CONSERV
                0097       IF (exactConserv) THEN
                0098
d2a11ab670 Jean*0099        facEmP = 0.
                0100        IF ( fluidIsWater.AND.useRealFreshWaterFlux ) facEmP=mass2rUnit
                0101        facMass = 0.
                0102        IF ( selectAddFluid.GE.1 ) facMass = mass2rUnit
                0103
fb5eaa30cd Jean*0104 C--   Compute the Divergence of The Barotropic Flow :
                0105
21b8df1d84 Jean*0106 C-    Initialise
e1fb02e8f0 Jean*0107        DO j=1-OLy,sNy+OLy
                0108         DO i=1-OLx,sNx+OLx
8dc89a3cca Jean*0109          hDivFlow(i,j)      = 0. _d 0
                0110          utrans(i,j)        = 0. _d 0
                0111          vtrans(i,j)        = 0. _d 0
d2a11ab670 Jean*0112         ENDDO
fb5eaa30cd Jean*0113        ENDDO
                0114
d2a11ab670 Jean*0115        DO k=1,Nr
21b8df1d84 Jean*0116
fb5eaa30cd Jean*0117 C-    Calculate velocity field "volume transports" through tracer cell faces
4606c28752 Jean*0118 C     anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport).
fb5eaa30cd Jean*0119         DO j=1,sNy+1
                0120          DO i=1,sNx+1
                0121           uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
4606c28752 Jean*0122      &                 *deepFacC(k)*rhoFacC(k)
fb5eaa30cd Jean*0123      &                 *drF(k)*_hFacW(i,j,k,bi,bj)
                0124           vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
4606c28752 Jean*0125      &                 *deepFacC(k)*rhoFacC(k)
fb5eaa30cd Jean*0126      &                 *drF(k)*_hFacS(i,j,k,bi,bj)
                0127          ENDDO
                0128         ENDDO
                0129
21b8df1d84 Jean*0130 C-    Integrate vertically the Horizontal Divergence
fb5eaa30cd Jean*0131         DO j=1,sNy
                0132          DO i=1,sNx
8dc89a3cca Jean*0133            hDivFlow(i,j) = hDivFlow(i,j)
fb5eaa30cd Jean*0134      &       +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
                0135      &                            +vTrans(i,j+1)-vTrans(i,j) )
d2a11ab670 Jean*0136 #ifdef ALLOW_ADDFLUID
                0137      &       -facMass*addMass(i,j,k,bi,bj)
                0138 #endif /* ALLOW_ADDFLUID */
fb5eaa30cd Jean*0139          ENDDO
                0140         ENDDO
                0141
                0142 C-    End DO k=1,Nr
d2a11ab670 Jean*0143        ENDDO
fb5eaa30cd Jean*0144
8dc89a3cca Jean*0145 C------------------------------------
4606c28752 Jean*0146 C note: deep-model not implemented for P-coordinate + realFreshWaterFlux ;
                0147 C       anelastic: always assumes that rhoFacF(1) = 1
bfdbc242ac Jean*0148        IF ( myIter.EQ.nIter0 .AND. myIter.NE.0
                0149      &    .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
615c650f5e Jean*0150
8dc89a3cca Jean*0151 C     needs previous time-step value of E-P-R, that has not been loaded
615c650f5e Jean*0152 C     and was not in old pickup-file ; try to use etaN & etaH instead.
                0153          IF ( usePickupBeforeC54 ) THEN
8dc89a3cca Jean*0154           DO j=1,sNy
                0155            DO i=1,sNx
615c650f5e Jean*0156             dEtaHdt(i,j,bi,bj) = (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
edb6656069 Mart*0157      &                         / (implicDiv2Dflow*deltaTFreeSurf)
615c650f5e Jean*0158            ENDDO
                0159           ENDDO
                0160          ENDIF
                0161
                0162          DO j=1,sNy
                0163           DO i=1,sNx
8dc89a3cca Jean*0164             PmEpR(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)
                0165      &                       + hDivFlow(i,j)*recip_rA(i,j,bi,bj)
4606c28752 Jean*0166      &                                      *recip_deepFac2F(1)
62fd6ae4e5 Jean*0167             PmEpR(i,j,bi,bj) = PmEpR(i,j,bi,bj)*rUnit2mass
8dc89a3cca Jean*0168           ENDDO
615c650f5e Jean*0169          ENDDO
bfdbc242ac Jean*0170        ELSEIF ( myIter.EQ.nIter0 ) THEN
8dc89a3cca Jean*0171          DO j=1,sNy
                0172           DO i=1,sNx
edb6656069 Mart*0173             ks = kSurfC(i,j,bi,bj)
8dc89a3cca Jean*0174             PmEpR(i,j,bi,bj) = 0. _d 0
                0175             dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
4606c28752 Jean*0176      &                                         *recip_deepFac2F(ks)
8dc89a3cca Jean*0177           ENDDO
21b8df1d84 Jean*0178          ENDDO
8dc89a3cca Jean*0179        ELSE
21b8df1d84 Jean*0180 C--    Needs to get valid PmEpR (for T,S forcing) also in overlap regions
                0181 C       (e.g., if using KPP) => set over full index range
                0182          DO j=1-OLy,sNy+OLy
                0183           DO i=1-OLx,sNx+OLx
                0184             PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
                0185           ENDDO
                0186          ENDDO
8dc89a3cca Jean*0187          DO j=1,sNy
                0188           DO i=1,sNx
2d3bf0eeda Jean*0189             ks = kSurfC(i,j,bi,bj)
8dc89a3cca Jean*0190             dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
4606c28752 Jean*0191      &                                         *recip_deepFac2F(ks)
8dc89a3cca Jean*0192      &                           -facEmP*EmPmR(i,j,bi,bj)
                0193           ENDDO
                0194          ENDDO
                0195        ENDIF
                0196 C------------------------------------
                0197
2d3bf0eeda Jean*0198 #ifdef ALLOW_OBCS
                0199 C--    reset dEtaHdt to zero outside the OB interior region
                0200        IF ( useOBCS ) THEN
                0201          DO j=1,sNy
                0202           DO i=1,sNx
                0203             dEtaHdt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)*maskInC(i,j,bi,bj)
                0204           ENDDO
                0205          ENDDO
                0206        ENDIF
                0207 #endif /* ALLOW_OBCS */
                0208
                0209 C-- end if exactConserv block
fb5eaa30cd Jean*0210       ENDIF
                0211
                0212 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0213
bfdbc242ac Jean*0214       IF ( exactConserv .AND. myIter.NE.nIter0 ) THEN
21b8df1d84 Jean*0215 C--   Update etaN at the end of the time step :
fb5eaa30cd Jean*0216 C     Incorporate the Implicit part of -Divergence(Barotropic_Flow)
                0217
                0218        IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
e1fb02e8f0 Jean*0219         DO j=1-OLy,sNy+OLy
                0220          DO i=1-OLx,sNx+OLx
21b8df1d84 Jean*0221            etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
fb5eaa30cd Jean*0222          ENDDO
                0223         ENDDO
                0224        ELSE
                0225         DO j=1,sNy
                0226          DO i=1,sNx
                0227            etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
edb6656069 Mart*0228      &      + implicDiv2Dflow*dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
fb5eaa30cd Jean*0229          ENDDO
                0230         ENDDO
                0231        ENDIF
                0232
b3d6da8289 Dimi*0233 #ifdef ALLOW_OBCS
2d3bf0eeda Jean*0234 C--    Was added on Dec 30, 2004 (to fix unrealistic etaN ?), but no longer
                0235 C      needed with proper masking in solver (matrix+cg2d_b,x) and masking
                0236 C      of dEtaHdt above. etaN next to OB does not enter present algorithm but
                0237 C      leave it commented out in case we would implement an aternative scheme.
                0238 c      IF ( useOBCS ) CALL OBCS_APPLY_ETA( bi, bj, etaN, myThid )
21b8df1d84 Jean*0239 #endif /* ALLOW_OBCS */
604cf5604a Alis*0240
bfdbc242ac Jean*0241 C-- end if exactConserv and not myIter=nIter0 block
fb5eaa30cd Jean*0242       ENDIF
                0243
                0244 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0245
cdc9f269ae Patr*0246 # ifdef NONLIN_FRSURF
aecc8b0f47 Mart*0247 #  if (defined ALLOW_AUTODIFF_TAMC) && \
                0248       (defined ALLOW_DEPTH_CONTROL) && (!defined DISABLE_RSTAR_CODE)
edb6656069 Mart*0249       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0250 CADJ STORE dEtaHdt(:,:,bi,bj) = comlev1_bibj, key = tkey, kind = isbyte
aecc8b0f47 Mart*0251 #  endif
00b29feb62 Jean*0252       IF (select_rStar .NE. 0) THEN
cdc9f269ae Patr*0253 #  ifndef DISABLE_RSTAR_CODE
65de7f3d8d Jean*0254 C-- note: rStarDhDt is similar to rStarDhCDt from S/R CALC_R_STAR
                0255 C         except for deep-factor and rho factor.
8dc89a3cca Jean*0256         DO j=1,sNy
                0257          DO i=1,sNx
2d3bf0eeda Jean*0258            ks = kSurfC(i,j,bi,bj)
65de7f3d8d Jean*0259            rStarDhDt(i,j) = dEtaHdt(i,j,bi,bj)
                0260      &                     *deepFac2F(ks)*rhoFacF(ks)
                0261      &                     *recip_Rcol(i,j,bi,bj)
00b29feb62 Jean*0262          ENDDO
cdc9f269ae Patr*0263         ENDDO
                0264 #  endif /* DISABLE_RSTAR_CODE */
00b29feb62 Jean*0265       ENDIF
aecc8b0f47 Mart*0266 #  if (defined ALLOW_AUTODIFF_TAMC) && \
                0267       (defined ALLOW_DEPTH_CONTROL) && (!defined DISABLE_RSTAR_CODE)
edb6656069 Mart*0268 CADJ STORE rStarDhDt = comlev1_bibj, key = tkey, kind = isbyte
aecc8b0f47 Mart*0269 #  endif
cdc9f269ae Patr*0270 # endif /* NONLIN_FRSURF */
                0271
ef891516b3 Jean*0272 #endif /* EXACT_CONSERV */
00b29feb62 Jean*0273
                0274 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0275
fb5eaa30cd Jean*0276       DO k=Nr,1,-1
                0277 C--    Integrate continuity vertically for vertical velocity
                0278
                0279        CALL INTEGRATE_FOR_W(
65de7f3d8d Jean*0280      I                       bi, bj, k, uFld, vFld,
                0281      I                       addMass, rStarDhDt,
fb5eaa30cd Jean*0282      O                       wVel,
                0283      I                       myThid )
21b8df1d84 Jean*0284
f4db9965d8 Jean*0285 #ifdef EXACT_CONSERV
bfdbc242ac Jean*0286        IF ( k.EQ.Nr .AND. myIter.NE.0 .AND. usingPCoords
                0287      &      .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
e0c3eb6fa1 Jean*0288
2d3bf0eeda Jean*0289          DO j=1,sNy
                0290           DO i=1,sNx
                0291             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
62fd6ae4e5 Jean*0292      &               +mass2rUnit*PmEpR(i,j,bi,bj)*maskC(i,j,k,bi,bj)
e0c3eb6fa1 Jean*0293           ENDDO
2d3bf0eeda Jean*0294          ENDDO
e0c3eb6fa1 Jean*0295
                0296        ENDIF
f4db9965d8 Jean*0297 #endif /* EXACT_CONSERV */
e0c3eb6fa1 Jean*0298
fb5eaa30cd Jean*0299 #ifdef ALLOW_OBCS
2d3bf0eeda Jean*0300 C--    reset W to zero outside the OB interior region
                0301        IF ( useOBCS ) THEN
                0302          DO j=1,sNy
                0303           DO i=1,sNx
                0304              wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
                0305           ENDDO
                0306          ENDDO
                0307        ENDIF
                0308 C--    Apply OBC to W (non-hydrostatic):
                0309        IF ( useOBCS.AND.nonHydrostatic )
                0310      &                CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
21b8df1d84 Jean*0311 #endif /* ALLOW_OBCS */
fb5eaa30cd Jean*0312
                0313 C-    End DO k=Nr,1,-1
                0314       ENDDO
                0315
                0316 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0317
bfdbc242ac Jean*0318 C--   End bi,bj loops
                0319        ENDDO
                0320       ENDDO
                0321
52dd931a85 Timo*0322 #ifdef ALLOW_AUTODIFF_MONITOR
                0323 C     In reverse, print adjoint variable for etaN
                0324       CALL DUMMY_FOR_ETAN( myTime, myIter, myThid )
                0325       _EXCH_XY_RL( etaN , myThid )
                0326 #else /* ALLOW_AUTODIFF_MONITOR */
bfdbc242ac Jean*0327       IF ( exactConserv .AND. myIter.NE.nIter0
                0328      &                  .AND. implicDiv2Dflow .NE. 0. _d 0 )
c1e2d059e3 Jean*0329      &    _EXCH_XY_RL( etaN , myThid )
52dd931a85 Timo*0330 #endif /* ALLOW_AUTODIFF_MONITOR */
                0331
bfdbc242ac Jean*0332       IF ( implicitIntGravWave .OR. myIter.EQ.nIter0 )
c1e2d059e3 Jean*0333      &    _EXCH_XYZ_RL( wVel , myThid )
bfdbc242ac Jean*0334
                0335 #ifdef EXACT_CONSERV
                0336 C--   Update etaH (from etaN and dEtaHdt):
                0337       IF (exactConserv) THEN
c1e2d059e3 Jean*0338 c       IF ( useRealFreshWaterFlux .AND. myIter.EQ.nIter0 )
                0339 c    &    _EXCH_XY_RL( PmEpR, myThid )
bfdbc242ac Jean*0340 #ifdef ALLOW_DEBUG
                0341         IF (debugMode) CALL DEBUG_CALL('UPDATE_ETAH',myThid)
                0342 #endif
                0343         CALL UPDATE_ETAH( myTime, myIter, myThid )
                0344       ENDIF
c1e2d059e3 Jean*0345
                0346 #ifdef NONLIN_FRSURF
                0347 # ifndef DISABLE_SIGMA_CODE
                0348       IF ( nonlinFreeSurf.GT.0 .AND. selectSigmaCoord.NE.0 ) THEN
                0349         CALL EXCH_XY_RL( dEtaHdt, myThid )
                0350         CALL UPDATE_ETAWS( myTime, myIter, myThid )
                0351       ENDIF
                0352 # endif /* DISABLE_SIGMA_CODE */
                0353 #endif /* NONLIN_FRSURF */
                0354
bfdbc242ac Jean*0355 #endif /* EXACT_CONSERV */
                0356
fb5eaa30cd Jean*0357       RETURN
                0358       END