/[MITgcm]/MITgcm/pkg/ecco/the_main_loop.F
ViewVC logotype

Contents of /MITgcm/pkg/ecco/the_main_loop.F

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


Revision 1.2 - (show annotations) (download)
Thu Nov 13 19:35:10 2003 UTC (20 years, 7 months ago) by heimbach
Branch: MAIN
Changes since 1.1: +2 -2 lines
added missing cost_averages_bar_directives.h
(plus name change)

1 C $Header: /u/gcmpack/MITgcm/pkg/ecco/the_main_loop.F,v 1.1 2003/11/06 22:10:08 heimbach Exp $
2
3 #include "PACKAGES_CONFIG.h"
4 #include "CPP_OPTIONS.h"
5
6 #ifdef ALLOW_OBCS
7 # include "OBCS_OPTIONS.h"
8 #endif
9 #ifdef ALLOW_SEAICE
10 # include "SEAICE_OPTIONS.h"
11 #endif
12
13 subroutine the_main_loop( myTime, myIter, mythid )
14
15 c ==================================================================
16 c SUBROUTINE the_main_loop
17 c ==================================================================
18 c
19 c o Run the ocean model and evaluate the specified cost function.
20 c
21 c *the_main_loop* is the top-level routine for the Tangent Linear and
22 c Adjoint Model Compiler (TAMC). For this purpose, the initialization
23 c of the model was split into two parts. Those parameters that do
24 c not depend on a specific model run are set in *initialise_fixed*,
25 c whereas those that do depend on the specific realization are
26 c initialized in *initialise_varia*. In order to do a so called
27 c checkpointing during the adjoint calculation and to account for the
28 c typical data involved in oceanographic applications a call tree
29 c that is divided into yearly, monthly, daily, and step parts can
30 c be used.
31 c
32 c This routine is to be used in conjuction with the MITgcmuv release
33 c checkpoint 24.
34 c
35 c started: Christian Eckert eckert@mit.edu 30-Jun-1999
36 c
37 c changed: Christian Eckert eckert@mit.edu 14-Jul-1999
38 c
39 c - The call to mapping was moved to initialise_varia,
40 c since this routine has to be called before
41 c ini_predictor.
42 c
43 c Christian Eckert eckert@mit.edu 11-Feb-2000
44 c
45 c - Restructured the code in order to create a package
46 c for the MITgcmUV.
47 c
48 c Patrick Heimbach heimbach@mit.edu 3-Jun-2000
49 c - corrected computation of ikey_dynamics and
50 c added computation of ikey_dynamics for the case
51 c undef ALLOW_TAMC_CHECKPOINTING
52 c
53 c Patrick Heimbach heimbach@mit.edu 6-Jun-2000
54 c - corrected initialisation of comlev1 common blocks
55 c
56 c Dimitris Menemenlis menemenlis@jpl.nasa.gov 26-Feb-2003
57 c - modifications for pkg/seaice
58 c
59 c ==================================================================
60 c SUBROUTINE the_main_loop
61 c ==================================================================
62
63 implicit none
64
65 c == global variables ==
66
67 #include "SIZE.h"
68 #include "EEPARAMS.h"
69 #include "PARAMS.h"
70
71 c**************************************
72 #ifdef ALLOW_AUTODIFF_TAMC
73
74 c These includes are needed for
75 c AD-checkpointing.
76 c They provide the fields to be stored.
77
78 # include "GRID.h"
79 # include "DYNVARS.h"
80 # include "FFIELDS.h"
81 # include "EOS.h"
82 # include "GAD.h"
83
84 # ifdef ALLOW_CD_CODE
85 # include "CD_CODE_VARS.h"
86 # endif
87 # ifdef ALLOW_PASSIVE_TRACER
88 # include "TR1.h"
89 # endif
90 # ifdef ALLOW_PTRACERS
91 # include "PTRACERS.h"
92 # endif
93 # ifdef ALLOW_NONHYDROSTATIC
94 # include "CG3D.h"
95 # endif
96 # ifdef EXACT_CONSERV
97 # include "SURFACE.h"
98 # endif
99 # ifdef ALLOW_OBCS
100 # include "OBCS.h"
101 # endif
102 # ifdef ALLOW_EXF
103 # include "exf_fields.h"
104 # include "exf_clim_fields.h"
105 # ifdef ALLOW_BULKFORMULAE
106 # include "exf_constants.h"
107 # endif
108 # endif /* ALLOW_EXF */
109 # ifdef ALLOW_SEAICE
110 # include "SEAICE.h"
111 # endif
112 # ifdef ALLOW_DIVIDED_ADJOINT_MPI
113 # include "mpif.h"
114 # endif
115
116 # include "tamc.h"
117 # include "ctrl.h"
118 # include "ctrl_dummy.h"
119 # include "cost.h"
120
121 #endif /* ALLOW_AUTODIFF_TAMC */
122 c**************************************
123
124 c == routine arguments ==
125 c note: under the multi-threaded model myiter and
126 c mytime are local variables passed around as routine
127 c arguments. Although this is fiddly it saves the need to
128 c impose additional synchronisation points when they are
129 c updated.
130 c myiter - iteration counter for this thread
131 c mytime - time counter for this thread
132 c mythid - thread number for this instance of the routine.
133 integer mythid
134 integer myiter
135 _RL mytime
136
137 c == local variables ==
138
139 integer bi,bj
140 integer iloop
141 integer mydate(4)
142 #ifdef ALLOW_SNAPSHOTS
143 character yprefix*3
144 #endif
145
146 #ifdef ALLOW_TAMC_CHECKPOINTING
147 integer ilev_1
148 integer ilev_2
149 integer ilev_3
150 integer max_lev2
151 integer max_lev3
152 #endif
153
154 c-- == end of interface ==
155
156 #ifndef DISABLE_DEBUGMODE
157 IF ( debugLevel .GE. debLevB )
158 & CALL DEBUG_ENTER('THE_MAIN_LOOP',myThid)
159 #endif
160
161 #ifdef ALLOW_AUTODIFF_TAMC
162 c-- Initialize storage for the initialisations.
163 CADJ INIT tapelev3 = USER
164 c-- Some more initialisations to please TAMC
165 CADJ INIT tapelev_ini_bibj_k = USER
166 # ifdef ALLOW_DIVIDED_ADJOINT
167 CADJ INIT onetape = user
168 cphCADJ INIT onetape = common, 1
169 cph We want to avoid common blocks except in the inner loop.
170 cph Reason: the active write and consecutive read may occur
171 cph in separate model executions for which the info
172 cph in common blocks are lost.
173 cph Thus, we can only store real values (no integers)
174 cph because we only have active file handling to real available.
175 # endif
176 # ifdef ALLOW_TAMC_CHECKPOINTING
177 nIter0 = INT( startTime/deltaTClock )
178 ikey_dynamics = 1
179 # endif
180 #endif /* ALLOW_AUTODIFF_TAMC */
181
182 myTime = startTime
183 myIter = nIter0
184
185 CALL TIMER_START('ECCO SPIN-UP', mythid)
186
187 c-- Get the current date.
188 call CAL_TIMESTAMP( myiter, mytime, mydate, mythid )
189
190 C-- Set initial conditions (variable arrays)
191 #ifndef DISABLE_DEBUGMODE
192 IF ( debugLevel .GE. debLevB )
193 & CALL DEBUG_CALL('INITIALISE_VARIA',myThid)
194 #endif
195 CALL TIMER_START('INITIALISE_VARIA [THE_MAIN_LOOP]', mythid)
196 CALL INITIALISE_VARIA( mythid )
197 CALL TIMER_STOP ('INITIALISE_VARIA [THE_MAIN_LOOP]', mythid)
198
199 c-- Dump for start state.
200 #ifndef DISABLE_DEBUGMODE
201 IF ( debugLevel .GE. debLevB )
202 & CALL DEBUG_CALL('WRITE_STATE',myThid)
203 #endif
204 CALL TIMER_START('WRITE_STATE [THE_MAIN_LOOP]', mythid)
205 CALL WRITE_STATE( mytime, myiter, mythid )
206 CALL TIMER_STOP ('WRITE_STATE [THE_MAIN_LOOP]', mythid)
207
208 #ifdef ALLOW_MONITOR
209 C-- Check status of solution (statistics, cfl, etc...)
210 #ifndef DISABLE_DEBUGMODE
211 IF ( debugLevel .GE. debLevB )
212 & CALL DEBUG_CALL('MONITOR',myThid)
213 #endif
214 CALL TIMER_START('MONITOR [THE_MAIN_LOOP]', mythid)
215 CALL MONITOR( myIter, myTime, myThid )
216 CALL TIMER_STOP ('MONITOR [THE_MAIN_LOOP]', mythid)
217 #endif /* ALLOW_MONITOR */
218
219
220 #ifdef ALLOW_COST
221 c-- Compute the cost function contribution of the boundary forcing,
222 c-- i.e. heat flux, salt flux, zonal and meridional wind stress.
223 call timer_start('COST_FORCING [ECCO SPIN-UP]', mythid)
224 call cost_forcing( myiter, mytime, mythid )
225 call timer_stop ('COST_FORCING [ECCO SPIN-UP]', mythid)
226
227 # ifdef ALLOW_OBCS_COST_CONTRIBUTION
228 call timer_start('cost_obcs [ECCO SPIN-UP]', mythid)
229 call cost_obcs( myiter, mytime, mythid )
230 call timer_stop ('cost_obcs [ECCO SPIN-UP]', mythid)
231 # endif
232
233 #endif /* ALLOW_COST */
234
235 call timer_stop ('ECCO SPIN-UP', mythid)
236 _BARRIER
237
238 c-- Do the model integration.
239 call timer_start('ECCO MAIN LOOP',mythid)
240
241 c >>>>>>>>>>>>>>>>>>>>>>>>>>> LOOP <<<<<<<<<<<<<<<<<<<<<<<<<<<<
242 c >>>>>>>>>>>>>>>>>>>>>>>>>>> STARTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<
243
244 #ifdef ALLOW_AUTODIFF_TAMC
245 #ifdef ALLOW_TAMC_CHECKPOINTING
246 c-- Implement a three level checkpointing. For a two level
247 c-- checkpointing delete the middle loop; for n levels (n > 3)
248 c-- insert more loops.
249
250 c-- Check the choice of the checkpointing parameters in relation
251 c-- to nTimeSteps: (nchklev_1*nchklev_2*nchklev_3 .ge. nTimeSteps)
252 if (nchklev_1*nchklev_2*nchklev_3 .lt. nTimeSteps) then
253 print*
254 print*, ' the_main_loop: TAMC checkpointing parameters'
255 print*, ' nchklev_1*nchklev_2*nchklev_3 = ',
256 & nchklev_1*nchklev_2*nchklev_3
257 print*, ' are not consistent with nTimeSteps = ',
258 & nTimeSteps
259 stop ' ... stopped in the_main_loop.'
260 endif
261 max_lev3=nTimeSteps/(nchklev_1*nchklev_2)+1
262 max_lev2=nTimeSteps/nchklev_1+1
263
264 c**************************************
265 #ifdef ALLOW_DIVIDED_ADJOINT
266 CADJ loop = divided
267 #endif
268 c**************************************
269
270 do ilev_3 = 1,nchklev_3
271 if(ilev_3.le.max_lev3) then
272 c**************************************
273 #include "checkpoint_lev3_directives.h"
274 c**************************************
275
276 c-- Initialise storage for the middle loop.
277 CADJ INIT tapelev2 = USER
278
279 do ilev_2 = 1,nchklev_2
280 if(ilev_2.le.max_lev2) then
281 c**************************************
282 #include "checkpoint_lev2_directives.h"
283 c**************************************
284
285
286 c**************************************
287 #ifdef ALLOW_AUTODIFF_TAMC
288 c-- Initialize storage for the innermost loop.
289 c-- Always check common block sizes for the checkpointing!
290 c--
291 CADJ INIT comlev1 = COMMON,nchklev_1
292 CADJ INIT comlev1_bibj = COMMON,nchklev_1*nsx*nsy*nthreads_chkpt
293 CADJ INIT comlev1_bibj_k = COMMON,nchklev_1*nsx*nsy*nr*nthreads_chkpt
294 c--
295 # ifdef ALLOW_KPP
296 CADJ INIT comlev1_kpp = COMMON,nchklev_1*nsx*nsy
297 # endif /* ALLOW_KPP */
298 c--
299 # ifdef ALLOW_GMREDI
300 CADJ INIT comlev1_gmredi_k_gad
301 CADJ & = COMMON,nchklev_1*nsx*nsy*nr*nthreads_chkpt*maxpass
302 # endif /* ALLOW_GMREDI */
303 c--
304 # ifdef ALLOW_PTRACERS
305 CADJ INIT comlev1_bibj_ptracers = COMMON,
306 CADJ & nchklev_1*nsx*nsy*nthreads_chkpt*NUMBER_OF_PTRACERS
307 # endif /* ALLOW_PTRACERS */
308 c--
309 # ifndef DISABLE_MULTIDIM_ADVECTION
310 CADJ INIT comlev1_bibj_k_gad
311 CADJ & = COMMON,nchklev_1*nsx*nsy*nr*nthreads_chkpt*maxpass
312 CADJ INIT comlev1_bibj_k_gad_pass
313 CADJ & = COMMON,nchklev_1*nsx*nsy*nr*nthreads_chkpt*maxpass*maxcube
314 # endif /* DISABLE_MULTIDIM_ADVECTION */
315 c--
316 # if (defined (ALLOW_EXF) && defined (ALLOW_BULKFORMULAE))
317 CADJ INIT comlev1_exf_1
318 CADJ & = COMMON,nchklev_1*snx*nsx*sny*nsy*nthreads_chkpt
319 CADJ INIT comlev1_exf_2
320 CADJ & = COMMON,niter_bulk*nchklev_1*snx*nsx*sny*nsy*nthreads_chkpt
321 # endif
322 c--
323 # ifdef ALLOW_SEAICE
324 # ifdef SEAICE_ALLOW_DYNAMICS
325 CADJ INIT comlev1_lsr = COMMON,nchklev_1*2
326 # endif
327 # ifdef SEAICE_MULTILEVEL
328 CADJ INIT comlev1_multdim
329 CADJ & = COMMON,nchklev_1*nsx*nsy*nthreads_chkpt*multdim
330 # endif
331 # endif /* ALLOW_SEAICE */
332 c--
333 #endif /* ALLOW_AUTODIFF_TAMC */
334 c**************************************
335
336 do ilev_1 = 1,nchklev_1
337
338 c-- The if-statement below introduces a some flexibility in the
339 c-- choice of the 3-tupel ( nchklev_1, nchklev_2, nchklev_3 ).
340 c--
341 c-- Requirement: nchklev_1*nchklev_2*nchklev_3 .ge. nTimeSteps .
342
343 iloop = (ilev_3 - 1)*nchklev_2*nchklev_1 +
344 & (ilev_2 - 1)*nchklev_1 + ilev_1
345
346 if ( iloop .le. nTimeSteps ) then
347
348 #else /* ALLOW_TAMC_CHECKPOINTING undefined */
349 c-- Initialise storage for the reference trajectory without TAMC check-
350 c-- pointing.
351 CADJ INIT history = USER
352 CADJ INIT comlev1_bibj = COMMON,nchklev_0*nsx*nsy*nthreads_chkpt
353 CADJ INIT comlev1_bibj_k = COMMON,nchklev_0*nsx*nsy*nr*nthreads_chkpt
354 CADJ INIT comlev1_kpp = COMMON,nchklev_0*nsx*nsy
355
356 c-- Check the choice of the checkpointing parameters in relation
357 c-- to nTimeSteps: (nchklev_0 .ge. nTimeSteps)
358 if (nchklev_0 .lt. nTimeSteps) then
359 print*
360 print*, ' the_main_loop: ',
361 & 'TAMC checkpointing parameter nchklev_0 = ',
362 & nchklev_0
363 print*, ' is not consistent with nTimeSteps = ',
364 & nTimeSteps
365 stop ' ... stopped in the_main_loop.'
366 endif
367
368 do iloop = 1, nTimeSteps
369
370 #endif /* ALLOW_TAMC_CHECKPOINTING */
371
372 #else /* ALLOW_AUTODIFF_TAMC undefined */
373 c-- Start the main loop of ecco_Objfunc. Automatic differentiation is
374 c-- NOT enabled.
375 do iloop = 1, nTimeSteps
376 #endif /* ALLOW_AUTODIFF_TAMC */
377
378 #ifdef ALLOW_TAMC_CHECKPOINTING
379 nIter0 = INT( startTime/deltaTClock )
380 ikey_dynamics = ilev_1
381 #endif
382
383 c-- Set the model iteration counter and the model time.
384 myiter = nIter0 + (iloop-1)
385 mytime = startTime + float(iloop-1)*deltaTclock
386
387 #ifdef ALLOW_COST
388
389 c-- Accumulate time averages of temperature, salinity, and SSH.
390 call timer_start('COST_AVERAGESFIELDS [ECCO MAIN]', mythid)
391 call cost_averagesFields( mytime, mythid )
392 call timer_stop ('COST_AVERAGESFIELDS [ECCO MAIN]', mythid)
393 #ifdef ALLOW_COST_ATLANTIC
394 c-- Compute meridional heat transport
395 call timer_start('cost_atlantic [ECCO MAIN]', mythid)
396 call cost_atlantic( mytime, myiter,mythid )
397 call timer_stop ('cost_atlantic [ECCO MAIN]', mythid)
398 #endif
399 #endif /* ALLOW_COST */
400
401 #ifdef EXACT_CONSERV
402 IF (exactConserv) THEN
403 C-- Update etaH(n+1) :
404 CALL TIMER_START('UPDATE_ETAH [FORWARD_STEP]',mythid)
405 CALL UPDATE_ETAH( myTime, myIter, myThid )
406 CALL TIMER_STOP ('UPDATE_ETAH [FORWARD_STEP]',mythid)
407 ENDIF
408 #endif /* EXACT_CONSERV */
409
410 #ifdef NONLIN_FRSURF
411 IF ( select_rStar.NE.0 ) THEN
412 C-- r* : compute the future level thickness according to etaH(n+1)
413 CALL TIMER_START('CALC_R_STAR [FORWARD_STEP]',mythid)
414 CALL CALC_R_STAR(etaH, myTime, myIter, myThid )
415 CALL TIMER_STOP ('CALC_R_STAR [FORWARD_STEP]',mythid)
416 ELSEIF ( nonlinFreeSurf.GT.0) THEN
417 C-- compute the future surface level thickness according to etaH(n+1)
418 CALL TIMER_START('CALC_SURF_DR [FORWARD_STEP]',mythid)
419 CALL CALC_SURF_DR(etaH, myTime, myIter, myThid )
420 CALL TIMER_STOP ('CALC_SURF_DR [FORWARD_STEP]',mythid)
421 ENDIF
422 #endif /* NONLIN_FRSURF */
423
424 #ifdef ALLOW_AUTODIFF_TAMC
425 c**************************************
426 #include "checkpoint_lev1_directives.h"
427 c**************************************
428 #endif
429
430 #ifndef DISABLE_DEBUGMODE
431 IF ( debugLevel .GE. debLevB )
432 & CALL DEBUG_CALL('EXF_GETFORCING',myThid)
433 #endif
434 CALL TIMER_START('EXF_GETFORCING [FORWARD_STEP]',mythid)
435 CALL EXF_GETFORCING( mytime, myiter, mythid )
436 CALL TIMER_STOP ('EXF_GETFORCING [FORWARD_STEP]',mythid)
437
438 #ifdef ALLOW_SEAICE
439 C-- Call sea ice model to compute forcing/external data fields. In
440 C addition to computing prognostic sea-ice variables and diagnosing the
441 C forcing/external data fields that drive the ocean model, SEAICE_MODEL
442 C also sets theta to the freezing point under sea-ice. The implied
443 C surface heat flux is then stored in variable surfaceTendencyTice,
444 C which is needed by KPP package (kpp_calc.F and kpp_transport_t.F)
445 C to diagnose surface buoyancy fluxes and for the non-local transport
446 C term. Because this call precedes model thermodynamics, temperature
447 C under sea-ice may not be "exactly" at the freezing point by the time
448 C theta is dumped or time-averaged.
449 cph this simple runtime flag causes a lot of recomp.
450 cph IF ( useSEAICE ) THEN
451 #ifndef DISABLE_DEBUGMODE
452 IF ( debugLevel .GE. debLevB )
453 & CALL DEBUG_CALL('SEAICE_MODEL',myThid)
454 #endif
455 CALL TIMER_START('SEAICE_MODEL [FORWARD_STEP]',myThid)
456 CALL SEAICE_MODEL( myTime, myIter, myThid )
457 CALL TIMER_STOP ('SEAICE_MODEL [FORWARD_STEP]',myThid)
458 #ifdef ALLOW_COST_ICE
459 CALL COST_ICE ( myTime, myIter, myThid )
460 #endif
461 cph ENDIF
462 #endif /* ALLOW_SEAICE */
463
464 #if (defined (ALLOW_AUTODIFF_TAMC) && \
465 defined (ALLOW_AUTODIFF_MONITOR))
466 C Include call to a dummy routine. Its adjoint will be
467 C called at the proper place in the adjoint code.
468 C The adjoint routine will print out adjoint values
469 C if requested. The location of the call is important,
470 C it has to be after the adjoint of the exchanges
471 C (DO_GTERM_BLOCKING_EXCHANGES).
472 myiter = niter0 + iloop
473 mytime = starttime + float(iloop)*deltaTClock
474 CALL DUMMY_IN_STEPPING( myTime, myIter, myThid )
475 CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
476 myiter = nIter0 + (iloop-1)
477 mytime = startTime + float(iloop-1)*deltaTclock
478 #endif
479
480 cph(
481 #ifdef ALLOW_SNAPSHOTS
482 yprefix = 'kf_'
483 call ecco_check_exp( mythid, myiter, mytime, yprefix )
484 #endif
485 cph)
486
487 C-- Step forward fields and calculate time tendency terms.
488 #ifndef DISABLE_DEBUGMODE
489 IF ( debugLevel .GE. debLevB )
490 & CALL DEBUG_CALL('THERMODYNAMICS',myThid)
491 #endif
492 CALL TIMER_START('THERMODYNAMICS [FORWARD_STEP]',mythid)
493 CALL THERMODYNAMICS( myTime, myIter, myThid )
494 CALL TIMER_STOP ('THERMODYNAMICS [FORWARD_STEP]',mythid)
495
496 C-- do exchanges (needed for DYNAMICS) when using stagger time-step :
497 #ifndef DISABLE_DEBUGMODE
498 IF ( debugLevel .GE. debLevB )
499 & CALL DEBUG_CALL('DO_STAGGER_FIELDS_EXCH.',myThid)
500 #endif
501 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
502 CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
503 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
504
505 #ifdef ALLOW_SHAP_FILT
506 IF (useSHAP_FILT .AND.
507 & staggerTimeStep .AND. shap_filt_TrStagg ) THEN
508 #ifndef DISABLE_DEBUGMODE
509 IF ( debugLevel .GE. debLevB )
510 & CALL DEBUG_CALL('SHAP_FILT_APPLY_TS',myThid)
511 #endif
512 CALL TIMER_START('SHAP_FILT [FORWARD_STEP]',myThid)
513 CALL SHAP_FILT_APPLY_TS( gT, gS, myTime, myIter, myThid )
514 CALL TIMER_STOP ('SHAP_FILT [FORWARD_STEP]',myThid)
515 ENDIF
516 #endif
517
518 #ifdef ALLOW_ZONAL_FILT
519 IF (useZONAL_FILT .AND.
520 & staggerTimeStep .AND. zonal_filt_TrStagg ) THEN
521 #ifndef DISABLE_DEBUGMODE
522 IF ( debugLevel .GE. debLevB )
523 & CALL DEBUG_CALL('ZONAL_FILT_APPLY_TS',myThid)
524 #endif
525 CALL TIMER_START('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
526 CALL ZONAL_FILT_APPLY_TS( gT, gS, myThid )
527 CALL TIMER_STOP ('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
528 ENDIF
529 #endif
530
531 C-- Step forward fields and calculate time tendency terms.
532 #ifndef ALLOW_AUTODIFF_TAMC
533 IF ( momStepping ) THEN
534 #endif
535 #ifndef DISABLE_DEBUGMODE
536 IF ( debugLevel .GE. debLevB )
537 & CALL DEBUG_CALL('DYNAMICS',myThid)
538 #endif
539 CALL TIMER_START('DYNAMICS [FORWARD_STEP]',mythid)
540 CALL DYNAMICS( myTime, myIter, myThid )
541 CALL TIMER_STOP ('DYNAMICS [FORWARD_STEP]',mythid)
542 #ifndef ALLOW_AUTODIFF_TAMC
543 ENDIF
544 #endif
545
546 #ifdef ALLOW_NONHYDROSTATIC
547 C-- Step forward W field in N-H algorithm
548 IF ( momStepping .AND. nonHydrostatic ) THEN
549 #ifndef DISABLE_DEBUGMODE
550 IF ( debugLevel .GE. debLevB )
551 & CALL DEBUG_CALL('CALC_GW',myThid)
552 #endif
553 CALL TIMER_START('CALC_GW [FORWARD_STEP]',myThid)
554 CALL CALC_GW(myThid)
555 CALL TIMER_STOP ('CALC_GW [FORWARD_STEP]',myThid)
556 ENDIF
557 #endif
558
559 #ifdef NONLIN_FRSURF
560 C-- update hfacC,W,S and recip_hFac according to etaH(n+1) :
561 IF ( nonlinFreeSurf.GT.0) THEN
562 IF ( select_rStar.GT.0 ) THEN
563 CALL TIMER_START('UPDATE_R_STAR [FORWARD_STEP]',myThid)
564 CALL UPDATE_R_STAR( myTime, myIter, myThid )
565 CALL TIMER_STOP ('UPDATE_R_STAR [FORWARD_STEP]',myThid)
566 ELSE
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 C- update also CG2D matrix (and preconditioner)
573 IF ( momStepping .AND. nonlinFreeSurf.GT.2 ) THEN
574 CALL TIMER_START('UPDATE_CG2D [FORWARD_STEP]',myThid)
575 CALL UPDATE_CG2D( myTime, myIter, myThid )
576 CALL TIMER_STOP ('UPDATE_CG2D [FORWARD_STEP]',myThid)
577 ENDIF
578 #endif
579
580 C-- Apply Filters to u*,v* before SOLVE_FOR_PRESSURE
581 #ifdef ALLOW_SHAP_FILT
582 IF (useSHAP_FILT .AND. shap_filt_uvStar) THEN
583 CALL TIMER_START('SHAP_FILT [FORWARD_STEP]',myThid)
584 IF (implicDiv2Dflow.LT.1.) THEN
585 C-- Explicit+Implicit part of the Barotropic Flow Divergence
586 C => Filtering of uVel,vVel is necessary
587 CALL SHAP_FILT_APPLY_UV( uVel,vVel,
588 & myTime+deltaT, myIter+1, myThid )
589 ENDIF
590 CALL SHAP_FILT_APPLY_UV( gU,gV,myTime+deltaT,myIter+1,myThid)
591 CALL TIMER_STOP ('SHAP_FILT [FORWARD_STEP]',myThid)
592 ENDIF
593 #endif
594
595 #ifdef ALLOW_ZONAL_FILT
596 IF (useZONAL_FILT .AND. zonal_filt_uvStar) THEN
597 CALL TIMER_START('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
598 IF (implicDiv2Dflow.LT.1.) THEN
599 C-- Explicit+Implicit part of the Barotropic Flow Divergence
600 C => Filtering of uVel,vVel is necessary
601 CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
602 ENDIF
603 CALL ZONAL_FILT_APPLY_UV( gU, gV, myThid )
604 CALL TIMER_STOP ('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
605 ENDIF
606 #endif
607
608 C-- Solve elliptic equation(s).
609 C Two-dimensional only for conventional hydrostatic or
610 C three-dimensional for non-hydrostatic and/or IGW scheme.
611 IF ( momStepping ) THEN
612 CALL TIMER_START('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
613 CALL SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
614 CALL TIMER_STOP ('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
615 ENDIF
616
617 #ifdef ALLOW_AUTODIFF_TAMC
618 cph This is needed because convective_adjustment calls
619 cph find_rho which may use pressure()
620 CADJ STORE totphihyd = comlev1, key = ikey_dynamics
621 #endif
622 C-- Correct divergence in flow field and cycle time-stepping
623 C arrays (for all fields) ; update time-counter
624 myIter = nIter0 + iLoop
625 myTime = startTime + deltaTClock * float(iLoop)
626 CALL TIMER_START('THE_CORRECTION_STEP [FORWARD_STEP]',myThid)
627 CALL THE_CORRECTION_STEP(myTime, myIter, myThid)
628 CALL TIMER_STOP ('THE_CORRECTION_STEP [FORWARD_STEP]',myThid)
629
630 C-- Do "blocking" sends and receives for tendency "overlap" terms
631 c CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
632 c CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )
633 c CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
634
635 C-- Do "blocking" sends and receives for field "overlap" terms
636 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
637 CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
638 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
639
640 #ifdef ALLOW_FLT
641 C-- Calculate float trajectories
642 IF (useFLT) THEN
643 CALL TIMER_START('FLOATS [FORWARD_STEP]',myThid)
644 CALL FLT_MAIN(myIter,myTime, myThid)
645 CALL TIMER_STOP ('FLOATS [FORWARD_STEP]',myThid)
646 ENDIF
647 #endif
648
649 #ifdef ALLOW_MONITOR
650 C-- Check status of solution (statistics, cfl, etc...)
651 CALL TIMER_START('MONITOR [FORWARD_STEP]',myThid)
652 CALL MONITOR( myIter, myTime, myThid )
653 CALL TIMER_STOP ('MONITOR [FORWARD_STEP]',myThid)
654 #endif /* ALLOW_MONITOR */
655
656 C-- Do IO if needed.
657 CALL TIMER_START('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
658 CALL DO_THE_MODEL_IO( myTime, myIter, myThid )
659 CALL TIMER_STOP ('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
660
661 C-- Save state for restarts
662 C Note: (jmc: is it still the case after ckp35 ?)
663 C =====
664 C Because of the ordering of the timestepping code and
665 C tendency term code at end of loop model arrays hold
666 C U,V,T,S at "time-level" N but gu, gv, gs, gt, guNM1,...
667 C at "time-level" N+1/2 (guNM1 at "time-level" N+1/2 is
668 C gu at "time-level" N-1/2) and etaN at "time-level" N+1/2.
669 C where N = I+timeLevBase-1
670 C Thus a checkpoint contains U.0000000000, GU.0000000001 and
671 C etaN.0000000001 in the indexing scheme used for the model
672 C "state" files. This example is referred to as a checkpoint
673 C at time level 1
674 CALL TIMER_START('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
675 CALL WRITE_CHECKPOINT(
676 & .FALSE., myTime, myIter, myThid )
677 CALL TIMER_STOP ('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
678
679 ctest myiter = niter0 + iloop
680 ctest mytime = starttime + float(iloop)*deltaTClock
681 c Calculate current date.
682 ctest call cal_TimeStamp( myiter, mytime, mydate, mythid )
683
684 #ifdef ALLOW_AUTODIFF_TAMC
685 #ifdef ALLOW_TAMC_CHECKPOINTING
686 endif
687 enddo
688 endif
689 enddo
690 endif
691 enddo
692 #else
693 enddo
694 #endif
695
696 #else
697 enddo
698 #endif /* ALLOW_AUTODIFF_TAMC */
699
700 _BARRIER
701 call timer_stop ('ECCO MAIN LOOP', mythid)
702
703 call timer_start('ECCO SPIN-DOWN', mythid)
704
705 #ifdef ALLOW_COST
706
707 #ifdef ALLOW_DIVIDED_ADJOINT
708 CADJ STORE mytime = onetape
709 #endif
710 c-- Accumulate time averages of temperature, salinity, and SSH.
711 #ifndef DISABLE_DEBUGMODE
712 IF ( debugLevel .GE. debLevB )
713 & CALL DEBUG_CALL('cost_averagesfields',myThid)
714 #endif
715 call timer_start('cost_averagesfields [ECCO SPIN-DOWN]', mythid)
716 call cost_averagesfields( mytime, mythid )
717 call timer_stop ('cost_averagesfields [ECCO SPIN-DOWN]', mythid)
718 #ifdef ALLOW_DIVIDED_ADJOINT
719 c**************************************
720 #include "cost_averages_bar_directives.h"
721 c**************************************
722 #endif
723
724 #ifdef ALLOW_COST_ATLANTIC
725 c-- Compute meridional heat transport
726 #ifndef DISABLE_DEBUGMODE
727 IF ( debugLevel .GE. debLevB )
728 & CALL DEBUG_CALL('cost_atlantic',myThid)
729 #endif
730 call timer_start('cost_atlantic [ECCO SPIN-DOWN]', mythid)
731 call cost_atlantic( mytime, myiter,mythid )
732 call timer_stop ('cost_atlantic [ECCO SPIN-DOWN]', mythid)
733 #endif
734
735 c-- Compute cost function contribution of Temperature and Salinity.
736 #ifndef DISABLE_DEBUGMODE
737 IF ( debugLevel .GE. debLevB )
738 & CALL DEBUG_CALL('cost_hyd',myThid)
739 #endif
740 call timer_start('cost_hyd [ECCO SPIN-DOWN]', mythid)
741 call cost_hyd( myiter, mytime, mythid )
742 call timer_stop ('cost_hyd [ECCO SPIN-DOWN]', mythid)
743
744 #ifdef ALLOW_CURMTR_COST_CONTRIBUTION
745 #ifndef DISABLE_DEBUGMODE
746 IF ( debugLevel .GE. debLevB )
747 & CALL DEBUG_CALL('cost_curmtr',myThid)
748 #endif
749 call timer_start('cost_curmtr [ECCO SPIN-DOWN]', mythid)
750 call cost_curmtr( myiter, mytime, mythid )
751 call timer_stop ('cost_curmtr [ECCO SPIN-DOWN]', mythid)
752 #endif
753
754 c-- Compute cost function contribution of SSH.
755 #ifdef ALLOW_SSH_COST_CONTRIBUTION
756 #ifndef DISABLE_DEBUGMODE
757 IF ( debugLevel .GE. debLevB )
758 & CALL DEBUG_CALL('cost_ssh',myThid)
759 #endif
760 call timer_start('cost_ssh [ECCO SPIN-DOWN]', mythid)
761 call cost_ssh( myiter, mytime, mythid )
762 call timer_stop ('cost_ssh [ECCO SPIN-DOWN]', mythid)
763 #endif
764
765 c-- Compute cost function contribution of SSS.
766 #ifdef ALLOW_SST_COST_CONTRIBUTION
767 #ifndef DISABLE_DEBUGMODE
768 IF ( debugLevel .GE. debLevB )
769 & CALL DEBUG_CALL('cost_sst',myThid)
770 #endif
771 call timer_start('cost_sst [ECCO SPIN-DOWN]', mythid)
772 call cost_sst( myiter, mytime, mythid )
773 call timer_stop ('cost_sst [ECCO SPIN-DOWN]', mythid)
774 #endif
775
776 c-- Compute cost function contribution of SSS.
777 #ifdef ALLOW_SSS_COST_CONTRIBUTION
778 #ifndef DISABLE_DEBUGMODE
779 IF ( debugLevel .GE. debLevB )
780 & CALL DEBUG_CALL('cost_sss',myThid)
781 #endif
782 call timer_start('cost_sss [ECCO SPIN-DOWN]', mythid)
783 call cost_sss( myiter, mytime, mythid )
784 call timer_stop ('cost_sss [ECCO SPIN-DOWN]', mythid)
785 #endif
786
787 c-- Compute cost function contribution of drifter's velocities.
788 #ifdef ALLOW_DRIFTER_COST_CONTRIBUTION
789 #ifndef DISABLE_DEBUGMODE
790 IF ( debugLevel .GE. debLevB )
791 & CALL DEBUG_CALL('cost_drifter',myThid)
792 #endif
793 call timer_start('cost_drifter [ECCO SPIN-DOWN]', mythid)
794 call cost_drifter( myiter, mytime, mythid )
795 call timer_stop ('cost_drifter [ECCO SPIN-DOWN]', mythid)
796 #endif
797
798 c-- Compute cost function contribution of wind stress observations.
799 #ifdef ALLOW_SCAT_COST_CONTRIBUTION
800 #ifndef DISABLE_DEBUGMODE
801 IF ( debugLevel .GE. debLevB )
802 & CALL DEBUG_CALL('cost_scat',myThid)
803 #endif
804 call timer_start('cost_scat [ECCO SPIN-DOWN]', mythid)
805 call cost_scat( myiter, mytime, mythid )
806 call timer_stop ('cost_scat [ECCO SPIN-DOWN]', mythid)
807 #endif
808
809 c-- Compute cost function contribution of wind stress observations.
810 #ifdef ALLOW_MEAN_HFLUX_COST_CONTRIBUTION
811 call timer_start('cost_mean_heatflux [ECCO SPIN-DOWN]', mythid)
812 call cost_mean_heatflux( myiter, mytime, mythid )
813 call timer_stop ('cost_mean_heatflux [ECCO SPIN-DOWN]', mythid)
814 #endif
815
816 c-- Compute cost function contribution of wind stress observations.
817 #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
818 call timer_start('cost_mean_saltflux [ECCO SPIN-DOWN]', mythid)
819 call cost_mean_saltflux( myiter, mytime, mythid )
820 call timer_stop ('cost_mean_saltflux [ECCO SPIN-DOWN]', mythid)
821 #endif
822
823 c-- Compute cost function contribution of drift between the first
824 c and the last year.
825 #ifdef ALLOW_DRIFT_COST_CONTRIBUTION
826 #ifndef DISABLE_DEBUGMODE
827 IF ( debugLevel .GE. debLevB )
828 & CALL DEBUG_CALL('cost_drift',myThid)
829 #endif
830 call timer_start('cost_drift [ECCO SPIN-DOWN]', mythid)
831 call cost_drift( myiter, mytime, mythid )
832 call timer_stop ('cost_drift [ECCO SPIN-DOWN]', mythid)
833 #endif
834 #ifdef ALLOW_DRIFTW_COST_CONTRIBUTION
835 #ifndef DISABLE_DEBUGMODE
836 IF ( debugLevel .GE. debLevB )
837 & CALL DEBUG_CALL('cost_driftw',myThid)
838 #endif
839 call timer_start('cost_driftw [ECCO SPIN-DOWN]', mythid)
840 call cost_driftw( myiter, mytime, mythid )
841 call timer_stop ('cost_driftw [ECCO SPIN-DOWN]', mythid)
842 #endif
843 _BARRIER
844
845 c-- Compute initial vs. final T/S deviation
846 #ifdef ALLOW_COST_INI_FIN
847 call timer_start('cost_ini_fin [ECCO SPIN-DOWN]', mythid)
848 call cost_theta_ini_fin( myiter, mytime, mythid )
849 call cost_salt_ini_fin( myiter, mytime, mythid )
850 call timer_stop ('cost_ini_fin [ECCO SPIN-DOWN]', mythid)
851 #endif
852 _BARRIER
853
854 c-- Sum all cost function contributions.
855 #ifndef DISABLE_DEBUGMODE
856 IF ( debugLevel .GE. debLevB )
857 & CALL DEBUG_CALL('cost_final',myThid)
858 #endif
859 call timer_start('COST_FINAL [ECCO SPIN-DOWN]', mythid)
860 call cost_final( mythid )
861 call timer_stop ('COST_FINAL [ECCO SPIN-DOWN]', mythid)
862
863 #endif /* ALLOW_COST */
864
865 call timer_stop ('ECCO SPIN-DOWN', mythid)
866
867 #ifndef DISABLE_DEBUGMODE
868 IF ( debugLevel .GE. debLevB )
869 & CALL DEBUG_LEAVE('THE_MAIN_LOOP',myThid)
870 #endif
871
872 return
873 end
874

  ViewVC Help
Powered by ViewVC 1.1.22