/[MITgcm]/MITgcm_contrib/AITCZ/code/forward_step.F
ViewVC logotype

Contents of /MITgcm_contrib/AITCZ/code/forward_step.F

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


Revision 1.1 - (show annotations) (download)
Wed Aug 20 15:24:59 2003 UTC (21 years, 11 months ago) by czaja
Branch: MAIN
CVS Tags: HEAD
Initial creation of Arnaud's simple coupled simulation.

1 C $Header: /u/u0/gcmpack/models/MITgcmUV/model/src/forward_step.F,v 1.23 2001/09/27 20:12:10 heimbach Exp $
2 C $Name: release1_beta1 $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: FORWARD_STEP
8 C !INTERFACE:
9 SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid )
10
11 C !DESCRIPTION: \bv
12 C *==================================================================
13 C | SUBROUTINE forward_step
14 C | o Run the ocean model and, optionally, evaluate a cost function.
15 C *==================================================================
16 C |
17 C | THE_MAIN_LOOP is the toplevel routine for the Tangent Linear and
18 C | Adjoint Model Compiler (TAMC). For this purpose the initialization
19 C | of the model was split into two parts. Those parameters that do
20 C | not depend on a specific model run are set in INITIALISE_FIXED,
21 C | whereas those that do depend on the specific realization are
22 C | initialized in INITIALISE_VARIA.
23 C |
24 C *==================================================================
25 C \ev
26
27 C !USES:
28 IMPLICIT NONE
29 C == Global variables ==
30 #include "SIZE.h"
31 #include "EEPARAMS.h"
32 #include "PARAMS.h"
33 #include "DYNVARS.h"
34 #include "FFIELDS.h"
35
36 #ifdef ALLOW_NONHYDROSTATIC
37 #include "CG3D.h"
38 #endif
39
40 #ifdef ALLOW_SHAP_FILT
41 #include "SHAP_FILT.h"
42 #endif
43
44
45 #ifdef ALLOW_AUTODIFF_TAMC
46 #include "tamc.h"
47 #include "ctrl.h"
48 #include "ctrl_dummy.h"
49 #include "cost.h"
50 #endif
51
52
53 C !LOCAL VARIABLES:
54 C == Routine arguments ==
55 C note: under the multi-threaded model myiter and
56 C mytime are local variables passed around as routine
57 C arguments. Although this is fiddly it saves the need to
58 C impose additional synchronisation points when they are
59 C updated.
60 C myiter - iteration counter for this thread
61 C mytime - time counter for this thread
62 C mythid - thread number for this instance of the routine.
63 integer iloop
64 integer mythid
65 integer myiter
66 _RL mytime
67 INTEGER bi,bj
68
69 CEOP
70
71
72 #ifdef ALLOW_AUTODIFF_TAMC
73 C-- Reset the model iteration counter and the model time.
74 myiter = nIter0 + (iloop-1)
75 mytime = startTime + float(iloop-1)*deltaTclock
76 #endif
77
78 #ifdef ALLOW_AUTODIFF_TAMC
79 C Include call to a dummy routine. Its adjoint will be
80 C called at the proper place in the adjoint code.
81 C The adjoint routine will print out adjoint values
82 C if requested. The location of the call is important,
83 C it has to be after the adjoint of the exchanges
84 C (DO_GTERM_BLOCKING_EXCHANGES).
85 CALL DUMMY_IN_STEPPING( myTime, myIter, myThid )
86 #endif
87
88 #ifdef EXACT_CONSERV
89 IF (exactConserv) THEN
90 C-- Update etaH(n+1) :
91 DO bj=myByLo(myThid),myByHi(myThid)
92 DO bi=myBxLo(myThid),myBxHi(myThid)
93 CALL CALC_EXACT_ETA( .FALSE., bi,bj, uVel,vVel,
94 I startTime, nIter0, myThid )
95 ENDDO
96 ENDDO
97 IF (implicDiv2Dflow .NE. 1. _d 0 )
98 & _EXCH_XY_R8(etaH, myThid )
99 ENDIF
100 #endif /* EXACT_CONSERV */
101
102 #ifdef NONLIN_FRSURF
103 C-- compute the future surface level thickness
104 C according to etaH(n+1)
105 IF ( nonlinFreeSurf.GT.0) THEN
106 CALL CALC_SURF_DR(etaH, myTime, myIter, myThid )
107 ENDIF
108 #endif /* NONLIN_FRSURF */
109
110 C-- Load forcing/external data fields.
111 #ifdef INCLUDE_EXTERNAL_FORCING_PACKAGE
112 C NOTE, that although the exf package is part of the
113 C distribution, it is not currently maintained, i.e.
114 C exf is disabled by default in genmake.
115 #ifdef ALLOW_BULKFORMULAE
116 CADJ STORE theta = comlev1, key = ikey_dynamics
117 #endif
118 CALL EXF_GETFORCING( mytime, myiter, mythid )
119 #else
120 CALL TIMER_START('EXTERNAL_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)
121 CALL EXTERNAL_FIELDS_LOAD( mytime, myiter, mythid )
122 CALL TIMER_STOP ('EXTERNAL_FIELDS_LOAD[THE_MAIN_LOOP]',mythid)
123 #endif
124
125 C-- Step forward fields and calculate time tendency terms.
126 CALL TIMER_START('THERMODYNAMICS [THE_MAIN_LOOP]',mythid)
127 CALL THERMODYNAMICS( myTime, myIter, myThid )
128 CALL TIMER_STOP ('THERMODYNAMICS [THE_MAIN_LOOP]',mythid)
129
130 #ifdef ALLOW_SHAP_FILT
131 IF (useSHAP_FILT .AND.
132 & staggerTimeStep .AND. shap_filt_TrStagg ) THEN
133 CALL TIMER_START('SHAP_FILT [THE_MAIN_LOOP]',myThid)
134 CALL SHAP_FILT_APPLY_TS( gT, gS, myTime, myIter, myThid )
135 CALL TIMER_STOP ('SHAP_FILT [THE_MAIN_LOOP]',myThid)
136 ENDIF
137 #endif
138
139
140 C-- Step forward fields and calculate time tendency terms.
141 IF ( momStepping ) THEN
142 CALL TIMER_START('DYNAMICS [THE_MAIN_LOOP]',mythid)
143 CALL DYNAMICS( myTime, myIter, myThid )
144 CALL TIMER_STOP ('DYNAMICS [THE_MAIN_LOOP]',mythid)
145 ENDIF
146
147 #ifndef ALLOW_AUTODIFF_TAMC
148 #ifdef ALLOW_NONHYDROSTATIC
149 C-- Step forward W field in N-H algorithm
150 IF ( momStepping .AND. nonHydrostatic ) THEN
151 CALL TIMER_START('CALC_GW [THE_MAIN_LOOP]',myThid)
152 CALL CALC_GW(myThid)
153 CALL TIMER_STOP ('CALC_GW [THE_MAIN_LOOP]',myThid)
154 ENDIF
155 #endif
156 #endif /* ALLOW_AUTODIFF_TAMC */
157
158 #ifdef NONLIN_FRSURF
159 C-- update hfacC,W,S and recip_hFac according to etaH(n+1) :
160 IF ( momStepping ) THEN
161 IF ( nonlinFreeSurf.GT.0) THEN
162 CALL UPDATE_SURF_DR( myTime, myIter, myThid )
163 ENDIF
164 C- update also CG2D matrix (and preconditioner)
165 IF ( nonlinFreeSurf.GT.2) THEN
166 CALL UPDATE_CG2D( myTime, myIter, myThid )
167 ENDIF
168 ENDIF
169 #endif
170
171 C-- This block has been moved to the_correction_step(). We have
172 C left this code here to indicate where the filters used to be
173 C in the algorithm before JMC moved them to after the pressure
174 C solver.
175 C
176 C modif omp: zonal filter before elliptic solver
177 #ifdef ALLOW_SHAP_FILT
178 IF (useSHAP_FILT .AND. shap_filt_uvStar) THEN
179 CALL TIMER_START('SHAP_FILT [THE_MAIN_LOOP]',myThid)
180 CALL SHAP_FILT_APPLY_UV( gUnm1,gVnm1, myTime,myIter,myThid )
181 IF (implicDiv2Dflow.LT.1.) THEN
182 C-- Explicit+Implicit part of the Barotropic Flow Divergence
183 C => Filtering of uVel,vVel is necessary
184 CALL SHAP_FILT_APPLY_UV( uVel,vVel, myTime,myIter,myThid )
185 ENDIF
186 CALL TIMER_STOP ('SHAP_FILT [THE_MAIN_LOOP]',myThid)
187 ENDIF
188 #endif
189
190 C
191 C modif omp: zonal filter before elliptic solver
192 #ifdef ALLOW_ZONAL_FILT
193 IF (zonal_filt_lat.LT.90.) THEN
194 CALL TIMER_START('ZONAL_FILT_APPLY [THE_MAIN_LOOP]',myThid)
195 CALL ZONAL_FILT_APPLY(
196 U gUnm1, gVnm1, gT, gS,
197 I myThid )
198 CALL TIMER_STOP ('ZONAL_FILT_APPLY [THE_MAIN_LOOP]',myThid)
199 ENDIF
200 #endif
201
202 C-- Solve elliptic equation(s).
203 C Two-dimensional only for conventional hydrostatic or
204 C three-dimensional for non-hydrostatic and/or IGW scheme.
205 IF ( momStepping ) THEN
206 CALL TIMER_START('SOLVE_FOR_PRESSURE [THE_MAIN_LOOP]',myThid)
207 CALL SOLVE_FOR_PRESSURE( myThid )
208 CALL TIMER_STOP ('SOLVE_FOR_PRESSURE [THE_MAIN_LOOP]',myThid)
209 ENDIF
210
211 C-- Correct divergence in flow field and cycle time-stepping
212 C arrays (for all fields) ; update time-counter
213 myIter = nIter0 + iLoop
214 myTime = startTime + deltaTClock * float(iLoop)
215 CALL TIMER_START('THE_CORRECTION_STEP [THE_MAIN_LOOP]',myThid)
216 CALL THE_CORRECTION_STEP(myTime, myIter, myThid)
217 CALL TIMER_STOP ('THE_CORRECTION_STEP [THE_MAIN_LOOP]',myThid)
218
219 C-- Do "blocking" sends and receives for tendency "overlap" terms
220 c CALL TIMER_START('BLOCKING_EXCHANGES [THE_MAIN_LOOP]',myThid)
221 c CALL DO_GTERM_BLOCKING_EXCHANGES( myThid )
222 c CALL TIMER_STOP ('BLOCKING_EXCHANGES [THE_MAIN_LOOP]',myThid)
223
224 C-- Do "blocking" sends and receives for field "overlap" terms
225 CALL TIMER_START('BLOCKING_EXCHANGES [THE_MAIN_LOOP]',myThid)
226 CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid )
227 CALL TIMER_STOP ('BLOCKING_EXCHANGES [THE_MAIN_LOOP]',myThid)
228
229 #ifdef ALLOW_FLT
230 C-- Calculate float trajectories
231 IF (useFLT) THEN
232 CALL TIMER_START('FLOATS [THE_MAIN_LOOP]',myThid)
233 CALL FLT_MAIN(myIter,myTime, myThid)
234 CALL TIMER_STOP ('FLOATS [THE_MAIN_LOOP]',myThid)
235 ENDIF
236 #endif
237
238 #ifndef EXCLUDE_MONITOR
239 C-- Check status of solution (statistics, cfl, etc...)
240 CALL MONITOR( myIter, myTime, myThid )
241 #endif /* EXCLUDE_MONITOR */
242
243 #ifndef ALLOW_AUTODIFF_TAMC
244 C-- Do IO if needed.
245 CALL TIMER_START('DO_THE_MODEL_IO [THE_MAIN_LOOP]',myThid)
246 CALL DO_THE_MODEL_IO( myTime, myIter, myThid )
247 CALL TIMER_STOP ('DO_THE_MODEL_IO [THE_MAIN_LOOP]',myThid)
248
249 C-- Save state for restarts
250 C Note: (jmc: is it still the case after ckp35 ?)
251 C =====
252 C Because of the ordering of the timestepping code and
253 C tendency term code at end of loop model arrays hold
254 C U,V,T,S at "time-level" N but gu, gv, gs, gt, guNM1,...
255 C at "time-level" N+1/2 (guNM1 at "time-level" N+1/2 is
256 C gu at "time-level" N-1/2) and etaN at "time-level" N+1/2.
257 C where N = I+timeLevBase-1
258 C Thus a checkpoint contains U.0000000000, GU.0000000001 and
259 C etaN.0000000001 in the indexing scheme used for the model
260 C "state" files. This example is referred to as a checkpoint
261 C at time level 1
262 CALL TIMER_START('WRITE_CHECKPOINT [THE_MAIN_LOOP]',myThid)
263 CALL WRITE_CHECKPOINT(
264 & .FALSE., myTime, myIter, myThid )
265 CALL TIMER_STOP ('WRITE_CHECKPOINT [THE_MAIN_LOOP]',myThid)
266
267 #endif /* ALLOW_AUTODIFF_TAMC */
268
269 END

  ViewVC Help
Powered by ViewVC 1.1.22