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