Skip to content

Commit

Permalink
Introduced a new atmospher constructor based on a toml config file.
Browse files Browse the repository at this point in the history
  • Loading branch information
@lsoucasse
lsoucasse committed Jun 12, 2024
1 parent 528661d commit 045c178
Show file tree
Hide file tree
Showing 5 changed files with 151 additions and 103 deletions.
17 changes: 17 additions & 0 deletions src/janus/data/tests/config_instellation.toml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
[planet]
time = 0.0
pl_radius = 6.371e6
pl_mass = 5.972e24

[star]
time = 100e+6 # yr
star_mass = 1.0 # M_sun, mass of star

[atmos]
T_surf = 2800.0
P_surf = 2.8e7
P_top = 1.0
req_levels = 150
albedo_pl = 0.1
albedo_s = 0.1
zenith_angle = 48.19 # cronin+14 (also for scaling by a factor of 3/8 ^^)
14 changes: 14 additions & 0 deletions src/janus/data/tests/config_runaway.toml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
[planet]
time = 0.0
pl_radius = 6.371e6
pl_mass = 5.972e24

[star]
time = 100e+6 # yr
star_mass = 1.0 # M_sun, mass of star
mean_distance = 0.5 # au, orbital distance

[atmos]
P_surf = 3.0e7
P_top = 1.0
trppT = 0.0
72 changes: 65 additions & 7 deletions src/janus/utils/atmosphere_column.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# atmosphere_column.py
# class for atmospheric column data

import toml
import numpy as np
import netCDF4 as nc
from janus.utils import phys
Expand All @@ -11,10 +12,12 @@
class atmos:

def __init__(self, T_surf: float, P_surf: float, P_top: float, pl_radius: float, pl_mass: float,
band_edges:list,
vol_mixing: dict = {}, vol_partial: dict = {},
band_edges:list, vol_mixing: dict = {}, vol_partial: dict = {},
req_levels: int = 100, water_lookup: bool=False, alpha_cloud:float=0.0,
trppT: float = 290.0, minT: float = 1.0, maxT: float = 9000.0, do_cloud: bool=False, re: float=0., lwm: float=0., clfr: float=0.):
trppT: float = 290.0, minT: float = 1.0, maxT: float = 9000.0, do_cloud: bool=False,
re: float=0., lwm: float=0., clfr: float=0.,
albedo_s: float=0.0, albedo_pl: float=0.175, zenith_angle: float=54.74
):

"""Atmosphere class
Expand Down Expand Up @@ -62,7 +65,12 @@ def __init__(self, T_surf: float, P_surf: float, P_top: float, pl_radius: float,
Liquid water mass fraction [kg/kg]
clfr : float
Water cloud fraction [adimensional]
albedo_s : float
Surface albedo
albedo_pl : float
Bond albedo (scattering) applied to self.toa_heating in socrates.py
zenith_angle : float
solar zenith angle, Hamano+15 (arccos(1/sqrt(3) = 54.74), Wordsworth+10: 48.19
"""

# Parse volatiles
Expand Down Expand Up @@ -128,9 +136,9 @@ def __init__(self, T_surf: float, P_surf: float, P_top: float, pl_radius: float,
self.toa_heating = 0. # ASF

self.inst_sf = 3.0/8.0 # Scale factor applied to instellation (see Cronin+14 for definitions)
self.albedo_s = 0.0 # surface albedo
self.albedo_pl = 0.175 # Bond albedo (scattering) applied to self.toa_heating in socrates.py
self.zenith_angle = 54.74 # solar zenith angle, Hamano+15 (arccos(1/sqrt(3) = 54.74), Wordsworth+10: 48.19 (arccos(2/3)), see Cronin+14 for definitions
self.albedo_s = albedo_s # surface albedo
self.albedo_pl = albedo_pl # Bond albedo (scattering) applied to self.toa_heating in socrates.py
self.zenith_angle = zenith_angle # solar zenith angle, Hamano+15 (arccos(1/sqrt(3) = 54.74), Wordsworth+10: 48.19 (arccos(2/3)), see Cronin+14 for definitions

self.planet_mass = pl_mass
self.planet_radius = pl_radius
Expand Down Expand Up @@ -250,6 +258,56 @@ def __init__(self, T_surf: float, P_surf: float, P_top: float, pl_radius: float,
self.liquid_water_fraction = 0.0
self.cloud_fraction = 0.0

#New contructor based on toml file
@classmethod
def from_file(atm, file: str, band_edges: list, vol_mixing: dict = {}, vol_partial: dict = {}):

with open(file, 'r') as f:
cfg = toml.load(f)

# Set parameters according to toml file if key present
# Otherwise set default values
T_surf = cfg['atmos']['T_surf'] if "T_surf" in cfg['atmos'] else 0.0
P_surf = cfg['atmos']['P_surf'] if "P_surf" in cfg['atmos'] else 1.0e5
P_top = cfg['atmos']['P_top'] if "P_top" in cfg['atmos'] else 1.0
pl_radius = cfg['planet']['pl_radius'] if "pl_radius" in cfg['planet'] else 6.371e6
pl_mass = cfg['planet']['pl_mass'] if "pl_mass" in cfg['planet'] else 5.972e24
#optional arguments
req_levels = cfg['atmos']['req_levels'] if "req_levels" in cfg['atmos'] else 100
alpha = cfg['atmos']['alpha'] if "alpha" in cfg['atmos'] else 0.
trppT = cfg['atmos']['trppT'] if "trppT" in cfg['atmos'] else 290.0
do_cloud = cfg['atmos']['do_cloud'] if "do_cloud" in cfg['atmos'] else False
re = cfg['atmos']['re'] if "re" in cfg['atmos'] else 0.
lwm = cfg['atmos']['lwm'] if "lwm" in cfg['atmos'] else 0.
clfr = cfg['atmos']['clfr'] if "clfr" in cfg['atmos'] else 0.
albedo_pl = cfg['atmos']['albedo_pl'] if "albedo_pl" in cfg['atmos'] else 0.175
albedo_s = cfg['atmos']['albedo_s'] if "albedo_s" in cfg['atmos'] else 0.
zenith_angle = cfg['atmos']['zenith_angle'] if "zenith_angle" in cfg['atmos'] else 54.74

return atm(T_surf, P_surf, P_top, pl_radius, pl_mass,
band_edges, vol_mixing, vol_partial,
#optional arguments
req_levels=req_levels,
alpha_cloud=alpha,
trppT=trppT,
do_cloud=do_cloud,
re=re,
lwm=lwm,
clfr=clfr,
albedo_pl=albedo_pl,
albedo_s=albedo_s,
zenith_angle=zenith_angle)

def setSurfaceTemperature(self, Tsurf: float):

self.ts = Tsurf
self.tmp[0] = Tsurf
self.tmpl[0] = Tsurf

def setTropopauseTemperature(self):

T_eqm = (self.instellation * self.inst_sf * (1.0 - self.albedo_pl) /phys.sigma)**(1.0/4.0)
self.trppT = T_eqm * (0.5**0.25) # radiative skin temperature

def write_PT(self,filename: str="output/PT.tsv", punit:str = "Pa"):
"""Write PT profile to file, with descending pressure.
Expand Down
80 changes: 26 additions & 54 deletions tests/test_instellation.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#!/usr/bin/env python3

from importlib.resources import files
import os, shutil
import os, shutil, toml
import numpy as np

from janus.modules.stellar_luminosity import InterpolateStellarLuminosity
Expand All @@ -12,54 +12,6 @@
import janus.utils.phys as phys
from janus.utils.ReadSpectralFile import ReadBandEdges

def run_once(sep, dirs, T_magma, P_surf, skin_d, band_edges):

# Planet
time = { "planet": 0., "star": 100e+6 } # yr,
star_mass = 1.0 # M_sun, mass of star
pl_radius = 6.371e6 # m, planet radius
pl_mass = 5.972e24 # kg, planet mass

# Boundary conditions for pressure & temperature
P_top = 1.0 # Pa

# Define volatiles by mole fractions
vol_mixing = {
"H2O" : 1.0,
"CO2" : 0.0,
"N2" : 0.0
}

rscatter = True
A_B = 0.1 # bond albedo
A_S = 0.1
inst_sf = 3.0/8.0

##### Function calls

S_0 = InterpolateStellarLuminosity(star_mass, time, sep)

zenith_angle = 48.19 # cronin+14 (also for scaling by a factor of 3/8 ^^)

T_eqm = (S_0 * inst_sf * (1.0 - A_B) /phys.sigma)**(1.0/4.0)
T_trpp = T_eqm * (0.5**0.25) # radiative skin temperature

# Create atmosphere object
atm = atmos(T_magma, P_surf * 1e5, P_top, pl_radius, pl_mass, band_edges,
vol_mixing=vol_mixing, trppT=T_trpp, req_levels=150)
atm.albedo_pl = A_B
atm.albedo_s = A_S
atm.inst_sf = inst_sf
atm.zenith_angle = zenith_angle
atm.instellation = S_0
atm.skin_d = skin_d
atm.tmp_magma = T_magma

# Do rad trans
atm = MCPA_CBL(dirs, atm, False, rscatter, T_surf_max=9.0e99, T_surf_guess = T_trpp+100)

return [atm.SW_flux_down[0], atm.LW_flux_up[0], atm.net_flux[0], atm.ts, T_trpp]

def test_instellation():

# Set up dirs
Expand All @@ -84,19 +36,39 @@ def test_instellation():
)
band_edges = ReadBandEdges(dirs["output"]+"star.sf")

# Set parameters
P_surf = 280.0 # surface pressure [bar]
T_magma = 3000.0 # magma temperature [K]
skin_d = 1e-2 # conductive skin thickness [m]
# Open config file
cfg_file = dirs["janus"]+"data/tests/config_instellation.toml"
with open(cfg_file, 'r') as f:
cfg = toml.load(f)

# Planet
time = { "planet": cfg['planet']['time'], "star": cfg['star']['time']}
star_mass = cfg['star']['star_mass']

# Define volatiles by mole fractions
vol_mixing = {
"H2O" : 1.0,
"CO2" : 0.0,
"N2" : 0.0
}

#Get reference data
ref = np.loadtxt(dirs["janus"]+"data/tests/data_instellation.csv",
dtype=float, skiprows=1, delimiter=',')

# Create atmosphere object
atm = atmos.from_file(cfg_file, band_edges, vol_mixing=vol_mixing, vol_partial={})

r_arr = np.linspace(0.3, 1.4, 7) # orbital distance range [AU]
for i in range(7):
print("Orbital separation = %.2f AU" % r_arr[i])
out = run_once(r_arr[i], dirs, T_magma, P_surf, skin_d, band_edges)

atm.instellation = InterpolateStellarLuminosity(star_mass, time, r_arr[i])
atmos.setTropopauseTemperature(atm)

atm = MCPA_CBL(dirs, atm, False, rscatter = True, T_surf_max=9.0e99, T_surf_guess = atm.trppT+100)

out = [atm.SW_flux_down[0], atm.LW_flux_up[0], atm.net_flux[0], atm.ts, atm.trppT]
print(out)
print(ref[i][1:6])
np.testing.assert_allclose(out, ref[i][1:6], rtol=1e-5, atol=0)
Expand Down
71 changes: 29 additions & 42 deletions tests/test_runaway_greenhouse.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/env python3

import os, shutil
import os, shutil, toml
import numpy as np
from matplotlib.ticker import MultipleLocator
from importlib.resources import files
Expand All @@ -13,46 +13,6 @@
import janus.utils.StellarSpectrum as StellarSpectrum
from janus.utils.ReadSpectralFile import ReadBandEdges


def run_once(T_surf, dirs, band_edges):

# Planet
time = { "planet": 0., "star": 4.5e9 } # yr,
star_mass = 1.0 # M_sun, mass of star
mean_distance = 0.5 # au, orbital distance
pl_radius = 6.371e6 # m, planet radius
pl_mass = 5.972e24 # kg, planet mass

# Boundary conditions for pressure & temperature
P_top = 1.0 # Pa

# Define volatiles by mole fractions
P_surf = 300 * 1e5
vol_mixing = {
"H2O" : 1.0,
"CO2" : 0.0,
"N2" : 0.0
}

# Rayleigh scattering on/off
rscatter = False

# Tropopause calculation
trppT = 0.0 # Fixed tropopause value if not calculated dynamically

##### Function calls

# Create atmosphere object
atm = atmos(T_surf, P_surf, P_top, pl_radius, pl_mass, band_edges, vol_mixing=vol_mixing, trppT=trppT)

# Compute stellar heating
atm.instellation = InterpolateStellarLuminosity(star_mass, time, mean_distance)

# Do rad trans
_, atm_moist = RadConvEqm(dirs, time, atm, standalone=True, cp_dry=False, trppD=False, rscatter=rscatter)

return [atm_moist.LW_flux_up[0]]

def test_runaway_greenhouse():

# Set up dirs
Expand All @@ -77,14 +37,41 @@ def test_runaway_greenhouse():
)
band_edges = ReadBandEdges(dirs["output"]+"star.sf")

# Open config file
cfg_file = dirs["janus"]+"data/tests/config_runaway.toml"
with open(cfg_file, 'r') as f:
cfg = toml.load(f)

# Planet
time = { "planet": cfg['planet']['time'], "star": cfg['star']['time']}
star_mass = cfg['star']['star_mass']
mean_distance = cfg['star']['mean_distance']

vol_mixing = {
"H2O" : 1.0,
"CO2" : 0.0,
"N2" : 0.0
}

#Get reference values
OLR_ref = np.loadtxt(dirs["janus"]+"data/tests/data_runaway_greenhouse.csv",
dtype=float, skiprows=1, delimiter=',')

# Create atmosphere object
atm = atmos.from_file(cfg_file, band_edges, vol_mixing=vol_mixing, vol_partial={})

# Compute stellar heating
atm.instellation = InterpolateStellarLuminosity(star_mass, time, mean_distance)

#Run Janus
Ts_arr = np.linspace(200, 2800, 20)
for i in range(20):
out = run_once(Ts_arr[i], dirs, band_edges)

atmos.setSurfaceTemperature(atm, Ts_arr[i])

_, atm_moist = RadConvEqm(dirs, time, atm, standalone=True, cp_dry=False, trppD=False, rscatter=False)

out = [atm_moist.LW_flux_up[0]]
print("Output %s; Reference %s" % (out, OLR_ref[i][1]))
np.testing.assert_allclose(out, OLR_ref[i][1], rtol=1e-5, atol=0)

Expand Down

0 comments on commit 045c178

Please sign in to comment.