* This file contains subroutines related to reading the
* absorption cross sections of gases that contribute to atmospheric transmission:
* Some of these subroutines are also called from rxn.f when loading photolysis cross sections
* for these same gases. It is possible to have different cross sections for 
* transmission and for photolysis, e.g. for ozone, Bass et al. could be used
* for transmission while Molina and Molina could be used for photolysis.  
* This flexibility can be useful but users should be aware.
* For xsections that are temperature dependent, caution should be used in passing the proper 
* temperature to the data routines.  Usually, transmission is for layers, TLAY(NZ-1), while
* photolysis is at levels, T(NZ).
* The following subroutines are her: 
*     rdo3xs
*       o3_mol
*       o3_rei
*       o3_bas
*       o3_wmo
*       o3_jpl
*     rdo2xs
*     rdno2xs
*       no2xs_d
*       no2xs_jpl94
*       no2xs_har
*     rdso2xs
*=============================================================================*

      SUBROUTINE rdo3xs(mabs, nz,t,nw,wl, xs)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read ozone molecular absorption cross section.  Re-grid data to match    =*
*=  specified wavelength working grid. Interpolate in temperature as needed  =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  MABS   - INTEGER, option for splicing different combinations of       (I)=*
*=           absorption cross secttions                                      =*
*=  NZ     - INTEGER, number of altitude levels or layers                 (I)=*
*=  T      - REAL, temperature of levels or layers                        (I)=*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid. In vacuum, nm                          =*
*=  XS     - REAL, molecular absoprtion cross section (cm^2) of O3 at     (O)=*
*=           each specified wavelength (WMO value at 273)                    =*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

* input: (altitude working grid)

      INTEGER iw, nw
      REAL wl(kw)

      INTEGER iz, nz
      REAL t(kz)

* internal

      INTEGER mabs
      REAL tc

* output:
* ozone absorption cross sections interpolated to 
*   working wavelength grid (iw)
*   working altitude grid (iz) for temperature of layer or level (specified in call)
* Units are cm2 molecule-1 in vacuum

      REAL xs(kz,kw)

* wavelength-interpolated values from different O3 data sources 
* also values of significant wavelengths, converted to vacuum

      REAL rei218(kw), rei228(kw), rei243(kw), rei295(kw)
      REAL v195, v345, v830

      REAL wmo203(kw), wmo273(kw)
      REAL v176, v850

      REAL jpl295(kw), jpl218(kw)
      REAL v186, v825

      REAL mol226(kw), mol263(kw), mol298(kw)
      REAL v185, v240, v350

      REAL c0(kw), c1(kw), c2(kw)
      REAL vb245, vb342

*_______________________________________________________________________
* read data from different sources
* rei = Reims group (Malicet et al., Brion et al.)
* jpl = JPL 2006 evaluation
* wmo = WMO 1985 O3 assessment
* mol = Molina and Molina
* bas = Bass et al.

      CALL o3_rei(nw,wl, rei218,rei228,rei243,rei295, v195,v345,v830)

      CALL o3_jpl(nw,wl, jpl218,jpl295, v186,v825)

      CALL o3_wmo(nw,wl, wmo203,wmo273, v176,v850)

      CALL o3_mol(nw,wl, mol226,mol263,mol298, v185,v240,v350)

      CALL o3_bas(nw,wl, c0,c1,c2, vb245,vb342)

****** option 1:

      IF(mabs. EQ. 1) THEN

* assign according to wavelength range:
*  175.439 - 185.185  1985WMO (203, 273 K)
*  185.185 - 195.00   2006JPL_O3 (218, 295 K)
*  195.00  - 345.00   Reims group (218, 228, 243, 295 K)
*  345.00  - 830.00   Reims group (295 K)
*  no extrapolations in temperature allowed

         DO 10 iw = 1, nw-1
         DO 20 iz = 1, nz

         IF(wl(iw) .LT. v185) THEN
            xs(iz,iw) = wmo203(iw) + 
     $           (wmo273(iw) - wmo203(iw))*(t(iz) - 203.)/(273. - 203.)
            IF (t(iz) .LE. 203.) xs(iz,iw) = wmo203(iw)
            IF (t(iz) .GE. 273.) xs(iz,iw) = wmo273(iw)
         ENDIF

         IF(wl(iw) .GE. v185 .AND. wl(iw) .LE. v195) THEN
            xs(iz,iw) = jpl218(iw) + 
     $           (jpl295(iw) - jpl218(iw))*(t(iz) - 218.)/(295. - 218.)
            IF (t(iz) .LE. 218.) xs(iz,iw) = jpl218(iw)
            IF (t(iz) .GE. 295.) xs(iz,iw) = jpl295(iw)
         ENDIF

         IF(wl(iw) .GE. v195 .AND. wl(iw) .LT. v345) THEN
            IF (t(iz) .GE. 218. .AND. t(iz) .LT. 228.) THEN
               xs(iz,iw) = rei218(iw) + 
     $              (t(iz)-218.)*(rei228(iw)-rei218(iw))/(228.-218.)
            ELSEIF (t(iz) .GE. 228. .AND. t(iz) .LT. 243.) THEN
               xs(iz,iw) = rei228(iw) +
     $              (t(iz)-228.)*(rei243(iw)-rei228(iw))/(243.-228.)
            ELSEIF (t(iz) .GE. 243. .AND. t(iz) .LT. 295.) THEN
               xs(iz,iw) = rei243(iw) +
     $              (t(iz)-243.)*(rei295(iw)-rei243(iw))/(295.-243.)
            ENDIF
            IF (t(iz) .LT. 218.) xs(iz,iw) = rei218(iw)
            IF (t(iz) .GE. 295.) xs(iz,iw) = rei295(iw)
         ENDIF

         IF(wl(iw) .GE. v345) THEN
            xs(iz,iw) = rei295(iw)
         ENDIF

 20      CONTINUE
 10      CONTINUE

      ELSEIF(mabs .EQ. 2) THEN

* use exclusively JPL-2006

         DO iw = 1, nw-1
         DO iz = 1, nz

            xs(iz,iw) = jpl218(iw) + 
     $           (jpl295(iw) - jpl218(iw))*(t(iz) - 218.)/(295. - 218.)
            IF (t(iz) .LE. 218.) xs(iz,iw) = jpl218(iw)
            IF (t(iz) .GE. 295.) xs(iz,iw) = jpl295(iw)

         ENDDO
         ENDDO

      ELSEIF(mabs .EQ. 3) THEN

* use exclusively Molina and Molina

         DO iw = 1, nw-1
         DO iz = 1, nz
            
            IF(wl(iw) .LT. v240) THEN
               xs(iz,iw) = mol226(iw) + 
     $              (t(iz)-226.)*(mol298(iw)-mol226(iw))/(298.-226.)
            ELSE
               IF(t(iz) .LT. 263.) THEN
                  xs(iz,iw) = mol226(iw) + 
     $                 (t(iz)-226.)*(mol263(iw)-mol226(iw))/(263.-226.)
               ELSE
                  xs(iz,iw) = mol263(iw) + 
     $                 (t(iz)-263.)*(mol298(iw)-mol263(iw))/(298.-263.)
               ENDIF
            ENDIF
            IF (t(iz) .LE. 226.) xs(iz,iw) = mol226(iw)
            IF (t(iz) .GE. 298.) xs(iz,iw) = mol298(iw)

         ENDDO
         ENDDO

      ELSEIF(mabs .EQ. 4) THEN

* use exclusively Bass et al.
* note limited wavelength range 245-342

         DO iw = 1, nw-1
         DO iz = 1, nz

            tc = t(iz) - 273.15
            xs(iz,iw) = c0(iw) + c1(iw)*tc + c2(iw)*tc*tc

         ENDDO
         ENDDO

      ELSE
         STOP 'mabs not set in rdxs.f'
      ENDIF

      RETURN
      END

*=============================================================================*

      SUBROUTINE o3_rei(nw,wl, 
     $     rei218,rei228,rei243,rei295, v195,v345,v830)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read and interpolate the O3 cross section from Reims group               =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  REI218 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
*=  REI228 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
*=  REI243 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
*=  REI295 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
*=  V195   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*=              e.g. start, stop, or other change                            =*
*=  V345   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*=  V830   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

*  input

      INTEGER nw, iw
      REAL wl(kw)

** internal

      INTEGER kdata
      PARAMETER (kdata = 70000)

      INTEGER n1, n2, n3, n4
      REAL x1(kdata), x2(kdata), x3(kdata), x4(kdata)
      REAL y1(kdata), y2(kdata), y3(kdata), y4(kdata)

      INTEGER i
      INTEGER ierr

* used for air-to-vacuum wavelength conversion

      REAL refrac, ri(kdata)
      EXTERNAL refrac

* output:

      REAL rei218(kw), rei228(kw), rei243(kw), rei295(kw)
      REAL v195, v345, v830

* data from the Reims group:
*=  For Hartley and Huggins bands, use temperature-dependent values from     =*
*=  Malicet et al., J. Atmos. Chem.  v.21, pp.263-273, 1995.                 =*
*=  over 345.01 - 830.00, use values from Brion, room temperature only

      OPEN(UNIT=kin,FILE='DATAE1/O3/1995Malicet_O3.txt',STATUS='old')
      DO i = 1, 2
         READ(kin,*)
      ENDDO
      n1 = 15001
      n2 = 15001
      n3 = 15001
      n4 = 15001
      DO i = 1, n1
         READ(kin,*) x1(i), y1(i), y2(i), y3(i), y4(i)
         x2(i) = x1(i)
         x3(i) = x1(i)
         x4(i) = x1(i)
      ENDDO
      CLOSE (kin)

*=  over 345.01 - 830.00, use values from Brion, room temperature only
* skip datum at 345.00 because already read in from 1995Malicet

      OPEN(UNIT=kin,FILE='DATAE1/O3/1998Brion_295.txt',STATUS='old')
      DO i = 1, 15
         READ(kin,*)
      ENDDO
      DO i = 1, 48515-15
         n1 = n1 + 1
         READ(kin,*) x1(n1), y1(n1)
      ENDDO
      CLOSE (kin)

      DO i = 1, n1
         ri(i) = refrac(x1(i), 2.45E19)
      ENDDO
      DO i = 1, n1
         x1(i) = x1(i) * ri(i)
      ENDDO

      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,               0.,0.)
      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,            1.e+38,0.)
      CALL inter2(nw,wl,rei295,n1,x1,y1,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - Reims 295K'
         STOP
      ENDIF

      DO i = 1, n2
         ri(i) = refrac(x2(i), 2.45E19)
      ENDDO
      DO i = 1, n2
         x2(i) = x2(i) * ri(i)
         x3(i) = x2(i)
         x4(i) = x2(i)
      ENDDO

      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,               0.,0.)
      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,            1.e+38,0.)
      CALL inter2(nw,wl,rei243,n2,x2,y2,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - Reims 243K'
         STOP
      ENDIF

      CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
      CALL addpnt(x3,y3,kdata,n3,               0.,0.)
      CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
      CALL addpnt(x3,y3,kdata,n3,            1.e+38,0.)
      CALL inter2(nw,wl,rei228,n3,x3,y3,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - Reims 228K'
         STOP
      ENDIF

      CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.)
      CALL addpnt(x4,y4,kdata,n4,               0.,0.)
      CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.)
      CALL addpnt(x4,y4,kdata,n4,            1.e+38,0.)
      CALL inter2(nw,wl,rei218,n4,x4,y4,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - Reims 218K'
         STOP
      ENDIF

* wavelength breaks must be converted to vacuum:

      v195 = 195.00 * refrac(195.00, 2.45E19)
      v345 = 345.00 * refrac(345.00, 2.45E19)
      v830 = 830.00 * refrac(830.00, 2.45E19)

      RETURN
      END

*=============================================================================*

      SUBROUTINE o3_wmo(nw,wl, wmo203,wmo273, v176,v850)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read and interpolate the O3 cross section                                =*
*=  data from WMO 85 Ozone Assessment                                        =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  WMO203 - REAL, cross section (cm^2) for O3 at 203K                    (O)=*
*=  WMO273 - REAL, cross section (cm^2) for O3 at 273K                    (O)=*
*=  V176   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*=              e.g. start, stop, or other change                            =*
*=  V850   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

*  input

      INTEGER nw, iw
      REAL wl(kw)

* internal

      INTEGER kdata
      PARAMETER (kdata = 200)

      INTEGER n1, n2
      REAL x1(kdata), x2(kdata)
      REAL y1(kdata), y2(kdata)

      INTEGER i, idum
      REAL a1, a2, dum
      INTEGER ierr

* used for air-to-vacuum wavelength conversion

      REAL refrac, ri(kdata)
      EXTERNAL refrac

* output

      REAL wmo203(kw), wmo273(kw)
      REAL v176, v850

*----------------------------------------------------------
* cross sections from WMO 1985 Ozone Assessment
* from 175.439 to 847.500 nm

      OPEN(UNIT=kin,FILE='DATAE1/wmo85',STATUS='old')
      DO i = 1, 3
         read(kin,*)
      ENDDO
      n1 = 158
      n2 = 158
      DO i = 1, n1
         READ(kin,*) idum, a1, a2, dum, dum, dum, y1(i), y2(i)
         x1(i) = (a1+a2)/2.
         x2(i) = (a1+a2)/2.
      ENDDO
      CLOSE (kin)

* convert wavelengths to vacuum

      DO i = 1, n1
         ri(i) = refrac(x1(i), 2.45E19)
      ENDDO
      DO i = 1, n1
         x1(i) = x1(i) * ri(i)
         x2(i) = x2(i) * ri(i)
      ENDDO

      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,               0.,0.)
      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,           1.e+38,0.)
      CALL inter2(nw,wl,wmo203,n1,x1,y1,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 cross section - WMO - 203K'
         STOP
      ENDIF

      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,               0.,0.)
      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,           1.e+38,0.)
      CALL inter2(nw,wl,wmo273,n2,x2,y2,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 cross section - WMO - 273K'
         STOP
      ENDIF

* wavelength breaks must be converted to vacuum:
      
      a1 = (175.438 + 176.991) / 2.
      v176 = a1 * refrac(a1,2.45E19)

      a1 = (847.5 + 852.5) / 2.
      v850 = a1 * refrac(a1, 2.45E19)

      RETURN
      END

*=============================================================================*

      SUBROUTINE o3_jpl(nw,wl, jpl218,jpl295, v186,v825)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read and interpolate the O3 cross section from JPL 2006                  =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  JPL218 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
*=  JPL295 - REAL, cross section (cm^2) for O3 at 295K                    (O)=*
*=  V186   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*=              e.g. start, stop, or other change                            =*
*=  V825   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

*  input

      INTEGER nw, iw
      REAL wl(kw)

* internal

      INTEGER kdata
      PARAMETER (kdata = 200)

      INTEGER n1, n2
      REAL x1(kdata), x2(kdata)
      REAL y1(kdata), y2(kdata)

      INTEGER i
      REAL dum
      INTEGER ierr

* used for air-to-vacuum wavelength conversion

      REAL refrac, ri(kdata)
      EXTERNAL refrac

* output

      REAL jpl295(kw), jpl218(kw)
      REAL v186, v825

***********

      OPEN(UNIT=kin,FILE='DATAE1/O3/2006JPL_O3.txt',STATUS='old')
      DO i = 1, 2
         read(kin,*)
      ENDDO
      n1 = 167
      n2 = 167
      DO i = 1, n1
         READ(kin,*) dum, dum, x1(i), y1(i), y2(i)
         y1(i) = y1(i) * 1.e-20
         y2(i) = y2(i) * 1.e-20
      ENDDO
      CLOSE (kin)

* convert wavelengths to vacuum

      DO i = 1, n1
         ri(i) = refrac(x1(i), 2.45E19)
      ENDDO
      DO i = 1, n1
         x1(i) = x1(i) * ri(i)
         x2(i) = x1(i)
      ENDDO

      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,               0.,0.)
      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,           1.e+38,0.)
      CALL inter2(nw,wl,jpl295,n1,x1,y1,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 cross section - WMO - 295K'
         STOP
      ENDIF

      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,               0.,0.)
      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,           1.e+38,0.)
      CALL inter2(nw,wl,jpl218,n2,x2,y2,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 cross section - WMO - 218K'
         STOP
      ENDIF

* wavelength breaks must be converted to vacuum:

      v186 = 186.051 * refrac(186.051, 2.45E19)
      v825 = 825.    * refrac(825.   , 2.45E19)


      RETURN
      END


*=============================================================================*

      SUBROUTINE o3_mol(nw,wl, mol226,mol263,mol298, v185,v240,v350)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read and interpolate the O3 cross section from Molina and Molina 1986    =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  MOL226 - REAL, cross section (cm^2) for O3 at 226 K                   (O)=*
*=  MOL263 - REAL, cross section (cm^2) for O3 at 263 K                   (O)=*
*=  MOL298 - REAL, cross section (cm^2) for O3 at 298 K                   (O)=*
*=  V185   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*=              e.g. start, stop, or other change                            =*
*=  V240   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*=  V350   - REAL, exact wavelength in vacuum for data breaks             (O)=*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

*  input

      INTEGER nw, iw
      REAL wl(kw)

* internal

      INTEGER i
      INTEGER ierr

      INTEGER kdata
      PARAMETER (kdata = 335)
      INTEGER n1, n2, n3
      REAL x1(kdata), x2(kdata), x3(kdata)
      REAL y1(kdata), y2(kdata), y3(kdata)

* used for air-to-vacuum wavelength conversion

      REAL refrac, ri(kdata)
      EXTERNAL refrac

* output

      REAL mol226(kw), mol263(kw), mol298(kw)
      REAL v185, v240, v350

*----------------------------------------------------------

      OPEN(UNIT=kin,FILE='DATAE1/O3/1986Molina.txt',STATUS='old')
      DO i = 1, 10
         READ(kin,*)
      ENDDO
      n1 = 0
      n2 = 0
      n3 = 0
      DO i = 1, 121-10
         n1 = n1 + 1
         n3 = n3 + 1
         READ(kin,*) x1(n1), y1(n1),  y3(n3)
         x3(n3) = x1(n1)
      ENDDO
      DO i = 1, 341-122
         n1 = n1 + 1
         n2 = n2 + 1
         n3 = n3 + 1
         READ(kin,*) x1(n1), y1(n1), y2(n2), y3(n3)
         x2(n2) = x1(n1)
         x3(n3) = x1(n1)
      ENDDO
      CLOSE (kin)

* convert all wavelengths from air to vacuum

      DO i = 1, n1
         ri(i) = refrac(x1(i), 2.45E19)
      ENDDO
      DO i = 1, n1
         x1(i) = x1(i) * ri(i)
      ENDDO

      DO i = 1, n2
         ri(i) = refrac(x2(i), 2.45E19)
      ENDDO
      DO i = 1, n2
         x2(i) = x2(i) * ri(i)
      ENDDO

      DO i = 1, n3
         ri(i) = refrac(x3(i), 2.45E19)
      ENDDO
      DO i = 1, n3
         x3(i) = x3(i) * ri(i)
      ENDDO

* convert wavelength breaks from air to vacuum

      v185 = 185.  * refrac(185. , 2.45E19)
      v240 = 240.5 * refrac(240.5, 2.45E19)
      v350 = 350.  * refrac(350. , 2.45E19)

* interpolate to working grid

      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,               0.,0.)
      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,            1.e+38,0.)
      CALL inter2(nw,wl,mol226,n1,x1,y1,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - 226K Molina'
         STOP
      ENDIF

      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,               0.,0.)
      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,            1.e+38,0.)
      CALL inter2(nw,wl,mol263,n2,x2,y2,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - 263K Molina'
         STOP
      ENDIF

      CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
      CALL addpnt(x3,y3,kdata,n3,               0.,0.)
      CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
      CALL addpnt(x3,y3,kdata,n3,            1.e+38,0.)
      CALL inter2(nw,wl,mol298,n3,x3,y3,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - 298K Molina'
         STOP
      ENDIF

      RETURN
      END

*=============================================================================*

      SUBROUTINE o3_bas(nw,wl, c0,c1,c2, vb245,vb342)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read and interpolate the O3 cross section from Bass 1985                 =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  c0     - REAL, coefficint for polynomial fit to cross section (cm^2)  (O)=*
*=  c1     - REAL, coefficint for polynomial fit to cross section (cm^2)  (O)=*
*=  c2     - REAL, coefficint for polynomial fit to cross section (cm^2)  (O)=*
*=  Vb245   - REAL, exact wavelength in vacuum for data breaks            (O)=*
*=              e.g. start, stop, or other change                            =*
*=  Vb342   - REAL, exact wavelength in vacuum for data breaks            (O)=*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

* input:

      INTEGER nw, iw
      REAL wl(kw)

* internal:

      INTEGER kdata
      PARAMETER (kdata = 2000)

      INTEGER i
      INTEGER ierr

      INTEGER n1, n2, n3
      REAL x1(kdata), x2(kdata), x3(kdata)
      REAL y1(kdata), y2(kdata), y3(kdata)

* used for air-to-vacuum wavelength conversion

      REAL refrac, ri(kdata)
      EXTERNAL refrac

* output:

      REAL c0(kw), c1(kw), c2(kw)
      REAL vb245, vb342

*******************

      OPEN(UNIT=kin,FILE='DATAE1/O3/1985Bass_O3.txt',STATUS='old')
      DO i = 1, 8
         READ(kin,*)
      ENDDO
      n1 = 1915
      n2 = 1915
      n3 = 1915
      DO i = 1, n1
         READ(kin,*) x1(i), y1(i), y2(i), y3(i)
         y1(i) = 1.e-20 * y1(i)
         y2(i) = 1.e-20 * y2(i)
         y3(i) = 1.e-20 * y3(i)
      ENDDO
      CLOSE (kin)

* convert all wavelengths from air to vacuum

      DO i = 1, n1
         ri(i) = refrac(x1(i), 2.45E19)
      ENDDO
      DO i = 1, n1
         x1(i) = x1(i) * ri(i)
         x2(i) = x1(i)
         x3(i) = x1(i)
      ENDDO

* convert wavelength breaks to vacuum

      vb245 = 245.018 * refrac(245.018, 2.45E19)
      vb342 = 341.981 * refrac(341.981, 2.45E19)

* interpolate to working grid

      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,               0.,0.)
      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,            1.e+38,0.)
      CALL inter2(nw,wl,c0,n1,x1,y1,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - c0 Bass'
         STOP
      ENDIF

      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,               0.,0.)
      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,            1.e+38,0.)
      CALL inter2(nw,wl,c1,n2,x2,y2,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - c1 Bass'
         STOP
      ENDIF

      CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
      CALL addpnt(x3,y3,kdata,n3,               0.,0.)
      CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
      CALL addpnt(x3,y3,kdata,n3,            1.e+38,0.)
      CALL inter2(nw,wl,c2,n3,x3,y3,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O3 xsect - c2 Bass'
         STOP
      ENDIF

      RETURN
      END

*=============================================================================*

      SUBROUTINE rdo2xs(nw,wl,o2xs1)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Compute equivalent O2 cross section, except                              =*
*=  the SR bands and the Lyman-alpha line.                                   =*
*-----------------------------------------------------------------------------* 
*=  PARAMETERS:                                   
*=  NW      - INTEGER, number of specified intervals + 1 in working       (I)=*
*=            wavelength grid                                                =*
*=  WL      - REAL, vector of lower limits of wavelength intervals in     (I)=*
*=            working wavelength grid           
*=            vertical layer at each specified wavelength                    =*
*=  O2XS1   - REAL, O2 molecular absorption cross section                    =*
*=
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

* Input

      INTEGER nw
      REAL wl(kw)

* Output O2 xsect, temporary, will be over-written in Lyman-alpha and 
*   Schumann-Runge wavelength bands.

      REAL o2xs1(kw)

* Internal

      INTEGER i, n, kdata
      PARAMETER (kdata = 200)
      REAL x1(kdata), y1(kdata)
      REAL x, y
      INTEGER ierr

*-----------------------------------------------------

* Read O2 absorption cross section data:
*  116.65 to 203.05 nm = from Brasseur and Solomon 1986
*  205 to 240 nm = Yoshino et al. 1988

* Note that subroutine la_srb.f will over-write values in the spectral regions
*   corresponding to:
* - Lyman-alpha (LA: 121.4-121.9 nm, Chabrillat and Kockarts parameterization) 
* - Schumann-Runge bands (SRB: 174.4-205.8 nm, Koppers parameteriaztion)

      n = 0

      OPEN(UNIT=kin,FILE='DATAE1/O2/O2_brasseur.abs')
      DO i = 1, 7
         READ(kin,*)
      ENDDO
      DO i = 1, 78
         READ(kin,*) x, y
         IF (x .LE. 204.) THEN
            n = n + 1
            x1(n) = x
            y1(n) = y
         ENDIF
      ENDDO
      CLOSE(kin)

      OPEN(UNIT=kin,FILE='DATAE1/O2/O2_yoshino.abs',STATUS='old')
      DO i = 1, 8
         READ(kin,*)
      ENDDO
      DO i = 1, 36
         n = n + 1
         READ(kin,*) x, y
         y1(n) = y*1.E-24
         x1(n) = x
      END DO
      CLOSE (kin)

* Add termination points and interpolate onto the 
*  user grid (set in subroutine gridw):

      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),y1(1))
      CALL addpnt(x1,y1,kdata,n,0.               ,y1(1))
      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n,              1.E+38,0.)
      CALL inter2(nw,wl,o2xs1, n,x1,y1, ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, 'O2 -> O + O'
         STOP
      ENDIF

*------------------------------------------------------

      RETURN
      END

*=============================================================================*

      SUBROUTINE rdno2xs(nz,tlay,nw,wl,no2xs)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read NO2 molecular absorption cross section.  Re-grid data to match      =*
*=  specified wavelength working grid.                                       =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  NO2XS  - REAL, molecular absoprtion cross section (cm^2) of NO2 at    (O)=*
*=           each specified wavelength                                       =*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

* input: (altitude working grid)

      INTEGER nz
      REAL tlay(kz)

      INTEGER nw
      REAL wl(kw)

      INTEGER mabs

* output:

      REAL no2xs(kz,kw)

*_______________________________________________________________________

* options for NO2 cross section:
* 1 = Davidson et al. (1988), indepedent of T
* 2 = JPL 1994 (same as JPL 1997, JPL 2002)
* 3 = Harder et al.
* 4 = JPL 2006, interpolating between midpoints of bins
* 5 = JPL 2006, bin-to-bin interpolation

      mabs = 4

      IF (mabs. EQ. 1) CALL no2xs_d(nz,tlay,nw,wl, no2xs)
      IF (mabs .EQ. 2) CALL no2xs_jpl94(nz,tlay,nw,wl, no2xs)
      IF (mabs .EQ. 3) CALL no2xs_har(nz,tlay,nw,wl, no2xs)
      IF (mabs .EQ. 4) CALL no2xs_jpl06a(nz,tlay,nw,wl, no2xs)
      IF (mabs .EQ. 5) CALL no2xs_jpl06b(nz,tlay,nw,wl, no2xs)

*_______________________________________________________________________

      RETURN
      END

*=============================================================================*

      SUBROUTINE no2xs_d(nz,t,nw,wl, no2xs)

      IMPLICIT NONE
      INCLUDE 'params'

* read and interpolate NO2 xs from Davidson et al. (1988).

* input:

      INTEGER nz, iz
      REAL t(nz)

      INTEGER nw, iw
      REAL wl(nw)

* output:

      REAL no2xs(kz,kw)

* local:

      INTEGER kdata
      PARAMETER (kdata=1000)
      REAL x1(kdata)
      REAL y1(kdata)
      REAL yg(kw)
      REAL dum
      INTEGER ierr
      INTEGER i, n, idum
      CHARACTER*40 fil

************* NO2 absorption cross sections
*     measurements by:
* Davidson, J. A., C. A. Cantrell, A. H. McDaniel, R. E. Shetter,
* S. Madronich, and J. G. Calvert, Visible-ultraviolet absorption
* cross sections for NO2 as a function of temperature, J. Geophys.
* Res., 93, 7105-7112, 1988.
*  Values at 273K from 263.8 to 648.8 nm in approximately 0.5 nm intervals

      fil = 'DATAE1/NO2/NO2_ncar_00.abs'
      OPEN(UNIT=kin,FILE=fil,STATUS='old')
      n = 750
      DO i = 1, n
         READ(kin,*) x1(i), y1(i), dum, dum, idum
      ENDDO
      CLOSE(kin)

* interpolate to wavelength grid

      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n,          0.,0.)
      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, fil
         STOP
      ENDIF

* assign, same at all altitudes (no temperature effect)     

      DO iz = 1, nz
         DO iw = 1, nw-1
            no2xs(iz,iw) = yg(iw)
         ENDDO
      ENDDO

*_______________________________________________________________________

      RETURN
      END

*=============================================================================*

      SUBROUTINE no2xs_jpl94(nz,t,nw,wl, no2xs)

      IMPLICIT NONE
      INCLUDE 'params'

* read and interpolate NO2 xs from JPL 1994

* input:

      INTEGER nz, iz
      REAL t(nz)

      INTEGER nw, iw
      REAL wl(nw)

* output:

      REAL no2xs(kz,kw)

* local:

      INTEGER kdata
      PARAMETER (kdata=100)
      INTEGER i, idum, n, n1
      REAL x1(kdata), x2(kdata), x3(kdata)
      REAL y1(kdata), y2(kdata), y3(kdata)
      REAL dum
      REAL yg1(nw), yg2(nw)

* cross section data from JPL 94 recommendation
* JPL 97 and JPL 2002 recommendations are identical

      OPEN(UNIT=kin,FILE='DATAE1/NO2/NO2_jpl94.abs',STATUS='old')
      READ(kin,*) idum, n
      DO i = 1, idum-2
         READ(kin,*)
      ENDDO 

* read in wavelength bins, cross section at T0 and temperature correction
* coefficient a;  see input file for details.
* data need to be scaled to total area per bin so that they can be used with
* inter3

      DO i = 1, n
         READ(kin,*) x1(i), x3(i), y1(i), dum, y2(i)
         y1(i) = (x3(i)-x1(i)) * y1(i)*1.E-20
         y2(i) = (x3(i)-x1(i)) * y2(i)*1.E-22
         x2(i) = x1(i) 
      ENDDO
      CLOSE(kin)

      x1(n+1) = x3(n)
      x2(n+1) = x3(n)
      n = n+1
      n1 = n

      CALL inter3(nw,wl,yg1,n,x1,y1,0)
      CALL inter3(nw,wl,yg2,n1,x2,y2,0)

* yg1, yg2 are per nm, so rescale by bin widths

      DO iw = 1, nw-1
         yg1(iw) = yg1(iw)/(wl(iw+1)-wl(iw))
         yg2(iw) = yg2(iw)/(wl(iw+1)-wl(iw))
      ENDDO

      DO iw = 1, nw-1
         DO iz = 1, nz
            no2xs(iz,iw) = yg1(iw) + yg2(iw)*(t(iz)-273.15)
         ENDDO
      ENDDO 

      RETURN
      END

*=============================================================================*

      SUBROUTINE no2xs_har(nz,t,nw,wl, no2xs)

      IMPLICIT NONE
      INCLUDE 'params'

* read and interpolate NO2 xs from Harder et al.

* input:

      INTEGER nz, iz
      REAL t(nz)

      INTEGER nw, iw
      REAL wl(nw)

* output:

      REAL no2xs(kz,kw)

* local:

      INTEGER kdata
      PARAMETER (kdata=150)
      INTEGER i, n, idum, ierr
      REAL x1(kdata), y1(kdata)
      REAL yg1(kw)

***

      OPEN(UNIT=kin,FILE='DATAE1/NO2/NO2_Har.abs',status='old')
      DO i = 1, 9
         READ(kin,*)
      ENDDO
      n = 135
      DO i = 1, n
         READ(kin,*) idum, y1(i)
         x1(i) = FLOAT(idum)
      ENDDO

      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),y1(1))
      CALL addpnt(x1,y1,kdata,n,               0.,y1(1))
      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),   0.)
      CALL addpnt(x1,y1,kdata,n,           1.e+38,   0.)
      CALL inter2(nw,wl,yg1,n,x1,y1,ierr)

      DO iw = 1, nw-1
         DO i = 1, nz
            no2xs(i,iw) = yg1(iw)
         ENDDO
      ENDDO 

      RETURN
      END

*=============================================================================*

      SUBROUTINE no2xs_jpl06a(nz,t,nw,wl, no2xs)

* read and interpolate NO2 xs from JPL2006

      IMPLICIT NONE
      INCLUDE 'params'

* input:

      INTEGER nz, iz
      REAL t(nz)

      INTEGER nw, iw
      REAL wl(nw)

* output:

      REAL no2xs(kz,kw)

* local

      INTEGER kdata
      PARAMETER (kdata=100)
      INTEGER i, n1, n2, ierr
      REAL x1(kdata), x2(kdata), y1(kdata), y2(kdata)
      REAL dum1, dum2
      REAL yg1(kw), yg2(kw)

* NO2 absorption cross section from JPL2006
* with interpolation of bin midpoints

      OPEN(UNIT=kin,FILE='DATAE1/NO2/NO2_jpl2006.abs',STATUS='old')
      DO i = 1, 3
         READ(kin,*)
      ENDDO 
      n1 = 81
      DO i = 1, n1
         READ(kin,*) dum1, dum2, y1(i), y2(i)
         x1(i) = 0.5 * (dum1 + dum2)
         x2(i) = x1(i) 
         y1(i) = y1(i)*1.E-20
         y2(i) = y2(i)*1.E-20
      ENDDO
      CLOSE(kin)
      n2 = n1

      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n1,               0.,0.)
      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),   0.)
      CALL addpnt(x1,y1,kdata,n1,            1.e+38,   0.)
      CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
      
      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
      CALL addpnt(x2,y2,kdata,n2,               0.,0.)
      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),   0.)
      CALL addpnt(x2,y2,kdata,n2,            1.e+38,   0.)
      CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
      
      DO iw = 1, nw-1
         DO iz = 1, nz
            no2xs(iz,iw) = yg1(iw) + 
     $           (yg2(iw)-yg1(iw))*(t(iz)-220.)/74.
         ENDDO
      ENDDO 

      RETURN
      END

*=============================================================================*

      SUBROUTINE no2xs_jpl06b(nz,t,nw,wl, no2xs)

* read and interpolate NO2 xs from Harder et al.

      IMPLICIT NONE
      INCLUDE 'params'

* input:

      INTEGER nz, iz
      REAL t(nz)

      INTEGER nw, iw
      REAL wl(nw)

* output:

      REAL no2xs(kz,kw)

* local

      INTEGER kdata
      PARAMETER (kdata=100)
      INTEGER i, n, n1, n2, ierr
      REAL x1(kdata), x2(kdata), y1(kdata), y2(kdata)
      REAL x3(kdata), y3(kdata)
      REAL dum1, dum2
      REAL yg1(kw), yg2(kw)

      OPEN(UNIT=kin,FILE='DATAE1/NO2/NO2_jpl2006.abs',STATUS='old')
      DO i = 1, 3
         READ(kin,*)
      ENDDO 
      n = 81
      do i = 1, n
         read(kin,*) x1(i), x3(i), y1(i), y2(i)
         y1(i) = (x3(i)-x1(i)) * y1(i)*1.E-20
         y2(i) = (x3(i)-x1(i)) * y2(i)*1.E-20
         x2(i) = x1(i) 
      ENDDO
      CLOSE(kin)
         
      x1(n+1) = x3(n)
      x2(n+1) = x3(n)
      n = n+1
      n1 = n

      CALL inter3(nw,wl,yg1,n,x1,y1,0)
      CALL inter3(nw,wl,yg2,n1,x2,y2,0)

* yg1, yg2 are per nm, so rescale by bin widths

      DO iw = 1, nw-1
         yg1(iw) = yg1(iw)/(wl(iw+1)-wl(iw))
         yg2(iw) = yg2(iw)/(wl(iw+1)-wl(iw))
      ENDDO

      DO iw = 1, nw-1
         DO iz = 1, nz
            no2xs(iz,iw) = yg1(iw) + 
     $           (yg2(iw)-yg1(iw))*(t(iz)-220.)/74.
         ENDDO

      ENDDO 

      RETURN
      END

*=============================================================================*

      SUBROUTINE rdso2xs(nw,wl,so2xs)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Read SO2 molecular absorption cross section.  Re-grid data to match      =*
*=  specified wavelength working grid.                                       =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  SO2XS  - REAL, molecular absoprtion cross section (cm^2) of SO2 at    (O)=*
*=           each specified wavelength                                       =*
*-----------------------------------------------------------------------------*
*=  EDIT HISTORY:                                                            =*
*=  02/97  Changed offset for grid-end interpolation to relative number      =*
*=         (x * (1 +- deltax)                                                =*
*-----------------------------------------------------------------------------*
*= 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                                            =*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

      INTEGER kdata
      PARAMETER(kdata=1000)

* input: (altitude working grid)
      INTEGER nw
      REAL wl(kw)

* output:

      REAL so2xs(kw)

* local:
      REAL x1(kdata)
      REAL y1(kdata)
      REAL yg(kw)
      REAL dum
      INTEGER ierr
      INTEGER i, l, n, idum
      CHARACTER*40 fil
*_______________________________________________________________________

************* absorption cross sections:
* SO2 absorption cross sections from J. Quant. Spectrosc. Radiat. Transfer
* 37, 165-182, 1987, T. J. McGee and J. Burris Jr.
* Angstrom vs. cm2/molecule, value at 221 K

      fil = 'DATA/McGee87'
      OPEN(UNIT=kin,FILE='DATAE1/SO2/SO2xs.all',STATUS='old')
      DO 11, i = 1,3 
         read(kin,*)
   11 CONTINUE
c      n = 681 
      n = 704 
      DO 12, i = 1, n
         READ(kin,*) x1(i), y1(i)
         x1(i) = x1(i)/10.
   12 CONTINUE
      CLOSE (kin)

      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
      CALL addpnt(x1,y1,kdata,n,          0.,0.)
      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, fil
         STOP
      ENDIF
      
      DO 13, l = 1, nw-1
         so2xs(l) = yg(l)
   13 CONTINUE

*_______________________________________________________________________

      RETURN
      END

*=============================================================================*

