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

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

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


Revision 1.5 - (hide annotations) (download)
Wed Mar 24 21:49:53 2004 UTC (20 years, 3 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint52n_post, checkpoint52m_post, checkpoint53
Changes since 1.4: +5 -3 lines
packages_write_pickup added to timestepping loop.

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

  ViewVC Help
Powered by ViewVC 1.1.22