       SUBROUTINE adaptative_tau(tphys,tau,u,v,cc,p,div,nx,ny)
      include 'undefined.h'
      integer nx,ny,i,j,dfetch
      DOUBLE PRECISION u(nx,ny),v(nx,ny),cc(nx,ny),p(nx,ny),div(nx,ny)
      DOUBLE PRECISION mineps,maxeps,epsilon,velmax,fmodmax1,tphys
     * ,findninf,tau,h,normmean,hi2

      h=1.0d0/(nx-2)
      hi2 = 1.d0/(h*h)
      do j=2,ny-1
         do i=2,nx-1
            div(i,j) = ((i-1)*u(i+1,j)-(i-2)*u(i,j)+(i-1.5d0)*(v(i,j+1)
     *           - v(i,j)))*hi2
         enddo
      enddo

      epsilon = normmean(div,nx,ny)
      velmax = findninf(u,v,nx,ny)
      i=dfetch("mineps",mineps)
      i=dfetch("maxeps",maxeps)

      if ((epsilon .gt. mineps) .or. (velmax .gt. 0.23)) then
         call rescale(u,v,p,0.5d0,nx,ny)
         tau = tau*0.5d0
         write(6,*)"tphys,newtau =",tphys,tau         
      else if ((epsilon .lt. maxeps) .and. (velmax .lt. 0.09d0)) then
         call rescale(u,v,p,2.0d0,nx,ny)
         tau = tau*2.0d0
         write(6,*)"tphys,newtau =",tphys,tau
      endif

      return
      end

      SUBROUTINE rescale(u,v,p,tscale,nx,ny)
      include 'undefined.h'
      integer nx,ny,i,j
      DOUBLE PRECISION u(nx,ny),v(nx,ny),p(nx,ny),tscale,tscale2

      tscale2 = tscale**2
      do j=1,ny
         do i=1,nx
            u(i,j)=tscale*u(i,j)
            v(i,j)=tscale*v(i,j)
            p(i,j)=tscale2*p(i,j)
         enddo
      enddo

      return
      end

