c     ============================================
      subroutine b4step1(maxmx,mbc,mx,meqn,q,
     &            xlower,dx,t,dt,maux,aux)
c     ============================================
c
c     # called from claw1 before each call to step1.
c     # used to print out shock location, h, u, dh each time step.
c
c     
      implicit double precision (a-h,o-z)
      dimension q(1-mbc:maxmx+mbc, meqn)
      dimension aux(1-mbc:maxmx+mbc, *)

      integer midpoint, local_mthbc(2), step
      double precision x0, x1, t1, midx, frac
      logical regrid
      integer icount, iold, inewleft, inewright
      double precision qiold, qioldm1, qioldp1
      double precision leftslope, rightslope, myslope
      double precision npoly, nstar
      logical treset
      common /ceos/ grav, gamma, gamma1, sK, npoly, nstar, Crho, ieos
      common /comxt/ dtcom,dxcom,tcom
      common /comxt/ treset


c     if (xlower.gt.-7.d0) then
c	 write(6,*) 'Stop in b4step1'
c        stop
c	 endif

c     # check for NANs in solution:
      call check4nans(maxmx,meqn,mbc,mx,q,t,1)


c     # Check for negative depth:
      do i=1-mbc,mx+mbc
         if (q(i,1) .lt. 0.d0) then
            write(6,*) 'Resetting q(i,1) to zero, i = ',i
            write(6,*) '   q(i,1) = ',q(i,1)
            q(i,1) = 0.d0
            stop
            endif
         enddo

      regrid = .true.
      frac = 0.1d0
      x0 = -10.d0
      x1 = -10.d0
      t1 = 4.d0
      step = 30

      startx = x0 + t/t1*(x1-x0)
      startpoint = int((startx-xlower)/dx)
      if (startpoint .gt. mx) startpoint = mx
      if (startpoint .lt. 1+3+step) startpoint = 1+3+step
      write(30,*) t,startx,startpoint

      dhmax = 0.d0
      do i=startpoint,mx
         dh =  dabs(q(i+1,1)- q(i,1))
	 if (dh.gt.dhmax) then
	    dhmax = dh
	    imax = i
         endif
      enddo

c     # Modification by RJL to look at where dv/dx is most negative:
      dvmax = 0.d0
      imax = mx
      do i=mx,1,-1
         dv =  dabs(q(i,2)/q(i,1)- q(i+1,2)/q(i+1,1))
	 if (dv.gt.dvmax) then
	    dvmax = dv
	    imax = i
         endif
      enddo

      xs = xlower + imax*dx
c     write(6,*) '+++ imax = ',imax, '   dvmax = ',dvmax
      hs = q(imax-3,1)
      h1 = (hs - q(imax-3-step,1)) / (step*dx)
      us = q(imax-3,2)/hs
      u1 = (us - q(imax-3-step,2)/q(imax-3-step,1)) / (step*dx)
      dh = q(imax-3,1) - q(imax+3,1)
         
c     # Changed to write t and xs at end only after t is reset.
c     write(29,291) t,xs,hs,us,dh,h1,u1
c     write(39,*) t,imax,xs,dvmax

 291  format(7e25.15)

c     Regridding
      if (regrid .and. xs .gt. frac*xlower) then
         if (xs.gt.1.d-20) then
            write(6,*) '***Aborting at t = ',t,
     &                 ' since xs is too large, xs =',  xs
            stop
            endif

         write(6,*) 'regridding at t, xs =', t, xs
         write(6,*) 'xlower from', xlower, 'to', 0.5d0*xlower
         xlower = 0.5d0*xlower
         dx = 0.5d0*dx
         dt = 0.5d0*dt
         dxcom = dx      !# used in pres_psi
         dtcom = dt

C$$$ccccccccccccccccccccccccccccccccccccc
c        do i=1,mx
c           write(501,*) i, q(i,1)
c           write(502,*) i, q(i,2)/q(i,1)
c        end do
c        write(501,*) imax
c        write(502,*) imax
C$$$cccccccccccccccccccccccccccccccccccccc         

c     Regrid the conserved quantities. 
c     Rewritten by RJL (12/09) to assume mx even, which simplifies things
c     from CG's original version of minmod slopes assuming mx odd.
c     q(i,:) is the cell average in cell i.
c     Say we have cells 1, 2, 3, 4, 5, 6 (mx is always even).
c     New cell 1 is the left half of old cell 4. 
c     New cell 6 is the right half of old cell 6. 

         mx2 = mx/2
         do m=1,3
            do iold=mx2+1, mx
               qiold = q(iold,m)
               
c              # one-sided slopes:
               slopeleft  = qiold - q(iold-1,m)
               sloperight = q(iold+1,m) - qiold

c              # minmod limiter:
               if (dabs(slopeleft).lt.dabs(sloperight)) then
                  slope = slopeleft
                else
                  slope = sloperight
                endif
               if (slopeleft*sloperight .le. 0.d0) then
                  slope = 0.d0
                endif

c              # Or, use zero slopes:
c              slope = 0.d0  

               ileft = 2*(iold - mx2 - 1) + 1
               iright = ileft + 1
               q(ileft,m) = qiold - 0.25d0*slope
               q(iright,m) = qiold + 0.25d0*slope
               enddo
            enddo


c     In front of the shock, nothing should be moving, but when the
c     shock gets very close to the beach, numerical error accumulates.
c     To avoid this, overwrite with the rest solution.

c     Is this still needed with wb method?  
c     Yes, since otherwise averaging densities, etc onto finer grid gives
c     states that are no longer exact equilibrium states.

         do i=imax+3,mx
 	    xc = xlower + (i-0.5d0)*dx
            q(i,1) = Crho*(-xc)**nstar
            q(i,2) = 0.d0
            pressure = grav/(nstar+1.d0)*Crho * (-xc)**(nstar+1.d0)
            q(i,3) = pressure / gamma1 

c           +++ old:
c           q(i,1) = (-grav*xc / (npoly*gamma*sK))**npoly
c           q(i,2) = 0.d0
c           q(i,3) = sK * q(i,1)**gamma / gamma1 +
c    &                   0.5d0*q(i,2)**2/q(i,1)
c           +++ debug:
c           write(51,501) i,xc,q(i,1),q(i,3)
c 501       format( i5,3e20.10)
         end do

c     The parameter mthbc has not been passed to us, but we know what it is:
         local_mthbc(1)=1
         local_mthbc(2)=3
         call bc1(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt,
     $        local_mthbc)

c        +++ debug:
c        do i=mx-2,mx+2
c           write(51,511) i,xc,q(i,1),q(i,2),q(i,3)
c 511       format( i5,4e16.6)
c          end do

C$$$ccccccccccccccccccccccccccccccccccccc
c        do i=1,mx
c           write(501,*) i, q(i,1)
c           write(502,*) i, q(i,2)/q(i,1)
c        end do
c        write(501,*) imax
c        write(502,*) imax
c        stop 'stop for debugging after first regridding'
C$$$cccccccccccccccccccccccccccccccccccccc         

      endif

c     # always reset to equilibrium soln near right boundary
c     # see if this works for n=1.5, nstar=3. case?
      do i=mx-20,mx
 	    xc = xlower + (i-0.5d0)*dx
            write(51,501) i,xc,q(i,1),q(i,3)
            q(i,1) = Crho*(-xc)**nstar
            q(i,2) = 0.d0
            pressure = grav/(nstar+1.d0)*Crho * (-xc)**(nstar+1.d0)
            q(i,3) = pressure / gamma1 
c           +++ debug:
            write(51,501) i,xc,q(i,1),q(i,3)
  501       format( i5,3e20.10)
         end do

c     xreset = -1.d-5  ! used for Euler
      xreset = -1.d-4
      if (xlower.gt.xreset .and. (.not. treset)) then
         write(6,*) '==== Resetting t to 0 at time t = ',t
         write(26,261) t
  261    format(e26.16,'   time tc reset to 0')
         t = 0.d0
         treset = .true.
         endif

      if (treset) then
         write(29,291) t,xs
        else
         write(28,291) t,xs
        endif


      return
      end
