subroutine convec(uf,vf,uv,vh,pf,nx,ny)
*     **********************************************************
*
*     nonsplit second-order upwind scheme for the linear advection
*     in the backfacing step channel:
*
*          phi_t + u phi_x + v phi_y = 0
*
*     velocity field given at cell centres and cell edges.
*
*     **********************************************************
      integer i,j,nx,ny,nxm,nym,mxym
      include 'gridsz'
      double precision dth,slimit
      double precision uf(nx,ny),vf(nx,ny),
     +                 uv(0:nx,0:ny),vh(0:nx,0:ny),
     +                 pf(0:nx+1,0:ny+1)
      double precision px(nxm,nym),py(nxm,nym),
     +                 pv(0:nxm,0:nym),ph(0:nxm,0:nym)
      common/timspc/dth
*     construct 1st derivatives with limiter in the x (normal) direction
*     and upwind in the y (transverse) direction
      do 5 j=1,ny
      do 5 i=1,nx
         px(i,j) = slimit(pf(i,j)-pf(i-1,j),pf(i+1,j)-pf(i,j))
         if (vf(i,j) .gt. 0.D0) then
             py(i,j) = pf(i,j) - pf(i,j-1)
         else
             py(i,j) = pf(i,j+1) - pf(i,j)
         end if
 5    continue
*     approximate vertical edge values.
      do 10 j=1,ny
      do 10 i=1,nx-1
         if (uv(i,j) .gt. 0.D0) then
             pv(i,j) = pf(i,j) + 0.5D0*( (1.D0 -
     +                 dth*uf(i,j) )*px(i,j) -
     +                 dth*vf(i,j)  *py(i,j) )
         else
             pv(i,j) = pf(i+1,j) - 0.5D0*( (1.D0 +
     +                 dth*uf(i+1,j) )*px(i+1,j) +
     +                 dth*vf(i+1,j)  *py(i+1,j) )
         end if
 10   continue
      do 22 j=1,ny
         pv(0,j) = pv(1,j)
         pv(nx,j) = pv(nx-1,j)
 22   continue
*
*     construct 1st derivatives with limiter in the y (normal) direction
*     and upwind in the x (transverse) direction
      do 25 j=1,ny
      do 25 i=1,ny
         py(i,j) = slimit(pf(i,j)-pf(i,j-1),pf(i,j+1)-pf(i,j))
         if (uf(i,j) .gt. 0.D0) then
             px(i,j) = pf(i,j) - pf(i-1,j)
         else
             px(i,j) = pf(i+1,j) - pf(i,j)
         end if
 25   continue
      do 30 j=1,ny-1
      do 30 i=1,nx
         if (vh(i,j) .gt. 0.D0) then
             ph(i,j) = pf(i,j) + 0.5D+00*( (1.0D+00 -
     +                 dth*vf(i,j) )*py(i,j) -
     +                 dth*uf(i,j)  *px(i,j) )
         else
             ph(i,j) = pf(i,j+1) - 0.5D+00*( (1.0D+00 +
     +                 dth*vf(i,j+1) )*py(i,j+1) +
     +                 dth*uf(i,j+1)  *px(i,j+1) )
         end if
 30   continue
      do 42 i=1,nx
         ph(i,0) = ph(i,1)
         ph(i,ny) = ph(i,ny-1)
 42   continue
*
*     calculate fluxes and update pf.
      do 50 j=1,ny
      do 50 i=1,nx
         pf(i,j) = pf(i,j) + dth*(
     +             uv(i-1,j)*pv(i-1,j) - uv(i,j)*pv(i,j) +
     +             vh(i,j-1)*ph(i,j-1) - vh(i,j)*ph(i,j) )
 50   continue
*     reflection boundary condition.
      call reflec(pf,nx,ny)
      end