/[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.16 by heimbach, Mon Aug 13 23:28:40 2001 UTC revision 1.125 by edhill, Sat Dec 3 08:30:32 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    #ifdef ALLOW_OFFLINE
8    # include "OFFLINE_OPTIONS.h"
9    #endif
10    #ifdef ALLOW_GMREDI
11    # include "GMREDI_OPTIONS.h"
12    #endif
13    
14    CBOP
15    C     !ROUTINE: FORWARD_STEP
16    C     !INTERFACE:
17        SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid )        SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid )
       IMPLICIT NONE  
18    
19  c     ==================================================================  C     !DESCRIPTION: \bv
20  c     SUBROUTINE forward_step  C     *==================================================================
21  c     ==================================================================  C     | SUBROUTINE forward_step
22  c  C     | o Run the ocean model and, optionally, evaluate a cost function.
23  c     o Run the ocean model and evaluate the specified cost function.  C     *==================================================================
24  c  C     |
25  c     *the_main_loop* is the toplevel routine for the Tangent Linear and  C     | THE_MAIN_LOOP is the toplevel routine for the Tangent Linear and
26  c     Adjoint Model Compiler (TAMC). For this purpose the initialization  C     | Adjoint Model Compiler (TAMC). For this purpose the initialization
27  c     of the model was split into two parts. Those parameters that do  C     | of the model was split into two parts. Those parameters that do
28  c     not depend on a specific model run are set in *initialise_fixed*,  C     | not depend on a specific model run are set in INITIALISE_FIXED,  
29  c     whereas those that do depend on the specific realization are  C     | whereas those that do depend on the specific realization are
30  c     initialized in *initialise_varia*.  C     | initialized in INITIALISE_VARIA.  
31    C     |
32    C     *==================================================================
33    C     \ev
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"
 #include "FFIELDS.h"  
 #include "TR1.h"  
42    
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  #endif
48    
49  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef HAVE_SIGREG
50  #include "tamc.h"  #include "SIGREG.h"
51  #include "ctrl.h"  #endif
52  #include "ctrl_dummy.h"  
53  #include "cost.h"  #ifdef ALLOW_SHAP_FILT
54    # include "SHAP_FILT.h"
55    #endif
56    #ifdef ALLOW_ZONAL_FILT
57    # include "ZONAL_FILT.h"
58  #endif  #endif
59    #ifdef COMPONENT_MODULE
60    # include "CPL_PARAMS.h"
61    #endif
62    
63    #ifdef ALLOW_AUTODIFF_TAMC
64    # include "FFIELDS.h"
65    
66  C     == routine arguments ==  # 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 ==
104  C     note: under the multi-threaded model myiter and  C     note: under the multi-threaded model myiter and
105  C           mytime are local variables passed around as routine  C           mytime are local variables passed around as routine
106  C           arguments. Although this is fiddly it saves the need to  C           arguments. Although this is fiddly it saves the need to
107  C           impose additional synchronisation points when they are  C           impose additional synchronisation points when they are
108  C           updated.  C           updated.
109  C     myiter - iteration counter for this thread  C     myIter - iteration counter for this thread
110  C     mytime - time counter for this thread  C     myTime - time counter for this thread
111  C     mythid - thread number for this instance of the routine.  C     myThid - thread number for this instance of the routine.
112        integer iloop        INTEGER iloop
113        integer mythid        INTEGER myThid
114        integer myiter        INTEGER myIter
115        _RL     mytime        _RL     myTime
116    
117    C     == Local variables ==
118    #ifdef COMPONENT_MODULE
119          INTEGER myItP1
120    #endif
121    CEOP
122    
123  C     == local variables ==  #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  #ifdef ALLOW_AUTODIFF_TAMC
135  C--     Reset the model iteration counter and the model time.  c**************************************
136          myiter = nIter0 + (iloop-1)  #include "checkpoint_lev1_directives.h"
137          mytime = startTime + float(iloop-1)*deltaTclock  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  #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  #ifdef ALLOW_AUTODIFF_TAMC
245  C       Include call to a dummy routine. Its adjoint will be  # ifdef ALLOW_PTRACERS
246  C       called at the proper place in the adjoint code.  cph this replaces _bibj storing of ptracer within thermodynamics
247  C       The adjoint routine will print out adjoint values  CADJ STORE ptracer  = comlev1, key = ikey_dynamics
248  C       if requested. The location of the call is important,  # endif
 C       it has to be after the adjoint of the exchanges  
 C       (DO_GTERM_BLOCKING_EXCHANGES).  
         CALL DUMMY_IN_STEPPING( myTime, myIter, myThid )  
249  #endif  #endif
250    
251  C--     Load forcing/external data fields.  #ifdef ALLOW_OFFLINE
252  #ifdef INCLUDE_EXTERNAL_FORCING_PACKAGE          call OFFLINE_FIELDS_LOAD( myTime, myIter, myThid )
253  C       NOTE, that although the exf package is part of the  #endif
254  C       distribution, it is not currently maintained, i.e.  
255  C       exf is disabled by default in genmake.  #ifdef ALLOW_PTRACERS
256          CALL EXF_GETFORCING( mytime, myiter, mythid )  # ifdef ALLOW_GCHEM
257  #else          IF ( useGCHEM ) THEN
258          CALL TIMER_START('EXTERNAL_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)  #ifdef ALLOW_DEBUG
259          CALL EXTERNAL_FIELDS_LOAD( mytime, myiter, mythid )           IF ( debugLevel .GE. debLevB )
260          CALL TIMER_STOP ('EXTERNAL_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)       &        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 COMPONENT_MODULE
268           IF ( useCoupler .AND. cpl_earlyExpImpCall ) THEN
269    C      Post coupling data that I export.
270    C      Read in coupling data that I import.
271             CALL TIMER_START('CPL_EXPORT-IMPORT  [FORWARD_STEP]',myThid)
272             CALL CPL_EXPORT_MY_DATA(       myIter, myTime, myThid )
273             CALL CPL_IMPORT_EXTERNAL_DATA( myIter, myTime, myThid )
274             CALL TIMER_STOP ('CPL_EXPORT-IMPORT  [FORWARD_STEP]',myThid)
275           ENDIF
276    #endif /* COMPONENT_MODULE */
277    
278    #ifdef ALLOW_EBM
279            IF ( useEBM ) THEN
280    # ifdef ALLOW_DEBUG
281             IF ( debugLevel .GE. debLevB )
282         &    CALL DEBUG_CALL('EBM',myThid)
283    # endif
284             CALL TIMER_START('EBM                [FORWARD_STEP]',mythid)
285             CALL EBM_DRIVER ( myTime, myIter, myThid )
286             CALL TIMER_STOP ('EBM                [FORWARD_STEP]',mythid)
287            ENDIF
288  #endif  #endif
289    
290  C--     Step forward fields and calculate time tendency terms.  C--     Step forward fields and calculate time tendency terms.
291          CALL TIMER_START('THERMODYNAMICS      [THE_MAIN_LOOP]',mythid)  
292    #ifdef ALLOW_DEBUG
293           IF ( debugLevel .GE. debLevB )
294         &    CALL DEBUG_CALL('DO_ATMOSPHERIC_PHYS',myThid)
295    #endif
296           CALL TIMER_START('DO_ATMOSPHERIC_PHYS [FORWARD_STEP]',mythid)
297           CALL DO_ATMOSPHERIC_PHYS( myTime, myIter, myThid )
298           CALL TIMER_STOP ('DO_ATMOSPHERIC_PHYS [FORWARD_STEP]',mythid)
299    
300    #ifdef ALLOW_AUTODIFF_TAMC
301    CADJ STORE theta              = comlev1, key = ikey_dynamics
302    CADJ STORE salt               = comlev1, key = ikey_dynamics
303    CADJ STORE totphihyd          = comlev1, key = ikey_dynamics
304    CADJ STORE surfaceforcingtice = comlev1, key = ikey_dynamics
305    # ifdef ALLOW_KPP
306    CADJ STORE uvel               = comlev1, key = ikey_dynamics
307    CADJ STORE vvel               = comlev1, key = ikey_dynamics
308    # endif
309    # ifdef EXACT_CONSERV
310    CADJ STORE empmr              = comlev1, key = ikey_dynamics
311    CADJ STORE pmepr              = comlev1, key = ikey_dynamics
312    # endif
313    #endif /* ALLOW_AUTODIFF_TAMC */
314    
315    #ifndef ALLOW_OFFLINE
316    #ifdef ALLOW_DEBUG
317           IF ( debugLevel .GE. debLevB )
318         &    CALL DEBUG_CALL('DO_OCEANIC_PHYS',myThid)
319    #endif
320           CALL TIMER_START('DO_OCEANIC_PHYS     [FORWARD_STEP]',mythid)
321           CALL DO_OCEANIC_PHYS( myTime, myIter, myThid )
322           CALL TIMER_STOP ('DO_OCEANIC_PHYS     [FORWARD_STEP]',mythid)
323    #endif
324    
325    #ifdef ALLOW_GCHEM
326    C     GCHEM package is an interface for any bio-geochemical or
327    C     ecosystem model you would like to include.
328    C     If GCHEM_SEPARATE_FORCING is not defined, you are
329    C     responsible for computing tendency terms for passive
330    C     tracers and storing them on a 3DxNumPtracers-array called
331    C     gchemTendency in GCHEM_CALC_TENDENCY. This tendency is then added
332    C     to gPtr in ptracers_forcing later-on.
333    C     If GCHEM_SEPARATE_FORCING is defined, you are reponsible for
334    C     UPDATING ptracers directly in GCHEM_FORCING_SEP. This amounts
335    C     to a completely separate time step that you have to implement
336    C     yourself (Eulerian seems to be fine in most cases).
337    CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
338    C     CAVEAT: Up to now, when GCHEM is turned on the field ptracerForcingSurf,
339    C     which is needed for KPP is not set properly. ptracerForcingSurf must
340    C     be treated differently depending on whether GCHEM_SEPARATE_FORCING
341    C     is define or not. TBD.
342    CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
343           IF ( useGCHEM ) THEN
344    #ifdef ALLOW_DEBUG
345            IF ( debugLevel .GE. debLevB )
346         &       CALL DEBUG_CALL('GCHEM_CALC_TENDENCY',myThid)
347    #endif
348            CALL TIMER_START('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid)
349            CALL GCHEM_CALC_TENDENCY( myTime, myIter, myThid )
350            CALL TIMER_STOP ('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid)
351           ENDIF
352    #endif /* ALLOW_GCHEM */
353    
354    #ifdef ALLOW_AUTODIFF_TAMC
355    cph needed to be moved here from do_oceanic_physics
356    cph to be visible down the road
357    c
358    CADJ STORE surfaceForcingS    = comlev1, key = ikey_dynamics
359    CADJ STORE surfaceForcingT    = comlev1, key = ikey_dynamics
360    CADJ STORE surfaceForcingTice = comlev1, key = ikey_dynamics
361    ctest(
362    CADJ STORE IVDConvCount       = comlev1, key = ikey_dynamics
363    ctest)
364    # ifdef ALLOW_PTRACERS
365    CADJ STORE surfaceForcingPtr  = comlev1, key = ikey_dynamics
366    # endif
367    c
368    # ifdef ALLOW_GMREDI
369    CADJ STORE Kwx                = comlev1, key = ikey_dynamics
370    CADJ STORE Kwy                = comlev1, key = ikey_dynamics
371    CADJ STORE Kwz                = comlev1, key = ikey_dynamics
372    #  ifdef GM_BOLUS_ADVEC
373    CADJ STORE GM_PsiX            = comlev1, key = ikey_dynamics
374    CADJ STORE GM_PsiY            = comlev1, key = ikey_dynamics
375    #  endif
376    # endif
377    c
378    # ifdef ALLOW_KPP
379    CADJ STORE KPPghat            = comlev1, key = ikey_dynamics
380    CADJ STORE KPPfrac            = comlev1, key = ikey_dynamics
381    CADJ STORE KPPdiffKzS         = comlev1, key = ikey_dynamics
382    CADJ STORE KPPdiffKzT         = comlev1, key = ikey_dynamics
383    # endif
384    #endif /* ALLOW_AUTODIFF_TAMC */
385    
386          IF ( .NOT.staggerTimeStep ) THEN
387    #ifdef ALLOW_DEBUG
388            IF ( debugLevel .GE. debLevB )
389         &    CALL DEBUG_CALL('THERMODYNAMICS',myThid)
390    #endif
391            CALL TIMER_START('THERMODYNAMICS      [FORWARD_STEP]',mythid)
392          CALL THERMODYNAMICS( myTime, myIter, myThid )          CALL THERMODYNAMICS( myTime, myIter, myThid )
393          CALL TIMER_STOP ('THERMODYNAMICS      [THE_MAIN_LOOP]',mythid)          CALL TIMER_STOP ('THERMODYNAMICS      [FORWARD_STEP]',mythid)
394    C--    if not staggerTimeStep: end
395          ENDIF
396    
397    #ifdef COMPONENT_MODULE
398           IF ( useCoupler .AND. .NOT.cpl_earlyExpImpCall ) THEN
399    C      Post coupling data that I export.
400    C      Read in coupling data that I import.
401             myItP1 = myIter + 1
402             CALL TIMER_START('CPL_EXPORT-IMPORT  [FORWARD_STEP]',myThid)
403             CALL CPL_EXPORT_MY_DATA(       myItP1, myTime, myThid )
404             CALL CPL_IMPORT_EXTERNAL_DATA( myItP1, myTime, myThid )
405             CALL TIMER_STOP ('CPL_EXPORT-IMPORT  [FORWARD_STEP]',myThid)
406    # ifndef ALLOW_AIM
407            IF ( useRealFreshWaterFlux ) THEN
408             CALL OCN_APPLY_IMPORT( .FALSE., myTime, myIter, myThid )
409            ENDIF
410    # endif
411           ENDIF
412    #endif /* COMPONENT_MODULE */
413    
414  C--     Step forward fields and calculate time tendency terms.  C--   Step forward fields and calculate time tendency terms.
415          CALL TIMER_START('DYNAMICS            [THE_MAIN_LOOP]',mythid)  #ifndef ALLOW_OFFLINE
416    #ifndef ALLOW_AUTODIFF_TAMC
417          IF ( momStepping ) THEN
418    #endif
419    #ifdef ALLOW_DEBUG
420            IF ( debugLevel .GE. debLevB )
421         &    CALL DEBUG_CALL('DYNAMICS',myThid)
422    #endif
423            CALL TIMER_START('DYNAMICS            [FORWARD_STEP]',mythid)
424          CALL DYNAMICS( myTime, myIter, myThid )          CALL DYNAMICS( myTime, myIter, myThid )
425          CALL TIMER_STOP ('DYNAMICS            [THE_MAIN_LOOP]',mythid)          CALL TIMER_STOP ('DYNAMICS            [FORWARD_STEP]',mythid)
   
426  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
427  #ifdef ALLOW_NONHYDROSTATIC        ENDIF
428  C--   Step forward W field in N-H algorithm  #endif
         IF ( nonHydrostatic ) THEN  
           CALL TIMER_START('CALC_GW          [THE_MAIN_LOOP]',myThid)  
           CALL CALC_GW(myThid)  
           CALL TIMER_STOP ('CALC_GW          [THE_MAIN_LOOP]',myThid)  
         ENDIF  
429  #endif  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
430    
431  C-- This block has been moved to the_correction_step(). We have  C--   Update time-counter
432  C   left this code here to indicate where the filters used to be        myIter = nIter0 + iLoop
433  C   in the algorithm before JMC moved them to after the pressure        myTime = startTime + deltaTClock * float(iLoop)
434  C   solver.  
435  C  #ifdef ALLOW_MNC
436  C#ifdef ALLOW_SHAP_FILT  C     Update the default next iter for MNC
437  CC--   Step forward all tiles, filter and exchange.        IF ( useMNC ) THEN
438  C      CALL TIMER_START('SHAP_FILT           [THE_MAIN_LOOP]',myThid)           CALL MNC_CW_CITER_SETG( 1, 1, -1, myIter , myThid )
439  C      CALL SHAP_FILT_APPLY(  
440  C     I                     gUnm1, gVnm1, gTnm1, gSnm1,  C        TODO: Logic should be added here so that users can specify, on
441  C     I                     myTime, myIter, myThid )  C        a per-citer-group basis, when it is time to update the
442  C      IF (implicDiv2Dflow.LT.1.) THEN  C        "current" (and not just the "next") iteration
443  CC--   Explicit+Implicit part of the Barotropic Flow Divergence  
444  CC      => Filtering of uVel,vVel is necessary  C        TODO: the following is just a temporary band-aid (mostly, for
445  C         CALL SHAP_FILT_UV( uVel, vVel, myTime, myThid )  C        Baylor) until someone writes a routine that better handles time
446  C      ENDIF  C        boundaries such as weeks, months, years, etc.
447  C      CALL TIMER_STOP ('SHAP_FILT           [THE_MAIN_LOOP]',myThid)           IF ( mnc_filefreq .GT. 0 ) THEN
448  C#endif             IF (DIFFERENT_MULTIPLE(mnc_filefreq,myTime,deltaTClock))
449  C       &          THEN
450  C#ifdef ALLOW_ZONAL_FILT               CALL MNC_CW_CITER_SETG( 1, 1, myIter, -1 , myThid )
451  C      IF (zonal_filt_lat.LT.90.) THEN             ENDIF
452  C        CALL TIMER_START('ZONAL_FILT_APPLY    [THE_MAIN_LOOP]',myThid)           ENDIF
453  C        CALL ZONAL_FILT_APPLY(         ENDIF
454  C     U           gUnm1, gVnm1, gTnm1, gSnm1,  #endif
455  C     I           myThid )  
456  C        CALL TIMER_STOP ('ZONAL_FILT_APPLY    [THE_MAIN_LOOP]',myThid)  C--   Update geometric factors:
457  C      ENDIF  #ifdef NONLIN_FRSURF
458  C#endif  C-    update hfacC,W,S and recip_hFac according to etaH(n+1) :
459          IF ( nonlinFreeSurf.GT.0) THEN
460           IF ( select_rStar.GT.0 ) THEN
461            CALL TIMER_START('UPDATE_R_STAR      [FORWARD_STEP]',myThid)
462            CALL UPDATE_R_STAR( myTime, myIter, myThid )
463            CALL TIMER_STOP ('UPDATE_R_STAR      [FORWARD_STEP]',myThid)
464           ELSE
465            CALL TIMER_START('UPDATE_SURF_DR     [FORWARD_STEP]',myThid)
466            CALL UPDATE_SURF_DR( myTime, myIter, myThid )
467            CALL TIMER_STOP ('UPDATE_SURF_DR     [FORWARD_STEP]',myThid)
468           ENDIF
469          ENDIF
470    C-    update also CG2D matrix (and preconditioner)
471          IF ( momStepping .AND. nonlinFreeSurf.GT.2 ) THEN
472            CALL TIMER_START('UPDATE_CG2D        [FORWARD_STEP]',myThid)
473            CALL UPDATE_CG2D( myTime, myIter, myThid )
474            CALL TIMER_STOP ('UPDATE_CG2D        [FORWARD_STEP]',myThid)
475          ENDIF
476    #endif
477    
478    C--   Apply Filters to u*,v* before SOLVE_FOR_PRESSURE
479    #ifdef ALLOW_SHAP_FILT
480          IF (useSHAP_FILT .AND. shap_filt_uvStar) THEN
481            CALL TIMER_START('SHAP_FILT_UV        [FORWARD_STEP]',myThid)
482            IF (implicDiv2Dflow.LT.1.) THEN
483    C--   Explicit+Implicit part of the Barotropic Flow Divergence
484    C      => Filtering of uVel,vVel is necessary
485              CALL SHAP_FILT_APPLY_UV( uVel,vVel,
486         &                             myTime, myIter, myThid )
487            ENDIF
488            CALL SHAP_FILT_APPLY_UV( gU,gV,myTime,myIter,myThid)
489            CALL TIMER_STOP ('SHAP_FILT_UV        [FORWARD_STEP]',myThid)
490          ENDIF
491    #endif
492    #ifdef ALLOW_ZONAL_FILT
493          IF (useZONAL_FILT .AND. zonal_filt_uvStar) THEN
494            CALL TIMER_START('ZONAL_FILT_UV       [FORWARD_STEP]',myThid)
495            IF (implicDiv2Dflow.LT.1.) THEN
496    C--   Explicit+Implicit part of the Barotropic Flow Divergence
497    C      => Filtering of uVel,vVel is necessary
498              CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
499            ENDIF
500            CALL ZONAL_FILT_APPLY_UV( gU, gV, myThid )
501            CALL TIMER_STOP ('ZONAL_FILT_UV       [FORWARD_STEP]',myThid)
502          ENDIF
503    #endif  
504    
505  C--   Solve elliptic equation(s).  C--   Solve elliptic equation(s).
506  C     Two-dimensional only for conventional hydrostatic or  C     Two-dimensional only for conventional hydrostatic or
507  C     three-dimensional for non-hydrostatic and/or IGW scheme.  C     three-dimensional for non-hydrostatic and/or IGW scheme.
508        CALL TIMER_START('SOLVE_FOR_PRESSURE  [THE_MAIN_LOOP]',myThid)  #ifndef ALLOW_OFFLINE
509        CALL SOLVE_FOR_PRESSURE( myThid )        IF ( momStepping ) THEN
510        CALL TIMER_STOP ('SOLVE_FOR_PRESSURE  [THE_MAIN_LOOP]',myThid)          CALL TIMER_START('SOLVE_FOR_PRESSURE  [FORWARD_STEP]',myThid)
511            CALL SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
512            CALL TIMER_STOP ('SOLVE_FOR_PRESSURE  [FORWARD_STEP]',myThid)
513          ENDIF
514    #endif
515    
516  C--   Correct divergence in flow field and cycle time-stepping  C--   Correct divergence in flow field and cycle time-stepping momentum
517  C     arrays (for all fields) ; update time-counter  c     IF ( momStepping ) THEN
518        myIter = nIter0 + iLoop  #ifndef ALLOW_OFFLINE
519        myTime = startTime + deltaTClock * float(iLoop)          CALL TIMER_START('MOM_CORRECTION_STEP [FORWARD_STEP]',myThid)
520        CALL TIMER_START('THE_CORRECTION_STEP [THE_MAIN_LOOP]',myThid)          CALL MOMENTUM_CORRECTION_STEP(myTime, myIter, myThid)
521        CALL THE_CORRECTION_STEP(myTime, myIter, myThid)          CALL TIMER_STOP ('MOM_CORRECTION_STEP [FORWARD_STEP]',myThid)
522        CALL TIMER_STOP ('THE_CORRECTION_STEP [THE_MAIN_LOOP]',myThid)  #endif
523    c     ENDIF
524    
525    #ifdef EXACT_CONSERV
526          IF (exactConserv) THEN
527    C--   Update etaH(n+1) :
528            CALL TIMER_START('UPDATE_ETAH         [FORWARD_STEP]',mythid)
529            CALL UPDATE_ETAH( myTime, myIter, myThid )
530            CALL TIMER_STOP ('UPDATE_ETAH         [FORWARD_STEP]',mythid)
531          ENDIF
532    #endif /* EXACT_CONSERV */
533    
534    #ifdef NONLIN_FRSURF
535          IF ( select_rStar.NE.0 ) THEN
536    C--   r* : compute the future level thickness according to etaH(n+1)
537            CALL TIMER_START('CALC_R_STAR       [FORWARD_STEP]',mythid)
538            CALL CALC_R_STAR(etaH, myTime, myIter, myThid )
539            CALL TIMER_STOP ('CALC_R_STAR       [FORWARD_STEP]',mythid)
540          ELSEIF ( nonlinFreeSurf.GT.0) THEN
541    C--   compute the future surface level thickness according to etaH(n+1)
542            CALL TIMER_START('CALC_SURF_DR      [FORWARD_STEP]',mythid)
543            CALL CALC_SURF_DR(etaH, myTime, myIter, myThid )
544            CALL TIMER_STOP ('CALC_SURF_DR      [FORWARD_STEP]',mythid)
545          ENDIF
546    #endif /* NONLIN_FRSURF */
547    
548    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
549          IF ( staggerTimeStep ) THEN
550    C--   do exchanges of U,V (needed for multiDim) when using stagger time-step :
551    #ifdef ALLOW_DEBUG
552            IF ( debugLevel .GE. debLevB )
553         &    CALL DEBUG_CALL('DO_STAGGER_FIELDS_EXCH.',myThid)
554    #endif
555            CALL TIMER_START('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
556            CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
557            CALL TIMER_STOP ('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
558    
559    C--   State-variables diagnostics
560           IF ( usediagnostics ) THEN
561            CALL TIMER_START('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
562            CALL DO_STATEVARS_DIAGS( myTime, 1, myIter, myThid )
563            CALL TIMER_STOP ('DO_STATEVARS_DIAGS  [FORWARD_STEP]',myThid)
564           ENDIF
565    
566    #ifdef ALLOW_DEBUG
567            IF ( debugLevel .GE. debLevB )
568         &    CALL DEBUG_CALL('THERMODYNAMICS',myThid)
569    #endif
570            CALL TIMER_START('THERMODYNAMICS      [FORWARD_STEP]',mythid)
571            CALL THERMODYNAMICS( myTime, myIter, myThid )
572            CALL TIMER_STOP ('THERMODYNAMICS      [FORWARD_STEP]',mythid)
573    
574    C--    if staggerTimeStep: end
575          ENDIF
576    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
577    
578    #ifdef ALLOW_AUTODIFF_TAMC
579    cph This is needed because convective_adjustment calls
580    cph find_rho which may use pressure()
581    CADJ STORE totphihyd  = comlev1, key = ikey_dynamics
582    #endif
583    C--   Cycle time-stepping Tracers arrays (T,S,+pTracers)
584            CALL TIMER_START('TRC_CORRECTION_STEP [FORWARD_STEP]',myThid)
585            CALL TRACERS_CORRECTION_STEP(myTime, myIter, myThid)
586            CALL TIMER_STOP ('TRC_CORRECTION_STEP [FORWARD_STEP]',myThid)
587    
588    #ifdef ALLOW_GCHEM
589    C     Add separate timestepping of chemical/biological/forcing
590    C     of ptracers here in GCHEM_FORCING_SEP
591            IF ( useGCHEM ) THEN
592    #ifdef ALLOW_DEBUG
593             IF ( debugLevel .GE. debLevB )
594         &        CALL DEBUG_CALL('GCHEM_FORCING_SEP',myThid)
595    #endif /* ALLOW_DEBUG */
596             CALL TIMER_START('GCHEM_FORCING_SEP  [FORWARD_STEP]',myThid)
597             CALL GCHEM_FORCING_SEP( myTime,myIter,myThid )
598             CALL TIMER_STOP ('GCHEM_FORCING_SEP  [FORWARD_STEP]',myThid)
599            ENDIF  
600    #endif /* ALLOW_GCHEM */
601    
602  C--   Do "blocking" sends and receives for tendency "overlap" terms  C--   Do "blocking" sends and receives for tendency "overlap" terms
603  c     CALL TIMER_START('BLOCKING_EXCHANGES  [THE_MAIN_LOOP]',myThid)  c     CALL TIMER_START('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
604  c     CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )  c     CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )
605  c     CALL TIMER_STOP ('BLOCKING_EXCHANGES  [THE_MAIN_LOOP]',myThid)  c     CALL TIMER_STOP ('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
606    
607  C--   Do "blocking" sends and receives for field "overlap" terms  C--   Do "blocking" sends and receives for field "overlap" terms
608        CALL TIMER_START('BLOCKING_EXCHANGES  [THE_MAIN_LOOP]',myThid)        CALL TIMER_START('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
609        CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )        CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
610        CALL TIMER_STOP ('BLOCKING_EXCHANGES  [THE_MAIN_LOOP]',myThid)        CALL TIMER_STOP ('BLOCKING_EXCHANGES  [FORWARD_STEP]',myThid)
611    
612    
613    C AMM
614    #ifdef ALLOW_GRIDALT
615            if (useGRIDALT) then
616             CALL GRIDALT_UPDATE(myThid)
617            endif
618    #endif
619    C AMM
620    
621    C AMM
622    #ifdef ALLOW_FIZHI
623            if( useFIZHI) then
624             CALL TIMER_START('FIZHI               [FORWARD_STEP]',mythid)
625             CALL STEP_FIZHI_CORR ( myTime, myIter, myThid )
626             CALL TIMER_STOP('FIZHI               [FORWARD_STEP]',mythid)
627            endif
628    #endif
629    C AMM
630    
631    #ifdef ALLOW_FLT
632    C--   Calculate float trajectories
633          IF (useFLT) THEN
634            CALL TIMER_START('FLOATS            [FORWARD_STEP]',myThid)
635            CALL FLT_MAIN(myIter,myTime, myThid)
636            CALL TIMER_STOP ('FLOATS            [FORWARD_STEP]',myThid)
637          ENDIF
638    #endif
639    
640    C--   State-variables time-averaging
641          CALL TIMER_START('DO_STATEVARS_TAVE   [FORWARD_STEP]',myThid)
642          CALL DO_STATEVARS_TAVE( myTime, myIter, myThid )
643          CALL TIMER_STOP ('DO_STATEVARS_TAVE   [FORWARD_STEP]',myThid)
644    
645  #ifndef EXCLUDE_MONITOR  #ifndef ALLOW_OFFLINE
646    #ifdef ALLOW_MONITOR
647  C--   Check status of solution (statistics, cfl, etc...)  C--   Check status of solution (statistics, cfl, etc...)
648          CALL TIMER_START('MONITOR             [FORWARD_STEP]',myThid)
649        CALL MONITOR( myIter, myTime, myThid )        CALL MONITOR( myIter, myTime, myThid )
650  #endif /* EXCLUDE_MONITOR */        CALL TIMER_STOP ('MONITOR             [FORWARD_STEP]',myThid)
651    #endif /* ALLOW_MONITOR */
652    #endif
653    
654    #ifdef ALLOW_COST
655    C--     compare model with data and compute cost function
656    C--     this is done after exchanges to allow interpolation
657          CALL TIMER_START('COST_TILE           [FORWARD_STEP]',myThid)
658          CALL COST_TILE  ( mytime, myiter, myThid )
659          CALL TIMER_STOP ('COST_TILE           [FORWARD_STEP]',myThid)
660    #endif
661    
 #ifndef ALLOW_AUTODIFF_TAMC  
662  C--   Do IO if needed.  C--   Do IO if needed.
663        CALL TIMER_START('DO_THE_MODEL_IO     [THE_MAIN_LOOP]',myThid)  #ifdef ALLOW_OFFLINE
664          CALL TIMER_START('OFFLINE_MODEL_IO    [FORWARD_STEP]',myThid)
665          CALL OFFLINE_MODEL_IO( myTime, myIter, myThid )
666          CALL TIMER_STOP ('OFFLINE_MODEL_IO    [FORWARD_STEP]',myThid)
667    #else
668          CALL TIMER_START('DO_THE_MODEL_IO     [FORWARD_STEP]',myThid)
669        CALL DO_THE_MODEL_IO( myTime, myIter, myThid )        CALL DO_THE_MODEL_IO( myTime, myIter, myThid )
670        CALL TIMER_STOP ('DO_THE_MODEL_IO     [THE_MAIN_LOOP]',myThid)        CALL TIMER_STOP ('DO_THE_MODEL_IO     [FORWARD_STEP]',myThid)
671    #endif
672    
673    #ifdef HAVE_SIGREG
674          IF ( useSIGREG ) THEN
675            IF ( i_got_signal .GT. 0 ) THEN
676              CALL PACKAGES_WRITE_PICKUP(
677         I         .TRUE., myTime, myIter, myThid )
678    #ifndef ALLOW_OFFLINE
679              CALL WRITE_CHECKPOINT(
680         I         .TRUE., myTime, myIter, myThid )  
681    #endif
682              STOP 'Checkpoint completed -- killed by signal handler'
683            ENDIF
684          ENDIF
685    #endif
686    
687  C--   Save state for restarts  C--   Save state for restarts
688  C     Note:    (jmc: is it still the case after ckp35 ?)        CALL TIMER_START('WRITE_CHECKPOINT    [FORWARD_STEP]',myThid)
689  C     =====        CALL PACKAGES_WRITE_PICKUP(
690  C     Because of the ordering of the timestepping code and       I               .FALSE., myTime, myIter, myThid )
691  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('WRITE_CHECKPOINT    [THE_MAIN_LOOP]',myThid)  
692        CALL WRITE_CHECKPOINT(        CALL WRITE_CHECKPOINT(
693       &        .FALSE., myTime, myIter, myThid )       I               .FALSE., myTime, myIter, myThid )  
694        CALL TIMER_STOP ('WRITE_CHECKPOINT    [THE_MAIN_LOOP]',myThid)  #endif
695          CALL TIMER_STOP ('WRITE_CHECKPOINT    [FORWARD_STEP]',myThid)
696    
697  #endif /* ALLOW_AUTODIFF_TAMC */  #ifdef ALLOW_DEBUG
698          IF ( debugLevel .GE. debLevB )
699         &    CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
700    #endif
701    
702          RETURN
703        END        END
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.125
  ViewVC Help
Powered by ViewVC 1.1.22