;------------------------------------------------------------------
pro defutvertdata,utvert,data
;this procedure reads in the data that is needed for the utvert plot
;(mtime,lev,lat,lon), it also computes the zonal and global mean, compared 
;with fortran tgcm code there is a difference to the tenths place on the 
;min and max of the data for both means
;written by Tom Freestone 9/3/03

;
; Define utvert.data(ntime,nlev) from processed var in data.
;
;
;  If plotting with altitude vertical axis, make sure data and levels are on zp since
;  levels on altitude won't work in deflondata
;
IF utvert.plotz GE 1 THEN BEGIN

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

field = utvert.field
ntime = field.ntime & nlev = field.nlev
;
; Check dimensions (return if field is only 2d lat x utvert)
;
dims = size(data,/dimensions)
ndims = n_elements(dims)
checkdims=size(dims)
error = 0
if (dims[0] ne field.nlon) then begin
  print,'>>> WARNING defutvertdata: 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 defutvertdata: 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 defutvertdata: 3rd dimension is not nlevs: dims[2]=',$
    dims[2],' but nlev=',field.nlev
  error = 1 
endif
if checkdims[1] gt 3 then begin
if (dims[3] ne field.ntime) then begin
  print,'>>> WARNING defutvertdata: 4th dimension is not ntime: dims[3]=',$
    dims[3],' but ntime=',field.ntime
  error = 1 
endif
endif
if nlev le 0 then begin
  print,'>>> defutvertdata: nlev=',nlev,' -- cannot make vertical slice.'
  error = 1 
endif
if error gt 0 then begin
  size = size(utvert.data)
  utvertdata = fltarr(size[1],size[2])
  utvertdata[*,*] = float(utvert.missing_value)
  utvert.data = ptr_new(utvertdata)
  return
endif
;
lons = *field.lons
lats = *field.lats
file = utvert.file
mtime = file.mtimes
nlon = field.nlon
utmin = utvert.utmin & utmax = utvert.utmax
interval = utvert.utint
dims=size(mtime)
for i=0,file.ntimes-1 do begin
  if file.mtimes[0,i] eq utmin[0] and file.mtimes[1,i] eq utmin[1] and $
    file.mtimes[2,i] eq utmin[2] then begin
    imtime_min=i
  endif
  if file.mtimes[0,i] eq utmax[0] and file.mtimes[1,i] eq utmax[1] and $
    file.mtimes[2,i] eq utmax[2] then begin
    imtime_max=i
  endif
endfor
imax=imtime_max
imin=imtime_min
if dims[0] eq 0 then begin
  nut =1
  times=float(mtime[0]*24.*60.+mtime[1]*60.+mtime[2])
endif else begin 
  nut =ceil(float(dims[2])/float(interval))
  times=findgen(nut)
  for x=0,nut-1 do begin
    times[x]=float(mtime[0,x]*24.*60.+mtime[1,x]*60.+mtime[2,x])
  endfor
endelse
timesindices = indgen(ntime)
fieldlevs = *field.levs
nfieldlevs = n_elements(fieldlevs)
;
; If axis limits levmin,max are set, reset levs accordingly:
;
if utvert.levmin lt utvert.levmax then begin ; set axes limits levs:
  nlev = 0
  for k=0,nfieldlevs-1 do $
    if fieldlevs[k] ge utvert.levmin-0.001 and fieldlevs[k] le utvert.levmax+0.001 then nlev = nlev+1

  if nlev lt 2 then print,'>>> WARNING pltutvert: 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 utvert.levmin-0.001 and fieldlevs[k] le utvert.levmax+0.001 then begin
      levs[kk] = fieldlevs[k]
      kk = kk+1
    endif
endif else levs = fieldlevs
nlev = n_elements(levs)
;
;determines the time array that the data should be plotted against
;
uts=mtime
utmin = utvert.utmin & utmax = utvert.utmax
utmin_mins = utmin[0]*60.*24. + utmin[1]*60. + utmin[2]
utmax_mins = utmax[0]*60.*24. + utmax[1]*60. + utmax[2]
itimesmin=ixfind_nearest(times,utmin_mins)
if utmin_mins gt times[0] or utmax_mins lt times[nut-1] then begin
  nutsnew = 0
  for j=0,nut-1 do begin
    if times[j] ge utmin_mins and times[j] le utmax_mins then begin
      nutsnew = nutsnew+1
    endif
  endfor
  if nutsnew le 1 then begin
    print,'>>> WARNING defutvertdata: utmin,max=',utmin,utmax,' nutsnew=', $
      nutsnew
    print,'    Cannot plot this utmin,max -- resorting to full UT ',$
      'grid..'
    nutsnew = nut
  endif else begin
    newuts = intarr(3,nutsnew)
    utindnew = intarr(nutsnew)
    jj = 0
    for j=0,nut-1 do begin
      if times[j] ge utmin_mins and times[j] le utmax_mins then begin
        newuts[*,jj] = mtime[*,j]
        utindnew[jj] = j
        jj=jj+1
      endif
    endfor
    uts = newuts
    nut = nutsnew
    utindices = utindnew
endelse
endif
;
; Define (nut,nlev) at utvert.slon and utvert.slat:
;
utvertdata = fltarr(nut,nlev)
utvertdata[*,*] = 0.

ilat = ixfind_nearest(lats,utvert.slat)
ilon = ixfind_nearest(lons,utvert.slon)
utvert.slat = lats[ilat]
utvert.slon = lons[ilon]
nlat = n_elements(lats)
nlon = n_elements(lons)

;
if utvert.levmin lt utvert.levmax then begin 
  if utvert.globalmean eq 1 then print,'Starting Global Mean, Please Wait...'
  if utvert.globalmean eq 1 then widget_control,/hourglass
  kk = 0
  for k=0,nfieldlevs-1 do $
    if fieldlevs[k] ge utvert.levmin-0.001 and fieldlevs[k] le utvert.levmax+0.001 then begin
        if utvert.zonalmean eq 0 and utvert.globalmean eq 0 then begin
          for j=0,nut-1 do begin
            utvertdata[j,kk] = float(data[ilon,ilat,k,j*interval+itimesmin])
          endfor
        endif else begin ; plot zonal mean
	  for j=0,nut-1 do begin
	    if utvert.zonalmean eq 1 then begin
              utvertdata[j,kk] = float(mean(data[*,ilat,k,j*interval+itimesmin]))
            endif
            if utvert.globalmean eq 1 then $
              utvertdata[j,kk] = glbmean(data[*,*,k,j*interval+itimesmin], $
		nlon,nlat,lats,utvert.missing_value)
          endfor
        endelse
      kk = kk+1
    endif
endif else begin
    if utvert.zonalmean eq 0 and utvert.globalmean eq 0 then begin
      for j=0,nut-1 do begin
        utvertdata[j,*] = float(data[ilon,ilat,*,j*interval+itimesmin])
      endfor
    endif else begin ; zonal mean
      for j=0,nut-1 do begin
        for k=0,nlev-1 do begin
	  if utvert.zonalmean eq 1 then $
            utvertdata[j,k] = float(mean(data[*,ilat,k,j*interval+itimesmin])) 
          if utvert.globalmean eq 1 then $
            utvertdata[j,k] = glbmean(data[*,*,k,j*interval+itimesmin],nlon, $
	      nlat,lats,utvert.missing_value)
        endfor
      endfor
    endelse ; zonal mean
endelse
if utvert.globalmean eq 1 then print,'Finished Global Mean'

;
; Take log10 of densities:
;
if (utvert.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
   (utvert.log10 eq 'this frame only') or $
   (utvert.log10 eq 'all fields') then begin
  missing_value = float(utvert.missing_value)
  log10f,utvertdata,missing_value
endif
;
; Update utvert structure:
;
*utvert.data = utvertdata
*utvert.levs = levs
*utvert.uts = uts
;
;  Interpolate WACCM data to 'standard' log pressure levels
;
IF utvert.ftype EQ 'WACCM' THEN BEGIN

  if utvert.plotz eq 0 then interpwdata, utvert

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