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