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

You can look at the approach used there and port it to Fortran.

NOTE: To the best of my knowledge there is no known algorithm for finding the optimal horner factorisation of a multivariate polynomial. The best one can do is to use a clever heuristic to find "good" factorisations. I implemented a greedy heuristic similar to the one described in the paper "[Greedy Algorithms for Optimizing Multivariate Horner Schemes][2]".

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]: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.330.7430&rep=rep1&type=pdf