/[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.103 - (hide annotations) (download)
Thu Apr 20 03:29:35 2006 UTC (18 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58f_post, checkpoint58d_post, checkpoint58e_post, checkpoint58g_post
Changes since 1.102: +4 -2 lines
Bracket Kuz, Kvz stores

1 heimbach 1.103 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.102 2006/04/04 14:52:43 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 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 jmc 1.21 #endif /* ALLOW_AUTODIFF_TAMC */
15 adcroft 1.1
16 cnh 1.9 CBOP
17     C !ROUTINE: THERMODYNAMICS
18     C !INTERFACE:
19 adcroft 1.1 SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
20 cnh 1.9 C !DESCRIPTION: \bv
21     C *==========================================================*
22 jmc 1.101 C | SUBROUTINE THERMODYNAMICS
23     C | o Controlling routine for the prognostic part of the
24     C | thermo-dynamics.
25 cnh 1.9 C *===========================================================
26     C | The algorithm...
27     C |
28     C | "Correction Step"
29     C | =================
30     C | Here we update the horizontal velocities with the surface
31     C | pressure such that the resulting flow is either consistent
32     C | with the free-surface evolution or the rigid-lid:
33     C | U[n] = U* + dt x d/dx P
34     C | V[n] = V* + dt x d/dy P
35     C |
36     C | "Calculation of Gs"
37     C | ===================
38     C | This is where all the accelerations and tendencies (ie.
39     C | physics, parameterizations etc...) are calculated
40     C | rho = rho ( theta[n], salt[n] )
41     C | b = b(rho, theta)
42     C | K31 = K31 ( rho )
43     C | Gu[n] = Gu( u[n], v[n], wVel, b, ... )
44     C | Gv[n] = Gv( u[n], v[n], wVel, b, ... )
45     C | Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
46     C | Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
47     C |
48     C | "Time-stepping" or "Prediction"
49     C | ================================
50     C | The models variables are stepped forward with the appropriate
51     C | time-stepping scheme (currently we use Adams-Bashforth II)
52     C | - For momentum, the result is always *only* a "prediction"
53     C | in that the flow may be divergent and will be "corrected"
54     C | later with a surface pressure gradient.
55     C | - Normally for tracers the result is the new field at time
56     C | level [n+1} *BUT* in the case of implicit diffusion the result
57     C | is also *only* a prediction.
58     C | - We denote "predictors" with an asterisk (*).
59     C | U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
60     C | V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
61     C | theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
62     C | salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63     C | With implicit diffusion:
64     C | theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
65     C | salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
66     C | (1 + dt * K * d_zz) theta[n] = theta*
67     C | (1 + dt * K * d_zz) salt[n] = salt*
68     C |
69     C *==========================================================*
70     C \ev
71    
72     C !USES:
73 adcroft 1.1 IMPLICIT NONE
74     C == Global variables ===
75     #include "SIZE.h"
76     #include "EEPARAMS.h"
77     #include "PARAMS.h"
78     #include "DYNVARS.h"
79     #include "GRID.h"
80 jmc 1.83 #ifdef ALLOW_GENERIC_ADVDIFF
81 adcroft 1.4 #include "GAD.h"
82 jmc 1.83 #endif
83 jmc 1.45 #ifdef ALLOW_PTRACERS
84 jmc 1.74 #include "PTRACERS_SIZE.h"
85 jmc 1.45 #include "PTRACERS.h"
86     #endif
87 heimbach 1.42 #ifdef ALLOW_TIMEAVE
88     #include "TIMEAVE_STATV.h"
89     #endif
90    
91 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
92     # include "tamc.h"
93     # include "tamc_keys.h"
94     # include "FFIELDS.h"
95 heimbach 1.30 # include "EOS.h"
96 adcroft 1.1 # ifdef ALLOW_KPP
97     # include "KPP.h"
98 heimbach 1.67 # endif
99     # ifdef ALLOW_GMREDI
100     # include "GMREDI.h"
101 adcroft 1.1 # endif
102 heimbach 1.68 # ifdef ALLOW_EBM
103     # include "EBM.h"
104 adcroft 1.1 # 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 jmc 1.83 #ifdef ALLOW_GENERIC_ADVDIFF
117 cnh 1.9 C !LOCAL VARIABLES:
118 adcroft 1.1 C == Local variables
119     C xA, yA - Per block temporaries holding face areas
120     C uTrans, vTrans, rTrans - Per block temporaries holding flow
121     C transport
122     C o uTrans: Zonal transport
123     C o vTrans: Meridional transport
124     C o rTrans: Vertical transport
125 jmc 1.64 C rTransKp1 o vertical volume transp. at interface k+1
126 adcroft 1.1 C maskUp o maskUp: land/water mask for W points
127     C fVer[STUV] o fVer: Vertical flux term - note fVer
128     C is "pipelined" in the vertical
129     C so we need an fVer for each
130     C variable.
131 jmc 1.81 C kappaRT, - Total diffusion in vertical at level k, for T and S
132     C kappaRS (background + spatially varying, isopycnal term).
133     C kappaRTr - Total diffusion in vertical at level k,
134     C for each passive Tracer
135     C kappaRk - Total diffusion in vertical, all levels, 1 tracer
136 jmc 1.39 C useVariableK = T when vertical diffusion is not constant
137 adcroft 1.1 C iMin, iMax - Ranges and sub-block indices on which calculations
138     C jMin, jMax are applied.
139     C bi, bj
140     C k, kup, - Index for layer above and below. kup and kDown
141     C kDown, km1 are switched with layer to be the appropriate
142     C index into fVerTerm.
143     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148 jmc 1.64 _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149 adcroft 1.1 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
152 jmc 1.81 _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
153     _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
154 heimbach 1.55 #ifdef ALLOW_PTRACERS
155     _RL fVerP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
156 jmc 1.81 _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
157 heimbach 1.55 #endif
158 jmc 1.81 _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
159 adcroft 1.1 INTEGER iMin, iMax
160     INTEGER jMin, jMax
161     INTEGER bi, bj
162     INTEGER i, j
163     INTEGER k, km1, kup, kDown
164 jmc 1.95 #ifdef ALLOW_ADAMSBASHFORTH_3
165     INTEGER iterNb, m1, m2
166     #endif
167     #ifdef ALLOW_TIMEAVE
168     LOGICAL useVariableK
169     #endif
170     #ifdef ALLOW_PTRACERS
171 heimbach 1.55 INTEGER iTracer, ip
172 jmc 1.95 #endif
173 adcroft 1.1
174 cnh 1.9 CEOP
175 adcroft 1.40
176 edhill 1.56 #ifdef ALLOW_DEBUG
177 heimbach 1.43 IF ( debugLevel .GE. debLevB )
178 jmc 1.63 & CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
179 adcroft 1.40 #endif
180 jmc 1.101
181 adcroft 1.1 #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 jmc 1.81 CHPF$& ,kappaRT,kappaRS
199 adcroft 1.1 CHPF$& )
200 heimbach 1.102 # ifdef ALLOW_PTRACERS
201     CHPF$ INDEPENDENT, NEW (fVerP,kappaRTr)
202     # endif
203 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
204    
205     DO bi=myBxLo(myThid),myBxHi(myThid)
206    
207     #ifdef ALLOW_AUTODIFF_TAMC
208     act1 = bi - myBxLo(myThid)
209     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
210     act2 = bj - myByLo(myThid)
211     max2 = myByHi(myThid) - myByLo(myThid) + 1
212     act3 = myThid - 1
213     max3 = nTx*nTy
214     act4 = ikey_dynamics - 1
215 heimbach 1.30 itdkey = (act1 + 1) + act2*max1
216 adcroft 1.1 & + act3*max1*max2
217     & + act4*max1*max2*max3
218     #endif /* ALLOW_AUTODIFF_TAMC */
219    
220 heimbach 1.41 C-- Set up work arrays with valid (i.e. not NaN) values
221     C These inital values do not alter the numerical results. They
222     C just ensure that all memory references are to valid floating
223     C point numbers. This prevents spurious hardware signals due to
224     C uninitialised but inert locations.
225    
226 adcroft 1.1 DO j=1-OLy,sNy+OLy
227     DO i=1-OLx,sNx+OLx
228 heimbach 1.41 xA(i,j) = 0. _d 0
229     yA(i,j) = 0. _d 0
230     uTrans(i,j) = 0. _d 0
231     vTrans(i,j) = 0. _d 0
232 adcroft 1.1 rTrans (i,j) = 0. _d 0
233 jmc 1.64 rTransKp1(i,j) = 0. _d 0
234 adcroft 1.1 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 jmc 1.81 kappaRT(i,j) = 0. _d 0
239     kappaRS(i,j) = 0. _d 0
240 adcroft 1.1 ENDDO
241     ENDDO
242    
243     DO k=1,Nr
244     DO j=1-OLy,sNy+OLy
245     DO i=1-OLx,sNx+OLx
246     C This is currently also used by IVDC and Diagnostics
247 jmc 1.81 kappaRk(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 jmc 1.97 ENDDO
252     ENDDO
253     ENDDO
254    
255     #ifdef ALLOW_PTRACERS
256     IF ( usePTRACERS ) THEN
257     DO ip=1,PTRACERS_num
258     DO j=1-OLy,sNy+OLy
259     DO i=1-OLx,sNx+OLx
260     fVerP (i,j,1,ip) = 0. _d 0
261     fVerP (i,j,2,ip) = 0. _d 0
262     kappaRTr(i,j,ip) = 0. _d 0
263     ENDDO
264     ENDDO
265     ENDDO
266     C- set tracer tendency to zero:
267 heimbach 1.102 DO iTracer=1,PTRACERS_num
268 jmc 1.97 DO k=1,Nr
269     DO j=1-OLy,sNy+OLy
270     DO i=1-OLx,sNx+OLx
271     gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
272     ENDDO
273 heimbach 1.42 ENDDO
274 adcroft 1.1 ENDDO
275     ENDDO
276 jmc 1.97 ENDIF
277     #endif
278 adcroft 1.1
279 jmc 1.95 #ifdef ALLOW_ADAMSBASHFORTH_3
280     C- Apply AB on T,S :
281     iterNb = myIter
282     IF (staggerTimeStep) iterNb = myIter - 1
283     m1 = 1 + MOD(iterNb+1,2)
284     m2 = 1 + MOD( iterNb ,2)
285     C compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
286     IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
287     I bi, bj, 0,
288     U theta, gtNm,
289     I tempStartAB, iterNb, myThid )
290     C compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
291     IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
292     I bi, bj, 0,
293     U salt, gsNm,
294     I saltStartAB, iterNb, myThid )
295     #endif /* ALLOW_ADAMSBASHFORTH_3 */
296    
297 jmc 1.72 c iMin = 1-OLx
298     c iMax = sNx+OLx
299     c jMin = 1-OLy
300     c jMax = sNy+OLy
301 adcroft 1.1
302     #ifdef ALLOW_AUTODIFF_TAMC
303     cph avoids recomputation of integrate_for_w
304 heimbach 1.30 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
305 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
306    
307 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
308     C-- MOST of THERMODYNAMICS will be disabled
309     #ifndef SINGLE_LAYER_MODE
310    
311 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
312 heimbach 1.30 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
313     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
314     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
315     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
316 heimbach 1.99 cph-test
317     cphCADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318     cphCADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
319 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
320    
321 adcroft 1.12 #ifndef DISABLE_MULTIDIM_ADVECTION
322 adcroft 1.4 C-- Some advection schemes are better calculated using a multi-dimensional
323     C method in the absence of any other terms and, if used, is done here.
324 adcroft 1.13 C
325     C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
326     C The default is to use multi-dimensinal advection for non-linear advection
327     C schemes. However, for the sake of efficiency of the adjoint it is necessary
328     C to be able to exclude this scheme to avoid excessive storage and
329     C recomputation. It *is* differentiable, if you need it.
330     C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
331     C disable this section of code.
332 jmc 1.24 IF (tempMultiDimAdvec) THEN
333 edhill 1.56 #ifdef ALLOW_DEBUG
334 heimbach 1.43 IF ( debugLevel .GE. debLevB )
335     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
336 adcroft 1.40 #endif
337 jmc 1.63 CALL GAD_ADVECTION(
338 jmc 1.69 I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
339     I GAD_TEMPERATURE,
340 jmc 1.63 I uVel, vVel, wVel, theta,
341     O gT,
342     I bi,bj,myTime,myIter,myThid)
343 jmc 1.23 ENDIF
344 jmc 1.24 IF (saltMultiDimAdvec) THEN
345 edhill 1.56 #ifdef ALLOW_DEBUG
346 heimbach 1.43 IF ( debugLevel .GE. debLevB )
347     & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
348 adcroft 1.40 #endif
349 jmc 1.63 CALL GAD_ADVECTION(
350 jmc 1.69 I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
351     I GAD_SALINITY,
352 jmc 1.63 I uVel, vVel, wVel, salt,
353     O gS,
354     I bi,bj,myTime,myIter,myThid)
355 adcroft 1.6 ENDIF
356 jmc 1.101
357 adcroft 1.17 C Since passive tracers are configurable separately from T,S we
358     C call the multi-dimensional method for PTRACERS regardless
359     C of whether multiDimAdvection is set or not.
360     #ifdef ALLOW_PTRACERS
361     IF ( usePTRACERS ) THEN
362 edhill 1.56 #ifdef ALLOW_DEBUG
363 heimbach 1.43 IF ( debugLevel .GE. debLevB )
364     & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
365 adcroft 1.40 #endif
366 adcroft 1.17 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
367     ENDIF
368     #endif /* ALLOW_PTRACERS */
369 adcroft 1.12 #endif /* DISABLE_MULTIDIM_ADVECTION */
370 adcroft 1.1
371 edhill 1.56 #ifdef ALLOW_DEBUG
372 jmc 1.63 IF ( debugLevel .GE. debLevB )
373 heimbach 1.43 & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
374 adcroft 1.40 #endif
375    
376 adcroft 1.1 C-- Start of thermodynamics loop
377     DO k=Nr,1,-1
378     #ifdef ALLOW_AUTODIFF_TAMC
379     C? Patrick Is this formula correct?
380     cph Yes, but I rewrote it.
381 jmc 1.81 cph Also, the kappaR? need the index and subscript k!
382 heimbach 1.30 kkey = (itdkey-1)*Nr + k
383 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
384    
385     C-- km1 Points to level above k (=k-1)
386     C-- kup Cycles through 1,2 to point to layer above
387     C-- kDown Cycles through 2,1 to point to current layer
388    
389     km1 = MAX(1,k-1)
390     kup = 1+MOD(k+1,2)
391     kDown= 1+MOD(k,2)
392    
393     iMin = 1-OLx
394     iMax = sNx+OLx
395     jMin = 1-OLy
396     jMax = sNy+OLy
397    
398 jmc 1.95 IF (k.EQ.Nr) THEN
399 jmc 1.64 DO j=1-Oly,sNy+Oly
400     DO i=1-Olx,sNx+Olx
401 jmc 1.95 rTransKp1(i,j) = 0. _d 0
402 jmc 1.64 ENDDO
403     ENDDO
404 jmc 1.95 ELSE
405     DO j=1-Oly,sNy+Oly
406     DO i=1-Olx,sNx+Olx
407     rTransKp1(i,j) = rTrans(i,j)
408     ENDDO
409     ENDDO
410     ENDIF
411 heimbach 1.66 #ifdef ALLOW_AUTODIFF_TAMC
412     CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
413     #endif
414 jmc 1.64
415 adcroft 1.1 C-- Get temporary terms used by tendency routines
416     CALL CALC_COMMON_FACTORS (
417     I bi,bj,iMin,iMax,jMin,jMax,k,
418     O xA,yA,uTrans,vTrans,rTrans,maskUp,
419     I myThid)
420 jmc 1.19
421 jmc 1.64 IF (k.EQ.1) THEN
422     C- Surface interface :
423     DO j=1-Oly,sNy+Oly
424     DO i=1-Olx,sNx+Olx
425     rTrans(i,j) = 0.
426     ENDDO
427     ENDDO
428     ELSE
429     C- Interior interface :
430     DO j=1-Oly,sNy+Oly
431     DO i=1-Olx,sNx+Olx
432     rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
433     ENDDO
434     ENDDO
435     ENDIF
436    
437 jmc 1.19 #ifdef ALLOW_GMREDI
438 heimbach 1.35
439 jmc 1.19 C-- Residual transp = Bolus transp + Eulerian transp
440     IF (useGMRedi) THEN
441     CALL GMREDI_CALC_UVFLOW(
442     & uTrans, vTrans, bi, bj, k, myThid)
443     IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
444     & rTrans, bi, bj, k, myThid)
445     ENDIF
446 heimbach 1.35
447 heimbach 1.66 #ifdef ALLOW_AUTODIFF_TAMC
448     CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
449 heimbach 1.35 #ifdef GM_BOLUS_ADVEC
450     CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
451     CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
452     #endif
453     #endif /* ALLOW_AUTODIFF_TAMC */
454    
455 jmc 1.19 #endif /* ALLOW_GMREDI */
456 adcroft 1.1
457     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
458     C-- Calculate the total vertical diffusivity
459 jmc 1.81 IF ( .NOT.implicitDiffusion ) THEN
460     CALL CALC_DIFFUSIVITY(
461     I bi,bj,iMin,iMax,jMin,jMax,k,
462     I maskUp,
463     O kappaRT,kappaRS,
464     I myThid)
465     ENDIF
466 heimbach 1.52 # ifdef ALLOW_AUTODIFF_TAMC
467 jmc 1.81 CADJ STORE kappaRT(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
468     CADJ STORE kappaRS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
469 heimbach 1.52 # endif /* ALLOW_AUTODIFF_TAMC */
470 adcroft 1.1 #endif
471    
472     iMin = 1-OLx+2
473     iMax = sNx+OLx-1
474     jMin = 1-OLy+2
475     jMax = sNy+OLy-1
476    
477     C-- Calculate active tracer tendencies (gT,gS,...)
478 jmc 1.90 C and step forward storing result in gT, gS, etc.
479 heimbach 1.100 C--
480     # ifdef ALLOW_AUTODIFF_TAMC
481     # if (defined (NONLIN_FRSURF) && defined (ALLOW_GMREDI))
482     CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
483 heimbach 1.103 CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
484     # ifdef GM_EXTRA_DIAGONAL
485 heimbach 1.100 CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
486     CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
487 heimbach 1.103 # endif
488 heimbach 1.100 # endif
489     # endif /* ALLOW_AUTODIFF_TAMC */
490 adcroft 1.1 IF ( tempStepping ) THEN
491 heimbach 1.99 #ifdef ALLOW_AUTODIFF_TAMC
492     CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
493     # ifdef NONLIN_FRSURF
494 heimbach 1.96 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
495 heimbach 1.99 # endif
496     #endif /* ALLOW_AUTODIFF_TAMC */
497 adcroft 1.1 CALL CALC_GT(
498     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
499 jmc 1.64 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
500 jmc 1.81 I kappaRT,
501 adcroft 1.1 U fVerT,
502 adcroft 1.7 I myTime,myIter,myThid)
503 jmc 1.95 #ifdef ALLOW_ADAMSBASHFORTH_3
504     IF ( AdamsBashforth_T ) THEN
505 adcroft 1.1 CALL TIMESTEP_TRACER(
506 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
507 jmc 1.95 I gtNm(1-Olx,1-Oly,1,1,1,m2),
508     U gT,
509 adcroft 1.1 I myIter, myThid)
510 jmc 1.95 ELSE
511     #endif
512     CALL TIMESTEP_TRACER(
513     I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
514     I theta,
515     U gT,
516     I myIter, myThid)
517     #ifdef ALLOW_ADAMSBASHFORTH_3
518     ENDIF
519     #endif
520 adcroft 1.1 ENDIF
521 jmc 1.44
522 adcroft 1.1 IF ( saltStepping ) THEN
523 heimbach 1.99 #ifdef ALLOW_AUTODIFF_TAMC
524     CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
525     # ifdef NONLIN_FRSURF
526 heimbach 1.96 CADJ STORE gS(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
527 heimbach 1.99 # endif
528     #endif /* ALLOW_AUTODIFF_TAMC */
529 adcroft 1.1 CALL CALC_GS(
530     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
531 jmc 1.64 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
532 jmc 1.81 I kappaRS,
533 adcroft 1.1 U fVerS,
534 adcroft 1.7 I myTime,myIter,myThid)
535 jmc 1.95 #ifdef ALLOW_ADAMSBASHFORTH_3
536     IF ( AdamsBashforth_S ) THEN
537 adcroft 1.1 CALL TIMESTEP_TRACER(
538 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
539 jmc 1.95 I gsNm(1-Olx,1-Oly,1,1,1,m2),
540     U gS,
541 adcroft 1.1 I myIter, myThid)
542 jmc 1.95 ELSE
543     #endif
544     CALL TIMESTEP_TRACER(
545     I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
546     I salt,
547     U gS,
548     I myIter, myThid)
549     #ifdef ALLOW_ADAMSBASHFORTH_3
550     ENDIF
551     #endif
552 adcroft 1.1 ENDIF
553 jmc 1.101
554 adcroft 1.17 #ifdef ALLOW_PTRACERS
555     IF ( usePTRACERS ) THEN
556 jmc 1.81 IF ( .NOT.implicitDiffusion ) THEN
557     CALL PTRACERS_CALC_DIFF(
558     I bi,bj,iMin,iMax,jMin,jMax,k,
559     I maskUp,
560     O kappaRTr,
561     I myThid)
562     ENDIF
563 heimbach 1.89 # ifdef ALLOW_AUTODIFF_TAMC
564     CADJ STORE kappaRTr(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
565     # endif /* ALLOW_AUTODIFF_TAMC */
566 heimbach 1.42 CALL PTRACERS_INTEGRATE(
567 adcroft 1.17 I bi,bj,k,
568 jmc 1.64 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
569 jmc 1.80 U fVerP,
570 jmc 1.81 I kappaRTr,
571 adcroft 1.17 I myIter,myTime,myThid)
572     ENDIF
573     #endif /* ALLOW_PTRACERS */
574 adcroft 1.1
575     #ifdef ALLOW_OBCS
576     C-- Apply open boundary conditions
577     IF (useOBCS) THEN
578 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
579 adcroft 1.1 END IF
580     #endif /* ALLOW_OBCS */
581 edhill 1.54
582 jmc 1.59 C-- Freeze water
583     C this bit of code is left here for backward compatibility.
584     C freezing at surface level has been moved to FORWARD_STEP
585     IF ( useOldFreezing .AND. .NOT. useSEAICE
586 jmc 1.61 & .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
587 jmc 1.59 #ifdef ALLOW_AUTODIFF_TAMC
588     CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
589     CADJ & , key = kkey, byte = isbyte
590     #endif /* ALLOW_AUTODIFF_TAMC */
591     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
592     ENDIF
593 adcroft 1.1
594     C-- end of thermodynamic k loop (Nr:1)
595     ENDDO
596 cheisey 1.31
597 spk 1.91 C All explicit advection/diffusion/sources should now be
598     C done. The updated tracer field is in gPtr. Accumalate
599     C explicit tendency and also reset gPtr to initial tracer
600     C field for implicit matrix calculation
601    
602     #ifdef ALLOW_MATRIX
603     IF (useMATRIX)
604 edhill 1.92 & CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
605 spk 1.91 #endif
606    
607 jmc 1.81 iMin = 1
608     iMax = sNx
609     jMin = 1
610     jMax = sNy
611 adcroft 1.1
612 jmc 1.63 C-- Implicit vertical advection & diffusion
613 jmc 1.81 IF ( tempStepping .AND. implicitDiffusion ) THEN
614     CALL CALC_3D_DIFFUSIVITY(
615     I bi,bj,iMin,iMax,jMin,jMax,
616     I GAD_TEMPERATURE, useGMredi, useKPP,
617     O kappaRk,
618     I myThid)
619     ENDIF
620 jmc 1.63 #ifdef INCLUDE_IMPLVERTADV_CODE
621 jmc 1.87 IF ( tempImplVertAdv ) THEN
622 jmc 1.63 CALL GAD_IMPLICIT_R(
623     I tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
624 jmc 1.81 I kappaRk, wVel, theta,
625 jmc 1.63 U gT,
626     I bi, bj, myTime, myIter, myThid )
627     ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
628     #else /* INCLUDE_IMPLVERTADV_CODE */
629     IF ( tempStepping .AND. implicitDiffusion ) THEN
630     #endif /* INCLUDE_IMPLVERTADV_CODE */
631 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
632 jmc 1.81 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
633 heimbach 1.30 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
634 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
635 jmc 1.63 CALL IMPLDIFF(
636 adcroft 1.1 I bi, bj, iMin, iMax, jMin, jMax,
637 jmc 1.87 I GAD_TEMPERATURE, kappaRk, recip_hFacC,
638 adcroft 1.7 U gT,
639 adcroft 1.1 I myThid )
640 jmc 1.63 ENDIF
641 jmc 1.81
642     #ifdef ALLOW_TIMEAVE
643     useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
644     & .OR. useGMredi .OR. ivdc_kappa.NE.0.
645     IF (taveFreq.GT.0. .AND. useVariableK ) THEN
646     IF (implicitDiffusion) THEN
647     CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
648     I Nr, 3, deltaTclock, bi, bj, myThid)
649     c ELSE
650     c CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
651     c I Nr, 3, deltaTclock, bi, bj, myThid)
652     ENDIF
653     ENDIF
654     #endif /* ALLOW_TIMEAVE */
655 adcroft 1.1
656 jmc 1.81 IF ( saltStepping .AND. implicitDiffusion ) THEN
657     CALL CALC_3D_DIFFUSIVITY(
658     I bi,bj,iMin,iMax,jMin,jMax,
659     I GAD_SALINITY, useGMredi, useKPP,
660     O kappaRk,
661     I myThid)
662     ENDIF
663    
664 jmc 1.63 #ifdef INCLUDE_IMPLVERTADV_CODE
665 jmc 1.87 IF ( saltImplVertAdv ) THEN
666 jmc 1.63 CALL GAD_IMPLICIT_R(
667     I saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
668 jmc 1.81 I kappaRk, wVel, salt,
669 jmc 1.63 U gS,
670     I bi, bj, myTime, myIter, myThid )
671     ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
672     #else /* INCLUDE_IMPLVERTADV_CODE */
673     IF ( saltStepping .AND. implicitDiffusion ) THEN
674     #endif /* INCLUDE_IMPLVERTADV_CODE */
675 adcroft 1.1 #ifdef ALLOW_AUTODIFF_TAMC
676 jmc 1.81 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
677 heimbach 1.30 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
678 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
679 jmc 1.63 CALL IMPLDIFF(
680 adcroft 1.1 I bi, bj, iMin, iMax, jMin, jMax,
681 jmc 1.87 I GAD_SALINITY, kappaRk, recip_hFacC,
682 adcroft 1.7 U gS,
683 adcroft 1.1 I myThid )
684 jmc 1.63 ENDIF
685 adcroft 1.1
686 adcroft 1.17 #ifdef ALLOW_PTRACERS
687 jmc 1.86 IF ( usePTRACERS ) THEN
688     C-- Vertical advection/diffusion (implicit) for passive tracers
689     CALL PTRACERS_IMPLICIT(
690     U kappaRk,
691     I bi, bj, myTime, myIter, myThid )
692 jmc 1.63 ENDIF
693 adcroft 1.17 #endif /* ALLOW_PTRACERS */
694    
695 adcroft 1.1 #ifdef ALLOW_OBCS
696     C-- Apply open boundary conditions
697 jmc 1.63 IF ( ( implicitDiffusion
698     & .OR. tempImplVertAdv
699     & .OR. saltImplVertAdv
700     & ) .AND. useOBCS ) THEN
701 adcroft 1.1 DO K=1,Nr
702 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
703 adcroft 1.1 ENDDO
704 jmc 1.63 ENDIF
705 adcroft 1.1 #endif /* ALLOW_OBCS */
706    
707 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
708 adcroft 1.1
709 jmc 1.39 C-- end bi,bj loops.
710 adcroft 1.1 ENDDO
711     ENDDO
712 adcroft 1.17
713 edhill 1.56 #ifdef ALLOW_DEBUG
714 adcroft 1.17 If (debugMode) THEN
715     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
716     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
717     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
718     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
719     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
720 jmc 1.90 CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
721     CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
722     #ifndef ALLOW_ADAMSBASHFORTH_3
723     CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
724     CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
725     #endif
726 adcroft 1.18 #ifdef ALLOW_PTRACERS
727     IF ( usePTRACERS ) THEN
728     CALL PTRACERS_DEBUG(myThid)
729     ENDIF
730     #endif /* ALLOW_PTRACERS */
731 adcroft 1.17 ENDIF
732 jmc 1.90 #endif /* ALLOW_DEBUG */
733 adcroft 1.40
734 edhill 1.56 #ifdef ALLOW_DEBUG
735 heimbach 1.43 IF ( debugLevel .GE. debLevB )
736 jmc 1.63 & CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
737 adcroft 1.17 #endif
738 adcroft 1.1
739 jmc 1.83 #endif /* ALLOW_GENERIC_ADVDIFF */
740    
741 adcroft 1.1 RETURN
742     END

  ViewVC Help
Powered by ViewVC 1.1.22