program concentration

  implicit none

  real :: L,T,D,C0,C1
  integer :: N,Nt,j
  real, dimension(:), allocatable :: x_reg, C

  call recup_donnees(L,T,D,C0,C1,N,Nt)
  allocate(x_reg(N),C(N))
  call mesh(N,L,x_reg)
  do j=1,Nt
    call calcul_c(N,L,real(j)*T/real(Nt),D,C0,C1,x_reg,C)
    call ecriture(N,x_reg,C)
  end do
  deallocate(x_reg,C)

end program concentration


!-------------------------
function conc(x,t,D,C0,C1)

  implicit none
  real, intent(in) :: x,t,D,C0,C1
  real :: conc

  conc = C0 + (C1-C0)*(1.-erf(x/(2.*sqrt(D*t))))

end function conc


!--------------------------------------
subroutine recup_donnees(L,T,D,C0,C1,N,Nt)

        implicit none
        integer, intent(out) :: N,Nt
        real, intent(out) :: L,T,D,C0,C1

        print*, "Donner la valeur de la longueur L:"
        read*, L
        print*, "Donner la valeur du temps T:"
        read*, T
        print*, "Donner la valeur de diffusivite D:"
        read*, D
        print*, "Donner la valeur de la concentration initiale C0:"
        read*, C0
        print*, "Donner la valeur de la concentration a gauche C1:"
        read*, C1
        print*, "Donner la valeur du nombre de points N:"
        read*, N
        print*, "Donner la valeur d'echantillons temporels Nt:"
        read*, Nt

end subroutine recup_donnees


!-------------------------
subroutine mesh(N,L,x_reg)

        implicit none
        ! Variables en arguments
        integer, intent(in) :: N
        real, intent(in) :: L
        real, dimension(N), intent(out) :: x_reg
        ! Variables locales
        integer :: k

        do k=1,N
          x_reg(k) = L / real(N-1) * real(k-1)
        end do

end subroutine mesh

!-----------------------------------------
subroutine calcul_c(N,L,t,D,C0,C1,x_reg,C)

        implicit none
        ! Variables en arguments
        integer, intent(in) :: N
        real, intent(in) :: L,t,D,C0,C1
        real, dimension(N), intent(in) :: x_reg
        real, dimension(N), intent(out) :: C
        ! Variables locales
        integer :: k
        real :: conc

        do k=1,N
           C(k) = conc(x_reg(k),t,D,C0,C1)
        end do

end subroutine calcul_c

!-----------------------------
subroutine ecriture(N,x_reg,C)

        implicit none
        ! Variables en arguments
        integer, intent(in) :: N
        real, dimension(N), intent(in) :: x_reg,C
        ! Variables locales
        integer :: k

        open(10,file='resultats.dat',position='append')
        do k=1,N
           write(10,*) x_reg(k),C(k)
        end do
        !write(10,*) " "
        close(10)

end subroutine ecriture