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

Annotation 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 - (hide annotations) (download)
Tue Aug 21 15:40:11 2001 UTC (22 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
FILE REMOVED
no longer needed

1 jmc 1.2 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 cnh 1.1 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