module VARGAS
				  
use INPUT
use MALHA

use MSFLIB


contains


subroutine principal_vargas
 
!     main (gabinet simulation)
!     Feb 28, 1998
 
      implicit real *8 (a-h,o-z)
 
!     nelmax must reserve 5 positions for the 5 outside wall temperatures
 
      external fcn,rkqc
      common /prop/ rhoair,gban,cvair,rkair,alfa,pr,cpair,rkwall
      common /geom/ thick
      common /temp/ tpoute,tpoutw,tpouts,tpoutn,tpoutt,  &
           tground
      common /const/ ht2
      common /rad/ tpair,alf,epsi,sig,riavg,rhoal,esp,cal
      common /hwl/ a,b,c
      common /heat/ ct,cq
 
!     insulation thermal conductivity of the 5 walls
!     1-east; 2-west; 3-south; 4-north and 5-top
 
      common /sol2/ rkiso(5)  
      common /rad1/ iwe,iww,iwn,iws
      common /mesh/ nel,nlx,nly,nlz,neleqp
      common /pert/ epsilon

      integer :: iconta    ! contador do conta
 
      pw1=(pw(2)-pw(1))/(tw(2)-tw(1))
      pwn=(pw(n2)-pw(n2-1))/(tw(n2)-tw(n2-1))
      call spline(tw,pw,n2,pw1,pwn,pw2)
 
      call gmesh

!	initial data
      tpout(1)=tpoute
      tpout(2)=tpoutw
      tpout(3)=tpouts
      tpout(4)=tpoutn
      tpout(5)=tpoutt
      iout=1

      if ( irest == 0 ) then
         do i=1,nel+5
            if(i.le.nel) then
	          tp(i)=tp0
      	    else
      	   	  tp(nel+iout)=tpout(iout)
      	   	  iout=iout+1
            endif
         enddo
      endif
 
!     compute vapor pressure
 
      temp=tp0-273.15d0
      call splint(tw,pw,pw2,n2,temp,pg)
      pv=hum0*pg


!     steady-state solution

      if (isteady.eq.1) then

 
         call mnewt(maxit,tp,nto,tol2,tol3,xl,nto)
         do l=1,nel
            if (itypel(l).eq.0) then
               temp=tp(l)-273.15d0
               call splint(tw,pw,pw2,n2,temp,pg)
               hum(l)=pv/pg
            else
               hum(l)=0.d0
            endif
         enddo

 
!     transient solution

      else

 
      time=0.d0
      nt=0
      do 100 j=1,ntime
      	do 110 i=1,nel+5
      	  t0(i)=tp(i)
 110    continue
        tend=time+htime
        nt=nt+1   
        call odeint(tp,nto,time,tend,tol,ht2,ht3,id1,id2   &
     			,xl,nto,fcn,rkqc)
        do l=1,nel
           if (itypel(l).eq.0) then
              temp=tp(l)-273.15d0
              call splint(tw,pw,pw2,n2,temp,pg)
              hum(l)=pv/pg
           else
              hum(l)=0.d0
           endif
        enddo
        t0norm=rnorm2(nto,t0,nto)
        do 120 l=1,nel+5
          t0(l)=tp(l)-t0(l)
 120    continue
        dtnorm=rnorm2(nto,t0,nto)

        if ( mod(j,fres) == 0 ) then
           open (9,file='itera')
           do iconta = 1, conta
             read(9,*)
           end do
           write(9,901) j, dtnorm/htime, timef()-tcpui
           901 format(/,i5,5x,2(1pe10.2,2x))
           close(9)
           conta = conta + 2
        end if

        if(dtnorm.lt.tol1*htime) then
        goto 300
        endif
        time=tend
 100  continue

      open (9,file='itera')
      do iconta = 1, conta
         read(9,*)
      end do
      write(9,*) 'O estado permanente não foi atingido!'
      close(9)


      endif


 300  continue
2000  format(5x,i10,2x,f14.7)


      end subroutine principal_vargas
       
!----------------------------------------------------------------

      function rnorm2(n,x,nd)
 
!	compute euclidean norm of a vector
 
      implicit real *8 (a-h,o-z)
      dimension x(nd)
      sum=0.d0
      do i=1,n
      	sum=sum+x(i)*x(i)
      enddo
      aux=sqrt(sum)
      rnorm2=aux
      return
      end function rnorm2

!------------------------------------------------------------

      function rninf(n,x,nd)
 
!	compute infinity norm of a vector
 
      implicit real *8 (a-h,o-z)
      dimension x(nd)
      rmax=0.d0
      do i=1,n
      	if(abs(x(i)).gt.rmax) then
           rmax=abs(x(i))
        endif
      enddo
      rninf=rmax
      return
      end function rninf

!---------------------------------------------------------------

      subroutine fore(n,fcn,time,fi,tend,xl,nelmax) 
  
      implicit real *8 (a-h,o-z)
      dimension fi(nelmax),xl(nelmax,8,3)

      common /const/ ht2 
      external fcn 

      real*8,dimension(:),allocatable :: f

      allocate ( f(nelmax) )

      k=0 
 50   k=k+1 
      time=min(time+ht2,tend) 
      call fcn(n,fi,f,xl,nelmax) 
      do 100 i=1,n 
         fi(i)=fi(i)+ht2*f(i) 
 100  continue 
      if (time.lt.tend) goto 50 

      deallocate ( f )

      return 
      end subroutine fore

!------------------------------------------------------------------ 

      subroutine usrfun(u,n,np,g,a,xl,nelmax)
 
!	compute the values of g for a given u and the jacobian
 
      implicit real *8 (a-h,o-z)
      dimension u(nelmax),g(np),a(np,np)
      dimension xl(nelmax,8,3)

      common /pert/ epsilon

      real*8,dimension(:),allocatable :: ge

      allocate ( ge(nelmax) )

      call fcn(n,u,g,xl,nelmax)
 
!     build the jacobian numerically
 
      do j=1,n
         u(j)=u(j)+epsilon
         call fcn(n,u,ge,xl,nelmax)
         u(j)=u(j)-epsilon
         do i=1,n
            a(i,j)=(ge(i)-g(i))/epsilon
         enddo
      enddo

      deallocate ( ge )

      return
      end subroutine usrfun

!-----------------------------------------------------

      SUBROUTINE mnewt(ntrial,x,n,tolx,tolf,xl,nelmax)

      implicit real *8 (a-h,o-z)
      dimension x(nelmax),xl(nelmax,8,3)

      integer :: iconta    ! contador do conta
      real*8 rnorma

      real*8,dimension(:),allocatable :: indx,fvec,p
      real*8,dimension(:,:),allocatable :: fjac

      allocate ( indx(n), fvec(n), p(n), fjac(n,n) )

      do 14  k=1,ntrial
        call usrfun(x,n,n,fvec,fjac,xl,nelmax)
        errf=0.
        do 11 i=1,n
          errf=errf+abs(fvec(i))
11      continue

        if ( mod(k,fres) == 0 ) then
           open (9,file='itera')
           do iconta = 1, conta
             read(9,*)
           end do
           rnorma = errf
           write(9,901) k, rnorma, timef()-tcpui
           close(9)
           conta = conta + 2
		end if

        if(errf.le.tolf) then
          deallocate ( indx, fvec, p, fjac )
          return
        end if

        do 12 i=1,n
          p(i)=-fvec(i)
12      continue

        call ludcmp(fjac,n,N,indx,d)
        call lubksb(fjac,n,N,indx,p)
        errx=0.
        do 13 i=1,n
          errx=errx+abs(p(i))
          x(i)=x(i)+p(i)
13      continue
        enorm=rnorm2(n,p,n) 

        if ( mod(k,fres) == 0 ) then
           if ( errf > errx .or. errf > enorm) then
              conta = conta - 2
              open (9,file='itera')
              do iconta = 1, conta
                read(9,*)
              end do
              rnorma = errx
              if ( enorm < errx ) rnorma = enorm
              write(9,901) k, rnorma, timef()-tcpui
              901 format(/,i5,5x,2(1pe10.2,2x))
              close(9)
              conta = conta + 2
           end if
        end if

        if(errx.le.tolx.or.enorm.le.tolx) then
          deallocate ( indx, fvec, p, fjac )
           return
        endif
14    continue

      deallocate ( indx, fvec, p, fjac )

      return
      END subroutine mnewt


end module VARGAS