/[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.140 - (show annotations) (download)
Fri Apr 28 22:53:14 2006 UTC (18 years, 1 month ago) by heimbach
Branch: MAIN
Changes since 1.139: +2 -1 lines
o SEAICE_CGRID adjoint, part 2.
  (all stores seem to be sorted out, but NANs in adjoint).

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

  ViewVC Help
Powered by ViewVC 1.1.22