function msis_utlat,field,slon,slev,lats,years,days,mtimes,zpoten,out
;
; Called from pltutlat:
;   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)
;
; zpoten = dblarr(nlats,ntime_file)
;
; Geopotential is in zpoten[nlats,ntimes]
; Output array out[nlats,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)
  nlats = n_elements(lats)
; print,'msis_utlat: field=',field,' slon=',slon,' slev=',slev,$
;   ' nlats=',nlats,' ntimes=',ntimes

  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_utlat: itime=',i4,' mtimes[*,itime]=',3i4,' sec=',i12)",$
;     itime,mtimes[*,itime],sec

    alt = fltarr(1)
    glon = slon

    found_field = 1
    if (field eq 'TN') then begin
      for j=0,nlats-1 do begin
        glat = lats[j]
        alt[0] = zpoten[j,itime]*1.e-5 ; cm to km
        msis00,date,sec,glat,glon,alt,t_alt,t_exo
; Average for 2002: ft7p=ft7a=179.4, ap=20.
;       msis00,date,sec,glat,glon,alt,t_alt,t_exo,nd,md,$
;         ft7p=179.4,ft7a=179.4,apd=20.
        out[j,itime] = t_alt
      endfor
    endif else begin
      for j=0,nlats-1 do begin
        glat = lats[j]
        alt[0] = zpoten[j,itime]*1.e-5 ; cm to km
;       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[j,itime] = nd.o2
          'O1' : out[j,itime] = nd.o
          'N'  : out[j,itime] = nd.n
          'N2' : out[j,itime] = nd.n2
          'AR' : out[j,itime] = nd.ar
          'HE' : out[j,itime] = nd.he
          'H'  : out[j,itime] = nd.h
          'RHO': out[j,itime] = md     ; mass density
          else: found_field = 0
        endcase
      endfor ; j=0,nlats-1
    endelse ; not TN
  endfor ; itime=0,ntimes-1

; if (found_field) then print,format=$
;   "('msis_utlat: slon=',f8.2,' slev=',f8.2,' ntimes=',i4,' nlats=',i4,' field ',a,': min,max=',2e12.4)",$
;   slon,slev,ntimes,nlats,field,min(out),max(out)

  return,found_field
end