!      program microalgae growth
use msflib ! biblioteca que contem o comando systemqq

logical chamada
	parameter (nelmax=100)
	dimension tp(nelmax)
	external fcn,rkqc
      common /const/ ht2 
	common /param/ ri0,p,rka,rmu_max,rn,rik,rm

      open(1,file='inpaula.txt')
	open(2,file='out-read.txt')
	open(4,file='out-C.txt')
!
!	integrador de EDO's com passo adaptativo - RK 4a/5a ordem
!	integrador de EDO's com passo fixo - RK 4a ordem  e For. Euler
!
	read(1,*)iflag
	write(2,*)'iflag=',iflag
	read(1,*)n,tau0
	write(2,*)'n=',n,'tau0=',tau0
	read(1,*) tend,dtau,teta0
	write(2,*)'tend=',tend,'dtau=',dtau,'teta0=',teta0
	read(1,*) rmu_max
	write(2,*)'rmu_max=',rmu_max
	read(1,*) rm
	write(2,*)'rm=',rm
	read(1,*) rik
	write(2,*)'rik=',rik
	read(1,*) rn
	write(2,*)'rn=',rn
	read(1,*) rka
	write(2,*)'rka=',rka
	read(1,*) p
	write(2,*)'p=',p
	read(1,*) ri0
	write(2,*)'ri0=',ri0

!	stop
!
!	initial values
!
	time=tau0
	do i=1,n
		tp(i)=teta0
	enddo
!	teta=teta0
!	time=tau0
	k=0 
	write(*,*) '     Table of results'
	write(*,*)'-----------------------------'
	write(*,*)'  Passo Nr     time       T(i)    '
	write(*,*)'-----------------------------'
	write(*,*)k,time,(tp(l),l=1,n)
!	write(2,*)time,tp(1)
!	write(3,*)time,teta
	write(4,*)time/86400.d0,tp(1)
	ht2=dtau/10
!
!	beginning of time loop
!
 50	k=k+1 
	tendi=time+dtau
      	write(*,*)'-------------time=',tendi
	if(iflag.eq.0) then	! RK-adaptativo
     call odeint(tp,n,time,tendi,1.e-6,ht2,1.e-20,id1,id2,nelmax,fcn,rkqc)
	endif
	if(iflag.eq.1) then	 ! RK-passo fixo
		call  rk4ord(tp,n,time,dtau,fcn,nelmax)
	endif
	if(iflag.eq.2) then	  ! Forward Euler
	    call fore(n,fcn,time,tp,tendi,nelmax)
	endif
	write(*,*)k,tendi/86400.d0,(tp(l),l=1,n)
	write(4,*)tendi/86400.d0,tp(1)
     	if (tendi.lt.tend) then
			time=tendi
			goto 50 
		endif
	close(2)
	close(4)

chamada = systemqq('notepad out-read.txt')       ! dados lidos
chamada = systemqq('notepad out-C.txt')       ! listagem da variavel
 ! chamada = systemqq('notepad outteo.txt')       ! listagem resultado teorico


chamada = systemqq('wgnuplot dados.gnu')  ! gráfico	geral
	stop
	end
!**************************************
	subroutine fcn(n,t,fi,f,nelmax)
	dimension fi(nelmax),f(nelmax)
	
	common /param/ ri0,p,rka,rmu_max,rn,rik,rm

		ri=0.4*ri0*exp(-p*rka*fi(1))  ! *375.5) ! 40% (photo-active radiation)
	!	rmu=rmu_max*ri**rn/(ri**rn+rik**rn)
           rmu=1.75e-5*ri**rn/(46.65d0**rn+ri**rn)
		f(1)=fi(1)*(rmu-rm)
	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 
!------------------------------------------------------------------ 





