pro plot_map_gpm_dbz,pnum,dbz,sfc_idx,lat,lon,colon,grid_colon,grid_lat,$
  dmin_radar,dmax_radar,top_color,llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp,$
  pxdim,pydim,mytable,pos,lon_cir,lat_cir,lon_sub,lat_sub,idx_ovp,height_ku_reverse,$
  mycbtable,cbpos,fs1,p0,p1,outfilename
  
;*****************************
;  Plot map - Ku or Ka - dbz
;*****************************

;  pull out surface
var=reform(dbz[sfc_idx,*,*])  ;MS 1

;  Get rid of the -9999.90 lats and lons for mask interpolation
r=where(lat ne -9999.90 and lon ne -9999.90)
good_colon=colon[r]
good_lat=lat[r]
good_var=var[r]
;  Interpolate mask
qhull,good_colon,good_lat,triangles,/delaunay,sphere=s
grid_mask=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=18,sectors=8,empty_sectors=3,missing=-8888)
;min_points=fewer than this, set to missing
;sectors=grid around point divided into pieces 1 to 8
;empty_sectors=max number of sectors that may be empty for point value ;4 to 3
  
;  Get rid of the -9999.90 lats and lons and var for cloud interpolation
r=where(lat ne -9999.90 and lon ne -9999.90 and var gt -9999,c)
if c le 10 then r=where(lat ne -9999.90 and lon ne -9999.90)
good_colon=colon[r]
good_lat=lat[r]
good_var=var[r]
;  Triangulate the cloud data
qhull,good_colon,good_lat,triangles,/delaunay,sphere=s
grid_var=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=8,sectors=8,empty_sectors=3,missing=-9999)
;  minpoints 8 to 5
;  Bytscal the data
data_image=bytscl(grid_var,top=top_color,min=dmin_radar,max=dmax_radar)
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_radar 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

p1=image(data_image,grid_colon,grid_lat,$
  limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$
  center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,$
  /current,dimensions=[pxdim,pydim],grid_units='degrees',$
  /box_axes,$
  ;label_format='MapGrid_Labels',$
  rgb_table=mytable,map_projection='Mercator',position=pos[pnum,*],font_size=fs1)

p2=mapcontinents(color='black',/hires)
;p1.mapgrid.box_axes=1
p1.mapgrid.box_color='black'
p1.mapgrid.color='black'
p1.mapgrid.label_color='black'
p1.mapgrid.linestyle='dot'
p1.mapgrid.label_position=0

for r=0,3 do begin
  p3=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot,color='black')
endfor

p4=plot(lon_sub[idx_ovp],lat_sub[idx_ovp],/overplot,color='black',linestyle=2,thick=2)

c0=colorbar(rgb_table=mycbtable,range=[dmin_radar,dmax_radar],$
  orientation=0,position=reform(cbpos[pnum,*]),font_size=10,$
  tickdir=1,border_on=1)
  
;  Write this out to a netcdf file
idx_ovp_gpm=where(lat ge llat_ovp and lat le ulat_ovp and $
                  lon ge lcolon_ovp and lon le rcolon_ovp,count_ovp_gpm)
data_var=var[idx_ovp_gpm]
lon_var=colon[idx_ovp_gpm]
lat_var=lat[idx_ovp_gpm]

s=size(data_var,/dimensions)
print,outfilename
fid=ncdf_create(outfilename,/clobber)
num_did=ncdf_dimdef(fid,'num',s[0])
vid=ncdf_vardef(fid,'grid_lat',num_did)
vid=ncdf_vardef(fid,'grid_lon',num_did)
vid=ncdf_vardef(fid,'grid_var',num_did)
vid=ncdf_vardef(fid,'surface_height')
ncdf_control,fid,/endef
ncdf_varput,fid,'grid_lat',lat_var
ncdf_varput,fid,'grid_lon',lon_var
ncdf_varput,fid,'grid_var',data_var
ncdf_varput,fid,'surface_height',height_ku_reverse[sfc_idx]

ncdf_close,fid



end