pro task20,iout

; itask=20 Look at MODTRAN6 PLT files to estimate sensitivity to
;          the number of streams in the calculation. 3 files

  printf,iout,' itask=20 Look at MODTRAN6 PLT files to estimate sensitivity to'
  printf,iout,'          the number of streams in the calculation. 3 files'

; ******************************
; Specify the ascii file listing
   flist='/ur/massie/libera/json/modtran6_streamscalc/z.list'
    noprr=1
   readflist,noprr,iout,flist,nfiles,files

   nstreams=fltarr(nfiles)
   nstreams=[4,8,16,32]

; ******************************
; The files are for the 400 -1700 nm range
; Will sum values (mW ster-1 cm-2 nm-1) over approx VIIRS bandpasses

; For the 400 - 1700 nm range

; See Documents Libera VIIRS bands and bandwidths.pdf
     nviirs=6
    vnm=fltarr(nviirs)
    vbandpass=fltarr(nviirs)
    vnm=[0.488, 0.55, 0.64, 0.86, 1.24, 1.61]
    vbandpass=[0.02, 0.02, 0.08, 0.039, 0.020, 0.06]  

; 3.74 has a bandpass of 0.18 microns

; ******************************
; Specify the nm range of the viirs values
     vnm1=fltarr(nviirs)
     vnm2=fltarr(nviirs)
     printf,iout,'  '
     printf,iout,'  task20: nviirs ',nviirs
     printf,iout,'  task20: viirs filter range in nm'
     printf,iout,'  task20: i,vnm1(i),vnm2(i)'
    for i=0,nviirs-1 do begin
      diff=0.5*vbandpass(i)*1000.0
      a1=1000.0*vnm(i)
     vnm1(i)=a1-diff
     vnm2(i)=a1+diff
     printf,iout,'  ',i,vnm1(i),vnm2(i)
    endfor

; ******************************
; Read in the data

    jfiles=nfiles
    ndat=1145

    wavenm=fltarr(ndat,jfiles)
    dwavenm=fltarr(ndat,jfiles)
    spectra=fltarr(ndat,jfiles)

; Loop over the 4 files
   for j=0,jfiles-1 do begin

     idat=50
     openr,idat,files(j)

     i1=0
     i2=ndat-1
     iskip=-1
    for i=i2,i1,iskip do begin
     readf,idat,a1,a2
     wavenm(i,j)=a1
; spectra values are in mW ster-1 cm-2 nm-1
     spectra(i,j)=a2
    endfor

    close,idat

    print,'  read in from ',files(j)

   endfor

; ******************************
; Calculate dwavenm
   for j=0,nfiles-1 do begin
     i1=0
     i2=ndat-1
     ilast=i2
    for i=i1,i2-1 do begin
      i3=i+1
     diff=wavenm(i3,j)-wavenm(i,j)
     dwavenm(i,j)=diff
    endfor
     dwavenm(ilast,j)=dwavenm(ilast-1,j)
   endfor

; ******************************
; Sum of the spectra points to obtain radiance in mW ster-1 cm-2
     viirsrad=fltarr(nviirs,jfiles)
     viirsradn=fltarr(nviirs,jfiles)

    for j=0,jfiles-1 do begin
    for i=i1,i2-1 do begin

; For this wavelength, determine if you are in range of a viirs filter
      kval=-999
     for k=0,nviirs-1 do begin
      if ((wavenm(i,j) ge vnm1(k)) and (wavenm(i,j) le vnm2(k))) then begin 
       kval=k
       goto,jump10
      endif
     endfor
     jump10:mn=0

     if (kval ge 0) then begin
      a2=dwavenm(i,j)*spectra(i,j)
      viirsrad(kval,j)=viirsrad(kval,j)+a2
      viirsradn(kval,j)=viirsradn(kval,j)+1.00
     endif

    endfor
    endfor

    printf,iout,'  '
    printf,iout,'  task20: min and max wavenm ',min(wavenm),max(wavenm)
    printf,iout,'  task20: min and max dwavenm ',min(dwavenm),max(dwavenm)
    printf,iout,'  task20: min and max spectra ',min(spectra),max(spectra)
    printf,iout,'  '
    printf,iout,'  task20: min and max viirsrad ',min(viirsrad),max(viirsrad)
    printf,iout,'  task20: min and max viirsradn ',min(viirsradn),max(viirsradn)
 
; ******************************
; Write out values
      printf,iout,'  '
      printf,iout,'  task20: viirsrad(i,j) is sum of product dwavenm(i,j) * spectra(i,j) '
      printf,iout,'  task20: viirsrad(i,j) for i=0,nfiles-1'
      printf,iout,'  task20: viirsradn(i,j) is number of spectral points added for each filter'
      printf,iout,'  task20: perc(i)=100.0*(a0-a32)/a32 where a0=viirsrad(i,j), a32=viirsrad(i,3)'
      printf,iout,'  task20: viirsradn(i,j) for i=0,nfiles-1'

    for j=0,jfiles-1 do begin

      printf,iout,'  '
      printf,iout,'  task20: j ',j
      printf,iout,'  task20: nstreams(j) ',nstreams(j)

       vec=fltarr(nviirs)
       vecn=fltarr(nviirs)
       perc=fltarr(nviirs)
     for i=0,nviirs-1 do begin
      vec(i)=viirsrad(i,j)
      a32=viirsrad(i,3)
      a0=viirsrad(i,j)
      perc(i)=100.0*(a0-a32)/a32
      vecn(i)=viirsradn(i,j)
     endfor

     printf,iout,'  '
     printf,iout,format='(2x,6(1x,f12.6))',vec
     printf,iout,format='(2x,6(1x,f10.4))',perc
     printf,iout,format='(2x,6(1x,i4))',vecn

     endfor

; ******************************
   return
   end