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