compute_surface_gradient init m_h_x_no_model m_h_y_no_model SIAFD_Regional ~SIAFD_Regional

compute_surface_gradient() void pism::stressbalance::SIAFD_Regional::compute_surface_gradient ( const Inputs &  inputs, array::Staggered1 &  h_x, array::Staggered1 &  h_y  ) privatevirtual 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 from pism::stressbalance::SIAFD. Definition at line 42 of file SIAFD_Regional.cc. References pism::array::Array2D< T >::box(), pism::Geometry::cell_type, pism::stressbalance::SIAFD::compute_surface_gradient(), pism::stressbalance::Inputs::geometry, pism::Component::m_grid, m_h_x_no_model, m_h_y_no_model, pism::stressbalance::Inputs::no_model_mask, pism::stressbalance::Inputs::no_model_surface_elevation, and pism::stressbalance::SIAFD::surface_gradient_haseloff().