c file:                  solirr.f
c
c external routines:     zensun,solirr
c
c internal routines:     none
c
c internal common:       sundat

c=======================================================================


      function solirr(w,nf) 1
c
c input:
c     w          wavelength in microns
c
c     nf         solar spectrum database switch
c               -1 = read solar.dat, replace data in sundat common block
c                1 = 5s 
c                2 = lowtran7 
c                3 = modtran3 (degraded to 20 cm-1 resolution)
c
c               note: the 5s model uses a 1/w**4 power law for w > 4 um
c output:
c                estraterstrial solar flux in w/m2/micron
c
      implicit none

      real ws1,ws2,dws,wl1,wl2,wl3,dwl1,dwl2,wm1,wm2,dwm,dw,wmin,wmax,
     &     wave,w,solirr,f,sun1,sun2a,sun2b,sun3
      integer i,nw,ns,nla,nlb,nm,nf

      parameter (ws1=.25, ws2=4., dws=.005)
      parameter (wl1=0., wl2=28400., wl3=57490., dwl1=20., dwl2=10)
      parameter (wm1=100., wm2=49960, dwm=20.)

      parameter (ns=751, nla=1440, nlb=2910, nm=2494)

      common /sundat/ sun1(ns),sun2a(nla),sun2b(nlb),sun3(nm)
      data nw,wmin,wmax,dw/0,0.,0.,0./
      save nw,wmin,wmax,dw


      if(nf.eq.-1) then                          ! user defined 

        if(nw.eq.0) then                         ! read solar.dat
          open(13,file='solar.dat',status='old')
          read(13,*) nw,wmin,wmax
          if(nw.gt.ns+nla+nlb+nm) stop 'Too many points in solar.dat'
          read(13,*) (sun1(i),i=1,nw)
          dw=(wmax-wmin)/nw
        endif

        wave=1.e4/w
        i=(wave-wmin)/dw+1.00001
        solirr=0.
        if(i.lt.1 .or. i.gt.nw-1) return
        f=(wave-wmin-dw*(i-1))/dw
        solirr=(1.-f)*sun1(i)+f*sun1(i+1)

      elseif(nf.eq.1) then                       ! 5s model

        i=(w-ws1)/dws+1.00001
        f=(w-ws1-dws*(i-1))/dws
        
        if(i.le.ns-1) then 
          solirr=(1.-f)*sun1(i)+f*sun1(i+1)
        else
          solirr=sun1(ns)*(ws2/w)**4
        endif

      elseif(nf.eq.2) then                       ! lowtran7 model

        wave=1.e4/w
        if(wave.lt.wl2) then
          i=(wave-wl1)/dwl1+1.00001
          f=(wave-wl1-dwl1*(i-1))/dwl1
          solirr=(1.-f)*sun2a(i)+f*sun2a(i+1)
        elseif(wave.lt.wl3) then
          i=(wave-wl2)/dwl2+1.00001
          f=(wave-wl2-dwl2*(i-1))/dwl2
          solirr=(1.-f)*sun2b(i)+f*sun2b(i+1)
        endif

      elseif(nf.eq.3) then                       ! modtran3 model

        wave=1.e4/w
        i=(wave-wm1)/dwm+1.00001
        solirr=0.
        if(i.lt.1 .or. i.gt.nm-1) return
        f=(wave-wm1-dwm*(i-1))/dwm
        solirr=(1.-f)*sun3(i)+f*sun3(i+1)

      endif
      return
      end