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

      subroutine copy_faces

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

c---------------------------------------------------------------------
c     
c This function copies the face values of a variable defined on a set 
c of cells to the overlap locations of the adjacent sets of cells. 
c Because a set of cells interfaces in each direction with exactly one 
c other set, we only need to fill six different buffers. We could try to 
c overlap communication with computation, by computing
c some internal values while communicating boundary values, but this
c adds so much overhead that it's not clearly useful. 
c---------------------------------------------------------------------

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

      integer i, j, k, c, m, requests(0:11), p0, p1, pp,ks,ke,is,ie,je,
     >     p2, p3, p4, p5, b_size(0:5), ss(0:5), js, add,
     >     sr(0:5), error, statuses(MPI_STATUS_SIZE, 0:11)

c---------------------------------------------------------------------
c     exit immediately if there are no faces to be copied           
c---------------------------------------------------------------------
      if (no_nodes .eq. 1) then
         call compute_rhs
         return
      endif

      ss(0) = start_send_east
      ss(1) = start_send_west
      ss(2) = start_send_north
      ss(3) = start_send_south
      ss(4) = start_send_top
      ss(5) = start_send_bottom

      sr(0) = start_recv_east
      sr(1) = start_recv_west
      sr(2) = start_recv_north
      sr(3) = start_recv_south
      sr(4) = start_recv_top
      sr(5) = start_recv_bottom

      b_size(0) = east_size   
      b_size(1) = west_size   
      b_size(2) = north_size  
      b_size(3) = south_size  
      b_size(4) = top_size    
      b_size(5) = bottom_size 

c---------------------------------------------------------------------
c     because the difference stencil for the diagonalized scheme is 
c     orthogonal, we do not have to perform the staged copying of faces, 
c     but can send all face information simultaneously to the neighboring 
c     cells in all directions          
c---------------------------------------------------------------------
      if (timeron) call timer_start(t_bpack)
      p0 = 0
      p1 = 0
      p2 = 0
      p3 = 0
      p4 = 0
      p5 = 0

      do  c = 1, ncells

c---------------------------------------------------------------------
c     fill the buffer to be sent to eastern neighbors (i-dir)
c---------------------------------------------------------------------
         if (cell_coord(1,c) .ne. ncells) then
		    ke=cell_size(3,c)-1
			je=cell_size(2,c)-1
			is=cell_size(1,c)-2
			ie=cell_size(1,c)-1
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))
            do   k = 0, ke
               do   j = 0, je
                  do   i = is, ie
				     pp = p0+k*(je+1)*2*5+j*2*5+(i-is)*5
                     do   m = 1, 5
                        out_buffer(ss(0)+pp+(m-1)) = u(m,i,j,k,c)                        
                     end do					 
                  end do                  
               end do
            end do
!DVM$ end region						
            p0 = p0+(ke+1)*(je+1)*(ie-is+1)*5
         endif

c---------------------------------------------------------------------
c     fill the buffer to be sent to western neighbors 
c---------------------------------------------------------------------
         if (cell_coord(1,c) .ne. 1) then
		    ke=cell_size(3,c)-1
			je=cell_size(2,c)-1
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))
            do   k = 0, ke
               do   j = 0,je
                  do   i = 0,1
				     pp = p1+k*(je+1)*2*5+j*2*5+(i-0)*5				  
                     do   m = 1, 5
                        out_buffer(ss(1)+pp+m-1) = u(m,i,j,k,c)
                     end do
                  end do
               end do
            end do
!DVM$ end region			
			p1=p1+(ke+1)*(je+1)*2*5
         endif

c---------------------------------------------------------------------
c     fill the buffer to be sent to northern neighbors (j_dir)
c---------------------------------------------------------------------
         if (cell_coord(2,c) .ne. ncells) then
		    ke=cell_size(3,c)-1						
			ie=cell_size(1,c)-1
			js=cell_size(2,c)-2
			je=cell_size(2,c)-1
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))			
            do   k = 0,ke
               do   j = js,je
                  do   i = 0,ie
		             pp = p2+k*(ie+1)*2*5+(j-js)*(ie+1)*5+i*5
                     do   m = 1, 5
                        out_buffer(ss(2)+pp+(m-1)) = u(m,i,j,k,c)
                     end do
                  end do
               end do
            end do
!DVM$ end region
			p2=p2+(ke+1)*(je-js+1)*(ie+1)*5
         endif

c---------------------------------------------------------------------
c     fill the buffer to be sent to southern neighbors 
c---------------------------------------------------------------------
         if (cell_coord(2,c).ne. 1) then
		    ke=cell_size(3,c)-1						
			ie=cell_size(1,c)-1			
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))						
            do   k = 0, ke
               do   j = 0, 1
                  do   i = 0, ie
		             pp = p3+k*(ie+1)*2*5+(j-0)*(ie+1)*5 + i*5
                     do   m = 1, 5
                        out_buffer(ss(3)+pp+(m-1)) = u(m,i,j,k,c)
                     end do
                  end do
               end do
            end do
!DVM$ end region	
            p3=p3+(ke+1)*2*(ie+1)*5
         endif

c---------------------------------------------------------------------
c     fill the buffer to be sent to top neighbors (k-dir)
c---------------------------------------------------------------------
         if (cell_coord(3,c) .ne. ncells) then
		    ks=cell_size(3,c)-2
		    ke=cell_size(3,c)-1						
			je=cell_size(2,c)-1		 
			ie=cell_size(1,c)-1		
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))									
            do   k = ks, ke
               do   j = 0, je
                  do   i = 0, ie
				     pp = p4+(k-ks)*(je+1)*(ie+1)*5+j*(ie+1)*5+i*5				  
                     do   m = 1, 5
                        out_buffer(ss(4)+pp+(m-1)) = u(m,i,j,k,c)
                     end do
                  end do
               end do
            end do
!DVM$ end region
            p4=p4+(ke-ks+1)*(je+1)*(ie+1)*5
         endif

c---------------------------------------------------------------------
c     fill the buffer to be sent to bottom neighbors
c---------------------------------------------------------------------
         if (cell_coord(3,c).ne. 1) then
			je=cell_size(2,c)-1		 
			ie=cell_size(1,c)-1		
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))
            do    k=0,1
               do   j = 0, je
                  do   i = 0, ie
				     pp = p5+(k-0)*(je+1)*(ie+1)*5+j*(ie+1)*5+i*5					  
                     do   m = 1, 5
                        out_buffer(ss(5)+pp+(m-1)) = u(m,i,j,k,c)
                     end do
                  end do
               end do
            end do
!DVM$ end region	
            p5=p5+2*(je+1)*(ie+1)*5
         endif

c---------------------------------------------------------------------
c     cell loop
c---------------------------------------------------------------------
      end do
      if (timeron) call timer_stop(t_bpack)

      if (timeron) call timer_start(t_exch)
!DVM$ get_actual(out_buffer)
	  
      call mpi_irecv(in_buffer(sr(0)), b_size(0), 
     >     dp_type, successor(1), WEST,  
     >     comm_rhs, requests(0), error)
      call mpi_irecv(in_buffer(sr(1)), b_size(1), 
     >     dp_type, predecessor(1), EAST,  
     >     comm_rhs, requests(1), error)
      call mpi_irecv(in_buffer(sr(2)), b_size(2), 
     >     dp_type, successor(2), SOUTH, 
     >     comm_rhs, requests(2), error)
      call mpi_irecv(in_buffer(sr(3)), b_size(3), 
     >     dp_type, predecessor(2), NORTH, 
     >     comm_rhs, requests(3), error)
      call mpi_irecv(in_buffer(sr(4)), b_size(4), 
     >     dp_type, successor(3), BOTTOM,
     >     comm_rhs, requests(4), error)
      call mpi_irecv(in_buffer(sr(5)), b_size(5), 
     >     dp_type, predecessor(3), TOP,   
     >     comm_rhs, requests(5), error)

      call mpi_isend(out_buffer(ss(0)), b_size(0), 
     >     dp_type, successor(1),   EAST, 
     >     comm_rhs, requests(6), error)
      call mpi_isend(out_buffer(ss(1)), b_size(1), 
     >     dp_type, predecessor(1), WEST, 
     >     comm_rhs, requests(7), error)
      call mpi_isend(out_buffer(ss(2)), b_size(2), 
     >     dp_type,successor(2),   NORTH, 
     >     comm_rhs, requests(8), error)
      call mpi_isend(out_buffer(ss(3)), b_size(3), 
     >     dp_type,predecessor(2), SOUTH, 
     >     comm_rhs, requests(9), error)
      call mpi_isend(out_buffer(ss(4)), b_size(4), 
     >     dp_type,successor(3),   TOP, 
     >     comm_rhs,   requests(10), error)
      call mpi_isend(out_buffer(ss(5)), b_size(5), 
     >     dp_type,predecessor(3), BOTTOM, 
     >     comm_rhs,requests(11), error)


      call mpi_waitall(12, requests, statuses, error)
      if (timeron) call timer_stop(t_exch)

c---------------------------------------------------------------------
c     unpack the data that has just been received;             
c---------------------------------------------------------------------
      if (timeron) call timer_start(t_bpack)
      p0 = 0
      p1 = 0
      p2 = 0
      p3 = 0
      p4 = 0
      p5 = 0
!DVM$ actual(in_buffer)

      do   c = 1, ncells

         if (cell_coord(1,c) .ne. 1) then
		    ke=cell_size(3,c)-1
			je=cell_size(2,c)-1
			
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))					 
            do   k = 0, ke
               do   j = 0, je
                  do   i = -2, -1
				     pp = p0+k*(je+1)*2*5+j*2*5+(i+2)*5
                     do   m = 1, 5
                        u(m,i,j,k,c) = in_buffer(sr(1)+pp+(m-1))
                     end do
                  end do
               end do
            end do
!DVM$ end region
            p0=p0+(ke+1)*(je+1)*2*5
         endif

         if (cell_coord(1,c) .ne. ncells) then
		    ke=cell_size(3,c)-1
			je=cell_size(2,c)-1
			ie=cell_size(1,c)+1
			is=cell_size(1,c)
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))					 
            do  k = 0, ke
               do  j = 0, je
                  do  i = is, ie
				     pp = p1+k*(je+1)*2*5+j*2*5+(i-is)*5
                     do   m = 1, 5
                        u(m,i,j,k,c) = in_buffer(sr(0)+pp+(m-1))
                     end do
                  end do
               end do
            end do
!DVM$ end region	
            p1=p1+(ke+1)*(je+1)*2*5
         end if
            
         if (cell_coord(2,c) .ne. 1) then
		    ke=cell_size(3,c)-1						
			ie=cell_size(1,c)-1
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))									
            do  k = 0, ke
               do   j = -2, -1
                  do  i = 0, ie
		             pp = p2+k*(ie+1)*2*5+(j+2)*(ie+1)*5+i*5
                     do   m = 1, 5
                        u(m,i,j,k,c) = in_buffer(sr(3)+pp+(m-1))
                     end do
                  end do
               end do
            end do
!DVM$ end region
            p2=p2+(ke+1)*2*(ie+1)*5
         endif
            
         if (cell_coord(2,c) .ne. ncells) then
		    ke=cell_size(3,c)-1						
			ie=cell_size(1,c)-1			
			js=cell_size(2,c)
			je=cell_size(2,c)+1
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))
            do  k = 0, ke
               do   j = js, je
                  do  i = 0, ie
		             pp = p3+k*(ie+1)*2*5+(j-js)*(ie+1)*5+i*5
                     do   m = 1, 5
                        u(m,i,j,k,c) = in_buffer(sr(2)+pp+(m-1))
                     end do
                  end do
               end do
            end do
!DVM$ end region	
            p3=p3+(ke+1)*2*(ie+1)*5
         endif

         if (cell_coord(3,c) .ne. 1) then
			je=cell_size(2,c)-1		 
			ie=cell_size(1,c)-1	
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))				
            do  k = -2, -1
               do  j = 0, je
                  do  i = 0, ie
					pp = p4+(k+2)*(je+1)*(ie+1)*5+j*(ie+1)*5+i*5	
                     do   m = 1, 5
                        u(m,i,j,k,c) = in_buffer(sr(5)+pp+(m-1))
                     end do
                  end do
               end do
            end do
!DVM$ end region	
            p4=p4+2*(je+1)*(ie+1)*5
         endif

         if (cell_coord(3,c) .ne. ncells) then
			je=cell_size(2,c)-1		 
			ie=cell_size(1,c)-1	
            ks=cell_size(3,c)
			ke=cell_size(3,c)+1
!DVM$ region			
!DVM$ PARALLEL(k,j,i), PRIVATE(pp,m), TIE(u(*,i,j,k,*))			
            do  k = ks, ke
               do  j = 0, je
                  do  i = 0, ie
				     pp=p5+(k-ks)*(je+1)*(ie+1)*5+j*(ie+1)*5+i*5					  
                     do   m = 1, 5
                        u(m,i,j,k,c) = in_buffer(sr(4)+pp+(m-1))
                     end do
                  end do
               end do
            end do
!DVM$ end region
            p5=p5+2*(je+1)*(ie+1)*5
         endif

c---------------------------------------------------------------------
c     cells loop
c---------------------------------------------------------------------
      end do
      if (timeron) call timer_stop(t_bpack)

c---------------------------------------------------------------------
c     do the rest of the rhs that uses the copied face values          
c---------------------------------------------------------------------
      call compute_rhs

      return
      end
