model/src/timestep.F

C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.36 2004/11/10 03:02:00 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: TIMESTEP
C     !INTERFACE:
      SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
     I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
     I                     guDissip, gvDissip,
     I                     myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R TIMESTEP                                              
C     | o Step model fields forward in time                       
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
C     phiSurfY ::          or geopotential (atmos) in X and Y direction
C     guDissip :: dissipation tendency (all explicit terms), u component
C     gvDissip :: dissipation tendency (all explicit terms), v component

      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER k
      _RL     dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL     dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     myTime
      INTEGER myIter, myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
      LOGICAL momForcing_In_AB
      LOGICAL momDissip_In_AB
      INTEGER i,j
      _RL ab15,ab05
      _RL phxFac,phyFac, psFac
      _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef ALLOW_CD_CODE
      _RL     guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif
CEOP

C     Adams-Bashforth timestepping weights
      IF (myIter .EQ. 0) THEN
       ab15=1.0
       ab05=0.0
      ELSE
       ab15=1.5+abeps
       ab05=-0.5-abeps
      ENDIF

C-- explicit part of the surface potential gradient is added in this S/R
      psFac = pfFacMom*(1. _d 0 - implicSurfPress)

C-- including or excluding momentum forcing from Adams-Bashforth:
      momForcing_In_AB = forcing_In_AB
      momForcing_In_AB = .TRUE.
      momDissip_In_AB  = .TRUE.

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C-    Initialize local arrays (not really necessary but safer)
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        gUtmp(i,j) = 0. _d 0
        gVtmp(i,j) = 0. _d 0
#ifdef ALLOW_CD_CODE
        guCor(i,j) = 0. _d 0
        gvCor(i,j) = 0. _d 0
#endif
       ENDDO
      ENDDO

C--   Stagger time step: grad Phi_Hyp will be added later
      IF (staggerTimeStep) THEN
        phxFac = pfFacMom
        phyFac = pfFacMom
      ELSE
C--   Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
C     note: already done in S/R mom_vecinv and mom_fluxform but would be better 
C     to add it to gU,gV here.
c       DO j=jMin,jMax
c        DO i=iMin,iMax
c         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - pfFacMom*dPhiHydX(i,j)
c         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - pfFacMom*dPhiHydY(i,j)
c        ENDDO
c       ENDDO
        phxFac = 0.
        phyFac = 0.
      ENDIF 

#ifdef ALLOW_MOM_VECINV
C--   Dissipation term inside the Adams-Bashforth:
C     note: already in gU,gV if using fluxform
      IF ( momViscosity .AND. momDissip_In_AB
     &          .AND. vectorInvariantMomentum ) THEN
        DO j=jMin,jMax
         DO i=iMin,iMax
          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
         ENDDO
        ENDDO
      ENDIF
#endif

C--   Forcing term inside the Adams-Bashforth:
      IF (momForcing .AND. momForcing_In_AB) THEN
        CALL EXTERNAL_FORCING_U(
     I     iMin,iMax,jMin,jMax,bi,bj,k,
     I     myTime,myThid)
        CALL EXTERNAL_FORCING_V(
     I     iMin,iMax,jMin,jMax,bi,bj,k,
     I     myTime,myThid)
      ENDIF

      IF (useCDscheme) THEN
C-     for CD-scheme, store gU,Vtmp = gU,V^n + forcing 
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gU(i,j,k,bi,bj)
          gVtmp(i,j) = gV(i,j,k,bi,bj)
         ENDDO
        ENDDO
      ENDIF

C-    Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
C     and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
#ifdef ALLOW_ADAMSBASHFORTH_3
      CALL ADAMS_BASHFORTH3(
     I                        bi, bj, k,
     U                        gU, guNm,
     I                        myIter, myThid )
      CALL ADAMS_BASHFORTH3(
     I                        bi, bj, k,
     U                        gV, gvNm,
     I                        myIter, myThid )
#else /* ALLOW_ADAMSBASHFORTH_3 */
      CALL ADAMS_BASHFORTH2(
     I                        bi, bj, k,
     U                        gU, guNm1,
     I                        myIter, myThid )
      CALL ADAMS_BASHFORTH2(
     I                        bi, bj, k,
     U                        gV, gvNm1,
     I                        myIter, myThid )
#endif /* ALLOW_ADAMSBASHFORTH_3 */
      
C--   Forcing term outside the Adams-Bashforth:
C     (not recommanded with CD-scheme ON)
      IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
       IF (useCDscheme) THEN
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gUtmp(i,j) - gU(i,j,k,bi,bj)
          gVtmp(i,j) = gVtmp(i,j) - gV(i,j,k,bi,bj)
         ENDDO
        ENDDO
       ENDIF
        CALL EXTERNAL_FORCING_U(
     I     iMin,iMax,jMin,jMax,bi,bj,k,
     I     myTime,myThid)
        CALL EXTERNAL_FORCING_V(
     I     iMin,iMax,jMin,jMax,bi,bj,k,
     I     myTime,myThid)

C-     for CD-scheme, compute gU,Vtmp = gU,V^n + forcing 
       IF (useCDscheme) THEN
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gUtmp(i,j) + gU(i,j,k,bi,bj)
          gVtmp(i,j) = gVtmp(i,j) + gV(i,j,k,bi,bj)
         ENDDO
        ENDDO
       ENDIF
      ENDIF

#ifdef ALLOW_CD_CODE
      IF (useCDscheme) THEN
C-     Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
C      and return coriolis terms on C-grid (guCor,gvCor)
        CALL CD_CODE_SCHEME(
     I                  bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
     O                  guCor,gvCor,
     I                  myTime, myIter, myThid)
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gU(i,j,k,bi,bj)
     &               + guCor(i,j)
          gVtmp(i,j) = gV(i,j,k,bi,bj)
     &               + gvCor(i,j)
         ENDDO
        ENDDO
      ELSE
#endif /* ALLOW_CD_CODE */
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gU(i,j,k,bi,bj)
          gVtmp(i,j) = gV(i,j,k,bi,bj)
         ENDDO
        ENDDO
#ifdef ALLOW_CD_CODE
      ENDIF
#endif

#ifdef NONLIN_FRSURF
      IF (.NOT. vectorInvariantMomentum 
     &    .AND. nonlinFreeSurf.GT.1) THEN
       IF (select_rStar.GT.0) THEN
        DO j=jMin,jMax
         DO i=iMin,iMax
           gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
           gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
         ENDDO
        ENDDO
       ELSE
        DO j=jMin,jMax
         DO i=iMin,iMax
          IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
           gUtmp(i,j) = gUtmp(i,j)
     &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
          ENDIF
          IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
           gVtmp(i,j) = gVtmp(i,j)
     &         *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
          ENDIF
         ENDDO
        ENDDO
       ENDIF
      ENDIF
#endif /* NONLIN_FRSURF */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_MOM_VECINV
C--   Dissipation term outside the Adams-Bashforth:
C     note: only implemented with vecinv formulation
      IF ( momViscosity .AND. .NOT.momDissip_In_AB
     &               .AND. vectorInvariantMomentum ) THEN
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
          gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
         ENDDO
        ENDDO
      ENDIF
#endif

C     Step forward zonal velocity (store in Gu)
      DO j=jMin,jMax
        DO i=iMin,iMax
          gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj) 
     &     +deltaTmom*( 
     &         gUtmp(i,j)
     &       - psFac*phiSurfX(i,j)
     &       - phxFac*dPhiHydX(i,j)
     &        )*_maskW(i,j,k,bi,bj)
        ENDDO
      ENDDO

C     Step forward meridional velocity (store in Gv)
      DO j=jMin,jMax
        DO i=iMin,iMax
          gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
     &     +deltaTmom*(
     &         gVtmp(i,j)
     &       - psFac*phiSurfY(i,j)
     &       - phyFac*dPhiHydY(i,j)
     &        )*_maskS(i,j,k,bi,bj)
        ENDDO
      ENDDO

      RETURN
      END
1	jmc	1.37	C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.36 2004/11/10 03:02:00 jmc Exp $
2	jmc	1.21	C $Name: $
3	cnh	1.1
4	edhill	1.34	#include "PACKAGES_CONFIG.h"
5	adcroft	1.10	#include "CPP_OPTIONS.h"
6	cnh	1.1
7	cnh	1.27	CBOP
8			C !ROUTINE: TIMESTEP
9			C !INTERFACE:
10	jmc	1.32	SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
11	jmc	1.31	I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
12	jmc	1.36	I guDissip, gvDissip,
13	jmc	1.32	I myTime, myIter, myThid )
14	cnh	1.27	C !DESCRIPTION: \bv
15			C ==========================================================
16			C \| S/R TIMESTEP
17			C \| o Step model fields forward in time
18			C ==========================================================
19			C \ev
20
21			C !USES:
22	jmc	1.21	IMPLICIT NONE
23	heimbach	1.18	C == Global variables ==
24	cnh	1.1	#include "SIZE.h"
25			#include "DYNVARS.h"
26	cnh	1.11	#include "EEPARAMS.h"
27	cnh	1.1	#include "PARAMS.h"
28			#include "GRID.h"
29	jmc	1.26	#include "SURFACE.h"
30	heimbach	1.18
31	cnh	1.27	C !INPUT/OUTPUT PARAMETERS:
32	cnh	1.1	C == Routine Arguments ==
33	jmc	1.29	C dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
34			C phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
35			C phiSurfY :: or geopotential (atmos) in X and Y direction
36	jmc	1.36	C guDissip :: dissipation tendency (all explicit terms), u component
37			C gvDissip :: dissipation tendency (all explicit terms), v component
38
39	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
40	jmc	1.32	INTEGER k
41	jmc	1.29	_RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
42			_RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43	jmc	1.21	_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45	jmc	1.36	_RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46			_RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47	jmc	1.32	_RL myTime
48	adcroft	1.17	INTEGER myIter, myThid
49	heimbach	1.18
50	cnh	1.27	C !LOCAL VARIABLES:
51	cnh	1.1	C == Local variables ==
52	jmc	1.32	LOGICAL momForcing_In_AB
53	jmc	1.36	LOGICAL momDissip_In_AB
54	adcroft	1.4	INTEGER i,j
55	adcroft	1.7	_RL ab15,ab05
56	jmc	1.21	_RL phxFac,phyFac, psFac
57	jmc	1.26	_RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58			_RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59	edhill	1.34	#ifdef ALLOW_CD_CODE
60	jmc	1.32	_RL guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61			_RL gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62			#endif
63	cnh	1.27	CEOP
64	cnh	1.1
65			C Adams-Bashforth timestepping weights
66	adcroft	1.25	IF (myIter .EQ. 0) THEN
67			ab15=1.0
68			ab05=0.0
69			ELSE
70	adcroft	1.17	ab15=1.5+abeps
71			ab05=-0.5-abeps
72	adcroft	1.25	ENDIF
73	cnh	1.1
74	jmc	1.29	C-- explicit part of the surface potential gradient is added in this S/R
75			psFac = pfFacMom*(1. _d 0 - implicSurfPress)
76
77	jmc	1.32	C-- including or excluding momentum forcing from Adams-Bashforth:
78			momForcing_In_AB = forcing_In_AB
79			momForcing_In_AB = .TRUE.
80	jmc	1.36	momDissip_In_AB = .TRUE.
81	jmc	1.32
82	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
83	jmc	1.32
84			C- Initialize local arrays (not really necessary but safer)
85			DO j=1-Oly,sNy+Oly
86			DO i=1-Olx,sNx+Olx
87			gUtmp(i,j) = 0. _d 0
88			gVtmp(i,j) = 0. _d 0
89	edhill	1.34	#ifdef ALLOW_CD_CODE
90	jmc	1.32	guCor(i,j) = 0. _d 0
91			gvCor(i,j) = 0. _d 0
92			#endif
93			ENDDO
94			ENDDO
95
96	jmc	1.36	C-- Stagger time step: grad Phi_Hyp will be added later
97			IF (staggerTimeStep) THEN
98			phxFac = pfFacMom
99			phyFac = pfFacMom
100			ELSE
101			C-- Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
102			C note: already done in S/R mom_vecinv and mom_fluxform but would be better
103			C to add it to gU,gV here.
104			c DO j=jMin,jMax
105			c DO i=iMin,iMax
106			c gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - pfFacMom*dPhiHydX(i,j)
107			c gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - pfFacMom*dPhiHydY(i,j)
108			c ENDDO
109			c ENDDO
110			phxFac = 0.
111			phyFac = 0.
112			ENDIF
113
114			#ifdef ALLOW_MOM_VECINV
115			C-- Dissipation term inside the Adams-Bashforth:
116			C note: already in gU,gV if using fluxform
117			IF ( momViscosity .AND. momDissip_In_AB
118			& .AND. vectorInvariantMomentum ) THEN
119			DO j=jMin,jMax
120			DO i=iMin,iMax
121			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
122			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
123			ENDDO
124			ENDDO
125			ENDIF
126			#endif
127
128	jmc	1.32	C-- Forcing term inside the Adams-Bashforth:
129			IF (momForcing .AND. momForcing_In_AB) THEN
130			CALL EXTERNAL_FORCING_U(
131			I iMin,iMax,jMin,jMax,bi,bj,k,
132			I myTime,myThid)
133			CALL EXTERNAL_FORCING_V(
134			I iMin,iMax,jMin,jMax,bi,bj,k,
135			I myTime,myThid)
136			ENDIF
137
138	jmc	1.37	IF (useCDscheme) THEN
139			C- for CD-scheme, store gU,Vtmp = gU,V^n + forcing
140			DO j=jMin,jMax
141			DO i=iMin,iMax
142			gUtmp(i,j) = gU(i,j,k,bi,bj)
143			gVtmp(i,j) = gV(i,j,k,bi,bj)
144			ENDDO
145			ENDDO
146			ENDIF
147
148	jmc	1.32	C- Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
149			C and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
150	jmc	1.37	#ifdef ALLOW_ADAMSBASHFORTH_3
151			CALL ADAMS_BASHFORTH3(
152			I bi, bj, k,
153			U gU, guNm,
154			I myIter, myThid )
155			CALL ADAMS_BASHFORTH3(
156			I bi, bj, k,
157			U gV, gvNm,
158			I myIter, myThid )
159			#else /* ALLOW_ADAMSBASHFORTH_3 */
160			CALL ADAMS_BASHFORTH2(
161			I bi, bj, k,
162			U gU, guNm1,
163			I myIter, myThid )
164			CALL ADAMS_BASHFORTH2(
165			I bi, bj, k,
166			U gV, gvNm1,
167			I myIter, myThid )
168			#endif /* ALLOW_ADAMSBASHFORTH_3 */
169	jmc	1.26
170	jmc	1.32	C-- Forcing term outside the Adams-Bashforth:
171			C (not recommanded with CD-scheme ON)
172			IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
173	jmc	1.37	IF (useCDscheme) THEN
174			DO j=jMin,jMax
175			DO i=iMin,iMax
176			gUtmp(i,j) = gUtmp(i,j) - gU(i,j,k,bi,bj)
177			gVtmp(i,j) = gVtmp(i,j) - gV(i,j,k,bi,bj)
178			ENDDO
179			ENDDO
180			ENDIF
181	jmc	1.32	CALL EXTERNAL_FORCING_U(
182			I iMin,iMax,jMin,jMax,bi,bj,k,
183			I myTime,myThid)
184			CALL EXTERNAL_FORCING_V(
185			I iMin,iMax,jMin,jMax,bi,bj,k,
186			I myTime,myThid)
187	jmc	1.37
188			C- for CD-scheme, compute gU,Vtmp = gU,V^n + forcing
189	jmc	1.32	IF (useCDscheme) THEN
190	jmc	1.28	DO j=jMin,jMax
191			DO i=iMin,iMax
192	jmc	1.37	gUtmp(i,j) = gUtmp(i,j) + gU(i,j,k,bi,bj)
193			gVtmp(i,j) = gVtmp(i,j) + gV(i,j,k,bi,bj)
194	jmc	1.28	ENDDO
195	jmc	1.26	ENDDO
196	jmc	1.28	ENDIF
197	jmc	1.26	ENDIF
198
199	edhill	1.34	#ifdef ALLOW_CD_CODE
200	jmc	1.32	IF (useCDscheme) THEN
201			C- Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
202			C and return coriolis terms on C-grid (guCor,gvCor)
203	heimbach	1.35	CALL CD_CODE_SCHEME(
204	edhill	1.34	I bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
205			O guCor,gvCor,
206			I myTime, myIter, myThid)
207	jmc	1.32	DO j=jMin,jMax
208			DO i=iMin,iMax
209			gUtmp(i,j) = gU(i,j,k,bi,bj)
210			& + guCor(i,j)
211			gVtmp(i,j) = gV(i,j,k,bi,bj)
212			& + gvCor(i,j)
213			ENDDO
214	cnh	1.1	ENDDO
215	jmc	1.37	ELSE
216			#endif /* ALLOW_CD_CODE */
217			DO j=jMin,jMax
218			DO i=iMin,iMax
219			gUtmp(i,j) = gU(i,j,k,bi,bj)
220			gVtmp(i,j) = gV(i,j,k,bi,bj)
221			ENDDO
222			ENDDO
223			#ifdef ALLOW_CD_CODE
224	jmc	1.32	ENDIF
225	jmc	1.37	#endif
226	adcroft	1.19
227	jmc	1.26	#ifdef NONLIN_FRSURF
228			IF (.NOT. vectorInvariantMomentum
229			& .AND. nonlinFreeSurf.GT.1) THEN
230	jmc	1.28	IF (select_rStar.GT.0) THEN
231			DO j=jMin,jMax
232			DO i=iMin,iMax
233	jmc	1.32	gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
234	jmc	1.28	gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
235			ENDDO
236			ENDDO
237			ELSE
238			DO j=jMin,jMax
239			DO i=iMin,iMax
240	jmc	1.32	IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
241			gUtmp(i,j) = gUtmp(i,j)
242			& *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
243			ENDIF
244	jmc	1.28	IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
245	jmc	1.26	gVtmp(i,j) = gVtmp(i,j)
246			& *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
247	jmc	1.28	ENDIF
248			ENDDO
249	jmc	1.26	ENDDO
250	jmc	1.28	ENDIF
251	jmc	1.26	ENDIF
252	jmc	1.32	#endif /* NONLIN_FRSURF */
253
254			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
255
256	jmc	1.36	#ifdef ALLOW_MOM_VECINV
257			C-- Dissipation term outside the Adams-Bashforth:
258			C note: only implemented with vecinv formulation
259			IF ( momViscosity .AND. .NOT.momDissip_In_AB
260			& .AND. vectorInvariantMomentum ) THEN
261			DO j=jMin,jMax
262			DO i=iMin,iMax
263			gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
264			gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
265			ENDDO
266			ENDDO
267			ENDIF
268			#endif
269
270	jmc	1.32	C Step forward zonal velocity (store in Gu)
271			DO j=jMin,jMax
272			DO i=iMin,iMax
273			gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
274			& +deltaTmom*(
275			& gUtmp(i,j)
276			& - psFac*phiSurfX(i,j)
277			& - phxFac*dPhiHydX(i,j)
278			& )*_maskW(i,j,k,bi,bj)
279			ENDDO
280			ENDDO
281	jmc	1.26
282	adcroft	1.7	C Step forward meridional velocity (store in Gv)
283	adcroft	1.19	DO j=jMin,jMax
284	cnh	1.1	DO i=iMin,iMax
285	jmc	1.32	gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
286	jmc	1.26	& +deltaTmom*(
287			& gVtmp(i,j)
288			& - psFac*phiSurfY(i,j)
289	jmc	1.29	& - phyFac*dPhiHydY(i,j)
290	adcroft	1.7	& )*_maskS(i,j,k,bi,bj)
291	cnh	1.1	ENDDO
292	adcroft	1.19	ENDDO
293	cnh	1.1
294			RETURN
295			END