proteus.outgas.wrapper¶

def calc_target_elemental_inventories(dirs: dict, config: Config, hf_row: dict):
    """
    Calculate total amount of volatile elements in the planet.

    This function pre-populates ``hf_row['O_kg_total']`` from the user's O_mode/
    O_budget settings (unless O_mode == 'ic_chemistry', in which case the
    first outgas call sets it). The atmosphere+dissolved O thus enters
    the M_ele aggregation and the Zalmoxis dry-mass subtraction, closing
    the bookkeeping asymmetry that can let M_atm exceed M_planet.
    """

    # zero by default, in case not included
    for e in element_list:
        hf_row[e + '_kg_total'] = 0.0

    # Calculate target elemental inventories from delivery config.
    # CALLIOPE's calc_target_masses computes element budgets (H, C, N, S)
    # from the delivery parameters. This is needed by ALL outgas modules,
    # not just CALLIOPE, because the element budgets drive the mass balance.
    from proteus.outgas.calliope import calc_target_masses

    calc_target_masses(dirs, config, hf_row)

    # Pre-populate O_kg_total from user budget (whole-planet O accounting,
    # issue #677). Stashes the user-supplied value as O_kg_user so the IC
    # consistency check can compare it against CALLIOPE's first-call
    # equilibrium output and flag a >50% divergence (see check_ic_oxygen_budget).
    o_kg_user = _resolve_oxygen_budget(config, hf_row)
    if o_kg_user is not None:
        hf_row['O_kg_total'] = o_kg_user
        # Stash the user budget as the one-shot IC-check baseline, but do not
        # re-arm the check once it has fired. check_ic_oxygen_budget resets
        # this to the -1.0 sentinel after it runs; re-stashing on every init
        # iteration would otherwise make the check fire repeatedly.
        if float(hf_row.get('O_kg_user_ic', 0.0)) >= 0.0:
            hf_row['O_kg_user_ic'] = o_kg_user
    else:
        # ic_chemistry mode: leave O_kg_total at 0; first CALLIOPE call writes it.
        hf_row['O_kg_user_ic'] = -1.0  # sentinel: no user budget supplied

    # Update total mass of tracked volatiles+noble gases. This includes oxygen and
    # therefore reflects the atmospheric+dissolved O mass that CALLIOPE
    # produces from the fO2 buffer. Does not include rock vapours.
    hf_row['M_ele'] = 0.0
    for e in vol_element_list + noble_gases:
        hf_row['M_ele'] += float(hf_row.get(e + '_kg_total', 0.0))

def check_desiccation(config: Config, hf_row: dict) -> bool:
    """
    Check if the planet has desiccated. This is done by checking if all volatile masses
    are below a threshold AND verifying that the loss is consistent with cumulative
    escape.

    Parameters
    ----------
    config : Config
        Configuration object
    hf_row : dict
        Dictionary of helpfile variables, at this iteration only

    Returns
    -------
    bool
        True if desiccation occurred, False otherwise

    Notes
    -----
    The escape-balance gate guards against the failure mode where an
    upstream AGNI or outgas failure wipes the atmosphere as a side
    effect (e.g. NaN volume mixing ratios that get zeroed downstream),
    causing all `*_kg_total` to collapse below `mass_thresh` in a single
    iteration without `esc_kg_cumulative` having moved. The CHILI sweep
    R7/R21 cascades fired this exact pattern: AGNI failed, the
    atmosphere was wiped, and the old check_desiccation reported success
    instead of letting the deadlock detector catch the upstream failure.

    The gate compares (M_vol_initial - current M_ele) against
    `1.5 * esc_kg_cumulative + 1e3 kg`. The 1.5x slack absorbs rounding
    and threshold-truncation noise; the 1e3 kg floor prevents the gate
    from firing on pathologically tiny inventories. If the gate is
    inactive (no baseline tracked yet, e.g. resuming an old CSV without
    `M_vol_initial`), the function falls back to the old threshold-only
    behaviour. The current total skips a non-finite element, as run_desiccated
    does, and mass an earlier desiccation removed, ``M_desiccated``, is left out
    of the loss the gate compares, so a desiccated row passes again on resume.

    Scope of "desiccated": the threshold loop below considers the volatile
    elements (`vol_element_list`: H, O, C, N, S) together with the noble gases,
    so a planet with a noble-gas dominated atmosphere is not desiccated.
    """

    # Threshold check: refuse desiccation while ANY volatile or noble gas is
    # still above the per-element mass threshold. Issue #677 fix: O is included
    # now. Under D1A the user's choice was to require O_kg_total below threshold
    # as well, on the grounds that an atmosphere with substantial O is not
    # meaningfully "desiccated" even if H/C/N/S are depleted. In practice
    # CALLIOPE drives O_kg_total to near-zero once H/C/N/S vanish, so this
    # change rarely affects the desiccation timing, but it keeps the
    # semantics honest under whole-planet O accounting.
    unreadable = [
        e for e in element_list if not np.isfinite(float(hf_row.get(f'{e}_kg_total', 0.0)))
    ]
    if unreadable:
        log.error(
            'Desiccation check refused: the total of %s is not finite; an upstream '
            'step left it unreadable and needs checking.',
            ', '.join(unreadable),
        )
        return False
    for e in vol_element_list + noble_gases:
        # Mass trapped in the solid, or held in condensates kept from re-vapourising,
        # can neither escape nor outgas, so only the escapable inventory counts;
        # otherwise it would block desiccation forever.
        reachable = escapable_inventory(hf_row, e, config)
        if reachable > config.outgas.mass_thresh:
            log.info('Not desiccated, %s = %.2e kg reachable' % (e, reachable))
            return False  # return, and allow run_outgassing to proceed

    # Escape-balance gate. Only enforced when a baseline has been
    # snapshotted by `escape.wrapper.run_escape` (first escape call).
    m_init_raw = hf_row.get('M_vol_initial', None)
    try:
        m_init = float(m_init_raw) if m_init_raw is not None else 0.0
    except (TypeError, ValueError):
        m_init = 0.0
    if not np.isfinite(m_init) or m_init <= 0.0:
        # No baseline -> fall back to old threshold-only behaviour.
        return True

    # cur_m_ele must sum over the same element set as the M_vol_initial
    # baseline (escape.wrapper.run_escape sums over the full element_list),
    # or the (m_init - cur_m_ele) comparison mixes element sets and
    # overstates "lost" mass by any noble-gas/rock-vapour inventory present
    # at baseline time.
    removed = float(hf_row.get('M_desiccated', 0.0))
    if not np.isfinite(removed):
        removed = 0.0
    cur_m_ele = readable_total(hf_row)
    lost = m_init - cur_m_ele - removed
    esc_cum = float(hf_row.get('esc_kg_cumulative', 0.0))

    # Allow 1.5x scaling slack plus a 1 t absolute floor for noise.
    if lost > 1.5 * esc_cum + 1.0e3:
        log.error(
            'Desiccation check refused: %.2e kg of volatile mass loss is '
            'unexplained by cumulative escape (%.2e kg). Initial=%.2e kg, '
            'current=%.2e kg. This usually indicates an AGNI or outgas-side '
            'failure that wiped the atmosphere as a side effect, not real '
            'escape. The deadlock detector should fire on the next iteration.',
            lost - esc_cum,
            esc_cum,
            m_init,
            cur_m_ele,
        )
        return False

    return True

def check_ic_oxygen_budget(
    config: Config,
    hf_row: dict,
    threshold_frac: float = 0.5,
) -> None:
    """Verify the user-supplied IC oxygen budget against CALLIOPE's chemistry.

    After the first ``run_outgassing`` call at IC, CALLIOPE (or
    atmodeller) has written ``hf_row['O_kg_total']`` from the fO2-
    buffered equilibrium. We compare that against the user-supplied
    O_budget that ``calc_target_elemental_inventories`` stashed earlier
    in ``hf_row['O_kg_user_ic']``. A large divergence usually means
    either:

    1. The user under-specified O_budget so the planet cannot actually
       support the equilibrium atmosphere CALLIOPE wants. The implied
       mantle FeO reservoir is being over-drawn.
    2. The user over-specified O_budget so M_planet now carries phantom
       oxygen mass that has no physical home (not in the atmosphere,
       not in the mantle FeO that PALEOS density accounts for).

    Both cases are best caught loudly at IC rather than silently
    corrupting the trajectory.

    Skipped when:
      - ``planet.fO2_source != "user_constant"``. When
        fO2_source = 'from_O_budget' the user O budget is the
        *authoritative* input that drives the chemistry,
        so any "divergence" between it and the derived O inventory is
        zero by construction. The check is meaningful only when fO2 is
        the input and O is the output.
      - ``O_kg_user_ic`` is the sentinel (-1.0), which marks
        ``O_mode == 'ic_chemistry'`` (the user opted into chemistry-
        derived O so any discrepancy is, by definition, accepted).

    Parameters
    ----------
    config : Config
        Configuration object. Reads ``config.planet.fO2_source`` to
        decide whether the check applies.
    hf_row : dict
        Helpfile row after the first outgas call at IC.
    threshold_frac : float
        Maximum allowed relative divergence between user O_budget and
        CALLIOPE's IC equilibrium value. Default 0.5 (50 percent).

    Raises
    ------
    ValueError
        If divergence exceeds ``threshold_frac``.
    """
    # Gate on fO2_source. The 'from_O_budget' source makes user O
    # authoritative so there is no divergence to flag.
    if config.planet.fO2_source != 'user_constant':
        return

    user_O = float(hf_row.get('O_kg_user_ic', -1.0))
    if user_O < 0.0:
        # Sentinel: O_mode == 'ic_chemistry' or check already fired.
        return

    chem_O = float(hf_row.get('O_kg_total', 0.0))
    if chem_O <= 0.0:
        # Degenerate state; nothing to compare against. Skip.
        return

    rel_div = abs(user_O - chem_O) / chem_O
    log.info(
        'IC oxygen budget check: user=%.3e kg, CALLIOPE-equilibrium=%.3e kg, '
        'relative divergence=%.1f%% (threshold %.0f%%)',
        user_O,
        chem_O,
        rel_div * 100,
        threshold_frac * 100,
    )

    if rel_div > threshold_frac:
        raise ValueError(
            f'IC oxygen budget mismatch (issue #677 consistency check): '
            f'user O_budget implies {user_O:.3e} kg, CALLIOPE equilibrium at IC '
            f'gives {chem_O:.3e} kg (relative divergence {rel_div * 100:.1f}%, '
            f'threshold {threshold_frac * 100:.0f}%). Either: '
            f'(a) adjust planet.elements.O_budget closer to the chemistry value; '
            f'(b) switch to O_mode="ic_chemistry" to defer to CALLIOPE; '
            f'(c) change outgas.fO2_shift_IW to bring chemistry into line with '
            f'the budget. Persistent disagreement usually indicates the implied '
            f'mantle FeO reservoir is being over- or under-drawn.'
        )

    # Mark check as done so subsequent outgas calls in init_stage don't re-fire.
    hf_row['O_kg_user_ic'] = -1.0

@contextmanager
def condensed_mass_withheld(config: Config, hf_row: dict, *, initial: bool):
    """Start the chemistry from no condensed mass, keeping held condensates out of it.

    Every backend starts from empty ``{e}_kg_condensed`` and ``{X}_cd_kg``
    columns; only atmodeller with condensates enabled fills them, so they stay 0
    otherwise. With ``outgas.prevent_desaturation`` on, the condensed mass so far
    is also taken out of each element total for the solve, so the chemistry cannot
    re-vapourise it, and on exit it is added back to the totals, to
    ``{e}_kg_condensed`` and to each condensate's ``{X}_cd_kg``, on top of what the
    solve newly condensed. The reservoirs then only grow. The exit runs also when
    the solve raises, so the row then holds the held mass, no more and no less.
    ``{X}_cd_activity`` stays the solver's value: it is the saturation of the
    inventory left in the solve, not of the held layer, so it can be below 1 while
    a layer is held.

    In the init stage (``initial``) nothing is held back: each init iteration
    resets the element totals to the delivered budget and solves the same starting
    state again, so condensate from an earlier trial solve is not kept.
    """
    held = {e: 0.0 for e in vol_element_list}
    if not initial:
        held = {e: locked_condensed_mass(hf_row, e, config) for e in vol_element_list}
    held_cd = {}
    if any(held.values()):
        for c in condensate_list:
            mass = float(hf_row.get(f'{c}_cd_kg', 0.0))
            held_cd[c] = max(0.0, mass) if np.isfinite(mass) else 0.0
    for c in condensate_list:
        hf_row[f'{c}_cd_kg'] = 0.0
    for e in vol_element_list:
        hf_row[f'{e}_kg_condensed'] = 0.0
        if held[e] > 0.0:
            hf_row[f'{e}_kg_total'] = float(hf_row.get(f'{e}_kg_total', 0.0)) - held[e]
    try:
        yield
    finally:
        for e in vol_element_list:
            if held[e] > 0.0:
                hf_row[f'{e}_kg_total'] = float(hf_row.get(f'{e}_kg_total', 0.0)) + held[e]
                hf_row[f'{e}_kg_condensed'] = (
                    float(hf_row.get(f'{e}_kg_condensed', 0.0)) + held[e]
                )
        for c, mass in held_cd.items():
            hf_row[f'{c}_cd_kg'] = float(hf_row.get(f'{c}_cd_kg', 0.0)) + mass

def desiccated_after_trapping(config: Config, hf_row: dict, desiccated: bool, step) -> bool:
    """Whether a desiccated planet is still desiccated after this step's trapping.

    A remelting step returns buried mass to the melt, where :func:`run_desiccated`
    keeps it while trapping is on. Once the mass the melt holds brings an element
    above ``mass_thresh`` the planet holds volatiles again, and the flag clears so
    that escape and outgassing take that mass from there.

    Parameters
    ----------
    config : Config
        Configuration object
    hf_row : dict
        Dictionary of helpfile variables, at this iteration only
    desiccated : bool
        Whether the planet was desiccated before this step
    step : TrappingStep or None
        What this step's trapping did, as returned by ``run_trapping``

    Returns
    -------
    bool
        Whether the planet is desiccated now.
    """
    if not desiccated or step is None or not step.remelted or step.total_trapped >= 0.0:
        return desiccated
    still = check_desiccation(config, hf_row)
    if not still:
        log.info(
            'Remelting returned %.3e kg to the melt; no longer desiccated', -step.total_trapped
        )
    return still

def outgassing_derives_o_kg_total(config: Config) -> bool:
    """Whether each outgassing call derives ``O_kg_total`` from the melt fO2.

    CALLIOPE and atmodeller derive the O budget from the fixed fO2 under
    ``planet.fO2_source = 'user_constant'``; ``O_mode = 'ic_chemistry'``
    requires that source at config load. The dummy outgassing derives an empty
    ``O_kg_total`` only in the init stage and otherwise keeps it as given, and
    ``from_O_budget`` takes the O budget as input.
    """
    return (
        config.outgas.module in ('calliope', 'atmodeller')
        and config.planet.fO2_source == 'user_constant'
    )

def run_crystallized(config: Config, hf_row: dict, dt: float):
    """Handle a crystallized mantle: volatile exchange frozen, but escape
    still erodes the atmosphere.

    After the mantle solidifies (Phi_global < phi_crit), volatiles can no
    longer exchange between mantle and atmosphere: dissolved volatiles are
    trapped in the solid and outgassing no longer replenishes the atmosphere.

    Escape, however, continues. ``run_escape`` runs earlier in the main loop
    with ``atmosphere_only=True`` in this regime, so it debits each element
    total (``*_kg_total``) by its share of the escaping mass, its atmospheric
    mass over the sum of the element masses in the atmosphere. The atmospheric
    reservoirs are scaled here by the fraction of that same sum that remains,
    so each element's atmosphere loses exactly what escape took from its total
    and the per-element totals and reservoirs stay consistent. ``M_atm``, the
    sum of the species masses, can differ from the element sum where the
    chemistry's molar masses do (by 1.8e-5 with CALLIOPE); scaling it by
    ``(M_atm - esc) / M_atm`` instead would leave that fraction of every step's
    loss in the atmosphere. The scaling is uniform, so partial-pressure ratios,
    VMRs and the mean molecular weight are unchanged.

    Parameters
    ----------
    config : Config
        Configuration object. Used to confirm the escape model is unfractionated.
    hf_row : dict
        Dictionary of helpfile variables, at this iteration only.
    dt : float
        Length of the current step [yr], used to size the escaped mass.
    """
    # Uniform scaling is valid only while escape removes atmospheric species in
    # proportion to their abundance. A fractionating model would change the
    # composition and require per-species debiting, so refuse it explicitly at
    # the point that relies on the assumption.

    reservoir = getattr(config.escape, 'reservoir', 'outgas')
    if reservoir not in ('outgas', 'bulk'):
        raise NotImplementedError(
            f'run_crystallized assumes unfractionated escape; escape.reservoir='
            f"'{reservoir}' is not supported once the mantle has solidified."
        )

    m_atm = float(hf_row.get('M_atm', 0.0))
    esc_rate = float(hf_row.get('esc_rate_total', 0.0))
    # The element masses in the atmosphere, the reservoir escape shared its loss over;
    # a row that records none has only the species sum.
    m_elements = sum(float(hf_row.get(f'{e}_kg_atm', 0.0)) for e in element_list) or m_atm

    if dt <= 0.0:
        # A non-positive step is a coupling error, not a benign no-op: surface it
        # so a regression that stalls the clock here does not look like "escape
        # silently stopped removing atmosphere".
        log.warning(
            'Crystallized mantle: non-positive dt=%.3e yr; skipping atmospheric '
            'escape debit this step.',
            dt,
        )
        return

    if m_atm <= 0.0 or esc_rate <= 0.0 or not m_elements > 0.0:
        # No atmosphere, no active escape, or an unreadable element atmosphere:
        # reservoirs stay as-is.
        log.info('Crystallized mantle: volatile exchange frozen, reservoirs preserved')
        return

    # Mass lost to escape over this step. With the mantle frozen it comes
    # entirely out of the atmosphere. Prefer the loss the escape step actually
    # applied: sizing this from the raw rate instead would empty the atmosphere
    # on an overshoot while the elemental totals kept all but the capped share.
    esc_step_kg = hf_row.get('esc_step_kg')
    if esc_step_kg is None or not np.isfinite(float(esc_step_kg)):
        esc_step_kg = esc_rate * secs_per_year * dt
    esc_step_kg = float(esc_step_kg)

    if esc_step_kg > m_elements:
        # Escape took more from the elemental totals than the column holds, so
        # the loss was sized from a reservoir the frozen mantle no longer
        # supplies and the two records of this step disagree by the excess.
        log.warning(
            'Escape removed %.3e kg from the elemental totals while the '
            'atmosphere holds only %.3e kg; the column empties and the excess '
            '%.3e kg leaves no reservoir that tracks it.',
            esc_step_kg,
            m_elements,
            esc_step_kg - m_elements,
        )
    retained = max(0.0, (m_elements - esc_step_kg) / m_elements)

    # Scale the atmospheric reservoirs by the retained fraction. Uniform
    # scaling preserves composition, so `*_vmr` and `atm_kg_per_mol` (mmw)
    # are left unchanged.
    # Scale every gas reservoir by the retained fraction when the frozen
    # atmosphere escapes. The noble gases are members of gas_list, so their
    # atmospheric mass shrinks with the reactive volatiles; otherwise a later
    # escape step that debits the noble total would leave the noble
    # atmospheric mass stale and larger than the total.
    for s in gas_list:
        hf_row[f'{s}_kg_atm'] = hf_row.get(f'{s}_kg_atm', 0.0) * retained
        hf_row[f'{s}_bar'] = hf_row.get(f'{s}_bar', 0.0) * retained
    # Element reservoirs whose symbol is not also a species name. The loop above
    # already scales the rock formers and noble gases through the species they
    # share a symbol with; H, O, C, N and S have no such twin, and escape sizes
    # the next step's loss from these keys, so leaving them would hold the
    # escapable reservoir at its pre-solidification value as the column empties.
    for e in element_list:
        if e in gas_list:
            continue
        hf_row[f'{e}_kg_atm'] = hf_row.get(f'{e}_kg_atm', 0.0) * retained
    hf_row['P_surf'] = hf_row.get('P_surf', 0.0) * retained
    # Keep the volatile/vapour partial-pressure split consistent with the
    # rescaled total (same uniform, composition-preserving scaling as above).
    hf_row['P_vol'] = hf_row.get('P_vol', 0.0) * retained
    hf_row['P_vap'] = hf_row.get('P_vap', 0.0) * retained
    hf_row['M_atm'] = m_atm * retained

    # The rock-vapour mass follows the same uniform scaling as its partial
    # pressure. Left unscaled it would freeze at its last pre-solidification
    # value while the vapour it describes is eroded by escape.
    hf_row['M_vaps'] = hf_row.get('M_vaps', 0.0) * retained
    # Track the atmospheric volatile mass separately from M_atm, which also
    # carries the rock vapour when vapourisation is enabled. Summing
    # vol_gas_list rather than scaling M_atm keeps the column equal to the
    # species masses it is defined as, which assert_mass_conservation checks.
    hf_row['M_vol_atm'] = sum(float(hf_row.get(s + '_kg_atm', 0.0)) for s in vol_gas_list)

    log.info(
        'Crystallized mantle: volatile exchange frozen; escape removed '
        '%.3e kg from the atmosphere (retained fraction %.4f)',
        m_atm - hf_row['M_atm'],
        retained,
    )

def run_desiccated(dirs: dict, config: Config, hf_row: dict, first_iter: bool):
    """
    Handle desiccation of the planet. This substitutes for run_outgassing when the planet
    has lost its entire volatile inventory. The atmosphere and the melt are emptied,
    noble gases included; the solid mantle keeps its share, which becomes each
    species and element total. With solid-phase trapping on, the melt is kept as
    well: a remelt returns buried mass to it, and emptying it every step would make
    what the planet keeps depend on how many steps the remelt takes. The mass
    removed is booked in ``M_desiccated`` and the Zalmoxis target.

    The kept share is the sum of the kept columns limited to [0, total before the
    call], each column scaled alike; a non-finite column keeps nothing and is
    logged. The kept columns, kg and mol, are set to their share. No total rises,
    so the mass removed is not negative and a row that passed
    ``check_desiccation`` passes it again. A non-finite total stays as it is, with
    its mol total set to NaN. A non-finite ``M_desiccated`` reads as 0, as in the
    check.

    Parameters
    ----------
    dirs : dict
        Dictionary of directory paths
    config : Config
        Configuration object
    hf_row : dict
        Dictionary of helpfile variables, at this iteration only
    first_iter : bool
        Whether the run is in its init stage.
    """

    # if desiccated, set all gas masses to zero
    log.info('Desiccation has occurred - no volatiles remaining')

    # With trapping on the melt is kept as well, so mass a remelt returns stays in it.
    keep_melt = trapping_active(config)
    held = ('_solid', '_liquid') if keep_melt else ('_solid',)
    totals = {
        n: float(hf_row.get(f'{n}_kg_total', 0.0))
        for n in dict.fromkeys(element_list + gas_list)
    }
    removed = readable_total(hf_row)
    for k in expected_keys():
        if k not in EMPTY_ATMOSPHERE_KEPT_KEYS and not k.endswith(held):
            hf_row[k] = 0.0
    for n, total in totals.items():
        if not np.isfinite(total):
            hf_row[f'{n}_kg_total'] = total
            if n in gas_list:
                hf_row[f'{n}_mol_total'] = np.nan
            continue
        solid = _kept_column(float(hf_row.get(f'{n}_kg_solid', 0.0)), n, 'solid')
        liquid = 0.0
        if keep_melt:
            liquid = _kept_column(float(hf_row.get(f'{n}_kg_liquid', 0.0)), n, 'melt')
        # Condensates (e.g. graphite) stay on the surface when the atmosphere goes.
        condensed = 0.0
        if n in vol_element_list:
            condensed = _kept_column(
                float(hf_row.get(f'{n}_kg_condensed', 0.0)), n, 'condensed'
            )
        mantle = solid + liquid
        held_mass = mantle + condensed
        kept = min(held_mass, total)
        share = kept / held_mass if held_mass > 0.0 else 0.0
        # A negative total keeps its value, in the solid column.
        hf_row[f'{n}_kg_solid'] = solid * share if held_mass > 0.0 else kept
        if keep_melt:
            hf_row[f'{n}_kg_liquid'] = liquid * share
        if n in vol_element_list:
            hf_row[f'{n}_kg_condensed'] = condensed * share
        hf_row[f'{n}_kg_total'] = kept
        if n in gas_list:
            mol = _kept_column(float(hf_row.get(f'{n}_mol_solid', 0.0)) * share, n, 'solid')
            hf_row[f'{n}_mol_solid'] = mol
            if keep_melt:
                mol_melt = float(hf_row.get(f'{n}_mol_liquid', 0.0)) * share
                hf_row[f'{n}_mol_liquid'] = _kept_column(mol_melt, n, 'melt')
                mol += hf_row[f'{n}_mol_liquid']
            hf_row[f'{n}_mol_total'] = mol
    removed -= readable_total(hf_row)
    booked = float(hf_row.get('M_desiccated', 0.0))
    hf_row['M_desiccated'] = (booked if np.isfinite(booked) else 0.0) + removed
    record_volatile_change(config, hf_row, -removed)

    # Vapourisation of refractories, under the same crystallised gate as
    # volatile outgassing path.
    if config.outgas.vapourise:
        if hf_row['Phi_global'] <= config.params.stop.solid.phi_crit:
            log.info('Skipping rock vapourisation for crystallised mantle')
        else:
            log.info('Calculating rock vapourisation at surface')
            run_vapourisation(dirs, config, hf_row, first_iter)

def run_outgassing(dirs: dict, config: Config, hf_row: dict, *, initial: bool):
    """
    Run outgassing model to get new volatile surface pressures

    Parameters
    ----------
    dirs : dict
        Dictionary of directory paths
    config : Config
        Configuration object
    hf_row : dict
        Dictionary of helpfile variables, at this iteration only
    initial : bool
        Whether this is an init-stage iteration, where the dummy outgassing
        derives an empty O budget from the outgassed species.
    """

    log.info('Calculating volatile outgassing at surface')

    # planet.fO2_source dispatch. Two runtime paths are wired today:
    # 'user_constant' (fO2 buffered to the configured IW offset) for
    # every backend, and 'from_O_budget' (authoritative-O chemistry)
    # for the CALLIOPE
    # and atmodeller backends. The Config-level validator
    # (planet_fO2_source_compat) rejects 'from_mantle_redox' and the
    # from_O_budget + dummy combo at config-load, so this guard is
    # unreachable under a normally-loaded Config. It remains as defence
    # in depth for programmatic Config construction (tests, scripted
    # runs) that bypasses the validators.
    fO2_source = config.planet.fO2_source
    if fO2_source not in ('user_constant', 'from_O_budget'):
        raise NotImplementedError(
            f'planet.fO2_source = "{fO2_source}" is recognised by the '
            'config schema but its runtime path is not yet wired into '
            'run_outgassing.'
        )

    if fO2_source == 'from_O_budget' and config.outgas.module == 'dummy':
        raise NotImplementedError(
            'planet.fO2_source = "from_O_budget" requires chemistry to '
            'invert against; outgas.module = "dummy" has none.'
        )

    # Default the derived-fO2 helpfile columns to the user-configured
    # buffer offset before the backend dispatch. Each backend overrides
    # the default as appropriate for its chemistry. Under fO2_source =
    # "user_constant" the pre-seeded value is the offset the chemistry
    # equilibrated to (because the fugacity constraint set it); leaving
    # it untouched is the right behaviour. Under fO2_source =
    # "from_O_budget" CALLIOPE and atmodeller overwrite the pre-seed
    # with their solver-derived value. The pre-seed exists primarily so
    # the dummy backend and any solve-skipped path land at the
    # physically meaningful configured value rather than the
    # ZeroHelpfileRow default 0.0, which downstream consumers would
    # otherwise mis-interpret as "the chemistry ran at exactly IW".
    hf_row['fO2_shift_IW_derived'] = float(config.outgas.fO2_shift_IW)
    hf_row['O_res'] = 0.0

    # The chemistry partitions a whole-planet inventory and rewrites every
    # `_kg_solid` column (0.0 in CALLIOPE and atmodeller; atmodeller writes its
    # condensates to `_kg_condensed`), so the trapped mass is hidden for the solve
    # and put back after it. Condensed mass is handled the same way when it is kept
    # from re-vapourising.
    with (
        condensed_mass_withheld(config, hf_row, initial=initial),
        trapped_mass_withheld(hf_row),
    ):
        _solve_chemistry(dirs, config, hf_row, initial=initial)

    # P_surf here is the volatile+noble gas total
    hf_row['P_vol'] = hf_row['P_surf']
    hf_row['P_vap'] = 0.0

    # calculate total atmosphere mass
    hf_row['M_atm'] = 0.0
    for s in gas_list:
        hf_row['M_atm'] += float(hf_row.get(s + '_kg_atm', 0.0))
    # Volatile-only atmospheric mass (rock vapour excluded). The caller
    # zeroes the rock-vapour species before this function runs, so the two
    # agree here; setting the column explicitly keeps the
    # assert_mass_conservation bookkeeping check live from the first
    # iteration rather than only after the mantle crystallises.
    hf_row['M_vol_atm'] = sum(float(hf_row.get(s + '_kg_atm', 0.0)) for s in vol_gas_list)

    # Derive element mass ratios in atmosphere
    for e1 in element_list:
        for e2 in element_list:
            key = f'{e2}/{e1}_atm'  # key to be set in helpfile
            if key not in hf_row:
                continue
            if min(hf_row[f'{e1}_kg_atm'], hf_row[f'{e2}_kg_atm']) < 1e-30:
                hf_row[key] = 0.0
            else:
                hf_row[key] = hf_row[f'{e2}_kg_atm'] / hf_row[f'{e1}_kg_atm']

    # print outgassed partial pressures (in order of descending abundance)
    species = [s for s in gas_list if s not in vap_list]
    vmrs = np.array([hf_row[s + '_vmr'] for s in species])

    for i in np.argsort(vmrs)[::-1]:
        s = species[i]
        _p = hf_row[s + '_bar']
        _x = hf_row[s + '_vmr']
        _s = '    %-6s     = %-9.2f bar (%.2e VMR)' % (s, _p, _x)
        if _p > 0.01:
            log.info(_s)
        else:
            # don't spam log with species of negligible abundance
            log.debug(_s)

    # print total pressure and mmw
    log.info('    total      = %-9.2f bar' % hf_row['P_surf'])
    log.info('    mmw        = %-9.5f g mol-1' % (hf_row['atm_kg_per_mol'] * 1e3))

def run_outgassing_and_vapourisation(
    dirs: dict, config: Config, hf_row: dict, first_iter: bool
):
    """Runs volatile outgassing and rock vapourisation together and combines the results.

    Two step process:
     - outgassing calculated from solubility+thermochemical equilibrium
     - vapourisation of refractory species from the surface.

    Vapourisation re-equilibrates the combined gas, allowing for
    computation of melt outgassing and dissolution.

    Parameters
    ----------
    dirs : dict
        Dictionary of directory paths
    config : Config
        Configuration object
    hf_row : dict
        Dictionary of helpfile variables, at this iteration only
    first_iter : bool
        Whether the run is in its init stage; the dummy outgassing derives an
        empty O budget only then.
    """

    # reset all rock-vapour masses to zero:
    hf_row['M_vaps'] = 0.0
    for s in gas_list:
        if s not in vap_list:
            continue
        else:
            hf_row[s + '_bar'] = 0.0
            hf_row[s + '_vmr'] = 0.0
            hf_row[s + '_kg_atm'] = 0.0
            hf_row[s + '_kg_total'] = 0.0

    # Only the rock-forming elements are reset: the volatile and noble totals are
    # owned by the outgassing and escape steps respectively.
    for e in element_list:
        if e in vol_element_list or e in noble_gases:
            continue
        hf_row[e + '_kg_atm'] = 0.0
        hf_row[e + '_kg_total'] = 0.0

    # Volatile outgassing
    run_outgassing(dirs, config, hf_row, initial=first_iter)

    # Vapourisation of refractories
    if config.outgas.vapourise:
        if hf_row['Phi_global'] <= config.params.stop.solid.phi_crit:
            log.info('Skipping rock vapourisation for crystallised mantle')
        else:
            log.info('Calculating rock vapourisation at surface')
            run_vapourisation(dirs, config, hf_row, first_iter)