[讨论] 计算方法的作业——请网友点评一下

!/*************************************************************************
!   > File Name: kt0301_2.f90
!   > Author: 青羽100P 
!   > Mail: [url=mailto:3443988515@qq.com]3443988515@qq.com[/url] 
!   > Created Time: 2017年10月11日 星期三 18时03分27秒
! ************************************************************************/

      program kt0301_2
!拉格朗日主程序代码
         implicit none
         integer::pa(17),pa_1(4),i,z,error
         real::tem(17),tem_1(4),chads
         !使用数组存放数值
         open(10,file='sample_unformatted_sequential.dat',f&
             &orm='unformatted',access='sequential',status='old')
         i=1
         do
           read(10,iostat=error),tem(i),pa(i)
           if(error/=0)exit
           i=i+1
         end do
         close(10)
         pa(1)=1000
         pa_1=(/pa(1),pa(4),pa(7),pa(10)/)
         tem_1=(/tem(1),tem(4),tem(7),tem(10)/)
         !挑出1000，700，400，200处的数值作为已知点
         !-------------------------
         open(11,file="chazhishuju.txt")
         do z=1,17
             call lar(4,pa_1,tem_1,pa(z),chads)
             write(11,"(i2,f10.4,i5,f10.4,f10.4)"),z,tem(z),pa(z),chads,chads-tem(z)
         end do
      end program
      include 'lar.f90'

!/************************************************************************
 !  > File Name: lar.f90
 !  > Author: 青羽100P 
 !  > Mail: [url=mailto:3443988515@qq.com]3443988515@qq.com[/url] 
 !  > Created Time: 2017年10月10日 星期二 21时59分40秒
 ! *************拉格朗日子程序代码***********************************************************/

      subroutine lar(n,x,y,cha,ds)
         implicit none
          integer::n,x(n),cha,i,j
          !n表示数组长度，x横坐标，cha表示要差值x的位置
          real::y(n),larlc,ds
          ds=0
          larlc=1
          do j=1,n
          larlc=1
          !保证不受上一次运算结果影响
              do i=1,n
                 if(i==j)then
                     larlc=larlc*1.0
                 else
                     larlc=larlc*((cha-x(i)*1.0)/(x(j)*1.0-x(i)*1.0))
                 end if
                 !此处if 块，实现连乘算式，当算到自己是等于1就可以了，也可以直接用cycle
                     !每个算式后 *1.0 保证实型运算，否则会被省略，结果全为零 
              end do
              ds=ds+larlc*y(j)
          end do
          print*,ds
      end subroutine