/[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.24 - (hide annotations) (download)
Wed Jul 3 20:22:39 2002 UTC (21 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint45d_post
Changes since 1.23: +13 -10 lines 
Oly must be > 3 when using an advection scheme other than 2nd O  and
 staggerTimeStep & multiDimAdvection & GM_AdvForm & not.GM_AdvSeparate
Add an exchange of gT,gS in this case.
1 jmc 1.24 C $Header: /u/gcmpack/MITgcm/model/src/thermodynamics.F,v 1.23 2002/06/15 03:28:39 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 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    
629     #ifdef ALLOW_AUTODIFF_TAMC
630     C? Patrick? What about this one?
631 adcroft 1.10 cph Keys iikey and idkey dont seem to be needed
632 adcroft 1.1 cph since storing occurs on different tape for each
633     cph impldiff call anyways.
634 adcroft 1.10 cph Thus, common block comlev1_impl isnt needed either.
635 adcroft 1.1 cph Storing below needed in the case useGMREDI.
636     iikey = (ikey-1)*maximpl
637     #endif /* ALLOW_AUTODIFF_TAMC */
638    
639     C-- Implicit diffusion
640     IF (implicitDiffusion) THEN
641    
642     IF (tempStepping) THEN
643     #ifdef ALLOW_AUTODIFF_TAMC
644     idkey = iikey + 1
645 heimbach 1.11 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
646 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
647     CALL IMPLDIFF(
648     I bi, bj, iMin, iMax, jMin, jMax,
649     I deltaTtracer, KappaRT, recip_HFacC,
650 adcroft 1.7 U gT,
651 adcroft 1.1 I myThid )
652     ENDIF
653    
654     IF (saltStepping) THEN
655     #ifdef ALLOW_AUTODIFF_TAMC
656     idkey = iikey + 2
657 heimbach 1.11 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
658 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
659     CALL IMPLDIFF(
660     I bi, bj, iMin, iMax, jMin, jMax,
661     I deltaTtracer, KappaRS, recip_HFacC,
662 adcroft 1.7 U gS,
663 adcroft 1.1 I myThid )
664     ENDIF
665    
666     #ifdef ALLOW_PASSIVE_TRACER
667     IF (tr1Stepping) THEN
668     #ifdef ALLOW_AUTODIFF_TAMC
669 heimbach 1.11 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
670 adcroft 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
671     CALL IMPLDIFF(
672     I bi, bj, iMin, iMax, jMin, jMax,
673     I deltaTtracer, KappaRT, recip_HFacC,
674 adcroft 1.7 U gTr1,
675 adcroft 1.1 I myThid )
676     ENDIF
677     #endif
678    
679 adcroft 1.17 #ifdef ALLOW_PTRACERS
680     C Vertical diffusion (implicit) for passive tracers
681     IF ( usePTRACERS ) THEN
682     CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
683     ENDIF
684     #endif /* ALLOW_PTRACERS */
685    
686 adcroft 1.1 #ifdef ALLOW_OBCS
687     C-- Apply open boundary conditions
688     IF (useOBCS) THEN
689     DO K=1,Nr
690 adcroft 1.7 CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
691 adcroft 1.1 ENDDO
692     END IF
693     #endif /* ALLOW_OBCS */
694    
695     C-- End If implicitDiffusion
696     ENDIF
697 heimbach 1.22
698     #endif /* SINGLE_LAYER_MODE */
699 adcroft 1.1
700     Ccs-
701     ENDDO
702     ENDDO
703    
704     #ifdef ALLOW_AIM
705     IF ( useAIM ) THEN
706     CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
707     ENDIF
708     #endif /* ALLOW_AIM */
709 jmc 1.24 IF ( staggerTimeStep ) THEN
710     IF ( useAIM .OR. useCubedSphereExchange ) THEN
711     IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
712     IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
713     ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
714     c .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
715     IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
716     IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
717     ENDIF
718     ENDIF
719 adcroft 1.17
720     #ifndef DISABLE_DEBUGMODE
721     If (debugMode) THEN
722     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
723     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
724     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
725     CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
726     CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
727     CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
728     CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
729     CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
730     CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
731 adcroft 1.18 #ifdef ALLOW_PTRACERS
732     IF ( usePTRACERS ) THEN
733     CALL PTRACERS_DEBUG(myThid)
734     ENDIF
735     #endif /* ALLOW_PTRACERS */
736 adcroft 1.17 ENDIF
737     #endif
738 adcroft 1.1
739     RETURN
740     END
  ViewVC Help
Powered by ViewVC 1.1.22