Debug terminator test in cam-sima

[K20shores]

If you plot the terminator with foo data made from cam-sima, timestep 0 looks correct, and everything after that looks wrong for Cl2. The data has been plotted two different ways from @nusbaume, and my much more complicated script which produces the below plots.

Both indicate that something went wrong with the terminator and foo chemistry in cam-sima. We don't know if this is something wrong with the data collection, generation, micm configuration, micm api usage, of if it's something wrong in cam-sima. Investigate

Acceptance criteria

• terminator with foo produces stable results for Cl, Cl2 and foo
  • we should still likely see a sunspot pattern for Cl2 and no hotspots outside of the sunspot (I think...)

Ideas

• Recreate a 24 hour simulation with a 3d field where each box is independent (so just a bunch of box models) where photolysis is driven by tuvx and solving with micm. This should produce good values for terminator + foo and be nearly identical if not identical to what we get from cam-sima (assuming temperature, pressure, air density, concentrations have the exact same values)
  • Base this off of the python chapman example and include it in our list of python examples because it will still be useful
• Recreate the cam-sima example. Probably, we can get the environmental fields and update the above script to use the same values. Compare the differences are start debugging

import xarray as xr
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import numpy as np
import cartopy.crs as ccrs

ds = xr.load_dataset('sima-terminator/chemistry-demo-foo/chemistry-demo.cam.h1i.0001-01-01-00000.nc')

Cl = ds["Cl"].values      # (time, lev, ncol)
Cl2 = ds["Cl2"].values    # (time, lev, ncol)
Foo = ds["Foo"].values    # (time, lev, ncol)
lat = ds["lat"].values
lon = ds["lon"].values
lev = ds["lev"].values

t_idx = 5
lev_idx = 19

# Compute change
# Cl2_change = Cl2[t_idx] - Cl2[0] if t_idx > 0 else Cl2[t_idx] - 0
# Cl2_change = Cl2[t_idx] - 1.95e-12
Cl2_change = Cl2[t_idx]


# ---------------------------
# 1. Wrap longitude to [-180, 180)
# ---------------------------
lon_wrapped = ((lon + 180) % 360) - 180
idx = np.argsort(lon_wrapped)
lon_wrapped = lon_wrapped[idx]
lat_sorted  = lat[idx]

Cl_sorted  = Cl[t_idx, lev_idx, idx]
Cl2_sorted = Cl2_change[lev_idx, idx]
Foo_sorted = Foo[t_idx, lev_idx, idx]

# ---------------------------
# 2. Make data cyclic
# ---------------------------
lon_cyclic = np.concatenate([lon_wrapped - 360, lon_wrapped, lon_wrapped + 360])
lat_cyclic = np.concatenate([lat_sorted, lat_sorted, lat_sorted])

Cl_cyclic  = np.concatenate([Cl_sorted,  Cl_sorted,  Cl_sorted])
Cl2_cyclic = np.concatenate([Cl2_sorted, Cl2_sorted, Cl2_sorted])
Foo_cyclic = np.concatenate([Foo_sorted, Foo_sorted, Foo_sorted])

# ---------------------------
# 3. Regular lat/lon grid
# ---------------------------
lon_grid = np.linspace(-180, 180, 360)
lat_grid = np.linspace(lat.min(), lat.max(), 180)
lon_mesh, lat_mesh = np.meshgrid(lon_grid, lat_grid)

# ---------------------------
# 4. Interpolate
# ---------------------------
Cl_grid = griddata((lon_cyclic, lat_cyclic), Cl_cyclic, (lon_mesh, lat_mesh), method="linear")
Cl2_grid = griddata((lon_cyclic, lat_cyclic), Cl2_cyclic, (lon_mesh, lat_mesh), method="linear")
Foo_grid = griddata((lon_cyclic, lat_cyclic), Foo_cyclic, (lon_mesh, lat_mesh), method="linear")

# ---------------------------
# 5. Plot
# ---------------------------
fig, axes = plt.subplots(
    1, 3, 
    figsize=(12, 5), 
    constrained_layout=True, 
    subplot_kw=dict(projection=ccrs.PlateCarree(central_longitude=-75)),
    dpi=300
)

hour = f'Hour {ds.time[t_idx].item().strftime("%Y-%m-%d %H:%M:%S")}'
level = f'at {ds.ilev[lev_idx].item():0.2f} hPa'
label = f'{hour} {level}'
fig.suptitle(f"{label}")

im1 = axes[0].imshow(
    Cl_grid, extent=(lon.min(), lon.max(), lat.min(), lat.max()),
    origin="lower", aspect="auto", cmap="viridis",
    transform=ccrs.PlateCarree()
)
im2 = axes[1].imshow(
    Cl2_grid, extent=(lon.min(), lon.max(), lat.min(), lat.max()),
    origin="lower", aspect="auto", cmap="plasma",
    transform=ccrs.PlateCarree()
)
im3 = axes[2].imshow(
    Foo_grid, extent=(lon.min(), lon.max(), lat.min(), lat.max()),
    origin="lower", aspect="auto", cmap="plasma",
    transform=ccrs.PlateCarree()
)

axes[0].set_title("Cl")
axes[1].set_title("$\\Delta \\mathrm{Cl}_2$")
axes[2].set_title("Foo")

for ax in axes:
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.coastlines(color='grey')

cb1 = plt.colorbar(im1, ax=axes[0], location='bottom')
cb2 = plt.colorbar(im2, ax=axes[1], location='bottom')
cb3 = plt.colorbar(im3, ax=axes[2], location='bottom')

reaction_text = (
    r"$\mathrm{Cl_2} \rightarrow 2\,\mathrm{Cl}$" "\n"
    r"$\mathrm{Cl} + \mathrm{Cl} \rightarrow \mathrm{Cl_2}$" "\n"
    r"$\mathrm{Cl_2} \rightarrow \mathrm{Cl} + \mathrm{foo}$" "\n"
    r"$\mathrm{Cl} + \mathrm{foo} \rightarrow \mathrm{Cl_2}$" "\n"
    r"$\mathrm{foo} + \mathrm{foo} \rightarrow \mathrm{Cl_2}$"
)

fig.text(
    0.94, 0.2, reaction_text,
    ha='center', va='bottom',
    fontsize=10,
    bbox=dict(facecolor='white', alpha=0.8, boxstyle='round,pad=0.3')
)

fig.savefig(f'cl_deltacl2_foo_t={t_idx}.jpg')