************************************************************************
*
*  NOTE:  Every startup routine must initialize dx, dy and dz.  Stagger
*  routines for interpolation may be used freely, derivative routines
*  only after calling init_derivs().
*
      subroutine start (case, init)
*
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      common /ccase/casec
      character*32 casec
*
      character *(*) case
      logical init
      real p1, p2
      integer ip1
*-----------------------------------------------------------------------
*
      casec = case

      write (*,*) 'start: case=',case
      if (case(1: 6) .eq. 'excite') then
        call excite (init)
      else
        print *,'WARNING: Unknown start parameter:',case
        nx = mx
        ny = my
        nz = mz
      endif
*
      end
************************************************************************
      subroutine excite (init)
      logical init
*
*  Initialize values for stratification: polytropic index poly
*
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      include 'stagger.ifc'
*
      real kx, ky, kz
*-----------------------------------------------------------------------
*
      nx = mx
      ny = my
      nz = mz

      dx = (x4 - x1) / (ub-lb)
      dy = widthy / my
      dz = widthz / mz
      x = x1  + (findgen(mx)-(lb-1))*dx

      if (tcool .eq. 0.0) then
        call profx (x,x2,x3,dnu1,dnu2,0.0,0.0,0.0,fcool)
      else
        call profx (x,x2,x3,dnu1,dnu2,-1./tcool,0.0,-1./tcool,fcool)
      endif

      if (eetop .le. 0.0) eetop = x(lb)
      if (eebot .le. 0.0) eebot = x(ub)

      poly = 1./(gamma-1.)
      grav = (poly+1.)/poly

      print *,'eetop,eebot,grav =',eetop,eebot,grav

      if (.not. init) return

      do i=1,mx
        ee(i,:,:) = x(i)
        r(i,:,:) = x(i)**poly
      enddo
      e = r*ee

      px = 0.
      py = 0.
      pz = 0.

      kx=pi/(x4-x1)
      ky=2.*pi/my
      kz=2.*pi/mz
      rand=12345.6789
      x0 = x1+0.5*dx
      do k=1,mz
        zk=kz*k
        do j=1,my
          yk=ky*j
          planf=0.
          do l=1,4
            planf=planf+cos(yk*l+rand*l)*cos(zk*l+rand*2*l)/l
          enddo
          do i=1,mx
            xk=kx*(x(i)-x0)
            px(i,j,k) = sin(xk)*planf
          enddo
        enddo
      enddo
        cs=sqrt((gamma*(gamma-1.))*ee)
      px = r*cs*drive*px/maxval(px)

      ppbot = (gamma-1.)*e(ub,1,1)
      ssbot = ppbot/r(ub,1,1)**gamma

      end
*-----------------------------------------------------------------------
      subroutine profx (x,x1,x2,dnu1,dnu2,visc1,visc2,visc3,tmp)
      include 'cparam.inc'
      dimension tmp(mx), x(mx)
*
*  piecewice parabolic profiles
*
      do 10 i=1,mx
        if (x(i).lt.x1-dnu1) then
          tmp(i)=visc1
        elseif (x(i).lt.x1+dnu1) then
          tmp(i)=visc1+(visc2-visc1)*prof(x(i),x1,dnu1)
        elseif (x(i).lt.x2-dnu2) then
          tmp(i)=visc2
        elseif (x(i).lt.x2+dnu2) then
          tmp(i)=visc2+(visc3-visc2)*prof(x(i),x2,dnu2)
        else
          tmp(i)=visc3
        endif
10    continue
      end
*-----------------------------------------------------------------------
      function prof(x,x1,d1)
*
*  Piecewice 3rd order parabolic profiles
*
      xi=(x-x1)/d1
      prof=.5+.25*xi*(3.-xi**2)
      end