radiation_aerosols.f

Go to the documentation of this file.

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


!========================================!
      module module_radiation_aerosols   !
!........................................!
!
      use physparam,only : iaermdl, iaerflg, lavoflg, lalwflg, laswflg, &
     &                     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:
      integer, parameter, public :: nf_aesw = 3     ! num of output fields for sw rad
      integer, parameter, public :: nf_aelw = 3     ! num of output fields for lw rad
      integer, parameter, public :: nlwstr  = 1     ! starting band number in ir region
      integer, parameter, public :: nspc    = 5     ! num of species for output aod (opnl)
      integer, parameter, public :: nspc1   = nspc + 1   ! total + species

      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"
      integer, save :: nswbnd  = nbdsw       ! number of actual bands for sw aerosols
                                             ! calculated according to laswflg setting
      integer, save :: nlwbnd  = nbdlw       ! number of actual bands for lw aerosols
                                             ! calculated according to lalwflg and lalw1bd settings
      integer, save :: nswlwbd = nbdsw+nbdlw ! total number of bands for sw+lw aerosols

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

!  ---  parameter constants:
      integer, parameter, public :: nwvsol  = 151   ! num of wvnum regions where solar
                                                    ! flux is constant
      integer, parameter, public :: nwvtot  = 57600 ! total num of wvnum included
      integer, parameter, public :: nwvtir  = 4000  ! total num of wvnum in ir range

!  ---  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 (w/m**2) 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  /

! --------------------------------------------------------------------- !
!   section-2 : module variables for stratospheric volcanic aerosols    !
!               from historical data (sato et al. 1993)                 !
! --------------------------------------------------------------------- !

!  ---  parameter constants:
      integer, parameter :: minvyr = 1850    ! lower lim (year) data available
      integer, parameter :: maxvyr = 1999    ! upper lim (year) data available

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

!  ---  static control variables:
      integer :: kyrstr, kyrend, kyrsav, kmonsav

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

!  ---  parameters and constants:
      integer, parameter :: nxc = 5    ! num of max componets in a profile
      integer, parameter :: nae = 7    ! num of aerosols profile structures
      integer, parameter :: ndm = 5    ! num of atmos aerosols domains
      integer, parameter :: imxae = 72 ! num of lon-points in glb aeros data set
      integer, parameter :: jmxae = 37 ! num of lat-points in glb aeros data set
      integer, parameter :: naerbnd=61 ! num of bands for clim aer data (opac)
      integer, parameter :: nrhlev =8  ! num of rh levels for rh-dep components
      integer, parameter :: ncm1 = 6   ! num of rh independent aeros species
      integer, parameter :: ncm2 = 4   ! num of rh dependent aeros species
      integer, parameter :: ncm  = ncm1+ncm2

      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)

      real (kind=kind_phys), save, dimension(NDM,NAE) :: haer, prsref,  &
     &                                                   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

      real (kind=kind_phys), dimension(NXC,IMXAE,JMXAE), save :: cmixg
      real (kind=kind_phys), dimension( 2 ,IMXAE,JMXAE), save :: denng
      integer,               dimension(NXC,IMXAE,JMXAE), save :: idxcg
      integer,               dimension(    IMXAE,JMXAE), save :: kprfg

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

!  ---  parameters and constants:
!   - KCM, KCM1, KCM2 are determined from subroutine 'set_aerspc'
      integer, parameter :: kaerbnd=61 ! num of bands for aer data (gocart)
      integer, parameter :: krhlev =36 ! num of rh levels for rh-dep components
!*    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
      integer, save      :: kcm1 = 0   ! num of rh indep aerosols (set in subr set_aerspc)
      integer, save      :: kcm2 = 0   ! num of rh dep aerosols   (set in subr set_aerspc)
      integer, save      :: kcm        ! =KCM1+KCM2               (set in subr set_aerspc)

      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

      integer,               allocatable, dimension(:) :: iendwv_grt
      real (kind=kind_phys), allocatable, dimension(:,:)  ::            &
     &                       rhidext0_grt, rhidssa0_grt, rhidasy0_grt
      real (kind=kind_phys), allocatable, dimension(:,:,:)::            &
     &                       rhdpext0_grt, rhdpssa0_grt, 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

      real (kind=kind_phys), allocatable, save, dimension(:,:)   ::     &
     &       extrhi_grt, ssarhi_grt, asyrhi_grt
      real (kind=kind_phys), allocatable, save, dimension(:,:,:) ::     &
     &       extrhd_grt, ssarhd_grt, asyrhd_grt

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

      integer, parameter :: imxg = 144 ! num of lon-points in geos dataset
      integer, parameter :: jmxg = 91  ! num of lat-points in geos dataset
      integer, parameter :: kmxg = 30  ! num of vertical layers in geos dataset
!*    integer, parameter :: NMXG = 12  ! num of gocart aer spec for opt calc
      integer, save      :: nmxg       ! to be determined by set_aerspc

      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

      real (kind=kind_phys),allocatable, save:: psclmg(:,:,:),          &
     &                                          dmclmg(:,:,:,:)

!   - geos-gocart lat/lon arrays
      real (kind=kind_phys), allocatable, save, dimension(:)   ::       &
     &                       geos_rlon, 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 or clim gocart aerosol fields
      logical, save :: get_fcst = .true.    ! option to get fcst field
      logical, save :: get_clim = .true.    ! option to get clim field

!  ------  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)
      type  gocart_index_type
         integer :: dust1, dust2, dust3, dust4, dust5,                  &   ! dust
     &              ssam,  sscm,                                        &   ! sea salt
     &              suso,                                               &   ! sulfate
     &              waso_phobic, waso_philic,                           &   ! oc
     &              soot_phobic, soot_philic                                ! bc
      endtype
      type(gocart_index_type), save :: dm_indx

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

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

!  --- default full-package setting
      integer, parameter          :: max_num_gridcomp = 5
      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
!
      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)
!
      real (kind=kind_phys), parameter :: wvn550 = 1.0e4/0.55
      integer, save      :: nv_aod = 1

!  ---  public interfaces

      public aer_init, aer_update, setaer


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

!-----------------------------------
      subroutine aer_init                                               &
!...................................
!  ---  inputs:
     &     ( nlay, me )
!  ---  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     !
!     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

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

      if ( iaerflg /= 100 ) then

!  --- ...  set up spectral one wavenumber solar/ir fluxes

        call set_spectrum
!  ---  inputs:   (module constants)
!  ---  outputs:  (in-scope variables)

!  --- ...  invoke tropospheric aerosol initialization

        if ( iaermdl == 0 ) 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

!  --- ...  invoke stratosperic volcanic aerosol initialization

      if ( lavoflg ) then

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

      endif    ! end if_lavoflg_block


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

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

!--------------------------------
      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 (w/m2)!
!   eirfwv(NWVTIR)   - ir flux(273k) for each individual wavenum (w/m2)!
!                                                                      !
!  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 = 2.0 * con_pi * con_plnk * (con_c**2)
      tmp2 = con_plnk * con_c / (con_boltz * con_t0c)

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

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

!-----------------------------------
      subroutine clim_aerinit                                           &
!...................................
!  ---  inputs:
     &     ( solfwv, eirfwv, me                                         &
!  ---  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
!  ---  inputs:   (in-scope variables, module constants)
!  ---  outputs:  (module variables)


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

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

      real (kind=kind_phys) :: sumsol, sumir

      logical :: file_exist
      character :: cline*80
!
!===>  ...  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',               &
     &        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

      read(niaercm,21) cline   ! ending wave num for 61 aeros spectral bands
  21  format(a80)
      read(niaercm,22) iendwv(:)
  22  format(13i6)

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

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

      read(niaercm,21) cline   ! rh indep ext coef for 61 bands, 6 species
      read(niaercm,28) rhidext0(:,:)
  28  format(8e10.3)

      read(niaercm,21) cline   ! rh indep sca coef for 61 bands, 6 species
      read(niaercm,28) rhidsca0(:,:)

      read(niaercm,21) cline   ! rh indep ssa coef for 61 bands, 6 species
      read(niaercm,28) rhidssa0(:,:)

      read(niaercm,21) cline   ! rh indep asy coef for 61 bands, 6 species
      read(niaercm,28) rhidasy0(:,:)

      read(niaercm,21) cline   ! rh dep ext coef for 61 bands, 8 rh lev, 4 species
      read(niaercm,28) rhdpext0(:,:,:)

      read(niaercm,21) cline   ! rh dep sca coef for 61 bands, 8 rh lev, 4 species
      read(niaercm,28) rhdpsca0(:,:,:)

      read(niaercm,21) cline   ! rh dep ssa coef for 61 bands, 8 rh lev, 4 species
      read(niaercm,28) rhdpssa0(:,:,:)

      read(niaercm,21) cline   ! rh dep asy coef for 61 bands, 8 rh lev, 4 species
      read(niaercm,28) rhdpasy0(:,:,:)

      read(niaercm,21) cline   ! stratospheric background aeros for 61 bands
      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
        solwaer(:,:) = f_zero

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

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

          lab_swdowhile : do while ( iw1 > iendwv(ii) )
            if ( ii == naerbnd ) exit lab_swdowhile
            ii = ii + 1
          enddo  lab_swdowhile

          sumsol = f_zero
          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

          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
        eirwaer(:,:) = f_zero

        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(wvnlw1(mb))
            iw2 = nint(wvnlw2(mb))
          endif

          lab_lwdowhile : do while ( iw1 > iendwv(ii) )
            if ( ii == naerbnd ) exit lab_lwdowhile
            ii = ii + 1
          enddo  lab_lwdowhile

          sumir = f_zero
          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

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

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

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

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

!  ---  for rh independent aerosol species

          do nc = 1, ncm1
            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

!  ---  for rh dependent aerosols species

          do nc = 1, ncm2
            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

        do nb = 1, nlwbnd

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

!  ---  for rh independent aerosol species

          do nc = 1, ncm1
            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

!  ---  for rh dependent aerosols species

          do nc = 1, ncm2
            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
!-----------------------------------


!-----------------------------------
      subroutine aer_update                                             &
!...................................
!  ---  inputs:
     &     ( iyear, imon, me )
!  ---  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

      if ( lalwflg .or. laswflg ) then ! update monthly tropspheric aerosol data
        call trop_update
      endif

      if ( lavoflg ) then              ! update yearly stratospheric volcanic aerosol data
        call volc_update
      endif


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

!--------------------------------
      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',                &
     &        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

      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

      do j = 1, jmxae
        do i = 1, imxae
          do m = 1, 2
            denng(m,i,j) = f_zero
          enddo
        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
!--------------------------------

!--------------------------------
      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',          &
     &            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
!-----------------------------------


!-----------------------------------
      subroutine setaer                                                 &
!...................................

!  ---  inputs:
     &     ( prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,xlon,xlat,         &
     &       imax,nlay,nlp1, lsswr,lslwr,                               &
!  ---  outputs:
     &       aerosw,aerolw                                              &
     &,      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

!SARAH
!         if ( iaerflg == 1 ) then      ! use opac aerosol climatology
          if ( iaermdl == 0 ) 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

!  ---  ...  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     ( wvnsw1(mb) > 20000 ) then   ! range of wvlth < 0.5mu
                tmp2 = 0.74
              elseif ( wvnsw2(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 * (wvnsw2(mb)+wvnsw1(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 * (wvnlw2(m) + wvnlw1(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     ( wvnsw1(mb) > 20000 ) then   ! range of wvlth < 0.5mu
                tmp2 = 0.74
              elseif ( wvnsw2(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 * (wvnsw2(mb)+wvnsw1(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 * (wvnlw2(m) + wvnlw1(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
!-----------------------------------

!-----------------------------------
      subroutine aer_property                                           &
!...................................

!  ---  inputs:
     &     ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,                   &
     &       alon,alat,slmsk, laersw,laerlw,                            &
     &       imax,nlay,nlp1,                                            &
!    &       IMAX,NLAY,NLP1,NSPC,                                       &
!  ---  outputs:
     &       aerosw,aerolw,aerodp                                       &
     &     )

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

      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

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

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


!-----------------------------------
      subroutine gocart_init                                            &
!...................................
!  ---  inputs:
     &     ( nwvtot,solfwv,soltot,nwvtir,eirfwv,                        &
     &       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

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

      integer :: i, mb, ib, ii, iw, iw1, iw2

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

        do ib = 1, nbdsw
          mb = ib + nswstr - 1
          ii = 1
          iw1 = nint(wvnsw1(mb))
          iw2 = nint(wvnsw2(mb))
!
! ---  locate the spectral band for 550nm (for aod diag)
!
          if (10000./iw1 >= 0.55 .and.                                  &
     &        10000./iw2 <= 0.55 )  then
              nv_aod =  ib
          endif

          lab_swdowhile : do while ( iw1 > iendwv_grt(ii) )
            if ( ii == kaerbnd ) exit lab_swdowhile
            ii = ii + 1
          enddo  lab_swdowhile

          sumsol = f_zero
          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

          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(wvnsw1(mb))
          iw2 = nint(wvnsw2(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

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

          lab_lwdowhile : do while ( iw1 > iendwv_grt(ii) )
            if ( ii == kaerbnd ) exit lab_lwdowhile
            ii = ii + 1
          enddo  lab_lwdowhile

          sumir = f_zero
          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(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,wvnsw1(mb),wvnsw2(mb)
           print *,'RAD -lamda1,lamda2: ',ib,10000./wvnsw1(mb),         &
     &                                   10000./wvnsw2(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,wvnlw1(ib),wvnlw2(ib)
           print *,'RAD -lamda1,lamda2: ',ib,10000./wvnlw1(ib),         &
     &                                   10000./wvnlw2(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
! =================
!-----------------------------
      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  
           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                   
           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

!-----------------------------------
!-----------------------------
      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',              &
     &        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    
          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    
          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
!-----------------------------------
! 
!-----------------------------
      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
!--------------------------------
!
!-----------------------------------
      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),                          &
     &       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),                         &
     &        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),                       &
     &         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),                      &
     &           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
!-----------------------------------

!-----------------------------------
      subroutine setgocartaer                                           &
!...................................

!  ---  inputs:
     &     ( alon,alat,prslk,rhlay,dz,hz,nswlwbd,                       &
     &       prsl,tvly,trcly,                                           &
     &       imax,nlay,nlp1, ivflip, lsswr,lslwr,                       &
!  ---  outputs:
     &       aerosw,aerolw                                              &
     &     )


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

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

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

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