subroutine dtcon
C
C set the timestep using various constraints.
C NB - additional constraints may be required depending on the problem!
C
C
      include 'global.h'
      include 'zone.h'
      integer i, j
      real  ridt, dtx, dt3, xvel, yvel, zvel, svel
     &     ,width, widthy, widthz
C      
C------------------------------------------------------------------------
C
C   Hydro constraint on timestep.  Use R*d(theta) if y geometry is angular
      ridt = 0.
C
      if(ndim .eq. 1) then
         j = 1
         k = 1
         do 10 i = 1, imax
            svel = sqrt(gam*zp(i,j,k)/zrho(i,j,k))/zdx(i)
            xvel = abs(zux(i,j,k)) / zdx(i)
            ridt = max(xvel,svel,ridt)
 10      continue
      else if(ndim .eq. 2) then
         k = 1
         do 20 j = 1, jmax
         do 20 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*zp(i,j,k)/zrho(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)
 20      continue
      else if(ndim .eq. 3) then 
         do 30 k = 1, kmax
         do 30 j = 1, jmax
         do 30 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*zp(i,j,k)/zrho(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)
 30      continue
      endif
C
C global time constraint for given courant parameter
      dtx = courant / ridt
C                                        
C limiting constraint on rate of increase of dt
      dt3 = 1.1 * dt
C                                        
C use smallest required timestep                                        
      dt = min( dt3, dtx )
C
C if timestep becomes too small, stop the program!
      if (time/dt .gt. 1.e6) then
         write(*,*) 'Timestep has become too small: dt = ',dt
         write(*,*) '                             time = ',time
         stop
      endif
C
      return
      end