/[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.52 - (show annotations) (download)
Thu Jun 29 18:49:50 2000 UTC (23 years, 10 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint30
Changes since 1.51: +4 -3 lines
The array ConvectCount(...) needs to always be declared becuase
it is used by IVDC and the time-averaging package. We should ultimately
move this into a common block specific to those routines that use it.

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

  ViewVC Help
Powered by ViewVC 1.1.22