c
c
c =========================================================
      subroutine src1(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt)
c =========================================================
      implicit double precision (a-h,o-z)
      dimension q(1-mbc:maxmx+mbc, meqn)
c
c
c     # solve the diffusion equation q_t = q_{xx} using Crank-Nicolson
c     # The LAPACK tridiagonal solver dgtsv is used, which is in tridiag.f
c
c
c     # local storage:
      parameter (maxldb = 5000)
      dimension b(maxldb,1), d(maxldb), dl(maxldb), du(maxldb)
      common /comsrc/ eps

      if (mx .gt. maxldb) then
          write(6,*) 'ERROR:  increase maxldb in src1.f'
          endif

      ldb = maxldb
      nrhs = 1
      dtdx2 = eps * dt / (2.d0*dx*dx) 
c     
c     # set coefficients in tridiagonal matrix and RHS:
      do i=1,mx
        dl(i) = -dtdx2
        d(i) = 1.d0 + 2.d0*dtdx2
        du(i) = -dtdx2
        b(i,1) = q(i,1) + dtdx2 * (q(i-1,1) - 2.d0*q(i,1) + q(i+1,1))
        enddo
c

c     # Boundary conditions:  At inflow use q=1:
      q0 = 1.d0

c     # At the outflow boundary, estimate q(mx+1,1) at the end of the
c     # time step by quadratic extrapolation of the current data:
      qmx1 = 2.d0*q(mx,1) - q(mx-1,1) 

c     # modify RHS for boundary values:
      b(1,1) = b(1,1) + dtdx2*q0
      b(mx,1) = b(mx,1) + dtdx2*qmx1
c
c
c     # Instead of the BC's above, could try to specify 
c     # no-flux boundary conditions for diffusion step:
c     # Adjust matrix entries to use q(0,1)=q(1,1) and q(mx+1,1)=q(mx,1)
c     # at end of time step:  (comment out above modifications to b if
c     # this is used)
c     d(1) = d(1) - dtdx2
c     d(mx) = d(mx) - dtdx2
      
      
c
c     # solve the tridiagonal system:
      call dgtsv(mx,nrhs,dl,d,du,b,ldb,info)

      if (info .ne. 0) then
         write(6,*) 'ERROR in src1 from call to dgtsv... info = ',info
         stop
         endif

      do i=1,mx
         q(i,1) = b(i,1)
	 enddo
c
      return
      end