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             : !! This routine evaluates particle fall velocities, vf(k) [cm s^-1]
       6             : !! and reynolds' numbers based on fall velocities, re(j,i,k) [dimensionless].
       7             : !! indices correspond to vertical level <k>, bin index <i>, and aerosol
       8             : !! group <j>.
       9             : !!
      10             : !! Method: Use the routined from Heymsfield and Westbrook [2010], which is
      11             : !! designed only for ice particles. Thus this routine uses the dry mass and
      12             : !! radius, not the wet mass and radius. The area ration (Ar) is determined
      13             : !! based upon the formulation of Schmitt and Heymsfield [JAS, 2009].
      14             : !!
      15             : !! @author  Chuck Bardeen
      16             : !! @version Mar-2010
      17           0 : subroutine setupvf_heymsfield2010(carma, cstate, j, rc)
      18             : 
      19             :   ! types
      20             :   use carma_precision_mod
      21             :   use carma_enums_mod
      22             :   use carma_constants_mod
      23             :   use carma_types_mod
      24             :   use carmastate_mod
      25             :   use carma_mod
      26             : 
      27             :   implicit none
      28             : 
      29             :   type(carma_type), intent(in)         :: carma   !! the carma object
      30             :   type(carmastate_type), intent(inout) :: cstate  !! the carma state object
      31             :   integer, intent(in)                  :: j       !! group index
      32             :   integer, intent(inout)               :: rc      !! return code, negative indicates failure
      33             : 
      34             :   ! Local declarations
      35             :   integer                  :: i, k
      36             :   real(kind=f)             :: rhoa_cgs, vg, rmfp, rkn, expon, x
      37             :   real(kind=f), parameter  :: c0      = 0.35_f
      38             :   real(kind=f), parameter  :: delta0  = 8.0_f
      39             : 
      40             :   real(kind=f)             :: dmax                ! maximum diameter
      41             : 
      42             : 
      43             :   ! Loop over all atltitudes.
      44           0 :   do k = 1, NZ
      45             : 
      46             :     ! This is <rhoa> in cartesian coordinates (good old cgs units)
      47           0 :     rhoa_cgs = rhoa(k) / zmet(k)
      48             : 
      49             :     ! <vg> is mean thermal velocity of air molecules [cm/s]
      50           0 :     vg = sqrt(8._f / PI * R_AIR * t(k))
      51             : 
      52             :     ! <rmfp> is mean free path of air molecules [cm]
      53           0 :     rmfp = 2._f * rmu(k) / (rhoa_cgs * vg)
      54             : 
      55             :     ! Loop over particle size bins.
      56           0 :     do i = 1,NBIN
      57             : 
      58             :       ! <rkn> is knudsen number
      59             : !      rkn = rmfp / r(i,j)
      60           0 :       rkn = rmfp / (r_wet(k,i,j) * rrat(i,j))
      61             : 
      62             :       ! <bpm> is the slip correction factor, the correction term for
      63             :       ! non-continuum effects.  Also used to calculate coagulation kernels
      64             :       ! and diffusion coefficients.
      65           0 :       expon = -.87_f / rkn
      66           0 :       expon = max(-POWMAX, expon)
      67           0 :       bpm(k,i,j) = 1._f + (1.246_f*rkn + 0.42_f*rkn*exp(expon))
      68             : 
      69           0 :       dmax = 2._f * r_wet(k,i,j) * rrat(i,j)
      70             : 
      71           0 :       x = (rhoa_cgs / (rmu(k)**2)) * &
      72           0 :            ((8._f * rmass(i,j) * GRAV) / (PI * (arat(i,j)**0.5_f)))
      73             : 
      74             :       ! Apply the slip correction factor. This is not included in the formulation
      75             :       ! from Heymsfield and Westbrook [2010].
      76             :       !
      77             :       ! NOTE: This is applied according to eq 8.46 and surrounding discussion in
      78             :       ! Seinfeld and Pandis [1998].
      79           0 :       x = x * bpm(k,i,j)
      80             : 
      81           0 :       re(k,i,j) = ((delta0**2) / 4._f) * (sqrt(1._f + (4._f * sqrt(x) / (delta0**2 * sqrt(c0)))) - 1._f)**2
      82             : 
      83             : 
      84           0 :       vf(k,i,j) = rmu(k) * re(k,i,j) / (rhoa_cgs * dmax)
      85             :     enddo      ! <i=1,NBIN>
      86             :   enddo      ! <k=1,NZ>
      87             : 
      88             :   ! Return to caller with particle fall velocities evaluated.
      89           0 :   return
      90           0 : end