Initial support for aerosols in PROTEUS, improve chemistry plot, and fix broken links in readme

README.md

@@ -32,11 +32,11 @@
 
 More information can be found on the [documentation](https://proteus-framework.org/PROTEUS/) pages:
 
-* [model description](https://proteus-framework.org/PROTEUS/model.html)
-* [installation guide](https://proteus-framework.org/PROTEUS/installation.html)
-* [usage guide](https://proteus-framework.org/PROTEUS/usage.html)
-* [contributing guide](https://proteus-framework.org/PROTEUS/CONTRIBUTING.html)
+* [model description](https://proteus-framework.org/PROTEUS/Explanations/model.html)
+* [installation guide](https://proteus-framework.org/PROTEUS/How-to/installation.html)
+* [usage guide](https://proteus-framework.org/PROTEUS/How-to/usage.html)
+* [contributing guide](https://proteus-framework.org/PROTEUS/Community/CONTRIBUTING.html)
 
-You can find help on the [discussions page](https://github.com/orgs/FormingWorlds/discussions) or by [contacting the developers](https://proteus-framework.org/PROTEUS/contact.html) directly.
+You can find help on the [discussions page](https://github.com/orgs/FormingWorlds/discussions) or by [contacting the developers](https://proteus-framework.org/PROTEUS/Community/contact.html) directly.
 
-If you make use of PROTEUS, please reference the papers outlined in the [bibliography](https://proteus-framework.org/PROTEUS/bibliography.html).
+If you make use of PROTEUS, please reference the papers outlined in the [bibliography](https://proteus-framework.org/PROTEUS/Reference/bibliography.html).

input/all_options.toml

@@ -207,6 +207,7 @@ version = "2.0"
     prevent_warming = true      # do not allow the planet to heat up
     surface_d       = 0.01      # m, conductive skin thickness
     surface_k       = 2.0       # W m-1 K-1, conductive skin thermal conductivity
+    aerosols_enabled = false    # enable aerosol radiative effects
     cloud_enabled   = false     # enable water cloud radiative effects
     cloud_alpha     = 0.0       # condensate retention fraction (1 -> fully retained)
     surf_state      = "fixed"   # surface scheme: "mixed_layer" | "fixed" | "skin"
@@ -233,6 +234,9 @@ version = "2.0"
         overlap_method  = "ee"          # gas overlap method
         surf_roughness  = 1e-3          # characteristic surface roughness [m]
         surf_windspeed  = 2.0           # characteristic surface wind speed [m/s].
+        cloud_enabled   = false         # include water cloud radiative effects?
+        cloud_alpha     = 0.5           # water condensate retention fraction
+        aerosols_enabled = false        # include aerosol radiative effects?
         rainout         = true          # include volatile condensation/evaporation aloft
         latent_heat     = false         # include latent heat release when `rainout=true`?
         oceans          = true          # form liquid oceans at planet surface?

input/demos/agni.toml

@@ -88,6 +88,9 @@ version = "2.0"
     albedo_pl       = 0.2   	# Bond albedo (scattering)
     module  = "agni"           # Which atmosphere module to use
     rayleigh = false
+    cloud_enabled   = false         # include water cloud radiative effects?
+    cloud_alpha     = 0.5           # water condensate retention fraction
+    aerosols_enabled = true        # include aerosol radiative effects?
 
     [atmos_clim.agni]
         verbosity       = 1             # output verbosity for agni (0:none, 1:info, 2:debug)

src/proteus/atmos_clim/agni.py

@@ -118,6 +118,36 @@ def _determine_condensates(vol_list: list):
     return [v for v in vol_list if v not in ALWAYS_DRY]
 
 
+def _determine_aerosols(dirs: dict) -> list:
+    """
+    Determine which aerosols are available.
+
+    Parameters
+    ----------
+        dirs : dict
+            Dictionary containing paths to directories
+
+    Returns
+    ----------
+        aerosols : list
+            List of available aerosols
+    """
+
+    scattering_dir = os.path.join(dirs['fwl'], 'scattering', 'scattering')
+    if not os.path.isdir(scattering_dir):
+        log.warning(f'Scattering data directory not found: {scattering_dir}')
+        return []
+
+    aerosols = []
+    for f in os.listdir(scattering_dir):
+        if f.endswith('.mon'):
+            aerosols.append(f.replace('.mon', ''))
+    aerosols = sorted(aerosols)
+
+    log.debug(f'Available aerosols: {aerosols}')
+    return aerosols
+
+
 def init_agni_atmos(dirs: dict, config: Config, hf_row: dict):
     """Initialise atmosphere struct for use by AGNI.
 
@@ -197,6 +227,13 @@ def init_agni_atmos(dirs: dict, config: Config, hf_row: dict):
     p_top = config.atmos_clim.agni.p_top
     p_surf = max(p_surf, p_top * 1.1)  # this will happen if the atmosphere is stripped
 
+    # Aerosol species dictionary (set MMR to zero initially)
+    aerosol_species = {}
+    if config.atmos_clim.aerosols_enabled:
+        aerosol_species = {a: 0.0 for a in _determine_aerosols(dirs)}
+        if len(aerosol_species) == 0:
+            log.warning('No data found for aerosol species')
+
     # Setup struct
     succ = jl.AGNI.atmosphere.setup_b(
         atmos,
@@ -218,6 +255,8 @@ def init_agni_atmos(dirs: dict, config: Config, hf_row: dict):
         IO_DIR=io_dir,
         flag_rayleigh=config.atmos_clim.rayleigh,
         flag_cloud=config.atmos_clim.cloud_enabled,
+        flag_aerosol=config.atmos_clim.aerosols_enabled,
+        aerosol_species=convert(jl.Dict, aerosol_species),
         overlap_method=config.atmos_clim.agni.overlap_method,
         albedo_s=config.atmos_clim.surf_greyalbedo,
         surface_material=surface_material,

src/proteus/atmos_clim/common.py

@@ -141,7 +141,15 @@ def read_ncdf_profile(nc_fpath: str, extra_keys: list = [], combine_edges: bool
             continue
 
         # Reading composition
-        if key == 'x_gas':
+        if key == 'gases':
+            gas_l = ds.variables['gases'][:]  # names (bytes matrix)
+            gases = []
+            for igas, gas in enumerate(gas_l):
+                gas_lbl = ''.join([c.decode(encoding='utf-8') for c in gas]).strip()
+                gases.append(gas_lbl)
+            out['gases'] = gases
+
+        elif key == 'x_gas':
             gas_l = ds.variables['gases'][:]  # names (bytes matrix)
             gas_x = ds.variables['x_gas'][:]  # vmrs (float matrix)
 
@@ -150,6 +158,27 @@ def read_ncdf_profile(nc_fpath: str, extra_keys: list = [], combine_edges: bool
                 gas_lbl = ''.join([c.decode(encoding='utf-8') for c in gas]).strip()
                 out[gas_lbl + '_vmr'] = np.array(gas_x[:, igas])
 
+        elif key == 'aerosols':
+            if 'aerosols' in ds.variables.keys():
+                aer_l = ds.variables['aerosols'][:]  # names (bytes matrix)
+                aerosols = []
+                for iaer, aer in enumerate(aer_l):
+                    aer_lbl = ''.join([c.decode(encoding='utf-8') for c in aer]).strip()
+                    if len(aer_lbl) > 0:
+                        aerosols.append(aer_lbl)
+                out['aerosols'] = aerosols
+
+        elif key == 'aer_mmr':
+            if 'aer_mmr' in ds.variables.keys():
+                aer_l = ds.variables['aerosols'][:]  # names (bytes matrix)
+                aer_x = ds.variables['aer_mmr'][:]  # mmrs (float matrix)
+
+                # get data for each aerosol
+                for iaer, aer in enumerate(aer_l):
+                    aer_lbl = ''.join([c.decode(encoding='utf-8') for c in aer]).strip()
+                    if len(aer_lbl) > 0:
+                        out[aer_lbl + '_mmr'] = np.array(aer_x[:, iaer])
+
         else:
             out[key] = np.array(ds.variables[key][:])
 
@@ -158,7 +187,10 @@ def read_ncdf_profile(nc_fpath: str, extra_keys: list = [], combine_edges: bool
 
     # convert to np arrays
     for key in out.keys():
-        out[key] = np.array(out[key], dtype=float)
+        try:
+            out[key] = np.array(out[key], dtype=float)
+        except (AttributeError, TypeError, ValueError):
+            out[key] = np.array(out[key])
 
     return out
 

src/proteus/cli.py

@@ -409,6 +409,14 @@ def surfaces():
     download_surface_albedos()
 
 
+@click.command()
+def scattering():
+    """Get scattering radiative properties"""
+    from .utils.data import download_scattering
+
+    download_scattering()
+
+
 @click.command()
 def reference():
     """Get reference data (exoplanet populations, mass-radius curves, etc.)"""
@@ -462,6 +470,7 @@ def spider():
 cli.add_command(get)
 get.add_command(spectral)
 get.add_command(surfaces)
+get.add_command(scattering)
 get.add_command(muscles)
 get.add_command(phoenix)
 get.add_command(solar)

src/proteus/config/_atmos_clim.py

@@ -18,7 +18,7 @@ def tmp_max_bigger_than_tmp_min(instance, attribute, value):
 
 def warn_if_dummy(instance, attribute, value):
     if (instance.module == 'dummy') and value:
-        raise ValueError('Dummy atmos_clim module is incompatible with Rayleigh scattering')
+        raise ValueError(f'Dummy atmos_clim module is incompatible with {attribute.name}=True')
 
 
 def check_overlap(instance, attribute, value):
@@ -310,6 +310,8 @@ class AtmosClim:
         Conductive skin thickness [m],
     surface_k: float
         Conductive skin thermal conductivity [W m-1 K-1].
+    aerosols_enabled: bool
+        Enable aerosol radiative effects.
     cloud_enabled: bool
         Enable water cloud radiative effects.
     cloud_alpha: float
@@ -346,7 +348,8 @@ class AtmosClim:
     prevent_warming: bool = field(default=False)
     surface_d: float = field(default=0.01, validator=gt(0))
     surface_k: float = field(default=2.0, validator=gt(0))
-    cloud_enabled: bool = field(default=False)
+    aerosols_enabled: bool = field(default=False, validator=warn_if_dummy)
+    cloud_enabled: bool = field(default=False, validator=warn_if_dummy)
     cloud_alpha: float = field(default=0.0, validator=(ge(0), le(1)))
     surf_greyalbedo: float = field(default=0.2, validator=(ge(0), le(1)))
     albedo_pl = field(default=0.0, validator=valid_albedo)

src/proteus/plot/cpl_chem_atmosphere.py

@@ -91,9 +91,12 @@ def plot_chem_atmosphere(
         log.warning('No atmosphere NetCDF files found in output folder')
         return
     nc_fpath = natural_sort(files)[-1]
-    atm_profile = read_ncdf_profile(nc_fpath, extra_keys=['pl', 'tmpl', 'x_gas'])
+    atm_profile = read_ncdf_profile(
+        nc_fpath, extra_keys=['pl', 'tmpl', 'x_gas', 'cloud_mmr', 'aer_mmr', 'aerosols']
+    )
 
     parr = atm_profile['pl'] * 1e-5  # convert to bar
+    tarr = atm_profile['tmpl']  # temperature profile
 
     # Get year
     year = float(nc_fpath.split('/')[-1].split('_atm')[0])
@@ -107,10 +110,13 @@ def plot_chem_atmosphere(
     else:
         has_offchem = False
 
-    # init plot
-    scale = 1.2
-    fig, ax = plt.subplots(1, 1, figsize=(6 * scale, 5 * scale))
+    # init plot with two panels
+    scale = 1.15
+    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5 * scale, 6 * scale))
 
+    # ====================
+    # Panel 1: Gas Species
+    # ====================
     # plot species profiles
     lw = 0.9
     al = 0.8
@@ -131,46 +137,127 @@ def plot_chem_atmosphere(
             xarr = [xarr[0]] + xarr
             if np.amax(xarr) >= xmin:
                 vmr = float(xarr[-1])
-                ax.plot(xarr, parr, ls='dashed', color=col, lw=_lw, alpha=al)
+                ax1.plot(xarr, parr, ls='dashed', color=col, lw=_lw, alpha=al)
 
         # plot from offline chemistry, if available (solid lines)
         if has_offchem and (gas in atm_offchem.keys()):
             xarr = list(atm_offchem[gas].values)
             if np.amax(xarr) >= xmin:
                 vmr = float(xarr[-1])  # prefer vmr from offline chemistry
-                ax.plot(xarr, parr, ls='solid', color=col, lw=_lw, alpha=al)
+                ax1.plot(xarr, parr, ls='solid', color=col, lw=_lw, alpha=al)
 
         # create legend entry and store surface vmr
         if vmr > 0.0:
-            ax.plot(1e30, 1e30, ls='solid', color=col, lw=_lw, alpha=al, label=lbl)
+            ax1.plot([], [], ls='solid', color=col, lw=_lw, alpha=al, label=lbl)
             vmr_surf.append(vmr)
 
-    # Decorate
-    ax.set_xscale('log')
-    ax.set_xlim(left=xmin, right=1.1)
-    ax.set_xlabel('Volume mixing ratio, at t=%.2e yr' % year)
-    ax.xaxis.set_major_locator(LogLocator(numticks=1000))
+    # Add temperature profile to Panel 1 (secondary x-axis at top)
+    ax1_temp = ax1.twiny()
+    ax1_temp.plot(tarr, parr, 'k', lw=1.5, alpha=0.6, label='Temperature')
+    ax1_temp.set_xlabel('Temperature [K]')
+    ax1_temp.set_xlim(left=0)
+
+    # Decorate Panel 1
+    ax1.set_xscale('symlog', linthresh=xmin)
+    ax1.set_xlim(left=0, right=1.1)
+    ax1.set_xlabel('Gas volume mixing ratio')
 
-    ax.set_ylabel('Pressure [bar]')
-    ax.set_yscale('log')
-    ax.set_ylim(bottom=np.amax(parr), top=np.amin(parr))
-    ax.yaxis.set_major_locator(LogLocator(numticks=1000))
+    ax1.set_ylabel('Pressure [bar]')
+    ax1.set_yscale('log')
+    ax1.set_ylim(bottom=np.amax(parr), top=np.amin(parr))
+    ax1.yaxis.set_major_locator(LogLocator(numticks=1000))
 
     # Legend (handles sorted by surface vmr)
-    handles, labels = ax.get_legend_handles_labels()
+    handles, labels = ax1.get_legend_handles_labels()
     order = np.argsort(vmr_surf)[::-1]
-    ax.legend(
+    ax1.legend(
         [handles[idx] for idx in order],
         [labels[idx] for idx in order],
-        loc='lower center',
-        bbox_to_anchor=(0.5, 1),
-        ncols=11,
-        fontsize=9,
-        borderpad=0.4,
+        loc='center left',
+        bbox_to_anchor=(1.00, 0.5),
+        ncols=max(1, len(order) // 18 + 1),
+        fontsize=8,
+        borderpad=0.0,
+        labelspacing=0.3,
+        columnspacing=1.0,
+        handlelength=0.9,
+        handletextpad=0.3,
+        frameon=False,
+    )
+
+    # =============================
+    # Panel 2: Aerosols and Clouds
+    # =============================
+    # Cloud profiles
+    if 'cloud_mmr' in atm_profile.keys():
+        cloud_mmr = atm_profile['cloud_mmr']
+        ax2.plot(
+            [cloud_mmr[0]] + list(cloud_mmr),
+            parr,
+            ls='solid',
+            color=get_colour('cloud'),
+            lw=1.5,
+            alpha=0.7,
+            label='cloud',
+        )
+
+    # Aerosol profiles
+    # Check for individual aerosol species
+    num_aerosols = 0
+    if 'aerosols' in atm_profile.keys():
+        for aer_name in atm_profile['aerosols']:
+            num_aerosols += 1
+            aer_mmr = atm_profile[f'{aer_name}_mmr']
+            ax2.plot(
+                [aer_mmr[0]] + list(aer_mmr),
+                parr,
+                ls='solid',
+                lw=1.5,
+                alpha=0.7,
+                color=get_colour(aer_name),
+                label=aer_name,
+            )
+
+    # Add temperature profile to Panel 2 (secondary x-axis at top)
+    ax2_temp = ax2.twiny()
+    ax2_temp.plot(tarr, parr, 'k', lw=1.5, alpha=0.6)
+    ax2_temp.set_xticks([])  # hide x-tick labels on secondary axis
+    ax2_temp.set_xlim(left=0)
+
+    # Decorate Panel 2
+    ax2.set_xscale('symlog', linthresh=xmin)
+    ax2.set_xlim(left=0, right=1.1)
+    ax2.set_xlabel('Aerosol mass mixing ratio')
+
+    ax2.set_ylabel('Pressure [bar]')
+    ax2.set_yscale('log')
+    ax2.set_ylim(bottom=np.amax(parr), top=np.amin(parr))
+    ax2.yaxis.set_major_locator(LogLocator(numticks=1000))
+
+    ax2.legend(
+        loc='center left',
+        bbox_to_anchor=(1.00, 0.5),
+        ncols=max(1, num_aerosols // 18 + 1),
+        fontsize=8,
+        borderpad=0.0,
         labelspacing=0.3,
         columnspacing=1.0,
-        handlelength=1.5,
+        handlelength=0.9,
         handletextpad=0.3,
+        frameon=False,
+    )
+
+    # Add time annotation
+    ax1.text(
+        0.02,
+        0.98,
+        f't = {year:.2e} yr',
+        transform=ax1.transAxes,
+        ha='left',
+        va='top',
+        fontsize=11,
+        weight='bold',
+        zorder=999,
     )
 
     # Save file