/[MITgcm]/MITgcm/verification/aim.5l_cs/code/dynamics.F
ViewVC logotype

Contents of /MITgcm/verification/aim.5l_cs/code/dynamics.F

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


Revision 1.2 - (show annotations) (download)
Tue Aug 21 15:40:11 2001 UTC (22 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
FILE REMOVED
no longer needed

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

  ViewVC Help
Powered by ViewVC 1.1.22