diff --git a/discontinuity_analysis/.gitignore b/discontinuity_analysis/.gitignore
new file mode 100644
index 0000000..cfb19a2
--- /dev/null
+++ b/discontinuity_analysis/.gitignore
@@ -0,0 +1,8 @@
+# latexmk / TeX intermediates
+*.aux
+*.log
+*.out
+*.toc
+*.fls
+*.fdb_latexmk
+*.synctex.gz
diff --git a/discontinuity_analysis/README.md b/discontinuity_analysis/README.md
new file mode 100644
index 0000000..850f0aa
--- /dev/null
+++ b/discontinuity_analysis/README.md
@@ -0,0 +1,49 @@
+# Discontinuity analysis of the tcpyPI pressure solver (issue #77 / PR #77)
+
+**The definitive writeup.** `discontinuity_analysis.pdf` (build:
+`latexmk -pdf`) is a self-contained mathematical account of the period-2
+convergence failure in the `pcmin` fixed-point iteration, written for a
+reader with no atmospheric-physics background. Every headline number was
+additionally reproduced by a second, independently written implementation of
+the map and CAPE routines.
+
+Summary of findings:
+
+- The failing sounding is the environment of **Hurricane Francine** (18 UTC
+  2024-09-11, Gulf of Mexico, hours before landfall; NHC best track confirms
+  the position) under a −1.6 K counterfactual SST.
+- The fixed-point map `g` is discontinuous; its only diagonal crossing falls
+  inside a 0.837 hPa jump at x_d ≈ 950.8289947515 hPa, so **no fixed point
+  exists** and the iteration locks into a superattracting, bit-exact 2-cycle.
+- The jump comes from one **discontinuous selection rule** in `cape()`
+  (integrate signed buoyancy to the topmost positive grid level, clamp at
+  zero): a one-bit sign decision at 350 hPa — on a feature worth
+  ~1e-2 J/kg, and ballistically unreachable (W = −61 J/kg there) — toggles
+  287 J/kg of negative area and flips eyewall CAPE between 0 and 224 J/kg.
+  Not quadrature, not grid resolution, not interpolation smoothness (each
+  tested and refuted), not floating point.
+- Smooth toy models separate the genuine discontinuities of the underlying
+  variational problem (argmax jumps with continuous value — Berge; reachable
+  energy at barrier tangencies, which is free-convection physics, not PI's
+  forced eyewall ascent) from the code's artificial one.
+- Jumps of `g` are generic (20% of 1,200 sampled real columns); the fatal
+  diagonal-in-gap placement is rare (~0.04 K of SST for this column) but is
+  found systematically by exactly the attribution-style parameter sweeps that
+  produced the case.
+- Repair: CAPE = max_q W(q) (value function; same quadrature and candidates,
+  one changed selection) + residual-checked bracketed root-finding with the
+  explicit non-hypercane endpoint hypothesis. Exact corrected roots:
+  x* = 950.3136889226 hPa (VMAX 83.9065 m/s, PMIN 920.7517 hPa) for #77 and
+  x* = 960.7687 hPa (VMAX 74.4149 m/s) for the 1995 case. On the
+  7,093-column benchmark: zero flag changes, 98.1% of VMAX within 0.1 m/s,
+  1.9% fragile-tail changes up to 7.9 m/s. The period-2 "rescue" becomes
+  dead code (its midpoint was never a root).
+
+Layout:
+
+- `discontinuity_analysis.tex` / `.pdf` — the document (build: `latexmk -pdf`)
+- `figures/` — all figures (PDF used by the doc, PNG previews)
+- `scripts/` — the analysis scripts behind every number and figure; data
+  files (`*.npz`) are regenerated by running `datagen.py`, `humphunt.py`,
+  `refine.py`, `fixdata.py`, `validate.py` against branch
+  `upstream/modernize_eg`
diff --git a/discontinuity_analysis/discontinuity_analysis.pdf b/discontinuity_analysis/discontinuity_analysis.pdf
new file mode 100644
index 0000000..b04270d
Binary files /dev/null and b/discontinuity_analysis/discontinuity_analysis.pdf differ
diff --git a/discontinuity_analysis/discontinuity_analysis.tex b/discontinuity_analysis/discontinuity_analysis.tex
new file mode 100644
index 0000000..f527633
--- /dev/null
+++ b/discontinuity_analysis/discontinuity_analysis.tex
@@ -0,0 +1,1180 @@
+\documentclass[11pt,a4paper]{article}
+
+\usepackage[margin=2.7cm]{geometry}
+\usepackage[T1]{fontenc}
+\usepackage{lmodern}
+\usepackage{amsmath,amssymb,amsthm}
+\usepackage{graphicx}
+\usepackage{booktabs}
+\usepackage{array}
+\usepackage[dvipsnames]{xcolor}
+\usepackage[colorlinks=true,linkcolor=MidnightBlue,citecolor=MidnightBlue,urlcolor=MidnightBlue]{hyperref}
+% TeX Live 2025's microtype 3.2a warns `Command \showhyphens has changed'
+% against this kernel even in an empty document.  Mute the kernel warning
+% during microtype's begin-document check only (hooks run in registration
+% order: mute -> microtype's check -> restore).
+\makeatletter
+\AddToHook{begindocument/before}{\let\ma@saved@warning\@latex@warning@no@line
+  \let\@latex@warning@no@line\@gobble}
+\makeatother
+\usepackage{microtype}
+\makeatletter
+\AddToHook{begindocument}{\let\@latex@warning@no@line\ma@saved@warning}
+\makeatother
+\emergencystretch=1.5em
+
+\graphicspath{{figures/}}
+
+\newtheorem{proposition}{Proposition}
+\newtheorem{observation}{Observation}
+\newtheorem{remark}{Remark}
+\theoremstyle{definition}
+\newtheorem{definition}{Definition}
+\newtheorem{convention}{Convention}
+
+\newcommand{\Etop}{E_{\mathrm{top}}}
+\newcommand{\Emax}{E_{\mathrm{max}}}
+\newcommand{\Ereach}{E_{\mathrm{reach}}}
+\newcommand{\hPa}{\,\mathrm{hPa}}
+\newcommand{\Jkg}{\,\mathrm{J\,kg^{-1}}}
+\newcommand{\K}{\,\mathrm{K}}
+\newcommand{\ms}{\,\mathrm{m\,s^{-1}}}
+
+\title{An artificial discontinuity in the tcpyPI pressure solver:\\
+a mathematical anatomy of issue \#77}
+\author{Analysis prepared from the \texttt{tcpyPI} source,\\
+the PR~\#77 regression data, and the 2024 sample benchmark}
+\date{July 12, 2026}
+
+\begin{document}
+\maketitle
+
+\begin{abstract}
+\noindent
+The potential-intensity solver in \texttt{tcpyPI} (a Python port of Emanuel's
+\texttt{pcmin}) computes a central pressure by fixed-point iteration
+$x_{n+1}=g(x_n)$ and, for the sounding of PR~\#77 --- Hurricane Francine's
+pre-landfall environment in a counterfactual-SST sensitivity member --- locks
+into a bit-exact period-2 cycle instead of converging.  We show that the map
+$g$ is \emph{discontinuous}, that its only diagonal crossing falls inside a
+downward jump of $0.84\hPa$ at $x_d\simeq950.8289947515\hPa$ --- so
+\emph{no fixed point exists} --- and we locate the source precisely: a
+discontinuous \emph{selection rule} inside the CAPE routine (``integrate the
+signed buoyancy up to the topmost grid level where buoyancy is positive, then
+clamp at zero'').  A one-bit sign decision at the $350\hPa$ level, taken on a
+marginal buoyancy feature worth $O(10^{-2})\Jkg$, toggles the inclusion of
+$287\Jkg$ of negative area and moves the eyewall CAPE between $0$ and
+$224\Jkg$.  Smooth toy models prove the pathology is a property of the
+selection functional, not of grid resolution, quadrature, or interpolation,
+and separate it from the \emph{genuine} discontinuities of the underlying
+variational problem (the maximizer, by Berge's maximum theorem; the reachable
+energy at barrier tangencies).  The mechanism is generic --- $20\%$ of
+sampled real columns carry at least one jump of $g$ --- while the fatal
+alignment with the diagonal is rare ($\sim0.04\K$ of SST for this column),
+which is exactly why parameter sweeps find it.  The repair is to define CAPE
+as the \emph{value function} $\max_q W(q)$ of the running signed work --- the
+convention that also matches the physics of forced eyewall ascent --- and to
+solve the pressure equation with a residual-checked bracketed method.  On
+both known failing soundings this converges immediately (exact roots
+$950.3137$ and $960.7687\hPa$); on a 7{,}093-column benchmark it changes no
+convergence flags and leaves $98.1\%$ of intensities within $0.1\ms$, and it
+renders the period-2 ``rescue'' logic dead code.
+\end{abstract}
+
+\tableofcontents
+
+%======================================================================
+\section{Conclusions first}\label{sec:conclusions}
+%======================================================================
+
+\begin{enumerate}
+\item \textbf{The failure is exactly reproducible and is not slow
+convergence.}  The iteration alternates bit-exactly between
+$A=950.653345443832$ and $B=951.279007983853\hPa$; the separation
+($0.626\hPa$) exceeds the $0.5\hPa$ stopping tolerance forever.  Raising the
+iteration cap cannot help.
+\item \textbf{No fixed point exists.}  The map $g$ jumps downward by
+$0.837\hPa$ at $x_d\simeq950.829\hPa$; the residual $g(x)-x$ passes
+from $+0.507$ to $-0.331\hPa$ across the jump without taking the value zero.
+A sign-change search reports $x_d$, but that is the location of a jump, not
+a root.  The jump's size is no accident: it is $\sim\!200\Jkg$ of clamped
+eyewall CAPE times the pressure sensitivity
+$g\,c_1/(R_d\overline T_v)\approx3.7\times10^{-3}\hPa$ per $\Jkg$ ---
+generically $0.7$--$0.9\hPa$, just above the $0.5\hPa$ tolerance
+(\S\ref{sec:jumpmag}).
+\item \textbf{The cause is one discontinuous selection.}  In \texttt{cape()},
+the level of neutral buoyancy is taken to be the \emph{highest sampled level
+with positive buoyancy}.  For the eyewall parcel this one-bit decision moves
+the terminal level between $750$ and $350\hPa$, dragging $287\Jkg$ of
+intervening negative area into a signed integral that is then clamped at
+zero: eyewall CAPE flips between $224.23\Jkg$ and $0$.  The saturated CAPE,
+outflow temperature, and outflow level are smooth there --- the discontinuity
+is not in the outflow calculation and not a quadrature error.  The rule is
+inherited verbatim from Emanuel's \texttt{pc\_min.m}: this is a defect of
+the reference algorithm itself, not of the Python port.
+\item \textbf{The same mechanism explains the historical 1995 failure}
+(cycle $961.017\leftrightarrow961.671\hPa$, triggered by the $400\hPa$
+sample; eyewall CAPE flips $0\leftrightarrow203\Jkg$), while an adjacent-day
+control column has an ordinary continuous root and converges.
+\item \textbf{The repair is one changed selection plus one changed solver.}
+Define CAPE as the maximum of the \emph{same} running signed work over
+candidate terminal levels ($\Emax=\max_qW(q)$; Berge's maximum theorem makes
+the value continuous even when the maximizing level jumps), and solve
+$g(x)=x$ with residual-checked bracketing.  Negative-buoyancy layers still
+count against the work; nothing is skipped or free.  Exact corrected results:
+$x_*=950.3136889226\hPa$, $V_{\max}=83.9065\ms$, $P_{\min}=920.7517\hPa$ for
+\#77, and $x_*=960.7687\hPa$, $V_{\max}=74.4149\ms$ for the 1995 case, each
+confirmed by two independent implementations of the map.  The upstream
+period-2 ``rescue''
+(midpoint collapse) is then dead code; its midpoint $\approx950.966\hPa$ was
+never a root of the stated equation.
+\end{enumerate}
+
+%======================================================================
+\section{The computation, stripped of meteorology}
+%======================================================================
+
+\noindent\fbox{\begin{minipage}{0.965\linewidth}
+\small
+\textbf{Glossary for the mathematician}\par\smallskip
+\textbf{Pressure as vertical coordinate.}  In atmospheric science the
+vertical coordinate is pressure $p$, measured in hectopascals (hPa); $p$
+decreases monotonically with altitude, from $\sim 1000\hPa$ at sea level to
+$0$ at the top of the atmosphere.  ``Above'' means ``smaller $p$''.  Data
+comes on a finite grid of \emph{pressure levels} $p_0 > p_1 > \dots > p_N$.
+
+\textbf{Parcel.}  An idealized blob of air lifted adiabatically from a
+starting level through the environment.  Its temperature along the ascent is
+determined by thermodynamics (conservation of entropy), \emph{not} by the
+environment; the environment only enters through comparison.
+
+\textbf{Buoyancy.}  $b(p)$ = (density temperature of the parcel at level $p$)
+$-$ (density temperature of the environment at $p$), in Kelvin.  Positive
+$b$ pushes the parcel up.  Density temperature is a smooth algebraic function
+of temperature, moisture and pressure.
+
+\textbf{CAPE} (convective available potential energy, $\Jkg$): the work the
+buoyancy force can do on the parcel, an integral of $b$ described below.
+
+\textbf{LNB} (level of neutral buoyancy): the level where an ascending
+parcel's buoyancy returns to zero --- where it ``detrains'' (exits).  The
+temperature there, the \emph{outflow temperature} $T_O$, sets the Carnot
+efficiency of the storm-as-heat-engine.
+
+\textbf{SST}: sea-surface temperature.  \textbf{MSL}: sea-level pressure far
+from the storm.
+\end{minipage}}
+
+\subsection{The work integral and the CAPE conventions}
+
+Fix a column: levels $p_0>\dots>p_N$ (the code retains levels with
+$p_j>50\hPa$; $N+1=28$ for the soundings here), environmental temperatures
+$T_j$ and moisture $r_j$.  For a parcel specified by its starting pressure,
+temperature and moisture, thermodynamics produces a buoyancy value $b_j$ at
+every level: below its condensation level by an explicit smooth formula,
+above it by inverting a smooth, strictly monotone entropy relation (a
+guarded Newton iteration in the code).
+
+\begin{observation}\label{obs:cont}
+Each $b_j$ is a continuous --- indeed, away from one benign switch,
+real-analytic --- function of the parcel parameters.  (The Newton stopping
+tolerance of $10^{-3}\K$ and a formula switch at the condensation level
+introduce micro-jumps of order $\lesssim10^{-3}\K$ in $b_j$; these propagate
+to $\lesssim10^{-2}\hPa$ in the pressure map, $50\times$ below the solver
+tolerance, and are accounted for separately in \S\ref{sec:micro}.)
+\end{observation}
+
+Define the \emph{running work integral} from the parcel's starting pressure
+$q_0$ up to a trial ``outflow'' pressure $q\le q_0$:
+\begin{equation}\label{eq:W}
+W(q) \;=\; R_d \int_{q}^{q_0} b(p)\, \mathrm{d}\ln p ,
+\end{equation}
+where $R_d=287$ in SI units and $b(p)$ is the piecewise-linear interpolant of
+the $b_j$.  The code evaluates \eqref{eq:W} by the trapezoid rule
+$R_d(b_j+b_{j-1})(p_{j-1}-p_j)/(p_{j-1}+p_j)$ on grid intervals, plus a
+partial term up to an interpolated zero of $b$ at the top --- i.e., exactly
+the integral of the piecewise-linear reconstruction, with the terminal point
+located \emph{sub-grid} by linear interpolation.  $W$ is a $C^1$-in-$q$,
+piecewise-quadratic function whose critical points are the zeros of $b$:
+local maxima at down-crossings ($b:+\to-$, candidate LNBs) and local minima
+at up-crossings.  Figure~\ref{fig:anatomy} shows the generic anatomy: one
+positive region of $b$, one hump of $W$, an unambiguous LNB at the argmax.
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.86\linewidth]{fig_healthy_anatomy.pdf}
+\caption{A generic strong-PI column (January, eastern tropical Pacific,
+$10^{\circ}$N $108^{\circ}$W; the saturated parcel of \S\ref{sec:map},
+$V_{\max}=72.5\ms$).  (a)~Buoyancy of the lifted parcel at the 28 retained
+grid levels: positive through the troposphere, crossing zero between the
+$100$ and $70\hPa$ levels.  (b)~The running work integral $W(q)$: a single
+dominant hump; CAPE ($7764\Jkg$) is its maximum, attained at the
+interpolated LNB ($91.2\hPa$).  For such profiles every convention below
+coincides.}
+\label{fig:anatomy}
+\end{figure}
+
+\begin{remark}[the $p_{\mathrm{top}}$ truncation can clip the LNB]
+\label{rem:ptop}
+The retained column ends at the last level above $p_{\mathrm{top}}=50\hPa$
+--- here $70\hPa$.  For the most intense columns the saturated parcel is
+\emph{still buoyant at $70\hPa$} (its true zero crossing lies above the
+cutoff); the code then sets $\mathrm{INB}$ to the top retained level and
+reports LNB $=70\hPa$, $T_O=T(70\hPa)$.  In the 2024 benchmark this affects
+$23/7093$ converging columns ($0.3\%$) --- precisely the extreme tail, every
+one with $V_{\max}>86\ms$.  The bias has two competing parts: the true LNB
+sits higher, where the stratospheric environment is \emph{warmer}, so the
+clipped $T_O$ is too cold (efficiency, hence $V_{\max}$, overestimated),
+while the uncaptured positive area above $70\hPa$ underestimates CAPE.  It
+cannot be quantified from the shipped sample, whose data itself stops at
+$50\hPa$; assessing it requires source data with higher levels (available
+in ERA5, and present in the PR~\#77 profile) plus a $p_{\mathrm{top}}$
+sensitivity test --- bearing in mind that the reversible parcel model
+itself becomes questionable in the stratosphere, which is presumably why
+\texttt{pcmin} truncates at all, and that under the legacy selector a
+deeper column would add \emph{more} marginal pockets to flicker (under
+max-work it is harmless).  This truncation bias is orthogonal to the
+discontinuity studied here, but a robust API should flag a clipped LNB
+($\mathrm{INB}$ at the profile top) rather than report it as an ordinary
+outflow level.
+\end{remark}
+
+\begin{remark}[the two negative tails of $W$]\label{rem:tails}
+In the $W$ figures the curve dips negative at \emph{both} ends, for
+different reasons.  Above the LNB, $b<0$ and $W$ falls away from its
+maximum --- mildly at first, then steeply in the stratosphere (to
+$\approx-8000\Jkg$ by $70\hPa$ in the \#77 sounding; the figures clip this)
+--- which is exactly what makes the LNB the argmax.  Near $1000\hPa$ the
+negativity is bookkeeping, not ascent physics: $W$ is anchored at the
+parcel's launch pressure $q_0=\min(x,1000)$, which for the storm-core
+parcels lies \emph{above} the lowest data level because the trial central
+pressure $x$ is below $1000\hPa$.  For $q>q_0$ (below the launch point)
+$W(q)=-R_d\int_q^{q_0}b\,\mathrm d\ln p<0$ when the parcel is buoyant at the
+bottom: displacing it downward costs work.  The code carries this segment as
+a signed ``partial area between the parcel pressure and the lowest level''
+correction (present identically in \texttt{pcmin.m}); the parcel never
+visits it.  The anchor itself, $W(q_0)=0$, is a candidate terminal point ---
+which is why $\Emax\ge0$ below holds with no clamp despite both tails.
+\end{remark}
+
+Turning $W$ into a single number requires choosing a terminal point, and
+here conventions diverge.  Three natural choices:
+
+\begin{convention}[value function]\label{conv:max}
+$\displaystyle \Emax=\max_{q\in[p_N,\,q_0]} W(q)$, the value of the
+parametric maximization; the LNB is the argmax.  Since $q=q_0$ is a
+candidate with $W(q_0)=0$, automatically $\Emax\ge0$; no clamping is needed.
+\end{convention}
+
+\begin{convention}[ballistic / reachable]\label{conv:reach}
+$\Ereach=\max\{W(q):q\in\mathcal R\}$, where $\mathcal R$ is the
+connected component of $\{W>0\}$ containing $q_0$.  A parcel starting
+\emph{at rest} converts work
+to kinetic energy and stalls where $W$ first returns to zero; $\Ereach$ is
+the best energy among levels such a parcel can reach unaided.
+\end{convention}
+
+\begin{convention}[topmost-positive-level --- the code's rule]\label{conv:top}
+Let $\mathrm{INB}=\max\{\,j : b_j>0\,\}$ (a \emph{grid} index).  The
+code integrates the signed work up to level $\mathrm{INB}$, adds the
+partial area up to the interpolated zero above it, and clamps the result at
+zero:
+\[
+\Etop=\max\bigl(W(z_{\mathrm{top}}),\,0\bigr),
+\]
+where $z_{\mathrm{top}}$ is the topmost down-crossing of the
+piecewise-linear $b$ (or the profile top if $b_N>0$).  The clamp
+$\max(\,\cdot\,,0)$ turns a possibly negative signed integral into $0$; it
+\emph{is} reachable here and pins one branch of the issue-77 flip at
+exactly $0$ --- though the discontinuity itself lives in the selector, not
+the clamp, which in fact attenuates the jump (\S\ref{sec:jumpmag}).
+\end{convention}
+
+To be clear about provenance: Convention~\ref{conv:top} is \emph{not} a
+porting divergence.  The selector, the clamp, and every formula above are
+inherited line-for-line from Emanuel's \texttt{pc\_min.m} (whose comment
+reads ``Find maximum level of positive buoyancy, INB''), translated
+verbatim in the very first \texttt{pyPI} commit (April 2020) and never
+altered since.  Neither the CAPE definition of Emanuel (1994,
+Eq.~6.3.6) nor Bister and Emanuel (2002) specifies a rule for profiles
+with multiple buoyant layers \cite{emanuel1994,bisteremanuel2002}; the
+topmost-positive-node selection is an undocumented implementation choice
+of the reference code, which the Python port faithfully reproduces ---
+including its discontinuity.
+
+Pointwise, $\Etop \le \Emax$ and $\Ereach\le\Emax$:
+Convention~\ref{conv:top} evaluates $W$ at one particular candidate,
+Convention~\ref{conv:max} at the best one.  Figure~\ref{fig:conventions}
+shows the four archetypal profiles that span all the ways the three can
+relate.  (a)~For a single-hump $W$ --- the overwhelmingly common case ---
+all three agree exactly.  They separate only when the buoyancy has a second
+positive layer aloft, and then in three distinct ways.  (b)~If the upper
+cell is \emph{dominated} (its local maximum of $W$ is lower) but no energy
+barrier intervenes ($W>0$ throughout), then $\Emax=\Ereach$ at the lower
+cell while $\Etop$, bound to the \emph{topmost} crossing, returns the
+smaller upper value: it undercounts by charging the full traversal cost of
+the valley for a destination that was not worth reaching.  (c)~If the upper
+feature is marginal \emph{and} sits behind a barrier ($W<0$ in between) ---
+the issue-77 configuration --- the signed integral to the top is negative
+and the clamp returns $\Etop=0$, while $\Emax=\Ereach$ sit unperturbed at
+the lower cell; worse, as \S\ref{sec:toy} shows, $\Etop$ here flips
+discontinuously between the two regimes as the marginal feature's sign
+flickers.  (d)~If the upper cell is \emph{dominant} but behind a barrier,
+$\Emax=\Etop$ at the top while the ballistic $\Ereach$ stalls at the lower
+cell: the one configuration where the reachability convention parts company
+with the other two.
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_conventions.pdf}
+\caption{The three CAPE conventions on four archetypal smooth profiles
+(shaded bands: layers of positive buoyancy; curves: the running work
+integral $W$).  Markers show each convention's terminal point and value:
+$\Emax$ (dot), $\Ereach$ (open square), $\Etop$ (cross).  (a)~One positive
+region: all agree.  (b)~Dominated upper cell, no barrier:
+$\Etop<\Emax=\Ereach$.  (c)~Marginal bump behind a barrier (the issue-77
+shape): $\Etop$ is clamped to $0$ from a negative signed integral;
+$\Emax=\Ereach$ ignore the bump.  (d)~Dominant upper cell behind a barrier:
+$\Ereach<\Emax=\Etop$.}
+\label{fig:conventions}
+\end{figure}
+
+\paragraph{Which convention does the physics of PI want?}
+The word ``CAPE'' is misleading here.  In an ordinary sounding it is read as
+kinetic energy available to a freely rising parcel, for which the ballistic
+Convention~\ref{conv:reach} is the natural semantics.  The PI quantity is
+instead \emph{hurricane CAPE}: in the Bister--Emanuel derivation the CAPE
+difference is an alternative expression for the thermodynamic work around
+the idealized storm cycle \cite{bisteremanuel2002,gilford2021}, and eyewall
+ascent in a mature vortex is driven primarily by frictional convergence ---
+organized, forced, predominantly hydrostatic --- not by a parcel coasting on
+its own buoyancy \cite{garner2015}.  An intervening negatively buoyant layer
+must still be \emph{charged against the work} (it reduces $W$), but it is
+not a hard reachability barrier for the forced secondary circulation.  That
+is precisely Convention~\ref{conv:max}: every layer costs what it costs,
+and the outflow level is the one maximizing extractable work.  This
+distinction drives the recommendation in \S\ref{sec:fix}.
+
+\subsection{The fixed-point map}\label{sec:map}
+
+The solver seeks the storm's central pressure $x$ (hPa) at the radius of
+maximum winds.  Given a trial $x$, it computes three CAPE values with the
+environment held fixed: $E_A$ (an $x$-independent ambient surface parcel),
+$E_M(x)$ (the ambient surface parcel displaced to pressure $x$, which
+enriches its moisture --- the ``eyewall parcel''), and $E_{MS}(x)$ (a parcel
+saturated at the SST at pressure $x$ --- the storm-core parcel, whose LNB
+also furnishes the outflow temperature $T_O(x)$).  These combine into
+\begin{equation}\label{eq:g}
+\begin{aligned}
+g(x) \;&=\; \mathrm{MSL}\cdot
+\exp\!\Bigl(-\,\mathrm{CAT}(x)\,/\,\bigl(R_d\,\overline{T}_v(x)\bigr)\Bigr),\\
+\mathrm{CAT}(x)\;&=\;\max\Bigl((E_M{-}E_A)+\tfrac12\,\kappa\,
+\tfrac{\mathrm{SST}}{T_O}\,(E_{MS}{-}E_M),\ 0\Bigr),
+\end{aligned}
+\end{equation}
+with $\kappa=C_k/C_D=0.9$ and $\overline{T}_v$ a smooth positive function of
+$x$.  The released code runs Picard iteration $x_{n+1}=g(x_n)$ from
+$x_0=970$, stopping when consecutive iterates satisfy
+$|x_{n+1}-x_{n}|\le 0.5\hPa$, capped at 200 steps.  (A two-back comparison
+$|x_{n+1}-x_{n-1}|$ appears only in the recently added cycle-detection
+patch, not in the stopping rule.)  All CAPE evaluations use
+Convention~\ref{conv:top}.  Note that even for continuous $g$, Picard
+iteration requires local contractivity $|g'(x_*)|<1$; existence of a fixed
+point alone is not enough --- and for discontinuous $g$, neither existence
+nor convergence is available.
+
+Away from the pathologies below, $g$ is a mild contraction: on the failing
+sounding its slope off the jump is $0.34$--$0.46$, and generic columns
+converge geometrically (Fig.~\ref{fig:healthymap}).
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.5\linewidth]{fig_healthy_map.pdf}
+\caption{The map $g$ for the generic column of Fig.~\ref{fig:anatomy}, with
+the cobweb of the iteration from $x_0=970$: contraction to a fixed point at
+$944.6\hPa$ in seven steps.}
+\label{fig:healthymap}
+\end{figure}
+
+%======================================================================
+\section{A smooth toy model: what is genuinely discontinuous}
+\label{sec:toy}
+%======================================================================
+
+Everything essential can be exhibited with $C^\infty$ data, which settles at
+the outset that grid resolution, quadrature order, and interpolation
+smoothness are \emph{not} the cause (a claim tested directly on the real
+sounding in \S\ref{sec:notcause}).  Take a smooth family of buoyancy
+profiles on a height coordinate $t$,
+\begin{equation}\label{eq:toy}
+b(t)\;=\;\underbrace{0.6\,e^{-\left(\frac{t-0.22}{0.10}\right)^2}}_{\text{main
+convective cell}}\;-\;\underbrace{\delta\,
+e^{-\left(\frac{t-0.5}{0.12}\right)^2}}_{\text{stable ``valley''}}\;+\;
+\underbrace{\varepsilon\,
+e^{-\left(\frac{t-0.8}{0.06}\right)^2}}_{\text{upper feature}},
+\end{equation}
+with $W(t)=\int_0^t b$, and evaluate the three conventions along two sweeps
+(Fig.~\ref{fig:toy}).  The sweeps animate the archetypes of
+Fig.~\ref{fig:conventions}: sweep 1 carries the profile into and out of
+archetype~(c) as the marginal bump's sign flickers; sweep 2 deforms
+archetype~(b) into archetype~(d) as the valley deepens and the barrier
+closes.
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_toy.pdf}
+\caption{The smooth toy model \eqref{eq:toy}.
+\textbf{Sweep 1} (top row; $\delta=1.1$ fixed): the amplitude $\varepsilon$ of
+a \emph{marginal} upper bump crosses zero.  $\Etop$ jumps from $0.095$ to $0$
+at $\varepsilon=0^+$: the instant the bump's peak becomes positive, the
+topmost zero of $b$ teleports above it and the terminal value
+$W\approx-0.13<0$ is clamped to $0$.  $\Emax$ and $\Ereach$ do not react at
+all --- the bump is dominated \emph{and} unreachable.  This is the issue-77
+mechanism in pure form.
+\textbf{Sweep 2} (bottom row; two genuinely competing cells, valley depth
+$\delta$ varies): $\Emax$ is continuous with a kink where the argmax
+exchanges humps ($\delta=0.881$), while $\Ereach$ has a \emph{genuine}
+downward jump at $\delta=0.534$, where the energy barrier closes and the
+upper cell becomes ballistically unreachable.}
+\label{fig:toy}
+\end{figure}
+
+The toy isolates three distinct behaviors:
+
+\begin{enumerate}
+\item \textbf{$\Emax$ is continuous, always} --- Berge's maximum theorem in
+its simplest form: $W(\cdot;\lambda)$ is a jointly continuous family over a
+fixed compact interval, so $\lambda\mapsto\max_q W(q;\lambda)$ is
+continuous \emph{even where the argmax jumps}.  The argmax correspondence is
+only upper hemicontinuous: at a dominance exchange the maximizer teleports
+(by $455\hPa$ in the real-data example of \S\ref{sec:physical}), but the
+value passes through with a kink, not a jump.
+\item \textbf{$\Ereach$ can jump for real physical reasons}, but only at a
+codimension-one tangency: the local minimum of $W$ in the valley touches
+zero, the set $\{W>0\}$ changes topology, and the reachable energy jumps.
+A bona-fide first-order transition --- but of the \emph{free-convection}
+problem, which is not the PI problem (\S\ref{sec:physical}).
+\item \textbf{$\Etop$ jumps at events with no physical content}: the sign of
+$b$ at its topmost zero is toggled by an $O(\varepsilon)$ feature carrying
+$O(\varepsilon)$ energy, dominated by orders of magnitude, and (in sweep 1)
+sitting behind a closed barrier.  The discontinuity locus of the
+``topmost positive point'' selector strictly contains the physical loci of
+items 1--2 and is enormously larger.  The same holds for smooth fields in
+general: an upper positive pocket is born at a fold ($b=0$, $\partial_t b=0$
+simultaneously), and the uppermost-root diagnostic jumps there even though
+$b$ is $C^\infty$ and the value function is not affected unless the new
+pocket actually dominates.
+\end{enumerate}
+
+\subsection{Why the symptom is a bit-exact 2-cycle}
+
+Insert a downward jump into an otherwise flat contraction so that the gap
+straddles the diagonal: $g(x)-x>0$ on the left branch, $<0$ on the right.
+No fixed point exists, every orbit must cross the jump forever, and because
+$g$ is nearly constant on each side (slope $\ll1$), the branch values
+$g(A)\approx B$ and $g(B)\approx A$ form a \emph{superattracting} 2-cycle:
+$g\circ g$ contracts strongly onto each branch point, so the alternation
+converges bitwise (Fig.~\ref{fig:toymap}).  The observed 16-digit lock is
+the canonical dynamics of a gap-straddling contraction --- an attracting
+border-collision cycle, not a curiosity of the sounding.
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.4\linewidth]{fig_toymap.pdf}
+\caption{Abstract gap map: piecewise-linear contraction (slope $0.25$) with a
+downward jump straddling the diagonal.  Orbits fall into the unique 2-cycle
+$\{g(B),g(A)\}$; the cycle is superattracting because $|g'(A)g'(B)|\ll1$.}
+\label{fig:toymap}
+\end{figure}
+
+%======================================================================
+\section{The failure in the wild: Hurricane Francine, issue \#77}
+\label{sec:pr77}
+%======================================================================
+
+\subsection{Provenance}
+
+The failing sounding is the environment of an \emph{actual tropical
+cyclone}.  PR~\#77 records $28.6^{\circ}$N, $92.1^{\circ}$W at 18~UTC on
+September 11, 2024 --- and the National Hurricane Center best track places
+Hurricane Francine at exactly that position and time, at 80~kt a few hours
+before its Louisiana landfall \cite{francine,pr77}.  The thermodynamic
+column is a \emph{counterfactual-climate sensitivity member} associated with
+the storm (SST perturbed by $-1.6\K$, mapped to a 1994/day-263 reference),
+not an in-situ measurement --- a distinction that matters physically but not
+to the numerical diagnosis.  The failure thus arises while asking the
+operational question this code exists to answer: \emph{how strong could this
+storm have been in a cooler climate?}  For that member the solver returns
+``missing''.  Column parameters: SST $=28.203^{\circ}$C, MSL
+$=1014.97\hPa$, 37 levels (28 retained).
+
+\subsection{The map has no fixed point}
+
+Scanning $g$ at $0.002\hPa$ resolution over $[946,958]$
+(Fig.~\ref{fig:pr77map}) yields:
+
+\begin{observation}\label{obs:pr77}
+On $[946,958]$ the map $g$ has exactly one discontinuity, a downward jump of
+$0.837\hPa$ at $x_d\simeq950.8289947515\hPa$.  The residual $F(x)=g(x)-x$
+satisfies $F>0$ for $x<x_d$ and $F<0$ for $x>x_d$ (Fig.~\ref{fig:capem}b):
+the only sign change occurs across the jump, so $g$ has no fixed point.
+Away from $x_d$, $g$ is smooth with slope $0.34$--$0.46$.  The iterates
+converge to the 2-cycle $\{A,B\}$ straddling $x_d$.
+\end{observation}
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_pr77_map.pdf}
+\caption{The issue-77 map (Hurricane Francine's environment, $-1.6\K$
+counterfactual).  (a)~$g$ on $[946,958]$ with the iteration cobweb from
+$x_0=970$; the red box marks the zoom.  (b)~Zoom: two smooth branches
+separated by a $0.837\hPa$ gap; the diagonal passes through the gap; the
+iteration locks into the bit-exact 2-cycle (red squares).}
+\label{fig:pr77map}
+\end{figure}
+
+\subsection{One-sided anatomy of the jump}
+
+Evaluating all internal quantities $10^{-9}\hPa$ on either side of $x_d$
+isolates the discontinuous term:
+
+\begin{center}
+\small
+\begin{tabular}{lrr}
+\toprule
+quantity & $x_d-10^{-9}$ & $x_d+10^{-9}$ \\
+\midrule
+buoyancy sample at $350\hPa$ (K) & $+7.58\times10^{-11}$ & $-7.99\times10^{-11}$ \\
+highest positive index & 18 ($350\hPa$) & 10 ($750\hPa$) \\
+diagnosed LNB (hPa) & $350.0000$ & $741.4664$ \\
+eyewall CAPE $E_M$ ($\Jkg$) & $0$ & $224.23146$ \\
+saturated CAPE $E_{MS}$ ($\Jkg$) & $8565.00966$ & $8565.00966$ \\
+saturated outflow $T_O$ (K) & $206.1096884$ & $206.1096884$ \\
+$g(x)$ (hPa) & $951.3363394$ & $950.4982353$ \\
+$F(x)=g(x)-x$ (hPa) & $+0.5073447$ & $-0.3307594$ \\
+\bottomrule
+\end{tabular}
+\end{center}
+
+\noindent
+The saturated terms --- and hence the outflow temperature entering
+\eqref{eq:g} --- are smooth; the discontinuity is entirely in the eyewall
+parcel's $E_M$, through the selector of Convention~\ref{conv:top}.
+(The $10^{-11}\K$ values demonstrate the one-sided limit; they are not the
+amplitude of the atmospheric feature.  At the two cycle phases the $350\hPa$
+buoyancy is $+0.0137$ and $-0.0350\K$: marginal and energetically
+negligible, but not floating-point noise.)
+
+Figure~\ref{fig:mech} shows the geometry.  The eyewall parcel's buoyancy has
+a healthy positive region below $750\hPa$ (the main hump of $W$, worth
+$\approx222$--$224\Jkg$ depending on the trial pressure), a deep negative
+valley from $750$ to $400\hPa$ ($-287\Jkg$), and the marginal excursion at
+the $350\hPa$ node.  When that node's sign flickers positive, the selector
+teleports its terminal point from $\approx745$ to $350\hPa$, the signed
+integral becomes $226.6+0.01-287.5=-60.9\Jkg$ (at cycle phase $A$;
+$-65.90\Jkg$ at $x_d^-$), and the final clamp $\max(\,\cdot\,,0)$ maps it
+to $0$.  Composing the three steps,
+\[
+b\;\longmapsto\;\mathrm{INB}=\max\{j: b_j>0\}
+\;\longmapsto\; W(z_{\mathrm{top}})
+\;\longmapsto\; \max(\,\cdot\,,0),
+\]
+the first arrow is the discontinuous one; the clamp pins one side at
+exactly $0$ --- and, perhaps counter-intuitively, \emph{attenuates} the
+jump rather than enlarging it (\S\ref{sec:jumpmag}).  Note also
+(Fig.~\ref{fig:mech}b)
+that $W(350\hPa)=-61\Jkg<0$: even a ballistic parcel could never reach the
+level whose sign decides between $224$ and $0$ --- it would stall inside the
+valley.  Under \emph{every} available physical semantics ($\Emax$,
+$\Ereach$) the answer on both sides of the flicker is the same
+$\approx222\Jkg$.
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_pr77_mechanism.pdf}
+\caption{The issue-77 mechanism.  (a)~Eyewall-parcel buoyancy at the two
+cycle phases (curves nearly coincide); inset: the marginal excursion at
+$350\hPa$ whose sign differs between phases by $\pm14$~mK.  (b)~The
+corresponding work integrals: the argmax ($222\Jkg$ at $745\hPa$, green) is
+identical in both phases; the code's terminal point at phase $A$ (red cross)
+sits at $350\hPa$ where $W=-61\Jkg$, clamped to $0$.}
+\label{fig:mech}
+\end{figure}
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_pr77_capem.pdf}
+\caption{(a)~$E_M$ as a function of the trial pressure $x$: piecewise smooth
+with the cliff at $x_d$ (the structure at $x<948$ is a second, harmless
+flicker of the same kind).  (b)~The residual $F=g-x$ near $x_d$: the only
+sign change is across the gap; no root exists.}
+\label{fig:capem}
+\end{figure}
+
+\subsection{The magnitude of the jump, and the true role of the clamp}
+\label{sec:jumpmag}
+
+Several distinct magnitudes are in play, and they should not be conflated:
+
+\begin{center}
+\small
+\begin{tabular}{@{}lll@{}}
+\toprule
+quantity & value & what it is\\
+\midrule
+perturbation that flips the outcome & infinitesimal & a sign crossing
+($\pm8\times10^{-11}\K$ one-sided)\\
+the deciding feature itself & $\sim10^{-2}\Jkg$ & the marginal pocket
+whose sign is read\\
+jump of $E_M$ & $224.23\Jkg$ & the discontinuity, in the CAPE
+functional\\
+jump of $g$ & $0.838\hPa$ & the \emph{same} discontinuity, attenuated
+via \eqref{eq:g}\\
+\bottomrule
+\end{tabular}
+\end{center}
+
+\noindent
+The trigger being infinitesimal is not a separate scale --- it is what
+``discontinuity'' means: a finite response to an arbitrarily small
+perturbation.  The two quantitative facts worth explaining are that the
+finite response bears no relation to the deciding feature's own energy
+($224\Jkg$ decided by a $10^{-2}\Jkg$ pocket --- the selector's doing,
+\S\ref{sec:pr77}), and that the fixed-point map inherits the jump three
+orders of magnitude smaller ($0.838\hPa$ --- the thermodynamics of
+\eqref{eq:g}).  The attenuation is worth making exact, because the
+result's size relative to the stopping tolerance is what makes the failure
+silent.  Across $x_d$ every term of \eqref{eq:g} except $E_M$ is smooth, so
+the jump is $E_M$'s jump pushed through $g$.  Since the exponent moves by
+only $8.8\times10^{-4}$, linearization is essentially exact:
+\begin{equation}\label{eq:jumpmag}
+|\Delta g| \;\approx\; g\cdot\frac{\Delta E_M}{R_d\,\overline T_v}
+\,\Bigl(1-\tfrac12\,\kappa\,\tfrac{\mathrm{SST}}{T_O}\Bigr).
+\end{equation}
+Three factors attenuate.  \emph{Partial cancellation}:
+$c_1=1-\tfrac12\kappa\,\mathrm{SST}/T_O=0.342$, because $E_M$ enters
+$\mathrm{CAT}$ twice with opposite signs ($+E_M$ directly, $-0.658\,E_M$
+through the dissipative-heating term), so two-thirds of any $E_M$ change
+cancels before reaching the pressure.  \emph{The thermal-energy
+denominator}: $R_d\overline T_v=87.0\,\mathrm{kJ\,kg^{-1}}$ --- pressure
+responds by a factor $e$ only to that much $\mathrm{CAT}$, and $224\Jkg$ is
+$0.26\%$ of it.  \emph{The level itself}: $g\approx951\hPa$ multiplies the
+relative change.  Product: $\partial g/\partial E_M=3.74\times10^{-3}\hPa$
+per $\Jkg$; with $\Delta E_M=224.23\Jkg$, \eqref{eq:jumpmag} predicts
+$0.838\hPa$ against the measured $0.8381\hPa$ --- three significant
+figures.  The same formula with the 1995 column's own factors
+($\Delta E_M=203.0\Jkg$, $T_O=212.8\K$, $c_1=0.369$) predicts $0.840\hPa$
+against its measured $0.839\hPa$.
+
+That second agreement exposes the \emph{generic scale} of these failures:
+a marginal column's dominated lower hump is worth a couple hundred
+$\Jkg$, and a couple hundred $\Jkg$ times
+$\sim3.7\times10^{-3}\hPa/\Jkg$ is $0.7$--$0.9\hPa$ --- sitting just above
+the solver's $0.5\hPa$ tolerance.  A substantially larger jump would
+require implausibly large clamped CAPE; a jump below tolerance is silently
+absorbed as ``convergence'' (the iterates land within $0.5\hPa$ of each
+other across the gap and the loop exits normally).  The window
+$0.5$--$1\hPa$ is precisely the silent-failure regime, and both known
+failures sit in it.  This is not a coincidence to be shrugged at: the two
+columns reach $\approx0.84\hPa$ from \emph{different} ingredients
+($224\times0.342$ vs.\ $203\times0.369$), because neither factor varies
+much over warm marginal columns.
+
+Finally, the clamp is \emph{not} what makes the map discontinuous --- and
+it does not enlarge the jump.  Disabling $\max(\,\cdot\,,0)$, so that $E_M$
+flips $-65.90\leftrightarrow224.23\Jkg$ instead of
+$0\leftrightarrow224.23\Jkg$, leaves the discontinuity fully intact and
+\emph{increases} the measured jump to $1.0845\hPa$
+(\eqref{eq:jumpmag} with $\Delta E_M=290.13\Jkg$: $1.084$), because the raw
+signed value on the flicker side lies farther from $224$ than the clamped
+$0$ does.  The clamp's actual roles are cosmetic and attenuating: it pins
+one branch at exactly $0$ --- the clean $0\leftrightarrow224$ signature ---
+and shaves a quarter off the jump.  The discontinuity lives in the selector
+alone.
+
+\subsection{The second instance, and the control}
+
+The 1995 sounding of the repository's regression tests fails by the
+identical geometry (Fig.~\ref{fig:1995}): its $400\hPa$ sample flickers,
+$E_M$ flips $0\leftrightarrow203.05\Jkg$, the residual jumps
+$+0.615\to-0.225\hPa$, and the iteration cycles at
+$961.017\leftrightarrow961.671\hPa$.  An adjacent-day
+control profile has an ordinary continuous crossing and converges --- the
+failure is a property of marginal columns, not of the data source.  Both
+failing columns are \emph{marginal} in the precise sense that the ambient
+CAPE $E_A$ is tiny ($1.19\Jkg$ for \#77; exactly $0$ for both 1995 profiles
+under the released quadrature): selector topology matters most near such
+boundaries, where the clamp is live and humps are shallow.
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.48\linewidth]{fig_1995_map.pdf}
+\caption{The second known failing sounding (real 1995 profile): the same
+gap-straddling geometry and 2-cycle.}
+\label{fig:1995}
+\end{figure}
+
+\subsection{Why the upstream ``rescue'' is not a solution}
+
+A recent upstream patch detects the 2-cycle and collapses it to its
+midpoint ($\approx950.966\hPa$ here), exiting through the normal
+convergence path.  But the one-sided residuals at the jump are $+0.507$ and
+$-0.331\hPa$: \emph{there is no root there to approximate}.  The midpoint is
+a convention for returning a number, not a solution of $g(x)=x$; moreover
+the thermodynamic fields ultimately reported are evaluated at one cycle
+phase, an arbitrary parity-dependent choice worth $\sim1.2\ms$ in
+$V_{\max}$.  Cycle detection remains useful as telemetry and as a guard
+against burning 200 iterations, but it treats the symptom.  The disease is
+the discontinuity.
+
+%======================================================================
+\section{What is generic, what is unusual}
+\label{sec:generic}
+%======================================================================
+
+Two separate questions: how common are jumps of $g$, and how common is the
+\emph{fatal configuration} (diagonal crossing inside a jump)?
+
+We scanned $g$ at $0.1\hPa$ resolution over $x\in[900,1005]$ for 1{,}200
+randomly sampled converging columns of the repository's 2024 North-Atlantic
+sample (Fig.~\ref{fig:survey}a).  \textbf{Jumps are generic}: $239/1200$
+columns ($20\%$) carry at least one jump, with sizes up to $\sim1\hPa$ ---
+unsurprising, since the tropical upper troposphere hovers near moist
+neutrality ($b\approx0$ over deep layers), so marginal sign flickers under
+the $x$-dependent parcel moisture are everywhere.  \textbf{The fatal
+configuration is rare}: in zero of the 1{,}200 columns does the diagonal
+crossing fall inside a jump, consistent with the full 7{,}093-column sample
+containing no convergence failures.
+
+How close is ``close''?  Perturbing the Francine column's SST further
+(Fig.~\ref{fig:survey}b), the no-fixed-point window is $\sim0.04\K$ wide out
+of $6\K$ scanned.  This explains where the failure surfaced: an
+SST-attribution sweep around a real storm samples SST densely and will
+eventually land inside some column's window.  Single climatological runs
+rarely hit it; parameter sweeps over live cases --- where a silently missing
+value is most costly --- will.  The alignment is also less improbable than
+it looks: the iteration actively ferries orbits from both sides back across
+the jump, so once the crossing is inside the gap the cycle is
+\emph{attracting}, and a jump elsewhere passes unnoticed (either it does not
+straddle zero, or the induced cycle amplitude is below the $0.5\hPa$
+tolerance and is silently accepted as convergence).
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_survey.pdf}
+\caption{(a)~Largest jump of $g$ per column over 1{,}200 sampled real columns
+(only the 239 with a detected jump shown).  Jumps are common; none coincides
+with its column's diagonal crossing.  (b)~The Francine column under further
+SST perturbation: the fixed point moves smoothly except in a $\sim0.04\K$
+band (red) where it does not exist; the issue-77 member sits inside it.}
+\label{fig:survey}
+\end{figure}
+
+%======================================================================
+\section{To what extent is the discontinuity physical?}
+\label{sec:physical}
+%======================================================================
+
+Argmax-type constructions \emph{do} have genuine discontinuities, and one
+must not repair those away.  The precise accounting:
+
+\paragraph{The maximizer may genuinely jump; the value may not.}
+CAPE-as-energy is the \emph{value} of $\max_q W(q;\lambda)$; the LNB and
+outflow temperature are the \emph{maximizer}.  By the maximum theorem the
+value is continuous in the data; the maximizer is only upper hemicontinuous
+and jumps at dominance exchanges.  Figure~\ref{fig:bergefig} shows this on a
+real column: a $0.4\%$ moisture change teleports the argmax
+$870\to410\to180\hPa$ while the value crosses each exchange with a kink.
+Physically the maximizer's jump is a real bifurcation --- the preferred
+detrainment layer switches --- and any code reporting $T_O$ must live with
+it.  (In \eqref{eq:g}, $T_O$ comes from the \emph{saturated} parcel, whose
+$W$ is single-humped with a $\sim40\times$ dominance margin for
+storm-relevant SSTs, so this genuine discontinuity is dormant in practice;
+see \S\ref{sec:micro}.)  The maximizing level is in any case metadata rather
+than a prediction: where real eyewall air detrains depends on dynamics and
+mixing absent from a reversible parcel calculation.
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_berge.pdf}
+\caption{Berge's theorem on a real column (June, $28.5^{\circ}$C SST).
+(a)~The eyewall-parcel work integral has two competing humps $450\hPa$
+apart.  (b,c)~Sweeping boundary-layer moisture by $\pm6\%$: CAPE
+$=\max_qW$ is continuous through both dominance exchanges (kinks), while
+the argmax jumps.  Value and maximizer have different continuity classes;
+only the maximizer's jumps are physical.}
+\label{fig:bergefig}
+\end{figure}
+
+\paragraph{Reachable energy can genuinely jump --- in a different problem.}
+Under ballistic semantics (Convention~\ref{conv:reach}) the value itself
+jumps when an energy barrier closes at a tangency (toy sweep 2).  That is
+legitimate physics --- of \emph{spontaneous free convection}: an isolated
+parcel launched from rest does not teleport through a barrier.  It is not
+the physics of hurricane PI, whose eyewall ascent is forced by the storm's
+own frictionally convergent circulation \cite{garner2015,bisteremanuel2002}.
+A reachability-restricted CAPE answers a different question and should at
+most be a separately named diagnostic, never the PI default
+(\S\ref{sec:fix}).
+
+\paragraph{The code's discontinuity is neither of these.}
+Convention~\ref{conv:top} jumps whenever the topmost sign of $b$ flickers
+--- a strict superset of the physical loci above.  At the issue-77 event
+specifically: the deciding feature carries $O(10^{-2})\Jkg$ against a
+$224\Jkg$ decision (no energetic semantics); the selected terminal level has
+$W=-61\Jkg$, unreachable even ballistically (no reachability semantics);
+and both $\Emax$ and $\Ereach$ are constant through the event.  If the
+``topmost level'' rule encodes an intention --- ``outflow happens at the
+highest level convection can reach'' --- then $\Emax$'s argmax honors that
+intention on every profile where the top level is attainable and worth
+attaining; where the rules differ, $\Etop$ is charging the parcel an entry
+fee for a room it never enters.
+
+\begin{center}
+\small
+\begin{tabular}{@{}>{\raggedright\arraybackslash}p{0.30\linewidth}>{\raggedright\arraybackslash}p{0.24\linewidth}>{\raggedright\arraybackslash}p{0.36\linewidth}@{}}
+\toprule
+feature & status & reason \\
+\midrule
+multiple buoyant layers & physical possibility & stratification can create
+several neutral levels \\
+birth of a remote positive pocket & genuine mathematics & a smooth field can
+undergo a fold/tangency \\
+maximizing LNB jumping & genuine (variational) & work maxima exchange
+dominance; the common value stays continuous (Berge) \\
+reachable energy jumping & genuine, wrong problem & barrier closure is
+free-convection physics, not forced-ascent PI \\
+uppermost-positive-node jumping & genuine mathematics, wrong value rule &
+that node need not maximize work nor be reachable \\
+$E_M$ jumping $0\leftrightarrow224\Jkg$ & artifact & selector $+$ signed
+integration $+$ clamp amplify a one-bit sign \\
+pressure two-cycle & solver artifact & a discontinuous map is iterated as
+though contractive \\
+\bottomrule
+\end{tabular}
+\end{center}
+
+\subsection{What the cause is not}\label{sec:notcause}
+
+Three natural alternative hypotheses were tested and refuted, and deserve
+one paragraph so nobody re-litigates them.  \emph{Quadrature order}: the
+trapezoid rule with sub-grid terminal crossings integrates the
+piecewise-linear reconstruction exactly and handles simple crossings
+continuously; the jump survives with the same quadrature under
+Convention~\ref{conv:max}, and disappears with the same quadrature under a
+changed \emph{selection} --- the quadrature is innocent.  \emph{Grid
+resolution / interpolation smoothness}: rebuilding the Francine environment
+with a $C^1$ monotone-cubic interpolant of $T,r$ in $\ln p$ resampled at
+$5\hPa$ resolves the $350\hPa$ feature (it is real, and $4\times$ taller
+than the coarse grid shows) and happens to converge --- but the map still
+has an $0.83\hPa$ jump, merely relocated to $951.4\hPa$, off the diagonal
+\emph{by luck}; the same resampling with linear interpolation leaves the
+jump straddling the diagonal and still fails.  No reconstruction, of any
+smoothness class, can remove the discontinuity of $\Etop$: \S\ref{sec:toy}
+exhibits it for $C^\infty$ data.  \emph{Floating point}: the one-sided
+$10^{-11}\K$ values are ordinary continuity of $b_{18}(x)$ through zero, not
+roundoff; the feature itself is tens of millikelvin at the cycle phases.
+Discretization's real role is secondary: it quantizes \emph{where} flickers
+can fire (at nodes) and lends a $\sim100\hPa$-wide lever arm to a narrow
+feature, but the selection rule decides \emph{that} they fire.
+
+%======================================================================
+\section{The repair}
+\label{sec:fix}
+%======================================================================
+
+\subsection{Maximum-work CAPE: one changed selection}
+
+Replace Convention~\ref{conv:top} by Convention~\ref{conv:max} in all three
+CAPE evaluations of \eqref{eq:g}: accumulate the \emph{same} running signed
+work $W$ over the \emph{same} piecewise-linear reconstruction with the
+\emph{same} sub-grid crossing candidates, and take the maximum over
+candidate terminal points instead of the value at the topmost crossing.
+Return the argmax as the LNB (with an explicit tie rule).  No clamp is
+needed.  Nothing is skipped and nothing is free: every negative-buoyancy
+layer reduces $W$ along the way.  The change is one selection, not new
+machinery --- and for single-hump profiles it reproduces the legacy value
+exactly.
+
+This is also the physically correct convention for PI
+(\S\ref{sec:physical}): hurricane CAPE is thermodynamic work around a
+forced cycle, so the outflow level is the one maximizing extractable work,
+with traversal costs charged but no ballistic reachability requirement.
+The ballistic (first-connected-lobe) rule was also implemented as a
+proof-of-mechanism: it too removes the failures, but it answers the
+free-convection question ---
+and it perturbs the healthy control column by $1.2\hPa$ where max-work
+moves it by $0.0035\hPa$.  Maximum signed work is the minimal-change,
+correct-semantics repair.
+
+\subsection{Existence, and solving for a root}
+
+\begin{proposition}\label{prop:cont}
+With $\Emax$ in place of $\Etop$ in all three CAPE evaluations, and modulo
+the micro-jumps of \S\ref{sec:micro}, the map $g$ of \eqref{eq:g} is
+continuous on $[400,\mathrm{MSL}]$ wherever the saturated parcel's argmax is
+unique, and $g\le\mathrm{MSL}$ so that $F(x)=g(x)-x$ satisfies
+$F(\mathrm{MSL})\le0$.  If in addition the column satisfies the
+non-hypercane endpoint condition $g(400)>400$, then $F(400)>0$ and the
+intermediate value theorem gives a fixed point in $(400,\mathrm{MSL}]$.
+\end{proposition}
+
+\begin{proof}[Proof sketch]
+Each $b_j$ is continuous in $x$ (Observation~\ref{obs:cont}); $W$ is jointly
+continuous on a compact set; $\Emax=\max_qW$ is continuous (maximum
+theorem); $E_A$ is constant and $\overline T_v$ smooth and positive;
+$T_O(x)$ is continuous wherever the saturated argmax is unique;
+$\max(\,\cdot\,,0)$ and $\exp$ preserve continuity.  $\mathrm{CAT}\ge0$
+gives $g\le\mathrm{MSL}$.
+\end{proof}
+
+Both hypotheses are explicit for a reason.  If $g(400)\le400$ the physical
+floor has been reached and a distinct domain flag should be returned, not a
+root.  And because residual non-smoothness survives at the $10^{-2}\hPa$
+scale (\S\ref{sec:micro}), a bracketed solver (bisection or safeguarded
+Brent on $F$) should declare success only after an explicit residual check
+$|F(x_*)|\le\varepsilon_F$ --- bisection on a residual with an undetected
+jump converges to the jump, exactly the trap of Observation~\ref{obs:pr77}.
+With those guards, bracketing terminates deterministically in $\sim40$
+evaluations, needs no contractivity assumption, and cannot enter a
+period-two orbit.  (The legacy Picard loop also converges here, since
+$|g'|\le0.47$ off the former jump --- but bracketing turns an empirical
+slope into a guarantee.)
+
+\subsection{Results}
+
+\begin{itemize}
+\item \textbf{Failing cases (exact roots).}  \#77/Francine: $x_*=950.3136889226\hPa$,
+$V_{\max}=83.9065\ms$, $P_{\min}=920.7517\hPa$, $T_O=206.088\K$.  1995
+marginal: $x_*=960.7687\hPa$, $V_{\max}=74.4149\ms$.  The legacy Picard loop
+with max-work CAPE also satisfies its operational $0.5\hPa$ rule after five
+updates; note its fifth iterate ($950.4115$, residual $-0.064\hPa$) is the
+operational answer, not the root --- the values above come from bracketing.
+For the 1995 case the corrected answer sits $1.1\ms$ from what the upstream
+rescue reports ($75.47$, the physics of an arbitrary cycle phase).
+\item \textbf{Benchmark population} (patched compiled kernel, rescue
+deleted; 7{,}093 converging columns of the 2024 sample): \textbf{zero} flag
+changes; $39\%$ of $V_{\max}$ values bit-identical; $98.1\%$ within
+$0.1\ms$ (the sub-$10^{-3}$ bulk is float summation-order noise laundered
+through the $0.5\hPa$ stopping rule, Fig.~\ref{fig:fixfig}b); a $1.9\%$ tail
+(138 columns) moves by up to $7.9\ms$.  Fixed-pressure dissection shows
+every tail column is a dominated-hump or near-surface-sliver configuration
+--- precisely the profiles on which the legacy convention was fragile by
+construction (a millikelvin perturbation could move its answer by hundreds
+of $\Jkg$).
+\item \textbf{Reference parity and process.}  Agreement with the MATLAB
+\texttt{pcmin} reference is untouched on the $98\%$ and diverges only on
+the fragile tail --- necessarily, since the discontinuity is \emph{in} the
+reference convention.  Maximum-work CAPE is an intentional, argued
+scientific convention change: it should be reviewed upstream, versioned,
+and (if compatibility demands) initially exposed behind an explicit
+\texttt{cape\_definition} option.  Bit-identity on every column is not an
+appropriate acceptance criterion for an intentional correction; identity of
+flags, boundedness of the audited tail, and the residual check are.
+\end{itemize}
+
+\begin{figure}[tb]
+\centering
+\includegraphics[width=0.98\linewidth]{fig_fix.pdf}
+\caption{The repair: the same map \eqref{eq:g} evaluated with the two CAPE
+conventions on the Francine column.  (a)~With the legacy $\Etop$ (red,
+solid) the map has the gap and no fixed point.  With $\Emax$ (teal, dashed)
+it is continuous, crosses the diagonal, and the iteration converges (dot).
+On the right branch the two coincide to line width --- red shows through
+the dash gaps --- because there the topmost-positive level already selects
+the dominant lower cell; on the left branch they differ by the full effect
+of the $222\Jkg$ clamp.  (b)~Cumulative distribution of
+$|\Delta V_{\max}|$ between conventions over the benchmark: $39\%$
+bit-identical (left edge, clipped to $10^{-12}$), $98.1\%$ below $0.1\ms$,
+and a $1.9\%$ tail of discontinuity-adjacent columns.}
+\label{fig:fixfig}
+\end{figure}
+
+\subsection{Residual non-smoothness}\label{sec:micro}
+
+What remains after the repair, all quantified and benign:
+\begin{enumerate}
+\item \emph{Newton stopping rule} in the entropy inversion (tolerance
+$10^{-3}\K$): the returned temperature is a lagged iterate, so $b_j$ has
+micro-jumps up to $O(10^{-3}\K)$; measured step-to-step roughness of the
+repaired map is $\le7\times10^{-3}\hPa$ --- $50\times$ below the solver
+tolerance, and covered by the residual check.  A literal continuity theorem
+applies to the converged thermodynamic relation, not to the last bit of a
+stopped Newton iterate.
+\item \emph{Condensation-level formula switch} at grid nodes: same order.
+\item \emph{Saturated-parcel argmax exchange} (the genuine Berge jump of
+$T_O$): dormant for storm-relevant profiles (dominance margin
+$\sim40\times$), but a robust API should detect ties/near-ties in the
+saturated argmax and return a branch-ambiguity flag rather than silently
+jumping; if it ever fires, it is physics, and residual-checked bracketing
+still terminates.
+\end{enumerate}
+
+\subsection{Implementation and validation plan}
+
+Staged so numerical and scientific changes are separable:
+
+\begin{description}
+\item[A --- diagnostics, no changed outputs.]  Expose a scalar evaluator for
+$F(x)$ and CAPE diagnostics (all buoyancy roots, selected root, cumulative
+work at each candidate, one-sided values); add automatic jump and
+period-$k$ scans.
+\item[B --- max-work kernel.]  Implement $\Emax$ in a transparent scalar
+reference, then in the compiled kernel without arithmetic reassociation;
+return value plus argmax metadata; classify every material difference from
+legacy on the full sample before choosing the default/compatibility policy.
+\item[C --- outflow semantics.]  Specify the tie rule and the relation
+between the maximizing level and PI's outflow temperature; detect saturated
+argmax exchanges; flag LNBs clipped by the $p_{\mathrm{top}}$ truncation
+(Remark~\ref{rem:ptop}).  If ballistic reachability is wanted for other
+applications, expose it as a distinctly named free-convection diagnostic,
+never as the PI default.
+\item[D --- bracketed solver.]  Brent/bisection with cached CAPE
+evaluations, endpoint-sign check first, distinct flags for invalid input,
+entropy-solver failure, no-bracket/physical floor, and branch ambiguity;
+success requires $|F(x_*)|\le\varepsilon_F$.
+\end{description}
+
+Required regression families: the exact \#77 column with one-sided values
+pinned; the 1995 marginal and control pair; perturbations of SST, surface
+humidity, MSL and each temperature level around both marginal cases;
+analytic one-lobe/two-lobe/tangency profiles with known roots; and the full
+sample with output deltas \emph{and} LNB-branch-identity changes reported.
+On every successful column, verify $|F(x_*)|\le\varepsilon_F$, that the
+returned CAPE equals the maximum candidate work, and that the LNB obeys the
+declared tie convention.
+
+%======================================================================
+\section{Final interpretation}
+%======================================================================
+
+The causal chain, in full:
+\begin{align*}
+\text{marginal remote buoyancy pocket}
+&\;\longrightarrow\; \text{hard highest-positive-node selector}\\
+&\;\longrightarrow\; \text{finite CAPE jump (helped by the clamp)}\\
+&\;\longrightarrow\; \text{pressure-map jump straddling the diagonal}\\
+&\;\longrightarrow\; \text{attracting, bit-exact 2-cycle}\\
+&\;\longrightarrow\; \text{``missing'' PI for a real hurricane's
+counterfactual.}
+\end{align*}
+What is physical: multiple neutral levels, dominance exchanges between work
+maxima (maximizer jumps, value continuous), and barrier transitions in the
+free-convection problem.  What is conventional: calling the uppermost root
+``the'' LNB.  What is artificial: letting a marginal, energetically
+negligible, dynamically unreachable, effectively one-node pocket erase
+$\sim224\Jkg$ of lower-tropospheric work and thereby delete the fixed point
+of the pressure equation.  The repair --- maximum signed work as the CAPE
+value, argmax as reported metadata, residual-checked bracketing as the
+solver --- is one changed selection and one changed loop, restores an
+existence proof under explicit hypotheses, and makes the two-cycle chase
+disappear for the right reason.
+
+\appendix
+
+%======================================================================
+\section{Key numbers for the failing sounding}
+%======================================================================
+
+\begin{center}
+\footnotesize
+\begin{tabular}{@{}ll@{}}
+\toprule
+quantity & value\\
+\midrule
+provenance & Francine env., 18 UTC 2024-09-11, $28.6^{\circ}$N
+$92.1^{\circ}$W, $\Delta\mathrm{SST}=-1.6\K$\\
+SST, MSL & $28.20263672^{\circ}$C, $1014.96545\hPa$\\
+levels (retained $>50\hPa$) & 37 (28)\\
+cycle $A$, $B$ & $950.6533454438$, $951.2790079839\hPa$
+($B{-}A=0.6257$)\\
+jump location, size & $x_d\simeq950.8289947515\hPa$, $0.837\hPa$ downward\\
+one-sided residuals at $x_d$ & $+0.5073447$ / $-0.3307594\hPa$\\
+slope of $g$ off the jump & $0.34$--$0.46$ (on $[946,958]$)\\
+$b$ at $350\hPa$, cycle phases & $+0.0137$ / $-0.0350\K$\\
+$b$ at $350\hPa$, at $x_d^{\mp}$ & $+7.6\times10^{-11}$ /
+$-8.0\times10^{-11}\K$\\
+$E_M$ across $x_d$ & $0 \leftrightarrow 224.23146\Jkg$\\
+raw signed integral before clamp & $-65.90\Jkg$ at $x_d^-$
+($-60.86$ at phase $A$)\\
+pressure sensitivity $\partial g/\partial E_M$ &
+$3.74\times10^{-3}\hPa/\Jkg$ ($=g\,c_1/(R_d\overline T_v)$, $c_1=0.342$)\\
+jump with clamp disabled & $1.0845\hPa$ (measured; predicted $1.084$)\\
+$E_{MS}$, $T_O$ at $x_d$ (smooth) & $8565.00966\Jkg$, $206.1096884\K$\\
+ambient CAPE $E_A$ & $1.19\Jkg$ (\#77); $0$ (both 1995 profiles)\\
+corrected root (max-work, bracketed) & $x_*=950.3136889226\hPa$\\
+corrected outputs & $V_{\max}=83.9065\ms$, $P_{\min}=920.7517\hPa$\\
+1995 marginal: cycle & $961.017\leftrightarrow961.671\hPa$\\
+1995 marginal: corrected root & $x_*=960.7687\hPa$, $V_{\max}=74.4149\ms$\\
+1995 control (max-work vs legacy) & fixed point moves $0.0035\hPa$
+(ballistic rule: $1.2\hPa$)\\
+\bottomrule
+\end{tabular}
+\end{center}
+
+%======================================================================
+\section{Reproduction}\label{app:repro}
+%======================================================================
+
+All computations were run against branch \texttt{upstream/modernize\_eg}
+(commit \texttt{abc4b4e}) of the \texttt{tcpyPI} repository, using the
+committed 2024 ERA5 North-Atlantic sample and its committed outputs.  The
+analysis scripts live in \texttt{discontinuity\_analysis/scripts/}:
+instrumented map evaluation (\texttt{gmap.py}), fine scans
+(\texttt{scan.py}), the max-work prototype and full-sample validation
+(\texttt{maxw.py}, \texttt{validate.py}, \texttt{compare.py},
+\texttt{dissect.py}), the survey and SST sweep (\texttt{datagen.py}), the
+dominance-exchange hunt (\texttt{humphunt.py}), the reconstruction
+refutation (\texttt{refine.py}), and figure generators (\texttt{plots.py},
+\texttt{plots2.py}).  Every headline quantity (cycle values, jump location,
+one-sided values, corrected roots, $V_{\max}$, $P_{\min}$, benchmark
+statistics) was additionally reproduced by a second, independently written
+implementation of the map and CAPE routines, and the two agree to the
+digits shown.
+The \#77 profile is the one committed in PR~\#77's failing test; the 1995
+profiles are those of \texttt{tests/test\_oscillation\_rescue.py}.
+
+\begin{thebibliography}{9}
+
+\bibitem{bisteremanuel2002}
+M. Bister and K. A. Emanuel,
+``Low frequency variability of tropical cyclone potential intensity. 1.
+Interannual to interdecadal variability,''
+\emph{Journal of Geophysical Research}, 107(D24), 4801, 2002.
+\url{https://doi.org/10.1029/2001JD000776}.
+
+\bibitem{gilford2021}
+D. M. Gilford,
+``pyPI (v1.3): Tropical Cyclone Potential Intensity Calculations in Python,''
+\emph{Geoscientific Model Development}, 14, 2351--2369, 2021.
+\url{https://doi.org/10.5194/gmd-14-2351-2021}.
+
+\bibitem{garner2015}
+S. T. Garner,
+``The Relationship between Hurricane Potential Intensity and CAPE,''
+\emph{Journal of the Atmospheric Sciences}, 72, 141--163, 2015.
+\url{https://doi.org/10.1175/JAS-D-14-0008.1}.
+
+\bibitem{emanuel1994}
+K. A. Emanuel, \emph{Atmospheric Convection}, Oxford University Press, 1994.
+
+\bibitem{smith2022}
+R. K. Smith,
+``Effective buoyancy and CAPE: Some implications for tropical cyclones,''
+\emph{Quarterly Journal of the Royal Meteorological Society}, 2022.
+\url{https://doi.org/10.1002/qj.4294}.
+
+\bibitem{pr77}
+tcpyPI pull request \#77, ``Add failing test to illustrate failure of
+convergence,'' including the exact regression column.
+\url{https://github.com/dgilford/tcpyPI/pull/77}.
+
+\bibitem{francine}
+National Hurricane Center,
+\emph{Tropical Cyclone Report: Hurricane Francine (AL062024), 9--12
+September 2024}, best-track Table 1.
+\url{https://www.nhc.noaa.gov/data/tcr/AL062024_Francine.pdf}.
+
+\end{thebibliography}
+
+\end{document}
diff --git a/discontinuity_analysis/figures/fig_1995_map.pdf b/discontinuity_analysis/figures/fig_1995_map.pdf
new file mode 100644
index 0000000..55a18ec
Binary files /dev/null and b/discontinuity_analysis/figures/fig_1995_map.pdf differ
diff --git a/discontinuity_analysis/figures/fig_1995_map.png b/discontinuity_analysis/figures/fig_1995_map.png
new file mode 100644
index 0000000..21ac3fa
Binary files /dev/null and b/discontinuity_analysis/figures/fig_1995_map.png differ
diff --git a/discontinuity_analysis/figures/fig_berge.pdf b/discontinuity_analysis/figures/fig_berge.pdf
new file mode 100644
index 0000000..95727ba
Binary files /dev/null and b/discontinuity_analysis/figures/fig_berge.pdf differ
diff --git a/discontinuity_analysis/figures/fig_berge.png b/discontinuity_analysis/figures/fig_berge.png
new file mode 100644
index 0000000..55a323a
Binary files /dev/null and b/discontinuity_analysis/figures/fig_berge.png differ
diff --git a/discontinuity_analysis/figures/fig_conventions.pdf b/discontinuity_analysis/figures/fig_conventions.pdf
new file mode 100644
index 0000000..cf9b2e2
Binary files /dev/null and b/discontinuity_analysis/figures/fig_conventions.pdf differ
diff --git a/discontinuity_analysis/figures/fig_conventions.png b/discontinuity_analysis/figures/fig_conventions.png
new file mode 100644
index 0000000..d91266d
Binary files /dev/null and b/discontinuity_analysis/figures/fig_conventions.png differ
diff --git a/discontinuity_analysis/figures/fig_discretization.pdf b/discontinuity_analysis/figures/fig_discretization.pdf
new file mode 100644
index 0000000..f4abe48
Binary files /dev/null and b/discontinuity_analysis/figures/fig_discretization.pdf differ
diff --git a/discontinuity_analysis/figures/fig_discretization.png b/discontinuity_analysis/figures/fig_discretization.png
new file mode 100644
index 0000000..638879e
Binary files /dev/null and b/discontinuity_analysis/figures/fig_discretization.png differ
diff --git a/discontinuity_analysis/figures/fig_fix.pdf b/discontinuity_analysis/figures/fig_fix.pdf
new file mode 100644
index 0000000..6f1f811
Binary files /dev/null and b/discontinuity_analysis/figures/fig_fix.pdf differ
diff --git a/discontinuity_analysis/figures/fig_fix.png b/discontinuity_analysis/figures/fig_fix.png
new file mode 100644
index 0000000..9aef73a
Binary files /dev/null and b/discontinuity_analysis/figures/fig_fix.png differ
diff --git a/discontinuity_analysis/figures/fig_healthy_anatomy.pdf b/discontinuity_analysis/figures/fig_healthy_anatomy.pdf
new file mode 100644
index 0000000..358cd65
Binary files /dev/null and b/discontinuity_analysis/figures/fig_healthy_anatomy.pdf differ
diff --git a/discontinuity_analysis/figures/fig_healthy_anatomy.png b/discontinuity_analysis/figures/fig_healthy_anatomy.png
new file mode 100644
index 0000000..a913cb1
Binary files /dev/null and b/discontinuity_analysis/figures/fig_healthy_anatomy.png differ
diff --git a/discontinuity_analysis/figures/fig_healthy_map.pdf b/discontinuity_analysis/figures/fig_healthy_map.pdf
new file mode 100644
index 0000000..f4f644b
Binary files /dev/null and b/discontinuity_analysis/figures/fig_healthy_map.pdf differ
diff --git a/discontinuity_analysis/figures/fig_healthy_map.png b/discontinuity_analysis/figures/fig_healthy_map.png
new file mode 100644
index 0000000..6e24224
Binary files /dev/null and b/discontinuity_analysis/figures/fig_healthy_map.png differ
diff --git a/discontinuity_analysis/figures/fig_pr77_capem.pdf b/discontinuity_analysis/figures/fig_pr77_capem.pdf
new file mode 100644
index 0000000..873284a
Binary files /dev/null and b/discontinuity_analysis/figures/fig_pr77_capem.pdf differ
diff --git a/discontinuity_analysis/figures/fig_pr77_capem.png b/discontinuity_analysis/figures/fig_pr77_capem.png
new file mode 100644
index 0000000..3ee03b9
Binary files /dev/null and b/discontinuity_analysis/figures/fig_pr77_capem.png differ
diff --git a/discontinuity_analysis/figures/fig_pr77_map.pdf b/discontinuity_analysis/figures/fig_pr77_map.pdf
new file mode 100644
index 0000000..cfb3b8a
Binary files /dev/null and b/discontinuity_analysis/figures/fig_pr77_map.pdf differ
diff --git a/discontinuity_analysis/figures/fig_pr77_map.png b/discontinuity_analysis/figures/fig_pr77_map.png
new file mode 100644
index 0000000..401bc3c
Binary files /dev/null and b/discontinuity_analysis/figures/fig_pr77_map.png differ
diff --git a/discontinuity_analysis/figures/fig_pr77_mechanism.pdf b/discontinuity_analysis/figures/fig_pr77_mechanism.pdf
new file mode 100644
index 0000000..45ff2d5
Binary files /dev/null and b/discontinuity_analysis/figures/fig_pr77_mechanism.pdf differ
diff --git a/discontinuity_analysis/figures/fig_pr77_mechanism.png b/discontinuity_analysis/figures/fig_pr77_mechanism.png
new file mode 100644
index 0000000..fbb8526
Binary files /dev/null and b/discontinuity_analysis/figures/fig_pr77_mechanism.png differ
diff --git a/discontinuity_analysis/figures/fig_survey.pdf b/discontinuity_analysis/figures/fig_survey.pdf
new file mode 100644
index 0000000..fc9d044
Binary files /dev/null and b/discontinuity_analysis/figures/fig_survey.pdf differ
diff --git a/discontinuity_analysis/figures/fig_survey.png b/discontinuity_analysis/figures/fig_survey.png
new file mode 100644
index 0000000..6f22843
Binary files /dev/null and b/discontinuity_analysis/figures/fig_survey.png differ
diff --git a/discontinuity_analysis/figures/fig_toy.pdf b/discontinuity_analysis/figures/fig_toy.pdf
new file mode 100644
index 0000000..3e4bd78
Binary files /dev/null and b/discontinuity_analysis/figures/fig_toy.pdf differ
diff --git a/discontinuity_analysis/figures/fig_toy.png b/discontinuity_analysis/figures/fig_toy.png
new file mode 100644
index 0000000..a54f069
Binary files /dev/null and b/discontinuity_analysis/figures/fig_toy.png differ
diff --git a/discontinuity_analysis/figures/fig_toymap.pdf b/discontinuity_analysis/figures/fig_toymap.pdf
new file mode 100644
index 0000000..385cb1c
Binary files /dev/null and b/discontinuity_analysis/figures/fig_toymap.pdf differ
diff --git a/discontinuity_analysis/figures/fig_toymap.png b/discontinuity_analysis/figures/fig_toymap.png
new file mode 100644
index 0000000..15f5a49
Binary files /dev/null and b/discontinuity_analysis/figures/fig_toymap.png differ
diff --git a/discontinuity_analysis/scripts/caseC2.py b/discontinuity_analysis/scripts/caseC2.py
new file mode 100644
index 0000000..ae7cc53
--- /dev/null
+++ b/discontinuity_analysis/scripts/caseC2.py
@@ -0,0 +1,55 @@
+"""Regenerate the 'healthy column' case with a representative interior-LNB
+column (m=0, j=20, i=1). Saves healthy2.npz."""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.join(SCRATCH, "branch", "src"))
+sys.path.insert(0, SCRATCH)
+
+import numpy as np
+import xarray as xr
+from tcpyPI import utilities
+
+from datagen import g_scan, iterate_raw, W_curve
+
+M, J, I = 0, 20, 1
+
+ds = xr.open_dataset(os.path.join(SCRATCH, "sample_data.nc"))
+p_s = ds["p"].values.astype(float)
+sst = float(ds["sst"].values[M, J, I])
+msl = float(ds["msl"].values[M, J, I])
+T = utilities.T_Ctok(ds["t"].values[M, :, J, I].astype(float))
+R = ds["r"].values[M, :, J, I].astype(float) * 0.001
+R[np.isnan(R)] = 0.0
+lat = float(ds["lat"].values[J])
+lon = float(ds["lon"].values[I])
+month = int(ds["month"].values[M]) if "month" in ds else M + 1
+print(f"column: month={month} lat={lat} lon={lon} SST={sst:.2f}")
+
+SSTK = utilities.T_Ctok(sst)
+ES0 = utilities.es_cc(sst)
+
+save = {}
+pms = np.arange(850.0, 1005.0, 0.02)
+sc = g_scan(pms, SSTK, msl, p_s, T, R, ES0)
+save["c_pm"] = pms
+save["c_g"] = sc[:, 0]
+st, xs = iterate_raw(SSTK, msl, p_s, T, R, ES0, nmax=60)
+save["c_iter"] = xs
+print(f"iteration: status={st}, fixed point ~ {xs[-1]:.4f}, n={len(xs)}")
+
+pmc = xs[-1]
+PP = min(pmc, 1000.0)
+pts, Ws, b, Pt = W_curve(SSTK, utilities.rv(ES0, PP), PP, T, R, p_s)
+save["c_W_p"] = pts
+save["c_W_W"] = Ws
+save["c_b"] = b
+save["c_blevels"] = p_s[: int(np.count_nonzero(p_s > 50))]
+kmax = int(np.argmax(Ws))
+print(f"argmax LNB: {pts[kmax]:.2f} hPa, CAPE={Ws[kmax]:.1f}; "
+      f"b sign change: b[{save['c_blevels'][-2]:.0f}]={b[-2]:+.2f}, "
+      f"b[{save['c_blevels'][-1]:.0f}]={b[-1]:+.2f}")
+np.savez(os.path.join(SCRATCH, "healthy2.npz"), **save)
+print("saved healthy2.npz")
diff --git a/discontinuity_analysis/scripts/compare.py b/discontinuity_analysis/scripts/compare.py
new file mode 100644
index 0000000..c8f36f1
--- /dev/null
+++ b/discontinuity_analysis/scripts/compare.py
@@ -0,0 +1,52 @@
+"""Compare branch (pcmin CAPE + rescue) vs branch_fix (max-W CAPE, no rescue)."""
+
+import numpy as np
+import xarray as xr
+
+a = np.load("branch_pi.npz")      # baseline: branch as merged
+b = np.load("branch_fix_pi.npz")  # max-W CAPE, rescue deleted
+
+both1 = (a["IFL"] == 1) & (b["IFL"] == 1)
+n_both = int(both1.sum())
+
+# flag agreement
+flag_diff = np.argwhere(a["IFL"] != b["IFL"])
+print(f"columns with IFL=1 in both: {n_both}")
+print(f"columns where IFL differs: {len(flag_diff)}")
+for idx in flag_diff[:20]:
+    m, j, i = idx
+    print(f"  (m={m},j={j},i={i}): IFL {a['IFL'][m,j,i]} -> {b['IFL'][m,j,i]}, "
+          f"VMAX {a['VMAX'][m,j,i]:.4f} -> {b['VMAX'][m,j,i]:.4f}")
+
+for var in ["VMAX", "PMIN", "TO", "OTL"]:
+    x, y = a[var], b[var]
+    d = np.abs(x - y)[both1]
+    ident = int((d == 0).sum())
+    print(f"{var}: bit-identical {ident}/{n_both} "
+          f"({100*ident/n_both:.2f}%), max|d|={np.nanmax(d):.6g}, "
+          f"n(|d|>1e-8)={int((d>1e-8).sum())}, n(|d|>0.1)={int((d>0.1).sum())}")
+
+# which columns changed materially, and were they rescue cases?
+d_vmax = np.abs(a["VMAX"] - b["VMAX"])
+chg = np.argwhere(both1 & (d_vmax > 1e-8))
+print(f"\ncolumns with VMAX change > 1e-8: {len(chg)}")
+for idx in chg[:30]:
+    m, j, i = idx
+    print(f"  (m={m},j={j},i={i}): VMAX {a['VMAX'][m,j,i]:.6f}->{b['VMAX'][m,j,i]:.6f} "
+          f"TO {a['TO'][m,j,i]:.3f}->{b['TO'][m,j,i]:.3f} "
+          f"OTL {a['OTL'][m,j,i]:.2f}->{b['OTL'][m,j,i]:.2f}")
+
+# MATLAB reference comparison for both
+ml = xr.open_dataset("matlab_ref.nc")
+print("\nMATLAB ref vars:", list(ml.data_vars))
+mv = ml["Vmax"].values
+mp = ml["Pmin"].values
+mfl = ml["PI_flag"].values if "PI_flag" in ml else None
+print("MATLAB shapes:", mv.shape)
+
+for name, r in [("branch", a), ("fix", b)]:
+    ok = (r["IFL"] == 1) & np.isfinite(mv) & (mfl == 1 if mfl is not None else True)
+    dv = np.abs(r["VMAX"] - mv)[ok]
+    dp = np.abs(r["PMIN"] - mp)[ok]
+    print(f"{name}: vs MATLAB over {int(ok.sum())} cols: "
+          f"VMAX max|d|={np.nanmax(dv):.6f} PMIN max|d|={np.nanmax(dp):.6f}")
diff --git a/discontinuity_analysis/scripts/datagen.py b/discontinuity_analysis/scripts/datagen.py
new file mode 100644
index 0000000..359032c
--- /dev/null
+++ b/discontinuity_analysis/scripts/datagen.py
@@ -0,0 +1,271 @@
+"""Generate all data for the discontinuity analysis document.
+
+Uses the COMPILED cape() from the branch tree for speed (survey, sweeps), and
+the instrumented pure-python cape for per-case internals (TVRDIF, W curves).
+
+Outputs: docdata.npz in this directory.
+"""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.join(SCRATCH, "branch", "src"))
+sys.path.insert(0, SCRATCH)
+
+import numpy as np
+import xarray as xr
+from tcpyPI import constants, utilities
+from tcpyPI.pi import cape as cape_jit  # compiled pcmin-convention cape
+
+from gmap import cape_instr  # instrumented pure-python (same convention)
+
+EPS = constants.EPS
+RD = constants.RD
+
+
+def g_eval(PM, SSTK, MSL, P, T, R, ES0, CKCD=0.9, diss_flag=1):
+    """One iteration body PM -> PNEW using compiled cape. Returns PNEW & parts."""
+    PP = min(PM, 1000.0)
+    RP = EPS * R[0] * MSL / (PP * (EPS + R[0]) - R[0] * MSL)
+    CAPEM, _, _, flM = cape_jit(T[0], RP, PP, T, R, P, 0, 50, 1)
+    RPS = utilities.rv(ES0, PP)
+    CAPEMS, TOMS, LNBS, flS = cape_jit(SSTK, RPS, PP, T, R, P, 0, 50, 1)
+    TO = TOMS
+    RAT = SSTK / TO if diss_flag else 1.0
+    TV0 = utilities.Trho(T[0], R[0], R[0])
+    TVSST = utilities.Trho(SSTK, RPS, RPS)
+    TVAV = 0.5 * (TV0 + TVSST)
+    CAPEA, _, _, _ = cape_jit(T[0], R[0], P[0], T, R, P, 0, 50, 1)
+    CAT = max((CAPEM - CAPEA) + 0.5 * CKCD * RAT * (CAPEMS - CAPEM), 0.0)
+    return MSL * np.exp(-CAT / (RD * TVAV)), CAPEM, CAPEMS, TO
+
+
+def g_scan(pms, SSTK, MSL, P, T, R, ES0):
+    out = np.empty((len(pms), 4))
+    for k, pm in enumerate(pms):
+        out[k] = g_eval(pm, SSTK, MSL, P, T, R, ES0)
+    return out
+
+
+def iterate_raw(SSTK, MSL, P, T, R, ES0, x0=970.0, nmax=200, tol=0.5):
+    """Raw pcmin iteration without rescue. Returns (status, iterates).
+    status: 1 converged, 2 cycle/cap, 0 blowdown(PM<400)."""
+    PM, PMOLD, PNEW, NP = x0, x0, 0.0, 0
+    xs = []
+    while abs(PNEW - PMOLD) > tol:
+        PNEW = g_eval(PM, SSTK, MSL, P, T, R, ES0)[0]
+        xs.append(PM)
+        PMOLD, PM = PM, PNEW
+        NP += 1
+        if NP > nmax or PM < 400:
+            return (2 if PM >= 400 else 0), np.array(xs + [PM])
+    xs.append(PM)
+    return 1, np.array(xs)
+
+
+def W_curve(TP, RP, PP, T, R, P, ptop=50):
+    """Running signed work integral W(p_t): nodes + crossing candidates.
+    Returns arrays (p_t values, W values) densified, plus candidate list."""
+    base = cape_instr(TP, RP, PP, T, R, P, 0, ptop)
+    N = int(np.count_nonzero(P > ptop))
+    Pt, Tt = P[:N], T[:N]
+    b = base["TVRDIF"]
+    pts = [PP]
+    Ws = [0.0]
+    W = RD * (PP - Pt[0]) / (PP + Pt[0]) * b[0]
+    pts.append(Pt[0])
+    Ws.append(W)
+    for j in range(1, N):
+        b0, b1 = b[j - 1], b[j]
+        if b0 * b1 < 0.0:
+            Pc = (Pt[j] * b0 - Pt[j - 1] * b1) / (b0 - b1)
+            Wc = W + RD * b0 * (Pt[j - 1] - Pc) / (Pt[j - 1] + Pc)
+            pts.append(Pc)
+            Ws.append(Wc)
+        W += RD * (b1 + b0) * (Pt[j - 1] - Pt[j]) / (Pt[j] + Pt[j - 1])
+        pts.append(Pt[j])
+        Ws.append(W)
+    return np.array(pts), np.array(Ws), b, Pt
+
+
+# ---------------- profiles ----------------
+P77 = np.array([1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650,
+                600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125,
+                100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1], dtype=float)
+TC77 = np.array([25.260956, 23.078949, 20.881561, 18.652649, 16.615143, 14.63147,
+                 12.794586, 11.79306, 11.01236, 10.847565, 10.354492, 7.807007,
+                 5.473297, 2.5278625, -1.532135, -7.3142395, -13.635345,
+                 -20.613602, -28.588928, -37.270096, -45.40825, -50.047455,
+                 -54.84575, -59.173737, -62.878662, -65.78009, -69.30669,
+                 -65.512024, -60.76361, -55.100464, -51.38333, -43.95462,
+                 -39.709717, -33.193268, -21.576843, -15.402496, -13.186676])
+R77 = np.array([1.0783693e01, 1.0704287e01, 1.0680210e01, 1.0617845e01,
+                1.0320683e01, 9.8811483e00, 9.1884289e00, 7.1884680e00,
+                5.6963191e00, 3.5568204e00, 1.5912720e00, 1.0433695e00,
+                5.9723043e-01, 4.3974420e-01, 4.8722979e-01, 5.8590513e-01,
+                4.5599860e-01, 3.1293562e-01, 1.9222400e-01, 9.7611703e-02,
+                3.3851895e-02, 2.4188591e-02, 1.9636340e-02, 1.3214327e-02,
+                7.2453087e-03, 4.3027173e-03, 3.7014042e-03, 2.9414182e-03,
+                2.7806845e-03, 2.8306348e-03, 2.9053832e-03, 2.9848625e-03,
+                3.0779461e-03, 3.1315640e-03, 3.2939769e-03, 3.3872949e-03,
+                3.7360021e-03])
+SST77, MSL77 = 28.20263671875, 1014.9654541015625
+
+P95 = np.array([1000.0, 975.0, 950.0, 925.0, 900.0, 875.0, 850.0, 825.0, 800.0,
+                775.0, 750.0, 700.0, 650.0, 600.0, 550.0, 500.0, 450.0, 400.0,
+                350.0, 300.0, 250.0, 225.0, 200.0, 175.0, 150.0, 125.0, 100.0,
+                70.0, 50.0, 30.0, 20.0, 10.0, 7.0, 5.0, 3.0, 2.0, 1.0])
+TC95 = np.array([20.266, 18.11, 15.978, 13.94, 12.576, 11.213, 9.669, 8.084,
+                 6.516, 4.993, 4.374, 4.93, 1.852, -2.047, -6.917, -13.015,
+                 -19.698, -26.881, -34.178, -42.822, -51.573, -54.68, -57.207,
+                 -58.943, -60.751, -63.255, -63.277, -61.5, -58.453, -53.327,
+                 -50.156, -46.602, -42.283, -35.389, -22.349, -14.944, -11.886])
+R95 = np.array([9.65, 9.517, 9.459, 9.205, 7.784, 6.561, 5.952, 5.517, 4.976,
+                4.332, 3.102, 1.207, 1.521, 1.263, 0.857, 0.556, 0.33, 0.211,
+                0.104, 0.073, 0.037, 0.027, 0.018, 0.01, 0.005, 0.004, 0.003,
+                0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.004, 0.004])
+SST95, MSL95 = 25.267, 1012.9
+
+save = {}
+
+
+def prep(SSTC, MSL, TCp, Rp):
+    SSTK = utilities.T_Ctok(SSTC)
+    T = utilities.T_Ctok(TCp)
+    R = Rp * 0.001
+    ES0 = utilities.es_cc(SSTC)
+    return SSTK, MSL, T, R, ES0
+
+
+# ---------- case A: PR77 ----------
+print("case A: PR77 ...")
+S77 = prep(SST77, MSL77, TC77, R77)
+pms = np.arange(946.0, 958.0, 0.002)
+sc = g_scan(pms, S77[0], S77[1], P77, S77[2], S77[3], S77[4])
+save["a_pm"] = pms
+save["a_g"] = sc[:, 0]
+save["a_capem"] = sc[:, 1]
+save["a_capems"] = sc[:, 2]
+st, xs = iterate_raw(S77[0], S77[1], P77, S77[2], S77[3], S77[4], nmax=60)
+save["a_iter_status"] = st
+save["a_iter"] = xs
+# buoyancy + W at the two cycle phases
+for tag, pm in [("lo", 950.6533454438), ("hi", 951.2790079839)]:
+    PP = min(pm, 1000.0)
+    RP = EPS * S77[3][0] * MSL77 / (PP * (EPS + S77[3][0]) - S77[3][0] * MSL77)
+    pts, Ws, b, Pt = W_curve(S77[2][0], RP, PP, S77[2], S77[3], P77)
+    save[f"a_W_{tag}_p"] = pts
+    save[f"a_W_{tag}_W"] = Ws
+    save[f"a_b_{tag}"] = b
+save["a_blevels"] = P77[: int(np.count_nonzero(P77 > 50))]
+
+# ---------- case B: 1995 problematic ----------
+print("case B: 1995 ...")
+S95 = prep(SST95, MSL95, TC95, R95)
+pms = np.arange(955.0, 968.0, 0.002)
+sc = g_scan(pms, S95[0], S95[1], P95, S95[2], S95[3], S95[4])
+save["b_pm"] = pms
+save["b_g"] = sc[:, 0]
+save["b_capem"] = sc[:, 1]
+st, xs = iterate_raw(S95[0], S95[1], P95, S95[2], S95[3], S95[4], nmax=60)
+save["b_iter_status"] = st
+save["b_iter"] = xs
+
+# ---------- sample columns ----------
+print("loading sample ...")
+ds = xr.open_dataset(os.path.join(SCRATCH, "sample_data.nc"))
+p_s = ds["p"].values.astype(float)
+sst_s, msl_s = ds["sst"].values, ds["msl"].values
+t_s, r_s = ds["t"].values, ds["r"].values
+a_ref = np.load(os.path.join(SCRATCH, "branch_pi.npz"))
+both1 = a_ref["IFL"] == 1
+
+
+def col(m, j, i):
+    T = utilities.T_Ctok(t_s[m, :, j, i].astype(float))
+    R = r_s[m, :, j, i].astype(float) * 0.001
+    R[np.isnan(R)] = 0.0
+    SSTK = utilities.T_Ctok(sst_s[m, j, i])
+    ES0 = utilities.es_cc(sst_s[m, j, i])
+    return SSTK, float(msl_s[m, j, i]), T, R, ES0
+
+
+# ---------- case C: healthy column (strong PI, mid-domain) ----------
+# pick the column with max VMAX (deep tropics, textbook profile)
+mm, jj, ii = np.unravel_index(np.nanargmax(np.where(both1, a_ref["VMAX"], -1)),
+                              a_ref["VMAX"].shape)
+print(f"case C: healthy col (m={mm},j={jj},i={ii}) VMAX={a_ref['VMAX'][mm,jj,ii]:.1f}")
+SC = col(mm, jj, ii)
+pms = np.arange(850.0, 1005.0, 0.02)
+sc = g_scan(pms, SC[0], SC[1], p_s, SC[2], SC[3], SC[4])
+save["c_pm"] = pms
+save["c_g"] = sc[:, 0]
+st, xs = iterate_raw(SC[0], SC[1], p_s, SC[2], SC[3], SC[4], nmax=60)
+save["c_iter"] = xs
+save["c_idx"] = np.array([mm, jj, ii])
+# W curve + buoyancy for the saturated parcel at converged PM (textbook single hump)
+pmc = xs[-1]
+PP = min(pmc, 1000.0)
+pts, Ws, b, Pt = W_curve(SC[0], utilities.rv(SC[4], PP), PP, SC[2], SC[3], p_s)
+save["c_W_p"] = pts
+save["c_W_W"] = Ws
+save["c_b"] = b
+save["c_blevels"] = p_s[: int(np.count_nonzero(p_s > 50))]
+
+# ---------- survey: jump statistics over converging columns ----------
+print("survey ...")
+idxs = np.argwhere(both1)
+rng = np.random.default_rng(42)
+sel = idxs[rng.choice(len(idxs), size=min(1200, len(idxs)), replace=False)]
+pmgrid = np.arange(900.0, 1005.0, 0.1)
+maxjump = np.full(len(sel), np.nan)
+njumps = np.zeros(len(sel))
+fp_in_jump = np.zeros(len(sel), bool)
+dist_fp_jump = np.full(len(sel), np.nan)
+SMOOTH = 0.06  # smooth |dg| bound for 0.1 step given slope<=0.5
+for k, (m, j, i) in enumerate(sel):
+    S = col(m, j, i)
+    gs = np.empty(len(pmgrid))
+    for q, pm in enumerate(pmgrid):
+        gs[q] = g_eval(pm, S[0], S[1], p_s, S[2], S[3], S[4])[0]
+    dg = np.diff(gs)
+    jmask = np.abs(dg) > SMOOTH
+    njumps[k] = int(jmask.sum())
+    if jmask.any():
+        maxjump[k] = np.max(np.abs(dg[jmask]))
+    # diagonal crossing
+    F = gs - pmgrid
+    scg = np.where(np.diff(np.sign(F)) != 0)[0]
+    if len(scg):
+        q = scg[0]
+        fp_in_jump[k] = jmask[q]
+        if jmask.any():
+            jpos = pmgrid[:-1][jmask] + 0.05
+            dist_fp_jump[k] = np.min(np.abs(jpos - pmgrid[q]))
+save["s_maxjump"] = maxjump
+save["s_njumps"] = njumps
+save["s_fp_in_jump"] = fp_in_jump
+save["s_dist"] = dist_fp_jump
+save["s_n"] = len(sel)
+print(f"  columns with >=1 jump: {(njumps>0).sum()}/{len(sel)}, "
+      f"fp inside jump: {fp_in_jump.sum()}")
+
+# ---------- case D: SST sweep on PR77 column ----------
+print("SST sweep ...")
+dss = np.arange(-3.0, 3.0001, 0.01)
+sweep_status = np.zeros(len(dss))
+sweep_pm = np.full(len(dss), np.nan)
+for k, d in enumerate(dss):
+    S = prep(SST77 + d, MSL77, TC77, R77)
+    st, xs = iterate_raw(S[0], S[1], P77, S[2], S[3], S[4], nmax=200)
+    sweep_status[k] = st
+    if st == 1:
+        sweep_pm[k] = xs[-1]
+save["d_dss"] = dss
+save["d_status"] = sweep_status
+save["d_pm"] = sweep_pm
+print(f"  failures: {(sweep_status==2).sum()}/{len(dss)}")
+
+np.savez(os.path.join(SCRATCH, "docdata.npz"), **save)
+print("saved docdata.npz")
diff --git a/discontinuity_analysis/scripts/dissect.py b/discontinuity_analysis/scripts/dissect.py
new file mode 100644
index 0000000..1526945
--- /dev/null
+++ b/discontinuity_analysis/scripts/dissect.py
@@ -0,0 +1,132 @@
+"""Classify branch-vs-fix differences: evaluate BOTH cape definitions at
+IDENTICAL inputs (same PM), removing the iteration-path effect.
+
+For each changed column, at the branch's own converged PM:
+ - compare cape_pcmin vs cape_maxw for the three calls (A, M, S)
+ - classify: identical (path noise) / surface-sliver / dominated-hump
+"""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.join(SCRATCH, "branch", "src"))  # pcmin cape via cape_instr
+
+import numpy as np
+import xarray as xr
+from tcpyPI import constants, utilities
+
+sys.path.insert(0, SCRATCH)
+from gmap import cape_instr  # pcmin-convention instrumented cape
+
+
+def cape_maxw_generic(TP, RP, PP, T, R, P, ascent_flag=0, ptop=50):
+    base = cape_instr(TP, RP, PP, T, R, P, ascent_flag, ptop)
+    if base["INB"] == -1:
+        return base
+    N = int(np.count_nonzero(P > ptop))
+    Pt, Tt = P[:N], T[:N]
+    TVRDIF = base["TVRDIF"]
+    RD = constants.RD
+    W = RD * (PP - Pt[0]) / (PP + Pt[0]) * TVRDIF[0]
+    BESTW, BESTP, BESTT = 0.0, 0.0, Tt[0]
+    if W > BESTW:
+        BESTW, BESTP, BESTT = W, Pt[0], Tt[0]
+    for j in range(1, N):
+        b0, b1 = TVRDIF[j - 1], TVRDIF[j]
+        if b0 * b1 < 0.0:
+            Pc = (Pt[j] * b0 - Pt[j - 1] * b1) / (b0 - b1)
+            Wc = W + RD * b0 * (Pt[j - 1] - Pc) / (Pt[j - 1] + Pc)
+            if Wc > BESTW:
+                Tc = (Tt[j - 1] * (Pc - Pt[j]) + Tt[j] * (Pt[j - 1] - Pc)) / (Pt[j - 1] - Pt[j])
+                BESTW, BESTP, BESTT = Wc, Pc, Tc
+        W += RD * (b1 + b0) * (Pt[j - 1] - Pt[j]) / (Pt[j] + Pt[j - 1])
+        if W > BESTW:
+            BESTW, BESTP, BESTT = W, Pt[j], Tt[j]
+    if BESTW <= 0.0:
+        return dict(CAPED=0.0, TOB=Tt[0], LNB=0.0, IFLAG=1)
+    return dict(CAPED=BESTW, TOB=BESTT, LNB=BESTP, IFLAG=1)
+
+
+ds = xr.open_dataset(os.path.join(SCRATCH, "sample_data.nc"))
+p = ds["p"].values.astype(float)
+sst, msl = ds["sst"].values, ds["msl"].values
+t4, r4 = ds["t"].values, ds["r"].values
+
+a = np.load(os.path.join(SCRATCH, "branch_pi.npz"))
+b = np.load(os.path.join(SCRATCH, "branch_fix_pi.npz"))
+both1 = (a["IFL"] == 1) & (b["IFL"] == 1)
+d_vmax = np.abs(a["VMAX"] - b["VMAX"])
+
+chg = np.argwhere(both1 & (d_vmax > 0.1))
+print(f"dissecting {len(chg)} columns with |dVMAX| > 0.1 m/s")
+
+cls = {"A": 0, "M": 0, "S": 0, "surface_sliver": 0, "none@fixedPM": 0}
+details = []
+for idx in chg:
+    m, j, i = idx
+    T = utilities.T_Ctok(t4[m, :, j, i].astype(float))
+    R = r4[m, :, j, i].astype(float) * 0.001
+    R[np.isnan(R)] = 0.0
+    SSTK = utilities.T_Ctok(sst[m, j, i])
+    MSL = msl[m, j, i]
+    ES0 = utilities.es_cc(sst[m, j, i])
+    # reconstruct branch converged PM from its PMIN? PM isn't stored; instead
+    # evaluate both definitions at the same probe PM = branch PMIN's implied
+    # iterate is unknown, so use PM=970 start AND branch PMIN as two probes.
+    diffs = set()
+    sliver = False
+    for PM in (970.0, max(a["PMIN"][m, j, i], 401.0)):
+        PP = min(PM, 1000.0)
+        RP = constants.EPS * R[0] * MSL / (PP * (constants.EPS + R[0]) - R[0] * MSL)
+        for tag, (TPq, RPq, PPq) in {
+            "A": (T[0], R[0], p[0]),
+            "M": (T[0], RP, PP),
+            "S": (SSTK, utilities.rv(ES0, PP), PP),
+        }.items():
+            o1 = cape_instr(TPq, RPq, PPq, T, R, p, 0, 50)
+            o2 = cape_maxw_generic(TPq, RPq, PPq, T, R, p, 0, 50)
+            if abs(o1["CAPED"] - o2["CAPED"]) > 1e-9 or abs((o1["LNB"] or 0) - o2["LNB"]) > 0.5:
+                diffs.add(tag)
+                if o1["INB"] == 0 and o2["CAPED"] > 0:
+                    sliver = True
+    if not diffs:
+        cls["none@fixedPM"] += 1
+    elif sliver:
+        cls["surface_sliver"] += 1
+    else:
+        for tag in diffs:
+            cls[tag] += 1
+    details.append((tuple(idx), sorted(diffs), a["VMAX"][m, j, i], b["VMAX"][m, j, i],
+                    a["TO"][m, j, i], b["TO"][m, j, i]))
+
+print("classification counts:", cls)
+print()
+for d in details[:25]:
+    print(f"{d[0]} calls-differing={d[1]} VMAX {d[2]:.3f}->{d[3]:.3f} TO {d[4]:.2f}->{d[5]:.2f}")
+
+# Also: how many of ALL both1 columns have identical cape at fixed PM=970 for all 3 calls?
+print()
+print("=== sampling 300 unchanged columns to confirm exact equality at fixed PM ===")
+ok_cols = np.argwhere(both1 & (d_vmax <= 1e-8))
+rng = np.random.default_rng(0)
+sel = ok_cols[rng.choice(len(ok_cols), size=300, replace=False)]
+neq = 0
+for idx in sel:
+    m, j, i = idx
+    T = utilities.T_Ctok(t4[m, :, j, i].astype(float))
+    R = r4[m, :, j, i].astype(float) * 0.001
+    R[np.isnan(R)] = 0.0
+    SSTK = utilities.T_Ctok(sst[m, j, i])
+    MSL = msl[m, j, i]
+    ES0 = utilities.es_cc(sst[m, j, i])
+    PM = 970.0
+    PP = min(PM, 1000.0)
+    RP = constants.EPS * R[0] * MSL / (PP * (constants.EPS + R[0]) - R[0] * MSL)
+    for TPq, RPq, PPq in [(T[0], R[0], p[0]), (T[0], RP, PP), (SSTK, utilities.rv(ES0, PP), PP)]:
+        o1 = cape_instr(TPq, RPq, PPq, T, R, p, 0, 50)
+        o2 = cape_maxw_generic(TPq, RPq, PPq, T, R, p, 0, 50)
+        if o1["CAPED"] != o2["CAPED"]:
+            neq += 1
+            break
+print(f"columns (of 300 sampled unchanged) with ANY cape difference at PM=970: {neq}")
diff --git a/discontinuity_analysis/scripts/final_check.py b/discontinuity_analysis/scripts/final_check.py
new file mode 100644
index 0000000..ee29610
--- /dev/null
+++ b/discontinuity_analysis/scripts/final_check.py
@@ -0,0 +1,75 @@
+"""Final confirmation: PR#77 + 1995 problematic/control profiles through the
+patched numba kernel (max-W CAPE, rescue deleted)."""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+tree = sys.argv[1]  # 'branch' or 'branch_fix'
+sys.path.insert(0, os.path.join(SCRATCH, tree, "src"))
+
+import numpy as np
+import tcpyPI
+
+OPT = dict(CKCD=0.9, ascent_flag=0, diss_flag=1, ptop=50, miss_handle=1)
+
+# --- PR #77 profile ---
+P77 = np.array([1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650,
+                600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125,
+                100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1], dtype=float)
+TC77 = np.array([25.260956, 23.078949, 20.881561, 18.652649, 16.615143, 14.63147,
+                 12.794586, 11.79306, 11.01236, 10.847565, 10.354492, 7.807007,
+                 5.473297, 2.5278625, -1.532135, -7.3142395, -13.635345,
+                 -20.613602, -28.588928, -37.270096, -45.40825, -50.047455,
+                 -54.84575, -59.173737, -62.878662, -65.78009, -69.30669,
+                 -65.512024, -60.76361, -55.100464, -51.38333, -43.95462,
+                 -39.709717, -33.193268, -21.576843, -15.402496, -13.186676])
+R77 = np.array([1.0783693e01, 1.0704287e01, 1.0680210e01, 1.0617845e01,
+                1.0320683e01, 9.8811483e00, 9.1884289e00, 7.1884680e00,
+                5.6963191e00, 3.5568204e00, 1.5912720e00, 1.0433695e00,
+                5.9723043e-01, 4.3974420e-01, 4.8722979e-01, 5.8590513e-01,
+                4.5599860e-01, 3.1293562e-01, 1.9222400e-01, 9.7611703e-02,
+                3.3851895e-02, 2.4188591e-02, 1.9636340e-02, 1.3214327e-02,
+                7.2453087e-03, 4.3027173e-03, 3.7014042e-03, 2.9414182e-03,
+                2.7806845e-03, 2.8306348e-03, 2.9053832e-03, 2.9848625e-03,
+                3.0779461e-03, 3.1315640e-03, 3.2939769e-03, 3.3872949e-03,
+                3.7360021e-03])
+out = tcpyPI.pi(28.20263671875, 1014.9654541015625, P77, TC77, R77, **OPT)
+print(f"[{tree}] PR#77 profile:  VMAX={out[0]:.6f} PMIN={out[1]:.6f} IFL={out[2]} TO={out[3]:.4f} OTL={out[4]:.4f}")
+
+# --- 1995 problematic + control (from branch test) ---
+PROB = np.array([
+    [1000.0, 20.266, 9.65], [975.0, 18.11, 9.517], [950.0, 15.978, 9.459],
+    [925.0, 13.94, 9.205], [900.0, 12.576, 7.784], [875.0, 11.213, 6.561],
+    [850.0, 9.669, 5.952], [825.0, 8.084, 5.517], [800.0, 6.516, 4.976],
+    [775.0, 4.993, 4.332], [750.0, 4.374, 3.102], [700.0, 4.93, 1.207],
+    [650.0, 1.852, 1.521], [600.0, -2.047, 1.263], [550.0, -6.917, 0.857],
+    [500.0, -13.015, 0.556], [450.0, -19.698, 0.33], [400.0, -26.881, 0.211],
+    [350.0, -34.178, 0.104], [300.0, -42.822, 0.073], [250.0, -51.573, 0.037],
+    [225.0, -54.68, 0.027], [200.0, -57.207, 0.018], [175.0, -58.943, 0.01],
+    [150.0, -60.751, 0.005], [125.0, -63.255, 0.004], [100.0, -63.277, 0.003],
+    [70.0, -61.5, 0.003], [50.0, -58.453, 0.003], [30.0, -53.327, 0.003],
+    [20.0, -50.156, 0.003], [10.0, -46.602, 0.003], [7.0, -42.283, 0.003],
+    [5.0, -35.389, 0.003], [3.0, -22.349, 0.003], [2.0, -14.944, 0.004],
+    [1.0, -11.886, 0.004]])
+out = tcpyPI.pi(25.267, 1012.9, PROB[:, 0], PROB[:, 1], PROB[:, 2], **OPT)
+print(f"[{tree}] 1995 problematic: VMAX={out[0]:.6f} PMIN={out[1]:.6f} IFL={out[2]} TO={out[3]:.4f} OTL={out[4]:.4f}")
+
+CTRL = np.array([
+    [1000.0, 21.326, 11.215], [975.0, 19.192, 11.011], [950.0, 17.315, 10.556],
+    [925.0, 15.676, 9.931], [900.0, 14.286, 8.639], [875.0, 12.957, 7.58],
+    [850.0, 11.495, 6.97], [825.0, 10.019, 6.439], [800.0, 8.575, 5.734],
+    [775.0, 6.993, 5.198], [750.0, 5.373, 4.752], [700.0, 3.13, 3.122],
+    [650.0, 0.534, 2.562], [600.0, -2.573, 1.442], [550.0, -6.616, 0.865],
+    [500.0, -12.233, 0.61], [450.0, -18.954, 0.42], [400.0, -26.098, 0.292],
+    [350.0, -33.337, 0.184], [300.0, -40.639, 0.112], [250.0, -49.373, 0.057],
+    [225.0, -52.236, 0.03], [200.0, -55.706, 0.017], [175.0, -57.323, 0.01],
+    [150.0, -60.084, 0.005], [125.0, -65.224, 0.004], [100.0, -65.651, 0.003],
+    [70.0, -63.457, 0.003], [50.0, -59.522, 0.003], [30.0, -53.504, 0.003],
+    [20.0, -49.926, 0.003], [10.0, -45.925, 0.003], [7.0, -41.476, 0.003],
+    [5.0, -34.68, 0.003], [3.0, -23.495, 0.003], [2.0, -16.575, 0.004],
+    [1.0, -12.898, 0.004]])
+out = tcpyPI.pi(25.292, 1011.4, CTRL[:, 0], CTRL[:, 1], CTRL[:, 2], **OPT)
+ref = (68.98260512391965, 943.0252390748838, 214.15121722308447, 159.82630589536842)
+print(f"[{tree}] 1995 control:     VMAX={out[0]:.6f} PMIN={out[1]:.6f} IFL={out[2]} TO={out[3]:.4f} OTL={out[4]:.4f}")
+print(f"        control vs pinned: dVMAX={out[0]-ref[0]:+.3e} dPMIN={out[1]-ref[1]:+.3e} dTO={out[3]-ref[2]:+.3e} dOTL={out[4]-ref[3]:+.3e}")
diff --git a/discontinuity_analysis/scripts/fixdata.py b/discontinuity_analysis/scripts/fixdata.py
new file mode 100644
index 0000000..d4819a7
--- /dev/null
+++ b/discontinuity_analysis/scripts/fixdata.py
@@ -0,0 +1,33 @@
+"""g_fix(PM) scan for PR77 (max-W CAPE) + iterates; append-style save."""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, SCRATCH)
+
+import numpy as np
+from maxw import gmap_fix
+
+pms = np.arange(946.0, 958.0, 0.01)
+gs = np.array([gmap_fix(x)["PNEW"] for x in pms])
+
+# iterates
+PM, PMOLD, PNEW, NP = 970.0, 970.0, 0.0, 0
+xs = []
+while abs(PNEW - PMOLD) > 0.5 and NP <= 60:
+    PNEW = gmap_fix(PM)["PNEW"]
+    xs.append(PM)
+    PMOLD, PM = PM, PNEW
+    NP += 1
+xs.append(PM)
+
+np.savez(os.path.join(SCRATCH, "fixdata.npz"), pm=pms, g=gs, iters=np.array(xs))
+print(f"g_fix scan saved; converged in {NP} iterations to {PM:.6f}")
+
+# sanity: humpdata cape_pc continuity check
+h = np.load(os.path.join(SCRATCH, "humpdata.npz"))
+dpc = np.max(np.abs(np.diff(h["cape_pc"])))
+dmw = np.max(np.abs(np.diff(h["cape_mw"])))
+dlnb = np.nanmax(np.abs(np.diff(h["lnb_mw"])))
+print(f"hump sweep: max step cape_pcmin={dpc:.3f} cape_maxW={dmw:.3f} lnb_maxW={dlnb:.1f}")
diff --git a/discontinuity_analysis/scripts/gmap.py b/discontinuity_analysis/scripts/gmap.py
new file mode 100644
index 0000000..9913684
--- /dev/null
+++ b/discontinuity_analysis/scripts/gmap.py
@@ -0,0 +1,203 @@
+"""Instrumented reconstruction of the pcmin fixed-point map g(PM) for PR #77.
+
+Replicates _pi_numba's iteration body (branch upstream/modernize_eg) exactly,
+but exposes every internal quantity so discontinuities can be located and
+classified. Pure-python (numba disabled) for introspection.
+"""
+
+import os
+import sys
+
+os.environ["TCPYPI_DISABLE_NUMBA"] = "1"
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.join(SCRATCH, "branch", "src"))
+
+import numpy as np
+from tcpyPI import constants, utilities
+from tcpyPI.pi import solve_temperature_from_entropy
+
+# ---- PR #77 profile ----
+SSTC = 28.20263671875
+MSL = 1014.9654541015625
+P_FULL = np.array(
+    [1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750]
+    + [700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225]
+    + [200, 175, 150, 125, 100, 70, 50, 30, 20, 10, 7]
+    + [5, 3, 2, 1],
+    dtype=float,
+)
+TC_FULL = np.array(
+    [25.260956, 23.078949, 20.881561, 18.652649, 16.615143]
+    + [14.63147, 12.794586, 11.79306, 11.01236, 10.847565]
+    + [10.354492, 7.807007, 5.473297, 2.5278625, -1.532135]
+    + [-7.3142395, -13.635345, -20.613602, -28.588928, -37.270096]
+    + [-45.40825, -50.047455, -54.84575, -59.173737, -62.878662]
+    + [-65.78009, -69.30669, -65.512024, -60.76361, -55.100464]
+    + [-51.38333, -43.95462, -39.709717, -33.193268, -21.576843]
+    + [-15.402496, -13.186676]
+)
+R_FULL = np.array(
+    [1.0783693e01, 1.0704287e01, 1.0680210e01, 1.0617845e01]
+    + [1.0320683e01, 9.8811483e00, 9.1884289e00, 7.1884680e00]
+    + [5.6963191e00, 3.5568204e00, 1.5912720e00, 1.0433695e00]
+    + [5.9723043e-01, 4.3974420e-01, 4.8722979e-01, 5.8590513e-01]
+    + [4.5599860e-01, 3.1293562e-01, 1.9222400e-01, 9.7611703e-02]
+    + [3.3851895e-02, 2.4188591e-02, 1.9636340e-02, 1.3214327e-02]
+    + [7.2453087e-03, 4.3027173e-03, 3.7014042e-03, 2.9414182e-03]
+    + [2.7806845e-03, 2.8306348e-03, 2.9053832e-03, 2.9848625e-03]
+    + [3.0779461e-03, 3.1315640e-03, 3.2939769e-03, 3.3872949e-03]
+    + [3.7360021e-03]
+)
+
+CKCD = 0.9
+ASCENT_FLAG = 0
+PTOP = 50
+
+SSTK = utilities.T_Ctok(SSTC)
+T = utilities.T_Ctok(TC_FULL)
+R = R_FULL * 0.001
+ES0 = utilities.es_cc(SSTC)
+NK = 0
+
+
+def cape_instr(TP, RP, PP, T, R, P, ascent_flag=0, ptop=50):
+    """Copy of branch cape() with full instrumentation. Returns dict."""
+    N = int(np.count_nonzero(P > ptop))
+    P = P[:N]
+    T = T[:N]
+    R = R[:N]
+    nlvl = len(P)
+    TVRDIF = np.zeros((nlvl,))
+
+    if (RP < 1e-6) or (TP < 200):
+        return dict(CAPED=0.0, TOB=np.nan, LNB=np.nan, IFLAG=0, INB=-1,
+                    PLCL=np.nan, TVRDIF=TVRDIF, PA=0.0, NA=0.0, PAT=0.0)
+
+    TPC = utilities.T_ktoC(TP)
+    ESP = utilities.es_cc(TPC)
+    EVP = utilities.ev(RP, PP)
+    RH = min(EVP / ESP, 1.0)
+    S = utilities.entropy_S(TP, RP, PP)
+    PLCL = utilities.e_pLCL(TP, RH, PP)
+
+    CAPED = 0
+    TOB = T[0]
+    IFLAG = 1
+    jmin = int(1e6)
+    ncsteps = []
+
+    for j in range(nlvl):
+        jmin = int(min(jmin, j))
+        if P[j] >= PLCL:
+            TG = TP * (P[j] / PP) ** (constants.RD / constants.CPD)
+            RG = RP
+            TLVR = utilities.Trho(TG, RG, RG)
+            TVENV = utilities.Trho(T[j], R[j], R[j])
+            TVRDIF[j] = TLVR - TVENV
+        else:
+            TG, RG, IFLAG_N = solve_temperature_from_entropy(
+                S=S, P=P[j], RP=RP, T_initial=T[j]
+            )
+            if IFLAG_N == 2:
+                return dict(CAPED=0.0, TOB=T[0], LNB=P[0], IFLAG=2, INB=-1,
+                            PLCL=PLCL, TVRDIF=TVRDIF, PA=0.0, NA=0.0, PAT=0.0)
+            RMEAN = ascent_flag * RG + (1 - ascent_flag) * RP
+            TLVR = utilities.Trho(TG, RMEAN, RG)
+            TENV = utilities.Trho(T[j], R[j], R[j])
+            TVRDIF[j] = TLVR - TENV
+
+    NA = 0.0
+    PA = 0.0
+    INB = 0
+    for j in range(nlvl - 1, jmin, -1):
+        if TVRDIF[j] > 0:
+            INB = max(INB, j)
+
+    if INB == 0:
+        return dict(CAPED=0.0, TOB=T[0], LNB=0.0, IFLAG=IFLAG, INB=0,
+                    PLCL=PLCL, TVRDIF=TVRDIF, PA=0.0, NA=0.0, PAT=0.0)
+
+    for j in range(jmin + 1, INB + 1, 1):
+        PFAC = (
+            constants.RD * (TVRDIF[j] + TVRDIF[j - 1]) * (P[j - 1] - P[j])
+            / (P[j] + P[j - 1])
+        )
+        PA = PA + max(PFAC, 0.0)
+        NA = NA - min(PFAC, 0.0)
+
+    PMA = PP + P[jmin]
+    PFAC = constants.RD * (PP - P[jmin]) / PMA
+    PA = PA + PFAC * max(TVRDIF[jmin], 0.0)
+    NA = NA - PFAC * min(TVRDIF[jmin], 0.0)
+
+    PAT = 0.0
+    TOB = T[INB]
+    LNB = P[INB]
+    if INB < nlvl - 1:
+        PINB = (P[INB + 1] * TVRDIF[INB] - P[INB] * TVRDIF[INB + 1]) / (
+            TVRDIF[INB] - TVRDIF[INB + 1]
+        )
+        LNB = PINB
+        PAT = constants.RD * TVRDIF[INB] * (P[INB] - PINB) / (P[INB] + PINB)
+        TOB = (T[INB] * (PINB - P[INB + 1]) + T[INB + 1] * (P[INB] - PINB)) / (
+            P[INB] - P[INB + 1]
+        )
+
+    CAPED = max(PA + PAT - NA, 0.0)
+    return dict(CAPED=CAPED, TOB=TOB, LNB=LNB, IFLAG=1, INB=INB,
+                PLCL=PLCL, TVRDIF=TVRDIF, PA=PA, NA=NA, PAT=PAT)
+
+
+# Environmental CAPE (PM-independent)
+_env = cape_instr(T[NK], R[NK], P_FULL[NK], T, R, P_FULL, ASCENT_FLAG, PTOP)
+CAPEA = _env["CAPED"]
+
+
+def gmap(PM, diss_flag=1):
+    """One iteration body: PM -> PNEW, with internals."""
+    PP = min(PM, 1000.0)
+    RP_M = constants.EPS * R[NK] * MSL / (
+        PP * (constants.EPS + R[NK]) - R[NK] * MSL
+    )
+    m = cape_instr(T[NK], RP_M, PP, T, R, P_FULL, ASCENT_FLAG, PTOP)
+
+    RP_S = utilities.rv(ES0, PP)
+    s = cape_instr(SSTK, RP_S, PP, T, R, P_FULL, ASCENT_FLAG, PTOP)
+
+    TO = s["TOB"]
+    RAT = SSTK / TO if diss_flag else 1.0
+
+    TV0 = utilities.Trho(T[NK], R[NK], R[NK])
+    TVSST = utilities.Trho(SSTK, RP_S, RP_S)
+    TVAV = 0.5 * (TV0 + TVSST)
+    CAT_raw = (m["CAPED"] - CAPEA) + 0.5 * CKCD * RAT * (s["CAPED"] - m["CAPED"])
+    CAT = max(CAT_raw, 0.0)
+    PNEW = MSL * np.exp(-CAT / (constants.RD * TVAV))
+    return dict(PNEW=PNEW, CAPEM=m["CAPED"], CAPEMS=s["CAPED"], TO=TO,
+                RAT=RAT, CAT_raw=CAT_raw, CAT=CAT,
+                INB_M=m["INB"], INB_S=s["INB"],
+                PLCL_M=m["PLCL"], PLCL_S=s["PLCL"],
+                LNB_M=m["LNB"], LNB_S=s["LNB"],
+                PA_M=m["PA"], NA_M=m["NA"], PAT_M=m["PAT"],
+                PA_S=s["PA"], NA_S=s["NA"], PAT_S=s["PAT"],
+                TVRDIF_M=m["TVRDIF"], TVRDIF_S=s["TVRDIF"])
+
+
+if __name__ == "__main__":
+    # 1) Reproduce the raw (no-rescue) iteration to find the 2-cycle
+    print("=== Raw fixed-point iteration (no rescue), x0=970 ===")
+    PM, PMOLD, PNEW, NP = 970.0, 970.0, 0.0, 0
+    hist = []
+    while abs(PNEW - PMOLD) > 0.5 and NP <= 30:
+        out = gmap(PM)
+        PNEW = out["PNEW"]
+        hist.append((NP, PM, PNEW, out["CAPEM"], out["CAPEMS"], out["TO"],
+                     out["INB_M"], out["INB_S"], out["CAT_raw"]))
+        PMOLD, PM = PM, PNEW
+        NP += 1
+    for h in hist[:30]:
+        print(
+            f"n={h[0]:3d} PM={h[1]:.10f} -> PNEW={h[2]:.10f} "
+            f"CAPEM={h[3]:9.4f} CAPEMS={h[4]:9.4f} TO={h[5]:8.3f} "
+            f"INB_M={h[6]:2d} INB_S={h[7]:2d} CAT={h[8]:8.4f}"
+        )
diff --git a/discontinuity_analysis/scripts/humphunt.py b/discontinuity_analysis/scripts/humphunt.py
new file mode 100644
index 0000000..5dd7ab5
--- /dev/null
+++ b/discontinuity_analysis/scripts/humphunt.py
@@ -0,0 +1,123 @@
+"""Find a column with two comparable local maxima of W(p_t) (saturated parcel),
+and produce data showing: CAPE_maxW continuous vs SST while argmax (TO) jumps.
+"""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.join(SCRATCH, "branch", "src"))
+sys.path.insert(0, SCRATCH)
+
+import numpy as np
+import xarray as xr
+from tcpyPI import constants, utilities
+
+from gmap import cape_instr
+
+RD = constants.RD
+
+
+def W_all(TP, RP, PP, T, R, P, ptop=50):
+    base = cape_instr(TP, RP, PP, T, R, P, 0, ptop)
+    N = int(np.count_nonzero(P > ptop))
+    Pt, Tt = P[:N], T[:N]
+    b = base["TVRDIF"]
+    pts, Ws, Ts = [PP], [0.0], [TP]
+    W = RD * (PP - Pt[0]) / (PP + Pt[0]) * b[0]
+    pts.append(Pt[0]); Ws.append(W); Ts.append(Tt[0])
+    for j in range(1, N):
+        b0, b1 = b[j - 1], b[j]
+        if b0 * b1 < 0.0:
+            Pc = (Pt[j] * b0 - Pt[j - 1] * b1) / (b0 - b1)
+            Wc = W + RD * b0 * (Pt[j - 1] - Pc) / (Pt[j - 1] + Pc)
+            Tc = (Tt[j-1] * (Pc - Pt[j]) + Tt[j] * (Pt[j-1] - Pc)) / (Pt[j-1] - Pt[j])
+            pts.append(Pc); Ws.append(Wc); Ts.append(Tc)
+        W += RD * (b1 + b0) * (Pt[j - 1] - Pt[j]) / (Pt[j] + Pt[j - 1])
+        pts.append(Pt[j]); Ws.append(W); Ts.append(Tt[j])
+    return np.array(pts), np.array(Ws), np.array(Ts), b, Pt
+
+
+def local_maxima(pts, Ws):
+    out = []
+    for k in range(1, len(Ws) - 1):
+        if Ws[k] >= Ws[k - 1] and Ws[k] >= Ws[k + 1] and Ws[k] > 0:
+            out.append((pts[k], Ws[k]))
+    return out
+
+
+ds = xr.open_dataset(os.path.join(SCRATCH, "sample_data.nc"))
+p_s = ds["p"].values.astype(float)
+sst_s, msl_s = ds["sst"].values, ds["msl"].values
+t_s, r_s = ds["t"].values, ds["r"].values
+a_ref = np.load(os.path.join(SCRATCH, "branch_pi.npz"))
+both1 = a_ref["IFL"] == 1
+
+# hunt over ENVIRONMENTAL-parcel W curves (the M call is where #77-type marginal
+# humps live); require two maxima with ratio in [0.3, 1.0] and separation > 150 hPa
+idxs = np.argwhere(both1)
+rng = np.random.default_rng(1)
+sel = idxs[rng.choice(len(idxs), size=800, replace=False)]
+best = None
+for m, j, i in sel:
+    T = utilities.T_Ctok(t_s[m, :, j, i].astype(float))
+    R = r_s[m, :, j, i].astype(float) * 0.001
+    R[np.isnan(R)] = 0.0
+    MSL = float(msl_s[m, j, i])
+    PM = a_ref["PMIN"][m, j, i] + 20  # near converged radius pressure
+    PP = min(PM, 1000.0)
+    RP = constants.EPS * R[0] * MSL / (PP * (constants.EPS + R[0]) - R[0] * MSL)
+    if RP < 1e-6 or np.isnan(RP):
+        continue
+    pts, Ws, Ts, b, Pt = W_all(T[0], RP, PP, T, R, p_s)
+    lm = local_maxima(pts, Ws)
+    if len(lm) >= 2:
+        lm = sorted(lm, key=lambda x: -x[1])[:2]
+        ratio = lm[1][1] / lm[0][1]
+        sep = abs(lm[0][0] - lm[1][0])
+        if 0.25 < ratio and sep > 150:
+            score = ratio * sep
+            if best is None or score > best[0]:
+                best = (score, (m, j, i), ratio, sep)
+
+print("best two-hump column:", best)
+
+if best is not None:
+    m, j, i = best[1]
+    T = utilities.T_Ctok(t_s[m, :, j, i].astype(float))
+    R = r_s[m, :, j, i].astype(float) * 0.001
+    R[np.isnan(R)] = 0.0
+    MSL = float(msl_s[m, j, i])
+    SSTC0 = float(sst_s[m, j, i])
+    PM = a_ref["PMIN"][m, j, i] + 20
+    PP = min(PM, 1000.0)
+
+    # sweep a parameter (parcel moisture scaling) to force a dominance exchange:
+    # scale R[0] (boundary-layer moisture) by f in [0.9, 1.1]
+    fs = np.linspace(0.94, 1.06, 601)
+    cape_pc, cape_mw, lnb_mw, lnb_pc = [], [], [], []
+    for f in fs:
+        Rk = R.copy()
+        Rk[0] = R[0] * f
+        RP = constants.EPS * Rk[0] * MSL / (PP * (constants.EPS + Rk[0]) - Rk[0] * MSL)
+        pts, Ws, Ts, b, Pt = W_all(T[0], RP, PP, T, Rk, p_s)
+        kmax = int(np.argmax(Ws))
+        cape_mw.append(max(Ws[kmax], 0.0))
+        lnb_mw.append(pts[kmax] if Ws[kmax] > 0 else np.nan)
+        o = cape_instr(T[0], RP, PP, T, Rk, p_s, 0, 50)
+        cape_pc.append(o["CAPED"])
+        lnb_pc.append(o["LNB"])
+    np.savez(os.path.join(SCRATCH, "humpdata.npz"),
+             fs=fs, cape_pc=np.array(cape_pc), cape_mw=np.array(cape_mw),
+             lnb_mw=np.array(lnb_mw), lnb_pc=np.array(lnb_pc),
+             idx=np.array([m, j, i]), sst=SSTC0, pm=PM)
+    # W curves at three f values around the swap
+    swap_k = int(np.argmax(np.abs(np.diff(np.array(lnb_mw))))) if len(fs) > 1 else 300
+    for tag, f in [("lo", fs[max(swap_k - 60, 0)]), ("at", fs[swap_k]),
+                   ("hi", fs[min(swap_k + 60, len(fs) - 1)])]:
+        Rk = R.copy()
+        Rk[0] = R[0] * f
+        RP = constants.EPS * Rk[0] * MSL / (PP * (constants.EPS + Rk[0]) - Rk[0] * MSL)
+        pts, Ws, Ts, b, Pt = W_all(T[0], RP, PP, T, Rk, p_s)
+        np.savez(os.path.join(SCRATCH, f"humpW_{tag}.npz"), p=pts, W=Ws, f=f)
+    print(f"saved humpdata.npz; swap near f={fs[swap_k]:.4f}; col (m,j,i)={(m,j,i)} SST={SSTC0:.2f}")
diff --git a/discontinuity_analysis/scripts/maxw.py b/discontinuity_analysis/scripts/maxw.py
new file mode 100644
index 0000000..58b61b5
--- /dev/null
+++ b/discontinuity_analysis/scripts/maxw.py
@@ -0,0 +1,104 @@
+"""max-W CAPE: CAPE = max over outflow level p_t of the running signed
+buoyancy integral W(p_t), evaluated on the same piecewise-linear interpolant
+(trapezoid + interpolated-crossing partial areas) the current code uses.
+
+Candidates for the max: every grid node above jmin, and every interpolated
+sign-change crossing. Continuous in the profile/parcel parameters; equals the
+pcmin value exactly for single-crossing (typical) profiles.
+"""
+
+import numpy as np
+from tcpyPI import constants, utilities
+
+from gmap import (CAPEA as _unused, P_FULL, T, R, MSL, SSTK, ES0, NK, CKCD,
+                  ASCENT_FLAG, PTOP, cape_instr, gmap)
+
+
+def cape_maxw(TP, RP, PP, T, R, P, ascent_flag=0, ptop=50):
+    """cape() with the max-W integration; reuses cape_instr for TVRDIF."""
+    base = cape_instr(TP, RP, PP, T, R, P, ascent_flag, ptop)
+    if base["IFLAG"] != 1 and base["INB"] == -1:
+        return base  # improper parcel / non-convergence: unchanged
+    N = int(np.count_nonzero(P > ptop))
+    Pt, Tt = P[:N], T[:N]
+    TVRDIF = base["TVRDIF"]
+    nlvl = N
+    jmin = 0
+
+    RD = constants.RD
+    W = RD * (PP - Pt[jmin]) / (PP + Pt[jmin]) * TVRDIF[jmin]
+    BESTW, BESTP, BESTT = 0.0, 0.0, Tt[0]
+    if W > BESTW:
+        BESTW, BESTP, BESTT = W, Pt[jmin], Tt[jmin]
+    for j in range(jmin + 1, nlvl):
+        b0, b1 = TVRDIF[j - 1], TVRDIF[j]
+        if b0 * b1 < 0.0:
+            # interpolated crossing between j-1 and j (same formula as PINB/PAT)
+            Pc = (Pt[j] * b0 - Pt[j - 1] * b1) / (b0 - b1)
+            Wc = W + RD * b0 * (Pt[j - 1] - Pc) / (Pt[j - 1] + Pc)
+            if Wc > BESTW:
+                Tc = (Tt[j - 1] * (Pc - Pt[j]) + Tt[j] * (Pt[j - 1] - Pc)) / (
+                    Pt[j - 1] - Pt[j]
+                )
+                BESTW, BESTP, BESTT = Wc, Pc, Tc
+        W += RD * (b1 + b0) * (Pt[j - 1] - Pt[j]) / (Pt[j] + Pt[j - 1])
+        if W > BESTW:
+            BESTW, BESTP, BESTT = W, Pt[j], Tt[j]
+
+    if BESTW <= 0.0:
+        return dict(CAPED=0.0, TOB=Tt[0], LNB=0.0, IFLAG=1, INB=0,
+                    PLCL=base["PLCL"], TVRDIF=TVRDIF, PA=0.0, NA=0.0, PAT=0.0)
+    return dict(CAPED=BESTW, TOB=BESTT, LNB=BESTP, IFLAG=1, INB=-99,
+                PLCL=base["PLCL"], TVRDIF=TVRDIF, PA=np.nan, NA=np.nan,
+                PAT=np.nan)
+
+
+def gmap_fix(PM, diss_flag=1):
+    PP = min(PM, 1000.0)
+    RP_M = constants.EPS * R[NK] * MSL / (PP * (constants.EPS + R[NK]) - R[NK] * MSL)
+    m = cape_maxw(T[NK], RP_M, PP, T, R, P_FULL, ASCENT_FLAG, PTOP)
+    RP_S = utilities.rv(ES0, PP)
+    s = cape_maxw(SSTK, RP_S, PP, T, R, P_FULL, ASCENT_FLAG, PTOP)
+    env = cape_maxw(T[NK], R[NK], P_FULL[NK], T, R, P_FULL, ASCENT_FLAG, PTOP)
+    TO = s["TOB"]
+    RAT = SSTK / TO if diss_flag else 1.0
+    TV0 = utilities.Trho(T[NK], R[NK], R[NK])
+    TVSST = utilities.Trho(SSTK, RP_S, RP_S)
+    TVAV = 0.5 * (TV0 + TVSST)
+    CAT = max((m["CAPED"] - env["CAPED"]) + 0.5 * CKCD * RAT * (s["CAPED"] - m["CAPED"]), 0.0)
+    PNEW = MSL * np.exp(-CAT / (constants.RD * TVAV))
+    return dict(PNEW=PNEW, CAPEM=m["CAPED"], CAPEMS=s["CAPED"], TO=TO, RAT=RAT,
+                CAPEA=env["CAPED"], TVAV=TVAV)
+
+
+if __name__ == "__main__":
+    print("=== fixed-point iteration with max-W CAPE, x0=970 ===")
+    PM, PMOLD, PNEW, NP = 970.0, 970.0, 0.0, 0
+    while abs(PNEW - PMOLD) > 0.5 and NP <= 200:
+        o = gmap_fix(PM)
+        PNEW = o["PNEW"]
+        print(f"n={NP:3d} PM={PM:.10f} -> PNEW={PNEW:.10f}  CAPEM={o['CAPEM']:9.4f} TO={o['TO']:8.3f}")
+        PMOLD, PM = PM, PNEW
+        NP += 1
+    print(f"converged: NP={NP}, PM={PM:.8f}")
+
+    # final outputs (mirror _pi_numba post-loop)
+    o = gmap_fix(PM)
+    CATFAC = 0.5 * (1.0 + 1 / constants.b)
+    CAT = max((o["CAPEM"] - o["CAPEA"]) + CKCD * o["RAT"] * CATFAC * (o["CAPEMS"] - o["CAPEM"]), 0.0)
+    PMIN = MSL * np.exp(-CAT / (constants.RD * o["TVAV"]))
+    FAC = max(0.0, o["CAPEMS"] - o["CAPEM"])
+    VMAX = 0.8 * np.sqrt(CKCD * o["RAT"] * FAC)
+    print(f"VMAX={VMAX:.6f} m/s  PMIN={PMIN:.6f} hPa  TO={o['TO']:.4f} K")
+
+    # continuity check: rescan the old jump neighborhood
+    print()
+    print("=== g_fix(PM) over the old jump (950.7..950.95, step 2e-3) ===")
+    pms = np.arange(950.70, 950.95, 0.002)
+    gs = np.array([gmap_fix(x)["PNEW"] for x in pms])
+    dg = np.abs(np.diff(gs))
+    print(f"max |step-to-step change| = {dg.max():.3e} hPa (was 0.837 at the jump)")
+
+    # fixed-point location vs the rescue's midpoint answer
+    print()
+    print("comparison: branch rescue returns midpoint-ish ~950.966; PR cycle was (950.653, 951.279)")
diff --git a/discontinuity_analysis/scripts/plots.py b/discontinuity_analysis/scripts/plots.py
new file mode 100644
index 0000000..72d5e02
--- /dev/null
+++ b/discontinuity_analysis/scripts/plots.py
@@ -0,0 +1,412 @@
+"""All figures for the discontinuity analysis document."""
+
+import os
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+FIG = "/home/mares/repos/tcpyPI/discontinuity_analysis/figures"
+os.makedirs(FIG, exist_ok=True)
+
+# ---- palette (dataviz reference, light mode) ----
+BLUE = "#2a78d6"
+AQUA = "#1baf7a"
+YELLOW = "#eda100"
+GREEN = "#008300"
+VIOLET = "#4a3aa7"
+RED = "#e34948"
+INK = "#0b0b0b"
+INK2 = "#52514e"
+GRID = "#e5e4e0"
+
+plt.rcParams.update({
+    "font.size": 9,
+    "font.family": "serif",
+    "axes.edgecolor": INK2,
+    "axes.labelcolor": INK,
+    "axes.titlesize": 9.5,
+    "axes.titleweight": "bold",
+    "axes.spines.top": False,
+    "axes.spines.right": False,
+    "axes.grid": True,
+    "grid.color": GRID,
+    "grid.linewidth": 0.6,
+    "xtick.color": INK2,
+    "ytick.color": INK2,
+    "legend.frameon": False,
+    "legend.fontsize": 8,
+    "lines.linewidth": 1.4,
+    "figure.dpi": 150,
+    "savefig.bbox": "tight",
+})
+
+D = np.load(os.path.join(SCRATCH, "docdata.npz"))
+FX = np.load(os.path.join(SCRATCH, "fixdata.npz"))
+H = np.load(os.path.join(SCRATCH, "humpdata.npz"))
+HW = {t: np.load(os.path.join(SCRATCH, f"humpW_{t}.npz")) for t in ("lo", "at", "hi")}
+H2 = np.load(os.path.join(SCRATCH, "healthy2.npz"))
+A = np.load(os.path.join(SCRATCH, "branch_pi.npz"))
+B = np.load(os.path.join(SCRATCH, "branch_fix_pi.npz"))
+
+
+def cobweb(ax, xs, color=INK2, lw=0.8, zorder=5, arrow_last=False):
+    """Draw cobweb path for iterates xs on an (x, g(x)) axes."""
+    for k in range(len(xs) - 1):
+        x0, x1 = xs[k], xs[k + 1]
+        ax.plot([x0, x0], [x0, x1], color=color, lw=lw, zorder=zorder, alpha=0.85)
+        ax.plot([x0, x1], [x1, x1], color=color, lw=lw, zorder=zorder, alpha=0.85)
+
+
+def plot_map_with_jumps(ax, pm, g, jump_thresh=0.1, color=BLUE, **kw):
+    """Plot g with NaN breaks at jumps so no vertical connector is drawn."""
+    gg = g.copy()
+    br = np.where(np.abs(np.diff(g)) > jump_thresh)[0]
+    gplot = gg.astype(float)
+    pmp = pm.astype(float)
+    for q in br:
+        gplot = np.insert(gplot, q + 1, np.nan)
+        pmp = np.insert(pmp, q + 1, np.nan)
+    ax.plot(pmp, gplot, color=color, lw=kw.pop("lw", 1.4),
+            zorder=kw.pop("zorder", 4), **kw)
+    return br
+
+
+# =====================================================================
+# FIG toy map: abstract gap map with superattracting 2-cycle
+# =====================================================================
+fig, ax = plt.subplots(figsize=(3.4, 3.2))
+xstar = 0.5
+s = 0.25
+a1, delta = 0.62, 0.28
+xl = np.linspace(0.05, xstar, 200)
+xr = np.linspace(xstar, 0.95, 200)
+gl = a1 - s * (xl - xstar)
+gr = a1 - delta - s * (xr - xstar)
+ax.plot(xl, gl, color=BLUE, lw=1.6)
+ax.plot(xr, gr, color=BLUE, lw=1.6)
+ax.plot([xstar, xstar], [a1 - delta, a1], color=BLUE, lw=0.8, ls=":", alpha=0.6)
+lim = [0.05, 0.95]
+ax.plot(lim, lim, color=INK2, lw=0.8, ls="--")
+# iterate from 0.15
+xs = [0.12]
+for _ in range(9):
+    x = xs[-1]
+    xs.append(a1 - s * (x - xstar) if x <= xstar else a1 - delta - s * (x - xstar))
+cobweb(ax, xs, color=INK, lw=0.9)
+ax.scatter(xs[-2:], xs[-1:] + [xs[-1]], s=0)  # noop keep
+ax.set_xlim(*lim)
+ax.set_ylim(*lim)
+ax.set_xlabel(r"$x$")
+ax.set_ylabel(r"$g(x)$")
+ax.set_title("gap map: no fixed point,\nsuperattracting 2-cycle")
+ax.annotate("gap straddles\nthe diagonal", xy=(xstar, a1 - delta / 2),
+            xytext=(0.6, 0.35), fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+fig.savefig(f"{FIG}/fig_toymap.pdf"); fig.savefig(f"{FIG}/fig_toymap.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 1: healthy anatomy — b(p) and W(p_t), saturated parcel
+# =====================================================================
+fig, axs = plt.subplots(1, 2, figsize=(7.0, 3.1))
+ax = axs[0]
+ax.plot(H2["c_b"], H2["c_blevels"], color=BLUE, marker="o", ms=2.5, lw=1.2)
+ax.axvline(0, color=INK2, lw=0.7)
+ax.set_ylim(1010, 40)
+ax.set_xlabel(r"buoyancy $b_j$  [K]")
+ax.set_ylabel("pressure [hPa]  (up = higher altitude)")
+ax.set_title("(a) buoyancy of the lifted parcel")
+ax = axs[1]
+ax.plot(H2["c_W_W"], H2["c_W_p"], color=BLUE, lw=1.4)
+kmax = int(np.argmax(H2["c_W_W"]))
+ax.plot(H2["c_W_W"][kmax], H2["c_W_p"][kmax], "o", color=AQUA, ms=6, zorder=6)
+ax.annotate("LNB / argmax\n(CAPE = max $W$)",
+            xy=(H2["c_W_W"][kmax], H2["c_W_p"][kmax]), xytext=(0.35, 0.45),
+            textcoords="axes fraction", fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+ax.axvline(0, color=INK2, lw=0.7)
+ax.set_ylim(1010, 40)
+ax.set_xlabel(r"$W(p_t)$  [J kg$^{-1}$]")
+ax.set_title("(b) running work integral (single hump)")
+fig.suptitle("Generic column (January, eastern tropical Pacific, 10\u00b0N 108\u00b0W): one positive region, unambiguous LNB",
+             fontsize=9, y=1.02)
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_healthy_anatomy.pdf"); fig.savefig(f"{FIG}/fig_healthy_anatomy.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 2: healthy map — g smooth at the crossing, iteration converges
+# =====================================================================
+fig, ax = plt.subplots(figsize=(4.6, 3.4))
+plot_map_with_jumps(ax, H2["c_pm"], H2["c_g"], jump_thresh=0.15)
+lim = [900, 1000]
+ax.plot(lim, lim, color=INK2, lw=0.8, ls="--", label="diagonal $g(x)=x$")
+cobweb(ax, list(H2["c_iter"]), color=INK, lw=0.9)
+xfp = H2["c_iter"][-1]
+ax.plot(xfp, xfp, "o", color=AQUA, ms=7, zorder=7)
+ax.set_xlim(*lim)
+ax.set_ylim(*lim)
+ax.set_xlabel(r"$x=$ trial central pressure $P_M$ [hPa]")
+ax.set_ylabel(r"$g(x)$ [hPa]")
+ax.set_title("Generic column: contraction to a fixed point")
+ax.legend(loc="lower right")
+fig.savefig(f"{FIG}/fig_healthy_map.pdf"); fig.savefig(f"{FIG}/fig_healthy_map.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 3: PR77 map, full + zoom with 2-cycle
+# =====================================================================
+fig, axs = plt.subplots(1, 2, figsize=(8.6, 3.4))
+ax = axs[0]
+br = plot_map_with_jumps(ax, D["a_pm"], D["a_g"], jump_thresh=0.3)
+lim = [946, 958]
+ax.plot(lim, lim, color=INK2, lw=0.8, ls="--")
+cobweb(ax, list(D["a_iter"][:14]), color=INK, lw=0.8)
+ax.set_xlim(*lim)
+ax.set_ylim(949.8, 953.2)
+ax.set_xlabel(r"$x = P_M$ [hPa]")
+ax.set_ylabel(r"$g(x)$ [hPa]")
+ax.set_title("(a) issue-77 column: the map $g$")
+zx0, zx1 = 950.45, 951.5
+ax.add_patch(plt.Rectangle((zx0, zx0), zx1 - zx0, zx1 - zx0, fill=False,
+                           edgecolor=RED, lw=0.9, zorder=7))
+ax = axs[1]
+m = (D["a_pm"] >= zx0) & (D["a_pm"] <= zx1)
+plot_map_with_jumps(ax, D["a_pm"][m], D["a_g"][m], jump_thresh=0.3)
+ax.plot([zx0, zx1], [zx0, zx1], color=INK2, lw=0.8, ls="--")
+cyc = [950.6533454438, 951.2790079839]
+cobweb(ax, list(D["a_iter"][6:]) + cyc * 3, color=INK, lw=0.9)
+ax.plot(cyc, cyc[::-1], "s", color=RED, ms=5, zorder=8, label="2-cycle $\\{A,B\\}$")
+jump_x = 950.829
+ax.annotate("gap: $g(x)-x$ changes sign\nhere, but $g$ has no value\non the diagonal",
+            xy=(jump_x, jump_x + 0.02), xytext=(950.95, 950.62), fontsize=8,
+            color=INK2, arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+ax.set_xlim(zx0, zx1)
+ax.set_ylim(zx0, zx1)
+ax.set_xlabel(r"$x = P_M$ [hPa]")
+ax.set_title("(b) zoom: the diagonal crosses inside the gap")
+ax.legend(loc="upper right")
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_pr77_map.pdf"); fig.savefig(f"{FIG}/fig_pr77_map.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 4: PR77 mechanism — buoyancy zoom + W curves at both phases
+# =====================================================================
+fig, axs = plt.subplots(1, 2, figsize=(8.2, 3.4))
+ax = axs[0]
+lv = D["a_blevels"]
+ax.plot(D["a_b_hi"], lv, color=YELLOW, marker="o", ms=3, lw=1.2,
+        label=r"$x=B=951.279$ (bump $<0$)")
+ax.plot(D["a_b_lo"], lv, color=BLUE, marker="o", ms=3, lw=1.2,
+        label=r"$x=A=950.653$ (bump $>0$)")
+ax.axvline(0, color=INK2, lw=0.7)
+ax.set_ylim(1010, 240)
+ax.set_xlim(-4.0, 2.6)
+ax.set_xlabel(r"buoyancy $b_j$ [K]")
+ax.set_ylabel("pressure [hPa]")
+ax.set_title("(a) eyewall-parcel buoyancy at the two phases")
+ax.legend(loc="lower left", fontsize=7.5)
+axz = ax.inset_axes([0.60, 0.52, 0.37, 0.40])
+axz.plot(D["a_b_hi"], lv, color=YELLOW, marker="o", ms=3, lw=1.0)
+axz.plot(D["a_b_lo"], lv, color=BLUE, marker="o", ms=3, lw=1.0)
+axz.axvline(0, color=INK2, lw=0.6)
+axz.set_xlim(-0.08, 0.05)
+axz.set_ylim(410, 290)
+axz.set_title("350 hPa bump: $\\pm$14 mK", fontsize=7)
+axz.tick_params(labelsize=6)
+axz.grid(True, color=GRID, lw=0.5)
+
+ax = axs[1]
+ax.plot(D["a_W_hi_W"], D["a_W_hi_p"], color=YELLOW, lw=1.4, label=r"$W$ at $x=B$")
+ax.plot(D["a_W_lo_W"], D["a_W_lo_p"], color=BLUE, lw=1.4, label=r"$W$ at $x=A$")
+ax.axvline(0, color=INK2, lw=0.7)
+# selections
+khi = int(np.argmax(D["a_W_hi_W"]))
+ax.plot(D["a_W_hi_W"][khi], D["a_W_hi_p"][khi], "o", color=AQUA, ms=6, zorder=6)
+# pcmin selection at A: topmost crossing ~ level 18/19 where W ~ -61
+plo = D["a_W_lo_p"]
+wlo = D["a_W_lo_W"]
+ksel = np.argmin(np.abs(plo - 349.0))
+ax.plot(wlo[ksel], plo[ksel], "X", color=RED, ms=8, zorder=7)
+ax.annotate("code's selection at $x=A$:\n$W=-61$ J/kg $\\to$ clamped to 0",
+            xy=(wlo[ksel], plo[ksel]), xytext=(0.25, 0.30),
+            textcoords="axes fraction",
+            fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+ax.annotate("argmax (both phases):\n$W\\approx 222$ J/kg at 745 hPa",
+            xy=(np.max(D["a_W_hi_W"]), D["a_W_hi_p"][int(np.argmax(D['a_W_hi_W']))]),
+            xytext=(0.05, 0.72), textcoords="axes fraction",
+            fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+ax.set_ylim(1010, 240)
+ax.set_xlim(-350, 330)
+ax.set_xlabel(r"$W(p_t)$ [J kg$^{-1}$]")
+ax.set_title("(b) work integral and the two selections")
+ax.legend(loc="lower left", fontsize=7.5)
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_pr77_mechanism.pdf"); fig.savefig(f"{FIG}/fig_pr77_mechanism.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 5: CAPEM(PM) staircase + zoom on g
+# =====================================================================
+fig, axs = plt.subplots(1, 2, figsize=(8.2, 3.1))
+ax = axs[0]
+cm = D["a_capem"].copy()
+pm = D["a_pm"]
+br = np.where(np.abs(np.diff(cm)) > 30)[0]
+cmp_, pmp = cm.astype(float), pm.astype(float)
+for q in br[::-1]:
+    cmp_ = np.insert(cmp_, q + 1, np.nan)
+    pmp = np.insert(pmp, q + 1, np.nan)
+ax.plot(pmp, cmp_, color=BLUE, lw=1.4)
+ax.set_xlabel(r"$x = P_M$ [hPa]")
+ax.set_ylabel(r"CAPE$_M(x)$ [J kg$^{-1}$]")
+ax.set_title("(a) the discontinuous coefficient: CAPE$_M$")
+ax.annotate("222 J/kg jump from a\n$10^{-4}$ K sign flicker",
+            xy=(950.83, 100), xytext=(952.5, 60), fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+ax.set_xlim(946, 958)
+
+ax = axs[1]
+mm = (pm > 950.5) & (pm < 951.2)
+ax.plot(pm[mm], D["a_g"][mm] - pm[mm], color=BLUE, lw=1.4)
+ax.axhline(0, color=INK2, lw=0.8, ls="--")
+ax.set_xlabel(r"$x = P_M$ [hPa]")
+ax.set_ylabel(r"$F(x) = g(x) - x$ [hPa]")
+ax.set_title("(b) residual changes sign only across the gap")
+ax.annotate("$F>0$", xy=(950.62, 0.45), fontsize=9, color=INK2)
+ax.annotate("$F<0$", xy=(951.0, -0.28), fontsize=9, color=INK2)
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_pr77_capem.pdf"); fig.savefig(f"{FIG}/fig_pr77_capem.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 6: 1995 case map
+# =====================================================================
+fig, ax = plt.subplots(figsize=(4.6, 3.4))
+plot_map_with_jumps(ax, D["b_pm"], D["b_g"], jump_thresh=0.3)
+lim = [955, 968]
+ax.plot(lim, lim, color=INK2, lw=0.8, ls="--")
+cobweb(ax, list(D["b_iter"][4:20]), color=INK, lw=0.8)
+ax.set_xlim(959, 964)
+ax.set_ylim(959, 964)
+ax.set_xlabel(r"$x = P_M$ [hPa]")
+ax.set_ylabel(r"$g(x)$ [hPa]")
+ax.set_title("Second instance (real 1995 sounding):\nsame gap-straddling geometry")
+fig.savefig(f"{FIG}/fig_1995_map.pdf"); fig.savefig(f"{FIG}/fig_1995_map.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 7: survey + SST sweep
+# =====================================================================
+fig, axs = plt.subplots(1, 2, figsize=(8.4, 3.1))
+ax = axs[0]
+mj = D["s_maxjump"]
+mj = mj[np.isfinite(mj)]
+ax.hist(np.log10(mj), bins=28, color=BLUE, edgecolor="white", lw=0.5)
+ax.axvline(np.log10(0.5), color=RED, lw=1.2, ls="--")
+ax.annotate("convergence\ntolerance 0.5 hPa", xy=(np.log10(0.5), ax.get_ylim()[1] * 0.75),
+            xytext=(0.62, 0.8), textcoords="axes fraction", fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+n_all = int(D["s_n"])
+ax.set_xlabel(r"$\log_{10}$ largest jump of $g$ in column [hPa]")
+ax.set_ylabel("columns")
+ax.set_title(f"(a) jumps are generic: {len(mj)}/{n_all} sampled\ncolumns have $\\geq$1 jump; 0 fail")
+ax = axs[1]
+st = D["d_status"]
+pmfp = D["d_pm"]
+dss = D["d_dss"]
+ax.plot(dss, pmfp, color=BLUE, lw=1.0)
+bad = st == 2
+ax.axvspan(dss[bad].min() - 0.005, dss[bad].max() + 0.005, color=RED,
+           alpha=0.75, lw=0)
+ax.annotate("the Francine (issue-77) member ($\\Delta=0$) sits\nin a $\\sim$0.04 K-wide no-fixed-point band",
+            xy=(0.0, 950.5), xytext=(-2.9, 925), fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=RED, lw=1.0))
+ax.set_xlabel(r"further SST perturbation $\Delta$ [K]")
+ax.set_ylabel(r"fixed point $x^\ast$ [hPa]")
+ax.set_title("(b) failure band under SST perturbation\n(red: no fixed point, 4/601 offsets)")
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_survey.pdf"); fig.savefig(f"{FIG}/fig_survey.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 8: genuine argmax bifurcation (Berge picture), real column
+# =====================================================================
+fig = plt.figure(figsize=(8.6, 3.5))
+gs = fig.add_gridspec(2, 2, width_ratios=[1, 1.15], hspace=0.12)
+ax = fig.add_subplot(gs[:, 0])
+for tag, cc in [("lo", BLUE), ("at", AQUA), ("hi", YELLOW)]:
+    d = HW[tag]
+    ax.plot(d["W"], d["p"], color=cc, lw=1.3, label=f"$f={float(d['f']):.3f}$")
+ax.axvline(0, color=INK2, lw=0.7)
+ax.set_ylim(1010, 130)
+ax.set_xlim(-260, 320)
+ax.set_xlabel(r"$W(p_t)$ [J kg$^{-1}$]")
+ax.set_ylabel("pressure [hPa]")
+ax.set_title("(a) two competing humps as moisture scales")
+ax.legend(loc="lower left", title="moisture scale $f$", fontsize=7.5,
+          title_fontsize=7.5)
+fs = H["fs"]
+axt = fig.add_subplot(gs[0, 1])
+axt.plot(fs, H["cape_mw"], color=AQUA, lw=1.8)
+axt.set_ylabel(r"CAPE [J kg$^{-1}$]")
+axt.set_title("(b) value continuous ...")
+axt.tick_params(labelbottom=False)
+axb = fig.add_subplot(gs[1, 1], sharex=axt)
+axb.plot(fs, H["lnb_mw"], color=VIOLET, lw=1.3)
+axb.invert_yaxis()
+axb.set_ylabel("argmax [hPa]")
+axb.set_xlabel(r"boundary-layer moisture scale $f$")
+axb.set_title("(c) ... while the maximizer jumps", pad=2)
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_berge.pdf"); fig.savefig(f"{FIG}/fig_berge.png", dpi=100)
+plt.close(fig)
+
+# =====================================================================
+# FIG 9: the repair — g_fix continuous + population effect
+# =====================================================================
+fig, axs = plt.subplots(1, 2, figsize=(8.4, 3.3))
+ax = axs[0]
+m = (D["a_pm"] > 949.5) & (D["a_pm"] < 953)
+plot_map_with_jumps(ax, D["a_pm"][m], D["a_g"][m], jump_thresh=0.3, color=RED,
+                    lw=1.6, label=r"legacy map ($E_{\rm top}$ CAPE): gap, no fixed point")
+mf = (FX["pm"] > 949.5) & (FX["pm"] < 953)
+ax.plot(FX["pm"][mf], FX["g"][mf], color=AQUA, lw=1.6, ls=(0, (4, 3)), zorder=6,
+        label=r"repaired map ($E_{\max}$ CAPE): continuous")
+ax.plot([949.5, 953], [949.5, 953], color=INK2, lw=0.8, ls="--")
+cobweb(ax, list(FX["iters"][2:]), color=INK, lw=0.9)
+xf = FX["iters"][-1]
+ax.plot(xf, xf, "o", color=AQUA, mec=INK, mew=0.6, ms=7, zorder=8)
+ax.set_xlim(949.9, 952.4)
+ax.set_ylim(949.9, 952.4)
+ax.set_xlabel(r"$x = P_M$ [hPa]")
+ax.set_ylabel(r"$g(x)$ [hPa]")
+ax.set_title("(a) same map, one changed ingredient:\n$E_{\max}$ CAPE closes the gap")
+ax.legend(loc="upper left", fontsize=7.5)
+
+ax = axs[1]
+both1 = (A["IFL"] == 1) & (B["IFL"] == 1)
+d = np.abs(A["VMAX"] - B["VMAX"])[both1]
+d = np.sort(d[np.isfinite(d)])
+dnz = np.clip(d, 1e-12, None)
+ecdf = np.arange(1, len(d) + 1) / len(d)
+ax.semilogx(dnz, 100 * ecdf, color=BLUE, lw=1.6)
+ax.axvline(0.1, color=RED, lw=1.0, ls="--")
+ax.annotate("98.1% below 0.1 m/s\n(loop-tolerance noise)", xy=(0.1, 55),
+            xytext=(1.5e-8, 40), fontsize=8, color=INK2,
+            arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
+ax.set_xlabel(r"$|\Delta V_{\max}|$ pcmin vs max-$W$  [m s$^{-1}$]")
+ax.set_ylabel("cumulative share of columns [%]")
+ax.set_title("(b) population effect over 7093 columns:\nflags unchanged, tail = fragile columns")
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_fix.pdf"); fig.savefig(f"{FIG}/fig_fix.png", dpi=100)
+plt.close(fig)
+
+print("all figures written to", FIG)
diff --git a/discontinuity_analysis/scripts/plots2.py b/discontinuity_analysis/scripts/plots2.py
new file mode 100644
index 0000000..9510618
--- /dev/null
+++ b/discontinuity_analysis/scripts/plots2.py
@@ -0,0 +1,213 @@
+"""Redesigned toy figure + discretization figure."""
+
+import os
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+FIG = "/home/mares/repos/tcpyPI/discontinuity_analysis/figures"
+
+BLUE = "#2a78d6"
+AQUA = "#1baf7a"
+YELLOW = "#eda100"
+VIOLET = "#4a3aa7"
+RED = "#e34948"
+INK = "#0b0b0b"
+INK2 = "#52514e"
+GRID = "#e5e4e0"
+
+plt.rcParams.update({
+    "font.size": 9, "font.family": "serif",
+    "axes.edgecolor": INK2, "axes.labelcolor": INK,
+    "axes.titlesize": 9.5, "axes.titleweight": "bold",
+    "axes.spines.top": False, "axes.spines.right": False,
+    "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6,
+    "xtick.color": INK2, "ytick.color": INK2,
+    "legend.frameon": False, "legend.fontsize": 8,
+    "lines.linewidth": 1.4, "figure.dpi": 150, "savefig.bbox": "tight",
+})
+
+tgrid = np.linspace(0, 1.1, 4401)
+dt = tgrid[1] - tgrid[0]
+
+
+def G(c, s):
+    return np.exp(-((tgrid - c) / s) ** 2)
+
+
+def b_toy(eps, delta, h2=1.0):
+    return 0.6 * G(0.22, 0.10) - delta * G(0.5, 0.12) + eps * G(0.8, 0.06) \
+        + (h2 - 1.0) * 0.0
+
+
+def b_two(delta, a2=1.0):
+    return 0.6 * G(0.22, 0.10) - delta * G(0.5, 0.12) + a2 * G(0.8, 0.10)
+
+
+def conv(b):
+    """(E_top, E_max, E_reach, t_argmax) for buoyancy array b on tgrid."""
+    W = np.concatenate([[0.0], np.cumsum(0.5 * (b[1:] + b[:-1]) * dt)])
+    sgn = np.sign(b)
+    cross = np.where(np.diff(sgn) != 0)[0]
+    down = [q for q in cross if b[q] > 0 >= b[q + 1]]
+    # pcmin convention: terminal point = crossing above the TOPMOST positive
+    # point; if b > 0 at the domain top, integrate to the top (INB = last level)
+    if b[-1] > 0:
+        E_top = max(W[-1], 0.0)
+    else:
+        E_top = max(W[down[-1]], 0.0) if down else 0.0
+    kmax = int(np.argmax(W))
+    E_max = max(W[kmax], 0.0)
+    pos = W > 1e-12
+    if not pos.any():
+        E_reach = 0.0
+    else:
+        k0 = int(np.argmax(pos))
+        k1 = k0
+        while k1 < len(W) - 1 and W[k1 + 1] > 0:
+            k1 += 1
+        E_reach = float(W[k0:k1 + 1].max())
+    return E_top, E_max, E_reach, tgrid[kmax]
+
+
+# ---------------- fig_toy: 2x2, two sweeps ----------------
+fig, axs = plt.subplots(2, 2, figsize=(8.6, 6.2))
+
+# (a) epsilon family, delta=1.1 fixed
+ax = axs[0, 0]
+for eps, cc in [(0.05, BLUE), (0.6, YELLOW), (1.3, RED)]:
+    ax.plot(tgrid, 0.6 * G(0.22, 0.10) - 1.1 * G(0.5, 0.12) + eps * G(0.8, 0.06),
+            color=cc, label=rf"$\varepsilon={eps}$")
+ax.axhline(0, color=INK2, lw=0.7)
+ax.set_xlabel(r"height $t$")
+ax.set_ylabel(r"$b(t)$")
+ax.set_title(r"(a) sweep 1: marginal upper bump, amplitude $\varepsilon$")
+ax.legend(loc="upper left")
+
+# (b) value functions vs epsilon
+ax = axs[0, 1]
+epss = np.linspace(-0.1, 1.5, 1601)
+V = np.array([conv(0.6 * G(0.22, 0.10) - 1.1 * G(0.5, 0.12) + e * G(0.8, 0.06))
+              for e in epss])
+ax.plot(epss, V[:, 0], color=RED, lw=1.5, label=r"$E_{\rm top}$ (code's rule)")
+ax.plot(epss, V[:, 1], color=AQUA, lw=1.9, label=r"$E_{\max}$")
+ax.plot(epss, V[:, 2], color=VIOLET, ls="--", lw=1.3, label=r"$E_{\rm reach}$")
+ax.set_xlabel(r"bump amplitude $\varepsilon$")
+ax.set_ylabel("available energy")
+ax.set_title(r"(b) $E_{\rm top}$ jumps at $\varepsilon\!\approx\!0^+$ (artificial)")
+ax.legend(loc="center right")
+jt = np.argmax(np.abs(np.diff(V[:, 0])))
+print(f"toy sweep1: E_top jump at eps={epss[jt]:.4f}, size={V[jt,0]-V[jt+1,0]:+.4f}")
+print(f"  E_max kink near eps where argmax jumps: "
+      f"{epss[np.argmax(np.abs(np.diff(V[:,3])))]:.3f}")
+
+# (c) delta family
+ax = axs[1, 0]
+for delta, cc in [(0.3, BLUE), (0.75, YELLOW), (1.3, RED)]:
+    ax.plot(tgrid, b_two(delta), color=cc, label=rf"$\delta={delta}$")
+ax.axhline(0, color=INK2, lw=0.7)
+ax.set_xlabel(r"height $t$")
+ax.set_ylabel(r"$b(t)$")
+ax.set_title(r"(c) sweep 2: two real cells, valley depth $\delta$")
+ax.legend(loc="lower left")
+
+# (d) value functions vs delta
+ax = axs[1, 1]
+deltas = np.linspace(0.2, 1.8, 1601)
+V2 = np.array([conv(b_two(d)) for d in deltas])
+ax.plot(deltas, V2[:, 0], color=RED, lw=1.5, label=r"$E_{\rm top}$")
+ax.plot(deltas, V2[:, 1], color=AQUA, lw=1.9, label=r"$E_{\max}$")
+ax.plot(deltas, V2[:, 2], color=VIOLET, ls="--", lw=1.3, label=r"$E_{\rm reach}$")
+ax.set_xlabel(r"valley depth $\delta$")
+ax.set_ylabel("available energy")
+ax.set_title(r"(d) $E_{\rm reach}$ jumps at barrier closure (genuine)")
+ax.legend(loc="upper right")
+jr = np.argmax(np.abs(np.diff(V2[:, 2])))
+print(f"toy sweep2: E_reach jump at delta={deltas[jr]:.4f}, size={V2[jr,2]-V2[jr+1,2]:+.4f}")
+jm = np.argmax(np.abs(np.diff(V2[:, 3])))
+print(f"  E_max argmax jump at delta={deltas[jm]:.4f} "
+      f"(t: {V2[jm,3]:.2f}->{V2[jm+1,3]:.2f}); E_max step {V2[jm,1]-V2[jm+1,1]:+.5f}")
+
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_toy.pdf"); fig.savefig(f"{FIG}/fig_toy.png", dpi=100)
+plt.close(fig)
+
+# ---------------- fig_discretization ----------------
+R = np.load(os.path.join(SCRATCH, "refinedata.npz"))
+
+fig, axs = plt.subplots(1, 3, figsize=(10.2, 3.2))
+
+# (a) toy: truth vs sampled/PWL (two node phases) vs smooth recovery
+ax = axs[0]
+tt = np.linspace(0, 1, 2001)
+
+
+def truth(t):
+    return (-0.06 - 0.03 * (t - 0.5)
+            + 0.10 * np.exp(-((t - 0.53) / 0.015) ** 2))
+
+
+ax.plot(tt, truth(tt), color=INK, lw=1.6, label="ground truth")
+for off, cc, lbl in [(0.00, BLUE, "nodes miss the bump"),
+                     (0.0355, YELLOW, "node lands on the bump")]:
+    nodes = np.arange(0.08, 1.0, 0.0833) + off
+    nodes = nodes[nodes <= 1]
+    ax.plot(nodes, truth(nodes), "o", color=cc, ms=4)
+    ax.plot(nodes, truth(nodes), color=cc, lw=1.1, ls="-", alpha=0.9, label=lbl)
+ax.axhline(0, color=INK2, lw=0.7)
+ax.set_xlabel(r"height $t$")
+ax.set_ylabel(r"$b(t)$")
+ax.set_title("(a) sampling butchers a narrow feature:\nvisibility depends on node phase")
+ax.legend(loc="upper left", fontsize=7.5)
+
+# (b) real column: native vs reconstructions, zoom
+ax = axs[1]
+ax.plot(R["native_b"], R["native_blev"], "o-", color=INK, ms=3.5, lw=1.3,
+        label="native 37-level (PWL)")
+ax.plot(R["linear_b"], R["linear_blev"], color=YELLOW, lw=1.1,
+        label="linear-in-$\\log p$, 5 hPa")
+ax.plot(R["pchip_b"], R["pchip_blev"], color=AQUA, lw=1.4,
+        label="monotone cubic, 5 hPa")
+ax.axvline(0, color=INK2, lw=0.7)
+ax.set_ylim(520, 240)
+ax.set_xlim(-0.6, 0.25)
+ax.set_xlabel(r"buoyancy $b$ [K]")
+ax.set_ylabel("pressure [hPa]")
+ax.set_title("(b) issue-77 column, reconstructed:\nthe bump is real but reconstruction-dependent")
+ax.legend(loc="lower left", fontsize=7.5)
+
+# (c) residuals F = g - x for native / linear / pchip
+ax = axs[2]
+Dn = np.load(os.path.join(SCRATCH, "docdata.npz"))
+mm = (Dn["a_pm"] > 950.3) & (Dn["a_pm"] < 952.0)
+
+
+def plot_res(ax, pm, g, color, label, jump_thresh=0.25):
+    F = (g - pm).astype(float)
+    pmp = pm.astype(float)
+    br = np.where(np.abs(np.diff(F)) > jump_thresh)[0]
+    for q in br[::-1]:
+        F = np.insert(F, q + 1, np.nan)
+        pmp = np.insert(pmp, q + 1, np.nan)
+    ax.plot(pmp, F, color=color, lw=1.4, label=label)
+
+
+plot_res(ax, Dn["a_pm"][mm], Dn["a_g"][mm], INK, "native (fails)")
+ml = (R["linear_pm"] > 950.3) & (R["linear_pm"] < 952.0)
+plot_res(ax, R["linear_pm"][ml], R["linear_g"][ml], YELLOW, "linear 5 hPa (fails)")
+mp = (R["pchip_pm"] > 950.3) & (R["pchip_pm"] < 952.0)
+plot_res(ax, R["pchip_pm"][mp], R["pchip_g"][mp], AQUA, "cubic 5 hPa (converges)")
+ax.axhline(0, color=INK2, lw=0.8, ls="--")
+ax.set_xlabel(r"$x = P_M$ [hPa]")
+ax.set_ylabel(r"$F(x) = g(x) - x$ [hPa]")
+ax.set_title("(c) recovery moves the gap but cannot\nremove it: zero crossing vs gap position")
+ax.legend(loc="upper right", fontsize=7.5)
+
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_discretization.pdf"); fig.savefig(f"{FIG}/fig_discretization.png", dpi=100)
+plt.close(fig)
+print("saved fig_toy.pdf (redesigned), fig_discretization.pdf")
diff --git a/discontinuity_analysis/scripts/plots3.py b/discontinuity_analysis/scripts/plots3.py
new file mode 100644
index 0000000..2893953
--- /dev/null
+++ b/discontinuity_analysis/scripts/plots3.py
@@ -0,0 +1,115 @@
+"""Four-archetype figure: the three CAPE conventions on concrete profiles."""
+
+import os
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+FIG = "/home/mares/repos/tcpyPI/discontinuity_analysis/figures"
+
+BLUE = "#2a78d6"
+AQUA = "#1baf7a"
+VIOLET = "#4a3aa7"
+RED = "#e34948"
+INK = "#0b0b0b"
+INK2 = "#52514e"
+GRID = "#e5e4e0"
+
+plt.rcParams.update({
+    "font.size": 9, "font.family": "serif",
+    "axes.edgecolor": INK2, "axes.labelcolor": INK,
+    "axes.titlesize": 9.5, "axes.titleweight": "bold",
+    "axes.spines.top": False, "axes.spines.right": False,
+    "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6,
+    "xtick.color": INK2, "ytick.color": INK2,
+    "legend.frameon": False, "legend.fontsize": 7.5,
+    "lines.linewidth": 1.4, "figure.dpi": 150, "savefig.bbox": "tight",
+})
+
+t = np.linspace(0, 1.15, 4601)
+dt = t[1] - t[0]
+
+
+def G(c, s):
+    return np.exp(-((t - c) / s) ** 2)
+
+
+def conventions(b):
+    """Return dict with values and terminal locations for E_top/E_max/E_reach."""
+    W = np.concatenate([[0.0], np.cumsum(0.5 * (b[1:] + b[:-1]) * dt)])
+    sgn = np.sign(b)
+    cross = np.where(np.diff(sgn) != 0)[0]
+    down = [q for q in cross if b[q] > 0 >= b[q + 1]]
+    if b[-1] > 0:
+        z_top, W_top = t[-1], W[-1]
+    elif down:
+        z_top, W_top = t[down[-1]], W[down[-1]]
+    else:
+        z_top, W_top = np.nan, 0.0
+    E_top = max(W_top, 0.0)
+    kmax = int(np.argmax(W))
+    E_max, z_max = max(W[kmax], 0.0), t[kmax]
+    pos = W > 1e-12
+    if not pos.any():
+        E_reach, z_reach = 0.0, np.nan
+    else:
+        k0 = int(np.argmax(pos))
+        k1 = k0
+        while k1 < len(W) - 1 and W[k1 + 1] > 0:
+            k1 += 1
+        kr = k0 + int(np.argmax(W[k0:k1 + 1]))
+        E_reach, z_reach = float(W[kr]), t[kr]
+    return dict(W=W, E_top=E_top, z_top=z_top, W_top=W_top,
+                E_max=E_max, z_max=z_max, E_reach=E_reach, z_reach=z_reach)
+
+
+CASES = [
+    ("(a) single positive region:\nall three agree",
+     1.0 * G(0.35, 0.18) - 0.7 * G(0.85, 0.12)),
+    ("(b) dominated upper cell, no barrier:\n$E_{\\rm top}$ undercounts",
+     1.0 * G(0.25, 0.10) - 0.55 * G(0.5, 0.09) + 0.35 * G(0.72, 0.08)
+     - 0.8 * G(1.0, 0.10)),
+    ("(c) marginal bump behind a barrier\n(the issue-77 shape): $E_{\\rm top}$ clamps to 0",
+     1.0 * G(0.25, 0.10) - 1.1 * G(0.5, 0.10) + 0.12 * G(0.78, 0.05)
+     - 0.5 * G(1.0, 0.08)),
+    ("(d) dominant upper cell behind a barrier:\n$E_{\\rm reach}$ stalls below the others",
+     0.6 * G(0.22, 0.09) - 1.0 * G(0.45, 0.10) + 1.3 * G(0.75, 0.12)
+     - 0.9 * G(1.05, 0.08)),
+]
+
+fig, axs = plt.subplots(2, 2, figsize=(8.8, 5.8))
+for ax, (title, b) in zip(axs.flat, CASES):
+    c = conventions(b)
+    W = c["W"]
+    ylo, yhi = W.min() - 0.03, W.max() + 0.055
+    ax.fill_between(t, ylo, yhi, where=b > 0, color=BLUE, alpha=0.07, lw=0)
+    ax.plot(t, W, color=BLUE, lw=1.6)
+    ax.axhline(0, color=INK2, lw=0.7)
+    # markers: E_max dot, E_reach open square, E_top X (at its terminal point)
+    ax.plot(c["z_max"], c["E_max"], "o", color=AQUA, ms=8, zorder=6,
+            label=rf"$E_{{\max}}={c['E_max']:.3f}$")
+    ax.plot(c["z_reach"], c["E_reach"], "s", mfc="none", mec=VIOLET, mew=1.6,
+            ms=12, zorder=7, label=rf"$E_{{\rm reach}}={c['E_reach']:.3f}$")
+    lbl = rf"$E_{{\rm top}}={c['E_top']:.3f}$"
+    if c["E_top"] != c["W_top"]:
+        lbl = rf"$E_{{\rm top}}=0$ (clamped from ${c['W_top']:.3f}$)"
+    ax.plot(c["z_top"], c["W_top"], "X", color=RED, ms=9, zorder=8, label=lbl)
+    ax.set_ylim(ylo, yhi)
+    ax.set_title(title)
+    ax.legend(loc="lower left", handletextpad=0.4)
+    ax.set_xlabel(r"trial outflow height $t$")
+    ax.set_ylabel(r"$W(t)$")
+
+fig.tight_layout()
+fig.savefig(f"{FIG}/fig_conventions.pdf")
+fig.savefig(f"{FIG}/fig_conventions.png", dpi=100)
+
+for title, b in CASES:
+    c = conventions(b)
+    print(title.split(':')[0].replace('\n', ' '),
+          f"E_top={c['E_top']:.4f} (W_top={c['W_top']:.4f}) "
+          f"E_max={c['E_max']:.4f} E_reach={c['E_reach']:.4f}")
diff --git a/discontinuity_analysis/scripts/refine.py b/discontinuity_analysis/scripts/refine.py
new file mode 100644
index 0000000..1d2bfc3
--- /dev/null
+++ b/discontinuity_analysis/scripts/refine.py
@@ -0,0 +1,115 @@
+"""Discretization/recovery experiment on the PR#77 column.
+
+Reconstruct the environment (T, R as functions of log p) with a monotone cubic
+(PCHIP / Fritsch-Carlson) interpolant, resample at 5 hPa from 1000 to 55 hPa
+(plus the original coarse levels above), and re-run:
+  - the buoyancy profile of the environmental parcel at cycle point A
+  - the fixed-point map g on the refined profile (scan + iteration)
+
+Also: the same with plain linear-in-logp resampling (the PWL 'reconstruction'
+the trapezoid scheme implicitly commits to), for contrast.
+"""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.join(SCRATCH, "branch", "src"))
+sys.path.insert(0, SCRATCH)
+
+import numpy as np
+from tcpyPI import constants, utilities
+
+from datagen import (P77, TC77, R77, SST77, MSL77, prep, g_scan, iterate_raw,
+                     g_eval)
+from gmap import cape_instr
+
+EPS = constants.EPS
+
+
+def pchip_slopes(x, y):
+    """Fritsch-Carlson monotone cubic slopes."""
+    h = np.diff(x)
+    d = np.diff(y) / h
+    n = len(x)
+    m = np.zeros(n)
+    m[0] = d[0]
+    m[-1] = d[-1]
+    for k in range(1, n - 1):
+        if d[k - 1] * d[k] <= 0:
+            m[k] = 0.0
+        else:
+            w1 = 2 * h[k] + h[k - 1]
+            w2 = h[k] + 2 * h[k - 1]
+            m[k] = (w1 + w2) / (w1 / d[k - 1] + w2 / d[k])
+    return m
+
+
+def pchip_eval(x, y, xq):
+    m = pchip_slopes(x, y)
+    idx = np.clip(np.searchsorted(x, xq) - 1, 0, len(x) - 2)
+    h = x[idx + 1] - x[idx]
+    t = (xq - x[idx]) / h
+    h00 = (1 + 2 * t) * (1 - t) ** 2
+    h10 = t * (1 - t) ** 2
+    h01 = t * t * (3 - 2 * t)
+    h11 = t * t * (t - 1)
+    return (h00 * y[idx] + h10 * h * m[idx] + h01 * y[idx + 1] + h11 * h * m[idx + 1])
+
+
+def build_refined(P, TC, R, method="pchip", dp=5.0, pmin_fine=55.0):
+    """Refined profile: fine levels from P[0] down to pmin_fine, then original."""
+    x = np.log(P[::-1])  # increasing in log p ... P descending -> reverse
+    tcv = TC[::-1]
+    rv_ = R[::-1]
+    fine = np.arange(P[0], pmin_fine - 1e-9, -dp)
+    fine = np.unique(np.concatenate([fine, P[P >= pmin_fine]]))[::-1]
+    keep_above = P[P < pmin_fine]
+    Pn = np.concatenate([fine, keep_above])
+    xq = np.log(Pn[::-1])
+    if method == "pchip":
+        tq = pchip_eval(x, tcv, xq)[::-1]
+        rq = pchip_eval(x, rv_, xq)[::-1]
+    else:
+        tq = np.interp(xq, x, tcv)[::-1]
+        rq = np.interp(xq, x, rv_)[::-1]
+    return Pn, tq, rq
+
+
+save = {}
+for method in ("pchip", "linear"):
+    Pn, TCn, Rn = build_refined(P77, TC77, R77, method=method)
+    S = prep(SST77, MSL77, TCn, Rn)
+    pms = np.arange(946.0, 958.0, 0.005)
+    sc = g_scan(pms, S[0], S[1], Pn, S[2], S[3], S[4])
+    st, xs = iterate_raw(S[0], S[1], Pn, S[2], S[3], S[4], nmax=200)
+    save[f"{method}_pm"] = pms
+    save[f"{method}_g"] = sc[:, 0]
+    save[f"{method}_capem"] = sc[:, 1]
+    save[f"{method}_status"] = st
+    save[f"{method}_iters"] = xs
+    # env-parcel buoyancy at cycle point A (native cycle value)
+    pmA = 950.6533454438
+    PP = min(pmA, 1000.0)
+    RP = EPS * S[3][0] * S[1] / (PP * (EPS + S[3][0]) - S[3][0] * S[1])
+    o = cape_instr(S[2][0], RP, PP, S[2], S[3], Pn, 0, 50)
+    N = int(np.count_nonzero(Pn > 50))
+    save[f"{method}_b"] = o["TVRDIF"]
+    save[f"{method}_blev"] = Pn[:N]
+    save[f"{method}_capemA"] = o["CAPED"]
+    print(f"{method}: status={st} n_it={len(xs)} final={xs[-1]:.4f} "
+          f"CAPEM(A)={o['CAPED']:.3f} nlev={len(Pn)}")
+
+# native for reference
+S = prep(SST77, MSL77, TC77, R77)
+pmA = 950.6533454438
+PP = min(pmA, 1000.0)
+RP = EPS * S[3][0] * S[1] / (PP * (EPS + S[3][0]) - S[3][0] * S[1])
+o = cape_instr(S[2][0], RP, PP, S[2], S[3], P77, 0, 50)
+N = int(np.count_nonzero(P77 > 50))
+save["native_b"] = o["TVRDIF"]
+save["native_blev"] = P77[:N]
+print(f"native: CAPEM(A)={o['CAPED']:.3f}")
+
+np.savez(os.path.join(SCRATCH, "refinedata.npz"), **save)
+print("saved refinedata.npz")
diff --git a/discontinuity_analysis/scripts/scan.py b/discontinuity_analysis/scripts/scan.py
new file mode 100644
index 0000000..24c2c8c
--- /dev/null
+++ b/discontinuity_analysis/scripts/scan.py
@@ -0,0 +1,63 @@
+"""Fine scan of g(PM) to locate and classify discontinuities."""
+
+import numpy as np
+from gmap import gmap, CAPEA, P_FULL
+
+# coarse-to-fine: scan a wide window first
+pm_grid = np.arange(945.0, 958.0, 0.002)
+rows = []
+for pm in pm_grid:
+    o = gmap(pm)
+    rows.append((pm, o["PNEW"], o["CAPEM"], o["CAPEMS"], o["INB_M"], o["INB_S"],
+                 o["TO"], o["CAT_raw"], o["TVRDIF_M"][18], o["TVRDIF_M"][10],
+                 o["TVRDIF_M"][11], o["PLCL_M"]))
+arr = np.array(rows)
+
+pm, g, capem, capems, inbm, inbs = arr[:, 0], arr[:, 1], arr[:, 2], arr[:, 3], arr[:, 4], arr[:, 5]
+to, catraw, tv18, tv10, tv11, plclm = arr[:, 6], arr[:, 7], arr[:, 8], arr[:, 9], arr[:, 10], arr[:, 11]
+
+# find jumps in g
+dg = np.abs(np.diff(g))
+jump_idx = np.where(dg > 10 * np.median(dg) + 1e-6)[0]
+print("=== discontinuities of g on [945, 958] (grid step 0.002) ===")
+for i in jump_idx:
+    print(
+        f"PM in ({pm[i]:.4f}, {pm[i+1]:.4f}): g jumps {g[i]:.6f} -> {g[i+1]:.6f} "
+        f"(|dg|={dg[i]:.4f});  INB_M {int(inbm[i])}->{int(inbm[i+1])}, "
+        f"INB_S {int(inbs[i])}->{int(inbs[i+1])}, "
+        f"CAPEM {capem[i]:.4f}->{capem[i+1]:.4f}, "
+        f"TVRDIF_M[18] {tv18[i]:+.2e}->{tv18[i+1]:+.2e}, "
+        f"PLCL {plclm[i]:.4f}->{plclm[i+1]:.4f}"
+    )
+
+print()
+print("=== slope of g away from jumps ===")
+mask = np.ones(len(pm), bool)
+for i in jump_idx:
+    mask[max(0, i - 2) : i + 3] = False
+sl = np.diff(g)[mask[:-1] & mask[1:]] / 0.002
+print(f"median slope {np.median(sl):.6f}, max |slope| {np.max(np.abs(sl)):.6f}")
+
+# where does the diagonal fall?
+print()
+print("=== fixed-point analysis ===")
+F = g - pm
+sign_changes = np.where(np.diff(np.sign(F)) != 0)[0]
+for i in sign_changes:
+    kind = "JUMP (no root)" if i in jump_idx else "continuous crossing (root)"
+    print(f"F=g(PM)-PM changes sign in ({pm[i]:.4f},{pm[i+1]:.4f}): F {F[i]:+.4f}->{F[i+1]:+.4f}  [{kind}]")
+
+# raw signed integral in the INB=18 phase: is the clamp active?
+o = gmap(950.6533454438)
+print()
+print("=== at cycle point A=950.6533 (INB_M=18 phase) ===")
+print(f"PA={o['PA_M']:.4f} NA={o['NA_M']:.4f} PAT={o['PAT_M']:.4f} "
+      f"raw signed integral={o['PA_M'] + o['PAT_M'] - o['NA_M']:.4f} -> CAPEM={o['CAPEM']:.4f}")
+print("TVRDIF_M profile (levels 8..20):")
+for j in range(8, 21):
+    print(f"  j={j:2d} P={P_FULL[j]:6.1f} hPa  TVRDIF={o['TVRDIF_M'][j]:+.6f} K")
+o2 = gmap(951.2790079839)
+print()
+print("=== at cycle point B=951.2790 (INB_M=10 phase) ===")
+print(f"PA={o2['PA_M']:.4f} NA={o2['NA_M']:.4f} PAT={o2['PAT_M']:.4f} CAPEM={o2['CAPEM']:.4f}")
+print(f"TVRDIF_M[18] = {o2['TVRDIF_M'][18]:+.6f} K")
diff --git a/discontinuity_analysis/scripts/validate.py b/discontinuity_analysis/scripts/validate.py
new file mode 100644
index 0000000..09c8848
--- /dev/null
+++ b/discontinuity_analysis/scripts/validate.py
@@ -0,0 +1,53 @@
+"""Full-sample validation of max-W CAPE (rescue removed) vs committed outputs.
+
+Usage: python validate.py <tree>   where <tree> is 'branch' or 'branch_fix'.
+Runs pi() over every column/month of the committed 2024 sample and writes
+<tree>_pi.npz.
+"""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+tree = sys.argv[1]
+sys.path.insert(0, os.path.join(SCRATCH, tree, "src"))
+
+import numpy as np
+import xarray as xr
+from tcpyPI import pi
+
+ds = xr.open_dataset(os.path.join(SCRATCH, "sample_data.nc"))
+p = ds["p"].values.astype(float)
+
+sst = ds["sst"].values  # (month, lat, lon) or similar
+msl = ds["msl"].values
+t = ds["t"].values
+r = ds["r"].values
+
+print("dims:", ds["sst"].dims, sst.shape, "t:", ds["t"].dims, t.shape)
+
+nm, ny, nx = sst.shape
+VMAX = np.full((nm, ny, nx), np.nan)
+PMIN = np.full((nm, ny, nx), np.nan)
+IFL = np.zeros((nm, ny, nx), dtype=np.int64)
+TO = np.full((nm, ny, nx), np.nan)
+OTL = np.full((nm, ny, nx), np.nan)
+
+for m in range(nm):
+    for j in range(ny):
+        for i in range(nx):
+            if np.isnan(sst[m, j, i]) or np.isnan(msl[m, j, i]):
+                IFL[m, j, i] = 3
+                continue
+            out = pi(
+                sst[m, j, i], msl[m, j, i], p,
+                t[m, :, j, i].astype(float), r[m, :, j, i].astype(float),
+                CKCD=0.9, ascent_flag=0, diss_flag=1, V_reduc=0.8,
+                ptop=50, miss_handle=1,
+            )
+            VMAX[m, j, i], PMIN[m, j, i], IFL[m, j, i], TO[m, j, i], OTL[m, j, i] = out
+
+np.savez(os.path.join(SCRATCH, f"{tree}_pi.npz"),
+         VMAX=VMAX, PMIN=PMIN, IFL=IFL, TO=TO, OTL=OTL)
+print(f"saved {tree}_pi.npz; IFL counts:",
+      {int(k): int(v) for k, v in zip(*np.unique(IFL, return_counts=True))})