1 | c======================================================================= |
---|
2 | subroutine massflowrateco2(P,T,Sat,Radius,Matm,Ic) |
---|
3 | c |
---|
4 | c Determination of the mass transfer rate for CO2 condensation & |
---|
5 | c sublimation |
---|
6 | c |
---|
7 | c inputs: Pressure (P), Temperature (T), saturation ratio (Sat), |
---|
8 | c particle radius (Radius), molecular mass of the atm (Matm) |
---|
9 | c output: MASS FLUX Ic |
---|
10 | c |
---|
11 | c Authors: C. Listowski (2014) then J. Audouard (2016-2017) |
---|
12 | c |
---|
13 | c |
---|
14 | c Updates: |
---|
15 | c -------- |
---|
16 | c December 2017 - C. Listowski - Simplification of the derivation of |
---|
17 | c massflowrate by using explicit formula for surface temperature, |
---|
18 | c No Newton-Raphson routine anymore- see comment at relevant line |
---|
19 | c======================================================================= |
---|
20 | USE comcstfi_h, ONLY: pi |
---|
21 | |
---|
22 | implicit none |
---|
23 | |
---|
24 | include "microphys.h" |
---|
25 | |
---|
26 | c arguments: INPUT |
---|
27 | c ---------- |
---|
28 | REAL,INTENT(in) :: T,Matm |
---|
29 | REAL*8,INTENT(in) :: SAT |
---|
30 | REAL,INTENT(in) :: P |
---|
31 | DOUBLE PRECISION,INTENT(in) :: Radius |
---|
32 | c arguments: OUTPUT |
---|
33 | c ---------- |
---|
34 | DOUBLE PRECISION,INTENT(out) :: Ic |
---|
35 | c Local Variables |
---|
36 | c ---------- |
---|
37 | DOUBLE PRECISION Tsurf |
---|
38 | DOUBLE PRECISION C0,C1,C2 |
---|
39 | DOUBLE PRECISION kmix,Lsub,cond |
---|
40 | DOUBLE PRECISION Ak |
---|
41 | |
---|
42 | call coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2,Ak) |
---|
43 | |
---|
44 | Tsurf = 1./C1*dlog(Sat/Ak) + T |
---|
45 | |
---|
46 | !Note by CL - dec 2017 (see also technical note) |
---|
47 | !The above is a simplified version of Tsurf |
---|
48 | !compared to the one used by Listowski et al. 2014 (Ta), where a |
---|
49 | !Newton-Raphson routine must be used. Approximations made by |
---|
50 | !considering the orders of magnitude of the different factors lead to |
---|
51 | !simplification of the equation 5 of Listowski et al. (2014). |
---|
52 | !The error compared to the exact value determined by NR iterations |
---|
53 | !is less than 0.6% for all sizes, pressures, supersaturations |
---|
54 | !relevant to present Mars. Should also be ok for most conditions |
---|
55 | !in ancient Mars (However, needs to be double-cheked, as explained in |
---|
56 | !(Listowski et al. 2013,JGR) |
---|
57 | |
---|
58 | cond = 4.*pi*Radius*kmix |
---|
59 | Ic = cond*(Tsurf-T)/Lsub |
---|
60 | |
---|
61 | END |
---|
62 | |
---|
63 | c******************************************************************************** |
---|
64 | |
---|
65 | subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2,Ak) |
---|
66 | |
---|
67 | c******************************************************************************** |
---|
68 | c defini la fonction eq 6 papier 2014 (Listowski et al., 2014) |
---|
69 | use tracer_mod, only: rho_ice_co2 |
---|
70 | USE comcstfi_h, ONLY: pi |
---|
71 | |
---|
72 | implicit none |
---|
73 | include "microphys.h" |
---|
74 | |
---|
75 | c arguments: INPUT |
---|
76 | c ---------------- |
---|
77 | REAL,INTENT(in) :: P |
---|
78 | REAL,INTENT(in) :: T |
---|
79 | REAL*8,INTENT(in) :: S |
---|
80 | DOUBLE PRECISION,INTENT(in) :: rc |
---|
81 | REAL,INTENT(in) :: Matm !g.mol-1 ( = mmean(ig,l) ) |
---|
82 | c arguments: OUTPUT |
---|
83 | c ---------- |
---|
84 | DOUBLE PRECISION,INTENT(out) :: C0,C1,C2 |
---|
85 | DOUBLE PRECISION,INTENT(out) :: kmix,Lsub |
---|
86 | |
---|
87 | c local: |
---|
88 | c ------ |
---|
89 | DOUBLE PRECISION Cpatm,Cpn2,Cpco2 |
---|
90 | DOUBLE PRECISION psat, xinf, pco2 |
---|
91 | DOUBLE PRECISION Dv |
---|
92 | DOUBLE PRECISION l0,l1,l2,l3,l4 |
---|
93 | DOUBLE PRECISION knudsen, a, lambda ! F and S correction |
---|
94 | DOUBLE PRECISION Ak ! kelvin factor |
---|
95 | DOUBLE PRECISION vthatm,lpmt,rhoatm, vthco2 ! for Kn,th |
---|
96 | |
---|
97 | |
---|
98 | c DEFINE heat cap. J.kg-1.K-1 and To |
---|
99 | |
---|
100 | data Cpco2/0.7e3/ |
---|
101 | data Cpn2/1e3/ |
---|
102 | |
---|
103 | kmix = 0d0 |
---|
104 | Lsub = 0d0 |
---|
105 | |
---|
106 | C0 = 0d0 |
---|
107 | C1 = 0d0 |
---|
108 | C2 = 0d0 |
---|
109 | |
---|
110 | c Equilibirum pressure over a flat surface |
---|
111 | psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T)) ! (Pa) |
---|
112 | c Compute transport coefficient |
---|
113 | pco2 = psat * dble(S) |
---|
114 | c Latent heat of sublimation if CO2 co2 (J.kg-1) |
---|
115 | c version Azreg_Ainou (J/kg) : |
---|
116 | l0=595594. |
---|
117 | l1=903.111 |
---|
118 | l2=-11.5959 |
---|
119 | l3=0.0528288 |
---|
120 | l4=-0.000103183 |
---|
121 | Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 * |
---|
122 | & dble(T)**3 + l4 * dble(T)**4 ! J/kg |
---|
123 | c atmospheric density |
---|
124 | rhoatm = dble(P*Matm)/(rgp*dble(T)) ! g.m-3 |
---|
125 | rhoatm = rhoatm * 1.00e-3 !kg.m-3 |
---|
126 | call KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2 |
---|
127 | call Diffcoeff(P, T, Dv) |
---|
128 | |
---|
129 | Dv = Dv * 1.00e-4 !!! cm2.s-1 to m2.s-1 |
---|
130 | |
---|
131 | c ----- FS correction for Diff |
---|
132 | vthco2 = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1 |
---|
133 | knudsen = 3*Dv / (vthco2 * rc) |
---|
134 | lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black) |
---|
135 | Dv = Dv / (1. + lambda * knudsen) |
---|
136 | c ----- FS correction for Kth |
---|
137 | vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg |
---|
138 | Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1 |
---|
139 | lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/ |
---|
140 | & (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer |
---|
141 | knudsen = lpmt / rc |
---|
142 | lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black) |
---|
143 | kmix = kmix / (1. + lambda * knudsen) |
---|
144 | c --------------------- ASSIGN coeff values for FUNCTION |
---|
145 | xinf = dble(S) * psat / dble(P) |
---|
146 | Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) )) |
---|
147 | C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/ |
---|
148 | & (rgp*dble(T))) |
---|
149 | C1 = Lsub*mco2/(rgp*dble(T)**2) |
---|
150 | C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T)) |
---|
151 | |
---|
152 | END |
---|
153 | |
---|
154 | |
---|
155 | c====================================================================== |
---|
156 | subroutine Diffcoeff(P, T, Diff) |
---|
157 | c Compute diffusion coefficient CO2/N2 |
---|
158 | c cited in Ilona's lecture - from Reid et al. 1987 |
---|
159 | c====================================================================== |
---|
160 | IMPLICIT NONE |
---|
161 | |
---|
162 | include "microphys.h" |
---|
163 | |
---|
164 | c arguments |
---|
165 | c ----------- |
---|
166 | |
---|
167 | REAL,INTENT(in) :: P |
---|
168 | REAL,INTENT(in) :: T |
---|
169 | |
---|
170 | c output |
---|
171 | c ----------- |
---|
172 | |
---|
173 | DOUBLE PRECISION,INTENT(out) :: Diff |
---|
174 | |
---|
175 | c local |
---|
176 | c ----------- |
---|
177 | |
---|
178 | REAL Pbar !!! has to be in bar for the formula |
---|
179 | DOUBLE PRECISION dva, dvb, Mab ! Mab has to be in g.mol-1 |
---|
180 | |
---|
181 | Pbar = P * 1d-5 |
---|
182 | |
---|
183 | Mab = 2. / ( 1./mn2 + 1./mco2 ) * 1000. |
---|
184 | |
---|
185 | dva = 26.9 ! diffusion volume of CO2, Reid et al. 1987 (cited in Ilona's lecture) |
---|
186 | dvb = 18.5 ! diffusion volume of N2 |
---|
187 | |
---|
188 | Diff = 0.00143 * dble(T)**(1.75) / (dble(Pbar) * sqrt(Mab) |
---|
189 | & * (dble(dva)**(1./3.) + dble(dvb)**(1./3.))**2.) !!! in cm2.s-1 |
---|
190 | |
---|
191 | RETURN |
---|
192 | |
---|
193 | END |
---|
194 | |
---|
195 | |
---|
196 | c====================================================================== |
---|
197 | |
---|
198 | subroutine KthMixNEW(Kthmix,T,x,rho) |
---|
199 | |
---|
200 | c Compute thermal conductivity of CO2/N2 mixture |
---|
201 | c (***WITHOUT*** USE OF VISCOSITY) |
---|
202 | |
---|
203 | c (Mason & Saxena, 1958 - Wassiljeva 1904) |
---|
204 | c====================================================================== |
---|
205 | |
---|
206 | implicit none |
---|
207 | |
---|
208 | include "microphys.h" |
---|
209 | c arguments |
---|
210 | c ----------- |
---|
211 | |
---|
212 | REAL,INTENT(in) :: T |
---|
213 | DOUBLE PRECISION,INTENT(in) :: x |
---|
214 | DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3 |
---|
215 | |
---|
216 | c outputs |
---|
217 | c ----------- |
---|
218 | |
---|
219 | DOUBLE PRECISION,INTENT(out) :: Kthmix |
---|
220 | |
---|
221 | c local |
---|
222 | c ------------ |
---|
223 | |
---|
224 | DOUBLE PRECISION x1,x2 |
---|
225 | |
---|
226 | DOUBLE PRECISION Tc1, Tc2, Pc1, Pc2 |
---|
227 | |
---|
228 | DOUBLE PRECISION A12, A11, A22, A21 |
---|
229 | |
---|
230 | DOUBLE PRECISION Gamma1, Gamma2, M1, M2 |
---|
231 | DOUBLE PRECISION lambda_trans1, lambda_trans2,epsilon |
---|
232 | |
---|
233 | DOUBLE PRECISION kco2, kn2 |
---|
234 | |
---|
235 | x1 = x |
---|
236 | x2 = 1d0 - x |
---|
237 | |
---|
238 | M1 = mco2 |
---|
239 | M2 = mn2 |
---|
240 | |
---|
241 | Tc1 = 304.1282 !(Scalabrin et al. 2006) |
---|
242 | Tc2 = 126.192 ! (Lemmon & Jacobsen 2003) |
---|
243 | |
---|
244 | Pc1 = 73.773 ! (bars) |
---|
245 | Pc2 = 33.958 ! (bars) |
---|
246 | |
---|
247 | Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.) |
---|
248 | Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.) |
---|
249 | |
---|
250 | c Translational conductivities |
---|
251 | |
---|
252 | lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) ) |
---|
253 | & /Gamma1 |
---|
254 | |
---|
255 | lambda_trans2 = ( exp(0.0464 * T/Tc2) - exp(-0.2412 * T/Tc2) ) |
---|
256 | & /Gamma2 |
---|
257 | |
---|
258 | c Coefficient of Mason and Saxena |
---|
259 | epsilon = 1. |
---|
260 | |
---|
261 | A11 = 1. |
---|
262 | |
---|
263 | A22 = 1. |
---|
264 | |
---|
265 | A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)* |
---|
266 | & (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2)) |
---|
267 | |
---|
268 | A21 = epsilon * (1. + sqrt(lambda_trans2/lambda_trans1)* |
---|
269 | & (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1)) |
---|
270 | |
---|
271 | c INDIVIDUAL COND. |
---|
272 | |
---|
273 | call KthCO2Scalab(kco2,T,rho) |
---|
274 | call KthN2LemJac(kn2,T,rho) |
---|
275 | |
---|
276 | c MIXTURE COND. |
---|
277 | Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22) |
---|
278 | Kthmix = Kthmix*1e-3 ! from mW.m-1.K-1 to W.m-1.K-1 |
---|
279 | |
---|
280 | END |
---|
281 | |
---|
282 | c====================================================================== |
---|
283 | subroutine KthN2LemJac(kthn2,T,rho) |
---|
284 | c Compute thermal cond of N2 (Lemmon and Jacobsen, 2003) |
---|
285 | cWITH viscosity |
---|
286 | c====================================================================== |
---|
287 | |
---|
288 | implicit none |
---|
289 | |
---|
290 | include "microphys.h" |
---|
291 | c include "microphysCO2.h" |
---|
292 | |
---|
293 | |
---|
294 | c arguments |
---|
295 | c ----------- |
---|
296 | |
---|
297 | REAL,INTENT(in) :: T |
---|
298 | DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3 |
---|
299 | |
---|
300 | c outputs |
---|
301 | c ----------- |
---|
302 | |
---|
303 | DOUBLE PRECISION,INTENT(out) :: kthn2 |
---|
304 | |
---|
305 | c local |
---|
306 | c ------------ |
---|
307 | |
---|
308 | DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10 |
---|
309 | DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10 |
---|
310 | DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10 |
---|
311 | DOUBLE PRECISION d4,d5,d6,d7,d8,d9 |
---|
312 | DOUBLE PRECISION l4,l5,l6,l7,l8,l9 |
---|
313 | DOUBLE PRECISION t2,t3,t4,t5,t6,t7,t8,t9 |
---|
314 | DOUBLE PRECISION gamma4,gamma5,gamma6,gamma7,gamma8,gamma9 |
---|
315 | |
---|
316 | DOUBLE PRECISION Tc,rhoc |
---|
317 | |
---|
318 | DOUBLE PRECISION tau, delta |
---|
319 | |
---|
320 | DOUBLE PRECISION visco |
---|
321 | |
---|
322 | DOUBLE PRECISION k1, k2 |
---|
323 | |
---|
324 | N1 = 1.511d0 |
---|
325 | N2 = 2.117d0 |
---|
326 | N3 = -3.332d0 |
---|
327 | |
---|
328 | N4 = 8.862 |
---|
329 | N5 = 31.11 |
---|
330 | N6 = -73.13 |
---|
331 | N7 = 20.03 |
---|
332 | N8 = -0.7096 |
---|
333 | N9 = 0.2672 |
---|
334 | |
---|
335 | t2 = -1.0d0 |
---|
336 | t3 = -0.7d0 |
---|
337 | t4 = 0.0d0 |
---|
338 | t5 = 0.03 |
---|
339 | t6 = 0.2 |
---|
340 | t7 = 0.8 |
---|
341 | t8 = 0.6 |
---|
342 | t9 = 1.9 |
---|
343 | |
---|
344 | d4 = 1. |
---|
345 | d5 = 2. |
---|
346 | d6 = 3. |
---|
347 | d7 = 4. |
---|
348 | d8 = 8. |
---|
349 | d9 = 10. |
---|
350 | |
---|
351 | l4 = 0. |
---|
352 | gamma4 = 0. |
---|
353 | |
---|
354 | l5 = 0. |
---|
355 | gamma5 = 0. |
---|
356 | |
---|
357 | l6 = 1. |
---|
358 | gamma6 = 1. |
---|
359 | |
---|
360 | l7 = 2. |
---|
361 | gamma7 = 1. |
---|
362 | |
---|
363 | l8 = 2. |
---|
364 | gamma8 = 1. |
---|
365 | |
---|
366 | l9 = 2. |
---|
367 | gamma9 = 1. |
---|
368 | |
---|
369 | c---------------------------------------------------------------------- |
---|
370 | call viscoN2(T,visco) !! v given in microPa.s |
---|
371 | |
---|
372 | Tc = 126.192d0 |
---|
373 | rhoc = 11.1839 * 1000 * mn2 !!!from mol.dm-3 to kg.m-3 |
---|
374 | |
---|
375 | tau = Tc / T |
---|
376 | delta = rho/rhoc |
---|
377 | |
---|
378 | k1 = N1 * visco + N2 * tau**t2 + N3 * tau**t3 !!! mW m-1 K-1 |
---|
379 | c--------- residual thermal conductivity |
---|
380 | |
---|
381 | k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4) |
---|
382 | & + N5 * tau**t5 * delta**d5 * exp(-gamma5*delta**l5) |
---|
383 | & + N6 * tau**t6 * delta**d6 * exp(-gamma6*delta**l6) |
---|
384 | & + N7 * tau**t7 * delta**d7 * exp(-gamma7*delta**l7) |
---|
385 | & + N8 * tau**t8 * delta**d8 * exp(-gamma8*delta**l8) |
---|
386 | & + N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) |
---|
387 | |
---|
388 | kthn2 = k1 + k2 |
---|
389 | |
---|
390 | END |
---|
391 | |
---|
392 | |
---|
393 | c====================================================================== |
---|
394 | |
---|
395 | subroutine viscoN2(T,visco) |
---|
396 | |
---|
397 | c Compute viscosity of N2 (Lemmon and Jacobsen, 2003) |
---|
398 | |
---|
399 | c====================================================================== |
---|
400 | |
---|
401 | implicit none |
---|
402 | |
---|
403 | include "microphys.h" |
---|
404 | c include "microphysCO2.h" |
---|
405 | c arguments |
---|
406 | c ----------- |
---|
407 | |
---|
408 | REAL,INTENT(in) :: T |
---|
409 | |
---|
410 | c outputs |
---|
411 | c ----------- |
---|
412 | |
---|
413 | DOUBLE PRECISION,INTENT(out) :: visco |
---|
414 | |
---|
415 | |
---|
416 | c local |
---|
417 | c ------------ |
---|
418 | |
---|
419 | DOUBLE PRECISION a0,a1,a2,a3,a4 |
---|
420 | DOUBLE PRECISION Tstar,factor,sigma,M2 |
---|
421 | DOUBLE PRECISION RGCS |
---|
422 | |
---|
423 | |
---|
424 | c---------------------------------------------------------------------- |
---|
425 | |
---|
426 | |
---|
427 | factor = 98.94 ! (K) |
---|
428 | |
---|
429 | sigma = 0.3656 ! (nm) |
---|
430 | |
---|
431 | a0 = 0.431 |
---|
432 | a1 = -0.4623 |
---|
433 | a2 = 0.08406 |
---|
434 | a3 = 0.005341 |
---|
435 | a4 = -0.00331 |
---|
436 | |
---|
437 | M2 = mn2 * 1.00e3 !!! to g.mol-1 |
---|
438 | |
---|
439 | Tstar = T*1./factor |
---|
440 | |
---|
441 | RGCS = exp( a0 + a1 * log(Tstar) + a2 * (log(Tstar))**2. + |
---|
442 | & a3 * (log(Tstar))**3. + a4 * (log(Tstar))**4. ) |
---|
443 | |
---|
444 | |
---|
445 | visco = 0.0266958 * sqrt(M2*T) / ( sigma**2. * RGCS ) !!! microPa.s |
---|
446 | |
---|
447 | |
---|
448 | RETURN |
---|
449 | |
---|
450 | END |
---|
451 | |
---|
452 | |
---|
453 | c====================================================================== |
---|
454 | |
---|
455 | subroutine KthCO2Scalab(kthco2,T,rho) |
---|
456 | |
---|
457 | c Compute thermal cond of CO2 (Scalabrin et al. 2006) |
---|
458 | |
---|
459 | c====================================================================== |
---|
460 | |
---|
461 | implicit none |
---|
462 | |
---|
463 | |
---|
464 | |
---|
465 | c arguments |
---|
466 | c ----------- |
---|
467 | |
---|
468 | REAL,INTENT(in) :: T |
---|
469 | DOUBLE PRECISION,INTENT(in) :: rho |
---|
470 | |
---|
471 | c outputs |
---|
472 | c ----------- |
---|
473 | |
---|
474 | DOUBLE PRECISION,INTENT(out) :: kthco2 |
---|
475 | |
---|
476 | c LOCAL |
---|
477 | c ----------- |
---|
478 | |
---|
479 | DOUBLE PRECISION Tc,Pc,rhoc, Lambdac |
---|
480 | |
---|
481 | DOUBLE PRECISION Tr, rhor, k1, k2 |
---|
482 | |
---|
483 | DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10 |
---|
484 | DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10 |
---|
485 | DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10 |
---|
486 | |
---|
487 | Tc = 304.1282 !(K) |
---|
488 | Pc = 7.3773e6 !(MPa) |
---|
489 | rhoc = 467.6 !(kg.m-3) |
---|
490 | Lambdac = 4.81384 !(mW.m-1K-1) |
---|
491 | |
---|
492 | g1 = 0. |
---|
493 | g2 = 0. |
---|
494 | g3 = 1.5 |
---|
495 | g4 = 0.0 |
---|
496 | g5 = 1.0 |
---|
497 | g6 = 1.5 |
---|
498 | g7 = 1.5 |
---|
499 | g8 = 1.5 |
---|
500 | g9 = 3.5 |
---|
501 | g10 = 5.5 |
---|
502 | |
---|
503 | h1 = 1. |
---|
504 | h2 = 5. |
---|
505 | h3 = 1. |
---|
506 | h4 = 1. |
---|
507 | h5 = 2. |
---|
508 | h6 = 0. |
---|
509 | h7 = 5.0 |
---|
510 | h8 = 9.0 |
---|
511 | h9 = 0. |
---|
512 | h10 = 0. |
---|
513 | |
---|
514 | n1 = 7.69857587 |
---|
515 | n2 = 0.159885811 |
---|
516 | n3 = 1.56918621 |
---|
517 | n4 = -6.73400790 |
---|
518 | n5 = 16.3890156 |
---|
519 | n6 = 3.69415242 |
---|
520 | n7 = 22.3205514 |
---|
521 | n8 = 66.1420950 |
---|
522 | n9 = -0.171779133 |
---|
523 | n10 = 0.00433043347 |
---|
524 | |
---|
525 | Tr = T/Tc |
---|
526 | rhor = rho/rhoc |
---|
527 | |
---|
528 | k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) |
---|
529 | & + n3*Tr**(g3)*rhor**(h3) |
---|
530 | |
---|
531 | k2 = n4*Tr**(g4)*rhor**(h4) + n5*Tr**(g5)*rhor**(h5) |
---|
532 | & + n6*Tr**(g6)*rhor**(h6) + n7*Tr**(g7)*rhor**(h7) |
---|
533 | & + n8*Tr**(g8)*rhor**(h8) + n9*Tr**(g9)*rhor**(h9) |
---|
534 | & + n10*Tr**(g10)*rhor**(h10) |
---|
535 | |
---|
536 | k2 = exp(-5.*rhor**(2.)) * k2 |
---|
537 | |
---|
538 | kthco2 = (k1 + k2) * Lambdac ! mW |
---|
539 | |
---|
540 | END |
---|