module dpres
  integer, parameter :: dp = kind(1.D0)
end module dpres

subroutine decoupe(iCPU,NCPU,N,NP,dx,x0)
use dpres
implicit none
integer :: iCPU,NCPU,N,NP,reste
real(kind=dp) :: dx,x0

NP = N/NCPU
reste = mod(N,NCPU)
if(iCPU<reste) then
  NP = NP + 1
  x0 = -1.d0 + dfloat(NP*iCPU)*dx
else
  x0 = -1.d0 + dfloat((NP+1)*reste + NP*(iCPU-reste))*dx
endif
end subroutine decoupe

!routine d'initialisation des champs
subroutine init(N,dx,x0,x,u)
use dpres
implicit none
integer :: N,i
real(kind=dp) :: dx, pi = 4.d0*datan(1.d0),x0
real(kind=dp), dimension(0:N-1) :: x,u

do i=0,N-1 !boucle sur les points de calculs
  x(i) =  x0 + dfloat(i)*dx ! x € [-1,1[
  u(i) =  1.5d0 + dcos(2.d0*pi*x(i)) ! u : trigonometrique
enddo
end subroutine init
! integration temporelle suivant le shema d'Euler
! Equation du/dt = f(u)
! resolue comme u(t+dt) = u(t) + dt.f(u(t))
subroutine integre(N,dt,v,u,du)
use dpres
implicit none
integer :: N,i
real(kind=dp) :: dt,v
real(kind=dp), dimension(0:N-1) :: u,du

do i=0,N-1 !boucle sur les points de calculs
  u(i) = u(i) - dt*v*du(i); !shema euler
enddo
end subroutine integre
! calcul des derivees spatiales
! shema aux differences finies decentrees gauche
! du/dx = (u(x) - u(x-dx))/dx
subroutine derive(iCPU,NCPU,N,dx,u,du)
use dpres
implicit none
include "mpif.h"
integer :: N,i
real(kind=dp) :: dx
real(kind=dp), dimension(0:N-1) :: u,du
! variables paralleles
integer :: iCPU,NCPU,erreur,source,dest,flag
real(kind=dp) :: ug
integer :: statu(MPI_STATUS_SIZE)

source = iCPU - 1
dest   = iCPU + 1
if(iCPU == 0  ) source = NCPU - 1
if(iCPU == NCPU-1) dest   = 0

if(mod(iCPU,2)==0)then
  call MPI_SEND(u(N-1),1,MPI_DOUBLE_PRECISION,dest,100 + iCPU,MPI_COMM_WORLD,erreur)
  call MPI_RECV(ug,1,MPI_DOUBLE_PRECISION,source,100+source,MPI_COMM_WORLD,statu,erreur)
else
  call MPI_RECV(ug,1,MPI_DOUBLE_PRECISION,source,100+source,MPI_COMM_WORLD,statu,erreur)
  call MPI_SEND(u(N-1),1,MPI_DOUBLE_PRECISION,dest,100 + iCPU,MPI_COMM_WORLD,erreur)
endif

do i=1,N-1 ! boucle sur les points de calculs
  du(i) = (u(i) - u(i-1))/dx !diff. finies
enddo
du(0) = (u(0) - ug)/dx ! condition aux limites periodiques
end subroutine derive      ! u(-1) = u(1)
!subroutine de sauvegarde dans ficihier
subroutine ecrit(ifile, N,x,u)
use dpres
implicit none
integer :: ifile,N,i
real(kind=dp), dimension(0:N-1) :: x,u

do i =0,N-1 !boucle sur les points de calculs
  write(ifile,*) x(i),u(i) !ecriture
enddo
close(ifile)
end subroutine ecrit
!program principal
program convection
use dpres
implicit none
include "mpif.h"
integer :: N,it,Nt ! N: nombre de points de calculs, Nt: nombre de pas de temps
real(kind=dp) :: v = 3.42d0,dx,dt ! vitesse de convection, pas d'espace et de temps
real(kind=dp), dimension(:), allocatable :: x,u,du

integer :: myID, nbCPU, erreur, NP;
real(kind=dp) :: x0
character :: f_res*13,f_init*13,conc_CPU*3

call MPI_INIT(erreur);
call MPI_COMM_SIZE(MPI_COMM_WORLD,nbCPU,erreur)
call MPI_COMM_RANK(MPI_COMM_WORLD,myID ,erreur)


N = 10000 ! maillage d'espace
Nt= 10000 ! nombre de pas de temps
dx = 2.d0/dfloat(N) ! pas d'espace
dt = 0.1d0*dx/v   ! pas de temps de stabilite

call decoupe(myID,nbCPU,N,NP,dx,x0)
allocate(x(0:NP-1),u(0:NP-1),du(0:NP-1)) ! allocation dynamique
call init(NP,dx,x0,x,u) ! initialisation
! ouverture des fichiers et sauvegarde
write(conc_CPU,'(i3.3)') myID
f_init = 'initf_'//conc_CPU//'.dat'
open(unit=11, file=f_init, status="unknown")
call ecrit(11,NP,x,u)

do it = 0,Nt-1 ! boucle sur les pas de temps
  call derive(myID,nbCPU,NP,dx,u,du) ! calcul de gradients
  call integre(NP,dt,v,u,du) ! integration temporelle 
enddo
!sauvegarde du resultat
f_res = 'resf_'//conc_CPU//'.dat'
open(unit=12, file=f_res, status="unknown")
call ecrit(12,NP,x,u)
!desallocation
deallocate(x,u,du)
call MPI_Finalize(erreur)
end program convection
