;
pro getmapzdata,info,map,data 
;
;  This routine takes the requested altitude for mapping and interpolates the data to that 
;  altitude
; 
;  Get z coordinate data index
;
  field = map.field
  ixtime = map.imtime
  mapalt = map.alt

  fields = *info.fields

  ixz = -1
  for i=0,n_elements(fields)-1 do 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
  endfor
;  if ixz eq -1 then begin
;    print,'>>> WARNING altyax: Could not find Z field for mapping'
;    return
;  endif
;
; 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 map.plotz EQ 2 THEN BEGIN

    zdatatemp = zdata
    zdata = zdatatemp * (1. + zdatatemp/6370.0)
 
  ENDIF
 
  dims = size(zdata,/dimensions)
  ndims = n_elements(dims)
  
  nLon = dims[0]
  nLat = dims[1]

  mapdata = fltarr(nLon,nLat)
  
  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 (mapalt)
;
      iLev = IXFIND_NEAREST(zdata1timemean, mapalt)
    
      IF mapalt GT zdata1timemean[iLev] THEN BEGIN
  
        iLevLower = iLev
        iLevUpper = iLev + 1
    
      ENDIF ELSE IF mapalt LT zdata1timemean[iLev] THEN BEGIN
  
        iLevLower = iLev -1
        iLevUpper = iLev
    
      ENDIF

      IF mapalt NE zdata1timemean[iLev] THEN BEGIN
  
        zdata1tmL = zdata1timemean[iLevLower]
        zdata1tmU = zdata1timemean[iLevUpper]

        levDiffLower = mapalt - zdata1tmL
        levDiffUpper = zdata1tmU - mapalt
        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)

        FOR iLn = 0, nLon-1 DO BEGIN
	
	  FOR iLt = 0, nLat-1 DO BEGIN

            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(map.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
 	
  	      mapdata[iLn,iLt] = float(dataLower[iLn,iLt])*levRatUpper + float(dataUpper[iLn,iLt])*levRatLower

              IF isdensity(map.field.name) THEN mapdata[iLn,iLt] = 10^(mapdata[iLn,iLt])
	
	    ENDIF ELSE BEGIN; One of both data values missing so set output at this lat/lon to missing
	    
	      mapdata[iLn,iLt] = float(map.missing_value)
	      
	    ENDELSE
	    
	  ENDFOR
	
	ENDFOR
 	
      ENDIF ELSE BEGIN ;If chosen altitude equal to mean altitude of level
      
        mapdata = data[*,*,iLev,ixtime]
	
      ENDELSE		

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

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

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

        mapdata = float(data[*,*])
        map.plotz = 0 
      
    end
    else: begin
        print,'>>> defmapdata: var data has ',ndims,' dimensions --',$
  	' cannot make map data.'
        mapdata = fltarr(nlon,nlat)
        mapdata = 0.
        map.plotz = 0 
    end
  endcase

  map.data = ptr_new(mapdata)

return

end