/[MITgcm]/MITgcm_contrib/shelfice_remeshing/MANUAL/code/timestep.F
ViewVC logotype

Contents of /MITgcm_contrib/shelfice_remeshing/MANUAL/code/timestep.F

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


Revision 1.2 - (show annotations) (download)
Tue Oct 13 16:01:50 2015 UTC (9 years, 9 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
FILE REMOVED
*** empty log message ***

1 C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/MANUAL/code/timestep.F,v 1.1 2015/09/10 14:56:38 dgoldberg 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 C guExt :: forcing tendency, u component
60 C gvExt :: forcing tendency, v component
61 C gu_AB :: tendency increment from Adams-Bashforth, u component
62 C gv_AB :: tendency increment from Adams-Bashforth, v component
63 INTEGER i,j
64 _RL phxFac,phyFac, psFac
65 _RL guExt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66 _RL gvExt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 _RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 _RL gu_AB(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RL gv_AB(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 #ifdef ALLOW_NONHYDROSTATIC
72 _RL nhFac
73 #endif
74 #ifdef ALLOW_CD_CODE
75 _RL guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 #endif
78 CEOP
79
80 C-- explicit part of the surface potential gradient is added in this S/R
81 psFac = pfFacMom*(1. _d 0 - implicSurfPress)
82 & *recip_deepFacC(k)*recip_rhoFacC(k)
83
84 C-- factors for gradient (X & Y directions) of Hydrostatic Potential
85 phxFac = pfFacMom
86 phyFac = pfFacMom
87
88 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
89
90 C- Initialize local arrays (not really necessary for all but safer)
91 DO j=1-OLy,sNy+OLy
92 DO i=1-OLx,sNx+OLx
93 guExt(i,j) = 0. _d 0
94 gvExt(i,j) = 0. _d 0
95 gUtmp(i,j) = 0. _d 0
96 gVtmp(i,j) = 0. _d 0
97 #ifdef ALLOW_CD_CODE
98 guCor(i,j) = 0. _d 0
99 gvCor(i,j) = 0. _d 0
100 #endif
101 ENDDO
102 ENDDO
103
104 IF ( momForcing ) THEN
105 C-- Collect forcing term in local array guExt,gvExt:
106 CALL APPLY_FORCING_U(
107 U guExt,
108 I iMin,iMax,jMin,jMax, k, bi,bj,
109 I myTime, myIter, myThid )
110 CALL APPLY_FORCING_V(
111 U gvExt,
112 I iMin,iMax,jMin,jMax, k, bi,bj,
113 I myTime, myIter, myThid )
114 ENDIF
115 print *, 'JJgV1', gV(1,70,49,1,1)
116
117 IF ( .NOT.staggerTimeStep .AND. .NOT. implicitIntGravWave ) THEN
118 C-- Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
119 DO j=jMin,jMax
120 DO i=iMin,iMax
121 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
122 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
123 ENDDO
124 ENDDO
125 phxFac = 0.
126 phyFac = 0.
127 c ELSE
128 C-- Stagger time step: grad Phi_Hyp will be added later
129 ENDIF
130 print *, 'JJgV2', gV(1,70,49,1,1)
131
132 C-- Dissipation term inside the Adams-Bashforth:
133 IF ( momViscosity .AND. momDissip_In_AB) THEN
134 DO j=jMin,jMax
135 DO i=iMin,iMax
136 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
137 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
138 ENDDO
139 ENDDO
140 ENDIF
141 print *, 'JJgV3', gV(1,70,49,1,1)
142
143 C-- Forcing term inside the Adams-Bashforth:
144 IF ( momForcing .AND. momForcingOutAB.NE.1 ) THEN
145 DO j=jMin,jMax
146 DO i=iMin,iMax
147 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guExt(i,j)
148 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvExt(i,j)
149 ENDDO
150 ENDDO
151 ENDIF
152
153 #ifdef CD_CODE_NO_AB_MOMENTUM
154 IF ( useCDscheme ) THEN
155 C- CD-scheme, before doing AB, store gU,Vtmp = gU,V^n (+dissip. +forcing)
156 DO j=jMin,jMax
157 DO i=iMin,iMax
158 gUtmp(i,j) = gU(i,j,k,bi,bj)
159 gVtmp(i,j) = gV(i,j,k,bi,bj)
160 ENDDO
161 ENDDO
162 ENDIF
163 #endif /* CD_CODE_NO_AB_MOMENTUM */
164
165 print *, 'JJgV', gV(1,70,49,1,1)
166
167
168
169 C- Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
170 C and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
171 #ifdef ALLOW_ADAMSBASHFORTH_3
172 CALL ADAMS_BASHFORTH3(
173 I bi, bj, k, Nr,
174 U gU(1-OLx,1-OLy,1,bi,bj), guNm,
175 O gu_AB,
176 I mom_StartAB, myIter, myThid )
177 CALL ADAMS_BASHFORTH3(
178 I bi, bj, k, Nr,
179 U gV(1-OLx,1-OLy,1,bi,bj), gvNm,
180 O gv_AB,
181 I mom_StartAB, myIter, myThid )
182 #else /* ALLOW_ADAMSBASHFORTH_3 */
183 CALL ADAMS_BASHFORTH2(
184 I bi, bj, k, Nr,
185 U gU(1-OLx,1-OLy,1,bi,bj),
186 U guNm1(1-OLx,1-OLy,1,bi,bj),
187 O gu_AB,
188 I mom_StartAB, myIter, myThid )
189 CALL ADAMS_BASHFORTH2(
190 I bi, bj, k, Nr,
191 U gV(1-OLx,1-OLy,1,bi,bj),
192 U gvNm1(1-OLx,1-OLy,1,bi,bj),
193 O gv_AB,
194 I mom_StartAB, myIter, myThid )
195 #endif /* ALLOW_ADAMSBASHFORTH_3 */
196 #ifdef ALLOW_DIAGNOSTICS
197 IF ( useDiagnostics ) THEN
198 CALL DIAGNOSTICS_FILL(gu_AB,'AB_gU ',k,1,2,bi,bj,myThid)
199 CALL DIAGNOSTICS_FILL(gv_AB,'AB_gV ',k,1,2,bi,bj,myThid)
200 ENDIF
201 #endif /* ALLOW_DIAGNOSTICS */
202
203 C- Make a local copy in gU,Vtmp of gU,V^n+1/2 (+dissip. +forcing)
204 #ifdef CD_CODE_NO_AB_MOMENTUM
205 IF ( .NOT.useCDscheme ) THEN
206 #endif
207 DO j=jMin,jMax
208 DO i=iMin,iMax
209 gUtmp(i,j) = gU(i,j,k,bi,bj)
210 gVtmp(i,j) = gV(i,j,k,bi,bj)
211 ENDDO
212 ENDDO
213 #ifdef CD_CODE_NO_AB_MOMENTUM
214 ENDIF
215 #endif
216
217 C-- Forcing term outside the Adams-Bashforth:
218 IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
219 DO j=jMin,jMax
220 DO i=iMin,iMax
221 gUtmp(i,j) = gUtmp(i,j) + guExt(i,j)
222 gVtmp(i,j) = gVtmp(i,j) + gvExt(i,j)
223 ENDDO
224 ENDDO
225 ENDIF
226
227 C-- Dissipation term outside the Adams-Bashforth:
228 IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
229 DO j=jMin,jMax
230 DO i=iMin,iMax
231 gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
232 gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
233 ENDDO
234 ENDDO
235 ENDIF
236
237 #ifdef ALLOW_CD_CODE
238 IF ( useCDscheme ) THEN
239 C- Step forward D-grid velocity using C-grid gU,Vtmp = gU,V +dissip +forcing
240 C and return coriolis terms on C-grid (guCor,gvCor)
241 CALL CD_CODE_SCHEME(
242 I bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
243 O guCor,gvCor,
244 I myTime, myIter, myThid)
245
246 #ifdef CD_CODE_NO_AB_MOMENTUM
247 IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
248 DO j=jMin,jMax
249 DO i=iMin,iMax
250 gUtmp(i,j) = ( gU(i,j,k,bi,bj) + guExt(i,j) ) + guCor(i,j)
251 gVtmp(i,j) = ( gV(i,j,k,bi,bj) + gvExt(i,j) ) + gvCor(i,j)
252 ENDDO
253 ENDDO
254 ELSE
255 DO j=jMin,jMax
256 DO i=iMin,iMax
257 gUtmp(i,j) = gU(i,j,k,bi,bj) + guCor(i,j)
258 gVtmp(i,j) = gV(i,j,k,bi,bj) + gvCor(i,j)
259 ENDDO
260 ENDDO
261 ENDIF
262 IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
263 DO j=jMin,jMax
264 DO i=iMin,iMax
265 gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
266 gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
267 ENDDO
268 ENDDO
269 ENDIF
270 #else /* CD_CODE_NO_AB_MOMENTUM */
271 DO j=jMin,jMax
272 DO i=iMin,iMax
273 gUtmp(i,j) = gUtmp(i,j) + guCor(i,j)
274 gVtmp(i,j) = gVtmp(i,j) + gvCor(i,j)
275 ENDDO
276 ENDDO
277 #endif /* CD_CODE_NO_AB_MOMENTUM */
278 ENDIF
279 #endif /* ALLOW_CD_CODE */
280
281 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
282
283 #ifdef NONLIN_FRSURF
284 IF ( .NOT. vectorInvariantMomentum
285 & .AND. nonlinFreeSurf.GT.1 ) THEN
286 IF ( select_rStar.GT.0 ) THEN
287 # ifndef DISABLE_RSTAR_CODE
288 DO j=jMin,jMax
289 DO i=iMin,iMax
290 gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
291 gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
292 ENDDO
293 ENDDO
294 # endif /* DISABLE_RSTAR_CODE */
295 ELSEIF ( selectSigmaCoord.NE.0 ) THEN
296 # ifndef DISABLE_SIGMA_CODE
297 DO j=jMin,jMax
298 DO i=iMin,iMax
299 gUtmp(i,j) = gUtmp(i,j)
300 & /( 1. _d 0 + dEtaWdt(i,j,bi,bj)*deltaTFreeSurf
301 & *dBHybSigF(k)*recip_drF(k)
302 & *recip_hFacW(i,j,k,bi,bj)
303 & )
304 gVtmp(i,j) = gVtmp(i,j)
305 & /( 1. _d 0 + dEtaSdt(i,j,bi,bj)*deltaTFreeSurf
306 & *dBHybSigF(k)*recip_drF(k)
307 & *recip_hFacS(i,j,k,bi,bj)
308 & )
309 ENDDO
310 ENDDO
311 # endif /* DISABLE_SIGMA_CODE */
312 ELSE
313 DO j=jMin,jMax
314 DO i=iMin,iMax
315 IF ( k.EQ.kSurfW(i,j,bi,bj) ) THEN
316 gUtmp(i,j) = gUtmp(i,j)
317 & *_hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
318 ENDIF
319 IF ( k.EQ.kSurfS(i,j,bi,bj) ) THEN
320 gVtmp(i,j) = gVtmp(i,j)
321 & *_hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
322 ENDIF
323 ENDDO
324 ENDDO
325 ENDIF
326 ENDIF
327 #endif /* NONLIN_FRSURF */
328
329 #ifdef ALLOW_NONHYDROSTATIC
330 C-- explicit part of the NH pressure gradient is added here
331 IF ( use3Dsolver .AND. implicitNHPress.NE.1. _d 0 ) THEN
332 nhFac = pfFacMom*(1. _d 0 - implicitNHPress)
333 & *recip_deepFacC(k)*recip_rhoFacC(k)
334 IF ( exactConserv ) THEN
335 DO j=jMin,jMax
336 DO i=iMin,iMax
337 gUtmp(i,j) = gUtmp(i,j)
338 & - nhFac*_recip_dxC(i,j,bi,bj)*
339 & ( (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
340 & -( dPhiNH(i,j,bi,bj) - dPhiNH(i-1,j,bi,bj) )
341 & )
342 gVtmp(i,j) = gVtmp(i,j)
343 & - nhFac*_recip_dyC(i,j,bi,bj)*
344 & ( (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
345 & -( dPhiNH(i,j,bi,bj) - dPhiNH(i,j-1,bi,bj) )
346 & )
347 ENDDO
348 ENDDO
349 ELSE
350 DO j=jMin,jMax
351 DO i=iMin,iMax
352 gUtmp(i,j) = gUtmp(i,j)
353 & - nhFac*_recip_dxC(i,j,bi,bj)*
354 & (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
355 gVtmp(i,j) = gVtmp(i,j)
356 & - nhFac*_recip_dyC(i,j,bi,bj)*
357 & (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
358 ENDDO
359 ENDDO
360 ENDIF
361 ENDIF
362 #endif /* ALLOW_NONHYDROSTATIC */
363
364 C Step forward zonal velocity (store in Gu)
365 DO j=jMin,jMax
366 DO i=iMin,iMax
367 gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
368 & +deltaTMom*(
369 & gUtmp(i,j)
370 & - psFac*phiSurfX(i,j)
371 & - phxFac*dPhiHydX(i,j)
372 & )*_maskW(i,j,k,bi,bj)
373 ENDDO
374 ENDDO
375
376 C Step forward meridional velocity (store in Gv)
377 DO j=jMin,jMax
378 DO i=iMin,iMax
379 gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
380 & +deltaTMom*(
381 & gVtmp(i,j)
382 & - psFac*phiSurfY(i,j)
383 & - phyFac*dPhiHydY(i,j)
384 & )*_maskS(i,j,k,bi,bj)
385 ENDDO
386 ENDDO
387
388 #ifdef ALLOW_DIAGNOSTICS
389 IF ( momViscosity .AND. useDiagnostics ) THEN
390 CALL DIAGNOSTICS_FILL( guDissip,'Um_Diss ',k,1,2,bi,bj,myThid )
391 CALL DIAGNOSTICS_FILL( gvDissip,'Vm_Diss ',k,1,2,bi,bj,myThid )
392 ENDIF
393 IF ( momForcing .AND. useDiagnostics ) THEN
394 CALL DIAGNOSTICS_FILL( guExt,'Um_Ext ',k,1,2,bi,bj,myThid )
395 CALL DIAGNOSTICS_FILL( gvExt,'Vm_Ext ',k,1,2,bi,bj,myThid )
396 ENDIF
397 #ifdef ALLOW_CD_CODE
398 IF ( useCDscheme .AND. useDiagnostics ) THEN
399 CALL DIAGNOSTICS_FILL( guCor,'Um_Cori ',k,1,2,bi,bj,myThid )
400 CALL DIAGNOSTICS_FILL( gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid )
401 ENDIF
402 #endif /* ALLOW_CD_CODE */
403 #endif /* ALLOW_DIAGNOSTICS */
404
405 RETURN
406 END

  ViewVC Help
Powered by ViewVC 1.1.22