Skip to content

Snapshot gravity column is not consistent with the heights written beside it #74

Description

@timlichtenberg

Problem description

Atmosphere.write_ncdf recomputes the height arrays immediately before writing, but leaves grav_z as the adiabat solve left it. The gravity and z columns in the written file therefore come from two different hydrostatic integrations of the same column.

The sequence:

  1. GeneralAdiabat fills atm.z and atm.grav_z together, deriving grav_z[idx+1] from z[idx+1] (src/janus/utils/GeneralAdiabat.py:189). At this point the two agree by construction.
  2. The resample onto the radiative grid trims and interpolates both with the same operations (src/janus/utils/GeneralAdiabat.py:874 and :883), so they still agree.
  3. write_ncdf calls integrate_heights and reassigns self.z, self.zl from it (src/janus/utils/atmosphere_column.py:394). integrate_heights does compute a local gravity inside its loop and uses it for the hydrostatic step (src/janus/utils/height.py:29), but that is a loop variable and is never written back to atm.grav_z.
  4. var_grav[:] = self.grav_z[:] (src/janus/utils/atmosphere_column.py:534) then writes the stale array next to the fresh heights.

Expected behavior

The gravity column should be the local gravity at the heights written in z beside it, that is grav_s * R**2 / (R + z)**2 evaluated on the z the file actually contains.

Evidence

This is traced through the source rather than measured from a run, so I can point at the mechanism but not yet at a magnitude. The two integrations use different formulations, since integrate_heights builds dz from R_gas * T / (mu * g * p) while the adiabat carries its own density, so how far apart they land will depend on the case. Happy to attach a measured profile from a coupled run if that is useful before anyone picks this up.

Why it matters

PROTEUS consumes this column. proteus/atmos_clim/common.py requires gravity to be present in the file and reads it, and proteus/observe/petitRADTRANS.py takes the value at the layer nearest the configured reference pressure as the reference gravity for the synthetic spectrum.

Until FormingWorlds/PROTEUS#800 no JANUS run produced a final snapshot at all, so this column only ever appeared in the per-iteration files. It is reachable now, which is why it is worth fixing rather than leaving.

My computer

macOS 15 (Darwin 25.3.0), Python 3.12, JANUS at 26.07.13.

Additional notes

integrate_heights already computes exactly the quantity that is wanted and discards it. Writing it into atm.grav_z as the loop advances, alongside z, would make the two consistent by construction instead of leaving them to agree by coincidence.

Relevant people

@nichollsh

Metadata

Metadata

Assignees

No one assigned

    Labels

    InterraPriority 4: tbdPriority level 4: nice to have features and/or has some timebugSomething isn't working

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions