Line data    Source code

       1             : ! Include shortname defintions, so that the F77 code does not have to be modified to
       2             : ! reference the CARMA structure.
       3             : #include "carma_globaer.h"
       4             : 
       5             : !!  Calculates particle production rates due to heterogeneous
       6             : !!  nucleation <rnuclg>:
       7             : !!
       8             : !!  This was moved from sulfnuc to make the code more manageable.
       9             : !!
      10             : !!  @author Mike Mills, Chuck Bardeen
      11             : !!  @version Jun-2013
      12           0 : subroutine sulfhetnucrate(carma, cstate, iz, igroup, nucbin, h2o, h2so4, beta1, beta2, nucrate, rc)
      13             :   use carma_precision_mod
      14             :   use carma_enums_mod
      15             :   use carma_constants_mod
      16             :   use carma_types_mod
      17             :   use carmastate_mod
      18             :   use carma_mod
      19             :   use sulfate_utils
      20             : 
      21             :   implicit none
      22             : 
      23             :   type(carma_type), intent(in)         :: carma       !! the carma object
      24             :   type(carmastate_type), intent(inout) :: cstate      !! the carma state object
      25             :   integer, intent(in)                  :: iz          !! level index
      26             :   integer, intent(in)                  :: igroup      !! group index
      27             :   integer, intent(in)                  :: nucbin      !! bin in which nucleation occurs
      28             :   real(kind=f), intent(in)             :: h2o         !! H2O concentrations in molec/cm3
      29             :   real(kind=f), intent(in)             :: h2so4       !! H2SO4 concentrations in molec/cm3
      30             :   real(kind=f), intent(in)             :: beta1
      31             :   real(kind=f), intent(in)             :: beta2
      32             :   real(kind=f), intent(out)            :: nucrate     !! nucleation rate #/x/y/z/s
      33             :   integer, intent(inout)               :: rc          !! return code, negative indicates failure
      34             : 
      35             :   !  Local declarations
      36             :   real(kind=f)                         :: cnucl
      37             :   real(kind=f)                         :: chom
      38             :   real(kind=f)                         :: expc
      39             :   real(kind=f)                         :: chet
      40             :   real(kind=f)                         :: xm
      41             :   real(kind=f)                         :: xm1
      42             :   real(kind=f)                         :: fxm
      43             :   real(kind=f)                         :: fv2
      44             :   real(kind=f)                         :: fu2
      45             :   real(kind=f)                         :: fv3
      46             :   real(kind=f)                         :: fv4
      47             :   real(kind=f)                         :: v1
      48             :   real(kind=f)                         :: fv1
      49             :   real(kind=f)                         :: ftry
      50             :   real(kind=f)                         :: ftry1
      51             :   real(kind=f)                         :: rarea
      52             :   real(kind=f)                         :: gg
      53             :   real(kind=f)                         :: FM = cos(50._f * DEG2RAD)  ! cos(contact angle)
      54             :   real(kind=f)                         :: h2so4_cgs        ! H2SO4 densities in g/cm3
      55             :   real(kind=f)                         :: h2o_cgs          ! H2O densities in g/cm3
      56             :   real(kind=f)                         :: mass_cluster_dry ! dry mass of the cluster
      57             :   real(kind=f)                         :: nucrate_cgs      ! binary nucleation rate, j (# cm-3 s-1)
      58             :   real(kind=f)                         :: rh               ! relative humidity (0-1)
      59             :   real(kind=f)                         :: rstar       !! critical radius (cm)
      60             : 
      61             :   ! Zhao heterogeneous nucleation rate depends on calculations for ftry and rstar made in Zhao homogeneous nucleation.
      62             : 
      63             :   ! Compute H2SO4 densities in g/cm3
      64           0 :   h2so4_cgs = gc(iz, igash2so4) / zmet(iz)
      65             : 
      66             :   ! Compute H2O densities in g/cm3
      67           0 :   h2o_cgs = gc(iz, igash2o) / zmet(iz)
      68             : 
      69             :   ! Compute relative humidity of water wrt liquid water
      70           0 :   rh = (supsatl(iz, igash2o) + 1._f)
      71             : 
      72           0 :   call binary_nuc_zhao1995( carma, cstate, t(iz), wtpct(iz), rh, h2so4, h2so4_cgs, h2o, h2o_cgs, beta1, &
      73           0 :                 nucrate_cgs, mass_cluster_dry, rstar, ftry, rc )
      74             : 
      75           0 :   if (rstar > 0._f) then
      76             :     ! Heterogeneous nucleation which depends on r
      77           0 :     cnucl = 4._f * PI * rstar**(2._f)
      78             :     chom  = h2so4 * h2o * beta1 * cnucl
      79           0 :     expc  = 2.4e-16_f * exp(4.51872e+11_f / RGAS / t(iz))
      80           0 :     chet  = chom * expc * beta2
      81           0 : 
      82             :     xm    = r(nucbin, igroup) / rstar
      83           0 : 
      84             :     if (xm .lt. 1._f) then
      85           0 :       fxm = sqrt(1._f - 2._f * FM * xm + xm**(2._f))
      86           0 :       fv2 = (xm - FM) / fxm
      87           0 :       fu2 = (1._f - xm * FM) / fxm
      88           0 :       fv3 = (2._f + fv2) * xm**3._f  * (fv2 - 1._f)**(2._f)
      89           0 :       fv4 = 3._f * FM    * xm**2._f  * (fv2 - 1._f)
      90           0 :     else
      91             :       xm1 = 1._f / xm
      92           0 :       fxm = sqrt(1._f - 2._f * FM * xm1 + xm1**2._f)
      93           0 :       fu2 = (xm1 - FM) / fxm
      94           0 :       fv2 = (1._f - xm1 * FM) / fxm
      95           0 :       v1  = (FM**(2._f) - 1._f) / (fv2 + 1._f) / fxm**(2._f)
      96           0 :       fv3 = (2._f + fv2) * xm1 * v1**2._f
      97           0 :       fv4 = 3._f * FM * v1
      98           0 :     endif
      99             : 
     100             :     fv1   = 0.5_f * (1._f + fu2**3._f + fv3 + fv4)
     101           0 : 
     102             :     ftry1 = ftry * fv1
     103           0 :     if (ftry1 .lt. -1000._f) then
     104           0 :       nucrate = 0._f
     105           0 :     else
     106             : 
     107             :       rarea = 4._f * PI * r(nucbin, igroup)**2._f ! surface area per nucleus
     108           0 :       gg    = exp(ftry1)
     109           0 : 
     110             :       ! Calculate heterogeneous nucleation rate [embryos/s]
     111             :       ! NOTE: for [embryos/gridpoint/s], multipy rnuclg by pc [nuclei/gridpoint]
     112             :       nucrate = chet * gg * rarea  ! embryos/s
     113           0 :     end if
     114             :   end if
     115             : 
     116             :   return
     117           0 : end subroutine sulfhetnucrate