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

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

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

Revision 1.69 - (hide annotations) (download)
Wed Jun 6 14:55:45 2001 UTC (22 years, 11 months ago) by adcroft
Branch: MAIN
Changes since 1.68: +17 -1 lines 
Added a debugMode that uses same statistics stuff as monitor.F
Can be disabled with -DEXCLUDE_DEBUGMODE. Turn on at run-time
with debugMode=.true.  Default is enabled but off.
1 adcroft 1.69 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.68 2001/05/29 14:01:37 adcroft Exp $
2 adcroft 1.68 C $Name: $
3 cnh 1.1
4 adcroft 1.24 #include "CPP_OPTIONS.h"
5 cnh 1.1
6 cnh 1.8 SUBROUTINE DYNAMICS(myTime, myIter, myThid)
7 cnh 1.1 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 adcroft 1.40 IMPLICIT NONE
25 cnh 1.1
26     C == Global variables ===
27     #include "SIZE.h"
28     #include "EEPARAMS.h"
29 adcroft 1.6 #include "PARAMS.h"
30 adcroft 1.3 #include "DYNVARS.h"
31 adcroft 1.42 #include "GRID.h"
32 heimbach 1.49
33     #ifdef ALLOW_AUTODIFF_TAMC
34 heimbach 1.53 # include "tamc.h"
35     # include "tamc_keys.h"
36 heimbach 1.67 # include "FFIELDS.h"
37     # ifdef ALLOW_KPP
38     # include "KPP.h"
39     # endif
40     # ifdef ALLOW_GMREDI
41     # include "GMREDI.h"
42     # endif
43 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
44 heimbach 1.49
45 jmc 1.64 #ifdef ALLOW_TIMEAVE
46     #include "TIMEAVE_STATV.h"
47 jmc 1.62 #endif
48    
49 cnh 1.1 C == Routine arguments ==
50 cnh 1.8 C myTime - Current time in simulation
51     C myIter - Current iteration number in simulation
52 cnh 1.1 C myThid - Thread number for this instance of the routine.
53 cnh 1.8 _RL myTime
54     INTEGER myIter
55 adcroft 1.47 INTEGER myThid
56 cnh 1.1
57     C == Local variables
58     C xA, yA - Per block temporaries holding face areas
59 cnh 1.38 C uTrans, vTrans, rTrans - Per block temporaries holding flow
60     C transport
61 jmc 1.61 C o uTrans: Zonal transport
62 cnh 1.1 C o vTrans: Meridional transport
63 cnh 1.30 C o rTrans: Vertical transport
64 adcroft 1.68 C maskUp o maskUp: land/water mask for W points
65 adcroft 1.58 C fVer[STUV] o fVer: Vertical flux term - note fVer
66 cnh 1.1 C is "pipelined" in the vertical
67     C so we need an fVer for each
68     C variable.
69 adcroft 1.58 C rhoK, rhoKM1 - Density at current level, and level above
70 cnh 1.31 C phiHyd - Hydrostatic part of the potential phiHydi.
71 cnh 1.38 C In z coords phiHydiHyd is the hydrostatic
72 jmc 1.65 C Potential (=pressure/rho0) anomaly
73 cnh 1.38 C In p coords phiHydiHyd is the geopotential
74 jmc 1.65 C surface height anomaly.
75 jmc 1.63 C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
76     C phiSurfY or geopotentiel (atmos) in X and Y direction
77 cnh 1.30 C KappaRT, - Total diffusion in vertical for T and S.
78 cnh 1.38 C KappaRS (background + spatially varying, isopycnal term).
79 cnh 1.30 C iMin, iMax - Ranges and sub-block indices on which calculations
80     C jMin, jMax are applied.
81 cnh 1.1 C bi, bj
82 heimbach 1.53 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 cnh 1.38 C index into fVerTerm.
85 adcroft 1.68 C tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.
86 cnh 1.30 _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 cnh 1.31 _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97 cnh 1.30 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99 jmc 1.63 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100     _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101 cnh 1.31 _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
102     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
103 adcroft 1.42 _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
104     _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
105 adcroft 1.50 _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 adcroft 1.68 _RL tauAB
109 adcroft 1.12
110 jmc 1.62 C This is currently used by IVDC and Diagnostics
111 adcroft 1.45 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112    
113 cnh 1.1 INTEGER iMin, iMax
114     INTEGER jMin, jMax
115     INTEGER bi, bj
116     INTEGER i, j
117 heimbach 1.53 INTEGER k, km1, kup, kDown
118 cnh 1.1
119 jmc 1.62 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 adcroft 1.11 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 heimbach 1.53 C physics, parameterizations etc...) are calculated
139 adcroft 1.11 C rho = rho ( theta[n], salt[n] )
140 cnh 1.27 C b = b(rho, theta)
141 adcroft 1.11 C K31 = K31 ( rho )
142 jmc 1.61 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 adcroft 1.11 C
147 adcroft 1.12 C "Time-stepping" or "Prediction"
148 adcroft 1.11 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 adcroft 1.12 C With implicit diffusion:
163 adcroft 1.11 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 adcroft 1.12 C (1 + dt * K * d_zz) theta[n] = theta*
166     C (1 + dt * K * d_zz) salt[n] = salt*
167 adcroft 1.11 C---
168    
169 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
170     C-- dummy statement to end declaration part
171     ikey = 1
172 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
173 heimbach 1.49
174 cnh 1.1 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 adcroft 1.5 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 heimbach 1.53 DO k=1,Nr
186 adcroft 1.58 phiHyd(i,j,k) = 0. _d 0
187 adcroft 1.45 KappaRU(i,j,k) = 0. _d 0
188     KappaRV(i,j,k) = 0. _d 0
189 adcroft 1.50 sigmaX(i,j,k) = 0. _d 0
190     sigmaY(i,j,k) = 0. _d 0
191     sigmaR(i,j,k) = 0. _d 0
192 cnh 1.1 ENDDO
193 cnh 1.30 rhoKM1 (i,j) = 0. _d 0
194     rhok (i,j) = 0. _d 0
195 jmc 1.63 phiSurfX(i,j) = 0. _d 0
196     phiSurfY(i,j) = 0. _d 0
197 cnh 1.1 ENDDO
198     ENDDO
199    
200 cnh 1.35
201 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
202     C-- HPF directive to help TAMC
203 heimbach 1.53 CHPF$ INDEPENDENT
204     #endif /* ALLOW_AUTODIFF_TAMC */
205 heimbach 1.49
206 cnh 1.1 DO bj=myByLo(myThid),myByHi(myThid)
207 heimbach 1.49
208     #ifdef ALLOW_AUTODIFF_TAMC
209     C-- HPF directive to help TAMC
210 jmc 1.61 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
211 adcroft 1.68 CHPF$& ,phiHyd,utrans,vtrans,xA,yA
212 heimbach 1.53 CHPF$& ,KappaRT,KappaRS,KappaRU,KappaRV
213     CHPF$& )
214     #endif /* ALLOW_AUTODIFF_TAMC */
215 heimbach 1.49
216 cnh 1.1 DO bi=myBxLo(myThid),myBxHi(myThid)
217    
218 heimbach 1.49 #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 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
234 heimbach 1.49
235 cnh 1.7 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 cnh 1.27 rTrans(i,j) = 0. _d 0
239 cnh 1.30 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 cnh 1.7 ENDDO
248     ENDDO
249    
250 adcroft 1.45 DO k=1,Nr
251     DO j=1-OLy,sNy+OLy
252     DO i=1-OLx,sNx+OLx
253 jmc 1.62 C This is currently also used by IVDC and Diagnostics
254 adcroft 1.45 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 cnh 1.1 iMin = 1-OLx+1
262     iMax = sNx+OLx
263     jMin = 1-OLy+1
264     jMax = sNy+OLy
265 cnh 1.35
266 adcroft 1.5
267 heimbach 1.66 #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 heimbach 1.67 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
271 heimbach 1.66 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
272     #endif /* ALLOW_AUTODIFF_TAMC */
273    
274 adcroft 1.58 C-- Start of diagnostic loop
275     DO k=Nr,1,-1
276 heimbach 1.49
277     #ifdef ALLOW_AUTODIFF_TAMC
278 adcroft 1.58 C? Patrick, is this formula correct now that we change the loop range?
279     C? Do we still need this?
280 heimbach 1.66 cph kkey formula corrected.
281     cph Needed for rhok, rhokm1, in the case useGMREDI.
282     kkey = (ikey-1)*Nr + k
283 heimbach 1.67 CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
284     CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
285 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
286 heimbach 1.49
287 adcroft 1.58 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 adcroft 1.60 C-- Apply OBC to W if in N-H mode
296 adcroft 1.58 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 heimbach 1.67 #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 adcroft 1.58 CALL FIND_RHO(
311     I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
312     I theta, salt,
313     O rhoK,
314 cnh 1.30 I myThid )
315 heimbach 1.67 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 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
322 adcroft 1.58 I theta, salt,
323     O rhoKm1,
324 cnh 1.30 I myThid )
325 heimbach 1.67 ENDIF
326 adcroft 1.58 CALL GRAD_SIGMA(
327 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k,
328 adcroft 1.58 I rhoK, rhoKm1, rhoK,
329 adcroft 1.50 O sigmaX, sigmaY, sigmaR,
330     I myThid )
331 adcroft 1.58 ENDIF
332 heimbach 1.49
333 adcroft 1.58 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 jmc 1.62 ENDIF
342 heimbach 1.53
343 adcroft 1.58 C-- end of diagnostic k loop (Nr:1)
344 heimbach 1.49 ENDDO
345    
346 heimbach 1.67 #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 adcroft 1.58 #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 heimbach 1.54 I bi, bj, iMin, iMax, jMin, jMax,
364     I myThid )
365 heimbach 1.67 #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 heimbach 1.54
377 adcroft 1.58 #ifdef ALLOW_GMREDI
378 heimbach 1.67
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 adcroft 1.58 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
385 heimbach 1.53 IF (useGMRedi) THEN
386 adcroft 1.58 DO k=1,Nr
387 heimbach 1.53 CALL GMREDI_CALC_TENSOR(
388     I bi, bj, iMin, iMax, jMin, jMax, k,
389 adcroft 1.50 I sigmaX, sigmaY, sigmaR,
390     I myThid )
391 heimbach 1.53 ENDDO
392 heimbach 1.55 #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 heimbach 1.53 ENDIF
402 heimbach 1.67
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 adcroft 1.58 #endif /* ALLOW_GMREDI */
410 heimbach 1.53
411 adcroft 1.58 #ifdef ALLOW_KPP
412     C-- Compute KPP mixing coefficients
413 heimbach 1.53 IF (useKPP) THEN
414     CALL KPP_CALC(
415 heimbach 1.54 I bi, bj, myTime, myThid )
416 heimbach 1.66 #ifdef ALLOW_AUTODIFF_TAMC
417     ELSE
418 heimbach 1.67 CALL KPP_CALC_DUMMY(
419     I bi, bj, myTime, myThid )
420 heimbach 1.66 #endif /* ALLOW_AUTODIFF_TAMC */
421 adcroft 1.58 ENDIF
422 heimbach 1.66
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 adcroft 1.58 #endif /* ALLOW_KPP */
433 heimbach 1.53
434     #ifdef ALLOW_AUTODIFF_TAMC
435 adcroft 1.58 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, myTime, myThid )
449     CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
450 heimbach 1.53 ENDIF
451 adcroft 1.58 #endif /* ALLOW_AIM */
452    
453 heimbach 1.53
454 adcroft 1.58 C-- Start of thermodynamics loop
455     DO k=Nr,1,-1
456 heimbach 1.67 #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 adcroft 1.58
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 cnh 1.1
476     C-- Get temporary terms used by tendency routines
477     CALL CALC_COMMON_FACTORS (
478 adcroft 1.68 I bi,bj,iMin,iMax,jMin,jMax,k,
479     O xA,yA,uTrans,vTrans,rTrans,maskUp,
480 cnh 1.1 I myThid)
481 heimbach 1.67
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 heimbach 1.49
487 cnh 1.38 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
488 adcroft 1.12 C-- Calculate the total vertical diffusivity
489     CALL CALC_DIFFUSIVITY(
490 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,
491 adcroft 1.68 I maskUp,
492 adcroft 1.42 O KappaRT,KappaRS,KappaRU,KappaRV,
493 adcroft 1.12 I myThid)
494 cnh 1.38 #endif
495 adcroft 1.58
496     C-- Calculate active tracer tendencies (gT,gS,...)
497     C and step forward storing result in gTnm1, gSnm1, etc.
498 cnh 1.9 IF ( tempStepping ) THEN
499 adcroft 1.58 CALL CALC_GT(
500 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
501 adcroft 1.68 I xA,yA,uTrans,vTrans,rTrans,maskUp,
502 adcroft 1.50 I KappaRT,
503 adcroft 1.58 U fVerT,
504 cnh 1.37 I myTime, myThid)
505 adcroft 1.68 tauAB = 0.5d0 + abEps
506 adcroft 1.58 CALL TIMESTEP_TRACER(
507 adcroft 1.68 I bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
508 adcroft 1.58 I theta, gT,
509     U gTnm1,
510     I myIter, myThid)
511 cnh 1.9 ENDIF
512 adcroft 1.18 IF ( saltStepping ) THEN
513 adcroft 1.58 CALL CALC_GS(
514 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
515 adcroft 1.68 I xA,yA,uTrans,vTrans,rTrans,maskUp,
516 adcroft 1.50 I KappaRS,
517 adcroft 1.58 U fVerS,
518 cnh 1.37 I myTime, myThid)
519 adcroft 1.68 tauAB = 0.5d0 + abEps
520 adcroft 1.58 CALL TIMESTEP_TRACER(
521 adcroft 1.68 I bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
522 adcroft 1.58 I salt, gS,
523     U gSnm1,
524     I myIter, myThid)
525 adcroft 1.18 ENDIF
526 adcroft 1.58
527     #ifdef ALLOW_OBCS
528 adcroft 1.41 C-- Apply open boundary conditions
529 adcroft 1.58 IF (useOBCS) THEN
530     CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
531     END IF
532     #endif /* ALLOW_OBCS */
533 heimbach 1.54
534 adcroft 1.41 C-- Freeze water
535 heimbach 1.49 IF (allowFreezing) THEN
536     #ifdef ALLOW_AUTODIFF_TAMC
537 adcroft 1.58 CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
538     CADJ & , key = kkey, byte = isbyte
539 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
540     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
541 heimbach 1.49 END IF
542 adcroft 1.48
543 adcroft 1.58 C-- end of thermodynamic k loop (Nr:1)
544     ENDDO
545 adcroft 1.45
546 adcroft 1.11
547 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
548 heimbach 1.66 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 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
556 heimbach 1.51
557     C-- Implicit diffusion
558     IF (implicitDiffusion) THEN
559 heimbach 1.49
560 jmc 1.62 IF (tempStepping) THEN
561 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
562     idkey = iikey + 1
563 heimbach 1.66 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
564 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
565 heimbach 1.49 CALL IMPLDIFF(
566 adcroft 1.42 I bi, bj, iMin, iMax, jMin, jMax,
567 adcroft 1.58 I deltaTtracer, KappaRT, recip_HFacC,
568 adcroft 1.42 U gTNm1,
569     I myThid )
570 adcroft 1.58 ENDIF
571 heimbach 1.49
572     IF (saltStepping) THEN
573     #ifdef ALLOW_AUTODIFF_TAMC
574     idkey = iikey + 2
575 heimbach 1.66 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
576 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
577 heimbach 1.49 CALL IMPLDIFF(
578 adcroft 1.42 I bi, bj, iMin, iMax, jMin, jMax,
579 adcroft 1.58 I deltaTtracer, KappaRS, recip_HFacC,
580 adcroft 1.42 U gSNm1,
581     I myThid )
582 adcroft 1.58 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 heimbach 1.49 END IF
591 adcroft 1.58 #endif /* ALLOW_OBCS */
592 heimbach 1.49
593 adcroft 1.58 C-- End If implicitDiffusion
594 heimbach 1.53 ENDIF
595 heimbach 1.49
596 jmc 1.63 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 jmc 1.65 C-- Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
603 jmc 1.63 C (note: this loop will be replaced by CALL CALC_GRAD_ETA)
604     IF (implicSurfPress.NE.1.) THEN
605 jmc 1.65 CALL CALC_GRAD_PHI_SURF(
606     I bi,bj,iMin,iMax,jMin,jMax,
607     I etaN,
608     O phiSurfX,phiSurfY,
609     I myThid )
610 jmc 1.63 ENDIF
611 adcroft 1.58
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 jmc 1.63 I bi,bj,iMin,iMax,jMin,jMax,k,
650     I phiHyd, phiSurfX, phiSurfY,
651 adcroft 1.58 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 adcroft 1.44 C-- Implicit viscosity
680 adcroft 1.58 IF (implicitViscosity.AND.momStepping) THEN
681     #ifdef ALLOW_AUTODIFF_TAMC
682     idkey = iikey + 3
683 heimbach 1.66 CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
684 adcroft 1.58 #endif /* ALLOW_AUTODIFF_TAMC */
685 adcroft 1.42 CALL IMPLDIFF(
686     I bi, bj, iMin, iMax, jMin, jMax,
687     I deltaTmom, KappaRU,recip_HFacW,
688     U gUNm1,
689     I myThid )
690 adcroft 1.58 #ifdef ALLOW_AUTODIFF_TAMC
691     idkey = iikey + 4
692 heimbach 1.66 CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
693 adcroft 1.58 #endif /* ALLOW_AUTODIFF_TAMC */
694 adcroft 1.42 CALL IMPLDIFF(
695     I bi, bj, iMin, iMax, jMin, jMax,
696     I deltaTmom, KappaRV,recip_HFacS,
697     U gVNm1,
698     I myThid )
699 heimbach 1.49
700 adcroft 1.58 #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 heimbach 1.49
709 adcroft 1.58 #ifdef INCLUDE_CD_CODE
710     #ifdef ALLOW_AUTODIFF_TAMC
711     idkey = iikey + 5
712 heimbach 1.66 CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
713 adcroft 1.58 #endif /* ALLOW_AUTODIFF_TAMC */
714 adcroft 1.42 CALL IMPLDIFF(
715     I bi, bj, iMin, iMax, jMin, jMax,
716     I deltaTmom, KappaRU,recip_HFacW,
717     U vVelD,
718     I myThid )
719 adcroft 1.58 #ifdef ALLOW_AUTODIFF_TAMC
720     idkey = iikey + 6
721 heimbach 1.66 CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
722 adcroft 1.58 #endif /* ALLOW_AUTODIFF_TAMC */
723 adcroft 1.42 CALL IMPLDIFF(
724     I bi, bj, iMin, iMax, jMin, jMax,
725     I deltaTmom, KappaRV,recip_HFacS,
726     U uVelD,
727     I myThid )
728 adcroft 1.58 #endif /* INCLUDE_CD_CODE */
729     C-- End If implicitViscosity.AND.momStepping
730 heimbach 1.53 ENDIF
731 cnh 1.1
732 jmc 1.62 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 jmc 1.64 #ifdef ALLOW_TIMEAVE
741 jmc 1.62 IF (taveFreq.GT.0.) THEN
742 adcroft 1.68 CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
743 jmc 1.64 I deltaTclock, bi, bj, myThid)
744 jmc 1.62 IF (ivdc_kappa.NE.0.) THEN
745 jmc 1.64 CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
746     I deltaTclock, bi, bj, myThid)
747 jmc 1.62 ENDIF
748     ENDIF
749 jmc 1.64 #endif /* ALLOW_TIMEAVE */
750 jmc 1.62
751 cnh 1.1 ENDDO
752     ENDDO
753 adcroft 1.69
754     #ifndef EXCLUDE_DEBUGMODE
755     CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
756     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
757     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
758     CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
759     CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
760     CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
761     CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
762     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
763     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
764     CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
765     CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
766     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
767     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
768     #endif
769 cnh 1.1
770     RETURN
771     END
  ViewVC Help
Powered by ViewVC 1.1.22