[子程序] function可以返回一个数组吗

!
! PURPOSE: Entry point for the console application.
!
!****************************************************************************
module ini
implicit none
real*8,parameter::T=10.,wx=6.,wy=6.,wz=11., pi = 3.141592653589793239
integer,parameter::N=3
!T是温度，单位为微K，N是原子数，6，6，11刻画了原子的构型
contains

function ran() !returns random number between 0 - 1 
implicit none 
integer , save :: flag = 0
double precision :: ran 
if(flag==0) then 
call random_seed()
flag = 1 
endif 
call random_number(ran) ! built in fortran 90 random number function 
end function ran

function normal(mean,sigma) 
implicit none 
integer :: flag 
double precision, parameter :: pi = 3.141592653589793239 
double precision :: u1, u2, y1, y2, normal, mean, sigma 
save flag 
data flag /0/ 
u1 = ran(); u2 = ran() 
if (flag.eq.0) then 
y1 = sqrt(-2.0d0*log(u1))*cos(2.0d0*pi*u2) 
normal = mean + sigma*y1 
flag = 1 
else 
y2 = sqrt(-2.0d0*log(u1))*sin(2.0d0*pi*u2) 
normal = mean + sigma*y2 
flag = 0 
endif 
end function normal 
!The above codes are made in Fortran 90 language, if you have any question, you may write to [url=mailto:sealin2008@hotmail.com]sealin2008@hotmail.com[/url] 

subroutine position0(R0)
implicit none
real*8::R0(N,3),x,y
integer ::i
do i=1,N
do
x=ran()
y=ran()
if((2*x-1)**2+(2*y-1)**2<=1) exit
end do
R0(i,1)=wx*(2.0*x-1.0)
R0(i,2)=wy*(2.0*y-1.0)
end do
do i=1,N
R0(i,3)=normal(0.D0,wz/2.0)
end do
end subroutine position0

subroutine Velocity(V)
implicit none
real*8::V(N,3),r1,r2,r3
integer::i
do i=1,N
r1=ran()*2.*pi
r2=ran()*2.-1.0
r3=normal(0.d0,dsqrt(1.3806504E4*T/(87*1.65053878)))
V(i,1)=r3*sqrt(1-r2*r2)*cos(r1)
V(i,2)=r3*sqrt(1-r2*r2)*sin(r1)
V(i,3)=r3*r2
end do

end subroutine Velocity

function position(i,t,R0,V)
implicit none
real*8::R0(N,3),V(N,3),t,position(3)
integer::i
position(1)=R0(i,1)+t*V(i,1)
position(2)=R0(i,2)+t*V(i,2)
position(3)=R0(i,3)+t*V(i,3)
end function position
end module ini


program main
use ini
implicit none
real*8::R0(N,3),V(N,3)
integer ::i
call position0(R0)
call Velocity(V)
Open( 12 , File = 'out.txt' )
Do i = 1 , N
write( 12 , * ) R0(i,1),R0(i,2),R0(i,3)
End Do
Do i = 1 , N
write( 12 , * ) V(i,1),V(i,2),V(i,3)
End Do
Close( 12 )


end program main

  R0 = position0()
  V = Velocity()

 Function position0() result(R0)
    implicit none
    real*8::R0(N,3),x,y
    integer ::i
    do i=1,N
      do
        x=ran()
        y=ran()
        if((2*x-1)**2+(2*y-1)**2<=1) exit
      end do
      R0(i,1)=wx*(2.0*x-1.0)
      R0(i,2)=wy*(2.0*y-1.0)
    end do
    do i=1,N
      R0(i,3)=normal(0.D0,wz/2.0)
    end do
  end Function position0

  Function Velocity() result(V)
    implicit none
    real*8::V(N,3),r1,r2,r3
    integer::i
    do i=1,N
      r1=ran()*2.*pi
      r2=ran()*2.-1.0
      r3=normal(0.d0,dsqrt(1.3806504E4*T/(87*1.65053878)))
      V(i,1)=r3*sqrt(1-r2*r2)*cos(r1)
      V(i,2)=r3*sqrt(1-r2*r2)*sin(r1)
      V(i,3)=r3*r2
    end do
  end Function Velocity