PRO getzlatlonut, info, slice, slicetype, fields, zheights
;
;  This routine reads the height field from the model history input file and forms a 
;  vertical altitude grid for longitude and latitude slices and interpolates the data
;  to the new vertical grid 
;   
; Get heights from Z (or Z3 from WACCM) field:
;
  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 or Z3 field -- extra right-hand',$
      ' axis NOT added.'
    return
  endif
;
; Read Z field if necessary:
; 12/8/05 btf: Currently, altyax is not called for mag fields (see rhyaxis
;              in pltxxx routines). However, ZMAG should be available on
;              histories now, so this routine should check if current field
;              is mag, and if it is, use ZMAG instead of Z.
;
  if not ptr_valid(fields[ixz].data) then begin ; field has not been read
    print,'Reading Z field to get heights for lat/lon..'
    widget_control,/hourglass
;
; Read Z data from file:
;
    varget,info,fields[ixz],zdata
    field = fields[ixz]
;
; Process Z field:
;
    procfield,info,zdata,field,info.z_hist
    fields[ixz] = field
;
;  Update and report minimum and maximum Z field value
;    
    fmin = min(*fields[ixz].data) & fmax = max(*fields[ixz].data)
    rpt_minmax,info,fields[ixz],fmin,fmax
    *info.fields = fields
  endif
;
;  Get the data to use and save any variables  
;
  fielddata = *fields[ixz].data
;
;  Need to modify altitude scale for fields at interface levels
;    
  fieldixz = fields[ixz]
  nLevsZP = slice.field.nlev
  nLevsIn = fieldixz.nlev
  fileInfo = *info.pfile
  iLevs = fileInfo.ilevs
  nILevs = N_ELEMENTS(iLevs)
;
;  Handle interface variables in WACCM
;  
  IF (info.ftype EQ 'WACCM' AND slice.field.levstype EQ 'interface') THEN BEGIN 
  
    inAlts = fielddata

    altsIZPlus1 = inAlts[0:fieldixz.nlon-1,0:fieldixz.nlat-1,1:nLevsIn-1,0:fieldixz.ntime-1]
    altsDiffs = altsIZPlus1 - inAlts[*,*,0:nLevsIn-2,*]
    altsInterface = altsDiffs/2.
    altsILevs = FLTARR(fieldixz.nlon,fieldixz.nlat,nILevs,fieldixz.ntime)
    altsILevs[*,*,0:nILevs-3,*] = inAlts[*,*,0:nLevsIn-2,*] - altsInterface
    altsILevs[*,*,nILevs-2,*] = inAlts[*,*,nLevsIn-1,*] - altsInterface[*,*,nLevsIn-2,*]
    altsILevs[*,*,nILevs-1,*] = inAlts[*,*,nLevsIn-1,*] + altsInterface[*,*,nLevsIn-2,*]
 
    fielddata = altsILevs
    
    print, 'Altitudes converted from midpoint to interface levels '

  ENDIF
;
;  Get current data and dimensions from z field
;
  slicedata = *slice.data

  dims = SIZE(fielddata)
  field = slice.field
  lons = *field.lons
  lats = *field.lats
  mtimes = slice.mtime
;
;  Two plot types handled here, longitude and latitude slices
;
  IF slicetype EQ 'lon' THEN BEGIN
  
    dims = SIZE(slicedata)
    nZH = dims[2]
    nLats = dims[1]
    zheights2D = FLTARR(nLats,nZH)
    zheights = FLTARR(nZH)

    ilon = ixfind_nearest(lons,slice.slon)
;
;  Handle if user set minimum and maximum x axis latitude
;
    iBegLat = 0
    iEndLat = nLats - 1
    nLatsF = N_ELEMENTS(lats)
    IF nLatsF NE nLats THEN BEGIN
      latsS = *slice.lats
      iBegLat = ixfind_nearest(lats,latsS[0])
      iEndLat = ixfind_nearest(lats,latsS[N_ELEMENTS(latsS)-1])
    ENDIF
;
;  Get z field 
;
    zheights2D[0:nLats-1,*] = REFORM(fielddata[ilon,iBegLat:iEndLat,*,slice.imtime])
;
;  WACCM handled differently than TGCM.  Vertical grid is the same size as input data not zp grid. 
;  zp grid in WACCM is not from input file but formed elsewhere?
;
    IF info.ftype EQ 'WACCM' THEN BEGIN 
    
      minzheights2D = MIN(zheights2D(*,0))
      maxzheights2D = MAX(zheights2D(*,nZH-1))
      delzheights2D = (maxzheights2D - minzheights2D) / nZH
      zheightsGrid = FLTARR(nZH)
      FOR izhG = 0, nZH-1 DO zheightsGrid[izhG] = minzheights2D + izhG * delzheights2D
      zheightsGrid = zheightsGrid
      nLevsZGrid = nZH

    ENDIF ELSE BEGIN
           
      nLevsZGrid = FIX(ABS(FIX(MAX(zheights2D(*,nZH-1))) - FIX(MIN(zheights2D(*,0)))) / 2.)
      zheightsGrid = FINDGEN(nLevsZGrid) * 2. + FIX(MIN(zheights2D(*,0)))
  
    ENDELSE
;
;  Interpolate data to new z grid
;  
    newData = FLTARR(nLats,nLevsZGrid)

    FOR iLat = 0, nLats-1 DO BEGIN
  
      newData(iLat,*) = interpol(slicedata(iLat,*), zheights2D(iLat,*), zheightsGrid)
  
    ENDFOR
;
;  Get vertical z grid for latitude slices
;    
  ENDIF ELSE IF slicetype EQ 'lat' THEN BEGIN 
  
    dims = SIZE(slicedata)
    nZH = dims[2]
    nLons = dims[1]
    zheights2D = FLTARR(nLons,nZH)
    zheights = FLTARR(nZH)

    ilat = ixfind_nearest(lats,slice.slat)
;
;  Handle if user set minimum and maximum x axis latitude
;
    iBegLon = 0
    iEndLon = nLons - 1
    nLonsF = N_ELEMENTS(lons)
    IF nLonsF NE nLons THEN BEGIN
      lonsS = *slice.lons
      iBegLon = ixfind_nearest(lons,lonsS[0])
      iEndLon = ixfind_nearest(lons,lonsS[N_ELEMENTS(lonsS)-1])
    ENDIF
 
    zheights2D[0:nLons-1,*] = REFORM(fielddata[iBegLon:iEndLon,ilat,*,slice.imtime])
;
;  WACCM handled differently than TGCM.  Vertical grid is the same size as input data not zp grid. 
;  zp grid in WACCM is not from input file but formed elsewhere?
;
    IF info.ftype EQ 'WACCM' THEN BEGIN 
    
      minzheights2D = MIN(zheights2D(*,0))
      maxzheights2D = MAX(zheights2D(*,nZH-1))
      delzheights2D = (maxzheights2D - minzheights2D) / nZH
      zheightsGrid = FLTARR(nZH)
      FOR izhG = 0, nZH-1 DO zheightsGrid[izhG] = minzheights2D + izhG * delzheights2D
      zheightsGrid = zheightsGrid
      nLevsZGrid = nZH

    ENDIF ELSE BEGIN
           
      nLevsZGrid = FIX(ABS(FIX(MAX(zheights2D(*,nZH-1))) - FIX(MIN(zheights2D(*,0)))) / 2.)
      zheightsGrid = FINDGEN(nLevsZGrid) * 2. + FIX(MIN(zheights2D(*,0)))
  
    ENDELSE
;
;  Interpolate data to new z grid
;      
    newData = FLTARR(nLons,nLevsZGrid)

    FOR iLon = 0, nLons-1 DO BEGIN
  
      newData(iLon,*) = interpol(slicedata(iLon,*), zheights2D(iLon,*), zheightsGrid)
  
    ENDFOR    
;
;  Get vertical z grid for latitude slices
;    
  ENDIF ELSE IF slicetype EQ 'ut' THEN BEGIN 
  
    dims = SIZE(slicedata)
    nZH = dims[2]
    nTimes = dims[1]
    zheights2D = FLTARR(nTimes,nZH)
    zheights = FLTARR(nZH)

    ilat = ixfind_nearest(lats,slice.slat)
    ilon = ixfind_nearest(lons,slice.slon)
;
;  Handle if user set minimum and maximum x axis latitude
;
    iBegTime = 0
    iEndTime = nTimes - 1
;    nMTimesF = N_ELEMENTS(mtimes[0,*])
;    IF nMTimesF NE nTimes THEN BEGIN
;      mTimesS0 = mtimes[0,0]*24.*60. + mtimes[1,0]*60. + mtimes[2,0]
;      iBegTime = ixfind_nearest(lons,m)
;      iEndTime = ixfind_nearest(lons,lonsS[N_ELEMENTS(lonsS)-1])
;    ENDIF
 
    for iTime = 0, nTimes-1 do zheights2D[iTime,*] = REFORM(fielddata[ilon,ilat,*,iTime])
;
;  WACCM handled differently than TGCM.  Vertical grid is the same size as input data not zp grid. 
;  zp grid in WACCM is not from input file but formed elsewhere?
;
    IF info.ftype EQ 'WACCM' THEN BEGIN 
    
      minzheights2D = MIN(zheights2D(*,0))
      maxzheights2D = MAX(zheights2D(*,nZH-1))
      delzheights2D = (maxzheights2D - minzheights2D) / nZH
      zheightsGrid = FLTARR(nZH)
      FOR izhG = 0, nZH-1 DO zheightsGrid[izhG] = minzheights2D + izhG * delzheights2D
      zheightsGrid = zheightsGrid
      nLevsZGrid = nZH

    ENDIF ELSE BEGIN
           
      nLevsZGrid = FIX(ABS(FIX(MAX(zheights2D(*,nZH-1))) - FIX(MIN(zheights2D(*,0)))) / 2.)
      zheightsGrid = FINDGEN(nLevsZGrid) * 2. + FIX(MIN(zheights2D(*,0)))
  
    ENDELSE
;
;  Interpolate data to new z grid
;      
    newData = FLTARR(nTimes,nLevsZGrid)

    FOR iTime = 0, nTimes-1 DO BEGIN
  
      newData(iTime,*) = interpol(slicedata(iTime,*), zheights2D(iTime,*), zheightsGrid)
  
    ENDFOR    

  ENDIF ELSE BEGIN
  
    print, 'Unsupported plot type in getzlatlonut - Returning '
    return
    
  ENDELSE
;
;  Convert geopotential to geometric height if requested   
;   
  IF slice.plotz EQ 2 THEN BEGIN

    zheightsGridTemp = zheightsGrid
    zheightsGrid = zheightsGridTemp * (1. + zheightsGridTemp/6370.0)
 
  ENDIF
;
;  Handle if user set minimum and maximum y axis vertical levels
; 
  IF slice.levminz GT zheightsGrid[0] OR $
        (slice.levmaxz LT zheightsGrid[nLevsZGrid-1] AND slice.levmaxz GT 0) THEN BEGIN
    iBegLev = ixfind_nearest(zheightsGrid,slice.levminz)
    iEndLev = ixfind_nearest(zheightsGrid,slice.levmaxz)
    nLevsZG = iEndLev - iBegLev + 1
    newDataTemp = newData
    IF slicetype EQ 'lon' THEN newData = FLTARR(nLats,nLevsZG)
    IF slicetype EQ 'lat' THEN newData = FLTARR(nLons,nLevsZG)
    IF slicetype EQ 'ut' THEN newData = FLTARR(nTimes,nLevsZG)
    newData[*,0:nLevsZG-1] = newDataTemp[*,iBegLev:iEndLev]
    zheightsGridTemp = zheightsGrid
    zheightsGrid = FLTARR(NLevsZG)
    zheightsGrid = zheightsGridTemp[iBegLev:iEndLev]
  ENDIF
 
  *slice.data = newData
  zheights = zheightsGrid
;    
;  Put data in slice structure
;    
  *slice.levs = zheights
  slice.levminz = MIN(zheights)
  slice.levmaxz = MAX(zheights)
  slice.fmin = MIN(newData)
  slice.fmax = MAX(newData)

  RETURN

END