Skip to content

Commit

Permalink
Merge pull request #16 from NCAR/lipscomb/gradient_fix
Browse files Browse the repository at this point in the history 
Revert to previous hybrid gradient scheme

In a recent commit (89f6875), the hybrid gradient scheme was modified
to compute gradients at edges where an ice-covered land cell sits above
ice-free ocean. This commit reverts to the earlier rule that gradients
are set to zero where an ice-covered land sits above ice-free
ocean. Instead of computing a nonzero gradient, we compute a lateral
spreading term at these edges, assuming that the margin is a
marine-based cliff. Answer changes are modest, but generally reduce max
speeds at coastal margins and promote greater stability.
  • Loading branch information
@billsacks
billsacks authored May 14, 2018
2 parents 0ba5870 + 14f5f8d commit c713c09
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 c713c09

Please sign in to comment.