model/src/timestep.F

C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.27 2001/09/26 18:09:16 cnh Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: TIMESTEP
C     !INTERFACE:
      SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
     I                     phiHyd, 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: geopotentiel)
C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
C     phiSurfY             or geopotentiel (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     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---+----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)
      psFac = pfFacMom*(1. _d 0 - implicSurfPress)
      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)
     &        )*_maskW(i,j,k,bi,bj)
        ENDDO
      ENDDO

      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)
      psFac = pfFacMom*(1. _d 0 - implicSurfPress)
      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)
     &        )*_maskS(i,j,k,bi,bj)
        ENDDO
      ENDDO

      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.28	C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.27 2001/09/26 18:09:16 cnh 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			I phiHyd, 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.23	C phiHyd - Hydrostatic Potential (ocean: pressure/rho)
32			C (atmos: geopotentiel)
33	jmc	1.21	C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
34			C phiSurfY or geopotentiel (atmos) in X and Y direction
35	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
36	adcroft	1.4	INTEGER K
37	adcroft	1.19	_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
38	jmc	1.21	_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40	adcroft	1.17	INTEGER myIter, myThid
41	heimbach	1.18
42	cnh	1.27	C !LOCAL VARIABLES:
43	cnh	1.1	C == Local variables ==
44	adcroft	1.4	INTEGER i,j
45	adcroft	1.7	_RL ab15,ab05
46	jmc	1.21	_RL phxFac,phyFac, psFac
47	jmc	1.26	_RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48			_RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49	cnh	1.27	CEOP
50	cnh	1.1
51			C Adams-Bashforth timestepping weights
52	adcroft	1.25	IF (myIter .EQ. 0) THEN
53			ab15=1.0
54			ab05=0.0
55			ELSE
56	adcroft	1.17	ab15=1.5+abeps
57			ab05=-0.5-abeps
58	adcroft	1.25	ENDIF
59	cnh	1.1
60	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
61			C- Compute effective gU term (including Adams-Bashforth weights) :
62			DO j=jMin,jMax
63			DO i=iMin,iMax
64			gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
65			& + ab05*gUNm1(i,j,k,bi,bj)
66			#ifdef INCLUDE_CD_CODE
67			& + guCD(i,j,k,bi,bj)
68			#endif
69			ENDDO
70			ENDDO
71
72			#ifdef NONLIN_FRSURF
73			IF (.NOT. vectorInvariantMomentum
74			& .AND. nonlinFreeSurf.GT.1) THEN
75	jmc	1.28	IF (select_rStar.GT.0) THEN
76			DO j=jMin,jMax
77			DO i=iMin,iMax
78			gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
79			ENDDO
80			ENDDO
81			ELSE
82			DO j=jMin,jMax
83			DO i=iMin,iMax
84			IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
85	jmc	1.26	gUtmp(i,j) = gUtmp(i,j)
86			& *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
87	jmc	1.28	ENDIF
88			ENDDO
89	jmc	1.26	ENDDO
90	jmc	1.28	ENDIF
91	jmc	1.26	ENDIF
92			#endif
93
94	adcroft	1.7	C Step forward zonal velocity (store in Gu)
95	jmc	1.21	psFac = pfFacMom*(1. _d 0 - implicSurfPress)
96	adcroft	1.19	DO j=jMin,jMax
97	cnh	1.1	DO i=iMin,iMax
98	jmc	1.26	gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
99			& +deltaTmom*(
100			& gUtmp(i,j)
101			& - psFac*phiSurfX(i,j)
102	adcroft	1.7	& )*_maskW(i,j,k,bi,bj)
103	cnh	1.1	ENDDO
104	adcroft	1.19	ENDDO
105
106			IF (staggerTimeStep) THEN
107			C-- -grad Phi_Hyd has not been incorporated to gU and is added here:
108	jmc	1.23	phxFac = pfFacMom*deltaTmom
109	adcroft	1.19	DO j=jMin,jMax
110			DO i=iMin,iMax
111			gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)
112			& - _recip_dxC(i,j,bi,bj)
113			& (phiHyd(i,j,k)-phiHyd(i-1,j,k))phxFac
114			& *_maskW(i,j,k,bi,bj)
115			ENDDO
116			ENDDO
117			ENDIF
118	adcroft	1.12
119	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
120			C- Compute effective gV term (including Adams-Bashforth weights) :
121			DO j=jMin,jMax
122			DO i=iMin,iMax
123			gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
124			& + ab05*gVNm1(i,j,k,bi,bj)
125			#ifdef INCLUDE_CD_CODE
126			& + gvCD(i,j,k,bi,bj)
127			#endif
128			ENDDO
129			ENDDO
130
131			#ifdef NONLIN_FRSURF
132			IF (.NOT. vectorInvariantMomentum
133			& .AND. nonlinFreeSurf.GT.1) THEN
134	jmc	1.28	IF (select_rStar.GT.0) THEN
135			DO j=jMin,jMax
136			DO i=iMin,iMax
137			gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
138			ENDDO
139			ENDDO
140			ELSE
141			DO j=jMin,jMax
142			DO i=iMin,iMax
143			IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
144	jmc	1.26	gVtmp(i,j) = gVtmp(i,j)
145			& *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
146	jmc	1.28	ENDIF
147			ENDDO
148	jmc	1.26	ENDDO
149	jmc	1.28	ENDIF
150	jmc	1.26	ENDIF
151			#endif
152
153	adcroft	1.7	C Step forward meridional velocity (store in Gv)
154	jmc	1.21	psFac = pfFacMom*(1. _d 0 - implicSurfPress)
155	adcroft	1.19	DO j=jMin,jMax
156	cnh	1.1	DO i=iMin,iMax
157	jmc	1.26	gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
158			& +deltaTmom*(
159			& gVtmp(i,j)
160			& - psFac*phiSurfY(i,j)
161	adcroft	1.7	& )*_maskS(i,j,k,bi,bj)
162	cnh	1.1	ENDDO
163	adcroft	1.19	ENDDO
164	adcroft	1.12
165	adcroft	1.19	IF (staggerTimeStep) THEN
166			C-- -grad Phi_Hyd has not been incorporated to gV and is added here:
167	jmc	1.23	phyFac = pfFacMom*deltaTmom
168	adcroft	1.19	DO j=jMin,jMax
169			DO i=iMin,iMax
170			gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)
171			& - _recip_dyC(i,j,bi,bj)
172	adcroft	1.24	& (phiHyd(i,j,k)-phiHyd(i,j-1,k))phyFac
173	adcroft	1.19	& *_maskS(i,j,k,bi,bj)
174			ENDDO
175	adcroft	1.12	ENDDO
176			ENDIF
177	cnh	1.1
178			RETURN
179			END