PROGRAM main
IMPLICIT NONE
INTEGER,PARAMETER :: Nt=100000
INTEGER :: n
REAL :: v,t,tau,v0,g
REAL :: Dt,v_old,v_new
REAL :: rhs,t_max,t_min
REAL :: error,err_max

error   = 0.
err_max = 0.

CALL input(t_min,t_max,tau,g,v0)
Dt=(t_max-t_min)/float(Nt)
v_old=v0

CALL write_sol(t_min,v0,v0,error)

DO n=1,Nt
t=t_min+Dt*n
CALL exact_sol(v,t,v0,tau,g)

v_new = v_old + rhs(Dt,v_old,tau,g)


error = abs(v_new-v) 
err_max=max(err_max,error)

CALL write_sol(t,v,v_new,error)

v_old = v_new
ENDDO

WRITE(*,*) 'ERRORE MAX',Dt,err_max 
WRITE(*,*) 'Fine della simulazione' 

STOP
END PROGRAM main
!_____________________________________________________________
!_____________________________________________________________
SUBROUTINE exact_sol(v,t,v0,tau,g)
IMPLICIT NONE
REAL,INTENT(IN):: v0,t,tau,g
REAL,INTENT(OUT):: v

v = (v0-g*tau)*exp(-t/tau)+g*tau

RETURN
END SUBROUTINE exact_sol
!_____________________________________________________________
SUBROUTINE input(t_min,t_max,tau,g,v0)
IMPLICIT NONE
REAL,INTENT(OUT):: v0,tau,g,t_min,t_max

OPEN(UNIT=1,FILE='input.dat',STATUS='old',ACTION='read')
READ(1,*) t_min
READ(1,*) t_max
READ(1,*) v0
READ(1,*) tau
READ(1,*) g
CLOSE(1)
RETURN
END SUBROUTINE input
!_____________________________________________________________
FUNCTION rhs(dt,v,tau,g)
IMPLICIT NONE
REAL,INTENT(IN) :: dt,tau,g,v
REAL :: rhs

rhs = (-v/tau + g) * dt

RETURN
END FUNCTION rhs
!_____________________________________________________________
SUBROUTINE write_sol(x,ye,ya,err)
IMPLICIT NONE
REAL,INTENT(IN)::x,ye,ya,err

OPEN(UNIT=1,FILE='output.dat',STATUS='unknown',ACTION='write',POSITION='append')
WRITE(1,*) x,ye,ya,err
CLOSE(1)

RETURN
END SUBROUTINE write_sol
!_____________________________________________________________