/[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.52 - (show annotations) (download)
Sat Sep 11 21:27:13 2010 UTC (13 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62k, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62w, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint63, checkpoint63d, checkpoint63e, checkpoint63a, checkpoint63b, checkpoint63c
Changes since 1.51: +23 -4 lines
sigma (and hybrid-sigma) coordinate code for non-linear free-surface

1 C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.51 2010/01/23 00:04:03 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6 #ifdef ALLOW_CD_CODE
7 #include "CD_CODE_OPTIONS.h"
8 #endif
9
10 CBOP
11 C !ROUTINE: TIMESTEP
12 C !INTERFACE:
13 SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
14 I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
15 I guDissip, gvDissip,
16 I myTime, myIter, myThid )
17 C !DESCRIPTION: \bv
18 C *==========================================================*
19 C | S/R TIMESTEP
20 C | o Step model fields forward in time
21 C *==========================================================*
22 C \ev
23
24 C !USES:
25 IMPLICIT NONE
26 C == Global variables ==
27 #include "SIZE.h"
28 #include "EEPARAMS.h"
29 #include "PARAMS.h"
30 #include "GRID.h"
31 #include "SURFACE.h"
32 #include "RESTART.h"
33 #include "DYNVARS.h"
34 #ifdef ALLOW_NONHYDROSTATIC
35 #include "NH_VARS.h"
36 #endif
37
38 C !INPUT/OUTPUT PARAMETERS:
39 C == Routine Arguments ==
40 C dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
41 C phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
42 C phiSurfY :: or geopotential (atmos) in X and Y direction
43 C guDissip :: dissipation tendency (all explicit terms), u component
44 C gvDissip :: dissipation tendency (all explicit terms), v component
45
46 INTEGER bi,bj,iMin,iMax,jMin,jMax
47 INTEGER k
48 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49 _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51 _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52 _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53 _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54 _RL myTime
55 INTEGER myIter, myThid
56
57 C !LOCAL VARIABLES:
58 C == Local variables ==
59 INTEGER i,j
60 _RL phxFac,phyFac, psFac
61 _RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62 _RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63 #ifdef ALLOW_NONHYDROSTATIC
64 _RL nhFac
65 #endif
66 #ifdef ALLOW_DIAGNOSTICS
67 C Allow diagnosis of external forcing
68 LOGICAL DIAGNOSTICS_IS_ON
69 EXTERNAL DIAGNOSTICS_IS_ON
70 #endif
71 LOGICAL externForcDiagIsOn
72 #if (defined ALLOW_DIAGNOSTICS) || (defined CD_CODE_NO_AB_MOMENTUM)
73 _RL gUext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 _RL gVext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 #endif
76 #ifdef ALLOW_CD_CODE
77 _RL guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78 _RL gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 #endif
80 CEOP
81
82 C-- explicit part of the surface potential gradient is added in this S/R
83 psFac = pfFacMom*(1. _d 0 - implicSurfPress)
84 & *recip_deepFacC(k)*recip_rhoFacC(k)
85
86 C-- factors for gradient (X & Y directions) of Hydrostatic Potential
87 phxFac = pfFacMom
88 phyFac = pfFacMom
89
90 externForcDiagIsOn = useDiagnostics .AND. momForcing
91 #ifdef ALLOW_DIAGNOSTICS
92 IF ( externForcDiagIsOn ) THEN
93 externForcDiagIsOn = DIAGNOSTICS_IS_ON('Um_Ext ',myThid)
94 & .OR. DIAGNOSTICS_IS_ON('Vm_Ext ',myThid)
95 ENDIF
96 #endif /* ALLOW_DIAGNOSTICS */
97
98 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
99
100 C- Initialize local arrays (not really necessary but safer)
101 DO j=1-Oly,sNy+Oly
102 DO i=1-Olx,sNx+Olx
103 gUtmp(i,j) = 0. _d 0
104 gVtmp(i,j) = 0. _d 0
105 #ifdef ALLOW_CD_CODE
106 guCor(i,j) = 0. _d 0
107 gvCor(i,j) = 0. _d 0
108 #endif
109 ENDDO
110 ENDDO
111
112 IF ( .NOT.staggerTimeStep .AND. .NOT. implicitIntGravWave ) THEN
113 C-- Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
114 DO j=jMin,jMax
115 DO i=iMin,iMax
116 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
117 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
118 ENDDO
119 ENDDO
120 phxFac = 0.
121 phyFac = 0.
122 c ELSE
123 C-- Stagger time step: grad Phi_Hyp will be added later
124 ENDIF
125
126 C-- Dissipation term inside the Adams-Bashforth:
127 IF ( momViscosity .AND. momDissip_In_AB) THEN
128 DO j=jMin,jMax
129 DO i=iMin,iMax
130 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
131 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
132 ENDDO
133 ENDDO
134 ENDIF
135
136 C-- Forcing term inside the Adams-Bashforth:
137 IF ( momForcing .AND. momForcingOutAB.NE.1 ) THEN
138 #ifdef ALLOW_DIAGNOSTICS
139 IF ( externForcDiagIsOn ) THEN
140 DO j=1,sNy+1
141 DO i=1,sNx+1
142 gUext(i,j) = gU(i,j,k,bi,bj)
143 gVext(i,j) = gV(i,j,k,bi,bj)
144 ENDDO
145 ENDDO
146 ENDIF
147 #endif /* ALLOW_DIAGNOSTICS */
148
149 CALL EXTERNAL_FORCING_U(
150 I iMin,iMax,jMin,jMax,bi,bj,k,
151 I myTime,myThid)
152 CALL EXTERNAL_FORCING_V(
153 I iMin,iMax,jMin,jMax,bi,bj,k,
154 I myTime,myThid)
155
156 #ifdef ALLOW_DIAGNOSTICS
157 IF ( externForcDiagIsOn ) THEN
158 DO j=1,sNy+1
159 DO i=1,sNx+1
160 gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
161 gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
162 ENDDO
163 ENDDO
164 ENDIF
165 #endif /* ALLOW_DIAGNOSTICS */
166 ENDIF
167
168 #ifdef CD_CODE_NO_AB_MOMENTUM
169 IF (useCDscheme) THEN
170 C- CD-scheme, before doing AB, store gU,Vtmp = gU,V^n (+dissip. +forcing)
171 DO j=jMin,jMax
172 DO i=iMin,iMax
173 gUtmp(i,j) = gU(i,j,k,bi,bj)
174 gVtmp(i,j) = gV(i,j,k,bi,bj)
175 ENDDO
176 ENDDO
177 ENDIF
178 #endif /* CD_CODE_NO_AB_MOMENTUM */
179
180 C- Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
181 C and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
182 #ifdef ALLOW_ADAMSBASHFORTH_3
183 CALL ADAMS_BASHFORTH3(
184 I bi, bj, k,
185 U gU, guNm,
186 I mom_StartAB, myIter, myThid )
187 CALL ADAMS_BASHFORTH3(
188 I bi, bj, k,
189 U gV, gvNm,
190 I mom_StartAB, myIter, myThid )
191 #else /* ALLOW_ADAMSBASHFORTH_3 */
192 CALL ADAMS_BASHFORTH2(
193 I bi, bj, k,
194 U gU, guNm1,
195 I mom_StartAB, myIter, myThid )
196 CALL ADAMS_BASHFORTH2(
197 I bi, bj, k,
198 U gV, gvNm1,
199 I mom_StartAB, myIter, myThid )
200 #endif /* ALLOW_ADAMSBASHFORTH_3 */
201
202 C-- Forcing term outside the Adams-Bashforth:
203 IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
204
205 #if (defined ALLOW_DIAGNOSTICS) || (defined CD_CODE_NO_AB_MOMENTUM)
206 IF ( externForcDiagIsOn .OR. useCDscheme ) THEN
207 DO j=jMin,jMax
208 DO i=iMin,iMax
209 gUext(i,j) = gU(i,j,k,bi,bj)
210 gVext(i,j) = gV(i,j,k,bi,bj)
211 ENDDO
212 ENDDO
213 ENDIF
214 #endif /* ALLOW_DIAGNOSTICS or CD_CODE_NO_AB_MOMENTUM */
215
216 CALL EXTERNAL_FORCING_U(
217 I iMin,iMax,jMin,jMax,bi,bj,k,
218 I myTime,myThid)
219 CALL EXTERNAL_FORCING_V(
220 I iMin,iMax,jMin,jMax,bi,bj,k,
221 I myTime,myThid)
222
223 #if (defined ALLOW_DIAGNOSTICS) || (defined CD_CODE_NO_AB_MOMENTUM)
224 IF ( externForcDiagIsOn .OR. useCDscheme ) THEN
225 DO j=jMin,jMax
226 DO i=iMin,iMax
227 gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
228 gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
229 ENDDO
230 ENDDO
231 ENDIF
232 #endif /* ALLOW_DIAGNOSTICS or CD_CODE_NO_AB_MOMENTUM */
233
234 C-- Forcing term outside the Adams-Bashforth: end of if block
235 ENDIF
236
237 #ifdef ALLOW_CD_CODE
238 IF (useCDscheme) THEN
239 #ifdef CD_CODE_NO_AB_MOMENTUM
240 C- CD-scheme, add forcing to gU,Vtmp = gU,V^n (+ dissip.)
241 IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
242 DO j=jMin,jMax
243 DO i=iMin,iMax
244 gUtmp(i,j) = gUtmp(i,j) + guExt(i,j)
245 gVtmp(i,j) = gVtmp(i,j) + gvExt(i,j)
246 ENDDO
247 ENDDO
248 ENDIF
249 #else /* CD_CODE_NO_AB_MOMENTUM */
250 C- CD-scheme, store gU,Vtmp = gU,V^n+1/2 (+dissip.) +forcing
251 DO j=jMin,jMax
252 DO i=iMin,iMax
253 gUtmp(i,j) = gU(i,j,k,bi,bj)
254 gVtmp(i,j) = gV(i,j,k,bi,bj)
255 ENDDO
256 ENDDO
257 #endif /* CD_CODE_NO_AB_MOMENTUM */
258 IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
259 C- CD-scheme, add dissip. to gU,Vtmp = gU,V + forcing
260 DO j=jMin,jMax
261 DO i=iMin,iMax
262 gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
263 gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
264 ENDDO
265 ENDDO
266 ENDIF
267 C- Step forward D-grid velocity using C-grid gU,Vtmp = gU,V +dissip +forcing
268 C and return coriolis terms on C-grid (guCor,gvCor)
269 CALL CD_CODE_SCHEME(
270 I bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
271 O guCor,gvCor,
272 I myTime, myIter, myThid)
273 DO j=jMin,jMax
274 DO i=iMin,iMax
275 gUtmp(i,j) = gU(i,j,k,bi,bj)
276 & + guCor(i,j)
277 gVtmp(i,j) = gV(i,j,k,bi,bj)
278 & + gvCor(i,j)
279 ENDDO
280 ENDDO
281 ELSE
282 #endif /* ALLOW_CD_CODE */
283 DO j=jMin,jMax
284 DO i=iMin,iMax
285 gUtmp(i,j) = gU(i,j,k,bi,bj)
286 gVtmp(i,j) = gV(i,j,k,bi,bj)
287 ENDDO
288 ENDDO
289 #ifdef ALLOW_CD_CODE
290 ENDIF
291 #endif /* ALLOW_CD_CODE */
292
293 #ifdef NONLIN_FRSURF
294 IF ( .NOT. vectorInvariantMomentum
295 & .AND. nonlinFreeSurf.GT.1 ) THEN
296 IF ( select_rStar.GT.0 ) THEN
297 # ifndef DISABLE_RSTAR_CODE
298 DO j=jMin,jMax
299 DO i=iMin,iMax
300 gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
301 gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
302 ENDDO
303 ENDDO
304 # endif /* DISABLE_RSTAR_CODE */
305 ELSEIF ( selectSigmaCoord.NE.0 ) THEN
306 # ifndef DISABLE_SIGMA_CODE
307 DO j=jMin,jMax
308 DO i=iMin,iMax
309 gUtmp(i,j) = gUtmp(i,j)
310 & /( 1. _d 0 + dEtaWdt(i,j,bi,bj)*deltaTfreesurf
311 & *dBHybSigF(k)*recip_drF(k)
312 & *recip_hFacW(i,j,k,bi,bj)
313 & )
314 gVtmp(i,j) = gVtmp(i,j)
315 & /( 1. _d 0 + dEtaSdt(i,j,bi,bj)*deltaTfreesurf
316 & *dBHybSigF(k)*recip_drF(k)
317 & *recip_hFacS(i,j,k,bi,bj)
318 & )
319 ENDDO
320 ENDDO
321 # endif /* DISABLE_SIGMA_CODE */
322 ELSE
323 DO j=jMin,jMax
324 DO i=iMin,iMax
325 IF ( k.EQ.kSurfW(i,j,bi,bj) ) THEN
326 gUtmp(i,j) = gUtmp(i,j)
327 & *_hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
328 ENDIF
329 IF ( k.EQ.kSurfS(i,j,bi,bj) ) THEN
330 gVtmp(i,j) = gVtmp(i,j)
331 & *_hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
332 ENDIF
333 ENDDO
334 ENDDO
335 ENDIF
336 ENDIF
337 #endif /* NONLIN_FRSURF */
338
339 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
340
341 C-- Dissipation term outside the Adams-Bashforth:
342 IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
343 DO j=jMin,jMax
344 DO i=iMin,iMax
345 gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
346 gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
347 ENDDO
348 ENDDO
349 ENDIF
350
351 #ifdef ALLOW_NONHYDROSTATIC
352 C-- explicit part of the NH pressure gradient is added here
353 IF ( use3Dsolver .AND. implicitNHPress.NE.1. _d 0 ) THEN
354 nhFac = pfFacMom*(1. _d 0 - implicitNHPress)
355 & *recip_deepFacC(k)*recip_rhoFacC(k)
356 IF ( exactConserv ) THEN
357 DO j=jMin,jMax
358 DO i=iMin,iMax
359 gUtmp(i,j) = gUtmp(i,j)
360 & - nhFac*_recip_dxC(i,j,bi,bj)*
361 & ( (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
362 & -( dPhiNH(i,j,bi,bj) - dPhiNH(i-1,j,bi,bj) )
363 & )
364 gVtmp(i,j) = gVtmp(i,j)
365 & - nhFac*_recip_dyC(i,j,bi,bj)*
366 & ( (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
367 & -( dPhiNH(i,j,bi,bj) - dPhiNH(i,j-1,bi,bj) )
368 & )
369 ENDDO
370 ENDDO
371 ELSE
372 DO j=jMin,jMax
373 DO i=iMin,iMax
374 gUtmp(i,j) = gUtmp(i,j)
375 & - nhFac*_recip_dxC(i,j,bi,bj)*
376 & (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
377 gVtmp(i,j) = gVtmp(i,j)
378 & - nhFac*_recip_dyC(i,j,bi,bj)*
379 & (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
380 ENDDO
381 ENDDO
382 ENDIF
383 ENDIF
384 #endif /* ALLOW_NONHYDROSTATIC */
385
386 C Step forward zonal velocity (store in Gu)
387 DO j=jMin,jMax
388 DO i=iMin,iMax
389 gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
390 & +deltaTmom*(
391 & gUtmp(i,j)
392 & - psFac*phiSurfX(i,j)
393 & - phxFac*dPhiHydX(i,j)
394 & )*_maskW(i,j,k,bi,bj)
395 ENDDO
396 ENDDO
397
398 C Step forward meridional velocity (store in Gv)
399 DO j=jMin,jMax
400 DO i=iMin,iMax
401 gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
402 & +deltaTmom*(
403 & gVtmp(i,j)
404 & - psFac*phiSurfY(i,j)
405 & - phyFac*dPhiHydY(i,j)
406 & )*_maskS(i,j,k,bi,bj)
407 ENDDO
408 ENDDO
409
410 #ifdef ALLOW_DIAGNOSTICS
411 IF ( externForcDiagIsOn ) THEN
412 CALL DIAGNOSTICS_FILL( gUext,'Um_Ext ',k,1,2,bi,bj,myThid )
413 CALL DIAGNOSTICS_FILL( gVext,'Vm_Ext ',k,1,2,bi,bj,myThid )
414 ENDIF
415 #ifdef ALLOW_CD_CODE
416 IF ( useCDscheme .AND. useDiagnostics ) THEN
417 CALL DIAGNOSTICS_FILL( guCor,'Um_Cori ',k,1,2,bi,bj,myThid )
418 CALL DIAGNOSTICS_FILL( gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid )
419 ENDIF
420 #endif /* ALLOW_CD_CODE */
421 #endif /* ALLOW_DIAGNOSTICS */
422
423 RETURN
424 END

  ViewVC Help
Powered by ViewVC 1.1.22