Ignore:
Timestamp:
Sep 20, 2024, 12:32:04 PM (5 days ago)
Author:
Laurent Fairhead
Message:

Updating cirrus branch to trunk revision 5171

Location:
LMDZ6/branches/cirrus
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/cirrus

  • LMDZ6/branches/cirrus/libf/phylmd/StratAer/sulfate_aer_mod.F90

    r4750 r5202  
    77
    88!*******************************************************************
    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)
    1926   
    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
    2334   
    24     IMPLICIT NONE
     35      IMPLICIT NONE
    2536   
    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
    2952   
    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
    34148   
    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
    78150!********************************************************************
    79151    SUBROUTINE STRACOMP(sh,t_seri,pplay)
    … …  
    544616
    545617    END SUBROUTINE
    546 
    547 !****************************************************************
    548     SUBROUTINE DENH2SA_TABA(t_seri)
    549 
    550 !   AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T
    551 !   from Tabazadeh et al. (1994) abaques
    552 !   ---------------------------------------------
    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 domain
    559 !   klev: number of altitude bands in the model domain
    560 !   for IFS: perhaps add another dimension for longitude
    561 !
    562 !   OUTPUT:
    563 !   DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)
    564 !   
    565     USE dimphy, ONLY : klon,klev
    566     USE phys_local_var_mod, ONLY: R2SO4, DENSO4
    567    
    568     IMPLICIT NONE
    569    
    570     REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
    571        
    572     INTEGER i,j
    573    
    574 !----------------------------------------------------------------------
    575 !       ... Local variables
    576 !----------------------------------------------------------------------
    577       real, parameter :: a9 = -268.2616e4, a10 = 576.4288e3
    578      
    579       real :: a0, a1, a2, a3, a4, a5, a6, a7 ,a8
    580       real :: c1, c2, c3, c4, w
    581      
    582      
    583 !   Loop on model domain (2 dimension for UPMC model; 3 for IFS)
    584     DO i=1,klon
    585        DO j=1,klev
    586 !----------------------------------------------------------------------
    587 !       ... Temperature variables
    588 !----------------------------------------------------------------------
    589           c1 = t_seri(I,J)- 273.15
    590           c2 = c1**2
    591           c3 = c1*c2
    592           c4 = c1*c3
    593 !----------------------------------------------------------------------
    594 !       Polynomial Coefficients
    595 !----------------------------------------------------------------------
    596           a0 = 999.8426 + 334.5402e-4*c1 - 569.1304e-5*c2
    597           a1 = 547.2659 - 530.0445e-2*c1 + 118.7671e-4*c2 + 599.0008e-6*c3
    598           a2 = 526.295e1 + 372.0445e-1*c1 + 120.1909e-3*c2 - 414.8594e-5*c3 + 119.7973e-7*c4
    599           a3 = -621.3958e2 - 287.7670*c1 - 406.4638e-3*c2 + 111.9488e-4*c3 + 360.7768e-7*c4
    600           a4 = 409.0293e3 + 127.0854e1*c1 + 326.9710e-3*c2 - 137.7435e-4*c3 - 263.3585e-7*c4
    601           a5 = -159.6989e4 - 306.2836e1*c1 + 136.6499e-3*c2 + 637.3031e-5*c3
    602           a6 = 385.7411e4 + 408.3717e1*c1 - 192.7785e-3*c2
    603           a7 = -580.8064e4 - 284.4401e1*c1
    604           a8 = 530.1976e4 + 809.1053*c1
    605 !----------------------------------------------------------------------
    606 !       ... Summation
    607 !----------------------------------------------------------------------
    608 !     w : H2SO4 Weight fraction
    609           w=r2SO4(i,j)*0.01
    610           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        ENDDO
    615     ENDDO
    616 
    617   END SUBROUTINE DENH2SA_TABA
    618618 
    619619!****************************************************************
    … …  
    764764       RETURN
    765765       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!
    7671062END MODULE sulfate_aer_mod
Note: See TracChangeset for help on using the changeset viewer.