/[MITgcm]/MITgcm_contrib/jscott/code_rafmod/thermodynamics.F
ViewVC logotype

Contents of /MITgcm_contrib/jscott/code_rafmod/thermodynamics.F

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


Revision 1.1 - (show annotations) (download)
Tue Aug 21 16:34:17 2007 UTC (17 years, 10 months ago) by jscott
Branch: MAIN
code directory for crude ML horizontal mixing scheme

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

  ViewVC Help
Powered by ViewVC 1.1.22