;
pro getaltzdata,info,slice,data 
;
;  This routine takes the requested altitude for mapping and interpolates the data to that 
;  altitude
; 
;  Get variables from map info structures
;
  field = slice.field
  ixtime = slice.imtime
  slicealt = slice.alt

  fields = *info.fields
; 
;  Get z coordinate data index
;
  ixz = -1
  for i=0,n_elements(fields)-1 do begin
;
;  So far 3 types of altitude fields so check according to map.plotz flag
;
    IF slice.plotz EQ 1 OR slice.plotz EQ 2 THEN BEGIN
     
      if strcompress(fields[i].name,/remove_all) eq 'Z' then ixz = i    
      if strcompress(fields[i].name,/remove_all) eq 'Z3' then ixz = i

    ENDIF ELSE IF slice.plotz EQ 3 THEN BEGIN
 
      if strcompress(fields[i].name,/remove_all) eq 'ZG' then ixz = i   

    ENDIF
  
  endfor
;
; Read Z field if necessary(not read before):
;
  if not ptr_valid(fields[ixz].data) then begin ; field has not been read
;    print,'Reading Z field for altitude mapping..'
    widget_control,/hourglass
;
; Read Z data from file:
;
    varget,info,fields[ixz],zdata
;
;  Adjust altitude to kilometers
;
    IF info.ftype EQ 'WACCM' THEN BEGIN ;Geopotential altitude is in meters in WACCM history files

      zdata = zdata * 1.e-3  ; m to km

    ENDIF ELSE BEGIN; Centimeters in TGCM history files 

      zdata = zdata * 1.e-5  ; cm to km

    ENDELSE
  
  ENDIF ELSE BEGIN ;Field already read so grab it
  
    zdata = *fields[ixz].data
    
  ENDELSE  
;
;  Convert geopotential to geometric height if requested   
;   
  IF slice.plotz EQ 2 THEN BEGIN

    zdatatemp = zdata
    zdata = zdatatemp * (1. + zdatatemp/6370.0)
 
  ENDIF
;
;  Get latitude/longitude dimensions from zdata and set up output array
; 
  dims = size(zdata,/dimensions)
  ndims = n_elements(dims)
  
  nLon = dims[0]
  nLat = dims[1]

  slicedata = fltarr(nLon,nLat)
;
;  Check for up to 4 dimensional zdata.  Only processes data for 4 dimensions and less than 4 cancels altitude
;  plotting and returns
;  
  case ndims of
    4: begin
      if dims[2] ne field.nlev then $
  	print,'>>> WARNING defmapdata: 3rd dimension of 4d data ',$
  	  'should be nlev=',field.nlev,' but is ',dims[2] 
      if dims[3] ne field.ntime and field.ntime ne 1 then $
  	print,'>>> WARNING defmapdata: 4th dimension of data ',$
  	  'should be ntime=',field.ntime,' but is ',dims[3] 
;
;  Pick out altitude for current time.  Assume dimensions lon,lat,lev,time
;
      zdata1time = zdata[*,*,*,ixtime]
    
      nzd = SIZE(zdata1time)

      nLev = nzd[3]
;    
;  Need to average over longitudes and latitudes to get one altitude profile and then choose
;  level to map
; 
      zdata1timemean = FLTARR(nLev)
         
      FOR iLev = 0, nLev-1 DO BEGIN
    
        zdata1timemean[iLev] = MEAN(zdata[*,*,iLev])
      
      ENDFOR              
;
;  Need to interpolate data to user requested altitude.
;  First find altitude values above and below user selected altitude (slicealt)
;
      iLev = IXFIND_NEAREST(zdata1timemean, slicealt)
    
      IF slicealt GT zdata1timemean[iLev] THEN BEGIN
  
        iLevLower = iLev
        iLevUpper = iLev + 1
    
      ENDIF ELSE IF slicealt LT zdata1timemean[iLev] THEN BEGIN
  
        iLevLower = iLev -1
        iLevUpper = iLev
    
      ENDIF
;
;  Do interpolation where slicealt is not equal to closest level found above
;
      IF slicealt NE zdata1timemean[iLev] THEN BEGIN
;
;  Find zdata values above and below user selected altitude, calculate weights for interpolation, and 
;  use only non missing data values
;  
        zdata1tmL = zdata1timemean[iLevLower]
        zdata1tmU = zdata1timemean[iLevUpper]

        levDiffLower = slicealt - zdata1tmL
        levDiffUpper = zdata1tmU - slicealt
        levDiffAll =  zdata1tmU - zdata1tmL
        levRatLower = levDiffLower / levDiffAll
        levRatUpper = levDiffUpper / levDiffAll   

        goodDataLower = WHERE(data[*,*,iLevLower,ixtime] NE field.missing_value, ngood)

        dataLower = fltarr(nLon,nLat)
        dataUpper = fltarr(nLon,nLat)
;
;  Loop to interpolate data vertically for each latitude/longitude
;
        FOR iLn = 0, nLon-1 DO BEGIN
	
	  FOR iLt = 0, nLat-1 DO BEGIN
;
;  Get surrounding values and check for missing.  If both valid then take log if isdensity
;  
            dataTempLower = data[iLn,iLt,iLevLower,ixtime]
            dataTempUpper = data[iLn,iLt,iLevUpper,ixtime]

  	    IF dataTempLower NE field.missing_value AND dataTempUpper NE field.missing_value THEN BEGIN

  	      IF isdensity(slice.field.name) THEN BEGIN

  	  	dataLower[iLn,iLt] = ALOG10(dataTempLower)
  	  	dataUpper[iLn,iLt] = ALOG10(dataTempUpper)
;  	  	 dataLower = ALOG10(data[*,*,iLevLower,ixtime])
;  	  	 dataUpper = ALOG10(data[*,*,iLevUpper,ixtime])
  	  
  	      ENDIF ELSE BEGIN

  	  	dataLower[iLn,iLt] = dataTempLower
  	  	dataUpper[iLn,iLt] = dataTempUpper
;  	  	 dataLower = data[*,*,iLevLower,ixtime]
;  	  	 dataUpper = data[*,*,iLevUpper,ixtime]

  	      ENDELSE
;
;  Interpolate and take inverse log if needed
; 	
  	      slicedata[iLn,iLt] = float(dataLower[iLn,iLt])*levRatUpper + float(dataUpper[iLn,iLt])*levRatLower

              IF isdensity(slice.field.name) THEN slicedata[iLn,iLt] = 10^(slicedata[iLn,iLt])
	
	    ENDIF ELSE BEGIN; One or both data values missing so set output at this lat/lon to missing
	    
	      slicedata[iLn,iLt] = float(slice.missing_value)
	      
	    ENDELSE
	    
	  ENDFOR
	
	ENDFOR
 	
      ENDIF ELSE BEGIN ;If chosen altitude equal to mean altitude of level
      
        slicedata = data[*,*,iLev,ixtime]
	
      ENDELSE
      
      IF slice.plottype EQ 'meridional' THEN BEGIN
      
        IF slice.zonalmean GT 0 THEN BEGIN  ; zonal mean

          FOR iLt = 0, nLat-1 DO BEGIN
    
            slicedata[0:nLon-1,iLt] =  float(mean(slicedata[*,iLt]))		
 
          ENDFOR
	
        ENDIF 
	
      ENDIF
      
      IF slice.plottype EQ 'zonal' THEN BEGIN
 
        IF slice.meridionalmean GT 0 THEN BEGIN  ; meridional mean

          FOR iLn = 0, nLon-1 DO BEGIN
    
            slicedata[iLn,0:nLat-1] =  float(mean(slicedata[iLn,*]))		
 
          ENDFOR
	
        ENDIF 
	
      ENDIF 
 
    end
;
;  For dimensions less than 4, assign pressure level data and quit altitude interpolation
;
    3: begin
        print,'>>> WARNING defmapzdata: No vertical dimension in input data. No interpolation done'

        slicedata = float(data[*,*,ixtime])
        slice.plotz = 0 
      
    end
    2: begin

        print,'>>> WARNING defmapzdata: No vertical dimension in input data. No interpolation done'

        slicedata = float(data[*,*])
        slice.plotz = 0 
      
    end
    else: begin
        print,'>>> getaltzdata: var data has ',ndims,' dimensions --',$
  	' cannot make slice data.'
        slicedata = fltarr(nlon,nlat)
        slicedata = 0.
        slice.plotz = 0 
    end
  endcase
;
;  Handle wave amplitude calculation here
;      
  IF slice.plottype EQ 'meridional' OR slice.plottype EQ 'zonal' THEN BEGIN
  
    IF slice.amplitude NE 0 AND field.name NE 'YAXZ' THEN BEGIN  ;Get wavenumber amplitude  
;
; Convert data that is log10 to linear:
;
      iwvn = slice.wavenumber
      ampdata = fltarr(nlon)
      slicelatdata = fltarr(nlat)
      slicewvndata = fltarr(nlon,nlat)
      ampmin = 1.0E+24
      for j=0,nlat-1 do begin
    
        ampdata[*] = slicedata[*,j]
        ampdata[*] = ABS(FFT(ampdata[*], -1)) * 2.
        if ampdata[iwvn] LT ampmin AND ampdata[iwvn] GT 0.0 then ampmin = ampdata[iwvn]
        slicelatdata[j] = ampdata[iwvn]
      
        IF slice.plottype EQ 'zonal' THEN BEGIN

          slicewvndata[*,j] = ampdata
        
	  zeroind = WHERE(slicewvndata[*,j] EQ 0.0)
          if zeroind[0] NE -1 then slicewvndata[zeroind,j] = ampmin
;        help, zeroind

          IF slice.amplitude EQ 2 THEN slicewvndata[*,j] = ALOG10(slicewvndata[*,j])
    
        ENDIF

      endfor
  
      IF slice.plottype EQ 'meridional' THEN BEGIN

        zeroind = WHERE(slicelatdata[*] EQ 0.0)
        if zeroind[0] NE -1 then slicelatdata[zeroind] = ampmin

        IF slice.amplitude EQ 2 THEN slicelatdata[*] = ALOG10(slicelatdata[*])
      
        FOR i=0,nlon-1 DO slicedata[i,*] = slicelatdata[*]

      ENDIF
    
      IF slice.plottype EQ 'zonal' THEN slicedata[*,*] = slicewvndata[*,*]
;
;  Update min/max for data
; 
      slice.fmin = MIN(slicedata)
      slice.fmax = MAX(slicedata) 

    ENDIF

  ENDIF ELSE BEGIN
;
; Take log10 of densities:
;
    IF slice.plottype EQ 'meridional' OR slice.plottype EQ 'zonal' THEN BEGIN
      if (slice.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
         (slice.log10 eq 'this frame only') or $
         (slice.log10 eq 'all fields') then begin
        missing_value = float(slice.missing_value)
        log10f,slicedata,missing_value
      endif    
  
    ENDIF

  ENDELSE ; wavenumber calculation or not
;
;  Update data in slice structure
;
  slice.data = ptr_new(slicedata)

return

end