subroutine timrhs(uf,vf,px,py,ud,vd,ux,vx,uy,vy,fx,fy,ab,nx,ny)
c*******************************************************************
c
c     evaluate the forcing terms for the time stepping
c
c          u + dt*( 1/2/rey grad^2 u - (u grad u) - grad p )
c                                                                 
c
c*******************************************************************
      integer i,j,nx,ny,nxm,nym,mxym
      include 'gridsz'
      double precision fx(nx,ny),fy(nx,ny),
     +                 uf(nx,ny),vf(nx,ny),
     +                 ud(nx,ny),vd(nx,ny),
     +                 px(nx,ny),py(nx,ny),
     +                 ux(0:nx,0:ny),vx(0:nx,0:ny),
     +                 uy(0:nx,0:ny),vy(0:nx,0:ny),
     +                 dt,dth,hdtrey,rey,hf,ufx,vfx,uav,vav,
     +                 ab,uby,vby,utop,ubot,fadb,fadt,roa,grav
      common/timest/dt
      common/timspc/dth
      common/reynod/rey
      common/uvboud/uby(nym),vby(nym),utop,ubot
      common/gravity/grav
      common/density/roa(0:nxm+1,0:nym+1)
      parameter (hf=0.5D0)
      hdtrey = hf*dt/rey
c     ........................................................
      do 80 j=1,ny
      do 80 i=1,nx
         uav = hf*( ux(i-1,j) + ux(i,j) )
         vav = hf*( vy(i,j-1) + vy(i,j) )
         ufx = uav*( ux(i,j) - ux(i-1,j) ) +
     +         vav*( uy(i,j) - uy(i,j-1) )
         vfx = uav*( vx(i,j) - vx(i-1,j) ) +
     +         vav*( vy(i,j) - vy(i,j-1) )
         fx(i,j) = uf(i,j) + (hdtrey*ud(i,j) - 
     +             dt*px(i,j))/roa(i,j) - dth*ufx
         fy(i,j) = vf(i,j) + (hdtrey*vd(i,j) - 
     +             dt*py(i,j))/roa(i,j) - dth*vfx -
     +             dt*grav
   80 continue

      do 90 j=1,ny
         fx(1,j) = fx(1,j) + 3.2D0*ab*uby(j)
         fx(nx,j) = fx(nx,j) + 3.2D0*ab*uby(j)
 90   continue
         fadb = 3.2D0*ab*ubot
         fadt = 3.2D0*ab*utop
      do 95 i=1,nx
         fx(i,1) = fx(i,1) + fadb
         fx(i,ny) = fx(i,ny) + fadt
 95   continue

      end