OCN to GLC thermal forcing coupling

[matthewhoffman]

This pull request introduces coupling from OCN to GLC to pass thermal forcing at a prescribed ocean depth (300 m) from MPAS-Ocean to MALI, where MALI uses it to calculate grounded marine melting through an existing 'facemelting' parameterization. Facemelting as a function of ocean thermal state is a critical OCN/GLC coupling for Greenland, but it is relevant to Antarctica as well. This OCN-GLC thermal forcing coupling is activated whenever OCN is present and GLC is active.

This PR has pieces in MPAS-Ocean, the coupler, and MALI:

• In MPAS-Ocean, a new avgThermalForcingAtCritDepth field is added and calculated. The depth at which to calculate the thermal forcing is namelist configurable. The default is 300 m.
• In MALI, the field passed from the coupler is collected in ismip6_2dThermalForcing and used in the existing facemelting parameterization. Namelist and streams are updated to use facemelting in all simulations and output relevant variables. For compsets where TF is not being coupled, it will have values of 0 and result in no facemelting being applied; it is safe to have this option generally turned on.
• In the coupler, a new remapping object is created that handles the TF mapping independently of any other fields passed from OCN to GLC. New mapping files are added. Any coupling variables related to the existing ice-shelf coupling now include the 'shelf' suffix and all coupling variables related to this new thermal-forcing coupling have the 'tf' suffix.
• A new compset is created that is a standard G-case but with active MALI. New mapping files are added between the relevant grids for the new mapper.
• The new OCN2GLC_TF_SMAPNAME mapping should be a nearest source to destination mapping with the MPAS-Ocean source mesh masked to only include grid cells that are deeper than the critical depth to avoid extrapolating undefined values. A tool to generate that map file is here: Add script for adding map of sufficiently deep ocean bathymetry to a scrip file MPAS-Dev/MPAS-Tools#578
• Two new model_grids are added that are compatible with this feature and include the special mapping file required: TL319_IcoswISC30E3r5_gis1to10kmR2 and TL319_oQU240wLI_gis20
• The new coupling is controlled by a new namelist option in MPAS-Ocean: config_glc_thermal_forcing_coupling_mode. It defaults to 'off' and there are no compsets that explicitly set it for now, making this a stealth feature.

[matthewhoffman]

@cbegeman & @jonbob , I tagged you both on this PR, though we may want to bring in some other people to review/contribute, as well. The PR functions so in some sense is complete, but we may want to discuss some of the implementation choices and revise. I opened the PR now so you have some concrete to look at as we discuss.

[cbegeman]

The 300m thermal forcing field seems right on the ocean side. Given my limited coupler development experience, I think this looks good. Would it be helpful for me to look over some of the cpl and output files from your recent test? Or run this test on another machine TL319_IcoswISC30E3r5_gis20?

[matthewhoffman]

Testing

Set up and ran case with the following. Note that a custom mapping file is required that is not yet on the inputdata server.

export CASE_ROOT=/home/$USER/e3sm-gis/e3sm_cases
export COMPSET=GMPAS-JRA1p5-MALI
export GRID=TL319_IcoswISC30E3r5_gis20
export COMPILER=gnu
export E3SM_CASE=`date +"%Y%m%d"`.${COMPSET}.${GRID}.GIS-ocn-TFcpl-test.${MACHINE}

cd $E3SM_ROOT/cime/scripts
./create_newcase \
-case $CASE_ROOT/$E3SM_CASE \
-mach $MACHINE \
-compset $COMPSET \
-res $GRID \
--compiler $COMPILER

cd $CASE_ROOT/$E3SM_CASE
./case.setup
./xmlchange OCN2GLC_TF_SMAPNAME="/home/ac.mhoffman/e3sm-gis/mapping_files/map_IcoswISC30E3r5_to_gis20km_esmfneareststod.20240422.deeperThan300m.nc"
./case.build

./xmlchange STOP_OPTION=nmonths
./xmlchange STOP_N=1
./xmlchange REST_OPTION=nmonths
./xmlchange REST_N=1
./xmlchange JOB_QUEUE=compute
./xmlchange JOB_WALLCLOCK_TIME="01:00:00"

./case.submit

This panel on the left shows the thermal forcing on the MALI mesh after nearest neighbor remapping from MPAS-Ocean. IThe panel on the right shows the time-series of thermal forcing for an arbitrary grid point, showing that it is indeed changing in time (x-axis units are days, which is the coupling interval).

This figure shows the resulting facemelting thickness calculated in MALI based off of the thermal forcing from MPAS-Ocean. Blue are ocean cells and white contour is ice-sheet margin. You can see that we are only getting facemelting in a small number of cells. This is because the 20 km mesh is too coarse to resolve fjords and there are not a lot of cells where the ice-sheet margin terminates against the ocean.

[matthewhoffman]

@cbegeman , thank for your initial review. I went ahead and made the depth at which to calculate thermal forcing a namelist option and added information to the option description that a custom mapping file needs to be generated to be consistent with that depth. I also did some additional cleanup to better differentiate the coupling variables for this TF coupling from the existing ice-shelf coupling. If you want to try running this on Perlmutter, that would be good. We are restricted to machines where Albany exists, so Perlmutter (with gnu) would be the best place to try. You will also need the custom mapping file that is currently in my home directory on Chrysalis. Once we are happy with the branch, I'll coordinate with Jon about the best way to name it.

[matthewhoffman]

@jonbob and @xylar , everything is working as I would like and this is ready for review. Once the two of you and Carolyn are happy with it, I can open a formal PR on the main E3SM repo and ask Rob to review it. I'm leaving for two weeks of vacation at the end of next week, so if we're able to get it to that stage before I leave, that would be great, but it's ok if not.

It would be good to get #95 merged so we can use it to confirm that this PR doesn't break any of the shelf coupling. (Or I can test that with a temporary merge of the two branches)

[cbegeman]

@matthewhoffman Great! I can test it on perlmutter. I'll get in touch with you to get that path when I'm ready to set it up (in the next day or two).

[matthewhoffman]

@cbegeman , sounds good! (That path is buried in my notes above for setting up the case: /home/ac.mhoffman/e3sm-gis/mapping_files/map_IcoswISC30E3r5_to_gis20km_esmfneareststod.20240422.deeperThan300m.nc)

[matthewhoffman]

Testing with Icos ocean mesh coupled to 1km GIS mesh yields expected results. Here is a zoom-in on Humboldt Glacier, the widest marine-terminating glacier in Greenland:

and here is the northwest region of Greenland, which contains many small marine outlet glaciers:

[cbegeman]

These changes look good to me! I ran the test on perlmutter cpu at /pscratch/sd/c/cbegeman/e3sm_scratch/pm-cpu/20240529.GMPAS-JRA1p5-MALI.TL319_IcoswISC30E3r5_gis20.GIS-ocn-TFcpl-test.pm-cpu and verified that it ran and produced non-zero faceMeltingThickness.

[xylar]

A few comments for now.

[matthewhoffman]

@xylar and @jonbob , I've addressed the two comments from Xylar's first pass of this PR, and it's ready for further review when you both have time. I'd be happy to meet and walk you through the implementation if that would help.

[xylar]

@matthewhoffman, I think this looks good and is ready to go to E3SM once @jonbob agrees.

The one question that came to mind when I looked through the code was whether we want a flag and/or package for these changes. Right now as it's written, we'll be computing the thermal forcing at the given depth in all E3SM simulations with active ocean regardless of whether we're coupling to MALI or not. I think that's okay but wanted to make sure that's the approach we want to take. I don't believe the TF will be written out to any stream by default so the information is just being lost if it doesn't go to MALI. There's some complexity in figuring out if we're coupling to glc in order to set a namelist option to turn this feature on or off, so we may not want to go down that road. I just want us to have weighed the pros and cons.

[matthewhoffman]

@xylar , that's a good point. We could so something analogous to the ice-shelf basal melting in the coupler, which is activated when config_land_ice_flux_mode = 'coupled' in the ocean namelist. We would eventually want a compset that enables this feature, but for now we could put it in a testmod usernl like we do here: https://github.com/E3SM-Project/E3SM/blob/master/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glcshelf/user_nl_mpaso

What do you think? That might be better than leaving it fully exposed as it is now.

Looking closer at what I had done, the coupler operations for the TF coupling are only activated when ocean is present and glc is prognostic:

@@ -1768,6 +1771,7 @@ subroutine cime_init()
    if (ocn_present) then
       if (atm_prognostic) ocn_c2_atm = .true.
       if (atm_present   ) ocn_c2_atm = .true. ! needed for aoflux calc if aoflux=atm
       if (glc_prognostic) ocn_c2_glctf = .true.

So it would be nice to have the ocean calculations also only occur when these conditions are met. I'm not sure offhand if it's possible for the information to flow that direction. @jonbob might know or I will take a closer look.

[xylar]

Yes, I think that would be my preference. For now, this would be a stealth feature only activated manually with a user namelist option.

[matthewhoffman]

Ok, yes, let's add that grid. I'm sure it'll be useful in other contexts as we get these features running more regularly so might as well take care of it now. Let me know if I can do anything to help.

[matthewhoffman]

PR moved to main E3SM repo at: E3SM-Project#6632