Let's solve the $1 \times n$ case first. We need to find the number of sequences of $n$ terms from $1, 2, \ldots p-1$ whose sum is not divisible by $p$. Let $a_n$ be the number of such sequences. We claim $a_n = (p-2)a_{n-1} + (p-1)a_{n-2}$. Suppose we have an $n$-term sequence. If the sum of the first $n-1$ terms of this sequence is $s \pmod{p}$, then we can pick the first $n-1$ terms in $a_{n-1}$ ways and the last term can be any residue mod $p$ but $0$ and $-s$, for $(p-2)a_{n-1}$ ways. If the sum of the first $n-1$ terms is divisible by $p$, and the $n-1$th term is $t$, then the sum of the first $n-2$ terms is $-t \pmod{p}$ and across all $t$ there are $a_{n-2}$ ways to pick them. Then the last term can be any nonzero residue, for $p-1$ choices. Summing these two cases up gives the desired recurrence. It is trivial to check that $a_1 = p-1$, $a_2 = (p-1)(p-2)$. Combining this with the recurrence, which we can solve with a number of methods (the characteristic polynomial factors with roots $-1$ and $p-1$, so it works out nice), we can figure out that \[a_n = \frac{p-1}{p} \left[ (p-1)^n - (-1)^n \right].\] Now suppose we have an $m \times n$ board. Arbitrarily pick all entries except those in the last row and column. There are $p^{(m-1)(n-1)}$ ways to do this. We will now form a bijection between the number of ways to fill in the remaining $m+n-1$ entries and the number of ways to solve the problem for a $1 \times (m+n-1)$ matrix. Take an arbitrary working $1 \times (m+n-1)$ matrix, and suppose its entries are $x_1, x_2, \ldots x_{m-1}, x^*, x'_1, x'_2, \ldots x'_{n-1}$. Suppose in our partially filled $m \times n$ board that the sum of the $n-1$ picked entries in the $i$th row is $s_i$. Set the last entry of the $i$th row to be whatever is equivalent to $x_i - s_i$ mod $p$. This is guaranteed to make the row not divisible by $p$ because its sum is $x_i$. Now suppose in our partially filled $m \times n$ board that the sum of the $m-1$ picked entries in the $j$th column is $s'_j$. Set the last entry of the $j$th column to be whatever is equivalent to $-x'_j - s'_j$ mod $p$. This is guaranteed to make the column not divisible by $p$ because its sum is $-x'_j$. Finally, set the final entry in the corner to have the residue $x^* + x'_1 + x'_2 + \ldots + x'_{n-1} + s_1 + s_2 + \ldots + s_{m-1}$ mod $p$. This makes the sum of the final row equivalent to $x^*$ mod $p$ and the final column equivalent to $x_1 + x_2 + \ldots + x_{m-1} + x^* + x'_1 + x'_2 + \ldots + x'_{n-1}$ mod $p$, both of which we know are nonzero. (We are using the obvious fact that $s_1 + s_2 + \ldots + s_{m-1} = s'_1 + s'_2 + \ldots + s'_{n-1}$.) It is easy to check that this procedure is invertible. So the total number of ways to do an $m \times n$ board is $p^{(m-1)(n-1)} a_{m+n-1}$, and the answer is \[p^{mn-m-n} (p-1) \left[ (p-1)^{m+n-1} + (-1)^{m+n} \right].\]

We work modulo $p$. Let the row sums be $a_1$, $a_2$, $\ldots$, $a_m$, let the column sums be $b_1$, $b_2$, $\ldots$, $b_n$, and let the sum of all numbers in the grid be $s$, where $s$ and all of the $a_i$ and $b_j$ are remainders modulo $p$. So $\sum a_i \equiv \sum b_j \equiv s \pmod p$. We fill in our grid as follows: First we choose $s$. Then we choose the $a_i$ and $b_j$, subject to the condition that all of them are nonzero and $\sum a_i \equiv \sum b_j \equiv s \pmod p$. We leave the lowermost row and the rightmost column blank, and we fill in the remaining $(m - 1)(n - 1)$ cells arbitrarily. There are $p^{(m - 1)(n - 1)}$ ways to do this. After that, everything else in the grid is determined uniquely because the row and column sums are already fixed modulo $p$. Let $f_m(s)$ denote the number of ways to choose the $a_i$ subject to the condition that all of them are nonzero and $\sum a_i \equiv s \pmod p$. Similarly, let $g_n(s)$ denote the number of ways to choose the $b_j$ subject to the condition that all of them are nonzero and $\sum b_j \equiv s \pmod p$. Then, summing over all possible choices of $s$, we see that the number of valid grids is \[p^{(m - 1)(n - 1)} \sum_s f_m(s)g_n(s).\] We go on to calculate $f_m(s)$. The calculation for $g_n(s)$ will be analogous. Consider the polynomial \[Q(x) = (x + x^2 + \cdots + x^{p - 1})^m.\] The coefficient before $x^k$ in $Q$ is the number of ordered $m$-tuples of elements of the set $\{1, 2, \ldots, p - 1\}$ with sum $k$. Thus $f_m(s)$ is the sum of the coefficients of $Q$ before all $x^k$ with $k \equiv s \pmod p$. To extract this sum, we apply a roots of unity filter. Let \[R(x) = x^{p - s}Q(x).\] Then $f_m(s)$ is the sum of the coefficients of $R$ before all $x^k$ such that $p$ divides $k$. Consequently, \[f_m(s) = \frac{1}{p} \sum_\omega R(\omega),\]where the sum on the right-hand side is taken over all $p$-th complex roots of unity $\omega$. Observe that $R(1) = (p - 1)^m$ and $R(\omega) = (-1)^m\omega^{p - s}$ when $\omega \neq 1$. It follows that \[f_m(s) = \frac{1}{p}\left((p - 1)^m + (-1)^m \sum_{\omega \neq 1} \omega^{p - s}\right).\] When $s = 0$, the above simplifies to \[f_m(0) = \frac{1}{p}\Big((p - 1)^m + (-1)^m(p - 1)\Big).\] Otherwise, when $s \in \{1, 2, \ldots, p - 1\}$, we obtain instead \[f_m(s) = \frac{1}{p}\Big((p - 1)^m + (-1)^{m + 1}\Big).\] (At this point, as a quick reality check we can verify that $\sum_s f_m(s) = (p - 1)^m$.) Of course, analogous formulae hold for $g_n(s)$ as well. Putting it all together, we find that \begin{align*} \sum_s f_m(s)g_n(s) &= \frac{1}{p^2} \bigg( \Big((p - 1)^m + (-1)^m(p - 1)\Big)\Big((p - 1)^n + (-1)^n(p - 1)\Big)\\ &\phantom{= \frac{1}{p^2} \bigg(} + (p - 1) \cdot \Big((p - 1)^m + (-1)^{m + 1}\Big)\Big((p - 1)^n + (-1)^{n + 1}\Big) \bigg)\\ &= \frac{1}{p} \Big((p - 1)^{m + n} + (-1)^{m + n}(p - 1)\Big). \end{align*} Therefore, the total number of valid grids works out to \[p^{mn - m - n} \Big((p - 1)^{m + n} + (-1)^{m + n}(p - 1)\Big),\]and the solution is complete.