bed_smoother compute_3d_horizontal_velocity compute_diffusive_flux compute_diffusivity compute_I compute_surface_gradient diagnostics_impl diffusivity init interglacial m_bed_smoother m_D m_delta_0 m_delta_1 m_e_factor m_e_factor_interglacial m_eemian_end m_eemian_start m_event_sia m_h_x m_h_y m_holocene_start m_seconds_per_year m_stencil_width m_work_2d_0 m_work_2d_1 m_work_3d_0 m_work_3d_1 SIAFD surface_gradient_eta surface_gradient_haseloff surface_gradient_mahaffy surface_gradient_x surface_gradient_y update ~SIAFD

surface_gradient_haseloff() void pism::stressbalance::SIAFD::surface_gradient_haseloff ( const array::Scalar2 &  ice_surface_elevation, const array::CellType2 &  cell_type, array::Staggered1 &  h_x, array::Staggered1 &  h_y  ) protectedvirtual Compute the ice surface gradient using a modification of Marianne Haseloff's approach. The original code deals correctly with adjacent ice-free points with bed elevations which are above the surface of the ice nearby. This is done by setting surface gradient at the margin to zero at such locations. This code also deals with shelf fronts: sharp surface elevation change at the ice shelf front would otherwise cause abnormally high diffusivity values, which forces PISM to take shorter time-steps than necessary. (Note that the mass continuity code does not use SIA fluxes in floating areas.) This is done by assuming that the ice surface near shelf fronts is horizontal (i.e. here the surface gradient is set to zero also). The code below uses an interpretation of the standard Mahaffy scheme. We compute components of the surface gradient at staggered grid locations. The field h_x stores the x-component on the i-offset and j-offset grids, h_y — the y-component. The Mahaffy scheme for the x-component at grid points on the i-offset grid (offset in the x-direction) is just the centered finite difference using adjacent regular-grid points. (Similarly for the y-component at j-offset locations.) Mahaffy's prescription for computing the y-component on the i-offset can be interpreted as: compute the y-component at four surrounding j-offset staggered grid locations, compute the average of these four. The code below does just that. The first double for-loop computes x-components at i-offset locations and y-components at j-offset locations. Each computed number is assigned a weight (w_i and w_j) that is used below The second double for-loop computes x-components at j-offset locations and y-components at i-offset locations as averages of quantities computed earlier. The weight are used to keep track of the number of values used in the averaging process. This method communicates ghost values of h_x and h_y. They cannot be computed locally because the first loop uses width=2 stencil of surface, mask, and bed to compute values at all grid points including width=1 ghosts, then the second loop uses width=1 stencil to compute local values. (In other words, a purely local computation would require width=3 stencil of surface, mask, and bed fields.) Definition at line 360 of file SIAFD.cc. References pism::Component::m_grid, m_work_2d_0, m_work_2d_1, pism::array::Array::stencil_width(), and pism::array::Array::update_ghosts(). Referenced by pism::stressbalance::SIAFD_Regional::compute_surface_gradient(), and compute_surface_gradient().