!      program blasius
use msflib !  systemqq

logical calling
	parameter (nmax=20000)
!	dimension tp(nelmax)
!	external fcn,rkqc
	common /param/ eta0,deta,etaend,f10,f20
	common /noeq/ n	
	common /noeqt/ nt	
	common /div/ ndiv	
	common /prandtl/ pr

	common /etavec/ eta(nmax)	
	common /f1vec/ f1(nmax)	
	common /param1/ t10

	external ff,tf

 !     common /const/ ht2 
      open(1,file='inpclass.txt')
	open(2,file='f1.txt')
	open(3,file='f2.txt')
	open(4,file='f3.txt')
	open(5,file='t1.txt')
	open(6,file='t2.txt')
	open(7,file='f.txt')
	open(8,file='t.txt')

!	Solution of a Boundary Value Problem - BVP
!
!	data input
!
	read(1,*)n
	write(*,*)'n=',n
	read(1,*)eta0
	write(*,*)'eta0=',eta0
	read(1,*)etaend
	write(*,*)'etaend=',etaend
	read(1,*)deta
	write(*,*)'deta=',deta
	read(1,*)f10
	write(*,*)'f10=',f10
	read(1,*)f20
	write(*,*)'f20=',f20
	read(1,*)f30a
	write(*,*)'f30a=',f30a
	read(1,*)f30b
	write(*,*)'f30b=',f30b
	read(1,*)maxit
	write(*,*)'maxit=',maxit
	read(1,*)tol
	write(*,*)'tol=',tol
	read(1,*)nt
	write(*,*)'nt=',nt
	read(1,*)pr
	write(*,*)'pr=',pr
	read(1,*)t10
	write(*,*)'t10=',t10
	read(1,*)t20a
	write(*,*)'t20a=',t20a
	read(1,*)t20b
	write(*,*)'t20b=',t20b

	! Blasius solution

	call secante(f30a,f30b,maxit,tol,c,ff)

	! Energy equation

	ndiv=etaend/deta
	write(*,*)'ndiv=',ndiv

	call secante(t20a,t20b,maxit,tol,d,tf)

	write(*,*) 'f30=',c
	write(*,*) 't20=',d

	close(2)
	close(3)
	close(4)
	close(5)
	close(6)
	close(7)
	close(8)

calling = systemqq('notepad f1.txt')       ! list of data

calling = systemqq('notepad f2.txt')       ! list of data
calling = systemqq('notepad f3.txt')       ! list of data

calling = systemqq('notepad t1.txt')       ! list of data

calling = systemqq('notepad t2.txt')       ! list of data

calling = systemqq('notepad f.txt')       ! list of f-table

calling = systemqq('notepad t.txt')       ! list of t-table

calling = systemqq('wgnuplot data.gnu')  ! graph - velocities
calling = systemqq('wgnuplot datat.gnu')  ! graph - temperatures

	stop
	end
!**************************************
	subroutine fcn(n,t,fi,f,nelmax)
	dimension fi(nelmax),f(nelmax)
!
	f(1)=fi(2)
	f(2)=fi(3)
	f(3)=-fi(1)*fi(3)/2
!
	return
	end
!***************************************
!**************************************
	subroutine tfcn(n,t,fi,f,nelmax)
	dimension fi(nelmax),f(nelmax)
	common /etavec/ eta1(1)	
	common /f1vec/ f1(1)	
	common /div/ ndiv	
	common /param/ eta0,deta,etaend,f10,f20
	common /prandtl/ pr
!
!	identify appropriate position of eta in f-table
!
!	write(*,*)'ndiv=',ndiv

	ipos=t/deta+1

!	 write (*,*) ipos,t,eta1(ipos),f1(ipos)
!
	f(1)=fi(2)
	f(2)=-pr*f1(ipos)*fi(2)/2
!
	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 
!------------------------------------------------------------------ 
	function ff(x)
	parameter (nmax=100)
	dimension tp(nmax)
	external fcn,rkqc
	common /noeq/ n
	common /flag/ iflag
	common /etavec/ eta1(1)	
	common /f1vec/ f1(1)	


	common /param/ eta0,deta,etaend,f10,f20
!
!	initial values
!
	eta=eta0
	tp(1)=f10
	tp(2)=f20
	tp(3)=x
!
	k=0 
	if (iflag.eq.1) then
	eta1(k+1)=eta
	f1(k+1)=tp(1)
	write(2,*)eta,tp(1)
	write(3,*)eta,tp(2)
	write(4,*)eta,tp(3)
	write(7,*)'--------------------------------------------'

	write(7,*)'Nr do passo   eta     f1     f2      f3'

	write(7,*)'--------------------------------------------'

	write(7,*)k,eta,(tp(l),l=1,n)
	endif

	ht2=deta/1
	nd=nmax
!
!	beginning of time loop
!
 50	k=k+1 
	tendi=eta+deta
!      	write(*,*)'-------------eta=',tendi
!
!     call odeint(tp,n,time,tendi,1.e-6,ht2,1.e-20,id1,id2,nelmax,fcn,rkqc)
	call  rk4ord(tp,n,eta,deta,fcn,nd)
!    call fore(n,fcn,time,tp,tendi,nmax)
	if (iflag.eq.1) then
	if (k.eq.100.or.k.eq.200.or.k.eq.300.or.k.eq.400.or.k.eq.500.or.k.eq.600.or.k.eq.700.or.k.eq.800.or.k.eq.900.or.k.eq.1000) write(7,*)k,tendi,(tp(l),l=1,n)
	eta1(k+1)=tendi
	f1(k+1)=tp(1)
	write(2,*)tendi,tp(1)
	write(3,*)tendi,tp(2)
	write(4,*)tendi,tp(3)
	endif
!
     	if (tendi.lt.etaend) then
			eta=tendi
			goto 50 
		endif

	aux=tp(2)-1.
	ff=aux


      return 
      end 
!------------------------------------------------------


