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