/[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.6 by cnh, Sun Feb 4 14:38:47 2001 UTC revision 1.132 by heimbach, Thu Mar 2 23:41:10 2006 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                    doHalfStep, 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  |    |  |--DO_THE_MODEL_IO  
 C /|\   |  |   o Write model state  
 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  |    |  |--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  #include "CG2D.h"  
43  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_MNC
44  #include "CG3D.h"  #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 NONLIN_FRSURF
72    #  include "GRID.h"
73    # endif
74    # ifdef ALLOW_EXF
75    #  include "exf_fields.h"
76    #  include "exf_clim_fields.h"
77    #  ifdef ALLOW_BULKFORMULAE
78    #   include "exf_constants.h"
79    #  endif
80    # endif
81    # ifdef ALLOW_OBCS
82    #  include "OBCS.h"
83    # endif
84    # ifdef ALLOW_PTRACERS
85    #  include "PTRACERS_SIZE.h"
86    #  include "PTRACERS.h"
87    # endif
88    # ifdef ALLOW_CD_CODE
89    #  include "CD_CODE_VARS.h"
90    # endif
91    # ifdef ALLOW_EBM
92    #  include "EBM.h"
93    # endif
94    # ifdef EXACT_CONSERV
95    #  include "SURFACE.h"
96    # endif
97    # ifdef ALLOW_KPP
98    #  include "KPP.h"
99    # endif
100    # ifdef ALLOW_GMREDI
101    #  include "GMREDI.h"
102    # endif
103    #endif /* ALLOW_AUTODIFF_TAMC */
104    
105    C     !LOCAL VARIABLES:
106  C     == Routine arguments ==  C     == Routine arguments ==
107  C     doHalfStep - If .TRUE. then this is the last half step  C     note: under the multi-threaded model myiter and
108  C     iLoop         - Invocation count (counter in THE_MAIN_LOOP)  C           mytime are local variables passed around as routine
109  C     myCurrentIter - Iteration counter for this thread  C           arguments. Although this is fiddly it saves the need to
110  C     myCurrentTime - Time counter for this thread  C           impose additional synchronisation points when they are
111  C     myThid - Thread number for this instance of the routine.  C           updated.
112        LOGICAL doHalfStep  C     myIter - iteration counter for this thread
113        INTEGER iLoop  C     myTime - time counter for this thread
114        INTEGER myCurrentIter  C     myThid - thread number for this instance of the routine.
115        _RL     myCurrentTime        INTEGER iloop
116        INTEGER myThid          INTEGER myThid
117          INTEGER myIter
118          _RL     myTime
119    
120  C     == Local variables ==  C     == Local variables ==
121    #ifdef COMPONENT_MODULE
122          INTEGER myItP1
123    #endif
124    CEOP
125    
126    #ifdef ALLOW_DEBUG
127          IF ( debugLevel .GE. debLevB )
128         &    CALL DEBUG_ENTER('FORWARD_STEP',myThid)
129    #endif
130    
131    #ifdef ALLOW_AUTODIFF_TAMC
132    C--   Reset the model iteration counter and the model time.
133          myIter = nIter0 + (iloop-1)
134          myTime = startTime + float(iloop-1)*deltaTclock
135    #endif
136    
137  C--   Load forcing/external data fields  #ifdef ALLOW_AUTODIFF_TAMC
138        CALL TIMER_START('I/O (READ)         [FORWARD_STEP]',myThid)  c**************************************
139        CALL EXTERNAL_FIELDS_LOAD( myCurrentTime, myCurrentIter, myThid )  #include "checkpoint_lev1_directives.h"
140        CALL TIMER_STOP ('I/O (READ)         [FORWARD_STEP]',myThid)  c**************************************
141    #endif
142  C--   Step forward fields and calculate time tendency terms  
143        CALL TIMER_START('DYNAMICS           [FORWARD_STEP]',myThid)  C--   Switch on/off diagnostics for snap-shot output:
144        CALL DYNAMICS( myCurrentTime, myCurrentIter, myThid )  #ifdef ALLOW_DIAGNOSTICS
145        CALL TIMER_STOP ('DYNAMICS           [FORWARD_STEP]',myThid)        IF ( useDiagnostics ) THEN
146            CALL DIAGNOSTICS_SWITCH_ONOFF( myTime, myIter, myThid )
 #ifdef ALLOW_NONHYDROSTATIC  
 C--   Step forward W field in N-H algorithm  
       IF ( nonHydrostatic ) THEN  
        CALL TIMER_START('CALC_GW            [FORWARD_STEP]',myThid)  
        CALL CALC_GW( myThid)  
        CALL TIMER_STOP ('CALC_GW            [FORWARD_STEP]',myThid)  
147        ENDIF        ENDIF
148  #endif  #endif
149    
150  #ifdef ALLOW_SHAP_FILT  C--   State-variables diagnostics
151  C--   Step forward all tiles, filter and exchange.        IF ( usediagnostics ) THEN
152        CALL TIMER_START('SHAP_FILT          [FORWARD_STEP]',myThid)          CALL TIMER_START('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
153        CALL SHAP_FILT_APPLY(          CALL DO_STATEVARS_DIAGS( myTime, 0, myIter, myThid )
154       I                     gUnm1, gVnm1, gTnm1, gSnm1,          CALL TIMER_STOP ('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
155       I                     myCurrentTime, myCurrentIter, myThid )        ENDIF
156        CALL TIMER_STOP ('SHAP_FILT          [FORWARD_STEP]',myThid)  
157    C--   Call Bulk-Formulae forcing package
158    #ifdef ALLOW_BULK_FORCE
159          IF ( useBulkForce ) THEN
160    #ifdef ALLOW_DEBUG
161           IF ( debugLevel .GE. debLevB )
162         &    CALL DEBUG_CALL('BULKF_FIELDS_LOAD',myThid)
163  #endif  #endif
164           CALL TIMER_START('BULKF_FORCING      [FORWARD_STEP]',mythid)
165    C-    load all forcing fields at current time
166           CALL BULKF_FIELDS_LOAD( myTime, myIter, myThid )
167    C-    calculate qnet and empmr (and wind stress)
168           CALL BULKF_FORCING( myTime, myIter, myThid )
169           CALL TIMER_STOP ('BULKF_FORCING      [FORWARD_STEP]',mythid)
170          ENDIF
171    #endif /* ALLOW_BULK_FORCE */
172    
173  #ifdef ALLOW_ZONAL_FILT  C--   Call external forcing package
174        CALL ZONAL_FILT_APPLY(  # ifdef ALLOW_EXF
175       U           gUnm1, gVnm1, gTnm1, gSnm1,  #  ifdef ALLOW_DEBUG
176       I           myThid )        IF ( debugLevel .GE. debLevB )
177         &    CALL DEBUG_CALL('EXF_GETFORCING',myThid)
178    #  endif
179          CALL TIMER_START('EXF_GETFORCING     [FORWARD_STEP]',mythid)
180          CALL EXF_GETFORCING( mytime, myiter, mythid )
181          CALL TIMER_STOP ('EXF_GETFORCING     [FORWARD_STEP]',mythid)
182    # else /* ALLOW_EXF undef */
183    cph The following IF-statement creates an additional dependency
184    cph for the forcing fields requiring additional storing.
185    cph Therefore, the IF-statement will be put between CPP-OPTIONS,
186    cph assuming that ALLOW_SEAICE has not yet been differentiated.
187    #  if (defined (ALLOW_SEAICE) || defined (ALLOW_EBM))
188          IF ( .NOT. useSEAICE .AND. .NOT. useEBM ) THEN
189    #  endif
190    #ifdef ALLOW_DEBUG
191           IF ( debugLevel .GE. debLevB )
192         &    CALL DEBUG_CALL('EXTERNAL_FIELDS_LOAD',myThid)
193  #endif  #endif
194           CALL TIMER_START('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
195           CALL EXTERNAL_FIELDS_LOAD( mytime, myiter, mythid )
196           CALL TIMER_STOP ('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
197    #  if (defined (ALLOW_SEAICE) || defined (ALLOW_EBM))
198          ENDIF
199    #  endif
200    # endif /* ALLOW_EXF */
201    
202  C--   Do IO if needed.  #ifdef ALLOW_AUTODIFF
203  C     Note:  c--   Add control vector for forcing and parameter fields
204  C     =====        if ( myiter .EQ. nIter0 )
205  C     At this point model arrays hold U,V,T,S  at "time-level" N       &     CALL CTRL_MAP_FORCING (mythid)
206  C     and cg2d_x at "time-level" N-1/2 where N = I+timeLevBase-1.  #endif
207  C     By convention this is taken to be the model "state".  
208        CALL TIMER_START('I/O (WRITE)        [FORWARD_STEP]',myThid)  #if (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_AUTODIFF_MONITOR))
209        CALL DO_THE_MODEL_IO(  C     Include call to a dummy routine. Its adjoint will be
210       &        doHalfStep, myCurrentTime, myCurrentIter, myThid )  C     called at the proper place in the adjoint code.
211        CALL TIMER_STOP ('I/O (WRITE)        [FORWARD_STEP]',myThid)  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    # ifdef NONLIN_FRSURF
320    cph-test
321    CADJ STORE hFac_surfC         = comlev1, key = ikey_dynamics
322    CADJ STORE hfac_surfs         = comlev1, key = ikey_dynamics
323    CADJ STORE hfac_surfw         = comlev1, key = ikey_dynamics
324    CADJ STORE hFacC, hFacS, hFacW
325    CADJ &     = comlev1, key = ikey_dynamics
326    CADJ STORE recip_hFacC, recip_hFacS, recip_hFacW
327    CADJ &     = comlev1, key = ikey_dynamics
328    # endif
329    #endif /* ALLOW_AUTODIFF_TAMC */
330    
331    #ifndef ALLOW_OFFLINE
332    #ifdef ALLOW_DEBUG
333           IF ( debugLevel .GE. debLevB )
334         &    CALL DEBUG_CALL('DO_OCEANIC_PHYS',myThid)
335    #endif
336           CALL TIMER_START('DO_OCEANIC_PHYS     [FORWARD_STEP]',mythid)
337           CALL DO_OCEANIC_PHYS( myTime, myIter, myThid )
338           CALL TIMER_STOP ('DO_OCEANIC_PHYS     [FORWARD_STEP]',mythid)
339    #ifdef ALLOW_AUTODIFF_TAMC
340    CADJ STORE EmPmR    = comlev1, key = ikey_dynamics
341    #endif
342    #endif
343    
344        IF (.NOT. doHalfStep) THEN  #ifdef ALLOW_GCHEM
345    C     GCHEM package is an interface for any bio-geochemical or
346    C     ecosystem model you would like to include.
347    C     If GCHEM_SEPARATE_FORCING is not defined, you are
348    C     responsible for computing tendency terms for passive
349    C     tracers and storing them on a 3DxNumPtracers-array called
350    C     gchemTendency in GCHEM_CALC_TENDENCY. This tendency is then added
351    C     to gPtr in ptracers_forcing later-on.
352    C     If GCHEM_SEPARATE_FORCING is defined, you are reponsible for
353    C     UPDATING ptracers directly in GCHEM_FORCING_SEP. This amounts
354    C     to a completely separate time step that you have to implement
355    C     yourself (Eulerian seems to be fine in most cases).
356    CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
357    C     CAVEAT: Up to now, when GCHEM is turned on the field ptracerForcingSurf,
358    C     which is needed for KPP is not set properly. ptracerForcingSurf must
359    C     be treated differently depending on whether GCHEM_SEPARATE_FORCING
360    C     is define or not. TBD.
361    CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
362           IF ( useGCHEM ) THEN
363    #ifdef ALLOW_DEBUG
364            IF ( debugLevel .GE. debLevB )
365         &       CALL DEBUG_CALL('GCHEM_CALC_TENDENCY',myThid)
366    #endif
367            CALL TIMER_START('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid)
368            CALL GCHEM_CALC_TENDENCY( myTime, myIter, myThid )
369            CALL TIMER_STOP ('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid)
370           ENDIF
371    #endif /* ALLOW_GCHEM */
372    
373    #ifdef ALLOW_AUTODIFF_TAMC
374    cph needed to be moved here from do_oceanic_physics
375    cph to be visible down the road
376    c
377    CADJ STORE surfaceForcingS    = comlev1, key = ikey_dynamics
378    CADJ STORE surfaceForcingT    = comlev1, key = ikey_dynamics
379    CADJ STORE surfaceForcingTice = comlev1, key = ikey_dynamics
380    ctest(
381    CADJ STORE IVDConvCount       = comlev1, key = ikey_dynamics
382    ctest)
383    # ifdef ALLOW_PTRACERS
384    CADJ STORE surfaceForcingPtr  = comlev1, key = ikey_dynamics
385    # endif
386    c
387    # ifdef ALLOW_GMREDI
388    CADJ STORE Kwx                = comlev1, key = ikey_dynamics
389    CADJ STORE Kwy                = comlev1, key = ikey_dynamics
390    CADJ STORE Kwz                = comlev1, key = ikey_dynamics
391    #  ifdef GM_BOLUS_ADVEC
392    CADJ STORE GM_PsiX            = comlev1, key = ikey_dynamics
393    CADJ STORE GM_PsiY            = comlev1, key = ikey_dynamics
394    #  endif
395    # endif
396    c
397    # ifdef ALLOW_KPP
398    CADJ STORE KPPghat            = comlev1, key = ikey_dynamics
399    CADJ STORE KPPfrac            = comlev1, key = ikey_dynamics
400    CADJ STORE KPPdiffKzS         = comlev1, key = ikey_dynamics
401    CADJ STORE KPPdiffKzT         = comlev1, key = ikey_dynamics
402    # endif
403    c
404    # ifdef NONLIN_FRSURF
405    cph-test
406    CADJ STORE etaN               = comlev1, key = ikey_dynamics
407    CADJ STORE etaH               = comlev1, key = ikey_dynamics
408    CADJ STORE uVel               = comlev1, key = ikey_dynamics
409    CADJ STORE vVel               = comlev1, key = ikey_dynamics
410    #  ifdef ALLOW_CD_CODE
411    CADJ STORE etanm1             = comlev1, key = ikey_dynamics
412    #  endif
413    # endif
414    #endif /* ALLOW_AUTODIFF_TAMC */
415    
416          IF ( .NOT.staggerTimeStep ) THEN
417    #ifdef ALLOW_DEBUG
418            IF ( debugLevel .GE. debLevB )
419         &    CALL DEBUG_CALL('THERMODYNAMICS',myThid)
420    #endif
421            CALL TIMER_START('THERMODYNAMICS      [FORWARD_STEP]',mythid)
422            CALL THERMODYNAMICS( myTime, myIter, myThid )
423            CALL TIMER_STOP ('THERMODYNAMICS      [FORWARD_STEP]',mythid)
424    C--    if not staggerTimeStep: end
425          ENDIF
426    c #ifdef ALLOW_NONHYDROSTATIC
427          IF ( implicitIntGravWave ) THEN
428            CALL TIMER_START('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
429            CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
430            CALL TIMER_STOP ('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
431          ENDIF
432    c #endif
433    
434    #ifdef COMPONENT_MODULE
435           IF ( useCoupler .AND. .NOT.cpl_earlyExpImpCall ) THEN
436    C      Post coupling data that I export.
437    C      Read in coupling data that I import.
438             myItP1 = myIter + 1
439             CALL TIMER_START('CPL_EXPORT-IMPORT  [FORWARD_STEP]',myThid)
440             CALL CPL_EXPORT_MY_DATA(       myItP1, myTime, myThid )
441             CALL CPL_IMPORT_EXTERNAL_DATA( myItP1, myTime, myThid )
442             CALL TIMER_STOP ('CPL_EXPORT-IMPORT  [FORWARD_STEP]',myThid)
443    # ifndef ALLOW_AIM
444            IF ( useRealFreshWaterFlux ) THEN
445             CALL OCN_APPLY_IMPORT( .FALSE., myTime, myIter, myThid )
446            ENDIF
447    # endif
448           ENDIF
449    #endif /* COMPONENT_MODULE */
450    
451    C--   Step forward fields and calculate time tendency terms.
452    #ifndef ALLOW_OFFLINE
453    #ifndef ALLOW_AUTODIFF_TAMC
454          IF ( momStepping ) THEN
455    #endif
456    #ifdef ALLOW_DEBUG
457            IF ( debugLevel .GE. debLevB )
458         &    CALL DEBUG_CALL('DYNAMICS',myThid)
459    #endif
460            CALL TIMER_START('DYNAMICS            [FORWARD_STEP]',mythid)
461            CALL DYNAMICS( myTime, myIter, myThid )
462            CALL TIMER_STOP ('DYNAMICS            [FORWARD_STEP]',mythid)
463    #ifndef ALLOW_AUTODIFF_TAMC
464          ENDIF
465    #endif
466    #endif
467    
468    C--   Update time-counter
469          myIter = nIter0 + iLoop
470          myTime = startTime + deltaTClock * float(iLoop)
471    
472    #ifdef ALLOW_MNC
473    C     Update the default next iter for MNC
474          IF ( useMNC ) THEN
475             CALL MNC_CW_CITER_SETG( 1, 1, -1, myIter , myThid )
476    
477    C        TODO: Logic should be added here so that users can specify, on
478    C        a per-citer-group basis, when it is time to update the
479    C        "current" (and not just the "next") iteration
480    
481    C        TODO: the following is just a temporary band-aid (mostly, for
482    C        Baylor) until someone writes a routine that better handles time
483    C        boundaries such as weeks, months, years, etc.
484             IF ( mnc_filefreq .GT. 0 ) THEN
485               IF (DIFFERENT_MULTIPLE(mnc_filefreq,myTime,deltaTClock))
486         &          THEN
487                 CALL MNC_CW_CITER_SETG( 1, 1, myIter, -1 , myThid )
488               ENDIF
489             ENDIF
490           ENDIF
491    #endif
492    
493    C--   Update geometric factors:
494    #ifdef NONLIN_FRSURF
495    C-    update hfacC,W,S and recip_hFac according to etaH(n+1) :
496          IF ( nonlinFreeSurf.GT.0) THEN
497           IF ( select_rStar.GT.0 ) THEN
498    # ifndef DISABLE_RSTAR_CODE
499            CALL TIMER_START('UPDATE_R_STAR      [FORWARD_STEP]',myThid)
500            CALL UPDATE_R_STAR( myTime, myIter, myThid )
501            CALL TIMER_STOP ('UPDATE_R_STAR      [FORWARD_STEP]',myThid)
502    # endif /* DISABLE_RSTAR_CODE */
503           ELSE
504    #ifdef ALLOW_AUTODIFF_TAMC
505    cph-test
506    CADJ STORE hFac_surfC, hFac_surfS, hFac_surfW
507    CADJ &     = comlev1, key = ikey_dynamics
508    #endif
509            CALL TIMER_START('UPDATE_SURF_DR     [FORWARD_STEP]',myThid)
510            CALL UPDATE_SURF_DR( myTime, myIter, myThid )
511            CALL TIMER_STOP ('UPDATE_SURF_DR     [FORWARD_STEP]',myThid)
512           ENDIF
513          ENDIF
514    C-    update also CG2D matrix (and preconditioner)
515          IF ( momStepping .AND. nonlinFreeSurf.GT.2 ) THEN
516            CALL TIMER_START('UPDATE_CG2D        [FORWARD_STEP]',myThid)
517            CALL UPDATE_CG2D( myTime, myIter, myThid )
518            CALL TIMER_STOP ('UPDATE_CG2D        [FORWARD_STEP]',myThid)
519          ENDIF
520    #endif /* NONLIN_FRSURF */
521    
522    C--   Apply Filters to u*,v* before SOLVE_FOR_PRESSURE
523    #ifdef ALLOW_SHAP_FILT
524          IF (useSHAP_FILT .AND. shap_filt_uvStar) THEN
525            CALL TIMER_START('SHAP_FILT_UV        [FORWARD_STEP]',myThid)
526            IF (implicDiv2Dflow.LT.1.) THEN
527    C--   Explicit+Implicit part of the Barotropic Flow Divergence
528    C      => Filtering of uVel,vVel is necessary
529              CALL SHAP_FILT_APPLY_UV( uVel,vVel,
530         &                             myTime, myIter, myThid )
531            ENDIF
532            CALL SHAP_FILT_APPLY_UV( gU,gV,myTime,myIter,myThid)
533            CALL TIMER_STOP ('SHAP_FILT_UV        [FORWARD_STEP]',myThid)
534          ENDIF
535    #endif
536    #ifdef ALLOW_ZONAL_FILT
537          IF (useZONAL_FILT .AND. zonal_filt_uvStar) THEN
538            CALL TIMER_START('ZONAL_FILT_UV       [FORWARD_STEP]',myThid)
539            IF (implicDiv2Dflow.LT.1.) THEN
540    C--   Explicit+Implicit part of the Barotropic Flow Divergence
541    C      => Filtering of uVel,vVel is necessary
542              CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
543            ENDIF
544            CALL ZONAL_FILT_APPLY_UV( gU, gV, myThid )
545            CALL TIMER_STOP ('ZONAL_FILT_UV       [FORWARD_STEP]',myThid)
546          ENDIF
547    #endif  
548    
549  C--   Solve elliptic equation(s).  C--   Solve elliptic equation(s).
550  C     Two-dimensional only for conventional hydrostatic or  C     Two-dimensional only for conventional hydrostatic or
551  C     three-dimensional for non-hydrostatic and/or IGW scheme.  C     three-dimensional for non-hydrostatic and/or IGW scheme.
552        CALL TIMER_START('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)  #ifndef ALLOW_OFFLINE
553        CALL SOLVE_FOR_PRESSURE( myThid )        IF ( momStepping ) THEN
554        CALL TIMER_STOP ('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)          CALL TIMER_START('SOLVE_FOR_PRESSURE  [FORWARD_STEP]',myThid)
555            CALL SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
556  C--   Correct divergence in flow field and cycle time-stepping          CALL TIMER_STOP ('SOLVE_FOR_PRESSURE  [FORWARD_STEP]',myThid)
557  C     arrays (for all fields)        ENDIF
558        CALL THE_CORRECTION_STEP(myCurrentTime, myCurrentIter, myThid)  #endif
559    
560    C--   Correct divergence in flow field and cycle time-stepping momentum
561    c     IF ( momStepping ) THEN
562    #ifndef ALLOW_OFFLINE
563            CALL TIMER_START('MOM_CORRECTION_STEP [FORWARD_STEP]',myThid)
564            CALL MOMENTUM_CORRECTION_STEP(myTime, myIter, myThid)
565            CALL TIMER_STOP ('MOM_CORRECTION_STEP [FORWARD_STEP]',myThid)
566    #endif
567    c     ENDIF
568    
569    #ifdef EXACT_CONSERV
570          IF (exactConserv) THEN
571    C--   Update etaH(n+1) :
572            CALL TIMER_START('UPDATE_ETAH         [FORWARD_STEP]',mythid)
573            CALL UPDATE_ETAH( myTime, myIter, myThid )
574            CALL TIMER_STOP ('UPDATE_ETAH         [FORWARD_STEP]',mythid)
575          ENDIF
576    #endif /* EXACT_CONSERV */
577    
578    #ifdef NONLIN_FRSURF
579          IF ( select_rStar.NE.0 ) THEN
580    # ifndef DISABLE_RSTAR_CODE
581    C--   r* : compute the future level thickness according to etaH(n+1)
582            CALL TIMER_START('CALC_R_STAR       [FORWARD_STEP]',mythid)
583            CALL CALC_R_STAR(etaH, myTime, myIter, myThid )
584            CALL TIMER_STOP ('CALC_R_STAR       [FORWARD_STEP]',mythid)
585    # endif /* DISABLE_RSTAR_CODE */
586          ELSEIF ( nonlinFreeSurf.GT.0) THEN
587    C--   compute the future surface level thickness according to etaH(n+1)
588    #ifdef ALLOW_AUTODIFF_TAMC
589    CADJ STORE etaH          = comlev1, key = ikey_dynamics
590    #endif
591            CALL TIMER_START('CALC_SURF_DR      [FORWARD_STEP]',mythid)
592            CALL CALC_SURF_DR(etaH, myTime, myIter, myThid )
593            CALL TIMER_STOP ('CALC_SURF_DR      [FORWARD_STEP]',mythid)
594          ENDIF
595    #endif /* NONLIN_FRSURF */
596    
597    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
598          IF ( staggerTimeStep ) THEN
599    C--   do exchanges of U,V (needed for multiDim) when using stagger time-step :
600    #ifdef ALLOW_DEBUG
601            IF ( debugLevel .GE. debLevB )
602         &    CALL DEBUG_CALL('DO_STAGGER_FIELDS_EXCH.',myThid)
603    #endif
604            CALL TIMER_START('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
605            CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
606            CALL TIMER_STOP ('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
607    
608    C--   State-variables diagnostics
609           IF ( usediagnostics ) THEN
610            CALL TIMER_START('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
611            CALL DO_STATEVARS_DIAGS( myTime, 1, myIter, myThid )
612            CALL TIMER_STOP ('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
613           ENDIF
614    
615    #ifdef ALLOW_DEBUG
616            IF ( debugLevel .GE. debLevB )
617         &    CALL DEBUG_CALL('THERMODYNAMICS',myThid)
618    #endif
619            CALL TIMER_START('THERMODYNAMICS      [FORWARD_STEP]',mythid)
620            CALL THERMODYNAMICS( myTime, myIter, myThid )
621            CALL TIMER_STOP ('THERMODYNAMICS      [FORWARD_STEP]',mythid)
622    
623    C--    if staggerTimeStep: end
624          ENDIF
625    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
626    
627    #ifdef ALLOW_AUTODIFF_TAMC
628    cph This is needed because convective_adjustment calls
629    cph find_rho which may use pressure()
630    CADJ STORE totphihyd  = comlev1, key = ikey_dynamics
631    #endif
632    C--   Cycle time-stepping Tracers arrays (T,S,+pTracers)
633            CALL TIMER_START('TRC_CORRECTION_STEP [FORWARD_STEP]',myThid)
634            CALL TRACERS_CORRECTION_STEP(myTime, myIter, myThid)
635            CALL TIMER_STOP ('TRC_CORRECTION_STEP [FORWARD_STEP]',myThid)
636    
637    #ifdef ALLOW_GCHEM
638    C     Add separate timestepping of chemical/biological/forcing
639    C     of ptracers here in GCHEM_FORCING_SEP
640            IF ( useGCHEM ) THEN
641    #ifdef ALLOW_DEBUG
642             IF ( debugLevel .GE. debLevB )
643         &        CALL DEBUG_CALL('GCHEM_FORCING_SEP',myThid)
644    #endif /* ALLOW_DEBUG */
645             CALL TIMER_START('GCHEM_FORCING_SEP  [FORWARD_STEP]',myThid)
646             CALL GCHEM_FORCING_SEP( myTime,myIter,myThid )
647             CALL TIMER_STOP ('GCHEM_FORCING_SEP  [FORWARD_STEP]',myThid)
648            ENDIF  
649    #endif /* ALLOW_GCHEM */
650    
651  C--   Do "blocking" sends and receives for tendency "overlap" terms  C--   Do "blocking" sends and receives for tendency "overlap" terms
652  C     CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)  c     CALL TIMER_START('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
653  C     CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )  c     CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )
654  C     CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)  c     CALL TIMER_STOP ('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
655    
656  C--   Do "blocking" sends and receives for field "overlap" terms  C--   Do "blocking" sends and receives for field "overlap" terms
657        CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)        CALL TIMER_START('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
658        CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )        CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
659        CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)        CALL TIMER_STOP ('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
660    
661          IF ( usediagnostics ) THEN
662           CALL TIMER_START('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
663           CALL DO_STATEVARS_DIAGS( myTime, 2, myIter, myThid )
664           CALL TIMER_STOP ('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
665          ENDIF
666    
667  c     myCurrentIter = myCurrentIter + 1  C AMM
668  c     myCurrentTime = myCurrentTime + deltaTClock  #ifdef ALLOW_GRIDALT
669        myCurrentIter = nIter0 + iLoop          if (useGRIDALT) then
670        myCurrentTime = startTime + deltaTClock * float(iLoop)           CALL GRIDALT_UPDATE(myThid)
671            endif
672    #endif
673    C AMM
674    
675    C AMM
676    #ifdef ALLOW_FIZHI
677            if( useFIZHI) then
678             CALL TIMER_START('FIZHI               [FORWARD_STEP]',mythid)
679             CALL STEP_FIZHI_CORR ( myTime, myIter, myThid, dTtracerLev(1) )
680             CALL TIMER_STOP('FIZHI               [FORWARD_STEP]',mythid)
681            endif
682    #endif
683    C AMM
684    
685    #ifdef ALLOW_FLT
686    C--   Calculate float trajectories
687          IF (useFLT) THEN
688            CALL TIMER_START('FLOATS            [FORWARD_STEP]',myThid)
689            CALL FLT_MAIN(myIter,myTime, myThid)
690            CALL TIMER_STOP ('FLOATS            [FORWARD_STEP]',myThid)
691          ENDIF
692    #endif
693    
694    C--   State-variables time-averaging
695          CALL TIMER_START('DO_STATEVARS_TAVE   [FORWARD_STEP]',myThid)
696          CALL DO_STATEVARS_TAVE( myTime, myIter, myThid )
697          CALL TIMER_STOP ('DO_STATEVARS_TAVE   [FORWARD_STEP]',myThid)
698    
699    #ifndef ALLOW_OFFLINE
700    #ifdef ALLOW_MONITOR
701    C--   Check status of solution (statistics, cfl, etc...)
702          CALL TIMER_START('MONITOR             [FORWARD_STEP]',myThid)
703          CALL MONITOR( myIter, myTime, myThid )
704          CALL TIMER_STOP ('MONITOR             [FORWARD_STEP]',myThid)
705    #endif /* ALLOW_MONITOR */
706    #endif
707    
708    #ifdef ALLOW_COST
709    C--     compare model with data and compute cost function
710    C--     this is done after exchanges to allow interpolation
711          CALL TIMER_START('COST_TILE           [FORWARD_STEP]',myThid)
712          CALL COST_TILE  ( mytime, myiter, myThid )
713          CALL TIMER_STOP ('COST_TILE           [FORWARD_STEP]',myThid)
714    #endif
715    
716    C--   Do IO if needed.
717    #ifdef ALLOW_OFFLINE
718          CALL TIMER_START('OFFLINE_MODEL_IO    [FORWARD_STEP]',myThid)
719          CALL OFFLINE_MODEL_IO( myTime, myIter, myThid )
720          CALL TIMER_STOP ('OFFLINE_MODEL_IO    [FORWARD_STEP]',myThid)
721    #else
722          CALL TIMER_START('DO_THE_MODEL_IO     [FORWARD_STEP]',myThid)
723          CALL DO_THE_MODEL_IO( myTime, myIter, myThid )
724          CALL TIMER_STOP ('DO_THE_MODEL_IO     [FORWARD_STEP]',myThid)
725    #endif
726    
727    #ifdef HAVE_SIGREG
728          IF ( useSIGREG ) THEN
729            IF ( i_got_signal .GT. 0 ) THEN
730              CALL PACKAGES_WRITE_PICKUP(
731         I         .TRUE., myTime, myIter, myThid )
732    #ifndef ALLOW_OFFLINE
733              CALL WRITE_CHECKPOINT(
734         I         .TRUE., myTime, myIter, myThid )  
735    #endif
736              STOP 'Checkpoint completed -- killed by signal handler'
737            ENDIF
738          ENDIF
739    #endif
740    
741  C--   Save state for restarts  C--   Save state for restarts
742  C     Note:        CALL TIMER_START('WRITE_CHECKPOINT    [FORWARD_STEP]',myThid)
743  C     =====        CALL PACKAGES_WRITE_PICKUP(
744  C     Because of the ordering of the timestepping code and       I               .FALSE., myTime, myIter, myThid )
745  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 cg2d_x at "time-level" N+1/2.  
 C      where N = I+timeLevBase-1  
 C     Thus a checkpoint contains U.0000000000, GU.0000000001 and  
 C     cg2d_x.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)  
746        CALL WRITE_CHECKPOINT(        CALL WRITE_CHECKPOINT(
747       &        .FALSE., myCurrentTime, myCurrentIter, myThid )       I               .FALSE., myTime, myIter, myThid )  
748        CALL TIMER_STOP ('I/O (WRITE)        [FORWARD_STEP]',myThid)  #endif
749          CALL TIMER_STOP ('WRITE_CHECKPOINT    [FORWARD_STEP]',myThid)
750    
751        ENDIF  #ifdef ALLOW_DEBUG
752          IF ( debugLevel .GE. debLevB )
753         &    CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
754    #endif
755    
756        RETURN        RETURN
757        END        END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.132
  ViewVC Help
Powered by ViewVC 1.1.22