/[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.36 - (show annotations) (download)
Tue Jan 21 19:19:45 2003 UTC (21 years, 4 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint48b_post, checkpoint48c_pre, checkpoint48a_post, checkpoint47j_post, checkpoint48c_post, checkpoint48
Changes since 1.35: +4 -1 lines
Dummy init. of visbeckK for TAF.

1 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.35 2003/01/10 23:41:15 heimbach Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5 #ifdef ALLOW_AUTODIFF_TAMC
6 # ifdef ALLOW_GMREDI
7 # include "GMREDI_OPTIONS.h"
8 # endif
9 # ifdef ALLOW_KPP
10 # include "KPP_OPTIONS.h"
11 # endif
12 cswdice --- add ----
13 #ifdef ALLOW_THERM_SEAICE
14 #include "ICE.h"
15 #endif
16 cswdice ------
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 #include "GAD.h"
84 #ifdef ALLOW_PASSIVE_TRACER
85 #include "TR1.h"
86 #endif
87 #ifdef ALLOW_AUTODIFF_TAMC
88 # include "tamc.h"
89 # include "tamc_keys.h"
90 # include "FFIELDS.h"
91 # include "EOS.h"
92 # ifdef ALLOW_KPP
93 # include "KPP.h"
94 # endif
95 # ifdef ALLOW_GMREDI
96 # include "GMREDI.h"
97 # endif
98 #endif /* ALLOW_AUTODIFF_TAMC */
99 #ifdef ALLOW_TIMEAVE
100 #include "TIMEAVE_STATV.h"
101 #endif
102
103 C !INPUT/OUTPUT PARAMETERS:
104 C == Routine arguments ==
105 C myTime - Current time in simulation
106 C myIter - Current iteration number in simulation
107 C myThid - Thread number for this instance of the routine.
108 _RL myTime
109 INTEGER myIter
110 INTEGER myThid
111
112 C !LOCAL VARIABLES:
113 C == Local variables
114 C xA, yA - Per block temporaries holding face areas
115 C uTrans, vTrans, rTrans - Per block temporaries holding flow
116 C transport
117 C o uTrans: Zonal transport
118 C o vTrans: Meridional transport
119 C o rTrans: Vertical transport
120 C maskUp o maskUp: land/water mask for W points
121 C fVer[STUV] o fVer: Vertical flux term - note fVer
122 C is "pipelined" in the vertical
123 C so we need an fVer for each
124 C variable.
125 C rhoK, rhoKM1 - Density at current level, and level above
126 C phiHyd - Hydrostatic part of the potential phiHydi.
127 C In z coords phiHydiHyd is the hydrostatic
128 C Potential (=pressure/rho0) anomaly
129 C In p coords phiHydiHyd is the geopotential
130 C surface height anomaly.
131 C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
132 C phiSurfY or geopotentiel (atmos) in X and Y direction
133 C KappaRT, - Total diffusion in vertical for T and S.
134 C KappaRS (background + spatially varying, isopycnal term).
135 C iMin, iMax - Ranges and sub-block indices on which calculations
136 C jMin, jMax are applied.
137 C bi, bj
138 C k, kup, - Index for layer above and below. kup and kDown
139 C kDown, km1 are switched with layer to be the appropriate
140 C index into fVerTerm.
141 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
148 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
149 _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
150 _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
151 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154 _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155 _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
156 _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
157 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
158 _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
159 _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
160 C This is currently used by IVDC and Diagnostics
161 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
162 INTEGER iMin, iMax
163 INTEGER jMin, jMax
164 INTEGER bi, bj
165 INTEGER i, j
166 INTEGER k, km1, kup, kDown
167
168 CEOP
169
170 #ifdef ALLOW_AUTODIFF_TAMC
171 C-- dummy statement to end declaration part
172 ikey = 1
173 itdkey = 1
174 #endif /* ALLOW_AUTODIFF_TAMC */
175
176 C-- Set up work arrays with valid (i.e. not NaN) values
177 C These inital values do not alter the numerical results. They
178 C just ensure that all memory references are to valid floating
179 C point numbers. This prevents spurious hardware signals due to
180 C uninitialised but inert locations.
181 DO j=1-OLy,sNy+OLy
182 DO i=1-OLx,sNx+OLx
183 xA(i,j) = 0. _d 0
184 yA(i,j) = 0. _d 0
185 uTrans(i,j) = 0. _d 0
186 vTrans(i,j) = 0. _d 0
187 rhok (i,j) = 0. _d 0
188 phiSurfX(i,j) = 0. _d 0
189 phiSurfY(i,j) = 0. _d 0
190 ENDDO
191 ENDDO
192
193
194 #ifdef ALLOW_AUTODIFF_TAMC
195 C-- HPF directive to help TAMC
196 CHPF$ INDEPENDENT
197 #endif /* ALLOW_AUTODIFF_TAMC */
198
199 DO bj=myByLo(myThid),myByHi(myThid)
200
201 #ifdef ALLOW_AUTODIFF_TAMC
202 C-- HPF directive to help TAMC
203 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
204 CHPF$& ,phiHyd,utrans,vtrans,xA,yA
205 CHPF$& ,KappaRT,KappaRS
206 CHPF$& )
207 #endif /* ALLOW_AUTODIFF_TAMC */
208
209 DO bi=myBxLo(myThid),myBxHi(myThid)
210
211 #ifdef ALLOW_AUTODIFF_TAMC
212 act1 = bi - myBxLo(myThid)
213 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
214 act2 = bj - myByLo(myThid)
215 max2 = myByHi(myThid) - myByLo(myThid) + 1
216 act3 = myThid - 1
217 max3 = nTx*nTy
218 act4 = ikey_dynamics - 1
219 itdkey = (act1 + 1) + act2*max1
220 & + act3*max1*max2
221 & + act4*max1*max2*max3
222 #endif /* ALLOW_AUTODIFF_TAMC */
223
224 C-- Set up work arrays that need valid initial values
225 DO j=1-OLy,sNy+OLy
226 DO i=1-OLx,sNx+OLx
227 rTrans (i,j) = 0. _d 0
228 fVerT (i,j,1) = 0. _d 0
229 fVerT (i,j,2) = 0. _d 0
230 fVerS (i,j,1) = 0. _d 0
231 fVerS (i,j,2) = 0. _d 0
232 fVerTr1(i,j,1) = 0. _d 0
233 fVerTr1(i,j,2) = 0. _d 0
234 rhoKM1 (i,j) = 0. _d 0
235 ENDDO
236 ENDDO
237
238 DO k=1,Nr
239 DO j=1-OLy,sNy+OLy
240 DO i=1-OLx,sNx+OLx
241 C This is currently also used by IVDC and Diagnostics
242 phiHyd(i,j,k) = 0. _d 0
243 sigmaX(i,j,k) = 0. _d 0
244 sigmaY(i,j,k) = 0. _d 0
245 sigmaR(i,j,k) = 0. _d 0
246 ConvectCount(i,j,k) = 0.
247 KappaRT(i,j,k) = 0. _d 0
248 KappaRS(i,j,k) = 0. _d 0
249 #ifdef ALLOW_AUTODIFF_TAMC
250 cph all the following init. are necessary for TAF
251 cph although some of these are re-initialised later.
252 gT(i,j,k,bi,bj) = 0. _d 0
253 gS(i,j,k,bi,bj) = 0. _d 0
254 # ifdef ALLOW_PASSIVE_TRACER
255 gTr1(i,j,k,bi,bj) = 0. _d 0
256 # endif
257 # ifdef ALLOW_GMREDI
258 Kwx(i,j,k,bi,bj) = 0. _d 0
259 Kwy(i,j,k,bi,bj) = 0. _d 0
260 Kwz(i,j,k,bi,bj) = 0. _d 0
261 # ifdef GM_NON_UNITY_DIAGONAL
262 Kux(i,j,k,bi,bj) = 0. _d 0
263 Kvy(i,j,k,bi,bj) = 0. _d 0
264 # endif
265 # ifdef GM_EXTRA_DIAGONAL
266 Kuz(i,j,k,bi,bj) = 0. _d 0
267 Kvz(i,j,k,bi,bj) = 0. _d 0
268 # endif
269 # ifdef GM_BOLUS_ADVEC
270 GM_PsiX(i,j,k,bi,bj) = 0. _d 0
271 GM_PsiY(i,j,k,bi,bj) = 0. _d 0
272 # endif
273 # ifdef GM_VISBECK_VARIABLE_K
274 VisbeckK(i,j,bi,bj) = 0. _d 0
275 # endif
276 # endif /* ALLOW_GMREDI */
277 #endif /* ALLOW_AUTODIFF_TAMC */
278 ENDDO
279 ENDDO
280 ENDDO
281
282 iMin = 1-OLx
283 iMax = sNx+OLx
284 jMin = 1-OLy
285 jMax = sNy+OLy
286
287
288 #ifdef ALLOW_AUTODIFF_TAMC
289 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
290 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
291 #ifdef ALLOW_KPP
292 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
293 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
294 #endif
295 #endif /* ALLOW_AUTODIFF_TAMC */
296
297 C-- Start of diagnostic loop
298 DO k=Nr,1,-1
299
300 #ifdef ALLOW_AUTODIFF_TAMC
301 C? Patrick, is this formula correct now that we change the loop range?
302 C? Do we still need this?
303 cph kkey formula corrected.
304 cph Needed for rhok, rhokm1, in the case useGMREDI.
305 kkey = (itdkey-1)*Nr + k
306 #endif /* ALLOW_AUTODIFF_TAMC */
307
308 C-- Integrate continuity vertically for vertical velocity
309 c CALL INTEGRATE_FOR_W(
310 c I bi, bj, k, uVel, vVel,
311 c O wVel,
312 c I myThid )
313
314 #ifdef ALLOW_OBCS
315 #ifdef ALLOW_NONHYDROSTATIC
316 C-- Apply OBC to W if in N-H mode
317 c IF (useOBCS.AND.nonHydrostatic) THEN
318 c CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
319 c ENDIF
320 #endif /* ALLOW_NONHYDROSTATIC */
321 #endif /* ALLOW_OBCS */
322
323 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
324 C-- MOST of THERMODYNAMICS will be disabled
325 #ifndef SINGLE_LAYER_MODE
326
327 C-- Calculate gradients of potential density for isoneutral
328 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
329 c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
330 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
331 #ifdef ALLOW_AUTODIFF_TAMC
332 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
333 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
334 CADJ STORE pressure(:,:,k,bi,bj) =
335 CADJ & comlev1_bibj_k, key=kkey, byte=isbyte
336 #endif /* ALLOW_AUTODIFF_TAMC */
337 CALL FIND_RHO(
338 I bi, bj, iMin, iMax, jMin, jMax, k, k,
339 I theta, salt,
340 O rhoK,
341 I myThid )
342
343 IF (k.GT.1) THEN
344 #ifdef ALLOW_AUTODIFF_TAMC
345 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
346 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
347 CADJ STORE pressure(:,:,k-1,bi,bj) =
348 CADJ & comlev1_bibj_k, key=kkey, byte=isbyte
349 #endif /* ALLOW_AUTODIFF_TAMC */
350 CALL FIND_RHO(
351 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
352 I theta, salt,
353 O rhoKm1,
354 I myThid )
355 ENDIF
356 CALL GRAD_SIGMA(
357 I bi, bj, iMin, iMax, jMin, jMax, k,
358 I rhoK, rhoKm1, rhoK,
359 O sigmaX, sigmaY, sigmaR,
360 I myThid )
361 ENDIF
362
363 #ifdef ALLOW_AUTODIFF_TAMC
364 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
365 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
366 #endif /* ALLOW_AUTODIFF_TAMC */
367 C-- Implicit Vertical Diffusion for Convection
368 c ==> should use sigmaR !!!
369 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
370 CALL CALC_IVDC(
371 I bi, bj, iMin, iMax, jMin, jMax, k,
372 I rhoKm1, rhoK,
373 U ConvectCount, KappaRT, KappaRS,
374 I myTime, myIter, myThid)
375 ENDIF
376
377 #endif /* SINGLE_LAYER_MODE */
378
379 C-- end of diagnostic k loop (Nr:1)
380 ENDDO
381
382 #ifdef ALLOW_AUTODIFF_TAMC
383 cph avoids recomputation of integrate_for_w
384 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
385 CADJ STORE pressure (:,:,:,bi,bj) =
386 CADJ & comlev1_bibj, key=itdkey, byte=isbyte
387 #endif /* ALLOW_AUTODIFF_TAMC */
388
389 #ifdef ALLOW_OBCS
390 C-- Calculate future values on open boundaries
391 IF (useOBCS) THEN
392 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
393 I uVel, vVel, wVel, theta, salt,
394 I myThid )
395 ENDIF
396 #endif /* ALLOW_OBCS */
397
398
399 c********************************************
400 cswdice --- add ---
401 #ifdef ALLOW_THERM_SEAICE
402 C-- Determines forcing terms based on external fields
403 c-- including effects from ice
404 CALL ICE_FORCING(
405 I bi, bj, iMin, iMax, jMin, jMax,
406 I myThid )
407 #else
408
409 cswdice --- end add ---
410
411 C-- Determines forcing terms based on external fields
412 C relaxation terms, etc.
413 CALL EXTERNAL_FORCING_SURF(
414 I bi, bj, iMin, iMax, jMin, jMax,
415 I myThid )
416 cswdice --- add ----
417 #endif
418 cswdice --- end add ---
419 c******************************************
420
421
422
423
424 #ifdef ALLOW_AUTODIFF_TAMC
425 cph needed for KPP
426 CADJ STORE surfacetendencyU(:,:,bi,bj)
427 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
428 CADJ STORE surfacetendencyV(:,:,bi,bj)
429 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
430 CADJ STORE surfacetendencyS(:,:,bi,bj)
431 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
432 CADJ STORE surfacetendencyT(:,:,bi,bj)
433 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
434 #endif /* ALLOW_AUTODIFF_TAMC */
435
436 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
437 C-- MOST of THERMODYNAMICS will be disabled
438 #ifndef SINGLE_LAYER_MODE
439
440 #ifdef ALLOW_GMREDI
441
442 #ifdef ALLOW_AUTODIFF_TAMC
443 cph storing here is needed only for one GMREDI_OPTIONS:
444 cph define GM_BOLUS_ADVEC
445 cph but I've avoided the #ifdef for now, in case more things change
446 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
447 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
448 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
449 #endif /* ALLOW_AUTODIFF_TAMC */
450
451 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
452 IF (useGMRedi) THEN
453 CALL GMREDI_CALC_TENSOR(
454 I bi, bj, iMin, iMax, jMin, jMax,
455 I sigmaX, sigmaY, sigmaR,
456 I myThid )
457 #ifdef ALLOW_AUTODIFF_TAMC
458 ELSE
459 CALL GMREDI_CALC_TENSOR_DUMMY(
460 I bi, bj, iMin, iMax, jMin, jMax,
461 I sigmaX, sigmaY, sigmaR,
462 I myThid )
463 #endif /* ALLOW_AUTODIFF_TAMC */
464 ENDIF
465
466 #ifdef ALLOW_AUTODIFF_TAMC
467 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
468 CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
469 CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
470 #endif /* ALLOW_AUTODIFF_TAMC */
471
472 #endif /* ALLOW_GMREDI */
473
474 #ifdef ALLOW_KPP
475 C-- Compute KPP mixing coefficients
476 IF (useKPP) THEN
477 CALL KPP_CALC(
478 I bi, bj, myTime, myThid )
479 #ifdef ALLOW_AUTODIFF_TAMC
480 ELSE
481 CALL KPP_CALC_DUMMY(
482 I bi, bj, myTime, myThid )
483 #endif /* ALLOW_AUTODIFF_TAMC */
484 ENDIF
485
486 #ifdef ALLOW_AUTODIFF_TAMC
487 CADJ STORE KPPghat (:,:,:,bi,bj)
488 CADJ & , KPPdiffKzT(:,:,:,bi,bj)
489 CADJ & , KPPdiffKzS(:,:,:,bi,bj)
490 CADJ & , KPPfrac (:,: ,bi,bj)
491 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
492 #endif /* ALLOW_AUTODIFF_TAMC */
493
494 #endif /* ALLOW_KPP */
495
496 #ifdef ALLOW_AUTODIFF_TAMC
497 CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
498 CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
499 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
500 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
501 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
502 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
503 #ifdef ALLOW_PASSIVE_TRACER
504 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
505 #endif
506 #endif /* ALLOW_AUTODIFF_TAMC */
507
508 #ifdef ALLOW_AIM
509 C AIM - atmospheric intermediate model, physics package code.
510 IF ( useAIM ) THEN
511 CALL TIMER_START('AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
512 CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
513 CALL TIMER_STOP( 'AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
514 ENDIF
515 #endif /* ALLOW_AIM */
516
517 #ifdef ALLOW_TIMEAVE
518 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
519 CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
520 I deltaTclock, bi, bj, myThid)
521 ENDIF
522 #endif /* ALLOW_TIMEAVE */
523
524 #ifndef DISABLE_MULTIDIM_ADVECTION
525 C-- Some advection schemes are better calculated using a multi-dimensional
526 C method in the absence of any other terms and, if used, is done here.
527 C
528 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
529 C The default is to use multi-dimensinal advection for non-linear advection
530 C schemes. However, for the sake of efficiency of the adjoint it is necessary
531 C to be able to exclude this scheme to avoid excessive storage and
532 C recomputation. It *is* differentiable, if you need it.
533 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
534 C disable this section of code.
535 IF (tempMultiDimAdvec) THEN
536 CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
537 U theta,gT,
538 I myTime,myIter,myThid)
539 ENDIF
540 IF (saltMultiDimAdvec) THEN
541 CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
542 U salt,gS,
543 I myTime,myIter,myThid)
544 ENDIF
545 C Since passive tracers are configurable separately from T,S we
546 C call the multi-dimensional method for PTRACERS regardless
547 C of whether multiDimAdvection is set or not.
548 #ifdef ALLOW_PTRACERS
549 IF ( usePTRACERS ) THEN
550 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
551 ENDIF
552 #endif /* ALLOW_PTRACERS */
553 #endif /* DISABLE_MULTIDIM_ADVECTION */
554
555 C-- Start of thermodynamics loop
556 DO k=Nr,1,-1
557 #ifdef ALLOW_AUTODIFF_TAMC
558 C? Patrick Is this formula correct?
559 cph Yes, but I rewrote it.
560 cph Also, the KappaR? need the index and subscript k!
561 kkey = (itdkey-1)*Nr + k
562 #endif /* ALLOW_AUTODIFF_TAMC */
563
564 C-- km1 Points to level above k (=k-1)
565 C-- kup Cycles through 1,2 to point to layer above
566 C-- kDown Cycles through 2,1 to point to current layer
567
568 km1 = MAX(1,k-1)
569 kup = 1+MOD(k+1,2)
570 kDown= 1+MOD(k,2)
571
572 iMin = 1-OLx
573 iMax = sNx+OLx
574 jMin = 1-OLy
575 jMax = sNy+OLy
576
577 C-- Get temporary terms used by tendency routines
578 CALL CALC_COMMON_FACTORS (
579 I bi,bj,iMin,iMax,jMin,jMax,k,
580 O xA,yA,uTrans,vTrans,rTrans,maskUp,
581 I myThid)
582
583 #ifdef ALLOW_GMREDI
584
585 C-- Residual transp = Bolus transp + Eulerian transp
586 IF (useGMRedi) THEN
587 CALL GMREDI_CALC_UVFLOW(
588 & uTrans, vTrans, bi, bj, k, myThid)
589 IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
590 & rTrans, bi, bj, k, myThid)
591 ENDIF
592
593 #ifdef ALLOW_AUTODIFF_TAMC
594 #ifdef GM_BOLUS_ADVEC
595 CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
596 CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
597 CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
598 #endif
599 #endif /* ALLOW_AUTODIFF_TAMC */
600
601 #endif /* ALLOW_GMREDI */
602
603 #ifdef ALLOW_AUTODIFF_TAMC
604 CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
605 CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
606 #endif /* ALLOW_AUTODIFF_TAMC */
607
608 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
609 C-- Calculate the total vertical diffusivity
610 CALL CALC_DIFFUSIVITY(
611 I bi,bj,iMin,iMax,jMin,jMax,k,
612 I maskUp,
613 O KappaRT,KappaRS,
614 I myThid)
615 #endif
616
617 iMin = 1-OLx+2
618 iMax = sNx+OLx-1
619 jMin = 1-OLy+2
620 jMax = sNy+OLy-1
621
622 C-- Calculate active tracer tendencies (gT,gS,...)
623 C and step forward storing result in gTnm1, gSnm1, etc.
624 IF ( tempStepping ) THEN
625 CALL CALC_GT(
626 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
627 I xA,yA,uTrans,vTrans,rTrans,maskUp,
628 I KappaRT,
629 U fVerT,
630 I myTime,myIter,myThid)
631 CALL TIMESTEP_TRACER(
632 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
633 I theta, gT,
634 I myIter, myThid)
635 ENDIF
636 cswdice ---- add ---
637 #ifdef ALLOW_THERM_SEAICE
638 if (k.eq.1) then
639 call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )
640 endif
641 #endif
642 cswdice -- end add ---
643 IF ( saltStepping ) THEN
644 CALL CALC_GS(
645 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
646 I xA,yA,uTrans,vTrans,rTrans,maskUp,
647 I KappaRS,
648 U fVerS,
649 I myTime,myIter,myThid)
650 CALL TIMESTEP_TRACER(
651 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
652 I salt, gS,
653 I myIter, myThid)
654 ENDIF
655 #ifdef ALLOW_PASSIVE_TRACER
656 IF ( tr1Stepping ) THEN
657 CALL CALC_GTR1(
658 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
659 I xA,yA,uTrans,vTrans,rTrans,maskUp,
660 I KappaRT,
661 U fVerTr1,
662 I myTime,myIter,myThid)
663 CALL TIMESTEP_TRACER(
664 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
665 I Tr1, gTr1,
666 I myIter,myThid)
667 ENDIF
668 #endif
669 #ifdef ALLOW_PTRACERS
670 IF ( usePTRACERS ) THEN
671 CALL PTRACERS_INTEGERATE(
672 I bi,bj,k,
673 I xA,yA,uTrans,vTrans,rTrans,maskUp,
674 X KappaRS,
675 I myIter,myTime,myThid)
676 ENDIF
677 #endif /* ALLOW_PTRACERS */
678
679 #ifdef ALLOW_OBCS
680 C-- Apply open boundary conditions
681 IF (useOBCS) THEN
682 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
683 END IF
684 #endif /* ALLOW_OBCS */
685
686 C-- Freeze water
687 IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN
688 #ifdef ALLOW_AUTODIFF_TAMC
689 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
690 CADJ & , key = kkey, byte = isbyte
691 #endif /* ALLOW_AUTODIFF_TAMC */
692 CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
693 END IF
694
695 C-- end of thermodynamic k loop (Nr:1)
696 ENDDO
697
698 cswdice -- add ---
699 #ifdef ALLOW_THERM_SEAICE
700 c timeaveraging for ice model values
701 CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
702 #endif
703 cswdice --- end add ---
704
705
706
707
708 C-- Implicit diffusion
709 IF (implicitDiffusion) THEN
710
711 IF (tempStepping) THEN
712 #ifdef ALLOW_AUTODIFF_TAMC
713 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
714 #endif /* ALLOW_AUTODIFF_TAMC */
715 CALL IMPLDIFF(
716 I bi, bj, iMin, iMax, jMin, jMax,
717 I deltaTtracer, KappaRT, recip_HFacC,
718 U gT,
719 I myThid )
720 ENDIF
721
722 IF (saltStepping) THEN
723 #ifdef ALLOW_AUTODIFF_TAMC
724 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
725 #endif /* ALLOW_AUTODIFF_TAMC */
726 CALL IMPLDIFF(
727 I bi, bj, iMin, iMax, jMin, jMax,
728 I deltaTtracer, KappaRS, recip_HFacC,
729 U gS,
730 I myThid )
731 ENDIF
732
733 #ifdef ALLOW_PASSIVE_TRACER
734 IF (tr1Stepping) THEN
735 #ifdef ALLOW_AUTODIFF_TAMC
736 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
737 #endif /* ALLOW_AUTODIFF_TAMC */
738 CALL IMPLDIFF(
739 I bi, bj, iMin, iMax, jMin, jMax,
740 I deltaTtracer, KappaRT, recip_HFacC,
741 U gTr1,
742 I myThid )
743 ENDIF
744 #endif
745
746 #ifdef ALLOW_PTRACERS
747 C Vertical diffusion (implicit) for passive tracers
748 IF ( usePTRACERS ) THEN
749 CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
750 ENDIF
751 #endif /* ALLOW_PTRACERS */
752
753 #ifdef ALLOW_OBCS
754 C-- Apply open boundary conditions
755 IF (useOBCS) THEN
756 DO K=1,Nr
757 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
758 ENDDO
759 END IF
760 #endif /* ALLOW_OBCS */
761
762 C-- End If implicitDiffusion
763 ENDIF
764
765 #endif /* SINGLE_LAYER_MODE */
766
767 Ccs-
768 ENDDO
769 ENDDO
770
771 #ifdef ALLOW_AIM
772 c IF ( useAIM ) THEN
773 c CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
774 c ENDIF
775 #endif /* ALLOW_AIM */
776 c IF ( staggerTimeStep ) THEN
777 c IF ( useAIM .OR. useCubedSphereExchange ) THEN
778 c IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
779 c IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
780 c ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
781 cc .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
782 c IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
783 c IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
784 c ENDIF
785 c ENDIF
786
787 #ifndef DISABLE_DEBUGMODE
788 If (debugMode) THEN
789 CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
790 CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
791 CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
792 CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
793 CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
794 CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
795 CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
796 CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
797 CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
798 #ifdef ALLOW_PTRACERS
799 IF ( usePTRACERS ) THEN
800 CALL PTRACERS_DEBUG(myThid)
801 ENDIF
802 #endif /* ALLOW_PTRACERS */
803 ENDIF
804 #endif
805
806 RETURN
807 END

  ViewVC Help
Powered by ViewVC 1.1.22