! 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 ! SAMcKeen 2/17/13 Update 296.309 to 305nm flux from SAO (K Chance et al., 2010) ! 2.5% increase in JO1D at 17.34deg SZA, ~1% at 58.8deg SZA, both at 500m AGL ! SAM 2/15/13 Dobson value set to 325., OMI - 5 days of eval. in May 2010 over LA ! 14.2% dobson reduction leads to 25% increase in surf. JO1D at 17deg SZA, ~50% at 80deg SZA ! ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation ! O3->O1D+O2 quantum yields taken from FTUV function fo3qy, based on ! Matsumi, Y., F. J. Comes, G. Hancock, A. Hofzumahaus, A. J. Hynes, ! M. Kawasaki, and A. R. Ravishankara, Quantum yields for production of O(1D) ! in the ultraviolet photolysis of ozone: Recommendation based on evaluation ! of laboratory data, J. Geophys. Res., 107, 10.1029/2001JD000510, 2002. ! Revised numbers give 23.3% JO1D increase at 17.34deg SZA, ~50% at 58.8deg, both at 500m AGL. ! ----------------------------------------------------------------------- ! 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 !liqy. change nreakj from 23 to 27, adding 4 species. INTEGER, PARAMETER :: nabv = 10, nj = 200, nreakj = 27 !liqy-20140908 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), jpl295(130), jpl218(130), sra(11,9), & !liqy-change xqy from xqy(130,23) to xqy(130,nreakj) srb(11,5), tabv(nabv), tstd(51), txs(130,nreakj),wl(130),xqy(130,nreakj), & !liqy-20140908 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./ ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation DATA (jpl295(kl),kl=1,130)/ 62.2, 57.5844, 52.5359, 47.7, 42.8638, 38.4591, & 34.9, 32.4255, 31.5, 32.6204, 36.3, 42.2149, 53.9, 69.3001, 90.3, 118., & 154., 199., 255., 322., 401., 490., 590., 693., 802., 908., 1001., 1080., & 1125., 1148., 1122., 1064., 968., 840., 698.001, 547., 406., 282., 184., & 113., 65.1, 37.2611, 26.2, 23.4, 20.1, 17.9, 15.5, 13.5, 12.2, 10.2, 9.24, & 7.95, 6.91, 6.25, 4.66, 2.82778, 1.378, 0.697, 0.32, 0.146, 0.0779, 0.0306, & 0.0136, 0.00694, 0.00305, 0.0013, 0.00085, 0.000572, 0.000542, 0.000668, & 0.000956, 0.00115, 0.00158, 0.00258, 0.00295, 0.00393, 0.00656, 0.00697, & 0.00882, 0.0137, 0.0165, 0.0185, 0.0218, 0.0366, 0.0367, 0.041, 0.0481, & 0.0754, 0.0813, 0.0816, 0.0908, 0.121, 0.16, 0.158, 0.166, 0.183, 0.219, & 0.267, 0.287, 0.295, 0.319, 0.337, 0.358, 0.398, 0.439, 0.467, 0.481, 0.464, & 0.446, 0.447, 0.476, 0.513, 0.514, 0.478, 0.438, 0.406, 0.382, 0.356, 0.327, & 0.297, 0.271, 0.245684, 0.21025, 0.17025, 0.13775, 0.112725, 0.087725, 0.07355, & 0.062075, 0.048525/ DATA (jpl218(kl),kl=1,130)/ 62.2, 57.5844, 52.5359, 47.7, 42.8638, 38.4535, & 34.4, 32.0245, 31.2, 32.4203, 36.2, 42.1149, 54.2, 69.6001, 90.6, 119., & 155., 201., 256., 323., 403., 492., 589., 692., 799., 905., 995., 1074., & 1116., 1136., 1105., 1047., 952., 823., 681.001, 531., 391., 271., 175., & 105., 59.4, 33.3104, 22.9, 20.6, 17.3, 15.6, 13.3, 11.5, 10.4, 8.5, 7.76, & 6.53, 5.62, 5.05, 3.67, 2.15, 0.9852, 0.483, 0.204, 0.0797, 0.0779, 0.0306, & 0.0136, 0.00694, 0.00305, 0.0013, 0.00085, 0.000572, 0.000542, 0.000668, & 0.000956, 0.00115, 0.00158, 0.00258, 0.00295, 0.00393, 0.00656, 0.00697, & 0.00882, 0.0137, 0.0165, 0.0185, 0.0218, 0.0366, 0.0367, 0.041, 0.0481, & 0.0754, 0.0813, 0.0816, 0.0908, 0.121, 0.16, 0.158, 0.166, 0.183, 0.219, & 0.267, 0.287, 0.295, 0.319, 0.337, 0.358, 0.398, 0.439, 0.467, 0.481, 0.464, & 0.446, 0.447, 0.476, 0.513, 0.514, 0.478, 0.438, 0.406, 0.382, 0.356, 0.327, & 0.297, 0.271, 0.245684, 0.21025, 0.17025, 0.13775, 0.112725, 0.087725, 0.07355, & 0.062075, 0.048525/ 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.790E+14, 2.990E+14, 1.003E+14, 9.294E+13, 1.024E+14, 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/ !liqy !!CL2 DATA (xs(kl,24),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, 2.930E-21, 5.100E-21,7.270E-21,& 1.212E-20, 1.870E-20, 2.564E-20, 3.932E-20, 5.408E-20,7.340E-20,& 9.791E-20, 1.223E-19, 1.388E-19, 1.454E-19, 1.520E-19,1.586E-19,& 1.652E-19, 1.718E-19, 1.784E-19, 1.850E-19, 1.902E-19,1.954E-19,& 2.006E-19, 2.058E-19, 2.162E-19, 2.370E-19, 2.460E-19,2.550E-19,& 2.450E-19, 2.350E-19, 2.115E-19, 1.880E-19, 1.600E-19,1.320E-19,& 1.080E-19, 8.400E-20, 6.700E-20, 5.000E-20, 3.950E-20,2.900E-20,& 2.350E-20, 1.800E-20, 1.550E-20, 1.300E-20, 1.130E-20,9.600E-21,& 8.450E-21, 7.300E-21, 6.350E-21, 5.400E-21, 4.600E-21,3.800E-21,& 3.200E-21, 2.600E-21, 2.100E-21, 1.600E-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/ DATA (xqy(kl,24),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, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 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.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/ !!HOCL DATA (xs(kl,25),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, 6.785E-20,6.071E-20,& 5.600E-20, 5.402E-20, 5.489E-20, 5.914E-20, 6.645E-20,7.748E-20,& 9.227E-20, 1.090E-19, 1.280E-19, 1.470E-19, 1.659E-19,1.828E-19,& 1.960E-19, 2.031E-19, 2.058E-19, 2.018E-19, 1.924E-19,1.783E-19,& 1.604E-19, 1.408E-19, 1.198E-19, 1.002E-19, 8.215E-20,6.768E-20,& 5.650E-20, 4.958E-20, 4.650E-20, 4.671E-20, 4.934E-20,5.330E-20,& 5.733E-20, 6.013E-20, 6.100E-20, 6.120E-20, 6.120E-20,6.120E-20,& 6.095E-20, 6.070E-20, 6.020E-20, 5.970E-20, 5.905E-20,5.840E-20,& 5.750E-20, 5.660E-20, 5.450E-20, 4.950E-20, 4.235E-20,3.500E-20,& 2.810E-20, 2.220E-20, 1.765E-20, 1.430E-20, 1.205E-20,1.060E-20,& 9.680E-21, 8.880E-21, 8.040E-21, 7.080E-21, 6.020E-21,4.910E-21,& 3.845E-21, 2.880E-21, 2.080E-21, 1.440E-21, 9.700E-22,6.300E-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/ DATA (xqy(kl,25),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 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.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/ !!CLNO2 DATA (xs(kl,26),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,& 3.105E-18, 3.173E-18, 3.242E-18, 3.330E-18, 3.418E-18,3.413E-18,& 3.314E-18, 3.171E-18, 2.956E-18, 2.706E-18, 2.392E-18,2.090E-18,& 1.814E-18, 1.591E-18, 1.385E-18, 1.236E-18, 1.100E-18,9.845E-19,& 8.750E-19, 7.679E-19, 6.597E-19, 5.606E-19, 4.626E-19,3.846E-19,& 3.159E-19, 2.629E-19, 2.268E-19, 1.976E-19, 1.782E-19,1.657E-19,& 1.556E-19, 1.462E-19, 1.399E-19, 1.373E-19, 1.348E-19,1.319E-19,& 1.290E-19, 1.261E-19, 1.232E-19, 1.203E-19, 1.171E-19,1.139E-19,& 1.106E-19, 1.074E-19, 1.010E-19, 8.809E-20, 7.246E-20,5.831E-20,& 4.600E-20, 3.572E-20, 2.734E-20, 2.077E-20, 1.564E-20,1.170E-20,& 8.726E-21, 6.501E-21, 4.836E-21, 3.593E-21, 2.687E-21,2.008E-21,& 1.498E-21, 1.125E-21, 8.471E-22, 6.419E-22, 4.858E-22,3.627E-22,& 2.777E-22, 2.168E-22, 1.668E-22, 1.234E-22, 9.395E-23,7.482E-23,& 6.059E-23, 4.981E-23, 3.794E-23, 3.169E-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/ DATA (xqy(kl,26),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/ !!FMCL DATA (xs(kl,27),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, 3.922E-20, 4.513E-20, 5.021E-20, 5.371E-20,5.541E-20,& 5.452E-20, 5.817E-20, 5.800E-20, 5.645E-20, 5.236E-20,4.550E-20,& 4.033E-20, 3.512E-20, 2.308E-20, 2.046E-20, 8.900E-21,8.214E-21,& 3.511E-21, 2.221E-21, 1.498E-21, 1.181E-21, 8.634E-22,6.538E-22,& 4.710E-22, 2.883E-22, 2.250E-22, 2.061E-22, 2.006E-22,1.790E-22,& 1.557E-22, 1.323E-22, 9.346E-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,0.000E+00,& 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00/ DATA (xqy(kl,27),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, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, & 1.000000, 1.000000, 1.000000, 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.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/ !liqy-20140911 ! 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, & !liqy ph_cl2,ph_hocl,ph_clno2,ph_fmcl, & !liqy-20140911 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, & !liqy ph_cl2,ph_hocl,ph_clno2,ph_fmcl, & !liqy-20140911 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 (config_flags%cu_diag==1) 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. dobsi=325. ! 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) !liqy ph_cl2(i,k,j) = phot1(23,k) ph_hocl(i,k,j) = phot1(24,k) ph_clno2(i,k,j) = phot1(25,k) ph_fmcl(i,k,j) = phot1(26,k) !liqy-20140911 ! 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)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 (zetanj) 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: !orig tau = zt(lev) - 230. !orig a = 0.9*(0.369+2.85E-4*tau+1.28E-5*tau*tau+2.57E-8*tau*tau*tau) !orig b = -0.575 + 5.59E-3*tau - 1.439E-5*tau*tau - 3.27E-8*tau*tau*tau !orig c = 0.9*(0.518+9.87E-4*tau-3.94E-5*tau*tau+3.91E-7*tau*tau*tau) !orig xl0 = 308.20 + 4.4871E-2*tau + 6.9380E-5*tau*tau - & !orig 2.5452E-6*tau*tau*tau DO 60 kl = kl0, kl1 xl = wl(kl) !orig qy(lev,kl,2) = a*atan(b*(xl-xl0)) + c !orig IF (qy(lev,kl,2)<0.) qy(lev,kl,2) = 0.0 !orig IF (qy(lev,kl,2)>0.9) qy(lev,kl,2) = 0.9 qy(lev,kl,2) = 0. kt = 0.695 * zt(lev) q1 = 1. q2 = exp( -825.518/kt ) if( xl <= 305. ) then qy(lev,kl,2) = .90 else if( xl > 305. .and. xl <= 328. ) then tau = zt(lev)/300. qy(lev,kl,2) = 0.0765 & + 0.8036* (q1/(q1+q2))*exp( -((304.225 - xl)/5.576)**4 ) & + 8.9061*(tau)**2 *(q2/(q1+q2))*exp( -((314.957 - xl)/6.601)**2 ) & + 0.1192*(tau)**1.5 *exp( -((310.737 - xl)/2.187)**2 ) else if( xl > 328. .and. xl <= 340. ) then qy(lev,kl,2) = 0.08 else if( xl > 340. ) then qy(lev,kl,2) = 0. end if 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 !orig DO 90 kl = kl0, kl1 DO 90 kl = 1, 130 DO 80 lev = 1, nlevel !orig IF (kl<33 .OR. kl>61) GO TO 70 !orig tdiffx = zt(lev) - 230. !orig s(lev,kl,2) = (so3tx(kl,1)+so3tx(kl,2)*tdiffx+so3tx(kl,3)*tdiffx* & !orig tdiffx)*1.0E-18 ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation s(lev,kl,2) = jpl218(kl) + & (jpl295(kl) - jpl218(kl))*(zt(lev) - 218.)/(295. - 218.) IF (zt(lev) .LE. 218.) s(lev,kl,2) = jpl218(kl) IF (zt(lev) .GE. 295.) s(lev,kl,2) = jpl295(kl) s(lev,kl,2)=s(lev,kl,2)*1.E-20 s(lev,kl,3) = s(lev,kl,2) !orig 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. !orig DO 160 kl = kl0, kl1 DO 160 kl = 1, 130 DO 150 lay = 1, nlayer !orig tdiffx = vt(lay) - 230. !orig ao3(lay,kl) = xs(kl,2) !orig IF (kl>=33 .AND. kl<=61) ao3(lay,kl) = (so3tx(kl,1)+so3tx(kl,2)* & !orig tdiffx+so3tx(kl,3)*tdiffx*tdiffx)*1.0E-18 ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation ao3(lay,kl) = jpl218(kl) + & (jpl295(kl) - jpl218(kl))*(vt(lay) - 218.)/(295. - 218.) IF (vt(lay) .LE. 218.) ao3(lay,kl) = jpl218(kl) IF (vt(lay) .GE. 295.) ao3(lay,kl) = jpl295(kl) ao3(lay,kl)=ao3(lay,kl)*1.E-20 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(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