[求助] 新手入门，求助一个求解费米能级的问题

PROGRAM TEST
DOUBLE PRECISION DM, DN, D, EB, T, F
T=300.0D0
DM=0.067D0
DN=1.0D14
EB=8.0D-3
D=1.0D17
F=-0.5D0
T=T/1.1604D4
CALL FERMI(DN,DM,D, EB,T,F)
PRINT *, DN, DM,T,F
STOP
END
SUBROUTINE FERMI(DN,DM, D, EB,T,F)
DOUBLE PRECISION DM, DN, D, EB, F, T,S,X,Y,Z
X=DN/(1.0D0+DEXP(EB/T))/((T*DM)**1.5D0*6.8122D21)
IF(DN.GE.D) X=DN/((T*DM)**1.5D0*6.8122D21)
S=1.0D-1
F=F/T
DO 100 I=1,10
S=S*0.5D0
DO 110 J=1,5000
F=F+S
Y=0.0D0
IF (-F.GT.5.0D2) GOTO 110
Y=F**4+33.6*F*(1.0-0.68D0*DEXP(-0.17*(1.0+F)**2))
Y=1.0/(DEXP(-F)+1.32934038675/(50.0+Y)**0.375)-X
IF (I+J.EQ.2) Z=Y
110 IF (Y*Z.LE.0.0) GOTO 100
100 F=F-S
F=(F+0.5D0*S)*T
RETURN
END