[数值问题] 我这个小程序哪儿有问题啊，谢谢

program tyt20
implicit none
integer i,j,tmp
real*8::Z,Zi,u1,u2,Zf1,Zf2,R
real,parameter::pi=3.14159
real*8 A,B,C,D
real*8 A2,A4,A01,A02,b1,b2,n1,n2
real*8 E
real*8 F
real*8 H1,H2
real*8 t         !!反射系数的平方用t表示
character(len=80)::filename="data.dat"
integer,parameter::fileid=10 
open(10,file='data.dat')

R=2
H1=80
do i=1,int(40/0.01)
    H1=H1+0.01
    !n2=n2+0.001
    !R=R+1   
H2=80
do j=1,int(40/0.01)
     H2=H2+0.01
u1=1+1200.0/H1
u2=1+1200.0/H2
  n1=0.45
  n2=0.3
  A01=0.5
  A02=0.3
  Zi=-3                 
  Z=50.0                 
  A2=0.375    
  A4=0.35       
  Zf1=sqrt(u1/14.0)  
  Zf2=sqrt(u2/14.0)  
b1=(Zf1*tan(sqrt(u1*14.0)*n1*2.0*pi*A01)+Z*tan((1-n1)*2.0*pi*A01))/(Z-Zf1*tan(sqrt(u1*14.0)*n1*2.0*pi*A01)*tan((1-n1)*2.0*pi*A01))
b2=(Zf2*tan(sqrt(u2*14.0)*n2*2.0*pi*A02)+Z*tan((1-n2)*2.0*pi*A02))/(Z-Zf2*tan(sqrt(u2*14.0)*n2*2.0*pi*A02)*tan((1-n2)*2.0*pi*A02))
A=(cos(2*pi*A2)+sin(2*pi*A2)/b2)*cos(2*pi*A4)-sin(2*pi*A2)*sin(2*pi*A4)
B=((cos(2*pi*A2)+sin(2*pi*A2)/b2)*sin(2*pi*A4)+sin(2*pi*A2)*cos(2*pi*A4))*Z
C=(((-1/b1)*cos(2*pi*A2)+sin(2*pi*A2)-((1/b1)*sin(2*pi*A2)+cos(2*pi*A2))/b2)*cos(2*pi*A4)+((1/b1)*sin(2*pi*A2)+cos(2*pi*A2))*sin(2*pi*A4))/Z
D=((1/b1)*cos(2*pi*A2)-sin(2*pi*A2)+((1/b1)*sin(2*pi*A2)+cos(2*pi*A2))/b2)*sin(2*pi*A4)+((1/b1)*sin(2*pi*A2)+cos(2*pi*A2))*cos(2*pi*A4)
E=(A*R*(D-C*Zi)+C*R*(B+A*Zi))/(((D-C*Zi)**2+(R*C)**2)*Z)
F=((B+A*Zi)*(D-C*Zi)-(A*C)*(R**2))/(((D-C*Zi)**2+(R*C)**2)*Z)
t=100*((E**2-1+F**2)**2+4*(F**2))/(((E+1)**2+F**2)**2)
  tmp=t
if (t(i,j)<tmp)then
    tmp=t(i,j)
end if
end do
write(*,"('H1:'F10.4,'H2:'F10.4,'t:'F11.4)")H1,H2,t
write(10,"(F10.4,' 'F10.4,''F11.4)")H1,H2,t
end do
close(10)
stop
end