LCOV - code coverage report
Current view: | top level - chemistry/utils - modal_aero_wateruptake.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 442 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 9 | 0.0 %

              Line data    Source code
           1             : module modal_aero_wateruptake
           2             :
           3             : !   RCE 07.04.13:  Adapted from MIRAGE2 code
           4             :
           5             : use shr_kind_mod,     only: r8 => shr_kind_r8
           6             : use physconst,        only: pi, rhoh2o, rair
           7             : use ppgrid,           only: pcols, pver
           8             : use physics_types,    only: physics_state
           9             : use physics_buffer,   only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field
          10             :
          11             : use wv_saturation,    only: qsat_water
          12             : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props, &
          13             :                             rad_cnst_get_mode_props
          14             : use cam_history,      only: addfld, add_default, outfld, horiz_only
          15             : use cam_logfile,      only: iulog
          16             : use ref_pres,         only: top_lev => clim_modal_aero_top_lev
          17             : use phys_control,     only: phys_getopts
          18             : use cam_abortutils,   only: endrun
          19             :
          20             : implicit none
          21             : private
          22             : save
          23             :
          24             : public :: &
          25             :    modal_aero_wateruptake_init, &
          26             :    modal_aero_wateruptake_dr,   &
          27             :    modal_aero_wateruptake_sub,  &
          28             :    modal_aero_kohler
          29             :
          30             : public :: modal_aero_wateruptake_reg
          31             :
          32             : real(r8), parameter :: third = 1._r8/3._r8
          33             : real(r8), parameter :: pi43  = pi*4.0_r8/3.0_r8
          34             :
          35             :
          36             : ! Physics buffer indices
          37             : integer :: cld_idx        = 0
          38             : integer :: dgnum_idx      = 0
          39             : integer :: dgnumwet_idx   = 0
          40             : integer :: sulfeq_idx     = 0
          41             : integer :: wetdens_ap_idx = 0
          42             : integer :: qaerwat_idx    = 0
          43             : integer :: hygro_idx      = 0
          44             : integer :: dryvol_idx     = 0
          45             : integer :: dryrad_idx     = 0
          46             : integer :: drymass_idx    = 0
          47             : integer :: so4dryvol_idx  = 0
          48             : integer :: naer_idx       = 0
          49             :
          50             :
          51             : logical, public :: modal_strat_sulfate = .false.   ! If .true. then MAM sulfate surface area density used in stratospheric heterogeneous chemistry
          52             :
          53             : !===============================================================================
          54             : contains
          55             : !===============================================================================
          56             :
          57           0 : subroutine modal_aero_wateruptake_reg()
          58             :
          59             :   use physics_buffer,   only: pbuf_add_field, dtype_r8
          60             :   use rad_constituents, only: rad_cnst_get_info
          61             :
          62             :    integer :: nmodes
          63             :
          64           0 :    call rad_cnst_get_info(0, nmodes=nmodes)
          65           0 :    call pbuf_add_field('DGNUMWET',   'global',  dtype_r8, (/pcols, pver, nmodes/), dgnumwet_idx)
          66           0 :    call pbuf_add_field('WETDENS_AP', 'physpkg', dtype_r8, (/pcols, pver, nmodes/), wetdens_ap_idx)
          67             :
          68             :    ! 1st order rate for direct conversion of strat. cloud water to precip (1/s)
          69           0 :    call pbuf_add_field('QAERWAT',    'physpkg', dtype_r8, (/pcols, pver, nmodes/), qaerwat_idx)
          70             :
          71           0 :    if (modal_strat_sulfate) then
          72           0 :       call pbuf_add_field('MAMH2SO4EQ', 'global',  dtype_r8, (/pcols, pver, nmodes/), sulfeq_idx)
          73             :    end if
          74             :
          75             :
          76           0 : end subroutine modal_aero_wateruptake_reg
          77             :
          78             : !===============================================================================
          79             : !===============================================================================
          80             :
          81           0 : subroutine modal_aero_wateruptake_init(pbuf2d)
          82           0 :    use time_manager,  only: is_first_step
          83             :    use physics_buffer,only: pbuf_set_field
          84             :    use infnan,       only : nan, assignment(=)
          85             :
          86             :    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
          87             :    real(r8) :: real_nan
          88             :
          89             :    integer :: m, nmodes
          90             :    logical :: history_aerosol      ! Output the MAM aerosol variables and tendencies
          91             :
          92             :    character(len=3) :: trnum       ! used to hold mode number (as characters)
          93             :    !----------------------------------------------------------------------------
          94             :
          95           0 :    real_nan = nan
          96             :
          97           0 :    cld_idx        = pbuf_get_index('CLD')
          98           0 :    dgnum_idx      = pbuf_get_index('DGNUM')
          99             :
         100           0 :    hygro_idx      = pbuf_get_index('HYGRO')
         101           0 :    dryvol_idx     = pbuf_get_index('DRYVOL')
         102           0 :    dryrad_idx     = pbuf_get_index('DRYRAD')
         103           0 :    drymass_idx    = pbuf_get_index('DRYMASS')
         104           0 :    so4dryvol_idx  = pbuf_get_index('SO4DRYVOL')
         105           0 :    naer_idx       = pbuf_get_index('NAER')
         106             :
         107             :    ! assume for now that will compute wateruptake for climate list modes only
         108             :
         109           0 :    call rad_cnst_get_info(0, nmodes=nmodes)
         110             :
         111           0 :    do m = 1, nmodes
         112           0 :       write(trnum, '(i3.3)') m
         113             :       call addfld('dgnd_a'//trnum(2:3), (/ 'lev' /), 'A', 'm', &
         114           0 :          'dry dgnum, interstitial, mode '//trnum(2:3))
         115             :       call addfld('dgnw_a'//trnum(2:3), (/ 'lev' /), 'A', 'm', &
         116           0 :          'wet dgnum, interstitial, mode '//trnum(2:3))
         117             :       call addfld('wat_a'//trnum(3:3),  (/ 'lev' /), 'A', 'm', &
         118           0 :          'aerosol water, interstitial, mode '//trnum(2:3))
         119             :
         120             :       ! determine default variables
         121           0 :       call phys_getopts(history_aerosol_out = history_aerosol)
         122             :
         123           0 :       if (history_aerosol) then
         124           0 :          call add_default('dgnd_a'//trnum(2:3), 1, ' ')
         125           0 :          call add_default('dgnw_a'//trnum(2:3), 1, ' ')
         126           0 :          call add_default('wat_a'//trnum(3:3),  1, ' ')
         127             :       endif
         128             :
         129             :    end do
         130             :
         131           0 :    call addfld('PM25',     (/ 'lev' /), 'A', 'kg/m3', 'PM2.5 concentration')
         132           0 :    call addfld('PM25_SRF', horiz_only,  'A', 'kg/m3', 'surface PM2.5 concentration')
         133             :
         134           0 :    if (is_first_step()) then
         135             :       ! initialize fields in physics buffer
         136           0 :       call pbuf_set_field(pbuf2d, dgnumwet_idx, 0.0_r8)
         137           0 :       if (modal_strat_sulfate) then
         138             :       ! initialize fields in physics buffer to NaN (not a number)
         139             :       ! so model will crash if used before initialization
         140           0 :          call pbuf_set_field(pbuf2d, sulfeq_idx, real_nan)
         141             :       endif
         142             :    endif
         143             :
         144           0 : end subroutine modal_aero_wateruptake_init
         145             :
         146             : !===============================================================================
         147             :
         148             :
         149           0 : subroutine modal_aero_wateruptake_dr(state, pbuf, list_idx_in, dgnumdry_m, dgnumwet_m, &
         150             :                                      qaerwat_m, wetdens_m, hygro_m, dryvol_m, dryrad_m, drymass_m,&
         151             :                                      so4dryvol_m, naer_m)
         152             : !-----------------------------------------------------------------------
         153             : !
         154             : ! CAM specific driver for modal aerosol water uptake code.
         155             : !
         156             : ! *** N.B. *** The calculation has been enabled for diagnostic mode lists
         157             : !              via optional arguments.  If the list_idx arg is present then
         158             : !              all the optional args must be present.
         159             : !
         160             : !-----------------------------------------------------------------------
         161             :
         162           0 :    use time_manager,  only: is_first_step
         163             :    use cam_history,   only: outfld, fieldname_len
         164             :    use tropopause,    only: tropopause_find_cam, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
         165             :    ! Arguments
         166             :    type(physics_state), target, intent(in)    :: state          ! Physics state variables
         167             :    type(physics_buffer_desc),   pointer       :: pbuf(:)        ! physics buffer
         168             :
         169             :    integer,  optional,          intent(in)    :: list_idx_in
         170             :    real(r8), optional,          pointer       :: dgnumdry_m(:,:,:)
         171             :    real(r8), optional,          pointer       :: dgnumwet_m(:,:,:)
         172             :    real(r8), optional,          pointer       :: qaerwat_m(:,:,:)
         173             :    real(r8), optional,          pointer       :: wetdens_m(:,:,:)
         174             :    real(r8), optional,          pointer       :: hygro_m(:,:,:)
         175             :    real(r8), optional,          pointer       :: dryvol_m(:,:,:)
         176             :    real(r8), optional,          pointer       :: dryrad_m(:,:,:)
         177             :    real(r8), optional,          pointer       :: drymass_m(:,:,:)
         178             :    real(r8), optional,          pointer       :: so4dryvol_m(:,:,:)
         179             :    real(r8), optional,          pointer       :: naer_m(:,:,:)
         180             :
         181             :    ! local variables
         182             :
         183             :    integer  :: lchnk              ! chunk index
         184             :    integer  :: ncol               ! number of columns
         185             :    integer  :: list_idx           ! radiative constituents list index
         186             :    integer  :: stat
         187             :
         188             :    integer :: i, k, l, m
         189             :    integer :: itim_old
         190             :    integer :: nmodes
         191             :    integer :: nspec
         192             :    integer :: tropLev(pcols)
         193             :
         194             :    character(len=fieldname_len+3) :: fieldname
         195             :
         196           0 :    real(r8), pointer :: h2ommr(:,:) ! specific humidity
         197           0 :    real(r8), pointer :: t(:,:)      ! temperatures (K)
         198           0 :    real(r8), pointer :: pmid(:,:)   ! layer pressure (Pa)
         199             :
         200           0 :    real(r8), pointer :: cldn(:,:)      ! layer cloud fraction (0-1)
         201           0 :    real(r8), pointer :: dgncur_a(:,:,:)
         202           0 :    real(r8), pointer :: dgncur_awet(:,:,:)
         203           0 :    real(r8), pointer :: wetdens(:,:,:)
         204           0 :    real(r8), pointer :: qaerwat(:,:,:)
         205             :
         206           0 :    real(r8), pointer :: raer(:,:)        ! aerosol species MRs (kg/kg)
         207           0 :    real(r8), pointer :: maer(:,:,:)      ! accumulated aerosol mode MRs
         208           0 :    real(r8), pointer :: hygro(:,:,:)     ! volume-weighted mean hygroscopicity (--)
         209           0 :    real(r8), pointer :: naer(:,:,:)      ! aerosol number MR (bounded!) (#/kg-air)
         210           0 :    real(r8), pointer :: dryvol(:,:,:)    ! single-particle-mean dry volume (m3)
         211           0 :    real(r8), pointer :: so4dryvol(:,:,:) ! single-particle-mean so4 dry volume (m3)
         212           0 :    real(r8), pointer :: drymass(:,:,:)   ! single-particle-mean dry mass  (kg)
         213           0 :    real(r8), pointer :: dryrad(:,:,:)    ! dry volume mean radius of aerosol (m)
         214             :
         215           0 :    real(r8), allocatable :: wetrad(:,:,:)    ! wet radius of aerosol (m)
         216           0 :    real(r8), allocatable :: wetvol(:,:,:)    ! single-particle-mean wet volume (m3)
         217           0 :    real(r8), allocatable :: wtrvol(:,:,:)    ! single-particle-mean water volume in wet aerosol (m3)
         218             :
         219           0 :    real(r8), allocatable :: rhcrystal(:)
         220           0 :    real(r8), allocatable :: rhdeliques(:)
         221           0 :    real(r8), allocatable :: specdens_1(:)
         222             :
         223           0 :    real(r8), pointer :: sulfeq(:,:,:) ! H2SO4 equilibrium mixing ratios over particles (mol/mol)
         224           0 :    real(r8), allocatable :: wtpct(:,:,:)  ! sulfate aerosol composition, weight % H2SO4
         225           0 :    real(r8), allocatable :: sulden(:,:,:) ! sulfate aerosol mass density (g/cm3)
         226             :
         227             :    real(r8) :: specdens, so4specdens
         228             :    real(r8) :: sigmag
         229           0 :    real(r8), allocatable :: alnsg(:)
         230             :    real(r8) :: rh(pcols,pver)        ! relative humidity (0-1)
         231             :    real(r8) :: dmean, qh2so4_equilib, wtpct_mode, sulden_mode
         232             :
         233             :    real(r8) :: es(pcols)             ! saturation vapor pressure
         234             :    real(r8) :: qs(pcols)             ! saturation specific humidity
         235             :
         236             :    real(r8) :: pm25(pcols,pver)      ! PM2.5 diagnostics
         237             :    real(r8) :: rhoair(pcols,pver)
         238             :
         239             :    character(len=3) :: trnum       ! used to hold mode number (as characters)
         240             :    character(len=32) :: spectype
         241             :    !-----------------------------------------------------------------------
         242             :
         243           0 :    lchnk = state%lchnk
         244           0 :    ncol = state%ncol
         245             :
         246           0 :    list_idx = 0
         247           0 :    if (present(list_idx_in)) then
         248           0 :       list_idx = list_idx_in
         249             :
         250             :       ! check that all optional args are present
         251             :       if (.not. present(dgnumdry_m) .or. .not. present(dgnumwet_m) .or. &
         252           0 :           .not. present(qaerwat_m)  .or. .not. present(wetdens_m)) then
         253             :          call endrun('modal_aero_wateruptake_dr called for'// &
         254           0 :                      'diagnostic list but required args not present')
         255             :       end if
         256             :
         257             :       ! arrays for diagnostic calculations must be associated
         258             :       if (.not. associated(dgnumdry_m) .or. .not. associated(dgnumwet_m) .or. &
         259           0 :           .not. associated(qaerwat_m)  .or. .not. associated(wetdens_m)) then
         260             :          call endrun('modal_aero_wateruptake_dr called for'// &
         261           0 :                      'diagnostic list but required args not associated')
         262             :       end if
         263             :
         264           0 :       if (modal_strat_sulfate) then
         265             :          call endrun('modal_aero_wateruptake_dr cannot be called with optional arguments and'// &
         266           0 :                      ' having modal_strat_sulfate set to true')
         267             :       end if
         268             :    end if
         269             :
         270             :    ! loop over all aerosol modes
         271           0 :    call rad_cnst_get_info(list_idx, nmodes=nmodes)
         272             :
         273           0 :    if (modal_strat_sulfate) then
         274           0 :      call pbuf_get_field(pbuf,  sulfeq_idx, sulfeq )
         275             :    endif
         276             :
         277             :    allocate( &
         278           0 :       wetrad(pcols,pver,nmodes),   &
         279           0 :       wetvol(pcols,pver,nmodes),   &
         280           0 :       wtrvol(pcols,pver,nmodes),   &
         281           0 :       wtpct(pcols,pver,nmodes),    &
         282           0 :       sulden(pcols,pver,nmodes),   &
         283           0 :       rhcrystal(nmodes),           &
         284           0 :       rhdeliques(nmodes),          &
         285             :       specdens_1(nmodes),          &
         286           0 :       alnsg(nmodes)           )
         287             :
         288           0 :    wtpct(:,:,:)     = 75._r8
         289           0 :    sulden(:,:,:)    = 1.923_r8
         290             :
         291           0 :    if (list_idx == 0) then
         292           0 :       call pbuf_get_field(pbuf, dgnum_idx,      dgncur_a )
         293           0 :       call pbuf_get_field(pbuf, dgnumwet_idx,   dgncur_awet )
         294           0 :       call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens)
         295           0 :       call pbuf_get_field(pbuf, qaerwat_idx,    qaerwat)
         296           0 :       call pbuf_get_field(pbuf, hygro_idx,       hygro)
         297           0 :       call pbuf_get_field(pbuf, dryvol_idx,      dryvol)
         298           0 :       call pbuf_get_field(pbuf, dryrad_idx,      dryrad)
         299           0 :       call pbuf_get_field(pbuf, drymass_idx,     drymass)
         300           0 :       call pbuf_get_field(pbuf, so4dryvol_idx,   so4dryvol)
         301           0 :       call pbuf_get_field(pbuf, naer_idx,        naer)
         302             :
         303           0 :       if (is_first_step()) then
         304           0 :          dgncur_awet(:,:,:) = dgncur_a(:,:,:)
         305             :       end if
         306             :    else
         307           0 :       dgncur_a    => dgnumdry_m
         308           0 :       dgncur_awet => dgnumwet_m
         309           0 :       qaerwat     => qaerwat_m
         310           0 :       wetdens     => wetdens_m
         311           0 :       hygro       => hygro_m
         312           0 :       dryvol      => dryvol_m
         313           0 :       dryrad      => dryrad_m
         314           0 :       drymass     => drymass_m
         315           0 :       so4dryvol   => so4dryvol_m
         316           0 :       naer        => naer_m
         317             :    end if
         318             :
         319           0 :    if (modal_strat_sulfate) then
         320             :       ! get tropopause level
         321             :       !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
         322           0 :       tropLev(:) = 0
         323             :       !REMOVECAM_END
         324           0 :       call tropopause_find_cam(state, tropLev, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE)
         325             :    endif
         326             :
         327           0 :    h2ommr => state%q(:,:,1)
         328           0 :    t      => state%t
         329           0 :    pmid   => state%pmid
         330             :
         331           0 :    allocate( maer(pcols,pver,nmodes))
         332           0 :    maer(:,:,:) = 0._r8
         333             :
         334           0 :    do m = 1, nmodes
         335             :
         336             :       call rad_cnst_get_mode_props(list_idx, m, sigmag=sigmag, &
         337           0 :          rhcrystal=rhcrystal(m), rhdeliques=rhdeliques(m))
         338             :
         339             :       ! get mode info
         340           0 :       call rad_cnst_get_info(list_idx, m, nspec=nspec)
         341             :
         342           0 :       do l = 1, nspec
         343             :
         344             :          ! accumulate the aerosol masses of each mode
         345           0 :          call rad_cnst_get_aer_mmr(0, m, l, 'a', state, pbuf, raer)
         346           0 :          maer(:ncol,:,m)= maer(:ncol,:,m) + raer(:ncol,:)
         347             :
         348             :          ! get species interstitial mixing ratio ('a')
         349             :          call rad_cnst_get_aer_props(list_idx, m, l, density_aer=specdens, &
         350           0 :                                      spectype=spectype)
         351             :
         352           0 :          if (modal_strat_sulfate .and. (trim(spectype).eq.'sulfate')) then
         353           0 :             so4specdens=specdens
         354             :          end if
         355             :
         356           0 :          if (l == 1) then
         357             :             ! save off these values to be used as defaults
         358           0 :             specdens_1(m) = specdens
         359             :          end if
         360             :
         361             :       end do
         362             :
         363           0 :       alnsg(m) = log(sigmag)
         364             :
         365           0 :       if (modal_strat_sulfate) then
         366           0 :          do k = top_lev, pver
         367           0 :             do i = 1, ncol
         368           0 :                dmean = dgncur_awet(i,k,m)*exp(1.5_r8*alnsg(m)**2)
         369           0 :                call calc_h2so4_equilib_mixrat( t(i,k), pmid(i,k), h2ommr(i,k), dmean, &
         370           0 :                                                qh2so4_equilib, wtpct_mode, sulden_mode )
         371           0 :                sulfeq(i,k,m)  = qh2so4_equilib
         372           0 :                wtpct(i,k,m)   = wtpct_mode
         373           0 :                sulden(i,k,m)  = sulden_mode
         374             :             end do    ! i = 1, ncol
         375             :          end do    ! k = top_lev, pver
         376             :
         377           0 :           fieldname = ' '
         378           0 :           write(fieldname,fmt='(a,i1)') 'wtpct_a',m
         379           0 :           call outfld(fieldname,wtpct(1:ncol,1:pver,m), ncol, lchnk )
         380             :
         381           0 :           fieldname = ' '
         382           0 :           write(fieldname,fmt='(a,i1)') 'sulfeq_a',m
         383           0 :           call outfld(fieldname,sulfeq(1:ncol,1:pver,m), ncol, lchnk )
         384             :
         385           0 :           fieldname = ' '
         386           0 :           write(fieldname,fmt='(a,i1)') 'sulden_a',m
         387           0 :           call outfld(fieldname,sulden(1:ncol,1:pver,m), ncol, lchnk )
         388             :
         389             :        end if
         390             :
         391             :     end do    ! m = 1, nmodes
         392             :
         393             :    ! relative humidity calc
         394             :
         395           0 :    itim_old    =  pbuf_old_tim_idx()
         396           0 :    call pbuf_get_field(pbuf, cld_idx, cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         397             :
         398           0 :    do k = top_lev, pver
         399           0 :       call qsat_water(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol), qs(1:ncol), ncol)
         400           0 :       do i = 1, ncol
         401           0 :          if (qs(i) > h2ommr(i,k)) then
         402           0 :             rh(i,k) = h2ommr(i,k)/qs(i)
         403             :          else
         404           0 :             rh(i,k) = 0.98_r8
         405             :          endif
         406           0 :          rh(i,k) = max(rh(i,k), 0.0_r8)
         407           0 :          rh(i,k) = min(rh(i,k), 0.98_r8)
         408           0 :          if (cldn(i,k) .lt. 1.0_r8) then
         409           0 :             rh(i,k) = (rh(i,k) - cldn(i,k)) / (1.0_r8 - cldn(i,k))  ! clear portion
         410             :          end if
         411           0 :          rh(i,k) = max(rh(i,k), 0.0_r8)
         412           0 :          rhoair(i,k) = pmid(i,k)/(rair*t(i,k))
         413             :       end do
         414             :    end do
         415             :
         416             :    call modal_aero_wateruptake_sub( &
         417             :       ncol, nmodes, rhcrystal, rhdeliques, dryrad, &
         418             :       hygro, rh, dryvol, so4dryvol, so4specdens, tropLev, &
         419           0 :       wetrad, wetvol, wtrvol, sulden, wtpct)
         420             :
         421           0 :    qaerwat = 0.0_r8
         422             :
         423           0 :    do m = 1, nmodes
         424             :
         425           0 :       do k = top_lev, pver
         426           0 :          do i = 1, ncol
         427             :
         428           0 :             dgncur_awet(i,k,m) = dgncur_a(i,k,m) * (wetrad(i,k,m)/dryrad(i,k,m))
         429           0 :             qaerwat(i,k,m)     = rhoh2o*naer(i,k,m)*wtrvol(i,k,m)
         430             :
         431             :             ! compute aerosol wet density (kg/m3)
         432           0 :             if (wetvol(i,k,m) > 1.0e-30_r8) then
         433           0 :                wetdens(i,k,m) = (drymass(i,k,m) + rhoh2o*wtrvol(i,k,m))/wetvol(i,k,m)
         434             :             else
         435           0 :                wetdens(i,k,m) = specdens_1(m)
         436             :             end if
         437             :          end do
         438             :       end do
         439             :
         440             :    end do    ! modes
         441             :
         442           0 :    if (list_idx == 0) then
         443             :
         444           0 :       pm25(:,:)=0._r8
         445             :
         446           0 :       do m = 1, nmodes
         447             :          ! output to history
         448           0 :          write( trnum, '(i3.3)' ) m
         449           0 :          call outfld( 'wat_a'//trnum(3:3),  qaerwat(:,:,m),     pcols, lchnk)
         450           0 :          call outfld( 'dgnd_a'//trnum(2:3), dgncur_a(:,:,m),    pcols, lchnk)
         451           0 :          call outfld( 'dgnw_a'//trnum(2:3), dgncur_awet(:,:,m), pcols, lchnk)
         452             :
         453             :          ! calculate PM2.5 diagnostics -- dgncur_a is zero above top_lev
         454           0 :          do k = top_lev, pver
         455           0 :             do i=1,ncol
         456           0 :                pm25(i,k) = pm25(i,k)+maer(i,k,m)*(1._r8-(0.5_r8 - 0.5_r8*erf(log(2.5e-6_r8/dgncur_a(i,k,m))/ &
         457           0 :                                                  (2._r8**0.5_r8*alnsg(m)))))*rhoair(i,k)
         458             :             end do
         459             :          end do
         460             :       end do
         461             :
         462           0 :       call outfld('PM25',     pm25(:,:),    pcols, lchnk)
         463           0 :       call outfld('PM25_SRF', pm25(:,pver), pcols, lchnk)
         464             :
         465             :    end if
         466             :
         467           0 :    deallocate(maer, alnsg)
         468           0 :    deallocate( &
         469           0 :       wetrad, wetvol, wtrvol, wtpct, sulden, rhcrystal, rhdeliques, specdens_1 )
         470           0 : end subroutine modal_aero_wateruptake_dr
         471             :
         472             : !===============================================================================
         473             :
         474           0 : subroutine modal_aero_wateruptake_sub( &
         475           0 :    ncol, nmodes, rhcrystal, rhdeliques, dryrad, &
         476           0 :    hygro, rh, dryvol, so4dryvol, so4specdens, troplev, &
         477           0 :    wetrad, wetvol, wtrvol, sulden, wtpct)
         478             :
         479             : !-----------------------------------------------------------------------
         480             : !
         481             : ! Purpose: Compute aerosol wet radius
         482             : !
         483             : ! Method:  Kohler theory
         484             : !
         485             : ! Author:  S. Ghan
         486             : !
         487             : !-----------------------------------------------------------------------
         488             :
         489             :
         490             :    ! Arguments
         491             :    integer, intent(in)  :: ncol                    ! number of columns
         492             :    integer, intent(in)  :: nmodes
         493             :    integer, intent(in)  :: troplev(:)
         494             :
         495             :    real(r8), intent(in) :: rhcrystal(:)
         496             :    real(r8), intent(in) :: rhdeliques(:)
         497             :    real(r8), intent(in) :: dryrad(:,:,:)         ! dry volume mean radius of aerosol (m)
         498             :    real(r8), intent(in) :: hygro(:,:,:)          ! volume-weighted mean hygroscopicity (--)
         499             :    real(r8), intent(in) :: rh(:,:)               ! relative humidity (0-1)
         500             :    real(r8), intent(in) :: dryvol(:,:,:)         ! dry volume of single aerosol (m3)
         501             :    real(r8), intent(in) :: so4dryvol(:,:,:)      ! dry volume of sulfate in single aerosol (m3)
         502             :    real(r8), intent(in) :: so4specdens           ! mass density sulfate in single aerosol (kg/m3)
         503             :    real(r8), intent(in) :: wtpct(:,:,:)          ! sulfate aerosol composition, weight % H2SO4
         504             :    real(r8), intent(in) :: sulden(:,:,:)         ! sulfate aerosol mass density (g/cm3)
         505             :
         506             :    real(r8), intent(out) :: wetrad(:,:,:)        ! wet radius of aerosol (m)
         507             :    real(r8), intent(out) :: wetvol(:,:,:)        ! single-particle-mean wet volume (m3)
         508             :    real(r8), intent(out) :: wtrvol(:,:,:)        ! single-particle-mean water volume in wet aerosol (m3)
         509             :
         510             :    ! local variables
         511             :
         512             :    integer :: i, k, m
         513             :
         514             :    real(r8) :: hystfac                ! working variable for hysteresis
         515             :    !-----------------------------------------------------------------------
         516             :
         517             :
         518             :    ! loop over all aerosol modes
         519           0 :    do m = 1, nmodes
         520             :
         521           0 :       hystfac = 1.0_r8 / max(1.0e-5_r8, (rhdeliques(m) - rhcrystal(m)))
         522             :
         523           0 :       do k = top_lev, pver
         524           0 :          do i = 1, ncol
         525             :
         526           0 :             if ( modal_strat_sulfate .and. (k<troplev(i))) then
         527           0 :                wetvol(i,k,m) = dryvol(i,k,m)-so4dryvol(i,k,m)
         528           0 :                wetvol(i,k,m) = wetvol(i,k,m)+so4dryvol(i,k,m)*so4specdens/sulden(i,k,m)/wtpct(i,k,m)/10._r8
         529           0 :                wetvol(i,k,m) = max(wetvol(i,k,m), dryvol(i,k,m))
         530           0 :                wetrad(i,k,m) = (wetvol(i,k,m)/pi43)**third
         531           0 :                wetrad(i,k,m) = max(wetrad(i,k,m), dryrad(i,k,m))
         532           0 :                wtrvol(i,k,m) = wetvol(i,k,m) - dryvol(i,k,m)
         533           0 :                wtrvol(i,k,m) = max(wtrvol(i,k,m), 0.0_r8)
         534             :             else
         535             :               ! compute wet radius for each mode
         536           0 :               call modal_aero_kohler(dryrad(i:i,k,m), hygro(i:i,k,m), rh(i:i,k), wetrad(i:i,k,m), 1)
         537             :
         538           0 :               wetrad(i,k,m) = max(wetrad(i,k,m), dryrad(i,k,m))
         539           0 :               wetvol(i,k,m) = pi43*wetrad(i,k,m)**3
         540           0 :               wetvol(i,k,m) = max(wetvol(i,k,m), dryvol(i,k,m))
         541           0 :               wtrvol(i,k,m) = wetvol(i,k,m) - dryvol(i,k,m)
         542           0 :               wtrvol(i,k,m) = max(wtrvol(i,k,m), 0.0_r8)
         543             :
         544             :               ! apply simple treatment of deliquesence/crystallization hysteresis
         545             :               ! for rhcrystal < rh < rhdeliques, aerosol water is a fraction of
         546             :               ! the "upper curve" value, and the fraction is a linear function of rh
         547           0 :               if (rh(i,k) < rhcrystal(m)) then
         548           0 :                  wetrad(i,k,m) = dryrad(i,k,m)
         549           0 :                  wetvol(i,k,m) = dryvol(i,k,m)
         550           0 :                  wtrvol(i,k,m) = 0.0_r8
         551           0 :               else if (rh(i,k) < rhdeliques(m)) then
         552           0 :                  wtrvol(i,k,m) = wtrvol(i,k,m)*hystfac*(rh(i,k) - rhcrystal(m))
         553           0 :                  wtrvol(i,k,m) = max(wtrvol(i,k,m), 0.0_r8)
         554           0 :                  wetvol(i,k,m) = dryvol(i,k,m) + wtrvol(i,k,m)
         555           0 :                  wetrad(i,k,m) = (wetvol(i,k,m)/pi43)**third
         556             :               end if
         557             :             end if
         558             :
         559             :          end do  ! columns
         560             :       end do     ! levels
         561             :
         562             :    end do ! modes
         563             :
         564           0 : end subroutine modal_aero_wateruptake_sub
         565             :
         566             : !-----------------------------------------------------------------------
         567           0 :       subroutine modal_aero_kohler(   &
         568           0 :           rdry_in, hygro, s, rwet_out, im )
         569             :
         570             : ! calculates equlibrium radius r of haze droplets as function of
         571             : ! dry particle mass and relative humidity s using kohler solution
         572             : ! given in pruppacher and klett (eqn 6-35)
         573             :
         574             : ! for multiple aerosol types, assumes an internal mixture of aerosols
         575             :
         576             :       implicit none
         577             :
         578             : ! arguments
         579             :       integer :: im         ! number of grid points to be processed
         580             :       real(r8) :: rdry_in(:)    ! aerosol dry radius (m)
         581             :       real(r8) :: hygro(:)      ! aerosol volume-mean hygroscopicity (--)
         582             :       real(r8) :: s(:)          ! relative humidity (1 = saturated)
         583             :       real(r8) :: rwet_out(:)   ! aerosol wet radius (m)
         584             :
         585             : ! local variables
         586             :       integer, parameter :: imax=200
         587             :       integer :: i, n, nsol
         588             :
         589             :       real(r8) :: a, b
         590             :       real(r8) :: p40(imax),p41(imax),p42(imax),p43(imax) ! coefficients of polynomial
         591             :       real(r8) :: p30(imax),p31(imax),p32(imax) ! coefficients of polynomial
         592             :       real(r8) :: p
         593             :       real(r8) :: r3, r4
         594           0 :       real(r8) :: r(im)         ! wet radius (microns)
         595             :       real(r8) :: rdry(imax)    ! radius of dry particle (microns)
         596             :       real(r8) :: ss            ! relative humidity (1 = saturated)
         597             :       real(r8) :: slog(imax)    ! log relative humidity
         598             :       real(r8) :: vol(imax)     ! total volume of particle (microns**3)
         599             :       real(r8) :: xi, xr
         600             :
         601             :       complex(r8) :: cx4(4,imax),cx3(3,imax)
         602             :
         603             :       real(r8), parameter :: eps = 1.e-4_r8
         604             :       real(r8), parameter :: mw = 18._r8
         605             :       real(r8), parameter :: pi = 3.14159_r8
         606             :       real(r8), parameter :: rhow = 1._r8
         607             :       real(r8), parameter :: surften = 76._r8
         608             :       real(r8), parameter :: tair = 273._r8
         609             :       real(r8), parameter :: third = 1._r8/3._r8
         610             :       real(r8), parameter :: ugascon = 8.3e7_r8
         611             :
         612             :
         613             : !     effect of organics on surface tension is neglected
         614           0 :       a=2.e4_r8*mw*surften/(ugascon*tair*rhow)
         615             :
         616           0 :       do i=1,im
         617           0 :            rdry(i) = rdry_in(i)*1.0e6_r8   ! convert (m) to (microns)
         618           0 :            vol(i) = rdry(i)**3          ! vol is r**3, not volume
         619           0 :            b = vol(i)*hygro(i)
         620             :
         621             : !          quartic
         622           0 :            ss=min(s(i),1._r8-eps)
         623           0 :            ss=max(ss,1.e-10_r8)
         624           0 :            slog(i)=log(ss)
         625           0 :            p43(i)=-a/slog(i)
         626           0 :            p42(i)=0._r8
         627           0 :            p41(i)=b/slog(i)-vol(i)
         628           0 :            p40(i)=a*vol(i)/slog(i)
         629             : !          cubic for rh=1
         630           0 :            p32(i)=0._r8
         631           0 :            p31(i)=-b/a
         632           0 :            p30(i)=-vol(i)
         633             :       end do
         634             :
         635             :
         636           0 :        do 100 i=1,im
         637             :
         638             : !       if(vol(i).le.1.e-20)then
         639           0 :         if(vol(i).le.1.e-12_r8)then
         640           0 :            r(i)=rdry(i)
         641           0 :            go to 100
         642             :         endif
         643             :
         644           0 :         p=abs(p31(i))/(rdry(i)*rdry(i))
         645           0 :         if(p.lt.eps)then
         646             : !          approximate solution for small particles
         647           0 :            r(i)=rdry(i)*(1._r8+p*third/(1._r8-slog(i)*rdry(i)/a))
         648             :         else
         649           0 :            call makoh_quartic(cx4(1,i),p43(i),p42(i),p41(i),p40(i),1)
         650             : !          find smallest real(r8) solution
         651           0 :            r(i)=1000._r8*rdry(i)
         652           0 :            nsol=0
         653           0 :            do n=1,4
         654           0 :               xr=real(cx4(n,i))
         655           0 :               xi=aimag(cx4(n,i))
         656           0 :               if(abs(xi).gt.abs(xr)*eps) cycle
         657           0 :               if(xr.gt.r(i)) cycle
         658           0 :               if(xr.lt.rdry(i)*(1._r8-eps)) cycle
         659           0 :               if(xr.ne.xr) cycle
         660           0 :               r(i)=xr
         661           0 :               nsol=n
         662             :            end do
         663           0 :            if(nsol.eq.0)then
         664             :               write(iulog,*)   &
         665           0 :                'ccm kohlerc - no real(r8) solution found (quartic)'
         666           0 :               write(iulog,*)'roots =', (cx4(n,i),n=1,4)
         667           0 :               write(iulog,*)'p0-p3 =', p40(i), p41(i), p42(i), p43(i)
         668           0 :               write(iulog,*)'rh=',s(i)
         669           0 :               write(iulog,*)'setting radius to dry radius=',rdry(i)
         670           0 :               r(i)=rdry(i)
         671             : !             stop
         672             :            endif
         673             :         endif
         674             :
         675           0 :         if(s(i).gt.1._r8-eps)then
         676             : !          save quartic solution at s=1-eps
         677           0 :            r4=r(i)
         678             : !          cubic for rh=1
         679           0 :            p=abs(p31(i))/(rdry(i)*rdry(i))
         680           0 :            if(p.lt.eps)then
         681           0 :               r(i)=rdry(i)*(1._r8+p*third)
         682             :            else
         683           0 :               call makoh_cubic(cx3,p32,p31,p30,im)
         684             : !             find smallest real(r8) solution
         685           0 :               r(i)=1000._r8*rdry(i)
         686           0 :               nsol=0
         687           0 :               do n=1,3
         688           0 :                  xr=real(cx3(n,i))
         689           0 :                  xi=aimag(cx3(n,i))
         690           0 :                  if(abs(xi).gt.abs(xr)*eps) cycle
         691           0 :                  if(xr.gt.r(i)) cycle
         692           0 :                  if(xr.lt.rdry(i)*(1._r8-eps)) cycle
         693           0 :                  if(xr.ne.xr) cycle
         694           0 :                  r(i)=xr
         695           0 :                  nsol=n
         696             :               end do
         697           0 :               if(nsol.eq.0)then
         698             :                  write(iulog,*)   &
         699           0 :                   'ccm kohlerc - no real(r8) solution found (cubic)'
         700           0 :                  write(iulog,*)'roots =', (cx3(n,i),n=1,3)
         701           0 :                  write(iulog,*)'p0-p2 =', p30(i), p31(i), p32(i)
         702           0 :                  write(iulog,*)'rh=',s(i)
         703           0 :                  write(iulog,*)'setting radius to dry radius=',rdry(i)
         704           0 :                  r(i)=rdry(i)
         705             : !                stop
         706             :               endif
         707             :            endif
         708           0 :            r3=r(i)
         709             : !          now interpolate between quartic, cubic solutions
         710           0 :            r(i)=(r4*(1._r8-s(i))+r3*(s(i)-1._r8+eps))/eps
         711             :         endif
         712             :
         713           0 :   100 continue
         714             :
         715             : ! bound and convert from microns to m
         716           0 :       do i=1,im
         717           0 :          r(i) = min(r(i),30._r8) ! upper bound based on 1 day lifetime
         718           0 :          rwet_out(i) = r(i)*1.e-6_r8
         719             :       end do
         720             :
         721           0 :       return
         722             :       end subroutine modal_aero_kohler
         723             :
         724             :
         725             : !-----------------------------------------------------------------------
         726           0 :       subroutine makoh_cubic( cx, p2, p1, p0, im )
         727             : !
         728             : !     solves  x**3 + p2 x**2 + p1 x + p0 = 0
         729             : !     where p0, p1, p2 are real
         730             : !
         731             :       integer, parameter :: imx=200
         732             :       integer :: im
         733             :       real(r8) :: p0(imx), p1(imx), p2(imx)
         734             :       complex(r8) :: cx(3,imx)
         735             :
         736             :       integer :: i
         737             :       real(r8) :: eps, q(imx), r(imx), sqrt3, third
         738             :       complex(r8) :: ci, cq, crad(imx), cw, cwsq, cy(imx), cz(imx)
         739             :
         740             :       save eps
         741             :       data eps/1.e-20_r8/
         742             :
         743           0 :       third=1._r8/3._r8
         744           0 :       ci=cmplx(0._r8,1._r8,r8)
         745           0 :       sqrt3=sqrt(3._r8)
         746           0 :       cw=0.5_r8*(-1+ci*sqrt3)
         747           0 :       cwsq=0.5_r8*(-1-ci*sqrt3)
         748             :
         749           0 :       do i=1,im
         750           0 :       if(p1(i).eq.0._r8)then
         751             : !        completely insoluble particle
         752           0 :          cx(1,i)=(-p0(i))**third
         753           0 :          cx(2,i)=cx(1,i)
         754           0 :          cx(3,i)=cx(1,i)
         755             :       else
         756           0 :          q(i)=p1(i)/3._r8
         757           0 :          r(i)=p0(i)/2._r8
         758           0 :          crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
         759           0 :          crad(i)=sqrt(crad(i))
         760             :
         761           0 :          cy(i)=r(i)-crad(i)
         762           0 :          if (abs(cy(i)).gt.eps) cy(i)=cy(i)**third
         763           0 :          cq=q(i)
         764           0 :          cz(i)=-cq/cy(i)
         765             :
         766           0 :          cx(1,i)=-cy(i)-cz(i)
         767           0 :          cx(2,i)=-cw*cy(i)-cwsq*cz(i)
         768           0 :          cx(3,i)=-cwsq*cy(i)-cw*cz(i)
         769             :       endif
         770             :       enddo
         771             :
         772           0 :       return
         773             :       end subroutine makoh_cubic
         774             :
         775             :
         776             : !-----------------------------------------------------------------------
         777           0 :       subroutine makoh_quartic( cx, p3, p2, p1, p0, im )
         778             :
         779             : !     solves x**4 + p3 x**3 + p2 x**2 + p1 x + p0 = 0
         780             : !     where p0, p1, p2, p3 are real
         781             : !
         782             :       integer, parameter :: imx=200
         783             :       integer :: im
         784             :       real(r8) :: p0(imx), p1(imx), p2(imx), p3(imx)
         785             :       complex(r8) :: cx(4,imx)
         786             :
         787             :       integer :: i
         788             :       real(r8) :: third, q(imx), r(imx)
         789             :       complex(r8) :: cb(imx), cb0(imx), cb1(imx),   &
         790             :                      crad(imx), cy(imx), czero
         791             :
         792             :
         793           0 :       czero=cmplx(0.0_r8,0.0_r8,r8)
         794           0 :       third=1._r8/3._r8
         795             :
         796           0 :       do 10 i=1,im
         797             :
         798           0 :       q(i)=-p2(i)*p2(i)/36._r8+(p3(i)*p1(i)-4*p0(i))/12._r8
         799             :       r(i)=-(p2(i)/6)**3+p2(i)*(p3(i)*p1(i)-4*p0(i))/48._r8   &
         800           0 :        +(4*p0(i)*p2(i)-p0(i)*p3(i)*p3(i)-p1(i)*p1(i))/16
         801             :
         802           0 :       crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
         803           0 :       crad(i)=sqrt(crad(i))
         804             :
         805           0 :       cb(i)=r(i)-crad(i)
         806           0 :       if(cb(i).eq.czero)then
         807             : !        insoluble particle
         808           0 :          cx(1,i)=(-p1(i))**third
         809           0 :          cx(2,i)=cx(1,i)
         810           0 :          cx(3,i)=cx(1,i)
         811           0 :          cx(4,i)=cx(1,i)
         812             :       else
         813           0 :          cb(i)=cb(i)**third
         814             :
         815           0 :          cy(i)=-cb(i)+q(i)/cb(i)+p2(i)/6
         816             :
         817           0 :          cb0(i)=sqrt(cy(i)*cy(i)-p0(i))
         818           0 :          cb1(i)=(p3(i)*cy(i)-p1(i))/(2*cb0(i))
         819             :
         820           0 :          cb(i)=p3(i)/2+cb1(i)
         821           0 :          crad(i)=cb(i)*cb(i)-4*(cy(i)+cb0(i))
         822           0 :          crad(i)=sqrt(crad(i))
         823           0 :          cx(1,i)=(-cb(i)+crad(i))/2._r8
         824           0 :          cx(2,i)=(-cb(i)-crad(i))/2._r8
         825             :
         826           0 :          cb(i)=p3(i)/2-cb1(i)
         827           0 :          crad(i)=cb(i)*cb(i)-4*(cy(i)-cb0(i))
         828           0 :          crad(i)=sqrt(crad(i))
         829           0 :          cx(3,i)=(-cb(i)+crad(i))/2._r8
         830           0 :          cx(4,i)=(-cb(i)-crad(i))/2._r8
         831             :       endif
         832           0 :    10 continue
         833             :
         834           0 :       return
         835             :       end subroutine makoh_quartic
         836             :
         837             : !----------------------------------------------------------------------
         838           0 :       subroutine calc_h2so4_equilib_mixrat( temp, pres, qh2o, dmean, &
         839             :                                             qh2so4_equilib, wtpct, sulden )
         840             :
         841             :       implicit none
         842             :
         843             :       real(r8), intent(in)  :: temp            ! temperature (K)
         844             :       real(r8), intent(in)  :: pres            ! pressure (Pa)
         845             :       real(r8), intent(in)  :: qh2o            ! water vapor specific humidity (kg/kg)
         846             :       real(r8), intent(in)  :: dmean           ! mean diameter of particles in each mode
         847             :       real(r8), intent(out) :: qh2so4_equilib  ! h2so4 saturation mixing ratios over the particles (mol/mol)
         848             :       real(r8), intent(out) :: wtpct           ! sulfate composition, weight % H2SO4
         849             :       real(r8), intent(out) :: sulden          ! sulfate aerosol mass density (g/cm3)
         850             :
         851             :       ! Local declarations
         852             :       real(r8)            :: qh2o_kelvin ! water vapor specific humidity adjusted for Kelvin effect (kg/kg)
         853             :       real(r8)            :: wtpct_flat  ! sulfate composition over a flat surface, weight % H2SO4
         854             :       real(r8)            :: fk1, fk4, fk4_1, fk4_2
         855             :       real(r8)            :: factor_kulm                     ! Kulmala correction terms
         856             :       real(r8)            :: en, t, sig1, sig2, frac, surf_tens, surf_tens_mode, dsigma_dwt
         857             :       real(r8)            :: den1, den2, sulfate_density, drho_dwt
         858             :       real(r8)            :: akelvin, expon, akas, rkelvinH2O, rkelvinH2O_a, rkelvinH2O_b
         859             :       real(r8)            :: sulfequil, r
         860             :       real(r8), parameter :: t0_kulm     = 340._r8           !  T0 set in the low end of the Ayers measurement range (338-445K)
         861             :       real(r8), parameter :: t_crit_kulm = 905._r8           !  Critical temperature = 1.5 * Boiling point
         862             :       real(r8), parameter :: fk0         = -10156._r8 / t0_kulm + 16.259_r8   !  Log(Kulmala correction factor)
         863             :       real(r8), parameter :: fk2         = 1._r8 / t0_kulm
         864             :       real(r8), parameter :: fk3         = 0.38_r8 / (t_crit_kulm - t0_kulm)
         865             :       real(r8), parameter :: RGAS = 8.31430e+07_r8           ! ideal gas constant (erg/mol/K)
         866             :       real(r8), parameter :: wtmol_h2so4 =  98.078479_r8     ! molecular weight of sulphuric acid
         867             :       real(r8), parameter :: wtmol_h2o   =  18.015280_r8     ! molecular weight of water vapor
         868             :       real(r8), parameter :: wtmol_air   =  28.9644_r8       ! molecular weight of dry air
         869             :       real(r8)            :: gv_wt_abcd_en(6,95), gvh2ovp(95)
         870             :       real(r8)            :: dnwtp(46), dnc0(46), dnc1(46)
         871             :       real(r8)            :: stwtp(15), stc0(15), stc1(15)
         872             :       integer             :: i, k
         873             :
         874             :
         875             :       data stwtp/0._r8, 23.8141_r8, 38.0279_r8, 40.6856_r8, 45.335_r8, 52.9305_r8, &
         876             :          56.2735_r8, 59.8557_r8, 66.2364_r8, 73.103_r8, 79.432_r8, 85.9195_r8, &
         877             :          91.7444_r8, 97.6687_r8, 100._r8/
         878             :
         879             :       data stc0/117.564_r8, 103.303_r8, 101.796_r8, 100.42_r8, 98.4993_r8, 91.8866_r8, &
         880             :          88.3033_r8, 86.5546_r8, 84.471_r8, 81.2939_r8, 79.3556_r8, 75.608_r8, &
         881             :          70.0777_r8, 63.7412_r8, 61.4591_r8 /
         882             :
         883             :       data stc1/-0.153641_r8, -0.0982007_r8, -0.0872379_r8, -0.0818509_r8,           &
         884             :          -0.0746702_r8, -0.0522399_r8, -0.0407773_r8, -0.0357946_r8, -0.0317062_r8,   &
         885             :          -0.025825_r8, -0.0267212_r8, -0.0269204_r8, -0.0276187_r8, -0.0302094_r8,    &
         886             :          -0.0303081_r8 /
         887             :
         888             :
         889             :       data dnwtp / 0._r8, 1._r8, 5._r8, 10._r8, 20._r8, 25._r8, 30._r8, 35._r8, 40._r8,  &
         890             :          41._r8, 45._r8, 50._r8, 53._r8, 55._r8, 56._r8, 60._r8, 65._r8, 66._r8, 70._r8, &
         891             :          72._r8, 73._r8, 74._r8, 75._r8, 76._r8, 78._r8, 79._r8, 80._r8, 81._r8, 82._r8, &
         892             :          83._r8, 84._r8, 85._r8, 86._r8, 87._r8, 88._r8, 89._r8, 90._r8, 91._r8, 92._r8, &
         893             :          93._r8, 94._r8, 95._r8, 96._r8, 97._r8, 98._r8, 100._r8 /
         894             :
         895             :       data dnc0 / 1._r8, 1.13185_r8, 1.17171_r8, 1.22164_r8, 1.3219_r8, 1.37209_r8,         &
         896             :          1.42185_r8, 1.4705_r8, 1.51767_r8, 1.52731_r8, 1.56584_r8, 1.61834_r8, 1.65191_r8, &
         897             :          1.6752_r8, 1.68708_r8, 1.7356_r8, 1.7997_r8, 1.81271_r8, 1.86696_r8, 1.89491_r8,   &
         898             :          1.9092_r8, 1.92395_r8, 1.93904_r8, 1.95438_r8, 1.98574_r8, 2.00151_r8, 2.01703_r8, &
         899             :          2.03234_r8, 2.04716_r8, 2.06082_r8, 2.07363_r8, 2.08461_r8, 2.09386_r8, 2.10143_r8,&
         900             :          2.10764_r8, 2.11283_r8, 2.11671_r8, 2.11938_r8, 2.12125_r8, 2.1219_r8, 2.12723_r8, &
         901             :          2.12654_r8, 2.12621_r8, 2.12561_r8, 2.12494_r8, 2.12093_r8 /
         902             :
         903             :       data dnc1 / 0._r8,  -0.000435022_r8, -0.000479481_r8, -0.000531558_r8, -0.000622448_r8, &
         904             :          -0.000660866_r8, -0.000693492_r8, -0.000718251_r8, -0.000732869_r8, -0.000735755_r8, &
         905             :          -0.000744294_r8, -0.000761493_r8, -0.000774238_r8, -0.00078392_r8, -0.000788939_r8,  &
         906             :          -0.00080946_r8, -0.000839848_r8, -0.000845825_r8, -0.000874337_r8, -0.000890074_r8,  &
         907             :          -0.00089873_r8, -0.000908778_r8, -0.000920012_r8, -0.000932184_r8, -0.000959514_r8,  &
         908             :          -0.000974043_r8, -0.000988264_r8, -0.00100258_r8, -0.00101634_r8, -0.00102762_r8,    &
         909             :          -0.00103757_r8, -0.00104337_r8, -0.00104563_r8, -0.00104458_r8, -0.00104144_r8,      &
         910             :          -0.00103719_r8, -0.00103089_r8, -0.00102262_r8, -0.00101355_r8, -0.00100249_r8,      &
         911             :          -0.00100934_r8, -0.000998299_r8, -0.000990961_r8, -0.000985845_r8, -0.000984529_r8,  &
         912             :          -0.000989315_r8 /
         913             :
         914             :       ! Saturation vapor pressure of sulfuric acid
         915             :       !
         916             :       ! Limit extrapolation at extreme temperatures
         917           0 :       t=min(max(temp,140._r8),450._r8)
         918             :
         919             :       !!  Calculate the weight % H2SO4 composition of sulfate
         920           0 :       call calc_h2so4_wtpct(t, pres, qh2o, wtpct_flat)
         921             :
         922             :       !!  Calculate surface tension (erg/cm2) of sulfate of
         923             :       !!  different compositions as a linear function of temperature.
         924           0 :       i = 2 ! minimum wt% is 29.517
         925           0 :       do while (wtpct_flat.gt.stwtp(i))
         926           0 :        i = i + 1
         927             :       end do
         928           0 :       sig1 = stc0(i-1) + stc1(i-1) * t
         929           0 :       sig2 = stc0(i)   + stc1(i) * t
         930             :       ! calculate derivative needed later for kelvin factor for h2o
         931           0 :       dsigma_dwt = (sig2-sig1) / (stwtp(i)-stwtp(i-1))
         932           0 :       surf_tens = sig1 + dsigma_dwt*(wtpct_flat-stwtp(i))
         933             :
         934             :       !!  Calculate density (g/cm3) of sulfate of
         935             :       !!  different compositions as a linear function of temperature.
         936           0 :       i = 6 ! minimum wt% is 29.517
         937           0 :       do while (wtpct_flat .gt. dnwtp(i))
         938           0 :         i = i + 1
         939             :       end do
         940           0 :       den1 = dnc0(i-1) + dnc1(i-1) * t
         941           0 :       den2 = dnc0(i)   + dnc1(i) * t
         942             :       ! Calculate derivative needed later for Kelvin factor for H2O
         943           0 :       drho_dwt = (den2-den1) / (dnwtp(i)-dnwtp(i-1))
         944           0 :       sulfate_density = den1 + drho_dwt * (wtpct_flat-dnwtp(i-1))
         945             :
         946           0 :       r=dmean*100._r8/2._r8 ! calcuate mode radius (cm) from diameter (m)
         947             :
         948             :       ! Adjust for Kelvin effect for water
         949             :       rkelvinH2O_b = 1 + wtpct_flat * drho_dwt / &
         950           0 :          sulfate_density - 3._r8 * wtpct_flat * dsigma_dwt / (2._r8*surf_tens)
         951             :
         952             :       rkelvinH2O_a = 2._r8 * wtmol_h2so4 * surf_tens / &
         953           0 :            (sulfate_density * RGAS * t * r)
         954             :
         955           0 :       rkelvinH2O = exp (rkelvinH2O_a*rkelvinH2O_b)
         956             :
         957           0 :       qh2o_kelvin = qh2o/rkelvinH2O
         958           0 :       call calc_h2so4_wtpct(t, pres, qh2o_kelvin, wtpct)
         959             :
         960             :
         961           0 :       wtpct=max(wtpct,wtpct_flat)
         962             :
         963             :       ! Parameterized fit to Giauque's (1959) enthalpies v. wt %:
         964           0 :       en = 4.184_r8 * (23624.8_r8 - 1.14208e8_r8 / ((wtpct - 105.318_r8)**2 + 4798.69_r8))
         965           0 :       en = max(en, 0.0_r8)
         966             :
         967             :       !!  Calculate surface tension (erg/cm2) of sulfate of
         968             :       !!  different compositions as a linear function of temperature.
         969           0 :       i = 2 ! minimum wt% is 29.517
         970           0 :       do while (wtpct.gt.stwtp(i))
         971           0 :        i=i+1
         972             :       end do
         973           0 :       sig1=stc0(i-1)+stc1(i-1)*t
         974           0 :       sig2=stc0(i)+stc1(i)*t
         975           0 :       frac=(stwtp(i)-wtpct)/(stwtp(i)-stwtp(i-1))
         976           0 :       surf_tens_mode=sig1*frac+sig2*(1.0_r8-frac)
         977             :
         978             :       !!  Calculate density (g/cm3) of sulfate of
         979             :       !!  different compositions as a linear function of temperature.
         980           0 :       i = 6 ! minimum wt% is 29.517
         981           0 :       do while (wtpct .gt. dnwtp(i))
         982           0 :         i=i+1
         983             :       end do
         984           0 :       den1=dnc0(i-1)+dnc1(i-1)*t
         985           0 :       den2=dnc0(i)+dnc1(i)*t
         986           0 :       frac=(dnwtp(i)-wtpct)/(dnwtp(i)-dnwtp(i-1))
         987           0 :       sulden=den1*frac+den2*(1.0_r8-frac)
         988             :
         989             :       ! Ayers' (1980) fit to sulfuric acid equilibrium vapor pressure:
         990             :       ! (Remember this is the log)
         991             :       ! SULFEQ=-10156/Temp+16.259-En/(8.3143*Temp)
         992             :       !
         993             :       ! Kulmala correction (J. CHEM. PHYS. V.93, No.1, 1 July 1990)
         994           0 :       fk1   = -1._r8 / t
         995           0 :       fk4_1 = log(t0_kulm / t)
         996           0 :       fk4_2 = t0_kulm / t
         997           0 :       fk4   = 1.0_r8 + fk4_1 - fk4_2
         998           0 :       factor_kulm = fk1 + fk2 + fk3 * fk4
         999             :
        1000             :       ! This is for pure H2SO4
        1001           0 :       sulfequil = fk0 + 10156._r8 * factor_kulm
        1002             :
        1003             :       ! Adjust for WTPCT composition:
        1004           0 :       sulfequil = sulfequil - en / (8.3143_r8 * t)
        1005             :
        1006             :       ! Take the exponential:
        1007           0 :       sulfequil = exp(sulfequil)
        1008             :
        1009             :       ! Convert atmospheres ==> Pa
        1010           0 :       sulfequil = sulfequil * 1.01325e5_r8
        1011             :
        1012             :       ! Convert Pa ==> mol/mol
        1013           0 :       sulfequil = sulfequil / pres
        1014             :
        1015             :       ! Calculate Kelvin curvature factor for H2SO4 interactively with temperature:
        1016             :       ! (g/mol)*(erg/cm2)/(K * g/cm3 * erg/mol/K) = cm
        1017           0 :       akelvin = 2._r8 * wtmol_h2so4 * surf_tens_mode / (t * sulden * RGAS)
        1018             :
        1019           0 :       expon = akelvin / r  ! divide by mode radius (cm)
        1020           0 :       expon = max(-100._r8, expon)
        1021           0 :       expon = min(100._r8, expon)
        1022           0 :       akas = exp( expon )
        1023           0 :       qh2so4_equilib = sulfequil * akas ! reduce H2SO4 equilibrium mixing ratio by Kelvin curvature factor
        1024             :
        1025           0 :       return
        1026             :       end subroutine calc_h2so4_equilib_mixrat
        1027             :
        1028             :
        1029             : !----------------------------------------------------------------------
        1030           0 :       subroutine calc_h2so4_wtpct( temp, pres, qh2o, wtpct )
        1031             :
        1032             :   !!  This function calculates the weight % H2SO4 composition of
        1033             :   !!  sulfate aerosol, using Tabazadeh et. al. (GRL, 1931, 1997).
        1034             :   !!  Rated for T=185-260K, activity=0.01-1.0
        1035             :   !!
        1036             :   !!  Argument list input:
        1037             :   !!    temp = temperature (K)
        1038             :   !!    pres = atmospheric pressure (Pa)
        1039             :   !!    qh2o = water specific humidity (kg/kg)
        1040             :   !!
        1041             :   !!  Output:
        1042             :   !!    wtpct = weight % H2SO4 in H2O/H2SO4 particle (0-100)
        1043             :   !!
        1044             :   !! @author Mike Mills
        1045             :   !! @ version October 2013
        1046             :
        1047             :       use wv_saturation, only: qsat_water
        1048             :
        1049             :       implicit none
        1050             :
        1051             :       real(r8), intent(in)  :: temp  ! temperature (K)
        1052             :       real(r8), intent(in)  :: pres  ! pressure (Pa)
        1053             :       real(r8), intent(in)  :: qh2o  ! water vapor specific humidity (kg/kg)
        1054             :       real(r8), intent(out) :: wtpct ! sulfate weight % H2SO4 composition
        1055             :
        1056             :       ! Local declarations
        1057             :       real(r8)            :: atab1,btab1,ctab1,dtab1,atab2,btab2,ctab2,dtab2
        1058             :       real(r8)            :: contl, conth, contt, conwtp
        1059             :       real(r8)            :: activ
        1060             :       real(r8)            :: es ! saturation vapor pressure over water (Pa) (dummy)
        1061             :       real(r8)            :: qs ! saturation specific humidity over water (kg/kg)
        1062             :
        1063             :       ! calculate saturation specific humidity over pure water, qs (kg/kg)
        1064           0 :       call qsat_water(temp, pres, es, qs)
        1065             :
        1066             :       !  Activity = water specific humidity (kg/kg) / equilibrium water (kg/kg)
        1067           0 :       activ = qh2o/qs
        1068             :
        1069           0 :       if (activ.lt.0.05_r8) then
        1070           0 :         activ = max(activ,1.e-6_r8)    ! restrict minimum activity
        1071           0 :         atab1   = 12.37208932_r8
        1072           0 :         btab1   = -0.16125516114_r8
        1073           0 :         ctab1   = -30.490657554_r8
        1074           0 :         dtab1   = -2.1133114241_r8
        1075           0 :         atab2   = 13.455394705_r8
        1076           0 :         btab2   = -0.1921312255_r8
        1077           0 :         ctab2   = -34.285174607_r8
        1078           0 :         dtab2   = -1.7620073078_r8
        1079           0 :       elseif (activ.ge.0.05_r8.and.activ.le.0.85_r8) then
        1080             :         atab1   = 11.820654354_r8
        1081             :         btab1   = -0.20786404244_r8
        1082             :         ctab1   = -4.807306373_r8
        1083             :         dtab1   = -5.1727540348_r8
        1084             :         atab2   = 12.891938068_r8
        1085             :         btab2   = -0.23233847708_r8
        1086             :         ctab2   = -6.4261237757_r8
        1087             :         dtab2   = -4.9005471319_r8
        1088           0 :       elseif (activ.gt.0.85_r8) then
        1089           0 :         activ = min(activ,1._r8)      ! restrict maximum activity
        1090           0 :         atab1   = -180.06541028_r8
        1091           0 :         btab1   = -0.38601102592_r8
        1092           0 :         ctab1   = -93.317846778_r8
        1093           0 :         dtab1   = 273.88132245_r8
        1094           0 :         atab2   = -176.95814097_r8
        1095           0 :         btab2   = -0.36257048154_r8
        1096           0 :         ctab2   = -90.469744201_r8
        1097           0 :         dtab2   = 267.45509988_r8
        1098             :       else
        1099           0 :         write(iulog,*) 'calc_h2so4_wtpct: invalid activity: activ,qh2o,qs,temp,pres=',activ,qh2o,qs,temp,pres
        1100           0 :         call endrun( 'calc_h2so4_wtpct error' )
        1101           0 :         return
        1102             :       endif
        1103             :
        1104           0 :       contl = atab1*(activ**btab1)+ctab1*activ+dtab1
        1105           0 :       conth = atab2*(activ**btab2)+ctab2*activ+dtab2
        1106             :
        1107           0 :       contt = contl + (conth-contl) * ((temp -190._r8)/70._r8)
        1108           0 :       conwtp = (contt*98._r8) + 1000._r8
        1109             :
        1110           0 :       wtpct = (100._r8*contt*98._r8)/conwtp
        1111           0 :       wtpct = min(max(wtpct,25._r8),100._r8) ! restrict between 1 and 100 %
        1112             :
        1113           0 :       return
        1114             :       end subroutine calc_h2so4_wtpct
        1115             :
        1116             :
        1117             : !----------------------------------------------------------------------
        1118             :
        1119             :    end module modal_aero_wateruptake
        1120             :
        1121             :

Generated by: LCOV version 1.14