Line data    Source code

       1             : !=============================================================================
       2             : ! Common dust module
       3             : !=============================================================================
       4             : module dust_common
       5             : 
       6             :   use shr_kind_mod,     only: r8 => shr_kind_r8, cl => shr_kind_cl
       7             :   use cam_abortutils,   only: endrun
       8             :   use cam_logfile,      only: iulog
       9             : 
      10             :   implicit none
      11             :   private
      12             : 
      13             :   public :: dust_set_params
      14             : 
      15             : contains
      16             : 
      17             :   !=============================================================================
      18             :   !
      19             :   ! !DESCRIPTION: 
      20             :   !
      21             :   ! Compute source efficiency factor from topography
      22             :   ! Initialize other variables used in subroutine Dust:
      23             :   ! ovr_src_snk_mss(m,n) and tmp1.
      24             :   ! Define particle diameter and density needed by atm model
      25             :   ! as well as by dry dep model
      26             :   ! Source: Paul Ginoux (for source efficiency factor)
      27             :   ! Modifications by C. Zender and later by S. Levis
      28             :   ! Rest of subroutine from C. Zender's dust model
      29             :   !=============================================================================
      30           0 :   subroutine dust_set_params( nbin, dmt_grd, dmt_vwr, stk_crc )
      31             : 
      32             :     !
      33             :     ! !USES
      34             :     !
      35             :     use physconst,     only: pi,rair, gravit
      36             :     use mo_constants,  only: dust_density
      37             :     use infnan,        only: nan, assignment(=)
      38             : 
      39             :     !
      40             :     ! !ARGUMENTS:
      41             :     !
      42             :     integer, intent(in)  :: nbin
      43             :     real(r8),intent(in)  :: dmt_grd(:)
      44             :     real(r8),intent(out) :: dmt_vwr(:)
      45             :     real(r8),intent(out) :: stk_crc(:)
      46             : 
      47             :     !
      48             :     ! !REVISION HISTORY
      49             :     ! Created by Samual Levis
      50             :     ! Revised for CAM by Natalie Mahowald
      51             :     !EOP
      52             :     !------------------------------------------------------------------------
      53             : 
      54             :     !------------------------------------------------------------------------
      55             :     !Local Variables
      56             :     integer, parameter:: dst_src_nbr =3
      57             :     integer, parameter:: sz_nbr =200
      58             : 
      59             :     integer  :: m,n                     !indices
      60           0 :     real(r8) :: dmt_min(nbin)      ![m] Size grid minimum
      61           0 :     real(r8) :: dmt_max(nbin)      ![m] Size grid maximum
      62           0 :     real(r8) :: dmt_ctr(nbin)      ![m] Diameter at bin center
      63           0 :     real(r8) :: dmt_dlt(nbin)      ![m] Width of size bin
      64           0 :     real(r8) :: slp_crc(nbin)      ![frc] Slip correction factor
      65           0 :     real(r8) :: vlm_rsl(nbin)      ![m3 m-3] Volume concentration resolved
      66           0 :     real(r8) :: vlc_stk(nbin)      ![m s-1] Stokes settling velocity
      67           0 :     real(r8) :: vlc_grv(nbin)      ![m s-1] Settling velocity
      68           0 :     real(r8) :: ryn_nbr_grv(nbin)  ![frc] Reynolds number at terminal velocity
      69           0 :     real(r8) :: cff_drg_grv(nbin)  ![frc] Drag coefficient at terminal velocity
      70             :     real(r8) :: tmp                     !temporary 
      71             :     real(r8) :: ln_gsd                  ![frc] ln(gsd)
      72             :     real(r8) :: gsd_anl                 ![frc] Geometric standard deviation
      73             :     real(r8) :: dmt_vma                 ![m] Mass median diameter analytic She84 p.75 Tabl.1
      74             :     real(r8) :: dmt_nma                 ![m] Number median particle diameter
      75             :     real(r8) :: lgn_dst                 !Lognormal distribution at sz_ctr
      76             :     real(r8) :: eps_max                 ![frc] Relative accuracy for convergence
      77             :     real(r8) :: eps_crr                 ![frc] Current relative accuracy
      78             :     real(r8) :: itr_idx                 ![idx] Counting index
      79             :     real(r8) :: dns_mdp                 ![kg m-3] Midlayer density
      80             :     real(r8) :: mfp_atm                 ![m] Mean free path of air
      81             :     real(r8) :: vsc_dyn_atm             ![kg m-1 s-1] Dynamic viscosity of air
      82             :     real(r8) :: vsc_knm_atm             ![kg m-1 s-1] Kinematic viscosity of air
      83             :     real(r8) :: vlc_grv_old             ![m s-1] Previous gravitational settling velocity
      84             :     real(r8) :: series_ratio            !Factor for logarithmic grid
      85             :     real(r8) :: lngsdsqrttwopi_rcp      !Factor in lognormal distribution
      86             :     real(r8) :: sz_min(sz_nbr)          ![m] Size Bin minima
      87             :     real(r8) :: sz_max(sz_nbr)          ![m] Size Bin maxima
      88             :     real(r8) :: sz_ctr(sz_nbr)          ![m] Size Bin centers
      89             :     real(r8) :: sz_dlt(sz_nbr)          ![m] Size Bin widths
      90             : 
      91           0 :     stk_crc(:) = nan
      92           0 :     dmt_vwr(:) = nan
      93             : 
      94             :     ! Introducing particle diameter. Needed by atm model and by dry dep model.
      95             :     ! Taken from Charlie Zender's subroutines dst_psd_ini, dst_sz_rsl,
      96             :     ! grd_mk (dstpsd.F90) and subroutine lgn_evl (psdlgn.F90)
      97             : 
      98             :     ! Charlie allows logarithmic or linear option for size distribution
      99             :     ! however, he hardwires the distribution to logarithmic in his code
     100             :     ! therefore, I take his logarithmic code only
     101             :     ! furthermore, if dst_nbr == 4, he overrides the automatic grid calculation
     102             :     ! he currently works with dst_nbr = 4, so I only take the relevant code
     103             :     ! if dust_number ever becomes different from 4, must add call grd_mk (dstpsd.F90)
     104             :     ! as done in subroutine dst_psd_ini
     105             :     ! note that here dust_number = dst_nbr
     106             : 
     107             :     ! Override automatic grid with preset grid if available
     108           0 :     do n = 1, nbin
     109           0 :        dmt_min(n) = dmt_grd(n)                       ![m] Max diameter in bin
     110           0 :        dmt_max(n) = dmt_grd(n+1)                     ![m] Min diameter in bin
     111           0 :        dmt_ctr(n) = 0.5_r8 * (dmt_min(n)+dmt_max(n)) ![m] Diameter at bin ctr
     112           0 :        dmt_dlt(n) = dmt_max(n)-dmt_min(n)            ![m] Width of size bin
     113             :     end do
     114             : 
     115             :  ! sets dust_dmt_vwr ....
     116             : 
     117             :     ! Bin physical properties
     118             :     gsd_anl = 2.0_r8      ! [frc] Geometric std dev PaG77 p. 2080 Table1
     119             :     ln_gsd = log(gsd_anl)
     120             : 
     121             :     ! Set a fundamental statistic for each bin
     122             :     dmt_vma = 2.524e-6_r8 ! [m] Mass median diameter analytic She84 p.75 Table1
     123             :     dmt_vma = 3.5e-6_r8
     124             :     ! Compute analytic size statistics
     125             :     ! Convert mass median diameter to number median diameter (call vma2nma)
     126             :     dmt_nma = dmt_vma * exp(-3.0_r8*ln_gsd*ln_gsd) ! [m]
     127             :     ! Compute resolved size statistics for each size distribution
     128             :     ! In C. Zender's code call dst_sz_rsl
     129           0 :     do n = 1, nbin
     130           0 :        series_ratio = (dmt_max(n)/dmt_min(n))**(1.0_r8/sz_nbr)
     131           0 :        sz_min(1) = dmt_min(n)
     132           0 :        do m = 2, sz_nbr                            ! Loop starts at 2
     133           0 :           sz_min(m) = sz_min(m-1) * series_ratio
     134             :        end do
     135             : 
     136             :        ! Derived grid values
     137           0 :        do m = 1, sz_nbr-1                          ! Loop ends at sz_nbr-1
     138           0 :           sz_max(m) = sz_min(m+1)                  ! [m]
     139             :        end do
     140           0 :        sz_max(sz_nbr) = dmt_max(n)                 ! [m]
     141             : 
     142             :        ! Final derived grid values
     143           0 :        do m = 1, sz_nbr
     144           0 :           sz_ctr(m) = 0.5_r8 * (sz_min(m)+sz_max(m))
     145           0 :           sz_dlt(m) = sz_max(m)-sz_min(m)
     146             :        end do
     147           0 :        lngsdsqrttwopi_rcp = 1.0_r8 / (ln_gsd*sqrt(2.0_r8*pi))
     148           0 :        dmt_vwr(n) = 0.0_r8 ! [m] Mass wgted diameter resolved
     149           0 :        vlm_rsl(n) = 0.0_r8 ! [m3 m-3] Volume concentration resolved
     150           0 :        do m = 1, sz_nbr
     151             :           ! Evaluate lognormal distribution for these sizes (call lgn_evl)
     152           0 :           tmp = log(sz_ctr(m)/dmt_nma) / ln_gsd
     153           0 :           lgn_dst = lngsdsqrttwopi_rcp * exp(-0.5_r8*tmp*tmp) / sz_ctr(m)
     154             :           ! Integrate moments of size distribution
     155             :           dmt_vwr(n) = dmt_vwr(n) + sz_ctr(m) *                    &
     156             :                pi / 6.0_r8 * (sz_ctr(m)**3.0_r8) * & ![m3] Volume
     157           0 :                lgn_dst * sz_dlt(m)                ![# m-3] Number concentrn
     158             :           vlm_rsl(n) = vlm_rsl(n) +                                &
     159             :                pi / 6.0_r8 * (sz_ctr(m)**3.0_r8) * & ![m3] Volume
     160           0 :                lgn_dst * sz_dlt(m)                ![# m-3] Number concentrn
     161             :        end do
     162           0 :        dmt_vwr(n) = dmt_vwr(n) / vlm_rsl(n) ![m] Mass weighted diameter resolved
     163             :     end do
     164             : 
     165             :   ! sets stk_crc ...
     166             : 
     167             :     ! calculate correction to Stokes' settling velocity (subroutine stk_crc_get)
     168           0 :     eps_max = 1.0e-4_r8
     169           0 :     dns_mdp = 100000._r8 / (295.0_r8*rair) ![kg m-3] const prs_mdp & tpt_vrt
     170             :     ! Size-independent thermokinetic properties
     171             :     vsc_dyn_atm = 1.72e-5_r8 * ((295.0_r8/273.0_r8)**1.5_r8) * 393.0_r8 / &
     172           0 :          (295.0_r8+120.0_r8)      ![kg m-1 s-1] RoY94 p.102 tpt_mdp=295.0
     173             :     mfp_atm = 2.0_r8 * vsc_dyn_atm / &  !SeP97 p. 455 constant prs_mdp, tpt_mdp
     174           0 :          (100000._r8*sqrt(8.0_r8/(pi*rair*295.0_r8)))
     175           0 :     vsc_knm_atm = vsc_dyn_atm / dns_mdp ![m2 s-1] Kinematic viscosity of air
     176             : 
     177           0 :     do m = 1, nbin
     178           0 :        slp_crc(m) = 1.0_r8 + 2.0_r8 * mfp_atm *                      &
     179           0 :             (1.257_r8+0.4_r8*exp(-1.1_r8*dmt_vwr(m)/(2.0_r8*mfp_atm))) / &
     180           0 :             dmt_vwr(m)                      ! [frc] Slip correction factor SeP97 p.464
     181             :        vlc_stk(m) = (1.0_r8/18.0_r8) * dmt_vwr(m) * dmt_vwr(m) * dust_density * &
     182           0 :             gravit * slp_crc(m) / vsc_dyn_atm ! [m s-1] SeP97 p.466
     183             :     end do
     184             : 
     185             :     ! For Reynolds number flows Re < 0.1 Stokes' velocity is valid for
     186             :     ! vlc_grv SeP97 p. 466 (8.42). For larger Re, inertial effects become
     187             :     ! important and empirical drag coefficients must be employed
     188             :     ! Implicit equation for Re, Cd, and Vt is SeP97 p. 467 (8.44)
     189             :     ! Using Stokes' velocity rather than iterative solution with empirical
     190             :     ! drag coefficient causes 60% errors for D = 200 um SeP97 p. 468
     191             : 
     192             :     ! Iterative solution for drag coefficient, Reynolds number, and terminal veloc
     193           0 :     do m = 1, nbin
     194             : 
     195             :        ! Initialize accuracy and counter
     196           0 :        eps_crr = eps_max + 1.0_r8  ![frc] Current relative accuracy
     197           0 :        itr_idx = 0                 ![idx] Counting index
     198             : 
     199             :        ! Initial guess for vlc_grv is exact for Re < 0.1
     200           0 :        vlc_grv(m) = vlc_stk(m)     ![m s-1]
     201           0 :        eps_loop: do while(eps_crr > eps_max)
     202             : 
     203             :           ! Save terminal velocity for convergence test
     204           0 :           vlc_grv_old = vlc_grv(m) ![m s-1]
     205           0 :           ryn_nbr_grv(m) = vlc_grv(m) * dmt_vwr(m) / vsc_knm_atm !SeP97 p.460
     206             : 
     207             :           ! Update drag coefficient based on new Reynolds number
     208           0 :           if (ryn_nbr_grv(m) < 0.1_r8) then
     209           0 :              cff_drg_grv(m) = 24.0_r8 / ryn_nbr_grv(m) !Stokes' law Sep97 p.463 (8.32)
     210           0 :           else if (ryn_nbr_grv(m) < 2.0_r8) then
     211             :              cff_drg_grv(m) = (24.0_r8/ryn_nbr_grv(m)) *    &
     212             :                   (1.0_r8 + 3.0_r8*ryn_nbr_grv(m)/16.0_r8 + &
     213             :                   9.0_r8*ryn_nbr_grv(m)*ryn_nbr_grv(m)*     &
     214           0 :                   log(2.0_r8*ryn_nbr_grv(m))/160.0_r8)        !Sep97 p.463 (8.32)
     215           0 :           else if (ryn_nbr_grv(m) < 500.0_r8) then
     216             :              cff_drg_grv(m) = (24.0_r8/ryn_nbr_grv(m)) * &
     217           0 :                   (1.0_r8 + 0.15_r8*ryn_nbr_grv(m)**0.687_r8) !Sep97 p.463 (8.32)
     218           0 :           else if (ryn_nbr_grv(m) < 2.0e5_r8) then
     219           0 :              cff_drg_grv(m) = 0.44_r8                         !Sep97 p.463 (8.32)
     220             :           else
     221           0 :              write(iulog,'(a,es9.2)') "ryn_nbr_grv(m) = ",ryn_nbr_grv(m)
     222           0 :              call endrun ('Dustini error: Reynolds number too large in stk_crc_get()')
     223             :           endif
     224             : 
     225             :           ! Update terminal velocity based on new Reynolds number and drag coeff
     226             :           ! [m s-1] Terminal veloc SeP97 p.467 (8.44)
     227             :           vlc_grv(m) = sqrt(4.0_r8 * gravit * dmt_vwr(m) * slp_crc(m) * dust_density / &
     228           0 :                (3.0_r8*cff_drg_grv(m)*dns_mdp))   
     229           0 :           eps_crr = abs((vlc_grv(m)-vlc_grv_old)/vlc_grv(m)) !Relative convergence
     230           0 :           if (itr_idx == 12) then
     231             :              ! Numerical pingpong may occur when Re = 0.1, 2.0, or 500.0
     232             :              ! due to discontinuities in derivative of drag coefficient
     233           0 :              vlc_grv(m) = 0.5_r8 * (vlc_grv(m)+vlc_grv_old)  ! [m s-1]
     234             :           endif
     235           0 :           if (itr_idx > 20) then
     236           0 :              write(iulog,*) 'Dustini error: Terminal velocity not converging ',&
     237           0 :                   ' in stk_crc_get(), breaking loop...'
     238             :              ! to next iteration
     239           0 :              exit eps_loop
     240             :           endif
     241           0 :           itr_idx = itr_idx + 1
     242             :        end do eps_loop  !end while
     243             :     end do   !end loop over size
     244             : 
     245             :     ! Compute factors to convert Stokes' settling velocities to
     246             :     ! actual settling velocities
     247           0 :     do m = 1, nbin
     248           0 :        stk_crc(m) = vlc_grv(m) / vlc_stk(m)
     249             :     end do
     250             :     
     251           0 :   end subroutine dust_set_params
     252             : 
     253             : end module dust_common