;-----------------------------------------------------------------------
;
; 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)
  logdenmin = min(den) & logdenmax = max(den)
  title = 'MAVEN Density: min,max=' + strcompress(string(denmin)) + $
           strcompress(string(denmax)) + ' (kg/km^3)'
  xll = .15 & yll = .15
  xur = .90 & yur = .85
  pos = [xll,yll,xur,yur]
  contour,den,lat,alt,/irregular,nlevels=15,/follow,/nodata,position=pos,$
    xtitle='Latitude (deg)',ytitle='Altitude (km)',title=title,charsize=1.
  plots,lat,alt,psym=3,color=bytscl(den,top=!d.n_colors-1)
;
; 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
;-----------------------------------------------------------------------
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']
 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
     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_f = f

end
;-----------------------------------------------------------------------
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
;
; 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 ; 1st and last lons are probably equal, so skip last
nlats     = n_elements(*state.mtgcm_lat)
nlevs     = n_elements(*state.mtgcm_lev)
fnames    = *state.mtgcm_fnames

lons = fltarr(nlons)
lons = sltlons[0:nlons-1]
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
;-----------------------------------------------------------------------
pro mtgcm_maven,ps=ps
;
; Define state info:
;
state = {                        $
;
; Sat 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_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 
;
; 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
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

plot_sltslices,state,slt
if (keyword_set(ps)) then device,/close_file

end