I implemented this in Python: [multivar_horner][1]

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][2]. In most cases however it suffices to use a heuristic to find a single "good" factorisation. [multivar_horner][1] implements the greedy heuristic described in "[Greedy Algorithms for Optimizing Multivariate Horner Schemes][3]".

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.

  [1]: https://github.com/MrMinimal64/multivar_horner
  [2]: https://multivar-horner.readthedocs.io/en/latest/3_optimal_factorisation.html
  [3]: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.330.7430&rep=rep1&type=pdf