/[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.5 - (show annotations) (download)
Fri Dec 27 15:46:46 2013 UTC (10 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64s
Changes since 1.4: +120 -39 lines
- move calls to GAD_IMPLICIT_R & IMPLDIFF + DWNSLP_APPLY from thermodynamics.F
  to inside temp_integrate.F, salt_integrate.F & ptracers_integrate
- add fZon & fMer as output argument of S/R GAD_CALC_RHS

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

  ViewVC Help
Powered by ViewVC 1.1.22