diff --git a/.gitignore b/.gitignore
index ab7aae65..733dc353 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,6 +6,5 @@
 !README.md
 !version.txt
 !src/**.f90
-!ctem/setup.py
-!ctem/ctem/*.py
+!ctem/**.py
 !cmake/Modules/*.cmake
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index fd259936..eb4fdb32 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -29,7 +29,7 @@ ENDIF ()
 OPTION(USE_COARRAY "Use Coarray Fortran for parallelization of test particles" OFF)
 OPTION(USE_OPENMP "Use OpenMP for parallelization" ON)
 OPTION(USE_SIMD "Use SIMD vectorization" ON)
-OPTION(BUILD_SHARED_LIBS "Build using shared libraries" OFF)
+OPTION(BUILD_SHARED_LIBS "Build using shared libraries" ON)
 
 # The following section is modified from Numpy f2py documentation
 IF(PROJECT_SOURCE_DIR STREQUAL PROJECT_BINARY_DIR)
diff --git a/ctem/ctem/driver.py b/ctem/ctem/driver.py
index 0c9bd611..0249c19d 100644
--- a/ctem/ctem/driver.py
+++ b/ctem/ctem/driver.py
@@ -52,7 +52,7 @@ def __init__(self, param_file="ctem.in", isnew=True):
         }
 
         # Get the location of the CTEM executable
-        self.ctem_executable = Path(_pyfile).parent.parent.parent.parent / "build" / "src" / "CTEM"
+        self.ctem_executable = Path(_pyfile).parent.parent.parent / "bin" / "CTEM"
         if not self.ctem_executable.exists():
             print(f"CTEM driver not found at {self.ctem_executable}. Trying current directory.")
             self.ctem_executable = Path(currentdir) / "CTEM"
diff --git a/examples/global-lunar-bombardment/craterproduction.py b/examples/global-lunar-bombardment/craterproduction.py
deleted file mode 100644
index 4d81b687..00000000
--- a/examples/global-lunar-bombardment/craterproduction.py
+++ /dev/null
@@ -1,216 +0,0 @@
-import numpy as np
-from scipy import optimize
-
-# The Neukum production function for the Moon and Mars
-#Lunar PF from: Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86
-#Mars PF from: Ivanov (2001) SSR v. 96 pp. 87-104
-sfd_coef = {
-        "NPF_Moon" : [-3.0876,-3.557528,+0.781027,+1.021521,-0.156012,-0.444058,+0.019977,+0.086850,-0.005874,-0.006809,+8.25e-4, +5.54e-5],
-        "NPF_Mars" : [-3.384, -3.197,   +1.257,   +0.7915,  -0.4861,  -0.3630,  +0.1016,  +6.756e-2,-1.181e-2,-4.753e-3,+6.233e-4,+5.805e-5]
-    }
-sfd_range = {
-    "NPF_Moon" : [0.01,1000],
-    "NPF_Mars" : [0.015,362]
-}
-#Exponential time constant (Ga)
-tau = 6.93
-Nexp = 5.44e-14
-
-time_coef = {
-    "NPF_Moon" : Nexp,
-    "NPF_Mars" : Nexp * 10**(sfd_coef.get("NPF_Mars")[0]) / 10**(sfd_coef.get("NPF_Moon")[0])
-}
-
-def N1(T, pfmodel):
-    """
-    Return the N(1) value as a function of time for a particular production function model
-
-    Parameters
-    ----------
-    T : numpy array
-        Time in units of Ga
-    pfmodel : string
-        The production function model to use
-        Currently options are 'NPF_Moon' or 'NPF_Mars'
-
-    Returns
-    -------
-    N1 : numpy array
-        The number of craters per square kilometer greater than 1 km in diameter
-    """
-    T_valid_range = np.where((T <= 4.5) & (T >= 0.0), T, np.nan)
-    return time_coef.get(pfmodel) * (np.exp(tau * T_valid_range) - 1.0) + 10 ** (sfd_coef.get(pfmodel)[0]) * T_valid_range
-
-def CSFD(Dkm,planet):
-    """
-    Return the cumulative size-frequency distribution for a particular production function model
-
-    Parameters
-    ----------
-    Dkm : numpy array
-        Diameters in units of km
-    pfmodel : string
-        The production function model to use
-        Currently options are 'NPF_Moon' or 'NPF_Mars'
-
-    Returns
-    -------
-    CSFD : numpy array
-        The number of craters per square kilometer greater than Dkm in diameter at T=1 Ga
-    """
-    logCSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(sfd_coef.get(planet)))
-    return 10**(logCSFD)
-
-def DSFD(Dkm,planet):
-    """
-    Return the differential size-frequency distribution for a particular production function model
-
-    Parameters
-    ----------
-    Dkm : numpy array
-        Diameters in units of km
-    pfmodel : string
-        The production function model to use
-        Currently options are 'NPF_Moon' or 'NPF_Mars'
-
-    Returns
-    -------
-    DSFD : numpy array
-        The differential number of craters (dN/dD) per square kilometer greater than Dkm in diameter at T = 1 Ga
-    """
-    dcoef = sfd_coef.get(planet)[1:]
-    logDSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(dcoef))
-    return 10**(logDSFD) * CSFD(Dkm,planet) / Dkm
-
-def Tscale(T,planet):
-    """
-    Return the number density of craters at time T relative to time T = 1 Ga
-
-    Parameters
-    ----------
-    T : numpy array
-        Time in units of Ga
-    pfmodel : string
-        The production function model to use
-        Currently options are 'NPF_Moon' or 'NPF_Mars'
-
-    Returns
-    -------
-    Tscale : numpy array
-        N1(T) / CSFD(Dkm = 1.0)
-    """
-    return N1(T,planet) / CSFD(1.0,planet)
-
-def pf_csfd(T,Dkm,planet):
-    """
-    Return the cumulative size-frequency distribution for a particular production function model as a function of Time
-
-    Parameters
-    ----------
-    T : numpy array
-        Time in units of Ga
-    Dkm : numpy array
-        Diameters in units of km
-    pfmodel : string
-        The production function model to use
-        Currently options are 'NPF_Moon' or 'NPF_Mars'
-
-    Returns
-    -------
-    pf_csfd : numpy array
-        The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T
-    """
-    D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]),Dkm,np.nan)
-    return CSFD(D_valid_range,planet) * Tscale(T,planet)
-
-def pf_dsfd(T,Dkm,planet):
-    """
-    Return the differential size-frequency distribution for a particular production function model as a function of Time
-
-    Parameters
-    ----------
-    T : numpy array
-        Time in units of Ga
-    Dkm : numpy array
-        Diameters in units of km
-    pfmodel : string
-        The production function model to use
-        Currently options are 'NPF_Moon' or 'NPF_Mars'
-
-    Returns
-    -------
-    pf_dsfd : numpy array
-        The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T
-    """
-    D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]), Dkm, np.nan)
-    return DSFD(D_valid_range, planet) * Tscale(T, planet)
-
-def T_from_scale(TS,planet):
-    """
-    Return the time in Ga for the given number density of craters relative to that at 1 Ga.
-    This is the inverse of Tscale
-
-    Parameters
-    ----------
-    TS : numpy array
-        number density of craters relative to that at 1 Ga
-    pfmodel : string
-        The production function model to use
-        Currently options are 'NPF_Moon' or 'NPF_Mars'
-
-    Returns
-    -------
-    T_from_scale : numpy array
-        The time in Ga
-    """
-    def func(S):
-        return Tscale(S,planet) - TS
-    return optimize.fsolve(func, np.full_like(TS,4.4),xtol=1e-10)
-
-
-if __name__ == "__main__":
-    import matplotlib.pyplot as plt
-    import matplotlib.ticker as ticker
-    print("Tests go here")
-    print(f"T = 1 Ga, N(1) = {pf_csfd(1.0,1.00,'NPF_Moon')}")
-    print(f"T = 4.2 Ga, N(1) = {pf_csfd(4.2,1.00,'NPF_Moon')}")
-    print("Tscale test: Should return all 1s")
-    Ttest = np.logspace(-4,np.log10(4.4),num=100)
-    Tres = T_from_scale(Tscale(Ttest,'NPF_Mars'),'NPF_Mars')
-    print(Ttest / Tres)
-    #for i,t in enumerate(Ttest):
-    #    print(t,Tscale(t,'Moon'),Tres[i])
-
-    CSFDfig = plt.figure(1, figsize=(8, 7))
-    ax = {'NPF_Moon': CSFDfig.add_subplot(121),
-          'NPF_Mars': CSFDfig.add_subplot(122)}
-
-    tvals = [0.01,1.0,4.0]
-    x_min = 1e-3
-    x_max = 1e3
-    y_min = 1e-9
-    y_max = 1e3
-    Dvals = np.logspace(np.log10(x_min), np.log10(x_max))
-    for key in ax:
-        ax[key].title.set_text(key)
-        ax[key].set_xscale('log')
-        ax[key].set_yscale('log')
-        ax[key].set_ylabel('$\mathregular{N_{>D} (km^{-2})}$')
-        ax[key].set_xlabel('Diameter (km)')
-        ax[key].set_xlim(x_min, x_max)
-        ax[key].set_ylim(y_min, y_max)
-        ax[key].yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20))
-        ax[key].yaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100))
-        ax[key].xaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20))
-        ax[key].xaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100))
-        ax[key].grid(True,which="minor",ls="-",lw=0.5,zorder=5)
-        ax[key].grid(True,which="major",ls="-",lw=1,zorder=10)
-        for t in tvals:
-            prod = pf_csfd(t,Dvals,key)
-            ax[key].plot(Dvals, prod, '-', color='black', linewidth=1.0, zorder=50)
-            labeli = 15
-            ax[key].text(Dvals[labeli],prod[labeli],f"{t:.2f} Ga", ha="left", va="top",rotation=-72)
-
-    plt.tick_params(axis='y', which='minor')
-    plt.tight_layout()
-    plt.show()
diff --git a/examples/global-lunar-bombardment/ctem_driver.py b/examples/global-lunar-bombardment/ctem_driver.py
index e1c3650c..be37a458 100755
--- a/examples/global-lunar-bombardment/ctem_driver.py
+++ b/examples/global-lunar-bombardment/ctem_driver.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
-import ctem
-sim = ctem.Simulation()
+#!/usr/bin/env python3
+import ctem
+sim = ctem.Simulation()
 sim.run()
\ No newline at end of file
diff --git a/examples/global-lunar-bombardment/ctem_io_readers.py b/examples/global-lunar-bombardment/ctem_io_readers.py
deleted file mode 100644
index c825fdc9..00000000
--- a/examples/global-lunar-bombardment/ctem_io_readers.py
+++ /dev/null
@@ -1,146 +0,0 @@
-#!/usr/local/bin/python
-#
-#Cratered Terrain Evolution Model file reading utilities
-#
-#Original IDL design: Jim Richardson, Arecibo Observatory
-#Revised IDL design: David Minton, Purdue University
-#Re-engineered design and Python implementation: Matthew Route, Purdue University
-#August 2016
-#
-#Known issues for operation with CTEM
-#1) ctem.in has 1.0d0 value for maxcrat which is not readable by Python
-
-import numpy
-
-def read_ctemin(parameters,notset):
-    #Read and parse ctem.in file
-   inputfile = parameters['workingdir'] + parameters['ctemfile']
-
-   #Read ctem.in file
-   print('Reading input file '+ parameters['ctemfile'])
-   fp = open(inputfile,'r')
-   lines = fp.readlines()
-   fp.close()
-   
-   #Process file text
-   for line in lines:
-      fields = line.split()
-      if len(fields) > 0:
-         if ('pix' == fields[0].lower()): parameters['pix']=float(fields[1])
-         if ('gridsize' == fields[0].lower()): parameters['gridsize']=int(fields[1])
-         if ('seed' == fields[0].lower()): parameters['seed']=int(fields[1])
-         if ('sfdfile' == fields[0].lower()): parameters['sfdfile']=fields[1]
-         if ('impfile' == fields[0].lower()): parameters['impfile']=fields[1]
-         if ('maxcrat' == fields[0].lower()): parameters['maxcrat']=float(fields[1])
-         if ('sfdcompare' == fields[0].lower()): parameters['sfdcompare']=fields[1]
-         if ('interval' == fields[0].lower()): parameters['interval']=float(fields[1])
-         if ('numintervals' == fields[0].lower()): parameters['numintervals']=int(fields[1])
-         if ('popupconsole' == fields[0].lower()): parameters['popupconsole']=fields[1]
-         if ('saveshaded' == fields[0].lower()): parameters['saveshaded']=fields[1]
-         if ('saverego' == fields[0].lower()): parameters['saverego']=fields[1]
-         if ('savepres' == fields[0].lower()): parameters['savepres']=fields[1]
-         if ('savetruelist' == fields[0].lower()): parameters['savetruelist']=fields[1]
-         if ('runtype' == fields[0].lower()): parameters['runtype']=fields[1]
-         if ('restart' == fields[0].lower()): parameters['restart']=fields[1]
-         if ('quasimc' == fields[0].lower()): parameters['quasimc']=fields[1]
-         if ('realcraterlist' == fields[0].lower()): parameters['realcraterlist']=fields[1]
-         if ('shadedminh' == fields[0].lower()):
-            parameters['shadedminh'] = float(fields[1])
-            parameters['shadedminhdefault'] = 0
-         if ('shadedmaxh' == fields[0].lower()):
-            parameters['shadedmaxh'] = float(fields[1])
-            parameters['shadedmaxhdefault'] = 0
-       
-   #Test values for further processing
-   if (parameters['interval'] <= 0.0):
-       print('Invalid value for or missing variable INTERVAL in '+ inputfile)
-   if (parameters['numintervals'] <= 0):
-       print('Invalid value for or missing variable NUMINTERVALS in '+ inputfile)
-   if (parameters['pix'] <= 0.0):
-       print('Invalid value for or missing variable PIX in '+ inputfile)
-   if (parameters['gridsize'] <= 0):
-       print('Invalid value for or missing variable GRIDSIZE in '+ inputfile)
-   if (parameters['seed'] == 0):
-       print('Invalid value for or missing variable SEED in '+ inputfile)
-   if (parameters['sfdfile'] == notset):
-       print('Invalid value for or missing variable SFDFILE in '+ inputfile)
-   if (parameters['impfile'] == notset):
-       print('Invalid value for or missing variable IMPFILE in '+ inputfile)
-   if (parameters['popupconsole'] == notset):
-       print('Invalid value for or missing variable POPUPCONSOLE in '+ inputfile)
-   if (parameters['saveshaded'] == notset):
-       print('Invalid value for or missing variable SAVESHADED in '+ inputfile)
-   if (parameters['saverego'] == notset):
-       print('Invalid value for or missing variable SAVEREGO in '+ inputfile)
-   if (parameters['savepres'] == notset):
-       print('Invalid value for or missing variable SAVEPRES in '+ inputfile)
-   if (parameters['savetruelist'] == notset):
-       print('Invalid value for or missing variable SAVETRUELIST in '+ inputfile)
-   if (parameters['runtype'] == notset):
-       print('Invalid value for or missing variable RUNTYPE in '+ inputfile)
-   if (parameters['restart'] == notset):
-       print('Invalid value for or missing variable RESTART in '+ inputfile)
-       
-   return
-
-def read_formatted_ascii(filename, skip_lines):
-   #Generalized ascii text reader
-   #For use with sfdcompare, production, odist, tdist, pdist data files
-   data = numpy.genfromtxt(filename, skip_header = skip_lines)
-   return data
-
-def read_unformatted_binary(filename, gridsize):
-   #Read unformatted binary files created by Fortran
-   #For use with surface ejecta and surface dem data files
-   dt = numpy.float
-   data = numpy.fromfile(filename, dtype = dt)
-   data.shape = (gridsize,gridsize)
-   
-   return data
-
-def read_ctemdat(parameters, seedarr):
-    #Read and parse ctem.dat file
-    datfile = parameters['workingdir'] + 'ctem.dat'
- 
-    #Read ctem.dat file    
-    print('Reading input file '+ parameters['datfile'])
-    fp = open(datfile,'r')
-    lines = fp.readlines()
-    fp.close()
-   
-    #Parse file lines and update parameter fields
-    fields = lines[0].split()
-    if len(fields) > 0:
-      parameters['totalimpacts'] = float(fields[0])
-      parameters['ncount'] = int(fields[1])
-      parameters['curyear'] = float(fields[2])
-      parameters['restart'] = fields[3]
-      parameters['fracdone'] = float(fields[4])
-      parameters['masstot'] = float(fields[5])
-
-    #Parse remainder of file to build seed array
-    nlines = len(lines)
-    index = 1
-    while (index < nlines):
-        fields = lines[index].split()
-        seedarr[index - 1] = float(fields[0])
-        index += 1
-
-    parameters['seedn'] = index - 1
-
-    return
-
-def read_impact_mass(filename):
-   #Read impact mass file
-   
-   fp=open(filename,'r')
-   line=fp.readlines()
-   fp.close()
-   
-   fields = line[0].split()
-   if (len(fields) > 0):
-      mass = float(fields[0])
-   else:
-      mass = 0
-      
-   return mass
\ No newline at end of file
diff --git a/examples/global-lunar-bombardment/ctem_io_writers.py b/examples/global-lunar-bombardment/ctem_io_writers.py
deleted file mode 100644
index e7a2192c..00000000
--- a/examples/global-lunar-bombardment/ctem_io_writers.py
+++ /dev/null
@@ -1,282 +0,0 @@
-#!/usr/local/bin/python
-#
-#Cratered Terrain Evolution Model file and graphical writing utilities
-#
-#Original IDL design: Jim Richardson, Arecibo Observatory
-#Revised IDL design: David Minton, Purdue University
-#Re-engineered design and Python implementation: Matthew Route, Purdue University
-#August 2016
-#
-
-import matplotlib
-from matplotlib import pyplot
-import numpy
-import os
-import shutil
-import scipy
-from scipy import signal
-
-#Set pixel scaling common for image writing, at 1 pixel/ array element
-dpi = 72.0
-
-#Write production function to file production.dat
-#This file format does not exactly match that generated from IDL. Does it work?
-def write_production(parameters, production):
-    filename = parameters['workingdir'] + parameters['sfdfile']
-    numpy.savetxt(filename, production, fmt='%1.8e', delimiter='   ')
-    
-    return
-
-def create_dir_structure(parameters):
-    #Create directories for various output files if they do not already exist
-    directories=['dist','misc','rego','rplot','surf','shaded']
-    
-    for directory in directories:
-        dir_test = parameters['workingdir'] + directory
-        if not os.path.isdir(dir_test):
-            os.makedirs(dir_test)
-    
-    return
-
-def write_ctemdat(parameters, seedarr):
-    #Write various parameters and random number seeds into ctem.dat file
-    filename = parameters['workingdir'] + parameters['datfile']
-    fp = open(filename,'w')
-    
-    template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n"
-    fp.write(template % parameters)
-    
-    #Write random number seeds to the file
-    for index in range(parameters['seedn']):
-       fp.write("%12d\n" % seedarr[index])
-       
-    fp.close()   
-    
-    return
-
-def copy_dists(parameters):
-    #Save copies of distribution files
-        
-    orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution']
-    dest_list = ['odist', 'ocum', 'pdist', 'tdist']       
-        
-    for index in range(len(orig_list)):  
-        forig = parameters['workingdir'] + orig_list[index] + '.dat'
-        fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount']
-        shutil.copy2(forig, fdest)
-    
-    forig = parameters['workingdir'] + 'impactmass.dat'
-    fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount']
-    shutil.copy2(forig, fdest)   
-    
-    if (parameters['savetruelist'].upper() == 'T'):
-        forig = parameters['workingdir'] + 'tcumulative.dat'
-        fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount']
-        shutil.copy2(forig, fdest)
-    
-    return
-
-#Possible references
-#http://nbviewer.jupyter.org/github/ThomasLecocq/geophysique.be/blob/master/2014-02-25%20Shaded%20Relief%20Map%20in%20Python.ipynb   
-   
-def image_dem(parameters, surface_dem):
-    
-    #Create surface dem map
-    solar_angle = 20.0
-    dem_map = numpy.copy(surface_dem) - numpy.roll(surface_dem, 1, 0)
-    dem_map = (0.5 * numpy.pi) + numpy.arctan2(dem_map, parameters['pix'])
-    dem_map = dem_map - numpy.radians(solar_angle) * (0.5 * numpy.pi)
-    numpy.place(dem_map, dem_map > (0.5 * numpy.pi), 0.5  *numpy.pi)
-    dem_map = numpy.absolute(dem_map)
-    dem_map = 254.0 * numpy.cos(dem_map)
-    
-    #Save image to file
-    filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount']
-    height = parameters['gridsize'] / dpi
-    width = height
-    fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi)
-    fig.figimage(dem_map, cmap = matplotlib.cm.gray, origin = 'lower')
-    matplotlib.pyplot.savefig(filename)    
-    
-    return   
-   
-def image_regolith(parameters, regolith):
-    
-    #Create scaled regolith image
-    minref = parameters['pix'] * 1.0e-4
-    maxreg = numpy.amax(regolith)
-    minreg = numpy.amin(regolith)
-    if (minreg < minref): minreg = minref
-    if (maxreg < minref): maxreg = (minref + 1.0e3)
-    regolith_scaled = numpy.copy(regolith)
-    numpy.place(regolith_scaled, regolith_scaled < minref, minref)
-    regolith_scaled = 254.0 * ((numpy.log(regolith_scaled) - numpy.log(minreg)) / (numpy.log(maxreg) - numpy.log(minreg)))    
-    
-    #Save image to file
-    filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount']
-    height = parameters['gridsize'] / dpi
-    width = height
-    fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi)
-    fig.figimage(regolith_scaled, cmap = matplotlib.cm.nipy_spectral, origin = 'lower')
-    matplotlib.pyplot.savefig(filename)   
-   
-    return
-    
-def image_shaded_relief(parameters, surface_dem):
-    #The color scale and appearance of this do not quite match the IDL version
-
-    #Create image by convolving DEM with 3x3 illumination matrix
-    light = numpy.array([[1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [-1.0, -1.0, -1.0]])
-    convolved_map = scipy.signal.convolve2d(surface_dem, light, mode = 'same')
-    
-    #Adjust output to resemble IDL (north slopes illuminated, south in shadow)
-    convolved_map = convolved_map * -1.0
-    convolved_map[0,:]=0.0
-    convolved_map[-1,:] =0.0
-    convolved_map[:,0]=0.0
-    convolved_map[:,-1]=0.0    
-
-    #If no shadedmin/max parameters are read in from ctem.dat, determine the values from the data
-    if (parameters['shadedminhdefault'] == 1): shadedminh = numpy.amin(surface_dem)
-    if (parameters['shadedmaxhdefault'] == 1): shadedmaxh = numpy.amax(surface_dem)
-
-    #If min and max appear to be reversed, then fix them
-    if (parameters['shadedminh'] > parameters['shadedmaxh']):
-        temp = parameters['shadedminh']
-        shadedminh = parameters['shadedmaxh']
-        shadedmaxh = temp
-    else:
-        shadedminh = parameters['shadedminh']
-        shadedmaxh = parameters['shadedmaxh']
-            
-    #If dynamic range is valid, construct a shaded DEM
-    dynamic_range = shadedmaxh - shadedminh
-    if  (dynamic_range != 0):
-        dem_scaled = numpy.copy(surface_dem) - shadedminh
-        numpy.place(dem_scaled, dem_scaled < 0.0, 0.0)
-        numpy.place(dem_scaled, dem_scaled > dynamic_range, dynamic_range)
-        dem_scaled = dem_scaled / dynamic_range
-    else:
-        dem_scaled = numpy.copy(surface_dem) * 0.0
-        
-    #Generate shaded depth map with surface_dem color scaling (RGBA)
-    shaded = scipy.misc.bytescale(convolved_map)
-    if numpy.amax(shaded) == 0:  shaded=255
-    shadedscl = shaded / 255.0
-    
-    #shaded_imagearr = matplotlib.cm.jet(dem_scaled)
-    #print dem_scaled[0:4,0:4]
-    #shaded_imagearr[:,:,0] = shaded_imagearr[:,:,0] * shadedscl
-    #shaded_imagearr[:,:,1] = shaded_imagearr[:,:,1] * shadedscl
-    #shaded_imagearr[:,:,2] = shaded_imagearr[:,:,2] * shadedscl
-    #shaded_imagearr[:,:,3] = shaded_imagearr[:,:,3] * shadedscl
-    
-    #Delivers nearly proper coloring, but no shaded relief    
-    shaded_imagearr = dem_scaled * shadedscl
-    shaded_imagearr = matplotlib.cm.jet(shaded_imagearr)
-    shaded_imagearr = numpy.around(shaded_imagearr, decimals = 1)
-        
-    #Save image to file
-    filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount']
-    height = parameters['gridsize'] / dpi
-    width = height
-    fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi)
-    fig.figimage(shaded_imagearr, cmap = matplotlib.cm.jet, origin = 'lower')
-    matplotlib.pyplot.savefig(filename)    
-    return
-    
-def create_rplot(parameters,odist,pdist,tdist,ph1):
-    #Parameters: empirical saturation limit and dfrac
-    satlimit = 3.12636
-    dfrac = 2**(1./4) * 1.0e-3
-    
-    #Calculate geometric saturation
-    minx = (parameters['pix'] / 3.0) * 1.0e-3 
-    maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3     
-    geomem = numpy.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]])
-    geomep = numpy.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]])
-    
-    #Create distribution arrays without zeros for plotting on log scale
-    idx = numpy.nonzero(odist[:,5])
-    odistnz = odist[idx]
-    odistnz = odistnz[:,[2,5]]
-    
-    idx = numpy.nonzero(pdist[:,5])
-    pdistnz = pdist[idx]
-    pdistnz = pdistnz[:,[2,5]]
-    
-    idx = numpy.nonzero(tdist[:,5])
-    tdistnz = tdist[idx]
-    tdistnz = tdistnz[:,[2,5]]
-    
-    #Correct pdist
-    pdistnz[:,1] = pdistnz[:,1] * parameters['curyear'] / parameters['interval']    
-    
-    #Create sdist bin factors, which contain one crater per bin
-    area = (parameters['gridsize'] * parameters['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 = numpy.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" % parameters['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 = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount']
-    height = parameters['gridsize'] / dpi
-    width = height
-    fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi)
-    
-    #Alter background color to be black, and change axis colors accordingly
-    matplotlib.pyplot.style.use('dark_background')
-    matplotlib.pyplot.rcParams['axes.prop_cycle']
-    
-    #Plot data    
-    matplotlib.pyplot.plot(odistnz[:,0]*1.0e-3, odistnz[:,1], linewidth=3.0, color = 'blue')
-    matplotlib.pyplot.plot(pdistnz[:,0]*1.0e-3, pdistnz[:,1], linewidth=2.0, linestyle='dashdot', color = 'white')
-    matplotlib.pyplot.plot(tdistnz[:,0]*1.0e-3, tdistnz[:,1], linewidth=2.0, color = 'red')
-    
-    matplotlib.pyplot.plot(geomem[:,0], geomem[:,1], linewidth=2.0, linestyle =':', color = 'yellow')
-    matplotlib.pyplot.plot(geomep[:,0], geomep[:,1], linewidth=2.0, linestyle =':', color = 'yellow')
-    matplotlib.pyplot.plot(sdist[:,0], sdist[:,1], linewidth=2.0, linestyle =':', color = 'yellow')
-    
-    matplotlib.pyplot.plot(ph1[:,0] * dfrac, ph1[:,1], 'wo')    
-    
-    #Create plot labels    
-    matplotlib.pyplot.title('Crater Distribution R-Plot',fontsize=22)
-    matplotlib.pyplot.xlim([minx, maxx])
-    matplotlib.pyplot.xscale('log')
-    matplotlib.pyplot.xlabel('Crater Diameter (km)',fontsize=18)
-    matplotlib.pyplot.ylim([5.0e-4, 5.0])    
-    matplotlib.pyplot.yscale('log')
-    matplotlib.pyplot.ylabel('R Value', fontsize=18)
-    matplotlib.pyplot.text(1.0e-2, 1.0, timelabel, fontsize=18)
-    
-    matplotlib.pyplot.tick_params(axis='both', which='major', labelsize=14)
-    matplotlib.pyplot.tick_params(axis='both', which='minor', labelsize=12)
-    
-    matplotlib.pyplot.savefig(filename)   
-   
-    return
-
-def write_realcraters(parameters, realcraters):
-    """writes file of real craters for use in quasi-MC runs"""
-
-    filename = parameters['workingdir'] + 'craterlist.dat'
-    numpy.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t')
-
-    return
\ No newline at end of file
diff --git a/examples/morphology_test_cases/global/ctem_io_writers.py b/examples/morphology_test_cases/global/ctem_io_writers.py
index e7a2192c..d07d7d54 100644
--- a/examples/morphology_test_cases/global/ctem_io_writers.py
+++ b/examples/morphology_test_cases/global/ctem_io_writers.py
@@ -273,10 +273,3 @@ def create_rplot(parameters,odist,pdist,tdist,ph1):
    
     return
 
-def write_realcraters(parameters, realcraters):
-    """writes file of real craters for use in quasi-MC runs"""
-
-    filename = parameters['workingdir'] + 'craterlist.dat'
-    numpy.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t')
-
-    return
\ No newline at end of file
diff --git a/src/io/io_updatePbar.f90 b/src/io/io_updatePbar.f90
index 7ba679d9..eab51800 100644
--- a/src/io/io_updatePbar.f90
+++ b/src/io/io_updatePbar.f90
@@ -42,17 +42,17 @@ subroutine io_updatePbar(message)
 do k = 1,endpoint !- 1 
    bar(6+k:6+k)=pbarchar(k:k)
 end do
-!select case(flip)
-!case(1)
-!   bar(6+endpoint:6+endpoint)="/"
-!case(2)
-!   bar(6+endpoint:6+endpoint)="-"
-!case(3)
-!   bar(6+endpoint:6+endpoint)="\"
-!case(4)
-!   bar(6+endpoint:6+endpoint)="|"
-!end select
-!flip = flip + 1
+select case(flip)
+case(1)
+  bar(6+endpoint:6+endpoint)="/"
+case(2)
+  bar(6+endpoint:6+endpoint)="-"
+case(3)
+  bar(6+endpoint:6+endpoint)="\"
+case(4)
+  bar(6+endpoint:6+endpoint)="|"
+end select
+flip = flip + 1
 if (flip > 4) flip = 1
 ! print the progress bar.  
 write(fmtlabel,'("(A1,A",I2.2,",1X,A",I2.2,",$)")') PBARSIZE+7,MESSAGESIZE