Saves height_error to output netcdf file

src/janus/utils/atmosphere_column.py

@@ -389,13 +389,7 @@ def write_ncdf(self, fpath:str):
 
         # ----------------------
         # Calculate gravity and height (in case it hasn't been done already)
-        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)):
-            self.zl[i] = 0.5 * (self.z[i-1] + self.z[i])
-        self.zl[0] = 2*self.z[0] - self.zl[1] # estimate TOA height
+        self.z, self.zl, self.height_error = integrate_heights(self, self.planet_mass, self.planet_radius)
 
         # ----------------------
         # Prepare NetCDF
@@ -455,6 +449,7 @@ def write_ncdf(self, fpath:str):
         var_albsurf = ds.createVariable("surf_albedo",   'f4');                            # Surface albedo
         var_sknd    = ds.createVariable("cond_skin_d"   ,'f4');  var_sknd.units = "m"      # Conductive skin thickness
         var_sknk    = ds.createVariable("cond_skin_k"   ,'f4');  var_sknk.units = "W m-1 K-1"    # Conductive skin thermal conductivity
+        var_zok     = ds.createVariable("height_ok"     ,'S1'); # Hydrostatic integration is ok
 
         #     Store data
         var_tstar.assignValue(self.ts)
@@ -472,6 +467,10 @@ def write_ncdf(self, fpath:str):
         var_albsurf.assignValue(self.albedo_s)
         var_sknd.assignValue(self.skin_d)
         var_sknk.assignValue(self.skin_k)
+        if self.height_error:
+            var_zok[0] = 'n'
+        else:
+            var_zok[0] = 'y'
 
 
         # ----------------------

src/janus/utils/height.py

@@ -11,7 +11,7 @@ def gravity( m, r ):
 
 def integrate_heights(atm, m_planet, r_planet):
 
-    z_profile       = np.zeros(len(atm.p))
+    z       = np.zeros(len(atm.p))
     grav_s          = gravity( m_planet, r_planet )
 
     # Reverse arrays to go from high to low pressure
@@ -20,11 +20,11 @@ def integrate_heights(atm, m_planet, r_planet):
     for vol in atm.vol_list.keys():
         atm.x_gas[vol] = atm.x_gas[vol][::-1]
 
-    ok = True
-    for n in range(0, len(z_profile)-1):
+    height_error = False
+    for n in range(0, len(z)-1):
 
         # Gravity with height
-        grav_z = grav_s * ((r_planet)**2) / ((r_planet + z_profile[n])**2)
+        grav_z = grav_s * ((r_planet)**2) / ((r_planet + z[n])**2)
 
         # Mean molar mass depending on mixing ratio
         mean_molar_mass = 0
@@ -35,24 +35,30 @@ def integrate_heights(atm, m_planet, r_planet):
         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
+        z[n+1] = z[n] + dz
 
         # Check if heights are very large.
         # This implies that the hydrostatic/gravity integration failed.
-        if (z_profile[n+1] > 1.0e8) or (dz > 1e8):
-            ok = False
+        if (z[n+1] > 1.0e8) or (dz > 1e8):
+            height_error = True
             log.error("Hydrostatic integration blew up. Setting dummy values for height")
             break
 
     # Set dummy values
-    if not ok:
-        z_profile = np.linspace(0.0, 1000.0, len(atm.p))
+    if height_error:
+        z = 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]
     for vol in atm.vol_list.keys():
         atm.x_gas[vol] = atm.x_gas[vol][::-1]
-    z_profile = z_profile[::-1]
+    z = z[::-1]
 
-    return z_profile, ok
+    # Set cell edge values
+    zl = np.zeros(len(z)+1)
+    for i in range(1,len(z)):
+        zl[i] = 0.5 * (z[i-1] + z[i])
+    zl[0] = 2*z[0] - zl[1] # estimate TOA height
+
+    return z, zl, height_error