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

Last change on this file since 3700 was 3008, checked in by emillour, 2 years ago

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