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

compute_surface_gradient() void pism::stressbalance::SIAFD::compute_surface_gradient ( const Inputs &  inputs, array::Staggered1 &  h_x, array::Staggered1 &  h_y  ) protectedvirtual Compute the ice surface gradient for the SIA. There are three methods for computing the surface gradient. Which method is controlled by configuration parameter sia.surface_gradient_method which can have values haseloff, mahaffy, or eta. The most traditional method is to directly differentiate the surface elevation $h$ by the Mahaffy method [Mahaffy]. The haseloff method, suggested by Marianne Haseloff, modifies the Mahaffy method only where ice-free adjacent bedrock points are above the ice surface, and in those cases the returned gradient component is zero. The alternative method, when sia.surface_gradient_method = eta, transforms the thickness to something more regular and differentiates that. We get back to the gradient of the surface by applying the chain rule. In particular, as shown in [Vazquezetal2003] for the flat bed and $n=3$ case, if we define $$\eta = H^{(2n+2)/n}$$ then $\eta$ is more regular near the margin than $H$. So we compute the surface gradient by $$\nabla h = \frac{n}{(2n+2)} \eta^{(-n-2)/(2n+2)} \nabla \eta + \nabla b,$$ recalling that $h = H + b$. This method is only applied when $\eta > 0$ at a given point; otherwise $\nabla h = \nabla b$. In all cases we are computing the gradient by finite differences onto a staggered grid. In the method with $\eta$ we apply centered differences using (roughly) the same method for $\eta$ and $b$ that applies directly to the surface elevation $h$ in the mahaffy and haseloff methods. Parameters [out] h_x the X-component of the surface gradient, on the staggered grid [out] h_y the Y-component of the surface gradient, on the staggered grid Reimplemented in pism::stressbalance::SIAFD_Regional. Definition at line 188 of file SIAFD.cc. References pism::Geometry::bed_elevation, pism::Geometry::cell_type, pism::RuntimeError::formatted(), pism::stressbalance::Inputs::geometry, pism::Geometry::ice_surface_elevation, pism::Geometry::ice_thickness, pism::Component::m_config, PISM_ERROR_LOCATION, surface_gradient_eta(), surface_gradient_haseloff(), and surface_gradient_mahaffy(). Referenced by pism::stressbalance::SIAFD_Regional::compute_surface_gradient(), and update().