/[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.12 - (show annotations) (download)
Wed Jul 30 03:31:35 2014 UTC (9 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.11: +3 -3 lines
-  add new pkg "kl10" for mixing due to internal wave breaking
  ( http://www.sciencedirect.com/science/article/pii/S1463500310000144 )

1 C $Header: /u/gcmpack/MITgcm/model/src/temp_integrate.F,v 1.11 2014/07/22 12:04:09 jmc 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 APPLY_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 #ifdef ALLOW_DIAGNOSTICS
94 C !FUNCTIONS:
95 LOGICAL DIAGNOSTICS_IS_ON
96 EXTERNAL DIAGNOSTICS_IS_ON
97 #endif /* ALLOW_DIAGNOSTICS */
98
99 C !LOCAL VARIABLES:
100 C iMin, iMax :: 1rst index loop range
101 C jMin, jMax :: 2nd index loop range
102 C k :: vertical index
103 C kM1 :: =k-1 for k>1, =1 for k=1
104 C kUp :: index into 2 1/2D array, toggles between 1|2
105 C kDown :: index into 2 1/2D array, toggles between 2|1
106 C xA :: Tracer cell face area normal to X
107 C yA :: Tracer cell face area normal to X
108 C maskUp :: Land/water mask for Wvel points (interface k)
109 C uTrans :: Zonal volume transport through cell face
110 C vTrans :: Meridional volume transport through cell face
111 C rTrans :: Vertical volume transport at interface k
112 C rTransKp :: Vertical volume transport at inteface k+1
113 C fZon :: Flux of temperature (T) in the zonal direction
114 C fMer :: Flux of temperature (T) in the meridional direction
115 C fVer :: Flux of temperature (T) in the vertical direction
116 C at the upper(U) and lower(D) faces of a cell.
117 C gtForc :: Temperature forcing tendency
118 C gt_AB :: Adams-Bashforth temperature tendency increment
119 C useVariableK :: T when vertical diffusion is not constant
120 INTEGER iMin, iMax, jMin, jMax
121 INTEGER i, j, k
122 INTEGER kUp, kDown, kM1
123 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 _RL fVer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
133 _RL gtForc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134 _RL gt_AB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135 #ifdef ALLOW_DIAGNOSTICS
136 LOGICAL diagForcing, diagAB_tend
137 #endif
138 LOGICAL calcAdvection
139 INTEGER iterNb
140 #ifdef ALLOW_ADAMSBASHFORTH_3
141 INTEGER m2
142 #endif
143 #ifdef ALLOW_TIMEAVE
144 LOGICAL useVariableK
145 #endif
146
147 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
148
149 C- Loop ranges for daughter routines
150 iMin = 1-OLx+2
151 iMax = sNx+OLx-1
152 jMin = 1-OLy+2
153 jMax = sNy+OLy-1
154
155 iterNb = myIter
156 IF (staggerTimeStep) iterNb = myIter - 1
157
158 #ifdef ALLOW_DIAGNOSTICS
159 diagForcing = .FALSE.
160 diagAB_tend = .FALSE.
161 IF ( useDiagnostics .AND. tempForcing )
162 & diagForcing = DIAGNOSTICS_IS_ON( 'gT_Forc ', myThid )
163 IF ( useDiagnostics .AND. AdamsBashforthGt )
164 & diagAB_tend = DIAGNOSTICS_IS_ON( 'AB_gT ', myThid )
165 #endif
166
167 #ifdef ALLOW_AUTODIFF_TAMC
168 act1 = bi - myBxLo(myThid)
169 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
170 act2 = bj - myByLo(myThid)
171 max2 = myByHi(myThid) - myByLo(myThid) + 1
172 act3 = myThid - 1
173 max3 = nTx*nTy
174 act4 = ikey_dynamics - 1
175 itdkey = (act1 + 1) + act2*max1
176 & + act3*max1*max2
177 & + act4*max1*max2*max3
178 #endif /* ALLOW_AUTODIFF_TAMC */
179
180 C- Tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
181 DO k=1,Nr
182 DO j=1-OLy,sNy+OLy
183 DO i=1-OLx,sNx+OLx
184 gT(i,j,k,bi,bj) = 0. _d 0
185 ENDDO
186 ENDDO
187 ENDDO
188 DO j=1-OLy,sNy+OLy
189 DO i=1-OLx,sNx+OLx
190 fVer(i,j,1) = 0. _d 0
191 fVer(i,j,2) = 0. _d 0
192 ENDDO
193 ENDDO
194 #ifdef ALLOW_AUTODIFF
195 DO k=1,Nr
196 DO j=1-OLy,sNy+OLy
197 DO i=1-OLx,sNx+OLx
198 kappaRk(i,j,k) = 0. _d 0
199 ENDDO
200 ENDDO
201 ENDDO
202 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
203 CADJ STORE wFld(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
204 #endif /* ALLOW_AUTODIFF */
205
206 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
207 CALL CALC_3D_DIFFUSIVITY(
208 I bi, bj, iMin, iMax, jMin, jMax,
209 I GAD_TEMPERATURE, useGMredi, useKPP,
210 O kappaRk,
211 I myThid )
212 #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
213
214 #ifndef DISABLE_MULTIDIM_ADVECTION
215 C-- Some advection schemes are better calculated using a multi-dimensional
216 C method in the absence of any other terms and, if used, is done here.
217 C
218 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
219 C The default is to use multi-dimensinal advection for non-linear advection
220 C schemes. However, for the sake of efficiency of the adjoint it is necessary
221 C to be able to exclude this scheme to avoid excessive storage and
222 C recomputation. It *is* differentiable, if you need it.
223 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
224 C disable this section of code.
225 #ifdef GAD_ALLOW_TS_SOM_ADV
226 # ifdef ALLOW_AUTODIFF_TAMC
227 CADJ STORE som_T = comlev1_bibj, key=itdkey, byte=isbyte
228 # endif
229 IF ( tempSOM_Advection ) THEN
230 # ifdef ALLOW_DEBUG
231 IF (debugMode) CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
232 # endif
233 CALL GAD_SOM_ADVECT(
234 I tempImplVertAdv,
235 I tempAdvScheme, tempVertAdvScheme, GAD_TEMPERATURE,
236 I dTtracerLev, uFld, vFld, wFld, theta,
237 U som_T,
238 O gT,
239 I bi, bj, myTime, myIter, myThid )
240 ELSEIF (tempMultiDimAdvec) THEN
241 #else /* GAD_ALLOW_TS_SOM_ADV */
242 IF (tempMultiDimAdvec) THEN
243 #endif /* GAD_ALLOW_TS_SOM_ADV */
244 # ifdef ALLOW_DEBUG
245 IF (debugMode) CALL DEBUG_CALL('GAD_ADVECTION',myThid)
246 # endif
247 CALL GAD_ADVECTION(
248 I tempImplVertAdv,
249 I tempAdvScheme, tempVertAdvScheme, GAD_TEMPERATURE,
250 I dTtracerLev, uFld, vFld, wFld, theta,
251 O gT,
252 I bi, bj, myTime, myIter, myThid )
253 ENDIF
254 #endif /* DISABLE_MULTIDIM_ADVECTION */
255
256 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
257
258 C- Start vertical index (k) loop (Nr:1)
259 calcAdvection = tempAdvection .AND. .NOT.tempMultiDimAdvec
260 DO k=Nr,1,-1
261 #ifdef ALLOW_AUTODIFF_TAMC
262 kkey = (itdkey-1)*Nr + k
263 #endif /* ALLOW_AUTODIFF_TAMC */
264 kM1 = MAX(1,k-1)
265 kUp = 1+MOD(k+1,2)
266 kDown= 1+MOD(k,2)
267
268 #ifdef ALLOW_AUTODIFF_TAMC
269 CADJ STORE fVer(:,:,:) = comlev1_bibj_k, key=kkey,
270 CADJ & byte=isbyte, kind = isbyte
271 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
272 CADJ & byte=isbyte, kind = isbyte
273 # ifdef ALLOW_ADAMSBASHFORTH_3
274 CADJ STORE gtNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey,
275 CADJ & byte=isbyte, kind = isbyte
276 CADJ STORE gtNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey,
277 CADJ & byte=isbyte, kind = isbyte
278 # else
279 CADJ STORE gtNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
280 CADJ & byte=isbyte, kind = isbyte
281 # endif
282 #endif /* ALLOW_AUTODIFF_TAMC */
283 CALL CALC_ADV_FLOW(
284 I uFld, vFld, wFld,
285 U rTrans,
286 O uTrans, vTrans, rTransKp,
287 O maskUp, xA, yA,
288 I k, bi, bj, myThid )
289
290 C-- Collect forcing term in local array gtForc:
291 DO j=1-OLy,sNy+OLy
292 DO i=1-OLx,sNx+OLx
293 gtForc(i,j) = 0. _d 0
294 ENDDO
295 ENDDO
296 IF ( tempForcing ) THEN
297 CALL APPLY_FORCING_T(
298 U gtForc,
299 I iMin,iMax,jMin,jMax, k, bi,bj,
300 I myTime, myIter, myThid )
301 #ifdef ALLOW_DIAGNOSTICS
302 IF ( diagForcing ) THEN
303 CALL DIAGNOSTICS_FILL(gtForc,'gT_Forc ',k,1,2,bi,bj,myThid)
304 ENDIF
305 #endif /* ALLOW_DIAGNOSTICS */
306 ENDIF
307
308 #ifdef ALLOW_ADAMSBASHFORTH_3
309 c m1 = 1 + MOD(iterNb+1,2)
310 m2 = 1 + MOD( iterNb ,2)
311 CALL GAD_CALC_RHS(
312 I bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
313 I xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
314 I vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
315 I uTrans, vTrans, rTrans, rTransKp,
316 I diffKhT, diffK4T, KappaRk(1-OLx,1-OLy,k), diffKr4T,
317 I theta, gtNm(1-OLx,1-OLy,1,1,1,m2), dTtracerLev,
318 I GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
319 I calcAdvection, tempImplVertAdv, AdamsBashforth_T,
320 I tempVertDiff4, useGMRedi, useKPP,
321 O fZon, fMer,
322 U fVer, gT,
323 I myTime, myIter, myThid )
324 #else /* ALLOW_ADAMSBASHFORTH_3 */
325 CALL GAD_CALC_RHS(
326 I bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
327 I xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
328 I vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
329 I uTrans, vTrans, rTrans, rTransKp,
330 I diffKhT, diffK4T, KappaRk(1-OLx,1-OLy,k), diffKr4T,
331 I theta, gtNm1, dTtracerLev,
332 I GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
333 I calcAdvection, tempImplVertAdv, AdamsBashforth_T,
334 I tempVertDiff4, useGMRedi, useKPP,
335 O fZon, fMer,
336 U fVer, gT,
337 I myTime, myIter, myThid )
338 #endif
339
340 C-- External thermal forcing term(s) inside Adams-Bashforth:
341 IF ( tempForcing .AND. tracForcingOutAB.NE.1 ) THEN
342 DO j=1-OLy,sNy+OLy
343 DO i=1-OLx,sNx+OLx
344 gT(i,j,k,bi,bj) = gT(i,j,k,bi,bj) + gtForc(i,j)
345 ENDDO
346 ENDDO
347 ENDIF
348
349 IF ( AdamsBashforthGt ) THEN
350 #ifdef ALLOW_ADAMSBASHFORTH_3
351 CALL ADAMS_BASHFORTH3(
352 I bi, bj, k, Nr,
353 U gT, gtNm, gt_AB,
354 I tempStartAB, iterNb, myThid )
355 #else
356 CALL ADAMS_BASHFORTH2(
357 I bi, bj, k, Nr,
358 U gT, gtNm1, gt_AB,
359 I tempStartAB, iterNb, myThid )
360 #endif
361 #ifdef ALLOW_DIAGNOSTICS
362 IF ( diagAB_tend ) THEN
363 CALL DIAGNOSTICS_FILL(gt_AB,'AB_gT ',k,1,2,bi,bj,myThid)
364 ENDIF
365 #endif /* ALLOW_DIAGNOSTICS */
366 ENDIF
367
368 C-- External thermal forcing term(s) outside Adams-Bashforth:
369 IF ( tempForcing .AND. tracForcingOutAB.EQ.1 ) THEN
370 DO j=1-OLy,sNy+OLy
371 DO i=1-OLx,sNx+OLx
372 gT(i,j,k,bi,bj) = gT(i,j,k,bi,bj) + gtForc(i,j)
373 ENDDO
374 ENDDO
375 ENDIF
376
377 #ifdef NONLIN_FRSURF
378 IF (nonlinFreeSurf.GT.0) THEN
379 CALL FREESURF_RESCALE_G(
380 I bi, bj, k,
381 U gT,
382 I myThid )
383 IF ( AdamsBashforthGt ) THEN
384 #ifdef ALLOW_ADAMSBASHFORTH_3
385 # ifdef ALLOW_AUTODIFF_TAMC
386 CADJ STORE gtNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey,
387 CADJ & byte=isbyte, kind = isbyte
388 CADJ STORE gtNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey,
389 CADJ & byte=isbyte, kind = isbyte
390 # endif
391 CALL FREESURF_RESCALE_G(
392 I bi, bj, k,
393 U gtNm(1-OLx,1-OLy,1,1,1,1),
394 I myThid )
395 CALL FREESURF_RESCALE_G(
396 I bi, bj, k,
397 U gtNm(1-OLx,1-OLy,1,1,1,2),
398 I myThid )
399 #else
400 CALL FREESURF_RESCALE_G(
401 I bi, bj, k,
402 U gtNm1,
403 I myThid )
404 #endif
405 ENDIF
406 ENDIF
407 #endif /* NONLIN_FRSURF */
408
409 CALL TIMESTEP_TRACER(
410 I bi, bj, k, dTtracerLev(k),
411 I theta,
412 U gT,
413 I myIter, myThid )
414
415 C- end of vertical index (k) loop (Nr:1)
416 ENDDO
417
418 #ifdef ALLOW_DOWN_SLOPE
419 IF ( useDOWN_SLOPE ) THEN
420 IF ( usingPCoords ) THEN
421 CALL DWNSLP_APPLY(
422 I GAD_TEMPERATURE, bi, bj, kSurfC,
423 I recip_drF, recip_hFac, recip_rA,
424 I dTtracerLev,
425 I theta,
426 U gT,
427 I myTime, myIter, myThid )
428 ELSE
429 CALL DWNSLP_APPLY(
430 I GAD_TEMPERATURE, bi, bj, kLowC,
431 I recip_drF, recip_hFac, recip_rA,
432 I dTtracerLev,
433 I theta,
434 U gT,
435 I myTime, myIter, myThid )
436 ENDIF
437 ENDIF
438 #endif /* ALLOW_DOWN_SLOPE */
439
440 iMin = 0
441 iMax = sNx+1
442 jMin = 0
443 jMax = sNy+1
444
445 C-- Implicit vertical advection & diffusion
446
447 #ifdef INCLUDE_IMPLVERTADV_CODE
448 IF ( tempImplVertAdv ) THEN
449 #ifdef ALLOW_AUTODIFF_TAMC
450 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
451 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
452 CADJ STORE wFld(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
453 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
454 CADJ STORE recip_hFac(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
455 #endif /* ALLOW_AUTODIFF_TAMC */
456 CALL GAD_IMPLICIT_R(
457 I tempImplVertAdv, tempVertAdvScheme, GAD_TEMPERATURE,
458 I dTtracerLev,
459 I kappaRk, recip_hFac, wFld, theta,
460 U gT,
461 I bi, bj, myTime, myIter, myThid )
462 ELSEIF ( implicitDiffusion ) THEN
463 #else /* INCLUDE_IMPLVERTADV_CODE */
464 IF ( implicitDiffusion ) THEN
465 #endif /* INCLUDE_IMPLVERTADV_CODE */
466 #ifdef ALLOW_AUTODIFF_TAMC
467 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
468 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
469 #endif /* ALLOW_AUTODIFF_TAMC */
470 CALL IMPLDIFF(
471 I bi, bj, iMin, iMax, jMin, jMax,
472 I GAD_TEMPERATURE, kappaRk, recip_hFac,
473 U gT,
474 I myThid )
475 ENDIF
476
477 #ifdef ALLOW_TIMEAVE
478 useVariableK = useKPP .OR. usePP81 .OR. useKL10 .OR. useMY82
479 & .OR. useGGL90 .OR. useGMredi .OR. ivdc_kappa.NE.0.
480 IF ( taveFreq.GT.0. .AND. useVariableK
481 & .AND.implicitDiffusion ) THEN
482 CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
483 I Nr, 3, deltaTClock, bi, bj, myThid)
484 ENDIF
485 #endif /* ALLOW_TIMEAVE */
486
487 #endif /* ALLOW_GENERIC_ADVDIFF */
488
489 RETURN
490 END

  ViewVC Help
Powered by ViewVC 1.1.22