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
                0064       INTEGER i, j
b924705ae1 Matt*0065       INTEGER ks
e1fb02e8f0 Jean*0066       _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL hDivFlow(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d2a11ab670 Jean*0069       _RL facEmP, facMass
                0070 #ifndef ALLOW_ADDFLUID
                0071       _RL addMass(1)
                0072 #endif /* ndef ALLOW_ADDFLUID */
65de7f3d8d Jean*0073 #if (defined NONLIN_FRSURF) && !(defined DISABLE_RSTAR_CODE)
e1fb02e8f0 Jean*0074       _RL rStarDhDt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
aecc8b0f47 Mart*0075 # if (defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_DEPTH_CONTROL)
edb6656069 Mart*0076 C     tkey :: tape key (tile dependent)
                0077       INTEGER tkey
aecc8b0f47 Mart*0078 # endif
65de7f3d8d Jean*0079 #else
                0080       _RL rStarDhDt(1)
                0081 #endif
fb5eaa30cd Jean*0082 CEOP
                0083
bfdbc242ac Jean*0084 C--   Start bi,bj loops
                0085       DO bj=myByLo(myThid),myByHi(myThid)
                0086        DO bi=myBxLo(myThid),myBxHi(myThid)
                0087
fb5eaa30cd Jean*0088 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0089
79b5d5775c Jean*0090       IF ( exactConserv ) THEN
fb5eaa30cd Jean*0091
d2a11ab670 Jean*0092        facEmP = 0.
                0093        IF ( fluidIsWater.AND.useRealFreshWaterFlux ) facEmP=mass2rUnit
                0094        facMass = 0.
79b5d5775c Jean*0095        IF ( selectAddFluid.GE.1 .AND. myIter.NE.0 ) facMass = mass2rUnit
d2a11ab670 Jean*0096
fb5eaa30cd Jean*0097 C--   Compute the Divergence of The Barotropic Flow :
                0098
21b8df1d84 Jean*0099 C-    Initialise
e1fb02e8f0 Jean*0100        DO j=1-OLy,sNy+OLy
                0101         DO i=1-OLx,sNx+OLx
8dc89a3cca Jean*0102          hDivFlow(i,j)      = 0. _d 0
                0103          utrans(i,j)        = 0. _d 0
                0104          vtrans(i,j)        = 0. _d 0
d2a11ab670 Jean*0105         ENDDO
fb5eaa30cd Jean*0106        ENDDO
                0107
d2a11ab670 Jean*0108        DO k=1,Nr
21b8df1d84 Jean*0109
fb5eaa30cd Jean*0110 C-    Calculate velocity field "volume transports" through tracer cell faces
4606c28752 Jean*0111 C     anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport).
fb5eaa30cd Jean*0112         DO j=1,sNy+1
                0113          DO i=1,sNx+1
                0114           uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
4606c28752 Jean*0115      &                 *deepFacC(k)*rhoFacC(k)
fb5eaa30cd Jean*0116      &                 *drF(k)*_hFacW(i,j,k,bi,bj)
                0117           vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
4606c28752 Jean*0118      &                 *deepFacC(k)*rhoFacC(k)
fb5eaa30cd Jean*0119      &                 *drF(k)*_hFacS(i,j,k,bi,bj)
                0120          ENDDO
                0121         ENDDO
                0122
21b8df1d84 Jean*0123 C-    Integrate vertically the Horizontal Divergence
fb5eaa30cd Jean*0124         DO j=1,sNy
                0125          DO i=1,sNx
8dc89a3cca Jean*0126            hDivFlow(i,j) = hDivFlow(i,j)
fb5eaa30cd Jean*0127      &       +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
                0128      &                            +vTrans(i,j+1)-vTrans(i,j) )
d2a11ab670 Jean*0129 #ifdef ALLOW_ADDFLUID
                0130      &       -facMass*addMass(i,j,k,bi,bj)
                0131 #endif /* ALLOW_ADDFLUID */
fb5eaa30cd Jean*0132          ENDDO
                0133         ENDDO
                0134
                0135 C-    End DO k=1,Nr
d2a11ab670 Jean*0136        ENDDO
fb5eaa30cd Jean*0137
8dc89a3cca Jean*0138 C------------------------------------
4606c28752 Jean*0139 C note: deep-model not implemented for P-coordinate + realFreshWaterFlux ;
                0140 C       anelastic: always assumes that rhoFacF(1) = 1
bfdbc242ac Jean*0141        IF ( myIter.EQ.nIter0 .AND. myIter.NE.0
                0142      &    .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
615c650f5e Jean*0143
8dc89a3cca Jean*0144 C     needs previous time-step value of E-P-R, that has not been loaded
615c650f5e Jean*0145 C     and was not in old pickup-file ; try to use etaN & etaH instead.
                0146          IF ( usePickupBeforeC54 ) THEN
8dc89a3cca Jean*0147           DO j=1,sNy
                0148            DO i=1,sNx
615c650f5e Jean*0149             dEtaHdt(i,j,bi,bj) = (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
edb6656069 Mart*0150      &                         / (implicDiv2Dflow*deltaTFreeSurf)
615c650f5e Jean*0151            ENDDO
                0152           ENDDO
                0153          ENDIF
                0154
                0155          DO j=1,sNy
                0156           DO i=1,sNx
8dc89a3cca Jean*0157             PmEpR(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)
                0158      &                       + hDivFlow(i,j)*recip_rA(i,j,bi,bj)
4606c28752 Jean*0159      &                                      *recip_deepFac2F(1)
62fd6ae4e5 Jean*0160             PmEpR(i,j,bi,bj) = PmEpR(i,j,bi,bj)*rUnit2mass
8dc89a3cca Jean*0161           ENDDO
615c650f5e Jean*0162          ENDDO
bfdbc242ac Jean*0163        ELSEIF ( myIter.EQ.nIter0 ) THEN
8dc89a3cca Jean*0164          DO j=1,sNy
                0165           DO i=1,sNx
edb6656069 Mart*0166             ks = kSurfC(i,j,bi,bj)
8dc89a3cca Jean*0167             PmEpR(i,j,bi,bj) = 0. _d 0
                0168             dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
4606c28752 Jean*0169      &                                         *recip_deepFac2F(ks)
8dc89a3cca Jean*0170           ENDDO
21b8df1d84 Jean*0171          ENDDO
8dc89a3cca Jean*0172        ELSE
21b8df1d84 Jean*0173 C--    Needs to get valid PmEpR (for T,S forcing) also in overlap regions
                0174 C       (e.g., if using KPP) => set over full index range
                0175          DO j=1-OLy,sNy+OLy
                0176           DO i=1-OLx,sNx+OLx
                0177             PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
                0178           ENDDO
                0179          ENDDO
8dc89a3cca Jean*0180          DO j=1,sNy
                0181           DO i=1,sNx
2d3bf0eeda Jean*0182             ks = kSurfC(i,j,bi,bj)
8dc89a3cca Jean*0183             dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
4606c28752 Jean*0184      &                                         *recip_deepFac2F(ks)
8dc89a3cca Jean*0185      &                           -facEmP*EmPmR(i,j,bi,bj)
                0186           ENDDO
                0187          ENDDO
                0188        ENDIF
                0189 C------------------------------------
                0190
2d3bf0eeda Jean*0191 #ifdef ALLOW_OBCS
                0192 C--    reset dEtaHdt to zero outside the OB interior region
                0193        IF ( useOBCS ) THEN
                0194          DO j=1,sNy
                0195           DO i=1,sNx
                0196             dEtaHdt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)*maskInC(i,j,bi,bj)
                0197           ENDDO
                0198          ENDDO
                0199        ENDIF
                0200 #endif /* ALLOW_OBCS */
                0201
                0202 C-- end if exactConserv block
fb5eaa30cd Jean*0203       ENDIF
                0204
                0205 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0206
bfdbc242ac Jean*0207       IF ( exactConserv .AND. myIter.NE.nIter0 ) THEN
21b8df1d84 Jean*0208 C--   Update etaN at the end of the time step :
fb5eaa30cd Jean*0209 C     Incorporate the Implicit part of -Divergence(Barotropic_Flow)
                0210
79b5d5775c Jean*0211        IF ( implicDiv2Dflow.EQ.zeroRL ) THEN
e1fb02e8f0 Jean*0212         DO j=1-OLy,sNy+OLy
                0213          DO i=1-OLx,sNx+OLx
21b8df1d84 Jean*0214            etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
fb5eaa30cd Jean*0215          ENDDO
                0216         ENDDO
                0217        ELSE
                0218         DO j=1,sNy
                0219          DO i=1,sNx
                0220            etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
edb6656069 Mart*0221      &      + implicDiv2Dflow*dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
fb5eaa30cd Jean*0222          ENDDO
                0223         ENDDO
                0224        ENDIF
                0225
b3d6da8289 Dimi*0226 #ifdef ALLOW_OBCS
2d3bf0eeda Jean*0227 C--    Was added on Dec 30, 2004 (to fix unrealistic etaN ?), but no longer
                0228 C      needed with proper masking in solver (matrix+cg2d_b,x) and masking
                0229 C      of dEtaHdt above. etaN next to OB does not enter present algorithm but
                0230 C      leave it commented out in case we would implement an aternative scheme.
                0231 c      IF ( useOBCS ) CALL OBCS_APPLY_ETA( bi, bj, etaN, myThid )
21b8df1d84 Jean*0232 #endif /* ALLOW_OBCS */
604cf5604a Alis*0233
bfdbc242ac Jean*0234 C-- end if exactConserv and not myIter=nIter0 block
fb5eaa30cd Jean*0235       ENDIF
                0236
                0237 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0238
cdc9f269ae Patr*0239 # ifdef NONLIN_FRSURF
aecc8b0f47 Mart*0240 #  if (defined ALLOW_AUTODIFF_TAMC) && \
                0241       (defined ALLOW_DEPTH_CONTROL) && (!defined DISABLE_RSTAR_CODE)
edb6656069 Mart*0242       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0243 CADJ STORE dEtaHdt(:,:,bi,bj) = comlev1_bibj, key = tkey, kind = isbyte
aecc8b0f47 Mart*0244 #  endif
79b5d5775c Jean*0245       IF ( select_rStar.NE.0 ) THEN
cdc9f269ae Patr*0246 #  ifndef DISABLE_RSTAR_CODE
65de7f3d8d Jean*0247 C-- note: rStarDhDt is similar to rStarDhCDt from S/R CALC_R_STAR
                0248 C         except for deep-factor and rho factor.
8dc89a3cca Jean*0249         DO j=1,sNy
                0250          DO i=1,sNx
2d3bf0eeda Jean*0251            ks = kSurfC(i,j,bi,bj)
65de7f3d8d Jean*0252            rStarDhDt(i,j) = dEtaHdt(i,j,bi,bj)
                0253      &                     *deepFac2F(ks)*rhoFacF(ks)
                0254      &                     *recip_Rcol(i,j,bi,bj)
00b29feb62 Jean*0255          ENDDO
cdc9f269ae Patr*0256         ENDDO
                0257 #  endif /* DISABLE_RSTAR_CODE */
00b29feb62 Jean*0258       ENDIF
aecc8b0f47 Mart*0259 #  if (defined ALLOW_AUTODIFF_TAMC) && \
                0260       (defined ALLOW_DEPTH_CONTROL) && (!defined DISABLE_RSTAR_CODE)
edb6656069 Mart*0261 CADJ STORE rStarDhDt = comlev1_bibj, key = tkey, kind = isbyte
aecc8b0f47 Mart*0262 #  endif
cdc9f269ae Patr*0263 # endif /* NONLIN_FRSURF */
                0264
00b29feb62 Jean*0265 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0266
fb5eaa30cd Jean*0267       DO k=Nr,1,-1
                0268 C--    Integrate continuity vertically for vertical velocity
                0269
                0270        CALL INTEGRATE_FOR_W(
65de7f3d8d Jean*0271      I                       bi, bj, k, uFld, vFld,
                0272      I                       addMass, rStarDhDt,
fb5eaa30cd Jean*0273      O                       wVel,
79b5d5775c Jean*0274      I                       myIter, myThid )
21b8df1d84 Jean*0275
bfdbc242ac Jean*0276        IF ( k.EQ.Nr .AND. myIter.NE.0 .AND. usingPCoords
                0277      &      .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
e0c3eb6fa1 Jean*0278
2d3bf0eeda Jean*0279          DO j=1,sNy
                0280           DO i=1,sNx
                0281             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
62fd6ae4e5 Jean*0282      &               +mass2rUnit*PmEpR(i,j,bi,bj)*maskC(i,j,k,bi,bj)
e0c3eb6fa1 Jean*0283           ENDDO
2d3bf0eeda Jean*0284          ENDDO
e0c3eb6fa1 Jean*0285
                0286        ENDIF
                0287
fb5eaa30cd Jean*0288 #ifdef ALLOW_OBCS
2d3bf0eeda Jean*0289 C--    reset W to zero outside the OB interior region
                0290        IF ( useOBCS ) THEN
                0291          DO j=1,sNy
                0292           DO i=1,sNx
                0293              wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
                0294           ENDDO
                0295          ENDDO
                0296        ENDIF
                0297 C--    Apply OBC to W (non-hydrostatic):
                0298        IF ( useOBCS.AND.nonHydrostatic )
                0299      &                CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
21b8df1d84 Jean*0300 #endif /* ALLOW_OBCS */
fb5eaa30cd Jean*0301
                0302 C-    End DO k=Nr,1,-1
                0303       ENDDO
                0304
                0305 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0306
bfdbc242ac Jean*0307 C--   End bi,bj loops
                0308        ENDDO
                0309       ENDDO
                0310
52dd931a85 Timo*0311 #ifdef ALLOW_AUTODIFF_MONITOR
                0312 C     In reverse, print adjoint variable for etaN
                0313       CALL DUMMY_FOR_ETAN( myTime, myIter, myThid )
                0314       _EXCH_XY_RL( etaN , myThid )
                0315 #else /* ALLOW_AUTODIFF_MONITOR */
bfdbc242ac Jean*0316       IF ( exactConserv .AND. myIter.NE.nIter0
                0317      &                  .AND. implicDiv2Dflow .NE. 0. _d 0 )
c1e2d059e3 Jean*0318      &    _EXCH_XY_RL( etaN , myThid )
52dd931a85 Timo*0319 #endif /* ALLOW_AUTODIFF_MONITOR */
                0320
bfdbc242ac Jean*0321       IF ( implicitIntGravWave .OR. myIter.EQ.nIter0 )
c1e2d059e3 Jean*0322      &    _EXCH_XYZ_RL( wVel , myThid )
bfdbc242ac Jean*0323
                0324 C--   Update etaH (from etaN and dEtaHdt):
79b5d5775c Jean*0325       IF ( exactConserv ) THEN
c1e2d059e3 Jean*0326 c       IF ( useRealFreshWaterFlux .AND. myIter.EQ.nIter0 )
                0327 c    &    _EXCH_XY_RL( PmEpR, myThid )
bfdbc242ac Jean*0328 #ifdef ALLOW_DEBUG
79b5d5775c Jean*0329         IF ( debugMode ) CALL DEBUG_CALL('UPDATE_ETAH',myThid)
bfdbc242ac Jean*0330 #endif
                0331         CALL UPDATE_ETAH( myTime, myIter, myThid )
                0332       ENDIF
c1e2d059e3 Jean*0333
                0334 #ifdef NONLIN_FRSURF
                0335 # ifndef DISABLE_SIGMA_CODE
                0336       IF ( nonlinFreeSurf.GT.0 .AND. selectSigmaCoord.NE.0 ) THEN
                0337         CALL EXCH_XY_RL( dEtaHdt, myThid )
                0338         CALL UPDATE_ETAWS( myTime, myIter, myThid )
                0339       ENDIF
                0340 # endif /* DISABLE_SIGMA_CODE */
                0341 #endif /* NONLIN_FRSURF */
                0342
fb5eaa30cd Jean*0343       RETURN
                0344       END