!---------------------------------------------------------------------------
! pmlib_mesh_map_interpolation.f
!   version: 2D single-core
!---------------------------------------------------------------------------
!                  __    __
!    ___  _ _     / /__ / /
!   | _ \/ ` \ _ / // // _ \
!   | ,_/_,_,_\  \_/\_/\___/
!   |_|                      Mads M. Hejlesen* and Jens H. Walther*^
!                            * Technical University of Denmark
!                            ^ ETH-Zurich
!
!---------------------------------------------------------------------------
! Copyright (c) 2014 Technical University of Denmark, 
!                    Department of Mechanical Engineering
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program; if not, see <http://www.gnu.org/licenses/>
!
! Contact: mmhej@mek.dtu.dk or jhwa@mek.dtu.dk
!---------------------------------------------------------------------------
SUBROUTINE pmlib_mesh_map_interpolation(topo,dir,nvar_send,ierr)

!---------------------------------------------------------------------------
! Sets up the mesh boundary communication for superposition 
! of the interpolated mesh
!  - Finds the overlapping mesh indices
!  - Finds the communication sequence
!  - Allocates the send and recieve buffers
!
! The communication is calculated on processor 0 only
! after which it is broadcasted to the remaining processors.
!---------------------------------------------------------------------------

USE pmlib_mod_communication
USE pmlib_mod_topology

IMPLICIT NONE
!---------------------------------------------------------------------------
! Arguments
!---------------------------------------------------------------------------
  TYPE(class_topology)                      :: topo
  INTEGER, INTENT(IN)                       :: dir
  INTEGER, INTENT(IN)                       :: nvar_send
  INTEGER, INTENT(OUT)                      :: ierr

!---------------------------------------------------------------------------
! Local variables
!---------------------------------------------------------------------------
  CHARACTER(LEN=22)                     :: caller = 'pmlib_comm_boundaries'
  INTEGER                               :: ncomm, icomm
  INTEGER                               :: iproc, jproc
  INTEGER                               :: i,j,k

  INTEGER                               :: ineigh, neigh
  REAL(MK), DIMENSION(pmlib_ndim)       :: xmin_send, xmax_send
  REAL(MK), DIMENSION(pmlib_ndim)       :: xmin_recieve, xmax_recieve
  REAL(MK), DIMENSION(pmlib_ndim)       :: dx

  INTEGER,DIMENSION(2*pmlib_ndim)       :: nghost_recieve, nghost_send
  INTEGER, DIMENSION(pmlib_ndim)        :: ncell_send, ncell_recieve
  INTEGER                               :: nghost, ncell

!---------------------------------------------------------------------------
! Initiate subroutine
!---------------------------------------------------------------------------
  ierr = 0

!---------------------------------------------------------------------------
! Find number of communications
!---------------------------------------------------------------------------
  icomm = 0
  DO ineigh = 2*dir-1,2*dir
    
    neigh = topo%ineigh(ineigh)

    IF( neigh .NE. -1 )THEN
      nghost_send    = topo%nghost
      nghost_recieve = topo%nghost

      IF( ( nghost_send(ineigh) .GT. 0 ) .OR. &
        & ( ineigh .EQ. 2*dir-1 .AND.  nghost_recieve(2*dir) .GT. 0 ) .OR. &
        & ( ineigh .EQ. 2*dir   .AND. nghost_recieve(2*dir-1) .GT. 0 ) )THEN
        icomm = icomm + 1
      END IF
    END IF ! neighbours

  END DO !ineigh

!---------------------------------------------------------------------------
! Allocate communication lists
!---------------------------------------------------------------------------
  ncomm               = icomm
  communication%ncomm = ncomm
  communication%ctype = 2

  ALLOCATE( communication%imin_send(pmlib_ndim,ncomm) )
  ALLOCATE( communication%imax_send(pmlib_ndim,ncomm) )
  ALLOCATE( communication%imin_recieve(pmlib_ndim,ncomm) )
  ALLOCATE( communication%imax_recieve(pmlib_ndim,ncomm) )

  ALLOCATE( communication%ncell_send(pmlib_ndim,ncomm) )
  ALLOCATE( communication%ncell_recieve(pmlib_ndim,ncomm) )

!---------------------------------------------------------------------------
! Loop through neighbours
!---------------------------------------------------------------------------
  icomm = 0
  DO ineigh = 2*dir-1,2*dir
    
    neigh = topo%ineigh(ineigh)

    IF( neigh .NE. -1 )THEN
      nghost_send    = topo%nghost
      nghost_recieve = topo%nghost

      IF( ( nghost_send(ineigh) .GT. 0 ) .OR. &
        & ( ineigh .EQ. 2*dir-1 .AND.   nghost_recieve(2*dir) .GT. 0 ) .OR. &
        & ( ineigh .EQ. 2*dir   .AND. nghost_recieve(2*dir-1) .GT. 0 ) )THEN

!---------------------------------------------------------------------------
! Count communications
!---------------------------------------------------------------------------
        icomm = icomm + 1

!---------------------------------------------------------------------------
! Store number of mesh cells in the topologies
!---------------------------------------------------------------------------
        ncell_send     = topo%ncell
        ncell_recieve  = topo%ncell

!---------------------------------------------------------------------------
! Find indices of ghost cells
!---------------------------------------------------------------------------
        DO i = 1,pmlib_ndim

          IF(2*i-1 .EQ. ineigh)THEN

            communication%imin_send(i,icomm)    = 1 - nghost_send(2*i-1)
            communication%imin_recieve(i,icomm) = ncell_recieve(i) + 1 &
                                              & - nghost_send(2*i-1)

            communication%imax_send(i,icomm)    = nghost_recieve(2*i)
            communication%imax_recieve(i,icomm) = ncell_recieve(i) &
                                              & + nghost_recieve(2*i)

          ELSEIF(2*i .EQ. ineigh)THEN

            communication%imin_send(i,icomm)    = ncell_send(i) &
                                              & - nghost_recieve(2*i-1) + 1
            communication%imin_recieve(i,icomm) = 1 - nghost_recieve(2*i-1)

            communication%imax_send(i,icomm)    = ncell_send(i) &
                                              & + nghost_send(2*i)
            communication%imax_recieve(i,icomm) = nghost_send(2*i)

          ELSE
            ncell  = MIN(ncell_send(i),ncell_recieve(i))

            nghost = MIN(nghost_send(2*i-1),nghost_recieve(2*i-1))
            communication%imin_send(i,icomm)    = 1 - nghost
            communication%imin_recieve(i,icomm) = 1 - nghost

            nghost = MIN(nghost_send(2*i),nghost_recieve(2*i))
            communication%imax_send(i,icomm)    = ncell + nghost
            communication%imax_recieve(i,icomm) = ncell + nghost
          END IF

          communication%ncell_send(i,icomm) = &
                   & communication%imax_send(i,icomm) - &
                   &  communication%imin_send(i,icomm) + 1
      
          communication%ncell_recieve(i,icomm) = &
                   & communication%imax_recieve(i,icomm) - &
                   &  communication%imin_recieve(i,icomm) + 1

! check if the number of sended and recieved cells coinside
          IF(  communication%ncell_send(i,icomm) .NE. &
            &  communication%ncell_recieve(i,icomm))THEN
            ierr = -1
            CALL pmlib_write( caller, &
                & 'Number of sended and recieved cells does not correspond')
            GOTO 9999
          END IF

        END DO

      END IF ! nghost
    END IF ! neighbours

  END DO !ineigh

!---------------------------------------------------------------------------
! Check number of communications
!---------------------------------------------------------------------------
  IF(icomm .NE. ncomm)THEN
    ierr = -1
    CALL pmlib_write(caller,'Number of communications are inconsistent')
    GOTO 9999
  ENDIF

!---------------------------------------------------------------------------
! Allocate communication buffer
!---------------------------------------------------------------------------
  IF( ASSOCIATED(comm_buffer) ) THEN
    DEALLOCATE(comm_buffer)
    IF (ierr .NE. 0) THEN
      CALL pmlib_write(caller,'Failed to deallocate array.')
      GOTO 9999
    ENDIF
  END IF 
  ALLOCATE( comm_buffer(ncomm) )

  DO j = 1,ncomm

    ncell_send(1) = communication%ncell_send(1,j)  
    ncell_send(2) = communication%ncell_send(2,j)
 
! Allocate send array for communication j
    ALLOCATE( comm_buffer(j)%mesh_send(nvar_send, &
            & ncell_send(1),ncell_send(2)) )
    comm_buffer(j)%nsend = nvar_send*ncell_send(1)*ncell_send(2)
    comm_buffer(j)%ncell_send = ncell_send

! Allocate recieve array for communication j
    ALLOCATE( comm_buffer(j)%mesh_recieve(nvar_send, &
            & ncell_send(1),ncell_send(2)) )
    comm_buffer(j)%nrecieve = nvar_send*ncell_send(1)*ncell_send(2)
    comm_buffer(j)%ncell_recieve = ncell_send

  END DO

!---------------------------------------------------------------------------
! set counter for number of packed variables
!---------------------------------------------------------------------------
  communication%npacked = 0

!---------------------------------------------------------------------------
! Return
!---------------------------------------------------------------------------
 9999 CONTINUE
RETURN
END SUBROUTINE pmlib_mesh_map_interpolation

