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