model/src/timestep.F

C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.30 2003/02/09 02:00:50 jmc Exp $
C $Name:  $

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

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

      RETURN
      END
1	jmc	1.31	C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.30 2003/02/09 02:00:50 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.31	I 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 dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
32			C phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
33			C phiSurfY :: or geopotential (atmos) in X and Y direction
34	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
35	adcroft	1.4	INTEGER K
36	jmc	1.29	_RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37			_RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
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.29	C-- stagger time step: grad Phi_Hyp is not in gU,gV => add it in this S/R
61	jmc	1.30	IF (staggerTimeStep) THEN
62	jmc	1.29	phxFac = pfFacMom
63			phyFac = pfFacMom
64			ELSE
65			phxFac = 0.
66			phyFac = 0.
67			ENDIF
68
69			C-- explicit part of the surface potential gradient is added in this S/R
70			psFac = pfFacMom*(1. _d 0 - implicSurfPress)
71
72	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
73			C- Compute effective gU term (including Adams-Bashforth weights) :
74			DO j=jMin,jMax
75			DO i=iMin,iMax
76			gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
77			& + ab05*gUNm1(i,j,k,bi,bj)
78			#ifdef INCLUDE_CD_CODE
79			& + guCD(i,j,k,bi,bj)
80			#endif
81			ENDDO
82			ENDDO
83
84			#ifdef NONLIN_FRSURF
85			IF (.NOT. vectorInvariantMomentum
86			& .AND. nonlinFreeSurf.GT.1) THEN
87	jmc	1.28	IF (select_rStar.GT.0) THEN
88			DO j=jMin,jMax
89			DO i=iMin,iMax
90			gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
91			ENDDO
92			ENDDO
93			ELSE
94			DO j=jMin,jMax
95			DO i=iMin,iMax
96			IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
97	jmc	1.26	gUtmp(i,j) = gUtmp(i,j)
98			& *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
99	jmc	1.28	ENDIF
100			ENDDO
101	jmc	1.26	ENDDO
102	jmc	1.28	ENDIF
103	jmc	1.26	ENDIF
104			#endif
105
106	adcroft	1.7	C Step forward zonal velocity (store in Gu)
107	adcroft	1.19	DO j=jMin,jMax
108	cnh	1.1	DO i=iMin,iMax
109	jmc	1.26	gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
110			& +deltaTmom*(
111			& gUtmp(i,j)
112			& - psFac*phiSurfX(i,j)
113	jmc	1.29	& - phxFac*dPhiHydX(i,j)
114	adcroft	1.7	& )*_maskW(i,j,k,bi,bj)
115	cnh	1.1	ENDDO
116	adcroft	1.19	ENDDO
117
118	jmc	1.26	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
119			C- Compute effective gV term (including Adams-Bashforth weights) :
120			DO j=jMin,jMax
121			DO i=iMin,iMax
122			gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
123			& + ab05*gVNm1(i,j,k,bi,bj)
124			#ifdef INCLUDE_CD_CODE
125			& + gvCD(i,j,k,bi,bj)
126			#endif
127			ENDDO
128			ENDDO
129
130			#ifdef NONLIN_FRSURF
131			IF (.NOT. vectorInvariantMomentum
132			& .AND. nonlinFreeSurf.GT.1) THEN
133	jmc	1.28	IF (select_rStar.GT.0) THEN
134			DO j=jMin,jMax
135			DO i=iMin,iMax
136			gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
137			ENDDO
138			ENDDO
139			ELSE
140			DO j=jMin,jMax
141			DO i=iMin,iMax
142			IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
143	jmc	1.26	gVtmp(i,j) = gVtmp(i,j)
144			& *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
145	jmc	1.28	ENDIF
146			ENDDO
147	jmc	1.26	ENDDO
148	jmc	1.28	ENDIF
149	jmc	1.26	ENDIF
150			#endif
151
152	adcroft	1.7	C Step forward meridional velocity (store in Gv)
153	adcroft	1.19	DO j=jMin,jMax
154	cnh	1.1	DO i=iMin,iMax
155	jmc	1.26	gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
156			& +deltaTmom*(
157			& gVtmp(i,j)
158			& - psFac*phiSurfY(i,j)
159	jmc	1.29	& - phyFac*dPhiHydY(i,j)
160	adcroft	1.7	& )*_maskS(i,j,k,bi,bj)
161	cnh	1.1	ENDDO
162	adcroft	1.19	ENDDO
163	cnh	1.1
164			RETURN
165			END