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/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 diff --git a/src/janus/utils/atmosphere_column.py b/src/janus/utils/atmosphere_column.py index b66a2a1..41bd57a 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,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' # ---------------------- diff --git a/src/janus/utils/height.py b/src/janus/utils/height.py index 47cc339..c87eb7c 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,9 +9,9 @@ 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,11 +20,11 @@ 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 @@ -32,27 +32,33 @@ 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 + 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 diff --git a/src/janus/utils/observed_rho.py b/src/janus/utils/observed_rho.py new file mode 100644 index 0000000..5bfec19 --- /dev/null +++ b/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 +