public subroutine ply_sampling_var_compute_elemdev(var, threshold, min_mean)

This routine computes for each element whether the solution in it is
considered to be deviating from the mean above the given threshold or
not. The logical result is stored in var%deviates for each element.

The total number of deviating elements is stored in var%nDeviating.

The variation is computed by the sum of the absolute values of all higher
modes divided by the first mode.
As we are using series of Legendre polynomials this also is a
bounding estimation for the maximal (relative) deviation from the mean in
this element.

A variation of 0.01 for example would imply that the state in the element
is guaranteed to nowhere deviate from the mean by more than 1 percent.
However, this is a very rough estimation and the actual maximal deviation
from the mean is probably much lower (at least for sufficiently high
polynomial degrees).

If the mean is too close to 0, we use epsilon for the normalization
instead of the actual mean.

The computation is done for the current data found in var%dat, any
previous computations of these flags will be discarded by this routine.

Nodes of different colours represent the following:

Solid arrows point from a procedure to one which it calls. Dashed
arrows point from an interface to procedures which implement that interface.
This could include the module procedures in a generic interface or the
implementation in a submodule of an interface in a parent module.