;
; User input (called from main pro utproc below):
;
pro input,dirs,fileroots,file_beg,file_end,form,fields,loc_lons,loc_lats,$
  zpbot,zptop,mtime_beg,mtime_end,image_dir
;
; timegcm1.2_y2c:
;
;dirs = ["/aim/d/tgcm/timegcm1.2_y2c/sy2c"]
;fileroots = ["ROBLE.timegcm1.2.sy2c"]
;file_beg = [1]  & file_end = [30] ; 1,0,0 to 59,23,0 by 60
;file_beg = [1]  & file_end = [1]
;form = "(i03)"
;mtime_beg = [1,0,0]  & mtime_end  = [7,0,0]
;mtime_beg = [-1,0,0]  & mtime_end  = [-1,0,0]
;fields = ["TN","UN","VN","W","O2","O1","NO","N4S","NE","TE","TI","O2P","Z"]
;
;dirs = ["/aim/d/tgcm/tiegcm1.81_yr02/py2s"]
;fileroots = ["FOSTER.tiegcm1.81.py2s"]
 dirs = ["/aim/d/tgcm/tiegcm1.81_yr02/sy2s"]
 fileroots = ["FOSTER.tiegcm1.81.sy2s"]
 file_beg = [1]  & file_end = [1]
 mtime_beg = [-1,0,0]  & mtime_end  = [-1,0,0]
;file_beg = [1]  & file_end = [16]
;mtime_beg = [1,0,0]  & mtime_end  = [21,0,0]
;file_beg = [42]  & file_end = [52]
;mtime_beg = [70,0,0]  & mtime_end  = [90,0,0]
;file_beg = [93]  & file_end = [103]
;mtime_beg = [170,0,0]  & mtime_end  = [190,0,0]
;file_beg = [143]  & file_end = [153]
;mtime_beg = [270,0,0]  & mtime_end  = [290,0,0]
 fields = ["TN","UN","VN","OMEGA","O2","O1","NO","N4S","NE","TE","TI","O2P","Z"]
 form = "(i03)"
;
;
 loc_lons = [0.,0.,0.,0.,0.,0.,0.,0.,0.]
 loc_lats = [  -80. , -70., -60., -50., 0., 50.,  60.,  70.,  80.]
;loc_lons = [0.] & loc_lats = [60.]
;
;zpbot = -17. & zptop = 5.
;zpbot = -17. & zptop = -7.
 zpbot = -999. & zptop = 999.
;zpbot = -999. & zptop = 5.
;
; Selected lons and zp levels for ut vs lat:
;
 slons = [-180.,-120.,-60.,0.,60.,120.,180.]
;slevs = [-17., -10., -4., 0., 2., 4.]
 slevs = [-7., -4., 0., 2., 4.]
;
;image_dir = '.'
;image_dir = '/oreo/d/foster/utproc'
 image_dir = '/toshi/ftp/pub/foster/333roxbury/tgcm/timegcm1.2_sy2c/images'
end ; pro input
;-----------------------------------------------------------------------
function ixfind, z,val,del
;
; Find nearest index to val in z, with increment del
; If not found, return -1
;
nz = n_elements(z)
if del gt 0 then begin
  for i=0,nz-1 do begin
    if (val ge z(i)-0.5*del and val le z(i)+0.5*del) then return,i
  endfor
endif else begin
  for i=0,nz-1 do begin
    if (val eq z(i)) then return,i
  endfor
endelse
return,-1
end
;-----------------------------------------------------------------------
pro dayofyear_to_monthday, dayofyear,imonth,iday
;
; Given day of year (1-365), return month and day
;
iday = 0
imonth = 0
if (dayofyear eq 0) then dayofyear = 1 ; to allow for model day of 0
if (dayofyear lt 1 or dayofyear gt 365) then begin
  print,'>>> bad dayofyear=',dayofyear
  return
endif
;                J, F, M, A, M, J, J, A, S, O, N, D
dayspermonth = [31,28,31,30,31,30,31,31,30,31,30,31]
nmonths = n_elements(dayspermonth)
;
id = 0
for m=0,nmonths-1 do begin
  for i=1,dayspermonth[m] do begin
    id = id+1
    if (id eq dayofyear) then begin
      imonth = m+1
      iday = i
      break  
    endif
  endfor
  if (imonth ne 0) then break ; day was found above
endfor
;print,'dayofyear_to_monthday: dayofyear=',dayofyear,' imonth=',imonth,' iday=',iday
end
;-----------------------------------------------------------------------
function monthday_to_dayofyear, imonth,iday
;                J, F, M, A, M, J, J, A, S, O, N, D
dayspermonth = [31,28,31,30,31,30,31,31,30,31,30,31]
day = 0
for i=0,11 do begin
  if (i+1 eq imonth) then begin
    day = day+iday
    break
  endif
  day = day+dayspermonth[i]
endfor
;print,'monthday_to_dayofyear: imonth=',imonth,' iday=',iday,' dayofyear=',day
return,day
end
;-----------------------------------------------------------------------
function mtime_to_mins,mtime ; model time to integer minutes
return,mtime(0)*24*60+mtime(1)*60+mtime(2)
end
;-----------------------------------------------------------------------
function mtime_to_ddays,mtime ; model time to decimal days
return,float(mtime(0))+float(mtime(1))/24.+float(mtime(2))/(24.*60.)
end
;-----------------------------------------------------------------------
function timeaxis,years,mtimes,deltamin,xtitle,nxticks
;
; Get julian day of first time:
;
dayofyear_to_monthday,mtimes(0,0),imonth0,iday0
julian0 = julday(imonth0,iday0,years[0],mtimes(1,0),mtimes(2,0),0)
;
; Generate x-axis coords in julian days (step size in minutes):
; (this array is returned as the function value)
;
ntimes = n_elements(years)
time = timegen(ntimes,start=julian0,step_size=deltamin,units="Minutes")
;
; Add 0.5 to times since the julian day starts at 12 pm noon 
; (s.a., caldat call below)
;
time = time+0.5 
;
; Construct x-axis label giving beginning,ending model times,
; delta time in miniutes, and total time in days:
;
ddays0 = mtime_to_ddays(mtimes[*,0])
ddays1 = mtime_to_ddays(mtimes[*,ntimes-1])
total_ddays = ddays1-ddays0 ; total number of decimal days 
nxticks = 0 ; let idl select
;if (total_ddays ge 364) then nxticks = 7 ; does not appear to work
xtitle = string(format=$
  "('Model time (yyyy:ddd:hh:mm) ',i4,':',i3.3,':',i2.2,':',i2.2,' to ')",$
  years[0],mtimes[*,0])
xtitle = xtitle + string(format=$
  "(i4,':',i3.3,':',i2.2,':',i2.2,' by ',i6,' minutes ')",$
  years[ntimes-1],mtimes[*,ntimes-1],deltamin)
xtitle = xtitle + string(format="('(',f7.2,' days)')",total_ddays)
return,time
end
;-----------------------------------------------------------------------
; Main:
;
pro utproc,ps=ps,jpg=jpg
input,dirs,fileroots,file_beg,file_end,form,fields,loc_lons,loc_lats,$
  zpbot,zptop,mtime_beg,mtime_end,image_dir
ndirs = n_elements(dirs)
nfiles = 0
for idir=0,ndirs-1 do begin
  nfiles_dir = file_end[idir]-file_beg[idir]+1
  nfiles = nfiles+nfiles_dir
endfor
files = strarr(nfiles)
ifiles = 0
for idir=0,ndirs-1 do begin
  nfiles_dir = file_end[idir]-file_beg[idir]+1
  for ifile=0,nfiles_dir-1 do begin
    i = file_beg[idir]+ifile
    strnum = string(format=form,i)
    if (idir eq 0 and ifile eq 0) then file0 = fileroots[idir]+strnum+".nc"
    if (idir eq ndirs-1 and ifile eq nfiles_dir-1) then $
      file1 = fileroots[idir]+strnum+".nc"
    files[ifiles] = dirs[idir] + "/" + fileroots[idir] + strnum + ".nc"
    ifiles = ifiles+1 
  endfor
endfor
;
nfields = n_elements(fields)
f_units = strarr(nfields)
ncids = intarr(nfiles)
if (n_elements(loc_lons) ne n_elements(loc_lats)) then begin
  print,'>>> Number of loc_lons must equal number of loc_lats:'
  print,'    nloc_lons=',n_elements(loc_lons),' nloc_lats=',$
    n_elements(loc_lats)
  return
endif
nlocs = n_elements(loc_lons)
print,format="('nfiles=',i4,' nfields=',i4,' nlocs=',i4)",$
  nfiles,nfields,nlocs
;
mtime0 = mtime_beg
mtime1 = mtime_end
;
; Search for start and end times:
;
ifile0 = -1 & ifile1 = -1
ihist0 = -1 & ihist1 = -1
for ifile=0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  id = ncdf_varid(ncids[ifile],'mtime')
  if (id eq -1) then begin
    print,'>>> could not find id of mtime var on file ',files[ifile]
    return
  endif
  ncdf_varget,ncids[ifile],id,mtime
  ntime = n_elements(mtime)/3
  for i=0,ntime-1 do begin
;   print,format=$
;     "('Searching for mtime0=',3i4,' or mtime1=',3i4,' Found ',3i4,' on file ',a)",$
;     mtime0,mtime1,mtime[*,i],file_basename(files[ifile])
;
; If start day is < 0, use first history on first file:
    if (mtime0[0] lt 0 and ifile eq 0 and i eq 0) then begin
      ifile0 = ifile
      ihist0 = i
      print,format="('Using start mtime ',3i4,' (first history of first file ',a,')')",$
        mtime[*,ihist0],file_basename(files[ifile0])
    endif else begin
      if (mtime[0,i] eq mtime0[0] and mtime[1,i] eq mtime0[1] and $
          mtime[2,i] eq mtime0[2]) then begin
        ifile0 = ifile
        ihist0 = i
        print,format="('Found start mtime ',3i4,' on file ',a,' (history ',i3,')')",$
          mtime0,file_basename(files[ifile]),ihist0
      endif
    endelse
;
; If end day is < 0 or > 366, use last history of last file:
    if ((mtime1[0] lt 0 or mtime1[0] gt 366) and $
         ifile eq nfiles-1 and i eq ntime-1) then begin
      ifile1 = nfiles-1
      ihist1 = ntime-1
      print,format="('Using ending mtime ',3i4,' (last history on last file ',a,')')",$
        mtime[*,ihist1],file_basename(files[ifile1])
    endif else begin
      if (mtime[0,i] eq mtime1[0] and mtime[1,i] eq mtime1[1] and $
          mtime[2,i] eq mtime1[2]) then begin
        ifile1 = ifile
        ihist1 = i
        print,format="('Found end   mtime ',3i4,' on file ',a,' (history ',i3,')')",$
          mtime1,file_basename(files[ifile]),ihist1
      endif
    endelse
    if (ifile0 gt -1 and ifile1 gt -1) then break
  endfor ; i=0,ntime-1
  ncdf_close,ncids[ifile]
  if (ifile0 gt -1 and ifile1 gt -1) then break
endfor ; mtimes search: ifile=0,nfiles-1
if (ifile0 eq -1) then begin
  print,format="('>>> Could not find start mtime ',3i4,' on any file.')",mtime0
  return
endif
if (ifile1 eq -1) then begin
  print,format="('>>> Could not find end mtime ',3i4,' on any file.')",mtime1
  return
endif
;
; First file loop: get coord arrays and total number of times:
;
ntimes = 0
dtime1 = 0.
for ifile=ifile0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  finfo = ncdf_inquire(ncids[ifile])
  ndims = finfo.ndims
;
; Get spatial coord array dims from first file:
  if (ifile eq ifile0) then begin
    info1 = finfo
    ndims1 = info1.ndims
    dims1 = intarr(ndims1)
    for i=0,ndims1-1 do begin
      ncdf_diminq,ncids[ifile],i,name,len
      if (name eq 'lon') then nlon1=len
      if (name eq 'lat') then nlat1=len
      if (name eq 'lev') then nlev1=len
    endfor
  endif
;
; Coord array dims must be same as first file:
  for i=0,ndims-1 do begin
    ncdf_diminq,ncids[ifile],i,name,len
    for i=0,ndims-1 do begin
      ncdf_diminq,ncids[ifile],i,name,len
      if (name eq 'lon' and len ne nlon1) then begin
        print,'>>> lon dim (file 1)=',nlon1,' but lon dim (file ',ifile,')=',len
        return
      endif
      if (name eq 'lat' and len ne nlat1) then begin
        print,'>>> lat dim (file 1)=',nlat1,' but lat dim (file ',ifile,')=',len
        return
      endif
      if (name eq 'lev' and len ne nlev1) then begin
        print,'>>> lev dim (file 1)=',nlev1,' but lev dim (file ',ifile,')=',len
        return
      endif
    endfor
  endfor 
  for ivar = 0,finfo.nvars-1 do begin
    vinfo = ncdf_varinq(ncids[ifile],ivar)
    case vinfo.name of
      "mtime": ncdf_varget,ncids[ifile],ivar,mtime
      "lon"  : ncdf_varget,ncids[ifile],ivar,lon
      "lat"  : ncdf_varget,ncids[ifile],ivar,lat
      "lev"  : begin
        ncdf_varget,ncids[ifile],ivar,lev
        ncdf_attget,ncids[ifile],ivar,'units',units
        lev_units = string(units)
      end
      "time" : begin
        ncdf_varget,ncids[ifile],ivar,time
        ntime = n_elements(time)
        if (ifile eq ifile0) then ntime = ntime-ihist0
        if (ifile eq ifile1) then ntime = ihist1+1
        if (ifile0 eq ifile1 and ifile eq ifile0) then ntime = ihist1-ihist0+1
      end
      else:
    endcase
  endfor ; for ivar=0,nvars-1
  nlon = n_elements(lon)
  nlat = n_elements(lat)
  nlev = n_elements(lev)
  ntimes = ntimes+ntime
;
; Check for constant delta-time:
;
; if (ntime gt 0 and dtime1 eq 0.) then dtime1 = time[1]-time[0]
; for i=1,n_elements(time)-1 do begin
;   dtime = time[i]-time[i-1]
;   if (dtime ne dtime1) then begin
;     print,'>>> Delta-time has changed: dtime1=',dtime1,' current i=',i,$
;       ' time[i]=',time[i],' time[i-1]=',time[i-1],' dtime=',dtime
;     print,'Current file: ',files[ifile],' (ifile=',ifile,')'
;     return
;   endif
; endfor
  ncdf_attget,ncids[ifile],'missing_value',missing,/global
  ncdf_close,ncids[ifile]
  if (ifile eq ifile1) then break
endfor ; ifile=ifile0,nfiles-1
print,format="('nlon=',i4,' nlat=',i4,' nlev=',i4,' total ntimes=',i4)",$
  nlon,nlat,nlev,ntimes
;
;--------------------------------------------------------------------------
; For ut vs zp:
;
; Use nearest grid-point to each requested lat,lon location:
;
ilocs = intarr(2,nlocs) ; indices to nearest grid point for each location
gridlocs = fltarr(2,nlocs) ; nearest grid-point lat,lon for each location
for iloc=0,nlocs-1 do begin
  ilat = ixfind(lat,loc_lats[iloc],lat[1]-lat[0])
  ilon = ixfind(lon,loc_lons[iloc],lon[1]-lon[0])
  ilocs[0,iloc] = ilat
  ilocs[1,iloc] = ilon
  gridlocs(0,iloc) = [lat[ilat]] ; actual lat grid point
  gridlocs(1,iloc) = [lon[ilon]] ; actual lon grid point
endfor
;
; Set up zp array for y-axis:
;
dzp = lev[1]-lev[0] ; assume uniform vertical grid
if (zpbot eq zptop) then begin ; default to full zp grid
  zplev = lev
  nzp = nlev
  k0 = 0
  k1 = nlev-1
endif else begin
  if (zpbot lt lev[0]) then begin ; default to bottom of grid
    k0 = 0
  endif else begin
    k0 = ixfind(lev,zpbot,dzp)
    if (k0 eq -1) then begin
      print,'>>> Could not find grid point for zpbot=',zpbot
      print,'levels = ',lev
      return
    endif
  endelse
  if (zptop gt lev[nlev-1]) then begin ; default to top of grid
    k1 = nlev-1
  endif else begin
    k1 = ixfind(lev,zptop,dzp)
    if (k1 eq -1) then begin
      print,'>>> Could not find grid point for zptop=',zptop
      print,'levels = ',lev
      return
    endif
  endelse
  nzp = k1-k0+1
  zplev = fltarr(nzp)
  zplev[0:nzp-1] = lev[k0:k1]
endelse
print,format="('nzp=',i3,' k0,k1=',i3,',',i3,' zplev=',f9.3,' to ',f9.3)",$
  nzp,k0,k1,zplev[0],zplev[nzp-1]
;
; futv[nlocs,nzp,ntimes,nfields] for ut vs zp
; futlat[nslons,nslevs,ntimes,nfields] for ut vs lat
;
futv = dblarr(nlocs,nzp,ntimes,nfields)

nslons = n_elements(slons)
nslevs = n_elements(slevs)
;futlat = dblarr(nslons,nslevs,ntimes,nfields)

mtimes = intarr(3,ntimes)
years = intarr(ntimes)    ; integer years
decyears = fltarr(ntimes) ; decimal years
itime0 = 0
nyrbound = 0
iyrbound = intarr(20) ; assume max of 20 years
;
; File loop to read vars:
;
for ifile=ifile0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  finfo = ncdf_inquire(ncids[ifile])
;
; Get number of times (histories) on current file and update: 
;   ntime: Number of times in current file.
;   time0,time1: Beginning and ending time indices into 
;     total-time arrays (0->ntimes-1) for current file.
;
  ndims = finfo.ndims
  for idim=0,ndims-1 do begin
    ncdf_diminq,ncids[ifile],idim,dimname,dimsize
    if (dimname eq 'time') then begin
      ntime = dimsize 
;
; Number of histories to read from this file:
      if (ifile eq ifile0) then ntime = ntime-ihist0
      if (ifile eq ifile1) then ntime = ihist1+1
      if (ifile0 eq ifile1 and ifile eq ifile0) then ntime = ihist1-ihist0+1
      itime1 = itime0+ntime-1 
    endif
  endfor
;
; First,last histories to read from this file:
;
  ihistbeg = 0 & ihistend = ntime-1
  if (ifile eq ifile0) then begin
    ihistbeg = ihist0 & ihistend = ihist0+ntime-1
  endif else if (ifile eq ifile1) then begin
    ihistbeg = 0 & ihistend = ihist1
  endif
; print,format="(/,'itime0,1=',2i4,' ihistbeg,end=',2i4,' ntime=',i4)",$
;   itime0,itime1,ihistbeg,ihistend,ntime
;
; Get years:
;
  id = ncdf_varid(ncids[ifile],'year')
  if (id eq -1) then print,'>>> WARNING: could not find year variable'
  ncdf_varget,ncids[ifile],id,var
; print,' ' & print,'itime0=',itime0,' itime1=',itime1,' ntime=',ntime
  years[itime0:itime1] = fix(var[ihistbeg:ihistend])
;
; Get model times:
;
  id = ncdf_varid(ncids[ifile],'mtime')
  if (id eq -1) then print,'>>> WARNING: could not find mtime variable'
  ncdf_varget,ncids[ifile],id,var
  mtimes[*,itime0:itime1] = var[*,ihistbeg:ihistend]
;
; Check for year-boundary:
;
  for i=itime0,itime1 do begin
    decyears[i] = float(years[i])+float(mtimes[0,i])/365.
    if (decyears[i] eq float(fix(decyears[i]))) then begin
      print,'Found year-boundary at i=',i,' decyears=',decyears[i]
      iyrbound[nyrbound] = i
      nyrbound = nyrbound+1
    endif
  endfor
;
; Print time info for current file:
;
  print,format="(/,'File ',a,' (',i4,' of ',i4,'):')",$
    file_basename(files[ifile]),ifile-ifile0+1,ifile1-ifile0+1
  print,format=$
    "('Reading ',i3,' model times: ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3)",$
    ntime,mtimes[*,itime0],mtimes[*,itime1]
;
; Read requested 4d fields:
;
  found = intarr(nfields)
  found[*] = 0
  for ivar = 0,finfo.nvars-1 do begin
    vinfo = ncdf_varinq(ncids[ifile],ivar)
    for ifld=0,nfields-1 do begin
      if (vinfo.name eq fields[ifld]) then begin
        found[ifld] = 1
        if (vinfo.ndims ne 4) then begin ; should also verify dim lens
          print,'>>> Field ',fields[ifld],' is not 4-d: skipping..'
          futv[*,*,itime0:itime1,ifld] = 0.
          break
        endif
        ncdf_varget,ncids[ifile],ivar,var ; (lon,lat,lev,time)
        print,format="('Field ',a,': 4d min,max=',2e12.4)",fields[ifld],$
          min(var[*,*,k0:k1,ihistbeg:ihistend]),max(var[*,*,k0:k1,ihistbeg:ihistend])
        for iloc=0,nlocs-1 do begin
;
; f = fltarr(nlocs,nlev,ntimes,nfields)
          futv[iloc,*,itime0:itime1,ifld] = $
            var[ilocs[1,iloc],ilocs[0,iloc],k0:k1,ihistbeg:ihistend]
          if (k1 eq nlev-1) then $
            futv[iloc,nzp-1,itime0:itime1,ifld] = futv[iloc,nzp-2,itime0:itime1,ifld]
        endfor
;
; Check for units attribute of current variable (may not exist for all fields):
        f_units[ifld] = ""
        for iatt=0,vinfo.natts-1 do begin
          attname = ncdf_attname(ncids[ifile],ivar,iatt)
          if (attname eq "units") then begin
            ncdf_attget,ncids[ifile],ivar,attname,units
            f_units[ifld]=string(units)
          endif
        endfor
      endif ; is a requested field
    endfor ; requested field loop
  endfor ; ivar=0,nvars-1
  for ivar=0,nfields-1 do begin
    if (found[ivar] eq 0) then print,format=$
      "('(Field ',a,' was not found)')",fields[ivar]
  endfor
  itime0 = itime0+ntime
  ncdf_close,ncids[ifile]
  if (ifile eq ifile1) then break
endfor ; ifile=ifile0,nfiles-1
if (ntimes le 1) then begin
  print,'>>> ntimes=',ntimes,': must have more than 1 time.'  
  return
endif
;
; Insure consistent delta-time (minutes) between all histories.
;
; If a deltamin between 2 histories is zero (duplicate history), 
; then the arrays are shifted to overwrite the 1st history of the 
; duplicate pair, and ntimes is decremented.  Duplicate secondary 
; histories can occur between the end of one file and the beginning 
; of the next when the model is restarted.  If deltamin is inconsistent 
; w/ deltamin0 and deltamin is non-zero, then this is a fatal
; (the history write frequency was apparently changed mid-run).
;
deltamin0 = mtime_to_mins(mtimes[*,1])-mtime_to_mins(mtimes[*,0])
dtimeloop:
for i=1,ntimes-1 do begin
  deltamin = mtime_to_mins(mtimes[*,i])-mtime_to_mins(mtimes[*,i-1])
  if (deltamin ne deltamin0) then begin
    print,format="(/,'Warning: deltamin is inconsistent: deltamin0=',i5,' deltamin=',i5)",$
      deltamin0,deltamin
    print,format="('      i=',i5,' mtimes[i]=',3i4,' mtimes[i-1]=',3i4)",$
      i,mtimes[*,i],mtimes[*,i-1]
    if (deltamin ne 0) then return
    print,format="('      (Duplicate histories will be corrected by shifting the arrays)')"
    mtimes[*,i:ntimes-1] = shift(mtimes[*,i:ntimes-1],0,-1)
    years[i:ntimes-1] = shift(years[i:ntimes-1],-1)
    decyears[i:ntimes-1] = shift(decyears[i:ntimes-1],-1)
;
; f = dblarr(nlocs,nzp,ntimes,nfields)
; The shift function does not work if there are any degenerate dimensions
;   (i.e., dimensioned 1). Since nfields and nlocs can be degenerate
;   (only one chosen by the user), loop over these to avoid problems 
;   w/ shift function:
;
    for ifld=0,nfields-1 do begin
      for iloc=0,nlocs-1 do begin
        futv[iloc,*,i:ntimes-1,ifld] = shift(futv[iloc,*,i:ntimes-1,ifld],0,0,-1)
      endfor
    endfor
;
; This results in an error in the number of arguments to shift function:
;   futv[*,*,i:ntimes-1,*] = shift(futv[*,*,i:ntimes-1,*],0,0,-1,0)
    ntimes = ntimes-1
    goto,dtimeloop ; check for another one
  endif
endfor
;
print,format=$
  "(/,'Total ntimes=',i5,': model times ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3,' by ',i6,' minutes')",$
  ntimes,mtimes[*,0],mtimes[*,ntimes-1],deltamin
;
; Graphics device is X or ps:
;
if (keyword_set(ps)) then begin
  set_plot,'PS'
; device,/color,/landscape            ; default width=9.5" (prints on 8.5x11)
; device,/color,/landscape,/inches,xsize=11.,ysize=8.  ; 11" wide
; device,/color,/landscape,/inches,xsize=14.,ysize=8.  ; fits 1 yr w/ scrolling
  device,/color,/landscape,/inches,xsize=25.,ysize=8.  ; best for > 2-3 years
; device,/color,/landscape,/inches,xsize=35.,ysize=8.
endif else begin
  set_plot,'X'
; Could see not difference between default and TrueType font
  device,decomposed=0,retain=2,set_font='Times',/tt_font ; TrueType font
; device,decomposed=0,retain=2                           ; for TrueColor display
  window,xsize=650,ysize=450,retain=2
endelse

 nlevels = 12 ; number of contour levels
 loadct,33,ncolors=nlevels,bottom=3 ; Blue-Red
;
; Field plot loop:
;
iframe = 0
ijpg = 0
for ifld=0,nfields-1 do begin
  if (found[ifld] eq 0) then begin
    print,'Skipping field ',fields[ifld],' because it was not found.'
    continue ; field was not found
  endif
  for iloc=0,nlocs-1 do begin
    utv = futv[iloc,*,0:ntimes-1,ifld] ; [nlev,ntimes]
    utv = transpose(utv)   ; [ntimes,nlev]
    funits = f_units[ifld]
    if (f_units[ifld] eq 'CM') then begin
      utv = utv/1.e5          ; cm to km (e.g., Z)
      funits = 'KM'
    endif
    if (f_units[ifld] eq 'CM/S') then begin
      utv = utv/100.
      funits = 'M/S'
    endif
;
; Main title and y-axis title:
;
    title = fields[ifld]+'('+funits+'):  Lat,Lon = '+ $
      string(format="(f6.2)",gridlocs[0,iloc])+', '+$
      string(format="(f7.2)",gridlocs[1,iloc])
    ytitle = "Log Pressure ("+lev_units+")"
;
; Use IDL procedures timegen, julday, and label_date 
; (setting nxticks does not appear to work w/ label_date)
;
    time = timeaxis(years[0:ntimes-1],mtimes[*,0:ntimes-1],deltamin,xtitle,$
      nxticks)
    filelabel = "Files " + file0 + " to " + file1
;
; Contour:
    utcontour,utv,time,zplev,loc_lons[iloc],nlevels,xtitle,ytitle,title,filelabel,missing,$
      jpg=jpg,ps=ps
    iframe = iframe+1
    print,format="('Frame ',i3,': ',a)",iframe,title
;
; Jpegs work, but only one image per file.
    if keyword_set(jpg) then begin
      image = tvrd(true=3)
      ijpg = ijpg+1
;     jpg_outfile = string(format="(a,'/image',i3.3,'.jpg')",image_dir,ijpg)
      jpg_outfile = strcompress(string(format=$
        "(a,'/Frame',i4.4,'_',a,'_utzp_lat',f6.2,'_lon',f7.2,'.jpg')",$
        image_dir,ijpg,fields[ifld],loc_lats[iloc],loc_lons[iloc]),/remove_all)
      print,format="('Write jpeg file ',a)",jpg_outfile
      help,image
      write_jpeg,jpg_outfile,image,true=3,quality=100
    endif

  endfor ; iloc
endfor ; ifld
if (keyword_set(ps)) then begin
  device,/close
endif
end