/[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.59 - (show annotations) (download)
Thu Nov 13 21:54:11 2003 UTC (20 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint52a_post
Changes since 1.58: +12 -12 lines 
additional changes for FREEZE:
 - new S/R FREEZE_SURFACE only apllied to surface level.
 - add run-time parameter "useOldFreezing" to use the old version "FREEZE"
1 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.58 2003/11/13 06:35:14 dimitri 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 #ifdef ALLOW_DEBUG
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 #ifdef ALLOW_DEBUG
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 #ifdef ALLOW_DEBUG
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 #ifdef ALLOW_DEBUG
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 #ifdef ALLOW_DEBUG
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 #ifdef ALLOW_DEBUG
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 #ifdef ALLOW_DEBUG
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 #ifdef ALLOW_DEBUG
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 # ifdef ALLOW_SEAICE
475 CADJ STORE surfacetendencyTice(:,:,bi,bj)
476 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
477 # endif
478 #endif /* ALLOW_AUTODIFF_TAMC */
479
480 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
481 C-- MOST of THERMODYNAMICS will be disabled
482 #ifndef SINGLE_LAYER_MODE
483
484 #ifdef ALLOW_GMREDI
485
486 #ifdef ALLOW_AUTODIFF_TAMC
487 cph storing here is needed only for one GMREDI_OPTIONS:
488 cph define GM_BOLUS_ADVEC
489 cph but I've avoided the #ifdef for now, in case more things change
490 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
491 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
492 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
493 #endif /* ALLOW_AUTODIFF_TAMC */
494
495 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
496 IF (useGMRedi) THEN
497 #ifdef ALLOW_DEBUG
498 IF ( debugLevel .GE. debLevB )
499 & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
500 #endif
501 CALL GMREDI_CALC_TENSOR(
502 I bi, bj, iMin, iMax, jMin, jMax,
503 I sigmaX, sigmaY, sigmaR,
504 I myThid )
505 #ifdef ALLOW_AUTODIFF_TAMC
506 ELSE
507 CALL GMREDI_CALC_TENSOR_DUMMY(
508 I bi, bj, iMin, iMax, jMin, jMax,
509 I sigmaX, sigmaY, sigmaR,
510 I myThid )
511 #endif /* ALLOW_AUTODIFF_TAMC */
512 ENDIF
513
514 #ifdef ALLOW_AUTODIFF_TAMC
515 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
516 CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
517 CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
518 #endif /* ALLOW_AUTODIFF_TAMC */
519
520 #endif /* ALLOW_GMREDI */
521
522 #ifdef ALLOW_KPP
523 C-- Compute KPP mixing coefficients
524 IF (useKPP) THEN
525 #ifdef ALLOW_DEBUG
526 IF ( debugLevel .GE. debLevB )
527 & CALL DEBUG_CALL('KPP_CALC',myThid)
528 #endif
529 CALL KPP_CALC(
530 I bi, bj, myTime, myThid )
531 #ifdef ALLOW_AUTODIFF_TAMC
532 ELSE
533 CALL KPP_CALC_DUMMY(
534 I bi, bj, myTime, myThid )
535 #endif /* ALLOW_AUTODIFF_TAMC */
536 ENDIF
537
538 #ifdef ALLOW_AUTODIFF_TAMC
539 CADJ STORE KPPghat (:,:,:,bi,bj)
540 CADJ & , KPPdiffKzT(:,:,:,bi,bj)
541 CADJ & , KPPdiffKzS(:,:,:,bi,bj)
542 CADJ & , KPPfrac (:,: ,bi,bj)
543 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
544 #endif /* ALLOW_AUTODIFF_TAMC */
545
546 #endif /* ALLOW_KPP */
547
548 #ifdef ALLOW_AUTODIFF_TAMC
549 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
550 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
551 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
552 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
553 #ifdef ALLOW_PASSIVE_TRACER
554 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
555 #endif
556 #ifdef ALLOW_PTRACERS
557 cph-- moved to forward_step to avoid key computation
558 cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
559 cphCADJ & key=itdkey, byte=isbyte
560 #endif
561 #endif /* ALLOW_AUTODIFF_TAMC */
562
563 #ifdef ALLOW_AIM
564 C AIM - atmospheric intermediate model, physics package code.
565 IF ( useAIM ) THEN
566 #ifdef ALLOW_DEBUG
567 IF ( debugLevel .GE. debLevB )
568 & CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
569 #endif
570 CALL TIMER_START('AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
571 CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
572 CALL TIMER_STOP( 'AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
573 ENDIF
574 #endif /* ALLOW_AIM */
575
576 #ifndef DISABLE_MULTIDIM_ADVECTION
577 C-- Some advection schemes are better calculated using a multi-dimensional
578 C method in the absence of any other terms and, if used, is done here.
579 C
580 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
581 C The default is to use multi-dimensinal advection for non-linear advection
582 C schemes. However, for the sake of efficiency of the adjoint it is necessary
583 C to be able to exclude this scheme to avoid excessive storage and
584 C recomputation. It *is* differentiable, if you need it.
585 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
586 C disable this section of code.
587 IF (tempMultiDimAdvec) THEN
588 #ifdef ALLOW_DEBUG
589 IF ( debugLevel .GE. debLevB )
590 & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
591 #endif
592 CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
593 U theta,gT,
594 I myTime,myIter,myThid)
595 ENDIF
596 IF (saltMultiDimAdvec) THEN
597 #ifdef ALLOW_DEBUG
598 IF ( debugLevel .GE. debLevB )
599 & CALL DEBUG_CALL('GAD_ADVECTION',myThid)
600 #endif
601 CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
602 U salt,gS,
603 I myTime,myIter,myThid)
604 ENDIF
605 C Since passive tracers are configurable separately from T,S we
606 C call the multi-dimensional method for PTRACERS regardless
607 C of whether multiDimAdvection is set or not.
608 #ifdef ALLOW_PTRACERS
609 IF ( usePTRACERS ) THEN
610 #ifdef ALLOW_DEBUG
611 IF ( debugLevel .GE. debLevB )
612 & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
613 #endif
614 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
615 ENDIF
616 #endif /* ALLOW_PTRACERS */
617 #endif /* DISABLE_MULTIDIM_ADVECTION */
618
619 #ifdef ALLOW_DEBUG
620 IF ( debugLevel .GE. debLevB )
621 & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
622 #endif
623
624 C-- Start of thermodynamics loop
625 DO k=Nr,1,-1
626 #ifdef ALLOW_AUTODIFF_TAMC
627 C? Patrick Is this formula correct?
628 cph Yes, but I rewrote it.
629 cph Also, the KappaR? need the index and subscript k!
630 kkey = (itdkey-1)*Nr + k
631 #endif /* ALLOW_AUTODIFF_TAMC */
632
633 C-- km1 Points to level above k (=k-1)
634 C-- kup Cycles through 1,2 to point to layer above
635 C-- kDown Cycles through 2,1 to point to current layer
636
637 km1 = MAX(1,k-1)
638 kup = 1+MOD(k+1,2)
639 kDown= 1+MOD(k,2)
640
641 iMin = 1-OLx
642 iMax = sNx+OLx
643 jMin = 1-OLy
644 jMax = sNy+OLy
645
646 C-- Get temporary terms used by tendency routines
647 CALL CALC_COMMON_FACTORS (
648 I bi,bj,iMin,iMax,jMin,jMax,k,
649 O xA,yA,uTrans,vTrans,rTrans,maskUp,
650 I myThid)
651
652 #ifdef ALLOW_GMREDI
653
654 C-- Residual transp = Bolus transp + Eulerian transp
655 IF (useGMRedi) THEN
656 CALL GMREDI_CALC_UVFLOW(
657 & uTrans, vTrans, bi, bj, k, myThid)
658 IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
659 & rTrans, bi, bj, k, myThid)
660 ENDIF
661
662 #ifdef ALLOW_AUTODIFF_TAMC
663 #ifdef GM_BOLUS_ADVEC
664 CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
665 CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
666 CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
667 #endif
668 #endif /* ALLOW_AUTODIFF_TAMC */
669
670 #endif /* ALLOW_GMREDI */
671
672 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
673 C-- Calculate the total vertical diffusivity
674 CALL CALC_DIFFUSIVITY(
675 I bi,bj,iMin,iMax,jMin,jMax,k,
676 I maskUp,
677 O KappaRT,KappaRS,
678 I myThid)
679 # ifdef ALLOW_AUTODIFF_TAMC
680 CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
681 CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
682 # endif /* ALLOW_AUTODIFF_TAMC */
683 #endif
684
685 iMin = 1-OLx+2
686 iMax = sNx+OLx-1
687 jMin = 1-OLy+2
688 jMax = sNy+OLy-1
689
690 C-- Calculate active tracer tendencies (gT,gS,...)
691 C and step forward storing result in gTnm1, gSnm1, etc.
692 IF ( tempStepping ) THEN
693 CALL CALC_GT(
694 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
695 I xA,yA,uTrans,vTrans,rTrans,maskUp,
696 I KappaRT,
697 U fVerT,
698 I myTime,myIter,myThid)
699 CALL TIMESTEP_TRACER(
700 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
701 I theta, gT,
702 I myIter, myThid)
703 ENDIF
704
705 #ifdef ALLOW_THERM_SEAICE
706 IF (useThermSeaIce .AND. k.EQ.1) THEN
707 CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )
708 ENDIF
709 #endif
710
711 IF ( saltStepping ) THEN
712 CALL CALC_GS(
713 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
714 I xA,yA,uTrans,vTrans,rTrans,maskUp,
715 I KappaRS,
716 U fVerS,
717 I myTime,myIter,myThid)
718 CALL TIMESTEP_TRACER(
719 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
720 I salt, gS,
721 I myIter, myThid)
722 ENDIF
723 #ifdef ALLOW_PASSIVE_TRACER
724 ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
725 IF ( tr1Stepping ) THEN
726 CALL CALC_GTR1(
727 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
728 I xA,yA,uTrans,vTrans,rTrans,maskUp,
729 I KappaRT,
730 U fVerTr1,
731 I myTime,myIter,myThid)
732 CALL TIMESTEP_TRACER(
733 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
734 I Tr1, gTr1,
735 I myIter,myThid)
736 ENDIF
737 #endif
738 #ifdef ALLOW_PTRACERS
739 IF ( usePTRACERS ) THEN
740 CALL PTRACERS_INTEGRATE(
741 I bi,bj,k,
742 I xA,yA,uTrans,vTrans,rTrans,maskUp,
743 X fVerP, KappaRS,
744 I myIter,myTime,myThid)
745 ENDIF
746 #endif /* ALLOW_PTRACERS */
747
748 #ifdef ALLOW_OBCS
749 C-- Apply open boundary conditions
750 IF (useOBCS) THEN
751 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
752 END IF
753 #endif /* ALLOW_OBCS */
754
755 C-- Freeze water
756 C this bit of code is left here for backward compatibility.
757 C freezing at surface level has been moved to FORWARD_STEP
758 IF ( useOldFreezing .AND. .NOT. useSEAICE
759 & .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN
760 #ifdef ALLOW_AUTODIFF_TAMC
761 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
762 CADJ & , key = kkey, byte = isbyte
763 #endif /* ALLOW_AUTODIFF_TAMC */
764 CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
765 ENDIF
766
767 C-- end of thermodynamic k loop (Nr:1)
768 ENDDO
769
770 cswdice -- add ---
771 #ifdef ALLOW_THERM_SEAICE
772 c timeaveraging for ice model values
773 ceh3 This should be wrapped in an IF ( useThermSeaIce ) THEN
774 CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
775 #endif
776 cswdice --- end add ---
777
778
779
780
781 C-- Implicit diffusion
782 IF (implicitDiffusion) THEN
783
784 IF (tempStepping) THEN
785 #ifdef ALLOW_AUTODIFF_TAMC
786 CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
787 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
788 #endif /* ALLOW_AUTODIFF_TAMC */
789 CALL IMPLDIFF(
790 I bi, bj, iMin, iMax, jMin, jMax,
791 I deltaTtracer, KappaRT, recip_HFacC,
792 U gT,
793 I myThid )
794 ENDIF
795
796 IF (saltStepping) THEN
797 #ifdef ALLOW_AUTODIFF_TAMC
798 CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
799 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
800 #endif /* ALLOW_AUTODIFF_TAMC */
801 CALL IMPLDIFF(
802 I bi, bj, iMin, iMax, jMin, jMax,
803 I deltaTtracer, KappaRS, recip_HFacC,
804 U gS,
805 I myThid )
806 ENDIF
807
808 #ifdef ALLOW_PASSIVE_TRACER
809 IF (tr1Stepping) THEN
810 #ifdef ALLOW_AUTODIFF_TAMC
811 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
812 #endif /* ALLOW_AUTODIFF_TAMC */
813 CALL IMPLDIFF(
814 I bi, bj, iMin, iMax, jMin, jMax,
815 I deltaTtracer, KappaRT, recip_HFacC,
816 U gTr1,
817 I myThid )
818 ENDIF
819 #endif
820
821 #ifdef ALLOW_PTRACERS
822 C Vertical diffusion (implicit) for passive tracers
823 IF ( usePTRACERS ) THEN
824 CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
825 ENDIF
826 #endif /* ALLOW_PTRACERS */
827
828 #ifdef ALLOW_OBCS
829 C-- Apply open boundary conditions
830 IF (useOBCS) THEN
831 DO K=1,Nr
832 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
833 ENDDO
834 END IF
835 #endif /* ALLOW_OBCS */
836
837 C-- End If implicitDiffusion
838 ENDIF
839
840 #ifdef ALLOW_TIMEAVE
841 ceh3 needs an IF ( useTIMEAVE ) THEN
842 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
843 CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
844 I Nr, deltaTclock, bi, bj, myThid)
845 ENDIF
846 useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
847 IF (taveFreq.GT.0. .AND. useVariableK ) THEN
848 IF (implicitDiffusion) THEN
849 CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
850 I Nr, 3, deltaTclock, bi, bj, myThid)
851 ELSE
852 CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
853 I Nr, 3, deltaTclock, bi, bj, myThid)
854 ENDIF
855 ENDIF
856 #endif /* ALLOW_TIMEAVE */
857
858 #endif /* SINGLE_LAYER_MODE */
859
860 C-- end bi,bj loops.
861 ENDDO
862 ENDDO
863
864 #ifdef ALLOW_DEBUG
865 If (debugMode) THEN
866 CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
867 CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
868 CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
869 CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
870 CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
871 CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
872 CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
873 CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
874 CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
875 #ifdef ALLOW_PTRACERS
876 IF ( usePTRACERS ) THEN
877 CALL PTRACERS_DEBUG(myThid)
878 ENDIF
879 #endif /* ALLOW_PTRACERS */
880 ENDIF
881 #endif
882
883 #ifdef ALLOW_DEBUG
884 IF ( debugLevel .GE. debLevB )
885 & CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
886 #endif
887
888 RETURN
889 END
  ViewVC Help
Powered by ViewVC 1.1.22