      !        convert zsurf to km
      !         fix params 
      ! SM NOTES, Aug 2008
      ! deleted subroutine delted
      ! deleted subroutine nlayde
      ! deleted subroutine leqt1b
      ! added subroutines ps2str, tridiag
      ! added subroutines sphers, airmas, for future 
      ! use with better Schumman-Runge parameterization
      ! should be able do delete: nsurf, use 1 instead
    MODULE module_phot_mad
      ! preliminary fixed values for T, p, O3 and caer
      ! use values in correct units!!
      ! so3t(kl,i) - ozone cross sect. temperature dependence coefficients
      ! wl(kl) - array of nominal center wavelengths of spectral intervals
      ! f(kl) - extraterrestrial solar irradiance
      ! xs(kl,nr) - cross sections for nr species.
      ! xqy(kl,nr) - quantum yields.  some are read, others are computed.
      ! schumann-runge - part a
      ! schumann-runge parameters - part b
      ! ozone temperature coefficients for xsection
      ! -----------------------------------------------------------------------
      ! read stand profiles
      ! aerosol
      ! close (incvol)
      ! wavelengths and extraterrestrial irradiance
      ! -----------------------------------------------------
      ! wave_length
      ! WAVE LENGTHS USED BY PHOTOLYSIS PROGRAMS
      ! FILE CREATED AUGUST 19, 1994
      ! FROM MADRONICH 1989 DATA FILE
      ! et_flux
      ! EXTRA-TERRESTIAL IRRADIANCE
      ! FILE CREATED AUGUST 19, 1994
      ! FROM MADRONICH 1989 DATA FILE
      ! -----------------------------------------------------------------------
      ! current jindex assignments - if calculation order changes
      ! these change - subroutine yield must be changed!
      ! process       jindex               notes
      ! wavelength     none           center values, modified wmo grid
      ! e.t.irradian   none           photons per bin, not per nm
      ! absorption cross sections:
      ! o2 absorp       1             schumman-runge corrected in srband
      ! o3 -> 1d        2             at 275 k. correct t-dep in subgrid
      ! o3 -> 3p        3             at 275 k. correct t-dep in subgrid
      ! no2             4
      ! no3 -> no+o2    5
      ! no3 -> no2+o    6
      ! hno2            7
      ! hno3            8
      ! hno4            9
      ! h2o2           10
      ! ch2o -> rad    11
      ! ch2o -> mol    12
      ! ch3cho         13
      ! ch3coch3       14
      ! ch3coc2h5      15
      ! hcocho proc a  16
      ! ch3cocho       17
      ! hcoch=chcho    18              estimate, no reliable measurement
      ! ch3o2h         19
      ! ch3coo2h       20              actually use 0.28*(h2o2 value)
      ! ch3ono2        21
      ! hcocho proc b  22
      ! quantum yields:
      ! no2             4
      ! ch2o -> rad    11
      ! ch2o -> mol    12
      ! ch3cho         13
      ! hcoch=chcho    18              (energetic threshold)
      ! cross section and quantum yield data.  this section assigns
      ! yields for ntp air.  yields are read from data file
      ! some yields are temperature and/or pressure dependent:ch3cho, ch2o_b,
      ! o3.
      ! for ch2o_b and ch3cho, the data read above are stp values.  these will
      ! be
      ! corrected in subgrid for t & ad dependence, after the altitude
      ! dependent
      ! values of t and ad are computed at each layer or level as appropriate.
      ! the o3->o(1d) yield is also calculated later, in subgrid.
      ! -----------------------------------------------------
      ! o2 cross section
      ! -----------------------------------------------------
      ! o3 -> o1d cross section
      ! -----------------------------------------------------
      ! o3 -> o3p cross section
      ! -----------------------------------------------------
      ! no2 cross section
      ! no2 quantum yield
      ! -----------------------------------------------------
      ! no3 -> no + o2 cross section
      ! no3 -> no + o2 quantum yield
      ! -----------------------------------------------------
      ! no3 -> no2 + o cross section
      ! no3 -> no2 + o quantum yield
      ! -----------------------------------------------------
      ! hono cross section
      ! hono cross quantum yield
      ! -----------------------------------------------------
      ! hno3 cross section
      ! hno3 cross quantum yield
      ! HNO3 CROSS SECTION TEMPERATURE DEPENDENCE
      ! -----------------------------------------------------
      ! hno4 cross section
      ! hno4 cross quantum yield
      ! -----------------------------------------------------
      ! h2o2 cross section
      ! h2o2 cross quantum yield
      ! -----------------------------------------------------
      ! hcho -> ho2 cross section
      ! hcho -> ho2 quantum yield
      ! -----------------------------------------------------
      ! hcho -> h2 cross section
      ! hcho -> h2 quantum yield
      ! -----------------------------------------------------
      ! ch3cho cross section
      ! ch3cho (ntp) quantum yield
      ! -----------------------------------------------------
      ! ch3coch3 cross section
      ! ch3coch3 quantum yield
      ! -----------------------------------------------------
      ! ch3coc2h5 cross section
      ! ch3coc2h5 quantum yield
      ! -----------------------------------------------------
      ! hcocho proc a cross section
      ! hcocho a quantum yield
      ! -----------------------------------------------------
      ! ch3cocho cross section
      ! ch3cocho quantum yield
      ! -----------------------------------------------------
      ! dcb cross section
      ! dcb quantum yield
      ! -----------------------------------------------------
      ! ch3o2h cross section
      ! ch3o2h cross quantum yield
      ! -----------------------------------------------------
      ! ch3coo2h cross section
      ! ch3coo2h cross quantum yield
      ! -----------------------------------------------------
      ! ch3ono2 cross section
      ! ch3ono2 cross quantum yield
      ! -----------------------------------------------------
      ! hcocho proc b cross section
      ! hcocho b quantum yield
      ! macr cross section
      ! macr quantum yield
      ! .. Parameters ..
      INTEGER, PARAMETER :: kl0 = 30, kl1 = 130
      INTEGER, PARAMETER :: kldif = (kl1-kl0+1)*3
      INTEGER, PARAMETER :: nabv = 10, nj = 200, nreakj = 23
      INTEGER, PARAMETER :: mj = 2*nj - 2
      ! .. Local Scalars ..
      INTEGER :: ip, kl
      ! .. Local Arrays ..
      REAL :: aerstd(51), airstd(51), albedoph(130), caabv(nabv), fext(130), &
        o3abv(nabv), o3std(51), pabv(nabv), so3tx(70,3), sra(11,9), srb(11,5), &
        tabv(nabv), tstd(51), txs(130,nreakj), wl(130), xqy(130,23), &
        xs(130,nreakj), zabv(nabv), zstd(51)
      ! .. Data Statements ..
      DATA zabv/21., 22., 23., 24., 25., 30., 35., 40., 45., 50./
      DATA tabv/215.19, 215.19, 215.19, 215.19, 215.19, 217.39, 227.80, &
        243.19, 258.50, 265.70/
      DATA pabv/1.57E18, 1.34E18, 1.14E18, 9.76E17, 8.33E17, 3.83E17, 1.76E17, &
        8.31E16, 4.09E16, 2.14E16/
      DATA o3abv/4.88E12, 4.86E12, 4.73E12, 4.54E12, 4.32E12, 2.52E12, &
        1.40E12, 6.07E11, 2.03E11, 6.61E10/
      DATA caabv/1.64E-3, 1.23E-3, 9.45E-4, 7.49E-4, 6.30E-4, 1.90E-4, &
        5.00E-5, 1.32E-5, 3.46E-6, 9.14E-7/
      DATA zstd/0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., &
        14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., &
        28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., &
        42., 43., 44., 45., 46., 47., 48., 49., 50./
      DATA ((sra(kl,ip),ip=1,9),kl=1,11)/ -2.158311E+01, -4.164652E-01, &
        5.266362E-02, 1.655877E-02, 0., 0., 0., 0., 0., -2.184813E+01, &
        -4.753880E-01, 4.519945E-02, 3.228313E-02, 3.079373E-03, 0., 0., 0., &
        0., -2.200507E+01, -4.628729E-01, -5.022541E-02, 2.545036E-02, &
        5.791406E-02, 1.179966E-02, -8.296876E-03, -3.238368E-03, &
        -3.069686E-04, -2.205527E+01, -4.400848E-01, -5.687308E-03, &
        3.712279E-02, 6.025527E-03, 0., 0., 0., 0., -2.205261E+01, &
        -5.707936E-01, -3.330207E-02, 5.959032E-02, 1.510540E-02, &
        1.000376E-03, 0., 0., 0., -2.228000E+01, -3.960759E-01, -2.995798E-02, &
        4.918104E-02, 9.269080E-03, -1.173411E-03, -2.599386E-04, 0., 0., &
        -2.275796E+01, -2.054719E-01, -1.094205E-02, 2.079595E-02, &
        3.769638E-03, 0., 0., 0., 0., -2.297610E+01, -5.823677E-02, &
        -1.007612E-01, 2.404666E-02, 4.761876E-02, 4.169606E-03, &
        -7.126663E-03, -2.263652E-03, -1.971653E-04, -2.506084E+01, &
        3.442774E-02, -2.212047E-04, 6.186041E-07, -6.284394E-10, 0., 0., 0., &
        0., -2.313436E+01, 1.177283E-04, 0., 0., 0., 0., 0., 0., 0., &
        -2.312205E+01, 0., 0., 0., 0., 0., 0., 0., 0./
      DATA ((srb(kl,ip),ip=1,5),kl=1,11)/ -2.431640E+03, 4.729722E+02, &
        -3.452121E+01, 1.120677E+00, -1.365618E-02, -3.701955E+01, &
        3.623290E+00, -8.929223E-02, 0., 0., -1.086239E+03, 1.981847E+02, &
        -1.359057E+01, 4.155845E-01, -4.788462E-03, -1.213108E+03, &
        2.277459E+02, -1.612207E+01, 5.101389E-01, -6.090518E-03, &
        -8.334575E+01, 7.944254E+00, -1.898894E-01, 0., 0., -2.139117E+02, &
        2.612729E+01, -1.036749E+00, 1.317695E-02, 0., -3.281301E+02, &
        4.307004E+01, -1.870019E+00, 2.674331E-02, 0., 3.033416E+03, &
        -5.978911E+02, 4.370384E+01, -1.406715E+00, 1.683967E-02, &
        -2.535815E+00, 0., 0., 0., 0., -4.474937E+00, 0., 0., 0., 0., &
        -2.996639E+00, 0., 0., 0., 0./
      DATA ((so3tx(kl,ip),kl=33,61),ip=1,3)/9.630E+00, 8.320E+00, 6.880E+00, &
        5.370E+00, 3.960E+00, 2.710E+00, 1.750E+00, 1.060E+00, 5.960E-01, &
        3.330E-01, 2.400E-01, 2.100E-01, 1.800E-01, 1.600E-01, 1.400E-01, &
        1.200E-01, 1.050E-01, 9.000E-02, 8.000E-02, 7.000E-02, 6.000E-02, &
        5.500E-02, 4.000E-02, 2.190E-02, 1.010E-02, 5.080E-03, 2.120E-03, &
        8.290E-04, 2.940E-04, 1.190E-03, 3.640E-04, 2.460E-04, 1.030E-03, &
        1.690E-03, 1.450E-03, 8.940E-04, 7.830E-04, 4.940E-04, 3.550E-04, &
        2.950E-04, 2.750E-04, 2.500E-04, 2.300E-04, 2.080E-04, 1.860E-04, &
        1.640E-04, 1.450E-04, 1.280E-04, 1.121E-04, 1.000E-04, 9.200E-05, &
        7.500E-05, 4.830E-05, 3.430E-05, 1.820E-05, 8.850E-06, 4.270E-06, &
        5.300E-06, -1.740E-05, 2.470E-06, 1.170E-05, 1.260E-06, -6.860E-06, &
        -2.890E-06, 3.590E-06, 2.000E-06, 3.660E-06, 2.600E-06, 2.170E-06, &
        1.950E-06, 1.380E-06, 1.650E-06, 1.550E-06, 1.460E-06, 1.340E-06, &
        1.210E-06, 1.130E-06, 1.060E-06, 9.400E-07, 8.700E-07, 7.500E-07, &
        5.200E-07, 2.660E-07, 1.630E-07, 1.260E-07, 8.710E-08, 3.500E-08/
      DATA aerstd/2.40E-1, 1.06E-1, 4.56E-2, 1.91E-2, 1.01E-2, 7.63E-3, &
        5.38E-3, 5.00E-3, 5.15E-3, 4.94E-3, 4.82E-3, 4.51E-3, 4.74E-3, &
        4.37E-3, 4.28E-3, 4.03E-3, 3.83E-3, 3.78E-3, 3.88E-3, 3.08E-3, &
        2.26E-3, 1.64E-3, 1.23E-3, 9.45E-4, 7.49E-4, 6.30E-4, 5.50E-4, &
        4.21E-4, 3.22E-4, 2.48E-4, 1.90E-4, 1.45E-4, 1.11E-4, 8.51E-5, &
        6.52E-5, 5.00E-5, 3.83E-5, 2.93E-5, 2.25E-5, 1.72E-5, 1.32E-5, &
        1.01E-5, 7.72E-6, 5.91E-6, 4.53E-6, 3.46E-6, 2.66E-6, 2.04E-6, &
        1.56E-6, 1.19E-6, 9.14E-7/
      DATA (wl(kl),kl=1,130)/1.861E+02, 1.878E+02, 1.896E+02, 1.914E+02, &
        1.933E+02, 1.952E+02, 1.971E+02, 1.990E+02, 2.010E+02, 2.031E+02, &
        2.052E+02, 2.073E+02, 2.094E+02, 2.117E+02, 2.139E+02, 2.162E+02, &
        2.186E+02, 2.210E+02, 2.235E+02, 2.260E+02, 2.286E+02, 2.313E+02, &
        2.340E+02, 2.367E+02, 2.396E+02, 2.425E+02, 2.454E+02, 2.485E+02, &
        2.516E+02, 2.548E+02, 2.582E+02, 2.615E+02, 2.650E+02, 2.685E+02, &
        2.722E+02, 2.759E+02, 2.798E+02, 2.837E+02, 2.878E+02, 2.920E+02, &
        2.963E+02, 3.005E+02, 3.030E+02, 3.040E+02, 3.050E+02, 3.060E+02, &
        3.070E+02, 3.080E+02, 3.090E+02, 3.100E+02, 3.110E+02, 3.120E+02, &
        3.130E+02, 3.140E+02, 3.160E+02, 3.200E+02, 3.250E+02, 3.300E+02, &
        3.350E+02, 3.400E+02, 3.450E+02, 3.500E+02, 3.550E+02, 3.600E+02, &
        3.650E+02, 3.700E+02, 3.750E+02, 3.800E+02, 3.850E+02, 3.900E+02, &
        3.950E+02, 4.000E+02, 4.050E+02, 4.100E+02, 4.150E+02, 4.200E+02, &
        4.250E+02, 4.300E+02, 4.350E+02, 4.400E+02, 4.450E+02, 4.500E+02, &
        4.550E+02, 4.600E+02, 4.650E+02, 4.700E+02, 4.750E+02, 4.800E+02, &
        4.850E+02, 4.900E+02, 4.950E+02, 5.000E+02, 5.050E+02, 5.100E+02, &
        5.150E+02, 5.200E+02, 5.250E+02, 5.300E+02, 5.350E+02, 5.400E+02, &
        5.450E+02, 5.500E+02, 5.550E+02, 5.600E+02, 5.650E+02, 5.700E+02, &
        5.750E+02, 5.800E+02, 5.850E+02, 5.900E+02, 5.950E+02, 6.000E+02, &
        6.050E+02, 6.100E+02, 6.150E+02, 6.200E+02, 6.250E+02, 6.300E+02, &
        6.350E+02, 6.400E+02, 6.448E+02, 6.510E+02, 6.600E+02, 6.700E+02, &
        6.800E+02, 6.900E+02, 7.000E+02, 7.100E+02, 7.200E+02, 7.300E+02/
      DATA (fext(kl),kl=1,130)/3.620E+11, 4.730E+11, 5.610E+11, 6.630E+11, &
        6.900E+11, 9.560E+11, 1.150E+12, 1.270E+12, 1.520E+12, 1.780E+12, &
        2.200E+12, 2.690E+12, 4.540E+12, 7.140E+12, 8.350E+12, 8.390E+12, &
        1.080E+13, 1.180E+13, 1.600E+13, 1.340E+13, 1.410E+13, 1.570E+13, &
        1.380E+13, 1.600E+13, 1.450E+13, 2.200E+13, 1.990E+13, 1.970E+13, &
        1.940E+13, 2.910E+13, 4.950E+13, 4.530E+13, 1.070E+14, 1.200E+14, &
        1.100E+14, 1.040E+14, 8.240E+13, 1.520E+14, 2.150E+14, 3.480E+14, &
        3.396E+14, 2.730E+14, 9.109E+13, 8.745E+13, 9.577E+13, 8.507E+13, &
        9.383E+13, 1.030E+14, 9.722E+13, 7.751E+13, 1.277E+14, 1.087E+14, &
        1.102E+14, 1.184E+14, 3.153E+14, 5.930E+14, 6.950E+14, 8.150E+14, &
        7.810E+14, 8.350E+14, 8.140E+14, 8.530E+14, 9.170E+14, 8.380E+14, &
        1.040E+15, 1.100E+15, 9.790E+14, 1.130E+15, 8.890E+14, 1.140E+15, &
        9.170E+14, 1.690E+15, 1.700E+15, 1.840E+15, 1.870E+15, 1.950E+15, &
        1.810E+15, 1.670E+15, 1.980E+15, 2.020E+15, 2.180E+15, 2.360E+15, &
        2.310E+15, 2.390E+15, 2.380E+15, 2.390E+15, 2.440E+15, 2.510E+15, &
        2.300E+15, 2.390E+15, 2.480E+15, 2.400E+15, 2.460E+15, 2.490E+15, &
        2.320E+15, 2.390E+15, 2.420E+15, 2.550E+15, 2.510E+15, 2.490E+15, &
        2.550E+15, 2.530E+15, 2.540E+15, 2.500E+15, 2.570E+15, 2.580E+15, &
        2.670E+15, 2.670E+15, 2.700E+15, 2.620E+15, 2.690E+15, 2.630E+15, &
        2.680E+15, 2.660E+15, 2.590E+15, 2.690E+15, 2.610E+15, 2.620E+15, &
        2.620E+15, 2.630E+15, 2.392E+15, 3.998E+15, 5.115E+15, 5.225E+15, &
        5.215E+15, 5.105E+15, 5.140E+15, 5.010E+15, 4.930E+15, 4.895E+15/
      DATA (xs(kl,1),kl=1,130)/7.040E-24, 7.360E-24, 7.640E-24, 7.870E-24, &
        8.040E-24, 8.140E-24, 8.170E-24, 8.130E-24, 8.010E-24, 7.840E-24, &
        7.630E-24, 7.330E-24, 6.990E-24, 6.450E-24, 5.810E-24, 5.230E-24, &
        4.710E-24, 4.260E-24, 3.800E-24, 3.350E-24, 2.900E-24, 2.450E-24, &
        2.050E-24, 1.690E-24, 1.300E-24, 9.300E-25, 0., 0., 0., 0., 0., 0., &
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0./
      DATA (xs(kl,2),kl=1,130)/0.622E-18, 0.576E-18, 0.526E-18, 0.476E-18, &
        0.428E-18, 0.383E-18, 0.347E-18, 0.323E-18, 0.314E-18, 0.326E-18, &
        0.364E-18, 0.434E-18, 0.542E-18, 0.699E-18, 0.921E-18, 0.119E-17, &
        0.155E-17, 0.199E-17, 0.256E-17, 0.323E-17, 0.400E-17, 0.483E-17, &
        0.579E-17, 0.686E-17, 0.797E-17, 0.900E-17, 0.100E-16, 0.108E-16, &
        0.113E-16, 0.115E-16, 0.112E-16, 0.106E-16, 0.963E-17, 0.836E-17, &
        0.695E-17, 0.545E-17, 0.404E-17, 0.280E-17, 0.183E-17, 0.112E-17, &
        0.647E-18, 0.369E-18, 0.270E-18, 0.238E-18, 0.203E-18, 0.183E-18, &
        0.161E-18, 0.139E-18, 0.122E-18, 0.105E-18, 0.939E-19, 0.825E-19, &
        0.711E-19, 0.653E-19, 0.486E-19, 0.276E-19, 0.137E-19, 0.707E-20, &
        0.330E-20, 0.152E-20, 0.816E-21, 0.266E-21, 0.109E-21, 0.549E-22, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.291E-22, 0.314E-22, 0.399E-22, &
        0.654E-22, 0.683E-22, 0.866E-22, 0.125E-21, 0.149E-21, 0.171E-21, &
        0.212E-21, 0.357E-21, 0.368E-21, 0.406E-21, 0.489E-21, 0.711E-21, &
        0.843E-21, 0.828E-21, 0.909E-21, 0.122E-20, 0.162E-20, 0.158E-20, &
        0.160E-20, 0.178E-20, 0.207E-20, 0.255E-20, 0.274E-20, 0.288E-20, &
        0.307E-20, 0.317E-20, 0.336E-20, 0.388E-20, 0.431E-20, 0.467E-20, &
        0.475E-20, 0.455E-20, 0.435E-20, 0.442E-20, 0.461E-20, 0.489E-20, &
        0.484E-20, 0.454E-20, 0.424E-20, 0.390E-20, 0.360E-20, 0.343E-20, &
        0.317E-20, 0.274E-20, 0.261E-20, 0.240E-20, 0.207E-20, 0.172E-20, &
        0.137E-20, 0.111E-20, 0.913E-21, 0.793E-21, 0.640E-21, 0.514E-21/
      DATA (xs(kl,3),kl=1,130)/0.622E-18, 0.576E-18, 0.526E-18, 0.476E-18, &
        0.428E-18, 0.383E-18, 0.347E-18, 0.323E-18, 0.314E-18, 0.326E-18, &
        0.364E-18, 0.434E-18, 0.542E-18, 0.699E-18, 0.921E-18, 0.119E-17, &
        0.155E-17, 0.199E-17, 0.256E-17, 0.323E-17, 0.400E-17, 0.483E-17, &
        0.579E-17, 0.686E-17, 0.797E-17, 0.900E-17, 0.100E-16, 0.108E-16, &
        0.113E-16, 0.115E-16, 0.112E-16, 0.106E-16, 0.963E-17, 0.836E-17, &
        0.695E-17, 0.545E-17, 0.404E-17, 0.280E-17, 0.183E-17, 0.112E-17, &
        0.647E-18, 0.369E-18, 0.270E-18, 0.238E-18, 0.203E-18, 0.183E-18, &
        0.161E-18, 0.139E-18, 0.122E-18, 0.105E-18, 0.939E-19, 0.825E-19, &
        0.711E-19, 0.653E-19, 0.486E-19, 0.276E-19, 0.137E-19, 0.707E-20, &
        0.330E-20, 0.152E-20, 0.816E-21, 0.266E-21, 0.109E-21, 0.549E-22, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.291E-22, 0.314E-22, 0.399E-22, &
        0.654E-22, 0.683E-22, 0.866E-22, 0.125E-21, 0.149E-21, 0.171E-21, &
        0.212E-21, 0.357E-21, 0.368E-21, 0.406E-21, 0.489E-21, 0.711E-21, &
        0.843E-21, 0.828E-21, 0.909E-21, 0.122E-20, 0.162E-20, 0.158E-20, &
        0.160E-20, 0.178E-20, 0.207E-20, 0.255E-20, 0.274E-20, 0.288E-20, &
        0.307E-20, 0.317E-20, 0.336E-20, 0.388E-20, 0.431E-20, 0.467E-20, &
        0.475E-20, 0.455E-20, 0.435E-20, 0.442E-20, 0.461E-20, 0.489E-20, &
        0.484E-20, 0.454E-20, 0.424E-20, 0.390E-20, 0.360E-20, 0.343E-20, &
        0.317E-20, 0.274E-20, 0.261E-20, 0.240E-20, 0.207E-20, 0.172E-20, &
        0.137E-20, 0.111E-20, 0.913E-21, 0.793E-21, 0.640E-21, 0.514E-21/
      DATA (xs(kl,4),kl=1,130)/0.259E-18, 0.272E-18, 0.286E-18, 0.273E-18, &
        0.251E-18, 0.244E-18, 0.246E-18, 0.246E-18, 0.282E-18, 0.415E-18, &
        0.448E-18, 0.445E-18, 0.464E-18, 0.487E-18, 0.482E-18, 0.502E-18, &
        0.444E-18, 0.471E-18, 0.377E-18, 0.393E-18, 0.274E-18, 0.278E-18, &
        0.169E-18, 0.162E-18, 0.882E-19, 0.747E-19, 0.391E-19, 0.275E-19, &
        0.201E-19, 0.197E-19, 0.211E-19, 0.236E-19, 0.270E-19, 0.325E-19, &
        0.379E-19, 0.503E-19, 0.588E-19, 0.700E-19, 0.815E-19, 0.972E-19, &
        0.115E-18, 0.128E-18, 0.154E-18, 0.159E-18, 0.158E-18, 0.156E-18, &
        0.164E-18, 0.166E-18, 0.182E-18, 0.184E-18, 0.192E-18, 0.204E-18, &
        0.204E-18, 0.202E-18, 0.224E-18, 0.248E-18, 0.281E-18, 0.313E-18, &
        0.343E-18, 0.380E-18, 0.407E-18, 0.431E-18, 0.472E-18, 0.483E-18, &
        0.517E-18, 0.532E-18, 0.551E-18, 0.564E-18, 0.576E-18, 0.593E-18, &
        0.585E-18, 0.602E-18, 0.578E-18, 0.600E-18, 0.565E-18, 0.581E-18, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA ((xqy(kl,ip),kl=kl0,kl1),ip=1,3)/kldif*1./
      DATA (xqy(kl,4),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        0.999000, 0.998000, 0.997000, 0.996500, 0.996000, 0.996000, 0.996000, &
        0.996000, 0.995000, 0.995000, 0.995000, 0.995000, 0.995000, 0.994000, &
        0.994000, 0.994000, 0.993000, 0.992000, 0.991000, 0.990000, 0.989000, &
        0.988000, 0.987000, 0.986000, 0.984000, 0.983000, 0.981000, 0.979000, &
        0.975000, 0.969000, 0.960000, 0.927000, 0.694000, 0.355000, 0.134000, &
        0.060000, 0.018000, 0.000900, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,5),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.160E-19, 0.240E-19, 0.520E-19, 0.760E-19, &
        0.110E-18, 0.136E-18, 0.170E-18, 0.198E-18, 0.220E-18, 0.288E-18, &
        0.358E-18, 0.396E-18, 0.502E-18, 0.598E-18, 0.694E-18, 0.716E-18, &
        0.828E-18, 0.984E-18, 0.110E-17, 0.114E-17, 0.125E-17, 0.153E-17, &
        0.156E-17, 0.168E-17, 0.169E-17, 0.217E-17, 0.229E-17, 0.208E-17, &
        0.213E-17, 0.261E-17, 0.299E-17, 0.329E-17, 0.278E-17, 0.281E-17, &
        0.307E-17, 0.334E-17, 0.322E-17, 0.554E-17, 0.441E-17, 0.314E-17, &
        0.365E-17, 0.188E-17, 0.233E-17, 0.473E-17, 0.100E-16, 0.584E-17, &
        0.180E-17, 0.135E-17, 0.822E-18, 0.640E-18, 0.777E-17, 0.134E-17, &
        0.337E-18, 0.175E-19, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,5),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.038000, &
        0.191000, 0.326000, 0.311000, 0.272000, 0.233000, 0.194000, 0.156000, &
        0.117000, 0.078000, 0.039000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,6),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.160E-19, 0.240E-19, 0.520E-19, 0.760E-19, &
        0.110E-18, 0.136E-18, 0.170E-18, 0.198E-18, 0.220E-18, 0.288E-18, &
        0.358E-18, 0.396E-18, 0.502E-18, 0.598E-18, 0.694E-18, 0.716E-18, &
        0.828E-18, 0.984E-18, 0.110E-17, 0.114E-17, 0.125E-17, 0.153E-17, &
        0.156E-17, 0.168E-17, 0.169E-17, 0.217E-17, 0.229E-17, 0.208E-17, &
        0.213E-17, 0.261E-17, 0.299E-17, 0.329E-17, 0.278E-17, 0.281E-17, &
        0.307E-17, 0.334E-17, 0.322E-17, 0.554E-17, 0.441E-17, 0.314E-17, &
        0.365E-17, 0.188E-17, 0.233E-17, 0.473E-17, 0.100E-16, 0.584E-17, &
        0.180E-17, 0.135E-17, 0.822E-18, 0.640E-18, 0.777E-17, 0.134E-17, &
        0.337E-18, 0.175E-19, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,6),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.962000, &
        0.809000, 0.661000, 0.578000, 0.506000, 0.433000, 0.361000, 0.289000, &
        0.217000, 0.144000, 0.072000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,7),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.130E-19, 0.190E-19, 0.280E-19, &
        0.220E-19, 0.360E-19, 0.340E-19, 0.536E-19, 0.534E-19, 0.111E-18, &
        0.786E-19, 0.189E-18, 0.116E-18, 0.130E-18, 0.279E-18, 0.954E-19, &
        0.179E-18, 0.260E-18, 0.590E-19, 0.101E-18, 0.176E-18, 0.304E-19, &
        0.775E-20, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,7),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
      DATA (xs(kl,8),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.127E-16, &
        0.114E-16, 0.100E-16, 0.847E-17, 0.679E-17, 0.518E-17, 0.382E-17, &
        0.270E-17, 0.182E-17, 0.120E-17, 0.730E-18, 0.451E-18, 0.283E-18, &
        0.195E-18, 0.134E-18, 0.102E-18, 0.802E-19, 0.650E-19, 0.518E-19, &
        0.414E-19, 0.321E-19, 0.265E-19, 0.230E-19, 0.209E-19, 0.199E-19, &
        0.196E-19, 0.195E-19, 0.193E-19, 0.188E-19, 0.180E-19, 0.170E-19, &
        0.152E-19, 0.134E-19, 0.113E-19, 0.924E-20, 0.719E-20, 0.532E-20, &
        0.371E-20, 0.249E-20, 0.188E-20, 0.167E-20, 0.150E-20, 0.133E-20, &
        0.119E-20, 0.105E-20, 0.932E-21, 0.814E-21, 0.721E-21, 0.628E-21, &
        0.547E-21, 0.465E-21, 0.362E-21, 0.197E-21, 0.975E-22, 0.452E-22, &
        0.222E-22, 0.110E-22, 0.604E-23, 0.420E-23, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,8),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
      DATA (txs(kl,8),kl=1,130)/0.000000, 0.000000, 0.000000, 1.700000, &
        1.650000, 1.660000, 1.690000, 1.740000, 1.770000, 1.850000, 1.970000, &
        2.080000, 2.170000, 2.170000, 2.210000, 2.150000, 2.060000, 1.960000, &
        1.840000, 1.780000, 1.800000, 1.860000, 1.900000, 1.970000, 1.970000, &
        1.970000, 1.880000, 1.750000, 1.610000, 1.440000, 1.340000, 1.230000, &
        1.180000, 1.140000, 1.120000, 1.140000, 1.140000, 1.180000, 1.220000, &
        1.250000, 1.450000, 1.490000, 1.560000, 1.640000, 1.690000, 1.780000, &
        1.870000, 1.940000, 2.040000, 2.150000, 2.270000, 2.380000, 2.620000, &
        2.700000, 2.920000, 3.100000, 3.240000, 3.520000, 3.770000, 3.910000, &
        4.230000, 4.700000, 5.150000, 5.250000, 5.740000, 6.450000, 6.700000, &
        7.160000, 7.550000, 8.160000, 9.750000, 9.930000, 9.600000, 10.50000, &
        10.80000, 11.80000, 11.80000, 9.300000, 12.10000, 11.90000, 9.300000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,9),kl=1,130)/0.000E+00, 0.000E+00, 0.100E-16, 0.980E-17, &
        0.918E-17, 0.828E-17, 0.723E-17, 0.618E-17, 0.517E-17, 0.426E-17, &
        0.352E-17, 0.292E-17, 0.245E-17, 0.205E-17, 0.175E-17, 0.151E-17, &
        0.131E-17, 0.115E-17, 0.102E-17, 0.916E-18, 0.827E-18, 0.752E-18, &
        0.687E-18, 0.631E-18, 0.578E-18, 0.529E-18, 0.484E-18, 0.439E-18, &
        0.396E-18, 0.353E-18, 0.311E-18, 0.271E-18, 0.231E-18, 0.194E-18, &
        0.158E-18, 0.125E-18, 0.946E-19, 0.694E-19, 0.485E-19, 0.325E-19, &
        0.210E-19, 0.135E-19, 0.104E-19, 0.938E-20, 0.846E-20, 0.763E-20, &
        0.689E-20, 0.623E-20, 0.564E-20, 0.511E-20, 0.463E-20, 0.420E-20, &
        0.382E-20, 0.347E-20, 0.287E-20, 0.191E-20, 0.101E-20, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,9),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
      DATA (xs(kl,10),kl=1,130)/0.000E+00, 0.000E+00, 0.695E-18, 0.635E-18, &
        0.586E-18, 0.547E-18, 0.515E-18, 0.488E-18, 0.462E-18, 0.436E-18, &
        0.412E-18, 0.388E-18, 0.365E-18, 0.341E-18, 0.318E-18, 0.295E-18, &
        0.272E-18, 0.250E-18, 0.229E-18, 0.209E-18, 0.190E-18, 0.172E-18, &
        0.155E-18, 0.140E-18, 0.126E-18, 0.112E-18, 0.999E-19, 0.881E-19, &
        0.774E-19, 0.675E-19, 0.582E-19, 0.500E-19, 0.425E-19, 0.358E-19, &
        0.298E-19, 0.247E-19, 0.201E-19, 0.164E-19, 0.131E-19, 0.105E-19, &
        0.828E-20, 0.658E-20, 0.573E-20, 0.543E-20, 0.514E-20, 0.486E-20, &
        0.460E-20, 0.435E-20, 0.412E-20, 0.390E-20, 0.369E-20, 0.349E-20, &
        0.330E-20, 0.312E-20, 0.279E-20, 0.220E-20, 0.160E-20, 0.130E-20, &
        0.100E-20, 0.700E-21, 0.500E-21, 0.400E-21, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,10),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
      DATA (xs(kl,11),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.293E-21, 0.342E-21, 0.102E-20, 0.456E-21, 0.527E-21, 0.537E-21, &
        0.347E-21, 0.759E-21, 0.628E-21, 0.974E-21, 0.104E-20, 0.219E-20, &
        0.228E-20, 0.357E-20, 0.374E-20, 0.584E-20, 0.651E-20, 0.102E-19, &
        0.114E-19, 0.176E-19, 0.180E-19, 0.259E-19, 0.227E-19, 0.275E-19, &
        0.318E-19, 0.160E-19, 0.245E-19, 0.637E-19, 0.426E-19, 0.399E-19, &
        0.186E-19, 0.131E-19, 0.310E-19, 0.182E-19, 0.596E-20, 0.111E-19, &
        0.911E-20, 0.457E-19, 0.423E-19, 0.142E-19, 0.243E-19, 0.178E-19, &
        0.129E-20, 0.213E-19, 0.661E-20, 0.139E-20, 0.827E-20, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,11),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.264091, &
        0.288940, 0.297038, 0.295059, 0.289384, 0.285483, 0.288087, 0.301000, &
        0.326819, 0.366764, 0.420506, 0.485961, 0.559106, 0.633887, 0.702103, &
        0.733457, 0.740762, 0.747845, 0.753000, 0.754000, 0.754800, 0.754000, &
        0.753000, 0.752000, 0.751000, 0.749500, 0.745000, 0.739600, 0.731700, &
        0.723300, 0.690300, 0.593100, 0.458100, 0.305000, 0.122300, 0.003429, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,12),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.293E-21, 0.342E-21, 0.102E-20, 0.456E-21, 0.527E-21, 0.537E-21, &
        0.347E-21, 0.759E-21, 0.628E-21, 0.974E-21, 0.104E-20, 0.219E-20, &
        0.228E-20, 0.357E-20, 0.374E-20, 0.584E-20, 0.651E-20, 0.102E-19, &
        0.114E-19, 0.176E-19, 0.180E-19, 0.259E-19, 0.227E-19, 0.275E-19, &
        0.318E-19, 0.160E-19, 0.245E-19, 0.637E-19, 0.426E-19, 0.399E-19, &
        0.186E-19, 0.131E-19, 0.310E-19, 0.182E-19, 0.596E-20, 0.111E-19, &
        0.911E-20, 0.457E-19, 0.423E-19, 0.142E-19, 0.243E-19, 0.178E-19, &
        0.129E-20, 0.213E-19, 0.661E-20, 0.139E-20, 0.827E-20, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,12),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.492085, &
        0.483681, 0.483325, 0.487471, 0.492514, 0.495532, 0.493893, 0.485473, &
        0.468839, 0.443373, 0.409405, 0.368400, 0.323132, 0.307820, 0.294564, &
        0.280920, 0.266885, 0.253277, 0.249000, 0.247000, 0.245600, 0.248000, &
        0.251000, 0.254000, 0.257000, 0.260200, 0.264500, 0.269000, 0.273500, &
        0.278900, 0.310300, 0.394100, 0.508100, 0.676100, 0.759300, 0.636100, &
        0.501500, 0.373400, 0.229000, 0.103600, 0.005906, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,13),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.593E-21, 0.548E-21, &
        0.552E-21, 0.513E-21, 0.480E-21, 0.482E-21, 0.482E-21, 0.482E-21, &
        0.536E-21, 0.593E-21, 0.734E-21, 0.948E-21, 0.125E-20, 0.171E-20, &
        0.234E-20, 0.321E-20, 0.434E-20, 0.582E-20, 0.770E-20, 0.991E-20, &
        0.127E-19, 0.159E-19, 0.200E-19, 0.237E-19, 0.287E-19, 0.326E-19, &
        0.376E-19, 0.408E-19, 0.444E-19, 0.463E-19, 0.466E-19, 0.465E-19, &
        0.432E-19, 0.406E-19, 0.372E-19, 0.348E-19, 0.342E-19, 0.342E-19, &
        0.336E-19, 0.333E-19, 0.314E-19, 0.293E-19, 0.276E-19, 0.253E-19, &
        0.247E-19, 0.243E-19, 0.210E-19, 0.169E-19, 0.108E-19, 0.651E-20, &
        0.314E-20, 0.138E-20, 0.224E-21, 0.947E-22, 0.425E-22, 0.229E-22, &
        0.275E-23, 0.192E-23, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,13),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.310000, &
        0.349300, 0.377900, 0.430300, 0.501600, 0.566200, 0.561500, 0.541100, &
        0.512100, 0.473200, 0.430000, 0.392000, 0.376000, 0.360000, 0.344000, &
        0.328000, 0.312000, 0.296000, 0.279800, 0.262500, 0.245000, 0.227500, &
        0.210000, 0.175000, 0.109300, 0.051880, 0.006266, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,14),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.003E-20, 1.873E-21, &
        1.408E-21, 1.084E-21, 1.030E-21, 1.071E-21, 1.176E-21, 1.371E-21, &
        1.690E-21, 2.153E-21, 2.779E-21, 3.561E-21, 4.575E-21, 5.906E-21, &
        7.585E-21, 9.622E-21, 1.212E-20, 1.516E-20, 1.863E-20, 2.252E-20, &
        2.676E-20, 3.143E-20, 3.616E-20, 4.058E-20, 4.475E-20, 4.774E-20, &
        5.029E-20, 5.065E-20, 5.049E-20, 4.790E-20, 4.415E-20, 3.940E-20, &
        3.308E-20, 2.717E-20, 2.371E-20, 2.244E-20, 2.106E-20, 1.953E-20, &
        1.801E-20, 1.663E-20, 1.538E-20, 1.408E-20, 1.277E-20, 1.173E-20, &
        1.081E-20, 9.675E-21, 7.783E-21, 4.712E-21, 2.078E-21, 7.065E-22, &
        1.973E-22, 5.745E-23, 2.137E-23, 8.235E-24, 5.686E-24, 2.157E-24, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,14),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.766400, 0.779200, 0.792800, 0.776000, &
        0.720000, 0.672000, 0.620200, 0.586900, 0.551800, 0.457500, 0.355000, &
        0.270000, 0.205500, 0.145000, 0.120000, 0.110000, 0.100000, 0.090000, &
        0.080000, 0.070000, 0.060000, 0.050000, 0.047800, 0.045600, 0.043400, &
        0.041200, 0.036800, 0.028000, 0.030500, 0.033000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,15),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.234E-19, 0.685E-20, &
        0.215E-20, 0.168E-20, 0.159E-20, 0.164E-20, 0.181E-20, 0.203E-20, &
        0.230E-20, 0.268E-20, 0.320E-20, 0.387E-20, 0.471E-20, 0.582E-20, &
        0.728E-20, 0.913E-20, 0.115E-19, 0.145E-19, 0.180E-19, 0.222E-19, &
        0.267E-19, 0.321E-19, 0.374E-19, 0.430E-19, 0.479E-19, 0.525E-19, &
        0.555E-19, 0.576E-19, 0.575E-19, 0.555E-19, 0.518E-19, 0.460E-19, &
        0.388E-19, 0.319E-19, 0.269E-19, 0.251E-19, 0.233E-19, 0.217E-19, &
        0.202E-19, 0.188E-19, 0.173E-19, 0.158E-19, 0.142E-19, 0.128E-19, &
        0.114E-19, 0.101E-19, 0.796E-20, 0.463E-20, 0.196E-20, 0.705E-21, &
        0.207E-21, 0.545E-22, 0.145E-22, 0.431E-23, 0.471E-23, 0.157E-23, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,15),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.766400, 0.779200, 0.792800, 0.776000, &
        0.720000, 0.672000, 0.620200, 0.586900, 0.551800, 0.457500, 0.355000, &
        0.270000, 0.205500, 0.145000, 0.120000, 0.110000, 0.100000, 0.090000, &
        0.080000, 0.070000, 0.060000, 0.050000, 0.047800, 0.045600, 0.043400, &
        0.041200, 0.036800, 0.028000, 0.030500, 0.033000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,16),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.340E-20, 0.401E-20, 0.486E-20, 0.573E-20, 0.662E-20, 0.754E-20, &
        0.906E-20, 0.112E-19, 0.133E-19, 0.162E-19, 0.202E-19, 0.224E-19, &
        0.248E-19, 0.276E-19, 0.289E-19, 0.318E-19, 0.322E-19, 0.321E-19, &
        0.338E-19, 0.343E-19, 0.307E-19, 0.290E-19, 0.275E-19, 0.272E-19, &
        0.272E-19, 0.272E-19, 0.272E-19, 0.273E-19, 0.280E-19, 0.283E-19, &
        0.268E-19, 0.249E-19, 0.212E-19, 0.151E-19, 0.127E-19, 0.142E-19, &
        0.229E-19, 0.358E-20, 0.000E+00, 0.000E+00, 0.286E-21, 0.208E-20, &
        0.344E-20, 0.764E-20, 0.107E-19, 0.159E-19, 0.166E-19, 0.303E-19, &
        0.263E-19, 0.336E-19, 0.366E-19, 0.456E-19, 0.643E-19, 0.546E-19, &
        0.922E-19, 0.677E-19, 0.599E-19, 0.117E-18, 0.715E-19, 0.730E-19, &
        0.201E-18, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,16),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
        0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000/
      DATA (xs(kl,17),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.131E-19, 0.142E-19, 0.156E-19, &
        0.174E-19, 0.189E-19, 0.205E-19, 0.219E-19, 0.233E-19, 0.252E-19, &
        0.269E-19, 0.285E-19, 0.313E-19, 0.338E-19, 0.362E-19, 0.394E-19, &
        0.427E-19, 0.450E-19, 0.486E-19, 0.476E-19, 0.479E-19, 0.465E-19, &
        0.420E-19, 0.371E-19, 0.352E-19, 0.344E-19, 0.336E-19, 0.316E-19, &
        0.296E-19, 0.276E-19, 0.256E-19, 0.237E-19, 0.227E-19, 0.218E-19, &
        0.208E-19, 0.199E-19, 0.182E-19, 0.151E-19, 0.938E-20, 0.652E-20, &
        0.482E-20, 0.323E-20, 0.300E-20, 0.394E-20, 0.560E-20, 0.695E-20, &
        0.108E-19, 0.148E-19, 0.191E-19, 0.243E-19, 0.322E-19, 0.403E-19, &
        0.473E-19, 0.566E-19, 0.692E-19, 0.846E-19, 0.968E-19, 0.103E-18, &
        0.102E-18, 0.101E-18, 0.106E-18, 0.104E-18, 0.994E-19, 0.813E-19, &
        0.395E-19, 0.109E-19, 0.327E-20, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,17),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 0.990000, 0.990000, 0.990000, &
        0.980000, 0.980000, 0.980000, 0.980000, 0.970000, 0.970000, 0.970000, &
        0.960000, 0.960000, 0.940000, 0.920000, 0.880000, 0.825000, 0.750000, &
        0.660000, 0.560000, 0.480000, 0.400000, 0.320000, 0.250000, 0.200000, &
        0.150000, 0.120000, 0.100000, 0.080000, 0.060000, 0.050000, 0.040000, &
        0.030000, 0.020000, 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, &
        0.005000, 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,18),kl=1,130)/0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,18),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 0.500000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,19),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.320E-18, 0.268E-18, 0.226E-18, 0.193E-18, &
        0.167E-18, 0.147E-18, 0.129E-18, 0.115E-18, 0.102E-18, 0.899E-19, &
        0.797E-19, 0.708E-19, 0.623E-19, 0.548E-19, 0.483E-19, 0.422E-19, &
        0.369E-19, 0.321E-19, 0.278E-19, 0.242E-19, 0.209E-19, 0.180E-19, &
        0.154E-19, 0.131E-19, 0.111E-19, 0.925E-20, 0.763E-20, 0.622E-20, &
        0.501E-20, 0.402E-20, 0.352E-20, 0.333E-20, 0.316E-20, 0.299E-20, &
        0.283E-20, 0.268E-20, 0.254E-20, 0.240E-20, 0.227E-20, 0.215E-20, &
        0.204E-20, 0.193E-20, 0.172E-20, 0.138E-20, 0.105E-20, 0.801E-21, &
        0.612E-21, 0.467E-21, 0.356E-21, 0.270E-21, 0.206E-21, 0.160E-21, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,19),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
      DATA (xs(kl,20),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.898E-19, &
        0.841E-19, 0.784E-19, 0.148E-18, 0.138E-18, 0.129E-18, 0.122E-18, &
        0.114E-18, 0.107E-18, 0.998E-19, 0.932E-19, 0.868E-19, 0.807E-19, &
        0.748E-19, 0.690E-19, 0.633E-19, 0.579E-19, 0.529E-19, 0.480E-19, &
        0.433E-19, 0.390E-19, 0.349E-19, 0.313E-19, 0.278E-19, 0.247E-19, &
        0.217E-19, 0.190E-19, 0.162E-19, 0.138E-19, 0.118E-19, 0.991E-20, &
        0.829E-20, 0.690E-20, 0.566E-20, 0.452E-20, 0.361E-20, 0.288E-20, &
        0.229E-20, 0.181E-20, 0.157E-20, 0.147E-20, 0.138E-20, 0.131E-20, &
        0.125E-20, 0.118E-20, 0.112E-20, 0.106E-20, 0.100E-20, 0.947E-21, &
        0.893E-21, 0.838E-21, 0.743E-21, 0.581E-21, 0.428E-21, 0.332E-21, &
        0.262E-21, 0.192E-21, 0.141E-21, 0.910E-22, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,20),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
      DATA (xs(kl,21),kl=1,130)/0.180E-16, 0.181E-16, 0.179E-16, 0.174E-16, &
        0.167E-16, 0.159E-16, 0.146E-16, 0.133E-16, 0.118E-16, 0.101E-16, &
        0.850E-17, 0.695E-17, 0.540E-17, 0.411E-17, 0.301E-17, 0.215E-17, &
        0.163E-17, 0.105E-17, 0.754E-18, 0.524E-18, 0.382E-18, 0.272E-18, &
        0.202E-18, 0.152E-18, 0.112E-18, 0.876E-19, 0.677E-19, 0.596E-19, &
        0.541E-19, 0.510E-19, 0.489E-19, 0.469E-19, 0.448E-19, 0.419E-19, &
        0.377E-19, 0.336E-19, 0.285E-19, 0.238E-19, 0.190E-19, 0.145E-19, &
        0.106E-19, 0.752E-20, 0.610E-20, 0.553E-20, 0.496E-20, 0.457E-20, &
        0.418E-20, 0.380E-20, 0.341E-20, 0.302E-20, 0.279E-20, 0.256E-20, &
        0.232E-20, 0.209E-20, 0.171E-20, 0.110E-20, 0.644E-21, 0.416E-21, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,21),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
        1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
      DATA (xs(kl,22),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.340E-20, 0.401E-20, 0.486E-20, 0.573E-20, 0.662E-20, 0.754E-20, &
        0.906E-20, 0.112E-19, 0.133E-19, 0.162E-19, 0.202E-19, 0.224E-19, &
        0.248E-19, 0.276E-19, 0.289E-19, 0.318E-19, 0.322E-19, 0.321E-19, &
        0.338E-19, 0.343E-19, 0.307E-19, 0.290E-19, 0.275E-19, 0.272E-19, &
        0.272E-19, 0.272E-19, 0.272E-19, 0.273E-19, 0.280E-19, 0.283E-19, &
        0.268E-19, 0.249E-19, 0.212E-19, 0.151E-19, 0.127E-19, 0.142E-19, &
        0.229E-19, 0.358E-20, 0.000E+00, 0.000E+00, 0.286E-21, 0.208E-20, &
        0.344E-20, 0.764E-20, 0.107E-19, 0.159E-19, 0.166E-19, 0.303E-19, &
        0.263E-19, 0.336E-19, 0.366E-19, 0.456E-19, 0.643E-19, 0.546E-19, &
        0.922E-19, 0.677E-19, 0.599E-19, 0.117E-18, 0.715E-19, 0.730E-19, &
        0.201E-18, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,22),kl=1,130)/0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      DATA (xs(kl,23),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.164E-20, &
        0.211E-20, 0.227E-20, 0.264E-20, 0.340E-20, 0.465E-20, 0.598E-20, &
        0.803E-20, 0.986E-20, 0.118E-19, 0.137E-19, 0.160E-19, 0.189E-19, &
        0.228E-19, 0.275E-19, 0.307E-19, 0.321E-19, 0.335E-19, 0.349E-19, &
        0.363E-19, 0.378E-19, 0.393E-19, 0.407E-19, 0.421E-19, 0.435E-19, &
        0.449E-19, 0.462E-19, 0.487E-19, 0.527E-19, 0.564E-19, 0.589E-19, &
        0.616E-19, 0.556E-19, 0.553E-19, 0.543E-19, 0.365E-19, 0.318E-19, &
        0.316E-19, 0.156E-19, 0.428E-20, 0.113E-20, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
        0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
      DATA (xqy(kl,23),kl=1,130)/0.019894, 0.019716, 0.019528, 0.019339, &
        0.019140, 0.018941, 0.018742, 0.018543, 0.018333, 0.018113, 0.017893, &
        0.017673, 0.017453, 0.017212, 0.016982, 0.016741, 0.016531, 0.016238, &
        0.015976, 0.015714, 0.015442, 0.015159, 0.014876, 0.014593, 0.014290, &
        0.013986, 0.013682, 0.013357, 0.013032, 0.012697, 0.012341, 0.011995, &
        0.011629, 0.011314, 0.010874, 0.010487, 0.010078, 0.009670, 0.009240, &
        0.008800, 0.008350, 0.007910, 0.007648, 0.007543, 0.007438, 0.007333, &
        0.007229, 0.007124, 0.007019, 0.006914, 0.006810, 0.006705, 0.006600, &
        0.006495, 0.006286, 0.005867, 0.005343, 0.004819, 0.004295, 0.003772, &
        0.003248, 0.002724, 0.002200, 0.001676, 0.001153, 0.000629, 0.000105, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
      ! END MODULE module_data_photmad

    CONTAINS
      subroutine madronich1_driver(id,curr_secs,ktau,config_flags,haveaer,&
               gmt,julday,t_phy,moist,aerwrf,p8w,t8w,p_phy,               &
               chem,rho_phy,dz8w,                                         &
               xlat,xlong,z_at_w,gd_cloud,gd_cloud2,                      &
               ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,   &
               ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,       &
               ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,            &
               ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,   &
               pm2_5_dry,pm2_5_water,uvrad,                               &
               ids,ide, jds,jde, kds,kde,                                 &
               ims,ime, jms,jme, kms,kme,                                 &
               its,ite, jts,jte, kts,kte                                  )
   USE module_configure
   USE module_state_description
   USE module_model_constants
   USE module_data_radm2
   implicit none
   INTEGER,      INTENT(IN   ) :: id,julday,                    &
                                  ids,ide, jds,jde, kds,kde,    &
                                  ims,ime, jms,jme, kms,kme,    &
                                  its,ite, jts,jte, kts,kte
   INTEGER,      INTENT(IN   ) :: ktau
   REAL(KIND=8), INTENT(IN   ) :: curr_secs
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),        &
         INTENT(IN ) ::                                   moist
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                   &
         INTENT(INOUT ) ::                                         &
               pm2_5_dry,pm2_5_water
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                   &
         OPTIONAL,                                                 &
         INTENT(INOUT ) ::                                         &
               gd_cloud,gd_cloud2
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                   &
         INTENT(INOUT ) ::                                         &
           ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,&
           ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,    &
           ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,         &
           ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),         &
         INTENT(IN ) ::                                   chem
   REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,    &
          INTENT(IN   ) ::                                      &
                                                      t_phy,    &
                                                      p_phy,    &
                                                      dz8w,     &
                                              t8w,p8w,z_at_w ,  &     
                                                      aerwrf ,  &     
                                                    rho_phy           
   REAL,  DIMENSION( ims:ime , jms:jme )                   ,    &     
          INTENT(IN   ) ::                                      &     
                                                     xlat,      &
                                                     xlong
   REAL,  DIMENSION( ims:ime , jms:jme )                   ,    &     
          INTENT(INOUT   ) ::                               uvrad
      REAL,      INTENT(IN   ) ::                      gmt

   TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags

   LOGICAL, INTENT(IN) :: haveaer
         
                                                    
! LOCAL  VAR

   INTEGER(KIND=8) :: ixhour
   INTEGER :: ki,i,j,k,n,iprt
!   photolysis input

      real tt(kts:kte+1),o33(kts:kte+1),rhoa(kts:kte+1),aerext(kts:kte+1),qll(kts:kte+1),    &
           phizz(kts:kte+1),phot1(nreakj-1,kts:kte+1)
      real(kind=8) :: xtime,xhour
      real :: xmin,gmtp,uvb_dd1,uvb_du1,uvb_dir1
      real :: zenith,zenita,azimuth,dobsi
      real :: bext340,bexth2o,ctr
      integer :: naerspec
      real :: zsurf

!     print *,'gmt,julday in madronich1= ',gmt,julday
      xtime=curr_secs/60._8   !min since simulation start
      ixhour=int(gmt+.01,8)+int(xtime/60._8,8)
      xhour=real(ixhour,8)
      xmin=60.*gmt+real(xtime-xhour*60._8,8)
      gmtp=mod(xhour,24._8)
      gmtp=gmtp+xmin/60.
!     print *,'gmtp = ',gmtp,xhour,xmin
      do 100 j=jts,jte
      do 100 i=its,ite
!     write(0,*)i,j
        do k=kts,kte+1
        do n=1,nreakj-1
          phot1(n,k)=0.
        END DO
        END DO
        iprt = 0
        zenith=0.
        zenita=0.
        azimuth=0.
        call calc_zenith(xlat(i,j),-xlong(i,j),julday,gmtp,azimuth,zenith)
! if nighttime, skip radiative transfer calculation
        if(zenith.eq.90.) zenith = 89.9
        if(zenith.ge.90.) go to 199
        zenita = cos(zenith*pi/180.)
        if(zenith.gt.75.)   zenita=1./chap(zenith)
        zsurf = z_at_w(i,1,j)*0.001

! photmad berechnet photolysefrequenzen nach Madronich in folgender Reihenfolge
        
! o2 absorp       1             schumman-runge corrected in srband
! o3 -> 1d        2             at 275 k. correct t-dep in subgrid
! o3 -> 3p        3             at 275 k. correct t-dep in subgrid
! no2             4
! no3 -> no+o2    5
! no3 -> no2+o    6
! hno2            7
! hno3            8
! hno4            9
! h2o2           10
! ch2o -> rad    11
! ch2o -> mol    12
! ch3cho         13
! ch3coch3       14
! ch3coc2h5      15
! hcocho proc a  16
! ch3cocho       17
! hcoch=chcho    18              estimate, no reliable measurement
! ch3o2h         19
! ch3coo2h       20              actually use 0.28*(h2o2 value)
! ch3ono2        21
! hcocho proc b  22
        do k=kts,kte
          aerext(k+1)=aerwrf(i,k+1,j)
!
!--- if you have aerosols
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          bext340=5.E-6
          bexth2o=5.E-6
          if(haveaer.and.ktau.gt.1)then
 
! dry aerosol mass
!rf check ki (or ki+1 or ki-1 ?)
                aerext(k)=pm2_5_dry(i,k,j)*bext340+ &
                   pm2_5_water(i,k,j)*bexth2o
                aerext(k)=aerext(k)*1.E3
!                if(i.eq.70.and.j.eq.70) write(06,*) 'aerext',k,aerext(k)

          endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          qll(k+1)=0.
          tt(k+1) = t_phy(i,k,j)
          rhoa(k+1) = rho_phy(i,k,j)
          o33(k+1) = max(1.e-3,chem(i,k,j,p_o3))
         if(present(gd_cloud) .and. present(gd_cloud2))then
          qll(k+1) = 1.e3*(moist(i,k,j,p_qc)+moist(i,k,j,p_qi)+      &
                           gd_cloud(i,k,j)+gd_cloud2(i,k,j))         &
                     *rho_phy(i,k,j)
         else
          qll(k+1) = 1.e3*(moist(i,k,j,p_qc)+moist(i,k,j,p_qi))      &
                     *rho_phy(i,k,j)
         endif
          if(qll(k+1).lt.1.e-5)qll(k+1) = 0.
          phizz(k+1) = z_at_w(i,k+1,j)*.001-z_at_w(i,1,j)*.001
!         if((i.eq.1.and.j.eq.17))then
!            write(0,*)k+1,phizz(k+1),qll(k+1),moist(i,k,j,p_qc),moist(i,k,j,p_qi),rqccuten(i,k,j),rqicuten(i,k,j),rho_phy(i,k,j)
!         write(0,*)k+1,z_at_w(i,k+1,j),z_at_w(i,1,j),tt(k+1),o33(k+1),qll(k+1),rhoa(k+1)
!         endif
       END DO
        tt(1)=t8w(i,kts,j)
        o33(1)=max(1.e-3,chem(i,kts,j,p_o3))
        qll(1)=0.
        phizz(1)=0.
        aerext(1)=aerwrf(i,1,j)
        k=0
!         write(0,*)k+1,z_at_w(i,k+1,j),z_at_w(i,1,j),tt(k+1),o33(k+1),qll(k+1),rhoa(k+1)
!
! if you have aerosols....
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
        if(haveaer.and.ktau.gt.1)then
              aerext(1)=pm2_5_dry(i,1,j)*bext340+ &
              pm2_5_water(i,1,j)*bexth2o
              aerext(1)=aerext(1)*1.e3
        endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        rhoa(1)=p8w(i,kts,j)/t8w(i,kts,j)/r_d
        dobsi=350.
!         if((i.eq.87.and.j.eq.66).or.(i.eq.105.and.j.eq.70))then
!            print *,'before photmad, i,j = ',i,j,pm2_5_dry(i,1,j),pm2_5_water(i,1,j)
!            print *,k,rhoa(1),phizz(1),qll(1),aerext(1),o33(1),tt(1),zenita
!         endif
!write(0,*)'calling photolysis_mad --------',i,j
        call photolysis_mad(kte,zenita,zenith,zsurf, &
             phizz,tt,rhoa,o33,aerext,qll,dobsi,phot1, &
             nreacj,iprt,uvb_dd1,uvb_du1,uvb_dir1)
!write(0,*)'back from photolysis_mad ---------------------------- '
!       print *,'after photmad, i,j = ',i,j
        uvrad(i,j)=uvb_dd1+uvb_dir1-uvb_du1
        do k=kts,kte
        do n=1,nreakj-1
          phot1(n,k)=60.*phot1(n,k)
        END DO
        END DO
 199    continue
        do k=kts,kte
!
!
          ph_o31d(i,k,j) = phot1(1,k)
          ph_o33p(i,k,j) = phot1(2,k)
          ph_no2(i,k,j) = phot1(3,k)
          ph_no3o2(i,k,j) = phot1(4,k)
          ph_no3o(i,k,j) = phot1(5,k)
          ph_hno2(i,k,j) = phot1(6,k)
          ph_hno3(i,k,j) = phot1(7,k)
          ph_hno4(i,k,j) = phot1(8,k)
          ph_h2o2(i,k,j) = phot1(9,k)
          ph_ch2or(i,k,j) = phot1(10,k)
          ph_ch2om(i,k,j) = phot1(11,k)
          ph_ch3cho(i,k,j) = phot1(12,k)                              
          ph_ch3coch3(i,k,j) = phot1(13,k)
          ph_ch3coc2h5(i,k,j) = phot1(14,k)
          ph_hcocho(i,k,j) = phot1(15,k)
          ph_ch3cocho(i,k,j) = phot1(16,k)
          ph_hcochest(i,k,j) = phot1(17,k)
          ph_ch3o2h(i,k,j) = phot1(18,k)
          ph_ch3coo2h(i,k,j) = phot1(19,k)
          ph_ch3ono2(i,k,j) = phot1(20,k)
          ph_hcochob(i,k,j) = phot1(21,k)
          ph_macr(i,k,j) = phot1(22,k)
!         if(i.eq.5.and.j.eq.5)print *,i,j,k,phot1(3,k),phot1(4,k),phot1(17,k)
        END DO
 100    continue

END SUBROUTINE madronich1_driver

      SUBROUTINE photolysis_mad(mkxcc,a,zenith,zsurf,         &
          zmm5,tmm5,pmm5,o3mm5,aerext,wlmm5,                  &
          dobsnew,phot1,nrtest,iprt,uvb_dd1,uvb_du1,uvb_dir1)
        ! Note that nreakj can differ from nreacj in csolve1!!!!!!
        ! ----- INPUT ------------------------------------------------
        ! Input from MM5
        ! ------ OUTPUT ----------------------------------------------
        ! -----------------------------------------------------------------------
        ! This program calculates J-values for selected atmospheric molecules
        ! the program and subroutine structure is as follows
        ! runph  -main program and menu.
        ! calc_zenith  -zenith angle calc.
        ! in addition, the chapman function ch(zenith) is used between
        ! readd  -reads spectra and standard profiles
        ! photmat -main subroutine (calls the routines below)
        ! o3scal  -rescales ozone profile to user-selected new dobson
        ! subgrid    -regrids the altitude profiles (
        ! srband -computes effective ozone cross sections in the
        ! schumamkxcc+1+nabv-runge region (if necessary)
        ! trapez -interpolation subroutine
        ! optics    -computes optical parameters and calls delta-eddington
        ! delted  -delta-eddington code of wiscombe
        ! nlayde(ib)  -multylayer calc
        ! leqt1b(a,n,ncl,nuc,ia,b,m,ib,ijob,kl)  -solves
        ! pentadiag.system
        ! Original program written by S. Madronich, NCAR. The code was
        ! last modified by him on 18 Aug. 1987.
        ! The driver for the chemistry has been heavily modified by
        ! W. R. Stockwell at IFU Germany.  (The modifications allow
        ! the simulation conditions to be directly modified
        ! without recompiling the program [this does not hold for the present
        ! version here any more]).  The extremely large data base
        ! has bee replaced with separate files to allow cross sections and
        ! photolysis rates to be more easily updated.  However, the order
        ! of the input files should be modified with extreme caution.
        ! The present version is made for the use of computed
        ! fields of T, P, O3, wl, and aerosol (if available)
        ! Modifications were made by Renate Forkel in Aug./Sept 1995
        ! The following modifications were made:
        ! - No climatological input any more
        ! - Arbitrary vertical grid
        ! - Several seperate cloud layers are possible. Number of additional
        ! layers within the clouds depends on the optical depth
        ! - Vertically inhogeneous clouds are possible now
        ! - No variable values in commons any more
        ! -----------------------------------------------------------------------
        ! xnkg: Molecules per kg air
        ! .. Scalar Arguments ..
        REAL :: a, dobsnew,uvb_dd1,uvb_du1,uvb_dir1
        INTEGER :: iprt, mkxcc, nrtest
        ! .. Local Scalars ..
        REAL :: aeru, airu, bextn, df, dj, dz, ff, gaer, gcld, gray, haer, &
          hair, ho3, o3u, omaer, omcld, omray, reff, tu, wmicron, xnkg, xx, &
          znorm, zu
        INTEGER :: i, j, k, kk, kl, lev, nlayer, nlevel, nn, nr, nsurf
        ! .. Array Arguments ..
        REAL :: aerext(mkxcc+1), o3mm5(mkxcc+1), phot1(nreakj-1,mkxcc+1), &
          pmm5(mkxcc+1), tmm5(mkxcc+1), wlmm5(mkxcc+1), zmm5(mkxcc+1)
        ! .. Local Arrays ..
        REAL :: aaer(130), aer(mkxcc+1+nabv), air(mkxcc+1+nabv), ao2(nj,130), &
          ao3(nj,130), arayl(130), cloud(mkxcc+1+nabv), cvo2(nj), &
          d(nreakj,nj), endir(nj,130), endn(nj,130), enup(nj,130), &
          hilf1(mkxcc+1), hilfd(nj), o3(mkxcc+1+nabv), qy(nj,130,nreakj), &
          s(nj,130,nreakj), t(mkxcc+1+nabv), vaer(nj), vair(nj), vcld(nj), &
          vo3(nj), vt(nj), z(nj), zkm(mkxcc+1), zmid(nj), zz(mkxcc+1+nabv)
 
          REAL :: fldir(nj,130), fldn(nj,130), flup(nj,130)

        ! .. Data Statements ..
        DATA xnkg/2.143E25/

        REAL :: zenith, zsurf
        REAL zasl(nj)

        ! EXTERNAL sphers, airmas
        INTEGER, PARAMETER :: kzmax = 200
        REAL :: dsdh(0:kzmax,kzmax)
        INTEGER :: nid(0:kzmax)
        INTEGER :: ii
        REAL vcol(nj), scol(nj)

        ! EXTERNAL o3scal, optics, srband, subgrid, trapez
        ! wave length range !!! If photolysis is also desired for levels above
        ! 2
        ! kl0 should be set equal to 1 again!!!!!!!!!!!!!!
        i = 1
        j = 1
        ! albedoph ************************  specify ground albedoph
        ! use best estimate albedoph of demerjian et al.,
        ! adv.env.sci.tech.,v.10,p.369, (1980)
        DO kl = kl0, kl1

          IF (wl(kl)<400.) albedoph(kl) = 0.05

          IF ((wl(kl)>=400.) .AND. (wl(kl)<450.)) albedoph(kl) = 0.06

          IF ((wl(kl)>=450.) .AND. (wl(kl)<500.)) albedoph(kl) = 0.08

          IF ((wl(kl)>=500.) .AND. (wl(kl)<550.)) albedoph(kl) = 0.10

          IF ((wl(kl)>=550.) .AND. (wl(kl)<600.)) albedoph(kl) = 0.11

          IF ((wl(kl)>=600.) .AND. (wl(kl)<640.)) albedoph(kl) = 0.12

          IF ((wl(kl)>=640.) .AND. (wl(kl)<660.)) albedoph(kl) = 0.135

          IF (wl(kl)>=660.) albedoph(kl) = 0.15
        END DO



        nn = mkxcc + 1 + nabv

        ! omray = single scattering albedoph, rayleigh.  use 1.00
        ! gray =  asymetry factor for rayleigh scattering.  use 0.0
        ! arayl(kl) = rayleigh scattering cross section, from
        ! frohlich and shaw, appl.opt. v.11, p.1773 (1980).
        ! overrides tabulation of jdata.base
        hair = 8.05
        omray = 1.0
        gray = 0.0

        DO 10 kl = kl0, kl1
          wmicron = wl(kl)/1.E3
          xx = 3.916 + 0.074*wmicron + 0.050/wmicron
          arayl(kl) = 3.90E-28/(wmicron)**xx
10      CONTINUE

        ! aerosol***********************  specify aerosols
        ! aaer(kl) =  aerosol total vertical optical depth variation with
        ! wavelength.  estimated from elterman (1968)
        ! aer(i) = attenuation (per km) profile from elterman (1968).
        ! given in data statement in beginning of code, for 340 nm (kl=6
        ! same vertical shape at all wavelengths.
        ! normalized later (in subroutine subgrid) to total vertical dep
        ! this wavelength.
        ! omaer = aerosol single scattering albedoph.  use 0.99 for now.
        ! gaer = aerosol asymetry factor.  use 0.61 (hansen and travis 197
        ! (these are assuming particles of about 0.1 micron radius
        ! index of refraction of about 1.65 + 0.002i.
        ! haer = the aerosol scale height at top of atmosphere
        ! use equal to air (8.05 km)

        DO 20 kl = kl0, kl1
          aaer(kl) = 0.379*(340./wl(kl))
20      CONTINUE
        omaer = 0.990
        gaer = 0.610
        haer = 8.05


        omcld = 1.000
        gcld = 0.860

        ho3 = 4.50

        nsurf = 1

        ! Vertikales Gitter ohne Wolken: Niveaus zz(i), i=1,mkxcc+1+nabv


        ! Transformation of MM5-values
        ! bextn: cloud extinction coeff per g/m**3 [1/km]
        DO k = 1, mkxcc + 1
          zz(k) = zmm5(k)
!write(0,*)' here7a zz = zmm5 ',k,zmm5(k)
          t(k) = tmm5(k)
          air(k) = xnkg*pmm5(k)*1.E-6

          ! falls o3mm5 in ppm

          o3(k) = o3mm5(k)*1.E-6*air(k)
          aer(k) = aerext(k)

          ! bextn: cloud extinction coeff per g/m**3 [1/km]
          ! **** warning: parameterization for bextn only good
          ! for continental clouds

          IF (wlmm5(k)>0.) THEN

            ! vereinfachte Version, falls kein
            ! Sulfat uebergeben wird.

            reff = 9.6*wlmm5(k)**0.333
            bextn = (0.0275+1.3/reff)*1000.
            cloud(k) = wlmm5(k)*bextn
          ELSE
            cloud(k) = 0.
          END IF
          ! if(iprt.eq.1)write(6,'(i3,e12.4)')k,cloud(k)
        END DO

        znorm = (50.-zz(mkxcc+1))/(50.-20.)

        DO k = 1, nabv
          zz(mkxcc+1+k) = 50. - znorm*(50.-zabv(k))
!write(0,*)' here8a ',k,' zz(',mkxcc+1+k,') ',zz(mkxcc+1+k),znorm,zabv(k)
        END DO

        zu = zz(mkxcc+1)
        tu = t(mkxcc+1)
        o3u = o3(mkxcc+1)
        airu = air(mkxcc+1)
        aeru = aer(mkxcc+1)
        kk = 1

        ! Zufuegen von Werten oberhalb von MM5

        ! Die 'abv'-Werte sind bereits in den richtigen Einheiten
        DO k = mkxcc + 1 + 1, mkxcc + 1 + nabv
!           write(0,'(2i3,5e12.3)')  k,kk,zz(k),zabv(kk)
30        IF (zz(k)<=zabv(kk)) THEN
            dz = zz(k) - zu
            ff = dz/(zabv(kk)-zu)
            t(k) = tu + ff*(tabv(kk)-tu)
            o3(k) = o3u + ff*(o3abv(kk)-o3u)
            air(k) = airu + ff*(pabv(kk)-airu)
            aer(k) = aeru + ff*(caabv(kk)-aeru)
            cloud(k) = 0.
            ! if(iprt.eq.1)then
!            write(0,'(2i3,5e12.3)')  k,kk,zz(k),zabv(kk),ff,tabv(kk),air(k)
            ! endif
          ELSE
40          zu = zabv(kk)
            tu = tabv(kk)
            o3u = o3abv(kk)
            airu = pabv(kk)
            aeru = caabv(kk)
            kk = kk + 1

            IF (zabv(kk)<zz(k)) GO TO 40
            GO TO 30
          END IF

        END DO

        IF (iprt==1) THEN
          !WRITE (06,'('k, zz(k), t(k),  air(k), o3(k), aer(k) cloud(k)')')

          !DO k = 1, mkxcc + 1 + nabv
          !  WRITE (06,'(i3,6e12.4)') k, zz(k), t(k), air(k), o3(k), aer(k), &
          !    cloud(k)
          !END DO

        END IF
        ! if(iprt.eq.1)then
        ! do k=1,mkxcc+1
        ! write(06,'(i3,6e12.4)') k,zz(k),t(k),air(k),o3mm5(k), &
        ! aerext(k),cloud(k)
        ! enddo
        ! endif
        ! Scaling of O3-column with specified Dobson units
        ! -----------------------------------------------------------------------
        CALL o3scal(dobsnew,ho3,zz,o3,mkxcc+1+nabv)
        ! stop
        ! -----------------------------------------------------------------------

        ! Einfuegen zusaetzlicher Schichten in Wolken,
        ! T-abh. Absorptionsquerschnitte, etc

        ! -----------------------------------------------------------------------
        nlevel = 0

        IF (iprt==1) PRINT *, 'nlevel = ', nlevel

        CALL subgrid(zz,t,air,aer,o3,cloud,hair,haer,z,zmid,vair,vo3,vaer, &
          vcld,vt,ao3,nlevel,nlayer,s,qy,cvo2,omray,mkxcc+1,mkxcc+1+nabv)

        IF (iprt==1) PRINT *, 'nlevel = ', nlevel
        ! -----------------------------------------------------------------------

! compute altitudes relative to sea level
        DO ii = 1, nlevel
           zasl(ii) = z(ii) + zsurf
        ENDDO
! compute spherical increments
        CALL sphers(kzmax, nlevel, z, zenith, dsdh, nid)
! here can also call airmass to compute air mass factors for Lyman alpha and 
! Schmumann-Runge parameterizations
!     CALL airmas(nlevel, dsdh, nid, vair, vcol, scol)

        ! subroutine srband computes effective ozone cross sections in the
        ! schumann-runge region, using the parameterization of allen and
        ! freder

        ! ----------------------------------------------------------------------
        CALL srband(a,zmid,cvo2,ao2,vt,nlevel)
        ! ----------------------------------------------------------------------


        ! subroutine optics is the driver for the flux calculation.  output is
        ! endir(lev,kl)  -irradiance of direct solar beam
        ! endn(lev,kl)   -irradiance of down-welling diffuse light
        ! enup(lev,kl)   -irradiance of up-welling diffuse light
        ! for each level lev and wavelength bin kl. conversion to actinic
        ! flux is in loop 90

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

        CALL optics(iprt, &
          zenith, a, dsdh, nid, &
          vair,arayl,gray,omray,ao2,vo3,ao3,vcld,gcld,omcld, &
          vaer,aaer,gaer,omaer,nlevel, kzmax, &
          endir,endn,enup, fldir,fldn,flup)

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

        DO nr = 1, nreakj

          DO k = 1, nlevel
            d(nr,k) = 0.
          END DO

        END DO

        DO 50 kl = kl0, kl1

          DO 50 lev = 1, nlevel

! calculate actinic flux, direct + diffuse-dn + diffuse-up
            df = fext(kl)*(fldir(lev,kl) + fldn(lev,kl) + flup(lev,kl))

! calculate contribution to j-vals, sum over wavelengths

            DO 50 nr = 1, nreakj
              dj = df*s(lev,kl,nr)*qy(lev,kl,nr)
              d(nr,lev) = d(nr,lev) + dj
50      CONTINUE

! calculate spectral irradiance at surface:
         uvb_dd1=0.
         uvb_du1=0.
         uvb_dir1=0.   
        do  kl = 38,54 
!        radfakt=fext(kl) ! photons/cm**2/s/bin
!                         ! UV in Photonen/cm**2/s
         radfakt=fext(kl)*1e4 * 6.63e-34 * 2.9979e8 / (wl(kl)*1e-9)
                          ! UV in W/m**2
         uvb_dd1=uvb_dd1+radfakt*endn(1,kl)
         uvb_du1=uvb_du1+radfakt*enup(1,kl)
         uvb_dir1=uvb_dir1+radfakt*endir(1,kl)
        enddo

        ! print *,'d(4,1) = ',d(4,1)

! interpolate j-values to new altitudes:

        DO nr = 1, nreakj - 1
          ! write(06,'(' nr',i3)') nr
          DO lev = 1, mkxcc + 1
            zkm(lev) = zmm5(lev)
          END DO

          DO lev = 1, nlevel
            hilfd(lev) = d(nr+1,lev)
            ! if(iprt.eq.1.and.nr.eq.3)
            ! 1      write(06,'(i3,2e12.4)') lev,z(lev),d(4,lev)
          END DO
          ! hilfd muss noch auf mm5-hoehen zurueckinterpoliert werden
          CALL trapez(z,hilfd,1,nlevel,zkm,hilf1,1,mkxcc+1,nj,nj,mkxcc+1, &
            mkxcc+1)

          IF (iprt==1 .AND. nr==3) PRINT *, 'nach trapez '

          IF (iprt==1 .AND. nr==3) WRITE (06,'(5e12.6)') (hilf1(k),k=1, &
            mkxcc+1)

          IF (iprt==1 .AND. nr==3) PRINT *, 'nach trapez '

          DO lev = 1, mkxcc + 1
            phot1(nr,lev) = hilf1(lev)
          END DO

        END DO

        RETURN

      END SUBROUTINE photolysis_mad

      ! #################################################################

      SUBROUTINE srband(a,zmid,cvo2,ao2,vt,nlevel)
        ! this subroutine calculates effective o2 cross sections in the
        ! schumann-runge band, using the formulae from allen and frederick,
        ! j.atmos.sci., v39, p2066 (1982).  the coefficients are those modifi
        ! wuebbles (lll report, 1982).
        ! .. Scalar Arguments ..
        REAL :: a
        INTEGER :: nlevel
        ! .. Local Scalars ..
        REAL :: ao20, ao20lg, c, clog, e10, x1, x2, x3, x4, x5, x6, x7, x8, &
          y1, y2, y3, y4, zendep
        INTEGER :: lay, n20, nlayer
        ! .. Intrinsic Functions ..
        INTRINSIC alog
        ! .. Array Arguments ..
        REAL :: ao2(nj,130), cvo2(nj), vt(nj), zmid(nj)

        nlayer = nlevel - 1
        ! initialize cross sections:
        DO 10 lay = 1, nlayer

          DO 10 kl = kl0, kl1
10      ao2(lay,kl) = xs(kl,1)


        ! correct as needed
        ! use formula for 20-50 km.
        ! below 20 km, use 20 km value
        ! find layer near 20 km
        DO 20 lay = 1, nlayer

          IF (zmid(lay)>20.) THEN
            n20 = lay
            GO TO 30
          END IF

20      CONTINUE
30      CONTINUE

        e10 = alog(10.)

        DO 60 kl = kl0, kl1

          IF (wl(kl)>205.) RETURN

          DO 40 lay = n20, nlayer
            x1 = alog(4.696E-23*cvo2(lay)/0.2095)/e10

            IF (wl(kl)>=200.) x1 = vt(lay)
            x2 = x1*x1
            x3 = x2*x1
            x4 = x3*x1
            x5 = x4*x1
            x6 = x5*x1
            x7 = x6*x1
            x8 = x7*x1
            ao20lg = sra(kl,1) + sra(kl,2)*x1 + sra(kl,3)*x2 + sra(kl,4)*x3 + &
              sra(kl,5)*x4 + sra(kl,6)*x5 + sra(kl,7)*x6 + sra(kl,8)*x7 + &
              sra(kl,9)*x8
            ao20 = 10.**ao20lg

            y1 = alog(cvo2(lay))/e10
            y2 = y1*y1
            y3 = y2*y1
            y4 = y3*y1
            clog = srb(kl,1) + srb(kl,2)*y1 + srb(kl,3)*y2 + srb(kl,4)*y3 + &
              srb(kl,5)*y4
            c = 10.**clog
            zendep = a**c

            ao2(lay,kl) = ao20*zendep
40        CONTINUE
          ! assign values below 20 km
          DO 50 lay = 1, n20 - 1
            ao2(lay,kl) = ao2(n20,kl)
50        CONTINUE
60      CONTINUE

        RETURN

      END SUBROUTINE srband

      ! ######################################################################

      SUBROUTINE o3scal(dobsnew,ho3,zz,o3,nn)
        ! adjustment of o3 profiles to a user-selected dobson value.
        ! select value of dobnew in main program
        ! if don't want to use, don't call this subroutine
        ! .. Scalar Arguments ..
        REAL :: dobsnew, ho3
        INTEGER :: nn
        ! .. Local Scalars ..
        REAL :: dobsref
        INTEGER :: i
        ! .. Intrinsic Functions ..
        INTRINSIC max, min
        ! .. Array Arguments ..
        REAL :: o3(nn), zz(nn)
        ! write(6,*) o3
        dobsref = o3(nn)*1.E5*ho3
        ! write(06,'('nn: dobsref,dobsnew',2e12.4)')
        ! &   dobsref/2.687e16,dobsnew
        DO 10 i = 1, nn
10      dobsref = dobsref + o3(i)*0.5*(zz(min(i+1,nn))-zz(max(i-1,1)))*1.E5
        dobsref = dobsref/2.687E16
        ! write(06,'('dobsref,dobsnew',2e12.4)') dobsref,dobsnew
        DO 20 i = 1, nn
          o3(i) = o3(i)*dobsnew/dobsref
20      CONTINUE
        ! write(06,*) o3

        RETURN

      END SUBROUTINE o3scal

      ! #######################################################################

      FUNCTION chap(zeta)
        ! chapman function is used when the solar zenith angle exceeds
        ! 75 deg.
        ! interpolates between values given in, e.g., mccartney (1976).
        ! .. Scalar Arguments ..
        REAL :: zeta
        ! .. Local Scalars ..
        REAL :: rm
        INTEGER :: i
        ! .. Local Arrays ..
        REAL :: y(22)
        ! .. Function Return Value ..
        REAL :: chap
        ! .. Data Statements ..
        DATA (y(i),i=1,22)/3.800, 4.055, 4.348, 4.687, 5.083, 5.551, 6.113, &
          6.799, 7.650, 8.732, 10.144, 12.051, 14.730, 18.686, 24.905, 35.466, &
          55.211, 96.753, 197., 485., 1476., 9999./

        DO 10 i = 75, 96
          rm = i

          IF (zeta<rm) GO TO 20
10      CONTINUE
20      chap = y(i-75) + (y(i-74)-y(i-75))*(zeta-(rm-1.))
      END FUNCTION chap


      ! #################################################################
      SUBROUTINE subgrid(zz,t,air,aer,o3,cloud,hair,haer,z,zmid,vair,vo3,vaer, &
          vcld,vt,ao3,nlevel,nlayer,s,qy,cvo2,omray,nz,nn)
        ! Output
        ! this subroutine computes all altitude dependent quantities over the
        ! sub-divided grid
        ! the following quantities are computed at each level (altitude)
        ! z(lev) = altitude (km) of each level
        ! zair(lev) = air concentration at each level
        ! zt(lev)   = temperature at each level
        ! s(lev,kl,nr) = abs. cross sect at each altitude (t,p corrected)
        ! qy(lev,kl,nr) = quantum yield at each altitude (t,p corrected)
        ! the following quantities are computed at each layer (thickness)
        ! zmid(lay) = altitude of midpoint of layer
        ! vair(lay) = air column in layer (vertical)
        ! vo3(lay)  = ozone column   '       '
        ! vaer(lay) = aerosol '      '       '
        ! vcld(lay) = cloud   '      '       '
        ! vt(lay)   = average temperature of column
        ! ao3(lay,kl) = average o3 cross sect in layer (with ave. layer te
        ! .. Scalar Arguments ..
        REAL :: haer, hair, omray
        INTEGER :: nlayer, nlevel, nn, nz
        ! .. Local Scalars ..
        REAL :: a, ak300, akt, b, c, dz, dzo, dzt, dzu, fnum, hlocal, phi1, &
          phi2, phi20, sum, tau, tdiffx, x0, xl, xl0
        INTEGER :: i, idt, ii, il, lay, lev, llbt, nr
        ! .. Intrinsic Functions ..
        ! EXTERNAL trapez
        INTRINSIC alog, atan, exp, float, ifix, max
        ! .. Array Arguments ..
        REAL :: aer(nn), air(nn), ao3(nj,130), cloud(nn), cvo2(nj), o3(nn), &
          qy(nj,130,nreakj), s(nj,130,nreakj), t(nn), vaer(nj), vair(nj), &
          vcld(nj), vo3(nj), vt(nj), z(nj), zmid(nj), zz(nn)
        ! .. Local Arrays ..
        REAL :: zair(nj), zt(nj)
        ! ------------------------------------------------------  levels
        ! fnum: factor for calculation of additional levels in clouds
        fnum = 0.02

        DO i = 1, nj
          vcld(i) = 0.
        END DO

        lev = 0
        llbt = 0

        DO 30 i = 1, nn
          ! write(06,'('cloud',i3,2e12.4)') i,zz(i),cloud(i)
          IF (cloud(i)<=0.) THEN
            lev = lev + 1
            z(lev) = zz(i)
            zair(lev) = air(i)
            zt(lev) = t(i)
            ! write(06,'('z,t,air ',i3,3e12.4)')
            ! lev,z(lev),zt(lev),zair(lev)
          ELSE

            IF (i==1 .AND. lev==0) THEN
              llbt = llbt + 1
              ! lbase(llbt)=1
              lev = 1
              z(lev) = zz(i)
              zair(lev) = air(i)
              zt(lev) = t(i)
            END IF

            IF (i==1) GO TO 10

            IF (cloud(i-1)<=0.) THEN
              lev = lev + 1
              llbt = llbt + 1
              ! lbase(llbt)=lev
              z(lev) = 0.5*(zz(i)+zz(i-1))
              zt(lev) = 0.5*(t(i)+t(i-1))
              if(abs(air(i)-air(i-1)).lt.air(i-1)/1.e5)then
                 zair(lev) = air(i-1)
              else
                 hlocal = 1./alog(air(i-1)/air(i))
                 x0 = (z(lev)-zz(i-1))/(zz(i)+zz(i-1))
                 zair(lev) = air(i-1)*exp(-x0/hlocal)
              endif

              ! write(06,'('u:z,t,air ',i3,3e12.4)')
              ! lev,z(lev),zt(lev),zair(lev)
            END IF

            dzt = zz(i) - z(lev)
            ! number of layers depents on optical depth
            idt = max(ifix(cloud(i)*dzt*fnum),1)
            dzu = dzt/float(idt)
            vcld(lev) = cloud(i)

            DO il = 1, idt
              lev = lev + 1

              IF (lev>nj) THEN
                CALL wrf_error_fatal ( 'LEV > NJ, NJ GROESSER WAEHLEN')
              END IF

              z(lev) = z(lev-1) + dzu
              zt(lev) = t(i-1) + (z(lev)-zz(i-1))/(zz(i)-zz(i-1))*(t(i)-t(i-1) &
                )

              if(abs(air(i)-air(i-1)).lt.air(i-1)/1.e5)then
                 zair(lev) = air(i-1)
              else
                 hlocal = 1./alog(air(i-1)/air(i))
                 x0 = (z(lev)-zz(i-1))/(zz(i)+zz(i-1))
                 zair(lev) = air(i-1)*exp(-x0/hlocal)
              endif
              vcld(lev) = cloud(i)
              ! write(06,'('u:z,t,air ',i3,3e12.4)')
              ! lev,z(lev),zt(lev),zair(lev)
            END DO

10          CONTINUE

            IF (i==nn) GO TO 20
            ! if(lev.ne.1) vcld(lev)=cloud(i)
            dzt = (zz(i+1)-zz(i))*.5
            ! number of layers depents on optical depth
            idt = max(ifix(cloud(i)*dzt*fnum),1)
            dzo = dzt/float(idt)

            DO il = 1, idt
              lev = lev + 1

              IF (lev>nj) THEN
                CALL wrf_error_fatal ( 'LEV > NJ, NJ GROESSER WAEHLEN')
              END IF

              z(lev) = z(lev-1) + dzo
              zt(lev) = t(i) + (z(lev)-zz(i))/(zz(i+1)-zz(i))*(t(i+1)-t(i))
              if(abs(air(i)-air(i+1)).lt.air(i)/1.e5)then
                 zair(lev) = air(i)
              else
                 hlocal = 1./alog(air(i)/air(i+1))
                 x0 = (z(lev)-zz(i))/(zz(i+1)+zz(i))
                 zair(lev) = air(i)*exp(-x0/hlocal)
              endif
              vcld(lev-1) = cloud(i)
              ! write(06,'('o:z,t,air ',i3,3e12.4)')
              ! lev,z(lev),zt(lev),zair(lev)
            END DO

20          CONTINUE
          END IF
          ! write(06,'('o:z,t,air ',i3,3e12.4)') lev,z(lev),zt(lev),zair(lev)
30      CONTINUE

        ! number of levels including additional cloud levels

        nlevel = lev

        IF (nlevel>nj) print *, ' NLEVEL > NJ, NJ GROESSER WAEHLEN ', &
          nlevel

        ! write(06,'(' nlevel',i3)') nlevel


        ! assign default yields
        DO 40 nr = 1, nreakj

          DO 40 kl = kl0, kl1

            DO 40 lev = 1, nlevel
40      qy(lev,kl,nr) = xqy(kl,nr)
        ! assign default absorption cross sections
        DO 50 kl = kl0, kl1

          DO 50 lev = 1, nlevel

            DO 50 nr = 1, nreakj
50      s(lev,kl,nr) = xs(kl,nr)
        ! --------------------------------------------------------------------
        ! re-calculate altitude dependent quantum yields.  this currently
        ! applies to
        ! 2=o3->o(1d)
        ! 12=ch2o->h2+co
        ! 13=ch3cho->ch3+cho
        ! 14=ch3coch3
        ! 15=ch3coch2ch3
        ! 16=hcocho -> 0.13 hcho + 1.87 co process a
        ! 17=ch3cocho
        ! 22=hcocho -> 0.45 hcho + 1.55 co + 0.80 ho2 process b

        ! o3 and ketones yield is calculated from fit equations,
        ! while for ch3cho and the dicarbonyls yields are calculated
        ! from the ntp yield by linear adjustment to 1/yield.  the ch2o yield
        ! recalculated only for wavelengths longer than 329 nm. the yields for
        ! 1=o3->o(3p) are calculated as (1.- singlet d yield).

        DO 60 lev = 1, nlevel
          ! o3 ozone:
          tau = zt(lev) - 230.
          a = 0.9*(0.369+2.85E-4*tau+1.28E-5*tau*tau+2.57E-8*tau*tau*tau)
          b = -0.575 + 5.59E-3*tau - 1.439E-5*tau*tau - 3.27E-8*tau*tau*tau
          c = 0.9*(0.518+9.87E-4*tau-3.94E-5*tau*tau+3.91E-7*tau*tau*tau)
          xl0 = 308.20 + 4.4871E-2*tau + 6.9380E-5*tau*tau - &
            2.5452E-6*tau*tau*tau

          DO 60 kl = kl0, kl1
            xl = wl(kl)
            qy(lev,kl,2) = a*atan(b*(xl-xl0)) + c

            IF (qy(lev,kl,2)<0.) qy(lev,kl,2) = 0.0

            IF (qy(lev,kl,2)>0.9) qy(lev,kl,2) = 0.9
            qy(lev,kl,3) = 1.0 - qy(lev,kl,2)
            ! ch2o formaldehyde:
            IF ((xl>=330.) .AND. qy(lev,kl,12)>0.) THEN
              phi1 = qy(lev,kl,11)
              phi2 = qy(lev,kl,12)
              phi20 = 1. - phi1
              ak300 = ((1./phi2)-(1./phi20))/2.54E+19
              akt = ak300*(1.+61.69*(1.-zt(lev)/300.)*(xl/329.-1.))
              qy(lev,kl,12) = 1./((1./phi20)+zair(lev)*akt)
            END IF

            IF (qy(lev,kl,12)>1.) qy(lev,kl,12) = 1.0

            IF (qy(lev,kl,12)<0.) qy(lev,kl,12) = 0.0
            ! ch3cho acetaldehyde:
            IF (xqy(kl,13)/=0.) THEN
              qy(lev,kl,13) = 1./(1.+(1./xqy(kl,13)-1.)*zair(lev)/2.465E19)
            END IF
            ! ch3coch3  acetone:
            qy(lev,kl,14) = 0.0766 + 0.09415*exp(-zair(lev)/3.222E18)
            ! ch3coch2ch3  methyl ethyl ketone:
            qy(lev,kl,15) = qy(lev,kl,14)
            ! hcocho  glyoxal process a:
            IF (xqy(kl,16)/=0.) THEN
              qy(lev,kl,16) = 1./(1.+(1./xqy(kl,16)-1.)*zair(lev)/2.465E19)
            END IF
            ! ch3cocho  methylglyoxal:
            IF (xqy(kl,17)/=0.) THEN
              qy(lev,kl,17) = 1./(1.+(1./xqy(kl,17)-1.)*zair(lev)/2.465E19)
            END IF
            ! hcocho  glyoxal process b:
            IF (xqy(kl,22)/=0.) THEN
              qy(lev,kl,22) = 1./(1.+(1./xqy(kl,22)-1.)*zair(lev)/2.465E19)
            END IF

60      CONTINUE
        ! _______________________________________________________________________
        ! correct absorption cross sections for t and p dep.  for now, do
        ! 2=ozone
        DO 90 kl = kl0, kl1

          DO 80 lev = 1, nlevel

            IF (kl<33 .OR. kl>61) GO TO 70
            tdiffx = zt(lev) - 230.
            s(lev,kl,2) = (so3tx(kl,1)+so3tx(kl,2)*tdiffx+so3tx(kl,3)*tdiffx* &
              tdiffx)*1.0E-18
            s(lev,kl,3) = s(lev,kl,2)
70          CONTINUE
80        CONTINUE
90      CONTINUE

        ! ----------------------------------------------  layers
        nlayer = nlevel - 1
        ! write(06,'('  Layers    ',i3)') nlayer
        lay = 0

        DO 100 i = 1, nlayer
          lay = lay + 1
          dz = z(i+1) - z(i)
          zmid(lay) = z(i) + 0.5*dz
          vt(lay) = (zt(lay+1)+zt(lay))/2.
          vair(lay) = dz*1.E5*(zair(i+1)+zair(i))/2.
100     CONTINUE

        ! vo3(lay) = dz*1.e5*(o3(i+1) + o3(i))/2.  ! umr. dz in cm
        ! vcld(lay) = 0. *dz
        ! vaer(lay) = (aer(i+1)-aer(i))/alog(aer(i+1)/aer(i)) *dz ! bei

        CALL trapez(zz,o3,1,nn,zmid,vo3,1,nlayer,nn,nn,nj,nj)


        CALL trapez(zz,aer,1,nn,zmid,vaer,1,nlayer,nn,nn,nj,nj)

        DO 110 i = 1, nlayer
          lay = i
          dz = z(i+1) - z(i)
          ! write(06,'('layer ',i3,6e12.4)') lay,zmid(lay),vt(lay),
          ! & vair(lay)/dz,vo3(lay),vaer(lay),vcld(lay)
          dz = z(i+1) - z(i)

          ! umr. dz in cm

          vo3(i) = dz*1.E5*vo3(i)
          vcld(i) = vcld(i)*dz
110     vaer(i) = vaer(i)*dz


        ! normalize aerosol optical depth to unity sum
        sum = 0.

        DO 120 lay = 1, nlayer
          sum = sum + vaer(lay)
120     CONTINUE

        DO 130 lay = 1, nlayer
          vaer(lay) = vaer(lay)/sum
130     CONTINUE

        ! calculated vertical column of o2 above the midpoint of each layer:
        ! want to use this for computing the average schumann-runge cross
        ! secti
        ! in each layer.
        ! so use half of current layer and half of previous higher layer


        cvo2(nlayer) = 0.2095*vair(nlayer)/2.

        DO 140 ii = 2, nlayer
          lay = nlayer - ii + 1
          cvo2(lay) = cvo2(lay+1) + 0.2095*(vair(lay)+vair(lay+1))/2.
140     CONTINUE

        ! correct attenuation coefficients for pressure and/or temperature
        ! dep.
        ! for now do only ozone absorption.
        DO 160 kl = kl0, kl1

          DO 150 lay = 1, nlayer
            tdiffx = vt(lay) - 230.
            ao3(lay,kl) = xs(kl,2)

            IF (kl>=33 .AND. kl<=61) ao3(lay,kl) = (so3tx(kl,1)+so3tx(kl,2)* &
              tdiffx+so3tx(kl,3)*tdiffx*tdiffx)*1.0E-18

150       CONTINUE
160     CONTINUE

        ! write(06,'(' z        ')')
        ! write(06,'(5e12.4     )') z
        ! write(06,'(' zmid     ')')
        ! write(06,'(5e12.4     )') zmid
        ! write(06,'(' zt       ')')
        ! write(06,'(5e12.4     )') zt
        ! write(06,'(' vt       ')')
        ! write(06,'(5e12.4     )') vt
        ! write(06,'(' vo3      ')')
        ! write(06,'(5e12.4     )') vo3
        ! write(06,'(' vair     ')')
        ! write(06,'(5e12.4     )') vair
        ! write(06,'(' vaer     ')')
        ! write(06,'(5e12.4     )') vaer


        RETURN

      END SUBROUTINE subgrid

      SUBROUTINE optics(iprt, &
          zenith, a, dsdh, nid, &
          vair,arayl,gray,omray,ao2,vo3,ao3,vcld,gcld, &
          omcld,vaer,aaer,gaer,omaer,nlevel, kzmax, &
          endir,endn,enup, fldir,fldn,flup)

        ! sm this subroutine prepares the data needed for the flux
        ! calculation,
        ! sm     calls the scattering subroutine ps2str. it returns values of
        ! sm     flux(lev,kl) for altitude lev-1, wavelength kl.
        ! sm it calculates the optical depths (vertical)
        ! sm     for each layer, from the vertical profiles of o2, o3,
        ! sm     air, cloud, and aerosol, and from the associated 'cross
        ! sections

        ! .. Scalar Arguments ..

        INTEGER :: iprt, nlevel
        REAL :: zenith
        REAL :: gaer, gcld, gray, omaer, omcld, omray

        ! .. Local Scalars ..
        REAL :: dtabs, dtaer, dtair, dtcld, dto2, dto3, dtscat
        INTEGER :: ii, lay, lev

        ! .. Intrinsic Functions ..
        INTRINSIC amax1

        ! .. Array Arguments ..
        
        INTEGER :: kzmax
        REAL :: dsdh(0:kzmax,kzmax)
        INTEGER :: nid(0:kzmax)
        REAL :: aaer(130), ao2(nj,130), ao3(nj,130), arayl(130), &
          vaer(nj), vair(nj), vcld(nj), vo3(nj)
        REAL :: fldir(nj,130), fldn(nj,130), flup(nj,130)
        REAL :: endir(nj,130), endn(nj,130), enup(nj,130)
        LOGICAL, PARAMETER :: delta = .TRUE.

        ! .. Local Arrays ..
        REAL :: alb, dtau(nj), g(nj), om(nj)
        REAL :: fdr(nj), fup(nj), fdn(nj), edr(nj), eup(nj), edn(nj)
        REAL, PARAMETER :: largest = 1.E+36



        ! loop over wavelengths
        DO 30 kl = kl0, kl1

          ! calculate optical depths for all layers (including cloud
          ! sublayers)

          DO 10 lay = 1, nlevel - 1
            ii = nlevel - lay

            dtair = vair(lay)*arayl(kl)
            dto2 = 0.2095*vair(lay)*ao2(lay,kl)
            dto3 = vo3(lay)*ao3(lay,kl)
            dtcld = vcld(lay)
            dtaer = vaer(lay)*aaer(kl)

            dtscat = dtair + dtcld + dtaer
            dtabs = dto2 + dto3

            dtscat = AMAX1(dtscat, 1./largest)
            dtabs = AMAX1(dtabs, 1./largest)

            dtau(ii) = dtabs + dtscat
            om(ii) = dtscat/dtau(ii)
            om(ii) = AMAX1(om(ii), 1./largest)
            
            g(ii) = (gray*dtair+gcld*dtcld+gaer*dtaer)/dtscat

10        CONTINUE

          alb = albedoph(kl)

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

! call pseudo-spherical 2-stream, set to delta-Eddington in subroutine. 

          CALL ps2str(kzmax,nlevel,zenith,a,alb,dtau,om,g, &
                     dsdh, nid, delta, &
                     fdr, fup, fdn, edr, eup, edn)

          ! return to upright grid
          DO 20 ii = 1, nlevel
            lev = nlevel - ii + 1

            endir(lev,kl) = edr(ii)
            endn(lev,kl) = edn(ii)
            enup(lev,kl) = eup(ii)
            fldir(lev,kl) = fdr(ii)
            fldn(lev,kl) = fdn(ii)
            flup(lev,kl) = fup(ii)

20        CONTINUE

30      CONTINUE

        RETURN

      END SUBROUTINE optics

      ! #######################################################################

      SUBROUTINE calc_zenith(lat,long,ijd,gmt,azimuth,zenith)
        ! this subroutine calculates solar zenith and azimuth angles for a
        ! part
        ! time and location.  must specify:
        ! input:
        ! lat - latitude in decimal degrees
        ! long - longitude in decimal degrees
        ! gmt  - greenwich mean time - decimal military eg.
        ! 22.75 = 45 min after ten pm gmt
        ! output
        ! zenith
        ! azimuth
        ! .. Scalar Arguments ..
        REAL :: azimuth, gmt, lat, long, zenith
        INTEGER :: ijd
        ! .. Local Scalars ..
        REAL :: caz, csz, cw, d, decl, dr, ec, epsi, eqt, eyt, feqt, feqt1, &
          feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi, &
          pi, ra, raz, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr, &
          yt, zpt, zr
        INTEGER :: jd
        ! .. Intrinsic Functions ..
        INTRINSIC acos, atan, cos, min, sin, tan
        ! convert to radians
        pi = 3.1415926535590
        dr = pi/180.
        rlt = lat*dr
        rphi = long*dr

        ! print julian days current 'ijd'

        ! ???? + (yr - yref)

        jd = ijd



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

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

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

        IF (lzgmt>0) azimuth = azimuth + (2*(180.-azimuth))
        ! 200     format(' ',f7.2,2(12x,f7.2))
        RETURN

      END SUBROUTINE calc_zenith

      SUBROUTINE trapez(x,y,ianfa,ienda,u,v,ianfn,iendn,ix,iy,iu,iv)
        ! *  implemented 1992 by ansgar ruggaber, university of munich, frg
        ! *  funded by the german minister of research and technology (bmft)
        ! *  under contract no. 521-4007-07eu-738 8
        ! lineare interpolation of referencefield (x(i),y(i))
        ! to (u(i),v(i))
        ! save
        ! .. Scalar Arguments ..
        INTEGER :: ianfa, ianfn, ienda, iendn, iu, iv, ix, iy
        ! .. Array Arguments ..
        REAL :: u(iu), v(iv), x(ix), y(iy)
        ! .. Local Scalars ..
        REAL :: uumord, vumord, xumord, yumord
        INTEGER :: i, ianf, ianfnn, idrehu, idrehx, iendnn, iordu, iordx, j

        idrehx = 0

        IF (x(ianfa)>=x(ienda)) THEN
          ! das x-feld wird ansteigend geordnet, das y-feld entsprechend
          iordx = (ienda-ianfa+1)/2

          DO i = ianfa, iordx
            xumord = x(i)
            x(i) = x(ienda+1-i)
            x(ienda+1-i) = xumord
            yumord = y(i)
            y(i) = y(ienda+1-i)
            y(ienda+1-i) = yumord
          END DO

          idrehx = 1
        END IF

        idrehu = 0

        IF (u(ianfn)>=u(iendn)) THEN
          ! u-field increasing
          iordu = (iendn-ianfn+1)/2

          DO i = ianfn, iordu
            uumord = u(i)
            u(i) = u(iendn+1-i)
            u(iendn+1-i) = uumord
          END DO

          idrehu = 1
        END IF

        ianfnn = ianfn
10      CONTINUE

        IF (u(ianfnn)<x(ianfa)) THEN
          ! no extrapolation at x(ianfa)
          v(ianfnn) = 1.0E-12
          ianfnn = ianfnn + 1
          GO TO 10
        END IF

        iendnn = iendn
20      CONTINUE

        IF ((1.-1.0E-10)*u(iendnn)>x(ienda)) THEN
          ! no extrapolation at x(ienda)
          v(iendnn) = 1.0E-12
          iendnn = iendnn - 1
          GO TO 20
        END IF

        ianf = ianfa

        DO j = ianfnn, iendnn

          DO i = ianf, ienda

            IF (x(i)-u(j)) 30, 50, 40
30        END DO

          GO TO 70
40        v(j) = y(i-1) + (y(i)-y(i-1))/(x(i)-x(i-1))*(u(j)-x(i-1))
          GO TO 60
50        v(j) = y(i)
60        ianf = i
70      END DO

        IF (idrehx/=0) THEN
          ! x- und y-field in starting position
          DO i = ianfa, iordx
            xumord = x(i)
            x(i) = x(ienda+1-i)
            x(ienda+1-i) = xumord
            yumord = y(i)
            y(i) = y(ienda+1-i)
            y(ienda+1-i) = yumord
          END DO

        END IF

        IF (idrehu/=0) THEN

          DO i = ianfn, iordu
            uumord = u(i)
            u(i) = u(iendn+1-i)
            u(iendn+1-i) = uumord
            vumord = v(i)
            v(i) = v(iendn+1-i)
            v(iendn+1-i) = vumord
          END DO

        END IF

      END SUBROUTINE trapez

      SUBROUTINE photmad_init(z_at_w,aerwrf,g,ids,ide,jds,jde,kds,kde,ims,ime, &
          jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
        ! local stuff
        ! .. Scalar Arguments ..
        REAL, INTENT (IN) :: g
        INTEGER, INTENT (IN) :: ide, ids, ime, ims, ite, its, jde, jds, jme, &
          jms, jte, jts, kde, kds, kme, kms, kte, kts
        ! .. Array Arguments ..
        REAL, INTENT (INOUT) :: aerwrf(ims:ime,kms:kme,jms:jme)
        REAL, INTENT (IN) :: z_at_w(ims:ime,kms:kme,jms:jme)
        ! .. Local Scalars ..
        REAL :: z1
        INTEGER :: i, j, k
        ! .. Local Arrays ..
        REAL :: aerext(kts:kte), phizz(kts:kte), z(kts:kte)

        DO j = jts, jte

          IF (j>jde-1) GO TO 20

          DO i = its, ite

            IF (i>ide-1) GO TO 10

            ! z at w points

            z1 = z_at_w(i,kts,j)
            z(kts)=0.

            DO k = kts+1, kte
              z(k) = z_at_w(i,k,j) - z1
            END DO

            DO k = kts, kte
              phizz(k) = .001*z(k)
              aerext(k) = 0.
!               if(i.eq.its.and.j.eq.jts)print *,phizz(k),aerstd(k),kts,kte
!               print *,phizz(k),kts,kte,ite,jte
            END DO

!           IF (phizz(kte-1)>20.) THEN
!             CALL wrf_error_fatal ( 'phizz(kte-1)>20., set kl0 to 1')
!           END IF

            CALL trapez(zstd,aerstd,1,51,phizz,aerext,kts,kte,51,51,kte,kte)

            DO k = kts, kte
              aerwrf(i,k,j) = aerext(k)
!               if(i.eq.its)print *,k,i,j,aerext(k),phizz(k)
!               print *,k,i,j,aerext(k),phizz(k)
            END DO

10          CONTINUE
          END DO

20        CONTINUE
        END DO

      END SUBROUTINE photmad_init

! =========================================================================== 

      SUBROUTINE ps2str(kzmax, kz,zen,mu,rsfc,tauu,omu,gu, &
           dsdh, nid, delta,                               &
           fdr, fup, fdn, edr, eup, edn)

! ---------------------------------------------------------------------------- 
!   PURPOSE:                                                                  
!   Solve two-stream equations for multiple layers.  The subroutine is based  
!   on equations from:  Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989. 
!   It contains 9 two-stream methods to choose from.  A pseudo-spherical      
!   correction has also been added.                                           
! ---------------------------------------------------------------------------- 
!   PARAMETERS:                                                               
!   KZ      - INTEGER, number of specified altitude levels in the working (I) 
!             grid                                                            
!   ZEN     - REAL, solar zenith angle (degrees)                          (I) 
!   RSFC    - REAL, surface albedo at current wavelength                  (I) 
!   TAUU    - REAL, unscaled optical depth of each layer                  (I) 
!   OMU     - REAL, unscaled single scattering albedo of each layer       (I) 
!   GU      - REAL, unscaled asymmetry parameter of each layer            (I) 
!   DSDH    - REAL, slant path of direct beam through each layer crossed  (I) 
!             when travelling from the top of the atmosphere to layer i;      
!             DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                             
!   NID     - INTEGER, number of layers crossed by the direct beam when   (I) 
!             travelling from the top of the atmosphere to layer i;           
!             NID(i), i = 0..NZ-1                                             
!   DELTA   - LOGICAL, switch to use delta-scaling                        (I) 
!             .TRUE. -> apply delta-scaling                                   
!             .FALSE.-> do not apply delta-scaling                            
!   FDR     - REAL, contribution of the direct component to the total     (O) 
!             actinic flux at each altitude level                             
!   FUP     - REAL, contribution of the diffuse upwelling component to    (O) 
!             the total actinic flux at each altitude level                   
!   FDN     - REAL, contribution of the diffuse downwelling component to  (O) 
!             the total actinic flux at each altitude level                   
!   EDR     - REAL, contribution of the direct component to the total     (O) 
!             spectral irradiance at each altitude level                      
!   EUP     - REAL, contribution of the diffuse upwelling component to    (O) 
!             the total spectral irradiance at each altitude level            
!   EDN     - REAL, contribution of the diffuse downwelling component to  (O) 
!             the total spectral irradiance at each altitude level            
! ---------------------------------------------------------------------------- 

      IMPLICIT NONE

!******
! input:
!******
      INTEGER kz
      REAL zen, rsfc
      REAL tauu(kz), omu(kz), gu(kz)
      INTEGER kzmax
      REAL dsdh(0:kzmax,kzmax)
      INTEGER nid(0:kzmax)
      LOGICAL delta

!******
! output:
!******
      REAL fup(kz),fdn(kz),fdr(kz)
      REAL eup(kz),edn(kz),edr(kz)

!******
! local:
!******
      INTEGER nlevel
      REAL tausla(0:kz), tauc(0:kz)
      REAL mu2(0:kz), mu, sum

! internal coefficients and matrix

      REAL lam(kz),taun(kz),bgam(kz)
      REAL e1(kz),e2(kz),e3(kz),e4(kz)
      REAL cup(kz),cdn(kz),cuptn(kz),cdntn(kz)
      REAL mu1(kz)
      INTEGER row
      REAL a(2*kz),b(2*kz),d(2*kz),e(2*kz),y(2*kz)

!******
! other:
!******
      REAL pifs, fdn0
      REAL gi(kz), omi(kz), tempg
      REAL f, g, om
      REAL gam1, gam2, gam3, gam4

! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
! in delta-function, modified quadrature, hemispheric constant,
! Hybrid modified Eddington-delta function metods, p633,Table1.
! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440, 
! uncomment the following two lines and the appropriate statements further
! down.
!     REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
!    >     BETA1, BETAn, amu1, subd

      REAL expon, expon0, expon1, divisr, temp, up, dn
      REAL ssfc
      INTEGER nlayer, mrows, lev

      INTEGER i, j

! Some additional program constants:

      REAL pi
      PARAMETER (pi=3.1415926535898)
      REAL largest
      PARAMETER(largest=1.E+36)
      REAL precis
      PARAMETER(precis = 1.e-7)
      REAL eps
      PARAMETER (eps = 1.E-3)
!_______________________________________________________________________

! MU = cosine of solar zenith angle
! RSFC = surface albedo
! TAUU =  unscaled optical depth of each layer
! OMU  =  unscaled single scattering albedo
! GU   =  unscaled asymmetry factor
! KLEV = max dimension of number of layers in atmosphere
! NLAYER = number of layers in the atmosphere
! NLEVEL = nlayer + 1 = number of levels

! initial conditions:  pi*solar flux = 1;  diffuse incidence = 0

      nlevel = kz
      pifs = 1.      
      fdn0 = 0.
      nlayer = nlevel - 1
!      mu = COS(zen*pi/180.)

!************* compute coefficients for each layer:
! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
! EXPON0 = calculation of e when TAU is zero
! EXPON1 = calculation of e when TAU is TAUN
! CUP and CDN = calculation when TAU is zero
! CUPTN and CDNTN = calc. when TAU is TAUN
! DIVISR = prevents division by zero

        DO j = 0, kz
           tauc(j) = 0.
           tausla(j) = 0.
           mu2(j) = 1./SQRT(largest)
        END DO

       IF( .NOT. delta ) THEN
         DO i = 1, nlayer
           gi(i) = gu(i)
           omi(i) = omu(i)
           taun(i) = tauu(i)
         ENDDO
       ELSE 

! delta-scaling. Has to be done for delta-Eddington approximation, 
! delta discrete ordinate, Practical Improved Flux Method, delta function,
! and Hybrid modified Eddington-delta function methods approximations

         DO i = 1, nlayer
           f = gu(i)*gu(i)
           gi(i) = (gu(i) - f)/(1 - f)
           omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f)       
           taun(i) = (1 - omu(i)*f)*tauu(i)
         ENDDO
        END IF


! calculate slant optical depth at the top of the atmosphere when zen>90.
! in this case, higher altitude of the top layer is recommended which can 
! be easily changed in gridz.f.

         IF(zen .GT. 90.0) THEN
           IF(nid(0) .LT. 0) THEN
             tausla(0) = largest
           ELSE
             sum = 0.0
             DO j = 1, nid(0)
              sum = sum + 2.*taun(j)*dsdh(0,j)
             END DO
             tausla(0) = sum 
           END IF
         END IF
  

        DO 11, i = 1, nlayer

         g = gi(i)
         om = omi(i)
         tauc(i) = tauc(i-1) + taun(i)

! stay away from 1 by precision.  For g, also stay away from -1

         tempg = AMIN1(abs(g),1. - precis)
         g = SIGN(tempg,g)
         om = AMIN1(om,1.-precis)


! calculate slant optical depth
!              
          IF(nid(i) .LT. 0) THEN
            tausla(i) = largest
          ELSE
            sum = 0.0
            DO j = 1, MIN(nid(i),i)
               sum = sum + taun(j)*dsdh(i,j)
            ENDDO
            DO j = MIN(nid(i),i)+1,nid(i)
               sum = sum + 2.*taun(j)*dsdh(i,j)
            ENDDO
            tausla(i) = sum 
            IF(tausla(i) .EQ. tausla(i-1)) THEN
              mu2(i) = SQRT(largest)
            ELSE
              mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
              mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)), &
                           mu2(i) )
            END IF
          END IF

!** the following gamma equations are from pg 16,289, Table 1
!** save mu1 for each approx. for use in converting irradiance to actinic flux

! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):

        gam1 =  (7. - om*(4. + 3.*g))/4.
        gam2 = -(1. - om*(4. - 3.*g))/4.
        gam3 = (2. - 3.*g*mu)/4.
        gam4 = 1. - gam3
        mu1(i) = 0.5

! quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):

!          gam1 = 1.7320508*(2. - om*(1. + g))/2.
!          gam2 = 1.7320508*om*(1. - g)/2.
!          gam3 = (1. - 1.7320508*g*mu)/2.
!          gam4 = 1. - gam3
!          mu1(i) = 1./sqrt(3.)
         
! hemispheric mean (Toon et al., 1089, JGR, 94, 16287):

!          gam1 = 2. - om*(1. + g)
!          gam2 = om*(1. - g)
!          gam3 = (2. - g*mu)/4.
!          gam4 = 1. - gam3
!          mu1(i) = 0.5

! PIFM  (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
!         GAM1 = 0.25*(8. - OM*(5. + 3.*G))
!         GAM2 = 0.75*OM*(1.-G)
!         GAM3 = 0.25*(2.-3.*G*MU)
!         GAM4 = 1. - GAM3
!         mu1(i) = 0.5

! delta discrete ordinates  (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
!         GAM1 = 0.5*1.7320508*(2. - OM*(1. + G))
!         GAM2 = 0.5*1.7320508*OM*(1.-G)
!         GAM3 = 0.5*(1.-1.7320508*G*MU)
!         GAM4 = 1. - GAM3
!         mu1(i) = 1./sqrt(3.)

! Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
! in delta-function, modified quadrature, hemispheric constant,
! Hybrid modified Eddington-delta function metods, p633,Table1.
! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
!      YLM0 = 2.
!      YLM2 = -3.*G*MU
!      YLM4 = 0.875*G**3*MU*(5.*MU**2-3.)
!      YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4)
!     YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4
!    *+429.*MU**6)
!     YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4
!    *-25740.*MU**6+12155.*MU**8)
!     YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4
!    *+218790.*MU**6-230945.*MU**8+88179.*MU**10)
!      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
!      YLMS=0.25*YLMS
!      BETA0 = YLMS

!         amu1=1./1.7320508
!      YLM0 = 2.
!      YLM2 = -3.*G*amu1
!      YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.)
!      YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4)
!     YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4
!    *+429.*amu1**6)
!     YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4
!    *-25740.*amu1**6+12155.*amu1**8)
!     YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4
!    *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10)
!      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
!      YLMS=0.25*YLMS
!      BETA1 = YLMS

!         BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5
!    *-0.045776*G**7)


! Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
!         subd=4.*(1.-G*G*(1.-MU))
!         GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd
!         GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd
!         GAM3 = BETA0
!         GAM4 = 1. - GAM3
!         mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g)

!****
! delta function  (Meador, and Weaver, 1980, JAS, 37, 630):
!         GAM1 = (1. - OM*(1. - beta0))/MU
!         GAM2 = OM*BETA0/MU
!         GAM3 = BETA0
!         GAM4 = 1. - GAM3
!         mu1(i) = mu
!****
! modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
!         GAM1 = 1.7320508*(1. - OM*(1. - beta1))
!         GAM2 = 1.7320508*OM*beta1
!         GAM3 = BETA0
!         GAM4 = 1. - GAM3
!         mu1(i) = 1./sqrt(3.)

! hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
!         GAM1 = 2.*(1. - OM*(1. - betan))
!         GAM2 = 2.*OM*BETAn
!         GAM3 = BETA0
!         GAM4 = 1. - GAM3
!         mu1(i) = 0.5

!****

! lambda = pg 16,290 equation 21
! big gamma = pg 16,290 equation 22
! if gam2 = 0., then bgam = 0. 

         lam(i) = sqrt(gam1*gam1 - gam2*gam2)

         IF( gam2 .NE. 0.) THEN
            bgam(i) = (gam1 - lam(i))/gam2
         ELSE
            bgam(i) = 0.
         ENDIF

         expon = EXP(-lam(i)*taun(i))

! e1 - e4 = pg 16,292 equation 44
         
         e1(i) = 1. + bgam(i)*expon
         e2(i) = 1. - bgam(i)*expon
         e3(i) = bgam(i) + expon
         e4(i) = bgam(i) - expon

! the following sets up for the C equations 23, and 24
! found on page 16,290
! prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
! which is approx equiv to shifting MU by 0.5*EPS* (MU)**3

         expon0 = EXP(-tausla(i-1))
         expon1 = EXP(-tausla(i))
          
         divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i))
         temp = AMAX1(eps,abs(divisr))
         divisr = SIGN(temp,divisr)

         up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr
         dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr
         
! cup and cdn are when tau is equal to zero
! cuptn and cdntn are when tau is equal to taun

         cup(i) = up*expon0
         cdn(i) = dn*expon0
         cuptn(i) = up*expon1
         cdntn(i) = dn*expon1
 
   11 CONTINUE

!**************** set up matrix ******
! ssfc = pg 16,292 equation 37  where pi Fs is one (unity).

      ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs

! MROWS = the number of rows in the matrix

      mrows = 2*nlayer     
      
! the following are from pg 16,292  equations 39 - 43.
! set up first row of matrix:

      i = 1
      a(1) = 0.
      b(1) = e1(i)
      d(1) = -e2(i)
      e(1) = fdn0 - cdn(i)

      row=1

! set up odd rows 3 thru (MROWS - 1):
      i = 0
      DO 20, row = 3, mrows - 1, 2
         i = i + 1
         a(row) = e2(i)*e3(i) - e4(i)*e1(i)
         b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1)
         d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1)
         e(row) = e3(i)*(cup(i + 1) - cuptn(i)) + &
              e1(i)*(cdntn(i) - cdn(i + 1))
   20 CONTINUE

! set up even rows 2 thru (MROWS - 2): 

      i = 0
      DO 30, row = 2, mrows - 2, 2
         i = i + 1
         a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1)
         b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1)
         d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1)
         e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) - &
              (cdn(i + 1) - cdntn(i))*e4(i + 1)
   30 CONTINUE

! set up last row of matrix at MROWS:

      row = mrows
      i = nlayer
      
      a(row) = e1(i) - rsfc*e3(i)
      b(row) = e2(i) - rsfc*e4(i)
      d(row) = 0.
      e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)

! solve tri-diagonal matrix:

      CALL tridag(a, b, d, e, y, mrows)

!*** unfold solution of matrix, compute output fluxes:

      row = 1 
      lev = 1
      j = 1
      
! the following equations are from pg 16,291  equations 31 & 32

      fdr(lev) = EXP( -tausla(0) )
      edr(lev) = mu * fdr(lev)
      edn(lev) = fdn0
      eup(lev) =  y(row)*e3(j) - y(row + 1)*e4(j) + cup(j)
      fdn(lev) = edn(lev)/mu1(lev)
      fup(lev) = eup(lev)/mu1(lev)

      DO 60, lev = 2, nlayer + 1
         fdr(lev) = EXP(-tausla(lev-1))
         edr(lev) =  mu *fdr(lev)
         edn(lev) =  y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j)
         eup(lev) =  y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j)
         fdn(lev) = edn(lev)/mu1(j)
         fup(lev) = eup(lev)/mu1(j)

         row = row + 2
         j = j + 1
   60 CONTINUE
!_______________________________________________________________________

      END SUBROUTINE ps2str

! =========================================================================== 

      SUBROUTINE tridag(a,b,c,r,u,n)

!_______________________________________________________________________
! solves tridiagonal system.  From Numerical Recipies, p. 40
!_______________________________________________________________________

      IMPLICIT NONE

! input:
      INTEGER n
      REAL a, b, c, r
      DIMENSION a(n),b(n),c(n),r(n)

! output:
      REAL u
      DIMENSION u(n)

! local:
      INTEGER j

      REAL bet, gam
      DIMENSION gam(2*n)
!_______________________________________________________________________

      IF (b(1) .EQ. 0.) STOP 'subroutine tridag failed, 1001'
      bet   = b(1)
      u(1) = r(1)/bet
      DO 11, j = 2, n   
         gam(j) = c(j - 1)/bet
         bet = b(j) - a(j)*gam(j)
         IF (bet .EQ. 0.) STOP 'subroutine tridag failed, 2002'
         u(j) = (r(j) - a(j)*u(j - 1))/bet
   11 CONTINUE
      DO 12, j = n - 1, 1, -1  
         u(j) = u(j) - gam(j + 1)*u(j + 1)
   12 CONTINUE
!_______________________________________________________________________

      END SUBROUTINE tridag

! =========================================================================== 

      SUBROUTINE sphers(kz, nz, z, zen, dsdh, nid)

! ---------------------------------------------------------------------------- 
!   PURPOSE:                                                                  
!   Calculate slant path over vertical depth ds/dh in spherical geometry.     
!   Calculation is based on:  A.Dahlback, and K.Stamnes, A new spheric model  
!   for computing the radiation field available for photolysis and heating    
!   at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B)   
! ---------------------------------------------------------------------------- 
!   PARAMETERS:                                                               
!   NZ      - INTEGER, number of specified altitude levels in the working (I) 
!             grid                                                            
!   Z       - REAL, specified altitude working grid (km)                  (I) 
!   ZEN     - REAL, solar zenith angle (degrees)                          (I) 
!   DSDH    - REAL, slant path of direct beam through each layer crossed  (O) 
!             when travelling from the top of the atmosphere to layer i;      
!             DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                             
!   NID     - INTEGER, number of layers crossed by the direct beam when   (O) 
!             travelling from the top of the atmosphere to layer i;           
!             NID(i), i = 0..NZ-1                                             
! ---------------------------------------------------------------------------- 
!   EDIT HISTORY:                                                             
!   double precision fix for shallow layers - Julia Lee-Taylor Dec 2000       
! ---------------------------------------------------------------------------- 
!  This program is free software;  you can redistribute it and/or modify      
!  it under the terms of the GNU General Public License as published by the   
!  Free Software Foundation;  either version 2 of the license, or (at your    
!  option) any later version.                                                 
!  The TUV package is distributed in the hope that it will be useful, but     
!  WITHOUT ANY WARRANTY;  without even the implied warranty of MERCHANTIBI-   
!  LITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public      
!  License for more details.                                                  
!  To obtain a copy of the GNU General Public License, write to:              
!  Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.    
! ---------------------------------------------------------------------------- 
!  To contact the authors, please mail to:                                    
!  Sasha Madronich, NCAR/ACD, P.O.Box 3000, Boulder, CO, 80307-3000, USA  or  
!  send email to:  sasha@ucar.edu                                             
! ---------------------------------------------------------------------------- 
!  Copyright (C) 1994-2008   University Corporation for Atmospheric Research  
! ---------------------------------------------------------------------------- 

      IMPLICIT NONE

! input
      INTEGER kz, nz
      REAL zen, z(kz)

! output
      INTEGER nid(0:kz)
      REAL dsdh(0:kz,kz)

! more program constants
! pi:
      REAL pi
      PARAMETER(pi=3.1415926535898)
      REAL  dr
      PARAMETER ( dr = pi/180.)

! radius of the earth:
      REAL re
      REAL radius
      PARAMETER(radius=6.371E+3)

      REAL ze(kz)

! local 

      DOUBLE PRECISION zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
      INTEGER i, j, k
      INTEGER id

      INTEGER nlayer
      REAL zd(0:kz-1)

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

      zenrad = zen*dr

! number of layers:
      nlayer = nz - 1

! include the elevation above sea level to the radius of the earth:
      re = radius + z(1)
! correspondingly z changed to the elevation above earth surface:
      DO k = 1, nz
         ze(k) = z(k) - z(1)
      END DO

! inverse coordinate of z
      zd(0) = ze(nz)
      DO k = 1, nlayer
        zd(k) = ze(nz - k)
      END DO

! initialize dsdh(i,j), nid(i)
      DO i = 0, kz
       nid(i) = 0
       DO j = 1, kz
        dsdh(i,j) = 0.
       END DO
      END DO

! calculate ds/dh of every layer
      DO 100 i = 0, nlayer

        rpsinz = (re + zd(i)) * SIN(zenrad)
 
        IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
           nid(i) = -1
        ELSE


! Find index of layer in which the screening height lies
      
           id = i 
           IF( zen .GT. 90.0 ) THEN
              DO 10 j = 1, nlayer
                 IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND. &
                     (rpsinz .GE. ( zd(j) + re )) ) id = j
 10           CONTINUE
           END IF
 
           DO 20 j = 1, id

             sm = 1.0
             IF(j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0) &
                sm = -1.0
 
             rj = re + zd(j-1)
             rjp1 = re + zd(j)
 
             dhj = zd(j-1) - zd(j)
 
             ga = rj*rj - rpsinz*rpsinz
             gb = rjp1*rjp1 - rpsinz*rpsinz
             IF (ga .LT. 0.0) ga = 0.0
             IF (gb .LT. 0.0) gb = 0.0
 
             IF(id.GT.i .AND. j.EQ.id) THEN
                dsj = SQRT( ga )
             ELSE
                dsj = SQRT( ga ) - sm*SQRT( gb )
             END IF
             dsdh(i,j) = dsj / dhj
 20        CONTINUE
 
           nid(i) = id
 
        END IF

 100  CONTINUE

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

      END SUBROUTINE sphers

! =========================================================================== 

      SUBROUTINE airmas(kz, nz, dsdh, nid, cz, &
            vcol, scol)

! ---------------------------------------------------------------------------- 
!   PURPOSE:                                                                  
!   Calculate vertical and slant air columns, in spherical geometry, as a     
!   function of altitude.                                                     
! ---------------------------------------------------------------------------- 
!   PARAMETERS:                                                               
!   kZ      - INTEGER, number of specified altitude levels in the working (I) 
!             grid                                                            
!   DSDH    - REAL, slant path of direct beam through each layer crossed  (O) 
!             when travelling from the top of the atmosphere to layer i;      
!             DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                             
!   NID     - INTEGER, number of layers crossed by the direct beam when   (O) 
!             travelling from the top of the atmosphere to layer i;           
!             NID(i), i = 0..NZ-1                                             
!   VCOL    - REAL, output, vertical air column, molec cm-2, above level iz   
!   SCOL    - REAL, output, slant air column in direction of sun, above iz    
!             also in molec cm-2                                              
! ---------------------------------------------------------------------------- 

      IMPLICIT NONE

! Input:

      INTEGER kz, nz
      INTEGER nid(0:kz)
      REAL dsdh(0:kz,kz)
      REAL cz(nz)

! output: 

      REAL vcol(nz), scol(nz)

! internal:

      INTEGER id, j
      REAL sum, vsum
      REAL largest
      PARAMETER(largest=1.E+36)

! calculate vertical and slant column from each level:
! work downward

      vsum = 0.
      DO id = 0, nz - 1
         vsum = vsum + cz(nz-id)
         vcol(nz-id) = vsum
         sum = 0.
         IF(nid(id) .LT. 0) THEN
            sum = largest
         ELSE

! single pass layers:

            DO j = 1, MIN(nid(id), id)
               sum = sum + cz(nz-j)*dsdh(id,j)
            ENDDO

! double pass layers:

            DO j = MIN(nid(id),id)+1, nid(id)
               sum = sum + 2.*cz(nz-j)*dsdh(id,j)
            ENDDO

         ENDIF
         scol(nz - id) = sum

      ENDDO

      END SUBROUTINE airmas

    END MODULE module_phot_mad
