Lipscomb/calvingmip.2026.merge#84
Conversation
This commit supports the final CISM submission for CalvingMIP. Diagnostic additions: * CISM now computes the radial distance from the origin, the effective thickness, and the ice speed at the calving front for each of 2 axes for Experiments 2 and 4. For the circular domain, the axes are the y-axis and the line y = x (in the NE quadrant). For the Thule domain, the axes are Caprona A and Halbrane A. These 4 axes are the ones plotted in the CalvingMIP paper by Jordan et al. (2026). I added the scalars cf_radius1, cf_radius2, cf_thck1, cf_thck2, cf_speed1 and cf_speed2 to glide_types and glide_vars.def. This means that users will be able to compare results directly to the published values from CISM and other models without doing any post-processing. Still to do is to extend these diagnostics to all 8 axes per experiment. Algorithmic changes: * In full CF cells, H_eff(i,j) is always computed based on the interior neighbor thickness instead of capping H_eff(i,j) at H(i,j), so we can have H > H_eff temporarily. * If a CF cell has no interior edge neighbor, it now computes H_eff from an interior corner neighbor. * For calvingMIP experiments, CISM now calls subroutine advance_calving_front after applying the calving_based thickness change. This ensures H <= H_eff at the CF when the algoritm is finished. I modified this subroutine so that ice advances only into edge neighbors and not corner neighbors. This prevents the creation of thin interior cells when the CF advances. * I fixed a bug in the computation of fluxes from corner neighbors in glissade_input_fluxes. * The protected CF cells (i.e., cells that are allowed to accumulate ice rather than have that ice shifted back upstream) now include cells with two partial CF neighbors (which typically have a full interior diagonal neighbor). This allows ice to accumulate in cells at CF corners. * I modified remove_icebergs so that any ice-covered cell, including partial CF cells, can spread the fill to edge neighbors. With this change, we don't need to call the glissade_fill_with_buffer subroutine (which was written to account for partial CF cells that were dynamically inactive), but just the standard glissade_fill. * In glissade_utils, I added subroutines that (1) compute global sums over the four calvingMIP quadrants and (2) compute the value of a field at a point in a bounding box, given the values at the corners. With these changes, the CF advance/retreat rates are similar to before (as desired), and oscillations in H_eff and speed at the CF are significantly reduced, especially when the CF is advancing.
The 25-point Laplacian stencil contained a '26' instead of '36' in the mask computation. Fixed the typo.
This change makes CISM's no-advance calving mask more flexible. Previously, this mask was an integer, either 0 or 1. Now the mask can be either (1) an integer array (still called calving_mask) if which_ho_calving_front = 0, or (2) a real array (called subgrid_calving_mask) if which_ho_calving_front = 1. When using the real array, CISM will thin ice in partial CF cells so that areafrac = thck/thck_effective = (1 - subgrid_calving_mask). The subgrid calving mask can be used when which_ho_calving_front = 1 (i.e., the subgrid CF scheme). The mask is computed at initialization. It can be based on an optional config parameter called calving_front_radius. If calving_front_radius /= 0, subgrid_calving_mask is computed to enforce this radius. This option is now used for CalvingMIP spinups. By default, the mask is computed based on the initial ice thickness, with a value of 0.0 or 1.0 everywhere. I applied the new subgrid mask to the two calvingMIP spinups (circular and Thule). Previously, the final CF radius was very close to 750 km along the x- and y-axes, but around 752 km along the diagonals because of the binary mask. Now the CF radius along diagonals is very close to 750 km (within < 0.1 km). I combined three mask operations (apply binary mask, apply subgrid mask, and apply forced retreat) in one subroutine, glissade_apply_calving_mask, which is now called from glissade_calving_solve after doing the main calving. Also, I improved the CalvingMIP diagnostics. Previously, the CF radius, speed, and thickness were computed only for axes 1 and 2 of each experiment. Now they are computed for all 8 axes. There are output fields called cf_radius, cf_speed, and cf_thck, each with dimension 8. I created a new dimension, naxis = 8, to support this output. I made the logic more consistent so that in general, the CF location is found by interpolating between two points, one with a_eff > 0.5 and the other with a_eff < 0.5. I put some of the complicated Caprona logic in subroutines to allow code reuse. Other changes: * The value of thck_effective now cannot exceed the flotation thickness. * I reverted an earlier change: the iceberg removal subroutine again calls glissade_fill_with_buffer instead of glissade_fill. Using the buffer allows ice to accumulate, as desired, in ice-free ocean cells at the CF. * For calving masks, thck_effective is not set to 0.0 in all cells with ice_mask = 0, and thck_effective is not allowed to exceed the flotation thickness. * I wrote a subroutine, glissade_cleanup_tiny_thickness, which is now called after the calving, at the end of the prognostic part of the timestep. This subroutine removes ice from cells with H < 1.e-11 m, adding it to the calving flux. This avoids starting the next timestep with tiny but nonzero amounts of ice in cells near the calving front. This subroutine, along with glissade_cleanup_icefree_cells, is in the glissade_utils module. * The 'thklim' argument in some calls to glissade_get_masks now depends on the CF numerics. With a subgrid CF, the threshold is 1.e-11 (so partial cells with H < 1 m can have ice_mask = 1). Without a subgrid CF, the threshold is model%numerics%thklim (typically 1 m). * I reduced the small_dthck parameter in subroutine advance_calving_front to 0.1 m, to avoid a too-steep gradient at the CF. * In subroutine add_surface_and_basal_mass_balance, most computations are now skipped in cells with effective_areafrac = 0. * Quadrant sums are now reproducible for any processor count. With these changes, all the CalvingMIP spinups and experiments are performing well.
This commit contains changes in the calvingMIP test directory, specifcally the files calvingMIP.config.template, calvingMIP.Setup.py, and README.calvingMIP. These are the files used to set up the final CISM runs for the CalvingMIP paper. Among the changes: * Both spin-ups now use marine_margin = 6, calving ice at the margin based on subgrid_calving_mask. The calving amplitude is set to a large positive value to allow free advance until reaching the margin. * The spin-ups and Experiment 4 have apply_calving_mask = True. These runs need a mask to limit CF advance when no other calving is taking place. * All runs now have which_ho_flotation_function = 2 (which is the same as the old 3). * Experiments 1 to 4 write cf_radius, cf_thck and cf_speed to the scalar output for each of 8 axes for each experiment. * In the input files, topg is computed in a way that enforces symmetry. * The input files no longer contain calving_mask. Rather, subgrid_calving_mask is computed at startup. * The prescribed value of dthck_dx_cf is now 1.e-4 instead of 5.e-4. The smaller value gives slightly more accurate CF retreat.
This commit changes the Caprona CalvingMIP diagnostics to be more symmetric. I replaced bounding_box calculations with more accurate interpolations along the Caprona axes, and removed subroutine glissade_bounding_box. I added four new diagnostic fields: * cf_locx and cf_locy, the x and y coordinates of the calving front along each axis * cf_uvel and vf_vvel, the ice speed at the CF along each axis I removed cf_speed, which is redundant given cf_uvel and cf_vvel. Also, I fixed a bug in the computation of thck_effective just before the velocity solve. I was computing thck_effective = 50 m (= thck_effective_min) at the calving front, instead of the correct value. (Elsewhere in the code, thck_effective was computed correctly.) This fix causes modest changes in ice velocity at the CF.
Fixed two conflicts along the way: one in glide_nc_custom.F90 and one in glissade.F90
This commit removes a couple of unnecessary deallocation statements that I overlooked when merging in the latest version of main.
|
@Katetc, here is the calvingmip.2026 PR. This morning I merged in the latest changes from main (i.e., the recent reprosum changes), renaming the branch to lipscomb/calvingmip.2026.merge. It's ready for testing. One thing to consider is whether we want to add any CalvingMIP tests to the test suite. It would be handy to have these tests to make sure future code changes don't break CalvingMIP, but I'm not sure how hard it would be to add standalone CISM tests on idealized grids. The README file and setup scripts are in .../tests/calvingMIP. |
Katetc
left a comment
Katetc
left a comment
There was a problem hiding this comment.
Pretty low-key comments here. I'm not going to wait for these to start testing. I'll approve the PR now in case you and Gunter decide to push it through while I'm gone.
| !! (floating_mask(i,j-1) == 1 .and. thck(i,j-1) < thck_effective_min) .or. & | ||
| !! (floating_mask(i,j+1) == 1 .and. thck(i,j+1) < thck_effective_min) ) then | ||
| !! calving_front_mask(i,j) = 1 | ||
| !! endif |
There was a problem hiding this comment.
I'm not a fan of leaving commented out code hanging around. We can recover this code through GitHub very easily, if needed.
| !! call parallel_globalindex(i, j, ig, jg, parallel) | ||
| !! write(iulog,*) 'No interior neighbor:', ig, jg, thck(i,j) | ||
| !! write(iulog,*) ' New H_eff:', thck_effective(i,j) | ||
| !TODO - Look at cases with no interior neighbors |
There was a problem hiding this comment.
Is this a real "TODO" or a left over comment? It seems vague
| ! This can be useful in idealized experiments like CalvingMIP to check for | ||
| ! violations of reflectional or rotational symmetry. | ||
| ! Note: These sums are not independent of processor count | ||
| ! TODO: Make them reproducible, using quadrant masks? |
There was a problem hiding this comment.
So this code won't work with repro_sums? I think that's fine, but maybe we should add an issue in CISM github so we can track stuff like this.
| subroutine glissade_cleanup_icefree_cells(model) | ||
|
|
||
| ! Clean up prognostic variables in ice-free cells. | ||
| ! This means seting most tracers to zero (or min(artm,0) for the case of temperature). |
There was a problem hiding this comment.
Does this impact conservation of tracers? (thinking about isotopes here)
| model%geometry%topg, & ! m | ||
| model%climate%eus, & ! m | ||
| model%numerics%thklim, & ! m | ||
| this_thklim, & ! thklim (m) |
There was a problem hiding this comment.
This unit is a little odd. I think it should still be just 'm', we can see that it's thklim in the name of the variable...
| ! are alse removed. | ||
| ! | ||
| ! Note: Option 2 is now part of subroutine glissade_apply_calving_mask. | ||
| ! Consider whether the following logic could go in the same subroutine, or if it is still needed. |
There was a problem hiding this comment.
See, this one doesn't even have a "TODO". A good idea doomed to be lost forever... unless it was put into a github issue
| public :: verbose_calving | ||
|
|
||
| logical, parameter :: verbose_calving = .false. | ||
| !! logical, parameter :: verbose_calving = .false. |
There was a problem hiding this comment.
I disagree with the need for this commented out setting.
| call parallel_halo(subgrid_calving_mask, parallel) | ||
|
|
||
| deallocate(ice_mask) | ||
| deallocate(ocean_mask) |
There was a problem hiding this comment.
Still not a bad idea to add an "if allocated" for stuff like these.
| endif | ||
| enddo | ||
|
|
||
| y_axis_thru_edges = .false. | ||
| do i = nhalo+1, nx-nhalo | ||
| if (x0(i) == 0.0d0) then | ||
| y_axis_thru_edges = .true. | ||
| exit |
There was a problem hiding this comment.
Ok, these exit calls make me a little nervous. "Exit" can also be "stop my program right here" in some compilers (though called slightly differently usually). I just think that your number of values in a nhalo+1 to nx-nhalo sized loop is not very big, so exiting a bit early is not saving you much.
Or, could we do this with an "any" call like
if( any(x0(nhalo+1:nx-nhalo) == 0.0d0)) then y_axis_thru_edges = .true.
| do i = nhalo+1, nx-nhalo | ||
| if (x1(i) == 0.0d0) then | ||
| y_axis_thru_centers = .true. | ||
| exit |
There was a problem hiding this comment.
As I commented below:
Ok, these exit calls make me a little nervous. "Exit" can also be "stop my program right here" in some compilers (though called slightly differently usually). I just think that your number of values in a nhalo+1 to nx-nhalo sized loop is not very big, so exiting a bit early is not saving you much.
Or, could we do this with an "any" call like
if( any(x0(nhalo+1:nx-nhalo) == 0.0d0)) then y_axis_thru_edges = .true.
|
I ran the PR through Claude to see what the AI would suggest. I saved it as a patch, so it refers to each commit as a "patch." I'll post the review below: ✅ What's good
🔴 Issues to fix before merge
This should almost certainly be flipped back to .false. before merging — verbose output at every timestep will flood production logs and slow runs. The commented-out line suggests it was left on accidentally during development. This redistributes small_dthck = 5.0 m of ice even when total_flux == 0 (no downstream ocean neighbors). In that case the block is guarded by if (total_flux > 0.0d0), so the ice would be removed from the CF cell but not redistributed anywhere — a mass conservation bug. It's a corner case (CF cell with no ocean edge-neighbor), but worth a check.
🟡 Concerns / suggestions The commit message says this fixes corner-neighbor fluxes being missed, which makes sense — those corner cells could be grounded. But removing the floating_mask check allows calving to propagate thickness changes into grounded cells too. Is that intended? A comment explaining the reasoning would help reviewers and future maintainers. Fortran doesn't zero-initialize derived type members by default. These will only be set when which_calving == CF_ADVANCE_RETREAT_RATE, so they'll be uninitialized garbage in the NetCDF output for all other calving options. Add = 0.0d0 initializers. !TODO - Can I get the angle from the flow direction at the CF? (line ~395) These are fine to leave in as long as they're tracked in the issue tracker. If not, they risk being forgotten. |
2 participants
This commit contains code changes that improve CISM results in CalvingMIP experiments. This is the code version used for our final CalvingMIP submission in March 2026. Changes include:
There is one new config option:
The branch originated from the reprosum branch several months ago. Several commits were added to the reprosum branch before it was merged to main. I then merged the current version of main to this calvingmip branch, bringing along several newer commits. I verified that this merger does not change CalvingMIP results.