pro pltutlat,stateptr,utlatptr
  utlat = *utlatptr
  state = *stateptr
  sfields = *state.sfields
  nfields = n_elements(sfields)
  f_units = strarr(nfields)
  f_units[*] = ""

  lats = *state.lats
  lons = *state.lons
; levs = *state.levs
  levs = *state.ilevs

  zplevs = *utlat.zplevs
  nzplev = n_elements(zplevs)
  print,'pltutlat: nzplev=',nzplev,' zplevs=' & print,zplevs

  slons = *utlat.slons
  nslon = n_elements(slons)
  print,'pltutlat: nslon=',nslon,' slons=' & print,slons

  ilat0 = ixfind_simple(lats,utlat.lat0)
  ilat1 = ixfind_simple(lats,utlat.lat1)
  ylats = lats[ilat0:ilat1]
  nlats = n_elements(ylats)
;
; Indices into mtimes of mtime_beg and mtime_end, and
; number of times to plot:
;
  mtimes = *state.mtimes
  ntimes_total = state.ntimes
  for i=0,ntimes_total-1 do begin
    if (array_equal(mtimes[*,i],state.mtime_beg)) then begin
      it0 = i
      break
    endif
  endfor
  for i=0,ntimes_total-1 do begin
    if (array_equal(mtimes[*,i],state.mtime_end)) then begin
      it1 = i
      break
    endif
  endfor
  ntime_plot = it1-it0+1
  mtime_plot = intarr(3,ntime_plot)
  for it=it0,it1 do mtime_plot[*,it-it0] = mtimes[*,it]
  for it=1,ntime_plot-1 do begin
    deltamin = mtime_to_mins(mtime_plot[*,it])-mtime_to_mins(mtime_plot[*,it-1])
    if deltamin ne 0 then break
  endfor
  print,format=$
    "('Reading ',i5,' times: ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3,' by ',i4)",$
    ntime_plot,mtime_plot[*,0],mtime_plot[*,ntime_plot-1],deltamin
;
; Allocate utlat array:
;
  futlat   = dblarr(nlats,nslon,nzplev,ntime_plot,nfields)
  years    = intarr(ntime_plot)  ; integer years
  days     = intarr(ntime_plot)  ; integer year-days
  ut       = dblarr(ntime_plot)  ; ut
  decyears = dblarr(ntime_plot)  ; decimal years
;
; Read from files:
;
  files = *state.files
  ifile0 = state.ifile_beg
  ifile1 = state.ifile_end
  itime0_plot = 0
  msis_version = "NRL MSIS-2000"
  for ifile=ifile0,ifile1 do begin
    ncid = ncdf_open(files[ifile],/nowrite)

    id = ncdf_varid(ncid,'mtime')
    ncdf_varget,ncid,id,mtime_file
    ntime_file = n_elements(mtime_file)/3

    ncdf_attget,ncid,'missing_value',missing,/global
    ncdf_attget,ncid,'model_version',tgcm_version,/global
    tgcm_version = string(tgcm_version)
    print,'tgcm_version=',tgcm_version

    itime0_file = 0 & itime1_file = ntime_file-1
    for i=0,ntime_file-1 do begin
      if (array_equal(mtime_file[*,i],state.mtime_beg)) then itime0_file=i
      if (array_equal(mtime_file[*,i],state.mtime_end)) then itime1_file=i
    endfor
    ntime_file = itime1_file-itime0_file+1
    itime1_plot = itime0_plot+ntime_file-1

;   print,format=$
;     "(/,'File ',a,' (',i3,' of ',i3,')')",$
;     file_basename(files[ifile]),ifile+1,ifile1-ifile0+1
;   print,format=$
;     "('Utlat reading ',i4,' mtimes: ',3(i3,','),' to ',3(i3,','),' by ',i5,')')",$
;     ntime_file,mtime_plot[*,itime0_plot],mtime_plot[*,itime1_plot],deltamin

    text = '-------------------------------------------------------------------------------------------------'
    myprint,text,/stdout,outfile=state.outfile,/blankline

    text=string(format=$
      "(/,'File ',a,' (',i3,' of ',i3,') Resolution=',f4.1,' (deg)')",$
      file_basename(files[ifile]),ifile+1,ifile1-ifile0+1,state.horiz_resolution)
    myprint,text,/stdout,outfile=state.outfile,/blankline

    text=string(format=$
      "('Utlat reading ',i4,' mtimes: ',3(i3,','),' to ',3(i3,','),' by ',i5,')')",$
      ntime_file,mtime_plot[*,itime0_plot],mtime_plot[*,itime1_plot],deltamin)
    myprint,text,/stdout,outfile=state.outfile,/blankline
;
; Get years:
    id = ncdf_varid(ncid,'year') 
    if (id eq -1) then print,'>>> WARNING: could not find year variable'
    ncdf_varget,ncid,id,var
    years[itime0_plot:itime1_plot] = fix(var[itime0_file:itime1_file])
;
; Get days:
    id = ncdf_varid(ncid,'day')
    if (id eq -1) then print,'>>> WARNING: could not find day variable'
    ncdf_varget,ncid,id,var
    days[itime0_plot:itime1_plot] = fix(var[itime0_file:itime1_file])
;
; Get ut:
    id = ncdf_varid(ncid,'ut')
    if (id eq -1) then print,'>>> WARNING: could not find ut variable'
    ncdf_varget,ncid,id,var
    ut[itime0_plot:itime1_plot] = fix(var[itime0_file:itime1_file])
;
; Calculate decimal years:
    for i=itime0_plot,itime1_plot do begin
      decyears[i] = double(years[i])+double(days[i]-1)/365.+ut[i]/(365.*24.)
    endfor

    density_conversion = 0
    density_log10 = 0
;
; Loop over variables in current file:
    finfo = ncdf_inquire(ncid)
;
; Get tgcm Z geopotential and tgcm fields for density conversion 
; if calling msis.
; futlat   = dblarr(nlats,nslon,nzplev,ntime_plot,nfields)
;
    density_conversion = 0
    density_log10 = 0
    model = ''
    for i=0,nfields-1 do begin
      if (sfields[i].model eq 'MSIS') then model = 'MSIS'
    endfor
    if (model eq 'MSIS') then begin ; there is at least 1 msis field
      ier = get_tgcm_var(ncid,'Z',zpoten_read) 
      zpoten = dblarr(nlats,ntime_file)
;
; Check if unit conversion of tgcm density fields is necessary
; (i.e., if any density species is requested for both models)
; Note that fields_msis[*].density was been set in utproc.
;
      density_conversion = 0
      density_log10 = 0
      for i=0,nfields-1 do begin
        if (sfields[i].model eq 'MSIS' and sfields[i].density eq 1) then begin
          msisname = sfields[i].name
          for i=0,nfields-1 do begin
            if (sfields[i].model eq 'TGCM' and sfields[i].name eq msisname) then $
              density_conversion = 1
          endfor
        endif
      endfor
;
; Get tgcm fields necessary for density unit conversion to cm-3:
;
      if (density_conversion) then begin
        ixt  = get_tgcm_var(ncid,'TN',tn_read)
        ixo2 = get_tgcm_var(ncid,'O2',o2_read)
        ixo  = get_tgcm_var(ncid,'O1',o1_read)
        if (ixt lt 0 or ixo2 lt 0 or ixo lt 0) then begin
          print,'>>> WARNING: could not find all tgcm fields necessary for ',$
            'density conversion: ixt,ixo2,ixo=',ixt,' ',ixo2,' ',ixo
          density_conversion = 0
        endif
        ixp0 = get_tgcm_var(ncid,'p0_model',p0)
        if (ixp0 lt 0) then ixp0 = get_tgcm_var(ncid,'p0',p0) 
        if (ixp0 lt 0) then begin
          print,'>>> WARNING: could not find p0.'
          density_conversion = 0 
        endif
        if (density_conversion) then begin
          tn = dblarr(nlats,ntime_file)
          o2 = dblarr(nlats,ntime_file)
          o1 = dblarr(nlats,ntime_file)
          density_log10 = 1
        endif
      endif ; density conversion needed
      msis_data = dblarr(nlats,ntime_file)

      print,'pltutlat: zpoten_read min,max=',min(zpoten_read),', ',$
        max(zpoten_read),' density_conversion=',density_conversion

    endif ; plotting at least one msis field

    for ivar=0,finfo.nvars-1 do begin
      vinfo = ncdf_varinq(ncid,ivar)
;
; Check for selected field: 
      for ifld=0,nfields-1 do begin
        if (sfields[ifld].model ne 'TGCM') then continue ; want tgcm only
        if (vinfo.name eq sfields[ifld].name) then begin
          for iatt=0,vinfo.natts-1 do begin
            attname = ncdf_attname(ncid,ivar,iatt)
            if (attname eq "units") then begin
              ncdf_attget,ncid,ivar,attname,units
              f_units[ifld]=string(units)
            endif
          endfor
          ncdf_varget,ncid,ivar,var ; (lon,lat,lev,time)

          if (f_units[ifld] eq 'CM/S' or f_units[ifld] eq 'cm/s') then begin
            var = var/100.
            f_units[ifld] = 'M/S'
          endif
;
; Location (lat,lon) loop:
; futlat = dblarr(nlats,nslon,nzplev,ntime_plot,nfields)
;
          for k=0,nzplev-1 do begin
            klev = ixfind_simple(levs,zplevs[k])
            for i=0,nslon-1 do begin
              ilon = ixfind_simple(lons,slons[i])
              futlat[*,i,k,itime0_plot:itime1_plot,ifld] = $
                var[ilon,ilat0:ilat1,klev,itime0_file:itime1_file]
            endfor ; i=0,nslon-1
          endfor ; k=0,nzplev-1

;         print,format=$
;           "('Read field ',a8,' min,max=',2e12.4)",$
;           sfields[ifld].name,min(futlat[*,*,*,itime0_plot:itime1_plot,ifld]),$
;                              max(futlat[*,*,*,itime0_plot:itime1_plot,ifld])

          text=string(format=$
            "('Read field ',a8,' min,max=',2e12.4)",$
            sfields[ifld].name,min(futlat[*,*,*,itime0_plot:itime1_plot,ifld]),$
                               max(futlat[*,*,*,itime0_plot:itime1_plot,ifld]))
          myprint,text,/stdout,outfile=state.outfile

        endif ; is selected field
      endfor ; ifld=0,nfields-1
    endfor ; ivar=0,finfo.nvars-1 for tgcm fields
;
; Get msis fields:
; zpoten = dblarr(nlats,ntime_file)
;
    for ifld=0,nfields-1 do begin ; all selected fields
      if (sfields[ifld].model ne 'MSIS') then continue
      for k=0,nzplev-1 do begin                 ; selected levels
        klev = ixfind_simple(levs,zplevs[k])
        for i=0,nslon-1 do begin                ; selected lons
          ilon = ixfind_simple(lons,slons[i])
          zpoten(*,*) = zpoten_read[ilon,*,klev,*]

; function msis_utlat,field,slon,slev,lats,years,days,mtimes,zpoten,out
;     msis_data = dblarr(nlats,ntime_file)

;         print,'pltutlat call msis: field ',sfields[ifld].name,' slons[i]=',$
;           slons[i],' zplevs[k]=',zplevs[k]

          found_field = msis_utlat(sfields[ifld].name,slons[i],zplevs[k],ylats,$
            years(itime0_plot:itime1_plot),days(itime0_plot:itime1_plot),$
            mtime_plot(*,itime0_plot:itime1_plot),zpoten,msis_data)

;         print,'pltutlat after call msis: field ',sfields[ifld].name,' min,max=',$
;           min(msis_data),max(msis_data)

          if (sfields[ifld].density and density_log10) then msis_data = alog10(msis_data)
          futlat[*,i,k,itime0_plot:itime1_plot,ifld] = msis_data[*,*]

; From pltutvert:
;       found_field = msis_utvert(sfields[ifld].name,locs[0:1,iloc],$
;         zplevs[izp0:izp1],years(itime0_plot:itime1_plot),$
;         days(itime0_plot:itime1_plot),mtime_plot(*,itime0_plot:itime1_plot),$
;         zpoten,msis_data)
;       if (sfields[ifld].density and density_log10) then msis_data = alog10(msis_data)
;       futv[iloc,*,itime0_plot:itime1_plot,ifld] = msis_data[*,*]

        endfor ; i=0,nslon-1
      endfor ; k=0,nzplev-1
    endfor ; ifld=0,nfields-1 for msis fields

    itime0_plot = itime0_plot+ntime_file
    ncdf_close,ncid
  endfor ; ifile=ifile0,ifile1
;
; 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).
;
  dtimeloop:
  for i=1,ntime_plot-1 do begin
    deltamin = mtime_to_mins(mtime_plot[*,i])-mtime_to_mins(mtime_plot[*,i-1])
    if (deltamin eq 0) then begin
      print,format=$
        "(/,'Warning: duplicate history: i=',i3,',mtime[*,i]=',3i4,' mtime[*,i-1]=',3i4)",$
        i,mtime_plot[*,i],mtime_plot[*,i-1]
      print,format=$
        "('(Duplicate histories are corrected by shifting arrays)')"
      mtime_plot[*,i:ntime_plot-1] = shift(mtime_plot[*,i:ntime_plot-1],0,-1)
      years[i:ntime_plot-1] = shift(years[i:ntime_plot-1],-1)
      decyears[i:ntime_plot-1] = shift(decyears[i:ntime_plot-1],-1)
;
; The shift function does not work if there are any degenerate dimensions
;   (i.e., dimensioned 1). Since nfields, nlons, and nzplevs can be degenerate
;   (only one chosen by the user), loop over these to avoid problems
;   w/ shift function:
; futlat = dblarr(nlats,nslon,nzplev,ntime_plot,nfields)
;
      for ifld=0,nfields-1 do begin
        for k=0,nzplev-1 do begin
          for islon=0,nslon-1 do begin
            futlat[*,islon,k,i:ntime_plot-1,ifld] = $
              shift(futlat[*,islon,k,i:ntime_plot-1,ifld],0,0,0,-1)
          endfor
        endfor
      endfor
      ntime_plot = ntime_plot-1
      goto,dtimeloop ; check for another one
    endif
  endfor
;
; Set output device (as of 12/4/07 only PS and X are implemented,
; (but see non-gui utproc.pro for jpg)
;
  print,format="(/,'Contouring ut vs latitude:')"
; set_output_device,stateptr
;
; Number of contour levels, and load default color table:
; nlevels = 12 ; number of contour levels
; loadct,33,ncolors=nlevels,bottom=3 ; Blue-Red

  iframe = 0
  ijpg = 1
  ipng = 1
; ps = 0 & if (state.device eq 'PS') then ps = 1
; jpg = 0 & if (state.device eq 'JPG') then jpg = 1
;
; Field plot loop:
; futlat = dblarr(nlats,nslon,nzplev,ntime_plot,nfields)
;
for ifld=0,nfields-1 do begin

  text = string(format=$
    "('-------------------- Utlat: Begin plotting field ',a,' -------------------')",$
    sfields[ifld].name)
  myprint,text,/stdout,outfile=state.outfile,/blankline

  model = sfields[ifld].model
  utlat = dblarr(ntime_plot,nlats)
  for k=0,nzplev-1 do begin
    for i=0,nslon-1 do begin
      for j=0,nlats-1 do begin
        for it=0,ntime_plot-1 do begin
          utlat[it,j] = futlat[j,i,k,it,ifld]
        endfor
      endfor
      funits = f_units[ifld]
      if (f_units[ifld] eq 'CM' or f_units[ifld] eq 'cm') then begin
        utlat = utlat/1.e5          ; cm to km (e.g., Z)
        funits = 'KM'
      endif
;
; Main title and y-axis title:
;
      if (model eq 'TGCM') then version = tgcm_version else version = msis_version
      if (sfields[ifld].density and density_log10) then $
        title_log10 = 'log10 ' else title_log10 = ''
      title = version+' '+sfields[ifld].name+' ('+funits+'):  ZP level = '+ $
        string(format="(f6.2)",zplevs[k])+' Longitude = '+ $
        string(format="(f7.2)",slons[i])+ $
        ' Resolution = '+string(format="(f3.1,' (deg)')",state.horiz_resolution)
      ytitle = "Latitude"
;
; Use IDL procedures timegen, julday, and label_date
; (setting nxticks does not appear to work w/ label_date)
;
      time = timeaxis(decyears[0:ntime_plot-1],years[0:ntime_plot-1],$
        mtime_plot[*,0:ntime_plot-1],deltamin,xtitle,nxticks)
      filelabel = "Files " + file_basename(files[ifile0]) + " to " + $
                             file_basename(files[ifile1])
      ps = 1
;     ps = 0 & if (state.device eq 'PS') then ps = 1
;     jpg = 0 & if (state.device eq 'JPG') then jpg = 1
;     cmin=0. & cmax=0. & cint=0. 

      cmin=sfields[ifld].cmin
      cmax=sfields[ifld].cmax
      cint=sfields[ifld].cint
;     print,'pltutlat: sfields cmin,cmax,cint=',cmin,',',cmax,',',cint

;     loadct,33,ncolors=nlevels,bottom=3 ; Blue-Red
;
; Contour:
      utcontour,stateptr,utlat,time,ylats,slons[i],cmin,cmax,cint,$
        xtitle,ytitle,title,filelabel,missing,jpg=jpg,ps=ps,png=ipng
      if (title ne '') then begin
        iframe = iframe+1

;       if (state.device eq 'PS') then begin
;         print,format="('Frame ',i4,': ',a,' (',a,')')",$
;           iframe,title,strcompress(state.ps_outfile,/remove_all)
;       endif else if (state.device eq 'X') then begin
;         print,format="('Frame ',i4,': ',a,' (X display)')",iframe,title
;       endif else begin
;         print,format="('Frame ',i4,': ',a,' (unknown output device)')",$
;           iframe,title
;       endelse

;       print,format="('Frame ',i4,': ',a,' (ps file ',a,')')",$
;         iframe,title,strcompress(state.ps_outfile,/remove_all)

        text=string(format="('Frame ',i4,': ',a,' (',a,')')",$
          iframe,title,strcompress(state.ps_outfile,/remove_all))
        myprint,text,/stdout,outfile=state.outfile

      endif
    endfor ; i=0,nslon-1
  endfor ; k=0,nzplev-1
endfor ; ifld=0,nfields-1
end