Changeset 4950 for LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90
- Timestamp:
- May 22, 2024, 3:16:36 PM (4 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90
r4750 r4950 7 7 8 8 !******************************************************************* 9 SUBROUTINE STRACOMP_BIN(sh,t_seri,pplay) 10 ! 11 ! Aerosol H2SO4 weight fraction as a function of PH2O and temperature 12 ! INPUT: 13 ! sh: VMR of H2O 14 ! t_seri: temperature (K) 15 ! pplay: middle layer pression (Pa) 16 ! 17 ! OUTPUT: 18 ! R2SO4: aerosol H2SO4 weight fraction (percent) 9 SUBROUTINE STRACOMP_KELVIN(sh,t_seri,pplay) 10 ! 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) 16 ! 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) 19 26 20 USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands 21 USE aerophys 22 USE phys_local_var_mod, ONLY: R2SO4 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 23 34 24 IMPLICIT NONE35 IMPLICIT NONE 25 36 26 REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature 27 REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression in the middle of each layer (Pa) 28 REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! specific humidity 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 29 52 30 REAL ks(7) 31 REAL t,qh2o,ptot,pw 32 REAL a,b,c,det 33 REAL xsb,msb 53 DO ilon=1,klon 54 DO ilev=1,klev 55 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-16*temp)/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-2*RRSI(IK)*f_r_wetB(ilon,ilev,IK) 97 else 98 ! H2SO4-H2O flat surface, only on the first timestep 99 radwet = 1.e-2*RRSI(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-3*surtens*1.e-6*mvh2o/ (radwet*rgas*temp) ) 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-2*RRSI(IK)*(dens_aer_dry/(denso4ik*1.e3)/ & 123 & (r2so4ik*1.e-2))**third 124 fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2o / (radwet*rgas*temp) ) 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 34 148 35 INTEGER ilon,ilev 36 DATA ks/-21.661,2724.2,51.81,-15732.0,47.004,-6969.0,-4.6183/ 37 38 !******************************************************************* 39 !*** liquid aerosols process 40 !******************************************************************* 41 ! BINARIES LIQUID AEROROLS: 42 43 DO ilon=1,klon 44 DO ilev=1,klev 45 46 t = max(t_seri(ilon,ilev),185.) 47 qh2o=sh(ilon,ilev)/18.*28.9 48 ptot=pplay(ilon,ilev)/100. 49 pw = qh2o*ptot/1013.0 50 pw = min(pw,2.e-3/1013.) 51 pw = max(pw,2.e-5/1013.) 52 53 !******************************************************************* 54 !*** binaries aerosols h2so4/h2o 55 !******************************************************************* 56 a = ks(3) + ks(4)/t 57 b = ks(1) + ks(2)/t 58 c = ks(5) + ks(6)/t + ks(7)*log(t) - log(pw) 59 60 det = b**2 - 4.*a*c 61 62 IF (det > 0.) THEN 63 xsb = (-b - sqrt(det))/(2.*a) 64 msb = 55.51*xsb/(1.0 - xsb) 65 ELSE 66 msb = 0. 67 ENDIF 68 R2SO4(ilon,ilev) = 100*msb*0.098076/(1.0 + msb*0.098076) 69 70 ! H2SO4 min dilution: 0.5% 71 R2SO4(ilon,ilev) = max( R2SO4(ilon,ilev), 0.005 ) 72 ENDDO 73 ENDDO 74 100 RETURN 75 76 END SUBROUTINE STRACOMP_BIN 77 149 END SUBROUTINE STRACOMP_KELVIN 78 150 !******************************************************************** 79 151 SUBROUTINE STRACOMP(sh,t_seri,pplay) … … 544 616 545 617 END SUBROUTINE 546 547 !****************************************************************548 SUBROUTINE DENH2SA_TABA(t_seri)549 550 ! AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T551 ! from Tabazadeh et al. (1994) abaques552 ! ---------------------------------------------553 554 !555 ! INPUT:556 ! R2SO4: aerosol H2SO4 weight fraction (percent)557 ! t_seri: temperature (K)558 ! klon: number of latitude bands in the model domain559 ! klev: number of altitude bands in the model domain560 ! for IFS: perhaps add another dimension for longitude561 !562 ! OUTPUT:563 ! DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)564 !565 USE dimphy, ONLY : klon,klev566 USE phys_local_var_mod, ONLY: R2SO4, DENSO4567 568 IMPLICIT NONE569 570 REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature571 572 INTEGER i,j573 574 !----------------------------------------------------------------------575 ! ... Local variables576 !----------------------------------------------------------------------577 real, parameter :: a9 = -268.2616e4, a10 = 576.4288e3578 579 real :: a0, a1, a2, a3, a4, a5, a6, a7 ,a8580 real :: c1, c2, c3, c4, w581 582 583 ! Loop on model domain (2 dimension for UPMC model; 3 for IFS)584 DO i=1,klon585 DO j=1,klev586 !----------------------------------------------------------------------587 ! ... Temperature variables588 !----------------------------------------------------------------------589 c1 = t_seri(I,J)- 273.15590 c2 = c1**2591 c3 = c1*c2592 c4 = c1*c3593 !----------------------------------------------------------------------594 ! Polynomial Coefficients595 !----------------------------------------------------------------------596 a0 = 999.8426 + 334.5402e-4*c1 - 569.1304e-5*c2597 a1 = 547.2659 - 530.0445e-2*c1 + 118.7671e-4*c2 + 599.0008e-6*c3598 a2 = 526.295e1 + 372.0445e-1*c1 + 120.1909e-3*c2 - 414.8594e-5*c3 + 119.7973e-7*c4599 a3 = -621.3958e2 - 287.7670*c1 - 406.4638e-3*c2 + 111.9488e-4*c3 + 360.7768e-7*c4600 a4 = 409.0293e3 + 127.0854e1*c1 + 326.9710e-3*c2 - 137.7435e-4*c3 - 263.3585e-7*c4601 a5 = -159.6989e4 - 306.2836e1*c1 + 136.6499e-3*c2 + 637.3031e-5*c3602 a6 = 385.7411e4 + 408.3717e1*c1 - 192.7785e-3*c2603 a7 = -580.8064e4 - 284.4401e1*c1604 a8 = 530.1976e4 + 809.1053*c1605 !----------------------------------------------------------------------606 ! ... Summation607 !----------------------------------------------------------------------608 ! w : H2SO4 Weight fraction609 w=r2SO4(i,j)*0.01610 DENSO4(i,j) = 0.001*(a0 + w*(a1 + w*(a2 + w*(a3 + w*(a4 + &611 w*(a5 + w*(a6 + w*(a7 + w*(a8 + w*(a9 + w*a10))))))))))612 DENSO4(i,j) = max (0.0, DENSO4(i,j) )613 614 ENDDO615 ENDDO616 617 END SUBROUTINE DENH2SA_TABA618 618 619 619 !**************************************************************** … … 764 764 RETURN 765 765 END SUBROUTINE 766 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) 774 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=b1*exp( -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=presat*10./(akb*T) 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.2*log(T) + 5.7133e-3*T ) 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/(T*T*T) 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+a2*so4mfrac+a3*so4m2+a4*so4m3 & 878 & +a5*so4m4+a6*so4m5+a7*so4m6 879 b=+b1+b2*so4mfrac+b3*so4m2+b4*so4m3 & 880 & +b5*so4m4+b6*so4m5+b7*so4m6 881 c=+c1+c2*so4mfrac+c3*so4m2+c4*so4m3 & 882 & +c5*so4m4+c6*so4m5+c7*so4m6 883 density_out=(a+b*T+c*T*T) ! 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=Ma*so4frac/( Ma*so4frac+Mw*(1.-so4frac) ) 915 so4m2=so4mfrac*so4mfrac 916 so4m3=so4mfrac*so4m2 917 so4m4=so4mfrac*so4m3 918 so4m5=so4mfrac*so4m4 919 920 a=+a1+a2*so4mfrac+a3*so4m2+a4*so4m3+a5*so4m4+a6*so4m5 921 b=+b1+b2*so4mfrac+b3*so4m2+b4*so4m3+b5*so4m4+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.372089320*aw**(-0.16125516114) & 959 & -30.490657554*aw -2.1133114241 960 y2=13.455394705*aw**(-0.19213122550) & 961 & -34.285174607*aw -1.7620073078 962 else if(aw .le. 0.85 .and. aw .gt. 0.05) then 963 y1=11.820654354*aw**(-0.20786404244) & 964 & -4.8073063730*aw -5.1727540348 965 y2=12.891938068*aw**(-0.23233847708) & 966 & -6.4261237757*aw -4.9005471319 967 else 968 y1=-180.06541028*aw**(-0.38601102592) & 969 & -93.317846778*aw +273.88132245 970 y2=-176.95814097*aw**(-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.e3*mh2so4/(Mh2so4*sulfmolal) 979 ! h2so4_mass_fraction = mfh2so4 = mh2so4/(mh2o + mh2so4) 980 ! mh2o=mh2so4*(1-mfh2so4)/mfh2so4 981 ! combining the 2 equations 982 ! 1.e3*mh2so4/(Mh2so4*sulfmolal) = mh2so4*(1-mfh2so4)/mfh2so4 983 ! 1.e3/(Mh2so4*sulfmolal) = (1-mfh2so4)/mfh2so4 984 ! 1000*mfh2so4 = (1-mfh2so4)*Mh2so4*sulfmolal 985 ! mfh2so4*(1000.+Mh2so4*sulfmolal) = Mh2so4*sulfmolal 986 ! mfh2so4 = Mh2so4*sulfmolal / (1000.+Mh2so4*sulfmolal) 987 ! wph2so4 (% mass fraction)= 100.*Mh2so4*sulfmolal / (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_i*X_i*P_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.**(a12*xw*xw/(xw+b12*xa)**2/T) 1019 ! h2so4 saturation number density over H2SO4/H2O solution (molec/cm3) 1020 solh2so4_out=cacta*xa*presat 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)*T*T+(x(8)+x(9)*T+x(10)*T*T)*w & 1041 +(x(11)+x(12)*T+x(13)*T*T)*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+x2*T+x3*T*T)*1000. 1058 1059 return 1060 end function rmvh2o 1061 ! 767 1062 END MODULE sulfate_aer_mod
Note: See TracChangeset
for help on using the changeset viewer.