subroutine init

C Kelvin-Helmholtz Instability

C  top and bottom have opposite X velocity producing a shear layer
C  in the middle.  This contact is wiggled (a sine wave) in order to
C  excite the KH instability.  Density of the two fluids is slightly
C  different to make the instability visible.
C  2D Cartesian
C                                                          sep98 blondin
C=======================================================================
C
C GLOBALS
C
      include 'global.h'
      include 'zone.h'
C
C LOCALS
C
      integer i, j, k, amp
      real kappa
C
C Set up geometry and boundary conditions of grid
C
C Boundary condition flags : nleft, nright
C   = 0  :  reflecting boundary condition
C   = 1  :  inflow/outflow boundary condition
C   = 2  :  fixed inflow boundary condition
C   = 3  :  periodic
C Geometry flag : ngeom                         |  Cartesian:
C   = 0  :  planar                              |    gx = 0, gy = 0, gz = 0
C   = 1  :  cylindrical radial                  |  Cylindrical:
C   = 2  :  spherical   radial             3D= {     gx = 0, gy = 1, gz = 3
C   = 3  :  cylindrical angle                   |
C   = 4  :  spherical polar angle (theta)       |  Spherical:
C   = 5  :  spherical azimu angle (phi)         |    gx = 2, gy = 4, gz = 5

C Define the problem...

      ndim   = 2
      ngeomx = 0
      ngeomy = 0
      ngeomz = 0
      nleftx = 3
      nrightx= 3
      nlefty = 1
      nrighty= 1
      nleftz = 0
      nrightz= 0
      xmin   = 0.
      xmax   = 1.0
      ymin   = 0.
      ymax   = 1.0
      zmin   = 0.
      zmax   = 1.

      gam    = 5. / 3.
      time   = 0.0
      courant= 0.6

      pi   = 2.0 * asin(1.0)
      gamm = gam - 1.0

C If any dimension is angular, multiply coordinates by pi...

      if(ngeomy.ge.3) then
         ymin = ymin * pi
         ymax = ymax * pi
      endif
      if(ngeomz.ge.3) then
         zmin = zmin * pi
         zmax = zmax * pi
      endif

C Set up grid coordinates and parabolic coefficients

      call grid(imax,xmin,xmax,zxa,zxc,zdx,zparax)
      call grid(jmax,ymin,ymax,zya,zyc,zdy,zparay)
      call grid(kmax,zmin,zmax,zza,zzc,zdz,zparaz)

C======================================================================
C Set up parameters from the problem (Sod shock tube)

      pamb = 1.0 
      dtop = 0.9 
      dbot = 1.1 
      vel  = 0.1 
      amp  = 2
      kappa = 3.0

C=======================================================================
C Log parameters of problem in history file

      write (8,*) 
      write (8,*) 'Kelvin-Helmholtz instability'
      write (8,*) 
      write (8,*) 'Adiabatic index, gamma = ', gam
      write (8,*) 'Mach number = ',vel/sqrt(pamb*gam/dbot)
      write (8,*) 'Density ratio is ', dtop/dbot
      write (8,*) 'Perturbation amplitude = ',amp,' zones'
      write (8,*) 'Perturbation wavenumber = ',kappa

C initialize grid:

      kappa = kappa * (2.*pi/(xmax-xmin))
      do k = 1, kmax
        do i = 1, imax
          jmid = jmax/2 + amp*sin(kappa*zxa(i))
          do j = 1, jmid
            zro(i,j,k) = dbot
            zpr(i,j,k) = pamb
            zux(i,j,k) = +vel
            zuy(i,j,k) = 0.
            zuz(i,j,k) = 0.
          enddo
          do j = jmid+1, jmax
            zro(i,j,k) = dtop
            zpr(i,j,k) = pamb
            zux(i,j,k) = -vel
            zuy(i,j,k) = 0.
            zuz(i,j,k) = 0.
          enddo
        enddo
      enddo

C set value of cmin & cmax for normalizing density in movie images
      cmax = 1.3
      cmin = 0.7

C########################################################################
C Compute Courant-limited timestep
      ridt = 0.
      j = 1
      k = 1
      if(ndim .eq. 1) then
         do  i = 1, imax
            svel = sqrt(gam*zpr(i,j,k)/zro(i,j,k))/zdx(i)
            xvel = abs(zux(i,j,k)) / zdx(i)
            ridt = max(xvel,svel,ridt)
         enddo
      else if(ndim .eq. 2) then
         do j = 1, jmax
         do i = 1, imax
            widthy = zdy(j)
            if(ngeomy.gt.2) widthy = widthy*(zxa(i)+0.5*zdx(i))
            width  = min(zdx(i),widthy)
            svel = sqrt(gam*zpr(i,j,k)/zro(i,j,k))/width
            xvel = abs(zux(i,j,k)) / zdx(i)
            yvel = abs(zuy(i,j,k)) / widthy
            ridt = max(xvel,yvel,svel,ridt)
         enddo
         enddo
      else if(ndim .eq. 3) then 
         do k = 1, kmax
         do j = 1, jmax
         do i = 1, imax
            widthy = zdy(j)
            widthz = zdz(k)
            if(ngeomy.gt.2) widthy = widthy*(zxa(i)+0.5*zdx(i))
            if(ngeomz.gt.2) widthz = widthz*(zxa(i)+0.5*zdx(i))
            if(ngeomz.eq.5) widthz = widthz*sin(zya(j)+.5*zdy(j))
            width  = min(zdx(i),widthy,widthz)
            svel = sqrt(gam*zpr(i,j,k)/zro(i,j,k))/width
            xvel = abs(zux(i,j,k)) / zdx(i)
            yvel = abs(zuy(i,j,k)) / widthy
            zvel = abs(zuz(i,j,k)) / widthz
            ridt = max(xvel,yvel,zvel,svel,ridt)
         enddo
         enddo
         enddo
      endif

      dt = courant / ridt

      return
      end