c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine coffs(wvl_c, delwvl, a0, a1, a2, b, p0, p2, w, nc, ns)
c
c     -- Coeff/Optical coefficients and weights for SBUV sub-channels
c        /N. Nath/April '01 (Equidistant position version)
c
c     -- Function:  This is an initial set-up routine that assigns 
c        component wavelengths within each SBUV channel, determines 
c        ozone absorption coefficients, Rayleigh scattering 
c        coefficients, and phase function coefficients at those 
c        wavelengths, and also determines the weights of the components
c        for computation of the average channel intensity.  The 
c        components/sub-channels are chosen to be equidistant (with 
c        intrinsic equal weights), and the weights are products of the 
c        triangular slit response function and the variation in the 
c        solar flux.
c
            implicit none
c
c     -- Input & Output Variables:
c
            real        wvl_c(nc), delwvl
                              !center wavelengths (nm), half-width of 
                              !  triangular slit function (nm) (input)

            real        a0(nc,ns), a1(nc,ns), a2(nc,ns), b(nc,ns), 
     1                  p0(nc,ns), p2(nc,ns), w(nc,ns)
                              !a0-a2: o3 absorption coefficients (a0 in 
                              !  (atm-cm)**-1), b: Rayleigh scattering
                              !  coefficient (atm**-1), p0-p2: phase 
                              !  function coefficients, w: weights 
                              !  (output)

            integer     nc, ns!number of channels & subchannels (input)
                              !  (nc max = 20, ns max = 16)  
c
c     -- Local Variables:
c
            integer     nref, ncl, nsl
            parameter   (nref = 1410, ncl = 20, nsl = 16)
c
            real        wvl_r(nref), a0_r(nref), a1_r(nref), a2_r(nref),
     1                  b_r(nref), rho_r(nref), flux_r(nref) 
                              !reference data: wave lengths, o3 
                              !  absorption coeffs, Rayleigh scattering 
                              !  coeffs, molecular anisotropy, & solar 
                              !  flux

            real        wvl(ncl,nsl), flux(ncl,nsl), rho
                              !subchannel wavelength, flux & anisotropy

            real        xs(nsl)!relative sub-channel positions

            real        sum, c, c1

            integer     i, ic, is, l, u, m, lun_alfa
c
      include 'flnams.com'
      character*60 cofile
      equivalence (fname(15),cofile)
c
c           Read tabulated reference data
c
      call get_lun(lun_alfa)
      if (lun_alfa .lt. 0) then
        print*,'error in getting a logical unit number for the sbuv',
     1          ' coefficients data set'
        stop
      endif
c
c     lun_alfa=67
c     cofile='/misc/cgw/source/v8alpha/solar_bass.dat'
      print*,'opening ',cofile,lun_alfa
      open (unit=lun_alfa,file=cofile,status='old',form='formatted')
      do i = 1,nref
         read (lun_alfa, *) wvl_r(i), a0_r(i), a1_r(i), a2_r(i), b_r(i),
     1     rho_r(i), flux_r(i)
         wvl_r(i) = wvl_r(i)/10.0
                              !convert wvl in nm
      end do
      close (lun_alfa)
c
      do is = 1,ns
         xs(is) = -1.0 + (1.0/ns)*(2*is - 1)
      end do 
c
c           Determine sub-channel wavelengths
c
      do ic = 1,nc
         do is = 1,ns
            wvl(ic,is) = wvl_c(ic) + xs(is)*delwvl
         end do
      end do           
c
c           Interpolate optical coefficients at sub-channel wavelengths
c
      do ic = 1,nc
         do is = 1,ns
            l = 1             !use binary search for location of 
                              !  wvl(ic,is)
            u = nref
            m = (l + u)/2
            do while (m .gt. l)
               if (wvl(ic,is) .lt. wvl_r(m)) then
                     u = m
                  else
                     l = m
               end if
               m = (l + u)/2
            end do
            c = (wvl_r(l+1) - wvl(ic,is))/(wvl_r(l+1) - wvl_r(l))
            c1 = 1.0 - c
            a0(ic,is) = c*a0_r(l) + c1*a0_r(l+1)
            a1(ic,is) = c*a1_r(l) + c1*a1_r(l+1)
            a2(ic,is) = c*a2_r(l) + c1*a2_r(l+1)
                             !abs coeff = a0 + a1*(t-ts) + a2*(t-ts)**2
            b(ic,is) = c*b_r(l) + c1*b_r(l+1)
            rho = c*rho_r(l) + c1*rho_r(l+1)
            p0(ic,is) = 0.75*(1.0 + rho)/(1.0 + rho/2.0)
            p2(ic,is) = 0.75*(1.0 - rho)/(1.0 + rho/2.0)
                             !phase function = p0 + p2*cos(sc angle)**2
            flux(ic,is) = c*flux_r(l) + c1*flux_r(l+1)
         end do
      end do
c            
c           Determine weights
c           
      do ic = 1,nc
         sum = 0.0
         do is = 1,ns
            w(ic,is) = (1.0 - abs(xs(is)))*flux(ic,is)
            sum = sum + w(ic,is)
         end do
         do is = 1,ns
            w(ic,is) = w(ic,is)/sum
         end do
      end do
c
      return
      end