/[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.59 - (hide annotations) (download)
Thu Nov 13 21:54:11 2003 UTC (20 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint52a_post
Changes since 1.58: +12 -12 lines 
additional changes for FREEZE:
 - new S/R FREEZE_SURFACE only apllied to surface level.
 - add run-time parameter "useOldFreezing" to use the old version "FREEZE"
1 jmc 1.59 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.58 2003/11/13 06:35:14 dimitri 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 edhill 1.51
7 jmc 1.21 #ifdef ALLOW_AUTODIFF_TAMC
8     # ifdef ALLOW_GMREDI
9     # include "GMREDI_OPTIONS.h"
10     # endif
11     # ifdef ALLOW_KPP
12     # include "KPP_OPTIONS.h"
13 heimbach 1.42 # endif
14 dimitri 1.48 #ifdef ALLOW_PTRACERS
15     # include "PTRACERS_OPTIONS.h"
16     #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     # endif
102     # ifdef ALLOW_GMREDI
103     # include "GMREDI.h"
104     # endif
105     #endif /* ALLOW_AUTODIFF_TAMC */
106    
107 cnh 1.9 C !INPUT/OUTPUT PARAMETERS:
108 adcroft 1.1 C == Routine arguments ==
109     C myTime - Current time in simulation
110     C myIter - Current iteration number in simulation
111     C myThid - Thread number for this instance of the routine.
112     _RL myTime
113     INTEGER myIter
114     INTEGER myThid
115    
116 cnh 1.9 C !LOCAL VARIABLES:
117 adcroft 1.1 C == Local variables
118     C xA, yA - Per block temporaries holding face areas
119     C uTrans, vTrans, rTrans - Per block temporaries holding flow
120     C transport
121     C o uTrans: Zonal transport
122     C o vTrans: Meridional transport
123     C o rTrans: Vertical transport
124     C maskUp o maskUp: land/water mask for W points
125     C fVer[STUV] o fVer: Vertical flux term - note fVer
126     C is "pipelined" in the vertical
127     C so we need an fVer for each
128     C variable.
129     C rhoK, rhoKM1 - Density at current level, and level above
130     C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
131     C phiSurfY or geopotentiel (atmos) in X and Y direction
132     C KappaRT, - Total diffusion in vertical for T and S.
133     C KappaRS (background + spatially varying, isopycnal term).
134 jmc 1.39 C useVariableK = T when vertical diffusion is not constant
135 adcroft 1.1 C iMin, iMax - Ranges and sub-block indices on which calculations
136     C jMin, jMax are applied.
137     C bi, bj
138     C k, kup, - Index for layer above and below. kup and kDown
139     C kDown, km1 are switched with layer to be the appropriate
140     C index into fVerTerm.
141     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146     _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
148     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
149 heimbach 1.55 #ifdef ALLOW_PASSIVE_TRACER
150 adcroft 1.1 _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151 heimbach 1.55 #endif
152     #ifdef ALLOW_PTRACERS
153     _RL fVerP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
154     #endif
155 adcroft 1.1 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157     _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158     _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159     _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
160     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
161     _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
162     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
163     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
164 cnh 1.9 C This is currently used by IVDC and Diagnostics
165 adcroft 1.1 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
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     & CALL DEBUG_ENTER('FORWARD_STEP',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     phiSurfX(i,j) = 0. _d 0
232     phiSurfY(i,j) = 0. _d 0
233 adcroft 1.1 rTrans (i,j) = 0. _d 0
234     fVerT (i,j,1) = 0. _d 0
235     fVerT (i,j,2) = 0. _d 0
236     fVerS (i,j,1) = 0. _d 0
237     fVerS (i,j,2) = 0. _d 0
238 heimbach 1.55 #ifdef ALLOW_PASSIVE_TRACER
239 adcroft 1.1 fVerTr1(i,j,1) = 0. _d 0
240     fVerTr1(i,j,2) = 0. _d 0
241 heimbach 1.55 #endif
242     #ifdef ALLOW_PTRACERS
243     DO ip=1,PTRACERS_num
244     fVerP (i,j,1,ip) = 0. _d 0
245     fVerP (i,j,2,ip) = 0. _d 0
246     ENDDO
247     #endif
248 adcroft 1.1 ENDDO
249     ENDDO
250    
251     DO k=1,Nr
252     DO j=1-OLy,sNy+OLy
253     DO i=1-OLx,sNx+OLx
254     C This is currently also used by IVDC and Diagnostics
255 heimbach 1.20 sigmaX(i,j,k) = 0. _d 0
256     sigmaY(i,j,k) = 0. _d 0
257     sigmaR(i,j,k) = 0. _d 0
258 adcroft 1.1 ConvectCount(i,j,k) = 0.
259 heimbach 1.30 KappaRT(i,j,k) = 0. _d 0
260     KappaRS(i,j,k) = 0. _d 0
261 jmc 1.45 C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
262 heimbach 1.30 gT(i,j,k,bi,bj) = 0. _d 0
263     gS(i,j,k,bi,bj) = 0. _d 0
264     # ifdef ALLOW_PASSIVE_TRACER
265 edhill 1.51 ceh3 needs an IF ( use PASSIVE_TRACER) THEN
266 heimbach 1.5 gTr1(i,j,k,bi,bj) = 0. _d 0
267 heimbach 1.30 # endif
268 heimbach 1.42 # ifdef ALLOW_PTRACERS
269 edhill 1.51 ceh3 this should have an IF ( usePTRACERS ) THEN
270 heimbach 1.42 DO iTracer=1,PTRACERS_numInUse
271     gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
272     ENDDO
273     # endif
274 jmc 1.45 #ifdef ALLOW_AUTODIFF_TAMC
275     cph all the following init. are necessary for TAF
276     cph although some of these are re-initialised later.
277 heimbach 1.30 # ifdef ALLOW_GMREDI
278     Kwx(i,j,k,bi,bj) = 0. _d 0
279     Kwy(i,j,k,bi,bj) = 0. _d 0
280     Kwz(i,j,k,bi,bj) = 0. _d 0
281     # ifdef GM_NON_UNITY_DIAGONAL
282     Kux(i,j,k,bi,bj) = 0. _d 0
283     Kvy(i,j,k,bi,bj) = 0. _d 0
284     # endif
285     # ifdef GM_EXTRA_DIAGONAL
286     Kuz(i,j,k,bi,bj) = 0. _d 0
287     Kvz(i,j,k,bi,bj) = 0. _d 0
288     # endif
289     # ifdef GM_BOLUS_ADVEC
290     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
291     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
292 heimbach 1.36 # endif
293     # ifdef GM_VISBECK_VARIABLE_K
294     VisbeckK(i,j,bi,bj) = 0. _d 0
295 heimbach 1.30 # endif
296     # endif /* ALLOW_GMREDI */
297     #endif /* ALLOW_AUTODIFF_TAMC */
298 adcroft 1.1 ENDDO
299     ENDDO
300     ENDDO
301    
302 jmc 1.21 iMin = 1-OLx
303 adcroft 1.1 iMax = sNx+OLx
304 jmc 1.21 jMin = 1-OLy
305 adcroft 1.1 jMax = sNy+OLy
306    
307     #ifdef ALLOW_AUTODIFF_TAMC
308 heimbach 1.30 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
309     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
310 heimbach 1.38 CADJ STORE totphihyd
311     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
312 heimbach 1.11 #ifdef ALLOW_KPP
313 heimbach 1.30 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
314     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
315 heimbach 1.11 #endif
316 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
317    
318 edhill 1.56 #ifdef ALLOW_DEBUG
319 heimbach 1.43 IF ( debugLevel .GE. debLevB )
320     & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
321 adcroft 1.40 #endif
322    
323 adcroft 1.1 C-- Start of diagnostic loop
324     DO k=Nr,1,-1
325    
326     #ifdef ALLOW_AUTODIFF_TAMC
327     C? Patrick, is this formula correct now that we change the loop range?
328     C? Do we still need this?
329     cph kkey formula corrected.
330     cph Needed for rhok, rhokm1, in the case useGMREDI.
331 heimbach 1.30 kkey = (itdkey-1)*Nr + k
332 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
333    
334     C-- Integrate continuity vertically for vertical velocity
335 jmc 1.27 c CALL INTEGRATE_FOR_W(
336     c I bi, bj, k, uVel, vVel,
337     c O wVel,
338     c I myThid )
339 adcroft 1.1
340     #ifdef ALLOW_OBCS
341     #ifdef ALLOW_NONHYDROSTATIC
342     C-- Apply OBC to W if in N-H mode
343 jmc 1.27 c IF (useOBCS.AND.nonHydrostatic) THEN
344     c CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
345     c ENDIF
346 adcroft 1.1 #endif /* ALLOW_NONHYDROSTATIC */
347     #endif /* ALLOW_OBCS */
348    
349 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
350     C-- MOST of THERMODYNAMICS will be disabled
351     #ifndef SINGLE_LAYER_MODE
352    
353 adcroft 1.1 C-- Calculate gradients of potential density for isoneutral
354     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
355     c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
356     IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
357 edhill 1.56 #ifdef ALLOW_DEBUG
358 heimbach 1.43 IF ( debugLevel .GE. debLevB )
359     & CALL DEBUG_CALL('FIND_RHO',myThid)
360 adcroft 1.40 #endif
361 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
362     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
363     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
364     #endif /* ALLOW_AUTODIFF_TAMC */
365     CALL FIND_RHO(
366 mlosch 1.26 I bi, bj, iMin, iMax, jMin, jMax, k, k,
367 adcroft 1.1 I theta, salt,
368     O rhoK,
369     I myThid )
370 heimbach 1.30
371 adcroft 1.1 IF (k.GT.1) THEN
372     #ifdef ALLOW_AUTODIFF_TAMC
373     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
374     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
375     #endif /* ALLOW_AUTODIFF_TAMC */
376     CALL FIND_RHO(
377 mlosch 1.26 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
378 adcroft 1.1 I theta, salt,
379     O rhoKm1,
380     I myThid )
381     ENDIF
382 edhill 1.56 #ifdef ALLOW_DEBUG
383 heimbach 1.43 IF ( debugLevel .GE. debLevB )
384     & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
385 adcroft 1.40 #endif
386 adcroft 1.1 CALL GRAD_SIGMA(
387     I bi, bj, iMin, iMax, jMin, jMax, k,
388     I rhoK, rhoKm1, rhoK,
389     O sigmaX, sigmaY, sigmaR,
390     I myThid )
391     ENDIF
392    
393 heimbach 1.30 #ifdef ALLOW_AUTODIFF_TAMC
394     CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
395     CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
396     #endif /* ALLOW_AUTODIFF_TAMC */
397 adcroft 1.1 C-- Implicit Vertical Diffusion for Convection
398     c ==> should use sigmaR !!!
399     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
400 edhill 1.56 #ifdef ALLOW_DEBUG
401 heimbach 1.43 IF ( debugLevel .GE. debLevB )
402     & CALL DEBUG_CALL('CALC_IVDC',myThid)
403 adcroft 1.40 #endif
404 adcroft 1.1 CALL CALC_IVDC(
405     I bi, bj, iMin, iMax, jMin, jMax, k,
406     I rhoKm1, rhoK,
407     U ConvectCount, KappaRT, KappaRS,
408     I myTime, myIter, myThid)
409     ENDIF
410    
411 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
412    
413 adcroft 1.1 C-- end of diagnostic k loop (Nr:1)
414     ENDDO
415    
416     #ifdef ALLOW_AUTODIFF_TAMC
417     cph avoids recomputation of integrate_for_w
418 heimbach 1.30 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
419 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
420    
421     #ifdef ALLOW_OBCS
422     C-- Calculate future values on open boundaries
423     IF (useOBCS) THEN
424 edhill 1.56 #ifdef ALLOW_DEBUG
425 heimbach 1.43 IF ( debugLevel .GE. debLevB )
426     & CALL DEBUG_CALL('OBCS_CALC',myThid)
427 adcroft 1.40 #endif
428 jmc 1.16 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
429 adcroft 1.1 I uVel, vVel, wVel, theta, salt,
430     I myThid )
431     ENDIF
432     #endif /* ALLOW_OBCS */
433    
434 edhill 1.54
435 cheisey 1.32 #ifdef ALLOW_THERM_SEAICE
436 jmc 1.44 IF (useThermSeaIce) THEN
437 edhill 1.56 #ifdef ALLOW_DEBUG
438 heimbach 1.43 IF ( debugLevel .GE. debLevB )
439     & CALL DEBUG_CALL('ICE_FORCING',myThid)
440 adcroft 1.40 #endif
441 cheisey 1.31 C-- Determines forcing terms based on external fields
442 jmc 1.44 C including effects from ice
443 cheisey 1.31 CALL ICE_FORCING(
444     I bi, bj, iMin, iMax, jMin, jMax,
445     I myThid )
446 edhill 1.54 ELSE
447     #endif /* ALLOW_THERM_SEAICE */
448    
449     C-- Determines forcing terms based on external fields
450     C relaxation terms, etc.
451 edhill 1.56 #ifdef ALLOW_DEBUG
452 edhill 1.54 IF ( debugLevel .GE. debLevB )
453     & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
454     #endif
455     CALL EXTERNAL_FORCING_SURF(
456     I bi, bj, iMin, iMax, jMin, jMax,
457     I myTime, myIter, myThid )
458    
459     #ifdef ALLOW_THERM_SEAICE
460     C-- end of if/else block useThermSeaIce --
461 jmc 1.50 ENDIF
462     #endif /* ALLOW_THERM_SEAICE */
463 cheisey 1.31
464 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
465     cph needed for KPP
466     CADJ STORE surfacetendencyU(:,:,bi,bj)
467 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
468 adcroft 1.1 CADJ STORE surfacetendencyV(:,:,bi,bj)
469 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
470 adcroft 1.1 CADJ STORE surfacetendencyS(:,:,bi,bj)
471 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
472 adcroft 1.1 CADJ STORE surfacetendencyT(:,:,bi,bj)
473 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
474 heimbach 1.57 # ifdef ALLOW_SEAICE
475     CADJ STORE surfacetendencyTice(:,:,bi,bj)
476     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
477     # endif
478 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
479    
480 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
481     C-- MOST of THERMODYNAMICS will be disabled
482     #ifndef SINGLE_LAYER_MODE
483    
484 adcroft 1.1 #ifdef ALLOW_GMREDI
485    
486 heimbach 1.22 #ifdef ALLOW_AUTODIFF_TAMC
487 heimbach 1.34 cph storing here is needed only for one GMREDI_OPTIONS:
488     cph define GM_BOLUS_ADVEC
489     cph but I've avoided the #ifdef for now, in case more things change
490     CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
491     CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
492     CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
493 heimbach 1.22 #endif /* ALLOW_AUTODIFF_TAMC */
494 heimbach 1.35
495 adcroft 1.1 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
496     IF (useGMRedi) THEN
497 edhill 1.56 #ifdef ALLOW_DEBUG
498 heimbach 1.43 IF ( debugLevel .GE. debLevB )
499     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
500 adcroft 1.40 #endif
501 jmc 1.15 CALL GMREDI_CALC_TENSOR(
502     I bi, bj, iMin, iMax, jMin, jMax,
503 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
504     I myThid )
505     #ifdef ALLOW_AUTODIFF_TAMC
506     ELSE
507 jmc 1.15 CALL GMREDI_CALC_TENSOR_DUMMY(
508     I bi, bj, iMin, iMax, jMin, jMax,
509 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
510     I myThid )
511     #endif /* ALLOW_AUTODIFF_TAMC */
512     ENDIF
513    
514     #ifdef ALLOW_AUTODIFF_TAMC
515 heimbach 1.30 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
516     CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
517     CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
518 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
519    
520     #endif /* ALLOW_GMREDI */
521    
522     #ifdef ALLOW_KPP
523     C-- Compute KPP mixing coefficients
524     IF (useKPP) THEN
525 edhill 1.56 #ifdef ALLOW_DEBUG
526 heimbach 1.43 IF ( debugLevel .GE. debLevB )
527     & CALL DEBUG_CALL('KPP_CALC',myThid)
528 adcroft 1.40 #endif
529 adcroft 1.1 CALL KPP_CALC(
530     I bi, bj, myTime, myThid )
531     #ifdef ALLOW_AUTODIFF_TAMC
532     ELSE
533     CALL KPP_CALC_DUMMY(
534     I bi, bj, myTime, myThid )
535     #endif /* ALLOW_AUTODIFF_TAMC */
536     ENDIF
537    
538     #ifdef ALLOW_AUTODIFF_TAMC
539     CADJ STORE KPPghat (:,:,:,bi,bj)
540     CADJ & , KPPdiffKzT(:,:,:,bi,bj)
541     CADJ & , KPPdiffKzS(:,:,:,bi,bj)
542     CADJ & , KPPfrac (:,: ,bi,bj)
543 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
544 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
545    
546     #endif /* ALLOW_KPP */
547    
548     #ifdef ALLOW_AUTODIFF_TAMC
549 heimbach 1.30 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
550     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
551     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
552     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
553 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
554 heimbach 1.30 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
555 adcroft 1.1 #endif
556 heimbach 1.42 #ifdef ALLOW_PTRACERS
557     cph-- moved to forward_step to avoid key computation
558     cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
559     cphCADJ & key=itdkey, byte=isbyte
560     #endif
561 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
562    
563     #ifdef ALLOW_AIM
564     C AIM - atmospheric intermediate model, physics package code.
565     IF ( useAIM ) THEN
566 edhill 1.56 #ifdef ALLOW_DEBUG
567 heimbach 1.43 IF ( debugLevel .GE. debLevB )
568     & CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
569 adcroft 1.40 #endif
570 jmc 1.33 CALL TIMER_START('AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
571     CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
572     CALL TIMER_STOP( 'AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
573 adcroft 1.1 ENDIF
574     #endif /* ALLOW_AIM */
575 jmc 1.14
576 adcroft 1.12 #ifndef DISABLE_MULTIDIM_ADVECTION
577 adcroft 1.4 C-- Some advection schemes are better calculated using a multi-dimensional
578     C method in the absence of any other terms and, if used, is done here.
579 adcroft 1.13 C
580     C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
581     C The default is to use multi-dimensinal advection for non-linear advection
582     C schemes. However, for the sake of efficiency of the adjoint it is necessary
583     C to be able to exclude this scheme to avoid excessive storage and
584     C recomputation. It *is* differentiable, if you need it.
585     C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
586     C disable this section of code.
587 jmc 1.24 IF (tempMultiDimAdvec) THEN
588 edhill 1.56 #ifdef ALLOW_DEBUG
589 heimbach 1.43 IF ( debugLevel .GE. debLevB )
590     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
591 adcroft 1.40 #endif
592 heimbach 1.11 CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
593     U theta,gT,
594 adcroft 1.6 I myTime,myIter,myThid)
595 jmc 1.23 ENDIF
596 jmc 1.24 IF (saltMultiDimAdvec) THEN
597 edhill 1.56 #ifdef ALLOW_DEBUG
598 heimbach 1.43 IF ( debugLevel .GE. debLevB )
599     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
600 adcroft 1.40 #endif
601 heimbach 1.11 CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
602     U salt,gS,
603 adcroft 1.6 I myTime,myIter,myThid)
604     ENDIF
605 adcroft 1.17 C Since passive tracers are configurable separately from T,S we
606     C call the multi-dimensional method for PTRACERS regardless
607     C of whether multiDimAdvection is set or not.
608     #ifdef ALLOW_PTRACERS
609     IF ( usePTRACERS ) THEN
610 edhill 1.56 #ifdef ALLOW_DEBUG
611 heimbach 1.43 IF ( debugLevel .GE. debLevB )
612     & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
613 adcroft 1.40 #endif
614 adcroft 1.17 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
615     ENDIF
616     #endif /* ALLOW_PTRACERS */
617 adcroft 1.12 #endif /* DISABLE_MULTIDIM_ADVECTION */
618 adcroft 1.1
619 edhill 1.56 #ifdef ALLOW_DEBUG
620 heimbach 1.43 IF ( debugLevel .GE. debLevB )
621     & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
622 adcroft 1.40 #endif
623    
624 adcroft 1.1 C-- Start of thermodynamics loop
625     DO k=Nr,1,-1
626     #ifdef ALLOW_AUTODIFF_TAMC
627     C? Patrick Is this formula correct?
628     cph Yes, but I rewrote it.
629     cph Also, the KappaR? need the index and subscript k!
630 heimbach 1.30 kkey = (itdkey-1)*Nr + k
631 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
632    
633     C-- km1 Points to level above k (=k-1)
634     C-- kup Cycles through 1,2 to point to layer above
635     C-- kDown Cycles through 2,1 to point to current layer
636    
637     km1 = MAX(1,k-1)
638     kup = 1+MOD(k+1,2)
639     kDown= 1+MOD(k,2)
640    
641     iMin = 1-OLx
642     iMax = sNx+OLx
643     jMin = 1-OLy
644     jMax = sNy+OLy
645    
646     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     #ifdef ALLOW_GMREDI
653 heimbach 1.35
654 jmc 1.19 C-- Residual transp = Bolus transp + Eulerian transp
655     IF (useGMRedi) THEN
656     CALL GMREDI_CALC_UVFLOW(
657     & uTrans, vTrans, bi, bj, k, myThid)
658     IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
659     & rTrans, bi, bj, k, myThid)
660     ENDIF
661 heimbach 1.35
662     #ifdef ALLOW_AUTODIFF_TAMC
663     #ifdef GM_BOLUS_ADVEC
664     CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
665     CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
666     CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
667     #endif
668     #endif /* ALLOW_AUTODIFF_TAMC */
669    
670 jmc 1.19 #endif /* ALLOW_GMREDI */
671 adcroft 1.1
672     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
673     C-- Calculate the total vertical diffusivity
674     CALL CALC_DIFFUSIVITY(
675     I bi,bj,iMin,iMax,jMin,jMax,k,
676     I maskUp,
677 heimbach 1.2 O KappaRT,KappaRS,
678 adcroft 1.1 I myThid)
679 heimbach 1.52 # ifdef ALLOW_AUTODIFF_TAMC
680     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
681     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
682     # endif /* ALLOW_AUTODIFF_TAMC */
683 adcroft 1.1 #endif
684    
685     iMin = 1-OLx+2
686     iMax = sNx+OLx-1
687     jMin = 1-OLy+2
688     jMax = sNy+OLy-1
689    
690     C-- Calculate active tracer tendencies (gT,gS,...)
691     C and step forward storing result in gTnm1, gSnm1, etc.
692     IF ( tempStepping ) THEN
693     CALL CALC_GT(
694     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
695     I xA,yA,uTrans,vTrans,rTrans,maskUp,
696     I KappaRT,
697     U fVerT,
698 adcroft 1.7 I myTime,myIter,myThid)
699 adcroft 1.1 CALL TIMESTEP_TRACER(
700 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
701 adcroft 1.1 I theta, gT,
702     I myIter, myThid)
703     ENDIF
704 jmc 1.44
705 cheisey 1.32 #ifdef ALLOW_THERM_SEAICE
706 jmc 1.44 IF (useThermSeaIce .AND. k.EQ.1) THEN
707     CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )
708     ENDIF
709 cheisey 1.32 #endif
710 jmc 1.44
711 adcroft 1.1 IF ( saltStepping ) THEN
712     CALL CALC_GS(
713     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
714     I xA,yA,uTrans,vTrans,rTrans,maskUp,
715     I KappaRS,
716     U fVerS,
717 adcroft 1.7 I myTime,myIter,myThid)
718 adcroft 1.1 CALL TIMESTEP_TRACER(
719 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
720 adcroft 1.1 I salt, gS,
721     I myIter, myThid)
722     ENDIF
723     #ifdef ALLOW_PASSIVE_TRACER
724 edhill 1.51 ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
725 adcroft 1.1 IF ( tr1Stepping ) THEN
726     CALL CALC_GTR1(
727     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
728     I xA,yA,uTrans,vTrans,rTrans,maskUp,
729     I KappaRT,
730     U fVerTr1,
731 heimbach 1.8 I myTime,myIter,myThid)
732 adcroft 1.1 CALL TIMESTEP_TRACER(
733 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
734 adcroft 1.1 I Tr1, gTr1,
735 heimbach 1.8 I myIter,myThid)
736 adcroft 1.1 ENDIF
737     #endif
738 adcroft 1.17 #ifdef ALLOW_PTRACERS
739     IF ( usePTRACERS ) THEN
740 heimbach 1.42 CALL PTRACERS_INTEGRATE(
741 adcroft 1.17 I bi,bj,k,
742     I xA,yA,uTrans,vTrans,rTrans,maskUp,
743 heimbach 1.55 X fVerP, KappaRS,
744 adcroft 1.17 I myIter,myTime,myThid)
745     ENDIF
746     #endif /* ALLOW_PTRACERS */
747 adcroft 1.1
748     #ifdef ALLOW_OBCS
749     C-- Apply open boundary conditions
750     IF (useOBCS) THEN
751 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
752 adcroft 1.1 END IF
753     #endif /* ALLOW_OBCS */
754 edhill 1.54
755 jmc 1.59 C-- Freeze water
756     C this bit of code is left here for backward compatibility.
757     C freezing at surface level has been moved to FORWARD_STEP
758     IF ( useOldFreezing .AND. .NOT. useSEAICE
759     & .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN
760     #ifdef ALLOW_AUTODIFF_TAMC
761     CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
762     CADJ & , key = kkey, byte = isbyte
763     #endif /* ALLOW_AUTODIFF_TAMC */
764     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
765     ENDIF
766 adcroft 1.1
767     C-- end of thermodynamic k loop (Nr:1)
768     ENDDO
769 cheisey 1.31
770     cswdice -- add ---
771 cheisey 1.32 #ifdef ALLOW_THERM_SEAICE
772 cheisey 1.31 c timeaveraging for ice model values
773 edhill 1.51 ceh3 This should be wrapped in an IF ( useThermSeaIce ) THEN
774 cheisey 1.31 CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
775     #endif
776     cswdice --- end add ---
777    
778    
779    
780 adcroft 1.1
781     C-- Implicit diffusion
782     IF (implicitDiffusion) THEN
783    
784     IF (tempStepping) THEN
785     #ifdef ALLOW_AUTODIFF_TAMC
786 heimbach 1.52 CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
787 heimbach 1.30 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
788 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
789     CALL IMPLDIFF(
790     I bi, bj, iMin, iMax, jMin, jMax,
791     I deltaTtracer, KappaRT, recip_HFacC,
792 adcroft 1.7 U gT,
793 adcroft 1.1 I myThid )
794     ENDIF
795    
796     IF (saltStepping) THEN
797     #ifdef ALLOW_AUTODIFF_TAMC
798 heimbach 1.52 CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
799 heimbach 1.30 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
800 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
801     CALL IMPLDIFF(
802     I bi, bj, iMin, iMax, jMin, jMax,
803     I deltaTtracer, KappaRS, recip_HFacC,
804 adcroft 1.7 U gS,
805 adcroft 1.1 I myThid )
806     ENDIF
807    
808     #ifdef ALLOW_PASSIVE_TRACER
809     IF (tr1Stepping) THEN
810     #ifdef ALLOW_AUTODIFF_TAMC
811 heimbach 1.30 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
812 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
813     CALL IMPLDIFF(
814     I bi, bj, iMin, iMax, jMin, jMax,
815     I deltaTtracer, KappaRT, recip_HFacC,
816 adcroft 1.7 U gTr1,
817 adcroft 1.1 I myThid )
818     ENDIF
819     #endif
820    
821 adcroft 1.17 #ifdef ALLOW_PTRACERS
822     C Vertical diffusion (implicit) for passive tracers
823     IF ( usePTRACERS ) THEN
824     CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
825     ENDIF
826     #endif /* ALLOW_PTRACERS */
827    
828 adcroft 1.1 #ifdef ALLOW_OBCS
829     C-- Apply open boundary conditions
830     IF (useOBCS) THEN
831     DO K=1,Nr
832 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
833 adcroft 1.1 ENDDO
834     END IF
835     #endif /* ALLOW_OBCS */
836    
837     C-- End If implicitDiffusion
838     ENDIF
839 heimbach 1.22
840 jmc 1.39 #ifdef ALLOW_TIMEAVE
841 edhill 1.51 ceh3 needs an IF ( useTIMEAVE ) THEN
842 jmc 1.39 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
843     CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
844     I Nr, deltaTclock, bi, bj, myThid)
845     ENDIF
846     useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
847     IF (taveFreq.GT.0. .AND. useVariableK ) THEN
848     IF (implicitDiffusion) THEN
849     CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
850     I Nr, 3, deltaTclock, bi, bj, myThid)
851     ELSE
852     CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
853     I Nr, 3, deltaTclock, bi, bj, myThid)
854     ENDIF
855     ENDIF
856     #endif /* ALLOW_TIMEAVE */
857    
858 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
859 adcroft 1.1
860 jmc 1.39 C-- end bi,bj loops.
861 adcroft 1.1 ENDDO
862     ENDDO
863 adcroft 1.17
864 edhill 1.56 #ifdef ALLOW_DEBUG
865 adcroft 1.17 If (debugMode) THEN
866     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
867     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
868     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
869     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
870     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
871     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
872     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
873     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
874     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
875 adcroft 1.18 #ifdef ALLOW_PTRACERS
876     IF ( usePTRACERS ) THEN
877     CALL PTRACERS_DEBUG(myThid)
878     ENDIF
879     #endif /* ALLOW_PTRACERS */
880 adcroft 1.17 ENDIF
881 adcroft 1.40 #endif
882    
883 edhill 1.56 #ifdef ALLOW_DEBUG
884 heimbach 1.43 IF ( debugLevel .GE. debLevB )
885     & CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
886 adcroft 1.17 #endif
887 adcroft 1.1
888     RETURN
889     END
  ViewVC Help
Powered by ViewVC 1.1.22