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

       subroutine y_solve

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

c---------------------------------------------------------------------
c this function performs the solution of the approximate factorization
c step in the y-direction for all five matrix components
c simultaneously. The Thomas algorithm is employed to solve the
c systems for the y-lines. Boundary conditions are non-periodic
c---------------------------------------------------------------------

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

       integer i, j, k, stage, ip, kp, n, isize, jend, ksize, j1, j2,
     >         buffer_size, c, m, p, jstart, error,
     >         requests(2), statuses(MPI_STATUS_SIZE, 2)
       double precision  r1, r2, d, e, s(5), sm1, sm2,
     >                   fac1, fac2


c---------------------------------------------------------------------
c now do a sweep on a layer-by-layer basis, i.e. sweeping through cells
c on this node in the direction of increasing i for the forward sweep,
c and after that reversing the direction for the backsubstitution  
c---------------------------------------------------------------------

       if (timeron) call timer_start(t_ysolve)
c---------------------------------------------------------------------
c                          FORWARD ELIMINATION  
c---------------------------------------------------------------------
       do    stage = 1, ncells
          c      = slice(2,stage)

          jstart = 0
          jend   = cell_size(2,c)-1

          isize     = cell_size(1,c)
          ksize     = cell_size(3,c)
          ip        = cell_coord(1,c)-1
          kp        = cell_coord(3,c)-1

          buffer_size = (isize-start(1,c)-end(1,c)) * 
     >                  (ksize-start(3,c)-end(3,c))

          if ( stage .ne. 1) then

c---------------------------------------------------------------------
c            if this is not the first processor in this row of cells, 
c            receive data from predecessor containing the right hand
c            sides and the upper diagonal elements of the previous two rows
c---------------------------------------------------------------------

             if (timeron) call timer_start(t_ycomm)
             call mpi_irecv(in_buffer, 22*buffer_size, 
     >                      dp_type, predecessor(2), 
     >                      DEFAULT_TAG, comm_solve, 
     >                      requests(1), error)
             if (timeron) call timer_stop(t_ycomm)

c---------------------------------------------------------------------
c            communication has already been started. 
c            compute the left hand side while waiting for the msg
c---------------------------------------------------------------------
             call lhsy(c)

c---------------------------------------------------------------------
c            wait for pending communication to complete
c            This waits on the current receive and on the send
c            from the previous stage. They always come in pairs. 
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_ycomm)
             call mpi_waitall(2, requests, statuses, error)
             if (timeron) call timer_stop(t_ycomm)

c---------------------------------------------------------------------
c            unpack the buffer                                 
c---------------------------------------------------------------------
             j  = jstart
             j1 = jstart + 1
             n = 0
c---------------------------------------------------------------------
c            create a running pointer
c---------------------------------------------------------------------
             p = 0
             do    k = start(3,c), ksize-end(3,c)-1
                do    i = start(1,c), isize-end(1,c)-1
                   lhs(i,j,k,n+2,c) = lhs(i,j,k,n+2,c) -
     >                       in_buffer(p+1) * lhs(i,j,k,n+1,c)
                   lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) -
     >                       in_buffer(p+2) * lhs(i,j,k,n+1,c)
                   do    m = 1, 3
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) -
     >                       in_buffer(p+2+m) * lhs(i,j,k,n+1,c)
                   end do
                   d            = in_buffer(p+6)
                   e            = in_buffer(p+7)
                   do    m = 1, 3
                      s(m) = in_buffer(p+7+m)
                   end do
                   r1 = lhs(i,j,k,n+2,c)
                   lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) - d * r1
                   lhs(i,j,k,n+4,c) = lhs(i,j,k,n+4,c) - e * r1
                   do    m = 1, 3
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - s(m) * r1
                   end do
                   r2 = lhs(i,j1,k,n+1,c)
                   lhs(i,j1,k,n+2,c) = lhs(i,j1,k,n+2,c) - d * r2
                   lhs(i,j1,k,n+3,c) = lhs(i,j1,k,n+3,c) - e * r2
                   do    m = 1, 3
                      rhs(i,j1,k,m,c) = rhs(i,j1,k,m,c) - s(m) * r2
                   end do
                   p = p + 10
                end do
             end do

             do    m = 4, 5
                n = (m-3)*5
                do    k = start(3,c), ksize-end(3,c)-1
                   do    i = start(1,c), isize-end(1,c)-1
                      lhs(i,j,k,n+2,c) = lhs(i,j,k,n+2,c) -
     >                          in_buffer(p+1) * lhs(i,j,k,n+1,c)
                      lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) -
     >                          in_buffer(p+2) * lhs(i,j,k,n+1,c)
                      rhs(i,j,k,m,c)   = rhs(i,j,k,m,c) -
     >                          in_buffer(p+3) * lhs(i,j,k,n+1,c)
                      d                = in_buffer(p+4)
                      e                = in_buffer(p+5)
                      s(m)             = in_buffer(p+6)
                      r1 = lhs(i,j,k,n+2,c)
                      lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) - d * r1
                      lhs(i,j,k,n+4,c) = lhs(i,j,k,n+4,c) - e * r1
                      rhs(i,j,k,m,c)   = rhs(i,j,k,m,c) - s(m) * r1
                      r2 = lhs(i,j1,k,n+1,c)
                      lhs(i,j1,k,n+2,c) = lhs(i,j1,k,n+2,c) - d * r2
                      lhs(i,j1,k,n+3,c) = lhs(i,j1,k,n+3,c) - e * r2
                      rhs(i,j1,k,m,c)   = rhs(i,j1,k,m,c) - s(m) * r2
                      p = p + 6
                   end do
                end do
             end do

          else            

c---------------------------------------------------------------------
c            if this IS the first cell, we still compute the lhs
c---------------------------------------------------------------------
             call lhsy(c)
          endif

c---------------------------------------------------------------------
c         perform the Thomas algorithm; first, FORWARD ELIMINATION     
c---------------------------------------------------------------------
          n = 0

          do    k = start(3,c), ksize-end(3,c)-1
             do    j = jstart, jend-2
                do    i = start(1,c), isize-end(1,c)-1
                   j1 = j  + 1
                   j2 = j  + 2
                   fac1               = 1.d0/lhs(i,j,k,n+3,c)
                   lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                   lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                   do    m = 1, 3
                      rhs(i,j,k,m,c) = fac1*rhs(i,j,k,m,c)
                   end do
                   lhs(i,j1,k,n+3,c) = lhs(i,j1,k,n+3,c) -
     >                         lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+4,c)
                   lhs(i,j1,k,n+4,c) = lhs(i,j1,k,n+4,c) -
     >                         lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+5,c)
                   do    m = 1, 3
                      rhs(i,j1,k,m,c) = rhs(i,j1,k,m,c) -
     >                         lhs(i,j1,k,n+2,c)*rhs(i,j,k,m,c)
                   end do
                   lhs(i,j2,k,n+2,c) = lhs(i,j2,k,n+2,c) -
     >                         lhs(i,j2,k,n+1,c)*lhs(i,j,k,n+4,c)
                   lhs(i,j2,k,n+3,c) = lhs(i,j2,k,n+3,c) -
     >                         lhs(i,j2,k,n+1,c)*lhs(i,j,k,n+5,c)
                   do    m = 1, 3
                      rhs(i,j2,k,m,c) = rhs(i,j2,k,m,c) -
     >                         lhs(i,j2,k,n+1,c)*rhs(i,j,k,m,c)
                   end do
                end do
             end do
          end do

c---------------------------------------------------------------------
c         The last two rows in this grid block are a bit different, 
c         since they do not have two more rows available for the
c         elimination of off-diagonal entries
c---------------------------------------------------------------------

          j  = jend - 1
          j1 = jend
          do    k = start(3,c), ksize-end(3,c)-1
             do    i = start(1,c), isize-end(1,c)-1
                fac1               = 1.d0/lhs(i,j,k,n+3,c)
                lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                do    m = 1, 3
                   rhs(i,j,k,m,c) = fac1*rhs(i,j,k,m,c)
                end do
                lhs(i,j1,k,n+3,c) = lhs(i,j1,k,n+3,c) -
     >                      lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+4,c)
                lhs(i,j1,k,n+4,c) = lhs(i,j1,k,n+4,c) -
     >                      lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+5,c)
                do    m = 1, 3
                   rhs(i,j1,k,m,c) = rhs(i,j1,k,m,c) -
     >                      lhs(i,j1,k,n+2,c)*rhs(i,j,k,m,c)
                end do
c---------------------------------------------------------------------
c               scale the last row immediately (some of this is
c               overkill in case this is the last cell)
c---------------------------------------------------------------------
                fac2               = 1.d0/lhs(i,j1,k,n+3,c)
                lhs(i,j1,k,n+4,c) = fac2*lhs(i,j1,k,n+4,c)
                lhs(i,j1,k,n+5,c) = fac2*lhs(i,j1,k,n+5,c)  
                do    m = 1, 3
                   rhs(i,j1,k,m,c) = fac2*rhs(i,j1,k,m,c)
                end do
             end do
          end do

c---------------------------------------------------------------------
c         do the u+c and the u-c factors                 
c---------------------------------------------------------------------
          do    m = 4, 5
             n = (m-3)*5
             do    k = start(3,c), ksize-end(3,c)-1
                do    j = jstart, jend-2
                   do    i = start(1,c), isize-end(1,c)-1
                   j1 = j  + 1
                   j2 = j  + 2
                   fac1               = 1.d0/lhs(i,j,k,n+3,c)
                   lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                   lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                   rhs(i,j,k,m,c) = fac1*rhs(i,j,k,m,c)
                   lhs(i,j1,k,n+3,c) = lhs(i,j1,k,n+3,c) -
     >                         lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+4,c)
                   lhs(i,j1,k,n+4,c) = lhs(i,j1,k,n+4,c) -
     >                         lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+5,c)
                   rhs(i,j1,k,m,c) = rhs(i,j1,k,m,c) -
     >                         lhs(i,j1,k,n+2,c)*rhs(i,j,k,m,c)
                   lhs(i,j2,k,n+2,c) = lhs(i,j2,k,n+2,c) -
     >                         lhs(i,j2,k,n+1,c)*lhs(i,j,k,n+4,c)
                   lhs(i,j2,k,n+3,c) = lhs(i,j2,k,n+3,c) -
     >                         lhs(i,j2,k,n+1,c)*lhs(i,j,k,n+5,c)
                   rhs(i,j2,k,m,c) = rhs(i,j2,k,m,c) -
     >                         lhs(i,j2,k,n+1,c)*rhs(i,j,k,m,c)
                end do
             end do
          end do

c---------------------------------------------------------------------
c            And again the last two rows separately
c---------------------------------------------------------------------
             j  = jend - 1
             j1 = jend
             do    k = start(3,c), ksize-end(3,c)-1
                do    i = start(1,c), isize-end(1,c)-1
                fac1               = 1.d0/lhs(i,j,k,n+3,c)
                lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                rhs(i,j,k,m,c)     = fac1*rhs(i,j,k,m,c)
                lhs(i,j1,k,n+3,c) = lhs(i,j1,k,n+3,c) -
     >                      lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+4,c)
                lhs(i,j1,k,n+4,c) = lhs(i,j1,k,n+4,c) -
     >                      lhs(i,j1,k,n+2,c)*lhs(i,j,k,n+5,c)
                rhs(i,j1,k,m,c)   = rhs(i,j1,k,m,c) -
     >                      lhs(i,j1,k,n+2,c)*rhs(i,j,k,m,c)
c---------------------------------------------------------------------
c               Scale the last row immediately (some of this is overkill
c               if this is the last cell)
c---------------------------------------------------------------------
                fac2               = 1.d0/lhs(i,j1,k,n+3,c)
                lhs(i,j1,k,n+4,c) = fac2*lhs(i,j1,k,n+4,c)
                lhs(i,j1,k,n+5,c) = fac2*lhs(i,j1,k,n+5,c)
                rhs(i,j1,k,m,c)   = fac2*rhs(i,j1,k,m,c)

             end do
          end do
       end do

c---------------------------------------------------------------------
c         send information to the next processor, except when this
c         is the last grid block;
c---------------------------------------------------------------------

          if (stage .ne. ncells) then

c---------------------------------------------------------------------
c            create a running pointer for the send buffer  
c---------------------------------------------------------------------
             p = 0
             n = 0
             do    k = start(3,c), ksize-end(3,c)-1
                do    i = start(1,c), isize-end(1,c)-1
                   do    j = jend-1, jend
                      out_buffer(p+1) = lhs(i,j,k,n+4,c)
                      out_buffer(p+2) = lhs(i,j,k,n+5,c)
                      do    m = 1, 3
                         out_buffer(p+2+m) = rhs(i,j,k,m,c)
                      end do
                      p = p+5
                   end do
                end do
             end do

             do    m = 4, 5
                n = (m-3)*5
                do    k = start(3,c), ksize-end(3,c)-1
                   do    i = start(1,c), isize-end(1,c)-1
                      do    j = jend-1, jend
                         out_buffer(p+1) = lhs(i,j,k,n+4,c)
                         out_buffer(p+2) = lhs(i,j,k,n+5,c)
                         out_buffer(p+3) = rhs(i,j,k,m,c)
                         p = p + 3
                      end do
                   end do
                end do
             end do

c---------------------------------------------------------------------
c            pack and send the buffer
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_ycomm)
             call mpi_isend(out_buffer, 22*buffer_size, 
     >                     dp_type, successor(2), 
     >                     DEFAULT_TAG, comm_solve, 
     >                     requests(2), error)
             if (timeron) call timer_stop(t_ycomm)

          endif
       end do

c---------------------------------------------------------------------
c      now go in the reverse direction                      
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c                         BACKSUBSTITUTION 
c---------------------------------------------------------------------
       do    stage = ncells, 1, -1
          c = slice(2,stage)

          jstart = 0
          jend   = cell_size(2,c)-1

          isize = cell_size(1,c)
          ksize = cell_size(3,c)
          ip    = cell_coord(1,c)-1
          kp    = cell_coord(3,c)-1

          buffer_size = (isize-start(1,c)-end(1,c)) * 
     >                  (ksize-start(3,c)-end(3,c))

          if (stage .ne. ncells) then

c---------------------------------------------------------------------
c            if this is not the starting cell in this row of cells, 
c            wait for a message to be received, containing the 
c            solution of the previous two stations     
c---------------------------------------------------------------------

             if (timeron) call timer_start(t_ycomm)
             call mpi_irecv(in_buffer, 10*buffer_size, 
     >                      dp_type, successor(2), 
     >                      DEFAULT_TAG, comm_solve, 
     >                      requests(1), error)
             if (timeron) call timer_stop(t_ycomm)


c---------------------------------------------------------------------
c            communication has already been started
c            while waiting, do the block-diagonal inversion for the 
c            cell that was just finished                
c---------------------------------------------------------------------

             call pinvr(slice(2,stage+1))

c---------------------------------------------------------------------
c            wait for pending communication to complete
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_ycomm)
             call mpi_waitall(2, requests, statuses, error)
             if (timeron) call timer_stop(t_ycomm)

c---------------------------------------------------------------------
c            unpack the buffer for the first three factors         
c---------------------------------------------------------------------
             n = 0
             p = 0
             j  = jend
             j1 = j - 1
             do    m = 1, 3
                do   k = start(3,c), ksize-end(3,c)-1
                   do   i = start(1,c), isize-end(1,c)-1
                      sm1 = in_buffer(p+1)
                      sm2 = in_buffer(p+2)
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                        lhs(i,j,k,n+4,c)*sm1 -
     >                        lhs(i,j,k,n+5,c)*sm2
                      rhs(i,j1,k,m,c) = rhs(i,j1,k,m,c) -
     >                        lhs(i,j1,k,n+4,c) * rhs(i,j,k,m,c) - 
     >                        lhs(i,j1,k,n+5,c) * sm1
                      p = p + 2
                   end do
                end do
             end do

c---------------------------------------------------------------------
c            now unpack the buffer for the remaining two factors
c---------------------------------------------------------------------
             do    m = 4, 5
                n = (m-3)*5
                do   k = start(3,c), ksize-end(3,c)-1
                   do   i = start(1,c), isize-end(1,c)-1
                      sm1 = in_buffer(p+1)
                      sm2 = in_buffer(p+2)
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                        lhs(i,j,k,n+4,c)*sm1 -
     >                        lhs(i,j,k,n+5,c)*sm2
                      rhs(i,j1,k,m,c) = rhs(i,j1,k,m,c) -
     >                        lhs(i,j1,k,n+4,c) * rhs(i,j,k,m,c) - 
     >                        lhs(i,j1,k,n+5,c) * sm1
                      p = p + 2
                   end do
                end do
             end do

          else
c---------------------------------------------------------------------
c            now we know this is the first grid block on the back sweep,
c            so we don't need a message to start the substitution. 
c---------------------------------------------------------------------

             j  = jend - 1
             j1 = jend
             n = 0
             do   m = 1, 3
                do   k = start(3,c), ksize-end(3,c)-1
                   do   i = start(1,c), isize-end(1,c)-1
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) -
     >                             lhs(i,j,k,n+4,c)*rhs(i,j1,k,m,c)
                   end do
                end do
             end do

             do    m = 4, 5
                n = (m-3)*5
                do   k = start(3,c), ksize-end(3,c)-1
                   do   i = start(1,c), isize-end(1,c)-1
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) -
     >                             lhs(i,j,k,n+4,c)*rhs(i,j1,k,m,c)
                   end do
                end do
             end do
          endif

c---------------------------------------------------------------------
c         Whether or not this is the last processor, we always have
c         to complete the back-substitution 
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c         The first three factors
c---------------------------------------------------------------------
          n = 0
          do   m = 1, 3
             do   k = start(3,c), ksize-end(3,c)-1
                do   j = jend-2, jstart, -1
                   do    i = start(1,c), isize-end(1,c)-1
                      j1 = j  + 1
                      j2 = j  + 2
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                          lhs(i,j,k,n+4,c)*rhs(i,j1,k,m,c) -
     >                          lhs(i,j,k,n+5,c)*rhs(i,j2,k,m,c)
                   end do
                end do
             end do
          end do

c---------------------------------------------------------------------
c         And the remaining two
c---------------------------------------------------------------------
          do    m = 4, 5
             n = (m-3)*5
             do   k = start(3,c), ksize-end(3,c)-1
                do   j = jend-2, jstart, -1
                   do    i = start(1,c), isize-end(1,c)-1
                      j1 = j  + 1
                      j2 = j1 + 1
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                          lhs(i,j,k,n+4,c)*rhs(i,j1,k,m,c) -
     >                          lhs(i,j,k,n+5,c)*rhs(i,j2,k,m,c)
                   end do
                end do
             end do
          end do

c---------------------------------------------------------------------
c         send on information to the previous processor, if needed
c---------------------------------------------------------------------
          if (stage .ne.  1) then
             j  = jstart
             j1 = jstart + 1
             p = 0
             do    m = 1, 5
                do    k = start(3,c), ksize-end(3,c)-1
                   do    i = start(1,c), isize-end(1,c)-1
                      out_buffer(p+1) = rhs(i,j,k,m,c)
                      out_buffer(p+2) = rhs(i,j1,k,m,c)
                      p = p + 2
                   end do
                end do
             end do

c---------------------------------------------------------------------
c            pack and send the buffer
c---------------------------------------------------------------------

             if (timeron) call timer_start(t_ycomm)
             call mpi_isend(out_buffer, 10*buffer_size, 
     >                     dp_type, predecessor(2), 
     >                     DEFAULT_TAG, comm_solve, 
     >                     requests(2), error)
             if (timeron) call timer_stop(t_ycomm)

          endif

c---------------------------------------------------------------------
c         If this was the last stage, do the block-diagonal inversion          
c---------------------------------------------------------------------
          if (stage .eq. 1) call pinvr(c)

       end do

       if (timeron) call timer_stop(t_ysolve)

       return
       end