In R, what does a probability density function compute? I understand that the probability density function, pdf, of a continuous random variable is the probability of the variable taking on a given value. I also thought that for a continuous random variable, the probability of it taking a specific value is always 0 i.e. P(X=x) = 0.
Therefore, in R, I am wondering what does 

dnorm(0)

compute?
I know it means the probability of obtaining a 0 from a normal distribution (which is one of the distributions a continuous random variable can take) with mean 0 and standard deviation 1. I had expected the output to be 0, since I am asking for the probability of a continuous random variable taking on a specific value. However, the output is 0.3989.
What does 0.3989 represent and why isn't the output 0?
 A: You shoud distinguish between probability mass and probability density.
A discrete random variable is like a set of little stones: each stone has its own mass (weight). For example, if you toss a regular coin the mass of the "head-stone" is 1/2. The probability "density" function of a discrete variable is actually a probability mass function.
A continuous random variable is like a heap of dust. You may think that a single speck of dust has no mass, but you can calculate the mass of a fistful of dust if you know its volume and its density: its mass is $volume\times density$.
dnorm(0) is just the density of a (standard) normal variable. When the volume is zero, i.e. when $Z=0$, the probability mass is zero just because in $volume\times density$ the first factor is zero. As soon as the volume increases, as soon as you look for the probability mass of a range of values (a fistful of dust, i.e. an interval), you can compute a probability mass. For example:
$$\begin{align}P(-0.1<Z<0.1)&=P(Z\in(-0.1,0.1))\\ &=\mathtt{pnorm(0.1)-pnorm(-0.1)}=0.07965567\end{align}$$
Each single value (each speck of dust) has a density, but only intervals like $(-0.1,0.1)$ (fistfuls of dust) have a "mass". So you can compute the probability mass of intervals, not of single values. But the probability density of single values does exist and make sense.
The density of a uniform random variable is constant, the density of a normal variable is not. However you can approximate the mass of small intervals. For example:
> pnorm(0.1)-pnorm(-0.1)
[1] 0.07965567
> density_at0 <- dnorm(0)
> volume_around0 <- (0.1) - (-0.1)
> volume_around0 * density_at0
[1] 0.07978846

A: Densities are not probabilities. 
$$
$$
The density $f$ of a continuous random variable $X$ is defined as
$$
f(x) = F'(x) = \lim_{\epsilon \rightarrow 0} \frac{\Pr\big(x < X \leq x + \epsilon\big)}{\epsilon}
$$
with $F$ the cumulative distribution function (cdf).
$$
$$


*

*$f(x)$ can be larger than $1$. 

*The area under the density equals $1$: $\displaystyle \int_{\text{dom}(f)} f(x) \, \text{d}x = 1$. 

*The area under the density in the interval $(a,b)$ equals the probability that $X \in (a,b)$.


That last property makes the bridge between densities and probabilities.
$$
$$
The density of a $N(\mu, \sigma^2)$ random variable is given by
$$
f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \Big\{ - \frac{1}{2\sigma^2} (x - \mu)^2 \Big\}
$$
In particular, with $\mu=0$ and $\sigma=1$, we have $f(0)= \dfrac{1}{\sqrt{2\pi}} = 0.3989$.

To address @garciaj's comment:
The fact that $f(0) = 0.3989$ indicates that (cf. the definition of $f(x)$ above)
$$
\Pr\big(0 < X \leq \epsilon\big) \approx 0.3989 \, \epsilon
$$
for $\epsilon$ small enough. For example, using $\epsilon = 0.1$, we have
\begin{align*}
\Pr\big(0 < X \leq 0.1\big) 
  & = \Pr\big(X \leq 0.1\big) - \Pr\big(X \leq 0\big) \\
  & = \text{pnorm(0.1)} - \text{pnorm(0)} \\ 
  & = 0.03982784 \approx 0.3989 \times 0.1
\end{align*}

