;  GPM launched 20140227
;  GMI data starts 20140403
;volcano_plot_merra_wind_sub
;position_plots
;read_gpm_1cgmi
pro volcano_gpm_ovp,sdate

;  Enter date and hour of the overpass
;  This comes from the overpass list
;  https://home.chpc.utah.edu/~u0079358/realm/volcano/kilauea_hawaii.ovp.2014.txt
;sdate='20140403' & shour='00'
;sdate='20140404' & shour='23'
;sdate='20140408' & shour='00'
sdate='20140410' & shour='23'    
;sdate='20140411' & shour='23'
;sdate='20140417' & shour='00'
;sdate='20140419' & shour='00' 
;sdate='20140427' & shour='23'
;sdate='20140501' & shour='00'
;sdate='20140505' & shour='00' 
;sdate='20140506' & shour='23' 
;sdate='20140513' & shour='23'

;****************************************************************
;  Constants
;****************************************************************

;  Volcano information
volcano_name='kilauea_hawaii
aerosol_box='downstream'

;  Directories
;path_prefix='/Volumes/'  ;imac
path_prefix='/uufs/chpc.utah.edu/common/home/'  ;chpc
if volcano_name eq 'kilauea_hawaii' then begin
  fdir=path_prefix+'mace-group4/cloudsat/subset_'+volcano_name+'/'
endif else begin
  fdir=path_prefix+'mace-group5/cloudsat/subset_'+volcano_name+'/'
endelse

;  GPM data directory
gdir=path_prefix+'mace-group5/gpm/GPM_L1C/GPM_1CGMI.05/'

;*****************************************************************
;  Read the volcano subset data
;*****************************************************************

;  Get the file that matches this overpass date and hour
files1=file_search(fdir+volcano_name+'.'+aerosol_box+'*.ovp.'+$
  sdate+'.'+shour+'*cdf',count=numfiles)
print,files1

;  Get the data from the subset file
fid=ncdf_open(files1[0])
vid=ncdf_varid(fid,'ulat') & ncdf_varget,fid,vid,ulat
vid=ncdf_varid(fid,'ulon') & ncdf_varget,fid,vid,ulon
vid=ncdf_varid(fid,'blat') & ncdf_varget,fid,vid,blat
vid=ncdf_varid(fid,'blon') & ncdf_varget,fid,vid,blon
vid=ncdf_varid(fid,'llat') & ncdf_varget,fid,vid,llat
vid=ncdf_varid(fid,'llon') & ncdf_varget,fid,vid,llon
vid=ncdf_varid(fid,'rlat') & ncdf_varget,fid,vid,rlat
vid=ncdf_varid(fid,'rlon') & ncdf_varget,fid,vid,rlon
vid=ncdf_varid(fid,'julian_day') & ncdf_varget,fid,vid,julian_day  ;fractional days
vid=ncdf_varid(fid,'xlat') & ncdf_varget,fid,vid,xlat  ;lat of volcano
vid=ncdf_varid(fid,'xlon') & ncdf_varget,fid,vid,xlon  ;lon of volcano
vid=ncdf_varid(fid,'lat') & ncdf_varget,fid,vid,lat  ;cloudsat satellite track
vid=ncdf_varid(fid,'lon') & ncdf_varget,fid,vid,lon  ;cloudsat satellite track
ncdf_close,fid

;  Calculate the mean lat and lon of the cloudsat track
lon_t=mean(lon)
lat_t=mean(lat)

;  The time stamp is the first time in the array
julian_day_t=julian_day[0]
caldat,julian_day_t,mm,dd,yy,hh,mi,ss
julian_date,yy,mm,dd,tdoy
;  Get strings of the overpass time
tyear_str=string(yy,format='(I4)')
tmonth_str=string(mm,format='(I02)')
tday_str=string(dd,format='(I02)')
thour_str=string(hh,format='(I02)')
tmin_str=string(mi,format='(I02)')
tsec_str=string(ss,format='(I02)')
tdoy_str=string(tdoy,format='(I03)')
;  yyyymmdd
date=tyear_str+tmonth_str+tday_str
;  hrmiss
time=thour_str+tmin_str+tsec_str
;  hr:mi:ss
time_str=thour_str+':'+tmin_str+':'+tsec_str

;  Calculate big box around the region
box_offset=3.0
ulat1=max(ulat)+box_offset & llat1=min(blat)-box_offset
;  This fixes the map box if the wind is blowing from the west.
slon=min(llon)
tlon=min(rlon)
if slon lt tlon then begin
  llon1=min(llon)-box_offset & rlon1=max(rlon)+box_offset
endif else if tlon lt slon then begin
  llon1=min(rlon)-box_offset & rlon1=max(llon)+box_offset
endif
co_llon1=llon1
if co_llon1 lt 0 then co_llon1=co_llon1+360.0
co_rlon1=rlon1
if co_rlon1 lt 0 then co_rlon1=co_rlon1+360.0

;  Get the merra winds
volcano_plot_merra_wind_sub,date,time,xlat,xlon,llat1,llon1,ulat1,rlon1,$
  u850_merra,v850_merra,lon_merra,lat_merra,uwind,vwind,lat_f,lon_f

;  Map limit vector
if rlon1 le 180 and llon1 lt 0 then begin
  limit_vector=[llat1,llon1,ulat1,rlon1]
endif else if rlon1 gt 180 or (llon1 gt 0 and rlon1 lt 0) then begin
  limit_vector=[llat1,co_llon1,ulat1,co_rlon1]
endif

;******************************************
;  Now look for gpm overpasses
;******************************************

;******************************************************************************************
;  List of gpm files 4 consecutive days centered on overpass
;******************************************************************************************
;  2 days before
julian_day_tb2=julian_day_t-2D
caldat,julian_day_tb2,tb2month,tb2day,tb2year,tb2hour,tb2min,tb2sec
julian_date,tb2year,tb2month,tb2day,tb2doy
tb2year_str=string(tb2year,format='(I4)')
tb2month_str=string(tb2month,format='(I02)')
tb2day_str=string(tb2day,format='(I02)')
tb2doy_str=string(tb2doy,format='(I03)')
tb2date_str=tb2year_str+tb2month_str+tb2day_str
;  Day before
julian_day_tb=julian_day_t-1D
caldat,julian_day_tb,tbmonth,tbday,tbyear,tbhour,tbmin,tbsec
julian_date,tbyear,tbmonth,tbday,tbdoy
tbyear_str=string(tbyear,format='(I4)')
tbmonth_str=string(tbmonth,format='(I02)')
tbday_str=string(tbday,format='(I02)')
tbdoy_str=string(tbdoy,format='(I03)')
tbdate_str=tbyear_str+tbmonth_str+tbday_str
;  Day after
julian_day_ta=julian_day_t+1D
caldat,julian_day_ta,tamonth,taday,tayear,tahour,tamin,tasec
julian_date,tayear,tamonth,taday,tadoy
tayear_str=string(tayear,format='(I4)')
tamonth_str=string(tamonth,format='(I02)')
taday_str=string(taday,format='(I02)')
tadoy_str=string(tadoy,format='(I03)')
tadate_str=tayear_str+tamonth_str+taday_str
;  2 days after
julian_day_ta2=julian_day_t+2D
caldat,julian_day_ta2,ta2month,ta2day,ta2year,ta2hour,ta2min,ta2sec
julian_date,ta2year,ta2month,ta2day,ta2doy
ta2year_str=string(ta2year,format='(I4)')
ta2month_str=string(ta2month,format='(I02)')
ta2day_str=string(ta2day,format='(I02)')
ta2doy_str=string(ta2doy,format='(I03)')
ta2date_str=ta2year_str+ta2month_str+ta2day_str
;  Array of the 4 consecutive days
year_str=[tb2year_str,tbyear_str,tyear_str,tayear_str,ta2year_str]
doy_str=[tb2doy_str,tbdoy_str,tdoy_str,tadoy_str,ta2doy_str]
date_str=[tb2date_str,tbdate_str,date,tadate_str,ta2date_str]

;  loop through the days and get a list of gpm files and their start times
gpm_jdays=!NULL & gpm_files=!NULL  ;make sure the array is clobbered before we start
prev_gpm_file=''

for i=0,n_elements(date_str)-1 do begin
  ;  Create the gpm file name string
  gstr='*.'+date_str[i]+'*.HDF5'
  ;  Get a list of gpm files
  gfiles=file_search(gdir+gstr,count=num_g)
  ;  This makes an array of files and an array of file starting times in julian days
  for j=0,num_g-1 do begin
    ;  Extract the gpm file date and time from filename
    parts = strsplit(file_basename(gfiles[j]),'.',/extract)
    year=fix(strmid(parts[4],0,4))
    month=fix(strmid(parts[4],4,2))
    day=fix(strmid(parts[4],6,2))
    hh=fix(strmid(parts[4],10,2))
    mm=fix(strmid(parts[4],12,2))
    ss=fix(strmid(parts[4],14,2))
    jday=julday(month,day,year,hh,mm,ss)
    ;print,file_basename(gfiles[j])
    ;print,year,month,day,hh,mm,ss
    ;  Accumulate these start times and filenames in an array
    if gpm_jdays eq !NULL then begin
      gpm_jdays=jday
      gpm_files=gfiles[j]
    endif else begin
      gpm_jdays = [gpm_jdays,jday]
      gpm_files = [gpm_files,gfiles[j]]
    endelse
  endfor   ;end of loop through files
endfor  ;end of loop through 3 consecutive days

;  sort the files by closest in time
abs_jday=abs(julian_day_t-gpm_jdays)
sidx=sort(abs_jday)
sort_gpm_files=gpm_files[sidx]
sort_gpm_times=gpm_jdays[sidx]

;  Loop through the files and find the one closest in time and in gray box
m=0
count_ovp=0

while count_ovp eq 0 and m lt n_elements(sort_gpm_files) do begin

  close_gpm_file=sort_gpm_files[m]
  close_gpm_time=sort_gpm_times[m]
  ;  Read the gpm data
  read_gpm_1cgmi,close_gpm_file,$
    julian_day_s1,tb_s1,colon_s1,lon_s1,lat_s1,$
    julian_day_s2,tb_s2,colon_s2,lon_s2,lat_s2
  ;  See of there is any gpm data in the gray box
  cidx=where(lat_s1 ge llat1 and lat_s1 le ulat1 and $
    lon_s1 ge llon1 and lon_s1 le rlon1,count_ovp)
    ;  This is the sub-satellite point line of gpm.  Is it in the box
  jidx=where(lat_s1[*,100] ge llat1 and lat_s1[*,100] le ulat1 and $
    lon_s1[*,100] ge llon1 and lon_s1[*,100] le rlon1,cjidx)
  
  print,m,' ',count_ovp,' ',close_gpm_file
  if count_ovp eq 0 then m=m+1
endwhile

if close_gpm_file ne prev_gpm_file and count_ovp gt 0 then begin
  prev_gpm_file=close_gpm_file
  print,'plot this one'
  
  ;  GPM overpass time
  parts=strsplit(file_basename(close_gpm_file),'.',/extract)
  orbit=parts[5]
  ovp_time=julian_day_s1[jidx[0]]
  caldat,ovp_time,omonth,oday,oyear,ohour,omin,osec
  oyear_str=string(oyear,format='(I4)')
  omonth_str=string(omonth,format='(I02)')
  oday_str=string(oday,format='(I02)')
  ohour_str=string(ohour,format='(I02)')
  omin_str=string(omin,format='(I02)')
  osec_str=string(osec,format='(I02)')
  ovp_str=oyear_str+omonth_str+oday_str+' '+ohour_str+':'+omin_str+' '+string(orbit)
  
  ;  Cloudsat volcano overpass date and time in a string
  sym_str=tyear_str+tmonth_str+tday_str+' '+thour_str+':'+tmin_str
  
  ;  Colortable
  ct=['red','orange','green','blue','saddle brown','dodger blue','spring green',$
      'gray','crimson','olive','purple','magenta']
    
  ;  Top is the last color to scale 256 colors, 0-255
  top_color=252
  ;  Colortable  0-252  253=white
  mytable=colortable(39,ncolors=254)
  ;254=hot pink               ;gray=255
  mytable=[mytable,transpose([238,18,137]),transpose([220,220,220])]
  mycbtable=mytable[0:top_color,*]

  ;  Set up the plot
  pxdim=800 & pydim=800
  xl=0.1 & xr=0.90
  yb=0.2 & yt=0.80
  sx=0.05
  sy=0.09
  numplots_x=1
  numplots_y=1
  position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
  ;  Colorbar position
  cbpos=pos
  ;cbpos[*,0]=pos[*,2]+0.08  ;vertical
  ;cbpos[*,2]=cbpos[*,0]+0.01
  cbpos[*,1]=pos[*,3]-0.05  ;horizontal
  cbpos[*,3]=cbpos[*,1]+0.01
  cbpos=[xl,0.85,xr,0.87]
  ;  Clear plot space
  pnum=0
  p0=plot([0,1],[0,1],position=pos[pnum,*],axis_style=4,/nodata,/buffer,$
    dimensions=[pxdim,pydim])
  d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
  isx=(d[0,1,0]-d[0,0,0])
  isy=(d[1,1,0]-d[1,0,0])
  
  chan89_gmi=tb_s1[*,*,8]
  lon_gmi=lon_s1
  lat_gmi=lat_s1
  ;dmin_gmi=min(chan89_gmi[cidx]) & dmax_gmi=max(chan89_gmi[cidx])
  dmin_gmi=200 & dmax_gmi=288  ;standard max and min across all plots
  p0=image(chan89_gmi[cidx],lon_gmi[cidx],lat_gmi[cidx],$
     limit=limit_vector,$
     grid_units='degrees',clip=0,$
     box_axes=1,label_position=0,linestyle='dotted',font_size=12,$
     label_format='MapGrid_Labels',$
     rgb_table=mytable,min_value=dmin_gmi,max_value=dmax_gmi,$
     map_projection='Mercator',position=pos[pnum,*],/current)
  ;  This plots the wind vectors   
  v0=vector(u850_merra,v850_merra,lon_merra,lat_merra,$
    vector_style=1,/overplot,$
    x_subsample=2,y_subsample=2,$
    color='blue',grid_units='degrees')
  ;a1=arrow([xlon,lon_f],[xlat,lat_f],color='green',/data,/current)
  ;  Plots the hawaii outline
  p1=mapcontinents(/countries,/hires)
  ; Plot the cloudsat track lat lon points in the region
  p2=plot(lon,lat,/overplot,color='pink',thick=2)
  ;  circle the volcano
  s2=plot([xlon,xlon],[xlat,xlat],/overplot,symbol='o',color='deep sky blue',$
    sym_thick=1,sym_size=2,linestyle=6)
  ;  Draw the gray box
  p1=plot(blon,blat,/overplot,color='gray',thick=3)
  p1=plot(ulon,ulat,/overplot,color='gray',thick=3)
  p1=plot(llon,llat,/overplot,color='gray',thick=3)
  p1=plot(rlon,rlat,/overplot,color='gray',thick=3)
  ;  Star is center of cloudsat track
  p1=symbol(lon_t,lat_t,'star',/data,sym_color='red',$;label_string=sym_str,$
    label_font_size=10,label_color='red')
  ;  Puts up a colorbar
  c0=colorbar(range=[dmin_gmi,dmax_gmi],$  
     orientation=0,position=cbpos,font_size=12,rgb_table=mycbtable,$
     tickdir=0,border_on=1,title='GMI 89H Temp (K)',textpos=1)
  t1=text(0.1,0.95,'Cloudsat Overpass !C'+date+' '+time_str,$
    font_size=12)
  t1=text(0.4,0.95,'GPM Overpass !C'+ovp_str,font_size=12)
  ;  Save the image
  p0.save,volcano_name+'.'+aerosol_box+'.'+date+'.gpm_ovp.png',height=pydim
endif
stop
end