/[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.70 - (hide annotations) (download)
Thu Jul 1 21:45:11 2004 UTC (19 years, 11 months ago) by jmc
Branch: MAIN
Changes since 1.69: +2 -12 lines 
prepare splitting of thermodynamics: store convective counter in common block
1 jmc 1.70 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.69 2004/06/26 02:38:09 jmc Exp $
2 adcroft 1.1 C $Name: $
3    
4 edhill 1.51 #include "PACKAGES_CONFIG.h"
5 adcroft 1.1 #include "CPP_OPTIONS.h"
6 jmc 1.60 #ifdef ALLOW_PTRACERS
7     # include "PTRACERS_OPTIONS.h"
8     #endif
9 edhill 1.51
10 jmc 1.21 #ifdef ALLOW_AUTODIFF_TAMC
11     # ifdef ALLOW_GMREDI
12     # include "GMREDI_OPTIONS.h"
13     # endif
14     # ifdef ALLOW_KPP
15     # include "KPP_OPTIONS.h"
16 heimbach 1.42 # endif
17 jmc 1.21 #endif /* ALLOW_AUTODIFF_TAMC */
18 adcroft 1.1
19 cnh 1.9 CBOP
20     C !ROUTINE: THERMODYNAMICS
21     C !INTERFACE:
22 adcroft 1.1 SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
23 cnh 1.9 C !DESCRIPTION: \bv
24     C *==========================================================*
25     C | SUBROUTINE THERMODYNAMICS
26     C | o Controlling routine for the prognostic part of the
27     C | thermo-dynamics.
28     C *===========================================================
29     C | The algorithm...
30     C |
31     C | "Correction Step"
32     C | =================
33     C | Here we update the horizontal velocities with the surface
34     C | pressure such that the resulting flow is either consistent
35     C | with the free-surface evolution or the rigid-lid:
36     C | U[n] = U* + dt x d/dx P
37     C | V[n] = V* + dt x d/dy P
38     C |
39     C | "Calculation of Gs"
40     C | ===================
41     C | This is where all the accelerations and tendencies (ie.
42     C | physics, parameterizations etc...) are calculated
43     C | rho = rho ( theta[n], salt[n] )
44     C | b = b(rho, theta)
45     C | K31 = K31 ( rho )
46     C | Gu[n] = Gu( u[n], v[n], wVel, b, ... )
47     C | Gv[n] = Gv( u[n], v[n], wVel, b, ... )
48     C | Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
49     C | Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
50     C |
51     C | "Time-stepping" or "Prediction"
52     C | ================================
53     C | The models variables are stepped forward with the appropriate
54     C | time-stepping scheme (currently we use Adams-Bashforth II)
55     C | - For momentum, the result is always *only* a "prediction"
56     C | in that the flow may be divergent and will be "corrected"
57     C | later with a surface pressure gradient.
58     C | - Normally for tracers the result is the new field at time
59     C | level [n+1} *BUT* in the case of implicit diffusion the result
60     C | is also *only* a prediction.
61     C | - We denote "predictors" with an asterisk (*).
62     C | U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
63     C | V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
64     C | theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
65     C | salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
66     C | With implicit diffusion:
67     C | theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
68     C | salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
69     C | (1 + dt * K * d_zz) theta[n] = theta*
70     C | (1 + dt * K * d_zz) salt[n] = salt*
71     C |
72     C *==========================================================*
73     C \ev
74    
75     C !USES:
76 adcroft 1.1 IMPLICIT NONE
77     C == Global variables ===
78     #include "SIZE.h"
79     #include "EEPARAMS.h"
80     #include "PARAMS.h"
81     #include "DYNVARS.h"
82     #include "GRID.h"
83 adcroft 1.4 #include "GAD.h"
84 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
85     #include "TR1.h"
86     #endif
87 jmc 1.45 #ifdef ALLOW_PTRACERS
88     #include "PTRACERS.h"
89     #endif
90 heimbach 1.42 #ifdef ALLOW_TIMEAVE
91     #include "TIMEAVE_STATV.h"
92     #endif
93    
94 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
95     # include "tamc.h"
96     # include "tamc_keys.h"
97     # include "FFIELDS.h"
98 heimbach 1.30 # include "EOS.h"
99 adcroft 1.1 # ifdef ALLOW_KPP
100     # include "KPP.h"
101 heimbach 1.67 # endif
102     # ifdef ALLOW_GMREDI
103     # include "GMREDI.h"
104 adcroft 1.1 # endif
105 heimbach 1.68 # ifdef ALLOW_EBM
106     # include "EBM.h"
107 adcroft 1.1 # endif
108     #endif /* ALLOW_AUTODIFF_TAMC */
109    
110 cnh 1.9 C !INPUT/OUTPUT PARAMETERS:
111 adcroft 1.1 C == Routine arguments ==
112     C myTime - Current time in simulation
113     C myIter - Current iteration number in simulation
114     C myThid - Thread number for this instance of the routine.
115     _RL myTime
116     INTEGER myIter
117     INTEGER myThid
118    
119 cnh 1.9 C !LOCAL VARIABLES:
120 adcroft 1.1 C == Local variables
121     C xA, yA - Per block temporaries holding face areas
122     C uTrans, vTrans, rTrans - Per block temporaries holding flow
123     C transport
124     C o uTrans: Zonal transport
125     C o vTrans: Meridional transport
126     C o rTrans: Vertical transport
127 jmc 1.64 C rTransKp1 o vertical volume transp. at interface k+1
128 adcroft 1.1 C maskUp o maskUp: land/water mask for W points
129     C fVer[STUV] o fVer: Vertical flux term - note fVer
130     C is "pipelined" in the vertical
131     C so we need an fVer for each
132     C variable.
133     C rhoK, rhoKM1 - Density at current level, and level above
134     C KappaRT, - Total diffusion in vertical for T and S.
135     C KappaRS (background + spatially varying, isopycnal term).
136 jmc 1.39 C useVariableK = T when vertical diffusion is not constant
137 adcroft 1.1 C iMin, iMax - Ranges and sub-block indices on which calculations
138     C jMin, jMax are applied.
139     C bi, bj
140     C k, kup, - Index for layer above and below. kup and kDown
141     C kDown, km1 are switched with layer to be the appropriate
142     C index into fVerTerm.
143     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148 jmc 1.64 _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149 adcroft 1.1 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
152 heimbach 1.55 #ifdef ALLOW_PASSIVE_TRACER
153 adcroft 1.1 _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
154 heimbach 1.55 #endif
155     #ifdef ALLOW_PTRACERS
156     _RL fVerP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
157     #endif
158 adcroft 1.1 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
160     _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
161     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
162     _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
163     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
164     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
165 jmc 1.64 _RL kp1Msk
166 jmc 1.39 LOGICAL useVariableK
167 adcroft 1.1 INTEGER iMin, iMax
168     INTEGER jMin, jMax
169     INTEGER bi, bj
170     INTEGER i, j
171     INTEGER k, km1, kup, kDown
172 heimbach 1.55 INTEGER iTracer, ip
173 adcroft 1.1
174 cnh 1.9 CEOP
175 adcroft 1.40
176 edhill 1.56 #ifdef ALLOW_DEBUG
177 heimbach 1.43 IF ( debugLevel .GE. debLevB )
178 jmc 1.63 & CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
179 adcroft 1.40 #endif
180 adcroft 1.1
181     #ifdef ALLOW_AUTODIFF_TAMC
182     C-- dummy statement to end declaration part
183     ikey = 1
184 heimbach 1.30 itdkey = 1
185 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
186    
187     #ifdef ALLOW_AUTODIFF_TAMC
188     C-- HPF directive to help TAMC
189     CHPF$ INDEPENDENT
190     #endif /* ALLOW_AUTODIFF_TAMC */
191    
192     DO bj=myByLo(myThid),myByHi(myThid)
193    
194     #ifdef ALLOW_AUTODIFF_TAMC
195     C-- HPF directive to help TAMC
196 heimbach 1.2 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
197 jmc 1.37 CHPF$& ,utrans,vtrans,xA,yA
198 heimbach 1.2 CHPF$& ,KappaRT,KappaRS
199 adcroft 1.1 CHPF$& )
200     #endif /* ALLOW_AUTODIFF_TAMC */
201    
202     DO bi=myBxLo(myThid),myBxHi(myThid)
203    
204     #ifdef ALLOW_AUTODIFF_TAMC
205     act1 = bi - myBxLo(myThid)
206     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
207     act2 = bj - myByLo(myThid)
208     max2 = myByHi(myThid) - myByLo(myThid) + 1
209     act3 = myThid - 1
210     max3 = nTx*nTy
211     act4 = ikey_dynamics - 1
212 heimbach 1.30 itdkey = (act1 + 1) + act2*max1
213 adcroft 1.1 & + act3*max1*max2
214     & + act4*max1*max2*max3
215     #endif /* ALLOW_AUTODIFF_TAMC */
216    
217 heimbach 1.41 C-- Set up work arrays with valid (i.e. not NaN) values
218     C These inital values do not alter the numerical results. They
219     C just ensure that all memory references are to valid floating
220     C point numbers. This prevents spurious hardware signals due to
221     C uninitialised but inert locations.
222    
223 adcroft 1.1 DO j=1-OLy,sNy+OLy
224     DO i=1-OLx,sNx+OLx
225 heimbach 1.41 xA(i,j) = 0. _d 0
226     yA(i,j) = 0. _d 0
227     uTrans(i,j) = 0. _d 0
228     vTrans(i,j) = 0. _d 0
229     rhok (i,j) = 0. _d 0
230     rhoKM1 (i,j) = 0. _d 0
231 adcroft 1.1 rTrans (i,j) = 0. _d 0
232 jmc 1.64 rTransKp1(i,j) = 0. _d 0
233 adcroft 1.1 fVerT (i,j,1) = 0. _d 0
234     fVerT (i,j,2) = 0. _d 0
235     fVerS (i,j,1) = 0. _d 0
236     fVerS (i,j,2) = 0. _d 0
237 heimbach 1.55 #ifdef ALLOW_PASSIVE_TRACER
238 adcroft 1.1 fVerTr1(i,j,1) = 0. _d 0
239     fVerTr1(i,j,2) = 0. _d 0
240 heimbach 1.55 #endif
241     #ifdef ALLOW_PTRACERS
242     DO ip=1,PTRACERS_num
243     fVerP (i,j,1,ip) = 0. _d 0
244     fVerP (i,j,2,ip) = 0. _d 0
245     ENDDO
246     #endif
247 adcroft 1.1 ENDDO
248     ENDDO
249    
250     DO k=1,Nr
251     DO j=1-OLy,sNy+OLy
252     DO i=1-OLx,sNx+OLx
253     C This is currently also used by IVDC and Diagnostics
254 heimbach 1.20 sigmaX(i,j,k) = 0. _d 0
255     sigmaY(i,j,k) = 0. _d 0
256     sigmaR(i,j,k) = 0. _d 0
257 heimbach 1.30 KappaRT(i,j,k) = 0. _d 0
258     KappaRS(i,j,k) = 0. _d 0
259 jmc 1.45 C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
260 heimbach 1.30 gT(i,j,k,bi,bj) = 0. _d 0
261     gS(i,j,k,bi,bj) = 0. _d 0
262     # ifdef ALLOW_PASSIVE_TRACER
263 edhill 1.51 ceh3 needs an IF ( use PASSIVE_TRACER) THEN
264 heimbach 1.5 gTr1(i,j,k,bi,bj) = 0. _d 0
265 heimbach 1.30 # endif
266 heimbach 1.42 # ifdef ALLOW_PTRACERS
267 edhill 1.51 ceh3 this should have an IF ( usePTRACERS ) THEN
268 heimbach 1.42 DO iTracer=1,PTRACERS_numInUse
269     gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
270     ENDDO
271     # endif
272 jmc 1.45 #ifdef ALLOW_AUTODIFF_TAMC
273     cph all the following init. are necessary for TAF
274     cph although some of these are re-initialised later.
275 heimbach 1.30 # ifdef ALLOW_GMREDI
276     Kwx(i,j,k,bi,bj) = 0. _d 0
277     Kwy(i,j,k,bi,bj) = 0. _d 0
278     Kwz(i,j,k,bi,bj) = 0. _d 0
279     # ifdef GM_NON_UNITY_DIAGONAL
280     Kux(i,j,k,bi,bj) = 0. _d 0
281     Kvy(i,j,k,bi,bj) = 0. _d 0
282     # endif
283     # ifdef GM_EXTRA_DIAGONAL
284     Kuz(i,j,k,bi,bj) = 0. _d 0
285     Kvz(i,j,k,bi,bj) = 0. _d 0
286     # endif
287     # ifdef GM_BOLUS_ADVEC
288     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
289     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
290 heimbach 1.36 # endif
291     # ifdef GM_VISBECK_VARIABLE_K
292     VisbeckK(i,j,bi,bj) = 0. _d 0
293 heimbach 1.30 # endif
294     # endif /* ALLOW_GMREDI */
295     #endif /* ALLOW_AUTODIFF_TAMC */
296 adcroft 1.1 ENDDO
297     ENDDO
298     ENDDO
299    
300 jmc 1.21 iMin = 1-OLx
301 adcroft 1.1 iMax = sNx+OLx
302 jmc 1.21 jMin = 1-OLy
303 adcroft 1.1 jMax = sNy+OLy
304    
305     #ifdef ALLOW_AUTODIFF_TAMC
306 heimbach 1.30 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
307     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
308 heimbach 1.38 CADJ STORE totphihyd
309     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
310 heimbach 1.11 #ifdef ALLOW_KPP
311 heimbach 1.30 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
312     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
313 heimbach 1.11 #endif
314 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
315    
316 edhill 1.56 #ifdef ALLOW_DEBUG
317 heimbach 1.43 IF ( debugLevel .GE. debLevB )
318     & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
319 adcroft 1.40 #endif
320    
321 adcroft 1.1 C-- Start of diagnostic loop
322     DO k=Nr,1,-1
323    
324     #ifdef ALLOW_AUTODIFF_TAMC
325     C? Patrick, is this formula correct now that we change the loop range?
326     C? Do we still need this?
327     cph kkey formula corrected.
328     cph Needed for rhok, rhokm1, in the case useGMREDI.
329 heimbach 1.30 kkey = (itdkey-1)*Nr + k
330 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
331    
332     C-- Integrate continuity vertically for vertical velocity
333 jmc 1.27 c CALL INTEGRATE_FOR_W(
334     c I bi, bj, k, uVel, vVel,
335     c O wVel,
336     c I myThid )
337 adcroft 1.1
338     #ifdef ALLOW_OBCS
339     #ifdef ALLOW_NONHYDROSTATIC
340     C-- Apply OBC to W if in N-H mode
341 jmc 1.27 c IF (useOBCS.AND.nonHydrostatic) THEN
342     c CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
343     c ENDIF
344 adcroft 1.1 #endif /* ALLOW_NONHYDROSTATIC */
345     #endif /* ALLOW_OBCS */
346    
347 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
348     C-- MOST of THERMODYNAMICS will be disabled
349     #ifndef SINGLE_LAYER_MODE
350    
351 adcroft 1.1 C-- Calculate gradients of potential density for isoneutral
352     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
353     c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
354     IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
355 edhill 1.56 #ifdef ALLOW_DEBUG
356 heimbach 1.43 IF ( debugLevel .GE. debLevB )
357     & CALL DEBUG_CALL('FIND_RHO',myThid)
358 adcroft 1.40 #endif
359 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
360     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
361     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
362     #endif /* ALLOW_AUTODIFF_TAMC */
363     CALL FIND_RHO(
364 mlosch 1.26 I bi, bj, iMin, iMax, jMin, jMax, k, k,
365 adcroft 1.1 I theta, salt,
366     O rhoK,
367     I myThid )
368 heimbach 1.30
369 adcroft 1.1 IF (k.GT.1) THEN
370     #ifdef ALLOW_AUTODIFF_TAMC
371     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
372     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
373     #endif /* ALLOW_AUTODIFF_TAMC */
374     CALL FIND_RHO(
375 mlosch 1.26 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
376 adcroft 1.1 I theta, salt,
377     O rhoKm1,
378     I myThid )
379     ENDIF
380 edhill 1.56 #ifdef ALLOW_DEBUG
381 heimbach 1.43 IF ( debugLevel .GE. debLevB )
382     & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
383 adcroft 1.40 #endif
384 adcroft 1.1 CALL GRAD_SIGMA(
385     I bi, bj, iMin, iMax, jMin, jMax, k,
386     I rhoK, rhoKm1, rhoK,
387     O sigmaX, sigmaY, sigmaR,
388     I myThid )
389     ENDIF
390    
391 heimbach 1.30 #ifdef ALLOW_AUTODIFF_TAMC
392     CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
393     CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
394     #endif /* ALLOW_AUTODIFF_TAMC */
395 adcroft 1.1 C-- Implicit Vertical Diffusion for Convection
396     c ==> should use sigmaR !!!
397     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
398 edhill 1.56 #ifdef ALLOW_DEBUG
399 heimbach 1.43 IF ( debugLevel .GE. debLevB )
400     & CALL DEBUG_CALL('CALC_IVDC',myThid)
401 adcroft 1.40 #endif
402 adcroft 1.1 CALL CALC_IVDC(
403     I bi, bj, iMin, iMax, jMin, jMax, k,
404     I rhoKm1, rhoK,
405     I myTime, myIter, myThid)
406     ENDIF
407    
408 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
409    
410 adcroft 1.1 C-- end of diagnostic k loop (Nr:1)
411     ENDDO
412    
413     #ifdef ALLOW_AUTODIFF_TAMC
414     cph avoids recomputation of integrate_for_w
415 heimbach 1.30 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
416 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
417    
418     #ifdef ALLOW_OBCS
419     C-- Calculate future values on open boundaries
420     IF (useOBCS) THEN
421 edhill 1.56 #ifdef ALLOW_DEBUG
422 heimbach 1.43 IF ( debugLevel .GE. debLevB )
423     & CALL DEBUG_CALL('OBCS_CALC',myThid)
424 adcroft 1.40 #endif
425 jmc 1.16 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
426 adcroft 1.1 I uVel, vVel, wVel, theta, salt,
427     I myThid )
428     ENDIF
429     #endif /* ALLOW_OBCS */
430    
431 heimbach 1.66 #ifndef ALLOW_AUTODIFF_TAMC
432 jmc 1.62 IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
433 heimbach 1.66 #endif
434 edhill 1.54 C-- Determines forcing terms based on external fields
435     C relaxation terms, etc.
436 edhill 1.56 #ifdef ALLOW_DEBUG
437 edhill 1.54 IF ( debugLevel .GE. debLevB )
438     & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
439     #endif
440 heimbach 1.68 CALL EXTERNAL_FORCING_SURF(
441 edhill 1.54 I bi, bj, iMin, iMax, jMin, jMax,
442     I myTime, myIter, myThid )
443 heimbach 1.66 #ifndef ALLOW_AUTODIFF_TAMC
444     ENDIF
445     #endif
446 edhill 1.54
447 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
448     cph needed for KPP
449     CADJ STORE surfacetendencyU(:,:,bi,bj)
450 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
451 adcroft 1.1 CADJ STORE surfacetendencyV(:,:,bi,bj)
452 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
453 adcroft 1.1 CADJ STORE surfacetendencyS(:,:,bi,bj)
454 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
455 adcroft 1.1 CADJ STORE surfacetendencyT(:,:,bi,bj)
456 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
457 heimbach 1.57 # ifdef ALLOW_SEAICE
458     CADJ STORE surfacetendencyTice(:,:,bi,bj)
459     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
460     # endif
461 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
462    
463 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
464     C-- MOST of THERMODYNAMICS will be disabled
465     #ifndef SINGLE_LAYER_MODE
466    
467 adcroft 1.1 #ifdef ALLOW_GMREDI
468    
469 heimbach 1.22 #ifdef ALLOW_AUTODIFF_TAMC
470 heimbach 1.34 cph storing here is needed only for one GMREDI_OPTIONS:
471     cph define GM_BOLUS_ADVEC
472     cph but I've avoided the #ifdef for now, in case more things change
473     CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
474     CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
475     CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
476 heimbach 1.22 #endif /* ALLOW_AUTODIFF_TAMC */
477 heimbach 1.35
478 adcroft 1.1 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
479     IF (useGMRedi) THEN
480 edhill 1.56 #ifdef ALLOW_DEBUG
481 heimbach 1.43 IF ( debugLevel .GE. debLevB )
482     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
483 adcroft 1.40 #endif
484 jmc 1.15 CALL GMREDI_CALC_TENSOR(
485     I bi, bj, iMin, iMax, jMin, jMax,
486 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
487     I myThid )
488     #ifdef ALLOW_AUTODIFF_TAMC
489     ELSE
490 jmc 1.15 CALL GMREDI_CALC_TENSOR_DUMMY(
491     I bi, bj, iMin, iMax, jMin, jMax,
492 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
493     I myThid )
494     #endif /* ALLOW_AUTODIFF_TAMC */
495     ENDIF
496    
497     #ifdef ALLOW_AUTODIFF_TAMC
498 heimbach 1.30 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
499     CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
500     CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
501 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
502    
503     #endif /* ALLOW_GMREDI */
504    
505     #ifdef ALLOW_KPP
506     C-- Compute KPP mixing coefficients
507     IF (useKPP) THEN
508 edhill 1.56 #ifdef ALLOW_DEBUG
509 heimbach 1.43 IF ( debugLevel .GE. debLevB )
510     & CALL DEBUG_CALL('KPP_CALC',myThid)
511 adcroft 1.40 #endif
512 adcroft 1.1 CALL KPP_CALC(
513     I bi, bj, myTime, myThid )
514     #ifdef ALLOW_AUTODIFF_TAMC
515     ELSE
516     CALL KPP_CALC_DUMMY(
517     I bi, bj, myTime, myThid )
518     #endif /* ALLOW_AUTODIFF_TAMC */
519     ENDIF
520    
521     #ifdef ALLOW_AUTODIFF_TAMC
522     CADJ STORE KPPghat (:,:,:,bi,bj)
523     CADJ & , KPPdiffKzT(:,:,:,bi,bj)
524     CADJ & , KPPdiffKzS(:,:,:,bi,bj)
525     CADJ & , KPPfrac (:,: ,bi,bj)
526 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
527 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
528    
529     #endif /* ALLOW_KPP */
530    
531     #ifdef ALLOW_AUTODIFF_TAMC
532 heimbach 1.30 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
533     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
534     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
535     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
536 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
537 heimbach 1.30 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
538 adcroft 1.1 #endif
539 heimbach 1.42 #ifdef ALLOW_PTRACERS
540     cph-- moved to forward_step to avoid key computation
541     cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
542     cphCADJ & key=itdkey, byte=isbyte
543     #endif
544 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
545    
546     #ifdef ALLOW_AIM
547     C AIM - atmospheric intermediate model, physics package code.
548     IF ( useAIM ) THEN
549 edhill 1.56 #ifdef ALLOW_DEBUG
550 heimbach 1.43 IF ( debugLevel .GE. debLevB )
551     & CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
552 adcroft 1.40 #endif
553 jmc 1.33 CALL TIMER_START('AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
554     CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
555     CALL TIMER_STOP( 'AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
556 adcroft 1.1 ENDIF
557     #endif /* ALLOW_AIM */
558 jmc 1.14
559 adcroft 1.12 #ifndef DISABLE_MULTIDIM_ADVECTION
560 adcroft 1.4 C-- Some advection schemes are better calculated using a multi-dimensional
561     C method in the absence of any other terms and, if used, is done here.
562 adcroft 1.13 C
563     C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
564     C The default is to use multi-dimensinal advection for non-linear advection
565     C schemes. However, for the sake of efficiency of the adjoint it is necessary
566     C to be able to exclude this scheme to avoid excessive storage and
567     C recomputation. It *is* differentiable, if you need it.
568     C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
569     C disable this section of code.
570 jmc 1.24 IF (tempMultiDimAdvec) THEN
571 edhill 1.56 #ifdef ALLOW_DEBUG
572 heimbach 1.43 IF ( debugLevel .GE. debLevB )
573     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
574 adcroft 1.40 #endif
575 jmc 1.63 CALL GAD_ADVECTION(
576 jmc 1.69 I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
577     I GAD_TEMPERATURE,
578 jmc 1.63 I uVel, vVel, wVel, theta,
579     O gT,
580     I bi,bj,myTime,myIter,myThid)
581 jmc 1.23 ENDIF
582 jmc 1.24 IF (saltMultiDimAdvec) THEN
583 edhill 1.56 #ifdef ALLOW_DEBUG
584 heimbach 1.43 IF ( debugLevel .GE. debLevB )
585     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
586 adcroft 1.40 #endif
587 jmc 1.63 CALL GAD_ADVECTION(
588 jmc 1.69 I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
589     I GAD_SALINITY,
590 jmc 1.63 I uVel, vVel, wVel, salt,
591     O gS,
592     I bi,bj,myTime,myIter,myThid)
593 adcroft 1.6 ENDIF
594 adcroft 1.17 C Since passive tracers are configurable separately from T,S we
595     C call the multi-dimensional method for PTRACERS regardless
596     C of whether multiDimAdvection is set or not.
597     #ifdef ALLOW_PTRACERS
598     IF ( usePTRACERS ) THEN
599 edhill 1.56 #ifdef ALLOW_DEBUG
600 heimbach 1.43 IF ( debugLevel .GE. debLevB )
601     & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
602 adcroft 1.40 #endif
603 adcroft 1.17 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
604     ENDIF
605     #endif /* ALLOW_PTRACERS */
606 adcroft 1.12 #endif /* DISABLE_MULTIDIM_ADVECTION */
607 adcroft 1.1
608 edhill 1.56 #ifdef ALLOW_DEBUG
609 jmc 1.63 IF ( debugLevel .GE. debLevB )
610 heimbach 1.43 & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
611 adcroft 1.40 #endif
612    
613 adcroft 1.1 C-- Start of thermodynamics loop
614     DO k=Nr,1,-1
615     #ifdef ALLOW_AUTODIFF_TAMC
616     C? Patrick Is this formula correct?
617     cph Yes, but I rewrote it.
618     cph Also, the KappaR? need the index and subscript k!
619 heimbach 1.30 kkey = (itdkey-1)*Nr + k
620 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
621    
622     C-- km1 Points to level above k (=k-1)
623     C-- kup Cycles through 1,2 to point to layer above
624     C-- kDown Cycles through 2,1 to point to current layer
625    
626     km1 = MAX(1,k-1)
627     kup = 1+MOD(k+1,2)
628     kDown= 1+MOD(k,2)
629    
630     iMin = 1-OLx
631     iMax = sNx+OLx
632     jMin = 1-OLy
633     jMax = sNy+OLy
634    
635 jmc 1.64 kp1Msk=1.
636     IF (k.EQ.Nr) kp1Msk=0.
637     DO j=1-Oly,sNy+Oly
638     DO i=1-Olx,sNx+Olx
639     rTransKp1(i,j) = kp1Msk*rTrans(i,j)
640     ENDDO
641     ENDDO
642 heimbach 1.66 #ifdef ALLOW_AUTODIFF_TAMC
643     CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
644     #endif
645 jmc 1.64
646 adcroft 1.1 C-- Get temporary terms used by tendency routines
647     CALL CALC_COMMON_FACTORS (
648     I bi,bj,iMin,iMax,jMin,jMax,k,
649     O xA,yA,uTrans,vTrans,rTrans,maskUp,
650     I myThid)
651 jmc 1.19
652 jmc 1.64 IF (k.EQ.1) THEN
653     C- Surface interface :
654     DO j=1-Oly,sNy+Oly
655     DO i=1-Olx,sNx+Olx
656     rTrans(i,j) = 0.
657     ENDDO
658     ENDDO
659     ELSE
660     C- Interior interface :
661     DO j=1-Oly,sNy+Oly
662     DO i=1-Olx,sNx+Olx
663     rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
664     ENDDO
665     ENDDO
666     ENDIF
667    
668 jmc 1.19 #ifdef ALLOW_GMREDI
669 heimbach 1.35
670 jmc 1.19 C-- Residual transp = Bolus transp + Eulerian transp
671     IF (useGMRedi) THEN
672     CALL GMREDI_CALC_UVFLOW(
673     & uTrans, vTrans, bi, bj, k, myThid)
674     IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
675     & rTrans, bi, bj, k, myThid)
676     ENDIF
677 heimbach 1.35
678 heimbach 1.66 #ifdef ALLOW_AUTODIFF_TAMC
679     CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
680 heimbach 1.35 #ifdef GM_BOLUS_ADVEC
681     CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
682     CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
683     #endif
684     #endif /* ALLOW_AUTODIFF_TAMC */
685    
686 jmc 1.19 #endif /* ALLOW_GMREDI */
687 adcroft 1.1
688     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
689     C-- Calculate the total vertical diffusivity
690     CALL CALC_DIFFUSIVITY(
691     I bi,bj,iMin,iMax,jMin,jMax,k,
692     I maskUp,
693 heimbach 1.2 O KappaRT,KappaRS,
694 adcroft 1.1 I myThid)
695 heimbach 1.52 # ifdef ALLOW_AUTODIFF_TAMC
696     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
697     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
698     # endif /* ALLOW_AUTODIFF_TAMC */
699 adcroft 1.1 #endif
700    
701     iMin = 1-OLx+2
702     iMax = sNx+OLx-1
703     jMin = 1-OLy+2
704     jMax = sNy+OLy-1
705    
706     C-- Calculate active tracer tendencies (gT,gS,...)
707     C and step forward storing result in gTnm1, gSnm1, etc.
708     IF ( tempStepping ) THEN
709     CALL CALC_GT(
710     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
711 jmc 1.64 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
712 adcroft 1.1 I KappaRT,
713     U fVerT,
714 adcroft 1.7 I myTime,myIter,myThid)
715 adcroft 1.1 CALL TIMESTEP_TRACER(
716 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
717 adcroft 1.1 I theta, gT,
718     I myIter, myThid)
719     ENDIF
720 jmc 1.44
721 adcroft 1.1 IF ( saltStepping ) THEN
722     CALL CALC_GS(
723     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
724 jmc 1.64 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
725 adcroft 1.1 I KappaRS,
726     U fVerS,
727 adcroft 1.7 I myTime,myIter,myThid)
728 adcroft 1.1 CALL TIMESTEP_TRACER(
729 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
730 adcroft 1.1 I salt, gS,
731     I myIter, myThid)
732     ENDIF
733     #ifdef ALLOW_PASSIVE_TRACER
734 edhill 1.51 ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
735 adcroft 1.1 IF ( tr1Stepping ) THEN
736     CALL CALC_GTR1(
737     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
738 jmc 1.64 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
739 adcroft 1.1 I KappaRT,
740     U fVerTr1,
741 heimbach 1.8 I myTime,myIter,myThid)
742 adcroft 1.1 CALL TIMESTEP_TRACER(
743 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
744 adcroft 1.1 I Tr1, gTr1,
745 heimbach 1.8 I myIter,myThid)
746 adcroft 1.1 ENDIF
747     #endif
748 adcroft 1.17 #ifdef ALLOW_PTRACERS
749     IF ( usePTRACERS ) THEN
750 heimbach 1.42 CALL PTRACERS_INTEGRATE(
751 adcroft 1.17 I bi,bj,k,
752 jmc 1.64 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
753 heimbach 1.55 X fVerP, KappaRS,
754 adcroft 1.17 I myIter,myTime,myThid)
755     ENDIF
756     #endif /* ALLOW_PTRACERS */
757 adcroft 1.1
758     #ifdef ALLOW_OBCS
759     C-- Apply open boundary conditions
760     IF (useOBCS) THEN
761 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
762 adcroft 1.1 END IF
763     #endif /* ALLOW_OBCS */
764 edhill 1.54
765 jmc 1.59 C-- Freeze water
766     C this bit of code is left here for backward compatibility.
767     C freezing at surface level has been moved to FORWARD_STEP
768     IF ( useOldFreezing .AND. .NOT. useSEAICE
769 jmc 1.61 & .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
770 jmc 1.59 #ifdef ALLOW_AUTODIFF_TAMC
771     CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
772     CADJ & , key = kkey, byte = isbyte
773     #endif /* ALLOW_AUTODIFF_TAMC */
774     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
775     ENDIF
776 adcroft 1.1
777     C-- end of thermodynamic k loop (Nr:1)
778     ENDDO
779 cheisey 1.31
780 adcroft 1.1
781 jmc 1.63 C-- Implicit vertical advection & diffusion
782     #ifdef INCLUDE_IMPLVERTADV_CODE
783     IF ( tempImplVertAdv ) THEN
784     CALL GAD_IMPLICIT_R(
785     I tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
786     I kappaRT, wVel, theta,
787     U gT,
788     I bi, bj, myTime, myIter, myThid )
789     ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
790     #else /* INCLUDE_IMPLVERTADV_CODE */
791     IF ( tempStepping .AND. implicitDiffusion ) THEN
792     #endif /* INCLUDE_IMPLVERTADV_CODE */
793 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
794 heimbach 1.52 CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
795 heimbach 1.30 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
796 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
797 jmc 1.63 CALL IMPLDIFF(
798 adcroft 1.1 I bi, bj, iMin, iMax, jMin, jMax,
799     I deltaTtracer, KappaRT, recip_HFacC,
800 adcroft 1.7 U gT,
801 adcroft 1.1 I myThid )
802 jmc 1.63 ENDIF
803 adcroft 1.1
804 jmc 1.63 #ifdef INCLUDE_IMPLVERTADV_CODE
805     IF ( saltImplVertAdv ) THEN
806     CALL GAD_IMPLICIT_R(
807     I saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
808     I kappaRS, wVel, salt,
809     U gS,
810     I bi, bj, myTime, myIter, myThid )
811     ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
812     #else /* INCLUDE_IMPLVERTADV_CODE */
813     IF ( saltStepping .AND. implicitDiffusion ) THEN
814     #endif /* INCLUDE_IMPLVERTADV_CODE */
815 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
816 heimbach 1.52 CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
817 heimbach 1.30 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
818 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
819 jmc 1.63 CALL IMPLDIFF(
820 adcroft 1.1 I bi, bj, iMin, iMax, jMin, jMax,
821     I deltaTtracer, KappaRS, recip_HFacC,
822 adcroft 1.7 U gS,
823 adcroft 1.1 I myThid )
824 jmc 1.63 ENDIF
825 adcroft 1.1
826     #ifdef ALLOW_PASSIVE_TRACER
827 jmc 1.63 IF ( tr1Stepping .AND. implicitDiffusion ) THEN
828 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
829 heimbach 1.30 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
830 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
831     CALL IMPLDIFF(
832     I bi, bj, iMin, iMax, jMin, jMax,
833     I deltaTtracer, KappaRT, recip_HFacC,
834 adcroft 1.7 U gTr1,
835 adcroft 1.1 I myThid )
836 jmc 1.63 ENDIF
837 adcroft 1.1 #endif
838    
839 adcroft 1.17 #ifdef ALLOW_PTRACERS
840 jmc 1.63 c #ifdef INCLUDE_IMPLVERTADV_CODE
841     c IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN
842     c ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN
843     c #else
844     IF ( usePTRACERS .AND. implicitDiffusion ) THEN
845     C-- Vertical diffusion (implicit) for passive tracers
846 adcroft 1.17 CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
847 jmc 1.63 ENDIF
848 adcroft 1.17 #endif /* ALLOW_PTRACERS */
849    
850 adcroft 1.1 #ifdef ALLOW_OBCS
851     C-- Apply open boundary conditions
852 jmc 1.63 IF ( ( implicitDiffusion
853     & .OR. tempImplVertAdv
854     & .OR. saltImplVertAdv
855     & ) .AND. useOBCS ) THEN
856 adcroft 1.1 DO K=1,Nr
857 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
858 adcroft 1.1 ENDDO
859 jmc 1.63 ENDIF
860 adcroft 1.1 #endif /* ALLOW_OBCS */
861    
862 jmc 1.39 #ifdef ALLOW_TIMEAVE
863 dimitri 1.65 #ifndef HRCUBE
864 jmc 1.39 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
865 jmc 1.70 CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
866 jmc 1.39 I Nr, deltaTclock, bi, bj, myThid)
867     ENDIF
868     useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
869     IF (taveFreq.GT.0. .AND. useVariableK ) THEN
870     IF (implicitDiffusion) THEN
871     CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
872     I Nr, 3, deltaTclock, bi, bj, myThid)
873     ELSE
874     CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
875     I Nr, 3, deltaTclock, bi, bj, myThid)
876     ENDIF
877     ENDIF
878 dimitri 1.65 #endif /* ndef HRCUBE */
879 jmc 1.39 #endif /* ALLOW_TIMEAVE */
880    
881 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
882 adcroft 1.1
883 jmc 1.39 C-- end bi,bj loops.
884 adcroft 1.1 ENDDO
885     ENDDO
886 adcroft 1.17
887 edhill 1.56 #ifdef ALLOW_DEBUG
888 adcroft 1.17 If (debugMode) THEN
889     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
890     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
891     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
892     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
893     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
894     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
895     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
896     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
897     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
898 adcroft 1.18 #ifdef ALLOW_PTRACERS
899     IF ( usePTRACERS ) THEN
900     CALL PTRACERS_DEBUG(myThid)
901     ENDIF
902     #endif /* ALLOW_PTRACERS */
903 adcroft 1.17 ENDIF
904 adcroft 1.40 #endif
905    
906 edhill 1.56 #ifdef ALLOW_DEBUG
907 heimbach 1.43 IF ( debugLevel .GE. debLevB )
908 jmc 1.63 & CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
909 adcroft 1.17 #endif
910 adcroft 1.1
911     RETURN
912     END
  ViewVC Help
Powered by ViewVC 1.1.22