I implemented this in Python: multivar_horner
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".
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.