Add two-way land-ice coupling to driver-moab

[rljacob]

Port GLC (MALI land ice) coupling to driver-moab so it works the same way as
driver-mct. This first version covers two-way lnd-glc coupling: per-elevation-class surface
mass balance and surface temperature downscaled to the ice sheet, and
per-elevation-class ice cover, topography and ice mask upscaled back to land.
It also adds glc mesh migration to the coupler, fractions, exchange, coupler
history, restart and budgets. Planar MPAS meshes are now supported in the MOAB
mesh instance, which the MALI Greenland mesh requires.

Ocean-glc thermal forcing and ice shelf coupling, glc-to-ocean and glc-to-ice
runoff mapping, and x2gacc restart are deferred to a follow-up pull request.
Their MCT code is left in place and inert.

Validated on chrysalis with gnu and the real Albany first-order velocity
solver, and on anlgce-ub22 with the SIA solver. The moab coupler history
matches driver-mct bit for bit for every field sent to the ice sheet.

[BFB]

---

What this adds

Two-way lnd<->glc coupling with elevation classes in driver-moab:

• Mesh migration. The MALI polygon mesh is sent to the coupler with
  iMOAB_SendMesh/ReceiveMesh, building on @iulian787 old branch
  iulian787/landice_moab_driver. Planar MALI meshes (sphere_radius = 0) get
  unit-sphere coordinates from latVertex/lonVertex, and areas are converted
  to steradians so they are consistent with every other coupler mesh.
• App ids. mbglid (MALI on component PEs), mbgxid (glc on coupler,
  COUPLE_MPASLNDICE), and mbintxlg/mbintxgl as read-map holder apps
  (LND_GLC_COU / GLC_LND_COU). Weight ids scalar_l2g, flux_l2g,
  flux_g2l.
• Mapping. Offline map files, the same seq_maps.rc entries driver-mct
  uses, loaded with iMOAB_LoadMapFile. Only GLOBAL_IDs (indexToCellID)
  matter, so existing lnd2glc/glc2lnd map files work unchanged. glc aream
  comes from the lnd2glc fmap area_b (arearead=2) as in driver-mct, with an
  aream=area fallback when land is absent.
• Component side. An ocean-style refactor: glc_import_mct/glc_export_mct
  now take raw real(:,:) arrays, and the moab path moves data through tags on
  mbglid with transpose twins. New glc_domain_moab sets lat/lon/area/mask/frac.
• lnd->glc. Per-EC land fields accumulate in l2gacc_lm following the
  prep_rof pattern. At the glc coupling step the averages are staged into the
  land-mesh tags for the duration of one batched lfrac-weighted normalized map
  (instantaneous values are saved and restored around it), then elevation-class
  vertical interpolation runs on arrays, followed by the qice area/aream
  pre-adjustment and global SMB renormalization. Results are written directly
  into the x2g tag names, which makes the merge a no-op as it is in mct.
• glc->lnd. Because per-EC weights differ per class, map_glc2lnd_ec_moab
  stages pre-multiplied numerator tags plus Sg_icemsk_num, does one raw map,
  and divides on the land side with reciprocal-multiplies matching mct's
  seq_map_avNormArr.

Inherited MCT glc mapper inits (seq_map_init_rcfile) segfault in driver-moab
because the coupler gsmaps are unpopulated, so they are replaced with
seq_map_mapinit stubs — the same treatment every other ported mapper received.

Matching driver-mct bit for bit

Three changes were needed to make the coupler history agree with driver-mct:

1. The per-elevation-class arrays in the vertical interpolation are stored in
   default real, replicating what mct's map_ice_covered does today (the
   demotion itself is documented in driver-mct: lnd->glc elevation-class downscaling demotes fields to single precision #8657). This is what makes x2g_Sl_tsrf
   agree, and it demonstrates that iMOAB's ApplyScalarProjectionWeights
   matches mct_sMat_avMult exactly for the same weights.
2. The averaged lnd->glc accumulation is no longer set back into the coupler
   land mesh tags, so the coupler history keeps showing instantaneous l2x
   fields the way driver-mct does. seq_hist_writeaux gains an optional
   matrix argument so the l2x1yrg auxiliary file, whose data lives in the
   accumulator rather than the tags, can still be written.
3. The glc->lnd elevation-class normalization is done as a reciprocal multiply
   with a zero denominator passing through as multiplication by zero, matching
   mct's seq_map_avNormArr idiom. This removes the ~1e-15 differences in the
   per-class frac and topo fields.

Validation

chrysalis, gnu + openmpi, Albany first-order solver (IGELM_MLI uses plain
MALI, so MALI_USE_ALBANY defaults to TRUE):

Test	Result
SMS.ne30pg2_r05_IcoswISC30E3r5_gis4to40.IGELM_MLI.mali-landiceIG	PASS
ERS_Ld3 same	PASS including bit-for-bit COMPARE_base_rest
SMS_Vmct / ERS_Vmct_Ld3 same	PASS
SMS_D (debug)	PASS

Coupler history, moab vs mct, same machine and compiler. With
glc_renormalize_smb='off', every x2g field is exactly bit for bit —
including x2g_Flgl_qice and x2g_Sl_tsrf. That is the substantive result:
the mapping delivers bit-identical forcing to the ice sheet.

Seven fields still differ at 1.2e-15 or below, all traceable to a single root
cause outside this change. MALI's own history, written after
li_core_initial_solve and before any coupling, is bit for bit on 15 of 16
fields; the exception is surfaceSpeed at 2.4e-15. So Albany's first velocity
solve returns a marginally different answer in a MOAB-linked executable, and
the residual g2x and topo01 differences are downstream of that. Two
identical moab runs are bit for bit across all 286 fields, so this is
reproducible rather than solver noise.

With the default glc_renormalize_smb='on_if_glc_coupled_fluxes', nine fields
differ instead of seven; the two extra are x2g_Flgl_qice and x2g_Sl_tsrf,
from the SMB renormalization global sums adding identical terms in different
orders under the moab RCB partition versus the mct gsmap decomposition. This is
the same reason driver-mct is not reproducible across PE layouts for these sums.

histaux_l2x1yrg, a 1-year IGELM_MLI run, exercises the new matrix
argument — the only path in this branch not covered by the tests above. It
writes 33 fields (11 elevation classes of Sl_tsrf, Sl_topo, Flgl_qice),
time_bnds = 0, 365, with physical values: Sl_tsrf capped exactly at 273.15 K,
Sl_topo monotonic across classes, and Flgl_qice exactly zero for the bare
land class. Note that GLC_AVG_PERIOD=yearly is required for this feature; with
the default the driver aborts at the first averaging step. That is inherited
driver-mct behaviour, not new here, but it is easy to trip over.

anlgce-ub22, gnu, SIA solver: ERS.ne30pg2_r05_IcoswISC30E3r5_gis4to40.MALISIA
and ERS_Ld5 IGELM_MLI both pass with bit-for-bit restart, and the MALISIA
coupler history is identical to driver-mct across all 35 fields.

Note on a pre-existing MOAB/Trilinos conflict

Running MALI with Albany under driver-moab currently segfaults on chrysalis
because MOAB ships a static Zoltan 3.901 whose symbols preempt the Zoltan 15.1.1
in Trilinos' libzoltan.so that Albany was compiled against. This predates this
branch and is unrelated to it — the crash happens in MALI's initial velocity
solve, before any coupling exchange, and needs only MOAB and Albany in the same
executable. It is reported as #8660 and will be fixed by rebuilding MOAB against
Trilinos' Zoltan. The Albany validation above was done with a local
-Wl,--exclude-libs,libzoltan.a workaround that is deliberately not included
in this pull request.

Not included

Ocean-glc thermal forcing and ice shelf coupling, glc-to-ocean and glc-to-ice
runoff mapping, and x2gacc restart are deferred to a follow-up. The prep_ocn
and prep_ice glc mapper inits still need their moab context there. The existing
MCT code for these is left in place and inert, and the tests covering that scope
(tests.py:154-155, :353) remain pinned to Vmct.

Bit-for-bit status

Expected to be BFB for driver-mct: the only shared code touched is
glc_comp_mct.F (an import/export refactor to raw arrays, with the moab path
behind HAVE_MOAB) and mpas_moabmesh.F (planar mesh support, which only runs
when a MOAB mesh instance is created). The Vmct tests above all pass. This has
not yet been confirmed against master baselines, which should be done before
merge.

[bishtgautam]

@rljacob, This isn't touching anything in ELM for me to review.

[rljacob]

Fun with MOAB h5m files. The partition of Greenland in MALI on 16 procs.

[rljacob]

This was needed to build on GCE and is the subject of a separate PR #8655

[rljacob]

I'll wait for #8655 to be resolved and then rebase this or just remove this change.

[rljacob]

Reviewers who aren't coupler developers should just check the changes in glc_comp_mct.F90 and maybe test that other cases currently run with driver-mct are not broken (I'll run the land-ice-developer suite on this).

Also see the comments from Claude about slight difference in Surface speed between the mct and moab cases. Results are BFB identical between the two cases otherwise.

[matthewhoffman]

I inspected the change to glc_comp_mct.F and, other than the one item I was uncertain about, I didn't see anything that concerned me - all of the MOAB functionality seems cleanly gated off.

Two other items:

1. As for the need to support planar meshes, if I understand correctly, that's because MOAB no longer requires pre-built mapping files and needs a mesh object at runtime. How is this handling distortion from map projection? That was a complication that had to be handled carefully when creating mapping files.

2. The slight difference in surfaceSpeed is curious, but, as noted, does not seem to be related to the moab coupling. Still, I would be more comfortable if we knew where this originated in this case. Are the runs being compared run with a different number of processors? Albany is not BFB across processor count. Other than that, I would expect surfaceSpeed to be BFB, so if that isn't the explanation it would be good for use to identify where it's coming from. (cc: @mperego )

I'm approving now but I'm happy to re-review if the above items lead to further changes.

[rljacob]

./create_test e3sm_landice_developer --compiler gnu --driver mct --pesfile ../../cime_config/testmods_dirs/config_pes_tests.xml still passes on this branch on chrysalis.

[rljacob]

The slight difference in surfaceSpeed is curious, but, as noted, does not seem to be related to the moab coupling. Still, I would be more comfortable if we knew where this originated in this case. Are the runs being compared run with a different number of processors?

No they were on the same number of processors. This is what Claude thinks:

Albany produces a different answer with the same code, same partition, and same inputs in a MOAB-linked process. That points at something in the process environment the solve is sensitive to — remaining candidates are duplicate symbol resolution somewhere in the 112 extra exports that I haven't attributed (the HDF5 H5T_NATIVE_* globals are the only non-STL ones), or MOAB's Zoltan instance perturbing shared state despite being a separate copy.

We could also try linking with debug builds of Albany and MOAB (debug build of E3SM isn't enough).

[matthewhoffman]

Thanks, @rljacob . Those explanations seem plausible. @mperego , do you think those possible explanations are sufficient to accept this difference without worrying about it?

[mperego]

Yes, I would not worry about it. I also checked with @bartgol and he agrees. Thanks!