/[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.28 - (hide annotations) (download)
Wed Oct 9 00:57:57 2002 UTC (21 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint46l_post, checkpoint46l_pre, checkpoint46k_post, checkpoint46m_post
Changes since 1.27: +2 -2 lines 
* S/R aim_initialise.F replace S/R aim_init from file aim_do_inphys.F:
 - read AIM physics parameters from a file (data.aimphys)
 - set defaults values = F.Molteni paper (Clim.Dyn., 2002)
 - phiHyd removed from arg. list of S/R aim_do_atmos_physics
1 jmc 1.28 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.27 2002/10/07 16:24:45 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 jmc 1.27 c CALL INTEGRATE_FOR_W(
292     c I bi, bj, k, uVel, vVel,
293     c O wVel,
294     c I myThid )
295 adcroft 1.1
296     #ifdef ALLOW_OBCS
297     #ifdef ALLOW_NONHYDROSTATIC
298     C-- Apply OBC to W if in N-H mode
299 jmc 1.27 c IF (useOBCS.AND.nonHydrostatic) THEN
300     c CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
301     c ENDIF
302 adcroft 1.1 #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 mlosch 1.26 I bi, bj, iMin, iMax, jMin, jMax, k, k,
319 adcroft 1.1 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 mlosch 1.26 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
329 adcroft 1.1 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 jmc 1.28 CALL AIM_DO_ATMOS_PHYSICS( bi, bj, myTime, myThid )
460 adcroft 1.1 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 jmc 1.24 IF (tempMultiDimAdvec) THEN
483 heimbach 1.11 CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
484     U theta,gT,
485 adcroft 1.6 I myTime,myIter,myThid)
486 jmc 1.23 ENDIF
487 jmc 1.24 IF (saltMultiDimAdvec) THEN
488 heimbach 1.11 CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
489     U salt,gS,
490 adcroft 1.6 I myTime,myIter,myThid)
491     ENDIF
492 adcroft 1.17 C Since passive tracers are configurable separately from T,S we
493     C call the multi-dimensional method for PTRACERS regardless
494     C of whether multiDimAdvection is set or not.
495     #ifdef ALLOW_PTRACERS
496     IF ( usePTRACERS ) THEN
497     CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
498     ENDIF
499     #endif /* ALLOW_PTRACERS */
500 adcroft 1.12 #endif /* DISABLE_MULTIDIM_ADVECTION */
501 adcroft 1.1
502     C-- Start of thermodynamics loop
503     DO k=Nr,1,-1
504     #ifdef ALLOW_AUTODIFF_TAMC
505     C? Patrick Is this formula correct?
506     cph Yes, but I rewrote it.
507     cph Also, the KappaR? need the index and subscript k!
508     kkey = (ikey-1)*Nr + k
509     #endif /* ALLOW_AUTODIFF_TAMC */
510    
511     C-- km1 Points to level above k (=k-1)
512     C-- kup Cycles through 1,2 to point to layer above
513     C-- kDown Cycles through 2,1 to point to current layer
514    
515     km1 = MAX(1,k-1)
516     kup = 1+MOD(k+1,2)
517     kDown= 1+MOD(k,2)
518    
519     iMin = 1-OLx
520     iMax = sNx+OLx
521     jMin = 1-OLy
522     jMax = sNy+OLy
523    
524     C-- Get temporary terms used by tendency routines
525     CALL CALC_COMMON_FACTORS (
526     I bi,bj,iMin,iMax,jMin,jMax,k,
527     O xA,yA,uTrans,vTrans,rTrans,maskUp,
528     I myThid)
529 jmc 1.19
530     #ifdef ALLOW_GMREDI
531     C-- Residual transp = Bolus transp + Eulerian transp
532     IF (useGMRedi) THEN
533     CALL GMREDI_CALC_UVFLOW(
534     & uTrans, vTrans, bi, bj, k, myThid)
535     IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
536     & rTrans, bi, bj, k, myThid)
537     ENDIF
538     #endif /* ALLOW_GMREDI */
539 adcroft 1.1
540     #ifdef ALLOW_AUTODIFF_TAMC
541     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
542     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
543     #endif /* ALLOW_AUTODIFF_TAMC */
544    
545     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
546     C-- Calculate the total vertical diffusivity
547     CALL CALC_DIFFUSIVITY(
548     I bi,bj,iMin,iMax,jMin,jMax,k,
549     I maskUp,
550 heimbach 1.2 O KappaRT,KappaRS,
551 adcroft 1.1 I myThid)
552     #endif
553    
554     iMin = 1-OLx+2
555     iMax = sNx+OLx-1
556     jMin = 1-OLy+2
557     jMax = sNy+OLy-1
558    
559     C-- Calculate active tracer tendencies (gT,gS,...)
560     C and step forward storing result in gTnm1, gSnm1, etc.
561     IF ( tempStepping ) THEN
562     CALL CALC_GT(
563     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
564     I xA,yA,uTrans,vTrans,rTrans,maskUp,
565     I KappaRT,
566     U fVerT,
567 adcroft 1.7 I myTime,myIter,myThid)
568 adcroft 1.1 CALL TIMESTEP_TRACER(
569 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
570 adcroft 1.1 I theta, gT,
571     I myIter, myThid)
572     ENDIF
573     IF ( saltStepping ) THEN
574     CALL CALC_GS(
575     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
576     I xA,yA,uTrans,vTrans,rTrans,maskUp,
577     I KappaRS,
578     U fVerS,
579 adcroft 1.7 I myTime,myIter,myThid)
580 adcroft 1.1 CALL TIMESTEP_TRACER(
581 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
582 adcroft 1.1 I salt, gS,
583     I myIter, myThid)
584     ENDIF
585     #ifdef ALLOW_PASSIVE_TRACER
586     IF ( tr1Stepping ) THEN
587     CALL CALC_GTR1(
588     I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
589     I xA,yA,uTrans,vTrans,rTrans,maskUp,
590     I KappaRT,
591     U fVerTr1,
592 heimbach 1.8 I myTime,myIter,myThid)
593 adcroft 1.1 CALL TIMESTEP_TRACER(
594 adcroft 1.3 I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
595 adcroft 1.1 I Tr1, gTr1,
596 heimbach 1.8 I myIter,myThid)
597 adcroft 1.1 ENDIF
598     #endif
599 adcroft 1.17 #ifdef ALLOW_PTRACERS
600     IF ( usePTRACERS ) THEN
601     CALL PTRACERS_INTEGERATE(
602     I bi,bj,k,
603     I xA,yA,uTrans,vTrans,rTrans,maskUp,
604     X KappaRS,
605     I myIter,myTime,myThid)
606     ENDIF
607     #endif /* ALLOW_PTRACERS */
608 adcroft 1.1
609     #ifdef ALLOW_OBCS
610     C-- Apply open boundary conditions
611     IF (useOBCS) THEN
612 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
613 adcroft 1.1 END IF
614     #endif /* ALLOW_OBCS */
615    
616     C-- Freeze water
617     IF (allowFreezing) THEN
618     #ifdef ALLOW_AUTODIFF_TAMC
619 heimbach 1.11 CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
620 adcroft 1.1 CADJ & , key = kkey, byte = isbyte
621     #endif /* ALLOW_AUTODIFF_TAMC */
622     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
623     END IF
624    
625     C-- end of thermodynamic k loop (Nr:1)
626     ENDDO
627    
628     C-- Implicit diffusion
629     IF (implicitDiffusion) THEN
630    
631     IF (tempStepping) THEN
632     #ifdef ALLOW_AUTODIFF_TAMC
633 heimbach 1.11 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
634 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
635     CALL IMPLDIFF(
636     I bi, bj, iMin, iMax, jMin, jMax,
637     I deltaTtracer, KappaRT, recip_HFacC,
638 adcroft 1.7 U gT,
639 adcroft 1.1 I myThid )
640     ENDIF
641    
642     IF (saltStepping) THEN
643     #ifdef ALLOW_AUTODIFF_TAMC
644 heimbach 1.11 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
645 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
646     CALL IMPLDIFF(
647     I bi, bj, iMin, iMax, jMin, jMax,
648     I deltaTtracer, KappaRS, recip_HFacC,
649 adcroft 1.7 U gS,
650 adcroft 1.1 I myThid )
651     ENDIF
652    
653     #ifdef ALLOW_PASSIVE_TRACER
654     IF (tr1Stepping) THEN
655     #ifdef ALLOW_AUTODIFF_TAMC
656 heimbach 1.11 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
657 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
658     CALL IMPLDIFF(
659     I bi, bj, iMin, iMax, jMin, jMax,
660     I deltaTtracer, KappaRT, recip_HFacC,
661 adcroft 1.7 U gTr1,
662 adcroft 1.1 I myThid )
663     ENDIF
664     #endif
665    
666 adcroft 1.17 #ifdef ALLOW_PTRACERS
667     C Vertical diffusion (implicit) for passive tracers
668     IF ( usePTRACERS ) THEN
669     CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
670     ENDIF
671     #endif /* ALLOW_PTRACERS */
672    
673 adcroft 1.1 #ifdef ALLOW_OBCS
674     C-- Apply open boundary conditions
675     IF (useOBCS) THEN
676     DO K=1,Nr
677 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
678 adcroft 1.1 ENDDO
679     END IF
680     #endif /* ALLOW_OBCS */
681    
682     C-- End If implicitDiffusion
683     ENDIF
684 heimbach 1.22
685     #endif /* SINGLE_LAYER_MODE */
686 adcroft 1.1
687     Ccs-
688     ENDDO
689     ENDDO
690    
691     #ifdef ALLOW_AIM
692 jmc 1.27 c IF ( useAIM ) THEN
693     c CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
694     c ENDIF
695 adcroft 1.1 #endif /* ALLOW_AIM */
696 jmc 1.27 c IF ( staggerTimeStep ) THEN
697     c IF ( useAIM .OR. useCubedSphereExchange ) THEN
698     c IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
699     c IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
700     c ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
701     cc .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
702     c IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
703     c IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
704     c ENDIF
705     c ENDIF
706 adcroft 1.17
707     #ifndef DISABLE_DEBUGMODE
708     If (debugMode) THEN
709     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
710     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
711     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
712     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
713     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
714     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
715     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
716     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
717     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
718 adcroft 1.18 #ifdef ALLOW_PTRACERS
719     IF ( usePTRACERS ) THEN
720     CALL PTRACERS_DEBUG(myThid)
721     ENDIF
722     #endif /* ALLOW_PTRACERS */
723 adcroft 1.17 ENDIF
724     #endif
725 adcroft 1.1
726     RETURN
727     END
  ViewVC Help
Powered by ViewVC 1.1.22