subroutine init

C Hawley-Zabusky Shock-tube 

C a planar shock of Mach number 1.2 impinges an oblique
C discontinuity separating unshocked gas with a density
C contrast of dratio (=3) across the discontinuity.
C In order to see the long term evolution (including the
C roll-up of the KH instability), imax should be made 
C much bigger (ie, make the box long and skinny).
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
      real mach,xshock,pamb,damb,cangle,xcon,dratio
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 = 1
      nrightx= 1
      nlefty = 0
      nrighty= 0
      nleftz = 0
      nrightz= 0
      xmin   = 0.
      ymin   = 0.
      ymax   = 1.0
      xmax   = (ymax-ymin)*real(imax)/real(jmax)+xmin
      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)

      data mach   / 1.2  /  !  mach number of incoming shock
      data ishock / 5    /  !  initial location on grid of incoming shock
      data pamb   / 1.0  /  !  ambient gas pressure
      data damb   / 1.0  /  !  ambient gas density
      data cangle / 30.  /  !  angle of contact discontinuity in degrees
      data icon   / 10   /  !  initial location where contact discontinuit
      data dratio / 3.0  /  !  density ratio across discontinuity

C=======================================================================
C initialize grid:

      xcon   = zxa(icon)
      gas1   = (gam+1.)/(gam-1.)
      gas2   = 2.0*gam/(gam+1.)
      pratio = 1.+gas2*(mach*mach-1.)
      pshock = pamb * pratio
      dshock = damb * (1.+gas1*pratio)/(gas1+pratio)
      ushock = sqrt(pamb/damb/gam)*(pratio-1.0)
     &         *sqrt(gas2/(pratio+1.0/gas1))
      angle  = cangle * pi/180.
      tang   = tan(angle)

      do k = 1, kmax
       do j = 1, jmax 
        do i = 1, ishock      
         zuy(i,j,k) = 0.
         zuz(i,j,k) = 0.
         zro(i,j,k) = dshock
         zpr(i,j,k) = pshock
         zux(i,j,k) = ushock
        enddo
        zyap = zya(j) + zdy(j)
        xlocat1 = xcon + (zya(j)-zya(1))/tang
        xlocat2 = xcon + (zyap  -zya(1))/tang
        do i = ishock+1, imax  
        zxap = zxa(i) + zdx(i)
         if (zxa(i) .gt. xlocat2) then
           density = dratio * damb
         else if (zxap .lt. xlocat1) then
           density = damb
         else if ((zxa(i).le.xlocat1).and.(zxap.le.xlocat2) ) then
           frac = .5*(zxap-xlocat1)**2*tang/(zdx(i)*zdy(j))
           density = (1-frac)*damb + frac*damb*dratio
         else if ((zxa(i).le.xlocat1).and.(zxap.gt.xlocat2)) then
           frac = ((zxap-xlocat2) + .5*(xlocat2-xlocat1))/zdx(i)
           density = (1-frac)*damb + frac*damb*dratio
         else if ((zxa(i).ge.xlocat1).and.(zxap.gt.xlocat2)) then
           ypt = ((zxa(i)-xcon)*tang + zya(1))
           frac = .5*(zyap-ypt)*(xlocat2-zxa(i))/(zdx(i)*zdy(j))
           density = frac*damb + (1-frac)*damb*dratio
         else if ((zxa(i).ge.xlocat1).and.(zxap.lt.xlocat2)) then
           ypt = ((zxa(i)-xcon)*tang + zya(1))
           ypt2= ((zxap-xcon)*tang + zya(1))
           frac = ((ypt-zya(j)) +.5*(ypt2-ypt))/zdy(j)
           density = (1-frac)*damb + frac*damb*dratio
         endif
         zro(i,j,k) = density
         zpr(i,j,k) = pamb
         zuy(i,j,k) = 0.
         zuz(i,j,k) = 0.
         zux(i,j,k) = 0.
        enddo
       enddo
      enddo

C set value of cmin & cmax for normalizing density in movie images
      cmax = 5.0
      cmin = 0.0

C Log parameters of problem in history file
      write (8,*) 
      write (8,*) 'Hawley-Zabusky shock tube'
      write (8,*) 
      write (8,*) 'Shock mach number ', mach
      write (8,*) 'Adiabatic index, gamma = ', gam
      write (8,*) 'Ambient shock tube pressure is ', pamb
      write (8,*) 'Ambient shock tube density  is ', damb
      write (8,*) 'Post - shock pressure is ', pshock
      write (8,*) 'Post - shock density is ', dshock
      write (8,*) 'Post - shock velocity is ', ushock
      write (8,*) 'Density ratio at contact is ', dratio
      write (8,*) 'Angle of contact is ', cangle

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