From c1319da6067d20c3cbdf418e8f66228e272e6ad0 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Wed, 26 Aug 2020 18:50:09 -0600 Subject: [PATCH] New options to force retreat, remove isthmuses, and expand calving mask I added several options to support ISMIP6-based Antarctic experiments with ice-shelf collapse. First, I modified the force_retreat option, which uses a real mask with values in the range [0,1] to force ice to thin or retreat. Previously, force_retreat was a logical option (T or F). Now, it is an integer option (0, 1, or 2). Option 0 => no forced retreat Option 1 => Ice is thinned or removed wherever ice_fraction_retreat_mask > 0. Option 2 => Floating ice is removed in cells with ice_fraction_retreat_mask exceeding a threshold value (0.01 by default), provided these cells are connected to the ocean through other cells with ice_fraction_retreat_mask above the threshold. Grounded ice is unaffected. Option 1 is the same as the old option force_retreat = T, used for ISMIP6 Greenland experiments with forced retreat of the calving front. However, the old convention was to remove ice from cells with a mask value < 1; the new convention is to remove ice from cells with a mask value > 0. The new option 2 is appropriate for ISMIP6 Antarctic experiments in which a hydrofracture mask is read in as a forcing file. Initial experiments using a retreat mask from CESM2 21st century simulations showed that several shelves become unstable, for example when a shelf is divided into two regions connected by a narrow isthmus one cell wide. I addressed this problem by adding an option to remove ice isthmuses. An isthmus is define as a floating or weakly grounded cell with either ice-free ocean or thin floating ice on each side: e.g. if cell(i,j) is bordered by open ocean in cells (i-1,j) and (i+1,j). When isthmus cells are removed after applying the retreat mask, a potentially unstable shelf is divided into two separate regions. Subroutine remove_icebergs then removes the unstable region downstream, leaving the more stable region upstream. Also, I added a new logical option 'expand_calving_mask', to be used in conjunction with the calving_mask option. If the new option is set to true, then calving_mask is expanded at model initialization to include all floating ice in selected basins. This floating ice will then calve immediately, leading to acceleration of grounded ice. By default, this option is applied to 13 of the 16 ISMIP6 Antarctic basins. It is not applied to the three basins that include most of the Ross, Filchner-Ronne and Amery ice shelves. To change the basins where floating ice is calved, the user must insert basin numbers by hand. Note that this option differs from the ABUMIP float-kill option, in that all ice is calved in cells that are initially filled with floating ice, but not in cells where the ice subsequently floats. In multi-century transient experiments at 8 km and 4 km, I found that when forcing retreat using the ISMIP6 shelf-collapse mask, velocities can reach several tens of km/yr, but usually return to reasonble values of < 10 km/yr after a time step or two. In experiments using the expanded calving mask, ice speeds in the first few years can exceed 100 km/yr. These usually return to realistic values within a few simulation years, but sometimes the simulation can crash. More work would be needed to make this option truly robust. This commit is BFB except when using the new options. For earlier experiments with force_retreat = T in the config file, we will need to set force_retreat = 1 and replace ice_fraction_retreat_mask with (1 - ice_fraction_retreat_mask) to reproduce the results. --- libglide/glide_setup.F90 | 38 +++- libglide/glide_types.F90 | 15 +- libglissade/glissade.F90 | 304 ++++++++++++++++++++++++++++--- libglissade/glissade_calving.F90 | 95 +++++++++- libglissade/glissade_masks.F90 | 191 ++++++++++++++++++- 5 files changed, 599 insertions(+), 44 deletions(-) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 635c3080..3b9063ad 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -728,6 +728,8 @@ subroutine handle_options(section, model) call GetValue(section,'calving_domain',model%options%calving_domain) call GetValue(section,'apply_calving_mask', model%options%apply_calving_mask) call GetValue(section,'remove_icebergs', model%options%remove_icebergs) + call GetValue(section,'remove_isthmuses', model%options%remove_isthmuses) + call GetValue(section,'expand_calving_mask', model%options%expand_calving_mask) call GetValue(section,'limit_marine_cliffs', model%options%limit_marine_cliffs) call GetValue(section,'cull_calving_front', model%options%cull_calving_front) call GetValue(section,'adjust_input_thickness', model%options%adjust_input_thickness) @@ -1347,7 +1349,24 @@ subroutine print_options(model) else call write_log(' Icebergs will not be removed') endif + + if (model%options%remove_isthmuses) then + if (.not.model%options%remove_icebergs) then + model%options%remove_icebergs = .true. + write(message,*) ' Setting remove_icebergs = T for stability when remove_isthmuses = T' + call write_log(message) + endif + call write_log(' Isthmuses will be removed') + endif + if (model%options%expand_calving_mask) then + if (model%options%whichcalving == CALVING_GRID_MASK .or. model%options%apply_calving_mask) then + call write_log(' The calving mask will be expanded to include floating ice in select basins') + else + call write_log(' Not using a calving_mask; expand_calving_mask = T will be ignored') + endif + endif + if (model%options%limit_marine_cliffs) then call write_log(' The thickness of marine ice cliffs will be limited') call write_log(message) @@ -1883,20 +1902,23 @@ subroutine print_options(model) if (model%options%block_inception) then write(message,*) 'Inception outside the main ice sheet will be blocked' - else - write(message,*) 'Inception outside the main ice sheet is allowed' + call write_log(message) endif if (model%options%remove_ice_caps) then write(message,*) 'Ice caps will be removed and added to the calving flux' - else - write(message,*) 'Ice caps will not be removed' + call write_log(message) endif - if (model%options%force_retreat) then + if (model%options%force_retreat == FORCE_RETREAT_ALL_ICE) then write(message,*) 'Ice retreat will be forced using ice_fraction_retreat_mask' - else - write(message,*) 'Ice retreat will not be forced' + call write_log(message) + elseif (model%options%force_retreat == FORCE_RETREAT_FLOATING_ICE) then + write(message,*) 'Floating ice retreat will be forced using ice_fraction_retreat_mask' + call write_log(message) + if (.not.model%options%remove_isthmuses) then + call write_log(' Warning: Can be unstable when remove_isthmuses = F') + endif endif write(message,*) 'ho_whichice_age : ',model%options%which_ho_ice_age, & @@ -3238,7 +3260,7 @@ subroutine define_glide_restart_variables(options) ! If forcing ice retreat, then we need ice_fraction_retreat_mask (which specifies the cells where retreat is forced) ! and reference_thck (which sets up an upper thickness limit for partly retreating cells) - if (options%force_retreat) then + if (options%force_retreat /= FORCE_RETREAT_NONE) then call glide_add_to_restart_variable_list('ice_fraction_retreat_mask') call glide_add_to_restart_variable_list('reference_thck') endif diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index fa320e16..ffa3e9b4 100644 --- a/libglide/glide_types.F90 +++ b/libglide/glide_types.F90 @@ -193,6 +193,10 @@ module glide_types integer, parameter :: CALVING_DOMAIN_OCEAN_EDGE = 0 integer, parameter :: CALVING_DOMAIN_EVERYWHERE = 1 + integer, parameter :: FORCE_RETREAT_NONE = 0 + integer, parameter :: FORCE_RETREAT_ALL_ICE = 1 + integer, parameter :: FORCE_RETREAT_FLOATING_ICE = 2 + integer, parameter :: VERTINT_STANDARD = 0 integer, parameter :: VERTINT_KINEMATIC_BC = 1 @@ -634,6 +638,15 @@ module glide_types !> These are connected regions with zero basal traction and no connection to grounded ice. !> Safer to make it true, but not necessary for all applications + logical :: remove_isthmuses = .false. + !> if true, then identify and remove ice isthmuses after calving + !> These are narrow bridges connecting two regions of floating ice. + !> False by default, but may need to be true for the FORCE_RETREAT_FLOATING_ICE option. + + logical :: expand_calving_mask = .false. + !> if true, then expand the calving mask to include all ice that is floating at initialization + !> Note: By default, this is done for a hardwired set of ISMIP6 basins, excluding large shelves + logical :: limit_marine_cliffs = .false. !> if true, then thin marine-based cliffs based on a thickness threshold @@ -1011,7 +1024,7 @@ module glide_types logical :: remove_ice_caps = .false. !> Flag that indicates whether ice caps are removed and added to the calving flux - logical :: force_retreat = .false. + integer :: force_retreat = 0 !> Flag that indicates whether retreat is forced using ice_fraction_retreat_mask integer :: which_ho_ice_age = 1 diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 7784c07f..c5355cca 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -957,7 +957,7 @@ subroutine glissade_initialise(model, evolve_ice) ! If using a mask to force ice retreat, then set the reference thickness (if not already read in). - if (model%options%force_retreat) then + if (model%options%force_retreat /= FORCE_RETREAT_NONE) then ! Set reference_thck ! This field is loaded if present in the input file. Otherwise, reference_thck is set to the initial thickness. @@ -972,10 +972,20 @@ subroutine glissade_initialise(model, evolve_ice) call write_log(trim(message)) endif + ! Check whether a nonzero retreat mask has been read in. If not, then write a message. + ! Note: This is not necessarily an error. If reading the retreat mask from a forcing file, + ! the first nonzero mask may not be read until the model time is greater than tstart. + local_maxval = maxval(model%geometry%ice_fraction_retreat_mask) + global_maxval = parallel_reduce_max(local_maxval) + if (global_maxval < eps11) then + call write_log('Initial ice_fraction_retreat_mask = 0 everywhere') + endif + !WHL - debug if (this_rank == rtest) then i = itest j = jtest + print*, 'force_retreat option =', model%options%force_retreat print*, ' ' print*, 'reference_thck (m):' do j = jtest+3, jtest-3, -1 @@ -984,6 +994,14 @@ subroutine glissade_initialise(model, evolve_ice) enddo write(6,*) ' ' enddo + print*, ' ' + print*, 'ice_fraction_retreat_mask:' + do j = jtest+3, jtest-3, -1 + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') model%geometry%ice_fraction_retreat_mask(i,j) + enddo + write(6,*) ' ' + enddo endif endif ! force_retreat @@ -1043,7 +1061,7 @@ subroutine glissade_initialise(model, evolve_ice) model%calving%calving_front_x, model%calving%calving_front_y, & model%calving%calving_mask) - endif + endif ! calving grid mask ! Note: The DIVA solver needs a halo update for effective viscosity. ! This is done at the end of glissade_diagnostic_variable_solve, which in most cases is sufficient. @@ -2885,8 +2903,9 @@ subroutine glissade_calving_solve(model, init_calving) use glimmer_paramets, only: thk0, tim0, len0 use glissade_calving, only: glissade_calve_ice, glissade_cull_calving_front, & - glissade_remove_icebergs, glissade_limit_cliffs, verbose_calving - use glissade_masks, only: glissade_get_masks, glissade_calving_front_mask + glissade_remove_icebergs, glissade_remove_isthmuses, glissade_limit_cliffs, verbose_calving + use glissade_masks, only: glissade_get_masks, glissade_calving_front_mask, & + glissade_ocean_connection_mask use glissade_grounding_line, only: glissade_grounded_fraction implicit none @@ -2908,6 +2927,11 @@ subroutine glissade_calving_solve(model, init_calving) active_ice_mask, & ! = 1 if ice is present and dynamically active calving_front_mask ! = 1 for calving-front cells + integer, dimension(model%general%ewn, model%general%nsn) :: & + ocean_connection_mask, & ! = 1 for cells that are masked for retreat and are connected to the ocean + ! through other cells that are masked for retreat + retreat_mask ! local version of ice_fraction_retreat_mask; excludes grounded cells + real(dp) :: & maxthck, & ! max thickness of retreating ice dthck ! thickness loss for retreating ice @@ -2919,10 +2943,20 @@ subroutine glissade_calving_solve(model, init_calving) integer :: nx, ny ! horizontal grid dimensions integer :: itest, jtest, rtest ! coordinates of diagnostic point + real(dp), parameter :: & + retreat_mask_threshold = 0.01d0 ! threshold value for removing cells based on ice_fraction_retreat_mask; + ! set to a low value by default + ! Could make this a config parameter + + ! variables to expand the calving mask at initialization + logical, dimension(16) :: mask_basin ! true for basins whose floating ice is added to the calving mask + ! currently hardwired to 16 for ISMIP6 + integer :: bn ! basin number + !WHL - debug - real(dp) :: local_maxval, global_maxval logical, parameter :: verbose_retreat = .true. + nx = model%general%ewn ny = model%general%nsn @@ -2943,25 +2977,25 @@ subroutine glissade_calving_solve(model, init_calving) ! Note: This option is similar to apply_calving_mask. It is different in that the mask ! is a real number in the range [0,1], allowing thinning instead of complete removal. ! Do not thin or remove ice if this is the initial calving call; force retreat only during runtime. - ! This option is applied before calling glissade_calve_ice, so that ice thinned by the retreat mask + ! There are two forced retreat options: + ! Option 1: Thin or remove ice wherever ice_fraction_retreat_mask > 0. + ! Option 2: Remove floating ice (but not grounded ice) where ice_fraction_retreat_mask > 0. + ! + ! Option 1 is done now, before calling glissade_calve_ice, so that ice thinned by the retreat mask ! can undergo further thinning or removal by the calving scheme. + ! Option 2 is done after the main calving solve, after thin ice at the calving front has been removed + ! by other mechanisms. - if (model%options%force_retreat .and. .not.init_calving) then - if (main_task) print*, 'Forcing retreat using ice_fraction_retreat_mask, time =', model%numerics%time - - !WHL - debug - ! Check that a nonzero retreat mask has been read in. Otherwise, all ice would erroneously be removed. - local_maxval = maxval(model%geometry%ice_fraction_retreat_mask) - global_maxval = parallel_reduce_max(local_maxval) - if (global_maxval < eps11) then - call write_log('Error: Failed to read in a nonzero retreat mask with force_retreat = T', GM_FATAL) + if (model%options%force_retreat == FORCE_RETREAT_ALL_ICE .and. .not.init_calving) then + if (this_rank == rtest) then + print*, 'Forcing retreat using ice_fraction_retreat_mask, time =', model%numerics%time endif if (verbose_retreat .and. this_rank == rtest) then i = itest j = jtest print*, ' ' - print*, 'Before front retreat, thck (m):' + print*, 'Before forced retreat, thck (m):' do j = jtest+3, jtest-3, -1 do i = itest-3, itest+3 write(6,'(f10.3)',advance='no') model%geometry%thck(i,j)*thk0 @@ -2969,7 +3003,7 @@ subroutine glissade_calving_solve(model, init_calving) write(6,*) ' ' enddo print*, ' ' - print*, 'Apply front retreat, retreat_mask:' + print*, 'ice_fraction_retreat_mask:' do j = jtest+3, jtest-3, -1 do i = itest-3, itest+3 write(6,'(f10.3)',advance='no') model%geometry%ice_fraction_retreat_mask(i,j) @@ -2980,12 +3014,8 @@ subroutine glissade_calving_solve(model, init_calving) print*, 'maxthck (m):' do j = jtest+3, jtest-3, -1 do i = itest-3, itest+3 - if (model%geometry%ice_fraction_retreat_mask(i,j) < 1.0d0) then - write(6,'(f10.3)',advance='no') & - model%geometry%reference_thck(i,j)*thk0 * model%geometry%ice_fraction_retreat_mask(i,j) - else - write(6,'(f10.3)',advance='no') model%geometry%thck(i,j)*thk0 - endif + write(6,'(f10.3)',advance='no') model%geometry%reference_thck(i,j)*thk0 & + * (1.0d0 - model%geometry%ice_fraction_retreat_mask(i,j)) enddo write(6,*) ' ' enddo @@ -2993,8 +3023,9 @@ subroutine glissade_calving_solve(model, init_calving) do j = 1, model%general%nsn do i = 1, model%general%ewn - if (model%geometry%ice_fraction_retreat_mask(i,j) < 1.0d0) then - maxthck = model%geometry%reference_thck(i,j) * model%geometry%ice_fraction_retreat_mask(i,j) + if (model%geometry%ice_fraction_retreat_mask(i,j) > 0.0d0) then + maxthck = model%geometry%reference_thck(i,j) & + * (1.0d0 - model%geometry%ice_fraction_retreat_mask(i,j)) dthck = model%geometry%thck(i,j) - min(maxthck, model%geometry%thck(i,j)) model%geometry%thck(i,j) = model%geometry%thck(i,j) - dthck model%calving%calving_thck(i,j) = model%calving%calving_thck(i,j) + dthck @@ -3006,7 +3037,7 @@ subroutine glissade_calving_solve(model, init_calving) i = itest j = jtest print*, ' ' - print*, 'After front retreat, thck (m):' + print*, 'After forced retreat, thck (m):' do j = jtest+3, jtest-3, -1 do i = itest-3, itest+3 write(6,'(f10.3)',advance='no') model%geometry%thck(i,j)*thk0 @@ -3015,7 +3046,7 @@ subroutine glissade_calving_solve(model, init_calving) enddo endif - endif ! force_retreat + endif ! force_retreat_all_ice !TODO - Make sure no additional halo updates are needed before glissade_calve_ice @@ -3030,6 +3061,69 @@ subroutine glissade_calving_solve(model, init_calving) if (model%options%whichcalving == CALVING_GRID_MASK .or. model%options%apply_calving_mask) then + ! Optionally, expand the calving mask to include floating ice in select basins. + ! Note: Currently hardwired to include 13 of the 16 ISMIP6 basins. + ! Does not include the three largest shelves (Ross, Filchner-Ronne, Amery) + + if (init_calving .and. model%options%expand_calving_mask) then + + ! Identify basins whose floating ice will be added to the calving mask + ! Currently hardwired to the ISMIP6 basin numbers (1 to 16) + mask_basin(:) = .true. + mask_basin(2) = .false. ! Amery + mask_basin(7) = .false. ! Ross + mask_basin(14) = .false. ! Filchner-Ronne + + if (verbose_calving .and. this_rank==rtest) then + print*, 'Expanding the calving mask to ice shelves in select basins' + print*, 'basin number, mask_basin:' + do bn = 1, 16 + print*, bn, mask_basin(bn) + enddo + endif + + call glissade_get_masks(nx, ny, & + model%geometry%thck*thk0, model%geometry%topg*thk0, & + model%climate%eus*thk0, 0.0d0, & ! thklim = 0 + ice_mask, & + floating_mask = floating_mask, & + land_mask = land_mask) + + if (verbose_calving .and. this_rank==rtest) then + print*, ' ' + print*, 'initial calving_mask, itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(i10)',advance='no') model%calving%calving_mask(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'floating_mask, itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(i10)',advance='no') floating_mask(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! For basins with mask_basin = T, add floating ice to the calving mask. + do j = 1, model%general%nsn + do i = 1, model%general%ewn + bn = model%ocean_data%basin_number(i,j) + if (mask_basin(bn) .and. floating_mask(i,j) == 1) then + model%calving%calving_mask(i,j) = 1 + endif + enddo + enddo + + call parallel_halo(model%calving%calving_mask) + + endif ! expand_calving_mask + if (verbose_calving .and. this_rank==rtest) then print*, ' ' print*, 'Limit advance of calving front' @@ -3043,6 +3137,15 @@ subroutine glissade_calving_solve(model, init_calving) write(6,*) ' ' enddo print*, ' ' + print*, 'floating_mask, itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(i10)',advance='no') floating_mask(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' print*, 'calving_mask, itest, jtest, rank =', itest, jtest, rtest do j = jtest+3, jtest-3, -1 write(6,'(i6)',advance='no') j @@ -3114,6 +3217,153 @@ subroutine glissade_calving_solve(model, init_calving) endif + if (model%options%force_retreat == FORCE_RETREAT_FLOATING_ICE) then + + ! Remove floating ice based on ice_fraction_retreat_mask. + ! This is done after the main calving routine, to avoid complications + ! involving thin ice near the calving front that calves after transport. + ! The logic works as follows: + ! * Idenfity cells with ice_fraction_retreat_mask exceeding some threshold. + ! * Remove any such cells if they are adjacent to ocean cells, or are connected + ! to the ocean through other identified cells. + ! * Do not remove cells without a connection to the ocean. + ! In other words, do not hollow out ice shelves from the interior, since + ! this can be numerically unstable. + + ! Update masks + + call glissade_get_masks(nx, ny, & + thck_unscaled, model%geometry%topg*thk0, & + model%climate%eus*thk0, model%numerics%thklim*thk0, & + ice_mask, & + floating_mask = floating_mask, & + ocean_mask = ocean_mask) + + ! Identify floating cells with ice_fraction_retreat_mask exceeding a prescribed threshold. + + where (floating_mask == 1 .and. model%geometry%ice_fraction_retreat_mask > retreat_mask_threshold) + retreat_mask = 1 + elsewhere + retreat_mask = 0 + endwhere + + ! Identify cells that have retreat_mask = 1 and are either adjacent to ocean cells, + ! or are connected to the ocean through other cells with retreat_mask = 1. + + call glissade_ocean_connection_mask(nx, ny, & + itest, jtest, rtest, & + thck_unscaled, retreat_mask, & + ocean_mask, & + ocean_connection_mask) + + if (verbose_calving .and. this_rank==rtest) then + print*, ' ' + print*, 'Force floating ice retreat' + print*, ' ' + print*, 'Thickness before retreat, itest, jtest, rtest =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f7.1)',advance='no') thck_unscaled(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'floating_mask:' + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(i7)',advance='no') floating_mask(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'ocean_mask:' + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-7, itest+3 + write(6,'(i7)',advance='no') ocean_mask(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'ice_fraction_retreat_mask:' + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f7.2)',advance='no') model%geometry%ice_fraction_retreat_mask(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'ocean_connection_mask:' + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(i7)',advance='no') ocean_connection_mask(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! Remove ice from ocean-connected cells with retreat_mask = 1 + where (ocean_connection_mask == 1) + model%calving%calving_thck = model%calving%calving_thck + thck_unscaled + thck_unscaled = 0.0d0 + !TODO - Reset temperature and other tracers in cells where the ice calved? + endwhere + + endif ! force_retreat_floating_ice + + if (model%options%remove_isthmuses) then + + ! Optionally, remove isthmuses. + ! An isthmus is defined as a floating or weakly grounded grid cell with ice-free ocean + ! or thin floating ice on both sides. + ! When using a calving or retreat mask derived from an ESM or other model, + ! isthmuses may need to be removed to prevent unstable ice configurations, + ! e.g. a shelf split into two parts connected by a bridge one cell wide. + ! Isthmus removal should always be followed by iceberg removal. + + ! Update the masks + + call glissade_get_masks(nx, ny, & + thck_unscaled, model%geometry%topg*thk0, & + model%climate%eus*thk0, model%numerics%thklim*thk0, & + ice_mask, & + floating_mask = floating_mask, & + ocean_mask = ocean_mask, & + land_mask = land_mask) + + ! Compute f_ground_cell for isthmus removal + + call glissade_grounded_fraction(nx, ny, & + itest, jtest, rtest, & ! diagnostic only + thck_unscaled, & + model%geometry%topg*thk0, & + model%climate%eus*thk0, & + ice_mask, & + floating_mask, & + land_mask, & + model%options%which_ho_ground, & + model%options%which_ho_flotation_function, & + model%options%which_ho_fground_no_glp, & + model%geometry%f_flotation, & + model%geometry%f_ground, & + model%geometry%f_ground_cell, & + model%geometry%topg_stdev*thk0) + + call glissade_remove_isthmuses(& + nx, ny, & + itest, jtest, rtest, & + thck_unscaled, & + model%geometry%f_ground_cell, & + floating_mask, & + ocean_mask, & + model%calving%calving_thck) + + endif ! remove isthmuses + ! ------------------------------------------------------------------------ ! Remove any icebergs. ! For the velocity solver to be robust, we require that any floating cell diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index 9d1b7ef1..062dc646 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -41,7 +41,8 @@ module glissade_calving private public :: glissade_calving_mask_init, glissade_thck_calving_threshold_init, & glissade_calve_ice, glissade_cull_calving_front, & - glissade_remove_icebergs, glissade_limit_cliffs + glissade_remove_icebergs, glissade_remove_isthmuses, & + glissade_limit_cliffs public :: verbose_calving logical, parameter :: verbose_calving = .false. @@ -1570,6 +1571,98 @@ subroutine glissade_remove_icebergs(& end subroutine glissade_remove_icebergs +!--------------------------------------------------------------------------- + + subroutine glissade_remove_isthmuses(& + nx, ny, & + itest, jtest, rtest, & + thck, & + f_ground_cell, & + floating_mask, & + ocean_mask, & + calving_thck) + + ! Remove any ice isthmuses. + ! An isthmus is defined as a floating or weakly grounded grid cell with ice-free ocean + ! or thin floating ice (H < 10 m) on both sides: i.e., in cells (i-1,j) and (i+1,j), + ! or (i,j-1) and (i,j+1). + ! When using an ice_fraction_retreat_mask derived from a model (e.g., CESM), + ! it may be necessary to remove isthmuses to prevent unstable ice configurations, + ! e.g. a shelf split into two parts connected by a bridge one cell wide. + ! Isthmus removal should always be followed by iceberg removal. + + integer :: nx, ny !> horizontal grid dimensions + integer, intent(in) :: itest, jtest, rtest !> coordinates of diagnostic point + + real(dp), dimension(nx,ny), intent(inout) :: thck !> ice thickness (m) + real(dp), dimension(nx,ny), intent(in) :: f_ground_cell !> grounded fraction in each grid cell + integer, dimension(nx,ny), intent(in) :: floating_mask !> = 1 ice is present and floating, else = 0 + integer, dimension(nx,ny), intent(in) :: ocean_mask !> = 1 where topg is below sea level and ice is absent, else = 0 + real(dp), dimension(nx,ny), intent(inout) :: calving_thck !> thickness (m) lost due to calving in each grid cell; + !> on output, includes ice removed from isthmuses + + ! local variables + + integer :: i, j + + integer, dimension(nx,ny) :: & + ocean_plus_thin_ice_mask ! = 1 for ocean cells and cells with thin floating ice + + ! Both floating and weakly grounded cells can be identified as isthmuses and removed; + ! isthmuses_f_ground_threshold is used to identify weakly grounded cells. + real(dp), parameter :: & ! threshold for counting cells as grounded + isthmus_f_ground_threshold = 0.50d0 + + ! An isthmus cell has ice-free ocean or thin floating ice on each side: + ! isthmus_f_ground_threshold is used to identify thin floating ice. + real(dp), parameter :: & ! threshold (m) for counting floating ice as thin + isthmus_thck_threshold = 10.0d0 + + if (verbose_calving .and. this_rank == rtest) then + print*, ' ' + print*, 'In glissade_remove_isthmuses' + print*, ' ' + print*, 'Thickness before isthmus removal:' + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') thck(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ocean_plus_thin_ice_mask = ocean_mask + where (floating_mask == 1 .and. thck < isthmus_thck_threshold) + ocean_plus_thin_ice_mask = 1 + endwhere + + do j = 2, ny-1 + do i = 2, nx-1 + if (floating_mask(i,j) == 1 .or. f_ground_cell(i,j) < isthmus_f_ground_threshold) then + if ( (ocean_plus_thin_ice_mask(i-1,j) == 1 .and. ocean_plus_thin_ice_mask(i+1,j) == 1) .or. & + (ocean_plus_thin_ice_mask(i,j-1) == 1 .and. ocean_plus_thin_ice_mask(i,j+1) == 1) ) then + calving_thck(i,j) = calving_thck(i,j) + thck(i,j) + thck(i,j) = 0.0d0 + endif + endif + enddo + enddo + + if (verbose_calving .and. this_rank==rtest) then + print*, ' ' + print*, 'Thickness after isthmus removal' + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') thck(i,j) + enddo + write(6,*) ' ' + enddo + endif + + end subroutine glissade_remove_isthmuses + !--------------------------------------------------------------------------- subroutine glissade_limit_cliffs(& diff --git a/libglissade/glissade_masks.F90 b/libglissade/glissade_masks.F90 index 1a3baf2e..916e2179 100644 --- a/libglissade/glissade_masks.F90 +++ b/libglissade/glissade_masks.F90 @@ -47,10 +47,12 @@ module glissade_masks implicit none private - public :: glissade_get_masks, glissade_calving_front_mask, & + public :: glissade_get_masks, glissade_calving_front_mask, & glissade_marine_cliff_mask, glissade_ice_sheet_mask, & + glissade_ocean_connection_mask, & glissade_marine_connection_mask, glissade_lake_mask, & glissade_extend_mask, glissade_fill, glissade_fill_with_buffer + public :: initial_color, fill_color, boundary_color ! colors for fill subroutines @@ -752,18 +754,193 @@ subroutine glissade_ice_sheet_mask(nx, ny, & end subroutine glissade_ice_sheet_mask +!**************************************************************************** + + subroutine glissade_ocean_connection_mask(nx, ny, & + itest, jtest, rtest, & + thck, input_mask, & + ocean_mask, & + ocean_connection_mask) + + ! Create a masks consisting of cells with input_mask = 1, which either + ! are adjacent to the ocean, or are connected to the ocean through other cells + ! with input_mask = 1. + + ! The algorithm is as follows: + ! (1) Assign the initial color to cells with input_mask = 1. + ! Mark other cells with the boundary color. + ! (2) Seed the fill by giving the fill color to cells that have input_mask = 1 + ! and are adjacent to ocean cells (ocean_mask = 1). + ! (3) Recursively fill all cells that have input_mask = 1 and are connected + ! to the ocean by a path that passes through other cells with input_mask = 1. + ! (4) Repeat the recursion as necessary to spread the fill to adjacent processors. + ! (5) Once the fill is done, all cells with the fill color are assigned + ! to ocean_connection_mask. + ! + ! Note: The logic is general enough that we could replace ocean_mask and + ! ocean_connection_mask with generic input and output masks. + + integer, intent(in) :: nx, ny !> horizontal grid dimensions + + integer, intent(in) :: itest, jtest, rtest !> coordinates of diagnostic point + + real(dp), dimension(nx,ny), intent(in) :: & + thck !> ice thickness (m) + + integer, dimension(nx,ny), intent(in) :: & + input_mask, & !> = 1 for cells that meet some criterion specified elsewhere + ocean_mask !> = 1 for ice-free cells with topg below sea level + + integer, dimension(nx,ny), intent(out) :: & + ocean_connection_mask !> = 1 for cells with input_mask = 1, connected to the ocean + !> through other such cells + + ! local variables + + integer :: i, j, iter + + integer :: & + max_iter, & ! max(ewtasks, nstasks) + local_count, & ! local counter for filled values + global_count, & ! global counter for filled values + global_count_save ! globalcounter for filled values from previous iteration + + integer, dimension(nx,ny) :: & + color ! color variable for the fill + + logical, parameter :: verbose_ocean_connection_mask = .false. + + ! initialize + ! Note: Cells with input_mask = 1 receive the initial color, and other cells receive the boundary color. + + where (input_mask == 1) + color = initial_color + elsewhere + color = boundary_color + endwhere + + ! Identifying cells that have the initial color and are adjacent to the ocean. + ! Fill these cells and then recursively fill neighbors that have the initial color. + ! We may have to do this several times to incorporate connections between neighboring processors. + + max_iter = max(ewtasks,nstasks) + global_count_save = 0 + + do iter = 1, max_iter + + if (iter == 1) then ! identify cells adjacent to the ocean, which can seed the fill + + do j = 2, ny-1 + do i = 2, nx-1 + if (color(i,j) == initial_color .and. & + (ocean_mask(i-1,j) == 1 .or. ocean_mask(i+1,j) == 1 .or. & + ocean_mask(i,j-1) == 1 .or. ocean_mask(i,j+1) == 1) ) then + ! assign the fill color to this cell, and recursively fill neighbors with input_mask = 1 + call glissade_fill(nx, ny, & + i, j, & + color, input_mask) + endif + enddo + enddo + + else ! iter > 1 + + ! Check for halo cells that were just filled on neighbor processors + ! Note: In order for a halo cell to seed the fill on this processor, it must already have the fill color. + + call parallel_halo(color) + + ! west halo layer + i = nhalo + do j = 1, ny + if (color(i,j) == fill_color) then + call glissade_fill(nx, ny, & + i+1, j, & + color, input_mask) + endif + enddo + + ! east halo layers + i = nx - nhalo + 1 + do j = 1, ny + if (color(i,j) == fill_color) then + call glissade_fill(nx, ny, & + i-1, j, & + color, input_mask) + endif + enddo + + ! south halo layer + j = nhalo + do i = nhalo+1, nx-nhalo ! already checked halo corners above + if (color(i,j) == fill_color) then + call glissade_fill(nx, ny, & + i, j+1, & + color, input_mask) + endif + enddo + + ! north halo layer + j = ny-nhalo+1 + do i = nhalo+1, nx-nhalo ! already checked halo corners above + if (color(i,j) == fill_color) then + call glissade_fill(nx, ny, & + i, j-1, & + color, input_mask) + endif + enddo + + endif ! iter = 1 + + ! Count the number of filled cells. If converged, then exit the loop. + + local_count = 0 + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo + if (color(i,j) == fill_color) local_count = local_count + 1 + enddo + enddo + + global_count = parallel_reduce_sum(local_count) + + if (global_count == global_count_save) then + if (verbose_ocean_connection_mask .and. main_task) & + print*, 'ocean_connection_mask, fill converged: iter, global_count =', iter, global_count + exit + else + if (verbose_ocean_connection_mask .and. main_task) & + print*, 'ocean_connection_mask, convergence check: iter, global_count =', iter, global_count + global_count_save = global_count + endif + + enddo ! max_iter + + ! Any cells with the fill color are deemed to be ocean-connected. + + call parallel_halo(color) + + ocean_connection_mask(:,:) = 0 + + where (color == fill_color) + ocean_connection_mask = 1 + elsewhere + ocean_connection_mask = 0 + endwhere + + end subroutine glissade_ocean_connection_mask + !**************************************************************************** subroutine glissade_marine_connection_mask(nx, ny, & itest, jtest, rtest, & - thck, topg, & - eus, thklim, & + thck, topg, & + eus, thklim, & marine_connection_mask) - ! Identify cells that have a marine path to the ocean. - ! The path can include grounded marine-based ice. - ! The ocean is identified as ice-free cells with bed elevation below sea level, - ! or optionally with bed elevation below some threshold value. + ! Identify cells that have a marine path to the ocean. + ! The path can include grounded marine-based ice. + ! The ocean is identified as ice-free cells with bed elevation below sea level, + ! or optionally with bed elevation below some threshold value. ! subroutine arguments