model/src/timestep.F

C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.41 2005/09/30 00:22:37 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
      LOGICAL momStartAB
      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_DIAGNOSTICS
C    Allow diagnosis of external forcing
      LOGICAL externForcDiagIsOn
      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL DIAGNOSTICS_IS_ON
      _RL     gUext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     gVext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif
#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
      momStartAB = nIter0.EQ.0
      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--  factors for gradient (X & Y directions) of Hydrostatic Potential
      phxFac = pfFacMom
      phyFac = pfFacMom

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

#ifdef ALLOW_DIAGNOSTICS
      externForcDiagIsOn = useDiagnostics .AND. momForcing
      IF ( externForcDiagIsOn ) THEN
        externForcDiagIsOn = DIAGNOSTICS_IS_ON('Um_Ext  ',myThid)
     &                  .OR. DIAGNOSTICS_IS_ON('Vm_Ext  ',myThid)
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

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

      IF ( .NOT.staggerTimeStep ) THEN
C--   Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
        DO j=jMin,jMax
         DO i=iMin,iMax
          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
         ENDDO
        ENDDO
        phxFac = 0.
        phyFac = 0.
c     ELSE
C--   Stagger time step: grad Phi_Hyp will be added later
      ENDIF 

C--   Dissipation term inside the Adams-Bashforth:
      IF ( momViscosity .AND. momDissip_In_AB) 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

C--   Forcing term inside the Adams-Bashforth:
      IF (momForcing .AND. momForcing_In_AB) THEN
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
         DO j=1,sNy+1
          DO i=1,sNx+1
           gUext(i,j) = gU(i,j,k,bi,bj)
           gVext(i,j) = gV(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */

        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)

#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
         DO j=1,sNy+1
          DO i=1,sNx+1
           gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
           gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
          ENDDO
         ENDDO
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */
      ENDIF

      IF (useCDscheme) THEN
C-     for CD-scheme, store gU,Vtmp = gU,V^n + dissip. + forcing 
        IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
          DO j=jMin,jMax
           DO i=iMin,iMax
            gUtmp(i,j) = gU(i,j,k,bi,bj) + guDissip(i,j)
            gVtmp(i,j) = gV(i,j,k,bi,bj) + gvDissip(i,j)
           ENDDO
          ENDDO
        ELSE
          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
      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                        momStartAB, myIter, myThid )
      CALL ADAMS_BASHFORTH3(
     I                        bi, bj, k,
     U                        gV, gvNm,
     I                        momStartAB, 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 recommended 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
#ifdef ALLOW_DIAGNOSTICS
       IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
        DO j=1,sNy+1
         DO i=1,sNx+1
          gUext(i,j) = gU(i,j,k,bi,bj)
          gVext(i,j) = gV(i,j,k,bi,bj)
         ENDDO
        ENDDO
       ENDIF
#endif /* ALLOW_DIAGNOSTICS */

        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)

#ifdef ALLOW_DIAGNOSTICS
       IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
        DO j=1,sNy+1
         DO i=1,sNx+1
          gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
          gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
         ENDDO
        ENDDO
       ENDIF
#endif /* ALLOW_DIAGNOSTICS */

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-|--+----|

C--   Dissipation term outside the Adams-Bashforth:
      IF ( momViscosity .AND. .NOT.momDissip_In_AB ) 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

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

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
          CALL DIAGNOSTICS_FILL(gUext,'Um_Ext  ',k,1,2,bi,bj,myThid)
          CALL DIAGNOSTICS_FILL(gVext,'Vm_Ext  ',k,1,2,bi,bj,myThid)
      ENDIF
#ifdef ALLOW_CD_CODE
      IF ( useCDscheme .AND. useDiagnostics ) THEN
          CALL DIAGNOSTICS_FILL(guCor,'Um_Cori ',k,1,2,bi,bj,myThid)
          CALL DIAGNOSTICS_FILL(gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid)
      ENDIF
#endif /* ALLOW_CD_CODE */
#endif /* ALLOW_DIAGNOSTICS */

      RETURN
      END
1	jmc	1.42	C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.41 2005/09/30 00:22:37 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	jmc	1.42	LOGICAL momStartAB
55	adcroft	1.4	INTEGER i,j
56	adcroft	1.7	_RL ab15,ab05
57	jmc	1.21	_RL phxFac,phyFac, psFac
58	jmc	1.26	_RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59			_RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60	jmc	1.41	#ifdef ALLOW_DIAGNOSTICS
61	baylor	1.40	C Allow diagnosis of external forcing
62	jmc	1.41	LOGICAL externForcDiagIsOn
63			LOGICAL DIAGNOSTICS_IS_ON
64			EXTERNAL DIAGNOSTICS_IS_ON
65	baylor	1.40	_RL gUext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66			_RL gVext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67	jmc	1.41	#endif
68	edhill	1.34	#ifdef ALLOW_CD_CODE
69	jmc	1.32	_RL guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70			_RL gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71			#endif
72	cnh	1.27	CEOP
73	cnh	1.1
74			C Adams-Bashforth timestepping weights
75	jmc	1.42	momStartAB = nIter0.EQ.0
76	adcroft	1.25	IF (myIter .EQ. 0) THEN
77			ab15=1.0
78			ab05=0.0
79			ELSE
80	adcroft	1.17	ab15=1.5+abeps
81			ab05=-0.5-abeps
82	adcroft	1.25	ENDIF
83	cnh	1.1
84	jmc	1.29	C-- explicit part of the surface potential gradient is added in this S/R
85			psFac = pfFacMom*(1. _d 0 - implicSurfPress)
86
87	jmc	1.38	C-- factors for gradient (X & Y directions) of Hydrostatic Potential
88			phxFac = pfFacMom
89			phyFac = pfFacMom
90
91	jmc	1.32	C-- including or excluding momentum forcing from Adams-Bashforth:
92			momForcing_In_AB = forcing_In_AB
93			momForcing_In_AB = .TRUE.
94	jmc	1.36	momDissip_In_AB = .TRUE.
95	jmc	1.32
96	jmc	1.41	#ifdef ALLOW_DIAGNOSTICS
97			externForcDiagIsOn = useDiagnostics .AND. momForcing
98			IF ( externForcDiagIsOn ) THEN
99			externForcDiagIsOn = DIAGNOSTICS_IS_ON('Um_Ext ',myThid)
100			& .OR. DIAGNOSTICS_IS_ON('Vm_Ext ',myThid)
101			ENDIF
102			#endif /* ALLOW_DIAGNOSTICS */
103
104	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
105	jmc	1.32
106			C- Initialize local arrays (not really necessary but safer)
107			DO j=1-Oly,sNy+Oly
108			DO i=1-Olx,sNx+Olx
109			gUtmp(i,j) = 0. _d 0
110			gVtmp(i,j) = 0. _d 0
111	edhill	1.34	#ifdef ALLOW_CD_CODE
112	jmc	1.32	guCor(i,j) = 0. _d 0
113			gvCor(i,j) = 0. _d 0
114			#endif
115			ENDDO
116			ENDDO
117
118	jmc	1.38	IF ( .NOT.staggerTimeStep ) THEN
119	jmc	1.36	C-- Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
120	jmc	1.38	DO j=jMin,jMax
121			DO i=iMin,iMax
122			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
123			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
124			ENDDO
125			ENDDO
126	jmc	1.36	phxFac = 0.
127			phyFac = 0.
128	jmc	1.38	c ELSE
129			C-- Stagger time step: grad Phi_Hyp will be added later
130	jmc	1.36	ENDIF
131
132			C-- Dissipation term inside the Adams-Bashforth:
133	jmc	1.38	IF ( momViscosity .AND. momDissip_In_AB) THEN
134	jmc	1.36	DO j=jMin,jMax
135			DO i=iMin,iMax
136			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
137			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
138			ENDDO
139			ENDDO
140			ENDIF
141
142	jmc	1.32	C-- Forcing term inside the Adams-Bashforth:
143			IF (momForcing .AND. momForcing_In_AB) THEN
144	baylor	1.40	#ifdef ALLOW_DIAGNOSTICS
145	jmc	1.41	IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
146			DO j=1,sNy+1
147			DO i=1,sNx+1
148	baylor	1.40	gUext(i,j) = gU(i,j,k,bi,bj)
149			gVext(i,j) = gV(i,j,k,bi,bj)
150			ENDDO
151			ENDDO
152			ENDIF
153			#endif /* ALLOW_DIAGNOSTICS */
154
155	jmc	1.32	CALL EXTERNAL_FORCING_U(
156			I iMin,iMax,jMin,jMax,bi,bj,k,
157			I myTime,myThid)
158			CALL EXTERNAL_FORCING_V(
159			I iMin,iMax,jMin,jMax,bi,bj,k,
160			I myTime,myThid)
161	baylor	1.40
162			#ifdef ALLOW_DIAGNOSTICS
163	jmc	1.41	IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
164			DO j=1,sNy+1
165			DO i=1,sNx+1
166	baylor	1.40	gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
167			gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
168			ENDDO
169			ENDDO
170			ENDIF
171			#endif /* ALLOW_DIAGNOSTICS */
172	jmc	1.32	ENDIF
173
174	jmc	1.37	IF (useCDscheme) THEN
175	jmc	1.42	C- for CD-scheme, store gU,Vtmp = gU,V^n + dissip. + forcing
176			IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
177			DO j=jMin,jMax
178			DO i=iMin,iMax
179			gUtmp(i,j) = gU(i,j,k,bi,bj) + guDissip(i,j)
180			gVtmp(i,j) = gV(i,j,k,bi,bj) + gvDissip(i,j)
181			ENDDO
182			ENDDO
183			ELSE
184			DO j=jMin,jMax
185			DO i=iMin,iMax
186			gUtmp(i,j) = gU(i,j,k,bi,bj)
187			gVtmp(i,j) = gV(i,j,k,bi,bj)
188			ENDDO
189			ENDDO
190			ENDIF
191	jmc	1.37	ENDIF
192
193	jmc	1.32	C- Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
194			C and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
195	jmc	1.37	#ifdef ALLOW_ADAMSBASHFORTH_3
196			CALL ADAMS_BASHFORTH3(
197			I bi, bj, k,
198			U gU, guNm,
199	jmc	1.42	I momStartAB, myIter, myThid )
200	jmc	1.37	CALL ADAMS_BASHFORTH3(
201			I bi, bj, k,
202			U gV, gvNm,
203	jmc	1.42	I momStartAB, myIter, myThid )
204	jmc	1.37	#else /* ALLOW_ADAMSBASHFORTH_3 */
205			CALL ADAMS_BASHFORTH2(
206			I bi, bj, k,
207			U gU, guNm1,
208			I myIter, myThid )
209			CALL ADAMS_BASHFORTH2(
210			I bi, bj, k,
211			U gV, gvNm1,
212			I myIter, myThid )
213			#endif /* ALLOW_ADAMSBASHFORTH_3 */
214	jmc	1.26
215	jmc	1.32	C-- Forcing term outside the Adams-Bashforth:
216	baylor	1.40	C (not recommended with CD-scheme ON)
217	jmc	1.32	IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
218	jmc	1.37	IF (useCDscheme) THEN
219			DO j=jMin,jMax
220			DO i=iMin,iMax
221			gUtmp(i,j) = gUtmp(i,j) - gU(i,j,k,bi,bj)
222			gVtmp(i,j) = gVtmp(i,j) - gV(i,j,k,bi,bj)
223			ENDDO
224			ENDDO
225			ENDIF
226	baylor	1.40	#ifdef ALLOW_DIAGNOSTICS
227	jmc	1.41	IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
228			DO j=1,sNy+1
229			DO i=1,sNx+1
230	baylor	1.40	gUext(i,j) = gU(i,j,k,bi,bj)
231			gVext(i,j) = gV(i,j,k,bi,bj)
232			ENDDO
233			ENDDO
234			ENDIF
235			#endif /* ALLOW_DIAGNOSTICS */
236
237	jmc	1.32	CALL EXTERNAL_FORCING_U(
238			I iMin,iMax,jMin,jMax,bi,bj,k,
239			I myTime,myThid)
240			CALL EXTERNAL_FORCING_V(
241			I iMin,iMax,jMin,jMax,bi,bj,k,
242			I myTime,myThid)
243	jmc	1.37
244	baylor	1.40	#ifdef ALLOW_DIAGNOSTICS
245	jmc	1.41	IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
246			DO j=1,sNy+1
247			DO i=1,sNx+1
248	baylor	1.40	gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
249			gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
250			ENDDO
251			ENDDO
252			ENDIF
253			#endif /* ALLOW_DIAGNOSTICS */
254
255	jmc	1.37	C- for CD-scheme, compute gU,Vtmp = gU,V^n + forcing
256	jmc	1.32	IF (useCDscheme) THEN
257	jmc	1.28	DO j=jMin,jMax
258			DO i=iMin,iMax
259	jmc	1.37	gUtmp(i,j) = gUtmp(i,j) + gU(i,j,k,bi,bj)
260			gVtmp(i,j) = gVtmp(i,j) + gV(i,j,k,bi,bj)
261	jmc	1.28	ENDDO
262	jmc	1.26	ENDDO
263	jmc	1.28	ENDIF
264	jmc	1.26	ENDIF
265
266	edhill	1.34	#ifdef ALLOW_CD_CODE
267	jmc	1.32	IF (useCDscheme) THEN
268			C- Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
269			C and return coriolis terms on C-grid (guCor,gvCor)
270	heimbach	1.35	CALL CD_CODE_SCHEME(
271	edhill	1.34	I bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
272			O guCor,gvCor,
273			I myTime, myIter, myThid)
274	jmc	1.32	DO j=jMin,jMax
275			DO i=iMin,iMax
276			gUtmp(i,j) = gU(i,j,k,bi,bj)
277			& + guCor(i,j)
278			gVtmp(i,j) = gV(i,j,k,bi,bj)
279			& + gvCor(i,j)
280			ENDDO
281	cnh	1.1	ENDDO
282	jmc	1.37	ELSE
283			#endif /* ALLOW_CD_CODE */
284			DO j=jMin,jMax
285			DO i=iMin,iMax
286			gUtmp(i,j) = gU(i,j,k,bi,bj)
287			gVtmp(i,j) = gV(i,j,k,bi,bj)
288			ENDDO
289			ENDDO
290			#ifdef ALLOW_CD_CODE
291	jmc	1.32	ENDIF
292	jmc	1.37	#endif
293	adcroft	1.19
294	jmc	1.26	#ifdef NONLIN_FRSURF
295			IF (.NOT. vectorInvariantMomentum
296			& .AND. nonlinFreeSurf.GT.1) THEN
297	jmc	1.28	IF (select_rStar.GT.0) THEN
298			DO j=jMin,jMax
299			DO i=iMin,iMax
300	jmc	1.32	gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
301	jmc	1.28	gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
302			ENDDO
303			ENDDO
304			ELSE
305			DO j=jMin,jMax
306			DO i=iMin,iMax
307	jmc	1.32	IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
308			gUtmp(i,j) = gUtmp(i,j)
309			& *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
310			ENDIF
311	jmc	1.28	IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
312	jmc	1.26	gVtmp(i,j) = gVtmp(i,j)
313			& *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
314	jmc	1.28	ENDIF
315			ENDDO
316	jmc	1.26	ENDDO
317	jmc	1.28	ENDIF
318	jmc	1.26	ENDIF
319	jmc	1.32	#endif /* NONLIN_FRSURF */
320
321			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
322
323	jmc	1.36	C-- Dissipation term outside the Adams-Bashforth:
324	jmc	1.38	IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
325	jmc	1.36	DO j=jMin,jMax
326			DO i=iMin,iMax
327			gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
328			gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
329			ENDDO
330			ENDDO
331			ENDIF
332
333	jmc	1.32	C Step forward zonal velocity (store in Gu)
334			DO j=jMin,jMax
335			DO i=iMin,iMax
336			gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
337			& +deltaTmom*(
338			& gUtmp(i,j)
339			& - psFac*phiSurfX(i,j)
340			& - phxFac*dPhiHydX(i,j)
341			& )*_maskW(i,j,k,bi,bj)
342			ENDDO
343			ENDDO
344	jmc	1.26
345	adcroft	1.7	C Step forward meridional velocity (store in Gv)
346	adcroft	1.19	DO j=jMin,jMax
347	cnh	1.1	DO i=iMin,iMax
348	jmc	1.32	gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
349	jmc	1.26	& +deltaTmom*(
350			& gVtmp(i,j)
351			& - psFac*phiSurfY(i,j)
352	jmc	1.29	& - phyFac*dPhiHydY(i,j)
353	adcroft	1.7	& )*_maskS(i,j,k,bi,bj)
354	cnh	1.1	ENDDO
355	adcroft	1.19	ENDDO
356	cnh	1.1
357	jmc	1.39	#ifdef ALLOW_DIAGNOSTICS
358	jmc	1.41	IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
359			CALL DIAGNOSTICS_FILL(gUext,'Um_Ext ',k,1,2,bi,bj,myThid)
360			CALL DIAGNOSTICS_FILL(gVext,'Vm_Ext ',k,1,2,bi,bj,myThid)
361			ENDIF
362	jmc	1.39	#ifdef ALLOW_CD_CODE
363			IF ( useCDscheme .AND. useDiagnostics ) THEN
364			CALL DIAGNOSTICS_FILL(guCor,'Um_Cori ',k,1,2,bi,bj,myThid)
365			CALL DIAGNOSTICS_FILL(gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid)
366			ENDIF
367			#endif /* ALLOW_CD_CODE */
368			#endif /* ALLOW_DIAGNOSTICS */
369
370	cnh	1.1	RETURN
371			END