Lipscomb/hma glaciers4

[whlipscomb]

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.

[gunterl]

Note 1: When merging this branch, attention needs to be paid to not reverse the configuration modification in the test directory during PR#61 (such as powerlaw_c -> powerlaw_c_const).

Note2: When running the test cases MISMIP+, MISMIP3d, ismip-hom, and slab, results were BFB for thickness, masks,... but not for velocity fields (although the difference was on the order of 1e-11). In MISMIP+, MISMIP3d the differences span a few grid cells while it was field-wide for the slab test.