Context Navigation

massflowrateco2.F @ 2350

Last change on this file since 2350 was 2156, checked in by aslmd, 6 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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format