Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sulfate_aer_mod.f90 @ 5310

Last change on this file since 5310 was 5268, checked in by abarral, 2 weeks ago
.f90 <-> .F90 depending on cpp key use
Property svn:keywords set to `Id`
File size: 31.8 KB

Rev	Line
[2690]	1	MODULE sulfate_aer_mod
	2
	3	! microphysical routines based on UPMC aerosol model by Slimane Bekki
	4	! adapted for stratospheric sulfate aerosol in LMDZ by Christoph Kleinschmitt
	5
	6	CONTAINS
	7
[4750]	8	!*******************************************************************
[4950]	9	SUBROUTINE STRACOMP_KELVIN(sh,t_seri,pplay)
[4750]	10	!
[4950]	11	! Aerosol H2SO4 weight fraction as a function of PH2O and temperature
	12	! INPUT:
	13	! sh: MMR of H2O
	14	! t_seri: temperature (K)
	15	! pplay: middle layer pression (Pa)
[4750]	16	!
[4950]	17	! Modified in modules:
	18	! R2SO4: aerosol H2SO4 weight fraction (percent)
	19	! R2SO4B: aerosol H2SO4 weight fraction (percent) for each aerosol bin
	20	! DENSO4: aerosol density (gr/cm3)
	21	! DENSO4B: aerosol density (gr/cm3)for each aerosol bin
	22	! f_r_wet: factor for converting dry to wet radius
	23	! assuming 'flat surface' composition (does not depend on aerosol size)
	24	! f_r_wetB: factor for converting dry to wet radius
	25	! assuming 'curved surface' composition (depends on aerosol size)
[4750]	26
[4950]	27	USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands
	28	USE infotrac_phy, ONLY : nbtr_bin
	29	USE aerophys
	30	USE phys_local_var_mod, ONLY: R2SO4, R2SO4B, DENSO4, DENSO4B, f_r_wet, f_r_wetB
	31	USE strataer_local_var_mod, ONLY: RRSI
	32	! WARNING: in phys_local_var_mod R2SO4B, DENSO4B, f_r_wetB (klon,klev,nbtr_bin)
	33	! and dens_aer_dry must be declared somewhere
[4750]	34
[4950]	35	IMPLICIT NONE
[4750]	36
[4950]	37	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	38	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression in the middle of each layer (Pa)
	39	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! specific humidity (kg h2o/kg air)
	40
	41	! local variables
	42	integer :: ilon,ilev,ik
	43	real, parameter :: rath2oair = mAIRmol/mH2Omol
	44	real, parameter :: third = 1./3.
	45	real :: pph2ogas(klon,klev)
	46	real :: temp, wpp, xa, surtens, mvh2o, radwet, fkelvin, pph2okel, r2so4ik, denso4ik
	47	!----------------------------------------
	48
	49	! gas-phase h2o partial pressure (Pa)
	50	! vmr=sh*rath2oair
	51	pph2ogas(:,:) = pplay(:,:)sh(:,:)rath2oair
[4750]	52
[4950]	53	DO ilon=1,klon
	54	DO ilev=1,klev
[4750]	55
[4950]	56	temp = max(t_seri(ilon,ilev),190.)
	57	temp = min(temp,300.)
	58
	59	! * H2SO4-H2O flat surface *
	60	!! equilibrium H2O pressure over pure flat liquid water (Pa)
	61	!! pflath2o=psh2o(temp)
	62	! h2so4 weight percent(%) = f(P_h2o(Pa),T)
	63	R2SO4(ilon,ilev)=wph2so4(pph2ogas(ilon,ilev),temp)
	64	! h2so4 mass fraction (0<wpp<1)
	65	wpp=R2SO4(ilon,ilev)*1.e-2
	66	! mole fraction
	67	xa=18.wpp/(18.wpp+98.*(1.-wpp))
	68
	69	! CHECK:compare h2so4 sat/ pressure (see Marti et al., 97 & reef. therein)
	70	! R2SO4(ilon,ilev)=70. temp=298.15
	71	! equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
	72	! include conversion from molec/cm3 to Pa
	73	! ph2so4=solh2so4(temp,xa)(1.38065e-16temp)/10.
	74	! print*,' ph2so4=',ph2so4,temp,R2SO4(ilon,ilev)
	75	! good match with Martin, et Ayers, not with Gmitro (the famous 0.086)
	76
	77	! surface tension (mN/m=1.e-3.kg/s2) = f(T,h2so4 mole fraction)
	78	surtens=surftension(temp,xa)
	79	! molar volume of pure h2o (cm3.mol-1 =1.e-6.m3.mol-1)
	80	mvh2o= rmvh2o(temp)
	81	! aerosol density (gr/cm3) = f(T,h2so4 mass fraction)
	82	DENSO4(ilon,ilev)=density(temp,wpp)
	83	! ->x1000., to have it in kg/m3
	84	! factor for converting dry to wet radius
	85	f_r_wet(ilon,ilev) = (dens_aer_dry/(DENSO4(ilon,ilev)*1.e3)/ &
	86	& (R2SO4(ilon,ilev)1.e-2))*third
	87	! * End of H2SO4-H2O flat surface *
	88
	89
	90	! Loop on bin radius (RRSI in cm)
	91	DO IK=1,nbtr_bin
	92
	93	! * H2SO4-H2O curved surface - Kelvin effect factor *
	94	! wet radius (m) (RRSI(IK) in [cm])
	95	if (f_r_wetB(ilon,ilev,IK) .gt. 1.0) then
	96	radwet = 1.e-2RRSI(IK)f_r_wetB(ilon,ilev,IK)
	97	else
	98	! H2SO4-H2O flat surface, only on the first timestep
	99	radwet = 1.e-2RRSI(IK)f_r_wet(ilon,ilev)
	100	endif
	101	! Kelvin factor:
	102	! surface tension (mN/m=1.e-3.kg/s2)
	103	! molar volume of pure h2o (cm3.mol-1 =1.e-6.m3.mol-1)
	104	fkelvin=exp( 2.1.e-3surtens1.e-6mvh2o/ (radwetrgastemp) )
	105	! equilibrium: pph2o(gas) = pph2o(liq) = pph2o(liq_flat) * fkelvin
	106	! equilibrium: pph2o(liq_flat) = pph2o(gas) / fkelvin
	107	! h2o liquid partial pressure before Kelvin effect (Pa)
	108	pph2okel = pph2ogas(ilon,ilev) / fkelvin
	109	! h2so4 weight percent(%) = f(P_h2o(Pa),temp)
	110	r2so4ik=wph2so4(pph2okel,temp)
	111	! h2so4 mass fraction (0<wpp<1)
	112	wpp=r2so4ik*1.e-2
	113	! mole fraction
	114	xa=18.wpp/(18.wpp+98.*(1.-wpp))
	115	! aerosol density (gr/cm3) = f(T,h2so4 mass fraction)
	116	denso4ik=density(temp,wpp)
	117	!
	118	! recalculate Kelvin factor with surface tension and radwet
	119	! with new R2SO4B and DENSO4B
	120	surtens=surftension(temp,xa)
	121	! wet radius (m)
	122	radwet = 1.e-2RRSI(IK)(dens_aer_dry/(denso4ik*1.e3)/ &
	123	& (r2so4ik1.e-2))*third
	124	fkelvin=exp( 2.1.e-3surtens1.e-6mvh2o / (radwetrgastemp) )
	125	pph2okel=pph2ogas(ilon,ilev) / fkelvin
	126	! h2so4 weight percent(%) = f(P_h2o(Pa),temp)
	127	R2SO4B(ilon,ilev,IK)=wph2so4(pph2okel,temp)
	128	! h2so4 mass fraction (0<wpp<1)
	129	wpp=R2SO4B(ilon,ilev,IK)*1.e-2
	130	xa=18.wpp/(18.wpp+98.*(1.-wpp))
	131	! aerosol density (gr/cm3) = f(T,h2so4 mass fraction)
	132	DENSO4B(ilon,ilev,IK)=density(temp,wpp)
	133	! factor for converting dry to wet radius
	134	f_r_wetB(ilon,ilev,IK) = (dens_aer_dry/(DENSO4B(ilon,ilev,IK)*1.e3)/ &
	135	& (R2SO4B(ilon,ilev,IK)1.e-2))*third
	136	!
	137	! print*,'R,Rwet(m),kelvin,h2so4(%),ro=',RRSI(ik),radwet,fkelvin, &
	138	! & R2SO4B(ilon,ilev,IK),DENSO4B(ilon,ilev,IK)
	139	! print*,' equil.h2so4(molec/cm3), &
	140	! & sigma',solh2so4(temp,xa),surftension(temp,xa)
	141
	142	ENDDO
	143
	144	ENDDO
	145	ENDDO
	146
	147	RETURN
[4750]	148
[4950]	149	END SUBROUTINE STRACOMP_KELVIN
[2690]	150	!********************************************************************
	151	SUBROUTINE STRACOMP(sh,t_seri,pplay)
	152
	153	! AEROSOL H2SO4 WEIGHT FRACTION AS A FUNCTION OF PH2O AND TEMPERATURE
	154	! ----------------------------------------------------------------
	155	! INPUT:
	156	! H2O: VMR of H2O
	157	! t_seri: temperature (K)
	158	! PMB: pressure (mb)
	159	! klon: number of latitude bands in the model domain
	160	! klev: number of altitude bands in the model domain
	161	! for IFS: perhaps add another dimension for longitude
	162	!
	163	! OUTPUT:
	164	! R2SO4: aerosol H2SO4 weight fraction (percent)
	165
	166	USE dimphy, ONLY : klon,klev
	167	USE aerophys
	168	USE phys_local_var_mod, ONLY: R2SO4
	169
	170	IMPLICIT NONE
	171
	172	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	173	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
	174	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
	175
	176	REAL PMB(klon,klev), H2O(klon,klev)
	177	!
	178	! working variables
	179	INTEGER I,J,K
	180	REAL TP, PH2O, VAL, A, B
	181	! local variables to be saved on exit
	182	INTEGER INSTEP
	183	INTEGER, PARAMETER :: N=16, M=28
	184	DATA INSTEP/0/
	185	REAL F(N,M)
	186	REAL XC(N)
	187	REAL YC(M)
	188	REAL XC1, XC16, YC1, YC28
	189	!
	190	SAVE INSTEP,F,XC,YC,XC1,XC16,YC1,YC28
[3663]	191	!$OMP THREADPRIVATE(INSTEP,F,XC,YC,XC1,XC16,YC1,YC28)
[2690]	192
	193	! convert pplay (in Pa) to PMB (in mb)
	194	PMB(:,:)=pplay(:,:)/100.0
	195
	196	! convert specific humidity sh (in kg/kg) to VMR H2O
	197	H2O(:,:)=sh(:,:)*mAIRmol/mH2Omol
	198
	199	IF(INSTEP.EQ.0) THEN
	200
	201	INSTEP=1
	202	XC(1)=0.01
	203	XC(2)=0.1
	204	XC(3)=0.5
	205	XC(4)=1.0
	206	XC(5)=1.5
	207	XC(6)=2.0
	208	XC(7)=3.0
	209	XC(8)=5.0
	210	XC(9)=6.0
	211	XC(10)=8.0
	212	XC(11)=10.0
	213	XC(12)=12.0
	214	XC(13)=15.0
	215	XC(14)=20.0
	216	XC(15)=30.0
	217	XC(16)=100.0
	218	!
	219	YC(1)=175.0
	220	DO I=2,28
	221	YC(I)=YC(I-1)+5.0
	222	ENDDO
	223
	224	! CONVERSION mb IN 1.0E-4mB
	225	DO I=1,16
	226	XC(I)=XC(I)*1.0E-4
	227	ENDDO
	228	!
	229	XC1=XC(1)+1.E-10
	230	XC16=XC(16)-1.E-8
	231	YC1=YC(1)+1.E-5
	232	YC28=YC(28)-1.E-5
	233
	234	F(6,4)=43.45
	235	F(6,5)=53.96
	236	F(6,6)=60.62
	237	F(6,7)=65.57
	238	F(6,8)=69.42
	239	F(6,9)=72.56
	240	F(6,10)=75.17
	241	F(6,11)=77.38
	242	F(6,12)=79.3
	243	F(6,13)=80.99
	244	F(6,14)=82.5
	245	F(6,15)=83.92
	246	F(6,16)=85.32
	247	F(6,17)=86.79
	248	F(6,18)=88.32
	249	!
	250	! ADD FACTOR BECAUSE THE SLOP IS TOO IMPORTANT
	251	! NOT FOR THIS ONE BUT THE REST
	252	! LOG DOESN'T WORK
	253	A=(F(6,5)-F(6,4))/( (YC(5)-YC(4))*2.0)
	254	B=-A*YC(4) + F(6,4)
	255	F(6,1)=A*YC(1) + B
	256	F(6,2)=A*YC(2) + B
	257	F(6,3)=A*YC(3) + B
	258	!
	259	F(7,4)=37.02
	260	F(7,5)=49.46
	261	F(7,6)=57.51
	262	F(7,7)=63.12
	263	F(7,8)=67.42
	264	F(7,9)=70.85
	265	F(7,10)=73.70
	266	F(7,11)=76.09
	267	F(7,12)=78.15
	268	F(7,13)=79.96
	269	F(7,14)=81.56
	270	F(7,15)=83.02
	271	F(7,16)=84.43
	272	F(7,17)=85.85
	273	F(7,18)=87.33
	274	!
	275	A=(F(7,5)-F(7,4))/( (YC(5)-YC(4))*2.0)
	276	B=-A*YC(4) + F(7,4)
	277	F(7,1)=A*YC(1) + B
	278	F(7,2)=A*YC(2) + B
	279	F(7,3)=A*YC(3) + B
	280	!
	281	F(8,4)=25.85
	282	F(8,5)=42.26
	283	F(8,6)=52.78
	284	F(8,7)=59.55
	285	F(8,8)=64.55
	286	F(8,9)=68.45
	287	F(8,10)=71.63
	288	F(8,11)=74.29
	289	F(8,12)=76.56
	290	F(8,13)=78.53
	291	F(8,14)=80.27
	292	F(8,15)=81.83
	293	F(8,16)=83.27
	294	F(8,17)=84.67
	295	F(8,18)=86.10
	296	!
	297	A=(F(8,5)-F(8,4))/( (YC(5)-YC(4))*2.5 )
	298	B=-A*YC(4) + F(8,4)
	299	F(8,1)=A*YC(1) + B
	300	F(8,2)=A*YC(2) + B
	301	F(8,3)=A*YC(3) + B
	302	!
	303	F(9,4)=15.38
	304	F(9,5)=39.35
	305	F(9,6)=50.73
	306	F(9,7)=58.11
	307	F(9,8)=63.41
	308	F(9,9)=67.52
	309	F(9,10)=70.83
	310	F(9,11)=73.6
	311	F(9,12)=75.95
	312	F(9,13)=77.98
	313	F(9,14)=79.77
	314	F(9,15)=81.38
	315	F(9,16)=82.84
	316	F(9,17)=84.25
	317	F(9,18)=85.66
	318	!
	319	A=(F(9,5)-F(9,4))/( (YC(5)-YC(4))*7.0)
	320	B=-A*YC(4) + F(9,4)
	321	F(9,1)=A*YC(1) + B
	322	F(9,2)=A*YC(2) + B
	323	F(9,3)=A*YC(3) + B
	324	!
	325	F(10,4)=0.0
	326	F(10,5)=34.02
	327	F(10,6)=46.93
	328	F(10,7)=55.61
	329	F(10,8)=61.47
	330	F(10,9)=65.94
	331	F(10,10)=69.49
	332	F(10,11)=72.44
	333	F(10,12)=74.93
	334	F(10,13)=77.08
	335	F(10,14)=78.96
	336	F(10,15)=80.63
	337	F(10,16)=82.15
	338	F(10,17)=83.57
	339	F(10,18)=84.97
	340	!
	341	A=(F(10,6)-F(10,5))/( (YC(6)-YC(5))*1.5)
	342	B=-A*YC(5) + F(10,5)
	343	F(10,1)=A*YC(1) + B
	344	F(10,2)=A*YC(2) + B
	345	F(10,3)=A*YC(3) + B
	346	F(10,4)=A*YC(4) + B
	347	!
	348	F(11,4)=0.0
	349	F(11,5)=29.02
	350	F(11,6)=43.69
	351	F(11,7)=53.44
	352	F(11,8)=59.83
	353	F(11,9)=64.62
	354	F(11,10)=68.39
	355	F(11,11)=71.48
	356	F(11,12)=74.10
	357	F(11,13)=76.33
	358	F(11,14)=78.29
	359	F(11,15)=80.02
	360	F(11,16)=81.58
	361	F(11,17)=83.03
	362	F(11,18)=84.44
	363	!
	364	A=(F(11,6)-F(11,5))/( (YC(6)-YC(5))*2.5 )
	365	B=-A*YC(5) + F(11,5)
	366	F(11,1)=A*YC(1) + B
	367	F(11,2)=A*YC(2) + B
	368	F(11,3)=A*YC(3) + B
	369	F(11,4)=A*YC(4) + B
	370	!
	371	F(12,4)=0.0
	372	F(12,5)=23.13
	373	F(12,6)=40.86
	374	F(12,7)=51.44
	375	F(12,8)=58.38
	376	F(12,9)=63.47
	377	F(12,10)=67.43
	378	F(12,11)=70.66
	379	F(12,12)=73.38
	380	F(12,13)=75.70
	381	F(12,14)=77.72
	382	F(12,15)=79.51
	383	F(12,16)=81.11
	384	F(12,17)=82.58
	385	F(12,18)=83.99
	386	!
	387	A=(F(12,6)-F(12,5))/( (YC(6)-YC(5))*3.5 )
	388	B=-A*YC(5) + F(12,5)
	389	F(12,1)=A*YC(1) + B
	390	F(12,2)=A*YC(2) + B
	391	F(12,3)=A*YC(3) + B
	392	F(12,4)=A*YC(4) + B
	393	!
	394	F(13,4)=0.0
	395	F(13,5)=0.0
	396	F(13,6)=36.89
	397	F(13,7)=48.63
	398	F(13,8)=56.46
	399	F(13,9)=61.96
	400	F(13,10)=66.19
	401	F(13,11)=69.6
	402	F(13,12)=72.45
	403	F(13,13)=74.89
	404	F(13,14)=76.99
	405	F(13,15)=78.85
	406	F(13,16)=80.50
	407	F(13,17)=82.02
	408	F(13,18)=83.44
	409	!
	410	A=(F(13,7)-F(13,6))/( (YC(7)-YC(6))*2.0)
	411	B=-A*YC(6) + F(13,6)
	412	F(13,1)=A*YC(1) + B
	413	F(13,2)=A*YC(2) + B
	414	F(13,3)=A*YC(3) + B
	415	F(13,4)=A*YC(4) + B
	416	F(13,5)=A*YC(5) + B
	417	!
	418	F(14,4)=0.0
	419	F(14,5)=0.0
	420	F(14,6)=30.82
	421	F(14,7)=44.49
	422	F(14,8)=53.69
	423	F(14,9)=59.83
	424	F(14,10)=64.47
	425	F(14,11)=68.15
	426	F(14,12)=71.19
	427	F(14,13)=73.77
	428	F(14,14)=76.0
	429	F(14,15)=77.95
	430	F(14,16)=79.69
	431	F(14,17)=81.26
	432	F(14,18)=82.72
	433	!
	434	A=(F(14,7)-F(14,6))/( (YC(7)-YC(6))*2.5 )
	435	B=-A*YC(6) + F(14,6)
	436	F(14,1)=A*YC(1) + B
	437	F(14,2)=A*YC(2) + B
	438	F(14,3)=A*YC(3) + B
	439	F(14,4)=A*YC(4) + B
	440	F(14,5)=A*YC(5) + B
	441	!
	442	F(15,4)=0.0
	443	F(15,5)=0.0
	444	F(15,6)=0.0
	445	F(15,7)=37.71
	446	F(15,8)=48.49
	447	F(15,9)=56.40
	448	F(15,10)=61.75
	449	F(15,11)=65.89
	450	F(15,12)=69.25
	451	F(15,13)=72.07
	452	F(15,14)=74.49
	453	F(15,15)=76.59
	454	F(15,16)=78.45
	455	F(15,17)=80.12
	456	F(15,18)=81.64
	457	!
	458	A=(F(15,8)-F(15,7))/( (YC(8)-YC(7))*1.5)
	459	B=-A*YC(7) + F(15,7)
	460	F(15,1)=A*YC(1) + B
	461	F(15,2)=A*YC(2) + B
	462	F(15,3)=A*YC(3) + B
	463	F(15,4)=A*YC(4) + B
	464	F(15,5)=A*YC(5) + B
	465	F(15,6)=A*YC(6) + B
	466
	467	! SUPPOSE THAT AT GIVEN AND PH2O<2mB,
	468	! %H2SO4 = A *LOG(PH2O) +B
	469	! XC(1-5) :EXTENSION LEFT (LOW H2O)
	470	DO J=1,18
	471	A=(F(6,J)-F(7,J))/(LOG(XC(6))-LOG(XC(7)))
	472	B=-A*LOG(XC(6)) + F(6,J)
	473	DO K=1,5
	474	F(K,J)=A*LOG(XC(K)) + B
	475	ENDDO
	476	ENDDO
	477
	478	! XC(16) :EXTENSION RIGHT (HIGH H2O)
	479	DO J=1,18
	480	A=(F(15,J)-F(14,J))/(XC(15)-XC(14))
	481	B=-A*XC(15) + F(15,J)
	482	F(16,J)=A*XC(16) + B
	483	! F(16,2)=1.0
	484	ENDDO
	485
	486	! YC(16-25) :EXTENSION DOWN (HIGH T)
	487	DO I=1,16
	488	A=(F(I,18)-F(I,17))/(YC(18)-YC(17))
	489	B=-A*YC(18) + F(I,18)
	490	DO K=19,28
	491	F(I,K)=A*YC(K) + B
	492	ENDDO
	493	ENDDO
	494
	495	! MANUAL CORRECTIONS
	496	DO J=1,10
	497	F(1,J)=94.0
	498	ENDDO
	499
	500	DO J=1,6
	501	F(2,J)=77.0 +REAL(J)
	502	ENDDO
	503
	504	DO J=1,7
	505	F(16,J)=9.0
	506	ENDDO
	507
	508	DO I=1,16
	509	DO J=1,28
	510	IF (F(I,J).LT.9.0) F(I,J)=30.0
	511	IF (F(I,J).GT.99.99) F(I,J)=99.99
	512	ENDDO
	513	ENDDO
	514
	515	ENDIF
	516
	517	DO I=1,klon
	518	DO J=1,klev
	519	TP=t_seri(I,J)
	520	IF (TP.LT.175.1) TP=175.1
	521	! Partial pressure of H2O (mb)
	522	PH2O =PMB(I,J)*H2O(I,J)
	523	IF (PH2O.LT.XC1) THEN
	524	R2SO4(I,J)=99.99
	525	! PH2O=XC(1)+1.0E-10
	526	ELSE
	527	IF (PH2O.GT.XC16) PH2O=XC16
	528	! SIMPLE LINEAR INTERPOLATIONS
	529	CALL FIND(PH2O,TP,XC,YC,F,VAL,N,M)
	530	IF (PMB(I,J).GE.10.0.AND.VAL.LT.60.0) VAL=60.0
	531	R2SO4(I,J)=VAL
	532	ENDIF
	533	ENDDO
	534	ENDDO
	535
	536	END SUBROUTINE
	537
	538	!****************************************************************
	539	SUBROUTINE STRAACT(ACTSO4)
	540
	541	! H2SO4 ACTIVITY (GIAUQUE) AS A FUNCTION OF H2SO4 WP
	542	! ----------------------------------------
	543	! INPUT:
	544	! H2SO4: VMR of H2SO4
	545	! klon: number of latitude bands in the model domain
	546	! klev: number of altitude bands in the model domain
	547	! for IFS: perhaps add another dimension for longitude
	548	!
	549	! OUTPUT:
	550	! ACTSO4: H2SO4 activity (percent)
	551
	552	USE dimphy, ONLY : klon,klev
	553	USE phys_local_var_mod, ONLY: R2SO4
	554
	555	IMPLICIT NONE
	556
	557	REAL ACTSO4(klon,klev)
	558
	559	! Working variables
	560	INTEGER NN,I,J,JX,JX1
	561	REAL TC,TB,TA,XT
	562	PARAMETER (NN=109)
	563	REAL XC(NN), X(NN)
	564
	565	! H2SO4 activity
	566	DATA X/ &
	567	& 0.0,0.25,0.78,1.437,2.19,3.07,4.03,5.04,6.08 &
	568	& ,7.13,8.18,14.33,18.59,28.59,39.17,49.49 &
	569	& ,102.4,157.8,215.7,276.9,341.6,409.8,481.5,556.6 &
	570	& ,635.5,719.,808.,902.,1000.,1103.,1211.,1322.,1437.,1555. &
	571	& ,1677.,1800.,1926.,2054.,2183.,2312.,2442.,2572.,2701.,2829. &
	572	& ,2955.,3080.,3203.,3325.,3446.,3564.,3681.,3796.,3910.,4022. &
	573	& ,4134.,4351.,4564.,4771.,4974.,5171.,5364.,5551.,5732.,5908. &
	574	& ,6079.,6244.,6404.,6559.,6709.,6854.,6994.,7131.,7264.,7393. &
	575	& ,7520.,7821.,8105.,8373.,8627.,8867.,9093.,9308.,9511.,9703. &
	576	& ,9885.,10060.,10225.,10535.,10819.,11079.,11318.,11537. &
	577	& ,11740.,12097.,12407.,12676.,12915.,13126.,13564.,13910. &
	578	& ,14191.,14423.,14617.,14786.,10568.,15299.,15491.,15654. &
	579	& ,15811./
	580	! H2SO4 weight fraction (percent)
	581	DATA XC/ &
	582	& 100.0,99.982,99.963,99.945,99.927,99.908,99.890,99.872 &
	583	& ,99.853,99.835,99.817,99.725,99.634,99.452,99.270 &
	584	& ,99.090,98.196,97.319,96.457,95.610,94.777,93.959,93.156 &
	585	& ,92.365,91.588,90.824,90.073,89.334,88.607,87.892,87.188 &
	586	& ,86.495,85.814,85.143,84.482,83.832,83.191,82.560,81.939 &
	587	& ,81.327,80.724,80.130,79.545,78.968,78.399,77.839,77.286 &
	588	& ,76.741,76.204,75.675,75.152,74.637,74.129,73.628,73.133 &
	589	& ,72.164,71.220,70.300,69.404,68.530,67.678,66.847,66.037 &
	590	& ,65.245,64.472,63.718,62.981,62.261,61.557,60.868,60.195 &
	591	& ,59.537,58.893,58.263,57.646,56.159,54.747,53.405,52.126 &
	592	& ,50.908,49.745,48.634,47.572,46.555,45.580,44.646,43.749 &
	593	& ,42.059,40.495,39.043,37.691,36.430,35.251,33.107,31.209 &
	594	& ,29.517,27.999,26.629,23.728,21.397,19.482,17.882,16.525 &
	595	& ,15.360,13.461,11.980,10.792,9.819,8.932/
	596
	597	DO I=1,klon
	598	DO J=1,klev
	599	! HERE LINEAR INTERPOLATIONS
	600	XT=R2SO4(I,J)
	601	CALL POSACT(XT,XC,NN,JX)
	602	JX1=JX+1
	603	IF(JX.EQ.0) THEN
	604	ACTSO4(I,J)=0.0
	605	ELSE IF(JX.GE.NN) THEN
	606	ACTSO4(I,J)=15811.0
	607	ELSE
	608	TC=XT -XC(JX)
	609	TB=X(JX1) -X(JX)
	610	TA=XC(JX1) -XC(JX)
	611	TA=TB/TA
	612	ACTSO4(I,J)=X(JX) + TA*TC
	613	ENDIF
	614	ENDDO
	615	ENDDO
	616
	617	END SUBROUTINE
[4750]	618
	619	!****************************************************************
[2690]	620	SUBROUTINE DENH2SA(t_seri)
	621
	622	! AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T
	623	! ---------------------------------------------
	624	! VERY ROUGH APPROXIMATION (SEE FOR WATER IN HANDBOOK
	625	! LINEAR 2% FOR 30 DEGREES with RESPECT TO WATER)
	626	!
	627	! INPUT:
	628	! R2SO4: aerosol H2SO4 weight fraction (percent)
	629	! t_seri: temperature (K)
	630	! klon: number of latitude bands in the model domain
	631	! klev: number of altitude bands in the model domain
	632	! for IFS: perhaps add another dimension for longitude
	633	!
	634	! OUTPUT:
	635	! DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)
	636	!
	637	USE dimphy, ONLY : klon,klev
	638	USE phys_local_var_mod, ONLY: R2SO4, DENSO4
	639
	640	IMPLICIT NONE
	641
	642	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	643
	644	INTEGER I,J
	645
	646	! Loop on model domain (2 dimension for UPMC model; 3 for IFS)
	647	DO I=1,klon
	648	DO J=1,klev
	649	! RO AT 20C
	650	DENSO4(I,J)=0.78681252E-5R2SO4(I,J)R2SO4(I,J)+ 0.82185978E-2*R2SO4(I,J)+0.97968381
	651	DENSO4(I,J)=DENSO4(I,J)* ( 1.0 - (t_seri(I,J)-293.0)*0.02/30.0 )
	652	ENDDO
	653	ENDDO
	654
	655	END SUBROUTINE
	656
	657	!***********************************************************
	658	SUBROUTINE FIND(X,Y,XC,YC,F,VAL,N,M)
	659	!
	660	! BI-LINEAR INTERPOLATION
	661
	662	! INPUT:
	663	! X: Partial pressure of H2O (mb)
	664	! Y: temperature (K)
	665	! XC: Table partial pressure of H2O (mb)
	666	! YC: Table temperature (K)
	667	! F: Table aerosol H2SO4 weight fraction=f(XC,YC) (percent)
	668	!
	669	! OUTPUT:
	670	! VAL: aerosol H2SO4 weight fraction (percent)
	671
	672	IMPLICIT NONE
	673
	674	INTEGER N,M
	675	REAL X,Y,XC(N),YC(M),F(N,M),VAL
	676	!
	677	! working variables
	678	INTEGER IERX,IERY,JX,JY,JXP1,JYP1
	679	REAL SXY,SX1Y,SX1Y1,SXY1,TA,TB,T,UA,UB,U
	680
	681	IERX=0
	682	IERY=0
	683	CALL POSITION(XC,X,N,JX,IERX)
	684	CALL POSITION(YC,Y,M,JY,IERY)
	685
	686	IF(JX.EQ.0.OR.IERY.EQ.1) THEN
	687	VAL=99.99
	688	RETURN
	689	ENDIF
	690
	691	IF(JY.EQ.0.OR.IERX.EQ.1) THEN
	692	VAL=9.0
	693	RETURN
	694	ENDIF
	695
	696	JXP1=JX+1
	697	JYP1=JY+1
	698	SXY=F(JX, JY )
	699	SX1Y=F(JXP1,JY )
	700	SX1Y1=F(JXP1,JYP1)
	701	SXY1=F(JX, JYP1)
	702
	703	! x-slope.
	704	TA=X -XC(JX)
	705	TB=XC(JXP1)-XC(JX)
	706	T=TA/TB
	707
	708	! y-slope.
	709	UA=Y -YC(JY)
	710	UB=YC(JYP1)-YC(JY)
	711	U=UA/UB
	712
	713	! Use bilinear interpolation to determine function at point X,Y.
	714	VAL=(1.-T)(1.-U)SXY + T(1.0-U)SX1Y + TUSX1Y1 + (1.0-T)USXY1
	715
	716	IF(VAL.LT.9.0) VAL=9.0
	717	IF(VAL.GT.99.99) VAL=99.99
	718
	719	RETURN
	720	END SUBROUTINE
	721	!****************************************************************
	722	SUBROUTINE POSITION(XC,X,N,JX,IER)
	723
	724	IMPLICIT NONE
	725
	726	INTEGER N,JX,IER,I
	727	REAL X,XC(N)
	728
	729	IER=0
	730	IF(X.LT.XC(1)) THEN
	731	JX=0
	732	ELSE
	733	DO 10 I=1,N
	734	IF (X.LT.XC(I)) GO TO 20
	735	10 CONTINUE
	736	IER=1
	737	20 JX=I-1
	738	ENDIF
	739
	740	RETURN
	741	END SUBROUTINE
	742	!********************************************************************
	743	SUBROUTINE POSACT(XT,X,N,JX)
	744
	745	! POSITION OF XT IN THE ARRAY X
	746	! -----------------------------------------------
	747
	748	IMPLICIT NONE
	749
	750	INTEGER N
	751	REAL XT,X(N)
	752	! Working variables
	753	INTEGER JX,I
	754
	755	IF(XT.GT.X(1)) THEN
	756	JX=0
	757	ELSE
	758	DO 10 I=1,N
	759	IF (XT.GT.X(I)) GO TO 20
	760	10 CONTINUE
	761	20 JX=I
	762	ENDIF
	763
	764	RETURN
	765	END SUBROUTINE
[4950]	766	!********************************************************************
	767	!-----------------------------------------------------------------------
	768	real function psh2so4(T) result(psh2so4_out)
	769	! equilibrium H2SO4 pressure over pure H2SO4 solution (Pa)
	770	!
	771	!---->Ayers et.al. (1980), GRL (7) pp 433-436
	772	! plus corrections for lower temperatures by Kulmala and Laaksonen (1990)
	773	! and Noppel et al. (1990)
[2690]	774
[4950]	775	implicit none
	776	real, intent(in) :: T
	777	real, parameter :: &
	778	& b1=1.01325e5, &
	779	& b2=11.5, &
	780	& b3=1.0156e4, &
	781	& b4=0.38/545., &
	782	& tref=360.15
	783
	784	! saturation vapor pressure ( N/m2 = Pa = kg/(m.s2) )
	785	psh2so4_out=b1exp( -b2 +b3( 1./tref-1./T &
	786	& +b4*(1.+log(tref/T)-tref/T) ) )
	787
	788	return
	789	end function psh2so4
	790	!-----------------------------------------------------------------------
	791	real function ndsh2so4(T) result(ndsh2so4_out)
	792	! equilibrium H2SO4 number density over pure H2SO4 (molec/cm3)
	793
	794	implicit none
	795	real, intent(in) :: T
	796	real :: presat
	797
	798	! Boltzmann constant ( 1.38065e-23 J/K = m2⋅kg/(s2⋅K) )
	799	! akb idem in cm2⋅g/(s2⋅K)
	800	real, parameter :: akb=1.38065e-16
	801
	802	! pure h2so4 saturation vapor pressure (Pa)
	803	presat=psh2so4(T)
	804	! saturation number density (1/cm3) - (molec/cm3)
	805	ndsh2so4_out=presat10./(akbT)
	806
	807	return
	808	end function ndsh2so4
	809	!-----------------------------------------------------------------------
	810	real function psh2o(T) result(psh2o_out)
	811	! equilibrium H2O pressure over pure liquid water (Pa)
	812	!
	813	implicit none
	814	real, intent(in) :: T
	815
	816	if(T.gt.229.) then
	817	! Preining et al., 1981 (from Kulmala et al., 1998)
	818	! saturation vapor pressure (N/m2 = 1 Pa = 1 kg/(m·s2))
	819	psh2o_out=exp( 77.34491296 -7235.424651/T &
	820	& -8.2log(T) + 5.7133e-3T )
	821	else
	822	! Tabazadeh et al., 1997, parameterization for 185<T<260
	823	! saturation water vapor partial pressure (mb = hPa =1.E2 kg/(m·s2))
	824	! or from Clegg and Brimblecombe , J. Chem. Eng., p43, 1995.
	825	;
	826	psh2o_out=18.452406985 -3505.1578807/T &
	827	& -330918.55082/(T*T) &
	828	& +12725068.262/(TTT)
	829	! in Pa
	830	psh2o_out=100.*exp(psh2o_out)
	831	end if
	832	! print*,psh2o_out
	833
	834	return
	835	end function psh2o
	836	!-----------------------------------------------------------------------
	837	real function density(T,so4mfrac) result(density_out)
	838	! calculation of particle density (gr/cm3)
	839
	840	! requires Temperature (T) and acid mass fraction (so4mfrac)
	841	!---->Vehkamaeki et al. (2002)
	842
	843	implicit none
	844	real, intent(in) :: T, so4mfrac
	845	real, parameter :: &
	846	& a1= 0.7681724,&
	847	& a2= 2.184714, &
	848	& a3= 7.163002, &
	849	& a4=-44.31447, &
	850	& a5= 88.74606, &
	851	& a6=-75.73729, &
	852	& a7= 23.43228
	853	real, parameter :: &
	854	& b1= 1.808225e-3, &
	855	& b2=-9.294656e-3, &
	856	& b3=-3.742148e-2, &
	857	& b4= 2.565321e-1, &
	858	& b5=-5.362872e-1, &
	859	& b6= 4.857736e-1, &
	860	& b7=-1.629592e-1
	861	real, parameter :: &
	862	& c1=-3.478524e-6, &
	863	& c2= 1.335867e-5, &
	864	& c3= 5.195706e-5, &
	865	& c4=-3.717636e-4, &
	866	& c5= 7.990811e-4, &
	867	& c6=-7.458060e-4, &
	868	& c7= 2.581390e-4
	869	real :: a,b,c,so4m2,so4m3,so4m4,so4m5,so4m6
	870
	871	so4m2=so4mfrac*so4mfrac
	872	so4m3=so4mfrac*so4m2
	873	so4m4=so4mfrac*so4m3
	874	so4m5=so4mfrac*so4m4
	875	so4m6=so4mfrac*so4m5
	876
	877	a=+a1+a2so4mfrac+a3so4m2+a4*so4m3 &
	878	& +a5so4m4+a6so4m5+a7*so4m6
	879	b=+b1+b2so4mfrac+b3so4m2+b4*so4m3 &
	880	& +b5so4m4+b6so4m5+b7*so4m6
	881	c=+c1+c2so4mfrac+c3so4m2+c4*so4m3 &
	882	& +c5so4m4+c6so4m5+c7*so4m6
	883	density_out=(a+bT+cTT) ! units are gm/cm*3
	884
	885	return
	886	end function density
	887	!-----------------------------------------------------------------------
	888	real function surftension(T,so4frac) result(surftension_out)
	889	! calculation of surface tension (mN/meter)
	890	! requires Temperature (T) and acid mole fraction (so4frac)
	891	!---->Vehkamaeki et al. (2002)
	892
	893	implicit none
	894	real,intent(in) :: T, so4frac
	895	real :: a,b,so4mfrac,so4m2,so4m3,so4m4,so4m5,so4sig
	896	real, parameter :: &
	897	& a1= 0.11864, &
	898	& a2=-0.11651, &
	899	& a3= 0.76852, &
	900	& a4=-2.40909, &
	901	& a5= 2.95434, &
	902	& a6=-1.25852
	903	real, parameter :: &
	904	& b1=-1.5709e-4, &
	905	& b2= 4.0102e-4, &
	906	& b3=-2.3995e-3, &
	907	& b4= 7.611235e-3, &
	908	& b5=-9.37386e-3, &
	909	& b6= 3.89722e-3
	910	real, parameter :: convfac=1.e3 ! convert from newton/m to dyne/cm
	911	real, parameter :: Mw=18.01528, Ma=98.079
	912
	913	! so4 mass fraction
	914	so4mfrac=Maso4frac/( Maso4frac+Mw*(1.-so4frac) )
	915	so4m2=so4mfrac*so4mfrac
	916	so4m3=so4mfrac*so4m2
	917	so4m4=so4mfrac*so4m3
	918	so4m5=so4mfrac*so4m4
	919
	920	a=+a1+a2so4mfrac+a3so4m2+a4so4m3+a5so4m4+a6*so4m5
	921	b=+b1+b2so4mfrac+b3so4m2+b4so4m3+b5so4m4+b6*so4m5
	922	so4sig=a+b*T
	923	surftension_out=so4sig*convfac
	924
	925	return
	926	end function surftension
	927	!-----------------------------------------------------------------------
	928	real function wph2so4(pph2o,T) result(wph2so4_out)
	929	! Calculates the equilibrium composition of h2so4 aerosols
	930	! as a function of temperature and H2O pressure, using
	931	! the parameterization of Tabazadeh et al., GRL, p1931, 1997.
	932	!
	933	! Parameters
	934	!
	935	! input:
	936	! T.....temperature (K)
	937	! pph2o..... amhbiant 2o pressure (Pa)
	938	!
	939	! output:
	940	! wph2so4......sulfuric acid composition (weight percent wt % h2so4)
	941	! = h2so4 mass fraction*100.
	942	!
	943	implicit none
	944	real, intent(in) :: pph2o, T
	945
	946	real :: aw, rh, y1, y2, sulfmolal
	947
	948	! psh2o(T): equilibrium H2O pressure over pure liquid water (Pa)
	949	! relative humidity
	950	rh=pph2o/psh2o(T)
	951	! water activity
	952	! aw=min( 0.999,max(1.e-3,rh) )
	953	aw=min( 0.999999999,max(1.e-8,rh) )
	954
	955	! composition
	956	! calculation of h2so4 molality
	957	if(aw .le. 0.05 .and. aw .gt. 0.) then
	958	y1=12.372089320aw*(-0.16125516114) &
	959	& -30.490657554*aw -2.1133114241
	960	y2=13.455394705aw*(-0.19213122550) &
	961	& -34.285174607*aw -1.7620073078
	962	else if(aw .le. 0.85 .and. aw .gt. 0.05) then
	963	y1=11.820654354aw*(-0.20786404244) &
	964	& -4.8073063730*aw -5.1727540348
	965	y2=12.891938068aw*(-0.23233847708) &
	966	& -6.4261237757*aw -4.9005471319
	967	else
	968	y1=-180.06541028aw*(-0.38601102592) &
	969	& -93.317846778*aw +273.88132245
	970	y2=-176.95814097aw*(-0.36257048154) &
	971	& -90.469744201*aw +267.45509988
	972	end if
	973	! h2so4 molality (m=moles of h2so4 (solute)/ kg of h2o(solvent))
	974	sulfmolal = y1+((T-190.)*(y2-y1)/70.)
	975
	976	! for a solution containing mh2so4 and mh2o:
	977	! sulfmolal = (mh2so4(gr)/h2so4_molar_mass(gr/mole)) / (mh2o(gr)*1.e-3)
	978	! mh2o=1.e3(mh2so4/Mh2so4)/sulfmolal=1.e3mh2so4/(Mh2so4*sulfmolal)
	979	! h2so4_mass_fraction = mfh2so4 = mh2so4/(mh2o + mh2so4)
	980	! mh2o=mh2so4*(1-mfh2so4)/mfh2so4
	981	! combining the 2 equations
	982	! 1.e3mh2so4/(Mh2so4sulfmolal) = mh2so4*(1-mfh2so4)/mfh2so4
	983	! 1.e3/(Mh2so4*sulfmolal) = (1-mfh2so4)/mfh2so4
	984	! 1000mfh2so4 = (1-mfh2so4)Mh2so4*sulfmolal
	985	! mfh2so4(1000.+Mh2so4sulfmolal) = Mh2so4*sulfmolal
	986	! mfh2so4 = Mh2so4sulfmolal / (1000.+Mh2so4sulfmolal)
	987	! wph2so4 (% mass fraction)= 100.Mh2so4sulfmolal / (1000.+Mh2so4*sulfmolal)
	988	! recall activity of i = a_i = P_i/P_pure_i and
	989	! activity coefficient of i = gamma_i = a_i/X_i (X_i: mole fraction of i)
	990	! so P_i = gamma_iX_iP_pure_i
	991	! if ideal solution, gamma_i=1, P_i = X_i*P_pure_i
	992
	993	! h2so4 weight precent
	994	wph2so4_out = 9800.sulfmolal/(98.sulfmolal+1000.)
	995	! print*,rh,pph2o,psh2o(T),vpice(T)
	996	! print*,T,aw,sulfmolal,wph2so4_out
	997	wph2so4_out = max(wph2so4_out,15.)
	998	wph2so4_out = min(wph2so4_out,99.999)
	999
	1000	return
	1001	end function wph2so4
	1002	!-----------------------------------------------------------------------
	1003	real function solh2so4(T,xa) result(solh2so4_out)
	1004	! equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
	1005
	1006	implicit none
	1007	real, intent(in) :: T, xa ! T(K) xa(H2SO4 mass fraction)
	1008
	1009	real :: xw, a12,b12, cacta, presat
	1010
	1011	xw=1.0-xa
	1012
	1013	! pure h2so4 saturation number density (molec/cm3)
	1014	presat=ndsh2so4(T)
	1015	! compute activity of acid
	1016	a12=5.672E3 -4.074E6/T +4.421E8/(T*T)
	1017	b12=1./0.527
	1018	cacta=10.*(a12xwxw/(xw+b12xa)**2/T)
	1019	! h2so4 saturation number density over H2SO4/H2O solution (molec/cm3)
	1020	solh2so4_out=cactaxapresat
	1021
	1022	return
	1023	end function solh2so4
	1024	!-----------------------------------------------------------------------
	1025	real function rpmvh2so4(T,ws) result(rpmvh2so4_out)
	1026	! partial molar volume of h2so4 in h2so4/h2o solution (cm3/mole)
	1027
	1028	implicit none
	1029	real, intent(in) :: T, ws
	1030	real, dimension(22),parameter :: x=(/ &
	1031	& 2.393284E-02,-4.359335E-05,7.961181E-08,0.0,-0.198716351, &
	1032	& 1.39564574E-03,-2.020633E-06,0.51684706,-3.0539E-03,4.505475E-06, &
	1033	& -0.30119511,1.840408E-03,-2.7221253742E-06,-0.11331674116, &
	1034	& 8.47763E-04,-1.22336185E-06,0.3455282,-2.2111E-03,3.503768245E-06, &
	1035	& -0.2315332,1.60074E-03,-2.5827835E-06/)
	1036
	1037	real :: w
	1038
	1039	w=ws*0.01
	1040	rpmvh2so4_out=x(5)+x(6)T+x(7)TT+(x(8)+x(9)T+x(10)TT)*w &
	1041	+(x(11)+x(12)T+x(13)TT)w*w
	1042	! h2so4 partial molar volume in h2so4/h2o solution (cm3/mole)
	1043	rpmvh2so4_out=rpmvh2so4_out*1000.
	1044
	1045	return
	1046	end function rpmvh2so4
	1047	!-----------------------------------------------------------------------
	1048	real function rmvh2o(T) result(rmvh2o_out)
	1049	! molar volume of pure h2o (cm3/mole)
	1050
	1051	implicit none
	1052	real, intent(in) :: T
	1053	real, parameter :: x1=2.393284E-02,x2=-4.359335E-05,x3=7.961181E-08
	1054
	1055	! 1000: L/mole -> cm3/mole
	1056	! pure h2o molar volume (cm3/mole)
	1057	rmvh2o_out=(x1+x2T+x3TT)1000.
	1058
	1059	return
	1060	end function rmvh2o
	1061	!
[2690]	1062	END MODULE sulfate_aer_mod

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

Download in other formats:

Original Format