Katetc/hma glaciers4, Replace Lipscomb/hma_glaciers4 PR

[Katetc]

Moved commits from #63 to a branch with a less complicated merge history to get rid of duplicate commits. Original branch with appropriate dates and commit history found here: lipscomb/hma_glaciers4 This is the branch used to develop a glacier modeling scheme. This code was used for the runs we submitted to GlacierMIP3 for the European Alps. The glacier scheme will be described in full in a paper in prep by Minallah, Lipscomb, and Leguy. Briefly:

• CISM reads in glacier information for an RGI region or sub-region. Each glacier-covered cell is assigned an RGI integer ID, ice thickness, and basal topography. CISM converts the RGI IDs to local IDs, with glaciers numbered consecutively from 1 to N.
• CISM also reads in climate fields: typically, monthly means of surface air temperature and precipitation at a reference elevation. For the Alps, we had two input climates: a baseline climate from the 1980s with glaciers assumed to be in balance, and a warmer, more recent climate when glaciers are out of balance. The fields for the warmer climate are specified as anomalies relative to the baseline climate. Temperature is downscaled to the current glacier elevation based on a fixed lapse rate. Precip can be applied as either rain or snow based on the downscaled temperature.
• In addition, CISM reads in an observationally based SMB dataset; we have used the dataset from Hugonnet et al. (2021).
• The goal is to spin up the glaciers such that the ice thickness and extent are in good agreement with observationally-based targets for the baseline climate. The targets are typically the RGI glacier extents, along with the thickness that would be expected at a date (perhaps earlier than the RGI date) when the glaciers are in balance with the climate.
• During the spin-up, we adjust two parameters, alpha and mu, for each glacier to satisfy two equations: (1) SMB = 0 = alphaP - muTp over the target glacier extent for the baseline climate, and (2) SMB' = SMB_obs = alphaP' - muTp' over the target glacier extent for the recent warm climate. Here, P is the precip and Tp is proportional to the difference between the downscaled temperature and the freezing temperature.
• During the spin-up, the ice flow dynamics is the same as it would be for ice sheets, but at higher resolution (100–200 m for the Alps). For glaciers we've been using the DIVA velocity solver with power-law basal sliding. The parameter C_p (powerlaw_c) is adjusted to optimize the match with the thickness target. The spin-up typically takes several thousand years. Toward the end of the spin-up, we may turn off the C_p inversion.
• Various logic is applied during the spin-up to assign appropriate integer IDs as glaciers advance and retreat, and to discourage excessive advance and retreat.
• Once the glaciers are spun up, we can run forward with a new climate.

Much of the new code is in the new module glissade_glaciers.F90. There are some new glacier fields (some 1D, others 2D) in glide_vars.def, and there is a new 'glaciers' section of the config file.

There is a new restart option 2: "hybrid restart". For this kind of restart, CISM initializes itself using the restart file from a previous run, but with a new start and end date specified in the config file.

There are some new I/O capabilities. For example, it is possible to read in all the input climate data (e.g., monthly data for 20 different years, a total of 240 time slices) once at initialization and store the data for repeated application during the run (rather than read in the same monthly data multiple times).

Glacier-related diagnostics (e.g., total glacier area and volume in the region) are written to the log file along with the standard diagnostics.

[whlipscomb]

@Katetc: I've run several tests comparing this branch to the current main branch. With LIVVkit I ran the dome (s0, s1, s2), ismip-hom (a, c, f), shelf-circular, shelf-confined, and stream tests. These all have BFB matches. I also got BFB matches for a dome DIVA test and an ISMIP6-type Antarctic spin-up.

I then ran a short Alps 200 m spin-up, comparing this branch to hma_glaciers3 (the branch used for GlacierMIP3, before the rebase to main) and hma_glaciers4 (after the rebase). Compared to hma_glaciers3, results are not quite BFB but are the same within roundoff (which is good enough). Compared to hma_glaciers4, results are BFB as expected.

As we talked about, there are a few diagnostic print statements to put back in glide_setup.F90. Once those are done, please go ahead and do the merge to main. Thanks again for all your work to clean up the history.