source: trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F @ 2156

Last change on this file since 2156 was 2156, checked in by aslmd, 5 years ago

changed upper-case CO2 in lower-case co2 to be consistent everywhere. also Fortran is not case-sensitive so using a mix of upper case and lower case is not usually good practice.

  • Property svn:executable set to *
File size: 14.3 KB
Line 
1c=======================================================================
2      subroutine massflowrateco2(P,T,Sat,Radius,Matm,Ic)
3c
4c     Determination of the mass transfer rate for CO2 condensation &
5c     sublimation
6c
7c   inputs: Pressure (P), Temperature (T), saturation ratio (Sat),
8c           particle radius (Radius), molecular mass of the atm (Matm)
9c   output:  MASS FLUX Ic
10c
11c   Authors: C. Listowski (2014) then J. Audouard (2016-2017)
12
13c
14c   Updates:
15c   --------
16c   December 2017 - C. Listowski - Simplification of the derivation of
17c   massflowrate by using explicit formula for surface temperature,
18c   No Newton-Raphson routine anymore- see comment at relevant line
19c=======================================================================
20      USE comcstfi_h, ONLY: pi
21
22      implicit none
23
24      include "microphys.h"
25
26c   arguments: INPUT
27c   ----------
28      REAL,INTENT(in) :: T,Matm
29      REAL*8,INTENT(in) :: SAT
30      REAL,INTENT(in) :: P
31      DOUBLE PRECISION,INTENT(in) :: Radius
32c   arguments: OUTPUT
33c   ----------
34      DOUBLE PRECISION,INTENT(out) ::   Ic
35c   Local Variables
36c   ----------
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
63c********************************************************************************
64
65      subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2,Ak)
66
67c********************************************************************************
68c 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
75c   arguments: INPUT
76c   ----------------
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) )
82c   arguments: OUTPUT
83c   ----------
84      DOUBLE PRECISION,INTENT(out) ::  C0,C1,C2
85      DOUBLE PRECISION,INTENT(out) ::  kmix,Lsub
86
87c   local:
88c   ------
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
98c     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
110c     Equilibirum pressure over a flat surface
111      psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T))  ! (Pa)
112c     Compute transport coefficient
113      pco2 = psat * dble(S)
114c     Latent heat of sublimation if CO2  co2 (J.kg-1)
115c     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
123c     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
131c     ----- 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)
136c     ----- 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)
144c     --------------------- 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
155c======================================================================
156      subroutine Diffcoeff(P, T, Diff)
157c     Compute diffusion coefficient CO2/N2
158c     cited in Ilona's lecture - from Reid et al. 1987
159c======================================================================
160       IMPLICIT NONE
161
162       include "microphys.h"
163
164c      arguments
165c     -----------
166     
167      REAL,INTENT(in) :: P
168      REAL,INTENT(in) :: T
169     
170c     output
171c     -----------
172
173      DOUBLE PRECISION,INTENT(out) :: Diff
174   
175c      local
176c     -----------
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
196c======================================================================
197
198         subroutine KthMixNEW(Kthmix,T,x,rho)
199
200c        Compute thermal conductivity of CO2/N2 mixture
201c         (***WITHOUT*** USE OF VISCOSITY)
202
203c          (Mason & Saxena, 1958 - Wassiljeva 1904)
204c======================================================================
205
206       implicit none
207       
208       include "microphys.h"
209c      arguments
210c     -----------
211         
212         REAL,INTENT(in) :: T
213         DOUBLE PRECISION,INTENT(in) :: x
214         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
215
216c     outputs
217c     -----------
218
219         DOUBLE PRECISION,INTENT(out) :: Kthmix
220
221c     local
222c    ------------
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
250c 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     
258c     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
271c     INDIVIDUAL COND.
272
273         call KthCO2Scalab(kco2,T,rho)
274         call KthN2LemJac(kn2,T,rho)
275
276c     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
282c======================================================================
283         subroutine KthN2LemJac(kthn2,T,rho)
284c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
285cWITH viscosity
286c======================================================================
287
288       implicit none
289
290        include "microphys.h"
291c        include "microphysCO2.h"
292
293
294c      arguments
295c     -----------
296         
297         REAL,INTENT(in) :: T
298         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
299
300c     outputs
301c     -----------
302
303         DOUBLE PRECISION,INTENT(out) :: kthn2
304
305c     local
306c    ------------
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
369c----------------------------------------------------------------------           
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     
379c--------- 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
393c======================================================================
394
395         subroutine viscoN2(T,visco)
396
397c        Compute viscosity of N2 (Lemmon and Jacobsen, 2003)
398
399c======================================================================
400
401         implicit none
402
403       include "microphys.h"
404c       include "microphysCO2.h"
405c      arguments
406c     -----------
407         
408      REAL,INTENT(in) :: T
409
410c     outputs
411c     -----------
412
413      DOUBLE PRECISION,INTENT(out) :: visco
414
415
416c     local
417c    ------------
418
419      DOUBLE PRECISION a0,a1,a2,a3,a4
420      DOUBLE PRECISION Tstar,factor,sigma,M2
421      DOUBLE PRECISION RGCS
422     
423
424c---------------------------------------------------------------------- 
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
453c======================================================================
454
455         subroutine KthCO2Scalab(kthco2,T,rho)
456
457c        Compute thermal cond of CO2 (Scalabrin et al. 2006)
458
459c======================================================================
460
461         implicit none
462
463
464
465c      arguments
466c     -----------
467
468      REAL,INTENT(in) :: T
469      DOUBLE PRECISION,INTENT(in) :: rho
470
471c     outputs
472c     -----------
473
474      DOUBLE PRECISION,INTENT(out) :: kthco2
475
476c     LOCAL
477c     -----------
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
Note: See TracBrowser for help on using the repository browser.