;  GPM launched 20140227
;  GMI data starts 20140403

pro volcano_gpm_ovp_v2,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='20140512' & shour='00'
;sdate='20140513' & shour='23'
;sdate='20140515' & shour='00'   
;sdate='20140517' & shour='00'
;sdate='20140519' & shour='23'
;sdate='20140522' & shour='23'    
;sdate='20140526' & shour='00'
;sdate='20140528' & shour='00'
;sdate='20140529' & shour='23'  
;2014   05      31    00    40      06      43032      0.15      0.25     16.54   -165.92
;2014   06      02    00    28      53      43061      0.19      0.41     21.21   -163.89
;2014   06      06    00    03      14      43119      0.18      0.43     18.28   -157.00
;2014   06      07    00    45      05      43134      0.18      0.23     11.41   -166.27
;2014   06      09    00    34      26      43163      0.15      0.36     19.75   -165.03
;2014   06      13    00    08      15      43221      0.23      0.45     15.41   -157.86
;2014   06      14    23    57      12      43250      0.18      0.24     18.43   -155.46
;2014   06      16    00    38      27      43265      0.30      0.34     12.02   -164.85
;2014   06      18    00    27      07      43294      0.25      0.37     15.44   -162.53
;2014   06      20    00    14      23      43323      0.22      0.38     14.52   -159.25
;2014   06      25    00    33      11      43396      0.15      0.27     14.68   -163.94
;2014   06      29    00    07      52      43454      0.07      0.11     13.86   -157.60
;2014   06      30    23    57      29      43483      0.13      0.29     19.32   -155.74
;****************************************************************
;  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/'
idir=path_prefix+'mace-group5/gpm/GPM_L3/GPM_3IMERGHH.06/'

;*****************************************************************
;  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
vid=ncdf_varid(fid,'aod_mydatml2')
if vid ne -1 then begin
  ncdf_varget,fid,vid,aod
endif else begin
  aod=make_array(n_elements(lat),value=-9)
endelse
ncdf_close,fid

;  Calculate colongitude
colon=lon
r=where(lon lt 0,c)
if c gt 0 then colon[r]=colon[r]+360.0

;  Calculate the mean lat and lon of the cloudsat track
lon_t=mean(lon)
lat_t=mean(lat)
colon_t=mean(colon)
aod_mean=mean(aod)

;  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=5.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
lcolon1=llon1
if lcolon1 lt 0 then lcolon1=lcolon1+360.0
rcolon1=rlon1
if rcolon1 lt 0 then rcolon1=rcolon1+360.0

;  sample region aresa
ulat_st=max(ulat)
llat_st=min(blat)
llon_st=llon1+box_offset
rlon_st=rlon1-box_offset

;  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,volcano_name,$
  u500_merra,v500_merra
  colon_merra=lon_merra
  r=where(lon_merra lt 0,c)
  if c gt 0 then colon_merra[r]=colon_merra[r]+360.0

;  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,lcolon1,ulat1,rcolon1]
endif

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

;******************************************************************************************
;  List of 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]

;********************
;  Imerg
;********************

;  loop through the days and get a list of files and their start times
imerg_jdays=!NULL & imerg_files=!NULL  ;make sure the array is clobbered before we start

for i=0,n_elements(date_str)-1 do begin
  ;  Create the file name string
  gstr='*.'+date_str[i]+'*.HDF5'
  ;  Get a list of gpm files
  gfiles=file_search(idir+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 imerg_jdays eq !NULL then begin
      imerg_jdays=jday
      imerg_files=gfiles[j]
    endif else begin
      imerg_jdays = [imerg_jdays,jday]
      imerg_files = [imerg_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=min(abs(julian_day_t-imerg_jdays),close_idx)
imerg_file=imerg_files[close_idx]
read_gpm_3bimerg,imerg_file,julian_day_imerg,colon_imerg,lon_imerg,lat_imerg,precip_imerg

;********************************
;  GMI TB
;********************************

;  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 llat_st and lat_s1 le ulat_st and $
    lon_s1 ge llon_st and lon_s1 le rlon_st,count_ovp)
    ;  This is the sub-satellite point line of gpm.  Is it in the box
  jidx=where(lat_s1[*,100] ge llat_st and lat_s1[*,100] le ulat_st and $
    lon_s1[*,100] ge llon_st and lon_s1[*,100] le rlon_st,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)
  mytable=colortable(22,ncolors=253)
  ;254=hot pink               ;gray=255
  ;mytable=[mytable,transpose([238,18,137]),transpose([220,220,220])]
  mytable=[mytable,transpose([255,255,255]),transpose([238,18,137]),transpose([220,220,220])]
  mycbtable=mytable[0:top_color,*]

  ;  Set up the plot
  pxdim=1200 & pydim=900
  xl=0.10 & xr=0.95
  yb=0.10 & yt=0.80
  sx=0.1
  sy=0.05
  numplots_x=2
  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]=0.86;pos[*,3]-0.05  ;horizontal
  cbpos[*,3]=0.87;cbpos[*,1]+0.01
  ;cbpos=[xl,0.85,xr,0.87]
  dx=0.01 & dy=0.01
  ;  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])
 
  ;*********************
  ;  Grid to interpolate to
  ;*********************
  ;  Calculate regular grid array of the 2degree plotting overpass box
  delta=0.1
  ynum=fix((ulat1-llat1)/delta)
  grid_lat=findgen(ynum)*delta+llat1
  xnum=fix((rcolon1-lcolon1)/delta)
  grid_colon=findgen(xnum)*delta+lcolon1
  grid_lon=grid_colon
  result=where(grid_lon gt 180,count)
  if count gt 0 then grid_lon[result]=grid_lon[result]-360.0

  ;***************
  ;  GMI 89H
  ;***************
  
  var=tb_s1[*,*,8]
  var_lat=lat_s1
  var_lon=lon_s1
  var_colon=colon_s1
  dmin_gmi=200 & dmax_gmi=288  ;standard max and min across all plots
  
  ;  Interpolate mask
  r=where(var_lat ne -9999.90 and var_lon ne -9999.90)
  good_var=var[r]
  good_lat=var_lat[r]
  good_lon=var_lon[r]
  good_colon=var_colon[r]
  qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
  grid_mask=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,$
    /kriging,min_points=18,sectors=8,empty_sectors=3,missing=-8888)

  ;  Interpolate data
  r=where(var_lat ne -9999.90 and var_lon ne -9999.90 and var gt -9999,c)
  if c le 10 then r=where(var_lat ne -9999.90 and var_lon ne -9999.90)
  good_var=var[r]
  good_lat=var_lat[r]
  good_lon=var_lon[r]
  good_colon=var_colon[r]
  qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
  grid_var=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,/nearest_neighbor);,$
;    ;/kriging,min_points=26,sectors=8,empty_sectors=3,missing=-9999)
;  ;  minpoints 8 to 5
;  Bytscal the data
data_image=bytscl(grid_var,top=top_color,min=dmin_gmi,max=dmax_gmi)
;  ;result=where(grid_var eq -9999,count)
;  ;if count gt 0 then data_image[result]=255  ;gray
;  ;result=where(grid_var lt 2*dmin_radar and grid_var ne -9999,count)
;  ;result=where(grid_var lt dmin_gmi and grid_var ne -9999,count)
;  ;if count gt 0 then data_image[result]=255  ;gray
result=where(grid_mask eq -8888,count)
if count gt 0 then data_image[result]=253  ;white
;
;  pnum=0  
p0=image(data_image,grid_lon,grid_lat,/current,$
    map_projection='Mercator',limit=limit_vector,grid_units='degrees',$
;;;    center_longitude=(rcolon1-lcolon1)/2.0,$
;;;    /current,dimensions=[pxdim,pydim]
    /box_axes,clip=0,label_position=0,linestyle='dotted',font_size=12,$
     label_format='MapGrid_Labels',$
    rgb_table=mytable,$;map_projection='Mercator',
   position=pos[pnum,*])
;    
;;  Set up the map
;  p0=map('Mercator',$
;    limit=limit_vector,$
;    ;limit=[llat1,lcolon1,ulat1,rcolon1],$
;    font_size=12,$
;    box_axes=1,label_position=0,linestyle='dotted',$
;    position=pos[pnum,*],/buffer)   
;    
;    
;;p0.save,volcano_name+'.'+aerosol_box+'.'+date+'.gpm_ovp.26.png',height=pydim
;;stop
;;  ;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')
  v0=vector(u850_merra,v850_merra,lon_merra,lat_merra,$
  ;v0=vector(u850_merra,v850_merra,colon_merra,lat_merra,$
    vector_style=1,/overplot,$
;    ;/current,position=pos[pnum,*],$
    x_subsample=2,y_subsample=2,head_size=1.5,$
    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='red',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[pnum,*],font_size=12,rgb_table=mycbtable,$
   tickdir=0,border_on=1,title='GMI 89H Temp (K)',textpos=1)
  t1=text(0.05,0.95,'Cloudsat Overpass !C'+date+' '+time_str,$
    font_size=12)
  t1=text(0.2,0.95,'GPM Overpass !C'+ovp_str,font_size=12)
t1=text(0.05,0.92,'850 mb Wind',font_size=14)
t1=text(0.2,0.92,'AOD='+strcompress(aod_mean,/remove_all),font_size=14)
;  
;  
;**************************************************
;  500 mb winds
;**************************************************
if 1 eq 0 then begin
pnum=1
p0=image(data_image,grid_lon,grid_lat,/current,$
  map_projection='Mercator',limit=limit_vector,grid_units='degrees',$
  /box_axes,clip=0,label_position=0,linestyle='dotted',font_size=12,$
  ;label_format='MapGrid_Labels',$
  rgb_table=mytable,$
  position=pos[pnum,*]) 
  v0=vector(u500_merra,v500_merra,lon_merra,lat_merra,$
    vector_style=1,/overplot,$
    x_subsample=2,y_subsample=2,head_size=1.5,$
    color='blue',grid_units='degrees')
  p1=mapcontinents(/countries,/hires)
  p2=plot(lon,lat,/overplot,color='red',thick=2)
  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') 
 t1=text(pos[pnum,0],pos[pnum,3],'500 mb',font_size=14)
 endif
 
 ;********************************
 ;  Plot imerg
 ;********************************
 pnum=1
 r=where(precip_imerg gt -9999.9)
 dmin_precip=min(precip_imerg[r])
 dmax_precip=10;max(precip[r])
 var_image=bytscl(precip_imerg,max=dmax_precip,min=dmin_precip,top=top_color)
 r=where(precip_imerg eq -9999.9,c)
 if c gt 0 then var_image[r]=255
 r=where(precip_imerg eq 0,c)
 if c gt 0 then var_image[r]=253
 p0=image(var_image,lon_imerg,lat_imerg,$
   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_precip,max_value=dmax_precip,$
   map_projection='Mercator',position=pos[pnum,*],/current)
 ;  Plots the hawaii outline
 p1=mapcontinents(/countries,/usa,/hires)
 ;  Puts up a colorbar
 c0=colorbar(range=[dmin_precip,dmax_precip],$
   orientation=0,position=cbpos[pnum,*],font_size=12,rgb_table=mycbtable,$
   tickdir=0,border_on=1,title='IMERG Precip mm/hr',textpos=1)
 t1=text(0.5,0.95,file_basename(imerg_file),font_size=12)
 
  ;  Save the image
  p0.save,volcano_name+'.'+aerosol_box+'.'+date+'.gpm_ovp.n3.png',height=pydim
;endif
stop
end