!____. 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=100,Ng=2
!____.__________________________________________________________________
REAL,DIMENSION(1-Ng:Nx+Ng):: v_new,v_old,v_exact,x
REAL:: delta_t,delta_x,a_conv,rhs,rhs_uw,rhs_uw2
REAL:: error,error_max,time,v0
!____.__________________________________________________________________
!     metric definition
CALL griglia(Nx,Ng,x,delta_x)
!____.__________________________________________________________________
!     some initializations
DO j=1-ng,nx+ng
 v_old(j)=v0(x(j))
ENDDO
error=0.
error_max=-999999.9
!____.__________________________________________________________________
!     read on file eulero_stokes.dat
CALL input(itmin,itmax,itout,a_conv,delta_t)
!____.__________________________________________________________________
!     compute maximun number of iterations and open write unit
!____.__________________________________________________________________
write(*,*) 'cfl',delta_t*a_conv/delta_x
CALL write_disk(x,v_old,v_old,0,Nx,Ng)
DO it= itmin, itmax
  DO j=1,nx
   rhs = rhs_uw2(v_old(j),v_old(j-1),v_old(j-2),a_conv,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)=0.
  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,a_conv,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(*,*) 'solution stored in file velocity.dat'
WRITE(*,*) 'Simulation ended!'
STOP
END
!____.__________________________________________________________________
!____.__________________________________________________________________
!____.__________________________________________________________________
!____.__________________________________________________________________
FUNCTION rhs_uw2(vj,vjm,vjmm,a,Dx)
IMPLICIT NONE
REAL rhs_uw2,vj,vjm,vjmm,a,Dx

rhs_uw2=-a*(3.*vj-4.*vjm+vjmm)/(2.*Dx)

RETURN
END
!____.__________________________________________________________________
!____.__________________________________________________________________
FUNCTION v0(x)
IMPLICIT NONE
REAL v0,x,fx

fx=(.25-x)*(x-.5)

v0=max(0.,fx)

RETURN
END 
!____.__________________________________________________________________
!____.__________________________________________________________________
SUBROUTINE exact_sol(Nx,Ng,x,time,a_conv,v_exact)
 
IMPLICIT NONE 
INTEGER Nx,Ng,j
REAL time,x(1-ng:nx+ng),a_conv,v_exact(1-ng:nx+ng)
REAL v0, appo

DO j=1-ng,Nx+ng
appo=x(j)-a_conv*time
v_exact(j)=v0(appo)
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
  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)

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

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