LCOV - code coverage report
Current view: | top level - physics/cam - hetfrz_classnuc.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 168 | 174 | 96.6 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 5 | 5 | 100.0 %

              Line data    Source code
           1             : module hetfrz_classnuc
           2             :
           3             : !-----------------------------------------------------------------------
           4             : !
           5             : ! Purpose: Calculate heterogeneous freezing rates from classical nucleation theory
           6             : !
           7             : ! Public interfaces:
           8             : !
           9             : ! hetfrz_classnuc_init
          10             : ! hetfrz_classnuc_calc
          11             : !
          12             : ! Author:
          13             : !   Corinna Hoose, UiO, May 2009
          14             : !   Yong Wang and Xiaohong Liu, UWyo, 12/2012,
          15             : !   implement in CAM5 and constrain uncertain parameters using natural dust and
          16             : !   BC(soot) datasets.
          17             : !   Yong Wang and Xiaohong Liu, UWyo, 05/2013, implement the PDF-contact angle approach:
          18             : !   Y. Wang et al., Atmos. Chem. Phys., 2014. https://doi.org/10.5194/acp-14-10411-2014
          19             : !   Jack Chen, NCAR, 09/2015, modify calculation of dust activation fraction.
          20             : !
          21             : !-----------------------------------------------------------------------
          22             :
          23             : use shr_kind_mod,  only: r8 => shr_kind_r8
          24             : use wv_saturation, only: svp_water
          25             : use shr_spfn_mod,  only: erf => shr_spfn_erf
          26             :
          27             : use physconst,     only: pi, planck, boltz, mwso4, amu, pstd
          28             :
          29             : implicit none
          30             : private
          31             :
          32             : public :: hetfrz_classnuc_init, hetfrz_classnuc_calc
          33             :
          34             : real(r8) :: rair
          35             : real(r8) :: cpair
          36             : real(r8) :: rh2o
          37             : real(r8) :: rhoh2o
          38             : real(r8) :: mwh2o
          39             : real(r8) :: tmelt
          40             :
          41             : !*****************************************************************************
          42             : !                PDF theta model
          43             : !*****************************************************************************
          44             : ! some variables for PDF theta model
          45             : ! immersion freezing
          46             : !
          47             : ! With the original value of pdf_n_theta set to 101 the dust activation
          48             : ! fraction between -15 and 0 C could be overestimated.  This problem was
          49             : ! eliminated by increasing pdf_n_theta to 301.  To reduce the expense of
          50             : ! computing the dust activation fraction the integral is only evaluated
          51             : ! where dim_theta is non-zero.  This was determined to be between
          52             : ! dim_theta index values of 53 through 113.  These loop bounds are
          53             : ! hardcoded in the variables i1 and i2.
          54             :
          55             : logical, parameter :: pdf_imm_in = .true.
          56             : integer, parameter :: pdf_n_theta = 301
          57             : integer, parameter :: i1 = 53
          58             : integer, parameter :: i2 = 113
          59             :
          60             : real(r8) :: dim_theta(pdf_n_theta) = -huge(1._r8)
          61             : real(r8) :: pdf_imm_theta(pdf_n_theta) = 0.0_r8
          62             : real(r8) :: pdf_d_theta = -huge(1._r8)
          63             : real(r8) :: dim_f_imm(pdf_n_theta) = 0.0_r8
          64             :
          65             : integer  :: iulog
          66             :
          67             : real(r8), parameter :: n1 = 1.e19_r8           ! number of water molecules in contact with unit area of substrate [m-2]
          68             : real(r8), parameter :: rhplanck = 1._r8/planck
          69             : real(r8), parameter :: nus = 1.e13_r8          ! frequ. of vibration [s-1] higher freq. (as in P&K, consistent with Anupam's data)
          70             : real(r8), parameter :: rhwincloud = 0.98_r8    ! 98% RH in mixed-phase clouds (Korolev & Isaac, JAS 2006)
          71             :
          72             : logical, parameter :: tot_in = .false.
          73             :
          74             : real(r8) :: bc_limfac = -huge(1._r8)   ! soot ice nucleating fraction
          75             : real(r8) :: dust_limfac = -huge(1._r8) ! dust ice nucleating fraction
          76             :
          77             : !===================================================================================================
          78             : contains
          79             : !===================================================================================================
          80             :
          81        1536 : subroutine hetfrz_classnuc_init( &
          82             :    rair_in, cpair_in, rh2o_in, rhoh2o_in, mwh2o_in, &
          83             :    tmelt_in, iulog_in, bc_limfac_in, dust_limfac_in)
          84             :
          85             :    real(r8), intent(in) :: rair_in
          86             :    real(r8), intent(in) :: cpair_in
          87             :    real(r8), intent(in) :: rh2o_in
          88             :    real(r8), intent(in) :: rhoh2o_in
          89             :    real(r8), intent(in) :: mwh2o_in
          90             :    real(r8), intent(in) :: tmelt_in
          91             :    integer,  intent(in) :: iulog_in
          92             :    real(r8), intent(in) :: bc_limfac_in
          93             :    real(r8), intent(in) :: dust_limfac_in
          94             :
          95        1536 :    rair   = rair_in
          96        1536 :    cpair  = cpair_in
          97        1536 :    rh2o   = rh2o_in
          98        1536 :    rhoh2o = rhoh2o_in
          99        1536 :    mwh2o  = mwh2o_in
         100        1536 :    tmelt  = tmelt_in
         101        1536 :    iulog  = iulog_in
         102             :
         103        1536 :    bc_limfac = bc_limfac_in
         104        1536 :    dust_limfac = dust_limfac_in
         105             :
         106             :    ! Initialize all the PDF theta variables:
         107        1536 :    if (pdf_imm_in) then
         108             :       call hetfrz_classnuc_init_pdftheta()
         109             :    end if
         110             :
         111        1536 : end subroutine hetfrz_classnuc_init
         112             :
         113             : !===================================================================================================
         114             :
         115  1270492981 : subroutine hetfrz_classnuc_calc(ntypes, types,&
         116             :    deltat, T, p, supersatice,                 &
         117  1270492981 :    fn,                                        &
         118             :    r3lx, icnlx,                               &
         119             :    frzbcimm, frzduimm,                        &
         120             :    frzbccnt, frzducnt,                        &
         121             :    frzbcdep, frzdudep,                        &
         122  1270492981 :    hetraer, wact_factor, dstcoat,             &
         123  1270492981 :    total_aer_num, uncoated_aer_num,           &
         124  1270492981 :    total_interstitial_aer_num, total_cloudborne_aer_num, errstring)
         125             :
         126             :    integer, intent(in) :: ntypes
         127             :    character(len=*), intent(in) :: types(ntypes)
         128             :    real(r8), intent(in) :: deltat            ! timestep [s]
         129             :    real(r8), intent(in) :: T                 ! temperature [K]
         130             :    real(r8), intent(in) :: p                 ! pressure [Pa]
         131             :    real(r8), intent(in) :: supersatice       ! supersaturation ratio wrt ice at 100%rh over water [ ]
         132             :    real(r8), intent(in) :: r3lx              ! volume mean drop radius [m]
         133             :    real(r8), intent(in) :: icnlx             ! in-cloud droplet concentration [cm-3]
         134             :
         135             :    real(r8), intent(in) :: fn(ntypes)               ! fraction activated [ ] for cloud borne aerosol number
         136             :                                                     ! index values are 1:bc, 2:dust_a1, 3:dust_a3
         137             :    real(r8), intent(in) :: hetraer(ntypes)          ! bc and dust mass mean radius [m]
         138             :    real(r8), intent(in) :: wact_factor(ntypes)      ! water activity factor -- density*(1.-(OC+BC)/(OC+BC+SO4)) [mug m-3]
         139             :    real(r8), intent(in) :: dstcoat(ntypes)          ! coated fraction
         140             :    real(r8), intent(in) :: total_aer_num(ntypes)    ! total bc and dust number concentration(interstitial+cloudborne) [#/cm^3]
         141             :    real(r8), intent(in) :: uncoated_aer_num(ntypes) ! uncoated bc and dust number concentration(interstitial)
         142             :    real(r8), intent(in) :: total_interstitial_aer_num(ntypes) ! total bc and dust concentration(interstitial)
         143             :    real(r8), intent(in) :: total_cloudborne_aer_num(ntypes)   ! total bc and dust concentration(cloudborne)
         144             :
         145             :    real(r8), target, intent(out) :: frzbcimm           ! het. frz by BC immersion nucleation [cm-3 s-1]
         146             :    real(r8), target, intent(out) :: frzduimm           ! het. frz by dust immersion nucleation [cm-3 s-1]
         147             :    real(r8), target, intent(out) :: frzbccnt           ! het. frz by BC contact nucleation [cm-3 s-1]
         148             :    real(r8), target, intent(out) :: frzducnt           ! het. frz by dust contact nucleation [cm-3 s-1]
         149             :    real(r8), target, intent(out) :: frzbcdep           ! het. frz by BC deposition nucleation [cm-3 s-1]
         150             :    real(r8), target, intent(out) :: frzdudep           ! het. frz by dust deposition nucleation [cm-3 s-1]
         151             :
         152             :    character(len=*), intent(out) :: errstring
         153             :
         154             :    ! local variables
         155             :
         156             :    real(r8) :: tc
         157             :    real(r8) :: vwice
         158             :    real(r8) :: rhoice
         159             :    real(r8) :: sigma_iw                        ! [J/m2]
         160             :    real(r8) :: sigma_iv                        ! [J/m2]
         161             :    real(r8) :: eswtr                           ! [Pa]
         162             :
         163             :    real(r8) :: rgimm    ! critical germ size
         164             :    real(r8) :: rgdep
         165             :    real(r8) :: dg0dep   ! homogeneous energy of germ formation
         166             :    real(r8) :: dg0cnt
         167             :    real(r8) :: Adep     ! prefactors
         168             :    real(r8) :: Acnt
         169             :
         170             :    !********************************************************
         171             :    ! Hoose et al., 2010 fitting parameters
         172             :    !********************************************************
         173             :    !freezing parameters for immersion freezing
         174             :    !real(r8),parameter :: theta_imm_bc = 40.17         ! contact angle [deg], converted to rad later
         175             :    !real(r8),parameter :: dga_imm_bc = 14.4E-20        ! activation energy [J]
         176             :    !real(r8),parameter :: theta_imm_dust = 30.98       ! contact angle [deg], converted to rad later
         177             :    !real(r8),parameter :: dga_imm_dust = 15.7E-20      ! activation energy [J]
         178             :    !freezing parameters for deposition nucleation
         179             :    !real(r8),parameter :: theta_dep_dust = 12.7       ! contact angle [deg], converted to rad later !Zimmermann et al (2008), illite
         180             :    !real(r8),parameter :: dga_dep_dust = -6.21E-21    ! activation energy [J]
         181             :    !real(r8),parameter :: theta_dep_bc = 28.          ! contact angle [deg], converted to rad later !Moehler et al (2005), soot
         182             :    !real(r8),parameter :: dga_dep_bc = -2.E-19        ! activation energy [J]
         183             :    !********************************************************
         184             :    ! Wang et al., 2014 fitting parameters
         185             :    !********************************************************
         186             :    ! freezing parameters for immersion freezing
         187             :    real(r8),parameter :: theta_imm_bc = 48.0_r8            ! contact angle [deg], converted to rad later !DeMott et al (1990)
         188             :    real(r8),parameter :: dga_imm_bc = 14.15E-20_r8         ! activation energy [J]
         189             :    real(r8),parameter :: theta_imm_dust = 46.0_r8          ! contact angle [deg], converted to rad later !DeMott et al (2011) SD
         190             :    real(r8),parameter :: dga_imm_dust = 14.75E-20_r8       ! activation energy [J]
         191             :    ! freezing parameters for deposition nucleation
         192             :    real(r8),parameter :: theta_dep_dust = 20.0_r8          ! contact angle [deg], converted to rad later !Koehler et al (2010) SD
         193             :    real(r8),parameter :: dga_dep_dust = -8.1E-21_r8        ! activation energy [J]
         194             :    real(r8),parameter :: theta_dep_bc = 28._r8             ! contact angle [deg], converted to rad later !Moehler et al (2005), soot
         195             :    real(r8),parameter :: dga_dep_bc = -2.E-19_r8           ! activation energy [J]
         196             :
         197             :    ! form factor
         198             :    ! only consider flat surfaces due to uncertainty of curved surfaces
         199             :    real(r8),parameter :: m_depcnt_bc = COS(theta_dep_bc*pi/180._r8)
         200             :    real(r8),parameter :: f_depcnt_bc = (2+m_depcnt_bc)*(1-m_depcnt_bc)**2/4._r8
         201             :
         202             :    real(r8),parameter :: m_depcnt_dst = COS(theta_dep_dust*pi/180._r8)
         203             :    real(r8),parameter :: f_depcnt_dust = (2+m_depcnt_dst)*(1-m_depcnt_dst)**2/4._r8
         204             :
         205             :    real(r8),parameter :: m_imm_bc = COS(theta_imm_bc*pi/180._r8)
         206             :    real(r8),parameter :: f_imm_bc = (2+m_imm_bc)*(1-m_imm_bc)**2/4._r8
         207             :
         208             :    real(r8),parameter :: m_imm_dust = COS(theta_imm_dust*pi/180._r8)
         209             :    real(r8),parameter :: f_imm_dust = (2+m_imm_dust)*(1-m_imm_dust)**2/4._r8
         210             :
         211             :    real(r8) :: f_dep, f_cnt, f_imm
         212             :    real(r8) :: dga_dep, dga_imm
         213             :    real(r8) :: limfac
         214             :    real(r8) :: frzimm, frzcnt, frzdep
         215             :    real(r8), pointer :: frzimm_ptr, frzcnt_ptr, frzdep_ptr
         216             :
         217             :    logical :: pdf_imm
         218             :
         219             :    integer :: ispc
         220             :
         221  1270492981 :    real(r8) :: ktherm(ntypes), kcoll(ntypes)
         222             :
         223             :    real(r8), parameter :: Ktherm_bc = 4.2_r8   ! black carbon thermal conductivity [J/(m s K)]
         224             :    real(r8), parameter :: Ktherm_dst = 0.72_r8 ! clay thermal conductivity [J/(m s K)]
         225             :
         226             :    !------------------------------------------------------------------------------------------------
         227             :
         228  1270492981 :    errstring = ' '
         229             :
         230  1270492981 :    nullify(frzimm_ptr)
         231  1270492981 :    nullify(frzcnt_ptr)
         232  1270492981 :    nullify(frzdep_ptr)
         233             :
         234  1270492981 :    frzbcimm = 0._r8
         235  1270492981 :    frzbccnt= 0._r8
         236  1270492981 :    frzbcdep = 0._r8
         237  1270492981 :    frzduimm = 0._r8
         238  1270492981 :    frzducnt= 0._r8
         239  1270492981 :    frzdudep = 0._r8
         240             :
         241             :    ! get saturation vapor pressure
         242  1270492981 :    eswtr = svp_water(t)  ! 0 for liquid
         243             :
         244  1270492981 :    tc = T - tmelt
         245  1270492981 :    rhoice = 916.7_r8-0.175_r8*tc-5.e-4_r8*tc**2
         246  1270492981 :    vwice = mwh2o*amu/rhoice
         247  1270492981 :    sigma_iw = (28.5_r8+0.25_r8*tc)*1E-3_r8
         248  1270492981 :    sigma_iv = (76.1_r8-0.155_r8*tc + 28.5_r8+0.25_r8*tc)*1E-3_r8
         249             :
         250             :    ! critical germ size
         251  1270492981 :    rgimm = 2*vwice*sigma_iw/(boltz*T*LOG(supersatice))
         252             :
         253             :    ! critical germ size
         254             :    ! assume 98% RH in mixed-phase clouds (Korolev & Isaac, JAS 2006)
         255  1270492981 :    rgdep=2*vwice*sigma_iv/(boltz*T*LOG(rhwincloud*supersatice))
         256             :
         257             :    ! homogeneous energy of germ formation
         258  1270492981 :    dg0dep = 4*pi/3._r8*sigma_iv*rgdep**2
         259             :
         260             :    ! prefactor
         261             :    ! attention: division of small numbers
         262  1270492981 :    Adep = (rhwincloud*eswtr)**2*(vwice/(mwh2o*amu))/(boltz*T*nus)*SQRT(sigma_iv/(boltz*T))
         263             :
         264             :    ! homogeneous energy of germ formation
         265  1270492981 :    dg0cnt = 4*pi/3._r8*sigma_iv*rgimm**2
         266             :
         267             :    ! prefactor
         268             :    ! attention: division of small numbers
         269  1270492981 :    Acnt = rhwincloud*eswtr*4*pi/(nus*SQRT(2*pi*mwh2o*amu*boltz*T))
         270             :
         271  6352464905 :    do ispc = 1, ntypes
         272             :
         273 11434436829 :       select case (trim(types(ispc)))
         274             :       case ('black-c')
         275  2540985962 :          ktherm(ispc) = ktherm_bc
         276             :       case ('dust')
         277  2540985962 :          ktherm(ispc) = ktherm_dst
         278             :       case default
         279           0 :          errstring = 'hetfrz_classnuc_calc ERROR: unrecognized aerosol type: '//trim(types(ispc))
         280 10163943848 :          return
         281             :       end select
         282             :    end do
         283             :
         284  1270492981 :    call collkernel(T, p, eswtr, rhwincloud, r3lx, hetraer, Ktherm, Kcoll)
         285             :
         286  6352464905 :    do ispc = 1, ntypes
         287             :
         288 10163943848 :       select case (trim(types(ispc)))
         289             :       case ('black-c')
         290  2540985962 :          f_dep = f_depcnt_bc
         291  2540985962 :          f_cnt = f_depcnt_bc
         292  2540985962 :          f_imm = f_imm_bc
         293  2540985962 :          dga_dep = dga_dep_bc
         294  2540985962 :          dga_imm = dga_imm_bc
         295  2540985962 :          pdf_imm = .false.
         296  2540985962 :          limfac = bc_limfac
         297  2540985962 :          frzimm_ptr => frzbcimm
         298  2540985962 :          frzcnt_ptr => frzbccnt
         299  2540985962 :          frzdep_ptr => frzbcdep
         300             :       case ('dust')
         301  2540985962 :          f_dep = f_depcnt_dust
         302  2540985962 :          f_cnt = f_depcnt_dust
         303  2540985962 :          f_imm = f_imm_dust
         304  2540985962 :          dga_dep = dga_dep_dust
         305  2540985962 :          dga_imm = dga_imm_dust
         306  2540985962 :          pdf_imm = .true.
         307  2540985962 :          limfac = dust_limfac
         308  2540985962 :          frzimm_ptr => frzduimm
         309  2540985962 :          frzcnt_ptr => frzducnt
         310  2540985962 :          frzdep_ptr => frzdudep
         311             :       case default
         312           0 :          errstring = 'hetfrz_classnuc_calc ERROR: unrecognized aerosol type: '//trim(types(ispc))
         313 10163943848 :          return
         314             :       end select
         315             :
         316             :       call hetfrz_classnuc_calc_rates( f_dep, f_cnt, f_imm, dga_dep, dga_imm, pdf_imm, limfac, &
         317             :            kcoll(ispc), hetraer(ispc), icnlx, r3lx, T, supersatice, sigma_iw, &
         318             :            rgimm, rgdep, dg0dep, Adep, dg0cnt, Acnt, vwice, deltat, &
         319             :            fn(ispc), wact_factor(ispc), dstcoat(ispc), &
         320             :            total_aer_num(ispc), total_interstitial_aer_num(ispc), total_cloudborne_aer_num(ispc), uncoated_aer_num(ispc), &
         321  5081971924 :            frzimm, frzcnt, frzdep, errstring )
         322             :
         323             :       ! accumulate dust and bc frz rates
         324  5081971924 :       frzimm_ptr = frzimm_ptr + frzimm
         325  5081971924 :       frzcnt_ptr = frzcnt_ptr + frzcnt
         326  6352464905 :       frzdep_ptr = frzdep_ptr + frzdep
         327             :
         328             :    end do
         329             :
         330  1270492981 :  end subroutine  hetfrz_classnuc_calc
         331             :
         332  5081971924 :  subroutine  hetfrz_classnuc_calc_rates( f_dep, f_cnt, f_imm, dga_dep, dga_imm, pdf_imm, limfac, &
         333             :       kcoll, mradius, icnlx, r3lx, T, supersatice, sigma_iw, &
         334             :       rgimm, rgdep, dg0dep, Adep, dg0cnt, Acnt, vwice, deltat, &
         335             :       fn, wact_factor, dstcoat, &
         336             :       total_aer_num, total_interstitial_aer_num, total_cloudborne_aer_num, uncoated_aer_num, &
         337             :       frzimm, frzcnt, frzdep, errstring )
         338             :
         339             :    ! input
         340             :    real(r8), intent(in) :: f_dep                      ! deposition form factor
         341             :    real(r8), intent(in) :: f_cnt                      ! contact form factor
         342             :    real(r8), intent(in) :: f_imm                      ! immersion form factor
         343             :    real(r8), intent(in) :: dga_dep                    ! deposition activation energy [J]
         344             :    real(r8), intent(in) :: dga_imm                    ! immersion activation energy [J]
         345             :    logical,  intent(in) :: pdf_imm                    ! PDF theta model switch (TRUE for dust)
         346             :    real(r8), intent(in) :: limfac                     ! Limit to 1% of available potential IN (for BC), no limit for dust
         347             :    real(r8), intent(in) :: kcoll                      ! collision kernel [cm3 s-1]
         348             :    real(r8), intent(in) :: mradius                    ! mass mean radius [m]
         349             :    real(r8), intent(in) :: icnlx                      ! in-cloud droplet concentration [cm-3]
         350             :    real(r8), intent(in) :: r3lx                       ! volume mean drop radius [m]
         351             :    real(r8), intent(in) :: T                          ! temperature [K]
         352             :    real(r8), intent(in) :: supersatice                ! supersaturation ratio wrt ice at 100%rh over water [ ]
         353             :    real(r8), intent(in) :: sigma_iw                   ! [J/m2]
         354             :    real(r8), intent(in) :: rgimm                      ! critical germ size
         355             :    real(r8), intent(in) :: rgdep                      ! critical germ size
         356             :    real(r8), intent(in) :: dg0dep                     ! homogeneous energy of germ formation
         357             :    real(r8), intent(in) :: Adep                       ! deposition nucleation prefactor
         358             :    real(r8), intent(in) :: dg0cnt                     ! homogeneous energy of germ formation
         359             :    real(r8), intent(in) :: Acnt                       ! contact nucleation prefactor
         360             :
         361             :    real(r8), intent(in) :: vwice
         362             :    real(r8), intent(in) :: deltat                     ! timestep [s]
         363             :    real(r8), intent(in) :: fn                         ! fraction activated [ ] for cloud borne aerosol number
         364             :    real(r8), intent(in) :: wact_factor                ! water activity factor -- density*(1.-(OC+BC)/(OC+BC+SO4)) [mug m-3]
         365             :    real(r8), intent(in) :: dstcoat                    ! coated fraction
         366             :    real(r8), intent(in) :: total_aer_num              ! total bc and dust number concentration(interstitial+cloudborne) [#/cm^3]
         367             :    real(r8), intent(in) :: total_interstitial_aer_num ! total bc and dust concentration(interstitial)
         368             :    real(r8), intent(in) :: total_cloudborne_aer_num   ! total bc and dust concentration(cloudborne)
         369             :    real(r8), intent(in) :: uncoated_aer_num           ! uncoated bc and dust number concentration(interstitial)
         370             :    ! output
         371             :    real(r8), intent(out) :: frzimm ! het. frz by immersion nucleation [cm-3 s-1]
         372             :    real(r8), intent(out) :: frzcnt ! het. frz by contact nucleation [cm-3 s-1]
         373             :    real(r8), intent(out) :: frzdep ! het. frz by deposition nucleation [cm-3 s-1]
         374             :
         375             :    character(len=*), intent(out) :: errstring
         376             :
         377             :    ! local vars
         378             :    real(r8) :: aw                          ! water activity [ ]
         379             :    real(r8) :: molal                       ! molality [moles/kg]
         380             :
         381             :    real(r8) :: Aimm
         382             :    real(r8) :: Jdep
         383             :    real(r8) :: Jimm
         384             :    real(r8) :: Jcnt
         385             :    real(r8) :: dg0imm
         386             :    real(r8) :: rgimm_aer
         387             :    real(r8) :: sum_imm
         388             :    real(r8) :: dim_Jimm(pdf_n_theta)
         389             :
         390             :    logical :: do_frz
         391             :    logical :: do_imm
         392             :
         393             :    integer :: i
         394             :
         395             :    !*****************************************************************************
         396             :    !                take water activity into account
         397             :    !*****************************************************************************
         398             :    !   solute effect
         399  5081971924 :    aw = 1._r8
         400  5081971924 :    molal = 0._r8
         401             :
         402             :    ! The heterogeneous ice freezing temperatures of all IN generally decrease with
         403             :    ! increasing total solute mole fraction. Therefore, the large solution concentration
         404             :    ! will cause the freezing point depression and the ice freezing temperatures of all
         405             :    ! IN will get close to the homogeneous ice freezing temperatures. Since we take into
         406             :    ! account water activity for three heterogeneous freezing modes(immersion, deposition,
         407             :    ! and contact), we utilize interstitial aerosols(not cloudborne aerosols) to calculate
         408             :    ! water activity.
         409             :    ! If the index of IN is 0, it means three freezing modes of this aerosol are depressed.
         410             :
         411             :    !calculate molality
         412  5081971924 :    if ( total_interstitial_aer_num > 0._r8 ) then
         413             :       molal = (1.e-6_r8*wact_factor/(mwso4*total_interstitial_aer_num*1.e6_r8))/ &
         414  5081971924 :            (4*pi/3*rhoh2o*(MAX(r3lx,4.e-6_r8))**3)
         415  5081971924 :       aw = 1._r8/(1._r8+2.9244948e-2_r8*molal+2.3141243e-3_r8*molal**2+7.8184854e-7_r8*molal**3)
         416             :    end if
         417             :
         418             :    !*****************************************************************************
         419             :    !                immersion freezing begin
         420             :    !*****************************************************************************
         421             :
         422  5081971924 :    frzimm = 0._r8
         423  5081971924 :    frzcnt = 0._r8
         424  5081971924 :    frzdep = 0._r8
         425             :
         426             :    ! take solute effect into account
         427  5081971924 :    rgimm_aer = rgimm
         428             :
         429             :    ! if aw*Si<=1, the freezing point depression is strong enough to prevent freezing
         430             :
         431  5081971924 :    do_frz = aw*supersatice > 1._r8
         432  5081971924 :    if (do_frz) then
         433  4800189048 :       rgimm_aer = 2*vwice*sigma_iw/(boltz*T*LOG(aw*supersatice))
         434             :    else
         435   281782876 :       return
         436             :    endif
         437             :
         438  4800189048 :    do_imm = T <= 263.15_r8 ! temperature threshold for immersion freezing (-10 C)
         439             :
         440  4800189048 :    if (do_imm) then
         441             :       ! homogeneous energy of germ formation
         442  3927906824 :       dg0imm = 4*pi/3._r8*sigma_iw*rgimm_aer**2
         443             :
         444             :       ! prefactor
         445  3927906824 :       Aimm = n1*((vwice*rhplanck)/(rgimm_aer**3)*SQRT(3._r8/pi*boltz*T*dg0imm))
         446             :
         447             :       ! nucleation rate per particle
         448             :
         449  3927906824 :       if (pdf_imm) then
         450  1834215884 :          dim_Jimm(:) = 0._r8
         451 >11372*10^7 :          do i = i1,i2
         452             :             ! 1/sqrt(f)
         453 >11188*10^7 :             dim_Jimm(i) = Aimm*mradius**2/SQRT(dim_f_imm(i))*EXP((-dga_imm-dim_f_imm(i)*dg0imm)/(boltz*T))
         454 >11372*10^7 :             dim_Jimm(i) = max(dim_Jimm(i), 0._r8)
         455             :          end do
         456             :
         457             :          sum_imm  = 0._r8
         458 >11188*10^7 :          do i = i1,i2-1
         459 >11005*10^7 :             sum_imm = sum_imm + 0.5_r8*((pdf_imm_theta(i  )*exp(-dim_Jimm(i  )*deltat)+ &
         460 >22194*10^7 :                  pdf_imm_theta(i+1)*exp(-dim_Jimm(i+1)*deltat)))*pdf_d_theta
         461             :          end do
         462  1834215884 :          if (sum_imm > 0.99_r8) then
         463   724871142 :             sum_imm = 1.0_r8
         464             :          end if
         465             :       else
         466  2093690940 :          Jimm = Aimm*mradius**2/SQRT(f_imm)*EXP(( -dga_imm - f_imm*dg0imm )/(boltz*T))
         467  2093690940 :          sum_imm = exp(-Jimm*deltat)
         468             :       end if
         469             :    end if
         470             :
         471             :    !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         472             :    !   Deposition nucleation
         473             :    !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         474             :
         475             :    ! nucleation rate per particle
         476  4800189048 :    if (rgdep > 0) then
         477  4800189048 :       Jdep = Adep*mradius**2/SQRT(f_dep)*EXP((-dga_dep-f_dep*dg0dep)/(boltz*T))
         478             :    else
         479             :       Jdep = 0._r8
         480             :    end if
         481             :
         482             :    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         483             :    ! contact nucleation
         484             :    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         485             :
         486             :    ! nucleation rate per particle
         487  4800189048 :    Jcnt = Acnt*mradius**2*EXP((-dga_dep-f_cnt*dg0cnt)/(boltz*T))*Kcoll*icnlx
         488             :
         489             :    ! Limit to 1% of available potential IN (for BC), no limit for dust
         490             :    if (tot_in) then
         491             :       if (do_imm) then
         492             :          frzimm = MIN(limfac*fn*total_aer_num/deltat, fn*total_aer_num/deltat*(1._r8-sum_imm))
         493             :       end if
         494             :       frzdep = MIN(limfac*(1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat, &
         495             :                    (1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat*(1._r8-exp(-Jdep*deltat)))
         496             :       frzcnt = MIN(limfac*(1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat, &
         497             :                    (1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat*(1._r8-exp(-Jcnt*deltat)))
         498             :    else
         499  4800189048 :       if (do_imm) then
         500  3927906824 :          frzimm = MIN(limfac*total_cloudborne_aer_num /deltat, total_cloudborne_aer_num/deltat*(1._r8-sum_imm))
         501             :       end if
         502  4800189048 :       frzdep = MIN(limfac*uncoated_aer_num/deltat, uncoated_aer_num/deltat*(1._r8-exp(-Jdep*deltat)))
         503  4800189048 :       frzcnt = MIN(limfac*uncoated_aer_num/deltat, uncoated_aer_num/deltat*(1._r8-exp(-Jcnt*deltat)))
         504             :    end if
         505             :
         506  4800189048 :    if (frzcnt <= -1._r8) then
         507           0 :       write(iulog,*) 'hetfrz_classnuc_calc: frzcnt, Jcnt, Kcoll: ', frzcnt, Jcnt, Kcoll
         508           0 :       errstring = 'ERROR in hetfrz_classnuc_calc::frzducnt'
         509           0 :       return
         510             :    end if
         511             :
         512  5081971924 : end subroutine  hetfrz_classnuc_calc_rates
         513             :
         514             : !===================================================================================================
         515             :
         516             : !-----------------------------------------------------------------------
         517             : !
         518             : ! Purpose: calculate collision kernels as a function of environmental parameters and aerosol/droplet sizes
         519             : !
         520             : ! Author: Corinna Hoose, UiO, October 2009
         521             : !
         522             : ! Modifications: Yong Wang and Xiaohong Liu, UWyo, 12/2012
         523             : !
         524             : ! "Seinfeld & Pandis" referenced in several places in this routine is:
         525             : ! Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 3rd Edition
         526             : ! John H. Seinfeld, Spyros N. Pandis  ISBN: 978-1-118-94740-1
         527             : !-----------------------------------------------------------------------
         528             :
         529  1270492981 : subroutine collkernel( temp, pres, eswtr, rhwincloud, r3lx,  rad, Ktherm, Kcoll )
         530             :
         531             :    real(r8), intent(in) :: temp       ! temperature [K]
         532             :    real(r8), intent(in) :: pres       ! pressure [Pa]
         533             :    real(r8), intent(in) :: eswtr      ! saturation vapor pressure of water [Pa]
         534             :    real(r8), intent(in) :: r3lx       ! volume mean drop radius [m]
         535             :    real(r8), intent(in) :: rhwincloud ! in-cloud relative humidity over water [ ]
         536             :    real(r8), intent(in) :: rad(:)     ! aerosol radius [m]
         537             :    real(r8), intent(in) :: Ktherm(:)  ! thermal conductivity of aerosol [J/(m s K)]
         538             :    real(r8), intent(out) :: Kcoll(:)  ! collision kernel [cm3 s-1]
         539             :
         540             :    ! local variables
         541             :    real(r8) :: a, b, c, a_f, b_f, c_f, f
         542             :    real(r8) :: tc          ! temperature [deg C]
         543             :    real(r8) :: rho_air     ! air density [kg m-3]
         544             :    real(r8) :: viscos_air  ! dynamic viscosity of air [kg m-1 s-1]
         545             :    real(r8) :: Ktherm_air  ! thermal conductivity of air [J/(m s K)]
         546             :    real(r8) :: lambda      ! mean free path [m]
         547             :    real(r8) :: Kn          ! Knudsen number [ ]
         548             :    real(r8) :: Re          ! Reynolds number [ ]
         549             :    real(r8) :: Pr          ! Prandtl number [ ]
         550             :    real(r8) :: Sc          ! Schmidt number [ ]
         551             :    real(r8) :: vterm       ! terminal velocity [m s-1]
         552             :    real(r8) :: Dvap        ! water vapor diffusivity [m2 s-1]
         553             :    real(r8) :: Daer        ! aerosol diffusivity [m2 s-1]
         554             :    real(r8) :: latvap      ! latent heat of vaporization [J kg-1]
         555             :    real(r8) :: G           ! thermodynamic function in Cotton et al. [kg m-1 s-1]
         556             :    real(r8) :: f_t         ! factor by Waldmann & Schmidt [ ]
         557             :    real(r8) :: Q_heat      ! heat flux [J m-2 s-1]
         558             :    real(r8) :: Tdiff_cotton ! temperature difference between droplet and environment [K]
         559             :    real(r8) :: K_brownian,K_thermo_cotton,K_diffusio_cotton   ! collision kernels [m3 s-1]
         560             :
         561             :    integer :: ntot, idx
         562             :
         563             :    !------------------------------------------------------------------------------------------------
         564             :
         565  1270492981 :    ntot = size(ktherm)
         566             :
         567  6352464905 :    Kcoll(:) = 0._r8
         568             :
         569  1270492981 :    tc = temp - tmelt
         570             :
         571             :    ! air viscosity for tc<0, from depvel_part.F90
         572  1270492981 :    viscos_air = (1.718_r8+0.0049_r8*tc-1.2e-5_r8*tc*tc)*1.e-5_r8
         573             :    ! air density
         574  1270492981 :    rho_air = pres/(rair*temp)
         575             :    ! mean free path: Seinfeld & Pandis 8.6 (Book: ISBN: 978-1-118-94740-1)
         576  1270492981 :    lambda = 2*viscos_air/(pres*SQRT(8/(pi*rair*temp)))
         577             :    ! latent heat of vaporization, varies with T
         578  1270492981 :    latvap = 1000*(-0.0000614342_r8*tc**3 + 0.00158927_r8*tc**2 - 2.36418_r8*tc + 2500.79_r8)
         579             :    ! droplet terminal velocity after Chen & Liu, QJRMS 2004 (https://doi-org.cuucar.idm.oclc.org/10.1256/qj.03.41)
         580  1270492981 :    a = 8.8462e2_r8
         581  1270492981 :    b = 9.7593e7_r8
         582  1270492981 :    c = -3.4249e-11_r8
         583  1270492981 :    a_f = 3.1250e-1_r8
         584  1270492981 :    b_f = 1.0552e-3_r8
         585  1270492981 :    c_f = -2.4023_r8
         586  1270492981 :    f = EXP(EXP(a_f + b_f*(LOG(r3lx))**3 + c_f*rho_air**1.5_r8))
         587  1270492981 :    vterm = (a+ (b + c*r3lx)*r3lx)*r3lx*f
         588             :
         589             :    ! Reynolds number
         590  1270492981 :    Re = 2*vterm*r3lx*rho_air/viscos_air
         591             :    ! thermal conductivity of air: Seinfeld & Pandis eq. 15.75 (Book: ISBN: 978-1-118-94740-1)
         592  1270492981 :    Ktherm_air = 1.e-3_r8*(4.39_r8+0.071_r8*temp)  !J/(m s K)
         593             :    ! Prandtl number
         594  1270492981 :    Pr = viscos_air*cpair/Ktherm_air
         595             :    ! water vapor diffusivity: Pruppacher & Klett 13-3 (https://link.springer.com/book/10.1007/978-0-306-48100-0)
         596  1270492981 :    Dvap = 0.211e-4_r8*(temp/tmelt)*(pstd/pres)
         597             :    ! G-factor = rhoh2o*Xi in Rogers & Yau, p. 104
         598             :    G = rhoh2o/((latvap/(rh2o*temp) - 1)*latvap*rhoh2o/(Ktherm_air*temp) &
         599  1270492981 :        + rhoh2o*rh2o*temp/(Dvap*eswtr))
         600             :
         601  6352464905 :    do idx = 1,ntot
         602  6352464905 :       if (rad(idx)>0._r8) then
         603             :          ! Knudsen number (Seinfeld & Pandis 8.1) (Book: ISBN: 978-1-118-94740-1)
         604  5081971924 :          Kn = lambda/rad(idx)
         605             :          ! aerosol diffusivity
         606  5081971924 :          Daer = boltz*temp*(1 + Kn)/(6*pi*rad(idx)*viscos_air)
         607             :
         608             :          ! Schmidt number
         609  5081971924 :          Sc = viscos_air/(Daer*rho_air)
         610             :
         611             :          ! Young (1974) first equ. on page 771 (doi: 10.1175/1520-0469(1974)031<0768:TROCNI>2.0.CO;2)
         612  5081971924 :          K_brownian = 4*pi*r3lx*Daer*(1 + 0.3_r8*Re**0.5_r8*Sc**0.33_r8)
         613             :
         614             :          ! thermal conductivities from Seinfeld & Pandis, Table 8.6 (Book: ISBN: 978-1-118-94740-1)
         615             :          ! form factor
         616             :          f_t = 0.4_r8*(1._r8 + 1.45_r8*Kn + 0.4_r8*Kn*EXP(-1._r8/Kn))      &
         617  5081971924 :               *(Ktherm_air + 2.5_r8*Kn*Ktherm(idx))                      &
         618  5081971924 :               /((1._r8 + 3._r8*Kn)*(2._r8*Ktherm_air + 5._r8*Kn*Ktherm(idx)+Ktherm(idx)))
         619             :
         620             :          ! calculate T-Tc as in Cotton et al.
         621  5081971924 :          Tdiff_cotton = -G*(rhwincloud - 1._r8)*latvap/Ktherm_air
         622  5081971924 :          Q_heat = Ktherm_air/r3lx*(1._r8 + 0.3_r8*Re**0.5_r8*Pr**0.33_r8)*Tdiff_cotton
         623  5081971924 :          K_thermo_cotton = 4._r8*pi*r3lx*r3lx*f_t*Q_heat/pres
         624  5081971924 :          K_diffusio_cotton = -(1._r8/f_t)*(rh2o*temp/latvap)*K_thermo_cotton
         625  5081971924 :          Kcoll(idx) = 1.e6_r8*(K_brownian + K_thermo_cotton + K_diffusio_cotton)  ! convert m3/s -> cm3/s
         626             :
         627             :          ! set K to 0 if negative
         628  5081971924 :          if (Kcoll(idx) < 0._r8) Kcoll(idx) = 0._r8
         629             :       else
         630           0 :          Kcoll(idx) = 0._r8
         631             :       endif
         632             :    end do
         633             :
         634  1270492981 : end subroutine collkernel
         635             :
         636             : !===================================================================================================
         637             :
         638        1536 : subroutine hetfrz_classnuc_init_pdftheta()
         639             :
         640             :    ! Local variables:
         641             :    real(r8) :: theta_min, theta_max
         642             :    real(r8) :: x1_imm, x2_imm
         643             :    real(r8) :: norm_theta_imm
         644             :    real(r8) :: imm_dust_mean_theta
         645             :    real(r8) :: imm_dust_var_theta
         646             :    integer  :: i
         647             :    real(r8) :: m
         648             :    real(r8) :: temp
         649             :    !----------------------------------------------------------------------------
         650             :
         651        1536 :    theta_min           = pi/180._r8
         652        1536 :    theta_max           = 179._r8/180._r8*pi
         653        1536 :    imm_dust_mean_theta = 46.0_r8/180.0_r8*pi
         654        1536 :    imm_dust_var_theta  = 0.01_r8
         655             :
         656        1536 :    pdf_d_theta = (179._r8-1._r8)/180._r8*pi/(pdf_n_theta-1)
         657             :
         658        1536 :    x1_imm = (LOG(theta_min) - LOG(imm_dust_mean_theta))/(sqrt(2.0_r8)*imm_dust_var_theta)
         659        1536 :    x2_imm = (LOG(theta_max) - LOG(imm_dust_mean_theta))/(sqrt(2.0_r8)*imm_dust_var_theta)
         660        1536 :    norm_theta_imm = (ERF(x2_imm) - ERF(x1_imm))*0.5_r8
         661        1536 :    dim_theta      = 0.0_r8
         662        1536 :    pdf_imm_theta  = 0.0_r8
         663       95232 :    do i = i1, i2
         664       93696 :       dim_theta(i)     = 1._r8/180._r8*pi + (i-1)*pdf_d_theta
         665             :       pdf_imm_theta(i) = exp(-((LOG(dim_theta(i)) - LOG(imm_dust_mean_theta))**2._r8) / &
         666             :                               (2._r8*imm_dust_var_theta**2._r8) ) /                     &
         667       95232 :                          (dim_theta(i)*imm_dust_var_theta*SQRT(2*pi))/norm_theta_imm
         668             :    end do
         669             :
         670       95232 :    do i = i1, i2
         671       93696 :      m = cos(dim_theta(i))
         672       93696 :      temp = (2+m)*(1-m)**2/4._r8
         673       95232 :      dim_f_imm(i) = temp
         674             :    end do
         675             :
         676        1536 : end subroutine hetfrz_classnuc_init_pdftheta
         677             :
         678             : !===================================================================================================
         679             :
         680             : end module hetfrz_classnuc

Generated by: LCOV version 1.14