/[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.50 - (show annotations) (download)
Wed Jun 21 19:13:11 2000 UTC (23 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint29
Changes since 1.49: +41 -68 lines
Packaged KPP and GM/Redi:
 o moved temporary arrays such as K13, KapGM,... to private common blocks
 o calls to KPP_CALC and GMREDI_CALC_TENSOR replace KVMIX and CALC_ISOSLOPES.

Added call to gradSigma() to calculate density gradients.

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

  ViewVC Help
Powered by ViewVC 1.1.22