How do I calculate the delta term of a Convolutional Layer, given the delta terms and weights of the previous Convolutional Layer?
I am first deriving the error for a convolutional layer below for simplicity for a one dimensional array (input) which can easily be transferred to a multidimensional then:
We assume here that the $y^{l-1}$ of length $N$ are the inputs of the $l-1$-th conv. layer, $m$ is the kernel-size of weights $w$ denoting each weight by $w_i$ and the output is $x^l$.
Hence we can write (note the summation from zero): $$x_i^l = \sum\limits_{a=0}^{m-1} w_a y_{a+i}^{l-1}$$
where $y_i^l = f(x_i^l)$ and $f$ the activation function (e.g. sigmoidal).
With this at hand we can now consider some error function $E$ and the error function at the convolutional layer (the one of your previous layer) given by $\partial E / \partial y_i^l $. We now want to find out the dependency of the error in one the weights in the previous layer(s):
\begin{equation} \frac{\partial E}{\partial w_a} = \sum\limits_{a=0}^{N-m} \frac{\partial E}{\partial x_i^l} \frac{\partial x_i^l}{\partial w_a} = \sum\limits_{a=0}^{N-m}\frac{\partial E}{\partial w_a} y_{i+a}^{l-1} \end{equation}
where we have the sum over all expression in which $w_a$ occurs, which are $N-m$. Note also that we know the last term arises from the fact that $\frac{\partial x_i^l}{\partial w_a}= y_{i+a}^{l-1}$ which you can see from the first equation.
To compute the gradient we need to know the first term, which can be calculated by:
$$ \frac{\partial E}{\partial x_i^l} = \frac{\partial E}{\partial y_i^l} \frac{\partial y_i^l}{\partial x_i^l} = \frac{\partial E}{\partial y_i^l} \frac{\partial}{\partial x_i^l} f(x_i^{l})$$
where again the first term is the error in the previous layer and $f$ the nonlinear activation function.
Having all necessary entities we are now able to calculate the error and propagate it back efficiently to the precious layer:
$$ \delta^{l-1}_a = \frac{\partial E}{\partial y_i^{l-1} } = \sum\limits_{a=0}^{m-1} \frac{\partial E}{\partial x_{i-a}^l} \frac{\partial x_{i-a}^l}{\partial y_i^{l-1}} = \sum\limits_{a=0}^{m-1} \frac{\partial E}{\partial x^l_{i-a}} w_a^{flipped}$$
Note that the last step can be understood easy when writing down the $x_i^l$-s w.r.t. the $y_i^{l-1}$-s. The $flipped$ refers to a transposed weight maxtrix ($T$).
Therefore you can just calculate the error in the next layer by (now in vector notation):
$$\delta^{l} = (w^{l})^{T} \delta^{l+1} f'(x^{l})$$
which becomes for a convolutional and subsampling layer:
$$\delta^{l} = upsample((w^{l})^{T} \delta^{l+1}) f'(x^{l})$$
where the $upsample$ operation propagates the error through the max pooling layer.
Please feel free to add or correct me!
For references see:

Convolutional Neural Network by Stanford
Convolutional Neural Network by Andrew Gibiansky

and for a C++ implementation this (without requirement to install)