pro task86,iout

; itask=86 Graph 1d / 3d as a function of vertical optical depth and cloud distance
;          Graph 1d / 3d as a function of solar zenith angle and cloud distance
;          Graph 1d / 3d as a function of cloud height and cloud distance
;          Graph 1d / 3d as a function of surface reflectance and cloud distance

; ******************************
; Output ps file
  fileps5='task86_4panels.ps'

; Internal savesets used to calculate filesp5
  fsaves=strarr(4)
  fsaves(0)='task86tau.sav'
  fsaves(1)='task86refl.sav'
  fsaves(2)='task86sza.sav'
  fsaves(3)='task86cldheight.sav'

; Write out to ascii files
  fwrites=strarr(4)
  fwrites(0)='task86tau.dat'
  fwrites(1)='task86refl.dat'
  fwrites(2)='task86sza.dat'
  fwrites(3)='task86cldheight.dat'

; ******************************
; Graph all 4 panels on the same page

; Open the ps file
   set_plot, 'ps'
   device, /color, bits_per_pixel=8, file=fileps5,$
    /portrait,/inch,ysize=9.0,xsize=7.0,xoffset=0.6,yoffset=1.0

   titleg=' '

; One graph per page
   !P.MULTI=[0,2,2]

; Load the color table
   loadct,26

; Bold fonts
   !P.font=0

; Loop over the 4 panels
   for ngr=0,3 do begin

; **
   if (ngr eq 0) then begin
;  save,npts,xvec,yvec,veccolor,isymn,clrval,$
;  nlines,nvec,xvecn,yvecn,$
;  labelx,labely,title,$
;  xmin,xmax,ymin,ymax,$
;  x1a,y1a,x2a,y2a,$
;  ioplot,igraphpts,sizept,$
;  cs,cs2,cthick,$
;  ixlog,iylog,$
;  igrwords,ngrwords,grwords,xwords,ywords,wordc,$
;  file=fsaves(1)

; optical depth
    restore,fsaves(0)
    fwrite=fwrites(0)
    x1a=0.15
    x2a=0.45
    y1a=0.60
    y2a=0.90
    title='Gas Vert Optical Depth'
    xwords=fltarr(nlines)
    ywords=fltarr(nlines)
     xplace=0.5*(x1a+x2a)
     yw=0.95*y2a
     yeps=0.03
    for j=0,nlines-1 do begin
     xwords(j)=xplace
     yw=yw-yeps
     ywords(j)=yw
    endfor
   endif

   if (ngr eq 1) then begin
; sza
    restore,fsaves(2)
    fwrite=fwrites(2)
    x1a=0.60
    x2a=0.90
    y1a=0.60
    y2a=0.90
    title='Solar Zenith Angle'
     xplace=0.5*(x1a+x2a)
     yw=0.95*y2a
     yeps=0.03
    for j=0,nlines-1 do begin
     xwords(j)=xplace
     yw=yw-yeps
     ywords(j)=yw
    endfor
   endif

   if (ngr eq 2) then begin
; cloud height
    restore,fsaves(3)
    x1a=0.15
    x2a=0.45
    y1a=0.20
    y2a=0.50
    title='Cloud Height'
     xplace=0.5*(x1a+x2a)
     yw=0.8*y2a
     yeps=0.03
    for j=0,nlines-1 do begin
     xwords(j)=xplace
     yw=yw-yeps
     ywords(j)=yw
    endfor
   endif

   if (ngr eq 3) then begin
; surface reflectance
    restore,fsaves(1)
    fwrite=fwrites(1)
    x1a=0.60
    x2a=0.90
    y1a=0.20
    y2a=0.50
    title='Surface Reflectance'
     xplace=0.5*(x1a+x2a)
     yw=0.8*y2a
     yeps=0.03
    for j=0,nlines-1 do begin
     xwords(j)=xplace
     yw=yw-yeps
     ywords(j)=yw
    endfor
   endif

; ***
; NVEC            INT       =       21
; XVECN           FLOAT     = Array[21, 7]
; YVECN           FLOAT     = Array[21, 7]
; LABELX          STRING    = 'Distance from Cloud (km)'
; LABELY          STRING    = '1D / 3D intensity ratio'
; TITLE           STRING    = 'Dependence on Vertical Optical Depth'
   iwr=1
  if (iwr eq 1) then begin
    idat=60
   openw,idat,fwrite

   printf,idat,'  ',title
   printf,idat,'  labelx ',labelx
   printf,idat,'  labely ',labely
   printf,idat,'  nlines ',nlines
   printf,idat,'  grwords ',grwords
   printf,idat,'  i, cloud distance, data for nlines '

   printf,iout,' '
   printf,iout,'  ',title
   printf,iout,'  labelx ',labelx
   printf,iout,'  labely ',labely
   printf,iout,'  nlines ',nlines
   printf,iout,'  grwords ',grwords
   printf,iout,'  i, cloud distance, data for nlines '

    ywr=fltarr(nlines)

   for i=0,nvec-1 do begin
    for j=0,nlines-1 do begin
     ywr(j)=yvecn(i,j)
    endfor
    printf,idat,format='(2x,i2,8(1x,f12.6))',$
     i,xvecn(i,0),ywr
    printf,iout,format='(2x,i2,8(1x,f12.6))',$
     i,xvecn(i,0),ywr
   endfor
   close,idat
  endif

; ****
     nopr2=0
   graphlines100,nopr2,iout,$
   npts,xvec,yvec,veccolor,isymn,clrval,$
   nlines,nvec,xvecn,yvecn,$
   labelx,labely,title,$
   xmin,xmax,ymin,ymax,$
   x1a,y1a,x2a,y2a,$
   ioplot,igraphpts,sizept,$
   cs,cs2,cthick,$
   ixlog,iylog,$
   igrwords,ngrwords,grwords,xwords,ywords,wordc

   endfor
; llop over the ngr panels

; **
   device,/close
   print,fileps5

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