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

      subroutine rhs

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

c---------------------------------------------------------------------
c   compute the right hand sides
c---------------------------------------------------------------------

      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, k, m
      double precision  q
      double precision  tmp, utmp(6,isiz3), rtmp(5,isiz3)
      double precision  u21, u31, u41
      double precision  u21i, u31i, u41i, u51i
      double precision  u21j, u31j, u41j, u51j
      double precision  u21k, u31k, u41k, u51k
      double precision  u21im1, u31im1, u41im1, u51im1
      double precision  u21jm1, u31jm1, u41jm1, u51jm1
      double precision  u21km1, u31km1, u41km1, u51km1


      if (timeron) call timer_start(t_rhs)
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               do m = 1, 5
                  rsd(m,i,j,k) = - frct(m,i,j,k)
               end do
               tmp = 1.0d+00 / u(1,i,j,k)
               rho_i(i,j,k) = tmp
               qs(i,j,k) = 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
     >                         + u(3,i,j,k) * u(3,i,j,k)
     >                         + u(4,i,j,k) * u(4,i,j,k) )
     >                      * tmp
            end do
         end do
      end do

      if (timeron) call timer_start(t_rhsx)
c---------------------------------------------------------------------
c   xi-direction flux differences
c---------------------------------------------------------------------

      do k = 2, nz - 1
         do j = jst, jend
            do i = 1, nx
               flux(1,i) = u(2,i,j,k)
               u21 = u(2,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,i) = u(2,i,j,k) * u21 + c2 * 
     >                        ( u(5,i,j,k) - q )
               flux(3,i) = u(3,i,j,k) * u21
               flux(4,i) = u(4,i,j,k) * u21
               flux(5,i) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
            end do

            do i = ist, iend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                 - tx2 * ( flux(m,i+1) - flux(m,i-1) )
               end do
            end do

            do i = ist, nx
               tmp = rho_i(i,j,k)

               u21i = tmp * u(2,i,j,k)
               u31i = tmp * u(3,i,j,k)
               u41i = tmp * u(4,i,j,k)
               u51i = tmp * u(5,i,j,k)

               tmp = rho_i(i-1,j,k)

               u21im1 = tmp * u(2,i-1,j,k)
               u31im1 = tmp * u(3,i-1,j,k)
               u41im1 = tmp * u(4,i-1,j,k)
               u51im1 = tmp * u(5,i-1,j,k)

               flux(2,i) = (4.0d+00/3.0d+00) * tx3 * (u21i-u21im1)
               flux(3,i) = tx3 * ( u31i - u31im1 )
               flux(4,i) = tx3 * ( u41i - u41im1 )
               flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
     >                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tx3 * ( u21i**2 - u21im1**2 )
     >              + c1 * c5 * tx3 * ( u51i - u51im1 )
            end do

            do i = ist, iend
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dx1 * tx1 * (            u(1,i-1,j,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i+1,j,k) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) )
     >              + dx2 * tx1 * (            u(2,i-1,j,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i+1,j,k) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) )
     >              + dx3 * tx1 * (            u(3,i-1,j,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i+1,j,k) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) )
     >              + dx4 * tx1 * (            u(4,i-1,j,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i+1,j,k) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) )
     >              + dx5 * tx1 * (            u(5,i-1,j,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i+1,j,k) )
            end do

c---------------------------------------------------------------------
c   Fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,2,j,k) = rsd(m,2,j,k)
     >           - dssp * ( + 5.0d+00 * u(m,2,j,k)
     >                      - 4.0d+00 * u(m,3,j,k)
     >                      +           u(m,4,j,k) )
               rsd(m,3,j,k) = rsd(m,3,j,k)
     >           - dssp * ( - 4.0d+00 * u(m,2,j,k)
     >                      + 6.0d+00 * u(m,3,j,k)
     >                      - 4.0d+00 * u(m,4,j,k)
     >                      +           u(m,5,j,k) )
            end do

            do i = 4, nx - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i-2,j,k)
     >                        - 4.0d+00 * u(m,i-1,j,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i+1,j,k)
     >                        +           u(m,i+2,j,k) )
               end do
            end do


            do m = 1, 5
               rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
     >           - dssp * (             u(m,nx-4,j,k)
     >                      - 4.0d+00 * u(m,nx-3,j,k)
     >                      + 6.0d+00 * u(m,nx-2,j,k)
     >                      - 4.0d+00 * u(m,nx-1,j,k)  )
               rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
     >           - dssp * (             u(m,nx-3,j,k)
     >                      - 4.0d+00 * u(m,nx-2,j,k)
     >                      + 5.0d+00 * u(m,nx-1,j,k) )
            end do

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

      if (timeron) call timer_start(t_rhsy)
c---------------------------------------------------------------------
c   eta-direction flux differences
c---------------------------------------------------------------------
      do k = 2, nz - 1
         do i = ist, iend
            do j = 1, ny
               flux(1,j) = u(3,i,j,k)
               u31 = u(3,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,j) = u(2,i,j,k) * u31 
               flux(3,j) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q)
               flux(4,j) = u(4,i,j,k) * u31
               flux(5,j) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
            end do

            do j = jst, jend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                   - ty2 * ( flux(m,j+1) - flux(m,j-1) )
               end do
            end do

            do j = jst, ny
               tmp = rho_i(i,j,k)

               u21j = tmp * u(2,i,j,k)
               u31j = tmp * u(3,i,j,k)
               u41j = tmp * u(4,i,j,k)
               u51j = tmp * u(5,i,j,k)

               tmp = rho_i(i,j-1,k)
               u21jm1 = tmp * u(2,i,j-1,k)
               u31jm1 = tmp * u(3,i,j-1,k)
               u41jm1 = tmp * u(4,i,j-1,k)
               u51jm1 = tmp * u(5,i,j-1,k)

               flux(2,j) = ty3 * ( u21j - u21jm1 )
               flux(3,j) = (4.0d+00/3.0d+00) * ty3 * (u31j-u31jm1)
               flux(4,j) = ty3 * ( u41j - u41jm1 )
               flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
     >                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * ty3 * ( u31j**2 - u31jm1**2 )
     >              + c1 * c5 * ty3 * ( u51j - u51jm1 )
            end do

            do j = jst, jend

               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dy1 * ty1 * (            u(1,i,j-1,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j+1,k) )

               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) )
     >              + dy2 * ty1 * (            u(2,i,j-1,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j+1,k) )

               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) )
     >              + dy3 * ty1 * (            u(3,i,j-1,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j+1,k) )

               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) )
     >              + dy4 * ty1 * (            u(4,i,j-1,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j+1,k) )

               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) )
     >              + dy5 * ty1 * (            u(5,i,j-1,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j+1,k) )

            end do
         end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
         do i = ist, iend
            do m = 1, 5
               rsd(m,i,2,k) = rsd(m,i,2,k)
     >           - dssp * ( + 5.0d+00 * u(m,i,2,k)
     >                      - 4.0d+00 * u(m,i,3,k)
     >                      +           u(m,i,4,k) )
               rsd(m,i,3,k) = rsd(m,i,3,k)
     >           - dssp * ( - 4.0d+00 * u(m,i,2,k)
     >                      + 6.0d+00 * u(m,i,3,k)
     >                      - 4.0d+00 * u(m,i,4,k)
     >                      +           u(m,i,5,k) )
            end do
         end do

         do j = 4, ny - 3
            do i = ist, iend
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j-2,k)
     >                        - 4.0d+00 * u(m,i,j-1,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j+1,k)
     >                        +           u(m,i,j+2,k) )
               end do
            end do
         end do

         do i = ist, iend
            do m = 1, 5
               rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
     >           - dssp * (             u(m,i,ny-4,k)
     >                      - 4.0d+00 * u(m,i,ny-3,k)
     >                      + 6.0d+00 * u(m,i,ny-2,k)
     >                      - 4.0d+00 * u(m,i,ny-1,k)  )
               rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
     >           - dssp * (             u(m,i,ny-3,k)
     >                      - 4.0d+00 * u(m,i,ny-2,k)
     >                      + 5.0d+00 * u(m,i,ny-1,k) )
            end do
         end do

      end do
      if (timeron) call timer_stop(t_rhsy)

      if (timeron) call timer_start(t_rhsz)
c---------------------------------------------------------------------
c   zeta-direction flux differences
c---------------------------------------------------------------------
      do j = jst, jend
         do i = ist, iend
            do k = 1, nz
               utmp(1,k) = u(1,i,j,k)
               utmp(2,k) = u(2,i,j,k)
               utmp(3,k) = u(3,i,j,k)
               utmp(4,k) = u(4,i,j,k)
               utmp(5,k) = u(5,i,j,k)
               utmp(6,k) = rho_i(i,j,k)
            end do
            do k = 1, nz
               flux(1,k) = utmp(4,k)
               u41 = utmp(4,k) * utmp(6,k)

               q = qs(i,j,k)

               flux(2,k) = utmp(2,k) * u41 
               flux(3,k) = utmp(3,k) * u41 
               flux(4,k) = utmp(4,k) * u41 + c2 * (utmp(5,k)-q)
               flux(5,k) = ( c1 * utmp(5,k) - c2 * q ) * u41
            end do

            do k = 2, nz - 1
               do m = 1, 5
                  rtmp(m,k) =  rsd(m,i,j,k)
     >                - tz2 * ( flux(m,k+1) - flux(m,k-1) )
               end do
            end do

            do k = 2, nz
               tmp = utmp(6,k)

               u21k = tmp * utmp(2,k)
               u31k = tmp * utmp(3,k)
               u41k = tmp * utmp(4,k)
               u51k = tmp * utmp(5,k)

               tmp = utmp(6,k-1)

               u21km1 = tmp * utmp(2,k-1)
               u31km1 = tmp * utmp(3,k-1)
               u41km1 = tmp * utmp(4,k-1)
               u51km1 = tmp * utmp(5,k-1)

               flux(2,k) = tz3 * ( u21k - u21km1 )
               flux(3,k) = tz3 * ( u31k - u31km1 )
               flux(4,k) = (4.0d+00/3.0d+00) * tz3 * (u41k-u41km1)
               flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
     >                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tz3 * ( u41k**2 - u41km1**2 )
     >              + c1 * c5 * tz3 * ( u51k - u51km1 )
            end do

            do k = 2, nz - 1
               rtmp(1,k) = rtmp(1,k)
     >              + dz1 * tz1 * (            utmp(1,k-1)
     >                             - 2.0d+00 * utmp(1,k)
     >                             +           utmp(1,k+1) )
               rtmp(2,k) = rtmp(2,k)
     >          + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) )
     >              + dz2 * tz1 * (            utmp(2,k-1)
     >                             - 2.0d+00 * utmp(2,k)
     >                             +           utmp(2,k+1) )
               rtmp(3,k) = rtmp(3,k)
     >          + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) )
     >              + dz3 * tz1 * (            utmp(3,k-1)
     >                             - 2.0d+00 * utmp(3,k)
     >                             +           utmp(3,k+1) )
               rtmp(4,k) = rtmp(4,k)
     >          + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) )
     >              + dz4 * tz1 * (            utmp(4,k-1)
     >                             - 2.0d+00 * utmp(4,k)
     >                             +           utmp(4,k+1) )
               rtmp(5,k) = rtmp(5,k)
     >          + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) )
     >              + dz5 * tz1 * (            utmp(5,k-1)
     >                             - 2.0d+00 * utmp(5,k)
     >                             +           utmp(5,k+1) )
            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,j,2) = rtmp(m,2)
     >           - dssp * ( + 5.0d+00 * utmp(m,2)
     >                      - 4.0d+00 * utmp(m,3)
     >                      +           utmp(m,4) )
               rsd(m,i,j,3) = rtmp(m,3)
     >           - dssp * ( - 4.0d+00 * utmp(m,2)
     >                      + 6.0d+00 * utmp(m,3)
     >                      - 4.0d+00 * utmp(m,4)
     >                      +           utmp(m,5) )
            end do

            do k = 4, nz - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rtmp(m,k)
     >              - dssp * (            utmp(m,k-2)
     >                        - 4.0d+00 * utmp(m,k-1)
     >                        + 6.0d+00 * utmp(m,k)
     >                        - 4.0d+00 * utmp(m,k+1)
     >                        +           utmp(m,k+2) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,j,nz-2) = rtmp(m,nz-2)
     >           - dssp * (             utmp(m,nz-4)
     >                      - 4.0d+00 * utmp(m,nz-3)
     >                      + 6.0d+00 * utmp(m,nz-2)
     >                      - 4.0d+00 * utmp(m,nz-1)  )
               rsd(m,i,j,nz-1) = rtmp(m,nz-1)
     >           - dssp * (             utmp(m,nz-3)
     >                      - 4.0d+00 * utmp(m,nz-2)
     >                      + 5.0d+00 * utmp(m,nz-1) )
            end do
         end do
      end do
      if (timeron) call timer_stop(t_rhsz)
      if (timeron) call timer_stop(t_rhs)

      return
      end