最近师弟我在算rotation matrix elemetns相关的，其中有一条是closure的check。具体数学公式之类不详细解释，换成通俗语言问题是这样的:

文件A: 有160个数据，分别是x1(i),y1(i),z1(i)

文件B: 有160个数据，分别是x2(i),y2(i),z2(i)

经过简单运算，二文件数据之差是 x(i),y(i),z(i)

现有S=3(COS(X)^2-1)/2的计算公式，其中COS(X)^2是由x(i),y(i),z(i)运算得来。

现在写的程序是这样的:

      program main
	implicit none

	integer i,j
	integer ios
	integer m,n
	character(160) :: filename1,filename2,filename3,filename4,tmp 
	integer id
	integer n1,n2
	real x1,y1,z1
	real x2,y2,z2
	real x(10000),y(10000),z(10000)
	real v1,v2,v3
	real a,b,c
	real up,down,coss(10000),s1(10000),s2(10000),s3(10000),avg
      integer GetFileN,ii,jj
	character(len=10)::status


	print*, 'Running...'
      filename1='c6_'
      filename2='c8_'
      do m=1,1   ! number of files to read
      write(tmp,*)m    	 	
      open(10,file=trim(filename1)//trim(adjustl(tmp))
	1//'.dat',status='old' )
      open(20,file=trim(filename2)//trim(adjustl(tmp))
	1//'.dat',status='old' )
      do i=1,GetFileN(10) 
	read(10,*,iostat=ios) n1,x1,y1,z1   !read c6
	read(20,*,iostat=ios) n2,x2,y2,z2   !read c8
	if (ios /=0) then
	exit
	endif
      x(i)=x2-x1    !c8-c6 x
	y(i)=y2-y1    !c8-c6 y
	z(i)=z2-z1    !c8-c6 z
c	print*,x(i),y(i),z(i)
c      write(30,'(3f12.5)')x(i),y(i),z(i)
	a=0
	b=0
	c=1 ! labframe
	up=(x(i)*a+y(i)*b+z(i)*c)*(x(i)*a+y(i)*b+z(i)*c)
	down=(x(i)*x(i)+y(i)*y(i)+z(i)*z(i))*(a*a+b*b+c*c)
	coss(i)=up/down   ! cos(i)^2
c      print*, coss(i)   
	s1(i)=(3*coss(i)-1)/2  !n=0 d00
	s2(i)=-sqrt(1.5)*sqrt(1-coss(i))*sqrt(coss(i)) ! n=1 d10
	s3(i)=sqrt(0.375)*(1-coss(i)) !n=2 d20
c	print*, s1(i) 
c     print*, s2(i)
c	print*, s3(i)
	enddo
	enddo
	print*, 'Done.'
      close(10)
      close(20)
	end

      integer function GetFileN(iFileUnit)
      implicit none
      logical , parameter :: b = .True.
      integer , intent( IN ) :: iFileUnit
      character*(1) :: c
      GetFileN = 0
      rewind( iFileUnit )
	do while (b)
      read( iFileUnit , * ,end =999 ,Err = 999 )c
      GetFileN = GetFileN + 1
      end Do
999   rewind( iFileUnit )
      return
      end function GetFileN


请忽略n1和n2，不需要他们参与计算。并且对应A和B两种文件其实每种有上万个(所以用了
trim(filename1)//trim(adjustl(tmp))
	1//'.dat')

因为最后需要对这个求出的s1，s2,s3值求平均，用现在的写法计算起来超!慢。我想请教各位师兄师姐有没有别的办法可以更快的计算，谢谢您的指教!