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

Updating cirrus branch to trunk revision 5171

Location:
LMDZ6/branches/cirrus
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/cirrus

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

    r3677 r5202  
    1717!-----------------------------------------------------------------------
    1818
    19   USE phys_local_var_mod, ONLY: mdw, budg_sed_part, DENSO4, f_r_wet, vsed_aer
     19  USE phys_local_var_mod, ONLY: mdw, budg_sed_part, DENSO4, DENSO4B, f_r_wet, f_r_wetB, vsed_aer
     20  USE strataer_local_var_mod, ONLY: flag_new_strat_compo
    2021  USE dimphy, ONLY : klon,klev
    2122  USE infotrac_phy
    … …  
    8990
    9091      ! stokes-velocity with cunnigham slip- flow correction
    91       ZVAER(JL,JK,nb) = 2./9.*(DENSO4(JL,JK)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wet(JL,JK)*mdw(nb)/2.)**2.* &
    92          (1.+ 2.*zlair(JL,JK)/(f_r_wet(JL,JK)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wet(JL,JK)*mdw(nb)/zlair(JL,JK))))
    93 
     92      IF(flag_new_strat_compo) THEN
     93         ! stokes-velocity with cunnigham slip- flow correction
     94         ZVAER(JL,JK,nb) = 2./9.*(DENSO4B(JL,JK,nb)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wetB(JL,JK,nb)*mdw(nb)/2.)**2.* &
     95              (1.+ 2.*zlair(JL,JK)/(f_r_wetB(JL,JK,nb)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wetB(JL,JK,nb)*mdw(nb)/zlair(JL,JK))))
     96      ELSE
     97         ZVAER(JL,JK,nb) = 2./9.*(DENSO4(JL,JK)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wet(JL,JK)*mdw(nb)/2.)**2.* &
     98              (1.+ 2.*zlair(JL,JK)/(f_r_wet(JL,JK)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wet(JL,JK)*mdw(nb)/zlair(JL,JK))))
     99      ENDIF
     100     
    94101      ZSEDFLX(JL,nb)=ZVAER(JL,JK,nb)*ZRHO
    95102      ZSOLAERB(nb)=ZSOLAERB(nb)+ZDTGDP*ZSEDFLX(JL,nb)

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

    r4601 r5202  
    55  IMPLICIT NONE
    66!
    7   REAL,PARAMETER                         :: ropx=1500.0              ! default aerosol particle mass density [kg/m3]
    8   REAL,PARAMETER                         :: dens_aer_dry=1848.682308 ! dry aerosol particle mass density at T_0=293K[kg/m3]
    9   REAL,PARAMETER                         :: dens_aer_ref=1483.905336 ! aerosol particle mass density with 75% H2SO4 at T_0=293K[kg/m3]
    10   REAL,PARAMETER                         :: mdwmin=0.002e-6          ! dry diameter of smallest aerosol particles [m]
    11   REAL,PARAMETER                         :: V_rat=2.0                ! volume ratio of neighboring size bins
    12   REAL,PARAMETER                         :: mfrac_H2SO4=0.75         ! default mass fraction of H2SO4 in the aerosol
    13   REAL, PARAMETER                        :: mAIRmol=28.949*1.66E-27  ! Average mass of an air molecule [kg]
    14   REAL, PARAMETER                        :: mH2Omol=18.016*1.66E-27  ! Mass of an H2O molecule [kg]
    15   REAL, PARAMETER                        :: mH2SO4mol=98.082*1.66E-27! Mass of an H2SO4 molecule [kg]
    16   REAL, PARAMETER                        :: mSO2mol=64.06*1.66E-27   ! Mass of an SO2 molecule [kg]
    17   REAL, PARAMETER                        :: mSatom=32.06*1.66E-27    ! Mass of a S atom [kg]
    18   REAL, PARAMETER                        :: mOCSmol=60.07*1.66E-27   ! Mass of an OCS molecule [kg]
    19   REAL, PARAMETER                        :: mClatom=35.45*1.66E-27   ! Mass of an Cl atom [kg]
    20   REAL, PARAMETER                        :: mHClmol=36.46*1.66E-27   ! Mass of an HCl molecule [kg]
    21   REAL, PARAMETER                        :: mBratom=79.90*1.66E-27   ! Mass of an Br atom [kg]
    22   REAL, PARAMETER                        :: mHBrmol=80.92*1.66E-27   ! Mass of an HBr molecule [kg]
    23   REAL, PARAMETER                        :: mNOmol=30.01*1.66E-27    ! Mass of an NO molecule [kg]
    24   REAL, PARAMETER                        :: mNO2mol=46.01*1.66E-27   ! Mass of an NO2 molecule [kg]
    25   REAL, PARAMETER                        :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg]
     7  REAL,PARAMETER    :: ropx=1500.0              ! default aerosol particle mass density [kg/m3]
     8  REAL,PARAMETER    :: dens_aer_dry=1848.682308 ! dry aerosol particle mass density at T_0=293K[kg/m3]
     9  REAL,PARAMETER    :: dens_aer_ref=1483.905336 ! aerosol particle mass density with 75% H2SO4 at T_0=293K[kg/m3]
     10  REAL,PARAMETER    :: mdwmin=0.002e-6          ! dry diameter of smallest aerosol particles [m]
     11  REAL,PARAMETER    :: V_rat=2.0                ! volume ratio of neighboring size bins
     12  REAL,PARAMETER    :: mfrac_H2SO4=0.75         ! default mass fraction of H2SO4 in the aerosol
     13  REAL, PARAMETER   :: mAIRmol=28.949*1.66E-27  ! Average mass of an air molecule [kg]
     14  REAL, PARAMETER   :: mH2Omol=18.016*1.66E-27  ! Mass of an H2O molecule [kg]
     15  REAL, PARAMETER   :: mH2SO4mol=98.082*1.66E-27! Mass of an H2SO4 molecule [kg]
     16  REAL, PARAMETER   :: mSO2mol=64.06*1.66E-27   ! Mass of an SO2 molecule [kg]
     17  REAL, PARAMETER   :: mSatom=32.06*1.66E-27    ! Mass of a S atom [kg]
     18  REAL, PARAMETER   :: mOCSmol=60.07*1.66E-27   ! Mass of an OCS molecule [kg]
     19  REAL, PARAMETER   :: mClatom=35.45*1.66E-27   ! Mass of an Cl atom [kg]
     20  REAL, PARAMETER   :: mHClmol=36.46*1.66E-27   ! Mass of an HCl molecule [kg]
     21  REAL, PARAMETER   :: mBratom=79.90*1.66E-27   ! Mass of an Br atom [kg]
     22  REAL, PARAMETER   :: mHBrmol=80.92*1.66E-27   ! Mass of an HBr molecule [kg]
     23  REAL, PARAMETER   :: mNOmol=30.01*1.66E-27    ! Mass of an NO molecule [kg]
     24  REAL, PARAMETER   :: mNO2mol=46.01*1.66E-27   ! Mass of an NO2 molecule [kg]
     25  REAL, PARAMETER   :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg]
     26  REAL, PARAMETER   :: rgas=8.3145 ! molar gas cste (J⋅K−1⋅mol−1=m3⋅Pa⋅K−1⋅mol−1=kg⋅m2⋅s−2⋅K−1⋅mol−1)
     27  !
     28  REAL, PARAMETER   :: MH2O  =1000.*mH2Omol     ! Mass of 1 molec [g] (18.016*1.66E-24)
     29  REAL, PARAMETER   :: MH2SO4=1000.*mH2SO4mol   ! Mass of 1 molec [g] (98.082*1.66E-24)
     30  REAL, PARAMETER   :: BOLZ  =1.381E-16         ! Boltzmann constant [dyn.cm/K]
    2631!
    2732END MODULE aerophys

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

    r4762 r5202  
    2626  USE aerophys
    2727  USE infotrac_phy
    28   USE phys_local_var_mod, ONLY: DENSO4, f_r_wet
    29 
     28  USE phys_local_var_mod, ONLY: DENSO4, DENSO4B, f_r_wet, f_r_wetB
     29  USE strataer_local_var_mod, ONLY: flag_new_strat_compo
     30 
    3031  IMPLICIT NONE
    3132
    … …  
    4344  ! local variables in coagulation routine
    4445  INTEGER                                       :: i,j,k,nb,ilon,ilev
    45   REAL, DIMENSION(nbtr_bin)                     :: radius ! aerosol particle radius in each bin [m]
     46  REAL, DIMENSION(nbtr_bin)                     :: radiusdry ! dry aerosol particle radius in each bin [m]
     47  REAL, DIMENSION(nbtr_bin)                     :: radiuswet ! wet aerosol particle radius in each bin [m]
    4648  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_t ! Concentration Traceur at time t [U/KgA]
    4749  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
    4850  REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin)   :: ff   ! Volume fraction of intermediate particles
    49   REAL, DIMENSION(nbtr_bin)                     :: V    ! Volume of bins
     51  REAL, DIMENSION(nbtr_bin)                     :: Vdry ! Volume dry of bins
     52  REAL, DIMENSION(nbtr_bin)                     :: Vwet ! Volume wet of bins
    5053  REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: Vij  ! Volume sum of i and j
    5154  REAL                                          :: eta  ! Dynamic viscosity of air
    … …  
    8285  include "YOMCST.h"
    8386
    84   DO i=1, nbtr_bin
    85    radius(i)=mdw(i)/2.
    86    V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
    87   ENDDO
    88 
    89   DO j=1, nbtr_bin
    90   DO i=1, nbtr_bin
    91    Vij(i,j)= V(i)+V(j)
    92   ENDDO
    93   ENDDO
    94 
     87! ff(i,j,k): Volume fraction of Vi,j that is partitioned to each model bin k
     88! just need to be calculated in model initialization because mdw(:) size is fixed
     89! no need to recalculate radius, Vdry, Vij, and ff every timestep because it is for 
     90! dry aerosols
     91  DO i=1, nbtr_bin
     92     radiusdry(i)=mdw(i)/2.
     93     Vdry(i)=radiusdry(i)**3.  !neglecting factor 4*RPI/3
     94     Vwet(i)=0.0
     95  ENDDO
     96
     97  DO j=1, nbtr_bin
     98     DO i=1, nbtr_bin
     99        Vij(i,j)= Vdry(i)+Vdry(j)
     100     ENDDO
     101  ENDDO
     102 
    95103!--pre-compute the f(i,j,k) from Jacobson equation 13
    96104  ff=0.0
    … …  
    100108    IF (k.EQ.1) THEN
    101109      ff(i,j,k)= 0.0
    102     ELSEIF (k.GT.1.AND.V(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN
     110    ELSEIF (k.GT.1.AND.Vdry(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.Vdry(k)) THEN
    103111      ff(i,j,k)= 1.-ff(i,j,k-1)
    104112    ELSEIF (k.EQ.nbtr_bin) THEN
    105       IF (Vij(i,j).GE.v(k)) THEN
     113      IF (Vij(i,j).GE.Vdry(k)) THEN
    106114        ff(i,j,k)= 1.
    107115      ELSE
    108116        ff(i,j,k)= 0.0
    109117      ENDIF
    110     ELSEIF (k.LE.(nbtr_bin-1).AND.V(k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN
    111       ff(i,j,k)= V(k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))
     118    ELSEIF (k.LE.(nbtr_bin-1).AND.Vdry(k).LE.Vij(i,j).AND.Vij(i,j).LT.Vdry(k+1)) THEN
     119      ff(i,j,k)= Vdry(k)/Vij(i,j)*(Vdry(k+1)-Vij(i,j))/(Vdry(k+1)-Vdry(k))
    112120    ENDIF
    113121  ENDDO
    114122  ENDDO
    115123  ENDDO
    116 
     124! End of just need to be calculated at initialization because mdw(:) size is fixed
     125 
    117126  DO ilon=1, klon
    118127  DO ilev=1, klev
    … …  
    120129  IF (is_strato(ilon,ilev)) THEN
    121130  !compute actual wet particle radius & volume for every grid box
    122   DO i=1, nbtr_bin
    123    radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
    124    V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
    125   ENDDO
    126 
     131  IF(flag_new_strat_compo) THEN
     132     DO i=1, nbtr_bin
     133        radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2.
     134        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
     135!!      Vwet(i)= Vdry(i)*(f_r_wetB(ilon,ilev,i)**3)
     136     ENDDO
     137  ELSE
     138     DO i=1, nbtr_bin
     139        radiuswet(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
     140        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
     141!!      Vwet(i)= Vdry(i)*(f_r_wet(ilon,ilev)**3)
     142     ENDDO
     143  ENDIF
     144 
    127145!--Calculations for the coagulation kernel---------------------------------------------------------
    128146
    … …  
    150168  Di=0.0
    151169  DO i=1, nbtr_bin
    152    Kn(i)=mnfrpth/radius(i)
    153    Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radius(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
     170      Kn(i)=mnfrpth/radiuswet(i)
     171      Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radiuswet(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
    154172  ENDDO
    155173
    156174!--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
    157175  thvelpar=0.0
    158   DO i=1, nbtr_bin
    159    m_par(i)=4./3.*RPI*radius(i)**3.*DENSO4(ilon,ilev)*1000.
    160    thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
    161   ENDDO
     176  IF(flag_new_strat_compo) THEN
     177     DO i=1, nbtr_bin
     178        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4B(ilon,ilev,i)*1000.
     179        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
     180     ENDDO
     181  ELSE
     182     DO i=1, nbtr_bin
     183        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4(ilon,ilev)*1000.
     184        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
     185     ENDDO
     186  ENDIF
    162187
    163188!--pre-compute the particle mean free path mfppar(i) from equation 22
    … …  
    171196  delta=0.0
    172197  DO i=1, nbtr_bin
    173    delta(i)=((2.*radius(i)+mfppar(i))**3.-(4.*radius(i)**2.+mfppar(i)**2.)**1.5)/ &
    174            & (6.*radius(i)*mfppar(i))-2.*radius(i)
    175   ENDDO
    176 
     198      delta(i)=((2.*radiuswet(i)+mfppar(i))**3.-(4.*radiuswet(i)**2.+mfppar(i)**2.)**1.5)/ &
     199           & (6.*radiuswet(i)*mfppar(i))-2.*radiuswet(i)
     200  ENDDO
     201
     202!   beta(i,j): coagulation kernel (rate coefficient) of 2 colliding particles i,j
    177203!--pre-compute the beta(i,j) from equation 17 in Jacobson
    178204  num=0.0
    … …  
    180206  DO i=1, nbtr_bin
    181207!
    182    num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))
    183    denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
    184         & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))
    185    beta(i,j)=num/denom
     208     num=4.*RPI*(radiuswet(i)+radiuswet(j))*(Di(i)+Di(j))
     209     denom=(radiuswet(i)+radiuswet(j))/(radiuswet(i)+radiuswet(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
     210          & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radiuswet(i)+radiuswet(j)))
     211     beta(i,j)=num/denom
    186212!
    187213!--compute enhancement factor due to van der Waals forces
    188214   IF (ok_vdw .EQ. 0) THEN      !--no enhancement factor
    189      Evdw=1.0
     215      Evdw=1.0
    190216   ELSEIF (ok_vdw .EQ. 1) THEN  !--E(0) case
    191      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
    192      xvdW = LOG(1.+AvdWi)
    193      EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
     217      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
     218      xvdW = LOG(1.+AvdWi)
     219      EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
    194220   ELSEIF (ok_vdw .EQ. 2) THEN  !--E(infinity) case
    195      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
    196      xvdW = LOG(1.+AvdWi)
    197      EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
     221      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
     222      xvdW = LOG(1.+AvdWi)
     223      EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
    198224   ENDIF
    199225!
    … …  
    209235  denom=0.0
    210236  DO j=1, nbtr_bin
    211   denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
     237     !    fraction of coagulation of k and j that is not giving k
     238     denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
    212239  ENDDO
    213240
    … …  
    219246    num=0.0
    220247    DO j=1, k
    221     numi=0.0
    222     DO i=1, k-1
    223     numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
     248       numi=0.0
     249       DO i=1, k-1
     250!           
     251!           see Jacobson: " In order to conserve volume and volume concentration (which
     252!           coagulation physically does) while giving up some accuracy in number concentration"
     253!
     254!           Coagulation of i and j giving k
     255!           with V(i) and then V(j) because it considers i,j and j,i with the double loop
     256!
     257!           BUT WHY WET VOLUME V(i) in old STRATAER? tracers are already dry aerosols and coagulation
     258!           kernel beta(i,j) accounts for wet aerosols -> reply below
     259!
     260!             numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
     261            numi=numi+ff(i,j,k)*beta(i,j)*Vdry(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
     262       ENDDO
     263       num=num+numi
    224264    ENDDO
    225     num=num+numi
    226     ENDDO
    227265
    228266!--calculate new concentration of other bins
    229     tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/(1.+pdtcoag*denom)/V(k)
     267!      tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
     268    tr_tp1(ilon,ilev,k)=(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*Vdry(k) )
     269!
     270!       In constant composition (no dependency on aerosol size because no kelvin effect)
     271!       V(l)= (f_r_wet(ilon,ilev)**3)*((mdw(l)/2.)**3) = (f_r_wet(ilon,ilev)**3)*Vdry(i)
     272!       so numi and num are proportional (f_r_wet(ilon,ilev)**3)
     273!       and so
     274!        tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
     275!                     =(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num_dry)/( (1.+pdtcoag*denom)*Vdry(k) )
     276!          with num_dry=...beta(i,j)*Vdry(i)*....
     277!       so in old STRATAER (.not.flag_new_strat_compo), it was correct
    230278  ENDIF
    231279
    … …  
    234282!--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
    235283  DO i=1, nbtr_bin
    236    tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
     284     tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
    237285  ENDDO
    238286
    … …  
    240288  ENDDO !--end of loop klev
    241289  ENDDO !--end of loop klon
     290! *********************************************
    242291
    243292END SUBROUTINE COAGULATE

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

    r3677 r5202  
    99CONTAINS
    1010
     11      SUBROUTINE condens_evapor_rate_kelvin(R2SO4G,t_seri,pplay,R2SO4, &
     12          & DENSO4,f_r_wet,R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
     13!
     14!     INPUT:
     15!     R2SO4G: number density of gaseous H2SO4 [molecules/cm3]
     16!     t_seri: temperature (K)
     17!     pplay: pressure (Pa)
     18!     R2SO4: aerosol H2SO4 weight fraction (percent) - flat surface (does not depend on aerosol size)
     19!     DENSO4: aerosol density (gr/cm3)
     20!     f_r_wet: factor for converting dry to wet radius
     21!        assuming 'flat surface' composition (does not depend on aerosol size)
     22!     variables that depends on aerosol size because of Kelvin effect
     23!     R2SO4Gik: number density of gaseous H2SO4 [molecules/cm3] - depends on aerosol size
     24!     DENSO4ik: aerosol density (gr/cm3) - depends on aerosol size
     25!     f_r_wetik: factor for converting dry to wet radius - depends on aerosol size
     26!     RRSI: radius [cm]
     27
     28      USE aerophys
     29      USE infotrac_phy
     30      USE YOMCST, ONLY : RPI
     31      USE sulfate_aer_mod, ONLY : wph2so4, surftension, solh2so4, rpmvh2so4
     32      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI
     33     
     34      IMPLICIT NONE
     35     
     36      REAL, PARAMETER :: third=1./3.
     37     
     38      ! input variables
     39      REAL            :: R2SO4G !H2SO4 number density [molecules/cm3]
     40      REAL            :: t_seri
     41      REAL            :: pplay
     42      REAL            :: R2SO4
     43      REAL            :: DENSO4
     44      REAL            :: f_r_wet
     45      REAL            :: R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin)
     46     
     47      ! output variables
     48      REAL            :: FL(nbtr_bin)
     49      REAL            :: ASO4(nbtr_bin)
     50      REAL            :: DNDR(nbtr_bin)
     51     
     52      ! local variables
     53      INTEGER            :: IK
     54      REAL            :: ALPHA,CST
     55      REAL            :: WH2(nbtr_bin)
     56      REAL            :: RP,VTK,AA,FL1,RKNUD
     57      REAL            :: DND
     58      REAL            :: ATOT,AH2O
     59      REAL            :: RRSI_wet(nbtr_bin)
     60      REAL            :: FPATH, WPP, XA, FKELVIN
     61      REAL            :: surtens, mvh2so4, temp
     62     
     63! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
     64!    ------------------------------------------------------------------
     65!                                  EXCEPT CN
     66!       RK:H2SO4 WEIGHT PERCENT DOESN'T CHANGE
     67!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE
     68
     69!                   MOLECULAR ACCOMODATION OF H2SO4
     70!     H2SO4 accommodation  coefficient [condensation/evaporation]
     71      ALPHA = ALPH2SO4
     72!      FPLAIR=(2.281238E-5)*TAIR/PAIR
     73!     1.E2 (m to cm),
     74      CST=1.E2*2.281238E-5
     75!     same expression as in coagulate
     76!     in coagulate: mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
     77!     mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15)
     78!     mnfrpth=2.28E-5*t_seri/pplay
     79     
     80      temp = min( max(t_seri, 190.), 300.) ! 190K <= temp <= 300K
     81     
     82      RRSI_wet(:)=RRSI(:)*f_r_wetik(:)
     83
     84!     Pruppa and Klett
     85      FPATH=CST*t_seri/pplay
     86   
     87!     H2SO4 mass fraction in aerosol
     88      WH2(:)=R2SO4ik(:)*1.0E-2
     89
     90!                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
     91!                               AYERS ET AL (1980)
     92!                                  (MU-MU0)
     93!      RP=-10156.0/t_seri +16.259-(ACTSO4*4.184)/(8.31441*t_seri)
     94!                                  DROPLET H2SO4 PRESSURE IN DYN.CM-2
     95!      RP=EXP(RP)*1.01325E6/0.086
     96!!      RP=EXP(RP)*1.01325E6
     97!                                  H2SO4 NUMBER DENSITY NEAR DROPLET
     98
     99!      DND=RP*6.02E23/(8.31E7*t_seri)
     100
     101!                                 KELVIN EFFECT FACTOR
     102!CK 20160613: bug fix, removed factor 250 (from original code by S. Bekki)
     103!!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri*250.0)
     104!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri)
     105
     106!                                  MEAN KINETIC VELOCITY
     107!     DYN*CM*K/(K*GR)=(CM/SEC2)*CM
     108!                                  IN CM/SEC
     109      VTK=SQRT(8.0*BOLZ*t_seri/(RPI*MH2SO4))
     110!                                 KELVIN EFFECT FACTOR
     111
     112!     Loop on bin radius (RRSI in cm)
     113      DO IK=1,nbtr_bin
     114
     115      IF(R2SO4ik(IK) > 0.0) THEN
     116
     117!       h2so4 mass fraction (0<wpp<1)
     118        wpp=R2SO4ik(IK)*1.e-2   
     119        xa=18.*wpp/(18.*wpp+98.*(1.-wpp))
     120!       equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
     121        DND=solh2so4(t_seri,xa)
     122!          KELVIN EFFECT: 
     123!       surface tension (mN/m=1.e-3.kg/s2) = f(T,h2so4 mole fraction)
     124        surtens=surftension(temp,xa)
     125!       partial molar volume of h2so4 (cm3.mol-1 =1.e-6.m3.mol-1)
     126        mvh2so4= rpmvh2so4(temp,R2SO4ik(IK))
     127!       Kelvin factor (MKS)
     128        fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2so4/ (1.e-2*RRSI_wet(IK)*rgas*temp) )
     129!                             
     130        DNDR(IK) =DND*fkelvin
     131
     132        FL1=RPI*ALPHA*VTK*(R2SO4G-DNDR(IK))
     133
     134!       TURCO(1979) FOR HNO3:ALH2SO4 CONDENSATION= ALH2SO4 EVAPORATION
     135!       RPI*R2*VTK IS EQUIVALENT TO DIFFUSION COEFFICIENT
     136!       EXTENSION OF THE RELATION FOR DIFFUSION KINETICS
     137!       KNUDSEN NUMBER FPATH/RRSI
     138!       NEW VERSION (SEE NOTES)
     139        RKNUD=FPATH/RRSI_wet(IK)
     140!       SENFELD
     141        FL(IK)=FL1*RRSI_wet(IK)**2*( 1.0 +RKNUD ) &
     142     &     /( 1.0 +ALPHA/(2.0*RKNUD) +RKNUD )
     143!       TURCO
     144!        RL= (4.0/3.0 +0.71/RKNUD)/(1.0+1.0/RKNUD)
     145!     *         +4.0*(1.0-ALPHA)/(3.0*ALPHA)
     146!        FL=FL1*RRSI(IK)*RRSI(IK)
     147!     *         /( (3.0*ALPHA/4.0)*(1.0/RKNUD+RL*ALPHA) )
     148
     149!                         INITIAL NUMBER OF H2SO4 MOLEC OF 1 DROPLET
     150        ATOT=4.0*RPI*DENSO4ik(IK)*(RRSI_wet(IK)**3)/3.0 !attention: g and cm
     151        ASO4(IK)=WH2(IK)*ATOT/MH2SO4 !attention: g
     152!        ATOT=4.0*RPI*dens_aer(I,J)/1000.*(RRSI(IK)**3)/3.0
     153!        ASO4=mfrac_H2SO4*ATOT/MH2SO4
     154!                        INITIAL NUMBER OF H2O MOLEC OF 1 DROPLET
     155        AH2O=(1.0-WH2(IK))*ATOT/MH2O !attention: g
     156
     157!       CHANGE OF THE NUMBER OF H2SO4 MOLEC OF 1 DROPLET DURING DT
     158!       IT IS FOR KEM BUT THERE ARE OTHER WAYS
     159
     160      ENDIF 
     161
     162      ENDDO !loop over bins
     163
     164      END SUBROUTINE condens_evapor_rate_kelvin
     165     
     166!********************************************************************
    11167      SUBROUTINE condens_evapor_rate(R2SO4G,t_seri,pplay,ACTSO4,R2SO4, &
    12                    & DENSO4,f_r_wet,RRSI,Vbin,FL,ASO4,DNDR)
     168                   & DENSO4,f_r_wet,FL,ASO4,DNDR)
    13169!
    14170!     INPUT:
    … …  
    22178      USE infotrac_phy
    23179      USE YOMCST, ONLY : RPI
     180      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI
    24181
    25182      IMPLICIT NONE
    … …  
    33190      REAL DENSO4
    34191      REAL f_r_wet
    35       REAL RRSI(nbtr_bin)
    36       REAL Vbin(nbtr_bin)
    37 
     192     
    38193      ! output variables
    39194      REAL FL(nbtr_bin)
    … …  
    48203      REAL ATOT,AH2O
    49204      REAL RRSI_wet(nbtr_bin)
    50       REAL Vbin_wet(nbtr_bin)
    51       REAL MH2SO4,MH2O,BOLZ,FPATH
     205      REAL FPATH
    52206
    53207! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
    … …  
    57211!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE
    58212
    59 !                   WEIGHT OF 1 MOLEC IN G
    60       MH2O  =1000.*mH2Omol !18.016*1.66E-24
    61       MH2SO4=1000.*mH2SO4mol !98.082*1.66E-24
    62 !                   BOLTZMANN CONSTANTE IN DYN.CM/K
    63       BOLZ  =1.381E-16
    64213!                   MOLECULAR ACCOMODATION OF H2SO4
    65 !     raes and van dingen
    66       ALPHA =0.1   
     214!     H2SO4 accommodation coefficient [condensation/evaporation]
     215      ALPHA = ALPH2SO4
    67216!      FPLAIR=(2.281238E-5)*TAIR/PAIR
    68217!     1.E2 (m to cm),
    69218      CST=1.E2*2.281238E-5
    70219
    71       ! compute local wet particle radius and volume
     220      ! compute local wet particle radius [cm]
    72221      RRSI_wet(:)=RRSI(:)*f_r_wet
    73       Vbin_wet(:)=Vbin(:)*f_r_wet**3
    74 
     222     
    75223!     Pruppa and Klett
    76224      FPATH=CST*t_seri/pplay
    … …  
    138286
    139287!********************************************************************
    140       SUBROUTINE cond_evap_part(dt,FL,ASO4,f_r_wet,RRSI,Vbin,tr_seri)
     288      SUBROUTINE condens_evapor_part(dt,FL,ASO4,f_r_wet,tr_seri)
    141289
    142290      USE aerophys
    143291      USE infotrac_phy
    144292      USE YOMCST, ONLY : RPI
    145 
     293      USE strataer_local_var_mod, ONLY : RRSI,Vbin
     294     
    146295      IMPLICIT NONE
    147296
    … …  
    151300      REAL ASO4(nbtr_bin)
    152301      REAL f_r_wet
    153       REAL RRSI(nbtr_bin)
    154       REAL Vbin(nbtr_bin)
    155 
     302     
    156303      ! output variables
    157304      REAL tr_seri(nbtr)
    158 
     305     
    159306      ! local variables
    160307      REAL tr_seri_new(nbtr)
    … …  
    211358      tr_seri(:)=tr_seri_new(:)
    212359
    213       END SUBROUTINE cond_evap_part
     360      END SUBROUTINE condens_evapor_part
    214361
    215362END MODULE cond_evap_tstep_mod

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

    r4601 r5202  
    88  USE aerophys
    99  USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_H2SO4_strat
    10   USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, f_r_wet
     10  USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, &
     11       f_r_wet, R2SO4B, DENSO4B, f_r_wetB
    1112  USE nucleation_tstep_mod
    1213  USE cond_evap_tstep_mod
    … …  
    1415  USE YOMCST, ONLY : RPI, RD, RG
    1516  USE print_control_mod, ONLY: lunout
    16   USE strataer_local_var_mod
     17  USE strataer_local_var_mod ! contains also RRSI and Vbin
    1718 
    1819  IMPLICIT NONE
    … …  
    3536  REAL                      :: ntot !total number of molecules in the critical cluster (ntot>4)
    3637  REAL                      :: x    ! molefraction of H2SO4 in the critical cluster     
    37   REAL Vbin(nbtr_bin)
    3838  REAL a_xm, b_xm, c_xm
    3939  REAL PDT, dt
    4040  REAL H2SO4_init
    4141  REAL ACTSO4(klon,klev)
    42   REAL RRSI(nbtr_bin)
    4342  REAL nucl_rate
    4443  REAL cond_evap_rate
    … …  
    4847  REAL DNDR(nbtr_bin)
    4948  REAL H2SO4_sat
    50 
    51   DO it=1,nbtr_bin
    52     Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0
    53   ENDDO
    54 
     49  REAL R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin)
     50 
    5551  !coefficients for H2SO4 density parametrization used for nucleation if ntot<4
    5652  a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + &
    … …  
    6157       & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 )))))
    6258
    63   ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
    64   CALL STRAACT(ACTSO4)
    65 
    66   ! compute particle radius in cm RRSI from diameter in m
    67   DO it=1,nbtr_bin
    68     RRSI(it)=mdw(it)/2.*100.
    69   ENDDO
    70 
     59  IF(.not.flag_new_strat_compo) THEN
     60     ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
     61     CALL STRAACT(ACTSO4)
     62  ENDIF
     63 
    7164  DO ilon=1, klon
    7265!
    … …  
    10497      ENDIF
    10598      ! compute cond/evap rate in kg(H2SO4)/kgA/s
    106       CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
    107              & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
    108              & RRSI,Vbin,FL,ASO4,DNDR)
     99      IF(flag_new_strat_compo) THEN
     100         R2SO4ik(:)   = R2SO4B(ilon,ilev,:)
     101         DENSO4ik(:)  = DENSO4B(ilon,ilev,:)
     102         f_r_wetik(:) = f_r_wetB(ilon,ilev,:)
     103         CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
     104              & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
     105              & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
     106      ELSE
     107         CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
     108              & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
     109              & FL,ASO4,DNDR)
     110      ENDIF
    109111      ! Compute H2SO4 saturate vapor for big particules
    110112      H2SO4_sat = DNDR(nbtr_bin)/(pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol)
    … …  
    127129      tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt)
    128130      ! apply cond to bins
    129       CALL cond_evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
     131      CALL condens_evapor_part(dt,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:))
    130132      ! apply nucl. to bins
    131       CALL nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri(ilon,ilev,:))
     133      CALL nucleation_part(nucl_rate,ntot,x,dt,tr_seri(ilon,ilev,:))
    132134      ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
    133135      budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
    … …  
    142144        & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol
    143145    ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys)
    144     CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
    145            & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
    146            & RRSI,Vbin,FL,ASO4,DNDR)
     146    IF(flag_new_strat_compo) THEN
     147       CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
     148            & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
     149            & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
     150    ELSE
     151       CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
     152            & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
     153            & FL,ASO4,DNDR)
     154    ENDIF
    147155    ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet
    148156    DO it=1,nbtr_bin
    … …  
    159167    tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys)
    160168    ! apply evap to bins
    161     CALL cond_evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
     169    CALL condens_evapor_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:))
    162170    ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
    163171    budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &

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

    r3677 r5202  
    1616
    1717  USE phys_local_var_mod, ONLY: tr_seri, mdw, alpha_bin, piz_bin, cg_bin
    18   USE aerophys
     18  USE aerophys, ONLY: dens_aer_dry, dens_aer_ref, V_rat
    1919  USE aero_mod
    2020  USE infotrac_phy, ONLY : nbtr, nbtr_bin, nbtr_sulgas, id_SO2_strat
    … …  
    226226    40000.000,    0.2500,   1.48400,   1.0000E-08, &
    227227    50000.000,    0.2000,   1.49800,   1.0000E-08 /), (/nb_lambda_h2so4,4/), order=(/2,1/) )
    228 
    229   !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
    230     mdw(1)=mdwmin
    231     IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
    232       mdw(2)=mdw(1)*2.**(1./3.)
    233       DO it=3, nbtr_bin
    234         mdw(it)=mdw(it-1)*V_rat**(1./3.)
    235       ENDDO
    236     ELSE
    237       DO it=2, nbtr_bin
    238         mdw(it)=mdw(it-1)*V_rat**(1./3.)
    239       ENDDO
    240     ENDIF
    241     WRITE(lunout,*) 'init mdw=', mdw
    242 
     228     
    243229    !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K
    244230    DO bin_number=1, nbtr_bin

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

    r4912 r5202  
    7070!--------------------------------------------------------------------------------------------------
    7171
    72 SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri)
     72SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,tr_seri)
    7373
    7474  USE aerophys
    7575  USE infotrac_phy
    76 
     76  USE strataer_local_var_mod, ONLY : Vbin
     77 
    7778  IMPLICIT NONE
    7879
    … …  
    8283  REAL x    ! mole raction of H2SO4 in the critical cluster
    8384  REAL dt
    84   REAL Vbin(nbtr_bin)
    85 
     85 
    8686  ! output variable
    8787  REAL tr_seri(nbtr)

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

    r4767 r5202  
    5151 
    5252  !============= NUCLEATION VARS =============
     53  ! MOLECULAR ACCOMODATION OF H2SO4 (Raes and Van Dingen)
     54  REAL,SAVE    :: ALPH2SO4               ! H2SO4 accommodation  coefficient [condensation/evaporation]
     55  !$OMP THREADPRIVATE(ALPH2SO4)
     56 
    5357  ! flag to constraint nucleation rate in a lat/pres box
    5458  LOGICAL,SAVE :: flag_nuc_rate_box      ! Nucleation rate limit or not to a lat/pres
    … …  
    6468  INTEGER,SAVE :: flh2o  ! ds stratemit : flh2o =0 (tr_seri), flh2o=1 (dq)
    6569  !$OMP THREADPRIVATE(flh2o)
    66 !  REAL,ALLOCATABLE,SAVE    :: d_q_emiss(:,:)
    67 !  !$OMP THREADPRIVATE(d_q_emiss)
    6870 
    6971  REAL,ALLOCATABLE,SAVE    :: budg_emi(:,:)            !DIMENSION(klon,n)
    … …  
    144146  !$OMP THREADPRIVATE(day_emit_roc)
    145147 
     148  REAL,ALLOCATABLE,SAVE    :: RRSI(:) ! radius [cm] for each aerosol size
     149  REAL,ALLOCATABLE,SAVE    :: Vbin(:) ! volume [m3] for each aerosol size 
     150  !$OMP THREADPRIVATE(RRSI, Vbin)
    146151  REAL,SAVE    :: dlat, dlon             ! delta latitude and d longitude of grid in degree
    147152  !$OMP THREADPRIVATE(dlat, dlon)
    … …  
    153158    USE print_control_mod, ONLY : lunout
    154159    USE mod_phys_lmdz_para, ONLY : is_master
    155     USE infotrac_phy, ONLY: id_OCS_strat,id_SO2_strat,id_H2SO4_strat,nbtr_sulgas
     160    USE infotrac_phy, ONLY: id_OCS_strat,id_SO2_strat,id_H2SO4_strat,nbtr_sulgas,nbtr_bin
     161    USE phys_local_var_mod, ONLY : mdw
     162    USE aerophys, ONLY: mdwmin, V_rat
     163    USE YOMCST  , ONLY : RPI
     164   
     165    INTEGER :: it
    156166   
    157167    WRITE(lunout,*) 'IN STRATAER_LOCAL_VAR INIT WELCOME!'
    … …  
    185195   
    186196    ! nuc init
     197    ALPH2SO4 = 0.1
    187198    flag_nuc_rate_box = .FALSE.
    188199    nuclat_min=0  ; nuclat_max=0
    … …  
    238249    ENDIF ! if master
    239250   
     251    !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
     252    mdw(1)=mdwmin
     253    IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
     254       mdw(2)=mdw(1)*2.**(1./3.)
     255       DO it=3, nbtr_bin
     256          mdw(it)=mdw(it-1)*V_rat**(1./3.)
     257       ENDDO
     258    ELSE
     259       DO it=2, nbtr_bin
     260          mdw(it)=mdw(it-1)*V_rat**(1./3.)
     261       ENDDO
     262    ENDIF
     263    IF (is_master) WRITE(lunout,*) 'init mdw=', mdw
     264   
     265    !   compute particle radius RRSI [cm] and volume Vbin [m3] from diameter mdw [m]
     266    ALLOCATE(RRSI(nbtr_bin), Vbin(nbtr_bin))
     267   
     268    DO it=1,nbtr_bin
     269       !     [cm]
     270       RRSI(it)=mdw(it)/2.*100.
     271       !     [m3]
     272       Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0
     273    ENDDO
     274   
     275    IF (is_master) THEN
     276       WRITE(lunout,*) 'init RRSI=', RRSI
     277       WRITE(lunout,*) 'init Vbin=', Vbin
     278    ENDIF
     279   
    240280    WRITE(lunout,*) 'IN STRATAER INIT END'
    241281   

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

    r4601 r5202  
    1313    USE print_control_mod, ONLY : lunout
    1414    USE mod_phys_lmdz_para, ONLY : is_master
    15     USE strataer_local_var_mod, ONLY: flag_nuc_rate_box,nuclat_min,nuclat_max,nucpres_min,nucpres_max
     15    USE strataer_local_var_mod, ONLY: ALPH2SO4,flag_nuc_rate_box,nuclat_min,nuclat_max, &
     16         nucpres_min,nucpres_max
    1617   
    1718    !Config Key  = flag_nuc_rate_box
    … …  
    3031    CALL getin_p('nucpres_max',nucpres_max)
    3132   
     33    ! Read argument H2SO4 accommodation  coefficient [condensation/evaporation]
     34    CALL getin_p('alph2so4',ALPH2SO4)
     35   
    3236    !============= Print params =============
    3337    IF (is_master) THEN
     38       WRITE(lunout,*) 'IN STRATAER_NUC : ALPH2SO4 = ',alph2so4
    3439       WRITE(lunout,*) 'IN STRATAER_NUC : flag_nuc_rate_box = ',flag_nuc_rate_box
    3540       IF (flag_nuc_rate_box) THEN

  • 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

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

    r4769 r5202  
    99       presnivs, xlat, xlon, pphis, pphi, &
    1010       t_seri, pplay, paprs, sh, rh, tr_seri)
    11 
     11   
    1212    USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
    13         & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part
    14 
     13        & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
     14        & R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode, reff_sulfate
     15   
    1516    USE dimphy
    1617    USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
    … …  
    5657    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
    5758    REAL                                   :: m_aer               ! aerosol mass
    58     INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
     59    INTEGER                                :: it, k, i, j, ilon, ilev, itime, i_int, ieru
    5960    LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
    6061    REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
    … …  
    8283    INTEGER                                :: injdur_sai          ! injection duration for SAI case [days]
    8384    INTEGER                                :: yr, is_bissext
     85    REAL                                   :: samoment2, samoment3! 2nd and 3rd order moments of size distribution
    8486
    8587    IF (is_mpi_root .AND. flag_verbose_strataer) THEN
    … …  
    8890    ENDIF
    8991   
     92    !   radius [m]
    9093    DO it=1, nbtr_bin
    9194      r_bin(it)=mdw(it)/2.
    … …  
    117120
    118121    IF(flag_new_strat_compo) THEN
    119        IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Tabazadeh 1994', flag_new_strat_compo
    120        ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) : binary routine (from reprobus)
    121        ! H2SO4 mass fraction in aerosol (%) from Tabazadeh et al. (1994).
    122        CALL stracomp_bin(sh,t_seri,pplay)
    123        
    124        ! aerosol density (gr/cm3) - from Tabazadeh
    125        CALL denh2sa_taba(t_seri)
     122       IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
     123       ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
     124       CALL stracomp_kelvin(sh,t_seri,pplay)
    126125    ELSE
    127        IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Bekki 2D model', flag_new_strat_compo
     126       IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
    128127       ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
    129128       ! H2SO4 mass fraction in aerosol (%)
    … …  
    132131       ! aerosol density (gr/cm3)
    133132       CALL denh2sa(t_seri)
     133       
     134       ! compute factor for converting dry to wet radius (for every grid box)
     135       f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
    134136    ENDIF
    135137   
    136 ! compute factor for converting dry to wet radius (for every grid box)
    137     f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
    138 
    139138!--calculate mass of air in every grid box
    140139    DO ilon=1, klon
    … …  
    348347    ENDDO
    349348   
     349!--compute
     350!     sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
     351!     SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
     352!     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
     353!     nd_mode: DRY(?) particle concentration in different modes (part/m3)
     354     sulfmmr(:,:)=0.0
     355     SAD_sulfate(:,:)=0.0
     356     sulfmmr_mode(:,:,:)=0.0
     357     nd_mode(:,:,:)=0.0
     358     reff_sulfate(:,:)=0.0
     359     
     360     DO i=1,klon
     361        DO j=1,klev
     362           samoment2=0.0
     363           samoment3=0.0
     364           DO it=1, nbtr_bin
     365              !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
     366              !assume that particles consist of ammonium sulfate at the surface (132g/mol)
     367              !and are dry at T = 20 deg. C and 50 perc. humidity
     368             
     369              !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
     370              !     equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
     371              sulfmmr_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &        ! [DRY part/kgA in bin it]
     372                   &  *(4./3.)*RPI*(mdw(it)/2.)**3.   &                   ! [mdw: dry diameter in m]
     373                   &  *dens_aer_dry                                       ! [dry aerosol mass density in kg/m3]
     374             
     375              !     sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio)
     376              !     equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
     377              sulfmmr(i,j) = sulfmmr(i,j) + sulfmmr_mode(i,j,it)
     378             
     379              !     nd_mode: particle concentration in different modes (DRY part/m3)
     380              nd_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &             ! [DRY part/kgA in bin it]
     381                   & *pplay(i,j)/t_seri(i,j)/RD                           ! [air mass concentration in kg air /m3A]
     382             
     383              IF(flag_new_strat_compo) THEN
     384                 !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
     385                 SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
     386                      &  *4.*RPI*( mdw(it)*f_r_wetB(i,j,it)/2. )**2. &       ! [WET SA of part it in m2]
     387                      &  *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
     388!    samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
     389                 samoment2 = samoment2 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
     390                      &  *( mdw(it)*f_r_wetB(i,j,it)/2. )**2.                     ! [WET SA of part it in m2]
     391!    samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
     392                 samoment3 = samoment3 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
     393                      &  *( mdw(it)*f_r_wetB(i,j,it)/2. )**3.                     ! [WET SA of part it in m2]
     394              ELSE
     395!     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
     396                 SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
     397                      &  *4.*RPI*( mdw(it)*f_r_wet(i,j)/2. )**2. &           ! [WET SA of part it in m2]
     398                      &  *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
     399!    samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
     400                 samoment2 = samoment2 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
     401                      &  *( mdw(it)*f_r_wet(i,j)/2. )**2.                          ! [WET SA of part it in m2]
     402!    samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
     403                 samoment3 = samoment3 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
     404                      &  *( mdw(it)*f_r_wet(i,j)/2. )**3.                          ! [WET SA of part it in m2]
     405              ENDIF
     406           ENDDO
     407!     reff_sulfate: effective radius of WET sulfate aerosols (cm)
     408           reff_sulfate(i,j) = (samoment3 / samoment2) &
     409                & *1.e2                                              ! conversion from m to cm
     410        ENDDO
     411     ENDDO
     412     
    350413  END SUBROUTINE traccoag
    351414
Note: See TracChangeset for help on using the changeset viewer.