Add a regional sea-level projection capability to MALI

[hollyhan]

This PR adds a 1D gravitationally consistent sea-level model (Kendall et al., 2005; Gomez et al., 2010; Han et al., 2022; https://github.com/GomezGroup/1DSeaLevelModel_FWTW) as a submodule within MALI, adding the capability of regional sea-level projection associated with ice sheet changes solved by MALI. In this new addition, the sea-level model computes changes in the height of sea surface equipotential and viscoelastic solid Earth associated with surface (ice and ocean) loading changes, enabling MALI to take into account gravitational, rotational and deformational effects of the loading changes.

The implemented requirements and modified/added codes in the PR are as follow:

1. MALI to optionally build SLM from the top level directory (mpas-albany-landice): modified Makefile in MALI's top level directory such that SLM=True in make command builds the SLM.
2. MALI to initialize SLM at t=0 of a simulation: added wrapper subroutine 'slmodel_init' for initializing and the SLM in mpas_li_bedtopo.F in which sea-level model sub-modules are called.
3. MALI to call SL solver at every coupling time interval: added wrapper subroutine 'slmodel_solve' for calling the SLM in mpas_li_bedtopo.F in which sea-level model sub-modules are called.
4. MALI to call the SLM at every coupling time interval (i.e., SLM timestep): 1) Modified Registry.xml to include a namelist value of the coupling time, 2) modified mpas_li_core.F to add an alarm to the land ice core clock object for the coupling time interval, and 3) modified mpas_li_bedtopo.F to use mpas_timekeeping and check/reset the alarm for the coupling time interval.
5. MALI needs to accommodate SLM running on a single processor while running on multiple processors: Implement MPI scatter, gather and halo updates in mpas_li_bedtopo.F. MALI data (ice thickness) is gathered on the head node after which the SLM is called.
6. Interpolation between MALI and SLM native grids should be done: Scripfiles mpas_to_grid.nc and grid_to_mpas.nc for the regional, unstructured MALI mesh and the global Gaussian grid are generated using MPAS-mesh tool from_mpas.py and ncremap, respectively. Subroutines interpolate_init and interpolate and MPI scatter/gather functions are copied from the ocean SAL PR#4472 [https://github.com/Add in-line self attraction and loading for global tidal simulation E3SM-Project/E3SM#4472]
7. MALI and SLM must exchange model data: changes on the head node, MALI provides thickness as input to SLM through SLM subroutine sl_solver that calculates and outputs total bedtopography, which is used to update bedrock topography in mpas_li_bedtopo.F
8. Include the SLM namelist file such that SLM parameters can be accessed and modifed in MALI top-level directory and that the SLM does not need to be recompiled the parameters are modified.

[matthewhoffman]

@hollyhan , thanks for opening this PR! Can you update the PR title to something more descriptive, and expand the description to be more detailed? For a major new functionality like this PR, we want the PR description to first give a summary of the overall capability, then describe the implementation, which subroutine were modified or added, and any important algorithmic details (e.g. copying the remapping and gather/scatter code from the ocean SAL PR - you can include a link to that PR). I'll work on a more detailed line by line review starting today, but I thought I would mention that before I got started.

[matthewhoffman]

@hollyhan , overall this PR looks great and represents a lot of work! I made a lot of comments, but they are all little things - many of which I should have instructed you in before you started working. :) A few of the items might require some discussion, but most should be self explanatory. After the initial clean up pass, I'll take a closer look at a few details and might have more changes to request.

The procedure for making changes is to add new commits that fix the various issues. For the formatting and little changes, feel free to put many or all of the changes into a single commit. For more significant changes, keep them in a separate commit. After you've taken care of everything you can, push your updated branch to Github, and the PR will update automatically with your revisions. When you are ready for me to take another look, request another review.

[matthewhoffman]

We'll need to update this to the project's public fork before we can merge.

[hollyhan]

Thank you for the very insightful comments, @matthewhoffman ! I have resolved most of the comments and updated the PR.

[matthewhoffman]

@hollyhan , thanks for all your hard work on revising this PR last week! It is much easier to read and things look a lot more mature. After a second reading, I've found a smaller number of additional changes. I think with another round or two we will be close to done!

[matthewhoffman]

I'm seeing a number of whitespace deletions in this file that are not related to any code changes. It's best to remove those if possible. I can show you how to see where they are showing up.

[matthewhoffman]

I'm still seeing about dozen whitespace changes in this file. You might need to go to https://github.com/MALI-Dev/E3SM/pull/21/files to see what I mean.

[hollyhan]

Thank you so much for the second round of review, @matthewhoffman!
I've addressed all of the points of the conversation - I think it's ready for the next round. :)

[matthewhoffman]

@hollyhan , this is looking really good! There are just a handful of small things in this round. The one bigger thing is I realized we should add a finalize routine that deallocates all the allocatable arrays. Let me know if you info on the details of how to do that. Otherwise, I think we are very close being done! We'll need to have another discussion about the SLM repo and make sure all the details about that are sorted out before we can merge.

[matthewhoffman]

I'm still seeing about dozen whitespace changes in this file. You might need to go to https://github.com/MALI-Dev/E3SM/pull/21/files to see what I mean.

[hollyhan]

Thanks for the another round of review, @matthewhoffman! Most of the suggestions were easy to address. One thing I would need to discuss in more detail with you is adding the new li_bedtopo_finalize routine. I'll reach out to you!

[hollyhan]

Thanks again for the review, @matthewhoffman! I addressed all of the requested changes and updated the PR.

[matthewhoffman]

@hollyhan , I just have a couple small changes this time around. After that, because this branch has been around for so long, let's rebase it on the current version of MALI-Dev/develop and then also remove commit a2771ad, which won't be needed after the rebase. You can do both operations at once with an interactive rebase. I'm happy to walk you through it if needed.

[matthewhoffman]

fix indentation here.

[trhille]

@hollyhan, are these panel descriptions correct?
[Top panel: bedTopography at simulation year 50 on 8 procs] [Middle panel: bedTopography at simulation year 0 on 1 proc] [Bottom panel: bedTopography at simulation year 50 on 1 proc] ![Screen Shot 2022-03-30 at 1 43 43 PM]

[hollyhan]

@trhille , yes they are but in different colour scale. I re-uploaded the plots so they are now all in the same colour scale.

[trhille]

@hollyhan, can you explain how to set the length of the simulation? I'm using these settings in namelist.landice:

config_start_time = 2007-07-01_00:00:00
config_stop_time = '2107-07-01_00:00:00'

and these in namelist.sealevel:

&time_config
    itersl = 1
    starttime = 0
    dt1 = 5
...
&timewindow_config
    L_sim = 100
    dt1 = 5
    dt2 = 10
    dt3 = 10
    dt4 = 10
    Ldt1 = 20
    Ldt2 = 0
    Ldt3 = 0
    Ldt4 = 0

It ran for 55 years and then ended with:

At line 288 of file sl_init_mod.f90 (unit = 201)
Fortran runtime error: Cannot open file 'INPUT_FILES/icemodel/iceload5g11': No such file or directory

Does the number of iceload5g files limit the length of the run? And if so, how do I set up a run of arbitrary duration?

[trhille]

From my 55-yr runs, I can confirm that the issues with noisy or unstable bedTopography have been resolved by 4319b7b, that bedTopography is evolving, and that bedTopography is bit-for-bit between simulations on 1 processor and 8 processors.

Runs started in 2007. left: 1 processor; middle 8 processors; right: difference between the two:

1 processor, to show evolution:

[hollyhan]

Need to clean src/SeaLevelModel directory in these lines: https://github.com/MALI-Dev/E3SM/blob/hollyhan/mali/add_SLmodel/components/mpas-albany-landice/src/Makefile#L54-L63

@trhille, I realize that @matthewhoffman had previously made make clean work for src/SeaLevelModel' with make clean SLM=true. I've modified the makefile slightly so that simply doing make cleanwill also clean files insrc/SeaLevelModel` if the SLM submodule exists. https://github.com/MALI-Dev/E3SM/blob/hollyhan/mali/add_SLmodel/components/mpas-albany-landice/src/Makefile#L64-L68\

[hollyhan]

It seems that the variable lot (number of data vectors) is not defined, or is defined incorrectly, leading to that upper bound error I saw arising in L1720 of spharmt.f90. We should define this variable somewhere. Since this is a bug in legacy code and not in your coupling, this could be a separate PR if needed.

@trhille - thanks for catching this. This particular piece of code was adopted by the group that developed original form of the sea-level model that we are using. Leaving it as is has been no problem so far because the SLM has been compiled without -fbounds-check and its results were benchmarked to the ancestral version of sea-level model. I do also think that we need to address this issue at some point and apply the same debug flags as we do in compiling MALI for compiling the SLM. Since this code belongs to the SLM submodule, we'll make a PR directly in the SLM repository and update the submodule here once the issue has been addressed.

[hollyhan]

From my 55-yr runs, I can confirm that the issues with noisy or unstable bedTopography have been resolved by 4319b7b, that bedTopography is evolving, and that bedTopography is bit-for-bit between simulations on 1 processor and 8 processors.

Thanks for testing again and confirming this, @trhille!

[hollyhan]

The coupling interval is set in both MALI and SLM settings, and there is no check on consistency. This leads to an inconsistency in the model time between MALI and SLM if they are not set to the same value.

Thanks for noting this, @trhille. I've added a check to the code that checks the consistency between coupling interval setup by MALI and SLM. 07b6832

[hollyhan]

@hollyhan, can you explain how to set the length of the simulation? I'm using these settings in namelist.landice:

config_start_time = 2007-07-01_00:00:00
config_stop_time = '2107-07-01_00:00:00'

and these in namelist.sealevel:

&time_config
    itersl = 1
    starttime = 0
    dt1 = 5
...
&timewindow_config
    L_sim = 100
    dt1 = 5
    dt2 = 10
    dt3 = 10
    dt4 = 10
    Ldt1 = 20
    Ldt2 = 0
    Ldt3 = 0
    Ldt4 = 0

It ran for 55 years and then ended with:

At line 288 of file sl_init_mod.f90 (unit = 201)
Fortran runtime error: Cannot open file 'INPUT_FILES/icemodel/iceload5g11': No such file or directory

Does the number of iceload5g files limit the length of the run? And if so, how do I set up a run of arbitrary duration?

How you have defined the length of simulation is "almost" correct - you would also want to set Ldt1 same as the L_sim if you don't want to use the time window algorithm. I realize this is really not straight forward just by looking at the namelist because there's no description of the name of variables. You can get more detailed description of the sea level namelist variables here: https://github.com/MALI-Dev/1DSeaLevelModel_FWTW/blob/main/user_specs_mod.f90#L134-L135

The number of iceload files (e.g., iceload5g in inputfolder_ice in namelist.sealevel) does limit the length of the run in an indirect way. Those input ice files are globally-defined ice thickness outside (they can be from other ice model results/reconstructions), and they need to be pre-processed/generated before the simulation starts. They get merged with ice thickness passed by MALI at every coupling timestep when the SLM is called (i.e., iceload5g0 will get merged with ice thickness over the Thwaites region that MALI provides at the initial timestep, iceload5g1 will get merged with MALI ice thickness at model year 5 when the coupling interval is 5 years... iceload5g10 will be merged with model year 50.) So when MALI was calling the SLM at year 55, the SLM is looking for iceload5g11. I've only put 10 ice files in the inputfile for the test purposes, and that is why the run stopped. What you can do for now if you want to test the run for 100 years is to increase the coupling time interval to 10 years from 5 or copy one of the input files like iceload5g0 and create iceload5g11...iceload5g20 files (this should be fine because I prescribed the iceload5gX files to have fixed ice thickness so the sea-level change pattern calculated by the SLM is only associated with the thickness changes in the Thwaites region calculated by MALI.)

[matthewhoffman]

This looks like it is assuming the config_slm_coupling_interval option has zeros for month, day, hr, min, sec. I think we might want the option to run with coupling intervals of smaller than a year to do sensitivity tests to the coupling interval. We can use the MPAS timekeeping operations to parse the string more carefully.

[matthewhoffman]

Given that the SLM just uses decimal years for timekeeping, it would be really tricky to ensure an exact match in the time interval (due to MPAS having the option to track leap years, MPAS keeping track of differing days per month, etc.). I think a match in time interval within 1% would be a reasonable approximation to avoid potentially complicated accounting here.

The most consistent way to do the calculation would be to use mpas_set_timeInterval and mpas_get_timeInterval in components/mpas-framework/src/framework/mpas_timekeeping.F to return the MPAS time interval in seconds (the dt output argument), and compare that to the SLM time interval. Then throw an error if they differ by, say, >1%. One complication there is that for time intervals that require months and/or years, mpas_get_timeInterval needs to know the start time, because the result will depend on the start time for many time intervals (e.g., 1 month). For the purposes of this calculation, the start time is not important, but you'd have to provide something. It might be best to hard code a time, like 2001-01-01_00:00:00 (I picked an arbitrary start of a year that isn't a leap year.). That is probably better than using something like the model start time, which would result in possibly slightly different results depending on start year.

This ends up being rather complex for a simple comparison, but I don't see a good way around it given the two different timekeeping systems in the two models.

[matthewhoffman]

@hollyhan , when I compile and run in debug mode, I get the following error:

At line 1720 of file spharmt.f90
Fortran runtime error: Index '513' of dimension 1 of array 'a' above upper bound of 512

If disable debug mode, the error goes away. Have you seen this?

[hollyhan]

@hollyhan , when I compile and run in debug mode, I get the following error:

At line 1720 of file spharmt.f90
Fortran runtime error: Index '513' of dimension 1 of array 'a' above upper bound of 512

If disable debug mode, the error goes away. Have you seen this?

@matthewhoffman , Yes - that's the error message from the module code spharmt.f90 in SLM, which was adopted by the group in harvard when the new original ice-age version of the SLM was written. It's a bound error, and it doesn't affect the calculation, which is why I think no efforts were made to correct it by the original developers. Trevor ran into this issue during his check, and solved by make clean-ing the SLM directory (#21 (comment)), which is now done in the top level directory (1469bf4)

[matthewhoffman]

I've tested the COMPASS full_integration suite in three ways, all with ALBANY=true:

1. SLM=false DEBUG=true
2. SLM=true DEBUG=false (setting DEBUG to false because of the error I get when both are true)
3. SLM=true DEBUG=true (none of existing tests include running the SLM, so the out of bounds error does not occur)
   In all three cases, all existing tests pass. None of the existing tests test the new functionality of course, but I can confirm that the changes do not affect any existing functionality.

I also confirmed:

• the expected error message occurs if you try to compile with SLM=true but did not initialize the SLM git submodule.
• if you compile without SLM but try to run with it you get a useful error message

I will proceed to testing runs that actually use the SLM with the Thwaites configuration Holly set up.

[trhille]

ice-age version of the SLM was written. It's a bound error, and it doesn't affect the calculation, which is why I think no efforts were made to correct it by the original developers. Trevor ran into this issue during his check, and solved by make clean

@hollyhan That fix was when using DEBUG=false. So it seems that needs to be fixed for debug mode.

[trhille]

@hollyhan, great work on this! After discussion with @matthewhoffman, I'm approving this PR with the knowledge that a few clean-up items will be addressed in another PR.

[matthewhoffman]

@hollyhan, I finished testing with the Thwaites mesh and see the same results as you and @trhille. I agree this looks really good and works as desired. It also nicely avoid affecting model behavior when the sea level model is not in use. Thanks for all your hard work in all the rounds of revisions and unexpected issues that popped up. We discussed the few clean-up items that we can leave for a follow up PR, but this is ready to merge!