/[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.51 - (hide annotations) (download)
Thu Oct 9 04:19:18 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51i_post
Changes since 1.50: +8 -1 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

Y Y Y Y 13 16 16 16  0 16 16 16 16 16 16 16 16 13 12  0  0 pass  adjustment.128x64x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16  0  0 16 16  0  0 pass  adjustment.cs-32x32x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16 22  0 16 16 22  0 pass  adjust_nlfs.cs-32x32x1
Y Y Y Y -- 13 13 16 16 13 13 13 13 16 16 16 16 16 16 16 16 N/O   advect_cs
Y Y Y Y -- 22 16 16 16 16 16 16 13 16 16 16 16 16 16 16 16 N/O   advect_xy
Y Y Y Y -- 13 16 13 16 16 16 16 16 16 16 22 16 16 16 16 16 N/O   advect_xz
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  aim.5l_cs
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 16 16 16 16 13 16 pass  aim.5l_Equatorial_Channel
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 13 16 16 13 13 16 pass  aim.5l_LatLon
Y Y Y Y 13 16 16 16 16 16 16 16 16 16 13 12 13 13 16 13 16 pass  exp0
Y Y Y Y 14 16 16 16 16 16 16 16 22 16 16 16 13 16 16 22 16 pass  exp1
Y Y Y Y 13 13 16 13 16 16 16 16 16 13 13 16 16 13 13 13 13 pass  exp2
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  exp4
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 22 16 16 16 22 16 pass  exp5
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  front_relax
Y Y Y Y 14 16 16 13 13 16 16 13 13 16 13 13 16 12 13 13 16 pass  global_ocean.90x40x15
Y Y Y Y 10 16 16 13 13 16 13 16 16 13 13 13 13 16 16 13 16 FAIL  global_ocean.cs32x15
Y Y Y Y  6 11 12 13 13 12 13 16 13  9  9  9  9 10  9  9 11 FAIL  global_ocean_pressure
Y Y Y Y 14 16 16 13 16 16 16 13 13 13 13 13 16 12 16 13 16 pass  global_with_exf
Y Y Y Y 14 16 16 16 16 16 16 16 16 11 13 22 13 16 16  9 16 pass  hs94.128x64x5
Y Y Y Y 13 16 16 16 16 16 16 16 16 11 16 16 16 13 16 22 13 pass  hs94.1x64x5
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 13 13 16 16 22 13 pass  hs94.cs-32x32x5
Y Y Y Y 10 10 16 13 13 16 16 16 22 16 13 13 13 13 13 22 13 FAIL  ideal_2D_oce
Y Y Y Y  8 16 16 16 16 16 16 16 16 13 13  8 16 16 16 16 16 FAIL  internal_wave
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 13 22 13 13 13 22 16 pass  inverted_barometer
Y Y Y Y 12 16 16 16 16 16 16 16 16 16 13 12 13 13 13 13 13 FAIL  lab_sea
Y Y Y Y 11 16 16 16 16 16 16 16 13 13 13 12 13 16 13 12 13 FAIL  natl_box
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  plume_on_slope
Y Y Y Y 13 16 16 16 16 13 16 16 16 16 16 16 16 13 16 16 16 pass  solid-body.cs-32x32x1
1 edhill 1.51 C $Header: $
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     _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
150     _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152     _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153     _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154     _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
155     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
156     _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
157     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
158     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
159 cnh 1.9 C This is currently used by IVDC and Diagnostics
160 adcroft 1.1 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
161 jmc 1.39 LOGICAL useVariableK
162 adcroft 1.1 INTEGER iMin, iMax
163     INTEGER jMin, jMax
164     INTEGER bi, bj
165     INTEGER i, j
166     INTEGER k, km1, kup, kDown
167 heimbach 1.42 INTEGER iTracer
168 adcroft 1.1
169 cnh 1.9 CEOP
170 adcroft 1.40
171     #ifndef DISABLE_DEBUGMODE
172 heimbach 1.43 IF ( debugLevel .GE. debLevB )
173     & CALL DEBUG_ENTER('FORWARD_STEP',myThid)
174 adcroft 1.40 #endif
175 adcroft 1.1
176     #ifdef ALLOW_AUTODIFF_TAMC
177     C-- dummy statement to end declaration part
178     ikey = 1
179 heimbach 1.30 itdkey = 1
180 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
181    
182     #ifdef ALLOW_AUTODIFF_TAMC
183     C-- HPF directive to help TAMC
184     CHPF$ INDEPENDENT
185     #endif /* ALLOW_AUTODIFF_TAMC */
186    
187     DO bj=myByLo(myThid),myByHi(myThid)
188    
189     #ifdef ALLOW_AUTODIFF_TAMC
190     C-- HPF directive to help TAMC
191 heimbach 1.2 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
192 jmc 1.37 CHPF$& ,utrans,vtrans,xA,yA
193 heimbach 1.2 CHPF$& ,KappaRT,KappaRS
194 adcroft 1.1 CHPF$& )
195     #endif /* ALLOW_AUTODIFF_TAMC */
196    
197     DO bi=myBxLo(myThid),myBxHi(myThid)
198    
199     #ifdef ALLOW_AUTODIFF_TAMC
200     act1 = bi - myBxLo(myThid)
201     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
202     act2 = bj - myByLo(myThid)
203     max2 = myByHi(myThid) - myByLo(myThid) + 1
204     act3 = myThid - 1
205     max3 = nTx*nTy
206     act4 = ikey_dynamics - 1
207 heimbach 1.30 itdkey = (act1 + 1) + act2*max1
208 adcroft 1.1 & + act3*max1*max2
209     & + act4*max1*max2*max3
210     #endif /* ALLOW_AUTODIFF_TAMC */
211    
212 heimbach 1.41 C-- Set up work arrays with valid (i.e. not NaN) values
213     C These inital values do not alter the numerical results. They
214     C just ensure that all memory references are to valid floating
215     C point numbers. This prevents spurious hardware signals due to
216     C uninitialised but inert locations.
217    
218 adcroft 1.1 DO j=1-OLy,sNy+OLy
219     DO i=1-OLx,sNx+OLx
220 heimbach 1.41 xA(i,j) = 0. _d 0
221     yA(i,j) = 0. _d 0
222     uTrans(i,j) = 0. _d 0
223     vTrans(i,j) = 0. _d 0
224     rhok (i,j) = 0. _d 0
225     rhoKM1 (i,j) = 0. _d 0
226     phiSurfX(i,j) = 0. _d 0
227     phiSurfY(i,j) = 0. _d 0
228 adcroft 1.1 rTrans (i,j) = 0. _d 0
229     fVerT (i,j,1) = 0. _d 0
230     fVerT (i,j,2) = 0. _d 0
231     fVerS (i,j,1) = 0. _d 0
232     fVerS (i,j,2) = 0. _d 0
233     fVerTr1(i,j,1) = 0. _d 0
234     fVerTr1(i,j,2) = 0. _d 0
235     ENDDO
236     ENDDO
237    
238     DO k=1,Nr
239     DO j=1-OLy,sNy+OLy
240     DO i=1-OLx,sNx+OLx
241     C This is currently also used by IVDC and Diagnostics
242 heimbach 1.20 sigmaX(i,j,k) = 0. _d 0
243     sigmaY(i,j,k) = 0. _d 0
244     sigmaR(i,j,k) = 0. _d 0
245 adcroft 1.1 ConvectCount(i,j,k) = 0.
246 heimbach 1.30 KappaRT(i,j,k) = 0. _d 0
247     KappaRS(i,j,k) = 0. _d 0
248 jmc 1.45 C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
249 heimbach 1.30 gT(i,j,k,bi,bj) = 0. _d 0
250     gS(i,j,k,bi,bj) = 0. _d 0
251     # ifdef ALLOW_PASSIVE_TRACER
252 edhill 1.51 ceh3 needs an IF ( use PASSIVE_TRACER) THEN
253 heimbach 1.5 gTr1(i,j,k,bi,bj) = 0. _d 0
254 heimbach 1.30 # endif
255 heimbach 1.42 # ifdef ALLOW_PTRACERS
256 edhill 1.51 ceh3 this should have an IF ( usePTRACERS ) THEN
257 heimbach 1.42 DO iTracer=1,PTRACERS_numInUse
258     gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
259     ENDDO
260     # endif
261 jmc 1.45 #ifdef ALLOW_AUTODIFF_TAMC
262     cph all the following init. are necessary for TAF
263     cph although some of these are re-initialised later.
264 heimbach 1.30 # ifdef ALLOW_GMREDI
265     Kwx(i,j,k,bi,bj) = 0. _d 0
266     Kwy(i,j,k,bi,bj) = 0. _d 0
267     Kwz(i,j,k,bi,bj) = 0. _d 0
268     # ifdef GM_NON_UNITY_DIAGONAL
269     Kux(i,j,k,bi,bj) = 0. _d 0
270     Kvy(i,j,k,bi,bj) = 0. _d 0
271     # endif
272     # ifdef GM_EXTRA_DIAGONAL
273     Kuz(i,j,k,bi,bj) = 0. _d 0
274     Kvz(i,j,k,bi,bj) = 0. _d 0
275     # endif
276     # ifdef GM_BOLUS_ADVEC
277     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
278     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
279 heimbach 1.36 # endif
280     # ifdef GM_VISBECK_VARIABLE_K
281     VisbeckK(i,j,bi,bj) = 0. _d 0
282 heimbach 1.30 # endif
283     # endif /* ALLOW_GMREDI */
284     #endif /* ALLOW_AUTODIFF_TAMC */
285 adcroft 1.1 ENDDO
286     ENDDO
287     ENDDO
288    
289 jmc 1.21 iMin = 1-OLx
290 adcroft 1.1 iMax = sNx+OLx
291 jmc 1.21 jMin = 1-OLy
292 adcroft 1.1 jMax = sNy+OLy
293    
294     #ifdef ALLOW_AUTODIFF_TAMC
295 heimbach 1.30 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
296     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
297 heimbach 1.38 CADJ STORE totphihyd
298     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
299 heimbach 1.11 #ifdef ALLOW_KPP
300 heimbach 1.30 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
301     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
302 heimbach 1.11 #endif
303 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
304    
305 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
306 heimbach 1.43 IF ( debugLevel .GE. debLevB )
307     & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
308 adcroft 1.40 #endif
309    
310 adcroft 1.1 C-- Start of diagnostic loop
311     DO k=Nr,1,-1
312    
313     #ifdef ALLOW_AUTODIFF_TAMC
314     C? Patrick, is this formula correct now that we change the loop range?
315     C? Do we still need this?
316     cph kkey formula corrected.
317     cph Needed for rhok, rhokm1, in the case useGMREDI.
318 heimbach 1.30 kkey = (itdkey-1)*Nr + k
319 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
320    
321     C-- Integrate continuity vertically for vertical velocity
322 jmc 1.27 c CALL INTEGRATE_FOR_W(
323     c I bi, bj, k, uVel, vVel,
324     c O wVel,
325     c I myThid )
326 adcroft 1.1
327     #ifdef ALLOW_OBCS
328     #ifdef ALLOW_NONHYDROSTATIC
329     C-- Apply OBC to W if in N-H mode
330 jmc 1.27 c IF (useOBCS.AND.nonHydrostatic) THEN
331     c CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
332     c ENDIF
333 adcroft 1.1 #endif /* ALLOW_NONHYDROSTATIC */
334     #endif /* ALLOW_OBCS */
335    
336 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
337     C-- MOST of THERMODYNAMICS will be disabled
338     #ifndef SINGLE_LAYER_MODE
339    
340 adcroft 1.1 C-- Calculate gradients of potential density for isoneutral
341     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
342     c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
343     IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
344 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
345 heimbach 1.43 IF ( debugLevel .GE. debLevB )
346     & CALL DEBUG_CALL('FIND_RHO',myThid)
347 adcroft 1.40 #endif
348 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
349     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
350     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
351     #endif /* ALLOW_AUTODIFF_TAMC */
352     CALL FIND_RHO(
353 mlosch 1.26 I bi, bj, iMin, iMax, jMin, jMax, k, k,
354 adcroft 1.1 I theta, salt,
355     O rhoK,
356     I myThid )
357 heimbach 1.30
358 adcroft 1.1 IF (k.GT.1) THEN
359     #ifdef ALLOW_AUTODIFF_TAMC
360     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
361     CADJ STORE salt (:,:,k-1,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-1, k,
365 adcroft 1.1 I theta, salt,
366     O rhoKm1,
367     I myThid )
368     ENDIF
369 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
370 heimbach 1.43 IF ( debugLevel .GE. debLevB )
371     & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
372 adcroft 1.40 #endif
373 adcroft 1.1 CALL GRAD_SIGMA(
374     I bi, bj, iMin, iMax, jMin, jMax, k,
375     I rhoK, rhoKm1, rhoK,
376     O sigmaX, sigmaY, sigmaR,
377     I myThid )
378     ENDIF
379    
380 heimbach 1.30 #ifdef ALLOW_AUTODIFF_TAMC
381     CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
382     CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
383     #endif /* ALLOW_AUTODIFF_TAMC */
384 adcroft 1.1 C-- Implicit Vertical Diffusion for Convection
385     c ==> should use sigmaR !!!
386     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
387 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
388 heimbach 1.43 IF ( debugLevel .GE. debLevB )
389     & CALL DEBUG_CALL('CALC_IVDC',myThid)
390 adcroft 1.40 #endif
391 adcroft 1.1 CALL CALC_IVDC(
392     I bi, bj, iMin, iMax, jMin, jMax, k,
393     I rhoKm1, rhoK,
394     U ConvectCount, KappaRT, KappaRS,
395     I myTime, myIter, myThid)
396     ENDIF
397    
398 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
399    
400 adcroft 1.1 C-- end of diagnostic k loop (Nr:1)
401     ENDDO
402    
403     #ifdef ALLOW_AUTODIFF_TAMC
404     cph avoids recomputation of integrate_for_w
405 heimbach 1.30 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
406 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
407    
408     #ifdef ALLOW_OBCS
409     C-- Calculate future values on open boundaries
410     IF (useOBCS) THEN
411 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
412 heimbach 1.43 IF ( debugLevel .GE. debLevB )
413     & CALL DEBUG_CALL('OBCS_CALC',myThid)
414 adcroft 1.40 #endif
415 jmc 1.16 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
416 adcroft 1.1 I uVel, vVel, wVel, theta, salt,
417     I myThid )
418     ENDIF
419     #endif /* ALLOW_OBCS */
420    
421 cheisey 1.31
422 cheisey 1.32 #ifdef ALLOW_THERM_SEAICE
423 jmc 1.44 IF (useThermSeaIce) THEN
424 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
425 heimbach 1.43 IF ( debugLevel .GE. debLevB )
426     & CALL DEBUG_CALL('ICE_FORCING',myThid)
427 adcroft 1.40 #endif
428 cheisey 1.31 C-- Determines forcing terms based on external fields
429 jmc 1.44 C including effects from ice
430 cheisey 1.31 CALL ICE_FORCING(
431     I bi, bj, iMin, iMax, jMin, jMax,
432     I myThid )
433 jmc 1.50 ELSE
434 jmc 1.44 #endif /* ALLOW_THERM_SEAICE */
435 cheisey 1.31
436 adcroft 1.1 C-- Determines forcing terms based on external fields
437     C relaxation terms, etc.
438 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
439 heimbach 1.43 IF ( debugLevel .GE. debLevB )
440     & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
441 adcroft 1.40 #endif
442 cheisey 1.31 CALL EXTERNAL_FORCING_SURF(
443 adcroft 1.1 I bi, bj, iMin, iMax, jMin, jMax,
444 dimitri 1.47 I myTime, myIter, myThid )
445 jmc 1.50
446     #ifdef ALLOW_THERM_SEAICE
447     C-- end of if/else block useThermSeaIce --
448     ENDIF
449     #endif /* ALLOW_THERM_SEAICE */
450 cheisey 1.31
451 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
452     cph needed for KPP
453     CADJ STORE surfacetendencyU(:,:,bi,bj)
454 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
455 adcroft 1.1 CADJ STORE surfacetendencyV(:,:,bi,bj)
456 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
457 adcroft 1.1 CADJ STORE surfacetendencyS(:,:,bi,bj)
458 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
459 adcroft 1.1 CADJ STORE surfacetendencyT(:,:,bi,bj)
460 heimbach 1.30 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
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 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
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 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
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 KappaRT(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
533     CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
534     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
535     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
536     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
537     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
538 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
539 heimbach 1.30 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
540 adcroft 1.1 #endif
541 heimbach 1.42 #ifdef ALLOW_PTRACERS
542     cph-- moved to forward_step to avoid key computation
543     cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
544     cphCADJ & key=itdkey, byte=isbyte
545     #endif
546 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
547    
548     #ifdef ALLOW_AIM
549     C AIM - atmospheric intermediate model, physics package code.
550     IF ( useAIM ) THEN
551 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
552 heimbach 1.43 IF ( debugLevel .GE. debLevB )
553     & CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
554 adcroft 1.40 #endif
555 jmc 1.33 CALL TIMER_START('AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
556     CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
557     CALL TIMER_STOP( 'AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
558 adcroft 1.1 ENDIF
559     #endif /* ALLOW_AIM */
560 jmc 1.14
561 adcroft 1.12 #ifndef DISABLE_MULTIDIM_ADVECTION
562 adcroft 1.4 C-- Some advection schemes are better calculated using a multi-dimensional
563     C method in the absence of any other terms and, if used, is done here.
564 adcroft 1.13 C
565     C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
566     C The default is to use multi-dimensinal advection for non-linear advection
567     C schemes. However, for the sake of efficiency of the adjoint it is necessary
568     C to be able to exclude this scheme to avoid excessive storage and
569     C recomputation. It *is* differentiable, if you need it.
570     C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
571     C disable this section of code.
572 jmc 1.24 IF (tempMultiDimAdvec) THEN
573 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
574 heimbach 1.43 IF ( debugLevel .GE. debLevB )
575     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
576 adcroft 1.40 #endif
577 heimbach 1.11 CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
578     U theta,gT,
579 adcroft 1.6 I myTime,myIter,myThid)
580 jmc 1.23 ENDIF
581 jmc 1.24 IF (saltMultiDimAdvec) THEN
582 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
583 heimbach 1.43 IF ( debugLevel .GE. debLevB )
584     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
585 adcroft 1.40 #endif
586 heimbach 1.11 CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
587     U salt,gS,
588 adcroft 1.6 I myTime,myIter,myThid)
589     ENDIF
590 adcroft 1.17 C Since passive tracers are configurable separately from T,S we
591     C call the multi-dimensional method for PTRACERS regardless
592     C of whether multiDimAdvection is set or not.
593     #ifdef ALLOW_PTRACERS
594     IF ( usePTRACERS ) THEN
595 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
596 heimbach 1.43 IF ( debugLevel .GE. debLevB )
597     & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
598 adcroft 1.40 #endif
599 adcroft 1.17 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
600     ENDIF
601     #endif /* ALLOW_PTRACERS */
602 adcroft 1.12 #endif /* DISABLE_MULTIDIM_ADVECTION */
603 adcroft 1.1
604 adcroft 1.40 #ifndef DISABLE_DEBUGMODE
605 heimbach 1.43 IF ( debugLevel .GE. debLevB )
606     & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
607 adcroft 1.40 #endif
608    
609 adcroft 1.1 C-- Start of thermodynamics loop
610     DO k=Nr,1,-1
611     #ifdef ALLOW_AUTODIFF_TAMC
612     C? Patrick Is this formula correct?
613     cph Yes, but I rewrote it.
614     cph Also, the KappaR? need the index and subscript k!
615 heimbach 1.30 kkey = (itdkey-1)*Nr + k
616 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
617    
618     C-- km1 Points to level above k (=k-1)
619     C-- kup Cycles through 1,2 to point to layer above
620     C-- kDown Cycles through 2,1 to point to current layer
621    
622     km1 = MAX(1,k-1)
623     kup = 1+MOD(k+1,2)
624     kDown= 1+MOD(k,2)
625    
626     iMin = 1-OLx
627     iMax = sNx+OLx
628     jMin = 1-OLy
629     jMax = sNy+OLy
630    
631     C-- Get temporary terms used by tendency routines
632     CALL CALC_COMMON_FACTORS (
633     I bi,bj,iMin,iMax,jMin,jMax,k,
634     O xA,yA,uTrans,vTrans,rTrans,maskUp,
635     I myThid)
636 jmc 1.19
637     #ifdef ALLOW_GMREDI
638 heimbach 1.35
639 jmc 1.19 C-- Residual transp = Bolus transp + Eulerian transp
640     IF (useGMRedi) THEN
641     CALL GMREDI_CALC_UVFLOW(
642     & uTrans, vTrans, bi, bj, k, myThid)
643     IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
644     & rTrans, bi, bj, k, myThid)
645     ENDIF
646 heimbach 1.35
647     #ifdef ALLOW_AUTODIFF_TAMC
648     #ifdef GM_BOLUS_ADVEC
649     CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
650     CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
651     CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
652     #endif
653     #endif /* ALLOW_AUTODIFF_TAMC */
654    
655 jmc 1.19 #endif /* ALLOW_GMREDI */
656 adcroft 1.1
657     #ifdef ALLOW_AUTODIFF_TAMC
658     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
659     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
660     #endif /* ALLOW_AUTODIFF_TAMC */
661    
662     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
663     C-- Calculate the total vertical diffusivity
664     CALL CALC_DIFFUSIVITY(
665     I bi,bj,iMin,iMax,jMin,jMax,k,
666     I maskUp,
667 heimbach 1.2 O KappaRT,KappaRS,
668 adcroft 1.1 I myThid)
669     #endif
670    
671     iMin = 1-OLx+2
672     iMax = sNx+OLx-1
673     jMin = 1-OLy+2
674     jMax = sNy+OLy-1
675    
676     C-- Calculate active tracer tendencies (gT,gS,...)
677     C and step forward storing result in gTnm1, gSnm1, etc.
678     IF ( tempStepping ) THEN
679     CALL CALC_GT(
680     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
681     I xA,yA,uTrans,vTrans,rTrans,maskUp,
682     I KappaRT,
683     U fVerT,
684 adcroft 1.7 I myTime,myIter,myThid)
685 adcroft 1.1 CALL TIMESTEP_TRACER(
686 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
687 adcroft 1.1 I theta, gT,
688     I myIter, myThid)
689     ENDIF
690 jmc 1.44
691 cheisey 1.32 #ifdef ALLOW_THERM_SEAICE
692 jmc 1.44 IF (useThermSeaIce .AND. k.EQ.1) THEN
693     CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )
694     ENDIF
695 cheisey 1.32 #endif
696 jmc 1.44
697 adcroft 1.1 IF ( saltStepping ) THEN
698     CALL CALC_GS(
699     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
700     I xA,yA,uTrans,vTrans,rTrans,maskUp,
701     I KappaRS,
702     U fVerS,
703 adcroft 1.7 I myTime,myIter,myThid)
704 adcroft 1.1 CALL TIMESTEP_TRACER(
705 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
706 adcroft 1.1 I salt, gS,
707     I myIter, myThid)
708     ENDIF
709     #ifdef ALLOW_PASSIVE_TRACER
710 edhill 1.51 ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
711 adcroft 1.1 IF ( tr1Stepping ) THEN
712     CALL CALC_GTR1(
713     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
714     I xA,yA,uTrans,vTrans,rTrans,maskUp,
715     I KappaRT,
716     U fVerTr1,
717 heimbach 1.8 I myTime,myIter,myThid)
718 adcroft 1.1 CALL TIMESTEP_TRACER(
719 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
720 adcroft 1.1 I Tr1, gTr1,
721 heimbach 1.8 I myIter,myThid)
722 adcroft 1.1 ENDIF
723     #endif
724 adcroft 1.17 #ifdef ALLOW_PTRACERS
725     IF ( usePTRACERS ) THEN
726 heimbach 1.42 CALL PTRACERS_INTEGRATE(
727 adcroft 1.17 I bi,bj,k,
728     I xA,yA,uTrans,vTrans,rTrans,maskUp,
729     X KappaRS,
730     I myIter,myTime,myThid)
731     ENDIF
732     #endif /* ALLOW_PTRACERS */
733 adcroft 1.1
734     #ifdef ALLOW_OBCS
735     C-- Apply open boundary conditions
736     IF (useOBCS) THEN
737 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
738 adcroft 1.1 END IF
739     #endif /* ALLOW_OBCS */
740    
741     C-- Freeze water
742 jmc 1.44 IF ( allowFreezing .AND. .NOT. useSEAICE
743     & .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN
744 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
745 heimbach 1.11 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
746 adcroft 1.1 CADJ & , key = kkey, byte = isbyte
747     #endif /* ALLOW_AUTODIFF_TAMC */
748     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
749 jmc 1.44 ENDIF
750 adcroft 1.1
751     C-- end of thermodynamic k loop (Nr:1)
752     ENDDO
753 cheisey 1.31
754     cswdice -- add ---
755 cheisey 1.32 #ifdef ALLOW_THERM_SEAICE
756 cheisey 1.31 c timeaveraging for ice model values
757 edhill 1.51 ceh3 This should be wrapped in an IF ( useThermSeaIce ) THEN
758 cheisey 1.31 CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
759     #endif
760     cswdice --- end add ---
761    
762    
763    
764 adcroft 1.1
765     C-- Implicit diffusion
766     IF (implicitDiffusion) THEN
767    
768     IF (tempStepping) THEN
769     #ifdef ALLOW_AUTODIFF_TAMC
770 heimbach 1.30 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
771 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
772     CALL IMPLDIFF(
773     I bi, bj, iMin, iMax, jMin, jMax,
774     I deltaTtracer, KappaRT, recip_HFacC,
775 adcroft 1.7 U gT,
776 adcroft 1.1 I myThid )
777     ENDIF
778    
779     IF (saltStepping) THEN
780     #ifdef ALLOW_AUTODIFF_TAMC
781 heimbach 1.30 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
782 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
783     CALL IMPLDIFF(
784     I bi, bj, iMin, iMax, jMin, jMax,
785     I deltaTtracer, KappaRS, recip_HFacC,
786 adcroft 1.7 U gS,
787 adcroft 1.1 I myThid )
788     ENDIF
789    
790     #ifdef ALLOW_PASSIVE_TRACER
791     IF (tr1Stepping) THEN
792     #ifdef ALLOW_AUTODIFF_TAMC
793 heimbach 1.30 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
794 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
795     CALL IMPLDIFF(
796     I bi, bj, iMin, iMax, jMin, jMax,
797     I deltaTtracer, KappaRT, recip_HFacC,
798 adcroft 1.7 U gTr1,
799 adcroft 1.1 I myThid )
800     ENDIF
801     #endif
802    
803 adcroft 1.17 #ifdef ALLOW_PTRACERS
804     C Vertical diffusion (implicit) for passive tracers
805     IF ( usePTRACERS ) THEN
806     CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
807     ENDIF
808     #endif /* ALLOW_PTRACERS */
809    
810 adcroft 1.1 #ifdef ALLOW_OBCS
811     C-- Apply open boundary conditions
812     IF (useOBCS) THEN
813     DO K=1,Nr
814 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
815 adcroft 1.1 ENDDO
816     END IF
817     #endif /* ALLOW_OBCS */
818    
819     C-- End If implicitDiffusion
820     ENDIF
821 heimbach 1.22
822 jmc 1.39 #ifdef ALLOW_TIMEAVE
823 edhill 1.51 ceh3 needs an IF ( useTIMEAVE ) THEN
824 jmc 1.39 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
825     CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
826     I Nr, deltaTclock, bi, bj, myThid)
827     ENDIF
828     useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
829     IF (taveFreq.GT.0. .AND. useVariableK ) THEN
830     IF (implicitDiffusion) THEN
831     CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
832     I Nr, 3, deltaTclock, bi, bj, myThid)
833     ELSE
834     CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
835     I Nr, 3, deltaTclock, bi, bj, myThid)
836     ENDIF
837     ENDIF
838     #endif /* ALLOW_TIMEAVE */
839    
840 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
841 adcroft 1.1
842 jmc 1.39 C-- end bi,bj loops.
843 adcroft 1.1 ENDDO
844     ENDDO
845 adcroft 1.17
846     #ifndef DISABLE_DEBUGMODE
847     If (debugMode) THEN
848     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
849     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
850     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
851     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
852     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
853     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
854     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
855     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
856     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
857 adcroft 1.18 #ifdef ALLOW_PTRACERS
858     IF ( usePTRACERS ) THEN
859     CALL PTRACERS_DEBUG(myThid)
860     ENDIF
861     #endif /* ALLOW_PTRACERS */
862 adcroft 1.17 ENDIF
863 adcroft 1.40 #endif
864    
865     #ifndef DISABLE_DEBUGMODE
866 heimbach 1.43 IF ( debugLevel .GE. debLevB )
867     & CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
868 adcroft 1.17 #endif
869 adcroft 1.1
870     RETURN
871     END
  ViewVC Help
Powered by ViewVC 1.1.22