** Warning **

Issuing rollback() due to DESTROY without explicit disconnect() of DBD::mysql::db handle dbname=MITgcm at /usr/local/share/lxr/lib/LXR/Common.pm line 1224.

Last-Modified: Thu, 15 May 2024 05:11:27 GMT Content-Type: text/html; charset=utf-8 MITgcm/MITgcm/model/src/timestep.F
Back to home page

MITgcm

 
 

    


File indexing completed on 2019-08-15 05:10:18 UTC

view on githubraw file Latest commit df999eca on 2019-03-20 18:51:56 UTC
138482fdf6 Ed H*0001 #include "PACKAGES_CONFIG.h"
73470161ee Alis*0002 #include "CPP_OPTIONS.h"
da4112fbea Jean*0003 #ifdef ALLOW_CD_CODE
                0004 #include "CD_CODE_OPTIONS.h"
                0005 #endif
924557e60a Chri*0006 
9366854e02 Chri*0007 CBOP
                0008 C     !ROUTINE: TIMESTEP
                0009 C     !INTERFACE:
c0911a93b9 Jean*0010       SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
ecf17841ce Jean*0011      I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
4667e4066d Jean*0012      I                     guDissip, gvDissip,
c0911a93b9 Jean*0013      I                     myTime, myIter, myThid )
9366854e02 Chri*0014 C     !DESCRIPTION: \bv
                0015 C     *==========================================================*
da4112fbea Jean*0016 C     | S/R TIMESTEP
                0017 C     | o Step model fields forward in time
9366854e02 Chri*0018 C     *==========================================================*
                0019 C     \ev
262eb6133c Patr*0020 
9366854e02 Chri*0021 C     !USES:
                0022       IMPLICIT NONE
262eb6133c Patr*0023 C     == Global variables ==
924557e60a Chri*0024 #include "SIZE.h"
81bc00c2f0 Chri*0025 #include "EEPARAMS.h"
924557e60a Chri*0026 #include "PARAMS.h"
                0027 #include "GRID.h"
dce55d02b7 Jean*0028 #include "SURFACE.h"
16f5093311 Jean*0029 #include "RESTART.h"
                0030 #include "DYNVARS.h"
                0031 #ifdef ALLOW_NONHYDROSTATIC
                0032 #include "NH_VARS.h"
                0033 #endif
262eb6133c Patr*0034 
9366854e02 Chri*0035 C     !INPUT/OUTPUT PARAMETERS:
924557e60a Chri*0036 C     == Routine Arguments ==
7e1bd93d4f Jean*0037 C     dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
                0038 C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
                0039 C     phiSurfY ::          or geopotential (atmos) in X and Y direction
4667e4066d Jean*0040 C     guDissip :: dissipation tendency (all explicit terms), u component
                0041 C     gvDissip :: dissipation tendency (all explicit terms), v component
                0042 
924557e60a Chri*0043       INTEGER bi,bj,iMin,iMax,jMin,jMax
c0911a93b9 Jean*0044       INTEGER k
94e3f14b29 Jean*0045       _RL     dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL     dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
678051767f Jean*0047       _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0048       _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4667e4066d Jean*0049       _RL     guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0050       _RL     gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c0911a93b9 Jean*0051       _RL     myTime
f1ed2dfe29 Alis*0052       INTEGER myIter, myThid
262eb6133c Patr*0053 
9366854e02 Chri*0054 C     !LOCAL VARIABLES:
924557e60a Chri*0055 C     == Local variables ==
f6687d0ce7 Jean*0056 C     guExt    :: forcing tendency, u component
                0057 C     gvExt    :: forcing tendency, v component
                0058 C     gu_AB    :: tendency increment from Adams-Bashforth, u component
                0059 C     gv_AB    :: tendency increment from Adams-Bashforth, v component
7de45374e9 Alis*0060       INTEGER i,j
df999eca2c Jean*0061       _RL phFac, psFac
f6687d0ce7 Jean*0062       _RL     guExt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0063       _RL     gvExt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
dce55d02b7 Jean*0064       _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0065       _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b88f9ec5cb Jean*0066       _RL     gu_AB(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067       _RL     gv_AB(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
df999eca2c Jean*0068       _RL     gUdPx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL     gVdPy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
16f5093311 Jean*0070 #ifdef ALLOW_NONHYDROSTATIC
                0071       _RL     nhFac
                0072 #endif
138482fdf6 Ed H*0073 #ifdef ALLOW_CD_CODE
c0911a93b9 Jean*0074       _RL     guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075       _RL     gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0076 #endif
9366854e02 Chri*0077 CEOP
924557e60a Chri*0078 
7e1bd93d4f Jean*0079 C-- explicit part of the surface potential gradient is added in this S/R
                0080       psFac = pfFacMom*(1. _d 0 - implicSurfPress)
16f5093311 Jean*0081      &       *recip_deepFacC(k)*recip_rhoFacC(k)
7e1bd93d4f Jean*0082 
722a76e285 Jean*0083 C--  factors for gradient (X & Y directions) of Hydrostatic Potential
df999eca2c Jean*0084       phFac = pfFacMom
722a76e285 Jean*0085 
dce55d02b7 Jean*0086 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
c0911a93b9 Jean*0087 
f6687d0ce7 Jean*0088 C-    Initialize local arrays (not really necessary for all but safer)
94e3f14b29 Jean*0089       DO j=1-OLy,sNy+OLy
                0090         DO i=1-OLx,sNx+OLx
f6687d0ce7 Jean*0091           guExt(i,j) = 0. _d 0
                0092           gvExt(i,j) = 0. _d 0
da4112fbea Jean*0093           gUtmp(i,j) = 0. _d 0
                0094           gVtmp(i,j) = 0. _d 0
df999eca2c Jean*0095           gUdPx(i,j) = 0. _d 0
                0096           gVdPy(i,j) = 0. _d 0
138482fdf6 Ed H*0097 #ifdef ALLOW_CD_CODE
da4112fbea Jean*0098           guCor(i,j) = 0. _d 0
                0099           gvCor(i,j) = 0. _d 0
c0911a93b9 Jean*0100 #endif
da4112fbea Jean*0101         ENDDO
c0911a93b9 Jean*0102       ENDDO
                0103 
f6687d0ce7 Jean*0104       IF ( momForcing ) THEN
                0105 C--   Collect forcing term in local array guExt,gvExt:
                0106         CALL APPLY_FORCING_U(
                0107      U                        guExt,
                0108      I                        iMin,iMax,jMin,jMax, k, bi,bj,
                0109      I                        myTime, myIter, myThid )
                0110         CALL APPLY_FORCING_V(
                0111      U                        gvExt,
                0112      I                        iMin,iMax,jMin,jMax, k, bi,bj,
                0113      I                        myTime, myIter, myThid )
                0114       ENDIF
                0115 
df999eca2c Jean*0116       IF ( .NOT.staggerTimeStep .AND. .NOT.implicitIntGravWave ) THEN
4667e4066d Jean*0117 C--   Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
722a76e285 Jean*0118         DO j=jMin,jMax
                0119          DO i=iMin,iMax
df999eca2c Jean*0120            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phFac*dPhiHydX(i,j)
                0121            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phFac*dPhiHydY(i,j)
722a76e285 Jean*0122          ENDDO
                0123         ENDDO
                0124 c     ELSE
                0125 C--   Stagger time step: grad Phi_Hyp will be added later
da4112fbea Jean*0126       ENDIF
4667e4066d Jean*0127 
                0128 C--   Dissipation term inside the Adams-Bashforth:
722a76e285 Jean*0129       IF ( momViscosity .AND. momDissip_In_AB) THEN
4667e4066d Jean*0130         DO j=jMin,jMax
                0131          DO i=iMin,iMax
da4112fbea Jean*0132            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
                0133            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
4667e4066d Jean*0134          ENDDO
                0135         ENDDO
                0136       ENDIF
                0137 
c0911a93b9 Jean*0138 C--   Forcing term inside the Adams-Bashforth:
c2b6ed6bfd Jean*0139       IF ( momForcing .AND. momForcingOutAB.NE.1 ) THEN
f6687d0ce7 Jean*0140         DO j=jMin,jMax
                0141          DO i=iMin,iMax
                0142            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guExt(i,j)
                0143            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvExt(i,j)
ca16750d41 Bayl*0144          ENDDO
f6687d0ce7 Jean*0145         ENDDO
c0911a93b9 Jean*0146       ENDIF
                0147 
da4112fbea Jean*0148 #ifdef CD_CODE_NO_AB_MOMENTUM
f6687d0ce7 Jean*0149       IF ( useCDscheme ) THEN
da4112fbea Jean*0150 C-    CD-scheme, before doing AB, store gU,Vtmp = gU,V^n (+dissip. +forcing)
                0151         DO j=jMin,jMax
                0152          DO i=iMin,iMax
                0153            gUtmp(i,j) = gU(i,j,k,bi,bj)
                0154            gVtmp(i,j) = gV(i,j,k,bi,bj)
                0155          ENDDO
                0156         ENDDO
aeffad48e9 Jean*0157       ENDIF
da4112fbea Jean*0158 #endif /* CD_CODE_NO_AB_MOMENTUM */
aeffad48e9 Jean*0159 
c0911a93b9 Jean*0160 C-    Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
                0161 C     and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
aeffad48e9 Jean*0162 #ifdef ALLOW_ADAMSBASHFORTH_3
                0163       CALL ADAMS_BASHFORTH3(
94e3f14b29 Jean*0164      I                        bi, bj, k, Nr,
0b6eaf7dba Jean*0165      U                        gU(1-OLx,1-OLy,1,bi,bj), guNm,
                0166      O                        gu_AB,
117ee419f5 Jean*0167      I                        mom_StartAB, myIter, myThid )
aeffad48e9 Jean*0168       CALL ADAMS_BASHFORTH3(
94e3f14b29 Jean*0169      I                        bi, bj, k, Nr,
0b6eaf7dba Jean*0170      U                        gV(1-OLx,1-OLy,1,bi,bj), gvNm,
                0171      O                        gv_AB,
117ee419f5 Jean*0172      I                        mom_StartAB, myIter, myThid )
aeffad48e9 Jean*0173 #else /* ALLOW_ADAMSBASHFORTH_3 */
                0174       CALL ADAMS_BASHFORTH2(
94e3f14b29 Jean*0175      I                        bi, bj, k, Nr,
0b6eaf7dba Jean*0176      U                        gU(1-OLx,1-OLy,1,bi,bj),
                0177      U                        guNm1(1-OLx,1-OLy,1,bi,bj),
                0178      O                        gu_AB,
117ee419f5 Jean*0179      I                        mom_StartAB, myIter, myThid )
aeffad48e9 Jean*0180       CALL ADAMS_BASHFORTH2(
94e3f14b29 Jean*0181      I                        bi, bj, k, Nr,
0b6eaf7dba Jean*0182      U                        gV(1-OLx,1-OLy,1,bi,bj),
                0183      U                        gvNm1(1-OLx,1-OLy,1,bi,bj),
                0184      O                        gv_AB,
117ee419f5 Jean*0185      I                        mom_StartAB, myIter, myThid )
aeffad48e9 Jean*0186 #endif /* ALLOW_ADAMSBASHFORTH_3 */
b88f9ec5cb Jean*0187 #ifdef ALLOW_DIAGNOSTICS
                0188       IF ( useDiagnostics ) THEN
b6ca249bb1 Jean*0189         CALL DIAGNOSTICS_FILL(gu_AB,'AB_gU   ',k,1,2,bi,bj,myThid)
                0190         CALL DIAGNOSTICS_FILL(gv_AB,'AB_gV   ',k,1,2,bi,bj,myThid)
b88f9ec5cb Jean*0191       ENDIF
                0192 #endif /* ALLOW_DIAGNOSTICS */
da4112fbea Jean*0193 
f6687d0ce7 Jean*0194 C-     Make a local copy in gU,Vtmp of gU,V^n+1/2 (+dissip. +forcing)
                0195 #ifdef CD_CODE_NO_AB_MOMENTUM
                0196       IF ( .NOT.useCDscheme ) THEN
                0197 #endif
                0198         DO j=jMin,jMax
                0199          DO i=iMin,iMax
                0200            gUtmp(i,j) = gU(i,j,k,bi,bj)
                0201            gVtmp(i,j) = gV(i,j,k,bi,bj)
                0202          ENDDO
                0203         ENDDO
                0204 #ifdef CD_CODE_NO_AB_MOMENTUM
                0205       ENDIF
                0206 #endif
                0207 
c0911a93b9 Jean*0208 C--   Forcing term outside the Adams-Bashforth:
c2b6ed6bfd Jean*0209       IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
f6687d0ce7 Jean*0210         DO j=jMin,jMax
                0211          DO i=iMin,iMax
                0212            gUtmp(i,j) = gUtmp(i,j) + guExt(i,j)
                0213            gVtmp(i,j) = gVtmp(i,j) + gvExt(i,j)
ca16750d41 Bayl*0214          ENDDO
f6687d0ce7 Jean*0215         ENDDO
                0216       ENDIF
aeffad48e9 Jean*0217 
f6687d0ce7 Jean*0218 C--   Dissipation term outside the Adams-Bashforth:
                0219       IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
                0220         DO j=jMin,jMax
                0221          DO i=iMin,iMax
                0222            gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
                0223            gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
ca16750d41 Bayl*0224          ENDDO
f6687d0ce7 Jean*0225         ENDDO
dce55d02b7 Jean*0226       ENDIF
                0227 
138482fdf6 Ed H*0228 #ifdef ALLOW_CD_CODE
f6687d0ce7 Jean*0229       IF ( useCDscheme ) THEN
                0230 C-     Step forward D-grid velocity using C-grid gU,Vtmp = gU,V +dissip +forcing
                0231 C      and return coriolis terms on C-grid (guCor,gvCor)
                0232         CALL CD_CODE_SCHEME(
                0233      I                  bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
                0234      O                  guCor,gvCor,
                0235      I                  myTime, myIter, myThid)
                0236 
da4112fbea Jean*0237 #ifdef CD_CODE_NO_AB_MOMENTUM
                0238         IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
                0239          DO j=jMin,jMax
                0240           DO i=iMin,iMax
f6687d0ce7 Jean*0241            gUtmp(i,j) = ( gU(i,j,k,bi,bj) + guExt(i,j) ) + guCor(i,j)
                0242            gVtmp(i,j) = ( gV(i,j,k,bi,bj) + gvExt(i,j) ) + gvCor(i,j)
da4112fbea Jean*0243           ENDDO
                0244          ENDDO
f6687d0ce7 Jean*0245         ELSE
                0246          DO j=jMin,jMax
                0247           DO i=iMin,iMax
                0248            gUtmp(i,j) = gU(i,j,k,bi,bj) + guCor(i,j)
                0249            gVtmp(i,j) = gV(i,j,k,bi,bj) + gvCor(i,j)
                0250           ENDDO
da4112fbea Jean*0251          ENDDO
f6687d0ce7 Jean*0252         ENDIF
da4112fbea Jean*0253         IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
                0254          DO j=jMin,jMax
                0255           DO i=iMin,iMax
                0256            gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
                0257            gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
                0258           ENDDO
                0259          ENDDO
                0260         ENDIF
f6687d0ce7 Jean*0261 #else /* CD_CODE_NO_AB_MOMENTUM */
aeffad48e9 Jean*0262         DO j=jMin,jMax
                0263          DO i=iMin,iMax
f6687d0ce7 Jean*0264            gUtmp(i,j) = gUtmp(i,j) + guCor(i,j)
                0265            gVtmp(i,j) = gVtmp(i,j) + gvCor(i,j)
aeffad48e9 Jean*0266          ENDDO
                0267         ENDDO
f6687d0ce7 Jean*0268 #endif /* CD_CODE_NO_AB_MOMENTUM */
aeffad48e9 Jean*0269       ENDIF
da4112fbea Jean*0270 #endif /* ALLOW_CD_CODE */
fb481a83c2 Alis*0271 
7dcc33f8da Jean*0272 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0273 
dce55d02b7 Jean*0274 #ifdef NONLIN_FRSURF
da4112fbea Jean*0275       IF ( .NOT. vectorInvariantMomentum
                0276      &     .AND. nonlinFreeSurf.GT.1 ) THEN
a2a20dcddc Jean*0277        IF ( select_rStar.GT.0 ) THEN
                0278 # ifndef DISABLE_RSTAR_CODE
00b29feb62 Jean*0279         DO j=jMin,jMax
                0280          DO i=iMin,iMax
c0911a93b9 Jean*0281            gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
00b29feb62 Jean*0282            gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
                0283          ENDDO
                0284         ENDDO
a2a20dcddc Jean*0285 # endif /* DISABLE_RSTAR_CODE */
                0286        ELSEIF ( selectSigmaCoord.NE.0 ) THEN
                0287 # ifndef DISABLE_SIGMA_CODE
                0288         DO j=jMin,jMax
                0289          DO i=iMin,iMax
                0290            gUtmp(i,j) = gUtmp(i,j)
94e3f14b29 Jean*0291      &        /( 1. _d 0 + dEtaWdt(i,j,bi,bj)*deltaTFreeSurf
a2a20dcddc Jean*0292      &                    *dBHybSigF(k)*recip_drF(k)
                0293      &                    *recip_hFacW(i,j,k,bi,bj)
                0294      &         )
                0295            gVtmp(i,j) = gVtmp(i,j)
94e3f14b29 Jean*0296      &        /( 1. _d 0 + dEtaSdt(i,j,bi,bj)*deltaTFreeSurf
a2a20dcddc Jean*0297      &                    *dBHybSigF(k)*recip_drF(k)
                0298      &                    *recip_hFacS(i,j,k,bi,bj)
                0299      &         )
                0300          ENDDO
                0301         ENDDO
                0302 # endif /* DISABLE_SIGMA_CODE */
00b29feb62 Jean*0303        ELSE
                0304         DO j=jMin,jMax
                0305          DO i=iMin,iMax
a2a20dcddc Jean*0306           IF ( k.EQ.kSurfW(i,j,bi,bj) ) THEN
c0911a93b9 Jean*0307            gUtmp(i,j) = gUtmp(i,j)
616600b8d2 Patr*0308      &         *_hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
c0911a93b9 Jean*0309           ENDIF
a2a20dcddc Jean*0310           IF ( k.EQ.kSurfS(i,j,bi,bj) ) THEN
dce55d02b7 Jean*0311            gVtmp(i,j) = gVtmp(i,j)
616600b8d2 Patr*0312      &         *_hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
00b29feb62 Jean*0313           ENDIF
                0314          ENDDO
dce55d02b7 Jean*0315         ENDDO
00b29feb62 Jean*0316        ENDIF
dce55d02b7 Jean*0317       ENDIF
c0911a93b9 Jean*0318 #endif /* NONLIN_FRSURF */
                0319 
df999eca2c Jean*0320 C-- Explicit (remaining) part of the Hyd pressure gradient:
                0321       IF ( staggerTimeStep .OR. implicitIntGravWave ) THEN
                0322         DO j=jMin,jMax
                0323          DO i=iMin,iMax
                0324            gUdPx(i,j) = -phFac*dPhiHydX(i,j) - psFac*phiSurfX(i,j)
                0325            gVdPy(i,j) = -phFac*dPhiHydY(i,j) - psFac*phiSurfY(i,j)
                0326          ENDDO
                0327         ENDDO
                0328       ELSEIF ( implicSurfPress.NE.oneRL ) THEN
                0329         DO j=jMin,jMax
                0330          DO i=iMin,iMax
                0331            gUdPx(i,j) = -psFac*phiSurfX(i,j)
                0332            gVdPy(i,j) = -psFac*phiSurfY(i,j)
                0333          ENDDO
                0334         ENDDO
                0335       ENDIF
                0336 
16f5093311 Jean*0337 #ifdef ALLOW_NONHYDROSTATIC
                0338 C-- explicit part of the NH pressure gradient is added here
                0339       IF ( use3Dsolver .AND. implicitNHPress.NE.1. _d 0 ) THEN
                0340         nhFac = pfFacMom*(1. _d 0 - implicitNHPress)
                0341      &         *recip_deepFacC(k)*recip_rhoFacC(k)
982e105a17 Jean*0342        IF ( exactConserv ) THEN
16f5093311 Jean*0343         DO j=jMin,jMax
                0344          DO i=iMin,iMax
df999eca2c Jean*0345           gUdPx(i,j) = gUdPx(i,j)
                0346      &               - nhFac*_recip_dxC(i,j,bi,bj)
                0347      &               *( (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
982e105a17 Jean*0348      &                 -( dPhiNH(i,j,bi,bj) - dPhiNH(i-1,j,bi,bj) )
                0349      &                )
df999eca2c Jean*0350           gVdPy(i,j) = gVdPy(i,j)
                0351      &               - nhFac*_recip_dyC(i,j,bi,bj)
                0352      &               *( (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
982e105a17 Jean*0353      &                 -( dPhiNH(i,j,bi,bj) - dPhiNH(i,j-1,bi,bj) )
                0354      &                )
16f5093311 Jean*0355          ENDDO
                0356         ENDDO
982e105a17 Jean*0357        ELSE
                0358         DO j=jMin,jMax
                0359          DO i=iMin,iMax
df999eca2c Jean*0360           gUdPx(i,j) = gUdPx(i,j)
                0361      &               - nhFac*_recip_dxC(i,j,bi,bj)
                0362      &               *( phi_nh(i,j,k,bi,bj) - phi_nh(i-1,j,k,bi,bj) )
                0363           gVdPy(i,j) = gVdPy(i,j)
                0364      &               - nhFac*_recip_dyC(i,j,bi,bj)
                0365      &               *( phi_nh(i,j,k,bi,bj) - phi_nh(i,j-1,k,bi,bj) )
982e105a17 Jean*0366          ENDDO
                0367         ENDDO
                0368        ENDIF
16f5093311 Jean*0369       ENDIF
71350ce0df Jean*0370 #endif /* ALLOW_NONHYDROSTATIC */
16f5093311 Jean*0371 
c0911a93b9 Jean*0372 C     Step forward zonal velocity (store in Gu)
                0373       DO j=jMin,jMax
                0374         DO i=iMin,iMax
da4112fbea Jean*0375           gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
df999eca2c Jean*0376      &     + deltaTMom*( gUtmp(i,j) + gUdPx(i,j) )
                0377      &                *_maskW(i,j,k,bi,bj)
c0911a93b9 Jean*0378         ENDDO
                0379       ENDDO
dce55d02b7 Jean*0380 
0bb99fb476 Alis*0381 C     Step forward meridional velocity (store in Gv)
fb481a83c2 Alis*0382       DO j=jMin,jMax
924557e60a Chri*0383         DO i=iMin,iMax
c0911a93b9 Jean*0384           gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
df999eca2c Jean*0385      &     + deltaTMom*( gVtmp(i,j) + gVdPy(i,j) )
                0386      &                *_maskS(i,j,k,bi,bj)
924557e60a Chri*0387         ENDDO
fb481a83c2 Alis*0388       ENDDO
459620dd74 Alis*0389 
785cc16916 Jean*0390 #ifdef ALLOW_DIAGNOSTICS
df999eca2c Jean*0391       IF ( useDiagnostics ) THEN
                0392         IF ( staggerTimeStep .OR. implicitIntGravWave ) THEN
                0393          DO j=jMin,jMax
                0394           DO i=iMin,iMax
                0395             gUdPx(i,j) = gUdPx(i,j)*_maskW(i,j,k,bi,bj)
                0396             gVdPy(i,j) = gVdPy(i,j)*_maskS(i,j,k,bi,bj)
                0397           ENDDO
                0398          ENDDO
                0399         ELSE
                0400          DO j=jMin,jMax
                0401           DO i=iMin,iMax
                0402             gUdPx(i,j) = ( gUdPx(i,j) - phFac*dPhiHydX(i,j) )
                0403      &                 *_maskW(i,j,k,bi,bj)
                0404             gVdPy(i,j) = ( gVdPy(i,j) - phFac*dPhiHydY(i,j) )
                0405      &                 *_maskS(i,j,k,bi,bj)
                0406           ENDDO
                0407          ENDDO
                0408         ENDIF
                0409         CALL DIAGNOSTICS_FILL( gUdPx,'Um_dPhiX',k,1,2,bi,bj,myThid )
                0410         CALL DIAGNOSTICS_FILL( gVdPy,'Vm_dPhiY',k,1,2,bi,bj,myThid )
                0411       ENDIF
f564e199cc Jean*0412       IF ( momViscosity .AND. useDiagnostics ) THEN
                0413         CALL DIAGNOSTICS_FILL( guDissip,'Um_Diss ',k,1,2,bi,bj,myThid )
                0414         CALL DIAGNOSTICS_FILL( gvDissip,'Vm_Diss ',k,1,2,bi,bj,myThid )
                0415       ENDIF
f6687d0ce7 Jean*0416       IF ( momForcing .AND. useDiagnostics ) THEN
                0417         CALL DIAGNOSTICS_FILL( guExt,'Um_Ext  ',k,1,2,bi,bj,myThid )
                0418         CALL DIAGNOSTICS_FILL( gvExt,'Vm_Ext  ',k,1,2,bi,bj,myThid )
2de62a9c9a Jean*0419       ENDIF
785cc16916 Jean*0420 #ifdef ALLOW_CD_CODE
                0421       IF ( useCDscheme .AND. useDiagnostics ) THEN
da4112fbea Jean*0422         CALL DIAGNOSTICS_FILL( guCor,'Um_Cori ',k,1,2,bi,bj,myThid )
                0423         CALL DIAGNOSTICS_FILL( gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid )
785cc16916 Jean*0424       ENDIF
                0425 #endif /* ALLOW_CD_CODE */
                0426 #endif /* ALLOW_DIAGNOSTICS */
                0427 
924557e60a Chri*0428       RETURN
                0429       END