diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
index 8f565f1e..832eb5b0 100644
--- a/python/ctem/ctem/driver.py
+++ b/python/ctem/ctem/driver.py
@@ -22,7 +22,6 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'popupconsole': None,
             'saveshaded': None,
             'saverego': None,
-            'savepres': None,
             'savetruelist': None,
             'seedn': 1,
             'totalimpacts': 0,
@@ -76,7 +75,13 @@ def __init__(self, param_file="ctem.in", isnew=True):
 
         for k, v in self.output_filenames.items():
             self.output_filenames[k] = os.path.join(currentdir, v)
-        
+
+        self.directories = ['dist', 'misc', 'surf']
+        if self.user['saveshaded'].upper() == 'T':
+            self.directories.append('shaded')
+        if self.user['saverego'].upper() == 'T':
+            self.directories.append('rego')
+            
         # Set up data arrays
         self.seedarr = np.zeros(100, dtype=int)
         self.seedarr[0] = self.user['seed']
@@ -85,6 +90,7 @@ def __init__(self, param_file="ctem.in", isnew=True):
         self.tdist = np.zeros([1, 6])
         self.surface_dem = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
         self.surface_ejc = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
+        self.ph1 = None
         
         if self.user['sfdcompare'] is not None:
             # Read sfdcompare file
@@ -98,7 +104,7 @@ def __init__(self, param_file="ctem.in", isnew=True):
             if (self.user['restart'].upper() == 'F'):
                 print('Starting a new run')
 
-                util.create_dir_structure(self.user)
+                util.create_dir_structure(self.user, self.directories)
                 # Delete any old output files
                 for k, v in self.output_filenames.items():
                     if os.path.isfile(v):
@@ -260,7 +266,7 @@ def process_output(self):
         # Display results
         print(self.user['ncount'], '  Generating surface images and plots')
 
-        # Write surface dem, surface ejecta, shaded relief, and rplot data
+        # Write surface dem, surface ejecta and shaded relief data
         util.image_dem(self.user, self.surface_dem)
         if (self.user['saverego'].upper() == 'T'):
             util.image_regolith(self.user, self.surface_ejc)
@@ -268,8 +274,6 @@ def process_output(self):
             util.image_shaded_relief(self.user, self.surface_dem)
         
         if self.user['ncount'] > 0: # These aren't available yet from the initial conditions
-            if (self.user['savepres'].upper() == 'T'):
-                util.create_rplot(self.user, self.odist, self.pdist, self.tdist, self.ph1)
             
             # Save copy of crater distribution files
             # Update user: mass, curyear, regolith properties
@@ -318,7 +322,7 @@ def cleanup(self):
         """
         # This is a list of files generated by the main Fortran program
         print("Deleting all files generated by CTEM")
-        util.destroy_dir_structure(self.user)
+        util.destroy_dir_structure(self.user, self.directories)
         for key, filename in self.output_filenames.items():
             print(f"Deleting file {filename}")
             try:
diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py
index 8cab22e3..df0465a1 100644
--- a/python/ctem/ctem/util.py
+++ b/python/ctem/ctem/util.py
@@ -9,12 +9,12 @@
 dpi = 72.0
 
 # These are directories that are created by CTEM in order to store intermediate results
-directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded']
 
-def create_dir_structure(user):
+def create_dir_structure(user, directories):
     """
     Create directories for various output files if they do not already exist
     """
+
     for dirname in directories:
         dirpath = os.path.join(user['workingdir'], dirname)
         if not os.path.isdir(dirpath):
@@ -23,7 +23,7 @@ def create_dir_structure(user):
     
     return
 
-def destroy_dir_structure(user):
+def destroy_dir_structure(user, directories):
     """
     Deletes directories generated by create_dir_structure
     """
@@ -35,95 +35,6 @@ def destroy_dir_structure(user):
 
     return
 
-
-def create_rplot(user, odist, pdist, tdist, ph1):
-    # Parameters: empirical saturation limit and dfrac
-    satlimit = 3.12636
-    dfrac = 2 ** (1. / 4) * 1.0e-3
-    
-    # Calculate geometric saturation
-    minx = (user['pix'] / 3.0) * 1.0e-3
-    maxx = 3 * user['pix'] * user['gridsize'] * 1.0e-3
-    geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]])
-    geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]])
-    
-    # Create distribution arrays without zeros for plotting on log scale
-    idx = np.nonzero(odist[:, 5])
-    odistnz = odist[idx]
-    odistnz = odistnz[:, [2, 5]]
-    
-    idx = np.nonzero(pdist[:, 5])
-    pdistnz = pdist[idx]
-    pdistnz = pdistnz[:, [2, 5]]
-    
-    idx = np.nonzero(tdist[:, 5])
-    tdistnz = tdist[idx]
-    tdistnz = tdistnz[:, [2, 5]]
-    
-    # Correct pdist
-    pdistnz[:, 1] = pdistnz[:, 1] * user['curyear'] / user['interval']
-    
-    # Create sdist bin factors, which contain one crater per bin
-    area = (user['gridsize'] * user['pix'] * 1.0e-3) ** 2.
-    plo = 1
-    sq2 = 2 ** (1. / 2)
-    while (sq2 ** plo > minx):
-        plo = plo - 1
-    phi = plo + 1
-    while (sq2 ** phi < maxx):
-        phi = phi + 1
-    n = phi - plo + 1
-    sdist = np.zeros([n, 2])
-    p = plo
-    for index in range(n):
-        sdist[index, 0] = sq2 ** p
-        sdist[index, 1] = sq2 ** (2.0 * p + 1.5) / (area * (sq2 - 1))
-        p = p + 1
-    
-    # Create time label
-    tlabel = "%5.4e" % user['curyear']
-    tlabel = tlabel.split('e')
-    texp = str(int(tlabel[1]))
-    timelabel = 'Time = ' + r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs'
-    
-    # Save image to file
-    filename = os.path.join(user['workingdir'],'rplot',"rplot%06d.png" % user['ncount'])
-    height = user['gridsize'] / dpi
-    width = height
-    fig = plt.figure(figsize=(width, height), dpi=dpi)
-    
-    # Alter background color to be black, and change axis colors accordingly
-    plt.style.use('dark_background')
-    plt.rcParams['axes.prop_cycle']
-    
-    # Plot data
-    plt.plot(odistnz[:, 0] * 1.0e-3, odistnz[:, 1], linewidth=3.0, color='blue')
-    plt.plot(pdistnz[:, 0] * 1.0e-3, pdistnz[:, 1], linewidth=2.0, linestyle='dashdot', color='white')
-    plt.plot(tdistnz[:, 0] * 1.0e-3, tdistnz[:, 1], linewidth=2.0, color='red')
-    
-    plt.plot(geomem[:, 0], geomem[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    plt.plot(geomep[:, 0], geomep[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    plt.plot(sdist[:, 0], sdist[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    
-    plt.plot(ph1[:, 0] * dfrac, ph1[:, 1], 'wo')
-    
-    # Create plot labels
-    plt.title('Crater Distribution R-Plot', fontsize=22)
-    plt.xlim([minx, maxx])
-    plt.xscale('log')
-    plt.xlabel('Crater Diameter (km)', fontsize=18)
-    plt.ylim([5.0e-4, 5.0])
-    plt.yscale('log')
-    plt.ylabel('R Value', fontsize=18)
-    plt.text(1.0e-2, 1.0, timelabel, fontsize=18)
-    
-    plt.tick_params(axis='both', which='major', labelsize=14)
-    plt.tick_params(axis='both', which='minor', labelsize=12)
-    
-    plt.savefig(filename)
-    
-    return
-
 def image_dem(user, DEM):
     dpi = 300.0  # 72.0
     pix = user['pix']
@@ -166,7 +77,7 @@ def image_regolith(user, regolith):
     height = user['gridsize'] / dpi
     width = height
     fig = plt.figure(figsize=(width, height), dpi=dpi)
-    fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower')
+    fig.figimage(regolith_scaled, cmap=cm.cividis, origin='lower')
     plt.savefig(filename)
     
     return