/[MITgcm]/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/temp_integrate.F
ViewVC logotype

Contents of /MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/temp_integrate.F

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


Revision 1.3 - (show annotations) (download)
Fri May 2 05:35:22 2014 UTC (11 years, 3 months ago) by atn
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +2 -2 lines
bug fix. move #include back inside autodiff bracket.

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

  ViewVC Help
Powered by ViewVC 1.1.22