New physics options to support ISMIP7 simulations - #91
Open
whlipscomb wants to merge 43 commits into
Open
whlipscomb wants to merge 43 commits into
whlipscomb wants to merge 43 commits into
Conversation
This commit supports the final CISM submission for CalvingMIP. Diagnostic additions: * CISM now computes the radial distance from the origin, the effective thickness, and the ice speed at the calving front for each of 2 axes for Experiments 2 and 4. For the circular domain, the axes are the y-axis and the line y = x (in the NE quadrant). For the Thule domain, the axes are Caprona A and Halbrane A. These 4 axes are the ones plotted in the CalvingMIP paper by Jordan et al. (2026). I added the scalars cf_radius1, cf_radius2, cf_thck1, cf_thck2, cf_speed1 and cf_speed2 to glide_types and glide_vars.def. This means that users will be able to compare results directly to the published values from CISM and other models without doing any post-processing. Still to do is to extend these diagnostics to all 8 axes per experiment. Algorithmic changes: * In full CF cells, H_eff(i,j) is always computed based on the interior neighbor thickness instead of capping H_eff(i,j) at H(i,j), so we can have H > H_eff temporarily. * If a CF cell has no interior edge neighbor, it now computes H_eff from an interior corner neighbor. * For calvingMIP experiments, CISM now calls subroutine advance_calving_front after applying the calving_based thickness change. This ensures H <= H_eff at the CF when the algoritm is finished. I modified this subroutine so that ice advances only into edge neighbors and not corner neighbors. This prevents the creation of thin interior cells when the CF advances. * I fixed a bug in the computation of fluxes from corner neighbors in glissade_input_fluxes. * The protected CF cells (i.e., cells that are allowed to accumulate ice rather than have that ice shifted back upstream) now include cells with two partial CF neighbors (which typically have a full interior diagonal neighbor). This allows ice to accumulate in cells at CF corners. * I modified remove_icebergs so that any ice-covered cell, including partial CF cells, can spread the fill to edge neighbors. With this change, we don't need to call the glissade_fill_with_buffer subroutine (which was written to account for partial CF cells that were dynamically inactive), but just the standard glissade_fill. * In glissade_utils, I added subroutines that (1) compute global sums over the four calvingMIP quadrants and (2) compute the value of a field at a point in a bounding box, given the values at the corners. With these changes, the CF advance/retreat rates are similar to before (as desired), and oscillations in H_eff and speed at the CF are significantly reduced, especially when the CF is advancing.
The 25-point Laplacian stencil contained a '26' instead of '36' in the mask computation. Fixed the typo.
This change makes CISM's no-advance calving mask more flexible. Previously, this mask was an integer, either 0 or 1. Now the mask can be either (1) an integer array (still called calving_mask) if which_ho_calving_front = 0, or (2) a real array (called subgrid_calving_mask) if which_ho_calving_front = 1. When using the real array, CISM will thin ice in partial CF cells so that areafrac = thck/thck_effective = (1 - subgrid_calving_mask). The subgrid calving mask can be used when which_ho_calving_front = 1 (i.e., the subgrid CF scheme). The mask is computed at initialization. It can be based on an optional config parameter called calving_front_radius. If calving_front_radius /= 0, subgrid_calving_mask is computed to enforce this radius. This option is now used for CalvingMIP spinups. By default, the mask is computed based on the initial ice thickness, with a value of 0.0 or 1.0 everywhere. I applied the new subgrid mask to the two calvingMIP spinups (circular and Thule). Previously, the final CF radius was very close to 750 km along the x- and y-axes, but around 752 km along the diagonals because of the binary mask. Now the CF radius along diagonals is very close to 750 km (within < 0.1 km). I combined three mask operations (apply binary mask, apply subgrid mask, and apply forced retreat) in one subroutine, glissade_apply_calving_mask, which is now called from glissade_calving_solve after doing the main calving. Also, I improved the CalvingMIP diagnostics. Previously, the CF radius, speed, and thickness were computed only for axes 1 and 2 of each experiment. Now they are computed for all 8 axes. There are output fields called cf_radius, cf_speed, and cf_thck, each with dimension 8. I created a new dimension, naxis = 8, to support this output. I made the logic more consistent so that in general, the CF location is found by interpolating between two points, one with a_eff > 0.5 and the other with a_eff < 0.5. I put some of the complicated Caprona logic in subroutines to allow code reuse. Other changes: * The value of thck_effective now cannot exceed the flotation thickness. * I reverted an earlier change: the iceberg removal subroutine again calls glissade_fill_with_buffer instead of glissade_fill. Using the buffer allows ice to accumulate, as desired, in ice-free ocean cells at the CF. * For calving masks, thck_effective is not set to 0.0 in all cells with ice_mask = 0, and thck_effective is not allowed to exceed the flotation thickness. * I wrote a subroutine, glissade_cleanup_tiny_thickness, which is now called after the calving, at the end of the prognostic part of the timestep. This subroutine removes ice from cells with H < 1.e-11 m, adding it to the calving flux. This avoids starting the next timestep with tiny but nonzero amounts of ice in cells near the calving front. This subroutine, along with glissade_cleanup_icefree_cells, is in the glissade_utils module. * The 'thklim' argument in some calls to glissade_get_masks now depends on the CF numerics. With a subgrid CF, the threshold is 1.e-11 (so partial cells with H < 1 m can have ice_mask = 1). Without a subgrid CF, the threshold is model%numerics%thklim (typically 1 m). * I reduced the small_dthck parameter in subroutine advance_calving_front to 0.1 m, to avoid a too-steep gradient at the CF. * In subroutine add_surface_and_basal_mass_balance, most computations are now skipped in cells with effective_areafrac = 0. * Quadrant sums are now reproducible for any processor count. With these changes, all the CalvingMIP spinups and experiments are performing well.
This commit contains changes in the calvingMIP test directory, specifcally the files calvingMIP.config.template, calvingMIP.Setup.py, and README.calvingMIP. These are the files used to set up the final CISM runs for the CalvingMIP paper. Among the changes: * Both spin-ups now use marine_margin = 6, calving ice at the margin based on subgrid_calving_mask. The calving amplitude is set to a large positive value to allow free advance until reaching the margin. * The spin-ups and Experiment 4 have apply_calving_mask = True. These runs need a mask to limit CF advance when no other calving is taking place. * All runs now have which_ho_flotation_function = 2 (which is the same as the old 3). * Experiments 1 to 4 write cf_radius, cf_thck and cf_speed to the scalar output for each of 8 axes for each experiment. * In the input files, topg is computed in a way that enforces symmetry. * The input files no longer contain calving_mask. Rather, subgrid_calving_mask is computed at startup. * The prescribed value of dthck_dx_cf is now 1.e-4 instead of 5.e-4. The smaller value gives slightly more accurate CF retreat.
This commit changes the Caprona CalvingMIP diagnostics to be more symmetric. I replaced bounding_box calculations with more accurate interpolations along the Caprona axes, and removed subroutine glissade_bounding_box. I added four new diagnostic fields: * cf_locx and cf_locy, the x and y coordinates of the calving front along each axis * cf_uvel and vf_vvel, the ice speed at the CF along each axis I removed cf_speed, which is redundant given cf_uvel and cf_vvel. Also, I fixed a bug in the computation of thck_effective just before the velocity solve. I was computing thck_effective = 50 m (= thck_effective_min) at the calving front, instead of the correct value. (Elsewhere in the code, thck_effective was computed correctly.) This fix causes modest changes in ice velocity at the CF.
This commit enables a physically-based calving law that combines stress-based and thickness-based calving. To use this calving law, the user should set marine_margin = CALVING_STRESS = 8 and set calving_minthck to a positive value. Then during the calving, there will be a call to subroutine stress_based_calving following by a call to subroutine thickness_based calving. Each subroutine computes a thinning term, calving_dthck, in each CF cell. The model then applies whichever term has the greater magnitude. In this way, the calving rate increases as ice thins, even if the principal stresses do not exceed the threshold values. Other changes: * I added basin-scale diagnostics: iarea_basin, iareaf_basin, iareag_basin, ivol_basin, ivol_above_flotation_basin, imass_basin, imass_above_flotation_basin. These are computed in the same way as the global diagnostics iarea, iareaf, etc., but using basin masks so that each sum is over a single basin. I also added a global scalar called ivol_above_flotation. * I added a basin dimension in the glide_nc_custom module to support basin-scale output. * I added some logic so that when inverting for dT_ocn, there is no inversion beyond the initial calving front. Something like this was in place before but had been removed. The reason for this is that we want ice that advances beyond the initial CF to calve based on the calving law, without any help from the inversion. These changes are BFB except for stress-based calving with calving_minthck > 0. I did a number of AIS runs with stress-plus-thickness-based calving, with tau_eigenconstant2 = 1.0, tau_eigenconstant1 = 0.0, and calving_minthck = 250 m. The working hypothesis is that the first principal stress is much less important than the second, and that where the second principal stress is small, the CF will advance until the ice is too thin to remain intact. The spin-up looks good. There are some regions of CF retreat (e.g., Brunt, West and Ross) and advance (e.g., Ronne, PIG, and west of Thwaites), but these are small compared to the total extent of floating ice. In forward runs with warming, this calving law leads to faster Thwaites collapse than runs with a no-advance calving mask.
In glide_types, nbasin is now initialized to 1 instead of 0. Also, the basin_number array is now always allocated. This avoids some possible I/O errors and segmentation faults in runs without multiple basins.
This commit generalizes the Schoof basal sliding law, which can be written as:
tau_p * tau_c
tau_b = -------------------------
[tau_p^n + tau_c^n]^(1/n)
where tau_p = Cp * u_b^(1/m) is the basal shear stress from a Weertman-type power law,
tau_c = Cc * N is the basal shear stress from a Coulomb law, and tau_b is the resultant stress.
Typically, n = m, but the Schoof law has the same asymptotic behavior when n /= m.
The larger (smaller) the value of n, the sharper (more gradual) the transition
between powerlaw and Coulomb behavior.
With this commit, m and n can be different. The defaults are powerlaw_m = schoof_n = 3.0,
but other values can be set in the config file as powerlaw_m and schoof_n, respectively.
Other changes:
* The user can now invert for Cc and Cp simultaneously (but not independently)
for the Schoof and Tsai laws, assuming that Cp and Cc are related by Cp = gamma * Cc^p,
where gamma (= schoof_gamma in the code) and p (= schoof_p) can be set in the config file.
To invert for both, the user sets which_ho_coulomb_c = 1 and which_ho_powerlaw_c = 4
(a new option that implements the above relationship).
Alternatively, this can be done at basin scale with which_ho_coulomb_c = 3 and
which_ho_powerlaw_c = 4.
The default values are gamma = 1.e5 and p = 0.70. With these values, Cc_max = 1.0
maps to Cp_max = 1.0e5, and Cc_const = 0.1 maps to Cp_const = 20,000 (assuming m = 3).
The underlying assumption is that a bed with strong (weak) resistance to Coulomb sliding
will also have strong (weak) resistance to powerlaw sliding.
I'm not sure to what extent this assumption is supported by evidence.
* I renumbered some babc options. The recent modified Schoof option 12 is superseded
by the generalized Schoof law (option 11) with schoof_n = 1. The Tsai law is option 12 again.
* I removed the previous basin-scale inversion for Cc_hi and Cc_lo, since the
logic is complicated and the results weren't very satisfactory.
* I modified the effective pressure computation based on the parameter p_ocean_penetration.
Previously, this scheme computed N = overburden on land.
Now, N on land is capped at rhoi*grav*basal_physics%haf_threshold.
The threshold is 500 m by default.
* I removed ocean_p_timescale, which phased in the p_ocean-based value of N with a time lag.
This parameter was added a few years ago to prevent Thwaites collapse early in the run.
Thwaites is now more stable, so I'm no longer running with a nonzero timescale.
* The initial calving mask is now set to 0 in cells where each of the four edge neighbors
has calving_mask = 0. This prevents the model from creating one-cell holes in the center
of ice shelves. This avoids slow convergence in several grid cells.
(In the future, we might want to pre-process the thickness targets to avoid holes.)
I ran several AIS spin-ups with these changes:
* I compared schoof_n = 1, 2 and 3 for local CC/Cp inversion.
Convergence is fastest for schoof_n = 1: 21 nonlinear iterations per timestep
at the end of a 5 kyr AIS spin-up, compared to 24 for n = 2 and 26 for n = 3.
(These numbers could vary with different model settings.)
* The fit to observed thickness and velocity is good with the combined
Cc/Cp inversion, with rmse of ~26 m for thickness and 137 m/yr for velocity.
Revised some text written to the log file from glide_setup, based on the current basal friction options.
This commit modifies the calving mask computed from the initial ice extent. Until now, all cells that are initially ice-free ocean (ocean_mask = 1) were assigned calving_mask = 1 by default and kept ice-free during the run. This includes interior ice-free cells disconnected from the deep ocean. However, keeping these cells ice-free seems physically unrealistic and can lead to fluctuations that slow convergence. With this commit, the default behavior is to first set calving_mask = 1 in ice-free ocean cells and then do a flood fill to identify disconnected cells (ocean_connection_mask = 0). In these cells, calving_mask is set to 0. On the current AIS 8km input file, there are about 40 such cells. This commit is potentially answer-changing for any run with marine_margin = 5.
In a previous commit, I turned off the Laplacian term when inverting for coulomb_c or powerlaw_c, except for vertices that are currently inactive. At that time, turning off the Laplacian term seemed to improved convergence. With this commit, the Laplacian is computed at all vertices where inversion is occuring, active or inactive (provided inversion_babc_length_scale > 0). This gives a smoother Cc or Cp field that seems more realistic, and it doesn't slow convergence in current AIS spin-ups.
The various inversion schemes include a term proportional to dH/dt. This is a damping term that helps prevent the value of an inversion variable from oscillating around the solution. Until now, the damping term has included a hardwired factor of 2. For example, the Cp and Cc inversions include a term 2 * dH/dt / H0. With this commit, the user can replace the factor of 2 with a value set in the config file (inversion_damping_factor in the parameters section). The default value is 2.0, for backward compatibility. In AIS spin-ups, I experimented with different values, including schemes that start with a larger value and gradually transition to a smaller value. It turns out that a single time-independent value works well. I got good results with damping_factor = 1 or 2. Smaller values appear to be underdamped, giving a larger rmse for ice thickness. In current AIS spin-ups, the time to run 1000 yr on 256 cores is less than 50 minutes using the Schoof sliding law with schoof_n = 1. Before the latest round of commits, Schoof-law runs were taking over an hour.
The logic of the basin-scale inversion for deltaT_ocn assumes that deltaT_ocn is uniform across each basin. However, a recent commit added some logic to set deltaT_ocn = 0 during inversion in cells where calving_mask = 1. This logic makes sense for local (2D) inversion, but it breaks basin-scale inversion. This commit fixes the logic by applying it only for local inversion, not basin-scale. Also, the parameter 'damping_factor' is now passed into the subroutines for basin-scale inversion, not just local inversion. The goal is to use basin-scale rather than local inversion for ISMIP7, to reduce the number of tunable parameters.
This commit adds a new option, which_lateral_melt. It is based on a recent implementation by Heiko Goelzer in an older version of CISM. The choices are (default = 0): [0] No lateral melt [1] Constant lateral melt rate [2] ISMIP6 lateral melt rate [3] ISMIP6 lateral melt rate for coupled setup Main additions: * I defined a new derived type, glide_lateral_melt, with several arrays and scalars. * There are several new arrays in glide_vars.def. * There is a new optional section 'lateral_melt' in the config file. * There is a new module, glissade_lateral_melt, with subroutines for computing melt_thck (the thickness change due to lateral melt at the ice front) for options 1 and 2. * There is a new subroutine in glissade_masks which is like glissade_calving_front_mask except that it allows grounded ice (not just floating) at the marine margin. The new code compiles but isn't yet functional. I would like to support lateral melt with a subgrid front parameterization, similar to what we do for calving. The easiest way to do this is probably to add lateral melt as a forcing term (in addition to calving) in the calving module. This will be done in an upcoming commit.
This commit creates a new module called glissade_diagnostics. I moved three CalvingMIP diagnostic subroutines from glissade_calving to the new module. The goal is to reduce the length of the calving module while creating a suitable place for this and other diagnostic subroutines. More subroutines can be added later. For example, the subroutines now in glide_diagnostics could go here, as could some of the inline diagnostics in glissade.F90. This commit is BFB.
This commit moves subroutine glissade_calving_solve from subroutine glissade_tstep (in glissade.F90) to glissade_calving.F90. Several other subroutines in glissade_calving, which are called from glissade_calving_solve, are now private to the calving module. This is part of a general effort to shorten glissade.F90 by moving high-level drivers to the appropriate physics modules. I also moved subroutines glissade_quadrant_sum (now called sum_over_quadrants), glissade_stress_tensor_eigenvalues, and glissade_strain_rate_tensor_eigenvalues from glissade_calving to glissade_diagnostics. This commit is BFB.
This commit divides subroutine glissade_calve_ice into two shorter subroutines: calve_ice_basic and calve_ice_subgrid. The former handles the simple position-based calving options, most of which date back to Glide. The latter handles the more recent options using a subgrid CF. The commit is BFB.
This commit adds a new option, which_ho_calving_front = 2. The new option is a generalization of which_ho_calving_front = 1. The main difference is that option 1 allows only floating cells to be CF cells, whereas option 2 allows both floating and marine-grounded cells to be CF cells. The lateral melt scheme will use option 2. CF cells are still identified by subroutine glissade_calving_front_mask, but that subroutine now contains logic to support the new option. I removed some variables from the lateral_melt derived type, since lateral melt is going to be coordinated with subgrid calving. I also removed subroutine glissade_melt_front_mask, which was originally going to handle lateral melt at grounded fronts. Fronts with potential lateral melt will now be identified by the new subgrid CF option. I verified that answers are BFB when which_ho_calving_front = 1, but still need to test which_ho_calving_front = 2.
This commit adds subroutine glissade_calc_lsrf_usrf in glissade_utils. The new subroutine computes lsrf and usrf everywhere, given thck, topg and eus. It has the same functionality as calling glide_calclsrf followed by a usrf update, but I'm working toward removing glide dependencies in glissade. I also removed some unnecessary arguments in the scaling subroutines in glissade_velo_higher. The only necessary scalings are those converting between units of s and yr. This commit is BFB.
Ice caps are no longer removed in glissade_diagnostic_variable_solve, since state variables aren't supposed to be modified in this subroutine. I created a new subroutine, glissade_remove_ice_caps, in glissade_utils. This subroutine is now called during initialization (if not restarting) and near the end of each timestep (after calving). I also changed the block_inception logic so that the ice_sheet_mask array is updated just before applying the logic. This commit is BFB in the cases tested.
At the end of each timestep, CISM computes a mass budget, checking that the total change in ice mass is equal to the sum of the mass fluxes. Until now, the mass fluxes have been divided into three categories: SMB, BMB and calving. This commit adds two more: lateral melt and ice removal. If nonzero, the two new fluxes are now written to the diagnostic log file in addition to SMB, BMB and calving. I checked that mass is still conserved within roundoff. Ice removal includes the mass loss associated with (1) ice cap removal and (2) removal of tiny amounts of ice from nearly ice-free cells. Iceberg and isthmus removal are still counted as part of the calving flux. I put all the flux terms in a new derived type called glide_mass_flux. (Several were previously in the geometry derived type.) I added the necessary flux diagnostic computations in the glissade and glide_diagnostic modules, and added several new flux fields in glide_vars.def. I added two subroutines in glissade_diagnostics.F90: (1) glissade_mass_balance_diagnostics and (2) glissade_grounding_line_flux. Both are called from the last part of the diagnostic solve. The first consists of code previously in the glissade module, and the second was moved from the glissade_grounding_line module. I also added calving_front_mask to the calving derived type. This commit is BFB apart from the new diagnostics.
This commit adds diagnostic for computing the total area and volume of ice caps. These diagnostics are now written to the log file (unless ice cap removal is turned on, in which case both values are 0.) I renamed subroutine glissade_remove_ice_caps to glissade_handle_ice_caps. If ice cap removal is turned off, the subroutine computes ice-sheet and ice-cap masks for diagnostic output.
This commit adds some structure to the lateral melt calculations.
In module glissade_lateral_melt, there are now three subroutines:
(1) glissade_lateral_melt_constant, which computes a thinning term
given which_lateral_melt = 1.
(2) glissade_lateral_melt_ismip6, which computes a thinning term
given which_lateral_melt = 2.
(3) glissade_lateral_thermal_forcing_avg, which computes tforcing_2d,
the depth-averaged thermal forcing needed for (2).
Subroutine (3) takes ztop_tfavg and zbot_tfavg as inputs.
These variables specify the ocean depth over which TF is averaged.
Both are now part of the lateral_melt type and can be set in the config file.
The defaults are -200 m and -500 m, the standard ISMIP6 values.
These subroutines are called (depending on which_lateral_melt)
from subroutine calve_ice_subgrid in module glissade_calving.
They output a variable called latmelt_dthck, which is applied
immediately after applying calving_dthck.
The thinning logic is the same as for calving.
I verified that the commit is BFB for cases without lateral melt.
The new lateral melt code still needs to be tested.
I had put this subroutine in the glissade_utils module. Moved it to the glimmer_utils module to avoid a circular dependency.
Forced ice retreat is now part of the removal flux (which also includes ice cap removal) instead of the calving flux.
In the new calc_lsrf_usrf subroutine, the formula for lsrf for floating ice is now:
lsrf = eus - (rhoi/rhoo)*thck
The 'eus' term on the rhs was missing before.
This won't change answers (except maybe at roundoff level) with eus = 0,
but is an important correction when we start running with nonzero eus.
This commit fixes a bug in subroutine parallel_halo_integer_3d. The mpi_send and mpi_recv calls in this subroutine now pass the correct argument, mpi_integer, instead of passing mpi_real8. Thanks to Mariana and Heiko for spotting the bug.
This commit modifies the new subgrid calving front option, which_ho_calving_front = 2,
to improved stability.
The new option is similar to option 1, except that it allows CF cells to be either floating or
marine-grounded, with H_eff derived from interior neighbors that are floating or marine-grounded.
Option 1 requires CF cells to be floating.
This change can lead to instability if either (1) the new H_eff results in a large surface
elevation gradient at the CF, or (2) the new H_eff leads to a high cliff (e.g., > 100 m
above sea level) at the CF. The fix is to compute usrf_effective and apply two limiters:
(1) Limit the gradient (usrf_effective - usrf_neighbor)/dx to a prescribed max.
Note: usrf_effective is limited only if the CF cell has a surface higher than its interior neighbor.
(2) Limit the difference usrf_effective - usrf, where usrf is the state-variable elevation.
Given the new usrf_effective, it is straightforward to recompute thck_effective.
These limits introduce two new parameters: max_dusrf_dx and max_dusrf.
After some testing, I chose max_dusrf_dx = 0.001 and max_dusrf = 25 m.
With these values, typical AIS runs are stable, with no major changes in the maximum ice speed
or the number of nonlinear iterations needed to converge.
The parameters are hardwired for now, but they could be made user-configurable.
The logic for subgrid CF options 0 and 1 has not changed.
I also added subroutine glissade_subglacial_discharge in glissade_lateral_melt.F90.
The subroutine is called from the calving module when lateral melt is enabled.
The subglacial discharge, which is needed for the lateral melt parameterization,
is computed by assuming that all surface ablation (i.e., acab_applied, if negative)
reaches the bed and flows to the marine terminus. The discharge is then summed over basins.
This commit introduces two logical options for lateral melt: (1) subglacial_discharge_from_ablation (default = F) - if false, then read subroutine_discharge directly from file - if true, then compute submarine discharge from surface ablation (2) thermal_forcing_avg_3d_to_2d (default = F) - if false, then read 2d thermal forcing directly from file; - if true, then read 3d thermal forcing from file and average to 2d This gives the user the flexibility to either provide subglacial_discharge (in an input file or passed from the coupler) or compute it from ablation. Similarly, the user can provide either 3D or 2D thermal forcing. The code aborts if (1) subglacial_discharge_from_ablation = F and a subglacial_discharge field is not read in, (2) thermal_forcing_avg_3d_to_2d = F and 2d thermal_forcing is not read in, or (3) thermal_forcing_avg_3d_to_2d = T and 3d thermal_forcing is not read in. The 2D thermal forcing field is loadable and is now called thermal_forcing_2d. The 3D thermal forcing field (also loadable) is still called thermal_forcing. The lateral melt options are now: - LATERAL_MELT_NONE = 0 - LATERAL_MELT_CONSTANT = 1 - LATERAL_MELT_ISMIP = 2 I removed option 3 (COUPLED), since it is redundant given the new logical options above. I also fixed a minor initialization bug in subroutine glissade_calc_effecpress.
This commit adds supporting code for ISMIP7 Antarctic experiments with prescribed ocean thermal forcing and target sub-shelf melt rates. The protocols are described here: https://www.ismip.org/participants/focus-groups/ais-basal-melt https://drive.google.com/file/d/1SygMQte-7XgKj4e-Hpyj1tHi6XguShCP/view The protocol includes instructions for calibrating melt schemes. For the ISMIP6-based quadratic schemes supported by CISM (local, nonlocal and nonlocal-slope), the calibration works as follows: (1) Read in thermal forcing as usual; also read in target basal melt rates from observations. (2) For each value of gamma0 in an ensemble of possible values, compute deltaT_ocn_basin for each basin to match basin-average target melt rates. (3) Use an ISMIP7 toolbox to choose the best gamma0/deltaT_ocn_basin combination. The toolbox included additional TF and target melt rates to determine sensitivities. The main code changes are as follows: * Added a loadable 2D I/O field call bmb_obs. This is the observed basal mass balance (kg/m^2/yr), defined as positive for melting. Note: These units are equivalent to mm/yr w.e. (the units for smb, smb_obs, etc.) * Added a 2D field, bmlt_float_target, in the basal_melt derived type. This field is derived from bmb_obs, if present. For I/O, bmlt_float_target has units of m/yr ice, the same as bmlt_float. * Added a new option, which_ho_bmlt_float = HO_DELTAT_OCN_CALIBRATE_BASIN = 4. With this option enabled, the model calls a new subroutine, calibrate_deltaT_ocn_basin, at initialization. This subroutine does an iterative solve to find deltaT_ocn_basin in each basin such that the basin average of bmlt_float is equal to the average of bmlt_float_target in the domain where both values are defined. Note: bmlt_float_target derived from the ISMIP7 Paolo dataset has missing values in some regions where the CISM spin-up has floating ice present. These regions are not included in the average. * Changed which_ho_bmlt_float = HO_DELTAT_OCN_DTHCK_DT from option 4 to option 5. The logic to apply this option is now in glissade_bmlt_float_init instead of glissade_inversion. * Added thermal_forcing_mask to the basal_melt derived type so that it is available for diagnostics and I/O. * Added a basin-scale field, deltaT_ocn_basin, to the ocean_data derived type. * Moved several global scalars from the geometry derived type to a new 'scalars' derived type. Added total_bmlt_float, total_bmlt_float_target, bmlt_float_basin and bmlt_float_target_basin (all with units of kg/s) to this derived type. * Renamed glissade_bmlt_float_thermal_forcing_init to glissade_thermal_forcing_init. Moved some code there from glissade.F90 and glissade_inversion.F90. * Moved subroutine glissade_bmlt_float_solve from glissade.F90 to glissade_bmlt_float.F90. This makes it possible to call glissade_bmlt_float_solve from glissade_bmlt_float_init as needed for calibration. * Made sure that the units of bmlt_float and bmlt_float_target are consistently m/s within the code (converted to m/yr for I/O and diagnostics only). * Removed some custom parameters such as rhoi_ismip6; subroutines that used these parameters now use the standard parameters (rhoi, etc.) * Changed subroutine interpolate_thermal_forcing_to_lsrf to a generic subroutine called interpolate_3d_forcing_to_lsrf . * Removed the deprecated option which_ho_flotation_function = 3 and the associated field, model%geometry%topg_raised.
When running with ISMIP thermal forcing, bmlt_float is computed during each timestep, but until now has not been computed at initialization. With this commit, the user can compute bmlt_float at initialization (t = 0) by setting bmlt_float_init = .true. in the [options] section of the config file. (The default is .false.) This option is useful for ISMIP melt calibration, where we want to compute bmlt_float and write it to output files for many different input datasets and values of gamma0. It is convenient to set tend = 0 in these runs to minimize the number of output time slices. I confirmed that CISM can read in deltaT_ocn and other state variables from the restart file of an earlier calibration run (with restart option 2 = hybrid restart), while reading thermal_forcing from an ISMIP7 model output dataset (e.g., that of Mathiot). CISM then writes bmlt_float and other variables of interest to the initial output files. I also renamed 'bmb_obs' to 'bmb_float'. This is the basal melt rate in ISMIP7 input files with units of kg/m^2/yr. The 'bmb' suggests that it has the same units as smb, and the 'float' indicates that it applies to floating ice only. Often it will be from models rather than obs.
I found some significant energy conservation errors (> roundoff) in runs with subgrid calving
due to setting negative thicknesses to zero. I traced the negative thicknesses to two sources:
(1) Several cells have negative H in input files. Gunter will fix these.
(2) Sometimes, when calving extends into upstream interior cells, the calving thickness
exceeds the current ice thickness, resulting in negative H. I added a line of code
to limit the calving thickness to the current thickness. I also added a mask to prevent
calving in upstream interior cells that are not calving-eligible (e.g., grounded cells
in the case that only floating cells are allowed to calve).
I also modified subroutine glissade_cleanup_tiny_thickness to
(1) remove ice with tiny negative thickness (magnitude < eps11), which can arise from
roundoff errors in subgrid calving, and
(1) abort if any ice has negative thickness of larger magnitude (> eps11).
This should prevent any ice with negative thickness from remaining at the end
of a timestep.
Answers change slightly in runs with subgrid calving.
This commit removes some logic related to the parameter H0_float, which was introduced several years ago to reduce basal melting in cells with thin ice, under the assumption that these cells likely were partial cells at the calving front. This logic has been superseded by the subgrid CF parameterization. This changes the value of f_float, a variable used as a weighting function when computing basin averages of bmlt_float and related variables. It now depends on f_ground_cell but not H0_float. In general, f_float is equal to the floating fraction in ice-filled cells that are at least partly floating, and is zero for ice-free ocean. Modifying f_float changes answers for the basal melt calibration added recently. I also added some optional diagnostics in module glissade_bmlt_float.
This commit makes several changes in the isostasy code in libglide:
* When running with an elastic lithosphere, the load updates are done slightly differently.
E.g., suppose nlith = 100. Previously the load was computed at initialization and
updated at steps 1, 101, 201, etc. Now it is computed at initialization and updated
(more logically) at steps 100, 200, etc. This leads to small answer changes
when isostasy is turned on.
* The option ASTHENOSPHERE_RELAXING is now call ASTHENOSPHERE_RELAXING_CONST.
This is in preparation for introducting a new option in which the relaxation time
is not constant but laterally varying.
* The constant relaxed_tau is now called tau_relax_const. Its default value is 3000 yr
instead of 4000 yr (since 3000 yr is more typical in the literature for regions
with slow relaxation).
* I removed an unnecessary mpi gather in subroutine calc_elastic, which updates
the elastic load.
* I changed 'whichrelaxed' to 'which_relaxed' and added an exact restart capability
for option which_relaxed = 2. For this option, the input topg is assumed to be
in equilibrium with the load, and relx is computed as topg + load.
On standard restarts (is_restart = 1), relx and load should be read from the restart file
rather than recomputed.
Hybrid restarts (is_restart = 2) are tricky, since the run from which we restart
may or may not have had active isostasy.
- If the previous run had active isostasy, then relx and load will be in the restart file.
The user should set which_relaxed = 0 for the hybrid run.
- If the previous run lacked active isostasy, there is usually no relx field in the
restart file. Then the user should set which_restart = 2, so that relx is computed
at initialization as topg + load.
* I added some runtime isostasy diagnostics, e.g. some calls to point_diag.
This commit adds three new modules in libglissade: glissade_isostasy, glissade_isostasy_elastic, and glissade_isostasy_kelvin. These modules are similar to the glide modules isostasy, isostasy_elastic, and isostasy_kelvin. The biggest difference is that the new glissade modules support a laterally varying relaxation time, as discussed by van Calcar et al. (TC, 2026). There are now three asthenosphere options: * ASTHENOSPHERE_FLUID = 0 * ASTHENOSPHERE_RELAXING_CONST = 1 * ASTHENOSPHERE_RELAXING_LATVAR = 2 To use the new option, the user should set asthenosphere = 2 in the config file and add an input file containing a 2D field called tau_relax. This option is supported for the glissade dycore but not the glide dycore. To use the old option with spatially uniform relaxation time, the user should set asthenosphere = 1 (as before) and set tau_relax_const (which used to be called relaxed_tau) in the config file. Other minor changes: - I moved subroutine glissade_isostasy_solve from glissade to glissade_isostasy. - I made some minor changes in the isostasy subroutines of glide_setup. - I changed some subroutine names in the main isostasy module. I verified that the new glissade code gives the same answers (BFB) as the old glide code. (For testing, I inserted a hack to thin the ice before calling the isostasy; I commented out the hack before committing the code.) I did not remove the glide modules. In principle, it should still be possible to run the glide dycore using the glide isostasy modules, though I haven't tested this. The new asthenosphere option 2 is still to be tested.
Subroutine calc_elastic in glissade_isostasy_elastic computes the load for each grid cell
by taking a distance-weighted sum of local load factors in a region surrounding the cell.
Since this region can extend over many grid cells on multiple tasks, the sum
requires global arrays.
Some time ago, I parallelized calc_elastic in a simple but inefficient way:
(1) Gather the load factors to a global array (load_factors_global) on the main task.
(2) Compute the load in a global array on the main task. The sum includes
global indices only.
(3) Scatter the global load array to each task.
This calculation is slow because the sum (2) requires loops over all cells
on the main task while the other tasks are idle.
I realized there is a more efficient way to compute the load:
(1) Gather the load factors to load_factors_global on the main task.
(2) Broadcast load_factors_global to all tasks. This requires a new subroutine,
broadcast_real8_2d, in the broadcast interface of parallel_mpi.F90.
(3) On each task, compute the load by summing over local cells. The sum includes
a mix of local indices (i,j) and global indices(m,n). No scatter is needed.
This calculation is much faster because each task computes the load over only
its locally owned cells.
Note: An MPI gather followed by a broadcast is equivalent to an MPI allgather.
However, a gather followed by a broadcast was simpler to code than an allgather.
I tested the new method in short runs on Derecho (AIS runs on an 8km grid,
256 cores, with isostasy called at every time step).
* With the new method, a 100-year run takes ~7 minutes, compared to 6 minutes
with isostasy turned off. Thus, the cost of the isostasy is modest compared to the
rest of the model, even when it is called much more often than it would be in practice.
* With the old method, a 100-year run takes about 4 hours. Thus, the isostasy is
about 250 times slower than with the new method. This makes sense, given that the
old method does all load sums on one task, while the new method is embarrassingly parallel.
For now, I left the old method in the code for comparison.
The user can switch between methods by changing the value of new_load_sum.
A future commit will remove the old method.
The previous commit added a new, more efficient method of computing the load in each grid cell for an elastic lithosphere. This commit removes the old method. I changed the name of lithosphere_period to load_update_interval. Since the update is more efficient than it was before, I shortened the default value to 10 yr. I also updated the isostasy comments near the top of glissade_isostasy.F90.
This commit fixes an issue I missed when merging in the main branch.
14 tasks
Katetc
reviewed
Sep 11, 2026
Katetc
left a comment
Katetc
left a comment
Contributor
There was a problem hiding this comment.
Looks good, Bill! Please address these comments. There were a few bug-looking issues that we caught today. But, I can start running this code through the aux_cism test suite on Monday. I started an issue in cism_wrapper for this, and added the suggested namelist updates from Claude there. Please have a look and let me know if you have any other needed updates or changes: ESCOMP/CISM-wrapper#146
| tot_area_float, & ! total area of floating ice (m^2) | ||
| area_cell, & ! cell area | ||
| tot_volume, & ! total ice volume (m^3) | ||
| tot_volume_above_flotation, & ! total ice volume above flotation (kg) |
Contributor
There was a problem hiding this comment.
I think the units for your tot_volume_above_flotation should be (m^3) not kg
Contributor
Author
There was a problem hiding this comment.
Yes, good catch. Fixed.
| ! Optionally, write output to a specific basin with an applied thermal forcing anomaly | ||
| if (main_task) then | ||
| nb = model%ocean_data%thermal_forcing_anomaly_basin | ||
| if (nb >= 1) then |
Contributor
There was a problem hiding this comment.
Claude thinks you should add a top bounds check here as well. Something like
if (nb >= 1 .and. nb <= model%ocean_data%nbasin) then This is probably a good idea in places with nb >= 1 checks
Contributor
Author
There was a problem hiding this comment.
Good idea. Added a top bounds check here, and also in glissade_lateral_melt.F90 and glissade_inversion.F90.
| err_dmass_dt = tot_dmass_dt - & | ||
| (tot_smb_flux + tot_bmb_flux + tot_calving_flux + tot_latmelt_flux + tot_removal_flux) | ||
|
|
||
| ! uncomment to convert total fluxes from kg/s to Gt/yr |
Contributor
There was a problem hiding this comment.
Erm. Let's add a "TO DO" comment here that is along the lines of "TO DO: Make diagnostic units more flexible in debug output using namelist options rather than commented out code"
Contributor
Author
There was a problem hiding this comment.
Actually, I had already fixed the issue with output units, but I forgot to delete the commented-out code. It's gone now.
|
|
||
| err_dmass_dt = tot_dmass_dt - (tot_smb_flux + tot_bmb_flux + tot_calving_flux) | ||
| err_dmass_dt = tot_dmass_dt - & | ||
| (tot_smb_flux + tot_bmb_flux + tot_calving_flux + tot_latmelt_flux + tot_removal_flux) |
Contributor
There was a problem hiding this comment.
Claude and I are both confused about why you didn't include tot_gl_flux here. Is it already included in one of the terms? If so, maybe add a comment about that.
Contributor
Author
There was a problem hiding this comment.
Added the following comment:
! Note: The total mass budget consists of the five terms above.
! The GL flux is not part of the mass budget, since any ice fluxed across the GL
! is not added or lost, but simply changes from grounded to floating.
| private :: init_rbel, rbel_ow, rbel_iw | ||
|
|
||
| logical, parameter :: verbose_isostasy = .false. ! if true, print diagnostic messages | ||
| logical, parameter :: verbose_elastic = .true. |
Contributor
There was a problem hiding this comment.
Is this a code change needed to turn on verbose? Or is there a namelist setting associated with it? Either way, maybe it shouldn't default to True. If it requires a code change, I suggest a "TODO" comment that this should be a namelist option.
Contributor
Author
There was a problem hiding this comment.
Generally, the verbose options are compile-time settings rather than namelist settings and are false by default. I set verbose_elastic = .true. during testing and forget to turn it off. It is now false.
At some point, we could add a 'verbose' section to the namelist and make these runtime settings. But I'd first like to combine or remove some of them (especially in glissade_velo_higher.F90) so there aren't so many.
| private | ||
| public :: glissade_isostasy_init, glissade_isostasy_solve, verbose_isostasy | ||
|
|
||
| logical, parameter :: verbose_isostasy = .true. |
Contributor
There was a problem hiding this comment.
Should default to False and doesn't need to be a parameter.
Contributor
Author
There was a problem hiding this comment.
It is now false. I went through the other glissade modules, changing verbose variables to false and removing the 'parameter' specification as needed.
| ! the user should set which_relaxed = RELAXED_TOPO_STANDARD instead. | ||
| if (.not.parallel_is_zero(model%isostasy%relx)) then | ||
| call write_log ('Do not set which_relaxed = RELAXED_TOPO_COMPUTE if relx is in the input file') | ||
| call write_log ('Either remove relx or set which_relaxed = RELAXED_TOPO_STANDARD', GM_FATAL) |
Contributor
There was a problem hiding this comment.
Do you really want RELAXED_TOPO_STANDARD here? Claude is convinced that's not a real option.
Contributor
Author
There was a problem hiding this comment.
Claude is correct; this should be RELAXED_TOPO_DEFAULT. Fixed.
| ! The Schoof law supports inverting for Cc while assuming a fixed relationship to Cp: | ||
| ! Cp = gamma*Cc^p | ||
| ! Some sliding laws support inversion for both Cp and Cc at basin scale, | ||
| ! with two distinct thickness targets. |
Contributor
There was a problem hiding this comment.
Claude is concerned about the fact that the code in this module no longer has two distinct thickness targets. It looks like they have been collapsed together or the marine-grounded ice is now combined with the land-grounded ice. If that is the intent, then can we update this comment?
Contributor
Author
There was a problem hiding this comment.
Claude is correct. I removed the comment. I also added some code in glide_setup.F90 to abort the code if the user tries to invert for both Cp and Cc independently.
| real(dp) :: coulomb_c_min = 1.0d-3 !> min value of coulomb_c, unitless | ||
|
|
||
| ! parameters for Schoof basal friction law | ||
| ! Typically, schoof_n (from Eq. 2 in School 2005) has the same value as powerlaw_m, but this is not required. |
Contributor
Author
There was a problem hiding this comment.
Thanks, I try to spell names correctly because (among other reasons) the person in question might look at the code someday. Fixed.
| long_name: effective pressure factor from ocean_p | ||
| data: data%basal_physics%f_effecpress_ocean_p | ||
| load: 1 | ||
|
|
Contributor
There was a problem hiding this comment.
Claude is concerned that f_effecpress_ocean_p is removed but effecpress_ocean_p is not added into glide_vars.def. Is that on purpose or an oversight?
Contributor
Author
There was a problem hiding this comment.
effecpress_ocean_p is a temporary variable that may or may not be used to determine the value of effecpress. For that reason, I don't think it should be in glide_vars.def. For the same reason, I don't think it needs to be part of the basal_physics derived type, so I removed it from the derived type and made it local to the subroutine.
Contributor
Author
|
I've addressed all the comments above. It's very helpful to have Claude find issues that otherwise might have stayed in the code and confused CISM users. I pushed the changes to the lipscomb/ismip7 branch on escomp. |
The commit contains a number of small fixes suggested by Kate after a code review using Claude. I removed some obsolete comments and corrected some others. Also, all the verbose variables are now false by default, and they are no longer declared as parameters. These fixes are not answer-changing.
Mostly cleaning up comments; no answer changes
This commit introduces a new config parameter, thermal_forcing_basin_max, in the ocean_data type. If thermal_forcing_basin_max > 0, then the basin-average TF term in the ISMIP6 quadratic nonlocal and nonlocal-slope schemes cannot exceed thermal_forcing_basin_max. The limiter uses an exponential function, so that the basin-average TF term asymptotes to the max value smoothly and gradually instead of abruptly. The basin-average TF term is a proxy for the friction velocity. The physical justification for a limiter is that in very warm basins, the friction velocity will not increase without limit, but will eventually asymptote because of drag and other factors. The current motivation for a limiter is that the quadratic nonlocal scheme appears to be overly sensitive to ocean warming, except for very low gamma0. If the basal melt rate increases linearly instead of quadratically with TF as TF becomes large, then the sensitivity will decrease for given gamma0. To apply the limiter, the user should set thermal_forcing_basin_max to a positive value in the [parameters] section of the config file. The default value is 0.0, in which case no limiting is applied. This commit is answer-changing only for thermal_forcing_basin_max > 0.
This commit fixes a couple of minor bugs in some code used to compute thck_effective when running with the new subgrid calving front scheme which_ho_calving_front = HO_CALVING_FRONT_SUBGRID_FLOAT_GROUND = 2. I inserted some interior_mask variables where needed, and also inserted a missing 'eus' in a computation of thck_effective. These fixes will change answers only when which_ho_calving_front = 2.
This file contains hidden or 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
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
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 branch contains new physics options and other changes to support standalone CISM simulations for ISMIP7:
an improved elastic lithosphere/relaxing asthenosphere (ELRA) isostasy model. The new version supports laterally varying relaxation times and is fully parallel. It has been tested with a map of Antarctic relaxation times from van Calcar et al. (TC, 2026).
code for calibrating sub-shelf melt rates to match observations and ocean circulation models, as required for ISMIP7. During model initialization, CISM can compute temperature corrections in each basin such that the basin-average value of bmlt_float is equal to a prescribed target.
a generalized Schoof sliding law. The new law distinguishes the exponent m in the power law from an exponent n that regulates the transition between power law and Coulomb behavior. Also, users can now invert for Cp and Cc simultaneously by assuming a functional relation between the two coefficients.
work toward a lateral melt option for Greenland simulations, based on Heiko Goelzer's implementation for NorESM. The lateral melt rate at marine-terminating ice fronts is computed as a function of thermal forcing and subglacial discharge. When using a subgrid calving-front parameterization, the lateral melt calculation reuses code that was developed for subgrid calving. While most of the code is in place, the lateral melt option still needs some work that will be included in a future PR.
a new subgrid calving front option. The existing option 1 allows only floating cells to be calving-front cells. The new option 2 allows either floating or marine-grounded cells to be calving-front cells, as desired when lateral melt is turned on.
a calving option that combines stress-based calving with thickness-based calving. This option gives the best results to date for Antarctic-wide physically-based calving.
a new real-valued calving mask and some additional diagnostics to support calvingMIP experiments.
modified mass budget diagnostics. Previously, there were three mass budget terms: SMB, BMB and calving. Now there are two more terms: lateral melt and ice removal. Ice cap removal, if turned on, goes in the last category.
some code reorganization. For instance, I moved some high-level solver subroutines from glissade.F90 to the physically relevant module, and I restructured the calving, bmlt_float, and isostasy modules.