/[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.18 - (show annotations) (download)
Tue Mar 5 14:15:34 2002 UTC (22 years, 3 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint44f_post, checkpoint44g_post
Changes since 1.17: +6 -2 lines
Bug fix: Missing #ifdef's for when PTRACERS are not enabled.

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/thermodynamics.F,v 1.17 2002/03/04 17:26:41 adcroft Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: THERMODYNAMICS
8 C !INTERFACE:
9 SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
10 C !DESCRIPTION: \bv
11 C *==========================================================*
12 C | SUBROUTINE THERMODYNAMICS
13 C | o Controlling routine for the prognostic part of the
14 C | thermo-dynamics.
15 C *===========================================================
16 C | The algorithm...
17 C |
18 C | "Correction Step"
19 C | =================
20 C | Here we update the horizontal velocities with the surface
21 C | pressure such that the resulting flow is either consistent
22 C | with the free-surface evolution or the rigid-lid:
23 C | U[n] = U* + dt x d/dx P
24 C | V[n] = V* + dt x d/dy P
25 C |
26 C | "Calculation of Gs"
27 C | ===================
28 C | This is where all the accelerations and tendencies (ie.
29 C | physics, parameterizations etc...) are calculated
30 C | rho = rho ( theta[n], salt[n] )
31 C | b = b(rho, theta)
32 C | K31 = K31 ( rho )
33 C | Gu[n] = Gu( u[n], v[n], wVel, b, ... )
34 C | Gv[n] = Gv( u[n], v[n], wVel, b, ... )
35 C | Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
36 C | Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
37 C |
38 C | "Time-stepping" or "Prediction"
39 C | ================================
40 C | The models variables are stepped forward with the appropriate
41 C | time-stepping scheme (currently we use Adams-Bashforth II)
42 C | - For momentum, the result is always *only* a "prediction"
43 C | in that the flow may be divergent and will be "corrected"
44 C | later with a surface pressure gradient.
45 C | - Normally for tracers the result is the new field at time
46 C | level [n+1} *BUT* in the case of implicit diffusion the result
47 C | is also *only* a prediction.
48 C | - We denote "predictors" with an asterisk (*).
49 C | U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
50 C | V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
51 C | theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
52 C | salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
53 C | With implicit diffusion:
54 C | theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
55 C | salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
56 C | (1 + dt * K * d_zz) theta[n] = theta*
57 C | (1 + dt * K * d_zz) salt[n] = salt*
58 C |
59 C *==========================================================*
60 C \ev
61
62 C !USES:
63 IMPLICIT NONE
64 C == Global variables ===
65 #include "SIZE.h"
66 #include "EEPARAMS.h"
67 #include "PARAMS.h"
68 #include "DYNVARS.h"
69 #include "GRID.h"
70 #include "GAD.h"
71 #ifdef ALLOW_PASSIVE_TRACER
72 #include "TR1.h"
73 #endif
74 #ifdef ALLOW_AUTODIFF_TAMC
75 # include "tamc.h"
76 # include "tamc_keys.h"
77 # include "FFIELDS.h"
78 # ifdef ALLOW_KPP
79 # include "KPP.h"
80 # endif
81 # ifdef ALLOW_GMREDI
82 # include "GMREDI.h"
83 # endif
84 #endif /* ALLOW_AUTODIFF_TAMC */
85 #ifdef ALLOW_TIMEAVE
86 #include "TIMEAVE_STATV.h"
87 #endif
88
89 C !INPUT/OUTPUT PARAMETERS:
90 C == Routine arguments ==
91 C myTime - Current time in simulation
92 C myIter - Current iteration number in simulation
93 C myThid - Thread number for this instance of the routine.
94 _RL myTime
95 INTEGER myIter
96 INTEGER myThid
97
98 C !LOCAL VARIABLES:
99 C == Local variables
100 C xA, yA - Per block temporaries holding face areas
101 C uTrans, vTrans, rTrans - Per block temporaries holding flow
102 C transport
103 C o uTrans: Zonal transport
104 C o vTrans: Meridional transport
105 C o rTrans: Vertical transport
106 C maskUp o maskUp: land/water mask for W points
107 C fVer[STUV] o fVer: Vertical flux term - note fVer
108 C is "pipelined" in the vertical
109 C so we need an fVer for each
110 C variable.
111 C rhoK, rhoKM1 - Density at current level, and level above
112 C phiHyd - Hydrostatic part of the potential phiHydi.
113 C In z coords phiHydiHyd is the hydrostatic
114 C Potential (=pressure/rho0) anomaly
115 C In p coords phiHydiHyd is the geopotential
116 C surface height anomaly.
117 C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
118 C phiSurfY or geopotentiel (atmos) in X and Y direction
119 C KappaRT, - Total diffusion in vertical for T and S.
120 C KappaRS (background + spatially varying, isopycnal term).
121 C iMin, iMax - Ranges and sub-block indices on which calculations
122 C jMin, jMax are applied.
123 C bi, bj
124 C k, kup, - Index for layer above and below. kup and kDown
125 C kDown, km1 are switched with layer to be the appropriate
126 C index into fVerTerm.
127 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
134 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
135 _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
136 _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
137 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140 _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141 _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
142 _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
143 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
144 _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
145 _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
146 C This is currently used by IVDC and Diagnostics
147 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
148 INTEGER iMin, iMax
149 INTEGER jMin, jMax
150 INTEGER bi, bj
151 INTEGER i, j
152 INTEGER k, km1, kup, kDown
153
154 CEOP
155
156 #ifdef ALLOW_AUTODIFF_TAMC
157 C-- dummy statement to end declaration part
158 ikey = 1
159 #endif /* ALLOW_AUTODIFF_TAMC */
160
161 C-- Set up work arrays with valid (i.e. not NaN) values
162 C These inital values do not alter the numerical results. They
163 C just ensure that all memory references are to valid floating
164 C point numbers. This prevents spurious hardware signals due to
165 C uninitialised but inert locations.
166 DO j=1-OLy,sNy+OLy
167 DO i=1-OLx,sNx+OLx
168 xA(i,j) = 0. _d 0
169 yA(i,j) = 0. _d 0
170 uTrans(i,j) = 0. _d 0
171 vTrans(i,j) = 0. _d 0
172 DO k=1,Nr
173 phiHyd(i,j,k) = 0. _d 0
174 sigmaX(i,j,k) = 0. _d 0
175 sigmaY(i,j,k) = 0. _d 0
176 sigmaR(i,j,k) = 0. _d 0
177 ENDDO
178 rhoKM1 (i,j) = 0. _d 0
179 rhok (i,j) = 0. _d 0
180 phiSurfX(i,j) = 0. _d 0
181 phiSurfY(i,j) = 0. _d 0
182 ENDDO
183 ENDDO
184
185
186 #ifdef ALLOW_AUTODIFF_TAMC
187 C-- HPF directive to help TAMC
188 CHPF$ INDEPENDENT
189 #endif /* ALLOW_AUTODIFF_TAMC */
190
191 DO bj=myByLo(myThid),myByHi(myThid)
192
193 #ifdef ALLOW_AUTODIFF_TAMC
194 C-- HPF directive to help TAMC
195 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
196 CHPF$& ,phiHyd,utrans,vtrans,xA,yA
197 CHPF$& ,KappaRT,KappaRS
198 CHPF$& )
199 #endif /* ALLOW_AUTODIFF_TAMC */
200
201 DO bi=myBxLo(myThid),myBxHi(myThid)
202
203 #ifdef ALLOW_AUTODIFF_TAMC
204 act1 = bi - myBxLo(myThid)
205 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
206 act2 = bj - myByLo(myThid)
207 max2 = myByHi(myThid) - myByLo(myThid) + 1
208 act3 = myThid - 1
209 max3 = nTx*nTy
210 act4 = ikey_dynamics - 1
211 ikey = (act1 + 1) + act2*max1
212 & + act3*max1*max2
213 & + act4*max1*max2*max3
214 #endif /* ALLOW_AUTODIFF_TAMC */
215
216 C-- Set up work arrays that need valid initial values
217 DO j=1-OLy,sNy+OLy
218 DO i=1-OLx,sNx+OLx
219 rTrans (i,j) = 0. _d 0
220 fVerT (i,j,1) = 0. _d 0
221 fVerT (i,j,2) = 0. _d 0
222 fVerS (i,j,1) = 0. _d 0
223 fVerS (i,j,2) = 0. _d 0
224 fVerTr1(i,j,1) = 0. _d 0
225 fVerTr1(i,j,2) = 0. _d 0
226 ENDDO
227 ENDDO
228
229 DO k=1,Nr
230 DO j=1-OLy,sNy+OLy
231 DO i=1-OLx,sNx+OLx
232 C This is currently also used by IVDC and Diagnostics
233 ConvectCount(i,j,k) = 0.
234 KappaRT(i,j,k) = 0. _d 0
235 KappaRS(i,j,k) = 0. _d 0
236 #ifdef ALLOW_AUTODIFF_TAMC
237 gT(i,j,k,bi,bj) = 0. _d 0
238 gS(i,j,k,bi,bj) = 0. _d 0
239 #ifdef ALLOW_PASSIVE_TRACER
240 gTr1(i,j,k,bi,bj) = 0. _d 0
241 #endif
242 #endif
243 ENDDO
244 ENDDO
245 ENDDO
246
247 iMin = 1-OLx+1
248 iMax = sNx+OLx
249 jMin = 1-OLy+1
250 jMax = sNy+OLy
251
252
253 #ifdef ALLOW_AUTODIFF_TAMC
254 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
255 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
256 #ifdef ALLOW_KPP
257 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
258 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
259 #endif
260 #endif /* ALLOW_AUTODIFF_TAMC */
261
262 C-- Start of diagnostic loop
263 DO k=Nr,1,-1
264
265 #ifdef ALLOW_AUTODIFF_TAMC
266 C? Patrick, is this formula correct now that we change the loop range?
267 C? Do we still need this?
268 cph kkey formula corrected.
269 cph Needed for rhok, rhokm1, in the case useGMREDI.
270 kkey = (ikey-1)*Nr + k
271 CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
272 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
273 #endif /* ALLOW_AUTODIFF_TAMC */
274
275 C-- Integrate continuity vertically for vertical velocity
276 CALL INTEGRATE_FOR_W(
277 I bi, bj, k, uVel, vVel,
278 O wVel,
279 I myThid )
280
281 #ifdef ALLOW_OBCS
282 #ifdef ALLOW_NONHYDROSTATIC
283 C-- Apply OBC to W if in N-H mode
284 IF (useOBCS.AND.nonHydrostatic) THEN
285 CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
286 ENDIF
287 #endif /* ALLOW_NONHYDROSTATIC */
288 #endif /* ALLOW_OBCS */
289
290 C-- Calculate gradients of potential density for isoneutral
291 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
292 c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
293 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
294 #ifdef ALLOW_AUTODIFF_TAMC
295 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
296 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
297 #endif /* ALLOW_AUTODIFF_TAMC */
298 CALL FIND_RHO(
299 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
300 I theta, salt,
301 O rhoK,
302 I myThid )
303 IF (k.GT.1) THEN
304 #ifdef ALLOW_AUTODIFF_TAMC
305 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
306 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
307 #endif /* ALLOW_AUTODIFF_TAMC */
308 CALL FIND_RHO(
309 I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
310 I theta, salt,
311 O rhoKm1,
312 I myThid )
313 ENDIF
314 CALL GRAD_SIGMA(
315 I bi, bj, iMin, iMax, jMin, jMax, k,
316 I rhoK, rhoKm1, rhoK,
317 O sigmaX, sigmaY, sigmaR,
318 I myThid )
319 ENDIF
320
321 C-- Implicit Vertical Diffusion for Convection
322 c ==> should use sigmaR !!!
323 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
324 CALL CALC_IVDC(
325 I bi, bj, iMin, iMax, jMin, jMax, k,
326 I rhoKm1, rhoK,
327 U ConvectCount, KappaRT, KappaRS,
328 I myTime, myIter, myThid)
329 ENDIF
330
331 C-- end of diagnostic k loop (Nr:1)
332 ENDDO
333
334 #ifdef ALLOW_AUTODIFF_TAMC
335 cph avoids recomputation of integrate_for_w
336 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
337 #endif /* ALLOW_AUTODIFF_TAMC */
338
339 #ifdef ALLOW_OBCS
340 C-- Calculate future values on open boundaries
341 IF (useOBCS) THEN
342 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
343 I uVel, vVel, wVel, theta, salt,
344 I myThid )
345 ENDIF
346 #endif /* ALLOW_OBCS */
347
348 C-- Determines forcing terms based on external fields
349 C relaxation terms, etc.
350 CALL EXTERNAL_FORCING_SURF(
351 I bi, bj, iMin, iMax, jMin, jMax,
352 I myThid )
353 #ifdef ALLOW_AUTODIFF_TAMC
354 cph needed for KPP
355 CADJ STORE surfacetendencyU(:,:,bi,bj)
356 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
357 CADJ STORE surfacetendencyV(:,:,bi,bj)
358 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
359 CADJ STORE surfacetendencyS(:,:,bi,bj)
360 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
361 CADJ STORE surfacetendencyT(:,:,bi,bj)
362 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
363 #endif /* ALLOW_AUTODIFF_TAMC */
364
365 #ifdef ALLOW_GMREDI
366
367 #ifdef ALLOW_AUTODIFF_TAMC
368 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
369 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
370 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
371 #endif /* ALLOW_AUTODIFF_TAMC */
372 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
373 IF (useGMRedi) THEN
374 CALL GMREDI_CALC_TENSOR(
375 I bi, bj, iMin, iMax, jMin, jMax,
376 I sigmaX, sigmaY, sigmaR,
377 I myThid )
378 #ifdef ALLOW_AUTODIFF_TAMC
379 ELSE
380 CALL GMREDI_CALC_TENSOR_DUMMY(
381 I bi, bj, iMin, iMax, jMin, jMax,
382 I sigmaX, sigmaY, sigmaR,
383 I myThid )
384 #endif /* ALLOW_AUTODIFF_TAMC */
385 ENDIF
386
387 #ifdef ALLOW_AUTODIFF_TAMC
388 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
389 CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
390 CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
391 #endif /* ALLOW_AUTODIFF_TAMC */
392
393 #endif /* ALLOW_GMREDI */
394
395 #ifdef ALLOW_KPP
396 C-- Compute KPP mixing coefficients
397 IF (useKPP) THEN
398 CALL KPP_CALC(
399 I bi, bj, myTime, myThid )
400 #ifdef ALLOW_AUTODIFF_TAMC
401 ELSE
402 CALL KPP_CALC_DUMMY(
403 I bi, bj, myTime, myThid )
404 #endif /* ALLOW_AUTODIFF_TAMC */
405 ENDIF
406
407 #ifdef ALLOW_AUTODIFF_TAMC
408 CADJ STORE KPPghat (:,:,:,bi,bj)
409 CADJ & , KPPdiffKzT(:,:,:,bi,bj)
410 CADJ & , KPPdiffKzS(:,:,:,bi,bj)
411 CADJ & , KPPfrac (:,: ,bi,bj)
412 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
413 #endif /* ALLOW_AUTODIFF_TAMC */
414
415 #endif /* ALLOW_KPP */
416
417 #ifdef ALLOW_AUTODIFF_TAMC
418 CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
419 CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
420 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
421 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
422 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
423 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
424 #ifdef ALLOW_PASSIVE_TRACER
425 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
426 #endif
427 #endif /* ALLOW_AUTODIFF_TAMC */
428
429 #ifdef ALLOW_AIM
430 C AIM - atmospheric intermediate model, physics package code.
431 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
432 IF ( useAIM ) THEN
433 CALL TIMER_START('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
434 CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )
435 CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
436 ENDIF
437 #endif /* ALLOW_AIM */
438
439 #ifdef ALLOW_TIMEAVE
440 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
441 CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
442 I deltaTclock, bi, bj, myThid)
443 ENDIF
444 #endif /* ALLOW_TIMEAVE */
445
446 #ifndef DISABLE_MULTIDIM_ADVECTION
447 C-- Some advection schemes are better calculated using a multi-dimensional
448 C method in the absence of any other terms and, if used, is done here.
449 C
450 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
451 C The default is to use multi-dimensinal advection for non-linear advection
452 C schemes. However, for the sake of efficiency of the adjoint it is necessary
453 C to be able to exclude this scheme to avoid excessive storage and
454 C recomputation. It *is* differentiable, if you need it.
455 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
456 C disable this section of code.
457 IF (multiDimAdvection) THEN
458 IF (tempStepping .AND.
459 & tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.
460 & tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.
461 & tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN
462 CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
463 U theta,gT,
464 I myTime,myIter,myThid)
465 ENDIF
466 IF (saltStepping .AND.
467 & saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.
468 & saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.
469 & saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN
470 CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
471 U salt,gS,
472 I myTime,myIter,myThid)
473 ENDIF
474 ENDIF
475 C Since passive tracers are configurable separately from T,S we
476 C call the multi-dimensional method for PTRACERS regardless
477 C of whether multiDimAdvection is set or not.
478 #ifdef ALLOW_PTRACERS
479 IF ( usePTRACERS ) THEN
480 CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
481 ENDIF
482 #endif /* ALLOW_PTRACERS */
483 #endif /* DISABLE_MULTIDIM_ADVECTION */
484
485 C-- Start of thermodynamics loop
486 DO k=Nr,1,-1
487 #ifdef ALLOW_AUTODIFF_TAMC
488 C? Patrick Is this formula correct?
489 cph Yes, but I rewrote it.
490 cph Also, the KappaR? need the index and subscript k!
491 kkey = (ikey-1)*Nr + k
492 #endif /* ALLOW_AUTODIFF_TAMC */
493
494 C-- km1 Points to level above k (=k-1)
495 C-- kup Cycles through 1,2 to point to layer above
496 C-- kDown Cycles through 2,1 to point to current layer
497
498 km1 = MAX(1,k-1)
499 kup = 1+MOD(k+1,2)
500 kDown= 1+MOD(k,2)
501
502 iMin = 1-OLx
503 iMax = sNx+OLx
504 jMin = 1-OLy
505 jMax = sNy+OLy
506
507 C-- Get temporary terms used by tendency routines
508 CALL CALC_COMMON_FACTORS (
509 I bi,bj,iMin,iMax,jMin,jMax,k,
510 O xA,yA,uTrans,vTrans,rTrans,maskUp,
511 I myThid)
512
513 #ifdef ALLOW_AUTODIFF_TAMC
514 CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
515 CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
516 #endif /* ALLOW_AUTODIFF_TAMC */
517
518 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
519 C-- Calculate the total vertical diffusivity
520 CALL CALC_DIFFUSIVITY(
521 I bi,bj,iMin,iMax,jMin,jMax,k,
522 I maskUp,
523 O KappaRT,KappaRS,
524 I myThid)
525 #endif
526
527 iMin = 1-OLx+2
528 iMax = sNx+OLx-1
529 jMin = 1-OLy+2
530 jMax = sNy+OLy-1
531
532 C-- Calculate active tracer tendencies (gT,gS,...)
533 C and step forward storing result in gTnm1, gSnm1, etc.
534 IF ( tempStepping ) THEN
535 CALL CALC_GT(
536 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
537 I xA,yA,uTrans,vTrans,rTrans,maskUp,
538 I KappaRT,
539 U fVerT,
540 I myTime,myIter,myThid)
541 CALL TIMESTEP_TRACER(
542 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
543 I theta, gT,
544 I myIter, myThid)
545 ENDIF
546 IF ( saltStepping ) THEN
547 CALL CALC_GS(
548 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
549 I xA,yA,uTrans,vTrans,rTrans,maskUp,
550 I KappaRS,
551 U fVerS,
552 I myTime,myIter,myThid)
553 CALL TIMESTEP_TRACER(
554 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
555 I salt, gS,
556 I myIter, myThid)
557 ENDIF
558 #ifdef ALLOW_PASSIVE_TRACER
559 IF ( tr1Stepping ) THEN
560 CALL CALC_GTR1(
561 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
562 I xA,yA,uTrans,vTrans,rTrans,maskUp,
563 I KappaRT,
564 U fVerTr1,
565 I myTime,myIter,myThid)
566 CALL TIMESTEP_TRACER(
567 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
568 I Tr1, gTr1,
569 I myIter,myThid)
570 ENDIF
571 #endif
572 #ifdef ALLOW_PTRACERS
573 IF ( usePTRACERS ) THEN
574 CALL PTRACERS_INTEGERATE(
575 I bi,bj,k,
576 I xA,yA,uTrans,vTrans,rTrans,maskUp,
577 X KappaRS,
578 I myIter,myTime,myThid)
579 ENDIF
580 #endif /* ALLOW_PTRACERS */
581
582 #ifdef ALLOW_OBCS
583 C-- Apply open boundary conditions
584 IF (useOBCS) THEN
585 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
586 END IF
587 #endif /* ALLOW_OBCS */
588
589 C-- Freeze water
590 IF (allowFreezing) THEN
591 #ifdef ALLOW_AUTODIFF_TAMC
592 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
593 CADJ & , key = kkey, byte = isbyte
594 #endif /* ALLOW_AUTODIFF_TAMC */
595 CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
596 END IF
597
598 C-- end of thermodynamic k loop (Nr:1)
599 ENDDO
600
601
602 #ifdef ALLOW_AUTODIFF_TAMC
603 C? Patrick? What about this one?
604 cph Keys iikey and idkey dont seem to be needed
605 cph since storing occurs on different tape for each
606 cph impldiff call anyways.
607 cph Thus, common block comlev1_impl isnt needed either.
608 cph Storing below needed in the case useGMREDI.
609 iikey = (ikey-1)*maximpl
610 #endif /* ALLOW_AUTODIFF_TAMC */
611
612 C-- Implicit diffusion
613 IF (implicitDiffusion) THEN
614
615 IF (tempStepping) THEN
616 #ifdef ALLOW_AUTODIFF_TAMC
617 idkey = iikey + 1
618 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
619 #endif /* ALLOW_AUTODIFF_TAMC */
620 CALL IMPLDIFF(
621 I bi, bj, iMin, iMax, jMin, jMax,
622 I deltaTtracer, KappaRT, recip_HFacC,
623 U gT,
624 I myThid )
625 ENDIF
626
627 IF (saltStepping) THEN
628 #ifdef ALLOW_AUTODIFF_TAMC
629 idkey = iikey + 2
630 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
631 #endif /* ALLOW_AUTODIFF_TAMC */
632 CALL IMPLDIFF(
633 I bi, bj, iMin, iMax, jMin, jMax,
634 I deltaTtracer, KappaRS, recip_HFacC,
635 U gS,
636 I myThid )
637 ENDIF
638
639 #ifdef ALLOW_PASSIVE_TRACER
640 IF (tr1Stepping) THEN
641 #ifdef ALLOW_AUTODIFF_TAMC
642 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
643 #endif /* ALLOW_AUTODIFF_TAMC */
644 CALL IMPLDIFF(
645 I bi, bj, iMin, iMax, jMin, jMax,
646 I deltaTtracer, KappaRT, recip_HFacC,
647 U gTr1,
648 I myThid )
649 ENDIF
650 #endif
651
652 #ifdef ALLOW_PTRACERS
653 C Vertical diffusion (implicit) for passive tracers
654 IF ( usePTRACERS ) THEN
655 CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
656 ENDIF
657 #endif /* ALLOW_PTRACERS */
658
659 #ifdef ALLOW_OBCS
660 C-- Apply open boundary conditions
661 IF (useOBCS) THEN
662 DO K=1,Nr
663 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
664 ENDDO
665 END IF
666 #endif /* ALLOW_OBCS */
667
668 C-- End If implicitDiffusion
669 ENDIF
670
671 Ccs-
672 ENDDO
673 ENDDO
674
675 #ifdef ALLOW_AIM
676 IF ( useAIM ) THEN
677 CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
678 ENDIF
679 _EXCH_XYZ_R8(gT,myThid)
680 _EXCH_XYZ_R8(gS,myThid)
681 #else
682 IF (staggerTimeStep.AND.useCubedSphereExchange) THEN
683 _EXCH_XYZ_R8(gT,myThid)
684 _EXCH_XYZ_R8(gS,myThid)
685 ENDIF
686 #endif /* ALLOW_AIM */
687
688 #ifndef DISABLE_DEBUGMODE
689 If (debugMode) THEN
690 CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
691 CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
692 CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
693 CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
694 CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
695 CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
696 CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
697 CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
698 CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
699 #ifdef ALLOW_PTRACERS
700 IF ( usePTRACERS ) THEN
701 CALL PTRACERS_DEBUG(myThid)
702 ENDIF
703 #endif /* ALLOW_PTRACERS */
704 ENDIF
705 #endif
706
707 RETURN
708 END

  ViewVC Help
Powered by ViewVC 1.1.22