function denconv,name,varloc,tn,o2,o1,p0,levs
;
; Convert tgcm density field from mmr to cm3
;
; varloc[nzplev,ntime_file] is input density field at current loc,
; to be converted from mmr to cm3. tn,o2,o1 are inputs at same loc
; (o2,o1 inputs are in mmr)
;
dimsizes = size(varloc,/dimensions)
nlev = dimsizes[0]
ntime = dimsizes[1]
pkt  = dblarr(nlev,ntime)
barm = dblarr(nlev,ntime)
boltz = 1.3805e-16
for i=0,ntime-1 do begin
  for k=0,nlev-1 do begin
    barm[k,i] = 1.0/(o2[k,i]/32. + o1[k,i]/16. + $
      max([.00001,(1.-o2[k,i]-o1[k,i])])/28.)
    pkt[k,i] = p0*exp(-levs[k])/(boltz*tn[k,i])
  endfor
endfor
case name of
  'O2': wt = 32.
  'O1': wt = 16.
  'N' : wt = 12.
  'N2': wt = 24.
  else: begin
    print,'>>> WARNING denconv: do not know molecular weight of field ',name
    wt = 1.
  end
endcase
out = varloc * pkt * barm / wt
print,'denconv: field ',name,' cm3 min,max=',min(out),max(out)
return,out
end
;-----------------------------------------------------------------------
pro pltutvert,stateptr,utvstateptr
  utvstate = *utvstateptr
  state = *stateptr
  sfields = *state.sfields
  if (isnullfield(sfields)) then begin
    print,'>>> No fields selected.'
    return
  endif
  nfields = n_elements(sfields)
  f_units = strarr(nfields)
  lats = *state.lats
  lons = *state.lons
; levs = *state.levs
  levs = *state.ilevs
  nlevs = n_elements(levs)
;
; Get locs and indices into lat,lon coord arrays:
;
  nlocs = utvstate.nlocs
; print,format="('pltutvert: nlocs=',i3,' locs (lat,lon)=')",nlocs
  locs = utvstate.locs[*,0:nlocs-1]
  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(lats,locs[0,iloc],lats[1]-lats[0])
    ilon = ixfind(lons,locs[1,iloc],lons[1]-lons[0])
    ilocs[0,iloc] = ilat
    ilocs[1,iloc] = ilon
    gridlocs(0,iloc) = [lats[ilat]] ; actual lat grid point
    gridlocs(1,iloc) = [lons[ilon]] ; actual lon grid point
;   print,format="('(',2f8.2,')')",locs[*,iloc]
  endfor
;
; Get zp levels and start/end indices into levs:
;
  zplevs = *utvstate.zplevs
  nzplev = n_elements(zplevs)
  for k=0,nlevs-1 do begin
    if (levs[k] eq zplevs[0]) then izp0 = k
    if (levs[k] eq zplevs[nzplev-1]) then izp1 = k
  endfor 
; print,format="('pltutvert: izp0,1=',2i3,' nzplev=',i3,' zplevs=',/,(6f7.2))",$
;   izp0,izp1,nzplev,zplevs
;
; 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=$
;   "(/,'Utvert reading ',i5,' times: ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3,' by ',i4)",$
;   ntime_plot,mtime_plot[*,0],mtime_plot[*,ntime_plot-1],deltamin

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

  text=string(format=$
    "(/,'Utvert reading ',i5,' times: ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3,' by ',i4)",$
    ntime_plot,mtime_plot[*,0],mtime_plot[*,ntime_plot-1],deltamin)
  myprint,text,/stdout,outfile=state.outfile,/blankline

;
; Allocate utv array:
;
  futv  = dblarr(nlocs,nzplev,ntime_plot,nfields) ; plot array
  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=$
;     "('Reading ',i4,' mtimes: ',3(i3,','),' to ',3(i3,','),' by ',i5,')')",$
;     ntime_file,mtime_plot[*,itime0_plot],mtime_plot[*,itime1_plot],deltamin

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

    text=string(format=$
      "('Utvert 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

;
; 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
;
; Get info structure for current file:
    finfo = ncdf_inquire(ncid)
;
; Get tgcm Z geopotential and tgcm fields for density conversion 
; if calling msis.
;
    density_conversion = 0
    density_log10 = 1
    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(nzplev,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(nzplev,ntime_file)
          o2 = dblarr(nzplev,ntime_file)
          o1 = dblarr(nzplev,ntime_file)
;         density_log10 = 1
        endif
      endif ; density conversion needed
      msis_data = dblarr(nzplev,ntime_file)
    endif ; plotting at least one msis field
    varloc = dblarr(nzplev,ntime_file)
;
; Get tgcm fields:
;
    for ivar=0,finfo.nvars-1 do begin ; all variables in history file
      vinfo = ncdf_varinq(ncid,ivar)
;
; Check for selected field: 
      for ifld=0,nfields-1 do begin   ; all selected fields
        if (sfields[ifld].model ne 'TGCM') then continue ; want tgcm only
        if (vinfo.name eq sfields[ifld].name) then begin ; found selected field
          f_units[ifld] = ""
          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:
          for iloc=0,nlocs-1 do begin
            varloc[*,*] = var[ilocs[1,iloc],ilocs[0,iloc],izp0:izp1,itime0_file:itime1_file]
            if (density_conversion and sfields[ifld].density eq 1) then begin
              tn[*,*] = tn_read[ilocs[1,iloc],ilocs[0,iloc],izp0:izp1,itime0_file:itime1_file]
              o2[*,*] = o2_read[ilocs[1,iloc],ilocs[0,iloc],izp0:izp1,itime0_file:itime1_file]
              o1[*,*] = o1_read[ilocs[1,iloc],ilocs[0,iloc],izp0:izp1,itime0_file:itime1_file]
              varloc = denconv(sfields[ifld].name,varloc,tn,o2,o1,p0,levs[izp0:izp1])
              f_units[ifld]='cm-3'
;             if (density_log10) then varloc = alog10(varloc)
              if (density_log10) then log10f,varloc,missing
            endif
            futv[iloc,*,itime0_plot:itime1_plot,ifld] = varloc[*,*]
          endfor ; iloc=0,nlocs-1

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

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

        endif ; a selected field
      endfor ; ifld=0,nfields-1
    endfor ; ivar=0,nvars-1 ; for tgcm fields
;
; Get msis fields:
;
    for ifld=0,nfields-1 do begin ; all selected fields
      if (sfields[ifld].model ne 'MSIS') then continue
      for iloc=0,nlocs-1 do begin
;
; Get tgcm geopotential at this location from Z field, which was read into
; zpoten_read from tgcm history file (see above):
        zpoten(*,*) = zpoten_read(ilocs[1,iloc],ilocs[0,iloc],izp0:izp1,*)
;
; Procedure msis_utvert calls msis for current date and time at current location:
        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 log10f,msis_data,missing
        futv[iloc,*,itime0_plot:itime1_plot,ifld] = msis_data[*,*]
      endfor ; iloc=0,nlocs-1
      print,format=$
        "('Call ',a,' for field ',a8,' min,max=',2e12.4)",msis_version,$
        sfields[ifld].name,min(futv[*,*,itime0_plot:itime1_plot,ifld]),$
                           max(futv[*,*,itime0_plot:itime1_plot,ifld])
    endfor ; ifld=0,nfields-1 for msis fields
;
; Update start time for next file:
    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)
;
; f = dblarr(nlocs,nzp,ntime_plot,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:ntime_plot-1,ifld] = shift(futv[iloc,*,i:ntime_plot-1,ifld],0,0,-1)
        endfor
      endfor
      ntime_plot = ntime_plot-1
      goto,dtimeloop ; check for another one
    endif
  endfor
;
; Number of contour levels, and load default color table:
; nlevels = 12 ; number of contour levels
; loadct,33,ncolors=nlevels,bottom=3 ; Blue-Red
;
; Field plot loop:
;
print,format="(/,'Contouring ut vs pressure:')"
iframe = 0
ijpg = 0
for ifld=0,nfields-1 do begin

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

  model = sfields[ifld].model
  for iloc=0,nlocs-1 do begin
    utv = futv[iloc,*,0:ntime_plot-1,ifld] ; [nlev,ntimes]
    utv = transpose(utv)   ; [ntime_plot,nlev]
    funits = f_units[ifld]
    if (f_units[ifld] eq 'CM' or f_units[ifld] eq 'cm') then begin
      utv = utv/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+' '+title_log10+sfields[ifld].name+' ('+funits+'):  Lat,Lon = '+ $
      string(format="(f6.2)",gridlocs[0,iloc])+', '+$
      string(format="(f7.2)",gridlocs[1,iloc])

;   ytitle = "Log Pressure ("+lev_units+")" ; get this later
    ytitle = "Log Pressure"
;
; Use IDL procedures timegen, julday, and label_date
; (setting nxticks does not appear to work w/ label_date)
;
;   print,'call timeaxis: ntime_plot=',ntime_plot,' decyears[0:ntime_plot-1]=',$
;     decyears[0:ntime_plot-1]

    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,'pltutvert: sfields cmin,cmax,cint=',cmin,',',cmax,',',cint
;
; Contour:
    utcontour,stateptr,utv,time,zplevs,locs[1,iloc],cmin,cmax,cint,$
      xtitle,ytitle,title,filelabel,missing,jpg=jpg,ps=ps
    iframe = iframe+1
    latstring = strcompress(string(format="(f6.2)",gridlocs[0,iloc]),/remove_all)
    lonstring = strcompress(string(format="(f7.2)",gridlocs[1,iloc]),/remove_all)
    framename = string(format="(a,'_lat',a,'_lon',a)",sfields[ifld].name,$
      latstring,lonstring)
    text=string(format="('Frame ',i4,': Name=',a,' Title=',a)",$
      iframe,framename,title)
    myprint,text,/stdout,outfile=state.outfile
  endfor ; iloc
endfor ; ifld
end