/[MITgcm]/MITgcm/model/src/timestep.F
ViewVC logotype

Contents of /MITgcm/model/src/timestep.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.29 - (show annotations) (download)
Sat Feb 8 02:09:20 2003 UTC (21 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint48d_pre
Changes since 1.28: +28 -8 lines 
in preparation for r*:
 new S/R (calc_grad_phi_hyd.F) to compute Hydrostatic potential gradient.
 pass the 2 comp. of the grad. as arguments to momentum S/R.
for the moment, only used if it does not change the results.
1 C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.28 2003/01/26 21:06:11 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: TIMESTEP
8 C !INTERFACE:
9 SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
10 I phiHyd, dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
11 I myIter, myThid )
12 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 IMPLICIT NONE
21 C == Global variables ==
22 #include "SIZE.h"
23 #include "DYNVARS.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "GRID.h"
27 #include "SURFACE.h"
28
29 C !INPUT/OUTPUT PARAMETERS:
30 C == Routine Arguments ==
31 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 INTEGER bi,bj,iMin,iMax,jMin,jMax
37 INTEGER K
38 _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
39 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40 _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42 _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43 INTEGER myIter, myThid
44
45 C !LOCAL VARIABLES:
46 C == Local variables ==
47 INTEGER i,j
48 _RL ab15,ab05
49 _RL phxFac,phyFac, psFac
50 _RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51 _RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52 CEOP
53
54 C Adams-Bashforth timestepping weights
55 IF (myIter .EQ. 0) THEN
56 ab15=1.0
57 ab05=0.0
58 ELSE
59 ab15=1.5+abeps
60 ab05=-0.5-abeps
61 ENDIF
62
63 C-- stagger time step: grad Phi_Hyp is not in gU,gV => add it in this S/R
64 c IF (staggerTimeStep) THEN
65 IF (.FALSE.) THEN
66 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 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 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 gUtmp(i,j) = gUtmp(i,j)
102 & *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
103 ENDIF
104 ENDDO
105 ENDDO
106 ENDIF
107 ENDIF
108 #endif
109
110 C Step forward zonal velocity (store in Gu)
111 DO j=jMin,jMax
112 DO i=iMin,iMax
113 gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
114 & +deltaTmom*(
115 & gUtmp(i,j)
116 & - psFac*phiSurfX(i,j)
117 & - phxFac*dPhiHydX(i,j)
118 & )*_maskW(i,j,k,bi,bj)
119 ENDDO
120 ENDDO
121
122 c IF (.FALSE.) THEN
123 IF (staggerTimeStep) THEN
124 C-- -grad Phi_Hyd has not been incorporated to gU and is added here:
125 phxFac = pfFacMom*deltaTmom
126 DO j=jMin,jMax
127 DO i=iMin,iMax
128 gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)
129 c & - phxFac*dPhiHydX(i,j)*_maskW(i,j,k,bi,bj)
130 & - _recip_dxC(i,j,bi,bj)
131 & *(phiHyd(i,j,k)-phiHyd(i-1,j,k))*phxFac
132 & *_maskW(i,j,k,bi,bj)
133 ENDDO
134 ENDDO
135 ENDIF
136
137 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
138 C- Compute effective gV term (including Adams-Bashforth weights) :
139 DO j=jMin,jMax
140 DO i=iMin,iMax
141 gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
142 & + ab05*gVNm1(i,j,k,bi,bj)
143 #ifdef INCLUDE_CD_CODE
144 & + gvCD(i,j,k,bi,bj)
145 #endif
146 ENDDO
147 ENDDO
148
149 #ifdef NONLIN_FRSURF
150 IF (.NOT. vectorInvariantMomentum
151 & .AND. nonlinFreeSurf.GT.1) THEN
152 IF (select_rStar.GT.0) THEN
153 DO j=jMin,jMax
154 DO i=iMin,iMax
155 gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
156 ENDDO
157 ENDDO
158 ELSE
159 DO j=jMin,jMax
160 DO i=iMin,iMax
161 IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
162 gVtmp(i,j) = gVtmp(i,j)
163 & *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
164 ENDIF
165 ENDDO
166 ENDDO
167 ENDIF
168 ENDIF
169 #endif
170
171 C Step forward meridional velocity (store in Gv)
172 DO j=jMin,jMax
173 DO i=iMin,iMax
174 gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
175 & +deltaTmom*(
176 & gVtmp(i,j)
177 & - psFac*phiSurfY(i,j)
178 & - phyFac*dPhiHydY(i,j)
179 & )*_maskS(i,j,k,bi,bj)
180 ENDDO
181 ENDDO
182
183 c IF (.FALSE.) THEN
184 IF (staggerTimeStep) THEN
185 C-- -grad Phi_Hyd has not been incorporated to gV and is added here:
186 phyFac = pfFacMom*deltaTmom
187 DO j=jMin,jMax
188 DO i=iMin,iMax
189 gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)
190 c & - phyFac*dPhiHydY(i,j)*_maskS(i,j,k,bi,bj)
191 & - _recip_dyC(i,j,bi,bj)
192 & *(phiHyd(i,j,k)-phiHyd(i,j-1,k))*phyFac
193 & *_maskS(i,j,k,bi,bj)
194 ENDDO
195 ENDDO
196 ENDIF
197
198 RETURN
199 END
  ViewVC Help
Powered by ViewVC 1.1.22