diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml index d92eef1..83156f3 100644 --- a/.github/workflows/publish.yaml +++ b/.github/workflows/publish.yaml @@ -38,12 +38,12 @@ jobs: bump-my-version bump --new-version $GITHUB_REF_NAME patch git commit -am "Bump version to: $GITHUB_REF_NAME" - - name: Push back changes to main and tag + - name: Push back changes to master and tag run: | git tag --force $GITHUB_REF_NAME HEAD git push --force --tags - git switch -C main - git push --set-upstream -f origin main + git switch -C master + git push --set-upstream -f origin master deploy: needs: fix_release_deps diff --git a/examples/SocRadConv.py b/examples/SocRadConv.py index 5e647e7..42cb000 100755 --- a/examples/SocRadConv.py +++ b/examples/SocRadConv.py @@ -16,7 +16,7 @@ mpl.use('Agg') import time as t -import os, shutil +import os, shutil, toml import numpy as np from importlib.resources import files @@ -49,30 +49,16 @@ start = t.time() ##### Settings - - # Planet - time = { "planet": 0., "star": 4e+9 } # yr, - star_mass = 1.0 # M_sun, mass of star - mean_distance = 1.0 # au, orbital distance - pl_radius = 6.371e6 # m, planet radius - pl_mass = 5.972e24 # kg, planet mass - - # Boundary conditions for pressure & temperature - T_surf = 2800.0 # K - P_top = 1.0 # Pa - - # Define volatiles by mole fractions - # P_surf = 300.0 * 1e5 - # vol_partial = {} - # vol_mixing = { - # "CO2" : 0.0, - # "H2O" : 1.0, - # "N2" : 0.0, - # } - - # OR: + cfg_file = dirs["janus"]+"data/tests/config_janus.toml" + with open(cfg_file, 'r'): + cfg = toml.load(cfg_file) + + # Planet + time = { "planet": cfg['planet']['time'], "star": cfg['star']['time']} + star_mass = cfg['star']['star_mass'] + mean_distance = cfg['star']['mean_distance'] + # Define volatiles by partial pressures - P_surf = 0.0 vol_mixing = {} vol_partial = { "H2O" : 9.0e5, @@ -83,43 +69,6 @@ "H2" : 1.0e5, } - - # Stellar heating on/off - stellar_heating = True - - # Rayleigh scattering on/off - rscatter = False - - # Pure steam convective adjustment - pure_steam_adj = False - - # Tropopause calculation - trppD = False # Calculate dynamically? - trppT = 10.0 # Fixed tropopause value if not calculated dynamically - - # Water lookup tables enabled (e.g. for L vs T dependence) - water_lookup = False - - # Surface temperature time-stepping - surf_dt = False - cp_dry = False - # Options activated by surf_dt - cp_surf = 1e5 # Heat capacity of the ground [J.kg^-1.K^-1] - mix_coeff_atmos = 1e6 # mixing coefficient of the atmosphere [s] - mix_coeff_surf = 1e6 # mixing coefficient at the surface [s] - - # Cloud radiation - do_cloud = False - alpha = 1.0 - re = 1.0e-5 # Effective radius of the droplets [m] (drizzle forms above 20 microns) - lwm = 0.8 # Liquid water mass fraction [kg/kg] - how much liquid vs. gas is there upon cloud formation? 0 : saturated water vapor does not turn liquid ; 1 : the entire mass of the cell contributes to the cloud - clfr = 0.8 # Water cloud fraction - how much of the current cell turns into cloud? 0 : clear sky cell ; 1 : the cloud takes over the entire area of the cell (just leave at 1 for 1D runs) - - # Instellation scaling | 1.0 == no scaling - Sfrac = 1.0 - - ##### Function calls - # Tidy directory if os.path.exists(dirs["output"]): shutil.rmtree(dirs["output"]) @@ -129,7 +78,6 @@ print("Inserting stellar spectrum") StellarSpectrum.InsertStellarSpectrum( dirs["janus"]+"data/spectral_files/Dayspring/256/Dayspring.sf", - # dirs["janus"]+"data/spectral_files/Reach/Reach.sf", dirs["janus"]+"data/spectral_files/stellar_spectra/Sun_t4_4Ga_claire_12.txt", dirs["output"] ) @@ -137,25 +85,32 @@ band_edges = ReadBandEdges(dirs["output"]+"star.sf") # Create atmosphere object - atm = atmos(T_surf, P_surf, P_top, pl_radius, pl_mass, band_edges, alpha_cloud=alpha, - vol_mixing=vol_mixing, vol_partial=vol_partial, trppT=trppT, req_levels=100, water_lookup=water_lookup, do_cloud=do_cloud, re=re, lwm=lwm, clfr=clfr) + atm = atmos.from_file(cfg_file, band_edges, vol_mixing=vol_mixing, vol_partial=vol_partial) # Set stellar heating on or off - if stellar_heating == False: + if cfg['star']['stellar_heating'] == False: atm.instellation = 0. else: atm.instellation = InterpolateStellarLuminosity(star_mass, time, mean_distance) print("Instellation:", round(atm.instellation), "W/m^2") # Set up atmosphere with general adiabat - atm_dry, atm = RadConvEqm(dirs, time, atm, standalone=True, cp_dry=cp_dry, trppD=trppD, rscatter=rscatter, do_cloud=do_cloud, pure_steam_adj=pure_steam_adj, surf_dt=surf_dt, cp_surf=cp_surf, mix_coeff_atmos=mix_coeff_atmos, mix_coeff_surf=mix_coeff_surf) + atm_dry, atm = RadConvEqm(dirs, + time, + atm, + standalone=True, + cp_dry=False, + trppD=False, # Calculate dynamically? + rscatter=False, # Rayleigh scattering on/off + pure_steam_adj=False, # Pure steam convective adjustment + surf_dt=False, # Surface temperature time-stepping + # Options activated by surf_dt + cp_surf=1e5, # Heat capacity of the ground [J.kg^-1.K^-1], + mix_coeff_atmos=1e6, # mixing coefficient of the atmosphere [s] + mix_coeff_surf=1e6 # mixing coefficient at the surface [s] + ) # Plot abundances w/ TP structure - if (cp_dry): - ga.plot_adiabats(atm_dry,filename=dirs["output"]+"dry_ga.png") - atm_dry.write_PT(filename=dirs["output"]+"dry_pt.tsv") - plot_fluxes(atm_dry,filename=dirs["output"]+"dry_fluxes.png") - ga.plot_adiabats(atm,filename= dirs["output"]+"moist_ga.png") atm.write_PT(filename= dirs["output"]+"moist_pt.tsv") atm.write_ncdf( dirs["output"]+"moist_atm.nc") diff --git a/examples/demo_instellation.py b/examples/demo_instellation.py index 573b4d1..95d1797 100755 --- a/examples/demo_instellation.py +++ b/examples/demo_instellation.py @@ -8,7 +8,7 @@ from matplotlib.ticker import MultipleLocator from importlib.resources import files -import os, shutil +import os, shutil, toml import numpy as np from janus.modules.stellar_luminosity import InterpolateStellarLuminosity @@ -19,85 +19,10 @@ 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) - - # Plot case - plt.ioff() - fig,ax = plt.subplots(1,1) - ax.plot(atm.tmpl, atm.pl, color='black', lw=2) - ax.set_yscale("log"); ax.invert_yaxis() - ax.set_ylabel("Pressure [Pa]") - ax.set_xlabel("Temperature [K]") - ax.set_title("a = %g AU" % sep) - fig.savefig(dirs["output"]+"/recent.jpg",bbox_inches='tight', dpi=100) - plt.close() - - # Save netcdf - atm.write_ncdf(dirs["output"]+"/recent.nc") - - return [atm.SW_flux_down[0], atm.LW_flux_up[0], atm.net_flux[0], atm.ts, T_trpp] - - if __name__=='__main__': print("Start") - # Planet data - x_planets = { - "Earth": 1.0, - "Venus": 0.723, - "Mercury": 0.387 - } - c_planets = { - "Earth": 'deepskyblue', - "Venus": 'goldenrod', - "Mercury": 'violet' - } - # Set up dirs if os.environ.get('RAD_DIR') == None: raise Exception("Socrates environment variables not set! Have you installed Socrates and sourced set_rad_env?") @@ -105,7 +30,7 @@ def run_once(sep, dirs, T_magma, P_surf, skin_d, band_edges): "janus": str(files("janus"))+"/", "output": os.path.abspath(os.getcwd())+"/output/" } - + # Tidy directory if os.path.exists(dirs["output"]): shutil.rmtree(dirs["output"]) @@ -118,48 +43,76 @@ def run_once(sep, dirs, T_magma, P_surf, skin_d, band_edges): dirs["janus"]+"data/spectral_files/stellar_spectra/Sun_t4_4Ga_claire_12.txt", dirs["output"] ) - band_edges = ReadBandEdges(dirs["output"]+"star.sf") + # 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 + } + + # PARAMETERS - P_surf = 280.0 # surface pressure [bar] - T_magma = 3000.0 # magma temperature [K] - skin_d = 1e-2 # conductive skin thickness [m] - r_inner = 0.3 # inner orbital distane [AU] - r_outer = 1.4 # outer orbital distance [AU] - samples = 7 # number of samples - logx = False # log x-axis? legend = True # make legend? dx_tick = 0.1 # x-tick spacing (set to 0 for automatic) # /PARAMETERS + + # Create atmosphere object + atm = atmos.from_file(cfg_file, band_edges, vol_mixing=vol_mixing, vol_partial={}) + T_magma = atm.tmp_magma #get default value 3000 K but this could be a config option + + r_arr = np.linspace(0.3, 1.4, 7) # orbital distance range [AU] - # Run JANUS in a loop to generate data - if logx: - r_arr = np.logspace( np.log10(r_inner), np.log10(r_outer), samples) - else: - r_arr = np.linspace(r_inner, r_outer, samples) asf_arr = [] # ASF OLR_arr = [] # OLR net_arr = [] # net flux at TOA ts_arr = [] # surface temperature tr_arr = [] # tropopause temperature - for r in r_arr: - print("Orbital separation = %.2f AU" % r) - out = run_once(r, dirs, T_magma, P_surf, skin_d, band_edges) - asf_arr.append(out[0] * -1.0) # directed downwards => minus sign - OLR_arr.append(out[1]) - net_arr.append(out[2]) - ts_arr.append(out[3]) - tr_arr.append(out[4]) - print(" ") + for i in range(7): + print("Orbital separation = %.2f AU" % r_arr[i]) + + 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) + + asf_arr.append(atm.SW_flux_down[0] * -1.0) # directed downwards => minus sign + OLR_arr.append(atm.LW_flux_up[0]) + net_arr.append(atm.net_flux[0]) + ts_arr.append(atm.ts) + tr_arr.append(atm.trppT) + + # Plot case + plt.ioff() + fig,ax = plt.subplots(1,1) + ax.plot(atm.tmpl, atm.pl, color='black', lw=2) + ax.set_yscale("log"); ax.invert_yaxis() + ax.set_ylabel("Pressure [Pa]") + ax.set_xlabel("Temperature [K]") + ax.set_title("a = %.2f AU" %r_arr[i]) + fig.savefig(dirs["output"]+"/profile%.2fAU.jpg"%r_arr[i],bbox_inches='tight', dpi=100) + plt.close() + + # Save netcdf + atm.write_ncdf(dirs["output"]+"/profile%.2fAU.nc"%r_arr[i]) + + print(" ") save_arr = [r_arr, asf_arr, OLR_arr, net_arr, ts_arr, tr_arr] np.savetxt(dirs["output"]+"/data_%dK.csv"%T_magma, np.array(save_arr).T, fmt="%.5e", delimiter=",", header="r [AU], S_0 [W m-2], OLR [W m-2], net [W m-2], ts [K], tr[K] ") - # Setup plot - lw = 2.5 print("Making plots") plt.ioff() @@ -171,16 +124,13 @@ def run_once(sep, dirs, T_magma, P_surf, skin_d, band_edges): ax1.text(r_arr[0], -1.0, "Warming", verticalalignment='center', color='seagreen', weight='bold', zorder=8).set_bbox(dict(facecolor='white', alpha=0.7, linewidth=0)) ax1.axhline(y=0, color='black', lw=0.7, zorder=0) # zero flux - ax1.plot(r_arr, asf_arr, color='royalblue', lw=lw, label='ASF') # absorbed stellar flux (sw down) - ax1.plot(r_arr, OLR_arr, color='crimson', lw=lw, label='OLR') # outgoing longwave radiation (lw up) - ax1.plot(r_arr, net_arr, color='seagreen' , lw=lw, label='Net') # net upward-directed flux + ax1.plot(r_arr, asf_arr, color='royalblue', lw=2.5, label='ASF') # absorbed stellar flux (sw down) + ax1.plot(r_arr, OLR_arr, color='crimson', lw=2.5, label='OLR') # outgoing longwave radiation (lw up) + ax1.plot(r_arr, net_arr, color='seagreen' , lw=2.5, label='Net') # net upward-directed flux - if legend: - ax1.legend(loc='center right', framealpha=1.0) + ax1.legend(loc='center right', framealpha=1.0) ax1.set_yscale("symlog") ax1.set_ylabel("Upward flux [W m$^{-2}$]") - if logx: - ax1.set_xscale("log") ax1.set_xticklabels([]) if dx_tick > 1.0e-10: ax1.xaxis.set_minor_locator(MultipleLocator(dx_tick)) @@ -189,19 +139,27 @@ def run_once(sep, dirs, T_magma, P_surf, skin_d, band_edges): yticks = np.delete(yticks, np.argwhere( yticks == 0.0 )) ax1.set_yticks(yticks) + # Planet data + x_planets = { + "Earth": 1.0, + "Venus": 0.723, + "Mercury": 0.387 + } + c_planets = { + "Earth": 'deepskyblue', + "Venus": 'goldenrod', + "Mercury": 'violet' + } for p in x_planets.keys(): for ax in (ax1,ax2): - ax.axvline(x=x_planets[p], color=c_planets[p], lw=lw*0.5, linestyle='dashed', label=p , zorder=1) + ax.axvline(x=x_planets[p], color=c_planets[p], lw=3, linestyle='dashed', label=p , zorder=1) arr_magma = np.ones(len(r_arr))*T_magma - ax2.plot(r_arr, np.zeros(len(r_arr)), zorder=6, color='silver', lw=lw, label=r"$\tilde{T}_s$") # magma temperature - ax2.plot(r_arr, ts_arr - arr_magma, zorder=6, color='black', lw=lw, label=r"$T_s$") # surface solution temperature + ax2.plot(r_arr, np.zeros(len(r_arr)), zorder=6, color='silver', lw=2.5, label=r"$\tilde{T}_s$") # magma temperature + ax2.plot(r_arr, ts_arr - arr_magma, zorder=6, color='black', lw=2.5, label=r"$T_s$") # surface solution temperature - if legend: - ax2.legend(loc='center right', framealpha=1.0) + ax2.legend(loc='center right', framealpha=1.0) ax2.set_ylabel(r"$T - \tilde{T_s}$ [K]") - if logx: - ax2.set_xscale("log") if dx_tick > 1.0e-10: ax2.xaxis.set_minor_locator(MultipleLocator(dx_tick)) diff --git a/examples/demo_runaway_greenhouse.py b/examples/demo_runaway_greenhouse.py index ac58634..2bc0d64 100755 --- a/examples/demo_runaway_greenhouse.py +++ b/examples/demo_runaway_greenhouse.py @@ -5,7 +5,7 @@ mpl.rcParams.update({'font.size': 12}) import matplotlib.pyplot as plt -import os, shutil +import os, shutil, toml import numpy as np from matplotlib.ticker import MultipleLocator from importlib.resources import files @@ -18,47 +18,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 [T_surf, atm_moist.LW_flux_up[0]] - - if __name__=='__main__': print("Start") @@ -71,7 +30,7 @@ def run_once(T_surf, dirs, band_edges): "janus": str(files("janus"))+"/", "output": os.path.abspath(os.getcwd())+"/output/" } - + # Tidy directory if os.path.exists(dirs["output"]): shutil.rmtree(dirs["output"]) @@ -85,19 +44,44 @@ def run_once(T_surf, dirs, band_edges): dirs["output"] ) print(" ") - 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 + } + # 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 + # Run JANUS in a loop to generate runaway curve print("Running JANUS...") - Ts_arr = [] + Ts_arr = np.linspace(200, 2800, 20) OLR_arr = [] - for Ts in np.linspace(200, 2800, 20): - print("T_surf = %d K" % Ts) - out = run_once(Ts, dirs, band_edges) - Ts_arr.append(out[0]) - OLR_arr.append(out[1]) - print(" ") + for i in range(20): + print("T_surf = %d K" % Ts_arr[i]) + atmos.setSurfaceTemperature(atm, Ts_arr[i]) + + _, atm_moist = RadConvEqm(dirs, time, atm, standalone=True, cp_dry=False, trppD=False, rscatter=False) + + OLR_arr.append(atm_moist.LW_flux_up[0]) + print(" ") + OLR_arr = np.array(OLR_arr) Ts_arr = np.array(Ts_arr) diff --git a/pyproject.toml b/pyproject.toml index 34f4d07..a536520 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,15 +27,16 @@ classifiers = [ keywords = ["exoplanet", "atmosphere"] requires-python = '>=3.10' dependencies = [ - "matplotlib", - "natsort", - "netcdf4", - "numpy", - "pandas", - "scipy", - "seaborn", - "tomlkit", - "f90nml" + 'matplotlib', + 'natsort', + 'netcdf4', + 'numpy', + 'pandas', + 'scipy', + 'seaborn', + 'toml', + 'tomlkit', + 'f90nml' ] [project.urls] diff --git a/src/janus/data/tests/config_instellation.toml b/src/janus/data/tests/config_instellation.toml new file mode 100644 index 0000000..c095d71 --- /dev/null +++ b/src/janus/data/tests/config_instellation.toml @@ -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 ^^) diff --git a/src/janus/data/tests/config_janus.toml b/src/janus/data/tests/config_janus.toml new file mode 100644 index 0000000..00ba9c4 --- /dev/null +++ b/src/janus/data/tests/config_janus.toml @@ -0,0 +1,22 @@ +[planet] +time = 0.0 +pl_radius = 6.371e6 +pl_mass = 5.972e24 + +[star] +time = 4.0e+9 # yr +star_mass = 1.0 # M_sun, mass of star +mean_distance = 1.0 # au, orbital distance +stellar_heating = true + +[atmos] +T_surf = 2800.0 +P_surf = 0.0 +P_top = 1.0 +trppT = 10.0 +req_levels = 100 +do_cloud = false +alpha = 1.0 +re = 1.0e-5 # Effective radius of the droplets [m] (drizzle forms above 20 microns) +lwm = 0.8 # Liquid water mass fraction [kg/kg] - how much liquid vs. gas is there upon cloud formation? 0 : saturated water vapor does not turn liquid ; 1 : the entire mass of the cell contributes to the cloud +clfr = 0.8 # Water cloud fraction - how much of the current cell turns into cloud? 0 : clear sky cell ; 1 : the cloud takes over the entire area of the cell (just leave at 1 for 1D runs) diff --git a/src/janus/data/tests/config_runaway.toml b/src/janus/data/tests/config_runaway.toml new file mode 100644 index 0000000..da5a65a --- /dev/null +++ b/src/janus/data/tests/config_runaway.toml @@ -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 diff --git a/src/janus/modules/compute_moist_adiabat.py b/src/janus/modules/compute_moist_adiabat.py index c449139..db39119 100644 --- a/src/janus/modules/compute_moist_adiabat.py +++ b/src/janus/modules/compute_moist_adiabat.py @@ -18,7 +18,7 @@ import janus.utils.GeneralAdiabat as ga # Moist adiabat with multiple condensibles import janus.utils.socrates as socrates -def compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter=False, do_cloud=False): +def compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter=False): """Compute moist adiabat case Parameters @@ -33,8 +33,6 @@ def compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter=False, do_cloud Calculate tropopause dynamically? rscatter : bool Include Rayleigh scattering? - do_cloud : bool - Include water cloud radiation? """ @@ -45,11 +43,11 @@ def compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter=False, do_cloud atm_moist.rh = compute_Rh(atm_moist) - if do_cloud: + if atm.do_cloud: atm_moist = simple_cloud(atm_moist) # Before radiation, set up the cloud for Socrates using the current PT profile # Run SOCRATES - atm_moist = socrates.radCompSoc(atm_moist, dirs, recalc=False, rscatter=rscatter, do_cloud=do_cloud) + atm_moist = socrates.radCompSoc(atm_moist, dirs, recalc=False, rscatter=rscatter) if standalone == True: print("w/o stratosphere (net, OLR): " + str(round(atm_moist.net_flux[0], 3)) +" , "+str(round(atm_moist.LW_flux_up[0], 3)) + " W/m^2") @@ -63,10 +61,10 @@ def compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter=False, do_cloud # Reset stratosphere temperature and abundance levels atm_moist = set_stratosphere(atm_moist) - if do_cloud: + if atm.do_cloud: atm_moist = simple_cloud(atm_moist) # Update cloud location after previous PT changes # Recalculate fluxes w/ new atmosphere structure - atm_moist = socrates.radCompSoc(atm_moist, dirs, recalc=True, rscatter=rscatter, do_cloud=do_cloud) + atm_moist = socrates.radCompSoc(atm_moist, dirs, recalc=True, rscatter=rscatter) if standalone == True: print("w/ stratosphere (net, OLR): " + str(round(atm_moist.net_flux[0], 3)) + " , " + str(round(atm_moist.LW_flux_up[0], 3)) + " W/m^2") diff --git a/src/janus/modules/dry_adiabat_timestep.py b/src/janus/modules/dry_adiabat_timestep.py index 4c046a2..3080a4b 100644 --- a/src/janus/modules/dry_adiabat_timestep.py +++ b/src/janus/modules/dry_adiabat_timestep.py @@ -22,7 +22,7 @@ import janus.utils.phys as phys # Time integration for n steps -def compute_dry_adiabat(atm, dirs, standalone, rscatter=False, pure_steam_adj=False, surf_dt=False, cp_surf=1e5, mix_coeff_atmos=1e6, mix_coeff_surf=1e6, do_cloud=False): +def compute_dry_adiabat(atm, dirs, standalone, rscatter=False, pure_steam_adj=False, surf_dt=False, cp_surf=1e5, mix_coeff_atmos=1e6, mix_coeff_surf=1e6): # Dry adiabat settings rad_steps = 100 # Maximum number of radiation steps @@ -50,14 +50,14 @@ def compute_dry_adiabat(atm, dirs, standalone, rscatter=False, pure_steam_adj=Fa # Compute radiation, midpoint method time stepping try: - if do_cloud: + if atm.do_cloud: atm_dry = simple_cloud(atm_dry) # Before radiation, set up the cloud for Socrates using the current PT profile - atm_dry = socrates.radCompSoc(atm_dry, dirs, recalc=False, rscatter=rscatter, do_cloud=do_cloud) - dT_dry = atm_dry.net_heating * atm_dry.dt + atm_dry = socrates.radCompSoc(atm_dry, dirs, recalc=False, rscatter=rscatter) + dT_dry = atm_dry.net_heating * atm_dry.dt - # Limit the temperature change per step - dT_dry = np.where(dT_dry > dT_max, dT_max, dT_dry) - dT_dry = np.where(dT_dry < -dT_max, -dT_max, dT_dry) + # Limit the temperature change per step + dT_dry = np.where(dT_dry > dT_max, dT_max, dT_dry) + dT_dry = np.where(dT_dry < -dT_max, -dT_max, dT_dry) # Do the surface balance if surf_dt: diff --git a/src/janus/modules/solve_pt.py b/src/janus/modules/solve_pt.py index 0490fa7..7a94c94 100644 --- a/src/janus/modules/solve_pt.py +++ b/src/janus/modules/solve_pt.py @@ -14,7 +14,7 @@ from janus.modules.dry_adiabat_timestep import compute_dry_adiabat from janus.utils.atmosphere_column import atmos -def RadConvEqm(dirs, time, atm, standalone:bool, cp_dry:bool, trppD:bool, rscatter:bool, do_cloud:bool=False, +def RadConvEqm(dirs, time, atm, standalone:bool, cp_dry:bool, trppD:bool, rscatter:bool, pure_steam_adj=False, surf_dt=False, cp_surf=1e5, mix_coeff_atmos=1e6, mix_coeff_surf=1e6): """Sets the atmosphere to a temperature profile using the general adiabat. @@ -36,8 +36,6 @@ def RadConvEqm(dirs, time, atm, standalone:bool, cp_dry:bool, trppD:bool, rscatt Calculate tropopause dynamically? rscatter : bool Include rayleigh scattering? - do_cloud : bool - Include water cloud radiation? pure_steam_adj : bool Use pure steam adjustment? surf_dt : float @@ -53,13 +51,13 @@ def RadConvEqm(dirs, time, atm, standalone:bool, cp_dry:bool, trppD:bool, rscatt ### Moist/general adiabat - atm_moist = compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter, do_cloud) + atm_moist = compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter) ### Dry adiabat if cp_dry == True: # Compute dry adiabat w/ timestepping - atm_dry = compute_dry_adiabat(atm, dirs, standalone, rscatter, pure_steam_adj, surf_dt, cp_surf, mix_coeff_atmos, mix_coeff_surf, do_cloud) + atm_dry = compute_dry_adiabat(atm, dirs, standalone, rscatter, pure_steam_adj, surf_dt, cp_surf, mix_coeff_atmos, mix_coeff_surf) if standalone == True: print("Net, OLR => moist:", str(round(atm_moist.net_flux[0], 3)), str(round(atm_moist.LW_flux_up[0], 3)) + " W/m^2", end=" ") @@ -99,7 +97,7 @@ def MCPA(dirs, atm, standalone:bool, trppD:bool, rscatter:bool): """ ### Moist/general adiabat - return compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter, do_cloud=atm.do_cloud) + return compute_moist_adiabat(atm, dirs, standalone, trppD, rscatter) def MCPA_CBL(dirs, atm_inp, trppD:bool, rscatter:bool, atm_bc:int=0, T_surf_guess:float=-1, T_surf_max:float=-1, method:int=0, atol:float=1.0e-3): """Calculates the temperature profile using the multiple-condensible pseudoadiabat and steps T_surf to conserve energy. @@ -172,7 +170,7 @@ def ini_atm(Ts): def func(x): print("Evaluating at T_surf = %.1f K" % x) - atm_tmp = compute_moist_adiabat(ini_atm(x), dirs, False, trppD, rscatter, do_cloud=do_cloud) + atm_tmp = compute_moist_adiabat(ini_atm(x), dirs, False, trppD, rscatter) if atm_bc == 0: F_atm = atm_tmp.net_flux[0] @@ -227,7 +225,7 @@ def func(x): print(" Using T_surf = %g K" % T_surf) # Get atmosphere state from solution value - atm = compute_moist_adiabat(ini_atm(T_surf), dirs, False, trppD, rscatter, do_cloud=do_cloud) + atm = compute_moist_adiabat(ini_atm(T_surf), dirs, False, trppD, rscatter) atm.ts = T_surf if atm_bc == 0: diff --git a/src/janus/utils/atmosphere_column.py b/src/janus/utils/atmosphere_column.py index 03ccc37..77d82ef 100644 --- a/src/janus/utils/atmosphere_column.py +++ b/src/janus/utils/atmosphere_column.py @@ -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 @@ -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 @@ -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 @@ -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 @@ -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. diff --git a/src/janus/utils/socrates.py b/src/janus/utils/socrates.py index 3fa3a36..341edbf 100644 --- a/src/janus/utils/socrates.py +++ b/src/janus/utils/socrates.py @@ -19,7 +19,7 @@ def radCompSoc(atm, dirs, recalc, rscatter=False, - rewrite_cfg=True, rewrite_tmp=True, rewrite_gas=True, do_cloud=False): + rewrite_cfg=True, rewrite_tmp=True, rewrite_gas=True): """Runs SOCRATES to calculate fluxes and heating rates Parameters @@ -32,8 +32,6 @@ def radCompSoc(atm, dirs, recalc, rscatter=False, Is this function call a 'recalculation' case accounting for a tropopause? rscatter : bool Include Rayleigh scattering? - do_cloud : bool - Include water cloud radiation? rewrite_cfg : bool Re-write configuration values (e.g. TOA heating) rewrite_PT : bool @@ -144,7 +142,7 @@ def radCompSoc(atm, dirs, recalc, rscatter=False, if check_gas(fthis): nctools.ncout3d( fthis, 0, 0, atm.p, phys.molar_mass['H2O'] / atm.mu * atm.x_gas["H2O"], 'q', longname="q", units='kg/kg') # Clouds - if do_cloud: + if atm.do_cloud: fthis = basename+'.re' if check_tmp(fthis): nctools.ncout3d( fthis, 0, 0, atm.p, atm.re, 're', longname="Effective Radius", units='M') @@ -165,7 +163,7 @@ def radCompSoc(atm, dirs, recalc, rscatter=False, nctools.ncout3d(basename+'.'+vol_lower, 0, 0, atm.p, x_gas_this, vol_lower, longname=vol, units='kg/kg') # Call sequences for run SOCRATES + move data - if do_cloud: + if atm.do_cloud: seq_sw_ex = ["Cl_run_cdf","-B", basename,"-s", runspectral_file, "-R 1", str(atm.nbands), " -ch ", str(atm.nbands), " -S -g 2 -C ", str(atm.cloud_scheme), " -K ", str(atm.cloud_representation), " -d ", str(atm.droplet_type), " -v ", str(atm.solver), " -u", scatter_flag, " -o"] seq_sw_mv = ["fmove", basename,"currentsw"] diff --git a/tests/test_instellation.py b/tests/test_instellation.py index ea37dd6..94d233d 100755 --- a/tests/test_instellation.py +++ b/tests/test_instellation.py @@ -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 @@ -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 @@ -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) diff --git a/tests/test_runaway_greenhouse.py b/tests/test_runaway_greenhouse.py index b5e6a0a..46fbb4d 100755 --- a/tests/test_runaway_greenhouse.py +++ b/tests/test_runaway_greenhouse.py @@ -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 @@ -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 @@ -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)