!      program reservatorio
use msflib ! biblioteca que contem o comando systemqq

logical chamada
	parameter (nelmax=100)
	dimension tp(nelmax)
	external fcn,rkqc
      common /const/ ht2 
	common /param/ rps,vswept,r,rn,vs,cv,rar
	common /param2/ temp0,vres,p0
	common /param3/ vm

      open(1,file='inpaula.txt')
	open(2,file='outnum.txt')
	open(3,file='out-read.txt')
	open(10,file='out-m.txt')
	open(11,file='out-p.txt')
	open(12,file='out-vm.txt')
!
!	integrador de EDO's com passo fixo - RK 4a ordem
!
	read(1,*)n
	write(3,*)'n=',n
	read(1,*)tau0
	write(3,*)'tau0=',tau0
	read(1,*) tend
	write(3,*)'tend=',tend
	read(1,*) dtau
	write(3,*)'dtau=',dtau
	read(1,*) p0
	write(3,*)'p0=',p0
	read(1,*) iflag
	write(3,*)'iflag=',iflag
	read(1,*) ht2
	write(3,*)'ht2=',ht2
	read(1,*) rps
	write(3,*)'rps=',rps
	read(1,*) vswept
	write(3,*)'vswept=',vswept
	read(1,*) r
	write(3,*)'r=',r
	read(1,*) rn
	write(3,*)'rn=',rn
	read(1,*) rar
	write(3,*)'rar=',rar
	read(1,*) cv
	write(3,*)'cv=',cv
	read(1,*) vres
	write(3,*)'vres=',vres
	read(1,*) temp0
	write(3,*)'temp0=',temp0

!
!	compute constants
!
	rm0=p0*vres/rar/temp0
	write(3,*)'rm0=',rm0
	rhoar=rm0/vres
	write(3,*)'rhoar=',rhoar
	vs=1./rhoar
	write(3,*)'vs=',vs

!
!	initial values
!
	time=tau0
	tp(1)=rm0
	tp(2)=p0
!*****************
!
!
!	loop para apresentar resultados intermediarios 
!	em cada intervalo dtau durante a integracao das equacoes
!
	k=0 
	write(*,*) '     Table of results'
	write(*,*)'------------------------------------------'
	write(*,*)'  Passo Nr     time       x(i)     v(i)   '
	write(*,*)'------------------------------------------'
	write(*,*)k,time,(tp(l),l=1,n)
	write(2,*)time,(tp(l),l=1,n)
	write(10,*)time,tp(1)
	write(11,*)time,tp(2)
	write(12,*)time,vm
!
!
!	beginning of time loop
!
 50	k=k+1 
	tendi=time+dtau
      	write(*,*)'-------------time=',tendi
	if (iflag.eq.0) then ! resolve por RK adaptativo
	     call odeint(tp,n,time,tendi,1.e-6,ht2,1.e-20,id1,id2,nelmax,fcn,rkqc)
	else
		if (iflag.eq.1) then   ! resolve por RK passo fixo
			call  rk4ord(tp,n,time,dtau,fcn,nd)
		else	! resolve por Forward Euler
		    call fore(n,fcn,time,tp,tendi,nelmax)
		endif
	endif
	write(*,*)k,tendi,(tp(l),l=1,n)
	write(2,*)tendi,(tp(l),l=1,n)
	write(10,*)tendi,tp(1)
	write(11,*)tendi,tp(2)
	write(12,*)tendi,vm
!
     	if (tendi.lt.tend) then
			time=tendi
			goto 50 
		endif
	close(2)
	close(3)
	close(10)
	close(11)
	close(12)


chamada = systemqq('notepad outnum.txt')       ! listagem dos dados


chamada = systemqq('wgnuplot dados-m.gnu')  ! gráfico m

chamada = systemqq('wgnuplot dados-p.gnu')  ! gráfico p

chamada = systemqq('wgnuplot dados-vm.gnu')  ! gráfico vm

	stop
	end
!**************************************
	subroutine fcn(n,t,fi,f,nelmax)
	dimension fi(nelmax),f(nelmax)
	common /param/ rps,vswept,r,rn,vs,cv,rar
	common /param2/ temp0,vres,p0
	common /param3/ vm
!
	efv=1-r*((fi(2)/p0)**(1./rn)-1.)
	f(1)=vswept/vs*rps*efv*cv
	vm=f(1)
	f(2)=rar*temp0/vres*f(1)
	
	
	return
	end
!***************************************
!234567890123456789012345678901234567890123456789012345678901234567890 
      subroutine fore(n,fcn,time,fi,tend,nelmax) 
! 
!      implicit real *8 (a-h,o-z)
      parameter (nd1=100) 
      dimension fi(nelmax),f(nd1)
      common /const/ ht2 
      external fcn 
      k=0 
 50   k=k+1 
      time=min(time+ht2,tend) 
      call fcn(n,time,fi,f,nelmax) 
      do 100 i=1,n 
         fi(i)=fi(i)+ht2*f(i) 
 100  continue 
      if (time.lt.tend) goto 50 
      return 
      end 
!------------------------------------------------------------------ 





