/[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.74 - (show annotations) (download)
Mon Jul 30 20:37:45 2001 UTC (22 years, 9 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint40pre4, checkpoint40pre5
Changes since 1.73: +26 -11 lines
Extended iMin,jMin range for calc_common_factors, calc_diffusivity.

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.73 2001/07/20 19:16:28 adcroft Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 SUBROUTINE DYNAMICS(myTime, myIter, myThid)
7 C /==========================================================\
8 C | SUBROUTINE DYNAMICS |
9 C | o Controlling routine for the explicit part of the model |
10 C | dynamics. |
11 C |==========================================================|
12 C | This routine evaluates the "dynamics" terms for each |
13 C | block of ocean in turn. Because the blocks of ocean have |
14 C | overlap regions they are independent of one another. |
15 C | If terms involving lateral integrals are needed in this |
16 C | routine care will be needed. Similarly finite-difference |
17 C | operations with stencils wider than the overlap region |
18 C | require special consideration. |
19 C | Notes |
20 C | ===== |
21 C | C*P* comments indicating place holders for which code is |
22 C | presently being developed. |
23 C \==========================================================/
24 IMPLICIT NONE
25
26 C == Global variables ===
27 #include "SIZE.h"
28 #include "EEPARAMS.h"
29 #include "PARAMS.h"
30 #include "DYNVARS.h"
31 #include "GRID.h"
32 #ifdef ALLOW_PASSIVE_TRACER
33 #include "TR1.h"
34 #endif
35
36 #ifdef ALLOW_AUTODIFF_TAMC
37 # include "tamc.h"
38 # include "tamc_keys.h"
39 # include "FFIELDS.h"
40 # ifdef ALLOW_KPP
41 # include "KPP.h"
42 # endif
43 # ifdef ALLOW_GMREDI
44 # include "GMREDI.h"
45 # endif
46 #endif /* ALLOW_AUTODIFF_TAMC */
47
48 #ifdef ALLOW_TIMEAVE
49 #include "TIMEAVE_STATV.h"
50 #endif
51
52 C == Routine arguments ==
53 C myTime - Current time in simulation
54 C myIter - Current iteration number in simulation
55 C myThid - Thread number for this instance of the routine.
56 _RL myTime
57 INTEGER myIter
58 INTEGER myThid
59
60 C == Local variables
61 C xA, yA - Per block temporaries holding face areas
62 C uTrans, vTrans, rTrans - Per block temporaries holding flow
63 C transport
64 C o uTrans: Zonal transport
65 C o vTrans: Meridional transport
66 C o rTrans: Vertical transport
67 C maskUp o maskUp: land/water mask for W points
68 C fVer[STUV] o fVer: Vertical flux term - note fVer
69 C is "pipelined" in the vertical
70 C so we need an fVer for each
71 C variable.
72 C rhoK, rhoKM1 - Density at current level, and level above
73 C phiHyd - Hydrostatic part of the potential phiHydi.
74 C In z coords phiHydiHyd is the hydrostatic
75 C Potential (=pressure/rho0) anomaly
76 C In p coords phiHydiHyd is the geopotential
77 C surface height anomaly.
78 C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
79 C phiSurfY or geopotentiel (atmos) in X and Y direction
80 C KappaRT, - Total diffusion in vertical for T and S.
81 C KappaRS (background + spatially varying, isopycnal term).
82 C iMin, iMax - Ranges and sub-block indices on which calculations
83 C jMin, jMax are applied.
84 C bi, bj
85 C k, kup, - Index for layer above and below. kup and kDown
86 C kDown, km1 are switched with layer to be the appropriate
87 C index into fVerTerm.
88 C tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.
89 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
96 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
97 _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
98 _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
99 _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
100 _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105 _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
106 _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
107 _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
108 _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
109 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110 _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
111 _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112 _RL tauAB
113
114 C This is currently used by IVDC and Diagnostics
115 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116
117 INTEGER iMin, iMax
118 INTEGER jMin, jMax
119 INTEGER bi, bj
120 INTEGER i, j
121 INTEGER k, km1, kup, kDown
122
123 Cjmc : add for phiHyd output <- but not working if multi tile per CPU
124 c CHARACTER*(MAX_LEN_MBUF) suff
125 c LOGICAL DIFFERENT_MULTIPLE
126 c EXTERNAL DIFFERENT_MULTIPLE
127 Cjmc(end)
128
129 C--- The algorithm...
130 C
131 C "Correction Step"
132 C =================
133 C Here we update the horizontal velocities with the surface
134 C pressure such that the resulting flow is either consistent
135 C with the free-surface evolution or the rigid-lid:
136 C U[n] = U* + dt x d/dx P
137 C V[n] = V* + dt x d/dy P
138 C
139 C "Calculation of Gs"
140 C ===================
141 C This is where all the accelerations and tendencies (ie.
142 C physics, parameterizations etc...) are calculated
143 C rho = rho ( theta[n], salt[n] )
144 C b = b(rho, theta)
145 C K31 = K31 ( rho )
146 C Gu[n] = Gu( u[n], v[n], wVel, b, ... )
147 C Gv[n] = Gv( u[n], v[n], wVel, b, ... )
148 C Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
149 C Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
150 C
151 C "Time-stepping" or "Prediction"
152 C ================================
153 C The models variables are stepped forward with the appropriate
154 C time-stepping scheme (currently we use Adams-Bashforth II)
155 C - For momentum, the result is always *only* a "prediction"
156 C in that the flow may be divergent and will be "corrected"
157 C later with a surface pressure gradient.
158 C - Normally for tracers the result is the new field at time
159 C level [n+1} *BUT* in the case of implicit diffusion the result
160 C is also *only* a prediction.
161 C - We denote "predictors" with an asterisk (*).
162 C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
163 C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
164 C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
165 C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
166 C With implicit diffusion:
167 C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
168 C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
169 C (1 + dt * K * d_zz) theta[n] = theta*
170 C (1 + dt * K * d_zz) salt[n] = salt*
171 C---
172
173 #ifdef ALLOW_AUTODIFF_TAMC
174 C-- dummy statement to end declaration part
175 ikey = 1
176 #endif /* ALLOW_AUTODIFF_TAMC */
177
178 C-- Set up work arrays with valid (i.e. not NaN) values
179 C These inital values do not alter the numerical results. They
180 C just ensure that all memory references are to valid floating
181 C point numbers. This prevents spurious hardware signals due to
182 C uninitialised but inert locations.
183 DO j=1-OLy,sNy+OLy
184 DO i=1-OLx,sNx+OLx
185 xA(i,j) = 0. _d 0
186 yA(i,j) = 0. _d 0
187 uTrans(i,j) = 0. _d 0
188 vTrans(i,j) = 0. _d 0
189 DO k=1,Nr
190 phiHyd(i,j,k) = 0. _d 0
191 KappaRU(i,j,k) = 0. _d 0
192 KappaRV(i,j,k) = 0. _d 0
193 sigmaX(i,j,k) = 0. _d 0
194 sigmaY(i,j,k) = 0. _d 0
195 sigmaR(i,j,k) = 0. _d 0
196 ENDDO
197 rhoKM1 (i,j) = 0. _d 0
198 rhok (i,j) = 0. _d 0
199 phiSurfX(i,j) = 0. _d 0
200 phiSurfY(i,j) = 0. _d 0
201 ENDDO
202 ENDDO
203
204
205 #ifdef ALLOW_AUTODIFF_TAMC
206 C-- HPF directive to help TAMC
207 CHPF$ INDEPENDENT
208 #endif /* ALLOW_AUTODIFF_TAMC */
209
210 DO bj=myByLo(myThid),myByHi(myThid)
211
212 #ifdef ALLOW_AUTODIFF_TAMC
213 C-- HPF directive to help TAMC
214 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
215 CHPF$& ,phiHyd,utrans,vtrans,xA,yA
216 CHPF$& ,KappaRT,KappaRS,KappaRU,KappaRV
217 CHPF$& )
218 #endif /* ALLOW_AUTODIFF_TAMC */
219
220 DO bi=myBxLo(myThid),myBxHi(myThid)
221
222 #ifdef ALLOW_AUTODIFF_TAMC
223 act1 = bi - myBxLo(myThid)
224 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
225
226 act2 = bj - myByLo(myThid)
227 max2 = myByHi(myThid) - myByLo(myThid) + 1
228
229 act3 = myThid - 1
230 max3 = nTx*nTy
231
232 act4 = ikey_dynamics - 1
233
234 ikey = (act1 + 1) + act2*max1
235 & + act3*max1*max2
236 & + act4*max1*max2*max3
237 #endif /* ALLOW_AUTODIFF_TAMC */
238
239 C-- Set up work arrays that need valid initial values
240 DO j=1-OLy,sNy+OLy
241 DO i=1-OLx,sNx+OLx
242 rTrans (i,j) = 0. _d 0
243 fVerT (i,j,1) = 0. _d 0
244 fVerT (i,j,2) = 0. _d 0
245 fVerS (i,j,1) = 0. _d 0
246 fVerS (i,j,2) = 0. _d 0
247 fVerTr1(i,j,1) = 0. _d 0
248 fVerTr1(i,j,2) = 0. _d 0
249 fVerU (i,j,1) = 0. _d 0
250 fVerU (i,j,2) = 0. _d 0
251 fVerV (i,j,1) = 0. _d 0
252 fVerV (i,j,2) = 0. _d 0
253 ENDDO
254 ENDDO
255
256 DO k=1,Nr
257 DO j=1-OLy,sNy+OLy
258 DO i=1-OLx,sNx+OLx
259 C This is currently also used by IVDC and Diagnostics
260 ConvectCount(i,j,k) = 0.
261 KappaRT(i,j,k) = 0. _d 0
262 KappaRS(i,j,k) = 0. _d 0
263 ENDDO
264 ENDDO
265 ENDDO
266
267 iMin = 1-OLx+1
268 iMax = sNx+OLx
269 jMin = 1-OLy+1
270 jMax = sNy+OLy
271
272
273 #ifdef ALLOW_AUTODIFF_TAMC
274 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
275 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
276 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
277 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
278 #ifdef ALLOW_PASSIVE_TRACER
279 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
280 #endif
281 #endif /* ALLOW_AUTODIFF_TAMC */
282
283 C-- Start of diagnostic loop
284 DO k=Nr,1,-1
285
286 #ifdef ALLOW_AUTODIFF_TAMC
287 C? Patrick, is this formula correct now that we change the loop range?
288 C? Do we still need this?
289 cph kkey formula corrected.
290 cph Needed for rhok, rhokm1, in the case useGMREDI.
291 kkey = (ikey-1)*Nr + k
292 CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
293 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
294 #endif /* ALLOW_AUTODIFF_TAMC */
295
296 C-- Integrate continuity vertically for vertical velocity
297 CALL INTEGRATE_FOR_W(
298 I bi, bj, k, uVel, vVel,
299 O wVel,
300 I myThid )
301
302 #ifdef ALLOW_OBCS
303 #ifdef ALLOW_NONHYDROSTATIC
304 C-- Apply OBC to W if in N-H mode
305 IF (useOBCS.AND.nonHydrostatic) THEN
306 CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
307 ENDIF
308 #endif /* ALLOW_NONHYDROSTATIC */
309 #endif /* ALLOW_OBCS */
310
311 C-- Calculate gradients of potential density for isoneutral
312 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
313 c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
314 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
315 #ifdef ALLOW_AUTODIFF_TAMC
316 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
317 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
318 #endif /* ALLOW_AUTODIFF_TAMC */
319 CALL FIND_RHO(
320 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
321 I theta, salt,
322 O rhoK,
323 I myThid )
324 IF (k.GT.1) THEN
325 #ifdef ALLOW_AUTODIFF_TAMC
326 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
327 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
328 #endif /* ALLOW_AUTODIFF_TAMC */
329 CALL FIND_RHO(
330 I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
331 I theta, salt,
332 O rhoKm1,
333 I myThid )
334 ENDIF
335 CALL GRAD_SIGMA(
336 I bi, bj, iMin, iMax, jMin, jMax, k,
337 I rhoK, rhoKm1, rhoK,
338 O sigmaX, sigmaY, sigmaR,
339 I myThid )
340 ENDIF
341
342 C-- Implicit Vertical Diffusion for Convection
343 c ==> should use sigmaR !!!
344 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
345 CALL CALC_IVDC(
346 I bi, bj, iMin, iMax, jMin, jMax, k,
347 I rhoKm1, rhoK,
348 U ConvectCount, KappaRT, KappaRS,
349 I myTime, myIter, myThid)
350 ENDIF
351
352 C-- end of diagnostic k loop (Nr:1)
353 ENDDO
354
355 #ifdef ALLOW_AUTODIFF_TAMC
356 cph avoids recomputation of integrate_for_w
357 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
358 #endif /* ALLOW_AUTODIFF_TAMC */
359
360 #ifdef ALLOW_OBCS
361 C-- Calculate future values on open boundaries
362 IF (useOBCS) THEN
363 CALL OBCS_CALC( bi, bj, myTime+deltaT,
364 I uVel, vVel, wVel, theta, salt,
365 I myThid )
366 ENDIF
367 #endif /* ALLOW_OBCS */
368
369 C-- Determines forcing terms based on external fields
370 C relaxation terms, etc.
371 CALL EXTERNAL_FORCING_SURF(
372 I bi, bj, iMin, iMax, jMin, jMax,
373 I myThid )
374 #ifdef ALLOW_AUTODIFF_TAMC
375 cph needed for KPP
376 CADJ STORE surfacetendencyU(:,:,bi,bj)
377 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
378 CADJ STORE surfacetendencyV(:,:,bi,bj)
379 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
380 CADJ STORE surfacetendencyS(:,:,bi,bj)
381 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
382 CADJ STORE surfacetendencyT(:,:,bi,bj)
383 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
384 #endif /* ALLOW_AUTODIFF_TAMC */
385
386 #ifdef ALLOW_GMREDI
387
388 #ifdef ALLOW_AUTODIFF_TAMC
389 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
390 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
391 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
392 #endif /* ALLOW_AUTODIFF_TAMC */
393 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
394 IF (useGMRedi) THEN
395 DO k=1,Nr
396 CALL GMREDI_CALC_TENSOR(
397 I bi, bj, iMin, iMax, jMin, jMax, k,
398 I sigmaX, sigmaY, sigmaR,
399 I myThid )
400 ENDDO
401 #ifdef ALLOW_AUTODIFF_TAMC
402 ELSE
403 DO k=1, Nr
404 CALL GMREDI_CALC_TENSOR_DUMMY(
405 I bi, bj, iMin, iMax, jMin, jMax, k,
406 I sigmaX, sigmaY, sigmaR,
407 I myThid )
408 ENDDO
409 #endif /* ALLOW_AUTODIFF_TAMC */
410 ENDIF
411
412 #ifdef ALLOW_AUTODIFF_TAMC
413 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
414 CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
415 CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
416 #endif /* ALLOW_AUTODIFF_TAMC */
417
418 #endif /* ALLOW_GMREDI */
419
420 #ifdef ALLOW_KPP
421 C-- Compute KPP mixing coefficients
422 IF (useKPP) THEN
423 CALL KPP_CALC(
424 I bi, bj, myTime, myThid )
425 #ifdef ALLOW_AUTODIFF_TAMC
426 ELSE
427 CALL KPP_CALC_DUMMY(
428 I bi, bj, myTime, myThid )
429 #endif /* ALLOW_AUTODIFF_TAMC */
430 ENDIF
431
432 #ifdef ALLOW_AUTODIFF_TAMC
433 CADJ STORE KPPghat (:,:,:,bi,bj)
434 CADJ & , KPPviscAz (:,:,:,bi,bj)
435 CADJ & , KPPdiffKzT(:,:,:,bi,bj)
436 CADJ & , KPPdiffKzS(:,:,:,bi,bj)
437 CADJ & , KPPfrac (:,: ,bi,bj)
438 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
439 #endif /* ALLOW_AUTODIFF_TAMC */
440
441 #endif /* ALLOW_KPP */
442
443 #ifdef ALLOW_AUTODIFF_TAMC
444 CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key = ikey, byte = isbyte
445 CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key = ikey, byte = isbyte
446 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
447 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
448 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
449 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
450 #ifdef ALLOW_PASSIVE_TRACER
451 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
452 #endif
453 #endif /* ALLOW_AUTODIFF_TAMC */
454
455 #ifdef ALLOW_AIM
456 C AIM - atmospheric intermediate model, physics package code.
457 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
458 IF ( useAIM ) THEN
459 CALL TIMER_START('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
460 CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )
461 CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
462 ENDIF
463 #endif /* ALLOW_AIM */
464
465
466 C-- Start of thermodynamics loop
467 DO k=Nr,1,-1
468 #ifdef ALLOW_AUTODIFF_TAMC
469 C? Patrick Is this formula correct?
470 cph Yes, but I rewrote it.
471 cph Also, the KappaR? need the index and subscript k!
472 kkey = (ikey-1)*Nr + k
473 #endif /* ALLOW_AUTODIFF_TAMC */
474
475 C-- km1 Points to level above k (=k-1)
476 C-- kup Cycles through 1,2 to point to layer above
477 C-- kDown Cycles through 2,1 to point to current layer
478
479 km1 = MAX(1,k-1)
480 kup = 1+MOD(k+1,2)
481 kDown= 1+MOD(k,2)
482
483 iMin = 1-OLx
484 iMax = sNx+OLx
485 jMin = 1-OLy
486 jMax = sNy+OLy
487
488 C-- Get temporary terms used by tendency routines
489 CALL CALC_COMMON_FACTORS (
490 I bi,bj,iMin,iMax,jMin,jMax,k,
491 O xA,yA,uTrans,vTrans,rTrans,maskUp,
492 I myThid)
493
494 #ifdef ALLOW_AUTODIFF_TAMC
495 CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
496 CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
497 #endif /* ALLOW_AUTODIFF_TAMC */
498
499 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
500 C-- Calculate the total vertical diffusivity
501 CALL CALC_DIFFUSIVITY(
502 I bi,bj,iMin,iMax,jMin,jMax,k,
503 I maskUp,
504 O KappaRT,KappaRS,KappaRU,KappaRV,
505 I myThid)
506 #endif
507
508 iMin = 1-OLx+2
509 iMax = sNx+OLx-1
510 jMin = 1-OLy+2
511 jMax = sNy+OLy-1
512
513 C-- Calculate active tracer tendencies (gT,gS,...)
514 C and step forward storing result in gTnm1, gSnm1, etc.
515 IF ( tempStepping ) THEN
516 CALL CALC_GT(
517 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
518 I xA,yA,uTrans,vTrans,rTrans,maskUp,
519 I KappaRT,
520 U fVerT,
521 I myTime, myThid)
522 tauAB = 0.5d0 + abEps
523 CALL TIMESTEP_TRACER(
524 I bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
525 I theta, gT,
526 U gTnm1,
527 I myIter, myThid)
528 ENDIF
529 IF ( saltStepping ) THEN
530 CALL CALC_GS(
531 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
532 I xA,yA,uTrans,vTrans,rTrans,maskUp,
533 I KappaRS,
534 U fVerS,
535 I myTime, myThid)
536 tauAB = 0.5d0 + abEps
537 CALL TIMESTEP_TRACER(
538 I bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
539 I salt, gS,
540 U gSnm1,
541 I myIter, myThid)
542 ENDIF
543 #ifdef ALLOW_PASSIVE_TRACER
544 IF ( tr1Stepping ) THEN
545 CALL CALC_GTR1(
546 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
547 I xA,yA,uTrans,vTrans,rTrans,maskUp,
548 I KappaRT,
549 U fVerTr1,
550 I myTime, myThid)
551 tauAB = 0.5d0 + abEps
552 CALL TIMESTEP_TRACER(
553 I bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
554 I Tr1, gTr1,
555 U gTr1NM1,
556 I myIter, myThid)
557 ENDIF
558 #endif
559
560 #ifdef ALLOW_OBCS
561 C-- Apply open boundary conditions
562 IF (useOBCS) THEN
563 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
564 END IF
565 #endif /* ALLOW_OBCS */
566
567 C-- Freeze water
568 IF (allowFreezing) THEN
569 #ifdef ALLOW_AUTODIFF_TAMC
570 CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
571 CADJ & , key = kkey, byte = isbyte
572 #endif /* ALLOW_AUTODIFF_TAMC */
573 CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
574 END IF
575
576 C-- end of thermodynamic k loop (Nr:1)
577 ENDDO
578
579
580 #ifdef ALLOW_AUTODIFF_TAMC
581 C? Patrick? What about this one?
582 cph Keys iikey and idkey don't seem to be needed
583 cph since storing occurs on different tape for each
584 cph impldiff call anyways.
585 cph Thus, common block comlev1_impl isn't needed either.
586 cph Storing below needed in the case useGMREDI.
587 iikey = (ikey-1)*maximpl
588 #endif /* ALLOW_AUTODIFF_TAMC */
589
590 C-- Implicit diffusion
591 IF (implicitDiffusion) THEN
592
593 IF (tempStepping) THEN
594 #ifdef ALLOW_AUTODIFF_TAMC
595 idkey = iikey + 1
596 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
597 #endif /* ALLOW_AUTODIFF_TAMC */
598 CALL IMPLDIFF(
599 I bi, bj, iMin, iMax, jMin, jMax,
600 I deltaTtracer, KappaRT, recip_HFacC,
601 U gTNm1,
602 I myThid )
603 ENDIF
604
605 IF (saltStepping) THEN
606 #ifdef ALLOW_AUTODIFF_TAMC
607 idkey = iikey + 2
608 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
609 #endif /* ALLOW_AUTODIFF_TAMC */
610 CALL IMPLDIFF(
611 I bi, bj, iMin, iMax, jMin, jMax,
612 I deltaTtracer, KappaRS, recip_HFacC,
613 U gSNm1,
614 I myThid )
615 ENDIF
616
617 #ifdef ALLOW_PASSIVE_TRACER
618 IF (tr1Stepping) THEN
619 #ifdef ALLOW_AUTODIFF_TAMC
620 CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
621 #endif /* ALLOW_AUTODIFF_TAMC */
622 CALL IMPLDIFF(
623 I bi, bj, iMin, iMax, jMin, jMax,
624 I deltaTtracer, KappaRT, recip_HFacC,
625 U gTr1Nm1,
626 I myThid )
627 ENDIF
628 #endif
629
630 #ifdef ALLOW_OBCS
631 C-- Apply open boundary conditions
632 IF (useOBCS) THEN
633 DO K=1,Nr
634 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
635 ENDDO
636 END IF
637 #endif /* ALLOW_OBCS */
638
639 C-- End If implicitDiffusion
640 ENDIF
641
642 C-- Start computation of dynamics
643 iMin = 1-OLx+2
644 iMax = sNx+OLx-1
645 jMin = 1-OLy+2
646 jMax = sNy+OLy-1
647
648 C-- Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
649 C (note: this loop will be replaced by CALL CALC_GRAD_ETA)
650 IF (implicSurfPress.NE.1.) THEN
651 CALL CALC_GRAD_PHI_SURF(
652 I bi,bj,iMin,iMax,jMin,jMax,
653 I etaN,
654 O phiSurfX,phiSurfY,
655 I myThid )
656 ENDIF
657
658 C-- Start of dynamics loop
659 DO k=1,Nr
660
661 C-- km1 Points to level above k (=k-1)
662 C-- kup Cycles through 1,2 to point to layer above
663 C-- kDown Cycles through 2,1 to point to current layer
664
665 km1 = MAX(1,k-1)
666 kup = 1+MOD(k+1,2)
667 kDown= 1+MOD(k,2)
668
669 C-- Integrate hydrostatic balance for phiHyd with BC of
670 C phiHyd(z=0)=0
671 C distinguishe between Stagger and Non Stagger time stepping
672 IF (staggerTimeStep) THEN
673 CALL CALC_PHI_HYD(
674 I bi,bj,iMin,iMax,jMin,jMax,k,
675 I gTnm1, gSnm1,
676 U phiHyd,
677 I myThid )
678 ELSE
679 CALL CALC_PHI_HYD(
680 I bi,bj,iMin,iMax,jMin,jMax,k,
681 I theta, salt,
682 U phiHyd,
683 I myThid )
684 ENDIF
685
686 C-- Calculate accelerations in the momentum equations (gU, gV, ...)
687 C and step forward storing the result in gUnm1, gVnm1, etc...
688 IF ( momStepping ) THEN
689 CALL CALC_MOM_RHS(
690 I bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
691 I phiHyd,KappaRU,KappaRV,
692 U fVerU, fVerV,
693 I myTime, myThid)
694 CALL TIMESTEP(
695 I bi,bj,iMin,iMax,jMin,jMax,k,
696 I phiHyd, phiSurfX, phiSurfY,
697 I myIter, myThid)
698
699 #ifdef ALLOW_OBCS
700 C-- Apply open boundary conditions
701 IF (useOBCS) THEN
702 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
703 END IF
704 #endif /* ALLOW_OBCS */
705
706 #ifdef ALLOW_AUTODIFF_TAMC
707 #ifdef INCLUDE_CD_CODE
708 ELSE
709 DO j=1-OLy,sNy+OLy
710 DO i=1-OLx,sNx+OLx
711 guCD(i,j,k,bi,bj) = 0.0
712 gvCD(i,j,k,bi,bj) = 0.0
713 END DO
714 END DO
715 #endif /* INCLUDE_CD_CODE */
716 #endif /* ALLOW_AUTODIFF_TAMC */
717 ENDIF
718
719
720 C-- end of dynamics k loop (1:Nr)
721 ENDDO
722
723
724
725 C-- Implicit viscosity
726 IF (implicitViscosity.AND.momStepping) THEN
727 #ifdef ALLOW_AUTODIFF_TAMC
728 idkey = iikey + 3
729 CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
730 #endif /* ALLOW_AUTODIFF_TAMC */
731 CALL IMPLDIFF(
732 I bi, bj, iMin, iMax, jMin, jMax,
733 I deltaTmom, KappaRU,recip_HFacW,
734 U gUNm1,
735 I myThid )
736 #ifdef ALLOW_AUTODIFF_TAMC
737 idkey = iikey + 4
738 CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
739 #endif /* ALLOW_AUTODIFF_TAMC */
740 CALL IMPLDIFF(
741 I bi, bj, iMin, iMax, jMin, jMax,
742 I deltaTmom, KappaRV,recip_HFacS,
743 U gVNm1,
744 I myThid )
745
746 #ifdef ALLOW_OBCS
747 C-- Apply open boundary conditions
748 IF (useOBCS) THEN
749 DO K=1,Nr
750 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
751 ENDDO
752 END IF
753 #endif /* ALLOW_OBCS */
754
755 #ifdef INCLUDE_CD_CODE
756 #ifdef ALLOW_AUTODIFF_TAMC
757 idkey = iikey + 5
758 CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
759 #endif /* ALLOW_AUTODIFF_TAMC */
760 CALL IMPLDIFF(
761 I bi, bj, iMin, iMax, jMin, jMax,
762 I deltaTmom, KappaRU,recip_HFacW,
763 U vVelD,
764 I myThid )
765 #ifdef ALLOW_AUTODIFF_TAMC
766 idkey = iikey + 6
767 CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
768 #endif /* ALLOW_AUTODIFF_TAMC */
769 CALL IMPLDIFF(
770 I bi, bj, iMin, iMax, jMin, jMax,
771 I deltaTmom, KappaRV,recip_HFacS,
772 U uVelD,
773 I myThid )
774 #endif /* INCLUDE_CD_CODE */
775 C-- End If implicitViscosity.AND.momStepping
776 ENDIF
777
778 Cjmc : add for phiHyd output <- but not working if multi tile per CPU
779 c IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
780 c & .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
781 c WRITE(suff,'(I10.10)') myIter+1
782 c CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
783 c ENDIF
784 Cjmc(end)
785
786 #ifdef ALLOW_TIMEAVE
787 IF (taveFreq.GT.0.) THEN
788 CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
789 I deltaTclock, bi, bj, myThid)
790 IF (ivdc_kappa.NE.0.) THEN
791 CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
792 I deltaTclock, bi, bj, myThid)
793 ENDIF
794 ENDIF
795 #endif /* ALLOW_TIMEAVE */
796
797 ENDDO
798 ENDDO
799
800 #ifndef EXCLUDE_DEBUGMODE
801 If (debugMode) THEN
802 CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
803 CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
804 CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
805 CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
806 CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
807 CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
808 CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
809 CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
810 CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
811 CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
812 CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
813 CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
814 CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
815 CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
816 ENDIF
817 #endif
818
819 RETURN
820 END

  ViewVC Help
Powered by ViewVC 1.1.22