;
pro deflatdata,lat,data
;
; Define lat.data(nlon,nlev) from processed var in data.
;
;
;  If plotting with altitude vertical axis, make sure data and levels are on zp in
;  case data and levels are already on altitude which won't work in deflatdata
;
IF lat.plotz GE 1 THEN BEGIN

  *lat.data = *lat.datazp
  *lat.levs = *lat.levszp
  lat.levmin = MIN(*lat.levs)
  lat.levmax = MAX(*lat.levs)
;
;  When switching from zp to altitude with a user set vertical range need to reset
;  levmin and levmax so full range is used for altitude which is required by getzlatlonut
;  
  IF lat.ftype NE 'WACCM' THEN BEGIN
   lat.levmin = 0.0
   lat.levmax = 0.0
  ENDIF

ENDIF

field = lat.field
nlon = field.nlon & nlat = field.nlat & nlev = field.nlev
;
; Check dimensions (return if field is only 2d lat x lon)
;
dims = size(data,/dimensions)
ndims = n_elements(dims)
error = 0
if (dims[0] ne field.nlon) then begin
  print,'>>> WARNING deflatdata: 1st dimension is not nlons: dims[0]=',$
    dims[0],' but nlon=',field.nlon
  error = 1 
endif
if (dims[1] ne field.nlat) then begin
  print,'>>> WARNING deflatdata: 2nd dimension is not nlats: dims[1]=',$
    dims[1],' but nlat=',field.nlat
  error = 1 
endif
if (dims[2] ne field.nlev) then begin
  print,'>>> WARNING deflatdata: 3rd dimension is not nlats: dims[2]=',$
    dims[2],' but nlev=',field.nlev
  error = 1 
endif
if nlev le 0 then begin
  print,'>>> deflatdata: nlev=',nlev,' -- cannot make vertical slice.'
  error = 1 
endif
if error gt 0 then begin
  heap_free,lat.data
  return
endif
lons = *field.lons
lats = *field.lats
lonindices = indgen(nlon)
;
; Use limited x-axis range, if requested by user:
;
lonmin = lat.lonmin & lonmax = lat.lonmax
if lonmin gt lons[0] or lonmax lt lons[nlon-1] then begin
  nlonsnew = 0
  for j=0,nlon-1 do begin
    if lons[j] ge lonmin and lons[j] le lonmax then begin
      nlonsnew = nlonsnew+1 
    endif
  endfor 
  if nlonsnew le 2 then begin
    print,'>>> WARNING deflatdata: lonmin,max=',lonmin,lonmax,' nlonsnew=',$
      nlonsnew
    print,'    Cannot plot this lonmin,max -- resorting to full longitude ',$
      'grid..'
    nlonsnew = nlon
  endif else begin
    newlons = fltarr(nlonsnew)
    lonindnew = intarr(nlonsnew)
    jj = 0
    for j=0,nlon-1 do begin
      if lons[j] ge lonmin and lons[j] le lonmax then begin
        newlons[jj] = lons[j]
        lonindnew[jj] = j
        jj = jj+1
      endif
    endfor
    lons = newlons
    nlon = nlonsnew
    lonindices = lonindnew
  endelse
endif

fieldlevs = *field.levs
nfieldlevs = n_elements(fieldlevs)
;
;  Setting of min/max of data is done in interpolation routine so don't do here
;
IF lat.ftype EQ 'WACCM' THEN BEGIN

  levs = fieldlevs
  nlev = n_elements(levs)

ENDIF ELSE BEGIN
;
; If axis limits levmin,max are set, reset levs accordingly:
;
if lat.levmin lt lat.levmax then begin ; set axes limits levs:
  nlev = 0
  for k=0,nfieldlevs-1 do $
    if fieldlevs[k] ge lat.levmin and fieldlevs[k] le lat.levmax then nlev = nlev+1
  if nlev lt 2 then print,'>>> WARNING deflatdata: nlev point for y-axis must be',$
    ' at least 2: nlev=',nlev
  levs = fltarr(nlev)
  kk = 0
  for k=0,nfieldlevs-1 do $
    if fieldlevs[k] ge lat.levmin and fieldlevs[k] le lat.levmax then begin
      levs[kk] = fieldlevs[k]
      kk = kk+1
    endif
endif else levs = fieldlevs
nlev = n_elements(levs)

ENDELSE
;
; Define (nlon,nlev) at lat.slon:
;
latdata = fltarr(nlon,nlev)
latdata[*,*] = 0.
ilat = ixfind_nearest(lats,lat.slat)
if ilat lt 0 then begin
  print,'>>> deflatdata: cannot find latitude ',lat.slat,': nlat=',nlat,$
    ' ilat=',ilat,' lats=' & print,lats
  lat.data = ptr_new(latdata)
  return
endif
;
;  Setting of min/max of data is done in interpolation routine so don't do here for WACCM
;
IF lat.ftype EQ 'WACCM' THEN BEGIN

  for j=0,nlon-1 do begin
    latdata[j,*] = float(data[lonindices[j],ilat,*,lat.imtime])
  endfor

ENDIF ELSE BEGIN

  if lat.levmin lt lat.levmax then begin 
    kk = 0
    for k=0,nfieldlevs-1 do begin
      if fieldlevs[k] ge lat.levmin and fieldlevs[k] le lat.levmax then begin
        for j=0,nlon-1 do begin
          latdata[j,kk] = float(data[lonindices[j],ilat,k,lat.imtime])
        endfor
        kk = kk+1
      endif
    endfor
  endif else begin

    for j=0,nlon-1 do begin
      latdata[j,*] = float(data[lonindices[j],ilat,*,lat.imtime])
    endfor

  endelse

ENDELSE
;
; Take log10 of densities:
;
if (lat.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
   (lat.log10 eq 'this frame only') or $
   (lat.log10 eq 'all fields') then begin
  missing_value = float(lat.missing_value)
  log10f,latdata,missing_value
endif
;
; Update lon structure:
;
lat.slat = lats[ilat]
lat.data = ptr_new(latdata)
lat.lons = ptr_new(lons)
lat.levs = ptr_new(levs)
;
;  Interpolate WACCM data to 'standard' log pressure levels.  Already done if plotting on altitude
;  in vertical so skip in that case.
;
IF lat.ftype EQ 'WACCM' THEN BEGIN

  if lat.plotz eq 0 then interpwdata, lat

ENDIF
;
;  This is part of handling the case where the vertical range was set by the user then
;  the vertical grid was changed from zp to altitude
;
IF (lat.plotz GE 0 AND lat.levmin EQ 0 AND lat.levmax EQ 0 AND lat.ftype NE 'WACCM') THEN *lat.levszp = *lat.levs

end