Look up a padded Zalmoxis surface node at the target pressure

src/proteus/interior_struct/zalmoxis.py

@@ -693,7 +693,9 @@ def _resolve_zalmoxis_anchor_temperatures(
     P_cmb estimate on the very first call; the next structure iteration
     recomputes it against the converged Zalmoxis P_cmb. Zalmoxis anchors its
     upward integration at the CMB value and takes the surface value for its
-    first-guess profile, its zero-pressure nodes and its nodes below the CMB.
+    first-guess profile and its nodes below the CMB. Its zero-pressure nodes
+    take the adiabat continued to the target surface pressure once the adiabat
+    blend reaches 1, and the linear profile at the outer radius before that.
     The energetics IC is solved separately on the P-S tables
     (compute_initial_entropy).
 
@@ -3768,11 +3770,13 @@ def zalmoxis_solver(
             else (None, None)
         )
         _interp_cache: dict = {}
+        _p_surface = float(config_params['target_surface_pressure'])
         _rho_fixed = np.zeros(len(radii))
         _T_fixed = np.zeros(len(radii))
         for _i in range(len(radii)):
             _r = float(radii[_i])
-            _P = float(pressure[_i])
+            # a padded node (P = 0 past a surface stop) is the surface, at the target pressure
+            _P = float(pressure[_i]) if pressure[_i] != 0.0 else _p_surface
             # r-indexed T: below r_ref[0] (= first Aragog staggered node
             # near CMB), clamp to T_cmb; else interp along the Aragog grid.
             if _r <= float(_r_ref[0]):

tests/interior_struct/test_zalmoxis.py

@@ -2587,6 +2587,28 @@ def tf(r, P):
     assert 0.0 < cmb_radius < hf_row['R_int']
 
 
+def test_the_jax_rebuild_looks_up_a_padded_surface_node_at_the_target_pressure(
+    tmp_path, monkeypatch
+):
+    """A padded outer node (P = 0 past a surface stop) is the surface: the post-solve rebuild
+    evaluates its density at the target surface pressure, not at 0 Pa."""
+    model_results = _plausible_model_results()
+    model_results['pressure'][-1] = 0.0
+    r_arr, t_arr = _cooled_mantle_arrays(model_results)
+    _, rho_mock, _, hf_row, _, _, _ = _run_gate_solver(
+        tmp_path,
+        monkeypatch,
+        'PALEOS:MgSiO3',
+        (r_arr, t_arr),
+        lambda r, P: float(np.interp(r, r_arr, t_arr)),
+        main_results=model_results,
+    )
+    pressures = [c.args[0] for c in rho_mock.call_args_list]
+    assert len(pressures) == len(model_results['radii'])
+    assert pressures[-1] == pytest.approx(hf_row['P_surf'] * 1e5, rel=1e-12)
+    assert pressures[:-1] == pytest.approx(list(model_results['pressure'][:-1]), rel=1e-12)
+
+
 @pytest.mark.parametrize('mantle_eos', ['PALEOS-2phase:MgSiO3', 'PALEOS:MgSiO3'])
 def test_zalmoxis_solver_passes_callable_on_jax_viable_eos_without_jax(
     tmp_path, monkeypatch, mantle_eos