# Public copy: internal compute identifiers and local environment paths were removed or generalized. Scientific values are retained. Paths may need adapting and some original run inputs are not included. """Interactive 3D globe: FKESSLER baroclinic wave growth, days 0-10. Reads the 16 daily CAM h0 history files from the 15-day the compute environment run and builds an animated Plotly globe with jet ribbons, ascent markers, cloud-water markers placed at their true altitude, and surface precipitation. """ import glob import numpy as np import xarray as xr import plotly.graph_objects as go NDAY = 11 # days 0-10 inclusive CL_THRESHOLD_MG = 1000.0 # Kessler autoconversion threshold ~1 g/kg ALPHA = 0.35 # signed power stretch for the pressure anomaly R_JET = 1.16 HALFW = 0.030 # ---------------------------------------------------------------- load files = sorted(glob.glob("hpc/1f737b81-34c1-4b1b-ae39-9db1cbfd439a/*.cam.h0.*.nc")) ds = xr.concat([xr.open_dataset(f, decode_timedelta=False) for f in files], dim="time", data_vars="all") sub = ds.isel(time=slice(0, NDAY)) lat, lon, lev = sub.lat.values, sub.lon.values, sub.lev.values kjet = int(np.abs(lev - 230).argmin()) PS = sub.PS.values / 100.0 OM = sub.OMEGA.values CL = sub.CLDLIQ.values * 1e6 PRECT = sub.PRECT.values * 86400 * 1000.0 cl_max = CL.max(axis=1) cl_argk = CL.argmax(axis=1) psa = PS - PS.mean(axis=(1, 2), keepdims=True) psa_disp = np.sign(psa) * np.abs(psa) ** ALPHA DISP_MAX = float(np.abs(psa_disp).max()) true_amp = np.abs(psa).max(axis=(1, 2)) # ---------------------------------------------------------------- geometry def sph2cart(latd, lond, r=1.0): la, lo = np.deg2rad(latd), np.deg2rad(lond) return (r * np.cos(la) * np.cos(lo), r * np.cos(la) * np.sin(lo), r * np.sin(la)) def bilin(F, latq, lonq): li = np.interp(latq, lat, np.arange(lat.size)) i0 = np.clip(np.floor(li).astype(int), 0, lat.size - 2) fi = li - i0 lo = np.mod(lonq, 360.0) / 2.5 j0 = np.floor(lo).astype(int) % lon.size fj = lo - np.floor(lo) j1 = (j0 + 1) % lon.size return ((1 - fi) * ((1 - fj) * F[i0, j0] + fj * F[i0, j1]) + fi * ((1 - fj) * F[i0 + 1, j0] + fj * F[i0 + 1, j1])) def trace(U, V, lat0, lon0, nstep=46, dt=2400.0): """RK2 streamline of the (U,V) field on the sphere; returns degrees.""" R = 6.371e6 la, lo = lat0, lon0 path = [(la, lo)] for _ in range(nstep): u, v = bilin(U, la, lo), bilin(V, la, lo) cl_ = max(np.cos(np.deg2rad(la)), 0.15) la_t = la + np.rad2deg(v * dt / R) * 0.5 lo_t = lo + np.rad2deg(u * dt / (R * cl_)) * 0.5 u2, v2 = bilin(U, la_t, lo_t), bilin(V, la_t, lo_t) cl2 = max(np.cos(np.deg2rad(la_t)), 0.15) la = float(np.clip(la + np.rad2deg(0.5 * (v + v2) * dt / R), -88.5, 88.5)) lo = lo + np.rad2deg(0.5 * (u + u2) * dt / (R * 0.5 * (cl_ + cl2))) path.append((la, lo)) return np.array(path) def ribbon_mesh(paths, speeds): """Triangulate many streamlines into one ribbon mesh.""" X, Y, Z, I, J, K, C = [], [], [], [], [], [], [] base = 0 for p, sp in zip(paths, speeds): pts = np.stack(sph2cart(p[:, 0], p[:, 1], R_JET), axis=1) n = pts.shape[0] if n < 3: continue tang = np.gradient(pts, axis=0) tn = np.linalg.norm(tang, axis=1, keepdims=True); tn[tn == 0] = 1 tang /= tn rad = pts / np.linalg.norm(pts, axis=1, keepdims=True) binorm = np.cross(tang, rad) bn = np.linalg.norm(binorm, axis=1, keepdims=True); bn[bn == 0] = 1 binorm /= bn taper = np.minimum(1.0, np.minimum(np.arange(n), n - 1 - np.arange(n)) / 4.0)[:, None] left, right = pts + HALFW * binorm * taper, pts - HALFW * binorm * taper X += list(left[:, 0]) + list(right[:, 0]) Y += list(left[:, 1]) + list(right[:, 1]) Z += list(left[:, 2]) + list(right[:, 2]) C += list(sp) + list(sp) for a in range(n - 1): I += [base + a, base + a + 1] J += [base + n + a, base + n + a + 1] K += [base + n + a + 1, base + a] base += 2 * n return dict(x=X, y=Y, z=Z, i=I, j=J, k=K, intensity=C) seed_lats = np.array([-58, -52, -46, -40, -34, -28, 28, 34, 40, 46, 52, 58], float) seed_lons = np.array([10, 100, 190, 280], float) ribbons = [] for t in range(NDAY): U, V = sub.U.values[t, kjet], sub.V.values[t, kjet] spd = np.sqrt(U ** 2 + V ** 2) paths, sps = [], [] for la0 in seed_lats: for lo0 in seed_lons: p = trace(U, V, la0, lo0) paths.append(p) sps.append(bilin(spd, p[:, 0], p[:, 1])) ribbons.append(ribbon_mesh(paths, sps)) SPD_MAX = max(max(rb["intensity"]) for rb in ribbons) # ---------------------------------------------------------------- marker sets def altitude(k): return 1.02 + 0.34 * (1.0 - np.log10(np.clip(lev[k], 5, 1000)) / 3.0) asc, cld, rain = {}, {}, {} NTOP = 620 for t in range(NDAY): om = OM[t] band = (lev[:, None, None] > 150) & (lev[:, None, None] < 900) omm = np.where(band & np.ones_like(om, bool), om, 0.0) idx = np.argsort(omm.ravel())[:NTOP] kk, ii, jj = np.unravel_index(idx, om.shape) s = -omm[kk, ii, jj] keep = s > 1e-4 kk, ii, jj, s = kk[keep], ii[keep], jj[keep], s[keep] x, y, z = sph2cart(lat[ii], lon[jj], altitude(kk)) asc[t] = dict(x=x, y=y, z=z, s=s) cm, ka = cl_max[t], cl_argk[t] sig = cm > 50.0 bg = (cm > 0.5) & ~sig dec = np.zeros_like(bg); dec[::2, ::3] = True ii, jj = np.where(sig | (bg & dec)) kk = ka[ii, jj] x, y, z = sph2cart(lat[ii], lon[jj], altitude(kk)) cld[t] = dict(x=x, y=y, z=z, s=cm[ii, jj]) ii, jj = np.where(PRECT[t] > 0.05) x, y, z = sph2cart(lat[ii], lon[jj], 1.008) rain[t] = dict(x=x, y=y, z=z, s=PRECT[t][ii, jj]) ASC_MAX = max(float(asc[t]["s"].max()) for t in range(NDAY)) CL_MAXG = max(float(cld[t]["s"].max()) for t in range(NDAY)) PR_MAX = float(PRECT.max()) # ---------------------------------------------------------------- traces LA2, LO2 = np.meshgrid(lat, lon, indexing="ij") Xs, Ys, Zs = sph2cart(LA2, LO2, 1.0) R2 = lambda a, n=3: np.round(np.asarray(a, float), n) def msize(v, vmax, lo, hi, p=0.45): return lo + (hi - lo) * np.clip(np.asarray(v, float) / vmax, 0, 1) ** p def frame_traces(t): return [ go.Surface(x=R2(Xs, 4), y=R2(Ys, 4), z=R2(Zs, 4), surfacecolor=R2(psa_disp[t]), cmin=-DISP_MAX, cmax=DISP_MAX, colorscale="RdBu_r", showscale=(t == 0), colorbar=dict(title=dict(text="surface pressure
anomaly (hPa,
amplitude-stretched)", side="right"), x=0.02, len=0.42, y=0.26, thickness=12, tickvals=[-3.9, -2, 0, 2, 3.9], ticktext=["−50", "−7", "0", "+7", "+50"]), lighting=dict(ambient=0.85, diffuse=0.4, specular=0.05), hoverinfo="skip", name="wave"), go.Mesh3d(x=R2(ribbons[t]["x"], 4), y=R2(ribbons[t]["y"], 4), z=R2(ribbons[t]["z"], 4), i=ribbons[t]["i"], j=ribbons[t]["j"], k=ribbons[t]["k"], intensity=R2(ribbons[t]["intensity"], 2), cmin=8, cmax=SPD_MAX, colorscale=[[0, "#3d3f57"], [0.5, "#7f86b8"], [1, "#eaf0ff"]], showscale=False, opacity=0.92, hoverinfo="skip", name="jets"), go.Scatter3d(x=R2(asc[t]["x"], 4), y=R2(asc[t]["y"], 4), z=R2(asc[t]["z"], 4), mode="markers", marker=dict(size=msize(asc[t]["s"], ASC_MAX, 1.6, 5.2), symbol="diamond", color=R2(asc[t]["s"], 4), cmin=0, cmax=ASC_MAX, colorscale=[[0, "#ffe0b2"], [0.35, "#ff9b3d"], [1, "#c62828"]], showscale=False, opacity=0.85), hoverinfo="skip", name="rising air"), go.Scatter3d(x=R2(cld[t]["x"], 4), y=R2(cld[t]["y"], 4), z=R2(cld[t]["z"], 4), mode="markers", marker=dict(size=msize(cld[t]["s"], CL_MAXG, 1.3, 8.0, p=0.38), color=R2(cld[t]["s"], 2), cmin=0, cmax=CL_THRESHOLD_MG, colorscale=[[0, "#f2f6fa"], [0.25, "#dfe8f2"], [0.6, "#a8bdd4"], [0.999, "#4a6fa5"], [1, "#7b2ff7"]], showscale=(t == 0), colorbar=dict(title=dict(text="cloud liquid
water (mg kg⁻¹)", side="right"), x=0.98, len=0.42, y=0.72, thickness=12, tickvals=[0, 250, 500, 750, 1000], ticktext=["0", "250", "500", "750", "1000 = rain
threshold
(not reached)"]), opacity=0.55), hoverinfo="skip", name="cloud water"), go.Scatter3d(x=R2(rain[t]["x"], 4), y=R2(rain[t]["y"], 4), z=R2(rain[t]["z"], 4), mode="markers", marker=dict(size=msize(rain[t]["s"], PR_MAX, 3.0, 11.0, p=0.5), color="#0d47a1", opacity=0.9, line=dict(width=0.5, color="#82b1ff")), hoverinfo="skip", name="precipitation"), ] def readout(t): return ("Day %d peak wave |PS′| %.2f hPa peak ascent %.3f Pa s⁻¹
" "peak cloud water %.0f mg kg⁻¹ = %.0f%% of the 1000 mg kg⁻¹ rain threshold " "peak rain %.1f mm day⁻¹") % ( t, true_amp[t], float(asc[t]["s"].max()), float(cld[t]["s"].max()), 100 * float(cld[t]["s"].max()) / CL_THRESHOLD_MG, float(rain[t]["s"].max()) if rain[t]["s"].size else 0.0) def note(x, y, text, color, anchor="left"): return dict(x=x, y=y, xref="paper", yref="paper", text=text, showarrow=False, font=dict(size=11, color=color), align="left", xanchor=anchor, bgcolor="rgba(12,14,22,0.55)", borderpad=3) legend_notes = [ note(0.012, 0.965, "wind ribbons — streamlines traced in the
232 hPa flow; two midlatitude jets", "#cfd6ea"), note(0.012, 0.885, "rising air ◆ — 620 strongest ascending
cells; colour = true −ω (Pa s⁻¹)", "#ffab52"), note(0.988, 0.965, "condensation → cloud water ● at the
altitude of each column's maximum;
violet = nearing the rain threshold
(peak reaches 99%, never crosses)", "#a8c4e8", "right"), note(0.988, 0.855, "precipitation ● at the surface,
below the strongest ascent", "#82b1ff", "right"), ] frames = [go.Frame(data=frame_traces(t), name=str(t), layout=go.Layout(annotations=[dict(text=readout(t), x=0.5, y=0.055, xref="paper", yref="paper", showarrow=False, align="center", font=dict(size=12, color="#e8eaf0"))])) for t in range(NDAY)] steps = [dict(method="animate", label=str(t), args=[[str(t)], dict(mode="immediate", frame=dict(duration=520, redraw=True), transition=dict(duration=0))]) for t in range(NDAY)] fig = go.Figure(data=frame_traces(0), frames=frames) fig.update_layout( title=dict(text="FKESSLER f19_f19_mg17 — a baroclinic wave amplifies into cyclones and starts to rain
" "CESM 2.1.5 / CAM with Kessler warm-rain physics, days 0–10 of a 15-day run on the compute environment. " "Pressure anomaly is amplitude-stretched (signed |x|0.35) so the weak early wave is visible; " "all colour scales encode true values.", x=0.5, xanchor="center", y=0.965, font=dict(size=15, color="#f0f2f8")), paper_bgcolor="#0c0e16", plot_bgcolor="#0c0e16", scene=dict(xaxis=dict(visible=False, range=[-1.35, 1.35]), yaxis=dict(visible=False, range=[-1.35, 1.35]), zaxis=dict(visible=False, range=[-1.35, 1.35]), aspectmode="cube", bgcolor="#0c0e16", camera=dict(eye=dict(x=1.55, y=1.15, z=0.85))), margin=dict(l=0, r=0, t=92, b=96), annotations=legend_notes + [dict(text=readout(0), x=0.5, y=0.055, xref="paper", yref="paper", showarrow=False, align="center", font=dict(size=12, color="#e8eaf0"))], updatemenus=[dict(type="buttons", showactive=False, x=0.045, y=0.115, xanchor="left", direction="left", bgcolor="#232838", bordercolor="#3b4260", font=dict(color="#e8eaf0", size=11), buttons=[dict(label="▶ play days 0–10", method="animate", args=[None, dict(frame=dict(duration=760, redraw=True), fromcurrent=True, transition=dict(duration=0), mode="immediate")]), dict(label="❚❚ pause", method="animate", args=[[None], dict(frame=dict(duration=0, redraw=False), mode="immediate")])])], sliders=[dict(active=0, steps=steps, x=0.30, len=0.62, y=0.10, xanchor="left", currentvalue=dict(prefix="model day ", font=dict(color="#e8eaf0", size=13)), bgcolor="#232838", bordercolor="#3b4260", activebgcolor="#6f8fd6", font=dict(color="#b8c0d8", size=10), tickcolor="#5a6280")], showlegend=False) fig.write_html("fkessler_baroclinic_wave_globe.html", include_plotlyjs="cdn", auto_play=False) import sys if "--png" in sys.argv: for t in (1, 6, 10): f = go.Figure(data=frame_traces(t)) f.update_layout(paper_bgcolor="#0c0e16", margin=dict(l=0, r=0, t=34, b=0), showlegend=False, title=dict(text=f"day {t} — peak |PS′| {true_amp[t]:.2f} hPa, " f"peak cloud {float(cld[t]['s'].max()):.0f} mg/kg, " f"rain cells {rain[t]['x'].size}", x=0.5, font=dict(color="#f0f2f8", size=12)), scene=dict(xaxis=dict(visible=False), yaxis=dict(visible=False), zaxis=dict(visible=False), aspectmode="cube", bgcolor="#0c0e16", camera=dict(eye=dict(x=1.55, y=1.15, z=0.85)))) f.write_image(f"preview_day{t}.png", width=880, height=780, scale=1)