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

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

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


Revision 1.7 - (show annotations) (download)
Fri Apr 4 20:54:11 2014 UTC (10 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64w, checkpoint64v
Changes since 1.6: +12 -9 lines
- Start to include explicitly AUTODIFF_OPTIONS.h, COST_OPTIONS.h,
  and CTRL_OPTIONS.h in src files (to enable to skip the ECCO_CPPOPTIONS.h)
  For now, only in pkgs used in verification/hs94.1x64x5.
- Replace ALLOW_AUTODIFF_TAMC by ALLOW_AUTODIFF (except for tape/storage
  which are specific to TAF/TAMC).

1 C $Header: /u/gcmpack/MITgcm/model/src/temp_integrate.F,v 1.6 2014/01/13 20:49:54 gforget Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6 #ifdef ALLOW_AUTODIFF
7 # include "AUTODIFF_OPTIONS.h"
8 #endif
9 #ifdef ALLOW_GENERIC_ADVDIFF
10 # include "GAD_OPTIONS.h"
11 #endif
12
13 CBOP
14 C !ROUTINE: TEMP_INTEGRATE
15 C !INTERFACE:
16 SUBROUTINE TEMP_INTEGRATE(
17 I bi, bj, recip_hFac,
18 I uFld, vFld, wFld,
19 U KappaRk,
20 I myTime, myIter, myThid )
21 C !DESCRIPTION: \bv
22 C *==========================================================*
23 C | SUBROUTINE TEMP_INTEGRATE
24 C | o Calculate tendency for temperature
25 C | and integrates forward in time.
26 C *==========================================================*
27 C | A procedure called EXTERNAL_FORCING_T is called from
28 C | here. These procedures can be used to add per problem
29 C | heat flux source terms.
30 C | Note: Although it is slightly counter-intuitive the
31 C | EXTERNAL_FORCING routine is not the place to put
32 C | file I/O. Instead files that are required to
33 C | calculate the external source terms are generally
34 C | read during the model main loop. This makes the
35 C | logistics of multi-processing simpler and also
36 C | makes the adjoint generation simpler. It also
37 C | allows for I/O to overlap computation where that
38 C | is supported by hardware.
39 C | Aside from the problem specific term the code here
40 C | forms the tendency terms due to advection and mixing
41 C | The baseline implementation here uses a centered
42 C | difference form for the advection term and a tensorial
43 C | divergence of a flux form for the diffusive term. The
44 C | diffusive term is formulated so that isopycnal mixing
45 C | and GM-style subgrid-scale terms can be incorporated by
46 C | simply setting the diffusion tensor terms appropriately.
47 C *==========================================================*
48 C \ev
49
50 C !USES:
51 IMPLICIT NONE
52 C == GLobal variables ==
53 #include "SIZE.h"
54 #include "EEPARAMS.h"
55 #include "PARAMS.h"
56 #include "GRID.h"
57 #include "DYNVARS.h"
58 #include "RESTART.h"
59 #ifdef ALLOW_GENERIC_ADVDIFF
60 # include "GAD.h"
61 # include "GAD_SOM_VARS.h"
62 #endif
63 #ifdef ALLOW_TIMEAVE
64 # include "TIMEAVE_STATV.h"
65 #endif
66 #ifdef ALLOW_AUTODIFF
67 # include "tamc.h"
68 # include "tamc_keys.h"
69 #endif
70
71 C !INPUT/OUTPUT PARAMETERS:
72 C == Routine arguments ==
73 C bi, bj, :: tile indices
74 C recip_hFac :: reciprocal of cell open-depth factor (@ next iter)
75 C uFld,vFld :: Local copy of horizontal velocity field
76 C wFld :: Local copy of vertical velocity field
77 C KappaRk :: Vertical diffusion for Tempertature
78 C myTime :: current time
79 C myIter :: current iteration number
80 C myThid :: my Thread Id. number
81 INTEGER bi, bj
82 _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
84 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
85 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86 _RL KappaRk (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87 _RL myTime
88 INTEGER myIter
89 INTEGER myThid
90 CEOP
91
92 #ifdef ALLOW_GENERIC_ADVDIFF
93 C !LOCAL VARIABLES:
94 C iMin, iMax :: 1rst index loop range
95 C jMin, jMax :: 2nd index loop range
96 C k :: vertical index
97 C kM1 :: =k-1 for k>1, =1 for k=1
98 C kUp :: index into 2 1/2D array, toggles between 1|2
99 C kDown :: index into 2 1/2D array, toggles between 2|1
100 C xA :: Tracer cell face area normal to X
101 C yA :: Tracer cell face area normal to X
102 C maskUp :: Land/water mask for Wvel points (interface k)
103 C uTrans :: Zonal volume transport through cell face
104 C vTrans :: Meridional volume transport through cell face
105 C rTrans :: Vertical volume transport at interface k
106 C rTransKp :: Vertical volume transport at inteface k+1
107 C fZon :: Flux of temperature (T) in the zonal direction
108 C fMer :: Flux of temperature (T) in the meridional direction
109 C fVer :: Flux of temperature (T) in the vertical direction
110 C at the upper(U) and lower(D) faces of a cell.
111 C useVariableK :: T when vertical diffusion is not constant
112 INTEGER iMin, iMax, jMin, jMax
113 INTEGER i, j, k
114 INTEGER kUp, kDown, kM1
115 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121 _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124 _RL fVer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
125 _RL gt_AB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126 LOGICAL calcAdvection
127 INTEGER iterNb
128 #ifdef ALLOW_ADAMSBASHFORTH_3
129 INTEGER m1, m2
130 #endif
131 #ifdef ALLOW_TIMEAVE
132 LOGICAL useVariableK
133 #endif
134
135 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
136
137 C- Loop ranges for daughter routines
138 iMin = 1-OLx+2
139 iMax = sNx+OLx-1
140 jMin = 1-OLy+2
141 jMax = sNy+OLy-1
142
143 iterNb = myIter
144 IF (staggerTimeStep) iterNb = myIter - 1
145
146 #ifdef ALLOW_AUTODIFF_TAMC
147 act1 = bi - myBxLo(myThid)
148 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
149 act2 = bj - myByLo(myThid)
150 max2 = myByHi(myThid) - myByLo(myThid) + 1
151 act3 = myThid - 1
152 max3 = nTx*nTy
153 act4 = ikey_dynamics - 1
154 itdkey = (act1 + 1) + act2*max1
155 & + act3*max1*max2
156 & + act4*max1*max2*max3
157 #endif /* ALLOW_AUTODIFF_TAMC */
158
159 C- Tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
160 DO k=1,Nr
161 DO j=1-OLy,sNy+OLy
162 DO i=1-OLx,sNx+OLx
163 gT(i,j,k,bi,bj) = 0. _d 0
164 ENDDO
165 ENDDO
166 ENDDO
167 DO j=1-OLy,sNy+OLy
168 DO i=1-OLx,sNx+OLx
169 fVer(i,j,1) = 0. _d 0
170 fVer(i,j,2) = 0. _d 0
171 ENDDO
172 ENDDO
173 #ifdef ALLOW_AUTODIFF
174 DO k=1,Nr
175 DO j=1-OLy,sNy+OLy
176 DO i=1-OLx,sNx+OLx
177 kappaRk(i,j,k) = 0. _d 0
178 ENDDO
179 ENDDO
180 ENDDO
181 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
182 CADJ STORE wFld(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
183 #endif /* ALLOW_AUTODIFF */
184
185 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
186 CALL CALC_3D_DIFFUSIVITY(
187 I bi, bj, iMin, iMax, jMin, jMax,
188 I GAD_TEMPERATURE, useGMredi, useKPP,
189 O kappaRk,
190 I myThid )
191 #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
192
193 #ifndef DISABLE_MULTIDIM_ADVECTION
194 C-- Some advection schemes are better calculated using a multi-dimensional
195 C method in the absence of any other terms and, if used, is done here.
196 C
197 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
198 C The default is to use multi-dimensinal advection for non-linear advection
199 C schemes. However, for the sake of efficiency of the adjoint it is necessary
200 C to be able to exclude this scheme to avoid excessive storage and
201 C recomputation. It *is* differentiable, if you need it.
202 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
203 C disable this section of code.
204 #ifdef GAD_ALLOW_TS_SOM_ADV
205 # ifdef ALLOW_AUTODIFF_TAMC
206 CADJ STORE som_T = comlev1_bibj, key=itdkey, byte=isbyte
207 # endif
208 IF ( tempSOM_Advection ) THEN
209 # ifdef ALLOW_DEBUG
210 IF (debugMode) CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
211 # endif
212 CALL GAD_SOM_ADVECT(
213 I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
214 I GAD_TEMPERATURE, dTtracerLev,
215 I uFld, vFld, wFld, theta,
216 U som_T,
217 O gT,
218 I bi, bj, myTime, myIter, myThid )
219 ELSEIF (tempMultiDimAdvec) THEN
220 #else /* GAD_ALLOW_TS_SOM_ADV */
221 IF (tempMultiDimAdvec) THEN
222 #endif /* GAD_ALLOW_TS_SOM_ADV */
223 # ifdef ALLOW_DEBUG
224 IF (debugMode) CALL DEBUG_CALL('GAD_ADVECTION',myThid)
225 # endif
226 CALL GAD_ADVECTION(
227 I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
228 I GAD_TEMPERATURE, dTtracerLev,
229 I uFld, vFld, wFld, theta,
230 O gT,
231 I bi, bj, myTime, myIter, myThid )
232 ENDIF
233 #endif /* DISABLE_MULTIDIM_ADVECTION */
234
235 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
236
237 C- Start vertical index (k) loop (Nr:1)
238 calcAdvection = tempAdvection .AND. .NOT.tempMultiDimAdvec
239 DO k=Nr,1,-1
240 #ifdef ALLOW_AUTODIFF_TAMC
241 kkey = (itdkey-1)*Nr + k
242 #endif /* ALLOW_AUTODIFF_TAMC */
243 kM1 = MAX(1,k-1)
244 kUp = 1+MOD(k+1,2)
245 kDown= 1+MOD(k,2)
246
247 #ifdef ALLOW_AUTODIFF_TAMC
248 CADJ STORE fVer(:,:,:) = comlev1_bibj_k, key=kkey,
249 CADJ & byte=isbyte, kind = isbyte
250 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
251 CADJ & byte=isbyte, kind = isbyte
252 # ifdef ALLOW_ADAMSBASHFORTH_3
253 CADJ STORE gtNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey,
254 CADJ & byte=isbyte, kind = isbyte
255 CADJ STORE gtNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey,
256 CADJ & byte=isbyte, kind = isbyte
257 # else
258 CADJ STORE gtNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
259 CADJ & byte=isbyte, kind = isbyte
260 # endif
261 #endif /* ALLOW_AUTODIFF_TAMC */
262 CALL CALC_ADV_FLOW(
263 I uFld, vFld, wFld,
264 U rTrans,
265 O uTrans, vTrans, rTransKp,
266 O maskUp, xA, yA,
267 I k, bi, bj, myThid )
268
269 #ifdef ALLOW_ADAMSBASHFORTH_3
270 m1 = 1 + MOD(iterNb+1,2)
271 m2 = 1 + MOD( iterNb ,2)
272 CALL GAD_CALC_RHS(
273 I bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
274 I xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
275 I vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
276 I uTrans, vTrans, rTrans, rTransKp,
277 I diffKhT, diffK4T, KappaRk(1-OLx,1-OLy,k), diffKr4T,
278 I gtNm(1-OLx,1-OLy,1,1,1,m2), theta, dTtracerLev,
279 I GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
280 I calcAdvection, tempImplVertAdv, AdamsBashforth_T,
281 I tempVertDiff4, useGMRedi, useKPP,
282 O fZon, fMer,
283 U fVer, gT,
284 I myTime, myIter, myThid )
285 #else /* ALLOW_ADAMSBASHFORTH_3 */
286 CALL GAD_CALC_RHS(
287 I bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
288 I xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
289 I vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
290 I uTrans, vTrans, rTrans, rTransKp,
291 I diffKhT, diffK4T, KappaRk(1-OLx,1-OLy,k), diffKr4T,
292 I gtNm1, theta, dTtracerLev,
293 I GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
294 I calcAdvection, tempImplVertAdv, AdamsBashforth_T,
295 I tempVertDiff4, useGMRedi, useKPP,
296 O fZon, fMer,
297 U fVer, gT,
298 I myTime, myIter, myThid )
299 #endif
300
301 C-- External thermal forcing term(s) inside Adams-Bashforth:
302 IF ( tempForcing .AND. tracForcingOutAB.NE.1 )
303 & CALL EXTERNAL_FORCING_T(
304 I iMin, iMax, jMin, jMax, bi, bj, k,
305 I myTime, myThid )
306
307 IF ( AdamsBashforthGt ) THEN
308 #ifdef ALLOW_ADAMSBASHFORTH_3
309 CALL ADAMS_BASHFORTH3(
310 I bi, bj, k, Nr,
311 U gT, gtNm, gt_AB,
312 I tempStartAB, iterNb, myThid )
313 #else
314 CALL ADAMS_BASHFORTH2(
315 I bi, bj, k, Nr,
316 U gT, gtNm1, gt_AB,
317 I tempStartAB, iterNb, myThid )
318 #endif
319 #ifdef ALLOW_DIAGNOSTICS
320 IF ( useDiagnostics ) THEN
321 CALL DIAGNOSTICS_FILL(gt_AB,'AB_gT ',k,1,2,bi,bj,myThid)
322 ENDIF
323 #endif /* ALLOW_DIAGNOSTICS */
324 ENDIF
325
326 C-- External thermal forcing term(s) outside Adams-Bashforth:
327 IF ( tempForcing .AND. tracForcingOutAB.EQ.1 )
328 & CALL EXTERNAL_FORCING_T(
329 I iMin, iMax, jMin, jMax, bi, bj, k,
330 I myTime, myThid )
331
332 #ifdef NONLIN_FRSURF
333 IF (nonlinFreeSurf.GT.0) THEN
334 CALL FREESURF_RESCALE_G(
335 I bi, bj, k,
336 U gT,
337 I myThid )
338 IF ( AdamsBashforthGt ) THEN
339 #ifdef ALLOW_ADAMSBASHFORTH_3
340 # ifdef ALLOW_AUTODIFF_TAMC
341 CADJ STORE gtNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey,
342 CADJ & byte=isbyte, kind = isbyte
343 CADJ STORE gtNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey,
344 CADJ & byte=isbyte, kind = isbyte
345 # endif
346 CALL FREESURF_RESCALE_G(
347 I bi, bj, k,
348 U gtNm(1-OLx,1-OLy,1,1,1,1),
349 I myThid )
350 CALL FREESURF_RESCALE_G(
351 I bi, bj, k,
352 U gtNm(1-OLx,1-OLy,1,1,1,2),
353 I myThid )
354 #else
355 CALL FREESURF_RESCALE_G(
356 I bi, bj, k,
357 U gtNm1,
358 I myThid )
359 #endif
360 ENDIF
361 ENDIF
362 #endif /* NONLIN_FRSURF */
363
364 #ifdef ALLOW_ADAMSBASHFORTH_3
365 IF ( AdamsBashforth_T ) THEN
366 CALL TIMESTEP_TRACER(
367 I bi, bj, k, dTtracerLev(k),
368 I gtNm(1-OLx,1-OLy,1,1,1,m2),
369 U gT,
370 I myIter, myThid )
371 ELSE
372 #endif
373 CALL TIMESTEP_TRACER(
374 I bi, bj, k, dTtracerLev(k),
375 I theta,
376 U gT,
377 I myIter, myThid )
378 #ifdef ALLOW_ADAMSBASHFORTH_3
379 ENDIF
380 #endif
381
382 C- end of vertical index (k) loop (Nr:1)
383 ENDDO
384
385 #ifdef ALLOW_DOWN_SLOPE
386 IF ( useDOWN_SLOPE ) THEN
387 IF ( usingPCoords ) THEN
388 CALL DWNSLP_APPLY(
389 I GAD_TEMPERATURE, bi, bj, kSurfC,
390 I recip_drF, recip_hFacC, recip_rA,
391 I dTtracerLev,
392 I theta,
393 U gT,
394 I myTime, myIter, myThid )
395 ELSE
396 CALL DWNSLP_APPLY(
397 I GAD_TEMPERATURE, bi, bj, kLowC,
398 I recip_drF, recip_hFacC, recip_rA,
399 I dTtracerLev,
400 I theta,
401 U gT,
402 I myTime, myIter, myThid )
403 ENDIF
404 ENDIF
405 #endif /* ALLOW_DOWN_SLOPE */
406
407 iMin = 0
408 iMax = sNx+1
409 jMin = 0
410 jMax = sNy+1
411
412 C-- Implicit vertical advection & diffusion
413
414 #ifdef INCLUDE_IMPLVERTADV_CODE
415 IF ( tempImplVertAdv ) THEN
416 #ifdef ALLOW_AUTODIFF_TAMC
417 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
418 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
419 CADJ STORE wFld(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
420 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
421 CADJ STORE recip_hFac(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
422 #endif /* ALLOW_AUTODIFF_TAMC */
423 CALL GAD_IMPLICIT_R(
424 I tempImplVertAdv, tempVertAdvScheme, GAD_TEMPERATURE,
425 I dTtracerLev,
426 I kappaRk, recip_hFac, wFld, theta,
427 U gT,
428 I bi, bj, myTime, myIter, myThid )
429 ELSEIF ( implicitDiffusion ) THEN
430 #else /* INCLUDE_IMPLVERTADV_CODE */
431 IF ( implicitDiffusion ) THEN
432 #endif /* INCLUDE_IMPLVERTADV_CODE */
433 #ifdef ALLOW_AUTODIFF_TAMC
434 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
435 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
436 #endif /* ALLOW_AUTODIFF_TAMC */
437 CALL IMPLDIFF(
438 I bi, bj, iMin, iMax, jMin, jMax,
439 I GAD_TEMPERATURE, kappaRk, recip_hFac,
440 U gT,
441 I myThid )
442 ENDIF
443
444 #ifdef ALLOW_TIMEAVE
445 useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
446 & .OR. useGMredi .OR. ivdc_kappa.NE.0.
447 IF ( taveFreq.GT.0. .AND. useVariableK
448 & .AND.implicitDiffusion ) THEN
449 CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
450 I Nr, 3, deltaTClock, bi, bj, myThid)
451 ENDIF
452 #endif /* ALLOW_TIMEAVE */
453
454 #endif /* ALLOW_GENERIC_ADVDIFF */
455
456 RETURN
457 END

  ViewVC Help
Powered by ViewVC 1.1.22