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

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

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


Revision 1.67 - (show annotations) (download)
Fri Oct 24 05:29:35 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51o_pre
Changes since 1.66: +1 -16 lines
 o undid all of the cp51 checkin pending some ongoing code cleanups
   and discussion

1 C $Header: /u/u3/gcmpack/MITgcm/model/src/forward_step.F,v 1.65 2003/10/23 04:41:40 edhill Exp $
2 C $Name: $
3
4 #include "AD_CONFIG.h"
5 #include "PACKAGES_CONFIG.h"
6 #include "CPP_OPTIONS.h"
7
8 cswdptr -- add --
9 #ifdef ALLOW_GCHEM
10 # include "GCHEM_OPTIONS.h"
11 #endif
12 cswdptr -- end add ---
13
14 CBOP
15 C !ROUTINE: FORWARD_STEP
16 C !INTERFACE:
17 SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid )
18
19 C !DESCRIPTION: \bv
20 C *==================================================================
21 C | SUBROUTINE forward_step
22 C | o Run the ocean model and, optionally, evaluate a cost function.
23 C *==================================================================
24 C |
25 C | THE_MAIN_LOOP is the toplevel routine for the Tangent Linear and
26 C | Adjoint Model Compiler (TAMC). For this purpose the initialization
27 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,
29 C | whereas those that do depend on the specific realization are
30 C | initialized in INITIALISE_VARIA.
31 C |
32 C *==================================================================
33 C \ev
34
35 C !USES:
36 IMPLICIT NONE
37 C == Global variables ==
38 #include "SIZE.h"
39 #include "EEPARAMS.h"
40 #include "PARAMS.h"
41 #include "DYNVARS.h"
42 #include "FFIELDS.h"
43
44 #ifdef ALLOW_NONHYDROSTATIC
45 #include "CG3D.h"
46 #endif
47
48 #ifdef ALLOW_SHAP_FILT
49 #include "SHAP_FILT.h"
50 #endif
51 #ifdef ALLOW_ZONAL_FILT
52 #include "ZONAL_FILT.h"
53 #endif
54
55 #ifdef ALLOW_AUTODIFF_TAMC
56 # include "tamc.h"
57 # include "ctrl.h"
58 # include "ctrl_dummy.h"
59 # include "cost.h"
60 # include "EOS.h"
61 # ifdef INCLUDE_EXTERNAL_FORCING_PACKAGE
62 # include "exf_fields.h"
63 # ifdef ALLOW_BULKFORMULAE
64 # include "exf_constants.h"
65 # endif
66 # endif
67 # ifdef ALLOW_OBCS
68 # include "OBCS.h"
69 # endif
70 # ifdef ALLOW_PTRACERS
71 # include "PTRACERS.h"
72 # endif
73 #endif /* ALLOW_AUTODIFF_TAMC */
74
75 C !LOCAL VARIABLES:
76 C == Routine arguments ==
77 C note: under the multi-threaded model myiter and
78 C mytime are local variables passed around as routine
79 C arguments. Although this is fiddly it saves the need to
80 C impose additional synchronisation points when they are
81 C updated.
82 C myiter - iteration counter for this thread
83 C mytime - time counter for this thread
84 C mythid - thread number for this instance of the routine.
85 integer iloop
86 integer mythid
87 integer myiter
88 _RL mytime
89 #ifdef ALLOW_BULK_FORCE
90 INTEGER bi,bj
91 #endif
92
93 CEOP
94
95 #ifndef DISABLE_DEBUGMODE
96 IF ( debugLevel .GE. debLevB )
97 & CALL DEBUG_ENTER('FORWARD_STEP',myThid)
98 #endif
99
100 #ifdef ALLOW_AUTODIFF_TAMC
101 C-- Reset the model iteration counter and the model time.
102 myiter = nIter0 + (iloop-1)
103 mytime = startTime + float(iloop-1)*deltaTclock
104 #endif
105
106 #if (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_AUTODIFF_MONITOR))
107 C Include call to a dummy routine. Its adjoint will be
108 C called at the proper place in the adjoint code.
109 C The adjoint routine will print out adjoint values
110 C if requested. The location of the call is important,
111 C it has to be after the adjoint of the exchanges
112 C (DO_GTERM_BLOCKING_EXCHANGES).
113 CALL DUMMY_IN_STEPPING( myTime, myIter, myThid )
114 cph I've commented this line since it may conflict with MITgcm's adjoint
115 cph However, need to check whether that's still consistent
116 cph with the ecco-branch (it should).
117 cph CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
118 #endif
119
120 #ifdef EXACT_CONSERV
121 IF (exactConserv) THEN
122 C-- Update etaH(n+1) :
123 CALL TIMER_START('UPDATE_ETAH [FORWARD_STEP]',mythid)
124 CALL UPDATE_ETAH( myTime, myIter, myThid )
125 CALL TIMER_STOP ('UPDATE_ETAH [FORWARD_STEP]',mythid)
126 ENDIF
127 #endif /* EXACT_CONSERV */
128
129 #ifdef NONLIN_FRSURF
130 IF ( select_rStar.NE.0 ) THEN
131 C-- r* : compute the future level thickness according to etaH(n+1)
132 CALL TIMER_START('CALC_R_STAR [FORWARD_STEP]',mythid)
133 CALL CALC_R_STAR(etaH, myTime, myIter, myThid )
134 CALL TIMER_STOP ('CALC_R_STAR [FORWARD_STEP]',mythid)
135 ELSEIF ( nonlinFreeSurf.GT.0) THEN
136 C-- compute the future surface level thickness according to etaH(n+1)
137 CALL TIMER_START('CALC_SURF_DR [FORWARD_STEP]',mythid)
138 CALL CALC_SURF_DR(etaH, myTime, myIter, myThid )
139 CALL TIMER_STOP ('CALC_SURF_DR [FORWARD_STEP]',mythid)
140 ENDIF
141 #endif /* NONLIN_FRSURF */
142
143 C-- Load forcing/external data fields.
144 #ifdef ALLOW_AUTODIFF_TAMC
145 c**************************************
146 #include "checkpoint_lev1_directives.h"
147 c**************************************
148 #endif
149
150
151 C-- Call external forcing package
152 #ifdef ALLOW_BULK_FORCE
153 IF ( useBulkforce ) THEN
154 #ifndef DISABLE_DEBUGMODE
155 IF ( debugLevel .GE. debLevB )
156 & CALL DEBUG_CALL('BULKF_FIELDS_LOAD',myThid)
157 #endif
158 CALL TIMER_START('BULKF_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)
159 CALL BULKF_FIELDS_LOAD( mytime, myiter, mythid )
160 CALL TIMER_STOP ('BULKF_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)
161 c calculate qnet and empmr (and wind stress)
162 DO bj=myByLo(myThid),myByHi(myThid)
163 DO bi=myBxLo(myThid),myBxHi(myThid)
164 CALL BULKF_FORCING( bi,bj, mytime, myiter, mythid )
165 ENDDO
166 ENDDO
167 c Update the tile edges.
168 _EXCH_XY_R8(Qnet, mythid)
169 _EXCH_XY_R8(EmPmR, mythid)
170 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
171 C _EXCH_XY_R8(fu , mythid)
172 C _EXCH_XY_R8(fv , mythid)
173 ELSE
174 #endif /* ALLOW_BULK_FORCE */
175
176 # ifdef ALLOW_EXF
177 C NOTE, that although the exf package is part of the
178 C distribution, it is not currently maintained, i.e.
179 C exf is disabled by default in genmake.
180 #ifndef DISABLE_DEBUGMODE
181 IF ( debugLevel .GE. debLevB )
182 & CALL DEBUG_CALL('EXF_GETFORCING',myThid)
183 #endif
184 CALL TIMER_START('EXF_GETFORCING [FORWARD_STEP]',mythid)
185 CALL EXF_GETFORCING( mytime, myiter, mythid )
186 CALL TIMER_STOP ('EXF_GETFORCING [FORWARD_STEP]',mythid)
187 # else /* ALLOW_EXF undef */
188 cph The following IF-statement creates an additional dependency
189 cph for the forcing fields requiring additional storing.
190 cph Therefore, the IF-statement will be put between CPP-OPTIONS,
191 cph assuming that ALLOW_SEAICE has not yet been differentiated.
192 # ifdef ALLOW_SEAICE
193 IF ( .NOT. useSEAICE ) THEN
194 # endif
195 #ifndef DISABLE_DEBUGMODE
196 IF ( debugLevel .GE. debLevB )
197 & CALL DEBUG_CALL('EXTERNAL_FIELDS_LOAD',myThid)
198 #endif
199 CALL TIMER_START('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
200 CALL EXTERNAL_FIELDS_LOAD( mytime, myiter, mythid )
201 CALL TIMER_STOP ('EXTERNAL_FIELDS_LOAD[FORWARD_STEP]',mythid)
202 # ifdef ALLOW_SEAICE
203 ENDIF
204 # endif
205 # endif /* ALLOW_EXF */
206 #ifdef ALLOW_BULK_FORCE
207 C-- end of if/else block useBulfforce --
208 ENDIF
209 #endif /* ALLOW_BULK_FORCE */
210
211 #if (defined (ALLOW_ADJOINT_RUN) || defined (ALLOW_TANGENTLINEAR_RUN))
212 c-- Add control vector for forcing and parameter fields
213 if ( myiter .EQ. nIter0 )
214 & CALL CTRL_MAP_FORCING (mythid)
215 #endif
216
217 # ifdef ALLOW_SEAICE
218 C-- Call sea ice model to compute forcing/external data fields. In
219 C addition to computing prognostic sea-ice variables and diagnosing the
220 C forcing/external data fields that drive the ocean model, SEAICE_MODEL
221 C also sets theta to the freezing point under sea-ice. The implied
222 C surface heat flux is then stored in variable surfaceTendencyTice,
223 C which is needed by KPP package (kpp_calc.F and kpp_transport_t.F)
224 C to diagnose surface buoyancy fluxes and for the non-local transport
225 C term. Because this call precedes model thermodynamics, temperature
226 C under sea-ice may not be "exactly" at the freezing point by the time
227 C theta is dumped or time-averaged.
228 IF ( useSEAICE ) THEN
229 #ifndef DISABLE_DEBUGMODE
230 IF ( debugLevel .GE. debLevB )
231 & CALL DEBUG_CALL('SEAICE_MODEL',myThid)
232 #endif
233 CALL TIMER_START('SEAICE_MODEL [FORWARD_STEP]',myThid)
234 CALL SEAICE_MODEL( myTime, myIter, myThid )
235 CALL TIMER_STOP ('SEAICE_MODEL [FORWARD_STEP]',myThid)
236 ENDIF
237 # endif /* ALLOW_SEAICE */
238
239 #ifdef ALLOW_AUTODIFF_TAMC
240 # ifdef ALLOW_PTRACERS
241 cph this replaces _bibj storing of ptracer within thermodynamics
242 CADJ STORE ptracer = comlev1, key = ikey_dynamics
243 # endif
244 #endif
245
246 #ifdef ALLOW_PTRACERS
247 # ifdef ALLOW_GCHEM
248 CALL GCHEM_FIELDS_LOAD( mytime, myiter, mythid )
249 # endif
250 #endif
251
252 C-- Step forward fields and calculate time tendency terms.
253 #ifndef DISABLE_DEBUGMODE
254 IF ( debugLevel .GE. debLevB )
255 & CALL DEBUG_CALL('THERMODYNAMICS',myThid)
256 #endif
257 CALL TIMER_START('THERMODYNAMICS [FORWARD_STEP]',mythid)
258 CALL THERMODYNAMICS( myTime, myIter, myThid )
259 CALL TIMER_STOP ('THERMODYNAMICS [FORWARD_STEP]',mythid)
260
261 C-- do exchanges (needed for DYNAMICS) when using stagger time-step :
262 #ifndef DISABLE_DEBUGMODE
263 IF ( debugLevel .GE. debLevB )
264 & CALL DEBUG_CALL('DO_STAGGER_FIELDS_EXCH.',myThid)
265 #endif
266 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
267 CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid )
268 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
269
270 #ifdef ALLOW_SHAP_FILT
271 IF (useSHAP_FILT .AND.
272 & staggerTimeStep .AND. shap_filt_TrStagg ) THEN
273 #ifndef DISABLE_DEBUGMODE
274 IF ( debugLevel .GE. debLevB )
275 & CALL DEBUG_CALL('SHAP_FILT_APPLY_TS',myThid)
276 #endif
277 CALL TIMER_START('SHAP_FILT [FORWARD_STEP]',myThid)
278 CALL SHAP_FILT_APPLY_TS(gT,gS,myTime+deltaT,myIter+1,myThid)
279 CALL TIMER_STOP ('SHAP_FILT [FORWARD_STEP]',myThid)
280 ENDIF
281 #endif
282 #ifdef ALLOW_ZONAL_FILT
283 IF (useZONAL_FILT .AND.
284 & staggerTimeStep .AND. zonal_filt_TrStagg ) THEN
285 #ifndef DISABLE_DEBUGMODE
286 IF ( debugLevel .GE. debLevB )
287 & CALL DEBUG_CALL('ZONAL_FILT_APPLY_TS',myThid)
288 #endif
289 CALL TIMER_START('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
290 CALL ZONAL_FILT_APPLY_TS( gT, gS, myThid )
291 CALL TIMER_STOP ('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
292 ENDIF
293 #endif
294
295 C-- Step forward fields and calculate time tendency terms.
296 #ifndef ALLOW_AUTODIFF_TAMC
297 IF ( momStepping ) THEN
298 #endif
299 #ifndef DISABLE_DEBUGMODE
300 IF ( debugLevel .GE. debLevB )
301 & CALL DEBUG_CALL('DYNAMICS',myThid)
302 #endif
303 CALL TIMER_START('DYNAMICS [FORWARD_STEP]',mythid)
304 CALL DYNAMICS( myTime, myIter, myThid )
305 CALL TIMER_STOP ('DYNAMICS [FORWARD_STEP]',mythid)
306 #ifndef ALLOW_AUTODIFF_TAMC
307 ENDIF
308 #endif
309
310 #ifdef ALLOW_NONHYDROSTATIC
311 C-- Step forward W field in N-H algorithm
312 IF ( momStepping .AND. nonHydrostatic ) THEN
313 #ifndef DISABLE_DEBUGMODE
314 IF ( debugLevel .GE. debLevB )
315 & CALL DEBUG_CALL('CALC_GW',myThid)
316 #endif
317 CALL TIMER_START('CALC_GW [FORWARD_STEP]',myThid)
318 CALL CALC_GW(myThid)
319 CALL TIMER_STOP ('CALC_GW [FORWARD_STEP]',myThid)
320 ENDIF
321 #endif
322
323 #ifdef NONLIN_FRSURF
324 C-- update hfacC,W,S and recip_hFac according to etaH(n+1) :
325 IF ( nonlinFreeSurf.GT.0) THEN
326 IF ( select_rStar.GT.0 ) THEN
327 CALL TIMER_START('UPDATE_R_STAR [FORWARD_STEP]',myThid)
328 CALL UPDATE_R_STAR( myTime, myIter, myThid )
329 CALL TIMER_STOP ('UPDATE_R_STAR [FORWARD_STEP]',myThid)
330 ELSE
331 CALL TIMER_START('UPDATE_SURF_DR [FORWARD_STEP]',myThid)
332 CALL UPDATE_SURF_DR( myTime, myIter, myThid )
333 CALL TIMER_STOP ('UPDATE_SURF_DR [FORWARD_STEP]',myThid)
334 ENDIF
335 ENDIF
336 C- update also CG2D matrix (and preconditioner)
337 IF ( momStepping .AND. nonlinFreeSurf.GT.2 ) THEN
338 CALL TIMER_START('UPDATE_CG2D [FORWARD_STEP]',myThid)
339 CALL UPDATE_CG2D( myTime, myIter, myThid )
340 CALL TIMER_STOP ('UPDATE_CG2D [FORWARD_STEP]',myThid)
341 ENDIF
342 #endif
343
344 C-- Apply Filters to u*,v* before SOLVE_FOR_PRESSURE
345 #ifdef ALLOW_SHAP_FILT
346 IF (useSHAP_FILT .AND. shap_filt_uvStar) THEN
347 CALL TIMER_START('SHAP_FILT [FORWARD_STEP]',myThid)
348 IF (implicDiv2Dflow.LT.1.) THEN
349 C-- Explicit+Implicit part of the Barotropic Flow Divergence
350 C => Filtering of uVel,vVel is necessary
351 CALL SHAP_FILT_APPLY_UV( uVel,vVel,
352 & myTime+deltaT, myIter+1, myThid )
353 ENDIF
354 CALL SHAP_FILT_APPLY_UV( gU,gV,myTime+deltaT,myIter+1,myThid)
355 CALL TIMER_STOP ('SHAP_FILT [FORWARD_STEP]',myThid)
356 ENDIF
357 #endif
358 #ifdef ALLOW_ZONAL_FILT
359 IF (useZONAL_FILT .AND. zonal_filt_uvStar) THEN
360 CALL TIMER_START('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
361 IF (implicDiv2Dflow.LT.1.) THEN
362 C-- Explicit+Implicit part of the Barotropic Flow Divergence
363 C => Filtering of uVel,vVel is necessary
364 CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
365 ENDIF
366 CALL ZONAL_FILT_APPLY_UV( gU, gV, myThid )
367 CALL TIMER_STOP ('ZONAL_FILT_APPLY [FORWARD_STEP]',myThid)
368 ENDIF
369 #endif
370
371 C-- Solve elliptic equation(s).
372 C Two-dimensional only for conventional hydrostatic or
373 C three-dimensional for non-hydrostatic and/or IGW scheme.
374 IF ( momStepping ) THEN
375 CALL TIMER_START('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
376 CALL SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
377 CALL TIMER_STOP ('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid)
378 ENDIF
379
380 #ifdef ALLOW_AUTODIFF_TAMC
381 cph This is needed because convective_adjustment calls
382 cph find_rho which may use pressure()
383 CADJ STORE totphihyd = comlev1, key = ikey_dynamics
384 #endif
385 C-- Correct divergence in flow field and cycle time-stepping
386 C arrays (for all fields) ; update time-counter
387 myIter = nIter0 + iLoop
388 myTime = startTime + deltaTClock * float(iLoop)
389 CALL TIMER_START('THE_CORRECTION_STEP [FORWARD_STEP]',myThid)
390 CALL THE_CORRECTION_STEP(myTime, myIter, myThid)
391 CALL TIMER_STOP ('THE_CORRECTION_STEP [FORWARD_STEP]',myThid)
392
393 C-- Do "blocking" sends and receives for tendency "overlap" terms
394 c CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
395 c CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )
396 c CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
397
398 C-- Do "blocking" sends and receives for field "overlap" terms
399 CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
400 CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
401 CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid)
402
403 cswdptr -- add for seperate timestepping of chemical/biological/forcing
404 cswdptr of ptracers ---
405 #ifdef ALLOW_GCHEM
406 ceh3 This is broken -- this ifdef should not be visible!
407 #ifdef PTRACERS_SEPARATE_FORCING
408 ceh3 needs an IF ( use GCHEM ) THEN
409 call GCHEM_FORCING_SEP( myTime,myIter,myThid )
410 #endif /* PTRACERS_SEPARATE_FORCING */
411 #endif /* ALLOW_GCHEM */
412 cswdptr -- end add ---
413
414
415 #ifdef ALLOW_FLT
416 C-- Calculate float trajectories
417 IF (useFLT) THEN
418 CALL TIMER_START('FLOATS [FORWARD_STEP]',myThid)
419 CALL FLT_MAIN(myIter,myTime, myThid)
420 CALL TIMER_STOP ('FLOATS [FORWARD_STEP]',myThid)
421 ENDIF
422 #endif
423
424 #ifndef EXCLUDE_MONITOR
425 C-- Check status of solution (statistics, cfl, etc...)
426 CALL TIMER_START('MONITOR [FORWARD_STEP]',myThid)
427 CALL MONITOR( myIter, myTime, myThid )
428 CALL TIMER_STOP ('MONITOR [FORWARD_STEP]',myThid)
429 #endif /* EXCLUDE_MONITOR */
430
431 C-- Do IO if needed.
432 CALL TIMER_START('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
433 CALL DO_THE_MODEL_IO( myTime, myIter, myThid )
434 CALL TIMER_STOP ('DO_THE_MODEL_IO [FORWARD_STEP]',myThid)
435
436 C-- Save state for restarts
437 C Note: (jmc: is it still the case after ckp35 ?)
438 C =====
439 C Because of the ordering of the timestepping code and
440 C tendency term code at end of loop model arrays hold
441 C U,V,T,S at "time-level" N but gu, gv, gs, gt, guNM1,...
442 C at "time-level" N+1/2 (guNM1 at "time-level" N+1/2 is
443 C gu at "time-level" N-1/2) and etaN at "time-level" N+1/2.
444 C where N = I+timeLevBase-1
445 C Thus a checkpoint contains U.0000000000, GU.0000000001 and
446 C etaN.0000000001 in the indexing scheme used for the model
447 C "state" files. This example is referred to as a checkpoint
448 C at time level 1
449 CALL TIMER_START('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
450 CALL WRITE_CHECKPOINT(
451 & .FALSE., myTime, myIter, myThid )
452 CALL TIMER_STOP ('WRITE_CHECKPOINT [FORWARD_STEP]',myThid)
453
454 #ifndef DISABLE_DEBUGMODE
455 IF ( debugLevel .GE. debLevB )
456 & CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
457 #endif
458
459 END

  ViewVC Help
Powered by ViewVC 1.1.22