;-----------------------------------------------------------------------
;
; Main is mtgcm_maven below.
;
;-----------------------------------------------------------------------
function ixfind,z,val
nz = n_elements(z)
;print,'ixfind: nz=',nz,' z=' & print,z,val
if val ge 0.5*(z[nz-1]+z[nz-2]) then return,nz-1
if val le 0.5*(z[0]+z[1]) then return,0
for i=1,nz-2 do begin
  if val le 0.5*(z[i]+z[i+1]) and val ge 0.5*(z[i]+z[i-1]) then $
    return,i
endfor
return,-1
end
;-----------------------------------------------------------------------
function fslt, slt,ut,lon,request,mag=mag,verbose=verbose
;
; (slt and ut always in decimal hours, glon in -180 -> 180)
;
v = 0 & if keyword_set(verbose) then v = 1
;
if request ne "getlon" then glon = lon
case request of
  "getlt": begin ; return slt
       fslt = ut + glon/15.
       if (fslt lt 0.) then fslt = fslt+24.
       if (fslt ge 24.) then fslt = fslt-24.
       if v gt 0 then print,'fslt getlt: mag=',mag,' ut=',ut,$
         ' glon=',glon,' new lt=',fslt
      end
  "getut": begin ; return ut
       fslt = slt - glon/15.
       if (fslt lt 0.) then fslt = fslt+24.
       if (fslt gt 24.) then fslt = fslt-24.
       if v gt 0 then print,'fslt getut: mag=',mag,' slt=',slt,$
         ' glon=',glon,' new ut=',fslt
     end
  "getlon": begin ; return longitude
       fslt = (slt-ut) * 15.
       if keyword_set(mag) then fslt = fslt+71.
       if (fslt ge 180.) then fslt = fslt-360.
       if (fslt lt -180.) then fslt = fslt+360.
       if v gt 0 then print,'fslt getlon: mag=',mag,' ut=',ut,$
         ' slt=',slt,' new lon=',fslt
     end
  else: begin
;      print,'>>> fslt: unknown request=',request
       fslt = 0.
     end
endcase
return,fslt
end
;-----------------------------------------------------------------------
pro read_maven,state,file
;
; Read MAVEN satellite ascii data file, and scatter-plot density.
; At end, update state for main.
;
; 4/16/12 btf: read new file w/ different format
;
  print,' ' & print,'Enter read_maven'
  openr,lu,file,/get_lun
  print,'Opened file ',file,' lu=',lu
;
; Read file name on 2nd line:
; fname = ''
; readf,lu,format="(/,a)",fname 
; print,'fname = ',fname
;
; Skip rest of header:
; line = ''
; for i=1,5 do readf,lu,line
;
; Read and echo 4-line header:
  line = ''
  for i=1,4 do begin
    readf,lu,line
    print,line
  endfor
;
  mxdata = 2000
  year = intarr(mxdata) ; 4-digit year
  day  = lonarr(mxdata) ; julian day
  ut   = fltarr(mxdata) ; ut (UTCG)
  slt  = fltarr(mxdata) ; local time
  lat  = fltarr(mxdata) ; Latitude (deg) (Mars body-fixed)
  lon  = fltarr(mxdata) ; Longitude (deg) (Mars body-fixed)
  slat = fltarr(mxdata) ; Latitude (deg) (Subsolar point)
  slon = fltarr(mxdata) ; Longitude (deg) (Subsolar point)
  den  = fltarr(mxdata) ; MAVEN_Mars_Gram (kg/km^3)
  alt  = fltarr(mxdata) ; MARS_Altitude (km)
  date = strarr(mxdata) ; MAVEN date string

  months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
  month = '' ; 3-chars (must be one of months above)
  dd=0 & yy=0 & hh=0 & mm=0 & ss=0.0 
  rdlat=0. & rdlon=0. & rdden=0. & rdalt=0.
  rdslat=0. & rdslon=0.
  norbit = 1 ; number of orbits (current orbit while reading)
  altmin = 1.e10 ; minimum altitude of current orbit
  aveslt = 0.
;
; Read data to EOF:

  ndata = 0
  while ~ eof(lu) do begin
    readf,lu,line
;
; Skip periodic 3-line labeling:
;   if (strtrim(line) eq '') then begin
;     for i=1,2 do readf,lu,line
;     print,'Skipped 3-line labeling at ndata ',ndata
;     continue 
;   endif
;
; Check for max data:
    ndata = ndata+1
    if (ndata gt mxdata) then begin
      print,'>>> ndata > mxdata =',ndata,' Please increase mxdata.'
      free_lun,lu
      return
    endif
;
; Read data from line string:
    date[ndata-1] = strmid(line,0,24)
;   reads,line,format="(i2,1x,a3,i5,i3,1x,i2,1x,f6.3,f35.3,f36.3,f34.6,f22.6)",$
;     dd,month,yy,hh,mm,ss,rdlat,rdlon,rdden,rdalt
    reads,line,format=$
    "(i2,1x,a3,i5,i3,1x,i2,1x,f7.4,4x,f11.3, f13.3, f0,  f12.3,f11.3,14x,f12.3,f11.3)",$
      dd,month,yy,hh,   mm,   ss,     rdslat,rdslon,rdden, rdalt1,rdalt2,rdlat,rdlon
;   print,' ndata=',ndata,' dd ',dd,' month=',month,' yy=',yy,' hh=',hh,' mm=',mm,' ss=',ss,$
;         ' rdslat=',rdslat,' rdslon=',rdslon,' rdden=',rdden,' rdalt=',rdalt,' rdlat,lon=',$
;         rdlat,',',rdlon 
    imon=0
    for i=0,11 do if (month eq months[i]) then imon = i+1
    if (imon eq 0) then begin
      print,'>>> Could not find month ',month,' in months array'
      return
    endif
    year[ndata-1]=yy
    day[ndata-1]=julday(imon,dd,-4713)
    ut[ndata-1] = float(hh)+float(mm)/60.+float(ss)/3600.
    lat[ndata-1]=rdlat
    lon[ndata-1]=rdlon
    slat[ndata-1]=rdslat
    slon[ndata-1]=rdslon
    den[ndata-1]=rdden
;
; Get local time at given longitude and ut:
    dummy = 0.
    slt[ndata-1] = fslt(dummy,ut[ndata-1],slon[ndata-1],"getlt")
;
; rdalt1 = PdeticAlt, rdalt2 = MolaAlt
; Use PdeticAlt for locating beginning of an orbit, and use MolaAlt 
; for finding periapsis for each orbit, and making plots and interpolating
; mtgcm data.
;
    alt[ndata-1]=rdalt2
;
; Check for minimum Mola altitude:
    if (rdalt2 lt altmin) then begin
      altmin = rdalt2
      ndatamin = ndata-1 ; data index at periapsis of current orbit
    endif
;
; End current orbit (data file begins w/ PdeticAlt ==500, so ignore that one):
    if (rdalt1 eq 500. and ndata ne 1) then begin
; slt at periapsis, assuming constant ut:
      sltperi = fslt(dummy,state.maven_utconstant,slon[ndatamin],"getlt") 
      aveslt = aveslt + sltperi
;     print,'End orbit ',norbit,' altmin=',altmin,' ndatamin=',ndatamin,$
;       ' slonperi=',slon[ndatamin],' ut=',state.maven_utconstant,' sltperi=',$ 
;       sltperi,' aveslt (sum)=',aveslt
; Update average slt over all orbits at periapsis:
      altmin = 1.e10        ; begin new search for periapsis
      norbit = norbit+1     ; update orbit number
    endif
  endwhile ; not eof
;
; Update aveslt for last orbit:
  sltperi = fslt(dummy,state.maven_utconstant,slon[ndatamin],"getlt")
  aveslt = aveslt + sltperi
; print,'End orbit ',norbit,' altmin=',altmin,' ndatamin=',ndatamin,$
;   ' slonperi=',slon[ndatamin],' ut=',state.maven_utconstant,' sltperi=',$ 
;   sltperi,' aveslt (sum)=',aveslt
;
; Calculate final average over orbits of slt at periapsis:
  aveslt = aveslt / float(norbit)
  print,'End read: norbit=',norbit,' aveslt=',aveslt
  aveslt = float(round(aveslt))
  print,'rounded aveslt=',aveslt
;
; Truncate arrays to number of data points: 
  year = year[0:ndata-1]
  day  = day [0:ndata-1]
  ut   = ut  [0:ndata-1]
  lat  = lat [0:ndata-1]
  lon  = lon [0:ndata-1]
  den  = den [0:ndata-1]
  alt  = alt [0:ndata-1]
  slt  = slt [0:ndata-1]
;
; Print to stdout:
; for i=0,ndata-1 do begin
;   print,format=$
;     "('i=',i6,' ndata=',i6,' year=',i4,' day=',i4,' ut=',f8.3,' slt=',f8.3,' slat=',f8.3,' slon=',f9.3,' den=',e12.4,' alt=',f12.6,' lat=',f8.3,' lon=',f8.3)",$
;     i,ndata,year[i],day[i],ut[i],slt[i],slat[i],slon[i],den[i],alt[i],lat[i],lon[i]
; endfor
;
; Close file:
  free_lun,lu

  if (not state.ps) then begin
    window,xsize=700,ysize=600,retain=2
    device,decomposed=0
  endif
;
; Plot points at lat vs altitude, and contour on irregular grid,
; with points colored according to log10 density:
;
  denmin = min(den) & denmax = max(den)
  maven_den = den   ; save density w/o log10
  den = alog10(den) ; convert den to log10(den)
; fmin = min(den) & fmax = max(den)
  fmin = denmin & fmax = denmax
  title = 'MAVEN Density: min,max=' + strcompress(string(denmin)) + $
           strcompress(string(denmax)) + ' (kg/km^3)'
  xll = .15 & yll = .25
  xur = .90 & yur = .90
  xmid = .5*(xll+xur)
  pos = [xll,yll,xur,yur]
  chsize = 1.1
;
; psym = 1(+), 2(*), 3(.), 4(diamond), 5(triangle), 6(square), 7(x)
  psym = 1
  symsize = 1.2
; contour,den,lat,alt,/irregular,nlevels=15,/follow,/nodata,position=pos,$
;   xtitle='Latitude (deg)',ytitle='Altitude (km)',title=title,charsize=chsize
; plots,lat,alt,psym=psym,symsize=symsize,color=bytscl(den,top=!d.n_colors-1)

  cblabels = [strcompress(string(fmin,format="(g10.4)"),/remove_all),$
    strcompress(string(0.5*(fmin+fmax),format="(g10.4)"),/remove_all),$
    strcompress(string(fmax,format="(g10.4)"),/remove_all)]
  barw = .90*(xur-xll) & barh = .03
  xbar = xmid-.5*barw & ybar = yll-barh-.13
; clrbar,barw,barh,xbar,ybar,cblabels,charsize=chsize,botcolor=1,$
;   topcolor=!d.table_size-2
;
; Update state with satellite data:
  state.maven_ndata   = ndata
  *state.maven_date   = date
  *state.maven_year   = year
  *state.maven_day    = day
  *state.maven_ut     = ut
  *state.maven_lat    = lat
  *state.maven_lon    = lon
  *state.maven_den    = maven_den ; w/o log10
  *state.maven_alt    = alt
  *state.maven_slt    = slt
  state.maven_aveslt = aveslt
end ; read_maven
;-----------------------------------------------------------------------
pro read_mtgcm,state,files
print,' ' & print,'Enter read_mtgcm'
 nfiles = n_elements(files)
;
; Get total ntimes and mtimes(3,ntimes):
;
 ntimes = 0     ; total ntimes for all files
 mxtimes = 100
 mtimes = intarr(3,mxtimes) ; model times for all files
 for ifile=0,nfiles-1 do begin
   ncid = ncdf_open(files[ifile])
   id = ncdf_dimid(ncid,'time')
   ncdf_diminq,ncid,id,name,ntime 
   ncdf_varget,ncid,'mtime',mtime
   for i=0,ntime-1 do begin
     ii = ntimes + i
     mtimes[*,ii] = mtime[*,i]
   endfor
   ntimes = ntimes+ntime
   if (ntimes gt mxtimes) then begin
     print,'>>> need to increase mxtimes'
     return
   endif
   ncdf_close,ncid
 endfor
 mtimes = mtimes[*,0:ntimes-1]
;print,'nfiles=',nfiles,' ntimes=',ntimes
;for i=0,ntimes-1 do print,'i=',i,' mtime=',mtimes[*,i]
 ntimes_total = ntimes
;
; Calculated fields not on history: CO2, RHO
; (all others are on history)
;
 fields = ['TN' ,'UN' ,'VN' ,'O1' ,'CO' ,'O2' ,'N2' ,'AR' ,'NE' ,'Z',  $
           'CO2','RHO']
 hist   = [ 1   ,1    ,1    ,1    ,1    ,1    ,1    ,1    ,1    ,1    ,$
            0   ,0]
 loglin = ['LIN','LIN','LIN','LOG','LOG','LOG','LOG','LOG','LOG','LIN',$
           'LOG','LOG']
;
; Units *after* pro mtgcm_units:
 units  = ['deg K','m/s','m/s','cm3','cm3','cm3','cm3','cm3','cm3','km',$
           'cm3','kg/km3']
 nf = n_elements(fields)
;
; File loop to read data:
; 
 ntimes = 0     ; total ntimes for all files
 for ifile=0,nfiles-1 do begin
   ncid = ncdf_open(files[ifile])
   print,format="('Opened file ',a,' (',i3,' of ',i3,')')",$
     files[ifile],ifile,nfiles
   id = ncdf_dimid(ncid,'time')
   ncdf_diminq,ncid,id,name,ntime 
;
; Print global file attributes:
;
;  fstruct = ncdf_inquire(ncid)
;  print,'Global File Attributes:'
;  for i=0,fstruct.ngatts-1 do begin
;    attname = ncdf_attname(ncid,i,/global)
;    ncdf_attget,ncid,attname,att,/global
;    att = string(att)
;    print,'  ',attname,' = ',att
;  endfor 
;
; Coordinates:
;
   ncdf_varget, ncid, 'lon', lon
   ncdf_varget, ncid, 'lat', lat
   ncdf_varget, ncid, 'lev', lev
;
; Grid size:
;
   nlon = n_elements(lon)
   nlat = n_elements(lat)
   nlev = n_elements(lev)
;  print,'ntime=',ntime,' nlon=',nlon,' nlat=',nlat,' nlev=',nlev
 
   if (ifile eq 0) then begin
     f = fltarr(nlon,nlat,nlev,ntimes_total,nf)
   endif

   for ifld=0,nf-1 do begin
     if (hist[ifld] ne 1) then continue ; skip if not on history
     ncdf_varget,ncid,fields[ifld],fld
     if (ntime eq 1) then begin
       ii = ntimes
       f[*,*,*,ii,ifld] = fld
     endif else begin
       for i=0,ntime-1 do begin
         ii = ntimes + i
;        print,format="('Define field ',a,' i=',i3,' ii=',i3)",$
;          fields[ifld],i,ii
         f[*,*,*,ii,ifld] = fld[*,*,*,i]
       endfor
     endelse
   endfor ; field loop
   ntimes = ntimes+ntime
   ncdf_close, ncid
 endfor ; file loop
;
; Save model ut:
 ut = fltarr(ntimes)
 for i=0,ntimes-1 do ut[i] = mtimes[1,i]+mtimes[2,i]/60. 

; Update state:
   state.mtgcm_ntimes = ntimes
  *state.mtgcm_mtimes = mtimes
  *state.mtgcm_ut = ut
  *state.mtgcm_lat = lat
  *state.mtgcm_lon = lon
  *state.mtgcm_lev = lev
   state.mtgcm_nflds = nf
  *state.mtgcm_fnames = fields
  *state.mtgcm_loglin = loglin
  *state.mtgcm_hist = hist
  *state.mtgcm_units = units
  *state.mtgcm_f = f

end ; read_mtgcm
;-----------------------------------------------------------------------
pro mtgcm_units,state
;
; Do units conversion (winds and densities), and calculate co2 and rho
; This is called after read_mtgcm.
;
  f = *state.mtgcm_f ; (nlon,nlat,nlev,ntimes,nf)
  fnames = *state.mtgcm_fnames
  nflds = state.mtgcm_nflds
  ixco = -1 & ixo1 = -1 & ixn2 = -1 & ixtn = -1
  ixo2 = -1 & ixar = -1 & ixrho = -1 & ixco2 = -1
  for i=0,nflds-1 do begin
    if (fnames[i] eq 'CO') then ixco = i
    if (fnames[i] eq 'O1') then ixo1 = i
    if (fnames[i] eq 'O2') then ixo2 = i
    if (fnames[i] eq 'N2') then ixn2 = i
    if (fnames[i] eq 'TN') then ixtn = i
    if (fnames[i] eq 'AR') then ixar = i
    if (fnames[i] eq 'CO2') then ixco2 = i
    if (fnames[i] eq 'RHO') then ixrho = i
  endfor
  if (ixco eq -1) then print,'>>> mtgcm_units: Need CO'
  if (ixo1 eq -1) then print,'>>> mtgcm_units: Need O1'
  if (ixn2 eq -1) then print,'>>> mtgcm_units: Need N2'
  if (ixtn eq -1) then print,'>>> mtgcm_units: Need TN'
  if (ixo2 eq -1) then print,'>>> mtgcm_units: Could not find O2'
  if (ixar eq -1) then print,'>>> mtgcm_units: Could not find AR'
  if (ixco2 eq -1) then print,'>>> mtgcm_units: Could not find CO2'
  if (ixrho eq -1) then print,'>>> mtgcm_units: Could not find RHO'

  nlon = n_elements(*state.mtgcm_lon)
  nlat = n_elements(*state.mtgcm_lat)
  nlev = n_elements(*state.mtgcm_lev)
  ntime = state.mtgcm_ntimes
  levs = *state.mtgcm_lev

  rho  = fltarr(nlon,nlat,nlev)
  co2  = fltarr(nlon,nlat,nlev)
  barm = fltarr(nlon,nlat,nlev)
  pkt  = fltarr(nlon,nlat,nlev)
  p0 = 1.2e-3
  boltz = 1.3805e-16
  dz = levs[1]-levs[0] ; assume monotonic
;
  for it=0,ntime-1 do begin
;
; co2 = 1-o-co-n2 (insure co2 > 0.)
    o1 = f[*,*,*,it,ixo1]
    o2 = f[*,*,*,it,ixo2]
    co = f[*,*,*,it,ixco]
    n2 = f[*,*,*,it,ixn2]
    tn = f[*,*,*,it,ixtn]
    ar = f[*,*,*,it,ixar]
    co2 = 1.-o1-co-n2
    co2[where(co2 < 0.)] = .00001
;
; barm = 1/(o/16+(co+n2)/28+co2/44)
    barm = 1./(o1/16.+(co+n2)/28.+co2/44.)
    for k=0,nlev-1 do begin
      pkt[*,*,k] = p0*exp(-(levs[0]+float(k)*dz))/(boltz*tn[*,*,k])
    endfor
;
; Calculate rho in gm/cm3, and convert to kg/km3:
    k = 1.66e-24
    rho = (co2*pkt*barm*k + n2*pkt*barm*k + $
           co *pkt*barm*k + o1*pkt*barm*k) * 1.e12
;
; Convert densities: cm3 = mmr * pkt * barm / w
    o1 = o1 * pkt * barm / 16.
    o2 = o2 * pkt * barm / 32.
    co = co * pkt * barm / 28.
    n2 = n2 * pkt * barm / 28.
    ar = ar * pkt * barm / 40.
    co2 = co2 * pkt * barm / 44.

    f[*,*,*,it,ixo1] = o1[*,*,*]
    f[*,*,*,it,ixo2] = o2[*,*,*]
    f[*,*,*,it,ixco] = co[*,*,*]
    f[*,*,*,it,ixn2] = n2[*,*,*]
    f[*,*,*,it,ixar] = ar[*,*,*]
    f[*,*,*,it,ixco2] = co2[*,*,*]
    f[*,*,*,it,ixrho] = rho[*,*,*]
  endfor ; it=0,ntime-1
;
; cm/s -> m/s for neutral winds:
  for i=0,nflds-1 do begin
    if (fnames[i] eq 'UN' or fnames[i] eq 'VN') then begin
      f[*,*,*,*,i] = f[*,*,*,*,i] * .01
      print,'Converted ',fnames[i],' from m/s to cm/s: min,max=',$
        min(f[*,*,*,*,i]),max(f[*,*,*,*,i])
    endif
    if (fnames[i] eq 'Z') then begin
      f[*,*,*,*,i] = f[*,*,*,*,i] * 1.e-5
      print,'Converted ',fnames[i],' from cm to km: min,max=',$
        min(f[*,*,*,*,i]),max(f[*,*,*,*,i])
    endif
  endfor
;
; Update state:
  *state.mtgcm_f = f
end ; mtgcm_units
;-----------------------------------------------------------------------
pro make_sltslices,state,slt
;
; Make model local time slices (nlat,nlev,ntime,nflds) at requested slt.
; Also return longitude coords used in 3rd dimension of the slices
; Model fields are defined in mtgcm_f[nlon,nlat,nlev,ntimes,nf]
;
ntime = state.mtgcm_ntimes
nlat = n_elements(*state.mtgcm_lat)
nlon = n_elements(*state.mtgcm_lon)
nlev = n_elements(*state.mtgcm_lev)
nfld = n_elements(*state.mtgcm_fnames)
fnames = *state.mtgcm_fnames
hist = *state.mtgcm_hist
f = *state.mtgcm_f
ut = *state.mtgcm_ut
lons = *state.mtgcm_lon
;print,'make_sltslices: ntime=',ntime,' nlat=',nlat,' nlon=',nlon,' nfld=',nfld
;print,'make_sltslices: ntime=',ntime,' nlon=',nlon,' mtgcm lons=' & print,lons
sltlons = fltarr(ntime) ; longitudes at requested slt
dummy = 0.
;
; Use average periapsis slt and model ut to find model longitudes:

sltslices = fltarr(nlat,nlev,ntime,nfld)
for ifld=0,nfld-1 do begin
  for it=0,ntime-2 do begin
    sltlons[it] = fslt(slt,ut[it],dummy,"getlon") ; longitude at slt, ut[it]
    ii = ixfind(lons,sltlons[it])                 ; find index to sltlon
    if (ii eq -1) then begin
      print,'>>> make_sltslices: could not find longitude ',sltlons[it],$
        ' in model longitudes:' & print,'    model longitudes=',lons
      sltslices[*,*,it,ifld] = 0.
      continue
    endif
    sltslices[*,*,it,ifld] = f[ii,*,*,it,ifld]
  endfor
  sltlons[ntime-1] = 180.
endfor
print,'make_sltslices: ntime=',ntime,' unsorted sltlons=' & print,sltlons
;
; Sort sltlons into ascending order and reorder model data in same way:
 indices = sort(sltlons)
 sltlons = sltlons[indices]
 print,'make_sltslices: ntime=',ntime,' sorted sltlons=' & print,sltlons
 sltslices = sltslices[*,*,indices,*] ; sort data same as longitudes
;
; Set data for last time: note this is *after* sorting data into ascending
; longitudes, so 0 and ntime-1 refer to -/+180 longitudes, not 0 and 24 ut.
;
 sltslices[*,*,ntime-1,*] = sltslices[*,*,0,*]
;
; Return slices and longitudes in state:
;
*state.mtgcm_sltslices = sltslices
*state.mtgcm_sltlons   = sltlons

end ; make_sltslices
;-----------------------------------------------------------------------
pro plot_sltslices,state,slt
;
; Plot model slt slices at a few arbitrary zp levels:
; sltslices(nlat,nlev,ntime,nfld), sltlons(ntime)
;
sltslices = *state.mtgcm_sltslices
sltlons   = *state.mtgcm_sltlons
nflds     = state.mtgcm_nflds
;nlons     = state.mtgcm_ntimes+1
 nlons     = n_elements(sltlons)
nlats     = n_elements(*state.mtgcm_lat)
nlevs     = n_elements(*state.mtgcm_lev)
fnames    = *state.mtgcm_fnames
hist      = *state.mtgcm_hist
units     = *state.mtgcm_units
loglin    = *state.mtgcm_loglin
;orbit_lat = *state.maven_lat
;orbit_lon = *state.maven_lon

lons = sltlons

print,'plot_sltslices: nlons=',nlons,' lons=' & print,lons

lats = *state.mtgcm_lat
levs = *state.mtgcm_lev
sltslice = fltarr(nlons,nlats)

xll = .15 & yll = .35
xur = .95 & yur = .85
xmid = .5*(xll+xur)
conpos = [xll,yll,xur,yur]
mappos = [xll,yll,xur,yur]
plat = 0. & plon=0. & rot=0.
xticlabs = ['-180','-120','-60','0','60','120','180']
yticlabs = ['-90','-60','-30','0','30','60','90']
xtitle = 'LONGITUDE (Constant Local Time)'
ytitle = 'LATITUDE (deg)'
chsize = 1.2
barw = .90*(xur-xll) & barh = .03
xbar = xmid-.5*barw & ybar = yll-barh-.15
;
; Get average global geopotential at each level for plot title:
aveht = fltarr(nlevs)
for ifld=0,nflds-1 do begin
  if (fnames[ifld] eq 'Z') then begin
    for k=0,nlevs-1 do begin 
      for i=0,nlons-1 do begin
        for j=0,nlats-1 do begin
          sltslice[i,j] = sltslices[j,k,i,ifld]
        endfor
      endfor
      aveht[k] = mean(sltslice)
    endfor
  endif
endfor
print,' '
delz = 10 ; arbitrarilly plot every 10th level
for ifld=0,nflds-1 do begin
  for k=0,nlevs-1,delz do begin 
    for i=0,nlons-1 do begin
      for j=0,nlats-1 do begin
        sltslice[i,j] = sltslices[j,k,i,ifld]
      endfor
    endfor
;
; Optionally plot log10 of the field:
    plot_log10 = 1 ; (set to 0 if you do *not* want to plot log10)
    if (loglin[ifld] eq 'LOG' and plot_log10 eq 1) then $
      sltslice = alog10(sltslice)

    if (not state.ps) then window,1,xsize=700,ysize=600,retain=2
    map_set,plat,plon,rot,/cylindrical,position=mappos,/noborder
    fmin = min(sltslice) & fmax = max(sltslice)
    image = bytscl(sltslice,top=!d.table_size-3,min=fmin,max=fmax)+1
    image = map_patch(image,lons,lats,xstart=startx,ystart=starty,xsize=xsize,ysize=ysize)
    tv,image,startx,starty,xsize=xsize,ysize=ysize
    map_continents,color=0
;   map_grid,/label,color=0
    cblabels = [strcompress(string(fmin,format="(g10.4)"),/remove_all),$
      strcompress(string(0.5*(fmin+fmax),format="(g10.4)"),/remove_all),$
      strcompress(string(fmax,format="(g10.4)"),/remove_all)]
    if (loglin[ifld] eq 'LOG' and plot_log10 eq 1) then $
      title = 'MTGCM ' + 'log10 ' + fnames[ifld] + ' ('+units[ifld]+')'+ $
          ' SLT=' + strcompress(string(format="(f5.2)",slt)) + $
          ' Zp=' + strcompress(string(format="(f6.2)",levs[k])) + $
          ' Ave Ht=' + strcompress(string(format="(f6.2)",aveht[k])) $
    else $
      title = 'MTGCM ' + fnames[ifld] + ' ('+units[ifld]+')'+ $
          ' SLT=' + strcompress(string(format="(f5.2)",slt)) + $
          ' Zp=' + strcompress(string(format="(f6.2)",levs[k])) + $
          ' Ave Ht=' + strcompress(string(format="(f6.2)",aveht[k]))

    print,'Plot: ',title,' min,max=',fmin,fmax
    contour,sltslice,lons,lats,/noerase,position=mappos,/nodata,$
      xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
      xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
      yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs,title=title
    clrbar,barw,barh,xbar,ybar,cblabels,charsize=chsize,botcolor=1,$
      topcolor=!d.table_size-2
;   map_set,plat,plon,rot,/cylindrical,position=mappos,/noborder
;   plots,orbit_lat,orbit_lon,psym=1 ; not working
    if (not state.ps) then cursor,x,y,/wait,/up

  endfor ; k=0,nlevs-1
endfor ; ifld=0,nflds-1

end ; plot_sltslices
;-----------------------------------------------------------------------
function bracket,xa,na,x,ix0,ix1
;
; Bracket x between 2 indices of vector xa(na), returning
; ix0 as the first index, and ix1 as the second.
; Return 0 if no error, non-zero with warning msg otherwise.
; In case of error, ix0 == ix1 == -1
;
error = 0
ix0 = -1 & ix1 = -1
for i=1,na-1 do begin
  if (xa[i-1] le x and xa[i] ge x) then begin
    ix1=i & ix0=i-1
    break
  endif
endfor 
if (ix0 eq -1 or ix1 eq -1) then begin
; print,'>>> WARNING bracket: could not find x=',x,' in xa:',$
;   ' na=',na,' xa min,max=',min(xa),max(xa)
  error = 1
endif
return,error
end ; function bracket
;-----------------------------------------------------------------------
pro plot_mtgcm_interp,state
;
; Make lat vs alt scatter plots of mtgcm fields interpolated
; to maven sat track.
;
  mtgcm_interp = *state.mtgcm_interp
  nflds = state.mtgcm_nflds
  fnames = *state.mtgcm_fnames
  loglin = *state.mtgcm_loglin
  hist   = *state.mtgcm_hist
  units  = *state.mtgcm_units
  data_lat = *state.maven_lat
  data_alt = *state.maven_alt

  xll = .15 & yll = .25
  xur = .90 & yur = .90
  xmid = .5*(xll+xur)
  pos = [xll,yll,xur,yur]
  chsize = 1.1
; psym = 1(+), 2(*), 3(.), 4(diamond), 5(triangle), 6(square), 7(x)
  psym = 1
  symsize = 1.2

  for i=0,nflds-1 do begin
    if (fnames[i] eq 'Z') then continue 
    f = mtgcm_interp[*,i]
    lat = data_lat[where(f ne 1.e36)]
    alt = data_alt[where(f ne 1.e36)]
    f = f[where(f ne 1.e36)]
    if (loglin[i] eq 'LOG') then begin
      f = alog10(f)
      fmin = min(f) & fmax = max(f)
      print,format="('plot_mtgcm_interp: log10 ',a,' min,max=',2e12.4)",$
        fnames[i],fmin,fmax
      title = 'MTGCM '+ 'log10 ' + fnames[i] + ' (' + units[i] + ')' + ' min,max=' + $
               strcompress(string(fmin)) + strcompress(string(fmax))
    endif else begin
      fmin = min(f) & fmax = max(f)
      print,format="('plot_mtgcm_interp: ',a,' min,max=',2e12.4)",$
        fnames[i],fmin,fmax
      title = 'MTGCM ' + fnames[i] + ' (' + units[i] + ')' + ' min,max=' + $
               strcompress(string(fmin)) + strcompress(string(fmax))
    endelse

    contour,f,lat,alt,/irregular,nlevels=15,/follow,/nodata,position=pos,$
      xtitle='Latitude (deg)',ytitle='Altitude (km)',title=title,charsize=chsize
    plots,lat,alt,psym=2,symsize=symsize,$
      color=bytscl(f,top=!d.table_size-3,min=fmin,max=fmax)+1

    cblabels = [strcompress(string(fmin,format="(g10.4)"),/remove_all),$
      strcompress(string(0.5*(fmin+fmax),format="(g10.4)"),/remove_all),$
      strcompress(string(fmax,format="(g10.4)"),/remove_all)]
    barw = .90*(xur-xll) & barh = .03
    xbar = xmid-.5*barw & ybar = yll-barh-.13
    clrbar,barw,barh,xbar,ybar,cblabels,charsize=chsize,botcolor=1,$
      topcolor=!d.table_size-2

    if (not state.ps) then cursor,x,y,/wait,/up
  endfor ; i=0,nflds-1

end ; plot_mtgcm_interp
;-----------------------------------------------------------------------
pro mtgcm_maven
;
; Define state info:
;
state = { $
;
; Maven state data:
;
  maven_ndata:0,                       $ ; number of data points in maven file 
  maven_date:ptr_new(/allocate_heap),  $ ; maven date string
  maven_year:ptr_new(/allocate_heap),  $ ; pointer to 4-digit integer year(ndata)
  maven_day:ptr_new(/allocate_heap),   $ ; pointer to integer day(ndata)
  maven_ut:ptr_new(/allocate_heap),    $ ; pointer to float ut(ndata)
  maven_lat:ptr_new(/allocate_heap),   $ ; pointer to float latitude(ndata)
  maven_lon:ptr_new(/allocate_heap),   $ ; pointer to float longitude(ndata)
  maven_den:ptr_new(/allocate_heap),   $ ; pointer to float est density(ndata)
  maven_alt:ptr_new(/allocate_heap),   $ ; pointer to float altitude(ndata)
  maven_slt:ptr_new(/allocate_heap),   $ ; pointer to float local time(ndata)
  maven_aveslt:0.                  ,   $ ; slt at periapsis, averaged over all orbits
  maven_utconstant:12.             ,   $ ; constant ut assumed for calc of ave slt
;
; Model state data:
;
  mtgcm_ntimes:0,                      $ ; number of model times
  mtgcm_mtimes:ptr_new(/allocate_heap),$ ; model times(3,ntimes)
  mtgcm_ut:ptr_new(/allocate_heap),    $ ; model ut(ntimes)
  mtgcm_lat:ptr_new(/allocate_heap),   $ ; model latitudes
  mtgcm_lon:ptr_new(/allocate_heap),   $ ; model longitudes
  mtgcm_lev:ptr_new(/allocate_heap),   $ ; model levels
  mtgcm_nflds:0,                       $ ; number of model fields
  mtgcm_fnames:ptr_new(/allocate_heap),$ ; model field names
  mtgcm_loglin:ptr_new(/allocate_heap),$ ; do log (LOG) or linear (LIN) interpolation for each field
  mtgcm_hist:ptr_new(/allocate_heap),  $ ; 1 if field is on history, 0 otherwise (i.e., calculated)
  mtgcm_units:ptr_new(/allocate_heap), $ ; field units (see conversions in pro mtgcm_units)
  mtgcm_f:ptr_new(/allocate_heap),     $ ; model data (nlon,nlat,nlev,ntimes,nf)
  mtgcm_sltslices:ptr_new(/allocate_heap), $ ; slices at slt=0 (nlat,nlev,ntime,nfld)
  mtgcm_sltlons:ptr_new(/allocate_heap),   $ ; longitudes in 3rd dim of sltslices 
  mtgcm_interp:ptr_new(/allocate_heap),    $ ; interpolated mtgcm (ndata,nflds)
;
; Set ps:0 for plots to come to the screen, or set ps:1 to make
;   a ps file "mtgcm_maven.ps" containing plots.
;
  ps:1                                 $ ; 1 if plotting to ps, 0 otherwise
}
stateptr = ptr_new(state)
state = *stateptr
;
; setclrtab can be called w/ any of these color tables:
; (if called w/ -1, will get custom blue-white-red)
; (calling loadct w/o arg will print above table and prompt for table)
;
; 0-        B-W LINEAR   14-             STEPS   28-         Hardcandy
; 1-        BLUE/WHITE   15-     STERN SPECIAL   29-            Nature
; 2-   GRN-RED-BLU-WHT   16-              Haze   30-             Ocean
; 3-   RED TEMPERATURE   17- Blue - Pastel - R   31-        Peppermint
; 4- BLUE/GREEN/RED/YE   18-           Pastels   32-            Plasma
; 5-      STD GAMMA-II   19- Hue Sat Lightness   33-          Blue-Red
; 6-             PRISM   20- Hue Sat Lightness   34-           Rainbow
; 7-        RED-PURPLE   21-   Hue Sat Value 1   35-        Blue Waves
; 8- GREEN/WHITE LINEA   22-   Hue Sat Value 2   36-           Volcano
; 9- GRN/WHT EXPONENTI   23- Purple-Red + Stri   37-             Waves
;10-        GREEN-PINK   24-             Beach   38-         Rainbow18
;11-          BLUE-RED   25-         Mac Style   39-   Rainbow + white
;12-          16 LEVEL   26-             Eos A   40-   Rainbow + black
;13-           RAINBOW   27-             Eos B

;colortab = -1 ; custom blue-white-red
;colortab = 11 ; this does not look like BLUE-RED to me
;colortab = 17 ; 
 colortab = 25 ; 
 
if (state.ps) then begin
  psfile = 'mtgcm_maven.ps'
  set_plot,'PS'
  device, filename=psfile, /color, bits_per_pixel=8
  setclrtab,colortab,/ps 
endif else begin
  set_plot,'X'
  setclrtab,colortab 
endelse
;
; Read MAVEN orbit text file, obtaining satellite attitude and density data: 
;
maven_file = 'DD3_SolarLatLong_Mar29toApr3__LT500km.Demcak.txt'

read_maven,state,maven_file

print,' ' & print,'After read_maven:'
print,'maven_ndata=',state.maven_ndata
print,'maven_day min,max=',min(*state.maven_day),max(*state.maven_day)
print,'maven_ut  min,max=',min(*state.maven_ut ),max(*state.maven_ut)
print,'maven_lat min,max=',min(*state.maven_lat),max(*state.maven_lat)
print,'maven_lon min,max=',min(*state.maven_lon),max(*state.maven_lon)
print,'maven_slt min,max=',min(*state.maven_slt),max(*state.maven_slt)
print,'maven_alt min,max=',min(*state.maven_alt),max(*state.maven_alt)
print,'maven_den min,max=',min(*state.maven_den),max(*state.maven_den)
;
; Read MTGCM history files:
;
;mtgcm_files = $
;['/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D1.FISM.TAU0p5.damp7.5.nc',$
;'/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D2.FISM.TAU0p5.damp7.5.nc',$
;'/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D3.FISM.TAU0p5.damp7.5.nc',$
;'/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D4.FISM.TAU0p5.damp7.5.nc']
;mtgcm_files = $
;['/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D1.FISM.TAU0p5.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D2.FISM.TAU0p5.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D3.FISM.TAU0p5.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D4.FISM.TAU0p5.damp7.5.nc']
;
;mtgcm_files = $
;['/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D1.FISM.TES1.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D2.FISM.TES1.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D3.FISM.TES1.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D4.FISM.TES1.damp7.5.nc']
;mtgcm_files = $
;['/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D1.FISM.TES1.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D2.FISM.TES1.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D3.FISM.TES1.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D4.FISM.TES1.damp7.5.nc']
;
;mtgcm_files = $
;['/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D1.FISM.TES2.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D2.FISM.TES2.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D3.FISM.TES2.damp7.5.nc',$
; '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F070L330f.D4.FISM.TES2.damp7.5.nc']
 mtgcm_files = $
 ['/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D1.FISM.TES2.damp7.5.nc',$
  '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D2.FISM.TES2.damp7.5.nc',$
  '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D3.FISM.TES2.damp7.5.nc',$
  '/hao/iris1/bougher/tgcmproc_idlproc/BOUGHER.SWBM12.p1F130L330f.D4.FISM.TES2.damp7.5.nc']

read_mtgcm,state,mtgcm_files

mtgcm_units,state

print,' ' & print,'After read_mtgcm and units conversion:'
print,'mtgcm_ntimes = ',state.mtgcm_ntimes
;print,'mtgcm_mtimes = ',*state.mtgcm_mtimes
;print,'mtgcm_ut = ',*state.mtgcm_ut
print,'mtgcm_lat min,max=',min(*state.mtgcm_lat),max(*state.mtgcm_lat)
print,'mtgcm_lon min,max=',min(*state.mtgcm_lon),max(*state.mtgcm_lon)
print,'mtgcm_lev min,max=',min(*state.mtgcm_lev),max(*state.mtgcm_lev)
print,'mtgcm_nflds = ',state.mtgcm_nflds
print,'mtgcm_fnames = ',*state.mtgcm_fnames

mtgcm_fnames = *state.mtgcm_fnames
mtgcm_hist   = *state.mtgcm_hist
mtgcm_f      = *state.mtgcm_f
for i=0,state.mtgcm_nflds-1 do begin
  print,format="('MTGCM Field ',a,' 4d min,max=',e12.4,',',e12.4)",$
    mtgcm_fnames[i],min(mtgcm_f[*,*,*,*,i]),max(mtgcm_f[*,*,*,*,i])
endfor
;
; Make model slt slices (nlat,nlev,ntime,nflds),at local midnight:
; (where lons at requested slt at available uts are in 3rd dimension)
; Pass state.maven_aveslt as the constant slt at which to make plots.
;
;state.maven_aveslt = 14.0 ; try different slts for debug
state.maven_aveslt = 14.0 
make_sltslices,state,state.maven_aveslt

sltslices = *state.mtgcm_sltslices   ; (nlat,nlev,ntime+1,nflds)
nflds     = state.mtgcm_nflds
print,' ' & print,'Model slt slices:'
for i=0,nflds-1 do begin
  if (mtgcm_hist[i] ne 1) then continue
  print,format="('MTGCM Field ',a,' sltslices min,max=',e12.4,',',e12.4)",$
    mtgcm_fnames[i],min(sltslices[*,*,*,i]),max(sltslices[*,*,*,i])
endfor
;
; Contour constant slt slices (lat,lon): On for testing for each new slt!! 
;plot_sltslices,state,state.maven_aveslt
;
; Dimensions of sltslices(nlats,nlevs,nlons,nflds)
nlats   = n_elements(*state.mtgcm_lat)
nlevs   = n_elements(*state.mtgcm_lev)
lats    = *state.mtgcm_lat      ; mtgcm latitudes (nlats)
sltlons = *state.mtgcm_sltlons  ; sltslices longitudes (nlons)
nlons   = n_elements(sltlons)   ; should be same as ntime
 
ndata = state.maven_ndata
data_lats = *state.maven_lat    ; data latitudes
data_lons = *state.maven_lon    ; data longitudes
data_alts = *state.maven_alt    ; data altitudes
;
; Find index to geopotential field:
ixz = -1
for i=0,nflds-1 do if (mtgcm_fnames[i] eq 'Z') then ixz = i
if (ixz eq -1) then begin
  print,'>>> Could not find field Z'
  return
endif
;
sltslice_z = sltslices[*,*,*,ixz]
print,'sltslice Z min,max (km) = ',min(sltslice_z),max(sltslice_z)
loglin = *state.mtgcm_loglin
mtgcm_interp = fltarr(ndata,nflds) ; final interpolated mtgcm
;
; Loop over maven data:
;
for idata=0,ndata-1 do begin
;
; Bracket data longitude in sltlons:
  i0 = -1 & i1 = -1  
  ier = bracket(sltlons,nlons,data_lons[idata],i0,i1)
  if (ier ne 0) then begin
    print,'>>> WARNING: error return from bracket sltlons:',$
      ' data_lons[idata]=',data_lons[idata]
    print,'nlons=',nlons,' sltlons=' & print,sltlons
    return
  endif
;
; Bracket data latitude in model grid:
  j0 = -1 & j1 = -1
  ier = bracket(lats,nlats,data_lats[idata],j0,j1)
  if (ier ne 0) then begin
    print,'>>> WARNING: error return from bracket lats:',$
      ' data_lats[idata]=',data_lats[idata]
    return
  endif
;
; Bracket data alt using the geopotential column of each
; of the 4 surrounding lat,lon points:
  z_i0j0 = sltslice_z[j0,*,i0]
  z_i0j1 = sltslice_z[j0,*,i1]
  z_i1j1 = sltslice_z[j1,*,i1]
  z_i1j0 = sltslice_z[j1,*,i0]
  k0_i0j0 = -1 & k1_i0j0 = -1
  k0_i0j1 = -1 & k1_i0j1 = -1
  k0_i1j1 = -1 & k1_i1j1 = -1
  k0_i1j0 = -1 & k1_i1j0 = -1
  kerr0 = bracket(z_i0j0,nlevs,data_alts[idata],k0_i0j0,k1_i0j0)  
  kerr1 = bracket(z_i0j1,nlevs,data_alts[idata],k0_i0j1,k1_i0j1)  
  kerr2 = bracket(z_i1j1,nlevs,data_alts[idata],k0_i1j1,k1_i1j1)  
  kerr3 = bracket(z_i1j0,nlevs,data_alts[idata],k0_i1j0,k1_i1j0)  
;
; If alt is out of range, do not interpolate at this data point:
  if (kerr0 ne 0 or kerr1 ne 0 or kerr2 ne 0 or kerr3 ne 0) then begin 
    mtgcm_interp[idata,*] = 1.e36 
    continue
  endif
;
; Fields loop:
;
  for i=0,nflds-1 do begin
    if (mtgcm_fnames[i] eq 'Z') then begin
      mtgcm_interp[idata,i] = 1.e36
      continue ; do not interp geopotential
    endif
;
; Linear or log interpolation to maven altitude:
; (using mean geopotential columns of the 4 points)
;
; Linear interpolation in ht of the 4 surrounding points:
    if (loglin[i] eq 'LIN') then begin ; do linear interp in ht
; i0,j0
      frac = (data_alts[idata]-z_i0j0[k0_i0j0]) / (z_i0j0[k1_i0j0]-z_i0j0[k0_i0j0])
      f_i0j0 = frac*sltslices[j0,k1_i0j0,i0,i]+(1.-frac)*sltslices[j0,k0_i0j0,i0,i]
; i1,j0
      frac = (data_alts[idata]-z_i1j0[k0_i1j0]) / (z_i1j0[k1_i1j0]-z_i1j0[k0_i1j0])
      f_i1j0 = frac*sltslices[j0,k1_i1j0,i1,i]+(1.-frac)*sltslices[j0,k0_i1j0,i1,i]
; i1,j1
      frac = (data_alts[idata]-z_i1j1[k0_i1j1]) / (z_i1j1[k1_i1j1]-z_i1j1[k0_i1j1])
      f_i1j1 = frac*sltslices[j1,k1_i1j1,i1,i]+(1.-frac)*sltslices[j1,k0_i1j1,i1,i]
; i0,j1
      frac = (data_alts[idata]-z_i0j1[k0_i0j1]) / (z_i0j1[k1_i0j1]-z_i0j1[k0_i0j1])
      f_i0j1 = frac*sltslices[j1,k1_i0j1,i0,i]+(1.-frac)*sltslices[j1,k0_i0j1,i0,i]
    endif else begin                      ; do log interp in ht
;
; Log interpolation in ht of the 4 surrounding points:
; i0,j0
      exparg = (alog(sltslices[j0,k1_i0j0,i0,i] / sltslices[j0,k0_i0j0,i0,i]) / $
        (z_i0j0[k1_i0j0]-z_i0j0[k0_i0j0])) * (data_alts[idata]-z_i0j0[k0_i0j0])
      f_i0j0 = sltslices[j0,k0_i0j0,i0,i] * exp(exparg)
; i1,j0
      exparg = (alog(sltslices[j0,k1_i1j0,i1,i] / sltslices[j0,k0_i1j0,i1,i]) / $
        (z_i1j0[k1_i1j0]-z_i1j0[k0_i1j0])) * (data_alts[idata]-z_i1j0[k0_i1j0])
      f_i1j0 = sltslices[j0,k0_i1j0,i1,i] * exp(exparg)
; i1,j1
      exparg = (alog(sltslices[j1,k1_i1j1,i1,i] / sltslices[j1,k0_i1j1,i1,i]) / $
        (z_i1j1[k1_i1j1]-z_i1j1[k0_i1j1])) * (data_alts[idata]-z_i1j1[k0_i1j1])
      f_i1j1 = sltslices[j1,k0_i1j1,i1,i] * exp(exparg)
; i0,j1
      exparg = (alog(sltslices[j1,k1_i0j1,i0,i] / sltslices[j1,k0_i0j1,i0,i]) / $
        (z_i0j1[k1_i0j1]-z_i0j1[k0_i0j1])) * (data_alts[idata]-z_i0j1[k0_i0j1])
      f_i0j1 = sltslices[j1,k0_i0j1,i0,i] * exp(exparg)
    endelse
;
; Finally, do bilinear interpolation in lat,lon:
;
    fracx = (data_lats[idata]-   lats[j0]) / abs(   lats[j1]-   lats[j0])
    fracy = (data_lons[idata]-sltlons[i0]) / abs(sltlons[i1]-sltlons[i0])
    mtgcm_interp[idata,i] = fracx*(fracy*f_i1j1+(1.-fracy)*fracx*f_i0j1)+$
      (1.-fracx)*(fracy*f_i1j0+(1.-fracy)*(1.-fracx)*f_i0j0)

  endfor ; i=0,nflds-1 (model fields)
endfor ; idata=0,ndata-1 (maven data)
*state.mtgcm_interp = mtgcm_interp
;
; Scatter plot interpolated mtgcm:
;
 print,' '
 plot_mtgcm_interp,state
;
; Make table of  interpolated mtgcm:
;
file = 'mtgcm_interp.txt'
openw,lu,file,/get_lun
print,'Opened output text file ',file
date = *state.maven_date
nfiles = n_elements(mtgcm_files)
hist_files = strarr(nfiles)
for i=0,nfiles-1 do begin
  len = strlen(mtgcm_files[i])
  n = strpos(mtgcm_files[i],'BOUGHER') ; assume this is beginning of file name
  hist_files[i] = strmid(mtgcm_files[i],n,len-n)
endfor
printf,lu,format="('MTGCM history files (nfiles=',i3,'):',/,(a))",nfiles,hist_files
printf,lu,format="(/,13x,'Time (UTCG)','    Latitude','   Longitude','    Altitude',12a12)",$
  mtgcm_fnames[0:nflds-1]
for i=0,ndata-1 do begin
  printf,lu,format="(a,3f12.3,12e12.3)",$
    date[i],data_lats[i],data_lons[i],data_alts[i],mtgcm_interp[i,0:nflds-1]
endfor
close,lu
print,'Closed output text file ',file
free_lun,lu

if (state.ps gt 0) then begin
  device,/close_file
  print,'Closed ps file ',psfile
endif
end