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

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

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


Revision 1.19 - (hide annotations) (download)
Wed Mar 6 01:59:02 2002 UTC (22 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint44h_pre
Changes since 1.18: +11 -1 lines
add GM-bolus transport to Eulerian transport to advect Tracers (T,S, ...)

1 jmc 1.19 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.18 2002/03/05 14:15:34 adcroft Exp $
2 adcroft 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6 cnh 1.9 CBOP
7     C !ROUTINE: THERMODYNAMICS
8     C !INTERFACE:
9 adcroft 1.1 SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
10 cnh 1.9 C !DESCRIPTION: \bv
11     C *==========================================================*
12     C | SUBROUTINE THERMODYNAMICS
13     C | o Controlling routine for the prognostic part of the
14     C | thermo-dynamics.
15     C *===========================================================
16     C | The algorithm...
17     C |
18     C | "Correction Step"
19     C | =================
20     C | Here we update the horizontal velocities with the surface
21     C | pressure such that the resulting flow is either consistent
22     C | with the free-surface evolution or the rigid-lid:
23     C | U[n] = U* + dt x d/dx P
24     C | V[n] = V* + dt x d/dy P
25     C |
26     C | "Calculation of Gs"
27     C | ===================
28     C | This is where all the accelerations and tendencies (ie.
29     C | physics, parameterizations etc...) are calculated
30     C | rho = rho ( theta[n], salt[n] )
31     C | b = b(rho, theta)
32     C | K31 = K31 ( rho )
33     C | Gu[n] = Gu( u[n], v[n], wVel, b, ... )
34     C | Gv[n] = Gv( u[n], v[n], wVel, b, ... )
35     C | Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
36     C | Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
37     C |
38     C | "Time-stepping" or "Prediction"
39     C | ================================
40     C | The models variables are stepped forward with the appropriate
41     C | time-stepping scheme (currently we use Adams-Bashforth II)
42     C | - For momentum, the result is always *only* a "prediction"
43     C | in that the flow may be divergent and will be "corrected"
44     C | later with a surface pressure gradient.
45     C | - Normally for tracers the result is the new field at time
46     C | level [n+1} *BUT* in the case of implicit diffusion the result
47     C | is also *only* a prediction.
48     C | - We denote "predictors" with an asterisk (*).
49     C | U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
50     C | V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
51     C | theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
52     C | salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
53     C | With implicit diffusion:
54     C | theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
55     C | salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
56     C | (1 + dt * K * d_zz) theta[n] = theta*
57     C | (1 + dt * K * d_zz) salt[n] = salt*
58     C |
59     C *==========================================================*
60     C \ev
61    
62     C !USES:
63 adcroft 1.1 IMPLICIT NONE
64     C == Global variables ===
65     #include "SIZE.h"
66     #include "EEPARAMS.h"
67     #include "PARAMS.h"
68     #include "DYNVARS.h"
69     #include "GRID.h"
70 adcroft 1.4 #include "GAD.h"
71 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
72     #include "TR1.h"
73     #endif
74     #ifdef ALLOW_AUTODIFF_TAMC
75     # include "tamc.h"
76     # include "tamc_keys.h"
77     # include "FFIELDS.h"
78     # ifdef ALLOW_KPP
79     # include "KPP.h"
80     # endif
81     # ifdef ALLOW_GMREDI
82     # include "GMREDI.h"
83     # endif
84     #endif /* ALLOW_AUTODIFF_TAMC */
85     #ifdef ALLOW_TIMEAVE
86     #include "TIMEAVE_STATV.h"
87     #endif
88    
89 cnh 1.9 C !INPUT/OUTPUT PARAMETERS:
90 adcroft 1.1 C == Routine arguments ==
91     C myTime - Current time in simulation
92     C myIter - Current iteration number in simulation
93     C myThid - Thread number for this instance of the routine.
94     _RL myTime
95     INTEGER myIter
96     INTEGER myThid
97    
98 cnh 1.9 C !LOCAL VARIABLES:
99 adcroft 1.1 C == Local variables
100     C xA, yA - Per block temporaries holding face areas
101     C uTrans, vTrans, rTrans - Per block temporaries holding flow
102     C transport
103     C o uTrans: Zonal transport
104     C o vTrans: Meridional transport
105     C o rTrans: Vertical transport
106     C maskUp o maskUp: land/water mask for W points
107     C fVer[STUV] o fVer: Vertical flux term - note fVer
108     C is "pipelined" in the vertical
109     C so we need an fVer for each
110     C variable.
111     C rhoK, rhoKM1 - Density at current level, and level above
112     C phiHyd - Hydrostatic part of the potential phiHydi.
113     C In z coords phiHydiHyd is the hydrostatic
114     C Potential (=pressure/rho0) anomaly
115     C In p coords phiHydiHyd is the geopotential
116     C surface height anomaly.
117     C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
118     C phiSurfY or geopotentiel (atmos) in X and Y direction
119     C KappaRT, - Total diffusion in vertical for T and S.
120     C KappaRS (background + spatially varying, isopycnal term).
121     C iMin, iMax - Ranges and sub-block indices on which calculations
122     C jMin, jMax are applied.
123     C bi, bj
124     C k, kup, - Index for layer above and below. kup and kDown
125     C kDown, km1 are switched with layer to be the appropriate
126     C index into fVerTerm.
127     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132     _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
134     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
135     _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
136     _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
137     _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139     _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140     _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141     _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
142     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
143     _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
144     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
145     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
146 cnh 1.9 C This is currently used by IVDC and Diagnostics
147 adcroft 1.1 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
148     INTEGER iMin, iMax
149     INTEGER jMin, jMax
150     INTEGER bi, bj
151     INTEGER i, j
152     INTEGER k, km1, kup, kDown
153    
154 cnh 1.9 CEOP
155 adcroft 1.1
156     #ifdef ALLOW_AUTODIFF_TAMC
157     C-- dummy statement to end declaration part
158     ikey = 1
159     #endif /* ALLOW_AUTODIFF_TAMC */
160    
161     C-- Set up work arrays with valid (i.e. not NaN) values
162     C These inital values do not alter the numerical results. They
163     C just ensure that all memory references are to valid floating
164     C point numbers. This prevents spurious hardware signals due to
165     C uninitialised but inert locations.
166     DO j=1-OLy,sNy+OLy
167     DO i=1-OLx,sNx+OLx
168     xA(i,j) = 0. _d 0
169     yA(i,j) = 0. _d 0
170     uTrans(i,j) = 0. _d 0
171     vTrans(i,j) = 0. _d 0
172     DO k=1,Nr
173     phiHyd(i,j,k) = 0. _d 0
174     sigmaX(i,j,k) = 0. _d 0
175     sigmaY(i,j,k) = 0. _d 0
176     sigmaR(i,j,k) = 0. _d 0
177     ENDDO
178     rhoKM1 (i,j) = 0. _d 0
179     rhok (i,j) = 0. _d 0
180     phiSurfX(i,j) = 0. _d 0
181     phiSurfY(i,j) = 0. _d 0
182     ENDDO
183     ENDDO
184    
185    
186     #ifdef ALLOW_AUTODIFF_TAMC
187     C-- HPF directive to help TAMC
188     CHPF$ INDEPENDENT
189     #endif /* ALLOW_AUTODIFF_TAMC */
190    
191     DO bj=myByLo(myThid),myByHi(myThid)
192    
193     #ifdef ALLOW_AUTODIFF_TAMC
194     C-- HPF directive to help TAMC
195 heimbach 1.2 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
196 adcroft 1.1 CHPF$& ,phiHyd,utrans,vtrans,xA,yA
197 heimbach 1.2 CHPF$& ,KappaRT,KappaRS
198 adcroft 1.1 CHPF$& )
199     #endif /* ALLOW_AUTODIFF_TAMC */
200    
201     DO bi=myBxLo(myThid),myBxHi(myThid)
202    
203     #ifdef ALLOW_AUTODIFF_TAMC
204     act1 = bi - myBxLo(myThid)
205     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
206     act2 = bj - myByLo(myThid)
207     max2 = myByHi(myThid) - myByLo(myThid) + 1
208     act3 = myThid - 1
209     max3 = nTx*nTy
210     act4 = ikey_dynamics - 1
211     ikey = (act1 + 1) + act2*max1
212     & + act3*max1*max2
213     & + act4*max1*max2*max3
214     #endif /* ALLOW_AUTODIFF_TAMC */
215    
216     C-- Set up work arrays that need valid initial values
217     DO j=1-OLy,sNy+OLy
218     DO i=1-OLx,sNx+OLx
219     rTrans (i,j) = 0. _d 0
220     fVerT (i,j,1) = 0. _d 0
221     fVerT (i,j,2) = 0. _d 0
222     fVerS (i,j,1) = 0. _d 0
223     fVerS (i,j,2) = 0. _d 0
224     fVerTr1(i,j,1) = 0. _d 0
225     fVerTr1(i,j,2) = 0. _d 0
226     ENDDO
227     ENDDO
228    
229     DO k=1,Nr
230     DO j=1-OLy,sNy+OLy
231     DO i=1-OLx,sNx+OLx
232     C This is currently also used by IVDC and Diagnostics
233     ConvectCount(i,j,k) = 0.
234     KappaRT(i,j,k) = 0. _d 0
235     KappaRS(i,j,k) = 0. _d 0
236 heimbach 1.5 #ifdef ALLOW_AUTODIFF_TAMC
237     gT(i,j,k,bi,bj) = 0. _d 0
238     gS(i,j,k,bi,bj) = 0. _d 0
239     #ifdef ALLOW_PASSIVE_TRACER
240     gTr1(i,j,k,bi,bj) = 0. _d 0
241     #endif
242     #endif
243 adcroft 1.1 ENDDO
244     ENDDO
245     ENDDO
246    
247     iMin = 1-OLx+1
248     iMax = sNx+OLx
249     jMin = 1-OLy+1
250     jMax = sNy+OLy
251    
252    
253     #ifdef ALLOW_AUTODIFF_TAMC
254 heimbach 1.11 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
255     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
256     #ifdef ALLOW_KPP
257     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
258     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
259     #endif
260 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
261    
262     C-- Start of diagnostic loop
263     DO k=Nr,1,-1
264    
265     #ifdef ALLOW_AUTODIFF_TAMC
266     C? Patrick, is this formula correct now that we change the loop range?
267     C? Do we still need this?
268     cph kkey formula corrected.
269     cph Needed for rhok, rhokm1, in the case useGMREDI.
270     kkey = (ikey-1)*Nr + k
271     CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
272     CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
273     #endif /* ALLOW_AUTODIFF_TAMC */
274    
275     C-- Integrate continuity vertically for vertical velocity
276     CALL INTEGRATE_FOR_W(
277     I bi, bj, k, uVel, vVel,
278     O wVel,
279     I myThid )
280    
281     #ifdef ALLOW_OBCS
282     #ifdef ALLOW_NONHYDROSTATIC
283     C-- Apply OBC to W if in N-H mode
284     IF (useOBCS.AND.nonHydrostatic) THEN
285     CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
286     ENDIF
287     #endif /* ALLOW_NONHYDROSTATIC */
288     #endif /* ALLOW_OBCS */
289    
290     C-- Calculate gradients of potential density for isoneutral
291     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
292     c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
293     IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
294     #ifdef ALLOW_AUTODIFF_TAMC
295     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
296     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
297     #endif /* ALLOW_AUTODIFF_TAMC */
298     CALL FIND_RHO(
299     I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
300     I theta, salt,
301     O rhoK,
302     I myThid )
303     IF (k.GT.1) THEN
304     #ifdef ALLOW_AUTODIFF_TAMC
305     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
306     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
307     #endif /* ALLOW_AUTODIFF_TAMC */
308     CALL FIND_RHO(
309     I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
310     I theta, salt,
311     O rhoKm1,
312     I myThid )
313     ENDIF
314     CALL GRAD_SIGMA(
315     I bi, bj, iMin, iMax, jMin, jMax, k,
316     I rhoK, rhoKm1, rhoK,
317     O sigmaX, sigmaY, sigmaR,
318     I myThid )
319     ENDIF
320    
321     C-- Implicit Vertical Diffusion for Convection
322     c ==> should use sigmaR !!!
323     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
324     CALL CALC_IVDC(
325     I bi, bj, iMin, iMax, jMin, jMax, k,
326     I rhoKm1, rhoK,
327     U ConvectCount, KappaRT, KappaRS,
328     I myTime, myIter, myThid)
329     ENDIF
330    
331     C-- end of diagnostic k loop (Nr:1)
332     ENDDO
333    
334     #ifdef ALLOW_AUTODIFF_TAMC
335     cph avoids recomputation of integrate_for_w
336 heimbach 1.11 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
337 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
338    
339     #ifdef ALLOW_OBCS
340     C-- Calculate future values on open boundaries
341     IF (useOBCS) THEN
342 jmc 1.16 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
343 adcroft 1.1 I uVel, vVel, wVel, theta, salt,
344     I myThid )
345     ENDIF
346     #endif /* ALLOW_OBCS */
347    
348     C-- Determines forcing terms based on external fields
349     C relaxation terms, etc.
350     CALL EXTERNAL_FORCING_SURF(
351     I bi, bj, iMin, iMax, jMin, jMax,
352     I myThid )
353     #ifdef ALLOW_AUTODIFF_TAMC
354     cph needed for KPP
355     CADJ STORE surfacetendencyU(:,:,bi,bj)
356     CADJ & = comlev1_bibj, key=ikey, byte=isbyte
357     CADJ STORE surfacetendencyV(:,:,bi,bj)
358     CADJ & = comlev1_bibj, key=ikey, byte=isbyte
359     CADJ STORE surfacetendencyS(:,:,bi,bj)
360     CADJ & = comlev1_bibj, key=ikey, byte=isbyte
361     CADJ STORE surfacetendencyT(:,:,bi,bj)
362     CADJ & = comlev1_bibj, key=ikey, byte=isbyte
363     #endif /* ALLOW_AUTODIFF_TAMC */
364    
365     #ifdef ALLOW_GMREDI
366    
367     #ifdef ALLOW_AUTODIFF_TAMC
368     CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
369     CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
370     CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
371     #endif /* ALLOW_AUTODIFF_TAMC */
372     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
373     IF (useGMRedi) THEN
374 jmc 1.15 CALL GMREDI_CALC_TENSOR(
375     I bi, bj, iMin, iMax, jMin, jMax,
376 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
377     I myThid )
378     #ifdef ALLOW_AUTODIFF_TAMC
379     ELSE
380 jmc 1.15 CALL GMREDI_CALC_TENSOR_DUMMY(
381     I bi, bj, iMin, iMax, jMin, jMax,
382 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
383     I myThid )
384     #endif /* ALLOW_AUTODIFF_TAMC */
385     ENDIF
386    
387     #ifdef ALLOW_AUTODIFF_TAMC
388     CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
389     CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
390     CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
391     #endif /* ALLOW_AUTODIFF_TAMC */
392    
393     #endif /* ALLOW_GMREDI */
394    
395     #ifdef ALLOW_KPP
396     C-- Compute KPP mixing coefficients
397     IF (useKPP) THEN
398     CALL KPP_CALC(
399     I bi, bj, myTime, myThid )
400     #ifdef ALLOW_AUTODIFF_TAMC
401     ELSE
402     CALL KPP_CALC_DUMMY(
403     I bi, bj, myTime, myThid )
404     #endif /* ALLOW_AUTODIFF_TAMC */
405     ENDIF
406    
407     #ifdef ALLOW_AUTODIFF_TAMC
408     CADJ STORE KPPghat (:,:,:,bi,bj)
409     CADJ & , KPPdiffKzT(:,:,:,bi,bj)
410     CADJ & , KPPdiffKzS(:,:,:,bi,bj)
411     CADJ & , KPPfrac (:,: ,bi,bj)
412     CADJ & = comlev1_bibj, key=ikey, byte=isbyte
413     #endif /* ALLOW_AUTODIFF_TAMC */
414    
415     #endif /* ALLOW_KPP */
416    
417     #ifdef ALLOW_AUTODIFF_TAMC
418 heimbach 1.11 CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
419     CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
420     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
421     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
422     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
423     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
424 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
425 heimbach 1.11 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
426 adcroft 1.1 #endif
427     #endif /* ALLOW_AUTODIFF_TAMC */
428    
429     #ifdef ALLOW_AIM
430     C AIM - atmospheric intermediate model, physics package code.
431     C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
432     IF ( useAIM ) THEN
433     CALL TIMER_START('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
434     CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )
435     CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
436     ENDIF
437     #endif /* ALLOW_AIM */
438 jmc 1.14
439     #ifdef ALLOW_TIMEAVE
440     IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
441     CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
442     I deltaTclock, bi, bj, myThid)
443     ENDIF
444     #endif /* ALLOW_TIMEAVE */
445 adcroft 1.4
446 adcroft 1.12 #ifndef DISABLE_MULTIDIM_ADVECTION
447 adcroft 1.4 C-- Some advection schemes are better calculated using a multi-dimensional
448     C method in the absence of any other terms and, if used, is done here.
449 adcroft 1.13 C
450     C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
451     C The default is to use multi-dimensinal advection for non-linear advection
452     C schemes. However, for the sake of efficiency of the adjoint it is necessary
453     C to be able to exclude this scheme to avoid excessive storage and
454     C recomputation. It *is* differentiable, if you need it.
455     C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
456     C disable this section of code.
457 adcroft 1.6 IF (multiDimAdvection) THEN
458     IF (tempStepping .AND.
459     & tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.
460     & tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.
461 heimbach 1.11 & tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN
462     CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
463     U theta,gT,
464 adcroft 1.6 I myTime,myIter,myThid)
465 heimbach 1.11 ENDIF
466 adcroft 1.6 IF (saltStepping .AND.
467     & saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.
468     & saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.
469 heimbach 1.11 & saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN
470     CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
471     U salt,gS,
472 adcroft 1.6 I myTime,myIter,myThid)
473 heimbach 1.11 ENDIF
474 adcroft 1.6 ENDIF
475 adcroft 1.17 C Since passive tracers are configurable separately from T,S we
476     C call the multi-dimensional method for PTRACERS regardless
477     C of whether multiDimAdvection is set or not.
478     #ifdef ALLOW_PTRACERS
479     IF ( usePTRACERS ) THEN
480     CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
481     ENDIF
482     #endif /* ALLOW_PTRACERS */
483 adcroft 1.12 #endif /* DISABLE_MULTIDIM_ADVECTION */
484 adcroft 1.1
485     C-- Start of thermodynamics loop
486     DO k=Nr,1,-1
487     #ifdef ALLOW_AUTODIFF_TAMC
488     C? Patrick Is this formula correct?
489     cph Yes, but I rewrote it.
490     cph Also, the KappaR? need the index and subscript k!
491     kkey = (ikey-1)*Nr + k
492     #endif /* ALLOW_AUTODIFF_TAMC */
493    
494     C-- km1 Points to level above k (=k-1)
495     C-- kup Cycles through 1,2 to point to layer above
496     C-- kDown Cycles through 2,1 to point to current layer
497    
498     km1 = MAX(1,k-1)
499     kup = 1+MOD(k+1,2)
500     kDown= 1+MOD(k,2)
501    
502     iMin = 1-OLx
503     iMax = sNx+OLx
504     jMin = 1-OLy
505     jMax = sNy+OLy
506    
507     C-- Get temporary terms used by tendency routines
508     CALL CALC_COMMON_FACTORS (
509     I bi,bj,iMin,iMax,jMin,jMax,k,
510     O xA,yA,uTrans,vTrans,rTrans,maskUp,
511     I myThid)
512 jmc 1.19
513     #ifdef ALLOW_GMREDI
514     C-- Residual transp = Bolus transp + Eulerian transp
515     IF (useGMRedi) THEN
516     CALL GMREDI_CALC_UVFLOW(
517     & uTrans, vTrans, bi, bj, k, myThid)
518     IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
519     & rTrans, bi, bj, k, myThid)
520     ENDIF
521     #endif /* ALLOW_GMREDI */
522 adcroft 1.1
523     #ifdef ALLOW_AUTODIFF_TAMC
524     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
525     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
526     #endif /* ALLOW_AUTODIFF_TAMC */
527    
528     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
529     C-- Calculate the total vertical diffusivity
530     CALL CALC_DIFFUSIVITY(
531     I bi,bj,iMin,iMax,jMin,jMax,k,
532     I maskUp,
533 heimbach 1.2 O KappaRT,KappaRS,
534 adcroft 1.1 I myThid)
535     #endif
536    
537     iMin = 1-OLx+2
538     iMax = sNx+OLx-1
539     jMin = 1-OLy+2
540     jMax = sNy+OLy-1
541    
542     C-- Calculate active tracer tendencies (gT,gS,...)
543     C and step forward storing result in gTnm1, gSnm1, etc.
544     IF ( tempStepping ) THEN
545     CALL CALC_GT(
546     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
547     I xA,yA,uTrans,vTrans,rTrans,maskUp,
548     I KappaRT,
549     U fVerT,
550 adcroft 1.7 I myTime,myIter,myThid)
551 adcroft 1.1 CALL TIMESTEP_TRACER(
552 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
553 adcroft 1.1 I theta, gT,
554     I myIter, myThid)
555     ENDIF
556     IF ( saltStepping ) THEN
557     CALL CALC_GS(
558     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
559     I xA,yA,uTrans,vTrans,rTrans,maskUp,
560     I KappaRS,
561     U fVerS,
562 adcroft 1.7 I myTime,myIter,myThid)
563 adcroft 1.1 CALL TIMESTEP_TRACER(
564 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
565 adcroft 1.1 I salt, gS,
566     I myIter, myThid)
567     ENDIF
568     #ifdef ALLOW_PASSIVE_TRACER
569     IF ( tr1Stepping ) THEN
570     CALL CALC_GTR1(
571     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
572     I xA,yA,uTrans,vTrans,rTrans,maskUp,
573     I KappaRT,
574     U fVerTr1,
575 heimbach 1.8 I myTime,myIter,myThid)
576 adcroft 1.1 CALL TIMESTEP_TRACER(
577 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
578 adcroft 1.1 I Tr1, gTr1,
579 heimbach 1.8 I myIter,myThid)
580 adcroft 1.1 ENDIF
581     #endif
582 adcroft 1.17 #ifdef ALLOW_PTRACERS
583     IF ( usePTRACERS ) THEN
584     CALL PTRACERS_INTEGERATE(
585     I bi,bj,k,
586     I xA,yA,uTrans,vTrans,rTrans,maskUp,
587     X KappaRS,
588     I myIter,myTime,myThid)
589     ENDIF
590     #endif /* ALLOW_PTRACERS */
591 adcroft 1.1
592     #ifdef ALLOW_OBCS
593     C-- Apply open boundary conditions
594     IF (useOBCS) THEN
595 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
596 adcroft 1.1 END IF
597     #endif /* ALLOW_OBCS */
598    
599     C-- Freeze water
600     IF (allowFreezing) THEN
601     #ifdef ALLOW_AUTODIFF_TAMC
602 heimbach 1.11 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
603 adcroft 1.1 CADJ & , key = kkey, byte = isbyte
604     #endif /* ALLOW_AUTODIFF_TAMC */
605     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
606     END IF
607    
608     C-- end of thermodynamic k loop (Nr:1)
609     ENDDO
610    
611    
612     #ifdef ALLOW_AUTODIFF_TAMC
613     C? Patrick? What about this one?
614 adcroft 1.10 cph Keys iikey and idkey dont seem to be needed
615 adcroft 1.1 cph since storing occurs on different tape for each
616     cph impldiff call anyways.
617 adcroft 1.10 cph Thus, common block comlev1_impl isnt needed either.
618 adcroft 1.1 cph Storing below needed in the case useGMREDI.
619     iikey = (ikey-1)*maximpl
620     #endif /* ALLOW_AUTODIFF_TAMC */
621    
622     C-- Implicit diffusion
623     IF (implicitDiffusion) THEN
624    
625     IF (tempStepping) THEN
626     #ifdef ALLOW_AUTODIFF_TAMC
627     idkey = iikey + 1
628 heimbach 1.11 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
629 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
630     CALL IMPLDIFF(
631     I bi, bj, iMin, iMax, jMin, jMax,
632     I deltaTtracer, KappaRT, recip_HFacC,
633 adcroft 1.7 U gT,
634 adcroft 1.1 I myThid )
635     ENDIF
636    
637     IF (saltStepping) THEN
638     #ifdef ALLOW_AUTODIFF_TAMC
639     idkey = iikey + 2
640 heimbach 1.11 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
641 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
642     CALL IMPLDIFF(
643     I bi, bj, iMin, iMax, jMin, jMax,
644     I deltaTtracer, KappaRS, recip_HFacC,
645 adcroft 1.7 U gS,
646 adcroft 1.1 I myThid )
647     ENDIF
648    
649     #ifdef ALLOW_PASSIVE_TRACER
650     IF (tr1Stepping) THEN
651     #ifdef ALLOW_AUTODIFF_TAMC
652 heimbach 1.11 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
653 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
654     CALL IMPLDIFF(
655     I bi, bj, iMin, iMax, jMin, jMax,
656     I deltaTtracer, KappaRT, recip_HFacC,
657 adcroft 1.7 U gTr1,
658 adcroft 1.1 I myThid )
659     ENDIF
660     #endif
661    
662 adcroft 1.17 #ifdef ALLOW_PTRACERS
663     C Vertical diffusion (implicit) for passive tracers
664     IF ( usePTRACERS ) THEN
665     CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
666     ENDIF
667     #endif /* ALLOW_PTRACERS */
668    
669 adcroft 1.1 #ifdef ALLOW_OBCS
670     C-- Apply open boundary conditions
671     IF (useOBCS) THEN
672     DO K=1,Nr
673 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
674 adcroft 1.1 ENDDO
675     END IF
676     #endif /* ALLOW_OBCS */
677    
678     C-- End If implicitDiffusion
679     ENDIF
680    
681     Ccs-
682     ENDDO
683     ENDDO
684    
685     #ifdef ALLOW_AIM
686     IF ( useAIM ) THEN
687     CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
688     ENDIF
689 adcroft 1.7 _EXCH_XYZ_R8(gT,myThid)
690     _EXCH_XYZ_R8(gS,myThid)
691 adcroft 1.1 #else
692     IF (staggerTimeStep.AND.useCubedSphereExchange) THEN
693 adcroft 1.7 _EXCH_XYZ_R8(gT,myThid)
694     _EXCH_XYZ_R8(gS,myThid)
695 adcroft 1.1 ENDIF
696     #endif /* ALLOW_AIM */
697 adcroft 1.17
698     #ifndef DISABLE_DEBUGMODE
699     If (debugMode) THEN
700     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
701     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
702     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
703     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
704     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
705     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
706     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
707     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
708     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
709 adcroft 1.18 #ifdef ALLOW_PTRACERS
710     IF ( usePTRACERS ) THEN
711     CALL PTRACERS_DEBUG(myThid)
712     ENDIF
713     #endif /* ALLOW_PTRACERS */
714 adcroft 1.17 ENDIF
715     #endif
716 adcroft 1.1
717     RETURN
718     END

  ViewVC Help
Powered by ViewVC 1.1.22