diff --git a/src/janus/utils/atmosphere_column.py b/src/janus/utils/atmosphere_column.py
index b66a2a1..d852535 100644
--- a/src/janus/utils/atmosphere_column.py
+++ b/src/janus/utils/atmosphere_column.py
@@ -5,7 +5,7 @@
 import numpy as np
 import netCDF4 as nc
 from janus.utils import phys
-from janus.utils.height import AtmosphericHeight
+from janus.utils.height import integrate_heights
 import os, copy, platform, shutil
 import pwd
 from datetime import datetime
@@ -172,7 +172,7 @@ def __init__(self, T_surf: float, P_surf: float, P_top: float, pl_radius: float,
         self.rho            = np.zeros(self.nlev) # Density of atmosphere at a given level
         self.ifatm 			= np.zeros(self.nlev) # Defines nth level to which atmosphere is calculated
         self.cp      		= np.zeros(self.nlev) # Mean heat capacity
-        self.height_error   = True # error when doing hydrostatic integration?
+        self.height_error   = False # error when doing hydrostatic integration?
 
         # Define T and P arrays from surface up
         self.tmp[0]         = self.ts         		# K
@@ -389,7 +389,8 @@ def write_ncdf(self, fpath:str):
 
         # ----------------------
         # Calculate gravity and height (in case it hasn't been done already)
-        self.z = AtmosphericHeight(self, self.planet_mass, self.planet_radius)
+        self.z, height_err = integrate_heights(self, self.planet_mass, self.planet_radius)
+        self.height_error = height_err
 
         self.zl = np.zeros(len(self.z)+1)
         for i in range(1,len(self.z)):
diff --git a/src/janus/utils/height.py b/src/janus/utils/height.py
index 47cc339..ad99187 100644
--- a/src/janus/utils/height.py
+++ b/src/janus/utils/height.py
@@ -1,7 +1,7 @@
 import numpy as np
 import janus.utils.phys as phys
 
-import logging 
+import logging
 log = logging.getLogger("fwl."+__name__)
 
 
@@ -9,7 +9,7 @@ def gravity( m, r ):
     g = phys.G*m/r**2
     return g
 
-def AtmosphericHeight(atm, m_planet, r_planet):
+def integrate_heights(atm, m_planet, r_planet):
 
     z_profile       = np.zeros(len(atm.p))
     grav_s          = gravity( m_planet, r_planet )
@@ -20,7 +20,7 @@ def AtmosphericHeight(atm, m_planet, r_planet):
     for vol in atm.vol_list.keys():
         atm.x_gas[vol] = atm.x_gas[vol][::-1]
 
-    atm.height_error = False 
+    ok = True
     for n in range(0, len(z_profile)-1):
 
         # Gravity with height
@@ -32,22 +32,22 @@ def AtmosphericHeight(atm, m_planet, r_planet):
             mean_molar_mass += phys.molar_mass[vol]*atm.x_gas[vol][n]
 
         # Use hydrostatic equation to get height difference
-        dz = phys.R_gas * atm.tmp[n] / (mean_molar_mass * grav_z * atm.p[n]) * (atm.p[n] - atm.p[n+1]) 
-        
+        dz = phys.R_gas * atm.tmp[n] / (mean_molar_mass * grav_z * atm.p[n]) * (atm.p[n] - atm.p[n+1])
+
         # Next height
         z_profile[n+1] = z_profile[n] + dz
 
         # Check if heights are very large.
         # This implies that the hydrostatic/gravity integration failed.
-        if z_profile[n+1] > 1.0e8:
-            atm.height_error = True 
-            log.warning("Hydrostatic integration blew up. Setting dummy values for height")
+        if (z_profile[n+1] > 1.0e8) or (dz > 1e8):
+            ok = False
+            log.error("Hydrostatic integration blew up. Setting dummy values for height")
             break
 
-    # Set dummy values 
-    if atm.height_error:
+    # Set dummy values
+    if not ok:
         z_profile = np.linspace(0.0, 1000.0, len(atm.p))
-        
+
     # Reverse arrays again back to normal
     atm.p   = atm.p[::-1]
     atm.tmp = atm.tmp[::-1]
@@ -55,4 +55,4 @@ def AtmosphericHeight(atm, m_planet, r_planet):
         atm.x_gas[vol] = atm.x_gas[vol][::-1]
     z_profile = z_profile[::-1]
 
-    return z_profile
+    return z_profile, ok