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

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

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


Revision 1.70 - (show annotations) (download)
Mon Nov 3 22:40:22 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
Changes since 1.69: +3 -3 lines
 o convert all the '#ifndef "EXCLUDE_MONITOR"'-isms to the now-standard
   '#ifdef "ALLOW_MONITOR"' format as is done with all the other packages
   - note that this also means removal of the monitor package from the
     packages.conf list for the AD examples

1 C $Header: /u/u3/gcmpack/MITgcm/model/src/forward_step.F,v 1.69 2003/11/01 04:50:02 edhill Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 cswdptr -- add --
8 #ifdef ALLOW_GCHEM
9 # include "GCHEM_OPTIONS.h"
10 #endif
11 cswdptr -- end add ---
12
13 CBOP
14 C !ROUTINE: FORWARD_STEP
15 C !INTERFACE:
16 SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid )
17
18 C !DESCRIPTION: \bv
19 C *==================================================================
20 C | SUBROUTINE forward_step
21 C | o Run the ocean model and, optionally, evaluate a cost function.
22 C *==================================================================
23 C |
24 C | THE_MAIN_LOOP is the toplevel routine for the Tangent Linear and
25 C | Adjoint Model Compiler (TAMC). For this purpose the initialization
26 C | of the model was split into two parts. Those parameters that do
27 C | not depend on a specific model run are set in INITIALISE_FIXED,
28 C | whereas those that do depend on the specific realization are
29 C | initialized in INITIALISE_VARIA.
30 C |
31 C *==================================================================
32 C \ev
33
34 C !USES:
35 IMPLICIT NONE
36 C == Global variables ==
37 #include "SIZE.h"
38 #include "EEPARAMS.h"
39 #include "PARAMS.h"
40 #include "DYNVARS.h"
41 #include "FFIELDS.h"
42
43 #ifdef ALLOW_NONHYDROSTATIC
44 #include "CG3D.h"
45 #endif
46
47 #ifdef ALLOW_SHAP_FILT
48 #include "SHAP_FILT.h"
49 #endif
50 #ifdef ALLOW_ZONAL_FILT
51 #include "ZONAL_FILT.h"
52 #endif
53
54 #ifdef ALLOW_AUTODIFF_TAMC
55 # include "tamc.h"
56 # include "ctrl.h"
57 # include "ctrl_dummy.h"
58 # include "cost.h"
59 # include "EOS.h"
60 # ifdef ALLOW_EXF
61 # include "exf_fields.h"
62 # ifdef ALLOW_BULKFORMULAE
63 # include "exf_constants.h"
64 # endif
65 # endif
66 # ifdef ALLOW_OBCS
67 # include "OBCS.h"
68 # endif
69 # ifdef ALLOW_PTRACERS
70 # include "PTRACERS.h"
71 # endif
72 #endif /* ALLOW_AUTODIFF_TAMC */
73
74 C !LOCAL VARIABLES:
75 C == Routine arguments ==
76 C note: under the multi-threaded model myiter and
77 C mytime are local variables passed around as routine
78 C arguments. Although this is fiddly it saves the need to
79 C impose additional synchronisation points when they are
80 C updated.
81 C myiter - iteration counter for this thread
82 C mytime - time counter for this thread
83 C mythid - thread number for this instance of the routine.
84 integer iloop
85 integer mythid
86 integer myiter
87 _RL mytime
88 #ifdef ALLOW_BULK_FORCE
89 INTEGER bi,bj
90 #endif
91
92 CEOP
93
94 #ifndef DISABLE_DEBUGMODE
95 IF ( debugLevel .GE. debLevB )
96 & CALL DEBUG_ENTER('FORWARD_STEP',myThid)
97 #endif
98
99 #ifdef ALLOW_AUTODIFF_TAMC
100 C-- Reset the model iteration counter and the model time.
101 myiter = nIter0 + (iloop-1)
102 mytime = startTime + float(iloop-1)*deltaTclock
103 #endif
104
105 #if (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_AUTODIFF_MONITOR))
106 C Include call to a dummy routine. Its adjoint will be
107 C called at the proper place in the adjoint code.
108 C The adjoint routine will print out adjoint values
109 C if requested. The location of the call is important,
110 C it has to be after the adjoint of the exchanges
111 C (DO_GTERM_BLOCKING_EXCHANGES).
112 CALL DUMMY_IN_STEPPING( myTime, myIter, myThid )
113 cph I've commented this line since it may conflict with MITgcm's adjoint
114 cph However, need to check whether that's still consistent
115 cph with the ecco-branch (it should).
116 cph CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
117 #endif
118
119 #ifdef EXACT_CONSERV
120 IF (exactConserv) THEN
121 C-- Update etaH(n+1) :
122 CALL TIMER_START('UPDATE_ETAH [FORWARD_STEP]',mythid)
123 CALL UPDATE_ETAH( myTime, myIter, myThid )
124 CALL TIMER_STOP ('UPDATE_ETAH [FORWARD_STEP]',mythid)
125 ENDIF
126 #endif /* EXACT_CONSERV */
127
128 #ifdef NONLIN_FRSURF
129 IF ( select_rStar.NE.0 ) THEN
130 C-- r* : compute the future level thickness according to etaH(n+1)
131 CALL TIMER_START('CALC_R_STAR [FORWARD_STEP]',mythid)
132 CALL CALC_R_STAR(etaH, myTime, myIter, myThid )
133 CALL TIMER_STOP ('CALC_R_STAR [FORWARD_STEP]',mythid)
134 ELSEIF ( nonlinFreeSurf.GT.0) THEN
135 C-- compute the future surface level thickness according to etaH(n+1)
136 CALL TIMER_START('CALC_SURF_DR [FORWARD_STEP]',mythid)
137 CALL CALC_SURF_DR(etaH, myTime, myIter, myThid )
138 CALL TIMER_STOP ('CALC_SURF_DR [FORWARD_STEP]',mythid)
139 ENDIF
140 #endif /* NONLIN_FRSURF */
141
142 C-- Load forcing/external data fields.
143 #ifdef ALLOW_AUTODIFF_TAMC
144 c**************************************
145 #include "checkpoint_lev1_directives.h"
146 c**************************************
147 #endif
148
149
150 C-- Call external forcing package
151 #ifdef ALLOW_BULK_FORCE
152 IF ( useBulkforce ) THEN
153 #ifndef DISABLE_DEBUGMODE
154 IF ( debugLevel .GE. debLevB )
155 & CALL DEBUG_CALL('BULKF_FIELDS_LOAD',myThid)
156 #endif
157 CALL TIMER_START('BULKF_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)
158 CALL BULKF_FIELDS_LOAD( mytime, myiter, mythid )
159 CALL TIMER_STOP ('BULKF_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)
160 c calculate qnet and empmr (and wind stress)
161 DO bj=myByLo(myThid),myByHi(myThid)
162 DO bi=myBxLo(myThid),myBxHi(myThid)
163 CALL BULKF_FORCING( bi,bj, mytime, myiter, mythid )
164 ENDDO
165 ENDDO
166 c Update the tile edges.
167 _EXCH_XY_R8(Qnet, mythid)
168 _EXCH_XY_R8(EmPmR, mythid)
169 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
170 C _EXCH_XY_R8(fu , mythid)
171 C _EXCH_XY_R8(fv , mythid)
172 ELSE
173 #endif /* ALLOW_BULK_FORCE */
174
175 # ifdef ALLOW_EXF
176 C NOTE, that although the exf package is part of the
177 C distribution, it is not currently maintained, i.e.
178 C exf is disabled by default in genmake.
179 #ifndef DISABLE_DEBUGMODE
180 IF ( debugLevel .GE. debLevB )
181 & CALL DEBUG_CALL('EXF_GETFORCING',myThid)
182 #endif
183 CALL TIMER_START('EXF_GETFORCING [FORWARD_STEP]',mythid)
184 CALL EXF_GETFORCING( mytime, myiter, mythid )
185 CALL TIMER_STOP ('EXF_GETFORCING [FORWARD_STEP]',mythid)
186 # else /* ALLOW_EXF undef */
187 cph The following IF-statement creates an additional dependency
188 cph for the forcing fields requiring additional storing.
189 cph Therefore, the IF-statement will be put between CPP-OPTIONS,
190 cph assuming that ALLOW_SEAICE has not yet been differentiated.
191 # ifdef ALLOW_SEAICE
192 IF ( .NOT. useSEAICE ) THEN
193 # endif
194 #ifndef DISABLE_DEBUGMODE
195 IF ( debugLevel .GE. debLevB )
196 & CALL DEBUG_CALL('EXTERNAL_FIELDS_LOAD',myThid)
197 #endif
198 CALL TIMER_START('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
199 CALL EXTERNAL_FIELDS_LOAD( mytime, myiter, mythid )
200 CALL TIMER_STOP ('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
201 # ifdef ALLOW_SEAICE
202 ENDIF
203 # endif
204 # endif /* ALLOW_EXF */
205 #ifdef ALLOW_BULK_FORCE
206 C-- end of if/else block useBulfforce --
207 ENDIF
208 #endif /* ALLOW_BULK_FORCE */
209
210 #ifdef ALLOW_AUTODIFF
211 c-- Add control vector for forcing and parameter fields
212 if ( myiter .EQ. nIter0 )
213 & CALL CTRL_MAP_FORCING (mythid)
214 #endif
215
216 # ifdef ALLOW_SEAICE
217 C-- Call sea ice model to compute forcing/external data fields. In
218 C addition to computing prognostic sea-ice variables and diagnosing the
219 C forcing/external data fields that drive the ocean model, SEAICE_MODEL
220 C also sets theta to the freezing point under sea-ice. The implied
221 C surface heat flux is then stored in variable surfaceTendencyTice,
222 C which is needed by KPP package (kpp_calc.F and kpp_transport_t.F)
223 C to diagnose surface buoyancy fluxes and for the non-local transport
224 C term. Because this call precedes model thermodynamics, temperature
225 C under sea-ice may not be "exactly" at the freezing point by the time
226 C theta is dumped or time-averaged.
227 IF ( useSEAICE ) THEN
228 #ifndef DISABLE_DEBUGMODE
229 IF ( debugLevel .GE. debLevB )
230 & CALL DEBUG_CALL('SEAICE_MODEL',myThid)
231 #endif
232 CALL TIMER_START('SEAICE_MODEL [FORWARD_STEP]',myThid)
233 CALL SEAICE_MODEL( myTime, myIter, myThid )
234 CALL TIMER_STOP ('SEAICE_MODEL [FORWARD_STEP]',myThid)
235 ENDIF
236 # endif /* ALLOW_SEAICE */
237
238 #ifdef ALLOW_AUTODIFF_TAMC
239 # ifdef ALLOW_PTRACERS
240 cph this replaces _bibj storing of ptracer within thermodynamics
241 CADJ STORE ptracer = comlev1, key = ikey_dynamics
242 # endif
243 #endif
244
245 #ifdef ALLOW_PTRACERS
246 # ifdef ALLOW_GCHEM
247 CALL GCHEM_FIELDS_LOAD( mytime, myiter, mythid )
248 # endif
249 #endif
250
251 C-- Step forward fields and calculate time tendency terms.
252 #ifndef DISABLE_DEBUGMODE
253 IF ( debugLevel .GE. debLevB )
254 & CALL DEBUG_CALL('THERMODYNAMICS',myThid)
255 #endif
256 CALL TIMER_START('THERMODYNAMICS [FORWARD_STEP]',mythid)
257 CALL THERMODYNAMICS( myTime, myIter, myThid )
258 CALL TIMER_STOP ('THERMODYNAMICS [FORWARD_STEP]',mythid)
259
260 C-- do exchanges (needed for DYNAMICS) when using stagger time-step :
261 #ifndef DISABLE_DEBUGMODE
262 IF ( debugLevel .GE. debLevB )
263 & CALL DEBUG_CALL('DO_STAGGER_FIELDS_EXCH.',myThid)
264 #endif
265 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
266 CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
267 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
268
269 #ifdef ALLOW_SHAP_FILT
270 IF (useSHAP_FILT .AND.
271 & staggerTimeStep .AND. shap_filt_TrStagg ) THEN
272 #ifndef DISABLE_DEBUGMODE
273 IF ( debugLevel .GE. debLevB )
274 & CALL DEBUG_CALL('SHAP_FILT_APPLY_TS',myThid)
275 #endif
276 CALL TIMER_START('SHAP_FILT [FORWARD_STEP]',myThid)
277 CALL SHAP_FILT_APPLY_TS(gT,gS,myTime+deltaT,myIter+1,myThid)
278 CALL TIMER_STOP ('SHAP_FILT [FORWARD_STEP]',myThid)
279 ENDIF
280 #endif
281 #ifdef ALLOW_ZONAL_FILT
282 IF (useZONAL_FILT .AND.
283 & staggerTimeStep .AND. zonal_filt_TrStagg ) THEN
284 #ifndef DISABLE_DEBUGMODE
285 IF ( debugLevel .GE. debLevB )
286 & CALL DEBUG_CALL('ZONAL_FILT_APPLY_TS',myThid)
287 #endif
288 CALL TIMER_START('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
289 CALL ZONAL_FILT_APPLY_TS( gT, gS, myThid )
290 CALL TIMER_STOP ('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
291 ENDIF
292 #endif
293
294 C-- Step forward fields and calculate time tendency terms.
295 #ifndef ALLOW_AUTODIFF_TAMC
296 IF ( momStepping ) THEN
297 #endif
298 #ifndef DISABLE_DEBUGMODE
299 IF ( debugLevel .GE. debLevB )
300 & CALL DEBUG_CALL('DYNAMICS',myThid)
301 #endif
302 CALL TIMER_START('DYNAMICS [FORWARD_STEP]',mythid)
303 CALL DYNAMICS( myTime, myIter, myThid )
304 CALL TIMER_STOP ('DYNAMICS [FORWARD_STEP]',mythid)
305 #ifndef ALLOW_AUTODIFF_TAMC
306 ENDIF
307 #endif
308
309 #ifdef ALLOW_NONHYDROSTATIC
310 C-- Step forward W field in N-H algorithm
311 IF ( momStepping .AND. nonHydrostatic ) THEN
312 #ifndef DISABLE_DEBUGMODE
313 IF ( debugLevel .GE. debLevB )
314 & CALL DEBUG_CALL('CALC_GW',myThid)
315 #endif
316 CALL TIMER_START('CALC_GW [FORWARD_STEP]',myThid)
317 CALL CALC_GW(myThid)
318 CALL TIMER_STOP ('CALC_GW [FORWARD_STEP]',myThid)
319 ENDIF
320 #endif
321
322 #ifdef NONLIN_FRSURF
323 C-- update hfacC,W,S and recip_hFac according to etaH(n+1) :
324 IF ( nonlinFreeSurf.GT.0) THEN
325 IF ( select_rStar.GT.0 ) THEN
326 CALL TIMER_START('UPDATE_R_STAR [FORWARD_STEP]',myThid)
327 CALL UPDATE_R_STAR( myTime, myIter, myThid )
328 CALL TIMER_STOP ('UPDATE_R_STAR [FORWARD_STEP]',myThid)
329 ELSE
330 CALL TIMER_START('UPDATE_SURF_DR [FORWARD_STEP]',myThid)
331 CALL UPDATE_SURF_DR( myTime, myIter, myThid )
332 CALL TIMER_STOP ('UPDATE_SURF_DR [FORWARD_STEP]',myThid)
333 ENDIF
334 ENDIF
335 C- update also CG2D matrix (and preconditioner)
336 IF ( momStepping .AND. nonlinFreeSurf.GT.2 ) THEN
337 CALL TIMER_START('UPDATE_CG2D [FORWARD_STEP]',myThid)
338 CALL UPDATE_CG2D( myTime, myIter, myThid )
339 CALL TIMER_STOP ('UPDATE_CG2D [FORWARD_STEP]',myThid)
340 ENDIF
341 #endif
342
343 C-- Apply Filters to u*,v* before SOLVE_FOR_PRESSURE
344 #ifdef ALLOW_SHAP_FILT
345 IF (useSHAP_FILT .AND. shap_filt_uvStar) THEN
346 CALL TIMER_START('SHAP_FILT [FORWARD_STEP]',myThid)
347 IF (implicDiv2Dflow.LT.1.) THEN
348 C-- Explicit+Implicit part of the Barotropic Flow Divergence
349 C => Filtering of uVel,vVel is necessary
350 CALL SHAP_FILT_APPLY_UV( uVel,vVel,
351 & myTime+deltaT, myIter+1, myThid )
352 ENDIF
353 CALL SHAP_FILT_APPLY_UV( gU,gV,myTime+deltaT,myIter+1,myThid)
354 CALL TIMER_STOP ('SHAP_FILT [FORWARD_STEP]',myThid)
355 ENDIF
356 #endif
357 #ifdef ALLOW_ZONAL_FILT
358 IF (useZONAL_FILT .AND. zonal_filt_uvStar) THEN
359 CALL TIMER_START('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
360 IF (implicDiv2Dflow.LT.1.) THEN
361 C-- Explicit+Implicit part of the Barotropic Flow Divergence
362 C => Filtering of uVel,vVel is necessary
363 CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
364 ENDIF
365 CALL ZONAL_FILT_APPLY_UV( gU, gV, myThid )
366 CALL TIMER_STOP ('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
367 ENDIF
368 #endif
369
370 C-- Solve elliptic equation(s).
371 C Two-dimensional only for conventional hydrostatic or
372 C three-dimensional for non-hydrostatic and/or IGW scheme.
373 IF ( momStepping ) THEN
374 CALL TIMER_START('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
375 CALL SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
376 CALL TIMER_STOP ('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
377 ENDIF
378
379 #ifdef ALLOW_AUTODIFF_TAMC
380 cph This is needed because convective_adjustment calls
381 cph find_rho which may use pressure()
382 CADJ STORE totphihyd = comlev1, key = ikey_dynamics
383 #endif
384 C-- Correct divergence in flow field and cycle time-stepping
385 C arrays (for all fields) ; update time-counter
386 myIter = nIter0 + iLoop
387 myTime = startTime + deltaTClock * float(iLoop)
388 CALL TIMER_START('THE_CORRECTION_STEP [FORWARD_STEP]',myThid)
389 CALL THE_CORRECTION_STEP(myTime, myIter, myThid)
390 CALL TIMER_STOP ('THE_CORRECTION_STEP [FORWARD_STEP]',myThid)
391
392 C-- Do "blocking" sends and receives for tendency "overlap" terms
393 c CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
394 c CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )
395 c CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
396
397 C-- Do "blocking" sends and receives for field "overlap" terms
398 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
399 CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
400 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
401
402 cswdptr -- add for seperate timestepping of chemical/biological/forcing
403 cswdptr of ptracers ---
404 #ifdef ALLOW_GCHEM
405 ceh3 This is broken -- this ifdef should not be visible!
406 #ifdef PTRACERS_SEPARATE_FORCING
407 ceh3 needs an IF ( use GCHEM ) THEN
408 call GCHEM_FORCING_SEP( myTime,myIter,myThid )
409 #endif /* PTRACERS_SEPARATE_FORCING */
410 #endif /* ALLOW_GCHEM */
411 cswdptr -- end add ---
412
413
414 #ifdef ALLOW_FLT
415 C-- Calculate float trajectories
416 IF (useFLT) THEN
417 CALL TIMER_START('FLOATS [FORWARD_STEP]',myThid)
418 CALL FLT_MAIN(myIter,myTime, myThid)
419 CALL TIMER_STOP ('FLOATS [FORWARD_STEP]',myThid)
420 ENDIF
421 #endif
422
423 #ifdef ALLOW_MONITOR
424 C-- Check status of solution (statistics, cfl, etc...)
425 CALL TIMER_START('MONITOR [FORWARD_STEP]',myThid)
426 CALL MONITOR( myIter, myTime, myThid )
427 CALL TIMER_STOP ('MONITOR [FORWARD_STEP]',myThid)
428 #endif /* ALLOW_MONITOR */
429
430 C-- Do IO if needed.
431 CALL TIMER_START('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
432 CALL DO_THE_MODEL_IO( myTime, myIter, myThid )
433 CALL TIMER_STOP ('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
434
435 C-- Save state for restarts
436 C Note: (jmc: is it still the case after ckp35 ?)
437 C =====
438 C Because of the ordering of the timestepping code and
439 C tendency term code at end of loop model arrays hold
440 C U,V,T,S at "time-level" N but gu, gv, gs, gt, guNM1,...
441 C at "time-level" N+1/2 (guNM1 at "time-level" N+1/2 is
442 C gu at "time-level" N-1/2) and etaN at "time-level" N+1/2.
443 C where N = I+timeLevBase-1
444 C Thus a checkpoint contains U.0000000000, GU.0000000001 and
445 C etaN.0000000001 in the indexing scheme used for the model
446 C "state" files. This example is referred to as a checkpoint
447 C at time level 1
448 CALL TIMER_START('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
449 CALL WRITE_CHECKPOINT(
450 & .FALSE., myTime, myIter, myThid )
451 CALL TIMER_STOP ('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
452
453 #ifndef DISABLE_DEBUGMODE
454 IF ( debugLevel .GE. debLevB )
455 & CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
456 #endif
457
458 END

  ViewVC Help
Powered by ViewVC 1.1.22