/[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.20 - (hide annotations) (download)
Sun Mar 24 02:36:39 2002 UTC (22 years, 3 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint44h_post, checkpoint45
Changes since 1.19: +21 -13 lines
o Modified initialisations to break adjoint dependencies
o removed some store directives
o added options files for KPP, GMREDI

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

  ViewVC Help
Powered by ViewVC 1.1.22