;
pro pltline,info,animate=animate,ps=ps,png=png
;
; Plot line type plots, profile(zp/field), meridional(field/latitude), and zonal(field/longitude)
;
;  Grab data arrays needed
;
file = *info.pfile
line = *info.pline
field = line.field
fields = *info.fields

IF line.plottype EQ 'time' THEN BEGIN

  fields1D = *info.fields1D
  field1D = line.field1D
  times = *field1D.times
  data = *line.data
  rhyaxis = 0
  ntimes = n_elements(times)
  
ENDIF ELSE BEGIN  

flats = *fields[line.ifield].lats
nflats = n_elements(flats)
data = *line.data
lats = *line.lats & nlat = n_elements(lats)
lons = *line.lons & nlon = n_elements(lons)
levs = *line.levs 
nlev = n_elements(levs)
rhyaxis = line.rhyaxis
;
;  Turn off right hand axis for difference file or magnetic grid
;
if rhyaxis gt 0 and (file.diffs gt 0 or field.grid_type eq 'magnetic') then $
  rhyaxis = 0
;
;  Check if latitude range is set
;  
latrangeset = 0
if line.latmin gt flats[0] or line.latmax lt flats[nflats-1] then $
  latrangeset = 1

ENDELSE ;time or not time plot
;
;  Check for NANs in data
;
find_nans,info,field.name,data,numnans
if numnans gt 0 then print,'>>NaNs found:',numnans
if numnans gt 0 then  begin
  rhyaxis = 0   
  line.rhyaxis = 0
endif

;print,'pltline: nflats = ',nflats,' field lats = ' & print,flats
;print,'pltline: latmin,max = ',line.latmin,line.latmax,' nlat=',nlat,$
;  ' lats = ' & print,lats

; Write binary file for debug:
;
;datfile='lineslice.dat'
;openw,lu,datfile,/get_lun
;help,nlat & help,lats
;writeu,lu,nlat
;writeu,lu,float(lats)
;writeu,lu,nlev
;writeu,lu,float(levs)
;writeu,lu,data
;close,lu
;free_lun,lu
;print,'Wrote binary file ',datfile
;print,'nlat=',nlat,' lats=' & print,lats
;print,'nlev=',nlev,' levs=' & print,levs
;print,'data min,max=',min(data),max(data)
;
; Write ascii data file for debug:
;
;datfile='lineslice.txt'
;openw,lu,datfile,/get_lun
;help,nlat & help,lats
;printf,lu,format="('nlat=',i4)",nlat
;printf,lu,format="('lats=',/,(6e12.4))",float(lats)
;printf,lu,format="('nlev=',i4)",nlev
;printf,lu,format="('levs=',/,(6e12.4))",float(levs)
;printf,lu,format="('data=',/,(6e12.4))",data
;close,lu
;free_lun,lu
;print,'Wrote ascii file ',datfile
;print,'nlat=',nlat,' lats=' & print,lats
;print,'nlev=',nlev,' levs=' & print,levs
;print,'data min,max=',min(data),max(data)
;
if not keyword_set(animate) and not keyword_set(ps) and $
   not keyword_set(png) then begin
  widget_control,line.draw,get_value=index
  wset,index
  erase
endif
;
; Check for all missing value:
;
nmissing = 1 ; do not init this to 0 because check_missing uses "if keyword_set"
missing_value = float(line.missing_value)
text = 1 & if keyword_set(animate) then text = 0
nodata = check_missing(data,missing_value,field.name,text=text,$
  nmissing=nmissing)
;if nodata eq 1 then return
if nmissing gt 0 then print,'Line slices: Field ',field.name,': Found ',$
  nmissing,' missing data points (out of ',n_elements(data),')'
;
; Field min,max:
; If line.fmin >  line.fmax, then leave as such for the animator but use 
;                              field extremes for plotting current frame.
; If line.fmin =  line.fmax, then reset to data extremes (for the slice).
; If line.fmin <  line.fmax, then leave as such as for plotting.
;
;fminmax,data,framemin,framemax,missing_value
if line.fmin ge line.fmax then begin  ; line.fmin >= line.fmax
  fminmax,data,fmin,fmax,missing_value
  if line.fmin eq line.fmax AND line.fmax GT 0 then begin
    line.fmin=fmin-fmin*0.01 & line.fmax=fmin+fmin*0.01
    print,'pltline reset line.fmin,max=',line.fmin,line.fmax
  endif
; endif else $                          ; line.fmin > line.fmax: leave for anim
;   print,'pltline did NOT reset min,max: line.fmin,max=',line.fmin,line.fmax
endif else begin                        ; lon.fmin < lon.fmax
  fmin = line.fmin & fmax = line.fmax
endelse
;
; Image position:
;
;xur = .90 & yur = .90
;xll = .15 & yll = .31
xur = 0.908 & yur = 0.955
xll = .0 & yll = .11
;if rhyaxis gt 0 then begin 
;  xur=xur-.025 & xll=xll-.025 
;endif
axespos  = [xll,yll,xur,yur]
xmid = .5*(xur+xll)
;
;  Plot top field label:
;
IF line.plottype EQ 'time' THEN BEGIN

  if strlen(field1D.long_name) gt 0 then fieldlab = field1D.long_name $
    else fieldlab = field1D.name
  if file.diffs le 0 then begin
    if strlen(field1D.units) gt 0 then $
     fieldlab = fieldlab + ' (' + strcompress(field1D.units,/remove_all) + ')'
  endif else fieldlab = fieldlab + ' (' + field1D.difftype + ' ' + $
   strcompress(field1D.units,/remove_all) + ')'
 

ENDIF ELSE BEGIN  

  if (line.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
     (line.log10 eq 'all fields') then begin
    if strlen(field.long_name) gt 0 then fieldlab = 'LOG10 '+field.long_name $
      else fieldlab = 'LOG10 '+field.name
  endif else begin
    if strlen(field.long_name) gt 0 then fieldlab = field.long_name $
      else fieldlab = field.name
  endelse
  if file.diffs le 0 then begin
    if strlen(field.units) gt 0 then $
     fieldlab = fieldlab + ' (' + strcompress(field.units,/remove_all) + ')'
  endif else fieldlab = fieldlab + ' (' + field.difftype + ' ' + $
   strcompress(field.units,/remove_all) + ')'

ENDELSE
;
; Top label :
;
uthrs = float(line.mtime[1])+float(line.mtime[2])/60.
toplab = string(format="(' DAY=',i3,' UT=',f5.2)",line.mtime[0],uthrs)
if line.zonalmean gt 0 then begin
  toplab = toplab + string(format="(' LAT=',f7.2)",line.slat) + ' ZONAL MEAN'
endif else if line.meridionalmean gt 0 then begin
  toplab = toplab + string(format="(' LST=',f5.2)",line.sslt)
  toplab = toplab + ' MERIDIONAL MEAN'
endif else if field.grid_type eq 'geographic' then begin
  toplab = toplab + string(format="(' LST=',f5.2)",line.sslt)
  toplab = toplab + string(format="(' LAT=',f7.2)",line.slat)
  toplab = toplab + string(format="(' LON=',f7.2)",line.slon)
endif else begin
  toplab = toplab + string(format="(' LST=',f5.2)",line.sslt)
  toplab = toplab + string(format="(' MLAT=',f7.2)",line.slat)
  toplab = toplab + string(format="(' MLON=',f7.2)",line.slon)
endelse

if line.amplitude GT 0 then begin 
  toplab = string(format="(' DAY=',i3,' UT=',f5.2)",line.mtime[0],uthrs) + $
           string(format="(' LAT=',f7.2)",line.slat) + $
	   string(format="(' WAVENUMBER=',i2)",line.wavenumber) + ' AMPLITUDE'
endif
;
; Bottom label:
;
 botlab1 = strcompress(field.ncfile,/remove_all)
;
; Contour levels were made by mkclevels.pro.
;
;print,'pltlon: nlevels=',nlevels,' clevels=' & print,clevels

chsize = 1.2     ; info label character size
;
; Contour will not draw y-axis if extra non-linear right-hand axis is 
; to be drawn by pro altyax:
;
yaxstyle = 1 & if rhyaxis gt 0 then yaxstyle = 9
;
; Display image for each line plot type:
;
case line.plottype of
;
; Display image only. Use contour,/nodata to draw axes, ticmarks, labels:
;
  'profile': begin
     xtitle = info.ftype + ' ' + fieldlab ; field
     ytitle = field.levname + ' (' + field.levunits + ')' ; zp or height
;
;  WACCM vertical coordinate converted to log pressure so change label
;
     if info.ftype eq 'WACCM' THEN ytitle = 'ZP log pressure level (ln(p0/p))' 
;
;  Set y-axis title if geopotential or geometric altitude
;
     IF line.plotz EQ 1 THEN ytitle = ' GEOPOTENTIAL ALTITUDE (km) '
     IF line.plotz EQ 2 THEN ytitle = ' GEOMETRIC ALTITUDE (km) '
;
;  Only plot good data (not missing data)
;
     indexGood = WHERE(data NE line.missing_value)
     dataP = data(indexGood)
     levsP = levs(indexGood)
;
;  For FFT output need to set labels and frequency axis
;     
     IF line.amplitude GT 0 THEN BEGIN

       IF line.amplitude EQ 2 THEN BEGIN
         xtitle = xtitle + ' WAVENUMBER ' + 'AMPLITUDE LOG10'
       ENDIF ELSE IF line.amplitude EQ 1 THEN BEGIN
         xtitle = xtitle + ' WAVENUMBER ' + 'AMPLITUDE'
       ENDIF
       
     ENDIF

     plot,dataP, levsP, $
         xrange=[line.fmin,line.fmax], xstyle=1, $ ; exact axis
         yrange=[line.levmin,line.levmax], ystyle=yaxstyle, $ ; exact axis
         xlog=0, ylog=0, $
         charsize=1.1, $
         title=toplab, $
         xtitle = xtitle, $
         ytitle = ytitle, $
         subtitle = botlab1, $
         psym = -2, $
         xmargin = [10,7], $
         ymargin = [5,2]

     if rhyaxis gt 0 and line.plotz GE 1 then presyax,info,line,line.plottype,xll,xur,yll,yur,chsize

     if rhyaxis gt 0 and numnans eq 0 and line.plotz eq 0 then begin
       altyax,info,line,line.plottype,xll,xur,yll,yur,chsize 
     endif

  end

  'meridional': begin
;
; Get heights from Z (or Z3 from WACCM) field:
;
    fields = *info.fields
    getzprof,info,line,'line',fields,fzheights
    
    vfnd, *field.levs, fzheights, *line.levs, zheights, cut=1, linear=1
;
; Titles, labels, tic marks, etc:
;
; Top label (under field label):
;
     uthrs = float(line.mtime[1])+float(line.mtime[2])/60.
     toplab = string(format="('DAY=',i3,' UT=',f5.2)",line.mtime[0],uthrs)
     if line.zonalmean gt 0 then begin
       if line.plotz eq 0 then toplab = toplab + ' ZONAL MEAN' + $
                                 string(format="(' LEV=',f6.2)",line.slev) + $
                                  string(format="('(',f6.2)",zheights[line.ilev]) + "KM)"  
       if line.plotz eq 1 then toplab = toplab + ' ZONAL MEAN' + string(format="(' GPALT LEV=',f6.2)",line.alt) + "KM"
       if line.plotz eq 2 then toplab = toplab + ' ZONAL MEAN' + string(format="(' GMALT LEV=',f6.2)",line.alt) + "KM"
       if line.plotz eq 3 then toplab = toplab + ' ZONAL MEAN' + string(format="(' ZGALT LEV=',f6.2)",line.alt) + "KM"
     endif else if line.amplitude GT 0 THEN BEGIN 
       toplab = toplab + ' WAVENUMBER AMPLITUDE' + $
       string(format="(' LEV=',f6.2)",line.slev) + $
        string(format="('(',f6.2)",zheights[line.ilev]) + "KM)"
     endif else if field.grid_type eq 'geographic' and line.plotz eq 0 then begin
       toplab = toplab + string(format="(' LST=',f5.2)",line.sslt)
       toplab = toplab + string(format="(' LON=',f7.2)",line.slon)
       toplab = toplab + string(format="(' LEV=',f6.2)",line.slev) + string( format="('(',f6.2)",zheights[line.ilev]) + "KM)"
     endif else if field.grid_type eq 'geographic' and line.plotz gt 0 then begin
       toplab = toplab + string(format="(' LST=',f5.2)",line.sslt)
       toplab = toplab + string(format="(' LON=',f7.2)",line.slon)
       if line.plotz eq 1 then toplab = toplab + string(format="(' GPALT LEV=',f6.2)",line.alt) + "KM"
       if line.plotz eq 2 then toplab = toplab + string(format="(' GMALT LEV=',f6.2)",line.alt) + "KM"
       if line.plotz eq 3 then toplab = toplab + string(format="(' ZGALT LEV=',f6.2)",line.alt) + "KM"
     endif else begin
       toplab = toplab + string(format="(' LST=',f5.2)",line.sslt)
       toplab = toplab + string(format="('  MLON = ',f7.2)",line.slon)
       toplab = toplab + string(format="('  MLEV = ',f7.2)",line.slev)
     endelse
     if field.grid_type eq 'geographic' then begin
       xtitle = 'LATITUDE (DEGREES GEOGRAPHIC)'
     endif else begin
       xtitle = 'LATITUDE (DEGREES MAGNETIC)'
     endelse	 
     ytitle = info.ftype + ' ' + fieldlab ; field
;
;  For FFT output need to set labels and frequency axis
;     
     IF line.amplitude GT 0 THEN BEGIN

       IF line.amplitude EQ 2 THEN BEGIN
         ytitle = ytitle + ' WAVENUMBER ' + ' AMPLITUDE LOG10'
       ENDIF ELSE IF line.amplitude EQ 1 THEN BEGIN
         ytitle = ytitle + ' WAVENUMBER ' + ' AMPLITUDE'
       ENDIF
       
     ENDIF

     plot, lats, data, $
         xrange=[line.latmin,line.latmax], xstyle=1, $ ; exact axis
         yrange=[line.fmin,line.fmax], ystyle=1, $ ; exact axis
         xlog=0, ylog=0, $
         charsize=1.1, $
         title=toplab, $
         xtitle = xtitle, $
         ytitle = ytitle, $
         subtitle = botlab1, $
         psym = -2, $
         xmargin = [10,7], $
         ymargin = [5,2]
;         ticklen=0, $
;         background = 255, $
;         color = 0
;         xtickformat='empty_string', ytickformat='empty_string'
;         min_value=NULL_VALUE
;         color=color
;         noerase=noerase
;         /nodata
;         position=self.plot1d_position, $

  end

  'zonal': begin
;
; Get heights from Z (or Z3 from WACCM) field:
;
     fields = *info.fields
     getzprof,info,line,'line',fields,fzheights
    
     vfnd, *field.levs, fzheights, *line.levs, zheights, cut=1, linear=1
;
; Titles, labels, tic marks, etc:
;
; Top label (under field label):
;
     uthrs = float(line.mtime[1])+float(line.mtime[2])/60.
     toplab = string(format="('DAY=',i3,' UT=',f5.2)",line.mtime[0],uthrs)
     if line.meridionalmean gt 0 then begin
       toplab = toplab + string(format="(' LST=',f5.2)",line.sslt)
       toplab = toplab + ' MERIDIONAL MEAN' + $
         string(format="(' LEV=',f6.2)",line.slev) + $
         string( format="('(',f6.2)",zheights[line.ilev]) + "KM)"
     endif else if field.grid_type eq 'geographic' and line.plotz eq 0 then begin
       toplab = toplab + string(format="(' LAT=',f7.2)",line.slat)
       toplab = toplab + string(format="(' LEV=',f6.2)",line.slev) + string( format="('(',f6.2)",zheights[line.ilev]) + "KM)"
     endif else if field.grid_type eq 'geographic' and line.plotz gt 0 then begin
       toplab = toplab + string(format="(' LAT=',f6.2)",line.slat)
       if line.plotz eq 1 then toplab = toplab + string(format="(' GPALT LEV=',f6.2)",line.alt) + "KM"
       if line.plotz eq 2 then toplab = toplab + string(format="(' GMALT LEV=',f6.2)",line.alt) + "KM"
       if line.plotz eq 3 then toplab = toplab + string(format="(' ZGALT LEV=',f6.2)",line.alt) + "KM"
     endif else begin
       toplab = toplab + string(format="('  MLAT = ',f7.2)",line.slat)
       toplab = toplab + string(format="('  MLEV = ',f7.2)",line.slev)
     endelse
     
     if field.grid_type eq 'geographic' then begin
       xtitle = 'LONGITUDE (DEGREES GEOGRAPHIC)'
     endif else begin
       xtitle = 'LONGITUDE (DEGREES MAGNETIC)'
     endelse

     ytitle = info.ftype + ' ' + fieldlab ; field
;
;  For FFT output need to set labels and frequency axis
;     
     IF line.amplitude GT 0 THEN BEGIN

       toplab = toplab + ' WAVENUMBER AMPLITUDE'
       xtitle = 'WAVENUMBER'

       IF line.amplitude EQ 2 THEN BEGIN
         ytitle = ytitle + ' WAVENUMBER ' + ' AMPLITUDE LOG10'
       ENDIF ELSE IF line.amplitude EQ 1 THEN BEGIN
         ytitle = ytitle + ' WAVENUMBER ' + ' AMPLITUDE'
       ENDIF
;
;  Get wavenumbers(lons here) and wavenumber amplitude in equal size arrays for plotting
; 
       lons = indgen(nlon/2)
       lonminspec = line.lonmin
       lonmaxspec = line.lonmax
       IF line.lonmax EQ nlon THEN line.lonmax = line.lonmax/2-1

       dataspec = data
       data = fltarr(nlon/2)
       data[line.lonmin:line.lonmax] = dataspec[line.lonmin:line.lonmax]
       line.fmin = MIN(data[line.lonmin:line.lonmax])
       line.fmax = MAX(data[line.lonmin:line.lonmax])

     ENDIF
     	 
     plot, lons, data, $
         xrange=[line.lonmin,line.lonmax], xstyle=1, $ ; exact axis
         yrange=[line.fmin,line.fmax], ystyle=1, $ ; exact axis
         xlog=0, ylog=0, $
         charsize=1.0, $
         title=toplab, $
         xtitle = xtitle, $
         ytitle = ytitle, $
         subtitle = botlab1, $
         psym = -2, $
         xmargin = [10,7], $
         ymargin = [5,2]
;         ticklen=0, $
;         background = 255, $
;         color = 0
;         xtickformat='empty_string', ytickformat='empty_string'
;         min_value=NULL_VALUE
;         color=color
;         noerase=noerase
;         /nodata
;         position=self.plot1d_position, $

     IF line.amplitude GT 0 THEN data = dataspec
   
  end
   
  'time': begin
;
; Titles, labels, tic marks, etc:
;     
;     xtitle = field1D.tunits
     xtitle = 'Fractional Model Day'
     ytitle = info.ftype + ' ' + fieldlab ; field

     file = *info.pfile
     mtime = file.mtimes         
;
; Top label (under field label):
;
     if strlen(field1D.long_name) gt 0 then fieldname = field1D.long_name $
                                                 else fieldname = field1D.name

     timediff = line.timemax - line.timemin
;     daymin = FIX(line.timemin / 24. /60.)
;     daymax = FIX(line.timemax / 24. /60.)
     daymin = line.timemin / 24. /60. + mtime[0,0]
     daymax = line.timemax / 24. /60. + mtime[0,0]
;     jmin = ixfind_nearest(times,line.timemin)
;     jmax = ixfind_nearest(times,line.timemax)
;     IF daymin NE mtime[0,jmin] THEN daymin = mtime[0,jmin]
;     IF daymax NE mtime[0,jmax] THEN daymax = mtime[0,jmax]
     daydiff = daymax - daymin
     daymindiff = 24. * 60.

     IF timediff LT daymindiff THEN BEGIN						 
       toplab = string(format="('DAY=',i3)",daymin) + ' ' + fieldname + ' Model Time Series'
     ENDIF ELSE BEGIN
       toplab = string(format="('DAY=',f6.2)",daymin) + ' To' + $
                  string(format="(' DAY=',f6.2)",daymin+daydiff) + ' ' + fieldname + ' Model Time Series'
     ENDELSE
    
     uthrs = float(line.mtime[1])+float(line.mtime[2])/60.
     
     mdaydiff = 0
     mdaydiff = mtime[0,0]
;     mtimetimediff = mtime[0,0] - FIX(times[0]/1440)
;     IF mtimetimediff GT 0 THEN mdaydiff = mtimetimediff
     timesplt = times/1440.0 + mdaydiff
     timeminplt = line.timemin/1440.0 + mdaydiff
     timemaxplt = line.timemax/1440.0 + mdaydiff
     	 
     plot, timesplt, data, $
         xrange=[timeminplt,timemaxplt], xstyle=1, $ ; exact axis
         yrange=[line.fmin,line.fmax], ystyle=1, $ ; exact axis
         xlog=0, ylog=0, $
         charsize=1.0, $
         title=toplab, $
         xtitle = xtitle, $
         ytitle = ytitle, $
;	 xticks = ntticks, $
;         xtickname=timelabs, $
         subtitle = botlab1, $
         psym = -2, $
         xmargin = [10,7], $
         ymargin = [5,2]
;         ticklen=0, $
;         background = 255, $
;         color = 0
;         xtickformat='empty_string', ytickformat='empty_string'
;         min_value=NULL_VALUE
;         color=color
;         noerase=noerase
;         /nodata
;         position=self.plot1d_position, $
   
  end

  else: print,'>>> pltline: unknown plot type ',line.plottype

endcase ; line plot type

*info.pline = line

end