Your equation with $\nabla \mathbf{P}$ makes sense if this is the
divergence $\mathbf{\nabla}\cdot \mathbf{P}$. That is, the Poynting
vector is the energy flux, i.e. the energy per unit are per unit time.
To calculate the absorption, you want the net electromagnetic energy
entering the volume. This is the surface integral of the Poynting vector,
with a negative sign since you want the energy entering, and the
area element points outward. Writing the surface integral and using
the divergence theorem:
\begin{equation}
A = -\int_S \mathbf{P}\cdot d\mathbf{S} = -\int_V \mathbf{\nabla}\cdot
\mathbf{P} dV
\end{equation}
Taking the volume to zero gives your $A/V$ equation.
For single frequency harmonic fields, the time average of $\mathbf{P}$
is $\mathbf{S}$, so the time averaged absorption per unit volume is
\begin{equation}
\langle A/V \rangle = \mathbf{\nabla}\cdot \mathbf{S} \,.
\end{equation}
If you have only $\mathbf{S}(x,y,z,w)$, then you can approximate the
divergence as a central difference
\begin{eqnarray}
\mathbf{\nabla} \cdot \mathbf{S}(x,y,z,w) \simeq
\frac{\mathbf{S}(x+\Delta x,y,z,w)
-\mathbf{S}(x-\Delta x,y,z,w)}{2\Delta x}
\nonumber\\
+\frac{\mathbf{S}(x,y+\Delta y,z,w)
-\mathbf{S}(x,y-\Delta y,z,w)}{2\Delta x}
\nonumber\\
+\frac{\mathbf{S}(x,y,z+\Delta z,w)
-\mathbf{S}(x,y,z-\Delta z,w)}{2\Delta x}
\end{eqnarray}
with whatever conveniently small $\Delta x$ your simulation gives you.
Added:
The $\frac{1}{2}$ and real part are the time average. That is
with $A(t) = {\rm Re} A e^{-i\omega t}$, $B(t)={\rm Re} B e^{-i\omega t}$ where the real part of $A$ is half the sum of $A$ and its
complex conjugate $A^*$. The product is
\begin{equation}
A(t)B(t) = \frac{1}{4} \left [ AB^* +A^*B + ABe^{-i2\omega t}+A^*B^*e^{i2\omega t} \right ]
\end{equation}
Time averaging gives just the first term
\begin{equation}
\langle A(t)B(t)\rangle = \frac{1}{4} \left [ AB^* +A^*B \right ] = \frac{1}{2} {\rm Re} AB^*
\end{equation}
More additions:
Most people do not write the explicit time
dependence when they write single frequency results. Here is the above time average result applied to the Poynting vector explicitly.
That is, your time varying electric and magnetic fields for an angular
frequency $\omega$ are
\begin{eqnarray}
\mathbf{E}(x,y,z,t) &=& {\rm Re} \left [
\mathbf{E}(x,y,z,\omega) e^{-i\omega t} \right ]
=\frac{1}{2}\left [\mathbf{E}(x,y,z,\omega)e^{-i\omega t}
+\mathbf{E^*}(x,y,z,\omega)e^{i\omega t} \right]
\nonumber\\
\mathbf{H}(x,y,z,t) &=& {\rm Re} \left[
\mathbf{H}(x,y,z,\omega) e^{-i\omega t} \right]
=\frac{1}{2}\left [\mathbf{H}(x,y,z,\omega)e^{-i\omega t}
+\mathbf{H^*}(x,y,z,\omega)e^{i\omega t} \right]
\nonumber\\
\end{eqnarray}
So the correct time varying Poynting vector is
\begin{eqnarray}
\mathbf{P}(x,y,z,t) &=& \mathbf{E}(x,y,z,t) \times \mathbf{H}(x,y,z,t)
\nonumber\\
&=&\frac{1}{4} \left [
\mathbf{E}(x,y,z,\omega)\times \mathbf{H^*}(x,y,z,\omega)
+\mathbf{E^*}(x,y,z,\omega)\times \mathbf{H}(x,y,z,\omega) \right]
\nonumber\\
&&
+\frac{1}{4} \left [
\mathbf{E}(x,y,z,\omega)\times\mathbf{H}(x,y,z,\omega) e^{-i2\omega t}
+\mathbf{E^*}(x,y,z,\omega)\times\mathbf{H^*}(x,y,z,\omega)
e^{i2\omega t}\right]
\nonumber\\
\end{eqnarray}
The time average of $e^{\pm i 2\omega t}$ is zero. If that is not obvious, you can just do the time average
\begin{equation}
\left \langle e^{\pm 2i\omega t}\right\rangle = \frac{\omega}{2\pi}
\int_0^{\frac{2\pi}{\omega}} dt e^{\pm i2\omega t} = \int_0^{2\pi} \frac{du}{2\pi} e^{\pm i2u} = 0
\end{equation}
The first
term of $\mathbf{P}(x,y,z,t)$ is independent of time so it equals its time average.
The time averaged Poynting vector is
\begin{eqnarray}
\mathbf{S}
&=&\frac{1}{4} \left [
\mathbf{E}(x,y,z,\omega)\times \mathbf{H^*}(x,y,z,\omega)
+\mathbf{E^*}(x,y,z,\omega)\times \mathbf{H}(x,y,z,\omega) \right]
\nonumber\\
&=& \frac{1}{2}
{\rm Re} \left [
\mathbf{E}(x,y,z,\omega)\times \mathbf{H^*}(x,y,z,\omega)\right]
\nonumber\\
&=& \frac{1}{2}
{\rm Re} \left [\mathbf{E^*}(x,y,z,\omega)\times \mathbf{H}(x,y,z,\omega)
\right ]
\end{eqnarray}
