model/src/timestep.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/timestep.F,v 1.25 2001/07/31 15:01:33 adcroft Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

C     /==========================================================\
C     | S/R TIMESTEP                                             |
C     | o Step model fields forward in time                      |
C     \==========================================================/
      SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
     I                     phiHyd, phiSurfX, phiSurfY,
     I                     myIter, myThid )
      IMPLICIT NONE

C     == Global variables ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"

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 ==
      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)

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

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

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.26	C $Header: /u/gcmpack/models/MITgcmUV/model/src/timestep.F,v 1.25 2001/07/31 15:01:33 adcroft Exp $
2	jmc	1.21	C $Name: $
3	cnh	1.1
4	adcroft	1.10	#include "CPP_OPTIONS.h"
5	cnh	1.1
6			C /==========================================================\
7			C \| S/R TIMESTEP \|
8			C \| o Step model fields forward in time \|
9			C \==========================================================/
10	jmc	1.21	SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
11			I phiHyd, phiSurfX, phiSurfY,
12	adcroft	1.17	I myIter, myThid )
13	jmc	1.21	IMPLICIT NONE
14	heimbach	1.18
15			C == Global variables ==
16	cnh	1.1	#include "SIZE.h"
17			#include "DYNVARS.h"
18	cnh	1.11	#include "EEPARAMS.h"
19	cnh	1.1	#include "PARAMS.h"
20			#include "GRID.h"
21	jmc	1.26	#include "SURFACE.h"
22	heimbach	1.18
23	cnh	1.1	C == Routine Arguments ==
24	jmc	1.23	C phiHyd - Hydrostatic Potential (ocean: pressure/rho)
25			C (atmos: geopotentiel)
26	jmc	1.21	C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
27			C phiSurfY or geopotentiel (atmos) in X and Y direction
28	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
29	adcroft	1.4	INTEGER K
30	adcroft	1.19	_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31	jmc	1.21	_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33	adcroft	1.17	INTEGER myIter, myThid
34	heimbach	1.18
35	cnh	1.1	C == Local variables ==
36	adcroft	1.4	INTEGER i,j
37	adcroft	1.7	_RL ab15,ab05
38	jmc	1.21	_RL phxFac,phyFac, psFac
39	jmc	1.26	_RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40			_RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41	cnh	1.1
42			C Adams-Bashforth timestepping weights
43	adcroft	1.25	IF (myIter .EQ. 0) THEN
44			ab15=1.0
45			ab05=0.0
46			ELSE
47	adcroft	1.17	ab15=1.5+abeps
48			ab05=-0.5-abeps
49	adcroft	1.25	ENDIF
50	cnh	1.1
51	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
52			C- Compute effective gU term (including Adams-Bashforth weights) :
53			DO j=jMin,jMax
54			DO i=iMin,iMax
55			gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
56			& + ab05*gUNm1(i,j,k,bi,bj)
57			#ifdef INCLUDE_CD_CODE
58			& + guCD(i,j,k,bi,bj)
59			#endif
60			ENDDO
61			ENDDO
62
63			#ifdef NONLIN_FRSURF
64			IF (.NOT. vectorInvariantMomentum
65			& .AND. nonlinFreeSurf.GT.1) THEN
66			DO j=jMin,jMax
67			DO i=iMin,iMax
68			IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
69			gUtmp(i,j) = gUtmp(i,j)
70			& *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
71			ENDIF
72			ENDDO
73			ENDDO
74			ENDIF
75			#endif
76
77	adcroft	1.7	C Step forward zonal velocity (store in Gu)
78	jmc	1.21	psFac = pfFacMom*(1. _d 0 - implicSurfPress)
79	adcroft	1.19	DO j=jMin,jMax
80	cnh	1.1	DO i=iMin,iMax
81	jmc	1.26	gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
82			& +deltaTmom*(
83			& gUtmp(i,j)
84			& - psFac*phiSurfX(i,j)
85	adcroft	1.7	& )*_maskW(i,j,k,bi,bj)
86	cnh	1.1	ENDDO
87	adcroft	1.19	ENDDO
88
89			IF (staggerTimeStep) THEN
90			C-- -grad Phi_Hyd has not been incorporated to gU and is added here:
91	jmc	1.23	phxFac = pfFacMom*deltaTmom
92	adcroft	1.19	DO j=jMin,jMax
93			DO i=iMin,iMax
94			gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)
95			& - _recip_dxC(i,j,bi,bj)
96			& (phiHyd(i,j,k)-phiHyd(i-1,j,k))phxFac
97			& *_maskW(i,j,k,bi,bj)
98			ENDDO
99			ENDDO
100			ENDIF
101	adcroft	1.12
102	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
103			C- Compute effective gV term (including Adams-Bashforth weights) :
104			DO j=jMin,jMax
105			DO i=iMin,iMax
106			gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
107			& + ab05*gVNm1(i,j,k,bi,bj)
108			#ifdef INCLUDE_CD_CODE
109			& + gvCD(i,j,k,bi,bj)
110			#endif
111			ENDDO
112			ENDDO
113
114			#ifdef NONLIN_FRSURF
115			IF (.NOT. vectorInvariantMomentum
116			& .AND. nonlinFreeSurf.GT.1) THEN
117			DO j=jMin,jMax
118			DO i=iMin,iMax
119			IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
120			gVtmp(i,j) = gVtmp(i,j)
121			& *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
122			ENDIF
123			ENDDO
124			ENDDO
125			ENDIF
126			#endif
127
128	adcroft	1.7	C Step forward meridional velocity (store in Gv)
129	jmc	1.21	psFac = pfFacMom*(1. _d 0 - implicSurfPress)
130	adcroft	1.19	DO j=jMin,jMax
131	cnh	1.1	DO i=iMin,iMax
132	jmc	1.26	gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
133			& +deltaTmom*(
134			& gVtmp(i,j)
135			& - psFac*phiSurfY(i,j)
136	adcroft	1.7	& )*_maskS(i,j,k,bi,bj)
137	cnh	1.1	ENDDO
138	adcroft	1.19	ENDDO
139	adcroft	1.12
140	adcroft	1.19	IF (staggerTimeStep) THEN
141			C-- -grad Phi_Hyd has not been incorporated to gV and is added here:
142	jmc	1.23	phyFac = pfFacMom*deltaTmom
143	adcroft	1.19	DO j=jMin,jMax
144			DO i=iMin,iMax
145			gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)
146			& - _recip_dyC(i,j,bi,bj)
147	adcroft	1.24	& (phiHyd(i,j,k)-phiHyd(i,j-1,k))phyFac
148	adcroft	1.19	& *_maskS(i,j,k,bi,bj)
149			ENDDO
150	adcroft	1.12	ENDDO
151			ENDIF
152	cnh	1.1
153			RETURN
154			END