Refactored soca -> cice - #1246
Conversation
Introduces src/soca/IO/soca_io_mod.F90 with soca_io_file_reader / soca_io_file_writer types (init/enqueue/commit/close), and routes all state, geometry, and balance-coefficient I/O through it. PE 0 calls nf90_get_var / nf90_put_var; mpp_broadcast and mpp_gather distribute / collect data on the geometry's pelist. Removes all fms_io_mod usages (register_restart_field, restore_state, save_restart, read_data, file_exist, field_exist, fms_io_init, fms_io_exit) from soca, including the global_soca_geom_counter shim. fms_init / fms_end are retained (subsystem init only, no I/O). Resolves the GDAS Tsnz_h shape mismatch by building nf90_get_var count from the file's actual dim sizes, not the caller's buffer rank. Per-domain state writes (ocn/sfc/ice/wav/bio) and reads now use the new module; the FMS code paths are removed in the same change rather than gated. Refs #1125. Obsoletes PR #1241.
Under FMS, the writer would implicitly add the .nc suffix, so test YAMLs and downstream readers reference 'ocn.<exp>.<typ>.<date>.nc'. The direct-netCDF writer in soca_io_mod is faithful to whatever string it's given, so files were landing without the extension and no consumer could find them.
Six soca_add_test() calls listed test_socs_parameters_diffusion as a dependency (should be test_soca_). Because the named test does not exist, ctest treated the dependency as void and ran the variational tests concurrently with parameters_diffusion under -j>1, causing races on the diffusion calibration outputs.
Drops the SOCA_IO_SERIAL / SOCA_IO_PARALLEL mode bifurcation (and the method= arg + io.method config knob) and rewrites the reader to do per-PE strided nf90_get_var, mirroring FMS mpp_io's domain-decomposed read pattern. Each PE opens the file once via NF90_NOWRITE and reads its own compute-domain tile via nf90_get_var(start, count) -- the same shape as fms_io's READ_RECORD_ for fileset.NE.MPP_SINGLE. Adds a module-level read cache mirroring fms_io's get_file_unit + files_read(i)%var(j) tables: - one nf90_open per (PE, filename) for the whole run - cached (varid, file_ndims, middle_dim sizes) per (file, var) - soca_io_close_all() flushes the cache; wired into soca_geom_end Caller updates: drop method= from reader/writer init() and the soca_io_method_from_config helper from soca_fields_mod / soca_geom_mod / soca_balance_mod. 1-deg 20-mem LETKF read phase (LocalEnsembleDA before solver ctor): baseline (FMS): 3.91 s before this commit (no cache): 13.92 s with this commit: 3.97 s Outputs are bit-identical.
…cache) Adds a second reader_commit implementation -- PE 0 nf90_get_var the global field + mpp_broadcast -- alongside the existing per-PE strided path. Both paths share the file-handle + var-metadata cache. Select via SOCA_IO_READ_MODE=broadcast|strided (default broadcast); the env var is latched on first reader_commit. Why broadcast as default: in DA cycling the page cache is always cold, since each cycle reads files that the previous cycle / model run just wrote. On cold cache the broadcast path (1 sequential read on PE 0, kernel readahead happy) beats strided (8 concurrent strided readers thrash the prefetcher and pay 8x the open-syscall cost). 1-deg 20-mem LETKF, rancor, taskset -c 0-7, cold cache (sample 1): broadcast: 13.5 s read strided: 30.5 s read Warm cache reverses the result (strided 3.9 s vs broadcast 13.3 s) but cycling never sees warm cache. On a parallel filesystem (Lustre/GPFS) or multi-node setup the strided path may win again -- toggle the env var to test.
Drops SOCA_IO_READ_MODE in favor of a yaml block on the geometry
config:
geometry:
io:
read mode: broadcast | strided
If absent the module-level default (broadcast) is kept.
soca_geom_init calls soca_io_read_mode_from_config(f_conf) once at
startup; the choice latches for every subsequent reader_commit.
Unrecognized values abort so typos surface immediately.
Annotates testinput/letkf.yml's geometry block with the two read-mode options so a contributor reading the canonical letkf example sees both choices without having to hunt for the soca_io_mod docstring.
The write path emits FMS-style auto-numbered dim names (xaxis_N / yaxis_N / zaxis_N / Time); the file-header docstring incorrectly claimed xh / yh / Time. Correct the comment -- no code change.
Writer: - enqueue holds pointers to caller buffers instead of allocating and copying (mirrors FMS register_restart_field contract). Compute-slice extraction moves from enqueue to commit, so peak per-writer memory drops from sum(var_bytes) to one (nx_c x ny_c) tile. - commit uses the 3D mpp_gather overload for 3D vars: one collective per var instead of nlevels, and no per-level gbuf3d(:,:,k) = gbuf2d memcpy on root. - Reuse gbuf3d / tile3 across 3D vars when nlevels matches. Reader: - commit_reader_strided hoists tile2 out of the per-var loop and reuses tile3 / tile4 when trailing dims match. - Dead tile3/tile4 zero-inits removed (read_var_strided fully fills the buffer). read_var_strided: fixed-size stack arrays for st / ct instead of per-call heap allocation. put_axis_coord_data: reuse idxbuf across axes. Cleanup: - Remove dead 'use mpi' and 'use fckit_configuration_module'. - Remove dead cartesian_axis parameter on writer_enqueue_1d (it was stored but never written to the netcdf output). - Drop unused nprocs out-param of mpi_pelist. - soca_io_close_all: use ncc on nf90_close instead of silently discarding the status. - commit_reader_scatter: matching dead-init cleanup; comment updated to flag it as pending parallel-ensemble I/O exercise. - Trim verbose / redundant comments.
Drop dangling references to the removed soca_io_read_mode_from_config / 'io.read mode' YAML key in soca_geom_mod and the letkf test input. The selector was already gone from soca_io_mod; without these the PR fails to compile.
- soca_io_mod: allocate gbuf2d / gbuf3d as 1x1 dummies on non-root so the actuals passed to mpp_gather are always allocated (assumed-shape dummy requires it). Pattern matches commit_reader_scatter. - soca_io_mod: refresh stale 'Data is copied in' comment to reflect the pointer-based enqueue semantics (caller buffer kept alive/unmutated through commit; actual must satisfy TARGET association rules). - soca_geom_mod: TARGET on the 'self' dummy in soca_geom_init, soca_geom_init_fieldset, and soca_geom_write so the allocatable components (self%lonh, self%lat, ..., self%mask2d*) are valid pointer targets for the soca_io reader/writer enqueue. TARGET on the local fieldData / fieldDataVars too. - soca_fields_mod: TARGET on h_common and on the local vars(:) wrapper array so vars(n)%data sections used in enqueue are valid pointer targets. - soca_balance_mod: TARGET on local kct for the same reason.
…e checks Three related fixes in the direct-netCDF reader path: 1. read_var_strided previously hardcoded "last Fortran dim is trailing time" and set ct(file_ndims)=1 unconditionally. For a file with spatial-only layout (Temp(z,y,x) with no leading Time), this read only level 1 of z into a multi-level destination, leaving the rest uninitialized -- showing up as garbage scale values in the vertical diffusion calibration. Now discriminate file_ndims == dst_rank (no time, fill every dim from the file) vs file_ndims > dst_rank (trailing time + middle squeeze, e.g. CICE Tsnz_h's nksnow=1). A total-element-count check catches silent partial-fills and dim-size mismatches (e.g. file z=75 vs destination z=25). 2. Drop the read_cache / cached_open / soca_io_close_all machinery. Each reader_commit nf90_opens, reads all enqueued vars (with inline inq_varid + inquire_variable + inquire_dimension; microseconds), and nf90_closes. Holding NetCDF4 handles open across commits was bloating LETKF per-task memory by ~MB-to-GB scale (HDF5 metadata + chunk caches per open file). Restores per-task max to match develop. 3. Add check_buf_1d / check_buf_2d helpers called from every reader/writer enqueue. Catches unallocated and wrong-sized caller buffers at the enqueue site instead of silently storing the pointer for a later get_var/put_var to mishandle.
# Conflicts: # src/soca/Fields/soca_fields_mod.F90
|
I'm going to run a few cycles with the ensembles using this postprocessing and report back here. In the meantime, this is mostly ready for review, I don't anticipate major changes. If it looks painful to review, just know that I've been working on this on and off for half a year, and I'm with you 😆 . |
CICE6's linear_itd aborts when per-cat thickness `hi = vicen/aicen` falls
outside `[hicat[k], hicat[k+1]]`. Diagnostic scan of the existing
production-grid restarts (datatoanalyze/, 7 dates) found ~5% of ice cells
in every restart had at least one populated cat with out-of-bin hi --
~5,000-16,000 cells per date. CICE6 was crashing on these.
Two bugs at play:
1. **Default `category bounds` was CICE5 layout** with bounded top
(`{0, 0.6445, 1.391, 2.470, 4.567, 9.334}`) whereas CICE6 GFSv17
uses `{0, 0.64, 1.39, 2.47, 4.57, 1000}` (open top). Silent
mismatch produced restarts CICE refused to read. Make
`category bounds` a RequiredParameter -- users must spell out the
layout matching their CICE version.
2. **The standalone min-vice cleanup block** in PostProcessIce.cc
(now removed) redistributed mass from dropped cats to surviving
cats proportional to surviving aicen, without respect to bins.
This shifted survivors' per-cat hi outside their bins, producing
the CICE6 abort cases. The dd Python reference (mod_cice6_utils.py
adjust_thkncats_aice) doesn't have a separate min-vice cleanup --
its scipy.optimize.minimize call uses the bin edges as variable
bounds, so survivors are guaranteed in-bin. We rely on the rebin's
existing `ainMin = 1e-8` clip + Step 7 renormalisation (which
preserves per-cat hi), and add a final per-cat bin clip as
defence in depth:
const double hClipped = std::min(std::max(h, hLo), hHi);
if (hClipped != h) vicen[k] = aicen[k] * hClipped;
Aicen is left alone (Σaicen preserved); a small amount of mass is
lost to the clip, much smaller than the test failures the bug
produced.
Drops the now-unused `min cat ice volume` yaml parameter.
Test fallout: the 4 soca2cice_new* testrefs shift slightly (min-vice
cleanup is gone). The freeboard test shifts hardest (hice 0.4153 ->
0.4044, because the cleanup was undoing some of freeboard's
ice-volume growth). Test yamls updated to set `category bounds`
explicitly. 9/9 postprocice tests pass.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Log::info() writes only from task 0 by oops convention, so all the
counter / max-magnitude diagnostics in runPostprocess were silently
reporting only rank 0's piece of the grid -- misleading on multi-task
runs.
allReduceInPlace each counter (sum) and each magnitude (max) on
geom_.getComm() before the log lines, so what gets printed is the
global total / global max. Covers:
- rebin_visited, rebin_aicen_mutations (sum)
- rebin_max_delta_aicen_all (max)
- rebin_failures, freeboard_failures (sum)
- bin_clip_slots (sum)
- bin_clip_max_abs_dh (max)
- seedNewIce fallbacks (sum, separately because
seedNewIce runs after the main loop)
Also adds the bin-clip diagnostic itself (was uncounted): per-cat slots
clipped + max |dh| in metres. Only printed when at least one slot was
clipped.
9/9 postprocice tests still pass.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Log::warning fans out to every MPI task, so the per-rank warnings were printing N identical copies of the global-summed counts. Switch the ITD-rebin-failure, freeboard-failure and noice->ice Tfrz-fallback messages to Log::info (rank-0-only). These are diagnostic counters about non-fatal events the code already handles, so info is the right level. Mention the bin clip in the rebin-failure message so the log explains how out-of-bin cells get repaired. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
- Replace the 'min-vice cleanup' step in the per-cell pass with the new 'per-cat bin clip' step (the cleanup is gone; the clip is the guarantee that the restart is CICE-readable). - Drop the 'min cat ice volume' configuration entry (parameter no longer exists). - Mark 'itd: category bounds' as required in both the schema block's inline comment and the closing 'required fields' bullet. Switch the example value to the CICE6 GFSv17 layout ([0, 0.64, 1.39, 2.47, 4.57, 1000]) and document the CICE5 alternative; explain that the previous silent default was a footgun. - Tidy two surrounding references to 'min-vice cleanup' that pointed at the dropped step. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
|
I realized that I made a mistake in my previous legacy vs this branch vs Python scripts comparison computations (I forgot again that hi_h and hs_h in the history files are effectively volumes, not thicknesses). I have now cleanly reran the comparisons, so I'm removing my previous comment and reposting the updated results here. Note that the results are now consistently good for Dmitry's scripts, without weird artifacts around ice edge in my previous incorrect comparisons. The conclusion is the same: no degradation compared to Dmitry's scripts, significantly better than legacy Fortran soca2cice for ice thickness and snow depth. PostProcessIce: production-grid validationSide-by-side comparison of the new C++/atlas The comparison was done on For each cycle we hold the same DA analysis (from GFSv17 with artificially increased ice thickness and snow depth to test insertion of those variables) and the same CICE background fixed, and compare the resulting CICE restart written by each of the three pipelines. Aggregate cell-level RMSE vs analysisMask: cells where
Per-category ICE2ICE perturbation vs backgroundMask: per-cat slots where both
PlotsThree production-grid comparisons for
All 7 dates exhibit the same qualitative behaviour - the new path is consistently close to the analysis on ice thickness and snow (significantly better than legacy). |
|
I pushed a few updates to this branch yesterday that fix some issues with ice thicknesses out of category bounds (thank you @DmitryDukhovskoy for you help in identifying the issues!). I have now successfully run a few cycles of a cycling experiment with an ensemble where the new soca2cice is used for updating CICE restart for the deterministic forecast after 3DVar and for the ensemble forecasts after ensemble recentering. |
|
There is a lot here! Some thoughts on reading through, please take as you wish. Dynamics are missing (velocities and stresses). Weighting donors and ice-->ocn as aice vanishes are some options. There can be melt ponds and snow on ice in the same cell. One end-to-end test may be that DA shouldn't change much if the background is close to the analysis. I don't know if you have tests like this, but it's good sanity check. It's one reason I mention an issue of zeroing ponds if snow changes. Some variables can have subtle meanings (CICE-Consortium/CICE#1096 (comment)) Maybe the self cell should also be included in the neighborhood of donor cells. It wouldn't change anything in the noice2ice branches |
|
@NickSzapiro-NOAA thank you for looking through this and for your comments! I'll think about that and get back to you in probably a couple of weeks (code sprinting next week). |
DmitryDukhovskoy
left a comment
DmitryDukhovskoy
left a comment
There was a problem hiding this comment.
I looked through the changes in the SOCA code. The insertion logic closely follows the Python scripts, except for the redistribution of sea-ice state variables across the ice thickness categories. This difference should have a small impact, provided it does not produce category-mean thicknesses (vicen/aicen) that fall outside the prescribed thickness bounds. I don't see any obvious issues.
|
Thanks for putting together such a thorough update. The overall workflow was easy to follow, and the validation tables and figures were especially helpful in understanding how the new approach compares with the legacy method and Dmitry's scripts. I had two small comments on PostProcessIce.cc. First, I may be missing something, but a_aice appears to be created as a writable Atlas view of the input "analysis" field and is then clamped in place. Since "analysis" is passed as a const State& , would it be safer to copy aice into a local vector, similar to a_hice_vec and a_hsno_vec, before applying the bounds? Also, in the ICE2ICE path, it looks like "vicen" can be redistributed while the existing "sice" profile remains unchanged. Is the idea that any added ice volume takes on the salinity of that category? I was also wondering whether total salt is conserved through this redistribution. |
Hi @mjkagnes123, the sea ice salinity profile and enthalpy are computed only when the ice column has essentially no existing salt (sice ~0) or enthalpy is not within the correct bounds. This is always the case for newly formed ice (noice --> ice), but it is usually not true for the ice2ice pathway. When sice ~ 0, the model has to estimate the salinity profile, since there is no existing salt distribution to work from. The estimate is based on the BL99 (Bitz & Lipscomb, 1999) approach, as implemented in CICE4. The direct insertion approach does not guarantee conservation of salt, mass, or energy. |
guillaumevernieres
left a comment
guillaumevernieres
left a comment
There was a problem hiding this comment.
I did not review the code but I'm testing this within the gfs/ufs. It does as advertized!
Thanks @shlyaeva
mjkagnes123
left a comment
mjkagnes123
left a comment
There was a problem hiding this comment.
I really appreciate the work you did on this. The tables and figures were very useful in showing how this is better than the way of doing things. The improvements over the legacy approach are clear. Dmitrys explanation also helped me understand the part about treating salinity. My comment about 'a_aice' is still something I think could be improved. It is not a big deal and it does not stop this from being approved. I am approving this after I looked at 'PostProcessIce.cc'. The improvements over the legacy approach, like the salinity treatment are good.
Dooruk
left a comment
Dooruk
left a comment
There was a problem hiding this comment.
@shlyaeva this is quite a comprehensive overhaul of the current approach and I really appreciate the README. It will take us a while to implement this into our workflow. Since the old approach is still available for a while, I don't want to hold this back.
Also, thank you @mjkagnes123 for your thorough review!
@mjkagnes123 thank you for catching this, I've addressed it in the latest commit. |
|
Thank you everyone for the reviews! I'm going to merge in its current form once the tests pass, and we can improve on it in follow-up PRs as needed. |
|
It's very possible that I face planted again when modifying the configuration (volume vs thickness) ... But I don't think so. I'm finding a linear growth of seaice volume for both poles when using the new soca2cice. @Dooruk , @mjkagnes123 , have you guys had time to test? |
Edit: It's not as simple as what I wrote above, but I'd be curious to compare notes once you test this with your system. |
7 participants
Description
Sea-ice analysis postprocessing in C++/atlas (
PostProcessIce): replaces the FortranSoca2Cicevariable change with a C++/atlas postprocessor that projects an aggregate SOCA ice analysis onto a per-category CICE restart. The new algorithm is largely based on @DmitryDukhovskoy's Python scripts to insert ice concentration, thickness and snow depth (https://github.com/DmitryDukhovskoy/RTOFS_utilities/blob/gaea_pub/prepare_cice6/insert_iconc_ithkn_cice6_restart.py and https://github.com/DmitryDukhovskoy/RTOFS_utilities/blob/gaea_pub/prepare_cice6/insert_hsnow_cice6_restart.py). Validated on 7 production GDAS cycles spanning different seasons.What's new
soca::PostProcessIcein src/soca/PostProcess/,with a single public entry point:
Reads the CICE background restart, applies the per-cell pass, writes the postprocessed restart in update mode, returns an aggregate-ice State
{aice, hice, hsno}matching what was written. Used by the standalonesoca_postproc.xapp, byAnalysisPostproc.h's ensemble loop, and by gdasapp's increment handler.Per-cell pass (see README for details): case dispatch (LAND / ICE2NOICE / NOICE2ICE / ICE2ICE), ITD rebin (on by default), aicen-weighted snow distribution, optional freeboard enforcement, mass-conserving min-vice cleanup. Pure column-physics helpers (BL99 enthalpy, BZ99 salinity profile, thickness-category solver, freeboard) live in
IcePhysics.h/.ccwith their own unit tests inTestIcePhysics.cc.New-ice thermo seeding: cells that transitioned from
bg_aicen=0tonew_aicen>0get Tsfcn from the area-weighted mean of a KDTree donor with any ice (global lat/lon tree, donor data gathered once via a sparse halo exchange); sub-surface qice/sice synthesized from CICE physics at the ocean freezing point.soca_postproc.xstandalone application(src/mains/Postproc.h): takes
background+increment, formsanalysis = bg + incr, runsPostProcessIce. Also accepts an explicitanalysisblock.Dedicated CICE restart writer on the soca::State side:
soca::State::writeCice(cfg)calls a new Fortran entrysoca_fields::write_cicethat wraps the update-mode CICE writer (byte-copy input restart → output, overwrite only modelled variables, ~40 unmodelled CICE vars pass through).fields_metadata.yml: new<LEVEL>placeholder for templated per-layer entries (CICE per-layer restart varsqice00N,sice00N,qsno00N); combines with<CATEGORY>to expandncat × ice_leventries automatically. Per-cat iceio nameupdated to restart naming (aice<N>,vice<N>,vsno<N>- dropping the_hsuffix which is the CICE history convention; restarts don't carry it).AnalysisPostproc.h: ensemble loop migrated to callPostProcessIce::postprocessper member.Soca2Cicevariable change no longer used here.Tests
Four end-to-end tests on the 72×35 grid (one per major code path) plusa pure unit test on the column physics:
test_soca_soca2cice_new- baseline (rebin on)test_soca_soca2cice_new_freeboard- freeboard enforcementtest_soca_soca2cice_new_seed- new-ice seeding + thermotest_soca_soca2cice_new_bgfallback- analysis variables resolver(only aice analysed, hice/hsno fall back to background)
test_soca_icephysics- pure C++ unit test of the column-physics helperstest_soca_postprocice_vs_soca2cice- Fortran-vs-C++ regression(
Soca2Ciceoutput compared bit-wise toPostProcessIceoutput;budget = 2 differing cells at 1e-10 tolerance)
test_soca_ensanpproc- ensemble postproc end-to-endAll pass. The legacy Fortran
Soca2Cicepath is still available throughsoca_convertstate.xfor the regression test; it can be deprecated in a follow-up PR.Configuration
A minimal yaml stanza:
Full schema with defaults in src/soca/PostProcess/README.md.
Checklist