I implemented this in Python: multivar_horner
You can look at the approach used there and port it to Fortran.
NOTE: In contrast to the one dimensional case there are multiple possible Horner factorisations of multivariate polynomials. One can allow a search over the possible factorisations to find a minimal representation as described HERE. In most cases however it suffices to use a heuristic to find a single "good" factorisation. multivar_horner implements the greedy heuristic described in "Greedy Algorithms for Optimizing Multivariate Horner Schemes".
The authors of some related publications (including the above mentioned) claim to have an implementation of their proposed algorithms, but I was not able to find any publicly available ones.