FORTRAN:虚拟参数

FORTRAN: dummy argument

我正在尝试使用RK-4阶显式方法来解决(平面)三体问题。我有一组 4*N 一阶 ODE,其中 N 是对象的数量。 4 是因为对于每个对象,我的位置向量等于 {x,y},速度向量等于 {vx,vy}。

不幸的是程序不能正常运行,它在派生子程序中给我一个警告,我无法理解:

A dummy argument with an explicit INTENT(OUT) declaration is not given an explicit value

我试图从子例程中删除所有 INTENT 声明,以尝试解决代码中的问题。在那之后没有出现警告,但程序仍然继续错误地工作。 这是代码:

MODULE constant
!----------------------------
!   Initial Condition
!----------------------------
REAL*8, PARAMETER :: m0 = 1.9891e30 ! mass sun (kg)
REAL*8, PARAMETER :: m1 = 5.9736e24 ! mass earth (kg)
REAL*8, PARAMETER :: m2 = 1e20 ! mass of test particle (kg)
REAL*8, PARAMETER :: G = 6.6738e-11 ! G (Nm**2/kg**2)

END MODULE constant


    PROGRAM bodies_rk
    use constant
    implicit none
    integer :: n_step,n_periods,k_max,k
    real*8 :: ome,dt,duepi,m_per,t_per
    real*8 :: x0,y0,z0,vx0,vy0,vz0
    real*8 :: x1,y1,z1,vx1,vy1,vz1
    real*8 :: x2,y2,z2,vx2,vy2,vz2
    real*8 :: a_step,h,t,phi,covx,covy,covz,comx,comy,comz
    real*8 :: r0,r1,r2
    real*8, DIMENSION(12):: x

!--------------------------------
!       Test run for Runge-Kutta
!--------------------------------
duepi=8.d0*DATAN(1.d0)
ome=1.d0
t_per=duepi/ome
n_step=1000                    ! step per period
n_periods=5000                 ! periods

!! position of the Sun !! (in m)
x0 = 0.d0 ; y0 = 0.d0 ; z0 = 0.d0
!! position of the Earth !! 
x1 = 1.496e11 ; y1 = 0.d0 ; z1 = 0.d0 !
!! position of the test-particle !! 
x2 = 3e12 ; y2 = 0.d0 ; z2 = 0.d0 

r0 = sqrt(x0**2.d0+y0**2.d0+z0**2.d0)
r1 = sqrt(x1**2.d0+y1**2.d0+z1**2.d0)
r2 = sqrt(x2**2.d0+y2**2.d0+z2**2.d0)

!! velocity of the Sun !! (in m/s)
vx0 =  0.d0 ; vy0 =  0.d0; vz0 = 0.d0
!! velocity of the Earth !! 
vx1 = 0.d0; vy1 =  29.8e3 ; vz1 = 0.d0 !  
!! velocity of the test-particle !! 
vx2 = 0.d0; vy2 =  sqrt((G*(m0+m1+m2))/r2) ; vz2 = 0.d0 !  


!! time step size and number of time step !!
a_step=dfloat(n_step)
dt=t_per/a_step
k_max=n_periods*n_step

!! move to center of Mass and velocity !!
comx =  (m0*x0+m1*x1+m2*x2)/(m0+m1+m2)
comy =  (m0*y0+m1*y1+m2*y2)/(m0+m1+m2)
comz =  (m0*z0+m1*z1+m2*z2)/(m0+m1+m2)

covx =  (m0*vx0+m1*vx1+m2*vx2)/(m0+m1+m2)
covy =  (m0*vy0+m1*vy1+m2*vy2)/(m0+m1+m2)
covz =  (m0*vz0+m1*vz1+m2*vz2)/(m0+m1+m2)

x0 = x0 - comx ; x1 = x1 - comx ; x2 = x2 - comx
y0 = y0 - comy ; y1 = y1 - comy ; y2 = y2 - comy
z0 = z0 - comz ; z1 = z1 - comz ; y2 = y2 - comy

vx0 = vx0 - covx ; vx1 = vx1 - covx ; vx2 = vx2 - covx
vy0 = vy0 - covy ; vy1 = vy1 - covy ; vx2 = vx2 - covx 
vz0 = vz0 - covz ; vz1 = vz1 - covz ; vx2 = vx2 - covx


!! declare the vector of position for the objects !! 
x(1)=x0
x(2)=vx0
x(3)=y0
x(4)=vy0

x(5)=x1
x(6)=vx1
x(7)=y1
x(8)=vy1

x(9)=x2
x(10)=vx2
x(11)=y2
x(12)=vy2

t=0.d0

DO k=1,k_max    ! -------------------------->
  CALL kutta(x,t,dt)
  IF(MOD(k,n_step) == 0) THEN
    WRITE(23,*) x(5),x(7) ! test print: orbit of the Earth
  END IF 
END DO  !  <---------------------


END PROGRAM bodies_rk

SUBROUTINE kutta(x,t,dt)
implicit none
REAL*8, INTENT(OUT)                        :: x(12)
REAL*8, INTENT(IN OUT)                     :: t
!----------------------------------------------------------
!     Subroutine  Runge-Kutta 4th Explicit
!----------------------------------------------------------
REAL*8 :: k1(12), k2(12), k3(12), k4(12)
real*8 :: told,dt,h
integer :: j
REAL*8, DIMENSION(12) ::f

told=t
call phi(t,x(1),x(2),x(3),x(4),     &
&          x(5),x(6),x(7),x(8),     &
&          x(9),x(10),x(11),x(12),  &
&          f(1),f(2),f(3),f(4),     &
&          f(5),f(6),f(7),f(8),     &
&          f(9),f(10),f(11),f(12))

DO j=1,12
  k1(j)=f(j)
END DO

h=0.5d0*dt
call phi(t+h,x(1)+h*k1(1),x(2)+h*k1(2),x(3)+h*k1(3),x(4)+h*k1(4),  & 
&            x(5)+h*k1(5),x(6)+h*k1(6),x(7)+h*k1(7),x(8)+h*k1(8),  & 
&            x(9)+h*k1(9),x(10)+h*k1(10),x(11)+h*k1(11),x(12)+h*k1(12),  & 
&            f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10),f(11),f(12))

DO j=1,12
  k2(j)=f(j)
END DO

call phi(t+h,x(1)+h*k2(1),x(2)+h*k2(2),x(3)+h*k2(3),x(4)+h*k2(4),  & 
&            x(5)+h*k2(5),x(6)+h*k2(6),x(7)+h*k2(7),x(8)+h*k2(8),  & 
&            x(9)+h*k2(9),x(10)+h*k2(10),x(11)+h*k2(11),x(12)+h*k2(12),  & 
&            f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10),f(11),f(12))

DO j=1,12
  k3(j)=f(j)
END DO

call phi(t+dt,x(1)+dt*k3(1),x(2)+dt*k3(2),x(3)+dt*k3(3),x(4)+dt*k3(4),  &
&             x(5)+dt*k3(5),x(6)+dt*k3(6),x(7)+dt*k3(7),x(8)+dt*k3(8),  &
&             x(9)+dt*k3(9),x(10)+dt*k3(10),x(11)+dt*k3(11),x(12)+dt*k3(12),  &
&             f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10),f(11),f(12))

DO j=1,12
  k4(j)=f(j)
END DO

DO j=1,12
  x(j)=x(j)+(k1(j)+2.d0*k2(j)+2.d0*k3(j)+k4(j))*dt/6.d0
END DO

t=t+dt

RETURN
END SUBROUTINE kutta


subroutine phi(t,x0,vx0,y0,vy0,  &
&                fx0,fvx0,fy0,fvy0,   &
&                x1,vx1,y1,vy1,  &
&                fx1,fvx1,fy1,fvy1,   &
&                x2,vx2,y2,vy2,  &
&                fx2,fvx2,fy2,fvy2)
use constant
implicit none
REAL*8, INTENT(IN OUT)              :: t
REAL*8, INTENT(IN)                  :: x1,x0,x2
REAL*8, INTENT(IN)                  :: vx1,vx0,vx2
REAL*8, INTENT(IN)                  :: y1,y0,y2
REAL*8, INTENT(IN)                  :: vy1,vy0,vy2
REAL*8, INTENT(OUT)                 :: fx1,fx0,fx2
REAL*8, INTENT(OUT)                 :: fvx1,fvx0,fvx2
REAL*8, INTENT(OUT)                 :: fy1,fy0,fy2
REAL*8, INTENT(OUT)                 :: fvy1,fvy0,fvy2
real*8 :: r0,r1,r2
!----------
!       Define the components of the derivative of the vectorial field to solve:
!      dx/dt=fx(t,x,vx,y,vy)  dvx/dt=fvx(t,x,vx,y,vy)
!      dy/dt=fy(t,x,vx,y,vy)  dvy/dt=fvy(t,x,vx,y,vy)
!----------

r0=SQRT(((x2-x1)**2.d0+(y2-y1)**2.d0)**3.d0)
r1=SQRT(((x0-x2)**2.d0+(y0-y2)**2.d0)**3.d0)
r2=SQRT(((x1-x0)**2.d0+(y1-y0)**2.d0)**3.d0)

fx0=vx0
fy0=vy0
fvx0 = -G*((m1*(x0-x1))/r2**3.d0)-G*((m2*(x0-x2))/r1)
fvy0 = -G*((m1*(y0-y1))/r2**3.d0)-G*((m2*(y0-y2))/r1)

fx1=vx1
fy1=vy1
fvx0 = -G*((m2*(x1-x0))/r0**3.d0)-G*((m0*(x1-x2))/r2)
fvy0 = -G*((m2*(y1-y0))/r0**3.d0)-G*((m0*(y1-y2))/r2)

fx2=vx2
fy2=vy2
fvx2 = -G*((m0*(x2-x0))/r1**3.d0)-G*((m1*(x2-x1))/r0)
fvy2 = -G*((m0*(y2-y0))/r1**3.d0)-G*((m1*(y2-y1))/r0)

RETURN
END subroutine phi

在这个版本的代码中,我只是不删除子程序中的INTENT。也许 INTENT 警告与代码的错误输出无关,但至少我可以解决其中一个问题。

非常感谢。

fvx1fvy1 实际上从未在子例程 phi 中使用过——这就是编译器所抱怨的。如果您的代码依赖于这些值,它肯定会产生错误的结果。

顺便说一句:t+ht+dt 是常量表达式(你计算两个和),所以你不能将它用作意图为 out 或 [=16 的虚拟参数=],就像子例程 phi 中的 t 一样(尽管 t 也从未在该例程中使用过)。