source: NEMO/trunk/doc/si3_doc/Figures/scripts/Thermal_properties.pro @ 9996

Last change on this file since 9996 was 9982, checked in by vancop, 3 years ago

SI3 doc update

File size: 7.6 KB
Line 
1;
2;------------------------------------------------------------------------------
3;
4; Plot LIM thermal properties
5;
6;------------------------------------------------------------------------------
7
8; Model parameters
9; Sample experiments
10N_spl = 3 & Ts_spl = FLTARR(N_spl)
11Ts_spl[*] = [ -15., -10., -5. ] & Tw_spl = 2.
12colors = [ 50., 100., 150. ]
13
14;
15;==============================================================================
16; 1) User parameters
17;==============================================================================
18;
19
20out_dir=''
21file_out   = 'Thermal_properties'
22
23numplot_x = 2       ; number of horizontal plots
24numplot_y = 2       ; number of vertical plots
25ct = 13             ; colortable
26cs = 1.0            ; size of fonts on plots (1=normal)
27th = 2.             ; thickness of curves
28device = 'PS'       ; 'PS' or 'X'
29
30;
31;==============================================================================
32; 2) Initialize graphics
33;==============================================================================
34;
35
36dx = 6
37dy = 6.
38figuresize_x = numplot_x * dx; figure size on x direction
39figuresize_y = numplot_y * dy; figure size on y direction
40SET_PLOT, DEVICE
41;!P.FONT=1
42TVLCT, [0,255,0,0], [0,0,255,0], [0,0,0,255]
43IF ( device EQ 'PS' ) THEN BEGIN
44;   DEVICE, /COLOR, /LANDSCAPE, filename=out_dir+file_out+'.ps', $
45;           XSIZE=figuresize_x,YSIZE=figuresize_y,FONT_SIZE=12.0, SET_FONT='Helvetica', BITS = 16
46   DEVICE, /COLOR, /LANDSCAPE, filename=out_dir+file_out+'.ps', $
47           XSIZE=figuresize_x,YSIZE=figuresize_y,FONT_SIZE=9.0
48ENDIF
49
50IF ( device EQ 'X' ) THEN BEGIN
51   xsize = 1200
52   ysize = 800
53   colorkey = 'rd'
54   init_graphics_x, xsize, ysize, colorkey
55ENDIF
56
57!P.MULTI=[0,numplot_x,numplot_y]
58!X.MARGIN   = [6,3]
59!Y.MARGIN   = [4,3]
60
61;==============================================================================
62; 3) Computations and plots
63;==============================================================================
64
65; arrays
66N_T = 200
67N_S = 3
68miss_val = 9.99e32
69
70;--- Physical constants
71c0 = 2067.d0
72cw = 3991.d0
73L  = 334000.d0
74mu = 0.054d0
75rhoi = 917.d0
76rhow = 1026.d0
77rhos = 330.d0
78
79Sw = 34.d0
80
81beta1 = 0.09; pringle coefficients
82beta2 = 0.011;
83
84;--- Ice salinity
85;S = ( rhoi - rhos ) / rhoi * Sw
86S = FLTARR(N_S)
87Smin = 5.   
88Smax = 15.
89dS = ( Smax - Smin ) / FLOAT(N_S-1)
90S[0:N_S-1] = Smin + FINDGEN(N_S) * dS
91
92Tf = -mu * S
93
94;--- Temperatures
95T = FLTARR(N_T)
96Tmin = -20.
97Tf   = -mu * Smax
98Tmax = 0.
99dT = ( Tmax - Tmin ) / FLOAT(N_T-1)
100T[0:N_T-1] = Tmin + FINDGEN(N_T) * dT
101
102;--- Grids
103S_grid = FLTARR(N_S,N_T)
104T_grid = FLTARR(N_S,N_T)
105
106FOR j = 0, N_T - 1 DO BEGIN
107   S_grid(0:N_S-1,j) = S(0:N_S-1)
108ENDFOR
109
110FOR i = 0, N_S - 1 DO BEGIN
111   T_grid(i,0:N_T-1) = T(0:N_T-1)
112ENDFOR
113
114;--- Brine fraction
115phi = -mu * S_grid / T_grid
116
117;--- specific enthalpy
118q_i = c0 * ( T_grid + mu * S_grid ) - L * ( 1 + mu * S_grid / T_grid ) - cw*mu*S_grid
119
120;--- thermal conductivity (pringle)
121kpri = 2.11 + beta1 * S_grid / T_grid - beta2 * T_grid;
122
123;--- specific heat
124cp = c0 + L*mu*S_grid/(T_grid*T_grid)
125
126;--- PLOTS
127!P.MULTI=[0, numplot_x, numplot_y]
128LOADCT, 13
129
130PLOT, [Tmin, Tmax], [0., 50.], /NODATA, XTITLE='T (C)', YTITLE='phi (%)', CHARSIZE = cs
131FOR i = 0, N_S - 1 DO BEGIN
132   zaddr = where(phi LT 100)
133   OPLOT, T_grid(i,zaddr), phi(i,zaddr)*100., thick = th, color = i* 50
134ENDFOR
135
136PLOT, [Tmin, Tmax], [-4., 0.], /NODATA, XTITLE='T (C)', YTITLE='q_m (10^5 J/kg)', CHARSIZE = cs
137FOR i = 0, N_S - 1 DO BEGIN
138   zaddr = where(phi LT 100)
139   OPLOT, T_grid(i,zaddr), q_i(i,zaddr)/1.0e5, thick = th, color = i* 50
140ENDFOR
141
142PLOT, [Tmin, Tmax], [0., 3.], /NODATA, XTITLE='T (C)', YTITLE='k_i (W/m/K)', CHARSIZE = cs
143FOR i = 0, N_S - 1 DO BEGIN
144   zaddr = where(phi LT 100)
145   OPLOT, T_grid(i,zaddr), kpri(i,zaddr), thick = th, color = i* 50
146ENDFOR
147
148PLOT, [Tmin, Tmax], [0., 20.], /NODATA, XTITLE='T (C)', YTITLE='c_i (10^3 J/kg/K)', CHARSIZE = cs
149FOR i = 0, N_S - 1 DO BEGIN
150   zaddr = where(phi LT 100)
151   OPLOT, T_grid(i,zaddr), cp(i,zaddr)/1000., thick = th, color = i* 50
152ENDFOR
153
154
155;T[0:N_T-1] = FINDGEN(N_T) * dT + Tmin
156;dT = ( Tmax - Tf) / FLOAT(N/2-1)
157;T[N/2:N-1] = FINDGEN(N/2) * dT + Tf
158
159;;PLOT, T, TITLE = "T", CHARSIZE = cs
160;
161;; Specific heat
162;ci = FLTARR(N) & ci(*) = miss_val
163;cis = FLTARR(N) & cis(*) = miss_val
164;ci(0:N/2) = c0 + L*mu*S/T(0:N/2)/T(0:N/2)
165;cis(0:N/2) = mu*(c0 - L/T(0:N/2))
166;;PLOT, T(0:N/2), ci(0:N/2)/L, CHARSIZE = cs
167;;OPLOT, T(0:N/2), cis(0:N/2)/L, LINESTYLE = 1
168;
169;;--- Brine volume
170;e_i = FLTARR(N) & e_i(*) = 1.
171;i_addr = FINDGEN(N/2+1)
172;e_i(i_addr) = -mu * S / T(i_addr)
173;PLOT, T, e_i, CHARSIZE = cs, TITLE = "e"
174;OPLOT, T, REPLICATE(0.5, N), linestyle = 1
175;T_50 = -mu * S / 0.50 ; temperature at which brine volume is 50%
176;
177;; Notz (2005), Phd thesis, page 44
178;i_addr = WHERE(T LT Tf)
179;Sbr = - 21.4 * T(i_addr) - 0.886 * T(i_addr) * T(i_addr) - 0.0170 * T(i_addr) * T(i_addr) * T(i_addr)
180;zrS = S / Sbr(i_addr)
181;beta_ocs = 0.81
182;
183;zrhoi      = 916.8 - 0.1403 * T(i_addr)              ; pure ice density Pounder(1965)
184;zrhol      = rhow + beta_ocs * ( Sbr(i_addr) - Sw )
185;zrR        = zrhoi / zrhol
186;znum       = 1 - zrS
187;zden       = 1 + zrS * ( zrR - 1 )
188;e_i2 = FLTARR(N) & e_i2(*) = 1.
189;e_i2(i_addr) =  1. - znum / zden
190;;OPLOT, T, e_i2, linestyle = 3
191;
192;;OPLOT, [ T_50, T_50], [ 0., 1.], linestyle = 1
193;
194;;--- Enthalpies
195;Ei = FLTARR(N) & Ei(*) = miss_val
196;Ei2= FLTARR(N) & Ei2(*) = miss_val
197;Es = FLTARR(N) & Es(*) = miss_val
198;Es_wgt = FLTARR(N) & Es(*) = miss_val
199;Ew = FLTARR(N) & Ew(*) = miss_val
200;
201;i_addr = WHERE( T LE 0. )
202;Es(i_addr) = c0 * T(i_addr) - L
203;Es_wgt(i_addr) = Es(i_addr) * rhos / rhoi
204;i_addr = WHERE( T GT 0. )
205;Es(i_addr) = cw * T(i_addr)
206;Es_wgt(i_addr) = Es(i_addr) * ( rhoi - rhos ) / rhoi
207;
208;i_addr = WHERE( T GE Tf )
209;Ew(i_addr) = cw * T(i_addr)
210;i_addr = FINDGEN(N/2+1)
211;
212;Ei(i_addr) = c0 * ( T(i_addr) + mu * S ) - L * ( 1 + mu * S / T(i_addr) ) - cw*mu*S
213;Ei2(i_addr) = - ( 1 - e_i(i_addr) ) * L
214;
215;Emax = MAX(Ew(WHERE(Ew NE miss_val)))
216;
217;PLOT, [Tmin, Tmax], [ MIN(Ei)/1000., Emax/1000. ], /NODATA, TITLE = "Enthalpy", CHARSIZE = cs, XTITLE = "degC", YTITLE = "kJ/kg", YMINOR = 2
218;OPLOT, T, Ei/1000., MAX=0., THICK = 3.;, solid black
219;OPLOT, T, Ei2/1000., MAX=0., COLOR = 150, THICK = 3.; solid grey
220;;OPLOT, T, Es/1000., LINESTYLE = 1, THICK = 3.
221;OPLOT, T, Es_wgt/1000., COLOR = 100, THICK = 3
222;i_addr = WHERE( T GE 0. )
223;OPLOT, T(i_addr), Es_wgt(i_addr)/1000., COLOR = 200, THICK = 3
224;;OPLOT, T, Ew/1000., MAX=Emax/1000., LINESTYLE = 0, THICK = 3.
225;OPLOT, [ T_50, T_50], [ -400., 100.], COLOR = 100
226;OPLOT, [ Tmin, Tmax], [0., 0.], COLOR = 100
227;
228;;--- Sample experiment
229;Es_spl = c0 * Ts_spl - L
230;Ew_spl = cw * Tw_spl
231;
232;LOADCT, ct
233;FOR i = 0, N_spl - 1 DO BEGIN
234;   OPLOT, [ Ts_spl(i) ], [ rhos / rhoi * Es_spl(i) ]/1000., PSYM = 1, SYMSIZE = 2, THICK = 3, COLOR = COLORS(i)
235;ENDFOR
236;LOADCT, 0
237;OPLOT, [ Tw_spl ], [ (rhoi - rhos) / rhoi * Ew_spl ]/1000., PSYM = 1, SYMSIZE = 2, THICK = 3
238;
239;Ei_spl = Es_spl * rhos / rhoi + Ew_spl * ( rhoi - rhos ) / rhoi
240;
241;; Retrieve T
242;zA = ( L *( rhos/rhoi ) - (rhos/rhoi) * c0 * Ts_spl + (rhos-rhoi)/rhoi * cw * Tw_spl ) * rhoi
243;
244;aaa = c0
245;bbb = mu * S * ( c0 - cw ) - L + zA / rhoi
246;ccc = - L * mu * S
247;
248;discr = SQRT ( bbb * bbb - 4 * aaa * ccc )
249;
250;Ti_spl = ( - bbb - discr ) / 2.0 / aaa
251;
252;Ei_test = c0 * ( Ti_spl + mu * S ) - L * ( 1 + mu * S / Ti_spl ) - cw*mu*S
253;PRINT, "Ei_test : ", Ei_test
254;PRINT, "Ei_spl  : ", Ei_spl
255;
256;; Retrieve brine volume
257;e_i_spl = -mu*S/Ti_spl
258;PRINT, "e_i_spl = ", e_i_spl
259;
260;LOADCT, ct
261;FOR i = 0, N_spl - 1 DO BEGIN
262;   OPLOT, [ Ti_spl(i) ], [ Ei_spl(i) ]/1000., PSYM = 1, SYMSIZE = 2, COLOR = COLORS(i), THICK = 3
263;ENDFOR
264;LOADCT, 0
265
266;--------------------------------------------------------------------------------------------------
267DEVICE, /CLOSE
268SET_PLOT, 'X'
269END
Note: See TracBrowser for help on using the repository browser.