c * * * * * * * * * * * * * * * * * * * * * * * * *
c    This example is slightly modified version of example presented 
c    by Assyr Abdulle and described in the book of 
c    G.Wanner and E.Hairer Solving Ordinary Differential Equations II
c * * * * * * * * * * * * * * * * * * * * * * * * *
c
c    This driver shows how to use DUMKA3. It solves a
c    system of ODEs resulting from the 2-dimensional space 
c    discretization of the Brusselator equations (u=u(x,y,t),v=v(x,y,t)):
c     
c    u_t=1+u^2*v-4.4*u+0.1*(u_{xx}+u_{yy})+f(x,y,t)
c    v_t=3.4*u-u^2*v+0.1*(v_{xx}+v_{yy})     for t>=0, 0<= x <= 1, 0<= y <= 1
c
c    with initial conditions 
c
c    u(x,y,0)=22*y*(1-y)^{3/2}  v(x,y,0)=27*x*(1-x)^{3/2}
c
c    and periodic boundary conditions
c
c    u(x+1,y,t)=u(x,y,t),  v(x,y+1,t)=v(x,y,t).
c
c    The function f is defined by (inhomogeneity)
c
c    f(x,y,t)=5 if (x-0.3)^2+(y-0.6)^2<= 0.1^2 and t>=1.1
c            =0 else
c
c    We discretize the space variables with
c    x_i=i/(N+1), y_i=i/(N+1) for i=0,1,...,N,
c    with N=128. We obtain a system of 32768
c    equations. The spectral radius of the Jacobian 
c    can be estimated with the Gershgorin theorem   
c    (13200 is an estimation for it). Thus we 
c    provide an external function RHO, giving 
c    the spectral radius of the Jacobian matrix.
c    As output point we choose t_out=1.5 
c
c--------------------------------------------------------    
      include 'dumka3.f'
      implicit double precision (a-h,o-z)
      parameter (nsd=128,neqn=nsd*nsd*2)
      dimension y(nsd,nsd,2),z(nsd,nsd,2),z1(nsd,nsd,2),z2(nsd,nsd,2),
     & z3(nsd,nsd,2),z4(nsd,nsd,2),atol(1),rtol(1)
      integer iwork(12)
      external fbrus
c --- common parameters for the problem -----
      common/trans/alf,ns,nssq,nsm1sq
c ----- file for solution -----
       rewind 8
c ----- dimensions -----
      ns=nsd
      nssq=ns*ns
      n=2*nssq 
      alf=1.0d-1
      alph=alf
c--------------------------------------------------------
c     Initialize iwork: 
c      iwork(1)=0  DUMKA3 (!try to!) calculate spectral radius internally
c                    (if iwork(1)=1 - rho is used to determine spectral radius)
c      iwork(2)=0  The Jacobian is not constant
c                    (if iwork(2)=1 - rho is called only once)
c      iwork(3)=0  Calculate problem until t=tend 
c                    (if iwork(3)=1 - dumka exits after the first step)
c      iwork(4)=0  Atol and rtol are scalars.
c--------------------------------------------------------
      iwork(1)=0
      iwork(2)=0
      iwork(3)=0
      iwork(4)=0
c ----- initial and end point of integration -----
      t=0.0d0
      tend=1.5d0
c ----- initial values -----
      
      ans=ns
      do j=1,ns  
        yy=(j-1)/ans      
        do i=1,ns          
          xx=(i-1)/ans
          y(j,i,1)=22.d0*yy*(1.d0-yy)**(1.5d0)
          y(j,i,2)=27.d0*xx*(1.d0-xx)**(1.5d0)
        end do
      end do
c ----- required tolerance ----- 
      
      tol = 0.01d0
      rtol(1)=tol
      atol(1)=tol
c ----- initial step size -----
      h=1.0d-4 
c ----- integration -----
      write(6,*) 'integration of the 2-dim brusselator problem'  
      call dumka3(neqn,t,tend,h,fbrus,y,z,z1,z2,z3,
     &             z4,atol,rtol,iwork) 
      open(8,file='solDumkaU') 
      open(9,file='solDumkaV')
c ----- print solution -----
      do j=1,ns
        yy=(j-1)/ans 
        do i=1,ns
          xx=(i-1)/ans
          write (8,300) xx,yy,y(j,i,1)
          write (9,300) xx,yy,y(j,i,2)
        end do
      end do
      close (8)
      close(9)
   
c ----- print statistics -----
      write(6,*) 'Solution is tabulated in file sol.out'
      write(6,*) 'Max estimation of the spectral radius=',iwork(11)
      write(6,*) 'Min estimation of the spectral radius=',iwork(12)
      write(6,*) 'Max number of stages used=',iwork(10)
      write(6,*) 'Number of f eval. for the spectr. radius=',iwork(9)
      write (6,91) iwork(5),iwork(6),iwork(7),iwork(8)
 91   format(' Number of f evaluations=',i5,' steps=',i4,
     &        ' accpt=',i4,' rejct=',i3) 
300   format(1f22.16,' ',1f22.16,' ',1f22.16)
c--------------------------------------------------------
c     End of main program
c--------------------------------------------------------
      end      
c--------------------------------------------------------
c     The subroutine RHO gives an estimation of the spectral 
c     radius of the Jacobian matrix of the problem. This
c     is a bound for the whole interval and thus RHO is called
c     once.
c--------------------------------------------------------

      double precision function rho(neqn,t,y)
      implicit double precision (a-h,o-z)
      dimension y(ns,ns,2)
      common/trans/alf,ns,nssq,nsm1sq
        rho = 8.0d0*nssq*alf + 2.d0 
      return
      end 
c--------------------------------------------------------
c     The subroutine fun compute the value of f(x,y) and
c     has to be declared as external.
c--------------------------------------------------------
      subroutine fbrus(neqn,x,y,f)
c ----- brusselator with diffusion in 2 dim. space -----
      implicit double precision (a-h,o-z)
      dimension y(ns,ns,2),f(ns,ns,2)
      common/trans/alf,ns,nssq,nsm1sq
c ----- constants for inhomogenity -----
      ans=ns
      radsq=0.1d0**2
      if(x.ge.1.1d0)then
        bet=5.00d0
      else
        bet=0.00d0
      end if
c ----- big loop -----
      do j=1,ns
      do i=1,ns
c ----- left neighbour -----
         if(i.eq.1)then
            uleft=y(j,ns,1)
            vleft=y(j,ns,2)
         else
            uleft=y(j,i-1,1)
            vleft=y(j,i-1,2)
         end if
c ----- right neighbour -----
         if(i.eq.ns)then
            uright=y(j,1,1)
            vright=y(j,1,2)
         else
            uright=y(j,i+1,1)
            vright=y(j,i+1,2)
         end if
c ----- lower neighbour -----
         if(j.eq.1)then
            ulow=y(ns,i,1)
            vlow=y(ns,i,2)
         else
            ulow=y(j-1,i,1)
            vlow=y(j-1,i,2)
         end if
c ----- upper neighbour -----
         if(j.eq.ns)then
            uup=y(1,i,1)
            vup=y(1,i,2)
         else
            uup=y(j+1,i,1)
            vup=y(j+1,i,2)
         end if
c ----- the derivative -----
         uij=y(j,i,1)
         vij=y(j,i,2)
         f(j,i,1)=1.d0+uij*uij*vij-4.4d0*uij
     &        +alf*nssq*(uleft+uright+ulow+uup-4.d0*uij)
         f(j,i,2)=3.4d0*uij - uij*uij*vij
     &        +alf*nssq*(vleft+vright+vlow+vup-4.d0*vij)
c ----- inhomogenity -----

        yy=(j-1)/ans 
        xx=(i-1)/ans
        if(((xx-0.3d0)**2+(yy-0.6d0)**2).le.radsq)then
          f(j,i,1)=f(j,i,1)+bet
        end if
      end do
      end do
      return
      end