source: trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F90 @ 2800

Last change on this file since 2800 was 2389, checked in by cmathe, 4 years ago

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