/[MITgcm]/MITgcm/model/src/thermodynamics.F
ViewVC logotype

Contents of /MITgcm/model/src/thermodynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.119 - (show annotations) (download)
Tue Oct 23 16:55:46 2007 UTC (16 years, 8 months ago) by heimbach
Branch: MAIN
Changes since 1.118: +3 -1 lines
Add one ifdef.

1 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.118 2007/10/19 14:41:39 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6 #ifdef ALLOW_GENERIC_ADVDIFF
7 # include "GAD_OPTIONS.h"
8 #endif
9
10 #ifdef ALLOW_AUTODIFF_TAMC
11 # ifdef ALLOW_GMREDI
12 # include "GMREDI_OPTIONS.h"
13 # endif
14 # ifdef ALLOW_KPP
15 # include "KPP_OPTIONS.h"
16 # endif
17 #endif /* ALLOW_AUTODIFF_TAMC */
18
19 CBOP
20 C !ROUTINE: THERMODYNAMICS
21 C !INTERFACE:
22 SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
23 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 IMPLICIT NONE
77 C == Global variables ===
78 #include "SIZE.h"
79 #include "EEPARAMS.h"
80 #include "PARAMS.h"
81 #include "RESTART.h"
82 #include "DYNVARS.h"
83 #include "GRID.h"
84 #ifdef ALLOW_GENERIC_ADVDIFF
85 # include "GAD.h"
86 # include "GAD_SOM_VARS.h"
87 #endif
88 #ifdef ALLOW_PTRACERS
89 # include "PTRACERS_SIZE.h"
90 # include "PTRACERS.h"
91 #endif
92 #ifdef ALLOW_TIMEAVE
93 # include "TIMEAVE_STATV.h"
94 #endif
95
96 #ifdef ALLOW_AUTODIFF_TAMC
97 # include "tamc.h"
98 # include "tamc_keys.h"
99 # include "FFIELDS.h"
100 # include "SURFACE.h"
101 # include "EOS.h"
102 # ifdef ALLOW_KPP
103 # include "KPP.h"
104 # endif
105 # ifdef ALLOW_GMREDI
106 # include "GMREDI.h"
107 # endif
108 # ifdef ALLOW_EBM
109 # include "EBM.h"
110 # endif
111 #endif /* ALLOW_AUTODIFF_TAMC */
112
113 C !INPUT/OUTPUT PARAMETERS:
114 C == Routine arguments ==
115 C myTime - Current time in simulation
116 C myIter - Current iteration number in simulation
117 C myThid - Thread number for this instance of the routine.
118 _RL myTime
119 INTEGER myIter
120 INTEGER myThid
121
122 #ifdef ALLOW_GENERIC_ADVDIFF
123 C !LOCAL VARIABLES:
124 C == Local variables
125 C xA, yA - Per block temporaries holding face areas
126 C uFld, vFld, wFld - Local copy of velocity field (3 components)
127 C uTrans, vTrans, rTrans - Per block temporaries holding flow transport
128 C o uTrans: Zonal transport
129 C o vTrans: Meridional transport
130 C o rTrans: Vertical transport
131 C rTransKp1 o vertical volume transp. at interface k+1
132 C maskUp o maskUp: land/water mask for W points
133 C fVer[STUV] o fVer: Vertical flux term - note fVer
134 C is "pipelined" in the vertical
135 C so we need an fVer for each
136 C variable.
137 C kappaRT, - Total diffusion in vertical at level k, for T and S
138 C kappaRS (background + spatially varying, isopycnal term).
139 C kappaRTr - Total diffusion in vertical at level k,
140 C for each passive Tracer
141 C kappaRk - Total diffusion in vertical, all levels, 1 tracer
142 C useVariableK = T when vertical diffusion is not constant
143 C iMin, iMax - Ranges and sub-block indices on which calculations
144 C jMin, jMax are applied.
145 C bi, bj
146 C k, kup, - Index for layer above and below. kup and kDown
147 C kDown, km1 are switched with layer to be the appropriate
148 C index into fVerTerm.
149 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157 _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
160 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
161 _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
162 _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
163 #ifdef ALLOW_PTRACERS
164 _RL fVerP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
165 _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
166 #endif
167 _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
168 INTEGER iMin, iMax
169 INTEGER jMin, jMax
170 INTEGER bi, bj
171 INTEGER i, j
172 INTEGER k, km1, kup, kDown
173 #ifdef ALLOW_ADAMSBASHFORTH_3
174 INTEGER iterNb, m1, m2
175 #endif
176 #ifdef ALLOW_TIMEAVE
177 LOGICAL useVariableK
178 #endif
179 #ifdef ALLOW_PTRACERS
180 INTEGER iTracer, ip
181 #endif
182
183 CEOP
184
185 #ifdef ALLOW_DEBUG
186 IF ( debugLevel .GE. debLevB )
187 & CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
188 #endif
189
190 #ifdef ALLOW_AUTODIFF_TAMC
191 C-- dummy statement to end declaration part
192 ikey = 1
193 itdkey = 1
194 #endif /* ALLOW_AUTODIFF_TAMC */
195
196 #ifdef ALLOW_AUTODIFF_TAMC
197 C-- HPF directive to help TAMC
198 CHPF$ INDEPENDENT
199 #endif /* ALLOW_AUTODIFF_TAMC */
200
201 C-- Compute correction at the surface for Lin Free Surf.
202 #ifdef ALLOW_AUTODIFF_TAMC
203 TsurfCor = 0. _d 0
204 SsurfCor = 0. _d 0
205 CADJ STORE theta,salt,wvel = comlev1, key = ikey_dynamics
206 #endif
207 IF (linFSConserveTr) THEN
208 #ifdef ALLOW_AUTODIFF_TAMC
209 CADJ STORE theta,salt,wvel = comlev1, key = ikey_dynamics
210 #endif
211 CALL CALC_WSURF_TR(theta,salt,wVel,
212 & myTime,myIter,myThid)
213 ENDIF
214
215 DO bj=myByLo(myThid),myByHi(myThid)
216
217 #ifdef ALLOW_AUTODIFF_TAMC
218 C-- HPF directive to help TAMC
219 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
220 CHPF$& ,utrans,vtrans,xA,yA
221 CHPF$& ,kappaRT,kappaRS
222 CHPF$& )
223 # ifdef ALLOW_PTRACERS
224 CHPF$ INDEPENDENT, NEW (fVerP,kappaRTr)
225 # endif
226 #endif /* ALLOW_AUTODIFF_TAMC */
227
228 DO bi=myBxLo(myThid),myBxHi(myThid)
229
230 #ifdef ALLOW_AUTODIFF_TAMC
231 act1 = bi - myBxLo(myThid)
232 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
233 act2 = bj - myByLo(myThid)
234 max2 = myByHi(myThid) - myByLo(myThid) + 1
235 act3 = myThid - 1
236 max3 = nTx*nTy
237 act4 = ikey_dynamics - 1
238 itdkey = (act1 + 1) + act2*max1
239 & + act3*max1*max2
240 & + act4*max1*max2*max3
241 #endif /* ALLOW_AUTODIFF_TAMC */
242
243 C-- Set up work arrays with valid (i.e. not NaN) values
244 C These inital values do not alter the numerical results. They
245 C just ensure that all memory references are to valid floating
246 C point numbers. This prevents spurious hardware signals due to
247 C uninitialised but inert locations.
248
249 DO j=1-OLy,sNy+OLy
250 DO i=1-OLx,sNx+OLx
251 xA(i,j) = 0. _d 0
252 yA(i,j) = 0. _d 0
253 uTrans(i,j) = 0. _d 0
254 vTrans(i,j) = 0. _d 0
255 rTrans (i,j) = 0. _d 0
256 rTransKp1(i,j) = 0. _d 0
257 fVerT (i,j,1) = 0. _d 0
258 fVerT (i,j,2) = 0. _d 0
259 fVerS (i,j,1) = 0. _d 0
260 fVerS (i,j,2) = 0. _d 0
261 kappaRT(i,j) = 0. _d 0
262 kappaRS(i,j) = 0. _d 0
263 ENDDO
264 ENDDO
265
266 DO k=1,Nr
267 DO j=1-OLy,sNy+OLy
268 DO i=1-OLx,sNx+OLx
269 C This is currently also used by IVDC and Diagnostics
270 kappaRk(i,j,k) = 0. _d 0
271 C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
272 gT(i,j,k,bi,bj) = 0. _d 0
273 gS(i,j,k,bi,bj) = 0. _d 0
274 ENDDO
275 ENDDO
276 ENDDO
277
278 #ifdef ALLOW_PTRACERS
279 IF ( usePTRACERS ) THEN
280 DO ip=1,PTRACERS_num
281 DO j=1-OLy,sNy+OLy
282 DO i=1-OLx,sNx+OLx
283 fVerP (i,j,1,ip) = 0. _d 0
284 fVerP (i,j,2,ip) = 0. _d 0
285 kappaRTr(i,j,ip) = 0. _d 0
286 ENDDO
287 ENDDO
288 ENDDO
289 C- set tracer tendency to zero:
290 DO iTracer=1,PTRACERS_num
291 DO k=1,Nr
292 DO j=1-OLy,sNy+OLy
293 DO i=1-OLx,sNx+OLx
294 gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
295 ENDDO
296 ENDDO
297 ENDDO
298 ENDDO
299 ENDIF
300 #endif
301
302 #ifdef ALLOW_ADAMSBASHFORTH_3
303 C- Apply AB on T,S :
304 iterNb = myIter
305 IF (staggerTimeStep) iterNb = myIter - 1
306 m1 = 1 + MOD(iterNb+1,2)
307 m2 = 1 + MOD( iterNb ,2)
308 C compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
309 IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
310 I bi, bj, 0,
311 U theta, gtNm,
312 I tempStartAB, iterNb, myThid )
313 C compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
314 IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
315 I bi, bj, 0,
316 U salt, gsNm,
317 I saltStartAB, iterNb, myThid )
318 #endif /* ALLOW_ADAMSBASHFORTH_3 */
319
320 c iMin = 1-OLx
321 c iMax = sNx+OLx
322 c jMin = 1-OLy
323 c jMax = sNy+OLy
324
325 #ifdef ALLOW_AUTODIFF_TAMC
326 cph avoids recomputation of integrate_for_w
327 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
328 #endif /* ALLOW_AUTODIFF_TAMC */
329
330 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
331 C-- MOST of THERMODYNAMICS will be disabled
332 #ifndef SINGLE_LAYER_MODE
333
334 #ifdef ALLOW_AUTODIFF_TAMC
335 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
336 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
337 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
338 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
339 # if ((defined ALLOW_DEPTH_CONTROL) || (defined NONLIN_FRSURF))
340 CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
341 CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
342 # endif
343 #endif /* ALLOW_AUTODIFF_TAMC */
344
345 #ifndef DISABLE_MULTIDIM_ADVECTION
346 C-- Some advection schemes are better calculated using a multi-dimensional
347 C method in the absence of any other terms and, if used, is done here.
348 C
349 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
350 C The default is to use multi-dimensinal advection for non-linear advection
351 C schemes. However, for the sake of efficiency of the adjoint it is necessary
352 C to be able to exclude this scheme to avoid excessive storage and
353 C recomputation. It *is* differentiable, if you need it.
354 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
355 C disable this section of code.
356 #ifdef GAD_ALLOW_SOM_ADVECT
357 IF ( tempSOM_Advection ) THEN
358 #ifdef ALLOW_DEBUG
359 IF ( debugLevel .GE. debLevB )
360 & CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
361 #endif
362 CALL GAD_SOM_ADVECT(
363 I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
364 I GAD_TEMPERATURE,
365 I uVel, vVel, wVel, theta,
366 U som_T,
367 O gT,
368 I bi,bj,myTime,myIter,myThid)
369 ELSEIF (tempMultiDimAdvec) THEN
370 #else /* GAD_ALLOW_SOM_ADVECT */
371 IF (tempMultiDimAdvec) THEN
372 #endif /* GAD_ALLOW_SOM_ADVECT */
373 #ifdef ALLOW_DEBUG
374 IF ( debugLevel .GE. debLevB )
375 & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
376 #endif
377 CALL GAD_ADVECTION(
378 I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
379 I GAD_TEMPERATURE,
380 I uVel, vVel, wVel, theta,
381 O gT,
382 I bi,bj,myTime,myIter,myThid)
383 ENDIF
384 #ifdef GAD_ALLOW_SOM_ADVECT
385 IF ( saltSOM_Advection ) THEN
386 #ifdef ALLOW_DEBUG
387 IF ( debugLevel .GE. debLevB )
388 & CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
389 #endif
390 CALL GAD_SOM_ADVECT(
391 I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
392 I GAD_SALINITY,
393 I uVel, vVel, wVel, salt,
394 U som_S,
395 O gS,
396 I bi,bj,myTime,myIter,myThid)
397 ELSEIF (saltMultiDimAdvec) THEN
398 #else /* GAD_ALLOW_SOM_ADVECT */
399 IF (saltMultiDimAdvec) THEN
400 #endif /* GAD_ALLOW_SOM_ADVECT */
401 #ifdef ALLOW_DEBUG
402 IF ( debugLevel .GE. debLevB )
403 & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
404 #endif
405 CALL GAD_ADVECTION(
406 I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
407 I GAD_SALINITY,
408 I uVel, vVel, wVel, salt,
409 O gS,
410 I bi,bj,myTime,myIter,myThid)
411 ENDIF
412
413 C Since passive tracers are configurable separately from T,S we
414 C call the multi-dimensional method for PTRACERS regardless
415 C of whether multiDimAdvection is set or not.
416 #ifdef ALLOW_PTRACERS
417 IF ( usePTRACERS ) THEN
418 #ifdef ALLOW_DEBUG
419 IF ( debugLevel .GE. debLevB )
420 & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
421 #endif
422 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
423 ENDIF
424 #endif /* ALLOW_PTRACERS */
425 #endif /* DISABLE_MULTIDIM_ADVECTION */
426
427 #ifdef ALLOW_DEBUG
428 IF ( debugLevel .GE. debLevB )
429 & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
430 #endif
431
432 C-- Start of thermodynamics loop
433 DO k=Nr,1,-1
434 #ifdef ALLOW_AUTODIFF_TAMC
435 C? Patrick Is this formula correct?
436 cph Yes, but I rewrote it.
437 cph Also, the kappaR? need the index and subscript k!
438 kkey = (itdkey-1)*Nr + k
439 #endif /* ALLOW_AUTODIFF_TAMC */
440
441 C-- km1 Points to level above k (=k-1)
442 C-- kup Cycles through 1,2 to point to layer above
443 C-- kDown Cycles through 2,1 to point to current layer
444
445 km1 = MAX(1,k-1)
446 kup = 1+MOD(k+1,2)
447 kDown= 1+MOD(k,2)
448
449 iMin = 1-OLx
450 iMax = sNx+OLx
451 jMin = 1-OLy
452 jMax = sNy+OLy
453
454 IF (k.EQ.Nr) THEN
455 DO j=1-Oly,sNy+Oly
456 DO i=1-Olx,sNx+Olx
457 rTransKp1(i,j) = 0. _d 0
458 ENDDO
459 ENDDO
460 ELSE
461 DO j=1-Oly,sNy+Oly
462 DO i=1-Olx,sNx+Olx
463 rTransKp1(i,j) = rTrans(i,j)
464 ENDDO
465 ENDDO
466 ENDIF
467 #ifdef ALLOW_AUTODIFF_TAMC
468 CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
469 #endif
470
471 C-- Get temporary terms used by tendency routines :
472 C- Calculate horizontal "volume transport" through tracer cell face
473 C anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport)
474 CALL CALC_COMMON_FACTORS (
475 I uVel, vVel,
476 O uFld, vFld, uTrans, vTrans, xA, yA,
477 I k,bi,bj, myThid )
478
479 C- Calculate vertical "volume transport" through tracer cell face
480 IF (k.EQ.1) THEN
481 C- Surface interface :
482 DO j=1-Oly,sNy+Oly
483 DO i=1-Olx,sNx+Olx
484 wFld(i,j) = 0. _d 0
485 maskUp(i,j) = 0. _d 0
486 rTrans(i,j) = 0. _d 0
487 ENDDO
488 ENDDO
489 ELSE
490 C- Interior interface :
491 C anelastic: rTrans is scaled by rhoFacF (~ mass transport)
492 DO j=1-Oly,sNy+Oly
493 DO i=1-Olx,sNx+Olx
494 wFld(i,j) = wVel(i,j,k,bi,bj)
495 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
496 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
497 & *deepFac2F(k)*rhoFacF(k)
498 ENDDO
499 ENDDO
500 ENDIF
501
502 #ifdef ALLOW_GMREDI
503 C-- Residual transp = Bolus transp + Eulerian transp
504 IF (useGMRedi) THEN
505 CALL GMREDI_CALC_UVFLOW(
506 U uFld, vFld, uTrans, vTrans,
507 I k, bi, bj, myThid )
508 IF (K.GE.2) THEN
509 CALL GMREDI_CALC_WFLOW(
510 U wFld, rTrans,
511 I k, bi, bj, myThid )
512 ENDIF
513 ENDIF
514 # ifdef ALLOW_AUTODIFF_TAMC
515 CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
516 CADJ STORE wfld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
517 # ifdef GM_BOLUS_ADVEC
518 CADJ STORE ufld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
519 CADJ STORE vfld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
520 CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
521 CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
522 # endif
523 # endif /* ALLOW_AUTODIFF_TAMC */
524 #endif /* ALLOW_GMREDI */
525
526 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
527 C-- Calculate the total vertical diffusivity
528 IF ( .NOT.implicitDiffusion ) THEN
529 CALL CALC_DIFFUSIVITY(
530 I bi,bj,iMin,iMax,jMin,jMax,k,
531 I maskUp,
532 O kappaRT,kappaRS,
533 I myThid)
534 ENDIF
535 # ifdef ALLOW_AUTODIFF_TAMC
536 CADJ STORE kappaRT(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
537 CADJ STORE kappaRS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
538 # endif /* ALLOW_AUTODIFF_TAMC */
539 #endif
540
541 iMin = 1-OLx+2
542 iMax = sNx+OLx-1
543 jMin = 1-OLy+2
544 jMax = sNy+OLy-1
545
546 C-- Calculate active tracer tendencies (gT,gS,...)
547 C and step forward storing result in gT, gS, etc.
548 C--
549 # ifdef ALLOW_AUTODIFF_TAMC
550 # if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
551 # ifdef GM_NON_UNITY_DIAGONAL
552 CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
553 CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
554 # endif
555 # ifdef GM_EXTRA_DIAGONAL
556 CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
557 CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
558 # endif
559 # endif
560 # endif /* ALLOW_AUTODIFF_TAMC */
561 C
562 #ifdef ALLOW_AUTODIFF_TAMC
563 # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
564 cph-test
565 CADJ STORE uFld(:,:), vFld(:,:), wFld(:,:)
566 CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
567 CADJ STORE uTrans(:,:), vTrans(:,:)
568 CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
569 CADJ STORE xA(:,:), yA(:,:)
570 CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
571 # endif
572 #endif /* ALLOW_AUTODIFF_TAMC */
573 C
574 IF ( tempStepping ) THEN
575 #ifdef ALLOW_AUTODIFF_TAMC
576 CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
577 # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
578 CADJ STORE gt(:,:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
579 CADJ STORE fvert(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
580 # endif
581 #endif /* ALLOW_AUTODIFF_TAMC */
582 CALL CALC_GT(
583 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
584 I xA, yA, maskUp, uFld, vFld, wFld,
585 I uTrans, vTrans, rTrans, rTransKp1,
586 I kappaRT,
587 U fVerT,
588 I myTime,myIter,myThid)
589 #ifdef ALLOW_ADAMSBASHFORTH_3
590 IF ( AdamsBashforth_T ) THEN
591 CALL TIMESTEP_TRACER(
592 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
593 I gtNm(1-Olx,1-Oly,1,1,1,m2),
594 U gT,
595 I myIter, myThid)
596 ELSE
597 #endif
598 CALL TIMESTEP_TRACER(
599 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
600 I theta,
601 U gT,
602 I myIter, myThid)
603 #ifdef ALLOW_ADAMSBASHFORTH_3
604 ENDIF
605 #endif
606 ENDIF
607
608 IF ( saltStepping ) THEN
609 #ifdef ALLOW_AUTODIFF_TAMC
610 CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
611 # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
612 CADJ STORE gs(:,:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
613 CADJ STORE fvers(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
614 # endif
615 #endif /* ALLOW_AUTODIFF_TAMC */
616 CALL CALC_GS(
617 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
618 I xA, yA, maskUp, uFld, vFld, wFld,
619 I uTrans, vTrans, rTrans, rTransKp1,
620 I kappaRS,
621 U fVerS,
622 I myTime,myIter,myThid)
623 #ifdef ALLOW_ADAMSBASHFORTH_3
624 IF ( AdamsBashforth_S ) THEN
625 CALL TIMESTEP_TRACER(
626 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
627 I gsNm(1-Olx,1-Oly,1,1,1,m2),
628 U gS,
629 I myIter, myThid)
630 ELSE
631 #endif
632 CALL TIMESTEP_TRACER(
633 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
634 I salt,
635 U gS,
636 I myIter, myThid)
637 #ifdef ALLOW_ADAMSBASHFORTH_3
638 ENDIF
639 #endif
640 ENDIF
641
642 #ifdef ALLOW_PTRACERS
643 IF ( usePTRACERS ) THEN
644 IF ( .NOT.implicitDiffusion ) THEN
645 CALL PTRACERS_CALC_DIFF(
646 I bi,bj,iMin,iMax,jMin,jMax,k,
647 I maskUp,
648 O kappaRTr,
649 I myThid)
650 ENDIF
651 # ifdef ALLOW_AUTODIFF_TAMC
652 CADJ STORE kappaRTr(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
653 # endif /* ALLOW_AUTODIFF_TAMC */
654 CALL PTRACERS_INTEGRATE(
655 I bi,bj,k,
656 I xA, yA, maskUp, uFld, vFld, wFld,
657 I uTrans, vTrans, rTrans, rTransKp1,
658 I kappaRTr,
659 U fVerP,
660 I myTime,myIter,myThid)
661 ENDIF
662 #endif /* ALLOW_PTRACERS */
663
664 #ifdef ALLOW_OBCS
665 C-- Apply open boundary conditions
666 IF (useOBCS) THEN
667 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
668 END IF
669 #endif /* ALLOW_OBCS */
670
671 C-- Freeze water
672 C this bit of code is left here for backward compatibility.
673 C freezing at surface level has been moved to FORWARD_STEP
674 IF ( useOldFreezing .AND. .NOT. useSEAICE
675 & .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
676 #ifdef ALLOW_AUTODIFF_TAMC
677 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
678 CADJ & , key = kkey, byte = isbyte
679 #endif /* ALLOW_AUTODIFF_TAMC */
680 CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
681 ENDIF
682
683 C-- end of thermodynamic k loop (Nr:1)
684 ENDDO
685
686 C All explicit advection/diffusion/sources should now be
687 C done. The updated tracer field is in gPtr. Accumalate
688 C explicit tendency and also reset gPtr to initial tracer
689 C field for implicit matrix calculation
690
691 #ifdef ALLOW_MATRIX
692 IF (useMATRIX)
693 & CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
694 #endif
695
696 iMin = 1
697 iMax = sNx
698 jMin = 1
699 jMax = sNy
700
701 C-- Implicit vertical advection & diffusion
702 IF ( tempStepping .AND. implicitDiffusion ) THEN
703 CALL CALC_3D_DIFFUSIVITY(
704 I bi,bj,iMin,iMax,jMin,jMax,
705 I GAD_TEMPERATURE, useGMredi, useKPP,
706 O kappaRk,
707 I myThid)
708 ENDIF
709 #ifdef INCLUDE_IMPLVERTADV_CODE
710 IF ( tempImplVertAdv ) THEN
711 CALL GAD_IMPLICIT_R(
712 I tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
713 I kappaRk, wVel, theta,
714 U gT,
715 I bi, bj, myTime, myIter, myThid )
716 ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
717 #else /* INCLUDE_IMPLVERTADV_CODE */
718 IF ( tempStepping .AND. implicitDiffusion ) THEN
719 #endif /* INCLUDE_IMPLVERTADV_CODE */
720 #ifdef ALLOW_AUTODIFF_TAMC
721 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
722 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
723 #endif /* ALLOW_AUTODIFF_TAMC */
724 CALL IMPLDIFF(
725 I bi, bj, iMin, iMax, jMin, jMax,
726 I GAD_TEMPERATURE, kappaRk, recip_hFacC,
727 U gT,
728 I myThid )
729 ENDIF
730
731 #ifdef ALLOW_TIMEAVE
732 useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
733 & .OR. useGMredi .OR. ivdc_kappa.NE.0.
734 IF (taveFreq.GT.0. .AND. useVariableK ) THEN
735 IF (implicitDiffusion) THEN
736 CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
737 I Nr, 3, deltaTclock, bi, bj, myThid)
738 c ELSE
739 c CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
740 c I Nr, 3, deltaTclock, bi, bj, myThid)
741 ENDIF
742 ENDIF
743 #endif /* ALLOW_TIMEAVE */
744
745 IF ( saltStepping .AND. implicitDiffusion ) THEN
746 CALL CALC_3D_DIFFUSIVITY(
747 I bi,bj,iMin,iMax,jMin,jMax,
748 I GAD_SALINITY, useGMredi, useKPP,
749 O kappaRk,
750 I myThid)
751 ENDIF
752
753 #ifdef INCLUDE_IMPLVERTADV_CODE
754 IF ( saltImplVertAdv ) THEN
755 CALL GAD_IMPLICIT_R(
756 I saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
757 I kappaRk, wVel, salt,
758 U gS,
759 I bi, bj, myTime, myIter, myThid )
760 ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
761 #else /* INCLUDE_IMPLVERTADV_CODE */
762 IF ( saltStepping .AND. implicitDiffusion ) THEN
763 #endif /* INCLUDE_IMPLVERTADV_CODE */
764 #ifdef ALLOW_AUTODIFF_TAMC
765 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
766 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
767 #endif /* ALLOW_AUTODIFF_TAMC */
768 CALL IMPLDIFF(
769 I bi, bj, iMin, iMax, jMin, jMax,
770 I GAD_SALINITY, kappaRk, recip_hFacC,
771 U gS,
772 I myThid )
773 ENDIF
774
775 #ifdef ALLOW_PTRACERS
776 IF ( usePTRACERS ) THEN
777 C-- Vertical advection/diffusion (implicit) for passive tracers
778 CALL PTRACERS_IMPLICIT(
779 U kappaRk,
780 I bi, bj, myTime, myIter, myThid )
781 ENDIF
782 #endif /* ALLOW_PTRACERS */
783
784 #ifdef ALLOW_OBCS
785 C-- Apply open boundary conditions
786 IF ( ( implicitDiffusion
787 & .OR. tempImplVertAdv
788 & .OR. saltImplVertAdv
789 & ) .AND. useOBCS ) THEN
790 DO K=1,Nr
791 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
792 ENDDO
793 ENDIF
794 #endif /* ALLOW_OBCS */
795
796 #endif /* SINGLE_LAYER_MODE */
797
798 C-- end bi,bj loops.
799 ENDDO
800 ENDDO
801
802 #ifdef ALLOW_DEBUG
803 If (debugMode) THEN
804 CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
805 CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
806 CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
807 CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
808 CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
809 CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
810 CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
811 #ifndef ALLOW_ADAMSBASHFORTH_3
812 CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
813 CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
814 #endif
815 #ifdef ALLOW_PTRACERS
816 IF ( usePTRACERS ) THEN
817 CALL PTRACERS_DEBUG(myThid)
818 ENDIF
819 #endif /* ALLOW_PTRACERS */
820 ENDIF
821 #endif /* ALLOW_DEBUG */
822
823 #ifdef ALLOW_DEBUG
824 IF ( debugLevel .GE. debLevB )
825 & CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
826 #endif
827
828 #endif /* ALLOW_GENERIC_ADVDIFF */
829
830 RETURN
831 END

  ViewVC Help
Powered by ViewVC 1.1.22