/[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.138 - (show annotations) (download)
Tue Apr 4 14:52:43 2006 UTC (18 years, 2 months ago) by heimbach
Branch: MAIN
Changes since 1.137: +3 -3 lines 
Fix ptracers adjoint
o avoid extensive recomputatations
o fix missing re-init. of gptr (missed by TAF)
1 C $Header: /u/gcmpack/MITgcm/model/src/forward_step.F,v 1.137 2006/03/30 02:33:05 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 #ifdef ALLOW_GMREDI
8 # include "GMREDI_OPTIONS.h"
9 #endif
10 #ifdef ALLOW_OBCS
11 # include "OBCS_OPTIONS.h"
12 #endif
13
14 CBOP
15 C !ROUTINE: FORWARD_STEP
16 C !INTERFACE:
17 SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid )
18
19 C !DESCRIPTION: \bv
20 C *==================================================================
21 C | SUBROUTINE forward_step
22 C | o Run the ocean model and, optionally, evaluate a cost function.
23 C *==================================================================
24 C |
25 C | THE_MAIN_LOOP is the toplevel routine for the Tangent Linear and
26 C | Adjoint Model Compiler (TAMC). For this purpose the initialization
27 C | of the model was split into two parts. Those parameters that do
28 C | not depend on a specific model run are set in INITIALISE_FIXED,
29 C | whereas those that do depend on the specific realization are
30 C | initialized in INITIALISE_VARIA.
31 C |
32 C *==================================================================
33 C \ev
34
35 C !USES:
36 IMPLICIT NONE
37 C == Global variables ==
38 #include "SIZE.h"
39 #include "EEPARAMS.h"
40 #include "PARAMS.h"
41 #include "DYNVARS.h"
42
43 #ifdef ALLOW_MNC
44 #include "MNC_PARAMS.h"
45 EXTERNAL DIFFERENT_MULTIPLE
46 LOGICAL DIFFERENT_MULTIPLE
47 #endif
48
49 #ifdef HAVE_SIGREG
50 #include "SIGREG.h"
51 #endif
52
53 #ifdef ALLOW_SHAP_FILT
54 # include "SHAP_FILT.h"
55 #endif
56 #ifdef ALLOW_ZONAL_FILT
57 # include "ZONAL_FILT.h"
58 #endif
59 #ifdef COMPONENT_MODULE
60 # include "CPL_PARAMS.h"
61 #endif
62
63 #ifdef ALLOW_AUTODIFF_TAMC
64 # include "FFIELDS.h"
65
66 # include "tamc.h"
67 # include "ctrl.h"
68 # include "ctrl_dummy.h"
69 # include "cost.h"
70 # include "EOS.h"
71 # ifdef NONLIN_FRSURF
72 # include "GRID.h"
73 # endif
74 # ifdef ALLOW_EXF
75 # include "exf_fields.h"
76 # include "exf_clim_fields.h"
77 # ifdef ALLOW_BULKFORMULAE
78 # include "exf_constants.h"
79 # endif
80 # endif
81 # ifdef ALLOW_PTRACERS
82 # include "PTRACERS_SIZE.h"
83 # include "PTRACERS.h"
84 # endif
85 # ifdef ALLOW_OBCS
86 # include "OBCS.h"
87 # ifdef ALLOW_PTRACERS
88 # include "OBCS_PTRACERS.h"
89 # endif
90 # endif
91 # ifdef ALLOW_CD_CODE
92 # include "CD_CODE_VARS.h"
93 # endif
94 # ifdef ALLOW_EBM
95 # include "EBM.h"
96 # endif
97 # ifdef EXACT_CONSERV
98 # include "SURFACE.h"
99 # endif
100 # ifdef ALLOW_KPP
101 # include "KPP.h"
102 # endif
103 # ifdef ALLOW_GMREDI
104 # include "GMREDI.h"
105 # endif
106 #endif /* ALLOW_AUTODIFF_TAMC */
107
108 C !LOCAL VARIABLES:
109 C == Routine arguments ==
110 C note: under the multi-threaded model myiter and
111 C mytime are local variables passed around as routine
112 C arguments. Although this is fiddly it saves the need to
113 C impose additional synchronisation points when they are
114 C updated.
115 C myIter - iteration counter for this thread
116 C myTime - time counter for this thread
117 C myThid - thread number for this instance of the routine.
118 INTEGER iloop
119 INTEGER myThid
120 INTEGER myIter
121 _RL myTime
122
123 C == Local variables ==
124 #ifdef COMPONENT_MODULE
125 INTEGER myItP1
126 #endif
127 CEOP
128
129 #ifdef ALLOW_DEBUG
130 IF ( debugLevel .GE. debLevB )
131 & CALL DEBUG_ENTER('FORWARD_STEP',myThid)
132 #endif
133
134 #ifdef ALLOW_AUTODIFF_TAMC
135 C-- Reset the model iteration counter and the model time.
136 myIter = nIter0 + (iloop-1)
137 myTime = startTime + float(iloop-1)*deltaTclock
138 #endif
139
140 #ifdef ALLOW_AUTODIFF_TAMC
141 c**************************************
142 #include "checkpoint_lev1_directives.h"
143 c**************************************
144 #endif
145
146 C-- Switch on/off diagnostics for snap-shot output:
147 #ifdef ALLOW_DIAGNOSTICS
148 IF ( useDiagnostics ) THEN
149 CALL DIAGNOSTICS_SWITCH_ONOFF( myTime, myIter, myThid )
150 C-- State-variables diagnostics
151 CALL TIMER_START('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid)
152 CALL DO_STATEVARS_DIAGS( myTime, 0, myIter, myThid )
153 CALL TIMER_STOP ('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid)
154 ENDIF
155 #endif
156
157 C-- Call Bulk-Formulae forcing package
158 #ifdef ALLOW_BULK_FORCE
159 IF ( useBulkForce ) THEN
160 #ifdef ALLOW_DEBUG
161 IF ( debugLevel .GE. debLevB )
162 & CALL DEBUG_CALL('BULKF_FIELDS_LOAD',myThid)
163 #endif
164 CALL TIMER_START('BULKF_FORCING [FORWARD_STEP]',mythid)
165 C- load all forcing fields at current time
166 CALL BULKF_FIELDS_LOAD( myTime, myIter, myThid )
167 C- calculate qnet and empmr (and wind stress)
168 CALL BULKF_FORCING( myTime, myIter, myThid )
169 CALL TIMER_STOP ('BULKF_FORCING [FORWARD_STEP]',mythid)
170 ENDIF
171 #endif /* ALLOW_BULK_FORCE */
172
173 C-- Call external forcing package
174 # ifdef ALLOW_EXF
175 # ifdef ALLOW_DEBUG
176 IF ( debugLevel .GE. debLevB )
177 & CALL DEBUG_CALL('EXF_GETFORCING',myThid)
178 # endif
179 CALL TIMER_START('EXF_GETFORCING [FORWARD_STEP]',mythid)
180 CALL EXF_GETFORCING( mytime, myiter, mythid )
181 CALL TIMER_STOP ('EXF_GETFORCING [FORWARD_STEP]',mythid)
182 # else /* ALLOW_EXF undef */
183 cph The following IF-statement creates an additional dependency
184 cph for the forcing fields requiring additional storing.
185 cph Therefore, the IF-statement will be put between CPP-OPTIONS,
186 cph assuming that ALLOW_SEAICE has not yet been differentiated.
187 # if (defined (ALLOW_SEAICE) || defined (ALLOW_EBM))
188 IF ( .NOT. useSEAICE .AND. .NOT. useEBM ) THEN
189 # endif
190 #ifdef ALLOW_DEBUG
191 IF ( debugLevel .GE. debLevB )
192 & CALL DEBUG_CALL('EXTERNAL_FIELDS_LOAD',myThid)
193 #endif
194 CALL TIMER_START('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
195 CALL EXTERNAL_FIELDS_LOAD( mytime, myiter, mythid )
196 CALL TIMER_STOP ('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
197 # ifdef NONLIN_FRSURF
198 CADJ STORE SST = comlev1, key = ikey_dynamics
199 CADJ STORE SSS = comlev1, key = ikey_dynamics
200 # ifdef SHORTWAVE_HEATING
201 CADJ STORE Qsw = comlev1, key = ikey_dynamics
202 # endif
203 # endif
204 # if (defined (ALLOW_SEAICE) || defined (ALLOW_EBM))
205 ENDIF
206 # endif
207 # endif /* ALLOW_EXF */
208
209 #ifdef ALLOW_AUTODIFF
210 c-- Add control vector for forcing and parameter fields
211 if ( myiter .EQ. nIter0 )
212 & CALL CTRL_MAP_FORCING (mythid)
213 #endif
214
215 #if (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_AUTODIFF_MONITOR))
216 C Include call to a dummy routine. Its adjoint will be
217 C called at the proper place in the adjoint code.
218 C The adjoint routine will print out adjoint values
219 C if requested. The location of the call is important,
220 C it has to be after the adjoint of the exchanges
221 C (DO_GTERM_BLOCKING_EXCHANGES).
222 CALL DUMMY_IN_STEPPING( myTime, myIter, myThid )
223 cph I've commented this line since it may conflict with MITgcm's adjoint
224 cph However, need to check whether that's still consistent
225 cph with the ecco-branch (it should).
226 cph CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
227 #endif
228
229 # ifdef ALLOW_SEAICE
230 C-- Call sea ice model to compute forcing/external data fields. In
231 C addition to computing prognostic sea-ice variables and diagnosing the
232 C forcing/external data fields that drive the ocean model, SEAICE_MODEL
233 C also sets theta to the freezing point under sea-ice. The implied
234 C surface heat flux is then stored in variable surfaceTendencyTice,
235 C which is needed by KPP package (kpp_calc.F and kpp_transport_t.F)
236 C to diagnose surface buoyancy fluxes and for the non-local transport
237 C term. Because this call precedes model thermodynamics, temperature
238 C under sea-ice may not be "exactly" at the freezing point by the time
239 C theta is dumped or time-averaged.
240 IF ( useSEAICE ) THEN
241 #ifdef ALLOW_DEBUG
242 IF ( debugLevel .GE. debLevB )
243 & CALL DEBUG_CALL('SEAICE_MODEL',myThid)
244 #endif
245 CALL TIMER_START('SEAICE_MODEL [FORWARD_STEP]',myThid)
246 CALL SEAICE_MODEL( myTime, myIter, myThid )
247 CALL TIMER_STOP ('SEAICE_MODEL [FORWARD_STEP]',myThid)
248 ENDIF
249 # endif /* ALLOW_SEAICE */
250
251 #ifdef ALLOW_AUTODIFF_TAMC
252 # ifdef ALLOW_PTRACERS
253 cph this is now also replaced by storing elsewhere
254 cphCADJ STORE ptracer = comlev1, key = ikey_dynamics
255 # endif
256 #endif
257
258 #ifdef ALLOW_OFFLINE
259 IF ( useOffLine ) THEN
260 CALL TIMER_START('OFFLINE_FIELDS_LOAD [FORWARD_STEP]',myThid)
261 #ifdef ALLOW_DEBUG
262 IF ( debugLevel .GE. debLevB )
263 & CALL DEBUG_CALL('OFFLINE_FIELDS_LOAD',myThid)
264 #endif /* ALLOW_DEBUG */
265 CALL OFFLINE_FIELDS_LOAD( myTime, myIter, myThid )
266 CALL TIMER_STOP ('OFFLINE_FIELDS_LOAD [FORWARD_STEP]',myThid)
267 ENDIF
268 #endif /* ALLOW_OFFLINE */
269
270 #ifdef ALLOW_GCHEM
271 IF ( useGCHEM ) THEN
272 #ifdef ALLOW_DEBUG
273 IF ( debugLevel .GE. debLevB )
274 & CALL DEBUG_CALL('GCHEM_FIELDS_LOAD',myThid)
275 #endif /* ALLOW_DEBUG */
276 CALL GCHEM_FIELDS_LOAD( mytime, myiter, mythid )
277 ENDIF
278 #endif /* ALLOW_GCHEM */
279
280 #ifdef ALLOW_RBCS
281 IF ( useRBCS ) THEN
282 CALL RBCS_FIELDS_LOAD( mytime, myiter, mythid )
283 ENDIF
284 #endif
285
286 #ifdef COMPONENT_MODULE
287 IF ( useCoupler .AND. cpl_earlyExpImpCall ) THEN
288 C Post coupling data that I export.
289 C Read in coupling data that I import.
290 CALL TIMER_START('CPL_EXPORT-IMPORT [FORWARD_STEP]',myThid)
291 CALL CPL_EXPORT_MY_DATA( myIter, myTime, myThid )
292 CALL CPL_IMPORT_EXTERNAL_DATA( myIter, myTime, myThid )
293 CALL TIMER_STOP ('CPL_EXPORT-IMPORT [FORWARD_STEP]',myThid)
294 ENDIF
295 #endif /* COMPONENT_MODULE */
296
297 #ifdef ALLOW_EBM
298 IF ( useEBM ) THEN
299 # ifdef ALLOW_DEBUG
300 IF ( debugLevel .GE. debLevB )
301 & CALL DEBUG_CALL('EBM',myThid)
302 # endif
303 CALL TIMER_START('EBM [FORWARD_STEP]',mythid)
304 CALL EBM_DRIVER ( myTime, myIter, myThid )
305 CALL TIMER_STOP ('EBM [FORWARD_STEP]',mythid)
306 ENDIF
307 #endif
308
309 C-- Step forward fields and calculate time tendency terms.
310
311 #ifdef ALLOW_DEBUG
312 IF ( debugLevel .GE. debLevB )
313 & CALL DEBUG_CALL('DO_ATMOSPHERIC_PHYS',myThid)
314 #endif
315 CALL TIMER_START('DO_ATMOSPHERIC_PHYS [FORWARD_STEP]',mythid)
316 CALL DO_ATMOSPHERIC_PHYS( myTime, myIter, myThid )
317 CALL TIMER_STOP ('DO_ATMOSPHERIC_PHYS [FORWARD_STEP]',mythid)
318
319 #ifdef ALLOW_AUTODIFF_TAMC
320 # ifdef ALLOW_KPP
321 CADJ STORE uvel = comlev1, key = ikey_dynamics
322 CADJ STORE vvel = comlev1, key = ikey_dynamics
323 # endif
324 # ifdef EXACT_CONSERV
325 cphCADJ STORE empmr = comlev1, key = ikey_dynamics
326 cphCADJ STORE pmepr = comlev1, key = ikey_dynamics
327 # endif
328 # ifdef NONLIN_FRSURF
329 cph-test
330 CADJ STORE hFacC = comlev1, key = ikey_dynamics
331 # endif
332 #endif /* ALLOW_AUTODIFF_TAMC */
333
334 #ifndef ALLOW_AUTODIFF_TAMC
335 IF ( .NOT. useOffLine ) THEN
336 #endif
337 #ifdef ALLOW_DEBUG
338 IF ( debugLevel .GE. debLevB )
339 & CALL DEBUG_CALL('DO_OCEANIC_PHYS',myThid)
340 #endif
341 CALL TIMER_START('DO_OCEANIC_PHYS [FORWARD_STEP]',mythid)
342 CALL DO_OCEANIC_PHYS( myTime, myIter, myThid )
343 CALL TIMER_STOP ('DO_OCEANIC_PHYS [FORWARD_STEP]',mythid)
344 #ifdef ALLOW_AUTODIFF_TAMC
345 CADJ STORE EmPmR = comlev1, key = ikey_dynamics
346 # ifdef EXACT_CONSERV
347 CADJ STORE pmepr = comlev1, key = ikey_dynamics
348 # endif
349 #else
350 ENDIF
351 #endif
352
353 #ifdef ALLOW_AUTODIFF_TAMC
354 # ifdef NONLIN_FRSURF
355 cph-test
356 CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics
357 CADJ STORE hfac_surfs = comlev1, key = ikey_dynamics
358 CADJ STORE hfac_surfw = comlev1, key = ikey_dynamics
359 CADJ STORE hFacC, hFacS, hFacW
360 CADJ & = comlev1, key = ikey_dynamics
361 CADJ STORE recip_hFacC, recip_hFacS, recip_hFacW
362 CADJ & = comlev1, key = ikey_dynamics
363 c
364 CADJ STORE surfaceforcingu = comlev1, key = ikey_dynamics
365 CADJ STORE surfaceforcingv = comlev1, key = ikey_dynamics
366 # endif
367 #endif /* ALLOW_AUTODIFF_TAMC */
368
369 #ifdef ALLOW_GCHEM
370 C GCHEM package is an interface for any bio-geochemical or
371 C ecosystem model you would like to include.
372 C If GCHEM_SEPARATE_FORCING is not defined, you are
373 C responsible for computing tendency terms for passive
374 C tracers and storing them on a 3DxNumPtracers-array called
375 C gchemTendency in GCHEM_CALC_TENDENCY. This tendency is then added
376 C to gPtr in ptracers_forcing later-on.
377 C If GCHEM_SEPARATE_FORCING is defined, you are reponsible for
378 C UPDATING ptracers directly in GCHEM_FORCING_SEP. This amounts
379 C to a completely separate time step that you have to implement
380 C yourself (Eulerian seems to be fine in most cases).
381 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
382 C CAVEAT: Up to now, when GCHEM is turned on the field ptracerForcingSurf,
383 C which is needed for KPP is not set properly. ptracerForcingSurf must
384 C be treated differently depending on whether GCHEM_SEPARATE_FORCING
385 C is define or not. TBD.
386 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
387 IF ( useGCHEM ) THEN
388 #ifdef ALLOW_DEBUG
389 IF ( debugLevel .GE. debLevB )
390 & CALL DEBUG_CALL('GCHEM_CALC_TENDENCY',myThid)
391 #endif
392 CALL TIMER_START('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid)
393 CALL GCHEM_CALC_TENDENCY( myTime, myIter, myThid )
394 CALL TIMER_STOP ('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid)
395 ENDIF
396 #endif /* ALLOW_GCHEM */
397
398 #ifdef ALLOW_AUTODIFF_TAMC
399 cph needed to be moved here from do_oceanic_physics
400 cph to be visible down the road
401 c
402 CADJ STORE surfaceForcingS = comlev1, key = ikey_dynamics
403 CADJ STORE surfaceForcingT = comlev1, key = ikey_dynamics
404 CADJ STORE surfaceForcingTice = comlev1, key = ikey_dynamics
405 ctest(
406 CADJ STORE IVDConvCount = comlev1, key = ikey_dynamics
407 ctest)
408 # ifdef ALLOW_PTRACERS
409 CADJ STORE surfaceForcingPtr = comlev1, key = ikey_dynamics
410 # endif
411 c
412 # ifdef ALLOW_GMREDI
413 CADJ STORE Kwx = comlev1, key = ikey_dynamics
414 CADJ STORE Kwy = comlev1, key = ikey_dynamics
415 CADJ STORE Kwz = comlev1, key = ikey_dynamics
416 # ifdef GM_BOLUS_ADVEC
417 CADJ STORE GM_PsiX = comlev1, key = ikey_dynamics
418 CADJ STORE GM_PsiY = comlev1, key = ikey_dynamics
419 # endif
420 # endif
421 c
422 # ifdef ALLOW_KPP
423 CADJ STORE KPPghat = comlev1, key = ikey_dynamics
424 CADJ STORE KPPfrac = comlev1, key = ikey_dynamics
425 CADJ STORE KPPdiffKzS = comlev1, key = ikey_dynamics
426 CADJ STORE KPPdiffKzT = comlev1, key = ikey_dynamics
427 # endif
428 c
429 # ifdef NONLIN_FRSURF
430 CADJ STORE etaH = comlev1, key = ikey_dynamics
431 # ifdef ALLOW_CD_CODE
432 CADJ STORE etanm1 = comlev1, key = ikey_dynamics
433 # endif
434 # endif
435 #endif /* ALLOW_AUTODIFF_TAMC */
436
437 IF ( .NOT.staggerTimeStep ) THEN
438 #ifdef ALLOW_DEBUG
439 IF ( debugLevel .GE. debLevB )
440 & CALL DEBUG_CALL('THERMODYNAMICS',myThid)
441 #endif
442 CALL TIMER_START('THERMODYNAMICS [FORWARD_STEP]',mythid)
443 CALL THERMODYNAMICS( myTime, myIter, myThid )
444 CALL TIMER_STOP ('THERMODYNAMICS [FORWARD_STEP]',mythid)
445 C-- if not staggerTimeStep: end
446 ENDIF
447 c #ifdef ALLOW_NONHYDROSTATIC
448 IF ( implicitIntGravWave ) THEN
449 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
450 CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
451 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
452 ENDIF
453 c #endif
454
455 #ifdef COMPONENT_MODULE
456 IF ( useCoupler .AND. .NOT.cpl_earlyExpImpCall ) THEN
457 C Post coupling data that I export.
458 C Read in coupling data that I import.
459 myItP1 = myIter + 1
460 CALL TIMER_START('CPL_EXPORT-IMPORT [FORWARD_STEP]',myThid)
461 CALL CPL_EXPORT_MY_DATA( myItP1, myTime, myThid )
462 CALL CPL_IMPORT_EXTERNAL_DATA( myItP1, myTime, myThid )
463 CALL TIMER_STOP ('CPL_EXPORT-IMPORT [FORWARD_STEP]',myThid)
464 # ifdef ALLOW_OCN_COMPON_INTERF
465 IF ( useRealFreshWaterFlux ) THEN
466 CALL OCN_APPLY_IMPORT( .FALSE., myTime, myIter, myThid )
467 ENDIF
468 # endif /* ALLOW_OCN_COMPON_INTERF */
469 ENDIF
470 #endif /* COMPONENT_MODULE */
471
472 #ifdef ALLOW_AUTODIFF_TAMC
473 # ifdef NONLIN_FRSURF
474 CADJ STORE hFacC = comlev1, key = ikey_dynamics
475 CADJ STORE hFacS = comlev1, key = ikey_dynamics
476 CADJ STORE hFacW = comlev1, key = ikey_dynamics
477 CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics
478 CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics
479 CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics
480 CADJ STORE etaN = comlev1, key = ikey_dynamics
481 # endif
482 #endif
483 C-- Step forward fields and calculate time tendency terms.
484 #ifndef ALLOW_AUTODIFF_TAMC
485 IF ( momStepping ) THEN
486 #endif
487 #ifdef ALLOW_DEBUG
488 IF ( debugLevel .GE. debLevB )
489 & CALL DEBUG_CALL('DYNAMICS',myThid)
490 #endif
491 CALL TIMER_START('DYNAMICS [FORWARD_STEP]',mythid)
492 CALL DYNAMICS( myTime, myIter, myThid )
493 CALL TIMER_STOP ('DYNAMICS [FORWARD_STEP]',mythid)
494 #ifndef ALLOW_AUTODIFF_TAMC
495 ENDIF
496 #endif
497
498 #ifdef ALLOW_AUTODIFF_TAMC
499 # ifdef NONLIN_FRSURF
500 cph-test
501 CADJ STORE gU, gV = comlev1, key = ikey_dynamics
502 # endif
503 #endif
504
505 C-- Update time-counter
506 myIter = nIter0 + iLoop
507 myTime = startTime + deltaTClock * float(iLoop)
508
509 #ifdef ALLOW_MNC
510 C Update the default next iter for MNC
511 IF ( useMNC ) THEN
512 CALL MNC_CW_CITER_SETG( 1, 1, -1, myIter , myThid )
513
514 C TODO: Logic should be added here so that users can specify, on
515 C a per-citer-group basis, when it is time to update the
516 C "current" (and not just the "next") iteration
517
518 C TODO: the following is just a temporary band-aid (mostly, for
519 C Baylor) until someone writes a routine that better handles time
520 C boundaries such as weeks, months, years, etc.
521 IF ( mnc_filefreq .GT. 0 ) THEN
522 IF (DIFFERENT_MULTIPLE(mnc_filefreq,myTime,deltaTClock))
523 & THEN
524 CALL MNC_CW_CITER_SETG( 1, 1, myIter, -1 , myThid )
525 ENDIF
526 ENDIF
527 ENDIF
528 #endif
529
530 C-- Update geometric factors:
531 #ifdef NONLIN_FRSURF
532 C- update hfacC,W,S and recip_hFac according to etaH(n+1) :
533 IF ( nonlinFreeSurf.GT.0) THEN
534 IF ( select_rStar.GT.0 ) THEN
535 # ifndef DISABLE_RSTAR_CODE
536 # ifdef ALLOW_AUTODIFF_TAMC
537 cph-test
538 CADJ STORE hFacC = comlev1, key = ikey_dynamics
539 CADJ STORE hFacS = comlev1, key = ikey_dynamics
540 CADJ STORE hFacW = comlev1, key = ikey_dynamics
541 CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics
542 CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics
543 CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics
544 # endif
545 CALL TIMER_START('UPDATE_R_STAR [FORWARD_STEP]',myThid)
546 CALL UPDATE_R_STAR( myTime, myIter, myThid )
547 CALL TIMER_STOP ('UPDATE_R_STAR [FORWARD_STEP]',myThid)
548 # ifdef ALLOW_AUTODIFF_TAMC
549 cph-test
550 CADJ STORE hFacC = comlev1, key = ikey_dynamics
551 CADJ STORE hFacS = comlev1, key = ikey_dynamics
552 CADJ STORE hFacW = comlev1, key = ikey_dynamics
553 CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics
554 CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics
555 CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics
556 # endif
557 # endif /* DISABLE_RSTAR_CODE */
558 ELSE
559 #ifdef ALLOW_AUTODIFF_TAMC
560 CADJ STORE hFac_surfC, hFac_surfS, hFac_surfW
561 CADJ & = comlev1, key = ikey_dynamics
562 #endif
563 CALL TIMER_START('UPDATE_SURF_DR [FORWARD_STEP]',myThid)
564 CALL UPDATE_SURF_DR( myTime, myIter, myThid )
565 CALL TIMER_STOP ('UPDATE_SURF_DR [FORWARD_STEP]',myThid)
566 ENDIF
567 ENDIF
568 # ifdef ALLOW_AUTODIFF_TAMC
569 cph-test
570 CADJ STORE hFacC = comlev1, key = ikey_dynamics
571 CADJ STORE hFacS = comlev1, key = ikey_dynamics
572 CADJ STORE hFacW = comlev1, key = ikey_dynamics
573 CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics
574 CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics
575 CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics
576 # endif
577 C- update also CG2D matrix (and preconditioner)
578 IF ( momStepping .AND. nonlinFreeSurf.GT.2 ) THEN
579 CALL TIMER_START('UPDATE_CG2D [FORWARD_STEP]',myThid)
580 CALL UPDATE_CG2D( myTime, myIter, myThid )
581 CALL TIMER_STOP ('UPDATE_CG2D [FORWARD_STEP]',myThid)
582 ENDIF
583 #endif /* NONLIN_FRSURF */
584
585 C-- Apply Filters to u*,v* before SOLVE_FOR_PRESSURE
586 #ifdef ALLOW_SHAP_FILT
587 IF (useSHAP_FILT .AND. shap_filt_uvStar) THEN
588 CALL TIMER_START('SHAP_FILT_UV [FORWARD_STEP]',myThid)
589 IF (implicDiv2Dflow.LT.1.) THEN
590 C-- Explicit+Implicit part of the Barotropic Flow Divergence
591 C => Filtering of uVel,vVel is necessary
592 CALL SHAP_FILT_APPLY_UV( uVel,vVel,
593 & myTime, myIter, myThid )
594 ENDIF
595 CALL SHAP_FILT_APPLY_UV( gU,gV,myTime,myIter,myThid)
596 CALL TIMER_STOP ('SHAP_FILT_UV [FORWARD_STEP]',myThid)
597 ENDIF
598 #endif
599 #ifdef ALLOW_ZONAL_FILT
600 IF (useZONAL_FILT .AND. zonal_filt_uvStar) THEN
601 CALL TIMER_START('ZONAL_FILT_UV [FORWARD_STEP]',myThid)
602 IF (implicDiv2Dflow.LT.1.) THEN
603 C-- Explicit+Implicit part of the Barotropic Flow Divergence
604 C => Filtering of uVel,vVel is necessary
605 CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
606 ENDIF
607 CALL ZONAL_FILT_APPLY_UV( gU, gV, myThid )
608 CALL TIMER_STOP ('ZONAL_FILT_UV [FORWARD_STEP]',myThid)
609 ENDIF
610 #endif
611
612 C-- Solve elliptic equation(s).
613 C Two-dimensional only for conventional hydrostatic or
614 C three-dimensional for non-hydrostatic and/or IGW scheme.
615 IF ( momStepping ) THEN
616 #ifdef ALLOW_AUTODIFF_TAMC
617 # ifdef NONLIN_FRSURF
618 CADJ STORE uvel, vvel
619 CADJ & = comlev1, key = ikey_dynamics
620 CADJ STORE empmr,hfacs,hfacw
621 CADJ & = comlev1, key = ikey_dynamics
622 # endif
623 #endif
624 CALL TIMER_START('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
625 CALL SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
626 CALL TIMER_STOP ('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
627 ENDIF
628
629 C-- Correct divergence in flow field and cycle time-stepping momentum
630 #ifndef ALLOW_AUTODIFF_TAMC
631 IF ( momStepping ) THEN
632 #endif
633 CALL TIMER_START('MOM_CORRECTION_STEP [FORWARD_STEP]',myThid)
634 CALL MOMENTUM_CORRECTION_STEP(myTime, myIter, myThid)
635 CALL TIMER_STOP ('MOM_CORRECTION_STEP [FORWARD_STEP]',myThid)
636 #ifndef ALLOW_AUTODIFF_TAMC
637 ENDIF
638 #endif
639
640 #ifdef EXACT_CONSERV
641 IF (exactConserv) THEN
642 #ifdef ALLOW_AUTODIFF_TAMC
643 cph-test
644 cphCADJ STORE etaH = comlev1, key = ikey_dynamics
645 #endif
646 C-- Update etaH(n+1) :
647 CALL TIMER_START('UPDATE_ETAH [FORWARD_STEP]',mythid)
648 CALL UPDATE_ETAH( myTime, myIter, myThid )
649 CALL TIMER_STOP ('UPDATE_ETAH [FORWARD_STEP]',mythid)
650 ENDIF
651 #endif /* EXACT_CONSERV */
652
653 #ifdef NONLIN_FRSURF
654 IF ( select_rStar.NE.0 ) THEN
655 # ifndef DISABLE_RSTAR_CODE
656 C-- r* : compute the future level thickness according to etaH(n+1)
657 CALL TIMER_START('CALC_R_STAR [FORWARD_STEP]',mythid)
658 CALL CALC_R_STAR(etaH, myTime, myIter, myThid )
659 CALL TIMER_STOP ('CALC_R_STAR [FORWARD_STEP]',mythid)
660 # endif /* DISABLE_RSTAR_CODE */
661 ELSEIF ( nonlinFreeSurf.GT.0) THEN
662 C-- compute the future surface level thickness according to etaH(n+1)
663 # ifdef ALLOW_AUTODIFF_TAMC
664 CADJ STORE etaH = comlev1, key = ikey_dynamics
665 # endif
666 CALL TIMER_START('CALC_SURF_DR [FORWARD_STEP]',mythid)
667 CALL CALC_SURF_DR(etaH, myTime, myIter, myThid )
668 CALL TIMER_STOP ('CALC_SURF_DR [FORWARD_STEP]',mythid)
669 ENDIF
670 # ifdef ALLOW_AUTODIFF_TAMC
671 cph-test
672 CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics
673 # endif
674 #endif /* NONLIN_FRSURF */
675
676 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
677 IF ( staggerTimeStep ) THEN
678 C-- do exchanges of U,V (needed for multiDim) when using stagger time-step :
679 #ifdef ALLOW_DEBUG
680 IF ( debugLevel .GE. debLevB )
681 & CALL DEBUG_CALL('DO_STAGGER_FIELDS_EXCH.',myThid)
682 #endif
683 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
684 CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
685 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
686
687 #ifdef ALLOW_DIAGNOSTICS
688 C-- State-variables diagnostics
689 IF ( useDiagnostics ) THEN
690 CALL TIMER_START('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid)
691 CALL DO_STATEVARS_DIAGS( myTime, 1, myIter, myThid )
692 CALL TIMER_STOP ('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid)
693 ENDIF
694 #endif
695
696 #ifdef ALLOW_DEBUG
697 IF ( debugLevel .GE. debLevB )
698 & CALL DEBUG_CALL('THERMODYNAMICS',myThid)
699 #endif
700 CALL TIMER_START('THERMODYNAMICS [FORWARD_STEP]',mythid)
701 CALL THERMODYNAMICS( myTime, myIter, myThid )
702 CALL TIMER_STOP ('THERMODYNAMICS [FORWARD_STEP]',mythid)
703
704 C-- if staggerTimeStep: end
705 ENDIF
706 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
707
708 #ifdef ALLOW_AUTODIFF_TAMC
709 cph This is needed because convective_adjustment calls
710 cph find_rho which may use pressure()
711 CADJ STORE totphihyd = comlev1, key = ikey_dynamics
712 #endif
713 C-- Cycle time-stepping Tracers arrays (T,S,+pTracers)
714 CALL TIMER_START('TRC_CORRECTION_STEP [FORWARD_STEP]',myThid)
715 CALL TRACERS_CORRECTION_STEP(myTime, myIter, myThid)
716 CALL TIMER_STOP ('TRC_CORRECTION_STEP [FORWARD_STEP]',myThid)
717
718 #ifdef ALLOW_GCHEM
719 C Add separate timestepping of chemical/biological/forcing
720 C of ptracers here in GCHEM_FORCING_SEP
721 IF ( useGCHEM ) THEN
722 #ifdef ALLOW_DEBUG
723 IF ( debugLevel .GE. debLevB )
724 & CALL DEBUG_CALL('GCHEM_FORCING_SEP',myThid)
725 #endif /* ALLOW_DEBUG */
726 CALL TIMER_START('GCHEM_FORCING_SEP [FORWARD_STEP]',myThid)
727 CALL GCHEM_FORCING_SEP( myTime,myIter,myThid )
728 CALL TIMER_STOP ('GCHEM_FORCING_SEP [FORWARD_STEP]',myThid)
729 ENDIF
730 #endif /* ALLOW_GCHEM */
731
732 C-- Do "blocking" sends and receives for tendency "overlap" terms
733 c CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
734 c CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )
735 c CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
736
737 C-- Do "blocking" sends and receives for field "overlap" terms
738 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
739 CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
740 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
741
742 #ifdef ALLOW_DIAGNOSTICS
743 IF ( useDiagnostics ) THEN
744 CALL TIMER_START('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid)
745 CALL DO_STATEVARS_DIAGS( myTime, 2, myIter, myThid )
746 CALL TIMER_STOP ('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid)
747 ENDIF
748 #endif
749
750 C AMM
751 #ifdef ALLOW_GRIDALT
752 if (useGRIDALT) then
753 CALL GRIDALT_UPDATE(myThid)
754 endif
755 #endif
756 C AMM
757
758 C AMM
759 #ifdef ALLOW_FIZHI
760 if( useFIZHI) then
761 CALL TIMER_START('FIZHI [FORWARD_STEP]',mythid)
762 CALL STEP_FIZHI_CORR ( myTime, myIter, myThid, dTtracerLev(1) )
763 CALL TIMER_STOP('FIZHI [FORWARD_STEP]',mythid)
764 endif
765 #endif
766 C AMM
767
768 #ifdef ALLOW_FLT
769 C-- Calculate float trajectories
770 IF (useFLT) THEN
771 CALL TIMER_START('FLOATS [FORWARD_STEP]',myThid)
772 CALL FLT_MAIN(myIter,myTime, myThid)
773 CALL TIMER_STOP ('FLOATS [FORWARD_STEP]',myThid)
774 ENDIF
775 #endif
776
777 #ifdef ALLOW_TIMEAVE
778 C-- State-variables time-averaging
779 CALL TIMER_START('DO_STATEVARS_TAVE [FORWARD_STEP]',myThid)
780 CALL DO_STATEVARS_TAVE( myTime, myIter, myThid )
781 CALL TIMER_STOP ('DO_STATEVARS_TAVE [FORWARD_STEP]',myThid)
782 #endif
783
784 #ifdef ALLOW_MONITOR
785 IF ( .NOT.useOffLine ) THEN
786 C-- Check status of solution (statistics, cfl, etc...)
787 CALL TIMER_START('MONITOR [FORWARD_STEP]',myThid)
788 CALL MONITOR( myIter, myTime, myThid )
789 CALL TIMER_STOP ('MONITOR [FORWARD_STEP]',myThid)
790 ENDIF
791 #endif /* ALLOW_MONITOR */
792
793 #ifdef ALLOW_COST
794 C-- compare model with data and compute cost function
795 C-- this is done after exchanges to allow interpolation
796 CALL TIMER_START('COST_TILE [FORWARD_STEP]',myThid)
797 CALL COST_TILE ( mytime, myiter, myThid )
798 CALL TIMER_STOP ('COST_TILE [FORWARD_STEP]',myThid)
799 #endif
800
801 C-- Do IO if needed.
802 CALL TIMER_START('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
803 CALL DO_THE_MODEL_IO( myTime, myIter, myThid )
804 CALL TIMER_STOP ('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
805
806 #ifdef HAVE_SIGREG
807 IF ( useSIGREG ) THEN
808 IF ( i_got_signal .GT. 0 ) THEN
809 CALL PACKAGES_WRITE_PICKUP(
810 I .TRUE., myTime, myIter, myThid )
811 CALL WRITE_CHECKPOINT(
812 I .TRUE., myTime, myIter, myThid )
813 STOP 'Checkpoint completed -- killed by signal handler'
814 ENDIF
815 ENDIF
816 #endif /* HAVE_SIGREG */
817
818 C-- Save state for restarts
819 CALL TIMER_START('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
820 CALL PACKAGES_WRITE_PICKUP(
821 I .FALSE., myTime, myIter, myThid )
822 CALL WRITE_CHECKPOINT(
823 I .FALSE., myTime, myIter, myThid )
824 CALL TIMER_STOP ('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
825
826 #ifdef ALLOW_DEBUG
827 IF ( debugLevel .GE. debLevB )
828 & CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
829 #endif
830
831 RETURN
832 END
  ViewVC Help
Powered by ViewVC 1.1.22