Index: trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F (revision 2388) +++ (revision ) @@ -1,540 +1,0 @@ -c======================================================================= - subroutine massflowrateco2(P,T,Sat,Radius,Matm,Ic) -c -c Determination of the mass transfer rate for CO2 condensation & -c sublimation -c -c inputs: Pressure (P), Temperature (T), saturation ratio (Sat), -c particle radius (Radius), molecular mass of the atm (Matm) -c output: MASS FLUX Ic -c -c Authors: C. Listowski (2014) then J. Audouard (2016-2017) -c -c -c Updates: -c -------- -c December 2017 - C. Listowski - Simplification of the derivation of -c massflowrate by using explicit formula for surface temperature, -c No Newton-Raphson routine anymore- see comment at relevant line -c======================================================================= - USE comcstfi_h, ONLY: pi - - implicit none - - include "microphys.h" - -c arguments: INPUT -c ---------- - REAL,INTENT(in) :: T,Matm - REAL*8,INTENT(in) :: SAT - REAL,INTENT(in) :: P - DOUBLE PRECISION,INTENT(in) :: Radius -c arguments: OUTPUT -c ---------- - DOUBLE PRECISION,INTENT(out) :: Ic -c Local Variables -c ---------- - DOUBLE PRECISION Tsurf - DOUBLE PRECISION C0,C1,C2 - DOUBLE PRECISION kmix,Lsub,cond - DOUBLE PRECISION Ak - - call coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2,Ak) - - Tsurf = 1./C1*dlog(Sat/Ak) + T - - !Note by CL - dec 2017 (see also technical note) - !The above is a simplified version of Tsurf - !compared to the one used by Listowski et al. 2014 (Ta), where a - !Newton-Raphson routine must be used. Approximations made by - !considering the orders of magnitude of the different factors lead to - !simplification of the equation 5 of Listowski et al. (2014). - !The error compared to the exact value determined by NR iterations - !is less than 0.6% for all sizes, pressures, supersaturations - !relevant to present Mars. Should also be ok for most conditions - !in ancient Mars (However, needs to be double-cheked, as explained in - !(Listowski et al. 2013,JGR) - - cond = 4.*pi*Radius*kmix - Ic = cond*(Tsurf-T)/Lsub - - END - -c******************************************************************************** - - subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2,Ak) - -c******************************************************************************** -c defini la fonction eq 6 papier 2014 (Listowski et al., 2014) - use tracer_mod, only: rho_ice_co2 - USE comcstfi_h, ONLY: pi - - implicit none - include "microphys.h" - -c arguments: INPUT -c ---------------- - REAL,INTENT(in) :: P - REAL,INTENT(in) :: T - REAL*8,INTENT(in) :: S - DOUBLE PRECISION,INTENT(in) :: rc - REAL,INTENT(in) :: Matm !g.mol-1 ( = mmean(ig,l) ) -c arguments: OUTPUT -c ---------- - DOUBLE PRECISION,INTENT(out) :: C0,C1,C2 - DOUBLE PRECISION,INTENT(out) :: kmix,Lsub - -c local: -c ------ - DOUBLE PRECISION Cpatm,Cpn2,Cpco2 - DOUBLE PRECISION psat, xinf, pco2 - DOUBLE PRECISION Dv - DOUBLE PRECISION l0,l1,l2,l3,l4 - DOUBLE PRECISION knudsen, a, lambda ! F and S correction - DOUBLE PRECISION Ak ! kelvin factor - DOUBLE PRECISION vthatm,lpmt,rhoatm, vthco2 ! for Kn,th - - -c DEFINE heat cap. J.kg-1.K-1 and To - - data Cpco2/0.7e3/ - data Cpn2/1e3/ - - kmix = 0d0 - Lsub = 0d0 - - C0 = 0d0 - C1 = 0d0 - C2 = 0d0 - -c Equilibirum pressure over a flat surface - psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T)) ! (Pa) -c Compute transport coefficient - pco2 = psat * dble(S) -c Latent heat of sublimation if CO2 co2 (J.kg-1) -c version Azreg_Ainou (J/kg) : - l0=595594. - l1=903.111 - l2=-11.5959 - l3=0.0528288 - l4=-0.000103183 - Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 * - & dble(T)**3 + l4 * dble(T)**4 ! J/kg -c atmospheric density - rhoatm = dble(P*Matm)/(rgp*dble(T)) ! g.m-3 - rhoatm = rhoatm * 1.00e-3 !kg.m-3 - call KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2 - call Diffcoeff(P, T, Dv) - - Dv = Dv * 1.00e-4 !!! cm2.s-1 to m2.s-1 - -c ----- FS correction for Diff - vthco2 = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1 - knudsen = 3*Dv / (vthco2 * rc) - lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black) - Dv = Dv / (1. + lambda * knudsen) -c ----- FS correction for Kth - vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg - Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1 - lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/ - & (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer - knudsen = lpmt / rc - lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black) - kmix = kmix / (1. + lambda * knudsen) -c --------------------- ASSIGN coeff values for FUNCTION - xinf = dble(S) * psat / dble(P) - Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) )) - C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/ - & (rgp*dble(T))) - C1 = Lsub*mco2/(rgp*dble(T)**2) - C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T)) - - END - - -c====================================================================== - subroutine Diffcoeff(P, T, Diff) -c Compute diffusion coefficient CO2/N2 -c cited in Ilona's lecture - from Reid et al. 1987 -c====================================================================== - IMPLICIT NONE - - include "microphys.h" - -c arguments -c ----------- - - REAL,INTENT(in) :: P - REAL,INTENT(in) :: T - -c output -c ----------- - - DOUBLE PRECISION,INTENT(out) :: Diff - -c local -c ----------- - - REAL Pbar !!! has to be in bar for the formula - DOUBLE PRECISION dva, dvb, Mab ! Mab has to be in g.mol-1 - - Pbar = P * 1d-5 - - Mab = 2. / ( 1./mn2 + 1./mco2 ) * 1000. - - dva = 26.9 ! diffusion volume of CO2, Reid et al. 1987 (cited in Ilona's lecture) - dvb = 18.5 ! diffusion volume of N2 - - Diff = 0.00143 * dble(T)**(1.75) / (dble(Pbar) * sqrt(Mab) - & * (dble(dva)**(1./3.) + dble(dvb)**(1./3.))**2.) !!! in cm2.s-1 - - RETURN - - END - - -c====================================================================== - - subroutine KthMixNEW(Kthmix,T,x,rho) - -c Compute thermal conductivity of CO2/N2 mixture -c (***WITHOUT*** USE OF VISCOSITY) - -c (Mason & Saxena, 1958 - Wassiljeva 1904) -c====================================================================== - - implicit none - - include "microphys.h" -c arguments -c ----------- - - REAL,INTENT(in) :: T - DOUBLE PRECISION,INTENT(in) :: x - DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3 - -c outputs -c ----------- - - DOUBLE PRECISION,INTENT(out) :: Kthmix - -c local -c ------------ - - DOUBLE PRECISION x1,x2 - - DOUBLE PRECISION Tc1, Tc2, Pc1, Pc2 - - DOUBLE PRECISION A12, A11, A22, A21 - - DOUBLE PRECISION Gamma1, Gamma2, M1, M2 - DOUBLE PRECISION lambda_trans1, lambda_trans2,epsilon - - DOUBLE PRECISION kco2, kn2 - - x1 = x - x2 = 1d0 - x - - M1 = mco2 - M2 = mn2 - - Tc1 = 304.1282 !(Scalabrin et al. 2006) - Tc2 = 126.192 ! (Lemmon & Jacobsen 2003) - - Pc1 = 73.773 ! (bars) - Pc2 = 33.958 ! (bars) - - Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.) - Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.) - -c Translational conductivities - - lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) ) - & /Gamma1 - - lambda_trans2 = ( exp(0.0464 * T/Tc2) - exp(-0.2412 * T/Tc2) ) - & /Gamma2 - -c Coefficient of Mason and Saxena - epsilon = 1. - - A11 = 1. - - A22 = 1. - - A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)* - & (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2)) - - A21 = epsilon * (1. + sqrt(lambda_trans2/lambda_trans1)* - & (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1)) - -c INDIVIDUAL COND. - - call KthCO2Scalab(kco2,T,rho) - call KthN2LemJac(kn2,T,rho) - -c MIXTURE COND. - Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22) - Kthmix = Kthmix*1e-3 ! from mW.m-1.K-1 to W.m-1.K-1 - - END - -c====================================================================== - subroutine KthN2LemJac(kthn2,T,rho) -c Compute thermal cond of N2 (Lemmon and Jacobsen, 2003) -cWITH viscosity -c====================================================================== - - implicit none - - include "microphys.h" -c include "microphysCO2.h" - - -c arguments -c ----------- - - REAL,INTENT(in) :: T - DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3 - -c outputs -c ----------- - - DOUBLE PRECISION,INTENT(out) :: kthn2 - -c local -c ------------ - - DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10 - DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10 - DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10 - DOUBLE PRECISION d4,d5,d6,d7,d8,d9 - DOUBLE PRECISION l4,l5,l6,l7,l8,l9 - DOUBLE PRECISION t2,t3,t4,t5,t6,t7,t8,t9 - DOUBLE PRECISION gamma4,gamma5,gamma6,gamma7,gamma8,gamma9 - - DOUBLE PRECISION Tc,rhoc - - DOUBLE PRECISION tau, delta - - DOUBLE PRECISION visco - - DOUBLE PRECISION k1, k2 - - N1 = 1.511d0 - N2 = 2.117d0 - N3 = -3.332d0 - - N4 = 8.862 - N5 = 31.11 - N6 = -73.13 - N7 = 20.03 - N8 = -0.7096 - N9 = 0.2672 - - t2 = -1.0d0 - t3 = -0.7d0 - t4 = 0.0d0 - t5 = 0.03 - t6 = 0.2 - t7 = 0.8 - t8 = 0.6 - t9 = 1.9 - - d4 = 1. - d5 = 2. - d6 = 3. - d7 = 4. - d8 = 8. - d9 = 10. - - l4 = 0. - gamma4 = 0. - - l5 = 0. - gamma5 = 0. - - l6 = 1. - gamma6 = 1. - - l7 = 2. - gamma7 = 1. - - l8 = 2. - gamma8 = 1. - - l9 = 2. - gamma9 = 1. - -c---------------------------------------------------------------------- - call viscoN2(T,visco) !! v given in microPa.s - - Tc = 126.192d0 - rhoc = 11.1839 * 1000 * mn2 !!!from mol.dm-3 to kg.m-3 - - tau = Tc / T - delta = rho/rhoc - - k1 = N1 * visco + N2 * tau**t2 + N3 * tau**t3 !!! mW m-1 K-1 -c--------- residual thermal conductivity - - k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4) - & + N5 * tau**t5 * delta**d5 * exp(-gamma5*delta**l5) - & + N6 * tau**t6 * delta**d6 * exp(-gamma6*delta**l6) - & + N7 * tau**t7 * delta**d7 * exp(-gamma7*delta**l7) - & + N8 * tau**t8 * delta**d8 * exp(-gamma8*delta**l8) - & + N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) - - kthn2 = k1 + k2 - - END - - -c====================================================================== - - subroutine viscoN2(T,visco) - -c Compute viscosity of N2 (Lemmon and Jacobsen, 2003) - -c====================================================================== - - implicit none - - include "microphys.h" -c include "microphysCO2.h" -c arguments -c ----------- - - REAL,INTENT(in) :: T - -c outputs -c ----------- - - DOUBLE PRECISION,INTENT(out) :: visco - - -c local -c ------------ - - DOUBLE PRECISION a0,a1,a2,a3,a4 - DOUBLE PRECISION Tstar,factor,sigma,M2 - DOUBLE PRECISION RGCS - - -c---------------------------------------------------------------------- - - - factor = 98.94 ! (K) - - sigma = 0.3656 ! (nm) - - a0 = 0.431 - a1 = -0.4623 - a2 = 0.08406 - a3 = 0.005341 - a4 = -0.00331 - - M2 = mn2 * 1.00e3 !!! to g.mol-1 - - Tstar = T*1./factor - - RGCS = exp( a0 + a1 * log(Tstar) + a2 * (log(Tstar))**2. + - & a3 * (log(Tstar))**3. + a4 * (log(Tstar))**4. ) - - - visco = 0.0266958 * sqrt(M2*T) / ( sigma**2. * RGCS ) !!! microPa.s - - - RETURN - - END - - -c====================================================================== - - subroutine KthCO2Scalab(kthco2,T,rho) - -c Compute thermal cond of CO2 (Scalabrin et al. 2006) - -c====================================================================== - - implicit none - - - -c arguments -c ----------- - - REAL,INTENT(in) :: T - DOUBLE PRECISION,INTENT(in) :: rho - -c outputs -c ----------- - - DOUBLE PRECISION,INTENT(out) :: kthco2 - -c LOCAL -c ----------- - - DOUBLE PRECISION Tc,Pc,rhoc, Lambdac - - DOUBLE PRECISION Tr, rhor, k1, k2 - - DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10 - DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10 - DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10 - - Tc = 304.1282 !(K) - Pc = 7.3773e6 !(MPa) - rhoc = 467.6 !(kg.m-3) - Lambdac = 4.81384 !(mW.m-1K-1) - - g1 = 0. - g2 = 0. - g3 = 1.5 - g4 = 0.0 - g5 = 1.0 - g6 = 1.5 - g7 = 1.5 - g8 = 1.5 - g9 = 3.5 - g10 = 5.5 - - h1 = 1. - h2 = 5. - h3 = 1. - h4 = 1. - h5 = 2. - h6 = 0. - h7 = 5.0 - h8 = 9.0 - h9 = 0. - h10 = 0. - - n1 = 7.69857587 - n2 = 0.159885811 - n3 = 1.56918621 - n4 = -6.73400790 - n5 = 16.3890156 - n6 = 3.69415242 - n7 = 22.3205514 - n8 = 66.1420950 - n9 = -0.171779133 - n10 = 0.00433043347 - - Tr = T/Tc - rhor = rho/rhoc - - k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) - & + n3*Tr**(g3)*rhor**(h3) - - k2 = n4*Tr**(g4)*rhor**(h4) + n5*Tr**(g5)*rhor**(h5) - & + n6*Tr**(g6)*rhor**(h6) + n7*Tr**(g7)*rhor**(h7) - & + n8*Tr**(g8)*rhor**(h8) + n9*Tr**(g9)*rhor**(h9) - & + n10*Tr**(g10)*rhor**(h10) - - k2 = exp(-5.*rhor**(2.)) * k2 - - kthco2 = (k1 + k2) * Lambdac ! mW - - END Index: trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F90 (revision 2389) +++ trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F90 (revision 2389) @@ -0,0 +1,618 @@ + subroutine massflowrateco2(P,T,Sat,Radius,Matm,Ic) +!======================================================================================================================! +! Determination of the mass transfer rate for CO2 condensation sublimation +! +! inputs: Pressure (P), Temperature (T), saturation ratio (Sat), particle radius (Radius), molecular mass of the atm +! (Matm) +! output: MASS FLUX Ic +! +! Authors: C. Listowski (2014), J. Audouard (2016-2017), Christophe Mathé (2020) +! +! Updates: +! -------- +! December 2017 - C. Listowski - Simplification of the derivation of massflowrate by using explicit formula for +! surface temperature, No Newton-Raphson routine anymore- see comment at relevant line +!======================================================================================================================! + use comcstfi_h, only: pi + + implicit none + + include "microphys.h" +!----------------------------------------------------------------------------------------------------------------------! +! VARIABLES DECLARATION +!----------------------------------------------------------------------------------------------------------------------! +! Input arguments: +!----------------- + real, intent(in) :: & + P, & + T, & + Matm + + real(kind=8), intent(in) :: & + SAT + + double precision, intent(in) :: & + Radius +!----------------------------------------------------------------------------------------------------------------------! +! Output arguments: +!------------------ + double precision, intent(out) :: & + Ic +!----------------------------------------------------------------------------------------------------------------------! +! Local: +!------- + double precision :: & + Ak, & + C0, & + C1, & + C2, & + kmix, & + Lsub, & + cond, & + Tsurf +!======================================================================================================================! +! BEGIN ===============================================================================================================! +!======================================================================================================================! + call coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2,Ak) + + Tsurf = 1./C1*dlog(Sat/Ak) + T + + ! Note by CL - dec 2017 (see also technical note) + ! The above is a simplified version of Tsurf compared to the one used by Listowski et al. 2014 (Ta), where a + ! Newton-Raphson routine must be used. Approximations made by considering the orders of magnitude of the different + ! factors lead to simplification of the equation 5 of Listowski et al. (2014). + ! The error compared to the exact value determined by NR iterations is less than 0.6% for all sizes, pressures, + ! supersaturations relevant to present Mars. Should also be ok for most conditions in ancient Mars (However, needs to + ! end subroutine massflowrateco2 be double-cheked, as explained in (Listowski et al. 2013,JGR) + + cond = 4.* pi * Radius * kmix + Ic = cond * (Tsurf-T) / Lsub +!======================================================================================================================! +! END =================================================================================================================! +!======================================================================================================================! + end subroutine massflowrateco2 + +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! + + + subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2,Ak) +!======================================================================================================================! +! Defini la fonction eq 6 papier 2014 (Listowski et al., 2014) +!======================================================================================================================! + use tracer_mod, only: rho_ice_co2 + use comcstfi_h, only: pi + + implicit none + include "microphys.h" +!----------------------------------------------------------------------------------------------------------------------! +! VARIABLES DECLARATION +!----------------------------------------------------------------------------------------------------------------------! +! Input arguments: +!----------------- + real, intent(in) :: & + P, & + T, & + Matm ! g.mol-1 ( = mmean(ig,l) ) + + real(kind=8), intent(in) :: & + S + + double precision, intent(in) :: & + rc +!----------------------------------------------------------------------------------------------------------------------! +! Output arguments: +!------------------ + double precision, intent(out) :: & + C0, & + C1, & + C2, & + kmix, & + Lsub +!----------------------------------------------------------------------------------------------------------------------! +! Local: +!------- + double precision :: & + Cpatm, & + Cpn2, & + Cpco2, & + Dv, & + psat, & + xinf, & + pco2, & + l0, & + l1, & + l2, & + l3, & + l4, & + Ak, & ! kelvin factor +!----F and S correction + knudsen,& + a, & + lambda, & +!----For Kn,th + vthatm, & + lpmt, & + rhoatm, & + vthco2 +!----------------------------------------------------------------------------------------------------------------------! +! DEFINE heat cap. J.kg-1.K-1 and To + data Cpco2/0.7e3/ + data Cpn2/1e3/ +!======================================================================================================================! +! BEGIN ===============================================================================================================! +!======================================================================================================================! + kmix = 0d0 + Lsub = 0d0 + + C0 = 0d0 + C1 = 0d0 + C2 = 0d0 + + ! Equilibirum pressure over a flat surface + psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T)) ! (Pa) + + ! Compute transport coefficient + pco2 = psat * dble(S) + + ! Latent heat of sublimation if CO2 co2 (J.kg-1) version Azreg_Ainou (J/kg) : + l0=595594. + l1=903.111 + l2=-11.5959 + l3=0.0528288 + l4=-0.000103183 + + Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 * dble(T)**3 + l4 * dble(T)**4 ! J/kg + + ! Atmospheric density + rhoatm = dble(P*Matm)/(rgp*dble(T)) ! g.m-3 + rhoatm = rhoatm * 1.00e-3 ! kg.m-3 + + ! Compute thermal cond of mixture co2/N2 + call KthMixNEW(kmix,T,pco2/dble(P),rhoatm) + + call Diffcoeff(P, T, Dv) + + Dv = Dv * 1.00e-4 ! cm2.s-1 to m2.s-1 + + ! ----- FS correction for Diff + vthco2 = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1 + knudsen = 3*Dv / (vthco2 * rc) + lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptee, Dahneke 1983? en fait si (Monschick&Black) + Dv = Dv / (1. + lambda * knudsen) + + ! ----- FS correction for Kth + vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg + Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1 + lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/(dble(Matm)*1.00e-3))) ! mean free path related to heat transfer + knudsen = lpmt / rc + lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptee, Dahneke 1983? en fait si (Monschick&Black) + kmix = kmix / (1. + lambda * knudsen) + + ! --------------------- ASSIGN coeff values for FUNCTION + xinf = dble(S) * psat / dble(P) + Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) )) + C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/(rgp*dble(T))) + C1 = Lsub*mco2/(rgp*dble(T)**2) + C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T)) +!======================================================================================================================! +! END =================================================================================================================! +!======================================================================================================================! + end subroutine coefffunc + +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! + + + subroutine Diffcoeff(P, T, Diff) +!======================================================================================================================! +! Compute diffusion coefficient CO2/N2 cited in Ilona's lecture - from Reid et al. 1987 +!======================================================================================================================! + implicit none + + include "microphys.h" +!----------------------------------------------------------------------------------------------------------------------! +! VARIABLES DECLARATION +!----------------------------------------------------------------------------------------------------------------------! +! Input arguments: +!----------------- + real, intent(in) :: & + P, & + T +!----------------------------------------------------------------------------------------------------------------------! +! Output arguments: +!------------------ + double precision, intent(out) :: & + Diff +!----------------------------------------------------------------------------------------------------------------------! +! Local: +!------- + real :: & + Pbar ! has to be in bar for the formula + + double precision :: & + dva, & + dvb, & + Mab ! Mab has to be in g.mol-1 +!======================================================================================================================! +! BEGIN ===============================================================================================================! +!======================================================================================================================! + Pbar = P * 1d-5 + + Mab = 2. / ( 1./mn2 + 1./mco2 ) * 1000. + + dva = 26.9 ! diffusion volume of CO2, Reid et al. 1987 (cited in Ilona's lecture) + dvb = 18.5 ! diffusion volume of N2 + + !!! Diff in cm2.s-1 + Diff = 0.00143 * dble(T)**(1.75) / (dble(Pbar) * sqrt(Mab) * (dble(dva)**(1./3.) + dble(dvb)**(1./3.))**2.) + + return +!======================================================================================================================! +! END =================================================================================================================! +!======================================================================================================================! + end subroutine Diffcoeff + + +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! + + + subroutine KthMixNEW(Kthmix,T,x,rho) +!======================================================================================================================! +! Compute thermal conductivity of CO2/N2 mixture (***WITHOUT*** USE OF VISCOSITY (Mason & Saxena 1958 - Wassiljeva 1904) +!======================================================================================================================! + implicit none + + include "microphys.h" +!----------------------------------------------------------------------------------------------------------------------! +! VARIABLES DECLARATION +!----------------------------------------------------------------------------------------------------------------------! +! Input arguments: +!----------------- + real, intent(in) :: & + T + + double precision, intent(in) :: & + x, & + rho ! kg.m-3 +!----------------------------------------------------------------------------------------------------------------------! +! Output arguments: +!------------------ + double precision, intent(out) :: & + Kthmix +!----------------------------------------------------------------------------------------------------------------------! +! Local: +!------- + double precision :: & + x1, & + x2, & + A12, & + A11, & + A22, & + A21, & + Tc1, & + Tc2, & + Pc1, & + Pc2, & + Gamma1, & + Gamma2, & + M1, & + M2, & + lambda_trans1, & + lambda_trans2, & + epsilon, & + kco2, & + kn2 +!======================================================================================================================! +! BEGIN ===============================================================================================================! +!======================================================================================================================! + x1 = x + x2 = 1d0 - x + + M1 = mco2 + M2 = mn2 + + Tc1 = 304.1282 ! (Scalabrin et al. 2006) + Tc2 = 126.192 ! (Lemmon & Jacobsen 2003) + + Pc1 = 73.773 ! (bars) + Pc2 = 33.958 ! (bars) + + Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.) + Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.) + + ! Translational conductivities + lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) ) / Gamma1 + lambda_trans2 = ( exp(0.0464 * T/Tc2) - exp(-0.2412 * T/Tc2) ) / Gamma2 + + ! Coefficient of Mason and Saxena + epsilon = 1. + + A11 = 1. + A22 = 1. + A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)* (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2)) + A21 = epsilon * (1. + sqrt(lambda_trans2/lambda_trans1)*(M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1)) + + ! INDIVIDUAL COND. + call KthCO2Scalab(kco2,T,rho) + call KthN2LemJac(kn2,T,rho) + + ! MIXTURE COND. + Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22) + Kthmix = Kthmix*1e-3 ! from mW.m-1.K-1 to W.m-1.K-1 +!======================================================================================================================! +! END =================================================================================================================! +!======================================================================================================================! + end subroutine KthMixNEW + + +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! + + + subroutine KthN2LemJac(kthn2,T,rho) +!======================================================================================================================! +! Compute thermal cond of N2 (Lemmon and Jacobsen, 2003) WITH viscosity +!======================================================================================================================! + implicit none + + include "microphys.h" +!----------------------------------------------------------------------------------------------------------------------! +! VARIABLES DECLARATION +!----------------------------------------------------------------------------------------------------------------------! +! Input arguments: +!----------------- + real, intent(in) :: T + + double precision, intent(in) :: rho ! kg.m-3 +!----------------------------------------------------------------------------------------------------------------------! +! Output argument: +!----------------- + double precision, intent(out) :: kthn2 +!----------------------------------------------------------------------------------------------------------------------! +! Local: +!------- + double precision :: & + g1, g2, g3, g4, g5, g6, g7, g8, g9, g10, & + h1, h2, h3, h4, h5, h6, h7, h8, h9, h10, & + n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, & + d4, d5, d6, d7, d8, d9, & + l4, l5, l6, l7, l8, l9, & + t2, t3, t4, t5, t6, t7, t8, t9, & + gamma4, gamma5, gamma6, gamma7, gamma8, gamma9, & + Tc, & + rhoc, & + tau, & + delta, & + visco, & + k1, & + k2 +!======================================================================================================================! +! BEGIN ===============================================================================================================! +!======================================================================================================================! + N1 = 1.511d0 + N2 = 2.117d0 + N3 = -3.332d0 + + N4 = 8.862 + N5 = 31.11 + N6 = -73.13 + N7 = 20.03 + N8 = -0.7096 + N9 = 0.2672 + + t2 = -1.0d0 + t3 = -0.7d0 + t4 = 0.0d0 + t5 = 0.03 + t6 = 0.2 + t7 = 0.8 + t8 = 0.6 + t9 = 1.9 + + d4 = 1. + d5 = 2. + d6 = 3. + d7 = 4. + d8 = 8. + d9 = 10. + + l4 = 0. + gamma4 = 0. + + l5 = 0. + gamma5 = 0. + + l6 = 1. + gamma6 = 1. + + l7 = 2. + gamma7 = 1. + + l8 = 2. + gamma8 = 1. + + l9 = 2. + gamma9 = 1. +!----------------------------------------------------------------------------------------------------------------------! + call viscoN2(T,visco) !! v given in microPa.s + + Tc = 126.192d0 + rhoc = 11.1839 * 1000 * mn2 ! from mol.dm-3 to kg.m-3 + + tau = Tc / T + delta = rho/rhoc + + k1 = N1 * visco + N2 * tau**t2 + N3 * tau**t3 !!! mW m-1 K-1 + + !--------- residual thermal conductivity + k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4) & + + N5 * tau**t5 * delta**d5 * exp(-gamma5*delta**l5) & + + N6 * tau**t6 * delta**d6 * exp(-gamma6*delta**l6) & + + N7 * tau**t7 * delta**d7 * exp(-gamma7*delta**l7) & + + N8 * tau**t8 * delta**d8 * exp(-gamma8*delta**l8) & + + N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) + + kthn2 = k1 + k2 +!======================================================================================================================! +! END =================================================================================================================! +!======================================================================================================================! + end subroutine KthN2LemJac + + +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! + + + subroutine viscoN2(T,visco) +!======================================================================================================================! +! Compute viscosity of N2 (Lemmon and Jacobsen, 2003) +!======================================================================================================================! + implicit none + + include "microphys.h" +!----------------------------------------------------------------------------------------------------------------------! +! VARIABLES DECLARATION +!----------------------------------------------------------------------------------------------------------------------! +! Input argument: +!---------------- + real, intent(in) :: T +!----------------------------------------------------------------------------------------------------------------------! +! Output argument: +!----------------- + double precision, intent(out) :: visco +!----------------------------------------------------------------------------------------------------------------------! +! Local: +!------- +!-----1) Parameters: +!------------------- + double precision, parameter :: & + factor = 98.94, & ! (K) + sigma = 0.3656, & ! (nm) + a0 = 0.431, & + a1 = -0.4623, & + a2 = 0.08406, & + a3 = 0.005341, & + a4 = -0.00331 +!----------------------------------------------------------------------------------------------------------------------! +!-----2) Variables: +!------------------ + double precision :: & + Tstar, & + M2, & + RGCS +!======================================================================================================================! +! BEGIN ===============================================================================================================! +!======================================================================================================================! + M2 = mn2 * 1.00e3 !!! to g.mol-1 + + Tstar = T*1./factor + + RGCS = exp( a0 + a1 * log(Tstar) + a2 * (log(Tstar))**2. + a3 * (log(Tstar))**3. + a4 * (log(Tstar))**4. ) + + visco = 0.0266958 * sqrt(M2*T) / ( sigma**2. * RGCS ) !!! microPa.s + + return +!======================================================================================================================! +! END =================================================================================================================! +!======================================================================================================================! + end subroutine viscoN2 + + +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! +!______________________________________________________________________________________________________________________! + + + subroutine KthCO2Scalab(kthco2,T,rho) +!======================================================================================================================! +! Compute thermal cond of CO2 (Scalabrin et al. 2006) +!======================================================================================================================! + implicit none +!----------------------------------------------------------------------------------------------------------------------! +! VARIABLES DECLARATION +!----------------------------------------------------------------------------------------------------------------------! +! Input arguments: +!----------------- + real, intent(in) :: T + + double precision, intent(in) :: rho +!----------------------------------------------------------------------------------------------------------------------! +! Output argument: +!----------------- + double precision, intent(out) :: kthco2 +!----------------------------------------------------------------------------------------------------------------------! +! Local: +!------- +!------- +!-----1) Parameters: +!------------------- + double precision, parameter :: & + Tc = 304.1282, & !(K) + Pc = 7.3773e6, & !(MPa) + rhoc = 467.6, & !(kg.m-3) + Lambdac = 4.81384,& !(mW.m-1K-1) + g1 = 0., & + g2 = 0., & + g3 = 1.5, & + g4 = 0.0, & + g5 = 1.0, & + g6 = 1.5, & + g7 = 1.5, & + g8 = 1.5, & + g9 = 3.5, & + g10 = 5.5, & + h1 = 1., & + h2 = 5., & + h3 = 1., & + h4 = 1., & + h5 = 2., & + h6 = 0., & + h7 = 5.0, & + h8 = 9.0, & + h9 = 0., & + h10 = 0., & + n1 = 7.69857587, & + n2 = 0.159885811, & + n3 = 1.56918621, & + n4 = -6.73400790, & + n5 = 16.3890156, & + n6 = 3.69415242, & + n7 = 22.3205514, & + n8 = 66.1420950, & + n9 = -0.171779133,& + n10 = 0.00433043347 +!----------------------------------------------------------------------------------------------------------------------! +!-----2) Variables: +!------------------ + double precision :: & + Tr, & + rhor,& + k1, & + k2 +!======================================================================================================================! +! BEGIN ===============================================================================================================! +!======================================================================================================================! + Tr = T/Tc + rhor = rho/rhoc + + k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) + n3*Tr**(g3)*rhor**(h3) + + k2 = n4*Tr**(g4)*rhor**(h4) + n5*Tr**(g5)*rhor**(h5) & + + n6*Tr**(g6)*rhor**(h6) + n7*Tr**(g7)*rhor**(h7) & + + n8*Tr**(g8)*rhor**(h8) + n9*Tr**(g9)*rhor**(h9) & + + n10*Tr**(g10)*rhor**(h10) + + k2 = k2 * exp(-5.*rhor**(2.)) + + kthco2 = (k1 + k2) * Lambdac ! mW +!======================================================================================================================! +! END =================================================================================================================! +!======================================================================================================================! + end subroutine KthCO2Scalab