Skip to content

Commit

Permalink
Reverted to previous hybrid gradient scheme
Browse files Browse the repository at this point in the history 
In a recent commit (89f6875), I modified the hybrid gradient scheme to compute gradients
at edges where an ice-covered land cell sits above ice-free ocean. On further study,
this seems not to be a good idea. When an ice-covered cell at the land margin sits
above ice-free ocean, surface elevation gradients can be large. Even when gradients
are limited by max_slope, the rhs forcing term can become very large when the coastal
ice is several hundred meters thick, possibly promoting instability.

With this commit, gradients are set to zero where an ice-covered land cell sits above ice-free ocean,
as before the recent change. Instead, we compute a lateral spreading term at these edges,
assuming that the margin is a vertical cliff whose height can be approximated by stagthck.
The lateral spreading term can be large, but tends to be smaller than the gradient term would be.

This commit is BFB for standard LIVV tests, but gives different answers for
Greenland simulations, especially where land promontories border ice-free fjords.
  • Loading branch information
@whlipscomb
whlipscomb committed May 13, 2018
1 parent 0ba5870 commit 14f5f8d
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 65 deletions.
20 changes: 3 additions & 17 deletions libglissade/glissade_grid_operators.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -657,10 +657,6 @@ subroutine glissade_edgemask_gradient_margin_hybrid(nx, ny, &
! At land-terminating margins the gradient is nonzero (except for nunataks), and at marine-terminating
! margins the gradient is zero.
!
! TODO: Update this subroutine to be consistent with the logic in glissade_surface_elevation_gradient.
! I.e., add the case of ice-covered land over ice-free ocean, and incorporate thck_gradient_ramp.
! For now, this subroutine is used only by the glissade SIA solver (which is not run with marine
! boundaries) and the L1L2 solver (which is no longer supported), so the differences should not matter.
!----------------------------------------------------------------

!----------------------------------------------------------------
Expand Down Expand Up @@ -764,15 +760,15 @@ subroutine glissade_surface_elevation_gradient(nx, ny, &
! for ice shelves with a sharp drop in ice thickness and surface elevation at the margin.
!
! HO_GRADIENT_MARGIN_HYBRID = 1: The gradient is computed at edges where either
!
! (1) Both adjacent cells are ice-covered.
! (2) One cell is ice-covered (land or marine-based) and lies above ice-free land.
! (3) One cell is ice-covered land and lies above ice-free ocean.
!
! This method sets the gradient to zero at edges where
! (1) An ice-covered marine-based cell (grounded or floating) lies above ice-free ocean.
! (1) An ice-covered cell (grounded or floating) lies above ice-free ocean.
! Note: Inactive calving-front cells are treated as ice-free ocean.
! (2) An ice-covered land cell lies below an ice-free land cell (i.e., a nunatak).
!

! The aim is to give a reasonable gradient at both land-terminating and marine-terminating margins.
! At land-terminating margins the gradient is nonzero (except for nunataks), and at marine-terminating
! margins the gradient is zero.
Expand Down Expand Up @@ -934,11 +930,6 @@ subroutine glissade_surface_elevation_gradient(nx, ny, &
! upper cell has active ice, and ice-free lower cell is land; compute the gradient
ds_dx_edge(i,j) = edge_factor * sign_factor * (usrf(iu,j) - usrf(il,j)) / dx

elseif (active_ice_mask(iu,j) == 1 .and. land_mask(iu,j) == 1 .and. land_mask(il,j) == 0) then

! upper cell has active ice on land, and ice-free lower cell is ocean; compute the gradient
ds_dx_edge(i,j) = edge_factor * sign_factor * (usrf(iu,j) - usrf(il,j)) / dx

endif ! both cells have ice

enddo ! i
Expand Down Expand Up @@ -979,11 +970,6 @@ subroutine glissade_surface_elevation_gradient(nx, ny, &
! upper cell has active ice, and ice-free lower cell is land; compute the gradient
ds_dy_edge(i,j) = edge_factor * sign_factor * (usrf(i,ju) - usrf(i,jl)) / dy

elseif (active_ice_mask(i,ju) == 1 .and. land_mask(i,ju) == 1 .and. land_mask(i,jl) == 0) then

! upper cell has active ice on land, and ice-free lower cell is ocean; compute the gradient
ds_dy_edge(i,j) = edge_factor * sign_factor * (usrf(i,ju) - usrf(i,jl)) / dy

endif ! both cells have ice

enddo ! i
Expand Down
70 changes: 22 additions & 48 deletions libglissade/glissade_velo_higher.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -812,8 +812,6 @@ subroutine glissade_velo_higher_solve(model, &
xVertex, yVertex, & ! x and y coordinates of each vertex (m)
stagusrf, & ! upper surface averaged to vertices, for active cells (m)
stagthck, & ! ice thickness averaged to vertices, for active cells (m)
stagusrf_marine, & ! upper surface averaged to vertices, for active marine cells only (m)
stagthck_marine, & ! ice thickness averaged to vertices, for active marine cells only (m)
dusrf_dx, dusrf_dy, & ! gradient of upper surface elevation (m/m)
ubas, vbas ! basal ice velocity (m/yr); input to calcbeta

Expand Down Expand Up @@ -1531,18 +1529,6 @@ subroutine glissade_velo_higher_solve(model, &
active_marine_mask = 0
endwhere

! Compute marine version of stagthck and stagusrf, used for lateral shelf BCs

call glissade_stagger(nx, ny, &
thck, stagthck_marine, &
active_marine_mask, &
stagger_margin_in = 1)

call glissade_stagger(nx, ny, &
usrf, stagusrf_marine, &
active_marine_mask, &
stagger_margin_in = 1)

if (verbose_gridop .and. this_rank == rtest) then
print*, ' '
print*, 'thck, itest, jtest, rank =', itest, jtest, rtest
Expand All @@ -1562,15 +1548,6 @@ subroutine glissade_velo_higher_solve(model, &
enddo
write(6,*) ' '
enddo
print*, ' '
print*, 'stagthck_marine, 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,'(f10.3)',advance='no') stagthck_marine(i,j)
enddo
write(6,*) ' '
enddo
endif

!------------------------------------------------------------------------------
Expand All @@ -1586,9 +1563,8 @@ subroutine glissade_velo_higher_solve(model, &
! if ice-free. This is what Glide does, but is not appropriate if we have ice-covered
! marine-based cells lying above ice-free ocean cells, because the gradient is too big.
! gradient_margin_in = 1 computes gradients at edges with
! (1) ice-covered cells on either side,
! (2) ice-covered cells (land or marine-based) above ice-free land, or
! (3) ice-covered land cells above ice-free ocean.
! (1) ice-covered cells on either side, or
! (2) ice-covered cell (land or marine-based) above ice-free land
! This option is designed for both land- and ocean-terminating boundaries. It is the default.
! gradient_margin_in = 2 computes gradients only at edges with ice-covered cells
! on each side. This is appropriate for problems with ice shelves, but is
Expand Down Expand Up @@ -2156,7 +2132,7 @@ subroutine glissade_velo_higher_solve(model, &
calving_front_mask, &
active_cell, &
xVertex, yVertex, &
stagusrf_marine, stagthck_marine, &
stagusrf, stagthck, &
loadu, loadv)
call t_stopf('glissade_load_vector_lateral_bc')

Expand Down Expand Up @@ -4604,7 +4580,7 @@ subroutine load_vector_lateral_bc(nx, ny, &
calving_front_mask, &
active_cell, &
xVertex, yVertex, &
stagusrf_marine, stagthck_marine, &
stagusrf, stagthck, &
loadu, loadv)

integer, intent(in) :: &
Expand All @@ -4629,8 +4605,8 @@ subroutine load_vector_lateral_bc(nx, ny, &
xVertex, yVertex ! x and y coordinates of vertices

real(dp), dimension(nx-1,ny-1), intent(in) :: &
stagusrf_marine, & ! upper surface elevation (m) on staggered grid, for marine-based cells only
stagthck_marine ! ice thickness (m) on staggered grid, for marine-based cells only
stagusrf, & ! upper surface elevation (m) on staggered grid
stagthck ! ice thickness (m) on staggered grid

real(dp), dimension(nz,nx-1,ny-1), intent(inout) :: &
loadu, loadv ! load vector, divided into u and v components
Expand All @@ -4644,8 +4620,7 @@ subroutine load_vector_lateral_bc(nx, ny, &
! Sum over elements in active cells
! Loop over cells that contain locally owned vertices

! Note: Lateral shelf BCs are applied to active cells (either floating cells or
! marine-based grounded cells) that border the ocean.
! Note: Lateral shelf BCs are applied to active cells (either floating or grounded) that border the ocean.
! Inactive calving_front cells are treated as if they were ocean cells.

do j = nhalo+1, ny-nhalo+1
Expand All @@ -4659,10 +4634,9 @@ subroutine load_vector_lateral_bc(nx, ny, &
print*, 'calving_front_mask (i-1:i,j-1)=', calving_front_mask(i-1:i, j-1)
endif

!WHL - Old method is to compute the spreading term only for active floating cells.
! New method is to compute the spreading term for all active marine-based cells that border the ocean.
!! if (active_cell(i,j) .and. floating_mask(i,j) == 1) then ! ice is active and floating
if (active_cell(i,j) .and. land_mask(i,j) == 0) then
! Compute the spreading term for all active cells that share an edge with an ice-free ocean cell.

if (active_cell(i,j)) then

if ( ocean_mask(i-1,j) == 1 .or. &
(calving_front_mask(i-1,j) == 1 .and. .not.active_cell(i-1,j)) ) then ! compute lateral BC for west face
Expand All @@ -4671,7 +4645,7 @@ subroutine load_vector_lateral_bc(nx, ny, &
nz, sigma, &
'west', &
i, j, &
stagusrf_marine, stagthck_marine, &
stagusrf, stagthck, &
xVertex, yVertex, &
loadu, loadv)
endif
Expand All @@ -4683,7 +4657,7 @@ subroutine load_vector_lateral_bc(nx, ny, &
nz, sigma, &
'east', &
i, j, &
stagusrf_marine, stagthck_marine, &
stagusrf, stagthck, &
xVertex, yVertex, &
loadu, loadv)
endif
Expand All @@ -4695,7 +4669,7 @@ subroutine load_vector_lateral_bc(nx, ny, &
nz, sigma, &
'south', &
i, j, &
stagusrf_marine, stagthck_marine, &
stagusrf, stagthck, &
xVertex, yVertex, &
loadu, loadv)
endif
Expand All @@ -4707,7 +4681,7 @@ subroutine load_vector_lateral_bc(nx, ny, &
nz, sigma, &
'north', &
i, j, &
stagusrf_marine, stagthck_marine, &
stagusrf, stagthck, &
xVertex, yVertex, &
loadu, loadv)
endif
Expand All @@ -4725,7 +4699,7 @@ subroutine lateral_shelf_bc(nx, ny, &
nz, sigma, &
face, &
iCell, jCell, &
stagusrf_marine, stagthck_marine, &
stagusrf, stagthck, &
xVertex, yVertex, &
loadu, loadv)

Expand Down Expand Up @@ -4779,8 +4753,8 @@ subroutine lateral_shelf_bc(nx, ny, &
xVertex, yVertex ! x and y coordinates of vertices

real(dp), dimension(nx-1,ny-1), intent(in) :: &
stagusrf_marine, & ! upper surface elevation (m) on staggered grid, for marine-based cells only
stagthck_marine ! ice thickness (m) on staggered grid (m), for marine-based cells only
stagusrf, & ! upper surface elevation (m) on staggered grid
stagthck ! ice thickness (m) on staggered grid (m)

real(dp), dimension(nz,nx-1,ny-1), intent(inout) :: &
loadu, loadv ! load vector, divided into u and v components
Expand Down Expand Up @@ -4894,13 +4868,13 @@ subroutine lateral_shelf_bc(nx, ny, &
x(3) = x(2)
x(4) = x(1)

s(1) = stagusrf_marine(iNode(1), jNode(1))
s(2) = stagusrf_marine(iNode(2), jNode(2))
s(1) = stagusrf(iNode(1), jNode(1))
s(2) = stagusrf(iNode(2), jNode(2))
s(3) = s(2)
s(4) = s(1)

h(1) = stagthck_marine(iNode(1), jNode(1))
h(2) = stagthck_marine(iNode(2), jNode(2))
h(1) = stagthck(iNode(1), jNode(1))
h(2) = stagthck(iNode(2), jNode(2))
h(3) = h(2)
h(4) = h(1)

Expand Down Expand Up @@ -6256,7 +6230,7 @@ subroutine compute_3d_velocity_L1L2(nx, ny, &
! Setting gradient_margin_in = 0 takes the gradient over both neighboring cells,
! including ice-free cells.
! Setting gradient_margin_in = 1 computes a gradient if both neighbor cells are
! either ice-covered cells or land cells; else gradient = 0.
! ice-covered, or an ice-covered cell sits above ice-free land; else gradient = 0
! Setting gradient_margin_in = 2 computes a gradient only if both neighbor cells
! are ice-covered.
! At a land margin, either 0 or 1 is appropriate, but 2 is inaccurate.
Expand Down

0 comments on commit 14f5f8d

Please sign in to comment.