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

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

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

Revision 1.46 - (hide annotations) (download)
Wed Jun 7 01:55:13 2006 UTC (17 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59, checkpoint58y_post, checkpoint58t_post, checkpoint58m_post, checkpoint58w_post, checkpoint58o_post, checkpoint58p_post, checkpoint58q_post, mitgcm_mapl_00, checkpoint58r_post, checkpoint58n_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59h, checkpoint58k_post, checkpoint58v_post, checkpoint58l_post, checkpoint58x_post, checkpoint58h_post, checkpoint58j_post, checkpoint58i_post, checkpoint58u_post, checkpoint58s_post
Changes since 1.45: +3 -3 lines 
Modifications for bottom topography control
o replace hFacC by _hFacC at various places
o replace ALLOW_HFACC_CONTROL by ALLOW_DEPTH_CONTROL
o add non-self-adjoint cg2d_nsa
o update autodiff support routines
o re-initialise hfac after ctrl_depth_ini
o works for 5x5 box, doesnt work for global_ocean.90x40x15
1 heimbach 1.46 C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.45 2006/03/07 15:28:02 jmc Exp $
2 jmc 1.21 C $Name: $
3 cnh 1.1
4 edhill 1.34 #include "PACKAGES_CONFIG.h"
5 adcroft 1.10 #include "CPP_OPTIONS.h"
6 cnh 1.1
7 cnh 1.27 CBOP
8     C !ROUTINE: TIMESTEP
9     C !INTERFACE:
10 jmc 1.32 SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
11 jmc 1.31 I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
12 jmc 1.36 I guDissip, gvDissip,
13 jmc 1.32 I myTime, myIter, myThid )
14 cnh 1.27 C !DESCRIPTION: \bv
15     C *==========================================================*
16     C | S/R TIMESTEP
17     C | o Step model fields forward in time
18     C *==========================================================*
19     C \ev
20    
21     C !USES:
22 jmc 1.21 IMPLICIT NONE
23 heimbach 1.18 C == Global variables ==
24 cnh 1.1 #include "SIZE.h"
25     #include "DYNVARS.h"
26 cnh 1.11 #include "EEPARAMS.h"
27 cnh 1.1 #include "PARAMS.h"
28     #include "GRID.h"
29 jmc 1.26 #include "SURFACE.h"
30 heimbach 1.18
31 cnh 1.27 C !INPUT/OUTPUT PARAMETERS:
32 cnh 1.1 C == Routine Arguments ==
33 jmc 1.29 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 jmc 1.36 C guDissip :: dissipation tendency (all explicit terms), u component
37     C gvDissip :: dissipation tendency (all explicit terms), v component
38    
39 cnh 1.1 INTEGER bi,bj,iMin,iMax,jMin,jMax
40 jmc 1.32 INTEGER k
41 jmc 1.29 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
42     _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43 jmc 1.21 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44     _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45 jmc 1.36 _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46     _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47 jmc 1.32 _RL myTime
48 adcroft 1.17 INTEGER myIter, myThid
49 heimbach 1.18
50 cnh 1.27 C !LOCAL VARIABLES:
51 cnh 1.1 C == Local variables ==
52 jmc 1.42 LOGICAL momStartAB
53 adcroft 1.4 INTEGER i,j
54 jmc 1.21 _RL phxFac,phyFac, psFac
55 jmc 1.26 _RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56     _RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57 jmc 1.41 #ifdef ALLOW_DIAGNOSTICS
58 baylor 1.40 C Allow diagnosis of external forcing
59 jmc 1.41 LOGICAL externForcDiagIsOn
60     LOGICAL DIAGNOSTICS_IS_ON
61     EXTERNAL DIAGNOSTICS_IS_ON
62 baylor 1.40 _RL gUext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63     _RL gVext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 jmc 1.41 #endif
65 edhill 1.34 #ifdef ALLOW_CD_CODE
66 jmc 1.32 _RL guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RL gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     #endif
69 cnh 1.27 CEOP
70 cnh 1.1
71 jmc 1.43 C Start AB with low-order timestepping weights
72 jmc 1.42 momStartAB = nIter0.EQ.0
73 cnh 1.1
74 jmc 1.29 C-- explicit part of the surface potential gradient is added in this S/R
75     psFac = pfFacMom*(1. _d 0 - implicSurfPress)
76    
77 jmc 1.38 C-- factors for gradient (X & Y directions) of Hydrostatic Potential
78     phxFac = pfFacMom
79     phyFac = pfFacMom
80    
81 jmc 1.41 #ifdef ALLOW_DIAGNOSTICS
82     externForcDiagIsOn = useDiagnostics .AND. momForcing
83     IF ( externForcDiagIsOn ) THEN
84     externForcDiagIsOn = DIAGNOSTICS_IS_ON('Um_Ext ',myThid)
85     & .OR. DIAGNOSTICS_IS_ON('Vm_Ext ',myThid)
86     ENDIF
87     #endif /* ALLOW_DIAGNOSTICS */
88    
89 jmc 1.26 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
90 jmc 1.32
91     C- Initialize local arrays (not really necessary but safer)
92     DO j=1-Oly,sNy+Oly
93     DO i=1-Olx,sNx+Olx
94     gUtmp(i,j) = 0. _d 0
95     gVtmp(i,j) = 0. _d 0
96 edhill 1.34 #ifdef ALLOW_CD_CODE
97 jmc 1.32 guCor(i,j) = 0. _d 0
98     gvCor(i,j) = 0. _d 0
99     #endif
100     ENDDO
101     ENDDO
102    
103 jmc 1.44 IF ( .NOT.staggerTimeStep .AND. .NOT. implicitIntGravWave ) THEN
104 jmc 1.36 C-- Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
105 jmc 1.38 DO j=jMin,jMax
106     DO i=iMin,iMax
107     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
108     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
109     ENDDO
110     ENDDO
111 jmc 1.36 phxFac = 0.
112     phyFac = 0.
113 jmc 1.38 c ELSE
114     C-- Stagger time step: grad Phi_Hyp will be added later
115 jmc 1.36 ENDIF
116    
117     C-- Dissipation term inside the Adams-Bashforth:
118 jmc 1.38 IF ( momViscosity .AND. momDissip_In_AB) THEN
119 jmc 1.36 DO j=jMin,jMax
120     DO i=iMin,iMax
121     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
122     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
123     ENDDO
124     ENDDO
125     ENDIF
126    
127 jmc 1.32 C-- Forcing term inside the Adams-Bashforth:
128 jmc 1.45 IF ( momForcing .AND. momForcingOutAB.NE.1 ) THEN
129 baylor 1.40 #ifdef ALLOW_DIAGNOSTICS
130 jmc 1.41 IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
131     DO j=1,sNy+1
132     DO i=1,sNx+1
133 baylor 1.40 gUext(i,j) = gU(i,j,k,bi,bj)
134     gVext(i,j) = gV(i,j,k,bi,bj)
135     ENDDO
136     ENDDO
137     ENDIF
138     #endif /* ALLOW_DIAGNOSTICS */
139    
140 jmc 1.32 CALL EXTERNAL_FORCING_U(
141     I iMin,iMax,jMin,jMax,bi,bj,k,
142     I myTime,myThid)
143     CALL EXTERNAL_FORCING_V(
144     I iMin,iMax,jMin,jMax,bi,bj,k,
145     I myTime,myThid)
146 baylor 1.40
147     #ifdef ALLOW_DIAGNOSTICS
148 jmc 1.41 IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
149     DO j=1,sNy+1
150     DO i=1,sNx+1
151 baylor 1.40 gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
152     gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
153     ENDDO
154     ENDDO
155     ENDIF
156     #endif /* ALLOW_DIAGNOSTICS */
157 jmc 1.32 ENDIF
158    
159 jmc 1.37 IF (useCDscheme) THEN
160 jmc 1.42 C- for CD-scheme, store gU,Vtmp = gU,V^n + dissip. + forcing
161     IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
162     DO j=jMin,jMax
163     DO i=iMin,iMax
164     gUtmp(i,j) = gU(i,j,k,bi,bj) + guDissip(i,j)
165     gVtmp(i,j) = gV(i,j,k,bi,bj) + gvDissip(i,j)
166     ENDDO
167     ENDDO
168     ELSE
169     DO j=jMin,jMax
170     DO i=iMin,iMax
171     gUtmp(i,j) = gU(i,j,k,bi,bj)
172     gVtmp(i,j) = gV(i,j,k,bi,bj)
173     ENDDO
174     ENDDO
175     ENDIF
176 jmc 1.37 ENDIF
177    
178 jmc 1.32 C- Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
179     C and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
180 jmc 1.37 #ifdef ALLOW_ADAMSBASHFORTH_3
181     CALL ADAMS_BASHFORTH3(
182     I bi, bj, k,
183     U gU, guNm,
184 jmc 1.42 I momStartAB, myIter, myThid )
185 jmc 1.37 CALL ADAMS_BASHFORTH3(
186     I bi, bj, k,
187     U gV, gvNm,
188 jmc 1.42 I momStartAB, myIter, myThid )
189 jmc 1.37 #else /* ALLOW_ADAMSBASHFORTH_3 */
190     CALL ADAMS_BASHFORTH2(
191     I bi, bj, k,
192     U gU, guNm1,
193     I myIter, myThid )
194     CALL ADAMS_BASHFORTH2(
195     I bi, bj, k,
196     U gV, gvNm1,
197     I myIter, myThid )
198     #endif /* ALLOW_ADAMSBASHFORTH_3 */
199 jmc 1.26
200 jmc 1.32 C-- Forcing term outside the Adams-Bashforth:
201 baylor 1.40 C (not recommended with CD-scheme ON)
202 jmc 1.45 IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
203 jmc 1.37 IF (useCDscheme) THEN
204     DO j=jMin,jMax
205     DO i=iMin,iMax
206     gUtmp(i,j) = gUtmp(i,j) - gU(i,j,k,bi,bj)
207     gVtmp(i,j) = gVtmp(i,j) - gV(i,j,k,bi,bj)
208     ENDDO
209     ENDDO
210     ENDIF
211 baylor 1.40 #ifdef ALLOW_DIAGNOSTICS
212 jmc 1.41 IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
213     DO j=1,sNy+1
214     DO i=1,sNx+1
215 baylor 1.40 gUext(i,j) = gU(i,j,k,bi,bj)
216     gVext(i,j) = gV(i,j,k,bi,bj)
217     ENDDO
218     ENDDO
219     ENDIF
220     #endif /* ALLOW_DIAGNOSTICS */
221    
222 jmc 1.32 CALL EXTERNAL_FORCING_U(
223     I iMin,iMax,jMin,jMax,bi,bj,k,
224     I myTime,myThid)
225     CALL EXTERNAL_FORCING_V(
226     I iMin,iMax,jMin,jMax,bi,bj,k,
227     I myTime,myThid)
228 jmc 1.37
229 baylor 1.40 #ifdef ALLOW_DIAGNOSTICS
230 jmc 1.41 IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
231     DO j=1,sNy+1
232     DO i=1,sNx+1
233 baylor 1.40 gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
234     gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
235     ENDDO
236     ENDDO
237     ENDIF
238     #endif /* ALLOW_DIAGNOSTICS */
239    
240 jmc 1.37 C- for CD-scheme, compute gU,Vtmp = gU,V^n + forcing
241 jmc 1.32 IF (useCDscheme) THEN
242 jmc 1.28 DO j=jMin,jMax
243     DO i=iMin,iMax
244 jmc 1.37 gUtmp(i,j) = gUtmp(i,j) + gU(i,j,k,bi,bj)
245     gVtmp(i,j) = gVtmp(i,j) + gV(i,j,k,bi,bj)
246 jmc 1.28 ENDDO
247 jmc 1.26 ENDDO
248 jmc 1.28 ENDIF
249 jmc 1.26 ENDIF
250    
251 edhill 1.34 #ifdef ALLOW_CD_CODE
252 jmc 1.32 IF (useCDscheme) THEN
253     C- Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
254     C and return coriolis terms on C-grid (guCor,gvCor)
255 heimbach 1.35 CALL CD_CODE_SCHEME(
256 edhill 1.34 I bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
257     O guCor,gvCor,
258     I myTime, myIter, myThid)
259 jmc 1.32 DO j=jMin,jMax
260     DO i=iMin,iMax
261     gUtmp(i,j) = gU(i,j,k,bi,bj)
262     & + guCor(i,j)
263     gVtmp(i,j) = gV(i,j,k,bi,bj)
264     & + gvCor(i,j)
265     ENDDO
266 cnh 1.1 ENDDO
267 jmc 1.37 ELSE
268     #endif /* ALLOW_CD_CODE */
269     DO j=jMin,jMax
270     DO i=iMin,iMax
271     gUtmp(i,j) = gU(i,j,k,bi,bj)
272     gVtmp(i,j) = gV(i,j,k,bi,bj)
273     ENDDO
274     ENDDO
275     #ifdef ALLOW_CD_CODE
276 jmc 1.32 ENDIF
277 jmc 1.37 #endif
278 adcroft 1.19
279 jmc 1.26 #ifdef NONLIN_FRSURF
280     IF (.NOT. vectorInvariantMomentum
281     & .AND. nonlinFreeSurf.GT.1) THEN
282 jmc 1.28 IF (select_rStar.GT.0) THEN
283     DO j=jMin,jMax
284     DO i=iMin,iMax
285 jmc 1.32 gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
286 jmc 1.28 gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
287     ENDDO
288     ENDDO
289     ELSE
290     DO j=jMin,jMax
291     DO i=iMin,iMax
292 jmc 1.32 IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
293     gUtmp(i,j) = gUtmp(i,j)
294 heimbach 1.46 & *_hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
295 jmc 1.32 ENDIF
296 jmc 1.28 IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
297 jmc 1.26 gVtmp(i,j) = gVtmp(i,j)
298 heimbach 1.46 & *_hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
299 jmc 1.28 ENDIF
300     ENDDO
301 jmc 1.26 ENDDO
302 jmc 1.28 ENDIF
303 jmc 1.26 ENDIF
304 jmc 1.32 #endif /* NONLIN_FRSURF */
305    
306     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
307    
308 jmc 1.36 C-- Dissipation term outside the Adams-Bashforth:
309 jmc 1.38 IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
310 jmc 1.36 DO j=jMin,jMax
311     DO i=iMin,iMax
312     gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
313     gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
314     ENDDO
315     ENDDO
316     ENDIF
317    
318 jmc 1.32 C Step forward zonal velocity (store in Gu)
319     DO j=jMin,jMax
320     DO i=iMin,iMax
321     gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
322     & +deltaTmom*(
323     & gUtmp(i,j)
324     & - psFac*phiSurfX(i,j)
325     & - phxFac*dPhiHydX(i,j)
326     & )*_maskW(i,j,k,bi,bj)
327     ENDDO
328     ENDDO
329 jmc 1.26
330 adcroft 1.7 C Step forward meridional velocity (store in Gv)
331 adcroft 1.19 DO j=jMin,jMax
332 cnh 1.1 DO i=iMin,iMax
333 jmc 1.32 gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
334 jmc 1.26 & +deltaTmom*(
335     & gVtmp(i,j)
336     & - psFac*phiSurfY(i,j)
337 jmc 1.29 & - phyFac*dPhiHydY(i,j)
338 adcroft 1.7 & )*_maskS(i,j,k,bi,bj)
339 cnh 1.1 ENDDO
340 adcroft 1.19 ENDDO
341 cnh 1.1
342 jmc 1.39 #ifdef ALLOW_DIAGNOSTICS
343 jmc 1.41 IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
344     CALL DIAGNOSTICS_FILL(gUext,'Um_Ext ',k,1,2,bi,bj,myThid)
345     CALL DIAGNOSTICS_FILL(gVext,'Vm_Ext ',k,1,2,bi,bj,myThid)
346     ENDIF
347 jmc 1.39 #ifdef ALLOW_CD_CODE
348     IF ( useCDscheme .AND. useDiagnostics ) THEN
349     CALL DIAGNOSTICS_FILL(guCor,'Um_Cori ',k,1,2,bi,bj,myThid)
350     CALL DIAGNOSTICS_FILL(gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid)
351     ENDIF
352     #endif /* ALLOW_CD_CODE */
353     #endif /* ALLOW_DIAGNOSTICS */
354    
355 cnh 1.1 RETURN
356     END
  ViewVC Help
Powered by ViewVC 1.1.22