;********************************************
;  This plots a plan view (horizontal slice on a map)
;  Constant Latitude slice, constant longitude slice, and line segment slice
;  4 plot.
;**************************************************


;function MapGrid_Labels, orientation, location, fractional, defaultlabel
;  ;if (location eq 0) then $
;  ;  return, orientation ? 'Equator' : 'Prime Meridian'
;  degree = '!M' + STRING(176b) ; Use the Math symbol
;  label = STRTRIM(ROUND(ABS(location)),2) + degree
;  label=locaton
;  suffix = orientation ? ((location lt 0) ? 'S' : 'N') : $
;    ((location lt 0) ? 'W' : 'E')
;return, label + suffix


;*****************************
;  Function to plot a circle
;*****************************
;function circle, xcenter,ycenter,radius
;  points = (2 * !PI / 99.0) * FINDGEN(100)
;  x = xcenter + radius * COS(points )
;  y = ycenter + radius * SIN(points )
;  RETURN, TRANSPOSE([[x],[y]])
;end

;*****************************
;  Program to plot cband radar
;*****************************
pro cband_4panel_1ht

;  Choose radar site, date and time
program_name='nosmn'
date_str='20200328' & time_str='*';'085004'  

;program_name='normet'
;date_str='20200327' & time_str='*';'085004'  

;program_name='in2018_v01'
;date_str='20180209' & time_str='13'
;date_str='20180125' & time_str='234004'

;program_name='olympex'
;date_str='20151204' & time_str='163710' 

rho_filter=0  ;no
;rho_filter=1  ;yes

;  Get year month day strings
syy=strmid(date_str,0,4)
smm=strmid(date_str,4,2)
sdd=strmid(date_str,6,2)

shh=strmid(time_str,0,2)
smi=strmid(time_str,2,2)
sss=strmid(time_str,4,2)
;date_str=syy+smm+sdd+shh+smi+sss
;ymd_str=syy+smm+sdd

;  chpc or imac
path_prefix='/Volumes/' 
;path_prefix='/uufs/chpc.utah.edu/common/home/'

;  directory
if program_name eq 'nosmn' then begin  ;comble
  filedir=path_prefix+'mace-group5/field_programs/comble/cband_somna/'
  filestr='nosmn.ss.vol-dbz*noblock.'+syy+smm+sdd+'T'+time_str+'*pyart.41.201.nc'
  radar_name='cband'
endif else if program_name eq 'in2018_v01' then begin
  ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/'
  ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/cfradial/'
  filedir=path_prefix+'mace-group5/field_programs/investigator/in2018_v01/cband/'+date_str+'/'
  filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart.81.401.nc'
  radar_name='cband'
endif else if program_name eq 'in2019_v01' then begin
  filedir=path_prefix+'mace-group5/field_programs/investigator/in2019_v01_gpm_passes/out_vol/'
  filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart*.nc'
  radar_name='cband'
endif else if program_name eq 'olympex' then begin
  filedir=path_prefix+'mace-group6/field_programs/olympex/sband/'+strmid(date_str,4,4)+'/ppi/'
  filestr='olympex_NPOL1_'+date_str+'_'+time_str+'*pyart*.nc'
  radar_name='sband'
  ;  Get the flight path
  nav_date_str='20151213'
  read_nav_olympex,path_prefix,nav_date_str,julian_day_plane,temp,tas,theta,pres,ppmvw,lat_plane,lon_plane,$
    alt_plane,wwind,uwind,vwind,cwc_kng,twc_nev,lwc_nev,iwc_nev,tot_s,cwc_s,rh,cwc_csi,tot_c,tot_p
endif else if program_name eq 'normet' then begin
  filedir=path_prefix+'mace-group5/field_programs/comble/cband/'+date_str+'/'
  filestr=syy+smm+sdd+time_str+'*pyart*.nc'  
  radar_name='cband'
endif

;  file name
filenames=file_search(filedir+filestr,count=numfiles)
print,'numfiles',numfiles

;  Loop through the cband files
for idx=0,numfiles-1 do begin

  filename=filenames[idx]
  print,filename
  if program_name eq 'olympex' then begin
    parts=strsplit(file_basename(filename),'_',/extract)
    shh=strmid(parts[3],0,2)
    smi=strmid(parts[3],2,2)
    sss=strmid(parts[3],4,2)
    time_str=shh+smi+sss
  endif else if program_name eq 'in2019_v01' or program_name eq 'in2018_v01' then begin
    parts=strsplit(file_basename(filename),'-',/extract)
    shh=strmid(parts[3],0,2)
    smi=strmid(parts[3],2,2)
    sss=strmid(parts[3],4,2)
    time_str=shh+smi+sss   
  endif else if program_name eq 'nosmn' then begin
    parts=strsplit(file_basename(filename),'.',/extract)
    shh=strmid(parts[4],9,2)
    smi=strmid(parts[4],11,2)
    sss=strmid(parts[4],13,2)
    time_str=shh+smi+sss
  endif else if program_name eq 'normet' then begin
    parts=strsplit(file_basename(filename),'.',/extract)
    shh=strmid(parts[0],8,2)
    smi=strmid(parts[0],10,2)
    sss=strmid(parts[0],12,2)
    time_str=shh+smi+sss
  endif
  time_tag=shh+':'+smi+':'+sss  ;for image labels
  
  print,'start reading data'
  if program_name eq 'nosmn' then begin
    read_nosmn_cband,filename,lat,lon,point_z,rho,kdp,zdr,dbz,$
      point_x,point_y,base_time,time_offset,lat_radar,lon_radar
  endif else begin  
    fid=ncdf_open(filename)
    vid=ncdf_varid(fid,'base_time') & ncdf_varget,fid,vid,base_time
    vid=ncdf_varid(fid,'time_offset') & ncdf_varget,fid,vid,time_offset  
    ;vid=ncdf_varid(fid,'radar_time') & ncdf_varget,fid,vid,radar_time 
    ;ncdf_attget,fid,vid,'units',astr & astr=string(astr)
    vid=ncdf_varid(fid,'radar_latitude') & ncdf_varget,fid,vid,lat_radar
    vid=ncdf_varid(fid,'radar_longitude') & ncdf_varget,fid,vid,lon_radar 
    vid=ncdf_varid(fid,'point_latitude') & ncdf_varget,fid,vid,lat
    vid=ncdf_varid(fid,'point_longitude') & ncdf_varget,fid,vid,lon      
    vid=ncdf_varid(fid,'point_altitude') & ncdf_varget,fid,vid,alt   
    vid=ncdf_varid(fid,'point_x') & ncdf_varget,fid,vid,point_x
    vid=ncdf_varid(fid,'point_y') & ncdf_varget,fid,vid,point_y
    vid=ncdf_varid(fid,'point_z') & ncdf_varget,fid,vid,point_z
    vid=ncdf_varid(fid,'cross_correlation_ratio') & if vid ne -1 then ncdf_varget,fid,vid,rho
    vid=ncdf_varid(fid,'specific_differential_phase')
    if vid eq -1 then vid=ncdf_varid(fid,'corrected_specific_differential_phase')
    if vid ne -1 then ncdf_varget,fid,vid,kdp
    vid=ncdf_varid(fid,'corrected_differential_reflectivity') & if vid ne -1 then ncdf_varget,fid,vid,zdr
    vid=ncdf_varid(fid,'corrected_reflectivity')  
    if vid eq -1 then vid=ncdf_varid(fid,'DBZH_gmm')
    if vid eq -1 then vid=ncdf_varid(fid,'reflectivity')
    if vid ne -1 then ncdf_varget,fid,vid,dbz
    ncdf_close,fid
  endelse
   
  ;  convert point_x,y,z from m to km
  point_x=point_x/1000.0
  point_y=point_y/1000.0
  point_z=point_z/1000.0
    
  ;  convert time to julian_day
  ;time_offset=0L
  day1=julday(1,1,1970,0,0,0)
  sec_per_day=long(24L*60L*60L)
  julian_day=double(day1+((base_time+time_offset)/sec_per_day))    
  caldat,julian_day,mm,dd,yy,hh,mi,ss
  print,yy,mm,dd,hh,mi,ss

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

  ;  Zoom in on radar, dist_km amount (max distance is 123)
  dist_km=150.0  ;75
  x1=point_x[*,0,0]
  xr=where(x1 ge -1*dist_km and x1 le dist_km,xc)
  y1=point_y[0,*,0]
  yr=where(y1 ge -1*dist_km and y1 le dist_km,yc)    
  
  if dbz ne !NULL then begin
    dbz=dbz[xr,*,*]
    dbz=dbz[*,yr,*]
    
    lat=lat[xr,*,*]
    lon=lon[xr,*,*]
    
    lat=lat[*,yr,*]
    lon=lon[*,yr,*]  
  endif
  
  if kdp ne !NULL then begin  
    kdp=kdp[xr,*,*]
    zdr=zdr[xr,*,*]
    rho=rho[xr,*,*]
 
    kdp=kdp[*,yr,*]
    zdr=zdr[*,yr,*]
    rho=rho[*,yr,*]  
  endif
  
  ;  Height array
  ;height_cband=alt[0,0,*]/1000.0  ;msl
  height_cband=point_z[0,0,*]  ;agl
  
  ;*************************************
  ;  Find flight track at this time
  ;*************************************
  if julian_day_plane ne !NULL then begin
    xtime=5D/(60D*24D)
    r=where(julian_day_plane ge julian_day-xtime and julian_day_plane le julian_day+xtime and $
          lat_plane ge min(lat) and lat_plane le max(lat) and $
          lon_plane ge min(lon) and lon_plane le max(lon) ,c)
    if c gt 1 then begin
      sub_lat_plane=lat_plane[r]
      sub_lon_plane=lon_plane[r]
      sub_alt_plane=alt_plane[r]/1000.0
      print,'found track'
    endif else begin
      print,'no flight track'
    endelse
  endif
  
  ;  Filetag
  parts=strsplit(file_basename(filename),'.',/extract)
  file_tag=parts[0]
  if program_name eq 'nosmn' then file_tag=parts[0]+'.'+parts[4]
  ;file_tag=strmid(file_basename(filename),0,19)   
     
  ;  Set up plot size
  pxdim=1000 & pydim=1000
  ;  Position the plots
  xl=0.10 & xr=0.90
  yb=0.10 & yt=0.92
  sx=0.08
  sy=0.02
  numplots_x=2
  numplots_y=2
  position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
  ;  Colorbar position
  cbpos=pos
  cbpos[*,0]=pos[*,2]+0.06  ;vertical
  cbpos[*,2]=cbpos[*,0]+0.01
  cbpos[*,1]=cbpos[*,1]+0.03
  cbpos[*,3]=cbpos[*,3]-0.03
  
  ;cbpos[*,1]=cbpos[*,3]+0.005 ;horizontal
  ;cbpos[*,3]=cbpos[*,1]+0.01
      
  ;  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(43,ncolors=253)
                   ;254=hot pink               ;gray=255
  mytable=[mytable,transpose([255,255,255]),transpose([238,18,137]),transpose([230,230,230])]
  mycbtable=mytable[0:top_color,*]
    
  mytable2=colortable(72,ncolors=253)
  ;254=hot pink               ;gray=255
  mytable2=[mytable2,transpose([255,255,255]),transpose([238,18,137]),transpose([230,230,230])]
  mycbtable2=mytable2[0:top_color,*]  
  
  ;mytable3=colortable(72,ncolors=10) ;white11             ;gray12
  ;mytable3=[mytable3,transpose([255,255,255]),transpose([230,230,230])]  
  ;                                   light or    blue yel
  mytable3=[[165,0,38],[241,102,64],[254,212,133],[235,247,228],[175,219,234],[108,164,204],[61,90,167],$
    [100,100,100],[255,255,255],[230,230,230]]
  mycbtable3=mytable3[*,0:7]  
  
  ;top_color=254
  ;mytable=colortable(39,ncolors=255)
  ;mytable[0,*]=[250,250,250];  first color =gray
  ;mycbtable=mytable[1:-1,*]    

  ;  Strings
  dx=0.01 & dy=0.01
  fs1=14

  ;Clear plot space
  pnum=0
  p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],$
    axis_style=4,/nodata)
  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])

  ;  radar min and max
  ;dmin_cband=-17 & dmax_cband=31
  ;dmin_cband=-5 & dmax_cband=25
  
  ;r=where(kdp ne -9999,c)
  ;dmin_cband=min(kdp[r])
  ;dmax_cband=max(kdp[r])
  
  ;r=where(reflect ne -9999,c)
  ;dmin_cband=min(reflect[r])
  ;dmax_cband=max(reflect[r])
  
  ;plot_var=reflect
  ;var_name='reflect'
  ;radar_name='SBAND'
  
  
  ;***********************************
  ;  DBZ at 1.5km
  ;  plot dbz where rhohv > 0.9 and dbz > 10.0
  ;***********************************
  ;r=where(dbz ne -9999,c)
  dmin_cband=5;min(dbz[r])
  dmax_cband=50;max(dbz[r])
  ;  1.5 km
  pnum=0
  r=where(height_cband ge 1.5,c)
  sidx=r[0]
  data_var=reform(dbz[*,*,sidx])
  lat_var=reform(lat[*,*,sidx])
  lon_var=reform(lon[*,*,sidx])
  if kdp ne !NULL then kdp_var=reform(kdp[*,*,sidx])
  if rho ne !NULL then rho_var=reform(rho[*,*,sidx])
  ;  Filter dbz for nosmn
  if program_name eq 'nosmn' then begin
    r=where((rho_var le 0.9 and rho_var ne -9999) or (data_var lt 10.0 and data_var ne -9999),c)
    ;r=where(data_var lt 10.0 and data_var ne -9999,c) 
    if c gt 0 then data_var[r]=-8888
    t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-1*dy,'dBZ where rhohv > 0.9 and dBZ > 10',font_size=fs1)
  endif
  data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband)
  r=where(data_var eq -9999,c)
  if c gt 0 then data_image[r]=253  ;white
  r=where(data_var eq -8888,c)
  if c gt 0 then data_image[r]=255  ;gray
  if kdp ne !NULL and program_name ne 'nosmn' then begin
    r=where(kdp_var eq -9999,c)
    if c gt 0 then data_image[r]=255  ;gray
  endif
  llat=min(lat) & ulat=max(lat)
  llon=min(lon) & rlon=max(lon)
  p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$
    /current,dimensions=[pxdim,pydim],grid_units='degrees',$
    rgb_table=mytable,linestyle='dot',$
    label_format='MapGrid_Labels',/box_axes,label_position=0,$
    label_show=1,irregular=0,$
    map_projection='Mercator',position=pos[pnum,*],font_size=fs1)
  p2=mapcontinents(color='black',/hires)
  if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot)
  ;  calculate range rings around ship location (25km)
  lon_cir=fltarr(4,360)
  lat_cir=fltarr(4,360)
  for r=1,4 do begin
    arc_dist=(float(r)*25000.0)/(6371D*1000D)  ;m
    for s=0,359 do begin
      result=ll_arc_distance([lon_radar,lat_radar],arc_dist,s,/degrees)
      lon_cir[r-1,s]=result[0]
      lat_cir[r-1,s]=result[1]
    endfor
  endfor
  for r=0,3 do begin
    p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot,color='gray')
  endfor
  ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,$
  t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-5*dy,$
    string(height_cband[sidx],format='(F4.2)')+' km',font_size=fs1)
  
  ;  Colorbar - horizontal
  ;cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$
  ;  orientation=0,position=[pos[1,0],pos[1,3]+0.5*dy,pos[2,2],pos[2,3]+1.5*dy],font_size=fs1,$
  ;  tickdir=1,border_on=1);,title='dBZ')   
  ;t1=text(pos[1,0]-5*dx,pos[1,3]-1*dy,'dBZ',font_size=fs1)
  ;  Colorbar - vertical
  cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$
    orientation=1,$
    ;position=[pos[pnum,2]+4.5*dx,pos[pnum,1],pos[pnum,2]+5.5*dx,pos[pnum,3]-3*dy],$
    position=cbpos[pnum,*],$
    font_size=fs1,tickdir=1,border_on=1);,title='dBZ')
  t1=text(pos[pnum,2]-5*dx,pos[pnum,3]-5*dy,'dBZ',font_size=fs1)

  ;***********************************
  ;  Zdr
  ;***********************************    
  if zdr ne !NULL then begin
  ;r=where(dbz ne -9999,c)
  dmin_zdr=-1;min(dbz[r])
  dmax_zdr=1;max(dbz[r])
  ;  1.5 km
  pnum=1
  r=where(height_cband ge 1.5,c)
  sidx=r[0]
  data_var=reform(zdr[*,*,sidx])
  lat_var=reform(lat[*,*,sidx])
  lon_var=reform(lon[*,*,sidx])
  kdp_var=reform(kdp[*,*,sidx])
  dbz_var=reform(dbz[*,*,sidx])
  rho_var=reform(rho[*,*,sidx])
  ;  Filter dbz for nosmn
  if program_name eq 'nosmn' then begin
    r=where((rho_var le 0.9 and rho_var ne -9999) or (dbz_var lt 10.0 and dbz_var ne -9999),c)
    if c gt 0 then data_var[r]=-8888
    t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-1*dy,'ZDR where rhohv > 0.9 and dBZ > 10',font_size=fs1)
  endif
  data_image=bytscl(data_var,top=top_color,min=dmin_zdr,max=dmax_zdr)
  r=where(data_var eq -9999,c)
  if c gt 0 then data_image[r]=253  ;white
  r=where(data_var eq -8888,c)
  if c gt 0 then data_image[r]=255  ;gray
  r=where(finite(data_var,/NAN),c)
  if c gt 0 then data_image[r]=253  ;white
  if program_name ne 'nosmn' then begin
    r=where(kdp_var eq -9999,c)
    if c gt 0 then data_image[r]=255  ;gray
  endif
  llat=min(lat) & ulat=max(lat)
  llon=min(lon) & rlon=max(lon)
  p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$
    /current,dimensions=[pxdim,pydim],grid_units='degrees',$
    rgb_table=mytable2,linestyle='dot',$
    label_format='MapGrid_Labels',/box_axes,label_position=0,$
    label_show=1,irregular=0,$
    map_projection='Mercator',position=pos[pnum,*],font_size=fs1)
  p1['Latitudes'].label_show=0   
  p2=mapcontinents(color='black',/hires)
  if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot)
  for r=0,3 do begin
    p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot,color='gray')
  endfor
  t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-5*dy,$  
    string(height_cband[sidx],format='(F4.2)')+' km',font_size=fs1)  
  ;  Colorbar
  ;cb0=colorbar(rgb_table=mycbtable2,range=[dmin_zdr,dmax_zdr],$
  ;  orientation=0,position=[pos[5,0],pos[5,3]+0.5*dy,pos[6,2],pos[6,3]+1.5*dy],font_size=fs1,$
  ;  tickdir=1,border_on=1);,title='dBZ')   
  ;t1=text(pos[5,0]-5*dx,pos[5,3]-1*dy,'ZDR',font_size=fs1)  
  ;  Colorbar - vertical
  cb0=colorbar(rgb_table=mycbtable2,range=[dmin_zdr,dmax_zdr],$
    orientation=1,$
    ;position=[pos[pnum,2]+4.5*dx,pos[pnum,1],pos[pnum,2]+5.5*dx,pos[pnum,3]-3*dy],$
    position=cbpos[pnum,*],$
    font_size=fs1,tickdir=1,border_on=1);,title='dBZ')
  t1=text(pos[pnum,2]-5*dx,pos[pnum,3]-5*dy,'ZDR',font_size=fs1)
  endif
  ;***********************************
  ;  KDP
  ;***********************************
  if kdp ne !NULL then begin
  ;                                   light or    blue yel
  ;mytable3=[[165,0,38],[241,102,64],[254,212,133],[235,247,228],[175,219,234],[108,164,204],[61,90,167],$
  ;  [255,255,255],[230,230,230]]
  kdp_int=[-100,-0.1,-0.05,0,0.05,0.1,1.0,100.0]
  lrho=0.95
  dmin_kdp=-0.2;min(dbz[r])
  dmax_kdp=0.2;max(dbz[r])
  ;  0.5 km
  pnum=2
  r=where(height_cband ge 1.5,c)
  sidx=r[0]
  data_var=reform(kdp[*,*,sidx])
  lat_var=reform(lat[*,*,sidx])
  lon_var=reform(lon[*,*,sidx])
  rho_var=reform(rho[*,*,sidx])
  dbz_var=reform(dbz[*,*,sidx])
  ;  Filter dbz for nosmn
  if program_name eq 'nosmn' then begin
    r=where((rho_var le 0.9 and rho_var ne -9999) or (dbz_var lt 10.0 and dbz_var ne -9999),c)
    if c gt 0 then data_var[r]=-8888
    t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-1*dy,'KDP where rhohv > 0.9 and dBZ > 10',font_size=fs1)
  endif
  ;data_image=bytscl(data_var,top=top_color,min=dmin_kdp,max=dmax_kdp)
  data_image=byte(data_var) & data_image[*]=0B
  for i=0,6 do begin
    r=where(data_var ge kdp_int[i] and data_var lt kdp_int[i+1],c)
    if c gt 0 then data_image[r]=i
    print,kdp_int[i],kdp_int[i+1],c
  endfor
  r=where(data_var eq -9999,c)
  if c gt 0 then data_image[r]=8;253  ;white
  r=where(data_var eq 0,c)
  if c gt 0 then data_image[r]=8;253  ;white
  r=where(data_var eq -8888,c)
  if c gt 0 then data_image[r]=9;255  ;gray
  if kdp ne !NULL and program_name ne 'nosmn' then begin
    r=where(data_var eq -9999,c)
    if c gt 0 then data_image[r]=9;255  ;gray
  endif
  if rho_filter eq 1 then begin
    r=where(rho_var lt lrho and data_var ne 0 and data_var ne -9999,c)
    if c gt 0 then data_image[r]=7;254  ;pink
  endif
  llat=min(lat) & ulat=max(lat)
  llon=min(lon) & rlon=max(lon)
  p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$
    /current,dimensions=[pxdim,pydim],grid_units='degrees',$
    rgb_table=mytable3,linestyle='dot',$
    label_format='MapGrid_Labels',/box_axes,label_position=0,$
    label_show=1,irregular=0,$
    map_projection='Mercator',position=pos[pnum,*],font_size=fs1)
  p1['Longitudes'].label_show=0  
  p2=mapcontinents(color='black',/hires)
  if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot)
  for r=0,3 do begin
    p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot,color='gray')
  endfor
  t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-5*dy,$
    string(height_cband[sidx],format='(F4.2)')+' km',font_size=fs1)
  ;cb0=colorbar(rgb_table=mycbtable,range=[dmin_kdp,dmax_kdp],$
  ;  orientation=0,position=[pos[9,0],pos[9,3]+0.5*dy,pos[10,2],pos[10,3]+1.5*dy],font_size=fs1,$
  ;  tickdir=1,border_on=1);,title='dBZ')
  ;t1=text(pos[9,0]-5*dx,pos[9,3]-1*dy,'KDP',font_size=fs1)
  ;  Colorbar - vertical
  if rho_filter eq 1 then begin
    kdp_int_str=['<-0.1','-0.1','-0.05','0','0.05','0.1','1','>1','rho<0.95']
  endif else begin
    kdp_int_str=['<-0.1','-0.1','-0.05','0','0.05','0.1','1','>1','']
  endelse
  cb0=colorbar(rgb_table=mycbtable3,$;range=[dmin_kdp,dmax_kdp],$  ;mytable2
    orientation=1,tickname=kdp_int_str,textpos=0,major=9,$
    ;position=[pos[pnum,2]+4.5*dx,pos[pnum,1],pos[pnum,2]+5.5*dx,pos[pnum,3]-3*dy],$
    position=cbpos[pnum,*],$
    font_size=fs1,$
    tickdir=1,border_on=1);,title='dBZ')
  t1=text(pos[pnum,2]-5*dx,pos[pnum,3]-5*dy,'KDP',font_size=fs1)
  endif
  
  ;***********************************
  ;  RHOHV
  ;***********************************
  if rho ne !NULL  and 1 eq 0 then begin
  r=where(rho ne -9999,c)
  dmin_rho=0.0;min(rho[r])
  dmax_rho=1.0;max(rho[r])
  ;  0.5 km
  pnum=3
  r=where(height_cband ge 1.5,c)
  sidx=r[0]
  data_var=reform(rho[*,*,sidx])
  lat_var=reform(lat[*,*,sidx])
  lon_var=reform(lon[*,*,sidx])
  llat=min(lat) & ulat=max(lat)
  llon=min(lon) & rlon=max(lon)
  data_image=bytscl(data_var,top=top_color,min=dmin_rho,max=dmax_rho)
  r=where(finite(data_var,/NAN),c)
  if c gt 0 then data_image[r]=253
  r=where(data_var eq -9999,c)
  if c gt 0 then data_image[r]=255
  p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$
    /current,dimensions=[pxdim,pydim],grid_units='degrees',$
    rgb_table=mytable,linestyle='dot',$
    label_format='MapGrid_Labels',/box_axes,label_position=0,$
    label_show=0,irregular=0,$
    map_projection='Mercator',position=pos[pnum,*],font_size=fs1)
  p2=mapcontinents(color='black',/hires)
  if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot)
  for r=0,3 do begin
    p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot,color='gray')
  endfor
  t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-5*dy,string(height_cband[sidx],format='(F4.2)')+' km',$
    font_size=fs1)
  
  ;  Colorbar - vertical
  cb0=colorbar(rgb_table=mycbtable,range=[dmin_rho,dmax_rho],$
    orientation=1,$
    ;position=[pos[pnum,2]+4.5*dx,pos[pnum,1],pos[pnum,2]+5.5*dx,pos[pnum,3]-3*dy],$
    position=cbpos[pnum,*],$
    font_size=fs1,$
    tickdir=1,border_on=1);,title='dBZ')
  t1=text(pos[pnum,2]-8*dx,pos[pnum,3]-5*dy,'RHOHV',font_size=fs1)
  endif
  
  ;***********************************
  ;  DBZ at 1.5km
  ;***********************************
  ;r=where(dbz ne -9999,c)
  dmin_cband=5;min(dbz[r])
  dmax_cband=50;max(dbz[r])
  ;  1.5 km
  pnum=3
  r=where(height_cband ge 1.5,c)
  sidx=r[0]
  data_var=reform(dbz[*,*,sidx])
  lat_var=reform(lat[*,*,sidx])
  lon_var=reform(lon[*,*,sidx])
  if kdp ne !NULL then kdp_var=reform(kdp[*,*,sidx])
  llat=min(lat) & ulat=max(lat)
  llon=min(lon) & rlon=max(lon)
  data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband)
  r=where(data_var eq -9999,c)
  if c gt 0 then data_image[r]=253  ;white
  if kdp ne !NULL and program_name ne 'nosmn' then begin
    r=where(kdp_var eq -9999,c)
    if c gt 0 then data_image[r]=255  ;gray
  endif
  if program_name eq 'nosmn' then begin
    t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-1*dy,'All dBZ values',font_size=fs1)
  endif
  p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$
    /current,dimensions=[pxdim,pydim],grid_units='degrees',$
    rgb_table=mytable,linestyle='dot',$
    label_format='MapGrid_Labels',/box_axes,label_position=0,$
    label_show=0,irregular=0,$
    map_projection='Mercator',position=pos[pnum,*],font_size=fs1)
  p2=mapcontinents(color='black',/hires)
  if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot)
  ;  calculate range rings around ship location (25km)
  lon_cir=fltarr(4,360)
  lat_cir=fltarr(4,360)
  for r=1,4 do begin
    arc_dist=(float(r)*25000.0)/(6371D*1000D)  ;m
    for s=0,359 do begin
      result=ll_arc_distance([lon_radar,lat_radar],arc_dist,s,/degrees)
      lon_cir[r-1,s]=result[0]
      lat_cir[r-1,s]=result[1]
    endfor
  endfor
  for r=0,3 do begin
    p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot,color='gray')
  endfor
  ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,$
  t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-5*dy,$
    string(height_cband[sidx],format='(F4.2)')+' km',font_size=fs1)

  ;  Colorbar - horizontal
  ;cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$
  ;  orientation=0,position=[pos[1,0],pos[1,3]+0.5*dy,pos[2,2],pos[2,3]+1.5*dy],font_size=fs1,$
  ;  tickdir=1,border_on=1);,title='dBZ')
  ;t1=text(pos[1,0]-5*dx,pos[1,3]-1*dy,'dBZ',font_size=fs1)
  ;  Colorbar - vertical
  cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$
    orientation=1,$
    ;position=[pos[pnum,2]+4.5*dx,pos[pnum,1],pos[pnum,2]+5.5*dx,pos[pnum,3]-3*dy],$
    position=cbpos[pnum,*],$
    font_size=fs1,tickdir=1,border_on=1);,title='dBZ')
  t1=text(pos[pnum,2]-5*dx,pos[pnum,3]-5*dy,'dBZ',font_size=fs1)
  
  
  
  
  
  t1=text(pos[3,0],pos[2,3]+2*dy,'25 km range rings',font_size=fs1)
  t1=text(xl,pos[2,3]+2*dy,program_name+' '+date_str+' '+time_tag,font_size=fs1)
  imagename=file_tag+'.4plot.1ht.rho_filter.png'
  p0.save,imagename
endfor




end