!>  \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:                                                     !
!           ( none )                                                   !
!                                                                      !
!      'aer_update' -- updating aerosol data                           !
!         inputs:                                                      !
!           ( iyear, imon, me )                                        !
!         outputs:                                                     !
!           ( none )                                                   !
!                                                                      !
!      'setaer'     -- mapping aeros profile, compute aeros opticals   !
!         inputs:                                                      !
!           (prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,xlon,xlat,        !
!            IMAX,NLAY,NLP1, lsswr,lslwr,                              !
!         outputs:                                                     !
!           (aerosw,aerolw,tau_gocart)                                 !
!!          (aerosw,aerolw,aerodp)                                     !
!                                                                      !
!                                                                      !
!   external modules referenced:                                       !
!       'module physparam'               in 'physparam.f'              !
!       'module physcons'                in 'physcons.f'               !
!       '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.    !
!                                                                      !
!   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                       !
!                                                                      !
!   references for stratosperic volcanical aerosols:                   !
!     sato et al. 1993 - jgr, v98, d12, 22987-22994                    !
!                                                                      !
!                                                                      !
!!!!!  ==========================================================  !!!!!
!!!!!                       end descriptions                       !!!!!
!!!!!  ==========================================================  !!!!!



!> \ingroup rad
!! \defgroup module_radiation_aerosols module_radiation_aerosols
!> @{
!!  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
!! - Stratospheric volcanical aerosols:
!! Sato et al. 1993 \cite sato_et_al_1993
!========================================!
      module module_radiation_aerosols   !
!........................................!
!
      use physparam,only : iaermdl, iaerflg, lalw1bd, aeros_file,       &
     &                     ivflip, kind_phys, kind_io4, kind_io8
      use physcons, only : con_pi, con_rd, con_g, con_t0c, con_c,       &
     &                     con_boltz, con_plnk, con_amd

      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 gfs_phy_tracer_config,        only : gfs_phy_tracer, trcindx
!
      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:
!   - KCM, KCM1, KCM2 are determined from subroutine 'set_aerspc'
!> num of bands for aer data (gocart)
      integer, parameter :: KAERBND=61
!> num of rh levels for rh-dep components
      integer, parameter :: KRHLEV =36
!*    integer, parameter :: KCM1 = 8   ! num of rh independent aer !species
!*    integer, parameter :: KCM2 = 5   ! num of rh dependent aer species
!*    integer, parameter :: KCM  = KCM1 + KCM2
!> num of rh indep aerosols (set in subr set_aerspc)
      integer, save      :: KCM1 = 0
!> num of rh dep aerosols (set in subr set_aerspc)
      integer, save      :: KCM2 = 0
!> =KCM1+KCM2 (set in subr set_aerspc)
      integer, save      :: KCM

      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 /

!  ---  the following arrays are allocate and setup in subr 'gocrt_aerinit'
!  ------  gocart aerosol specification    ------
!  =>  transported aerosol species:
!      DU (5-bins)
!      SS (4 bins for climo mode and 5 bins for fcst mode)
!      SU (dms, so2, so4, msa)
!      OC (phobic, philic) and BC (phobic, philic)
!  =>  species and lumped species for aerosol optical properties
!      DU (5-bins, with 4 sub-groups in the submicron bin )
!      SS (ssam for submicron, sscm for coarse mode)
!      SU (so4)
!      OC (phobic, philic) and BC (phobic, philic)
!  =>  specification used for aerosol optical properties luts
!      DU (8 bins)
!      SS (ssam, sscm)
!      SU (suso)
!      OC (waso) and BC (soot)
!
!   - spectral band structure:
!      iendwv_grt(KAERBND)      - ending wavenumber (cm**-1) for each band
!   - relative humidity independent aerosol optical properties:
!   ===> species : dust (8 bins)
!      rhidext0_grt(KAERBND,KCM1) - extinction coefficient
!      rhidssa0_grt(KAERBND,KCM1) - single scattering albedo
!      rhidasy0_grt(KAERBND,KCM1) - asymmetry parameter
!   - relative humidity dependent aerosol optical properties:
!   ===> species : soot, suso, waso, ssam, sscm
!      rhdpext0_grt(KAERBND,KRHLEV,KCM2) - extinction coefficient
!      rhdpssa0_grt(KAERBND,KRHLEV,KCM2) - single scattering albedo
!      rhdpasy0_grt(KAERBND,KRHLEV,KCM2) - asymmetry parameter

!> spectral band structure: ending wavenumber (\f$cm^-1\f$) for each band
      integer,               allocatable, dimension(:) :: iendwv_grt
! relative humidity independent aerosol optical properties:
!! species : dust (8 bins)

!> \name relative humidity independent aerosol optical properties:
!! species : dust (8 bins)

!> extinction coefficient
      real (kind=kind_phys),allocatable, dimension(:,:)  :: rhidext0_grt
!> single scattering albedo
      real (kind=kind_phys),allocatable, dimension(:,:)  :: rhidssa0_grt
!> asymmetry parameter
      real (kind=kind_phys), allocatable, dimension(:,:) :: rhidasy0_grt
!
! relative humidity dependent aerosol optical properties:
! species : soot, suso, waso, ssam, sscm

!> \name relative humidity dependent aerosol optical properties:
!! species : soot, suso, waso, ssam, sscm

!> extinction coefficient
      real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpext0_grt
!> single scattering albedo
      real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpssa0_grt
!> asymmetry parameter
      real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpasy0_grt

!   - relative humidity independent aerosol optical properties:
!      extrhi_grt(KCM1,NSWLWBD) - extinction coefficient for sw+lw spectral band
!      ssarhi_grt(KCM1,NSWLWBD) - single scattering albedo for sw+lw spectral band
!      asyrhi_grt(KCM1,NSWLWBD) - asymmetry parameter for sw+lw spectral band
!   - relative humidity dependent aerosol optical properties:
!      extrhd_grt(KRHLEV,KCM2,NSWLWBD) - extinction 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

!>\name relative humidity independent aerosol optical properties

!> extinction coefficient for SW+LW spectral band
      real (kind=kind_phys),allocatable,save,dimension(:,:) ::          &
     &      extrhi_grt
!> single scattering albedo for SW+LW spectral band
      real (kind=kind_phys),allocatable,save,dimension(:,:) ::          &
     &      ssarhi_grt
!> asymmetry parameter for SW+LW spectral band
      real (kind=kind_phys),allocatable,save,dimension(:,:) ::          &
     &      asyrhi_grt

!> \name relative humidity dependent aerosol optical properties

!> extinction coefficient for SW+LW spectral band
      real (kind=kind_phys),allocatable,save,dimension(:,:,:) ::        &
     &      extrhd_grt
!> single scattering albedo for SW+LW band
      real (kind=kind_phys),allocatable,save,dimension(:,:,:) ::        &
     &      ssarhd_grt
!> asymmetry parameter for SW+LW band
      real (kind=kind_phys),allocatable,save,dimension(:,:,:) ::        &
     &      asyrhd_grt

!> \name module variables for gocart aerosol clim data set

! --------------------------------------------------------------------- !
!   section-5 : module variables for gocart aerosol climo data set      !
! --------------------------------------------------------------------- !
!     This version only supports geos3-gocart data set (Jan 2010)
!     Modified to support geos4-gocart data set        (May 2010)
!
!  geos3-gocart vs geos4-gocart
!  (1) Use the same module variables
!      IMXG,JMXG,KMXG,NMXG,psclmg,dmclmg,geos_rlon,geos_rlat
!  (2) Similarity between geos3 and geos 4:
!      identical lat/lon grids and aerosol specification;
!      direction of vertical index is bottom-up (sfc to toa)
!  (3) Difference between geos3 and geos4
!      vertical coordinate (sigma for geos3/hybrid_sigma_pressure for geos4)
!      aerosol units (mass concentration for geos3/mixing ratio for geos4)

!> num of lon-points in geos dataset
      integer, parameter :: IMXG = 144 
!> num of lat-points in geos dataset
      integer, parameter :: JMXG = 91 
!> num of vertical layers in geos dataset
      integer, parameter :: KMXG = 30  
!*    integer, parameter :: NMXG = 12 
!> to be determined by set_aerspc
      integer, save      :: NMXG    

      real (kind=kind_phys), parameter :: dltx = 360.0 / float(IMXG)
      real (kind=kind_phys), parameter :: dlty = 180.0 / float(JMXG-1)

!  --- the following arrays are allocated and setup in 'rd_gocart_clim'
!   - geos-gocart climo data (input dataset)
!     psclmg  - pressure in cb                   IMXG*JMXG*KMXG
!     dmclmg  - aerosol dry mass in g/m3         IMXG*JMXG*KMXG*NMXG
!               or aerosol mixing ratio in mol/mol or Kg/Kg

!> pressure in cb
      real (kind=kind_phys),allocatable, save:: psclmg(:,:,:)
!> aerosol dry mass in g/m3 or aerosol mixing ration in mol/mol or Kg/Kg
      real (kind=kind_phys),allocatable, save:: dmclmg(:,:,:,:)

!   - geos-gocart lat/lon arrays
!> geos-gocart longitude arrays
      real (kind=kind_phys), allocatable, save, dimension(:):: geos_rlon
!> geos-gocart latitude arrays
      real (kind=kind_phys), allocatable, save, dimension(:):: geos_rlat

!> control flag for gocart climo data set: xxxx as default; ver3 for geos3;
!! ver4 for geos4; 0000 for unknown data
      character*4, save  :: gocart_climo = 'xxxx'

!> molecular wght of gocart aerosol species
      real (kind=kind_io4), allocatable :: molwgt(:)

! ---------------------------------------------------------------------
! !
!   section-6 : module variables for gocart aerosol scheme options
!   !
! ---------------------------------------------------------------------
! !

!> logical parameter for gocart initialization control
      logical, save :: lgrtint = .true.

!> logical parameter for gocart debug print control
!     logical, save :: lckprnt = .true.
      logical, save :: lckprnt = .false.

!  --- the following index/flag/weight are set up in 'set_aerspc'

!> merging coefficients for fcst/clim; determined from fdaer
      real (kind=kind_phys), save :: ctaer = f_zero   ! user specified wgt

!> option to get fcst gocart aerosol field
      logical, save :: get_fcst = .true.   
!> option to get clim gocart aerosol field
      logical, save :: get_clim = .true.  

!  ------  gocart aerosol specification    ------
!  =>  transported aerosol species:
!      DU (5-bins)
!      SS (4 bins for climo mode and 5 bins for fcst mode)
!      SU (dms, so2, so4, msa)
!      OC (phobic, philic) and BC (phobic, philic)
!  =>  species and lumped species for aerosol optical properties
!      DU (5-bins, with 4 sub-groups in the submicron bin )
!      SS (ssam for submicron, sscm for coarse mode)
!      SU (so4)
!      OC (phobic, philic) and BC (phobic, philic)
!  =>  specification used for aerosol optical properties luts
!      DU (8 bins)
!      SS (ssam, sscm)
!      SU (suso)
!      OC (waso) and BC (soot)
!

!> index for rh dependent aerosol optical properties (2nd 
!! dimension for extrhd_grt, ssarhd_grt, and asyrhd_grt)
      integer, save :: isoot, iwaso, isuso, issam, isscm

!   - index for rh independent aerosol optical properties (1st
!     dimension for extrhi_grt, ssarhi_grt, and asyrhi_grt) is
!     not needed ===> hardwired to 8-bin dust

!   - index for gocart aerosol species to be included in the
!     calculations of aerosol optical properties (ext, ssa, asy)
!>  index for gocart aerosol species to be included in the
!!  calculations of aerosol optical properties (ext, ssa, asy)
      type  gocart_index_type
! dust
         integer :: dust1, dust2, dust3, dust4, dust5
! sea salt
         integer :: ssam,  sscm
! sulfate
         integer :: suso
! oc
         integer :: waso_phobic, waso_philic
! bc
         integer :: soot_phobic, soot_philic
      endtype
      type (gocart_index_type), save :: dm_indx

!>  index for gocart aerosols from prognostic tracer fields
      type  tracer_index_type
! dust
         integer :: du001, du002, du003, du004, du005
! sea salt
         integer :: ss001, ss002, ss003, ss004, ss005
! sulfate
         integer :: so4
! oc
         integer :: ocphobic, ocphilic
! bc
         integer :: bcphobic, bcphilic
      endtype
      type (tracer_index_type), save :: dmfcs_indx

!   - grid components to be included in the aeropt calculations
!> number of aerosol grid components
      integer, save                  :: num_gridcomp = 0 
!> aerosol grid components
      character, allocatable , save  :: gridcomp(:)*2   

!> default full-package setting
      integer, parameter          :: max_num_gridcomp = 5
!> data max_gridcomp  /'DU', 'BC', 'OC', 'SU', 'SS'/
      character*2                 :: max_gridcomp(max_num_gridcomp)
      data max_gridcomp  /'DU', 'BC', 'OC', 'SU', 'SS'/

! GOCART code modification end here (Sarah Lu)
! ------------------------!
! =======================================================================

!! ---  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

!  ---  public interfaces

      public aer_init, aer_update, setaer


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

!> 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 ) !  ---  inputs
!  ---  outputs: ( to module variables )

!  ==================================================================  !
!                                                                      !
!  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                             !
!                                                                      !
!  outputs: (to module variables)                                      !
!                                                                      !
!  external module variables: (in physparam)                           !
!     iaermdl - tropospheric aerosol model scheme flag                 !
!               =0 opac-clim; =1 gocart-clim, =2 gocart-prognostic     !
!               =5 opac-clim new spectral mapping                      !
!     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                     !
!     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                    !
!                                                                      !
!  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, gcrt_aerinit,                    !
!                       wrt_aerlog, set_volcaer, set_spectrum,         !
!                                                                      !
!  ==================================================================  !

!  ---  inputs:
      integer,  intent(in) :: NLAY, me

!  ---  output: ( none )

!  ---  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
!
      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= (iaerflg >= 100)         ! control flag for stratospheric volcanic aeros

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

      if ( me == 0 ) then

        call wrt_aerlog      ! write aerosol param info to log file
!  ---  inputs:   (in scope variables)
!  ---  outputs:  ( none )

      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
!  ---  inputs:   (module constants)
!  ---  outputs:  (in-scope variables)

!> -# 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                                         &
!  ---  outputs:
     &     )

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

!         call gcrt_climinit

!       elseif ( iaermdl == 2 ) then                ! gocart-prognostic scheme

!         call gcrt_aerinit

        else
          if ( me == 0 ) then
            print *,'  !!! ERROR in aerosol model scheme selection',    &
     &              ' iaermdl =',iaermdl
            stop
          endif
        endif

      endif    ! end if_iaerflg_block

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

      if ( lavoflg ) then

        call set_volcaer
!  ---  inputs:  (module variables)
!  ---  outputs: (module variables)

      endif    ! end if_lavoflg_block


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

!> This subroutine writes aerosol parameter configuration to run log file.
!--------------------------------
      subroutine wrt_aerlog
!................................
!  ---  inputs:    (in scope variables)
!  ---  outputs:   ( none )

!  ==================================================================  !
!                                                                      !
!  subprogram : wrt_aerlog                                             !
!                                                                      !
!    write aerosol parameter configuration to run log file.            !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
!  external module variables:  (in physparam)                          !
!   iaermdl  - aerosol scheme flag: 0:opac-clm; 1:gocart-clim;         !
!              2:gocart-prog; 5:opac-clim+new mapping                  !
!   iaerflg  - aerosol effect control flag: 3-digits (volc,lw,sw)      !
!   lalwflg  - toposphere lw aerosol effect: =f:no; =t:yes             !
!   laswflg  - toposphere sw aerosol effect: =f:no; =t:yes             !
!   lavoflg  - stratospherer volcanic aeros effect: =f:no; =t:yes      !
!                                                                      !
!  outputs: ( none )                                                   !
!                                                                      !
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call wrt_aerlog                                           !
!                                                                      !
!  ==================================================================  !

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

!
!===>  ...  begin here
!
      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
        stop
      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.
!>\section gel_set_spec General Algorithm
!--------------------------------
      subroutine set_spectrum
!................................
!  ---  inputs:   (module constants)
!  ---  outputs:  (in-scope variables)

!  ==================================================================  !
!                                                                      !
!  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:  (module constants)                                 
!!  -   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$)
!                                                                      !
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call set_spectrum                                         !
!                                                                      !
!  ==================================================================  !

!  ---  inputs: (module constants)
!     integer :: NWVTOT, NWVTIR

!  ---  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

!  ---  locals:
      real (kind=kind_phys) :: soltot, tmp1, tmp2, tmp3

      integer :: nb, nw, nw1, nw2, nmax, nmin
!
!===>  ...  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
!.............................
!  ---  inputs:   ( none )
!  ---  outputs:  (module variables)

!  ==================================================================  !
!                                                                      !
!  subprogram : set_volcaer                                            !
!                                                                      !
!    this is the initialization progrmam for stratospheric volcanic    !
!    aerosols.                                                         !
!                                                                      !
!  subroutines called: none                                            !
!                                                                      !
!  usage:    call set_volcaer                                          !
!                                                                      !
!  ==================================================================  !

!  ---  inputs: (none)

!  ---  output: (module variables)
!     integer :: ivolae(:,:,:)

!  ---  locals:
!
!===>  ...  begin here
!
!  ---  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                                         &          ! ---  inputs
     &     )                                                           !  ---  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                      !
!                                                                      !
!  outputs: (to module variables)                                      !
!                                                                      !
!  external module variables: (in physparam)                           !
!     iaerflg - abc 3-digit integer aerosol flag (abc:volc,lw,sw)      !
!               a: =0 use background stratospheric aerosol             !
!                  =1 incl stratospheric vocanic aeros (MINVYR-MAXVYR) !
!               b: =0 no topospheric aerosol in lw radiation           !
!                  =1 include tropspheric aerosols for lw radiation    !
!               c: =0 no topospheric aerosol in sw radiation           !
!                  =1 include tropspheric aerosols for sw radiation    !
!     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                     !
!     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                    !
!                                                                      !
!  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

!  ---  output: ( none )

!  ---  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
!
!  --- ...  invoke tropospheric aerosol initialization

!> - call set_aercoef() to invoke tropospheric aerosol initialization.
      call set_aercoef
!  ---  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
!................................
!  ---  inputs:   (in-scope variables, module constants)
!  ---  outputs:  (module variables)

!  ==================================================================  !
!                                                                      !
!  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)                                  !
!                                                                      !
!  external module variables:  (in physparam)                          !
!   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 )
!  ---  output: ( none )

!  ---  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
!
!> -# 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 !!'
        stop
      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:  (in physparam)                           !
!   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 ) !  ---  inputs:
!  ---  outputs: ( to module variables )

!  ==================================================================  !
!                                                                      !
!  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                !
!                                                                      !
!  outputs: ( none )                                                   !
!                                                                      !
!  external module variables: (in physparam)                           !
!     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

!  ---  output: ( none )

!  ---  locals: ( none )
!
!===> ...  begin here
!
      if ( imon < 1 .or. imon > 12 ) then
        print *,' ***** ERROR in specifying requested month !!! ',      &
     &          'imon=', imon
        print *,' ***** STOPPED in subroutinte aer_update !!!'
        stop
      endif

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

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


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

!> This subroutine updates the monthly global distribution of aerosol
!! profiles in five degree horizontal resolution.
!--------------------------------
      subroutine trop_update
!................................
!  ---  inputs:    (in scope variables, module variables)
!  ---  outputs:   (module variables)

!  ==================================================================  !
!                                                                      !
!  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                     !
!                                                                      !
!  outputs: (module variables)                                         !
!                                                                      !
!  external module variables: (in physparam)                           !
!    aeros_file   - external aerosol data file name                    !
!                                                                      !
!  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: ( none )
!  ---  output: ( none )

!  ---  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
!
!  --- ...  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 !!'
        stop
      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
!................................
!  ---  inputs:    (in scope variables, module variables)
!  ---  outputs:   (module variables)

!  ==================================================================  !
!                                                                      !
!  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     !
!                                                                      !
!  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

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

      character :: cline*80, volcano_file*32
      data volcano_file / 'volcanic_aerosols_1850-1859.txt ' /
!
!===>  ...  begin here
!
      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 !!'
            stop
          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,xlon,xlat,         &   !  ---  inputs
     &       IMAX,NLAY,NLP1, lsswr,lslwr,                               &
     &       aerosw,aerolw                                              &   !  ---  outputs
     &,      aerodp                                                     &
     &     )

!  ==================================================================  !
!                                                                      !
!  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 !
!     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      !
!                                                                      !
!  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  !
!                                                                      !
!  external module variable: (in physparam)                            !
!     iaerflg - aerosol effect control flag (volc,lw,sw, 3-dig)        !
!     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   !
!     ivflip  - control flag for direction of vertical index           !
!               =0: index from toa to surface                          !
!               =1: index from surface to toa                          !
!                                                                      !
!  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

      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

      logical, intent(in) :: lsswr, lslwr


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

      real (kind=kind_phys), dimension(:,:)    , intent(out) :: aerodp

!  ---  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), parameter :: rdg  = 180.0 / con_pi
      real (kind=kind_phys), parameter :: rovg = 0.001 * con_rd / con_g

!===>  ...  begin here

      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


      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 (ivflip == 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).

!SARAH
!         if ( iaerflg == 1 ) then      ! use opac aerosol climatology
          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,                                            &
!    &         IMAX,NLAY,NLP1,NSPC1,                                      &
!  ---  outputs:
     &         aerosw,aerolw,aerodp                                       &
     &       )

!  ---  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
! SARAH
!        elseif ( iaerflg == 2 )   then    ! use gocart aerosol scheme
         elseif ( iaermdl == 1 )   then    ! use gocart aerosol scheme

          call setgocartaer                                               &

!  ---  inputs:
     &       ( alon,alat,prslk,rhlay,dz,hz,NSWLWBD,                       &
     &         prsl,tvly,tracer,                                          &
     &         IMAX,NLAY,NLP1, ivflip, lsswr,lslwr,                       &
!  ---  outputs:
     &         aerosw,aerolw                                              &
     &     )

        endif     ! end if_iaerflg_block

      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 ( ivflip == 0 ) 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_ivflip_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,                                            &   
     &       aerosw,aerolw,aerodp                                       &     !  ---  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      !
!                                                                      !
!  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 !
!                                                                      !
!  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                   !
!                                                                      !
!  external module variables: (in physparam)                           !
!     ivflip  - control flag for direction of vertical index           !
!               =0: index from toa to surface                          !
!               =1: index from surface to toa                          !
!                                                                      !
!  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

      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

!  ---  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
!
!> -# 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 = i1
        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
              stop
            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
              stop
            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 = j1
        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
              stop
            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
              stop
            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 (ivflip == 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)
            stop
          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
!  ---  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)
!         if ( lspcaod ) then
           do m = 1, NSPC
             aerodp(i,m+1) = spcodp(m)
           enddo
!         endif

        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
!................................

!  ---  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:
!  ---  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 ( ivflip == 0 ) 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
!-----------------------------------

!> @}
! =======================================================================
! GOCART code modification starts here (Sarah lu)  ---------------------!
!!
!! gocart_init : set_aerspc, rd_gocart_clim, rd_gocart_luts, optavg_grt
!! setgocartaer: aeropt_grt, map_aermr

!> The initialization program for gocart aerosols
!! - determine weight and index for aerosol composition/luts
!! - read in monthly global distribution of gocart aerosols
!! - read and map the tabulated aerosol optical spectral data onto 
!!   corresponding SW/LW radiation spectral bands.
!!
!>\param NWVTOT           total num of wave numbers used in sw spectrum   
!!\param solfwv           (NWVTOT), solar flux for each individual
!!                        wavenumber (w/m2)
!!\param soltot           total solar flux for the spectrual range (w/m2)
!!\param NWVTIR           total num of wave numbers used in the ir region 
!!\param eirfwv           (NWVTIR), ir flux(273k) for each individual
!!                        wavenumber (w/m2)
!!\param NBDSW            num of bands calculated for sw aeros opt prop   
!!\param NLWBND           num of bands calculated for lw aeros opt prop   
!!\param NSWLWBD          total num of bands calc for sw+lw aeros opt prop
!!\param imon             month of the year                               
!!\param me               print message control flag    
!!\param raddt
!!\param fdaer
!>\section gel_go_ini General Algorithm
!! @{
!-----------------------------------
      subroutine gocart_init                                            &     
     &     ( NWVTOT,solfwv,soltot,NWVTIR,eirfwv,                        &   !  ---  inputs:
     &       NBDSW,NLWBND,NSWLWBD,imon,me,raddt,fdaer                   &   !  ---  outputs: ( none )
     &     )

!  ==================================================================  !
!                                                                      !
!  subprogram : gocart_init                                            !
!                                                                      !
!    this is the initialization program for gocart aerosols            !
!                                                                      !
!    - determine weight and index for aerosol composition/luts         !
!    - read in monthly global distribution of gocart aerosols          !
!    - read and map the tabulated aerosol optical spectral data        !
!        onto corresponding sw/lw radiation spectral bands.            !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
!  inputs:                                                             !
!   NWVTOT           - total num of wave numbers used in sw spectrum   !
!   solfwv(NWVTOT)   - solar flux for each individual wavenumber (w/m2)!
!   soltot           - total solar flux for the spectrual range  (w/m2)!
!   NWVTIR           - total num of wave numbers used in the ir region !
!   eirfwv(NWVTIR)   - ir flux(273k) for each individual wavenum (w/m2)!
!   NBDSW            - num of bands calculated for sw aeros opt prop   !
!   NLWBND           - num of bands calculated for lw aeros opt prop   !
!   NSWLWBD          - total num of bands calc for sw+lw aeros opt prop!
!   imon             - month of the year                               !
!   me               - print message control flag                      !
!                                                                      !
!  outputs: (to the module variables)                                  !
!                                                                      !
!  module variables:                                                   !
!     NBDSW   - total number of sw spectral bands                      !
!     wvnum1,wvnum2 (NSWSTR:NSWEND)                                    !
!             - start/end wavenumbers for each of sw bands             !
!     NBDLW   - total number of lw spectral bands                      !
!     wvnlw1,wvnlw2 (NBDLW)                                            !
!             - start/end wavenumbers for each of lw bands             !
!     NSWLWBD - total number of sw+lw bands used in this version       !
!     extrhi_grt  - extinction coef for rh-indep aeros  KCM1*NSWLWBD   !
!     ssarhi_grt  - single-scat-alb for rh-indep aeros  KCM1*NSWLWBD   !
!     asyrhi_grt  - asymmetry factor for rh-indep aeros KCM1*NSWLWBD   !
!     extrhd_grt  - extinction coef for rh-dep aeros KRHLEV*KCM2*NSWLWBD!
!     ssarhd_grt  - single-scat-alb for rh-dep aeros KRHLEV*KCM2*NSWLWBD!
!     asyrhd_grt  - asymmetry factor for rh-dep aerosKRHLEV*KCM2*NSWLWBD!
!     ctaer       - merging coefficients for fcst/clim fields          !
!     get_fcst    - option to get fcst aerosol fields                  !
!     get_clim    - option to get clim aerosol fields                  !
!     dm_indx  - index for aer spec to be included in aeropt calculations  !
!     dmfcs_indx  - index for prognostic aerosol fields                !
!     psclmg      - geos3/4-gocart pressure            IMXG*JMXG*KMXG  !
!     dmclmg      - geos3-gocart aerosol dry mass   IMXG*JMXG*KMXG*NMXG!
!                   or geos4-gocart aerosol mixing ratio               !
!                                                                      !
!  usage:    call gocart_init                                          !
!                                                                      !
!  subprograms called:  set_aerspc, rd_gocart_clim,                    !
!                       rd_gocart_luts, optavg_grt                     !
!                                                                      !
!  ==================================================================  !

      implicit none

!  ---  inputs:
      integer, intent(in) :: NWVTOT,NWVTIR,NBDSW,NLWBND,NSWLWBD,imon,me

      real (kind=kind_phys), intent(in) :: raddt, fdaer

      real (kind=kind_phys), intent(in) :: solfwv(:),soltot, eirfwv(:)

!  ---  output: ( none )

!  ---  locals:

      real (kind=kind_phys), dimension(NBDSW,KAERBND)  :: solwaer
      real (kind=kind_phys), dimension(NBDSW)          :: solbnd
      real (kind=kind_phys), dimension(NLWBND,KAERBND) :: eirwaer
      real (kind=kind_phys), dimension(NLWBND)         :: eirbnd
      real (kind=kind_phys) :: sumsol, sumir, fac, tmp, wvs, wve

      integer, dimension(NBDSW)  :: nv1, nv2
      integer, dimension(NLWBND) :: nr1, nr2

      integer :: i, mb, ib, ii, iw, iw1, iw2, ik, ibs, ibe

!===>  ...  begin here

!--------------------------------------------------------------------------
!  (1) determine aerosol specification index and merging coefficients
!--------------------------------------------------------------------------

      if ( .not. lgrtint ) then

!  --- ...  already done aerspc initialization, continue

        continue

      else

!  --- ...  set aerosol specification index and merging coefficients

        call set_aerspc(raddt,fdaer)
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

      endif  ! end if_lgrtinit_block

!
!--------------------------------------------------------------------------
!  (2) read gocart climatological data
!--------------------------------------------------------------------------

!  --- ...  read gocart climatological data, if needed

      if ( get_clim ) then

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

      endif

!
!--------------------------------------------------------------------------
!  (3) read and map the tabulated aerosol optical spectral data
!           onto corresponding radiation spectral bands
!--------------------------------------------------------------------------

      if ( .not. lgrtint ) then

!  --- ...  already done optical property interpolation, exit

        return

      else

!  --- ...  reset lgrtint

        lgrtint = .false.

!  --- ...  read tabulated aerosol optical input data
        call rd_gocart_luts
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

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

        solbnd (:)   = f_zero
        solwaer(:,:) = f_zero

        nv_aod = 1

        ibs = 1
        ibe = 1
        wvs = wvn_sw1(1)
        wve = wvn_sw1(1)
        do ib = 2, NBDSW
          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

        do ib = 1, NBDSW
          mb = ib + NSWSTR - 1
          ii = 1
          iw1 = nint(wvn_sw1(mb))
          iw2 = nint(wvn_sw2(mb))

          Lab_swdowhile : do while ( iw1 > iendwv_grt(ii) )
            if ( ii == KAERBND ) 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_grt(ii) ) then
              solwaer(ib,ii) = sumsol

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

          if ( iw2 /= iendwv_grt(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

          if((me==0) .and. lckprnt) print *,'RAD-nv1,nv2:',             &
     &        ib,iw1,iw2,nv1(ib),iendwv_grt(nv1(ib)),                   &
     &        nv2(ib),iendwv_grt(nv2(ib)),                              &
     &        10000./iw1, 10000./iw2
        enddo     ! end do_ib_block for sw

! --- check the spectral range for the nv_550 band
        if((me==0) .and. lckprnt) then
          mb = nv_aod + NSWSTR - 1
          iw1 = nint(wvn_sw1(mb))
          iw2 = nint(wvn_sw2(mb))
           print *,'RAD-nv_aod:',                                       &
     &       nv_aod, iw1, iw2, 10000./iw1, 10000./iw2
        endif
!
!  --- ...  compute ir flux weights and interval indices for mapping
!           spectral bands between lw radiation and aerosol data

        eirbnd (:)   = f_zero
        eirwaer(:,:) = f_zero

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

        do ib = 1, NLWBND
          ii = 1
          if ( NLWBND == 1 ) then
            iw1 = 400                   ! corresponding 25 mu
            iw2 = 2500                  ! corresponding 4  mu
          else
            iw1 = nint(wvn_lw1(ib))
            iw2 = nint(wvn_lw2(ib))
          endif

          Lab_lwdowhile : do while ( iw1 > iendwv_grt(ii) )
            if ( ii == KAERBND ) 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_grt(ii) ) then
              eirwaer(ib,ii) = sumir

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

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

          nr2(ib) = ii

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

          if(me==0 .and. lckprnt) print *,'RAD-nr1,nr2:',                &
     &        ib,iw1,iw2,nr1(ib),iendwv_grt(nr1(ib)),                    &
     &        nr2(ib),iendwv_grt(nr2(ib)),                               &
     &        10000./iw1, 10000./iw2
        enddo     ! end do_ib_block for lw

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

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

        if(me==0 .and. lckprnt) then
          print *, 'RAD -After optavg_grt, sw band info'
          do ib = 1, NBDSW
           mb = ib + NSWSTR - 1
           print *,'RAD -wvnsw1,wvnsw2: ',ib,wvn_sw1(mb),wvn_sw2(mb)
           print *,'RAD -lamda1,lamda2: ',ib,10000./wvn_sw1(mb),        &
     &                                   10000./wvn_sw2(mb)
           print *,'RAD -extrhi_grt:', extrhi_grt(:,ib)
!          do i = 1, KRHLEV
           do i = 1, KRHLEV, 10
             print *, 'RAD -extrhd_grt:',i,rhlev_grt(i),                &
     &                                extrhd_grt(i,:,ib)
           enddo
          enddo
          print *, 'RAD -After optavg_grt, lw band info'
          do ib = 1, NLWBND
           ii = NBDSW + ib
           print *,'RAD -wvnlw1,wvnlw2: ',ib,wvn_lw1(ib),wvn_lw2(ib)
           print *,'RAD -lamda1,lamda2: ',ib,10000./wvn_lw1(ib),        &
     &                                   10000./wvn_lw2(ib)
           print *,'RAD -extrhi_grt:', extrhi_grt(:,ii)
!          do i = 1, KRHLEV
           do i = 1, KRHLEV, 10
             print *, 'RAD -extrhd_grt:',i,rhlev_grt(i),                &
     &                                extrhd_grt(i,:,ii)
           enddo
          enddo
        endif

!  --- ...  dealoocate input data arrays no longer needed
        deallocate ( iendwv_grt   )
        if ( allocated(rhidext0_grt) ) then
          deallocate ( rhidext0_grt )
          deallocate ( rhidssa0_grt )
          deallocate ( rhidasy0_grt )
        endif
        if ( allocated(rhdpext0_grt) ) then
          deallocate ( rhdpext0_grt )
          deallocate ( rhdpssa0_grt )
          deallocate ( rhdpasy0_grt )
        endif

      endif  ! end if_lgrtinit_block

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

!> This subroutine determines merging coefficients ctaer; setup aerosol
!!  specification.
!-----------------------------
      subroutine set_aerspc(raddt,fdaer)
!.............................
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

! ==================================================================== !
!                                                                      !
! subprogram: set_aerspc                                               !
!                                                                      !
! determine merging coefficients ctaer;                                !
! set up aerosol specification: num_gridcomp, gridcomp, dm_indx,       !
!                       dmfcs_indx, isoot, iwaso, isuso, issam, isscm  !
!                                                                      !
! Aerosol optical properties (ext, ssa, asy) are determined from       !
! NMGX (<=12) aerosol species                                          !
! ==> DU: dust1 (4 sub-micron bins), dust2, dust3, dust4, dust5        !
!     BC: soot_phobic, soot_philic                                     !
!     OC: waso_phobic, waso_philic                                     !
!     SU: suso (=so4)                                                  !
!     SS: ssam (accumulation mode), sscm (coarse mode)                 !
!                                                                      !
! The current version only supports prognostic aerosols (from GOCART   !
! in-line calculations) and climo aerosols (from GEOS-GOCART runs)    !
!                                                                      !
!  ==================================================================  !
!
      implicit none

!  ---  inputs:
      real (kind=kind_phys), intent(in) :: raddt, fdaer
!  ---  output:

!  ---  local:
!     real (kind=kind_phys)     :: raddt
      integer                   :: i, indxr
      character*2               :: tp, gridcomp_tmp(max_num_gridcomp)

!! ===> determine ctaer (user specified weight for fcst fields)
!     raddt = min(fhswr,fhlwr) / 24.
      if( fdaer >= 99999. ) ctaer = f_one
      if((fdaer>0.).and.(fdaer<99999.)) ctaer=exp(-raddt/fdaer)

      if(me==0 .and. lckprnt) then
        print *, 'RAD -raddt, fdaer,ctaer: ', raddt, fdaer, ctaer
        if (ctaer == f_one ) then
          print *, 'LU -aerosol fields determined from fcst'
        elseif (ctaer == f_zero) then
          print *, 'LU -aerosol fields determined from clim'
        else
          print *, 'LU -aerosol fields determined from fcst/clim'
        endif
      endif

!! ===> determine get_fcst and get_clim
!!    if fcst is chosen (ctaer == f_one ), set get_clim to F
!!    if clim is chosen (ctaer == f_zero), set get_fcst to F
      if ( ctaer == f_one  )  get_clim = .false.
      if ( ctaer == f_zero )  get_fcst = .false.

!! ===> determine aerosol species to be included in the calculations
!!      of aerosol optical properties (ext, ssa, asy)

!*  If climo option is chosen, the aerosol composition is hardwired
!*  to full package. If not, the composition is determined from
!*  tracer_config on-the-fly (full package or subset)
      lab_if_fcst : if ( get_fcst ) then

!!      use tracer_config to determine num_gridcomp and gridcomp
        if ( gfs_phy_tracer%doing_GOCART )  then
         if ( gfs_phy_tracer%doing_DU )  then
            num_gridcomp  =  num_gridcomp  + 1
            gridcomp_tmp(num_gridcomp) = 'DU'
         endif
         if ( gfs_phy_tracer%doing_SU ) then
            num_gridcomp  =  num_gridcomp  + 1
            gridcomp_tmp(num_gridcomp) = 'SU'
         endif
         if ( gfs_phy_tracer%doing_SS ) then
            num_gridcomp  =  num_gridcomp  + 1
            gridcomp_tmp(num_gridcomp) = 'SS'
         endif
         if ( gfs_phy_tracer%doing_OC ) then
            num_gridcomp  =  num_gridcomp  + 1
            gridcomp_tmp(num_gridcomp) = 'OC'
         endif
         if ( gfs_phy_tracer%doing_BC ) then
            num_gridcomp  =  num_gridcomp  + 1
            gridcomp_tmp(num_gridcomp) = 'BC'
         endif
!
         if ( num_gridcomp > 0 ) then
           allocate ( gridcomp(num_gridcomp) )
           gridcomp(1:num_gridcomp) = gridcomp_tmp(1:num_gridcomp)
         else
           print *,'ERROR: prognostic aerosols not found,abort',me
           stop 1000
         endif

        else      ! gfs_phy_tracer%doing_GOCART=F

         print *,'ERROR: prognostic aerosols option off, abort',me
        stop 1001

        endif     ! end_if_gfs_phy_tracer%doing_GOCART_if_

      else lab_if_fcst

!!      set to full package (max_num_gridcomp and max_gridcomp)
        num_gridcomp = max_num_gridcomp
        allocate ( gridcomp(num_gridcomp) )
        gridcomp(1:num_gridcomp) = max_gridcomp(1:num_gridcomp)

      endif lab_if_fcst

!!
!! Aerosol specification is determined as such:
!! A. For radiation-aerosol feedback, the specification is based on the aeropt
!!    routine from Mian Chin and Hongbin Yu (hydrophobic and hydrophilic for
!!    OC/BC; submicron and supermicron for SS, 8-bins (with 4 subgroups for the
!!    the submicron bin) for DU, and SO4 for SU)
!! B. For transport, the specification is determined from GOCART in-line module
!! C. For LUTS, (waso, soot, ssam, sscm, suso, dust) is used, based on the
!!    the OPAC climo aerosol scheme (implemented by Yu-Tai Hou)

!!=== <A> determine dm_indx and NMXG
      indxr = 0
      dm_indx%waso_phobic = -999         ! OC
      dm_indx%soot_phobic = -999         ! BC
      dm_indx%ssam = -999                ! SS
      dm_indx%suso = -999                ! SU
      dm_indx%dust1 = -999               ! DU
      do i = 1, num_gridcomp
         tp = gridcomp(i)
         select case ( tp )
         case ( 'OC')    ! consider hydrophobic and hydrophilic
           dm_indx%waso_phobic = indxr + 1
           dm_indx%waso_philic = indxr + 2
           indxr = indxr + 2
         case ( 'BC')    ! consider hydrophobic and hydrophilic
           dm_indx%soot_phobic = indxr + 1
           dm_indx%soot_philic = indxr + 2
           indxr = indxr + 2
         case ( 'SS')    ! consider submicron and supermicron
           dm_indx%ssam = indxr + 1
           dm_indx%sscm = indxr + 2
           indxr = indxr + 2
         case ( 'SU')    ! consider SO4 only
           dm_indx%suso = indxr + 1
           indxr = indxr + 1
         case ( 'DU')    ! consider all 5 bins
           dm_indx%dust1 = indxr + 1
           dm_indx%dust2 = indxr + 2
           dm_indx%dust3 = indxr + 3
           dm_indx%dust4 = indxr + 4
           dm_indx%dust5 = indxr + 5
           indxr = indxr + 5
         case default
           print *,'ERROR: aerosol species not supported, abort',me
           stop 1002
         end select
      enddo
!!
      NMXG       = indxr      ! num of gocart aer spec for opt cal
!!

!!=== <B> determine dmfcs_indx
!!    SS: 5-bins are considered for transport while only two groups
!!        (accumulation/coarse modes) are considered for radiation
!!    DU: 5-bins are considered for transport while 8 bins (with the
!!        submicorn bin exptended to 4 bins) are considered for radiation
!!    SU: DMS, SO2, and MSA are not considered for radiation

      if ( get_fcst ) then
         if ( gfs_phy_tracer%doing_OC )  then
            dmfcs_indx%ocphobic = trcindx ('ocphobic', gfs_phy_tracer)
            dmfcs_indx%ocphilic = trcindx ('ocphilic', gfs_phy_tracer)
         endif
         if ( gfs_phy_tracer%doing_BC )  then
            dmfcs_indx%bcphobic = trcindx ('bcphobic', gfs_phy_tracer)
            dmfcs_indx%bcphilic = trcindx ('bcphilic', gfs_phy_tracer)
         endif
         if ( gfs_phy_tracer%doing_SS )  then
            dmfcs_indx%ss001 = trcindx ('ss001', gfs_phy_tracer)
            dmfcs_indx%ss002 = trcindx ('ss002', gfs_phy_tracer)
            dmfcs_indx%ss003 = trcindx ('ss003', gfs_phy_tracer)
            dmfcs_indx%ss004 = trcindx ('ss004', gfs_phy_tracer)
            dmfcs_indx%ss005 = trcindx ('ss005', gfs_phy_tracer)
         endif
         if ( gfs_phy_tracer%doing_SU )  then
            dmfcs_indx%so4 = trcindx ('so4', gfs_phy_tracer)
         endif
         if ( gfs_phy_tracer%doing_DU )  then
            dmfcs_indx%du001 = trcindx ('du001', gfs_phy_tracer)
            dmfcs_indx%du002 = trcindx ('du002', gfs_phy_tracer)
            dmfcs_indx%du003 = trcindx ('du003', gfs_phy_tracer)
            dmfcs_indx%du004 = trcindx ('du004', gfs_phy_tracer)
            dmfcs_indx%du005 = trcindx ('du005', gfs_phy_tracer)
         endif
      endif

!!
!!=== <C> determin KCM, KCM1, KCM2
!!    DU: submicron bin (dust1) contains 4 sub-groups (e.g., hardwire
!!        8 bins for aerosol optical properties luts)
!!    OC/BC: while hydrophobic aerosols are rh-independent, the luts
!!        for hydrophilic aerosols are used (e.g., use the coeff
!!        corresponding to rh=0)
!!
      indxr = 1
      isoot = -999
      iwaso = -999
      isuso = -999
      issam = -999
      isscm = -999
      do i = 1, num_gridcomp
         tp = gridcomp(i)
         if ( tp /= 'DU' ) then  !<--- non-dust aerosols
           select case ( tp )
           case ( 'OC ')
             iwaso = indxr
           case ( 'BC ')
             isoot = indxr
           case ( 'SU ')
             isuso = indxr
           case ( 'SS ')
             issam = indxr
             isscm = indxr + 1
           end select
           if ( tp /= 'SS' ) then
             indxr = indxr + 1
           else
             indxr = indxr + 2
           endif
         else                   !<--- dust aerosols
           KCM1 =  8            ! num of rh independent aer species
         endif
      enddo
      KCM2 = indxr - 1          ! num of rh dependent aer species
      KCM  = KCM1 + KCM2        ! total num of aer species

!!
!! check print starts here
      if( me == 0 .and. lckprnt) then
       print *, 'RAD -num_gridcomp:', num_gridcomp
       print *, 'RAD -gridcomp    :', gridcomp(:)
       print *, 'RAD -NMXG:',  NMXG
       print *, 'RAD -dm_indx ===> '
       print *, 'RAD -aerspc: dust1=', dm_indx%dust1
       print *, 'RAD -aerspc: dust2=', dm_indx%dust2
       print *, 'RAD -aerspc: dust3=', dm_indx%dust3
       print *, 'RAD -aerspc: dust4=', dm_indx%dust4
       print *, 'RAD -aerspc: dust5=', dm_indx%dust5
       print *, 'RAD -aerspc: ssam=',  dm_indx%ssam
       print *, 'RAD -aerspc: sscm=',  dm_indx%sscm
       print *, 'RAD -aerspc: suso=',  dm_indx%suso
       print *, 'RAD -aerspc: waso_phobic=',dm_indx%waso_phobic
       print *, 'RAD -aerspc: waso_philic=',dm_indx%waso_philic
       print *, 'RAD -aerspc: soot_phobic=',dm_indx%soot_phobic
       print *, 'RAD -aerspc: soot_philic=',dm_indx%soot_philic

       print *, 'RAD -KCM1 =', KCM1
       print *, 'RAD -KCM2 =', KCM2
       print *, 'RAD -KCM  =', KCM
       if ( KCM2 > 0 ) then
         print *, 'RAD -aerspc: issam=', issam
         print *, 'RAD -aerspc: isscm=', isscm
         print *, 'RAD -aerspc: isuso=', isuso
         print *, 'RAD -aerspc: iwaso=', iwaso
         print *, 'RAD -aerspc: isoot=', isoot
       endif

       if ( get_fcst ) then
         print *, 'RAD -dmfcs_indx ===> '
         print *, 'RAD -trc_du001=',dmfcs_indx%du001
         print *, 'RAD -trc_du002=',dmfcs_indx%du002
         print *, 'RAD -trc_du003=',dmfcs_indx%du003
         print *, 'RAD -trc_du004=',dmfcs_indx%du004
         print *, 'RAD -trc_du005=',dmfcs_indx%du005
         print *, 'RAD -trc_so4  =',dmfcs_indx%so4
         print *, 'RAD -trc_ocphobic=',dmfcs_indx%ocphobic
         print *, 'RAD -trc_ocphilic=',dmfcs_indx%ocphilic
         print *, 'RAD -trc_bcphobic=',dmfcs_indx%bcphobic
         print *, 'RAD -trc_bcphilic=',dmfcs_indx%bcphilic
         print *, 'RAD -trc_ss001=',dmfcs_indx%ss001
         print *, 'RAD -trc_ss002=',dmfcs_indx%ss002
         print *, 'RAD -trc_ss003=',dmfcs_indx%ss003
         print *, 'RAD -trc_ss004=',dmfcs_indx%ss004
         print *, 'RAD -trc_ss005=',dmfcs_indx%ss005
       endif
      endif
!! check print ends here

      return
!                                                                      !
      end subroutine set_aerspc

!-----------------------------------
!> This subroutine reads input gocart aerosol optical data from Mie
!! code calculations.
!-----------------------------
      subroutine rd_gocart_luts
!.............................
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

! ==================================================================== !
! subprogram: rd_gocart_luts                                           !
!   read input gocart aerosol optical data from Mie code calculations  !
!                                                                      !
! Remarks (Quanhua (Mark) Liu, JCSDA, June 2008)                       !
!  The LUT is for NCEP selected 61 wave numbers and 6 aerosols         !
!  (dust, soot, suso, waso, ssam, and sscm) and 36 aerosol effective   !
!  size in microns.                                                    !
!                                                                      !
!  The LUT is computed using Mie code with a logorithm size            !
!  distribution for each of 36 effective sizes. The standard deviation !
!  sigma of the size, and min/max size follows Chin et al. 2000        !
!  For each effective size, it corresponds a relative humidity value.  !
!                                                                      !
!  The LUT contains the density, sigma, relative humidity, mean mode   !
!  radius, effective size, mass extinction coefficient, single         !
!  scattering albedo, asymmetry factor, and phase function             !
!                                                                      !
!  ==================================================================  !
!
      implicit none

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

!  ---  locals:
      INTEGER, PARAMETER :: NP = 100, NP2 = 2*NP, nWave=100,            &
     &                      nAero=6, n_p=36
      INTEGER :: NW, NS, nH, n_bin
      real (kind=kind_io8), Dimension( NP2 ) :: Angle, Cos_Angle,       &
     &                                          Cos_Weight
      real (kind=kind_io8), Dimension(n_p,nAero) :: RH, rm, reff
      real (kind=kind_io8), Dimension(nWave,n_p,nAero) ::               & 
     &                      ext0, sca0, asy0
      real (kind=kind_io8), Dimension(NP2,n_p,nWave,nAero) :: ph0
      real (kind=kind_io8) :: wavelength(nWave), density(nAero),        &
     &                        sigma(nAero), wave,n_fac,PI,t1,s1,g1
      CHARACTER(len=80) :: AerosolName(nAero)
      INTEGER    :: i, j, k, l, ij

      character  :: aerosol_file*30
      logical    :: file_exist
      integer    :: indx_dust(8)          ! map 36 dust bins to gocart size bins

      data aerosol_file  /"NCEP_AEROSOL.bin"/
      data AerosolName/ ' Dust ', ' Soot ', ' SUSO ', ' WASO ',         &
     &                  ' SSAM ', ' SSCM '/

!! 8 dust bins
!!  1       2       3       4       5       6       7       8
!! .1-.18, .18-.3, .3-.6, 0.6-1.0, 1.0-1.8, 1.8-3, 3-6,  6-10  <-- def
!!  0.1399  0.2399  0.4499 0.8000 1.3994  2.3964 4.4964  7.9887 <-- reff
      data indx_dust/4, 8, 12, 18, 21, 24, 30, 36/

!     PI = acos(-1.d0)

! -- allocate aerosol optical data
      if ( .not. allocated( iendwv_grt ) ) then
        allocate ( iendwv_grt (KAERBND) )
      endif
      if (.not. allocated(rhidext0_grt) .and. KCM1 > 0 ) then
        allocate ( rhidext0_grt(KAERBND,KCM1))
        allocate ( rhidssa0_grt(KAERBND,KCM1))
        allocate ( rhidasy0_grt(KAERBND,KCM1))
      endif
      if (.not. allocated(rhdpext0_grt) .and. KCM2 > 0 ) then
        allocate ( rhdpext0_grt(KAERBND,KRHLEV,KCM2))
        allocate ( rhdpssa0_grt(KAERBND,KRHLEV,KCM2))
        allocate ( rhdpasy0_grt(KAERBND,KRHLEV,KCM2))
      endif

! -- read luts
      inquire (file = aerosol_file, exist = file_exist)

      if ( file_exist ) then
        if(me==0 .and. lckprnt) print *,'RAD -open :',aerosol_file
        close (NIAERCM)
        open (unit=NIAERCM,file=aerosol_file,status='OLD',              &
     &        action='read',form='UNFORMATTED')
      else
        print *,'    Requested aerosol data file "',aerosol_file,       &
     &          '" not found!', me
        print *,'    *** Stopped in subroutine RD_GOCART_LUTS !!'
        stop 1003
      endif              ! end if_file_exist_block

      READ(NIAERCM) (Cos_Angle(i),i=1,NP)
      READ(NIAERCM) (Cos_Weight(i),i=1,NP)
      READ(NIAERCM)
      READ(NIAERCM)
      READ(NIAERCM) NW,NS
      READ(NIAERCM)
      READ(NIAERCM) (wavelength(i),i=1,NW)

! --- check nAero and NW
      if (NW /= KAERBND) then
        print *, "Incorrect spectral band, abort ", NW
        stop 1004
      endif

! --- convert wavelength to wavenumber
      do i = 1, KAERBND
       iendwv_grt(i) = 10000. / wavelength(i)
       if(me==0 .and. lckprnt) print *,'RAD -wn,lamda:',                &
     &           i,iendwv_grt(i),wavelength(i)
      enddo

      DO j = 1, nAero
        if(me==0 .and. lckprnt) print *,'RAD -read LUTs:',              &
     &            j,AerosolName(j)
        READ(NIAERCM)
        READ(NIAERCM)
        READ(NIAERCM) n_bin, density(j), sigma(j)
        READ(NIAERCM)
        READ(NIAERCM) (RH(i,j),i=1, n_bin)
        READ(NIAERCM)
        READ(NIAERCM) (rm(i,j),i=1, n_bin)
        READ(NIAERCM)
        READ(NIAERCM) (reff(i,j),i=1, n_bin)

! --- check n_bin
        if (n_bin /= KRHLEV ) then
          print *, "Incorrect rh levels, abort ", n_bin
          stop 1005
        endif

! --- read luts
        DO k = 1, NW
          READ(NIAERCM) wave,(ext0(k,L,j),L=1,n_bin)
          READ(NIAERCM) (sca0(k,L,j),L=1,n_bin)
          READ(NIAERCM) (asy0(k,L,j),L=1,n_bin)
          READ(NIAERCM) (ph0(1:NP2,L,k,j),L=1,n_bin)
        END DO

! --- map luts input to module variables
        if (AerosolName(j) == ' Dust ' ) then
         if ( KCM1 > 0) then    !<-- only if rh independent aerosols are needed
          do i = 1, KCM1
           rhidext0_grt(1:KAERBND,i)=ext0(1:KAERBND,indx_dust(i),j)
           rhidssa0_grt(1:KAERBND,i)=sca0(1:KAERBND,indx_dust(i),j)
           rhidasy0_grt(1:KAERBND,i)=asy0(1:KAERBND,indx_dust(i),j)
          enddo
         endif
        else
         if ( KCM2 > 0) then    !<-- only if rh dependent aerosols are needed
          if (AerosolName(j) == ' Soot ') ij = isoot
          if (AerosolName(j) == ' SUSO ') ij = isuso
          if (AerosolName(j) == ' WASO ') ij = iwaso
          if (AerosolName(j) == ' SSAM ') ij = issam
          if (AerosolName(j) == ' SSCM ') ij = isscm
          if ( ij .ne. -999 ) then
             rhdpext0_grt(1:KAERBND,1:KRHLEV,ij) =                      &
     &               ext0(1:KAERBND,1:KRHLEV,j)
             rhdpssa0_grt(1:KAERBND,1:KRHLEV,ij) =                      &
     &               sca0(1:KAERBND,1:KRHLEV,j)
             rhdpasy0_grt(1:KAERBND,1:KRHLEV,ij) =                      &
     &               asy0(1:KAERBND,1:KRHLEV,j)
          endif   ! if_ij
         endif    ! if_KCM2
        endif
      END DO

      return
!...................................
      end subroutine rd_gocart_luts
!-----------------------------------
!                                                                      !
!> This subroutine computes mean aerosols optical properties over each
!! SW/LW 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_grt
!.............................
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

! ==================================================================== !
!                                                                      !
! subprogram: optavg_grt                                               !
!                                                                      !
!   compute mean aerosols optical properties over each sw/lw 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  ===================  !
!                                                                      !
! input arguments:                                                     !
!   nv1,nv2 (NBDSW)  - 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 (NBDSW,KAERBND)                                            !
!                    - solar flux weight over each sw radiation band   !
!                      vs each aerosol data spectral band              !
!   eirwaer (NLWBND,KAERBND)                                            !
!                    - ir flux weight over each lw radiation band      !
!                      vs each aerosol data spectral band              !
!   solbnd  (NBDSW)  - solar flux weight over each sw radiation band   !
!   eirbnd  (NLWBND) - ir flux weight over each lw radiation band      !
!   NBDSW            - total number of sw spectral bands               !
!   NLWBND           - total number of lw spectral bands               !
!   NSWLWBD          - total number of sw+lw spectral bands            !
!                                                                      !
! output arguments: (to module variables)                              !
!                                                                      !
!  ==================================================================  !
!
      implicit none

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

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

      integer :: ib, nb, ni, nh, nc
!
!===> ...  begin here

!  --- ...  allocate aerosol optical data
      if (.not. allocated(extrhd_grt) .and. KCM2 > 0 ) then
        allocate ( extrhd_grt(KRHLEV,KCM2,NSWLWBD) )
        allocate ( ssarhd_grt(KRHLEV,KCM2,NSWLWBD) )
        allocate ( asyrhd_grt(KRHLEV,KCM2,NSWLWBD) )
      endif
      if (.not. allocated(extrhi_grt) .and. KCM1 > 0 ) then
        allocate ( extrhi_grt(KCM1,NSWLWBD) )
        allocate ( ssarhi_grt(KCM1,NSWLWBD) )
        allocate ( asyrhi_grt(KCM1,NSWLWBD) )
      endif
!
!  --- ...  loop for each sw radiation spectral band

      do nb = 1, NBDSW
        rsolbd = f_one / solbnd(nb)

!  ---  for rh independent aerosol species

        lab_rhi: if (KCM1 >  0 ) then
        do nc = 1, KCM1
          sumk    = f_zero
          sumok   = f_zero
          sumokg  = f_zero
          sumreft = f_zero

          do ni = nv1(nb), nv2(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(nb,ni)

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

          refb = sumreft * rsolbd

          extrhi_grt(nc,nb) = sumk   * 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
        endif lab_rhi

!  ---  for rh dependent aerosols species

        lab_rhd: if (KCM2 > 0 ) then
        do nc = 1, KCM2
          do nh = 1, KRHLEV
            sumk    = 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)
              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
            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
        endif lab_rhd

      enddo   !  end do_nb_block for sw

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

      do nb = 1, NLWBND

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

!  ---  for rh independent aerosol species

        lab_rhi_lw: if (KCM1 > 0 ) then
        do nc = 1, KCM1
          sumk    = f_zero
          sumok   = f_zero
          sumokg  = f_zero
          sumreft = f_zero

          do ni = nr1(nb), nr2(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(nb,ni)

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

          refb = sumreft * rirbd

          extrhi_grt(nc,ib) = sumk   * 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
        endif lab_rhi_lw

!  ---  for rh dependent aerosols species

        lab_rhd_lw: if (KCM2 > 0 ) then
        do nc = 1, KCM2
          do nh = 1, KRHLEV
            sumk    = 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)
              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
            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
        endif lab_rhd_lw

      enddo   !  end do_nb_block for lw

!
      return
!................................
      end subroutine optavg_grt
!--------------------------------
!
!> This subroutine:
!! - 1. read in aerosol dry mass and surface pressure from GEOS3-GOCART
!! C3.1 2000 monthly dataset or aerosol mixing ratio and surface
!! pressure from GEOS4-GOCART 2000-2007 averaged monthly data set.
!! - 2. compute goes lat/lon array (for horizontal mapping)
!-----------------------------------
      subroutine rd_gocart_clim
!...................................
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

!  ==================================================================  !
!                                                                      !
! subprogram: rd_gocart_clim                                           !
!                                                                      !
!   1. read in aerosol dry mass and surface pressure from GEOS3-GOCART !
!      C3.1 2000 monthly data set                                      !
!      or aerosol mixing ratio and surface pressure from GEOS4-GOCART  !
!      2000-2007 averaged monthly data set                             !
!   2. compute goes lat/lon array (for horizontal mapping)             !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
! inputs arguments:                                                    !
!     imon    - month of the year                                      !
!     me      - print message control flag                             !
!                                                                      !
! outputs arguments: (to the module variables)                         !
!     psclmg   - pressure (sfc to toa)    cb   IMXG*JMXG*KMXG          !
!     dmclmg   - aerosol dry mass/mixing ratio     IMXG*JMXG*KMXG*NMXG !
!     geos_rlon - goes longitude          deg      IMXG                !
!     geos_rlat - goes latitude           deg      JMXG                !
!                                                                      !
!  usage:    call rd_gocart_clim                                       !
!                                                                      !
!  program history:                                                    !
!    05/18/2010  ---  Lu    Add the option to read GEOS4-GOCART climo  !
!  ==================================================================  !
!
      implicit none

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

!  ---  locals:
      integer, parameter :: MAXSPC = 5
      real (kind=kind_io4), parameter  :: PINT = 0.01
      real (kind=kind_io4), parameter  :: EPSQ = 0.0

      integer         :: i, j, k, numspci, ii
      integer         :: icmp, nrecl, nt1, nt2, nn(MAXSPC)
      character       :: ymd*6, yr*4, mn*2, tp*2,                       &
     &                   fname*30, fin*30, aerosol_file*40
      logical         :: file_exist

      real (kind=kind_io4), dimension(KMXG)             :: sig
      real (kind=kind_io4), dimension(IMXG,JMXG)        :: ps
      real (kind=kind_io4), dimension(IMXG,JMXG,KMXG)   :: temp
      real (kind=kind_io4), dimension(IMXG,JMXG,KMXG,MAXSPC):: buff
      real (kind=kind_phys)   :: pstmp

!     Add the following variables for GEOS4-GOCART
      real (kind=kind_io4), dimension(KMXG):: hyam, hybm
      real (kind=kind_io4)                 :: p0

      data yr /'2000'/           !!<=== use 2000 as the climo proxy

!* sigma_coordinate for GEOS3-GOCART
!* P(i,j,k) = PINT + SIG(k) * (PS(i,j) - PINT)
      data SIG  /                                                       &
     &     9.98547E-01,9.94147E-01,9.86350E-01,9.74300E-01,9.56950E-01, &
     &     9.33150E-01,9.01750E-01,8.61500E-01,8.11000E-01,7.50600E-01, &
     &     6.82900E-01,6.10850E-01,5.37050E-01,4.63900E-01,3.93650E-01, &
     &     3.28275E-01,2.69500E-01,2.18295E-01,1.74820E-01,1.38840E-01, &
     &     1.09790E-01,8.66900E-02,6.84150E-02,5.39800E-02,4.25750E-02, &
     &     3.35700E-02,2.39900E-02,1.36775E-02,5.01750E-03,5.30000E-04 /

!* hybrid_sigma_pressure_coordinate for GEOS4-GOCART
!* p(i,j,k) = a(k)*p0 + b(k)*ps(i,j)
      data hyam/                                                        &
     &   0, 0.0062694, 0.02377049, 0.05011813, 0.08278809, 0.1186361,   &
     &   0.1540329, 0.1836373, 0.2043698, 0.2167788, 0.221193,          &
     &   0.217729, 0.2062951, 0.1865887, 0.1615213, 0.1372958,          &
     &   0.1167039, 0.09920014, 0.08432171, 0.06656809, 0.04765031,     &
     &   0.03382346, 0.0237648, 0.01435208, 0.00659734, 0.002826232,    &
     &   0.001118959, 0.0004086494, 0.0001368611, 3.750308e-05/

      data hybm /                                                       &
     &   0.992555, 0.9642, 0.90556, 0.816375, 0.703815, 0.576585,       &
     &   0.44445, 0.324385, 0.226815, 0.149165, 0.089375,               &
     &   0.045865, 0.017485, 0.00348, 0, 0, 0, 0, 0,                    &
     &   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 /

      data p0 /1013.25 /

!===>  ...  begin here

! --- allocate and initialize gocart climatological data
      if ( .not. allocated (dmclmg) ) then
        allocate ( dmclmg(IMXG,JMXG,KMXG,NMXG) )
        allocate ( psclmg(IMXG,JMXG,KMXG) )
        allocate ( molwgt(NMXG) )
      endif

      dmclmg(:,:,:,:)  = f_zero
      psclmg(:,:,:)    = f_zero
      molwgt(:)        = f_zero

! --- allocate and initialize geos lat and lon arrays
      if ( .not. allocated ( geos_rlon  )) then
          allocate (geos_rlon(IMXG))
          allocate (geos_rlat(JMXG))
      endif

      geos_rlon(:) = f_zero
      geos_rlat(:) = f_zero

! --- compute geos lat and lon arrays
      do i = 1, IMXG
         geos_rlon(i)     = -180. + (i-1)* dltx
      end do
      do j = 2, JMXG-1
         geos_rlat(j)     = -90. + (j-1)* dlty
      end do
      geos_rlat(1)      = -89.5
      geos_rlat(JMXG)   =  89.5

! --- determine whether GEOS3 or GEOS4 data set is provided
      if ( gocart_climo == 'xxxx' ) then
        gocart_climo='0000'
! check geos3-gocart climo
        aerosol_file = '200001.PS.avg'
        inquire (file = aerosol_file, exist = file_exist)
        if ( file_exist ) gocart_climo='ver3'
! check geos4-gocart climo
        aerosol_file = 'gocart_climo_2000x2007_ps_01.bin'
        inquire (file = aerosol_file, exist = file_exist)
        if ( file_exist ) gocart_climo='ver4'
      endif
!
!
! --- read ps (sfc pressure) and compute 3d pressure field (psclmg)
!
      write(mn,'(i2.2)') imon
      ymd = yr//mn
      aerosol_file = 'null'
      if ( gocart_climo == 'ver3' ) then
        aerosol_file = ymd//'.PS.avg'
      elseif ( gocart_climo == 'ver4' ) then
        aerosol_file = 'gocart_climo_2000x2007_ps_'//mn//'.bin'
      endif
!
      inquire (file = aerosol_file, exist = file_exist)
      lab_if_ps : if ( file_exist ) then

       close(NIAERCM)
       if ( gocart_climo == 'ver3' ) then
        nrecl = 4 * (IMXG * JMXG)
        open(NIAERCM, file=trim(aerosol_file),                          &
     &       action='read',access='direct',recl=nrecl)
        read(NIAERCM, rec=1) ps
        do j = 1, JMXG
          do i = 1, IMXG
            do k = 1, KMXG
              pstmp = pint + sig(k) * (ps(i,j) - pint)
              psclmg(i,j,k) = 0.1 * pstmp       ! convert mb to cb
            enddo
          enddo
        enddo

       elseif ( gocart_climo == 'ver4' ) then
         open(NIAERCM, file=trim(aerosol_file),                         &
     &        action='read',status='old', form='unformatted')
         read(NIAERCM) ps(:,:)
         do j = 1, JMXG
           do i = 1, IMXG
             do k = 1, KMXG
              pstmp = hyam(k)*p0 + hybm(k)*ps(i,j)
              psclmg(i,j,k) = 0.1 * pstmp       ! convert mb to cb
            enddo
          enddo
         enddo

       endif     !  ---- end if_gocart_climo

      else lab_if_ps

        print *,' *** Requested aerosol data file "',                    &
     &          trim(aerosol_file),  '" not found!'
        print *,' *** Stopped in RD_GOCART_CLIM ! ', me
        stop 1006
      endif   lab_if_ps
!
! --- read aerosol dry mass (g/m3) or mixing ratios (mol/mol,kg/kg)
!
      lab_do_icmp : do icmp = 1, num_gridcomp

         tp = gridcomp(icmp)

!        determine aerosol_file
         aerosol_file = 'null'
         if ( gocart_climo == 'ver3' ) then
           if(tp == 'DU')   fname='.DU.STD.tv20.g.avg'
           if(tp == 'SS')   fname='.SS.STD.tv17.g.avg'
           if(tp == 'SU')   fname='.SU.STD.tv15.g.avg'
           if(tp == 'OC')   fname='.CC.STD.tv15.g.avg'
           if(tp == 'BC')   fname='.CC.STD.tv15.g.avg'
           aerosol_file=ymd//trim(fname)
         elseif ( gocart_climo == 'ver4' ) then
           fin = 'gocart_climo_2000x2007_'
           if(tp == 'DU')   fname=trim(fin)//'du_'
           if(tp == 'SS')   fname=trim(fin)//'ss_'
           if(tp == 'SU')   fname=trim(fin)//'su_'
           if(tp == 'OC')   fname=trim(fin)//'cc_'
           if(tp == 'BC')   fname=trim(fin)//'cc_'
           aerosol_file=trim(fname)//mn//'.bin'
         endif

         numspci = 4
         if(tp == 'DU')   numspci = 5
         inquire (file=trim(aerosol_file), exist = file_exist)
         lab_if_aer: if ( file_exist ) then
!
          close(NIAERCM)
          if ( gocart_climo == 'ver3' ) then
          nrecl = 4 * numspci * (IMXG * JMXG * KMXG + 3)
          open (NIAERCM, file=trim(aerosol_file),                       &
     &         action='read',access='direct', recl=nrecl)
          read(NIAERCM,rec=1)(nt1,nt2,nn(i),buff(:,:,:,i),i=1,numspci)

          elseif ( gocart_climo == 'ver4' ) then
           open (NIAERCM, file=trim(aerosol_file),                      &
     &           action='read',status='old', form='unformatted')
           do i = 1, numspci
             do k = 1, KMXG
              read(NIAERCM) temp(:,:,k)
              buff(:,:,k,i) =  temp(:,:,k)
             enddo
           enddo
          endif

!!===> fill dmclmg with working array buff
          select case ( tp )

! fill in DU from DU: du1, du2, du3, du4, du5
          case ('DU' )
           if ( dm_indx%dust1 /= -999) then
            do ii = 1, 5
             dmclmg(:,:,:,dm_indx%dust1+ii-1) = buff(:,:,:,ii)
            enddo
           else
            print *, 'ERROR: invalid DU index, abort! ',me
            stop 1007
           endif

! fill in BC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic
          case ('BC' )
           if ( dm_indx%soot_phobic /= -999) then
            dmclmg(:,:,:,dm_indx%soot_phobic)=buff(:,:,:,1)
            dmclmg(:,:,:,dm_indx%soot_philic)=buff(:,:,:,3)
            molwgt(dm_indx%soot_phobic) = 12.
            molwgt(dm_indx%soot_philic) = 12.
           else
            print *, 'ERROR: invalid BC index, abort! ',me
            stop 1008
           endif

! fill in SU from SU: dms, so2, so4, msa
          case ('SU' )
           if ( dm_indx%suso /= -999) then
            dmclmg(:,:,:,dm_indx%suso) = buff(:,:,:,3)
            molwgt(dm_indx%suso) = 96.
           else
            print *, 'ERROR: invalid SU index, abort! ',me
            stop 1009
           endif

! fill in OC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic
          case ('OC' )
           if ( dm_indx%waso_phobic /= -999) then
            dmclmg(:,:,:,dm_indx%waso_phobic) = 1.4*buff(:,:,:,2)
            dmclmg(:,:,:,dm_indx%waso_philic) = 1.4*buff(:,:,:,4)
            molwgt(dm_indx%waso_phobic) = 12.
            molwgt(dm_indx%waso_philic) = 12.
           else
            print *, 'ERROR: invalid OC index, abort! ',me
            stop 1010
           endif

! fill in SS from SS: ss1, ss2, ss3, ss4
          case ('SS' )
           if ( dm_indx%ssam /= -999) then
            dmclmg(:,:,:,dm_indx%ssam) = buff(:,:,:,1)
            dmclmg(:,:,:,dm_indx%sscm) = buff(:,:,:,2) +                &
     &                     buff(:,:,:,3)+buff(:,:,:,4)
           else
            print *, 'ERROR: invalid SS index, abort! ',me
            stop  1011
           endif

          case default

            print *, 'ERROR: invalid aerosol species, abort ',tp
            stop  1012

          end select

         else   lab_if_aer
          print *,' *** Requested aerosol data file "',aerosol_file,    &
     &            '" not found!'
          print *,' *** Stopped in RD_GOCART_CLIM ! ', me
          stop 1013
         endif  lab_if_aer

       enddo lab_do_icmp

      return
!...................................
      end subroutine rd_gocart_clim
!-----------------------------------
!
!...................................
      end subroutine gocart_init
!-----------------------------------
!! @}

!> This subroutine computes SW + LW aerosol optical properties for
!! gocart aerosol species (merged from fcst and clim fields).
!!
!>\param alon      IMAX, longitude of given points in degree
!!\param alat      IMAX, latitude of given points in degree       
!!\param prslk     (IMAX,NLAY), pressure in cb  
!!\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 NSWLWBD   total number of sw+ir bands for aeros opt prop  
!!\param prsl      (IMAX,NLAY), layer mean pressure in mb 
!!\param tvly      (IMAX,NLAY), layer mean virtual temperature in K  
!!\param trcly     (IMAX,NLAY,NTRAC), layer mean specific tracer in g/g
!!\param IMAX      horizontal dimension of arrays              
!!\param NLAY,NLP1 vertical dimensions of arrays                   
!!\param ivflip    control flag for direction of vertical index  
!!\n               =0: index from toa to surface                
!!\n               =1: index from surface to toa               
!!\param lsswr,lslwr       logical flag 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                
!>\section gen_setgo General Algorithm
!!@{
!-----------------------------------
      subroutine setgocartaer                                           &           
     &     ( alon,alat,prslk,rhlay,dz,hz,NSWLWBD,                       &       !  ---  inputs:
     &       prsl,tvly,trcly,                                           &
     &       IMAX,NLAY,NLP1, ivflip, lsswr,lslwr,                       &
     &       aerosw,aerolw                                              &       !  ---  outputs:
     &     )


!  ==================================================================  !
!                                                                      !
!  setgocartaer computes sw + lw aerosol optical properties for gocart !
!  aerosol species (merged from fcst and clim fields)                  !
!                                                                      !
!  inputs:                                                             !
!     alon, alat                                             IMAX      !
!             - longitude and latitude of given points in degree       !
!     prslk   - pressure                           cb   IMAX*NLAY      !
!     rhlay   - layer mean relative humidity            IMAX*NLAY      !
!     dz      - layer thickness                    m    IMAX*NLAY      !
!     hz      - level high                         m    IMAX*NLP1      !
!     NSWLWBD - total number of sw+ir bands for aeros opt prop  1      !
!     prsl    - layer mean pressure                mb   IMAX*NLAY      !
!     tvly    - layer mean virtual temperature     k   IMAX*NLAY      !
!     trcly   - layer mean specific tracer         g/g  IMAX*NLAY*NTRAC!
!     IMAX    - horizontal dimension of arrays                  1      !
!     NLAY,NLP1-vertical dimensions of arrays                   1      !
!     ivflip   - control flag for direction of vertical index   1      !
!               =0: index from toa to surface                          !
!               =1: index from surface to toa                          !
!     lsswr,lslwr                                                      !
!             - logical flag for sw/lw radiation calls          1      !
!                                                                      !
!  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!
!                                                                      !
!  module parameters and constants:                                    !
!     NBDSW   - total number of sw bands for aeros opt prop     1      !
!     NLWBND  - total number of ir bands for aeros opt prop     1      !
!                                                                      !
!  module variable: (set by subroutine gocart_init)                    !
!     dmclmg  - aerosols dry mass/mixing ratios   IMXG*JMXG*KMXG*NMXG  !
!     psclmg  - pressure                cb        IMXG*JMXG*KMXG       !
!                                                                      !
!  usage:    call setgocartaer                                         !
!                                                                      !
!  subprograms called:  map_aermr, aeropt_grt                          !
!                                                                      !
!  ==================================================================  !
!
      implicit none

!  ---  inputs:
      integer, intent(in) :: IMAX,NLAY,NLP1,ivflip,NSWLWBD
      logical, intent(in) :: lsswr, lslwr

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

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

!  ---  locals:
      real (kind=kind_phys), dimension(NLAY) :: rh1, dz1
      real (kind=kind_phys), dimension(NLAY,NSWLWBD)::tauae,ssaae,asyae
      real (kind=kind_phys), dimension(NLAY,max_num_gridcomp) ::        & 
     &                       tauae_gocart

      real (kind=kind_phys) :: tmp1, tmp2

      integer               :: i, i1, i2, j1, j2, k, m, m1, kp

!     prognostic aerosols on gfs grids
      real (kind=kind_phys), dimension(:,:,:),allocatable:: aermr,dmfcs

! aerosol (dry mass) on gfs grids/levels
      real (kind=kind_phys), dimension(:,:), allocatable ::             &
     &  dmanl,dmclm, dmclmx
      real (kind=kind_phys), dimension(KMXG)     :: pstmp, pkstr
      real (kind=kind_phys) :: ptop, psfc, tem, plv, tv, rho

!  ---  conversion constants
      real (kind=kind_phys), parameter :: hdltx = 0.5 * dltx
      real (kind=kind_phys), parameter :: hdlty = 0.5 * dlty

!===>  ...  begin here
!
      if ( .not. allocated(dmanl) ) then
        allocate ( dmclmx(KMXG,NMXG) )
        allocate ( dmanl(NLAY,NMXG) )
        allocate ( dmclm(NLAY,NMXG) )

        allocate ( aermr(IMAX,NLAY,NMXG) )
        allocate ( dmfcs(IMAX,NLAY,NMXG) )
      endif
!
!> -# Call map_aermr() to map input tracer array (trcly) to local
!!    tracer array (aermr).
      dmfcs(:,:,:) = f_zero
      lab_if_fcst : if ( get_fcst ) then

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

      endif       lab_if_fcst
!
!> -# Map geos-gocart climo (dmclmg) to gfs grids (dmclm).
      lab_do_IMAX : do i = 1, IMAX

        dmclm(:,:) = f_zero

        lab_if_clim : if ( get_clim ) then
!  ---  map grid in longitude direction
          i2 = 1
          j2 = 1
          tmp1 = alon(i)
          if (tmp1 > 180.) tmp1 = tmp1 - 360.0
          lab_do_IMXG : do i1 = 1, IMXG
            tmp2 = geos_rlon(i1)
            if (tmp2 > 180.) tmp2 = tmp2 - 360.0
            if (abs(tmp1-tmp2) <= hdltx) then
              i2 = i1
              exit lab_do_IMXG
            endif
          enddo  lab_do_IMXG

!  ---  map grid in latitude direction
          lab_do_JMXG : do j1 = 1, JMXG
            if (abs(alat(i)-geos_rlat(j1)) <= hdlty) then
              j2 = j1
              exit lab_do_JMXG
            endif
          enddo  lab_do_JMXG

!  ---  update local arrays pstmp and dmclmx
          pstmp(:)= psclmg(i2,j2,:)*1000.0      ! cb to Pa
          dmclmx(:,:) = dmclmg(i2,j2,:,:)

!  ---  map geos-gocart climo (dmclmx) to gfs level (dmclm)
          pkstr(:)=fpkap(pstmp(:))
          psfc = pkstr(1)                       ! pressure at sfc
          ptop = pkstr(KMXG)                    ! pressure at toa

!  ---  map grid in verical direction (follow how ozone is mapped
!       in radiation_gases routine)
          do k = 1, NLAY
           kp = k                              ! from sfc to toa
           if(ivflip==0) kp = NLAY - k + 1     ! from toa to sfc
           tmp1 = prslk(i,kp)

           do m1 = 1, KMXG - 1                 ! from sfc to toa
             if(tmp1 > pkstr(m1+1) .and. tmp1 <= pkstr(m1)) then
               tmp2 = f_one / (pkstr(m1)-pkstr(m1+1))
               tem = (pkstr(m1) - tmp1) * tmp2
               dmclm(kp,:) = tem * dmclmx(m1+1,:)+                      &
     &                   (f_one-tem) * dmclmx(m1,:)
             endif
           enddo

!*         if(tmp1 > psfc) dmclm(kp,:) = dmclmx(1,:)
!*         if(tmp1 < ptop) dmclm(kp,:) = dmclmx(KMXG,:)

          enddo
        endif    lab_if_clim
!
!  ---  compute fcst/clim merged aerosol loading (dmanl) and the
!       radiation optical properties (aerosw, aerolw)
!
        do k = 1, NLAY

!  ---  map global to local arrays (rh1 and dz1)
          rh1(k) = rhlay(i,k)
          dz1(k) = dz   (i,k)

!  ---  convert from mixing ratio to dry mass (g/m3)
          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
          if ( get_fcst ) then
            do m = 1,  NMXG            ! mixing ratio (g/g)
            dmfcs(i,k,m) = max(1000.*(rho*aermr(i,k,m)),f_zero)
            enddo     ! m_do_loop
          endif
          if ( get_clim .and. (gocart_climo == 'ver4') ) then
            do m = 1,  NMXG
             dmclm(k,m)=1000.*dmclm(k,m)*rho !mixing ratio (g/g)
             if ( molwgt(m) /= 0. ) then     !mixing ratio (mol/mol)
              dmclm(k,m)=dmclm(k,m) * (molwgt(m)/con_amd)
             endif
            enddo     ! m_do_loop
          endif


!  ---  determine dmanl from dmclm and dmfcs
          do m = 1, NMXG
             dmanl(k,m)= ctaer*dmfcs(i,k,m) +                           &
     &                 ( f_one-ctaer)*dmclm(k,m)
          enddo
        enddo

!> -# Call aeropt_grt() to alculate sw/lw aerosol optical properties
!!    for the corresponding frequency bands.

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

        if ( lsswr ) then

          if ( laswflg ) 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

          else

            aerosw(:,:,:,:) = f_zero

          endif

        endif     ! end if_lsswr_block

        if ( lslwr ) then

          if ( lalwflg ) then

            if ( NLWBND == 1 ) then
              m1 = NBDSW + 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 = NBDSW + 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

          else

            aerolw(:,:,:,:) = f_zero

          endif
        endif     ! end if_lslwr_block

      enddo  lab_do_IMAX

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

!>\ingroup setaer
!> This subroutine maps input tracer fields (trcly) to local tracer
!! array (aermr).
!-----------------------------
      subroutine map_aermr
!.............................
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

! ==================================================================== !
!                                                                      !
! subprogram: map_aermr                                                !
!                                                                      !
!   map input tracer fields (trcly) to local tracer array (aermr)      !
!                                                                      !
!  ====================  defination of variables  ===================  !
!                                                                      !
! input arguments:                                                     !
!     IMAX    - horizontal dimension of arrays                  1      !
!     NLAY    - vertical dimensions of arrays                   1      !
!     trcly   - layer tracer mass mixing ratio     g/g  IMAX*NLAY*NTRAC!
! output arguments: (to module variables)                              !
!     aermr   - layer aerosol mass mixing ratio    g/g  IMAX*NLAY*NMXG !
!                                                                      !
! note:                                                                !
!  NTRAC is the number of tracers excluding water vapor                !
!  NMXG is the number of prognostic aerosol species                    !
!  ==================================================================  !
!
      implicit none

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

!  ---  local:
      integer    :: i, indx, ii
      character  :: tp*2

! initialize
      aermr(:,:,:) = f_zero
      ii = 1        !! <---- trcly does not contain q

! ==>  DU: du1 (submicron bins), du2, du3, du4, du5
       if( gfs_phy_tracer%doing_DU ) then
         aermr(:,:,dm_indx%dust1) = trcly(:,:,dmfcs_indx%du001-ii)
         aermr(:,:,dm_indx%dust2) = trcly(:,:,dmfcs_indx%du002-ii)
         aermr(:,:,dm_indx%dust3) = trcly(:,:,dmfcs_indx%du003-ii)
         aermr(:,:,dm_indx%dust4) = trcly(:,:,dmfcs_indx%du004-ii)
         aermr(:,:,dm_indx%dust5) = trcly(:,:,dmfcs_indx%du005-ii)
       endif

! ==>  OC: oc_phobic, oc_philic
       if( gfs_phy_tracer%doing_OC ) then
         aermr(:,:,dm_indx%waso_phobic) =                               &
     &                     trcly(:,:,dmfcs_indx%ocphobic-ii)
         aermr(:,:,dm_indx%waso_philic) =                               &
     &                     trcly(:,:,dmfcs_indx%ocphilic-ii)
       endif

! ==>  BC: bc_phobic, bc_philic
       if( gfs_phy_tracer%doing_BC ) then
         aermr(:,:,dm_indx%soot_phobic) =                               &
     &                     trcly(:,:,dmfcs_indx%bcphobic-ii)
         aermr(:,:,dm_indx%soot_philic) =                               &
     &                     trcly(:,:,dmfcs_indx%bcphilic-ii)
       endif

! ==>  SS: ss1, ss2 (submicron bins), ss3, ss4, ss5
       if( gfs_phy_tracer%doing_SS ) then
          aermr(:,:,dm_indx%ssam) = trcly(:,:,dmfcs_indx%ss001-ii)      &
     &                            + trcly(:,:,dmfcs_indx%ss002-ii)
          aermr(:,:,dm_indx%sscm) = trcly(:,:,dmfcs_indx%ss003-ii)      &
     &                            + trcly(:,:,dmfcs_indx%ss004-ii)      &
     &                            + trcly(:,:,dmfcs_indx%ss005-ii)
       endif

! ==>  SU: so4
       if( gfs_phy_tracer%doing_SU ) then
          aermr(:,:,dm_indx%suso) = trcly(:,:,dmfcs_indx%so4-ii)
       endif

      return
!...................................
      end subroutine map_aermr
!-----------------------------------


!>\ingroup setaer
!! This subroutine computes aerosols optical properties in NSWLWBD
!! SW/LW bands. Aerosol distribution at each grid point is composed
!! from up to NMXG aerosol species (from NUM_GRIDCOMP components).
!-----------------------------------
      subroutine aeropt_grt
!...................................
!  ---  inputs:  (in scope variables)
!  ---  outputs: (in scope variables)

!  ==================================================================  !
!                                                                      !
!  subprogram: aeropt_grt                                              !
!                                                                      !
!  compute aerosols optical properties in NSWLWBD sw/lw bands.         !
!  Aerosol distribution at each grid point is composed from up to      !
!  NMXG aerosol species (from NUM_GRIDCOMP components).                !
!                                                                      !
!  input variables:                                                    !
!     dmanl  - aerosol dry mass                     g/m3   NLAY*NMXG   !
!     rh1    - relative humidity                     %     NLAY        !
!     dz1    - layer thickness                       km    NLAY        !
!     NLAY   - vertical dimensions                   -     1           !
!     ivflip - control flag for direction of vertical index            !
!               =0: index from toa to surface                          !
!               =1: index from surface to toa                          !
!                                                                      !
!  output variables:                                                   !
!     tauae  - aerosol optical depth                 -   NLAY*NSWLWBD  !
!     ssaae  - aerosol single scattering albedo      -   NLAY*NSWLWBD  !
!     asyae  - aerosol asymmetry parameter           -   NLAY*NSWLWBD  !
!                                                                      !
!  ==================================================================  !
!
      implicit none

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

!  ---  locals:
      real (kind=kind_phys) :: aerdm
      real (kind=kind_phys) :: ext1, ssa1, asy1, ex00, ss00, as00,      &
     &                         ex01, ss01, as01, exint
      real (kind=kind_phys) :: tau, ssa, asy,                           &
     &                         sum_tau, sum_ssa, sum_asy

!  ---  subgroups for sub-micron dust
!  ---  corresponds to 0.1-0.18, 0.18-0.3, 0.3-0.6, 0.6-1.0 micron

      real (kind=kind_phys) ::   fd(4)
      data fd  / 0.01053,0.08421,0.25263,0.65263 /

      character    :: tp*2
      integer      :: icmp, n, kk, ib, ih2, ih1, ii, ij, ijk
      real (kind=kind_phys) :: drh0, drh1, rdrh

      real (kind=kind_phys) :: qmin  !<--lower bound for opt calc
      data qmin  / 1.e-20 /

!===>  ...  begin here

! --- initialize (assume no aerosols)
      tauae = f_zero
      ssaae = f_one
      asyae = f_zero

      tauae_gocart = f_zero

!===> ... loop over vertical layers
!
      lab_do_layer : do kk = 1, NLAY

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

        ih2 = 1
        do while ( rh1(kk) > 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(kk) - rhlev_grt(ih1)
        if ( ih1 == ih2 ) then
          rdrh = f_zero
        else
          rdrh = drh1 / drh0
        endif

! --- loop through sw/lw spectral bands

        lab_do_ib : do ib = 1, NSWLWBD
          sum_tau = f_zero
          sum_ssa = f_zero
          sum_asy = f_zero

! --- loop through aerosol grid components
          lab_do_icmp : do icmp = 1, NUM_GRIDCOMP
            ext1 = f_zero
            ssa1 = f_zero
            asy1 = f_zero

            tp = gridcomp(icmp)

            select case ( tp )

! -- dust aerosols: no humidification effect
            case ( 'DU')
              do n = 1, KCM1

                if (n <= 4) then
                  aerdm = dmanl(kk,dm_indx%dust1) * fd(n)
                else
                  aerdm = dmanl(kk,dm_indx%dust1+n-4 )
                endif

                if (aerdm < qmin) aerdm = f_zero
                ex00 = extrhi_grt(n,ib)*(1000.*dz1(kk))*aerdm
                ss00 = ssarhi_grt(n,ib)
                as00 = asyrhi_grt(n,ib)
                ext1 = ext1 + ex00
                ssa1 = ssa1 + ex00 * ss00
                asy1 = asy1 + ex00 * ss00 * as00

              enddo

! -- suso aerosols: with humidification effect
            case ( 'SU')
              ij = isuso
              exint = extrhd_grt(ih1,ij,ib)                             &
     &          + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
              ss00 = ssarhd_grt(ih1,ij,ib)                              &
     &          + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
              as00 = asyrhd_grt(ih1,ij,ib)                              &
     &          + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))

              aerdm = dmanl(kk, dm_indx%suso)
              if (aerdm < qmin) aerdm = f_zero
              ex00 = exint*(1000.*dz1(kk))*aerdm
              ext1 = ex00
              ssa1 = ex00 * ss00
              asy1 = ex00 * ss00 * as00

! -- seasalt aerosols: with humidification effect
            case ( 'SS')
              do n = 1, 2                  !<---- ssam, sscm
                ij = issam + (n-1)
                exint = extrhd_grt(ih1,ij,ib)                           &
     &          + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
                ss00 = ssarhd_grt(ih1,ij,ib)                            &
     &          + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
                as00 = asyrhd_grt(ih1,ij,ib)                            &
     &          + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))

                aerdm = dmanl(kk, dm_indx%ssam+n-1)
                if (aerdm < qmin) aerdm = f_zero
                ex00 = exint*(1000.*dz1(kk))*aerdm
                ext1 = ext1 + ex00
                ssa1 = ssa1 + ex00 * ss00
                asy1 = asy1 + ex00 * ss00 * as00

              enddo

! -- organic carbon/black carbon:
!    using 'waso' and 'soot' for hydrophilic OC and BC
!    using 'waso' and 'soot' at RH=0 for hydrophobic OC and BC
            case ( 'OC', 'BC')
              if(tp == 'OC') then
                 ii = dm_indx%waso_phobic
                 ij = iwaso
              else
                 ii = dm_indx%soot_phobic
                 ij = isoot
              endif

! --- hydrophobic
              aerdm = dmanl(kk, ii)
              if (aerdm < qmin) aerdm = f_zero
              ex00 = extrhd_grt(1,ij,ib)*(1000.*dz1(kk))*aerdm
              ss00 = ssarhd_grt(1,ij,ib)
              as00 = asyrhd_grt(1,ij,ib)
! --- hydrophilic
              aerdm = dmanl(kk, ii+1)
              if (aerdm < qmin) aerdm = f_zero
              exint = extrhd_grt(ih1,ij,ib)                             &
     &         + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
              ex01 = exint*(1000.*dz1(kk))*aerdm
              ss01 = ssarhd_grt(ih1,ij,ib)                              &
     &          + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
              as01 = asyrhd_grt(ih1,ij,ib)                              &
     &          + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))

              ext1 = ex00 + ex01
              ssa1 = (ex00 * ss00) + (ex01 * ss01)
              asy1 = (ex00 * ss00 * as00) + (ex01 * ss01 * as01)

            end select

! --- determine tau, ssa, asy for each grid component
            tau = ext1
            if (ext1 > f_zero) ssa=min(f_one,ssa1/ext1)
            if (ssa1 > f_zero) asy=min(f_one,asy1/ssa1)

! --- save tau at 550 nm for each grid component
            if ( ib == nv_aod ) then
              do ijk = 1, max_num_gridcomp
                if ( tp == max_gridcomp(ijk) )  then
                   tauae_gocart(kk,ijk) = tau
                endif
              enddo
            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 lab_do_icmp


! --- determine total tau, ssa, asy for aerosol mixture
          tauae(kk,ib) = sum_tau
          if (sum_tau > f_zero) ssaae(kk,ib) = sum_ssa / sum_tau
          if (sum_ssa > f_zero) asyae(kk,ib) = sum_asy / sum_ssa

        enddo   lab_do_ib

      enddo lab_do_layer

!
      return
!...................................
      end subroutine aeropt_grt
!--------------------------------

!................................
      end subroutine setgocartaer
!--------------------------------
!! @}
!
! GOCART code modification end here (Sarah Lu)  ------------------------!
! =======================================================================

!..........................................!
      end module module_radiation_aerosols !
!==========================================!
!> @}