
c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine  tzetar

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   block-diagonal matrix-vector multiplication                       
c---------------------------------------------------------------------

       include 'header.h'

       integer i, j, k
       double precision  t1, t2, t3, ac, xvel, yvel, zvel, r1, r2, r3, 
     >                   r4, r5, btuz, ac2u, uzik1


       if (timeron) call timer_start(t_tzetar)
!$omp parallel do default(shared)
!$omp& private(i,j,k,t1,t2,t3,ac,xvel,yvel,zvel,r1,r2,r3, 
!$omp&              r4,r5,btuz,ac2u,uzik1)
       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2

                xvel = us(i,j,k)
                yvel = vs(i,j,k)
                zvel = ws(i,j,k)
                ac   = speed(i,j,k)

                ac2u = ac*ac

                r1 = rhs(1,i,j,k)
                r2 = rhs(2,i,j,k)
                r3 = rhs(3,i,j,k)
                r4 = rhs(4,i,j,k)
                r5 = rhs(5,i,j,k)      

                uzik1 = u(1,i,j,k)
                btuz  = bt * uzik1

                t1 = btuz/ac * (r4 + r5)
                t2 = r3 + t1
                t3 = btuz * (r4 - r5)

                rhs(1,i,j,k) = t2
                rhs(2,i,j,k) = -uzik1*r2 + xvel*t2
                rhs(3,i,j,k) =  uzik1*r1 + yvel*t2
                rhs(4,i,j,k) =  zvel*t2  + t3
                rhs(5,i,j,k) =  uzik1*(-xvel*r2 + yvel*r1) + 
     >                    qs(i,j,k)*t2 + c2iv*ac2u*t1 + zvel*t3

             end do
          end do
       end do
       if (timeron) call timer_stop(t_tzetar)

       return
       end