Skip to content

Add emission and deposition tables to the budget - #457

Open
Behrooz-Roozitalab wants to merge 1 commit into
NCAR:mainfrom
Behrooz-Roozitalab:chemistry_tables
Open

Add emission and deposition tables to the budget#457
Behrooz-Roozitalab wants to merge 1 commit into
NCAR:mainfrom
Behrooz-Roozitalab:chemistry_tables

Conversation

@Behrooz-Roozitalab

Copy link
Copy Markdown
Contributor

Updates to the budget code including:

  • including 'tropospheric' or global as an option in the yaml file.
  • including an option to create and use annual mean files instead of using the raw data.

Also, two tables are added to this version (could be called in YAML file) to calculate the emissions and depositions for all the variables in the files.

…ons and depositions- for all variables, add a calculation option using annual mean files

@justin-richling justin-richling left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This PR is mostly good @Behrooz-Roozitalab, thank you for breaking this up logically! It looks like a lot of suggested changes, but it's really only related to commented out code and some indentation suggestions. I will merge this when you get to these, thanks again!

@@ -136,7 +153,7 @@ def aerosol_gas_tables(adfobj, trop_val=None, **kwargs):
# if True, calculate only Tropospheric values
# if False, all layers
# tropopause is defiend as either directly or indirectly. Look for tropopause to see the definition
Tropospheric = bres['Tropospheric']
#Tropospheric = bres['Tropospheric']

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for the late review on this @Behrooz-Roozitalab! The only question I have is if this is commented out, could we just remove it, or is this something that might get flagged later?

raise AdfError(errmsg)

# Initialize nicknames dictionary
#nicknames = {}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we remove this if it is not being used?

Comment on lines +191 to +196
# Calculated duration of time period in seconds?
#durations[case] = (end_period-start_period).days*86400 #+365*86400


# Get number of years for calculations
#num_yrs[case] = (int(end_year)-int(start_year)) #+1

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we also remove these commented out lines?

except subprocess.CalledProcessError as e:
print(f"NCO Error (Exit Code {e.returncode}): {e.stderr}")

#os.sys(f"ncra {Files} {output_location}/{case}_ANN.nc")

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this also safe to remove?


# Gather dictionary data for current case
# NOTE: The calculations can take a long time...
#Dic_scn_var_comp[case] = make_Dic_scn_var_comp_2D(adfobj, VARIABLES, data_dir, dic_SE, Files, ext1_SE, AEROSOLS)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we also remove this and line 273, #inside = Inside_SE_region(current_lat,current_lon,dir_shapefile)?

needed_vars = []
Dic_all_data={}

# all_data=[]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we remove this line?

if file == 0:
Dic_all_data[var]=[]

# Star gathering of variable data

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's indent this to match the formatting and fix typo to # Start gathering ...

data=data*delP
else:
data=data
# End if

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like this is also off on the formatting, could you indent/move it to the correct spot?

nickname = case

# Collect row data in a list of dictionaries
#durations[case]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove this commented out line.

@@ -136,7 +153,7 @@ def aerosol_gas_tables(adfobj, trop_val=None, **kwargs):
# if True, calculate only Tropospheric values
# if False, all layers
# tropopause is defiend as either directly or indirectly. Look for tropopause to see the definition
Tropospheric = bres['Tropospheric']
#Tropospheric = bres['Tropospheric']

### NOT WORKING FOR NOW

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this still not currently working, or can we remove this line too?

@brianpm

brianpm commented Aug 28, 2026

Copy link
Copy Markdown
Contributor

A note on this review. This is a new, AI-generated review process we're trying out, aimed at catching more errors and converging on consistent style and conventions across the ADF. It follows the rules in the new AGENTS.md file. If you disagree with any of the conventions it applies, think a finding is wrong, or have other comments on the process, please say so — that feedback is genuinely useful and will shape both the review process and AGENTS.md itself. Thanks!

Summary. Adds two new analysis scripts (emissions_table.py, depositions_table.py, ~1900 lines, both derived from aerosol_gas_tables.py) and reworks aerosol_gas_tables.py to add a Climate fast-path (pre-average history files with ncra) and to move the Tropospheric switch from YAML into a function kwarg. The feature is useful and the PR merges cleanly against main with zero conflicts against any other open PR. Not ready as-is — there is a guaranteed NameError, a copy/paste bug in a budget term, an ncra call that runs twice, and a broad except: → 0 pattern that turns missing model fields into zeros in published tables.

Classification: new/modified diagnostic scripts (scripts/analysis/) plus a defaults change. §4.1, §4.2, §4.3, §4.5, §6 apply. CI does not lint or test any of these files (§3), so everything below is un-caught by the workflows.

Blocking

1. depositions_table.py:248subprocess is never imported.
The file uses subprocess.run at lines 248 and 250 and subprocess.CalledProcessError at 253, but the import block (lines 1–11) has no import subprocessemissions_table.py:4 and aerosol_gas_tables.py:4 both have it, this one was dropped. The import in aerosol_gas_tables.py does not cover this: it binds the name in that module's globals only, and depositions_table imports nothing from it.
Failure: any config running depositions_table with {depositions_table: {kwargs: {Climate: True}}} dies with NameError: name 'subprocess' is not defined on the first case. Climate: False (the default) never reaches the line, which is why it hasn't shown up.

2. aerosol_gas_tables.py:1257 — the CHMP fallback assigns the wrong variable.

try:
    spc_chmp = Dic_scn_var_comp[current_var][current_var+'_CHMP']
except:
    spc_chml = 0                                    # <- should be spc_chmp
spc_chmp = np.where(np.isnan(trop), np.nan, spc_chmp)   # line 1258

Two bugs in three lines: the handler doesn't define spc_chmp, so line 1258 raises UnboundLocalError anyway (the try protects nothing); and it clobbers spc_chml, which was legitimately read at line 1240 and is still needed. Failure: any run whose history files lack <var>_CHMP — the exact case this try was added for.

3. ncra runs twice, and the guarded call is the second one.
aerosol_gas_tables.py:329-335, emissions_table.py:242-248, depositions_table.py:248-254:

subprocess.run(cmd, check=True)                                  # unguarded, runs first
try:
    subprocess.run(cmd, check=True, capture_output=True, text=True)   # same cmd again
except FileNotFoundError: ...
except subprocess.CalledProcessError as e: ...

The first line looks like a leftover. It doubles the cost of the most expensive step in the Climate path (an ncra over every history file in the period), and because it is unguarded it raises before the try ever executes — so the error handlers below are unreachable. Delete the bare call and keep the guarded one. Also §4.5: on a failure the code prints and then falls through to make_Dic_scn_var_comp* with a {case}_ANN.nc that was never written, producing a second, more confusing error — it should warn and fall back to the raw-file path, or continue.

4. The except: <term> = 0 pattern silently fabricates budget numbers (§5.2).
~15 new blocks in aerosol_gas_tables.py (1201, 1210, 1220, 1230, 1241, 1274, 1286, 1297, 1312, 1322, 1335, 1345, 1377, 1394) and the equivalents in the two new files replace a missing field with the scalar 0, which propagates to a rounded 0.0 in the published table. A missing <var>_DDF and a genuinely zero dry deposition are then indistinguishable to the reader, and any total computed from those terms is wrong with no indication.

This is the §4.2 "missing input files are normal" case, but the required behavior there is warn-and-skip, not warn-and-zero. The script already tracks missing_vars_tot — these terms should be omitted from chem_dict (so the column is absent or NaN), with the variable name accumulated into the existing missing-variable message. All of these are also bare except: (§6.3) and will swallow KeyboardInterrupt and any bug in the indexing expression; they should be except KeyError:.

5. New scripts are not registered in config_cam_baseline_example.yaml (§4.1).
emissions_table and depositions_table appear in no config file in the repo. The established pattern is the commented-out opt-in line — config_cam_baseline_example.yaml:515 already has #- aerosol_gas_tables. Without it there is nothing telling a user these exist or how to pass Climate/Tropospheric via the kwargs dict form. (Entry-point names themselves are correct: emissions_table.pyemissions_table(), depositions_table.pydepositions_table(), and neither basename collides on the flat scripts/ sys.path.)

6. Tropospheric removed from adf_variable_defaults.yaml with no migration path (§4.3).
The key is deleted from the budget_tables block and the read at aerosol_gas_tables.py:156 is commented out; it is now a kwarg defaulting to True. Anyone running with a custom variable_defaults file that sets Tropospheric: False will have that setting silently ignored and get tropospheric budgets instead of total-column — no error, just different numbers. The PR description says this option was added to the YAML file; it was in fact removed from it. Needs either a backward-compatible fallback (Tropospheric = kwarg if not None else bres.get('Tropospheric', True)) or a loud note in the PR description per §4.3.

Non-blocking

7. Removing the AREA fallback may break unstructured grids. aerosol_gas_tables.py:452-467 — the old code tried area, then AREA, then a rectilinear dlon/dlat computation; now AREA is gone and a KeyError on area falls straight to the rectilinear branch. For an SE/unstructured grid, lat/lon are 1-D ncol arrays, so np.meshgrid(lon, lat) builds an (ncol, ncol) field and lat[1]-lat[0] is meaningless. The stated motivation (CESM AREA is wrong in some versions) is fair, but the replacement should detect an unstructured grid and warn/skip rather than compute a rectilinear area for it.

8. ncra output lands in the website plot directory. All three scripts write {case}_ANN.nc to Path(adfobj.plot_location[0]). Two issues: it puts a potentially large intermediate data file in the directory that becomes the published website, and [0] is hard-coded regardless of the loop index i, so every case's intermediate goes into the first case's plot dir. A scratch/temp location (or plot_location[i]) would be better, plus cleanup after use.

9. Duration is inconsistent between the two branches. emissions_table.py:183 vs :188Climate: True hard-codes 365*86400, while the raw-file path uses (end_period - start_period).days * 86400 from datetime, i.e. the proleptic Gregorian calendar. For a period spanning leap years those disagree by ~0.07% per leap year, and CAM's default calendar is noleap, so the datetime value is the wrong one. Worth making both 365 * 86400 * n_years.

10. Docstrings do not follow §6.2. Both new files: no module docstring; the entry-point docstrings are free prose with no Parameters / Returns sections and no Notes block listing the adfobj attributes used (these scripts touch variable_defaults, get_cam_info, case_nicknames, climo_yrs, hist_string, plot_location, compare_obs, get_baseline_info, debug_log, add_website_data — exactly what that convention exists to surface). The heading underline at emissions_table.py:21 is also shorter than its title, which Sphinx will drop. §6.2 applies to new code even though the file it was copied from predates the convention.

11. black on new code (§6.1). The two new files are new code, so §6.1 applies in full — currently they carry the source file's spacing (Climate=None,**kwargs, an f-string with no placeholders at emissions_table.py:262). A darker --check --diff --revision main...HEAD pass would cover it. Low priority relative to everything above.

12. YAML formatting and a behavior question. adf_variable_defaults.yaml:2438-2443 — the AEROSOLS continuation lines mix 3- and 4-space indentation, and the diff drops two blank lines around ext1_SE for no reason. Also worth confirming the intent of adding the modal species (so4_a1dst_a3): membership in AEROSOLS selects which branch calc_budget_data takes (DDF+DDFC / WDF+WDFC vs plain DDF/WDF), so this changes results for any existing user who has those in VARIABLES.

Optional / follow-up (pre-existing, not this PR's job)

  • Get_files returns areas unbound when area=False (aerosol_gas_tables.py:446-470) — NameError if ever called that way.
  • Time averaging is unweighted throughout (np.nanmean(..., axis=0) at emissions_table.py:761, and ncra without -w), so monthly means are averaged without month-length weighting. Pre-existing and consistent between the old and new paths, so Climate does not make it worse — but it is a real ~0.1% bias in all of these tables and would make a good standalone issue.
  • env/conda_environment.yaml has no nco entry even though the framework already requires ncrcat/ncatted for time-series generation (lib/adf_diag.py:668). Pre-existing framework gap, worth its own issue.
  • Duplicate import numpy as np at lines 1 and 6 of all three files.

Not verified

  • Whether the Climate: True numbers actually match the raw-file numbers. This needs a real ADF run on chemistry history output; it was not run for this review. Given the duration/num_yrs differences between the branches (finding Mismatch in baseline climo calculation #9), a side-by-side of the two paths on the same case is probably the most valuable check.
  • Science correctness of the new emission/deposition term lists and the so4 → MW 32.066 assignment (aerosol_gas_tables.py:172-176) — a call for the chemistry owner, not a code review. Note the substring test is lowercase 'so4' in var, so it applies to the new modal species but not to SO4, which is presumably already in the MW dict; worth confirming that's deliberate.
  • Whether any downstream consumer depends on budget_tables: Tropospheric in a custom variable-defaults file (finding Should config file have ability to pass configuration arguments directly to scripts? #6) — only this repo was grepped.

Findings #1#3 are mechanical and quick. #4 (the silent zeros) has the real scientific consequence and is the biggest ask. Happy to help with any of these.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants