Skip to content

planetWRF support: modern MarsWRF/planetWRF output through the full CAP chain, plus a time-shift fix - #198

Merged
rurata merged 12 commits into
develfrom
aeolis
Aug 28, 2026
Merged

planetWRF support: modern MarsWRF/planetWRF output through the full CAP chain, plus a time-shift fix#198
rurata merged 12 commits into
develfrom
aeolis

Conversation

@mirfjc

@mirfjc mirfjc commented Aug 27, 2026

Copy link
Copy Markdown

Summary

Makes MarsFormat -gcm marswrf work on current planetWRF output, both full wrfout files and reduced, published subsets, so that planetWRF runs go through the whole CAP chain (MarsFormat, MarsInterp, MarsFiles, MarsVars, MarsPlot) unchanged. Along the way it fixes three bugs, one of which affects all models (see "Changes that affect Ames users").

Tested on:

  • a full planetWRF wrfout (641 variables, 36x18x100, 8 frames per sol, 2 sols),
  • a reduced published file (16 variables, no grid metadata, 4 frames per sol, 2 sols),
  • one full Mars year of the reduced file (7.4 GB, 2676 frames): MarsFormat 28 s, MarsInterp -t pstd 53 min, MarsFiles -bd -ba 10 6 s, MarsPlot time-lat pages. Peak memory about 3x the input file size (whole-file loading; same as before this PR).

pytest tests/test_marsformat.py tests/test_marsinterp.py tests/test_marsfiles.py tests/test_marsvars.py tests/test_marsplot.py: 60 passed, 4 skipped.

MarsWRF / planetWRF path (bin/MarsFormat.py)

Every planetWRF-specific quantity is now taken from the newest source the file offers, with fallbacks for older or reduced files. The chosen source is printed.

quantity first choice fallbacks
temperature T_PHY (T+T0)*((P+PB)/P0)^(R_D/CP)
pressure (pfull3D) P_PHY P+PB
time [sols] Times strings XTIME + SIMULATION_START_DATE, JULIAN, START_DATE + spacing from L_S
Ls (areo) L_S sol2ls(time)
eta, ak/bk ZNW, ZNU, P_TOP least-squares fit of P_PHY vs PSFC
winds U_PHY, V_PHY, W_PHY (mass points) 2-point destagger of U, V, W
heights (zfull) Z_PHY - HGT Z_PHY minus hydrostatic surface estimate; PH+PHB
local time LTST (true solar, as eot_offset) UT + lon/15
gas constant, cp rgas3D, cp3D from MW_AIR_3D, RCP_3D R_D, CP attributes
gravity g3D when DO_VARIABLE_GRAVITY = 1 G attribute

Bugs fixed in the existing marswrf code:

  • Temperature was derived from P_TOP + PB[0] (double-counting P_TOP, ignoring the perturbation pressure and using time 0). Against the model's own T_PHY this was 5.7 K rms, up to 53 K near topography; now 0.003 K.
  • Time axis was left in minutes: planetWRF writes XTIME with units = "".
  • Longitude was converted to 0-360 but not sorted (185..355, 5..175).
  • Variables on non-CAP dimensions (soil layers, dust bins, radiation layers, ZNU, C1H, ...) broke MarsInterp; they are now dropped with a message.

New option -stag / --keep_staggered: also write the winds on their native Arakawa-C locations (u_stag, v_stag, w_half with lon_u, lat_v, phalf axes) plus interface pressure phalf3D and height zhalf. w_half, phalf3D, zhalf go through MarsInterp; u_stag, v_stag are for users' own scripts.

Changes that affect Ames users

  • time_shift_calc (amescap/FV3_utils.py) was one sample interval off. The default branch computed the input index as n*dt_out - lon_shift - tod[0] + dt_in; the + dt_in is one full bin (and the explicit -t "3 15" branch omitted - tod[0]). On 3-hourly data every time-shifted field moved by 3 h (surface temperature peaked at 09:00 local). On the Ames 24-per-sol, bin-centred diurn axis (0.5, 1.5, ... h) the same error is 0.5 h. Both branches now use (UT - tod[0]) / dt_in with UT = LT - lon_shift, verified against hand linear interpolation at four longitudes. Existing 3 AM / 3 PM comparisons will move by half an hour.
  • MarsInterp -t zstd now falls back to zsurf inside the input file when no DDDDD.fixed.nc exists (previously an error exit).
  • MarsInterp's copy-through test 'pfull' not in dims or 'phalf' not in dims was always true; variables on the model vertical grid but not on lat/lon crashed the run. They are now skipped with a message.
  • MarsVars: div and curl use the exact C-grid forms (cgrid_div, cgrid_curl) when u_stag/v_stag are in the file, else the existing A-grid spherical_div/spherical_curl. Planet constants (rgas, Cp, g, radius, R/cp) are read from the file only when all four WRF attributes R_D, CP, G, RADIUS are present; Ames files keep every default, including the original theta exponent. fms_Z_calc and mass_stream gained keyword arguments with the old constants as defaults.
  • MarsFiles -t applies eot_offset (true local solar time) only when the file carries it.
  • MarsPlot: optional nlev, cticks, cfmt fields appended to Axis Options (last commit, independent of the rest; drop it if unwanted). Old templates render identically.

Notes

  • The branch was cut from main after PR Devel #195, so it carries the three main commits not yet in devel (the Devel #195 merge and the pypdf bump).
  • Not addressed here, could follow as a separate PR: vinterp copies the full array several times per output level (flatten() inside the level loop), find_n is a Python double loop, and MarsInterp/MarsFormat hold whole files in memory; time-chunking would cut the 21 GB peak to well under 1 GB.
  • Docs: new page docs/source/planetwrf.rst (quick start, source-priority table, -stag, constants, true solar time).

dependabot Bot and others added 12 commits July 26, 2026 18:17
Bumps [pypdf](https://github.com/py-pdf/pypdf) from 6.13.3 to 6.14.2.
- [Release notes](https://github.com/py-pdf/pypdf/releases)
- [Changelog](https://github.com/py-pdf/pypdf/blob/main/CHANGELOG.md)
- [Commits](py-pdf/pypdf@6.13.3...6.14.2)

---
updated-dependencies:
- dependency-name: pypdf
  dependency-version: 6.14.2
  dependency-type: direct:production
...

Signed-off-by: dependabot[bot] <support@github.com>
…pdf-6.14.2

Bump pypdf from 6.13.3 to 6.14.2
The marswrf path derived temperature from (P_TOP + PB[0]), which
double-counts P_TOP (PB already includes it) and drops the perturbation
pressure P. Against planetWRF's own T_PHY this gave 5.7 K rms and up to
53 K errors near topography.

Now: use the physics-grid T_PHY and P_PHY directly when present and
filled (they are zero on the first frame after a restart, where the
reconstruction is used instead); otherwise reconstruct with the full
pressure, temp = (T + T0) * ((P + PB) / P0) ** (R_D / CP). Both paths
reproduce T_PHY to ~0.01 K. The full 3D pressure is saved as pfull3D.
Tested on a 2-sol bite of a full planetWRF wrfout (641 variables,
36x18x100) and on a reduced, published file (16 variables, no grid
metadata). Both now go through MarsFormat, MarsInterp -t pstd and
MarsPlot cleanly.

- Time axis [sols] from the best clock in the file, in priority:
  'Times' date strings > XTIME + SIMULATION_START_DATE > JULIAN +
  MODEL_MARS_YEAR > START_DATE + uniform spacing inferred from L_S.
  planetWRF writes XTIME with units='' so the old unit sniffing left
  the time axis in minutes. Ls is computed with sol2ls() only if L_S
  is absent.
- Longitude converted to 0-360 AND sorted; the bare modulo left a
  non-monotonic coordinate (185..355, 5..175).
- Vertical coordinate: ak = P_TOP*(1-eta), bk = eta from ZNW/ZNU;
  when ZNU/ZNW/P_TOP are absent (reduced files) eta and P_TOP are
  recovered by a per-level least-squares fit of P_PHY against PSFC.
- zfull (layer height above the local surface) saved when PH/PHB/HGT
  are present, instead of being computed and discarded.
- U_PHY/V_PHY/W_PHY aliased to the wind names when the staggered
  U/V/W are absent, so ucomp/vcomp/w exist after renaming.
- Variables on non-CAP dimensions (soil layers, dust bins, radiation
  layers) and vertical-only metadata (ZNU, C1H, ...) are dropped with
  a message; MarsInterp could not copy them into interpolated files.
- 'Times' drop guarded for files that do not carry it.

tests/test_marsformat.py: 9 passed.
…MarsFormat

time_shift_calc (amescap/FV3_utils.py): the default branch computed the
input index as n*dt_out - lon_shift - tod[0] + dt_in, i.e. one full
sample interval too late, and the target_times branch omitted -tod[0].
On a planetWRF diurn file with 8 samples per sol this shifted every
field by 3 h: surface temperature peaked at 09:00 local after the
shift instead of 12:00. With bin-centred Ames diurn axes (0.5, 1.5,
... h) the same term is a 0.5 h error.

Both branches now use (UT - tod[0]) / dt_in with UT = LT - lon_shift.
Verified against hand linear interpolation of the UT series at four
longitudes for both the default and the explicit -t "12 15" forms
(agreement to 0.1 K); the diurnal peak is now at 12:00 local at every
longitude.

MarsFormat (marswrf): the non-CAP-dimension prune now also covers
coordinate variables. xarray promotes XLAT_U, XLONG_V, ... to coords
through the 'coordinates' attribute of U and V, and they survived to
break MarsFiles -t (tuple.index('lon')).

tests: test_marsfiles.py 14 passed, 4 skipped.
Recent planetWRF writes LTST, the true local solar time of the model
sun. Its departure from mean local time (UT + lon/15) is one scalar per
frame, the equation of time, which on Mars drifts by tens of minutes
over the year (-3 min at Ls 68 in the test file).

- MarsFormat (marswrf): when LTST is present, derive eot_offset(time)
  [hr] = median over the grid of LTST - (UT + lon/15), fill restart
  frames from neighbours, and save it. LOCAL_TIME is the naive clock
  and is not used.
- time_shift_calc: new keyword eot_offset [hr] added to the longitude
  shift (one value per call, i.e. per sol).
- MarsFiles -t: if the diurn file carries eot_offset, shift sol by sol
  with that offset and report that the output axis is true local solar
  time; otherwise unchanged (mean local time).

Verified against hand interpolation including the offset at four
longitudes for -t "12 15" (exact to 0.01 K).

tests: test_marsfiles.py 14 passed, 4 skipped; test_marsformat.py 9 passed.
…ck in MarsInterp

MarsFormat (marswrf):
- zfull (height above local surface) now comes first from Z_PHY, the
  model's own layer-centre geopotential height above the areoid,
  minus HGT. Without HGT (reduced files) the surface height is
  estimated hydrostatically from the bottom layer (Z_PHY, P_PHY, T_PHY,
  PSFC; 1.6 m rms against HGT on the test file) and saved as zsurf.
  PH/PHB/HGT remain the fallback for older files.
- ucomp, vcomp, w are taken from U_PHY, V_PHY, W_PHY whenever present:
  these are the model's own values averaged to mass-point layer
  centres, whereas U, V, W sit on the Arakawa-C staggered locations
  (W on layer interfaces). The 2-point destaggering is kept as the
  fallback for files that only carry U, V, W.

MarsInterp: for -t zstd, when no DDDDD.fixed.nc exists, use zsurf
carried in the input file itself (first frame if time-dependent),
instead of exiting. MarsFormat writes zsurf for converted models, so
zstd (height above the areoid, topography visible in zonal means)
now works on planetWRF output without a fixed file.

tests: test_marsinterp + test_marsfiles 27 passed, 4 skipped;
test_marsformat 9 passed.
…rp copy fix

New option -stag / --keep_staggered (marswrf only). In addition to the
mass-point ucomp, vcomp, w (from U_PHY, V_PHY, W_PHY), write the winds
on their native C-grid locations:
  u_stag(time, pfull, lat, lon_u)   lon_u from XLONG_U, 0-360, sorted
  v_stag(time, pfull, lat_v, lon)   lat_v from XLAT_V
  w_half(time, phalf, lat, lon)     W on layer interfaces
  phalf3D(time, phalf, lat, lon)    PF_PHY, interface pressure
  zhalf(time, phalf, lat, lon)      ZF_PHY minus the surface, if present
w_half, phalf3D and zhalf sit on CAP's phalf dimension and go through
MarsInterp like any interface field; u_stag and v_stag are for the
user's own scripts (CAP plots only the mass-point grid) and are pruned
from interpolated files. Verified against the raw U and V.

MarsInterp: the copy-through branch tested
  'pfull' not in dims or 'phalf' not in dims
which is always true, so a variable on pfull with non-standard
horizontal dims crashed the run (cannot find dimension pfull). Now
such variables are skipped with a message.

MarsFormat prune rule for vertical-only metadata now recognises lat_v
and lon_u as horizontal dimensions.

tests: test_marsformat + test_marsinterp + test_marsfiles 36 passed,
4 skipped.
FV3_utils:
- cgrid_div, cgrid_curl: divergence and relative vorticity formed on
  the Arakawa-C grid from u_stag/v_stag (MarsFormat -stag), i.e. as
  the model itself forms them. Divergence at mass points from the
  cell-face winds; vorticity at cell corners then averaged to mass
  points (pole corners excluded). Verified against a hand calculation
  on the raw WRF arrays to 1e-6 relative.
- fms_Z_calc(rgas=191., g=3.72): the gas constant and gravity are now
  arguments, scalar or arrays shaped like T, instead of hard-coded.
  Defaults reproduce the previous results exactly.
- mass_stream(g=3.72, a=3400e3): same, for the stream function.

MarsVars:
- set_planet_constants(f): rgas/Rd, Cp, g, planetary radius and
  R/cp from the file when it carries the planetWRF attribute set
  (R_D, CP, G, RADIUS, all four required); the 3D fields rgas3D, cp3D
  (composition-dependent) and g3D (variable gravity) take precedence
  when present. Files without them keep the Ames defaults and the
  original theta exponent R/(M_co2*Cp).
- div and curl use the C-grid winds when u_stag/v_stag are in the file
  (native-level files from MarsFormat -stag), else the A-grid
  spherical_div/curl as before. R = 3400 km replaced by the planetary
  radius from the file where known (3389.92 km for planetWRF).
- compute_zfull/zhalf and msf pass the file constants through.

MarsFormat (marswrf):
- rgas3D = 8314.46/MW_AIR_3D and cp3D = rgas3D/RCP_3D [J/kg/K] when
  the composition fields exist (8314.46/MW reproduces R_D exactly).
- g3D = G*(1 - Z_PHY/RADIUS)**2 when DO_VARIABLE_GRAVITY = 1
  (planetWRF g_of_phi with phi = G*Z_PHY); checked to 1e-4.

CAP's hydrostatic zfull with the file constants now tracks the model's
own Z_PHY - HGT within 0.4% up to 145 km on the test file.

tests: test_marsformat, test_marsinterp, test_marsfiles, test_marsvars
45 passed, 4 skipped.
Extra key = value fields after the positional Axis Options
(X, Y, cmap, scale, proj) are now collected for 2D plots:
  nlev   = 41              number of filled contour levels when
                           Cmin, Cmax is a pair (default 21)
  cticks = 120,160,200,240 colorbar tick values
  cfmt   = %.1f            colorbar tick label format
Existing templates are unaffected (fields were previously ignored;
defaults unchanged). Documented in the generated template.

tests: test_marsplot.py 15 passed.
Documents the MarsFormat -gcm marswrf path: quick start, the source
priority for every derived quantity (temperature, pressure, time, Ls,
eta, winds, heights, local time, gas constant, gravity), the -stag
option, planet constants in MarsVars, and true local solar time.
@mirfjc
mirfjc requested review from falconstryker and rurata August 27, 2026 01:51
@rurata
rurata merged commit bcc9aca into devel Aug 28, 2026
63 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants