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