Rpointers update source cism #2
Merged
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
This commit updates CISM to be python3 compatible and the test cases be bit for bit with the "python2" version of CISM with the exception of the ross test case which is not bit for bit for the creation in the netCDF input file. This is the last commit that will ensure bit for bit with code using python2. The reason the ross case is not bit for bit depends on how some computation is now performed in Python3 and updated package version. In particular, for this test I did the following: A. I made the number of decimal points be trimmed to 3 for the velocity variable in the python3 version B. I manually specified the option dtype=float to the call to numpy.sin in the python3 version Doing A and B lead to bit for bit results with python2 counterpart. However, doing the multiplication of the corrected instances in A and B (i.e., velocity*numpy.sin(azimuth, dtype=float)) did not lead to bit for bit results with the python2 counterpart.
This commit removes the bit-for-bit compatibility with python2. The reason is due to the precision python3 writes to the config file which now contains more digits than python2. This led to minor differences in the outputs. This commit should be used as a new baseline for future code development comparison.
This commit contains additional changes towards adding flexibility for future python release. In particular the option parser (which is deprecated since version 3.2) is replaced by argument parser. In the ISMIP-HOM test case, the plotISMIP_HOM.py had to be slightly modified: a. The "reduce" function now needs to be imported from the "functools" module. b. The "map" function now returns an object as opposed to a list. While this had no big influences in this code, it did in the "read" function defined in this script and it had to be replaced by a list comprehension on line 222. In addition it seems that the map function in python2 could skip empty lines when reading through a text file. It's no longer the case in python3 and I had to add a check on line 219 to avoid an error.
Leguy/update toward cism3
The Intel compiler has been taking a long time to build glide_io.F90.
Some testing showed that this is because 'use glide_types' and three other use statements
were repeated in each of a large number of accessor subroutines.
For some reason, this hasn't been an issue on the Gnu compiler.
I modified ../utils/generate_ncvars.py and ../libglimmer/ncdf_template.in
so that these 'use' statements appear only at the top of the glide_io module.
With this change, the Intel build time on derecho with 8 cores ('make -j 8')
decreases from about 4:40 to 1:10, a welcome improvement.
Removed redundant 'use' statements that slowed the Intel build
Ktc/fixin slap
Katec/derecho update
Cleaning up some comments, e.g., changing 'Glimmer' to 'CISM'
Until now, the exponent n in the Glen flow law has been set in glimmer_physcon.F90
as an integer parameter, gn = 3.
With this commit, n is renamed n_glen (a more greppable name) for use in Glissade.
It is declared in glimmer_physcon.F90 as real(dp) with default value 3.0d0.
Since n_glen is no longer a parameter, it can now be set in the config file like other
physical 'constants' (e.g., rhoi and rhoo) that are not truly constant, but can
take different values in different models or benchmarking experiments.
To avoid changing answers in old Glide code, I retained the integer parameter gn = 3 in glimmer_physcon.F90.
This parameter is now used only in the Glide dycore (e.g., glide_velo.F90).
In Glissade, I replaced gn with n_glen in several places:
(1) In subroutine glissade_interior_dissipation_sia, the dissipation factor includes n_glen.
Note: n_glen does not appear explicitly in the 1st-order dissipation, which is proportional to tau_eff^2/efvs.
(2) In glissade_velo_sia, n_glen appears in the vertical integral for the velocity.
(3) In glissade_velo_higher, flwafact = flwa**(-1./n_glen) where flwa = A.
(4) In glissade_velo_higher, the exponent p_effstr = (1.d0 - n_glen)/n_glen is used
to compute effective_viscosity for BP, DIVA, or L1L2.
(5) In glissade_velo_higher, subroutine compute_3d_velocity_l1l2 has a factor proportional to
tau**((n_glen-1.)/2.) in the vertical integral.
For (1) and (2), n_glen was previously assumed to be an integer. Switching it to real(dp)
is answer-changing at the machine roundoff level for runs using the local SIA solver
in glissade_velo_sia.F90.
For (3), (4) and (5), n_glen replaces the equivalent expression real(gn,dp).
For this reason, answers are BFB when using the BP, DIVA or L1L2 solver.
In glissade_basal_traction, n_glen appeared in the expression for beta in the Coulomb sliding option,
HO_BABC_COULOMB_FRICTION. Here, I replaced n_glen with powerlaw_m (also = 3.0d0 by default)
to be consistent with the expressions for beta in the Schoof and Tsai laws.
In glimmer_scales, Glen's n appears in the expressions for the scaling parameters vis0 and vis_scale,
Here, I defined a local integer parameter gn = 3 so that the scales are unchanged.
It is now possible for the user to specify arbitrary values of n_glen in tests such as the slab test.
Another minor change: I added some logic to the subroutine that computes L1L2 velocities.
For which_ho_efvs = 2 = HO_EFVS_NONLINEAR, the effective viscosity (efvs) is computed from
the effective strain rate using an equation from Perego et al. (2012).
But for option 0 (efvs = constant) and option 1 (efvs = multiple of flow factor),
this strain rate equation in the code does not apply. For these options, I added an
alternative that computes velocity in terms of the strain-rate-independent efvs.
This allows us to use L1L2 for problems with constant efvs (e.g., the slab test).
The code was calling subroutine init_isostasy with isostasy = 0 = ISOSTASY_NONE. This subroutine is now called only if isostasy = 1 = ISOSTASY_COMPUTE.
I modified glissade.F90 to abort cleanly with a call to glide_finalise after an advective CFL error. This is done only when the user does *not* specify adaptive subcycling. The clean abort allows the new slabStability script to keep going, launching a new run with a shorter timestep. In subroutine glissade_flow_factor of glissade_therm.F90, I removed the FLWA_INPUT option (option 3 of whichflwa). This option is redundant with option 0, FLWA_CONST_FLWA, in which the user can set default_flwa in the parameters section of the config file, and otherwise CISM uses default_flwa = 1.0e-16 Pa^-n yr^-1. We probably should rename default_flwa to constant_flwa, but leaving it for now to avoid breaking config files in test cases. Note: This option was added by Matt Hoffman in 2014. I am unaware of test cases that use this option (most have flow_law = 0 or 2), but if there are any, we will need to fix them by switching to whichflwa = 0. In subroutine glissade_therm_driver of glissade_therm.F90, I increased the threshold for column energy conservation errors from 1.0d-8 to 1.0d-7 W/m^2. For very small timesteps of ~1.0e-6 yr, the smaller threshold can be exceeded because of roundoff errors. In subroutine glissade_check_cfl of glissade_transport.F90, I modified the fatal abort for large CFL violations (advective CFL > 10). Now, CISM aborts for CFL > 10 only when adaptive_cfl_threshold > 0, i.e. transport subcycling is enabled. This prevents excessive subcycling for egregious CFL violations. If adaptive_cfl_threshold = 0. (the default), i.e. transport subcycling is not enabled, then the code exits cleanly (with a call to glide_finalise) in glissade.F90 when advective CFL > 1. I added a TODO note to set the default value of geot (the geothermal heat flux) to 0 instead of 0.05 W/m^2. With the default value, simple prognostic tests like the dome are not mass-conserving. Not changing just yet because answers will change for the dome test.
This commit modifies the run and plot scripts for the Dukowicz slab test case, as described in Section 5 of this paper: J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34, https://doi.org/10.5194/tc-6-21-2012. The test case consists of an ice slab of uniform thickness moving down an inclined plane by a combination of sliding and shearing. Analytic Stokes and first-order velocity solutions exist for all values of Glen's exponent n >= 1. The solutions for n = 1 are derived in Dukowicz (2012), and solutions for n > 1 are derived in an unpublished manuscript by Dukowicz (2013). These solutions can be compared to those simulated by CISM. The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman with support for n = 1. They came with warnings that the test is not supported. The test is now supported, and the scripts include some new features: * The user may specify any value of n >= 1 (not necessarily an integer). The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant). * Physics parameters are no longer hard-coded. The user can enter the ice thickness, beta, viscosity coefficient (mu_n), and slope angle (theta) on the command line. * The user can specify time parameters dt (the dynamic time step) and nt (number of steps). The previous version did not support transient runs. * The user can specify a small thickness perturbation dh, which is added to the initial uniform thickness via random sampling from a Gaussian distribution. The perturbation will grow or decay, depending on the solver stability for given dx and dt. For n = 1, the viscosity coefficient mu_1 has a default value of 1.e6 Pa yr in the relation mu = mu_1 * eps((1-n)/n), where eps is the effective strain rate. For n > 1, the user can specify a coefficient mu_n; otherwise the run script computes mu_n such that the basal and surface speeds are nearly the same as for an n = 1 case with the mu_1 = 1.e6 Pa yr and the same values of thickness, beta, and theta. (Note: There is a subtle difference between the Dukowicz and CISM definitions of the effective strain rate; the Dukowicz value is twice as large. Later, it might be helpful to make the Dukowicz convention consistent with CISM.) I modified the plotting script, plotSlab.py, to look in the config and output files for physics parameters that are no longer hardwired. I slightly modified the analytic formulas to give the correct solution for non-integer n. This script creates two plots. The first plot shows excellent agreement between higher-order CISM solutions and the analytic solution for small values of the slope angle theta. For steep slopes, the answers diverge as expected. For the second plot, the extent of the y-axis is wrong. This remains to be fixed. I also added a new script, stabilitySlab.py, to carry out stability tests as described in: Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance of depth-integrated ice-dynamics solvers, to be submitted to The Cryosphere. The idea is that for a given set of physics parameters and stress-balance approximation (DIVA, L1L2, etc.), the script launches multiple CISM runs at a range of grid resolutions. At each grid resolution, the script determines the maximum stable time step. A run is deemed stable when the standard deviation of an initial small thickness perturbation is reduced over the course of 100 time steps. A run is unstable if the standard deviation increases or if the model aborts (usually with a CFL violation). I have run the stability script for several solvers (DIVA, L1L2, SIA, SSA) for each of two physical cases: one with thick shearing ice and one with thin sliding ice. Each suite of experiments runs in a few minutes on 4 Macbook cores for solvers other than BP. As expected, DIVA and SSA are much more stable than L1L2 and SIA. This and the previous commit correspond to the CISM code and scripts used for the initial submission by Robinson et al. (2021).
Rewrote the slab README file to describe the new command line options for runSlab.py, and the new script stabilitySlab.py.
This commit introduces new algorithms for computing the 3D velocity field in the L1L2 solver. Also, I fixed a bug so that the L1L2 solver works correctly for n_glen = 1. The goal was to improve the stability of the L1L2 solver in slab tests. At fine grid resolution (< 200 m), the solver is unstable even for very small dt, instead of following the stability limits of an SIA solver (as is the case for Yelmo). The default method (method 1) for the 3D velocity computes most quantities on the staggered grid and is unchanged. There are two new methods: * Method 2 follows the Yelmo algorithm as closely as possible, with some quantities on cell edges and some on vertices. Unlike Yelmo, there is necessarily an averaging of uvel and vvel from edges to corners at the end, since CISM assumes B-grid velocities. * Method 3 moves all intermediate calculations to cell edges, with a final averaging to vertices at the end. To support the new methods, there is a new subroutine, glissade_average_to_edges, in module glissade_grid_operators. Both new methods give results similar to method 1. Neither improves stability at very high resolution. This suggests that the stability of Yelmo might result from the overall C-grid velocity scheme rather than details of the 3D velocity computation. For each method, there is now an option to exclude the membrane stresses in the tau_xz and tau_yz terms that appear in the vertical integration factor. When these stresses are excluded, the stability curve for L1L2 solver is close to that of an SIA solver, like Yelmo.
This commit includeds the following changes: * I fixed some typos in the README file for the slab case. * In runSlab.py, I added the option to use the local SIA solver (contained in module glissade_velo_sia.F90). * I added some code which can apply an initial sinusoidal perturbation of wavelength 2*dx, instead of a random Gaussian perturbation. This is useful for getting reproducible results.
This commit adds a hybrid velocity solver as described by Robinson et al. (TC, 2021). The solver first computes SIA velocities using the local SIA solver (which_ho_approx = -1) with zero basal slip, then computes SSA velocities using the Glissade higher-order SSA solver (which_ho_approx = 1), and finally adds the two solutions. The new logic is in module glissade_velo. To use the hybrid solver, set which_ho_approx = HO_APPROX_HYBRID = 5 in the config file. For the slab test, the hybrid solver has the expected stability properties, aligned with SSA at coarse resolution and SIA at fine resolution. Answers are as expected for a dome test.
This commit streamlines the 3D velocity subroutine for L1L2 and makes a small but important change in the algorithm. The change is to evaluate the membrane stresses in tau_xz and tau_yz at layer midpoints instead of lower layer boundaries, consistent with where the SIA terms are evaluated. I found that a small vertical offset can disrupt the balance between these two terms and make the L1L2 solver unstable for the slab problem at resolutions finer than ~200 m. With the fix, the L1L2 stability curve is parallel to the SIA stability curve at fine resolution, with L1L2 being slightly more stable than SIA. I also removed the alternate L1L2 discretization methods that were added in a recent commit. The alternate strategies evaluated some terms at edges instead of vertices, and did not change the results significantly. With this commit, all terms are evaluated at either cell centers or vertices. This is the code version used for the slab tests shown for CISM in Section 3 of the paper by Robinson, Goldberg & Lipscomb (2021, TC). For these tests, I temporarily changed the energy conservation criterion in glissade_therm.F90 to avoid false non-conservation alarms when running at very fine resolution. See the comments in that module. This commit is answer-changing for L1L2, in a good way.
This is the first of several commits to clean up the calving module, implement a new subgrid calving front scheme, and remove deprecated options. With this commit, the functionality of subroutine glissade_effective_calving_thck has been moved to subroutine glissade_calving_front_mask. In later commits, I will remove deprecated code from glissade_calving_front_mask.
This commit removes a complicated way of computing the threshold for thickness-based calving.
Some time ago, I added code to specify the threshold as a 2D field called thck_calving_threshold.
Now, the threshold is simply a scalar equal to calving_minthck.
The thinning rate at the CF if given by
dH_eff/dt = -(Hc_min - H_eff) / tau_c,
where Hc_min = calving_minthck, H_eff = thck_effective (constrained by H_eff >= H),
and tau_c = calving_timescale. Then dH_eff/dt is converted to dH/dt:
dH/dt = dH_eff/dt * H/H_eff
Since cells with H < Hc_min < H_eff do not calve, the CF can advance.
Calving options 6, 7, and 8 (thickness-based, eigencalving, and damage-based) all now
use the new method of computing H_eff. For each scheme, unprotected ice (i.e., cells
beyond a protected row of partical CF cells with H < H_eff) are removed after transport.
Calving is iterated until dt_calving = 0 in all CF cells.
A next step is to remove the deprecated subgrid CF scheme and replace it with the scheme
developed to support these calving options.
This commit removes the variables thck_calving_front and active_ice_mask, and in general all the logic associated with the old subgrid CF scheme. In that scheme, CF cells with thck < thck_calving_front were identified as dynamically inactive. This could lead to flickering and oscillations when a cell changed from inactive to active (or vice versa), with large changes in the local force balance. An upcoming commit will replace the old scheme with the new scheme developed for damage-based calving, eigencalving, and thickness-based calving. Most elements of the new scheme are now in place, but are controlled by the choice of calving scheme rather than the value of which_ho_calving_front.
This commit changes the subgrid calving front scheme associated with the
config option which_ho_calving_front = 1 = HO_CALVING_FRONT_SUBGRID.
In the old scheme, partially filled cells at the calving front were dynamically inactive.
This led to undesirable flickering near the CF.
In the new scheme, these cells are active, like all other cells with thck > thklim.
They have an effective thickness, thck_effective, which is now part of the calving derived type.
For most calving-front cells, thck_effective > thck.
The following calculations use thck_effective in place of thck in partial CF cells:
(1) several calving options: thickness-based calving, eigencalving, damage-based calving.
These options now require which_ho_calving_front = 1; else the code aborts.
(2) SMB and BMB calculations: An effective area fraction, thck/thck_effective,
is used to reduce SMB and BMB in partly filled CF cells.
(3) HO velocity solver: The solver replaces thck with thck_effective in partial CF cells,
for more accurate computation of gradients and lateral pressure at the CF.
This commit implements (1) only, so as to give the same answers (BFB) in tests
of new calving schemes with MISMIP+ geometry.
A follow-up commit will implement (2) and (3), which will change answers for most runs
using the new subgrid CF scheme.
With this commit, the SMB/BMB calculation uses effective_areafrac as computed using the new subgrid CF scheme, when that scheme is enabled. To make the calculation more accurate, I moved some code that removes unprotected ice beyond the subgrid CF. This ice used to be removed during the calving solve, but now is removed between the transport solve and the SMB/BMB solve. Accordingly, the field model%calving%calving_thck is now initialized to zero at the start of the timestep, instead of at the start of the calving solve. This is an answer-changing commit for runs with the new subgrid CF scheme.
For the new subgrid CF scheme, the velocity solver now sets H to the effective thickness (calving%thck_effective) instead of the state-variable thickness (geometry%thck) in calving-front cells. The goal is to make the velocity solver smoother and more accurate by applying a large lateral pressure term to partial CF cells, instead of a large pressure gradient term between full cells and partial CF cells. In practice, the velocities near the CF are not very different, at least for the MISMIP+ calving tests. This commit changes answers for runs with the new CF scheme.
Various cleanup for the calving module and related subroutines: * Moved some stress tensor and strain rate calculations from glissade_diagnostic_variable_solve to new subroutines in the calving module. * Created a config variable called dthck_dx_cf, replacing a hardwired constant in glissade_masks. This is the thickness gradient used to compute thck_effective at the subgrid calving front. The default value is still 0.002. * Removed the seldom-if-ever used Huybrechts calving option * Replaced some diagnostic print loops with calls to subroutine point_diag * Made indentation consistent in subroutine calls
Added a parallel_halo update for thck_effective in subroutine glissade_calving_front_mask. Without this update, there are many symmetry errors in the velocity solver. With the update, the results for multiple cores are the same (within roundoff) as the results on a single core.
This commit includes several changes to damage-based calving and eigencalving, to obtain more realistic whole-Antarctic calving simulations: * Changed the parallel logic for calving iterations using a new subroutine called parallel_reduce_log_or (true globally iff true on any processor). Added a companion subroutine, parallel_reduce_log_and (true globally iff true on every processor). * Changed the logic for protecting calving-front cells during transport when using the subgrid CF scheme. Cells are now protected if they have "full" neighbors. Cells are full if in the ice interior, or at the CF with a thickness close to that of the thickest interior neighbor. * Added calls to thickness-based calving as a follow-up to damage-based calving and eigencalving. Calving of thin ice helps prevent the formation of extensive thin shelves. * Simplified the lateral_rate_factor calculation. Now the factor is simply set to the number of ocean neighbors. A whole-Antarctic run with damage-based calving looks fairly good after 1 ka. Damage accumulation leads to reasonable CF locations for the Ross, Filchner, Amery, and Larsen C shelves. The Ronne CF is a bit too far advanced. Thickness-based calving (minthck = 100 m) largely prevents the formation of thin, extensive shelves, with the exception of Thwaites, Shackleton and Moscow-Totten. A whole-Antarctic run with eigencalving also looks fairly good after 1 ka.
This commit adds subroutines parallel_reduce_log_or and parallel_reduce_log_and in module parallel_slap. These subroutines were recently added in parallel_mpi. To support the serial build, I moved subroutine glissade_test_comm_row_col to module parallel_mpi, changing 'glissade' to 'parallel' in the name. A new subroutine of the same name in module parallel_slap prints an error message and aborts. By default, this subroutine is not called. It was written to test the row and column communicators, which are supported for parallel runs only. The serial code now builds cleanly on Cheyenne.
This commit adds a calving option, whichcalving = CF_ADVANCE_RETREAT_RATE = 9.
This option will be used for CalvingMIP:
https://github.com/JRowanJordan/CalvingMIP/wiki
Some of the CalvingMIP experiments require the calving front to be stationary
or to advance and retreat at a specified rate.
This work is still in progress. I've done the following so far:
* Introduced the new calving option, which is implemented with the subroutine
calving_front_advance_retreat in glissade_calving.F90. This subroutine
computes the ice volume that must be removed to shift the CF at the desired rate.
* Added another subroutine, glissade_redistribute_unprotected_ice, to support
the new option. This subroutine looks for unprotected ice which is out ahead
of the CF after the ice transport. Instead of calving unprotected cells
(which has been the default behavior), it moves this ice back upstream to the
partial cells (H < H_effective) at the pre-transport CF.
The ice in these cells can then calve according to the prescribed rate.
* Got the new option working in a MISMIP+ experiment with flow in the x direction.
This is simpler than the 2D flow in the calvingMIP experiments.
* Made several logic changes to get calvingMIP spinups to work. These include:
- When running with a uniform constant C_p or C_c, initialize the C_p and C_c
arrays on restart. These arrays were being erroneously set to zero on restart.
- Recompute ice and calving masks (including effective_areafrac) after the transport
ONLY if running with which_ho_calving_front = 0. When running with a subgrid CF,
the logic now requires the calving scheme to use the pre-transport calving masks.
- When using the new calving option 9, do NOT remove unprotected ice immediately
after the transport. This ice is removed by the calving scheme instead.
Later, this behavior may be extended to other subgrid CF options such as
eigencalving and damage-based calving.
* Changed the little-used Huybrechts calving scheme to option 10.
Future work:
* Make the CF advance/retreat rate a config parameter instead of hardwiring it.
* Modify the ice redistribution scheme to put ice in more than one upstream grid cell.
* Simplify the calving logic by removing the iterations; instead, just calve a certain ice volume
in the initial CF cell and (if necessary) its upstream neighbors.
* Rework the other subgrid CF schemes (thickness-based, damage-based, and eigencalving)
in a similar way.
With this commit, the new calving option (CF_ADVANCE_RETREAT_RATE = 9)
works for the circular and Thule domains.
To enable symmetric retreat with a 2D flow pattern (and not just 1D flow as
in MISMIP+), I changed the redistribution logic at the calving front.
Instead of moving ice to the thickest upstream neighbor, we distribute ice to each
upstream neighbor in proportion to the flux received from that neighbor.
Fluxes are estimated using subroutine glissade_input_fluxes.
These fluxes are not identical to those from incremental remapping, but are very close.
(I copied this subroutine from a glacier branch that will likely be merged to main
before this one. In case of conflicts, we should defer to this newer version.)
There is some similar post-calving logic. If a column melts entirely before
all the thinning is used up, the thinning is extended to upstream neighbors
that were full (H = H_eff) at the start of the timestep.
If cells at the CF are over-full (H > H_eff) after calving, then the excess ice
is distributed to downstream neighbers, advancing the CF.
With the new logic, it is no longer necessary to iterate the calving for option 9,
so I removed the iterations. Later, it might be possible to remove iteration
from other subgrid calving options (eigencalving, etc.).
To support the prescribed advance and retreat rates in the CalvingMIP experiments,
I added two config parameters: cf_advance_retreat_amplitude and cf_advance_retreat_period.
Let A = amplitude, tau = period, t = time elapsed since tstart.
Then the target rate W for the calving front is
W = A * sin(2*pi*t/tau).
For A > 0, the CF will advance until t = tau/2 and then retreat.
For A < 0, the CF will retreat until t = tau/2 and then advance.
I ran CalvingMIP experiments 1 and 2 for the circular shelf and experiments 3 and 4 for Thule.
The results look good except for the CF readvance in the second half of Expt. 4.
That will take a bit more coding.
Two changes in masks:
(1) I changed the way H_eff = thck_effective is computed in calving-front cells.
Previously, H_eff could be derived only from neighboring cells that are floating.
Now, H_eff can also be derived from neighboring cells that are marine-grounded.
This gives reasonable values for H_eff in CF cells adjacent to marine-grounded cells.
(2) The call to glissade_redistribute_unprotected_ice now follows the use of
the calving mask to remove ice from masked-out cells. This prevents the
redistribution from generating a sawtooth calving front.
I also changed the default dthck_dx_cf to 0.
This is a parameter used for the subgrid CF scheme; it might be deprecated later.
Some time ago, I wrote thickness-based, eigenvalue-based, and damage-based calving schemes
for testing with MISMIP+. These were the three calving schemes developed for use with
the subgrid CF option.
Recently, for CalvingMIP, I wrote a calving scheme with prescribed advance or retreat rates.
This scheme has some new logic that seems equally applicable for the other three schemes.
With this commit, the other three schemes have been rewritten. The logic of the four
calving schemes using a subgrid CF is now as follows:
(1) Compute some masks related to calving.
(2) Where ice has been transported downstream from a partial CF cell to cells beyond the CF,
move it back upstream.
(3) For CF cells, compute the thinning rate that will give the desired rate of CF advance or retreat. The thinning rate is computed from the lateral calving rate, which depends on the calving law.
(4) Apply the calving-derived thinning. If a full column is removed at the CF,
do additional thinning upstream.
(5) Where H > H_effective in CF cells, set H = H_effective and move the extra ice downstream,
advancing the CF.
Only step (3) depends on the specific calving law. Each of the four subgrid CF schemes
uses the same code for (1), (2), (4), and (5), with different subroutiness only to compute
the calving-derived thinning rate. As a result, some redundant code has been removed.
Other improvements: Calving no longer needs to be iterated multiple times when
a full column is removed at the CF. Also eigencalving and damage-based calving
no longer required a subsequent call to thickness-based calving for cleanup.
I tested each scheme using the radially symmetric geometry of CalvingMIP Experiment 2.
Results look generally reasonable, although in some cases the radial symmetry is lost.
(E.g., calving is favored or disfavored along grid diagonals in comparison to the x and y axes.)
I will follow up with tests using Experiment 4, the Thule domain.
This commit modifies the eigencalving scheme and introduces a new stress-based calving law. The eigencalving scheme (calving option 7) now is more similar to Levermann et al. (2012). The calving rate is proportional to the geometric mean of the eigenvalues of the horizontal strain rate tensor, provided that both eigenvalues are positive. In practice, this means that calving is restricted to regions where the second (across-flow) eigenvalue is positive, since the first eigenvalue is almost always positive. In partial CF cells, the strain rate is computed as a weighted average of the value in the partial cell and an upstream value in a full cell. For this scheme, the user supplies a config parameter called eigenconstant, which is analogous to K from Levermann et al. (2012) but is scaled differently. I got reasonable results for a value of ~5.e5 m. The stress-based scheme (new calving option 8) resembles the scheme of Morlighem et al. (2016) and borrows some features of the previous eigencalving scheme. The calving rate is proportional to |v| * stress_eff/stress_threshold, where |v| is the ice speed at the CF and stress_eff is an effective stress formed as a weighted sum of the two eigenvalues. The resulting CF is located where stress_eff ~ stress_threshold. Weighing the second eigenvalue more strongly gives results similar to eigencalving. For this scheme, the user supplies the stress threshold along with parameters tau_eigenconstant1 and tau_eigenconstant2, which are weights for the two eigenvalues. These two schemes closely resemble the EC and VM schemes investigated by Milner et al. (2023, TC) for steady-state Antarctic simulations. EC and VM were the top-performing schemes out of the four schemes studied. I fixed a calving bug in glissade_velo_higher.F90. For the subgrid CF scheme, the velocity solver replaces thck with thck_effective. The bug fix is to recompute usrf to be consistent with the replacement thickness. In subroutine glissade_add_smb of glissade_transport.F90, I added some logic to turn off surface ablation in cells with ocean_mask = 1. This fixes some column conservation errors. I'm not clear on why we didn't see these errors earlier. In subroutine glissade_calving_front_mask of glissade_masks.F90, I added an argument called thck_full. This field is used to determine when a CF cell has filled with ice, allowing the CF to advance. I added logic to protect ice-free ocean cells from ice removal if surrounded by three partial CF cells with ice. Previously, such cells tended to remain ice-free. I added effective_areafrac to the calving derived type. I made eps_eigen1 and eps_eigen2 loadable on restart. Still need to modify damage-based calving (now option 10) consistent with the eigencalving changes.
Got rid of some unneeded mask calls Added eps_eigen and tau_eigen to restart
CalvingMIP Experiments 2 and 4 have a prescribed, radially symmetric calving rate. In the previous code version, retreat was a bit faster along the N/S/E/W axes than along diagonals. This commit contains several changes that improve radial symmetry and give better agreement with the prescribed retreat and advance rates. In the calving module, there is a new subroutine, compute_calving_front_length, that estimates the CF length in each grid cell. The idea is to draw line segments between adjacent cells along the CF and compute the length of the segments. For example, if the CF lies along a row parallel to the x-axis, cf_length = dx, and if the CF lies along a column parallel to the y-axis, cf_length = dy. If the CF is parallel to a diagonal (e.g., y = x), then cf_length = sqrt(dx^2 + dy^2). For some cells, cf_length is the sum of an along-axis segment and a diagonal segment. There is logic to ensure that each CF cell is connected to at most 2 CF neighbors. I tried a new method for restricting flow from partial CF cells to unprotected ocean cells. The old method is based on subroutine redistribute_unprotected_ice; any ice that flows into unprotected cells is pushed back upstream to the cells where it came from, using a simplified method to estimate fluxes. I tried doing this in a more exact way, by passing edgemask_e and edgemask_n to glissade_horizontal_remap. These are integer masks with a value of 1 for all edges across which transport is allowed, and otherwise a value of 0. I set edgemask = 0 for certain edges of partial CF and ocean cells when using the subgrid CF scheme. I got the method partially working, but in more complex configurations (e.g., with retreat in CalvingMIP Experiment 4), I couldn't figure out how to avoid negative areas. The problem is that for some patterns of masking, we can have outflow across an open edge without having any inflow across the adjacent masked edge. So I reverted back to the old method. For now, there is a logical parameter 'use_edgemasks' that can be set to true to activate the new edgemask method. I may remove this option later, unless I somehow get it fully working. Minor changes: * Fixed some initialization logic so that Cp is not set to 0 on restart when running with constant Cc or Cp * Added effective_areafrac to I/O variables * Replaced some diagnostic prints with calls to subroutine point_diag With these changes, calvingMIP results look good for Experiments 1 to 4. It remains to test these changes for physically-based calving laws.
Katetc/hma glaciers4, Replace Lipscomb/hma_glaciers4 PR
This commit makes several changes in calving logic, mostly to improve whole-Antarctic runs with stress-based calving: * The user can set a minimum value for thck_effective by passing calving_minthck to subroutine glissade_calving_front_mask. This helps prevent the propagation of ice shelves with unrealistically thin ice and low calving rates. Recall that the rate of volume loss at the CF is proportional to thck_effective. * When computing thck_effective, the upstream thickness now can come from neighboring land-based cells, not just floating cells. Since thck_effective is capped at the flotation thickness, thick land-based ice will not result in huge values of thck_effective. * I modified the CF length computation to account for CF cells with 0 or 1 CF neighbors. In these cells, the CF length is now at least as large as dx. * The calving-front length is now passed into the subroutines for stress-based calving, damage-based calving, and eigencalving. The rate of volume loss is proportional to the CF length, consistent with CalvingMIP. Specifically, the rate of loss is equal to Cr * dt * H_eff * cf_length / (dx*dy), where Cr is the lateral calving rate (m/yr) computed from the stress or the strain rate. * I modified the damage-based calving scheme to be more consistent with stress-based calving. However, this scheme needs more work, since it tends to generate damage = 1 over a large region. * I removed the 'tooth factor' for the CalvingMIP scheme. The results are very similar when cf_length is computed in the standard way for protruding CF cells with 3 ocean neighbors, instead of increasing the effective cf_length. * Modified some CF location logic for the CalvingMIP axes With these changes, the computed CF locations look quite good for an AIS spin-up with stress-based calving. Spurious CF advance is fairly minimal. I was seeing too much retreat for the Amery shelf and some shelves in Queen Maud Land, but this is reduced by having a stress threshold of 75K instead of 50K. The higher threshold also reduces spurious retreat for Ross and Filchner-Ronne.
The option numbers for calving options 6 to 10 are now as follows: integer, parameter :: CF_ADVANCE_RETREAT_RATE = 6 integer, parameter :: CALVING_THCK_THRESHOLD = 7 integer, parameter :: CALVING_STRESS = 8 integer, parameter :: EIGEN_CALVING = 9 integer, parameter :: CALVING_DAMAGE = 10 CalvingMIP experiments 1 to 4 now use option 6. Option 8 for stress-based calving is unchanged. This is the code version used for the CalvingMIP experiments that Gunter submitted in late July 2024.
The calving-front length is now passed to the thickness-based calving scheme, consistent with the other calving laws that use the subgrid CF scheme. I also created a new config parameter called thck_effective_min; this is a lower limit on thck_effective for CF cells. Earlier, I called that parameter 'calving_minthck', but calving_minthck has a distinct meaning for thickness-based calving, so it's clearer to have two different parameters. Also removed thck_full as an output argument in subroutine glissade_calving_front_mask. For the CalvingMIP experiments submitted in late July 2024, this commit gives the same answers (BFB) as the previous commit, 'Changed option numbers...'
In an earlier commit on this branch, I added optional edge masks to the remapping transport. The goal was to prevent ice from leaving partial CF cells. However, the masks did not work, since blocking flow across some edges and not others leads to negative ice thickness in many simulations. This commit removes the edgemask logic. I also turned off the verbose_calving option. This commit is BFB for calvingMIP and other simulations with calving.
In an earlier commit, I added a variable 'time' to the calving derived type. This was used for CalvingMIP experiments with a prescribed advance/retreat rate, anticipating two different rates during a two-part run (e.g., Experiment 4). However, the variable was applied in a way that breaks exact restart for other CalvingMIP experiments. This commit removes it, replacing it with the regular time variable (model%numerics%time). This fixes exact restart. CalvingMIP results are BFB. In particular, Experiment 4 results are unchanged, since the advance rate during the second phase (after year 500) is time-independent.
The new python setup script for CalvingMIP (to be checked in soon) creates config files based on a template, using the configparser module to tailor each file to a given experiment. The configparser requires each new section to have a unique name. This means we cannot have multiple sections called 'CF output'. Instead, the parser adds sections called 'CF output1', 'CF output2', etc. I made a couple of changes in CISM to handle these sections correctly: * I added logic such that any section with a name that includes 'CF output' as a substring is treated as a CF output file. * I removed logic that requires 'CF restart' to be the final section of the config file. It can now be followed by sections added by the parser. With these changes, the setup script is able to set up config files that do not need any additional user edits.
I added a directory called calvingMIP in the tests directory. It includes three files: * calvingMIP.Setup.py, a Python setup script which creates the input netCDF files for the circular and Thule domains, and then creates seven subdirectories: SpinupCircular, SpinupThule, Experiment1, Experiment2, Experiment3, Experiment4, and Experiment5. Each subdirectory includes a config file suitable for the given experiment. * calvingMIP.config.template, a template file for creating the experiment-specific config files * a README file with instructions For more details, see the README file and the CalvingMIP wiki: https://github.com/JRowanJordan/CalvingMIP/wiki I've run multiple tests to make sure the experiments are set up correctly and yield (nearly) the same results as the recent CalvingMIP submission. Experiments 1 to 5 now start in hybrid restart mode (restart = 2) instead of standard restart (restart = 1). This leads to small differences in results (typically around the 8th sig fig after 1000 years).
This commit adds a few lines needed to run the mismip+ tests.
New calving schemes for CalvingMIP and Antarctic simulations
12 tasks
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.
Updates to CISM to use new rpointers in NorESM3_0_alpha01