;
pro pltlon,info,animate=animate,ps=ps,png=png
;
; Plot longitude slice (y:x = zp:lat):
;
file = *info.pfile
lon = *info.plon
field = lon.field
fields = *info.fields
flats = *fields[lon.ifield].lats
nflats = n_elements(flats)
data = *lon.data
lats = *lon.lats & nlat = n_elements(lats)
levs = *lon.levs 
nlev = n_elements(levs)
;
rhyaxis = lon.rhyaxis
if rhyaxis gt 0 and (file.diffs gt 0 or field.grid_type eq 'magnetic') then $
  rhyaxis = 0
latrangeset = 0
if lon.latmin gt flats[0] or lon.latmax lt flats[nflats-1] then $
  latrangeset = 1

find_nans,info,field.name,data,numnans
if numnans gt 0 then print,'>>NaNs found:',numnans
if numnans gt 0 then  begin
  rhyaxis = 0   
  lon.rhyaxis = 0
endif
;
if not keyword_set(animate) and not keyword_set(ps) and $
   not keyword_set(png) then begin
  widget_control,lon.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(lon.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,'Lon slices: Field ',field.name,': Found ',$
  nmissing,' missing data points (out of ',n_elements(data),')'
;
; Field min,max:
; If lon.fmin >  lon.fmax, then leave as such for the animator but use 
;                              field extremes for plotting current frame.
; If lon.fmin =  lon.fmax, then reset to data extremes (for the slice).
; If lon.fmin <  lon.fmax, then leave as such as for plotting.
;
fminmax,data,framemin,framemax,missing_value

if lon.fmin ge lon.fmax then begin  ; lon.fmin >= lon.fmax
  fminmax,data,fmin,fmax,missing_value
  if lon.fmin eq lon.fmax then begin
    lon.fmin=fmin & lon.fmax=fmax
    print,'pltlon reset lon.fmin,max=',lon.fmin,lon.fmax
  endif
endif else begin                        ; lon.fmin < lon.fmax
  fmin = lon.fmin & fmax = lon.fmax
endelse

if lon.setcontour eq 'off' and lon.plotz eq 0 then begin
  cmin=0. & cmax=0. & cint=0.
  deflondata,lon,*fields[lon.ifield].data
;
;  Set local levels and number of levels variable since changed in deflondata for WACCM
;
  IF info.ftype EQ 'WACCM' THEN BEGIN

    levs = *lon.levs 
    nlev = n_elements(levs)

  ENDIF
  
  fminmax,*lon.data,fmin,fmax,lon.missing_value
  levels = mkclevels(fmin,fmax,cmin,cmax,cint)
  lon.fmin=cmin & lon.fmax=cmax & lon.cint=cint
  *lon.clevels = levels
  *info.plon =lon 
  data = *lon.data
endif
;
;  For altitude vertical coordinate and 1D/2D wavenumber analysis, get min and max levels of data
;
if lon.plotz gt 0 and (lon.amplitude gt 0 or lon.zonalfreq gt 0) then begin
  cmin=0. & cmax=0. & cint=0.
  fminmax,*lon.data,fmin,fmax,lon.missing_value
  levels = mkclevels(fmin,fmax,cmin,cmax,cint)
  lon.fmin=cmin & lon.fmax=cmax & lon.cint=cint
  *lon.clevels = levels
  *info.plon =lon 
  data = *lon.data
endif
;
; Image position:
;
xur = .90 & yur = .90
xll = .15 & yll = .31
if rhyaxis gt 0 then begin 
  xur=xur-.025 & xll=xll-.025 
endif

axespos  = [xll,yll,xur,yur]
;
;  For magnetic fields need to move axes to fit contour
;
IF field.grid_type EQ 'magnetic' THEN axespos = [xll-.002,yll,xur-.002,yur]

xmid = .5*(xur+xll)
;
; Gap between pole and first latitude (>0 only at top latitude on geographic grid) :
;
polegap_south = (90.-abs(lats[0]))/180.*(xur-xll)
if lon.latmin gt flats[0] then polegap_south = 0.
polegap_north = (90.-abs(lats[nlat-1]))/180.*(xur-xll)
if lon.latmax lt flats[nflats-1] then polegap_north = 0.
imagepos = [xll+polegap_south,yll,xur-polegap_south,yur]
;print,'pltlon: polegap_south=',polegap_south,' polegap_north=',polegap_north
;
; Byte-scale to an image. Make any missing data black:
;
if (lon.plottype eq 'image_only') or (lon.plottype eq 'image+contours') then begin
  missing_indices = where(data eq missing_value,nmissing)
  if nmissing gt 0 then begin
    print,'Field ',field.name,' lon slice: Found ',nmissing,' missing data points.'
    image = bytscl(data,min=fmin,max=fmax,top=!d.table_size-3)+1
    image[missing_indices] = 0
  endif else begin
    image = bytscl(data,min=fmin,max=fmax,top=!d.table_size-3)+1
  endelse
  if keyword_set(ps) eq 0 then begin
    xyll = convert_coord([xll+polegap_south],[yll],/norm,/to_device)
    xll_dev = xyll(0,0) & yll_dev = xyll(1,0)
    xyur = convert_coord([xur-polegap_north],[yur],/norm,/to_device)
    xur_dev = xyur(0,0) & yur_dev = xyur(1,0)
    image=congrid(image,xur_dev-xll_dev+1,yur_dev-yll_dev+1,/interp,/minus_one)
    xsize=0. & ysize=0.
  endif else begin ; ps
    imscale = 8
    image=congrid(image,nlat*imscale,nlev*imscale,/interp,/minus_one)
    xsize=(xur-polegap_north)-(xll+polegap_south) & ysize=yur-yll
  endelse
endif ; making an image
;
; Set up for color bar:
;
if lon.plottype ne 'monochrome_contours' then begin
  barw = xur-xll-.05*(xur-xll) & barh = .03
  xbar = xmid-.5*barw
  ybar = yll-barh-.11
  baroff = .035
  cblabels = [strcompress(string(fmin,format="(g9.3)"),/remove_all),$
    strcompress(string(0.5*(fmin+fmax),format="(g9.3)"),/remove_all),$
    strcompress(string(fmax,format="(g9.3)"),/remove_all)]
endif
;
; Titles, labels, tic marks, etc:
;
if field.grid_type eq 'geographic' then begin
  xtitle = 'LATITUDE (DEGREES GEOGRAPHIC)'
endif else begin
  xtitle = 'LATITUDE (DEGREES MAGNETIC)'
endelse

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 lon.plotz EQ 1 THEN ytitle = ' GEOPOTENTIAL HEIGHT (km) '
IF lon.plotz EQ 2 THEN ytitle = ' GEOMETRIC HEIGHT (km) '

xticlabs = ['-90','-60','-30','0','30','60','90']
;
; Field label:
;
if (lon.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
   (lon.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 lon.amplitude gt 0 then begin
  IF lon.amplitude EQ 1 THEN BEGIN
    if strlen(field.long_name) gt 0 then fieldlab = field.long_name $
      else fieldlab = field.name
  ENDIF ELSE IF lon.amplitude EQ 2 THEN BEGIN
    if strlen(field.long_name) gt 0 then fieldlab = 'LOG10 '+field.long_name $
      else fieldlab = 'LOG10 '+field.name
  ENDIF
endif

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) + ')'
;
; Top label (under field label):
;
uthrs = float(lon.mtime[1])+float(lon.mtime[2])/60.
toplab = string(format="('DAY = ',i3,' UT = ',f5.2)",lon.mtime[0],uthrs)
if lon.zonalmean gt 0 then begin
  toplab = toplab + ' ZONAL MEAN'
endif else if lon.amplitude gt 0 then begin
  IF lon.amplitude EQ 1 THEN BEGIN
    toplab = toplab + ' WAVENUMBER ' + string(format="(i2)",lon.wavenumber) + ' AMPLITUDE'
  ENDIF ELSE IF lon.amplitude EQ 2 THEN BEGIN
      toplab = toplab + ' WAVENUMBER ' + string(format="(i2)",lon.wavenumber) + ' AMPLITUDE(LOG10)'
  ENDIF
endif else if lon.zonalfreq gt 0 then begin
  IF lon.zonalfreq EQ 1 THEN BEGIN
    toplab = toplab + ' WAVENUMBER ' + string(format="(i2)",lon.wavenumber) + ' FREQUENCY ' + $
                                                  string(format="(i3)",lon.frequency) + ' AMPLITUDE'
  ENDIF ELSE IF lon.zonalfreq EQ 2 THEN BEGIN
    toplab = toplab + ' WAVENUMBER ' + string(format="(i2)",lon.wavenumber) + ' FREQUENCY ' + $
                                                  string(format="(i3)",lon.frequency) + ' AMPLITUDE(LOG10)'
  ENDIF ELSE IF lon.zonalfreq EQ 3 THEN BEGIN
    toplab = toplab + ' WAVENUMBER ' + string(format="(i2)",lon.wavenumber) + ' FREQUENCY ' + $
                                                  string(format="(i3)",lon.frequency) + ' PHASE'
  ENDIF ELSE IF lon.zonalfreq EQ 4 THEN BEGIN
    toplab = toplab + ' WAVENUMBER ' + string(format="(i2)",lon.wavenumber) + ' FREQUENCY ' + $
                                                  string(format="(i3)",lon.frequency) + ' PHASE(LOG10)'
  ENDIF
endif else if field.grid_type eq 'geographic' then begin
  toplab = toplab + string(format="('  LON = ',f7.2)",lon.slon)
  toplab = toplab + string(format="('  SLT = ',f5.2)",lon.sslt)
endif else begin
  toplab = toplab + string(format="('  MLON = ',f7.2)",lon.slon)
  toplab = toplab + string(format="('  MLT = ',f5.2)",lon.sslt)
endelse
;
; Bottom label:
;
 botlab1 = strcompress(field.ncfile,/remove_all)
;botlab1 = strcompress(field.ncfile,/remove_all)+ $
;  string(format="(' (',i3,',',i2,',',i2,')')",lon.mtime[*])
;
; Min,max label:
;
minmaxlab = 'MIN,MAX = '+string(format="(g10.4)",framemin)+','+$
                         string(format="(g10.4)",framemax)
if lon.plottype ne 'image_only' then $
  minmaxlab = minmaxlab + string(format="(' INTERVAL = ',g10.4)",lon.cint)
;
; Contour levels were made by mkclevels.pro.
;
clevels = *lon.clevels
nlevels = n_elements(clevels)

clinestyle = intarr(nlevels)
clinestyle[*] = 0 ; solid
for i=0,nlevels-1 do if clevels[i] lt 0. then clinestyle[i] = 2 ; dashed
conlab = 'Contour from '+string(format="(g10.4)",lon.fmin)+ $
         ' to '         +string(format="(g10.4)",lon.fmax)+ $
         ' Interval = ' +string(format="(g10.4)",lon.cint)
c_charsize = 1.1 ; contour line label character size
chsize = 1.2     ; info label character size
irregular = 0 & if field.grid_type eq 'magnetic' then irregular = 1
;
; 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

;if lon.fmin eq lon.missing_value and lon.fmax eq lon.missing_value then begin
;  lon.plottype = 'image_only'
;endif 
;
; Display image and/or contour:
;
case lon.plottype of
;
; Display image only. Use contour,/nodata to draw axes, ticmarks, labels:
;
  'image_only': begin
    tv,image,xll+polegap_south,yll,/norm,xsize=xsize,ysize=ysize       ; ras image
    clrbar,barw,barh,xbar,ybar,cblabels,botcolor=1,topcolor=!d.table_size-2,$
      charsize=chsize,laboffset=baroff
    if latrangeset then $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$       ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        charsize=chsize $
    else $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$       ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,charsize=chsize
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.21,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.25,botlab1  ,/norm,align=0.5,charsize=chsize
  end
;
; Display image with contour overlay:
;
  'image+contours': begin
    tv,image,xll+polegap_south,yll,/norm,xsize=xsize,ysize=ysize 
    contour,data,lats,levs,/noerase,position=imagepos,$              ; contour
      /follow,xstyle=5,ystyle=5,c_charsize=c_charsize,$
      levels=clevels,c_linestyle=clinestyle,max_value=missing_value,$
      c_colors=[lon.clineclr]
    if latrangeset then $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$       ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        charsize=chsize $
    else $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$       ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,charsize=chsize
    clrbar,barw,barh,xbar,ybar,cblabels,botcolor=1,topcolor=!d.table_size-2,$
      charsize=chsize,laboffset=baroff
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.21,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.25,botlab1  ,/norm,align=0.5,charsize=chsize
  end
;
; Contour only (B&W):
;
  'monochrome_contours': begin
    if latrangeset then $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$       ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        charsize=chsize $
    else $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$       ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,charsize=chsize
    contour,data,lats,levs,/noerase,position=imagepos,$              ; contour
      /follow,xstyle=5,ystyle=5,c_charsize=c_charsize,$
      levels=clevels,c_linestyle=clinestyle,max_value=missing_value
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.12,conlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.16,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.20,botlab1,/norm,align=0.5,charsize=chsize
  end
;
; Color fill contour:
;
  'colorfill_contours': begin
    if latrangeset then $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$        ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        charsize=chsize $
    else $
      contour,data,lats,levs,/noerase,/nodata,position=axespos,$        ; axes
        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=yaxstyle,ticklen=-.015,$
        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,charsize=chsize

    contour,data,lats,levs,/noerase,position=imagepos,/follow,/fill,$   ; colorfill
      xstyle=5,ystyle=5,levels=clevels,max_value=missing_value

    contour,data,lats,levs,/noerase,position=imagepos,/follow,$         ; contour
      xstyle=5,ystyle=5,c_charsize=c_charsize,levels=clevels,$
      c_linestyle=clinestyle,max_value=missing_value,c_colors=[lon.clineclr]

    clrbar,barw,barh,xbar,ybar,cblabels,botcolor=1,topcolor=!d.table_size-2,$
      charsize=chsize,laboffset=baroff
    xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.21,conlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.24,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.27,botlab1,/norm,align=0.5,charsize=chsize
  end
  else: print,'>>> pltlon: unknown plot type ',lon.plottype
endcase ; plot type

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

if rhyaxis gt 0 and numnans eq 0 and lon.plotz eq 0 then begin
  altyax,info,lon,'lon',xll,xur,yll,yur,chsize 
endif
;
; lon.fmin,max may have been changed, so update state:
;
*info.plon = lon
end