;****************************************************************
;  Text file format located here:
;  https://www.ready.noaa.gov/hysplitusersguide/S263.htm
;****************************************************************

;****************************************************************
;  INPUT filename  name of ascii trajectory file
;  
;  OUTPUT julian_day_traj
;         age_traj
;         lat_traj
;         lon_traj
;         ht_traj
;         pres_traj
;****************************************************************

pro read_hysplit_traj,filename,julian_day_traj,age_traj,$
  lat_traj,lon_traj,ht_traj,pres_traj,lat_start,lon_start

;  File name
;path_prefix='/Volumes/'
;path_prefix='/uufs/chpc.utah.edu/common/home/'
;fdir=path_prefix+'mace-group4/modis/hysplit/marcus_trajectories/2018_1_5_7/'
;filename=fdir+'traj_-65.62_108.75jan1500summer2018010507'
;filename='hysplit_1-6-20_1000.txt'
;filename='hysplit_1-6-20_2000.txt'
;filename='0000_11_24.txt'

;  Open the file
openr,1,filename
;  Record #1 I8-Number of meteorological grids used in calculation
readf,1,num_grids,format='(I8)'
;  Records Loop #2 through the number of grids
;  A8 - Meteorological Model identification
;  5I6 - Data file starting Year, Month, Day, Hour, Forecast Hour
met_model=''
yy_grid=make_array(/int,num_grids)
mm_grid=make_array(/int,num_grids)
dd_grid=make_array(/int,num_grids)
hh_grid=make_array(/int,num_grids)
fh_grid=make_array(/int,num_grids)
for i=0,num_grids-1 do begin
  readf,1,met_model,yy,mm,dd,hh,fh,format='(A8,5(4X,I2))'
  yy_grid[i]=yy
  mm_grid[i]=mm
  dd_grid[i]=dd
  hh_grid[i]=hh
  fh_grid[i]=fh
endfor
;  Record #3 I6-number of different trajectories in file
;  1X,A8 - direction of trajectory calculation (forward,backward)
;  1X,A8 - vertical motion calculation method (omega, theta, ...)
dir_traj=''
vmotion=''
readf,1,num_traj,dir_traj,vmotion,format='(I6,1X,A8,1X,A8)'
;  Record Loop #4 through the number of different trajectories in file
;  4I6 - starting year, month, day, hour
;  2F9.3 - starting latitude, longitude
;  F8.1 - starting level above ground (meters)
yy_start=make_array(/int,num_traj)
mm_start=make_array(/int,num_traj)
dd_start=make_array(/int,num_traj)
hh_start=make_array(/int,num_traj)
lat_start=make_array(/float,num_traj)
lon_start=make_array(/float,num_traj)
ht_start=make_array(/float,num_traj)
for i=0,num_traj-1 do begin
  readf,1,yy,mm,dd,hh,lat,lon,ht,format='(4(I6),2(F9.3),F8.1)'
  yy_start[i]=yy
  mm_start[i]=mm
  dd_start[i]=dd
  hh_start[i]=hh
  lat_start[i]=lat
  lon_start[i]=lon  
  ht_start[i]=ht   
endfor
;  Record #5
;  I6 - number (n) of diagnostic output variables
;  n(1X,A8) - label identification of each variable (PRESSURE, THETA, ...)
readf,1,num_var,format='(I6)'
;Record Loop #6 through the number of hours in the simulation
;I6 - trajectory number
;I6 - meteorological grid number or antecedent trajectory number
;5I6 - year month day hour minute of the point
;I6 - forecast hour at point
;F8.1 - age of the trajectory in hours
;2F9.3 - position latitude and longitude
;1X,F8.1 - position height in meters above ground
;n(1X,F8.1) - n diagnostic output variables (1st to be output is always pressure)
num=10000
traj_id=make_array(/int,num)
grid_id=make_array(/int,num)
yy_traj=make_array(/int,num)
mm_traj=make_array(/int,num)
dd_traj=make_array(/int,num)
hh_traj=make_array(/int,num)
mi_traj=make_array(/float,num)
fh_traj=make_array(/float,num)
age_traj=make_array(/float,num)
lat_traj=make_array(/float,num)
lon_traj=make_array(/float,num)
ht_traj=make_array(/float,num)
pres_traj=make_array(/float,num)
i=0
while (~ EOF(1)) do begin
  readf,1,tid,gid,yy,mm,dd,hh,mi,fh,age,lat,lon,ht,pres,format='(8(I6),F8.1,2(F9.3),1X,F8.1,1X,F8.1)'
  traj_id[i]=tid
  grid_id[i]=gid
  yy_traj[i]=yy
  mm_traj[i]=mm
  dd_traj[i]=dd
  hh_traj[i]=hh
  mi_traj[i]=mi
  fh_traj[i]=fh
  age_traj[i]=age
  lat_traj[i]=lat
  lon_traj[i]=lon
  ht_traj[i]=ht
  pres_traj[i]=pres
  i=i+1
endwhile
traj_id=traj_id[0:i-1]
grid_id=grid_id[0:i-1]
yy_traj=yy_traj[0:i-1]
mm_traj=mm_traj[0:i-1]
dd_traj=dd_traj[0:i-1]
hh_traj=hh_traj[0:i-1]
mi_traj=mi_traj[0:i-1]
fh_traj=fh_traj[0:i-1]
age_traj=age_traj[0:i-1]
lat_traj=lat_traj[0:i-1]
lon_traj=lon_traj[0:i-1]
ht_traj=ht_traj[0:i-1]
pres_traj=pres_traj[0:i-1]
close,1
;print,traj_id[-1],grid_id[-1],yy_traj[-1],mm_traj[-1],dd_traj[-1],$
;hh_traj[-1],mi_traj[-1],fh_traj[-1],age_traj[-1],lat_traj[-1],$
;lon_traj[-1],ht_traj[-1],pres_traj[-1]

;  Convert these date-time strings to IDL julian day
yy_traj='20'+string(yy_traj,format='(I02)')
julian_day_traj=julday(mm_traj,dd_traj,yy_traj,hh_traj,mi_traj,0)

;  Turn the plot on or off
if 1 eq 0 then begin


;  Set up the map
pxdim=900  ;  pixel size of the image in the x direction
pydim=850   ;  pixel size of the image in the y direction
;  These values divide up the plot space using the subroutine
;  position plots.  The location of the plots will be defined by
;  pos
xl=0.10 & xr=0.95   
yb=0.60 & yt=0.97
sx=0.05
sy=0.07
numplots_x=1
numplots_y=2
position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
;;  Lat lon range of the map  ;antarctic
;ulat=-40.0 ;upper lat
;llat=-70.0 ;lower lat
;llon=70 ;left lon
;rlon=150  ;right lon

ulat=40;max(lat_traj)+5. ;upper lat
llat=10;min(lat_traj)-5. ;lower lat
llon=-179.0;max(lon_traj)+5.   ;left lon
rlon=-140.0;min(lon_traj)-5.  ;right lon



map_fontsize=10  ;map font size
continent_color='green'  ;color to draw the continents
fs1=7
;  This the the date-time string format for the x axis
dummy=label_date(date_format=['%N/%D/%Z!C%H:%I'])

;  Plot 0 will be the map
map_pos=[0.1,0.1,0.9,0.5]
m0=map('mercator',$
  limit=[llat,llon,ulat,rlon],/buffer,$
  ;center_longitude=[((rlon-llon)/2.0)+llon],$
  box_axes=1,label_position=0,linestyle='dotted',$
  position=map_pos,font_size=map_fontsize)
mc=mapcontinents(/countries,color=continent_color,thick=1.5)
mc=mapcontinents(/USA,color=continent_color,thick=1.5)
;m0.save,'test_map2.png',height=700

;  Plot the track projection
p1=plot(lon_traj,lat_traj,/overplot,color='blue')
;  Plot the track initial trajectory
p0=plot(lon_start,lat_start,/overplot,color='black',symbol='diamond',thick=1.5,sym_size=1)

;  Plot 1 will be a time series of altitude
pnum=0
p0=plot(julian_day_traj,ht_traj,position=pos[pnum,*],/current,font_size=fs1,$
        xtickformat='label_date',ytitle='Altitude',xtitle='Time (UTC)',color='blue',$
        xrange=[min(julian_day_traj),max(julian_day_traj)]);,$
        ;symbol='o',sym_size=0.5)
      
;  Plot 2 will be a time series of age        
pnum=1
p0=plot(julian_day_traj,age_traj,position=pos[pnum,*],/current,font_size=fs1,$
        xtickformat='label_date',ytitle='Age',color='blue',xrange=[min(julian_day_traj),max(julian_day_traj)])        

;  This is the name of the image file 
m0.save,'traj_plot.png',height=pydim  

endif  ;turn off plot

end