/[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.115 - (show annotations) (download)
Wed Jan 24 08:06:25 2007 UTC (17 years, 5 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59, checkpoint58y_post, checkpoint58w_post, checkpoint59e, checkpoint59d, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint58v_post, checkpoint58x_post
Changes since 1.114: +5 -3 lines
Updates for NLFS adjoint.

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

  ViewVC Help
Powered by ViewVC 1.1.22