;
pro deflondata,lon,data
;
; Define lon.data(nlat,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 deflondata
;
IF lon.plotz GE 1 THEN BEGIN

  *lon.data = *lon.datazp
  *lon.levs = *lon.levszp
  lon.levmin = MIN(*lon.levs)
  lon.levmax = MAX(*lon.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 lon.ftype NE 'WACCM' THEN BEGIN
   lon.levmin = 0.0
   lon.levmax = 0.0
  ENDIF
  
ENDIF

field = lon.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 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 nlevs: 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
  heap_free,lon.data
  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 
    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
  endelse
endif
;
;  Pull out vertical levels to check for user specified vertical range
;
fieldlevs = *field.levs
nfieldlevs = n_elements(fieldlevs)
;
;  Setting of min/max of data is done in interpolation routine so don't do here
;
IF lon.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 lon.levmin lt lon.levmax then begin ; set axes limits levs:
    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

ENDELSE
;
; 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
;
;  Load data and do zonal mean and wave amplitude calculations if requested
;  WACCM is handled in a different way than the others and check for subsetted
;  axes ranges done later in interpwdata for WACCM
;
IF lon.ftype EQ 'WACCM' THEN BEGIN
;
;  Handle wave amplitude calculation here
;
  IF lon.amplitude NE 0 AND field.name NE 'YAXZ' THEN BEGIN
  
    iwvn = lon.wavenumber
    ampdata = fltarr(nlon)
    ampmin = 1.0E+24
    for j=0,nlat-1 do begin
      for k=0,nlev-1 do begin
        
        ampdata[*] = float(data[*,latindices[j],k,lon.imtime])
        ampdata[*] = ABS(FFT(ampdata[*], -1)) * 2.
	if ampdata[iwvn] LT ampmin AND ampdata[iwvn] GT 0.0 then ampmin = ampdata[iwvn]
	londata[j,k] = ampdata[iwvn]
    
      endfor ;nlev
    endfor ;nlat

    IF lon.amplitude EQ 2 THEN BEGIN
;
;  Check for zero values and set to a small number then take log of data if requested
;
      for j=0,nlat-1 do begin

        zeroind = WHERE(londata[j,*] EQ 0.0)
        if zeroind[0] NE -1 then londata[j,zeroind] = ampmin
      
      endfor

      londata[*,*] = ALOG10(londata[*,*])  

    ENDIF ;Log of amplitude
;
;  Handle wave amplitude and phase 2D calculation here
;
  ENDIF ELSE IF lon.zonalfreq NE 0 AND field.name NE 'YAXZ' THEN BEGIN
;
;  If amplitude and phase already calculated, just load from saved arrays
; 
    IF lon.ampphs EQ 1 THEN BEGIN   
;
;  Handle user specified vertical range
;      
      for j = 0,nlat-1 do begin 	     
        IF lon.zonalfreq EQ 1 OR lon.zonalfreq EQ 2 THEN BEGIN
	  amp = *lon.amp
          londata[j,*] = REFORM(amp[latindices[j],*])
          IF lon.zonalfreq EQ 2 THEN londata[j,*] = ALOG10(londata[j,*])
        ENDIF ELSE IF lon.zonalfreq EQ 3 OR lon.zonalfreq EQ 4 THEN BEGIN
	  phs = *lon.phs
          londata[j,*] = REFORM(phs[latindices[j],*])
          IF lon.zonalfreq EQ 4 THEN londata[j,*] = ALOG10(londata[j,*])
        ENDIF 
      endfor 
 
    ENDIF ELSE BEGIN ;2D FFT not done for this field
    
      iwvn = lon.wavenumber
      ifreq = lon.frequency
      ntime = field.ntime
      ntime2 = ntime/2+1
      nshift = (nlon+1)/2-1
      uspcmin = 1.0E+24
      lonshift = lons[0] / (lons[1] - lons[0])
;
;  Need to pick out correct westward(negative)/eastward(positive) wavenumber since 0 is at center of array
;  and set frequency index to positive if specified as negative
;      
      iwvn = nshift - iwvn
;      IF ifreq LT 0 THEN iwvn = nshift + iwvn
;      IF ifreq GE 0 THEN iwvn = nshift - iwvn
;      IF ifreq LT 0 THEN ifreq = -ifreq
;
;  Do 2D FFT to get amplitude and phase
; 
      uspc = fltarr(nlon,ntime2,field.nlat,nlev)
      uphs = fltarr(nlon,ntime2,field.nlat,nlev)
      for k = 0,nlev-1 do begin
        for j = 0,nlat-1 do begin
          tmpa = fltarr(nlon,ntime)
          for it = 0,ntime-1 do begin
            for i = 0,nlon-1 do begin
              tmpa(i,it)=data(i,j,k,it)
            endfor
            tmpa(*,it) = shift(tmpa(*,it),lonshift) ;;this is to put 0 deg at the first point
          endfor
          fftmp = fft(tmpa)

          for it = 0,ntime2-1 do begin
            uspc(*,it,j,k) = shift(abs(fftmp(*,it)),nshift)*2.
            uphs(*,it,j,k) = -shift(atan(fftmp(*,it),/phase),nshift)*12./!pi
	    if min(uspc(*,it,j,k)) LT uspcmin AND min(uspc(*,it,j,k)) GT 0.0 then uspcmin = min(uspc(*,it,j,k))
          endfor

        endfor  	  
      endfor

      uphs(where(uphs lt 0))=uphs(where(uphs lt 0))+24.
      
      uspc(*,0,*,*)=uspc(*,0,*,*)/2.

      uphs(where(uphs le 0)) = min(uphs(where(uphs gt 0)))
;
;  Fill output data array
;     
      IF lon.zonalfreq EQ 1 OR lon.zonalfreq EQ 2 THEN BEGIN

        for j=0,nlat-1 do begin
          londata[j,*] = REFORM(uspc[iwvn,ifreq,j,*])
        endfor

	IF lon.zonalfreq EQ 2 THEN BEGIN
;
;  Check for zero values and set to a small number then take log of data if requested
;
          for j=0,nlat-1 do begin

            zeroind = WHERE(londata[j,*] EQ 0.0)
            if zeroind[0] NE -1 then londata[j,zeroind] = uspcmin
      
          endfor

	  londata[*,*] = ALOG10(londata[*,*])

        ENDIF

      ENDIF ELSE IF lon.zonalfreq EQ 3 OR lon.zonalfreq EQ 4 THEN BEGIN

        for j=0,nlat-1 do begin
          londata[j,*] = REFORM(uphs[iwvn,ifreq,j,*])
        endfor
	IF lon.zonalfreq EQ 4 THEN londata[*,*] = ALOG10(londata[*,*])

      ENDIF ;amplitude or phase
;
;  Save amplitude and phase so don't have to recalculate and set flag
; 
      lon.amp = ptr_new(REFORM(uspc[iwvn,ifreq,*,*]))
      lon.phs = ptr_new(REFORM(uphs[iwvn,ifreq,*,*]))

      lon.ampphs = 1 ;Value of 1 means 2D FFT already calculated
 
    ENDELSE; Previously calculated 2D FFT for current field
;
;  Case where just need to load data into output array
;
  ENDIF ELSE BEGIN ;no spectral analysis

    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
;
; 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') AND (lon.amplitude EQ 0) then begin
      missing_value = float(lon.missing_value)
      log10f,londata,missing_value
    endif

  ENDELSE ;spectral analysis or not

ENDIF ELSE BEGIN; Not WACCM
;
;  Vertical level minimum/maximum values already set 
;
  if lon.levmin lt lon.levmax then begin
;
;  Handle wave amplitude calculation here
;  
    IF lon.amplitude NE 0 AND field.name NE 'YAXZ' THEN BEGIN
    
      kk = 0
      iwvn = lon.wavenumber
      ampdata = fltarr(nlon)
      ampmin = 1.0E+24
      for k=0,nfieldlevs-1 do begin
;
;  Handle user specified vertical range
;      
  	if fieldlevs[k] ge lon.levmin and fieldlevs[k] le lon.levmax then begin

  	  for j=0,nlat-1 do begin
  	  
  	    ampdata[*] = float(data[*,latindices[j],k,lon.imtime])
  	    ampdata[*] = ABS(FFT(ampdata[*], -1)) * 2.
  	    if ampdata[iwvn] LT ampmin AND ampdata[iwvn] GT 0.0 then ampmin = ampdata[iwvn]
  	    londata[j,kk] = ampdata[iwvn]
  	   
  	  endfor
  
  	  kk = kk+1
  	  
  	endif ;vertical range
      
      endfor
;
;  Check for zero values and set to a small number then take log of data if requested
; 
      for j=0,nlat-1 do begin
  
  	zeroind = WHERE(londata[j,*] EQ 0.0)
  	if zeroind[0] NE -1 then londata[j,zeroind] = ampmin
  	
      endfor

      IF lon.amplitude EQ 2 THEN londata[*,*] = ALOG10(londata[*,*])
;
;  Handle wave amplitude and phase 2D calculation here
;  
    ENDIF ELSE IF lon.zonalfreq NE 0 AND field.name NE 'YAXZ' THEN BEGIN
;
;  If amplitude and phase already calculated, just load from saved arrays
; 
      IF lon.ampphs EQ 1 THEN BEGIN ;2D FFT already done for this field
 
        kk = 0
        for k = 0,nfieldlevs-1 do begin      
;
;  Handle user specified vertical range
; 
  	  if fieldlevs[k] ge lon.levmin and fieldlevs[k] le lon.levmax then begin
            for j = 0,nlat-1 do begin              
              IF lon.zonalfreq EQ 1 OR lon.zonalfreq EQ 2 THEN BEGIN
	        amp = *lon.amp
                londata[j,kk] = REFORM(amp[latindices[j],k])
	        IF lon.zonalfreq EQ 2 THEN londata[j,kk] = ALOG10(londata[j,kk])
              ENDIF ELSE IF lon.zonalfreq EQ 3 OR lon.zonalfreq EQ 4 THEN BEGIN
	        phs = *lon.phs
                londata[j,kk] = REFORM(phs[latindices[j],k])
	        IF lon.zonalfreq EQ 4 THEN londata[j,kk] = ALOG10(londata[j,kk])
              ENDIF	
            endfor 
  	    kk = kk+1
          endif 
	  
	endfor            
 
      ENDIF ELSE BEGIN ;2D FFT not done for this field   
    
        iwvn = lon.wavenumber
        ifreq = lon.frequency
        ntime = field.ntime
        ntime2 = ntime/2+1
        nshift = (nlon+1)/2-1
        lonshift = lons[0] / (lons[1] - lons[0])
;
;  Need to pick out correct westward(negative)/eastward(positive) wavenumber since 0 is at center of array
;  and set frequency index to positive if specified as negative
;      
        iwvn = nshift - iwvn
;        IF ifreq LT 0 THEN iwvn = nshift + iwvn
;        IF ifreq GE 0 THEN iwvn = nshift - iwvn
;        IF ifreq LT 0 THEN ifreq = -ifreq
;
;  Do 2D FFT to get amplitude and phase
;       
        uspc = fltarr(nlon,ntime2,field.nlat,nfieldlevs)
        uphs = fltarr(nlon,ntime2,field.nlat,nfieldlevs)
        for k = 0,nfieldlevs-1 do begin      

            for j = 0,nlat-1 do begin
              tmpa = fltarr(nlon,ntime)
              for it = 0,ntime-1 do begin
                for i = 0,nlon-1 do begin
                  tmpa(i,it)=data(i,j,k,it)
                endfor
                tmpa(*,it) = shift(tmpa(*,it),lonshift) ;;this is to put 0 deg at the first point
              endfor
              fftmp = fft(tmpa)

              for it = 0,ntime2-1 do begin
                uspc(*,it,j,k) = shift(abs(fftmp(*,it)),nshift)*2.
                uphs(*,it,j,k) = -shift(atan(fftmp(*,it),/phase),nshift)*12./!pi
              endfor
            endfor
  
        endfor

        uphs(where(uphs lt 0))=uphs(where(uphs lt 0))+24.

        uspc(*,0,*,*)=uspc(*,0,*,*)/2.

        uphs(where(uphs le 0)) = min(uphs(where(uphs gt 0)))
;
;  Fill output data array
;           
        kk = 0
        for k = 0,nfieldlevs-1 do begin      
;
;  Handle user specified vertical range
;      
 	  if fieldlevs[k] ge lon.levmin and fieldlevs[k] le lon.levmax then begin

            FOR j = 0, nlat-1 DO BEGIN

              IF lon.zonalfreq EQ 1 OR lon.zonalfreq EQ 2 THEN BEGIN
                londata[j,kk] = REFORM(uspc[iwvn,ifreq,j,k])
	        IF lon.zonalfreq EQ 2 THEN londata[j,kk] = ALOG10(londata[j,k])
              ENDIF ELSE IF lon.zonalfreq EQ 3 OR lon.zonalfreq EQ 4 THEN BEGIN
                londata[j,kk] = REFORM(uphs[iwvn,ifreq,j,k])
	        IF lon.zonalfreq EQ 4 THEN londata[j,kk] = ALOG10(londata[j,k])
              ENDIF	

            ENDFOR

  	    kk = kk+1
  	  
          endif
      
        endfor
;
;  Save amplitude and phase so don't have to recalculate and set flag
;  
        lon.amp = ptr_new(REFORM(uspc[iwvn,ifreq,*,*]))
        lon.phs = ptr_new(REFORM(uphs[iwvn,ifreq,*,*]))

        lon.ampphs = 1 ;Value of 1 means 2D FFT already calculated
 
      ENDELSE; Previously calculated 2D FFT for current field
;
;  Case where just need to load data into output array (no spectral analysis requested)
;
    ENDIF ELSE BEGIN ;No spectral analysis
 
      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
;
; 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') AND (lon.amplitude EQ 0) then begin
        missing_value = float(lon.missing_value)
        log10f,londata,missing_value
      endif

    ENDELSE ;Spectral analysis or not

  endif else begin  ;  Vertical level minimum/maximum values not set 
;
;  Handle wave amplitude calculation here
;  
    IF lon.amplitude NE 0 AND field.name NE 'YAXZ' THEN BEGIN
    
      iwvn = lon.wavenumber
      ampdata = fltarr(nlon)
      ampmin = 1.0E+24
      for j=0,nlat-1 do begin
  	for k=0,nlev-1 do begin
  	  
  	  ampdata[*] = float(data[*,latindices[j],k,lon.imtime])
  	  ampdata[*] = ABS(FFT(ampdata[*], -1)) * 2.
  	  if ampdata[iwvn] LT ampmin AND ampdata[iwvn] GT 0.0 then ampmin = ampdata[iwvn]
  	  londata[j,k] = ampdata[iwvn]
      
  	endfor
      endfor
  
      for j=0,nlat-1 do begin
  
  	zeroind = WHERE(londata[j,*] EQ 0.0)
  	if zeroind[0] NE -1 then londata[j,zeroind] = ampmin
  	
      endfor

      IF lon.amplitude EQ 2 THEN londata[*,*] = ALOG10(londata[*,*])
      
    ENDIF ELSE IF lon.zonalfreq NE 0 AND field.name NE 'YAXZ' THEN BEGIN
;
;  If amplitude and phase already calculated, just load from saved arrays
;      
      IF lon.ampphs EQ 1 THEN BEGIN
;
;  Handle user specified vertical range
;      
        for j = 0,nlat-1 do begin	       
          IF lon.zonalfreq EQ 1 OR lon.zonalfreq EQ 2 THEN BEGIN
	    amp = *lon.amp
            londata[j,*] = REFORM(amp[latindices[j],*])
            IF lon.zonalfreq EQ 2 THEN londata[j,*] = ALOG10(londata[j,*])
          ENDIF ELSE IF lon.zonalfreq EQ 3 OR lon.zonalfreq EQ 4 THEN BEGIN
	    phs = *lon.phs
            londata[j,*] = REFORM(phs[latindices[j],*])
            IF lon.zonalfreq EQ 4 THEN londata[j,*] = ALOG10(londata[j,*])
          ENDIF       
        endfor 

      ENDIF ELSE BEGIN
    
        iwvn = lon.wavenumber
        ifreq = lon.frequency
        ntime = field.ntime
        ntime2 = ntime/2+1
        nshift = (nlon+1)/2-1
        lonshift = lons[0] / (lons[1] - lons[0])
;
;  Need to pick out correct westward(negative)/eastward(positive) wavenumber since 0 is at center of array
;  and set frequency index to positive if specified as negative
;      
        iwvn = nshift - iwvn
;        IF ifreq LT 0 THEN iwvn = nshift + iwvn
;        IF ifreq GE 0 THEN iwvn = nshift - iwvn
;        IF ifreq LT 0 THEN ifreq = -ifreq
;
;  Do 2D FFT to get amplitude and phase
;       
        uspc = fltarr(nlon,ntime2,field.nlat,nlev)
        uphs = fltarr(nlon,ntime2,field.nlat,nlev)
        for k = 0,nlev-1 do begin
;
;  Handle user specified vertical range
;            
          for j = 0,nlat-1 do begin
            tmpa = fltarr(nlon,ntime)
            for it = 0,ntime-1 do begin
              for i = 0,nlon-1 do begin
                tmpa(i,it)=data(i,j,k,it)
              endfor
              tmpa(*,it) = shift(tmpa(*,it),lonshift) ;;this is to put 0 deg at the first point
            endfor
            fftmp = fft(tmpa)

            for it = 0,ntime2-1 do begin
              uspc(*,it,j,k) = shift(abs(fftmp(*,it)),nshift)*2.
              uphs(*,it,j,k) = -shift(atan(fftmp(*,it),/phase),nshift)*12./!pi
            endfor
          endfor
  
        endfor

        uphs(where(uphs lt 0))=uphs(where(uphs lt 0))+24.

        uspc(*,0,*,*)=uspc(*,0,*,*)/2.

        uphs(where(uphs le 0)) = min(uphs(where(uphs gt 0)))
;
;  Fill output data array
;                 
        IF lon.zonalfreq EQ 1 OR lon.zonalfreq EQ 2 THEN BEGIN
          londata = REFORM(uspc[iwvn,ifreq,*,*])
	  IF lon.zonalfreq EQ 2 THEN londata[*,*] = ALOG10(londata[*,*])
        ENDIF ELSE IF lon.zonalfreq EQ 3 OR lon.zonalfreq EQ 4 THEN BEGIN
          londata = REFORM(uphs[iwvn,ifreq,*,*])
	  IF lon.zonalfreq EQ 4 THEN londata[*,*] = ALOG10(londata[*,*])
        ENDIF	
;
;  Save amplitude and phase so don't have to recalculate and set flag
;   
        lon.amp = ptr_new(REFORM(uspc[iwvn,ifreq,*,*]))
        lon.phs = ptr_new(REFORM(uphs[iwvn,ifreq,*,*]))
 
        lon.ampphs = 1 ;Value of 1 means 2D FFT already calculated
 
      ENDELSE; Previously calculated 2D FFT for current field
;
;  Case where just need to load data into output array
;
    ENDIF ELSE BEGIN ;No spectral analysis

      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
;
; 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') AND (lon.amplitude EQ 0) then begin
        missing_value = float(lon.missing_value)
        log10f,londata,missing_value
      endif

    ENDELSE ;spectral analysis or not

  endelse ;vertical level range set or not

ENDELSE ;WACCM or TGCM
;
; Update lon structure:
;
lon.data = ptr_new(londata)
lon.lats = ptr_new(lats)
lon.levs = ptr_new(levs)
;
;  Interpolate WACCM data to 'standard' log pressure levels
;
IF lon.ftype EQ 'WACCM' THEN BEGIN

  if lon.plotz eq 0 then interpwdata, lon

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 (lon.plotz GE 0 AND lon.levmin EQ 0 AND lon.levmax EQ 0 AND lon.ftype NE 'WACCM') THEN *lon.levszp = *lon.levs

end