Merge pull request #59 from FormingWorlds/density

Calculates observed bulk density

.gitignore

@@ -34,4 +34,6 @@ build/
 fwl_data*/
 dist/
 .pytest_cache
+.ruff_cache
+
 

README.md

@@ -3,21 +3,27 @@
 
 ## JANUS (1D convective atmosphere model)
 
-Generates a temperature profile using the generalised moist pseudoadiabat and a prescribed stratosphere. Calculates radiative fluxes using SOCRATES.   
+Generates a temperature profile using the generalised moist pseudoadiabat and a prescribed stratosphere. Calculates radiative fluxes using SOCRATES.
 
 Pronounced *jan-us*. *Jan* as in "january", and *us* as in the collective pronoun.
 
 ### Documentation
-https://proteus-code.readthedocs.io
-
-### Contributors (abbreviations & email addresses):
-* TL - Tim Lichtenberg ([email protected])
-* MH - Mark Hammond ([email protected])
-* RB – Ryan Boukrouche ([email protected])
-* RJG – RJ Graham ([email protected])
-* HN - Harrison Nicholls ([email protected])
-* HII - Hamish Innes ([email protected])
-* LS - Laurent Soucasse ([email protected])
+https://fwl-proteus.readthedocs.io
+
+## Contributors
+
+| Name  | Email address |
+| -     | -             |
+Tim Lichtenberg         | tim.lichtenberg[at]rug.nl |
+Harrison Nicholls       | harrison.nicholls[at]physics.ox.ac.uk |
+Laurent Soucasse        | l.soucasse[at]esciencecenter.nl |
+Stef Smeets             | s.smeets[at]esciencecenter.nl |
+Mark Hammond            | mark.hammond[at]physics.ox.ac.uk |
+RJ Graham               | arejaygraham[at]uchicago.edu |
+Raymond Pierrehumbert   | raymond.pierrehumbert[at]physics.ox.ac.uk |
+Ryan Boukrouche         | ryan.boukrouche[at]astro.su.se |
+Hamish Innes            | hamish.innes[at]fu-berlin.de |
+
 
 ### Repository structure
 * `README.md`           - This file

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,12 +389,7 @@ 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.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
@@ -454,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)
@@ -471,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

@@ -1,17 +1,17 @@
 import numpy as np
 import janus.utils.phys as phys
 
-import logging 
+import logging
 log = logging.getLogger("fwl."+__name__)
 
 
 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))
+    z       = np.zeros(len(atm.p))
     grav_s          = gravity( m_planet, r_planet )
 
     # Reverse arrays to go from high to low pressure
@@ -20,39 +20,45 @@ 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 
-    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
         for vol in atm.vol_list.keys():
             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
+        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:
-            atm.height_error = True 
-            log.warning("Hydrostatic integration blew up. Setting dummy values for height")
+        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 atm.height_error:
-        z_profile = np.linspace(0.0, 1000.0, len(atm.p))
-        
+    # Set dummy values
+    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]
+
+    # 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_profile
+    return z, zl, height_error

src/janus/utils/observed_rho.py

@@ -0,0 +1,40 @@
+import numpy as np
+
+def calc_observed_rho(atm):
+    """Calculate the observed bulk density.
+
+    Copied from AGNI.
+
+    Parameters
+    ----------
+        atm : atmos
+            Atmosphere object from atmosphere_column.py
+
+    Returns
+    ----------
+        rho : float
+            Observed bulk density [kg m-3]
+    """
+
+    # transspec_r::Float64            # planet radius probed in transmission [m]
+    # transspec_m::Float64            # mass [kg] of atmosphere + interior
+    # transspec_rho::Float64          # bulk density [kg m-3] implied by r and m
+
+    # arguments
+    transspec_p:float = 1e2 # (INPUT) level probed in transmission [Pa]
+
+    # get the observed height
+    idx = int(np.argmin(np.abs(atm.p - transspec_p)))
+    transspec_r = atm.z[idx] + atm.planet_radius
+
+    # get mass of whole atmosphere, assuming hydrostatic
+    transspec_m = np.amax(atm.pl) * 4 * np.pi * atm.planet_radius**2 / atm.grav_s
+
+    # add mass of the interior component
+    transspec_m += atm.planet_mass
+
+    # get density of all enclosed by observed layer
+    transspec_rho = 3.0 * transspec_m / (4.0 * np.pi * transspec_r**3)
+
+    return transspec_rho
+