function msis_utvert,field,loc,zplevs,years,days,mtimes,zpoten,out
;
; zpoten = dblarr(nzplev,ntime_file)
;
; Geopotential is in zpoten[0,0,nzplevs,ntimes]
; Output array out[nzplevs,ntimes]
;
; pro msis00, date, sec, glat, glon, alt, t_alt, t_exo, nd, md, $
;             ft7p=ft7p_usr, ft7a=ft7a_usr, apd=ap_usr, $
;             SI=si_yesno, Ap_3hr=ap3_yesno, HELP=h_flag
;  t_alt (output) -- scalar or array of same type as alt, with the
;                    temperature (K) at each of the input altitudes
;  t_exo (output) -- scalar, exospheric temperature (K)
;
;     Note: if only temperatures are desired, omit the following 
;           two output arguments, so that time is not wasted
;           computing the number densities.
;
;  nd (output)    -- number density, an array of the same length
;                    as the input argument alt, each element of
;                    which is a structure with seven scalar fields,
;                    one for each of seven atmospheric species:
;                    nd.h  -- atomic hydrogen
;                    nd.he -- helium
;                    nd.n  -- atomic nitrogen
;                    nd.o  -- atomic oxygen
;                    nd.n2 -- molecular nitrogen
;                    nd.o2 -- molecular oxygen
;                    nd.ar -- argon
;                    nd.ao -- anomalous oxygen
;
;  md (output)    -- scalar or array of same type as alt, with the
;                    total mass density of all species at each of
;                    the requested altitudes.
;
  ntimes = n_elements(years)
  nzplevs = n_elements(zplevs)
; print,'msis_utvert: field=',field,' loc=',loc,' ntimes=',ntimes,$
;   ' nzplevs=',nzplevs,' zplevs=',zplevs
  glat = loc[0]
  glon = loc[1]

  for itime=0,ntimes-1 do begin
    date = long(years[itime])*1000+long(days[itime])
    sec = float(mtimes[1,itime])*3600.+float(mtimes[2,itime])*60.

;   print,format="('msis_utvert: itime=',i4,' mtimes[*,itime]=',3i4,' sec=',i12)",$
;     itime,mtimes[*,itime],sec

    alt = fltarr(nzplevs)
;   alt[*] = zpoten[0,0,*,itime]*1.e-5 ; cm to km
    alt[*] = zpoten[*,itime]*1.e-5 ; cm to km

    found_field = 1

; pro msis00, date, sec, glat, glon, alt, t_alt, t_exo, nd, md, $
;             ft7p=ft7p_usr, ft7a=ft7a_usr, apd=ap_usr, $
;             SI=si_yesno, Ap_3hr=ap3_yesno, HELP=h_flag

    if (field eq 'TN') then begin
;     msis00,date,sec,glat,glon,alt,t_alt,t_exo
; Average for 2002: ft7p=ft7a=179.4, ap=20.

;     if (itime eq 0) then begin
;       print,'msis_utvert call msis00 for TN:' 
;       help,date & help,sec & help,glat & help,glon & help,alt
;     endif

      msis00,date,sec,glat,glon,alt,t_alt,t_exo,nd,md,$
        ft7p=179.4,ft7a=179.4,apd=20.
      out[*,itime] = t_alt
    endif else begin
      msis00,date,sec,glat,glon,alt,t_alt,t_exo,nd,md
;     msis00,date,sec,glat,glon,alt,t_alt,t_exo,nd,md,$
;       ft7p=75.,ft7a=75.,apd=4.
      case field of
        'O2' : out[*,itime] = nd.o2
        'O1' : out[*,itime] = nd.o
        'N'  : out[*,itime] = nd.n
        'N2' : out[*,itime] = nd.n2
        'AR' : out[*,itime] = nd.ar
        'HE' : out[*,itime] = nd.he
        'H'  : out[*,itime] = nd.h
        'RHO': out[*,itime] = md     ; mass density
        else: found_field = 0
      endcase
    endelse ; not TN
  endfor ; itime=0,ntimes-1

; if (found_field) then print,format=$
;   "('msis_utvert: lat,lon=',2f8.2,' ntimes=',i4,' nzplevs=',i4,' field ',a,': min,max=',2e12.4)",$
;   glat,glon,ntimes,nzplevs,field,min(out),max(out)

  return,found_field
end