pro plot_model_inputs, date, datadir, imagedir
;  Code to read in and plot model forcing data from the DW ascii inpu
;  files as opposed to the model output as in plot_global_force.pro
;
;  This version works with KC model input

varstrin1='sst'
varstrin2='flx'

varstr_qh='qh'
varstr_qe='qe'
varstr_lw='lw'
varstr_swh='swh'
varstr_per='per'

;  Specify directory

flxfile=datadir+'/'+date+'_'+varstrin2+'.dat'
sstfile=datadir+'/'+date+'_'+varstrin1+'.dat'
plotdir=imagedir+'/'

string = 'wc -l ' + sstfile + ' | awk ''{printf "%s",$1}'''
spawn, string, lines
ngrid=long(lines[0])
print,ngrid

openr,unit,flxfile,/get_lun
openr,unit2,sstfile,/get_lun

;  Specify simulation/data file characteristics

tempres = 1     ; Number of output lines per hour
nhours = 48     ; Number of hours in simulation
ninputs = 9     ; Number of input flux files

minlat=-60.0
minlon=0.0
maxlat=60.0
maxlon=360.0

res = 0.5

nxgrid = (maxlon-minlon)/res
nygrid = (maxlat-minlat)/res + 1

gridarr=fltarr(6)
timearr=fltarr(13,ninputs)    ; for 15 min res w/o first point

qh=fltarr(nxgrid,nygrid,4)
qh(*,*,*) = -999.99

qe=fltarr(nxgrid,nygrid,4)
qe(*,*,*) = -999.99

lw=fltarr(nxgrid,nygrid,4)
lw(*,*,*) = -999.99

swh=fltarr(nxgrid,nygrid,4)
swh(*,*,*) = -999.99

per=fltarr(nxgrid,nygrid,4)
per(*,*,*) = -999.99

;  Loop through cells in the grid

for igrid = 0L, ngrid-1L do begin

   ;  Read in grid cell info

   readf,unit2,gridarr

   xindex = fix(gridarr[4]) - 1
   yindex = fix(gridarr[5]) - 1

   ;  Read in temporal information for grid cell

   readf,unit,timearr

   ;  Start loop through hours to be plotted

   for hr=0,3 do begin

      ;  Specify hour to be plotted

      outhr = 0 + hr*6
      ;indexout = fix(outhr*tempres) - 1    ; -1 is for idl index start at 0
      indexout = hr+4    ; -1 is for idl index start at 0

      lw(xindex,yindex,hr) = timearr(4,indexout)
      qe(xindex,yindex,hr) = timearr(5,indexout)
      qh(xindex,yindex,hr) = timearr(6,indexout)
      swh(xindex,yindex,hr) = timearr(7,indexout)
      per(xindex,yindex,hr) = timearr(8,indexout)

   endfor

endfor     ; end loop through grid

;  Loop through hours again for plotting

for hr=0,3 do begin

   outhr = 0 + hr*6

;  Specify hour to be plotted

   outfile_qe=string(plotdir,varstr_qe,outhr,format='(a,a,"_",i2.2,".png")')
   outfile_qh=string(plotdir,varstr_qh,outhr,format='(a,a,"_",i2.2,".png")')
   outfile_lw=string(plotdir,varstr_lw,outhr,format='(a,a,"_",i2.2,".png")')
   outfile_swh=string(plotdir,varstr_swh,outhr,format='(a,a,"_",i2.2,".png")')
   outfile_per=string(plotdir,varstr_per,outhr,format='(a,a,"_",i2.2,".png")')

   ;  Now plot

   outstr=string(date,outhr,format='("Latent flux, ",a," ",i2.2, "Z")')
   pimage, qe(*,*,hr), $
   /order, $
   minlat=minlat,maxlat=maxlat,minlon=minlon,maxlon=maxlon, $
   nxtick=6,nytick=4, /label, $
   ;mis_val=-999.99,mis_col=1,coast=2,/vertical, $
   ;mis_val=-999.99,mis_col=3,coast=2,ccolor=1,/vertical, $
   mis_val=-999.99,mis_col=3,coast=2,ccolor=1, $
   vmin=-100,vmax=500.0, nctick=6, $
   ;/pngout, pngfile=outfile_q,xws=800,yws=350, $ ; with vertical
   /pngout, pngfile=outfile_qe,xws=700,yws=400, $
   ;/eps_port, ps_file=outfile, $
   title=outstr, $
   units='Latent Flux (W/m2)'

   outstr=string(date,outhr,format='("Sensible flux, ",a," ",i2.2, "Z")')
   pimage, qh(*,*,hr), $
   /order, $
   minlat=minlat,maxlat=maxlat,minlon=minlon,maxlon=maxlon, $
   nxtick=6,nytick=4, /label, $
   ;mis_val=-999.99,mis_col=1,coast=2,/vertical, $
   ;mis_val=-999.99,mis_col=3,coast=2,ccolor=1,/vertical, $
   mis_val=-999.99,mis_col=3,coast=2,ccolor=1, $
   vmin=-20,vmax=100.0, nctick=6, $
   ;/pngout, pngfile=outfile_q,xws=800,yws=350, $ ; with vertical
   /pngout, pngfile=outfile_qh,xws=700,yws=400, $
   ;/eps_port, ps_file=outfile, $
   title=outstr, $
   units='Sensible Flux (W/m2)'

   outstr=string(date,outhr,format='("LW Flux, ",a," ",i2.2, "Z")')
   pimage, lw(*,*,hr), $
   /order, $
   minlat=minlat,maxlat=maxlat,minlon=minlon,maxlon=maxlon, $
   nxtick=6,nytick=4, /label, $
   ;mis_val=-999.99,mis_col=1,coast=2,/vertical, $
   ;mis_val=-999.99,mis_col=3,coast=2,ccolor=1,/vertical, $
   mis_val=-999.99,mis_col=3,coast=2,ccolor=1, $
   vmin=0,vmax=200.0, nctick=5, $
   ;/pngout, pngfile=outfile_u,xws=800,yws=350, $
   /pngout, pngfile=outfile_lw,xws=700,yws=400, $
   ;/eps_port, ps_file=outfile2, $
   title=outstr, $
   units='LW Flux (m/s)'

   outstr=string(date,outhr,format='("SWH, ",a," ",i2.2, "Z")')
   pimage, swh(*,*,hr), $
   /order, $
   minlat=minlat,maxlat=maxlat,minlon=minlon,maxlon=maxlon, $
   nxtick=6,nytick=4, /label, $
   ;mis_val=-999.99,mis_col=1,coast=2,/vertical, $
   ;mis_val=-999.99,mis_col=3,coast=2,ccolor=1,/vertical, $
   mis_val=-999.99,mis_col=3,coast=2,ccolor=1, $
   vmin=0,vmax=10.0, nctick=5, $
   ;/pngout, pngfile=outfile_u,xws=800,yws=350, $
   /pngout, pngfile=outfile_swh,xws=700,yws=400, $
   ;/eps_port, ps_file=outfile2, $
   title=outstr, $
   units='SWH (m)'

   outstr=string(date,outhr,format='("Wave Period, ",a," ",i2.2, "Z")')
   pimage, per(*,*,hr), $
   /order, $
   minlat=minlat,maxlat=maxlat,minlon=minlon,maxlon=maxlon, $
   nxtick=6,nytick=4, /label, $
   ;mis_val=-999.99,mis_col=1,coast=2,/vertical, $
   ;mis_val=-999.99,mis_col=3,coast=2,ccolor=1,/vertical, $
   mis_val=-999.99,mis_col=3,coast=2,ccolor=1, $
   vmin=0,vmax=30.0, nctick=6, $
   ;/pngout, pngfile=outfile_u,xws=800,yws=350, $
   /pngout, pngfile=outfile_per,xws=700,yws=400, $
   ;/eps_port, ps_file=outfile2, $
   title=outstr, $
   units='Period (s)'

endfor

end