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