/[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.8 - (hide 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 jmc 1.8 C $Header: /u/gcmpack/MITgcm/model/src/temp_integrate.F,v 1.7 2014/04/04 20:54:11 jmc Exp $
2 jmc 1.7 C $Name: $
3 jmc 1.1
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6 jmc 1.7 #ifdef ALLOW_AUTODIFF
7     # include "AUTODIFF_OPTIONS.h"
8     #endif
9 jmc 1.4 #ifdef ALLOW_GENERIC_ADVDIFF
10     # include "GAD_OPTIONS.h"
11     #endif
12 jmc 1.1
13     CBOP
14     C !ROUTINE: TEMP_INTEGRATE
15     C !INTERFACE:
16     SUBROUTINE TEMP_INTEGRATE(
17 jmc 1.5 I bi, bj, recip_hFac,
18 jmc 1.4 I uFld, vFld, wFld,
19     U KappaRk,
20 jmc 1.1 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 jmc 1.8 C | A procedure called APPLY_FORCING_T is called from
28 jmc 1.1 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 jmc 1.4 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 jmc 1.1 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 jmc 1.5 #include "GRID.h"
57     #include "DYNVARS.h"
58 jmc 1.1 #include "RESTART.h"
59     #ifdef ALLOW_GENERIC_ADVDIFF
60 jmc 1.4 # include "GAD.h"
61     # include "GAD_SOM_VARS.h"
62 jmc 1.1 #endif
63 jmc 1.5 #ifdef ALLOW_TIMEAVE
64     # include "TIMEAVE_STATV.h"
65     #endif
66 jmc 1.7 #ifdef ALLOW_AUTODIFF
67 jmc 1.1 # include "tamc.h"
68     # include "tamc_keys.h"
69     #endif
70    
71     C !INPUT/OUTPUT PARAMETERS:
72     C == Routine arguments ==
73 jmc 1.5 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 jmc 1.4 INTEGER bi, bj
82 jmc 1.5 _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 jmc 1.1 _RL myTime
88     INTEGER myIter
89     INTEGER myThid
90     CEOP
91    
92     #ifdef ALLOW_GENERIC_ADVDIFF
93     C !LOCAL VARIABLES:
94 jmc 1.5 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 jmc 1.8 C gt_AB :: Adams-Bashforth temperature tendency increment
112 jmc 1.5 C useVariableK :: T when vertical diffusion is not constant
113 jmc 1.4 INTEGER iMin, iMax, jMin, jMax
114 jmc 1.1 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 jmc 1.5 _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 jmc 1.1 _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 jmc 1.5 #ifdef ALLOW_TIMEAVE
133     LOGICAL useVariableK
134     #endif
135    
136 jmc 1.4 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 jmc 1.5 IF (staggerTimeStep) iterNb = myIter - 1
146 jmc 1.1
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 jmc 1.4 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 jmc 1.1 DO j=1-OLy,sNy+OLy
169     DO i=1-OLx,sNx+OLx
170 jmc 1.5 fVer(i,j,1) = 0. _d 0
171     fVer(i,j,2) = 0. _d 0
172 jmc 1.1 ENDDO
173     ENDDO
174 jmc 1.4 #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 jmc 1.5 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
183     CADJ STORE wFld(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
184 jmc 1.4 #endif /* ALLOW_AUTODIFF */
185    
186     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
187 jmc 1.5 CALL CALC_3D_DIFFUSIVITY(
188 jmc 1.4 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 jmc 1.8 I tempImplVertAdv,
215     I tempAdvScheme, tempVertAdvScheme, GAD_TEMPERATURE,
216     I dTtracerLev, uFld, vFld, wFld, theta,
217 jmc 1.4 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 jmc 1.8 I tempImplVertAdv,
229     I tempAdvScheme, tempVertAdvScheme, GAD_TEMPERATURE,
230     I dTtracerLev, uFld, vFld, wFld, theta,
231 jmc 1.4 O gT,
232     I bi, bj, myTime, myIter, myThid )
233     ENDIF
234     #endif /* DISABLE_MULTIDIM_ADVECTION */
235 jmc 1.1
236 jmc 1.4 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 jmc 1.1 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 jmc 1.5 CADJ STORE fVer(:,:,:) = comlev1_bibj_k, key=kkey,
250 jmc 1.1 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 jmc 1.5 O fZon, fMer,
284     U fVer, gT,
285 jmc 1.1 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 jmc 1.5 O fZon, fMer,
298     U fVer, gT,
299 jmc 1.1 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 jmc 1.8 & 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 jmc 1.1
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 jmc 1.8 & 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 jmc 1.1
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 jmc 1.5 #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 jmc 1.7 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 jmc 1.5 #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 jmc 1.7 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
438 jmc 1.5 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 jmc 1.1 #endif /* ALLOW_GENERIC_ADVDIFF */
458    
459     RETURN
460     END

  ViewVC Help
Powered by ViewVC 1.1.22