c
c
c  diagnostic main program:
c
c        common / sr / alb(751),wmin,wmax,dw
c        dimension sc(4)
c        write(*,1000) (.25+.005*(i-1),i=1,751)
c        call suralb(1,0,sc)
c        write(*,1000) (alb(i),i=1,751)
c        call suralb(2,0,sc)
c        write(*,1000) (alb(i),i=1,751)
c        call suralb(3,0,sc)
c        write(*,1000) (alb(i),i=1,751)
c        call suralb(4,0,sc)
c        write(*,1000) (alb(i),i=1,751)
c        call suralb(5,0,sc)
c        write(*,1000) (alb(i),i=1,751)
c        call suralb(6,0,sc)
c        write(*,1000) (alb(i),i=1,751)
c   1000 format(10f10.4)
c        end

c file:                  filter.f
c
c external routines:     setfilt,filter
c
c internal routines      meteo, goese,goesw,avhr81,avhr82,avhr91,avhr92,
c                        avhr101,avhr102,avhr111,avhr112,gtr1,gtr2,nm410
c                        nm936,usersat,mfrsr1,mfrsr2,mfrsr3,mfrsr4,mfrsr5,
c                        mfrsr6,avhr83,avhr84,avhr85,setlow
c
c internal common        fltblk
c=======================================================================

      subroutine setfilt(isat,wlinf,wlsup,wlinc,wl1,wl2,nwl) 1,26
      parameter (iwvmx=1000)
      common /fltblk/filfun(iwvmx),wmin,wmax,dw,nnf
      
c
c input   isat      filter selection flag
c 
c         wlinf     lower wavelength limit
c                      or central wavelength for isat=-2,-3
c         wlsup     upper wavelength limit
c                      or equivalent width for isat=-2,-3
c         wlinc     wavelength increment
c output
c         wl1       lower wavelength limit
c         wl2       upper wavelength limit
c         nwl       number of sample points
c
c NOTE: actual filter values are returned by the function "filter"
c
      if(wlsup.eq.0.and.isat.le.2) then
        write(*,*) 'Error -- must specify non-zero WLSUP when ISAT=',isat
        stop
      endif

      if(isat.eq.-4) then   ! center, equivalent width, gaussian
        sqpi=sqrt(acos(-1.))
        xc=2
        xlim=xc*sqpi
        wmin=wlinf-xc*wlsup
        wmax=wlinf+xc*wlsup
        nnf=1000
        dw=(wmax-wmin)/(nnf-1)
        do i=1,nnf
          xx=-xlim+(i-1)*(2*xlim)/(nnf-1)
          filfun(i)=exp(-xx**2)
        enddo
      endif
       
      if(isat.eq.-3) then  ! center, equivalent width, triangular
        wmin=wlinf-wlsup
        wmax=wlinf+wlsup
        filfun(1)=0.
        filfun(2)=1.
        filfun(3)=0.
        nnf=3
        dw=wlsup
      endif

      if(isat.eq.-2) then   ! center, width, flat top
        wmin=wlinf-.5*wlsup
        wmax=wlinf+.5*wlsup
        nnf=0
      endif

      if(isat.eq.0.) then ! wl_min, wl_max, flat  top
        wmin=wlinf
        wmax=wlsup
        nwl=1
        nnf=0
        if(wlinf.eq.wlsup) then
          wlinc=1.
          wl1=wmin
          wl2=wmax
          return
        endif
      endif

      if(isat.eq.-1) call usersat(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.1)  call meteo(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.2)  call goese(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.3)  call goesw(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.4)  call avhr81(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.5)  call avhr82(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.6)  call avhr91(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.7)  call avhr92(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.8)  call avhr101(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.9)  call avhr102(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.10) call avhr111(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.11) call avhr112(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.12) call gtr1(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.13) call gtr2(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.14) call nm410(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.15) call nm936(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.16) call mfrsr1(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.17) call mfrsr2(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.18) call mfrsr3(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.19) call mfrsr4(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.20) call mfrsr5(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.21) call mfrsr6(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.22) call avhr83(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.23) call avhr84(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.24) call avhr85(filfun,wmin,wmax,dw,nnf)
      if(isat.eq.25) call setlow(filfun,wmin,wmax,dw,nnf)

      if(wmin.lt.0.2) then
        write(*,*) 'Error in setfilt -- illegal wavelength limits '
        write(*,*) wmin,wmax
        stop
      endif

      if(wlinc.gt.1.) then
c       equal spacing in wavenumber
        nwl=((10000./wmin) - (10000./wmax))/wlinc + 1.
      elseif(wlinc.lt.0.) then
c       equal spacing in ln(wavelength) and ln(wavenumber)
        nwl=1+log(wmax/wmin)/abs(wlinc)
      else
c       use default wavelength increment
        if(wlinc.eq.0)
     &       wlinc=(wmax-wmin)/max(10,1+int((wmax-wmin)/0.005))
c       equal spacing in wavelength
        nwl=nint((wmax-wmin)/wlinc)+1
      endif

      wl1=wmin
      wl2=wmax

      return
      end