!
! $Id: cond_evap_tstep_mod.F90 5202 2024-09-20 10:32:04Z abarral $
!
MODULE cond_evap_tstep_mod

! based on UPMC aerosol model by Slimane Bekki
! adapted for stratospheric sulfate aerosol in LMDZ by Christoph Kleinschmitt 

CONTAINS

      SUBROUTINE condens_evapor_rate_kelvin(R2SO4G,t_seri,pplay,R2SO4, &
          & DENSO4,f_r_wet,R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
!
!     INPUT: 
!     R2SO4G: number density of gaseous H2SO4 [molecules/cm3]
!     t_seri: temperature (K)
!     pplay: pressure (Pa)
!     R2SO4: aerosol H2SO4 weight fraction (percent) - flat surface (does not depend on aerosol size)
!     DENSO4: aerosol density (gr/cm3) 
!     f_r_wet: factor for converting dry to wet radius 
!        assuming 'flat surface' composition (does not depend on aerosol size)
!     variables that depends on aerosol size because of Kelvin effect
!     R2SO4Gik: number density of gaseous H2SO4 [molecules/cm3] - depends on aerosol size
!     DENSO4ik: aerosol density (gr/cm3) - depends on aerosol size
!     f_r_wetik: factor for converting dry to wet radius - depends on aerosol size
!     RRSI: radius [cm]

      USE aerophys
      USE infotrac_phy
      USE YOMCST, ONLY : RPI
      USE sulfate_aer_mod, ONLY : wph2so4, surftension, solh2so4, rpmvh2so4
      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI
      
      IMPLICIT NONE
      
      REAL, PARAMETER :: third=1./3.
      
      ! input variables
      REAL            :: R2SO4G !H2SO4 number density [molecules/cm3]
      REAL            :: t_seri
      REAL            :: pplay
      REAL            :: R2SO4
      REAL            :: DENSO4
      REAL            :: f_r_wet
      REAL            :: R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin)
      
      ! output variables
      REAL            :: FL(nbtr_bin)
      REAL            :: ASO4(nbtr_bin)
      REAL            :: DNDR(nbtr_bin)
      
      ! local variables
      INTEGER            :: IK
      REAL            :: ALPHA,CST
      REAL            :: WH2(nbtr_bin)
      REAL            :: RP,VTK,AA,FL1,RKNUD
      REAL            :: DND
      REAL            :: ATOT,AH2O
      REAL            :: RRSI_wet(nbtr_bin)
      REAL            :: FPATH, WPP, XA, FKELVIN
      REAL            :: surtens, mvh2so4, temp
      
! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
!    ------------------------------------------------------------------
!                                  EXCEPT CN
!       RK:H2SO4 WEIGHT PERCENT DOESN'T CHANGE
!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE

!                   MOLECULAR ACCOMODATION OF H2SO4
!     H2SO4 accommodation  coefficient [condensation/evaporation]
      ALPHA = ALPH2SO4
!      FPLAIR=(2.281238E-5)*TAIR/PAIR
!     1.E2 (m to cm), 
      CST=1.E2*2.281238E-5
!     same expression as in coagulate
!     in coagulate: mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
!     mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15)
!     mnfrpth=2.28E-5*t_seri/pplay
      
      temp = min( max(t_seri, 190.), 300.) ! 190K <= temp <= 300K
      
      RRSI_wet(:)=RRSI(:)*f_r_wetik(:)

!     Pruppa and Klett
      FPATH=CST*t_seri/pplay
    
!     H2SO4 mass fraction in aerosol
      WH2(:)=R2SO4ik(:)*1.0E-2

!                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
!                               AYERS ET AL (1980)
!                                  (MU-MU0)
!      RP=-10156.0/t_seri +16.259-(ACTSO4*4.184)/(8.31441*t_seri)
!                                  DROPLET H2SO4 PRESSURE IN DYN.CM-2
!      RP=EXP(RP)*1.01325E6/0.086
!!      RP=EXP(RP)*1.01325E6
!                                  H2SO4 NUMBER DENSITY NEAR DROPLET

!      DND=RP*6.02E23/(8.31E7*t_seri)

!                                 KELVIN EFFECT FACTOR
!CK 20160613: bug fix, removed factor 250 (from original code by S. Bekki)
!!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri*250.0)
!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri)

!                                  MEAN KINETIC VELOCITY
!     DYN*CM*K/(K*GR)=(CM/SEC2)*CM
!                                  IN CM/SEC
      VTK=SQRT(8.0*BOLZ*t_seri/(RPI*MH2SO4))
!                                 KELVIN EFFECT FACTOR

!     Loop on bin radius (RRSI in cm)
      DO IK=1,nbtr_bin

      IF(R2SO4ik(IK) > 0.0) THEN

!       h2so4 mass fraction (0<wpp<1)
        wpp=R2SO4ik(IK)*1.e-2    
        xa=18.*wpp/(18.*wpp+98.*(1.-wpp)) 
!       equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
        DND=solh2so4(t_seri,xa)
!          KELVIN EFFECT:  
!       surface tension (mN/m=1.e-3.kg/s2) = f(T,h2so4 mole fraction)
        surtens=surftension(temp,xa)
!       partial molar volume of h2so4 (cm3.mol-1 =1.e-6.m3.mol-1)
        mvh2so4= rpmvh2so4(temp,R2SO4ik(IK))
!       Kelvin factor (MKS)
        fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2so4/ (1.e-2*RRSI_wet(IK)*rgas*temp) )
!                              
        DNDR(IK) =DND*fkelvin

        FL1=RPI*ALPHA*VTK*(R2SO4G-DNDR(IK))

!       TURCO(1979) FOR HNO3:ALH2SO4 CONDENSATION= ALH2SO4 EVAPORATION
!       RPI*R2*VTK IS EQUIVALENT TO DIFFUSION COEFFICIENT
!       EXTENSION OF THE RELATION FOR DIFFUSION KINETICS
!       KNUDSEN NUMBER FPATH/RRSI
!       NEW VERSION (SEE NOTES)
        RKNUD=FPATH/RRSI_wet(IK)
!       SENFELD
        FL(IK)=FL1*RRSI_wet(IK)**2*( 1.0 +RKNUD ) & 
     &     /( 1.0 +ALPHA/(2.0*RKNUD) +RKNUD )
!       TURCO
!        RL= (4.0/3.0 +0.71/RKNUD)/(1.0+1.0/RKNUD)
!     *         +4.0*(1.0-ALPHA)/(3.0*ALPHA)
!        FL=FL1*RRSI(IK)*RRSI(IK)
!     *         /( (3.0*ALPHA/4.0)*(1.0/RKNUD+RL*ALPHA) )

!                         INITIAL NUMBER OF H2SO4 MOLEC OF 1 DROPLET
        ATOT=4.0*RPI*DENSO4ik(IK)*(RRSI_wet(IK)**3)/3.0 !attention: g and cm
        ASO4(IK)=WH2(IK)*ATOT/MH2SO4 !attention: g
!        ATOT=4.0*RPI*dens_aer(I,J)/1000.*(RRSI(IK)**3)/3.0
!        ASO4=mfrac_H2SO4*ATOT/MH2SO4
!                        INITIAL NUMBER OF H2O MOLEC OF 1 DROPLET
        AH2O=(1.0-WH2(IK))*ATOT/MH2O !attention: g

!       CHANGE OF THE NUMBER OF H2SO4 MOLEC OF 1 DROPLET DURING DT
!       IT IS FOR KEM BUT THERE ARE OTHER WAYS

      ENDIF  

      ENDDO !loop over bins

      END SUBROUTINE condens_evapor_rate_kelvin
      
!********************************************************************
      SUBROUTINE condens_evapor_rate(R2SO4G,t_seri,pplay,ACTSO4,R2SO4, &
                   & DENSO4,f_r_wet,FL,ASO4,DNDR)
!
!     INPUT: 
!     R2SO4: aerosol H2SO4 weight fraction (percent)
!     ACTSO4: H2SO4 activity
!     R2SO4G: number density of gaseous H2SO4 [molecules/cm3]
!     t_seri: temperature (K)
!     DENSO4: aerosol density (gr/cm3)

      USE aerophys
      USE infotrac_phy
      USE YOMCST, ONLY : RPI
      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI

      IMPLICIT NONE

      ! input variables
      REAL R2SO4G !H2SO4 number density [molecules/cm3]
      REAL t_seri
      REAL pplay
      REAL ACTSO4
      REAL R2SO4
      REAL DENSO4
      REAL f_r_wet
      
      ! output variables
      REAL FL(nbtr_bin)
      REAL ASO4(nbtr_bin)
      REAL DNDR(nbtr_bin)

      ! local variables
      INTEGER IK
      REAL ALPHA,CST
      REAL WH2,RP,VTK,AA,FL1,RKNUD
      REAL DND
      REAL ATOT,AH2O
      REAL RRSI_wet(nbtr_bin)
      REAL FPATH

! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
!    ------------------------------------------------------------------
!                                  EXCEPT CN
!       RK:H2SO4 WEIGHT PERCENT DOESN'T CHANGE
!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE

!                   MOLECULAR ACCOMODATION OF H2SO4
!     H2SO4 accommodation coefficient [condensation/evaporation]
      ALPHA = ALPH2SO4
!      FPLAIR=(2.281238E-5)*TAIR/PAIR
!     1.E2 (m to cm), 
      CST=1.E2*2.281238E-5

      ! compute local wet particle radius [cm]
      RRSI_wet(:)=RRSI(:)*f_r_wet
      
!     Pruppa and Klett
      FPATH=CST*t_seri/pplay

      
!     H2SO4 mass fraction in aerosol
      WH2=R2SO4*1.0E-2
      IF(WH2.EQ.0.0) RETURN
!                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
!                               AYERS ET AL (1980)
!                                  (MU-MU0)
      RP=-10156.0/t_seri +16.259-(ACTSO4*4.184)/(8.31441*t_seri)
!                                  DROPLET H2SO4 PRESSURE IN DYN.CM-2
      RP=EXP(RP)*1.01325E6/0.086
!      RP=EXP(RP)*1.01325E6
!                                  H2SO4 NUMBER DENSITY NEAR DROPLET
!     R=8.31E7 DYN.CM.MOL-1*K-1
!                               R/AVOGADRO NUMBER=DYN.CM.MOLEC-1*K-1
      DND=RP*6.02E23/(8.31E7*t_seri)
!                                  MEAN KINETIC VELOCITY
!     DYN*CM*K/(K*GR)=(CM/SEC2)*CM
!                                  IN CM/SEC
      VTK=SQRT(8.0*BOLZ*t_seri/(RPI*MH2SO4))
!                                 KELVIN EFFECT FACTOR
!CK 20160613: bug fix, removed factor 250 (from original code by S. Bekki)
!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri*250.0)
      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri)

!     Loop on bin radius (RRSI in cm)
      DO IK=1,nbtr_bin
!                                 KELVIN EFFECT
        DNDR(IK) =DND*EXP(AA/RRSI_wet(IK))

        FL1=RPI*ALPHA*VTK*(R2SO4G-DNDR(IK))

!       TURCO(1979) FOR HNO3:ALH2SO4 CONDENSATION= ALH2SO4 EVAPORATION
!       RPI*R2*VTK IS EQUIVALENT TO DIFFUSION COEFFICIENT
!       EXTENSION OF THE RELATION FOR DIFFUSION KINETICS
!       KNUDSEN NUMBER FPATH/RRSI
!       NEW VERSION (SEE NOTES)
        RKNUD=FPATH/RRSI_wet(IK)
!       SENFELD
        FL(IK)=FL1*RRSI_wet(IK)**2*( 1.0 +RKNUD ) & 
     &     /( 1.0 +ALPHA/(2.0*RKNUD) +RKNUD )
!       TURCO
!        RL= (4.0/3.0 +0.71/RKNUD)/(1.0+1.0/RKNUD)
!     *         +4.0*(1.0-ALPHA)/(3.0*ALPHA)
!        FL=FL1*RRSI(IK)*RRSI(IK)
!     *         /( (3.0*ALPHA/4.0)*(1.0/RKNUD+RL*ALPHA) )

!                         INITIAL NUMBER OF H2SO4 MOLEC OF 1 DROPLET
        ATOT=4.0*RPI*DENSO4*(RRSI_wet(IK)**3)/3.0 !attention: g and cm
        ASO4(IK)=WH2*ATOT/MH2SO4 !attention: g
!        ATOT=4.0*RPI*dens_aer(I,J)/1000.*(RRSI(IK)**3)/3.0
!        ASO4=mfrac_H2SO4*ATOT/MH2SO4
!                        INITIAL NUMBER OF H2O MOLEC OF 1 DROPLET
        AH2O=(1.0-WH2)*ATOT/MH2O !attention: g

!       CHANGE OF THE NUMBER OF H2SO4 MOLEC OF 1 DROPLET DURING DT
!       IT IS FOR KEM BUT THERE ARE OTHER WAYS

      ENDDO !loop over bins

      END SUBROUTINE condens_evapor_rate

!********************************************************************
      SUBROUTINE condens_evapor_part(dt,FL,ASO4,f_r_wet,tr_seri)

      USE aerophys
      USE infotrac_phy
      USE YOMCST, ONLY : RPI
      USE strataer_local_var_mod, ONLY : RRSI,Vbin
      
      IMPLICIT NONE

      ! input variables
      REAL dt
      REAL FL(nbtr_bin)
      REAL ASO4(nbtr_bin)
      REAL f_r_wet
      
      ! output variables
      REAL tr_seri(nbtr)
      
      ! local variables
      REAL tr_seri_new(nbtr)
      INTEGER IK,JK,k
      REAL Vnew
      REAL RRSI_wet(nbtr_bin)
      REAL Vbin_wet(nbtr_bin)
      REAL sum_IK(nbtr_bin)
      REAL ff(nbtr_bin,nbtr_bin)

      tr_seri_new(:)=tr_seri(:)

      ! compute local wet particle radius and volume
      RRSI_wet(:)=RRSI(:)*f_r_wet
      Vbin_wet(:)=Vbin(:)*f_r_wet**3 *1.e6 !Vbin_wet in cm3 (as Vnew)

      ! compute distribution factor for particles of intermediate size (from Jacobson 1994, equation 13)
      DO IK=1,nbtr_bin
        Vnew=4.0*RPI*(RRSI_wet(IK)**3)/3.0*(1.+FL(IK)*dt/ASO4(IK))
        ff(IK,:)=0.0 
        DO k=1, nbtr_bin 
          IF (k.LE.(nbtr_bin-1)) THEN 
            IF (Vbin_wet(k).LE.Vnew.AND.Vnew.LT.Vbin_wet(k+1)) THEN
              ff(IK,k)= Vbin_wet(k)/Vnew*(Vbin_wet(k+1)-Vnew)/(Vbin_wet(k+1)-Vbin_wet(k))
            ENDIF
          ENDIF
          IF (k.EQ.1.AND.Vnew.LE.Vbin_wet(k)) THEN
            ff(IK,k)= 1.
          ENDIF
          IF (k.GT.1) THEN 
            IF (Vbin_wet(k-1).LT.Vnew.AND.Vnew.LT.Vbin_wet(k)) THEN
              ff(IK,k)= 1.-ff(IK,k-1)
            ENDIF
          ENDIF
          IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin_wet(k)) THEN
            ff(IK,k)= 1.
          ENDIF
        ENDDO
        ! correction of ff for volume conservation
        DO k=1, nbtr_bin         
          ff(IK,k)=ff(IK,k)*Vnew/Vbin_wet(k)
        ENDDO
      ENDDO !loop over bins

      DO IK=1, nbtr_bin
        sum_IK(IK)=0.0
        DO JK=1, nbtr_bin
          sum_IK(IK)=sum_IK(IK)+tr_seri(JK+nbtr_sulgas)*ff(JK,IK)
        ENDDO
        ! compute new particle concentrations
        tr_seri_new(IK+nbtr_sulgas)=sum_IK(IK)
      ENDDO

      tr_seri(:)=tr_seri_new(:)

      END SUBROUTINE condens_evapor_part

END MODULE cond_evap_tstep_mod