这是一个叶片前缘求导圆的，

dimension z(100),z0(100),ys(100),yp(100)
dimension z1(100),y1(100),y1_1(100),k(100),x0(100),y0(100),r1(100),r2(100),r12(100)

open(1,file='_1r_r1_100_50.dat')
open(2,file='2.dat')  



write(*,*) '出入截面的坐标点数'
read(*,*) N  
do i=1,N
  read(1,*) z0(i),yp(i),ys(i) !原始坐标
end do
write(*,*) '积叠线位置'
read(*,*) xjd
do i=1,N
  z(i)=z0(i)+xjd
  write(2,*) z(i)
end do


!*****在压力面插值函数曲线前两个节点间取99个等分点*
d1=(z(5)-z(2))/99     !!前缘压力面横坐标间距
k12=(ys(2)-ys(1))/(z(2)-z(1))
k1=-1/k12
a=ys(1)-k1*z(1)
do m=1,100
  z1(m)=z(2)+(m-1)*d1
  y1(m)=ys(3)*(z1(m)-z(4))*(z1(m)-z(5))/((z(3)-z(4))*(z(3)-z(5)))+ys(4)*(z1(m)-z(3))*(z1(m)-z(5))/((z(4)-z(3))*(z(4)-z(5)))+ys(5)*(z1(m)-z(3))*(z1(m)-z(4))/((z(5)-z(3))*(z(5)-z(4)))
  y1_1(m)=ys(3)*(z1(m)-z(4)+z1(m)-z(5))/((z(3)-z(4))*(z(3)-z(5)))+ys(4)*(z1(m)-z(3)+z1(m)-z(5))/((z(4)-z(3))*(z(4)-z(5)))+ys(5)*(z1(m)-z(3)+z1(m)-z(4))/((z(5)-z(3))*(z(5)-z(4)))
  k(m)=-1/y1_1(m)
  b=y1(m)-k(m)*z1(m)
  x0(m)=(a-b)/(k(m)-k1)
  y0(m)=k1*x0(m)+a
  r1(m)=sqrt((x0(m)-z(1))**2+(y0(m)-ys(1))**2)
  r2(m)=sqrt((x0(m)-z1(m))**2+(y0(m)-y1(m))**2)
  r12(m)=abs(r1(m)-r2(m))
end do
dth=0.01
do m=1,100
if(z(1)<x0(m)) then
  if(r12(m)<dth) then
    x0=x0(m)
    y0=y0(m)
    r0=r1(m)
    xt=z1(m)
    yt=y1(m)
  end if
end if
end do
write(*,*)'x0=',x0,'y0=',y0,'r0=',r0 
end
！！！运行后出现sqrt domain error 的错误，找不出错误的地方。
求助！！！！