Index: geo_match_3d_rr_or_z_comparisons.pro
===================================================================
--- geo_match_3d_rr_or_z_comparisons.pro	(revision 460)
+++ geo_match_3d_rr_or_z_comparisons.pro	(revision 502)
@@ -9,12 +9,12 @@
 ;
 ; DESCRIPTION
 ; -----------
-; Performs a case-by-case statistical analysis of geometry-matched PR and GR
+; Performs a case-by-case statistical analysis of geometry-matched (D)PR and GR
 ; rain rate from data contained in a geo-match netCDF file.  Rain rate for GR
 ; is taken from the geo-match file if this field is flagged as available,
 ; otherwise it is derived from the volume-averaged GR reflectivity using a Z-R 
-; relationship.  PR rainrate is the volume-averaged rain rate stored in the
-; netCDF file and previously derived from the 3-D rainrate in the 2A-25 product.
+; relationship.  (D)PR rainrate is the volume-averaged rain rate stored in the
+; netCDF file and previously derived from the 3-D rainrate in the 2A product.
 ;
 ; INTERNAL MODULES
 ; ----------------
@@ -23,22 +23,17 @@
 ;
 ; 2) geo_match_xxx_plots - Workhorse procedure to read data, compute statistics,
 ;                          create vertical profiles, histogram, scatter plots, 
-;                          and tablulations of PR-GR rainrate or reflectivity
-;                          differences, and display PR and GR reflectivity and
+;                          and tablulations of (D)PR-GR rainrate or reflectivity
+;                          differences, and display (D)PR and GR reflectivity and
 ;                          rainrate and GR dual-pol field PPI plots in an
 ;                          animation sequence.
 ;
-; 3) print_table_headers - Does what it says.
 ;
-;
 ; NON-SYSTEM ROUTINES CALLED
 ; --------------------------
 ; 1) fprep_geo_match_profiles() or fprep_dpr_geo_match_profiles()
-; 2) z_r_rainrate()
-; 3) calc_geo_pr_gv_meandiffs_wght_idx
-; 4) plot_scatter_by_bb_prox_ps
-; 5) plot_scatter_by_bb_prox
-; 6) plot_geo_match_ppi_anim_ps
+; 2) select_geomatch_subarea()
+; 3) render_rr_or_z_plots()
 ;
 ;
 ; HISTORY
@@ -116,6 +111,22 @@
 ;   NCPATH and/or the default netCDF file path on startup of the procedure, if
 ;   the RECALL_NCPATH keyword is set and user system variable is defined in the
 ;   IDL session.
+; 04/09/15  Morris/GPM GV/SAIC
+; - MAJOR REWRITE FOR INTERACTIVE SELECTION OF ANALYSIS SUB-AREAS.
+; - Extracted internal function plot_sweep_2_zbuf_4xsec from this file into a
+;   stand alone function module common_utils/plot_sweep_2_zbuf_4xsec.pro.
+; - Extracted major blocks of code and the internal function print_table_headers
+;   that do the statistics, plots, and PDF file output into a new stand alone
+;   function module display/render_rr_or_z_plots.pro.
+; - Calls new function select_geomatch_subarea() to clip the data arrays to be
+;   analyzed to the area around a user selected point.
+; 04/28/15  Morris/GPM GV/SAIC
+; - Fixed oversights in reading and processing RC and RP fields from DPR files.
+; - Added structure element rr_field_used to indicate which GR rain rate field
+;   is being processed.
+; - Added tests for existence of SAVE_DIR directory before passing it along from
+;   main routine to geo_match_xxx_plots.
+; - Handle empty string for SUBSET_METHOD keyword parameter.
 ;
 ;
 ; EMAIL QUESTIONS OR COMMENTS TO:
@@ -124,78 +135,6 @@
 ;-
 ;===============================================================================
 ;
-; MODULE 3:  print_table_headers
-;
-
-pro print_table_headers, var1, var2, intype, BB_LAYERS=bb_layers, $
-                         PS_UNIT=tempunit
-
-IF N_ELEMENTS(tempunit) EQ 1 THEN do_ps=1 ELSE do_ps=0
-
-; pad field type to 2 chars as needed
-CASE intype OF
-    'Z' : type=intype+' '
-   'RR' : type=intype
-   'ZR' : type=intype
-   ELSE : message, 'illegal value for intype, must be Z, ZR, or RR'
-ENDCASE
-
-; set up spacing based on lengths of var1 and var2
-CASE (STRLEN(var1)*10+STRLEN(var2)) OF
-   22 : BEGIN
-           diffvar = ' '+var1+'-'+var2
-           maxvars = ' '+ var1 + 'Max'+type+'   '+ var2 +'Max'+type
-        END
-   23 : BEGIN
-           diffvar = var1+'-'+var2
-           maxvars = ' '+ var1 + 'Max'+type+'  '+ var2 +'Max'+type
-        END
-   32 : BEGIN
-           diffvar = var1+'-'+var2
-           maxvars = var1 + 'Max'+type+'   '+ var2 +'Max'+type
-        END
-   ELSE : message, 'illegal string lengths for var1 and var2, must sum to 4 or 5'
-ENDCASE
-
-IF N_ELEMENTS(bb_layers) EQ 0 THEN BEGIN
-   ; print the header for stats broken out by CAPPI levels and rain type
-   print, '' & IF (do_ps EQ 1) THEN printf, tempunit, ''
-   textout = 'Statistics grouped by fixed height levels (km):'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-   print, '' & IF (do_ps EQ 1) THEN printf, tempunit, ''
-   textout = ' Vert. |   Any Rain Type  |    Stratiform    |' $
-             +'    Convective     |     Dataset Statistics      |'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-   textout = ' Layer |  '+diffvar+'   NumPts |  '+diffvar+'   NumPts |' $
-             +'  '+diffvar+'   NumPts  | AvgDist  '+maxvars+' |'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-   textout = ' ----- | -------   ------ | -------   ------ |' $
-             +' -------   ------  | -------  --------  -------- |'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-
-ENDIF ELSE BEGIN
-   ; print the header for stats broken out by BB proximity
-   print, '' & IF (do_ps EQ 1) THEN printf, tempunit, ''
-   print, '' & IF (do_ps EQ 1) THEN printf, tempunit, ''
-   textout = 'Statistics grouped by proximity to Bright Band:'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-   print, '' & IF (do_ps EQ 1) THEN printf, tempunit, ''
-   textout = 'Surface|   Any Rain Type  |    Stratiform    |' $
-             +'    Convective     |     Dataset Statistics      |'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-   textout = ' type  |  '+diffvar+'   NumPts |  '+diffvar+'   NumPts |' $
-             +'  '+diffvar+'   NumPts  | AvgDist  '+maxvars+' |'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-   textout = ' ----- | -------   ------ | -------   ------ |' $
-             +' -------   ------  | -------  --------  -------- |'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-
-ENDELSE
-
-end
-
-;===============================================================================
-;
 ; MODULE 2:  geo_match_xxx_plots
 ;
 ; DESCRIPTION
@@ -241,7 +180,9 @@
                               hide_rntype, pr_or_dpr, PS_DIR=ps_dir, B_W=b_w, $
                               S2KU=s2ku, ZR=zr_force, ALT_BB_HGT=alt_bb_hgt, $
                               GR_RR_FIELD=gr_rr_field, BATCH=batch, $
-                              MAX_RANGE=max_range
+                              MAX_RANGE=max_range, SUBSET_METHOD=submeth, $
+                              MIN_FOR_SUBSET=subthresh, SAVE_DIR=save_dir, $
+                              STEP_MANUAL=step_manual
 
 ; "include" file for read_geo_match_netcdf() structs returned
 @geo_match_nc_structs.inc
@@ -249,28 +190,6 @@
 ; "include" file for PR data constants
 @pr_params.inc
 
-bname = file_basename( ncfilepr )
-prlen = strlen( bname )
-pctString = STRTRIM(STRING(FIX(pctabvthresh)),2)
-parsed = STRSPLIT( bname, '.', /extract )
-site = parsed[1]
-yymmdd = parsed[2]
-orbit = parsed[3]
-version = parsed[4]
-IF pr_or_dpr EQ 'DPR' THEN BEGIN
-   swath=parsed[6]
-   instrument='2A'+parsed[5]
-ENDIF ELSE BEGIN
-   swath='NS'
-   instrument='Ku'
-  ; leave this here for now, expect PR V08x version labels soon, though
-   CASE version OF
-        '6' : version = 'V6'
-        '7' : version = 'V7'
-       ELSE : print, "Using PR version = ", version
-   ENDCASE
-ENDELSE
-
 ; set up pointers for each field to be returned from fprep_geo_match_profiles()
 ptr_geometa=ptr_new(/allocate_heap)
 ptr_sweepmeta=ptr_new(/allocate_heap)
@@ -281,14 +200,29 @@
 ; define pointer for GR rain rate only if not using Z-R rainrate
 IF KEYWORD_SET(zr_force) EQ 0 THEN BEGIN
    IF pr_or_dpr EQ 'DPR' THEN BEGIN
-      ; only define the pointer specific to the rain rate field to be used
+      ; only define the pointers specific to the rain rate field to be used
       CASE gr_rr_field OF
-         'RC' : ptr_gvrc=ptr_new(/allocate_heap)
-         'RP' : ptr_gvrp=ptr_new(/allocate_heap)
-         'RR' : ptr_gvrr=ptr_new(/allocate_heap)
-         ELSE : ptr_gvrr=ptr_new(/allocate_heap)
+         'RC' : BEGIN
+                   ptr_gvrc=ptr_new(/allocate_heap)
+                   ptr_pctgoodrcgv=ptr_new(/allocate_heap)
+                END
+         'RP' : BEGIN
+                   ptr_gvrp=ptr_new(/allocate_heap)
+                   ptr_pctgoodrpgv=ptr_new(/allocate_heap)
+                END
+         'RR' : BEGIN
+                   ptr_gvrr=ptr_new(/allocate_heap)
+                   ptr_pctgoodrrgv=ptr_new(/allocate_heap)
+                END
+         ELSE : BEGIN
+                   ptr_gvrr=ptr_new(/allocate_heap)
+                   ptr_pctgoodrrgv=ptr_new(/allocate_heap)
+                END
       ENDCASE
-   ENDIF ELSE ptr_gvrr=ptr_new(/allocate_heap)     ; new for Version 2.2 PR matchup
+   ENDIF ELSE BEGIN
+      ptr_gvrr=ptr_new(/allocate_heap)          ; new for Version 2.2 PR matchup
+      ptr_pctgoodrrgv=ptr_new(/allocate_heap)
+   ENDELSE
 ENDIF
 
 ptr_BestHID=ptr_new(/allocate_heap)       ; new for Version 2.3 matchup file
@@ -303,6 +237,7 @@
 ptr_botm=ptr_new(/allocate_heap)
 ptr_lat=ptr_new(/allocate_heap)
 ptr_lon=ptr_new(/allocate_heap)
+ptr_pia=ptr_new(/allocate_heap)
 ptr_nearSurfRain=ptr_new(/allocate_heap)
 ptr_nearSurfRain_Comb=ptr_new(/allocate_heap)
 ptr_rnFlag=ptr_new(/allocate_heap)
@@ -316,7 +251,6 @@
 ptr_pctgoodpr=ptr_new(/allocate_heap)
 ptr_pctgoodgv=ptr_new(/allocate_heap)
 ptr_pctgoodrain=ptr_new(/allocate_heap)
-ptr_pctgoodrrgv=ptr_new(/allocate_heap)   ; new for Version 2.2 matchup file
 
 ; structure to hold bright band variables
 BBparms = {meanBB : -99.99, BB_HgtLo : -99, BB_HgtHi : -99}
@@ -340,7 +274,7 @@
     PTRGVZMEAN=ptr_gvz, PTRzcor=ptr_zcor, PTRzraw=ptr_zraw, PTRrain3d=ptr_rain3, $
     PTRGVRRMEAN=ptr_gvrr, PTRGVMODEHID=ptr_BestHID, PTRGVDZEROMEAN=ptr_GR_DP_Dzero, $
     PTRGVZDRMEAN=ptr_GR_DP_Zdr, PTRGVKDPMEAN=ptr_GR_DP_Kdp, PTRGVRHOHVMEAN=ptr_GR_DP_RHOhv, $
-    PTRtop=ptr_top, PTRbotm=ptr_botm, PTRlat=ptr_lat, PTRlon=ptr_lon, $
+    PTRtop=ptr_top, PTRbotm=ptr_botm, PTRlat=ptr_lat, PTRlon=ptr_lon, PTRpia=ptr_pia, $
     PTRsfcrainpr=ptr_nearSurfRain, PTRsfcraincomb=ptr_nearSurfRain_Comb, $
     PTRrainflag_int=ptr_rnFlag, PTRraintype_int=ptr_rnType, $
     PTRxCorners=ptr_xCorner, PTRyCorners=ptr_yCorner, PTRbbProx=ptr_bbProx, $
@@ -355,16 +289,18 @@
     PCT_ABV_THRESH=pctAbvThresh, S2KU=s2ku, PTRfieldflags=ptr_fieldflags, $
     PTRgeometa=ptr_geometa, PTRsweepmeta=ptr_sweepmeta, PTRsitemeta=ptr_sitemeta, $
     PTRGVZMEAN=ptr_gvz, PTRzcor=ptr_zcor, PTRzraw=ptr_zraw, PTRrain3d=ptr_rain3, $
-    PTRGVRRMEAN=ptr_gvrr, PTRGVMODEHID=ptr_BestHID, PTRGVDZEROMEAN=ptr_GR_DP_Dzero, $
+    PTRGVRCMEAN=ptr_gvrc, PTRGVRPMEAN=ptr_gvrp, PTRGVRRMEAN=ptr_gvrr, $
+    PTRGVMODEHID=ptr_BestHID, PTRGVDZEROMEAN=ptr_GR_DP_Dzero, $
     PTRGVZDRMEAN=ptr_GR_DP_Zdr, PTRGVKDPMEAN=ptr_GR_DP_Kdp, PTRGVRHOHVMEAN=ptr_GR_DP_RHOhv, $
-    PTRtop=ptr_top, PTRbotm=ptr_botm, PTRlat=ptr_lat, PTRlon=ptr_lon, $
+    PTRtop=ptr_top, PTRbotm=ptr_botm, PTRlat=ptr_lat, PTRlon=ptr_lon, PTRpia=ptr_pia ,$
     PTRsfcrainpr=ptr_nearSurfRain, PTRsfcraincomb=ptr_nearSurfRain_Comb, $
     PTRrainflag_int=ptr_rnFlag, PTRraintype_int=ptr_rnType, $
     PTRxCorners=ptr_xCorner, PTRyCorners=ptr_yCorner, PTRbbProx=ptr_bbProx, $
     PTRhgtcat=ptr_hgtcat, PTRdist=ptr_dist, PTRpridx_long=ptr_pr_index, $
     PTRpctgoodpr=ptr_pctgoodpr, PTRpctgoodgv=ptr_pctgoodgv, $
-    PTRpctgoodrain=ptr_pctgoodrain, PTRpctgoodrrgv=ptr_pctgoodrrgv, BBPARMS=BBparms, $
-    ALT_BB_HGT=alt_bb_hgt )
+    PTRpctgoodrain=ptr_pctgoodrain, PTRpctgoodrcgv=ptr_pctgoodrcgv, $
+    PTRpctgoodrpgv=ptr_pctgoodrpgv, PTRpctgoodrrgv=ptr_pctgoodrrgv, $
+    BBPARMS=BBparms, ALT_BB_HGT=alt_bb_hgt )
  END
 ENDCASE
 
@@ -389,49 +325,59 @@
   myflags=*ptr_fieldflags
     ptr_free,ptr_fieldflags
   gvz=*ptr_gvz
-  gvz_in = gvz     ; for plotting as PPI
     ptr_free,ptr_gvz
   zraw=*ptr_zraw
-  zraw_in = zraw   ; for plotting as PPI
     ptr_free,ptr_zraw
   zcor=*ptr_zcor
-  zcor_in = zcor   ; for plotting as PPI
     ptr_free,ptr_zcor
   rain3=*ptr_rain3
-  rain3_in = rain3 ; for plotting as PPI
-    ptr_free,ptr_rain3
+  ptr_free,ptr_rain3
+  have_pia=0
+  IF ptr_valid(ptr_pia) THEN BEGIN
+     pia=*ptr_pia
+     ptr_free,ptr_pia
+     IF pr_or_dpr EQ 'DPR' THEN have_pia=myflags.have_piaFinal $
+     ELSE have_pia=myflags.have_pia
+  ENDIF ELSE pia = -1
 
  ; initialize flag as to source of GR rain rate to use to "compute Z-R"
   have_gvrr = 0
+  gvrr = -1
+  pctgoodrrgv = -1
+  rr_field_used = 'Z-R'
 
   IF pr_or_dpr EQ 'DPR' THEN BEGIN
      CASE gr_rr_field OF
         'RC' : IF ptr_valid(ptr_gvrc) THEN BEGIN
                  gvrr=*ptr_gvrc
-                 gvrr_in = gvrr
                  ptr_free,ptr_gvrc
                  have_gvrr=myflags.have_GR_RC_rainrate
+                 IF ptr_valid(ptr_pctgoodrcgv) THEN pctgoodrrgv=*ptr_pctgoodrcgv
+                 rr_field_used = 'RC'
                ENDIF
         'RP' : IF ptr_valid(ptr_gvrp) THEN BEGIN
                  gvrr=*ptr_gvrp
-                 gvrr_in = gvrr
                  ptr_free,ptr_gvrp
                  have_gvrr=myflags.have_GR_RP_rainrate
-              ENDIF
+                 IF ptr_valid(ptr_pctgoodrpgv) THEN pctgoodrrgv=*ptr_pctgoodrpgv
+                 rr_field_used = 'RP'
+               ENDIF
         ELSE : IF ptr_valid(ptr_gvrr) THEN BEGIN
                  gvrr=*ptr_gvrr
-                 gvrr_in = gvrr
                  ptr_free,ptr_gvrr
                  have_gvrr=myflags.have_GR_RR_rainrate
+                 IF ptr_valid(ptr_pctgoodrrgv) THEN pctgoodrrgv=*ptr_pctgoodrrgv
+                 rr_field_used = 'RR'
                ENDIF
      ENDCASE
   ENDIF ELSE BEGIN
-    IF ptr_valid(ptr_gvrr) THEN BEGIN
-       gvrr=*ptr_gvrr
-       gvrr_in = gvrr ; for plotting as PPI
-       ptr_free,ptr_gvrr
-       have_gvrr=myflags.have_GR_rainrate   ; should just be 0 for version<2.2
-    ENDIF
+     IF ptr_valid(ptr_gvrr) THEN BEGIN
+        gvrr=*ptr_gvrr
+        ptr_free,ptr_gvrr
+        have_gvrr=myflags.have_GR_rainrate   ; should just be 0 for version<2.2
+        IF ptr_valid(ptr_pctgoodrrgv) THEN pctgoodrrgv=*ptr_pctgoodrrgv
+        rr_field_used = 'RR'
+     ENDIF
   ENDELSE
 
   haveHID = 0                          ; first check myflags values for all these?
@@ -439,35 +385,35 @@
      HIDcat=*ptr_BestHID
      haveHID = 1
      ptr_free,ptr_BestHID
-  ENDIF
+  ENDIF ELSE HIDcat=-1
 
   haveD0 = 0
   IF ptr_valid(ptr_GR_DP_Dzero) THEN BEGIN
      Dzero=*ptr_GR_DP_Dzero
      haveD0 = 1
      ptr_free,ptr_GR_DP_Dzero
-  ENDIF
+  ENDIF ELSE Dzero=-1
 
   haveZdr = 0
   IF ptr_valid(ptr_GR_DP_Zdr) THEN BEGIN
      Zdr=*ptr_GR_DP_Zdr
      haveZdr = 1
      ptr_free,ptr_GR_DP_Zdr
-  ENDIF
+  ENDIF ELSE Zdr=-1
 
   haveKdp = 0
   IF ptr_valid(ptr_GR_DP_Kdp) THEN BEGIN
      Kdp=*ptr_GR_DP_Kdp
      haveKdp = 1
      ptr_free,ptr_GR_DP_Kdp
-  ENDIF
+  ENDIF ELSE Kdp=-1
 
   haveRHOhv = 0
   IF ptr_valid(ptr_GR_DP_RHOhv) THEN BEGIN
      RHOhv=*ptr_GR_DP_RHOhv
      haveRHOhv = 1
      ptr_free,ptr_GR_DP_RHOhv
-  ENDIF
+  ENDIF ELSE RHOhv=-1
 
   top=*ptr_top
   botm=*ptr_botm
@@ -501,959 +447,149 @@
     ptr_free,ptr_pctgoodpr
     ptr_free,ptr_pctgoodgv
     ptr_free,ptr_pctgoodrain
-  IF ptr_valid(ptr_pctgoodrrgv) THEN BEGIN
-     pctgoodrrgv=*ptr_pctgoodrrgv
-     ptr_free,ptr_pctgoodrrgv
-  ENDIF
+;  IF ptr_valid(ptr_pctgoodrrgv) THEN BEGIN
+;     pctgoodrrgv=*ptr_pctgoodrrgv
+;     ptr_free,ptr_pctgoodrrgv
+;  ENDIF ELSE pctgoodrrgv=-1
 
-show_ppis=1   ; initialize to ON, override if PS and BATCH both are set
+; stuff the flags, structs, and data arrays into structures to pass along
+; - at this point, they will all be copies of the originals and we can
+;   butcher them as we please
 
- ; open a file to hold output stats to be appended to the Postscript file,
- ; if Postscript output is indicated
-  IF KEYWORD_SET( ps_dir ) THEN BEGIN
-     do_ps = 1
-     temptext = ps_dir + '/dbzdiffstats_temp.txt'
-     OPENW, tempunit, temptext, /GET_LUN
-    ; figure out whether to plot PPIs and prompt, or just write Postscript and 
-    ; wrap up
-     IF KEYWORD_SET( batch ) THEN show_ppis=0     ; don't display PPIs/animation
-  ENDIF ELSE do_ps = 0
+haveIt = { have_gvrr : have_gvrr, $
+           haveHID : haveHID, $
+           haveD0 : haveD0, $
+           haveZdr : haveZdr, $
+           haveKdp : haveKdp, $
+           haveRHOhv : haveRHOhv, $
+           have_pia : have_pia }
 
-nfp = mygeometa.num_footprints
-nswp = mygeometa.num_sweeps
-site_lat = mysite.site_lat
-site_lon = mysite.site_lon
-siteID = string(mysite.site_id)
-nsweeps = mygeometa.num_sweeps
+dataStruc = { haveFlags : haveIt, $
+              mygeometa : mygeometa, $
+              mysite : mysite, $
+              mysweeps : mysweeps, $
+              gvz : gvz, $
+              zraw : zraw, $
+              zcor : zcor, $
+              rain3 : rain3, $
+              gvrr : gvrr, $
+              rr_field_used : rr_field_used, $
+              HIDcat : HIDcat, $
+              Dzero : Dzero, $
+              Zdr : Zdr, $
+              Kdp : Kdp, $
+              RHOhv : RHOhv, $
+              top : top, $
+              botm : botm, $
+              lat : lat, $
+              lon : lon, $
+              pia : pia, $
+              rnflag : rnflag, $
+              rntype : rntype, $
+              pr_index : pr_index, $
+              xcorner : xcorner, $
+              ycorner : ycorner, $
+              bbProx : bbProx, $
+              dist : dist, $
+              hgtcat : hgtcat, $
+              pctgoodpr : pctgoodpr, $
+              pctgoodgv : pctgoodgv, $
+              pctgoodrain : pctgoodrain, $
+              pctgoodrrgv : pctgoodrrgv, $
+              BBparms : BBparms, $
+              heights : heights, $
+              hgtinterval : hgtinterval, $
+              is_subset : 0 }
 
-; override flag if forcing Z-R usage
-IF KEYWORD_SET(zr_force) THEN have_gvrr = 0
-print, ''
-IF (have_gvrr) THEN BEGIN
-   gvrr_txt = "Using GR RR field."
-   var2='RR'
-ENDIF ELSE BEGIN
-   gvrr_txt = "Using GR RR from Z-R."
-   var2='ZR'
-ENDELSE
-print, ''
+; - - - - - - - - - - - - - - - - - - - - - - -
 
-; For each 'column' of data, find the maximum GR reflectivity value for the
-;  footprint, and use this value to define a GR match to the PR-indicated rain type.
-;  Using Default GR dBZ thresholds of >=35 for "GV Convective" and <=25 for 
-;  "GV Stratiform", or other GR dBZ thresholds provided as user parameters,
-;  set PR rain type to "other" (3) where: PR type is Convective and GR is not, or
-;  PR is Stratiform and GR indicates Convective.  For GR reflectivities between
-;  'gvstratiform' and 'gvconvective' thresholds, leave the PR rain type as-is.
+IF N_ELEMENTS(submeth) EQ 1 THEN BEGIN
+  ; bring up the PPI location selector and cut out the area of interest
+   dataStrucTrimmed = select_geomatch_subarea( ncfilepr, hide_rntype, pr_or_dpr, $
+                                               startelev, dataStruc, $
+                                               SUBSET_METHOD=submeth, $
+                                               RR_OR_Z=xxx, $
+                                               RANGE_MAX=subthresh )
 
-print, ''
-max_gvz_per_fp = MAX( gvz, DIMENSION=2)
-IF ( gvstratiform GT 0.0 ) THEN BEGIN
-   idx2other = WHERE( rnType[*,0] EQ 2 AND max_gvz_per_fp LE gvstratiform, count2other )
-   IF ( count2other GT 0 ) THEN rnType[idx2other,*] = 3
-   fmtstrng='("No. of footprints switched from Convective to Other = ",I0,",' $
-            +' based on Stratiform dBZ threshold = ",F0.1)'
-   print, FORMAT=fmtstrng, count2other, gvstratiform
-ENDIF ELSE BEGIN
-   print, "Leaving PR Convective Rain Type assignments unchanged."
-ENDELSE
-IF ( gvconvective GT 0.0 ) THEN BEGIN
-   idx2other = WHERE( rnType[*,0] EQ 1 AND max_gvz_per_fp GE gvconvective, count2other )
-   IF ( count2other GT 0 ) THEN rnType[idx2other,*] = 3
-   fmtstrng='("No. of footprints switched from Stratiform to Other = ",I0,",' $
-            +' based on Convective dBZ threshold = ",F0.1)'
-   print, FORMAT=fmtstrng, count2other, gvconvective
-ENDIF ELSE BEGIN
-   print, "Leaving PR Stratiform Rain Type assignments unchanged."
-ENDELSE
-
-; make a copy of the adjusted rain type field for use in PPI plots
-rntype4ppi = REFORM( rnType[*,0] )
-; if rain type "hiding" is on, set all samples to "Other" rain type
-hide_rntype = KEYWORD_SET( hide_rntype )
-IF hide_rntype THEN rnType4ppi[*] = 3
-
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-; optional data *clipping* based on percent completeness of the volume averages:
-; Decide which PR and GR points to include, based on percent of expected points
-; in bin-averaged results that were above dBZ thresholds set when the matchups
-; were done.  If unspecified or set to zero, include all points regardless of
-; 'completeness' of the volume averages.
-
-
-IF ( pctAbvThresh GT 0.0 ) THEN BEGIN
-    ; clip to the 'good' points, where 'pctAbvThresh' fraction of bins in average
-    ; were above threshold
-   if have_gvrr AND xxx EQ 'RR' then begin
-      thresh_msg = "Thresholding by rain rate cutoff.  " + gvrr_txt
-      idxgoodenuff = WHERE( pctgoodrrgv GE pctAbvThresh $
-                       AND  pctgoodrain GE pctAbvThresh, countgoodpct )
-   endif else begin
-      IF xxx EQ 'Z' THEN thresh_msg = "Thresholding by reflectivity cutoffs." $
-      ELSE thresh_msg = "Thresholding by reflectivity cutoffs.  " + gvrr_txt
-      idxgoodenuff = WHERE( pctgoodpr GE pctAbvThresh $
-                       AND  pctgoodgv GE pctAbvThresh, countgoodpct )
-   endelse
-   IF ( countgoodpct GT 0 ) THEN BEGIN
-       IF have_gvrr EQ 0 THEN gvrr = z_r_rainrate(gvz[idxgoodenuff]) $  ; override empty field
-          ELSE gvrr=gvrr[idxgoodenuff]   ; using scaled GR rainrate from matchup file
-       gvz = gvz[idxgoodenuff]
-       zraw = zraw[idxgoodenuff]
-       zcor = zcor[idxgoodenuff]
-       rain3 = rain3[idxgoodenuff]
-       top = top[idxgoodenuff]
-       botm = botm[idxgoodenuff]
-       lat = lat[idxgoodenuff]
-       lon = lon[idxgoodenuff]
-       rnFlag = rnFlag[idxgoodenuff]
-       rnType = rnType[idxgoodenuff]
-       dist = dist[idxgoodenuff]
-       bbProx = bbProx[idxgoodenuff]
-       hgtcat = hgtcat[idxgoodenuff]
-;       pr_index = pr_index[idxgoodenuff] : NO! don't clip - must be full array for PPIs
-       IF ( PPIbyThresh ) THEN BEGIN
-           idx2plot=idxgoodenuff  ;idxpractual2d[idxgoodenuff]
-           n2plot=countgoodpct
-       ENDIF
+   IF size(dataStrucTrimmed, /TYPE) NE 8 THEN BEGIN
+      status = 1
+      message, "Unable to run statistics for storm area, skipping case.",/info
    ENDIF ELSE BEGIN
-       print, "No complete-volume points, quitting case."
-       goto, errorExit
-   ENDELSE
-ENDIF ELSE BEGIN
-  ; pctAbvThresh is 0, take/plot ALL non-bogus points
-   IF have_gvrr EQ 0 THEN gvrr = z_r_rainrate(gvz) $    ; override empty field
-      ELSE gvrr=gvrr   ; using scaled GR rainrate from matchup file
-   IF ( PPIbyThresh ) THEN BEGIN
-      idx2plot=WHERE( pctgoodpr GE 0.0 AND  pctgoodgv GE 0.0, countactual2d )
-      n2plot=countactual2d
-   ENDIF
-ENDELSE
-
-; we only use unclipped arrays for PPIs, so we make full copies of these arrays
-gvz_in2 = gvz_in
-zcor_in2 = zcor_in
-rain3_in2 = rain3_in
-IF have_gvrr EQ 0 THEN BEGIN
-   gvrr_in2 = z_r_rainrate(gvz_in)
-   gvrr_in2 = REFORM( gvrr_in2, nfp, nswp)
-ENDIF ELSE gvrr_in2 = gvrr_in
-
-; optional data *blanking* based on percent completeness of the volume averages
-; for PPI plots, operating on the full arrays of gvz and zcor
-
-IF ( PPIbyThresh ) THEN BEGIN
-   idx3d = LONG( gvz_in )   ; make a copy
-  ; re-set this for our later use in PPI plotting
-   idx3d[*,*] = 0L       ; initialize all points to 0
-   idx3d[idx2plot] = 2L  ; tag the points to be plotted in post-threshold PPI
-   idx2blank = WHERE( idx3d EQ 0L, count2blank )
-   IF ( count2blank GT 0 ) THEN BEGIN
-     gvz_in2[idx2blank] = 0.0
-     zcor_in2[idx2blank] = 0.0
-     rain3_in2[idx2blank] = 0.0
-     gvrr_in2[idx2blank] = 0.0
-   ENDIF
-ENDIF
-
-; determine the non-missing points-in-common between PR and GR, data value-wise,
-; to make sure the same points are plotted on PR and GR full/post-threshold PPIs
-idx2blank2 = WHERE( (gvz_in2 LT 0.0) OR (zcor_in2 LE 0.0) $
-                    OR (rain3_in2 LT 0.0) OR (gvrr_in2 LT 0.0), count2blank2 )
-IF ( count2blank2 GT 0 ) THEN BEGIN
-   gvz_in2[idx2blank2] = 0.0
-   zcor_in2[idx2blank2] = 0.0
-   rain3_in2[idx2blank2] = 0.0
-   gvrr_in2[idx2blank2] = 0.0
-ENDIF
-
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-; build an array of range categories from the GV radar, using ranges previously
-; computed from lat and lon by fprep_geo_match_profiles():
-; - range categories: 0 for 0<=r<50, 1 for 50<=r<100, 2 for r>=100, etc.
-;   (for now we only have points roughly inside 100km, so limit to 2 categs.)
-distcat = ( FIX(dist) / 50 ) < 1
-
-; get info from array of height category for the fixed-height levels, for profiles
-nhgtcats = N_ELEMENTS(heights)
-num_in_hgt_cat = LONARR( nhgtcats )
-FOR i=0, nhgtcats-1 DO BEGIN
-   hgtstr =  string(heights[i], FORMAT='(f0.1)')
-   idxhgt = where(hgtcat EQ i, counthgts)
-   num_in_hgt_cat[i] = counthgts
-ENDFOR
-
-; get info from array of BB proximity
-num_in_BB_Cat = LONARR(4)
-idxabv = WHERE( bbProx EQ 3, countabv )
-num_in_BB_Cat[3] = countabv
-idxblo = WHERE( bbProx EQ 1, countblo )
-num_in_BB_Cat[1] = countblo
-idxin = WHERE( bbProx EQ 2, countin )
-num_in_BB_Cat[2] = countin
-idxnobb = WHERE( bbProx EQ 0, countnobb )
-num_in_BB_Cat[0] = countnobb
-
-; build an array of sample volume depth for weighting of the layer averages and
-; mean differences
-voldepth = (top-botm) > 0.0
-
-;minz4hist = 18.  ; not used, replaced with cutoff
-maxz4hist = 55.   ; used only for Z
-IF N_ELEMENTS(max_range) NE 1 THEN rangecut = 100. ELSE rangecut = max_range
-
-; define a structure to hold difference statistics computed within and returned
-; by the called function calc_geo_pr_gv_meandiffs_wght_idx()
-the_struc = { diffsset, $
-              meandiff: -99.999, meandist: -99.999, fullcount: 0L,  $
-              maxpr: -99.999, maxgv: -99.999, $
-              meandiffc: -99.999, countc: 0L, $
-              meandiffs: -99.999, counts: 0L, $
-              AvgDifByHist: -99.999 $
-            }
-
-print, ''
-
-CASE xxx OF
-   'RR' : BEGIN
-          difftext='-GR Rain Rate difference statistics (mm/h) - GR Site: '
-          cutoff = 0.1     ; cutoff RR value to include in mean diff. calcs.
-          yvar = rain3
-          xvar = gvrr
-          src1 = pr_or_dpr
-          src2 = 'GR'
-;          field = var2
-          END
-    'Z' : BEGIN
-          difftext='-GR Reflectivity difference statistics (dBZ) - GR Site: '
-          cutoff = 10.     ; cutoff dBZ value to include in mean diff. calcs.
-          bs = 2.0         ; fixed histogram bin size, else bins by 1.0, messy
-          src1 = pr_or_dpr
-          src2 = 'GR'
-          yvar = zcor
-          xvar = gvz
-;          field = xxx
-          END
-ENDCASE
-
-CASE pr_or_dpr OF
-  'PR' : textout = pr_or_dpr + difftext + siteID + '   Orbit: ' + orbit + $
-                   '  Version: '+version
-  'DPR' : BEGIN
-          textout = pr_or_dpr + ' ' + Instrument + difftext + siteID
-          print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-          textout = 'Orbit: '+orbit+'  Version: '+version+'  Swath Type: '+swath
-          END
-ENDCASE
-print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-
-textout = pr_or_dpr + ' time = ' + mygeometa.atimeNearestApproach + $
-          '   GR start time = ' + mysweeps[0].atimeSweepStart
-print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-textout = 'Required percent of above-threshold ' + pr_or_dpr + $
-          ' and GR bins in matched volumes >= '+pctString+"%"
-print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-IF ( pctAbvThresh GT 0.0 ) THEN BEGIN
-   print, thresh_msg & IF (do_ps EQ 1) THEN printf, tempunit, thresh_msg
-ENDIF
-IF ( s2ku ) THEN BEGIN
-   textout = 'GR reflectivity has S-to-Ku frequency adjustments applied.'
-   print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-ENDIF
-
-mnprarr = fltarr(3,nhgtcats)  ; each level, for raintype all, stratiform, convective
-mngvarr = fltarr(3,nhgtcats)  ; each level, for raintype all, stratiform, convective
-levhasdata = intarr(nhgtcats) & levhasdata[*] = 0
-levsdata = 0
-max_hgt_w_data = 0.0
-
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-;print_table_headers, src1, src2, field, PS_UNIT=tempunit
-print_table_headers, src1, src2, xxx, PS_UNIT=tempunit
-
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-; Compute a mean rainrate difference at each level
-
-for lev2get = 0, nhgtcats-1 do begin
-   havematch = 0
-   thishgt = (lev2get+1)*hgtinterval
-   IF ( num_in_hgt_cat[lev2get] GT 0 ) THEN BEGIN
-      flag = ''
-      idx4hist = lonarr(num_in_hgt_cat[lev2get])  ; array indices used for point-to-point mean diffs
-      idx4hist[*] = -1L
-      if (lev2get eq BBparms.BB_HgtLo OR lev2get eq BBparms.BB_HgtHi) then flag = ' @ BB'
-      diffstruc = the_struc
-      calc_geo_pr_gv_meandiffs_wght_idx, yvar, xvar, rnType, dist, distcat, $
-                                         hgtcat, lev2get, cutoff, rangecut, $
-                                         mnprarr, mngvarr, havematch, $
-                                         diffstruc, idx4hist, voldepth
-      if(havematch eq 1) then begin
-         levsdata = levsdata + 1
-         levhasdata[lev2get] = 1
-         max_hgt_w_data = thishgt
-        ; format level's stats for table output
-         FMT55 = '(3("    ",f7.3,"    ",i4),"  ",3("   ",f7.3))'
-         stats55 = STRING(diffstruc.meandiff, diffstruc.fullcount, $
-                          diffstruc.meandiffs, diffstruc.counts, $
-                          diffstruc.meandiffc, diffstruc.countc, $
-                          diffstruc.meandist, diffstruc.maxpr, diffstruc.maxgv, $
-                          FORMAT = FMT55 )
-        ; extract/format level's stats for graphic plots output
-         mndifhstr = string(diffstruc.AvgDifByHist, FORMAT='(f0.3)')
-         mndifstr = string(diffstruc.meandiff, diffstruc.fullcount, FORMAT='(f0.3," (",i0,")")')
-         mndifstrc = string(diffstruc.meandiffc, diffstruc.countc, FORMAT='(f0.3," (",i0,")")')
-         mndifstrs = string(diffstruc.meandiffs, diffstruc.counts, FORMAT='(f0.3," (",i0,")")')
-         idx4hist[*] = -1L
-         textout = STRING(heights[lev2get], stats55, flag, FORMAT='(" ",f4.1,a0," ",a0)')
-         print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-      endif else begin
-         textout = "No above-threshold points at height " + $
-                   STRING(heights[lev2get], FORMAT='(f0.3)')
-         print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-      endelse
-   ENDIF ELSE BEGIN
-      print, "No points at height " + string(heights[lev2get], FORMAT='(f0.3)')
-   ENDELSE
-
-endfor         ; lev2get = 0, nhgtcats-1
-
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-; Build the mean RR profile plot panel
-
-orig_device = !D.NAME
-CASE xxx OF
-   'RR' : IF (have_gvrr) THEN gr_rr_zr = ' DP RR' ELSE gr_rr_zr = ' Z-R RR'
-    'Z' : gr_rr_zr = ' Zc'
-ENDCASE
-
-IF ( do_ps EQ 1 ) THEN BEGIN
-  ; set up postscript plot params. and file path/name
-   cd, ps_dir
-   b_w = keyword_set(b_w)
-   IF ( s2ku ) THEN add2nm = '_S2Ku' ELSE add2nm = ''
-   PSFILEpdf = ps_dir+'/'+site+'.'+yymmdd+'.'+orbit+"."+version+'.'+pr_or_dpr+'_' $
-               +instrument+'_'+swath+'.Pct'+pctString+add2nm+'_'+xxx+'_PDF_SCATR.ps'
-   print, "Output sent to ", PSFILEpdf
-   set_plot,/copy,'ps'
-   device,filename=PSFILEpdf,/color,bits=8,/inches,xoffset=0.25,yoffset=2.55, $
-          xsize=8.,ysize=8.
-
-   ; Set up color table
-   ;
-   common colors, r_orig, g_orig, b_orig, r_curr, g_curr, b_curr ;load color table
-   IF ( b_w EQ 0) THEN  LOADCT, 6  ELSE  LOADCT, 33
-   ncolor=255
-   red=bytarr(256) & green=bytarr(256) & blue=bytarr(256)
-   red=r_curr & green=g_curr & blue=b_curr
-   red(0)=255 & green(0)=255 & blue(0)=255
-   red(1)=115 & green(1)=115 & blue(1)=115  ; gray for GR
-   red(ncolor)=0 & green(ncolor)=0 & blue(ncolor)=0 
-   tvlct,red,green,blue
-   !P.COLOR=0 ; make the title and axis annotation black
-   !X.THICK=2 ; make the ticks and borders thicker
-   !Y.THICK=2 ; ditto
-   !P.FONT=0 ; use the device fonts supplied by postscript
-
-   IF ( b_w EQ 0) THEN BEGIN
-     PR_COLR=200
-     GV_COLR=60
-     ST_LINE=1    ; dotted for stratiform
-     CO_LINE=2    ; dashed for convective
-   ENDIF ELSE BEGIN
-     PR_COLR=ncolor
-     GV_COLR=ncolor
-     ST_LINE=0    ; solid for stratiform
-     CO_LINE=1    ; dotted for convective
-   ENDELSE
-
-   CHARadj=0.75
-   THIKadjPR=1.5
-   THIKadjGV=1.5
-   ST_THK=1
-   CO_THK=1
-ENDIF ELSE BEGIN
-  ; set up x-window plot params.
-   device, decomposed = 0
-   LOADCT, 2
-
-   IF ( pctAbvThresh GT 0.0 ) THEN BEGIN
-      IF ( pctAbvThresh EQ 100.0 ) THEN gt_ge = ' ' ELSE gt_ge = ' >='
-      CASE xxx OF
-         'Z' : wintxt = "With" + gt_ge + pctString + $
-                        "% of averaged bins above dBZ thresholds"  
-        'RR' : wintxt = "With" + gt_ge + pctString + $
-                        "% of averaged bins above rainrate threshold"
-      ENDCASE
-   ENDIF ELSE BEGIN
-      wintxt = "With all non-missing "+pr_or_dpr+"/GR matched samples"
-   ENDELSE
-
-   Window, xsize=700, ysize=700, TITLE = site+gr_rr_zr+' vs. '+instrument+'.'+ $
-           swath+"."+version+"  --  "+wintxt, RETAIN=2
-   PR_COLR=30
-   GV_COLR=70
-   ST_LINE=1    ; dotted for stratiform
-   CO_LINE=2    ; dashed for convective
-   CHARadj=1.0
-   THIKadjPR=1.0
-   THIKadjGV=1.0
-   ST_THK=3
-   CO_THK=2
-ENDELSE
-
-
-!P.Multi=[0,2,2,0,0]
-
-if (levsdata eq 0) then begin
-   print, "No valid data levels found!"
-   nframes = 0
-   goto, nextFile
-endif
-
-idxlev2plot = WHERE( levhasdata EQ 1 )
-h2plot = heights[idxlev2plot]
-
-; figure out the y-axis range.  Use the greater of max_hgt_w_data*2.0
-; and meanbb*2 as the proposed range.  Cut off at 20 km if result>20.
-prop_max_y = max_hgt_w_data*2.0 > (FIX((BBparms.meanbb*2)/1.5) + 1) *1.5
-
-CASE xxx OF
-   'Z' : BEGIN
-           plot, [15,50], [0,20 < prop_max_y], /NODATA, COLOR=255, $
-             XSTYLE=1, YSTYLE=1, YTICKINTERVAL=hgtinterval, YMINOR=1, thick=1, $
-             XTITLE='Level Mean Reflectivity, dBZ', YTITLE='Height Level, km', $
-             CHARSIZE=1*CHARadj, BACKGROUND=0
-           xvals = [15,50]   ; x-endpoints of BB line, dBZ scale
-         END
-  'RR' : BEGIN
-           plot, [0.1,150], [0,20 < prop_max_y], /NODATA, COLOR=255, $
-             XSTYLE=1, YSTYLE=1, YTICKINTERVAL=hgtinterval, YMINOR=1, thick=1, $
-             XTITLE='Level Mean Rain Rate, mm/h', YTITLE='Height Level, km', $
-             CHARSIZE=1*CHARadj, BACKGROUND=0, /xlog
-           xvals = [0.1,150]   ; x-endpoints of BB line, RR scale
-         END
-ENDCASE
-
-IF (~ hideTotals) THEN BEGIN
-  ; plot the profile for all points regardless of rain type
-   prmnz2plot = mnprarr[0,*]
-   prmnz2plot = prmnz2plot[idxlev2plot]
-   gvmnz2plot = mngvarr[0,*]
-   gvmnz2plot = gvmnz2plot[idxlev2plot]
-   oplot, prmnz2plot, h2plot, COLOR=PR_COLR, thick=1*THIKadjPR
-   oplot, gvmnz2plot, h2plot, COLOR=GV_COLR, thick=1*THIKadjGV
-ENDIF
-
-; plot the profile for stratiform rain type points
-prmnz2plot = mnprarr[1,*]
-prmnz2plot = prmnz2plot[idxlev2plot]
-gvmnz2plot = mngvarr[1,*]
-gvmnz2plot = gvmnz2plot[idxlev2plot]
-idxhavezs = WHERE( prmnz2plot GT 0.0 and gvmnz2plot GT 0.0, counthavezs )
-IF ( counthavezs GT 0 ) THEN BEGIN
-   oplot, prmnz2plot[idxhavezs], h2plot[idxhavezs], COLOR=PR_COLR, LINESTYLE=ST_LINE, thick=3*THIKadjPR
-   oplot, gvmnz2plot[idxhavezs], h2plot[idxhavezs], COLOR=GV_COLR, LINESTYLE=ST_LINE, thick=3*THIKadjGV
-ENDIF
-
-; plot the profile for convective rain type points
-prmnz2plot = mnprarr[2,*]
-prmnz2plot = prmnz2plot[idxlev2plot]
-gvmnz2plot = mngvarr[2,*]
-gvmnz2plot = gvmnz2plot[idxlev2plot]
-idxhavezs = WHERE( prmnz2plot GT 0.0 and gvmnz2plot GT 0.0, counthavezs )
-IF ( counthavezs GT 0 ) THEN BEGIN
-   oplot, prmnz2plot[idxhavezs], h2plot[idxhavezs], COLOR=PR_COLR, LINESTYLE=CO_LINE, thick=3*THIKadjPR
-   oplot, gvmnz2plot[idxhavezs], h2plot[idxhavezs], COLOR=GV_COLR, LINESTYLE=CO_LINE, thick=3*THIKadjGV
-ENDIF
-
-; plot the mean BB indicator line and its legend info.
-yvalsbb = [BBparms.meanbb, BBparms.meanbb]
-plots, xvals, yvalsbb, COLOR=255, LINESTYLE=2;, THICK=3*THIKadjGV
-plots, [0.29,0.33], [0.805,0.805], COLOR=255, /NORMAL, LINESTYLE=2
-XYOutS, 0.34, 0.8, 'Mean BB Hgt', COLOR=255, CHARSIZE=1*CHARadj, /NORMAL
-
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-; Compute a mean rainrate difference at each BB proximity layer and plot PDFs
-
-mnprarrbb = fltarr(3,4)  ; each BB-relative level, for raintype all, stratiform, convective
-mngvarrbb = fltarr(3,4)  ; each BB-relative level, for raintype all, stratiform, convective
-levhasdatabb = intarr(4) & levhasdatabb[*] = 0
-levsdatabb = 0
-bblevstr = ['Unknown', ' Below', 'Within', ' Above']
-xoff = [0.0, 0.0, -0.5, 0.0 ]  ; for positioning legend in PDFs
-yoff = [0.0, 0.0, -0.5, -0.5 ]
-
-;print_table_headers, src1, src2, field, /BB, PS_UNIT=tempunit
-print_table_headers, src1, src2, xxx, /BB, PS_UNIT=tempunit
-
-; define a 2D array to capture array indices of values used for point-to-point
-; mean diffs, for each of the 3 bbProx levels, for plotting these points in the
-; scatter plots
-numZpts = N_ELEMENTS(yvar)
-idx4hist3 = lonarr(3,numZpts)
-idx4hist3[*,*] = -1L
-num4hist3 = lonarr(3)  ; track number of points used in each bbProx layer
-idx4hist = idx4hist3[0,*]  ; a 1-D array for passing to function in the layer loop
-
-; set up for known versus unknown BB proximity case
-IF bbparms.meanBB NE -99.99 THEN BEGIN
-   bblevBeg = 1
-   bblevEnd = 3
-   pmultifac = 4
-   pmultirows = 2
-ENDIF ELSE BEGIN
-   bblevBeg = 0
-   bblevEnd = 0
-   pmultifac = 3
-   pmultirows = 2
-ENDELSE
-
-for bblev2get = bblevBeg, bblevEnd do begin
-   havematch = 0
-   !P.Multi=[pmultifac-bblev2get,2,pmultirows,0,0]
-   IF ( num_in_BB_cat[bblev2get] GT 0 ) THEN BEGIN
-      flag = ''
-      if (bblev2get eq 2) then flag = ' @ BB'
-      diffstruc = the_struc
-      calc_geo_pr_gv_meandiffs_wght_idx, yvar, xvar, rnType, dist, distcat, $
-                                         bbProx, bblev2get, cutoff, rangecut, $
-                                         mnprarrbb, mngvarrbb, havematch, $
-                                         diffstruc, idx4hist, voldepth
-      if(havematch eq 1) then begin
-         levsdatabb = levsdatabb + 1
-         levhasdatabb[bblev2get] = 1
-        ; format level's stats for table output
-         FMT55='("  ",f7.3,"    ",i4,2("    ",f7.3,"    ",i4),"  ",3("   ",f7.3))'
-         stats55 = STRING(diffstruc.meandiff, diffstruc.fullcount, $
-                       diffstruc.meandiffs, diffstruc.counts, $
-                       diffstruc.meandiffc, diffstruc.countc, $
-                       diffstruc.meandist, diffstruc.maxpr, diffstruc.maxgv, $
-                       FORMAT = FMT55 )
-        ; capture points used, and format level's stats for graphic plots output
-         num4hist3[bblev2get-bblevBeg] = diffstruc.fullcount
-         idx4hist3[bblev2get-bblevBeg,*] = idx4hist
-         rr_pr2 = yvar[idx4hist[0:diffstruc.fullcount-1]]
-         rr_gv2 = xvar[idx4hist[0:diffstruc.fullcount-1]]
-         type2 = rnType[idx4hist[0:diffstruc.fullcount-1]]
-         bbProx2 = bbProx[idx4hist[0:diffstruc.fullcount-1]]
-         mndifhstr = string(diffstruc.AvgDifByHist, FORMAT='(f0.3)')
-         mndifstr = string(diffstruc.meandiff, diffstruc.fullcount, $
-                           FORMAT='(f0.3," (",i0,")")')
-         IF diffstruc.countc EQ 0 THEN mndifstrc = 'None' $
-         ELSE mndifstrc = string(diffstruc.meandiffc, diffstruc.countc, $
-                                 FORMAT='(f0.3," (",i0,")")')
-         IF diffstruc.counts EQ 0 THEN mndifstrs = 'None' $
-         ELSE mndifstrs = string(diffstruc.meandiffs, diffstruc.counts, $
-                                 FORMAT='(f0.3," (",i0,")")')
-         idx4hist[*] = -1
-         textout = STRING(bblevstr[bblev2get], stats55, flag, $
-                          FORMAT='(a0," ",a0," ",a0)')
-         print, textout & IF (do_ps EQ 1) THEN printf, tempunit, textout
-
-        ; Plot the PDF graph for this level
-         hgtline = 'Layer = ' + bblevstr[bblev2get] + " BB"
-
-        ; DO ANY/ALL RAINTYPE PDFS FIRST
-         CASE xxx OF
-         'RR' : BEGIN
-               ; define a set of "nlogcats" log-spaced interval boundaries
-               ; - yields nlogcats+1 rainrate categories
-                nlogcats = 16
-                logbins = 10^(findgen(nlogcats)/5.-1)
-               ; figure out index of interval where each point falls:
-               ; -- ranges from -1 (below lowest bound) to nlogcats-1 (above highest bound)
-                bin4pr = VALUE_LOCATE( logbins, rr_pr2 )
-                bin4gr = VALUE_LOCATE( logbins, rr_gv2 )  ; ditto for GR rainrate
-               ; compute histogram of log range category, ignoring the lowest (below 0.1 mm/h)
-                prhist = histogram( bin4pr, min=0, max=nlogcats-1,locations = prhiststart )
-                nxhist = histogram( bin4gr, min=0, max=nlogcats-1,locations = prhiststart )
-               ; will label every other interval start value on plot, no room for all
-                labelbins=['0.10','0.25','0.63','1.58','3.98','10.0','25.1','63','158']
-                plot, [0,MAX(prhiststart)],[0,FIX((MAX(prhist)>MAX(nxhist))*1.1)], $
-                      /NODATA, COLOR=255, CHARSIZE=1*CHARadj, $
-                      XTITLE=bblevstr[bblev2get]+' BB Rain Rate, mm/h', $
-                      YTITLE='Number of PR Footprints', $
-                      YRANGE=[ 0, FIX((MAX(prhist)>MAX(nxhist))*1.1) + 1 ], $
-                      BACKGROUND=0, xtickname=labelbins , $
-                      xtickinterval=2,xminor=2
-                END
-          'Z' : BEGIN
-                prhist = histogram(rr_pr2, min=cutoff, max=maxz4hist, binsize = bs, $
-                                   locations = prhiststart)
-                nxhist = histogram(rr_gv2, min=cutoff, max=maxz4hist, binsize = bs)
-                plot, [15,MAX(prhiststart)],[0,FIX((MAX(prhist)>MAX(nxhist))*1.1)], $
-                      /NODATA, COLOR=255, $
-                      XTITLE=bblevstr[bblev2get]+' BB Reflectivity, dBZ', $
-                      YTITLE='Number of DPR Footprints', $
-                      YRANGE=[ 0, FIX((MAX(prhist)>MAX(nxhist))*1.1) + 1 ], $
-                      BACKGROUND=0
-                END
-         ENDCASE
-         IF ( ~ hideTotals ) THEN BEGIN
-            oplot, prhiststart, prhist, COLOR=PR_COLR
-            oplot, prhiststart, nxhist, COLOR=GV_COLR
-            xyouts, 0.34, 0.95, pr_or_dpr+' (all)', COLOR=PR_COLR, /NORMAL, $
-                    CHARSIZE=1*CHARadj
-            plots, [0.29,0.33], [0.955,0.955], COLOR=PR_COLR, /NORMAL
-            xyouts, 0.34, 0.925, siteID+' (all)', COLOR=GV_COLR, /NORMAL, $
-                    CHARSIZE=1*CHARadj
-            plots, [0.29,0.33], [0.93,0.93], COLOR=GV_COLR, /NORMAL
+      status = render_rr_or_z_plots(ncfilepr, xxx, looprate, elevs2show, startelev, $
+                                    PPIorient, windowsize, pctabvthresh, PPIbyThresh, $
+                                    gvconvective, gvstratiform, hideTotals, $
+                                    hide_rntype, pr_or_dpr, dataStrucTrimmed, $
+                                    PS_DIR=ps_dir, B_W=b_w, S2KU=s2ku, ZR=zr_force, $
+                                    BATCH=batch, MAX_RANGE=max_range, $
+                                    STEP_MANUAL=step_manual )
+      wdelete, 1
+      saveIt=0
+      IF ( N_ELEMENTS(save_dir) EQ 1 ) THEN BEGIN
+         doodah = ""
+         PRINT, STRING(7B)  ; ring the terminal bell
+         WHILE (doodah NE 'Y' AND doodah NE 'N') DO BEGIN
+            READ, doodah, PROMPT='Save subset variables to file?  Enter Y or N : '
+            doodah = STRTRIM(STRUPCASE(doodah),2)
+            CASE doodah OF
+              'Y' : saveIt=1
+              'N' : saveIt=0
+             ELSE : BEGIN
+                      PRINT, STRING(7B)
+                      PRINT, "Illegal response, enter Y or N."
+                    END
+            ENDCASE
+         ENDWHILE
+      ENDIF
+      IF ( saveIt ) THEN BEGIN
+        ; set up IDL SAVE file path/name
+         bname = file_basename( ncfilepr )
+         pctString = STRTRIM(STRING(FIX(pctabvthresh)),2)
+         parsed = STRSPLIT( bname, '.', /extract )
+         site = parsed[1]
+         yymmdd = parsed[2]
+         orbit = parsed[3]
+         version = parsed[4]
+         IF pr_or_dpr EQ 'DPR' THEN BEGIN
+            swath=parsed[6]
+            instrument='2A'+parsed[5]
+         ENDIF ELSE BEGIN
+            swath='NS'
+            instrument='Ku'
+           ; leave this here for now, expect PR V08x version labels soon, though
+            CASE version OF
+                 '6' : version = 'V6'
+                 '7' : version = 'V7'
+                ELSE : print, "Using PR version = ", version
+            ENDCASE
+         ENDELSE
+         IF ( s2ku ) THEN add2nm = '_S2Ku' ELSE add2nm = ''
+         IF datastrucTrimmed.is_subset THEN BEGIN
+           ; format the storm lat/lon position into a string to be added to the PS name
+            IF datastrucTrimmed.storm_lat LT 0.0 THEN hemi='S' ELSE hemi='N'
+            IF datastrucTrimmed.storm_lon LT 0.0 THEN ew='W' ELSE ew='E'
+            addpos='_'+STRING(ABS(datastrucTrimmed.storm_lat),FORMAT='(f0.2)')+hemi+'_'+ $
+                   STRING(ABS(datastrucTrimmed.storm_lon),FORMAT='(f0.2)')+ew
+            add2nm = add2nm+addpos
          ENDIF
-
-         headline = pr_or_dpr+'-'+siteID+' Biases:'
-         xyouts, 0.775+xoff[bblev2get],0.925+yoff[bblev2get], headline, $
-                 COLOR=255, /NORMAL, CHARSIZE=1*CHARadj
-
-         mndifline = 'Any/All: ' + mndifstr
-         mndiflinec = 'Convective: ' + mndifstrc
-         mndiflines = 'Stratiform: ' + mndifstrs
-         mndifhline = 'By Area Mean: ' + mndifhstr
-         xyouts, 0.775+xoff[bblev2get],0.9+yoff[bblev2get], mndifline, $
-                 COLOR=255, /NORMAL, CHARSIZE=1*CHARadj
-         xyouts, 0.775+xoff[bblev2get],0.875+yoff[bblev2get], mndiflinec, $
-                 COLOR=255, /NORMAL, CHARSIZE=1*CHARadj
-         xyouts, 0.775+xoff[bblev2get],0.85+yoff[bblev2get], mndiflines, $
-                 COLOR=255, /NORMAL, CHARSIZE=1*CHARadj
-
-        ; OVERLAY CONVECTIVE RAINTYPE PDFS, IF ANY POINTS
-         idxconvhist= WHERE( type2 EQ RainType_convective, nconv )
-         IF ( nconv GT 0 ) THEN BEGIN
-           CASE xxx OF
-           'RR' : BEGIN
-                  bin4pr = VALUE_LOCATE( logbins, rr_pr2[idxconvhist] )
-                  bin4gr = VALUE_LOCATE( logbins, rr_gv2[idxconvhist] )
-                  prhist = histogram( bin4pr, min=0, max=18,locations = prhiststart )
-                  nxhist = histogram( bin4gr, min=0, max=18,locations = prhiststart )
-                  END
-            'Z' : BEGIN
-                  prhist = histogram(rr_pr2[idxconvhist], min=cutoff, max=maxz4hist, $
-                                     binsize = bs, locations = prhiststart)
-                  nxhist = histogram(rr_gv2[idxconvhist], min=cutoff, max=maxz4hist, $
-                                     binsize = bs)
-                  END
-           ENDCASE
-           oplot, prhiststart, prhist, COLOR=PR_COLR, LINESTYLE=CO_LINE, $
-                  thick=3*THIKadjPR
-           oplot, prhiststart, nxhist, COLOR=GV_COLR, LINESTYLE=CO_LINE, $
-                  thick=3*THIKadjGV
-           xyouts, 0.34, 0.85, pr_or_dpr+' (Conv)', COLOR=PR_COLR, /NORMAL, $
-                   CHARSIZE=1*CHARadj
-           xyouts, 0.34, 0.825, siteID+' (Conv)', COLOR=GV_COLR, /NORMAL, $
-                   CHARSIZE=1*CHARadj
-           plots, [0.29,0.33], [0.855,0.855], COLOR=PR_COLR, /NORMAL, $
-                  LINESTYLE=CO_LINE, thick=3*THIKadjPR
-           plots, [0.29,0.33], [0.83,0.83], COLOR=GV_COLR, /NORMAL, $
-                  LINESTYLE=CO_LINE, thick=3*THIKadjGV
-         ENDIF
-
-        ; OVERLAY STRATIFORM RAINTYPE PDFS, IF ANY POINTS
-         idxstrathist= WHERE( type2 EQ RainType_stratiform, nstrat )
-         IF ( nstrat GT 0 ) THEN BEGIN
-           CASE xxx OF
-           'RR' : BEGIN
-                  bin4pr = VALUE_LOCATE( logbins, rr_pr2[idxstrathist] )
-                  bin4gr = VALUE_LOCATE( logbins, rr_gv2[idxstrathist] )
-                  prhist = histogram( bin4pr, min=0, max=18,locations = prhiststart )
-                  nxhist = histogram( bin4gr, min=0, max=18,locations = prhiststart )
-                  END
-            'Z' : BEGIN
-                  prhist = histogram(rr_pr2[idxstrathist], min=cutoff, max=maxz4hist, $
-                                     binsize = bs, locations = prhiststart)
-                  nxhist = histogram(rr_gv2[idxstrathist], min=cutoff, max=maxz4hist, $
-                                     binsize = bs)
-                  END
-           ENDCASE
-           oplot, prhiststart, prhist, COLOR=PR_COLR, LINESTYLE=ST_LINE, $
-                  thick=3*THIKadjPR
-           oplot, prhiststart, nxhist, COLOR=GV_COLR, LINESTYLE=ST_LINE, $
-                  thick=3*THIKadjGV
-           xyouts, 0.34, 0.9, pr_or_dpr+' (Strat)', COLOR=PR_COLR, /NORMAL, $
-                   CHARSIZE=1*CHARadj
-           xyouts, 0.34, 0.875, siteID+' (Strat)', COLOR=GV_COLR, /NORMAL, $
-                   CHARSIZE=1*CHARadj
-           plots, [0.29,0.33], [0.905,0.905], COLOR=PR_COLR, /NORMAL, $
-                  LINESTYLE=ST_LINE, thick=3*THIKadjPR
-           plots, [0.29,0.33], [0.88,0.88], COLOR=GV_COLR, /NORMAL, $
-                  LINESTYLE=ST_LINE, thick=3*THIKadjGV
-         ENDIF
- 
-      endif else begin
-         print, "No above-threshold points ", bblevstr[bblev2get], " Bright Band"
-      endelse
-   ENDIF ELSE BEGIN
-      print, "No points at proximity = ", bblevstr[bblev2get], " Bright Band"
-      xyouts, 0.6+xoff[bblev2get],0.75+yoff[bblev2get], bblevstr[bblev2get] + $
-              " BB: NO POINTS", COLOR=255, /NORMAL, CHARSIZE=1.5
+         SAVFILE = save_dir+'/'+site+'.'+yymmdd+'.'+orbit+"."+version+'.'+pr_or_dpr+'_' $
+                     +instrument+'_'+swath+'.Pct'+pctString+add2nm+'_'+xxx+'.sav'
+         SAVE, ncfilepr, xxx, looprate, elevs2show, startelev, PPIorient, windowsize, $
+               pctabvthresh, PPIbyThresh, gvconvective, gvstratiform, hideTotals, $
+               hide_rntype, pr_or_dpr, datastrucTrimmed, FILE=SAVFILE
+         print, "Data saved to ", SAVFILE
+      ENDIF
    ENDELSE
-
-endfor         ; bblev2get = 1,3
-
-IF ( s2ku ) THEN xyouts, 0.29, 0.775, '('+siteID+' Ku-adjusted)', COLOR=GV_COLR, $
-                 /NORMAL, CHARSIZE=1*CHARadj
-
-; Write a data identification line at the bottom of the page below the PDF
-; plots for Postscript output.  This line also goes at the top of the scatter
-; plots, hence the name.
-
-;CASE xxx OF
-;   'RR' : IF (have_gvrr) THEN gr_rr_zr = ' DP RR' ELSE gr_rr_zr = ' Z-R RR'
-;    'Z' : gr_rr_zr = ' Zc'
-;ENDCASE
-IF ( s2ku ) THEN kutxt=' Ku-adjusted ' else kutxt=''
-IF ( pctAbvThresh GT 0.0 ) THEN BEGIN
-   IF ( pctAbvThresh EQ 100.0 ) THEN gt_ge = "     " ELSE gt_ge = "    >="
-   SCATITLE = site+kutxt+gr_rr_zr+' vs. '+pr_or_dpr+' '+instrument+'/'+swath+"/" $
-              +version+gt_ge+pctString+"% bins above threshold"
 ENDIF ELSE BEGIN
-   SCATITLE = site+kutxt+gr_rr_zr+' vs. '+pr_or_dpr+' '+instrument+'/'+swath+"/"+version $
-              +" -- All non-missing pairs"
+   ; call the routine to produce the graphics and output, just doing entire area
+   status = render_rr_or_z_plots(ncfilepr, xxx, looprate, elevs2show, startelev, $
+                                 PPIorient, windowsize, pctabvthresh, PPIbyThresh, $
+                                 gvconvective, gvstratiform, hideTotals, $
+                                 hide_rntype, pr_or_dpr, dataStruc, PS_DIR=ps_dir, $
+                                 B_W=b_w, S2KU=s2ku, ZR=zr_force, BATCH=batch, $
+                                 MAX_RANGE=max_range, STEP_MANUAL=step_manual )
 ENDELSE
 
-TITLE2 = "Orbit:  "+orbit+"  --  GR Start Time:  "+mysweeps[0].atimeSweepStart
-
-IF ( do_ps EQ 1 ) THEN BEGIN
-   xyouts, 0.5, -0.07, SCATITLE, alignment=0.5, color=255, /normal, $
-           charsize=1., Charthick=1.5
-   xyouts, 0.5, -0.10, TITLE2, alignment=0.5, color=255, /normal, $
-           charsize=1., Charthick=1.5
-ENDIF
-
-IF xxx EQ 'RR' THEN BEGIN
-   xlblstr=STRJOIN([STRING(logbins[0:nlogcats-2],FORMAT='(F0.2)'), $
-                    '>'+STRING(logbins[nlogcats-1],FORMAT='(F5.1)')], ', ', /single)
-   print, '' & print, 'Histogram bin lower bounds (mm/h):'
-   print, xlblstr & print, ''
-   IF ( do_ps EQ 1 ) THEN BEGIN
-      ; write the array of histogram intervals at the bottom of the PDF plot page
-      xyouts, 0.05, -0.15, '!11'+'Histogram bin lower bounds (mm/h):'+'!X', $
-              /NORMAL, COLOR=255, CHARSIZE=0.667, Charthick=1.5
-      xyouts, 0.05, -0.17, '!11'+xlblstr+'!X', /NORMAL, COLOR=255, CHARSIZE=0.667
-   ENDIF
-ENDIF
-
-IF ( do_ps EQ 1 ) THEN BEGIN
-   erase                 ; start a new page in the PS file for the stat tables
-;   device, /landscape   ; change page setup
-   FREE_LUN, tempunit    ; close the temp file for writing
-   OPENR, tempunit2, temptext, /GET_LUN  ; open the temp file for reading
-   statstr = ''
-   fmt='(a0)'
-   xtext = 0.05 & ytext = 0.95
-  ; write the stats tables out to the Postscript file
-   while (eof(tempunit2) ne 1) DO BEGIN
-     readf, tempunit2, statstr, format=fmt
-     xyouts, xtext, ytext, '!11'+statstr+'!X', /NORMAL, COLOR=255, CHARSIZE=0.667
-     ytext = ytext - 0.02
-   endwhile
-  ; close and delete the temp file
-   FREE_LUN, tempunit2
-   FILE_DELETE, temptext, VERBOSE=1
-ENDIF
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-; Build the Scatter Plots
-
-IF pr_or_dpr EQ 'PR' THEN sat_instr = pr_or_dpr $
-ELSE sat_instr = pr_or_dpr+'/'+instrument+'/'+swath
-IF bbparms.meanBB EQ -99.99 THEN skipBB=1 ELSE skipBB=0
-IF xxx EQ 'RR' THEN BEGIN
-   min_xy=0.5
-   max_xy=150.0
-   units='mm/h'
-ENDIF
-
-IF ( do_ps EQ 1 ) THEN BEGIN
-   erase
-   device,/inches,xoffset=0.25,yoffset=0.55,xsize=8.,ysize=10.,/portrait
-   samphgt = (top+botm)/2.0
-   plot_scatter_by_bb_prox_ps, PSFILEpdf, SCATITLE, siteID, yvar, xvar, $
-                            rnType, bbProx, num4hist3, idx4hist3, S2KU=s2ku, $
-                            MIN_XY=min_xy, MAX_XY=max_xy, UNITS=units, $
-                            SAT_INSTR=sat_instr, SKIP_BB=skipBB, HEIGHTS=samphgt
-ENDIF ELSE BEGIN
-   scatwinsize = windowsize > 375  ; constrain size to be 375 pixels or greater
-   plot_scatter_by_bb_prox, SCATITLE, siteID, yvar, xvar, rnType, bbProx, $
-                            num4hist3, idx4hist3, scatwinsize, S2KU=s2ku, $
-                            MIN_XY=min_xy, MAX_XY=max_xy, UNITS=units, $
-                            SAT_INSTR=sat_instr, SKIP_BB=skipBB
-ENDELSE
-
-SET_PLOT, orig_device
-
-; - - - - - - - - - - - - - - - - - - - - - - - -
-
-  ; Build the PPI animation loop.
-
-  ; Check that we have as many sweeps as (startelev+elevs2show); if not, adjust
-  ; elevs2show
-
-   IF (startelev LE mygeometa.num_sweeps ) THEN BEGIN
-      IF (elevs2show+startelev) LE mygeometa.num_sweeps THEN BEGIN
-           nframes = elevs2show
-      ENDIF ELSE BEGIN
-           nframes = mygeometa.num_sweeps - (startelev + 1)
-           print, "Number of sweeps present = ", mygeometa.num_sweeps
-           print, "First, Last sweep requested = ", startelev+1, ',', $
-                  startelev+elevs2show
-           print, "Number of sweeps to show (adjusted): ", nframes
-      ENDELSE
-   ENDIF ELSE BEGIN
-        elevs2show = 1
-        nframes = 1
-        startelev = mygeometa.num_sweeps - 1
-        print, "Number of sweeps present = ", mygeometa.num_sweeps
-        print, "First, Last sweep requested = ", startelev+1, ',', $
-               startelev+elevs2show
-        print, "Showing only sweep number: ", startelev+1
-   ENDELSE
-
-   IF ( elevs2show EQ 0 ) THEN GOTO, nextFile
-   do_pixmap=0
-   IF ( elevs2show GT 1 ) THEN BEGIN
-      do_pixmap=1
-      retain = 0
-      print, ""
-      print, "Please wait while PPI image animation is being prepared..."
-      print, ""
-   ENDIF ELSE retain = 2
-
-   !P.MULTI=[0,1,1]
-   IF ( N_ELEMENTS(windowsize) NE 1 ) THEN windowsize = 375
-   xsize = windowsize[0]
-   ysize = xsize
-   windownum = 2
-   title = ""
-
-   ppi_comn = { winSize : windowsize, $
-                winNum : windownum, $
-                winTitle : TITLE2, $
-                nframes : nframes, $
-                startelev : startelev, $
-                looprate : looprate, $
-                mysweeps : mysweeps, $
-                PPIorient : PPIorient, $
-                PPIbyThresh : PPIbyThresh, $
-                pctString : pctString, $
-                site_Lat : mysite.site_lat, $
-                site_Lon : mysite.site_lon, $
-                site_ID : siteID, $
-                xCorner : xCorner, $
-                yCorner : yCorner, $
-                pr_index : pr_index, $
-                mygeometa : mygeometa, $
-                rntype4ppi : rntype4ppi }
-
-   IF (have_gvrr) THEN gr_rr_zr = ' DP' ELSE gr_rr_zr = ' Z-R'
-   IF pctAbvThresh GT 0.0 AND PPIbythresh THEN sayPct = 1 ELSE sayPct = 0
-
-   IF (haveKdp and haveZdr and haveHID) THEN BEGIN
-      fieldData = ptrarr(3,3, /allocate_heap)
-      fieldIDs = [ ['CZ','CZ','DR'], $
-                   ['RR','RR','FH'], $
-                   ['KD','D0','RH'] ]
-      sources = [ [pr_or_dpr+'/'+instrument, siteID, siteID], $
-                  [pr_or_dpr+'/'+instrument, siteID+gr_rr_zr, siteID], $
-                  [siteID,siteID,siteID] ]
-      thresholded = [ [0,0,0], $
-                      [sayPct,sayPct,0], $
-                      [0,0,0] ]
-      *fieldData[0] = zcor_in
-      *fieldData[1] = gvz_in
-      *fieldData[2] = Zdr
-      *fieldData[3] = rain3_in2
-      *fieldData[4] = gvrr_in2
-      *fieldData[5] = HIDcat
-      *fieldData[6] = Kdp
-      *fieldData[7] = Dzero
-      *fieldData[8] = RHOhv
-   ENDIF ELSE BEGIN
-      fieldData = ptrarr(2,2, /allocate_heap)
-      fieldIDs = [['CZ','CZ'],['RR','RR']]
-      sources = [['PR',siteID],['PR',siteID+gr_rr_zr]]
-      thresholded = [[0,0],[sayPct,sayPct]]
-      *fieldData[0] = zcor_in
-      *fieldData[1] = gvz_in
-      *fieldData[2] = rain3_in2
-      *fieldData[3] = gvrr_in2
-   ENDELSE
-
-   plot_geo_match_ppi_anim_ps, fieldIDs, sources, fieldData, thresholded, $
-                               ppi_comn, DO_PS=do_ps, SHOW_PPIS=show_ppis
-
-   FOR nfieldptr = 0, N_ELEMENTS(fieldData)-1 DO ptr_free, fieldData[nfieldptr]
-
-
-nextFile:
-
-IF ( do_ps EQ 1 ) THEN BEGIN  ; wrap up the postscript file
-   set_plot,/copy,'ps'
-   device,/close
-   SET_PLOT, orig_device
-  ; try to convert it from PS to PDF, using ps2pdf utility
-   command1 = 'which ps2pdf'
-   spawn, command1, result, errout
-   IF result NE '' THEN BEGIN
-      print, 'Converting ', PSFILEpdf, ' to PDF format.'
-      command2 = 'ps2pdf '+PSFILEpdf
-      spawn, command2, result, errout
-      print, 'Removing Postscript version'
-      command3 = 'rm -v '+PSFILEpdf
-      spawn, command3, result, errout
-   ENDIF
-ENDIF
-
-; - - - - - - - - - - - - - - - - - - - - - - -
-
-
-something = ""
-IF nframes LT 2 AND show_ppis THEN BEGIN
-   print, ''
-   READ, something, PROMPT='Hit Return to proceed to next case, Q to Quit: '
-ENDIF
-
-catch, wdel_err
-IF wdel_err EQ 0 THEN BEGIN
-   IF ( elevs2show GT 0 AND nframes GT 0 AND show_ppis ) THEN WDELETE, 2
-ENDIF ELSE BEGIN
-   print, ""
-   print, !error_state.MSG
-   catch, /CANCEL
-   print, "Please do not close non-animating images window manually!  Continue..."
-ENDELSE
-
-errorExit2:
-
-if ( levsdata NE 0 ) THEN BEGIN
-   catch, wdel_err
-   IF wdel_err EQ 0 THEN BEGIN
-      if ( do_ps EQ 0 ) THEN WDELETE, 0
-   ENDIF ELSE BEGIN
-      print, ""
-      print, !error_state.MSG
-      catch, /CANCEL
-      print, "Please do not close profile/PDF window manually!  Continue..."
-   ENDELSE
-   catch, wdel_err
-   IF wdel_err EQ 0 THEN BEGIN
-      if ( do_ps EQ 0 ) THEN WDELETE, 3
-   ENDIF ELSE BEGIN
-      print, ""
-      print, !error_state.MSG
-      catch, /CANCEL
-      print, "Please do not close scatter plot window manually!  Continue..."
-   ENDELSE
-endif
-
-status = 0
-IF something EQ 'Q' OR something EQ 'q' THEN status = 2
-
 errorExit:
-
 return, status
 end
 
@@ -1571,6 +707,43 @@
 ;                 then the user variable will override the NCPATH parameter
 ;                 value on program startup.
 ;
+; subset_method - Method to use to select subset areas from the matchup data:
+;                 'D' = select an area within a cutoff distance (defined by the
+;                       'min_for_subset' parameter) from a user-selected point.
+;                 'V' = select an area of contiguous data values around the
+;                       user-selected start location that are at/above the
+;                       'min_for_subset' value.  The data value to be
+;                       thresholded is defined by the 'rr_or_z' parameter,
+;                       and the threshold applies to the highest data value
+;                       in the vertical column along the PR/DPR ray (e.g., to
+;                       the composite reflectivity).  If either the PR/DPR or
+;                       the matching ground radar value exceeds the threshold,
+;                       then the data for that ray will be included in the
+;                       subset area.
+;                  If subset_method is unspecified then the analysis will be
+;                  performed over the entire domain of the matchup dataset.
+;
+; min_for_subset - Threshold value to be used to define points to be included
+;                  in a user-selected subset area.  If subset_method is 'D',
+;                  then min_for_subset is a distance in km.  If subset_method
+;                  is 'V', then the parameter units are defined by the rr_or_z
+;                  parameter value.  This parameter is ignored if no value is
+;                  specified for subset_method.
+;
+; save_dir       - Optional directory specification to which the subsetted
+;                  variables in a structure will be saved in an IDL SAVE file if
+;                  the user chooses to save them.
+;
+; step_manual    - Flag and Rate value to toggle and control the alternative
+;                  method of animation of PPI images.  If unset, animation is
+;                  automated in an XINTERANIMATE window (default, legacy
+;                  behavior).  If set to a non-zero value, then the PPI images
+;                  will be stepped through under user control: either one at a
+;                  time in forward or reverse, or in an automatic forward
+;                  sequence where the pause, in seconds, between frames is
+;                  defined by the step_manual value.  The automated sequence
+;                  will only play one time in the latter mode, starting from
+;                  the currently-displayed frame.
 
 pro geo_match_3d_rr_or_z_comparisons, RR_OR_Z=rr_or_z, $
                                       INSTRUMENT=instrument, $
@@ -1595,7 +768,11 @@
                                       S2KU = s2ku, $
                                       USE_ZR = use_zr, $
                                       GR_RR_FIELD=gr_rr_field_in, $
-                                      RECALL_NCPATH=recall_ncpath
+                                      RECALL_NCPATH=recall_ncpath, $
+                                      SUBSET_METHOD=subset_method, $
+                                      MIN_FOR_SUBSET=min_for_subset, $
+                                      SAVE_DIR=save_dir, $
+                                      STEP_MANUAL=step_manual
 
 print
 print, "#############################################################"
@@ -1618,6 +795,45 @@
    ENDCASE
 ENDELSE
 
+IF N_ELEMENTS( subset_method ) EQ 1 THEN BEGIN
+   CASE STRUPCASE(subset_method) OF
+       'D' : BEGIN
+               submeth = 'D'
+               IF N_ELEMENTS(min_for_subset) NE 1 THEN BEGIN
+                  print, "Setting a default subset area radius of 20km"
+                  subthresh = 20.
+               ENDIF ELSE subthresh = min_for_subset
+             END
+       'V' : BEGIN
+               submeth = 'V'
+               IF N_ELEMENTS(min_for_subset) NE 1 THEN BEGIN
+                  CASE xxx OF
+                    'Z' : BEGIN
+                          print, "Setting a default subset threshold of 30 dBZ"
+                          subthresh = 30.
+                          END
+                   'RR' : BEGIN
+                          print, "Setting a default subset threshold of 1 mm/h"
+                          subthresh = 1.
+                          END
+                  ENDCASE
+               ENDIF ELSE subthresh = min_for_subset
+             END
+        '' : BREAK   ; silently ignore empty string
+      ELSE : message, "Only allowed values for SUBSET_METHOD are D "+ $
+                      "(distance) and V (value)"
+   ENDCASE
+ENDIF
+
+IF N_ELEMENTS(submeth) EQ 1 AND N_ELEMENTS(save_dir) EQ 1 THEN BEGIN
+  ; check for existence of save_dir, if not empty string
+   IF save_dir NE '' THEN BEGIN
+      IF FILE_TEST(save_dir, /DIRECTORY) THEN real_save_dir = save_dir $
+      ELSE MESSAGE, "SAVE_DIR directory "+save_dir+ $
+                     " does not exist, disabling SAVE files.", /INFO
+   ENDIF
+ENDIF
+
 ; set up the loop speed for xinteranimate, 0<= speed <= 100
 IF ( N_ELEMENTS(looprate) EQ 1 ) THEN BEGIN
   IF ( looprate LT 0 OR looprate GT 100 ) THEN looprate = 3
@@ -1678,8 +894,11 @@
 ENDELSE
 
 IF ( N_ELEMENTS(sitefilter) NE 1 ) THEN BEGIN
-   print, "Defaulting to * for file pattern."
-   ncfilepatt = '*'
+   CASE STRUPCASE(pr_or_dpr) OF
+      'PR' : ncfilepatt = 'GRtoPR.*'
+     'DPR' : ncfilepatt = 'GRtoDPR*'
+   ENDCASE
+   print, "Defaulting to "+ncfilepatt+" for file pattern."
 ENDIF ELSE ncfilepatt = '*'+sitefilter+'*'
 
 PPIorient = keyword_set(ppi_vertical)
@@ -1765,7 +984,7 @@
       print, "Setting GR_RR_FIELD value to RR (for DROPS)."
       gr_rr_field = 'RR'
    ENDELSE
-ENDIF
+ENDIF ELSE gr_rr_field = ''
 
 ; specify whether to skip graphical PPI output to screen in Postscript mode
 IF ( PS_DIR NE '' AND KEYWORD_SET(batch) ) THEN do_batch=1 ELSE do_batch=0
@@ -1808,7 +1027,9 @@
                                       hideRntype, pr_or_dpr, ALT_BB_HGT=alt_bb_hgt, $
                                       PS_DIR=ps_dir, B_W=b_w, S2KU=s2ku, ZR=zr_force, $
                                       GR_RR_FIELD=gr_rr_field, BATCH=do_batch, $
-                                      MAX_RANGE=max_range_in )
+                                      MAX_RANGE=max_range_in, SUBSET_METHOD=submeth, $
+                                      MIN_FOR_SUBSET=subthresh, SAVE_DIR=real_save_dir, $
+                                      STEP_MANUAL=step_manual )
 
          if (action EQ 2) then break    ; manual request to quit
       endfor
@@ -1823,15 +1044,18 @@
                                  hideRntype, pr_or_dpr, ALT_BB_HGT=alt_bb_hgt, $
                                  PS_DIR=ps_dir, B_W=b_w, S2KU=s2ku, ZR=zr_force, $
                                  GR_RR_FIELD=gr_rr_field, BATCH=do_batch, $
-                                 MAX_RANGE=max_range_in )
+                                 MAX_RANGE=max_range_in, SUBSET_METHOD=submeth, $
+                                 MIN_FOR_SUBSET=subthresh, SAVE_DIR=real_save_dir, $
+                                 STEP_MANUAL=step_manual )
 
-      if (action NE 0) then break    ; either an error, or manual request to quit
+;      if (action NE 0) then break    ; either an error, or manual request to quit
       newpathpr = FILE_DIRNAME(ncfilepr)  ; set the path to the last file's path
       IF KEYWORD_SET(recall_ncpath) THEN BEGIN
          ; define/assign new default path for session as user system variable
           IF (haveUserVar EQ 1) THEN !LAST_NCPATH = newpathpr $
           ELSE DEFSYSV, '!LAST_NCPATH', newpathpr
       ENDIF
+      PRINT, "Select next file to process or close File Selector to quit."
       ncfilepr = dialog_pickfile(path=newpathpr, filter = ncfilepatt)
    endwhile
 ENDELSE
