Lipscomb/hma glaciers4 #63
Open
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Until now, the exponent n in the Glen flow law has been set in glimmer_physcon.F90
as an integer parameter, gn = 3.
With this commit, n is renamed n_glen (a more greppable name) for use in Glissade.
It is declared in glimmer_physcon.F90 as real(dp) with default value 3.0d0.
Since n_glen is no longer a parameter, it can now be set in the config file like other
physical 'constants' (e.g., rhoi and rhoo) that are not truly constant, but can
take different values in different models or benchmarking experiments.
To avoid changing answers in old Glide code, I retained the integer parameter gn = 3 in glimmer_physcon.F90.
This parameter is now used only in the Glide dycore (e.g., glide_velo.F90).
In Glissade, I replaced gn with n_glen in several places:
(1) In subroutine glissade_interior_dissipation_sia, the dissipation factor includes n_glen.
Note: n_glen does not appear explicitly in the 1st-order dissipation, which is proportional to tau_eff^2/efvs.
(2) In glissade_velo_sia, n_glen appears in the vertical integral for the velocity.
(3) In glissade_velo_higher, flwafact = flwa**(-1./n_glen) where flwa = A.
(4) In glissade_velo_higher, the exponent p_effstr = (1.d0 - n_glen)/n_glen is used
to compute effective_viscosity for BP, DIVA, or L1L2.
(5) In glissade_velo_higher, subroutine compute_3d_velocity_l1l2 has a factor proportional to
tau**((n_glen-1.)/2.) in the vertical integral.
For (1) and (2), n_glen was previously assumed to be an integer. Switching it to real(dp)
is answer-changing at the machine roundoff level for runs using the local SIA solver
in glissade_velo_sia.F90.
For (3), (4) and (5), n_glen replaces the equivalent expression real(gn,dp).
For this reason, answers are BFB when using the BP, DIVA or L1L2 solver.
In glissade_basal_traction, n_glen appeared in the expression for beta in the Coulomb sliding option,
HO_BABC_COULOMB_FRICTION. Here, I replaced n_glen with powerlaw_m (also = 3.0d0 by default)
to be consistent with the expressions for beta in the Schoof and Tsai laws.
In glimmer_scales, Glen's n appears in the expressions for the scaling parameters vis0 and vis_scale,
Here, I defined a local integer parameter gn = 3 so that the scales are unchanged.
It is now possible for the user to specify arbitrary values of n_glen in tests such as the slab test.
Another minor change: I added some logic to the subroutine that computes L1L2 velocities.
For which_ho_efvs = 2 = HO_EFVS_NONLINEAR, the effective viscosity (efvs) is computed from
the effective strain rate using an equation from Perego et al. (2012).
But for option 0 (efvs = constant) and option 1 (efvs = multiple of flow factor),
this strain rate equation in the code does not apply. For these options, I added an
alternative that computes velocity in terms of the strain-rate-independent efvs.
This allows us to use L1L2 for problems with constant efvs (e.g., the slab test).
The code was calling subroutine init_isostasy with isostasy = 0 = ISOSTASY_NONE. This subroutine is now called only if isostasy = 1 = ISOSTASY_COMPUTE.
I modified glissade.F90 to abort cleanly with a call to glide_finalise after an advective CFL error. This is done only when the user does *not* specify adaptive subcycling. The clean abort allows the new slabStability script to keep going, launching a new run with a shorter timestep. In subroutine glissade_flow_factor of glissade_therm.F90, I removed the FLWA_INPUT option (option 3 of whichflwa). This option is redundant with option 0, FLWA_CONST_FLWA, in which the user can set default_flwa in the parameters section of the config file, and otherwise CISM uses default_flwa = 1.0e-16 Pa^-n yr^-1. We probably should rename default_flwa to constant_flwa, but leaving it for now to avoid breaking config files in test cases. Note: This option was added by Matt Hoffman in 2014. I am unaware of test cases that use this option (most have flow_law = 0 or 2), but if there are any, we will need to fix them by switching to whichflwa = 0. In subroutine glissade_therm_driver of glissade_therm.F90, I increased the threshold for column energy conservation errors from 1.0d-8 to 1.0d-7 W/m^2. For very small timesteps of ~1.0e-6 yr, the smaller threshold can be exceeded because of roundoff errors. In subroutine glissade_check_cfl of glissade_transport.F90, I modified the fatal abort for large CFL violations (advective CFL > 10). Now, CISM aborts for CFL > 10 only when adaptive_cfl_threshold > 0, i.e. transport subcycling is enabled. This prevents excessive subcycling for egregious CFL violations. If adaptive_cfl_threshold = 0. (the default), i.e. transport subcycling is not enabled, then the code exits cleanly (with a call to glide_finalise) in glissade.F90 when advective CFL > 1. I added a TODO note to set the default value of geot (the geothermal heat flux) to 0 instead of 0.05 W/m^2. With the default value, simple prognostic tests like the dome are not mass-conserving. Not changing just yet because answers will change for the dome test.
This commit modifies the run and plot scripts for the Dukowicz slab test case, as described in Section 5 of this paper: J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34, https://doi.org/10.5194/tc-6-21-2012. The test case consists of an ice slab of uniform thickness moving down an inclined plane by a combination of sliding and shearing. Analytic Stokes and first-order velocity solutions exist for all values of Glen's exponent n >= 1. The solutions for n = 1 are derived in Dukowicz (2012), and solutions for n > 1 are derived in an unpublished manuscript by Dukowicz (2013). These solutions can be compared to those simulated by CISM. The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman with support for n = 1. They came with warnings that the test is not supported. The test is now supported, and the scripts include some new features: * The user may specify any value of n >= 1 (not necessarily an integer). The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant). * Physics parameters are no longer hard-coded. The user can enter the ice thickness, beta, viscosity coefficient (mu_n), and slope angle (theta) on the command line. * The user can specify time parameters dt (the dynamic time step) and nt (number of steps). The previous version did not support transient runs. * The user can specify a small thickness perturbation dh, which is added to the initial uniform thickness via random sampling from a Gaussian distribution. The perturbation will grow or decay, depending on the solver stability for given dx and dt. For n = 1, the viscosity coefficient mu_1 has a default value of 1.e6 Pa yr in the relation mu = mu_1 * eps((1-n)/n), where eps is the effective strain rate. For n > 1, the user can specify a coefficient mu_n; otherwise the run script computes mu_n such that the basal and surface speeds are nearly the same as for an n = 1 case with the mu_1 = 1.e6 Pa yr and the same values of thickness, beta, and theta. (Note: There is a subtle difference between the Dukowicz and CISM definitions of the effective strain rate; the Dukowicz value is twice as large. Later, it might be helpful to make the Dukowicz convention consistent with CISM.) I modified the plotting script, plotSlab.py, to look in the config and output files for physics parameters that are no longer hardwired. I slightly modified the analytic formulas to give the correct solution for non-integer n. This script creates two plots. The first plot shows excellent agreement between higher-order CISM solutions and the analytic solution for small values of the slope angle theta. For steep slopes, the answers diverge as expected. For the second plot, the extent of the y-axis is wrong. This remains to be fixed. I also added a new script, stabilitySlab.py, to carry out stability tests as described in: Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance of depth-integrated ice-dynamics solvers, to be submitted to The Cryosphere. The idea is that for a given set of physics parameters and stress-balance approximation (DIVA, L1L2, etc.), the script launches multiple CISM runs at a range of grid resolutions. At each grid resolution, the script determines the maximum stable time step. A run is deemed stable when the standard deviation of an initial small thickness perturbation is reduced over the course of 100 time steps. A run is unstable if the standard deviation increases or if the model aborts (usually with a CFL violation). I have run the stability script for several solvers (DIVA, L1L2, SIA, SSA) for each of two physical cases: one with thick shearing ice and one with thin sliding ice. Each suite of experiments runs in a few minutes on 4 Macbook cores for solvers other than BP. As expected, DIVA and SSA are much more stable than L1L2 and SIA. This and the previous commit correspond to the CISM code and scripts used for the initial submission by Robinson et al. (2021).
Rewrote the slab README file to describe the new command line options for runSlab.py, and the new script stabilitySlab.py.
This commit introduces new algorithms for computing the 3D velocity field in the L1L2 solver. Also, I fixed a bug so that the L1L2 solver works correctly for n_glen = 1. The goal was to improve the stability of the L1L2 solver in slab tests. At fine grid resolution (< 200 m), the solver is unstable even for very small dt, instead of following the stability limits of an SIA solver (as is the case for Yelmo). The default method (method 1) for the 3D velocity computes most quantities on the staggered grid and is unchanged. There are two new methods: * Method 2 follows the Yelmo algorithm as closely as possible, with some quantities on cell edges and some on vertices. Unlike Yelmo, there is necessarily an averaging of uvel and vvel from edges to corners at the end, since CISM assumes B-grid velocities. * Method 3 moves all intermediate calculations to cell edges, with a final averaging to vertices at the end. To support the new methods, there is a new subroutine, glissade_average_to_edges, in module glissade_grid_operators. Both new methods give results similar to method 1. Neither improves stability at very high resolution. This suggests that the stability of Yelmo might result from the overall C-grid velocity scheme rather than details of the 3D velocity computation. For each method, there is now an option to exclude the membrane stresses in the tau_xz and tau_yz terms that appear in the vertical integration factor. When these stresses are excluded, the stability curve for L1L2 solver is close to that of an SIA solver, like Yelmo.
This commit includeds the following changes: * I fixed some typos in the README file for the slab case. * In runSlab.py, I added the option to use the local SIA solver (contained in module glissade_velo_sia.F90). * I added some code which can apply an initial sinusoidal perturbation of wavelength 2*dx, instead of a random Gaussian perturbation. This is useful for getting reproducible results.
This commit adds a hybrid velocity solver as described by Robinson et al. (TC, 2021). The solver first computes SIA velocities using the local SIA solver (which_ho_approx = -1) with zero basal slip, then computes SSA velocities using the Glissade higher-order SSA solver (which_ho_approx = 1), and finally adds the two solutions. The new logic is in module glissade_velo. To use the hybrid solver, set which_ho_approx = HO_APPROX_HYBRID = 5 in the config file. For the slab test, the hybrid solver has the expected stability properties, aligned with SSA at coarse resolution and SIA at fine resolution. Answers are as expected for a dome test.
This commit streamlines the 3D velocity subroutine for L1L1 and makes a small but important change in the algorithm. The change is to evaluate the membrane stresses in tau_xz and tau_yz at layer midpoints instead of lower layer boundaries, consistent with where the SIA terms are evaluated. I found that a small vertical offset can disrupt the balance between these two terms and make the L1L2 solver unstable for the slab problem at resolutions finer than ~200 m. With the fix, the L1L2 stability curve is parallel to the SIA stability curve at fine resolution, with L1L2 being slightly more stable than SIA. I also removed the alternate L1L2 discretization methods that were added in a recent commit. The alternate strategies evaluated some terms at edges instead of vertices, and did not change the results significantly. With this commit, all terms are evaluated at either cell centers or vertices. This is the code version used for the slab tests shown for CISM in Section 3 of the paper by Robinson, Goldberg & Lipscomb (2021, TC). For these tests, I temporarily changed the energy conservation criterion in glissade_therm.F90 to avoid false non-conservation alarms when running at very fine resolution. See the comments in that module. This commit is answer-changing for L1L2, in a good way.
This is the first of a series of commits toward an improved basal water scheme for Glissade.
It is a steady-state scheme, with basal water routed conservatively down a hydraulic gradient
from input cells (where bmlt > 0) to output cells. It is based on the serial code that was written
some years ago by Jesse Johnson for Glimmer.
Although this scheme is cruder than more recent, state-of-the-art schemes such as SHAKTI
(Sommers et al., 2018), which solve time-evolving equations for the hydraulic head,
the hope is that it will allow more realistic flow simply by putting basal water
in the right locations. The current local till scheme puts water only where there
is basal melting (not in downstream locations), and is not conservative.
To use the new scheme with Glissade, set which_ho_bwat = 3 = HO_BWAT_FLUX_ROUTING in the config file.
The driver is subroutine glissade_bwat_flux_routing, in module glissade_basal_water.F90.
The scheme has 4 steps:
(1) Compute the effective pressure N for each grounded grid cell.
For now, we assume N = 0, which is equivalent to assuming that local water pressure
balances the full ice overburden pressure.
(2) Compute the hydraulic head h for each grounded grid cell. This is given by
h = z_b + p_w / (rhow*g)
where z_b = bed elevation (m)
p_w = water pressure (Pa) = p_i - N
p_i = ice overburden pressure = rhoi*g*H
N = effective pressure (Pa) = part of overburden not supported by water
H = ice thickness (m)
(3) Route basal water down the gradient of hydraulic head h. This is done by
(a) filling any depressions in h(x,y)
(b) adding small increments to h so that all water has a downhill outlet
(c) sorting the grid cells in order from low to high h
(d) initializing F = bmlt*dx*dy in each cell
(e) looping through cells from high to low h, and for each cell, partitioning
the input bwatflx among one or more down-gradient cells
For step (e), there are three routing options:
(i) D8; water is routed to the down-gradient cell with the lowest h
(ii) Dinf; water is routed to the two down-gradient cells with the lowest h,
in proportion to the gradient
(iii) FD8; water is routed to all down-gradient cells in proportion to the gradient
The original Glimmer code contained only the FD8 option. I added the others since
they are less dispersive. The user can choose the method by setting a new
config parameter, ho_flux_routing_scheme (= 0 for D8, = 1 for Dinf, = 2 for FD8).
D8 is the default.
(4) Compute the steady-state water depth based on the simple expression:
F = q * dx
where F = total water flux (m^3/s), computed in step 3
q = water flux per unit width (m^2/s), a function of grad(h) and water depth b
dx = grid cell width
Following Sommers et al. (2018, Eq. 6), we assume the flux is given by
q = (b^3 * g) / [(12*nu)(1 + omega*Re)] * -grad(h)
where b = basal water depth (m)
nu = kinematic viscosity of water (m^2/s)
Re = Reynolds number (unitless)
omega = an empirical constant (unitless)
For small Re, the flow is laminar, and for large Re, the flow is turbulent.
For now, assume laminar flow (Re -> 0).
Given F from step 3 and grad(h) from step 2, it is straightforward to solve for b.
These equations are similar, but not identical, to what is assumed in the Glimmer/Glide version.
For this reason, I did not compare try to obtain the same answers as in Glide.
To make the routing scheme more versatile, I wrote two subroutines that were not
in the old Glimmer flux-routing scheme:
(1) fix_flats (step 3b above)
This subroutine uses the algorithm of Garbrecht & Mertz (1997, J. Hydrol.)
to increment the surface elevation in flat regions, ensuring that all cells
have a downhill outlet. It repeatedly calls subroutine find_flats.
See G&M (1997) for details.
(2) find_flats
This subroutine identifies all cells without a downslope gradient.
These are regions where the surface is flat, usually after filling depressions,
and the water has no downhill outlet.
I verified that fix_flats is working, first for the example problem in G&M (1997)
and then for a dome problem with a central depression in the hydraulic head.
For diagnostics, I added output fields head and bwatflx in a new basal_hydro derived type.
I moved bwat and stagbwat to this derived type, along with some parameters that are used
for the local till model. This led to minor code changes, replacing 'temper' or 'basal_physics'
with 'basal_hydro', in several modules. Other fields and parameters could be added later
to the basal_hydro type to support new basal hydrology models, such as SHAKTI.
I coded these equations and set up a simple dome test problem with a basal melting source
beneath the dome center, where h is high. I verified that water flows downhill conservatively;
i.e., the total output flux at the margin is equal to the input flux from bmlt.
I plotted F = bwatflx for the D8, Dinf and FD8 options. As expected, FD8 gives a fairly uniform,
diffuse flow, while D8 gives a sharper flow in several distinct streams.
Dinf seems not to work well for the dome geometry because many ties are broken asymmetrically.
Next steps:
* Implement a parallel scheme on multiple tasks.
This could be done simply using gathers and scatters, or scalably using halo updates and new logic.
* Try out in more realistic ice sheet problems.
This commit includes a number of changes to allow the new flux-routing scheme
to run on multiple processors. Thus requires passing the parallel derived type
to a number of subroutines and revising the flux-routing logic.
In several iterative calculations with a stopping criterion (e.g., finding depressions,
and steps 1, 2 and 3 of the fix_flats algorithm), I replaced local with global sums.
Now, the iteration ends when a global criterion is met (e.g., no more depressions).
I modified subroutine route_basal_water as follows:
- In the serial version, the code loops from high to low cells in one iteration.
For each cell, any incoming flux is routed to neighboring downslope cells.
When we get to the lowest cell, all the water flux has reached the margin.
- In the parallel version, there are multiple iterations.
In each iteration, we loop from high to low over the locally owned cells
(not halo cells) on each processor.
The first iteration includes any water flux from basal melting (bmlt > 0).
For each local cell, the flux is routed downhill. Two things can happen:
(1) The flux reaches the low-lying margin, In this case, we are done with it.
(2) The flux is routed to a downslope halo cell. In this case, the flux in the
halo cell is communicated to the neighboring processor, and then routed
downslope to the locally owned cell adjacent to the halo.
In the next iteration, halo fluxes computed on the previous iteration are routed downhill.
When all the water has reached the margin, the iteration halts.
The total water flux is accumulated on each iteration, so that when the iteration
is done, the outgoing flux reaching the global margin should be equal to the
incoming flux received from basal melting.
To support computations of global sums, I added a parallel_global_sum interface
in the parallel modules. This could be useful in other modules too.
I also added a new halo update subroutine, parallel_halo_real8_4d, to support
efficient halo updates for arrays with two dimensions other than ewn and nsn.
Answers with a serial build are the same as for an MPI built with np = 1.
I checked that depressions are filled and flats are fixed correctly on 4 processors.
I verified that answers are the same (within roundoff) on 1 versus 4 processors
for a dome problem with a simple hydraulic head field that has flow across processors.
It might be possible to make the routing algorithm more scalable, e.g. by reducing the
number of global sums. However, this might not be worth the effort, if the
flux-routing scheme remains much cheaper than the velocity solver.
This commit includes some bug fixes and other changes in the parallel flux-routing scheme:
* Fixed a bug in the Dinf routing option. Results now look sensible.
* Do not apply bmlt < 0 (i.e., refreezing) to the basal water flux.
To conserve water, refreezing will need to be handled separately.
* Changed the criteria for bwat_mask = 0. This is the mask that identifies cells
through which basal water can flow. We do not want to exclude ice-free cells
in the interior, but we want to count all water that leaves the ice sheet.
Now, bwat_mask is set to 1 in the following cells:
(a) floating or open ocean
(b) cells at the edge of the global domain. I added a subroutine,
parallel_global_edge_mask, to identify these cells at initialization.
(c) ice-free cells with overwrite_acab_mask = 1
* Added an option for overwrite_acab: OVERWRITE_ACAB_INPUT_MASK = 3.
With this option, cells where the SMB is overwritten (usually with a negative value)
are read from the input file. For the dome hydro problem, this mask can be used
to zero out basal water outside the original dome boundary.
* Increased count_max, the max number of iterations for the flux routing, to 50.
A new iteration is needed whenever there is a nonzero flux in one of more halo cells.
With Dinf or FD8, there can be multiple crossings of the boundary between processors
as water flows down the gradient. (Up to ~40 on a 4km Greenland grid.)
This commit contains several minor changes to support the flux-routing hydrology scheme in runs with Northern Hemisphere paleo ice sheets * Added a new effective pressure option: which_ho_effecpress = HO_EFFECPRESS_BWAT_RAMP = 5. This is similar to the existing option HO_EFFECPRESS_BWAT = 4. The only difference is that as bwat increases from 0 to bwat_till_max, the effective pressure ramps down linearly, unlike the more complex formulation of Bueler and van Pelt (2015). * Added a new logical option, smooth_input_usrf. When this option is set to true, the initial upper surface elevation field (usrf) is smoothed, and the thickness is then adjusted to be consistent. This is helpful for a 4-km N. Hemisphere simulation using the input file cism_USGS-huy3-S1D_4km.nc. In this file, parts of the GrIS have rough topography and fairly smooth thck, leading to rough usrf and large surface gradients that crash the solver in the first few time steps. With several smoothing passes, the code starts cleanly. The smoothing is done in subroutine glissade_smooth_usrf, in glissade_utils.F90. Note that usrf is smoothed only for grounded ice, and not for floating ice, ice-free ocean, or ice-free land. * Made the fill_depression routine more efficient by restructuring with an outer and inner loop. Halo updates are called only from the outer loop. Although more efficient than before, it is still very expensive to fill depressions with the current algorithm. The cost of the entire code more than doubles on a 4-km N. Hemisphere grid. In a future commit, I will try to implement a different algorithm that scales better. * Added the max value of bmlt_ground in the diagnostic log file
This commit adds an efficient algorithm for filling depressions in the head field
when running the flux-routing hydrology scheme.
The old scheme was very slow to converge on large grids such as the 4km
Northern Hemisphere grid.
The new scheme is based on the algorithm of Planchon and Darboux (2001).
The basic idea is:
* Initially, set phi = phi_in on the boundary, and set phi = a large number elsewhere
(where phi is the head field).
* Loop through the domain. For each cell c, with value phi(c) not yet fixed to a known value,
compute phi_min8(n), the current minimum of phi in the 8 neighbor cells.
- If phi_in(c) > phi_min8(n) + eps, then set phi(c) = phi_in(c) and mark that cell as having a known value,
since phi(c) cannot go any lower.
- If phi_in(c) < phi_min8(n) + eps, but phi(c) > phi_min8(c) + eps, set phi(c) = phi_min8(n) + eps.
Do not mark the cell as having a known value, because it might be lowered further.
* Continue until no further lowering of phi is possible. At that point, phi = phi_out.
Here, eps is a small number greater than zero, which ensures that there are no flat surfaces
when the depression-filling is done. Thus, it is no longer necessary to call fix_flats.
On the 4km N.H. grid, the number of depression-fill iterations is reduced from several hundred
per time step to ~10.
I ran the EISMINT-2 tests with Glide for the first time in recent memory
to see if they were still working. They were not. Experiment G, which
has substantial basal sliding, crashed with a checkerboard instability
in the thickness field.
The main problem was a sign error in glide_velo.F90, subroutine slipvelo.
The code originally contained this line:
model%velowk%fslip = rhograv * btrc
which I had replaced some time ago with this line:
model%velowk%fslip = rhoi * grav * btrc
However, the correct replacement includes a minus sign:
model%velowk%fslip = -rhoi * grav * btrc
There are a number of constants and fields in this module that are defined as negative.
I added some comments to clarify which expressions implicitly include minus signs.
A minor fix: I had added bmlt_ground as a diagnostic field without making sure
it is allocated for Glide runs. Now, bmlt_ground is always allocated.
For Glide, which has no melting beneath floating ice, bmlt_ground is simply
a copy of bmlt.
To each config file, I added diagnostics at the central dome point (31,31) every 5000 years.
This makes it possible to read relevant results (e.g., total area and volume, and
basal temperature at the divide) directly from the log file.
I changed bmlt_ground back to bmlt in the list of output fields.
I then ran the full suite of EISMINT-2 experiments with Glide. All runs finished cleanly.
I compared results from Experiments A, B, C, D, G and H with those in the paper
by Payne et al. (2000, J. Glaciol.).
Glide results are generally consistent with the published model means.
The paper does not say which specific model was Tony's early version of Glide,
but it probably was V, W or Z based on the spoke patterns in Fig. 11 for Expt H.
Glide has changed enough that we would not expect to duplicate his answers.
This commit is answer-changing (in a good way) for Glide runs with nonzero basal sliding.
This commit adds a directory with config files for running several EISMINT-2 experiments
(A, B, C, D, F, G and H) with the Glissade local shallow-ice solver.
These are the seven experiments discussed in the paper by Payne et al. (2000, J.Glaciol.).
Tests I, J, K and L are omitted.
The Glissade local SIA solver is similar to Glide, except that the continuity equation
is solved using an explicit upwind-weighted advection scheme (usually incremental remapping)
instead of an implicit diffusion scheme. See glissade_velo_sia.F90 for the solver algorithm.
The different settings, compared to the Glide tests, are as follows:
* dycore = 2 (Glissade)
* evolution = 3 (IR transport)
* vertical_evolution = 0 (not constrained by upper kinematic BC)
* which_ho_approx = -1 (Glissade local SIA solver)
* which_ho_sparse = 3
* ice_limit = 1 (minimum thickess for active ice = 1 m)
Note: The Glissade local SIA scheme does not need a sparse solver,
but the default value (which_ho_sparse = 0) triggers an error in parallel runs.
Tests are typically run as follows:
> ./cism_driver e2.a.config
where cism_driver is a symbolic link to the executable:
cism_driver@ -> ../../../builds/mac-gnu/cism_driver/cism_driver
The Glissade SIA solver is slower than Glide.
Unlike Glide, however, the Glissade solver can be run in parallel, e.g.
> mpirun -n 4 ./cism_driver e2.a.config
This commit modifies module eismint_forcing.F90 so that the global temperature
and mass balance forcing are computed correctly on multiple processors.
I ran experiments A and G to completion and compared to Glide results.
The results are different (e.g., Glissade gives a higher dome), but
are in the same ballpark as the results in Payne et al. (2000).
One config file, e2.a.config.diva, is included with typical higher-order settings
for the Glissade DIVA solver. Other config files could be modified in a similar way.
In some runs, a spurious nonzero tf_anomaly was appearing in subroutine glissade_bmlt_float_thermal_forcing, as a result of not being initialized to zero during each pass through the calling subroutine in glissade.F90. In the first few timesteps, the anomaly grew large enough to trigger a fatal error in the check for reasonable values of thermal forcing. This commit adds logic to initialize tf_anomaly and tf_anomaly_basin correctly. It's unclear why this bug didn't show up before.
This is an initial implementation of the sliding law suggested by Zoet & Iversion (2020, Science).
This law takes the form (see Eq. 3 in the ZI paper):
tau_b = C_c * N * [u_b/(u_b + u_t)]^(1/m), Eq. 3 in ZI(2020)
where C_c = a constant in the range [0,1], aka coulomb_C
N = effective pressure
u_t = threshold speed controlling the transition between powerlaw and Coulomb behavior
m = powerlaw exponent
Here, C_c and m correspond to tan(phi) and p, respectively, in ZI Eq. 3.
This law is implemented as which_ho_babc option HO_BABC_ZOET_IVERSON = 7,
replacing the unsupported option HO_BABC_YIELD_NEWTON.
C_c can be implemented either as a constant (default value = 0.42) or as a 2D field (coulomb_c_inversion).
For the latter, I implemented an inversion procedure, closely following the procedure for Cp inversion
with the Schoof and powerlaw options.
The threshold speed u_t is a new config parameter.
Tim van den Akker already added ZI on a branch created from main.
This commit replicates Tim's changes on a branch that is rebased onto the recent basal hydrology changes,
and adds code changes for C_c inversion.
Note: The new code has not yet been tested thoroughly.
This commit may be squashed with other commits after debugging.
This commit includes a number of changes in the options for computing the
yield stress in Coulomb-type friction laws. These are laws in which the
yield stress is written as N*C_c or N*tan(phi), where N is effective pressure,
C_c is a scalar in the range [0,1], and phi is a friction angle, such that
tan(phi) is in the range [0,1].
Note that C_c and tan(phi) are functionally equivalent, and CISM generally
uses the C_c terminology (coulomb_c in the code).
The goal is to allow mixing and matching of options for specifying or reducing
the yield stress: e.g. via ocean connection, or the presence of meltwater at the bed,
or permanently weak till.
The main changes are:
* There is a new set of options for which_ho_effecpress; some are removed and some are added.
* The ocean connection option is now applied simply by setting p_ocean_penetration > 0,
without needing to set which_ho_effecpress = 3.
* There is a new option called which_ho_coulomb_c for setting till strength.
* There is a new option which_ho_powerlaw_c, analogous to which_ho_coulomb_c.
* The options which_ho_cc_inversion and which_ho_cp_inversion are deprecated.
Instead, the user sets which_ho_coulomb_c/powerlaw_c = 1 to compute these
coefficients by inversion, and which_ho_coulomb_c/powerlaw_c = 2 to read these
coefficients from an external file (which may have been obtained by inversion).
* Added 2D fields coulomb_c_2d and powerlaw_c_2d to go with these new options,
superseding coulomb_c_inversion and powerlaw_c_inversion.
The new fields are part of the basal_physics derived type (instead of the inversion type).
With these changes, config files used for ISMIP6 experiments will not, in general,
run cleanly and give the same result with the new options. Several changes will be needed
in the CESM namelists.
More details are given below.
*****
The new options for which_ho_effecpress are:
[0] N is equal to overburden pressure, rhoi * g * H.
[1] N is reduced where the bed is at or near pressure melting point, with a linear ramp.
[2] N is reduced where basal water is present (bwat > 0), with a linear ramp.
[3] N is reduced where there is meltwater flux at the bed (bwatflx > 0), with a linear ramp.
[4] N is reduced where basal water is present (bwat > 0), following Bueler & van Pelt.
Options 0, 1, and 4 are the same as current options.
Option 2 is the same as the old option 5.
Option 3 is new, replacing the old ocean_penetration option.
The options for which_ho_coulomb_c are:
[0] C_c is equal to a spatially uniform constant, the parameter coulomb_c.
I changed the default value from 0.42 to 1.0. In other words, yield stress
is equal to overburden pressure by default, but can be overridden by
setting coulomb_c = 0.50, or a similar value, in the config file.
[1] C_c is found by inversion, nudging the simulated ice thickness toward a target thickness.
[2] C_c is a 2D field read from an external file. The values could be based on geology
or obtained from a previous inversion.
[3] C_c is a function of bed elevation, following the pseudo-plastic law as typically applied.
Note: Option 3 is not supported with this commit; will add support in an upcoming commit.
The options for which_ho_powerlaw_c are:
[0] C_p is equal to a spatially uniform constant, the parameter powerlaw_c.
The default value remains 1.0d4 Pa m^(-1/3) yr^(1/3), which is appropriate
for a power law with exponent m = 3.
[1] C_p is found by inversion, nudging the simulated ice thickness toward a target thickness.
[2] C_p is a 2D field read from an external data set. The values could be based on geology
or obtained from a previous inversion.
There is no analog to which_ho_coulomb_c = 3.
A conceptual point: The basal melt mechanisms that decrease N by increasing water pressure
could equally be viewed as mechanisms that decrease C_c by weakening the till.
Here, we view these mechanisms as decreasing N. Thus, N becomes a dynamic field
that can vary rapidly with changes in the basal water system or ice thickness, whereas
C_c is viewed as a slowly varying field embodying quasi-permanent properties of the till.
With which_ho_coulomb_c = 2 and isostasy turned on, C_c will vary with the bed elevation,
but this is generally a slow change.
If the user mistakenly sets which_ho_effecpress = 3 to specify ocean penetration with p > 0,
the code should get the same answer as before, provided bwatflx = 0.
But it is safer to set which_ho_effecpress = 0 in this case, while setting p > 0.
Another small change: Changed the units of bwatflx from m^3/yr to m/yr, so that this flux
is not proportional to the grid cell area.
I set up test cases for ISMIP6-style Antarctic spin-up runs and confirmed that the answers
are BFB with the appropriate choices of new options.
This commit implements option which_ho_coulomb_c = HO_COULOMB_C_ELEVATION = 3, with coulomb_c varying linearly between a max value at high bed elevations (>= bedmax) to a min value at low bed elevations (<= bedmin). This option generalizes the procedure used in the pseudo-plastic basal friction law to set tan(phi), which is equivalent to coulomb_c. For details, see the git log for the previous commit. To use the new option, users should set the config parameters coulomb_c_max, coulomb_c_min, coulomb_c_bedmax, and coulomb_c_bedmin, if not using the default values. The new coulomb_c option can be applied to any basal friction law with a coulomb_c term. Currently, these are friction laws 2 (new pseudo-plastic), 7 (Zoet-Iverson), 10 (basic Coulomb), 11 (Schoof), and 12 (Tsai). Previously, bed elevation dependence was hardwired for the pseudo-plastic friction law (the default in CESM Greenland runs), but was not supported for other friction laws. For backward compatiblity, CISM now supports two flavors of pseudo-plastic law: * HO_BABC_PSEUDO_PLASTIC_OLD = 3, the old option, with hardwired elevation dependence * HO_BABC_PSEUDO_PLASTIC = 2, the new option, with optional elevation dependence The older option can be removed when it is no longer the CESM default. To make room for the new option as #2, I set the little-used HO_BABC_YIELD_PICARD = 15. To streamline the code, powerlaw_c_2d and coulomb_c_2d are now computed near the start of subroutine calcbeta, instead of being computed separately for each friction law. Note: While adding the new option, I found that the old option has a minor bug: tanphi is a function of topg(ew,ns). However, tanphi should be colocated with beta on the staggered grid, whereas topg lives on the unstaggered grid. I left the bug in place for backward compatibility, given that this option will be deprecated going forward. Note: The default value of coulomb_c_min is 1.0e-3. For the pseudo-plastic law, a more appropriate value is coulomb_c_min ~ 0.10, which is close to tan(phimin) with phimin = 5 degrees. In runs with a pseudoplastic law, answers are BFB with which_ho_babc = 3, and are modestly different (as expected) with which_ho_babc = 2 combined with which_ho_coulomb_c = 3. In Antarctic runs with other friction laws, answers are BFB except when using the new options.
Previously, basal_physics%powerlaw_c and basal_physics%coulomb_c were constants, and basal_physics%powerlaw_c_2d and basal_physics%coulomb_c_2d were 2D fields. The usual CISM convention is not to put '2d' in the names of fields except to distinguish them from 3d fields (as in the case of uvel_2d, vvel_2d). Instead, we put 'const' or 'constant' in the names of parameters that need to be distinguished from 2D fields. With this commit, the constants are named basal_physics%powerlaw_c_const and basal_physics%coulomb_c_const. The 2D fields are simply basal_physics%powerlaw_c and basal_physics%coulomb_c. This commit is BFB. However, it requires changing config files. In the parameters section, powerlaw_c -> powerlaw_c_const and coulomb_c -> coulomb_c_const. In lists of output fields, powerlaw_c_2d (or the older powerlaw_c_inversion) becomes powerlaw_c, and similarly for coulomb_c.
When running ISMIP6 experiments, it was convenient to read in some parameters and fields from the input file. However, this has led to some confusion and errors. In particular, when gamma0 can either be read from an input or restart file or specified in the config file, there is a risk of overwriting the desired value. With this commit, gamma0 is no longer an I/O option. It is not listed in glide_vars.def or written to restart files. The only way to set gamma0 /= 0.0 (the default is 0.0) is to specify a nonzero gamma0 in the config file. Thus, gamma0 is now treated as most other parameters are treated. I also removed the bmlt_float_ismip6_magnitude option. This option was used to specify whether CISM was running with low-, median-, or high-sensitivity gamma0, as defined for ISMIP6, and thereby choose specific values of gamma0 and deltaT_basin. Six deltaT_basin variants (deltaT_basin_nonlocal_median, deltaT_basin_local_pct5, etc.) are no longer in the code. If present, deltaT_basin is still read from the input file, but the user must put the correct field in the input file, instead of letting CISM choose from among six possible fields. When inverting for deltaT_basin, there is no need to put it in the input file; it is initialized to zero everywhere. Finally, I removed the thermal_forcing_baseline field, which was no longer used but was simply copied to thermal_forcing. This commit is BFB, apart from removing some old options. Note: Some older Antarctic input files contain a field called thermal_forcing_baseline. This should be changed to thermal_forcing, else it will not be read in.
This commit removes most of the plume model that was developed several years ago but was never robust. The plume model was written to support option whichbmlt_float = BMLT_FLOAT_MISOMIP = 5, with the goal of computing sub-shelf melt rates given the prescribed MISOMIP temperature and salinity profiles. I left the option in the code, since it might be supported with a plume or reduced-order ocean model that is still to be developed. However, the code now aborts with a fatal error if the user selects this option. In glissade_bmlt_float.F90, I removed most of the plume code but kept a few subroutines that might be useful later.
This commit removes the old bmlt inversion option, which_ho_bmlt_inversion, which inverts for bmlt_float in floating grid cells, adjusting the melt rate to bring the thickness in each grid cell closer to an observational target. This scheme was used for the original ISMIP6 Antarctic runs but was never very robust. The melt rate is sensitive to the ice thickness, so it is hard to converge on stable values. The resulting melt rate field is noisier than real melt rates would be. The bmlt_basin inversion option (which_ho_bmlt_basin_inversion) remains. This scheme inverts for deltaT_basin in each basin and is more robust.
The steady-state flux-routing scheme computes the basal water flux (bwatflx) based on the assumption that the total water inflow from basal melting is equal to the outflow to the ocean. It then diagnoses bwat from bwatflx, given an expression for the water velocity in terms of grad(head). Previously, the effective pressure could be reduced in subroutine calc_effecpress based on either bwat or bwatflx. However, the bwat diagnosed in the flux-routing scheme is problematic for basal thermodynamics. The thermal solver assumes that the bed is at the melting point wherever bwat > 0. As a result, a very small bwat beneath frozen ice in the flux router can drive large, abrupt increases in basal temperature. The solution is to diagnose bwat locally in the flux router, without returning bwat to glissade. As a result, the option which_ho_effecpress = HO_EFFECPRESS_BWAT = 2 can no longer be used with the flux router. It should be used only with the local till model, which prognoses bwat. To reduce N when using the flux router, set which_ho_effecpress = HO_EFFECPRESS_BWATFLX = 3. For consistency, effecpress_bwat_threshold is now used for option 4 (Bueler-Van Pelt) as well as option 2 (linear ramp). The BvP scheme has a related parameter, bwat_till_max, which normally should be set to the same value as effecpress_bwat_threshold. Both parameters now have defaults of 2 m. I verified that a sample Greenland run with the pseudo-plastic law and local till is BFB. A minor change: I commented out the call to subroutine effective_pressure in the flux router, and added comments explaining that we assume N = 0 for the purpose of computing the head. With these changes, Antarctic simulations are more stable than before. However, with which_ho_effecpress = 3 and effecpress_bwat_threshold > 1 m/yr, only the largest Antarctic channels have significant reductions in N. We should consider replacing the linear ramp with other functional forms.
The flux-routing basal water scheme returns a 2D bwatflx field, which is used to reduce the effective pressure N. Previously, we used a linear ramp function to compute N from bwatflx. However, bwatflx can range over several orders of magnitude, and the marginal reduction in N per unit of added bwatflx could decrease gradually as bwatflx grows. With this commit, N is prognosed as a function of Fw/F0, where Fw = bwatflx and F0 is an empirical constant. More precisely, the prognosed quantity is called f_effecpress (or f_N for short), defined as the ratio of effective pressure to overburden pressure, N/N_o = N/(rhoi*g*H). That is, f_N is the fraction of overburden that is *not* balanced by water pressure. It lies in the range (0,1] and evolves as follows: df_N/dt = [1 - f_N*(Fw/F0)] / tau_N, where tau_N is a relaxation timescale. With 'f_N' in the numerator of the rhs, it becomes harder to reduce N as f_N approaches zero. The '1' in the numerator nudges f_N back toward 1 when Fw is small or zero. The asymptotic value attained with a given Fw is f_N = min(F0/Fw, 1). Note: Initially, I tried diagnosing N based on the instantaneous value of (Fw/F0), but this can lead to unstable oscillations in ice speed and thickness. For now, the default values of the new constants are tau_N = 100 yr and F0 = 0.01 m/yr. These can be set in the [parameters] section of the config file. I made a related change in the ZI law. This law has the form tau_b = N * C_c * [u_b/(u_b + u_t)]^(1/m). I replaced 'N' with min(N, N_max), where N is the actual effective pressure and N_max is another empirical constant. The idea is that for thick, slow-sliding ice in the power-law regime, the sliding coefficient (i.e., the term that multiplies u_b^(1/m)) can asymptote instead of continuing to increase with higher N. This makes the ZI law more like the Schoof law in the power-law limit. The default N_max = 1.e8 Pa, an overburden corresponding to an unrealistic ~10 km of ice. Thus, N is not capped in the ZI law by default, but only when the user sets a smaller N_max in the config file. Finally, I changed the way N is computed when it can be reduced by either basal melting (which_ho_effecpress >= 1) or by p > 0. Previously, N could be multiplied successively by two numbers < 1 when both effects were active. On further reflection, one reduction should not necessarily reinforce the other. The new computation is N = min(N(bwatflx), N(ocean_p)). That is, we choose the minimum value from the two calculations.
This commit fixes a bug in the computation of effective pressure from bwatflx. A halo update of bwatflx was missing, leading to an error in staggered effecpress and resulting ice speeds at processor boundaries. Thanks to Sarah Bradley for spotting the error. The fix is to add a halo update for bwatflx at the end of the flux-routing scheme, and/or a halo update for effecpress before computing effecpress_stag. Either update fixes the bug, but to be on the safe side, I added both updates. I also added and revised some basal water and effecpress diagnostic prints.
I added logic to prevent evaluating a term with thck_gradient_ramp in the denominator unless thck_gradient_ramp > 0. This turned up in recent debugging; not sure why it didn't turn up earlier.
Until now, the inversion for powerlaw_c and coulomb_c has used an
observational thickness target (or equivalently, a surface elevation target).
With this commit, it is possible to invert for powerlaw_c or coulomb_c based on
a thickness target, a velocity target (observed surface speed), or both.
The two targets are combined in one subroutine, with similar logic for each target.
The rate of change of C_p or C_c is proportional to:
(H - Hobs) / babc_thck_scale with a thickness target, and/or
(v - vobs) / babc_velo_scale with a velocity target.
The rate is damped by a term proportional to dH/dt for a thickness target,
but there is no damping proportional to dv/dt. (I tried adding a damping term
proportional to dv/dt, but it led to oscillations in coulomb_c.)
To support the velocity target, these parameters and fields have been added:
babc_velo_scale = scale for velocity differences
usfc_obs = x component of observed surface speed
vsfc_obs = y component of observed surface speed
velo_sfc_obs = observed surface speed = sqrt(usfc_obs^2 + vsfc_obs^2)
Note: babc_thck_scale = babc_velo_scale = 0.0 by default.
Previously (e.g., for ISMIP6 runs), babc_thck_scale was equal to 100 m by default.
To enable inversion based on thickness, velocity, or both at once,
the user can set either or both parameters to a positive value in the config file.
For runs with velocity inversion, I am using a value of 200 m/yr.
To reduce code duplication, I consolidated several subroutines in glissade_inversion.F90.
The same subroutines now handle both powerlaw_c and coulomb_c inversion,
with which_ho_powerlaw_c and which_ho_coulomb_c determining the appropriate arguments
to pass between subroutines.
I verified that answers are BFB with the consolidation.
I removed dthck_dt from the list of inversion restart variables.
It is not, in fact, needed for exact restart.
The scheme seems to be working as desired in a 3 ka Antarctic spin-up.
Still need to do a close analysis of runs with and without a velocity target.
One challenge in Antarctic spin-ups is to prevent the retreat of the East Thwaites grounding line early in the run. The ice can thin by several hundred meters and never recover. With p = 1, the entire Thwaites basin can collapse. With this commit, users can phase in the effect of ocean_p (aka p_ocean_penetration or simply 'p') on effective pressure N by specifying a relaxation timescale > 0. The timescale is a new parameter, basal_physics%ocean_p_timescale. When the timescale > 0, the effective pressure is relaxed using a new prognostic field, basal_physics%f_effecpress_ocean_p. This is a scalar field in the range [0,1] that specifies the fractional value of N relative to overburden. If p > 0, then N is multiplied by f_effecpress_ocean_p for each gridcell. Typically, the original N is the overburden pressure rhoi*g*H, unless N is also reduced by meltwater at the bed (e.g., with a hydrology model). With a timescale of zero (the default value), N = rhoi*g*H*(1 - Hf/H)^p. In this case, f_effecpress_ocean_p = (1 - Hf/H)^p (where Hf = flotation thickness). With a timescale > 0, we take (1 - Hf/H)^p as a target fraction, but relax toward that value over a period of ocean_p_timescale. As a result, N does not immediately drop to a low value in regions where Hf is close to H. This gives the ice time to thicken, so that Hf/H becomes smaller and the target fraction is larger. Thus, f_effecpress_ocean_p may never reach the low value that was the initial target. There was already a field called f_effecpress that acts in a similar way by phasing in the contribution of bwatflx to reducing N. This field is now call f_effecpress_bwat to avoid confusion with f_effecpress_ocean_p. In new Antarctic simulations, I set basal_physics%ocean_p_timescale = 500 yr, the same value as inversion%babc_timescale. This change inhibits retreat of the E. Thwaites GL (there is still some thinning, but much less than before) and generally reduces thickness biases for lightly grounded ice. In making this change, I reverted to an earlier method for compounding reductions in N. This method is multiplicative. Thus, if both basal water and ocean_p act to reduce N, they act in concert; we can multiply N by f_effecpress_bwat and again by f_effecpress_ocean_p. The more recent method was to choose the minimum of (1) N as reduced by basal water and (2) N as reduced by ocean_p > 0. Going back to the earlier method allows lower N and reduces the need for low C_c in some regions. Note: Setting ocean_p_timescale = 0 (the default) results in roundoff-level answer changes compared to previous code, because some operations have been reordered. Other changes: - Updated the ho_whicheffecpress descriptions in glide_setup. - Modified the interface for glide_define_restart_variables to pass in 'model' instead of 'model%options'. This allows restart variables to be added based on whether certain parameters (e.g., p_ocean_penetration) are zero or nonzero. - Changed a comment in glissade_bmlt_float
This commit adds a 2D scale factor called 'area_factor', to the glacier derived type, along with an option called 'scale_area'. When scale_area = .true., the area_factor is computed in each grid cell as cos^2(theta), where theta is the latitude. When scale_area = .false. (the default), area_factor is set to 1.0 everywhere. To use this option, simply add 'scale_area = .true.' to the [glaciers] section of the config file, and make sure that 'lat' (in degrees) is present in the input file. This commit is BFB except for diagnostic output. The glacier areas and volumes in the diagnostic log file are now smaller and should agree better with the true areas.
Subroutine glissade_glacier_update used to be responsible for inversion only, but now does other updates including glacier advance/retreat and an update of the smb_glacier_id mask, which determines where the SMB is applied during the following year. These updates require annual-average SMB and related fields, whether or not we are doing inversion. I found that some of these fields were being computed only when inversion was turned on. I changed the logic so that these fields are also computed with inversion off. This commit changes many lines of code, but most changes are simply changes in indentation, with some calculations moved outside of 'if inversion' loops. This commit is BFB for runs with inversion. The code now seems to be working correctly for runs without inversion.
This commit turns on verbose output for subroutine glide_forcing_read_once. It can take a long time to read in 240 time slices of GlacierMIP forcing (nearly 30 minutes for the full Alps at 200 m) and copy to local arrays. Logging progress to the cism output file allows the user to verify that the code isn't hanging.
Previously, the SMB for glacier cells was computed for each dynamic timestep based on the air temperature, precip (or snowfall), and surface elevation at that time. With this commit, the SMB is computed at the end of the year based on annual averages of air temperature and precip (or snowfall). The SMB is then applied uniformly during the following year. The differences are greatest for cells in the ablation zone. In these cells, it was typical to have some accumulation at the start of the year, then to melt all the ice during the summer (often having some melt potential left over), than add a bit of accumulation at the end of the year. These cells will now be ice-free at the end of the year, assuming the annual average SMB is negative enough to remove the advective inflow. I added model%climate%smb to the restart file for glacier runs, to preserve exact restart. In addition, I simplified the SMB masks: - I replaced smb_glacier_mask_init with cism_glacier_mask_init, which is held constant during the run. - I set smb_glacier_mask to 0 for cells with cism_glacier_id_init = cism_glacier_id = 0, instead of setting the mask to 1 in cells downstream of active cells. The goal is to have less flickering of ice thickness near the terminus. However, there is still some flickering. It might help if we go back to allowing melting in cells downstream of active cells. To be studied further. Also, I introduced an optional config parameter called precip_lapse. If precip_lapse > 0, then the precip rate increases in proportion to the difference between the ice surface elevation and the reference elevation. Huss and Hock (2015) introduced such a lapse rate, with values of 1.0 to 2.5e-4 (in units of fractional change per meter). The CISM default value is 0.0, but later we can test nonzero values, which will tend to reduce alpha_snow. This commit is answer-changing for all glacier runs.
This commit adds an elevation correction to the inversion for mu and alpha
in glacier calculations. The general effect is to reduce mu and alpha.
Recall the two equations solved when inverting for each glacier:
0 = alpha * snow - mu * Tpos,
smb_aux = alpha * snow_aux - mu * Tpos_aux,
where smb_aux, snow_aux, and Tpos_aux are the mass balance, snowfall, and
temperature surplus for the auxiliary climate (typically the climate
of the past two decades, for which we have geodetic mass balance estimates).
Both snow/snow_aux and Tpos/Tpos_aux are glacier-area averages.
We define Tpos = max(artm - Tmlt, 0) and Tpos_aux = max(artm_aux - Tmlt, 0).
Here, armt_ref_aux is computed from artm_aux by a lapse-rate correction,
and snow_aux is usually computed from precip_aux using a temperature threshold.
Suppose a glacier has been thinning at a rate of ~1 m/yr for the past two decades.
The direct cause is climate warming: artm increases at the reference elevation.
There is also an SMB-elevation feedback that grows over time. As the glacier thins,
the surface is lower and therefore warmer than it would be otherwise.
With the latest code changes, artm_aux is computed as follows:
artm_aux = artm_ref_aux + (usrf_aux - usrf_ref_aux) * T_lapse,
where usrf_aux, the effective surface elevation, is given by
usrf_aux = usrf + delta_usrf,
delta_usrf = (smb_aux - smb)*(rhow/rhoi)/1000 * dt_aux.
Here, delta_usrf is interpreted as the change in surface elevation during the transition
between the baseline climate and the auxiliary climate, due to the (usually negative)
SMB anomaly, and dt_aux is the length of the transition period.
(More precisely, dt_aux is twice the transition period, if the changes are linear and
delta_usrf represents the surface elevation halfway through the transition.)
Thus, in a warming climate, artm_aux is warmer than artm for two reasons:
(1) warming of the climate at the reference elevation
(2) warming of the ice surface due to loss of elevation.
Effect (2) is large enough to matter on decadal time scales.
If we ignore it (as we've done until now), we will estimate artm_aux to be too cool
and therefore mu to be too large. Including it, we lower mu (and alpha) and get a
lower sensitivity to a warming climate.
I added a new 2d field, smb_aux, that is written to the restart file,
and confirmed exact restart.
Another change: After changing the construction of SMB masks in the previous commit,
I reverted to the earlier treatment. With this treatment, the cells that can receive
a nonzero SMB include cells with cism_glacier_id_init = cism_glacier_id = 0,
provided these cells are just downstream of glaciated cells and have SMB < 0.
Removing these cells from the SMB mask resulted in excessive melting
(since mu must be larger if downstream cells with SMB < 0 aren't in the mask).
I ran a full-Alps commitment experiment with the changes.
The committed losses are still very high; more changes to follow.
…beta This commit includes several major glacier changes. First, I set up the inversion to distinguish between three different dates: * baseline_date = nominal date of the spin-up, when the glacier is in balance with the climate * rgi_date = date of RGI observations (ice extent and thickness), typically around 2003, when most glaciers were already out of balance * smbobs_date = central date of the SMB observations, currently Hugonnet et al. (2021), which provide the SMB term in the second inversion equation The inversion scheme tries to compute mu and alpha to give SMB = 0 at the baseline date and SMB = smb_obs for the recent climate. With this commit, CISM computes the SMB at both the baseline and smbobs date, and interpolates to get the SMB at the RGI date. This SMB can be averaged between the baseline and RGI dates to estimate the thickness change between these two dates. This thickness change is then used to correct the baseline thickness target. For most glaciers, the baseline thickness target exceeds the RGI thickness estimates. Thus, the spun-up baseline glaciers are thicker than the RGI glaciers. The model can then be run forward to the RGI date to a state that is a good match to the RGI thicknesses, and is thinner than the baseline state. This distinction is important for GlacierMIP3, which stipulates that models should be initialized to the RGI date, with ice losses are computed relative to this date. Initializing to an earlier date will overestimate the losses. Next, I streamlined the tuning procedure for glacier-specific parameters mu_star, alpha_snow, and beta_artm. These changes were motivated by the fact that for many glaciers, especially the smallest ones, the Hugonnet SMB estimates have large errors. So it is not worth taking extraordinary measures to match Hugonnet for all glaciers. Briefly, the procedure is now as follows: For glaciers with smb_obs of the right sign (i.e. smb_obs < 0, with D > 0), solve the 2-equation system. * If mu and alpha are in range, we keep them. * If not, we prescribe alpha within range and solve the 1-equation system. If mu is in range, we keep it. * If not, we increment beta, making the air temperature warmer or cooler (uniformly for the baseline and recent climate). We keep adjusting beta until mu is in range. For glaciers with smb_obs of the wrong sign, we ignore smb_obs and solve the 1-equation system as above. For glaciers with little or no baseline melting, we increase beta to induce some melting. I added some diagnostics to count and write out the glaciers that violate Eqs. 1 and/or 2. I removed the parameter beta_artm_aux, which is no longer used. I modified the criterion for computing smb_glacier_id in cells that border glaciated cells. Diagonal neighbors are now included, not just edge neighbors. This increases the baseline glacier volume, since there is a greater area with SMB < 0 that must be balanced by cells with SMB > 0. I removed subroutines reset_glacier_fields, accumulate_glacier_fields, and average_glacier fields. This code in now inlined. For simplicity, I removed the option to set an elevation lapse rate for precipitation. I changed some of the default glacier parameter options. In particular, the new default Tmlt = -1 C. I verified exact restart with the new inversion scheme. With these changes, the baseline area and volume for the full Alps are 10 to 15% greater than the RGI area and volume. This is with a baseline date of 1984, using a 1979-1988 climatology.
Inversion for glacier parameters requires climate data for two periods: a baseline period when glaciers are assumed to be in balance with the climate, and a recent period when glaciers are out of balance and we have SMB observations. We've been reading artm_ref, snow, and precip for these periods from two different forcing files. For the recent period, we've read in 'auxiliary' fields artm_ref_aux, snow_aux, precip_aux. With this commit, instead of reading in auxiliary fields, we read anomaly fields artm_ref_anomaly, snow_anomaly, and precip_anomaly. The reference air temperature for the recent period is given by artm_ref_recent = artm_ref + artm_ref_anomaly, and likewise for snow and precip. These could be read from separate files, but for now I put the baseline and anomaly fields in a single file. With this change, the answers are the same within rounding error, but not BFB. The dates for the baseline and recent climates, along with the RGI reference date, are now config parameters. I added the anomaly fields in glide_types, removed the aux fields, and changed some variable names for clarity and consistency. I also confirmed exact restart.
This commit fixes some minor issues from recent glacier runs:
(1) Fixed some log messages that go along with reading read_once files.
The log files were indicating that certain time slices were being read erroneously,
but in fact the files were read correctly.
(2) For the glacier diagnostics, I added a count of the number of glaciers with
nonzero area and volume.
In glacier runs, ice can form in a grid cell that is adjacent to two different glaciers. We then have to decide which glacier this grid cell belongs to. The old criterion was to choose the glacier with the more negative SMB. This prevents pirating of glaciers with high melting by glaciers with low melting. This commit introduces a new criterion based on the input ice fluxes. Given the thickness and velocity fields, CISM computes the ice volume flux across each of the four edges of a newly advanced glacier cell. These edges fluxes are computed in a new subroutine in glissade_utils.F90. If there are incoming fluxes from two different glaciers, the cell is assigned to the glacier providing the largest flux across an edge. This change prevents the Glacier des Bossons (cism_glacier_id = 3481) from advancing across the mouth of neighboring glacier (3482; I think it's called the Glacier de Taconnaz). However, it still allows the Bossons glacier to turn left and deliver ice to Taconnaz, resulting in unrealistic advance of Taconnaz and retreat of Bossons. I then introduced a dynamic change for glaciers. Based on cism_glacier_id_init, we compute a boundary mask at initialization. When a cell edge has different glaciers on each side, we set the mask to 1, which forces powerlaw_c to be held at its maximum value at the two vertices of the edge. This minimizes sliding (though internal deformation is still allowed), reducing flow between glaciers. This change reduces, but does not eliminate, the spurious flow from Bossons to Taconnaz. To further reduce this flow, we would probably need to resolve the topography better. In reality, a narrow ridge separates the two glaciers, preventing flow from one to the other.
In subroutine glissade_glacier_update, inversion-related calculations are now done only when actually inverting for mu_star, alpha_snow, and/or powerlaw_c. This prevents extraneous calculations in forward runs without inversion. I also reduced and cleaned up some diagnostic output. This commit is BFB for runs with inversion.
This commit adds some options to support anomaly forcing in forward glacier runs.
It is now fairly straightforward to do a historical run starting at the baseline date
and extending to the RGI date or recent date, with anomaly forcing ramped up linearly.
The anomaly fields for glaciers are artm_anomaly, snow_anomaly, and precip_anomaly.
(I renamed artm_ref_anomaly to artm_anomaly.) These fields can now be used in two ways:
(1) When inverting for mu_star and alpha, the anomaly fields are added to the baseline
climate fields to obtain values for the recent climate, which in turn are used
to compute the SMB for the recent climate.
(2) In forward runs, the anomaly fields are read in at initialization and then
can be ramped up over some timescale. We sometimes do this in ISMIP6 forward runs.
In case(1), we should have enable_artm_anomaly = enable_snow_anomaly =
enable_precip_anomaly = .false. (the default values).
This is because the anomalies are used only for inversion; they are not
part of the baseline climate.
In case (2), the user should set enable_artm_anomaly = enable_snow_anomaly =
enable_precip_anomaly = .true. in the config file.
Then the anomalies are added to the baseline fields (artm, snow, and precip) in glissade.F90.
To make the forcing more flexible, I added a config variable called artm_anomaly_tstart.
This is the time (in years) when we begin applying the anomaly.
The default is year 0, which previously was the only option.
I also changed the anomaly routines to increase the anomaly at each timestep
during the ramp-up period. The old default was to increase the anomaly only once per year,
following ISMIP6 protocols. The ISMIP6 behavior can be recreated by uncommenting one line.
This changes the answers slightly.
Spin-up answers also change slightly, because usrf is accessed earlier in the timestep
for the lapse-rate correction to artm.
Following a glacier spin-up, I did a historical run from the baseline date (1984)
to the RGI date (2003). The Alps lose about 15 km^3 of ice by the RGI date, which is
still about 6 km^3 above the RGI target value, even though the SMB values are close to
the values used to set the baseline targets. This might be because of slower flow,
or nonlinear decrease of the SMB with rising temperatures.
This commit contains several minor changes to allow the code to compile
(and do so efficiently) on the Derecho Intel compiler.
Notably, changes in generate_ncvars.py and ncdf_template.in will result
in many fewer 'use glide_types' and other use statements that appeared
in many subroutines of glide_io.F90. Now, these use statements appear
only at the top of the module. As a result, CISM compiles in just over
a minute on 8 cores ('make -j 8), compared to more than 4 minutes before.
(For some reason, this wasn't a problem on the gnu compiler.)
Also added a missing use statement (use glide_stop) in glide_initialise.F90
and a missing include statement (#include <ctype.h>) in writestats.c.
When the glacier branch is rebased to main, these changes will already have
been done on main, possibly leading to minor conflicts, but these
shouldn't be hard to resolve.
Our glacier grids are different from typical ice sheet grids in that the nominal resolution (e.g., dx = 200 m) is coarser than the true cell dimensions, which are given by dx * cos(latitude). For grid cells in the Alps, which lie near 45 N, a typical grid cell on a 200-m grid represents a region whose dimensions are roughly 140 x 140 m. This raises a couple of issues. First, to diagnose the true area of a grid cell, we need to scale the nominal area (say, 40000 m^2) by cos^2(lat). There was already some code to handle this in glacier area computations. With this commit, the adjustment is applied consistently across the code. We define a 2D field called cell_area, which can vary from cell to cell. By default, cell_area = dew*dns (where dew and dns equal the nominal resolution). When glaciers are enabled and scale_area = .true., cell_area(i,j) is multiplied by cos^2(lat(i,j)) for each cell. This value of cell_area is used only for diagnostics, not in the ice dynamics. Second, the gravitational driving force and internal ice stresses depend on the distance (dew or dns) between cell centers. To get these forces correct, we should use the true rather than nomimal dimensions. This commit introduces a new config option, length_scale_factor, that can be used to modify dew and dns. Since dew and dns are scalars (not 2D fields), the same scale factor is applied everywhere. We should choose a factor that corresponds to the average latitude in a region. For instance, we could set length_scale_factor = sqrt(2)/2 ~ 0.707 for a region at latitude 45 N. The default value is length_scale_factor = 1. When length_scale_factor /= 1, dew and dns are modified at initialization. This changes answers throughout the code. For now, the length scaling can be applied only when glaciers are enabled. We could potentially make it more general.
This commit changes the treatment of ice-free peripheral cells for purposes of SMB inversion and weighting. Here, 'ice-free' means ice-free at the end of the timestep, after applying a negative SMB. But these cells often have ice after advection, before applying the SMB. Recall that the two equations used for inversion are both summed over the initial glacier area. This includes all cells with cism_glacier_id_init > 0. The question is how to deal with adjacent cells where cism_glacier_id_init = 0, if ice can flow into these cells and melt. One choice is to ignore melting in these peripheral cells when doing the inversion. This typically leads to higher values of mu_star, since the estimated ablation is assumed to occur in fewer cells than actually have ice loss. The spun-up glaciers tend to be too small. Another choice is to include all such peripheral cells when doing the inversion. This typically leads to lower values of mu_star, since the estimated ablation has more cells to work with. The spun-up glaciers tend to be too big. A compromise is to assign these cells a weight between 0 and 1 when summing over glaciers for Eqs. 1 and 2. The weight w is computed as the ratio between the applied SMB and the potential SMB. For instance, an ice-free cell near a glacier terminus might have a (potential) computed SMB of -5 m/yr. But suppose the applied SMB is -2 m/yr, because it takes only 40% of the potential SMB to melt all the ice advected into the cell. Then we assign the cell a weight w = 0.4 when computing all-glacier sums. We can think of the cell as a partial cell that is only 40% ice-covered. This change has the desired effect. The spun-up glaciers, on average, have areas and volumes closer to their targets. Other changes: * Removed the field usrf_target_rgi. The target surface elevation field at the RGI date is now usrf_obs, the observed surface elevation. * Added a subroutine to compute the min and max SMB for each glacier * Added code to count the number of cells included in each mask for each glacier * Added code to sum area and volume over just the initial extent of each glacier, not including advanced cells * At startup, set glacier thickness targets in ice-filled cells to at least the dynamic minimum thickness * Added parallel_reduce_max and parallel_reduce_min subroutines for 1D arrays * Streamlined some diagnostics
The usual way of spinning up glaciers is to solve two equations for two tunable
parameters, mu_star and alpha_snow.
An earlier scheme, which is still supported, is to solve one equation for mu_star:
SMB = alpha_snow * snow - mu_star * Tpos,
where we assume that snow and Tpos are from a balanced climate, hence SMB = 0 and
mu_star = alpha_snow * snow / Tpos.
Instead of inverting for alpha_snow, we set alpha_snow to a prescribed constant.
This commit updates the 1-equation scheme to more closely follow the logic of the
2-equation scheme. For example, weights between 0 and 1 are applied to ice-free cells
adjacent to ice-covered cells at the glacier periphery. Also, the temperature parameter
beta_artm is now adjusted as needed to bring mu_star into a prescribed range.
For the Alps, adjusting beta_artm (with a max adjustment of 5 C) brings mu_star
into range for all but a handful of glaciers within 100 years.
The logic that computes RGI and recent climate SMBs during the inversion is now applied
only when inverting for both mu_star and alpha_snow, not for mu_star alone.
This commit is answer-changing for the 1-equation scheme but not the 2-equation scheme.
Until now, the config option 'restart' has had two possible values:
* restart = 0 if not a restart (i.e., read in a CF input file and initialize the ice state)
* restart = 1 if a restart (i.e., read in a CF restart file that includes the full ice state)
This commit adds a new option, restart = 2, called a 'hybrid restart' because of its
similarity to a CESM hybrid run. We now refer to option 1 as a 'standard restart'.
A hybrid restart works as follows:
- The run is initialized from a file in the [CF input] section, as for restart = 0.
However, this file has 'restart' or '.r.' in its name and includes the full ice state
as needed for exact restart. Typically, it is the final restart time slice from a spin-up.
- The initial time (model%numerics%time) is set to 'tstart' as specified in the config file.
This differs from a standard restart, which takes its initial time from the restart file.
- The initial tstep_count = 0.
This differs from a standard restart, which takes tstep_count from the restart file.
For glaciers, we can use the hybrid restart option for commitment runs and other
forward runs starting from a spun-up state. For instance, say we want to do a 2000-year
spin-up followed by a forward run from 2003–2100. The workflow is as follows:
- Do the spin-up and write a final restart file.
- Set up a directory for the forward run with the required input and forcing files,
including the restart file from the spin-up.
- In the config file:
* Set tstart = 2003, tend = 2100., and restart = 2.
* In the [CF input] section, set 'name' to the name of the restart file from the spin-up.
This should be different from the name of the [CF restart] file for the forward run.
E.g., one file could be spinup.restart.nc, and the other could be forward.restart.nc.
* Change other options changes as needed, e.g. change the inversion options.
- Launch the forward run.
It is no longer necessary to modify the model time or tstep_count by hand in the restart file
from the spin-up. Also, it is not necessary to use the same name for (1) the restart file
from the spin-up and (2) the restart file for the forward run.
If the [CF output] file for the forward run is configured with 'write_init = .true.',
then this output file will include the ice state at the start of the forward run.
For a standard restart, CISM does not write to the output file at the start of the forward run.
Testing the new restart option, I confirmed that a hybrid restart is exact, as expected,
apart from the new values of the model time and tstep_count.
I also confirmed that the standard restart (restart = 1) works as before.
This commit adds a model option called forcewrite_final. If true, the model is forced to write to each netCDF output file, including restart files, when a run finishes. This can be useful if we want to write output at regular intervals (e.g., frequency = 50 years) and also write output at the end of the run, even when the total number of years (tstart - tend) is not divisible by the frequency. For other runs (e.g., a series of short debugging runs), we might not want to write output when finishing each run. For this reason, the default is forcewrite_final = .false. To override the default, set forcewrite_final = .true. in the config file. Even if model%options%forcewrite_final = .false., it remains possible to force a final write by calling subroutine glide_finalise with an optional argument set to '.true.'. This argument used to be called 'crash'; now it is called 'forcewrite_arg', since it might be desirable to write final output regardless of whether the model has crashed. Note: As before, output is not written to netCDF when the model aborts with a fatal error. In that case, subroutine parallel_stop calls mpi_abort, which aborts the run without writing output.
This commit includes several tweaks to the strategies for limiting glacier advance and retreat.
I added a logical glacier option: 'redistributed_advanced_ice'.
The default is .false. for backward compatibility.
When the option is set to .true., advanced ice in the accumulation zone is thinned
at a prescribed rate ('thinning_rate_advanced_ice', a new config parameter).
The ice mass removed is redistributed uniformly across the initial extent of the glacier.
In runs with a thinning rate of 2 m/yr, this change reduces but does not eliminate advanced ice.
A higher thinning rate removes more ice, but with diminishing returns;
there is a tradeoff between a correct glacier margin and artificial high thinning rates.
I modified the way glacier IDs are assigned to advanced cells in subroutine glacier_advance_retreat.
Occasionally, such a cell is adjacent to two or more neighboring glaciers.
Previously, it was given the ID of the glacier contributing the greatest flux.
This leads to difficulties when no glacier contributes a positive flux.
In the new code, the cell is assigned the neighbor ID that results in the most negative SMB.
This is similar to the way peripheral cells in the ablation zone are assigned smb_glacier_id.
I introduced a new config parameter, smb_weight_advanced_ice, in the range [0,1].
This is the weight given in the inversion calculation to glacier-free cells in the ablation zone,
where ice can be advected and melt without ever being thick enough to glaciate the cell.
Trial and error has shown that w = 0 tends to drive high mu_star and spurious retreat,
whereas w = 1 results in low mu_star and spurious advance.
I tried setting w based on the ratio of applied to potential SMB, but this gives w = 0
in retreated cells, which encourages further retreat.
A value w = 0.5 seems a good compromise; this corresponds to the case that half the potential SMB
is used to melt ice and the other half is not used. This value is the default for now.
I modified the computation of smb_glacier_id. Previously, this ID was set to 0 in the accumulation zone.
Now, every glacier-free cell adjacent to a glacier cell is given smb_glacier_id > 0,
but any positive SMB is set to zero. The result is the same, but the logic is clearer.
I removed the deprecated subroutine 'remove_snowfields'.
In glissade_utils, I added a subroutine to estimate the input ice flux to each cell from each neighbor.
I decided not to use this subroutine for glaciers, but I left it in case it's useful.
This commit adds some diagnostic scalars and fields related to thickness inversion: * thck_target, a 2D field. This is the thickness target for inversion; it can now be written to output files. It is not needed in the restart file. * glacier%area_target and glacier%volume target. These are targets for each glacier, obtained by summing cell_area and thck_target over the glacier. * tot_glc_area_target and tot_glc_volume_target. These are computed by summing over all glaciers before writing to the diagnostic log file. * rmse_thck and rmse_thck_init_extent. These are root-mean-square errors of (thck - thck_target). With these new scalars and output fields, it is easier to determine which parameter settings come closest to matching the targets.
This commit adds three scalar glacier diagnostics: * glacier_total_area = area summed over all glaciers * glacier_total_volume = volume summed over all glaciers * nglacier_active = number of active glaciers (i.e., glaciers with nonzero area) Each is now part of the glacier derived type, and each can be added to the appropriate variable list in the config file. Within the code, total_area has units of m^2 and total_volume has units of m^3. However, the netcdf variables have scale factors to convert to km^2 and km^3.
This commit changes the inversion algorithm for Cp in advanced glacier cells (i.e., cells that are initially ice-free but become glaciated). The Cp evolution equation has three terms: * a term proportional to the thickness difference from the target * a term proportional to dH/dt * a relaxation term that nudges Cp toward Cp_const For advanced cells, the dH/dt term is now ignored. This means that Cp decreases smoothly toward Cp_min, without oscillations. I made this change to improve stability for the Lower Grindelwald Glacier at 100-m resolution. Larger Cp_const, Cp_min, and babc_relax_factor also improve stability.
This commit cleans up a rebase consistency issue. Subroutines glide_finalise_all and glide_finalise now have an optional argument called 'finalise_arg' rather than 'crash_arg'.
This commit cleans up a rebase issue. I added the 'model_id' argument in several calls to glide_add_to_restart_variable_list.
I added a note to the README file for the slab stability test.
This test runs the slab problem at multiple spatial resolutions and finds the
maximum stable timestep at each resolution.
The test suggested in the README file is this one:
python stabilitySlab.py -n 4 -a DIVA -theta 0.0375 -thk 1000. -mu 1.e5 -beta 1000. \
-dh 0.1 -nt 100 -nr 12 -rmin 10. -rmax 40000.
This test fails with an energy conservation error.
However, energy conservation is not really violated; we just have a poor choice
for the conservation threshold for this problem. Here is the fix:
Comment out this line in glissade_therm.F90:
if (abs((efinal-einit-delta_e)/dttem) > 1.0d-7) then
Uncomment this line:
if (abs((efinal-einit-delta_e)/(efinal)) > 1.0d-8) then
The README file now includes instructions for the fix.
|
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. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
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:
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.