SUBROUTINE compute_tstep(mx,my,dx,dy,x,y,aux,frside,itstep,dt_min)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,*) :: aux,frside

  umax = 0
  vmax = 0
  DO i = 1-mbc,mx+mbc-1
     DO j = 1-mbc,my+mbc-1
        SELECT CASE(itstep)
        CASE(1)
           !! Darcy velocities divided by kappa (\nu_{ij})
           u1 = MAX(ABS(aux(i,j,1)),ABS(aux(i+1,j,1)))
           v1 = MAX(ABS(aux(i,j,2)),ABS(aux(i,j+1,2)))
           umax = MAX(u1/aux(i,j,3),umax)
           vmax = MAX(v1/aux(i,j,3),vmax)
        CASE(2)
           !! Pointwise velocities only (\nu)
           fri =  MAX(frside(i,j,1),1e-5)
           frip = MAX(frside(i+1,j,1),1e-5)
           u1 = MAX(ABS(aux(i,j,1))/fri,ABS(aux(i+1,j,1))/frip)

           frj =  MAX(frside(i,j,2),1e-5)
           frjp = MAX(frside(i,j+1,2),1e-5)
           v1 = MAX(ABS(aux(i,j,2))/frj,ABS(aux(i,j+1,2))/frjp)

           umax = MAX(u1,umax)
           vmax = MAX(v1,vmax)
        END SELECT
     END DO
  END DO

  cfl = 0.9
  IF (umax /= 0 .AND. vmax /= 0) THEN
     dt_min = cfl*MIN(dx/umax, dy/vmax)
  ELSE IF (umax == 0.0 .AND. vmax /= 0.0) THEN
     dt_min = cfl*dy/vmax
  ELSE IF (vmax == 0.0 .AND. umax /= 0.0) THEN
     dt_min = cfl*dx/umax
  ENDIF

  RETURN
END SUBROUTINE compute_tstep
