/[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.31 - (hide annotations) (download)
Fri Nov 15 19:58:21 2002 UTC (21 years, 7 months ago) by cheisey
Branch: MAIN
Changes since 1.30: +38 -2 lines 
Two packages:  bulk_forcing (Bulk forcing)
and thermodynamic_seaice (adapted from LANL CICE.v2.0.2)
Earlier integration from Stephaine Dutkiewicz
and Patrick Heimbach.

Two ifdef statements for compile time,
ALLOW_TSEAICE and ALLOW_BULKFORMULA

Two switches in data.pkg to turn on at run-time:

cat data.pkg
# Packages
 &PACKAGES
 useBulkf=.TRUE.,
 useTSeaIce=.TRUE.,
 &

The bulk package requires an additional parameter file
with two namelists.

 cat data.blk
 &BULKF_PARM01
 RainFile=       'ncep_precip_m_cubed.bin',
 SolarFile=      'ncep_downsolr_cubed.bin',
 AirTempFile=    'ncep_tair_cubed.bin',
 AirhumidityFile='ncep_qair_g_cubed.bin',
 LongwaveFile=   'ncep_netlw_cubed.bin',
 UWindFile=      'ncep_uwind_cubed.bin',
 VWindFile=      'ncep_vwind_cubed.bin',
 WspeedFile=    ' ',
 RunoffFile=    ' ',
 QnetFile=       ' ',
 EmPFile=        'ncep_emp_calc_cubed.bin',
 CloudFile=      'ncep_totalcloud_cubed.bin',
 &

 &BULKF_PARM02
 qnet_off=0.0,
 empmr_off=0.0,
 conservcycle=311040000.,
 &



c ADAPTED FROM:
c LANL CICE.v2.0.2
c-----------------------------------------------------------------------
c.. thermodynamics (vertical physics) based on M. Winton 3-layer model
c.. See Bitz, C. M. and W. H. Lipscomb, 1999:  "An energy-conserving
c..       thermodynamic sea ice model for climate study."  J. Geophys.
c..       Res., 104, 15669 - 15677.
c..     Winton, M., 1999:  "A reformulated three-layer sea ice model."
c..       Submitted to J. Atmos. Ocean. Technol.

c.. authors Elizabeth C. Hunke and William Lipscomb
c..         Fluid Dynamics Group, Los Alamos National Laboratory
c-----------------------------------------------------------------------

WARNING:  useSEAICE and useTSEAICE are mutually exclusive.

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