! Copyright (C) 2016      David Tsiklauri <d.tsiklauri@qmul.ac.uk>
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details:
! <http://www.gnu.org/licenses/>.

MODULE constants  
INTEGER, PARAMETER :: nx=20000, ns=20
INTEGER, PARAMETER :: num = KIND(1.D0)
REAL(num),PARAMETER :: pi = 3.14159265358979_num
END MODULE constants

MODULE functions  
USE constants
IMPLICIT NONE
CONTAINS

PURE FUNCTION dfdx2(t,x,dx,f)  
REAL(num), INTENT(IN) :: t,dx
REAL(num), dimension (-1:nx+2),INTENT(IN) :: f,x
REAL(num), dimension (-1:nx+2) :: dfdx2
INTEGER :: i

do i=1,nx
dfdx2(i) = (-(1.0_num/12.0_num)*f(i+2)+(4.0_num/3.0_num)*f(i+1)-(5.0_num/2.0_num)*f(i)+(4.0_num/3.0_num)*f(i-1)-(1.0_num/12.0_num)*f(i-2))/dx**2

!two options: 1st is 2th order spatial derivative for diffusion equation and option 2 is 3rd order spatial derivative for linear KdV equation
!option 1 (-(1.0_num/12.0_num)*f(i+2)+(4.0_num/3.0_num)*f(i+1)-(5.0_num/2.0_num)*f(i)+(4.0_num/3.0_num)*f(i-1)-(1.0_num/12.0_num)*f(i-2))/dx**2
!option 2  -(0.5_num*f(i+2)-f(i+1)+f(i-1)-0.5_num*f(i-2))/dx**3
!remeber to alter power of dt in line with dt=0.1_num*dx**2 for option 1 and dt=0.1_num*dx**3 for option 2

end do
END FUNCTION dfdx2

SUBROUTINE bc(f)  
USE constants
REAL(num), dimension (-1:nx+2) :: f

f(-1)  =f(nx-2) 
f(0)   =f(nx-1)      
f(nx+1)=f(2)
f(nx+2)=f(3)

END SUBROUTINE bc
end module functions

PROGRAM kdv
USE constants
use functions

IMPLICIT NONE
CHARACTER (len = 6) :: str_ts 
REAL(num), DIMENSION(-1:nx+2) :: f,fp1,fp2,fp3,x
REAL(num) :: xs,xe,dx,t1,t2,t,dt

INTEGER :: cc,i,nt

xs=-250.0_num
xe=250.0_num

dx=(xe-xs)/(real(nx-1.0_num))

t1=0.0_num
t2=12.50_num

dt=0.1_num*dx**2

nt=int((t2-t1)/dt)
print*,dx
print*,dt

do i=-1,nx+2 
x(i)=xs+(i-1)*dx
enddo

f=0.0_num
fp1=0.0_num
fp2=0.0_num
fp3=0.0_num

do i=1,nx
f(i)=exp(-x(i)**2) 
!initial condition
enddo

t=t1

if (mod(cc,nt/ns)==0) then
write (str_ts, '(a, i2.2)') "",cc/(nt/ns)
open (unit = 82, file = ""//str_ts, form = 'unformatted', status = 'replace')
write(82) f,x,t
close(82)
print *,t,cc/(nt/ns),maxval(abs(f))
end if


do while (t .le. t2) 

fp1=f+0.5_num*dt*dfdx2(t,x,dx,f)
call bc(fp1)
 
fp2=f+0.5_num*dt*dfdx2(t+0.5_num*dt,x,dx,fp1)
call bc(fp2)

fp3=f+dt*dfdx2(t+0.5_num*dt,x,dx,fp2)
call bc(fp3)

f=f+dt/6.0_num*((fp1-f)/(0.5_num*dt)+2.0_num*(fp2-f)/(0.5_num*dt)+2.0_num*(fp3-f)/dt+dfdx2(t+dt,x,dx,fp3))
call bc(f)

t=t+dt
cc=cc+1

if (mod(cc,nt/ns)==0) then
write (str_ts, '(a, i2.2)') "",cc/(nt/ns)
open (unit = 82, file = ""//str_ts, form = 'unformatted', status = 'replace')
write(82) f,x,t
close(82)
print *,t,cc/(nt/ns),maxval(abs(f))
end if

enddo
END PROGRAM kdv