/[MITgcm]/MITgcm/model/src/dynamics.F
ViewVC logotype

Contents of /MITgcm/model/src/dynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.53 - (show annotations) (download)
Mon Sep 11 23:07:29 2000 UTC (23 years, 8 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint31
Changes since 1.52: +257 -206 lines
Various corrections and additions of store directives for TAMC.
Changes of interfaces to packaged GMRedi and KPP.
Tested for exp(0,2,4).

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.52 2000/06/29 18:49:50 adcroft Exp $
2
3 #include "CPP_OPTIONS.h"
4
5 SUBROUTINE DYNAMICS(myTime, myIter, myThid)
6 C /==========================================================\
7 C | SUBROUTINE DYNAMICS |
8 C | o Controlling routine for the explicit part of the model |
9 C | dynamics. |
10 C |==========================================================|
11 C | This routine evaluates the "dynamics" terms for each |
12 C | block of ocean in turn. Because the blocks of ocean have |
13 C | overlap regions they are independent of one another. |
14 C | If terms involving lateral integrals are needed in this |
15 C | routine care will be needed. Similarly finite-difference |
16 C | operations with stencils wider than the overlap region |
17 C | require special consideration. |
18 C | Notes |
19 C | ===== |
20 C | C*P* comments indicating place holders for which code is |
21 C | presently being developed. |
22 C \==========================================================/
23 IMPLICIT NONE
24
25 C == Global variables ===
26 #include "SIZE.h"
27 #include "EEPARAMS.h"
28 #include "CG2D.h"
29 #include "PARAMS.h"
30 #include "DYNVARS.h"
31 #include "GRID.h"
32
33 #ifdef ALLOW_AUTODIFF_TAMC
34 # include "tamc.h"
35 # include "tamc_keys.h"
36 # ifdef ALLOW_KPP
37 # include "KPP.h"
38 # endif
39 #endif /* ALLOW_AUTODIFF_TAMC */
40
41 C == Routine arguments ==
42 C myTime - Current time in simulation
43 C myIter - Current iteration number in simulation
44 C myThid - Thread number for this instance of the routine.
45 _RL myTime
46 INTEGER myIter
47 INTEGER myThid
48
49 C == Local variables
50 C xA, yA - Per block temporaries holding face areas
51 C uTrans, vTrans, rTrans - Per block temporaries holding flow
52 C transport
53 C rVel o uTrans: Zonal transport
54 C o vTrans: Meridional transport
55 C o rTrans: Vertical transport
56 C o rVel: Vertical velocity at upper and
57 C lower cell faces.
58 C maskC,maskUp o maskC: land/water mask for tracer cells
59 C o maskUp: land/water mask for W points
60 C aTerm, xTerm, cTerm - Work arrays for holding separate terms in
61 C mTerm, pTerm, tendency equations.
62 C fZon, fMer, fVer[STUV] o aTerm: Advection term
63 C o xTerm: Mixing term
64 C o cTerm: Coriolis term
65 C o mTerm: Metric term
66 C o pTerm: Pressure term
67 C o fZon: Zonal flux term
68 C o fMer: Meridional flux term
69 C o fVer: Vertical flux term - note fVer
70 C is "pipelined" in the vertical
71 C so we need an fVer for each
72 C variable.
73 C rhoK, rhoKM1 - Density at current level, level above and level
74 C below.
75 C rhoKP1
76 C buoyK, buoyKM1 - Buoyancy at current level and level above.
77 C phiHyd - Hydrostatic part of the potential phiHydi.
78 C In z coords phiHydiHyd is the hydrostatic
79 C pressure anomaly
80 C In p coords phiHydiHyd is the geopotential
81 C surface height
82 C anomaly.
83 C etaSurfX, - Holds surface elevation gradient in X and Y.
84 C etaSurfY
85 C KappaRT, - Total diffusion in vertical for T and S.
86 C KappaRS (background + spatially varying, isopycnal term).
87 C iMin, iMax - Ranges and sub-block indices on which calculations
88 C jMin, jMax are applied.
89 C bi, bj
90 C k, kup, - Index for layer above and below. kup and kDown
91 C kDown, km1 are switched with layer to be the appropriate
92 C index into fVerTerm.
93 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL rVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
99 _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101 _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103 _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105 _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
109 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
110 _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
111 _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
112 _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
113 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114 _RL rhokp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116 _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117 _RL buoyK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118 _RL rhotmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119 _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120 _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121 _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
122 _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
123 _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
124 _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
125 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
126 _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
127 _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
128
129 C This is currently also used by IVDC and Diagnostics
130 C #ifdef INCLUDE_CONVECT_CALL
131 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
132 C #endif
133
134 INTEGER iMin, iMax
135 INTEGER jMin, jMax
136 INTEGER bi, bj
137 INTEGER i, j
138 INTEGER k, km1, kup, kDown
139 LOGICAL BOTTOM_LAYER
140
141 #ifdef ALLOW_AUTODIFF_TAMC
142 INTEGER isbyte
143 PARAMETER( isbyte = 4 )
144
145 INTEGER act1, act2, act3, act4
146 INTEGER max1, max2, max3
147 INTEGER iikey, kkey
148 INTEGER maximpl
149 #endif /* ALLOW_AUTODIFF_TAMC */
150
151 C--- The algorithm...
152 C
153 C "Correction Step"
154 C =================
155 C Here we update the horizontal velocities with the surface
156 C pressure such that the resulting flow is either consistent
157 C with the free-surface evolution or the rigid-lid:
158 C U[n] = U* + dt x d/dx P
159 C V[n] = V* + dt x d/dy P
160 C
161 C "Calculation of Gs"
162 C ===================
163 C This is where all the accelerations and tendencies (ie.
164 C physics, parameterizations etc...) are calculated
165 C rVel = sum_r ( div. u[n] )
166 C rho = rho ( theta[n], salt[n] )
167 C b = b(rho, theta)
168 C K31 = K31 ( rho )
169 C Gu[n] = Gu( u[n], v[n], rVel, b, ... )
170 C Gv[n] = Gv( u[n], v[n], rVel, b, ... )
171 C Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )
172 C Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )
173 C
174 C "Time-stepping" or "Prediction"
175 C ================================
176 C The models variables are stepped forward with the appropriate
177 C time-stepping scheme (currently we use Adams-Bashforth II)
178 C - For momentum, the result is always *only* a "prediction"
179 C in that the flow may be divergent and will be "corrected"
180 C later with a surface pressure gradient.
181 C - Normally for tracers the result is the new field at time
182 C level [n+1} *BUT* in the case of implicit diffusion the result
183 C is also *only* a prediction.
184 C - We denote "predictors" with an asterisk (*).
185 C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
186 C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
187 C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
188 C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
189 C With implicit diffusion:
190 C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
191 C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
192 C (1 + dt * K * d_zz) theta[n] = theta*
193 C (1 + dt * K * d_zz) salt[n] = salt*
194 C---
195
196 #ifdef ALLOW_AUTODIFF_TAMC
197 C-- dummy statement to end declaration part
198 ikey = 1
199 #endif /* ALLOW_AUTODIFF_TAMC */
200
201 C-- Set up work arrays with valid (i.e. not NaN) values
202 C These inital values do not alter the numerical results. They
203 C just ensure that all memory references are to valid floating
204 C point numbers. This prevents spurious hardware signals due to
205 C uninitialised but inert locations.
206 DO j=1-OLy,sNy+OLy
207 DO i=1-OLx,sNx+OLx
208 xA(i,j) = 0. _d 0
209 yA(i,j) = 0. _d 0
210 uTrans(i,j) = 0. _d 0
211 vTrans(i,j) = 0. _d 0
212 aTerm(i,j) = 0. _d 0
213 xTerm(i,j) = 0. _d 0
214 cTerm(i,j) = 0. _d 0
215 mTerm(i,j) = 0. _d 0
216 pTerm(i,j) = 0. _d 0
217 fZon(i,j) = 0. _d 0
218 fMer(i,j) = 0. _d 0
219 DO k=1,Nr
220 phiHyd (i,j,k) = 0. _d 0
221 KappaRU(i,j,k) = 0. _d 0
222 KappaRV(i,j,k) = 0. _d 0
223 sigmaX(i,j,k) = 0. _d 0
224 sigmaY(i,j,k) = 0. _d 0
225 sigmaR(i,j,k) = 0. _d 0
226 ENDDO
227 rhoKM1 (i,j) = 0. _d 0
228 rhok (i,j) = 0. _d 0
229 rhoKP1 (i,j) = 0. _d 0
230 rhoTMP (i,j) = 0. _d 0
231 buoyKM1(i,j) = 0. _d 0
232 buoyK (i,j) = 0. _d 0
233 maskC (i,j) = 0. _d 0
234 ENDDO
235 ENDDO
236
237
238 #ifdef ALLOW_AUTODIFF_TAMC
239 C-- HPF directive to help TAMC
240 CHPF$ INDEPENDENT
241 #endif /* ALLOW_AUTODIFF_TAMC */
242
243 DO bj=myByLo(myThid),myByHi(myThid)
244
245 #ifdef ALLOW_AUTODIFF_TAMC
246 C-- HPF directive to help TAMC
247 CHPF$ INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV
248 CHPF$& ,phiHyd,utrans,vtrans,maskc,xA,yA
249 CHPF$& ,KappaRT,KappaRS,KappaRU,KappaRV
250 CHPF$& )
251 #endif /* ALLOW_AUTODIFF_TAMC */
252
253 DO bi=myBxLo(myThid),myBxHi(myThid)
254
255 #ifdef ALLOW_AUTODIFF_TAMC
256 act1 = bi - myBxLo(myThid)
257 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
258
259 act2 = bj - myByLo(myThid)
260 max2 = myByHi(myThid) - myByLo(myThid) + 1
261
262 act3 = myThid - 1
263 max3 = nTx*nTy
264
265 act4 = ikey_dynamics - 1
266
267 ikey = (act1 + 1) + act2*max1
268 & + act3*max1*max2
269 & + act4*max1*max2*max3
270 #endif /* ALLOW_AUTODIFF_TAMC */
271
272 C-- Set up work arrays that need valid initial values
273 DO j=1-OLy,sNy+OLy
274 DO i=1-OLx,sNx+OLx
275 rTrans(i,j) = 0. _d 0
276 rVel (i,j,1) = 0. _d 0
277 rVel (i,j,2) = 0. _d 0
278 fVerT (i,j,1) = 0. _d 0
279 fVerT (i,j,2) = 0. _d 0
280 fVerS (i,j,1) = 0. _d 0
281 fVerS (i,j,2) = 0. _d 0
282 fVerU (i,j,1) = 0. _d 0
283 fVerU (i,j,2) = 0. _d 0
284 fVerV (i,j,1) = 0. _d 0
285 fVerV (i,j,2) = 0. _d 0
286 phiHyd(i,j,1) = 0. _d 0
287 ENDDO
288 ENDDO
289
290 DO k=1,Nr
291 DO j=1-OLy,sNy+OLy
292 DO i=1-OLx,sNx+OLx
293 #ifdef INCLUDE_CONVECT_CALL
294 ConvectCount(i,j,k) = 0.
295 #endif
296 KappaRT(i,j,k) = 0. _d 0
297 KappaRS(i,j,k) = 0. _d 0
298 ENDDO
299 ENDDO
300 ENDDO
301
302 iMin = 1-OLx+1
303 iMax = sNx+OLx
304 jMin = 1-OLy+1
305 jMax = sNy+OLy
306
307
308 k = 1
309 BOTTOM_LAYER = k .EQ. Nr
310
311 #ifdef DO_PIPELINED_CORRECTION_STEP
312 C-- Calculate gradient of surface pressure
313 CALL CALC_GRAD_ETA_SURF(
314 I bi,bj,iMin,iMax,jMin,jMax,
315 O etaSurfX,etaSurfY,
316 I myThid)
317 C-- Update fields in top level according to tendency terms
318 CALL CORRECTION_STEP(
319 I bi,bj,iMin,iMax,jMin,jMax,k,
320 I etaSurfX,etaSurfY,myTime,myThid)
321
322 #ifdef ALLOW_OBCS
323 IF (openBoundaries) THEN
324 #ifdef ALLOW_AUTODIFF_TAMC
325 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
326 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
327 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
328 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
329 #endif /* ALLOW_AUTODIFF_TAMC */
330 CALL APPLY_OBCS1( bi, bj, k, myThid )
331 END IF
332 #endif
333
334 IF ( .NOT. BOTTOM_LAYER ) THEN
335 C-- Update fields in layer below according to tendency terms
336 CALL CORRECTION_STEP(
337 I bi,bj,iMin,iMax,jMin,jMax,k+1,
338 I etaSurfX,etaSurfY,myTime,myThid)
339 #ifdef ALLOW_OBCS
340 IF (openBoundaries) THEN
341 #ifdef ALLOW_AUTODIFF_TAMC
342 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
343 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
344 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
345 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
346 #endif /* ALLOW_AUTODIFF_TAMC */
347 CALL APPLY_OBCS1( bi, bj, k+1, myThid )
348 END IF
349 #endif
350 ENDIF
351 #endif
352 C-- Density of 1st level (below W(1)) reference to level 1
353 #ifdef INCLUDE_FIND_RHO_CALL
354 #ifdef ALLOW_AUTODIFF_TAMC
355 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
356 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
357 #endif /* ALLOW_AUTODIFF_TAMC */
358 CALL FIND_RHO(
359 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
360 O rhoKm1,
361 I myThid )
362 #endif
363
364 IF (.NOT. BOTTOM_LAYER) THEN
365
366 C-- Check static stability with layer below
367 C-- and mix as needed.
368 #ifdef INCLUDE_FIND_RHO_CALL
369 #ifdef ALLOW_AUTODIFF_TAMC
370 CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj
371 CADJ & , key = ikey, byte = isbyte
372 CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj
373 CADJ & , key = ikey, byte = isbyte
374 #endif /* ALLOW_AUTODIFF_TAMC */
375 CALL FIND_RHO(
376 I bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,
377 O rhoKp1,
378 I myThid )
379 #endif
380
381 #ifdef INCLUDE_CONVECT_CALL
382
383 #ifdef ALLOW_AUTODIFF_TAMC
384 CADJ STORE rhoKm1(:,:) = comlev1_bibj, key = ikey, byte = isbyte
385 CADJ STORE rhoKp1(:,:) = comlev1_bibj, key = ikey, byte = isbyte
386 #endif /* ALLOW_AUTODIFF_TAMC */
387 CALL CONVECT(
388 I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
389 U ConvectCount,
390 I myTime,myIter,myThid)
391
392 #ifdef ALLOW_AUTODIFF_TAMC
393 CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)
394 CADJ & = comlev1_bibj, key = ikey, byte = isbyte
395 CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)
396 CADJ & = comlev1_bibj, key = ikey, byte = isbyte
397 #endif /* ALLOW_AUTODIFF_TAMC */
398
399 #endif
400
401 C-- Implicit Vertical Diffusion for Convection
402 IF (ivdc_kappa.NE.0.) THEN
403 CALL CALC_IVDC(
404 I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
405 U ConvectCount, KappaRT, KappaRS,
406 I myTime,myIter,myThid)
407 ENDIF
408
409 C-- Recompute density after mixing
410 #ifdef INCLUDE_FIND_RHO_CALL
411 CALL FIND_RHO(
412 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
413 O rhoKm1,
414 I myThid )
415 #endif
416 ENDIF
417
418 C-- Calculate buoyancy
419 CALL CALC_BUOYANCY(
420 I bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,
421 O buoyKm1,
422 I myThid )
423
424 C-- Integrate hydrostatic balance for phiHyd with BC of
425 C-- phiHyd(z=0)=0
426 CALL CALC_PHI_HYD(
427 I bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyKm1,
428 U phiHyd,
429 I myThid )
430
431 #ifdef ALLOW_GMREDI
432 CALL GRAD_SIGMA(
433 I bi, bj, iMin, iMax, jMin, jMax, k,
434 I rhoKm1, rhoKm1, rhoKm1,
435 O sigmaX, sigmaY, sigmaR,
436 I myThid )
437 #endif
438
439 C-- Start of downward loop
440 DO k=2,Nr
441
442 #ifdef ALLOW_AUTODIFF_TAMC
443 kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
444 #endif /* ALLOW_AUTODIFF_TAMC */
445
446 BOTTOM_LAYER = k .EQ. Nr
447
448 #ifdef DO_PIPELINED_CORRECTION_STEP
449 IF ( .NOT. BOTTOM_LAYER ) THEN
450 C-- Update fields in layer below according to tendency terms
451 CALL CORRECTION_STEP(
452 I bi,bj,iMin,iMax,jMin,jMax,k+1,
453 I etaSurfX,etaSurfY,myTime,myThid)
454 #ifdef ALLOW_OBCS
455 IF (openBoundaries) THEN
456 #ifdef ALLOW_AUTODIFF_TAMC
457 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k
458 CADJ & , key = kkey, byte = isbyte
459 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k
460 CADJ & , key = kkey, byte = isbyte
461 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
462 CADJ & , key = kkey, byte = isbyte
463 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
464 CADJ & , key = kkey, byte = isbyte
465 #endif /* ALLOW_AUTODIFF_TAMC */
466 CALL APPLY_OBCS1( bi, bj, k+1, myThid )
467 END IF
468 #endif
469 ENDIF
470 #endif /* DO_PIPELINED_CORRECTION_STEP */
471
472 C-- Density of k level (below W(k)) reference to k level
473 #ifdef INCLUDE_FIND_RHO_CALL
474 #ifdef ALLOW_AUTODIFF_TAMC
475 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
476 CADJ & , key = kkey, byte = isbyte
477 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
478 CADJ & , key = kkey, byte = isbyte
479 #endif /* ALLOW_AUTODIFF_TAMC */
480 CALL FIND_RHO(
481 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
482 O rhoK,
483 I myThid )
484
485 #ifdef ALLOW_AUTODIFF_TAMC
486 CADJ STORE rhoK(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
487 #endif /* ALLOW_AUTODIFF_TAMC */
488 #endif
489
490 IF (.NOT. BOTTOM_LAYER) THEN
491
492 C-- Check static stability with layer below and mix as needed.
493 C-- Density of k+1 level (below W(k+1)) reference to k level.
494 #ifdef INCLUDE_FIND_RHO_CALL
495 #ifdef ALLOW_AUTODIFF_TAMC
496 CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj_k
497 CADJ & , key = kkey, byte = isbyte
498 CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj_k
499 CADJ & , key = kkey, byte = isbyte
500 #endif /* ALLOW_AUTODIFF_TAMC */
501 CALL FIND_RHO(
502 I bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,
503 O rhoKp1,
504 I myThid )
505 #ifdef ALLOW_AUTODIFF_TAMC
506 CADJ STORE rhoKp1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
507 #endif /* ALLOW_AUTODIFF_TAMC */
508 #endif
509
510 #ifdef INCLUDE_CONVECT_CALL
511 CALL CONVECT(
512 I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,
513 U ConvectCount,
514 I myTime,myIter,myThid)
515
516 #endif
517
518 C-- Implicit Vertical Diffusion for Convection
519 IF (ivdc_kappa.NE.0.) THEN
520 #ifdef ALLOW_AUTODIFF_TAMC
521 CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
522 #endif /* ALLOW_AUTODIFF_TAMC */
523 CALL CALC_IVDC(
524 I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
525 U ConvectCount, KappaRT, KappaRS,
526 I myTime,myIter,myThid)
527 END IF
528
529 C-- Recompute density after mixing
530 #ifdef INCLUDE_FIND_RHO_CALL
531 #ifdef ALLOW_AUTODIFF_TAMC
532 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
533 CADJ & , key = kkey, byte = isbyte
534 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
535 CADJ & , key = kkey, byte = isbyte
536 #endif /* ALLOW_AUTODIFF_TAMC */
537 CALL FIND_RHO(
538 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
539 O rhoK,
540 I myThid )
541 #endif
542
543 C-- IF (.NOT. BOTTOM_LAYER) ends here
544 ENDIF
545
546 C-- Calculate buoyancy
547 CALL CALC_BUOYANCY(
548 I bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
549 O buoyK,
550 I myThid )
551
552 C-- Integrate hydrostatic balance for phiHyd with BC of
553 C-- phiHyd(z=0)=0
554 CALL CALC_PHI_HYD(
555 I bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
556 U phiHyd,
557 I myThid )
558
559 #ifdef INCLUDE_FIND_RHO_CALL
560 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
561
562 #ifdef ALLOW_AUTODIFF_TAMC
563 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k
564 CADJ & , key = kkey, byte = isbyte
565 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k
566 CADJ & , key = kkey, byte = isbyte
567 #endif /* ALLOW_AUTODIFF_TAMC */
568
569 CALL FIND_RHO(
570 I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
571 O rhoTmp,
572 I myThid )
573 #endif
574
575
576 #ifdef ALLOW_GMREDI
577 CALL GRAD_SIGMA(
578 I bi, bj, iMin, iMax, jMin, jMax, k,
579 I rhoK, rhotmp, rhoK,
580 O sigmaX, sigmaY, sigmaR,
581 I myThid )
582 #endif
583
584 DO J=jMin,jMax
585 DO I=iMin,iMax
586 #ifdef INCLUDE_FIND_RHO_CALL
587 rhoKm1 (I,J) = rhoK(I,J)
588 #endif
589 buoyKm1(I,J) = buoyK(I,J)
590 ENDDO
591 ENDDO
592
593 C-- end of k loop
594 ENDDO
595
596 #ifdef ALLOW_GMREDI
597 IF (useGMRedi) THEN
598 DO k=1, Nr
599 CALL GMREDI_CALC_TENSOR(
600 I bi, bj, iMin, iMax, jMin, jMax, k,
601 I sigmaX, sigmaY, sigmaR,
602 I myThid )
603 ENDDO
604 ENDIF
605 #endif
606
607 #ifdef ALLOW_AUTODIFF_TAMC
608 CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key = ikey, byte = isbyte
609 CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key = ikey, byte = isbyte
610
611 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
612 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
613 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
614 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
615
616 C-- dummy initialization to break data flow because
617 C-- calc_div_ghat has a condition for initialization
618 DO J=jMin,jMax
619 DO I=iMin,iMax
620 cg2d_b(i,j,bi,bj) = 0.0
621 ENDDO
622 ENDDO
623 #endif /* ALLOW_AUTODIFF_TAMC */
624
625
626 #ifdef ALLOW_KPP
627 C-- Compute KPP mixing coefficients
628 IF (useKPP) THEN
629
630 CALL TIMER_START('KPP_CALC [DYNAMICS]', myThid)
631 CALL KPP_CALC(
632 I bi, bj, myTime, myThid )
633 CALL TIMER_STOP ('KPP_CALC [DYNAMICS]', myThid)
634
635 #ifdef ALLOW_AUTODIFF_TAMC
636 ELSE
637 DO j=1-OLy,sNy+OLy
638 DO i=1-OLx,sNx+OLx
639 KPPhbl (i,j,bi,bj) = 1.0
640 KPPfrac(i,j,bi,bj) = 0.0
641 DO k = 1,Nr
642 KPPghat (i,j,k,bi,bj) = 0.0
643 KPPviscAz (i,j,k,bi,bj) = viscAz
644 KPPdiffKzT(i,j,k,bi,bj) = diffKzT
645 KPPdiffKzS(i,j,k,bi,bj) = diffKzS
646 ENDDO
647 ENDDO
648 ENDDO
649 #endif /* ALLOW_AUTODIFF_TAMC */
650 ENDIF
651
652 #ifdef ALLOW_AUTODIFF_TAMC
653 CADJ STORE KPPghat (:,:,:,bi,bj)
654 CADJ & , KPPviscAz (:,:,:,bi,bj)
655 CADJ & , KPPdiffKzT(:,:,:,bi,bj)
656 CADJ & , KPPdiffKzS(:,:,:,bi,bj)
657 CADJ & , KPPfrac (:,: ,bi,bj)
658 CADJ & = comlev1_bibj, key = ikey, byte = isbyte
659 #endif /* ALLOW_AUTODIFF_TAMC */
660
661 #endif /* ALLOW_KPP */
662
663 C-- Start of upward loop
664 DO k = Nr, 1, -1
665
666 C-- km1 Points to level above k (=k-1)
667 C-- kup Cycles through 1,2 to point to layer above
668 C-- kDown Cycles through 2,1 to point to current layer
669
670 km1 =max(1,k-1)
671 kup =1+MOD(k+1,2)
672 kDown=1+MOD(k,2)
673
674 iMin = 1-OLx+2
675 iMax = sNx+OLx-1
676 jMin = 1-OLy+2
677 jMax = sNy+OLy-1
678
679 #ifdef ALLOW_AUTODIFF_TAMC
680 kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
681
682 CADJ STORE rvel (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte
683 CADJ STORE rTrans(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
684 CADJ STORE KappaRT(:,:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
685 CADJ STORE KappaRS(:,:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
686 #endif /* ALLOW_AUTODIFF_TAMC */
687
688 C-- Get temporary terms used by tendency routines
689 CALL CALC_COMMON_FACTORS (
690 I bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
691 O xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
692 I myThid)
693
694 #ifdef ALLOW_OBCS
695 IF (openBoundaries) THEN
696 CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )
697 ENDIF
698 #endif
699
700 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
701 C-- Calculate the total vertical diffusivity
702 CALL CALC_DIFFUSIVITY(
703 I bi,bj,iMin,iMax,jMin,jMax,k,
704 I maskC,maskUp,
705 O KappaRT,KappaRS,KappaRU,KappaRV,
706 I myThid)
707 #endif
708 C-- Calculate accelerations in the momentum equations
709 IF ( momStepping ) THEN
710 CALL CALC_MOM_RHS(
711 I bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
712 I xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
713 I phiHyd,KappaRU,KappaRV,
714 U aTerm,xTerm,cTerm,mTerm,pTerm,
715 U fZon, fMer, fVerU, fVerV,
716 I myTime, myThid)
717 #ifdef ALLOW_AUTODIFF_TAMC
718 #ifdef INCLUDE_CD_CODE
719 ELSE
720 DO j=1-OLy,sNy+OLy
721 DO i=1-OLx,sNx+OLx
722 guCD(i,j,k,bi,bj) = 0.0
723 gvCD(i,j,k,bi,bj) = 0.0
724 END DO
725 END DO
726 #endif
727 #endif /* ALLOW_AUTODIFF_TAMC */
728 ENDIF
729 C-- Calculate active tracer tendencies
730 IF ( tempStepping ) THEN
731 CALL CALC_GT(
732 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
733 I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
734 I KappaRT,
735 U aTerm,xTerm,fZon,fMer,fVerT,
736 I myTime, myThid)
737 ENDIF
738 IF ( saltStepping ) THEN
739 CALL CALC_GS(
740 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
741 I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
742 I KappaRS,
743 U aTerm,xTerm,fZon,fMer,fVerS,
744 I myTime, myThid)
745 ENDIF
746 #ifdef ALLOW_OBCS
747 C-- Calculate future values on open boundaries
748 IF (openBoundaries) THEN
749 Caja CALL CYCLE_OBCS( k, bi, bj, myThid )
750 CALL SET_OBCS( k, bi, bj, myTime+deltaTclock, myThid )
751 ENDIF
752 #endif
753 C-- Prediction step (step forward all model variables)
754 CALL TIMESTEP(
755 I bi,bj,iMin,iMax,jMin,jMax,k,
756 I myIter, myThid)
757 #ifdef ALLOW_OBCS
758 C-- Apply open boundary conditions
759 IF (openBoundaries) THEN
760 #ifdef ALLOW_AUTODIFF_TAMC
761 CADJ STORE gunm1(:,:,k,bi,bj) = comlev1_bibj_k
762 CADJ & , key = kkey, byte = isbyte
763 CADJ STORE gvnm1(:,:,k,bi,bj) = comlev1_bibj_k
764 CADJ & , key = kkey, byte = isbyte
765 CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k, key = kkey, byte = isbyte
766 #endif /* ALLOW_AUTODIFF_TAMC */
767 CADJ &
768 CALL APPLY_OBCS2( bi, bj, k, myThid )
769 END IF
770 #endif
771 C-- Freeze water
772 IF (allowFreezing) THEN
773 #ifdef ALLOW_AUTODIFF_TAMC
774 CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
775 CADJ & , key = kkey, byte = isbyte
776 #endif /* ALLOW_AUTODIFF_TAMC */
777 CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
778 END IF
779
780 #ifdef DIVG_IN_DYNAMICS
781 C-- Diagnose barotropic divergence of predicted fields
782 CALL CALC_DIV_GHAT(
783 I bi,bj,iMin,iMax,jMin,jMax,k,
784 I xA,yA,
785 I myThid)
786 #endif /* DIVG_IN_DYNAMICS */
787
788 C-- Cumulative diagnostic calculations (ie. time-averaging)
789 #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
790 IF (taveFreq.GT.0.) THEN
791 CALL DO_TIME_AVERAGES(
792 I myTime, myIter, bi, bj, k, kup, kDown,
793 I rVel, ConvectCount,
794 I myThid )
795 ENDIF
796 #endif
797
798
799 C-- k loop
800 ENDDO
801
802 #ifdef ALLOW_AUTODIFF_TAMC
803 maximpl = 6
804 iikey = (ikey-1)*maximpl
805 #endif /* ALLOW_AUTODIFF_TAMC */
806
807 C-- Implicit diffusion
808 IF (implicitDiffusion) THEN
809
810 IF (tempStepping) THEN
811 #ifdef ALLOW_AUTODIFF_TAMC
812 idkey = iikey + 1
813 #endif /* ALLOW_AUTODIFF_TAMC */
814 CALL IMPLDIFF(
815 I bi, bj, iMin, iMax, jMin, jMax,
816 I deltaTtracer, KappaRT,recip_HFacC,
817 U gTNm1,
818 I myThid )
819 END IF
820
821 IF (saltStepping) THEN
822 #ifdef ALLOW_AUTODIFF_TAMC
823 idkey = iikey + 2
824 #endif /* ALLOW_AUTODIFF_TAMC */
825 CALL IMPLDIFF(
826 I bi, bj, iMin, iMax, jMin, jMax,
827 I deltaTtracer, KappaRS,recip_HFacC,
828 U gSNm1,
829 I myThid )
830 END IF
831
832 C-- implicitDiffusion
833 ENDIF
834
835 C-- Implicit viscosity
836 IF (implicitViscosity) THEN
837
838 IF (momStepping) THEN
839 #ifdef ALLOW_AUTODIFF_TAMC
840 idkey = iikey + 3
841 #endif /* ALLOW_AUTODIFF_TAMC */
842 CALL IMPLDIFF(
843 I bi, bj, iMin, iMax, jMin, jMax,
844 I deltaTmom, KappaRU,recip_HFacW,
845 U gUNm1,
846 I myThid )
847 #ifdef ALLOW_AUTODIFF_TAMC
848 idkey = iikey + 4
849 #endif /* ALLOW_AUTODIFF_TAMC */
850 CALL IMPLDIFF(
851 I bi, bj, iMin, iMax, jMin, jMax,
852 I deltaTmom, KappaRV,recip_HFacS,
853 U gVNm1,
854 I myThid )
855
856 #ifdef INCLUDE_CD_CODE
857
858 #ifdef ALLOW_AUTODIFF_TAMC
859 idkey = iikey + 5
860 #endif /* ALLOW_AUTODIFF_TAMC */
861 CALL IMPLDIFF(
862 I bi, bj, iMin, iMax, jMin, jMax,
863 I deltaTmom, KappaRU,recip_HFacW,
864 U vVelD,
865 I myThid )
866 #ifdef ALLOW_AUTODIFF_TAMC
867 idkey = iikey + 6
868 #endif /* ALLOW_AUTODIFF_TAMC */
869 CALL IMPLDIFF(
870 I bi, bj, iMin, iMax, jMin, jMax,
871 I deltaTmom, KappaRV,recip_HFacS,
872 U uVelD,
873 I myThid )
874
875 #endif
876
877 C-- momStepping
878 ENDIF
879
880 C-- implicitViscosity
881 ENDIF
882
883 ENDDO
884 ENDDO
885
886
887 RETURN
888 END

  ViewVC Help
Powered by ViewVC 1.1.22