/[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.2 - (hide 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 heimbach 1.2 C $Header: /u/gcmpack/MITgcm/pkg/ecco/the_main_loop.F,v 1.1 2003/11/06 22:10:08 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    
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 heimbach 1.2 #include "cost_averages_bar_directives.h"
721 heimbach 1.1 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