◆ update_impl() [1/2]

void pism::energy::TemperatureModel::update_impl	(	double	t,
		double	dt,
		const Inputs &	inputs
	)

protectedvirtual

Takes a semi-implicit time-step for the temperature equation.

This method should be kept because it is worth having alternative physics, and so that older results can be reproduced. In particular, this method is documented by papers [BBL,BBssasliding]. The new browser page bombproofenth essentially documents the cold-ice-BOMBPROOF method here, but the newer enthalpy-based method is slightly different and (we hope) a superior implementation of the conservation of energy principle.

The conservation of energy equation written in terms of temperature is

\[ \rho c_p(T) \frac{dT}{dt} = k \frac{\partial^2 T}{\partial z^2} + \Sigma,\]

where \(T(t,x,y,z)\) is the temperature of the ice. This equation is the shallow approximation of the full 3D conservation of energy. Note \(dT/dt\) stands for the material derivative, so advection is included. Here \(\rho\) is the density of ice, \(c_p\) is its specific heat, and \(k\) is its conductivity. Also \(\Sigma\) is the volume strain heating (with SI units of \(J/(\text{s} \text{m}^3) = \text{W}\,\text{m}^{-3}\)).

We handle horizontal advection explicitly by first-order upwinding. We handle vertical advection implicitly by centered differencing when possible, and retreat to implicit first-order upwinding when necessary. There is a CFL condition for the horizontal explicit upwinding [MortonMayers]. We report any CFL violations, but they are designed to not occur.

The vertical conduction term is always handled implicitly (i.e. by backward Euler).

We work from the bottom of the column upward in building the system to solve (in the semi-implicit time-stepping scheme). The excess energy above pressure melting is converted to melt-water, and that a fraction of this melt water is transported to the ice base according to the scheme in excessToFromBasalMeltLayer().

The method uses equally-spaced calculation but the columnSystemCtx methods coarse_to_fine(), fine_to_coarse() interpolate back-and-forth from this equally-spaced computational grid to the (usually) non-equally spaced storage grid.

An instance of tempSystemCtx is used to solve the tridiagonal system set-up here.

In this procedure two scalar fields are modified: basal_melt_rate and m_work. But basal_melt_rate will never need to communicate ghosted values (horizontal stencil neighbors). The ghosted values for m_ice_temperature are updated from the values in m_work in the communication done by energy_step().

The (older) scheme cold-ice-BOMBPROOF, implemented here, is very reliable, but there is still an extreme and rare fjord situation which causes trouble. For example, it occurs in one column of ice in one fjord perhaps only once in a 200ka simulation of the whole sheet, in my (ELB) experience modeling the Greenland ice sheet. It causes the discretized advection bulge to give temperatures below that of the coldest ice anywhere, a continuum impossibility. So as a final protection there is a "bulge limiter" which sets the temperature to the surface temperature of the column minus the bulge maximum (15 K) if it is below that level. The number of times this occurs is reported as a "BPbulge" percentage.

Implements pism::energy::EnergyModel.

Definition at line 168 of file TemperatureModel.cc.