!____. file: poisson.f90
!____. This program apply the eulero integration of 
!____. an ordinary differential equation 
PROGRAM poisson
!____.__________________________________________________________________
IMPLICIT NONE
INTEGER:: i,j
INTEGER:: ntime,nout
INTEGER,PARAMETER:: Nx=100,Ng=1
!____.__________________________________________________________________
REAL,DIMENSION(1-ng:Nx+ng):: v_new, v_old, v_exact, x
REAL:: delta_t,delta_x,rhs,rhs_cent
REAL:: time,cfl
REAL:: error,error_max
!____.__________________________________________________________________
!     metric definition
call griglia(Nx,Ng,x,delta_x)
!____.__________________________________________________________________
!     some initializations
call exact_sol(Nx,Ng,x,0.,v_exact)

DO i=1-Ng,Nx+Ng
  v_old(i)=0.
ENDDO
error=0.
error_max=-999999.9
!____.__________________________________________________________________

!call jacobi(nx,ng,v_old,v_new,delta_x)
call gauss_seidel(nx,ng,v_old,v_new,delta_x)

error=0.
DO j=1,Nx
  error = error + abs(v_new(j)-v_exact(j))
ENDDO
error_max=max(error_max,error/float(Nx))

!____.  write on disk
call write_disk(x,v_new,v_exact,0,Nx,Ng)

WRITE(*,*) 'dt,   error=',delta_t,error_max
WRITE(*,*) 'solution stored in file velocity.dat'
WRITE(*,*) 'the end falks!'
STOP
!____.__________________________________________________________________
END
!____.__________________________________________________________________
!____.__________________________________________________________________
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE jacobi(nx,ng,v_old,v_new,delta_x)
!____.
IMPLICIT NONE
INTEGER:: nx,ng,i,j,m,itmax
REAL,DIMENSION(1-Ng:Nx+Ng):: v_new,v_old
REAL:: tol,error,appo,delta_x

tol =1.E-10
m=0

DO 
error=0.
 DO j=1,nx
   v_new(j)=(-delta_x**2+v_old(j-1)+v_old(j+1))/2.
   error = error + abs(v_new(j)-v_old(j))
 ENDDO
do j=1-ng,0
   v_new(j)=0.
   v_new(-j+nx+1)=0.
enddo
    if (error/float(Nx).lt.tol) then
       exit
    else 
       DO j=1,nx
        v_old(j)=v_new(j)
       ENDDO
    endif
   m = m + 1
ENDDO
write(*,*) 'it=',m,'error=',error
!____.
RETURN
END SUBROUTINE jacobi
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE gauss_seidel(nx,ng,v_old,v_new,delta_x)
!____.
IMPLICIT NONE
INTEGER:: nx,ng,i,j,m,itmax
REAL,DIMENSION(1-ng:nx+ng):: v_new, v_old
REAL:: tol,error,appo,delta_x

tol =1.E-10
itmax=40
m=0


write(*,*) 'sono qui'
DO 
error=0.
 DO j=1,nx
   v_new(j)=(-delta_x**2+v_new(j-1)+v_old(j+1))/2.
   error = error + abs(v_new(j)-v_old(j))
 ENDDO
    if (error.lt.tol) then
       exit
    else 
       DO j=1,nx
        v_old(j)=v_new(j)
       ENDDO
    endif
   m = m + 1
ENDDO
write(*,*) 'it=',m,'error=',error
!____.
RETURN
END
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE exact_sol(Nx,Ng,x,time,v_exact)
 
IMPLICIT NONE 
INTEGER:: Nx,Ny,Ng,i,j,WN,m,n
REAL,DIMENSION(1-ng:nx+ng)::x,v_exact
REAL:: appo, appo1, appo2, pi, time

WN=30.
pi = acos(-1.)

DO i=1-Ng,Nx+Ng
  v_exact(i)= 0.5*(x(i)**2-1.)
ENDDO

RETURN
ENDSUBROUTINE exact_sol
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE write_disk(x,v,u,it,nx,ng)
 
IMPLICIT NONE 
INTEGER:: it,nx,ny,ng,i,j
REAL,DIMENSION(1-ng:Nx+ng):: x,v,u
CHARACTER(32):: opfile

opfile='acl12345678.dat'
WRITE(opfile(4:11),1000) it
1000  FORMAT(i8.8)

OPEN(UNIT=13,FILE=opfile,STATUS='unknown',ACTION='write')

 DO i=1-Ng,Nx+ng
  WRITE(13,*) x(i),v(i),u(i)
 ENDDO

CLOSE(13)

RETURN
END SUBROUTINE write_disk
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE griglia(Nx,Ng,x,delta_x)
IMPLICIT NONE
INTEGER Nx,j,Ng
REAL x(1-ng:Nx+ng),delta_x
REAL L
L=2.
delta_x=L/float(Nx+1)


DO j=1-ng,Nx+ng
  x(j)=-1.+j*delta_x
ENDDO

RETURN
ENDSUBROUTINE griglia














