!____. file: eulero_upwind.f
!____. This program apply the eulero integration of 
!____. an ordinary differential equation
PROGRAM main
!____.__________________________________________________________________
IMPLICIT NONE
INTEGER:: it, j
INTEGER:: itmin,itmax,itout
INTEGER,PARAMETER:: Nx=20,Ng=1
!____.__________________________________________________________________
REAL,DIMENSION(1-Ng:Nx+Ng):: v_new,v_old,v_exact,x
REAL:: delta_t,delta_x,nu,rhs_cent,rhs,va,vb
REAL:: error,error_max,time,v0
!____.__________________________________________________________________
!     metric definition
CALL griglia(Nx,Ng,x,delta_x)
!____.__________________________________________________________________
!     read on file eulero_stokes.dat
CALL input(itmin,itmax,itout,nu,delta_t,va,vb)
!____.__________________________________________________________________
!     some initializations
call exact_sol(Nx,Ng,x,0.,nu,va,vb,v_exact)
DO j=1-ng,nx+ng
 v_old(j)=v_exact(j)
ENDDO
error=0.
error_max=-999999.9
!____.__________________________________________________________________
!     compute maximun number of iterations and open write unit
!____.__________________________________________________________________
write(*,*) 'cfl',delta_t*nu/delta_x**2
CALL write_disk(x,v_old,v_old,0,Nx,Ng)
DO it= itmin, itmax
  DO j=1,nx
   rhs = rhs_cent(v_old(j+1),v_old(j),v_old(j-1),nu,delta_x)
   v_new(j) = v_old(j) + delta_t * rhs
  ENDDO
  time = it*delta_t
!____.  boundary conditions
  DO j=1-ng,0
    v_new(     j) = va
    v_new(Nx+1+j) = vb
  ENDDO
!____.   update the solution 
  DO j=1-ng,Nx+ng
   v_old(j)= v_new(j)
  ENDDO
!____.   compute the error 
  CALL exact_sol(Nx,Ng,x,time,nu,va,vb,v_exact)

  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
  if (mod(it,itout).eq.0) CALL write_disk(x,v_new,v_exact,it,Nx,Ng)
ENDDO

WRITE(*,*) 'dt, error=',delta_t,error_max
WRITE(*,*) 'Simulation ended!'
STOP
END
!____.__________________________________________________________________
!____.__________________________________________________________________
!____.__________________________________________________________________
!____.__________________________________________________________________
FUNCTION rhs_cent(vjp,vj,vjm,a,Dx)
IMPLICIT NONE
REAL rhs_cent,vjp,vj,vjm,a,Dx

rhs_cent=a*(vjp-2.*vj+vjm)/Dx**2

RETURN
END
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE exact_sol(Nx,Ng,x,time,nu,va,vb,v_exact)
 
IMPLICIT NONE 
INTEGER Nx,Ng,j, MW, m
REAL time,x(1-ng:nx+ng),nu,v_exact(1-ng:nx+ng)
REAL appo, va, vb, vs, appo2
REAL,PARAMETER:: pi=acos(-1.)

MW = 1000

DO j=1-ng,Nx+ng
vs = va + (vb - va) * x(j)
appo2 = 0.
DO m = 1, MW
appo = (2 * m -1) * pi
appo2 = appo2 + 400./appo * sin(appo*x(j))  &
      * exp(-nu*appo**2*time)
ENDDO

v_exact(j) = - appo2 + vs
ENDDO

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

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

OPEN(UNIT=13,FILE=opfile,STATUS='unknown',ACTION='write')
 DO j=1-Ng,Nx+ng
  WRITE(13,*) x(j),v(j),u(j)
 ENDDO
CLOSE(13)

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

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

RETURN
END
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE input(itmin,itmax,itout,nu,delta_t,va,vb)
IMPLICIT NONE
INTEGER :: itmin,itmax,itout
REAL :: nu,delta_t,va,vb
OPEN(UNIT=1,FILE='input.dat',STATUS='old',ACTION='read')
READ(1,*)itmin
READ(1,*)itmax
READ(1,*)itout
READ(1,*)nu
READ(1,*)delta_t
READ(1,*)va
READ(1,*)vb
CLOSE(1)
RETURN
END SUBROUTINE input









