PISM ice-sheet grids and coupling links for OASIS, behind ice.enabled - #57
Merged
Conversation
The ISM-mapper is an OASIS component, so a PISM domain has to appear in grids/masks/areas next to the atmosphere and ocean grids. New PISM grid class alongside ORCA, registered in the factory, plus an ice section in the config and an append step in the writer. Off unless ice.enabled is true, so existing configs are unaffected. Centres come from the PISM file's own lat/lon; corners are inverse-projected from the cell corner (x,y) with the projection in its mapping variable. Areas are spherical excess over two triangles, not a plane formula: the domain wraps the pole, so cells on the +/-180 meridian would get a lon span of ~2*pi. Verified on antarct_cr 8 km (761x761): area total 36.03 against 36.28 x10^6 km2 from the plane area corrected by the polar stereographic scale factor. Cell areas run 67.0 km2 at the pole to 53.8 at the outer corner against a nominal 64, right for lat_ts=-71. All cells wind counterclockwise seen from outside the sphere, which SCRIP needs and the area check cannot detect. Also makes the FESOM grid export opt-in via ocean.write_oasis_grid, default off. FESOM writes feom itself at runtime in dist order; what we write here is mesh order, so the two disagree. NOTE this changes default output for every FESOM config: feom no longer appears unless the config asks for it.
Two links cover the whole ISM field set: a link is one (source, target, map), not one per field. Everything OIFS sends the ISM-mapper is atm -> ice, and plit is ice -> atm. method picks the ice -> atm map. bilinear reproduces EC-Earth. conserv is likely better for plit: it is a 0/1 mask on 8 km going to ~100 km cells, and what suorog thresholds against ECE_LANDICE_THRESH is an area fraction, which conservative gives directly and a centre-based map does not. Left as bilinear by default so the reference behaviour is what you get without asking.
Three links, not two: runoff is on the runoff-mapper's atmosphere grid (R096), not atma, so it needs its own. The SCRIPR letter is the SOURCE grid type. ismp is logically rectangular, so the ice -> atm link is LR, not the D it had. BILINEAR and CONSERV both abort on this grid, so the method profiles are now gauswgt and distwgt. 381 of the 761x761 cells have a corner longitude spread above 180 deg along the polar seam, and the pole cell's corners encircle the pole, which breaks SCRIP's enclosing-quadrilateral search and its per-cell longitude line integral. Consequence worth knowing: plit reaches suorog as an interpolated value, not the area fraction ECE_LANDICE_THRESH is written for. Tested on compute (3 nodes, one rank per link) against ocp-tool's own TCO95 CORE3 output: rmp_A096_to_ismp_BILINEAR 69 MB ok rmp_R096_to_ismp_GAUSWGT_25 122 MB ok rmp_ismp_to_A096_GAUSWGT_25 44 MB ok rmp_ismp_to_A096_DISTWGT_4 34 MB ok ismp->A096 BILINEAR, CONSERV OASIS abort rc=1
OASIS restarts hold the source field on the source grid and remap at read time, so this writes only <prefix>_plit on the ice grid. The atmosphere fluxes are safely zero-filled by NNOREST; a mask is not. With LAG the receiver takes its first value from the restart, and an empty mask strips the glacier tiles from the atmosphere for a coupling interval, which shows up as degenerate skin temperatures and then a floating invalid in ece_updclie_cpl. IceRegionConfig gains prefix, since the ISM-mapper field prefix (antar) is not the OASIS grid name (ismp). Defaults to the grid name. Verified on antarct_cr 8 km: 221566 ice cells, 14.37e6 km2 against the 13.7-14.2e6 of Antarctica including shelves.
The FESOM export is off by default because FESOM writes its own feom entry at runtime with the distribution reordering, so the pool build has no business writing one. Weight regeneration is the other case. The remapping weights for a new mesh are built against the ocean grid in the OASIS files, so if it is not there they cannot be built, and the leg starts on the previous mesh's weights and dies in MCT with a column-count mismatch. write_oasis_grid overrides the template so the caller that goes on to build weights can ask for it, rather than the regeneration path deciding on everyone's behalf.
For the ice sheet's basal melt and discharge. RnfA is a regular 512x256 lat-lon grid, so unlike the ice to atmosphere case the target is not the awkward side, but the source is still pole-centred, hence LR and a neighbour-based map.
This reverts commit 5d384dc.
…arine CCN
WHY. The shipped climate.v020/<res>/month_dms is Lana et al. (2011). DMS-Rev3
(Hulswar et al. 2022, ESSD 14, 2963-2987) rests on 865,109 observations against L11's
47,313 -- an ~18-fold increase -- with revised filtering, dynamic biogeochemical
province boundaries and a new smoothing algorithm. Globally it is only 4 % lower, but
the paper is explicit that the largest revisions are in the polar oceans, which is the
only region this campaign is looking at.
MEASURED ON THE MODEL GRID, SO 65-45S, L11 -> Rev3:
annual 1.93 -> 2.24 (+16 %)
DJF 4.23 -> 4.36 ( +3 %)
SON 1.34 -> 2.02 (+51 %)
November 2.18 -> 4.24 (+95 %)
global max 53.25 -> 26.56 (-50 %, the paper's "reduced patchiness")
So Rev3 does not raise the summer peak, it BROADENS it -- the bloom starts a month
earlier. The model's SO shortwave CRE error is worst in SON (-7.03 pp) and DJF (-6.94),
so the upgrade lands in the season the old climatology was weakest. ACP 23, 15181
(2023) independently finds the oceanic DMS dataset moves SO summertime emissions by
112 % against 50-60 % for the transfer velocity, i.e. the dataset is the larger lever.
UNITS WERE DERIVED, NOT READ. DMS_REV3.nc carries no attributes at all. Integrating
the shipped DMSflux globally gives 28.19 Tg S/yr (DMSfluxSI 28.36) against a published
27.1-28.1, which fixes the flux as mol m-2 d-1 and the concentration as nmol/L; every
other assumption is out by orders of magnitude. Only the concentration is ingested --
the shipped flux uses Rev3's own wind climatology, whereas the model must form the flux
online from its own wind, which is quadratic in speed. The shipped flux is kept as a
validation target.
METHOD. dms_rev3_to_grib.py transfers Rev3 onto the model's reduced Gaussian grid by
nearest VALID neighbour in 3-D (so the dateline and poles cannot give a spurious
neighbour). Nearest rather than bilinear because 1 deg source against ~0.9 deg target
is matched resolution -- bilinear would need valid-weight renormalisation at every coast
for no gain, and would bleed land NaNs into coastal ocean. Values are overwritten into
the SHIPPED month_dms messages rather than building GRIB from scratch, so grid, dataDate,
paramId and all headers stay byte-identical to what the reader already accepts. The land
sentinel comes from the template, preserving the model's own mask; a Rev3 hole at a
template ocean point would read as land, so L11 is retained there and the count reported.
For TCO95: zero holes, 12040 land points held at -99.
icmcl_dms.py then interleaves the result as the eighth field of each date group, which
ece_updclie_climr.F90 requires -- 84 -> 96 messages, 12 dates x 8 fields, verified.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Groundwork for coupling PISM to AWI-ESM3 through the EC-Earth4 ISM-mapper. The ISM-mapper is an OASIS component, so a PISM domain has to appear in
grids.nc/masks.nc/areas.ncand needs remap weights, and ocp-tool is where both already live.Draft: the
plitremap choice below is unresolved, and the esm_tools side that would exercise this is not written yet.What is here
ocp_tool/grids/pism.py— aPISMgrid class with the same five methods asORCA(cell_latitudes,cell_longitudes,cell_corners,cell_areas,cell_masks), registered ingrids/factory().ocp_tool/config.py—IceConfig/IceRegionConfig/_load_ice_config, gated exactly likepaleo.ocp_tool/oasis_writer.py—_append_pism_grid_to_oasis_files, alongside the existing runoff / FESOM / AMIP appends.ocp_tool/oasis_weights.py—awiesm3_ismp_links(). Two links cover the whole ISM field set, since a link is one(source, target, map)and not one per field.configs/TCO95_CORE3_ICE.yaml—TCO95_CORE3.yamlplus theicesection.Centres come from the PISM file's own
lat/lon. Corners are inverse-projected from the cell corner(x, y)using the projection in itsmappingvariable. Areas are spherical excess over two triangles rather than a plane formula: the domain wraps the pole, so cells on the ±180 meridian would otherwise get a longitude span of ~2π.Two things reviewers should not have to find
1. This changes default output for every FESOM config.
feomwas being appended to the OASIS files unconditionally. FESOM writes its ownfeomentry at runtime in domain-decomposition order, and what ocp-tool writes is mesh order, so the two disagree —permute_feom_to_runtime.pyin esm_tools exists because of that difference. It is now behindocean.write_oasis_grid, default off, sofeomno longer appears unless a config asks for it. Setwrite_oasis_grid: trueto restore the old behaviour.Everything else is additive and inert without
ice.enabled.2. The
plitremap is an open decision.awiesm3_ismp_links(method=...)defaults to"bilinear", reproducing EC-Earth. I think"conserv"is actually right for the ice→atm direction:plitis a 0/1 mask on 8 km cells landing on ~100 km cells, and what OIFSsuorogcompares againstECE_LANDICE_THRESHis an area fraction. Conservative gives that directly; a centre-based map gives an interpolated value that only resembles one. Thesuorogcomment already describes it as a "fractional ice sheet mask (plit, 0-1 after regridding)". Left as bilinear so the reference behaviour is what you get without asking, but worth settling before weights are generated in anger.Verification
Grid geometry checked three ways on
antarct_cr8 km, 761×761, since a wrong grid fails as silently-bad weights rather than a crash:lat_ts = -71, where the scale factor is 0.973 at the pole and 1.089 at the grid edge.A full TCO95/CORE3 run completes and puts
ismpin all three OASIS files at 761×761 withcrn_ismp = 4, bit-identical (max|diff| = 0onlon,lat,clo,cla,srf) to an independently generated reference, withfeomcorrectly absent.Noted, not fixed
_append_fesom_grid_to_oasis_filesguards its optional import withexcept ImportError, butimport pyfesom2raisesAttributeErrorwhen cmocean meets a matplotlib new enough to have droppedmatplotlib.cm.register_cmap, so the whole run dies instead of skipping. Left alone deliberately: widening the catch would turn a loud failure into OASIS files silently missing the ocean grid, which fails far more confusingly later. It is also moot on the default path now, since that function is no longer called unless asked for.