radiation_aerosols.f

!>  \file radiation_aerosols.f
!!  This file contains climatological atmospheric aerosol schemes for
!!  radiation computations

!  ==========================================================  !!!!!
!            'module_radiation_aerosols' description           !!!!!
!  ==========================================================  !!!!!
!                                                                      !
!   this module contains climatological atmospheric aerosol schemes for!
!   radiation computations.                                            !
!                                                                      !
!   in the module, the externally callable subroutines are :           !
!                                                                      !
!      'aer_init'   -- initialization                                  !
!         inputs:                                                      !
!           ( NLAY, me )                                               !
!         outputs:                                                     !
!           ( errflg, errmsg )                                         !
!                                                                      !
!      'aer_update' -- updating aerosol data                           !
!         inputs:                                                      !
!           ( iyear, imon, me )                                        !
!         outputs:                                                     !
!           ( errflg, errmsg )                                         !
!                                                                      !
!      'setaer'     -- mapping aeros profile, compute aeros opticals   !
!         inputs:                                                      !
!           (prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,aerfld,xlon,xlat, !
!            IMAX,NLAY,NLP1, lsswr,lslwr,                              !
!         outputs:                                                     !
!          (aerosw,aerolw,aerodp,errmsg,errflg)                        !
!                                                                      !
!                                                                      !
!   external modules referenced:                                       !
!       'module module_radsw_parameters' in 'radsw_xxxx#_param.f'      !
!       'module module_radlw_parameters' in 'radlw_xxxx#_param.f'      !
!       'module module_radlw_cntr_para'  in 'radsw_xxxx#_param.f'      !
!                                                                      !
!   output variable definitions:                                       !
!       aerosw(IMAX,NLAY,NBDSW,1) - aerosols optical depth for sw      !
!       aerosw(IMAX,NLAY,NBDSW,2) - aerosols single scat albedo for sw !
!       aerosw(IMAX,NLAY,NBDSW,3) - aerosols asymmetry parameter for sw!
!                                                                      !
!       aerolw(IMAX,NLAY,NBDLW,1) - aerosols optical depth for lw      !
!       aerolw(IMAX,NLAY,NBDLW,2) - aerosols single scattering albedo  !
!       aerolw(IMAX,NLAY,NBDLW,3) - aerosols asymmetry parameter       !
!                                                                      !
!                                                                      !
!   program history:                                                   !
!     apr     2003  ---  y.-t. hou     created                         !
!     nov 04, 2003  ---  y.-t. hou     modified version                !
!     apr 15, 2005  ---  y.-t. hou     modified module structure       !
!     jul     2006  ---  y.-t. hou     add volcanic forcing            !
!     feb     2007  ---  y.-t. hou     add generalized spectral band   !
!                   interpolation for sw aerosol optical properties    !
!     mar     2007  ---  y.-t. hou     add generalized spectral band   !
!                   interpolation for lw aerosol optical properties    !
!     aug     2007  ---  y.-t. hou     change clim-aer vert domain     !
!                   from pressure reference to sigma reference         !
!     sep     2007  ---  y.-t. hou     moving temporary allocatable    !
!                   module variable arrays to subroutine dynamically   !
!                   allocated arrays (eliminate deallocate calls)      !
!     jan     2010  ---  sarah lu      add gocart option               !
!     may     2010  ---  sarah lu      add geos4-gocart climo          !
!     jul     2010  --   s. moorthi - merged NEMS version with new GFS !
!                        version                                       !
!     oct 23, 2010  ---  Hsin-mu lin   modified subr setclimaer to     !
!        interpolate the 5 degree aerosol data to small domain based on!
!        the nearby 4 points instead of previous nearby assignment by  !
!        using the 5 degree data. this process will eliminate the dsw  !
!        jagged edges in the east conus where aerosol effect are lagre.!
!     dec     2010  ---  y.-t. hou     modified and optimized bi-linear!
!        horizontal interpolation in subr setclimaer. added safe guard !
!        measures in lat/lon indexing and added sea/land mask variable !
!        slmsk as input field to help aerosol profile selection.       !
!     jan     2011  ---  y.-t. hou     divided the program into two    !
!        separated interchangeable modules: a climatology aerosol      !
!        module, and a gocart aerosol scheme module. the stratospheric !
!        volcanic aerosol part is still within the two driver modules, !
!        and may also become a separate one in the further development.!
!        unified in/out argument list for both clim and gocart types of!
!        schemes and added vertically integrated aer-opt-dep, aerodp,  !
!        to replace tau_gocart as optional output for various species. !
!     aug     2012  ---  y.-t. hou     changed the initialization subr !
!        'aerinit' into two parts: 'aer_init' is called at the start   !
!        of run to set up module parameters; and 'aer_update' is       !
!        called within the time loop to check and update data sets.    !
!     nov     2012  ---  y.-t. hou     modified control parameters thru!
!        module 'physparam'.                                            !
!     jan     2013  ---  sarah lu and y.-t. hou   reintegrate both     !
!        opac-clim and gocart schemes into one module to make the      !
!        program best utilize common components. added aerosol model   !
!        scheme selection control variable iaer_mdl to the namelist.   !
!      Aug     2013  --- s. moorthi - merge sarah's gocart changes with!
!                                     yutai's changes                  !
!      13Feb2014  --- Sarah lu - compute aod at 550nm                  !
!      jun     2018  --- h-m lin and y-t hou   updated spectral band   !
!        mapping method for aerosol optical properties. controled by   !
!        internal variable lmap_new through namelist variable iaer.    !
!      may     2019  --- sarah lu, restore the gocart option, allowing !
!        aerosol ext, ssa, asy determined from MERRA2 monthly climo    !
!        with new spectral band mapping method                         !
!                                                                      !
!   references for opac climatological aerosols:                       !
!     hou et al. 2002  (ncep office note 441)                          !
!     hess et al. 1998 - bams v79 831-844                              !
!                                                                      !
!   references for gocart interactive aerosols:                        !
!     chin et al., 2000 - jgr, v105, 24671-24687                       !
!     colarco et al., 2010 - jgr, v115, D14207                         !
!                                                                      !
!   references for merra2 aerosol reanalysis:                          !
!     randles et al., 2017 - jclim, v30, 6823-6850                     !
!     buchard et al., 2017 - jclim, v30, 6851-6871                     !
!                                                                      !
!   references for stratosperic volcanical aerosols:                   !
!     sato et al. 1993 - jgr, v98, d12, 22987-22994                    !
!                                                                      !
!                                                                      !
!!!!!  ==========================================================  !!!!!
!!!!!                       end descriptions                       !!!!!
!!!!!  ==========================================================  !!!!!


!> This module contains climatological atmospheric aerosol schemes for
!! radiation computations.
      module module_radiation_aerosols
!
      use machine,  only : kind_phys, kind_io4, kind_io8
      use module_iounitdef,        only : NIAERCM
      use module_radsw_parameters, only : NBDSW,  wvnsw1=>wvnum1,       &
     &                                    NSWSTR, wvnsw2=>wvnum2
      use module_radlw_parameters, only : NBDLW,  wvnlw1, wvnlw2
!
      use funcphys,                     only : fpkap
      use aerclm_def,                   only : ntrcaerm

!
      implicit   none
!
      private

!  ---  version tag and last revision date
      character(40), parameter ::                                       &
     &   VTAGAER='NCEP-Radiation_aerosols  v5.2  Jan 2013 '
!    &   VTAGAER='NCEP-Radiation_aerosols  v5.1  Nov 2012 '
!    &   VTAGAER='NCEP-Radiation_aerosols  v5.0  Aug 2012 '

! --- general use parameter constants:
!> num of output fields for SW rad
      integer, parameter, public :: NF_AESW = 3
!> num of output fields for LW rad
      integer, parameter, public :: NF_AELW = 3
!> starting band number in ir region
      integer, parameter, public :: NLWSTR  = 1
!> num of species for output aod (opnl)
      integer, parameter, public :: NSPC    = 5
!> total+species
      integer, parameter, public :: NSPC1   = NSPC + 1

      real (kind=kind_phys), parameter :: f_zero = 0.0
      real (kind=kind_phys), parameter :: f_one  = 1.0

!  ---  module control parameters set in subroutine "aer_init"
!> number of actual bands for sw aerosols; calculated according to
!! laswflg setting
      integer, save :: NSWBND  = NBDSW
!> number of actual bands for lw aerosols; calculated according to
!! lalwflg and lalw1bd settings
      integer, save :: NLWBND  = NBDLW
!> total number of bands for sw+lw aerosols
      integer, save :: NSWLWBD = NBDSW+NBDLW
!  LW aerosols effect control flag
!    =.true.:aerosol effect is included in LW radiation
!    =.false.:aerosol effect is not included in LW radiation
      logical, save :: lalwflg = .true.
!  SW aerosols effect control flag
!    =.true.:aerosol effect is included in SW radiation
!    =.false.:aerosol effect is not included in SW radiation
      logical, save :: laswflg = .true.
!  stratospheric volcanic aerosol effect flag
!    =.true.:historical events of stratosphere volcanic aerosol effect
!            is included radiation (limited by data availability)
!    =.false.:volcanic aerosol effect is not included in radiation
      logical, save :: lavoflg = .true.

      logical, save :: lmap_new = .true.  ! use new mapping method (set in aer_init)

! --------------------------------------------------------------------- !
!   section-1 : module variables for spectral band interpolation        !
!               similar to gfdl-sw treatment (2000 version)             !
! --------------------------------------------------------------------- !

!  ---  parameter constants:
!> num of wvnum regions where solar flux is constant
      integer, parameter, public :: NWVSOL  = 151

!> total num of wvnum included
      integer, parameter, public :: NWVTOT  = 57600
!> total num of wvnum in ir range
      integer, parameter, public :: NWVTIR  = 4000

!> number of wavenumbers in each region where the solar flux is constant
      integer, dimension(NWVSOL), save :: nwvns0

      data nwvns0   / 100,  11,  14,  18,  24,  33,  50,  83,  12,  12, &
     &  13,  15,  15,  17,  18,  20,  21,  24,  26,  30,  32,  37,  42, &
     &  47,  55,  64,  76,  91, 111, 139, 179, 238, 333,  41,  42,  45, &
     &  46,  48,  51,  53,  55,  58,  61,  64,  68,  71,  75,  79,  84, &
     &  89,  95, 101, 107, 115, 123, 133, 142, 154, 167, 181, 197, 217, &
     & 238, 263, 293, 326, 368, 417, 476, 549, 641, 758, 909, 101, 103, &
     & 105, 108, 109, 112, 115, 117, 119, 122, 125, 128, 130, 134, 137, &
     & 140, 143, 147, 151, 154, 158, 163, 166, 171, 175, 181, 185, 190, &
     & 196, 201, 207, 213, 219, 227, 233, 240, 248, 256, 264, 274, 282, &
     & 292, 303, 313, 325, 337, 349, 363, 377, 392, 408, 425, 444, 462, &
     & 483, 505, 529, 554, 580, 610, 641, 675, 711, 751, 793, 841, 891, &
     & 947,1008,1075,1150,1231,1323,1425,1538,1667,1633,14300 /

!> solar flux \f$w/m^2\f$ in each wvnumb region where it is constant
      real (kind=kind_phys), dimension(NWVSOL), save :: s0intv

      data  s0intv(  1: 50)       /                                     &
     &     1.60000E-6, 2.88000E-5, 3.60000E-5, 4.59200E-5, 6.13200E-5,  &
     &     8.55000E-5, 1.28600E-4, 2.16000E-4, 2.90580E-4, 3.10184E-4,  &
     &     3.34152E-4, 3.58722E-4, 3.88050E-4, 4.20000E-4, 4.57056E-4,  &
     &     4.96892E-4, 5.45160E-4, 6.00600E-4, 6.53600E-4, 7.25040E-4,  &
     &     7.98660E-4, 9.11200E-4, 1.03680E-3, 1.18440E-3, 1.36682E-3,  &
     &     1.57560E-3, 1.87440E-3, 2.25500E-3, 2.74500E-3, 3.39840E-3,  &
     &     4.34000E-3, 5.75400E-3, 7.74000E-3, 9.53050E-3, 9.90192E-3,  &
     &     1.02874E-2, 1.06803E-2, 1.11366E-2, 1.15830E-2, 1.21088E-2,  &
     &     1.26420E-2, 1.32250E-2, 1.38088E-2, 1.44612E-2, 1.51164E-2,  &
     &     1.58878E-2, 1.66500E-2, 1.75140E-2, 1.84450E-2, 1.94106E-2 /
      data  s0intv( 51:100)       /                                     &
     &     2.04864E-2, 2.17248E-2, 2.30640E-2, 2.44470E-2, 2.59840E-2,  &
     &     2.75940E-2, 2.94138E-2, 3.13950E-2, 3.34800E-2, 3.57696E-2,  &
     &     3.84054E-2, 4.13490E-2, 4.46880E-2, 4.82220E-2, 5.22918E-2,  &
     &     5.70078E-2, 6.19888E-2, 6.54720E-2, 6.69060E-2, 6.81226E-2,  &
     &     6.97788E-2, 7.12668E-2, 7.27100E-2, 7.31610E-2, 7.33471E-2,  &
     &     7.34814E-2, 7.34717E-2, 7.35072E-2, 7.34939E-2, 7.35202E-2,  &
     &     7.33249E-2, 7.31713E-2, 7.35462E-2, 7.36920E-2, 7.23677E-2,  &
     &     7.25023E-2, 7.24258E-2, 7.20766E-2, 7.18284E-2, 7.32757E-2,  &
     &     7.31645E-2, 7.33277E-2, 7.36128E-2, 7.33752E-2, 7.28965E-2,  &
     &     7.24924E-2, 7.23307E-2, 7.21050E-2, 7.12620E-2, 7.10903E-2 /
      data  s0intv(101:151)       /                        7.12714E-2,  &
     &     7.08012E-2, 7.03752E-2, 7.00350E-2, 6.98639E-2, 6.90690E-2,  &
     &     6.87621E-2, 6.52080E-2, 6.65184E-2, 6.60038E-2, 6.47615E-2,  &
     &     6.44831E-2, 6.37206E-2, 6.24102E-2, 6.18698E-2, 6.06320E-2,  &
     &     5.83498E-2, 5.67028E-2, 5.51232E-2, 5.48645E-2, 5.12340E-2,  &
     &     4.85581E-2, 4.85010E-2, 4.79220E-2, 4.44058E-2, 4.48718E-2,  &
     &     4.29373E-2, 4.15242E-2, 3.81744E-2, 3.16342E-2, 2.99615E-2,  &
     &     2.92740E-2, 2.67484E-2, 1.76904E-2, 1.40049E-2, 1.46224E-2,  &
     &     1.39993E-2, 1.19574E-2, 1.06386E-2, 1.00980E-2, 8.63808E-3,  &
     &     6.52736E-3, 4.99410E-3, 4.39350E-3, 2.21676E-3, 1.33812E-3,  &
     &     1.12320E-3, 5.59000E-4, 3.60000E-4, 2.98080E-4, 7.46294E-5  /

      real (kind=kind_phys), dimension(NBDSW), save :: wvn_sw1, wvn_sw2
      real (kind=kind_phys), dimension(NBDLW), save :: wvn_lw1, wvn_lw2
! --------------------------------------------------------------------- !
!   section-2 : module variables for stratospheric volcanic aerosols    !
!               from historical data (sato et al. 1993)                 !
! --------------------------------------------------------------------- !

!  ---  parameter constants:
!> lower limit (year) data available
      integer, parameter :: MINVYR = 1850
!> upper limit (year) data available
      integer, parameter :: MAXVYR = 1999

!> monthly, 45-deg lat-zone aerosols data set in subroutine 'aer_init'
      integer, allocatable, save :: ivolae(:,:,:)

!  ---  static control variables:
!> starting year of data in the input file
      integer :: kyrstr
!> ending year of data in the input file
      integer :: kyrend
!> the year of data in use in the input file
      integer :: kyrsav
!> the month of data in use in the input file
      integer :: kmonsav

! --------------------------------------------------------------------- !
!   section-3 : module variables for opac climatological aerosols       !
!               optical properties (hess et al. 1989)                   !
! --------------------------------------------------------------------- !

!  ---  parameters and constants:
!> num of max componets in a profile
      integer, parameter :: NXC = 5
!> num of aerosols profile structures
      integer, parameter :: NAE = 7
!> num of atmos aerosols domains
      integer, parameter :: NDM = 5
!> num of lon-points in glb aeros data set
      integer, parameter :: IMXAE = 72
!> num of lat-points in glb aeros data set
      integer, parameter :: JMXAE = 37
!> num of bands for clim aer data (opac)
      integer, parameter :: NAERBND=61
!> num of rh levels for rh-dep components
      integer, parameter :: NRHLEV =8
!> num of rh independent aeros species
      integer, parameter :: NCM1 = 6
!> num of rh dependent aeros species
      integer, parameter :: NCM2 = 4
      integer, parameter :: NCM  = NCM1+NCM2

!> predefined relative humidity levels
      real (kind=kind_phys), dimension(NRHLEV), save :: rhlev
      data  rhlev (:) / 0.0, 0.5, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99 /

!  ---  the following arrays are for climatological data that are
!           allocated and read in subroutine 'clim_aerinit'.
!   - global aerosol distribution:
!      haer  (NDM,NAE)  - scale height of aerosols (km)
!      prsref(NDM,NAE)  - ref pressure lev (sfc to toa) in mb (100Pa)
!      sigref(NDM,NAE)  - ref sigma lev (sfc to toa)

!> scale height of aerosols (km)
      real (kind=kind_phys), save, dimension(NDM,NAE) :: haer
!> ref pressure lev (sfc to toa) in mb (100Pa)
      real (kind=kind_phys), save, dimension(NDM,NAE) :: prsref
!> ref sigma lev (sfc to toa)
      real (kind=kind_phys), save, dimension(NDM,NAE) :: sigref

!  ---  the following arrays are allocate and setup in subr 'clim_aerinit'
!   - for relative humidity independent aerosol optical properties:
!      species : insoluble        (inso); soot             (soot);
!                mineral nuc mode (minm); mineral acc mode (miam);
!                mineral coa mode (micm); mineral transport(mitr).
!      extrhi(NCM1,NSWLWBD) - extinction coefficient for sw+lw spectral band
!      scarhi(NCM1,NSWLWBD) - scattering coefficient for sw+lw spectral band
!      ssarhi(NCM1,NSWLWBD) - single scattering albedo for sw+lw spectral band
!      asyrhi(NCM1,NSWLWBD) - asymmetry parameter for sw+lw spectral band
!   - for relative humidity dependent aerosol optical properties:
!      species : water soluble    (waso); sea salt acc mode(ssam);
!                sea salt coa mode(sscm); sulfate droplets (suso).
!      rh level: 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99%
!      extrhd(NRHLEV,NCM2,NSWLWBD) - extinction coefficient for sw+lw band
!      scarhd(NRHLEV,NCM2,NSWLWBD) - scattering coefficient for sw+lw band
!      ssarhd(NRHLEV,NCM2,NSWLWBD) - single scattering albedo for sw+lw band
!      asyrhd(NRHLEV,NCM2,NSWLWBD) - asymmetry parameter for sw+lw band
!   - for stratospheric aerosols optical properties:
!      extstra(NSWLWBD)            - extinction coefficient for sw+lw band

      real (kind=kind_phys), allocatable, save, dimension(:,:)   ::     &
     &       extrhi, scarhi, ssarhi, asyrhi
      real (kind=kind_phys), allocatable, save, dimension(:,:,:) ::     &
     &       extrhd, scarhd, ssarhd, asyrhd
      real (kind=kind_phys), allocatable, save, dimension(:)     ::     &
     &       extstra

!  ---  the following arrays are calculated in subr 'clim_aerinit'
!   - for topospheric aerosol profile distibution:
!      kprfg (    IMXAE*JMXAE)   - aeros profile index
!      idxcg (NXC*IMXAE*JMXAE)   - aeros component index
!      cmixg (NXC*IMXAE*JMXAE)   - aeros component mixing ratio
!      denng ( 2 *IMXAE*JMXAE)   - aerosols number density

!> \name topospheric aerosol profile distribution

!> aeros component mixing ratio
      real (kind=kind_phys), dimension(NXC,IMXAE,JMXAE), save :: cmixg
!> aeros number density
      real (kind=kind_phys), dimension( 2 ,IMXAE,JMXAE), save :: denng
!> aeros component index
      integer,               dimension(NXC,IMXAE,JMXAE), save :: idxcg
!> aeros profile index
      integer,               dimension(    IMXAE,JMXAE), save :: kprfg

! --------------------------------------------------------------------- !
!   section-4 : module variables for gocart aerosol optical properties  !
! --------------------------------------------------------------------- !
!> \name module variables for gocart aerosol optical properties

!  ---  parameters and constants:
!> num of bands for aer data (gocart)
      integer, parameter :: KAERBNDD=61
      integer, parameter :: KAERBNDI=56
!> num of rh levels for rh-dep components
      integer, parameter :: KRHLEV =36
!> num of gocart rh indep aerosols
      integer, parameter :: KCM1 = 5
!> num of gocart rh dep aerosols
      integer, parameter :: KCM2 = 10
!> num of gocart aerosols
      integer, parameter :: KCM  = KCM1 + KCM2

      real (kind=kind_phys), dimension(KRHLEV) :: rhlev_grt             &
      data  rhlev_grt (:)/ .00, .05, .10, .15, .20, .25, .30, .35,      &
     &      .40, .45, .50, .55, .60, .65, .70, .75, .80, .81, .82,      &
     &      .83, .84, .85, .86, .87, .88, .89, .90, .91, .92, .93,      &
     &      .94, .95, .96, .97, .98, .99 /

!> \name relative humidity independent aerosol optical properties:
!! species: du001, du002, du003, du004, du005
!    extrhi_grt(KCM1,NSWLWBD) - extinction coefficient for sw+lw band
!    scarhi_grt(KCM1,NSWLWBD) - scattering coefficient for sw+lw band
!    ssarhi_grt(KCM1,NSWLWBD) - single scattering albedo for sw+lw band
!    asyrhi_grt(KCM1,NSWLWBD) - asymmetry parameter for sw+lw band
      real (kind=kind_phys),allocatable,save,dimension(:,:) ::          &
     &   extrhi_grt, extrhi_grt_550, scarhi_grt, ssarhi_grt, asyrhi_grt
!
!> \name relative humidity dependent aerosol optical properties:
!! species : ss001, ss002, ss003, ss004, ss005, so4,
!!           bcphobic, bcphilic, ocphobic, ocphilic
!    extrhd_grt(KRHLEV,KCM2,NSWLWBD) - extinction coefficient for sw+lw band
!    scarhd_grt(KRHLEV,KCM2,NSWLWBD) - scattering coefficient for sw+lw band
!    ssarhd_grt(KRHLEV,KCM2,NSWLWBD) - single scattering albedo for sw+lw band
!    asyrhd_grt(KRHLEV,KCM2,NSWLWBD) - asymmetry parameter for sw+lw band

!> extinction coefficient
      real (kind=kind_phys),allocatable,save,dimension(:,:,:) ::        &
     &   extrhd_grt, extrhd_grt_550, scarhd_grt, ssarhd_grt, asyrhd_grt

!> gocart species
      integer, parameter         :: num_gc = 5
      character*2                :: gridcomp(num_gc)
      integer, dimension (num_gc):: num_radius, radius_lower
      integer, dimension (num_gc):: trc_to_aod

      data gridcomp     /'DU', 'SS', 'SU', 'BC', 'OC'/
      data num_radius   /5, 5, 1, 2, 2 /
      data radius_lower /1, 6, 11, 12, 14 /
      data trc_to_aod   /1, 5, 4, 2, 3/   ! dust, soot, waso, suso, ssam

! =======================================================================
! --------------------------------------------------------------------- !
!   section-5 : module variables for aod diagnostic                     !
! --------------------------------------------------------------------- !
!! ---  the following are for diagnostic purpose to output aerosol optical depth
!       aod from 10 components are grouped into 5 major different species:
!      1:dust (inso,minm,miam,micm,mitr); 2:black carbon (soot)
!      3:water soluble (waso);            4:sulfate (suso);      5:sea salt (ssam,sscm)
!
!      idxspc (NCM)         - index conversion array
!      lspcaod              - logical flag for aod from individual species
!
!> index conversion array:data  idxspc / 1, 2, 1, 1, 1, 1, 3, 5, 5, 4 /
      integer, dimension(NCM) :: idxspc
      data  idxspc / 1, 2, 1, 1, 1, 1, 3, 5, 5, 4 /
!
!   - wvn550 is the wavenumber (1/cm) of wavelenth 550nm for diagnostic aod output
!     nv_aod is the sw spectral band covering wvn550 (comp in aer_init)
!
!> the wavenumber (\f$cm^-1\f$) of wavelength 550nm for diagnostic aod output
      real (kind=kind_phys), parameter :: wvn550 = 1.0e4/0.55
!> the sw spectral band covering wvn550 (comp in aer_init)
      integer, save      :: nv_aod = 1
      integer            :: i550,id550

!  ---  public interfaces

      public aer_init, aer_update, setaer

! =================
      contains
! =================

!> \defgroup mod_radiation_aerosols Radiation Aerosols Module
!>  This module contains climatological atmospheric aerosol schemes for
!!  radiation computations.
!!
!!\version NCEP-Radiation_aerosols  v5.2  Jan 2013
!!
!!\n This module has three externally callable subroutines:
!! - aer_init() -- initialization; called at the start of run to set up
!!                 module parameters
!! - aer_update() -- updating aerosol data; called within the time loop
!!                   to check and update data sets
!! - setaer() -- mapping aeros profile, compute aeros opticals
!!
!!\n References:
!! - OPAC climatological aerosols:
!! Hou et al. 2002 \cite hou_et_al_2002; Hess et al. 1998
!! \cite hess_et_al_1998
!! - GOCART interactive aerosols:
!! Chin et al., 2000 \cite chin_et_al_2000
!! Colarco et al., 2010 - jgr, v115, D14207\cite colarco_et_al_2010
!!
!! - MERRA2 aerosol reanalysis:
!! Randles et al., 2017 - jclim, v30, 6823-6850\cite randles_et_al_2017
!! Buchard et al., 2017 - jclim, v30, 6851-6871\cite buchard_et_al_2017
!!
!! - Stratospheric volcanical aerosols:
!! Sato et al. 1993 \cite sato_et_al_1993
!!
!> The initialization program is to set up necessary parameters and
!! working arrays.
!!
!>\param NLAY    number of model vertical layers (not used)
!>\param me      print message control flag
!>\section gen_al General Algorithm
!-----------------------------------
      subroutine aer_init                                               &
     &     ( NLAY, me, iaermdl, iaerflg, lalw1bd, aeros_file, con_pi,   &
     &     con_t0c, con_c, con_boltz, con_plnk, errflg, errmsg)

!  ==================================================================  !
!                                                                      !
!  aer_init is the initialization program to set up necessary          !
!    parameters and working arrays.                                    !
!                                                                      !
!  inputs:                                                             !
!     NLAY    - number of model vertical layers  (not used)            !
!     me      - print message control flag                             !
!     iaermdl - tropospheric aerosol model scheme flag                 !
!               =0 opac-clim; =1 gocart-clim, =2 gocart-prognostic     !
!               =5 opac-clim new spectral mapping                      !
!     lalw1bd = logical lw aeros propty 1 band vs multi-band cntl flag !
!               =t use 1 broad band optical property                   !
!               =f use multi bands optical property                    !
!                                                                      !
!  outputs: (CCPP error handling)                                      !
!     errmsg  - CCPP error message                                     !
!     errflg  - CCPP error flag                                        !
!                                                                      !
!  internal module variables:                                          !
!     lalwflg - logical lw aerosols effect control flag                !
!               =t compute lw aerosol optical prop                     !
!     laswflg - logical sw aerosols effect control flag                !
!               =t compute sw aerosol optical prop                     !
!     lavoflg - logical stratosphere volcanic aerosol control flag     !
!               =t include volcanic aerosol effect                     !
!                                                                      !
!  internal module constants:                                          !
!     NWVSOL  - num of wvnum regions where solar flux is constant      !
!     NWVTOT  - total num of wave numbers used in sw spectrum          !
!     NWVTIR  - total num of wave numbers used in the ir region        !
!     NSWBND  - total number of sw spectral bands                      !
!     NLWBND  - total number of lw spectral bands                      !
!                                                                      !
!  usage:    call aer_init                                             !
!                                                                      !
!  subprograms called:  clim_aerinit, gocart_aerinit,                  !
!                       wrt_aerlog, set_volcaer, set_spectrum,         !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
      integer,          intent(in) :: NLAY, me, iaermdl, iaerflg
      logical,          intent(in) :: lalw1bd
      character(len=26),intent(in) :: aeros_file
      real(kind_phys),  intent(in) :: con_pi,con_t0c, con_c, con_boltz, & 
     &     con_plnk
!  ---  output:
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      real (kind=kind_phys), dimension(NWVTOT) :: solfwv        ! one wvn sol flux
      real (kind=kind_phys), dimension(NWVTIR) :: eirfwv        ! one wvn ir flux
!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      kyrstr  = 1
      kyrend  = 1
      kyrsav  = 1
      kmonsav = 1

      laswflg= (mod(iaerflg,10) > 0)    ! control flag for sw tropospheric aerosol
      lalwflg= (mod(iaerflg/10,10) > 0) ! control flag for lw tropospheric aerosol
      lavoflg= (mod(iaerflg/100,10) >0) ! control flag for stratospheric volcanic aeros

!> -# Call wrt_aerlog to write aerosol parameter configuration to output logs.

      if ( me == 0 ) then

        call wrt_aerlog(iaermdl, iaerflg, lalw1bd, errflg, errmsg)      ! write aerosol param info to log file
!  ---  inputs:   (in scope variables)
!  ---  outputs:  (CCPP error handling)

      endif

      if ( iaerflg == 0 ) return      ! return without any aerosol calculations

!  --- ...  in sw, aerosols optical properties are computed for each radiation
!           spectral band; while in lw, optical properties can be calculated
!           for either only one broad band or for each of the lw radiation bands

      if ( laswflg ) then
        NSWBND = NBDSW
      else
        NSWBND = 0
      endif

      if ( lalwflg ) then
        if ( lalw1bd ) then
          NLWBND = 1
        else
          NLWBND = NBDLW
        endif
      else
        NLWBND = 0
      endif

      NSWLWBD = NSWBND + NLWBND

      wvn_sw1(:) = wvnsw1(:)
      wvn_sw2(:) = wvnsw2(:)
      wvn_lw1(:) = wvnlw1(:)
      wvn_lw2(:) = wvnlw2(:)

! note: for result consistency, the defalt opac-clim aeros setting still use
!       old spectral band mapping. use iaermdl=5 to use new mapping method

      if ( iaermdl == 0 ) then                    ! opac-climatology scheme
        lmap_new = .false.

        wvn_sw1(2:NBDSW-1) = wvn_sw1(2:NBDSW-1) + 1
        wvn_lw1(2:NBDLW) = wvn_lw1(2:NBDLW) + 1
      else
        lmap_new = .true.
      endif

      if ( iaerflg /= 100 ) then

!> -# Call set_spectrum to set up spectral one wavenumber solar/IR
!! fluxes.

        call set_spectrum(con_pi, con_t0c, con_c, con_boltz, con_plnk,  &
     &        errflg, errmsg)
!  ---  inputs:   (module constants)
!  ---  outputs:  (ccpp error handling)

!> -# Call clim_aerinit() to invoke tropospheric aerosol initialization.

        if ( iaermdl==0 .or. iaermdl==5 ) then      ! opac-climatology scheme
          call clim_aerinit                                             &
!  ---  inputs:
     &     ( solfwv, eirfwv, me, aeros_file,                            &
!  ---  outputs:
     &     errflg, errmsg)

        elseif ( iaermdl==1 .or. iaermdl==2 ) then  ! gocart clim/prog scheme

          call gocart_aerinit                                           &
!  ---  inputs:
     &     ( solfwv, eirfwv, me,                                        &
!  ---  outputs:
     &     errflg, errmsg)

        else
          if ( me == 0 ) then
            print *,'  !!! ERROR in aerosol model scheme selection',    &
     &              ' iaermdl =',iaermdl
            errflg = 1
            errmsg = 'ERROR(aer_init): aerosol model scheme selected'// &
     &           'is invalid'
            return
          endif
        endif

      endif    ! end if_iaerflg_block

!> -# Call set_volcaer to invoke stratospheric volcanic aerosol
!! initialization.

      if ( lavoflg ) then

        call set_volcaer(errflg, errmsg)
!  ---  inputs:  (module variables)
!  ---  outputs: (module variables: ccpp error handling)

      endif    ! end if_lavoflg_block


! =================
      contains
! =================

!> This subroutine writes aerosol parameter configuration to run log file.
!--------------------------------
      subroutine wrt_aerlog(iaermdl, iaerflg, lalw1bd, errflg, errmsg)
!  ==================================================================  !
!                                                                      !
!  subprogram : wrt_aerlog                                             !
!                                                                      !
!    write aerosol parameter configuration to run log file.            !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
!  internal module variables:                                          !
!   lalwflg  - toposphere lw aerosol effect: =f:no; =t:yes             !
!   laswflg  - toposphere sw aerosol effect: =f:no; =t:yes             !
!   lavoflg  - stratosphere volcanic aeros effect: =f:no; =t:yes       !
!                                                                      !
!  inputs:                                                             !
!   iaerflg  - aerosol effect control flag: 3-digits (volc,lw,sw)      !
!   iaermdl  - tropospheric aerosol model scheme flag                  !
!                                                                      !
!  outputs:                                                            !
!   errmsg   - CCPP error message                                      !
!   errflg   - CCPP error flag                                         !
!                                                                      !
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call wrt_aerlog                                           !
!                                                                      !
!  ==================================================================  !

!  ---  inputs: ()
      integer,          intent(in) :: iaermdl, iaerflg
      logical,          intent(in) :: lalw1bd
!  ---  output: (CCPP error handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg
!  ---  locals:

!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      print *, VTAGAER    ! print out version tag

      if ( iaermdl==0 .or. iaermdl==5 ) then
        print *,' - Using OPAC-seasonal climatology for tropospheric',  &
     &          ' aerosol effect'
      elseif ( iaermdl == 1 ) then
        print *,' - Using GOCART-climatology for tropospheric',         &
     &          ' aerosol effect'
      elseif ( iaermdl == 2 ) then
        print *,' - Using GOCART-prognostic aerosols for tropospheric', &
     &          ' aerosol effect'
      else
        print *,' !!! ERROR in selection of aerosol model scheme',      &
     &          ' IAER_MDL =',iaermdl
        errflg = 1
        errmsg = 'ERROR(wrt_aerlog): Selected aerosol model scheme is'//&
     &       'is invalid'
        return
      endif   ! end_if_iaermdl_block

      print *,'   IAER=',iaerflg,'  LW-trop-aer=',lalwflg,              &
     &        '  SW-trop-aer=',laswflg,'  Volc-aer=',lavoflg

      if ( iaerflg <= 0 ) then        ! turn off all aerosol effects
        print *,' - No tropospheric/volcanic aerosol effect included'
        print *,'      Input values of aerosol optical properties to'   &
     &         ,' both SW and LW radiations are set to zeros'
      else
        if ( iaerflg >= 100 ) then    ! incl stratospheric volcanic aerosols
          print *,' - Include stratospheric volcanic aerosol effect'
        else                       ! no stratospheric volcanic aerosols
          print *,' - No stratospheric volcanic aerosol effect'
        endif

        if ( laswflg ) then          ! chcek for sw effect
          print *,'   - Compute multi-band aerosol optical'             &
     &           ,' properties for SW input parameters'
        else
          print *,'   - No SW radiation aerosol effect, values of'      &
     &           ,' aerosol properties to SW input are set to zeros'
        endif

        if ( lalwflg ) then          ! check for lw effect
          if ( lalw1bd ) then
            print *,'   - Compute 1 broad-band aerosol optical'         &
     &           ,' properties for LW input parameters'
          else
            print *,'   - Compute multi-band aerosol optical'           &
     &           ,' properties for LW input parameters'
          endif
        else
          print *,'   - No LW radiation aerosol effect, values of'      &
     &           ,' aerosol properties to LW input are set to zeros'
        endif
      endif     ! end if_iaerflg_block
!
      return
!................................
      end subroutine wrt_aerlog
!--------------------------------

!> This subroutine defines the one wavenumber solar fluxes based on toa
!! solar spectral distribution, and define the one wavenumber IR fluxes
!! based on black-body emission distribution at a predefined temperature.
      subroutine set_spectrum(con_pi, con_t0c, con_c, con_boltz,        &
     &     con_plnk, errflg, errmsg)

!  ==================================================================  !
!                                                                      !
!  subprogram : set_spectrum                                           !
!                                                                      !
!    define the one wavenumber solar fluxes based on toa solar spectral!
!    distrobution, and define the one wavenumber ir fluxes based on    !
!    black-body emission distribution at a predefined temperature.     !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
!> -  inputs: (CCPP Interstitials)
!!  -   con_pi:  Physical constant (pi)
!!  -   con_t0c: Physical constant (temperature kelvin at zero celcius)
!!  -   con_c:   Physical constant (speed of light)
!!  -   con_boltz: Physical constant (Boltzmann constant)
!!  -   con_plnk: Physical constant (Planck constant)
!!
!> -  inputs: (in-scope variables)
!!  -   NWVTOT:  total num of wave numbers used in sw spectrum
!!  -   NWVTIR:  total num of wave numbers used in the ir region
!!
!> -  outputs: (in-scope variables)
!!  -   solfwv(NWVTOT):   solar flux for each individual wavenumber
!!                        (\f$W/m^2\f$)
!!  -   eirfwv(NWVTIR):   ir flux(273k) for each individual wavenumber
!!                        (\f$W/m^2\f$)
!!
!> -  outputs: (CCPP error-handling) 
!!  -   errflg:           CCPP error flag
!!  -   errmsg:           CCPP error message
!                                                                      !
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call set_spectrum                                         !
!                                                                      !
!  ==================================================================  !

!  ---  inputs: (module constants)
!     integer :: NWVTOT, NWVTIR
!  ---  inputs: (CCPP Interstitials)
      real(kind_phys),intent(in) :: con_pi, con_t0c, con_c, con_boltz,  &
     &     con_plnk

!  ---  output: (in-scope variables)
!     real (kind=kind_phys), dimension(NWVTOT) :: solfwv        ! one wvn sol flux
!     real (kind=kind_phys), dimension(NWVTIR) :: eirfwv        ! one wvn ir flux
!  ---  output: (CCPP error-handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg
!  ---  locals:
      real (kind=kind_phys) :: soltot, tmp1, tmp2, tmp3

      integer :: nb, nw, nw1, nw2, nmax, nmin

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0
!
!===>  ...  begin here
!
!     nmax = min( NWVTOT, nint( maxval(wvnsw2) ))
!     nmin = max( 1,      nint( minval(wvnsw1) ))

!  ---  check print
!     print *,' MINWVN, MAXWVN = ',nmin, nmax
!  --- ...  define the one wavenumber solar fluxes based on toa solar
!           spectral distribution

!     soltot1 = f_zero
!     soltot  = f_zero
      do nb = 1, NWVSOL
        if ( nb == 1 ) then
          nw1 = 1
        else
          nw1 = nw1 + nwvns0(nb-1)
        endif

        nw2 = nw1 + nwvns0(nb) - 1

        do nw = nw1, nw2
          solfwv(nw) = s0intv(nb)
!         soltot1 = soltot1 + s0intv(nb)
!         if ( nw >= nmin .and. nw <= nmax ) then
!           soltot = soltot + s0intv(nb)
!         endif
        enddo
      enddo

!  --- ...  define the one wavenumber ir fluxes based on black-body
!           emission distribution at a predefined temperature

      tmp1 = (con_pi + con_pi) * con_plnk * con_c* con_c
      tmp2 = con_plnk * con_c / (con_boltz * con_t0c)

!$omp parallel do private(nw,tmp3)
      do nw = 1, NWVTIR
        tmp3 = 100.0 * nw
        eirfwv(nw) = (tmp1 * tmp3**3) / (exp(tmp2*tmp3) - 1.0)
      enddo
!
      return
!................................
      end subroutine set_spectrum
!--------------------------------


!> The initialization program for stratospheric volcanic aerosols.
!-----------------------------
      subroutine set_volcaer(errflg, errmsg)
!.............................
!  ---  inputs:   ( none )                                             !
!  outputs: (CCPP error handling)                                      !
!   errflg           - CCPP error flag                                 !
!   errmsg           - CCPP error message                              ! 
!  ==================================================================  !
!                                                                      !
!  subprogram : set_volcaer                                            !
!                                                                      !
!    this is the initialization progrmam for stratospheric volcanic    !
!    aerosols.                                                         !
!                                                                      !
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call set_volcaer                                          !
!                                                                      !
!  ==================================================================  !

!  ---  inputs: (none)

!  ---  output: (CCPP error handling)
!     integer :: ivolae(:,:,:)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg
!  ---  locals:
!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

!  ---  allocate data space

      if ( .not. allocated(ivolae) ) then
        allocate ( ivolae(12,4,10) )   ! for 12-mon,4-lat_zone,10-year
      endif
!
      return
!................................
      end subroutine set_volcaer
!--------------------------------
!
!...................................
      end subroutine aer_init
!-----------------------------------


!> This subroutine is the opac-climatology aerosol initialization
!! program to set up necessary parameters and working arrays.
!>\param solfwv     (NWVTOT), solar flux for each individual wavenumber
!!                  \f$(w/m^2)\f$
!!\param eirfwv     (NWVTIR), IR flux(273k) for each individual wavenumber
!!                  \f$(w/m^2)\f$
!!\param me         print message control flag
!!
!!\section gen_clim_aerinit General Algorithm
      subroutine clim_aerinit                                           &
     &     ( solfwv, eirfwv, me, aeros_file,                            &          ! ---  inputs
     &     errflg, errmsg)                                                         !  ---  outputs

!  ==================================================================  !
!                                                                      !
!  clim_aerinit is the opac-climatology aerosol initialization program !
!  to set up necessary parameters and working arrays.                  !
!                                                                      !
!  inputs:                                                             !
!   solfwv(NWVTOT)   - solar flux for each individual wavenumber (w/m2)!
!   eirfwv(NWVTIR)   - ir flux(273k) for each individual wavenum (w/m2)!
!   me               - print message control flag                      !
!   aeros_file       - external aerosol data file name                 !
!                                                                      !
!  outputs: (CCPP error handling)                                      !
!   errflg           - CCPP error flag                                 !
!   errmsg           - CCPP error message                              !
!                                                                      !
!  internal module variables:                                          !
!     lalwflg - logical lw aerosols effect control flag                !
!               =t compute lw aerosol optical prop                     !
!     laswflg - logical sw aerosols effect control flag                !
!               =t compute sw aerosol optical prop                     !
!                                                                      !
!  module constants:                                                   !
!     NWVSOL  - num of wvnum regions where solar flux is constant      !
!     NWVTOT  - total num of wave numbers used in sw spectrum          !
!     NWVTIR  - total num of wave numbers used in the ir region        !
!     NSWBND  - total number of sw spectral bands                      !
!     NLWBND  - total number of lw spectral bands                      !
!     NAERBND - number of bands for climatology aerosol data           !
!     NCM1    - number of rh independent aeros species                 !
!     NCM2    - number of rh dependent aeros species                   !
!                                                                      !
!  usage:    call clim_aerinit                                         !
!                                                                      !
!  subprograms called:  set_aercoef, optavg                            !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
      real (kind=kind_phys), dimension(:) :: solfwv        ! one wvn sol flux
      real (kind=kind_phys), dimension(:) :: eirfwv        ! one wvn ir flux
      integer,  intent(in) :: me
      character(len=26), intent(in) :: aeros_file
!  ---  output: (CCPP error handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      real (kind=kind_phys), dimension(NAERBND,NCM1)       ::           &
     &       rhidext0, rhidsca0, rhidssa0, rhidasy0
      real (kind=kind_phys), dimension(NAERBND,NRHLEV,NCM2)::           &
     &       rhdpext0, rhdpsca0, rhdpssa0, rhdpasy0
      real (kind=kind_phys), dimension(NAERBND)            :: straext0

      real (kind=kind_phys), dimension(NSWBND,NAERBND) :: solwaer
      real (kind=kind_phys), dimension(NSWBND)         :: solbnd
      real (kind=kind_phys), dimension(NLWBND,NAERBND) :: eirwaer
      real (kind=kind_phys), dimension(NLWBND)         :: eirbnd

      integer, dimension(NSWBND) :: nv1, nv2
      integer, dimension(NLWBND) :: nr1, nr2
!
!===>  ...  begin here
!
! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

!  --- ...  invoke tropospheric aerosol initialization

!> - call set_aercoef() to invoke tropospheric aerosol initialization.
      call set_aercoef(aeros_file, errflg, errmsg)
!  ---  inputs:   (in-scope variables, module constants)
!  ---  outputs:  (module variables)


! =================
      contains
! =================

!> The initialization program for climatological aerosols. The program
!! reads and maps the pre-tabulated aerosol optical spectral data onto
!! corresponding SW radiation spectral bands.
!!\section det_set_aercoef General Algorithm
!--------------------------------
      subroutine set_aercoef(aeros_file,errflg, errmsg)
!................................
!  ---  inputs:   (in-scope variables, module constants)
!  ---  outputs:  (CCPP error handling)

!  ==================================================================  !
!                                                                      !
!  subprogram : set_aercoef                                            !
!                                                                      !
!    this is the initialization progrmam for climatological aerosols   !
!                                                                      !
!    the program reads and maps the pre-tabulated aerosol optical      !
!    spectral data onto corresponding sw radiation spectral bands.     !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
!  inputs:  (in-scope variables, module constants)                     !
!   solfwv(:)    - real, solar flux for individual wavenumber (w/m2)   !
!   eirfwv(:)    - real, lw flux(273k) for individual wavenum (w/m2)   !
!   me           - integer, select cpu number as print control flag    !
!                                                                      !
!  outputs: (to the module variables)                                  !
!  outputs: (CCPP error handling)                                      !
!   errflg       - CCPP error flag                                     !
!   errmsg       - CCPP error message                                  !
!                                                                      !
!  external module variables:                                          !
!   lalwflg   - module control flag for lw trop-aer: =f:no; =t:yes     !
!   laswflg   - module control flag for sw trop-aer: =f:no; =t:yes     !
!   aeros_file- external aerosol data file name                        !
!                                                                      !
!  internal module variables:                                          !
!     IMXAE   - number of longitude points in global aeros data set    !
!     JMXAE   - number of latitude points in global aeros data set     !
!     wvnsw1,wvnsw2 (NSWSTR:NSWEND)                                    !
!             - start/end wavenumbers for each of sw bands             !
!     wvnlw1,wvnlw2 (     1:NBDLW)                                     !
!             - start/end wavenumbers for each of lw bands             !
!     NSWLWBD - total num of bands (sw+lw) for aeros optical properties!
!     NSWBND  - number of sw spectral bands actually invloved          !
!     NLWBND  - number of lw spectral bands actually invloved          !
!     NIAERCM - unit number for reading input data set                 !
!     extrhi  - extinction coef for rh-indep aeros         NCM1*NSWLWBD!
!     scarhi  - scattering coef for rh-indep aeros         NCM1*NSWLWBD!
!     ssarhi  - single-scat-alb for rh-indep aeros         NCM1*NSWLWBD!
!     asyrhi  - asymmetry factor for rh-indep aeros        NCM1*NSWLWBD!
!     extrhd  - extinction coef for rh-dep aeros    NRHLEV*NCM2*NSWLWBD!
!     scarhd  - scattering coef for rh-dep aeros    NRHLEV*NCM2*NSWLWBD!
!     ssarhd  - single-scat-alb for rh-dep aeros    NRHLEV*NCM2*NSWLWBD!
!     asyrhd  - asymmetry factor for rh-dep aeros   NRHLEV*NCM2*NSWLWBD!
!                                                                      !
!  major local variables:                                              !
!   for handling spectral band structures                              !
!     iendwv   - ending wvnum (cm**-1) for each band  NAERBND          !
!   for handling optical properties of rh independent species (NCM1)   !
!         1. insoluble        (inso); 2. soot             (soot);      !
!         3. mineral nuc mode (minm); 4. mineral acc mode (miam);      !
!         5. mineral coa mode (micm); 6. mineral transport(mitr).      !
!     rhidext0 - extinction coefficient             NAERBND*NCM1       !
!     rhidsca0 - scattering coefficient             NAERBND*NCM1       !
!     rhidssa0 - single scattering albedo           NAERBND*NCM1       !
!     rhidasy0 - asymmetry parameter                NAERBND*NCM1       !
!   for handling optical properties of rh ndependent species (NCM2)    !
!         1. water soluble    (waso); 2. sea salt acc mode(ssam);      !
!         3. sea salt coa mode(sscm); 4. sulfate droplets (suso).      !
!         rh level (NRHLEV): 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99%    !
!     rhdpext0 - extinction coefficient             NAERBND,NRHLEV,NCM2!
!     rhdpsca0 - scattering coefficient             NAERBND,NRHLEV,NCM2!
!     rhdpssa0 - single scattering albedo           NAERBND,NRHLEV,NCM2!
!     rhdpasy0 - asymmetry parameter                NAERBND,NRHLEV,NCM2!
!   for handling optical properties of stratospheric bkgrnd aerosols   !
!     straext0 - extingction coefficients             NAERBND          !
!                                                                      !
!  usage:    call set_aercoef                                          !
!                                                                      !
!  subprograms called:  optavg                                         !
!                                                                      !
!  ==================================================================  !
!
!  ---  inputs:  ( none )
      character(len=26),intent(in) :: aeros_file
!  ---  output: (CCPP error handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      integer, dimension(NAERBND) :: iendwv

      integer :: i, j, k, m, mb, ib, ii, id, iw, iw1, iw2, ik, ibs, ibe

      real (kind=kind_phys) :: sumsol, sumir, fac, tmp, wvs, wve

      logical :: file_exist
      character :: cline*80
!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

!> -# Reading climatological aerosols optical data from aeros_file,
!! including:

      inquire (file=aeros_file, exist=file_exist)

      if ( file_exist ) then
        close (NIAERCM)
        open  (unit=NIAERCM,file=aeros_file,status='OLD',               &
     &        action='read',form='FORMATTED')
        rewind (NIAERCM)
      else
        print *,'    Requested aerosol data file "',aeros_file,         &
     &          '" not found!'
        print *,'    *** Stopped in subroutine aero_init !!'
        errflg = 1
        errmsg = 'ERROR(set_aercoef): Requested aerosol data file '//   &
     &       aeros_file//' not found'
        return
      endif     ! end if_file_exist_block

!  --- ...  skip monthly global distribution

      do m = 1, 12
        read (NIAERCM,12) cline
  12    format(a80/)

        do j = 1, JMXAE
          do i = 1, IMXAE
            read(NIAERCM,*) id
          enddo
        enddo
      enddo   ! end do_m_block

!  --- ...  aloocate and input aerosol optical data

      if ( .not. allocated( extrhi ) ) then
        allocate ( extrhi (       NCM1,NSWLWBD) )
        allocate ( scarhi (       NCM1,NSWLWBD) )
        allocate ( ssarhi (       NCM1,NSWLWBD) )
        allocate ( asyrhi (       NCM1,NSWLWBD) )
        allocate ( extrhd (NRHLEV,NCM2,NSWLWBD) )
        allocate ( scarhd (NRHLEV,NCM2,NSWLWBD) )
        allocate ( ssarhd (NRHLEV,NCM2,NSWLWBD) )
        allocate ( asyrhd (NRHLEV,NCM2,NSWLWBD) )
        allocate ( extstra(            NSWLWBD) )
      endif

!>  - ending wave num for 61 aerosol spectral bands
      read(NIAERCM,21) cline
  21  format(a80)
      read(NIAERCM,22) iendwv(:)
  22  format(13i6)

!>  - atmos scale height for 5 domains, 7 profs
      read(NIAERCM,21) cline
      read(NIAERCM,24) haer(:,:)
  24  format(20f4.1)

!>  - reference pressure for 5 domains, 7 profs
      read(NIAERCM,21) cline
      read(NIAERCM,26) prsref(:,:)
  26  format(10f7.2)

!>  - rh independent ext coef for 61 bands, 6 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhidext0(:,:)
  28  format(8e10.3)

!>  - rh independent sca coef for 61 bands, 6 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhidsca0(:,:)

!>  - rh independent ssa coef for 61 bands, 6 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhidssa0(:,:)

!>  - rh independent asy coef for 61 bands, 6 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhidasy0(:,:)

!>  - rh dependent ext coef for 61 bands, 8 rh lev, 4 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhdpext0(:,:,:)

!>  - rh dependent sca coef for 61 bands, 8 rh lev, 4 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhdpsca0(:,:,:)

!>  - rh dependent ssa coef for 61 bands, 8 rh lev, 4 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhdpssa0(:,:,:)

!>  - rh dependent asy coef for 61 bands, 8 rh lev, 4 species
      read(NIAERCM,21) cline
      read(NIAERCM,28) rhdpasy0(:,:,:)

!>  - stratospheric background aeros for 61 bands
      read(NIAERCM,21) cline
      read(NIAERCM,28) straext0(:)

      close (NIAERCM)

!> -# Convert pressure reference level (in mb) to sigma reference level
!!    assume an 1000mb reference surface pressure.

      sigref(:,:) = 0.001 * prsref(:,:)

!> -# Compute solar flux weights and interval indices for mapping
!!    spectral bands between SW radiation and aerosol data.

      if ( laswflg ) then
        solbnd (:)   = f_zero
!$omp parallel do private(i,j)
        do j=1,naerbnd
          do i=1,nswbnd
            solwaer(i,j) = f_zero
          enddo
        enddo

        ibs = 1
        ibe = 1
        wvs = wvn_sw1(1)
        wve = wvn_sw1(1)
        nv_aod = 1
        do ib = 2, NSWBND
          mb = ib + NSWSTR - 1
          if ( wvn_sw2(mb) >= wvn550 .and. wvn550 >= wvn_sw1(mb) ) then
            nv_aod = ib                  ! sw band number covering 550nm wavelenth
          endif

          if ( wvn_sw1(mb) < wvs ) then
            wvs = wvn_sw1(mb)
            ibs = ib
          endif
          if ( wvn_sw1(mb) > wve ) then
            wve = wvn_sw1(mb)
            ibe = ib
          endif
        enddo

!$omp parallel do private(ib,mb,ii,iw1,iw2,iw,sumsol,fac,tmp,ibs,ibe)
        do ib = 1, NSWBND
          mb = ib + NSWSTR - 1
          ii = 1
          iw1 = nint(wvn_sw1(mb))
          iw2 = nint(wvn_sw2(mb))

          Lab_swdowhile : do while ( iw1 > iendwv(ii) )
            if ( ii == NAERBND ) exit Lab_swdowhile
            ii = ii + 1
          enddo  Lab_swdowhile

          if ( lmap_new ) then
            if (ib == ibs) then
              sumsol = f_zero
            else
              sumsol = -0.5 * solfwv(iw1)
            endif
            if (ib == ibe) then
              fac = f_zero
            else
              fac = -0.5
            endif
            solbnd(ib) = sumsol
          else
            sumsol = f_zero
          endif
          nv1(ib) = ii

          do iw = iw1, iw2
            solbnd(ib) = solbnd(ib) + solfwv(iw)
            sumsol     = sumsol     + solfwv(iw)

            if ( iw == iendwv(ii) ) then
              solwaer(ib,ii) = sumsol

              if ( ii < NAERBND ) then
                sumsol = f_zero
                ii = ii + 1
              endif
            endif
          enddo

          if ( iw2 /= iendwv(ii) ) then
            solwaer(ib,ii) = sumsol
          endif

          if ( lmap_new ) then
            tmp = fac * solfwv(iw2)
            solwaer(ib,ii) = solwaer(ib,ii) + tmp
            solbnd(ib) = solbnd(ib) + tmp
          endif

          nv2(ib) = ii
!         frcbnd(ib) = solbnd(ib) / soltot
        enddo     ! end do_ib_block for sw
      endif    ! end if_laswflg_block

!> -# Compute LW flux weights and interval indices for mapping
!!    spectral bands between lw radiation and aerosol data.

      if ( lalwflg ) then
        eirbnd (:)   = f_zero
!$omp parallel do private(i,j)
        do j=1,naerbnd
          do i=1,nlwbnd
            eirwaer(i,j) = f_zero
          enddo
        enddo

        ibs = 1
        ibe = 1
        if (NLWBND > 1 ) then
          wvs = wvn_lw1(1)
          wve = wvn_lw1(1)
          do ib = 2, NLWBND
            mb = ib + NLWSTR - 1
            if ( wvn_lw1(mb) < wvs ) then
              wvs = wvn_lw1(mb)
              ibs = ib
            endif
            if ( wvn_lw1(mb) > wve ) then
              wve = wvn_lw1(mb)
              ibe = ib
            endif
          enddo
        endif

!$omp parallel do private(ib,ii,iw1,iw2,iw,mb,sumir,fac,tmp,ibs,ibe)
        do ib = 1, NLWBND
          ii = 1
          if ( NLWBND == 1 ) then
!           iw1 = 250                   ! corresponding 40 mu
            iw1 = 400                   ! corresponding 25 mu
            iw2 = 2500                  ! corresponding 4  mu
          else
            mb = ib + NLWSTR - 1
            iw1 = nint(wvn_lw1(mb))
            iw2 = nint(wvn_lw2(mb))
          endif

          Lab_lwdowhile : do while ( iw1 > iendwv(ii) )
            if ( ii == NAERBND ) exit Lab_lwdowhile
            ii = ii + 1
          enddo  Lab_lwdowhile

          if ( lmap_new ) then
            if (ib == ibs) then
              sumir = f_zero
            else
              sumir = -0.5 * eirfwv(iw1)
            endif
            if (ib == ibe) then
              fac = f_zero
            else
              fac = -0.5
            endif
            eirbnd(ib) = sumir
          else
            sumir = f_zero
          endif
          nr1(ib) = ii

          do iw = iw1, iw2
            eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
            sumir  = sumir  + eirfwv(iw)

            if ( iw == iendwv(ii) ) then
              eirwaer(ib,ii) = sumir

              if ( ii < NAERBND ) then
                sumir = f_zero
                ii = ii + 1
              endif
            endif
          enddo

          if ( iw2 /= iendwv(ii) ) then
            eirwaer(ib,ii) = sumir
          endif

          if ( lmap_new ) then
            tmp = fac * eirfwv(iw2)
            eirwaer(ib,ii) = eirwaer(ib,ii) + tmp
            eirbnd(ib) = eirbnd(ib) + tmp
          endif

          nr2(ib) = ii
        enddo     ! end do_ib_block for lw
      endif    ! end if_lalwflg_block

!> -# Call optavg() to compute spectral band mean properties for each
!! species.

      call optavg
!  ---  inputs:  (in-scope variables, module variables)
!  ---  outputs: (module variables)

!  ---  check print
!     do ib = 1, NSWBND
!       print *,' After optavg, for sw band:',ib
!       print *,'  extrhi:', extrhi(:,ib)
!       print *,'  scarhi:', scarhi(:,ib)
!       print *,'  ssarhi:', ssarhi(:,ib)
!       print *,'  asyrhi:', asyrhi(:,ib)
!       mb = ib + NSWSTR - 1
!       print *,'  wvnsw1,wvnsw2 :',wvnsw1(mb),wvnsw2(mb)
!       do i = 1, NRHLEV
!         print *,'  extrhd for rhlev:',i
!         print *,extrhd(i,:,ib)
!         print *,'  scarhd for rhlev:',i
!         print *,scarhd(i,:,ib)
!         print *,'  ssarhd for rhlev:',i
!         print *,ssarhd(i,:,ib)
!         print *,'  asyrhd for rhlev:',i
!         print *,asyrhd(i,:,ib)
!       enddo
!       print *,' extstra:', extstra(ib)
!     enddo
!     print *,'  wvnlw1 :',wvnlw1
!     print *,'  wvnlw2 :',wvnlw2
!     do ib = 1, NLWBND
!       ii = NSWBND + ib
!       print *,' After optavg, for lw band:',ib
!       print *,'  extrhi:', extrhi(:,ii)
!       print *,'  scarhi:', scarhi(:,ii)
!       print *,'  ssarhi:', ssarhi(:,ii)
!       print *,'  asyrhi:', asyrhi(:,ii)
!       do i = 1, NRHLEV
!         print *,'  extrhd for rhlev:',i
!         print *,extrhd(i,:,ii)
!         print *,'  scarhd for rhlev:',i
!         print *,scarhd(i,:,ii)
!         print *,'  ssarhd for rhlev:',i
!         print *,ssarhd(i,:,ii)
!         print *,'  asyrhd for rhlev:',i
!         print *,asyrhd(i,:,ii)
!       enddo
!       print *,' extstra:', extstra(ii)
!     enddo
!
      return
!................................
      end subroutine set_aercoef
!--------------------------------

!> This subroutine computes mean aerosols optical properties over each
!! SW radiation spectral band for each of the species components. This
!! program follows GFDL's approach for thick cloud optical property in
!! SW radiation scheme (2000).
!--------------------------------
      subroutine optavg
!................................
!  ---  inputs:  (in-scope variables, module variables
!  ---  outputs: (module variables)

! ==================================================================== !
!                                                                      !
! subprogram: optavg                                                   !
!                                                                      !
!   compute mean aerosols optical properties over each sw radiation    !
!   spectral band for each of the species components.  This program    !
!   follows gfdl's approach for thick cloud opertical property in      !
!   sw radiation scheme (2000).                                        !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
! major input variables:                                               !
!   nv1,nv2 (NSWBND) - start/end spectral band indices of aerosol data !
!                      for each sw radiation spectral band             !
!   nr1,nr2 (NLWBND) - start/end spectral band indices of aerosol data !
!                      for each ir radiation spectral band             !
!   solwaer (NSWBND,NAERBND)                                           !
!                    - solar flux weight over each sw radiation band   !
!                      vs each aerosol data spectral band              !
!   eirwaer (NLWBND,NAERBND)                                           !
!                    - ir flux weight over each lw radiation band      !
!                      vs each aerosol data spectral band              !
!   solbnd  (NSWBND) - solar flux weight over each sw radiation band   !
!   eirbnd  (NLWBND) - ir flux weight over each lw radiation band      !
!   NSWBND           - total number of sw spectral bands               !
!   NLWBND           - total number of lw spectral bands               !
!                                                                      !
! external module variables:                                           !
!   laswflg          - control flag for sw spectral region             !
!   lalwflg          - control flag for lw spectral region             !
!                                                                      !
! output variables: (to module variables)                              !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
!  ---  output:

!  ---  locals:
      real (kind=kind_phys) :: sumk, sums, sumok, sumokg, sumreft,      &
     &       sp, refb, reft, rsolbd, rirbd

      integer :: ib, nb, ni, nh, nc
!
!===> ...  begin here
!
!  --- ...  loop for each sw radiation spectral band

      if ( laswflg ) then

!$omp parallel do private(nb,nc,sumk,sums,sumok,sumokg,sumreft)
!$omp+private(ni,nh,sp,reft,refb,rsolbd)
        do nb = 1, NSWBND
          rsolbd = f_one / solbnd(nb)

!  ---  for rh independent aerosol species

          do nc = 1, NCM1        !  ---  for rh independent aerosol species
            sumk    = f_zero
            sums    = f_zero
            sumok   = f_zero
            sumokg  = f_zero
            sumreft = f_zero

            do ni = nv1(nb), nv2(nb)
              sp   = sqrt( (f_one - rhidssa0(ni,nc))                    &
     &             / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
              reft = (f_one - sp) / (f_one + sp)
              sumreft = sumreft + reft*solwaer(nb,ni)

              sumk    = sumk    + rhidext0(ni,nc)*solwaer(nb,ni)
              sums    = sums    + rhidsca0(ni,nc)*solwaer(nb,ni)
              sumok   = sumok   + rhidssa0(ni,nc)*solwaer(nb,ni)        &
     &                * rhidext0(ni,nc)
              sumokg  = sumokg  + rhidssa0(ni,nc)*solwaer(nb,ni)        &
     &                * rhidext0(ni,nc)*rhidasy0(ni,nc)
            enddo

            refb = sumreft * rsolbd

            extrhi(nc,nb) = sumk   * rsolbd
            scarhi(nc,nb) = sums   * rsolbd
            asyrhi(nc,nb) = sumokg / (sumok + 1.0e-10)
            ssarhi(nc,nb) = 4.0*refb                                    &
     &         / ( (f_one+refb)**2 - asyrhi(nc,nb)*(f_one-refb)**2 )
          enddo   ! end do_nc_block for rh-ind aeros


          do nc = 1, NCM2        !  ---  for rh dependent aerosols species
            do nh = 1, NRHLEV
              sumk    = f_zero
              sums    = f_zero
              sumok   = f_zero
              sumokg  = f_zero
              sumreft = f_zero

              do ni = nv1(nb), nv2(nb)
                sp   = sqrt( (f_one - rhdpssa0(ni,nh,nc))               &
     &               / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
                reft = (f_one - sp) / (f_one + sp)
                sumreft = sumreft + reft*solwaer(nb,ni)

                sumk    = sumk    + rhdpext0(ni,nh,nc)*solwaer(nb,ni)
                sums    = sums    + rhdpsca0(ni,nh,nc)*solwaer(nb,ni)
                sumok   = sumok   + rhdpssa0(ni,nh,nc)*solwaer(nb,ni)   &
     &                  * rhdpext0(ni,nh,nc)
                sumokg  = sumokg  + rhdpssa0(ni,nh,nc)*solwaer(nb,ni)   &
     &                  * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
              enddo

              refb = sumreft * rsolbd

              extrhd(nh,nc,nb) = sumk   * rsolbd
              scarhd(nh,nc,nb) = sums   * rsolbd
              asyrhd(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
              ssarhd(nh,nc,nb) = 4.0*refb                               &
     &         / ( (f_one+refb)**2 - asyrhd(nh,nc,nb)*(f_one-refb)**2 )
            enddo   ! end do_nh_block
          enddo   ! end do_nc_block for rh-dep aeros

!  ---  for stratospheric background aerosols

          sumk = f_zero
          do ni = nv1(nb), nv2(nb)
            sumk = sumk + straext0(ni)*solwaer(nb,ni)
          enddo

          extstra(nb) = sumk * rsolbd

!  ---  check print
!         if ( nb > 6 .and. nb < 10) then
!           print *,' in optavg for sw band',nb
!           print *,'  nv1, nv2:',nv1(nb),nv2(nb)
!           print *,'  solwaer:',solwaer(nb,nv1(nb):nv2(nb))
!           print *,'  extrhi:', extrhi(:,nb)
!           do i = 1, NRHLEV
!             print *,'  extrhd for rhlev:',i
!             print *,extrhd(i,:,nb)
!           enddo
!           print *,'  sumk, rsolbd, extstra:',sumk,rsolbd,extstra(nb)
!         endif

        enddo   !  end do_nb_block for sw
      endif   !  end if_laswflg_block

!  --- ...  loop for each lw radiation spectral band

      if ( lalwflg ) then

!$omp parallel do private(nb,ib,nc,rirbd,sumk,sums,sumok,sumokg,sumreft)
!$omp+private(ni,nh,sp,reft,refb,rsolbd)
        do nb = 1, NLWBND

          ib = NSWBND + nb
          rirbd = f_one / eirbnd(nb)

          do nc = 1, NCM1        !  ---  for rh independent aerosol species
            sumk    = f_zero
            sums    = f_zero
            sumok   = f_zero
            sumokg  = f_zero
            sumreft = f_zero

            do ni = nr1(nb), nr2(nb)
              sp   = sqrt( (f_one - rhidssa0(ni,nc))                    &
     &             / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
              reft = (f_one - sp) / (f_one + sp)
              sumreft = sumreft + reft*eirwaer(nb,ni)

              sumk    = sumk    + rhidext0(ni,nc)*eirwaer(nb,ni)
              sums    = sums    + rhidsca0(ni,nc)*eirwaer(nb,ni)
              sumok   = sumok   + rhidssa0(ni,nc)*eirwaer(nb,ni)        &
     &                * rhidext0(ni,nc)
              sumokg  = sumokg  + rhidssa0(ni,nc)*eirwaer(nb,ni)        &
     &                * rhidext0(ni,nc)*rhidasy0(ni,nc)
            enddo

            refb = sumreft * rirbd

            extrhi(nc,ib) = sumk   * rirbd
            scarhi(nc,ib) = sums   * rirbd
            asyrhi(nc,ib) = sumokg / (sumok + 1.0e-10)
            ssarhi(nc,ib) = 4.0*refb                                       &
     &         / ( (f_one+refb)**2 - asyrhi(nc,ib)*(f_one-refb)**2 )
          enddo   ! end do_nc_block for rh-ind aeros

          do nc = 1, NCM2        !  ---  for rh dependent aerosols species
            do nh = 1, NRHLEV
              sumk    = f_zero
              sums    = f_zero
              sumok   = f_zero
              sumokg  = f_zero
              sumreft = f_zero

              do ni = nr1(nb), nr2(nb)
                sp   = sqrt( (f_one - rhdpssa0(ni,nh,nc))               &
     &             / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
                reft = (f_one - sp) / (f_one + sp)
                sumreft = sumreft + reft*eirwaer(nb,ni)

                sumk    = sumk    + rhdpext0(ni,nh,nc)*eirwaer(nb,ni)
                sums    = sums    + rhdpsca0(ni,nh,nc)*eirwaer(nb,ni)
                sumok   = sumok   + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni)   &
     &                  * rhdpext0(ni,nh,nc)
                sumokg  = sumokg  + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni)   &
     &                  * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
              enddo

              refb = sumreft * rirbd

              extrhd(nh,nc,ib) = sumk   * rirbd
              scarhd(nh,nc,ib) = sums   * rirbd
              asyrhd(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
              ssarhd(nh,nc,ib) = 4.0*refb                               &
     &         / ( (f_one+refb)**2 - asyrhd(nh,nc,ib)*(f_one-refb)**2 )
            enddo   ! end do_nh_block
          enddo   ! end do_nc_block for rh-dep aeros

!  ---  for stratospheric background aerosols

          sumk = f_zero
          do ni = nr1(nb), nr2(nb)
            sumk = sumk + straext0(ni)*eirwaer(nb,ni)
          enddo

          extstra(ib) = sumk * rirbd

!  ---  check print
!         if ( nb >= 1 .and. nb < 5) then
!           print *,' in optavg for ir band:',nb
!           print *,'  nr1, nr2:',nr1(nb),nr2(nb)
!           print *,'  eirwaer:',eirwaer(nb,nr1(nb):nr2(nb))
!           print *,'  extrhi:', extrhi(:,ib)
!           do i = 1, NRHLEV
!             print *,'  extrhd for rhlev:',i
!             print *,extrhd(i,:,ib)
!           enddo
!           print *,'  sumk, rirbd, extstra:',sumk,rirbd,extstra(ib)
!         endif

        enddo   !  end do_nb_block for lw
      endif   !  end if_lalwflg_block
!
      return
!................................
      end subroutine optavg
!--------------------------------
!
!...................................
      end subroutine clim_aerinit
!-----------------------------------


!> This subroutine checks and updates time varying climatology aerosol
!! data sets.
!!
!>\param iyear    4-digit calender year
!!\param imon     month of the year
!!\param me       print message control flag
!>\section gen_aer_upd General Algorithm
!-----------------------------------
      subroutine aer_update                                             &
     &     ( iyear, imon, me, iaermdl, aeros_file, errflg, errmsg )

!  ==================================================================  !
!                                                                      !
!  aer_update checks and update time varying climatology aerosol       !
!    data sets.                                                        !
!                                                                      !
!  inputs:                                                  size       !
!     iyear      - 4-digit calender year                      1        !
!     imon       - month of the year                          1        !
!     me         - print message control flag                 1        !
!     iaermdl    - tropospheric aerosol model scheme flag     1        !
!     aeros_file - external aerosol data file name          len=26     !
!                                                                      !
!  outputs: (CCPP error handling)                           len=*      !
!     errmsg  - CCPP error message                            1        !
!     errflg  - CCPP error flag                                        !  
!                                                                      !
!  internal module variables:                                          !
!     lalwflg     - control flag for tropospheric lw aerosol           !
!     laswflg     - control flag for tropospheric sw aerosol           !
!     lavoflg     - control flag for stratospheric volcanic aerosol    !
!                                                                      !
!  usage:    call aero_update                                          !
!                                                                      !
!  subprograms called:  trop_update, volc_update                       !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
      integer,  intent(in) :: iyear, imon, me, iaermdl
      character(len=26),intent(in) :: aeros_file
!  ---  output: (CCPP error-handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg
!  ---  locals: ( none )
!
!===> ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      if ( imon < 1 .or. imon > 12 ) then
        print *,' ***** ERROR in specifying requested month !!! ',      &
     &          'imon=', imon
        print *,' ***** STOPPED in subroutinte aer_update !!!'
        errflg = 1
        errmsg = 'ERROR(aer_update): Requested month not valid'
        return
      endif

!> -# Call trop_update() to update monthly tropospheric aerosol data.
      if ( lalwflg .or. laswflg ) then

        if ( iaermdl == 0 .or. iaermdl==5 ) then    ! opac-climatology scheme
        call trop_update(aeros_file, errflg, errmsg)
        endif

      endif

!> -# Call volc_update() to update yearly stratospheric volcanic aerosol data.
      if ( lavoflg ) then
        call volc_update(errflg, errmsg)
      endif


! =================
      contains
! =================

!> This subroutine updates the monthly global distribution of aerosol
!! profiles in five degree horizontal resolution.
!--------------------------------
      subroutine trop_update(aeros_file, errflg, errmsg)

!  ==================================================================  !
!                                                                      !
!  subprogram : trop_update                                            !
!                                                                      !
!    updates the  monthly global distribution of aerosol profiles in   !
!    five degree horizontal resolution.                                !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
!  inputs:  (in-scope variables, module constants)                     !
!   imon     - integer, month of the year                              !
!   me       - integer, print message control flag                     !
!  inputs:  (CCPP Interstitials)                                       !
!   aeros_file   - external aerosol data file name                     !
!                                                                      !
!  outputs: (module variables)                                         !
! 
!  outputs: (CCPP error-handling)                                      !
!   errmsg  - Error message                                            !
!   errflg  - Error flag                                               !
!                                                                      !
!  internal module variables:                                          !
!    kprfg (    IMXAE*JMXAE)   - aeros profile index                   !
!    idxcg (NXC*IMXAE*JMXAE)   - aeros component index                 !
!    cmixg (NXC*IMXAE*JMXAE)   - aeros component mixing ratio          !
!    denng ( 2 *IMXAE*JMXAE)   - aerosols number density               !
!                                                                      !
!    NIAERCM      - unit number for input data set                     !
!                                                                      !
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call trop_update                                          !
!                                                                      !
!  ==================================================================  !

!  ---  inputs: (CCPP Interstitials)
      character(len=26),intent(in) :: aeros_file
!  ---  output: (CCPP error handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
!     real (kind=kind_io8)  :: cmix(NXC), denn, tem
      real (kind=kind_phys) :: cmix(NXC), denn, tem
      integer               :: idxc(NXC), kprf

      integer :: i, id, j, k, m, nc
      logical :: file_exist

      character :: cline*80, ctyp*3
!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

!  --- ...  reading climatological aerosols data

      inquire (file=aeros_file, exist=file_exist)

      if ( file_exist ) then
        close(NIAERCM)
        open (unit=NIAERCM,file=aeros_file,status='OLD',                &
     &        action='read',form='FORMATTED')
        rewind (NIAERCM)

        if ( me == 0 ) then
          print *,'   Opened aerosol data file: ',aeros_file
        endif
      else
        print *,'    Requested aerosol data file "',aeros_file,         &
     &          '" not found!'
        print *,'    *** Stopped in subroutine trop_update !!'
        errflg = 1
        errmsg = 'ERROR(trop_update):Requested aerosol data file '//    &
     &       aeros_file // ' not found.'
        return
      endif      ! end if_file_exist_block

!$omp parallel do private(i,j,m)
      do j = 1, JMXAE
        do i = 1, IMXAE
          do m = 1, NXC
            idxcg(m,i,j) = 0
            cmixg(m,i,j) = f_zero
          enddo
        enddo
      enddo

!$omp parallel do private(i,j)
      do j = 1, JMXAE
        do i = 1, IMXAE
          denng(1,i,j) = f_zero
          denng(2,i,j) = f_zero
        enddo
      enddo

!  --- ...  loop over 12 month global distribution

      Lab_do_12mon : do m = 1, 12

        read(NIAERCM,12) cline
  12    format(a80/)

        if ( m /= imon ) then
!         if ( me == 0 ) print *,'  *** Skipped ',cline

          do j = 1, JMXAE
            do i = 1, IMXAE
              read(NIAERCM,*) id
            enddo
          enddo
        else
          if ( me == 0 ) print *,'  --- Reading ',cline

          do j = 1, JMXAE
            do i = 1, IMXAE
              read(NIAERCM,14) (idxc(k),cmix(k),k=1,NXC),kprf,denn,nc,  &
     &                         ctyp
  14          format(5(i2,e11.4),i2,f8.2,i3,1x,a3)

              kprfg(i,j)     = kprf
              denng(1,i,j)   = denn       ! num density of 1st layer
              if ( kprf >= 6 ) then
                denng(2,i,j) = cmix(NXC)  ! num density of 2dn layer
              else
                denng(2,i,j) = f_zero
              endif

              tem = f_one
              do k = 1, NXC-1
                idxcg(k,i,j) = idxc(k)    ! component index
                cmixg(k,i,j) = cmix(k)    ! component mixing ratio
                tem          = tem - cmix(k)
              enddo
              idxcg(NXC,i,j) = idxc(NXC)
              cmixg(NXC,i,j) = tem        ! to make sure all add to 1.
            enddo
          enddo

          close (NIAERCM)
          exit  Lab_do_12mon
        endif     ! end if_m_block

      enddo  Lab_do_12mon

!  --  check print

!     print *,'  IDXCG :'
!     print 16,idxcg
! 16  format(40i3)
!     print *,'  CMIXG :'
!     print 17,cmixg
!     print *,'  DENNG :'
!     print 17,denng
!     print *,'  KPRFG :'
!     print 17,kprfg
! 17  format(8e16.9)
!
      return
!................................
      end subroutine trop_update
!--------------------------------


!> This subroutine searches historical volcanic data sets to find and
!! read in monthly 45-degree lat-zone band of optical depth.
!--------------------------------
      subroutine volc_update(errflg, errmsg)
!................................
!  ---  inputs:    (in scope variables, module variables)
!  ---  outputs:   (CCPP error handling)

!  ==================================================================  !
!                                                                      !
!  subprogram : volc_update                                            !
!                                                                      !
!    searches historical volcanic data sets to find and read in        !
!    monthly 45-degree lat-zone band data of optical depth.            !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
!  inputs:  (in-scope variables, module constants)                     !
!   iyear    - integer, 4-digit calender year                 1        !
!   imon     - integer, month of the year                     1        !
!   me       - integer, print message control flag            1        !
!   NIAERCM  - integer, unit number for input data set        1        !
!                                                                      !
!  outputs: (module variables)                                         !
!   ivolae   - integer, monthly, 45-deg lat-zone volc odp      12*4*10 !
!   kyrstr   - integer, starting year of data in the input file        !
!   kyrend   - integer, ending   year of data in the input file        !
!   kyrsav   - integer, the year of data in use in the input file      !
!   kmonsav  - integer, the month of data in use in the input file     !
!                                                                      !
!  outputs: (CCPP error-handling)                                      !
!   errmsg  - Error message                                            !
!   errflg  - Error flag                                               !
!                                                                      ! 
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call volc_aerinit                                         !
!                                                                      !
!  ==================================================================  !

!  ---  inputs: (in-scope variables, module constants)
!     integer :: iyear, imon, me, NIAERCM

!  ---  output: (module variables)
!     integer :: ivolae(:,:,:), kyrstr, kyrend, kyrsav, kmonsav
!  ---  output: (CCPP error-handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      integer :: i, j, k
      logical :: file_exist

      character :: cline*80, volcano_file*32
      data volcano_file / 'volcanic_aerosols_1850-1859.txt ' /
!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      kmonsav = imon

      if ( kyrstr<=iyear .and. iyear<=kyrend ) then   ! use previously input data
        kyrsav = iyear
        return
      else                                            ! need to input new data
        kyrsav = iyear
        kyrstr = iyear - mod(iyear,10)
        kyrend = kyrstr + 9

!  ---  check print
!       print *,'  kyrstr, kyrend, kyrsav, kmonsav =',                  &
!    &          kyrstr,kyrend,kyrsav,kmonsav

        if ( iyear < MINVYR .or. iyear > MAXVYR ) then
!         if ( .not. allocated(ivolae) ) then
!           allocate ( ivolae(12,4,10) )   ! for 12-mon,4-lat_zone,10-year
!         endif
          ivolae(:,:,:) = 1            ! set as lowest value
          if ( me == 0 ) then
            print *,'   Request volcanic date out of range,',           &
     &              ' optical depth set to lowest value'
          endif
        else
          write(volcano_file(19:27),60) kyrstr,kyrend
  60      format(i4.4,'-',i4.4)

          inquire (file=volcano_file, exist=file_exist)
          if ( file_exist ) then
            close(NIAERCM)
            open (unit=NIAERCM,file=volcano_file,status='OLD',          &
     &            action='read',form='FORMATTED')

            read(NIAERCM,62) cline
  62        format(a80)

!  ---  check print
            if ( me == 0 ) then
              print *,'   Opened volcanic data file: ',volcano_file
              print *, cline
            endif

            do k = 1, 10
              do j = 1, 4
                read(NIAERCM,64) (ivolae(i,j,k),i=1,12)
  64            format(12i5)
              enddo
            enddo

            close (NIAERCM)
          else
            print *,'   Requested volcanic data file "',                &
     &              volcano_file,'" not found!'
            print *,'   *** Stopped in subroutine VOLC_AERINIT !!'
            errflg = 1
            errmsg = 'ERROR(volc_update): Requested volcanic data '//   &
     &              'file '//volcano_file//' not found!'
            return
          endif   ! end if_file_exist_block

        endif   ! end if_iyear_block
      endif   ! end if_kyrstr_block

!  ---  check print
      if ( me == 0 ) then
        k = mod(kyrsav,10) + 1
        print *,' CHECK: Sample Volcanic data used for month, year:',   &
     &           imon, iyear
        print *,  ivolae(kmonsav,:,k)
      endif
!
      return
!................................
      end subroutine volc_update
!--------------------------------
!
!...................................
      end subroutine aer_update
!-----------------------------------

!> This subroutine computes aerosols optical properties.
!>\param prsi    (IMAX,NLP1), pressure at interface in mb
!!\param prsl    (IMAX,NLAY), layer mean pressure in mb
!!\param prslk   (IMAX,NLAY), exner function = \f$(p/p0)^{rocp}\f$
!!\param tvly    (IMAX,NLAY), layer virtual temperature in K
!!\param rhlay   (IMAX,NLAY), layer mean relative humidity
!!\param slmsk   (IMAX), sea/land mask (sea:0,land:1,sea-ice:2)
!!\param tracer  (IMAX,NLAY,NTRAC), aerosol tracer concentration
!!\param xlon    (IMAX), longitude of given points in radiance, ok for
!!               both 0->2pi or -pi->+pi ranges
!!\param xlat    (IMAX), latitude of given points in radiance, default
!!               to pi/2 -> -pi/2, otherwise see in-line comment
!!\param IMAX           1, horizontal dimension of arrays
!!\param NLAY,NLP1      1, vertical dimensions of arrays
!!\param lsswr,lslwr    logical flags for sw/lw radiation calls
!!\param aerosw    (IMAX,NLAY,NBDSW,NF_AESW), aeros opt properties for sw
!!\n                    (:,:,:,1): optical depth
!!\n                    (:,:,:,2): single scattering albedo
!!\n                    (:,:,:,3): asymmetry parameter
!!\param aerolw    (IMAX,NLAY,NBDLW,NF_AELW), aeros opt properties for lw
!!\n                    (:,:,:,1): optical depth
!!\n                    (:,:,:,2): single scattering albedo
!!\n                    (:,:,:,3): asymmetry parameter
!!\param aerodp    (IMAX,NSPC1), vertically integrated optical depth
!>\section general_setaer General Algorithm
!-----------------------------------
      subroutine setaer                                                 &
     &     ( prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,aerfld,xlon,xlat,  &   !  ---  inputs
     &       IMAX,NLAY,NLP1, lsswr,lslwr,iaermdl,iaerflg,top_at_1,      &
     &       con_pi,con_rd,con_g,aerosw,aerolw,                         &   !  ---  outputs
     &       aerodp, ext550, errflg, errmsg                             &
     &     )

!  ==================================================================  !
!                                                                      !
!  setaer computes aerosols optical properties                         !
!                                                                      !
!  inputs:                                                   size      !
!     prsi    - pressure at interface              mb      IMAX*NLP1   !
!     prsl    - layer mean pressure                mb      IMAX*NLAY   !
!     prslk   - exner function = (p/p0)**rocp              IMAX*NLAY   !
!     tvly    - layer virtual temperature          k       IMAX*NLAY   !
!     rhlay   - layer mean relative humidity               IMAX*NLAY   !
!     slmsk   - sea/land mask (sea:0,land:1,sea-ice:2)       IMAX      !
!     tracer  - aerosol tracer concentration           IMAX*NLAY*NTRAC !
!     aerfld - prescribed aerosol mixing rat          IMAX*NLAY*NTRCAER!
!     xlon    - longitude of given points in radiance        IMAX      !
!               ok for both 0->2pi or -pi->+pi ranges                  !
!     xlat    - latitude of given points in radiance         IMAX      !
!               default to pi/2 -> -pi/2, otherwise see in-line comment!
!     IMAX    - horizontal dimension of arrays                  1      !
!     NLAY,NLP1-vertical dimensions of arrays                   1      !
!     lsswr,lslwr                                                      !
!             - logical flags for sw/lw radiation calls         1      !
!     con_pi  - Physical constant (pi)                                 !
!     con_t0c - Physical constant (temperature kelvin at zero celcius) !
!     con_c   - Physical constant (speed of light)                     !
!     iaermdl - tropospheric aerosol model scheme flag                 !
!     iaerflg - aerosol effect control flag                            ! 
!     top_at_1 - Vertical ordering convection flag                     !
!                                                                      !
!  outputs:                                                            !
!     aerosw - aeros opt properties for sw      IMAX*NLAY*NBDSW*NF_AESW!
!               (:,:,:,1): optical depth                               !
!               (:,:,:,2): single scattering albedo                    !
!               (:,:,:,3): asymmetry parameter                         !
!     aerolw - aeros opt properties for lw      IMAX*NLAY*NBDLW*NF_AELW!
!               (:,:,:,1): optical depth                               !
!               (:,:,:,2): single scattering albedo                    !
!               (:,:,:,3): asymmetry parameter                         !
!     tau_gocart - 550nm aeros opt depth     IMAX*NLAY*MAX_NUM_GRIDCOMP!
!!    aerodp - vertically integrated optical depth         IMAX*NSPC1  !
!                                                                      !
!     errflg  - CCPP error flag                                        !
!     errmsg  - CCPP error message                                     !
!                                                                      !
!  internal module variable:                                           !
!     laswflg - tropospheric aerosol control flag for sw radiation     !
!               =f: no sw aeros calc.  =t: do sw aeros calc.           !
!     lalwflg - tropospheric aerosol control flag for lw radiation     !
!               =f: no lw aeros calc.  =t: do lw aeros calc.           !
!     lavoflg - control flag for stratospheric vocanic aerosols        !
!               =t: add volcanic aerosols to the background aerosols   !
!  internal module variable: (set by subroutine aer_init)              !
!     ivolae  - stratosphere volcanic aerosol optical depth (fac 1.e4) !
!                                                     12*4*10          !
!  usage:    call setaer                                               !
!                                                                      !
!  subprograms called:  aer_property                                   !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
      integer, intent(in) :: IMAX, NLAY, NLP1, iaermdl, iaerflg
      real (kind=kind_phys), intent(in) :: con_pi, con_rd, con_g
      real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl,  &
     &       prslk, tvly, rhlay
      real (kind=kind_phys), dimension(:),   intent(in) :: xlon, xlat,  &
     &       slmsk
      real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer
      real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld

      logical, intent(in) :: lsswr, lslwr, top_at_1


!  ---  outputs:
      real (kind=kind_phys), dimension(:,:,:,:), intent(out) ::         &
     &       aerosw, aerolw

      real (kind=kind_phys), dimension(:,:)    , intent(out) :: aerodp
      real (kind=kind_phys), dimension(:,:)    , intent(out) :: ext550
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      real (kind=kind_phys), parameter :: psrfh = 5.0    ! ref press (mb) for upper bound

      real (kind=kind_phys), dimension(IMAX) :: alon,alat,volcae,rdelp
!     real (kind=kind_phys), dimension(IMAX) :: sumodp
      real (kind=kind_phys) :: prsln(NLP1),hz(IMAX,NLP1),dz(IMAX,NLAY)
      real (kind=kind_phys) :: tmp1, tmp2, psrfl

      integer               :: kcutl(IMAX), kcuth(IMAX)
      integer               :: i, i1, j, k, m, mb, kh, kl

      logical               :: laddsw=.false.,  laersw=.false.
      logical               :: laddlw=.false.,  laerlw=.false.

!  ---  conversion constants
      real (kind=kind_phys) :: rdg
      real (kind=kind_phys) :: rovg

!===>  ...  begin here
      rdg  = 180._kind_phys / con_pi
      rovg = 0.001_kind_phys * con_rd / con_g

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      do m = 1, NF_AESW
        do j = 1, NBDSW
          do k = 1, NLAY
            do i = 1, IMAX
              aerosw(i,k,j,m) = f_zero
            enddo
          enddo
        enddo
      enddo

      do m = 1, NF_AELW
        do j = 1, NBDLW
          do k = 1, NLAY
            do i = 1, IMAX
              aerolw(i,k,j,m) = f_zero
            enddo
          enddo
        enddo
      enddo

!     sumodp = f_zero
      do i = 1, IMAX
       do k = 1, NSPC1
         aerodp(i,k) = f_zero
       enddo
      enddo
      ext550(:,:) = f_zero

      if ( .not. (lsswr .or. lslwr) ) then
        return
      endif

      if ( iaerflg <= 0 ) then
        return
      endif

      laersw = lsswr .and. laswflg
      laerlw = lslwr .and. lalwflg

!> -# Convert lat/lon from radiance to degree.

      do i = 1, IMAX
        alon(i) = xlon(i) * rdg
        if (alon(i) < f_zero) alon(i) = alon(i) + 360.0
        alat(i) = xlat(i) * rdg          ! if xlat in pi/2 -> -pi/2 range
!       alat(i) = 90.0 - xlat(i)*rdg     ! if xlat in 0 -> pi range
      enddo

!> -# Compute level height and layer thickness.

      if ( laswflg .or. lalwflg ) then

        lab_do_IMAX : do i = 1, IMAX

          lab_if_flip : if (.not. top_at_1) then       ! input from sfc to toa

            do k = 1, NLAY
              prsln(k) = log(prsi(i,k))
            enddo
            prsln(NLP1)= log(prsl(i,NLAY))

            do k = NLAY, 1, -1
              dz(i,k) = rovg * (prsln(k) - prsln(k+1)) * tvly(i,k)
            enddo
            dz(i,NLAY)  = 2.0 * dz(i,NLAY)

            hz(i,1) = f_zero
            do k = 1, NLAY
              hz(i,k+1) = hz(i,k) + dz(i,k)
            enddo

          else  lab_if_flip                         ! input from toa to sfc

            prsln(1) = log(prsl(i,1))
            do k = 2, NLP1
              prsln(k) = log(prsi(i,k))
            enddo

            do k = 1, NLAY
              dz(i,k) = rovg * (prsln(k+1) - prsln(k)) * tvly(i,k)
            enddo
            dz(i,1) = 2.0 * dz(i,1)

            hz(i,NLP1) = f_zero
            do k = NLAY, 1, -1
              hz(i,k) = hz(i,k+1) + dz(i,k)
            enddo

          endif  lab_if_flip

        enddo  lab_do_IMAX


!> -# Calculate SW aerosol optical properties for the corresponding
!!    frequency bands:
!!    - if opac aerosol climatology is used, call aer_property(): this
!!      subroutine maps the 5 degree global climatological aerosol data
!!      set onto model grids, and compute aerosol optical properties for
!!      SW and LW radiations.
!!    - if gocart aerosol scheme is used, call setgocartaer(): this
!!      subroutine computes sw + lw aerosol optical properties for gocart
!!      aerosol species (merged from fcst and clim fields).

          if ( iaermdl==0 .or. iaermdl==5 ) then  ! use opac aerosol climatology

          call aer_property                                               &
!  ---  inputs:
     &       ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,                   &
     &         alon,alat,slmsk, laersw,laerlw,                            &
     &         IMAX,NLAY,NLP1,top_at_1,                                   &
!    &         IMAX,NLAY,NLP1,NSPC1,                                      &
!  ---  outputs:
     &         aerosw,aerolw,aerodp,errflg,errmsg                         &
     &       )

!
          elseif ( iaermdl==1 .or. iaermdl==2) then ! use gocart aerosols

          call aer_property_gocart                                        &
!  ---  inputs:
     &       ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,aerfld,            &
     &         alon,alat,slmsk,laersw,laerlw,con_rd,                      &
     &         IMAX,NLAY,NLP1,                                            &
!  ---  outputs:
     &         aerosw,aerolw,aerodp,ext550,errflg,errmsg                  &
     &       )
        endif     ! end if_iaerflg_block


!  ---  check print
!       do m = 1, NBDSW
!         print *,'  ***  CHECK AEROSOLS PROPERTIES FOR SW BAND =',m,   &
!    &            ' ***'
!         do k = 1, 10
!           print *,'  LEVEL :',k
!           print *,'  TAUAER:',aerosw(:,k,m,1)
!           print *,'  SSAAER:',aerosw(:,k,m,2)
!           print *,'  ASYAER:',aerosw(:,k,m,3)
!         enddo
!       enddo
!       print *,'  ***  CHECK AEROSOLS OPTICAL DEPTH FOR 550nm REGION'
!       print *, aerodp(:,1)
!       if ( laod_out ) then
!         do m = 1, NSPC1
!           print *,'  ***  CHECK AEROSOLS OPTICAL DEPTH FOR SPECIES:', &
!    &              m
!           print *, aerodp(:,m)
!           sumodp(:) = sumodp(:) + aerodp(:,m)
!         enddo

!
!         print *,'  ***  CHECK AEROSOLS OPTICAL DEPTH FOR ALL SPECIES:'
!         print *, sumodp(:)
!       endif
!       do m = 1, NBDLW
!         print *,'  ***  CHECK AEROSOLS PROPERTIES FOR LW BAND =',m,   &
!    &            ' ***'
!         do k = 1, 10
!           print *,'  LEVEL :',k
!           print *,'  TAUAER:',aerolw(:,k,m,1)
!           print *,'  SSAAER:',aerolw(:,k,m,2)
!           print *,'  ASYAER:',aerolw(:,k,m,3)
!         enddo
!       enddo

      endif   ! end if_laswflg_or_lalwflg_block

!> -# Compute stratosphere volcanic forcing:
!!    - select data in 4 lat bands, interpolation at the boundaries
!!    - Find lower boundary of stratosphere: polar, fixed at 25000pa
!!      (250mb); tropic, fixed at 15000pa (150mb); mid-lat, interpolation
!!    - SW: add volcanic aerosol optical depth to the background value
!!    - Smoothing profile at boundary if needed
!!    - LW: add volcanic aerosol optical depth to the background value
!  ---  ...  stratosphere volcanic forcing

      if ( lavoflg ) then

        if ( iaerflg == 100 ) then
          laddsw = lsswr
          laddlw = lslwr
        else
          laddsw = lsswr .and. laswflg
          laddlw = lslwr .and. lalwflg
        endif

        i1 = mod(kyrsav, 10) + 1

!  ---  select data in 4 lat bands, interpolation at the boundaires

        do i = 1, IMAX
          if      ( alat(i) > 46.0 ) then
            volcae(i) = 1.0e-4 * ivolae(kmonsav,1,i1)
          else if ( alat(i) > 44.0 ) then
            volcae(i) = 5.0e-5                                          &
     &                * (ivolae(kmonsav,1,i1) + ivolae(kmonsav,2,i1))
          else if ( alat(i) >  1.0 ) then
            volcae(i) = 1.0e-4 * ivolae(kmonsav,2,i1)
          else if ( alat(i) > -1.0 ) then
            volcae(i) = 5.0e-5                                          &
     &                * (ivolae(kmonsav,2,i1) + ivolae(kmonsav,3,i1))
          else if ( alat(i) >-44.0 ) then
            volcae(i) = 1.0e-4 * ivolae(kmonsav,3,i1)
          else if ( alat(i) >-46.0 ) then
            volcae(i) = 5.0e-5                                          &
     &                * (ivolae(kmonsav,3,i1) + ivolae(kmonsav,4,i1))
          else
            volcae(i) = 1.0e-4 * ivolae(kmonsav,4,i1)
          endif
        enddo

        if (top_at_1) then         ! input data from toa to sfc

!  ---  find lower boundary of stratosphere

          do i = 1, IMAX

            tmp1 = abs( alat(i) )
            if ( tmp1 > 70.0 ) then          ! polar, fixed at 25000pa (250mb)
              psrfl = 250.0
            elseif ( tmp1 < 20.0 ) then      ! tropic, fixed at 15000pa (150mb)
              psrfl = 150.0
            else                             ! mid-lat, interpolation
              psrfl = 110.0 + 2.0*tmp1
            endif

            kcuth(i) = NLAY - 1
            kcutl(i) = 2
            rdelp(i) = f_one / prsi(i,2)

            lab_do_kcuth0 : do k = 2, NLAY-2
              if ( prsi(i,k) >= psrfh ) then
                kcuth(i) = k - 1
                exit lab_do_kcuth0
              endif
            enddo  lab_do_kcuth0

            lab_do_kcutl0 : do k = 2, NLAY-2
              if ( prsi(i,k) >= psrfl ) then
                kcutl(i) = k - 1
                rdelp(i) = f_one / (prsi(i,k) - prsi(i,kcuth(i)))
                exit lab_do_kcutl0
              endif
            enddo  lab_do_kcutl0
          enddo

!  ---  sw: add volcanic aerosol optical depth to the background value

          if ( laddsw ) then
            do m = 1, NBDSW
              mb = NSWSTR + m - 1

              if     ( wvn_sw1(mb) > 20000 ) then  ! range of wvlth < 0.5mu
                tmp2 = 0.74
              elseif ( wvn_sw2(mb) < 20000 ) then  ! range of wvlth > 0.5mu
                tmp2 = 1.14
              else                                 ! range of wvlth in btwn
                tmp2 = 0.94
              endif
              tmp1 = (0.275e-4 * (wvn_sw2(mb)+wvn_sw1(mb))) ** tmp2

              do i = 1, IMAX
                kh = kcuth(i)
                kl = kcutl(i)
                do k = kh, kl
                  tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) * rdelp(i))
                  aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
                enddo

!  ---  smoothing profile at boundary if needed

                if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl+1,m,1) ) then
                  tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl+1,m,1)
                  aerosw(i,kl  ,m,1) = 0.8 * tmp2
                  aerosw(i,kl+1,m,1) = 0.2 * tmp2
                endif
              enddo    ! end do_i_block
            enddo      ! end do_m_block

!  ---  check print
!           do i = 1, IMAX
!             print *,' LEV  PRESS      TAU      FOR PROFILE:',i,       &
!    &                '  KCUTH, KCUTL =',kcuth(i),kcutl(i)
!             kh = kcuth(i) - 1
!             kl = kcutl(i) + 10
!             do k = kh, kl
!               write(6,71) k, prsl(i,k), aerosw(i,k,1,1)
! 71            format(i3,2e11.4)
!             enddo
!           enddo

          endif        ! end if_laddsw_block

!  ---  lw: add volcanic aerosol optical depth to the background value

          if ( laddlw ) then
            if ( NLWBND == 1 ) then

              tmp1 = (0.55 / 11.0) ** 1.2
              do i = 1, IMAX
                kh = kcuth(i)
                kl = kcutl(i)
                do k = kh, kl
                  tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) * rdelp(i))  &
     &                 * volcae(i)
                  do m = 1, NBDLW
                    aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2
                  enddo
                enddo
              enddo    ! end do_i_block

            else

              do m = 1, NBDLW
                tmp1 = (0.275e-4 * (wvn_lw2(m) + wvn_lw1(m))) ** 1.2

                do i = 1, IMAX
                  kh = kcuth(i)
                  kl = kcutl(i)
                  do k = kh, kl
                    tmp2 = tmp1 * ((prsi(i,k+1)-prsi(i,k)) * rdelp(i))
                    aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
                  enddo
                enddo    ! end do_i_block
              enddo      ! end do_m_block

            endif      ! end if_NLWBND_block
          endif        ! end if_laddlw_block

        else                            ! input data from sfc to toa

!  ---  find lower boundary of stratosphere

          do i = 1, IMAX

            tmp1 = abs( alat(i) )
            if ( tmp1 > 70.0 ) then          ! polar, fixed at 25000pa (250mb)
              psrfl = 250.0
            elseif ( tmp1 < 20.0 ) then      ! tropic, fixed at 15000pa (150mb)
              psrfl = 150.0
            else                             ! mid-lat, interpolation
              psrfl = 110.0 + 2.0*tmp1
            endif

            kcuth(i) = 2
            kcutl(i) = NLAY - 1
            rdelp(i) = f_one / prsi(i,NLAY-1)

            lab_do_kcuth1 : do k = NLAY-1, 2, -1
              if ( prsi(i,k) >= psrfh ) then
                kcuth(i) = k
                exit lab_do_kcuth1
              endif
            enddo  lab_do_kcuth1

            lab_do_kcutl1 : do k = NLAY, 2, -1
              if ( prsi(i,k) >= psrfl ) then
                kcutl(i) = k
                rdelp(i) = f_one / (prsi(i,k) - prsi(i,kcuth(i)+1))
                exit lab_do_kcutl1
              endif
            enddo  lab_do_kcutl1
          enddo

!  ---  sw: add volcanic aerosol optical depth to the background value

          if ( laddsw ) then
            do m = 1, NBDSW
              mb = NSWSTR + m - 1

              if     ( wvn_sw1(mb) > 20000 ) then  ! range of wvlth < 0.5mu
                tmp2 = 0.74
              elseif ( wvn_sw2(mb) < 20000 ) then  ! range of wvlth > 0.5mu
                tmp2 = 1.14
              else                                 ! range of wvlth in btwn
                tmp2 = 0.94
              endif
              tmp1 = (0.275e-4 * (wvn_sw2(mb)+wvn_sw1(mb))) ** tmp2

              do i = 1, IMAX
                kh = kcuth(i)
                kl = kcutl(i)
                do k = kl, kh
                  tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) * rdelp(i))
                  aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
                enddo

!  ---  smoothing profile at boundary if needed

                if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl-1,m,1) ) then
                  tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl-1,m,1)
                  aerosw(i,kl  ,m,1) = 0.8 * tmp2
                  aerosw(i,kl-1,m,1) = 0.2 * tmp2
                endif
              enddo    ! end do_i_block
            enddo      ! end do_m_block

!  ---  check print
!           do i = 1, IMAX
!             print *,' LEV  PRESS      TAU      FOR PROFILE:',i,       &
!    &                '  KCUTH, KCUTL =',kcuth(i),kcutl(i)
!             kh = kcuth(i) + 1
!             kl = kcutl(i) - 10
!             do k = kh, kl, -1
!               write(6,71) NLP1-k,prsl(i,k),aerosw(i,k,1,1)
!             enddo
!           enddo

          endif        ! end if_laddsw_block

!  ---  lw: add volcanic aerosol optical depth to the background value

          if ( laddlw ) then
            if ( NLWBND == 1 ) then

              tmp1 = (0.55 / 11.0) ** 1.2
              do i = 1, IMAX
                kh = kcuth(i)
                kl = kcutl(i)
                do k = kl, kh
                  tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) * rdelp(i))  &
     &                 * volcae(i)
                  do m = 1, NBDLW
                    aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2
                  enddo
                enddo
              enddo    ! end do_i_block

            else

              do m = 1, NBDLW
                tmp1 = (0.275e-4 * (wvn_lw2(m) + wvn_lw1(m))) ** 1.2

                do i = 1, IMAX
                  kh = kcuth(i)
                  kl = kcutl(i)
                  do k = kl, kh
                    tmp2 = tmp1 * ((prsi(i,k)-prsi(i,k+1)) * rdelp(i))
                    aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
                  enddo
                enddo    ! end do_i_block
              enddo      ! end do_m_block

            endif      ! end if_NLWBND_block
          endif        ! end if_laddlw_block

        endif                           ! end if_top_at_1_block

      endif   ! end if_lavoflg_block
!
      return
!...................................
      end subroutine setaer
!-----------------------------------

!> This subroutine maps the 5 degree global climatological aerosol data
!! set onto model grids, and compute aerosol optical properties for SW
!! and LW radiations.
!!\param prsi           (IMAX,NLP1), pressure at interface in mb
!!\param prsl           (IMAX,NLAY), layer mean pressure(not used)
!!\param prslk          (IMAX,NLAY), exner function=\f$(p/p0)^{rocp}\f$ (not used)
!!\param tvly           (IMAX,NLAY), layer virtual temperature (not used)
!!\param rhlay          (IMAX,NLAY), layer mean relative humidity
!!\param dz             (IMAX,NLAY), layer thickness in m
!!\param hz             (IMAX,NLP1), level high in m
!!\param tracer         (IMAX,NLAY,NTRAC), aer tracer concentrations (not used)
!!\param alon, alat     (IMAX), longitude and latitude of given points in degree
!!\param slmsk          (IMAX), sea/land mask (sea:0,land:1,sea-ice:2)
!!\param laersw,laerlw  logical flag for sw/lw aerosol calculations
!!\param IMAX           horizontal dimension of arrays
!!\param NLAY,NLP1      vertical dimensions of arrays
!!\param NSPC           num of species for optional aod output fields
!!\param aerosw        (IMAX,NLAY,NBDSW,NF_AESW), aeros opt properties for sw
!!\n                              (:,:,:,1): optical depth
!!\n                              (:,:,:,2): single scattering albedo
!!\n                              (:,:,:,3): asymmetry parameter
!!\param aerolw        (IMAX,NLAY,NBDLW,NF_AELW), aeros opt properties for lw
!!\n                              (:,:,:,1): optical depth
!!\n                              (:,:,:,2): single scattering albedo
!!\n                              (:,:,:,3): asymmetry parameter
!!\param aerodp        (IMAX,NSPC+1), vertically integrated aer-opt-depth
!!\section gel_aer_pro General Algorithm
!-----------------------------------
      subroutine aer_property                                           &
     &     ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,                   &     !  ---  inputs:
     &       alon,alat,slmsk, laersw,laerlw,                            &
     &       IMAX,NLAY,NLP1,top_at_1,                                   &
     &       aerosw,aerolw,aerodp,errflg,errmsg                         &     !  ---  outputs:
     &     )

!  ==================================================================  !
!                                                                      !
!  aer_property maps the 5 degree global climatological aerosol data   !
!  set onto model grids, and compute aerosol optical properties for sw !
!  and lw radiations.                                                  !
!                                                                      !
!  inputs:                                                             !
!     prsi    - pressure at interface              mb      IMAX*NLP1   !
!     prsl    - layer mean pressure         (not used)     IMAX*NLAY   !
!     prslk   - exner function=(p/p0)**rocp (not used)     IMAX*NLAY   !
!     tvly    - layer virtual temperature   (not used)     IMAX*NLAY   !
!     rhlay   - layer mean relative humidity               IMAX*NLAY   !
!     dz      - layer thickness                    m       IMAX*NLAY   !
!     hz      - level high                         m       IMAX*NLP1   !
!     tracer  - aer tracer concentrations   (not used)  IMAX*NLAY*NTRAC!
!     alon, alat                                             IMAX      !
!             - longitude and latitude of given points in degree       !
!     slmsk   - sea/land mask (sea:0,land:1,sea-ice:2)       IMAX      !
!     laersw,laerlw                                             1      !
!             - logical flag for sw/lw aerosol calculations            !
!     IMAX    - horizontal dimension of arrays                  1      !
!     NLAY,NLP1-vertical dimensions of arrays                   1      !
!!    NSPC    - num of species for optional aod output fields   1      !
!     top_at_1 - vertical ordering flag                                !
!                                                                      !
!  outputs:                                                            !
!     aerosw - aeros opt properties for sw      IMAX*NLAY*NBDSW*NF_AESW!
!               (:,:,:,1): optical depth                               !
!               (:,:,:,2): single scattering albedo                    !
!               (:,:,:,3): asymmetry parameter                         !
!     aerolw - aeros opt properties for lw      IMAX*NLAY*NBDLW*NF_AELW!
!               (:,:,:,1): optical depth                               !
!               (:,:,:,2): single scattering albedo                    !
!               (:,:,:,3): asymmetry parameter                         !
!!    aerodp - vertically integrated aer-opt-depth         IMAX*NSPC+1 !
!                                                                      !
!     errflg  - CCPP error flag                                        !
!     errmsg  - CCPP error message                                     !
!                                                                      !
!  module parameters and constants:                                    !
!     NSWBND  - total number of actual sw spectral bands computed      !
!     NLWBND  - total number of actual lw spectral bands computed      !
!     NSWLWBD - total number of sw+lw bands computed                   !
!                                                                      !
!  module variable: (set by subroutine aer_init)                       !
!     kprfg   - aerosols profile index                IMXAE*JMXAE      !
!               1:ant  2:arc  3:cnt  4:mar  5:des  6:marme 7:cntme     !
!     idxcg   - aerosols component index              NXC*IMXAE*JMXAE  !
!               1:inso    2:soot    3:minm    4:miam    5:micm         !
!               6:mitr    7:waso    8:ssam    9:sscm   10:suso         !
!     cmixg   - aerosols component mixing ratio       NXC*IMXAE*JMXAE  !
!     denng   - aerosols number density                2 *IMXAE*JMXAE  !
!               1:for domain-1   2:domain-2 (prof marme/cntme only)    !
!                                                                      !
!  usage:    call aer_property                                         !
!                                                                      !
!  subprograms called:  radclimaer                                     !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
      integer, intent(in) :: IMAX, NLAY, NLP1
!     integer, intent(in) :: IMAX, NLAY, NLP1, NSPC
      logical, intent(in) :: laersw, laerlw, top_at_1

      real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl,  &
     &       prslk, tvly, rhlay, dz, hz
      real (kind=kind_phys), dimension(:),   intent(in) :: alon, alat,  &
     &       slmsk
      real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer

!  ---  outputs:
      real (kind=kind_phys), dimension(:,:,:,:), intent(out) ::         &
     &       aerosw, aerolw
      real (kind=kind_phys), dimension(:,:)    , intent(out) :: aerodp
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      real (kind=kind_phys), dimension(NCM) :: cmix
      real (kind=kind_phys), dimension(  2) :: denn
      real (kind=kind_phys), dimension(NSPC) :: spcodp

      real (kind=kind_phys), dimension(NLAY) :: delz, rh1, dz1
      integer,               dimension(NLAY) :: idmaer

      real (kind=kind_phys), dimension(NLAY,NSWLWBD):: tauae,ssaae,asyae
!test real (kind=kind_phys), dimension(IMAX,NLAY) :: aersav

      real (kind=kind_phys) :: tmp1, tmp2, rps, dtmp, h1
      real (kind=kind_phys) :: wi, wj, w11, w12, w21, w22

      integer :: i, ii, i1, i2, i3,  j1, j2, j3,  k, m, m1,             &
     &           kp, kpa, kpi, kpj

!  ---  conversion constants
      real (kind=kind_phys), parameter :: dltg = 360.0 / float(IMXAE)
      real (kind=kind_phys), parameter :: hdlt = 0.5 * dltg
      real (kind=kind_phys), parameter :: rdlt = 1.0 / dltg

!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

!> -# Map aerosol data to model grids
!!    - Map grid in longitude direction, lon from 0 to 355 deg resolution
!!    - Map grid in latitude direction, lat from 90n to 90s in 5 deg resolution

      i1 = 1
      i2 = 2
      j1 = 1
      j2 = 2

      lab_do_IMAX : do i = 1, IMAX

!  ---  map grid in longitude direction, lon from 0 to 355 deg resolution

!       print *,' Seeking lon index for point i =',i
        i3 = 1
        lab_do_IMXAE : do while ( i3 <= IMXAE )
          tmp1 = dltg * (i3 - 1)
          dtmp = alon(i) - tmp1
!         print *,'   alon, i3, tlon, dlon =',alon(i),i3,tmp1,dtmp

          if ( dtmp > dltg ) then
            i3 = i3 + 1
            if ( i3 > IMXAE ) then
              print *,' ERROR! In setclimaer alon>360. ipt =',i,        &
     &           ',  dltg,alon,tlon,dlon =',dltg,alon(i),tmp1,dtmp
              errflg = 1
              errmsg = 'ERROR(aer_property)'
              return
            endif
          elseif ( dtmp >= f_zero ) then
            i1 = i3
            i2 = mod(i3,IMXAE) + 1
            wi = dtmp * rdlt
            if ( dtmp <= hdlt ) then
              kpi = i3
            else
              kpi = i2
            endif
!           print *,'   found i1, i2, wi =',i1,i2,wi
            exit lab_do_IMXAE
          else
            i3 = i3 - 1
            if ( i3 < 1 ) then
              print *,' ERROR! In setclimaer alon< 0. ipt =',i,         &
     &           ',  dltg,alon,tlon,dlon =',dltg,alon(i),tmp1,dtmp
              errflg = 1
              errmsg = 'ERROR(aer_property)'
              return
            endif
          endif
        enddo  lab_do_IMXAE

!  ---  map grid in latitude direction, lat from 90n to 90s in 5 deg resolution

!       print *,' Seeking lat index for point i =',i
        j3 = 1
        lab_do_JMXAE : do while ( j3 <= JMXAE )
          tmp2 = 90.0 - dltg * (j3 - 1)
          dtmp = tmp2 - alat(i)
!         print *,'   alat, j3, tlat, dlat =',alat(i),j3,tmp2,dtmp

          if ( dtmp > dltg ) then
            j3 = j3 + 1
            if ( j3 >= JMXAE ) then
              print *,' ERROR! In setclimaer alat<-90. ipt =',i,        &
     &           ',  dltg,alat,tlat,dlat =',dltg,alat(i),tmp2,dtmp
              errflg = 1
              errmsg = 'ERROR(aer_property)'
              return
            endif
          elseif ( dtmp >= f_zero ) then
            j1 = j3
            j2 = j3 + 1
            wj = dtmp * rdlt
            if ( dtmp <= hdlt ) then
              kpj = j3
            else
              kpj = j2
            endif
!           print *,'   found j1, j2, wj =',j1,j2,wj
            exit lab_do_JMXAE
          else
            j3 = j3 - 1
            if ( j3 < 1 ) then
              print *,' ERROR! In setclimaer alat>90. ipt =',i,         &
     &           ',  dltg,alat,tlat,dlat =',dltg,alat(i),tmp2,dtmp
              errflg = 1
              errmsg = 'ERROR(aer_property)'
              return
            endif
          endif
        enddo  lab_do_JMXAE

!> -# Determin the type of aerosol profile (kp) and scale hight for
!!    domain 1 (h1) to be used at this grid point.

        kp = kprfg(kpi,kpj)                     ! nearest typical aeros profile as default
        kpa = max( kprfg(i1,j1),kprfg(i1,j2),kprfg(i2,j1),kprfg(i2,j2) )
        h1 = haer(1,kp)
        denn(2) = f_zero
        ii = 1

        if ( kp /= kpa ) then
          if ( kpa == 6 ) then                  ! if ocean prof with mineral aeros overlay
            ii = 2                              ! need 2 types of densities
            if ( slmsk(i) > f_zero ) then       ! but actually a land/sea-ice point
              kp = 7                            ! reset prof index to land
              h1 = 0.5*(haer(1,6) + haer(1,7))  ! use a transition scale hight
            else
              kp = kpa
              h1 = haer(1,6)
            endif
          elseif ( kpa == 7 ) then              ! if land prof with mineral aeros overlay
            ii = 2                              ! need 2 types of densities
            if ( slmsk(i) <= f_zero ) then      ! but actually an ocean point
              kp = 6                            ! reset prof index to ocean
              h1 = 0.5*(haer(1,6) + haer(1,7))  ! use a transition scale hight
            else
              kp = kpa
              h1 = haer(1,7)
            endif
          else                                  ! lower atmos without mineral aeros overlay
!           h1 = 0.5*(haer(1,kp) + haer(1,kpa)) ! use a transition scale hight
            h1 = haer(1,kpa)
            kp = kpa
          endif
        endif

!> -# Compute horizontal bi-linear interpolation weights

        w11 = (f_one-wi) * (f_one-wj)
        w12 = (f_one-wi) *       wj
        w21 =        wi  * (f_one-wj)
        w22 =        wi  * wj

!  ---  check print
!       print *,'  Grid pt', i,',   alon, alat =',alon(i),alat(i),      &
!    &                       ',   tlon, tlat =',tmp1,tmp2
!       print *,'   lon grid index i1, i2 =',i1,i2,',  weight wi =',wi
!       print *,'   lat grid index j1, j2 =',j1,j2,',  weight wj =',wj
!       print *,'   bi-linear weights w11,w21,w12,w22 =',w11,w21,w12,w22
!       print *,'   kp,kpa,slmsk,h1 =',kp,m1,slmsk(i),h1

!> -# Do horizontal bi-linear interpolation on aerosol partical density
!!   (denn)

        do m = 1, ii                            ! ii=1 for domain 1; =2 for domain 2.
          denn(m) = w11*denng(m,i1,j1) + w12*denng(m,i1,j2)             &
     &            + w21*denng(m,i2,j1) + w22*denng(m,i2,j2)
        enddo  ! end_do_m_loop

!> -# Do horizontal bi-linear interpolation on mixing ratios

        cmix(:) = f_zero
        do m = 1, NXC
          ii = idxcg(m,i1,j1)
          if ( ii > 0 ) then
            cmix(ii) = cmix(ii) + w11*cmixg(m,i1,j1)
          endif
          ii = idxcg(m,i1,j2)
          if ( ii > 0 ) then
            cmix(ii) = cmix(ii) + w12*cmixg(m,i1,j2)
          endif
          ii = idxcg(m,i2,j1)
          if ( ii > 0 ) then
            cmix(ii) = cmix(ii) + w21*cmixg(m,i2,j1)
          endif
          ii = idxcg(m,i2,j2)
          if ( ii > 0 ) then
            cmix(ii) = cmix(ii) + w22*cmixg(m,i2,j2)
          endif
        enddo  ! end_do_m_loop

!  ---  check print
!       print *,'   denn =',denn(:)
!       print *,'   cmix =',cmix(:)

!> -# Prepare to setup domain index array and effective layer thickness,
!!    also convert pressure level to sigma level to follow the terrain.

        do k = 1, NLAY
          rh1(k) = rhlay(i,k)
          dz1(k) = dz   (i,k)
        enddo

        lab_if_flip : if (.not. top_at_1) then       ! input from sfc to toa

          if ( prsi(i,1) > 100.0 ) then
            rps = f_one / prsi(i,1)
          else
            print *,' !!! (1) Error in subr radiation_aerosols:',       &
     &              ' unrealistic surface pressure =', i,prsi(i,1)
            errflg = 1
            errmsg = 'ERROR(aer_property): Unrealistic surface pressure'
            return
          endif

          ii = 1
          do k = 1, NLAY
            if (prsi(i,k+1)*rps < sigref(ii,kp)) then
              ii = ii + 1
              if (ii == 2 .and. prsref(2,kp) == prsref(3,kp)) then
                ii = 3
              endif
            endif
            idmaer(k) = ii

            if ( ii > 1 ) then
              tmp1 = haer(ii,kp)
            else
              tmp1 = h1
            endif

            if (tmp1 > f_zero) then
              tmp2 = f_one / tmp1
              delz(k) = tmp1 * (exp(-hz(i,k)*tmp2)-exp(-hz(i,k+1)*tmp2))
            else
              delz(k) = dz1(k)
            endif
          enddo

        else  lab_if_flip                         ! input from toa to sfc

          if ( prsi(i,NLP1) > 100.0 ) then
            rps =  1.0 / prsi(i,NLP1)
          else
            print *,' !!! (2) Error in subr radiation_aerosols:',       &
     &              ' unrealistic surface pressure =', i,prsi(i,NLP1)
          endif

          ii = 1
          do k = NLAY, 1, -1
            if (prsi(i,k)*rps < sigref(ii,kp)) then
              ii = ii + 1
              if (ii == 2 .and. prsref(2,kp) == prsref(3,kp)) then
                ii = 3
              endif
            endif
            idmaer(k) = ii

            if ( ii > 1 ) then
              tmp1 = haer(ii,kp)
            else
              tmp1 = h1
            endif

            if (tmp1 > f_zero) then
              tmp2   = f_one / tmp1
              delz(k) = tmp1 * (exp(-hz(i,k+1)*tmp2)-exp(-hz(i,k)*tmp2))
            else
              delz(k) = dz1(k)
            endif
          enddo

        endif  lab_if_flip

!  ---  check print

!       print *,' in setclimaer, profile:',i
!       print *,'  rh   :',rh1
!       print *,'  dz   :',dz1
!       print *,'  delz :',delz
!       print *,'  idmaer:',idmaer

!> -# Call radclimaer() to calculate SW/LW aerosol optical properties
!!    for the corresponding frequency bands.

        call radclimaer(top_at_1)
!  ---  inputs:  (in-scope variables)
!  ---  outputs: (in-scope variables)

        if ( laersw ) then

          do m = 1, NBDSW
            do k = 1, NLAY
              aerosw(i,k,m,1) = tauae(k,m)
              aerosw(i,k,m,2) = ssaae(k,m)
              aerosw(i,k,m,3) = asyae(k,m)
            enddo
          enddo

!  ---  total aod (optional)
         do k = 1, NLAY
           aerodp(i,1) = aerodp(i,1) + tauae(k,nv_aod)
         enddo

!  ---  for diagnostic output (optional)
         do m = 1, NSPC
           aerodp(i,m+1) = spcodp(m)
         enddo

        endif     ! end if_larsw_block

        if ( laerlw ) then

          if ( NLWBND == 1 ) then
            m1 = NSWBND + 1
            do m = 1, NBDLW
              do k = 1, NLAY
                aerolw(i,k,m,1) = tauae(k,m1)
                aerolw(i,k,m,2) = ssaae(k,m1)
                aerolw(i,k,m,3) = asyae(k,m1)
              enddo
            enddo
          else
            do m = 1, NBDLW
              m1 = NSWBND + m
              do k = 1, NLAY
                aerolw(i,k,m,1) = tauae(k,m1)
                aerolw(i,k,m,2) = ssaae(k,m1)
                aerolw(i,k,m,3) = asyae(k,m1)
              enddo
            enddo
          endif

        endif     ! end if_laerlw_block

      enddo  lab_do_IMAX

! =================
      contains
! =================

!> This subroutine computes aerosols optical properties in NSWLWBD
!! bands. there are seven different vertical profile structures. in the
!! troposphere, aerosol distribution at each grid point is composed
!! from up to six components out of ten different substances.
!--------------------------------
      subroutine radclimaer(top_at_1)
!................................

!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

!  ==================================================================  !
!                                                                      !
!  compute aerosols optical properties in NSWLWBD bands. there are     !
!  seven different vertical profile structures. in the troposphere,    !
!  aerosol distribution at each grid point is composed from up to      !
!  six components out of a total of ten different substances.          !
!                                                                      !
!  ref: wmo report wcp-112 (1986)                                      !
!                                                                      !
!  input variables:                                                    !
!     cmix   - mixing ratioes of aerosol components  -     NCM         !
!     denn   - aerosol number densities              -     2           !
!     rh1    - relative humidity                     -     NLAY        !
!     delz   - effective layer thickness             km    NLAY        !
!     idmaer - aerosol domain index                  -     NLAY        !
!     NXC    - number of different aerosol components-     1           !
!     NLAY   - vertical dimensions                   -     1           !
!                                                                      !
!  output variables:                                                   !
!     tauae  - optical depth                         -     NLAY*NSWLWBD!
!     ssaae  - single scattering albedo              -     NLAY*NSWLWBD!
!     asyae  - asymmetry parameter                   -     NLAY*NSWLWBD!
!!    aerodp - vertically integrated aer-opt-depth   -     IMAX*NSPC+1 !
!                                                                      !
!  ==================================================================  !
!
      real (kind=kind_phys) :: crt1, crt2
      parameter (crt1=30.0, crt2=0.03333)

!  ---  inputs:
      logical, intent(in) :: top_at_1
!  ---  outputs:

!  ---  locals:
      real (kind=kind_phys) :: cm, hd, hdi, sig0u, sig0l, ratio, tt0,   &
     &      ex00, sc00, ss00, as00, ex01, sc01, ss01, as01,     tt1,    &
     &      ex02, sc02, ss02, as02, ex03, sc03, ss03, as03,     tt2,    &
     &      ext1, sca1, ssa1, asy1, drh0, drh1, rdrh

      integer :: ih1, ih2, kk, idom, icmp, ib, ii, ic, ic1
      integer :: idx

!===> ...  begin here

       spcodp = f_zero

!===> ... loop over vertical layers from top to surface

      lab_do_layer : do kk = 1, NLAY

! --- linear interp coeffs for rh-dep species

        ih2 = 1
        do while ( rh1(kk) > rhlev(ih2) )
          ih2 = ih2 + 1
          if ( ih2 > NRHLEV ) exit
        enddo
        ih1 = max( 1, ih2-1 )
        ih2 = min( NRHLEV, ih2 )

        drh0 = rhlev(ih2) - rhlev(ih1)
        drh1 = rh1(kk) - rhlev(ih1)
        if ( ih1 == ih2 ) then
          rdrh = f_zero
        else
          rdrh = drh1 / drh0
        endif

! --- assign optical properties in each domain

        idom = idmaer(kk)

        lab_if_idom : if (idom == 5) then
! --- 5th domain - upper stratosphere assume no aerosol

          do ib = 1, NSWLWBD
            tauae(kk,ib) = f_zero
            if ( ib <= NSWBND ) then
              ssaae(kk,ib) = 0.99
              asyae(kk,ib) = 0.696
            else
              ssaae(kk,ib) = 0.5
              asyae(kk,ib) = 0.3
            endif
          enddo

        elseif (idom == 4) then    lab_if_idom
! --- 4th domain - stratospheric layers

          do ib = 1, NSWLWBD
            tauae(kk,ib) = extstra(ib) * delz(kk)
            if ( ib <= NSWBND ) then
              ssaae(kk,ib) = 0.99
              asyae(kk,ib) = 0.696
            else
              ssaae(kk,ib) = 0.5
              asyae(kk,ib) = 0.3
            endif
          enddo

! --- compute aod from individual species' contribution (optional)
          idx = idxspc(10)             ! for sulfate
          spcodp(idx) = spcodp(idx) + tauae(kk,nv_aod)

        elseif (idom == 3) then    lab_if_idom
! --- 3rd domain - free tropospheric layers
!   1:inso 0.17e-3; 2:soot 0.4; 7:waso 0.59983; n:730

          do ib = 1, NSWLWBD
            ex01 = extrhi(1,ib)
            sc01 = scarhi(1,ib)
            ss01 = ssarhi(1,ib)
            as01 = asyrhi(1,ib)

            ex02 = extrhi(2,ib)
            sc02 = scarhi(2,ib)
            ss02 = ssarhi(2,ib)
            as02 = asyrhi(2,ib)

            ex03 = extrhd(ih1,1,ib)                                     &
     &           + rdrh * (extrhd(ih2,1,ib) - extrhd(ih1,1,ib))
            sc03 = scarhd(ih1,1,ib)                                     &
     &           + rdrh * (scarhd(ih2,1,ib) - scarhd(ih1,1,ib))
            ss03 = ssarhd(ih1,1,ib)                                     &
     &           + rdrh * (ssarhd(ih2,1,ib) - ssarhd(ih1,1,ib))
            as03 = asyrhd(ih1,1,ib)                                     &
     &           + rdrh * (asyrhd(ih2,1,ib) - asyrhd(ih1,1,ib))

            ext1 = 0.17e-3*ex01 + 0.4*ex02 + 0.59983*ex03
            sca1 = 0.17e-3*sc01 + 0.4*sc02 + 0.59983*sc03
            ssa1 = 0.17e-3*ss01*ex01 + 0.4*ss02*ex02 + 0.59983*ss03*ex03
            asy1 = 0.17e-3*as01*sc01 + 0.4*as02*sc02 + 0.59983*as03*sc03

            tauae(kk,ib) = ext1 * 730.0 * delz(kk)
            ssaae(kk,ib) = min(f_one, ssa1/ext1)
            asyae(kk,ib) = min(f_one, asy1/sca1)

! --- compute aod from individual species' contribution (optional)
            if ( ib==nv_aod ) then
             spcodp(1) = spcodp(1) + 0.17e-3*ex01*730.0*delz(kk)   ! dust (inso)   #1
             spcodp(2) = spcodp(2) + 0.4    *ex02*730.0*delz(kk)   ! black carbon  #2
             spcodp(3) = spcodp(3) + 0.59983*ex03*730.0*delz(kk)   ! water soluble #7
            endif

          enddo

        elseif (idom == 1) then    lab_if_idom
! --- 1st domain - mixing layer

          lab_do_ib : do ib = 1, NSWLWBD
            ext1 = f_zero
            sca1 = f_zero
            ssa1 = f_zero
            asy1 = f_zero

            lab_do_icmp : do icmp = 1, NCM
              ic = icmp
              idx = idxspc(icmp)

              cm = cmix(icmp)
              lab_if_cm : if ( cm > f_zero ) then

                lab_if_ic : if ( ic <= NCM1 ) then        ! component withour rh dep
                  tt0  = cm * extrhi(ic,ib)
                  ext1 = ext1 + tt0
                  sca1 = sca1 + cm * scarhi(ic,ib)
                  ssa1 = ssa1 + cm * ssarhi(ic,ib) * extrhi(ic,ib)
                  asy1 = asy1 + cm * asyrhi(ic,ib) * scarhi(ic,ib)
                else  lab_if_ic                           ! component with rh dep
                  ic1 = ic - NCM1

                  ex00 = extrhd(ih1,ic1,ib)                             &
     &               + rdrh * (extrhd(ih2,ic1,ib) - extrhd(ih1,ic1,ib))
                  sc00 = scarhd(ih1,ic1,ib)                             &
     &               + rdrh * (scarhd(ih2,ic1,ib) - scarhd(ih1,ic1,ib))
                  ss00 = ssarhd(ih1,ic1,ib)                             &
     &               + rdrh * (ssarhd(ih2,ic1,ib) - ssarhd(ih1,ic1,ib))
                  as00 = asyrhd(ih1,ic1,ib)                             &
     &               + rdrh * (asyrhd(ih2,ic1,ib) - asyrhd(ih1,ic1,ib))

                  tt0  = cm * ex00
                  ext1 = ext1 + tt0
                  sca1 = sca1 + cm * sc00
                  ssa1 = ssa1 + cm * ss00 * ex00
                  asy1 = asy1 + cm * as00 * sc00
                endif  lab_if_ic

! --- compute aod from individual species' contribution (optional)
                if ( ib==nv_aod ) then
                 spcodp(idx) = spcodp(idx) + tt0*denn(1)*delz(kk)   ! idx for dif species
                endif

              endif  lab_if_cm
            enddo  lab_do_icmp

            tauae(kk,ib) = ext1 * denn(1) * delz(kk)
            ssaae(kk,ib) = min(f_one, ssa1/ext1)
            asyae(kk,ib) = min(f_one, asy1/sca1)
          enddo  lab_do_ib

        elseif (idom == 2) then    lab_if_idom
! --- 2nd domain - mineral transport layers

          do ib = 1, NSWLWBD
            tauae(kk,ib) = extrhi(6,ib) * denn(2) * delz(kk)
            ssaae(kk,ib) = ssarhi(6,ib)
            asyae(kk,ib) = asyrhi(6,ib)
          enddo

! --- compute aod from individual species' contribution (optional)
          spcodp(1) = spcodp(1) + tauae(kk,nv_aod)            ! dust

        else  lab_if_idom
! --- domain index out off range, assume no aerosol

          do ib = 1, NSWLWBD
            tauae(kk,ib) = f_zero
            ssaae(kk,ib) = f_one
            asyae(kk,ib) = f_zero
          enddo

!         write(6,19) kk,idom
! 19      format(/'  ***  ERROR in sub AEROS: domain index out'         &
!    &,            ' of range!  K, IDOM =',3i5,' ***')
!         stop 19

        endif  lab_if_idom

      enddo  lab_do_layer

!
!===> ... smooth profile at domain boundaries
!
      if (top_at_1) then    ! input from toa to sfc

        do ib = 1, NSWLWBD
        do kk = 2, NLAY
          if ( tauae(kk,ib) > f_zero ) then
            ratio = tauae(kk-1,ib) / tauae(kk,ib)
          else
            ratio = f_one
          endif

          tt0 = tauae(kk,ib) + tauae(kk-1,ib)
          tt1 = 0.2 * tt0
          tt2 = tt0 - tt1

          if ( ratio > crt1 ) then
            tauae(kk,ib)   = tt1
            tauae(kk-1,ib) = tt2
          endif

          if ( ratio < crt2 ) then
            tauae(kk,ib)   = tt2
            tauae(kk-1,ib) = tt1
          endif
        enddo   ! do_kk_loop
        enddo   ! do_ib_loop

      else                      ! input from sfc to toa

        do ib = 1, NSWLWBD
        do kk = NLAY-1, 1, -1
          if ( tauae(kk,ib) > f_zero ) then
            ratio = tauae(kk+1,ib) / tauae(kk,ib)
          else
            ratio = f_one
          endif

          tt0 = tauae(kk,ib) + tauae(kk+1,ib)
          tt1 = 0.2 * tt0
          tt2 = tt0 - tt1

          if ( ratio > crt1 ) then
            tauae(kk,ib)   = tt1
            tauae(kk+1,ib) = tt2
          endif

          if ( ratio < crt2 ) then
            tauae(kk,ib)   = tt2
            tauae(kk+1,ib) = tt1
          endif
        enddo   ! do_kk_loop
        enddo   ! do_ib_loop

      endif

!
      return
!................................
      end subroutine radclimaer
!--------------------------------
!
!...................................
      end subroutine aer_property
!-----------------------------------

!> This subroutine is the gocart aerosol initialization
!! program to set up necessary parameters and working arrays.
!>\param solfwv     (NWVTOT), solar flux for each individual wavenumber
!!                  \f$(w/m^2)\f$
!!\param eirfwv     (NWVTIR), IR flux(273k) for each individual wavenumber
!!                  \f$(w/m^2)\f$
!!\param me         print message control flag
!!
!>\section gel_go_ini General Algorithm
!-----------------------------------
      subroutine gocart_aerinit                                         &
     &     ( solfwv, eirfwv, me,                                        &
     &     errflg, errmsg)

!  ==================================================================  !
!                                                                      !
!  subprogram : gocart_aerinit                                         !
!                                                                      !
!  gocart_aerinit is the gocart aerosol initialization program         !
!  to set up necessary parameters and working arrays.                  !
!                                                                      !
!  inputs:                                                             !
!   solfwv(NWVTOT)   - solar flux for each individual wavenumber (w/m2)!
!   eirfwv(NWVTIR)   - ir flux(273k) for each individual wavenum (w/m2)!
!   me               - print message control flag                      !
!                                                                      !
!  outputs: (CCPP error handling)                                      !
!   errflg           - CCPP error flag                                 !
!   errmsg           - CCPP error message                              !
!                                                                      !
!  module variables:                                                   !
!     NWVSOL  - num of wvnum regions where solar flux is constant      !
!     NWVTOT  - total num of wave numbers used in sw spectrum          !
!     NWVTIR  - total num of wave numbers used in the ir region        !
!     NSWBND  - total number of sw spectral bands                      !
!     NLWBND  - total number of lw spectral bands                      !
!     NAERBND - number of bands for climatology aerosol data           !
!     KCM1    - number of rh independent aeros species                 !
!     KCM2    - number of rh dependent aeros species                   !
!                                                                      !
!  usage:    call gocart_init                                          !
!                                                                      !
!  subprograms called:  rd_gocart_luts, optavg_gocart                  !
!                                                                      !
!  ==================================================================  !

      implicit none

!  ---  inputs:
      real (kind=kind_phys), dimension(:) :: solfwv        ! one wvn sol flux
      real (kind=kind_phys), dimension(:) :: eirfwv        ! one wvn ir flux

      integer,  intent(in) :: me

!  ---  output: (CCPP error handling)
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      real (kind=kind_phys), dimension(kaerbndi,kcm1)       ::          &
     &       rhidext0_grt, rhidsca0_grt, rhidssa0_grt, rhidasy0_grt
      real (kind=kind_phys), dimension(kaerbndd,krhlev,kcm2)::          &
     &       rhdpext0_grt, rhdpsca0_grt, rhdpssa0_grt, rhdpasy0_grt

      real (kind=kind_phys), dimension(nswbnd,kaerbndd) :: solwaer
      real (kind=kind_phys), dimension(nswbnd)          :: solbnd
      real (kind=kind_phys), dimension(nlwbnd,kaerbndd) :: eirwaer
      real (kind=kind_phys), dimension(nlwbnd)          :: eirbnd

      real (kind=kind_phys), dimension(nswbnd,kaerbndi) :: solwaer_du
      real (kind=kind_phys), dimension(nswbnd)          :: solbnd_du
      real (kind=kind_phys), dimension(nlwbnd,kaerbndi) :: eirwaer_du
      real (kind=kind_phys), dimension(nlwbnd)          :: eirbnd_du

      integer, dimension(nswbnd) :: nv1, nv2, nv1_du, nv2_du
      integer, dimension(nlwbnd) :: nr1, nr2, nr1_du, nr2_du

      integer, dimension(kaerbndd) :: iendwv
      integer, dimension(kaerbndi) :: iendwv_du
      real (kind=kind_phys), dimension(kaerbndd) :: wavelength
      real (kind=kind_phys), dimension(kaerbndi) :: wavelength_du
      real (kind=kind_phys) :: sumsol, sumir, sumsol_du, sumir_du

      integer :: i, j, k, mb, ib, ii, iix, iw, iw1, iw2

!
!===>  ...  begin here

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

!
!  --- ...  invoke gocart aerosol initialization


      if (KCM /= ntrcaerm ) then
        print *, 'ERROR in # of gocart aer species',KCM
        errflg = 1
        errmsg = 'ERROR(gocart_init): Incorrect # of species'
        return
      endif

!  --- ...  aloocate and input aerosol optical data

      if ( .not. allocated( extrhi_grt ) ) then
        allocate ( extrhi_grt (       kcm1,nswlwbd) )
        allocate ( extrhi_grt_550 (       kcm1,1) )
        allocate ( scarhi_grt (       kcm1,nswlwbd) )
        allocate ( ssarhi_grt (       kcm1,nswlwbd) )
        allocate ( asyrhi_grt (       kcm1,nswlwbd) )
        allocate ( extrhd_grt (krhlev,kcm2,nswlwbd) )
        allocate ( extrhd_grt_550 (krhlev,kcm2,1) )
        allocate ( scarhd_grt (krhlev,kcm2,nswlwbd) )
        allocate ( ssarhd_grt (krhlev,kcm2,nswlwbd) )
        allocate ( asyrhd_grt (krhlev,kcm2,nswlwbd) )
      endif

!  --- ...  read tabulated GOCART aerosols optical data

      call rd_gocart_luts
!  ---  inputs:  (in scope variables, module variables)
!  ---  outputs: (in scope variables)

!  --- ...  convert wavelength to wavenumber
!           wavelength and wavelength_du are read-in by rd_gocart_luts

      do i = 1, kaerbndd
       iendwv(i) = int(10000. / wavelength(i))
      enddo

      do i = 1, kaerbndi
       iendwv_du(i) = int(10000. / wavelength_du(i))
      enddo

!  --- ...  compute solar flux weights and interval indices for mapping
!           spectral bands between sw radiation and aerosol data

      if ( laswflg ) then
        solbnd (:)   = f_zero
        solbnd_du (:)= f_zero
        do i=1,nswbnd
          do j=1,kaerbndd
            solwaer(i,j) = f_zero
          enddo
          do j=1,kaerbndi
            solwaer_du(i,j) = f_zero
          enddo
        enddo

        do ib = 1, nswbnd
          mb  = ib + nswstr - 1
          ii  = 1
          iix = 1
          iw1 = nint(wvnsw1(mb))
          iw2 = nint(wvnsw2(mb))

          if ( wvnsw2(mb)>=wvn550 .and. wvn550>=wvnsw1(mb) ) then
            nv_aod = ib                  ! sw band number covering 550nm wavelenth
          endif

! -- for rd-dependent
          do while ( iw1 > iendwv(ii) )
            if ( ii == kaerbndd ) exit
            ii = ii + 1
          enddo
          sumsol = f_zero
          nv1(ib) = ii

! -- for rd-independent
          do while ( iw1 > iendwv_du(iix) )
            if ( iix == kaerbndi ) exit
            iix = iix + 1
          enddo
          sumsol_du = f_zero
          nv1_du(ib) = iix

          do iw = iw1, iw2
! -- for rd-dependent
            solbnd(ib) = solbnd(ib) + solfwv(iw)
            sumsol = sumsol + solfwv(iw)

            if ( iw == iendwv(ii) ) then
              solwaer(ib,ii) = sumsol
              if ( ii < kaerbndd ) then
                sumsol = f_zero
                ii = ii + 1
              endif
            endif

! -- for rd-independent
            solbnd_du(ib) = solbnd_du(ib) + solfwv(iw)
            sumsol_du = sumsol_du + solfwv(iw)

            if ( iw == iendwv_du(iix) ) then
              solwaer_du(ib,iix) = sumsol_du
              if ( iix < kaerbndi ) then
                sumsol_du = f_zero
                iix = iix + 1
              endif
            endif
          enddo

          if ( iw2 /= iendwv(ii) ) then
            solwaer(ib,ii) = sumsol
          endif
          if ( iw2 /= iendwv_du(iix) ) then
            solwaer_du(ib,iix) = sumsol_du
          endif

          nv2(ib) = ii
          nv2_du(ib) = iix
        enddo     ! end do_ib_block for sw
      endif    ! end if_laswflg_block

!  --- ...  compute lw flux weights and interval indices for mapping
!           spectral bands between lw radiation and aerosol data

      if ( lalwflg ) then
        eirbnd    (:)   = f_zero
        eirbnd_du (:)   = f_zero
        do i=1,nlwbnd
          do j=1,kaerbndd
            eirwaer(i,j) = f_zero
          enddo
          do j=1,kaerbndi
            eirwaer_du(i,j) = f_zero
          enddo
        enddo

        do ib = 1, nlwbnd
          ii = 1
          iix = 1
          if ( nlwbnd == 1 ) then
            iw1 = 400                   ! corresponding 25 mu
            iw2 = 2500                  ! corresponding 4  mu
          else
            mb = ib + nlwstr - 1
            iw1 = nint(wvnlw1(mb))
            iw2 = nint(wvnlw2(mb))
          endif

! -- for rd-dependent
          do while ( iw1 > iendwv(ii) )
            if ( ii == kaerbndd ) exit
            ii = ii + 1
          enddo
          sumir = f_zero
          nr1(ib) = ii

! -- for rd-independent
          do while ( iw1 > iendwv_du(iix) )
            if ( iix == kaerbndi ) exit
            iix = iix + 1
          enddo
          sumir_du = f_zero
          nr1_du(ib) = iix

          do iw = iw1, iw2
! -- for rd-dependent
            eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
            sumir  = sumir  + eirfwv(iw)

            if ( iw == iendwv(ii) ) then
              eirwaer(ib,ii) = sumir

              if ( ii < kaerbndd ) then
                sumir = f_zero
                ii = ii + 1
              endif
            endif

! -- for rd-independent
            eirbnd_du(ib) = eirbnd_du(ib) + eirfwv(iw)
            sumir_du  = sumir_du  + eirfwv(iw)

            if ( iw == iendwv_du(iix) ) then
              eirwaer_du(ib,iix) = sumir_du

              if ( iix < kaerbndi ) then
                sumir_du = f_zero
                iix = iix + 1
              endif
            endif
          enddo

          if ( iw2 /= iendwv(ii) ) then
            eirwaer(ib,ii) = sumir
          endif
          if ( iw2 /= iendwv_du(iix) ) then
            eirwaer_du(ib,iix) = sumir_du
          endif

          nr2(ib) = ii
          nr2_du(ib) = iix
        enddo     ! end do_ib_block for lw
      endif    ! end if_lalwflg_block

!  ---  compute spectral band mean properties for each species

      call optavg_gocart
!  ---  inputs:  (in-scope variables, module variables)
!  ---  outputs: (module variables)


!  ---  check print
!      if (me == 0) then
!      do ib = 1, NSWBND
!        mb = ib + NSWSTR - 1
!        print *, '  wvnsw1,wvnsw2 :',wvnsw1(mb),wvnsw2(mb)
!        print *, ' After optavg_gocart, for sw band:',ib
!        print *, '  extrhi:', extrhi_grt(:,ib)
!        print *, '  scarhi:', scarhi_grt(:,ib)
!        print *, '  ssarhi:', ssarhi_grt(:,ib)
!        print *, '  asyrhi:', asyrhi_grt(:,ib)
!       do i = 1, KRHLEV
!        print *, '  extrhd for rhlev:',i
!        print *,  extrhd_grt(i,:,ib)
!        print *, '  scarhd for rhlev:',i
!        print *, scarhd_grt(i,:,ib)
!        print *, '  ssarhd for rhlev:',i
!        print *, ssarhd_grt(i,:,ib)
!        print *, '  asyrhd for rhlev:',i
!        print *, asyrhd_grt(i,:,ib)
!        enddo
!      enddo
!      print *, '  wvnlw1 :',wvnlw1
!      print *, '  wvnlw2 :',wvnlw2
!      do ib = 1, NLWBND
!        ii = NSWBND + ib
!        print *,' After optavg_gocart, for lw band:',ib
!        print *,'  extrhi_grt:', extrhi_grt(:,ii)
!        print *,'  scarhi_grt:', scarhi_grt(:,ii)
!        print *,'  ssarhi_grt:', ssarhi_grt(:,ii)
!        print *,'  asyrhi_grt:', asyrhi_grt(:,ii)
!       do i = 1, KRHLEV
!        print *,'  extrhd for rhlev:',i
!        print *, extrhd_grt(i,:,ib)
!        print *,'  scarhd for rhlev:',i
!        print *, scarhd_grt(i,:,ib)
!        print *,'  ssarhd for rhlev:',i
!        print *, ssarhd_grt(i,:,ib)
!        print *,'  asyrhd for rhlev:',i
!        print *, asyrhd_grt(i,:,ib)
!        enddo
!      enddo
!      endif

! =================
      contains
! =================

!-----------------------------
!>   read GMAO pre-tabultaed aerosol optical data for dust, seasalt,    !
!!   sulfate, black carbon, and organic carbon aerosols 
      subroutine rd_gocart_luts
!.............................
!  ---  inputs:  (in scope variables, module variables)
!  ---  outputs: (in scope variables)

! ==================================================================== !
!                                                                      !
! subprogram: rd_gocart_luts                                           !
!   read GMAO pre-tabultaed aerosol optical data for dust, seasalt,    !
!   sulfate, black carbon, and organic carbon aerosols                 !
!                                                                      !
!  major local variables:                                              !
!   for handling spectral band structures                              !
!     iendwv    - ending wvnum (cm**-1) for each band    kaerbndd      !
!     iendwv_du - ending wvnum (cm**-1) for each band    kaerbndi      !
!   for handling optical properties of rh independent species (kcm1)   !
!     1=du001, 2=du002, 3=du003, 4=du004, 5=du005                      !
!     rhidext0_grt - extinction coefficient          kaerbndi*kcm1     !
!     rhidsca0_grt - scattering coefficient          kaerbndi*kcm1     !
!     rhidssa0_grt - single scattering albedo        kaerbndi*kcm1     !
!     rhidasy0_grt - asymmetry parameter             kaerbndi*kcm1     !
!   for handling optical properties of rh ndependent species (kcm2)    !
!     1=ss001, 2=ss002, 3=ss003, 4=ss004, 5=ss005, 6=so4,              !
!     7=bcphobic, 8=bcphilic, 9=ocphobic, 10=ocphilic                  !
!     rhdpext0_grt - extinction coefficient        kaerbndd*krhlev*kcm2!
!     rhdpsca0_grt - scattering coefficient        kaerbndd*krhlev*kcm2!
!     rhdpssa0_grt - single scattering albedo      kaerbndd*krhlev*kcm2!
!     rhdpasy0_grt - asymmetry parameter           kaerbndd*krhlev*kcm2!
!                                                                      !
!  usage:    call rd_gocart_luts                                       !
!                                                                      !
!  ==================================================================  !
!
      implicit none

!  ---  inputs: (none)
!  ---  output: (none)

!  ---  locals:
       integer             :: iradius, ik, ibeg
       integer, parameter  :: numspc = 5        ! # of aerosol species

! - input tabulated aerosol optical spectral data from GSFC
       real, dimension(kaerbndd) :: lambda             ! wavelength (m) for non-dust
       real, dimension(kaerbndi) :: lambda_du          ! wavelength (m) for dust
       real, dimension(krhlev)   :: rh                 ! relative humidity (fraction)
       real, dimension(kaerbndd,krhlev,numspc) :: bext! extinction efficiency (m2/kg)
       real, dimension(kaerbndd,krhlev,numspc) :: bsca! scattering efficiency (m2/kg)
       real, dimension(kaerbndd,krhlev,numspc) :: g   ! asymmetry factor (dimensionless)
       real, dimension(kaerbndi,krhlev,numspc) :: bext_du! extinction efficiency (m2/kg)
       real, dimension(kaerbndi,krhlev,numspc) :: bsca_du! scattering efficiency (m2/kg)
       real, dimension(kaerbndi,krhlev,numspc) :: g_du   ! asymmetry factor (dimensionless)
!
       logical           :: file_exist
       character*50      :: fin, dummy

! ---  read LUTs for dust aerosols
       fin='optics_'//gridcomp(1)//'.dat'
       inquire (file=trim(fin), exist=file_exist)
       if ( file_exist ) then
         close(niaercm)
         open (unit=niaercm, file=fin, status='OLD')
         rewind(niaercm)
       else
         print *,' Requested luts file ',trim(fin),' not found'
         print *,' ** Stopped in rd_gocart_luts ** '
         errflg = 1
         errmsg = 'Requested luts file '//trim(fin)//' not found'
         return
       endif      ! end if_file_exist_block

       iradius = 5
! read lambda and compute mpwavelength (m)
       read(niaercm,'(a40)') dummy
       read(niaercm,*) (lambda_du(i), i=1, kaerbndi)
! read rh, relative humidity (fraction)
       read(niaercm,'(a40)') dummy
       read(niaercm,*) (rh(i), i=1, krhlev)
! read bext (m2 (kg dry mass)-1)
       do k = 1, iradius
         read(niaercm,'(a40)') dummy
         do j=1, krhlev
          read(niaercm,*) (bext_du(i,j,k), i=1,kaerbndi)
         enddo
       enddo
! read bsca (m2 (kg dry mass)-1)
       do k = 1, iradius
         read(niaercm,'(a40)') dummy
         do j=1, krhlev
          read(niaercm,*) (bsca_du(i,j,k), i=1, kaerbndi)
         enddo
       enddo
! read g (dimensionless)
       do k = 1, iradius
         read(niaercm,'(a40)') dummy
         do j=1, krhlev
          read(niaercm,*) (g_du(i,j,k), i=1, kaerbndi)
         enddo
       enddo

! fill rhidext0 local arrays for dust aerosols (flip i-index)
       do i = 1, kaerbndi         ! convert from m to micron
          j = kaerbndi -i + 1     ! flip i-index
          wavelength_du(j) = 1.e6 * lambda_du(i)
           if (int(wavelength_du(j)*100) == 55) then
              id550=j
           endif
       enddo
       do k = 1, iradius
          do i = 1, kaerbndi
           ii = kaerbndi -i + 1
           rhidext0_grt(ii,k) = bext_du(i,1,k)
           rhidsca0_grt(ii,k) = bsca_du(i,1,k)
           if ( bext_du(i,1,k) /= f_zero) then
            rhidssa0_grt(ii,k) = bsca_du(i,1,k)/bext_du(i,1,k)
           else
            rhidssa0_grt(ii,k) = f_one
           endif
           rhidasy0_grt(ii,k) = g_du(i,1,k)
          enddo
       enddo

! ---  read LUTs for non-dust aerosols
       do ib = 2, num_gc       ! loop thru SS, SU, BC, OC
        fin='optics_'//gridcomp(ib)//'.dat'
        inquire (file=trim(fin), exist=file_exist)
        if ( file_exist ) then
          close(niaercm)
          open (unit=niaercm, file=fin, status='OLD')
          rewind(niaercm)
        else
          print *,' Requested luts file ',trim(fin),' not found'
          print *,' ** Stopped in rd_gocart_luts ** '
          errflg = 1
          errmsg = 'Requested luts file '//trim(fin)//' not found'
          return
        endif      ! end if_file_exist_block

        ibeg  =  radius_lower(ib) - kcm1
        iradius = num_radius(ib)

! read lambda and compute mpwavelength (m)
        read(niaercm,'(a40)') dummy
        read(niaercm,*) (lambda(i), i=1, kaerbndd)
! read rh, relative humidity (fraction)
        read(niaercm,'(a40)') dummy
        read(niaercm,*) (rh(i), i=1, krhlev)
! read bext
        do k = 1, iradius
         read(niaercm,'(a40)') dummy
         do j=1, krhlev
          read(niaercm,*) (bext(i,j,k), i=1,kaerbndd)
         enddo
        enddo
! read bsca
        do k = 1, iradius
         read(niaercm,'(a40)') dummy
         do j=1, krhlev
          read(niaercm,*) (bsca(i,j,k), i=1, kaerbndd)
         enddo
        enddo
! read g
        do k = 1, iradius
         read(niaercm,'(a40)') dummy
         do j=1, krhlev
          read(niaercm,*) (g(i,j,k), i=1, kaerbndd)
         enddo
        enddo

! fill rhdpext0 local arrays for non-dust aerosols (flip i-index)
        do i = 1, kaerbndd         ! convert from m to micron
           j = kaerbndd -i + 1     ! flip i-index
           wavelength(j) = 1.e6 * lambda(i)
           if (int(wavelength(j)*100) == 55) then
              i550=j
           endif
        enddo
        do k = 1, iradius
          ik = ibeg + k - 1
          do i = 1, kaerbndd
          ii = kaerbndd -i + 1
          do j = 1, krhlev
           rhdpext0_grt(ii,j,ik) = bext(i,j,k)
           rhdpsca0_grt(ii,j,ik) = bsca(i,j,k)
           if ( bext(i,j,k) /= f_zero) then
            rhdpssa0_grt(ii,j,ik) = bsca(i,j,k)/bext(i,j,k)
           else
            rhdpssa0_grt(ii,j,ik) = f_one
           endif
           rhdpasy0_grt(ii,j,ik) = g(i,j,k)
          enddo
          enddo
        enddo

       enddo       !! ib-loop

      return
!...................................
      end subroutine rd_gocart_luts
!-----------------------------------

!--------------------------------
!>   compute mean aerosol optical properties over each sw radiation     
!!   spectral band for each of the species components.  This program    
!!   follows optavg routine (in turn follows gfdl's approach for thick  
!!   cloud opertical property in sw radiation scheme (2000).           
      subroutine optavg_gocart
!................................
!  ---  inputs:  (in-scope variables, module variables)
!  ---  outputs: (module variables)

! ==================================================================== !
!                                                                      !
! subprogram: optavg_gocart                                            !
!                                                                      !
!   compute mean aerosol optical properties over each sw radiation     !
!   spectral band for each of the species components.  This program    !
!   follows optavg routine (in turn follows gfdl's approach for thick  !
!   cloud opertical property in sw radiation scheme (2000).            !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
! major input variables:                                               !
!   nv1,nv2 (nswbnd) - start/end spectral band indices of aerosol data !
!                      for each sw radiation spectral band             !
!   nr1,nr2 (nlwbnd) - start/end spectral band indices of aerosol data !
!                      for each ir radiation spectral band             !
!   nv1_du,nv2_du(nswbnd) - start/end spectral band indices of aer data!
!                      for each sw radiation spectral band             !
!   nr1_du,nr2_du(nlwbnd) - start/end spectral band indices of aer data!
!                      for each ir radiation spectral band             !
!   solwaer (nswbnd,kaerbndd)                                          !
!                    - solar flux weight over each sw radiation band   !
!                      vs each aerosol data spectral band              !
!   eirwaer (nlwbnd,kaerbndd)                                          !
!                    - ir flux weight over each lw radiation band      !
!                      vs each aerosol data spectral band              !
!   solwaer_du (nswbnd,kaerbndi)                                       !
!                    - solar flux weight over each sw radiation band   !
!                      vs each aerosol data spectral band              !
!   eirwaer_du (nlwbnd,kaerbndi)                                       !
!                    - ir flux weight over each lw radiation band      !
!                      vs each aerosol data spectral band              !
!   solbnd  (nswbnd) - solar flux weight over each sw radiation band   !
!   eirbnd  (nlwbnd) - ir flux weight over each lw radiation band      !
!   solbnd_du(nswbnd) - solar flux weight over each sw radiation band  !
!   eirbnd_du(nlwbnd) - ir flux weight over each lw radiation band     !
!   nswbnd           - total number of sw spectral bands               !
!   nlwbnd           - total number of lw spectral bands               !
!                                                                      !
! external module variables:                                           !
!   laswflg          - control flag for sw spectral region             !
!   lalwflg          - control flag for lw spectral region             !
!                                                                      !
! output variables: (to module variables)                              !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
!  ---  output:

!  ---  locals:
      real (kind=kind_phys) :: sumk, sums, sumok, sumokg, sumreft,      &
     &       sp, refb, reft, rsolbd, rirbd

      integer :: ib, nb, ni, nh, nc
!
!===> ...  begin here
!
!  --- ...  loop for each sw radiation spectral band

      if ( laswflg ) then
        do nb = 1, nswbnd
          rsolbd = f_one / solbnd_du(nb)
          do nc = 1, kcm1         !  ---  for rh independent aerosol species
            sumk    = f_zero
            sums    = f_zero
            sumok   = f_zero
            sumokg  = f_zero
            sumreft = f_zero

            do ni = nv1_du(nb), nv2_du(nb)
              sp   = sqrt( (f_one - rhidssa0_grt(ni,nc))                &
     &             / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
              reft = (f_one - sp) / (f_one + sp)
              sumreft = sumreft + reft*solwaer_du(nb,ni)

              sumk    = sumk    + rhidext0_grt(ni,nc)*solwaer_du(nb,ni)
              sums    = sums    + rhidsca0_grt(ni,nc)*solwaer_du(nb,ni)
              sumok   = sumok   + rhidssa0_grt(ni,nc)*solwaer_du(nb,ni) &
     &                * rhidext0_grt(ni,nc)
              sumokg  = sumokg  + rhidssa0_grt(ni,nc)*solwaer_du(nb,ni) &
     &                * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
            enddo

            refb = sumreft * rsolbd

            extrhi_grt(nc,nb) = sumk   * rsolbd
            if (nb==nv_aod) then
              extrhi_grt_550(nc,1) = rhidext0_grt(id550,nc)
            endif
            scarhi_grt(nc,nb) = sums   * rsolbd
            asyrhi_grt(nc,nb) = sumokg / (sumok + 1.0e-10)
            ssarhi_grt(nc,nb) = 4.0*refb                                &
     &         / ( (f_one+refb)**2 - asyrhi_grt(nc,nb)*(f_one-refb)**2 )
          enddo   ! end do_nc_block for rh-ind aeros

          rsolbd = f_one / solbnd(nb)
          do nc = 1, kcm2         !  ---  for rh dependent aerosol species
            do nh = 1, krhlev
              sumk    = f_zero
              sums    = f_zero
              sumok   = f_zero
              sumokg  = f_zero
              sumreft = f_zero

              do ni = nv1(nb), nv2(nb)
                sp   = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc))           &
     &          /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)))
                reft = (f_one - sp) / (f_one + sp)
                sumreft = sumreft + reft*solwaer(nb,ni)

                sumk    = sumk  + rhdpext0_grt(ni,nh,nc)*solwaer(nb,ni)
                sums    = sums  + rhdpsca0_grt(ni,nh,nc)*solwaer(nb,ni)
                sumok   = sumok + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni) &
     &                  * rhdpext0_grt(ni,nh,nc)
                sumokg  = sumokg + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni)&
     &                  * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
              enddo

              refb = sumreft * rsolbd

              extrhd_grt(nh,nc,nb) = sumk   * rsolbd
              if (nb==nv_aod) then
                extrhd_grt_550(nh,nc,1) = rhdpext0_grt(i550,nh,nc)
              endif
              scarhd_grt(nh,nc,nb) = sums   * rsolbd
              asyrhd_grt(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
              ssarhd_grt(nh,nc,nb) = 4.0*refb                           &
     &         /((f_one+refb)**2 - asyrhd_grt(nh,nc,nb)*(f_one-refb)**2)

            enddo   ! end do_nh_block
          enddo   ! end do_nc_block for rh-dep aeros

       enddo   !  end do_nb_block for sw
      endif   !  end if_laswflg_block

!  --- ...  loop for each lw radiation spectral band

      if ( lalwflg ) then

        do nb = 1, nlwbnd

          ib = nswbnd + nb

          rirbd = f_one / eirbnd_du(nb)
          do nc = 1, kcm1         !  ---  for rh independent aerosol species
            sumk    = f_zero
            sums    = f_zero
            sumok   = f_zero
            sumokg  = f_zero
            sumreft = f_zero

            do ni = nr1_du(nb), nr2_du(nb)
              sp   = sqrt( (f_one - rhidssa0_grt(ni,nc))                &
     &             / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
              reft = (f_one - sp) / (f_one + sp)
              sumreft = sumreft + reft*eirwaer_du(nb,ni)

              sumk    = sumk    + rhidext0_grt(ni,nc)*eirwaer_du(nb,ni)
              sums    = sums    + rhidsca0_grt(ni,nc)*eirwaer_du(nb,ni)
              sumok   = sumok   + rhidssa0_grt(ni,nc)*eirwaer_du(nb,ni) &
     &                * rhidext0_grt(ni,nc)
              sumokg  = sumokg  + rhidssa0_grt(ni,nc)*eirwaer_du(nb,ni) &
     &                * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
            enddo

            refb = sumreft * rirbd

            extrhi_grt(nc,ib) = sumk   * rirbd
            scarhi_grt(nc,ib) = sums   * rirbd
            asyrhi_grt(nc,ib) = sumokg / (sumok + 1.0e-10)
            ssarhi_grt(nc,ib) = 4.0*refb                                &
     &         / ( (f_one+refb)**2 - asyrhi_grt(nc,ib)*(f_one-refb)**2 )

          enddo   ! end do_nc_block for rh-ind aeros

          rirbd = f_one / eirbnd(nb)
          do nc = 1, kcm2        !  ---  for rh dependent aerosol species
            do nh = 1, krhlev
              sumk    = f_zero
              sums    = f_zero
              sumok   = f_zero
              sumokg  = f_zero
              sumreft = f_zero

              do ni = nr1(nb), nr2(nb)
                sp   = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc))           &
     &          /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)))
                reft = (f_one - sp) / (f_one + sp)
                sumreft = sumreft + reft*eirwaer(nb,ni)

                sumk    = sumk  + rhdpext0_grt(ni,nh,nc)*eirwaer(nb,ni)
                sums    = sums  + rhdpsca0_grt(ni,nh,nc)*eirwaer(nb,ni)
                sumok   = sumok + rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) &
     &                  * rhdpext0_grt(ni,nh,nc)
                sumokg  = sumokg+ rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) &
     &                  * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
              enddo

              refb = sumreft * rirbd

              extrhd_grt(nh,nc,ib) = sumk   * rirbd
              scarhd_grt(nh,nc,ib) = sums   * rirbd
              asyrhd_grt(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
              ssarhd_grt(nh,nc,ib) = 4.0*refb                           &
     &         /((f_one+refb)**2 - asyrhd_grt(nh,nc,ib)*(f_one-refb)**2)
            enddo   ! end do_nh_block
          enddo   ! end do_nc_block for rh-dep aeros

        enddo   !  end do_nb_block for lw
      endif   !  end if_lalwflg_block
!
      return
      return
!...................................
      end subroutine optavg_gocart
!-----------------------------------

!...................................
      end subroutine gocart_aerinit
!-----------------------------------

!> This subroutine compute aerosol optical properties for SW
!! and LW radiations.
!!\param prsi           (IMAX,NLP1), pressure at interface in mb
!!\param prsl           (IMAX,NLAY), layer mean pressure(not used)
!!\param prslk          (IMAX,NLAY), exner function=\f$(p/p0)^{rocp}\f$ (not used)
!!\param tvly           (IMAX,NLAY), layer virtual temperature (not used)
!!\param rhlay          (IMAX,NLAY), layer mean relative humidity
!!\param dz             (IMAX,NLAY), layer thickness in m
!!\param hz             (IMAX,NLP1), level high in m
!!\param tracer         (IMAX,NLAY,NTRAC), aer tracer concentrations
!!\param aerfld         (IMAX,NLAY,NTRCAER), aer tracer concentrations
!!\param alon, alat     (IMAX), longitude and latitude of given points in degree
!!\param slmsk          (IMAX), sea/land mask (sea:0,land:1,sea-ice:2)
!!\param laersw,laerlw  logical flag for sw/lw aerosol calculations
!!\param IMAX           horizontal dimension of arrays
!!\param NLAY,NLP1      vertical dimensions of arrays
!!\param NSPC           num of species for optional aod output fields
!!\param aerosw        (IMAX,NLAY,NBDSW,NF_AESW), aeros opt properties for sw
!!\n                              (:,:,:,1): optical depth
!!\n                              (:,:,:,2): single scattering albedo
!!\n                              (:,:,:,3): asymmetry parameter
!!\param aerolw        (IMAX,NLAY,NBDLW,NF_AELW), aeros opt properties for lw
!!\n                              (:,:,:,1): optical depth
!!\n                              (:,:,:,2): single scattering albedo
!!\n                              (:,:,:,3): asymmetry parameter
!!\param aerodp        (IMAX,NSPC+1), vertically integrated aer-opt-depth
!!\section gel_go_aer_pro General Algorithm
!-----------------------------------
      subroutine aer_property_gocart                                    &
!...................................

!  ---  inputs:
     &     ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,aerfld,            &
     &       alon,alat,slmsk, laersw,laerlw,con_rd,                     &
     &       imax,nlay,nlp1,                                            &
!  ---  outputs:
     &       aerosw,aerolw,aerodp,ext550,errflg,errmsg                  &
     &     )

!  ==================================================================  !
!                                                                      !
!  aer_property_gocart maps prescribed gocart aerosol data set onto    !
!  model grids, and compute aerosol optical properties for sw and      !
!  lw radiations.                                                      !
!                                                                      !
!  inputs:                                                             !
!     prsi    - pressure at interface              mb      IMAX*NLP1   !
!     prsl    - layer mean pressure         (not used)     IMAX*NLAY   !
!     prslk   - exner function=(p/p0)**rocp (not used)     IMAX*NLAY   !
!     tvly    - layer virtual temperature   (not used)     IMAX*NLAY   !
!     rhlay   - layer mean relative humidity               IMAX*NLAY   !
!     dz      - layer thickness                    m       IMAX*NLAY   !
!     hz      - level high                         m       IMAX*NLP1   !
!     tracer  - aer tracer concentrations   (not used)  IMAX*NLAY*NTRAC!
!     aerfld  - prescribed aer tracer mixing ratios   IMAX*NLAY*NTRCAER!
!     alon, alat                                             IMAX      !
!             - longitude and latitude of given points in degree       !
!     slmsk   - sea/land mask (sea:0,land:1,sea-ice:2)       IMAX      !
!     laersw,laerlw                                             1      !
!             - logical flag for sw/lw aerosol calculations            !
!     IMAX    - horizontal dimension of arrays                  1      !
!     NLAY,NLP1-vertical dimensions of arrays                   1      !
!     con_rd  - Physical constant (gas constant for dry air)           !
!                                                                      !
!  outputs:                                                            !
!     aerosw - aeros opt properties for sw      IMAX*NLAY*NBDSW*NF_AESW!
!               (:,:,:,1): optical depth                               !
!               (:,:,:,2): single scattering albedo                    !
!               (:,:,:,3): asymmetry parameter                         !
!     aerolw - aeros opt properties for lw      IMAX*NLAY*NBDLW*NF_AELW!
!               (:,:,:,1): optical depth                               !
!               (:,:,:,2): single scattering albedo                    !
!               (:,:,:,3): asymmetry parameter                         !
!     aerodp - vertically integrated aer-opt-depth         IMAX*NSPC+1 !
!     errflg  - CCPP error flag                                        !
!     errmsg  - CCPP error message                                     !
!                                                                      !
!  module parameters and constants:                                    !
!     NSWBND  - total number of actual sw spectral bands computed      !
!     NLWBND  - total number of actual lw spectral bands computed      !
!     NSWLWBD - total number of sw+lw bands computed                   !
!                                                                      !
!  module variable: (set by subroutine aer_init)                       !
!                                                                      !
!  usage:    call aer_property_gocart                                  !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
      integer, intent(in) :: IMAX, NLAY, NLP1
      logical, intent(in) :: laersw, laerlw
      real (kind=kind_phys), intent(in) :: con_rd
      real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl,  &
     &       prslk, tvly, rhlay, dz, hz
      real (kind=kind_phys), dimension(:),   intent(in) :: alon, alat,  &
     &       slmsk
      real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer
      real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld

!  ---  outputs:
      real (kind=kind_phys), dimension(:,:,:,:), intent(out) ::         &
     &       aerosw, aerolw
      real (kind=kind_phys), dimension(:,:)    , intent(out) :: aerodp
      real (kind=kind_phys), dimension(:,:)    , intent(out) :: ext550
      integer,          intent(out) :: errflg
      character(len=*), intent(out) :: errmsg

!  ---  locals:
      real (kind=kind_phys), dimension(nlay,nswlwbd):: tauae,ssaae,asyae
      real (kind=kind_phys), dimension(nlay,1):: tauae_550
      real (kind=kind_phys), dimension(nlay,nspc)        :: spcodp

      real (kind=kind_phys),dimension(nlay,kcm) ::  aerms
      real (kind=kind_phys),dimension(nlay) :: dz1, rh1
      real (kind=kind_phys) :: plv, tv, rho
      integer               :: i, m, m1, k

!
!===>  ...  begin here
!

! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      lab_do_IMAXg : do i = 1, IMAX

! --- initialize tauae, ssaae, asyae
        do m = 1, NSWLWBD
          do k = 1, NLAY
           tauae(k,m) = f_zero
           if (m==nv_aod) then
             tauae_550(k,1) = f_zero
           endif
           ssaae(k,m) = f_one
           asyae(k,m) = f_zero
          enddo
        enddo

! --- set floor value for aerms (kg/m3)
        do k = 1, NLAY
          do m = 1, kcm
            aerms(k,m) = 1.e-15
          enddo
        enddo

        do k = 1, NLAY
          do m = 1, nspc
            spcodp(k,m) = f_zero
          enddo
        enddo

        do k = 1, NLAY
          rh1(k) = rhlay(i,k)          !
          dz1(k) = 1000.*dz (i,k)      ! thickness converted from km to m
          plv = 100.*prsl(i,k)         ! convert pressure from mb to Pa
          tv =  tvly(i,k)              ! virtual temp in K
          rho = plv / ( con_rd * tv)   ! air density in kg/m3

          do m = 1, KCM
           aerms(k,m) = aerfld(i,k,m)*rho  ! dry mass (kg/m3)
          enddo
!
! --- calculate sw/lw aerosol optical properties for the
!     corresponding frequency bands

          call aeropt
!   ---   inputs:  (in-scope variables)
!   ---   outputs: (in-scope variables)

        enddo  ! end_do_k_loop

! ----------------------------------------------------------------------

! --- update aerosw and aerolw arrays
        if ( laersw ) then

          do m = 1, NBDSW
            do k = 1, NLAY
              aerosw(i,k,m,1) = tauae(k,m)
              aerosw(i,k,m,2) = ssaae(k,m)
              aerosw(i,k,m,3) = asyae(k,m)
            enddo
          enddo

! --- update diagnostic aod arrays
          do k = 1, NLAY
            aerodp(i,1) = aerodp(i,1) + tauae_550(k,1)
            ext550(i,k) = tauae_550(k,1)
            do m = 1, NSPC
              aerodp(i,m+1) = aerodp(i,m+1)+spcodp(k,m)
            enddo
          enddo

        endif     ! end if_larsw_block

        if ( laerlw ) then

          if ( NLWBND == 1 ) then
            m1 = NSWBND + 1
            do m = 1, NBDLW
              do k = 1, NLAY
                aerolw(i,k,m,1) = tauae(k,m1)
                aerolw(i,k,m,2) = ssaae(k,m1)
                aerolw(i,k,m,3) = asyae(k,m1)
              enddo
            enddo
          else
            do m = 1, NBDLW
              m1 = NSWBND + m
              do k = 1, NLAY
                aerolw(i,k,m,1) = tauae(k,m1)
                aerolw(i,k,m,2) = ssaae(k,m1)
                aerolw(i,k,m,3) = asyae(k,m1)
              enddo
            enddo
          endif

        endif     ! end if_laerlw_block

      enddo  lab_do_IMAXg

! =================
      contains
! =================

!--------------------------------
!> compute aerosols optical properties in NSWLWBD bands for gocart
!!  aerosol species
      subroutine aeropt
!................................

!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

!  ==================================================================  !
!                                                                      !
!  compute aerosols optical properties in NSWLWBD bands for gocart     !
!  aerosol species                                                     !
!                                                                      !
!  input variables:                                                    !
!     rh1    - relative humidity                     %     NLAY        !
!     dz1    - layer thickness                       m     NLAY        !
!     aerms  - aerosol mass concentration           kg/m3  NLAY*KCM    !
!     NLAY   - vertical dimensions                   -     1           !
!                                                                      !
!  output variables:                                                   !
!     tauae  - optical depth                         -     NLAY*NSWLWBD!
!     ssaae  - single scattering albedo              -     NLAY*NSWLWBD!
!     asyae  - asymmetry parameter                   -     NLAY*NSWLWBD!
!     aerodp - vertically integrated aer-opt-depth   -     IMAX*NSPC+1 !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
!  ---  outputs:

!  ---  locals:
      real (kind=kind_phys) :: drh0, drh1, rdrh
      real (kind=kind_phys) :: cm, ext01, ext01_550, sca01,asy01,ssa01
      real (kind=kind_phys) :: ext1, ext1_550, asy1, ssa1,sca1,tau_550
      real (kind=kind_phys) :: sum_tau,sum_asy,sum_ssa,tau,asy,ssa
      real (kind=kind_phys) :: sum_tau_550
      integer               :: ih1, ih2, nbin, ib, ntrc, ktrc

! --- linear interp coeffs for rh-dep species
        ih2 = 1
        do while ( rh1(k) > rhlev_grt(ih2) )
          ih2 = ih2 + 1
          if ( ih2 > krhlev ) exit
        enddo
        ih1 = max( 1, ih2-1 )
        ih2 = min( krhlev, ih2 )

        drh0 = rhlev_grt(ih2) - rhlev_grt(ih1)
        drh1 = rh1(k) - rhlev_grt(ih1)
        if ( ih1 == ih2 ) then
         rdrh = f_zero
        else
         rdrh = drh1 / drh0
        endif

! --- compute optical properties for each spectral bands
        do ib = 1, nswlwbd

           sum_tau = f_zero
           if (ib == nv_aod ) then
             sum_tau_550 = f_zero
             ext1_550 = f_zero
           endif
           sum_ssa = f_zero
           sum_asy = f_zero

! --- determine tau, ssa, asy for dust aerosols
           ext1 = f_zero
           asy1 = f_zero
           sca1 = f_zero
           ssa1 = f_zero
           asy = f_zero
           ssa = f_zero
           do m = 1, kcm1
            cm =  max(aerms(k,m),0.0) * dz1(k)
            ext1 = ext1 + cm*extrhi_grt(m,ib)
            if (ib == nv_aod) then
              ext1_550 = ext1_550 + cm*extrhi_grt_550(m,1)
            endif
            sca1 = sca1 + cm*scarhi_grt(m,ib)
            ssa1 = ssa1 + cm*extrhi_grt(m,ib) * ssarhi_grt(m,ib)
            asy1 = asy1 + cm*scarhi_grt(m,ib) * asyrhi_grt(m,ib)
           enddo    ! m-loop
           tau = ext1
           if (ext1 > f_zero) ssa=min(f_one, ssa1/ext1)
           if (sca1 > f_zero) asy=min(f_one, asy1/sca1)

! --- update aod from individual species
           if ( ib==nv_aod ) then
             tau_550 = ext1_550
             spcodp(k,1) =  tau_550
             sum_tau_550 = sum_tau_550 + tau_550
           endif
! --- update sum_tau, sum_ssa, sum_asy
           sum_tau = sum_tau + tau
           sum_ssa = sum_ssa + tau * ssa
           sum_asy = sum_asy + tau * ssa * asy

! --- determine tau, ssa, asy for non-dust aerosols
           do ntrc = 2, nspc
            ext1 = f_zero
            if ( ib==nv_aod ) then
              ext1_550 = f_zero
            endif
            asy1 = f_zero
            sca1 = f_zero
            ssa1 = f_zero
            ktrc = trc_to_aod(ntrc)
            do nbin = 1, num_radius(ntrc)
             m1 =  radius_lower(ntrc) + nbin - 1
             m =   m1 - num_radius(1)                   ! exclude dust aerosols
             cm =  max(aerms(k,m1),0.0) * dz1(k)
             ext01 = extrhd_grt(ih1,m,ib) +                             &
     &              rdrh * (extrhd_grt(ih2,m,ib)-extrhd_grt(ih1,m,ib))
             if ( ib==nv_aod ) then
               ext01_550 = extrhd_grt_550(ih1,m,1) +                    &
     &         rdrh * (extrhd_grt_550(ih2,m,1)-extrhd_grt_550(ih1,m,1))
             endif
             sca01 = scarhd_grt(ih1,m,ib) +                             &
     &              rdrh * (scarhd_grt(ih2,m,ib)-scarhd_grt(ih1,m,ib))
             ssa01 = ssarhd_grt(ih1,m,ib) +                             &
     &              rdrh * (ssarhd_grt(ih2,m,ib)-ssarhd_grt(ih1,m,ib))
             asy01 = asyrhd_grt(ih1,m,ib) +                             &
     &              rdrh * (asyrhd_grt(ih2,m,ib)-asyrhd_grt(ih1,m,ib))
             ext1 = ext1 + cm*ext01
             if ( ib==nv_aod ) then
               ext1_550 = ext1_550 + cm*ext01_550
             endif
             sca1 = sca1 + cm*sca01
             ssa1 = ssa1 + cm*ext01 * ssa01
             asy1 = asy1 + cm*sca01 * asy01
            enddo        ! end_do_nbin_loop
            tau = ext1
            if (ext1 > f_zero) ssa=min(f_one, ssa1/ext1)
            if (sca1 > f_zero) asy=min(f_one, asy1/sca1)
! --- update aod from individual species
            if ( ib==nv_aod ) then
              tau_550 = ext1_550
              spcodp(k,ktrc) = tau_550
              sum_tau_550 = sum_tau_550 + tau_550
            endif
! --- update sum_tau, sum_ssa, sum_asy
            sum_tau = sum_tau + tau
            sum_ssa = sum_ssa + tau * ssa
            sum_asy = sum_asy + tau * ssa * asy
           enddo        ! end_do_ntrc_loop

! --- determine total tau, ssa, asy for aerosol mixture
           tauae(k,ib) = sum_tau
           if ( ib==nv_aod ) then
             tauae_550(k,1) = sum_tau_550
           endif
           if (sum_tau > f_zero) ssaae(k,ib) = sum_ssa / sum_tau
           if (sum_ssa > f_zero) asyae(k,ib) = sum_asy / sum_ssa

        enddo         ! end_do_ib_loop
!
      return
!................................
      end subroutine aeropt
!--------------------------------

!...................................
      end subroutine aer_property_gocart
!-----------------------------------
!
! =======================================================================
!..........................................!
      end module module_radiation_aerosols !
!==========================================!