Katetc/hma glaciers4, Replace Lipscomb/hma_glaciers4 PR #64
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
I added module glissade_glacier.F90 to support simulations with glacier regions. I added the glide_glacier derived type, an enable_glaciers config option, and several glacier I/O fields. The glacier module is still under construction. This version contains subroutine glissade_glacier_init, which performs the following tasks: * Read in the 2D array glacier_id, which associates each grid cell with a unique glacier ID, typically a Randolph Glacier Inventory (RGI) ID. * Count the number of glaciated grid cells. * Sort the glaciers in order of ascending ID, using a quicksort algorithm. * Count the number of glaciers (nglacier). * Create an array that maps CISM-specific glacier indices (1:nglacier) to the RGI indices. * Create a 2D array called glacier_id_cism, which is like glacier_id except that it associates each cell with a CISM-specific glacier index. * Allocate some other CISM-specific arrays of size(nglacier). * Compute the initial area and volume of each glacier. These will serve as inversion targets. I tested the initialization subroutine on an Everest-region 100m grid of size 1411x1061, with about 1200 glaciers and 150,000 glaciated cells. The logic seems to be working.
This commit adds infrastructure for glacier simulations, including:
* subroutine glissade_glacier_smb, which computes the SMB using the formula
SMB = snow - mu_star * max(artm - Tmlt, 0)
where snow = solid precip rate
artm = surface air temperature
Tmlt = temperature threshold for ablation (set to -1 C)
mu_star = tunable parameter with units of mm/yr w.e./deg
This formula is based on the SMB equation in OGGM (Maussion et al., 2019).
The subroutine is called from glissade_tstep during the SMB calculations.
I added a 'snow' array to the 'climate' derived type and glide_vars.def.
* subroutines glissade_glacier_inversion, glacier_invert_powerlaw_c, and
glacier_invert_mu_star.
These subroutines invert for mu_star based on a glacier area target,
and for powerlaw_c based on a glacier volume target.
Subroutine glissade_glacier_inversion is the inversion driver;
it is called from glissade_diagnostic_variable_solve.
Inversion for glacier-specific powerlaw_c is similar to inversion for
vertex-based powerlaw_c. There are two terms: one proportional to (V - V_target)
and one proportional to dV/dt. The powerlaw_c inversion seems to be working.
Inversion for mu_star uses just one term, proportional to (A - A_target).
This is because glacier area change will be discontinuous (i.e., the area changes
in increments of one gridcell), so dA/dt is not well defined.
I have not tested mu_star inversion, since glacier areas are still fixed.
* I added some new variables to the glacier derived type.
This type now includes 9 arrays of dimension(nglacier):
glacierid, cism_to_rgi_glacier_id, area, volume, area_target, volume_target,
dvolume_dt, mu_star, and powerlaw_c. Several were added to glide_vars.def.
In addition, there are two glacier arrays of size(ewn,nsn): rgi_glacier_id and
cism_glacier_id. The former is an input array using an RGI integer ID;
the latter is numbered from 1 to nglacier.
* I added config parameters mu_star_const, mu_star_min and mu_star_max.
The default values are likely to change.
* I added a netCDF dimension called glacierid with dimension(nglacier),
where nglacier is the total number of glaciers.
This allows CISM to read and write 1D output arrays with dimension(nglacier).
* I modified generate_ncvars.py to use parallel_put_var instead of
distributed_put_var to write 1D arrays such as the new glacier arrays.
I added an interface, parallel_put_var_integer_1d, in the cism_parallel modules.
* I moved the glissade_glacier_init call before the call to glide_io_createall,
so that nglacier and glacierid are correct when output files are created.
* In glide_diagnostics.F90, mass is now given in units of Gt when dm_dt_diag = 1.
Note that the loops from 1 to nglacier are done on all processors,
with identical answers on each processor. This will not scale well.
Since many glaciers straddle multiple processors, it is not obvious
how to distribute the glaciers across all processors.
I ran the model on the Everest grid for 10 years without obvious errors,
with CISM cycling through a climatological forcing file with 12 monthly slices.
Output fields look correct.
This commit includes changes to ensure exact restart when inverting for glacier properties. The following fields in the glacier derived type are written to the restart file: * rgi_glacier_id * cism_glacier_id * cism_to_rgi_glacier_id * glacier_mu_star * glacier_powerlaw_c * glacier_area_target * glacier_volume_target The first two fields are defined on the horizontal grid; the others have dimension(nglacier). Area and volume targets are needed for exact restart when inverting for mu_star and powerlaw_c, but not for forward runs. The following are recomputed on restart and are not needed in the restart file: nglacier, ngdiag, glacierid, glacier_area, glacier_volume. On restart, we need to know nglacier to read 1D glacier arrays from the restart file. We could add nglacier to the config file when restarting, but until now we have avoided making config files different for restart compared to start-up (apart from changing 'tend' and setting restart = 1). Instead, I added a call to a new subroutine, glimmer_nc_get_dimlength, in glimmer_ncio.F90. This subroutine parses the restart file for the length of dimension 'glacierid'. On the Everest grid, I confirmed exact restart for a short run with inversion. Note: I am running with global_bc = 3 (no ice in cells adjacent to the global boundary). With periodic or outflow BC, restart is not exact because there are glaciers along the global boundary on the Everest grid. Since the velocity grid is smaller than the ice grid, velocities are not correct when ice is present at the left and lower boundaries of the global domain.
This commit includes the following changes for glacier calculations:
* Modified the inversion for mu_star. In the previous version, there was a prognostic equation
for mu_star based on the difference between the current and target glacier areas.
This could be slow to converge. Since SMB is a linear function of snowfall and temperature,
I implemented a more direct method, computing the value of mu_star that gives SMB = 0
over the initial area of each glacier. This method uses the cism_glacier_id_init array
instead of the area_target array, which now is purely diagnostic.
In general, mu_star will be set to a value that prevents large advance or retreat.
An exception would be if a glacier has no ablation zone; i.e., the monthly mean air temperature
never exceeds Tmlt. I set Tmlt = -2 C (instead of -1 C, as suggested by Maussion et al. 2018)
to make it less likely that glaciers lack ablation zones.
The method would need to be modified for marine-terminating glaciers.
* Introduced a parameter inversion_time_interval, which controls how often the inversion
calculation is called. The interval must be an integer number of years, with default value 1 yr.
Inverting on shorter timescales would introduce unnecessary sub-annual variations.
Three fields - snow, Tpos, and dthck_dt - are accumulated and averaged over this interval
to support the inversion. Here, Tpos = max(artm - Tmlt, 0.0).
* Added a subroutine, glissade_glacier_advance_retreat, to support re-indexing
as glaciers advance and retreat. The rules are as follows:
- At start-up, glaciated cells have cism_glacier_id in the range [1, nglacier].
Other cells have cism_glacier_id = 0. The initial indices are saved as cism_glacier_id_init.
- When a cell has H < H_min and cism_glacier_id > 0, we set cism_glacier_id = 0.
It no longer contributes to glacier area or volume.
Here, H_min (= 5 m by default) is a threshold for counting ice as part of a glacier.
- When a cell has H >= H_min and cism_glacier_id = 0, we give it a nonzero ID:
either (1) cism_glacier_id_init, if the initial ID > 0,
or (2) the ID of an adjacent glaciated neighbor (the neighbor with
the highest surface elevation, if there is more than one).
Preference is given to (1), to preserve the original glacier outlines if possible.
- If H >= H_min in a cell with cism_glacier_id_init = 0 and no glaciated neighbors,
we do not give it a glacier ID. Instead, we set H = H_min and remove the excess ice.
There is no glacier inception; we only allow existing glaciers to advance.
* Put the tunable glacier parameters at the top of module glissade_glacier.
If desired, these could be added to the glacier derived type and set in the config file.
This commit adds a section called 'glacier' in the config file. Currently, two glacier options and two parameters can be set in this section: * set_mu_star (formerly glacier_mu_star in the 'ho_options' section) * set_powerlaw_c (formerly glacier_powerlaw_c in the 'ho_options' section) * minthck (min ice thickness counted as a glacier) * tmlt (min air temp at which ablation occurs) Later, I would like to group other sets of physics options and parameters in their own sections: e.g., 'calving', 'basal_physics'. I also removed some old, commented-out basal process options from glide_types.
When glaciers are enabled, CISM now writes some global glacier diagnostics (number of glaciers, total area and area target, total volume and volume target) and single-glacier diagnostics (area and area target, volume and volume target, mu_star and powerlaw_c for glacier ngdiag) to the log file.
This commit includes several fixes: * If constant (HO_POWERLAW_C_CONSTANT, HO_COULOMB_C_CONSTANT), powerlaw_c and coulomb_c are now initialized once at the start of the run instead of repeatedly in calcbeta. The previous logic was overriding the glacier-derived powerlaw_c in calcbeta. As a result, which_ho_powerlaw_c is no longer passed to calcbeta. We still pass which_ho_coulomb_c to handle the case of elevation-dependent coulomb_c. * Subroutine glissade_glacier_advance_retreat now has logic to update acab_applied in grid cells where the ice thickness is limited to block glacier inception. This subroutine is now called during each timestep instead of once a year. In the future, we might want to call it once a year (to save calculation), but then we might also want to keep track of the ice removal in a separate budget, to avoid having large negative acab corrections during one timestep per year. * I added a short subroutine, reset_glacier_fields, to set the accumulated fields back to zero after each glacier inversion calculation.
With this commit, glacier%minthck is no longer a config parameter. Instead, it is set to a value slightly less than model%numerics%thklim, which determines the threshold of dynamically active ice and usually is set to 1 m. Recall that any cell with a nonzero glacier ID is set to ng = 0 if it becomes thinner than glacier%minthck. Any cell with ng = 0 and H > glacier%minthck receives an ID > 0 if it is part of an initial glacier or adjacent to a cell with ng > 0; otherwise, H is set to glacier%minthck to keep it inactive. I added a subroutine, glacier_powerlaw_c_to_2d, that fills the 2D array model%basal_physics%powerlaw_c, given powerlaw_c for each glacier. CISM sets model%basal_physics%powerlaw_c = 0 at vertices that are not adjacent to any glacier cells. This setting could cause problems if non-glacier cells were dynamically active.
I changed the target for the volume inversion. Instead of trying to match the total glacier volume from the input file, CISM now tries to match the volume over the observed glacier footprint, i.e., the region covered by the initial glacier. This is a new variable in the glacier derived type, called volume_over_init_region. The goal is to avoid cancelling errors in the volume inversion. If a glacier advances, its area will increase, so if CISM is simply trying to match the total volume, it will tend to make the interior of the glacier too thin (and conversely for retreating glaciers). With the new method, CISM will generate glaciers that have the correct thickness over the observed footprint.
Until now, the surface air temperature artm has been read in as a 2D field and applied to glaciers as read. This commit allows the input artm to be corrected using a spatially uniform temperature lapse rate. For glacier runs, forcing files should now have fields called 'artm_ref' and 'usrf_ref', where usrf_ref is the reference surface elevation at which artm_ref is valid. Note: usrf_ref is a new name for what used to be called smb_reference_usrf. It has no scaling parameter (unlike usurf, whose scaling parameter thk0 will at some point be removed). The 'options' section of the config file should set artm_input_function = 3. This is a new option, ARTM_INPUT_FUNCTION_XY_LAPSE, which specifies that artm should be read in as a function of (x,y) and corrected with a uniform lapse rate. Cf. option 1, in which the correction is given by a 2D field, artm_gradz. The 'parameters' section of the config file should specify t_lapse, which I moved from the glacier type to the climate type. Its default value is 0. For HMA glacier runs, I am setting t_lapse = 0.005 degC/m. Also, I modified the volume inversion calculation to compute dthck_dt only over the initial glacier footprint defined by cism_glacier_id_init. I increased powerlaw_c_timescale from 10 yr to 25 yr. I added the 2D arrays snow_accum and Tpos_accum to glide_vars.def.
This commit includes a new method of inverting for powerlaw_c in glacier runs. Multi-century runs showed that it does not work well to have one value of C_p for an entire glacier. There are century-scale oscillations in C_p and mu_star in different parts of the glacier, especially large glaciers with a long residence time. For example, if a glacier advances too far, then mu_star will increase, which drives retreat and reduces the total volume (as well as the area), which drives C_p higher, which leads to thickening and re-advance. I wrote a new powerlaw_c inversion scheme that adjusts C_p in each grid cell based on the thickness bias (H - H_obs), using an equation similar to that for ice sheets. Thus, C_p for glaciers is now a 2D array. The main difference from the ice-sheet approach is that the single thickness scale (typically 100 m) is replaced by max(thck_obs, glacier_thck_scale), where glacier_thck_scale = 100 m. Thus, thickness errors are weighted more heavily for thin ice than for thick ice. To support this method, I added usrf_obs and powerlaw_c (the 2D field) to the restart file. In new multi-century simulations, C_p no longer has huge oscillations. Long glacier tongues are intact, and thickness biases are much reduced. I tried two values of T_mlt: -2 C and -4 C. The thickness biases in the two runs are similar. In the second run, values of mu_star are typically 50% or less of the values in the first run. I also tried a new inversion scheme for mu_star, prognosing mu_star using an equation of the same form as the equation for powerlaw_c, but with an area target instead of a volume target (and still with a single value per glacier). However, this method does not work well. In the observations, there are a number of ice-free grid cells in regions with SMB > 0, perhaps because of steep topography. Ice tends to advance into these cells, increasing the glacier area. As a result, mu_star increases, compensating for the area gain at high elevations with area loss at low elevations. Many glacier tongues are lost. The old scheme (setting mu_star so that SMB = 0 over the initial footprint) works better. For now, I left the new subroutine (glacier_invert_mu_star_alternate) in the code for reference. Minor changes: * Inserted a halo update for model%geometry%thck before the glacier initialization. This is a temporary hack that is needed because I am running on a reduced domain with glaciers at the global boundary. The no-ice global BCs automatically remove ice at the global boundary. The new halo update removes ice from these cells before creating the thickness targets, so that the inversion algorithm does not aim (in vain) for nonzero targets. * Changed mu_star_min and mu_star_max to 20 and 20000, respectively. * Changed the glacier_powerlaw_c_timescale to 100 yr, and inserted a thickness scale of 100 m. * Moved glissade_usrf_to thck and glissade_thck_to_usrf to the glissade_utils module.
With this commit, glacier advance and retreat (leading to glacier index changes) are computed at the end of each year, instead of every timestep. This prevents spurious winter advance and summer retreat associated with subannual thickness changes. Recall that the advance/retreat subroutine limits the ice thickness in non-glacier cells, and the limiting is treated as a negative contribution to acab. In the future, we might want to classify this limiting as part of a non-physical correction flux. I added a new output field, glacier_mu_star_2d, which is simply glacier_mu_star mapped onto the horizontal grid. I changed the reset timing for glacier Tpos_accum and snow_accum. These are now zeroed out at the start of a new year (after writing output) instead of at the end of the previous year (before writing output).
Exact restart was not working; needs glacier_mu_star.
This commit rearranges some calls in the glacier inversion subroutine, such that when running with glacier_mu_star and glacier_powerlaw_c inversion off, but reading these fields from external files, some standard glacier diagnostics (area and volume) are updated during the run.
This commit adds some functionality to the option enable_artm_anomaly. Previously, we set enable_artm_anomaly = T when reading a 2D artm_anomaly field from the restart file. With this commit, it is also possible to prescribe a spatially uniform anomaly, e.g. 1 deg C everywhere. To use the new option, simply leave out artm_anomaly from the input file, and set artm_anomaly_const in the [parameters] section of the config file. This will yield uniform warming of the desired value. While adding this option, I found a logic bug in glissade.F90 that prevented the correct application of both an elevation adjustment and an anomaly to artm. I rearranged some logic so that both adjustments can be applied in the same run. In GlacierMIP-style experiments, for example, we may want to prescribe a warming anomaly while also adjusting for elevation. This commit is not BFB. In changing the logic, I fixed a one-timestep lag in setting the ice surface temperature to artm.
This is a bug fix commit following a rebase to lipscomb/basal_physics3. I removed a 'public' statement for subroutine usrf_to_thck, which no longer is part of glissade_inversion.F90. I also changed some calls from usrf_to_thck to glissade_usrf_to_thck. Since I missed this during the rebase, some recent commits may not compile.
With this commit, CISM can read a 2D field consisting of observed glacier SMB, smb_obs.
The source data is assumed to be from Hugonnet et al. (2021) or a similar observational dataset.
It consists of the mean SMB per glacier over some period, typically one or two decades.
Before being read in, the per-glacier data must be mapped to the CISM grid, typically
with a single value over the entire area of a given glacier.
It would be possible to read in a 1D list of SMB values per glacier, but it is simpler for CISM
to work with 2D gridded values.
Once read in, SMB_obs can be used to compute mu_star for each glacier using the relation
SMB = P_s - mu_star * max(T - T_mlt, 0),
where P_s is solid precip, T is surface air temperature, and T_mlt is a temperature threshold for ablation,
with all values being monthly.
Summing this relation over each glacier and each month of the year, it is straightforward to find mu_star.
In earlier code versions, we assumed an equilibrium mass balance, i.e. SMB = 0 over the observed glacier area.
Using the Hugonnet data is likely to give a more accurate mu_star for present-day climate.
I also removed a deprecated subroutine, glacier_invert_mu_star_alternate.
With the previous commit, mu_star can be computed during inversion based on smb_obs,
the observed SMB from a dataset such as Hugonnet et al. (2021). In general, the resulting
SMB is not equal to 0; for most, SMB < 0.
This makes it tricky to run a long spin-up (e.g., long enough to invert for powerlaw_c)
while maintaining the present-day glacier footprint.
The solution adopted here is to compute a temperature adjustment, delta_artm, that results
in SMB ~ 0 over the initial glacier area. The adjustment is applied throughout the spin-up
to minimize glacier advance and retreat.
Recall the SMB formula:
SMB = P_s - mu_star * max(T - T_mlt, 0)
This can be modified to
SMB = P_s - mu_star * max(T + dT - T_mlt, 0),
where dT = delta_artm is the desired correction.
Setting SMB = 0, summing the remaining terms over an annual cycle, and ignoring the max operation,
we can solve for dT. Ignoring the 'max' yields an undershoot, but after a few annual iterations,
the SMB approaches zero as desired.
When inverting for powerlaw_c, delta_artm is automatically written to the restart file.
If switching from an inversion run (set_powerlaw_c = 1) to a forward run (set_powerlaw_c = 2)
and starting the forward run from the restart file, a nonzero delta_artm will be in the restart file.
In this case, delta_artm is automatically reset to zero as desired for the forward run.
I added a new 2D field, smb_obs, in the climate derived type. This is the field read from the input file.
In the glacier derived type, I removed the 2D smb_obs field and added two 1D fields: smb and smb_obs.
These are glacier-average fields that can be computed and output as diagnostics.
I also added delta_artm as a 1D glacier-average diagnostic field.
At the start of the run, the 2D smb_obs (aka climate%smb_obs) is read from the input file.
It is converted to the 1D glacier%smb_obs, which is passed repeatedly to the inversion subroutine
during spin-up. The 2D field is not used again.
On restart, the 1D glacier_smb_obs (aka glacier%smb_obs) is read from the restart file,
and the 2D field is not needed.
To compute the adjusted temperature, I added subroutine glacier_adjust_artm.
In several places, it is necessary to gather data over the 2D grid to compute 1D glacier average values,
or to scatter the 1D glacier-average values to the 2D grid. To consolidate these computations,
I added subroutines glacier_2d_to_1d and glacier_1d_to_2d.
In a recent commit, the inversion procedure for powerlaw_c (and coulomb_c) was modified to include a relaxation term. This term is a function of the ratio of powerlaw_c to a default value, and nudges powerlaw_c back toward that value. Thus, powerlaw_c is not continually pushed toward the extreme max or min value. Instead, the thickness error term is eventually balanced by the relaxation term. This commit implements similar logic for powerlaw_c inversion in glaciers. It requires a new parameter, glacier_powerlaw_c_relax_factor, which for now is declared at the top of glissade_glacier.F90 with a value of 0.05 (unitless). Later, we could add this and other parameters to the glacier derived type and make them config parameters. This commit is answer-changing for glacier spin-up with inversion for powerlaw_c.
In glacier runs to date, the snowfall rate has been read directly from the forcing file and applied without corrections, regardless of the downscaled surface temperature. This commit adds a new option in the glacier derived type: glacier_snow_calc. Option 0: Take the snowfall rate directly from the input field 'snow'. Option 1: Compute the snowfall rate from the precip and the downscaled artm. Option 1 is the default, anticipating that we will likely use this option going forward. The option can be set in the [glacier] section of the config file. Precip is assumed to fall entirely as snow at air temperatures below snow_threshold_min, and entirely as rain at temperatures above snow_threshold_max. At intermediate temperatures, the precip fraction that falls as snow follows a linear ramp. The two threshold values can be set in the [glacier] section of the config file. In the monthly climatological input file, I plotted the ratio snow/precip and compared it to contours of artm in different months. Based on this comparison, I chose default values of snow_threshold_min = -5 C and snow_threshold_max = 5 C. To be consistent, I changed the default of T_mlt to -5 C. I added a field snow_dartm_2d (analogous to Tpos_dartm_2d) to accumulate the monthly snowfall for the case that artm is adjusted. I also added a short subroutine, glacier_snow_calc, to compute the snowfall rate given the precip and surface air temperature. I also added glacier_delta_artm to glide_vars.def and added it to the list of restart fields when running glaciers with powerlaw_c inversion. This is needed for exact restart. In a few places, I replaced model%climate%artm with model%climate%artm_corrected. The latter is needed if running with a prescribed artm anomaly.
This commit adds a new way to read forcing files, in anticipation of GlacierMIP3. For GlacierMIP3, we will run to steady state with forcing from a given 20-year period, e.g., 2081-2100. Instead of cycling through these 20 years repeatedly, we will alternate years at random (actually, based on a list that was generated randomly). Thus, in a 2000-year run, we will use each year of forcing data about 100 times. Reading in a new forcing time slice every model month is expensive. An alternative is to read in all the forcing data just once, at initialization, and store it in a 3D array in which the third dimension is a time index. For GlacierMIP3, the time index runs from 1 to 240 (20 years * 12 months). To activate the new option, the user should set read_once = .true. in the [CF forcing] section of the config file. The default is read_once = .false, which gives the standard behavior. It is allowed to have two or more forcing files, with one or more read in the standard way, and one or more read in the new way. Any forcing file that can be read in the new way can also be read in the standard way, with results that are BFB. The user should check that the 2D fields to be read once are assigned 'read_once: 1' in glide_vars.def, and that the corresponding 3D fields are declared in glide_types. Currently, three fields have read_once = 1: precip, artm_ref, and snow. The associated 3D arrays are precip_read_once, artm_ref_read_once, and snow_read_once. These fields are used to compute glacier SMB. (Typically, either precip or snow is used, but not both.) The forcing file also contains the field usrf_ref, but this is the same for all time slices, so we don't need to save a 3D version. To enable this option, I added two subroutines to ncdf_template.F90.in: *_read_forcing_once, which reads in all time slices of the selected fields and stores them in 3D arrays *_retrieve_forcing, which copies data from the appropriate time slice to the standard 2D arrays Here, * = 'glide', 'glad', etc. So far, this option is used only for glide. The new subroutines are autogenerated in files glide_io.F90, glad_io.F90, etc. I modified generate_ncvars.py to insert the appropriate code for each variable with read_once = 1. The new subroutines are called from subroutines cism_init_dycore and cism_run_dycore, respectively, in cism_front_end.F90. The subroutine glide_read_forcing is called as before, to handle the forcing files with read_once = .false. A related change: In both glide_read_forcing and glide_read_forcing_once, I set the roundoff parameter eps = 1.d-3. This ensures that single-precision time values to the nearest month (e.g., 1979.0833) are interpreted correctly. The old value of 1.d-4 in glide_read_forcing allowed roundoff errors.
GlacierMIP3 specifies that future forcing should be read in based on a shuffled list of forcing years. For example, suppose the forcing is from years 2081-2100. Then the protocol states that the first few years of forcing data should be from years 2081, 2094, 2098, 2084, and 2090 (which were determined at random); and similarly until the end of the simulation. With this commit, a new attribute called shuffle_file can be specified in the [CF forcing] section of the config file. This attribute is a string containing the name of an ASCII file. The file, if present, is opened and read on each timestep, to determine whether it is time to read a new time slice and if so, the appropriate forcing year for the time slice. For the years above, we would first read 12 months of 2081 data, followed by 2094, and so on. If no shuffle file is provided, the code defaults to reading the forcing data corresponding to the model time. The ASCII file should consist of two columns of integers. The first is a consecutive list of years (0, 1, 2, 3, ...), and the second is the list of years. These are assumed to have Fortran format '(i6,i8)'. Note that the first column starts at year 0. I verified that the code works as expected for a sample GlacierMIP3 forcing file.
This commit introduces a new glacier config parameter called snow_reduction_factor. This parameter is a number between 0 and 1. It determines the fraction of incoming snowfall that is allowed to accumulate in grid cells outside the initial glacier mask. The initial mask is usually based on RGI observations. The motivation is as follows: * We want each glacier's steady-state area to be as close as possible to the initial area. Thus, we adjust either mu_star or delta_artm for each glacier so that the net SMB over the initial mask is close to zero. * Some glaciers will expand laterally past the initial boundary. To limit expansion, we apply the full computed ablation in these grid cells. * We still have the freedom to adjust the fraction of the computed snowfall applied in grid cells outside the initial mask. * If all the snowfall accumulates, we generally get too much expansion, and the steady-state ice volume is too high. But if no snowfall is allowed to accumulate, the steady-state ice volume is too low. * As a compromise, we allow a prescribed fraction of the input snowfall to accumulate. Some trial and error shows that snow_reduction_factor = 0.4 to 0.5 works well in many cases. The default value is 0.5. I also added a glacier parameter called diagnostic_minthck, which sets a thickness threshold for purposes of glacier area and volume diagnostics. For instance, a threshold of 10 m means that glacier ice thinner than 10 m does not contribute to the diagnosed glacier area. This makes it easier to match the nominal area and volume targets for each glacier. This parameter has no effect, however, on the dynamics. Typically, thklim and glacier%minthck are set to a smaller value of 1 m. For powerlaw_c inversion, the glacier code now uses several parameters that are part of the inversion derived type: babc_timescale, babc_thck_scale, and babc_relax_factor. Previously, these were hardwired parameters in the glacier module. Now we use the same parameters as are used for ice-sheet inversion. I changed the names of output fields glacier_area_target and glacier_volume_target to glacier_area_init and glacier_volume_init. I also added some useful diagnostic print statements for large glaciers.
The original way of computing mu_star was to adjust mu_star such that the modeled SMB = 0 when integrated over each glacier. A recent commit introduced the capability to adjust mu_star such that the modeled SMB matches the observed SMB, given the input temperature and snow forcing. At the same time, we compute a temperature correction delta_artm such that the modeled SMB over each glacier = 0 after the correction. The first method does not require smb_obs in the input data set, while the second method does. There were some logic issues when smb_obs was present but was not needed, or was needed but was not present. To more easily handle the logic, this commit introduces a new glacier config option, match_smb_obs, which is false for the first method and true for the second method. The default is false. With match_smb_obs = F, CISM will zero out smb_obs, if present in the input file. With match_smb_obs = T, CISM will throw a fatal error if smb_obs is missing in the input file. For match_smb_obs = F, the input temperature forcing should be appropriate for a period when the glacier was in balance with the climate. For match_smb_obs = T, the input temperature forcing should match the period of SMB observation. I also fixed a minor error in a diagnostic SMB calculation.
Until now, inversion of mu_star has been based on either of two criteria:
(1) Compute mu_star such that smb = 0, integrated over the initial glacier footprint
during a balanced climate (e.g., mid 20th century)
(2) Compute mu_star such that smb = smb_obs, integrated over the initial glacier footprint
during a period of SMB observations (e.g., Hugonnet).
In either case, we used the input snowfall without adjustment.
However, the snowfall for many glaciers is inaccurate, leading to inaccurate mu_star.
With this commit, we can invert for two glacier-specific parameters: mu_star and snow_factor,
where snow_factor is a scalar that multiplies the observation-based snowfall.
That is, the SMB is given by
SMB = snow_factor * snow - mu_star * max(T - Tmlt, 0).
We enforce both (1) and (2), resulting in a system of two equations and two unknowns.
This system is solved for mu_star and snow_factor for each glacier.
With the extra degree of freedom, spun-up glacier areas and volumes are generally
in better agreement with the initial values.
For some glaciers (usually small ones), mu_star must be adjusted to fall within an allowed range,
in which case (1) is satisfied but not (2).
To run the 2-parameter inversion scheme, the user should specify set_mu_star = 1 and set_snow_factor = 1
in the config file.
To run a 1-parameter scheme that enforces criterion (1) only, the user should specify set_mu_star = 1
and (optionally) set_snow_factor = 0. If set_snow_factor is not specified, it defaults to 0.
The 2-parameter scheme requires reading in two sets of forcing data, typically usrf_ref, artm_ref, snow,
and/or precip for both the balanced climate and unbalanced climate.
The input fields for the balanced climate are read as before.
The input fields for the unbalanced climate are read from a separate forcing file containing auxiliary fields
called usrf_ref_aux, artm_ref_aux, snow_aux, and precip_aux.
If both are specified in the config file under [CF forcing], CISM will read and handle them correctly.
The field smb_obs remains (at least for now) in the input file, not a forcing file.
The method of inverting for the single parameter mu_star with smb = smb_obs is no longer supported.
This method required a temperature correction delta_artm, often with unrealistically large corrections.
I removed delta_artm from the code.
A recent commit added the option to read forcing files just once at the start of a run, copy the date to arrays, and read data from these arrays at runtime, instead of repeatedly reading data from a netCDF file. This commit adds the ability to do this with multiple files. For example, during the spin-up we can have one forcing file with artm_ref, precip, and/or snow from a balanced climate, and another with auxiliary fields artm_ref_aux, precip_aux, and/or snow_aux from a recent unbalanced climate. Until now, the standard spin-up procedure has been to read both files at each monthly time step. But since there are only 12 monthly values in each file, it is more efficient to read and save all 12 values at start-up. With this commit, CISM can read multiple forcing files just once at initialization. All the fields to be read once should be listed in glide_vars.def with read_once = 1. Each file should have read_once = .true. in the [CF forcing] section of the config file. New logic in subroutines glide_read_forcing_once and glide_retrieve_forcing (both in the autogenerated glide_io.F90) will take care of the rest.
This commit introduces a more systematic way to compute SMB for
(1) advanced glacier cells (cism_glacier_id_init = 0 but cism_glacier_id = 0) and
(2) cells adjacent to glacier-covered cells (cism_glacier_id = 0 but cism_glacier_id > 0
for a neighbor).
We use masks to determine where to apply the computed SMB. There are two versions of the mask:
- smb_glacier_id_init, based on the initial glacier footprints (from cism_glacier_id_init).
This mask is used for inversion of mu_star and snow_factor.
- smb_glacier_id, based on the current glacier footprints (from cism_glacier_id).
This mask determines where the SMB is applied at runtime.
The rules for smb_glacier_id are as follows:
- Where cism_glacier_id_init > 0 and cism_glacier_id > 0, set smb_glacier_id(i,j) = cism_glacier_id(i,j)
and apply the SMB.
- In advanced grid cells (cism_glacier_id_init = 0 but cism_glacier_id > 0),
compute a potential SMB assuming smb_glacier_id(i,j) = cism_glacier_id(i,j).
Apply this SMB if negative; else set smb_glacier_id(i,j) = 0.
- In retreated grid cells (cism_glacier_id_init > 0 but cism_glacier_id = 0),
compute a potential SMB assuming smb_glacier_id(i,j) = cism_glacier_id_init(i,j).
Apply this SMB if positive; else set smb_glacier_id(i,j) = 0.
- In other glacier-free cells (cism_glacier_id_init = cism_glacier_id = 0), check
for glacier-covered edge neighbors (cism_glacier_id > 0). For each neighbor (ii,jj),
compute a potential SMB assuming smb_glacier_id(i,j) = cism_glacier_id(ii,jj).
Apply this SMB if negative; else set smb_glacier_id(i,j) = 0.
If there are neighbors with SMB < 0 from two or more glaciers, choose the glacier ID
that results in the lowest SMB.
The rules for smb_glacier_id_init are the same as for smb_glacier_id, except that
we assume cism_glacier_id = cism_glacier_id_init, so there are no advanced cells.
The biggest change, compared to previous code, is that we no longer allow accumulation
outside the initial glacier extent.
Other changes:
I added some code to prevent 'pirating' of one glacier by another.
This can happen when two glaciers are adjacent, but one has larger mu_star and hence
more ablation for a given climate. If a particular advanced cell winds up in the glacier
with less ablation, it can steal ice from the adjacent glacier, allowing the
slow-melting glacier to advance unrealistically.
The fix is to transfer advanced cells as needed to the glacier with more ablation.
I added a subroutine to remove "snowfields", defined as patches of isolated ice,
isolated from the initial glacier extent. The code to remove snowfields is similar
to the existing code for removing icebergs.
I made Tmlt a 2D field, anticipating that we may allow it to vary spatially in future commits.
For now, it is set everywhere to tmlt_const, a config parameter that is set to -4 C by default.
I renamed the main runtime glacier subroutine from glissade_glacier_inversion to glissade_glacier_update.
The advance_retreat subroutine is now called from inside this subroutine instead of the glissade module.
I added diagnostic subroutines to compute accumulation area ratios and areas of advance
and retreat for each glacier.
I removed the deprecated snow_reduction_factor.
The 2-parameter inversion scheme solves a pair of coupled equations:
p * S - mu * Tp = 0
p' * S' - mu * Tp' = B
where p = snow_factor, mu = mu_star, S = snow, Tp = max(T - Tmlt, 0),
B = observed SMB, and a prime denotes the auxiliary climate associated
with a (mostly) negative mass balance.
The solution can be written as
mu = -B * S / D
p = -B * Tp / D
D = S * Tp' - S' * Tp
Thus, D is a snowfall-weighted temperature difference; D > 0 for a warming climate,
provided S and S' are not too different.
For the majority of glaciers, mu_star and snow_factor fall within realistic ranges.
The ranges (for now) are specified as (200, 5000) for mu_star and (0.5, 3.0) for snow_factor.
However, many glaciers have (1) mu_star > mu_star_max or (2) mu_star < mu_star_min.
Some glaciers have mu_star < 0, meaning that B > 0 is associated with warming
or B < 0 is associated with cooling, which is unrealistic.
Until now, we've simply set mu_star = mu_star_max in case (1) or mu_star = mu_star_min in case (2).
With this commit, the code corrects the auxiliary temperature (artm_aux) to bring mu into the desired range.
With mu_star in the desired range, snow_factor is usually in the desired range also.
Strictly speaking, this is a 3rd parameter in the inversion, but this 3rd parameter
is only used to assist in finding sensible values of mu_star and snow_factor;
it is not used in subsequent forward runs.
The correction is a new glacier-specific variable called artm_aux_corr.
It is limited to be no greater than 3 C in either direction.
The adjustment stops when mu is in the prescribed range or artm_aux_corr reaches
its limit, whichever comes first.
Initially, I tried adjusting Tmlt on a glacier-by-glacier basis, but this turned out
to be numerically problematic. Tmlt is now the same for all glaciers, as it was earlier.
This commit is answer-changing for many glaciers.
When writing equations, we have been calling the multiplicative snow factor 'alpha', and the temperature correction factor 'beta'. This commit changes snow_factor to alpha_snow and changes artm_aux_corr to beta_artm_aux, consistent with this notation. This commit is BFB.
The following parameters can now be set by the user in the config file, instead of being hardwired in the glacier module: - mu_star_const, mu_star_min, mu_star_max - alpha_snow_const, alpha_snow_min, alpha_snow_max - beta_artm_aux_max, beta_artm_aux_increment
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.
|
@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. |
TimvandenAkker
pushed a commit
that referenced
this pull request
Jan 22, 2025
Katetc/hma glaciers4, Replace Lipscomb/hma_glaciers4 PR
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.
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:
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.