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

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

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


Revision 1.22 - (hide annotations) (download)
Thu May 30 02:31:01 2002 UTC (22 years ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint45b_post, checkpoint45c_post
Changes since 1.21: +18 -1 lines
Included CPP option SINGLE_LAYER_MODE
to configure barotropic setup (Martin Losch).

1 heimbach 1.22 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.21 2002/04/30 22:53:17 jmc Exp $
2 adcroft 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5 jmc 1.21 #ifdef ALLOW_AUTODIFF_TAMC
6     # ifdef ALLOW_GMREDI
7     # include "GMREDI_OPTIONS.h"
8     # endif
9     # ifdef ALLOW_KPP
10     # include "KPP_OPTIONS.h"
11     # endif
12     #endif /* ALLOW_AUTODIFF_TAMC */
13 adcroft 1.1
14 cnh 1.9 CBOP
15     C !ROUTINE: THERMODYNAMICS
16     C !INTERFACE:
17 adcroft 1.1 SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
18 cnh 1.9 C !DESCRIPTION: \bv
19     C *==========================================================*
20     C | SUBROUTINE THERMODYNAMICS
21     C | o Controlling routine for the prognostic part of the
22     C | thermo-dynamics.
23     C *===========================================================
24     C | The algorithm...
25     C |
26     C | "Correction Step"
27     C | =================
28     C | Here we update the horizontal velocities with the surface
29     C | pressure such that the resulting flow is either consistent
30     C | with the free-surface evolution or the rigid-lid:
31     C | U[n] = U* + dt x d/dx P
32     C | V[n] = V* + dt x d/dy P
33     C |
34     C | "Calculation of Gs"
35     C | ===================
36     C | This is where all the accelerations and tendencies (ie.
37     C | physics, parameterizations etc...) are calculated
38     C | rho = rho ( theta[n], salt[n] )
39     C | b = b(rho, theta)
40     C | K31 = K31 ( rho )
41     C | Gu[n] = Gu( u[n], v[n], wVel, b, ... )
42     C | Gv[n] = Gv( u[n], v[n], wVel, b, ... )
43     C | Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
44     C | Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
45     C |
46     C | "Time-stepping" or "Prediction"
47     C | ================================
48     C | The models variables are stepped forward with the appropriate
49     C | time-stepping scheme (currently we use Adams-Bashforth II)
50     C | - For momentum, the result is always *only* a "prediction"
51     C | in that the flow may be divergent and will be "corrected"
52     C | later with a surface pressure gradient.
53     C | - Normally for tracers the result is the new field at time
54     C | level [n+1} *BUT* in the case of implicit diffusion the result
55     C | is also *only* a prediction.
56     C | - We denote "predictors" with an asterisk (*).
57     C | U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
58     C | V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
59     C | theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
60     C | salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
61     C | With implicit diffusion:
62     C | theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63     C | salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
64     C | (1 + dt * K * d_zz) theta[n] = theta*
65     C | (1 + dt * K * d_zz) salt[n] = salt*
66     C |
67     C *==========================================================*
68     C \ev
69    
70     C !USES:
71 adcroft 1.1 IMPLICIT NONE
72     C == Global variables ===
73     #include "SIZE.h"
74     #include "EEPARAMS.h"
75     #include "PARAMS.h"
76     #include "DYNVARS.h"
77     #include "GRID.h"
78 adcroft 1.4 #include "GAD.h"
79 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
80     #include "TR1.h"
81     #endif
82     #ifdef ALLOW_AUTODIFF_TAMC
83     # include "tamc.h"
84     # include "tamc_keys.h"
85     # include "FFIELDS.h"
86     # ifdef ALLOW_KPP
87     # include "KPP.h"
88     # endif
89     # ifdef ALLOW_GMREDI
90     # include "GMREDI.h"
91     # endif
92     #endif /* ALLOW_AUTODIFF_TAMC */
93     #ifdef ALLOW_TIMEAVE
94     #include "TIMEAVE_STATV.h"
95     #endif
96    
97 cnh 1.9 C !INPUT/OUTPUT PARAMETERS:
98 adcroft 1.1 C == Routine arguments ==
99     C myTime - Current time in simulation
100     C myIter - Current iteration number in simulation
101     C myThid - Thread number for this instance of the routine.
102     _RL myTime
103     INTEGER myIter
104     INTEGER myThid
105    
106 cnh 1.9 C !LOCAL VARIABLES:
107 adcroft 1.1 C == Local variables
108     C xA, yA - Per block temporaries holding face areas
109     C uTrans, vTrans, rTrans - Per block temporaries holding flow
110     C transport
111     C o uTrans: Zonal transport
112     C o vTrans: Meridional transport
113     C o rTrans: Vertical transport
114     C maskUp o maskUp: land/water mask for W points
115     C fVer[STUV] o fVer: Vertical flux term - note fVer
116     C is "pipelined" in the vertical
117     C so we need an fVer for each
118     C variable.
119     C rhoK, rhoKM1 - Density at current level, and level above
120     C phiHyd - Hydrostatic part of the potential phiHydi.
121     C In z coords phiHydiHyd is the hydrostatic
122     C Potential (=pressure/rho0) anomaly
123     C In p coords phiHydiHyd is the geopotential
124     C surface height anomaly.
125     C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
126     C phiSurfY or geopotentiel (atmos) in X and Y direction
127     C KappaRT, - Total diffusion in vertical for T and S.
128     C KappaRS (background + spatially varying, isopycnal term).
129     C iMin, iMax - Ranges and sub-block indices on which calculations
130     C jMin, jMax are applied.
131     C bi, bj
132     C k, kup, - Index for layer above and below. kup and kDown
133     C kDown, km1 are switched with layer to be the appropriate
134     C index into fVerTerm.
135     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140     _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
142     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
143     _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
144     _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
145     _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147     _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148     _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149     _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
150     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
151     _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
152     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
153     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
154 cnh 1.9 C This is currently used by IVDC and Diagnostics
155 adcroft 1.1 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
156     INTEGER iMin, iMax
157     INTEGER jMin, jMax
158     INTEGER bi, bj
159     INTEGER i, j
160     INTEGER k, km1, kup, kDown
161    
162 cnh 1.9 CEOP
163 adcroft 1.1
164     #ifdef ALLOW_AUTODIFF_TAMC
165     C-- dummy statement to end declaration part
166     ikey = 1
167     #endif /* ALLOW_AUTODIFF_TAMC */
168    
169     C-- Set up work arrays with valid (i.e. not NaN) values
170     C These inital values do not alter the numerical results. They
171     C just ensure that all memory references are to valid floating
172     C point numbers. This prevents spurious hardware signals due to
173     C uninitialised but inert locations.
174     DO j=1-OLy,sNy+OLy
175     DO i=1-OLx,sNx+OLx
176     xA(i,j) = 0. _d 0
177     yA(i,j) = 0. _d 0
178     uTrans(i,j) = 0. _d 0
179     vTrans(i,j) = 0. _d 0
180     rhok (i,j) = 0. _d 0
181     phiSurfX(i,j) = 0. _d 0
182     phiSurfY(i,j) = 0. _d 0
183     ENDDO
184     ENDDO
185    
186    
187     #ifdef ALLOW_AUTODIFF_TAMC
188     C-- HPF directive to help TAMC
189     CHPF$ INDEPENDENT
190     #endif /* ALLOW_AUTODIFF_TAMC */
191    
192     DO bj=myByLo(myThid),myByHi(myThid)
193    
194     #ifdef ALLOW_AUTODIFF_TAMC
195     C-- HPF directive to help TAMC
196 heimbach 1.2 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
197 adcroft 1.1 CHPF$& ,phiHyd,utrans,vtrans,xA,yA
198 heimbach 1.2 CHPF$& ,KappaRT,KappaRS
199 adcroft 1.1 CHPF$& )
200     #endif /* ALLOW_AUTODIFF_TAMC */
201    
202     DO bi=myBxLo(myThid),myBxHi(myThid)
203    
204     #ifdef ALLOW_AUTODIFF_TAMC
205     act1 = bi - myBxLo(myThid)
206     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
207     act2 = bj - myByLo(myThid)
208     max2 = myByHi(myThid) - myByLo(myThid) + 1
209     act3 = myThid - 1
210     max3 = nTx*nTy
211     act4 = ikey_dynamics - 1
212     ikey = (act1 + 1) + act2*max1
213     & + act3*max1*max2
214     & + act4*max1*max2*max3
215     #endif /* ALLOW_AUTODIFF_TAMC */
216    
217     C-- Set up work arrays that need valid initial values
218     DO j=1-OLy,sNy+OLy
219     DO i=1-OLx,sNx+OLx
220     rTrans (i,j) = 0. _d 0
221     fVerT (i,j,1) = 0. _d 0
222     fVerT (i,j,2) = 0. _d 0
223     fVerS (i,j,1) = 0. _d 0
224     fVerS (i,j,2) = 0. _d 0
225     fVerTr1(i,j,1) = 0. _d 0
226     fVerTr1(i,j,2) = 0. _d 0
227 heimbach 1.20 rhoKM1 (i,j) = 0. _d 0
228 adcroft 1.1 ENDDO
229     ENDDO
230    
231     DO k=1,Nr
232     DO j=1-OLy,sNy+OLy
233     DO i=1-OLx,sNx+OLx
234     C This is currently also used by IVDC and Diagnostics
235 heimbach 1.20 phiHyd(i,j,k) = 0. _d 0
236     sigmaX(i,j,k) = 0. _d 0
237     sigmaY(i,j,k) = 0. _d 0
238     sigmaR(i,j,k) = 0. _d 0
239 adcroft 1.1 ConvectCount(i,j,k) = 0.
240     KappaRT(i,j,k) = 0. _d 0
241     KappaRS(i,j,k) = 0. _d 0
242 heimbach 1.5 #ifdef ALLOW_AUTODIFF_TAMC
243     gT(i,j,k,bi,bj) = 0. _d 0
244     gS(i,j,k,bi,bj) = 0. _d 0
245     #ifdef ALLOW_PASSIVE_TRACER
246     gTr1(i,j,k,bi,bj) = 0. _d 0
247     #endif
248 heimbach 1.20 #ifdef ALLOW_GMREDI
249     Kwx(i,j,k,bi,bj) = 0. _d 0
250     Kwy(i,j,k,bi,bj) = 0. _d 0
251     Kwz(i,j,k,bi,bj) = 0. _d 0
252     #ifdef GM_NON_UNITY_DIAGONAL
253     Kux(i,j,k,bi,bj) = 0. _d 0
254     Kvy(i,j,k,bi,bj) = 0. _d 0
255     #endif
256     #endif /* ALLOW_GMREDI */
257 heimbach 1.5 #endif
258 adcroft 1.1 ENDDO
259     ENDDO
260     ENDDO
261    
262 jmc 1.21 iMin = 1-OLx
263 adcroft 1.1 iMax = sNx+OLx
264 jmc 1.21 jMin = 1-OLy
265 adcroft 1.1 jMax = sNy+OLy
266    
267    
268     #ifdef ALLOW_AUTODIFF_TAMC
269 heimbach 1.11 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
270     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
271     #ifdef ALLOW_KPP
272     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
273     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
274     #endif
275 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
276    
277     C-- Start of diagnostic loop
278     DO k=Nr,1,-1
279    
280     #ifdef ALLOW_AUTODIFF_TAMC
281     C? Patrick, is this formula correct now that we change the loop range?
282     C? Do we still need this?
283     cph kkey formula corrected.
284     cph Needed for rhok, rhokm1, in the case useGMREDI.
285     kkey = (ikey-1)*Nr + k
286     CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
287     CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
288     #endif /* ALLOW_AUTODIFF_TAMC */
289    
290     C-- Integrate continuity vertically for vertical velocity
291     CALL INTEGRATE_FOR_W(
292     I bi, bj, k, uVel, vVel,
293     O wVel,
294     I myThid )
295    
296     #ifdef ALLOW_OBCS
297     #ifdef ALLOW_NONHYDROSTATIC
298     C-- Apply OBC to W if in N-H mode
299     IF (useOBCS.AND.nonHydrostatic) THEN
300     CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
301     ENDIF
302     #endif /* ALLOW_NONHYDROSTATIC */
303     #endif /* ALLOW_OBCS */
304    
305 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
306     C-- MOST of THERMODYNAMICS will be disabled
307     #ifndef SINGLE_LAYER_MODE
308    
309 adcroft 1.1 C-- Calculate gradients of potential density for isoneutral
310     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
311     c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
312     IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
313     #ifdef ALLOW_AUTODIFF_TAMC
314     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
315     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
316     #endif /* ALLOW_AUTODIFF_TAMC */
317     CALL FIND_RHO(
318     I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
319     I theta, salt,
320     O rhoK,
321     I myThid )
322     IF (k.GT.1) THEN
323     #ifdef ALLOW_AUTODIFF_TAMC
324     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
325     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
326     #endif /* ALLOW_AUTODIFF_TAMC */
327     CALL FIND_RHO(
328     I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
329     I theta, salt,
330     O rhoKm1,
331     I myThid )
332     ENDIF
333     CALL GRAD_SIGMA(
334     I bi, bj, iMin, iMax, jMin, jMax, k,
335     I rhoK, rhoKm1, rhoK,
336     O sigmaX, sigmaY, sigmaR,
337     I myThid )
338     ENDIF
339    
340     C-- Implicit Vertical Diffusion for Convection
341     c ==> should use sigmaR !!!
342     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
343     CALL CALC_IVDC(
344     I bi, bj, iMin, iMax, jMin, jMax, k,
345     I rhoKm1, rhoK,
346     U ConvectCount, KappaRT, KappaRS,
347     I myTime, myIter, myThid)
348     ENDIF
349    
350 heimbach 1.22 #endif /* SINGLE_LAYER_MODE */
351    
352 adcroft 1.1 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 heimbach 1.11 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
358 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
359    
360     #ifdef ALLOW_OBCS
361     C-- Calculate future values on open boundaries
362     IF (useOBCS) THEN
363 jmc 1.16 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
364 adcroft 1.1 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 heimbach 1.22 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
387     C-- MOST of THERMODYNAMICS will be disabled
388     #ifndef SINGLE_LAYER_MODE
389    
390 adcroft 1.1 #ifdef ALLOW_GMREDI
391    
392 heimbach 1.22 #ifdef ALLOW_AUTODIFF_TAMC
393     CADJ STORE sigmaX(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
394     CADJ STORE sigmaY(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
395     CADJ STORE sigmaR(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
396     #endif /* ALLOW_AUTODIFF_TAMC */
397 adcroft 1.1 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
398     IF (useGMRedi) THEN
399 jmc 1.15 CALL GMREDI_CALC_TENSOR(
400     I bi, bj, iMin, iMax, jMin, jMax,
401 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
402     I myThid )
403     #ifdef ALLOW_AUTODIFF_TAMC
404     ELSE
405 jmc 1.15 CALL GMREDI_CALC_TENSOR_DUMMY(
406     I bi, bj, iMin, iMax, jMin, jMax,
407 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
408     I myThid )
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 & , KPPdiffKzT(:,:,:,bi,bj)
435     CADJ & , KPPdiffKzS(:,:,:,bi,bj)
436     CADJ & , KPPfrac (:,: ,bi,bj)
437     CADJ & = comlev1_bibj, key=ikey, byte=isbyte
438     #endif /* ALLOW_AUTODIFF_TAMC */
439    
440     #endif /* ALLOW_KPP */
441    
442     #ifdef ALLOW_AUTODIFF_TAMC
443 heimbach 1.11 CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
444     CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
445     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
446     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
447     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
448     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
449 adcroft 1.1 #ifdef ALLOW_PASSIVE_TRACER
450 heimbach 1.11 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
451 adcroft 1.1 #endif
452     #endif /* ALLOW_AUTODIFF_TAMC */
453    
454     #ifdef ALLOW_AIM
455     C AIM - atmospheric intermediate model, physics package code.
456     C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
457     IF ( useAIM ) THEN
458     CALL TIMER_START('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
459     CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )
460     CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid)
461     ENDIF
462     #endif /* ALLOW_AIM */
463 jmc 1.14
464     #ifdef ALLOW_TIMEAVE
465     IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
466     CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
467     I deltaTclock, bi, bj, myThid)
468     ENDIF
469     #endif /* ALLOW_TIMEAVE */
470 adcroft 1.4
471 adcroft 1.12 #ifndef DISABLE_MULTIDIM_ADVECTION
472 adcroft 1.4 C-- Some advection schemes are better calculated using a multi-dimensional
473     C method in the absence of any other terms and, if used, is done here.
474 adcroft 1.13 C
475     C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
476     C The default is to use multi-dimensinal advection for non-linear advection
477     C schemes. However, for the sake of efficiency of the adjoint it is necessary
478     C to be able to exclude this scheme to avoid excessive storage and
479     C recomputation. It *is* differentiable, if you need it.
480     C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
481     C disable this section of code.
482 adcroft 1.6 IF (multiDimAdvection) THEN
483     IF (tempStepping .AND.
484     & tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.
485     & tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.
486 heimbach 1.11 & tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN
487     CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
488     U theta,gT,
489 adcroft 1.6 I myTime,myIter,myThid)
490 heimbach 1.11 ENDIF
491 adcroft 1.6 IF (saltStepping .AND.
492     & saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.
493     & saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.
494 heimbach 1.11 & saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN
495     CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
496     U salt,gS,
497 adcroft 1.6 I myTime,myIter,myThid)
498 heimbach 1.11 ENDIF
499 adcroft 1.6 ENDIF
500 adcroft 1.17 C Since passive tracers are configurable separately from T,S we
501     C call the multi-dimensional method for PTRACERS regardless
502     C of whether multiDimAdvection is set or not.
503     #ifdef ALLOW_PTRACERS
504     IF ( usePTRACERS ) THEN
505     CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
506     ENDIF
507     #endif /* ALLOW_PTRACERS */
508 adcroft 1.12 #endif /* DISABLE_MULTIDIM_ADVECTION */
509 adcroft 1.1
510     C-- Start of thermodynamics loop
511     DO k=Nr,1,-1
512     #ifdef ALLOW_AUTODIFF_TAMC
513     C? Patrick Is this formula correct?
514     cph Yes, but I rewrote it.
515     cph Also, the KappaR? need the index and subscript k!
516     kkey = (ikey-1)*Nr + k
517     #endif /* ALLOW_AUTODIFF_TAMC */
518    
519     C-- km1 Points to level above k (=k-1)
520     C-- kup Cycles through 1,2 to point to layer above
521     C-- kDown Cycles through 2,1 to point to current layer
522    
523     km1 = MAX(1,k-1)
524     kup = 1+MOD(k+1,2)
525     kDown= 1+MOD(k,2)
526    
527     iMin = 1-OLx
528     iMax = sNx+OLx
529     jMin = 1-OLy
530     jMax = sNy+OLy
531    
532     C-- Get temporary terms used by tendency routines
533     CALL CALC_COMMON_FACTORS (
534     I bi,bj,iMin,iMax,jMin,jMax,k,
535     O xA,yA,uTrans,vTrans,rTrans,maskUp,
536     I myThid)
537 jmc 1.19
538     #ifdef ALLOW_GMREDI
539     C-- Residual transp = Bolus transp + Eulerian transp
540     IF (useGMRedi) THEN
541     CALL GMREDI_CALC_UVFLOW(
542     & uTrans, vTrans, bi, bj, k, myThid)
543     IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
544     & rTrans, bi, bj, k, myThid)
545     ENDIF
546     #endif /* ALLOW_GMREDI */
547 adcroft 1.1
548     #ifdef ALLOW_AUTODIFF_TAMC
549     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
550     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
551     #endif /* ALLOW_AUTODIFF_TAMC */
552    
553     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
554     C-- Calculate the total vertical diffusivity
555     CALL CALC_DIFFUSIVITY(
556     I bi,bj,iMin,iMax,jMin,jMax,k,
557     I maskUp,
558 heimbach 1.2 O KappaRT,KappaRS,
559 adcroft 1.1 I myThid)
560     #endif
561    
562     iMin = 1-OLx+2
563     iMax = sNx+OLx-1
564     jMin = 1-OLy+2
565     jMax = sNy+OLy-1
566    
567     C-- Calculate active tracer tendencies (gT,gS,...)
568     C and step forward storing result in gTnm1, gSnm1, etc.
569     IF ( tempStepping ) THEN
570     CALL CALC_GT(
571     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
572     I xA,yA,uTrans,vTrans,rTrans,maskUp,
573     I KappaRT,
574     U fVerT,
575 adcroft 1.7 I myTime,myIter,myThid)
576 adcroft 1.1 CALL TIMESTEP_TRACER(
577 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
578 adcroft 1.1 I theta, gT,
579     I myIter, myThid)
580     ENDIF
581     IF ( saltStepping ) THEN
582     CALL CALC_GS(
583     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
584     I xA,yA,uTrans,vTrans,rTrans,maskUp,
585     I KappaRS,
586     U fVerS,
587 adcroft 1.7 I myTime,myIter,myThid)
588 adcroft 1.1 CALL TIMESTEP_TRACER(
589 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
590 adcroft 1.1 I salt, gS,
591     I myIter, myThid)
592     ENDIF
593     #ifdef ALLOW_PASSIVE_TRACER
594     IF ( tr1Stepping ) THEN
595     CALL CALC_GTR1(
596     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
597     I xA,yA,uTrans,vTrans,rTrans,maskUp,
598     I KappaRT,
599     U fVerTr1,
600 heimbach 1.8 I myTime,myIter,myThid)
601 adcroft 1.1 CALL TIMESTEP_TRACER(
602 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
603 adcroft 1.1 I Tr1, gTr1,
604 heimbach 1.8 I myIter,myThid)
605 adcroft 1.1 ENDIF
606     #endif
607 adcroft 1.17 #ifdef ALLOW_PTRACERS
608     IF ( usePTRACERS ) THEN
609     CALL PTRACERS_INTEGERATE(
610     I bi,bj,k,
611     I xA,yA,uTrans,vTrans,rTrans,maskUp,
612     X KappaRS,
613     I myIter,myTime,myThid)
614     ENDIF
615     #endif /* ALLOW_PTRACERS */
616 adcroft 1.1
617     #ifdef ALLOW_OBCS
618     C-- Apply open boundary conditions
619     IF (useOBCS) THEN
620 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
621 adcroft 1.1 END IF
622     #endif /* ALLOW_OBCS */
623    
624     C-- Freeze water
625     IF (allowFreezing) THEN
626     #ifdef ALLOW_AUTODIFF_TAMC
627 heimbach 1.11 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
628 adcroft 1.1 CADJ & , key = kkey, byte = isbyte
629     #endif /* ALLOW_AUTODIFF_TAMC */
630     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
631     END IF
632    
633     C-- end of thermodynamic k loop (Nr:1)
634     ENDDO
635    
636    
637     #ifdef ALLOW_AUTODIFF_TAMC
638     C? Patrick? What about this one?
639 adcroft 1.10 cph Keys iikey and idkey dont seem to be needed
640 adcroft 1.1 cph since storing occurs on different tape for each
641     cph impldiff call anyways.
642 adcroft 1.10 cph Thus, common block comlev1_impl isnt needed either.
643 adcroft 1.1 cph Storing below needed in the case useGMREDI.
644     iikey = (ikey-1)*maximpl
645     #endif /* ALLOW_AUTODIFF_TAMC */
646    
647     C-- Implicit diffusion
648     IF (implicitDiffusion) THEN
649    
650     IF (tempStepping) THEN
651     #ifdef ALLOW_AUTODIFF_TAMC
652     idkey = iikey + 1
653 heimbach 1.11 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
654 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
655     CALL IMPLDIFF(
656     I bi, bj, iMin, iMax, jMin, jMax,
657     I deltaTtracer, KappaRT, recip_HFacC,
658 adcroft 1.7 U gT,
659 adcroft 1.1 I myThid )
660     ENDIF
661    
662     IF (saltStepping) THEN
663     #ifdef ALLOW_AUTODIFF_TAMC
664     idkey = iikey + 2
665 heimbach 1.11 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
666 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
667     CALL IMPLDIFF(
668     I bi, bj, iMin, iMax, jMin, jMax,
669     I deltaTtracer, KappaRS, recip_HFacC,
670 adcroft 1.7 U gS,
671 adcroft 1.1 I myThid )
672     ENDIF
673    
674     #ifdef ALLOW_PASSIVE_TRACER
675     IF (tr1Stepping) THEN
676     #ifdef ALLOW_AUTODIFF_TAMC
677 heimbach 1.11 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
678 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
679     CALL IMPLDIFF(
680     I bi, bj, iMin, iMax, jMin, jMax,
681     I deltaTtracer, KappaRT, recip_HFacC,
682 adcroft 1.7 U gTr1,
683 adcroft 1.1 I myThid )
684     ENDIF
685     #endif
686    
687 adcroft 1.17 #ifdef ALLOW_PTRACERS
688     C Vertical diffusion (implicit) for passive tracers
689     IF ( usePTRACERS ) THEN
690     CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
691     ENDIF
692     #endif /* ALLOW_PTRACERS */
693    
694 adcroft 1.1 #ifdef ALLOW_OBCS
695     C-- Apply open boundary conditions
696     IF (useOBCS) THEN
697     DO K=1,Nr
698 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
699 adcroft 1.1 ENDDO
700     END IF
701     #endif /* ALLOW_OBCS */
702    
703     C-- End If implicitDiffusion
704     ENDIF
705 heimbach 1.22
706     #endif /* SINGLE_LAYER_MODE */
707 adcroft 1.1
708     Ccs-
709     ENDDO
710     ENDDO
711    
712     #ifdef ALLOW_AIM
713     IF ( useAIM ) THEN
714     CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
715     ENDIF
716 adcroft 1.7 _EXCH_XYZ_R8(gT,myThid)
717     _EXCH_XYZ_R8(gS,myThid)
718 adcroft 1.1 #else
719     IF (staggerTimeStep.AND.useCubedSphereExchange) THEN
720 adcroft 1.7 _EXCH_XYZ_R8(gT,myThid)
721     _EXCH_XYZ_R8(gS,myThid)
722 adcroft 1.1 ENDIF
723     #endif /* ALLOW_AIM */
724 adcroft 1.17
725     #ifndef DISABLE_DEBUGMODE
726     If (debugMode) THEN
727     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
728     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
729     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
730     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
731     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
732     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
733     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
734     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
735     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
736 adcroft 1.18 #ifdef ALLOW_PTRACERS
737     IF ( usePTRACERS ) THEN
738     CALL PTRACERS_DEBUG(myThid)
739     ENDIF
740     #endif /* ALLOW_PTRACERS */
741 adcroft 1.17 ENDIF
742     #endif
743 adcroft 1.1
744     RETURN
745     END

  ViewVC Help
Powered by ViewVC 1.1.22