/[MITgcm]/MITgcm/model/src/forward_step.F
ViewVC logotype

Diff of /MITgcm/model/src/forward_step.F

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

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