Pierre,
You have a number of syntax problems. To the best of my knowledge,
dt=1.0/(20.0E9)_2
should be
dt=1.0_2 / 20.0D9
It would certainly be preferred.
20.0E9 implies real4, while 20.0D9 is real8
1.0 is real4 while 1.0_2 or 1.0D0 is real8, although for 0.5,1,2,3 or 8 there is no difference, unlike 0.1 which is not a power of 2.
I don't think you can have ')_2' to imply kind=2 (real*8)
I have changed your program and use ..d0 for real*8 and also changed the write statements to put one value per line. This is just my preferred style, which you are free to ignore.
I also introduced 'k' to take array syntax out of EXP. Again just my style preference.
The program ran ok using either
ftn95 pierre /debug /link
or
ftn95 pierre /check /link
It takes about 400mb of memory to run.
PROGRAM test
!
integer j,N,time,l , k
Real8 c,Pi,dt,dx,a,b,si,ti ,t
parameter (c=299792458.0_2)
parameter (Pi=3.14159265358979323846d0)
parameter (time=8001) !you can try time=801
parameter (N=6004) ! you can try N=604
REAL8, DIMENSION(0:N-1,-1:time) :: Phi
!Phi(position,time)=Phi(i,j) time=0 ⇒ t=0s
REAL8, DIMENSION(0:N-1,1:5) :: Y
REAL8, DIMENSION(0:N-1) :: X
REAL*8, DIMENSION(0:N-1,1:5) :: An
REAL START,FINISH
! beginning of the calculus
CALL CLOCK@(START)
dt = 1.0d0 / 20.0d9
dx = c / 10.0d9
b = c2 * dt2 / dx**2
a = 2.0d0 * (1.0d0-b)
si = 3.0d0 / (2.0d0Pi1.0d9)
! analytical solution
An = 0
Do j=0,N-1
x(j) = jdx
end do
do l=1,5
ti = l dt1600.0d0 ! you can try ti=ldt*160.0_2
do k = 0,n-1
An(k,l)= exp(-1.0d0/2.0d0 * (((ti-x(k)/c)/si)-8.0d0)**2 )
end do
end do
Phi(:,-1:0) = 0
! Absorbing boundary condition.
Do j=1,time
t = jdt
Phi(0,j) = exp (-1.0d0/2.0d0((t/si)-8.0d0)*2)
do k = 1,n-2
Phi(k,j) = -Phi(k,j-2) + aPhi(k,j-1) + b(Phi(k+1,j-1) + Phi(k-1,j-1))
end do
Phi(N-1,j) = Phi(N-2,j-1) - 1.0d0/3.0d0(Phi(N-2,j) - Phi(N-1,j-1))
end do
!end of the calculus
CALL CLOCK@(FINISH)
PRINT *,'seconds used = ',FINISH-START
y = 0
y(:,1) = Phi(:,1600) ! you can try y(:,1)=Phi(:,160)
y(:,2) = Phi(:,3200) ! you can try y(:,2)=Phi(:,320)
y(:,3) = Phi(:,4800) ! you can try y(:,3)=Phi(:,480)
y(:,4) = Phi(:,6400) ! you can try y(:,4)=Phi(:,640)
y(:,5) = Phi(:,8000) ! you can try y(:,5)=Phi(:,800)
!output to a .txt file
OPEN (UNIT=5,FILE='X.txt')
do k = 0,n-1
write(5,*) k, x(k)
end do
Close(5)
!output to a .txt file
OPEN (UNIT=5,FILE='Y.txt')
do j = 1,5
do k = 0,n-1
write(5,*) j,k,y(k,j)
end do
end do
Close(5)
!output to a .txt file
OPEN (UNIT=5,FILE='An.txt')
do j = 1,5
do k = 0,n-1
write(5,*) j,k,An(k,j)
end do
end do
Close(5)
END PROGRAM test
