LCOV - code coverage report
Current view: | top level - physics/cam - phys_prop.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 10 | 487 | 2.1 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 2 | 20 | 10.0 %

              Line data    Source code
           1             : module phys_prop
           2             :
           3             : ! Properties of aerosols that are used by radiation and other parameterizations.
           4             :
           5             : ! *****N.B.*****
           6             : ! This module is a utility used by the rad_constituents module.  The properties stored
           7             : ! here are meant to be accessed via that module.  This module knows nothing about how
           8             : ! this data is associated with the constituents that are radiatively active or those that
           9             : ! are being used for diagnostic calculations.  That is the responsibility of the
          10             : ! rad_constituents module.
          11             :
          12             : use shr_kind_mod,   only: r8 => shr_kind_r8
          13             : use spmd_utils,     only: masterproc
          14             : use radconstants,   only: nlwbands, nswbands, idx_sw_diag
          15             : use ioFileMod,      only: getfil
          16             : use cam_pio_utils,  only: cam_pio_openfile
          17             : use pio,            only: file_desc_t, var_desc_t, pio_get_var, pio_inq_varid, &
          18             :                           pio_inq_dimlen, pio_inq_dimid , pio_nowrite, pio_closefile, &
          19             :                           pio_seterrorhandling, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR, PIO_NOERR
          20             :
          21             : use cam_logfile,    only: iulog
          22             : use cam_abortutils, only: endrun
          23             :
          24             : implicit none
          25             : private
          26             : save
          27             :
          28             : integer, parameter, public :: ot_length = 32
          29             :
          30             : public :: &
          31             :    physprop_accum_unique_files,  &! Make a list of the unique set of files that contain properties
          32             :                                   ! This is an initialization step that must be done before calling physprop_init
          33             :    physprop_init,                &! Initialization -- read the input datasets
          34             :    physprop_get_id,              &! Return ID used to access the property data from the input files
          35             :    physprop_get                   ! Return data for specified ID
          36             :
          37             : ! Data from one input dataset is stored in a structure of type(physprop_type).
          38             : type :: physprop_type
          39             :    character(len=256) :: sourcefile ! Absolute pathname of data file.
          40             :    character(len=ot_length)  :: opticsmethod ! one of {hygro,nonhygro}
          41             :
          42             :    ! for hygroscopic species of externally mixed aerosols
          43             :    real(r8), pointer :: sw_hygro_ext(:,:)
          44             :    real(r8), pointer :: sw_hygro_ssa(:,:)
          45             :    real(r8), pointer :: sw_hygro_asm(:,:)
          46             :    real(r8), pointer :: lw_hygro_abs(:,:)
          47             :
          48             :    ! for nonhygroscopic species of externally mixed aerosols
          49             :    real(r8), pointer :: sw_nonhygro_ext(:)
          50             :    real(r8), pointer :: sw_nonhygro_ssa(:)
          51             :    real(r8), pointer :: sw_nonhygro_asm(:)
          52             :    real(r8), pointer :: sw_nonhygro_scat(:)
          53             :    real(r8), pointer :: sw_nonhygro_ascat(:)
          54             :    real(r8), pointer :: lw_abs(:)
          55             :
          56             :    ! complex refractive index
          57             :    complex(r8), pointer :: refindex_aer_sw(:)
          58             :    complex(r8), pointer :: refindex_aer_lw(:)
          59             :
          60             :    ! for radius-dependent mass-specific quantities
          61             :    real(r8), pointer :: r_sw_ext(:,:)
          62             :    real(r8), pointer :: r_sw_scat(:,:)
          63             :    real(r8), pointer :: r_sw_ascat(:,:)
          64             :    real(r8), pointer :: r_lw_abs(:,:)
          65             :    real(r8), pointer :: mu(:)
          66             :
          67             :    ! for modal optics
          68             :    real(r8), pointer :: extpsw(:,:,:,:) ! specific extinction
          69             :    real(r8), pointer :: abspsw(:,:,:,:) ! specific absorption
          70             :    real(r8), pointer :: asmpsw(:,:,:,:) ! asymmetry factor
          71             :    real(r8), pointer :: absplw(:,:,:,:) ! specific absorption
          72             :    real(r8), pointer :: refrtabsw(:,:)  ! table of real refractive indices for aerosols visible
          73             :    real(r8), pointer :: refitabsw(:,:)  ! table of imag refractive indices for aerosols visible
          74             :    real(r8), pointer :: refrtablw(:,:)  ! table of real refractive indices for aerosols infrared
          75             :    real(r8), pointer :: refitablw(:,:)  ! table of imag refractive indices for aerosols infrared
          76             :
          77             :    ! microphysics parameters.
          78             :    character(len=32) :: aername ! for output of number concentration
          79             :    real(r8) :: density_aer      ! density of aerosol (kg/m3)
          80             :    real(r8) :: hygro_aer        ! hygroscopicity of aerosol
          81             :    real(r8) :: dryrad_aer       ! number mode radius (m) of aerosol size distribution
          82             :    real(r8) :: dispersion_aer   ! geometric standard deviation of aerosol size distribution
          83             :    real(r8) :: num_to_mass_aer  ! ratio of number concentration to mass concentration (#/kg)
          84             :                                 ! *** Is this actually (kg/#) ???
          85             :    ! mode parameters
          86             :    integer :: ncoef       ! number of Chebyshev coefficients
          87             :    integer :: prefr       ! dimension in table of real refractive indices
          88             :    integer :: prefi       ! dimension in table of imag refractive indices
          89             :    real(r8) :: sigmag     ! geometric standard deviation of the number distribution for aerosol mode
          90             :    real(r8) :: dgnum      ! geometric dry mean diameter of the number distribution for aerosol mode
          91             :    real(r8) :: dgnumlo    ! lower limit of dgnum
          92             :    real(r8) :: dgnumhi    ! upper limit of dgnum
          93             :    real(r8) :: rhcrystal  ! crystalization relative humidity for mode
          94             :    real(r8) :: rhdeliques ! deliquescence relative humidity for mode
          95             :
          96             : endtype physprop_type
          97             :
          98             : ! This module stores data in an array of physprop_type structures.  The way this data
          99             : ! is accessed outside the module is via a physprop ID, which is an index into the array.
         100             : integer :: numphysprops = 0 ! an incremental total across ALL clim and diag constituents
         101             : type (physprop_type), pointer :: physprop(:)
         102             :
         103             : ! Temporary storage location for filenames in namelist, and construction of dynamic index
         104             : ! to properties.  The unique filenames specified in the namelist are the identifiers of
         105             : ! the properties.  Searching the uniquefilenames array provides the index into the physprop
         106             : ! array.
         107             : character(len=256), allocatable :: uniquefilenames(:)
         108             :
         109             : ! Number of evenly spaced intervals in rh used in this module and in the aer_rad_props module
         110             : ! for calculations of aerosol hygroscopic growth.
         111             : integer, parameter, public :: nrh = 1000
         112             :
         113             : !================================================================================================
         114             : contains
         115             : !================================================================================================
         116             :
         117        1536 : subroutine physprop_accum_unique_files(radname, type)
         118             :
         119             :    ! Count number of aerosols in input radname array.  Aerosols are identified
         120             :    ! as strings with a ".nc" suffix.
         121             :    ! Construct a cumulative list of unique filenames containing physical property data.
         122             :
         123             :    character(len=*), intent(in)  :: radname(:)
         124             :    character(len=1), intent(in)  :: type(:)
         125             :
         126             :    integer :: ncnst, i
         127             :    character(len=*), parameter :: subname = 'physprop_accum_unique_files'
         128             :    !------------------------------------------------------------------------------------
         129             :
         130             :    ! Initial guess for number of files we need.
         131        1536 :    if (.not. allocated(uniquefilenames)) allocate(uniquefilenames(50))
         132             :
         133        1536 :    ncnst = ubound(radname, 1)
         134             :
         135       13824 :    do i = 1, ncnst
         136             :
         137             :       ! check if radname is either a bulk aerosol or a mode
         138       13824 :       if (type(i) == 'A' .or. type(i) == 'M') then
         139             :
         140             :          ! check if this filename has been used by another aerosol.  If not
         141             :          ! then add it to the list of unique names.
         142           0 :          if (physprop_get_id(radname(i)) < 0) then
         143           0 :             numphysprops = numphysprops + 1
         144           0 :             if (numphysprops > size(uniquefilenames)) then
         145           0 :                call double_capacity(uniquefilenames)
         146             :             end if
         147           0 :             uniquefilenames(numphysprops) = trim(radname(i))
         148             :          endif
         149             :
         150             :       endif
         151             :    enddo
         152             :
         153             :  contains
         154             :
         155             :    ! Simple routine to re-allocate an array with twice the size, but with
         156             :    ! the inital values being preserved.
         157           0 :    subroutine double_capacity(array)
         158             :      character(len=256), intent(inout), allocatable :: array(:)
         159             :
         160           0 :      character(len=256), allocatable :: tmp(:)
         161             :      integer :: ierr
         162             :
         163           0 :      allocate(tmp(size(array)*2), stat=ierr)
         164           0 :      if ( ierr /= 0 ) then
         165           0 :         call endrun('physprop_accum_unique_files: Allocation error.')
         166             :      end if
         167             :
         168           0 :      tmp(:size(array)) = array
         169             :
         170           0 :      deallocate(array, stat=ierr)
         171             :      if ( ierr /= 0 ) then
         172           0 :         call endrun('physprop_accum_unique_files: Deallocation error.')
         173             :      end if
         174             :
         175           0 :      call move_alloc(tmp, array)
         176             :
         177           0 :    end subroutine double_capacity
         178             :
         179             : end subroutine physprop_accum_unique_files
         180             :
         181             : !================================================================================================
         182             :
         183        1536 : subroutine physprop_init()
         184             :
         185             :    ! Read properties from the aerosol data files.
         186             :
         187             :    ! ***N.B.*** The calls to physprop_accum_unique_files must be made before calling
         188             :    !            this init routine.  physprop_accum_unique_files is responsible for building
         189             :    !            the list of files to be read here.
         190             :
         191             :    ! Local variables
         192             :    integer            :: fileindex
         193             :    type(file_desc_t)  :: nc_id ! index to netcdf file
         194             :    character(len=256) :: locfn ! path to actual file used
         195             :    character(len=32)  :: aername_str ! string read from netCDF file -- may contain trailing
         196             :                                      ! nulls which aren't dealt with by trim()
         197             :
         198             :    integer :: ierr ! error codes from mpi
         199             :
         200             :    !------------------------------------------------------------------------------------
         201             :
         202        3072 :    allocate(physprop(numphysprops))
         203             :
         204        1536 :    do fileindex = 1, numphysprops
         205           0 :       nullify(physprop(fileindex)%sw_hygro_ext)
         206           0 :       nullify(physprop(fileindex)%sw_hygro_ssa)
         207           0 :       nullify(physprop(fileindex)%sw_hygro_asm)
         208           0 :       nullify(physprop(fileindex)%lw_hygro_abs)
         209             :
         210           0 :       nullify(physprop(fileindex)%sw_nonhygro_ext)
         211           0 :       nullify(physprop(fileindex)%sw_nonhygro_ssa)
         212           0 :       nullify(physprop(fileindex)%sw_nonhygro_asm)
         213           0 :       nullify(physprop(fileindex)%sw_nonhygro_scat)
         214           0 :       nullify(physprop(fileindex)%sw_nonhygro_ascat)
         215           0 :       nullify(physprop(fileindex)%lw_abs)
         216             :
         217           0 :       nullify(physprop(fileindex)%refindex_aer_sw)
         218           0 :       nullify(physprop(fileindex)%refindex_aer_lw)
         219             :
         220           0 :       nullify(physprop(fileindex)%r_sw_ext)
         221           0 :       nullify(physprop(fileindex)%r_sw_scat)
         222           0 :       nullify(physprop(fileindex)%r_sw_ascat)
         223           0 :       nullify(physprop(fileindex)%r_lw_abs)
         224           0 :       nullify(physprop(fileindex)%mu)
         225             :
         226           0 :       nullify(physprop(fileindex)%extpsw)
         227           0 :       nullify(physprop(fileindex)%abspsw)
         228           0 :       nullify(physprop(fileindex)%asmpsw)
         229           0 :       nullify(physprop(fileindex)%absplw)
         230           0 :       nullify(physprop(fileindex)%refrtabsw)
         231           0 :       nullify(physprop(fileindex)%refitabsw)
         232           0 :       nullify(physprop(fileindex)%refrtablw)
         233           0 :       nullify(physprop(fileindex)%refitablw)
         234             :
         235           0 :       call getfil(uniquefilenames(fileindex), locfn, 0)
         236           0 :       physprop(fileindex)%sourcefile = locfn
         237             :
         238             :       ! Open the physprop file
         239           0 :       call cam_pio_openfile(nc_id, locfn, PIO_NOWRITE)
         240             :
         241           0 :       call aerosol_optics_init(physprop(fileindex), nc_id)
         242             :
         243             :       ! Close the physprop file
         244        1536 :       call pio_closefile(nc_id)
         245             :
         246             :    end do
         247        1536 : end subroutine physprop_init
         248             :
         249             : !================================================================================================
         250             :
         251           0 : integer function physprop_get_id(filename)
         252             :
         253             :    ! Look for filename in the global list of unique filenames (module data uniquefilenames).
         254             :    ! If found, return it's index in the list.  Otherwise return -1.
         255             :
         256             :    character(len=*), intent(in) :: filename
         257             :    integer iphysprop
         258             :
         259           0 :    physprop_get_id = -1
         260           0 :    do iphysprop = 1, numphysprops
         261           0 :      if(trim(uniquefilenames(iphysprop)) == trim(filename) ) then
         262           0 :        physprop_get_id = iphysprop
         263           0 :        return
         264             :      endif
         265             :    enddo
         266             :
         267             : end function physprop_get_id
         268             :
         269             : !================================================================================================
         270             :
         271           0 : subroutine physprop_get(id, sourcefile, opticstype, &
         272             :    sw_hygro_ext, sw_hygro_ssa, sw_hygro_asm, lw_hygro_abs, &
         273             :    sw_nonhygro_ext, sw_nonhygro_ssa, sw_nonhygro_asm, &
         274             :    sw_nonhygro_scat, sw_nonhygro_ascat, lw_abs, &
         275             :    refindex_aer_sw, refindex_aer_lw, &
         276             :    r_sw_ext, r_sw_scat, r_sw_ascat, r_lw_abs, mu, &
         277             :    extpsw, abspsw, asmpsw, absplw, refrtabsw, &
         278             :    refitabsw, refrtablw, refitablw, &
         279           0 :    aername, density_aer, hygro_aer, dryrad_aer, dispersion_aer, &
         280             :    num_to_mass_aer, ncoef, prefr, prefi, sigmag, &
         281             :    dgnum, dgnumlo, dgnumhi, rhcrystal, rhdeliques)
         282             :
         283             :    ! Return requested properties for specified ID.
         284             :
         285             :    ! Arguments
         286             :    integer,                            intent(in)  :: id
         287             :    character(len=256),       optional, intent(out) :: sourcefile ! Absolute pathname of data file.
         288             :    character(len=ot_length), optional, intent(out) :: opticstype
         289             :    real(r8),          optional, pointer     :: sw_hygro_ext(:,:)
         290             :    real(r8),          optional, pointer     :: sw_hygro_ssa(:,:)
         291             :    real(r8),          optional, pointer     :: sw_hygro_asm(:,:)
         292             :    real(r8),          optional, pointer     :: lw_hygro_abs(:,:)
         293             :    real(r8),          optional, pointer     :: sw_nonhygro_ext(:)
         294             :    real(r8),          optional, pointer     :: sw_nonhygro_ssa(:)
         295             :    real(r8),          optional, pointer     :: sw_nonhygro_asm(:)
         296             :    real(r8),          optional, pointer     :: sw_nonhygro_scat(:)
         297             :    real(r8),          optional, pointer     :: sw_nonhygro_ascat(:)
         298             :    real(r8),          optional, pointer     :: lw_abs(:)
         299             :    complex(r8),       optional, pointer     :: refindex_aer_sw(:)
         300             :    complex(r8),       optional, pointer     :: refindex_aer_lw(:)
         301             :    real(r8),          optional, pointer     :: r_sw_ext(:,:)
         302             :    real(r8),          optional, pointer     :: r_sw_scat(:,:)
         303             :    real(r8),          optional, pointer     :: r_sw_ascat(:,:)
         304             :    real(r8),          optional, pointer     :: r_lw_abs(:,:)
         305             :    real(r8),          optional, pointer     :: mu(:)
         306             :    real(r8),          optional, pointer     :: extpsw(:,:,:,:)
         307             :    real(r8),          optional, pointer     :: abspsw(:,:,:,:)
         308             :    real(r8),          optional, pointer     :: asmpsw(:,:,:,:)
         309             :    real(r8),          optional, pointer     :: absplw(:,:,:,:)
         310             :    real(r8),          optional, pointer     :: refrtabsw(:,:)
         311             :    real(r8),          optional, pointer     :: refitabsw(:,:)
         312             :    real(r8),          optional, pointer     :: refrtablw(:,:)
         313             :    real(r8),          optional, pointer     :: refitablw(:,:)
         314             :    character(len=20), optional, intent(out) :: aername
         315             :    real(r8),          optional, intent(out) :: density_aer
         316             :    real(r8),          optional, intent(out) :: hygro_aer
         317             :    real(r8),          optional, intent(out) :: dryrad_aer
         318             :    real(r8),          optional, intent(out) :: dispersion_aer
         319             :    real(r8),          optional, intent(out) :: num_to_mass_aer
         320             :    integer,           optional, intent(out) :: ncoef
         321             :    integer,           optional, intent(out) :: prefr
         322             :    integer,           optional, intent(out) :: prefi
         323             :    real(r8),          optional, intent(out) :: sigmag
         324             :    real(r8),          optional, intent(out) :: dgnum
         325             :    real(r8),          optional, intent(out) :: dgnumlo
         326             :    real(r8),          optional, intent(out) :: dgnumhi
         327             :    real(r8),          optional, intent(out) :: rhcrystal
         328             :    real(r8),          optional, intent(out) :: rhdeliques
         329             :
         330             :    ! Local variables
         331             :    character(len=*), parameter :: subname = 'physprop_get'
         332             :    !------------------------------------------------------------------------------------
         333             :
         334           0 :    if (id <= 0 .or. id > numphysprops) then
         335           0 :       write(iulog,*) subname//': illegal ID value: ', id
         336           0 :       call endrun('physprop_get: ID out of range')
         337             :    end if
         338             :
         339           0 :    if (present(sourcefile))        sourcefile        =  physprop(id)%sourcefile
         340           0 :    if (present(opticstype))        opticstype        =  physprop(id)%opticsmethod
         341           0 :    if (present(sw_hygro_ext))      sw_hygro_ext      => physprop(id)%sw_hygro_ext
         342           0 :    if (present(sw_hygro_ssa))      sw_hygro_ssa      => physprop(id)%sw_hygro_ssa
         343           0 :    if (present(sw_hygro_asm))      sw_hygro_asm      => physprop(id)%sw_hygro_asm
         344           0 :    if (present(lw_hygro_abs))      lw_hygro_abs      => physprop(id)%lw_hygro_abs
         345           0 :    if (present(sw_nonhygro_ext))   sw_nonhygro_ext   => physprop(id)%sw_nonhygro_ext
         346           0 :    if (present(sw_nonhygro_ssa))   sw_nonhygro_ssa   => physprop(id)%sw_nonhygro_ssa
         347           0 :    if (present(sw_nonhygro_asm))   sw_nonhygro_asm   => physprop(id)%sw_nonhygro_asm
         348           0 :    if (present(sw_nonhygro_scat))  sw_nonhygro_scat  => physprop(id)%sw_nonhygro_scat
         349           0 :    if (present(sw_nonhygro_ascat)) sw_nonhygro_ascat => physprop(id)%sw_nonhygro_ascat
         350           0 :    if (present(lw_abs))            lw_abs            => physprop(id)%lw_abs
         351             :
         352           0 :    if (present(refindex_aer_sw))   refindex_aer_sw   => physprop(id)%refindex_aer_sw
         353           0 :    if (present(refindex_aer_lw))   refindex_aer_lw   => physprop(id)%refindex_aer_lw
         354             :
         355           0 :    if (present(r_sw_ext))          r_sw_ext      => physprop(id)%r_sw_ext
         356           0 :    if (present(r_sw_scat))         r_sw_scat     => physprop(id)%r_sw_scat
         357           0 :    if (present(r_sw_ascat))        r_sw_ascat    => physprop(id)%r_sw_ascat
         358           0 :    if (present(r_lw_abs))          r_lw_abs      => physprop(id)%r_lw_abs
         359           0 :    if (present(mu))                mu            => physprop(id)%mu
         360             :
         361           0 :    if (present(extpsw))            extpsw        => physprop(id)%extpsw
         362           0 :    if (present(abspsw))            abspsw        => physprop(id)%abspsw
         363           0 :    if (present(asmpsw))            asmpsw        => physprop(id)%asmpsw
         364           0 :    if (present(absplw))            absplw        => physprop(id)%absplw
         365           0 :    if (present(refrtabsw))         refrtabsw     => physprop(id)%refrtabsw
         366           0 :    if (present(refitabsw))         refitabsw     => physprop(id)%refitabsw
         367           0 :    if (present(refrtablw))         refrtablw     => physprop(id)%refrtablw
         368           0 :    if (present(refitablw))         refitablw     => physprop(id)%refitablw
         369             :
         370           0 :    if (present(aername))         aername         =  physprop(id)%aername
         371           0 :    if (present(density_aer))     density_aer     =  physprop(id)%density_aer
         372           0 :    if (present(hygro_aer))       hygro_aer       =  physprop(id)%hygro_aer
         373           0 :    if (present(dryrad_aer))      dryrad_aer      =  physprop(id)%dryrad_aer
         374           0 :    if (present(dispersion_aer))  dispersion_aer  =  physprop(id)%dispersion_aer
         375           0 :    if (present(num_to_mass_aer)) num_to_mass_aer =  physprop(id)%num_to_mass_aer
         376             :
         377           0 :    if (present(ncoef))           ncoef           =  physprop(id)%ncoef
         378           0 :    if (present(prefr))           prefr           =  physprop(id)%prefr
         379           0 :    if (present(prefi))           prefi           =  physprop(id)%prefi
         380           0 :    if (present(sigmag))          sigmag          =  physprop(id)%sigmag
         381           0 :    if (present(dgnum))           dgnum           =  physprop(id)%dgnum
         382           0 :    if (present(dgnumlo))         dgnumlo         =  physprop(id)%dgnumlo
         383           0 :    if (present(dgnumhi))         dgnumhi         =  physprop(id)%dgnumhi
         384           0 :    if (present(rhcrystal))       rhcrystal       =  physprop(id)%rhcrystal
         385           0 :    if (present(rhdeliques))      rhdeliques      =  physprop(id)%rhdeliques
         386             :
         387           0 : end subroutine physprop_get
         388             :
         389             : !================================================================================================
         390             : ! Private methods
         391             : !================================================================================================
         392             :
         393           0 : subroutine aerosol_optics_init(phys_prop, nc_id)
         394             :
         395             :    ! Determine the opticstype, then call the
         396             :    ! appropriate routine to read the data.
         397             :
         398             :    type(physprop_type), intent(inout) :: phys_prop  ! data after interp onto cam rh mesh
         399             :    type(file_desc_t),   intent(inout) :: nc_id      ! indentifier for netcdf file
         400             :
         401             :    integer :: opticslength_id, opticslength
         402             :    type(var_desc_t) :: op_type_id
         403             :    integer :: ierr ! mpi error codes
         404             :    character(len=ot_length)  :: opticstype_str ! string read from netCDF file -- may contain trailing
         405             :                                         ! nulls which aren't dealt with by trim()
         406             :    !------------------------------------------------------------------------------------
         407             :
         408           0 :    ierr = pio_inq_dimid(nc_id, 'opticsmethod_len', opticslength_id)
         409           0 :    ierr = pio_inq_dimlen(nc_id, opticslength_id, opticslength)
         410           0 :    if ( opticslength .gt. ot_length ) then
         411           0 :       call endrun(" optics type length in "//phys_prop%sourcefile//" excedes maximum length of 32")
         412             :    endif
         413           0 :    ierr = pio_inq_varid(nc_id, 'opticsmethod', op_type_id)
         414           0 :    ierr = pio_get_var(nc_id, op_type_id,phys_prop%opticsmethod )
         415             :
         416           0 :    select case (phys_prop%opticsmethod)
         417             :    case ('zero')
         418           0 :       call zero_optics_init(phys_prop, nc_id)
         419             :
         420             :    case ('hygro')
         421           0 :       call hygro_optics_init(phys_prop, nc_id)
         422             :
         423             :    case ('hygroscopic')
         424           0 :       call hygroscopic_optics_init(phys_prop, nc_id)
         425             :
         426             :    case ('nonhygro')
         427           0 :       call nonhygro_optics_init(phys_prop, nc_id)
         428             :
         429             :    case ('insoluble')
         430           0 :       call insoluble_optics_init(phys_prop, nc_id)
         431             :
         432             :    case ('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3')
         433           0 :       call volcanic_radius_optics_init(phys_prop, nc_id)
         434             :
         435             :    case ('volcanic')
         436           0 :       call volcanic_optics_init(phys_prop, nc_id)
         437             :
         438             :    case ('modal')
         439           0 :       call modal_optics_init(phys_prop, nc_id)
         440             :
         441             :    ! other types of optics can be added here
         442             :
         443             :    case default
         444             :       call endrun('aerosol_optics_init: unsupported optics type '//&
         445           0 :          trim(phys_prop%opticsmethod)//' in file '//phys_prop%sourcefile)
         446             :    end select
         447             :
         448           0 : end subroutine aerosol_optics_init
         449             :
         450             : !================================================================================================
         451             :
         452           0 : subroutine hygro_optics_init(phys_prop, nc_id)
         453             :
         454             :    ! Read optics data of type 'hygro' and interpolate it to CAM's rh mesh.
         455             :
         456             :    type (physprop_type), intent(inout) :: phys_prop  ! data after interp onto cam rh mesh
         457             :    type (file_desc_t),   intent(inout) :: nc_id      ! indentifier for netcdf file
         458             :
         459             :    ! Local variables
         460             :    integer :: ierr ! error flag
         461             :
         462             :    integer :: rh_idx_id, lw_band_id, sw_band_id
         463             :    integer :: kbnd, krh
         464             :    integer :: rh_id, sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
         465             :    integer :: nbnd, swbands
         466             :
         467             :    ! temp data from hygroscopic file before interpolation onto cam-rh-mesh
         468             :    integer  :: nfilerh ! number of rh values in file
         469           0 :    real(r8), allocatable, dimension(:) :: frh
         470           0 :    real(r8), allocatable, dimension(:,:)  :: fsw_ext
         471           0 :    real(r8), allocatable, dimension(:,:)  :: fsw_ssa
         472           0 :    real(r8), allocatable, dimension(:,:)  :: fsw_asm
         473             :
         474             :    real(r8) :: rh ! real rh value on cam rh mesh (indexvalue)
         475             :    !------------------------------------------------------------------------------------
         476             :
         477           0 :    allocate(phys_prop%sw_hygro_ext(nrh,nswbands))
         478           0 :    allocate(phys_prop%sw_hygro_ssa(nrh,nswbands))
         479           0 :    allocate(phys_prop%sw_hygro_asm(nrh,nswbands))
         480           0 :    allocate(phys_prop%lw_abs(nlwbands))
         481             :
         482           0 :    ierr = pio_inq_dimid(nc_id, 'rh_idx', rh_idx_id)
         483             :
         484           0 :    ierr = pio_inq_dimlen(nc_id, rh_idx_id, nfilerh)
         485             :
         486           0 :    ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
         487             :
         488           0 :    ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
         489             :
         490           0 :    ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
         491             :
         492           0 :    if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
         493           0 :         ' has the wrong number of lwbands')
         494             :
         495           0 :    ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
         496             :
         497           0 :    if(swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
         498           0 :          ' has the wrong number of sw bands')
         499             :
         500           0 :    ierr = pio_inq_varid(nc_id, 'rh', rh_id)
         501             :
         502           0 :    ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
         503             :
         504           0 :    ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
         505             :
         506           0 :    ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
         507             :
         508           0 :    ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
         509             :
         510             :    ! specific optical properties on file's rh mesh
         511           0 :    allocate(fsw_ext(nfilerh,nswbands))
         512           0 :    allocate(fsw_asm(nfilerh,nswbands))
         513           0 :    allocate(fsw_ssa(nfilerh,nswbands))
         514           0 :    allocate(frh(nfilerh))
         515             :
         516           0 :    ierr = pio_get_var(nc_id, rh_id, frh)
         517             :
         518           0 :    ierr = pio_get_var(nc_id, sw_ext_id, fsw_ext)
         519             :
         520           0 :    ierr = pio_get_var(nc_id, sw_ssa_id, fsw_ssa)
         521             :
         522           0 :    ierr = pio_get_var(nc_id, sw_asm_id, fsw_asm)
         523             :
         524           0 :    ierr = pio_get_var(nc_id, lw_ext_id, phys_prop%lw_abs)
         525             :
         526             :    ! interpolate onto cam's rh mesh
         527           0 :    do kbnd = 1,nswbands
         528           0 :       do krh = 1, nrh
         529           0 :          rh = 1.0_r8 / nrh * (krh - 1)
         530           0 :          phys_prop%sw_hygro_ext(krh,kbnd) = &
         531           0 :             exp_interpol( frh, fsw_ext(:,kbnd) / fsw_ext(1,kbnd), rh ) &
         532           0 :             * fsw_ext(1, kbnd)
         533           0 :          phys_prop%sw_hygro_ssa(krh,kbnd) = &
         534           0 :             lin_interpol( frh, fsw_ssa(:,kbnd) / fsw_ssa(1,kbnd), rh ) &
         535           0 :             * fsw_ssa(1, kbnd)
         536           0 :          phys_prop%sw_hygro_asm(krh,kbnd) = &
         537           0 :             lin_interpol( frh, fsw_asm(:,kbnd) / fsw_asm(1,kbnd), rh ) &
         538           0 :             * fsw_asm(1, kbnd)
         539             :       enddo
         540             :    enddo
         541             :
         542           0 :    deallocate (fsw_ext, fsw_asm, fsw_ssa, frh)
         543             :
         544             :    ! read refractive index data if available
         545           0 :    call refindex_aer_init(phys_prop, nc_id)
         546             :
         547             :    ! read bulk aero props
         548           0 :    call bulk_props_init(phys_prop, nc_id)
         549             :
         550           0 : end subroutine hygro_optics_init
         551             :
         552             : !================================================================================================
         553             :
         554           0 : subroutine zero_optics_init(phys_prop, nc_id)
         555             :
         556             :    ! Read optics data of type 'nonhygro'
         557             :
         558             :    type (physprop_type), intent(inout) :: phys_prop  ! storage for file data
         559             :    type (file_desc_t),   intent(inout) :: nc_id      ! indentifier for netcdf file
         560             :
         561             :    ! Local variables
         562             :    integer :: lw_band_id, sw_band_id
         563             :    integer :: sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
         564             :    integer :: swbands, nbnd
         565             :    integer :: ierr ! error flag
         566             :    !------------------------------------------------------------------------------------
         567             :
         568             :    ! perhaps this doesn't even need allocated.
         569           0 :    allocate (phys_prop%sw_nonhygro_ext(nswbands))
         570           0 :    allocate (phys_prop%sw_nonhygro_ssa(nswbands))
         571           0 :    allocate (phys_prop%sw_nonhygro_asm(nswbands))
         572           0 :    allocate (phys_prop%lw_abs(nlwbands))
         573             :
         574           0 :    phys_prop%sw_nonhygro_ext = 0._r8
         575           0 :    phys_prop%sw_nonhygro_ssa = 0._r8
         576           0 :    phys_prop%sw_nonhygro_asm = 0._r8
         577           0 :    phys_prop%lw_abs = 0._r8
         578             :
         579           0 : end subroutine zero_optics_init
         580             :
         581             : !================================================================================================
         582             :
         583           0 : subroutine insoluble_optics_init(phys_prop, nc_id)
         584             :
         585             :    ! Read optics data of type 'nonhygro'
         586             :
         587             :    type (physprop_type), intent(inout) :: phys_prop  ! storage for file data
         588             :    type (file_desc_t),   intent(inout) :: nc_id      ! indentifier for netcdf file
         589             :
         590             :    ! Local variables
         591             :    integer :: lw_band_id, sw_band_id
         592             :    integer :: sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
         593             :    integer :: swbands, nbnd
         594             :    integer :: ierr ! error flag
         595             :    integer :: start(2), count(2)
         596             :    !------------------------------------------------------------------------------------
         597             :
         598           0 :    allocate (phys_prop%sw_nonhygro_ext(nswbands))
         599           0 :    allocate (phys_prop%sw_nonhygro_ssa(nswbands))
         600           0 :    allocate (phys_prop%sw_nonhygro_asm(nswbands))
         601           0 :    allocate (phys_prop%lw_abs(nlwbands))
         602             :
         603           0 :    ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
         604             :
         605           0 :    ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
         606             :
         607           0 :    ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
         608             :
         609           0 :    if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
         610           0 :         ' has the wrong number of lwbands')
         611             :
         612           0 :    ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
         613             :
         614           0 :    if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
         615           0 :         ' has the wrong number of sw bands')
         616             :
         617             :    ! read file data
         618           0 :    ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
         619           0 :    ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
         620           0 :    ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
         621           0 :    ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
         622             :
         623           0 :    start = 1
         624           0 :    count=(/1,swbands/)
         625             :
         626           0 :    ierr = pio_get_var(nc_id, sw_ext_id, start, count, phys_prop%sw_nonhygro_ext)
         627           0 :    ierr = pio_get_var(nc_id, sw_ssa_id, start, count, phys_prop%sw_nonhygro_ssa)
         628           0 :    ierr = pio_get_var(nc_id, sw_asm_id, start, count, phys_prop%sw_nonhygro_asm)
         629           0 :    count = (/1,nbnd/)
         630           0 :    ierr = pio_get_var(nc_id, lw_ext_id, start, count, phys_prop%lw_abs)
         631             :
         632             :    ! read refractive index data if available
         633           0 :    call refindex_aer_init(phys_prop, nc_id)
         634             :
         635             :    ! read bulk aero props
         636           0 :    call bulk_props_init(phys_prop, nc_id)
         637             :
         638           0 : end subroutine insoluble_optics_init
         639             :
         640             : !================================================================================================
         641             :
         642           0 : subroutine volcanic_radius_optics_init(phys_prop, nc_id)
         643             :
         644             :    ! Read optics data of type 'volcanic_radius'
         645             :
         646             :    type (physprop_type), intent(inout) :: phys_prop  ! storage for file data
         647             :    type (file_desc_t),   intent(inout) :: nc_id      ! indentifier for netcdf file
         648             :
         649             :    ! Local variables
         650             :    integer :: lw_band_id, sw_band_id, mu_id, mu_did
         651             :    integer :: sw_ext_id, sw_scat_id, sw_ascat_id, lw_abs_id
         652             :    integer :: swbands, nbnd, n_mu_samples
         653             :    integer :: ierr ! error flag
         654             :    !------------------------------------------------------------------------------------
         655             :
         656           0 :    ierr = pio_inq_dimid(nc_id, 'mu_samples', mu_did)
         657           0 :    ierr = pio_inq_dimlen(nc_id, mu_did, n_mu_samples)
         658             :
         659           0 :    allocate (phys_prop%r_sw_ext(nswbands,n_mu_samples))
         660           0 :    allocate (phys_prop%r_sw_scat(nswbands,n_mu_samples))
         661           0 :    allocate (phys_prop%r_sw_ascat(nswbands,n_mu_samples))
         662           0 :    allocate (phys_prop%r_lw_abs(nlwbands,n_mu_samples))
         663           0 :    allocate (phys_prop%mu(n_mu_samples))
         664             :
         665           0 :    ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
         666             :
         667           0 :    ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
         668             :
         669           0 :    ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
         670             :
         671           0 :    if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
         672           0 :         ' has the wrong number of lwbands')
         673             :
         674           0 :    ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
         675             :
         676           0 :    if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
         677           0 :         ' has the wrong number of sw bands')
         678             :
         679             :    ! read file data
         680           0 :    ierr = pio_inq_varid(nc_id, 'bext_sw', sw_ext_id)
         681           0 :    ierr = pio_inq_varid(nc_id, 'bsca_sw', sw_scat_id)
         682           0 :    ierr = pio_inq_varid(nc_id, 'basc_sw', sw_ascat_id)
         683           0 :    ierr = pio_inq_varid(nc_id, 'babs_lw', lw_abs_id)
         684           0 :    ierr = pio_inq_varid(nc_id, 'mu_samples', mu_id)
         685             :
         686           0 :    ierr = pio_get_var(nc_id, sw_ext_id, phys_prop%r_sw_ext)
         687           0 :    ierr = pio_get_var(nc_id, sw_scat_id, phys_prop%r_sw_scat)
         688           0 :    ierr = pio_get_var(nc_id, sw_ascat_id, phys_prop%r_sw_ascat)
         689           0 :    ierr = pio_get_var(nc_id, lw_abs_id, phys_prop%r_lw_abs)
         690           0 :    ierr = pio_get_var(nc_id, mu_id, phys_prop%mu)
         691             :
         692             :    ! read bulk aero props
         693           0 :    call bulk_props_init(phys_prop, nc_id)
         694             :
         695           0 : end subroutine volcanic_radius_optics_init
         696             :
         697             : !================================================================================================
         698             :
         699           0 : subroutine volcanic_optics_init(phys_prop, nc_id)
         700             :
         701             :    ! Read optics data of type 'volcanic'
         702             :
         703             :    type (physprop_type), intent(inout) :: phys_prop  ! storage for file data
         704             :    type (file_desc_t)  , intent(inout) :: nc_id      ! indentifier for netcdf file
         705             :
         706             :    ! Local variables
         707             :    integer :: lw_band_id, sw_band_id
         708             :    integer :: sw_ext_id, sw_scat_id, sw_ascat_id, lw_abs_id
         709             :    integer :: swbands, nbnd
         710             :    integer :: ierr ! error flag
         711             :    !------------------------------------------------------------------------------------
         712             :
         713           0 :    allocate (phys_prop%sw_nonhygro_ext(nswbands))
         714           0 :    allocate (phys_prop%sw_nonhygro_scat(nswbands))
         715           0 :    allocate (phys_prop%sw_nonhygro_ascat(nswbands))
         716           0 :    allocate (phys_prop%lw_abs(nlwbands))
         717             :
         718           0 :    ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
         719           0 :    ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
         720             :
         721           0 :    ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
         722             :
         723           0 :    if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
         724           0 :         ' has the wrong number of lwbands')
         725             :
         726           0 :    ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
         727           0 :    if(masterproc) write(iulog,*) 'swbands',swbands
         728             :
         729           0 :    if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
         730           0 :         ' has the wrong number of sw bands')
         731             :
         732             :    ! read file data
         733           0 :    ierr = pio_inq_varid(nc_id, 'bext_sw', sw_ext_id)
         734           0 :    ierr = pio_inq_varid(nc_id, 'bsca_sw', sw_scat_id)
         735           0 :    ierr = pio_inq_varid(nc_id, 'basc_sw', sw_ascat_id)
         736           0 :    ierr = pio_inq_varid(nc_id, 'babs_lw', lw_abs_id)
         737             :
         738           0 :    ierr = pio_get_var(nc_id, sw_ext_id, phys_prop%sw_nonhygro_ext)
         739           0 :    ierr = pio_get_var(nc_id, sw_scat_id, phys_prop%sw_nonhygro_scat)
         740           0 :    ierr = pio_get_var(nc_id, sw_ascat_id, phys_prop%sw_nonhygro_ascat)
         741           0 :    ierr = pio_get_var(nc_id, lw_abs_id, phys_prop%lw_abs)
         742             :
         743             :    ! read bulk aero props
         744           0 :    call bulk_props_init(phys_prop, nc_id)
         745             :
         746           0 : end subroutine volcanic_optics_init
         747             :
         748             : !================================================================================================
         749             :
         750           0 : subroutine hygroscopic_optics_init(phys_prop, nc_id)
         751             :
         752             :    ! Read optics data of type 'hygroscopic' and interpolate it to CAM's rh mesh.
         753             :
         754             :    type (physprop_type), intent(inout) :: phys_prop  ! data after interp onto cam rh mesh
         755             :    type (file_desc_T),   intent(inout) :: nc_id      ! indentifier for netcdf file
         756             :
         757             :    ! Local variables
         758             :    integer :: ierr ! error flag
         759             :
         760             :    integer :: rh_idx_id, lw_band_id, sw_band_id
         761             :    integer :: kbnd, krh
         762             :    integer :: rh_id, sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
         763             :    integer :: nbnd, swbands
         764             :
         765             :    ! temp data from hygroscopic file before interpolation onto cam-rh-mesh
         766             :    integer  :: nfilerh ! number of rh values in file
         767           0 :    real(r8), allocatable, dimension(:) :: frh
         768           0 :    real(r8), allocatable, dimension(:,:)  :: fsw_ext
         769           0 :    real(r8), allocatable, dimension(:,:)  :: fsw_ssa
         770           0 :    real(r8), allocatable, dimension(:,:)  :: fsw_asm
         771           0 :    real(r8), allocatable, dimension(:,:)  :: flw_abs
         772             :
         773             :    real(r8) :: rh ! real rh value on cam rh mesh (indexvalue)
         774             :    character(len=*), parameter :: sub = 'hygroscopic_optics_init'
         775             :    !------------------------------------------------------------------------------------
         776             :
         777           0 :    allocate(phys_prop%sw_hygro_ext(nrh,nswbands))
         778           0 :    allocate(phys_prop%sw_hygro_ssa(nrh,nswbands))
         779           0 :    allocate(phys_prop%sw_hygro_asm(nrh,nswbands))
         780           0 :    allocate(phys_prop%lw_hygro_abs(nrh,nlwbands))
         781             :
         782           0 :    ierr = pio_inq_dimid(nc_id, 'rh_idx', rh_idx_id)
         783           0 :    ierr = pio_inq_dimlen(nc_id, rh_idx_id, nfilerh)
         784             :
         785           0 :    ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
         786           0 :    ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
         787           0 :    if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
         788           0 :         ' has the wrong number of lwbands')
         789             :
         790           0 :    ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
         791           0 :    ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
         792           0 :    if(swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
         793           0 :         ' has the wrong number of sw bands')
         794             :
         795           0 :    ierr = pio_inq_varid(nc_id, 'rh', rh_id)
         796           0 :    ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
         797           0 :    ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
         798           0 :    ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
         799           0 :    ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
         800             :
         801             :    ! specific optical properties on file's rh mesh
         802           0 :    allocate(fsw_ext(nfilerh,nswbands))
         803           0 :    allocate(fsw_asm(nfilerh,nswbands))
         804           0 :    allocate(fsw_ssa(nfilerh,nswbands))
         805           0 :    allocate(flw_abs(nfilerh,nlwbands))
         806           0 :    allocate(frh(nfilerh))
         807             :
         808           0 :    ierr = pio_get_var(nc_id, rh_id, frh)
         809           0 :    ierr = pio_get_var(nc_id, sw_ext_id, fsw_ext)
         810           0 :    ierr = pio_get_var(nc_id, sw_ssa_id, fsw_ssa)
         811           0 :    ierr = pio_get_var(nc_id, sw_asm_id, fsw_asm)
         812           0 :    ierr = pio_get_var(nc_id, lw_ext_id, flw_abs)
         813             :
         814             :    ! interpolate onto cam's rh mesh
         815           0 :    do kbnd = 1,nswbands
         816           0 :       do krh = 1, nrh
         817           0 :          rh = 1.0_r8 / nrh * (krh - 1)
         818           0 :          phys_prop%sw_hygro_ext(krh,kbnd) = &
         819           0 :             exp_interpol( frh, fsw_ext(:,kbnd) / fsw_ext(1,kbnd), rh ) &
         820           0 :             * fsw_ext(1, kbnd)
         821           0 :          phys_prop%sw_hygro_ssa(krh,kbnd) = &
         822           0 :             lin_interpol( frh, fsw_ssa(:,kbnd) / fsw_ssa(1,kbnd), rh ) &
         823           0 :             * fsw_ssa(1, kbnd)
         824           0 :          phys_prop%sw_hygro_asm(krh,kbnd) = &
         825           0 :             lin_interpol( frh, fsw_asm(:,kbnd) / fsw_asm(1,kbnd), rh ) &
         826           0 :             * fsw_asm(1, kbnd)
         827             :       enddo
         828             :    enddo
         829           0 :    do kbnd = 1,nlwbands
         830           0 :       do krh = 1, nrh
         831           0 :          rh = 1.0_r8 / nrh * (krh - 1)
         832           0 :          phys_prop%lw_hygro_abs(krh,kbnd) = &
         833           0 :             exp_interpol( frh, flw_abs(:,kbnd) / flw_abs(1,kbnd), rh ) &
         834           0 :             * flw_abs(1, kbnd)
         835             :       enddo
         836             :    enddo
         837             :
         838           0 :    deallocate (fsw_ext, fsw_asm, fsw_ssa, flw_abs, frh)
         839             :
         840             :    ! read refractive index data if available
         841           0 :    call refindex_aer_init(phys_prop, nc_id)
         842             :
         843             :    ! read bulk aero props
         844           0 :    call bulk_props_init(phys_prop, nc_id)
         845             :
         846           0 : end subroutine hygroscopic_optics_init
         847             :
         848             : !================================================================================================
         849             :
         850           0 : subroutine nonhygro_optics_init(phys_prop, nc_id)
         851             :
         852             :    ! Read optics data of type 'nonhygro'
         853             :
         854             :    type (physprop_type), intent(inout) :: phys_prop  ! storage for file data
         855             :    type (file_desc_t)  , intent(inout) :: nc_id      ! indentifier for netcdf file
         856             :
         857             :    ! Local variables
         858             :    integer :: lw_band_id, sw_band_id
         859             :    integer :: sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
         860             :    integer :: swbands, nbnd
         861             :    integer :: ierr ! error flag
         862             :    !------------------------------------------------------------------------------------
         863             :
         864           0 :    allocate (phys_prop%sw_nonhygro_ext(nswbands))
         865           0 :    allocate (phys_prop%sw_nonhygro_ssa(nswbands))
         866           0 :    allocate (phys_prop%sw_nonhygro_asm(nswbands))
         867           0 :    allocate (phys_prop%lw_abs(nlwbands))
         868             :
         869           0 :    ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
         870           0 :    ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
         871             :
         872           0 :    ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
         873             :
         874           0 :    if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
         875           0 :         ' has the wrong number of lwbands')
         876             :
         877           0 :    ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
         878             :
         879           0 :    if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
         880           0 :         ' has the wrong number of sw bands')
         881             :
         882             :    ! read file data
         883           0 :    ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
         884           0 :    ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
         885           0 :    ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
         886           0 :    ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
         887             :
         888           0 :    ierr = pio_get_var(nc_id, sw_ext_id, phys_prop%sw_nonhygro_ext)
         889           0 :    ierr = pio_get_var(nc_id, sw_ssa_id, phys_prop%sw_nonhygro_ssa)
         890           0 :    ierr = pio_get_var(nc_id, sw_asm_id, phys_prop%sw_nonhygro_asm)
         891           0 :    ierr = pio_get_var(nc_id, lw_ext_id, phys_prop%lw_abs)
         892             :
         893             :    ! read refractive index data if available
         894           0 :    call refindex_aer_init(phys_prop, nc_id)
         895             :
         896             :    ! read bulk aero props
         897           0 :    call bulk_props_init(phys_prop, nc_id)
         898             :
         899           0 : end subroutine nonhygro_optics_init
         900             :
         901             : !================================================================================================
         902             :
         903           0 : subroutine refindex_aer_init(phys_prop, nc_id)
         904             :
         905             : !  Read refractive indices of aerosol
         906             :
         907             :    type (physprop_type), intent(inout) :: phys_prop  ! storage for file data
         908             :    type (file_desc_T),   intent(inout) :: nc_id      ! indentifier for netcdf file
         909             :
         910             :    ! Local variables
         911             :    integer :: i
         912             :    integer :: istat1, istat2, istat3     ! status flags
         913             :    integer :: vid_real, vid_im           ! variable ids
         914           0 :    real(r8), pointer :: ref_real(:), ref_im(:)  ! tmp storage for components of complex index
         915             :    character(len=*), parameter :: subname = 'refindex_aer_init'
         916             :    !------------------------------------------------------------------------------------
         917             :
         918             :    ! assume that the dimensions lw_band and sw_band have already been checked
         919             :    ! by the calling subroutine
         920             :
         921             :    ! Check that the variables are present before allocating storage and reading.
         922             :    ! Since we're setting complex data values, both the real and imaginary parts must
         923             :    ! be present or neither will be read.
         924             :
         925             :    ! set PIO to return control to the caller when variable not found
         926           0 :    call pio_seterrorhandling(nc_id, PIO_BCAST_ERROR)
         927             :
         928           0 :    istat1 = pio_inq_varid(nc_id, 'refindex_real_aer_sw', vid_real)
         929           0 :    istat2 = pio_inq_varid(nc_id, 'refindex_im_aer_sw',   vid_im)
         930             :
         931           0 :    if (istat1 == PIO_NOERR  .and. istat2 == PIO_NOERR) then
         932             :
         933           0 :       allocate(ref_real(nswbands), ref_im(nswbands))
         934             :
         935           0 :       istat3 = pio_get_var(nc_id, vid_real, ref_real)
         936           0 :       if (istat3 /= PIO_NOERR) then
         937           0 :          call endrun(subname//': ERROR reading refindex_real_aer_sw')
         938             :       end if
         939             :
         940           0 :       istat3 = pio_get_var(nc_id, vid_im, ref_im)
         941           0 :       if (istat3 /= PIO_NOERR) then
         942           0 :          call endrun(subname//': ERROR reading refindex_im_aer_sw')
         943             :       end if
         944             :
         945             :       ! successfully read refindex data -- set complex values in physprop object
         946           0 :       allocate(phys_prop%refindex_aer_sw(nswbands))
         947           0 :       do i = 1, nswbands
         948           0 :          phys_prop%refindex_aer_sw(i) = cmplx(ref_real(i), abs(ref_im(i)),&
         949           0 :               kind=r8)
         950             :       end do
         951             :
         952           0 :       deallocate(ref_real, ref_im)
         953             :
         954             :    end if
         955             :
         956           0 :    istat1 = pio_inq_varid(nc_id, 'refindex_real_aer_lw', vid_real)
         957           0 :    istat2 = pio_inq_varid(nc_id, 'refindex_im_aer_lw',   vid_im)
         958             :
         959           0 :    if (istat1 == PIO_NOERR  .and. istat2 == PIO_NOERR) then
         960             :
         961           0 :       allocate(ref_real(nlwbands), ref_im(nlwbands))
         962             :
         963           0 :       istat3 = pio_get_var(nc_id, vid_real, ref_real)
         964           0 :       if (istat3 /= PIO_NOERR) then
         965           0 :          call endrun(subname//': ERROR reading refindex_real_aer_lw')
         966             :       end if
         967             :
         968           0 :       istat3 = pio_get_var(nc_id, vid_im, ref_im)
         969           0 :       if (istat3 /= PIO_NOERR) then
         970           0 :          call endrun(subname//': ERROR reading refindex_im_aer_lw')
         971             :       end if
         972             :
         973             :       ! successfully read refindex data -- set complex value in physprop object
         974           0 :       allocate(phys_prop%refindex_aer_lw(nlwbands))
         975           0 :       do i = 1, nlwbands
         976           0 :          phys_prop%refindex_aer_lw(i) = cmplx(ref_real(i), abs(ref_im(i)),&
         977           0 :               kind=r8)
         978             :       end do
         979             :
         980           0 :       deallocate(ref_real, ref_im)
         981             :
         982             :    end if
         983             :
         984             :    ! reset PIO to handle errors internally
         985           0 :    call pio_seterrorhandling(nc_id, PIO_INTERNAL_ERROR)
         986             :
         987           0 : end subroutine refindex_aer_init
         988             :
         989             : !================================================================================================
         990             :
         991           0 : subroutine modal_optics_init(props, ncid)
         992             :
         993             : !  Read optics data for modal aerosols
         994             :
         995             :    type (physprop_type), intent(inout) :: props   ! storage for file data
         996             :    type (file_desc_T),   intent(inout) :: ncid    ! indentifier for netcdf file
         997             :
         998             :    ! Local variables
         999             :    integer :: ierr
        1000             :    integer :: did
        1001             :    integer :: ival
        1002             :    type(var_desc_t) :: vid
        1003           0 :    real(r8), pointer :: rval(:,:,:,:,:) ! temp array used to eliminate a singleton dimension
        1004             :
        1005             :    character(len=*), parameter :: subname = 'modal_optics_init'
        1006             :    !------------------------------------------------------------------------------------
        1007             :
        1008             :    ! Check dimensions for number of lw and sw bands
        1009             :
        1010           0 :    ierr = pio_inq_dimid(ncid, 'lw_band', did)
        1011           0 :    ierr = pio_inq_dimlen(ncid, did, ival)
        1012           0 :    if (ival .ne. nlwbands) call endrun(subname//':'//props%sourcefile// &
        1013           0 :         ' has the wrong number of lw bands')
        1014             :
        1015           0 :    ierr = pio_inq_dimid(ncid, 'sw_band', did)
        1016           0 :    ierr = pio_inq_dimlen(ncid, did, ival)
        1017           0 :    if (ival .ne. nswbands) call endrun(subname//':'//props%sourcefile// &
        1018           0 :         ' has the wrong number of sw bands')
        1019             :
        1020             :    ! Get other dimensions
        1021           0 :    ierr = pio_inq_dimid(ncid, 'coef_number', did)
        1022           0 :    ierr = pio_inq_dimlen(ncid, did, props%ncoef)
        1023             :
        1024           0 :    ierr = pio_inq_dimid(ncid, 'refindex_real', did)
        1025           0 :    ierr = pio_inq_dimlen(ncid, did, props%prefr)
        1026             :
        1027           0 :    ierr = pio_inq_dimid(ncid, 'refindex_im', did)
        1028           0 :    ierr = pio_inq_dimlen(ncid, did, props%prefi)
        1029             :
        1030             :    ! Allocate arrays
        1031             :    allocate( &
        1032             :       props%extpsw(props%ncoef,props%prefr,props%prefi,nswbands), &
        1033             :       props%abspsw(props%ncoef,props%prefr,props%prefi,nswbands), &
        1034             :       props%asmpsw(props%ncoef,props%prefr,props%prefi,nswbands), &
        1035             :       props%absplw(props%ncoef,props%prefr,props%prefi,nlwbands), &
        1036             :       props%refrtabsw(props%prefr,nswbands), &
        1037             :       props%refitabsw(props%prefi,nswbands), &
        1038             :       props%refrtablw(props%prefr,nlwbands), &
        1039           0 :       props%refitablw(props%prefi,nlwbands)  )
        1040             :
        1041             :
        1042             :    ! allocate temp to remove the mode dimension from the sw variables
        1043           0 :    allocate(rval(props%ncoef,props%prefr,props%prefi,1,nswbands))
        1044             :
        1045           0 :    ierr = pio_inq_varid(ncid, 'extpsw', vid)
        1046           0 :    ierr = pio_get_var(ncid, vid, rval)
        1047           0 :    props%extpsw = rval(:,:,:,1,:)
        1048             :
        1049           0 :    ierr = pio_inq_varid(ncid, 'abspsw', vid)
        1050           0 :    ierr = pio_get_var(ncid, vid, rval)
        1051           0 :    props%abspsw = rval(:,:,:,1,:)
        1052             :
        1053           0 :    ierr = pio_inq_varid(ncid, 'asmpsw', vid)
        1054           0 :    ierr = pio_get_var(ncid, vid, rval)
        1055           0 :    props%asmpsw = rval(:,:,:,1,:)
        1056             :
        1057           0 :    deallocate(rval)
        1058             :
        1059             :    ! allocate temp to remove the mode dimension from the lw variables
        1060           0 :    allocate(rval(props%ncoef,props%prefr,props%prefi,1,nlwbands))
        1061             :
        1062           0 :    ierr = pio_inq_varid(ncid, 'absplw', vid)
        1063           0 :    ierr = pio_get_var(ncid, vid, rval)
        1064           0 :    props%absplw = rval(:,:,:,1,:)
        1065             :
        1066           0 :    deallocate(rval)
        1067             :
        1068           0 :    ierr = pio_inq_varid(ncid, 'refindex_real_sw', vid)
        1069           0 :    ierr = pio_get_var(ncid, vid, props%refrtabsw)
        1070             :
        1071           0 :    ierr = pio_inq_varid(ncid, 'refindex_im_sw', vid)
        1072           0 :    ierr = pio_get_var(ncid, vid, props%refitabsw)
        1073             :
        1074           0 :    ierr = pio_inq_varid(ncid, 'refindex_real_lw', vid)
        1075           0 :    ierr = pio_get_var(ncid, vid, props%refrtablw)
        1076             :
        1077           0 :    ierr = pio_inq_varid(ncid, 'refindex_im_lw', vid)
        1078           0 :    ierr = pio_get_var(ncid, vid, props%refitablw)
        1079             :
        1080           0 :    ierr = pio_inq_varid(ncid, 'sigmag', vid)
        1081           0 :    ierr = pio_get_var(ncid, vid, props%sigmag)
        1082             :
        1083           0 :    ierr = pio_inq_varid(ncid, 'dgnum', vid)
        1084           0 :    ierr = pio_get_var(ncid, vid, props%dgnum)
        1085             :
        1086           0 :    ierr = pio_inq_varid(ncid, 'dgnumlo', vid)
        1087           0 :    ierr = pio_get_var(ncid, vid, props%dgnumlo)
        1088             :
        1089           0 :    ierr = pio_inq_varid(ncid, 'dgnumhi', vid)
        1090           0 :    ierr = pio_get_var(ncid, vid, props%dgnumhi)
        1091             :
        1092           0 :    ierr = pio_inq_varid(ncid, 'rhcrystal', vid)
        1093           0 :    ierr = pio_get_var(ncid, vid, props%rhcrystal)
        1094             :
        1095           0 :    ierr = pio_inq_varid(ncid, 'rhdeliques', vid)
        1096           0 :    ierr = pio_get_var(ncid, vid, props%rhdeliques)
        1097             :
        1098           0 : end subroutine modal_optics_init
        1099             :
        1100             : !================================================================================================
        1101             :
        1102           0 : subroutine bulk_props_init(physprop, nc_id)
        1103             :
        1104             : !  Read props for bulk aerosols
        1105             :
        1106             :    type (physprop_type), intent(inout) :: physprop ! storage for file data
        1107             :    type (file_desc_T),   intent(inout) :: nc_id    ! indentifier for netcdf file
        1108             :
        1109             :    ! Local variables
        1110             :    integer :: ierr
        1111             :
        1112             :    type(var_desc_T) :: vid
        1113             :
        1114             :    logical :: debug = .false.
        1115             :
        1116             :    character(len=*), parameter :: subname = 'bulk_props_init'
        1117             :    !------------------------------------------------------------------------------------
        1118             :
        1119             :    ! read microphys
        1120           0 :    ierr = pio_inq_varid(nc_id, 'name', vid)
        1121           0 :    ierr = pio_get_var(nc_id, vid, physprop%aername)
        1122             :
        1123             :    ! use GLC function to remove trailing nulls and blanks.
        1124             :    ! physprop%aername = aername_str(:GLC(aername_str))
        1125             :
        1126           0 :    ierr = pio_inq_varid(nc_id, 'density', vid)
        1127           0 :    ierr = pio_get_var(nc_id, vid, physprop%density_aer)
        1128             :
        1129           0 :    ierr = pio_inq_varid(nc_id, 'sigma_logr', vid)
        1130           0 :    ierr = pio_get_var(nc_id, vid, physprop%dispersion_aer)
        1131             :
        1132           0 :    ierr = pio_inq_varid(nc_id, 'dryrad', vid)
        1133           0 :    ierr = pio_get_var(nc_id, vid, physprop%dryrad_aer)
        1134             :
        1135           0 :    ierr = pio_inq_varid(nc_id, 'hygroscopicity', vid)
        1136           0 :    ierr = pio_get_var(nc_id, vid, physprop%hygro_aer)
        1137             :
        1138           0 :    ierr = pio_inq_varid(nc_id, 'num_to_mass_ratio', vid)
        1139           0 :    ierr = pio_get_var(nc_id, vid, physprop%num_to_mass_aer)
        1140             :
        1141             :    ! Output select data to log file
        1142           0 :    if (debug .and. masterproc .and. idx_sw_diag > 0) then
        1143           0 :       if (trim(physprop%aername) == 'SULFATE') then
        1144           0 :          write(iulog, '(2x, a)') '_______ hygroscopic growth in visible band _______'
        1145           0 :          call aer_optics_log_rh('SO4', physprop%sw_hygro_ext(:,idx_sw_diag), &
        1146           0 :             physprop%sw_hygro_ssa(:,idx_sw_diag), physprop%sw_hygro_asm(:,idx_sw_diag))
        1147             :       end if
        1148           0 :       write(iulog, *) subname//': finished for ', trim(physprop%aername)
        1149             :    end if
        1150             :
        1151           0 : end subroutine bulk_props_init
        1152             :
        1153             : !================================================================================================
        1154             :
        1155           0 : function exp_interpol(x, f, y) result(g)
        1156             : ! Purpose:
        1157             : !   interpolates f(x) to point y
        1158             : !   assuming f(x) = f(x0) exp a(x - x0)
        1159             : !   where a = ( ln f(x1) - ln f(x0) ) / (x1 - x0)
        1160             : !   x0 <= x <= x1
        1161             : !   assumes x is monotonically increasing
        1162             : ! Author: D. Fillmore
        1163             :
        1164             :    implicit none
        1165             :
        1166             :    real(r8), intent(in), dimension(:) :: x  ! grid points
        1167             :    real(r8), intent(in), dimension(:) :: f  ! grid function values
        1168             :    real(r8), intent(in) :: y                ! interpolation point
        1169             :    real(r8) :: g                            ! interpolated function value
        1170             :
        1171             :    integer :: k  ! interpolation point index
        1172             :    integer :: n  ! length of x
        1173             :    real(r8) :: a
        1174             :
        1175           0 :    n = size(x)
        1176             :
        1177             :    ! find k such that x(k) < y =< x(k+1)
        1178             :    ! set k = 1 if y <= x(1)  and  k = n-1 if y > x(n)
        1179             :
        1180           0 :    if (y <= x(1)) then
        1181             :      k = 1
        1182           0 :    else if (y >= x(n)) then
        1183           0 :      k = n - 1
        1184             :    else
        1185             :      k = 1
        1186           0 :      do while (y > x(k+1) .and. k < n)
        1187           0 :        k = k + 1
        1188             :      end do
        1189             :    end if
        1190             :
        1191             :    ! interpolate
        1192           0 :    a = (  log( f(k+1) / f(k) )  ) / ( x(k+1) - x(k) )
        1193           0 :    g = f(k) * exp( a * (y - x(k)) )
        1194             :    return
        1195             : end function exp_interpol
        1196             :
        1197             : !================================================================================================
        1198             :
        1199           0 : function lin_interpol(x, f, y) result(g)
        1200             : ! Purpose:
        1201             : !   interpolates f(x) to point y
        1202             : !   assuming f(x) = f(x0) + a * (x - x0)
        1203             : !   where a = ( f(x1) - f(x0) ) / (x1 - x0)
        1204             : !   x0 <= x <= x1
        1205             : !   assumes x is monotonically increasing
        1206             : ! Author: D. Fillmore
        1207             :
        1208             :    implicit none
        1209             :
        1210             :    real(r8), intent(in), dimension(:) :: x  ! grid points
        1211             :    real(r8), intent(in), dimension(:) :: f  ! grid function values
        1212             :    real(r8), intent(in) :: y                ! interpolation point
        1213             :    real(r8) :: g                            ! interpolated function value
        1214             :
        1215             :    integer :: k  ! interpolation point index
        1216             :    integer :: n  ! length of x
        1217             :    real(r8) :: a
        1218             :
        1219           0 :    n = size(x)
        1220             :
        1221             :    ! find k such that x(k) < y =< x(k+1)
        1222             :    ! set k = 1 if y <= x(1)  and  k = n-1 if y > x(n)
        1223             :
        1224           0 :    if (y <= x(1)) then
        1225             :      k = 1
        1226           0 :    else if (y >= x(n)) then
        1227           0 :      k = n - 1
        1228             :    else
        1229             :      k = 1
        1230           0 :      do while (y > x(k+1) .and. k < n)
        1231           0 :        k = k + 1
        1232             :      end do
        1233             :    end if
        1234             :
        1235             :    ! interpolate
        1236           0 :    a = (  f(k+1) - f(k) ) / ( x(k+1) - x(k) )
        1237           0 :    g = f(k) + a * (y - x(k))
        1238             :    return
        1239             : end function lin_interpol
        1240             :
        1241             : !================================================================================================
        1242             :
        1243             : subroutine aer_optics_log(name, ext, ssa, asm)
        1244             :
        1245             :    ! Purpose:
        1246             :    !   write aerosol optical constants to log file
        1247             :
        1248             :    ! Author: D. Fillmore
        1249             :
        1250             :    character(len=*), intent(in) :: name
        1251             :    real(r8), intent(in) :: ext(:)
        1252             :    real(r8), intent(in) :: ssa(:)
        1253             :    real(r8), intent(in) :: asm(:)
        1254             :
        1255             :    integer :: kbnd, nbnd
        1256             :    !------------------------------------------------------------------------------------
        1257             :
        1258             :    nbnd = ubound(ext, 1)
        1259             :
        1260             :    write(iulog, '(2x, a)') name
        1261             :    write(iulog, '(2x, a, 4x, a, 4x, a, 4x, a)') 'SW band', 'ext (m^2 kg^-1)', ' ssa', ' asm'
        1262             :    do kbnd = 1, nbnd
        1263             :       write(iulog, '(2x, i7, 4x, f13.2, 4x, f4.2, 4x, f4.2)') kbnd, ext(kbnd), ssa(kbnd), asm(kbnd)
        1264             :    end do
        1265             :
        1266             : end subroutine aer_optics_log
        1267             :
        1268             : !================================================================================================
        1269             :
        1270             :
        1271           0 : subroutine aer_optics_log_rh(name, ext, ssa, asm)
        1272             :
        1273             :    ! Purpose:
        1274             :    !   write out aerosol optical properties
        1275             :    !   for a set of test rh values
        1276             :    !   to test hygroscopic growth interpolation
        1277             :
        1278             :    ! Author: D. Fillmore
        1279             :
        1280             :    character(len=*), intent(in) :: name
        1281             :    real(r8), intent(in) :: ext(nrh)
        1282             :    real(r8), intent(in) :: ssa(nrh)
        1283             :    real(r8), intent(in) :: asm(nrh)
        1284             :
        1285             :    integer :: krh_test
        1286             :    integer, parameter :: nrh_test = 36
        1287             :    integer :: krh
        1288             :    real(r8) :: rh
        1289             :    real(r8) :: rh_test(nrh_test)
        1290             :    real(r8) :: exti
        1291             :    real(r8) :: ssai
        1292             :    real(r8) :: asmi
        1293             :    real(r8) :: wrh
        1294             :    !------------------------------------------------------------------------------------
        1295             :
        1296           0 :    do krh_test = 1, nrh_test
        1297           0 :       rh_test(krh_test) = sqrt(sqrt(sqrt(sqrt(((krh_test - 1.0_r8) / (nrh_test - 1))))))
        1298             :    enddo
        1299           0 :    write(iulog, '(2x, a)') name
        1300           0 :    write(iulog, '(2x, a, 4x, a, 4x, a, 4x, a)') '   rh', 'ext (m^2 kg^-1)', '  ssa', '  asm'
        1301             :
        1302             :    ! loop through test rh values
        1303           0 :    do krh_test = 1, nrh_test
        1304             :       ! find corresponding rh index
        1305           0 :       rh = rh_test(krh_test)
        1306           0 :       krh = min(floor( (rh) * nrh ) + 1, nrh - 1)
        1307           0 :       wrh = (rh) *nrh - krh
        1308           0 :       exti = ext(krh + 1) * (wrh + 1) - ext(krh) * wrh
        1309           0 :       ssai = ssa(krh + 1) * (wrh + 1) - ssa(krh) * wrh
        1310           0 :       asmi = asm(krh + 1) * (wrh + 1) - asm(krh) * wrh
        1311           0 :       write(iulog, '(2x, f5.3, 4x, f13.3, 4x, f5.3, 4x, f5.3)') rh_test(krh_test), exti, ssai, asmi
        1312             :    end do
        1313             :
        1314           0 : end subroutine aer_optics_log_rh
        1315             :
        1316             :
        1317             : !================================================================================================
        1318             :
        1319           0 : end module phys_prop

Generated by: LCOV version 1.14