Lipscomb/smb cleanup.v2 #75
Open
Conversation
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
This commit modifies an inversion option used to optimize floating ice thickness at the basin scale. The option was previously called which_ho_bmlt_basin, but now is called which_ho_deltaT_basin (which is more descriptive, since we are inverting for temperature corrections, not melt rates). There are four choices: 0: deltaT_ocn = 0 in each basin 1: invert for a single value of deltaT_ocn in each basin 2: read deltaT_ocn from an external file for each basin 3: prescribe deltaT_ocn in each basin using ISMIP6 values The inversion calculation now uses the same inversion parameters as deltaT_ocn inversion: deltaT_ocn_timescale, deltaT_ocn_thck_scale, deltaT_ocn_temp_scale, and deltaT_ocn_relax (Note: deltaT_ocn_relax = 0 in nearly all cases; it is now a parameter instead of a 2D field.) Thus, inversion for deltaT_basin is now logically equivalent to deltaT_ocn inversion, except that the thickness terms are basin-scale averages instead of local values, and we compute a basin-scale correction instead of a 2D field. The parameters bmlt_basin_timescale and dbmlt_dtemp_scale are no longer used. Also, in subroutine ismip6_bmlt float, I removed some logic to limit deltaT_ocn based on the current thermal forcing at the lower ice surface. When applied after the inversion was finished, this logic broke the assumption that deltaT_ocn should remain fixed. Removing the logic ensures that deltaT_ocn stays fixed in forward runs without inversion, and makes little difference in results.
This commit makes a number of changes to inversion methods. The biggest change is a new option to invert for flow_enhancement_factor E based on a surface speed target, as suggested by Tim van den Akker. I removed the option to invert for E based on a thickness target. The options to invert for C_c, C_p and deltaT_ocn based on a thickness target remain. I fixed a bug in subroutine ismip6_bmlt_float. This subroutine could change the value of deltaT_ocn even when that field was supposed to be held fixed. Thanks to Tim for pointing out the bug. For all thickness-based inversions, I changed the leading (dH/H0) term to (dH/H0)^n, where n is a user-prescribed exponent. The default is n = 1, but I have been running with n = 2 to give greater weight to larger errors dH. For the which_ho_flotation_function options, I removed the deprecated option 2. The new option 2 is the old option 3 (using a linear flotation function without extrapolating value to ice-free ocean cells). I created a new option 3, which computes the flotation function based on a field called topg_raised. This is a field in which the coarse-grid bed topography is higher (i.e., shallower) to better capture pinning points on the BedMachine grid. We are now testing various ways to compute topg_raised. I removed an old option to adjust topg on coarse grids based on the standard deviation of topg as computed on a finer grid. Other changes: * When inverting for C_c, we now initialize C_c to the same value everywhere (instead of starting with a low value for cells that are floating or ice-free ocean). * In the inversion for C_c and C_p, I removed two optional terms that weren't being used. * Changed the default timescales for all inversions to 200 years. * Replaced some verbose diagnostics with calls to the point_diag subroutine. * Added some parallel_globalindex diagnostics to the matrix checks in the velocity solver. * Set verbose_calving = F.
This commit includes more changes related to Antarctic spin-up. First, I added a config paramater called thermal_forcing_basin_min. The default value is zero. If set to a value greater than zero, this becomes the minimum value of thermal_forcing_basin applied in the ISMIP6 nonlocal and nonlocal-slope basal melting schemes. This parameter becomes important in cool basins where negative temperature corrections (in response to thin shelf ice) can create local values TF < 0. Assuming the basin-scale TF > 0, the product of local and basin-scale TF is negative, leading to basal freezing. In long spin-ups with inversion for deltaT_ocn based on a thickness target, many cells in the Ronne and Ross shelves benefit from modest basal freezing. This is consistent with observations (see, e.g., Fig. 5 of Lambert et al., 2023, TC) and reduces a persistent thin bias. Second, I added a subroutine called glissade_rms_error, which can be used to evaluate model biases relative to observations. In subroutine glide_write_diagnostics, this subroutine is now called to compute the RMSE for modeled ice thickness and surface speed. Current AIS spin-ups have an RMSE of about 35 m for thickness and 190 m/yr for speed. The thickness bias is quite low; the speed bias is larger because of substantial slow biases for PIG, Thwaites, and the Brunt ice shelf. Other changes: * I added an optional mask argument for the parallel_global_sum interface. This is useful for RMSE calculations. * I moved the parallel_global_sum_staggered interface from the PCG solver to parallel_mpi. * I added a function parallel_is_zero; this returns .true. when any value of a 2D field in the global domain is nonzero; else it returns .false. This is useful for checking whether a given field has been read in. * I added some point_diag calls in several modules as part of ongoing cleanup.
Added some diagnostics related to the ice age in (1) glissade_transport.F90, subroutine glissade_vertical_remap (2) glide_diagnostics.F90, subroutine glide_write_diag The ice_age is now incremented at the start of each time step, so that the maximum possible age at the end of the step is equal to the total model time.
This commit allows the user to add Laplacian smoothing terms when inverting
for basal friction coefficients (Cc or Cp), ocean temperature corrections (deltaT_ocn),
or the flow enhancement factor (E). The goal is to inhibit noisy fields
with large differences between adjacent grid cells.
The Laplacian operator for a scalar H at a cell center is discretized as
del^2(H) = d^2(H)/dx^2 + d^2(H)/dy^2,
where d^2(H)/dx^2 = (dH/dx|e - dh/dx|w) / dx, where 'e' is the east face,
'w' is the west face, dx is the cell width, and dH/dx is discretized in
the usual way. The y term is analogous.
The new Laplacian terms take the form del^2(H) * L0^2/tau0,
where tau0 is the usual timescale, and L0 is a new length scale,
typically of the same order of magnitude as dx and dy. Larger L0 leads to
more smoothing.
In glissade_grid_operators.F90, I added subroutines for computing the Laplacian
of both cell-centered and staggered variables.
The user can now specify the length scale for each inversion in the config file.
Thee default values are 0, i.e. no smoothing. On an 8km Antarctic grid,
a value L0 = 10 km noticeably smooths the fields with only a small increase
(a few m) in the rmse for H.
I hoped that this change would make the spun-up AIS more robust when inversion
is turned off, but the improvement is modest; more work is still needed.
I also changed the logic in inversion for Cc, Cp, and E.
For these parameters, we expect the physical changes to vary more in proportion
to the logs of these quantities than the absolute values.
The inversion subroutines now convert the Cc (and other quantities) to log10(Cc),
prognose changes in log10(Cc), and convert back to Cc at the end.
This is not very different from what was done before, but the logic is now a bit cleaner.
I set a default timescale of 200 yr for each inversion.
I changed the max and min values of E to 10.0 and 0.1.
At some point, these should be config parameters.
In subroutine glissade_inversion_basal_friction, I set f_ground_weight = F by default.
(Previously, it was true.) With this change, the staggering uses the unweighted thickness
in each cell, instead of weighting by f_ground. This means the inversion might
better discourage thinning of partly grounded cells. In practice, the difference
is small, although thickness biases are reduced in a few local areas.
This commit introduces a new config parameter, inversion_toggle_frequency. The goal is to make the spin-up more robust by buffering against undesired advance and retreat when inversion is turned off. If the toggle frequency = 0. (the default value), then inversion takes place based on the values of which_ho_coulomb_c, which_ho_deltaT_ocn, etc. (CISM does inversion when these options are set to 1.) If the toggle frequency > 0., then the inversion is toggled on and off at regular intervals. For example, if frequency = 1000 and CISM is configured to do inversion, then the inversion is turned on during the following time intervals (yr): [0,1000.], (2000.,3000.], (4000.,5000.], etc. At other invervals, the inversion is turned off: (1000.,2000.], (3000.,4000.0], (5000.,6000.], etc. (Note the parentheses and brackets to indicated open and closed intervals.) I also introduced a logical variable called inversion_relax_everywhere, which is true by default but false when toggling the inversion. * If true, then wherever the usual inversion conditions are not met (e.g., in floating cells when we are inverting for C_c), we relax toward the default values. The Laplacian smoothing term is also turned on, while the terms based on (H - H_obs) and dH/dt are set to 0. * If false, then the relaxation and Laplacian terms are turned off in cells that do not meet the usual inversion conditions, and we retain values computed previously when the inversion conditions were satisfied. Typically, the previous inverted values will discourage undesired advance and retreat. I removed the parameter deltaT_ocn_relax, which was always set to zero in previous runs. I added a parameter deltaT_basin_relax, which is the basin average of deltaT_ocn. This value is computed at runtime and is uniform within each basin. For grounded cells in a given basin, deltaT_ocn is now relaxed toward the basin-average value instead of toward zero.
Several commits ago, I introduced a parameter called thck_error_exponent with a default value of 1. Setting this parameter to a value of 1 modified the size of the inversion terms based on the thickness bias, raising the quantity (H - H_obs)/H0 to an arbitrary power. Tests with values > 1 showed that results did not generally improve; instead, the overall thickness bias was greater. With this commit, an exponent of 1 is assumed in all cases.
This commit moves several inversion subroutine calls from glissade_diagnostic_variable_solve (in module glissade.F90) to a new subroutine, glissade_inversion_solve (in glissade_inversion.F90). It is not quite BFB; there are roundoff-level changes in ice speed after a few timesteps, and these changes propagate to other variables. Answers agree to several decimal places. Future commits will move other *_solve subroutines out of glissade.F90. I also removed a logical parameter called f_ground_weight that was used to control the way thickness fields are staggered to vertices near the grounding line. This staggering is now always done in a simple way without weighting.
I changed the way inactive cells are treated during inversion. Here, 'inactive' means that an certain parameter is not physically relevant for a cell. For example, the value of coulomb_c makes no difference in fully floating cells, so floating cells are inactive when inverting for coulomb_c. Similarly, the value of deltaT_ocn makes no difference in fully grounded cells, so grounded cells are inactive when inverting for deltaT_ocn. In inactive cells, during inversion, we now zero out terms associated with the thickness error and dH/dt, but we keep the terms associated with relaxation and Laplacian smoothing. The values in active cells will then transition smoothly toward the values in inactive cells. Values in inactive cells will approach the relaxation value when we are far from any active cells. The goal is to make the spin-up more robust by removing limiting abrupt transitions. I removed the config option deltaT_ocn_extrapolate. This was used to set deltaT_ocn to the basin-average value in inactive cells, but is no longer needed if we continually relax to the basin-average value in these cells. I also added a config parameter, inversion_phaseout_timescale. Setting this parameter to a value > 0 increases each of the inversion timescales (babc_timescale, deltaT_ocn_timescale, flow_enhancement_timescale) by a factor exp(time/phaseout_period). For instance, if this parameter = 5000 yr, then each timescale increases by a factor of e every 5000 yr, eventually reaching a large value such that the changes due to inversion are negligible. The goal is to improve robustness by phasing out the inversion gradually instead of suddenly.
Replaced some messy diagnostics in glissade.F90 with calls to point_diag. Still need to do this in several other modules.
This commit adds a new module, glissade_mass_balance, in a file of the same name. This module will take care of calculations related to the surface mass balance (and some aspects of the basal mass balance). For now, it contains an initialization subroutine with some operations previously done in glissade_initialise. I will add more functionality in upcoming commits, moving code from the glissade and glissade_transport modules. This is part of a larger effort to clean up the glissade_module.
This commit moves several subroutines and blocks of code from the glissade and glissade_transport modules to the new glissade_mass_balance module. This is part of the overall goal of shortening and simplifying glissade. It touches a lot of lines, but the changes are BFB.
Until now, the code has had several redundant SMB fields. For example, the following pairs represent the same field with different units (mm/yr w.e. for smb, m/yr ice or scaled CISM units for acab): * smb_anomaly, acab_anomaly * smb_ref, acab_ref * smb_gradz, acab_gradz * smb_3d, acab_3d * smb_factor, acab_factor This commit removes the fields with the 'acab_' prefix. Option enable_acab_anomaly is now enable_smb_anomaly. This means that for the non-default SMB input options (SMB_INPUT_FUNCTION_XY_GRADZ, SMB_INPUT_FUNCTION_XYZ) and for SMB anomaly forcing, the user must supply input fields with units of mm/yr w.e. (smb_input = SMB_INPUT_MMYR_WE). This is already the case for most recent simulations. Note that the overwrite_acab option is unchanged and still uses acab units. In the glissade_mass_balance module, subroutine glissade_apply_smb downscales the SMB to the ice surface and adds anomalies as needed, then converts the SMB to acab units. If acab is read directly from the input file (smb_input = SMB_INPUT_MYR_ICE), no unit conversion is necessary. If acab is passed in from CESM via Glad, no conversion is necessary in Glissade, since Glad does the required conversions. These changes are BFB in the glacier and ice sheet tests I ran, but could yield roundoff-level differences when using some of these options.
This is a first step toward removing the confusing scale factors from CISM. For a long time, the variables acab and bmlt have had a scale factor called scale_acab = scyr * thk0 / tim0, used to convert acab and similar fields from unit of m/yr (in I/O) to dimensionless model units originally defined for the Glimmer code. I redefined scale_acab = scyr. This means that the conversion is from m/yr to m/s, without the extra tim0 and thk0 factors. This is consistent with the goal of having SI units for internal calculations, while using human-friendly time scales for I/O.. I removed the factor (tim0/thk0) from various modules where it appeared in conversions. In a few standard tests, results are either BFB or are changed at the machine roundoff level.
Following up on the previous commit, this seems a good time to remove the confusing scale factors from CISM. This commit removes tim0, which has units of time. As a result, the timestep dt now has units of seconds instead of being dimensionless. The old value was equal to the number of seconds in 400 yr. The method is as follows: * Comment out every line containing tim0, including the declaration in glimmer_paramets. * Replace each tim0 operation with the equivalent, assuming tim0 = 1. * Recompile the code to make sure all instances of tim0 were removed. * Run some tests to make sure the answers are the same within roundoff. In general, answers will not be BFB, since some multiplication and divisions will be modified or removed. This commit fixes a diagnostic mass conservation issue with the previous commit, which changed the units of bmlt; have to be careful adding bmlt_applied*dt to calving_thck when bmlt_applied lacks the thk0 scale but calving_thck still includes it. Later I'll delete the commented-out lines, but I'm leaving them for now to assist with debugging when we do more extensive testing. For this commit, the dome test and a glacier test are BFB, while an AIS test has answer changes within roundoff.
This commit removes len0 from the code, effectively setting it equal to 1 (instead of 2.e5 m, its previous value). The grid quantities model%numerics%dew and model%numerics%dns now have units of m. Answers change at the machine roundoff level.
In glissade_inversion_init, there is logic to remove cells with H < thklim at the start of the run. This is problematic if thklim = 1.0 m and if there are many cells in the input file with H = 1.0 m. We don't want to remove them, but they are removed if H is very slightly less than 1.0 m due to rounding. This commit fixes the issue by adding buffers against roundoff-level differences. In glissade_get_masks, the logic is that ice_mask = 1 if thck > thklim. It makes sense to have '>' rather than '>=' in case thklim = 0. But suppose thklim = 1.0 m, and many cells in the input file also have thck = 1.0 m. We would like to make these input cells as ice-covered. This will be the case if we reduce thklim by a very small amount, e.g. 1.0e-11. We now do this in glissade_initialise: model%numerics%thklim = max(model%numerics%thklim - eps11, 0.0d0) In glide_write_diag, there is a threshold parameter, diag_minthck, which is set to either 0 or thklim. The above adjustment to thklim ensures that cells with H = 1.0 m (exactly) are included in the diagnostics when the user sets thklim = 1.0 m in the config file. These fixes can lead to small answer changes.
This commit removes thk0 from the code, effectively setting it equal to 1 (instead of 2000 m, its previous value). Many basic quantities like model%geometry%thck now have units of m. This is a long-overdue code simplification. Compared to the previous commit (which resolves some rounding issues), this commit is answer-changing only at the machine roundoff level. The global sums of BMB and calving fluxes are different, but they compensate in a way that preserves global mass conservation.
This commit removes vel0 from the code, effectively setting it equal to 1 (instead of 500 m/yr * (1/scyr), its previous value). Velocities now have units of m/yr in I/O and m/s in much of the code. Thus, the scale factor for velocities in glide_vars.def is scyr. However, glissade_velo_higher still uses m/yr. I plan to switch to m/s for consistency with the rest of the code. Answers change at the machine roundoff level.
This commit removes acc0 from the code, effectively setting it to 1. Answers change at the machine roundoff level, at most.
This commit removes tau0, vis0, and evs0 from the code, effectively setting them to 1. I also removed a Glide scaling parameter called trc0. Answers change at the machine roundoff level.
This commit removes scale_acab, scale_uvel, scale_beta, and other scales previously defined in glimmer_scales.F90. These scales were formed from combinations of the scale factors thk0, tim0, len0, etc., along with scyr. Now that the scale factors have been removed, the only quantity left in scale_acab, etc. is scyr, which converts between seconds (the mks unit) and years (an easier unit to visualize for most glacier processes). Many of these scales were scale factors for the variables in glide_vars.def. Now that thk0, len0, vel0, etc. have been removed, I replaced these scales in glide_vars.def with factors of either scyr or 1.0/scyr. For instance, acab has a factor of scyr, since it has units of m/yr for output but m/s within the code. Note: For ease of debugging, I commented out most of the lines containing scale factors instead of removing them. The next step is to do some more testing and fix anything I've broken, before deleting the comments and tidying up. Answer changes are within machine roundoff.
For the past several commits, I have commented out lines containing references to thk0, tim0, scale_acab, and other scale factors, instead of removing these lines from the code. The reasoning was that if we need to debug answer differences later, it could be helpful to see the commented-out lines. This commit cleans up the code by removing the commented-out lines. It is the last step in scale factor removal. It touches many modules but is BFB.
This commit fixes a logic error that led to inexact restart for glacier runs. The problem was with the parallel_get_var functions used to read in 1D glacier fields such as glacier_mu_star and glacier_alpha_snow. (These fields are vectors with one element per glacier.) Previously, I modified several parallel_get_var functions to fetch these fields, but neglected to include a 'start' variable in the argument list. On restart, CISM was therefore reading these fields from the first time slice of the restart file instead of the most recent slice. The error was not obvious in the output, since the erroneous values are typically used for only one timestep before being recomputed, if inverting for mu_star and alpha_snow. The errors are large and persistent only if transitioning from inversion to no inversion and reading a restart file with multiple time slices. For the runs submitted to GlacierMIP, the errors were small, since the runs without inversion were read from restart files containing only the latest time slice. I found the error when switching from inversion to no inversion with a restart file containing multiple time slices. The fix is simply to call the parallel_get_var functions with a 'start' argument that identifies the correct time slice. I modified generate_ncvars.py so that glide_io.F90 and other autogenerated files are created with this argument added. I then re-ran the 10,000-year Alps spin-ups with the new code. Compared to the runs done last year for GlacierMIP3, the answer changes are small. For example: Total ice area = 2500.94 km^2 (old), 2501.10 km^2 (new) Total ice volume = 138.020 km^3 (old), 138.018 km^3 (new) Aletsch area = 88.79 km^2 (old), 88.80 km^2 (new) This commit should come after the commit 791597e on the main CISM branch.
Until now, the DIVA velocities have been too large in regions with steep slopes.
This can be understood as follows: The equations are derived with the assumption
that the flow is roughly parallel to the x-axis, which is not the case for steep slopes.
The d/dx(mu du/dx) term in the momentum balance should include a component
perpendicular to the bed, with a proportionality constant of sin^2(theta),
where theta is the angle between the bed and the x-axis.
In the standard DIVA equations, this component is missing, since the x derivative
of the vertically-averaged u is typically small. As a result, DIVA underestimates
the resistive interior stresses and overestimates the ice speed.
A solution is suggested by an analysis of the slab problem presented in Dukowicz (2012, TC).
For this problem, John derived analytic solutions for the velocity based on the Stokes
and BP equations. Studying these solutions, I noticed that the DIVA solution can be
modified to reproduce the analytic BP solution (at least for the slab problem)
by making the following changes:
(1) When assembling the basal BC terms in the stiffness matrix, integrate over the
actual area of the basal surface (in general, a curved 2D surface embedded in 3D)
rather than the projection of this surface onto the xy plane.
This means dividing the projected area by cos(theta), which increases the resistance
to sliding.
(2) Compute a slope correction term gamma = cos^2(theta) + 4 sin^2(theta).
Insert this term in the following places:
(a) Instead of setting beta_eff = beta/(1 + beta*F2),
set beta_eff = (gamma*beta) / (gamma + beta*F2*cos(theta)), thereby
increasing beta_eff.
Recall that beta_eff is defined such that beta_eff * u_mean = beta*u_b,
so a larger beta_eff implies a smaller u_mean.
(b) When computing the 3D DIVA velocity profile, multiply the vertical shear
by cos(theta)/gamma, thus reducing the shear.
(c) When computing du/dz, set du/dz = (tau_x*sigma)/(efvs*gamma), instead of just
(tau_x*sigma)/efvs. Here, sigma = (s - z)/H is the unitless vertical coordinate.
(d) When computing the effective strain rate, multiply the (du/dz)^2 and (dv/dz)^2 terms
by gamma, thus increasing the effective viscosity for Glen's n > 1.
With these changes, the DIVA speeds match the BP speeds almost exactly for the slab problem.
In general, these corrections work well for slab-like problems where the flow is dominated
by shear perpendicular to the bed. They do not work as well for cases with significant
along-flow tensile or lateral stresses (e.g., ISMIP-HOM test A).
The changes are needed in both the x and y momentum equations. So we compute slope correction terms
for both axes, leading to separate parameters beta_eff_x and beta_eff_y.
To implement the corrections, the user should set diva_slope_correction = .true. in the config file.
The default (at least for now) is .false., which reverts to the standard DIVA solver without slope corrections.
At some point we'll change the default to true.
I used these corrections for the runs in the CISM glacier paper (Minallah, Lipscomb et al., 2025, GMD).
This made it possible to increase the max_slope parameter to 1.0 (i.e., a 45-degree angle)
without getting excessively large flow speeds.
Minor changes and fixes:
* A fix to support periodic BCs in the north-south direction for the slab
* Don't allow Cp to be reset at restart when running with Cp inversion for glaciers
* Option to count the fraction of edges steeper than max_slope
I added a few features to the slab test case to test the DIVA slope correction option: * The user can now specify the angle phi between the direction of steepest descent and the x axis. Previously we had phi = 0 implicitly; now the user can set phi = 0, 90, 180, or 270 degrees. The periodic BCs do not allow diagonal directions. * I added a cos(theta) correction to the ice thickness so that the analytic solution is correct for steep slopes. * I added the option diva_slope_correction = .true. to the default config file. With these changes, I showed that with the slope correction, the CISM DIVA solution for the slab problem is nearly indistinguishable from the BP solution.
The first line in runISMIP_HOM.py contained an erroneous 'python2'. Changed to 'python' so it will run cleanly on newer python versions.
This commit consolidates two chunks of code which carry out operations (downscaling, adding anomalies, etc.) related to climate forcing: (1) code in glissade_thermal_solve to compute artm, precip, and/or snow (2) code in glissade_apply_smb to compute acab These calculations are now done in a new subroutine, glissade_prepare_climate_forcing, in module glissade_mass_balance.F90. The subroutine is called early in the timestep, before glissade_thermal_solve, since the thermal solver needs artm. The acab field is computed here and used later in the timestep, in glissade_apply_smb. I also removed some obsolete variables like thck_unscaled and topg_unscaled. Now that model%geometry%thck and topg have units of m, the alternate variables are not needed. This commit is BFB.
This commit includes several changes related to glacier spin-up, anomaly forcing, and diagnostics: * I modified subroutine glissade_invert_powerlaw_c to more closely resemble subroutine invert_basal_friction in module glissade_inversion, as recently modified for Antarctic spin-up. The glacier calculation now applies the prognostic terms to log(Cp). In ice-free regions, Cp relaxes toward its default value. This is answer-changing. * In glide_diagnostics, I replaced a couple of manual rms error calculations with calls to subroutine glissade_rms_error. The diagnostic results are unchanged. * I modified the anomaly calculations so that, by default, the entire anomaly for a given year is applied starting from Jan. 1 of that year. This is different from the initMIP convention of delaying until the following year. This leads to answer changes for the historical and commitment runs. I left code comments on how to implement other conventions. * Turned off some verbose diagnostics (but left verbose_glacier = .true. for now)
This commit simply changes the default setting to .true. for the new DIVA slope correction option. This means that older config files do not need to set the option explicitly. We can reproduce older results by setting the config option to .false.
When implementing the DIVA slope correction for glaciers, I changed the DIVA 3D velocity calculation in a way that leads to NaNs for marine-based ice. This commit fixes the problem. It leads to small answer changes for land-based ice, compared to the initial DIVA slope correction. With the correction turned off, the agreement with earlier code (before the slope correction) is now at roundoff level.
The goal is to give the ais_spinup branch the same history as the current version of main, up to the first ais_spinup commit. This should make it easier to bring the ais_spinup branch to main. Merge branch 'main' into lipscomb/ais_spinup.tmp
Preparing the smb_cleanup branch to be merged to main. Actual history was (1) Start the ais_spinup branch (2) From this branch, create the smb_cleanup branch (3) Continue the ais_spinup branch (mostly glacier commits) (4) Copy ais_spinup to ais_spinup.v2, and merge main (5) Copy smb_cleanup to smb_cleanup.v2, and merge ais_spinup.v2
When I changed the DIVA 3D velocity calculation to be correct for large slopes, I introduced a divzero for floating ice. Then when I fixed the divzero, I broke the calculation for large slopes. With this commit, I think the velocity is correct for both floating ice and large slopes.
When testing methods of making the inverted ice state more robust, I added options to (1) toggle the inversion on and off and (2) gradually phase out the inversion by lengthening the inversion timescales. Neither method seemed to improve robustness, so I am removing them.
This commit fixes a rare error that can occur when doing remapping transport for very thin ice layers. The error happens when, within machine precision, the total column thickness is greater than zero, but the thickness of individual layers is not. The fix is to call make_remap_mask for each layer with thck_layer as an input argument, instead of calling it just once with thck as an argument. Results are BFB except in the rare cases when errors used to be triggered. Thanks to Tim van den Akker for diagnosing the problem and figuring out the fix.
Some time ago, I commented out some logic setting f_ground = 1 for all land-based cells (including ice-free cells), but I can't remember why. It seems sensible to set f_ground = 1 on land, so I uncommented this code. This leads to small answer changes. Also did some minor code cleanup in glissade_inversion.
This commit replaces verbose print statements with calls to subroutine print_diag in several modules. Answers are BFB.
Like the previous commit, this commit replaces verbose print statements with calls to subroutine point_diag in several modules. Answers are BFB.
This commit replaces 'print*,' with 'write(6,*)' throughout the code. Later, might want to replace '6' with a unit variable.
This commit modifies the algorithm for extrapolating ocean thermal forcing
into sub-ice-shelf cavities. The new algorithm reduces the unphysical striping
in the previous version, giving smoother and more reasonable fields.
Thanks to Michele Petrini for testing and improving the algorithm.
The algorithm now works as follows:
* At startup, set TF to an unphysical value everywhere in sub-shelf cavities.
* Extrapolate off-shelf values into cavities, iterating until there are physical values
in all cells that levels that need them. Weigh adjacent cells using a Laplacian smoother.
- With a 1-level Laplacian smoother, the extrapolation gets stuck from time to time.
In that case we fill levels in the vertical direction as needed, and then continue the iteration.
- With a 3-level Laplacian smoother, the extrapolation fills nearly all cells and levels
without getting stuck. The only levels remaining are those which are separated by one or more
levels or cells from the nearest filled value.
- When the Laplacian smoother has gone as far as it can, fill the remaining levels
by averaging valid values at all vertical levels in adjacent grid cells.
Three stencil sizes are supported for the Laplacian smoother:
* 9 (3 x 3 at one level)
* 25 (5 x 5 at one level)
* 27 (3 x 3 at three levels)
The 27-point smoother converges in fewer iterations than the other two, so it is the default.
Each cell and level that starts with TF = unphys_val receives its initial value by taking
a smoothed average of neighboring values that have been filled.
This initial value can be smoothed one or more additional times by setting max_count_smoothing > 1.
The default is max_count_smooothing = 2, which gives one additional smoothing.
This additional smoothing reduces the striping along grid diagonals in large cavities, especially Ross.
Further smoothing does not make a big difference.
With max_count_smoothing = 2 or more, the 9-point stencil gives results similar to the 25-point stencil.
For now, max_count_smoothing is hardwired; it could be a config option if desired.
Optionally, the user can use a parameter called cavity_buffer to set TF to an unphysical value
in open-ocean cells that are near but not at the calving front (either 1 or 2 cells away).
The default buffer is 0, but it can be set to 1 or 2. It is hardwired; could be a config option if desired.
The full extrapolation, starting without any valid values in cavities, is done only at initialization.
On subsequent timesteps (and on restart), valid values in cavities are retained; the extrapolation is applied
only where the grounding line has retreated.
I compared the extrapolated TF at the lower ice surface to the values obtained from Xylar's ISMIP6 TF data set,
where TF is already extrapolated into cavities. A couple of differences:
* For Larsen C the CISM TF is notably larger; I'm not sure why.
* For Ross, the CISM TF is larger in the Byrd sector than in most of the cavity; this difference does not appear in Xylar's TF.
For most cavities, however, the two patterns look pretty similar.
This commit removes the use of, and calls to, subroutine point_diag from modules glissade_grid_operators and glissade_utils, where using the glide_diagnostics module creates circular dependencies. In the longer term, we might want to move point_diag to glissade_utils or some other low-level module.
An earlier commit added a Laplacian smoothing term to the equations for inversion of coulomb_c, powerlaw_c, deltaT_ocn, and flow_enhancement factor. The goal was to make the inversion more robust by having smoother transitions in ice properties near grounding lines. I found that in Antarctic spin-ups, applying the Laplacian term to coulomb_c can cause oscillations. Instead of converging in 25 to 30 iterations, the nonlinear solver fails to converge. This significantly slows down the code, from < 1 hr per 1000 model years to about 1:45 per 1000 model years on 256 Derecho cores. This commit sets the Laplacian term to zero in the C_c and C_p inversions for cells that are at least partly grounded. The Laplacian term is applied only in floating cells and ice-free ocean, where it is balanced by a relaxation term. Thus, we still have a fairly smooth transition in C_c or C_p when the GL advances downstream, but downstream floating cells no longer have a smoothing influence on upstream grounded cells. With this change, performance improves. The time to run 1000 yr is again < 1 hr. Applying the Laplacian term to inversion of deltaT_ocn and flow_enhancement_factor does not significantly slow down the code, so I left the Laplacian term in those equations (at least for now). I tried a 10 kyr spin-up with inversion, followed by 5 kyr without inversion. The spin-up is fairly robust. There is some localized thickening and thinning (e.g., some thickening and GL advance in Queen Maud Land around the Lazarev shelf), but overall, not much change in grounding lines in key regions such as the ASE, Filchner-Ronne, and Ross. I also fixed a minor error in the verbose diagnostics for Picard acceleration. This commit is answer-changing only for runs with Laplacian smoothing during inversion.
Until now, inversion for coulomb_c has been local; the value at each vertex is adjusted based on the thickness at that vertex. This commit introduces an alternative, which_ho_coulomb_c_basin: * The default is 0, implying no basin-scale inversion of coulomb_c. * When which_ho_coulomb_c_basin = 1, CISM inverts for coulomb_c at the basin scale based on the basin-mean thickness of grounded ice. * When which_ho_coulomb_c_basin = 2, CISM reads the basin-scale values from an external file (typically, these would be values from a previous inversion). The logic is similar to that for the which_ho_deltaT_basin option, but with tendency terms based on those for local coulomb_c inversion. The thickness targets are given by the initial thickness of grounded ice, interpolated to vertices. They are contained in a new field, grounded_thck_target, that can be written to and read from restart files. For summing and averaging vertex quantities over basins, I adopted the convention that vertex (i,j) belongs to basin nb if and only if cell (i,j) belongs to basin nb. In short test runs, the new inversion option seems to be working as desired. I verified exact restart.
…ions
This commit extends the basal_hydro derived type with options and parameters related to
subglacial hydrology. Most of these can be set by the user in the config file.
Some variables that were previously part of the basal_physics type are now in the basal_hydro type.
in order to separate hydrology variables from variables connected to basal sliding laws.
In glide_setup, there are two new subroutines, handle_basal_hydro and print_basal_hydro,
to deal with this derived type. In the future there will be similar changes for other sets
of closely related variables and parameters.
The which_ho_effecpress options have been modified. The options are now
HO_EFFECPRESS_OVERBURDEN = 0
HO_EFFECPRESS_BPMP = 1
HO_EFFECPRESS_BWAT = 2
HO_EFFECPRESS_CAVITY_SHEET = 3
HO_EFFECPRESS_BWAT_BVP = 4
* Option 0 (N = overburden = p_i = rhoi*g*H) is unchanged and is still the default.
* Option 1 is unchanged.
* Option 2 has been generalized; it now supports a macroporous-sheet option.
The term 'macroporous sheet' is taken from de Fleurian et al. (2018), the SHMIP paper.
The water pressure is given by
p_w = pi * min((bwat/bwat_threshold)**gamma, 1.0),
where bwat is the diagnosed basal water and gamma is an exponent.
The default values are bwat_threshold = 0.1 m and gamma = 7/2 (as in Flowers & Clarke, 2002).
With this scheme, p_w = p_i (i.e., N = 0) for bwat >= bwat_threshold.
To set a nonzero min value for N, the user should choose a nonzero value for effecpress_delta.
The default (as before) is effecpress_delta = 0.02, but this value can be modified in the config file.
* Option 3 is a new cavity-sheet parameterization, also described by de Fleurian et al. (2018).
With this scheme, N is found by setting the cavity opening rate equal to the closing rate.
There are two modes of opening: by sliding and by melting.
- Opening by sliding is given by v_o = u_b * max(bump_height - bwat, 0) / bump_wavelength,
where u_b is the sliding speed (either the modeled speed or a fixed parameter),
and bump_height and bump_wavelength are empirical constants.
The new option cavity_open_slide allows the user to choose either fixed or dynamic u_b,
or to keep this term turned off.
- Opening by melting is given by V_o = bmlt_ground + (qflx*grad(phi))/(rhoi*L),
where bmlt_ground is the basal melt rate based on geothermal, frictional, and conductive fluxes,
and qflx*grad(phi)) is the rate of energy production by dissipation at cavity walls.
The new option cavity_open_slide allows the user to choose either of these two opening terms,
or both, or neither.
Closing is given by v_c = c_close * A * bwat * N^n,
where n = Glen's n, A is a flow factor, and c_close is an empirical constant, set to 2/n^n by default.
With v_c = v_o + V_o , it is straightforward to solve for N.
* Option 4 is unchanged; it is used with the local basal till model.
I removed the old HO_EFFECPRESS_BWATFLX option.
CISM now distinguishes between two ways of computing the basal water depth:
(1) bwat is the depth prognosed in the local basal till model, used with option 4 above.
When bwat > 0 at the bed, CISM requires the water to refreeze before the bed temperature T_b
can drop below the pressure melting point, T_pmp.
(2) bwat_diag is the depth diagnosed from the steady-state basal water fluxes in the flow-routing scheme.
In this scheme, basal water can be present when T_b < T_pmp (if it arrives from upstream),
so bwat_diag > 0 does not require T_b = T_pmp.
Options 2 and 3 use bwat_diag to compute N.
In the glissade_basal_water module, I changed the logic for initializing phi_in
subroutine fill_depressions, to get clean results at boundaries.
I also did some cleanup in this module, removing deprecated subroutines
(including fix_flats and find_flats) and adding some point_diag calls.
I tested the new macroporous-sheet and cavity-sheet options with a simple dome geometry.
The answers look reasonable, at least qualitatively.
I checked that for an Antarctic run with N based on ocean_p but otherwise set to overburden,
the results are BFB as expected.
Still to do:
* Implement a SHAKTI-style transition from laminar to turbulent flow.
* Test the new options for a SHMIP test case.
Previously, there were separate options to invert for Cc locally (which_ho_coulomb_c = 1) and at basin scale (which_ho_coulomb_c_basin = 1). Similarly, there were separate options to invert for deltaT_ocn locally (which_ho_deltaT_ocn = 1) and at basin scale (which_ho_deltaT_basin = 1). This commit removes which_ho_coulomb_c_basin and which_ho_deltaT_basin as separate options. Basin-scale inversion is now folded into the options which_ho_coulomb_c and which_ho_deltaT_ocn. The user should set these options to 3 to do inversion at basin scale. Answers for basin-scale inversion are not changed; only the option numbers have changed. In removing which_ho_deltaT_basin as an variable, I removed the old option HO_DELTAT_BASIN_ISMIP6 = 3, which prescribed deltaT_ocn in each basin based on Nico Jourdain's values for ISMIP6. This option has not been used recently. Also, I added an option for basin-scale powerlaw_c inversion: which_ho_powerlaw_c = 3. This mirrors the recently added option for basin-scale coulomb_c inversion.
This commit makes two changes related to the effective pressure N:
(1) There are two distinct ways to compute N: based on subglacial hydrology
(e.g., N is a function of bwat) or based on the flotation thickness and ocean_p.
Previously, these were combined in such a way that if both had the effect
of reducing N, the combined effect was greater than the effect of either on its own.
With this commit, the two calculations are independent, and CISM takes the minimum
of the two values.
(2) The flux-routing hydrology scheme computes the basal water depth bwat.
Then for the option HO_EFFECPRESS_BWAT, it computes N as a function of bwat.
The standard macroporous-sheet formulation increases the water pressure in proportion
to (bwat/bwat_threshold)^gamma, where gamma >=1. I found that this formulation leads
to instability, since modest increases in bwat can cause a large, sudden reduction in N.
I implemented a more stable alternative, where the water pressure increases
in proportion to bwat / (bwat + bwat_threshold).
This commit moves subroutine point_diag from the glide_diagnostics module to the glimmer_utils module to avoid circular 'use' statements. The issue is that glide_diagnostics is not very low-level; it uses some glissade modules from which it could be convenient to call point_diag. In that case, point_diag needs to live elsewhere.
This commit adds an option smb_input_function = SMB_INPUT_FUNCTION_PDD = 3, which computes the SMB using a simple positive-degree scheme. The scheme works as follows: (1) Downscale artm_ref to the surface to get artm, if needed (based on artm_input_function). (2) Based on artm, partition the precip between rain and snow. (3) Based on artm, compute the ablation. (4) Computes SMB as snowfall minus ablation. Step (2) requires snow threshold temperatures, which I added as config variables. Step (3) requires a degree factor and a melt threshold, also added as config variables. We now have duplicates of several variables (tmlt, snow_threshold_min, snow_threshold_max): one variable in the climate derived type and another in the glacier derived type. In a future commit, these should be consolidated. The new option will be used for AIS paleo spin-ups where artm_ref is computed by interpolating between LIG, LGM, and PI climates based on a climate index.
This commit changes some coulomb_c options as follows: * changes the basin-scale inversion options for which_ho_coulomb_c * removes the option which_ho_coulomb_c_relax * adds the option elevation_based_coulomb_c. The new set of which_ho_coulomb_c options is HO_COULOMB_C_CONSTANT = 0 HO_COULOMB_C_INVERSION = 1 HO_COULOMB_C_EXTERNAL = 2 HO_COULOMB_C_INVERSION_BASIN = 3 HO_COULOMB_C_EXTERNAL_BASIN = 4 Options 0, 1 and 2 are as before. Option 3 is modified: * In the old version, the inversion was for a single value of Cc in each basin. A single value turns out to be unrealistic. This value is too low in high-elevation regions and/or too low in low-elevation regions. * In the new version, Cc is assumed to be elevation-dependent: Cc = coulomb_c_hi for z >= bed_hi Cc = coulomb_c_lo for z <= bed_lo log(Cc) is interpolated linearly between log(coulomb_c_lo) and log(coulomb_c_hi) for z between bed_lo and bed_hi This is the same elevation-based formulation CISM has used before, except that log(Cc) instead of Cc is now linearly interpolated. Given this elevation dependence, the inversion is for coulomb_c_lo in each basin. This makes it possible to invert for a low value in some basins (e.g., the Siple coast) and higher values in other basins. Unless Cc_lo is small, the Siple ice tends to be too thick. In the new code, both coulomb_c_lo and coulomb_c_hi are 2D fields, although only coulomb_c_lo is spatially varying; coulomb_c_hi is set to a constant everywhere. The elevation target is the mean elevation of marine-grounded_ice at initialization. Option 4 is new. It is related to option 3 in the same way option 2 is related to option 1. If doing basin-scale inversion, we use option 3 for the spin-up, and then run forward with option 4, reading the saved value of coulomb_c_lo from the restart file. I removed option which_ho_coulomb_c_relax, which has not been very useful. When inverting locally for coulomb_c, the relaxation is now toward coulomb_c_const rather than an elevation-based range of options. I added the logical config option elevation_based_coulomb_c. The default is F. This option must be set to F for which_ho_coulomb_c = 1 or 2, and to T for which_ho_coulomb_c = 3 or 4. To run with fixed elevation-based Cc, the user sets which_ho_coulomb_c = 0 and elevation_based_coulomb_c = T. This results in both coulomb_c_hi and coulomb_c_lo being constants, which are 0.5 and 0.1 by default but can be set in the config file. The defaults for bed_hi and bed_lo are 0 m and -500 m. The default bed_hi is lower than I've used before, to allow higher Cc for all ice grounded above sea level. These values can be set by the user in the config file. I tried some AIS spin-ups with the new which_ho_coulomb_c = 3, along with active subglacial hydrology. Results look reasonable overall, but there are some large thickness biases, e.g.: * PIG is too thick and Thwaites is too thin * The Siple coast ice is too thick near the GL and too thin farther upstream. * The Transantarctic Mountains and AP are too thick. The hope is that Thwaites will improve with a new target thickness dataset, which would improve PIG as well (since the Thwaites errors leak into the PIG basin). The Siple coast might improve with more work on the hydrology. The TA and AP might improve with either hydrology changes or finer resolution. Other changes: * Replaced the parallel_is_nonzero function with parallel_is_zero, and added several calls to this function to simplify some logic * Fixed the SMB units for one of the calls to check_fill_values (to check dthck_dt_obs) * Reduced the default for bwat_threshold to 1.e-3 I confirmed that this commit is BFB for several Antarctic and glacier applications that do not use the new options.
A key link between the active subglacial hydrology and the ice dynamics is the function that relates h = bwat (the subglacial water depth) to the effective pressure N. The previous function was N/overburden = 1 - h/(h + h0), where h0 = bwat_threshold. This function asymptotes to zero very slowly as h increases relative to h0. The new function is N/overburden = max[delta, exp(-h/h0)], where h0 = bwat_threshold is now an e-folding scale, and delta ~ 0.02. This function asymptotes to delta quickly as h increases relative to h0, but not so quickly as to make the dynamics unstable, at least in tests to date.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 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 began with some cleanup and changes to CISM's surface mass balance options, and then absorbed some separate work going on at the same time. Here's a summary of changes:
Created a new module, glissade_mass_balance, consolidating some code that was spread among several modules. The new module contains a subroutine to prepare SMB and artm forcing in the form CISM needs (based on the input climate fields and the choices for artm_input_function and smb_input_function), as well as subroutines to apply the SMB at runtime and then remap tracers to the prescribed sigma levels. There is a new option to compute the SMB using a PDD scheme.
Removed the scaling factors (thk0, tim0, len0, vel0, etc.) that had been in the code since Glimmer days. Most fields now use MKS units. The SMB-related fields use units of either (1) mm/yr water equivalent or (2) m/s of ice. Generally, any I/O field whose name includes 'smb' is in the first category, while the 'acab' and 'bmlt' fields used for internal calculations are in the second category.
Fixed a restart error for glacier runs. This error led to restarts that were close but not exact. The final runs for the glacier framework paper (Minallah et al., 2025, GMD) incorporate the fix.
Added a slope correction for the DIVA solver. I found a way to improve the accuracy of the DIVA equations for large slopes by adding some functions of theta (the angle formed by the ice surface with the horizontal) to the equations. The final runs for the glacier framework paper include these correction terms.
Modified the extrapolation of ocean thermal forcing to sub-shelf cavities. The new algorithm incorporates some ideas of Michele Petrini. The biggest change is that the extrapolation is now done using a 3-layer, 27-point Laplacian stencil, instead of alternating frequently between horizontal and vertical extrapolations. This greatly reduces unrealistic striping in the extrapolated thermal forcing.
Added optional Laplacian smoothing terms to the inversions for coulomb_c, powerlaw_c and deltaT_ocn. These terms make the inverted fields less noisy.
Added an option to invert for coulomb_c at the basin scale in spin-ups. This greatly reduces the number of free parameters, compared to inverting for coulomb_c at every vertex.
Modified the which_ho_effecpress option to connect the effective pressure N to the subglacial hydrology. The new option 2 relates N to basal water depth, and option 3 computes N based on simple equations for cavity opening and closing.
Replaced many diagnostic print statements with calls to subroutine point_diag, which now lives in the glimmer_utils module. Replaced many other print statements with write statements.