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

       subroutine  tzetar(c)

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

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

       include 'header.h'

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

c---------------------------------------------------------------------
c      treat only one cell                                             
c---------------------------------------------------------------------
       do    k = start(3,c), cell_size(3,c)-end(3,c)-1
          do    j = start(2,c), cell_size(2,c)-end(2,c)-1
             do    i = start(1,c), cell_size(1,c)-end(1,c)-1

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

                ac2u = ac*ac

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

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

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

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

             end do
          end do
       end do

       return
       end
