c===========================================================
c     Given input x-values x_j, j = 1 , n, returns gaussian 
c     profile, g_j and first two derivatives g'_j and 
c     g''_j.  Gaussian profile defined by
c
c     g(x;a,x0,del) = a * exp( - ((x-x0)/del)**2 )
c
c     Returns 0.0 for all three profiles if all inputs are 
c     not valid.
c===========================================================
      subroutine dvgaussian(g,dg,ddg,x,n,amp,x0,del)

         implicit        none

         integer         n
         real*8          g(n),     dg(n),    ddg(n),    x(n)
         real*8          amp,      x0,       del

         real*8          xmx0,     delm2,    mdelm2
         integer         j

         if( del .eq. 0.0d0  .or.  n .lt. 1 ) then
c-----------------------------------------------------------
c           Invalid input.
c-----------------------------------------------------------
            do j = 1 , n
               g(j)   = 0.0d0
               dg(j)  = 0.0d0
               ddg(j) = 0.0d0
            end do
         else
            delm2  = 1.0d0 / (del ** 2)
            mdelm2 = -2.0d0 * delm2
            do j = 1 , n
               xmx0 = x(j) - x0
               g(j)   = amp * exp(-delm2 * xmx0**2)
               dg(j)  = mdelm2 * xmx0 * g(j)
               ddg(j) = mdelm2 * (1.0d0 + mdelm2 * xmx0**2) 
     &                         * g(j)
            end do
         end if

         return

      end

c===========================================================
c     Output to standard out for subsequent plotting via 
c     gnuplot.
c===========================================================
      subroutine gnuout(u,x,nx,t,stride)
         implicit       none

         integer        nx,           stride
         real*8         u(nx),        x(nx),        t

         integer        j
         
c-----------------------------------------------------------
c        Basic gnuplot 'splot' format for parametric 
c        plotting is:
c 
c        <x>    <y>    <f(x,y)>
c-----------------------------------------------------------
         do j = 1 , nx , stride
            write(*,*) t, x(j), u(j)
         end do
c-----------------------------------------------------------
c        Empty line separates rows of data.
c-----------------------------------------------------------
         write(*,*)

         return

      end