Skip to content

Conversation

@matthewhoffman
Copy link
Contributor

@matthewhoffman matthewhoffman commented Apr 4, 2025

This pull request extends the existing thermal-forcing coupling functionality from a single depth to a predefined set of z-levels. The implementation follows that of multiple elevation classes in the LND-GLC coupling. A new xml variable is added, GLC_NZOC. It has a default value of 0, meaning TF coupling is disabled. The standard value for when the coupling is activated is 30, which follows the TF z-levels used by ISMIP6. Other values are supported, and an additional example with 4 levels is included. A new coupler module, glc_zocnclass_mod.F90, is added which is follows very closely from glc_elevclass_mod.F90. It contains subroutines and functions to initialize the z-level classes and to access the z-level values. It also has subroutines to provide the z-level indices as strings, which the coupling field initialization code and the component import/export routines uses to create GLC_NZOC fields, one for each layer.

This 3d TF coupling is used for ice-sheet outlet glacier facemelting (melting at vertical glacier margins terminating in the ocean), and so is used for a process that we do not anticipate resolving explicitly in MPAS-Ocean (or Omega). As such, this coupling should be included anytime both OCN and GLC are active, and logic has been added so that this is the case based on compset definition. There are versions of the GLC_NZOC variable in both the MPAS-Ocean and MALI Registries so that it can be set by namelist settings controlled by CIME/build operations.

The MPAS-Ocean mpas_ocn_time_average_coupled.F module has been modified to loop over the z-levels and calculate TF at each one. Then, the time-averaged 3d TF values are passed to the coupler, where they are remapped to the GLC mesh and passed to MALI. MALI assigns the values to its existing 3d TF field, and the uses its bathymetry-aware extrapolation routine to extrapolate the TF values from the edge of the MPAS-Ocean domain to wherever the ice-sheet margin is on the MALI mesh. From there, the TF is used to force the facemelting parameterization. Mapping files have been updated to not include extrapolation in the mapping so as to allow the bathymetry-aware extrapolation in MALI to operate, where appropriate.

The facemelting flux in MALI is exported to the previously unused Fogg_rofi flux, which already gets combined in the coupler with ice runoff to be imported in MPAS-Ocean in the Foxx_rofi flux. Long term, we may want to separate these runoff fluxes so that the horizontal and vertical distribution of facemelting can be handled differently than other solid runoff, but this provides a reasonable first approximation.

Note that this 3d TF can also be used to force the ISMIP6 ice-shelf basal melt parameterization in MALI, which provides a less sophisticated alternative to using ice-shelf basal melt rates calculated in the coupler (or MPAS-Ocean). This simpler approach will be used for the few, small ice shelves in Greenland, which are not included in the MPAS-Ocean meshes (and would not be resolved even if they were). This approach can also be used for Antarctica as an alternative that avoids grounding line migration issues in MPAS-Ocean by taking advantage of the TF extrapolation in MALI. These different methods for ice-shelf basal melt fluxes are handled through a new XML variable OCN_GLC_ISMF_COUPLING, which specifies how ice-shelf melt fluxes are handled in both MPAS-Ocean and MALI.

The previous 2d thermal forcing coupling has been removed. Tests have been updated to test this new 3d TF capability for grids with both GIS and AIS.

[NML]
[non-BFB] only for configurations with active MALI and ocn

@matthewhoffman matthewhoffman added MPAS-ocean Concerning the MPAS-ocean model coupled to E3SM. MPAS-albany-landice Concerning the MPAS-Albany land ice model Coupler Related to code in driver-mct or driver-moab or component connections to the coupler. Stealth PR has feature which, if turned on, could change climate. fka FCC labels Apr 4, 2025
@matthewhoffman matthewhoffman requested review from jonbob and xylar April 4, 2025 20:17
@matthewhoffman
Copy link
Contributor Author

Discussion and testing documented at E3SM-Ocean-Discussion#121

Copy link
Contributor

@xylar xylar left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Approving based on a careful inspection of the code (with a particular focus on the MPAS-Ocean changes) and my own testing, including adding debug fields that reconstructed the vertical coordinate in addition to the thermal forcing on the z-levels that glc accepts. Also approving based on other testing reported in E3SM-Ocean-Discussion#121

@xylar
Copy link
Contributor

xylar commented Apr 7, 2025

Assuming #7195 goes in before this, there are a few things we'll need to fix after rebasing on to that merge:

  • fix one more ocn_ismf
  • update the mapping files between ocn and glc for SORRMr4 to following this conversion.

@jonbob jonbob requested a review from rljacob April 8, 2025 14:45
@jonbob
Copy link
Contributor

jonbob commented Apr 8, 2025

@rljacob -- we've asked you to review this as well since it touches the coupler

@matthewhoffman
Copy link
Contributor Author

@rljacob , I'd be happy to meet and walk you through it if that would be easier for you.

@matthewhoffman
Copy link
Contributor Author

matthewhoffman commented Apr 10, 2025

FYI, I plan to rebase this to deal with known conflicts with #7195, now that 7195 is on master. I hope to do that tomorrow.

@xylar xylar self-requested a review April 10, 2025 15:42
@xylar
Copy link
Contributor

xylar commented Apr 10, 2025

I'm running a GMPAS-JRA1p5-DIB-PISMF-SIS.TL319_oQU240wLI_ais8to30 case and seeing:

 84:  0 0x0000000000012cf0 __funlockfile()  :0
 84:  1 0x0000000001e787bb glc_comp_mct_mp_init_ocean_z_levels_.A()  /lcrc/group/e3sm/ac.xylar/scratch/chrys/20250410.GMPAS-JRA1p5-DIB-PISMF-SIS.TL319_oQU240wLI_ais8to30.chrysalis/mpas-albany-landice/driver/glc_comp_mct.f90:1829
 84:  2 0x0000000001e7c376 glc_comp_mct_mp_glc_init_mct_()  /lcrc/group/e3sm/ac.xylar/scratch/chrys/20250410.GMPAS-JRA1p5-DIB-PISMF-SIS.TL319_oQU240wLI_ais8to30.chrysalis/mpas-albany-landice/driver/glc_comp_mct.f90:706
 84:  3 0x0000000000452c04 component_mod_mp_component_init_cc_()  /gpfs/fs1/home/ac.xylar/e3sm_work/E3SM/matthewhoffman/ocn-glc/3d-tf-ocn-glc/driver-mct/main/component_mod.F90:257
 84:  4 0x000000000043f1dc cime_comp_mod_mp_cime_init_()  /gpfs/fs1/home/ac.xylar/e3sm_work/E3SM/matthewhoffman/ocn-glc/3d-tf-ocn-glc/driver-mct/main/cime_comp_mod.F90:1531
 84:  5 0x000000000044f89a MAIN__()  /gpfs/fs1/home/ac.xylar/e3sm_work/E3SM/matthewhoffman/ocn-glc/3d-tf-ocn-glc/driver-mct/main/cime_driver.F90:122
 84:  6 0x000000000041b9a2 main()  ???:0
 84:  7 0x000000000003ad85 __libc_start_main()  ???:0
 84:  8 0x000000000041b8ae _start()  ???:0
 84: =================================

I'm looking into this but wanted to give you a heads up that a fix is likely needed.

@xylar
Copy link
Contributor

xylar commented Apr 10, 2025

@matthewhoffman, I think maybe you'd better take this one on. It seems like maybe the ismip6Melt package isn't active but the code things there are 30 vertical z-levels and is trying to print them out even though the fields aren't allocated.

https://github.com/E3SM-Project/E3SM/blob/matthewhoffman/ocn-glc/3d-tf-ocn-glc/components/mpas-albany-landice/src/mode_forward/mpas_li_core_interface.F#L170-L179

https://github.com/E3SM-Project/E3SM/blob/matthewhoffman/ocn-glc/3d-tf-ocn-glc/components/mpas-albany-landice/driver/glc_comp_mct.F#L1832-L1834

@xylar
Copy link
Contributor

xylar commented Apr 10, 2025

It seems like there will be 30 vertical levels if both MALI and MPAS-Ocean are active:
https://github.com/E3SM-Project/E3SM/blob/matthewhoffman/ocn-glc/3d-tf-ocn-glc/driver-mct/cime_config/config_component_e3sm.xml#L837
But that the package is not always on for this situation. So it seems like another check is needed to make sure the package is also on before accessing the ismip6shelfMelt* variables.

@matthewhoffman
Copy link
Contributor Author

Thanks, @xylar , that makes sense. I'll address that today.

@matthewhoffman matthewhoffman force-pushed the matthewhoffman/ocn-glc/3d-tf-ocn-glc branch from 1e3b2e1 to 10d062f Compare April 12, 2025 03:42
@matthewhoffman
Copy link
Contributor Author

@xylar , I've rebased and made some changes:

  • dealt with conflicts or issues needing updating during the rebase
  • updated MALI grid files to fix the dimension issue we diagnosed (and other minor update Trevor found)
  • fixed an out-of-bounds error that the other fixes exposed

I also touched base with Jon about needing some new mapping files for SOwISC12to30E3r4 to use TF coupling, now that that came in as part of the rebase.

I confirmed the following tests pass after the changes that could affect these grids:

  • SMS_D_Ld5.TL319_oQU240wLI_gis4to40.MPAS_FOLISIO_JRA1p5.chrysalis_gnu
  • SMS_D_Ld5.TL319_oQU240wLI_ais8to30.MPAS_FOLISIO_JRA1p5.chrysalis_gnu
  • SMS_D_Ld5.TL319_IcoswISC30E3r5_ais4to20.MPAS_FOLISIO_JRA1p5.chrysalis_gnu

However, SMS_D_Ld5.TL319_oQU240wLI_ais8to30.GMPAS-JRA1p5-DIB-PISMF-SIS.chrysalis_gnu still fails. The errors Xylar was seeing are now gone and the code gets through init for all components. However, it dies with this error in the e3sm log file:

 55: MCT::m_AttrVect::indexRA_:: FATAL--attribute not found: "norm8wt" Traceback:
 55: |X|MCT::m_AttrVect::indexRA_
 55: 037.MCT(MPEU)::die.: from MCT::m_AttrVect::indexRA_()

I will need to look into this further, but am stopping here for now.

logger.warning("WARNING: The specified compset is requesting ocean ICs spunup from a G-case")
logger.warning(" But no file available for this grid.")
if ocn_ismf == 'data':
if ocn_glc_ismf_coupling == 'data_mpaso':
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great, glad you caught this one, too. I needed this change in my test-merge branch.

@matthewhoffman
Copy link
Contributor Author

@xylar , I've pushed commits that will handle the situation where both OCN and GLC are active but we don't want TF coupling (e.g. MALI DATA and STATIC modes). With these changes, SMS_D_Ld5.TL319_oQU240wLI_ais8to30.GMPAS-JRA1p5-DIB-PISMF-SIS.chrysalis_gnu now runs without error.

@matthewhoffman
Copy link
Contributor Author

...and the compsets that do include TF coupling are now broken. I'll wait to update again until I can get all test cases working again.

@matthewhoffman
Copy link
Contributor Author

I fixed a bug introduced in the previous commit and reran these tests to confirm things are working as intended:

  • SMS_D_Ld5.TL319_oQU240wLI_gis4to40.MPAS_FOLISIO_JRA1p5.chrysalis_gnu PASS
  • SMS_D_Ld5.TL319_oQU240wLI_ais8to30.MPAS_FOLISIO_JRA1p5.chrysalis_gnu PASS
  • SMS_D_Ld5.TL319_oQU240wLI_ais8to30.GMPAS-JRA1p5-DIB-PISMF-SIS.chrysalis_gnu PASS
  • SMS_D_Ld5.TL319_IcoswISC30E3r5_ais4to20.MPAS_FOLISIO_JRA1p5.chrysalis_gnu PASS

Other than the needed mapping files for SOwISC12to30E3r4, everything is resolved that I am aware of.

@matthewhoffman
Copy link
Contributor Author

@jonbob , I added a commit that updates the SOwISC12to30E3r4 mapping files. I successfully ran a test with SMS_D_Ld5.TL319_SOwISC12to30E3r4_ais4to20.MPAS_FOLISIO_JRA1p5.chrysalis_gnu to ensure that they don't throw any errors.

There are no remaining revisions to make that I am aware of. @rljacob , this is at a good point for you to review if you would like to review it.

"ERS_Ld5_D.T62_oQU240.GMPAS-IAF.mpaso-conservation_check",
"ERS_Ld5_PS.ne30pg2_r05_IcoswISC30E3r5.CRYO1850-DISMF.mpaso-scaled_dib_dismf",
"ERS_Ld5.TL319_oQU240wLI_gis20.MPAS_LISIO_JRA1p5.mpaso-ocn_glc_tf_coupling",
# OCN/GLC 3d TF coupling GIS test:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should we run this e3sm_ocnice_stealth_features suite somewhere on a regular basis?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rljacob -- it runs every night on gcp12

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh I see. Its included in extra coverage.

<map name="OCN2GLC_SHELF_FMAPNAME">cpl/gridmaps/oQU240wLI/map_oQU240wLI_to_gis20km_esmfaave.20240919.nc</map>
<map name="OCN2GLC_SHELF_SMAPNAME">cpl/gridmaps/oQU240wLI/map_oQU240wLI_to_gis20km_esmfbilin.20240919.nc</map>
<map name="OCN2GLC_TF_SMAPNAME">cpl/gridmaps/oQU240wLI/map_oQU240wLI_to_gis20km_deeperThan300m.esmfneareststod.20240919.nc</map>
<map name="OCN2GLC_TF_SMAPNAME">cpl/gridmaps/oQU240wLI/map_oQU240wLI_to_gis20km_esmfbilin.20240919.nc</map>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did tempestremap not work?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still can't get it to work when both grids are non-global like this. I understand that mbtempest can do that, but I haven't really entrained that into my workflow

xylar and others added 14 commits May 23, 2025 23:53
'none' is set when the compset does not include the ocean
Both Greenland meshes currently have issues that cause problems when
using OCN2GLC_TF_SMAPNAME mapping files that don't have extrapolation.
For the gis4to40 mesh, there is an "inland sea" region that does not get
valid ocean thermal forcing, leading to MALI errors in facemelting.
Here, I've switched the map file to one that includes extrapolation to
work around this issue.  The inland sea problem is being fixed in a
different PR, and after that is merged, this map file can be switched
back to the version without extrapolation, so that the bathymetry-aware
extrapolation inside of MALI can be used instead.

The gis1to10kmR2 mesh has a different problem - the ocean "gutter" is
too narrow and only overlaps the IcoswISC30E3r5 grid in a few places.
Switching to a OCN2GLC_TF_SMAPNAME mapping file that includes
extrapolation similarly works around this problem.  A more thorough
solution will require creating a new high-res GIS mesh with a larger
gutter.
We want 3d TF coupling on for all compsets with both MPAS-Ocean and MALI
active because we want to use 3d TF for facemelting even if we use
ice-shelf basal melt fluxes from the coupler.  However, we do not want
facemelting when MALI is DATA or STATIC modes, so the TF coupling should
be inactive in those cases.
1. OCN_ISMF -> OCN_GLC_ISMF_COUPLING for landice freshwater flux logic
2. OCN_ISMF -> OCN_GLC_ISMF_COUPLING for sowisc12to30e3r4 mesh logic
* mpas.ais8to30km: removed nISMIP6OceanLayers dimension and assoc. vars
* mpas.ais4to20km: removed nISMIP6OceanLayers dimension and assoc. vars
* mpas.gis4to40km: added missing muFriction field
This distinction is needed for compsets where both OCN and GLC are
active but we do not want TF coupling (e.g. MALI STATIC and DATA modes).
@matthewhoffman matthewhoffman force-pushed the matthewhoffman/ocn-glc/3d-tf-ocn-glc branch from 4361fa2 to b14f9bf Compare May 24, 2025 05:56
@matthewhoffman
Copy link
Contributor Author

Force-pushed after rebase to deal with known conflicts from #7210 being merged to master.

@matthewhoffman
Copy link
Contributor Author

@jonbob , I've rebased this PR and added the needed commit to introduce the time-averaging for the melt fluxes to the coupler. After those updates, I retested with:

./create_test --wait --walltime 0:30:00  ERS_D_Ld5.TL319_oQU240wLI_ais8to30.MPAS_FOLISIO_JRA1p5.chrysalis_gnu

and the test passed.

This is ready for final review and merge.

@jonbob jonbob added the NML label May 27, 2025
"ERS_Ld5_PS.ne30pg2_r05_IcoswISC30E3r5.CRYO1850-DISMF.mpaso-scaled_dib_dismf",
"ERS_Ld5.TL319_oQU240wLI_gis20.MPAS_LISIO_JRA1p5.mpaso-ocn_glc_tf_coupling",
# OCN/GLC 3d TF coupling GIS test:
"ERS_Ld5.TL319_oQU240wLI_gis4to40.MPAS_FOLISIO_JRA1p5.mpaso-jra_1958",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@matthewhoffman -- it doesn't look like this stealth test is actually turning the feature on?

@jonbob jonbob added non-BFB PR makes roundoff changes to answers. and removed Stealth PR has feature which, if turned on, could change climate. fka FCC labels May 28, 2025
jonbob added a commit that referenced this pull request May 28, 2025
Add 3d thermal forcing to OCN-GLC coupling

This pull request extends the existing thermal-forcing coupling
functionality from a single depth to a predefined set of z-levels. The
implementation follows that of multiple elevation classes in the LND-GLC
coupling. A new xml variable is added, GLC_NZOC. It has a default value
of 0, meaning TF coupling is disabled. The standard value for when the
coupling is activated is 30, which follows the TF z-levels used by
ISMIP6. Other values are supported, and an additional example with 4
levels is included. A new coupler module, glc_zocnclass_mod.F90, is
added which is follows very closely from glc_elevclass_mod.F90. It
contains subroutines and functions to initialize the z-level classes and
to access the z-level values. It also has subroutines to provide the
z-level indices as strings, which the coupling field initialization code
and the component import/export routines uses to create GLC_NZOC fields,
one for each layer.

This 3d TF coupling is used for ice-sheet outlet glacier facemelting
(melting at vertical glacier margins terminating in the ocean), and so
is used for a process that we do not anticipate resolving explicitly in
MPAS-Ocean (or Omega). As such, this coupling should be included anytime
both OCN and GLC are active, and logic has been added so that this is
the case based on compset definition. There are versions of the GLC_NZOC
variable in both the MPAS-Ocean and MALI Registries so that it can be
set by namelist settings controlled by CIME/build operations.

The MPAS-Ocean mpas_ocn_time_average_coupled.F module has been modified
to loop over the z-levels and calculate TF at each one. Then, the time-
averaged 3d TF values are passed to the coupler, where they are remapped
to the GLC mesh and passed to MALI. MALI assigns the values to its
existing 3d TF field, and the uses its bathymetry-aware extrapolation
routine to extrapolate the TF values from the edge of the MPAS-Ocean
domain to wherever the ice-sheet margin is on the MALI mesh. From there,
the TF is used to force the facemelting parameterization. Mapping files
have been updated to not include extrapolation in the mapping so as to
allow the bathymetry-aware extrapolation in MALI to operate, where
appropriate.

The facemelting flux in MALI is exported to the previously unused
Fogg_rofi flux, which already gets combined in the coupler with ice
runoff to be imported in MPAS-Ocean in the Foxx_rofi flux. Long term, we
may want to separate these runoff fluxes so that the horizontal and
vertical distribution of facemelting can be handled differently than
other solid runoff, but this provides a reasonable first approximation.

Note that this 3d TF can also be used to force the ISMIP6 ice-shelf
basal melt parameterization in MALI, which provides a less sophisticated
alternative to using ice-shelf basal melt rates calculated in the
coupler (or MPAS-Ocean). This simpler approach will be used for the few,
small ice shelves in Greenland, which are not included in the MPAS-Ocean
meshes (and would not be resolved even if they were). This approach can
also be used for Antarctica as an alternative that avoids grounding line
migration issues in MPAS-Ocean by taking advantage of the TF
extrapolation in MALI. These different methods for ice-shelf basal melt
fluxes are handled through a new XML variable OCN_GLC_ISMF_COUPLING,
which specifies how ice-shelf melt fluxes are handled in both MPAS-Ocean
and MALI.

The previous 2d thermal forcing coupling has been removed. Tests have
been updated to test this new 3d TF capability for grids with both GIS
and AIS.

[NML]
[non-BFB] only for configurations with active MALI and ocn
@jonbob
Copy link
Contributor

jonbob commented May 28, 2025

Testing shows expected results:

  • NML DIFF only for SMS_D_Ld1.ne30pg2_r05_IcoswISC30E3r5.WCYCL1850.chrysalis_intel.allactive-wcprod
  • NML and baseline DIFFs for e3sm_landice_developer using chrysalis_gnu

Merged to next

@jonbob jonbob merged commit db098a3 into master May 29, 2025
7 checks passed
@jonbob jonbob deleted the matthewhoffman/ocn-glc/3d-tf-ocn-glc branch May 29, 2025 16:25
@xylar
Copy link
Contributor

xylar commented May 29, 2025

Thanks @jonbob!! And thanks @matthewhoffman for the hard work on this!

@jonbob
Copy link
Contributor

jonbob commented May 29, 2025

merged to master and expected baseline and NML DIFFs blessed

@xylar xylar mentioned this pull request Jul 4, 2025
54 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Coupler Related to code in driver-mct or driver-moab or component connections to the coupler. MPAS-albany-landice Concerning the MPAS-Albany land ice model MPAS-ocean Concerning the MPAS-ocean model coupled to E3SM. NML non-BFB PR makes roundoff changes to answers.

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants