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

Annotation 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 - (hide annotations) (download)
Fri May 2 05:35:22 2014 UTC (11 years, 2 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 atn 1.3 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 atn 1.1 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 atn 1.2 #ifdef ALLOW_SALT_PLUME
13     # include "SALT_PLUME_OPTIONS.h"
14     #endif
15 atn 1.1
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 atn 1.3 I theta,1,
396 atn 1.1 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