PRO getzprof, info, slice, slicetype, fields, zheights
;
;  This routine reads the height field from the model history input file and picks out the 
;  subset of the data needed for line plots.  Needs to be cleaned up and checked for accuracy 
;  with meridional and zonal line plots. 
;
;
; 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 getzprof: 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 right-hand axes..'
    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)
  
  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

  slicedata = *slice.data	; save
  log10save = slice.log10	; save
  slicefmin = slice.fmin        ; save(changed in getzprof)
  slicefmax = slice.fmax        ; save(changed in getzprof)
;
;  Meridional line plot type needs data in latitude dimension and zonal line plot type needs data in 
;  longitude dimension and rest need data in vertical dimension
;
  IF slicetype eq 'meridional' THEN BEGIN

    dims = SIZE(fielddata)
    nZH = dims[3]
    zheights = FLTARR(nZH)

    FOR iZH = 0, nZH-1 DO BEGIN
 
      zheights[iZH] = FLOAT(MEAN(fielddata[slice.ilon,*,iZH,slice.imtime]))

    ENDFOR
  
  ENDIF ELSE IF slicetype eq 'zonal' THEN BEGIN

    dims = SIZE(fielddata)
    nZH = dims[3]
    zheights = FLTARR(nZH)

    FOR iZH = 0, nZH-1 DO BEGIN
 
      zheights[iZH] = FLOAT(MEAN(fielddata[*,slice.ilat,iZH,slice.imtime]))

    ENDFOR

  ENDIF ELSE IF slicetype eq 'lon' THEN BEGIN

    dims = SIZE(fielddata)
    nZH = dims[3]
    zheights = FLTARR(nZH)
    field = slice.field
    lons = *field.lons
    ilon = ixfind_nearest(lons,slice.slon)

    FOR iZH = 0, nZH-1 DO BEGIN
 
      zheights[iZH] = FLOAT(MEAN(fielddata[ilon,*,iZH,slice.imtime]))

    ENDFOR

  ENDIF ELSE IF slicetype eq 'lat' THEN BEGIN

    dims = SIZE(fielddata)
    nZH = dims[3]
    zheights = FLTARR(nZH)
    field = slice.field
    lats = *field.lats
    ilat = ixfind_nearest(lats,slice.slat)

    FOR iZH = 0, nZH-1 DO BEGIN
 
      zheights[iZH] = FLOAT(MEAN(fielddata[*,ilat,iZH,slice.imtime]))

    ENDFOR

  ENDIF ELSE IF slicetype eq 'ut' THEN BEGIN

    dims = SIZE(fielddata)
    nZH = dims[3]
    zheights = FLTARR(nZH)
    field = slice.field
    lats = *field.lats
    ilat = ixfind_nearest(lats,slice.slat)
    lons = *field.lons
    ilon = ixfind_nearest(lons,slice.slon)

    FOR iZH = 0, nZH-1 DO BEGIN
 
      zheights[iZH] = FLOAT(MEAN(fielddata[ilon,ilat,iZH,*]))

    ENDFOR

  ENDIF ELSE BEGIN

    dims = SIZE(fielddata)
    nZH = dims[3]
    zheights = FLTARR(nZH)
    zheights[*] = FLOAT(fielddata[slice.ilon,slice.ilat,*,slice.imtime])

  ENDELSE

  *slice.data = slicedata ; restore
 
  slice.log10 = log10save ; restore
  slice.fmin = slicefmin  ; restore(changed in getzprof)
  slice.fmax = slicefmax  ; restore(changed in getzprof)
;
;  Check to see if vertical coordinate is altitude and, if so, which option
;
  IF slicetype EQ 'profile' THEN BEGIN

    IF info.ftype EQ 'WACCM' AND slice.plotz GT 0 THEN BEGIN
    
      levszporig = *slice.field.levs
      levszpnew = *slice.levszp
      oldzheights = zheights
      zheights = interpol(oldzheights, levszporig, levszpnew)
    
    ENDIF    
;
;  Convert geopotential to geometric height if requested   
;   
    IF slice.plotz EQ 2 THEN BEGIN

      zheightsTemp = zheights
      zheights = zheightsTemp * (1. + zheightsTemp/6370.0)
 
    ENDIF
;    
;  Put data in slice structure
;    
    IF slice.plotz GE 1 THEN BEGIN
  
      *slice.levs = zheights
      slice.levmin = MIN(zheights)
      slice.levmax = MAX(zheights)

    ENDIF

  ENDIF

  RETURN

END