c
c      =======================================================
       subroutine src2(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,
     &               dx,dy,q,maux,aux,t,dt)
c      =======================================================

       implicit double precision (a-h,o-z)

       dimension aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,maux)
       dimension q(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,meqn)

c
       integer ndim
       parameter (ndim = 2*300*150)
       integer info(4), idid
       double precision rtol, atol
       double precision y(ndim),dydt(ndim),work(8+5*ndim)
       integer neqn, i,j,l, m

       external f
      
       common /grid/ nx,ny
       
       common /rkcdid/  nfe,nsteps,naccpt,nrejct,nfesig,maxm

       
      if(ndim.lt.mx*my*2) then
         write(6,*) "need to change ndim in src2"
         write(6,*) "ndim = ", mx*my*2
         stop
      endif 

      do i=1, ndim
        y(i) = 0.d0
        dydt(i) = 0.d0
      enddo

c     specify parameters for rkc:
c  
      nx = mx
      ny = my
c
      tende = t+dt
      info(1) = 0  !# compute a solution at tend only
      info(2) = 0  !# 0 - rkc computes the spectral radius 
      info(3) = 0  !# the Jacobian may not be constant
      info(4) = 0  !# ATOL is a scalar
c
      rtol = 1.d-8
      atol = rtol
      idid = 0
      neqn = 2*nx*ny
    
      time = t

c     set y-vector, i.e. initial values of solution 
     
       
      do m=1, 2
        do j=1, my
          do i=1, mx
            l = (m-1)*mx*my+i+(j-1)*mx
            y(l) = q(i,j,m)  
          enddo
        enddo
      enddo
      
 10   continue
      
      call rkc(neqn,f,y,time,tende,rtol,atol,info,work,idid) 
c     # if not done yet, take another step
      if(idid.eq.2) goto 10 
c
      if(idid.gt.2) then
        write(6,*) "Error in rkc"
        write(6,*) 'Failed at t= ',t,' with idid = ',idid
        stop
      endif 
c     # output information for rkc solver
      write(6,*) "rkc - method:"
      write(6,*) "there were", nsteps, "steps (", nrejct," rejected)."     
      write(6,*) "the maximum number of stages used was", maxm
c
      do j=1, my
        do i=1, mx
          l = i+(j-1)*mx
          q(i,j,1) = y(l)
          q(i,j,2) = y(mx*my+l)
        enddo
      enddo 
c    
c
      return
      end
c
c=====================================================================
      subroutine f(neqn,time,y,dydt)
      implicit double precision (a-h,o-z)
      integer neqn
      dimension y(*),dydt(*)
      double precision t, dx, dy
      integer i,j,k,l,ibc,jbc
      double precision u(0:300+1,0:150+1,1:2)
      integer nx,ny
      common /grid/ nx,ny 
      common /dmetric/ dmetric(-1:300+2,-1:150+2,5)
      common /reaction/ D,del,alf,beta,gam,tau1,tau2

      if(nx.gt.300.or.ny.gt.150) then
        write(6,*) "need to increase size of u-array of subroutine f" 
        write(6,*)  "in src2"
        stop
      endif 
      dx = 4.d0/nx
      dy = 2.d0/ny
      do m = 1, 2
        do j=1, ny
          do i=1, nx
            u(i,j,m) = y((m-1)*nx*ny+i+(j-1)*nx)
          enddo
        enddo
      enddo  
c     set boundary condition for sphere on top and bottom
      do m=1, 2
          do i=1, nx
            i2 = nx+1-i
            u(i,0,m) = u(i2,1,m)
            u(i,ny+1,m) = u(i2,ny,m)  
          enddo 
      enddo
c     periodic boundary condition in xc-direction
      do m=1, 2
         do j=0, ny+1
            u(0,j,m) = u(nx,j,m)
            u(nx+1,j,m) = u(1,j,m)
         enddo  
      enddo  
c
      do m=1, 2
        if(m.eq.1) Dr = D*del
        if(m.eq.2) Dr = del
        do j=1, ny
          do i=1, nx
            l = (m-1)*nx*ny+i+(j-1)*nx
            q1 = (u(i,j+1,m)+u(i+1,j+1,m)-u(i,j-1,m)-u(i+1,j-1,m))/4.d0
            q2 = (u(i,j+1,m)+u(i-1,j+1,m)-u(i,j-1,m)-u(i-1,j-1,m))/4.d0
            q3 = (u(i+1,j,m)+u(i+1,j+1,m)-u(i-1,j,m)-u(i-1,j+1,m))/4.d0
            q4 = (u(i+1,j-1,m)+u(i+1,j,m)-u(i-1,j-1,m)-u(i-1,j,m))/4.d0
c 
            dydt(l)=Dr*((dsqrt(0.5d0*(dmetric(i,j,1)+dmetric(i+1,j,1)))
     &   *(0.5d0*(dmetric(i,j,2)+dmetric(i+1,j,2))*(u(i+1,j,m)
     &   -u(i,j,m))/dx + 0.5d0*(dmetric(i,j,3)+dmetric(i+1,j,3))*q1/dy)
     &   -dsqrt(0.5d0*(dmetric(i-1,j,1)+dmetric(i,j,1)))*
     &   (0.5d0*(dmetric(i-1,j,2)+dmetric(i,j,2))*(u(i,j,m)-u(i-1,j,m))
     &   /dx+0.5d0*(dmetric(i-1,j,3)+dmetric(i,j,3))*q2/dy))/
     &   (dsqrt(dmetric(i,j,1))*dx)+(dsqrt(0.5d0*(dmetric(i,j,1)
     &   +dmetric(i,j+1,1)))*(0.5d0*(dmetric(i,j+1,4)+
     &   dmetric(i,j,4))*q3/dx + 0.5d0*(dmetric(i,j,5)+dmetric(i,j+1,5))
     &   *(u(i,j+1,m)-u(i,j,m))/dy)-dsqrt(0.5*(dmetric(i,j,1)
     &   +dmetric(i,j-1,1)))*(0.5d0*(dmetric(i,j,4)+dmetric(i,j-1,4))
     &   *q4/dx + 0.5d0*(dmetric(i,j,5)+dmetric(i,j-1,5))*(u(i,j,m)
     &   -u(i,j-1,m))/dy)) /(dsqrt(dmetric(i,j,1))*dy) )
c
            if(m.eq.1) then
             dydt(l) = dydt(l)+alf*u(i,j,1)*(1.d0-tau1*(u(i,j,2)**2.d0))
     &                        +u(i,j,2)*(1.d0-tau2*u(i,j,1))         
            else
              dydt(l) = dydt(l)+beta*u(i,j,2)*
     &        (1.d0+(alf*tau1*u(i,j,1)*u(i,j,2))/beta) 
     &        + u(i,j,1)*(gam+tau2*u(i,j,2))
            endif
          enddo
        enddo
      enddo 

      return
      end

      double precision function spcrad(neqn,time,y)
c--------------------------
c  This is a dummy routine.
c--------------------------
      integer          neqn
      double precision time,y(*)

      spcrad = 0d0  !# for dummy routine
      return
      end

