model/src/timestep.F

C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.51 2010/01/23 00:04:03 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_CD_CODE
#include "CD_CODE_OPTIONS.h"
#endif

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 "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "RESTART.h"
#include "DYNVARS.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "NH_VARS.h"
#endif

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 ==
      INTEGER i,j
      _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_NONHYDROSTATIC
      _RL     nhFac
#endif
#ifdef ALLOW_DIAGNOSTICS
C    Allow diagnosis of external forcing
      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL DIAGNOSTICS_IS_ON
#endif
      LOGICAL externForcDiagIsOn
#if (defined ALLOW_DIAGNOSTICS) || (defined CD_CODE_NO_AB_MOMENTUM)
      _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-- explicit part of the surface potential gradient is added in this S/R
      psFac = pfFacMom*(1. _d 0 - implicSurfPress)
     &       *recip_deepFacC(k)*recip_rhoFacC(k)

C--  factors for gradient (X & Y directions) of Hydrostatic Potential
      phxFac = pfFacMom
      phyFac = pfFacMom

      externForcDiagIsOn = useDiagnostics .AND. momForcing
#ifdef ALLOW_DIAGNOSTICS
      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 .AND. .NOT. implicitIntGravWave ) 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. momForcingOutAB.NE.1 ) THEN
#ifdef ALLOW_DIAGNOSTICS
        IF ( 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 ( 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

#ifdef CD_CODE_NO_AB_MOMENTUM
      IF (useCDscheme) THEN
C-    CD-scheme, before doing AB, store gU,Vtmp = gU,V^n (+dissip. +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
#endif /* CD_CODE_NO_AB_MOMENTUM */

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                        mom_StartAB, myIter, myThid )
      CALL ADAMS_BASHFORTH3(
     I                        bi, bj, k,
     U                        gV, gvNm,
     I                        mom_StartAB, myIter, myThid )
#else /* ALLOW_ADAMSBASHFORTH_3 */
      CALL ADAMS_BASHFORTH2(
     I                        bi, bj, k,
     U                        gU, guNm1,
     I                        mom_StartAB, myIter, myThid )
      CALL ADAMS_BASHFORTH2(
     I                        bi, bj, k,
     U                        gV, gvNm1,
     I                        mom_StartAB, myIter, myThid )
#endif /* ALLOW_ADAMSBASHFORTH_3 */

C--   Forcing term outside the Adams-Bashforth:
      IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN

#if (defined ALLOW_DIAGNOSTICS) || (defined CD_CODE_NO_AB_MOMENTUM)
        IF ( externForcDiagIsOn .OR. useCDscheme ) THEN
         DO j=jMin,jMax
          DO i=iMin,iMax
           gUext(i,j) = gU(i,j,k,bi,bj)
           gVext(i,j) = gV(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDIF
#endif /* ALLOW_DIAGNOSTICS or CD_CODE_NO_AB_MOMENTUM */

        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)

#if (defined ALLOW_DIAGNOSTICS) || (defined CD_CODE_NO_AB_MOMENTUM)
        IF ( externForcDiagIsOn .OR. useCDscheme ) THEN
         DO j=jMin,jMax
          DO i=iMin,iMax
           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 or CD_CODE_NO_AB_MOMENTUM */

C--   Forcing term outside the Adams-Bashforth: end of if block
      ENDIF

#ifdef ALLOW_CD_CODE
      IF (useCDscheme) THEN
#ifdef CD_CODE_NO_AB_MOMENTUM
C-     CD-scheme, add forcing to gU,Vtmp = gU,V^n (+ dissip.)
        IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
         DO j=jMin,jMax
          DO i=iMin,iMax
           gUtmp(i,j) = gUtmp(i,j) + guExt(i,j)
           gVtmp(i,j) = gVtmp(i,j) + gvExt(i,j)
          ENDDO
         ENDDO
        ENDIF
#else /* CD_CODE_NO_AB_MOMENTUM */
C-     CD-scheme, store gU,Vtmp = gU,V^n+1/2 (+dissip.) +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 /* CD_CODE_NO_AB_MOMENTUM */
        IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
C-     CD-scheme, add dissip. to gU,Vtmp = gU,V  + forcing
         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 D-grid velocity using C-grid gU,Vtmp = gU,V +dissip +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 /* ALLOW_CD_CODE */

#ifdef NONLIN_FRSURF
      IF ( .NOT. vectorInvariantMomentum
     &     .AND. nonlinFreeSurf.GT.1 ) THEN
       IF ( select_rStar.GT.0 ) THEN
# ifndef DISABLE_RSTAR_CODE
        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
# endif /* DISABLE_RSTAR_CODE */
       ELSEIF ( selectSigmaCoord.NE.0 ) THEN
# ifndef DISABLE_SIGMA_CODE
        DO j=jMin,jMax
         DO i=iMin,iMax
           gUtmp(i,j) = gUtmp(i,j)
     &        /( 1. _d 0 + dEtaWdt(i,j,bi,bj)*deltaTfreesurf
     &                    *dBHybSigF(k)*recip_drF(k)
     &                    *recip_hFacW(i,j,k,bi,bj)
     &         )
           gVtmp(i,j) = gVtmp(i,j)
     &        /( 1. _d 0 + dEtaSdt(i,j,bi,bj)*deltaTfreesurf
     &                    *dBHybSigF(k)*recip_drF(k)
     &                    *recip_hFacS(i,j,k,bi,bj)
     &         )
         ENDDO
        ENDDO
# endif /* DISABLE_SIGMA_CODE */
       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

#ifdef ALLOW_NONHYDROSTATIC
C-- explicit part of the NH pressure gradient is added here
      IF ( use3Dsolver .AND. implicitNHPress.NE.1. _d 0 ) THEN
        nhFac = pfFacMom*(1. _d 0 - implicitNHPress)
     &         *recip_deepFacC(k)*recip_rhoFacC(k)
       IF ( exactConserv ) THEN
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gUtmp(i,j)
     &               - nhFac*_recip_dxC(i,j,bi,bj)*
     &                ( (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
     &                 -( dPhiNH(i,j,bi,bj) - dPhiNH(i-1,j,bi,bj) )
     &                )
          gVtmp(i,j) = gVtmp(i,j)
     &               - nhFac*_recip_dyC(i,j,bi,bj)*
     &                ( (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
     &                 -( dPhiNH(i,j,bi,bj) - dPhiNH(i,j-1,bi,bj) )
     &                )
         ENDDO
        ENDDO
       ELSE
        DO j=jMin,jMax
         DO i=iMin,iMax
          gUtmp(i,j) = gUtmp(i,j)
     &               - nhFac*_recip_dxC(i,j,bi,bj)*
     &                  (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
          gVtmp(i,j) = gVtmp(i,j)
     &               - nhFac*_recip_dyC(i,j,bi,bj)*
     &                  (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
         ENDDO
        ENDDO
       ENDIF
      ENDIF
#endif /* ALLOW_NONHYDROSTATIC */

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 ( 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.52	C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.51 2010/01/23 00:04:03 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	jmc	1.48	#ifdef ALLOW_CD_CODE
7			#include "CD_CODE_OPTIONS.h"
8			#endif
9	cnh	1.1
10	cnh	1.27	CBOP
11			C !ROUTINE: TIMESTEP
12			C !INTERFACE:
13	jmc	1.32	SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
14	jmc	1.31	I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
15	jmc	1.36	I guDissip, gvDissip,
16	jmc	1.32	I myTime, myIter, myThid )
17	cnh	1.27	C !DESCRIPTION: \bv
18			C ==========================================================
19	jmc	1.48	C \| S/R TIMESTEP
20			C \| o Step model fields forward in time
21	cnh	1.27	C ==========================================================
22			C \ev
23
24			C !USES:
25	jmc	1.21	IMPLICIT NONE
26	heimbach	1.18	C == Global variables ==
27	cnh	1.1	#include "SIZE.h"
28	cnh	1.11	#include "EEPARAMS.h"
29	cnh	1.1	#include "PARAMS.h"
30			#include "GRID.h"
31	jmc	1.49	#include "SURFACE.h"
32	jmc	1.47	#include "RESTART.h"
33	jmc	1.49	#include "DYNVARS.h"
34			#ifdef ALLOW_NONHYDROSTATIC
35			#include "NH_VARS.h"
36			#endif
37	heimbach	1.18
38	cnh	1.27	C !INPUT/OUTPUT PARAMETERS:
39	cnh	1.1	C == Routine Arguments ==
40	jmc	1.29	C dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
41			C phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
42			C phiSurfY :: or geopotential (atmos) in X and Y direction
43	jmc	1.36	C guDissip :: dissipation tendency (all explicit terms), u component
44			C gvDissip :: dissipation tendency (all explicit terms), v component
45
46	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
47	jmc	1.32	INTEGER k
48	jmc	1.29	_RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49			_RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50	jmc	1.21	_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52	jmc	1.36	_RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53			_RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54	jmc	1.32	_RL myTime
55	adcroft	1.17	INTEGER myIter, myThid
56	heimbach	1.18
57	cnh	1.27	C !LOCAL VARIABLES:
58	cnh	1.1	C == Local variables ==
59	adcroft	1.4	INTEGER i,j
60	jmc	1.21	_RL phxFac,phyFac, psFac
61	jmc	1.26	_RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62			_RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
64			_RL nhFac
65			#endif
66	jmc	1.41	#ifdef ALLOW_DIAGNOSTICS
67	baylor	1.40	C Allow diagnosis of external forcing
68	jmc	1.41	LOGICAL DIAGNOSTICS_IS_ON
69			EXTERNAL DIAGNOSTICS_IS_ON
70	jmc	1.48	#endif
71			LOGICAL externForcDiagIsOn
72			#if (defined ALLOW_DIAGNOSTICS) \|\| (defined CD_CODE_NO_AB_MOMENTUM)
73	baylor	1.40	_RL gUext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74			_RL gVext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75	jmc	1.41	#endif
76	edhill	1.34	#ifdef ALLOW_CD_CODE
77	jmc	1.32	_RL guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78			_RL gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79			#endif
80	cnh	1.27	CEOP
81	cnh	1.1
82	jmc	1.29	C-- explicit part of the surface potential gradient is added in this S/R
83			psFac = pfFacMom*(1. _d 0 - implicSurfPress)
84	jmc	1.49	& recip_deepFacC(k)recip_rhoFacC(k)
85	jmc	1.29
86	jmc	1.38	C-- factors for gradient (X & Y directions) of Hydrostatic Potential
87			phxFac = pfFacMom
88			phyFac = pfFacMom
89
90	jmc	1.48	externForcDiagIsOn = useDiagnostics .AND. momForcing
91	jmc	1.41	#ifdef ALLOW_DIAGNOSTICS
92			IF ( externForcDiagIsOn ) THEN
93			externForcDiagIsOn = DIAGNOSTICS_IS_ON('Um_Ext ',myThid)
94			& .OR. DIAGNOSTICS_IS_ON('Vm_Ext ',myThid)
95			ENDIF
96			#endif /* ALLOW_DIAGNOSTICS */
97
98	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
99	jmc	1.32
100			C- Initialize local arrays (not really necessary but safer)
101			DO j=1-Oly,sNy+Oly
102	jmc	1.48	DO i=1-Olx,sNx+Olx
103			gUtmp(i,j) = 0. _d 0
104			gVtmp(i,j) = 0. _d 0
105	edhill	1.34	#ifdef ALLOW_CD_CODE
106	jmc	1.48	guCor(i,j) = 0. _d 0
107			gvCor(i,j) = 0. _d 0
108	jmc	1.32	#endif
109	jmc	1.48	ENDDO
110	jmc	1.32	ENDDO
111
112	jmc	1.44	IF ( .NOT.staggerTimeStep .AND. .NOT. implicitIntGravWave ) THEN
113	jmc	1.36	C-- Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
114	jmc	1.38	DO j=jMin,jMax
115			DO i=iMin,iMax
116	jmc	1.48	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
117			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
118	jmc	1.38	ENDDO
119			ENDDO
120	jmc	1.36	phxFac = 0.
121			phyFac = 0.
122	jmc	1.38	c ELSE
123			C-- Stagger time step: grad Phi_Hyp will be added later
124	jmc	1.48	ENDIF
125	jmc	1.36
126			C-- Dissipation term inside the Adams-Bashforth:
127	jmc	1.38	IF ( momViscosity .AND. momDissip_In_AB) THEN
128	jmc	1.36	DO j=jMin,jMax
129			DO i=iMin,iMax
130	jmc	1.48	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
131			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
132	jmc	1.36	ENDDO
133			ENDDO
134			ENDIF
135
136	jmc	1.32	C-- Forcing term inside the Adams-Bashforth:
137	jmc	1.45	IF ( momForcing .AND. momForcingOutAB.NE.1 ) THEN
138	baylor	1.40	#ifdef ALLOW_DIAGNOSTICS
139	jmc	1.48	IF ( externForcDiagIsOn ) THEN
140	jmc	1.41	DO j=1,sNy+1
141			DO i=1,sNx+1
142	baylor	1.40	gUext(i,j) = gU(i,j,k,bi,bj)
143			gVext(i,j) = gV(i,j,k,bi,bj)
144			ENDDO
145			ENDDO
146			ENDIF
147			#endif /* ALLOW_DIAGNOSTICS */
148
149	jmc	1.32	CALL EXTERNAL_FORCING_U(
150			I iMin,iMax,jMin,jMax,bi,bj,k,
151			I myTime,myThid)
152			CALL EXTERNAL_FORCING_V(
153			I iMin,iMax,jMin,jMax,bi,bj,k,
154			I myTime,myThid)
155	baylor	1.40
156			#ifdef ALLOW_DIAGNOSTICS
157	jmc	1.48	IF ( externForcDiagIsOn ) THEN
158	jmc	1.41	DO j=1,sNy+1
159			DO i=1,sNx+1
160	baylor	1.40	gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
161			gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
162			ENDDO
163			ENDDO
164			ENDIF
165			#endif /* ALLOW_DIAGNOSTICS */
166	jmc	1.32	ENDIF
167
168	jmc	1.48	#ifdef CD_CODE_NO_AB_MOMENTUM
169	jmc	1.37	IF (useCDscheme) THEN
170	jmc	1.48	C- CD-scheme, before doing AB, store gU,Vtmp = gU,V^n (+dissip. +forcing)
171			DO j=jMin,jMax
172			DO i=iMin,iMax
173			gUtmp(i,j) = gU(i,j,k,bi,bj)
174			gVtmp(i,j) = gV(i,j,k,bi,bj)
175			ENDDO
176			ENDDO
177	jmc	1.37	ENDIF
178	jmc	1.48	#endif /* CD_CODE_NO_AB_MOMENTUM */
179	jmc	1.37
180	jmc	1.32	C- Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
181			C and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
182	jmc	1.37	#ifdef ALLOW_ADAMSBASHFORTH_3
183			CALL ADAMS_BASHFORTH3(
184			I bi, bj, k,
185			U gU, guNm,
186	jmc	1.47	I mom_StartAB, myIter, myThid )
187	jmc	1.37	CALL ADAMS_BASHFORTH3(
188			I bi, bj, k,
189			U gV, gvNm,
190	jmc	1.47	I mom_StartAB, myIter, myThid )
191	jmc	1.37	#else /* ALLOW_ADAMSBASHFORTH_3 */
192			CALL ADAMS_BASHFORTH2(
193			I bi, bj, k,
194			U gU, guNm1,
195	jmc	1.47	I mom_StartAB, myIter, myThid )
196	jmc	1.37	CALL ADAMS_BASHFORTH2(
197			I bi, bj, k,
198			U gV, gvNm1,
199	jmc	1.47	I mom_StartAB, myIter, myThid )
200	jmc	1.37	#endif /* ALLOW_ADAMSBASHFORTH_3 */
201	jmc	1.48
202	jmc	1.32	C-- Forcing term outside the Adams-Bashforth:
203	jmc	1.45	IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
204	jmc	1.48
205			#if (defined ALLOW_DIAGNOSTICS) \|\| (defined CD_CODE_NO_AB_MOMENTUM)
206			IF ( externForcDiagIsOn .OR. useCDscheme ) THEN
207			DO j=jMin,jMax
208			DO i=iMin,iMax
209			gUext(i,j) = gU(i,j,k,bi,bj)
210			gVext(i,j) = gV(i,j,k,bi,bj)
211			ENDDO
212	jmc	1.37	ENDDO
213	jmc	1.48	ENDIF
214			#endif /* ALLOW_DIAGNOSTICS or CD_CODE_NO_AB_MOMENTUM */
215	baylor	1.40
216	jmc	1.32	CALL EXTERNAL_FORCING_U(
217			I iMin,iMax,jMin,jMax,bi,bj,k,
218			I myTime,myThid)
219			CALL EXTERNAL_FORCING_V(
220			I iMin,iMax,jMin,jMax,bi,bj,k,
221			I myTime,myThid)
222	jmc	1.37
223	jmc	1.48	#if (defined ALLOW_DIAGNOSTICS) \|\| (defined CD_CODE_NO_AB_MOMENTUM)
224			IF ( externForcDiagIsOn .OR. useCDscheme ) THEN
225			DO j=jMin,jMax
226			DO i=iMin,iMax
227			gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
228			gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
229			ENDDO
230	baylor	1.40	ENDDO
231	jmc	1.48	ENDIF
232			#endif /* ALLOW_DIAGNOSTICS or CD_CODE_NO_AB_MOMENTUM */
233
234			C-- Forcing term outside the Adams-Bashforth: end of if block
235			ENDIF
236	baylor	1.40
237	jmc	1.48	#ifdef ALLOW_CD_CODE
238			IF (useCDscheme) THEN
239			#ifdef CD_CODE_NO_AB_MOMENTUM
240			C- CD-scheme, add forcing to gU,Vtmp = gU,V^n (+ dissip.)
241			IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
242			DO j=jMin,jMax
243			DO i=iMin,iMax
244			gUtmp(i,j) = gUtmp(i,j) + guExt(i,j)
245			gVtmp(i,j) = gVtmp(i,j) + gvExt(i,j)
246			ENDDO
247			ENDDO
248			ENDIF
249			#else /* CD_CODE_NO_AB_MOMENTUM */
250			C- CD-scheme, store gU,Vtmp = gU,V^n+1/2 (+dissip.) +forcing
251	jmc	1.28	DO j=jMin,jMax
252			DO i=iMin,iMax
253	jmc	1.48	gUtmp(i,j) = gU(i,j,k,bi,bj)
254			gVtmp(i,j) = gV(i,j,k,bi,bj)
255	jmc	1.28	ENDDO
256	jmc	1.26	ENDDO
257	jmc	1.48	#endif /* CD_CODE_NO_AB_MOMENTUM */
258			IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
259			C- CD-scheme, add dissip. to gU,Vtmp = gU,V + forcing
260			DO j=jMin,jMax
261			DO i=iMin,iMax
262			gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
263			gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
264			ENDDO
265			ENDDO
266			ENDIF
267			C- Step forward D-grid velocity using C-grid gU,Vtmp = gU,V +dissip +forcing
268	jmc	1.32	C and return coriolis terms on C-grid (guCor,gvCor)
269	heimbach	1.35	CALL CD_CODE_SCHEME(
270	edhill	1.34	I bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
271			O guCor,gvCor,
272			I myTime, myIter, myThid)
273	jmc	1.32	DO j=jMin,jMax
274			DO i=iMin,iMax
275	jmc	1.48	gUtmp(i,j) = gU(i,j,k,bi,bj)
276			& + guCor(i,j)
277			gVtmp(i,j) = gV(i,j,k,bi,bj)
278			& + gvCor(i,j)
279	jmc	1.32	ENDDO
280	cnh	1.1	ENDDO
281	jmc	1.37	ELSE
282			#endif /* ALLOW_CD_CODE */
283			DO j=jMin,jMax
284			DO i=iMin,iMax
285	jmc	1.48	gUtmp(i,j) = gU(i,j,k,bi,bj)
286			gVtmp(i,j) = gV(i,j,k,bi,bj)
287	jmc	1.37	ENDDO
288			ENDDO
289			#ifdef ALLOW_CD_CODE
290	jmc	1.32	ENDIF
291	jmc	1.48	#endif /* ALLOW_CD_CODE */
292	adcroft	1.19
293	jmc	1.26	#ifdef NONLIN_FRSURF
294	jmc	1.48	IF ( .NOT. vectorInvariantMomentum
295			& .AND. nonlinFreeSurf.GT.1 ) THEN
296	jmc	1.52	IF ( select_rStar.GT.0 ) THEN
297			# ifndef DISABLE_RSTAR_CODE
298	jmc	1.28	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	jmc	1.52	# endif /* DISABLE_RSTAR_CODE */
305			ELSEIF ( selectSigmaCoord.NE.0 ) THEN
306			# ifndef DISABLE_SIGMA_CODE
307			DO j=jMin,jMax
308			DO i=iMin,iMax
309			gUtmp(i,j) = gUtmp(i,j)
310			& /( 1. _d 0 + dEtaWdt(i,j,bi,bj)*deltaTfreesurf
311			& dBHybSigF(k)recip_drF(k)
312			& *recip_hFacW(i,j,k,bi,bj)
313			& )
314			gVtmp(i,j) = gVtmp(i,j)
315			& /( 1. _d 0 + dEtaSdt(i,j,bi,bj)*deltaTfreesurf
316			& dBHybSigF(k)recip_drF(k)
317			& *recip_hFacS(i,j,k,bi,bj)
318			& )
319			ENDDO
320			ENDDO
321			# endif /* DISABLE_SIGMA_CODE */
322	jmc	1.28	ELSE
323			DO j=jMin,jMax
324			DO i=iMin,iMax
325	jmc	1.52	IF ( k.EQ.kSurfW(i,j,bi,bj) ) THEN
326	jmc	1.32	gUtmp(i,j) = gUtmp(i,j)
327	heimbach	1.46	& *_hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
328	jmc	1.32	ENDIF
329	jmc	1.52	IF ( k.EQ.kSurfS(i,j,bi,bj) ) THEN
330	jmc	1.26	gVtmp(i,j) = gVtmp(i,j)
331	heimbach	1.46	& *_hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
332	jmc	1.28	ENDIF
333			ENDDO
334	jmc	1.26	ENDDO
335	jmc	1.28	ENDIF
336	jmc	1.26	ENDIF
337	jmc	1.32	#endif /* NONLIN_FRSURF */
338
339			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
340
341	jmc	1.36	C-- Dissipation term outside the Adams-Bashforth:
342	jmc	1.38	IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
343	jmc	1.36	DO j=jMin,jMax
344			DO i=iMin,iMax
345			gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
346			gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
347			ENDDO
348			ENDDO
349			ENDIF
350
351	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
352			C-- explicit part of the NH pressure gradient is added here
353			IF ( use3Dsolver .AND. implicitNHPress.NE.1. _d 0 ) THEN
354			nhFac = pfFacMom*(1. _d 0 - implicitNHPress)
355			& recip_deepFacC(k)recip_rhoFacC(k)
356	jmc	1.51	IF ( exactConserv ) THEN
357	jmc	1.49	DO j=jMin,jMax
358			DO i=iMin,iMax
359			gUtmp(i,j) = gUtmp(i,j)
360			& - nhFac_recip_dxC(i,j,bi,bj)
361	jmc	1.51	& ( (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
362			& -( dPhiNH(i,j,bi,bj) - dPhiNH(i-1,j,bi,bj) )
363			& )
364	jmc	1.49	gVtmp(i,j) = gVtmp(i,j)
365			& - nhFac_recip_dyC(i,j,bi,bj)
366	jmc	1.51	& ( (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
367			& -( dPhiNH(i,j,bi,bj) - dPhiNH(i,j-1,bi,bj) )
368			& )
369	jmc	1.49	ENDDO
370			ENDDO
371	jmc	1.51	ELSE
372			DO j=jMin,jMax
373			DO i=iMin,iMax
374			gUtmp(i,j) = gUtmp(i,j)
375			& - nhFac_recip_dxC(i,j,bi,bj)
376			& (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
377			gVtmp(i,j) = gVtmp(i,j)
378			& - nhFac_recip_dyC(i,j,bi,bj)
379			& (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
380			ENDDO
381			ENDDO
382			ENDIF
383	jmc	1.49	ENDIF
384	jmc	1.50	#endif /* ALLOW_NONHYDROSTATIC */
385	jmc	1.49
386	jmc	1.32	C Step forward zonal velocity (store in Gu)
387			DO j=jMin,jMax
388			DO i=iMin,iMax
389	jmc	1.48	gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
390			& +deltaTmom*(
391	jmc	1.32	& gUtmp(i,j)
392			& - psFac*phiSurfX(i,j)
393			& - phxFac*dPhiHydX(i,j)
394			& )*_maskW(i,j,k,bi,bj)
395			ENDDO
396			ENDDO
397	jmc	1.26
398	adcroft	1.7	C Step forward meridional velocity (store in Gv)
399	adcroft	1.19	DO j=jMin,jMax
400	cnh	1.1	DO i=iMin,iMax
401	jmc	1.32	gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
402	jmc	1.26	& +deltaTmom*(
403			& gVtmp(i,j)
404			& - psFac*phiSurfY(i,j)
405	jmc	1.29	& - phyFac*dPhiHydY(i,j)
406	adcroft	1.7	& )*_maskS(i,j,k,bi,bj)
407	cnh	1.1	ENDDO
408	adcroft	1.19	ENDDO
409	cnh	1.1
410	jmc	1.39	#ifdef ALLOW_DIAGNOSTICS
411	jmc	1.48	IF ( externForcDiagIsOn ) THEN
412			CALL DIAGNOSTICS_FILL( gUext,'Um_Ext ',k,1,2,bi,bj,myThid )
413			CALL DIAGNOSTICS_FILL( gVext,'Vm_Ext ',k,1,2,bi,bj,myThid )
414	jmc	1.41	ENDIF
415	jmc	1.39	#ifdef ALLOW_CD_CODE
416			IF ( useCDscheme .AND. useDiagnostics ) THEN
417	jmc	1.48	CALL DIAGNOSTICS_FILL( guCor,'Um_Cori ',k,1,2,bi,bj,myThid )
418			CALL DIAGNOSTICS_FILL( gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid )
419	jmc	1.39	ENDIF
420			#endif /* ALLOW_CD_CODE */
421			#endif /* ALLOW_DIAGNOSTICS */
422
423	cnh	1.1	RETURN
424			END