From 14f5f8df2d224ba9cb3f6a8bec9a3eaaf7cd0610 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Fri, 11 May 2018 22:35:38 -0600 Subject: [PATCH] Reverted to previous hybrid gradient scheme In a recent commit (89f68751), 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. --- libglissade/glissade_grid_operators.F90 | 20 ++----- libglissade/glissade_velo_higher.F90 | 70 ++++++++----------------- 2 files changed, 25 insertions(+), 65 deletions(-) diff --git a/libglissade/glissade_grid_operators.F90 b/libglissade/glissade_grid_operators.F90 index eedc0137..ddfcf24a 100644 --- a/libglissade/glissade_grid_operators.F90 +++ b/libglissade/glissade_grid_operators.F90 @@ -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. !---------------------------------------------------------------- !---------------------------------------------------------------- @@ -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. @@ -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 @@ -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 diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90 index 0e5f0279..b0d8b00b 100644 --- a/libglissade/glissade_velo_higher.F90 +++ b/libglissade/glissade_velo_higher.F90 @@ -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 @@ -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 @@ -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 !------------------------------------------------------------------------------ @@ -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 @@ -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') @@ -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) :: & @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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) @@ -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.