MODULE module_ftuv_subs

      use module_wave_data, only : nw

      implicit none

      integer, parameter :: dp = selected_real_kind(14,300)
!-----------------------------------------------------------------------------
!     data for inter2, inter3 and interi
!-----------------------------------------------------------------------------
      integer, private, save :: nintervals
      integer, private, allocatable, save :: xi(:), xcnt(:)
      real(dp), private, allocatable, save :: xfrac(:,:)

      real(dp), private, parameter :: km2cm = 1.e5_dp

!-----------------------------------------------------------------------------
!     data for lyman
!-----------------------------------------------------------------------------
      integer, private, parameter :: nla = 2
      integer, private, parameter :: ngast = 17

!-----------------------------------------------------------------------------
!     data for schu_inti, schu
!-----------------------------------------------------------------------------
      integer, private, parameter  :: tdim = 501
      !BSINGH(PNNL)-Lahey compiler doesn't allow use of 'real' in variable declaration
      real(dp), private, parameter  :: tdim_real = 501.0
      !real(dp), private, parameter :: t_del = 5._dp/real(tdim-1,kind=dp)!BSINGH(PNNL)-commented out
      !real(dp), private, parameter :: t_fac = real(tdim-1,kind=dp)/5._dp!BSINGH(PNNL)-commented out
      real(dp), private, parameter :: t_del = 5._dp/(tdim_real-1._dp)!BSINGH(PNNL)-redefined t_del and t_fac
      real(dp), private, parameter :: t_fac = (tdim_real-1._dp)/5._dp

      integer, private :: ii, jj
      real(dp), private, save :: d_table(0:tdim,ngast-1)
      real(dp), private, save :: x_table(0:tdim,ngast-1)
      real(dp), private, save :: o2_table(tdim)
      real(dp), private, dimension(12,ngast-1), save :: a_schu, b_schu
!-----------------------------------------------------------------------------
! a_schu(16,12) coefficients for rj(m) (table 1 in kockarts 1994)
! b_schu(16,12) rj(o2)(table 2 in kockarts 1994)
! rjm attenuation coefficients rj(m)
! rjo2 rj(o2)
!-----------------------------------------------------------------------------
      data ((a_schu(jj,ii),jj=1,12),ii=1,ngast-1) /                             &
!a 57000-56500.5 cm-1
       1.13402e-01,1.00088e-20,3.48747e-01,2.76282e-20,3.47322e-01,1.01267e-19, &
       1.67351e-01,5.63588e-19,2.31433e-02,1.68267e-18,0.00000e+00,0.00000e+00, &
!a 56500-56000.5 cm-1
       2.55268e-03,1.64489e-21,1.85483e-01,2.03591e-21,2.60603e-01,4.62276e-21, &
       2.50337e-01,1.45106e-20,1.92340e-01,7.57381e-20,1.06363e-01,7.89634e-19, &
!a 56000-55500.5 cm-1
       4.21594e-03,8.46639e-22,8.91886e-02,1.12935e-21,2.21334e-01,1.67868e-21, &
       2.84446e-01,3.94782e-21,2.33442e-01,1.91554e-20,1.63433e-01,2.25346e-19, &
!a 55500-55000.5 cm-1
       3.93529e-03,6.79660e-22,4.46906e-02,9.00358e-22,1.33060e-01,1.55952e-21, &
       3.25506e-01,3.43763e-21,2.79405e-01,1.62086e-20,2.10316e-01,1.53883e-19, &
!a 55000-54500.5 cm-1
       2.60939e-03,2.33791e-22,2.08101e-02,3.21734e-22,1.67186e-01,5.77191e-22, &
       2.80694e-01,1.33362e-21,3.26867e-01,6.10533e-21,1.96539e-01,7.83142e-20, &
!a 54500-54000.5 cm-1
       9.33711e-03,1.32897e-22,3.63980e-02,1.78786e-22,1.46182e-01,3.38285e-22, &
       3.81762e-01,8.93773e-22,2.58549e-01,4.28115e-21,1.64773e-01,4.67537e-20, &
!a 54000-53500.5 cm-1
       9.51799e-03,1.00252e-22,3.26320e-02,1.33766e-22,1.45962e-01,2.64831e-22, &
       4.49823e-01,6.42879e-22,2.14207e-01,3.19594e-21,1.45616e-01,2.77182e-20, &
!a 53500-53000.5 cm-1
       7.87331e-03,3.38291e-23,6.91451e-02,4.77708e-23,1.29786e-01,8.30805e-23, &
       3.05103e-01,2.36167e-22,3.35007e-01,8.59109e-22,1.49766e-01,9.63516e-21, &
!a 53000-52500.5 cm-1
       6.92175e-02,1.56323e-23,1.44403e-01,3.03795e-23,2.94489e-01,1.13219e-22, &
       3.34773e-01,3.48121e-22,9.73632e-02,2.10693e-21,5.94308e-02,1.26195e-20, &
!a 52500-52000.5 cm-1
       1.47873e-01,8.62033e-24,3.15881e-01,3.51859e-23,4.08077e-01,1.90524e-22, &
       8.08029e-02,9.93062e-22,3.90399e-02,6.38738e-21,8.13330e-03,9.93644e-22, &
!a 52000-51500.5 cm-1
       1.50269e-01,1.02621e-23,2.39823e-01,3.48120e-23,3.56408e-01,1.69494e-22, &
       1.61277e-01,6.59294e-22,8.89713e-02,2.94571e-21,3.25063e-03,1.25548e-20, &
!a 51500-51000.5 cm-1
       2.55746e-01,8.49877e-24,2.94733e-01,2.06878e-23,2.86382e-01,9.30992e-23, &
       1.21011e-01,3.66239e-22,4.21105e-02,1.75700e-21,0.00000e+00,0.00000e+00, &
!a 51000-50500.5 cm-1
       5.40111e-01,7.36085e-24,2.93263e-01,2.46742e-23,1.63417e-01,1.37832e-22, &
       3.23781e-03,2.15052e-21,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!a 50500-50000.5 cm-1
       8.18514e-01,7.17937e-24,1.82262e-01,4.17496e-23,0.00000e+00,0.00000e+00, &
       0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!a 50000-49500.5 cm-1
       8.73680e-01,7.13444e-24,1.25583e-01,2.77819e-23,0.00000e+00,0.00000e+00, &
       0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!a 49500-49000.5 cm-1
       3.32476e-04,7.00362e-24,9.89000e-01,6.99600e-24,0.00000e+00,0.00000e+00, &
       0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 /
      data ((b_schu(jj,ii),jj=1,12),ii=1,ngast-1) /                             &
! 57000-56500.5 cm-1
       1.07382e-21,9.95029e-21,7.19430e-21,2.48960e-20,2.53735e-20,7.54467e-20, &
       4.48987e-20,2.79981e-19,9.72535e-20,9.29745e-19,2.30892e-20,4.08009e-17, &
! 56500-56000.5 cm-1
       3.16903e-22,1.98251e-21,5.87326e-22,3.44057e-21,2.53094e-21,8.81484e-21, &
       8.82299e-21,4.17179e-20,2.64703e-20,2.43792e-19,8.73831e-20,1.46371e-18, &
! 56000-55500.5 cm-1
       1.64421e-23,9.26011e-22,2.73137e-22,1.33640e-21,9.79188e-22,2.99706e-21, &
       3.37768e-21,1.39438e-20,1.47898e-20,1.04322e-19,4.08014e-20,6.31023e-19, &
! 55500-55000.5 cm-1
       8.68729e-24,7.31056e-22,8.78313e-23,1.07173e-21,8.28170e-22,2.54986e-21, &
       2.57643e-21,9.42698e-21,9.92377e-21,5.21402e-20,3.34301e-20,2.91785e-19, &
! 55000-54500.5 cm-1
       1.20679e-24,2.44092e-22,2.64326e-23,4.03998e-22,2.53514e-22,8.53166e-22, &
       1.29834e-21,3.74482e-21,5.12103e-21,2.65798e-20,2.10948e-20,2.35315e-19, &
! 54500-54000.5 cm-1
       2.79656e-24,1.40820e-22,3.60824e-23,2.69510e-22,4.02850e-22,8.83735e-22, &
       1.77198e-21,6.60221e-21,9.60992e-21,8.13558e-20,4.95591e-21,1.22858e-17, &
! 54000-53500.5 cm-1
       2.36959e-24,1.07535e-22,2.83333e-23,2.16789e-22,3.35242e-22,6.42753e-22, &
       1.26395e-21,5.43183e-21,4.88083e-21,5.42670e-20,3.27481e-21,1.58264e-17, &
! 53500-53000.5 cm-1
       8.65018e-25,3.70310e-23,1.04351e-23,6.43574e-23,1.17431e-22,2.70904e-22, &
       4.88705e-22,1.65505e-21,2.19776e-21,2.71172e-20,2.65257e-21,2.13945e-17, &
! 53000-52500.5 cm-1
       9.63263e-25,1.54249e-23,4.78065e-24,2.97642e-23,6.40637e-23,1.46464e-22, &
       1.82634e-22,7.12786e-22,1.64805e-21,2.37376e-17,9.33059e-22,1.13741e-20, &
! 52500-52000.5 cm-1
       1.08414e-24,8.37560e-24,9.15550e-24,2.99295e-23,9.38405e-23,1.95845e-22, &
       2.84356e-22,3.39699e-21,1.94524e-22,2.72227e-19,1.18924e-21,3.20246e-17, &
! 52000-51500.5 cm-1
       1.52817e-24,1.01885e-23,1.22946e-23,4.16517e-23,9.01287e-23,2.34869e-22, &
       1.93510e-22,1.44956e-21,1.81051e-22,5.17773e-21,9.82059e-22,6.22768e-17, &
! 51500-51000.5 cm-1
       2.12813e-24,8.48035e-24,5.23338e-24,1.93052e-23,1.99464e-23,7.48997e-23, &
       4.96642e-22,6.15691e-17,4.47504e-23,2.76004e-22,8.26788e-23,1.65278e-21, &
! 51000-50500.5 cm-1
       3.81336e-24,7.32307e-24,5.60549e-24,2.04651e-23,3.36883e-22,6.15708e-17, &
       2.09877e-23,1.07474e-22,9.13562e-24,8.41252e-22,0.00000e+00,0.00000e+00, &
! 50500-50000.5 cm-1
       5.75373e-24,7.15986e-24,5.90031e-24,3.05375e-23,2.97196e-22,8.92000e-17, &
       8.55920e-24,1.66709e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
! 50000-49500.5 cm-1
       6.21281e-24,7.13108e-24,3.30780e-24,2.61196e-23,1.30783e-22,9.42550e-17, &
       2.69241e-24,1.46500e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
! 49500-49000.5 cm-1
       6.81118e-24,6.98767e-24,7.55667e-25,2.75124e-23,1.94044e-22,1.45019e-16, &
       1.92236e-24,3.73223e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 /

!-----------------------------------------------------------------------------
!     data for set_o2_xsect      
!-----------------------------------------------------------------------------
      integer, public,  parameter :: nwint = 105

      real(dp), public, save  :: wlint(nwint)
      real(dp), private, save :: xso2int(nwint)
      real(dp), private, save :: wlla(nla)
      real(dp), private, save :: wlgast(ngast)

!-----------------------------------------------------------------------------
! o2 parameters
!-----------------------------------------------------------------------------
      data (wlint(ii),ii=1,105)/  &
            0.0000E+00, 0.1166E+03, 0.1167E+03, 0.1173E+03, 0.1179E+03,  &
            0.1187E+03, 0.1194E+03, 0.1202E+03, 0.1208E+03, 0.1214E+03,  &
            0.1219E+03, 0.1223E+03, 0.1231E+03, 0.1238E+03, 0.1246E+03,  &
            0.1254E+03, 0.1262E+03, 0.1270E+03, 0.1278E+03, 0.1286E+03,  &
            0.1294E+03, 0.1303E+03, 0.1311E+03, 0.1320E+03, 0.1329E+03,  &
            0.1338E+03, 0.1346E+03, 0.1356E+03, 0.1365E+03, 0.1374E+03,  &
            0.1384E+03, 0.1399E+03, 0.1418E+03, 0.1439E+03, 0.1459E+03,  &
            0.1481E+03, 0.1504E+03, 0.1526E+03, 0.1550E+03, 0.1574E+03,  &
            0.1600E+03, 0.1626E+03, 0.1653E+03, 0.1681E+03, 0.1709E+03,  &
            0.1731E+03, 0.1746E+03, 0.1754E+03, 0.1770E+03, 0.1786E+03,  &
            0.1802E+03, 0.1818E+03, 0.1835E+03, 0.1852E+03, 0.1869E+03,  &
            0.1887E+03, 0.1905E+03, 0.1923E+03, 0.1942E+03, 0.1961E+03,  &
            0.1980E+03, 0.2000E+03, 0.2020E+03, 0.2041E+03, 0.2050E+03,  &
            0.2060E+03, 0.2070E+03, 0.2080E+03, 0.2090E+03, 0.2100E+03,  &
            0.2110E+03, 0.2120E+03, 0.2130E+03, 0.2140E+03, 0.2150E+03,  &
            0.2160E+03, 0.2170E+03, 0.2180E+03, 0.2190E+03, 0.2200E+03,  &
            0.2210E+03, 0.2220E+03, 0.2230E+03, 0.2240E+03, 0.2250E+03,  &
            0.2260E+03, 0.2270E+03, 0.2280E+03, 0.2290E+03, 0.2300E+03,  &
            0.2310E+03, 0.2320E+03, 0.2330E+03, 0.2340E+03, 0.2350E+03,  &
            0.2360E+03, 0.2370E+03, 0.2380E+03, 0.2390E+03, 0.2400E+03,  &
            0.2410E+03, 0.2420E+03, 0.2430E+03, 0.2430E+03, 0.1000E+39 /

!-----------------------------------------------------------------------------
! o2 parameters
!-----------------------------------------------------------------------------
      data (xso2int(ii),ii=1,105)/  &
            0.0000E+00, 0.1000E-19, 0.6350E-18, 0.7525E-18, 0.1425E-18,  &
            0.2025E-18, 0.2413E-17, 0.6400E-17, 0.5251E-17, 0.7329E-18,  &
            0.3445E-18, 0.3425E-18, 0.3925E-18, 0.8918E-17, 0.9197E-17,  &
            0.6625E-18, 0.2700E-18, 0.1575E-18, 0.3240E-18, 0.4990E-18,  &
            0.4875E-18, 0.5525E-18, 0.1068E-17, 0.1850E-17, 0.2275E-17,  &
            0.3425E-17, 0.5890E-17, 0.8365E-17, 0.1090E-16, 0.1275E-16,  &
            0.1340E-16, 0.1380E-16, 0.1440E-16, 0.1445E-16, 0.1350E-16,  &
            0.1220E-16, 0.1071E-16, 0.9075E-17, 0.7410E-17, 0.5775E-17,  &
            0.4210E-17, 0.2765E-17, 0.1655E-17, 0.9760E-18, 0.5900E-18,  &
            0.3660E-18, 0.2061E-18, 0.3889E-19, 0.4092E-20, 0.2273E-20,  &
            0.1256E-20, 0.6902E-21, 0.3750E-21, 0.2097E-21, 0.1226E-21,  &
            0.6508E-22, 0.4579E-22, 0.3220E-22, 0.2507E-22, 0.1941E-22,  &
            0.1609E-22, 0.1319E-22, 0.9907E-23, 0.7419E-23, 0.7240E-23,  &
            0.7090E-23, 0.6955E-23, 0.6770E-23, 0.6595E-23, 0.6375E-23,  &
            0.6145E-23, 0.5970E-23, 0.5805E-23, 0.5655E-23, 0.5470E-23,  &
            0.5240E-23, 0.5005E-23, 0.4760E-23, 0.4550E-23, 0.4360E-23,  &
            0.4175E-23, 0.3990E-23, 0.3780E-23, 0.3560E-23, 0.3330E-23,  &
            0.3095E-23, 0.2875E-23, 0.2875E-23, 0.2875E-23, 0.2700E-23,  &
            0.2530E-23, 0.2340E-23, 0.2175E-23, 0.2020E-23, 0.1860E-23,  &
            0.1705E-23, 0.1555E-23, 0.1410E-23, 0.1280E-23, 0.1160E-23,  &
            0.1055E-23, 0.9550E-24, 0.4500E-24, 0.0000E+00, 0.0000E+00 /

!-----------------------------------------------------------------------------
! the lyman-alpha grid
!-----------------------------------------------------------------------------
      data (wlla(ii),ii=1,2)/ 121.4_dp, 121.9_dp /

!-----------------------------------------------------------------------------
! o2 parameters
!-----------------------------------------------------------------------------
      data (wlgast(ii),ii=1,17)/  &
            0.1754E+03_dp, 0.1770E+03_dp, 0.1786E+03_dp, 0.1802E+03_dp, 0.1818E+03_dp,  &
            0.1835E+03_dp, 0.1852E+03_dp, 0.1869E+03_dp, 0.1887E+03_dp, 0.1905E+03_dp,  &
            0.1923E+03_dp, 0.1942E+03_dp, 0.1961E+03_dp, 0.1980E+03_dp, 0.2000E+03_dp,  &
            0.2020E+03_dp, 0.2041E+03_dp /

!-----------------------------------------------------------------------------
! setaer arrays
!-----------------------------------------------------------------------------
      real(dp) :: op_cb1w(nw-1), om_cb1w(nw-1), g_cb1w(nw-1)
      real(dp) :: op_cb2w(nw-1), om_cb2w(nw-1), g_cb2w(nw-1)
      real(dp) :: op_oc1w(nw-1), g_oc1w(nw-1)
      real(dp) :: op_oc2w(nw-1), om_oc2w(nw-1), g_oc2w(nw-1)
      real(dp) :: op_antw(nw-1), g_antw(nw-1)
      real(dp) :: op_so4w(nw-1)
      real(dp) :: op_salw(nw-1)

      CONTAINS

      subroutine photoin( chem_opt, nz, zen, o3toms, esfact,     &
                          o3top, o2top, albedo, z, tlev,         & 
                          tlay, airlev, rh, xlwc, o3,            &
                          acb1, acb2, aoc1, aoc2, aant,          & 
                          aso4, asal, tauaer300, tauaer400, tauaer600, &
                          tauaer999, waer300, waer400, waer600,  &
                          waer999, gaer300, gaer400, gaer600,    &
                          gaer999, aer_ra_feedback, prate, adjcoe, radfld )

      use module_wave_data, only : nw, tuv_jmax, deltaw, sflx,   &
                                   c20, c40, c60, c80, sq
!-----------------------------------------------
! ... inter-face routine
!-----------------------------------------------
      implicit none
!----------------------------------------------------------
! ... dummy arguments
!----------------------------------------------------------
      integer, intent(in)    :: chem_opt      
      integer, intent(in)    :: nz      
      real(dp), intent(in)    :: zen     
      real(dp), intent(in)    :: o3toms  
      real(dp), intent(in)    :: esfact  
      real(dp), intent(in)    :: o3top   
      real(dp), intent(in)    :: o2top  
      real(dp), intent(in)    :: albedo(nw-1)
      real(dp), intent(in)    :: z(nz)
      real(dp), intent(in)    :: tlev(nz)
      real(dp), intent(in)    :: tlay(nz-1)
      real(dp), intent(in)    :: airlev(nz)
      real(dp), intent(in)    :: rh(nz)
      real(dp), intent(in)    :: xlwc(nz)
      real(dp), intent(in)    :: o3(nz)
      real(dp), intent(in)    :: acb1(nz)
      real(dp), intent(in)    :: acb2(nz)
      real(dp), intent(in)    :: aoc1(nz)
      real(dp), intent(in)    :: aoc2(nz)
      real(dp), intent(in)    :: aant(nz)
      real(dp), intent(in)    :: aso4(nz)
      real(dp), intent(in)    :: asal(nz)
! rajesh: declare aerosol optical properties
      real(dp), intent(in)    :: tauaer300(nz-1)
      real(dp), intent(in)    :: tauaer400(nz-1)
      real(dp), intent(in)    :: tauaer600(nz-1)
      real(dp), intent(in)    :: tauaer999(nz-1)
      real(dp), intent(in)    :: waer300(nz-1)
      real(dp), intent(in)    :: waer400(nz-1)
      real(dp), intent(in)    :: waer600(nz-1)
      real(dp), intent(in)    :: waer999(nz-1)
      real(dp), intent(in)    :: gaer300(nz-1)
      real(dp), intent(in)    :: gaer400(nz-1)
      real(dp), intent(in)    :: gaer600(nz-1)
      real(dp), intent(in)    :: gaer999(nz-1)
      INTEGER, INTENT(IN)     :: aer_ra_feedback

      real(dp), intent(out)   :: prate(nz,tuv_jmax)
      real(dp), intent(out)   :: adjcoe(nz,tuv_jmax)
      real(dp), intent(out)   :: radfld(nz,nw-1)
!----------------------------------------------------------
! ... local variables
!----------------------------------------------------------
      integer :: i, j, k, iw, n
      real(dp) :: factor, delzint
      character(len=8 ) :: cdate(4)
      character(len=10) :: ctime(4)
!----------------------------------------------------------
! ... altitude grid
!----------------------------------------------------------
      integer  :: iz, izi
      real(dp) :: colinc(nz)
      real(dp) :: vcol(nz)
      real(dp) :: scol(nz)
      real(dp) :: to3(nz)
!----------------------------------------------------------
! ... fast-j adj
!----------------------------------------------------------
      real(dp), dimension(nz,tuv_jmax) :: adjcoe1, adjcoe2
!----------------------------------------------------------
! ... solar zenith angle
! slant pathlengths in spherical geometry
!----------------------------------------------------------
      integer  :: nid(0:nz-1)
      real(dp) :: dsdh(0:nz-1,nz-1)
!----------------------------------------------------------
! ... extra terrestrial solar flux and earth-sun distance ^-2
!----------------------------------------------------------
      real(dp), dimension(nw-1) :: etf, delw, xsec
!--------------------------------------------------------------
! ... o2 absorption cross section
!--------------------------------------------------------------
      real(dp) :: xso2(nw-1,nz)
!--------------------------------------------------------------
! ... atmospheric optical parameters:
!--------------------------------------------------------------
      real(dp), dimension(nz-1,nw-1) :: dtrl, dto3, dto2
      real(dp), dimension(nz-1,nw-1) :: dtcld, omcld, gcld
      real(dp), dimension(nz-1,nw-1) :: dtcb1, omcb1, gcb1
      real(dp), dimension(nz-1,nw-1) :: dtcb2, omcb2, gcb2
      real(dp), dimension(nz-1,nw-1) :: dtoc1, omoc1, goc1
      real(dp), dimension(nz-1,nw-1) :: dtoc2, omoc2, goc2
      real(dp), dimension(nz-1,nw-1) :: dtant, omant, gant
      real(dp), dimension(nz-1,nw-1) :: dtso4, omso4, gso4
      real(dp), dimension(nz-1,nw-1) :: dtsal, omsal, gsal
! rajesh: declare optical property arrays for aerosols
      real(dp), dimension(nz-1,nw-1) :: dtaer, omaer, gaer
!--------------------------------------------------------------
! ... spectral irradiance and actinic flux (scalar irradiance):
!--------------------------------------------------------------
      real(dp), dimension(nz,nw-1) :: radxx
!-------------------------------------------------------------
! ... j-values:
!-------------------------------------------------------------
      integer :: m
!-------------------------------------------------------------
! ... location and time
!-------------------------------------------------------------
      integer :: iyear, imonth, iday
      real(dp) :: dtime, ut0
!-------------------------------------------------------------
! ... earth-sun distance effect
!-------------------------------------------------------------
      etf(1:nw-1) = sflx(1:nw-1) * esfact 
!-------------------------------------------------------
! ... air profile and rayleigh optical depths (inter-face)
!-------------------------------------------------------
      call setair( nz, z, airlev, dtrl, colinc, o2top )
!-------------------------------------------------------------
! ... ozone optical depths (must give temperature) (inter-face)
!-------------------------------------------------------------
      call setozo( nz, z, tlay, dto3, to3, o3, airlev, o3top, o3toms )
!-------------------------------------------------------------
! ... cloud optical depths (must cloud water g/m3) (inter-face)
!-------------------------------------------------------------
      call setcld( nz, z, xlwc, dtcld, omcld, gcld )

!-------------------------------------------------------------
! ... aerosol optical depths ...
!-------------------------------------------------------------
!      call setaer( nz, z, airlev, rh, acb1,        &
!                   acb2, aoc1, aoc2, aant, aso4,   &
!                   asal,                           &
!                   dtcb1, omcb1, gcb1,             &
!                   dtcb2, omcb2, gcb2,             &
!                   dtoc1, omoc1, goc1,             &
!                   dtoc2, omoc2, goc2,             &
!                   dtant, omant, gant,             &
!                   dtso4, omso4, gso4,             &
!                   dtsal, omsal, gsal )

!------------------------------------------------------------
! rajesh: Conform Aerosol Optical Properties from 4 wavelengths to 
! the entire spectra of FTUV
!------------------------------------------------------------
! Initialize aerosol optical properties to zero
      dtaer(:,:) = 0._dp
      omaer(:,:) = 0._dp
       gaer(:,:) = 0._dp
      if(aer_ra_feedback == 1) then 
       call aer_wrf2ftuv(nz, z, tauaer300, tauaer400,    &
                         tauaer600, tauaer999, waer300, &
                         waer400, waer600, waer999,     &
                         gaer300, gaer400, gaer600,     &
                         gaer999, dtaer, omaer, gaer    )
      endif 

!------------------------------------------------------------
! ... photo-chemical and photo-biological weigting functions.
! for pchem, need to know temperature and pressure profiles.
! output:
! from pchem: sq(nw,nz,nj) - for each reaction jlabel(nj)
!-------------------------------------------------------------
      call pchem( chem_opt, nz, tlev, airlev )

!-------------------------------------------------------------
! ... slant path lengths for spherical geometry
!-------------------------------------------------------------
       call spheres( nz, z, zen, dsdh, nid )

       call airmas( nz, z, zen, dsdh, nid, colinc, vcol, scol )

!---------------------------------------------------------------
! ... modification of coefficent of j-vales function of to3 and zenith
!---------------------------------------------------------------
      if( zen <= 20.5_dp ) then
         call setz( nz, to3, tlev, c20, 20, adjcoe )
      else if( zen > 20.5_dp .and. zen < 40.5_dp ) then
         call setz( nz, to3, tlev, c20, 20, adjcoe1 )
         call setz( nz, to3, tlev, c40, 40, adjcoe2 )
         factor = (zen - 20.5_dp) / 20._dp
         adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) &
                               + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor
      else if( zen >= 40.5_dp .and. zen < 60.5_dp ) then
         call setz( nz, to3, tlev, c40, 40, adjcoe1 )
         call setz( nz, to3, tlev, c60, 60, adjcoe2 )
         factor = (zen - 40.5_dp) / 20._dp
         adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) &
                               + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor
      else if( zen >= 60.5_dp .and. zen < 80._dp ) then
         call setz( nz, to3, tlev, c60, 60, adjcoe1 )
         call setz( nz, to3, tlev, c80, 80, adjcoe2 )
         factor = (zen - 60.5_dp) / 19.5_dp
         adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) &
                               + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor
      else if( zen >= 80. ) then
         call setz( nz, to3, tlev, c80, 80, adjcoe )
      end if
!------------------------------------------------------------------
! ... effective o2 optical depth (sr bands, must know zenith angle!)
! assign o2 cross section to sq(*,*,1)
!------------------------------------------------------------------
       call set_o2_xsect( nz, z, colinc, vcol, scol, dto2, xso2 )
       sq(1:nw-1,1:nz,1) = xso2(1:nw-1,1:nz)

!-----------------------------------------------
! ... main wavelength loop
!-----------------------------------------------
      delw(1:nw-1) = deltaw(1:nw-1) * etf(1:nw-1)
!---------------------------------------------------
! ... monochromatic radiative transfer:
! outputs are fdir, fdn, fup
!---------------------------------------------------
      call rtlink( nz,                  &
                   nw,                  &
                   zen,                 &
                   albedo,              & 
                   z,                   &
                   nid,                 &
                   dsdh,                & 
                   dtrl,                & ! opt depth of air (rayleigh)
                   dto3,                & ! opt depth of o3
                   dto2,                & ! opt depth of o2
                   dtcld, omcld, gcld,  & ! opt depth of cloud
!                   dtcb1, omcb1, gcb1,  &
!                   dtcb2, omcb2, gcb2,  &
!                   dtoc1, omoc1, goc1,  &
!                   dtoc2, omoc2, goc2,  &
!                   dtant, omant, gant,  &
!                   dtso4, omso4, gso4,  &
!                   dtsal, omsal, gsal,  &
                   dtaer, omaer, gaer,  &  ! pass aerosol optical properties   
                   radfld )                           ! output of abs. scart flux.
!----------------------------------------------------------
! Interplation the top level
!----------------------------------------------------------
         radfld(:,:) = max( radfld(:,:),0._dp )
         delzint = (z(nz-2) - z(nz-3))/(z(nz-1) - z(nz-2))
         do iw = 1,nw-1
            radfld(1,iw) = radfld(2,iw) + (radfld(2,iw)-radfld(3,iw))*delzint
            radfld(1,iw) = max(radfld(1,iw),radfld(2,iw))
         end do

!----------------------------------------------------------
! ... j-val calculation
! spherical irradiance (actinic flux)
! as a function of altitude
! convert to quanta s-1 nm-1 cm-2
! (1.e-4 * (wc*1e-9) / (hc = 6.62e-34 * 2.998e8))
!----------------------------------------------------------
      do m = 1,tuv_jmax
        do iz = 1,nz
          izi = nz - iz + 1
          xsec(:nw-1) = sq(:nw-1,iz,m) * delw(:nw-1)
          prate(iz,m) = dot_product( radfld(izi,:nw-1), xsec(:nw-1) )
        end do
        prate(1:nz,m) = prate(1:nz,m) * adjcoe(1:nz,m)
      end do

      end subroutine photoin

!----------------------------------------------------------------------------
!rajesh: subroutine to convert aerosol optical properties from 4
!wavelengths 
! to entire spectra of FTUV
!-----------------------------------------------------------------------------
      subroutine aer_wrf2ftuv( nzlev, z, tauaer300, tauaer400,     &
                               tauaer600, tauaer999,               &
                               waer300, waer400, waer600, waer999, &
                               gaer300, gaer400, gaer600, gaer999, &
                               dtaer, omaer, gaer )
       use module_wave_data, only : nw, wc

!----------------------------------------------------------------------
! The routine is based on aerosol treatment in module_ra_rrtmg_sw.F
! INPUT: 
! nzlev: number of specified altitude levels in the working grid
! z: specified altitude working grid   
! Aerosol optical properties at 300, 400, 600 and 999 nm. 
!   tauaer300, tauaer400, tauaer600, tauaer999: Layer AODs
!   waer300, waer400, waer600, waer999: Layer SSAs
!   gaer300, gaer400, gaer600, gaer999: Layer asymmetry parameters

! OUTPUT:
! dtaer: Layer AOD at FTUV wavelengths
! omaer: Layer SSA at FTUV wavelengths
! gaer : Layer asymmetry parameters at FTUV wavelengths
!------------------------------------------------------------------------
      implicit none

!-----------------------------------------------------------------------------
!       ... Dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in)   :: nzlev
      real(dp), intent(in)  :: z(nzlev)
      real(dp), intent(in)  :: tauaer300(nzlev-1), tauaer400(nzlev-1),    &
                               tauaer600(nzlev-1), tauaer999(nzlev-1)
      real(dp), intent(in)  :: waer300(nzlev-1), waer400(nzlev-1),        &
                               waer600(nzlev-1), waer999(nzlev-1)
      real(dp), intent(in)  :: gaer300(nzlev-1), gaer400(nzlev-1),        &
                               gaer600(nzlev-1), gaer999(nzlev-1)
! Output arrays
      real(dp), intent(out) :: dtaer(nzlev-1,nw-1),                   &
                               omaer(nzlev-1,nw-1), gaer(nzlev-1,nw-1)

! Local Variables
      integer     :: k, wn, nloop
      real(dp)    :: ang, slope
      real(dp), parameter :: thresh = 1.e-9_dp

      ang   = 0._dp
      slope = 0._dp

! Start Calculation
      do wn = 1,nw-1       ! wavelength loop
       do k = 1,nzlev-1    ! level loop 

! use angstrom exponent to calculate aerosol optical depth; wc is in nm.  
        if(tauaer300(k) .gt. thresh .and. tauaer999(k) .gt. thresh) then
         ang = log(tauaer300(k)/tauaer999(k))/log(0.999_dp/0.3_dp)
         dtaer(k,wn) = tauaer400(k)*(0.4_dp/(wc(wn)*1.e-3_dp))**ang
!         print *, ang, dtaer(k,wn), tauaer600(k)*(600./wc(wn))**ang

! ssa - use linear interpolation/extrapolation
         slope = (waer600(k)-waer400(k))/0.2_dp
         omaer(k,wn) = slope*((wc(wn)*1.e-3_dp)-0.6_dp)+waer600(k)
         if(omaer(k,wn) .lt. 0.4_dp) omaer(k,wn)=0.4_dp
         if(omaer(k,wn) .ge. 1.0_dp) omaer(k,wn)=1.0_dp

! asymmetry parameter - use linear interpolation/extrapolation
         slope = (gaer600(k)-gaer400(k))/0.2_dp
         gaer(k,wn) = slope*((wc(wn)*1.e-3_dp)-0.6_dp)+gaer600(k)
         if(gaer(k,wn) .lt. 0.5_dp) gaer(k,wn) = 0.5_dp
         if(gaer(k,wn) .ge. 1.0_dp) gaer(k,wn) = 1.0_dp
        endif
       enddo  ! k 
      enddo   ! wn

      end subroutine aer_wrf2ftuv

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

      subroutine setaer( nzlev, z, airden, rh, acb1,    &
                         acb2, aoc1, aoc2, aant, aso4,  &
                         asal,                          &
                         op_cb1, om_cb1, g_cb1,         &
                         op_cb2, om_cb2, g_cb2,         &
                         op_oc1, om_oc1, g_oc1,         &
                         op_oc2, om_oc2, g_oc2,         &
                         op_ant, om_ant, g_ant,         &
                         op_so4, om_so4, g_so4,         &
                         op_sal, om_sal, g_sal )

      use module_wave_data, only : nw, wc

!-----------------------------------------------------------------------------
!=  PURPOSE:                                                                 
!
!=  Cal aer opt depth
!-----------------------------------------------------------------------------
!
!=  PARAMETERS:                                                              
!=  NZLEV   - INTEGER, number of specified altitude levels in the working (I)
!=            grid                                                           
!=  Z       - real(dp), specified altitude working grid (km)                  (I)
!=  axxx    - aerosol mix ratio                                           (I)
!
!=  op      - real(dp), optical depth (absorption)                            (O)
!=  om      - real(dp), single albedo             
!=  g       - asysmetric factor                                           (O)
!
!-----------------------------------------------------------------------------
!
!     VERTICAL DOMAIN is from bottom(1)  to TOP (TOP=nzlev)
!        CCM from top(1) to bottom(nzlev)
!-----------------------------------------------------------------------------
      implicit none

!-----------------------------------------------------------------------------
!	... Dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in)   :: nzlev
      real(dp), intent(in)  :: z(nzlev)
      real(dp), intent(in)  :: airden(nzlev)
      real(dp), intent(in)  :: rh(nzlev)
      real(dp), intent(in)  :: acb1(nzlev)
      real(dp), intent(in)  :: acb2(nzlev)
      real(dp), intent(in)  :: aoc1(nzlev)
      real(dp), intent(in)  :: aoc2(nzlev)
      real(dp), intent(in)  :: aant(nzlev)
      real(dp), intent(in)  :: aso4(nzlev)
      real(dp), intent(in)  :: asal(nzlev)
      real(dp), intent(out) :: op_cb1(nzlev-1,nw-1), om_cb1(nzlev-1,nw-1), g_cb1(nzlev-1,nw-1)
      real(dp), intent(out) :: op_cb2(nzlev-1,nw-1), om_cb2(nzlev-1,nw-1), g_cb2(nzlev-1,nw-1)
      real(dp), intent(out) :: op_oc1(nzlev-1,nw-1), om_oc1(nzlev-1,nw-1), g_oc1(nzlev-1,nw-1)
      real(dp), intent(out) :: op_oc2(nzlev-1,nw-1), om_oc2(nzlev-1,nw-1), g_oc2(nzlev-1,nw-1)
      real(dp), intent(out) :: op_ant(nzlev-1,nw-1), om_ant(nzlev-1,nw-1), g_ant(nzlev-1,nw-1)
      real(dp), intent(out) :: op_so4(nzlev-1,nw-1), om_so4(nzlev-1,nw-1), g_so4(nzlev-1,nw-1)
      real(dp), intent(out) :: op_sal(nzlev-1,nw-1), om_sal(nzlev-1,nw-1), g_sal(nzlev-1,nw-1)

!-----------------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------------
      integer     :: k, wn, nz
      real(dp)    :: wcen
      real(dp)    :: dz(nzlev)
      real(dp)    :: sig_oc2(nzlev-1), sig_cb2(nzlev-1)
      real(dp)    :: opt_oc1(nzlev-1), opt_oc2(nzlev-1), opt_cb1(nzlev-1), opt_cb2(nzlev-1)
      real(dp)    :: opt_ant(nzlev-1), opt_so4(nzlev-1), opt_sal(nzlev-1)  
      real(dp)    :: rw(nzlev-1), num(nzlev-1), wght(nzlev-1)


!-----------------------------------------------------------------------------
! 	... assign the constants
!-----------------------------------------------------------------------------
      real(dp), parameter :: rm_cb1 = 0.035_dp
      real(dp), parameter :: rm_cb2 = 0.035_dp
      real(dp), parameter :: rm_oc1 = 0.10_dp
      real(dp), parameter :: rm_oc2 = 0.10_dp
      real(dp), parameter :: rm_ant = 0.15_dp
      real(dp), parameter :: rm_so4 = 0.24_dp
      real(dp), parameter :: rm_sal = 1.30_dp

      real(dp), parameter :: de_cb1 = 1.0_dp
      real(dp), parameter :: de_cb2 = 1.0_dp
      real(dp), parameter :: de_oc1 = 1.5_dp 
      real(dp), parameter :: de_oc2 = 1.5_dp 
      real(dp), parameter :: de_ant = 1.7_dp  
      real(dp), parameter :: de_so4 = 1.7_dp
      real(dp), parameter :: de_sal = 2.2_dp

      real(dp), parameter :: ml_cb1 = 12._dp
      real(dp), parameter :: ml_cb2 = 12._dp
      real(dp), parameter :: ml_oc1 = 48._dp 
      real(dp), parameter :: ml_oc2 = 48._dp 
      real(dp), parameter :: ml_ant = 90._dp  
      real(dp), parameter :: ml_so4 = 96._dp
      real(dp), parameter :: ml_sal = 22._dp

      real(dp), parameter :: qx_cb1 = 2.138_dp
      real(dp), parameter :: qx_cb2 = 2.138_dp
      real(dp), parameter :: qw_cb1 = 0.510_dp
      real(dp), parameter :: qw_cb2 = 0.510_dp
      real(dp), parameter :: sx_cb2 = 0.120_dp
      real(dp), parameter :: sw_cb2 = 1.000_dp

      real(dp), parameter :: qx_oc1 = 0.675_dp
      real(dp), parameter :: qx_oc2 = 0.675_dp
      real(dp), parameter :: qw_oc1 = 1.118_dp
      real(dp), parameter :: qw_oc2 = 1.118_dp
      real(dp), parameter :: sx_oc2 = 0.980_dp
      real(dp), parameter :: sw_oc2 = 1.000_dp

      real(dp), parameter :: qx_ant = 0.726_dp
      real(dp), parameter :: qw_ant = 1.770_dp

      real(dp), parameter :: qx_so4 = 2.830_dp
      real(dp), parameter :: qw_so4 = 3.605_dp

      real(dp), parameter :: qx_sal = 2.600_dp
      real(dp), parameter :: qw_sal = 2.575_dp

      real(dp), parameter :: pi = 3.1416_dp

!-----------------------------------------------------------------------------
! 	... hydroscropic growth
!-----------------------------------------------------------------------------
      real(dp), parameter :: a1_cb2 = 1.000_dp
      real(dp), parameter :: a2_cb2 =-0.829_dp
      real(dp), parameter :: a3_cb2 = 1.350_dp

      real(dp), parameter :: a1_oc2 = 1.005_dp
      real(dp), parameter :: a2_oc2 =-0.0666_dp
      real(dp), parameter :: a3_oc2 = 0.7608_dp

      real(dp), parameter :: a1_ant = 1.000_dp
      real(dp), parameter :: a2_ant = 0.493_dp
      real(dp), parameter :: a3_ant = 0.527_dp

      real(dp), parameter :: a1_so4 = 1.000_dp
      real(dp), parameter :: a2_so4 = 0.493_dp
      real(dp), parameter :: a3_so4 = 0.528_dp

      real(dp), parameter :: a1_sal = 1.011_dp
      real(dp), parameter :: a2_sal = 0.457_dp
      real(dp), parameter :: a3_sal = 1.102_dp

      nz = nzlev - 1
!-----------------------------------------------------------------------------
! 	... calculate vertical interveral dz
!-----------------------------------------------------------------------------    
      dz(1:nz) = abs( z(2:nz+1) - z(1:nz) ) * km2cm * pi * 1.e-8_dp
!-----------------------------------------------------------------------------
!        -------- CB1 -------------
!-----------------------------------------------------------------------------
         rw(:nz) = rm_cb1
         call cnum( acb1, airden, de_cb1, ml_cb1, rm_cb1, &
                    rw, num, wght, nzlev )
      
         opt_cb1(:nz) = qx_cb1*wght(:nz) + (1._dp - wght(:nz))*qw_cb1     
         opt_cb1(:nz) = opt_cb1(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
!-----------------------------------------------------------------------------
!        -------- CB2 -------------
!-----------------------------------------------------------------------------
         rw(:nz) = rm_cb2*(a1_cb2 + rh(:nz)*(a2_cb2 + a3_cb2*rh(:nz)))
         call cnum( acb2, airden, de_cb2, ml_cb2, rm_cb2, &
                    rw, num, wght, nzlev )
      
         opt_cb2(:nz) = qx_cb2*wght(:nz) + (1._dp - wght(:nz))*qw_cb2     
         opt_cb2(:nz) = opt_cb2(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
         sig_cb2(:nz) = sx_cb2*wght(:nz) + (1._dp - wght(:nz))*sw_cb2   
!-----------------------------------------------------------------------------
!        -------- OC1 -------------
!-----------------------------------------------------------------------------
         rw(:nz) = rm_oc1
         call cnum( aoc1, airden, de_oc1, ml_oc1, rm_oc1, &
                    rw, num, wght, nzlev )
      
         opt_oc1(:nz) = qx_oc1*wght(:nz) + (1._dp - wght(:nz))*qw_oc1     
         opt_oc1(:nz) = opt_oc1(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
!-----------------------------------------------------------------------------
!        -------- OC2 -------------
!-----------------------------------------------------------------------------
         rw(:nz) = rm_oc2*(a1_oc2 + rh(:nz)*(a2_oc2 + a3_oc2*rh(:nz)))
         call cnum( aoc2, airden, de_oc2, ml_oc2, rm_oc2, &
                    rw, num, wght, nzlev )
      
         opt_oc2(:nz) = qx_oc2*wght(:nz) + (1._dp - wght(:nz))*qw_oc2     
         opt_oc2(:nz) = opt_oc2(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
         sig_oc2(:nz) = sx_oc2*wght(:nz) + (1._dp - wght(:nz))*sw_oc2   
!-----------------------------------------------------------------------------
!        -------- ANT-NH4NO3 -------------
!-----------------------------------------------------------------------------
         rw(:nz) = rm_ant*(a1_ant + rh(:nz)*(a2_ant + a3_ant*rh(:nz)))
         call cnum( aant, airden, de_ant, ml_ant, rm_ant, &
                    rw, num, wght, nzlev )
      
         opt_ant(:nz) = qx_ant*wght(:nz) + (1._dp - wght(:nz))*qw_ant     
         opt_ant(:nz) = opt_ant(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
!-----------------------------------------------------------------------------
!        -------- SO4 -------------
!-----------------------------------------------------------------------------
         rw(:nz) = rm_so4*(a1_so4 + rh(:nz)*(a2_so4 + a3_so4*rh(:nz)))
         call cnum( aso4, airden, de_so4, ml_so4, rm_so4, &
                    rw, num, wght, nzlev )
      
         opt_so4(:nz) = qx_so4*wght(:nz) + (1._dp - wght(:nz))*qw_so4     
         opt_so4(:nz) = opt_so4(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
!-----------------------------------------------------------------------------
!        -------- SALT -------------
!-----------------------------------------------------------------------------
         rw(:nz) = rm_sal*(a1_sal + rh(:nz)*(a2_sal + a3_sal*rh(:nz)))
         call cnum( asal, airden, de_sal, ml_sal, rm_sal, &
                    rw, num, wght, nzlev )
      
         opt_sal(:nz) = qx_sal*wght(:nz) + (1._dp - wght(:nz))*qw_sal     
         opt_sal(:nz) = opt_sal(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)

WAVE_LENGTH_LOOP : &
         do wn = 1,nw-1
            wcen = wc(wn)
            op_cb1(:nz,wn) = opt_cb1(:nz)*op_cb1w(wn)
            om_cb1(:nz,wn) = om_cb1w(wn)
             g_cb1(:nz,wn) = g_cb1w(wn)

            op_cb2(:nz,wn) = opt_cb2(:nz)*op_cb2w(wn)
            om_cb2(:nz,wn) = sig_cb2(:nz)*om_cb2w(wn)
             g_cb2(:nz,wn) = g_cb2w(wn)

            op_oc1(:nz,wn) = opt_oc1(:nz)*op_oc1w(wn)
            om_oc1(:nz,wn) = 0.98_dp
             g_oc1(:nz,wn) = g_oc1w(wn)

            op_oc2(:nz,wn) = opt_oc2(:nz)*op_oc2w(wn)
            om_oc2(:nz,wn) = sig_oc2(:nz)
             g_oc2(:nz,wn) = g_oc2w(wn)

            op_ant(:nz,wn) = opt_ant(:nz)*op_antw(wn)
            om_ant(:nz,wn) = 1.00_dp
             g_ant(:nz,wn) = g_antw(wn)

            op_so4(:nz,wn) = opt_so4(:nz)*op_so4w(wn)
            om_so4(:nz,wn) = 1.00_dp
             g_so4(:nz,wn) = 0.75_dp

            op_sal(:nz,wn) = opt_sal(:nz)*op_salw(wn)
            om_sal(:nz,wn) = 1.00_dp
             g_sal(:nz,wn) = 0.77_dp
         enddo WAVE_LENGTH_LOOP

      end subroutine setaer

      subroutine aer_init
!-----------------------------------------------------------------------------
!   ... initialize setaer
!-----------------------------------------------------------------------------

      use module_wave_data, only : nw, wc

      implicit none

!-----------------------------------------------------------------------
!     Assign the wave-length dep values for aerososl
!-----------------------------------------------------------------------
      real(dp), parameter :: d1_cb = 1.61892_dp
      real(dp), parameter :: d2_cb =-0.0006853_dp
      real(dp), parameter :: d3_cb =-1.07806e-6_dp

      real(dp), parameter :: o1_cb = 1.29185_dp
      real(dp), parameter :: o2_cb = 7.6863e-5_dp
      real(dp), parameter :: o3_cb =-1.28567e-6_dp

      real(dp), parameter :: g1_cb = 2.87326_dp
      real(dp), parameter :: g2_cb =-0.004482_dp
      real(dp), parameter :: g3_cb = 1.50361e-6_dp

      real(dp), parameter :: d1_oc = 10.5098_dp
      real(dp), parameter :: d2_oc =-0.0301455_dp
      real(dp), parameter :: d3_oc = 2.22837e-5_dp

      real(dp), parameter :: g1_oc = 1.67125_dp
      real(dp), parameter :: g2_oc =-0.0008341_dp
      real(dp), parameter :: g3_oc =-1.18801e-6_dp

      real(dp), parameter :: d1_ant = 9.44794_dp
      real(dp), parameter :: d2_ant =-0.0268029_dp
      real(dp), parameter :: d3_ant = 1.97106e-5_dp

      real(dp), parameter :: g1_ant = 1.34742_dp
      real(dp), parameter :: g2_ant =-1.85850e-5_dp
      real(dp), parameter :: g3_ant =-1.54297e-6_dp

      real(dp), parameter :: d1_so4 = 0.942351_dp
      real(dp), parameter :: d2_so4 = 0.00220168_dp
      real(dp), parameter :: d3_so4 =-3.9622e-6_dp

      real(dp), parameter :: d1_sal = 1.0_dp
      real(dp), parameter :: d2_sal = 0.0_dp
      real(dp), parameter :: d3_sal = 0.0_dp

      op_cb1w(:nw-1) = d1_cb + wc(:nw-1)*(d2_cb + d3_cb*wc(:nw-1))
      om_cb1w(:nw-1) = .12_dp*(o1_cb + wc(:nw-1)*(o2_cb + o3_cb*wc(:nw-1)))
      g_cb1w(:nw-1)  = .12_dp*(g1_cb + wc(:nw-1)*(g2_cb + g3_cb*wc(:nw-1)))

      op_cb2w(:nw-1) = op_cb1w(:nw-1)
      om_cb2w(:nw-1) = om_cb1w(:nw-1)
      g_cb2w(:nw-1)  = g_cb1w(:nw-1)

      op_oc1w(:nw-1) = d1_oc + wc(:nw-1)*(d2_oc + d3_oc*wc(:nw-1))
      g_oc1w(:nw-1)  = 0.50_dp*(g1_oc + wc(:nw-1)*(g2_oc + g3_oc*wc(:nw-1)))

      op_oc2w(:nw-1) = op_oc1w(:nw-1)
!     om_oc2w(:nw-1) = sig_oc2
      g_oc2w(:nw-1)  = g_oc1w(:nw-1)

      op_antw(:nw-1) = d1_ant + wc(:nw-1)*(d2_ant + d3_ant*wc(:nw-1))
      g_antw(:nw-1)  = 0.63_dp*(g1_ant + wc(:nw-1)*(g2_ant + g3_ant*wc(:nw-1)))

      op_so4w(:nw-1) = d1_so4 + wc(:nw-1)*(d2_so4 + d3_so4*wc(:nw-1))

      op_salw(:nw-1) = d1_sal + wc(:nw-1)*(d2_sal + d3_sal*wc(:nw-1))

      end subroutine aer_init

      subroutine cnum( mix, denx, cden, molw, rm, &
                       rw, num, weight, nzlev )
!-----------------------------------------------------------------------------
!   PURPOSE:                                                                 
!
!   Cal aer num and percent weigh
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
!	... Dummy arguments
!-----------------------------------------------------------------------------
      integer,  intent(in)  :: nzlev
      real(dp), intent(in)  :: molw, rm
      real(dp), intent(in)  :: cden
      real(dp), intent(in)  :: rw(nzlev-1)
      real(dp), intent(in)  :: mix(nzlev-1), denx(nzlev-1)
      real(dp), intent(out) :: num(nzlev-1), weight(nzlev-1)
!-----------------------------------------------------------------------------
!	... Local variables
!-----------------------------------------------------------------------------
      REAL(dp), PARAMETER :: PI = 3.14159265358979324
      REAL(dp), PARAMETER :: vol_fac = 4._dp * PI / 3._dp
      integer  :: k
      real(dp) :: rm_cubic, factor
      real(dp) :: mas(nzlev-1), masw(nzlev-1)

      rm_cubic = rm**3
      factor   = 1._dp/(vol_fac * cden * 1.e-12_dp * rm_cubic)
!-----------------------------------------------------------------------------
!     Convert volumn mix ratio to mass
!-----------------------------------------------------------------------------
      where( mix(:nzlev-1) >= 1.e-15_dp )
!===============================================================================
! XUEXI
!     UNLIKE MOZART mix = mix ratio
!                   mix = ug/m3
!     1e-12 for  ug/m3 --> gram/cm3 
!===============================================================================
         mas(:nzlev-1)    = mix(:nzlev-1) * 1.e-12_dp                              ! 1e-12 for  ug/m3 --> gram/cm3 
         num(:nzlev-1)    = mas(:nzlev-1) * factor
      elsewhere
        num(:nzlev-1)    = 0._dp
        weight(:nzlev-1) = 0._dp
      endwhere

      do k = 1,nzlev-1
         if( abs(rw(k) - rm) >= 1.e-6_dp ) then
           masw(k) = vol_fac*(rw(k)**3 - rm_cubic)*1.e-12_dp * num(k)   ! water mass gram/cm3
         else
           masw(k) = 0._dp
         end if
      end do

      where( mix(:nzlev-1) >= 1.e-15_dp )
         weight(:nzlev-1) = mas(:nzlev-1)/( mas(:nzlev-1) + masw(:nzlev-1) )        ! weight percentage (%*0.01)
         weight(:nzlev-1) = max( 0._dp, min( 1._dp,weight(:nzlev-1) ) )
      endwhere

      end subroutine cnum

      subroutine setair( nz, z, airlev, dtrl, &
                         cz, o2top )

      use module_wave_data, only : nw, wc
!-----------------------------------------------------------------------------
! purpose:
! set up an altitude profile of air molecules. subroutine includes a
! shape-conserving scaling method that allows scaling of the entire
! profile to a given sea-level pressure.
!-----------------------------------------------------------------------------
! parameters:
! pmbnew - real(dp), sea-level pressure (mb) to which profile should be
! scaled. if pmbnew < 0, no scaling is done
! nz - integer, number of specified altitude levels in the working (i)
! grid
! z - real(dp), specified altitude working grid (km) (i)
! nw - integer, number of specified intervals + 1 in working (i)
! wavelength grid
! wl - real(dp), vector of lower limits of wavelength intervals in (i)
! working wavelength grid
! airlev - real(dp), air density (molec/cc) at each specified altitude (o)
! dtrl - real(dp), rayleigh optical depth at each specified altitude (o)
! and each specified wavelength
! cz - real(dp), number of air molecules per cm^2 at each specified (o)
! altitude layer
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      real(dp), intent(in) :: o2top
      real(dp), intent(in) :: z(nz)
      real(dp), intent(in) :: airlev(nz)
!-----------------------------------------------------------------------------
! ... air density (molec cm-3) at each grid level
! rayleigh optical depths
!-----------------------------------------------------------------------------
      real(dp), intent(out) :: dtrl(nz-1,nw-1)
      real(dp), intent(out) :: cz(nz)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: i, ip1, iw
      real(dp) :: hscale
      real(dp) :: srayl
      real(dp) :: deltaz
      real(dp) :: wmicrn, xx
!-----------------------------------------------------
! ... compute column increments (logarithmic integrals)
!-----------------------------------------------------
      do i = 1,nz-1
         ip1 = i + 1
         deltaz = km2cm * (z(ip1) - z(i))
         cz(i) = (airlev(ip1) - airlev(i)) / log(airlev(ip1)/airlev(i)) * deltaz
      end do
!-----------------------------------------------------
! ... include exponential tail integral from infinity to 50 km,
! fold tail integral into top layer
! specify scale height near top of data. (scale height at 40km????)
!-----------------------------------------------------
! hscale = 8.05e5
      cz(nz) = o2top/.2095_dp
!-----------------------------------------------------
! ... compute rayleigh cross sections and depths:
!-----------------------------------------------------
      do iw = 1,nw-1
!-----------------------------------------------------
! ... rayleigh scattering cross section from wmo 1985 (originally from
! nicolet, m., on the molecular scattering in the terrestrial atmosphere:
! an empirical formula for its calculation in the homoshpere, planet.
! space sci., 32, 1467-1468, 1984.
!-----------------------------------------------------
         wmicrn = 1.e-3_dp*wc(iw)
         if( wmicrn <= .55_dp ) then
            xx = 3.6772_dp + 0.389_dp*wmicrn + 0.09426_dp/wmicrn
         else
            xx = 4.04_dp
         end if
         srayl = 4.02e-28_dp/(wmicrn)**xx
         dtrl(1:nz-2,iw) = cz(1:nz-2)*srayl
         dtrl(nz-1,iw) = (cz(nz-1) + cz(nz))*srayl
      end do

      end subroutine setair

      subroutine zenith( lat, long, idate, ut, zen )
!-----------------------------------------------------------------------------
! purpose:
! calculate solar zenith angle for a given time and location.
! calculation is based on equations given in: paltridge and platt, radia-
! tive processes in meteorology and climatology, elsevier, pp. 62,63, 1976.
! fourier coefficients originally from: spencer, j.w., 1971, fourier
! series representation of the position of the sun, search, 2:172.
! note: this approximate program does not account fro changes from year
! to year.
!-----------------------------------------------------------------------------
! parameters:
! lat - real(dp), latitude of location (degrees) (i)
! long - real(dp), longitude of location (degrees) (i)
! idate - integer, date in the form yymmdd (i)
! ut - real(dp), local time in decimal ut (e.g., 16.25 means 15 minutes (i)
! after 4 pm)
! zen - real(dp), solar zenith angle (degrees) (o)
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: idate
      real(dp), intent(in) :: lat,long
      real(dp), intent(in) :: ut
      real(dp), intent(out) :: zen
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      real(dp), parameter :: pi = 3.1415926
      real(dp), parameter :: d2r = pi/180.0
      real(dp), parameter :: r2d = 1.0/d2r
      integer :: i
      integer :: iiyear, imth, iday, ijd
      integer :: imn(12) = (/ 31,28,31,30,31,30,31,31,30,31,30,31 /)
      real(dp) :: lbut,lzut
      real(dp) :: rlt
      real(dp) :: d, tz, rdecl, eqr, eqh, zpt
      real(dp) :: csz, zr, caz, raz
      real(dp) :: sintz, costz, sin2tz, cos2tz, sin3tz, cos3tz
!-----------------------------------------------------------------------------
! ... convert to radians
!-----------------------------------------------------------------------------
      rlt = lat*d2r
!-----------------------------------------------------------------------------
! ... parse date
!-----------------------------------------------------------------------------
      iiyear = idate/10000
      imth = (idate - iiyear*10000)/100
      iday = idate - iiyear*10000 - imth*100
!-----------------------------------------------------------------------------
! ... identify and correct leap years
!-----------------------------------------------------------------------------
      if( mod(iiyear,4) == 0 ) then
         imn(2) = 29
      else
         imn(2) = 28
      end if
!-----------------------------------------------------------------------------
! ... compute current (julian) day of year ijd = 1 to 365
!-----------------------------------------------------------------------------
      ijd = 0
      do i = 1,imth-1
         ijd = ijd + imn(i)
      end do
      ijd = ijd + iday
!-----------------------------------------------------------------------------
! ... calculate decimal julian day from start of year:
!-----------------------------------------------------------------------------
      d = real(ijd-1,kind=dp) + ut/24._dp
!-----------------------------------------------------------------------------
! ... equation 3.8 for "day-angle"
!-----------------------------------------------------------------------------
      tz = 2._dp*pi*d/365._dp
!-----------------------------------------------------------------------------
! ... calculate sine and cosine from addition theoremes for
! better performance; the computation of sin2tz,
! sin3tz, cos2tz and cos3tz is about 5-6 times faster
! than the evaluation of the intrinsic functions
!
! it is sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y)
! and cos(x+y) = cos(x)*cos(y)-sin(x)*sin(y)
!
! sintz = sin(tz) costz = cos(tz)
! sin2tz = sin(2.*tz) cos2tz = sin(2.*tz)
! sin3tz = sin(3.*tz) cos3tz = cos(3.*tz)
!-----------------------------------------------------------------------------
      sintz = sin( tz )
      costz = cos( tz )
      sin2tz = 2._dp*sintz*costz
      cos2tz = costz*costz - sintz*sintz
      sin3tz = sintz*cos2tz + costz*sin2tz
      cos3tz = costz*cos2tz - sintz*sin2tz
!-----------------------------------------------------------------------------
! ... equation 3.7 for declination in radians
!-----------------------------------------------------------------------------
      rdecl = 0.006918_dp - 0.399912_dp*costz + 0.070257_dp*sintz &
                          - 0.006758_dp*cos2tz + 0.000907_dp*sin2tz &
                          - 0.002697_dp*cos3tz + 0.001480_dp*sin3tz
!-----------------------------------------------------------------------------
! ... equation 3.11 for equation of time in radians
!-----------------------------------------------------------------------------
      eqr = 0.000075_dp + 0.001868_dp*costz - 0.032077_dp*sintz &
                        - 0.014615_dp*cos2tz - 0.040849_dp*sin2tz
!-----------------------------------------------------------------------------
! ... convert equation of time to hours
!-----------------------------------------------------------------------------
      eqh = eqr*24._dp/(2._dp*pi)
!-----------------------------------------------------------------------------
! ... calculate local hour angle (hours):
!-----------------------------------------------------------------------------
      lbut = 12._dp - eqh - long*24._dp/360._dp
!-----------------------------------------------------------------------------
! ... convert to angle from ut
!-----------------------------------------------------------------------------
      lzut = 15._dp*(ut - lbut)
      zpt = lzut*d2r
!-----------------------------------------------------------------------------
! ... equation 2.4 for cosine of zenith angle
!-----------------------------------------------------------------------------
      csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
      zr = acos(csz)
      zen = zr*r2d

      end subroutine zenith

      subroutine calc_zenith( lat, long, ijd, gmt, azimuth, zenith )
!-------------------------------------------------------------------
! this subroutine calculates solar zenith and azimuth angles for a
! part  time and location.  must specify:
! input:
! lat - latitude in decimal degrees
! long - longitude in decimal degrees
! gmt  - greenwich mean time - decimal military eg.
! 22.75 = 45 min after ten pm gmt
! output:
! zenith
! azimuth
! .. Scalar Arguments ..
!-------------------------------------------------------------------
        integer, intent(in)   :: ijd
        real(dp), intent(in)  :: gmt, lat, long
        real(dp), intent(out) :: azimuth, zenith

!-------------------------------------------------------------------
! .. Local variables
!-------------------------------------------------------------------
      real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp
      real(dp), parameter :: r2d = 1.0_dp/d2r

      integer  :: jd
      real(dp) :: caz, csz, cw, d, decl, ec, epsi, eqt, eyt, feqt, feqt1,       &
          feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi,     &
          pi, ra, raz, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr,    &
          yt, zpt, zr
!-------------------------------------------------------------------
! .. Intrinsic Functions ..
!-------------------------------------------------------------------
        intrinsic acos, atan, cos, min, sin, tan

!-------------------------------------------------------------------
! convert to radians
!-------------------------------------------------------------------
        rlt = lat*d2r
        rphi = long*d2r

        jd = ijd

        d = jd + gmt/24.0_dp
!-------------------------------------------------------------------
! calc geom mean longitude
!-------------------------------------------------------------------
        ml = 279.2801988_dp + d*(.9856473354_dp + 2.267E-13_dp*d)
        rml = ml*d2r

!-------------------------------------------------------------------
! calc equation of time in sec
! w = mean long of perigee
! e = eccentricity
! epsi = mean obliquity of ecliptic
!-------------------------------------------------------------------
        w = 282.4932328_dp + d*(4.70684E-5_dp + 3.39E-13_dp*d)
        wr = w*d2r
        ec   = 1.6720041E-2_dp - d*(1.1444E-9_dp + 9.4E-17_dp*d)
        epsi = 23.44266511_dp - d*(3.5626E-7_dp + 1.23E-15_dp*d)
        pepsi = epsi*d2r
        yt = (tan(pepsi/2.0_dp))**2
        cw = cos(wr)
        sw = sin(wr)
        ssw = sin(2.0_dp*wr)
        eyt = 2._dp*ec*yt
        feqt1 = -sin(rml)*cw*(eyt + 2._dp*ec)
        feqt2 = cos(rml)*sw*(2._dp*ec - eyt)
        feqt3 = sin(2._dp*rml)*(yt - (5._dp*ec**2/4._dp)*(cw**2 - sw**2))
        feqt4 = cos(2._dp*rml)*(5._dp*ec**2*ssw/4._dp)
        feqt5 = sin(3._dp*rml)*(eyt*cw)
        feqt6 = -cos(3._dp*rml)*(eyt*sw)
        feqt7 = -sin(4._dp*rml)*(.5_dp*yt**2)
        feqt = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7
        eqt = feqt*13751.0_dp

!-------------------------------------------------------------------
! convert eq of time from sec to deg
!-------------------------------------------------------------------
        reqt = eqt/240._dp
!-------------------------------------------------------------------
! calc right ascension in rads
!-------------------------------------------------------------------
        ra = ml - reqt
        rra = ra*d2r
!-------------------------------------------------------------------
! calc declination in rads, deg
!-------------------------------------------------------------------
        tab = 0.43360_dp*sin(rra)
        rdecl = atan(tab)
        decl = rdecl/d2r
!-------------------------------------------------------------------
! calc local hour angle
!-------------------------------------------------------------------
        lbgmt = 12.0_dp - eqt/3600._dp + long*24._dp/360._dp
        lzgmt = 15.0_dp*(gmt-lbgmt)
        zpt = lzgmt*d2r
        csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
        if(csz.gt.1)print *,'calczen,csz ',csz
        csz = min( 1._dp,csz )
!       zr = acos(csz)
!       zenith = zr/d2r
        zr = acos(csz)
        zenith = zr/d2r
!-------------------------------------------------------------------
! calc local solar azimuth
!-------------------------------------------------------------------
        caz = (sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr))
        if(caz < -0.999999_dp)then
          azimuth=180._dp
        elseif(caz > 0.999999_dp)then
          azimuth=0._dp
        else
          raz = acos(caz)
          azimuth = raz/d2r
        endif
!       caz = min(1.,(sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr)))
!       if(caz.lt.-1)print *,calczen ,caz
!       caz = max(-1.,caz)
!       raz = acos(caz)
!       azimuth = raz/d2r

        if( lzgmt > 0._dp ) azimuth = azimuth + (2._dp*(180._dp - azimuth))

      end subroutine calc_zenith

      subroutine sundis( julday, esrm2 )
!-----------------------------------------------------------------------------
! purpose:
! calculate earth-sun distance variation for a given date. based on
! fourier coefficients originally from: spencer, j.w., 1971, fourier
! series representation of the position of the sun, search, 2:172
!-----------------------------------------------------------------------------
! parameters:
! idate - integer, specification of the date, from yymmdd (i)
! esrm2 - real(dp), variation of the earth-sun distance (o)
! esrm2 = (average e/s dist)^2 / (e/s dist on day idate)^2
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: julday
      real(dp), intent(out) :: esrm2
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      real(dp), parameter :: pi = 3.1415926_dp
      real(dp) :: dayn, thet0
      real(dp) :: sinth, costh, sin2th, cos2th
!-----------------------------------------------------------------------------
! ... parse date to find day number (julian day)
!-----------------------------------------------------------------------------
      dayn = real(julday - 1,kind=dp) + .5_dp
!-----------------------------------------------------------------------------
! ... define angular day number and compute esrm2:
!-----------------------------------------------------------------------------
      thet0 = 2._dp*pi*dayn/365._dp
!-----------------------------------------------------------------------------
! ... calculate sin(2*thet0), cos(2*thet0) from
! addition theorems for trig functions for better
! performance; the computation of sin2th, cos2th
! is about 5-6 times faster than the evaluation
! of the intrinsic functions sin and cos
!-----------------------------------------------------------------------------
      sinth  = sin( thet0 )
      costh  = cos( thet0 )
      sin2th = 2._dp*sinth*costh
      cos2th = costh*costh - sinth*sinth
      esrm2  = 1.000110_dp + .034221_dp*costh + .001280_dp*sinth + .000719_dp*cos2th + .000077_dp*sin2th

      end subroutine sundis

      subroutine spheres( nz, z, zen, dsdh, nid )
!-----------------------------------------------------------------------------
! purpose:
! calculate slant path over vertical depth ds/dh in spherical geometry.
! calculation is based on: a.dahlback, and k.stamnes, a new spheric model
! for computing the radiation field available for photolysis and heating
! at twilight, planet.space sci., v39, n5, pp. 671-683, 1991 (appendix b)
!-----------------------------------------------------------------------------
! parameters:
! nz - integer, number of specified altitude levels in the working (i)
! grid
! z - real(dp), specified altitude working grid (km) (i)
! zen - real(dp), solar zenith angle (degrees) (i)
! dsdh - real(dp), slant path of direct beam through each layer crossed (o)
! when travelling from the top of the atmosphere to layer i;
! dsdh(i,j), i = 0..nz-1, j = 1..nz-1
! nid - integer, number of layers crossed by the direct beam when (o)
! travelling from the top of the atmosphere to layer i;
! nid(i), i = 0..nz-1
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in)  :: nz
      integer, intent(out) :: nid(0:nz-1)
      real(dp), intent(in) :: zen
      real(dp), intent(in) :: z(nz)
      real(dp), intent(out) :: dsdh(0:nz-1,nz-1)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      real(dp), parameter :: radius = 6.371e+3_dp         ! radius earth (km)
      real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp
      integer :: i, j, k
      integer :: id
      real(dp) :: re, ze(nz)
      real(dp) :: zd(0:nz-1)
      real(dp) :: zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm

      zenrad = zen*d2r
!-----------------------------------------------------------------------------
! ... include the elevation above sea level to the radius of the earth:
!-----------------------------------------------------------------------------
      re = radius + z(1)
!-----------------------------------------------------------------------------
! ... correspondingly z changed to the elevation above earth surface:
!-----------------------------------------------------------------------------
      ze(1:nz) = z(1:nz) - z(1)
!-----------------------------------------------------------------------------
! ... inverse coordinate of z
!-----------------------------------------------------------------------------
      zd(0) = ze(nz)
      do k = 1,nz-1
        zd(k) = ze(nz - k)
      end do
!-----------------------------------------------------------------------------
! ... initialize dsdh(i,j), nid(i)
!-----------------------------------------------------------------------------
      dsdh(0:nz-1,1:nz-1) = 0._dp
      nid(0:nz-1) = 0
!-----------------------------------------------------------------------------
! ... calculate ds/dh of every layer
!-----------------------------------------------------------------------------
      do i = 0,nz-1
        rpsinz = (re + zd(i)) * sin(zenrad)
        if( zen > 90._dp .and. rpsinz < re ) then
           nid(i) = -1
        else
!-----------------------------------------------------------------------------
! ... find index of layer in which the screening height lies
!-----------------------------------------------------------------------------
           id = i
           if( zen > 90._dp ) then
              do j = 1,nz-1
                 if( (rpsinz < ( zd(j-1) + re ) ) .and. (rpsinz >= ( zd(j) + re )) ) then
                    id = j
                 end if
              end do
           end if
           do j = 1,id
             if( j == id .and. id == i .and. zen > 90._dp ) then
                sm = -1._dp
             else
                sm = 1._dp
             end if
             rj = re + zd(j-1)
             rjp1 = re + zd(j)
             dhj = zd(j-1) - zd(j)
             ga = max( 0._dp,rj*rj - rpsinz*rpsinz )
             gb = max( 0._dp,rjp1*rjp1 - rpsinz*rpsinz )
             if( id > i .and. j == id ) then
                dsj = sqrt( ga )
             else
                dsj = sqrt( ga ) - sm*sqrt( gb )
             end if
             dsdh(i,j) = dsj / dhj
           end do
           nid(i) = id
        end if
      end do

      end subroutine spheres

      subroutine setz( nz, cz, tlev, c, ndx, adjcoe )

      use module_wave_data, only : tuv_jmax
!-----------------------------------------------------------------------------
! adjcoe - real(dp), coross section adjust coefficients
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      integer, intent(in) :: ndx
      real(dp), intent(in) :: cz(nz)
      real(dp), intent(in) :: tlev(nz)
      real(dp), intent(in) :: c(5,tuv_jmax)
      real(dp), intent(inout) :: adjcoe(nz,tuv_jmax)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: k, m
      real(dp) :: tt, adjin
      real(dp) :: c0, c1, c2
      real(dp) :: xz(nz)
!-----------------------------------------------------------------------------
! 1 o2 + hv -> o + o
! 2 o3 -> o2 + o(1d)
! 3 o3 -> o2 + o(3p)
! 4 no2 -> no + o(3p)
! 5 no3 -> no + o2
! 6 no3 -> no2 + o(3p)
! 7 n2o5 -> no3 + no + o(3p)
! 8 n2o5 -> no3 + no2
! 9 n2o + hv -> n2 + o(1d)
! 10 ho2 + hv -> oh + o
! 11 h2o2 -> 2 oh
! 12 hno2 -> oh + no
! 13 hno3 -> oh + no2
! 14 hno4 -> ho2 + no2
! 15 ch2o -> h + hco
! 16 ch2o -> h2 + co
! 17 ch3cho -> ch3 + hco
! 18 ch3cho -> ch4 + co
! 19 ch3cho -> ch3co + h
! 20 c2h5cho -> c2h5 + hco
! 21 chocho -> products
! 22 ch3cocho -> products
! 23 ch3coch3
! 24 ch3ooh -> ch3o + oh
! 25 ch3ono2 -> ch3o+no2
! 26 pan + hv -> products
! 27 mvk + hv -> products
! 28 macr + hv -> products
! 29 hyac + hv -> products
! 30 glyald + hv -> products
!-----------------------------------------------------------------------------
      xz(1:nz) = cz(1:nz)*1.e-18_dp
      do k = 1,tuv_jmax
         adjcoe(1:nz,k) = 1._dp
      end do
      tt = tlev(1)/281._dp
      do m = 1,tuv_jmax
         adjin = 1._dp
         if( m == 2 ) then
!----------------------------------------------------------------------
! ... temperature modification
! t0.9 (1.05) t0.95(1.025) t1.0(1.0) t1.15(1.02) t1.1(1.04)
!----------------------------------------------------------------------
            select case( ndx )
               case( 20 )
                  c0 = 4.52372_dp ; c1 = -5.94317_dp ; c2 = 2.63156_dp
               case( 40 )
                  c0 = 4.99378_dp ; c1 = -7.92752_dp ; c2 = 3.94715_dp
               case( 60 )
                  c0 = .969867_dp ; c1 = -.841035_dp ; c2 = .878835_dp
               case( 80 )
                  c0 = 1.07801_dp ; c1 = -2.39580_dp ; c2 = 2.32632_dp
            end select
            adjin = c0 + tt*(c1 + c2*tt)
         else if( m == 11 ) then
!----------------------------------------------------------------------
! ... temperature modification
! t0.9 (1.05) t0.95(1.025) t1.0(1.0) t1.15(1.02) t1.1(1.04)
!----------------------------------------------------------------------
            select case( ndx )
               case( 20 )
                  c0 = 2.43360_dp ; c1 = -3.61363_dp ; c2 = 2.19018_dp
               case( 40 )
                  c0 = 3.98265_dp ; c1 = -6.90516_dp ; c2 = 3.93602_dp
               case( 60 )
                  c0 = 3.49843_dp ; c1 = -5.98839_dp ; c2 = 3.50262_dp
               case( 80 )
                  c0 = 3.06312_dp ; c1 = -5.26281_dp ; c2 = 3.20980_dp
            end select
            adjin = c0 + tt*(c1 + c2*tt)
         end if
         call calcoe( nz, c(1,m), xz, tt, adjin, adjcoe(1,m) )
      end do

      end subroutine setz

      subroutine setozo( nz, z, tlay, dto3, &
                         to3, o3, airlev, o3top, o3toms )

      use module_wave_data, only : nw, wl, xso3, s226, s263, s298
!-----------------------------------------------------------------------------
! purpose:
! set up an altitude profile of ozone, and corresponding absorption
! optical depths. subroutine includes a shape-conserving scaling method
! that allows scaling of the entire profile to a given overhead ozone
! column amount.
!-----------------------------------------------------------------------------
! parameters:
! dobnew - real(dp), overhead ozone column amount (du) to which profile (i)
! should be scaled. if dobnew < 0, no scaling is done
! nz - integer, number of specified altitude levels in the working (i)
! grid
! z - real(dp), specified altitude working grid (km)
! nw - integer, number of specified intervals + 1 in working (i)
! wavelength grid
! wl - real(dp), vector of lower limits of wavelength intervals in (i)
! working wavelength grid
! xso3 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
! each specified wavelength (wmo value at 273)
! s226 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
! each specified wavelength (value from molina and molina at 226k)
! s263 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
! each specified wavelength (value from molina and molina at 263k)
! s298 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
! each specified wavelength (value from molina and molina at 298k)
! tlay - real(dp), temperature (k) at each specified altitude layer (i)
! dto3 - real(dp), optical depth due to ozone absorption at each (o)
! specified altitude at each specified wavelength
! to3 - real(dp), totol ozone column (o)
!-----------------------------------------------------------------------------

      implicit none

!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      real(dp), intent(in) :: o3top
      real(dp), intent(in) :: o3toms
      real(dp), intent(in) :: z(nz)
      real(dp), intent(in) :: tlay(nz-1)
      real(dp), intent(in) :: airlev(nz)
      real(dp), intent(out) :: dto3(nz-1,nw-1)
      real(dp), intent(out) :: to3(nz)
      real(dp), intent(in) :: o3(nz)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer  :: k, iw, nd, ilat
      real(dp) :: so3
      real(dp) :: scale
      real(dp) :: div1, div2
      real(dp) :: o3den(nz)
      real(dp) :: cz(nz)

!-----------------------------------------------------------------------------
! ... compute column increments
!-----------------------------------------------------------------------------
!     o3(1:nz) = o3(1:nz)*airlev(1:nz)
!     cz(1:nz-1) = (o3(2:nz)) * 1.e5 * (z(2:nz) - z(1:nz-1))

      o3den(1:nz) = o3(1:nz)*airlev(1:nz)
      cz(1:nz-1)  = 0.5_dp*(o3den(2:nz) + o3den(1:nz-1))*km2cm*(z(2:nz) - z(1:nz-1))

      to3(nz)  = o3top
      do k = nz-1,1,-1
         to3(k) = to3(k+1) + cz(k)
      end do

! Do not double count o3top
      cz(nz-1)    = cz(nz-1) + o3top
!-----------------------------------------------------------------------------
! ... scale o3 using toms data
!-----------------------------------------------------------------------------
      if( o3toms > 0.0_dp ) then
        scale = o3toms/(to3(1)/2.687e16_dp)
        cz(1:nz)  = cz(1:nz)*scale
        to3(1:nz) = to3(1:nz)*scale
      endif
!-----------------------------------------------------------------------------
! ... calculate ozone optical depth for each layer, with temperature
! correction. output, dto3(kz,kw)
!-----------------------------------------------------------------------------
      div1 = 1._dp / ( 263._dp - 226._dp)
      div2 = 1._dp / ( 298._dp - 263._dp)
      do iw = 1,nw-1
         so3 = xso3(iw)
         do k = 1,nz-1
            if( wl(iw) > 240.5_dp .and. wl(iw+1) < 350._dp ) then
               if( tlay(k) < 263._dp ) then
                  so3 = s226(iw) + (s263(iw) - s226(iw)) * div1 * (tlay(k) - 226._dp)
               else
                  so3 = s263(iw) + (s298(iw) - s263(iw)) * div2 * (tlay(k) - 263._dp)
               end if
            end if
            dto3(k,iw) = cz(k)*so3
         end do
      end do

      end subroutine setozo

      subroutine set_o2_xsect( nz, z, cz, vcol, scol, dto2, xso2 )

      use module_wave_data, only : nw, wl
!-----------------------------------------------------------------------------
! purpose:
! compute equivalent optical depths for o2 absorption, parameterized in
! the sr bands and the lyman-alpha line.
!-----------------------------------------------------------------------------
! parameters:
! nz - integer, number of specified altitude levels in the working (i)
! grid
! z - real(dp), specified altitude working grid (km)
! nw - integer, number of specified intervals + 1 in working
! wavelength grid
! wl - real(dp), vector of lower limits of wavelength intervals in
! working wavelength grid
! cz - real(dp), number of air molecules per cm^2 at each specified
! altitude layer
! zen - real(dp), solar zenith angle
! dto2 - real(dp), optical depth due to o2 absorption at each specified (o)
! vertical layer at each specified wavelength
! xso2 - real(dp), molecular absorption cross section in sr bands at (o)
! each specified altitude and wavelength. includes herzberg
! continuum.
!-----------------------------------------------------------------------------
! edit history:
! 02/98 included lyman-alpha parameterization
! 03/97 fix dto2 problem at top level (nz)
! 02/97 changed offset for grid-end interpolation to relative number
! (x * (1 +- deltax))
! 08/96 modified for early exit, no redundant read of data and smaller
! internal grid if possible; internal grid uses user grid points
! whenever possible
! 07/96 modified to work on internal grid and interpolate final values
! onto the user-defined grid
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      real(dp), intent(in) :: cz(nz)
      real(dp), intent(in) :: z(nz)
      real(dp), intent(in) :: vcol(nz)
      real(dp), intent(in) :: scol(nz)
      real(dp), intent(out) :: dto2(nz-1,nw-1)
      real(dp), intent(out) :: xso2(nw-1,nz)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: iw, iz, igast
      real(dp) :: secchi(nz)
!-----------------------------------------------------------------------------
! ... o2 optical depth and equivalent cross section on kockarts grid
!-----------------------------------------------------------------------------
      real(dp) :: dto2k(nz-1,ngast-1)
      real(dp) :: xso2k(nz,ngast-1)
!-----------------------------------------------------------------------------
! ... o2 optical depth and equivalent cross section in the lyman-alpha region
!-----------------------------------------------------------------------------
      real(dp) :: dto2la(nz-1,nla-1)
      real(dp) :: xso2la(nz,nla-1)
!-----------------------------------------------------------------------------
! ... temporary one-dimensional storage for optical depth and cross section values
! xxtmp - on internal grid
! xxuser - on user defined grid
!-----------------------------------------------------------------------------
      real(dp), dimension(2*nwint) :: dttmp, xstmp    !	CHANGED by Guohui
      real(dp), dimension(nwint) :: dtuser, xsuser    !	CHANGED by Guohui
      real(dp) :: o2col(nz)
      real(dp) :: x, y
      real(dp) :: delo2
!-----------------------------------------------------------------------------
! ... check, whether user grid is in the o2 absorption band at all...
! if not, set cross section and optical depth values to zero and return
!-----------------------------------------------------------------------------
      dto2(:nz-1,:nw-1) = 0._dp
      xso2(:nw-1,:nz) = 0._dp
      if( wl(1) > 243._dp ) then
         return
      end if
!-----------------------------------------------------------------------------
! ... sec xhi or chapman calculation
! for zen > 95 degrees, use zen = 95. (this is only to compute effective o2
! cross sections. still, better than setting dto2 = 0. as was done up to
! version 4.0) sm 1/2000
! in future could replace with mu2(iz) (but mu2 is also wavelength-depenedent)
! or imporved chapman function
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! ... slant o2 column
!-----------------------------------------------------------------------------
      o2col(1:nz) = 0.2095_dp * scol(1:nz)
!-----------------------------------------------------------------------------
! ... effective secant of solar zenith angle. use 2.0 if no direct sun.
! for nz, use value at nz-1
!-----------------------------------------------------------------------------
      secchi(1:nz-1) = scol(1:nz-1)/vcol(1:nz-1)
      where( secchi(1:nz-1) == 0._dp )
         secchi(1:nz-1) = 2._dp
      endwhere
      secchi(nz) = secchi(nz-1)
!-----------------------------------------------------------------------------
! ... if necessary:
! kockarts parameterization of the sr bands, output values of o2
! optical depth and o2 equivalent cross section are on his grid
!-----------------------------------------------------------------------------
      if( wl(1) < wlgast(ngast) .and. wl(nw) > wlgast(1) ) then
           call schu( nz, o2col, secchi, dto2k, xso2k )
      else
         dto2k(:,:) = 0._dp
         xso2k(:,:) = 0._dp
      end if
!-----------------------------------------------------------------------------
! ... lyman-alpha parameterization, output values of o2 opticaldepth
! and o2 effective (equivalent) cross section
!-----------------------------------------------------------------------------
      if( wl(1) <= wlla(nla) .and. wl(nw) >= wlla(1) ) then
         call lymana( nz, o2col, secchi, dto2la, xso2la )
      else
         dto2la(:,:) = 0._dp
         xso2la(:,:) = 0._dp
      end if
!-----------------------------------------------------------------------------
! ... loop through the altitude levels
!-----------------------------------------------------------------------------
      do iz = 1,nz
         igast = 0
!-----------------------------------------------------------------------------
! ... loop through the internal wavelength grid
!-----------------------------------------------------------------------------
         do iw = 1,nwint-1
!-----------------------------------------------------------------------------
! ... if outside kockarts grid and outside lyman-alpha, use the
! jpl/brasseur+solomon data, if inside
! kockarts grid, use the parameterized values from the call to schu,
! if inside lyman-alpha, use the paraemterized values from call to lymana
!-----------------------------------------------------------------------------
            if( wlint(iw+1) <= wlgast(1) .or. wlint(iw) >= wlgast(ngast) ) then
              if( wlint(iw+1) <= wlla(1) .or. wlint(iw) >= wlla(nla) ) then
                 xstmp(iw) = xso2int(iw)
              else
                 xstmp(iw) = xso2la(iz,1)
              end if
            else
               igast = igast + 1
               xstmp(iw) = xso2k(iz,igast)
            end if
!-----------------------------------------------------------------------------
! ... compute the area in each bin (for correct interpolation purposes only!)
!-----------------------------------------------------------------------------
            xstmp(iw) = xstmp(iw) * (wlint(iw+1) - wlint(iw))
         end do
!-----------------------------------------------------------------------------
! ... interpolate o2 cross section from the internal grid onto the user grid
!-----------------------------------------------------------------------------
         call inter3( nw, wl, xsuser, nwint, wlint, xstmp )
         xso2(1:nw-1,iz) = xsuser(1:nw-1)/(wl(2:nw) - wl(1:nw-1))
      end do
      do iz = 1,nz-1
         igast = 0
         delo2 = .2095_dp * cz(iz) ! vertical o2 column
!-----------------------------------------------------------------------------
! ... loop through the internal wavelength grid
!-----------------------------------------------------------------------------
         do iw = 1,nwint-1
!-----------------------------------------------------------------------------
! ... if outside kockarts grid and outside lyman-alpha, use the
! jpl/brasseur+solomon data, if inside
! kockarts grid, use the parameterized values from the call to schu,
! if inside lyman-alpha, use the paraemterized values from call to lymana
!-----------------------------------------------------------------------------
            if( wlint(iw+1) <= wlgast(1) .or. wlint(iw) >= wlgast(ngast) ) then
              if( wlint(iw+1) <= wlla(1) .or. wlint(iw) >= wlla(nla) ) then
                 dttmp(iw) = xso2int(iw) * delo2
              else
                 dttmp(iw) = dto2la(iz,1)
              end if
            else
               igast = igast + 1
               dttmp(iw) = dto2k(iz,igast)
            end if
!-----------------------------------------------------------------------------
! ... compute the area in each bin (for correct interpolation purposes only!)
!-----------------------------------------------------------------------------
            dttmp(iw) = dttmp(iw) * (wlint(iw+1) - wlint(iw))
         end do
!-----------------------------------------------------------------------------
! ... interpolate o2 optical depth from the internal grid onto the user grid
!-----------------------------------------------------------------------------
         call inter3( nw, wl, dtuser, nwint, wlint, dttmp )
         dto2(iz,1:nw-1) = dtuser(1:nw-1)/(wl(2:nw) - wl(1:nw-1))
      end do

      end subroutine set_o2_xsect

      subroutine setcld( nz, z, xlwc, dtcld, omcld, gcld )

      use module_wave_data, only : nw
!-----------------------------------------------------------------------------
!= PURPOSE:
!= Set up cloud optical depth, single albedo and g
!-----------------------------------------------------------------------------
!= PARAMETERS:
!= PARAMETERS:
!= NZ - INTEGER, number of specified altitude levels in the working (I)
!= grid
!= Z - real(dp), specified altitude working grid (km) (I)
!= XLWC Cloud water content g/M3 (I)
!=
!= dtcld - cloud optical depth
!= omcld - cloud droplet single albedo
!= gcld  - g
!-----------------------------------------------------------------------------
!
! VERTICAL DOMAIN is from bottom(1) to TOP (TOP=nz)
! CCM from top(1) to bottom(nz)
!-----------------------------------------------------------------------------
!****************************************************************************
      implicit none
!-----------------------------------------------------------------------------
! ... Dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      real(dp), intent(in) :: z(nz)
      real(dp), intent(in) :: xlwc(nz)
      real(dp), intent(out) :: dtcld(nz-1,nw-1)
      real(dp), intent(out) :: omcld(nz-1,nw-1)
      real(dp), intent(out) :: gcld(nz-1,nw-1)
!-----------------------------------------------------------------------------
! ... Local variables
!-----------------------------------------------------------------------------
      real(dp), parameter :: wden = 1.0_dp * 1.e6_dp ! g/m3 (1 m3 water = 1e6 g water)
      real(dp), parameter :: re = 10.0_dp * 1.e-6_dp ! assuming cloud drop radius = 10 um to M

      integer  :: i                                  ! vertical index for each blk
      real(dp) :: dz(nz-1)
!-----------------------------------------------------------------------------
! ... calculate vertical interveral dz
!-----------------------------------------------------------------------------
      do i=1,nz-1
          dz(i) = (z(i+1)-z(i)) * 1.e3_dp ! dz in Meters
      end do
!----------------------------------------------------
! ....calculate cloud optical depth T
! following Liao et al. JGR, 104, 23697, 1999
!----------------------------------------------------
      do i = 1,nz-1
        dtcld(i,:) = 1.5_dp * xlwc(i)*dz(i)/ (wden * re)
        dtcld(i,:) = max( dtcld(i,:), 0._dp )
        omcld(i,:) = .9999_dp
        gcld (i,:) = .85_dp
      end do

      end subroutine setcld

      subroutine schu_inti
!-----------------------------------------------------------------------------
! ... initialize the tables
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: i, iw, k, j1, jp1, j
      real(dp) :: col
      real(dp), dimension(6) :: a0, a1, b0, b1
      do iw = 1,ngast-1
         x_table(0,iw) = sum( a_schu(1:11:2,iw) )
         d_table(0,iw) = sum( b_schu(1:11:2,iw) )
         do k = 1,tdim
            col = 22._dp + t_del*real(k-1,kind=dp)
            o2_table(k) = col
            col = 10._dp**col
            a1(:) = a_schu(2:12:2,iw)*col
            b1(:) = b_schu(2:12:2,iw)*col
            where( a1(:) < 500._dp )
               a0(:) = exp( -a1(:) )
            elsewhere
               a0(:) = 0._dp
            endwhere
            where( b1(:) < 500._dp )
               b0(:) = exp( -b1(:) )
            elsewhere
               b0(:) = 0._dp
            endwhere
            x_table(k,iw) = dot_product( a_schu(1:11:2,iw),a0(:) )
            d_table(k,iw) = dot_product( b_schu(1:11:2,iw),b0(:) )
         end do
      end do

      end subroutine schu_inti

      subroutine schu( nz, o2col, secchi, dto2, xscho2 )
!-----------------------------------------------------------------------------
! purpose:
! calculate the equivalent absorption cross section of o2 in the sr bands.
! the algorithm is based on: g.kockarts, penetration of solar radiation
! in the schumann-runge bands of molecular oxygen: a robust approximation,
! annales geophysicae, v12, n12, pp. 1207ff, dec 1994. calculation is
! done on the wavelength grid used by kockarts (1994). final values do
! include effects from the herzberg continuum.
!-----------------------------------------------------------------------------
! parameters:
! nz - integer, number of specified altitude levels in the working (i)
! grid
! o2col - real(dp), slant overhead o2 column (molec/cc) at each specified (i)
! altitude
! dto2 - real(dp), optical depth due to o2 absorption at each specified (o)
! vertical layer at each specified wavelength
! xscho2 - real(dp), molecular absorption cross section in sr bands at (o)
! each specified wavelength. includes herzberg continuum
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      real(dp), intent(inout) :: dto2(nz-1,ngast-1)
      real(dp), intent(inout) :: xscho2(nz,ngast-1)
      real(dp), intent(in) :: o2col(nz)
      real(dp), intent(in) :: secchi(nz)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: i, iw, j, j1, jp1, k, ki, kp1
      integer :: index(nz)
      integer :: minind(1), maxind(1)
      real(dp) :: a0, a1, b0, b1
      real(dp), dimension(6) :: ac, bc
      real(dp), dimension(nz,6) :: aa, bb
      real(dp), dimension(nz) :: rjm, rjo2
      real(dp), dimension(nz) :: rjmi, rjo2i, lo2col, dels
!-----------------------------------------------------------------------------
! ... initialize r(m)
!-----------------------------------------------------------------------------
      rjm(1:nz) = 0._dp
      rjo2(1:nz) = 0._dp
!-----------------------------------------------------------------------------
! ... initialize the table interpolation
!-----------------------------------------------------------------------------
      where( o2col(:) /= 0 )
         lo2col(:) = log10( o2col(:) )
      endwhere
      do ki = 1,nz
         if( o2col(ki) /= 0._dp ) then
            if( lo2col(ki) <= o2_table(1) ) then
               dels(ki) = 0._dp
               index(ki) = 1
            else if( lo2col(ki) >= o2_table(tdim) ) then
               dels(ki) = 1._dp
               index(ki) = tdim-1
            else
               do k = 2,tdim
                  if( lo2col(ki) <= o2_table(k) ) then
                     dels(ki) = t_fac*(lo2col(ki) - o2_table(k-1))
                     index(ki) = k-1
                     exit
                  end if
               end do
            end if
         else
            index(ki) = 0
            dels(ki) = 0._dp
         end if
      end do
!-----------------------------------------------------------------------------
! ... calculate sum of exponentials (eqs 7 and 8 of kockarts 1994)
!-----------------------------------------------------------------------------
      do iw = 1,ngast-1
         do k = 1,nz
            ki = index(k)
            rjm(k) = x_table(ki,iw) + dels(k)*(x_table(ki+1,iw) - x_table(ki,iw))
            rjo2(k) = d_table(ki,iw) + dels(k)*(d_table(ki+1,iw) - d_table(ki,iw))
         end do
         do k = 1,nz-1
            if( rjm(k) > 1.e-100_dp ) then
               kp1 = k + 1
               if( rjm(kp1) > 0._dp ) then
                  dto2(k,iw) = log( rjm(kp1) ) / secchi(kp1) - log( rjm(k) ) * secchi(k)
               else
                  dto2(k,iw) = 1000._dp
               end if
            else
               dto2(k,iw) = 1000._dp
            end if
         end do
         do k = 1,nz
            if( rjm(k) > 1.e-100_dp ) then
               if( rjo2(k) > 1.e-100_dp ) then
                  xscho2(k,iw) = rjo2(k)/rjm(k)
               else
                  xscho2(k,iw) = 0._dp
               end if
            else
               xscho2(k,iw) = 0._dp
            end if
         end do
      end do

      end subroutine schu

      subroutine rtlink( nz,                  &
                         nw,                  &
                         zen,                 &
                         albedo,              & 
                         z,                   &
                         nid,                 &
                         dsdh,                & 
                         dtrl,                &
                         dto3,                &
                         dto2,                &
                         dtcld, omcld, gcld,  &
!                         dtcb1, omcb1, gcb1,  &
!                         dtcb2, omcb2, gcb2,  &
!                         dtoc1, omoc1, goc1,  &
!                         dtoc2, omoc2, goc2,  &
!                         dtant, omant, gant,  &
!                         dtso4, omso4, gso4,  &
!                         dtsal, omsal, gsal,  &
! rajesh: pass aerosol optical properties
                         dtaer, omaer, gaer,  &  
                         radfld )

      implicit none
!-----------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) :: nz
      integer, intent(in) :: nw
      real(dp), intent(in)    :: z(nz)
      real(dp), intent(in)    :: zen
      real(dp), intent(in)    :: albedo(nw-1)
      integer, intent(in)    :: nid(0:nz-1)
      real(dp), intent(in)    :: dsdh(0:nz-1,nz-1)
      real(dp), intent(in)    :: dtrl(nz-1,nw-1)
      real(dp), intent(in)    :: dto3(nz-1,nw-1)
      real(dp), intent(in)    :: dto2(nz-1,nw-1)
      real(dp), intent(in)    :: dtcld(nz-1,nw-1), omcld(nz-1,nw-1), gcld(nz-1,nw-1)
!      real(dp), intent(in)    :: dtcb1(nz-1,nw-1), omcb1(nz-1,nw-1), gcb1(nz-1,nw-1)
!      real(dp), intent(in)    :: dtcb2(nz-1,nw-1), omcb2(nz-1,nw-1), gcb2(nz-1,nw-1)
!      real(dp), intent(in)    :: dtoc1(nz-1,nw-1), omoc1(nz-1,nw-1), goc1(nz-1,nw-1)
!      real(dp), intent(in)    :: dtoc2(nz-1,nw-1), omoc2(nz-1,nw-1), goc2(nz-1,nw-1)
!      real(dp), intent(in)    :: dtant(nz-1,nw-1), omant(nz-1,nw-1), gant(nz-1,nw-1)
!      real(dp), intent(in)    :: dtso4(nz-1,nw-1), omso4(nz-1,nw-1), gso4(nz-1,nw-1)
!      real(dp), intent(in)    :: dtsal(nz-1,nw-1), omsal(nz-1,nw-1), gsal(nz-1,nw-1)
! rajesh: declare dust optical properties
      real(dp), intent(in)    :: dtaer(nz-1,nw-1), omaer(nz-1,nw-1), gaer(nz-1,nw-1) 
      real(dp), intent(out)   :: radfld(nz,nw-1)           
!-----------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------
      real(dp), parameter :: smallest = tiny( 1.0 )
      integer :: i, ii, iw
      real(dp) :: dtsct, dtabs
      real(dp) :: dscld, dacld
!      real(dp) :: dscb1, dacb1
!      real(dp) :: dscb2, dacb2
!      real(dp) :: dsoc1, daoc1
!      real(dp) :: dsoc2, daoc2
!      real(dp) :: dsant, daant
!      real(dp) :: dsso4, daso4
!      real(dp) :: dssal, dasal
! rajesh: declare scattering and absorbing optical depths
      real(dp) :: dsaer, daaer
      real(dp), dimension(nz-1,nw-1) :: dt, om, g
!-----------------------------------------------------------------------
! ... set any coefficients specific to rt scheme
!-----------------------------------------------------------------------
      do iw = 1,nw-1
         do i = 1,nz-1

            dscld = dtcld(i,iw)*omcld(i,iw)
            dacld = dtcld(i,iw)*(1._dp - omcld(i,iw))

!            dscb1 = dtcb1(i,iw)*omcb1(i,iw)
!            dacb1 = dtcb1(i,iw)*(1._dp - omcb1(i,iw))

!            dscb2 = dtcb2(i,iw)*omcb2(i,iw)
!            dacb2 = dtcb2(i,iw)*(1._dp - omcb2(i,iw))

!            dsoc1 = dtoc1(i,iw)*omoc1(i,iw)
!            daoc1 = dtoc1(i,iw)*(1._dp - omoc1(i,iw))

!            dsoc2 = dtoc2(i,iw)*omoc2(i,iw)
!            daoc2 = dtoc2(i,iw)*(1._dp - omoc2(i,iw))

!            dsant = dtant(i,iw)*omant(i,iw)
!            daant = dtant(i,iw)*(1._dp - omant(i,iw))

!            dsso4 = dtso4(i,iw)*omso4(i,iw)
!            daso4 = dtso4(i,iw)*(1._dp - omso4(i,iw))

!            dssal = dtsal(i,iw)*omsal(i,iw)
!            dasal = dtsal(i,iw)*(1._dp - omsal(i,iw))

! rajesh: determine scattering and absorption optical depths for dust
            dsaer = dtaer(i,iw)*omaer(i,iw)
            daaer = dtaer(i,iw)*(1._dp - omaer(i,iw))
            dtsct = dtrl(i,iw) + dscld + dsaer 
!                    dscb1 + dscb2 + dsoc1 + dsoc2 + dsant + dsso4 + dssal
            dtabs = dto3(i,iw) + dto2(i,iw) + dacld + daaer
!                    dacb1 + dacb2 + daoc1 + daoc2 + daant + daso4 + dasal

            dtabs = max( dtabs,smallest )
            dtsct = max( dtsct,smallest )
!-----------------------------------------------------------------------
! ... invert z-coordinate
!-----------------------------------------------------------------------
            ii = nz - i

            dt(ii,iw) = dtsct + dtabs
            if( dtsct /= smallest ) then
               om(ii,iw) = dtsct/(dtsct + dtabs)
               g(ii,iw) = ( gcld(i,iw)*dscld + &
!                            gcb1(i,iw)*dscb1 + &
!                            gcb2(i,iw)*dscb2 + &
!                            goc1(i,iw)*dsoc1 + &
!                            goc2(i,iw)*dsoc2 + &
!                            gant(i,iw)*dsant + &
!                            gso4(i,iw)*dsso4 + &
! rajesh: add aerosol contribution
                            gaer(i,iw)*dsaer ) / dtsct
               g(ii,iw) = max( smallest, g(ii,iw) )
               g(ii,iw) = min( 1.d0, g(ii,iw) )
            else
               om(ii,iw) = smallest
               g(ii,iw)  = smallest
            end if
!            print *, dt(ii,iw), om(ii,iw), g(ii,iw)
         end do
      end do

!-----------------------------------------------------------------------
! ... call rt routine
!-----------------------------------------------------------------------
      call ps2str( nz, nw, zen, albedo, dt, om, &
                   g, dsdh, nid, radfld )

      end subroutine rtlink

      subroutine tridec( syscnt, order, lower, main, upper )
!--------------------------------------------------------------------
! dimension of lower(syscnt,order), main(syscnt,order), upper(syscnt,order)
! arguments
!
! latest revision june 1995
!
! purpose trdi computes the solution of the tridiagonal
! linear system,
! b(1)*x(1)+c(1)*x(2) = y(1)
! a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1) = y(i)
! i=2,3,...,n-1
! a(n)*x(n-1)+b(n)*x(n) = y(n)
!
! usage call trdi (n, a, b, c, x )
!
! arguments
!
! on input n
! the number of unknowns. this subroutine
! requires that n be greater than 2.
!
! a
! the subdiagonal of the matrix is stored in
! locations a(2) through a(n).
!
! b
! the diagonal of the matrix is stored in
! locations b(1) through b(n).
!
! c
! the super-diagonal of the matrix is stored in
! locations c(1) through c(n-1).
!
! x
! the right-hand side of the equations is
! stored in y(1) through y(n).
!
! on output x
! an array which contains the solution to the
! system of equations.
!
! special conditions this subroutine executes satisfactorily
! if the input matrix is diagonally dominant
! and non-singular. the diagonal elements
! are used to pivot, and no tests are made to
! determine singularity. if a singular, or
! nearly singular, matrix is used as input,
! a divide by zero or floating point overflow
! may result.
!
! precision compiler dependent
!
! language fortran
!
! history originally written by nancy werner at ncar
! in october, 1973. modified by s. walters at
! ncar in june 1989 to functionally replace
! the cal routine tridla.
!
! portability fortran 90
!
! algorithm an lu-decomposition is obtained using the
! algorithm described in the reference below.
!
! to avoid unnecessary divisions, the alpha
! values used in the routine are the
! reciprocals of the alpha values described
! in the reference below.
!
! accuracy every component of the residual of the linear
! system (i.e. the difference between y and
! the matrix applied to x) should be less in
! magnitude than ten times the machine precision
! times the matrix order times the maximum
! absolute component of the solution vector
! times the largest absolute row sum of the
! input matrix.
!
! timing the timing is roughly proportional to the
! order n of the linear system.
!
! references analysis of numerical methods by
! e. isaacson and h. keller
! (john wiley and sons, 1966) pp. 55-58.
!--------------------------------------------------------------------
      implicit none
!--------------------------------------------------------------------
! ... dummy args
!--------------------------------------------------------------------
      integer, intent(in) :: syscnt, order
      real(dp), intent(in), dimension(syscnt,order) :: lower
      real(dp), intent(inout), dimension(syscnt,order) :: main, upper
!--------------------------------------------------------------------
! ... local variables
!--------------------------------------------------------------------
      integer :: i
!----------------------------------------------------------------------
! ... lu-decomposition
!----------------------------------------------------------------------
      main(:,1) = 1._dp / main(:,1)
      upper(:,1) = upper(:,1)*main(:,1)
      do i = 2,order-1
         main(:,i) = 1._dp / (main(:,i) - lower(:,i)*upper(:,i-1))
         upper(:,i) = upper(:,i)*main(:,i)
      end do

      end subroutine tridec

      subroutine trislv( syscnt, order, lower, main, upper, x )

      implicit none
!----------------------------------------------------------------------
! ... dummy args
!----------------------------------------------------------------------
      integer, intent(in) :: syscnt, order
      real(dp), intent(in), dimension(syscnt,order) :: lower, &
                                                       main, &
                                                       upper
      real(dp), intent(inout), dimension(syscnt,order) :: x
!----------------------------------------------------------------------
! ... local variables
!----------------------------------------------------------------------
      integer :: i, im1, j, n, nm1
      nm1 = order - 1
      n = order
!----------------------------------------------------------------------
! ... solve the system
!----------------------------------------------------------------------
      x(:,1) = x(:,1)*main(:,1)
      do i = 2,nm1
         x(:,i) = (x(:,i) - lower(:,i)*x(:,i-1))*main(:,i)
      end do
      x(:,n) = (x(:,n) - lower(:,n)*x(:,nm1))/(main(:,n) - lower(:,n)*upper(:,nm1))
      do i = nm1,1,-1
         x(:,i) = x(:,i) - upper(:,i)*x(:,i+1)
      end do

      end subroutine trislv

      subroutine ps2str( nz, nw, zen, rsfc, tauu, omu, &
                         gu, dsdh, nid, radfld )
!-----------------------------------------------------------------------------
! purpose:
! solve two-stream equations for multiple layers. the subroutine is based
! on equations from: toon et al., j.geophys.res., v94 (d13), nov 20, 1989.
! it contains 9 two-stream methods to choose from. a pseudo-spherical
! correction has also been added.
!-----------------------------------------------------------------------------
! parameters:
! nz - integer, number of specified altitude levels in the working (i)
! grid
! zen - real(dp), solar zenith angle (degrees) (i)
! rsfc - real(dp), surface albedo at current wavelength (i)
! tauu - real(dp), unscaled optical depth of each layer (i)
! omu - real(dp), unscaled single scattering albedo of each layer (i)
! gu - real(dp), unscaled asymmetry parameter of each layer (i)
! dsdh - real(dp), slant path of direct beam through each layer crossed (i)
! when travelling from the top of the atmosphere to layer i;
! dsdh(i,j), i = 0..nz-1, j = 1..nz-1
! nid - integer, number of layers crossed by the direct beam when (i)
! travelling from the top of the atmosphere to layer i;
! nid(i), i = 0..nz-1
! delta - logical, switch to use delta-scaling (i)
! .true. -> apply delta-scaling
! .false.-> do not apply delta-scaling
! fdr - real(dp), contribution of the direct component to the total (o)
! actinic flux at each altitude level
! fup - real(dp), contribution of the diffuse upwelling component to (o)
! the total actinic flux at each altitude level
! fdn - real(dp), contribution of the diffuse downwelling component to (o)
! the total actinic flux at each altitude level
! edr - real(dp), contribution of the direct component to the total (o)
! spectral irradiance at each altitude level
! eup - real(dp), contribution of the diffuse upwelling component to (o)
! the total spectral irradiance at each altitude level
! edn - real(dp), contribution of the diffuse downwelling component to (o)
! the total spectral irradiance at each altitude level
!-----------------------------------------------------------------------------

      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz, nw
      integer, intent(in) :: nid(0:nz-1)
      real(dp), intent(in) :: zen
      real(dp), intent(in) :: rsfc(nw-1)
      real(dp), dimension(nz-1,nw-1), intent(in) :: tauu, omu, gu
      real(dp), intent(in) :: dsdh(0:nz-1,nz-1)
      real(dp), intent(out) :: radfld(nz,nw-1)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! ... mu = cosine of solar zenith angle
! rsfc = surface albedo
! tauu = unscaled optical depth of each layer
! omu = unscaled single scattering albedo
! gu = unscaled asymmetry factor
! klev = max dimension of number of layers in atmosphere
! nlayer = number of layers in the atmosphere
! nlevel = nlayer + 1 = number of levels
!-----------------------------------------------------------------------------
      real(dp), parameter :: smallest = tiny( 1.0_dp )
      real(dp), parameter :: largest  = huge( 1.0_dp )
      real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp
      integer :: mrows, nzm1, nzm2
      real(dp), parameter :: eps = 1.e-3_dp
      real(dp), parameter :: pifs = 1._dp
      real(dp), parameter :: fdn0 = 0._dp
      integer :: row
      integer :: lev
      integer :: i, ip1, iw
      integer :: j, jl, ju
      real(dp) :: precis, wrk
      real(dp) :: tempg
      real(dp) :: mu, suma
      real(dp) :: g, om
      real(dp) :: gam1, gam2, gam3, gam4
      real(dp), dimension(nz-1) :: f, gi, omi
      real(dp), dimension(0:nz) :: tauc, mu2
      real(dp), dimension(nz-1) :: lam, taun, bgam
      real(dp), dimension(nz-1) :: cdn
      real(dp), dimension(0:nz,nw-1) :: tausla
      real(dp), dimension(nz-1,nw-1) :: cup, cuptn, cdntn
      real(dp), dimension(nz-1,nw-1) :: e1, e2, e3, e4
      real(dp), dimension(2*(nz-1)) :: a, b, d, e
      real(dp), dimension(nw-1,2*(nz-1)) :: sub, main, super, y
!-----------------------------------------------------------------------------
! ... for calculations of associated legendre polynomials for gama1,2,3,4
! in delta-function, modified quadrature, hemispheric constant,
! hybrid modified eddington-delta function metods, p633,table1.
! w.e.meador and w.r.weaver, gas,1980,v37,p.630
! w.j.wiscombe and g.w. grams, gas,1976,v33,p2440,
! uncomment the following two lines and the appropriate statements
! further down.
!-----------------------------------------------------------------------------
      real(dp) :: expon, expon0, expon1, divisr, temp, up, dn
      real(dp) :: ssfc
!-----------------------------------------------------------------------------
! ... initial conditions: pi*solar flux = 1; diffuse incidence = 0
!-----------------------------------------------------------------------------
      nzm1 = nz - 1
      nzm2 = nz - 2
      mrows = 2*nzm1
      precis = epsilon( precis )
      mu = cos( zen*d2r )
wave_loop : &
      do iw = 1,nw-1
!-----------------------------------------------------------------------------
! ... compute coefficients for each layer:
! gam1 - gam4 = 2-stream coefficients, different for different approximations
! expon0 = calculation of e when tau is zero
! expon1 = calculation of e when tau is taun
! cup and cdn = calculation when tau is zero
! cuptn and cdntn = calc. when tau is taun
! divisr = prevents division by zero
!-----------------------------------------------------------------------------
         tauc(0:nz) = 0._dp
         tausla(0:nz,iw) = 0._dp
         mu2(0:nz) = sqrt( smallest )
!-----------------------------------------------------------------------------
! ... delta-scaling. have to be done for delta-eddington approximation,
! delta discrete ordinate, practical improved flux method, delta function,
! and hybrid modified eddington-delta function methods approximations
!-----------------------------------------------------------------------------
         f(1:nzm1) = gu(:,iw)*gu(:,iw)
         gi(1:nzm1) = (gu(:,iw) - f(1:nzm1))/(1._dp - f(1:nzm1))
         omi(1:nzm1) = (1._dp - f(1:nzm1))*omu(1:nzm1,iw)/(1._dp - omu(1:nzm1,iw)*f(1:nzm1))
         taun(1:nzm1) = (1._dp - omu(1:nzm1,iw)*f(1:nzm1))*tauu(1:nzm1,iw)
!-----------------------------------------------------------------------------
! ... calculate slant optical depth at the top of the atmosphere when zen>90.
! in this case, higher altitude of the top layer is recommended which can
! be easily changed in gridz.f.
!-----------------------------------------------------------------------------
         if( zen > 90._dp ) then
            if( nid(0) < 0 ) then
               tausla(0,iw) = largest
            else
               ju = nid(0)
               tausla(0,iw) = 2._dp*dot_product( taun(1:ju),dsdh(0,1:ju) )
            end if
         end if
level_loop : &
         do i = 1,nzm1
            g = gi(i)
            om = omi(i)
            tauc(i) = tauc(i-1) + taun(i)
!-----------------------------------------------------------------------------
! ... stay away from 1 by precision. for g, also stay away from -1
!-----------------------------------------------------------------------------
            tempg = min( abs(g),1._dp - precis )
            g = sign( tempg,g )
            om = min( om,1._dp - precis )
!-----------------------------------------------------------------------------
! ... calculate slant optical depth
!-----------------------------------------------------------------------------
            if( nid(i) < 0 ) then
               tausla(i,iw) = largest
            else
               ju = min( nid(i),i )
               suma = dot_product( taun(1:ju),dsdh(i,1:ju) )
               jl = min( nid(i),i ) + 1
               tausla(i,iw) = suma + 2._dp*dot_product( taun(jl:nid(i)),dsdh(i,jl:nid(i)) )
               if( abs(tausla(i,iw)-tausla(i-1,iw)) < 1.0e-30_dp ) then
                 mu2(i) = sqrt( largest )
               else
                 mu2(i) = (tauc(i) - tauc(i-1))/(tausla(i,iw) - tausla(i-1,iw))
                 mu2(i) = sign( max( abs(mu2(i)),sqrt(smallest) ),mu2(i) )
               end if
            end if
!-----------------------------------------------------------------------------
! ... the following gamma equations are from pg 16,289, table 1
! eddington approximation(joseph et al., 1976, jas, 33, 2452):
!-----------------------------------------------------------------------------
            gam1 = .25_dp*(7._dp - om*(4._dp + 3._dp*g))
            gam2 = -.25_dp*(1._dp - om*(4._dp - 3._dp*g))
            gam3 = .25_dp*(2._dp - 3._dp*g*mu)
            gam4 = 1._dp - gam3
!-----------------------------------------------------------------------------
! ... lambda = pg 16,290 equation 21
! big gamma = pg 16,290 equation 22
!-----------------------------------------------------------------------------
            lam(i) = sqrt( gam1*gam1 - gam2*gam2 )
            bgam(i) = (gam1 - lam(i))/gam2
            wrk = lam(i)*taun(i)
            if( wrk < 500._dp ) then
               expon = exp( -wrk )
            else
               expon = 0._dp
            end if
!-----------------------------------------------------------------------------
! ... e1 - e4 = pg 16,292 equation 44
!-----------------------------------------------------------------------------
            e1(i,iw) = 1._dp + bgam(i)*expon
            e2(i,iw) = 1._dp - bgam(i)*expon
            e3(i,iw) = bgam(i) + expon
            e4(i,iw) = bgam(i) - expon
!-----------------------------------------------------------------------------
! ... the following sets up for the c equations 23, and 24
! found on page 16,290
! prevent division by zero (if lambda=1/mu, shift 1/mu^2 by eps = 1.e-3
! which is approx equiv to shifting mu by 0.5*eps* (mu)**3
!-----------------------------------------------------------------------------
            if( tausla(i-1,iw) < 500._dp ) then
               expon0 = exp( -tausla(i-1,iw) )
            else
               expon0 = 0._dp
            end if
            if( tausla(i,iw) < 500._dp ) then
               expon1 = exp( -tausla(i,iw) )
            else
               expon1 = 0._dp
            end if

            divisr = lam(i)*lam(i) - 1._dp/(mu2(i)*mu2(i))
            temp = max( eps,abs(divisr) )
            divisr = 1._dp/sign( temp,divisr )
            up = om*pifs*((gam1 - 1._dp/mu2(i))*gam3 + gam4*gam2)*divisr
            dn = om*pifs*((gam1 + 1._dp/mu2(i))*gam4 + gam2*gam3)*divisr
!-----------------------------------------------------------------------------
! ... cup and cdn are when tau is equal to zero
! cuptn and cdntn are when tau is equal to taun
!-----------------------------------------------------------------------------
            cup(i,iw) = up*expon0
            cdn(i) = dn*expon0
            cuptn(i,iw) = up*expon1
            cdntn(i,iw) = dn*expon1

         end do level_loop
!-----------------------------------------------------------------------------
! ... set up matrix
! ssfc = pg 16,292 equation 37 where pi fs is one (unity).
!-----------------------------------------------------------------------------
        if( tausla(nzm1,iw) < 500._dp ) then
           ssfc = rsfc(iw)*mu*exp( -tausla(nzm1,iw) )*pifs
        else
           ssfc = 0._dp
        end if
!-----------------------------------------------------------------------------
! ... the following are from pg 16,292 equations 39 - 43.
! set up first row of matrix:
!-----------------------------------------------------------------------------
        a(1) = 0._dp
        b(1) = e1(1,iw)
        d(1) = -e2(1,iw)
        e(1) = fdn0 - cdn(1)
!-----------------------------------------------------------------------------
! ... set up odd rows 3 thru (mrows - 1):
!-----------------------------------------------------------------------------
        a(3:mrows-1:2) = e2(1:nzm2,iw)*e3(1:nzm2,iw) - e4(1:nzm2,iw)*e1(1:nzm2,iw)
        b(3:mrows-1:2) = e1(1:nzm2,iw)*e1(2:nzm1,iw) - e3(1:nzm2,iw)*e3(2:nzm1,iw)
        d(3:mrows-1:2) = e3(1:nzm2,iw)*e4(2:nzm1,iw) - e1(1:nzm2,iw)*e2(2:nzm1,iw)
        e(3:mrows-1:2) = e3(1:nzm2,iw)*(cup(2:nzm1,iw) - cuptn(1:nzm2,iw)) + e1(1:nzm2,iw)*(cdntn(1:nzm2,iw) - cdn(2:nzm1))
!-----------------------------------------------------------------------------
! ... set up even rows 2 thru (mrows - 2):
!-----------------------------------------------------------------------------
        a(2:mrows-2:2) = e2(2:nzm1,iw)*e1(1:nzm2,iw) - e3(1:nzm2,iw)*e4(2:nzm1,iw)
        b(2:mrows-2:2) = e2(1:nzm2,iw)*e2(2:nzm1,iw) - e4(1:nzm2,iw)*e4(2:nzm1,iw)
        d(2:mrows-2:2) = e1(2:nzm1,iw)*e4(2:nzm1,iw) - e2(2:nzm1,iw)*e3(2:nzm1,iw)
        e(2:mrows-2:2) = (cup(2:nzm1,iw) - cuptn(1:nzm2,iw))*e2(2:nzm1,iw) - (cdn(2:nzm1) - cdntn(1:nzm2,iw))*e4(2:nzm1,iw)
!-----------------------------------------------------------------------------
! ... set up last row of matrix at mrows:
!-----------------------------------------------------------------------------
        a(mrows) = e1(nzm1,iw) - rsfc(iw)*e3(nzm1,iw)
        b(mrows) = e2(nzm1,iw) - rsfc(iw)*e4(nzm1,iw)
        d(mrows) = 0._dp
        e(mrows) = ssfc - cuptn(nzm1,iw) + rsfc(iw)*cdntn(nzm1,iw)
        sub(iw,1:mrows) = a(1:mrows)
        main(iw,1:mrows) = b(1:mrows)
        super(iw,1:mrows) = d(1:mrows)
        y(iw,1:mrows) = e(1:mrows)
      end do wave_loop
!-----------------------------------------------------------------------------
! ... solve the system
!-----------------------------------------------------------------------------
      call tridec( nw-1, mrows, sub, main, super )
      call trislv( nw-1, mrows, sub, main, super, y )
!-----------------------------------------------------------------------------
! ... unfold solution of matrix, compute output fluxes
!-----------------------------------------------------------------------------
      do iw = 1,nw-1
!-----------------------------------------------------------------------------
! ... the following equations are from pg 16,291 equations 31 & 32
!-----------------------------------------------------------------------------
         e(:mrows) = y(iw,:mrows)
         if( tausla(0,iw) < 500._dp ) then
            radfld(1,iw) = 2._dp*(fdn0 + e(1)*e3(1,iw) - e(2)*e4(1,iw) + cup(1,iw)) + exp( -tausla(0,iw) )
         else
            radfld(1,iw) = 2._dp*(fdn0 + e(1)*e3(1,iw) - e(2)*e4(1,iw) + cup(1,iw))
         end if
         where( tausla(1:nzm1,iw) < 500._dp )
            cdn(1:nzm1) = exp( -tausla(1:nzm1,iw) )
         elsewhere
            cdn(1:nzm1) = 0._dp
         endwhere
         radfld(2:nz,iw) = 2._dp*(e(1:mrows-1:2)*(e3(1:nzm1,iw) + e1(1:nzm1,iw)) &
                            + e(2:mrows:2)*(e4(1:nzm1,iw) + e2(1:nzm1,iw)) &
                            + cdntn(1:nzm1,iw) + cuptn(1:nzm1,iw)) + cdn(1:nzm1)
      end do

      end subroutine ps2str


      subroutine pchem( chem_opt, nz, tlev, airlev )

      use module_state_description, only : MOZART_KPP, MOZCART_KPP, &
                                           MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP
      use module_wave_data, only : r01, r04, r44, r08, r06, r10,  &
                                   r11, r14, r15, r17, r18,       &
                                   xs_mvk, xs_macr, xs_hyac, xs_glyald, nj

!-----------------------------------------------------------------------------
! purpose:
! load various "weighting functions" (products of cross section and
! quantum yield at each altitude and for wavelength). the altitude
! dependence is necessary to ensure the consideration of pressure and
! temperature dependence of the cross sections or quantum yields.
! the actual reading, evaluation and interpolation is done is separate
! subroutines for ease of management and manipulation. please refer to
! the inline documentation of the specific subroutines for detail information.
!-----------------------------------------------------------------------------
! parameters:
! nw - integer, number of specified intervals + 1 in working (i)
! wavelength grid
! wl - real(dp), vector of lower limits of wavelength intervals in (i)
! working wavelength grid
! nz - integer, number of altitude levels in working altitude grid (i)
! tlev - real(dp), temperature (k) at each specified altitude level (i)
! airlev - real(dp), air density (molec/cc) at each altitude level (i)
! j - integer, counter for number of weighting functions defined (io)
! sq - real(dp), cross section x quantum yield (cm^2) for each (o)
! photolysis reaction defined, at each defined wavelength and
! at each defined altitude level
! jlabel - character*40, string identifier for each photolysis reaction (o)
! defined
!-----------------------------------------------------------------------------

      implicit none

!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in)  :: chem_opt
      integer, intent(in)  :: nz
      real(dp), intent(in) :: tlev(nz)
      real(dp), intent(in) :: airlev(nz)

!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: j
      character(len=40) :: jlabel(nj)
!-----------------------------------------------------------------------------
! ... (1) o2 + hv -> o + o
! reserve first position. cross section parameterization in
! schumman-runge and lyman-alpha regions are zenith-angle dependent
!-----------------------------------------------------------------------------
      j = 1
      jlabel(j) = 'o2 + hv -> o + o '
!-----------------------------------------------------------------------------
! ... (2,3) o3 + hv -> (both channels) (tmp dep)
!-----------------------------------------------------------------------------
      call r01( nz, tlev, airlev, j, jlabel )
!-----------------------------------------------------------------------------
! ... (4) no2 + hv -> no + o(3p)
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'no2 + hv -> no + o(3p) '
!-----------------------------------------------------------------------------
! ... (5) no3 + hv -> no2 + o(3p)
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'no3 + hv -> no2 + o(3p) '
!-----------------------------------------------------------------------------
! ... (6) no3 + hv -> no + o2
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'no3 + hv -> no + o2 '
!-----------------------------------------------------------------------------
! ... (7,8) n2o5 + hv -> (both channels) (tmp dep)
!-----------------------------------------------------------------------------
      call r04( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (9) n2o + hv -> n2 + o(1d) (tmp dep)
!-----------------------------------------------------------------------------
      call r44( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (10) ho2 + hv -> oh + o
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ho2 + hv -> oh + o '
!-----------------------------------------------------------------------------
! ... (11) h2o2 + hv -> 2 oh
!-----------------------------------------------------------------------------
      call r08( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (12) hno2 + hv -> oh + no
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'hno2 + hv -> oh + no '
!-----------------------------------------------------------------------------
! ... (13) hno3 + hv -> oh + no2
!-----------------------------------------------------------------------------
      call r06( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (14) hno4 + hv -> ho2 + no2
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'hno4 + hv -> ho2 + no2 '
!-----------------------------------------------------------------------------
! ... (15,16) ch2o + hv -> (both channels)
!-----------------------------------------------------------------------------
      call r10( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (17,18,19) ch3cho + hv -> (all three channels)
!-----------------------------------------------------------------------------
      call r11( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (20) c2h5cho + hv -> c2h5 + hco
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'c2h5cho + hv -> c2h5 + hco '
!-----------------------------------------------------------------------------
! ... (21) chocho + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'chocho + hv -> products '
!-----------------------------------------------------------------------------
! ... (22) ch3cocho + hv -> products
!-----------------------------------------------------------------------------
      call r14( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (23) ch3coch3 + hv -> products
!-----------------------------------------------------------------------------
      call r15( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (24) ch3ooh + hv -> ch3o + oh
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ch3ooh + hv -> ch3o + oh '
!-----------------------------------------------------------------------------
! ... (25) ch3ono2 + hv -> ch3o + no2
!-----------------------------------------------------------------------------
      call r17( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (26) pan + hv -> products
!-----------------------------------------------------------------------------
      call r18( nz,tlev,airlev,j,jlabel )
!-----------------------------------------------------------------------------
! ... (27) mvk + hv -> products
!     (28) macr + hv -> products
!     (29) hyac + hv -> products
!     (30) glyald + hv -> products
!-----------------------------------------------------------------------------
is_mozart : &
      if( chem_opt == MOZART_KPP .or. chem_opt == MOZCART_KPP .or. &
          chem_opt == MOZART_MOSAIC_4BIN_KPP .or. &
          chem_opt == MOZART_MOSAIC_4BIN_AQ_KPP) then
        call xs_mvk( nz, tlev, airlev )
        call xs_macr( nz, tlev, airlev )
        call xs_hyac( nz, tlev, airlev )
        call xs_glyald( nz, tlev, airlev )
      else is_mozart
!-----------------------------------------------------------------------------
! ... (27) cloo + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'cloo + hv -> products '
!-----------------------------------------------------------------------------
! ... (28,29) clono2 + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'clono2 + hv -> cl + no3'
      j = j + 1
      jlabel(j) = 'clono2 + hv -> clo + no2'
!-----------------------------------------------------------------------------
! ... (30) ch3cl + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ch3cl + hv -> products '
!-----------------------------------------------------------------------------
! ... (31) ccl2o + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ccl2o + hv -> products'
!-----------------------------------------------------------------------------
! ... (32) ccl4 + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ccl4 + hv -> products '
!-----------------------------------------------------------------------------
! ... (33) cclfo + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'cclfo + hv -> products '
!-----------------------------------------------------------------------------
! ... (34) ccf2o + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ccf2o + hv -> products '
!-----------------------------------------------------------------------------
! ... (35) cf2clcfcl2 (cfc-113) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'cf2clcfcl2 (cfc-113) + hv -> products'
!-----------------------------------------------------------------------------
! ... (36) cf2clcf2cl (cfc-114) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'cf2clcf2cl (cfc-114) + hv -> products '
!-----------------------------------------------------------------------------
! ... (37) cf3cf2cl (cfc-115) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'cf3cf2cl (cfc-115) + hv -> products '
!-----------------------------------------------------------------------------
! ... (38) ccl3f (cfc-111) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ccl3f (cfc-111) + hv -> products'
!-----------------------------------------------------------------------------
! ... (39) ccl2f2 (cfc-112) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j) = 'ccl2f2 (cfc-112) + hv -> products'
!-----------------------------------------------------------------------------
! (40) ch3ccl3 + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='ch3ccl3 + hv -> products'
!-----------------------------------------------------------------------------
! (41) cf3chcl2 (hcfc-123) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf3chcl2 (hcfc-123) + hv -> products'
!-----------------------------------------------------------------------------
! (42) cf3chfcl (hcfc-124) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf3chfcl (hcfc-124) + hv -> products'
!-----------------------------------------------------------------------------
! (43) ch3cfcl2 (hcfc-141b) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='ch3cfcl2 (hcfc-141b) + hv -> products'
!-----------------------------------------------------------------------------
! (44) ch3cf2cl (hcfc-142b) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='ch3cf2cl (hcfc-142b) + hv -> products'
!-----------------------------------------------------------------------------
! (45) cf3cf2chcl2 (hcfc-225ca) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf3cf2chcl2 (hcfc-225ca) + hv -> products'
!-----------------------------------------------------------------------------
! (46) cf2clcf2chfcl (hcfc-225cb) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf2clcf2chfcl (hcfc-225cb) + hv -> products'
!-----------------------------------------------------------------------------
! (47) chclf2 (hcfc-22) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='chclf2 (hcfc-22) + hv -> products'
!-----------------------------------------------------------------------------
! (48) brono2 + hv -> br + o + no2
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='brono2 + hv -> br + o + no2'
!-----------------------------------------------------------------------------
! (49) brono2 + hv -> bro + no2
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='brono2 + hv -> bro + no2'
!-----------------------------------------------------------------------------
! (50) ch3br + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)=' ch3br + hv -> products'
!-----------------------------------------------------------------------------
! (51) chbr3 + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='chbr3 + hv -> products'
!-----------------------------------------------------------------------------
! (52) cf3br (halon-1301) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf3br (halon-1301) + hv -> products'
!-----------------------------------------------------------------------------
! (53) cf2brcf2br (halon-2402) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf2brcf2br (halon-2402) + hv -> products'
!-----------------------------------------------------------------------------
! (54) cf2br2 (halon-1202) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf2br2 (halon-1202) + hv -> products'
!-----------------------------------------------------------------------------
! (55) cf2brcl (halon-1211) + hv -> products
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cf2brcl (halon-1211) + hv -> products '
!-----------------------------------------------------------------------------
! (56) cl2 + hv -> cl + cl
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='cl2 + hv -> cl + cl '
!-----------------------------------------------------------------------------
! (57) hocl + hv -> oh + cl
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='hocl + hv -> oh + cl'
!-----------------------------------------------------------------------------
! (58) fmcl + hv -> cl + co + ho2
!-----------------------------------------------------------------------------
      j = j + 1
      jlabel(j)='fmcl + hv -> cl + co + ho2'

      end if is_mozart


      end subroutine pchem

      subroutine lymana( nz, o2col, secchi, dto2la, xso2la )
!-----------------------------------------------------------------------------
! purpose:
! calculate the effective absorption cross section of o2 in the lyman-alpha
! bands and an effective o2 optical depth at all altitudes. parameterized
! after: chabrillat, s., and g. kockarts, simple parameterization of the
! absorption of the solar lyman-alpha line, geophysical research letters,
! vol.24, no.21, pp 2659-2662, 1997.
!-----------------------------------------------------------------------------
! parameters:
! nz - integer, number of specified altitude levels in the working (i)
! grid
! o2col - real(dp), slant overhead o2 column (molec/cc) at each specified (i)
! altitude
! dto2la - real(dp), optical depth due to o2 absorption at each specified (o)
! vertical layer
! xso2la - real(dp), molecular absorption cross section in la bands (o)
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      real(dp), intent(in) :: o2col(nz)
      real(dp), intent(in) :: secchi(nz)
      real(dp), intent(out) :: dto2la(nz-1,nla-1)
      real(dp), intent(out) :: xso2la(nz,nla-1)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: i, k, kp1
      real(dp), dimension(nz) :: rm, ro2
      real(dp), save :: b(3), c(3), d(3), e(3)
      data b / 6.8431e-01_dp, 2.29841e-01_dp, 8.65412e-02_dp/, &
           c /8.22114e-21_dp, 1.77556e-20_dp, 8.22112e-21_dp/, &
           d / 6.0073e-21_dp, 4.28569e-21_dp, 1.28059e-20_dp/, &
           e /8.21666e-21_dp, 1.63296e-20_dp, 4.85121e-17_dp/
!-----------------------------------------------------------------------------
! ... calculate reduction factors at every altitude
!-----------------------------------------------------------------------------
      rm(:) = 0._dp
      ro2(:) = 0._dp
      do k = 1,nz
        do i = 1,3
          rm(k) = rm(k) + b(i) * exp( -c(i)*o2col(k) )
          ro2(k) = ro2(k) + d(i) * exp( -e(i)*o2col(k) )
        end do
      end do
!-----------------------------------------------------------------------------
! ... calculate effective o2 optical depths and effective o2 cross sections
!-----------------------------------------------------------------------------
      do k = 1,nz-1
         if( rm(k) > 1.e-100_dp ) then
            kp1 = k + 1
            if( rm(kp1) > 0._dp ) then
               dto2la(k,1) = log( rm(kp1) )/secchi(kp1) - log( rm(k) )/secchi(k)
            else
               dto2la(k,1) = 1000._dp
            end if
         else
            dto2la(k,1) = 1000._dp
         end if
      end do
      do k = 1,nz
         if( rm(k) > 1.e-100_dp ) then
            if( ro2(k) > 1.e-100_dp ) then
               xso2la(k,1) = ro2(k)/rm(k)
            else
               xso2la(k,1) = 0._dp
            end if
         else
            xso2la(k,1) = 0._dp
         end if
      end do

      end subroutine lymana


      subroutine inter2( ng, xg, yg, n, x, y, ierr )
!-----------------------------------------------------------------------------
! purpose:
! map input data given on single, discrete points onto a set of target
! bins.
! the original input data are given on single, discrete points of an
! arbitrary grid and are being linearly interpolated onto a specified set
! of target bins. in general, this is the case for most of the weighting
! functions (action spectra, molecular cross section, and quantum yield
! data), which have to be matched onto the specified wavelength intervals.
! the average value in each target bin is found by averaging the trapezoi-
! dal area underneath the input data curve (constructed by linearly connec-
! ting the discrete input values).
! some caution should be used near the endpoints of the grids. if the
! input data set does not span the range of the target grid, an error
! message is printed and the execution is stopped, as extrapolation of the
! data is not permitted.
! if the input data does not encompass the target grid, use addpnt to
! expand the input array.
!-----------------------------------------------------------------------------
! parameters:
! ng - integer, number of bins + 1 in the target grid (i)
! xg - real(dp), target grid (e.g., wavelength grid); bin i is defined (i)
! as [xg(i),xg(i+1)] (i = 1..ng-1)
! yg - real(dp), y-data re-gridded onto xg, yg(i) specifies the value for (o)
! bin i (i = 1..ng-1)
! n - integer, number of points in input grid (i)
! x - real(dp), grid on which input data are defined (i)
! y - real(dp), input y-data (i)
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: ng, n
      integer, intent(out) :: ierr
      real(dp), intent(in) :: x(n)
      real(dp), intent(in) :: y(n)
      real(dp), intent(in) :: xg(ng)
      real(dp), intent(out) :: yg(ng)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: ngintv
      integer :: i, k, jstart
      real(dp) :: area, xgl, xgu
      real(dp) :: darea, slope
      real(dp) :: a1, a2, b1, b2
      ierr = 0
!-----------------------------------------------------------------------------
! ... test for correct ordering of data, by increasing value of x
!-----------------------------------------------------------------------------
      do i = 2,n
         if( x(i) <= x(i-1) ) then
            ierr = 1
            write(*,*) 'inter2: x coord not monotonically increasing'
            return
         end if
      end do
      do i = 2,ng
        if( xg(i) <= xg(i-1) ) then
           ierr = 2
           write(*,*) 'inter2: xg coord not monotonically increasing'
           return
        end if
      end do
!-----------------------------------------------------------------------------
! ... check for xg-values outside the x-range
!-----------------------------------------------------------------------------
      if( x(1) > xg(1) .or. x(n) < xg(ng) ) then
          write(*,*) 'inter2: data does not span grid'
          write(*,*) '        use addpnt to expand data and re-run'
! call endrun
      end if
!-----------------------------------------------------------------------------
! ... find the integral of each grid interval and use this to
! calculate the average y value for the interval
! xgl and xgu are the lower and upper limits of the grid interval
!-----------------------------------------------------------------------------
      jstart = 1
      ngintv = ng - 1
      do i = 1,ngintv
!-----------------------------------------------------------------------------
! ... initalize:
!-----------------------------------------------------------------------------
         area = 0._dp
         xgl = xg(i)
         xgu = xg(i+1)
!-----------------------------------------------------------------------------
! ... discard data before the first grid interval and after the
! last grid interval
! for internal grid intervals, start calculating area by interpolating
! between the last point which lies in the previous interval and the
! first point inside the current interval
!-----------------------------------------------------------------------------
         k = jstart
         if( k <= n-1 ) then
!-----------------------------------------------------------------------------
! ... if both points are before the first grid, go to the next point
!-----------------------------------------------------------------------------
            do
               if( x(k+1) <= xgl ) then
                  jstart = k - 1
                  k = k+1
                  if( k <= n-1 ) then
                     cycle
                  else
                     exit
                  end if
               else
                  exit
               end if
            end do
!-----------------------------------------------------------------------------
! ... if the last point is beyond the end of the grid,
! complete and go to the next grid
!-----------------------------------------------------------------------------
            do
               if( k <= n-1 .and. x(k) < xgu ) then
                  jstart = k-1
!-----------------------------------------------------------------------------
! ... compute x-coordinates of increment
!-----------------------------------------------------------------------------
                  a1 = max( x(k),xgl )
                  a2 = min( x(k+1),xgu )
!-----------------------------------------------------------------------------
! ... if points coincide, contribution is zero
!-----------------------------------------------------------------------------
                  if( x(k+1) == x(k) ) then
                     darea = 0._dp
                  else
                     slope = (y(k+1) - y(k))/(x(k+1) - x(k))
                     b1 = y(k) + slope*(a1 - x(k))
                     b2 = y(k) + slope*(a2 - x(k))
                     darea = .5*(a2 - a1)*(b2 + b1)
                  end if
!-----------------------------------------------------------------------------
! ... find the area under the trapezoid from a1 to a2
!-----------------------------------------------------------------------------
                  area = area + darea
                  k = k+1
                  cycle
               else
                  exit
               end if
            end do
         end if
!-----------------------------------------------------------------------------
! ... calculate the average y after summing the areas in the interval
!-----------------------------------------------------------------------------
         yg(i) = area/(xgu - xgl)
      end do

      end subroutine inter2

      subroutine inter_inti( ng, xg, n, x )
!-----------------------------------------------------------------------------
! ... initialization
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: ng, n
      real(dp), intent(in) :: xg(ng)
      real(dp), intent(in) :: x(n)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: i, ii, iil, astat
      integer :: ndim(1)
      allocate( xi(ng), xcnt(ng-1), stat=astat )
      if( astat /= 0 ) then
         write(*,*) 'inter_inti: failed to allocate wrk arrays; error = ',astat
         stop
      else
         xi(:) = 0
         xcnt(:) = 0
      end if
      iil = 1
      do i = 1,ng
         do ii = iil,n-1
            if( xg(i) < x(ii) ) then
               xi(i) = ii - 1
               iil = ii
               exit
            end if
         end do
      end do
      nintervals = count( xi(:) /= 0 )
      if( nintervals == 0 ) then
         write(*,*) 'inter_inti: wavelength grids do not overlap'
         stop
      else
         nintervals = nintervals - 1
      end if
      xcnt(1:nintervals) = xi(2:nintervals+1) - xi(1:nintervals) + 1
      ndim(:) = maxval( xcnt(1:nintervals) )
      allocate( xfrac(ndim(1),nintervals),stat=astat )
      if( astat /= 0 ) then
         write(*,*) 'inter_inti: failed to allocate wrk array; error = ',astat
         stop
      else
         xfrac(:,:) = 1._dp
      end if
      do i = 1,nintervals
        iil = xi(i)
        xfrac(1,i) = (min( x(iil+1),xg(i+1) ) - xg(i))/(x(iil+1) - x(iil))
        if( xcnt(i) > 1 ) then
           iil = xi(i) + xcnt(i) - 1
           xfrac(xcnt(i),i) = (xg(i+1) - x(iil))/(x(iil+1) - x(iil))
        end if
      end do
!     write(*,*) 'inter_inti: diagnostics; ng,n,nintervals = ',ng,n,nintervals
!     write(*,'('' xi'')')
!     write(*,'(10i4)') xi(1:nintervals+1)
!     write(*,'('' xcnt'')')
!     write(*,'(10i4)') xcnt(1:nintervals)
!     write(*,'('' xfrac'')')
!     do i = 1,nintervals
!        write(*,'(1p,5e21.13)') xfrac(1:xcnt(i),i)
!     end do

      end subroutine inter_inti

      subroutine inter3( ng, xg, yg, n, x, y )
!-----------------------------------------------------------------------------
! purpose:
! map input data given on a set of bins onto a different set of target
! bins.
! the input data are given on a set of bins (representing the integral
! of the input quantity over the range of each bin) and are being matched
! onto another set of bins (target grid). a typical example would be an
! input data set spcifying the extra-terrestrial flux on wavelength inter-
! vals, that has to be matched onto the working wavelength grid.
! the resulting area in a given bin of the target grid is calculated by
! simply adding all fractional areas of the input data that cover that
! particular target bin.
! some caution should be used near the endpoints of the grids. if the
! input data do not span the full range of the target grid, the area in
! the "missing" bins will be assumed to be zero. if the input data extend
! beyond the upper limit of the target grid, the user has the option to
! integrate the "overhang" data and fold the remaining area back into the
! last target bin. using this option is recommended when re-gridding
! vertical profiles that directly affect the total optical depth of the
! model atmosphere.
!-----------------------------------------------------------------------------
! parameters:
! ng - integer, number of bins + 1 in the target grid (i)
! xg - real(dp), target grid (e.g. working wavelength grid); bin i (i)
! is defined as [xg(i),xg(i+1)] (i = 1..ng-1)
! yg - real(dp), y-data re-gridded onto xg; yg(i) specifies the (o)
! y-value for bin i (i = 1..ng-1)
! n - integer, number of bins + 1 in the input grid (i)
! x - real(dp), input grid (e.g. data wavelength grid); bin i is (i)
! defined as [x(i),x(i+1)] (i = 1..n-1)
! y - real(dp), input y-data on grid x; y(i) specifies the (i)
! y-value for bin i (i = 1..n-1)
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: n, ng
      real(dp), intent(in) :: xg(ng)
      real(dp), intent(in) :: x(n)
      real(dp), intent(in) :: y(n)
      real(dp), intent(out) :: yg(ng)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: i, ii, iil
!-----------------------------------------------------------------------------
! ... do interpolation
!-----------------------------------------------------------------------------
      yg(:) = 0.
      do i = 1,nintervals
         iil = xi(i)
         ii = xcnt(i)
         if( ii == 1 ) then
            yg(i) = xfrac(1,i)*y(iil)
         else
            yg(i) = dot_product( xfrac(1:ii,i),y(iil:iil+ii-1) )
         end if
      end do

      end subroutine inter3

      subroutine calcoe( nz, c, xz, tt, adjin, adjcoe )
!-----------------------------------------------------------------------------
! parameters:
! adjcoe - real(dp), coross section adjust coefficients (in and out)
! c(5,28)-polynomal coef
! tt -nomarlized temperature
!-----------------------------------------------------------------------------*
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      real(dp), intent(in) :: adjin
      real(dp), intent(in) :: tt
      real(dp), intent(in) :: c(5)
      real(dp), intent(in) :: xz(nz)
      real(dp), intent(inout) :: adjcoe(nz)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: k
      real(dp) :: x
      do k = 1,nz
         x = xz(k)
         adjcoe(k) = adjin * (1._dp + .01_dp*(c(1) + x*(c(2) + x*(c(3) + x*(c(4) + x*c(5))))))
      end do

      end subroutine calcoe


      subroutine airmas( nz, z, zen, dsdh, nid, cz, &
                         vcol, scol )
!-----------------------------------------------------------------------------
! purpose:
! calculate vertical and slant air columns, in spherical geometry, as a
! function of altitude.
!-----------------------------------------------------------------------------
! parameters:
! nz - integer, number of specified altitude levels in the working (i)
! grid
! z - real(dp), specified altitude working grid (km) (i)
! zen - real(dp), solar zenith angle (degrees) (i)
! dsdh - real(dp), slant path of direct beam through each layer crossed (o)
! when travelling from the top of the atmosphere to layer i;
! dsdh(i,j), i = 0..nz-1, j = 1..nz-1
! nid - integer, number of layers crossed by the direct beam when (o)
! travelling from the top of the atmosphere to layer i;
! nid(i), i = 0..nz-1
! vcol - real(dp), output, vertical air column, molec cm-2, above level iz
! scol - real(dp), output, slant air column in direction of sun, above iz
! also in molec cm-2
!-----------------------------------------------------------------------------
      implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in) :: nz
      integer, intent(in) :: nid(0:nz-1)
      real(dp), intent(in) :: z(nz)
      real(dp), intent(in) :: zen
      real(dp), intent(in) :: dsdh(0:nz-1,nz-1)
      real(dp), intent(in) :: cz(nz)
      real(dp), intent(out) :: vcol(nz)
      real(dp), intent(out) :: scol(nz)
!-----------------------------------------------------------------------------
! ... parameters
!-----------------------------------------------------------------------------
      real(dp), parameter :: largest = huge(1.0_dp)
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
      integer :: id, j
      real(dp) :: sum, ssum, vsum, ratio
!-----------------------------------------------------------------------------
! ... calculate vertical and slant column from each level: work downward
!-----------------------------------------------------------------------------
      vsum = 0._dp
      ssum = 0._dp
      do id = 1,nz-1
         vsum = vsum + cz(nz-id)
         vcol(nz-id) = vsum
         sum = 0.
         if( nid(id) < 0 ) then
            sum = largest
         else
!-----------------------------------------------------------------------------
! ... single pass layers:
!-----------------------------------------------------------------------------
            do j = 1,min( nid(id),id )
               sum = sum + cz(nz-j)*dsdh(id,j)
            end do
!-----------------------------------------------------------------------------
! ... double pass layers:
!-----------------------------------------------------------------------------
            do j = min( nid(id),id )+1,nid(id)
               sum = sum + 2._dp*cz(nz-j)*dsdh(id,j)
            end do
         end if
         scol(nz - id) = sum
      end do
!-----------------------------------------------------------------------------
! ... special section to set scol(nz)
!-----------------------------------------------------------------------------
      if( scol(nz-2) /= 0._dp ) then
         ratio = scol(nz-1)/scol(nz-2)
         scol(nz) = ratio * scol(nz-1)
      else
         scol(nz) = 0._dp
      end if

      end subroutine airmas

      END MODULE module_ftuv_subs