Add 3d thermal forcing to OCN-GLC coupling

[matthewhoffman]

This pull request extends the existing thermal forcing functionality from a single depth to a predefined set of z-levels. The coupler implementation follows that of multiple elevation classes in the LND-GLC coupling. A new coupler 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 ocean 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.

The intention is that this 3d TF coupling will be used for 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. 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 the 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.

The facemelting flux in MALI is exported to the previously unused Fogg_rofl flux, which already gets combined in the coupler with river runoff to be imported in MPAS-Ocean in the Foxx_rofl 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 river 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 could 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.

This PR is now functional, but some additional work is needed before this is complete and ready for review at the main E3SM repo:

• new OCN2GLC_TF_SMAPNAME mapping files that use bilin method but exclude extrapolation need to be created for all the grids that currently support that mapping. We should wait until the gis4to40 mesh is merged and this branch rebased, so that these new mapping files can be added for that mesh too.
• a decision needs to be made whether to continue to support the previous 2d TF coupling in addition to this new 3d version. I don't see a use-case where the 2d version would ever be used over the 3d version, so I advocate for deprecating the 2d version.
• a mask field should be added to the coupling so that MALI knows what locations on its mesh were valid ocean locations in MPAS-Ocean. Right now this is handled in a hacky (but functional) way.
• ocean developers may want to modify the way 3d TF is being calculated in MPAS-Ocean
• there are some temporary commits to namelist and streams settings to ease testing that should be removed or adjusted.
• the existing 2d TF test (ERS_Ld5.TL319_oQU240wLI_gis20.MPAS_LISIO_JRA1p5.mpaso-ocn_glc_tf_coupling) should be replaced by a new test. We should wait until the gis4to40 mesh is merged and this branch rebased, so we can use that mesh for the updated test. The new test should be ERS_Ld5.TL319_oQU240wLI_gis4to40.MPAS_FOLISIO_JRA1p5.mpaso-jra_1958 and the ocn_glc_tf_coupling testmod should be removed from components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/
• a new AIS 3d TF test should be added: ERS_Ld5.TL319_oQU240wLI_ais8to30.MPAS_LISIO_JRA1p5.mpaso-jra_1958
• add a package for the 3d TF variable in MPAS-Ocean
• with this branch, TF-coupling for ice shelves (in addition to facemelting) is on for any compset that has both OCN and GLC active. Add the ability to use the ice-shelf melt rates in the coupler, so as not to set that development back. This might be tricky because it requires namelist settings changed in both MPAS-Ocean and MALI. The MPASO_ISMF should allow this to be handled fairly smoothly. cases with MPASO%.*PISMF will activate ISMF in the coupler. cases without a '*ISMF" modifier for MPASO will have no ISMF calculated in the coupler or MPAS-Ocean, so 3d tf shouldn't interfere with anything. However, nl mods in MALI will be needed. Therefore the ocn_glcshelf testmod should be moved from MPAS-O to MALI testmods.

Required mapping files

OCN2GLC_TF_SMAPNAME: bilinear without extrapolation

• oQU240wLI_ais8to30
• oQU240wLI_gis4to40
• IcoswISC30E3r5_ais8to30
• IcoswISC30E3r5_ais4to20
• IcoswISC30E3r5_gis1to10kmR2
• IcoswISC30E3r5_gis4to40

GLC2OCN_LIQ_RMAPNAME: nearest-neighbor from rof2ocn masked

• oQU240wLI_ais8to30
• oQU240wLI_gis4to40
• IcoswISC30E3r5_ais8to30
• IcoswISC30E3r5_ais4to20
• IcoswISC30E3r5_gis1to10kmR2
• IcoswISC30E3r5_gis4to40

GLC2OCN_ICE_RMAPNAME: nearest-neighbor from rof2ocn no-mask

• oQU240wLI_ais8to30
• oQU240wLI_gis4to40
• IcoswISC30E3r5_ais8to30
• IcoswISC30E3r5_ais4to20
• IcoswISC30E3r5_gis1to10kmR2
• IcoswISC30E3r5_gis4to40

[matthewhoffman]

Testing with:

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

Note this is smoketest with debug on for testing purposes. The ais8to30 mesh is chosen as a low-res testing mesh over the gis20 mesh because the gis20 mesh has very little ocean mesh overlap and very few marine-terminating ice-sheet margins. The ais8to30 mesh also allows evaluation of how the TF coupling works for ice-shelf basal melting.

Results look like things are functioning correctly:

TF on the MALI mesh at two different depths:
Surface:

Deeper:

TF cross-section:

Ice-shelf basal melt rate (m/yr) on the MALI mesh calculated from ISMIP6 parameterization inside MALI:
(calculated as floatingBasalMassBalApplied/-910*3.15e7)

[matthewhoffman]

@cbegeman , @xylar , I marked you both as reviewers. I think a review from one of you is sufficient, but I wanted this on both of your radars.

Other people who may want to look at this or weigh in on details (but I don't think a review is required): @irenavankova, @darincomeau , @alexolinhager, @trhille

[matthewhoffman]

@jonbob , see note above about new OCN2GLC_TF_SMAPNAME mapping files. Let me know if you want me to create those or if you prefer to do so.

[matthewhoffman]

GIS Testing

./create_test --wait --walltime 0:30:00  SMS_D_Ld5.TL319_IcoswISC30E3r5_gis1to10kmR2.MPAS_FOLISIO_JRA1p5.chrysalis_intel

Note: need to modify env_mach_pes in the test directory and set all -4 values to -8 in the <entry id="NTASKS"> section to get this case to run. (Followed by ./case.setup --reset; ./case.build; ./case.submit )

Results look wonky, presumably because there is insufficient overlap between the OCN and GLC meshes. The gis1to10kmR2 MALI mesh has a very narrow gutter.

The 20km GIS mesh is also not useful for testing because the coarse resolution fails to resolve most marine outlet glaciers. Proper GIS testing will require a medium to high res mesh with a larger gutter (like the 4km mesh currently being tested).

[matthewhoffman]

After talking with @darincomeau , it became clear that there already is an unused coupler variable for ice-sheet runoff flux to the ocean, Fogg_rofl. This gets combined with river runoff to be received by MPAS-Ocean in the Foxx_rofl variable, which is assigned to the riverRunoff variable. While in the long term, we might want to handle river and facemelt runoff fluxes differently in terms of horizontal smoothing and vertical profile, it was easy to attach facemelting to this existing variable and get the FWF passed to MPAS-Ocean. I've added a commit that does that. It includes logic to also include the ice-shelf basal melt flux in Fogg_rofl if MALI is calculating that flux itself using 3d TF in the ISMIP6 parameterization (as opposed to ISMF being calculate by the coupler or MPAS-Ocean).

Here is a figure showing that the facemeling flux in MALI is being received by MPAS-Ocean properly. You have to zoom far in to see the nonzero facemelt flux values on the 8km MALI mesh, so they are highlighted with red circles.

This run dies in sea ice, so more work needs to go into confirming the flux magnitude is reasonable (for reference, the value at the mouth of the Amazon is 2e-4) and if some sort of smoothing needs to be added (which might requires separating this flux from Foxx_rofl).

[matthewhoffman]

@cbegeman and @xylar , I tried changing the melt flux from MALI to MPAS-Ocean to occur through rofi instead of rofl. I ran into a problem in that ice-shelf basal melt flux from the Jourdain parameterization can sometimes be freeze-on. When that happens I hit this error in the MPAS-Ocean driver:
https://github.com/E3SM-Project/E3SM/blob/master/components/mpas-ocean/driver/ocn_comp_mct.F#L2363-L2365

           if(iceRunoffFlux(n) < 0.0_RKIND) then
               call shr_sys_abort ('Error: incoming rofi_F is negative')
           end if

For the intermediate step of this PR, do you think it is better to switch back to rofl and ignore the heat flux for now, or disable that check in the driver (I'm not sure how dangerous that would be)?

Also, are you sure that rofi / iceRunoffFlux has its corresponding heat flux automatically accounted for? I'm not seeing where that happens.

[cbegeman]

Also, are you sure that rofi / iceRunoffFlux has its corresponding heat flux automatically accounted for? I'm not seeing where that happens.

Yes, see here:
https://github.com/E3SM-Project/E3SM/blob/b03e6d10557bfe21ac7cff8fd9986e37f7f284fd/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F#L678-L683

[cbegeman]

For the intermediate step of this PR, do you think it is better to switch back to rofl and ignore the heat flux for now, or disable that check in the driver (I'm not sure how dangerous that would be)?

I'm fine with switching back to rofl for now, but I would want to come up with a better solution before we devote significant simulation time to this intermediate step. I don't think we would have any problem with disabling the check, it's more of a check for prescribed forcing fields, as I don't think there's a way to get negative ice runoff prognostically.

[xylar]

@matthewhoffman, I would feel best about using rofi and disabling the check. Doug Jacobsen put that in 10 years ago without explanation and I don't think it is necessary. It's okay to freeze and thereby warm up the ocean a little -- this seems less dangerous than the current problem we have where positive ice runoff can supercool the ocean, generating frazil.

[xylar]

@matthewhoffman, I have only partially looked through the code but wanted to flag a significant problem before I continue.

You are using a z-level approach to find the depth that each value in glcZLevels corresponds to. This is not really a safe assumption anywhere in MPAS-Ocean but it is particularly poor for the Antarctic application because the vertical coordinate is very far from z-level inside of and within about a 20-cell ring around ice-shelf cavities. Instead, we need to find zMid that is closest to a given entry in glcZLevels.

[xylar]

Just wanting to check that these are the only 3 supported values. It sounded in our discussion like the number was flexible.

[xylar]

I see, these are hard-coded in glc_zocnclass_init_default()

[matthewhoffman]

That's correct. I followed the implementation of land model multiple elevation classes. In principle any number of levels (and any array of level elevations) is supported, but in practice a small set of numbers of levels are hard-coded and each has a hard-coded array of level positions. This seems like a reasonable tradeoff between flexibility and complexity.

[xylar]

Same question as above: is this restricted to only these 3 values?

[matthewhoffman]

See my longer response to a related question below. I followed the implementation of the land multiple elevation classes, which hardcodes a small number of options for the number of levels, and each one corresponds to a preset array of levels. There is functionality (in land MECs and also my new ocean z-levels) to initialize to a custom array of levels, but there is not at present a way to read them in from file. That would need to be added in the future if we wanted to support that level of generality.

[xylar]

I don't think we want to make this change.

[matthewhoffman]

Thanks for flagging - there are a number of streams settings I've committed for convenience but I will need to remove once the draft version of this PR is complete. I've separated such changes in commits with obvious names.

[xylar]

Okay, I've looked through the code and what I've mentioned so far seem like my only concerns. I really appreciate what a heap of work this has been. It's really exciting!

[xylar]

I see, these are hard-coded in glc_zocnclass_init_default()

[xylar]

Any way to make this more general (at least in the longer term)?

[matthewhoffman]

Probably - and that is what I had originally imagined doing - but the way I implemented it is how the LND multiple elevation classes are handled: https://github.com/E3SM-Project/E3SM/blob/master/driver-mct/shr/glc_elevclass_mod.F90#L85-L109 and after thinking about it more, the likelihood of wanting to specify a custom set of z-levels seems pretty low.

That said, the land multiple elevation class module includes an override version of the init routine that supports passing a pre-defined array instead of using one of the hard-coded ones: https://github.com/E3SM-Project/E3SM/blob/master/driver-mct/shr/glc_elevclass_mod.F90#L33-L36
I've included that functionality in the new z-level module. The override functionality is not used anywhere, but if we ever wanted to add that, we would control that here: https://github.com/E3SM-Project/E3SM/blob/master/driver-mct/shr/seq_flds_mod.F90#L472 The more complicated part would be implementing a mechanism by which to input the desired list of z-levels (presumably in a case xml file).

[xylar]

A few more small things that I'm changing in the update I'm working on. Just want to point them out for posterity.

[xylar]

@matthewhoffman, I moved the 3d thermal forcing calculation to a subroutine of its own. I switched to finding the depth of the glc z levels using zMid from the ocean and linear interpolation rather than nearest point.

I made sure that the test you ran above also passed for me:

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

That's as far as I got so far. It would be really good if we could write out the time-averaged fields on both the MPAS-Ocean and the MALI sides and carefully look at them (before extrapolation) to make sure they make sense.

[xylar]

Sorry for the force-push. Bad practice. I had an uncommitted change that I didn't want to add in an extra fix-up commit but probably should have.

[xylar]

Assuming the remaining check boxes (mapping files and removing temp commits) get checked, I think this is ready to move to E3SM from my point of view.

[matthewhoffman]

@xylar , FYI, I removed the temporary changes to the MPAS-O streams file from the commit where they occurred and force-pushed.

I also made a list in the PR description of the mapping files that need to be updated. Once I've coordinated with Jon to get those in place, I'll move this PR over to the main repo for final review.

[jonbob]

@matthewhoffman -- the new files have been added to this branch, and staged on the inputdata repo. Please test as much as you think is necessary

[matthewhoffman]

Testing

Ultra low res ocean

• SMS_D_Ld5.TL319_oQU240wLI_ais8to30.MPAS_FOLISIO_JRA1p5.chrysalis_gnu: PASS
• SMS_D_Ld5.TL319_oQU240wLI_gis4to40.MPAS_FOLISIO_JRA1p5.chrysalis_gnu: FAIL - appears to be related to 'inland seas' extrapolation issue fixed in this Add temperature and salinity to ocean extrapolation MALI-Dev/E3SM#124

Low res ocean

• SMS_D_Ld5.TL319_IcoswISC30E3r5_ais8to30.MPAS_FOLISIO_JRA1p5.chrysalis_gnu: PASS
• SMS_D_Ld5.TL319_IcoswISC30E3r5_ais4to20.MPAS_FOLISIO_JRA1p5.chrysalis_gnu: PASS
• SMS_D_Ld5.TL319_IcoswISC30E3r5_gis1to10kmR2.MPAS_FOLISIO_JRA1p5.chrysalis_gnu: FAIL - gis1to10 gutter is too narrow to receive sufficient ocean data for extrapolation to work. We could switch mapping back to nearest neighbor with some loss of accuracy due to missing bathymetry evaluation
• SMS_D_Ld5.TL319_IcoswISC30E3r5_gis4to40.MPAS_FOLISIO_JRA1p5.chrysalis_gnu: FAIL - appears to be related to 'inland seas' extrapolation issue fixed in this Add temperature and salinity to ocean extrapolation MALI-Dev/E3SM#124

[matthewhoffman]

With the workarounds in the most recent commit (1e3b2e1), the 3 Greenland tests that failed now pass. The larger issues related to those failures are outside the scope of this PR, so I can consider this PR ready to merge. I will move this PR to the main repo and any further evaluation can occur there.

[matthewhoffman]

PR moved to E3SM-Project#7215