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

      subroutine jacu(k)

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

c---------------------------------------------------------------------
c   compute the upper triangular part of the jacobian matrix
c---------------------------------------------------------------------


      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      integer k

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j
      double precision  r43
      double precision  c1345
      double precision  c34
      double precision  tmp1, tmp2, tmp3


      r43 = ( 4.0d+00 / 3.0d+00 )
      c1345 = c1 * c3 * c4 * c5
      c34 = c3 * c4

!$omp do schedule(static)
         do j = jend, jst, -1
            do i = iend, ist, -1

c---------------------------------------------------------------------
c   form the block daigonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               du(1,1,i,j) =  1.0d+00
     >                       + dt * 2.0d+00 * (   tx1 * dx1
     >                                          + ty1 * dy1
     >                                          + tz1 * dz1 )
               du(1,2,i,j) =  0.0d+00
               du(1,3,i,j) =  0.0d+00
               du(1,4,i,j) =  0.0d+00
               du(1,5,i,j) =  0.0d+00

               du(2,1,i,j) =  dt * 2.0d+00
     >           * ( - tx1 * r43 - ty1 - tz1 )
     >           * ( c34 * tmp2 * u(2,i,j,k) )
               du(2,2,i,j) =  1.0d+00
     >          + dt * 2.0d+00 * c34 * tmp1 
     >          * (  tx1 * r43 + ty1 + tz1 )
     >          + dt * 2.0d+00 * (   tx1 * dx2
     >                             + ty1 * dy2
     >                             + tz1 * dz2  )
               du(2,3,i,j) = 0.0d+00
               du(2,4,i,j) = 0.0d+00
               du(2,5,i,j) = 0.0d+00

               du(3,1,i,j) = dt * 2.0d+00
     >           * ( - tx1 - ty1 * r43 - tz1 )
     >           * ( c34 * tmp2 * u(3,i,j,k) )
               du(3,2,i,j) = 0.0d+00
               du(3,3,i,j) = 1.0d+00
     >         + dt * 2.0d+00 * c34 * tmp1
     >              * (  tx1 + ty1 * r43 + tz1 )
     >         + dt * 2.0d+00 * (  tx1 * dx3
     >                           + ty1 * dy3
     >                           + tz1 * dz3 )
               du(3,4,i,j) = 0.0d+00
               du(3,5,i,j) = 0.0d+00

               du(4,1,i,j) = dt * 2.0d+00
     >           * ( - tx1 - ty1 - tz1 * r43 )
     >           * ( c34 * tmp2 * u(4,i,j,k) )
               du(4,2,i,j) = 0.0d+00
               du(4,3,i,j) = 0.0d+00
               du(4,4,i,j) = 1.0d+00
     >         + dt * 2.0d+00 * c34 * tmp1
     >              * (  tx1 + ty1 + tz1 * r43 )
     >         + dt * 2.0d+00 * (  tx1 * dx4
     >                           + ty1 * dy4
     >                           + tz1 * dz4 )
               du(4,5,i,j) = 0.0d+00

               du(5,1,i,j) = -dt * 2.0d+00
     >  * ( ( ( tx1 * ( r43*c34 - c1345 )
     >     + ty1 * ( c34 - c1345 )
     >     + tz1 * ( c34 - c1345 ) ) * ( u(2,i,j,k) ** 2 )
     >   + ( tx1 * ( c34 - c1345 )
     >     + ty1 * ( r43*c34 - c1345 )
     >     + tz1 * ( c34 - c1345 ) ) * ( u(3,i,j,k) ** 2 )
     >   + ( tx1 * ( c34 - c1345 )
     >     + ty1 * ( c34 - c1345 )
     >     + tz1 * ( r43*c34 - c1345 ) ) * ( u(4,i,j,k) ** 2 )
     >      ) * tmp3
     >   + ( tx1 + ty1 + tz1 ) * c1345 * tmp2 * u(5,i,j,k) )

               du(5,2,i,j) = dt * 2.0d+00
     > * ( tx1 * ( r43*c34 - c1345 )
     >   + ty1 * (     c34 - c1345 )
     >   + tz1 * (     c34 - c1345 ) ) * tmp2 * u(2,i,j,k)
               du(5,3,i,j) = dt * 2.0d+00
     > * ( tx1 * ( c34 - c1345 )
     >   + ty1 * ( r43*c34 -c1345 )
     >   + tz1 * ( c34 - c1345 ) ) * tmp2 * u(3,i,j,k)
               du(5,4,i,j) = dt * 2.0d+00
     > * ( tx1 * ( c34 - c1345 )
     >   + ty1 * ( c34 - c1345 )
     >   + tz1 * ( r43*c34 - c1345 ) ) * tmp2 * u(4,i,j,k)
               du(5,5,i,j) = 1.0d+00
     >   + dt * 2.0d+00 * ( tx1 + ty1 + tz1 ) * c1345 * tmp1
     >   + dt * 2.0d+00 * (  tx1 * dx5
     >                    +  ty1 * dy5
     >                    +  tz1 * dz5 )

c---------------------------------------------------------------------
c   form the first block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i+1,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               au(1,1,i,j) = - dt * tx1 * dx1
               au(1,2,i,j) =   dt * tx2
               au(1,3,i,j) =   0.0d+00
               au(1,4,i,j) =   0.0d+00
               au(1,5,i,j) =   0.0d+00

               au(2,1,i,j) =  dt * tx2
     >          * ( - ( u(2,i+1,j,k) * tmp1 ) ** 2
     >     + c2 * qs(i+1,j,k) * tmp1 )
     >          - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i+1,j,k) )
               au(2,2,i,j) =  dt * tx2
     >          * ( ( 2.0d+00 - c2 ) * ( u(2,i+1,j,k) * tmp1 ) )
     >          - dt * tx1 * ( r43 * c34 * tmp1 )
     >          - dt * tx1 * dx2
               au(2,3,i,j) =  dt * tx2
     >              * ( - c2 * ( u(3,i+1,j,k) * tmp1 ) )
               au(2,4,i,j) =  dt * tx2
     >              * ( - c2 * ( u(4,i+1,j,k) * tmp1 ) )
               au(2,5,i,j) =  dt * tx2 * c2 

               au(3,1,i,j) =  dt * tx2
     >              * ( - ( u(2,i+1,j,k) * u(3,i+1,j,k) ) * tmp2 )
     >         - dt * tx1 * ( - c34 * tmp2 * u(3,i+1,j,k) )
               au(3,2,i,j) =  dt * tx2 * ( u(3,i+1,j,k) * tmp1 )
               au(3,3,i,j) =  dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
     >          - dt * tx1 * ( c34 * tmp1 )
     >          - dt * tx1 * dx3
               au(3,4,i,j) = 0.0d+00
               au(3,5,i,j) = 0.0d+00

               au(4,1,i,j) = dt * tx2
     >          * ( - ( u(2,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 )
     >          - dt * tx1 * ( - c34 * tmp2 * u(4,i+1,j,k) )
               au(4,2,i,j) = dt * tx2 * ( u(4,i+1,j,k) * tmp1 )
               au(4,3,i,j) = 0.0d+00
               au(4,4,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
     >          - dt * tx1 * ( c34 * tmp1 )
     >          - dt * tx1 * dx4
               au(4,5,i,j) = 0.0d+00

               au(5,1,i,j) = dt * tx2
     >          * ( ( c2 * 2.0d0 * qs(i+1,j,k)
     >              - c1 * u(5,i+1,j,k) )
     >          * ( u(2,i+1,j,k) * tmp2 ) )
     >          - dt * tx1
     >          * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i+1,j,k)**2 )
     >              - (     c34 - c1345 ) * tmp3 * ( u(3,i+1,j,k)**2 )
     >              - (     c34 - c1345 ) * tmp3 * ( u(4,i+1,j,k)**2 )
     >              - c1345 * tmp2 * u(5,i+1,j,k) )
               au(5,2,i,j) = dt * tx2
     >          * ( c1 * ( u(5,i+1,j,k) * tmp1 )
     >             - c2
     >             * (  u(2,i+1,j,k)*u(2,i+1,j,k) * tmp2
     >                  + qs(i+1,j,k) * tmp1 ) )
     >           - dt * tx1
     >           * ( r43*c34 - c1345 ) * tmp2 * u(2,i+1,j,k)
               au(5,3,i,j) = dt * tx2
     >           * ( - c2 * ( u(3,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
     >           - dt * tx1
     >           * (  c34 - c1345 ) * tmp2 * u(3,i+1,j,k)
               au(5,4,i,j) = dt * tx2
     >           * ( - c2 * ( u(4,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
     >           - dt * tx1
     >           * (  c34 - c1345 ) * tmp2 * u(4,i+1,j,k)
               au(5,5,i,j) = dt * tx2
     >           * ( c1 * ( u(2,i+1,j,k) * tmp1 ) )
     >           - dt * tx1 * c1345 * tmp1
     >           - dt * tx1 * dx5

c---------------------------------------------------------------------
c   form the second block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j+1,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               bu(1,1,i,j) = - dt * ty1 * dy1
               bu(1,2,i,j) =   0.0d+00
               bu(1,3,i,j) =  dt * ty2
               bu(1,4,i,j) =   0.0d+00
               bu(1,5,i,j) =   0.0d+00

               bu(2,1,i,j) =  dt * ty2
     >           * ( - ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
     >           - dt * ty1 * ( - c34 * tmp2 * u(2,i,j+1,k) )
               bu(2,2,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
     >          - dt * ty1 * ( c34 * tmp1 )
     >          - dt * ty1 * dy2
               bu(2,3,i,j) =  dt * ty2 * ( u(2,i,j+1,k) * tmp1 )
               bu(2,4,i,j) = 0.0d+00
               bu(2,5,i,j) = 0.0d+00

               bu(3,1,i,j) =  dt * ty2
     >           * ( - ( u(3,i,j+1,k) * tmp1 ) ** 2
     >      + c2 * ( qs(i,j+1,k) * tmp1 ) )
     >       - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j+1,k) )
               bu(3,2,i,j) =  dt * ty2
     >                   * ( - c2 * ( u(2,i,j+1,k) * tmp1 ) )
               bu(3,3,i,j) =  dt * ty2 * ( ( 2.0d+00 - c2 )
     >                   * ( u(3,i,j+1,k) * tmp1 ) )
     >       - dt * ty1 * ( r43 * c34 * tmp1 )
     >       - dt * ty1 * dy3
               bu(3,4,i,j) =  dt * ty2
     >                   * ( - c2 * ( u(4,i,j+1,k) * tmp1 ) )
               bu(3,5,i,j) =  dt * ty2 * c2

               bu(4,1,i,j) =  dt * ty2
     >              * ( - ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
     >       - dt * ty1 * ( - c34 * tmp2 * u(4,i,j+1,k) )
               bu(4,2,i,j) = 0.0d+00
               bu(4,3,i,j) =  dt * ty2 * ( u(4,i,j+1,k) * tmp1 )
               bu(4,4,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
     >                        - dt * ty1 * ( c34 * tmp1 )
     >                        - dt * ty1 * dy4
               bu(4,5,i,j) = 0.0d+00

               bu(5,1,i,j) =  dt * ty2
     >          * ( ( c2 * 2.0d0 * qs(i,j+1,k)
     >               - c1 * u(5,i,j+1,k) )
     >          * ( u(3,i,j+1,k) * tmp2 ) )
     >          - dt * ty1
     >          * ( - (     c34 - c1345 )*tmp3*(u(2,i,j+1,k)**2)
     >              - ( r43*c34 - c1345 )*tmp3*(u(3,i,j+1,k)**2)
     >              - (     c34 - c1345 )*tmp3*(u(4,i,j+1,k)**2)
     >              - c1345*tmp2*u(5,i,j+1,k) )
               bu(5,2,i,j) =  dt * ty2
     >          * ( - c2 * ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
     >          - dt * ty1
     >          * ( c34 - c1345 ) * tmp2 * u(2,i,j+1,k)
               bu(5,3,i,j) =  dt * ty2
     >          * ( c1 * ( u(5,i,j+1,k) * tmp1 )
     >          - c2 
     >          * ( qs(i,j+1,k) * tmp1
     >               + u(3,i,j+1,k)*u(3,i,j+1,k) * tmp2 ) )
     >          - dt * ty1
     >          * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j+1,k)
               bu(5,4,i,j) =  dt * ty2
     >          * ( - c2 * ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
     >          - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j+1,k)
               bu(5,5,i,j) =  dt * ty2
     >          * ( c1 * ( u(3,i,j+1,k) * tmp1 ) )
     >          - dt * ty1 * c1345 * tmp1
     >          - dt * ty1 * dy5

c---------------------------------------------------------------------
c   form the third block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j,k+1)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               cu(1,1,i,j) = - dt * tz1 * dz1
               cu(1,2,i,j) =   0.0d+00
               cu(1,3,i,j) =   0.0d+00
               cu(1,4,i,j) = dt * tz2
               cu(1,5,i,j) =   0.0d+00

               cu(2,1,i,j) = dt * tz2
     >           * ( - ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >           - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k+1) )
               cu(2,2,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
     >           - dt * tz1 * c34 * tmp1
     >           - dt * tz1 * dz2 
               cu(2,3,i,j) = 0.0d+00
               cu(2,4,i,j) = dt * tz2 * ( u(2,i,j,k+1) * tmp1 )
               cu(2,5,i,j) = 0.0d+00

               cu(3,1,i,j) = dt * tz2
     >           * ( - ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >           - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k+1) )
               cu(3,2,i,j) = 0.0d+00
               cu(3,3,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
     >           - dt * tz1 * ( c34 * tmp1 )
     >           - dt * tz1 * dz3
               cu(3,4,i,j) = dt * tz2 * ( u(3,i,j,k+1) * tmp1 )
               cu(3,5,i,j) = 0.0d+00

               cu(4,1,i,j) = dt * tz2
     >        * ( - ( u(4,i,j,k+1) * tmp1 ) ** 2
     >            + c2 * ( qs(i,j,k+1) * tmp1 ) )
     >        - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k+1) )
               cu(4,2,i,j) = dt * tz2
     >             * ( - c2 * ( u(2,i,j,k+1) * tmp1 ) )
               cu(4,3,i,j) = dt * tz2
     >             * ( - c2 * ( u(3,i,j,k+1) * tmp1 ) )
               cu(4,4,i,j) = dt * tz2 * ( 2.0d+00 - c2 )
     >             * ( u(4,i,j,k+1) * tmp1 )
     >             - dt * tz1 * ( r43 * c34 * tmp1 )
     >             - dt * tz1 * dz4
               cu(4,5,i,j) = dt * tz2 * c2

               cu(5,1,i,j) = dt * tz2
     >     * ( ( c2 * 2.0d0 * qs(i,j,k+1)
     >       - c1 * u(5,i,j,k+1) )
     >            * ( u(4,i,j,k+1) * tmp2 ) )
     >       - dt * tz1
     >       * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k+1)**2)
     >           - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k+1)**2)
     >           - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k+1)**2)
     >          - c1345 * tmp2 * u(5,i,j,k+1) )
               cu(5,2,i,j) = dt * tz2
     >       * ( - c2 * ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k+1)
               cu(5,3,i,j) = dt * tz2
     >       * ( - c2 * ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k+1)
               cu(5,4,i,j) = dt * tz2
     >       * ( c1 * ( u(5,i,j,k+1) * tmp1 )
     >       - c2
     >       * ( qs(i,j,k+1) * tmp1
     >            + u(4,i,j,k+1)*u(4,i,j,k+1) * tmp2 ) )
     >       - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k+1)
               cu(5,5,i,j) = dt * tz2
     >       * ( c1 * ( u(4,i,j,k+1) * tmp1 ) )
     >       - dt * tz1 * c1345 * tmp1
     >       - dt * tz1 * dz5

            end do
         end do
!$omp end do nowait

      return
      end