! $Id: courant.f,v 1.8 2000/06/09 13:08:21 aake Exp $
!***********************************************************************
      subroutine Courant
!
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      real cx, cy, cz
      integer ijk
      common /cdf/ sawtooth_factor

!-----------------------------------------------------------------------
!  HD/MHD wave Courant number
!-----------------------------------------------------------------------
      cx  = (1./dx)
      cy  = (1./dy)
      cz  = (1./dz)
      cs  = (gamma*(gamma-1.))*e/r
      cs  = sqrt(cs)
      scr1  = amax1(
     &         cx*(abs(ux) + cs),
     &         cy*(abs(uy) + cs),
     &         cz*(abs(uz) + cs))
      Cu = maxval(scr1(lb+1:ub,:,:))
      nflop = nflop + 5 + 13

!-----------------------------------------------------------------------
!  Diffusive Courant number
!-----------------------------------------------------------------------

!  The sawtooth factor is the amplitude of ddxup(ddxdn(f)) for
!  a sawtooth function with amplitude 1 (f=[+1,-1,...]).  For the
!  stagger_6th.f, it is equal to 6.17.  The quenching function 
!  is an additional factor 2. on top, for a total of 12.3.

      sawtooth_factor = 6.2
!!      sawtooth_factor = 12.3           ! quench()
 
!  The safety factor aims at marginal stability, where the diffusive
!  term is twice (2.) the fluctuation, and with a staggered sawtooth
!  in three (3.) dimensions as the worst case.

      cdtmax = 0.75
      safety = sawtooth_factor*3./2.*Cdtmax

!  According to this estimate, one should be able to run with
!  dt*nu = 1/(6.2*1.5) = 1/(9.3) = 0.11 (for dx=1). In practice 
!  this depends somewhat on the time stepping method. With rk3_2n, 
!  the limit is around (dt=0.4)*(nu=0.35) = 0.14.
 
      cx = cx**2
      cy = cy**2
      cz = cz**2

      scr1 = amax1( cx*vx, cy*vy, cz*vz)
      Cv  = maxval(scr1(lb+1:ub,:,:))
      Cv = Cv*safety
      nflop = nflop + 7

!-----------------------------------------------------------------------
!  New time step
!-----------------------------------------------------------------------
      write (*,'(1x,a,2f5.2,1pe10.2,1pe10.3,0pf7.0)') 
     & 'Courant: Cu, Cv,  dt, t, kz/s = ', 
     &  Cu*dt, Cv*dt, dt, t, mw*1e-3/cpu

      dtold = dt
      if (nstep .gt. 3) then
        if (tcool .gt. 0.0) then
          dt = amin1(abs(tcool)*1.5, dtold*1.4, Cdt/amax1(Cu, Cv))
        else
          dt = amin1(dtold*1.4, Cdt/amax1(Cu, Cv))
        endif
      endif

      end