;-----------------------------------------------------------------------
;
; 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.
;
  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
;
  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) ; Solar_Lat_Plain Angle (deg)
  lon  = fltarr(mxdata) ; Solar_long_plain Angle (deg)
  den  = fltarr(mxdata) ; MAVEN_Mars_Gram (kg/km^3)
  alt  = fltarr(mxdata) ; MARS_Altitude (km)

  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.
;
; 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:
    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
;   print,' ndata=',ndata,' dd ',dd,' month=',month,' yy=',yy,' hh=',hh,' mm=',mm,' ss=',ss,$
;         ' rdlat=',rdlat,' rdlon=',rdlon,' rdden=',rdden,' rdalt=',rdalt 
    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
    den[ndata-1]=rdden
    alt[ndata-1]=rdalt
;
; Get local time at given longitude and ut:
    dummy = 0.
    slt[ndata-1] = fslt(dummy,ut[ndata-1],lon[ndata-1],"getlt")
  endwhile
;
; 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,' lat=',f8.3,' lon=',f9.3,' den=',e12.4,' alt=',f12.6)",$
;     i,ndata,year[i],day[i],ut[i],slt[i],lat[i],lon[i],den[i],alt[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
  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=3,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
;
; Update state with satellite data:
  state.maven_ndata = ndata
  *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
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

 fields = ['TN' ,'UN' ,'VN' ,'O1' ,'CO' ,'O2' ,'N2' ,'AR' ,'NE' ,'Z'  ]
 loglin = ['LIN','LIN','LIN','LOG','LOG','LOG','LOG','LOG','LOG','LIN']
 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

; Is this the best place to generate proper units for fields needed for prints/plots?
   for ifld=0,nf-1 do begin
     ncdf_varget,ncid,fields[ifld],fld

     if (ntime eq 1) then begin
       ii = ntimes
       f[*,*,*,ii,ifld] = fld
         if (ifld eq 1) then begin  
	     f[*,*,*,ii,ifld] = fld/100.
	 endif    
         if (ifld eq 2) then begin 
	     f[*,*,*,ii,ifld] = fld/100.
	 endif    
     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]
         if (ifld eq 1) then begin  
	     f[*,*,*,ii,ifld] = fld[*,*,*,i]/100.
	 endif    
         if (ifld eq 2) then begin 
	     f[*,*,*,ii,ifld] = fld[*,*,*,i]/100.
	 endif    
       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_f = f

end ; read_mtgcm
;-----------------------------------------------------------------------
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
f = *state.mtgcm_f
ut = *state.mtgcm_ut
lons = *state.mtgcm_lon
;print,'make_sltslices: ntime=',ntime,' nlat=',nlat,' nlon=',nlon,' nfld=',nfld
sltslices = fltarr(nlat,nlev,ntime,nfld)
sltlons = fltarr(ntime) ; longitudes at requested slt
dummy = 0.

for ifld=0,nfld-1 do begin
  for it=0,ntime-1 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
; print,'make_sltslices: sltlons = ' & print,sltlons
; print,'make_sltslices: field ',fnames[ifld],' min,max=',$
;   min(sltslices[*,*,*,ifld]),',',max(sltslices[*,*,*,ifld])
endfor
;
; Sort sltlons into ascending order and add periodic point:
 indices = sort(sltlons)
 sltlons = sltlons[indices]

 nlons = ntime+1
 tmp = fltarr(nlons)
 for i=0,nlons-2 do tmp[i]=sltlons[i]
 tmp[nlons-1] = 180.
 sltlons = tmp
;print,'nlons=',nlons,' sorted sltlons=',sltlons

 sltslices = sltslices[*,*,indices,*]
 tmp = fltarr(nlat,nlev,nlons,nfld)
 tmp[*,*,0:nlons-2,*] = sltslices[*,*,*,*]
 tmp[*,*,nlons-1,*] = tmp[*,*,0,*]
 sltslices = tmp
;
; 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
nlats     = n_elements(*state.mtgcm_lat)
nlevs     = n_elements(*state.mtgcm_lev)
fnames    = *state.mtgcm_fnames

lons = sltlons
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

print,' '
for ifld=0,nflds-1 do begin
  for k=0,nlevs-1,10 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

    if (not state.ps) then window,1,xsize=700,ysize=600,retain=2
    map_set,plat,plon,rot,/cylindrical,position=mappos,/noborder
    image = bytscl(sltslice,top=!d.table_size-3,min=min(sltslice),max=max(sltslice))+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
    fmin = min(sltslice) & fmax = max(sltslice)
    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)]
    title = 'MTGCM ' + fnames[ifld] + ' SLT=' + strcompress(string(format="(f5.2)",slt)) + $
        ' Zp=' + strcompress(string(format="(f6.2)",levs[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

    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 ; pro 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
  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

  for i=0,nflds-2 do begin ; skip Z, assuming it is in last place
    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)]
    fmin = min(f) & fmax = max(f)

    print,'plot_mtgcm_interp: i=',i,' field ',fnames[i],' f min,max=',$
      fmin,fmax,' lat min,max=',min(lat),max(lat),' alt min,max=',$
      min(alt),max(alt)

    title = 'MTGCM '+fnames[i]+': min,max=' + strcompress(string(fmin)) + $
             strcompress(string(fmax))
    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=3,thick=3,color=bytscl(f,top=!d.n_colors-1)
    plots,lat,alt,psym=2,color=bytscl(f,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

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

end ; plot_mtgcm_interp
;-----------------------------------------------------------------------
pro mtgcm_maven,ps=ps
;
; Define state info:
;
state = { $
;
; Maven state data:
;
  maven_ndata:0,                       $ ; number of data points in maven file 
  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)
;
; 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_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)
;
; Other:
;
  ps:0                                 $ ; 1 if plotting to ps, 0 otherwise
}
stateptr = ptr_new(state)
state = *stateptr
setclrtab,-1
state.ps = 0
if (keyword_set(ps)) then begin
  set_plot,'PS'
  device, filename='mtgcm_maven.ps', /color, bits_per_pixel=8
  state.ps = 1
endif else begin
  set_plot,'X'
endelse
;
; Read MAVEN orbit text file, obtaining satellite attitude and density data: 
;
maven_file = 'DD1_SolarLatLong_Dec26toDec31_LT500km.txt'

read_maven,state,maven_file

print,' ' & print,'After read_mtgcm:'
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_zmslt/BOUGHER.SWBM11.p1F130L270d.D1.FISM.TES1.damp7.5.nc',$
  '/hao/iris1/bougher/tgcmproc_zmslt/BOUGHER.SWBM11.p1F130L270d.D2.FISM.TES1.damp7.5.nc',$
  '/hao/iris1/bougher/tgcmproc_zmslt/BOUGHER.SWBM11.p1F130L270d.D3.FISM.TES1.damp7.5.nc',$
  '/hao/iris1/bougher/tgcmproc_zmslt/BOUGHER.SWBM11.p1F130L270d.D4.FISM.TES1.damp7.5.nc']

read_mtgcm,state,mtgcm_files

print,' ' & print,'After read_mtgcm:'
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_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)
;
slt = 0. 
make_sltslices,state,slt

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
  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):
;plot_sltslices,state,slt
;
; Dimensions of sltslices(nlats,nlevs,nlons,nflds)
nlats   = n_elements(*state.mtgcm_lat)
nlevs   = n_elements(*state.mtgcm_lev)
nlons   = state.mtgcm_ntimes+1  ; include periodic point set by make_sltslices
lats    = *state.mtgcm_lat      ; mtgcm latitudes (nlats)
sltlons = *state.mtgcm_sltlons  ; sltslices longitudes (nlons)
 
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
;
; Convert geopotential sltslice from cm to km:
sltslice_z = sltslices[*,*,*,ixz]*1.e-5
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 return
;
; Bracket data latitude in model grid:
  j0 = -1 & j1 = -1
  ier = bracket(lats,nlats,data_lats[idata],j0,j1)
  if (ier ne 0) then return
;
; Bracket data alt using average height of the 4 surrounding points:
  k0 = -1 & k1 = -1
  z = fltarr(4,nlevs)
  z[0,*] = sltslice_z[j0,*,i0] & z[1,*] = sltslice_z[j0,*,i1]
  z[2,*] = sltslice_z[j1,*,i1] & z[3,*] = sltslice_z[j1,*,i0]
  meanz = mean(z,dimension=1) ; meanz(nlevs)
  kerr = bracket(meanz,nlevs,data_alts[idata],k0,k1)
;
; If alt is out of range, do not interpolate at this data point:
  if (kerr ne 0) then begin
    mtgcm_interp[idata,*] = 1.e36 ; "missing mtgcm data"
;   print,'idata=',idata,' alt ',data_alts[idata],' is out of range'
    continue
  endif ;
; Fields loop:
;
  for i=0,nflds-1 do begin
    if (mtgcm_fnames[i] eq 'Z') then continue ; do not interp geopotential
;
; 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
      frac = (data_alts[idata]-meanz[k0]) / (meanz[k1]-meanz[k0])
      f_i0j0 = frac*sltslices[j0,k1,i0,i]+(1.-frac)*sltslices[j0,k0,i0,i]
      f_i1j0 = frac*sltslices[j0,k1,i1,i]+(1.-frac)*sltslices[j0,k0,i1,i]
      f_i1j1 = frac*sltslices[j1,k1,i1,i]+(1.-frac)*sltslices[j1,k0,i1,i]
      f_i0j1 = frac*sltslices[j1,k1,i0,i]+(1.-frac)*sltslices[j1,k0,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,i0,i] / sltslices[j0,k0,i0,i]) / $
        (meanz[k1]-meanz[k0])) * (data_alts[idata]-meanz[k0])
      f_i0j0 = sltslices[j0,k0,i0,i] * exp(exparg)
; i1,j0
      exparg = (alog(sltslices[j0,k1,i1,i] / sltslices[j0,k0,i1,i]) / $
        (meanz[k1]-meanz[k0])) * (data_alts[idata]-meanz[k0])
      f_i1j0 = sltslices[j0,k0,i1,i] * exp(exparg)
; i1,j1
      exparg = (alog(sltslices[j1,k1,i1,i] / sltslices[j1,k0,i1,i]) / $
        (meanz[k1]-meanz[k0])) * (data_alts[idata]-meanz[k0])
      f_i1j1 = sltslices[j1,k0,i1,i] * exp(exparg)
; i0,j1
      exparg = (alog(sltslices[j1,k1,i0,i] / sltslices[j1,k0,i0,i]) / $
        (meanz[k1]-meanz[k0])) * (data_alts[idata]-meanz[k0])
      f_i0j1 = sltslices[j1,k0,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)

; print,format="('idata=',i5,' mtgcm_interp[idata,*]=',9e12.4)",$
;   idata,mtgcm_interp[idata,0:nflds-2] ; assuming Z is in last place!

endfor ; idata=0,ndata-1 (maven data)
*state.mtgcm_interp = mtgcm_interp
;
; Scatter plot interpolated mtgcm:
;
plot_mtgcm_interp,state
;
; Make table of  interpolated mtgcm:
;
file = 'mtgcm_interp.txt'
openw,lu,file,/get_lun
printf,lu,format="(' idata','    Latitude','   Longitude','    Altitude',9a12)",$
  mtgcm_fnames[0:nflds-2]
for i=0,ndata-1 do begin
  printf,lu,format="(i6,3f12.3,9e12.3)",$
    i,data_lats[i],data_lons[i],data_alts[i],mtgcm_interp[i,0:nflds-2]
endfor
close,lu
free_lun,lu
if (keyword_set(ps)) then device,/close_file

end