/[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.18 - (hide annotations) (download)
Tue Mar 5 14:15:34 2002 UTC (22 years, 3 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint44f_post, checkpoint44g_post
Changes since 1.17: +6 -2 lines
Bug fix: Missing #ifdef's for when PTRACERS are not enabled.

1 adcroft 1.18 C $Header: /u/gcmpack/models/MITgcmUV/model/src/thermodynamics.F,v 1.17 2002/03/04 17:26:41 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    
513     #ifdef ALLOW_AUTODIFF_TAMC
514     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
515     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
516     #endif /* ALLOW_AUTODIFF_TAMC */
517    
518     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
519     C-- Calculate the total vertical diffusivity
520     CALL CALC_DIFFUSIVITY(
521     I bi,bj,iMin,iMax,jMin,jMax,k,
522     I maskUp,
523 heimbach 1.2 O KappaRT,KappaRS,
524 adcroft 1.1 I myThid)
525     #endif
526    
527     iMin = 1-OLx+2
528     iMax = sNx+OLx-1
529     jMin = 1-OLy+2
530     jMax = sNy+OLy-1
531    
532     C-- Calculate active tracer tendencies (gT,gS,...)
533     C and step forward storing result in gTnm1, gSnm1, etc.
534     IF ( tempStepping ) THEN
535     CALL CALC_GT(
536     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
537     I xA,yA,uTrans,vTrans,rTrans,maskUp,
538     I KappaRT,
539     U fVerT,
540 adcroft 1.7 I myTime,myIter,myThid)
541 adcroft 1.1 CALL TIMESTEP_TRACER(
542 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
543 adcroft 1.1 I theta, gT,
544     I myIter, myThid)
545     ENDIF
546     IF ( saltStepping ) THEN
547     CALL CALC_GS(
548     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
549     I xA,yA,uTrans,vTrans,rTrans,maskUp,
550     I KappaRS,
551     U fVerS,
552 adcroft 1.7 I myTime,myIter,myThid)
553 adcroft 1.1 CALL TIMESTEP_TRACER(
554 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
555 adcroft 1.1 I salt, gS,
556     I myIter, myThid)
557     ENDIF
558     #ifdef ALLOW_PASSIVE_TRACER
559     IF ( tr1Stepping ) THEN
560     CALL CALC_GTR1(
561     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
562     I xA,yA,uTrans,vTrans,rTrans,maskUp,
563     I KappaRT,
564     U fVerTr1,
565 heimbach 1.8 I myTime,myIter,myThid)
566 adcroft 1.1 CALL TIMESTEP_TRACER(
567 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
568 adcroft 1.1 I Tr1, gTr1,
569 heimbach 1.8 I myIter,myThid)
570 adcroft 1.1 ENDIF
571     #endif
572 adcroft 1.17 #ifdef ALLOW_PTRACERS
573     IF ( usePTRACERS ) THEN
574     CALL PTRACERS_INTEGERATE(
575     I bi,bj,k,
576     I xA,yA,uTrans,vTrans,rTrans,maskUp,
577     X KappaRS,
578     I myIter,myTime,myThid)
579     ENDIF
580     #endif /* ALLOW_PTRACERS */
581 adcroft 1.1
582     #ifdef ALLOW_OBCS
583     C-- Apply open boundary conditions
584     IF (useOBCS) THEN
585 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
586 adcroft 1.1 END IF
587     #endif /* ALLOW_OBCS */
588    
589     C-- Freeze water
590     IF (allowFreezing) THEN
591     #ifdef ALLOW_AUTODIFF_TAMC
592 heimbach 1.11 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
593 adcroft 1.1 CADJ & , key = kkey, byte = isbyte
594     #endif /* ALLOW_AUTODIFF_TAMC */
595     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
596     END IF
597    
598     C-- end of thermodynamic k loop (Nr:1)
599     ENDDO
600    
601    
602     #ifdef ALLOW_AUTODIFF_TAMC
603     C? Patrick? What about this one?
604 adcroft 1.10 cph Keys iikey and idkey dont seem to be needed
605 adcroft 1.1 cph since storing occurs on different tape for each
606     cph impldiff call anyways.
607 adcroft 1.10 cph Thus, common block comlev1_impl isnt needed either.
608 adcroft 1.1 cph Storing below needed in the case useGMREDI.
609     iikey = (ikey-1)*maximpl
610     #endif /* ALLOW_AUTODIFF_TAMC */
611    
612     C-- Implicit diffusion
613     IF (implicitDiffusion) THEN
614    
615     IF (tempStepping) THEN
616     #ifdef ALLOW_AUTODIFF_TAMC
617     idkey = iikey + 1
618 heimbach 1.11 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
619 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
620     CALL IMPLDIFF(
621     I bi, bj, iMin, iMax, jMin, jMax,
622     I deltaTtracer, KappaRT, recip_HFacC,
623 adcroft 1.7 U gT,
624 adcroft 1.1 I myThid )
625     ENDIF
626    
627     IF (saltStepping) THEN
628     #ifdef ALLOW_AUTODIFF_TAMC
629     idkey = iikey + 2
630 heimbach 1.11 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
631 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
632     CALL IMPLDIFF(
633     I bi, bj, iMin, iMax, jMin, jMax,
634     I deltaTtracer, KappaRS, recip_HFacC,
635 adcroft 1.7 U gS,
636 adcroft 1.1 I myThid )
637     ENDIF
638    
639     #ifdef ALLOW_PASSIVE_TRACER
640     IF (tr1Stepping) THEN
641     #ifdef ALLOW_AUTODIFF_TAMC
642 heimbach 1.11 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
643 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
644     CALL IMPLDIFF(
645     I bi, bj, iMin, iMax, jMin, jMax,
646     I deltaTtracer, KappaRT, recip_HFacC,
647 adcroft 1.7 U gTr1,
648 adcroft 1.1 I myThid )
649     ENDIF
650     #endif
651    
652 adcroft 1.17 #ifdef ALLOW_PTRACERS
653     C Vertical diffusion (implicit) for passive tracers
654     IF ( usePTRACERS ) THEN
655     CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
656     ENDIF
657     #endif /* ALLOW_PTRACERS */
658    
659 adcroft 1.1 #ifdef ALLOW_OBCS
660     C-- Apply open boundary conditions
661     IF (useOBCS) THEN
662     DO K=1,Nr
663 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
664 adcroft 1.1 ENDDO
665     END IF
666     #endif /* ALLOW_OBCS */
667    
668     C-- End If implicitDiffusion
669     ENDIF
670    
671     Ccs-
672     ENDDO
673     ENDDO
674    
675     #ifdef ALLOW_AIM
676     IF ( useAIM ) THEN
677     CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
678     ENDIF
679 adcroft 1.7 _EXCH_XYZ_R8(gT,myThid)
680     _EXCH_XYZ_R8(gS,myThid)
681 adcroft 1.1 #else
682     IF (staggerTimeStep.AND.useCubedSphereExchange) THEN
683 adcroft 1.7 _EXCH_XYZ_R8(gT,myThid)
684     _EXCH_XYZ_R8(gS,myThid)
685 adcroft 1.1 ENDIF
686     #endif /* ALLOW_AIM */
687 adcroft 1.17
688     #ifndef DISABLE_DEBUGMODE
689     If (debugMode) THEN
690     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
691     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
692     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
693     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
694     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
695     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
696     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
697     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
698     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
699 adcroft 1.18 #ifdef ALLOW_PTRACERS
700     IF ( usePTRACERS ) THEN
701     CALL PTRACERS_DEBUG(myThid)
702     ENDIF
703     #endif /* ALLOW_PTRACERS */
704 adcroft 1.17 ENDIF
705     #endif
706 adcroft 1.1
707     RETURN
708     END

  ViewVC Help
Powered by ViewVC 1.1.22