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

Annotation 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 - (hide 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 dgoldberg 1.2 C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/MANUAL/code/timestep.F,v 1.1 2015/09/10 14:56:38 dgoldberg Exp $
2 dgoldberg 1.1 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