;
pro deflondata,lon,data

;latmin = -30. & latmax = 30. ; lat limits for x-axis
;
; Define lon.data(nlat,nlev) from processed var in data.
;
field = lon.field
nlon = field.nlon & nlat = field.nlat & nlev = field.nlev
;print,' '
;print,'deflondata: Field ',field.name,' nlat=',nlat,' nlev=',nlev
;print,'  slon=',lon.slon,' sslt=',lon.sslt,' ut=',lon.ut
;
; 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 deflondata: 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 deflondata: 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 deflondata: 3rd dimension is not nlats: dims[2]=',$
    dims[2],' but nlev=',field.nlev
  error = 1 
endif
if nlev le 0 then begin
  print,'>>> deflondata: nlev=',nlev,' -- cannot make vertical slice.'
  error = 1 
endif
if error gt 0 then begin
  size = size(lon.data)
  londata = fltarr(size[1],size[2])
  londata[*,*] = float(lon.missing_value)
  lon.data = ptr_new(londata)
  return
endif
lons = *field.lons
lats = *field.lats
latindices = indgen(nlat)
;
; Use limited x-axis range, if requested by user:
;
latmin = lon.latmin & latmax = lon.latmax
if latmin gt lats[0] or latmax lt lats[nlat-1] then begin
  nlatsnew = 0
  for j=0,nlat-1 do begin
    if lats[j] ge latmin and lats[j] le latmax then begin
      nlatsnew = nlatsnew+1 
;     print,'deflondata: j=',j,' latmin,max=',latmin,latmax,' lats[j]=',lats[j],$
;       ' nlatsnew=',nlatsnew
    endif
  endfor 
  if nlatsnew le 2 then begin
    print,'>>> WARNING deflondata: latmin,max=',latmin,latmax,' nlatsnew=',$
      nlatsnew
    print,'    Cannot plot this latmin,max -- resorting to full latitude ',$
      'grid..'
    nlatsnew = nlat
  endif else begin
    newlats = fltarr(nlatsnew)
    latindnew = intarr(nlatsnew)
    jj = 0
    for j=0,nlat-1 do begin
      if lats[j] ge latmin and lats[j] le latmax then begin
        newlats[jj] = lats[j]
        latindnew[jj] = j
        jj = jj+1
      endif
    endfor
    lats = newlats
    nlat = nlatsnew
    latindices = latindnew
;   print,'deflondata limit latitude x-axis: nlat=',nlat,' lats=' & print,lats
;   print,'latindices=' & print,latindices
  endelse
endif
;
fieldlevs = *field.levs
nfieldlevs = n_elements(fieldlevs)
;print,'deflondata: nfieldlevs=',nfieldlevs,' fieldlevs=' & print,fieldlevs
;
; If axis limits levmin,max are set, reset levs accordingly:
;
if lon.levmin lt lon.levmax then begin ; set axes limits levs:
; print,'deflondata: lon.levmin,max=',lon.levmin,lon.levmax
  nlev = 0
  for k=0,nfieldlevs-1 do $
    if fieldlevs[k] ge lon.levmin and fieldlevs[k] le lon.levmax then nlev = nlev+1
  if nlev lt 2 then print,'>>> WARNING pltlon: 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 lon.levmin and fieldlevs[k] le lon.levmax then begin
      levs[kk] = fieldlevs[k]
      kk = kk+1
    endif
endif else levs = fieldlevs
nlev = n_elements(levs)
;
;print,'deflondata: nfieldlevs=',nfieldlevs,' nlev=',nlev,' levs=' & print,levs
;
; Define (nlat,nlev) at lon.slon:
;
londata = fltarr(nlat,nlev)
londata[*,*] = 0.
ilon = -1
for i=0,nlon-1 do begin
  if lons[i] le lon.slon+.01 and lons[i] ge lon.slon-.01 then begin
    ilon = i 
    goto,found
  endif
endfor
found:
if ilon lt 0 then begin
  print,'>>> deflondata: cannot find longitude ',lon.slon,': nlon=',nlon,$
    ' ilon=',ilon,' lons=' & print,lons
  lon.data = ptr_new(londata)
  return
endif
;
if lon.levmin lt lon.levmax then begin 
  kk = 0
  for k=0,nfieldlevs-1 do $
    if fieldlevs[k] ge lon.levmin and fieldlevs[k] le lon.levmax then begin
      if lon.zonalmean eq 0 then begin
        for j=0,nlat-1 do begin
          londata[j,kk] = float(data[ilon,latindices[j],k,lon.imtime])
        endfor
      endif else begin ; plot zonal mean
        for j=0,nlat-1 do begin
          londata[j,kk] = float(mean(data[*,latindices[j],k,lon.imtime]))
        endfor
      endelse
      kk = kk+1
    endif
endif else begin
  if lon.zonalmean eq 0 then begin
    for j=0,nlat-1 do begin
      londata[j,*] = float(data[ilon,latindices[j],*,lon.imtime])
    endfor
  endif else begin ; zonal mean
    for j=0,nlat-1 do begin
      for k=0,nlev-1 do begin
        londata[j,k] = float(mean(data[*,latindices[j],k,lon.imtime]))
      endfor
    endfor
  endelse ; zonal mean
endelse
;
; Take log10 of densities:
;
if (lon.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
   (lon.log10 eq 'this frame only') or $
   (lon.log10 eq 'all fields') then begin
; print,'deflondata: taking log10 of field ',field.name,'...'
  missing_value = float(lon.missing_value)
  log10f,londata,missing_value
endif
;
; Update lon structure:
;
lon.data = ptr_new(londata)
lon.lats = ptr_new(lats)
lon.levs = ptr_new(levs)
end