/[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.63 - (show annotations) (download)
Tue Oct 7 04:31:30 2003 UTC (20 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint51j_post, checkpoint51i_pre
Changes since 1.62: +8 -3 lines
fix Problem with bulk_force & therm_seaice.

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

  ViewVC Help
Powered by ViewVC 1.1.22