/[MITgcm]/MITgcm/model/src/temp_integrate.F
ViewVC logotype

Annotation of /MITgcm/model/src/temp_integrate.F

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


Revision 1.6 - (hide annotations) (download)
Mon Jan 13 20:49:54 2014 UTC (10 years, 4 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint64u, checkpoint64t
Changes since 1.5: +8 -3 lines
- re-include store directives that were removed on Dec 27th.
- note : this only changes the location of the calc_3d_diffusivity
  recomputation call, which does not change results unless packages
  used in forward (gmredi) are omitted in adjoint (common practice).

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

  ViewVC Help
Powered by ViewVC 1.1.22