/[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.8 - (show annotations) (download)
Fri Jul 11 18:43:41 2014 UTC (9 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.7: +17 -14 lines
- new file "apply_forcing.F" containing all the code previously in
  external_forcing.F, but with new argument list: pass, as new argument,
  the current level tendency array to update (instead of a direct update
  of the common bloc array). Change the corresponding calls.

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

  ViewVC Help
Powered by ViewVC 1.1.22