diff --git a/gnssrefl/sd_libs.py b/gnssrefl/sd_libs.py
index e99098641..9d5cd6324 100644
--- a/gnssrefl/sd_libs.py
+++ b/gnssrefl/sd_libs.py
@@ -939,7 +939,7 @@ def subdaily_resids_last_stage(station, year, th, biasCor_rh, spline_at_GPS, fs,
 
     return badpoints2
 
-def RH_ortho_plot2( station, H0, year,  txtdir, fs, time_rh, rh, gap_min_val,th,spline,delta_out):
+def RH_ortho_plot2( station, H0, year,  txtdir, fs, time_rh, rh, gap_min_val,th,spline,delta_out,csvfile):
     """
 
     Makes a plot of the final spline fit to the data. Output time interval controlled by the user.
@@ -970,6 +970,8 @@ def RH_ortho_plot2( station, H0, year,  txtdir, fs, time_rh, rh, gap_min_val,th,
         used for fitting 
     delta_out : int
         how often spline is printed, in seconds
+    csvfile : bool
+        print out csv instead of plain txt
     """
 
     firstpoint = float(th[0]); lastpoint = float(th[-1])
@@ -1006,8 +1008,12 @@ def RH_ortho_plot2( station, H0, year,  txtdir, fs, time_rh, rh, gap_min_val,th,
     fout.write('{0:1s}  {1:30s}  {2:8.3f} \n'.format('%','orthometric height minus RH, where Hortho (m) is ', H0  ))
     #fout.write('{0:1s}  {1:30s}  \n'.format('%','This assumes RH is measured relative to the L1 phase center.  '))
     #fout.write('{0:1s}  {1:30s}  \n'.format('%','MJD, RH(m), YY,MM,DD,HH,MM,SS, quasi-sea-level(m)'))
-    fout.write('{0:1s}  {1:30s}  \n'.format('%','MJD              RH(m)  YYYY  MM  DD  HH  MM  SS   quasi-sea-level(m)'))
-    fout.write('{0:1s}  {1:30s}  \n'.format('%','(1)               (2)   (3)  (4) (5)  (6) (7) (8)    (9)'))
+    if csvfile:
+        fout.write('{0:1s}  {1:30s}  \n'.format('%','MJD,             RH(m), YYYY, MM, DD, HH, MM, SS,  quasi-sea-level(m)'))
+        fout.write('{0:1s}  {1:30s}  \n'.format('%','(1),              (2),  (3), (4),(5), (6),(7),(8),   (9)'))
+    else:
+        fout.write('{0:1s}  {1:30s}  \n'.format('%','MJD              RH(m)  YYYY  MM  DD  HH  MM  SS   quasi-sea-level(m)'))
+        fout.write('{0:1s}  {1:30s}  \n'.format('%','(1)               (2)   (3)  (4) (5)  (6) (7) (8)    (9)'))
 
 
     # difference function to find time between all rh measurements
@@ -1050,8 +1056,10 @@ def RH_ortho_plot2( station, H0, year,  txtdir, fs, time_rh, rh, gap_min_val,th,
         for i in range(0,N_new):
             if not np.isnan(spline_new[i]):
                 rhout = spline_new[i]
-                fout.write('{0:15.7f}  {1:10.3f} {2:4.0f} {3:3.0f} {4:3.0f} {5:3.0f} {6:3.0f} {7:3.0f} {8:10.3f} \n'.format(
-                    mjd_new[i], rhout, theyear[i], xm[i], xd[i], xh[i], xmin[i], xs[i], H0-rhout))
+                if csvfile:
+                    fout.write('{0:15.7f}, {1:10.3f},{2:4.0f},{3:3.0f},{4:3.0f},{5:3.0f},{6:3.0f},{7:3.0f},{8:10.3f} \n'.format( mjd_new[i], rhout, theyear[i], xm[i], xd[i], xh[i], xmin[i], xs[i], H0-rhout))
+                else:
+                    fout.write('{0:15.7f}  {1:10.3f} {2:4.0f} {3:3.0f} {4:3.0f} {5:3.0f} {6:3.0f} {7:3.0f} {8:10.3f} \n'.format(mjd_new[i], rhout, theyear[i], xm[i], xd[i], xh[i], xmin[i], xs[i], H0-rhout))
 
 
     fout.close()
diff --git a/gnssrefl/subdaily.py b/gnssrefl/subdaily.py
index 5ce447aff..7e30cb206 100644
--- a/gnssrefl/subdaily.py
+++ b/gnssrefl/subdaily.py
@@ -771,6 +771,9 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
     else:
         apply_if_corr = False
 
+    csvfile_spline = kwargs.get('csvfile_spline',False)
+    print('csvfile spline ', csvfile_spline)
+
     apply_rhdot  = kwargs.get('apply_rhdot',True)
     if apply_rhdot:
         apply_rhdot_corr = True
@@ -1134,7 +1137,7 @@ def rhdot_correction2(station,fname,fname_new,pltit,outlierV,outlierV2,**kwargs)
                                                fs, strsig, hires_figs,txtdir, ii,jj,th_even, spline_whole_time)
     H0 = sd.find_ortho_height(station,extension)
     # this writes out spline file and makes plot .... 
-    sd.RH_ortho_plot2( station, H0, year, txtdir, fs, th[jj],biasCor_rh[jj],gap_min_val,th,spline,delta_out)
+    sd.RH_ortho_plot2( station, H0, year, txtdir, fs, th[jj],biasCor_rh[jj],gap_min_val,th,spline,delta_out,csvfile_spline)
     print('\nRMS with frequency biases and RHdot taken out (m) ', np.round(newsigma,3) , '\n' )
 
 
diff --git a/gnssrefl/subdaily_cl.py b/gnssrefl/subdaily_cl.py
index dc9a1d132..bac296017 100644
--- a/gnssrefl/subdaily_cl.py
+++ b/gnssrefl/subdaily_cl.py
@@ -241,6 +241,9 @@ def subdaily(station: str, year: int, txtfile_part1: str = '', txtfile_part2: st
     if csvfile:
         print('>>>> WARNING: csvfile option is currently turned off.  We are working to add it back.')
         csvfile = False
+        csvfile_spline = True
+    else:
+        csvfile_spline = False
 
     outputs = [] # this is for multiple years
 
@@ -314,10 +317,13 @@ def subdaily(station: str, year: int, txtfile_part1: str = '', txtfile_part2: st
 
     # not sure why tv and corr are being returned.
     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,gap_min_val=gap_min_val,year=year,extension=extension,knots2=knots2)
+       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,
+                                      gap_min_val=gap_min_val,year=year, extension=extension, 
+                                      knots2=knots2,csvfile_spline=csvfile_spline)
        if plt:
            mplt.show()