subroutine init

C Sedov-Taylor Blastwave

C   the central 4 zones are set to a very high pressure to
C   drive a blastwave into the surrounding medium.
C   2D cylindrical (R-Z)
C                                                         sep98 blondin
C=======================================================================
C
C GLOBALS
C
      include 'global.h'
      include 'zone.h'
C
C LOCALS
C
      integer i, j, k
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 = 1
      ngeomy = 0
      ngeomz = 0
      nleftx = 0
      nrightx= 1
      nlefty = 0
      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
      pst  = 1.0e+09 
      damb = gam 

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

      write (8,*) 
      write (8,*) 'Sedov-Taylor Blastwave'
      write (8,*) 
      write (8,*) 'Adiabatic index, gamma = ', gam
      write (8,*) 'Blast pressure = ',pst
      write (8,*) 'Outer x location = ',xmax
      write (8,*) 'Inner x location = ',xmin
      write (8,*) 'Outer y location = ',ymax
      write (8,*) 'Inner y location = ',ymin
      write (8,*) 'Outer z location = ',zmax
      write (8,*) 'Inner z location = ',zmin

C initialize grid:

      do k = 1, kmax
        do j = 1, jmax
          do i = 1, imax
            zro(i,j,k) = damb
            zpr(i,j,k) = pamb
            zux(i,j,k) = 0.0
            zuy(i,j,k) = 0.
            zuz(i,j,k) = 0.
          enddo
        enddo
      enddo

      do j = 1, 2
        do i = 1, 2
          zpr(i,j,1) = pst
        enddo
      enddo

C set value of cmin & cmax for normalizing density in movie images
      cmax = damb*4.5
      cmin = 0.0

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