1 | ; |
---|
2 | ;------------------------------------------------------------------------------ |
---|
3 | ; |
---|
4 | ; Plot LIM thermal properties |
---|
5 | ; |
---|
6 | ;------------------------------------------------------------------------------ |
---|
7 | |
---|
8 | ; Model parameters |
---|
9 | ; Sample experiments |
---|
10 | N_spl = 3 & Ts_spl = FLTARR(N_spl) |
---|
11 | Ts_spl[*] = [ -15., -10., -5. ] & Tw_spl = 2. |
---|
12 | colors = [ 50., 100., 150. ] |
---|
13 | |
---|
14 | ; |
---|
15 | ;============================================================================== |
---|
16 | ; 1) User parameters |
---|
17 | ;============================================================================== |
---|
18 | ; |
---|
19 | |
---|
20 | out_dir='' |
---|
21 | file_out = 'Thermal_properties' |
---|
22 | |
---|
23 | numplot_x = 2 ; number of horizontal plots |
---|
24 | numplot_y = 2 ; number of vertical plots |
---|
25 | ct = 13 ; colortable |
---|
26 | cs = 1.0 ; size of fonts on plots (1=normal) |
---|
27 | th = 2. ; thickness of curves |
---|
28 | device = 'PS' ; 'PS' or 'X' |
---|
29 | |
---|
30 | ; |
---|
31 | ;============================================================================== |
---|
32 | ; 2) Initialize graphics |
---|
33 | ;============================================================================== |
---|
34 | ; |
---|
35 | |
---|
36 | dx = 6 |
---|
37 | dy = 6. |
---|
38 | figuresize_x = numplot_x * dx; figure size on x direction |
---|
39 | figuresize_y = numplot_y * dy; figure size on y direction |
---|
40 | SET_PLOT, DEVICE |
---|
41 | ;!P.FONT=1 |
---|
42 | TVLCT, [0,255,0,0], [0,0,255,0], [0,0,0,255] |
---|
43 | IF ( 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 |
---|
48 | ENDIF |
---|
49 | |
---|
50 | IF ( device EQ 'X' ) THEN BEGIN |
---|
51 | xsize = 1200 |
---|
52 | ysize = 800 |
---|
53 | colorkey = 'rd' |
---|
54 | init_graphics_x, xsize, ysize, colorkey |
---|
55 | ENDIF |
---|
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 |
---|
66 | N_T = 200 |
---|
67 | N_S = 3 |
---|
68 | miss_val = 9.99e32 |
---|
69 | |
---|
70 | ;--- Physical constants |
---|
71 | c0 = 2067.d0 |
---|
72 | cw = 3991.d0 |
---|
73 | L = 334000.d0 |
---|
74 | mu = 0.054d0 |
---|
75 | rhoi = 917.d0 |
---|
76 | rhow = 1026.d0 |
---|
77 | rhos = 330.d0 |
---|
78 | |
---|
79 | Sw = 34.d0 |
---|
80 | |
---|
81 | beta1 = 0.09; pringle coefficients |
---|
82 | beta2 = 0.011; |
---|
83 | |
---|
84 | ;--- Ice salinity |
---|
85 | ;S = ( rhoi - rhos ) / rhoi * Sw |
---|
86 | S = FLTARR(N_S) |
---|
87 | Smin = 5. |
---|
88 | Smax = 15. |
---|
89 | dS = ( Smax - Smin ) / FLOAT(N_S-1) |
---|
90 | S[0:N_S-1] = Smin + FINDGEN(N_S) * dS |
---|
91 | |
---|
92 | Tf = -mu * S |
---|
93 | |
---|
94 | ;--- Temperatures |
---|
95 | T = FLTARR(N_T) |
---|
96 | Tmin = -20. |
---|
97 | Tf = -mu * Smax |
---|
98 | Tmax = 0. |
---|
99 | dT = ( Tmax - Tmin ) / FLOAT(N_T-1) |
---|
100 | T[0:N_T-1] = Tmin + FINDGEN(N_T) * dT |
---|
101 | |
---|
102 | ;--- Grids |
---|
103 | S_grid = FLTARR(N_S,N_T) |
---|
104 | T_grid = FLTARR(N_S,N_T) |
---|
105 | |
---|
106 | FOR j = 0, N_T - 1 DO BEGIN |
---|
107 | S_grid(0:N_S-1,j) = S(0:N_S-1) |
---|
108 | ENDFOR |
---|
109 | |
---|
110 | FOR i = 0, N_S - 1 DO BEGIN |
---|
111 | T_grid(i,0:N_T-1) = T(0:N_T-1) |
---|
112 | ENDFOR |
---|
113 | |
---|
114 | ;--- Brine fraction |
---|
115 | phi = -mu * S_grid / T_grid |
---|
116 | |
---|
117 | ;--- specific enthalpy |
---|
118 | q_i = c0 * ( T_grid + mu * S_grid ) - L * ( 1 + mu * S_grid / T_grid ) - cw*mu*S_grid |
---|
119 | |
---|
120 | ;--- thermal conductivity (pringle) |
---|
121 | kpri = 2.11 + beta1 * S_grid / T_grid - beta2 * T_grid; |
---|
122 | |
---|
123 | ;--- specific heat |
---|
124 | cp = c0 + L*mu*S_grid/(T_grid*T_grid) |
---|
125 | |
---|
126 | ;--- PLOTS |
---|
127 | !P.MULTI=[0, numplot_x, numplot_y] |
---|
128 | LOADCT, 13 |
---|
129 | |
---|
130 | PLOT, [Tmin, Tmax], [0., 50.], /NODATA, XTITLE='T (C)', YTITLE='phi (%)', CHARSIZE = cs |
---|
131 | FOR 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 |
---|
134 | ENDFOR |
---|
135 | |
---|
136 | PLOT, [Tmin, Tmax], [-4., 0.], /NODATA, XTITLE='T (C)', YTITLE='q_m (10^5 J/kg)', CHARSIZE = cs |
---|
137 | FOR 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 |
---|
140 | ENDFOR |
---|
141 | |
---|
142 | PLOT, [Tmin, Tmax], [0., 3.], /NODATA, XTITLE='T (C)', YTITLE='k_i (W/m/K)', CHARSIZE = cs |
---|
143 | FOR 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 |
---|
146 | ENDFOR |
---|
147 | |
---|
148 | PLOT, [Tmin, Tmax], [0., 20.], /NODATA, XTITLE='T (C)', YTITLE='c_i (10^3 J/kg/K)', CHARSIZE = cs |
---|
149 | FOR 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 |
---|
152 | ENDFOR |
---|
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 | ;-------------------------------------------------------------------------------------------------- |
---|
267 | DEVICE, /CLOSE |
---|
268 | SET_PLOT, 'X' |
---|
269 | END |
---|