model/src/timestep.F

C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.29 2003/02/08 02:09:20 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: TIMESTEP
C     !INTERFACE:
      SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
     I             phiHyd, dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
     I                     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     phiHyd     :: Hydrostatic Potential (ocean: pressure/rho) 
C                                         (atmos: geopotential)
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
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER K
      _RL     phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) 
      _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)
      INTEGER myIter, myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
      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)
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-- stagger time step: grad Phi_Hyp is not in gU,gV => add it in this S/R
      IF (staggerTimeStep) THEN
c     IF (.FALSE.) THEN
        phxFac = pfFacMom
        phyFac = pfFacMom
      ELSE
        phxFac = 0.
        phyFac = 0.
      ENDIF 

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

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-    Compute effective gU term (including Adams-Bashforth weights) :
      DO j=jMin,jMax
       DO i=iMin,iMax
        gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
     &             + ab05*gUNm1(i,j,k,bi,bj)  
#ifdef INCLUDE_CD_CODE
     &             + guCD(i,j,k,bi,bj)
#endif
       ENDDO
      ENDDO
      
#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)
         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
         ENDDO
        ENDDO
       ENDIF
      ENDIF
#endif 

C     Step forward zonal velocity (store in Gu)
      DO j=jMin,jMax
        DO i=iMin,iMax
          gUNm1(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

      IF (.FALSE.) THEN
c     IF (staggerTimeStep) THEN
C--   -grad Phi_Hyd has not been incorporated to gU and is added here:
        phxFac = pfFacMom*deltaTmom
        DO j=jMin,jMax
          DO i=iMin,iMax
            gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)
     &       - _recip_dxC(i,j,bi,bj)
     &         *(phiHyd(i,j,k)-phiHyd(i-1,j,k))*phxFac
     &         *_maskW(i,j,k,bi,bj)
          ENDDO
        ENDDO
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-    Compute effective gV term (including Adams-Bashforth weights) :
      DO j=jMin,jMax
       DO i=iMin,iMax
        gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
     &             + ab05*gVNm1(i,j,k,bi,bj)  
#ifdef INCLUDE_CD_CODE
     &             + gvCD(i,j,k,bi,bj)
#endif
       ENDDO
      ENDDO
      
#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
           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.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 

C     Step forward meridional velocity (store in Gv)
      DO j=jMin,jMax
        DO i=iMin,iMax
          gVNm1(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

      IF (.FALSE.) THEN
c     IF (staggerTimeStep) THEN
C--   -grad Phi_Hyd has not been incorporated to gV and is added here:
        phyFac = pfFacMom*deltaTmom
        DO j=jMin,jMax
          DO i=iMin,iMax
            gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)
     &       - _recip_dyC(i,j,bi,bj)
     &         *(phiHyd(i,j,k)-phiHyd(i,j-1,k))*phyFac
     &         *_maskS(i,j,k,bi,bj)
          ENDDO
        ENDDO
      ENDIF

      RETURN
      END
1	jmc	1.30	C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.29 2003/02/08 02:09:20 jmc Exp $
2	jmc	1.21	C $Name: $
3	cnh	1.1
4	adcroft	1.10	#include "CPP_OPTIONS.h"
5	cnh	1.1
6	cnh	1.27	CBOP
7			C !ROUTINE: TIMESTEP
8			C !INTERFACE:
9	jmc	1.21	SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
10	jmc	1.29	I phiHyd, dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
11	adcroft	1.17	I myIter, myThid )
12	cnh	1.27	C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| S/R TIMESTEP
15			C \| o Step model fields forward in time
16			C ==========================================================
17			C \ev
18
19			C !USES:
20	jmc	1.21	IMPLICIT NONE
21	heimbach	1.18	C == Global variables ==
22	cnh	1.1	#include "SIZE.h"
23			#include "DYNVARS.h"
24	cnh	1.11	#include "EEPARAMS.h"
25	cnh	1.1	#include "PARAMS.h"
26			#include "GRID.h"
27	jmc	1.26	#include "SURFACE.h"
28	heimbach	1.18
29	cnh	1.27	C !INPUT/OUTPUT PARAMETERS:
30	cnh	1.1	C == Routine Arguments ==
31	jmc	1.29	C phiHyd :: Hydrostatic Potential (ocean: pressure/rho)
32			C (atmos: geopotential)
33			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	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
37	adcroft	1.4	INTEGER K
38	adcroft	1.19	_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
39	jmc	1.29	_RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40			_RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41	jmc	1.21	_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43	adcroft	1.17	INTEGER myIter, myThid
44	heimbach	1.18
45	cnh	1.27	C !LOCAL VARIABLES:
46	cnh	1.1	C == Local variables ==
47	adcroft	1.4	INTEGER i,j
48	adcroft	1.7	_RL ab15,ab05
49	jmc	1.21	_RL phxFac,phyFac, psFac
50	jmc	1.26	_RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51			_RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52	cnh	1.27	CEOP
53	cnh	1.1
54			C Adams-Bashforth timestepping weights
55	adcroft	1.25	IF (myIter .EQ. 0) THEN
56			ab15=1.0
57			ab05=0.0
58			ELSE
59	adcroft	1.17	ab15=1.5+abeps
60			ab05=-0.5-abeps
61	adcroft	1.25	ENDIF
62	cnh	1.1
63	jmc	1.29	C-- stagger time step: grad Phi_Hyp is not in gU,gV => add it in this S/R
64	jmc	1.30	IF (staggerTimeStep) THEN
65			c IF (.FALSE.) THEN
66	jmc	1.29	phxFac = pfFacMom
67			phyFac = pfFacMom
68			ELSE
69			phxFac = 0.
70			phyFac = 0.
71			ENDIF
72
73			C-- explicit part of the surface potential gradient is added in this S/R
74			psFac = pfFacMom*(1. _d 0 - implicSurfPress)
75
76	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
77			C- Compute effective gU term (including Adams-Bashforth weights) :
78			DO j=jMin,jMax
79			DO i=iMin,iMax
80			gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
81			& + ab05*gUNm1(i,j,k,bi,bj)
82			#ifdef INCLUDE_CD_CODE
83			& + guCD(i,j,k,bi,bj)
84			#endif
85			ENDDO
86			ENDDO
87
88			#ifdef NONLIN_FRSURF
89			IF (.NOT. vectorInvariantMomentum
90			& .AND. nonlinFreeSurf.GT.1) THEN
91	jmc	1.28	IF (select_rStar.GT.0) THEN
92			DO j=jMin,jMax
93			DO i=iMin,iMax
94			gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
95			ENDDO
96			ENDDO
97			ELSE
98			DO j=jMin,jMax
99			DO i=iMin,iMax
100			IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
101	jmc	1.26	gUtmp(i,j) = gUtmp(i,j)
102			& *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
103	jmc	1.28	ENDIF
104			ENDDO
105	jmc	1.26	ENDDO
106	jmc	1.28	ENDIF
107	jmc	1.26	ENDIF
108			#endif
109
110	adcroft	1.7	C Step forward zonal velocity (store in Gu)
111	adcroft	1.19	DO j=jMin,jMax
112	cnh	1.1	DO i=iMin,iMax
113	jmc	1.26	gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
114			& +deltaTmom*(
115			& gUtmp(i,j)
116			& - psFac*phiSurfX(i,j)
117	jmc	1.29	& - phxFac*dPhiHydX(i,j)
118	adcroft	1.7	& )*_maskW(i,j,k,bi,bj)
119	cnh	1.1	ENDDO
120	adcroft	1.19	ENDDO
121
122	jmc	1.30	IF (.FALSE.) THEN
123			c IF (staggerTimeStep) THEN
124	adcroft	1.19	C-- -grad Phi_Hyd has not been incorporated to gU and is added here:
125	jmc	1.23	phxFac = pfFacMom*deltaTmom
126	adcroft	1.19	DO j=jMin,jMax
127			DO i=iMin,iMax
128			gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)
129			& - _recip_dxC(i,j,bi,bj)
130			& (phiHyd(i,j,k)-phiHyd(i-1,j,k))phxFac
131			& *_maskW(i,j,k,bi,bj)
132			ENDDO
133			ENDDO
134			ENDIF
135	adcroft	1.12
136	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
137			C- Compute effective gV term (including Adams-Bashforth weights) :
138			DO j=jMin,jMax
139			DO i=iMin,iMax
140			gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
141			& + ab05*gVNm1(i,j,k,bi,bj)
142			#ifdef INCLUDE_CD_CODE
143			& + gvCD(i,j,k,bi,bj)
144			#endif
145			ENDDO
146			ENDDO
147
148			#ifdef NONLIN_FRSURF
149			IF (.NOT. vectorInvariantMomentum
150			& .AND. nonlinFreeSurf.GT.1) THEN
151	jmc	1.28	IF (select_rStar.GT.0) THEN
152			DO j=jMin,jMax
153			DO i=iMin,iMax
154			gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
155			ENDDO
156			ENDDO
157			ELSE
158			DO j=jMin,jMax
159			DO i=iMin,iMax
160			IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
161	jmc	1.26	gVtmp(i,j) = gVtmp(i,j)
162			& *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
163	jmc	1.28	ENDIF
164			ENDDO
165	jmc	1.26	ENDDO
166	jmc	1.28	ENDIF
167	jmc	1.26	ENDIF
168			#endif
169
170	adcroft	1.7	C Step forward meridional velocity (store in Gv)
171	adcroft	1.19	DO j=jMin,jMax
172	cnh	1.1	DO i=iMin,iMax
173	jmc	1.26	gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
174			& +deltaTmom*(
175			& gVtmp(i,j)
176			& - psFac*phiSurfY(i,j)
177	jmc	1.29	& - phyFac*dPhiHydY(i,j)
178	adcroft	1.7	& )*_maskS(i,j,k,bi,bj)
179	cnh	1.1	ENDDO
180	adcroft	1.19	ENDDO
181	adcroft	1.12
182	jmc	1.30	IF (.FALSE.) THEN
183			c IF (staggerTimeStep) THEN
184	adcroft	1.19	C-- -grad Phi_Hyd has not been incorporated to gV and is added here:
185	jmc	1.23	phyFac = pfFacMom*deltaTmom
186	adcroft	1.19	DO j=jMin,jMax
187			DO i=iMin,iMax
188			gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)
189			& - _recip_dyC(i,j,bi,bj)
190	adcroft	1.24	& (phiHyd(i,j,k)-phiHyd(i,j-1,k))phyFac
191	adcroft	1.19	& *_maskS(i,j,k,bi,bj)
192			ENDDO
193	adcroft	1.12	ENDDO
194			ENDIF
195	cnh	1.1
196			RETURN
197			END