diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index 532e6fc..048b71f 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -74,7 +74,7 @@ jobs:
       - name: Test with pytest
         run: |
           export FWL_DATA="/home/runner/work/fwl_data"
-          export SOCRATES="/home/runner/work/JANUS/JANUS/SOCRATES"
+          export RAD_DIR="/home/runner/work/JANUS/JANUS/SOCRATES"
           coverage run -m pytest
 
       - name: Report coverage
diff --git a/.gitignore b/.gitignore
index 1d8e60d..b3575d8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,5 @@ utils/*.ipynb*
 build/
 fwl_data*/
 dist/
+.pytest_cache
+
diff --git a/docs/getting_started.md b/docs/getting_started.md
index de956ad..00230fc 100644
--- a/docs/getting_started.md
+++ b/docs/getting_started.md
@@ -41,7 +41,7 @@ If you install and compile [SOCRATES](https://github.com/nichollsh/SOCRATES) you
 you can override the path using the `SOCRATES` environment variable, e.g.
 
 ```console
-SOCRATES=/home/user/path/to/SOCRATES pytest
+RAD_DIR=/home/user/path/to/SOCRATES pytest
 ```
 
 ### `FWL_DATA`
diff --git a/examples/SocRadConv.py b/examples/SocRadConv.py
index be14184..8a5c780 100755
--- a/examples/SocRadConv.py
+++ b/examples/SocRadConv.py
@@ -20,6 +20,9 @@
 import numpy as np
 from importlib.resources import files
 
+from janus.utils.logs import setup_logger
+log = setup_logger()
+
 from janus.modules import RadConvEqm, plot_fluxes, plot_emission
 from janus.utils import atmos, CleanOutputDir, DownloadSpectralFiles, DownloadStellarSpectra, plot_adiabats, ReadBandEdges, StellarSpectrum
 import mors
@@ -33,7 +36,7 @@
 ####################################
 if __name__ == "__main__":
 
-    print("Start JANUS")
+    log.info("Start JANUS")
 
     # Set up dirs
     if os.environ.get('RAD_DIR') == None:
@@ -83,7 +86,7 @@
     StellarSpectrum.PrepareStellarSpectrum(spec.wl, spec.fl, socstar)
 
     # Move/prepare spectral file
-    print("Inserting stellar spectrum")
+    log.info("Inserting stellar spectrum")
     StellarSpectrum.InsertStellarSpectrum(
         str(FWL_DATA_DIR / 'spectral_files'/'Oak'/'318'/'Oak.sf'),
         socstar,
@@ -99,10 +102,10 @@
     if cfg['star']['stellar_heating'] == False:
         atm.instellation = 0.
     else:
-        mors.DownloadEvolutionTracks("/Baraffe")
+        mors.DownloadEvolutionTracks("Baraffe")
         baraffe = mors.BaraffeTrack(star_mass)
-        atm.instellation = baraffe.BaraffeSolarConstant(time['star'], mean_distance)
-        print("Instellation:", round(atm.instellation), "W/m^2")
+        atm.instellation = baraffe.BaraffeSolarConstant(time['star'], mean_distance) 
+        log.info("Instellation: %.2e W/m^2" % atm.instellation)
 
     # Set up atmosphere with general adiabat
     atm_dry, atm = RadConvEqm(dirs,
@@ -132,4 +135,5 @@
     CleanOutputDir(dirs['output'])
 
     end = t.time()
-    print("Runtime:", round(end - start,2), "s")
+    log.info("Runtime: %.2f s" % float(end-start))
+
diff --git a/examples/demo_instellation.py b/examples/demo_instellation.py
index ce6a565..6122173 100755
--- a/examples/demo_instellation.py
+++ b/examples/demo_instellation.py
@@ -11,6 +11,9 @@
 import os, shutil, toml
 import numpy as np
 
+from janus.utils.logs import setup_logger
+log = setup_logger()
+
 from janus.modules import MCPA_CBL
 from janus.utils import atmos, CleanOutputDir, DownloadSpectralFiles, DownloadStellarSpectra, ReadBandEdges, StellarSpectrum
 
@@ -20,7 +23,7 @@
 
 if __name__=='__main__':
 
-    print("Start")
+    log.info("Start")
 
     # Set up dirs
     if os.environ.get('RAD_DIR') == None:
@@ -37,7 +40,7 @@
     os.mkdir(dirs["output"])
 
     #Download required spectral files
-    DownloadSpectralFiles("/Oak")
+    DownloadSpectralFiles("Oak")
     DownloadStellarSpectra()
 
     # Read spectrum
@@ -50,7 +53,7 @@
 
 
     # Setup spectral file
-    print("Inserting stellar spectrum")
+    log.info("Inserting stellar spectrum")
     StellarSpectrum.InsertStellarSpectrum(
         str(FWL_DATA_DIR / 'spectral_files'/'Oak'/'318'/'Oak.sf'),
         socstar,
@@ -66,7 +69,7 @@
     # Star luminosity
     time = { "planet": cfg['planet']['time'], "star": cfg['star']['time']}
     star_mass = cfg['star']['star_mass']
-    mors.DownloadEvolutionTracks("/Baraffe")
+    mors.DownloadEvolutionTracks("Baraffe")
     baraffe = mors.BaraffeTrack(star_mass)
 
     # Define volatiles by mole fractions
@@ -94,7 +97,7 @@
     ts_arr  = []    # surface temperature
     tr_arr  = []    # tropopause temperature
     for i in range(7):
-      print("Orbital separation = %.2f AU" % r_arr[i])
+      log.info("Orbital separation = %.2f AU" % r_arr[i])
 
       atm.instellation = baraffe.BaraffeSolarConstant(time['star'], r_arr[i])
       atmos.setTropopauseTemperature(atm)
@@ -121,14 +124,14 @@
       # Save netcdf
       atm.write_ncdf(dirs["output"]+"/profile%.2fAU.nc"%r_arr[i])
 
-      print(" ")
+      log.info(" ")
 
     save_arr = [r_arr, asf_arr, OLR_arr, net_arr, ts_arr, tr_arr]
     np.savetxt(dirs["output"]+"/data_%dK.csv"%T_magma,
                np.array(save_arr).T, fmt="%.5e", delimiter=",",
                header="r [AU],  S_0 [W m-2], OLR [W m-2], net [W m-2], ts [K], tr[K] ")
-
-    print("Making plots")
+    
+    log.info("Making plots")
 
     plt.ioff()
     fig,(ax1, ax2) = plt.subplots(2,1, figsize=(6,6.4))
@@ -188,4 +191,5 @@
     CleanOutputDir(dirs['output'])
 
     # Done
-    print("Done!")
+    log.info("Done!")
+
diff --git a/examples/demo_runaway_greenhouse.py b/examples/demo_runaway_greenhouse.py
index e77f914..8fcdfe5 100755
--- a/examples/demo_runaway_greenhouse.py
+++ b/examples/demo_runaway_greenhouse.py
@@ -10,6 +10,9 @@
 from matplotlib.ticker import MultipleLocator
 from importlib.resources import files
 
+from janus.utils.logs import setup_logger
+log = setup_logger()
+
 from janus.modules import RadConvEqm
 from janus.utils import atmos, CleanOutputDir, DownloadSpectralFiles, DownloadStellarSpectra, ReadBandEdges, StellarSpectrum
 import mors
@@ -18,8 +21,8 @@
 
 if __name__=='__main__':
 
-    print("Start")
-    print(" ")
+    log.info("Start")
+    log.info(" ")
 
     # Set up dirs
     if os.environ.get('RAD_DIR') == None:
@@ -36,7 +39,7 @@
     os.mkdir(dirs["output"])
 
     #Download required spectral files
-    DownloadSpectralFiles("/Oak")
+    DownloadSpectralFiles("Oak")
     DownloadStellarSpectra()
 
     # Read spectrum
@@ -49,13 +52,13 @@
 
 
     # Setup spectral file
-    print("Inserting stellar spectrum")
+    log.info("Inserting stellar spectrum")
     StellarSpectrum.InsertStellarSpectrum(
         str(FWL_DATA_DIR / 'spectral_files'/'Oak'/'318'/'Oak.sf'),
         socstar,
         dirs["output"]
     )
-    print(" ")
+    log.info(" ")
     band_edges = ReadBandEdges(dirs["output"]+"star.sf")
 
     # Open config file
@@ -77,24 +80,24 @@
     atm = atmos.from_file(cfg_file, band_edges, vol_mixing=vol_mixing, vol_partial={})
 
     # Compute stellar heating
-    mors.DownloadEvolutionTracks("/Baraffe")
+    mors.DownloadEvolutionTracks("Baraffe")
     baraffe = mors.BaraffeTrack(star_mass)
     atm.instellation = baraffe.BaraffeSolarConstant(time['star'], mean_distance)
 
     #Run Janus
 
     # Run JANUS in a loop to generate runaway curve
-    print("Running JANUS...")
+    log.info("Running JANUS...")
     Ts_arr = np.linspace(200, 2800, 20)
     OLR_arr = []
     for i in range(20):
-      print("T_surf = %d K" % Ts_arr[i])
+      log.info("T_surf = %d K" % Ts_arr[i])
       atmos.setSurfaceTemperature(atm, Ts_arr[i])
 
       _, atm_moist = RadConvEqm(dirs, time, atm, standalone=True, cp_dry=False, trppD=False, rscatter=False)
 
       OLR_arr.append(atm_moist.LW_flux_up[0])
-      print(" ")
+      log.info(" ")
 
     OLR_arr = np.array(OLR_arr)
     Ts_arr  = np.array(Ts_arr)
@@ -110,7 +113,7 @@
                           dtype=float, skiprows=2, delimiter=',').T
 
     # Setup plot
-    print("Making plot")
+    log.info("Making plot")
     fig,ax = plt.subplots(1,1, figsize=(7,4))
 
     # Plot data
@@ -134,11 +137,12 @@
 
     fig.savefig(dirs["output"]+"runaway_demo.pdf", bbox_inches='tight')
     fig.savefig(dirs["output"]+"runaway_demo.png", bbox_inches='tight', dpi=190)
-    print(" ")
+    log.info(" ")
 
     # Tidy
     CleanOutputDir(os.getcwd())
     CleanOutputDir(dirs['output'])
 
     # Done
-    print("Done!")
+    log.info("Done!")
+
diff --git a/pyproject.toml b/pyproject.toml
index 0e5432d..6c556a8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -42,7 +42,6 @@ dependencies = [
   'tomlkit',
   'f90nml',
   'fwl-mors',
-  'datetime'
 ]
 
 [project.urls]
diff --git a/src/janus/modules/compute_moist_adiabat.py b/src/janus/modules/compute_moist_adiabat.py
index db39119..9442490 100644
--- a/src/janus/modules/compute_moist_adiabat.py
+++ b/src/janus/modules/compute_moist_adiabat.py
@@ -15,6 +15,9 @@
 from janus.modules.water_cloud import simple_cloud
 from janus.modules.relative_humidity import compute_Rh
 
+import logging
+log = logging.getLogger("fwl."+__name__)
+
 import janus.utils.GeneralAdiabat as ga # Moist adiabat with multiple condensibles
 import janus.utils.socrates as socrates
 
@@ -50,7 +53,7 @@ def compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter=False):
     atm_moist = socrates.radCompSoc(atm_moist, dirs, recalc=False, rscatter=rscatter)
 
     if standalone == True:
-        print("w/o stratosphere (net, OLR): " + str(round(atm_moist.net_flux[0], 3)) +" , "+str(round(atm_moist.LW_flux_up[0], 3)) + " W/m^2")
+        log.info("w/o stratosphere (net, OLR): %.3f, %.3f W/m^2" % ( atm_moist.net_flux[0] , atm_moist.LW_flux_up[0]))
 
     # Calculate tropopause
     if (trppD == True) or (atm_moist.trppT > atm_moist.minT):
@@ -67,6 +70,6 @@ def compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter=False):
         atm_moist = socrates.radCompSoc(atm_moist, dirs, recalc=True, rscatter=rscatter)
 
         if standalone == True:
-            print("w/ stratosphere (net, OLR): " + str(round(atm_moist.net_flux[0], 3)) + " , " + str(round(atm_moist.LW_flux_up[0], 3)) + " W/m^2")
+            log.info("w/ stratosphere (net, OLR): %.3f, %.3f W/m^2" % (atm_moist.net_flux[0] , atm_moist.LW_flux_up[0] )) 
 
     return atm_moist
diff --git a/src/janus/modules/find_tropopause.py b/src/janus/modules/find_tropopause.py
index 6bcb40e..6c09eb3 100644
--- a/src/janus/modules/find_tropopause.py
+++ b/src/janus/modules/find_tropopause.py
@@ -11,6 +11,10 @@
 
 import numpy as np
 
+import logging 
+log = logging.getLogger("fwl."+__name__)
+
+
 def find_tropopause(atm_moist, dynamic: bool, verbose=True):
     """Computes tropopause location via two methods: dynamically based on heating rates, or by a set temperature value.
 
@@ -26,8 +30,7 @@ def find_tropopause(atm_moist, dynamic: bool, verbose=True):
     # Heating criterion
     if dynamic:
 
-        if verbose:
-            print("TROPOPAUSE SET BY HEATING")
+        log.debug("TROPOPAUSE SET BY HEATING")
 
         # Find tropopause index
         trpp_idx   = 0 
@@ -44,7 +47,6 @@ def find_tropopause(atm_moist, dynamic: bool, verbose=True):
         if np.size(signchange_indices) > 0 and np.max(atm_moist.net_heating) > DeltaT_max_sign:
 
             # First guess: maximum heating index
-            # print(np.argmax(atm_moist.net_heating))
             max_heat_idx = np.argmax(atm_moist.net_heating)
 
             # Lower height while heating still significant
@@ -84,9 +86,6 @@ def find_tropopause(atm_moist, dynamic: bool, verbose=True):
         # If significant tropopause found or isothermal atmosphere from stellar heating
         if trpp_idx != 0:
 
-            # # Print tropopause index for debugging
-            # print("Tropopause @ (index, P/Pa, T/K):", trpp_idx, round(atm_moist.pl[trpp_idx],3), round(atm_moist.tmpl[trpp_idx],3))
-        
             atm_moist.trppidx   = trpp_idx                  # index
             atm_moist.trppP     = atm_moist.pl[trpp_idx]    # pressure 
             atm_moist.trppT     = atm_moist.tmpl[trpp_idx]  # temperature
@@ -94,8 +93,7 @@ def find_tropopause(atm_moist, dynamic: bool, verbose=True):
     # Temperature criterion
     else:
 
-        if verbose:
-            print("TROPOPAUSE SET BY CONTANT VALUE OF", atm_moist.trppT, "K")
+        log.debug("TROPOPAUSE SET BY CONTANT VALUE OF %.2f K"%atm_moist.trppT)
 
         # Tropopause is triggered
         if np.any(atm_moist.tmpl <= atm_moist.trppT) or np.any(atm_moist.tmp <= atm_moist.trppT):
diff --git a/src/janus/modules/relative_humidity.py b/src/janus/modules/relative_humidity.py
index 692c87b..cc1bd7c 100644
--- a/src/janus/modules/relative_humidity.py
+++ b/src/janus/modules/relative_humidity.py
@@ -41,23 +41,11 @@ def update_for_constant_RH(atm, R_h):
     """
 
     p_sat = [ga.p_sat('H2O', atm.tmp[i]) for i in range(len(atm.p))]
-    #print("update_for_constant_RH: p_sat = ", p_sat)
 
     i_lcl, lcl = find_intersection(p_sat, atm.p_vol['H2O'], tolerance=1e0)
-    #print("update_for_constant_RH: atm.p[i_lcl:] = ", atm.p[i_lcl:])
-    #print("update_for_constant_RH: atm.p_vol['H2O'][i_lcl:] = ", atm.p_vol['H2O'][i_lcl:])
-    #print("update_for_constant_RH: R_h[i_lcl:] = ", R_h[i_lcl:])
-    #print("update_for_constant_RH: atm.ps = ", atm.ps)
-    #print("update_for_constant_RH: atm.x_gas['H2O'][i_lcl:] = ", atm.x_gas['H2O'][i_lcl:])
 
-    #atm.p_vol['H2O'][i_lcl:] = R_h[i_lcl:] * p_sat[i_lcl:]  # Updating the partial pressure of H2O from the LCL down
     atm.p_vol['H2O'] = R_h * p_sat
-    #print("update_for_constant_RH: updated atm.p_vol['H2O'][i_lcl:] = ", atm.p_vol['H2O'][i_lcl:])
     atm.p = sum(atm.p_vol[molecule] for molecule in atm.vol_list.keys()) # Updating the total pressure
-    #print("update_for_constant_RH: updated atm.p[i_lcl:] = ", atm.p[i_lcl:])
     atm.ps = float(sum(atm.p_vol[molecule][-1] for molecule in atm.vol_list.keys())) # Not forgetting the surface pressure
-    #print("update_for_constant_RH: updated atm.ps = ", atm.ps)
-    #atm.x_gas['H2O'][i_lcl:] = R_h[i_lcl:] * p_sat[i_lcl:] / atm.p[i_lcl:] # Updating the mixing ratio of H2O
-    #print("atm.p_vol['H2O'][-1], atm.x_gas['H2O'], p_sat[-1] = ", atm.p_vol['H2O'][-1], atm.x_gas['H2O'][-1], p_sat[-1])
 
     return atm
diff --git a/src/janus/modules/solve_pt.py b/src/janus/modules/solve_pt.py
index 7a94c94..6d71478 100644
--- a/src/janus/modules/solve_pt.py
+++ b/src/janus/modules/solve_pt.py
@@ -10,6 +10,10 @@
 """
 
 import scipy.optimize as optimise
+
+import logging 
+log = logging.getLogger("fwl."+__name__)
+
 from janus.modules.compute_moist_adiabat import compute_moist_adiabat
 from janus.modules.dry_adiabat_timestep import compute_dry_adiabat
 from janus.utils.atmosphere_column import atmos
@@ -60,9 +64,9 @@ def RadConvEqm(dirs, time, atm, standalone:bool, cp_dry:bool, trppD:bool, rscatt
         atm_dry   = compute_dry_adiabat(atm, dirs, standalone, rscatter, pure_steam_adj, surf_dt, cp_surf, mix_coeff_atmos, mix_coeff_surf)
 
         if standalone == True:
-            print("Net, OLR => moist:", str(round(atm_moist.net_flux[0], 3)), str(round(atm_moist.LW_flux_up[0], 3)) + " W/m^2", end=" ")
-            print("| dry:", str(round(atm_dry.net_flux[0], 3)), str(round(atm_dry.LW_flux_up[0], 3)) + " W/m^2", end=" ")
-            print()
+            log.info("Net, OLR => moist: %.3f, %.3f W/m^2" % ( atm_moist.net_flux[0],atm_moist.LW_flux_up[0]))
+            log.info("^^       =>   dry: %.3f, %.3f W/m^2" % ( atm_dry.net_flux[0],  atm_dry.LW_flux_up[0]))
+            log.info("")
     else: 
         atm_dry = {}
     
@@ -169,7 +173,7 @@ def ini_atm(Ts):
     # We want to optimise this function (returns residual of F_atm and F_skn, given T_surf)
     def func(x):
 
-        print("Evaluating at T_surf = %.1f K" % x)
+        log.info("Evaluating at T_surf = %.1f K" % x)
         atm_tmp = compute_moist_adiabat(ini_atm(x), dirs, False, trppD, rscatter)
 
         if atm_bc == 0:
@@ -178,12 +182,12 @@ def func(x):
             F_atm = atm_tmp.net_flux[-1]  
 
         F_skn = skin(atm_tmp)
-        print("    F_atm = %+.2e W m-2      F_skn = %+.2e W m-2" % (F_atm,F_skn))
+        log.info("    F_atm = %+.2e W m-2      F_skn = %+.2e W m-2" % (F_atm,F_skn))
 
         del atm_tmp
         return float(F_skn - F_atm)
     
-    print("Solving for global energy balance with conductive lid (T_magma = %.1f K)" % attrs["tmp_magma"])
+    log.info("Solving for global energy balance with conductive lid (T_magma = %.1f K)" % attrs["tmp_magma"])
 
     # Use an 'initial guess' method
     if method == 0:
@@ -209,9 +213,9 @@ def func(x):
     succ  = bool(r.converged)
 
     if not succ:
-        print("WARNING: Did not find solution for surface skin balance")
+        log.warning("Did not find solution for surface skin balance")
     else:
-        print("Found surface solution")
+        log.info("Found surface solution")
         
     # Check bounds on T_surf
     T_surf = max(T_surf_sol, minT)
@@ -220,9 +224,9 @@ def func(x):
         T_surf = min(T_surf, T_surf_max)
 
     if T_surf != T_surf_sol:
-        print("T_surf limits activated")
-        print("    Found T_surf = %g K" % T_surf_sol)
-        print("    Using T_surf = %g K" % T_surf)
+        log.info("T_surf limits activated")
+        log.info("    Found T_surf = %g K" % T_surf_sol)
+        log.info("    Using T_surf = %g K" % T_surf)
 
     # Get atmosphere state from solution value
     atm = compute_moist_adiabat(ini_atm(T_surf), dirs, False, trppD, rscatter)
@@ -235,10 +239,10 @@ def func(x):
 
     F_olr = atm.LW_flux_up[0]
 
-    print("    T_surf = %g K"       % T_surf)
-    print("    F_atm  = %.4e W m-2" % F_atm)
-    print("    F_skn  = %.4e W m-2" % skin(atm))
-    print("    F_olr  = %.4e W m-2" % F_olr)
+    log.info("    T_surf = %g K"       % T_surf)
+    log.info("    F_atm  = %.4e W m-2" % F_atm)
+    log.info("    F_skn  = %.4e W m-2" % skin(atm))
+    log.info("    F_olr  = %.4e W m-2" % F_olr)
 
     return atm
 
diff --git a/src/janus/set_socrates_env.py b/src/janus/set_socrates_env.py
index 08a59db..86abdec 100644
--- a/src/janus/set_socrates_env.py
+++ b/src/janus/set_socrates_env.py
@@ -18,7 +18,7 @@
 import platformdirs
 import requests
 
-logger = logging.getLogger(__name__)
+log = logging.getLogger("fwl."+__name__)
 
 if not SOCRATES_DIR.exists():
     raise RuntimeError(f'Cannot find SOCRATES in this location: {SOCRATES_DIR}')
@@ -26,8 +26,8 @@
 with open(SOCRATES_DIR / 'version') as f:
     SOCRATES_VERSION = f.readline()
 
-logger.info(f'socrates location: %s', SOCRATES_DIR)
-logger.info('socrates version: %s', SOCRATES_VERSION)
+log.debug('socrates location: %s' % SOCRATES_DIR)
+log.debug('socrates version: %s' % SOCRATES_VERSION)
 
 sep = os.pathsep
 
diff --git a/src/janus/socrates.py b/src/janus/socrates.py
index 3471662..fa3bad2 100644
--- a/src/janus/socrates.py
+++ b/src/janus/socrates.py
@@ -10,7 +10,7 @@
 import requests
 
 SOCRATES_DATA_DIR = Path(platformdirs.user_data_dir('socrates'))
-SOCRATES_DIR = Path(os.environ.get('SOCRATES', SOCRATES_DATA_DIR / 'SOCRATES'))
+SOCRATES_DIR = Path(os.environ.get('RAD_DIR', SOCRATES_DATA_DIR / 'SOCRATES'))
 
 
 def _set_permissions(drc: Path):
@@ -64,4 +64,4 @@ def download_socrates(ref: str = 'main'):
     symlink.unlink(missing_ok=True)
     symlink.symlink_to(target_dir)
 
-    print(f'SOCRATES downloaded to {target_dir}')
+    click.echo(f'SOCRATES downloaded to {target_dir}')
diff --git a/src/janus/utils/GeneralAdiabat.py b/src/janus/utils/GeneralAdiabat.py
index 11fbec9..063d0f3 100644
--- a/src/janus/utils/GeneralAdiabat.py
+++ b/src/janus/utils/GeneralAdiabat.py
@@ -21,6 +21,9 @@
 import numpy as np
 import janus.utils.water_tables as wt
 
+import logging
+log = logging.getLogger("fwl."+__name__)
+
 from janus.utils.cp_funcs import *
 from janus.utils.ClimateUtilities import *
 from janus.utils.atmosphere_column import atmos
@@ -161,8 +164,6 @@ def atm_z(atm, idx):
     # # Use gas phase mean molar mass
     # atm.mu[idx] = atm.mu[idx]
 
-    # print(atm.grav_z[idx], atm.mu[idx], atm.p[idx], atm.p[idx+1])
-    
 
     mu_c = 0.
     if atm.xc[idx] > 0:
@@ -180,18 +181,15 @@ def atm_z(atm, idx):
     dz =  -dp / atm.rho[idx]
     atm.zl[idx+1] = atm.zl[idx] + dz
     if dz<0:
-        print("WARNING: dz  = %g < 0 " % dz)
-        print("         (dp  = %g)" % dp)
-        print("         (rho = %g)" % atm.rho[idx])
+        log.warning("dz  = %g < 0 " % dz)
+        log.warning("(dp  = %g)" % dp)
+        log.warning("(rho = %g)" % atm.rho[idx])
     
     # Next gravity
     atm.grav_z[idx+1] = atm.grav_s * ((atm.planet_radius)**2) / ((atm.planet_radius+atm.z[idx+1])**2)
 
     # Calculate scale height in km
     # H = 1e-3 * phys.R_gas * atm.tmp[idx] / (atm.mu[idx] * atm.grav_z[idx])
-    # print(H)
-
-    # print(idx, T_mean_down, dz, atm.z[idx], atm.z[idx+1], atm.grav_z[idx], atm.grav_z[idx+1])
 
     return atm
 
@@ -462,12 +460,12 @@ def slopeRay( logpa, logT ):
     pa      = math.exp(logpa)
     T       = math.exp(logT)
     qsat    = eps*(satvph2o(T)/pa)
-    print('qsat=%.3f'%qsat)
+    log.debug('qsat=%.3f'%qsat)
     num     = (1. + (L/(Ra*T))*qsat)*Ra
-    print('numerator=%.3f'%num)
+    log.debug('numerator=%.3f'%num)
     den     = cpa + (cpc + (L/(Rc*T) - 1.)*(L/T))*qsat
-    print('denominator=%.3f'%den)
-    print('dlnT/dlnPa=%.3f'%(num/den))
+    log.debug('denominator=%.3f'%den)
+    log.debug('dlnT/dlnPa=%.3f'%(num/den))
     return num/den
 
 
@@ -525,7 +523,6 @@ def dlnT_dlnP_d(lnP, lnT, atm):
         # Coefficients
         eta_vol     = atm.x_gas[vol][idx] / atm.xd[idx]
         eta_cond    = atm.x_cond[vol][idx] / atm.xd[idx]
-        #print(eta_vol)
         # sums for saturated comps
         # HII added, ensure below the critical point !
         if ((np.isclose(atm.p_vol[vol][idx] ,p_sat(vol,tmp, water_lookup=atm.water_lookup)) or
@@ -544,7 +541,6 @@ def dlnT_dlnP_d(lnP, lnT, atm):
             
         # sums for subsaturated comps
         else: 
-            #print(eta_vol)
             denom_sum += cpv(vol,tmp) * eta_vol #The eta_vol (x_vol/x_d) is there so that if there are multiple dry components the average dry specific heat is used
         
     
@@ -581,7 +577,6 @@ def moist_slope(lnP, lnT, atm):
         # Coefficients
         eta_vol     = atm.x_gas[vol][idx] / atm.xd[idx]
         
-        #print(eta_vol)
         # sums for saturated comps
         # Ensure below the critical point for volatile before we
         # assume condensation occurs
@@ -651,7 +646,6 @@ def condensation( atm, idx, wet_list, dry_list, prs_reset):
         # Calculate the individual partial pressures of dry species
         for vol in atm.vol_list:
             if vol in dry_list:
-                #print(atm.x_gas[vol][idx-1])
                 atm.p_vol[vol][idx] = p_dry_tot * (atm.x_gas[vol][idx-1]/dry_frac_sum)
                 
        
@@ -837,8 +831,6 @@ def general_adiabat( atm ):
         # Calculate next local gravity and height
         atm = atm_z(atm, idx)
 
-        # print("RK4 step, idx:", idx, round(atm.p[idx+1],5), round(atm.tmp[idx+1],5))
-
         # Set next level to calculate
         idx             += 1
         atm.ifatm[idx]  = idx
diff --git a/src/janus/utils/RayleighSpectrum.py b/src/janus/utils/RayleighSpectrum.py
index 1370ad2..f15be03 100755
--- a/src/janus/utils/RayleighSpectrum.py
+++ b/src/janus/utils/RayleighSpectrum.py
@@ -62,8 +62,6 @@ def cross_section(wavelength, info_dict):#wavelength in m
 def band_integrator(species_list, molar_mixing_ratio_list, wavelength1_list, wavelength2_list):
     total_coefficient = 0
     for n, species in enumerate(species_list):
-        #print(species)
-        #print(species_info(species))
         info_tup = species_info(species)
         
         #Coefficient w/o mixing ratio weight
@@ -146,7 +144,6 @@ def rayleigh_coeff_adder(species_list = ['co2'], mixing_ratio_list = [1.], spect
     
     # This loop puts the cross sections in the proper format and then adds them to the list line by line
     for n in range(cross_section_list.shape[0]):
-        #print('{0:5d}'.format(np.array(cross_section_list[n,0],dtype=int))+'\t\t %.9E\n'%(cross_section_list[n,1]))
         block3_list.append('{0:5d}'.format(np.array(cross_section_list[n,0],dtype=int))+'        %.9E\n'%(cross_section_list[n,1]))
     
     # Adding the end statement of block 3
diff --git a/src/janus/utils/StellarSpectrum.py b/src/janus/utils/StellarSpectrum.py
index 7404f98..15b6cf2 100644
--- a/src/janus/utils/StellarSpectrum.py
+++ b/src/janus/utils/StellarSpectrum.py
@@ -1,18 +1,13 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Jun 15 16:07:20 2023
-
-@author: nichollsh
-"""
-
 import numpy as np
 import shutil , os
 import subprocess
 from scipy.interpolate import PchipInterpolator
 
-from .. import set_socrates_env  # noqa
+import logging 
+log = logging.getLogger("fwl."+__name__)
+
 
+from .. import set_socrates_env 
 
 def PrepareStellarSpectrum(wl, fl, star_file, nbins_max=95000):
     """Write a stellar spectrum.
@@ -40,7 +35,7 @@ def PrepareStellarSpectrum(wl, fl, star_file, nbins_max=95000):
     if (len(wl) != len(fl)):
         raise Exception("Stellar wavelength and flux arrays have different lengths")
     if (len(wl) < 500):
-        print("WARNING: Loaded stellar spectrum is very short!")
+        log.warning("Loaded stellar spectrum is very short!")
 
     nbins_max = min(len(wl), nbins_max)
     
@@ -49,14 +44,14 @@ def PrepareStellarSpectrum(wl, fl, star_file, nbins_max=95000):
 
     # Down-sample spectrum when necessary or requested
     if len(wl) > socrates_nbins_max:
-        print("Rebinning stellar spectrum")
+        log.info("Rebinning stellar spectrum")
 
         # Store old wl,fl arrays
         wl_orig = wl
         fl_orig = fl
 
         if (nbins_max < 500):
-            print("WARNING: Requested number of bins is small (%d bins)" % nbins_max)
+            log.warning("Requested number of bins is small (%d bins)" % nbins_max)
 
         nbins_max = min( int(socrates_nbins_max), nbins_max) # Must be fewer than 100k
 
@@ -134,6 +129,6 @@ def InsertStellarSpectrum(orig_file:str, star_file:str, output_folder:str):
               ]
     p = subprocess.run(['prep_spec'], stdout=subprocess.PIPE, input='\n'.join(inputs), encoding='ascii')
     if (p.returncode != 0):
-        print("WARNING: prep_spec returned with code %d" % p.returncode)
+        log.warning("prep_spec returned with code %d" % p.returncode)
     
 
diff --git a/src/janus/utils/cp_funcs.py b/src/janus/utils/cp_funcs.py
index 2750c6d..40049e8 100644
--- a/src/janus/utils/cp_funcs.py
+++ b/src/janus/utils/cp_funcs.py
@@ -229,7 +229,6 @@ def cp_cond( vol, tmp, cp_mode='T-dependent'):
 
     match vol:
     
-        #print(vol)
         # https://webbook.nist.gov/cgi/fluid.cgi?TLow=274&THigh=647&TInc=20&Applet=on&Digits=5&ID=C7732185&Action=Load&Type=SatP&TUnit=K&PUnit=MPa&DUnit=mol%2Fl&HUnit=kJ%2Fmol&WUnit=m%2Fs&VisUnit=uPa*s&STUnit=N%2Fm&RefState=DEF
         case "H2O":
             
diff --git a/src/janus/utils/data.py b/src/janus/utils/data.py
index fb7b73a..7444009 100644
--- a/src/janus/utils/data.py
+++ b/src/janus/utils/data.py
@@ -5,11 +5,11 @@
 import platformdirs
 from osfclient.api import OSF
 
-logger = logging.getLogger(__name__)
+log = logging.getLogger("fwl."+__name__)
 
 FWL_DATA_DIR = Path(os.environ.get('FWL_DATA', platformdirs.user_data_dir('fwl_data')))
 
-logger.info(f'FWL data location: {FWL_DATA_DIR}')
+log.info(f'FWL data location: {FWL_DATA_DIR}')
 
 basic_list = (
         "Dayspring/256",
@@ -36,7 +36,7 @@ def download_folder(*, storage, folders: list[str], data_dir: Path):
             parts = file.path.split('/')[1:]
             target = Path(data_dir, *parts)
             target.parent.mkdir(parents=True, exist_ok=True)
-            print(f'Downloading {file.path}...')
+            log.info(f'Downloading {file.path}...')
             with open(target, 'wb') as f:
                 file.write_to(f)
             break
diff --git a/src/janus/utils/find_lcl.py b/src/janus/utils/find_lcl.py
index 5c50a00..465d608 100644
--- a/src/janus/utils/find_lcl.py
+++ b/src/janus/utils/find_lcl.py
@@ -1,4 +1,7 @@
 import numpy as np
+import logging
+
+log = logging.getLogger("fwl."+__name__)
 
 def find_intersection(curve1, curve2, tolerance=1e-1):
     # Starting from the TOA, find the last index where the absolute difference is within the tolerance
@@ -8,9 +11,7 @@ def find_intersection(curve1, curve2, tolerance=1e-1):
         first_intersection_index = intersection_indices[-1]
         first_intersection_value = curve1[first_intersection_index]
         
-        #print("LCL index:", first_intersection_index)
-        #print("LCL value:", first_intersection_value)
         return first_intersection_index, first_intersection_value
     else:
-        print("No LCL found")
+        log.error("No LCL found")
         return None
\ No newline at end of file
diff --git a/src/janus/utils/height.py b/src/janus/utils/height.py
index 2d0d1cd..47cc339 100644
--- a/src/janus/utils/height.py
+++ b/src/janus/utils/height.py
@@ -1,6 +1,10 @@
 import numpy as np
 import janus.utils.phys as phys
 
+import logging 
+log = logging.getLogger("fwl."+__name__)
+
+
 def gravity( m, r ):
     g = phys.G*m/r**2
     return g
@@ -37,7 +41,7 @@ def AtmosphericHeight(atm, m_planet, r_planet):
         # This implies that the hydrostatic/gravity integration failed.
         if z_profile[n+1] > 1.0e8:
             atm.height_error = True 
-            print("WARNING: Hydrostatic integration blew up. Setting dummy values for height")
+            log.warning("Hydrostatic integration blew up. Setting dummy values for height")
             break
 
     # Set dummy values 
diff --git a/src/janus/utils/logs.py b/src/janus/utils/logs.py
new file mode 100644
index 0000000..b8797cc
--- /dev/null
+++ b/src/janus/utils/logs.py
@@ -0,0 +1,31 @@
+import logging
+import sys
+
+# Simple terminal-logger instance 
+def setup_logger(level:str="INFO"):
+
+    custom_logger = logging.getLogger("fwl."+__name__)
+    custom_logger.handlers.clear()
+
+    level = str(level).strip().upper()
+    if level not in ["INFO", "DEBUG", "ERROR", "WARNING"]:
+        raise ValueError("Invalid log level '%s'"%level)
+    level_code = logging.getLevelName(level)
+
+    # Add terminal output to logger
+    sh = logging.StreamHandler(sys.stdout)
+    sh.setLevel(level_code)
+    custom_logger.addHandler(sh)
+    custom_logger.setLevel(level_code)
+
+    # Capture unhandled exceptions
+    # https://stackoverflow.com/a/16993115
+    def handle_exception(exc_type, exc_value, exc_traceback):
+        if issubclass(exc_type, KeyboardInterrupt):
+            custom_logger.error("KeyboardInterrupt")
+            sys.__excepthook__(exc_type, exc_value, exc_traceback)
+            return
+        custom_logger.critical("Uncaught exception", exc_info=(exc_type, exc_value, exc_traceback))
+    sys.excepthook = handle_exception
+    
+    return custom_logger
diff --git a/src/janus/utils/nctools.py b/src/janus/utils/nctools.py
index 5647fd0..bb17f79 100644
--- a/src/janus/utils/nctools.py
+++ b/src/janus/utils/nctools.py
@@ -47,8 +47,6 @@ def ncout_surf(file, lon, lat, basis, alb):
         raise RuntimeError(' Error in ncout_surf: arrays dont match', nvals
                             , n_lon * n_lat * levels)
 
-    #print('ncout_surf - file: ',file)
-
     ncdf_file = create_cdf(file)
 
     write_dim(ncdf_file, n_lon, lon, 'lon', 'f4', 'lon', 'degree', 'LONGITUDE')
@@ -83,8 +81,6 @@ def ncout_spectral_surf(file, lon, lat, bands, alb):
         raise RuntimeError(' Error in ncout_surf: arrays dont match', nvals
                            , n_lon * n_lat * bands)
 
-    #print('ncout_spectral_surf - file: ', file)
-
     ncdf_file = create_cdf(file)
     write_dim(ncdf_file, n_lon, lon, 'lon', 'f4', 'lon', 'degree', 'LONGITUDE')
     write_dim(ncdf_file, n_lat, lat, 'lat', 'f4', 'lat', 'degree', 'LATITUDE')
@@ -123,7 +119,6 @@ def ncout2d(file, lon, lat, val, name = None, longname = None
     if (type(name) is not str):
         patterns = re.compile(r'\\|\.')
         name = file[patterns.search(file).start() + 1:]
-    #print('ncout2d - file: ',file)
 
     ncdf_file = create_cdf(file)
 
@@ -173,8 +168,6 @@ def ncout3d(file, lon, lat, p, val, name = None
         patterns = re.compile(r'\\|\.')
         name = file[patterns.search(file).start() + 1:]
 
-    #print('ncout3d - file: ',file)
-
     ncdf_file = create_cdf(file)
 
     write_dim(ncdf_file, n_lon, lon, 'lon', 'f4', 'lon', 'degree', 'LONGITUDE')
@@ -256,17 +249,12 @@ def ncout_opt_prop(file, lon, lat, p, bands, absp, scat, phf):
         scat_vals = scat_vals[:, order, :, :]
         phf_vals = phf_vals[:, :, order, :, :]
 
-    #print('ncout_opt_prop - file: ', file)
-
     ncdf_file = create_cdf(file)
 
     write_dim(ncdf_file, n_lon, lon, 'lon', 'f4', 'lon', 'degree', 'LONGITUDE')
     write_dim(ncdf_file, n_lat, lat, 'lat', 'f4', 'lat', 'degree', 'LATITUDE')
     write_dim(ncdf_file, levels, p, 'plev', 'f4', 'plev', 'Pa', 'PRESSURE')
     write_dim(ncdf_file, 1, 1, 'mom', 'i2', 'mom', 'none', 'moment')
-    #print('bands=', bands)
-    #print('np.sum(bands)=', np.sum(bands))
-    #print('np.arange(bands)+1=', np.arange(bands)+1)
     write_dim(ncdf_file, bands, np.arange(bands)+1, 'band', 'i2', 'band',
               'none', 'band')
 
diff --git a/src/janus/utils/socrates.py b/src/janus/utils/socrates.py
index c81ef24..997b825 100644
--- a/src/janus/utils/socrates.py
+++ b/src/janus/utils/socrates.py
@@ -12,6 +12,9 @@
 import subprocess
 import f90nml
 
+import logging 
+log = logging.getLogger("fwl."+__name__)
+
 import janus.utils.nctools as nctools
 import janus.utils.RayleighSpectrum as RayleighSpectrum
 from janus.utils.atmosphere_column import atmos
@@ -60,14 +63,14 @@ def radCompSoc(atm, dirs, recalc, rscatter=False,
 
     # Check that atmosphere is okay
     if np.any(atm.cp <= 0):
-        print("ERROR: Negative heat capacity!")
+        log.error("Negative heat capacity!")
         exit(1)
     vals = []
     for arr in atm.x_gas.values():
         vals.append(arr)
     if np.any(np.array(vals).flatten() < 0):
-        print("ERROR: Negative mixing ratio(s)!")
-        print(atm.x_gas)
+        log.error("Negative mixing ratio(s)!")
+        log.error(str(atm.x_gas))
         exit(1)
     
     # Rayleigh scattering
@@ -329,9 +332,3 @@ def CleanOutputDir( output_dir ):
     for file in natural_sort(files_to_delete):
         os.remove(file)
 
-
-# Disable and enable print: https://stackoverflow.com/questions/8391411/suppress-calls-to-print-python
-def blockPrint():
-    sys.stdout = open(os.devnull, 'w')
-def enablePrint():
-    sys.stdout = sys.__stdout__
diff --git a/tests/helpers/__init__.py b/tests/helpers/__init__.py
index d744047..4dbd680 100644
--- a/tests/helpers/__init__.py
+++ b/tests/helpers/__init__.py
@@ -67,7 +67,7 @@ def get_atmosphere_config(*, band_edges, cfg_name: str, distance: float):
 
     atm = atmos.from_file(cfg_file, band_edges, vol_mixing=vol_mixing, vol_partial={})
 
-    mors.DownloadEvolutionTracks('/Baraffe')
+    mors.DownloadEvolutionTracks('Baraffe')
     baraffe = mors.BaraffeTrack(star_mass)
     atm.instellation = baraffe.BaraffeSolarConstant(star_time, distance)
 
diff --git a/tools/plotting_tools/plot_cff.py b/tools/plotting_tools/plot_cff.py
index 435c54e..47be5c3 100644
--- a/tools/plotting_tools/plot_cff.py
+++ b/tools/plotting_tools/plot_cff.py
@@ -76,8 +76,6 @@
 atm = pkl.load(atm_file_stream)
 atm_file_stream.close()
 
-# print(atm.p)
-
 ##### PLOT A
 
 # # Loop through volatiles, options: "H2O", "CO2", "H2", "N2", "CH4", "CO", "O2"