LCOV - code coverage report
Current view: | top level - chemistry/aerosol - dust_common.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 65 | 77 | 84.4 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 1 | 100.0 %

              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        1536 :   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        3072 :     real(r8) :: dmt_min(nbin)      ![m] Size grid minimum
          61        3072 :     real(r8) :: dmt_max(nbin)      ![m] Size grid maximum
          62        3072 :     real(r8) :: dmt_ctr(nbin)      ![m] Diameter at bin center
          63        3072 :     real(r8) :: dmt_dlt(nbin)      ![m] Width of size bin
          64        3072 :     real(r8) :: slp_crc(nbin)      ![frc] Slip correction factor
          65        3072 :     real(r8) :: vlm_rsl(nbin)      ![m3 m-3] Volume concentration resolved
          66        3072 :     real(r8) :: vlc_stk(nbin)      ![m s-1] Stokes settling velocity
          67        3072 :     real(r8) :: vlc_grv(nbin)      ![m s-1] Settling velocity
          68        3072 :     real(r8) :: ryn_nbr_grv(nbin)  ![frc] Reynolds number at terminal velocity
          69        3072 :     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        1536 :     stk_crc(:) = nan
          92        1536 :     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        6144 :     do n = 1, nbin
         109        4608 :        dmt_min(n) = dmt_grd(n)                       ![m] Max diameter in bin
         110        4608 :        dmt_max(n) = dmt_grd(n+1)                     ![m] Min diameter in bin
         111        4608 :        dmt_ctr(n) = 0.5_r8 * (dmt_min(n)+dmt_max(n)) ![m] Diameter at bin ctr
         112        6144 :        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        6144 :     do n = 1, nbin
         130        4608 :        series_ratio = (dmt_max(n)/dmt_min(n))**(1.0_r8/sz_nbr)
         131        4608 :        sz_min(1) = dmt_min(n)
         132      921600 :        do m = 2, sz_nbr                            ! Loop starts at 2
         133      921600 :           sz_min(m) = sz_min(m-1) * series_ratio
         134             :        end do
         135             :
         136             :        ! Derived grid values
         137      921600 :        do m = 1, sz_nbr-1                          ! Loop ends at sz_nbr-1
         138      921600 :           sz_max(m) = sz_min(m+1)                  ! [m]
         139             :        end do
         140        4608 :        sz_max(sz_nbr) = dmt_max(n)                 ! [m]
         141             :
         142             :        ! Final derived grid values
         143      926208 :        do m = 1, sz_nbr
         144      921600 :           sz_ctr(m) = 0.5_r8 * (sz_min(m)+sz_max(m))
         145      926208 :           sz_dlt(m) = sz_max(m)-sz_min(m)
         146             :        end do
         147        4608 :        lngsdsqrttwopi_rcp = 1.0_r8 / (ln_gsd*sqrt(2.0_r8*pi))
         148        4608 :        dmt_vwr(n) = 0.0_r8 ! [m] Mass wgted diameter resolved
         149        4608 :        vlm_rsl(n) = 0.0_r8 ! [m3 m-3] Volume concentration resolved
         150      926208 :        do m = 1, sz_nbr
         151             :           ! Evaluate lognormal distribution for these sizes (call lgn_evl)
         152      921600 :           tmp = log(sz_ctr(m)/dmt_nma) / ln_gsd
         153      921600 :           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      921600 :                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      926208 :                lgn_dst * sz_dlt(m)                ![# m-3] Number concentrn
         161             :        end do
         162        6144 :        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        1536 :     eps_max = 1.0e-4_r8
         169        1536 :     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        1536 :          (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        1536 :          (100000._r8*sqrt(8.0_r8/(pi*rair*295.0_r8)))
         175        1536 :     vsc_knm_atm = vsc_dyn_atm / dns_mdp ![m2 s-1] Kinematic viscosity of air
         176             :
         177        6144 :     do m = 1, nbin
         178        4608 :        slp_crc(m) = 1.0_r8 + 2.0_r8 * mfp_atm *                      &
         179        4608 :             (1.257_r8+0.4_r8*exp(-1.1_r8*dmt_vwr(m)/(2.0_r8*mfp_atm))) / &
         180        4608 :             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        6144 :             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        6144 :     do m = 1, nbin
         194             :
         195             :        ! Initialize accuracy and counter
         196        4608 :        eps_crr = eps_max + 1.0_r8  ![frc] Current relative accuracy
         197        4608 :        itr_idx = 0                 ![idx] Counting index
         198             :
         199             :        ! Initial guess for vlc_grv is exact for Re < 0.1
         200        4608 :        vlc_grv(m) = vlc_stk(m)     ![m s-1]
         201       10752 :        eps_loop: do while(eps_crr > eps_max)
         202             :
         203             :           ! Save terminal velocity for convergence test
         204        4608 :           vlc_grv_old = vlc_grv(m) ![m s-1]
         205        4608 :           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        4608 :           if (ryn_nbr_grv(m) < 0.1_r8) then
         209        4608 :              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        4608 :                (3.0_r8*cff_drg_grv(m)*dns_mdp))
         229        4608 :           eps_crr = abs((vlc_grv(m)-vlc_grv_old)/vlc_grv(m)) !Relative convergence
         230        4608 :           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        4608 :           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        4608 :           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        6144 :     do m = 1, nbin
         248        6144 :        stk_crc(m) = vlc_grv(m) / vlc_stk(m)
         249             :     end do
         250             :
         251        1536 :   end subroutine dust_set_params
         252             :
         253             : end module dust_common

Generated by: LCOV version 1.14