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.