Il problema sta proprio nella function che restituisce un integer!

Gli integer del fortran sono a 32bit e quindi possono rappresentare piu' o meno un numero nel range da -2 miliardi a +2 miliardi ( precisamente da -2^32 a +2^32 -1).

Poiche' il fattoriale di numeri maggiori di 13 e' piu' grande, da li subentrano i problemi.

La function fattoriale deve usare i DOUBLE PRECISION.