chem/module_ftuv_subs.F

   1
   2       MODULE module_ftuv_subs
   3
   4       use module_wave_data, only : nw
   5
   6       implicit none
   7
   8       integer, parameter :: dp = selected_real_kind(14,300)
   9 !-----------------------------------------------------------------------------
  10 !     data for inter2, inter3 and interi
  11 !-----------------------------------------------------------------------------
  12       integer, private, save :: nintervals
  13       integer, private, allocatable, save :: xi(:), xcnt(:)
  14       real(dp), private, allocatable, save :: xfrac(:,:)
  15
  16       real(dp), private, parameter :: km2cm = 1.e5_dp
  17
  18 !-----------------------------------------------------------------------------
  19 !     data for lyman
  20 !-----------------------------------------------------------------------------
  21       integer, private, parameter :: nla = 2
  22       integer, private, parameter :: ngast = 17
  23
  24 !-----------------------------------------------------------------------------
  25 !     data for schu_inti, schu
  26 !-----------------------------------------------------------------------------
  27       integer, private, parameter  :: tdim = 501
  28       !BSINGH(PNNL)-Lahey compiler doesn't allow use of 'real' in variable declaration
  29       real(dp), private, parameter  :: tdim_real = 501.0
  30       !real(dp), private, parameter :: t_del = 5._dp/real(tdim-1,kind=dp)!BSINGH(PNNL)-commented out
  31       !real(dp), private, parameter :: t_fac = real(tdim-1,kind=dp)/5._dp!BSINGH(PNNL)-commented out
  32       real(dp), private, parameter :: t_del = 5._dp/(tdim_real-1._dp)!BSINGH(PNNL)-redefined t_del and t_fac
  33       real(dp), private, parameter :: t_fac = (tdim_real-1._dp)/5._dp
  34
  35       integer, private :: ii, jj
  36       real(dp), private, save :: d_table(0:tdim,ngast-1)
  37       real(dp), private, save :: x_table(0:tdim,ngast-1)
  38       real(dp), private, save :: o2_table(tdim)
  39       real(dp), private, dimension(12,ngast-1), save :: a_schu, b_schu
  40 !-----------------------------------------------------------------------------
  41 ! a_schu(16,12) coefficients for rj(m) (table 1 in kockarts 1994)
  42 ! b_schu(16,12) rj(o2)(table 2 in kockarts 1994)
  43 ! rjm attenuation coefficients rj(m)
  44 ! rjo2 rj(o2)
  45 !-----------------------------------------------------------------------------
  46       data ((a_schu(jj,ii),jj=1,12),ii=1,ngast-1) /                             &
  47 !a 57000-56500.5 cm-1
  48        1.13402e-01,1.00088e-20,3.48747e-01,2.76282e-20,3.47322e-01,1.01267e-19, &
  49        1.67351e-01,5.63588e-19,2.31433e-02,1.68267e-18,0.00000e+00,0.00000e+00, &
  50 !a 56500-56000.5 cm-1
  51        2.55268e-03,1.64489e-21,1.85483e-01,2.03591e-21,2.60603e-01,4.62276e-21, &
  52        2.50337e-01,1.45106e-20,1.92340e-01,7.57381e-20,1.06363e-01,7.89634e-19, &
  53 !a 56000-55500.5 cm-1
  54        4.21594e-03,8.46639e-22,8.91886e-02,1.12935e-21,2.21334e-01,1.67868e-21, &
  55        2.84446e-01,3.94782e-21,2.33442e-01,1.91554e-20,1.63433e-01,2.25346e-19, &
  56 !a 55500-55000.5 cm-1
  57        3.93529e-03,6.79660e-22,4.46906e-02,9.00358e-22,1.33060e-01,1.55952e-21, &
  58        3.25506e-01,3.43763e-21,2.79405e-01,1.62086e-20,2.10316e-01,1.53883e-19, &
  59 !a 55000-54500.5 cm-1
  60        2.60939e-03,2.33791e-22,2.08101e-02,3.21734e-22,1.67186e-01,5.77191e-22, &
  61        2.80694e-01,1.33362e-21,3.26867e-01,6.10533e-21,1.96539e-01,7.83142e-20, &
  62 !a 54500-54000.5 cm-1
  63        9.33711e-03,1.32897e-22,3.63980e-02,1.78786e-22,1.46182e-01,3.38285e-22, &
  64        3.81762e-01,8.93773e-22,2.58549e-01,4.28115e-21,1.64773e-01,4.67537e-20, &
  65 !a 54000-53500.5 cm-1
  66        9.51799e-03,1.00252e-22,3.26320e-02,1.33766e-22,1.45962e-01,2.64831e-22, &
  67        4.49823e-01,6.42879e-22,2.14207e-01,3.19594e-21,1.45616e-01,2.77182e-20, &
  68 !a 53500-53000.5 cm-1
  69        7.87331e-03,3.38291e-23,6.91451e-02,4.77708e-23,1.29786e-01,8.30805e-23, &
  70        3.05103e-01,2.36167e-22,3.35007e-01,8.59109e-22,1.49766e-01,9.63516e-21, &
  71 !a 53000-52500.5 cm-1
  72        6.92175e-02,1.56323e-23,1.44403e-01,3.03795e-23,2.94489e-01,1.13219e-22, &
  73        3.34773e-01,3.48121e-22,9.73632e-02,2.10693e-21,5.94308e-02,1.26195e-20, &
  74 !a 52500-52000.5 cm-1
  75        1.47873e-01,8.62033e-24,3.15881e-01,3.51859e-23,4.08077e-01,1.90524e-22, &
  76        8.08029e-02,9.93062e-22,3.90399e-02,6.38738e-21,8.13330e-03,9.93644e-22, &
  77 !a 52000-51500.5 cm-1
  78        1.50269e-01,1.02621e-23,2.39823e-01,3.48120e-23,3.56408e-01,1.69494e-22, &
  79        1.61277e-01,6.59294e-22,8.89713e-02,2.94571e-21,3.25063e-03,1.25548e-20, &
  80 !a 51500-51000.5 cm-1
  81        2.55746e-01,8.49877e-24,2.94733e-01,2.06878e-23,2.86382e-01,9.30992e-23, &
  82        1.21011e-01,3.66239e-22,4.21105e-02,1.75700e-21,0.00000e+00,0.00000e+00, &
  83 !a 51000-50500.5 cm-1
  84        5.40111e-01,7.36085e-24,2.93263e-01,2.46742e-23,1.63417e-01,1.37832e-22, &
  85        3.23781e-03,2.15052e-21,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
  86 !a 50500-50000.5 cm-1
  87        8.18514e-01,7.17937e-24,1.82262e-01,4.17496e-23,0.00000e+00,0.00000e+00, &
  88        0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
  89 !a 50000-49500.5 cm-1
  90        8.73680e-01,7.13444e-24,1.25583e-01,2.77819e-23,0.00000e+00,0.00000e+00, &
  91        0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
  92 !a 49500-49000.5 cm-1
  93        3.32476e-04,7.00362e-24,9.89000e-01,6.99600e-24,0.00000e+00,0.00000e+00, &
  94        0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 /
  95       data ((b_schu(jj,ii),jj=1,12),ii=1,ngast-1) /                             &
  96 ! 57000-56500.5 cm-1
  97        1.07382e-21,9.95029e-21,7.19430e-21,2.48960e-20,2.53735e-20,7.54467e-20, &
  98        4.48987e-20,2.79981e-19,9.72535e-20,9.29745e-19,2.30892e-20,4.08009e-17, &
  99 ! 56500-56000.5 cm-1
 100        3.16903e-22,1.98251e-21,5.87326e-22,3.44057e-21,2.53094e-21,8.81484e-21, &
 101        8.82299e-21,4.17179e-20,2.64703e-20,2.43792e-19,8.73831e-20,1.46371e-18, &
 102 ! 56000-55500.5 cm-1
 103        1.64421e-23,9.26011e-22,2.73137e-22,1.33640e-21,9.79188e-22,2.99706e-21, &
 104        3.37768e-21,1.39438e-20,1.47898e-20,1.04322e-19,4.08014e-20,6.31023e-19, &
 105 ! 55500-55000.5 cm-1
 106        8.68729e-24,7.31056e-22,8.78313e-23,1.07173e-21,8.28170e-22,2.54986e-21, &
 107        2.57643e-21,9.42698e-21,9.92377e-21,5.21402e-20,3.34301e-20,2.91785e-19, &
 108 ! 55000-54500.5 cm-1
 109        1.20679e-24,2.44092e-22,2.64326e-23,4.03998e-22,2.53514e-22,8.53166e-22, &
 110        1.29834e-21,3.74482e-21,5.12103e-21,2.65798e-20,2.10948e-20,2.35315e-19, &
 111 ! 54500-54000.5 cm-1
 112        2.79656e-24,1.40820e-22,3.60824e-23,2.69510e-22,4.02850e-22,8.83735e-22, &
 113        1.77198e-21,6.60221e-21,9.60992e-21,8.13558e-20,4.95591e-21,1.22858e-17, &
 114 ! 54000-53500.5 cm-1
 115        2.36959e-24,1.07535e-22,2.83333e-23,2.16789e-22,3.35242e-22,6.42753e-22, &
 116        1.26395e-21,5.43183e-21,4.88083e-21,5.42670e-20,3.27481e-21,1.58264e-17, &
 117 ! 53500-53000.5 cm-1
 118        8.65018e-25,3.70310e-23,1.04351e-23,6.43574e-23,1.17431e-22,2.70904e-22, &
 119        4.88705e-22,1.65505e-21,2.19776e-21,2.71172e-20,2.65257e-21,2.13945e-17, &
 120 ! 53000-52500.5 cm-1
 121        9.63263e-25,1.54249e-23,4.78065e-24,2.97642e-23,6.40637e-23,1.46464e-22, &
 122        1.82634e-22,7.12786e-22,1.64805e-21,2.37376e-17,9.33059e-22,1.13741e-20, &
 123 ! 52500-52000.5 cm-1
 124        1.08414e-24,8.37560e-24,9.15550e-24,2.99295e-23,9.38405e-23,1.95845e-22, &
 125        2.84356e-22,3.39699e-21,1.94524e-22,2.72227e-19,1.18924e-21,3.20246e-17, &
 126 ! 52000-51500.5 cm-1
 127        1.52817e-24,1.01885e-23,1.22946e-23,4.16517e-23,9.01287e-23,2.34869e-22, &
 128        1.93510e-22,1.44956e-21,1.81051e-22,5.17773e-21,9.82059e-22,6.22768e-17, &
 129 ! 51500-51000.5 cm-1
 130        2.12813e-24,8.48035e-24,5.23338e-24,1.93052e-23,1.99464e-23,7.48997e-23, &
 131        4.96642e-22,6.15691e-17,4.47504e-23,2.76004e-22,8.26788e-23,1.65278e-21, &
 132 ! 51000-50500.5 cm-1
 133        3.81336e-24,7.32307e-24,5.60549e-24,2.04651e-23,3.36883e-22,6.15708e-17, &
 134        2.09877e-23,1.07474e-22,9.13562e-24,8.41252e-22,0.00000e+00,0.00000e+00, &
 135 ! 50500-50000.5 cm-1
 136        5.75373e-24,7.15986e-24,5.90031e-24,3.05375e-23,2.97196e-22,8.92000e-17, &
 137        8.55920e-24,1.66709e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
 138 ! 50000-49500.5 cm-1
 139        6.21281e-24,7.13108e-24,3.30780e-24,2.61196e-23,1.30783e-22,9.42550e-17, &
 140        2.69241e-24,1.46500e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
 141 ! 49500-49000.5 cm-1
 142        6.81118e-24,6.98767e-24,7.55667e-25,2.75124e-23,1.94044e-22,1.45019e-16, &
 143        1.92236e-24,3.73223e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 /
 144
 145 !-----------------------------------------------------------------------------
 146 !     data for set_o2_xsect
 147 !-----------------------------------------------------------------------------
 148       integer, public,  parameter :: nwint = 105
 149
 150       real(dp), public, save  :: wlint(nwint)
 151       real(dp), private, save :: xso2int(nwint)
 152       real(dp), private, save :: wlla(nla)
 153       real(dp), private, save :: wlgast(ngast)
 154
 155 !-----------------------------------------------------------------------------
 156 ! o2 parameters
 157 !-----------------------------------------------------------------------------
 158       data (wlint(ii),ii=1,105)/  &
 159             0.0000E+00, 0.1166E+03, 0.1167E+03, 0.1173E+03, 0.1179E+03,  &
 160             0.1187E+03, 0.1194E+03, 0.1202E+03, 0.1208E+03, 0.1214E+03,  &
 161             0.1219E+03, 0.1223E+03, 0.1231E+03, 0.1238E+03, 0.1246E+03,  &
 162             0.1254E+03, 0.1262E+03, 0.1270E+03, 0.1278E+03, 0.1286E+03,  &
 163             0.1294E+03, 0.1303E+03, 0.1311E+03, 0.1320E+03, 0.1329E+03,  &
 164             0.1338E+03, 0.1346E+03, 0.1356E+03, 0.1365E+03, 0.1374E+03,  &
 165             0.1384E+03, 0.1399E+03, 0.1418E+03, 0.1439E+03, 0.1459E+03,  &
 166             0.1481E+03, 0.1504E+03, 0.1526E+03, 0.1550E+03, 0.1574E+03,  &
 167             0.1600E+03, 0.1626E+03, 0.1653E+03, 0.1681E+03, 0.1709E+03,  &
 168             0.1731E+03, 0.1746E+03, 0.1754E+03, 0.1770E+03, 0.1786E+03,  &
 169             0.1802E+03, 0.1818E+03, 0.1835E+03, 0.1852E+03, 0.1869E+03,  &
 170             0.1887E+03, 0.1905E+03, 0.1923E+03, 0.1942E+03, 0.1961E+03,  &
 171             0.1980E+03, 0.2000E+03, 0.2020E+03, 0.2041E+03, 0.2050E+03,  &
 172             0.2060E+03, 0.2070E+03, 0.2080E+03, 0.2090E+03, 0.2100E+03,  &
 173             0.2110E+03, 0.2120E+03, 0.2130E+03, 0.2140E+03, 0.2150E+03,  &
 174             0.2160E+03, 0.2170E+03, 0.2180E+03, 0.2190E+03, 0.2200E+03,  &
 175             0.2210E+03, 0.2220E+03, 0.2230E+03, 0.2240E+03, 0.2250E+03,  &
 176             0.2260E+03, 0.2270E+03, 0.2280E+03, 0.2290E+03, 0.2300E+03,  &
 177             0.2310E+03, 0.2320E+03, 0.2330E+03, 0.2340E+03, 0.2350E+03,  &
 178             0.2360E+03, 0.2370E+03, 0.2380E+03, 0.2390E+03, 0.2400E+03,  &
 179             0.2410E+03, 0.2420E+03, 0.2430E+03, 0.2430E+03, 0.1000E+39 /
 180
 181 !-----------------------------------------------------------------------------
 182 ! o2 parameters
 183 !-----------------------------------------------------------------------------
 184       data (xso2int(ii),ii=1,105)/  &
 185             0.0000E+00, 0.1000E-19, 0.6350E-18, 0.7525E-18, 0.1425E-18,  &
 186             0.2025E-18, 0.2413E-17, 0.6400E-17, 0.5251E-17, 0.7329E-18,  &
 187             0.3445E-18, 0.3425E-18, 0.3925E-18, 0.8918E-17, 0.9197E-17,  &
 188             0.6625E-18, 0.2700E-18, 0.1575E-18, 0.3240E-18, 0.4990E-18,  &
 189             0.4875E-18, 0.5525E-18, 0.1068E-17, 0.1850E-17, 0.2275E-17,  &
 190             0.3425E-17, 0.5890E-17, 0.8365E-17, 0.1090E-16, 0.1275E-16,  &
 191             0.1340E-16, 0.1380E-16, 0.1440E-16, 0.1445E-16, 0.1350E-16,  &
 192             0.1220E-16, 0.1071E-16, 0.9075E-17, 0.7410E-17, 0.5775E-17,  &
 193             0.4210E-17, 0.2765E-17, 0.1655E-17, 0.9760E-18, 0.5900E-18,  &
 194             0.3660E-18, 0.2061E-18, 0.3889E-19, 0.4092E-20, 0.2273E-20,  &
 195             0.1256E-20, 0.6902E-21, 0.3750E-21, 0.2097E-21, 0.1226E-21,  &
 196             0.6508E-22, 0.4579E-22, 0.3220E-22, 0.2507E-22, 0.1941E-22,  &
 197             0.1609E-22, 0.1319E-22, 0.9907E-23, 0.7419E-23, 0.7240E-23,  &
 198             0.7090E-23, 0.6955E-23, 0.6770E-23, 0.6595E-23, 0.6375E-23,  &
 199             0.6145E-23, 0.5970E-23, 0.5805E-23, 0.5655E-23, 0.5470E-23,  &
 200             0.5240E-23, 0.5005E-23, 0.4760E-23, 0.4550E-23, 0.4360E-23,  &
 201             0.4175E-23, 0.3990E-23, 0.3780E-23, 0.3560E-23, 0.3330E-23,  &
 202             0.3095E-23, 0.2875E-23, 0.2875E-23, 0.2875E-23, 0.2700E-23,  &
 203             0.2530E-23, 0.2340E-23, 0.2175E-23, 0.2020E-23, 0.1860E-23,  &
 204             0.1705E-23, 0.1555E-23, 0.1410E-23, 0.1280E-23, 0.1160E-23,  &
 205             0.1055E-23, 0.9550E-24, 0.4500E-24, 0.0000E+00, 0.0000E+00 /
 206
 207 !-----------------------------------------------------------------------------
 208 ! the lyman-alpha grid
 209 !-----------------------------------------------------------------------------
 210       data (wlla(ii),ii=1,2)/ 121.4_dp, 121.9_dp /
 211
 212 !-----------------------------------------------------------------------------
 213 ! o2 parameters
 214 !-----------------------------------------------------------------------------
 215       data (wlgast(ii),ii=1,17)/  &
 216             0.1754E+03_dp, 0.1770E+03_dp, 0.1786E+03_dp, 0.1802E+03_dp, 0.1818E+03_dp,  &
 217             0.1835E+03_dp, 0.1852E+03_dp, 0.1869E+03_dp, 0.1887E+03_dp, 0.1905E+03_dp,  &
 218             0.1923E+03_dp, 0.1942E+03_dp, 0.1961E+03_dp, 0.1980E+03_dp, 0.2000E+03_dp,  &
 219             0.2020E+03_dp, 0.2041E+03_dp /
 220
 221 !-----------------------------------------------------------------------------
 222 ! setaer arrays
 223 !-----------------------------------------------------------------------------
 224       real(dp) :: op_cb1w(nw-1), om_cb1w(nw-1), g_cb1w(nw-1)
 225       real(dp) :: op_cb2w(nw-1), om_cb2w(nw-1), g_cb2w(nw-1)
 226       real(dp) :: op_oc1w(nw-1), g_oc1w(nw-1)
 227       real(dp) :: op_oc2w(nw-1), om_oc2w(nw-1), g_oc2w(nw-1)
 228       real(dp) :: op_antw(nw-1), g_antw(nw-1)
 229       real(dp) :: op_so4w(nw-1)
 230       real(dp) :: op_salw(nw-1)
 231
 232       CONTAINS
 233
 234       subroutine photoin( chem_opt, nz, zen, o3toms, esfact,     &
 235                           o3top, o2top, albedo, z, tlev,         &
 236                           tlay, airlev, rh, xlwc, o3,            &
 237                           acb1, acb2, aoc1, aoc2, aant,          &
 238                           aso4, asal, tauaer300, tauaer400, tauaer600, &
 239                           tauaer999, waer300, waer400, waer600,  &
 240                           waer999, gaer300, gaer400, gaer600,    &
 241                           gaer999, aer_ra_feedback, prate, adjcoe, radfld )
 242
 243       use module_wave_data, only : nw, tuv_jmax, deltaw, sflx,   &
 244                                    c20, c40, c60, c80, sq
 245 !-----------------------------------------------
 246 ! ... inter-face routine
 247 !-----------------------------------------------
 248       implicit none
 249 !----------------------------------------------------------
 250 ! ... dummy arguments
 251 !----------------------------------------------------------
 252       integer, intent(in)    :: chem_opt
 253       integer, intent(in)    :: nz
 254       real(dp), intent(in)    :: zen
 255       real(dp), intent(in)    :: o3toms
 256       real(dp), intent(in)    :: esfact
 257       real(dp), intent(in)    :: o3top
 258       real(dp), intent(in)    :: o2top
 259       real(dp), intent(in)    :: albedo(nw-1)
 260       real(dp), intent(in)    :: z(nz)
 261       real(dp), intent(in)    :: tlev(nz)
 262       real(dp), intent(in)    :: tlay(nz-1)
 263       real(dp), intent(in)    :: airlev(nz)
 264       real(dp), intent(in)    :: rh(nz)
 265       real(dp), intent(in)    :: xlwc(nz)
 266       real(dp), intent(in)    :: o3(nz)
 267       real(dp), intent(in)    :: acb1(nz)
 268       real(dp), intent(in)    :: acb2(nz)
 269       real(dp), intent(in)    :: aoc1(nz)
 270       real(dp), intent(in)    :: aoc2(nz)
 271       real(dp), intent(in)    :: aant(nz)
 272       real(dp), intent(in)    :: aso4(nz)
 273       real(dp), intent(in)    :: asal(nz)
 274 ! rajesh: declare aerosol optical properties
 275       real(dp), intent(in)    :: tauaer300(nz-1)
 276       real(dp), intent(in)    :: tauaer400(nz-1)
 277       real(dp), intent(in)    :: tauaer600(nz-1)
 278       real(dp), intent(in)    :: tauaer999(nz-1)
 279       real(dp), intent(in)    :: waer300(nz-1)
 280       real(dp), intent(in)    :: waer400(nz-1)
 281       real(dp), intent(in)    :: waer600(nz-1)
 282       real(dp), intent(in)    :: waer999(nz-1)
 283       real(dp), intent(in)    :: gaer300(nz-1)
 284       real(dp), intent(in)    :: gaer400(nz-1)
 285       real(dp), intent(in)    :: gaer600(nz-1)
 286       real(dp), intent(in)    :: gaer999(nz-1)
 287       INTEGER, INTENT(IN)     :: aer_ra_feedback
 288
 289       real(dp), intent(out)   :: prate(nz,tuv_jmax)
 290       real(dp), intent(out)   :: adjcoe(nz,tuv_jmax)
 291       real(dp), intent(out)   :: radfld(nz,nw-1)
 292 !----------------------------------------------------------
 293 ! ... local variables
 294 !----------------------------------------------------------
 295       integer :: i, j, k, iw, n
 296       real(dp) :: factor, delzint
 297       character(len=8 ) :: cdate(4)
 298       character(len=10) :: ctime(4)
 299 !----------------------------------------------------------
 300 ! ... altitude grid
 301 !----------------------------------------------------------
 302       integer  :: iz, izi
 303       real(dp) :: colinc(nz)
 304       real(dp) :: vcol(nz)
 305       real(dp) :: scol(nz)
 306       real(dp) :: to3(nz)
 307 !----------------------------------------------------------
 308 ! ... fast-j adj
 309 !----------------------------------------------------------
 310       real(dp), dimension(nz,tuv_jmax) :: adjcoe1, adjcoe2
 311 !----------------------------------------------------------
 312 ! ... solar zenith angle
 313 ! slant pathlengths in spherical geometry
 314 !----------------------------------------------------------
 315       integer  :: nid(0:nz-1)
 316       real(dp) :: dsdh(0:nz-1,nz-1)
 317 !----------------------------------------------------------
 318 ! ... extra terrestrial solar flux and earth-sun distance ^-2
 319 !----------------------------------------------------------
 320       real(dp), dimension(nw-1) :: etf, delw, xsec
 321 !--------------------------------------------------------------
 322 ! ... o2 absorption cross section
 323 !--------------------------------------------------------------
 324       real(dp) :: xso2(nw-1,nz)
 325 !--------------------------------------------------------------
 326 ! ... atmospheric optical parameters:
 327 !--------------------------------------------------------------
 328       real(dp), dimension(nz-1,nw-1) :: dtrl, dto3, dto2
 329       real(dp), dimension(nz-1,nw-1) :: dtcld, omcld, gcld
 330       real(dp), dimension(nz-1,nw-1) :: dtcb1, omcb1, gcb1
 331       real(dp), dimension(nz-1,nw-1) :: dtcb2, omcb2, gcb2
 332       real(dp), dimension(nz-1,nw-1) :: dtoc1, omoc1, goc1
 333       real(dp), dimension(nz-1,nw-1) :: dtoc2, omoc2, goc2
 334       real(dp), dimension(nz-1,nw-1) :: dtant, omant, gant
 335       real(dp), dimension(nz-1,nw-1) :: dtso4, omso4, gso4
 336       real(dp), dimension(nz-1,nw-1) :: dtsal, omsal, gsal
 337 ! rajesh: declare optical property arrays for aerosols
 338       real(dp), dimension(nz-1,nw-1) :: dtaer, omaer, gaer
 339 !--------------------------------------------------------------
 340 ! ... spectral irradiance and actinic flux (scalar irradiance):
 341 !--------------------------------------------------------------
 342       real(dp), dimension(nz,nw-1) :: radxx
 343 !-------------------------------------------------------------
 344 ! ... j-values:
 345 !-------------------------------------------------------------
 346       integer :: m
 347 !-------------------------------------------------------------
 348 ! ... location and time
 349 !-------------------------------------------------------------
 350       integer :: iyear, imonth, iday
 351       real(dp) :: dtime, ut0
 352 !-------------------------------------------------------------
 353 ! ... earth-sun distance effect
 354 !-------------------------------------------------------------
 355       etf(1:nw-1) = sflx(1:nw-1) * esfact
 356 !-------------------------------------------------------
 357 ! ... air profile and rayleigh optical depths (inter-face)
 358 !-------------------------------------------------------
 359       call setair( nz, z, airlev, dtrl, colinc, o2top )
 360 !-------------------------------------------------------------
 361 ! ... ozone optical depths (must give temperature) (inter-face)
 362 !-------------------------------------------------------------
 363       call setozo( nz, z, tlay, dto3, to3, o3, airlev, o3top, o3toms )
 364 !-------------------------------------------------------------
 365 ! ... cloud optical depths (must cloud water g/m3) (inter-face)
 366 !-------------------------------------------------------------
 367       call setcld( nz, z, xlwc, dtcld, omcld, gcld )
 368
 369 !-------------------------------------------------------------
 370 ! ... aerosol optical depths ...
 371 !-------------------------------------------------------------
 372 !      call setaer( nz, z, airlev, rh, acb1,        &
 373 !                   acb2, aoc1, aoc2, aant, aso4,   &
 374 !                   asal,                           &
 375 !                   dtcb1, omcb1, gcb1,             &
 376 !                   dtcb2, omcb2, gcb2,             &
 377 !                   dtoc1, omoc1, goc1,             &
 378 !                   dtoc2, omoc2, goc2,             &
 379 !                   dtant, omant, gant,             &
 380 !                   dtso4, omso4, gso4,             &
 381 !                   dtsal, omsal, gsal )
 382
 383 !------------------------------------------------------------
 384 ! rajesh: Conform Aerosol Optical Properties from 4 wavelengths to
 385 ! the entire spectra of FTUV
 386 !------------------------------------------------------------
 387 ! Initialize aerosol optical properties to zero
 388       dtaer(:,:) = 0._dp
 389       omaer(:,:) = 0._dp
 390        gaer(:,:) = 0._dp
 391       if(aer_ra_feedback == 1) then
 392        call aer_wrf2ftuv(nz, z, tauaer300, tauaer400,    &
 393                          tauaer600, tauaer999, waer300, &
 394                          waer400, waer600, waer999,     &
 395                          gaer300, gaer400, gaer600,     &
 396                          gaer999, dtaer, omaer, gaer    )
 397       endif
 398
 399 !------------------------------------------------------------
 400 ! ... photo-chemical and photo-biological weigting functions.
 401 ! for pchem, need to know temperature and pressure profiles.
 402 ! output:
 403 ! from pchem: sq(nw,nz,nj) - for each reaction jlabel(nj)
 404 !-------------------------------------------------------------
 405       call pchem( chem_opt, nz, tlev, airlev )
 406
 407 !-------------------------------------------------------------
 408 ! ... slant path lengths for spherical geometry
 409 !-------------------------------------------------------------
 410        call spheres( nz, z, zen, dsdh, nid )
 411
 412        call airmas( nz, z, zen, dsdh, nid, colinc, vcol, scol )
 413
 414 !---------------------------------------------------------------
 415 ! ... modification of coefficent of j-vales function of to3 and zenith
 416 !---------------------------------------------------------------
 417       if( zen <= 20.5_dp ) then
 418          call setz( nz, to3, tlev, c20, 20, adjcoe )
 419       else if( zen > 20.5_dp .and. zen < 40.5_dp ) then
 420          call setz( nz, to3, tlev, c20, 20, adjcoe1 )
 421          call setz( nz, to3, tlev, c40, 40, adjcoe2 )
 422          factor = (zen - 20.5_dp) / 20._dp
 423          adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) &
 424                                + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor
 425       else if( zen >= 40.5_dp .and. zen < 60.5_dp ) then
 426          call setz( nz, to3, tlev, c40, 40, adjcoe1 )
 427          call setz( nz, to3, tlev, c60, 60, adjcoe2 )
 428          factor = (zen - 40.5_dp) / 20._dp
 429          adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) &
 430                                + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor
 431       else if( zen >= 60.5_dp .and. zen < 80._dp ) then
 432          call setz( nz, to3, tlev, c60, 60, adjcoe1 )
 433          call setz( nz, to3, tlev, c80, 80, adjcoe2 )
 434          factor = (zen - 60.5_dp) / 19.5_dp
 435          adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) &
 436                                + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor
 437       else if( zen >= 80. ) then
 438          call setz( nz, to3, tlev, c80, 80, adjcoe )
 439       end if
 440 !------------------------------------------------------------------
 441 ! ... effective o2 optical depth (sr bands, must know zenith angle!)
 442 ! assign o2 cross section to sq(*,*,1)
 443 !------------------------------------------------------------------
 444        call set_o2_xsect( nz, z, colinc, vcol, scol, dto2, xso2 )
 445        sq(1:nw-1,1:nz,1) = xso2(1:nw-1,1:nz)
 446
 447 !-----------------------------------------------
 448 ! ... main wavelength loop
 449 !-----------------------------------------------
 450       delw(1:nw-1) = deltaw(1:nw-1) * etf(1:nw-1)
 451 !---------------------------------------------------
 452 ! ... monochromatic radiative transfer:
 453 ! outputs are fdir, fdn, fup
 454 !---------------------------------------------------
 455       call rtlink( nz,                  &
 456                    nw,                  &
 457                    zen,                 &
 458                    albedo,              &
 459                    z,                   &
 460                    nid,                 &
 461                    dsdh,                &
 462                    dtrl,                & ! opt depth of air (rayleigh)
 463                    dto3,                & ! opt depth of o3
 464                    dto2,                & ! opt depth of o2
 465                    dtcld, omcld, gcld,  & ! opt depth of cloud
 466 !                   dtcb1, omcb1, gcb1,  &
 467 !                   dtcb2, omcb2, gcb2,  &
 468 !                   dtoc1, omoc1, goc1,  &
 469 !                   dtoc2, omoc2, goc2,  &
 470 !                   dtant, omant, gant,  &
 471 !                   dtso4, omso4, gso4,  &
 472 !                   dtsal, omsal, gsal,  &
 473                    dtaer, omaer, gaer,  &  ! pass aerosol optical properties
 474                    radfld )                           ! output of abs. scart flux.
 475 !----------------------------------------------------------
 476 ! Interplation the top level
 477 !----------------------------------------------------------
 478          radfld(:,:) = max( radfld(:,:),0._dp )
 479          delzint = (z(nz-2) - z(nz-3))/(z(nz-1) - z(nz-2))
 480          do iw = 1,nw-1
 481             radfld(1,iw) = radfld(2,iw) + (radfld(2,iw)-radfld(3,iw))*delzint
 482             radfld(1,iw) = max(radfld(1,iw),radfld(2,iw))
 483          end do
 484
 485 !----------------------------------------------------------
 486 ! ... j-val calculation
 487 ! spherical irradiance (actinic flux)
 488 ! as a function of altitude
 489 ! convert to quanta s-1 nm-1 cm-2
 490 ! (1.e-4 * (wc*1e-9) / (hc = 6.62e-34 * 2.998e8))
 491 !----------------------------------------------------------
 492       do m = 1,tuv_jmax
 493         do iz = 1,nz
 494           izi = nz - iz + 1
 495           xsec(:nw-1) = sq(:nw-1,iz,m) * delw(:nw-1)
 496           prate(iz,m) = dot_product( radfld(izi,:nw-1), xsec(:nw-1) )
 497         end do
 498         prate(1:nz,m) = prate(1:nz,m) * adjcoe(1:nz,m)
 499       end do
 500
 501       end subroutine photoin
 502
 503 !----------------------------------------------------------------------------
 504 !rajesh: subroutine to convert aerosol optical properties from 4
 505 !wavelengths
 506 ! to entire spectra of FTUV
 507 !-----------------------------------------------------------------------------
 508       subroutine aer_wrf2ftuv( nzlev, z, tauaer300, tauaer400,     &
 509                                tauaer600, tauaer999,               &
 510                                waer300, waer400, waer600, waer999, &
 511                                gaer300, gaer400, gaer600, gaer999, &
 512                                dtaer, omaer, gaer )
 513        use module_wave_data, only : nw, wc
 514
 515 !----------------------------------------------------------------------
 516 ! The routine is based on aerosol treatment in module_ra_rrtmg_sw.F
 517 ! INPUT:
 518 ! nzlev: number of specified altitude levels in the working grid
 519 ! z: specified altitude working grid
 520 ! Aerosol optical properties at 300, 400, 600 and 999 nm.
 521 !   tauaer300, tauaer400, tauaer600, tauaer999: Layer AODs
 522 !   waer300, waer400, waer600, waer999: Layer SSAs
 523 !   gaer300, gaer400, gaer600, gaer999: Layer asymmetry parameters
 524
 525 ! OUTPUT:
 526 ! dtaer: Layer AOD at FTUV wavelengths
 527 ! omaer: Layer SSA at FTUV wavelengths
 528 ! gaer : Layer asymmetry parameters at FTUV wavelengths
 529 !------------------------------------------------------------------------
 530       implicit none
 531
 532 !-----------------------------------------------------------------------------
 533 !       ... Dummy arguments
 534 !-----------------------------------------------------------------------------
 535       integer, intent(in)   :: nzlev
 536       real(dp), intent(in)  :: z(nzlev)
 537       real(dp), intent(in)  :: tauaer300(nzlev-1), tauaer400(nzlev-1),    &
 538                                tauaer600(nzlev-1), tauaer999(nzlev-1)
 539       real(dp), intent(in)  :: waer300(nzlev-1), waer400(nzlev-1),        &
 540                                waer600(nzlev-1), waer999(nzlev-1)
 541       real(dp), intent(in)  :: gaer300(nzlev-1), gaer400(nzlev-1),        &
 542                                gaer600(nzlev-1), gaer999(nzlev-1)
 543 ! Output arrays
 544       real(dp), intent(out) :: dtaer(nzlev-1,nw-1),                   &
 545                                omaer(nzlev-1,nw-1), gaer(nzlev-1,nw-1)
 546
 547 ! Local Variables
 548       integer     :: k, wn, nloop
 549       real(dp)    :: ang, slope
 550       real(dp), parameter :: thresh = 1.e-9_dp
 551
 552       ang   = 0._dp
 553       slope = 0._dp
 554
 555 ! Start Calculation
 556       do wn = 1,nw-1       ! wavelength loop
 557        do k = 1,nzlev-1    ! level loop
 558
 559 ! use angstrom exponent to calculate aerosol optical depth; wc is in nm.
 560         if(tauaer300(k) .gt. thresh .and. tauaer999(k) .gt. thresh) then
 561          ang = log(tauaer300(k)/tauaer999(k))/log(0.999_dp/0.3_dp)
 562          dtaer(k,wn) = tauaer400(k)*(0.4_dp/(wc(wn)*1.e-3_dp))**ang
 563 !         print *, ang, dtaer(k,wn), tauaer600(k)*(600./wc(wn))**ang
 564
 565 ! ssa - use linear interpolation/extrapolation
 566          slope = (waer600(k)-waer400(k))/0.2_dp
 567          omaer(k,wn) = slope*((wc(wn)*1.e-3_dp)-0.6_dp)+waer600(k)
 568          if(omaer(k,wn) .lt. 0.4_dp) omaer(k,wn)=0.4_dp
 569          if(omaer(k,wn) .ge. 1.0_dp) omaer(k,wn)=1.0_dp
 570
 571 ! asymmetry parameter - use linear interpolation/extrapolation
 572          slope = (gaer600(k)-gaer400(k))/0.2_dp
 573          gaer(k,wn) = slope*((wc(wn)*1.e-3_dp)-0.6_dp)+gaer600(k)
 574          if(gaer(k,wn) .lt. 0.5_dp) gaer(k,wn) = 0.5_dp
 575          if(gaer(k,wn) .ge. 1.0_dp) gaer(k,wn) = 1.0_dp
 576         endif
 577        enddo  ! k
 578       enddo   ! wn
 579
 580       end subroutine aer_wrf2ftuv
 581
 582 !--------------------------------------------------------------------
 583
 584       subroutine setaer( nzlev, z, airden, rh, acb1,    &
 585                          acb2, aoc1, aoc2, aant, aso4,  &
 586                          asal,                          &
 587                          op_cb1, om_cb1, g_cb1,         &
 588                          op_cb2, om_cb2, g_cb2,         &
 589                          op_oc1, om_oc1, g_oc1,         &
 590                          op_oc2, om_oc2, g_oc2,         &
 591                          op_ant, om_ant, g_ant,         &
 592                          op_so4, om_so4, g_so4,         &
 593                          op_sal, om_sal, g_sal )
 594
 595       use module_wave_data, only : nw, wc
 596
 597 !-----------------------------------------------------------------------------
 598 !=  PURPOSE:
 599 !
 600 !=  Cal aer opt depth
 601 !-----------------------------------------------------------------------------
 602 !
 603 !=  PARAMETERS:
 604 !=  NZLEV   - INTEGER, number of specified altitude levels in the working (I)
 605 !=            grid
 606 !=  Z       - real(dp), specified altitude working grid (km)                  (I)
 607 !=  axxx    - aerosol mix ratio                                           (I)
 608 !
 609 !=  op      - real(dp), optical depth (absorption)                            (O)
 610 !=  om      - real(dp), single albedo
 611 !=  g       - asysmetric factor                                           (O)
 612 !
 613 !-----------------------------------------------------------------------------
 614 !
 615 !     VERTICAL DOMAIN is from bottom(1)  to TOP (TOP=nzlev)
 616 !        CCM from top(1) to bottom(nzlev)
 617 !-----------------------------------------------------------------------------
 618       implicit none
 619
 620 !-----------------------------------------------------------------------------
 621 !       ... Dummy arguments
 622 !-----------------------------------------------------------------------------
 623       integer, intent(in)   :: nzlev
 624       real(dp), intent(in)  :: z(nzlev)
 625       real(dp), intent(in)  :: airden(nzlev)
 626       real(dp), intent(in)  :: rh(nzlev)
 627       real(dp), intent(in)  :: acb1(nzlev)
 628       real(dp), intent(in)  :: acb2(nzlev)
 629       real(dp), intent(in)  :: aoc1(nzlev)
 630       real(dp), intent(in)  :: aoc2(nzlev)
 631       real(dp), intent(in)  :: aant(nzlev)
 632       real(dp), intent(in)  :: aso4(nzlev)
 633       real(dp), intent(in)  :: asal(nzlev)
 634       real(dp), intent(out) :: op_cb1(nzlev-1,nw-1), om_cb1(nzlev-1,nw-1), g_cb1(nzlev-1,nw-1)
 635       real(dp), intent(out) :: op_cb2(nzlev-1,nw-1), om_cb2(nzlev-1,nw-1), g_cb2(nzlev-1,nw-1)
 636       real(dp), intent(out) :: op_oc1(nzlev-1,nw-1), om_oc1(nzlev-1,nw-1), g_oc1(nzlev-1,nw-1)
 637       real(dp), intent(out) :: op_oc2(nzlev-1,nw-1), om_oc2(nzlev-1,nw-1), g_oc2(nzlev-1,nw-1)
 638       real(dp), intent(out) :: op_ant(nzlev-1,nw-1), om_ant(nzlev-1,nw-1), g_ant(nzlev-1,nw-1)
 639       real(dp), intent(out) :: op_so4(nzlev-1,nw-1), om_so4(nzlev-1,nw-1), g_so4(nzlev-1,nw-1)
 640       real(dp), intent(out) :: op_sal(nzlev-1,nw-1), om_sal(nzlev-1,nw-1), g_sal(nzlev-1,nw-1)
 641
 642 !-----------------------------------------------------------------------------
 643 !       ... local variables
 644 !-----------------------------------------------------------------------------
 645       integer     :: k, wn, nz
 646       real(dp)    :: wcen
 647       real(dp)    :: dz(nzlev)
 648       real(dp)    :: sig_oc2(nzlev-1), sig_cb2(nzlev-1)
 649       real(dp)    :: opt_oc1(nzlev-1), opt_oc2(nzlev-1), opt_cb1(nzlev-1), opt_cb2(nzlev-1)
 650       real(dp)    :: opt_ant(nzlev-1), opt_so4(nzlev-1), opt_sal(nzlev-1)
 651       real(dp)    :: rw(nzlev-1), num(nzlev-1), wght(nzlev-1)
 652
 653
 654 !-----------------------------------------------------------------------------
 655 !       ... assign the constants
 656 !-----------------------------------------------------------------------------
 657       real(dp), parameter :: rm_cb1 = 0.035_dp
 658       real(dp), parameter :: rm_cb2 = 0.035_dp
 659       real(dp), parameter :: rm_oc1 = 0.10_dp
 660       real(dp), parameter :: rm_oc2 = 0.10_dp
 661       real(dp), parameter :: rm_ant = 0.15_dp
 662       real(dp), parameter :: rm_so4 = 0.24_dp
 663       real(dp), parameter :: rm_sal = 1.30_dp
 664
 665       real(dp), parameter :: de_cb1 = 1.0_dp
 666       real(dp), parameter :: de_cb2 = 1.0_dp
 667       real(dp), parameter :: de_oc1 = 1.5_dp
 668       real(dp), parameter :: de_oc2 = 1.5_dp
 669       real(dp), parameter :: de_ant = 1.7_dp
 670       real(dp), parameter :: de_so4 = 1.7_dp
 671       real(dp), parameter :: de_sal = 2.2_dp
 672
 673       real(dp), parameter :: ml_cb1 = 12._dp
 674       real(dp), parameter :: ml_cb2 = 12._dp
 675       real(dp), parameter :: ml_oc1 = 48._dp
 676       real(dp), parameter :: ml_oc2 = 48._dp
 677       real(dp), parameter :: ml_ant = 90._dp
 678       real(dp), parameter :: ml_so4 = 96._dp
 679       real(dp), parameter :: ml_sal = 22._dp
 680
 681       real(dp), parameter :: qx_cb1 = 2.138_dp
 682       real(dp), parameter :: qx_cb2 = 2.138_dp
 683       real(dp), parameter :: qw_cb1 = 0.510_dp
 684       real(dp), parameter :: qw_cb2 = 0.510_dp
 685       real(dp), parameter :: sx_cb2 = 0.120_dp
 686       real(dp), parameter :: sw_cb2 = 1.000_dp
 687
 688       real(dp), parameter :: qx_oc1 = 0.675_dp
 689       real(dp), parameter :: qx_oc2 = 0.675_dp
 690       real(dp), parameter :: qw_oc1 = 1.118_dp
 691       real(dp), parameter :: qw_oc2 = 1.118_dp
 692       real(dp), parameter :: sx_oc2 = 0.980_dp
 693       real(dp), parameter :: sw_oc2 = 1.000_dp
 694
 695       real(dp), parameter :: qx_ant = 0.726_dp
 696       real(dp), parameter :: qw_ant = 1.770_dp
 697
 698       real(dp), parameter :: qx_so4 = 2.830_dp
 699       real(dp), parameter :: qw_so4 = 3.605_dp
 700
 701       real(dp), parameter :: qx_sal = 2.600_dp
 702       real(dp), parameter :: qw_sal = 2.575_dp
 703
 704       real(dp), parameter :: pi = 3.1416_dp
 705
 706 !-----------------------------------------------------------------------------
 707 !       ... hydroscropic growth
 708 !-----------------------------------------------------------------------------
 709       real(dp), parameter :: a1_cb2 = 1.000_dp
 710       real(dp), parameter :: a2_cb2 =-0.829_dp
 711       real(dp), parameter :: a3_cb2 = 1.350_dp
 712
 713       real(dp), parameter :: a1_oc2 = 1.005_dp
 714       real(dp), parameter :: a2_oc2 =-0.0666_dp
 715       real(dp), parameter :: a3_oc2 = 0.7608_dp
 716
 717       real(dp), parameter :: a1_ant = 1.000_dp
 718       real(dp), parameter :: a2_ant = 0.493_dp
 719       real(dp), parameter :: a3_ant = 0.527_dp
 720
 721       real(dp), parameter :: a1_so4 = 1.000_dp
 722       real(dp), parameter :: a2_so4 = 0.493_dp
 723       real(dp), parameter :: a3_so4 = 0.528_dp
 724
 725       real(dp), parameter :: a1_sal = 1.011_dp
 726       real(dp), parameter :: a2_sal = 0.457_dp
 727       real(dp), parameter :: a3_sal = 1.102_dp
 728
 729       nz = nzlev - 1
 730 !-----------------------------------------------------------------------------
 731 !       ... calculate vertical interveral dz
 732 !-----------------------------------------------------------------------------
 733       dz(1:nz) = abs( z(2:nz+1) - z(1:nz) ) * km2cm * pi * 1.e-8_dp
 734 !-----------------------------------------------------------------------------
 735 !        -------- CB1 -------------
 736 !-----------------------------------------------------------------------------
 737          rw(:nz) = rm_cb1
 738          call cnum( acb1, airden, de_cb1, ml_cb1, rm_cb1, &
 739                     rw, num, wght, nzlev )
 740
 741          opt_cb1(:nz) = qx_cb1*wght(:nz) + (1._dp - wght(:nz))*qw_cb1
 742          opt_cb1(:nz) = opt_cb1(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
 743 !-----------------------------------------------------------------------------
 744 !        -------- CB2 -------------
 745 !-----------------------------------------------------------------------------
 746          rw(:nz) = rm_cb2*(a1_cb2 + rh(:nz)*(a2_cb2 + a3_cb2*rh(:nz)))
 747          call cnum( acb2, airden, de_cb2, ml_cb2, rm_cb2, &
 748                     rw, num, wght, nzlev )
 749
 750          opt_cb2(:nz) = qx_cb2*wght(:nz) + (1._dp - wght(:nz))*qw_cb2
 751          opt_cb2(:nz) = opt_cb2(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
 752          sig_cb2(:nz) = sx_cb2*wght(:nz) + (1._dp - wght(:nz))*sw_cb2
 753 !-----------------------------------------------------------------------------
 754 !        -------- OC1 -------------
 755 !-----------------------------------------------------------------------------
 756          rw(:nz) = rm_oc1
 757          call cnum( aoc1, airden, de_oc1, ml_oc1, rm_oc1, &
 758                     rw, num, wght, nzlev )
 759
 760          opt_oc1(:nz) = qx_oc1*wght(:nz) + (1._dp - wght(:nz))*qw_oc1
 761          opt_oc1(:nz) = opt_oc1(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
 762 !-----------------------------------------------------------------------------
 763 !        -------- OC2 -------------
 764 !-----------------------------------------------------------------------------
 765          rw(:nz) = rm_oc2*(a1_oc2 + rh(:nz)*(a2_oc2 + a3_oc2*rh(:nz)))
 766          call cnum( aoc2, airden, de_oc2, ml_oc2, rm_oc2, &
 767                     rw, num, wght, nzlev )
 768
 769          opt_oc2(:nz) = qx_oc2*wght(:nz) + (1._dp - wght(:nz))*qw_oc2
 770          opt_oc2(:nz) = opt_oc2(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
 771          sig_oc2(:nz) = sx_oc2*wght(:nz) + (1._dp - wght(:nz))*sw_oc2
 772 !-----------------------------------------------------------------------------
 773 !        -------- ANT-NH4NO3 -------------
 774 !-----------------------------------------------------------------------------
 775          rw(:nz) = rm_ant*(a1_ant + rh(:nz)*(a2_ant + a3_ant*rh(:nz)))
 776          call cnum( aant, airden, de_ant, ml_ant, rm_ant, &
 777                     rw, num, wght, nzlev )
 778
 779          opt_ant(:nz) = qx_ant*wght(:nz) + (1._dp - wght(:nz))*qw_ant
 780          opt_ant(:nz) = opt_ant(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
 781 !-----------------------------------------------------------------------------
 782 !        -------- SO4 -------------
 783 !-----------------------------------------------------------------------------
 784          rw(:nz) = rm_so4*(a1_so4 + rh(:nz)*(a2_so4 + a3_so4*rh(:nz)))
 785          call cnum( aso4, airden, de_so4, ml_so4, rm_so4, &
 786                     rw, num, wght, nzlev )
 787
 788          opt_so4(:nz) = qx_so4*wght(:nz) + (1._dp - wght(:nz))*qw_so4
 789          opt_so4(:nz) = opt_so4(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
 790 !-----------------------------------------------------------------------------
 791 !        -------- SALT -------------
 792 !-----------------------------------------------------------------------------
 793          rw(:nz) = rm_sal*(a1_sal + rh(:nz)*(a2_sal + a3_sal*rh(:nz)))
 794          call cnum( asal, airden, de_sal, ml_sal, rm_sal, &
 795                     rw, num, wght, nzlev )
 796
 797          opt_sal(:nz) = qx_sal*wght(:nz) + (1._dp - wght(:nz))*qw_sal
 798          opt_sal(:nz) = opt_sal(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz)
 799
 800 WAVE_LENGTH_LOOP : &
 801          do wn = 1,nw-1
 802             wcen = wc(wn)
 803             op_cb1(:nz,wn) = opt_cb1(:nz)*op_cb1w(wn)
 804             om_cb1(:nz,wn) = om_cb1w(wn)
 805              g_cb1(:nz,wn) = g_cb1w(wn)
 806
 807             op_cb2(:nz,wn) = opt_cb2(:nz)*op_cb2w(wn)
 808             om_cb2(:nz,wn) = sig_cb2(:nz)*om_cb2w(wn)
 809              g_cb2(:nz,wn) = g_cb2w(wn)
 810
 811             op_oc1(:nz,wn) = opt_oc1(:nz)*op_oc1w(wn)
 812             om_oc1(:nz,wn) = 0.98_dp
 813              g_oc1(:nz,wn) = g_oc1w(wn)
 814
 815             op_oc2(:nz,wn) = opt_oc2(:nz)*op_oc2w(wn)
 816             om_oc2(:nz,wn) = sig_oc2(:nz)
 817              g_oc2(:nz,wn) = g_oc2w(wn)
 818
 819             op_ant(:nz,wn) = opt_ant(:nz)*op_antw(wn)
 820             om_ant(:nz,wn) = 1.00_dp
 821              g_ant(:nz,wn) = g_antw(wn)
 822
 823             op_so4(:nz,wn) = opt_so4(:nz)*op_so4w(wn)
 824             om_so4(:nz,wn) = 1.00_dp
 825              g_so4(:nz,wn) = 0.75_dp
 826
 827             op_sal(:nz,wn) = opt_sal(:nz)*op_salw(wn)
 828             om_sal(:nz,wn) = 1.00_dp
 829              g_sal(:nz,wn) = 0.77_dp
 830          enddo WAVE_LENGTH_LOOP
 831
 832       end subroutine setaer
 833
 834       subroutine aer_init
 835 !-----------------------------------------------------------------------------
 836 !   ... initialize setaer
 837 !-----------------------------------------------------------------------------
 838
 839       use module_wave_data, only : nw, wc
 840
 841       implicit none
 842
 843 !-----------------------------------------------------------------------
 844 !     Assign the wave-length dep values for aerososl
 845 !-----------------------------------------------------------------------
 846       real(dp), parameter :: d1_cb = 1.61892_dp
 847       real(dp), parameter :: d2_cb =-0.0006853_dp
 848       real(dp), parameter :: d3_cb =-1.07806e-6_dp
 849
 850       real(dp), parameter :: o1_cb = 1.29185_dp
 851       real(dp), parameter :: o2_cb = 7.6863e-5_dp
 852       real(dp), parameter :: o3_cb =-1.28567e-6_dp
 853
 854       real(dp), parameter :: g1_cb = 2.87326_dp
 855       real(dp), parameter :: g2_cb =-0.004482_dp
 856       real(dp), parameter :: g3_cb = 1.50361e-6_dp
 857
 858       real(dp), parameter :: d1_oc = 10.5098_dp
 859       real(dp), parameter :: d2_oc =-0.0301455_dp
 860       real(dp), parameter :: d3_oc = 2.22837e-5_dp
 861
 862       real(dp), parameter :: g1_oc = 1.67125_dp
 863       real(dp), parameter :: g2_oc =-0.0008341_dp
 864       real(dp), parameter :: g3_oc =-1.18801e-6_dp
 865
 866       real(dp), parameter :: d1_ant = 9.44794_dp
 867       real(dp), parameter :: d2_ant =-0.0268029_dp
 868       real(dp), parameter :: d3_ant = 1.97106e-5_dp
 869
 870       real(dp), parameter :: g1_ant = 1.34742_dp
 871       real(dp), parameter :: g2_ant =-1.85850e-5_dp
 872       real(dp), parameter :: g3_ant =-1.54297e-6_dp
 873
 874       real(dp), parameter :: d1_so4 = 0.942351_dp
 875       real(dp), parameter :: d2_so4 = 0.00220168_dp
 876       real(dp), parameter :: d3_so4 =-3.9622e-6_dp
 877
 878       real(dp), parameter :: d1_sal = 1.0_dp
 879       real(dp), parameter :: d2_sal = 0.0_dp
 880       real(dp), parameter :: d3_sal = 0.0_dp
 881
 882       op_cb1w(:nw-1) = d1_cb + wc(:nw-1)*(d2_cb + d3_cb*wc(:nw-1))
 883       om_cb1w(:nw-1) = .12_dp*(o1_cb + wc(:nw-1)*(o2_cb + o3_cb*wc(:nw-1)))
 884       g_cb1w(:nw-1)  = .12_dp*(g1_cb + wc(:nw-1)*(g2_cb + g3_cb*wc(:nw-1)))
 885
 886       op_cb2w(:nw-1) = op_cb1w(:nw-1)
 887       om_cb2w(:nw-1) = om_cb1w(:nw-1)
 888       g_cb2w(:nw-1)  = g_cb1w(:nw-1)
 889
 890       op_oc1w(:nw-1) = d1_oc + wc(:nw-1)*(d2_oc + d3_oc*wc(:nw-1))
 891       g_oc1w(:nw-1)  = 0.50_dp*(g1_oc + wc(:nw-1)*(g2_oc + g3_oc*wc(:nw-1)))
 892
 893       op_oc2w(:nw-1) = op_oc1w(:nw-1)
 894 !     om_oc2w(:nw-1) = sig_oc2
 895       g_oc2w(:nw-1)  = g_oc1w(:nw-1)
 896
 897       op_antw(:nw-1) = d1_ant + wc(:nw-1)*(d2_ant + d3_ant*wc(:nw-1))
 898       g_antw(:nw-1)  = 0.63_dp*(g1_ant + wc(:nw-1)*(g2_ant + g3_ant*wc(:nw-1)))
 899
 900       op_so4w(:nw-1) = d1_so4 + wc(:nw-1)*(d2_so4 + d3_so4*wc(:nw-1))
 901
 902       op_salw(:nw-1) = d1_sal + wc(:nw-1)*(d2_sal + d3_sal*wc(:nw-1))
 903
 904       end subroutine aer_init
 905
 906       subroutine cnum( mix, denx, cden, molw, rm, &
 907                        rw, num, weight, nzlev )
 908 !-----------------------------------------------------------------------------
 909 !   PURPOSE:
 910 !
 911 !   Cal aer num and percent weigh
 912 !-----------------------------------------------------------------------------
 913       implicit none
 914 !-----------------------------------------------------------------------------
 915 !       ... Dummy arguments
 916 !-----------------------------------------------------------------------------
 917       integer,  intent(in)  :: nzlev
 918       real(dp), intent(in)  :: molw, rm
 919       real(dp), intent(in)  :: cden
 920       real(dp), intent(in)  :: rw(nzlev-1)
 921       real(dp), intent(in)  :: mix(nzlev-1), denx(nzlev-1)
 922       real(dp), intent(out) :: num(nzlev-1), weight(nzlev-1)
 923 !-----------------------------------------------------------------------------
 924 !       ... Local variables
 925 !-----------------------------------------------------------------------------
 926       REAL(dp), PARAMETER :: PI = 3.14159265358979324
 927       REAL(dp), PARAMETER :: vol_fac = 4._dp * PI / 3._dp
 928       integer  :: k
 929       real(dp) :: rm_cubic, factor
 930       real(dp) :: mas(nzlev-1), masw(nzlev-1)
 931
 932       rm_cubic = rm**3
 933       factor   = 1._dp/(vol_fac * cden * 1.e-12_dp * rm_cubic)
 934 !-----------------------------------------------------------------------------
 935 !     Convert volumn mix ratio to mass
 936 !-----------------------------------------------------------------------------
 937       where( mix(:nzlev-1) >= 1.e-15_dp )
 938 !===============================================================================
 939 ! XUEXI
 940 !     UNLIKE MOZART mix = mix ratio
 941 !                   mix = ug/m3
 942 !     1e-12 for  ug/m3 --> gram/cm3
 943 !===============================================================================
 944          mas(:nzlev-1)    = mix(:nzlev-1) * 1.e-12_dp                              ! 1e-12 for  ug/m3 --> gram/cm3
 945          num(:nzlev-1)    = mas(:nzlev-1) * factor
 946       elsewhere
 947         num(:nzlev-1)    = 0._dp
 948         weight(:nzlev-1) = 0._dp
 949       endwhere
 950
 951       do k = 1,nzlev-1
 952          if( abs(rw(k) - rm) >= 1.e-6_dp ) then
 953            masw(k) = vol_fac*(rw(k)**3 - rm_cubic)*1.e-12_dp * num(k)   ! water mass gram/cm3
 954          else
 955            masw(k) = 0._dp
 956          end if
 957       end do
 958
 959       where( mix(:nzlev-1) >= 1.e-15_dp )
 960          weight(:nzlev-1) = mas(:nzlev-1)/( mas(:nzlev-1) + masw(:nzlev-1) )        ! weight percentage (%*0.01)
 961          weight(:nzlev-1) = max( 0._dp, min( 1._dp,weight(:nzlev-1) ) )
 962       endwhere
 963
 964       end subroutine cnum
 965
 966       subroutine setair( nz, z, airlev, dtrl, &
 967                          cz, o2top )
 968
 969       use module_wave_data, only : nw, wc
 970 !-----------------------------------------------------------------------------
 971 ! purpose:
 972 ! set up an altitude profile of air molecules. subroutine includes a
 973 ! shape-conserving scaling method that allows scaling of the entire
 974 ! profile to a given sea-level pressure.
 975 !-----------------------------------------------------------------------------
 976 ! parameters:
 977 ! pmbnew - real(dp), sea-level pressure (mb) to which profile should be
 978 ! scaled. if pmbnew < 0, no scaling is done
 979 ! nz - integer, number of specified altitude levels in the working (i)
 980 ! grid
 981 ! z - real(dp), specified altitude working grid (km) (i)
 982 ! nw - integer, number of specified intervals + 1 in working (i)
 983 ! wavelength grid
 984 ! wl - real(dp), vector of lower limits of wavelength intervals in (i)
 985 ! working wavelength grid
 986 ! airlev - real(dp), air density (molec/cc) at each specified altitude (o)
 987 ! dtrl - real(dp), rayleigh optical depth at each specified altitude (o)
 988 ! and each specified wavelength
 989 ! cz - real(dp), number of air molecules per cm^2 at each specified (o)
 990 ! altitude layer
 991 !-----------------------------------------------------------------------------
 992       implicit none
 993 !-----------------------------------------------------------------------------
 994 ! ... dummy arguments
 995 !-----------------------------------------------------------------------------
 996       integer, intent(in) :: nz
 997       real(dp), intent(in) :: o2top
 998       real(dp), intent(in) :: z(nz)
 999       real(dp), intent(in) :: airlev(nz)
1000 !-----------------------------------------------------------------------------
1001 ! ... air density (molec cm-3) at each grid level
1002 ! rayleigh optical depths
1003 !-----------------------------------------------------------------------------
1004       real(dp), intent(out) :: dtrl(nz-1,nw-1)
1005       real(dp), intent(out) :: cz(nz)
1006 !-----------------------------------------------------------------------------
1007 ! ... local variables
1008 !-----------------------------------------------------------------------------
1009       integer :: i, ip1, iw
1010       real(dp) :: hscale
1011       real(dp) :: srayl
1012       real(dp) :: deltaz
1013       real(dp) :: wmicrn, xx
1014 !-----------------------------------------------------
1015 ! ... compute column increments (logarithmic integrals)
1016 !-----------------------------------------------------
1017       do i = 1,nz-1
1018          ip1 = i + 1
1019          deltaz = km2cm * (z(ip1) - z(i))
1020          cz(i) = (airlev(ip1) - airlev(i)) / log(airlev(ip1)/airlev(i)) * deltaz
1021       end do
1022 !-----------------------------------------------------
1023 ! ... include exponential tail integral from infinity to 50 km,
1024 ! fold tail integral into top layer
1025 ! specify scale height near top of data. (scale height at 40km????)
1026 !-----------------------------------------------------
1027 ! hscale = 8.05e5
1028       cz(nz) = o2top/.2095_dp
1029 !-----------------------------------------------------
1030 ! ... compute rayleigh cross sections and depths:
1031 !-----------------------------------------------------
1032       do iw = 1,nw-1
1033 !-----------------------------------------------------
1034 ! ... rayleigh scattering cross section from wmo 1985 (originally from
1035 ! nicolet, m., on the molecular scattering in the terrestrial atmosphere:
1036 ! an empirical formula for its calculation in the homoshpere, planet.
1037 ! space sci., 32, 1467-1468, 1984.
1038 !-----------------------------------------------------
1039          wmicrn = 1.e-3_dp*wc(iw)
1040          if( wmicrn <= .55_dp ) then
1041             xx = 3.6772_dp + 0.389_dp*wmicrn + 0.09426_dp/wmicrn
1042          else
1043             xx = 4.04_dp
1044          end if
1045          srayl = 4.02e-28_dp/(wmicrn)**xx
1046          dtrl(1:nz-2,iw) = cz(1:nz-2)*srayl
1047          dtrl(nz-1,iw) = (cz(nz-1) + cz(nz))*srayl
1048       end do
1049
1050       end subroutine setair
1051
1052       subroutine zenith( lat, long, idate, ut, zen )
1053 !-----------------------------------------------------------------------------
1054 ! purpose:
1055 ! calculate solar zenith angle for a given time and location.
1056 ! calculation is based on equations given in: paltridge and platt, radia-
1057 ! tive processes in meteorology and climatology, elsevier, pp. 62,63, 1976.
1058 ! fourier coefficients originally from: spencer, j.w., 1971, fourier
1059 ! series representation of the position of the sun, search, 2:172.
1060 ! note: this approximate program does not account fro changes from year
1061 ! to year.
1062 !-----------------------------------------------------------------------------
1063 ! parameters:
1064 ! lat - real(dp), latitude of location (degrees) (i)
1065 ! long - real(dp), longitude of location (degrees) (i)
1066 ! idate - integer, date in the form yymmdd (i)
1067 ! ut - real(dp), local time in decimal ut (e.g., 16.25 means 15 minutes (i)
1068 ! after 4 pm)
1069 ! zen - real(dp), solar zenith angle (degrees) (o)
1070 !-----------------------------------------------------------------------------
1071       implicit none
1072 !-----------------------------------------------------------------------------
1073 ! ... dummy arguments
1074 !-----------------------------------------------------------------------------
1075       integer, intent(in) :: idate
1076       real(dp), intent(in) :: lat,long
1077       real(dp), intent(in) :: ut
1078       real(dp), intent(out) :: zen
1079 !-----------------------------------------------------------------------------
1080 ! ... local variables
1081 !-----------------------------------------------------------------------------
1082       real(dp), parameter :: pi = 3.1415926
1083       real(dp), parameter :: d2r = pi/180.0
1084       real(dp), parameter :: r2d = 1.0/d2r
1085       integer :: i
1086       integer :: iiyear, imth, iday, ijd
1087       integer :: imn(12) = (/ 31,28,31,30,31,30,31,31,30,31,30,31 /)
1088       real(dp) :: lbut,lzut
1089       real(dp) :: rlt
1090       real(dp) :: d, tz, rdecl, eqr, eqh, zpt
1091       real(dp) :: csz, zr, caz, raz
1092       real(dp) :: sintz, costz, sin2tz, cos2tz, sin3tz, cos3tz
1093 !-----------------------------------------------------------------------------
1094 ! ... convert to radians
1095 !-----------------------------------------------------------------------------
1096       rlt = lat*d2r
1097 !-----------------------------------------------------------------------------
1098 ! ... parse date
1099 !-----------------------------------------------------------------------------
1100       iiyear = idate/10000
1101       imth = (idate - iiyear*10000)/100
1102       iday = idate - iiyear*10000 - imth*100
1103 !-----------------------------------------------------------------------------
1104 ! ... identify and correct leap years
1105 !-----------------------------------------------------------------------------
1106       if( mod(iiyear,4) == 0 ) then
1107          imn(2) = 29
1108       else
1109          imn(2) = 28
1110       end if
1111 !-----------------------------------------------------------------------------
1112 ! ... compute current (julian) day of year ijd = 1 to 365
1113 !-----------------------------------------------------------------------------
1114       ijd = 0
1115       do i = 1,imth-1
1116          ijd = ijd + imn(i)
1117       end do
1118       ijd = ijd + iday
1119 !-----------------------------------------------------------------------------
1120 ! ... calculate decimal julian day from start of year:
1121 !-----------------------------------------------------------------------------
1122       d = real(ijd-1,kind=dp) + ut/24._dp
1123 !-----------------------------------------------------------------------------
1124 ! ... equation 3.8 for "day-angle"
1125 !-----------------------------------------------------------------------------
1126       tz = 2._dp*pi*d/365._dp
1127 !-----------------------------------------------------------------------------
1128 ! ... calculate sine and cosine from addition theoremes for
1129 ! better performance; the computation of sin2tz,
1130 ! sin3tz, cos2tz and cos3tz is about 5-6 times faster
1131 ! than the evaluation of the intrinsic functions
1132 !
1133 ! it is sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y)
1134 ! and cos(x+y) = cos(x)*cos(y)-sin(x)*sin(y)
1135 !
1136 ! sintz = sin(tz) costz = cos(tz)
1137 ! sin2tz = sin(2.*tz) cos2tz = sin(2.*tz)
1138 ! sin3tz = sin(3.*tz) cos3tz = cos(3.*tz)
1139 !-----------------------------------------------------------------------------
1140       sintz = sin( tz )
1141       costz = cos( tz )
1142       sin2tz = 2._dp*sintz*costz
1143       cos2tz = costz*costz - sintz*sintz
1144       sin3tz = sintz*cos2tz + costz*sin2tz
1145       cos3tz = costz*cos2tz - sintz*sin2tz
1146 !-----------------------------------------------------------------------------
1147 ! ... equation 3.7 for declination in radians
1148 !-----------------------------------------------------------------------------
1149       rdecl = 0.006918_dp - 0.399912_dp*costz + 0.070257_dp*sintz &
1150                           - 0.006758_dp*cos2tz + 0.000907_dp*sin2tz &
1151                           - 0.002697_dp*cos3tz + 0.001480_dp*sin3tz
1152 !-----------------------------------------------------------------------------
1153 ! ... equation 3.11 for equation of time in radians
1154 !-----------------------------------------------------------------------------
1155       eqr = 0.000075_dp + 0.001868_dp*costz - 0.032077_dp*sintz &
1156                         - 0.014615_dp*cos2tz - 0.040849_dp*sin2tz
1157 !-----------------------------------------------------------------------------
1158 ! ... convert equation of time to hours
1159 !-----------------------------------------------------------------------------
1160       eqh = eqr*24._dp/(2._dp*pi)
1161 !-----------------------------------------------------------------------------
1162 ! ... calculate local hour angle (hours):
1163 !-----------------------------------------------------------------------------
1164       lbut = 12._dp - eqh - long*24._dp/360._dp
1165 !-----------------------------------------------------------------------------
1166 ! ... convert to angle from ut
1167 !-----------------------------------------------------------------------------
1168       lzut = 15._dp*(ut - lbut)
1169       zpt = lzut*d2r
1170 !-----------------------------------------------------------------------------
1171 ! ... equation 2.4 for cosine of zenith angle
1172 !-----------------------------------------------------------------------------
1173       csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
1174       zr = acos(csz)
1175       zen = zr*r2d
1176
1177       end subroutine zenith
1178
1179       subroutine calc_zenith( lat, long, ijd, gmt, azimuth, zenith )
1180 !-------------------------------------------------------------------
1181 ! this subroutine calculates solar zenith and azimuth angles for a
1182 ! part  time and location.  must specify:
1183 ! input:
1184 ! lat - latitude in decimal degrees
1185 ! long - longitude in decimal degrees
1186 ! gmt  - greenwich mean time - decimal military eg.
1187 ! 22.75 = 45 min after ten pm gmt
1188 ! output:
1189 ! zenith
1190 ! azimuth
1191 ! .. Scalar Arguments ..
1192 !-------------------------------------------------------------------
1193         integer, intent(in)   :: ijd
1194         real(dp), intent(in)  :: gmt, lat, long
1195         real(dp), intent(out) :: azimuth, zenith
1196
1197 !-------------------------------------------------------------------
1198 ! .. Local variables
1199 !-------------------------------------------------------------------
1200       real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp
1201       real(dp), parameter :: r2d = 1.0_dp/d2r
1202
1203       integer  :: jd
1204       real(dp) :: caz, csz, cw, d, decl, ec, epsi, eqt, eyt, feqt, feqt1,       &
1205           feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi,     &
1206           pi, ra, raz, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr,    &
1207           yt, zpt, zr
1208 !-------------------------------------------------------------------
1209 ! .. Intrinsic Functions ..
1210 !-------------------------------------------------------------------
1211         intrinsic acos, atan, cos, min, sin, tan
1212
1213 !-------------------------------------------------------------------
1214 ! convert to radians
1215 !-------------------------------------------------------------------
1216         rlt = lat*d2r
1217         rphi = long*d2r
1218
1219         jd = ijd
1220
1221         d = jd + gmt/24.0_dp
1222 !-------------------------------------------------------------------
1223 ! calc geom mean longitude
1224 !-------------------------------------------------------------------
1225         ml = 279.2801988_dp + d*(.9856473354_dp + 2.267E-13_dp*d)
1226         rml = ml*d2r
1227
1228 !-------------------------------------------------------------------
1229 ! calc equation of time in sec
1230 ! w = mean long of perigee
1231 ! e = eccentricity
1232 ! epsi = mean obliquity of ecliptic
1233 !-------------------------------------------------------------------
1234         w = 282.4932328_dp + d*(4.70684E-5_dp + 3.39E-13_dp*d)
1235         wr = w*d2r
1236         ec   = 1.6720041E-2_dp - d*(1.1444E-9_dp + 9.4E-17_dp*d)
1237         epsi = 23.44266511_dp - d*(3.5626E-7_dp + 1.23E-15_dp*d)
1238         pepsi = epsi*d2r
1239         yt = (tan(pepsi/2.0_dp))**2
1240         cw = cos(wr)
1241         sw = sin(wr)
1242         ssw = sin(2.0_dp*wr)
1243         eyt = 2._dp*ec*yt
1244         feqt1 = -sin(rml)*cw*(eyt + 2._dp*ec)
1245         feqt2 = cos(rml)*sw*(2._dp*ec - eyt)
1246         feqt3 = sin(2._dp*rml)*(yt - (5._dp*ec**2/4._dp)*(cw**2 - sw**2))
1247         feqt4 = cos(2._dp*rml)*(5._dp*ec**2*ssw/4._dp)
1248         feqt5 = sin(3._dp*rml)*(eyt*cw)
1249         feqt6 = -cos(3._dp*rml)*(eyt*sw)
1250         feqt7 = -sin(4._dp*rml)*(.5_dp*yt**2)
1251         feqt = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7
1252         eqt = feqt*13751.0_dp
1253
1254 !-------------------------------------------------------------------
1255 ! convert eq of time from sec to deg
1256 !-------------------------------------------------------------------
1257         reqt = eqt/240._dp
1258 !-------------------------------------------------------------------
1259 ! calc right ascension in rads
1260 !-------------------------------------------------------------------
1261         ra = ml - reqt
1262         rra = ra*d2r
1263 !-------------------------------------------------------------------
1264 ! calc declination in rads, deg
1265 !-------------------------------------------------------------------
1266         tab = 0.43360_dp*sin(rra)
1267         rdecl = atan(tab)
1268         decl = rdecl/d2r
1269 !-------------------------------------------------------------------
1270 ! calc local hour angle
1271 !-------------------------------------------------------------------
1272         lbgmt = 12.0_dp - eqt/3600._dp + long*24._dp/360._dp
1273         lzgmt = 15.0_dp*(gmt-lbgmt)
1274         zpt = lzgmt*d2r
1275         csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
1276         if(csz.gt.1)print *,'calczen,csz ',csz
1277         csz = min( 1._dp,csz )
1278 !       zr = acos(csz)
1279 !       zenith = zr/d2r
1280         zr = acos(csz)
1281         zenith = zr/d2r
1282 !-------------------------------------------------------------------
1283 ! calc local solar azimuth
1284 !-------------------------------------------------------------------
1285         caz = (sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr))
1286         if(caz < -0.999999_dp)then
1287           azimuth=180._dp
1288         elseif(caz > 0.999999_dp)then
1289           azimuth=0._dp
1290         else
1291           raz = acos(caz)
1292           azimuth = raz/d2r
1293         endif
1294 !       caz = min(1.,(sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr)))
1295 !       if(caz.lt.-1)print *,calczen ,caz
1296 !       caz = max(-1.,caz)
1297 !       raz = acos(caz)
1298 !       azimuth = raz/d2r
1299
1300         if( lzgmt > 0._dp ) azimuth = azimuth + (2._dp*(180._dp - azimuth))
1301
1302       end subroutine calc_zenith
1303
1304       subroutine sundis( julday, esrm2 )
1305 !-----------------------------------------------------------------------------
1306 ! purpose:
1307 ! calculate earth-sun distance variation for a given date. based on
1308 ! fourier coefficients originally from: spencer, j.w., 1971, fourier
1309 ! series representation of the position of the sun, search, 2:172
1310 !-----------------------------------------------------------------------------
1311 ! parameters:
1312 ! idate - integer, specification of the date, from yymmdd (i)
1313 ! esrm2 - real(dp), variation of the earth-sun distance (o)
1314 ! esrm2 = (average e/s dist)^2 / (e/s dist on day idate)^2
1315 !-----------------------------------------------------------------------------
1316       implicit none
1317 !-----------------------------------------------------------------------------
1318 ! ... dummy arguments
1319 !-----------------------------------------------------------------------------
1320       integer, intent(in) :: julday
1321       real(dp), intent(out) :: esrm2
1322 !-----------------------------------------------------------------------------
1323 ! ... local variables
1324 !-----------------------------------------------------------------------------
1325       real(dp), parameter :: pi = 3.1415926_dp
1326       real(dp) :: dayn, thet0
1327       real(dp) :: sinth, costh, sin2th, cos2th
1328 !-----------------------------------------------------------------------------
1329 ! ... parse date to find day number (julian day)
1330 !-----------------------------------------------------------------------------
1331       dayn = real(julday - 1,kind=dp) + .5_dp
1332 !-----------------------------------------------------------------------------
1333 ! ... define angular day number and compute esrm2:
1334 !-----------------------------------------------------------------------------
1335       thet0 = 2._dp*pi*dayn/365._dp
1336 !-----------------------------------------------------------------------------
1337 ! ... calculate sin(2*thet0), cos(2*thet0) from
1338 ! addition theorems for trig functions for better
1339 ! performance; the computation of sin2th, cos2th
1340 ! is about 5-6 times faster than the evaluation
1341 ! of the intrinsic functions sin and cos
1342 !-----------------------------------------------------------------------------
1343       sinth  = sin( thet0 )
1344       costh  = cos( thet0 )
1345       sin2th = 2._dp*sinth*costh
1346       cos2th = costh*costh - sinth*sinth
1347       esrm2  = 1.000110_dp + .034221_dp*costh + .001280_dp*sinth + .000719_dp*cos2th + .000077_dp*sin2th
1348
1349       end subroutine sundis
1350
1351       subroutine spheres( nz, z, zen, dsdh, nid )
1352 !-----------------------------------------------------------------------------
1353 ! purpose:
1354 ! calculate slant path over vertical depth ds/dh in spherical geometry.
1355 ! calculation is based on: a.dahlback, and k.stamnes, a new spheric model
1356 ! for computing the radiation field available for photolysis and heating
1357 ! at twilight, planet.space sci., v39, n5, pp. 671-683, 1991 (appendix b)
1358 !-----------------------------------------------------------------------------
1359 ! parameters:
1360 ! nz - integer, number of specified altitude levels in the working (i)
1361 ! grid
1362 ! z - real(dp), specified altitude working grid (km) (i)
1363 ! zen - real(dp), solar zenith angle (degrees) (i)
1364 ! dsdh - real(dp), slant path of direct beam through each layer crossed (o)
1365 ! when travelling from the top of the atmosphere to layer i;
1366 ! dsdh(i,j), i = 0..nz-1, j = 1..nz-1
1367 ! nid - integer, number of layers crossed by the direct beam when (o)
1368 ! travelling from the top of the atmosphere to layer i;
1369 ! nid(i), i = 0..nz-1
1370 !-----------------------------------------------------------------------------
1371       implicit none
1372 !-----------------------------------------------------------------------------
1373 ! ... dummy arguments
1374 !-----------------------------------------------------------------------------
1375       integer, intent(in)  :: nz
1376       integer, intent(out) :: nid(0:nz-1)
1377       real(dp), intent(in) :: zen
1378       real(dp), intent(in) :: z(nz)
1379       real(dp), intent(out) :: dsdh(0:nz-1,nz-1)
1380 !-----------------------------------------------------------------------------
1381 ! ... local variables
1382 !-----------------------------------------------------------------------------
1383       real(dp), parameter :: radius = 6.371e+3_dp         ! radius earth (km)
1384       real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp
1385       integer :: i, j, k
1386       integer :: id
1387       real(dp) :: re, ze(nz)
1388       real(dp) :: zd(0:nz-1)
1389       real(dp) :: zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
1390
1391       zenrad = zen*d2r
1392 !-----------------------------------------------------------------------------
1393 ! ... include the elevation above sea level to the radius of the earth:
1394 !-----------------------------------------------------------------------------
1395       re = radius + z(1)
1396 !-----------------------------------------------------------------------------
1397 ! ... correspondingly z changed to the elevation above earth surface:
1398 !-----------------------------------------------------------------------------
1399       ze(1:nz) = z(1:nz) - z(1)
1400 !-----------------------------------------------------------------------------
1401 ! ... inverse coordinate of z
1402 !-----------------------------------------------------------------------------
1403       zd(0) = ze(nz)
1404       do k = 1,nz-1
1405         zd(k) = ze(nz - k)
1406       end do
1407 !-----------------------------------------------------------------------------
1408 ! ... initialize dsdh(i,j), nid(i)
1409 !-----------------------------------------------------------------------------
1410       dsdh(0:nz-1,1:nz-1) = 0._dp
1411       nid(0:nz-1) = 0
1412 !-----------------------------------------------------------------------------
1413 ! ... calculate ds/dh of every layer
1414 !-----------------------------------------------------------------------------
1415       do i = 0,nz-1
1416         rpsinz = (re + zd(i)) * sin(zenrad)
1417         if( zen > 90._dp .and. rpsinz < re ) then
1418            nid(i) = -1
1419         else
1420 !-----------------------------------------------------------------------------
1421 ! ... find index of layer in which the screening height lies
1422 !-----------------------------------------------------------------------------
1423            id = i
1424            if( zen > 90._dp ) then
1425               do j = 1,nz-1
1426                  if( (rpsinz < ( zd(j-1) + re ) ) .and. (rpsinz >= ( zd(j) + re )) ) then
1427                     id = j
1428                  end if
1429               end do
1430            end if
1431            do j = 1,id
1432              if( j == id .and. id == i .and. zen > 90._dp ) then
1433                 sm = -1._dp
1434              else
1435                 sm = 1._dp
1436              end if
1437              rj = re + zd(j-1)
1438              rjp1 = re + zd(j)
1439              dhj = zd(j-1) - zd(j)
1440              ga = max( 0._dp,rj*rj - rpsinz*rpsinz )
1441              gb = max( 0._dp,rjp1*rjp1 - rpsinz*rpsinz )
1442              if( id > i .and. j == id ) then
1443                 dsj = sqrt( ga )
1444              else
1445                 dsj = sqrt( ga ) - sm*sqrt( gb )
1446              end if
1447              dsdh(i,j) = dsj / dhj
1448            end do
1449            nid(i) = id
1450         end if
1451       end do
1452
1453       end subroutine spheres
1454
1455       subroutine setz( nz, cz, tlev, c, ndx, adjcoe )
1456
1457       use module_wave_data, only : tuv_jmax
1458 !-----------------------------------------------------------------------------
1459 ! adjcoe - real(dp), coross section adjust coefficients
1460 !-----------------------------------------------------------------------------
1461       implicit none
1462 !-----------------------------------------------------------------------------
1463 ! ... dummy arguments
1464 !-----------------------------------------------------------------------------
1465       integer, intent(in) :: nz
1466       integer, intent(in) :: ndx
1467       real(dp), intent(in) :: cz(nz)
1468       real(dp), intent(in) :: tlev(nz)
1469       real(dp), intent(in) :: c(5,tuv_jmax)
1470       real(dp), intent(inout) :: adjcoe(nz,tuv_jmax)
1471 !-----------------------------------------------------------------------------
1472 ! ... local variables
1473 !-----------------------------------------------------------------------------
1474       integer :: k, m
1475       real(dp) :: tt, adjin
1476       real(dp) :: c0, c1, c2
1477       real(dp) :: xz(nz)
1478 !-----------------------------------------------------------------------------
1479 ! 1 o2 + hv -> o + o
1480 ! 2 o3 -> o2 + o(1d)
1481 ! 3 o3 -> o2 + o(3p)
1482 ! 4 no2 -> no + o(3p)
1483 ! 5 no3 -> no + o2
1484 ! 6 no3 -> no2 + o(3p)
1485 ! 7 n2o5 -> no3 + no + o(3p)
1486 ! 8 n2o5 -> no3 + no2
1487 ! 9 n2o + hv -> n2 + o(1d)
1488 ! 10 ho2 + hv -> oh + o
1489 ! 11 h2o2 -> 2 oh
1490 ! 12 hno2 -> oh + no
1491 ! 13 hno3 -> oh + no2
1492 ! 14 hno4 -> ho2 + no2
1493 ! 15 ch2o -> h + hco
1494 ! 16 ch2o -> h2 + co
1495 ! 17 ch3cho -> ch3 + hco
1496 ! 18 ch3cho -> ch4 + co
1497 ! 19 ch3cho -> ch3co + h
1498 ! 20 c2h5cho -> c2h5 + hco
1499 ! 21 chocho -> products
1500 ! 22 ch3cocho -> products
1501 ! 23 ch3coch3
1502 ! 24 ch3ooh -> ch3o + oh
1503 ! 25 ch3ono2 -> ch3o+no2
1504 ! 26 pan + hv -> products
1505 ! 27 mvk + hv -> products
1506 ! 28 macr + hv -> products
1507 ! 29 hyac + hv -> products
1508 ! 30 glyald + hv -> products
1509 !-----------------------------------------------------------------------------
1510       xz(1:nz) = cz(1:nz)*1.e-18_dp
1511       do k = 1,tuv_jmax
1512          adjcoe(1:nz,k) = 1._dp
1513       end do
1514       tt = tlev(1)/281._dp
1515       do m = 1,tuv_jmax
1516          adjin = 1._dp
1517          if( m == 2 ) then
1518 !----------------------------------------------------------------------
1519 ! ... temperature modification
1520 ! t0.9 (1.05) t0.95(1.025) t1.0(1.0) t1.15(1.02) t1.1(1.04)
1521 !----------------------------------------------------------------------
1522             select case( ndx )
1523                case( 20 )
1524                   c0 = 4.52372_dp ; c1 = -5.94317_dp ; c2 = 2.63156_dp
1525                case( 40 )
1526                   c0 = 4.99378_dp ; c1 = -7.92752_dp ; c2 = 3.94715_dp
1527                case( 60 )
1528                   c0 = .969867_dp ; c1 = -.841035_dp ; c2 = .878835_dp
1529                case( 80 )
1530                   c0 = 1.07801_dp ; c1 = -2.39580_dp ; c2 = 2.32632_dp
1531             end select
1532             adjin = c0 + tt*(c1 + c2*tt)
1533          else if( m == 11 ) then
1534 !----------------------------------------------------------------------
1535 ! ... temperature modification
1536 ! t0.9 (1.05) t0.95(1.025) t1.0(1.0) t1.15(1.02) t1.1(1.04)
1537 !----------------------------------------------------------------------
1538             select case( ndx )
1539                case( 20 )
1540                   c0 = 2.43360_dp ; c1 = -3.61363_dp ; c2 = 2.19018_dp
1541                case( 40 )
1542                   c0 = 3.98265_dp ; c1 = -6.90516_dp ; c2 = 3.93602_dp
1543                case( 60 )
1544                   c0 = 3.49843_dp ; c1 = -5.98839_dp ; c2 = 3.50262_dp
1545                case( 80 )
1546                   c0 = 3.06312_dp ; c1 = -5.26281_dp ; c2 = 3.20980_dp
1547             end select
1548             adjin = c0 + tt*(c1 + c2*tt)
1549          end if
1550          call calcoe( nz, c(1,m), xz, tt, adjin, adjcoe(1,m) )
1551       end do
1552
1553       end subroutine setz
1554
1555       subroutine setozo( nz, z, tlay, dto3, &
1556                          to3, o3, airlev, o3top, o3toms )
1557
1558       use module_wave_data, only : nw, wl, xso3, s226, s263, s298
1559 !-----------------------------------------------------------------------------
1560 ! purpose:
1561 ! set up an altitude profile of ozone, and corresponding absorption
1562 ! optical depths. subroutine includes a shape-conserving scaling method
1563 ! that allows scaling of the entire profile to a given overhead ozone
1564 ! column amount.
1565 !-----------------------------------------------------------------------------
1566 ! parameters:
1567 ! dobnew - real(dp), overhead ozone column amount (du) to which profile (i)
1568 ! should be scaled. if dobnew < 0, no scaling is done
1569 ! nz - integer, number of specified altitude levels in the working (i)
1570 ! grid
1571 ! z - real(dp), specified altitude working grid (km)
1572 ! nw - integer, number of specified intervals + 1 in working (i)
1573 ! wavelength grid
1574 ! wl - real(dp), vector of lower limits of wavelength intervals in (i)
1575 ! working wavelength grid
1576 ! xso3 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
1577 ! each specified wavelength (wmo value at 273)
1578 ! s226 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
1579 ! each specified wavelength (value from molina and molina at 226k)
1580 ! s263 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
1581 ! each specified wavelength (value from molina and molina at 263k)
1582 ! s298 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i)
1583 ! each specified wavelength (value from molina and molina at 298k)
1584 ! tlay - real(dp), temperature (k) at each specified altitude layer (i)
1585 ! dto3 - real(dp), optical depth due to ozone absorption at each (o)
1586 ! specified altitude at each specified wavelength
1587 ! to3 - real(dp), totol ozone column (o)
1588 !-----------------------------------------------------------------------------
1589
1590       implicit none
1591
1592 !-----------------------------------------------------------------------------
1593 ! ... dummy arguments
1594 !-----------------------------------------------------------------------------
1595       integer, intent(in) :: nz
1596       real(dp), intent(in) :: o3top
1597       real(dp), intent(in) :: o3toms
1598       real(dp), intent(in) :: z(nz)
1599       real(dp), intent(in) :: tlay(nz-1)
1600       real(dp), intent(in) :: airlev(nz)
1601       real(dp), intent(out) :: dto3(nz-1,nw-1)
1602       real(dp), intent(out) :: to3(nz)
1603       real(dp), intent(in) :: o3(nz)
1604 !-----------------------------------------------------------------------------
1605 ! ... local variables
1606 !-----------------------------------------------------------------------------
1607       integer  :: k, iw, nd, ilat
1608       real(dp) :: so3
1609       real(dp) :: scale
1610       real(dp) :: div1, div2
1611       real(dp) :: o3den(nz)
1612       real(dp) :: cz(nz)
1613
1614 !-----------------------------------------------------------------------------
1615 ! ... compute column increments
1616 !-----------------------------------------------------------------------------
1617 !     o3(1:nz) = o3(1:nz)*airlev(1:nz)
1618 !     cz(1:nz-1) = (o3(2:nz)) * 1.e5 * (z(2:nz) - z(1:nz-1))
1619
1620       o3den(1:nz) = o3(1:nz)*airlev(1:nz)
1621       cz(1:nz-1)  = 0.5_dp*(o3den(2:nz) + o3den(1:nz-1))*km2cm*(z(2:nz) - z(1:nz-1))
1622
1623       to3(nz)  = o3top
1624       do k = nz-1,1,-1
1625          to3(k) = to3(k+1) + cz(k)
1626       end do
1627
1628 ! Do not double count o3top
1629       cz(nz-1)    = cz(nz-1) + o3top
1630 !-----------------------------------------------------------------------------
1631 ! ... scale o3 using toms data
1632 !-----------------------------------------------------------------------------
1633       if( o3toms > 0.0_dp ) then
1634         scale = o3toms/(to3(1)/2.687e16_dp)
1635         cz(1:nz)  = cz(1:nz)*scale
1636         to3(1:nz) = to3(1:nz)*scale
1637       endif
1638 !-----------------------------------------------------------------------------
1639 ! ... calculate ozone optical depth for each layer, with temperature
1640 ! correction. output, dto3(kz,kw)
1641 !-----------------------------------------------------------------------------
1642       div1 = 1._dp / ( 263._dp - 226._dp)
1643       div2 = 1._dp / ( 298._dp - 263._dp)
1644       do iw = 1,nw-1
1645          so3 = xso3(iw)
1646          do k = 1,nz-1
1647             if( wl(iw) > 240.5_dp .and. wl(iw+1) < 350._dp ) then
1648                if( tlay(k) < 263._dp ) then
1649                   so3 = s226(iw) + (s263(iw) - s226(iw)) * div1 * (tlay(k) - 226._dp)
1650                else
1651                   so3 = s263(iw) + (s298(iw) - s263(iw)) * div2 * (tlay(k) - 263._dp)
1652                end if
1653             end if
1654             dto3(k,iw) = cz(k)*so3
1655          end do
1656       end do
1657
1658       end subroutine setozo
1659
1660       subroutine set_o2_xsect( nz, z, cz, vcol, scol, dto2, xso2 )
1661
1662       use module_wave_data, only : nw, wl
1663 !-----------------------------------------------------------------------------
1664 ! purpose:
1665 ! compute equivalent optical depths for o2 absorption, parameterized in
1666 ! the sr bands and the lyman-alpha line.
1667 !-----------------------------------------------------------------------------
1668 ! parameters:
1669 ! nz - integer, number of specified altitude levels in the working (i)
1670 ! grid
1671 ! z - real(dp), specified altitude working grid (km)
1672 ! nw - integer, number of specified intervals + 1 in working
1673 ! wavelength grid
1674 ! wl - real(dp), vector of lower limits of wavelength intervals in
1675 ! working wavelength grid
1676 ! cz - real(dp), number of air molecules per cm^2 at each specified
1677 ! altitude layer
1678 ! zen - real(dp), solar zenith angle
1679 ! dto2 - real(dp), optical depth due to o2 absorption at each specified (o)
1680 ! vertical layer at each specified wavelength
1681 ! xso2 - real(dp), molecular absorption cross section in sr bands at (o)
1682 ! each specified altitude and wavelength. includes herzberg
1683 ! continuum.
1684 !-----------------------------------------------------------------------------
1685 ! edit history:
1686 ! 02/98 included lyman-alpha parameterization
1687 ! 03/97 fix dto2 problem at top level (nz)
1688 ! 02/97 changed offset for grid-end interpolation to relative number
1689 ! (x * (1 +- deltax))
1690 ! 08/96 modified for early exit, no redundant read of data and smaller
1691 ! internal grid if possible; internal grid uses user grid points
1692 ! whenever possible
1693 ! 07/96 modified to work on internal grid and interpolate final values
1694 ! onto the user-defined grid
1695 !-----------------------------------------------------------------------------
1696       implicit none
1697 !-----------------------------------------------------------------------------
1698 ! ... dummy arguments
1699 !-----------------------------------------------------------------------------
1700       integer, intent(in) :: nz
1701       real(dp), intent(in) :: cz(nz)
1702       real(dp), intent(in) :: z(nz)
1703       real(dp), intent(in) :: vcol(nz)
1704       real(dp), intent(in) :: scol(nz)
1705       real(dp), intent(out) :: dto2(nz-1,nw-1)
1706       real(dp), intent(out) :: xso2(nw-1,nz)
1707 !-----------------------------------------------------------------------------
1708 ! ... local variables
1709 !-----------------------------------------------------------------------------
1710       integer :: iw, iz, igast
1711       real(dp) :: secchi(nz)
1712 !-----------------------------------------------------------------------------
1713 ! ... o2 optical depth and equivalent cross section on kockarts grid
1714 !-----------------------------------------------------------------------------
1715       real(dp) :: dto2k(nz-1,ngast-1)
1716       real(dp) :: xso2k(nz,ngast-1)
1717 !-----------------------------------------------------------------------------
1718 ! ... o2 optical depth and equivalent cross section in the lyman-alpha region
1719 !-----------------------------------------------------------------------------
1720       real(dp) :: dto2la(nz-1,nla-1)
1721       real(dp) :: xso2la(nz,nla-1)
1722 !-----------------------------------------------------------------------------
1723 ! ... temporary one-dimensional storage for optical depth and cross section values
1724 ! xxtmp - on internal grid
1725 ! xxuser - on user defined grid
1726 !-----------------------------------------------------------------------------
1727       real(dp), dimension(2*nwint) :: dttmp, xstmp    ! CHANGED by Guohui
1728       real(dp), dimension(nwint) :: dtuser, xsuser    ! CHANGED by Guohui
1729       real(dp) :: o2col(nz)
1730       real(dp) :: x, y
1731       real(dp) :: delo2
1732 !-----------------------------------------------------------------------------
1733 ! ... check, whether user grid is in the o2 absorption band at all...
1734 ! if not, set cross section and optical depth values to zero and return
1735 !-----------------------------------------------------------------------------
1736       dto2(:nz-1,:nw-1) = 0._dp
1737       xso2(:nw-1,:nz) = 0._dp
1738       if( wl(1) > 243._dp ) then
1739          return
1740       end if
1741 !-----------------------------------------------------------------------------
1742 ! ... sec xhi or chapman calculation
1743 ! for zen > 95 degrees, use zen = 95. (this is only to compute effective o2
1744 ! cross sections. still, better than setting dto2 = 0. as was done up to
1745 ! version 4.0) sm 1/2000
1746 ! in future could replace with mu2(iz) (but mu2 is also wavelength-depenedent)
1747 ! or imporved chapman function
1748 !-----------------------------------------------------------------------------
1749 !-----------------------------------------------------------------------------
1750 ! ... slant o2 column
1751 !-----------------------------------------------------------------------------
1752       o2col(1:nz) = 0.2095_dp * scol(1:nz)
1753 !-----------------------------------------------------------------------------
1754 ! ... effective secant of solar zenith angle. use 2.0 if no direct sun.
1755 ! for nz, use value at nz-1
1756 !-----------------------------------------------------------------------------
1757       secchi(1:nz-1) = scol(1:nz-1)/vcol(1:nz-1)
1758       where( secchi(1:nz-1) == 0._dp )
1759          secchi(1:nz-1) = 2._dp
1760       endwhere
1761       secchi(nz) = secchi(nz-1)
1762 !-----------------------------------------------------------------------------
1763 ! ... if necessary:
1764 ! kockarts parameterization of the sr bands, output values of o2
1765 ! optical depth and o2 equivalent cross section are on his grid
1766 !-----------------------------------------------------------------------------
1767       if( wl(1) < wlgast(ngast) .and. wl(nw) > wlgast(1) ) then
1768            call schu( nz, o2col, secchi, dto2k, xso2k )
1769       else
1770          dto2k(:,:) = 0._dp
1771          xso2k(:,:) = 0._dp
1772       end if
1773 !-----------------------------------------------------------------------------
1774 ! ... lyman-alpha parameterization, output values of o2 opticaldepth
1775 ! and o2 effective (equivalent) cross section
1776 !-----------------------------------------------------------------------------
1777       if( wl(1) <= wlla(nla) .and. wl(nw) >= wlla(1) ) then
1778          call lymana( nz, o2col, secchi, dto2la, xso2la )
1779       else
1780          dto2la(:,:) = 0._dp
1781          xso2la(:,:) = 0._dp
1782       end if
1783 !-----------------------------------------------------------------------------
1784 ! ... loop through the altitude levels
1785 !-----------------------------------------------------------------------------
1786       do iz = 1,nz
1787          igast = 0
1788 !-----------------------------------------------------------------------------
1789 ! ... loop through the internal wavelength grid
1790 !-----------------------------------------------------------------------------
1791          do iw = 1,nwint-1
1792 !-----------------------------------------------------------------------------
1793 ! ... if outside kockarts grid and outside lyman-alpha, use the
1794 ! jpl/brasseur+solomon data, if inside
1795 ! kockarts grid, use the parameterized values from the call to schu,
1796 ! if inside lyman-alpha, use the paraemterized values from call to lymana
1797 !-----------------------------------------------------------------------------
1798             if( wlint(iw+1) <= wlgast(1) .or. wlint(iw) >= wlgast(ngast) ) then
1799               if( wlint(iw+1) <= wlla(1) .or. wlint(iw) >= wlla(nla) ) then
1800                  xstmp(iw) = xso2int(iw)
1801               else
1802                  xstmp(iw) = xso2la(iz,1)
1803               end if
1804             else
1805                igast = igast + 1
1806                xstmp(iw) = xso2k(iz,igast)
1807             end if
1808 !-----------------------------------------------------------------------------
1809 ! ... compute the area in each bin (for correct interpolation purposes only!)
1810 !-----------------------------------------------------------------------------
1811             xstmp(iw) = xstmp(iw) * (wlint(iw+1) - wlint(iw))
1812          end do
1813 !-----------------------------------------------------------------------------
1814 ! ... interpolate o2 cross section from the internal grid onto the user grid
1815 !-----------------------------------------------------------------------------
1816          call inter3( nw, wl, xsuser, nwint, wlint, xstmp )
1817          xso2(1:nw-1,iz) = xsuser(1:nw-1)/(wl(2:nw) - wl(1:nw-1))
1818       end do
1819       do iz = 1,nz-1
1820          igast = 0
1821          delo2 = .2095_dp * cz(iz) ! vertical o2 column
1822 !-----------------------------------------------------------------------------
1823 ! ... loop through the internal wavelength grid
1824 !-----------------------------------------------------------------------------
1825          do iw = 1,nwint-1
1826 !-----------------------------------------------------------------------------
1827 ! ... if outside kockarts grid and outside lyman-alpha, use the
1828 ! jpl/brasseur+solomon data, if inside
1829 ! kockarts grid, use the parameterized values from the call to schu,
1830 ! if inside lyman-alpha, use the paraemterized values from call to lymana
1831 !-----------------------------------------------------------------------------
1832             if( wlint(iw+1) <= wlgast(1) .or. wlint(iw) >= wlgast(ngast) ) then
1833               if( wlint(iw+1) <= wlla(1) .or. wlint(iw) >= wlla(nla) ) then
1834                  dttmp(iw) = xso2int(iw) * delo2
1835               else
1836                  dttmp(iw) = dto2la(iz,1)
1837               end if
1838             else
1839                igast = igast + 1
1840                dttmp(iw) = dto2k(iz,igast)
1841             end if
1842 !-----------------------------------------------------------------------------
1843 ! ... compute the area in each bin (for correct interpolation purposes only!)
1844 !-----------------------------------------------------------------------------
1845             dttmp(iw) = dttmp(iw) * (wlint(iw+1) - wlint(iw))
1846          end do
1847 !-----------------------------------------------------------------------------
1848 ! ... interpolate o2 optical depth from the internal grid onto the user grid
1849 !-----------------------------------------------------------------------------
1850          call inter3( nw, wl, dtuser, nwint, wlint, dttmp )
1851          dto2(iz,1:nw-1) = dtuser(1:nw-1)/(wl(2:nw) - wl(1:nw-1))
1852       end do
1853
1854       end subroutine set_o2_xsect
1855
1856       subroutine setcld( nz, z, xlwc, dtcld, omcld, gcld )
1857
1858       use module_wave_data, only : nw
1859 !-----------------------------------------------------------------------------
1860 != PURPOSE:
1861 != Set up cloud optical depth, single albedo and g
1862 !-----------------------------------------------------------------------------
1863 != PARAMETERS:
1864 != PARAMETERS:
1865 != NZ - INTEGER, number of specified altitude levels in the working (I)
1866 != grid
1867 != Z - real(dp), specified altitude working grid (km) (I)
1868 != XLWC Cloud water content g/M3 (I)
1869 !=
1870 != dtcld - cloud optical depth
1871 != omcld - cloud droplet single albedo
1872 != gcld  - g
1873 !-----------------------------------------------------------------------------
1874 !
1875 ! VERTICAL DOMAIN is from bottom(1) to TOP (TOP=nz)
1876 ! CCM from top(1) to bottom(nz)
1877 !-----------------------------------------------------------------------------
1878 !****************************************************************************
1879       implicit none
1880 !-----------------------------------------------------------------------------
1881 ! ... Dummy arguments
1882 !-----------------------------------------------------------------------------
1883       integer, intent(in) :: nz
1884       real(dp), intent(in) :: z(nz)
1885       real(dp), intent(in) :: xlwc(nz)
1886       real(dp), intent(out) :: dtcld(nz-1,nw-1)
1887       real(dp), intent(out) :: omcld(nz-1,nw-1)
1888       real(dp), intent(out) :: gcld(nz-1,nw-1)
1889 !-----------------------------------------------------------------------------
1890 ! ... Local variables
1891 !-----------------------------------------------------------------------------
1892       real(dp), parameter :: wden = 1.0_dp * 1.e6_dp ! g/m3 (1 m3 water = 1e6 g water)
1893       real(dp), parameter :: re = 10.0_dp * 1.e-6_dp ! assuming cloud drop radius = 10 um to M
1894
1895       integer  :: i                                  ! vertical index for each blk
1896       real(dp) :: dz(nz-1)
1897 !-----------------------------------------------------------------------------
1898 ! ... calculate vertical interveral dz
1899 !-----------------------------------------------------------------------------
1900       do i=1,nz-1
1901           dz(i) = (z(i+1)-z(i)) * 1.e3_dp ! dz in Meters
1902       end do
1903 !----------------------------------------------------
1904 ! ....calculate cloud optical depth T
1905 ! following Liao et al. JGR, 104, 23697, 1999
1906 !----------------------------------------------------
1907       do i = 1,nz-1
1908         dtcld(i,:) = 1.5_dp * xlwc(i)*dz(i)/ (wden * re)
1909         dtcld(i,:) = max( dtcld(i,:), 0._dp )
1910         omcld(i,:) = .9999_dp
1911         gcld (i,:) = .85_dp
1912       end do
1913
1914       end subroutine setcld
1915
1916       subroutine schu_inti
1917 !-----------------------------------------------------------------------------
1918 ! ... initialize the tables
1919 !-----------------------------------------------------------------------------
1920       implicit none
1921 !-----------------------------------------------------------------------------
1922 ! ... local variables
1923 !-----------------------------------------------------------------------------
1924       integer :: i, iw, k, j1, jp1, j
1925       real(dp) :: col
1926       real(dp), dimension(6) :: a0, a1, b0, b1
1927       do iw = 1,ngast-1
1928          x_table(0,iw) = sum( a_schu(1:11:2,iw) )
1929          d_table(0,iw) = sum( b_schu(1:11:2,iw) )
1930          do k = 1,tdim
1931             col = 22._dp + t_del*real(k-1,kind=dp)
1932             o2_table(k) = col
1933             col = 10._dp**col
1934             a1(:) = a_schu(2:12:2,iw)*col
1935             b1(:) = b_schu(2:12:2,iw)*col
1936             where( a1(:) < 500._dp )
1937                a0(:) = exp( -a1(:) )
1938             elsewhere
1939                a0(:) = 0._dp
1940             endwhere
1941             where( b1(:) < 500._dp )
1942                b0(:) = exp( -b1(:) )
1943             elsewhere
1944                b0(:) = 0._dp
1945             endwhere
1946             x_table(k,iw) = dot_product( a_schu(1:11:2,iw),a0(:) )
1947             d_table(k,iw) = dot_product( b_schu(1:11:2,iw),b0(:) )
1948          end do
1949       end do
1950
1951       end subroutine schu_inti
1952
1953       subroutine schu( nz, o2col, secchi, dto2, xscho2 )
1954 !-----------------------------------------------------------------------------
1955 ! purpose:
1956 ! calculate the equivalent absorption cross section of o2 in the sr bands.
1957 ! the algorithm is based on: g.kockarts, penetration of solar radiation
1958 ! in the schumann-runge bands of molecular oxygen: a robust approximation,
1959 ! annales geophysicae, v12, n12, pp. 1207ff, dec 1994. calculation is
1960 ! done on the wavelength grid used by kockarts (1994). final values do
1961 ! include effects from the herzberg continuum.
1962 !-----------------------------------------------------------------------------
1963 ! parameters:
1964 ! nz - integer, number of specified altitude levels in the working (i)
1965 ! grid
1966 ! o2col - real(dp), slant overhead o2 column (molec/cc) at each specified (i)
1967 ! altitude
1968 ! dto2 - real(dp), optical depth due to o2 absorption at each specified (o)
1969 ! vertical layer at each specified wavelength
1970 ! xscho2 - real(dp), molecular absorption cross section in sr bands at (o)
1971 ! each specified wavelength. includes herzberg continuum
1972 !-----------------------------------------------------------------------------
1973       implicit none
1974 !-----------------------------------------------------------------------------
1975 ! ... dummy arguments
1976 !-----------------------------------------------------------------------------
1977       integer, intent(in) :: nz
1978       real(dp), intent(inout) :: dto2(nz-1,ngast-1)
1979       real(dp), intent(inout) :: xscho2(nz,ngast-1)
1980       real(dp), intent(in) :: o2col(nz)
1981       real(dp), intent(in) :: secchi(nz)
1982 !-----------------------------------------------------------------------------
1983 ! ... local variables
1984 !-----------------------------------------------------------------------------
1985       integer :: i, iw, j, j1, jp1, k, ki, kp1
1986       integer :: index(nz)
1987       integer :: minind(1), maxind(1)
1988       real(dp) :: a0, a1, b0, b1
1989       real(dp), dimension(6) :: ac, bc
1990       real(dp), dimension(nz,6) :: aa, bb
1991       real(dp), dimension(nz) :: rjm, rjo2
1992       real(dp), dimension(nz) :: rjmi, rjo2i, lo2col, dels
1993 !-----------------------------------------------------------------------------
1994 ! ... initialize r(m)
1995 !-----------------------------------------------------------------------------
1996       rjm(1:nz) = 0._dp
1997       rjo2(1:nz) = 0._dp
1998 !-----------------------------------------------------------------------------
1999 ! ... initialize the table interpolation
2000 !-----------------------------------------------------------------------------
2001       where( o2col(:) /= 0 )
2002          lo2col(:) = log10( o2col(:) )
2003       endwhere
2004       do ki = 1,nz
2005          if( o2col(ki) /= 0._dp ) then
2006             if( lo2col(ki) <= o2_table(1) ) then
2007                dels(ki) = 0._dp
2008                index(ki) = 1
2009             else if( lo2col(ki) >= o2_table(tdim) ) then
2010                dels(ki) = 1._dp
2011                index(ki) = tdim-1
2012             else
2013                do k = 2,tdim
2014                   if( lo2col(ki) <= o2_table(k) ) then
2015                      dels(ki) = t_fac*(lo2col(ki) - o2_table(k-1))
2016                      index(ki) = k-1
2017                      exit
2018                   end if
2019                end do
2020             end if
2021          else
2022             index(ki) = 0
2023             dels(ki) = 0._dp
2024          end if
2025       end do
2026 !-----------------------------------------------------------------------------
2027 ! ... calculate sum of exponentials (eqs 7 and 8 of kockarts 1994)
2028 !-----------------------------------------------------------------------------
2029       do iw = 1,ngast-1
2030          do k = 1,nz
2031             ki = index(k)
2032             rjm(k) = x_table(ki,iw) + dels(k)*(x_table(ki+1,iw) - x_table(ki,iw))
2033             rjo2(k) = d_table(ki,iw) + dels(k)*(d_table(ki+1,iw) - d_table(ki,iw))
2034          end do
2035          do k = 1,nz-1
2036             if( rjm(k) > 1.e-100_dp ) then
2037                kp1 = k + 1
2038                if( rjm(kp1) > 0._dp ) then
2039                   dto2(k,iw) = log( rjm(kp1) ) / secchi(kp1) - log( rjm(k) ) * secchi(k)
2040                else
2041                   dto2(k,iw) = 1000._dp
2042                end if
2043             else
2044                dto2(k,iw) = 1000._dp
2045             end if
2046          end do
2047          do k = 1,nz
2048             if( rjm(k) > 1.e-100_dp ) then
2049                if( rjo2(k) > 1.e-100_dp ) then
2050                   xscho2(k,iw) = rjo2(k)/rjm(k)
2051                else
2052                   xscho2(k,iw) = 0._dp
2053                end if
2054             else
2055                xscho2(k,iw) = 0._dp
2056             end if
2057          end do
2058       end do
2059
2060       end subroutine schu
2061
2062       subroutine rtlink( nz,                  &
2063                          nw,                  &
2064                          zen,                 &
2065                          albedo,              &
2066                          z,                   &
2067                          nid,                 &
2068                          dsdh,                &
2069                          dtrl,                &
2070                          dto3,                &
2071                          dto2,                &
2072                          dtcld, omcld, gcld,  &
2073 !                         dtcb1, omcb1, gcb1,  &
2074 !                         dtcb2, omcb2, gcb2,  &
2075 !                         dtoc1, omoc1, goc1,  &
2076 !                         dtoc2, omoc2, goc2,  &
2077 !                         dtant, omant, gant,  &
2078 !                         dtso4, omso4, gso4,  &
2079 !                         dtsal, omsal, gsal,  &
2080 ! rajesh: pass aerosol optical properties
2081                          dtaer, omaer, gaer,  &
2082                          radfld )
2083
2084       implicit none
2085 !-----------------------------------------------------------------------
2086 ! ... dummy arguments
2087 !-----------------------------------------------------------------------
2088       integer, intent(in) :: nz
2089       integer, intent(in) :: nw
2090       real(dp), intent(in)    :: z(nz)
2091       real(dp), intent(in)    :: zen
2092       real(dp), intent(in)    :: albedo(nw-1)
2093       integer, intent(in)    :: nid(0:nz-1)
2094       real(dp), intent(in)    :: dsdh(0:nz-1,nz-1)
2095       real(dp), intent(in)    :: dtrl(nz-1,nw-1)
2096       real(dp), intent(in)    :: dto3(nz-1,nw-1)
2097       real(dp), intent(in)    :: dto2(nz-1,nw-1)
2098       real(dp), intent(in)    :: dtcld(nz-1,nw-1), omcld(nz-1,nw-1), gcld(nz-1,nw-1)
2099 !      real(dp), intent(in)    :: dtcb1(nz-1,nw-1), omcb1(nz-1,nw-1), gcb1(nz-1,nw-1)
2100 !      real(dp), intent(in)    :: dtcb2(nz-1,nw-1), omcb2(nz-1,nw-1), gcb2(nz-1,nw-1)
2101 !      real(dp), intent(in)    :: dtoc1(nz-1,nw-1), omoc1(nz-1,nw-1), goc1(nz-1,nw-1)
2102 !      real(dp), intent(in)    :: dtoc2(nz-1,nw-1), omoc2(nz-1,nw-1), goc2(nz-1,nw-1)
2103 !      real(dp), intent(in)    :: dtant(nz-1,nw-1), omant(nz-1,nw-1), gant(nz-1,nw-1)
2104 !      real(dp), intent(in)    :: dtso4(nz-1,nw-1), omso4(nz-1,nw-1), gso4(nz-1,nw-1)
2105 !      real(dp), intent(in)    :: dtsal(nz-1,nw-1), omsal(nz-1,nw-1), gsal(nz-1,nw-1)
2106 ! rajesh: declare dust optical properties
2107       real(dp), intent(in)    :: dtaer(nz-1,nw-1), omaer(nz-1,nw-1), gaer(nz-1,nw-1)
2108       real(dp), intent(out)   :: radfld(nz,nw-1)
2109 !-----------------------------------------------------------------------
2110 ! ... local variables
2111 !-----------------------------------------------------------------------
2112       real(dp), parameter :: smallest = tiny( 1.0 )
2113       integer :: i, ii, iw
2114       real(dp) :: dtsct, dtabs
2115       real(dp) :: dscld, dacld
2116 !      real(dp) :: dscb1, dacb1
2117 !      real(dp) :: dscb2, dacb2
2118 !      real(dp) :: dsoc1, daoc1
2119 !      real(dp) :: dsoc2, daoc2
2120 !      real(dp) :: dsant, daant
2121 !      real(dp) :: dsso4, daso4
2122 !      real(dp) :: dssal, dasal
2123 ! rajesh: declare scattering and absorbing optical depths
2124       real(dp) :: dsaer, daaer
2125       real(dp), dimension(nz-1,nw-1) :: dt, om, g
2126 !-----------------------------------------------------------------------
2127 ! ... set any coefficients specific to rt scheme
2128 !-----------------------------------------------------------------------
2129       do iw = 1,nw-1
2130          do i = 1,nz-1
2131
2132             dscld = dtcld(i,iw)*omcld(i,iw)
2133             dacld = dtcld(i,iw)*(1._dp - omcld(i,iw))
2134
2135 !            dscb1 = dtcb1(i,iw)*omcb1(i,iw)
2136 !            dacb1 = dtcb1(i,iw)*(1._dp - omcb1(i,iw))
2137
2138 !            dscb2 = dtcb2(i,iw)*omcb2(i,iw)
2139 !            dacb2 = dtcb2(i,iw)*(1._dp - omcb2(i,iw))
2140
2141 !            dsoc1 = dtoc1(i,iw)*omoc1(i,iw)
2142 !            daoc1 = dtoc1(i,iw)*(1._dp - omoc1(i,iw))
2143
2144 !            dsoc2 = dtoc2(i,iw)*omoc2(i,iw)
2145 !            daoc2 = dtoc2(i,iw)*(1._dp - omoc2(i,iw))
2146
2147 !            dsant = dtant(i,iw)*omant(i,iw)
2148 !            daant = dtant(i,iw)*(1._dp - omant(i,iw))
2149
2150 !            dsso4 = dtso4(i,iw)*omso4(i,iw)
2151 !            daso4 = dtso4(i,iw)*(1._dp - omso4(i,iw))
2152
2153 !            dssal = dtsal(i,iw)*omsal(i,iw)
2154 !            dasal = dtsal(i,iw)*(1._dp - omsal(i,iw))
2155
2156 ! rajesh: determine scattering and absorption optical depths for dust
2157             dsaer = dtaer(i,iw)*omaer(i,iw)
2158             daaer = dtaer(i,iw)*(1._dp - omaer(i,iw))
2159             dtsct = dtrl(i,iw) + dscld + dsaer
2160 !                    dscb1 + dscb2 + dsoc1 + dsoc2 + dsant + dsso4 + dssal
2161             dtabs = dto3(i,iw) + dto2(i,iw) + dacld + daaer
2162 !                    dacb1 + dacb2 + daoc1 + daoc2 + daant + daso4 + dasal
2163
2164             dtabs = max( dtabs,smallest )
2165             dtsct = max( dtsct,smallest )
2166 !-----------------------------------------------------------------------
2167 ! ... invert z-coordinate
2168 !-----------------------------------------------------------------------
2169             ii = nz - i
2170
2171             dt(ii,iw) = dtsct + dtabs
2172             if( dtsct /= smallest ) then
2173                om(ii,iw) = dtsct/(dtsct + dtabs)
2174                g(ii,iw) = ( gcld(i,iw)*dscld + &
2175 !                            gcb1(i,iw)*dscb1 + &
2176 !                            gcb2(i,iw)*dscb2 + &
2177 !                            goc1(i,iw)*dsoc1 + &
2178 !                            goc2(i,iw)*dsoc2 + &
2179 !                            gant(i,iw)*dsant + &
2180 !                            gso4(i,iw)*dsso4 + &
2181 ! rajesh: add aerosol contribution
2182                             gaer(i,iw)*dsaer ) / dtsct
2183                g(ii,iw) = max( smallest, g(ii,iw) )
2184                g(ii,iw) = min( 1.d0, g(ii,iw) )
2185             else
2186                om(ii,iw) = smallest
2187                g(ii,iw)  = smallest
2188             end if
2189 !            print *, dt(ii,iw), om(ii,iw), g(ii,iw)
2190          end do
2191       end do
2192
2193 !-----------------------------------------------------------------------
2194 ! ... call rt routine
2195 !-----------------------------------------------------------------------
2196       call ps2str( nz, nw, zen, albedo, dt, om, &
2197                    g, dsdh, nid, radfld )
2198
2199       end subroutine rtlink
2200
2201       subroutine tridec( syscnt, order, lower, main, upper )
2202 !--------------------------------------------------------------------
2203 ! dimension of lower(syscnt,order), main(syscnt,order), upper(syscnt,order)
2204 ! arguments
2205 !
2206 ! latest revision june 1995
2207 !
2208 ! purpose trdi computes the solution of the tridiagonal
2209 ! linear system,
2210 ! b(1)*x(1)+c(1)*x(2) = y(1)
2211 ! a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1) = y(i)
2212 ! i=2,3,...,n-1
2213 ! a(n)*x(n-1)+b(n)*x(n) = y(n)
2214 !
2215 ! usage call trdi (n, a, b, c, x )
2216 !
2217 ! arguments
2218 !
2219 ! on input n
2220 ! the number of unknowns. this subroutine
2221 ! requires that n be greater than 2.
2222 !
2223 ! a
2224 ! the subdiagonal of the matrix is stored in
2225 ! locations a(2) through a(n).
2226 !
2227 ! b
2228 ! the diagonal of the matrix is stored in
2229 ! locations b(1) through b(n).
2230 !
2231 ! c
2232 ! the super-diagonal of the matrix is stored in
2233 ! locations c(1) through c(n-1).
2234 !
2235 ! x
2236 ! the right-hand side of the equations is
2237 ! stored in y(1) through y(n).
2238 !
2239 ! on output x
2240 ! an array which contains the solution to the
2241 ! system of equations.
2242 !
2243 ! special conditions this subroutine executes satisfactorily
2244 ! if the input matrix is diagonally dominant
2245 ! and non-singular. the diagonal elements
2246 ! are used to pivot, and no tests are made to
2247 ! determine singularity. if a singular, or
2248 ! nearly singular, matrix is used as input,
2249 ! a divide by zero or floating point overflow
2250 ! may result.
2251 !
2252 ! precision compiler dependent
2253 !
2254 ! language fortran
2255 !
2256 ! history originally written by nancy werner at ncar
2257 ! in october, 1973. modified by s. walters at
2258 ! ncar in june 1989 to functionally replace
2259 ! the cal routine tridla.
2260 !
2261 ! portability fortran 90
2262 !
2263 ! algorithm an lu-decomposition is obtained using the
2264 ! algorithm described in the reference below.
2265 !
2266 ! to avoid unnecessary divisions, the alpha
2267 ! values used in the routine are the
2268 ! reciprocals of the alpha values described
2269 ! in the reference below.
2270 !
2271 ! accuracy every component of the residual of the linear
2272 ! system (i.e. the difference between y and
2273 ! the matrix applied to x) should be less in
2274 ! magnitude than ten times the machine precision
2275 ! times the matrix order times the maximum
2276 ! absolute component of the solution vector
2277 ! times the largest absolute row sum of the
2278 ! input matrix.
2279 !
2280 ! timing the timing is roughly proportional to the
2281 ! order n of the linear system.
2282 !
2283 ! references analysis of numerical methods by
2284 ! e. isaacson and h. keller
2285 ! (john wiley and sons, 1966) pp. 55-58.
2286 !--------------------------------------------------------------------
2287       implicit none
2288 !--------------------------------------------------------------------
2289 ! ... dummy args
2290 !--------------------------------------------------------------------
2291       integer, intent(in) :: syscnt, order
2292       real(dp), intent(in), dimension(syscnt,order) :: lower
2293       real(dp), intent(inout), dimension(syscnt,order) :: main, upper
2294 !--------------------------------------------------------------------
2295 ! ... local variables
2296 !--------------------------------------------------------------------
2297       integer :: i
2298 !----------------------------------------------------------------------
2299 ! ... lu-decomposition
2300 !----------------------------------------------------------------------
2301       main(:,1) = 1._dp / main(:,1)
2302       upper(:,1) = upper(:,1)*main(:,1)
2303       do i = 2,order-1
2304          main(:,i) = 1._dp / (main(:,i) - lower(:,i)*upper(:,i-1))
2305          upper(:,i) = upper(:,i)*main(:,i)
2306       end do
2307
2308       end subroutine tridec
2309
2310       subroutine trislv( syscnt, order, lower, main, upper, x )
2311
2312       implicit none
2313 !----------------------------------------------------------------------
2314 ! ... dummy args
2315 !----------------------------------------------------------------------
2316       integer, intent(in) :: syscnt, order
2317       real(dp), intent(in), dimension(syscnt,order) :: lower, &
2318                                                        main, &
2319                                                        upper
2320       real(dp), intent(inout), dimension(syscnt,order) :: x
2321 !----------------------------------------------------------------------
2322 ! ... local variables
2323 !----------------------------------------------------------------------
2324       integer :: i, im1, j, n, nm1
2325       nm1 = order - 1
2326       n = order
2327 !----------------------------------------------------------------------
2328 ! ... solve the system
2329 !----------------------------------------------------------------------
2330       x(:,1) = x(:,1)*main(:,1)
2331       do i = 2,nm1
2332          x(:,i) = (x(:,i) - lower(:,i)*x(:,i-1))*main(:,i)
2333       end do
2334       x(:,n) = (x(:,n) - lower(:,n)*x(:,nm1))/(main(:,n) - lower(:,n)*upper(:,nm1))
2335       do i = nm1,1,-1
2336          x(:,i) = x(:,i) - upper(:,i)*x(:,i+1)
2337       end do
2338
2339       end subroutine trislv
2340
2341       subroutine ps2str( nz, nw, zen, rsfc, tauu, omu, &
2342                          gu, dsdh, nid, radfld )
2343 !-----------------------------------------------------------------------------
2344 ! purpose:
2345 ! solve two-stream equations for multiple layers. the subroutine is based
2346 ! on equations from: toon et al., j.geophys.res., v94 (d13), nov 20, 1989.
2347 ! it contains 9 two-stream methods to choose from. a pseudo-spherical
2348 ! correction has also been added.
2349 !-----------------------------------------------------------------------------
2350 ! parameters:
2351 ! nz - integer, number of specified altitude levels in the working (i)
2352 ! grid
2353 ! zen - real(dp), solar zenith angle (degrees) (i)
2354 ! rsfc - real(dp), surface albedo at current wavelength (i)
2355 ! tauu - real(dp), unscaled optical depth of each layer (i)
2356 ! omu - real(dp), unscaled single scattering albedo of each layer (i)
2357 ! gu - real(dp), unscaled asymmetry parameter of each layer (i)
2358 ! dsdh - real(dp), slant path of direct beam through each layer crossed (i)
2359 ! when travelling from the top of the atmosphere to layer i;
2360 ! dsdh(i,j), i = 0..nz-1, j = 1..nz-1
2361 ! nid - integer, number of layers crossed by the direct beam when (i)
2362 ! travelling from the top of the atmosphere to layer i;
2363 ! nid(i), i = 0..nz-1
2364 ! delta - logical, switch to use delta-scaling (i)
2365 ! .true. -> apply delta-scaling
2366 ! .false.-> do not apply delta-scaling
2367 ! fdr - real(dp), contribution of the direct component to the total (o)
2368 ! actinic flux at each altitude level
2369 ! fup - real(dp), contribution of the diffuse upwelling component to (o)
2370 ! the total actinic flux at each altitude level
2371 ! fdn - real(dp), contribution of the diffuse downwelling component to (o)
2372 ! the total actinic flux at each altitude level
2373 ! edr - real(dp), contribution of the direct component to the total (o)
2374 ! spectral irradiance at each altitude level
2375 ! eup - real(dp), contribution of the diffuse upwelling component to (o)
2376 ! the total spectral irradiance at each altitude level
2377 ! edn - real(dp), contribution of the diffuse downwelling component to (o)
2378 ! the total spectral irradiance at each altitude level
2379 !-----------------------------------------------------------------------------
2380
2381       implicit none
2382 !-----------------------------------------------------------------------------
2383 ! ... dummy arguments
2384 !-----------------------------------------------------------------------------
2385       integer, intent(in) :: nz, nw
2386       integer, intent(in) :: nid(0:nz-1)
2387       real(dp), intent(in) :: zen
2388       real(dp), intent(in) :: rsfc(nw-1)
2389       real(dp), dimension(nz-1,nw-1), intent(in) :: tauu, omu, gu
2390       real(dp), intent(in) :: dsdh(0:nz-1,nz-1)
2391       real(dp), intent(out) :: radfld(nz,nw-1)
2392 !-----------------------------------------------------------------------------
2393 ! ... local variables
2394 !-----------------------------------------------------------------------------
2395 !-----------------------------------------------------------------------------
2396 ! ... mu = cosine of solar zenith angle
2397 ! rsfc = surface albedo
2398 ! tauu = unscaled optical depth of each layer
2399 ! omu = unscaled single scattering albedo
2400 ! gu = unscaled asymmetry factor
2401 ! klev = max dimension of number of layers in atmosphere
2402 ! nlayer = number of layers in the atmosphere
2403 ! nlevel = nlayer + 1 = number of levels
2404 !-----------------------------------------------------------------------------
2405       real(dp), parameter :: smallest = tiny( 1.0_dp )
2406       real(dp), parameter :: largest  = huge( 1.0_dp )
2407       real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp
2408       integer :: mrows, nzm1, nzm2
2409       real(dp), parameter :: eps = 1.e-3_dp
2410       real(dp), parameter :: pifs = 1._dp
2411       real(dp), parameter :: fdn0 = 0._dp
2412       integer :: row
2413       integer :: lev
2414       integer :: i, ip1, iw
2415       integer :: j, jl, ju
2416       real(dp) :: precis, wrk
2417       real(dp) :: tempg
2418       real(dp) :: mu, suma
2419       real(dp) :: g, om
2420       real(dp) :: gam1, gam2, gam3, gam4
2421       real(dp), dimension(nz-1) :: f, gi, omi
2422       real(dp), dimension(0:nz) :: tauc, mu2
2423       real(dp), dimension(nz-1) :: lam, taun, bgam
2424       real(dp), dimension(nz-1) :: cdn
2425       real(dp), dimension(0:nz,nw-1) :: tausla
2426       real(dp), dimension(nz-1,nw-1) :: cup, cuptn, cdntn
2427       real(dp), dimension(nz-1,nw-1) :: e1, e2, e3, e4
2428       real(dp), dimension(2*(nz-1)) :: a, b, d, e
2429       real(dp), dimension(nw-1,2*(nz-1)) :: sub, main, super, y
2430 !-----------------------------------------------------------------------------
2431 ! ... for calculations of associated legendre polynomials for gama1,2,3,4
2432 ! in delta-function, modified quadrature, hemispheric constant,
2433 ! hybrid modified eddington-delta function metods, p633,table1.
2434 ! w.e.meador and w.r.weaver, gas,1980,v37,p.630
2435 ! w.j.wiscombe and g.w. grams, gas,1976,v33,p2440,
2436 ! uncomment the following two lines and the appropriate statements
2437 ! further down.
2438 !-----------------------------------------------------------------------------
2439       real(dp) :: expon, expon0, expon1, divisr, temp, up, dn
2440       real(dp) :: ssfc
2441 !-----------------------------------------------------------------------------
2442 ! ... initial conditions: pi*solar flux = 1; diffuse incidence = 0
2443 !-----------------------------------------------------------------------------
2444       nzm1 = nz - 1
2445       nzm2 = nz - 2
2446       mrows = 2*nzm1
2447       precis = epsilon( precis )
2448       mu = cos( zen*d2r )
2449 wave_loop : &
2450       do iw = 1,nw-1
2451 !-----------------------------------------------------------------------------
2452 ! ... compute coefficients for each layer:
2453 ! gam1 - gam4 = 2-stream coefficients, different for different approximations
2454 ! expon0 = calculation of e when tau is zero
2455 ! expon1 = calculation of e when tau is taun
2456 ! cup and cdn = calculation when tau is zero
2457 ! cuptn and cdntn = calc. when tau is taun
2458 ! divisr = prevents division by zero
2459 !-----------------------------------------------------------------------------
2460          tauc(0:nz) = 0._dp
2461          tausla(0:nz,iw) = 0._dp
2462          mu2(0:nz) = sqrt( smallest )
2463 !-----------------------------------------------------------------------------
2464 ! ... delta-scaling. have to be done for delta-eddington approximation,
2465 ! delta discrete ordinate, practical improved flux method, delta function,
2466 ! and hybrid modified eddington-delta function methods approximations
2467 !-----------------------------------------------------------------------------
2468          f(1:nzm1) = gu(:,iw)*gu(:,iw)
2469          gi(1:nzm1) = (gu(:,iw) - f(1:nzm1))/(1._dp - f(1:nzm1))
2470          omi(1:nzm1) = (1._dp - f(1:nzm1))*omu(1:nzm1,iw)/(1._dp - omu(1:nzm1,iw)*f(1:nzm1))
2471          taun(1:nzm1) = (1._dp - omu(1:nzm1,iw)*f(1:nzm1))*tauu(1:nzm1,iw)
2472 !-----------------------------------------------------------------------------
2473 ! ... calculate slant optical depth at the top of the atmosphere when zen>90.
2474 ! in this case, higher altitude of the top layer is recommended which can
2475 ! be easily changed in gridz.f.
2476 !-----------------------------------------------------------------------------
2477          if( zen > 90._dp ) then
2478             if( nid(0) < 0 ) then
2479                tausla(0,iw) = largest
2480             else
2481                ju = nid(0)
2482                tausla(0,iw) = 2._dp*dot_product( taun(1:ju),dsdh(0,1:ju) )
2483             end if
2484          end if
2485 level_loop : &
2486          do i = 1,nzm1
2487             g = gi(i)
2488             om = omi(i)
2489             tauc(i) = tauc(i-1) + taun(i)
2490 !-----------------------------------------------------------------------------
2491 ! ... stay away from 1 by precision. for g, also stay away from -1
2492 !-----------------------------------------------------------------------------
2493             tempg = min( abs(g),1._dp - precis )
2494             g = sign( tempg,g )
2495             om = min( om,1._dp - precis )
2496 !-----------------------------------------------------------------------------
2497 ! ... calculate slant optical depth
2498 !-----------------------------------------------------------------------------
2499             if( nid(i) < 0 ) then
2500                tausla(i,iw) = largest
2501             else
2502                ju = min( nid(i),i )
2503                suma = dot_product( taun(1:ju),dsdh(i,1:ju) )
2504                jl = min( nid(i),i ) + 1
2505                tausla(i,iw) = suma + 2._dp*dot_product( taun(jl:nid(i)),dsdh(i,jl:nid(i)) )
2506                if( abs(tausla(i,iw)-tausla(i-1,iw)) < 1.0e-30_dp ) then
2507                  mu2(i) = sqrt( largest )
2508                else
2509                  mu2(i) = (tauc(i) - tauc(i-1))/(tausla(i,iw) - tausla(i-1,iw))
2510                  mu2(i) = sign( max( abs(mu2(i)),sqrt(smallest) ),mu2(i) )
2511                end if
2512             end if
2513 !-----------------------------------------------------------------------------
2514 ! ... the following gamma equations are from pg 16,289, table 1
2515 ! eddington approximation(joseph et al., 1976, jas, 33, 2452):
2516 !-----------------------------------------------------------------------------
2517             gam1 = .25_dp*(7._dp - om*(4._dp + 3._dp*g))
2518             gam2 = -.25_dp*(1._dp - om*(4._dp - 3._dp*g))
2519             gam3 = .25_dp*(2._dp - 3._dp*g*mu)
2520             gam4 = 1._dp - gam3
2521 !-----------------------------------------------------------------------------
2522 ! ... lambda = pg 16,290 equation 21
2523 ! big gamma = pg 16,290 equation 22
2524 !-----------------------------------------------------------------------------
2525             lam(i) = sqrt( gam1*gam1 - gam2*gam2 )
2526             bgam(i) = (gam1 - lam(i))/gam2
2527             wrk = lam(i)*taun(i)
2528             if( wrk < 500._dp ) then
2529                expon = exp( -wrk )
2530             else
2531                expon = 0._dp
2532             end if
2533 !-----------------------------------------------------------------------------
2534 ! ... e1 - e4 = pg 16,292 equation 44
2535 !-----------------------------------------------------------------------------
2536             e1(i,iw) = 1._dp + bgam(i)*expon
2537             e2(i,iw) = 1._dp - bgam(i)*expon
2538             e3(i,iw) = bgam(i) + expon
2539             e4(i,iw) = bgam(i) - expon
2540 !-----------------------------------------------------------------------------
2541 ! ... the following sets up for the c equations 23, and 24
2542 ! found on page 16,290
2543 ! prevent division by zero (if lambda=1/mu, shift 1/mu^2 by eps = 1.e-3
2544 ! which is approx equiv to shifting mu by 0.5*eps* (mu)**3
2545 !-----------------------------------------------------------------------------
2546             if( tausla(i-1,iw) < 500._dp ) then
2547                expon0 = exp( -tausla(i-1,iw) )
2548             else
2549                expon0 = 0._dp
2550             end if
2551             if( tausla(i,iw) < 500._dp ) then
2552                expon1 = exp( -tausla(i,iw) )
2553             else
2554                expon1 = 0._dp
2555             end if
2556
2557             divisr = lam(i)*lam(i) - 1._dp/(mu2(i)*mu2(i))
2558             temp = max( eps,abs(divisr) )
2559             divisr = 1._dp/sign( temp,divisr )
2560             up = om*pifs*((gam1 - 1._dp/mu2(i))*gam3 + gam4*gam2)*divisr
2561             dn = om*pifs*((gam1 + 1._dp/mu2(i))*gam4 + gam2*gam3)*divisr
2562 !-----------------------------------------------------------------------------
2563 ! ... cup and cdn are when tau is equal to zero
2564 ! cuptn and cdntn are when tau is equal to taun
2565 !-----------------------------------------------------------------------------
2566             cup(i,iw) = up*expon0
2567             cdn(i) = dn*expon0
2568             cuptn(i,iw) = up*expon1
2569             cdntn(i,iw) = dn*expon1
2570
2571          end do level_loop
2572 !-----------------------------------------------------------------------------
2573 ! ... set up matrix
2574 ! ssfc = pg 16,292 equation 37 where pi fs is one (unity).
2575 !-----------------------------------------------------------------------------
2576         if( tausla(nzm1,iw) < 500._dp ) then
2577            ssfc = rsfc(iw)*mu*exp( -tausla(nzm1,iw) )*pifs
2578         else
2579            ssfc = 0._dp
2580         end if
2581 !-----------------------------------------------------------------------------
2582 ! ... the following are from pg 16,292 equations 39 - 43.
2583 ! set up first row of matrix:
2584 !-----------------------------------------------------------------------------
2585         a(1) = 0._dp
2586         b(1) = e1(1,iw)
2587         d(1) = -e2(1,iw)
2588         e(1) = fdn0 - cdn(1)
2589 !-----------------------------------------------------------------------------
2590 ! ... set up odd rows 3 thru (mrows - 1):
2591 !-----------------------------------------------------------------------------
2592         a(3:mrows-1:2) = e2(1:nzm2,iw)*e3(1:nzm2,iw) - e4(1:nzm2,iw)*e1(1:nzm2,iw)
2593         b(3:mrows-1:2) = e1(1:nzm2,iw)*e1(2:nzm1,iw) - e3(1:nzm2,iw)*e3(2:nzm1,iw)
2594         d(3:mrows-1:2) = e3(1:nzm2,iw)*e4(2:nzm1,iw) - e1(1:nzm2,iw)*e2(2:nzm1,iw)
2595         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))
2596 !-----------------------------------------------------------------------------
2597 ! ... set up even rows 2 thru (mrows - 2):
2598 !-----------------------------------------------------------------------------
2599         a(2:mrows-2:2) = e2(2:nzm1,iw)*e1(1:nzm2,iw) - e3(1:nzm2,iw)*e4(2:nzm1,iw)
2600         b(2:mrows-2:2) = e2(1:nzm2,iw)*e2(2:nzm1,iw) - e4(1:nzm2,iw)*e4(2:nzm1,iw)
2601         d(2:mrows-2:2) = e1(2:nzm1,iw)*e4(2:nzm1,iw) - e2(2:nzm1,iw)*e3(2:nzm1,iw)
2602         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)
2603 !-----------------------------------------------------------------------------
2604 ! ... set up last row of matrix at mrows:
2605 !-----------------------------------------------------------------------------
2606         a(mrows) = e1(nzm1,iw) - rsfc(iw)*e3(nzm1,iw)
2607         b(mrows) = e2(nzm1,iw) - rsfc(iw)*e4(nzm1,iw)
2608         d(mrows) = 0._dp
2609         e(mrows) = ssfc - cuptn(nzm1,iw) + rsfc(iw)*cdntn(nzm1,iw)
2610         sub(iw,1:mrows) = a(1:mrows)
2611         main(iw,1:mrows) = b(1:mrows)
2612         super(iw,1:mrows) = d(1:mrows)
2613         y(iw,1:mrows) = e(1:mrows)
2614       end do wave_loop
2615 !-----------------------------------------------------------------------------
2616 ! ... solve the system
2617 !-----------------------------------------------------------------------------
2618       call tridec( nw-1, mrows, sub, main, super )
2619       call trislv( nw-1, mrows, sub, main, super, y )
2620 !-----------------------------------------------------------------------------
2621 ! ... unfold solution of matrix, compute output fluxes
2622 !-----------------------------------------------------------------------------
2623       do iw = 1,nw-1
2624 !-----------------------------------------------------------------------------
2625 ! ... the following equations are from pg 16,291 equations 31 & 32
2626 !-----------------------------------------------------------------------------
2627          e(:mrows) = y(iw,:mrows)
2628          if( tausla(0,iw) < 500._dp ) then
2629             radfld(1,iw) = 2._dp*(fdn0 + e(1)*e3(1,iw) - e(2)*e4(1,iw) + cup(1,iw)) + exp( -tausla(0,iw) )
2630          else
2631             radfld(1,iw) = 2._dp*(fdn0 + e(1)*e3(1,iw) - e(2)*e4(1,iw) + cup(1,iw))
2632          end if
2633          where( tausla(1:nzm1,iw) < 500._dp )
2634             cdn(1:nzm1) = exp( -tausla(1:nzm1,iw) )
2635          elsewhere
2636             cdn(1:nzm1) = 0._dp
2637          endwhere
2638          radfld(2:nz,iw) = 2._dp*(e(1:mrows-1:2)*(e3(1:nzm1,iw) + e1(1:nzm1,iw)) &
2639                             + e(2:mrows:2)*(e4(1:nzm1,iw) + e2(1:nzm1,iw)) &
2640                             + cdntn(1:nzm1,iw) + cuptn(1:nzm1,iw)) + cdn(1:nzm1)
2641       end do
2642
2643       end subroutine ps2str
2644
2645
2646       subroutine pchem( chem_opt, nz, tlev, airlev )
2647
2648       use module_state_description, only : MOZART_KPP, MOZCART_KPP, T1_MOZCART_KPP, &
2649                                            MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP
2650       use module_wave_data, only : r01, r04, r44, r08, r06, r10,  &
2651                                    r11, r14, r15, r17, r18,       &
2652                                    xs_mvk, xs_macr, xs_hyac, xs_glyald, nj
2653
2654 !-----------------------------------------------------------------------------
2655 ! purpose:
2656 ! load various "weighting functions" (products of cross section and
2657 ! quantum yield at each altitude and for wavelength). the altitude
2658 ! dependence is necessary to ensure the consideration of pressure and
2659 ! temperature dependence of the cross sections or quantum yields.
2660 ! the actual reading, evaluation and interpolation is done is separate
2661 ! subroutines for ease of management and manipulation. please refer to
2662 ! the inline documentation of the specific subroutines for detail information.
2663 !-----------------------------------------------------------------------------
2664 ! parameters:
2665 ! nw - integer, number of specified intervals + 1 in working (i)
2666 ! wavelength grid
2667 ! wl - real(dp), vector of lower limits of wavelength intervals in (i)
2668 ! working wavelength grid
2669 ! nz - integer, number of altitude levels in working altitude grid (i)
2670 ! tlev - real(dp), temperature (k) at each specified altitude level (i)
2671 ! airlev - real(dp), air density (molec/cc) at each altitude level (i)
2672 ! j - integer, counter for number of weighting functions defined (io)
2673 ! sq - real(dp), cross section x quantum yield (cm^2) for each (o)
2674 ! photolysis reaction defined, at each defined wavelength and
2675 ! at each defined altitude level
2676 ! jlabel - character*40, string identifier for each photolysis reaction (o)
2677 ! defined
2678 !-----------------------------------------------------------------------------
2679
2680       implicit none
2681
2682 !-----------------------------------------------------------------------------
2683 ! ... dummy arguments
2684 !-----------------------------------------------------------------------------
2685       integer, intent(in)  :: chem_opt
2686       integer, intent(in)  :: nz
2687       real(dp), intent(in) :: tlev(nz)
2688       real(dp), intent(in) :: airlev(nz)
2689
2690 !-----------------------------------------------------------------------------
2691 ! ... local variables
2692 !-----------------------------------------------------------------------------
2693       integer :: j
2694       character(len=40) :: jlabel(nj)
2695 !-----------------------------------------------------------------------------
2696 ! ... (1) o2 + hv -> o + o
2697 ! reserve first position. cross section parameterization in
2698 ! schumman-runge and lyman-alpha regions are zenith-angle dependent
2699 !-----------------------------------------------------------------------------
2700       j = 1
2701       jlabel(j) = 'o2 + hv -> o + o '
2702 !-----------------------------------------------------------------------------
2703 ! ... (2,3) o3 + hv -> (both channels) (tmp dep)
2704 !-----------------------------------------------------------------------------
2705       call r01( nz, tlev, airlev, j, jlabel )
2706 !-----------------------------------------------------------------------------
2707 ! ... (4) no2 + hv -> no + o(3p)
2708 !-----------------------------------------------------------------------------
2709       j = j + 1
2710       jlabel(j) = 'no2 + hv -> no + o(3p) '
2711 !-----------------------------------------------------------------------------
2712 ! ... (5) no3 + hv -> no2 + o(3p)
2713 !-----------------------------------------------------------------------------
2714       j = j + 1
2715       jlabel(j) = 'no3 + hv -> no2 + o(3p) '
2716 !-----------------------------------------------------------------------------
2717 ! ... (6) no3 + hv -> no + o2
2718 !-----------------------------------------------------------------------------
2719       j = j + 1
2720       jlabel(j) = 'no3 + hv -> no + o2 '
2721 !-----------------------------------------------------------------------------
2722 ! ... (7,8) n2o5 + hv -> (both channels) (tmp dep)
2723 !-----------------------------------------------------------------------------
2724       call r04( nz,tlev,airlev,j,jlabel )
2725 !-----------------------------------------------------------------------------
2726 ! ... (9) n2o + hv -> n2 + o(1d) (tmp dep)
2727 !-----------------------------------------------------------------------------
2728       call r44( nz,tlev,airlev,j,jlabel )
2729 !-----------------------------------------------------------------------------
2730 ! ... (10) ho2 + hv -> oh + o
2731 !-----------------------------------------------------------------------------
2732       j = j + 1
2733       jlabel(j) = 'ho2 + hv -> oh + o '
2734 !-----------------------------------------------------------------------------
2735 ! ... (11) h2o2 + hv -> 2 oh
2736 !-----------------------------------------------------------------------------
2737       call r08( nz,tlev,airlev,j,jlabel )
2738 !-----------------------------------------------------------------------------
2739 ! ... (12) hno2 + hv -> oh + no
2740 !-----------------------------------------------------------------------------
2741       j = j + 1
2742       jlabel(j) = 'hno2 + hv -> oh + no '
2743 !-----------------------------------------------------------------------------
2744 ! ... (13) hno3 + hv -> oh + no2
2745 !-----------------------------------------------------------------------------
2746       call r06( nz,tlev,airlev,j,jlabel )
2747 !-----------------------------------------------------------------------------
2748 ! ... (14) hno4 + hv -> ho2 + no2
2749 !-----------------------------------------------------------------------------
2750       j = j + 1
2751       jlabel(j) = 'hno4 + hv -> ho2 + no2 '
2752 !-----------------------------------------------------------------------------
2753 ! ... (15,16) ch2o + hv -> (both channels)
2754 !-----------------------------------------------------------------------------
2755       call r10( nz,tlev,airlev,j,jlabel )
2756 !-----------------------------------------------------------------------------
2757 ! ... (17,18,19) ch3cho + hv -> (all three channels)
2758 !-----------------------------------------------------------------------------
2759       call r11( nz,tlev,airlev,j,jlabel )
2760 !-----------------------------------------------------------------------------
2761 ! ... (20) c2h5cho + hv -> c2h5 + hco
2762 !-----------------------------------------------------------------------------
2763       j = j + 1
2764       jlabel(j) = 'c2h5cho + hv -> c2h5 + hco '
2765 !-----------------------------------------------------------------------------
2766 ! ... (21) chocho + hv -> products
2767 !-----------------------------------------------------------------------------
2768       j = j + 1
2769       jlabel(j) = 'chocho + hv -> products '
2770 !-----------------------------------------------------------------------------
2771 ! ... (22) ch3cocho + hv -> products
2772 !-----------------------------------------------------------------------------
2773       call r14( nz,tlev,airlev,j,jlabel )
2774 !-----------------------------------------------------------------------------
2775 ! ... (23) ch3coch3 + hv -> products
2776 !-----------------------------------------------------------------------------
2777       call r15( nz,tlev,airlev,j,jlabel )
2778 !-----------------------------------------------------------------------------
2779 ! ... (24) ch3ooh + hv -> ch3o + oh
2780 !-----------------------------------------------------------------------------
2781       j = j + 1
2782       jlabel(j) = 'ch3ooh + hv -> ch3o + oh '
2783 !-----------------------------------------------------------------------------
2784 ! ... (25) ch3ono2 + hv -> ch3o + no2
2785 !-----------------------------------------------------------------------------
2786       call r17( nz,tlev,airlev,j,jlabel )
2787 !-----------------------------------------------------------------------------
2788 ! ... (26) pan + hv -> products
2789 !-----------------------------------------------------------------------------
2790       call r18( nz,tlev,airlev,j,jlabel )
2791 !-----------------------------------------------------------------------------
2792 ! ... (27) mvk + hv -> products
2793 !     (28) macr + hv -> products
2794 !     (29) hyac + hv -> products
2795 !     (30) glyald + hv -> products
2796 !-----------------------------------------------------------------------------
2797 is_mozart : &
2798       if( chem_opt == MOZART_KPP .or. chem_opt == MOZCART_KPP .or. &
2799           chem_opt == T1_MOZCART_KPP .or. &
2800           chem_opt == MOZART_MOSAIC_4BIN_KPP .or. &
2801           chem_opt == MOZART_MOSAIC_4BIN_AQ_KPP) then
2802         call xs_mvk( nz, tlev, airlev )
2803         call xs_macr( nz, tlev, airlev )
2804         call xs_hyac( nz, tlev, airlev )
2805         call xs_glyald( nz, tlev, airlev )
2806       else is_mozart
2807 !-----------------------------------------------------------------------------
2808 ! ... (27) cloo + hv -> products
2809 !-----------------------------------------------------------------------------
2810       j = j + 1
2811       jlabel(j) = 'cloo + hv -> products '
2812 !-----------------------------------------------------------------------------
2813 ! ... (28,29) clono2 + hv -> products
2814 !-----------------------------------------------------------------------------
2815       j = j + 1
2816       jlabel(j) = 'clono2 + hv -> cl + no3'
2817       j = j + 1
2818       jlabel(j) = 'clono2 + hv -> clo + no2'
2819 !-----------------------------------------------------------------------------
2820 ! ... (30) ch3cl + hv -> products
2821 !-----------------------------------------------------------------------------
2822       j = j + 1
2823       jlabel(j) = 'ch3cl + hv -> products '
2824 !-----------------------------------------------------------------------------
2825 ! ... (31) ccl2o + hv -> products
2826 !-----------------------------------------------------------------------------
2827       j = j + 1
2828       jlabel(j) = 'ccl2o + hv -> products'
2829 !-----------------------------------------------------------------------------
2830 ! ... (32) ccl4 + hv -> products
2831 !-----------------------------------------------------------------------------
2832       j = j + 1
2833       jlabel(j) = 'ccl4 + hv -> products '
2834 !-----------------------------------------------------------------------------
2835 ! ... (33) cclfo + hv -> products
2836 !-----------------------------------------------------------------------------
2837       j = j + 1
2838       jlabel(j) = 'cclfo + hv -> products '
2839 !-----------------------------------------------------------------------------
2840 ! ... (34) ccf2o + hv -> products
2841 !-----------------------------------------------------------------------------
2842       j = j + 1
2843       jlabel(j) = 'ccf2o + hv -> products '
2844 !-----------------------------------------------------------------------------
2845 ! ... (35) cf2clcfcl2 (cfc-113) + hv -> products
2846 !-----------------------------------------------------------------------------
2847       j = j + 1
2848       jlabel(j) = 'cf2clcfcl2 (cfc-113) + hv -> products'
2849 !-----------------------------------------------------------------------------
2850 ! ... (36) cf2clcf2cl (cfc-114) + hv -> products
2851 !-----------------------------------------------------------------------------
2852       j = j + 1
2853       jlabel(j) = 'cf2clcf2cl (cfc-114) + hv -> products '
2854 !-----------------------------------------------------------------------------
2855 ! ... (37) cf3cf2cl (cfc-115) + hv -> products
2856 !-----------------------------------------------------------------------------
2857       j = j + 1
2858       jlabel(j) = 'cf3cf2cl (cfc-115) + hv -> products '
2859 !-----------------------------------------------------------------------------
2860 ! ... (38) ccl3f (cfc-111) + hv -> products
2861 !-----------------------------------------------------------------------------
2862       j = j + 1
2863       jlabel(j) = 'ccl3f (cfc-111) + hv -> products'
2864 !-----------------------------------------------------------------------------
2865 ! ... (39) ccl2f2 (cfc-112) + hv -> products
2866 !-----------------------------------------------------------------------------
2867       j = j + 1
2868       jlabel(j) = 'ccl2f2 (cfc-112) + hv -> products'
2869 !-----------------------------------------------------------------------------
2870 ! (40) ch3ccl3 + hv -> products
2871 !-----------------------------------------------------------------------------
2872       j = j + 1
2873       jlabel(j)='ch3ccl3 + hv -> products'
2874 !-----------------------------------------------------------------------------
2875 ! (41) cf3chcl2 (hcfc-123) + hv -> products
2876 !-----------------------------------------------------------------------------
2877       j = j + 1
2878       jlabel(j)='cf3chcl2 (hcfc-123) + hv -> products'
2879 !-----------------------------------------------------------------------------
2880 ! (42) cf3chfcl (hcfc-124) + hv -> products
2881 !-----------------------------------------------------------------------------
2882       j = j + 1
2883       jlabel(j)='cf3chfcl (hcfc-124) + hv -> products'
2884 !-----------------------------------------------------------------------------
2885 ! (43) ch3cfcl2 (hcfc-141b) + hv -> products
2886 !-----------------------------------------------------------------------------
2887       j = j + 1
2888       jlabel(j)='ch3cfcl2 (hcfc-141b) + hv -> products'
2889 !-----------------------------------------------------------------------------
2890 ! (44) ch3cf2cl (hcfc-142b) + hv -> products
2891 !-----------------------------------------------------------------------------
2892       j = j + 1
2893       jlabel(j)='ch3cf2cl (hcfc-142b) + hv -> products'
2894 !-----------------------------------------------------------------------------
2895 ! (45) cf3cf2chcl2 (hcfc-225ca) + hv -> products
2896 !-----------------------------------------------------------------------------
2897       j = j + 1
2898       jlabel(j)='cf3cf2chcl2 (hcfc-225ca) + hv -> products'
2899 !-----------------------------------------------------------------------------
2900 ! (46) cf2clcf2chfcl (hcfc-225cb) + hv -> products
2901 !-----------------------------------------------------------------------------
2902       j = j + 1
2903       jlabel(j)='cf2clcf2chfcl (hcfc-225cb) + hv -> products'
2904 !-----------------------------------------------------------------------------
2905 ! (47) chclf2 (hcfc-22) + hv -> products
2906 !-----------------------------------------------------------------------------
2907       j = j + 1
2908       jlabel(j)='chclf2 (hcfc-22) + hv -> products'
2909 !-----------------------------------------------------------------------------
2910 ! (48) brono2 + hv -> br + o + no2
2911 !-----------------------------------------------------------------------------
2912       j = j + 1
2913       jlabel(j)='brono2 + hv -> br + o + no2'
2914 !-----------------------------------------------------------------------------
2915 ! (49) brono2 + hv -> bro + no2
2916 !-----------------------------------------------------------------------------
2917       j = j + 1
2918       jlabel(j)='brono2 + hv -> bro + no2'
2919 !-----------------------------------------------------------------------------
2920 ! (50) ch3br + hv -> products
2921 !-----------------------------------------------------------------------------
2922       j = j + 1
2923       jlabel(j)=' ch3br + hv -> products'
2924 !-----------------------------------------------------------------------------
2925 ! (51) chbr3 + hv -> products
2926 !-----------------------------------------------------------------------------
2927       j = j + 1
2928       jlabel(j)='chbr3 + hv -> products'
2929 !-----------------------------------------------------------------------------
2930 ! (52) cf3br (halon-1301) + hv -> products
2931 !-----------------------------------------------------------------------------
2932       j = j + 1
2933       jlabel(j)='cf3br (halon-1301) + hv -> products'
2934 !-----------------------------------------------------------------------------
2935 ! (53) cf2brcf2br (halon-2402) + hv -> products
2936 !-----------------------------------------------------------------------------
2937       j = j + 1
2938       jlabel(j)='cf2brcf2br (halon-2402) + hv -> products'
2939 !-----------------------------------------------------------------------------
2940 ! (54) cf2br2 (halon-1202) + hv -> products
2941 !-----------------------------------------------------------------------------
2942       j = j + 1
2943       jlabel(j)='cf2br2 (halon-1202) + hv -> products'
2944 !-----------------------------------------------------------------------------
2945 ! (55) cf2brcl (halon-1211) + hv -> products
2946 !-----------------------------------------------------------------------------
2947       j = j + 1
2948       jlabel(j)='cf2brcl (halon-1211) + hv -> products '
2949 !-----------------------------------------------------------------------------
2950 ! (56) cl2 + hv -> cl + cl
2951 !-----------------------------------------------------------------------------
2952       j = j + 1
2953       jlabel(j)='cl2 + hv -> cl + cl '
2954 !-----------------------------------------------------------------------------
2955 ! (57) hocl + hv -> oh + cl
2956 !-----------------------------------------------------------------------------
2957       j = j + 1
2958       jlabel(j)='hocl + hv -> oh + cl'
2959 !-----------------------------------------------------------------------------
2960 ! (58) fmcl + hv -> cl + co + ho2
2961 !-----------------------------------------------------------------------------
2962       j = j + 1
2963       jlabel(j)='fmcl + hv -> cl + co + ho2'
2964
2965       end if is_mozart
2966
2967
2968       end subroutine pchem
2969
2970       subroutine lymana( nz, o2col, secchi, dto2la, xso2la )
2971 !-----------------------------------------------------------------------------
2972 ! purpose:
2973 ! calculate the effective absorption cross section of o2 in the lyman-alpha
2974 ! bands and an effective o2 optical depth at all altitudes. parameterized
2975 ! after: chabrillat, s., and g. kockarts, simple parameterization of the
2976 ! absorption of the solar lyman-alpha line, geophysical research letters,
2977 ! vol.24, no.21, pp 2659-2662, 1997.
2978 !-----------------------------------------------------------------------------
2979 ! parameters:
2980 ! nz - integer, number of specified altitude levels in the working (i)
2981 ! grid
2982 ! o2col - real(dp), slant overhead o2 column (molec/cc) at each specified (i)
2983 ! altitude
2984 ! dto2la - real(dp), optical depth due to o2 absorption at each specified (o)
2985 ! vertical layer
2986 ! xso2la - real(dp), molecular absorption cross section in la bands (o)
2987 !-----------------------------------------------------------------------------
2988       implicit none
2989 !-----------------------------------------------------------------------------
2990 ! ... dummy arguments
2991 !-----------------------------------------------------------------------------
2992       integer, intent(in) :: nz
2993       real(dp), intent(in) :: o2col(nz)
2994       real(dp), intent(in) :: secchi(nz)
2995       real(dp), intent(out) :: dto2la(nz-1,nla-1)
2996       real(dp), intent(out) :: xso2la(nz,nla-1)
2997 !-----------------------------------------------------------------------------
2998 ! ... local variables
2999 !-----------------------------------------------------------------------------
3000       integer :: i, k, kp1
3001       real(dp), dimension(nz) :: rm, ro2
3002       real(dp), save :: b(3), c(3), d(3), e(3)
3003       data b / 6.8431e-01_dp, 2.29841e-01_dp, 8.65412e-02_dp/, &
3004            c /8.22114e-21_dp, 1.77556e-20_dp, 8.22112e-21_dp/, &
3005            d / 6.0073e-21_dp, 4.28569e-21_dp, 1.28059e-20_dp/, &
3006            e /8.21666e-21_dp, 1.63296e-20_dp, 4.85121e-17_dp/
3007 !-----------------------------------------------------------------------------
3008 ! ... calculate reduction factors at every altitude
3009 !-----------------------------------------------------------------------------
3010       rm(:) = 0._dp
3011       ro2(:) = 0._dp
3012       do k = 1,nz
3013         do i = 1,3
3014           rm(k) = rm(k) + b(i) * exp( -c(i)*o2col(k) )
3015           ro2(k) = ro2(k) + d(i) * exp( -e(i)*o2col(k) )
3016         end do
3017       end do
3018 !-----------------------------------------------------------------------------
3019 ! ... calculate effective o2 optical depths and effective o2 cross sections
3020 !-----------------------------------------------------------------------------
3021       do k = 1,nz-1
3022          if( rm(k) > 1.e-100_dp ) then
3023             kp1 = k + 1
3024             if( rm(kp1) > 0._dp ) then
3025                dto2la(k,1) = log( rm(kp1) )/secchi(kp1) - log( rm(k) )/secchi(k)
3026             else
3027                dto2la(k,1) = 1000._dp
3028             end if
3029          else
3030             dto2la(k,1) = 1000._dp
3031          end if
3032       end do
3033       do k = 1,nz
3034          if( rm(k) > 1.e-100_dp ) then
3035             if( ro2(k) > 1.e-100_dp ) then
3036                xso2la(k,1) = ro2(k)/rm(k)
3037             else
3038                xso2la(k,1) = 0._dp
3039             end if
3040          else
3041             xso2la(k,1) = 0._dp
3042          end if
3043       end do
3044
3045       end subroutine lymana
3046
3047
3048       subroutine inter2( ng, xg, yg, n, x, y, ierr )
3049 !-----------------------------------------------------------------------------
3050 ! purpose:
3051 ! map input data given on single, discrete points onto a set of target
3052 ! bins.
3053 ! the original input data are given on single, discrete points of an
3054 ! arbitrary grid and are being linearly interpolated onto a specified set
3055 ! of target bins. in general, this is the case for most of the weighting
3056 ! functions (action spectra, molecular cross section, and quantum yield
3057 ! data), which have to be matched onto the specified wavelength intervals.
3058 ! the average value in each target bin is found by averaging the trapezoi-
3059 ! dal area underneath the input data curve (constructed by linearly connec-
3060 ! ting the discrete input values).
3061 ! some caution should be used near the endpoints of the grids. if the
3062 ! input data set does not span the range of the target grid, an error
3063 ! message is printed and the execution is stopped, as extrapolation of the
3064 ! data is not permitted.
3065 ! if the input data does not encompass the target grid, use addpnt to
3066 ! expand the input array.
3067 !-----------------------------------------------------------------------------
3068 ! parameters:
3069 ! ng - integer, number of bins + 1 in the target grid (i)
3070 ! xg - real(dp), target grid (e.g., wavelength grid); bin i is defined (i)
3071 ! as [xg(i),xg(i+1)] (i = 1..ng-1)
3072 ! yg - real(dp), y-data re-gridded onto xg, yg(i) specifies the value for (o)
3073 ! bin i (i = 1..ng-1)
3074 ! n - integer, number of points in input grid (i)
3075 ! x - real(dp), grid on which input data are defined (i)
3076 ! y - real(dp), input y-data (i)
3077 !-----------------------------------------------------------------------------
3078       implicit none
3079 !-----------------------------------------------------------------------------
3080 ! ... dummy arguments
3081 !-----------------------------------------------------------------------------
3082       integer, intent(in) :: ng, n
3083       integer, intent(out) :: ierr
3084       real(dp), intent(in) :: x(n)
3085       real(dp), intent(in) :: y(n)
3086       real(dp), intent(in) :: xg(ng)
3087       real(dp), intent(out) :: yg(ng)
3088 !-----------------------------------------------------------------------------
3089 ! ... local variables
3090 !-----------------------------------------------------------------------------
3091       integer :: ngintv
3092       integer :: i, k, jstart
3093       real(dp) :: area, xgl, xgu
3094       real(dp) :: darea, slope
3095       real(dp) :: a1, a2, b1, b2
3096       ierr = 0
3097 !-----------------------------------------------------------------------------
3098 ! ... test for correct ordering of data, by increasing value of x
3099 !-----------------------------------------------------------------------------
3100       do i = 2,n
3101          if( x(i) <= x(i-1) ) then
3102             ierr = 1
3103             write(*,*) 'inter2: x coord not monotonically increasing'
3104             return
3105          end if
3106       end do
3107       do i = 2,ng
3108         if( xg(i) <= xg(i-1) ) then
3109            ierr = 2
3110            write(*,*) 'inter2: xg coord not monotonically increasing'
3111            return
3112         end if
3113       end do
3114 !-----------------------------------------------------------------------------
3115 ! ... check for xg-values outside the x-range
3116 !-----------------------------------------------------------------------------
3117       if( x(1) > xg(1) .or. x(n) < xg(ng) ) then
3118           write(*,*) 'inter2: data does not span grid'
3119           write(*,*) '        use addpnt to expand data and re-run'
3120 ! call endrun
3121       end if
3122 !-----------------------------------------------------------------------------
3123 ! ... find the integral of each grid interval and use this to
3124 ! calculate the average y value for the interval
3125 ! xgl and xgu are the lower and upper limits of the grid interval
3126 !-----------------------------------------------------------------------------
3127       jstart = 1
3128       ngintv = ng - 1
3129       do i = 1,ngintv
3130 !-----------------------------------------------------------------------------
3131 ! ... initalize:
3132 !-----------------------------------------------------------------------------
3133          area = 0._dp
3134          xgl = xg(i)
3135          xgu = xg(i+1)
3136 !-----------------------------------------------------------------------------
3137 ! ... discard data before the first grid interval and after the
3138 ! last grid interval
3139 ! for internal grid intervals, start calculating area by interpolating
3140 ! between the last point which lies in the previous interval and the
3141 ! first point inside the current interval
3142 !-----------------------------------------------------------------------------
3143          k = jstart
3144          if( k <= n-1 ) then
3145 !-----------------------------------------------------------------------------
3146 ! ... if both points are before the first grid, go to the next point
3147 !-----------------------------------------------------------------------------
3148             do
3149                if( x(k+1) <= xgl ) then
3150                   jstart = k - 1
3151                   k = k+1
3152                   if( k <= n-1 ) then
3153                      cycle
3154                   else
3155                      exit
3156                   end if
3157                else
3158                   exit
3159                end if
3160             end do
3161 !-----------------------------------------------------------------------------
3162 ! ... if the last point is beyond the end of the grid,
3163 ! complete and go to the next grid
3164 !-----------------------------------------------------------------------------
3165             do
3166                if( k <= n-1 .and. x(k) < xgu ) then
3167                   jstart = k-1
3168 !-----------------------------------------------------------------------------
3169 ! ... compute x-coordinates of increment
3170 !-----------------------------------------------------------------------------
3171                   a1 = max( x(k),xgl )
3172                   a2 = min( x(k+1),xgu )
3173 !-----------------------------------------------------------------------------
3174 ! ... if points coincide, contribution is zero
3175 !-----------------------------------------------------------------------------
3176                   if( x(k+1) == x(k) ) then
3177                      darea = 0._dp
3178                   else
3179                      slope = (y(k+1) - y(k))/(x(k+1) - x(k))
3180                      b1 = y(k) + slope*(a1 - x(k))
3181                      b2 = y(k) + slope*(a2 - x(k))
3182                      darea = .5*(a2 - a1)*(b2 + b1)
3183                   end if
3184 !-----------------------------------------------------------------------------
3185 ! ... find the area under the trapezoid from a1 to a2
3186 !-----------------------------------------------------------------------------
3187                   area = area + darea
3188                   k = k+1
3189                   cycle
3190                else
3191                   exit
3192                end if
3193             end do
3194          end if
3195 !-----------------------------------------------------------------------------
3196 ! ... calculate the average y after summing the areas in the interval
3197 !-----------------------------------------------------------------------------
3198          yg(i) = area/(xgu - xgl)
3199       end do
3200
3201       end subroutine inter2
3202
3203       subroutine inter_inti( ng, xg, n, x )
3204 !-----------------------------------------------------------------------------
3205 ! ... initialization
3206 !-----------------------------------------------------------------------------
3207       implicit none
3208 !-----------------------------------------------------------------------------
3209 ! ... dummy arguments
3210 !-----------------------------------------------------------------------------
3211       integer, intent(in) :: ng, n
3212       real(dp), intent(in) :: xg(ng)
3213       real(dp), intent(in) :: x(n)
3214 !-----------------------------------------------------------------------------
3215 ! ... local variables
3216 !-----------------------------------------------------------------------------
3217       integer :: i, ii, iil, astat
3218       integer :: ndim(1)
3219       allocate( xi(ng), xcnt(ng-1), stat=astat )
3220       if( astat /= 0 ) then
3221          write(*,*) 'inter_inti: failed to allocate wrk arrays; error = ',astat
3222          stop
3223       else
3224          xi(:) = 0
3225          xcnt(:) = 0
3226       end if
3227       iil = 1
3228       do i = 1,ng
3229          do ii = iil,n-1
3230             if( xg(i) < x(ii) ) then
3231                xi(i) = ii - 1
3232                iil = ii
3233                exit
3234             end if
3235          end do
3236       end do
3237       nintervals = count( xi(:) /= 0 )
3238       if( nintervals == 0 ) then
3239          write(*,*) 'inter_inti: wavelength grids do not overlap'
3240          stop
3241       else
3242          nintervals = nintervals - 1
3243       end if
3244       xcnt(1:nintervals) = xi(2:nintervals+1) - xi(1:nintervals) + 1
3245       ndim(:) = maxval( xcnt(1:nintervals) )
3246       allocate( xfrac(ndim(1),nintervals),stat=astat )
3247       if( astat /= 0 ) then
3248          write(*,*) 'inter_inti: failed to allocate wrk array; error = ',astat
3249          stop
3250       else
3251          xfrac(:,:) = 1._dp
3252       end if
3253       do i = 1,nintervals
3254         iil = xi(i)
3255         xfrac(1,i) = (min( x(iil+1),xg(i+1) ) - xg(i))/(x(iil+1) - x(iil))
3256         if( xcnt(i) > 1 ) then
3257            iil = xi(i) + xcnt(i) - 1
3258            xfrac(xcnt(i),i) = (xg(i+1) - x(iil))/(x(iil+1) - x(iil))
3259         end if
3260       end do
3261 !     write(*,*) 'inter_inti: diagnostics; ng,n,nintervals = ',ng,n,nintervals
3262 !     write(*,'('' xi'')')
3263 !     write(*,'(10i4)') xi(1:nintervals+1)
3264 !     write(*,'('' xcnt'')')
3265 !     write(*,'(10i4)') xcnt(1:nintervals)
3266 !     write(*,'('' xfrac'')')
3267 !     do i = 1,nintervals
3268 !        write(*,'(1p,5e21.13)') xfrac(1:xcnt(i),i)
3269 !     end do
3270
3271       end subroutine inter_inti
3272
3273       subroutine inter3( ng, xg, yg, n, x, y )
3274 !-----------------------------------------------------------------------------
3275 ! purpose:
3276 ! map input data given on a set of bins onto a different set of target
3277 ! bins.
3278 ! the input data are given on a set of bins (representing the integral
3279 ! of the input quantity over the range of each bin) and are being matched
3280 ! onto another set of bins (target grid). a typical example would be an
3281 ! input data set spcifying the extra-terrestrial flux on wavelength inter-
3282 ! vals, that has to be matched onto the working wavelength grid.
3283 ! the resulting area in a given bin of the target grid is calculated by
3284 ! simply adding all fractional areas of the input data that cover that
3285 ! particular target bin.
3286 ! some caution should be used near the endpoints of the grids. if the
3287 ! input data do not span the full range of the target grid, the area in
3288 ! the "missing" bins will be assumed to be zero. if the input data extend
3289 ! beyond the upper limit of the target grid, the user has the option to
3290 ! integrate the "overhang" data and fold the remaining area back into the
3291 ! last target bin. using this option is recommended when re-gridding
3292 ! vertical profiles that directly affect the total optical depth of the
3293 ! model atmosphere.
3294 !-----------------------------------------------------------------------------
3295 ! parameters:
3296 ! ng - integer, number of bins + 1 in the target grid (i)
3297 ! xg - real(dp), target grid (e.g. working wavelength grid); bin i (i)
3298 ! is defined as [xg(i),xg(i+1)] (i = 1..ng-1)
3299 ! yg - real(dp), y-data re-gridded onto xg; yg(i) specifies the (o)
3300 ! y-value for bin i (i = 1..ng-1)
3301 ! n - integer, number of bins + 1 in the input grid (i)
3302 ! x - real(dp), input grid (e.g. data wavelength grid); bin i is (i)
3303 ! defined as [x(i),x(i+1)] (i = 1..n-1)
3304 ! y - real(dp), input y-data on grid x; y(i) specifies the (i)
3305 ! y-value for bin i (i = 1..n-1)
3306 !-----------------------------------------------------------------------------
3307       implicit none
3308 !-----------------------------------------------------------------------------
3309 ! ... dummy arguments
3310 !-----------------------------------------------------------------------------
3311       integer, intent(in) :: n, ng
3312       real(dp), intent(in) :: xg(ng)
3313       real(dp), intent(in) :: x(n)
3314       real(dp), intent(in) :: y(n)
3315       real(dp), intent(out) :: yg(ng)
3316 !-----------------------------------------------------------------------------
3317 ! ... local variables
3318 !-----------------------------------------------------------------------------
3319       integer :: i, ii, iil
3320 !-----------------------------------------------------------------------------
3321 ! ... do interpolation
3322 !-----------------------------------------------------------------------------
3323       yg(:) = 0.
3324       do i = 1,nintervals
3325          iil = xi(i)
3326          ii = xcnt(i)
3327          if( ii == 1 ) then
3328             yg(i) = xfrac(1,i)*y(iil)
3329          else
3330             yg(i) = dot_product( xfrac(1:ii,i),y(iil:iil+ii-1) )
3331          end if
3332       end do
3333
3334       end subroutine inter3
3335
3336       subroutine calcoe( nz, c, xz, tt, adjin, adjcoe )
3337 !-----------------------------------------------------------------------------
3338 ! parameters:
3339 ! adjcoe - real(dp), coross section adjust coefficients (in and out)
3340 ! c(5,28)-polynomal coef
3341 ! tt -nomarlized temperature
3342 !-----------------------------------------------------------------------------*
3343       implicit none
3344 !-----------------------------------------------------------------------------
3345 ! ... dummy arguments
3346 !-----------------------------------------------------------------------------
3347       integer, intent(in) :: nz
3348       real(dp), intent(in) :: adjin
3349       real(dp), intent(in) :: tt
3350       real(dp), intent(in) :: c(5)
3351       real(dp), intent(in) :: xz(nz)
3352       real(dp), intent(inout) :: adjcoe(nz)
3353 !-----------------------------------------------------------------------------
3354 ! ... local variables
3355 !-----------------------------------------------------------------------------
3356       integer :: k
3357       real(dp) :: x
3358       do k = 1,nz
3359          x = xz(k)
3360          adjcoe(k) = adjin * (1._dp + .01_dp*(c(1) + x*(c(2) + x*(c(3) + x*(c(4) + x*c(5))))))
3361       end do
3362
3363       end subroutine calcoe
3364
3365
3366       subroutine airmas( nz, z, zen, dsdh, nid, cz, &
3367                          vcol, scol )
3368 !-----------------------------------------------------------------------------
3369 ! purpose:
3370 ! calculate vertical and slant air columns, in spherical geometry, as a
3371 ! function of altitude.
3372 !-----------------------------------------------------------------------------
3373 ! parameters:
3374 ! nz - integer, number of specified altitude levels in the working (i)
3375 ! grid
3376 ! z - real(dp), specified altitude working grid (km) (i)
3377 ! zen - real(dp), solar zenith angle (degrees) (i)
3378 ! dsdh - real(dp), slant path of direct beam through each layer crossed (o)
3379 ! when travelling from the top of the atmosphere to layer i;
3380 ! dsdh(i,j), i = 0..nz-1, j = 1..nz-1
3381 ! nid - integer, number of layers crossed by the direct beam when (o)
3382 ! travelling from the top of the atmosphere to layer i;
3383 ! nid(i), i = 0..nz-1
3384 ! vcol - real(dp), output, vertical air column, molec cm-2, above level iz
3385 ! scol - real(dp), output, slant air column in direction of sun, above iz
3386 ! also in molec cm-2
3387 !-----------------------------------------------------------------------------
3388       implicit none
3389 !-----------------------------------------------------------------------------
3390 ! ... dummy arguments
3391 !-----------------------------------------------------------------------------
3392       integer, intent(in) :: nz
3393       integer, intent(in) :: nid(0:nz-1)
3394       real(dp), intent(in) :: z(nz)
3395       real(dp), intent(in) :: zen
3396       real(dp), intent(in) :: dsdh(0:nz-1,nz-1)
3397       real(dp), intent(in) :: cz(nz)
3398       real(dp), intent(out) :: vcol(nz)
3399       real(dp), intent(out) :: scol(nz)
3400 !-----------------------------------------------------------------------------
3401 ! ... parameters
3402 !-----------------------------------------------------------------------------
3403       real(dp), parameter :: largest = huge(1.0_dp)
3404 !-----------------------------------------------------------------------------
3405 ! ... local variables
3406 !-----------------------------------------------------------------------------
3407       integer :: id, j
3408       real(dp) :: sum, ssum, vsum, ratio
3409 !-----------------------------------------------------------------------------
3410 ! ... calculate vertical and slant column from each level: work downward
3411 !-----------------------------------------------------------------------------
3412       vsum = 0._dp
3413       ssum = 0._dp
3414       do id = 1,nz-1
3415          vsum = vsum + cz(nz-id)
3416          vcol(nz-id) = vsum
3417          sum = 0.
3418          if( nid(id) < 0 ) then
3419             sum = largest
3420          else
3421 !-----------------------------------------------------------------------------
3422 ! ... single pass layers:
3423 !-----------------------------------------------------------------------------
3424             do j = 1,min( nid(id),id )
3425                sum = sum + cz(nz-j)*dsdh(id,j)
3426             end do
3427 !-----------------------------------------------------------------------------
3428 ! ... double pass layers:
3429 !-----------------------------------------------------------------------------
3430             do j = min( nid(id),id )+1,nid(id)
3431                sum = sum + 2._dp*cz(nz-j)*dsdh(id,j)
3432             end do
3433          end if
3434          scol(nz - id) = sum
3435       end do
3436 !-----------------------------------------------------------------------------
3437 ! ... special section to set scol(nz)
3438 !-----------------------------------------------------------------------------
3439       if( scol(nz-2) /= 0._dp ) then
3440          ratio = scol(nz-1)/scol(nz-2)
3441          scol(nz) = ratio * scol(nz-1)
3442       else
3443          scol(nz) = 0._dp
3444       end if
3445
3446       end subroutine airmas
3447
3448       END MODULE module_ftuv_subs