subroutine init

C Standing Strong Shock - the Achilles heel of PPM
C  1D Cartesian
C                                                          sep98 blondin
C=======================================================================

C GLOBALS
      include 'global.h'
      include 'zone.h'

C LOCALS
      integer i, j, k
      real mach,xshock,pamb,damb

C Set up geometry and boundary conditions of grid

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

      ! Define the problem...

      ndim   = 1
      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

      ! If only 1D, then force jmax=kmax=1

      if(ndim.eq.1) jmax = 1
      if(ndim.lt.3) kmax = 1

      ! 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

      ! 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======================================================================
      ! Set up parameters from the problem (Sod shock tube)

      data mach   / 10.  /  !  mach number of incoming shock
      data pamb   / 1.0  /  !  ambient gas pressure
      data damb   / 1.0  /  !  ambient gas density

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

      ishock = imax / 2
      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))
      upwind = mach * sqrt(pamb/damb/gam)

      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 - upwind
        enddo
        do i = ishock + 1, imax
         zro(i,j,k) = damb
         zpr(i,j,k) = pamb
         zuy(i,j,k) = 0.
         zuz(i,j,k) = 0.
         zux(i,j,k) = -upwind
        enddo
       enddo
      enddo

      ! set value of cmin & cmax for normalizing density in movie images
      cmax = 6.0
      cmin = 0.0

      ! Log parameters of problem in history file
      write (8,*) 
      write (8,*) 'Strong Standing Shock'
      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

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