From 1655b6d37cc22cc0f136931f9459496b7f2d4c0f Mon Sep 17 00:00:00 2001 From: Harrison Nicholls Date: Wed, 11 Sep 2024 17:03:38 +0100 Subject: [PATCH 1/6] Routine to calculate observed bulk density --- src/janus/utils/observed_rho.py | 43 +++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 src/janus/utils/observed_rho.py diff --git a/src/janus/utils/observed_rho.py b/src/janus/utils/observed_rho.py new file mode 100644 index 0000000..8ea4583 --- /dev/null +++ b/src/janus/utils/observed_rho.py @@ -0,0 +1,43 @@ +import numpy as np +import logging + +log = logging.getLogger("fwl."+__name__) + +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 * atm.transspec_m / (4.0 * np.pi * transspec_r**3) + + return transspec_rho + From e37c965b9315c05433ec12b68375791293f37074 Mon Sep 17 00:00:00 2001 From: Harrison Nicholls Date: Wed, 11 Sep 2024 17:09:37 +0100 Subject: [PATCH 2/6] Updated gitignore and tidied imports --- .gitignore | 2 ++ src/janus/utils/observed_rho.py | 3 --- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index b3575d8..3185e03 100644 --- a/.gitignore +++ b/.gitignore @@ -34,4 +34,6 @@ build/ fwl_data*/ dist/ .pytest_cache +.ruff_cache + diff --git a/src/janus/utils/observed_rho.py b/src/janus/utils/observed_rho.py index 8ea4583..709603c 100644 --- a/src/janus/utils/observed_rho.py +++ b/src/janus/utils/observed_rho.py @@ -1,7 +1,4 @@ import numpy as np -import logging - -log = logging.getLogger("fwl."+__name__) def calc_observed_rho(atm): """Calculate the observed bulk density. From 582e96ae34ba6c3adfc34e58761dda391ecf8218 Mon Sep 17 00:00:00 2001 From: Harrison Nicholls Date: Wed, 11 Sep 2024 17:15:25 +0100 Subject: [PATCH 3/6] Fixed typo --- src/janus/utils/observed_rho.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/janus/utils/observed_rho.py b/src/janus/utils/observed_rho.py index 709603c..5bfec19 100644 --- a/src/janus/utils/observed_rho.py +++ b/src/janus/utils/observed_rho.py @@ -34,7 +34,7 @@ def calc_observed_rho(atm): transspec_m += atm.planet_mass # get density of all enclosed by observed layer - transspec_rho = 3.0 * atm.transspec_m / (4.0 * np.pi * transspec_r**3) + transspec_rho = 3.0 * transspec_m / (4.0 * np.pi * transspec_r**3) return transspec_rho From d473120924c7755c718cde30b9dba6d7b9d83bb7 Mon Sep 17 00:00:00 2001 From: Harrison Nicholls Date: Wed, 11 Sep 2024 17:22:35 +0100 Subject: [PATCH 4/6] Updated README --- README.md | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 405d63c..45f51a7 100644 --- a/README.md +++ b/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 (tim.lichtenberg@rug.nl) -* MH - Mark Hammond (mark.hammond@physics.ox.ac.uk) -* RB – Ryan Boukrouche (ryan.boukrouche@astro.su.se) -* RJG – RJ Graham (arejaygraham@uchicago.edu) -* HN - Harrison Nicholls (harrison.nicholls@physics.ox.ac.uk) -* HII - Hamish Innes (hamish.innes@physics.ox.ac.uk) -* LS - Laurent Soucasse (l.soucasse@esciencecenter.nl) +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 From 0e2071869f121f54c4ab668cbf135b5241dd342a Mon Sep 17 00:00:00 2001 From: Harrison Nicholls Date: Wed, 11 Sep 2024 20:05:13 +0100 Subject: [PATCH 5/6] Improvements for check with hydrostatic integration fails --- src/janus/utils/atmosphere_column.py | 7 ++++--- src/janus/utils/height.py | 24 ++++++++++++------------ 2 files changed, 16 insertions(+), 15 deletions(-) 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 From f428906727fa608e62c92fd19ddf8e4f495b5ac4 Mon Sep 17 00:00:00 2001 From: Harrison Nicholls Date: Wed, 11 Sep 2024 20:16:00 +0100 Subject: [PATCH 6/6] Saves height_error to output netcdf file --- src/janus/utils/atmosphere_column.py | 13 ++++++------- src/janus/utils/height.py | 28 +++++++++++++++++----------- 2 files changed, 23 insertions(+), 18 deletions(-) diff --git a/src/janus/utils/atmosphere_column.py b/src/janus/utils/atmosphere_column.py index d852535..41bd57a 100644 --- a/src/janus/utils/atmosphere_column.py +++ b/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' # ---------------------- diff --git a/src/janus/utils/height.py b/src/janus/utils/height.py index ad99187..c87eb7c 100644 --- a/src/janus/utils/height.py +++ b/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