diff --git a/CHANGELOG.md b/CHANGELOG.md
index a4db64af6..f491a6422 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,7 +4,7 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 ## 1.9.6
-November 15, 2023
+November 18, 2023
 
 Added highrate (1-sec) GNSS archive from Spain, IGN ES.  RINEX 3 only. 
 For this archive, use ignes and make sure to say -rate high -samplerate 1 . 
@@ -15,6 +15,16 @@ output file where frequency is being written in column 11.
 
 Updated sc02 usecase and notebook install instructions.
 
+Provided more feedback to people using nolook option in rinex2snr. Changed "local directory"
+printed to the screen to the actual directory.
+
+Substantially changed plots in subdaily to use datetime instead of day of year.  Ultimately 
+all of subdaily should be changed to MJD so as to allow datasets crossing year boundaries.
+It also removes gaps from the last spline output plot, but it has not been implemented in the output
+txt file as yet (nor in the second to last plot)
+
+
+
 ## 1.9.5
 November 13, 2023
 
diff --git a/README.md b/README.md
index 83bc43f17..1a97f6c3f 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # gnssrefl
 
-**github version: 1.9.5** [![PyPI Version](https://img.shields.io/pypi/v/gnssrefl.svg)](https://pypi.python.org/pypi/gnssrefl) [![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.5601495.svg)](http://dx.doi.org/10.5281/zenodo.5601495) [![Documentation Status](https://readthedocs.org/projects/gnssrefl/badge/?version=latest)](https://gnssrefl.readthedocs.io/en/latest/?badge=latest)
+**github version: 1.9.6** [![PyPI Version](https://img.shields.io/pypi/v/gnssrefl.svg)](https://pypi.python.org/pypi/gnssrefl) [![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.5601495.svg)](http://dx.doi.org/10.5281/zenodo.5601495) [![Documentation Status](https://readthedocs.org/projects/gnssrefl/badge/?version=latest)](https://gnssrefl.readthedocs.io/en/latest/?badge=latest)
 
 **For those of you trying to convert RINEX 3 files to RINEX 2.11, be careful. There are a lot of things going
 on in that translation.  It is better to let gnssrefl to do that conversion for you.  rinex2snr allows RINEX 3
diff --git a/gnssrefl/rinex2snr.py b/gnssrefl/rinex2snr.py
index 3f2438802..7f0dc6cd7 100755
--- a/gnssrefl/rinex2snr.py
+++ b/gnssrefl/rinex2snr.py
@@ -206,12 +206,16 @@ def run_rinex2snr(station, year_list, doy_list, isnr, orbtype, rate,dec_rate,arc
                 rgz = station + cdoy + '0.' + cyy + 'o.gz'
                 localpath2 =  os.environ['REFL_CODE'] + '/' + cyyyy + '/rinex/' + station + '/'
                 if nol:
-                    print('Will first assume RINEX file ', station, ' year:', year, ' doy:', doy, 'is in the local directory')
+                    print('Will first assume RINEX file ', station, ' year:', year, ' doy:', doy, 'is in this directory:')
+                    current_local = os.getcwd()
+                    print(current_local)
                     # this assumes RINEX file is in local directory or "nearby"
                     if version == 2:
+
                         if mk:
                             the_makan_option(station,cyyyy,cyy,cdoy) # looks everywhere in your local directories
                         if not os.path.exists(r):
+                            print('Unsuccessful, so how trying other directories')
                             # could try this way? - look for file in localpath2. gunzip if necessary
                             allgood = get_local_rinexfile(r,localpath2)
                         if os.path.exists(r):
@@ -1341,7 +1345,8 @@ def get_local_rinexfile(rfile,localpath2):
                 allgood = True
                 subprocess.call(['rm',rd])
 
-    # now check in $REFL_CODE/YYYY/rinex/ssss
+    # now check in 
+    print('Checking for the file in : ', localpath2)
     if not allgood:
         r = localpath2 + rfile
         # hatanaka version in REFL_CODE
diff --git a/gnssrefl/sd_libs.py b/gnssrefl/sd_libs.py
index 222a10cf3..8cb9d5aba 100644
--- a/gnssrefl/sd_libs.py
+++ b/gnssrefl/sd_libs.py
@@ -3,10 +3,32 @@
 import math
 import matplotlib.pyplot as plt
 import numpy as np
+from astropy.time import Time
 
+# my code
 import gnssrefl.gps as g
 import gnssrefl.gnssir_v2 as guts2
 
+def mjd_to_obstimes(mjd):
+    """
+    takes mjd array and converts to datetime for plotting.
+
+    Parameters
+    ----------
+    mjd : numpy array of floats
+        mod julian date
+
+    Returns
+    -------
+    dt : numpy array of datetime objects
+
+    """
+
+    dt = Time(mjd,format='mjd').utc.datetime; 
+
+    return dt
+
+
 def write_spline_output(splineout, iyear, th, spline, delta_out, station, txtdir,Hortho):
     """
     Writing the output of the spline fit to the final RH time series.
@@ -394,7 +416,7 @@ def stack_two_more(otimes,tv,ii,jj,stats, station, txtdir, sigma,kplt,hires_figs
     fs = 10
     otimesarray = np.asarray(otimes)
     # new plot 
-    plt.figure()
+    plt.figure(figsize=(10,5))
     plt.plot(tv[:,5],tv[:,2], '.',markersize=4,label='obs')
     plt.plot(tv[ii,5],tv[ii,2], 'ro',markersize=4,label='outliers')
     plt.xlabel('Azimuth (degrees)')
@@ -624,9 +646,9 @@ def writejsonfile(ntv,station, outfile):
 
     return True
 
-def rhdot_plots(th,correction,rhdot_at_th, tvel,yvel,fs,station,txtdir,hires_figs):
+def rhdot_plots(th,correction,rhdot_at_th, tvel,yvel,fs,station,txtdir,hires_figs,year):
     """
-    make rhdot correction plots
+    makes the rhdot correction plots
 
     Parameters
     ----------
@@ -648,8 +670,15 @@ def rhdot_plots(th,correction,rhdot_at_th, tvel,yvel,fs,station,txtdir,hires_fig
         file directory for output
     hires_figs : bool
         whether you want eps instead of png 
+    year : int
+        calendar year
 
     """
+    mjd0 = g.fdoy2mjd(year, th[0] ) # MJD of first point
+    # not being used
+    #th_obs = mjd_to_obstimes(th+mjd0)
+
+
     fig=plt.figure(figsize=(10,6))
     plt.subplot(2,1,1)
     plt.plot(th, correction,'b.')
@@ -657,32 +686,45 @@ def rhdot_plots(th,correction,rhdot_at_th, tvel,yvel,fs,station,txtdir,hires_fig
     plt.xlim((np.min(th), np.max(th)))
     plt.grid()
     plt.title('The RH correction from RHdot effect ')
-    #plt.xlim((plot_begin, plot_end))
+    #plt.xlim((th_obs[0], th_obs[-1])) #fig.autofmt_xdate()
 
     plt.subplot(2,1,2)
-    A1 = np.min(th) ; A2 = np.max(th)
+    #mjd1 = g.fdoy2mjd(year, tvel[0] ) # MJD of first point
+    #tvel_obs = mjd_to_obstimes(tvel+mjd1)
+
+    # ???
+    A1 = np.min(th)  ; A2 = np.max(th) 
     jj = (tvel >= A1) & (tvel <= A2)
-    plt.plot(th, rhdot_at_th,'bo',label='at GNSS obs')
+
+    plt.plot(th, rhdot_at_th,'b.',label='at GNSS obs')
+    # disaster - did not work
+    #plt.plot(tvel_obs, yvel,'c-', label='spline fit')
     plt.plot(tvel[jj], yvel[jj],'c-', label='spline fit')
     plt.legend(loc="upper left")
     plt.grid()
     plt.title('surface velocity')
-    plt.ylabel('meters/hour'); plt.xlabel('days of the year')
+    plt.ylabel('meters/hour'); plt.xlabel('days of the year') 
     plt.xlim((A1,A2))
+    #plt.xlim((th[0], th[-1])) #fig.autofmt_xdate()
+
     if hires_figs:
         g.save_plot(txtdir + '/' + station + '_rhdot3.eps')
     else:
         g.save_plot(txtdir + '/' + station + '_rhdot3.png')
 
 
-def RH_ortho_plot( station, H0, th_even, spline_whole_time,txtdir, fs):
+def RH_ortho_plot( station, H0, year, th_even, spline_whole_time,txtdir, fs, time_rh, rh, gap_min_val):
     """
+    Makes a plot of the final spline fit to the data. Output time interval controlled by the user.
+
     Parameters
     ----------
     station : str
         name of station, 4 ch
     H0 : float
         datum correction (orthometric height) to convert RH to MSL data, in meters
+    year : int
+        year of the time series (ultimately should not be needed)
     th_even : numpy of floats
         time in fractional days of of year, I think
     spline_whole_time : numpy of floats
@@ -691,18 +733,55 @@ def RH_ortho_plot( station, H0, th_even, spline_whole_time,txtdir, fs):
        location of plot
     fs : int
         fontsize
+    time_rh : numpy of floats
+        time of rh values, in fractional doy I believe
+    rh : numpy of floats
+        refl hgt in meters
+    gap_min_val : float
+        minimum length gap allowed, in day of year units
     """
 
+
+    # difference function to find time between all rh measurements
+    gdiff = np.diff( time_rh )
+    mjdx = g.fdoy2mjd(year, time_rh[0] ) # 
+    mjdt = mjdx + (time_rh - time_rh[0])
+    mjd = mjdt[0:-1]
+
+    # MJD of the gdiff values
+    # mjd = time_rh[0:-1] + mjdx
+    #print('base value ', year, time_rh[0], mjdx, mjd[-1])
+
+    # find indices of gaps  that are larger than a certain value
+    ii = (gdiff > gap_min_val)
+    N = len(mjd[ii])
+    Ngdiff = len(gdiff)
+
+    # figure out mjd of the first point
+    mjd0 = g.fdoy2mjd(year, th_even[0])
+    mjd_even = mjd0 + (th_even - th_even[0])  
+    mjd_even_obstimes = mjd_to_obstimes(mjd_even)
+
+
+    if (Ngdiff > 0):
+        for i in range(0,Ngdiff):
+            if ii[i]:
+                e0 = mjd[i] ; e1 = mjd[i+1] # beginning and ending of the gap
+                #newl = [ e0, e1 ]
+                #print('gap', mjd[i],mjd[i+1])
+                bad_indices = (mjd_even > e0) & (mjd_even < e1 )
+                #mjd_even[bad_indices] = np.nan
+                mjd_even_obstimes[bad_indices] = np.datetime64("NaT")
+
+
     fig=plt.figure(figsize=(10,5))
-    #H0 = find_ortho_height(station,'')
-    plt.plot(th_even, H0 - spline_whole_time, 'c-')
+    plt.plot(mjd_even_obstimes, H0 -spline_whole_time, 'c')
+    #plt.plot(mjd_even, H0 -spline_whole_time, 'c')
     plt.grid()
-    plt.xlabel('Time (days)',fontsize=fs)
-    # would like to add year, but it is not currently sent to this code, alas
-    # 
-    #plt.xlabel('Time (days)/' + str(year) ,fontsize=fs)
     plt.ylabel('meters',fontsize=fs)
     plt.title(station.upper() + ' Water Level ', fontsize=fs)
+    fig.autofmt_xdate()
+
     pfile = txtdir + '/' + station + '_H0.png'
     g.save_plot(pfile)
 
@@ -812,3 +891,93 @@ def quickTr(year, doy,frachours):
 
 
     return datestring
+
+
+def subdaily_resids_last_stage(station, year, th, biasCor_rh, spline_at_GPS, fs, strsig, hires_figs, 
+                               txtdir, ii, jj, th_even,spline_whole_time):
+    """
+    Makes the final residual plot for subdaily (after RHdot and IF correction made).
+    Returns the bad points ...
+
+
+    Parameters
+    ----------
+    station : str
+        4 character station name
+    year : int
+        calendar year
+    th : numpy array of ??
+        time variable of some kind, fractional day of year ?
+    biasCor_rh : numpy array of floats
+        refl hgts that have been corrected for RHdot and IF
+    spline_at_GPS : numpy array of floats
+        RH derived From the spline fit and calculated at GPS time tags
+    fs : int
+        font size
+    strsig : str
+        sigma string to go on the legend
+    hires_figs : bool
+        whether to save the plots with better resolution
+    txtdir : str
+        directory where the plot will be saved
+    ii : numpy array
+        indices of the outliers?
+    jj : numpy array
+        indices of the values to keep?
+    th_even : numpy array
+        evenly spaced time values, day of year
+    spline_whole_time : numpy array of flots
+        splinefit for ???
+
+    Returns
+    -------
+    badpoints2 : numpy array of floats
+         RH residuals
+    """
+    # convert to mjd and then obstimes
+    mjd0 = g.fdoy2mjd(year, th[0])
+    th_obs = mjd_to_obstimes(mjd0 + th - th[0])
+    mjd1 = g.fdoy2mjd(year, th_even[0])
+    th_even_obs = mjd_to_obstimes(mjd1 + th_even - th_even[0])
+
+
+    fig=plt.figure(figsize=(10,5))
+    plt.plot(th_obs, biasCor_rh, 'b.', label='RH ' + strsig)
+    plt.plot(th_even_obs, spline_whole_time, 'c-',label='spline')
+    plt.plot(th_obs[ii], biasCor_rh[ii], 'rx', label='outliers')
+
+    plt.legend(loc="best",fontsize=fs)
+    plt.grid() ; plt.gca().invert_yaxis()
+    plt.ylabel('meters',fontsize=fs); 
+    #plt.xlabel('days of the year',fontsize=fs)
+    plt.title(station.upper() + ' RH with RHdot and InterFrequency Corrections Applied',fontsize=fs)
+    fig.autofmt_xdate()
+
+    # put hires_figs boolean here
+    if hires_figs:
+        g.save_plot(txtdir + '/' + station + '_rhdot4.eps')
+    else:
+        g.save_plot(txtdir + '/' + station + '_rhdot4.png')
+
+    fig=plt.figure(figsize=(10,5))
+    plt.plot(th_obs, biasCor_rh - spline_at_GPS, 'b.',label='all residuals')
+    plt.title('Station:' + station + ' Residuals to new spline fit',fontsize=fs)
+    plt.grid()
+    plt.ylabel('meters',fontsize=fs)
+    #plt.xlabel('days of the year',fontsize=fs)
+
+    plt.plot(th_obs[ii], (biasCor_rh -spline_at_GPS)[ii], 'r.',label='outliers')
+    # will write these residauls out to a file
+    badpoints2 =  (biasCor_rh -spline_at_GPS)[ii]
+    plt.legend(loc="upper left",fontsize=fs)
+    fig.autofmt_xdate()
+
+    # print('\n RMS with frequency biases and RHdot taken out (m) ', np.round(newsigma,3) , '\n' )
+    if hires_figs:
+        g.save_plot(txtdir + '/' + station + '_rhdot5.eps')
+    else:
+        g.save_plot(txtdir + '/' + station + '_rhdot5.png')
+
+    plt.close() # dont send this figure to the screen
+
+    return badpoints2
diff --git a/gnssrefl/subdaily.py b/gnssrefl/subdaily.py
index 88a5a5a68..139a94f68 100644
--- a/gnssrefl/subdaily.py
+++ b/gnssrefl/subdaily.py
@@ -693,16 +693,14 @@ def flipit(tvd,col):
     tnew = tnew[ii]
     ynew = ynew[ii]
 
-    # this is what should be used for the spline.
-    #plt.figure()
-    #plt.plot(tnew,ynew,'.')
-    #plt.show()
     return tnew, ynew
 
 
 def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs):
     """
-    Improved code to compute rhdot correction and bias correction for subdaily
+    Improved code to compute rhdot correction and interfrequency bias correction for subdaily
+    this assumes you have at least removed crude outliers in the previous section of the 
+    subdaily code.
 
     Parameters
     ----------
@@ -730,17 +728,26 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     apply_rhdot : bool, optional
         whether you want to apply the rhdot correction
         default is true
+    gap_min_val : float, optional
+        gap allowed in last spline, in hours
+    year : int
+        hopefully this will go away ...
+
 
     """
     # output will go to REFL_CODE/Files unless txtdir provided
     xdir = os.environ['REFL_CODE']
+    gap_min_val = kwargs.get('gap_min_val',6.0)
+    gap_min_val = gap_min_val/24 # change to DOY units
 
     val = kwargs.get('txtdir',[])
+
     if len(val) == 0:
         txtdir = xdir + '/Files/'
     else:
         txtdir = val
 
+    year = kwargs.get('year',0)
 
     delta_out = kwargs.get('delta_out',0)
     if (delta_out  == 0):
@@ -791,6 +798,7 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     print('\nComputes rhdot correction and interfrequency bias correction for subdaily')
     print('\nInput filename:', fname)
     print('\nOutput filename: ', fname_new)
+    print('\nMinimum gap allowed in spline output, in day units', gap_min_val)
 
     # read in the lomb scargle values which are the output of gnssir
     # i.e. the reflector heights
@@ -829,8 +837,8 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     print('Average num of RH obs per day', '{0:5.1f} '.format (len(h)/Ndays) )
     print('Knots per day: ', knots_per_day, ' Number of knots: ', numKnots)
     # currently using 3 sigma
-    print('Outlier criterion for the first splinefit (m):', outlierV)
-    print('Outlier criterion for the second splinefit (m):', outlierV2)
+    print('Outlier criterion provided by user for the first splinefit (m):', outlierV)
+    print('Outlier criterion provided by user for the second splinefit (m):', outlierV2)
 
     firstKnot_in_minutes = 15
     t1 = tnew.min()+firstKnot_in_minutes/60/24
@@ -855,9 +863,6 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
 
     obsPerHour= perday/24
 
-    fig=plt.figure(figsize=(10,6))
-
-    plt.subplot(2,1,1)
 
     tvel = spl_x[1:N]
     yvel = obsPerHour*np.diff(spl_y)
@@ -883,11 +888,21 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     label1 = 'w/o RHdot ' + str(np.round(sigmaBefore,3)) + 'm'
     label2 = 'w RHdot ' + str(np.round(sigmaAfter,3)) + 'm'
 
-    plt.plot(th, h, 'b.', label=label1,markersize=4)
+    # start the figure, convert from doy to MJD to obstime
+    mjd0 = g.fdoy2mjd(year,th[0])
+    th_obs = sd.mjd_to_obstimes(mjd0 + th-th[0])
+
+    fig=plt.figure(figsize=(10,6))
+    plt.subplot(2,1,1)
+    plt.plot(th_obs, h, 'b.', label=label1,markersize=4)
     iw = (spl_x > th[0]) & (spl_x < th[-1])
-    plt.plot(spl_x[iw], spl_y[iw], 'c--', label='spline') # otherwise wonky spline makes a goofy plot
 
-    plt.plot(th,correctedRH,'m.',label=label2,markersize=4)
+    # change to datetime
+    tm_mjd = g.fdoy2mjd(year, spl_x[iw][0])
+    spl_x_obs = sd.mjd_to_obstimes(tm_mjd + spl_x[iw] - spl_x[iw][0] )
+    plt.plot(spl_x_obs, spl_y[iw], 'c--', label='spline') # otherwise wonky spline makes a goofy plot
+
+    plt.plot(th_obs,correctedRH,'m.',label=label2,markersize=4)
 
     print('\nRMS no RHdot correction (m)', '{0:6.3f}'.format ( sigmaBefore))
     print('RMS w/ RHdot correction (m)', '{0:6.3f} \n'.format ( sigmaAfter ))
@@ -895,15 +910,16 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     plt.gca().invert_yaxis()
     plt.legend(loc="upper left")
     plt.ylabel('meters',fontsize=fs)
-    plt.title( station.upper() + ' Reflector Heights')
+    plt.title( station.upper() + ' Reflector Heights',fontsize=fs)
     outlierstring = str(outlierV) + '(m) outliers'
     plt.yticks(fontsize=fs); plt.xticks(fontsize=fs)
     plt.grid()
-    plt.xlim((plot_begin, plot_end))
+    fig.autofmt_xdate()
+    #plt.xlim((plot_begin, plot_end))
 
 
     plt.subplot(2,1,2)
-    plt.plot(th, residual_after,'r.',label='all pts')
+    plt.plot(th_obs, residual_after,'r.',label='all pts')
 
     if outlierV is None:
     # use 3 sigma
@@ -931,13 +947,15 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     correctedRH_new = correctedRH[ii]
     NV = len(correctedRH)
 
-    plt.plot(th[ii], residual_after[ii],'b.',label='kept pts')
+    plt.plot(th_obs[ii], residual_after[ii],'b.',label='kept pts')
     plt.grid()
-    plt.title('Residuals to the spline fit',fontsize=fs)
-    plt.xlabel('days of the year',fontsize=fs)
+    #plt.title('Residuals to the spline fit',fontsize=fs)
+    #plt.xlabel('days of the year',fontsize=fs)
     plt.ylabel('meters',fontsize=fs)
     plt.legend(loc="upper left",fontsize=fs)
-    plt.xlim((plot_begin, plot_end))
+    plt.yticks(fontsize=fs); plt.xticks(fontsize=fs)
+    #plt.xlim((plot_begin, plot_end))
+    fig.autofmt_xdate()
 
     if hires_figs:
         g.save_plot(txtdir + '/' + station + '_rhdot2.eps')
@@ -949,7 +967,7 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     residual_after = residual_after[ii]
 
     # make the RHdot plot as well
-    sd.rhdot_plots(th,correction,rhdot_at_th, tvel,yvel,fs,station,txtdir,hires_figs)
+    sd.rhdot_plots(th,correction,rhdot_at_th, tvel,yvel,fs,station,txtdir,hires_figs,year)
 
     writecsv = False ; extraline = ''
     # write out the new solutions with RHdot and without 3 sigma outliers
@@ -962,6 +980,7 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     # this is horrible code and needs to be fixed
     biasCorrected_RH = correctedRH_new
 
+    print('----------------------------------------------------------\n')
     print('Check inter-frequency biases for GPS, Glonass, and Galileo\n')
 
     if 1 not in tvd_new[:,10]:
@@ -1003,6 +1022,7 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
         #if pltit:
         #    plt.show()
         return tvd_new, correction
+    print('----------------------------------------------------------\n')
 
     # now try to write the bias corrected values
     # first attempt was wrong because i forgot to sort the corrected column in biasCorrected_RH
@@ -1044,9 +1064,6 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     newsigma = np.std(biasCor_rh-spline_at_GPS)
     strsig = str(round(newsigma,3)) + '(m)'
 
-    fig=plt.figure(figsize=(10,5))
-    plt.plot(th, biasCor_rh, 'b.', label='RH ' + strsig)
-    plt.plot(th_even, spline_whole_time, 'c-',label='spline')
 
     if outlierV2 is None:
         # use 3 sigma to find outliers
@@ -1062,55 +1079,23 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
         # points to keep
         jj = np.abs(biasCor_rh -spline_at_GPS) <= OutlierLimit
 
-
-    plt.plot(th[ii], biasCor_rh[ii], 'rx', label='outliers')
-
-    plt.legend(loc="best",fontsize=fs)
-    plt.grid()
-    plt.gca().invert_yaxis()
-    plt.ylabel('meters',fontsize=fs)
-    plt.title(station.upper() + ' RH with RHdot and InterFrequency Corrections Applied',fontsize=fs)
-    plt.xlabel('days of the year',fontsize=fs)
-    # put hires_figs boolean here
-    if hires_figs:
-        g.save_plot(txtdir + '/' + station + '_rhdot4.eps')
-    else:
-        g.save_plot(txtdir + '/' + station + '_rhdot4.png')
-
-
-    H0 = sd.find_ortho_height(station,'')
-    sd.RH_ortho_plot( station, H0, th_even, spline_whole_time,txtdir, fs)
     
+    # make the plot externally now
+    badpoints2 = sd.subdaily_resids_last_stage(station, year, th, biasCor_rh, spline_at_GPS, 
+                                               fs, strsig, hires_figs,txtdir, ii,jj,th_even, spline_whole_time)
+    H0 = sd.find_ortho_height(station,'')
+    sd.RH_ortho_plot( station, H0, year, th_even, spline_whole_time,txtdir, fs, th[jj],biasCor_rh[jj],gap_min_val)
+    print('\nRMS with frequency biases and RHdot taken out (m) ', np.round(newsigma,3) , '\n' )
 
-    fig=plt.figure(figsize=(10,5))
-    #plt.subplot(2,1,1)
-    plt.plot(th, biasCor_rh - spline_at_GPS, 'b.',label='all residuals')
-    plt.title('Station:' + station + ' Residuals to new spline fit',fontsize=fs)
-    plt.grid()
-    plt.ylabel('meters',fontsize=fs)
-    plt.xlabel('days of the year',fontsize=fs)
-
-    plt.plot(th[ii], (biasCor_rh -spline_at_GPS)[ii], 'r.',label='outliers')
-    # will write these residauls out to a file
-    badpoints2 =  (biasCor_rh -spline_at_GPS)[ii]
-    plt.legend(loc="upper left",fontsize=fs)
-
-    print('\n RMS with frequency biases and RHdot taken out (m) ', np.round(newsigma,3) , '\n' )
-    if hires_figs:
-        g.save_plot(txtdir + '/' + station + '_rhdot5.eps')
-    else:
-        g.save_plot(txtdir + '/' + station + '_rhdot5.png')
-
-    plt.close() # dont send this figure to the screen
 
-    # bias corrected - and again without 3 sigma outliers
+    # write out the files with RH dot and IF bias corrected - and again without 3 sigma outliers
     bias_corrected_filename = fname_new + 'IF'; extraline = ''; writecsv = False
     biasCor_rh = tvd_new[jj,24]
     write_subdaily(bias_corrected_filename,station,tvd_new[jj,:], writecsv,extraline, newRH_IF=biasCor_rh)
 
     new_outliers = tvd_new[ii,:]
 
-    # write outliers again ... 
+    # write outliers to a file ... again ... 
     sd.writeout_spline_outliers(new_outliers,txtdir,badpoints2,'outliers.spline2.txt')
 
     # I was looking at issue of delT being too big
@@ -1118,13 +1103,6 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     year = int(tvd[0,0]);
     sd.write_spline_output(splineout, year, th, spline, delta_out,station,txtdir,H0)
 
-    if  False:
-        plt.figure()
-        res = (h-spline_at_GPS)
-        plt.plot( tvd[ii,14],  res[ii], 'o', label='outlier')
-        plt.plot( tvd[:,14], res, '.',label='residuals')
-        plt.xlabel('delta Time (minutes)')
-
 
     return tvd, correction
 
diff --git a/gnssrefl/subdaily_cl.py b/gnssrefl/subdaily_cl.py
index c1cbd55cf..3e73f479d 100644
--- a/gnssrefl/subdaily_cl.py
+++ b/gnssrefl/subdaily_cl.py
@@ -42,8 +42,9 @@ def parse_arguments():
     parser.add_argument("-knots_test", default=None, type=int, help="test knots")
     parser.add_argument("-hires_figs", default=None, type=str, help="hi-resolution eps figures, default is False")
     parser.add_argument("-apply_rhdot", default=None, type=str, help="apply rhdot, default is True")
-    parser.add_argument("-fs", default=None, type=int, help="fontsize for figures. default is 12")
+    parser.add_argument("-fs", default=None, type=int, help="fontsize for figures. default is 10")
     parser.add_argument("-alt_sigma", default=None, type=str, help="boolean test for alternate sigma definition. default is False")
+    parser.add_argument("-gap_min_val", default=None, type=float, help="min gap allowed in splinefit output file. default is 6 hours")
 
     args = parser.parse_args().__dict__
 
@@ -61,7 +62,7 @@ def subdaily(station: str, year: int, txtfile_part1: str = '', txtfile_part2: st
         doy2: int = 366, testing: bool = True, ampl: float = 0, h1: float=0.4, h2: float=300.0, 
         azim1: int=0, azim2: int = 360, peak2noise: float = 0, kplt: bool = False, 
         subdir: str = None, delta_out : int = 1800, if_corr: bool = True, knots_test: int = 0, 
-             hires_figs : bool=False, apply_rhdot : bool=True, fs: int = 12, alt_sigma: bool= False):
+             hires_figs : bool=False, apply_rhdot : bool=True, fs: int = 10, alt_sigma: bool= False, gap_min_val: float=6.0):
     """
     Subdaily combines multiple day gnssir solutions and applies relevant corrections. 
     It only works for one year at a time; you can restricts time periods within a year with -doy1 and -doy2
@@ -175,10 +176,14 @@ def subdaily(station: str, year: int, txtfile_part1: str = '', txtfile_part2: st
         whether you want the RH dot correction applied
         for a lake or river you would not want it to be.
     fs : int, optional
-        fontsize for Figures. default is 12 for now.
+        fontsize for Figures. default is 10 for now.
     alt_sigma : bool, optional
         whether you want to use Nievinski definition for outlier criterion.
-        will change to kwargs when I get a chance.
+        in part 1 of the code (the crude outlier detector)
+
+    gap_min_val : float, optional
+        removes splinefit values from output txt and plot for gaps 
+        bigger than this value, in hours
 
     """
 
@@ -248,7 +253,8 @@ def subdaily(station: str, year: int, txtfile_part1: str = '', txtfile_part2: st
     if rhdot:
        tv, corr = t.rhdot_correction2(station, input2spline, output4spline, plt, spline_outlier1, spline_outlier2, 
                    knots=knots,txtdir=txtdir,testing=testing,delta_out=delta_out,
-                   if_corr=if_corr,knots_test=knots_test,hires_figs=hires_figs,apply_rhdot=apply_rhdot,fs=fs)
+                   if_corr=if_corr,knots_test=knots_test,hires_figs=hires_figs,
+                   apply_rhdot=apply_rhdot,fs=fs,gap_min_val=gap_min_val,year=year)
        if plt:
            mplt.show()
 
diff --git a/setup.py b/setup.py
index a55f4cf3b..71ecb652e 100644
--- a/setup.py
+++ b/setup.py
@@ -42,7 +42,7 @@
 ]
 setup(
     name="gnssrefl",
-    version="1.9.5",
+    version="1.9.6",
     author="Kristine Larson",
     author_email="kristinem.larson@gmail.com",
     description="A GNSS reflectometry software package ",