/[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.1 - (hide annotations) (download)
Tue Apr 22 04:32:26 2014 UTC (11 years, 2 months ago) by atn
Branch: MAIN
in progress:
1. sort out bi,bj loop
2. skip remove saltplumeflux in external_forcing_surf, (thus) skip kpp
3. move saltplume outside k-loop in [salt,temp]_integrate.F
4. add diagnostics

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

  ViewVC Help
Powered by ViewVC 1.1.22