subroutine upwind(ul,ur,vl,vr,u,v)
c     ****************************************************************
c
c     Godunov upwinding
c
c     ****************************************************************
      double precision ul,ur,vl,vr,u,v
      if (ul .ge. 0.D+00 .and. ul+ur .ge. 0.D+00) then
       u = ul
      elseif ( ul .le. 0.D+00 .and. ur .ge. 0.D+00) then
       u = 0.D+00
      else
       u = ur
      endif
      if (u) 10,20,30
   10    v = vr
      return
   20    v = 0.5D+00*(vl+vr)
      return
   30    v = vl
      return
      end

      subroutine upwsec(u,ul,ur,v,vl,vr,du,dul,dur,dv,dvl,dvr,
     +                  duup,dvup)
c     ****************************************************************
c
c     second order upwinding of the derivative
c
c     ****************************************************************
      double precision u,ul,ur,v,vl,vr,du,dul,dur,dv,dvl,dvr,
     +                 duup,dvup,dth,duc
      common/timspc/dth
      if (u .ge. 0.D0) then
         duc = 0.5D0*( 1.D0 - u*dth )
         duup = u - ul + duc*( du - dul )
         dvup = v - vl + duc*( dv - dvl )
      else
         duc = -0.5D0*( 1.D0 + u*dth )
         duup = ur - u + duc*( dur - du )
         dvup = vr - v + duc*( dvr - dv )
      endif
      end

      function dm(dr,dl)
*     ********************************************************
*
*     limiter for the second-order hj upwind scheme
*
*     ********************************************************
      double precision dm,dr,dl
      if( dabs(dr) .le. dabs(dl) ) then
          dm = dr
      else
          dm = dl
      endif
      end

      function slimit(ul,ur)
c     ***********************************************************
c
c     limiter definition in the interior
c
c     ul,ur         ->         values on the left and right sides
c     slimit        <-         resolved value
c
c     ***********************************************************
      double precision ul,ur,dut,slimit
      if ( ul*ur .le. 0.D+00 ) then
         slimit = 0.D+00
      else
         dut = 0.5D0*( ul + ur )
         slimit = dsign( 1.D0, dut )*dmin1( 
     +            dabs(dut), 2.D0*dabs(ul), 2.D0*dabs(ur) )
      end if
c      slimit = 0.5D0*(ul+ur)
      end

      function blimit(ub,u1,u2)
c     ***********************************************************
c
c     limiter definition on the boundaries
c
c     ub,u1,u2      ->         values spaced by h/2,h and h
c     blimit        <-         resolved value
c
c     ***********************************************************
      double precision ub,u1,u2,blimit,umax,umin,ubt
      ubt = 0.25D0*( 3.D0*u1 - u2 + 2.D0*ub )
      umax = dmax1( u1, ub )
      umin = dmin1( u1, ub )
      if( ubt .gt. umax ) then
         blimit = 2.D0*( u1 - umax )
      else if( ubt .lt. umin ) then
         blimit = 2.D0*( u1 - umin )
      else
         blimit = 0.5D0*( u2 + u1 - 2.D0*ub )
      end if
      end