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

      subroutine z_solve

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

c---------------------------------------------------------------------
c     Performs line solves in Z direction by first factoring
c     the block-tridiagonal matrix into an upper triangular matrix, 
c     and then performing back substitution to solve for the unknow
c     vectors of each line.  
c     
c     Make sure we treat elements zero to cell_size in the direction
c     of the sweep.
c---------------------------------------------------------------------

      include 'header.h'
      include 'mpinpb.h'

      integer c, kstart, stage,
     >     first, last, recv_id, error, r_status(MPI_STATUS_SIZE),
     >     isize,jsize,ksize,send_id

      kstart = 0

      if (timeron) call timer_start(t_zsolve)
c---------------------------------------------------------------------
c     in our terminology stage is the number of the cell in the y-direct
c     i.e. stage = 1 means the start of the line stage=ncells means end
c---------------------------------------------------------------------
      do stage = 1,ncells
         c = slice(3,stage)
         isize = cell_size(1,c) - 1
         jsize = cell_size(2,c) - 1
         ksize = cell_size(3,c) - 1
c---------------------------------------------------------------------
c     set last-cell flag
c---------------------------------------------------------------------
         if (stage .eq. ncells) then
            last = 1
         else
            last = 0
         endif

         if (stage .eq. 1) then
c---------------------------------------------------------------------
c     This is the first cell, so solve without receiving data
c---------------------------------------------------------------------
            first = 1
c            call lhsz(c)
            call z_solve_cell(first,last,c)
         else
c---------------------------------------------------------------------
c     Not the first cell of this line, so receive info from
c     processor working on preceeding cell
c---------------------------------------------------------------------
            first = 0
            if (timeron) call timer_start(t_zcomm)
            call z_receive_solve_info(recv_id,c)
c---------------------------------------------------------------------
c     overlap computations and communications
c---------------------------------------------------------------------
c            call lhsz(c)
c---------------------------------------------------------------------
c     wait for completion
c---------------------------------------------------------------------
            call mpi_wait(send_id,r_status,error)
            call mpi_wait(recv_id,r_status,error)
            if (timeron) call timer_stop(t_zcomm)
c---------------------------------------------------------------------
c     install C'(kstart+1) and rhs'(kstart+1) to be used in this cell
c---------------------------------------------------------------------
            call z_unpack_solve_info(c)
            call z_solve_cell(first,last,c)
         endif

         if (last .eq. 0) call z_send_solve_info(send_id,c)
      enddo

c---------------------------------------------------------------------
c     now perform backsubstitution in reverse direction
c---------------------------------------------------------------------
      do stage = ncells, 1, -1
         c = slice(3,stage)
         first = 0
         last = 0
         if (stage .eq. 1) first = 1
         if (stage .eq. ncells) then
            last = 1
c---------------------------------------------------------------------
c     last cell, so perform back substitute without waiting
c---------------------------------------------------------------------
            call z_backsubstitute(first, last,c)
         else
            if (timeron) call timer_start(t_zcomm)
            call z_receive_backsub_info(recv_id,c)
            call mpi_wait(send_id,r_status,error)
            call mpi_wait(recv_id,r_status,error)
            if (timeron) call timer_stop(t_zcomm)
            call z_unpack_backsub_info(c)
            call z_backsubstitute(first,last,c)
         endif
         if (first .eq. 0) call z_send_backsub_info(send_id,c)
      enddo

      if (timeron) call timer_stop(t_zsolve)

      return
      end
      
c---------------------------------------------------------------------
c---------------------------------------------------------------------
      
      subroutine z_unpack_solve_info(c)
c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     unpack C'(-1) and rhs'(-1) for
c     all i and j
c---------------------------------------------------------------------

      include 'header.h'

      integer i,j,m,n,ptr,c,kstart 

      kstart = 0
      ptr = 0
      do j=0,JMAX-1
         do i=0,IMAX-1
            do m=1,BLOCK_SIZE
               do n=1,BLOCK_SIZE
                  lhsc(m,n,i,j,kstart-1,c) = out_buffer(ptr+n)
               enddo
               ptr = ptr+BLOCK_SIZE
            enddo
            do n=1,BLOCK_SIZE
               rhs(n,i,j,kstart-1,c) = out_buffer(ptr+n)
            enddo
            ptr = ptr+BLOCK_SIZE
         enddo
      enddo

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------
      
      subroutine z_send_solve_info(send_id,c)

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

c---------------------------------------------------------------------
c     pack up and send C'(kend) and rhs'(kend) for
c     all i and j
c---------------------------------------------------------------------

      include 'header.h'
      include 'mpinpb.h'

      integer i,j,m,n,ksize,ptr,c,ip,jp
      integer error,send_id,buffer_size

      ksize = cell_size(3,c)-1
      ip = cell_coord(1,c) - 1
      jp = cell_coord(2,c) - 1
      buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*
     >     (BLOCK_SIZE*BLOCK_SIZE + BLOCK_SIZE)

c---------------------------------------------------------------------
c     pack up buffer
c---------------------------------------------------------------------
      ptr = 0
      do j=0,JMAX-1
         do i=0,IMAX-1
            do m=1,BLOCK_SIZE
               do n=1,BLOCK_SIZE
                  in_buffer(ptr+n) = lhsc(m,n,i,j,ksize,c)
               enddo
               ptr = ptr+BLOCK_SIZE
            enddo
            do n=1,BLOCK_SIZE
               in_buffer(ptr+n) = rhs(n,i,j,ksize,c)
            enddo
            ptr = ptr+BLOCK_SIZE
         enddo
      enddo

c---------------------------------------------------------------------
c     send buffer 
c---------------------------------------------------------------------
      if (timeron) call timer_start(t_zcomm)
      call mpi_isend(in_buffer, buffer_size,
     >     dp_type, successor(3),
     >     BOTTOM+ip+jp*NCELLS, comm_solve,
     >     send_id,error)
      if (timeron) call timer_stop(t_zcomm)

      return
      end

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

      subroutine z_send_backsub_info(send_id,c)

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

c---------------------------------------------------------------------
c     pack up and send U(jstart) for all i and j
c---------------------------------------------------------------------

      include 'header.h'
      include 'mpinpb.h'

      integer i,j,n,ptr,c,kstart,ip,jp
      integer error,send_id,buffer_size

c---------------------------------------------------------------------
c     Send element 0 to previous processor
c---------------------------------------------------------------------
      kstart = 0
      ip = cell_coord(1,c)-1
      jp = cell_coord(2,c)-1
      buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*BLOCK_SIZE
      ptr = 0
      do j=0,JMAX-1
         do i=0,IMAX-1
            do n=1,BLOCK_SIZE
               in_buffer(ptr+n) = rhs(n,i,j,kstart,c)
            enddo
            ptr = ptr+BLOCK_SIZE
         enddo
      enddo

      if (timeron) call timer_start(t_zcomm)
      call mpi_isend(in_buffer, buffer_size,
     >     dp_type, predecessor(3), 
     >     TOP+ip+jp*NCELLS, comm_solve, 
     >     send_id,error)
      if (timeron) call timer_stop(t_zcomm)

      return
      end

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

      subroutine z_unpack_backsub_info(c)

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

c---------------------------------------------------------------------
c     unpack U(ksize) for all i and j
c---------------------------------------------------------------------

      include 'header.h'

      integer i,j,n,ptr,c

      ptr = 0
      do j=0,JMAX-1
         do i=0,IMAX-1
            do n=1,BLOCK_SIZE
               backsub_info(n,i,j,c) = out_buffer(ptr+n)
            enddo
            ptr = ptr+BLOCK_SIZE
         enddo
      enddo

      return
      end


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

      subroutine z_receive_backsub_info(recv_id,c)

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

c---------------------------------------------------------------------
c     post mpi receives
c---------------------------------------------------------------------

      include 'header.h'
      include 'mpinpb.h'

      integer error,recv_id,ip,jp,c,buffer_size
      ip = cell_coord(1,c) - 1
      jp = cell_coord(2,c) - 1
      buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*BLOCK_SIZE
      call mpi_irecv(out_buffer, buffer_size,
     >     dp_type, successor(3), 
     >     TOP+ip+jp*NCELLS, comm_solve, 
     >     recv_id, error)

      return
      end

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

      subroutine z_receive_solve_info(recv_id,c)

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

c---------------------------------------------------------------------
c     post mpi receives 
c---------------------------------------------------------------------

      include 'header.h'
      include 'mpinpb.h'

      integer ip,jp,recv_id,error,c,buffer_size
      ip = cell_coord(1,c) - 1
      jp = cell_coord(2,c) - 1
      buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*
     >     (BLOCK_SIZE*BLOCK_SIZE + BLOCK_SIZE)
      call mpi_irecv(out_buffer, buffer_size,
     >     dp_type, predecessor(3), 
     >     BOTTOM+ip+jp*NCELLS, comm_solve,
     >     recv_id, error)

      return
      end
      
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine z_backsubstitute(first, last, c)

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

c---------------------------------------------------------------------
c     back solve: if last cell, then generate U(ksize)=rhs(ksize)
c     else assume U(ksize) is loaded in un pack backsub_info
c     so just use it
c     after call u(kstart) will be sent to next cell
c---------------------------------------------------------------------

      include 'header.h'

      integer first, last, c, i, k
      integer m,n,j,jsize,isize,ksize,kstart
      
      kstart = 0
      isize = cell_size(1,c)-end(1,c)-1      
      jsize = cell_size(2,c)-end(2,c)-1
      ksize = cell_size(3,c)-1
      if (last .eq. 0) then
         do j=start(2,c),jsize
            do i=start(1,c),isize
c---------------------------------------------------------------------
c     U(jsize) uses info from previous cell if not last cell
c---------------------------------------------------------------------
               do m=1,BLOCK_SIZE
                  do n=1,BLOCK_SIZE
                     rhs(m,i,j,ksize,c) = rhs(m,i,j,ksize,c) 
     >                    - lhsc(m,n,i,j,ksize,c)*
     >                    backsub_info(n,i,j,c)
                  enddo
               enddo
            enddo
         enddo
      endif
      do k=ksize-1,kstart,-1
         do j=start(2,c),jsize
            do i=start(1,c),isize
               do m=1,BLOCK_SIZE
                  do n=1,BLOCK_SIZE
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) 
     >                    - lhsc(m,n,i,j,k,c)*rhs(n,i,j,k+1,c)
                  enddo
               enddo
            enddo
         enddo
      enddo

      return
      end

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

      subroutine z_solve_cell(first,last,c)

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

c---------------------------------------------------------------------
c     performs guaussian elimination on this cell.
c     
c     assumes that unpacking routines for non-first cells 
c     preload C' and rhs' from previous cell.
c     
c     assumed send happens outside this routine, but that
c     c'(KMAX) and rhs'(KMAX) will be sent to next cell.
c---------------------------------------------------------------------

      include 'header.h'
      include 'work_lhs_vec.h'

      integer first,last,c
      integer i,j,k,m,n,isize,ksize,jsize,kstart

      kstart = 0
      isize = cell_size(1,c)-end(1,c)-1
      jsize = cell_size(2,c)-end(2,c)-1
      ksize = cell_size(3,c)-1

c---------------------------------------------------------------------
c     zero the left hand side for starters
c     set diagonal values to 1. This is overkill, but convenient
c---------------------------------------------------------------------
      do i = 0, isize
         do m = 1, 5
            do n = 1, 5
               lhsa(m,n,i,0) = 0.0d0
               lhsb(m,n,i,0) = 0.0d0
               lhsa(m,n,i,ksize) = 0.0d0
               lhsb(m,n,i,ksize) = 0.0d0
            enddo
            lhsb(m,m,i,0) = 1.0d0
            lhsb(m,m,i,ksize) = 1.0d0
         enddo
      enddo

      do j=start(2,c),jsize 

c---------------------------------------------------------------------
c     This function computes the left hand side for the three z-factors 
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     Compute the indices for storing the block-diagonal matrix;
c     determine c (labeled f) and s jacobians for cell c
c---------------------------------------------------------------------

         do k = start(3,c)-1, cell_size(3,c)-end(3,c)
            do i=start(1,c),isize

               tmp1 = 1.0d0 / u(1,i,j,k,c)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               fjac(1,1,i,k) = 0.0d+00
               fjac(1,2,i,k) = 0.0d+00
               fjac(1,3,i,k) = 0.0d+00
               fjac(1,4,i,k) = 1.0d+00
               fjac(1,5,i,k) = 0.0d+00

               fjac(2,1,i,k) = - ( u(2,i,j,k,c)*u(4,i,j,k,c) ) 
     >              * tmp2 
               fjac(2,2,i,k) = u(4,i,j,k,c) * tmp1
               fjac(2,3,i,k) = 0.0d+00
               fjac(2,4,i,k) = u(2,i,j,k,c) * tmp1
               fjac(2,5,i,k) = 0.0d+00

               fjac(3,1,i,k) = - ( u(3,i,j,k,c)*u(4,i,j,k,c) )
     >              * tmp2 
               fjac(3,2,i,k) = 0.0d+00
               fjac(3,3,i,k) = u(4,i,j,k,c) * tmp1
               fjac(3,4,i,k) = u(3,i,j,k,c) * tmp1
               fjac(3,5,i,k) = 0.0d+00

               fjac(4,1,i,k) = - (u(4,i,j,k,c)*u(4,i,j,k,c) * tmp2 ) 
     >              + c2 * qs(i,j,k,c)
               fjac(4,2,i,k) = - c2 *  u(2,i,j,k,c) * tmp1 
               fjac(4,3,i,k) = - c2 *  u(3,i,j,k,c) * tmp1
               fjac(4,4,i,k) = ( 2.0d+00 - c2 )
     >              *  u(4,i,j,k,c) * tmp1 
               fjac(4,5,i,k) = c2

               fjac(5,1,i,k) = ( c2 * 2.0d0 * qs(i,j,k,c)
     >              - c1 * ( u(5,i,j,k,c) * tmp1 ) )
     >              * ( u(4,i,j,k,c) * tmp1 )
               fjac(5,2,i,k) = - c2 * ( u(2,i,j,k,c)*u(4,i,j,k,c) )
     >              * tmp2 
               fjac(5,3,i,k) = - c2 * ( u(3,i,j,k,c)*u(4,i,j,k,c) )
     >              * tmp2
               fjac(5,4,i,k) = c1 * ( u(5,i,j,k,c) * tmp1 )
     >              - c2 * ( qs(i,j,k,c)
     >              + u(4,i,j,k,c)*u(4,i,j,k,c) * tmp2 )
               fjac(5,5,i,k) = c1 * u(4,i,j,k,c) * tmp1

               njac(1,1,i,k) = 0.0d+00
               njac(1,2,i,k) = 0.0d+00
               njac(1,3,i,k) = 0.0d+00
               njac(1,4,i,k) = 0.0d+00
               njac(1,5,i,k) = 0.0d+00

               njac(2,1,i,k) = - c3c4 * tmp2 * u(2,i,j,k,c)
               njac(2,2,i,k) =   c3c4 * tmp1
               njac(2,3,i,k) =   0.0d+00
               njac(2,4,i,k) =   0.0d+00
               njac(2,5,i,k) =   0.0d+00

               njac(3,1,i,k) = - c3c4 * tmp2 * u(3,i,j,k,c)
               njac(3,2,i,k) =   0.0d+00
               njac(3,3,i,k) =   c3c4 * tmp1
               njac(3,4,i,k) =   0.0d+00
               njac(3,5,i,k) =   0.0d+00

               njac(4,1,i,k) = - con43 * c3c4 * tmp2 * u(4,i,j,k,c)
               njac(4,2,i,k) =   0.0d+00
               njac(4,3,i,k) =   0.0d+00
               njac(4,4,i,k) =   con43 * c3 * c4 * tmp1
               njac(4,5,i,k) =   0.0d+00

               njac(5,1,i,k) = - (  c3c4
     >              - c1345 ) * tmp3 * (u(2,i,j,k,c)**2)
     >              - ( c3c4 - c1345 ) * tmp3 * (u(3,i,j,k,c)**2)
     >              - ( con43 * c3c4
     >              - c1345 ) * tmp3 * (u(4,i,j,k,c)**2)
     >              - c1345 * tmp2 * u(5,i,j,k,c)

               njac(5,2,i,k) = (  c3c4 - c1345 ) * tmp2 * u(2,i,j,k,c)
               njac(5,3,i,k) = (  c3c4 - c1345 ) * tmp2 * u(3,i,j,k,c)
               njac(5,4,i,k) = ( con43 * c3c4
     >              - c1345 ) * tmp2 * u(4,i,j,k,c)
               njac(5,5,i,k) = ( c1345 )* tmp1


            enddo
         enddo

c---------------------------------------------------------------------
c     now joacobians set, so form left hand side in z direction
c---------------------------------------------------------------------
         do k = start(3,c), ksize-end(3,c)
            do i=start(1,c),isize

               tmp1 = dt * tz1
               tmp2 = dt * tz2

               lhsa(1,1,i,k) = - tmp2 * fjac(1,1,i,k-1)
     >              - tmp1 * njac(1,1,i,k-1)
     >              - tmp1 * dz1 
               lhsa(1,2,i,k) = - tmp2 * fjac(1,2,i,k-1)
     >              - tmp1 * njac(1,2,i,k-1)
               lhsa(1,3,i,k) = - tmp2 * fjac(1,3,i,k-1)
     >              - tmp1 * njac(1,3,i,k-1)
               lhsa(1,4,i,k) = - tmp2 * fjac(1,4,i,k-1)
     >              - tmp1 * njac(1,4,i,k-1)
               lhsa(1,5,i,k) = - tmp2 * fjac(1,5,i,k-1)
     >              - tmp1 * njac(1,5,i,k-1)

               lhsa(2,1,i,k) = - tmp2 * fjac(2,1,i,k-1)
     >              - tmp1 * njac(2,1,i,k-1)
               lhsa(2,2,i,k) = - tmp2 * fjac(2,2,i,k-1)
     >              - tmp1 * njac(2,2,i,k-1)
     >              - tmp1 * dz2
               lhsa(2,3,i,k) = - tmp2 * fjac(2,3,i,k-1)
     >              - tmp1 * njac(2,3,i,k-1)
               lhsa(2,4,i,k) = - tmp2 * fjac(2,4,i,k-1)
     >              - tmp1 * njac(2,4,i,k-1)
               lhsa(2,5,i,k) = - tmp2 * fjac(2,5,i,k-1)
     >              - tmp1 * njac(2,5,i,k-1)

               lhsa(3,1,i,k) = - tmp2 * fjac(3,1,i,k-1)
     >              - tmp1 * njac(3,1,i,k-1)
               lhsa(3,2,i,k) = - tmp2 * fjac(3,2,i,k-1)
     >              - tmp1 * njac(3,2,i,k-1)
               lhsa(3,3,i,k) = - tmp2 * fjac(3,3,i,k-1)
     >              - tmp1 * njac(3,3,i,k-1)
     >              - tmp1 * dz3 
               lhsa(3,4,i,k) = - tmp2 * fjac(3,4,i,k-1)
     >              - tmp1 * njac(3,4,i,k-1)
               lhsa(3,5,i,k) = - tmp2 * fjac(3,5,i,k-1)
     >              - tmp1 * njac(3,5,i,k-1)

               lhsa(4,1,i,k) = - tmp2 * fjac(4,1,i,k-1)
     >              - tmp1 * njac(4,1,i,k-1)
               lhsa(4,2,i,k) = - tmp2 * fjac(4,2,i,k-1)
     >              - tmp1 * njac(4,2,i,k-1)
               lhsa(4,3,i,k) = - tmp2 * fjac(4,3,i,k-1)
     >              - tmp1 * njac(4,3,i,k-1)
               lhsa(4,4,i,k) = - tmp2 * fjac(4,4,i,k-1)
     >              - tmp1 * njac(4,4,i,k-1)
     >              - tmp1 * dz4
               lhsa(4,5,i,k) = - tmp2 * fjac(4,5,i,k-1)
     >              - tmp1 * njac(4,5,i,k-1)

               lhsa(5,1,i,k) = - tmp2 * fjac(5,1,i,k-1)
     >              - tmp1 * njac(5,1,i,k-1)
               lhsa(5,2,i,k) = - tmp2 * fjac(5,2,i,k-1)
     >              - tmp1 * njac(5,2,i,k-1)
               lhsa(5,3,i,k) = - tmp2 * fjac(5,3,i,k-1)
     >              - tmp1 * njac(5,3,i,k-1)
               lhsa(5,4,i,k) = - tmp2 * fjac(5,4,i,k-1)
     >              - tmp1 * njac(5,4,i,k-1)
               lhsa(5,5,i,k) = - tmp2 * fjac(5,5,i,k-1)
     >              - tmp1 * njac(5,5,i,k-1)
     >              - tmp1 * dz5

               lhsb(1,1,i,k) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(1,1,i,k)
     >              + tmp1 * 2.0d+00 * dz1
               lhsb(1,2,i,k) = tmp1 * 2.0d+00 * njac(1,2,i,k)
               lhsb(1,3,i,k) = tmp1 * 2.0d+00 * njac(1,3,i,k)
               lhsb(1,4,i,k) = tmp1 * 2.0d+00 * njac(1,4,i,k)
               lhsb(1,5,i,k) = tmp1 * 2.0d+00 * njac(1,5,i,k)

               lhsb(2,1,i,k) = tmp1 * 2.0d+00 * njac(2,1,i,k)
               lhsb(2,2,i,k) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(2,2,i,k)
     >              + tmp1 * 2.0d+00 * dz2
               lhsb(2,3,i,k) = tmp1 * 2.0d+00 * njac(2,3,i,k)
               lhsb(2,4,i,k) = tmp1 * 2.0d+00 * njac(2,4,i,k)
               lhsb(2,5,i,k) = tmp1 * 2.0d+00 * njac(2,5,i,k)

               lhsb(3,1,i,k) = tmp1 * 2.0d+00 * njac(3,1,i,k)
               lhsb(3,2,i,k) = tmp1 * 2.0d+00 * njac(3,2,i,k)
               lhsb(3,3,i,k) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(3,3,i,k)
     >              + tmp1 * 2.0d+00 * dz3
               lhsb(3,4,i,k) = tmp1 * 2.0d+00 * njac(3,4,i,k)
               lhsb(3,5,i,k) = tmp1 * 2.0d+00 * njac(3,5,i,k)

               lhsb(4,1,i,k) = tmp1 * 2.0d+00 * njac(4,1,i,k)
               lhsb(4,2,i,k) = tmp1 * 2.0d+00 * njac(4,2,i,k)
               lhsb(4,3,i,k) = tmp1 * 2.0d+00 * njac(4,3,i,k)
               lhsb(4,4,i,k) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(4,4,i,k)
     >              + tmp1 * 2.0d+00 * dz4
               lhsb(4,5,i,k) = tmp1 * 2.0d+00 * njac(4,5,i,k)

               lhsb(5,1,i,k) = tmp1 * 2.0d+00 * njac(5,1,i,k)
               lhsb(5,2,i,k) = tmp1 * 2.0d+00 * njac(5,2,i,k)
               lhsb(5,3,i,k) = tmp1 * 2.0d+00 * njac(5,3,i,k)
               lhsb(5,4,i,k) = tmp1 * 2.0d+00 * njac(5,4,i,k)
               lhsb(5,5,i,k) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(5,5,i,k) 
     >              + tmp1 * 2.0d+00 * dz5

               lhsc(1,1,i,j,k,c) =  tmp2 * fjac(1,1,i,k+1)
     >              - tmp1 * njac(1,1,i,k+1)
     >              - tmp1 * dz1
               lhsc(1,2,i,j,k,c) =  tmp2 * fjac(1,2,i,k+1)
     >              - tmp1 * njac(1,2,i,k+1)
               lhsc(1,3,i,j,k,c) =  tmp2 * fjac(1,3,i,k+1)
     >              - tmp1 * njac(1,3,i,k+1)
               lhsc(1,4,i,j,k,c) =  tmp2 * fjac(1,4,i,k+1)
     >              - tmp1 * njac(1,4,i,k+1)
               lhsc(1,5,i,j,k,c) =  tmp2 * fjac(1,5,i,k+1)
     >              - tmp1 * njac(1,5,i,k+1)

               lhsc(2,1,i,j,k,c) =  tmp2 * fjac(2,1,i,k+1)
     >              - tmp1 * njac(2,1,i,k+1)
               lhsc(2,2,i,j,k,c) =  tmp2 * fjac(2,2,i,k+1)
     >              - tmp1 * njac(2,2,i,k+1)
     >              - tmp1 * dz2
               lhsc(2,3,i,j,k,c) =  tmp2 * fjac(2,3,i,k+1)
     >              - tmp1 * njac(2,3,i,k+1)
               lhsc(2,4,i,j,k,c) =  tmp2 * fjac(2,4,i,k+1)
     >              - tmp1 * njac(2,4,i,k+1)
               lhsc(2,5,i,j,k,c) =  tmp2 * fjac(2,5,i,k+1)
     >              - tmp1 * njac(2,5,i,k+1)

               lhsc(3,1,i,j,k,c) =  tmp2 * fjac(3,1,i,k+1)
     >              - tmp1 * njac(3,1,i,k+1)
               lhsc(3,2,i,j,k,c) =  tmp2 * fjac(3,2,i,k+1)
     >              - tmp1 * njac(3,2,i,k+1)
               lhsc(3,3,i,j,k,c) =  tmp2 * fjac(3,3,i,k+1)
     >              - tmp1 * njac(3,3,i,k+1)
     >              - tmp1 * dz3
               lhsc(3,4,i,j,k,c) =  tmp2 * fjac(3,4,i,k+1)
     >              - tmp1 * njac(3,4,i,k+1)
               lhsc(3,5,i,j,k,c) =  tmp2 * fjac(3,5,i,k+1)
     >              - tmp1 * njac(3,5,i,k+1)

               lhsc(4,1,i,j,k,c) =  tmp2 * fjac(4,1,i,k+1)
     >              - tmp1 * njac(4,1,i,k+1)
               lhsc(4,2,i,j,k,c) =  tmp2 * fjac(4,2,i,k+1)
     >              - tmp1 * njac(4,2,i,k+1)
               lhsc(4,3,i,j,k,c) =  tmp2 * fjac(4,3,i,k+1)
     >              - tmp1 * njac(4,3,i,k+1)
               lhsc(4,4,i,j,k,c) =  tmp2 * fjac(4,4,i,k+1)
     >              - tmp1 * njac(4,4,i,k+1)
     >              - tmp1 * dz4
               lhsc(4,5,i,j,k,c) =  tmp2 * fjac(4,5,i,k+1)
     >              - tmp1 * njac(4,5,i,k+1)

               lhsc(5,1,i,j,k,c) =  tmp2 * fjac(5,1,i,k+1)
     >              - tmp1 * njac(5,1,i,k+1)
               lhsc(5,2,i,j,k,c) =  tmp2 * fjac(5,2,i,k+1)
     >              - tmp1 * njac(5,2,i,k+1)
               lhsc(5,3,i,j,k,c) =  tmp2 * fjac(5,3,i,k+1)
     >              - tmp1 * njac(5,3,i,k+1)
               lhsc(5,4,i,j,k,c) =  tmp2 * fjac(5,4,i,k+1)
     >              - tmp1 * njac(5,4,i,k+1)
               lhsc(5,5,i,j,k,c) =  tmp2 * fjac(5,5,i,k+1)
     >              - tmp1 * njac(5,5,i,k+1)
     >              - tmp1 * dz5

            enddo
         enddo


c---------------------------------------------------------------------
c     outer most do loops - sweeping in i direction
c---------------------------------------------------------------------
         if (first .eq. 1) then 

c---------------------------------------------------------------------
c     multiply c(i,j,kstart) by b_inverse and copy back to c
c     multiply rhs(kstart) by b_inverse(kstart) and copy to rhs
c---------------------------------------------------------------------
!dir$ ivdep
            do i=start(1,c),isize
               call binvcrhs( lhsb(1,1,i,kstart),
     >                        lhsc(1,1,i,j,kstart,c),
     >                        rhs(1,i,j,kstart,c) )
            enddo

         endif

c---------------------------------------------------------------------
c     begin inner most do loop
c     do all the elements of the cell unless last 
c---------------------------------------------------------------------
         do k=kstart+first,ksize-last
!dir$ ivdep
            do i=start(1,c),isize

c---------------------------------------------------------------------
c     subtract A*lhs_vector(k-1) from lhs_vector(k)
c     
c     rhs(k) = rhs(k) - A*rhs(k-1)
c---------------------------------------------------------------------
               call matvec_sub(lhsa(1,1,i,k),
     >                         rhs(1,i,j,k-1,c),rhs(1,i,j,k,c))

c---------------------------------------------------------------------
c     B(k) = B(k) - C(k-1)*A(k)
c     call matmul_sub(aa,i,j,k,c,cc,i,j,k-1,c,bb,i,j,k,c)
c---------------------------------------------------------------------
               call matmul_sub(lhsa(1,1,i,k),
     >                         lhsc(1,1,i,j,k-1,c),
     >                         lhsb(1,1,i,k))

c---------------------------------------------------------------------
c     multiply c(i,j,k) by b_inverse and copy back to c
c     multiply rhs(i,j,1) by b_inverse(i,j,1) and copy to rhs
c---------------------------------------------------------------------
               call binvcrhs( lhsb(1,1,i,k),
     >                        lhsc(1,1,i,j,k,c),
     >                        rhs(1,i,j,k,c) )

            enddo
         enddo

c---------------------------------------------------------------------
c     Now finish up special cases for last cell
c---------------------------------------------------------------------
         if (last .eq. 1) then

!dir$ ivdep
            do i=start(1,c),isize
c---------------------------------------------------------------------
c     rhs(ksize) = rhs(ksize) - A*rhs(ksize-1)
c---------------------------------------------------------------------
               call matvec_sub(lhsa(1,1,i,ksize),
     >                         rhs(1,i,j,ksize-1,c),rhs(1,i,j,ksize,c))

c---------------------------------------------------------------------
c     B(ksize) = B(ksize) - C(ksize-1)*A(ksize)
c     call matmul_sub(aa,i,j,ksize,c,
c     $              cc,i,j,ksize-1,c,bb,i,j,ksize,c)
c---------------------------------------------------------------------
               call matmul_sub(lhsa(1,1,i,ksize),
     >                         lhsc(1,1,i,j,ksize-1,c),
     >                         lhsb(1,1,i,ksize))

c---------------------------------------------------------------------
c     multiply rhs(ksize) by b_inverse(ksize) and copy to rhs
c---------------------------------------------------------------------
               call binvrhs( lhsb(1,1,i,ksize),
     >                       rhs(1,i,j,ksize,c) )
            enddo

         endif
      enddo


      return
      end