[数值问题] 求解运行时的问题。daltA1没算进去

program tu5
 implicit none
 real*8 A1,A2,A4
 real*8 Z0,Zr,Zi
 real*8 A,B,C,D
 real*8 E
 real*8 F
 real*8 ref,Rf
 real*8 C1,D1,E1,F1,refs,t
 real,parameter::pi=3.14159
 real*8::A3=0.0
 integer i
 real*8::daltA1=0.0
 integer k
 character(len=80)::filename="data.txt"
 integer,parameter::fileid=10

 !!!!!!!!!!
 A2=0.25
 A4=0.425
 A1=0.625
 Z0=50
 Zr=2
 Zi=-3
 !!!!!!!!!!************************************************
 open(10,file="data.txt")
    do k=1,100
       daltA1=daltA1+0.001
       do i=1,int(0.2/0.001)
          A3=A3+0.001

    call Get_Rf(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,E,F,E1,F1,ref,refs,Rf)
   write(*,*)t
 write(*,"('A3:'F8.5,'Rf:'F8.5,'daltA1:'F8.5)")A3,Rf,daltA1
 write(10,"('A3:'F8.5,'Rf:'F8.5,'daltA1:'F8.5)")A3,Rf,daltA1

 end do
 end do
 close(fileid)
 stop
 end
 !!!!!!!!***************************************************

 !!!!!!!!!!!!!***********************************************求ref（A1）*****************
subroutine A3_to_A(A3,A2,A4,A)
 implicit none
 real*8 A2,A3,A4,A
 real,parameter::pi=3.14159

 A=cos(2.0*pi*A2)*cos(2.0*pi*A4)-sin(2.0*pi*A2)*(-cos(2.0*pi*A4)/tan(2.0*pi*A3)+sin(2.0*pi*A4))              !!A的函数式
!write(*,*)A
 return
 end

 subroutine A3_to_B(A3,A2,A4,Z0,B)
 implicit none
 real*8 A2,A4,B,Z0,A3
 real,parameter::pi=3.14159


 B=Z0*(cos(2.0*pi*A2)*sin(2.0*pi*A4)+sin(2.0*pi*A2)*(sin(2.0*pi*A4)/tan(2.0*pi*A3)+cos(2.0*pi*A4)))            !!B的函数式
!write(*,*)B
 return
 end
 subroutine A3_to_C(A3,A1,A2,A4,Z0,C)
 implicit none
 real*8 A1,A2,A3,A4,C,Z0
 real,parameter::pi=3.14159


 C=1.0/Z0*(cos(2.0*pi*A4)*((-cos(2.0*pi*A2))/tan(2.0*pi*A1)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*A1)+cos(2.0*pi*A2))*((-cos(2.0*pi*A4))/tan(2.0*pi*A3)+sin(2.0*pi*A4)))
 !write(*,*)C
 return
 end

 subroutine A3_to_D(A3,A1,A2,A4,D)
 implicit none
 real*8 A1,A2,A3,A4,D
 real,parameter::pi=3.14159

 D=-sin(2.0*pi*A4)*(-cos(2.0*pi*A2)/tan(2.0*pi*A1)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*A1)+cos(2.0*pi*A2))*(sin(2.0*pi*A4)/tan(2.0*pi*A3)+cos(2.0*pi*A4))
 !write(*,*)D
 return
 end
 !!!!!!!!!!!***************************************E,F,ref*******
 subroutine Get_E(Z0,Zi,Zr,A3,A1,A2,A4,E)
 implicit none
 real*8 A,B,C,D,Z0,Zi,Zr,A3,E,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call A3_to_C(A3,A1,A2,A4,Z0,C)
 call A3_to_D(A3,A1,A2,A4,D)


 E=(A*Zr*(D-C*Zi)+C*Zr*(B+A*Zi))/(((D-C*Zi)**2+Zr**2*C**2)*Z0)
 !write(*,*)E
 return
 end

 subroutine Get_F(A3,Z0,Zi,Zr,A1,A2,A4,F)
 implicit none
 real*8 A,B,C,D,Zi,Zr,Z0,F,A3,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call A3_to_C(A3,A1,A2,A4,Z0,C)
 call A3_to_D(A3,A1,A2,A4,D)

 F=((B+A*Zi)*(D-C*Zi)-A*C*Zr**2)/(((D-C*Zi)**2+Zr**2*C**2)*Z0)
 !write(*,*)F
 return
 end

 subroutine Get_ref(A3,A1,A2,A4,Z0,Zi,Zr,E,F,ref)
 implicit none
 real*8 E,F,A3,ref,Z0,Zi,Zr,A1,A2,A4
 call Get_E(Z0,Zi,Zr,A3,A1,A2,A4,E)
 call Get_F(A3,Z0,Zi,Zr,A1,A2,A4,F)

 ref=((E**2-1+F**2)**2+4*F**2)/(((E+1)**2+F**2)**2)
 !write(*,*)ref
 return
 end
 !!!!!!!!!!!!!!!!!!**********************************************************
 !!!!!!!!!!!!!!!!!!!*******************求ref（A1+daltA1）********************
subroutine daltA1_t0_t(daltA1,A1,t)
 implicit none
 real*8 daltA1,A1,t

 t=daltA1+A1

 !write(*,*)t
 return
 end

 subroutine Get_C1(A3,daltA1,A1,A2,A4,Z0,t,C1)
 implicit none
 real*8 A1,A2,A3,A4,C1,Z0,daltA1,t
 real,parameter::pi=3.14159

 call daltA1_t0_t(daltA1,A1,t)

 C1=1.0/Z0*(cos(2.0*pi*A4)*((-cos(2.0*pi*A2))/tan(2.0*pi*t)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*t)+cos(2.0*pi*A2))*((-cos(2.0*pi*A4))/tan(2.0*pi*A3)+sin(2.0*pi*A4)))
 !write(*,*)C1
 return
 end

 subroutine Get_D1(A3,daltA1,A1,A2,A4,t,D1)
 implicit none
 real*8 A1,A2,A3,A4,daltA1,t,D1
 real,parameter::pi=3.14159

 call daltA1_t0_t(daltA1,A1,t)

 D1=-sin(2.0*pi*A4)*(-cos(2.0*pi*A2)/tan(2.0*pi*t)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*t)+cos(2.0*pi*A2))*(sin(2.0*pi*A4)/tan(2.0*pi*A3)+cos(2.0*pi*A4))
 !write(*,*)D1
 return
 end
 !!!!!!!!!!******************E1,F1,refs**********************
 subroutine Get_E1(A3,daltA1,A1,A2,A4,Z0,Zi,Zr,t,E1)
 implicit none
 real*8 A,B,C1,D1,Z0,Zi,Zr,A3,E1,daltA1,t,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call Get_C1(A3,daltA1,A1,A2,A4,Z0,t,C1)
 call Get_D1(A3,daltA1,A1,A2,A4,t,D1)


 E1=(A*Zr*(D1-C1*Zi)+C1*Zr*(B+A*Zi))/(((D1-C1*Zi)**2+Zr**2*C1**2)*Z0)

 return
 end
 subroutine Get_F1(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,F1)
 implicit none
 real*8 A,B,C1,D1,Zi,Zr,Z0,F1,A3,daltA1,t,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call Get_C1(A3,daltA1,A1,A2,A4,Z0,t,C1)
 call Get_D1(A3,daltA1,A1,A2,A4,t,D1)


 F1=((B+A*Zi)*(D1-C1*Zi)-A*C1*Zr**2)/(((D1-C1*Zi)**2+Zr**2*C1**2)*Z0)

 return
 end

 subroutine Get_refs(A3,A1,A2,A4,Zi,Zr,Z0,E1,F1,t,refs)
 implicit none
 real*8 E1,F1,A3,refs,daltA1,A1,A2,A4,Z0,Zi,Zr,t
 call Get_E1(A3,daltA1,A1,A2,A4,Z0,Zi,Zr,t,E1)
 call Get_F1(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,F1)

 refs=((E1**2-1+F1**2)**2+4*F1**2)/(((E1+1)**2+F1**2)**2)
 !write(*,*)refs
 return
 end
 !!!!!!!!!!!!!!!!!!!!!!!!*********************************************
 subroutine Get_Rf(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,E,F,E1,F1,ref,refs,Rf)
 implicit none
 real*8 ref,refs,Rf
 real*8 A3,daltA1,E,F,E1,F1,A1,A2,A4,Z0,Zi,Zr,t
 call Get_ref(A3,A1,A2,A4,Z0,Zi,Zr,E,F,ref)
 call Get_refs(A3,A1,A2,A4,Zi,Zr,Z0,E1,F1,t,refs)

 Rf=refs/ref

 return
 end