couple_mod.F

Go to the documentation of this file.

      MODULE couple_mod
!******************************************************************
! THAT MODULE CONTAINS THE SUBROUTINES NEEDED TO EXCHANGE VALUES OF VARIABLE
! BETWEEN TWO MODULES (TILL NOW TELEMAC2D AND TOMAWAC) WHEN THEY DEAL
! WITH DIFFERENT MESH.
! ASSUMPTIONS:
! SAME NUMBER OF PROCESSES ON BOTH COUPLED MODELS
!******************************************************************

      USE bief
      USE declarations_special
      USE declarations_parallel
      USE interface_parallel

      TYPE couple_data_type
        ! RECEIVER DATA

        ! SENDER DATA
        INTEGER, POINTER, DIMENSION(:,:) :: send_list
        INTEGER, POINTER, DIMENSION(:)   :: npoin_all

        ! LIST WITH FOR EACH PROCESS
        ! LOCAL NODE ID WHERE TO GET DATA FROM FOR SENDING
        INTEGER, POINTER, DIMENSION(:,:)   :: send_map

        ! LIST WITH FOR EACH PROCESS
        ! LOCAL NODE ID WHERE TO PUT DATA THAT IS RECEIVED
        ! THIRD DIMENSION IS USED TO HAVE INFLUENCE OF MULTIPLE POINTS.
        ! E.G. FOR INTERPOLATION
        INTEGER, POINTER, DIMENSION(:,:,:)   :: recv_map

        !NUMBER OF POINTS TO SEND TO OTHER PROCESSES
        INTEGER, POINTER, DIMENSION(:)   :: poin2send
        !NUMBER OF POINTS TO SRECEIVE FROM OTHER PROCESSES
        INTEGER, POINTER, DIMENSION(:)   :: poin2recv

        ! MAXIMUM NUMBER OF POINTS ON ANY PROCESS
        INTEGER:: maxpoin
        ! NUMBER OF POINTS ON EACH PROCESS SENDING DATA
        INTEGER, POINTER, DIMENSION(:)   :: send_npoin
        !GLOBAL NODE NUMBER OF THE SENDER
        INTEGER, POINTER, DIMENSION(:)   :: send_knolg

        !NUMBER OF PROCCESS TO SEND TO (INCLUDING LOCAL)
        INTEGER :: np_send
        !NUMBRE OF PROCESSES TO RECEIVE FROM (INCLUDING LOCAL)
        INTEGER :: np_recv
        !NUMBER OF PROCCESS TO SEND TO (EXCLUDING LOCAL)
        INTEGER :: np_send_noloc
        !NUMBRE OF PROCESSES TO RECEIVE FROM (EXCLUDING LOCAL)
        INTEGER :: np_recv_noloc

        !BUFFER FOR SENDING DATA; FIRST DIMENSION: ID; SECOND DIMENSION: PROCESS
        DOUBLE PRECISION, POINTER, DIMENSION(:,:)   :: send_buf
        !BUFFER FOR RECEIVING DATA; FIRST DIMENSION: ID; SECOND DIMENSION: PROCESS
        DOUBLE PRECISION, POINTER, DIMENSION(:,:)   :: recv_buf
        !LIST OF MPI_REQUEST FOR MPI WAIT_ALL
        INTEGER, POINTER, DIMENSION(:) :: reqlist

        !MAXIMUM NUMBER OF NODES USED TO CALCULATE VALUES (USING INTERPOLATION)
        INTEGER :: npp
        !INTERPOLATION WEIGHTS
        DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: weights

      END TYPE

      !RECV_LIST: NPX2 LIST WITH IN THE FIRST COLUMN THE LOCAL POINT NUMBER
      ! ON THE SENDER AND ON THE SECOND COLUMN THE PROCESS ID

      type(couple_data_type), ALLOCATABLE, DIMENSION(:) :: cpld
      INTEGER :: nrcouple
! MAXIMUM OF MODULES COUPLING 2 ARE NEEDED FOR THE MOMENT
      INTEGER, PARAMETER :: maxcouple = 4
! MAXIMUM OF VARIABLES TO EXCHANGE
      INTEGER, PARAMETER :: maxvar    = 10
! MAXIMUM OF POINTS TO INTERPOLATE ON A POINT ON THE OTHER MESH.
      INTEGER, PARAMETER :: maxnpp    = 20

      LOGICAL :: do_couple = .false.

      SAVE
      PRIVATE

      PUBLIC :: init_couple,add_sender, add_receiver,
     &    send_couple, receive_couple, end_couple, do_couple

      CONTAINS

!                   **********************
                    SUBROUTINE init_couple
!                   **********************
     &()
!
!***********************************************************************
!BRIEF    CONSTRUCTOR FOR COUPLING ROUTINE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      INTEGER :: ICOUPLE, IERR
!
!-----------------------------------------------------------------------
!
!INITIALIZE COUPLING DATA

      do_couple = .true.
      nrcouple =  0
      ALLOCATE(cpld(maxcouple),stat=ierr)
      CALL check_allocate(ierr,'CPLD')

      !TODO CHECK IF NEEDED
      DO icouple=1,maxcouple
        ALLOCATE(cpld(icouple)%SEND_NPOIN(max(ncsize,1)),stat=ierr)
        CALL check_allocate(ierr,'SEND_NPOIN')
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE

!                   **********************
                    SUBROUTINE  add_sender
!                   **********************
     &(mesh_send,cid)
!
!***********************************************************************
!
!BRIEF    .
!
!NOTE: SENDER MUST BE ADDED BEFORE RECIEVER IS ADDED
!

!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CID          |-->| ID OF THE COUPLING (STARTS AT 1)
!| MESH_SEND    |-->| MESH (LOCAL) OF THE SENDER
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)   :: CID
      type(bief_mesh), INTENT(IN)   :: mesh_send
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

      INTEGER :: NPOIN, IERR , IPOIN, IPROC, ITAG, IREQ
      INTEGER , DIMENSION(1) :: SEND_BUF
      INTEGER , DIMENSION(:), ALLOCATABLE :: RECV_BUF
!
!-----------------------------------------------------------------------
!
      ALLOCATE(cpld(cid)%REQLIST(max(ncsize,1)-1),stat=ierr)
      CALL check_allocate(ierr,'REQLIST')

      ALLOCATE (recv_buf(max(ncsize,1)),stat=ierr)
      CALL check_allocate(ierr,'RECV_BUF')
      recv_buf = 0
      !COUNT NUMBER OF COUPLINGS
      nrcouple = nrcouple +1
      npoin = mesh_send%NPOIN

      !ADD NUMBER OF POINTS OF SENDER TO ALL PROCESSES
      send_buf(1) = npoin

#if defined (HAVE_MPI)
      IF (ncsize.GT.1) THEN
        CALL mpi_allgather(send_buf, 1, mpi_integer, recv_buf, 1,
     &   mpi_integer, comm, ierr)
      ELSE
          recv_buf(1) = send_buf(1)
      ENDIF
#else
      recv_buf(1) = send_buf(1)
#endif

      IF (ierr.NE.0) THEN
        WRITE (lu,*) 'ALL TO ALL PROBLEM IN ADD_SENDER'
        CALL plante(1)
      ENDIF


      IF (ncsize.GT.1) THEN
        CALL p_sync
      ENDIF

      cpld(cid)%SEND_NPOIN(1:max(ncsize,1)) = recv_buf(1:max(ncsize,1))
      DEALLOCATE (recv_buf)
      cpld(cid)%MAXPOIN = maxval(cpld(cid)%SEND_NPOIN(1:max(ncsize,1)))

      ! SAVE GLOBAL POINT NUMBERS

      ALLOCATE(cpld(cid)%SEND_KNOLG(npoin),stat=ierr)
      CALL check_allocate(ierr,'SEND_KNOLG')
      IF (ncsize.GT.1) THEN
        DO ipoin = 1,npoin
          cpld(cid)%SEND_KNOLG(ipoin)=mesh_send%KNOLG%I(ipoin)
        ENDDO
      ELSE
        ! NOT SURE IF THIS IS NEEDED
        DO ipoin = 1,npoin
          cpld(cid)%SEND_KNOLG(ipoin)= ipoin
        ENDDO
      ENDIF

      ! SEND GLOBAL POINT NUMBERS TO ALL OTHER PROCS
      !  ALREADY WORKS ON 1 PROCESSOR
      ireq = 0
      DO iproc = 0,max(ncsize,1)-1
        IF (iproc.EQ.ipid) THEN
          cycle
        ENDIF
        itag = ipid+iproc*ncsize
        ireq = ireq+1
        CALL p_iwriti(cpld(cid)%SEND_KNOLG,npoin,
     &      iproc,itag,cpld(cid)%REQLIST(ireq))
      ENDDO

      RETURN
      END SUBROUTINE

!-----------------------------------------------------------------------
!
!                   ***********************
                    SUBROUTINE read_recv
!                   ***********************
     &        (file_recv,recv_var,cid,npoin,send_node,npp)
!
!***********************************************************************
!
!BRIEF    READS THE VARIABLE(S) IN THE SELAFIN CONTAINING THE MAPPING
!         BETWEEN TWO COUPLED MODELS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         ,  INTENT(IN) :: CID, NPOIN
      TYPE(bief_file) ,  INTENT(IN) :: FILE_RECV
      CHARACTER(LEN=16), INTENT(IN) :: RECV_VAR
      DOUBLE PRECISION,  INTENT(OUT), DIMENSION(:,:), ALLOCATABLE ::
     & send_node
      INTEGER,  INTENT(INOUT) :: NPP
      CHARACTER(LEN=2)  :: STRI
      CHARACTER(LEN=16), ALLOCATABLE :: VAR_LIST(:)
      CHARACTER(LEN=16), ALLOCATABLE :: UNIT_LIST(:)
      CHARACTER(LEN=16) :: VARNAME
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      INTEGER :: IERR,  NVAR, IVAR, INUM, I, INDVAR, IPOIN, NRFOUND
      DOUBLE PRECISION :: TMP

      ! DETERMINE NUMBER OF POINTS USED IN INTERPOLATIONS   (FROM FIND_VARIABLE)

      ! READ VARIABLE NAMES FROM SLF FILE
      CALL get_data_nvar(file_recv%FMT,file_recv%LU,nvar,ierr)
      CALL check_call(ierr, 'FIND_VARIABLE:GET_DATA_NVAR')
      ALLOCATE(var_list(nvar),stat=ierr)
      ALLOCATE(unit_list(nvar),stat=ierr)
      CALL check_allocate(ierr, 'VAR_LIST')
      CALL get_data_var_list(file_recv%FMT,file_recv%LU,nvar,
     &                   var_list,unit_list,ierr)
      CALL check_call(ierr, 'FIND_VARIABLE:GET_DATA_VAR_LIST')

      ! LOOK FOR SPACE IN VARIABLE NAME
      inum = 0
      DO i=1,16
        IF (recv_var(i:i).EQ.' ') THEN
            inum = i
            EXIT
        ENDIF
      ENDDO
      IF (inum.EQ.0) THEN
        WRITE (lu,*) 'INVALID VARIABLE NAME FOR COUPLING: ', recv_var
        CALL plante (1)
      ENDIF
      nrfound = 0
      DO ivar = 1,nvar
        IF (var_list(ivar)(1:inum-1).EQ.recv_var(1:inum-1).AND.
     &    var_list(ivar)(inum:inum).NE. 'W')  THEN
          nrfound = nrfound + 1
        ENDIF
      ENDDO
      IF (nrfound.EQ.0) THEN
        WRITE (lu,*) 'ERROR: COUPLE_MOD NO COUPLING INFO FOR ',
     &        recv_var
        CALL plante(1)
      ENDIF


      ! LOOK FOR THE HIGHEST VARIABLE NUMBER.
      ! IT IS ASSUMED THAT VARIABLES ARE NUMBERED CONTINUOUSLY
      ! THIS IS CHECKED LATER WHEN READING THE VARIABLES
      npp = 0
      DO ivar = 1,nvar
        IF (var_list(ivar)(1:inum-1).EQ.recv_var(1:inum-1).AND.
     &     var_list(ivar)(inum:inum).NE. 'W')  THEN
          READ (var_list(ivar)(inum:inum+1),*) indvar
          IF (indvar.GT.maxnpp .OR. indvar.LE.0) THEN
            WRITE (lu,*) 'INVALID NUMBER FOR ',recv_var,': ',indvar
            CALL plante (1)
          ENDIF
          npp = max(npp,indvar)
        ENDIF
      ENDDO

      DEALLOCATE(var_list)
      DEALLOCATE(unit_list)

!----------------------------------------------------------------------
! ALLOCATE MEMORY FOR WEIGHTS AND SEND NODES

      ALLOCATE(send_node(npoin,npp),stat=ierr)
      CALL check_allocate(ierr,'SEND_NODE')

      ALLOCATE(cpld(cid)%WEIGHTS(npoin,npp),stat=ierr)
      CALL check_allocate(ierr,'WEIGHTS')


      ! LOOP OVER ALL NEEDED VARIABLES
      DO i=1,npp
        WRITE(stri,'(I2.2)') i
        varname = recv_var(1:inum-1)//stri
      ! READ  ID FROM WHERE TO RECEIVE DATA (GLOBAL ID OF THE SENDER)
        CALL find_variable(file_recv%FMT,file_recv%LU,varname,
     &      send_node(:,i),npoin,ierr)
        IF (ierr.NE.0) THEN
          WRITE (lu,*) 'VARIABLE ',varname, ' CANNOT BE READ '
          CALL plante(1)
        ENDIF

        ! READ INTERPOLATION WEIGHTS
        varname = recv_var(1:inum-1)//'WTS'//stri
        CALL find_variable(file_recv%FMT,file_recv%LU,varname,
     &      cpld(cid)%WEIGHTS(:,i),npoin,ierr)
        IF (ierr.NE.0) THEN
          WRITE (lu,*) 'VARIABLE ',varname, ' CANNOT BE READ '
          CALL plante(1)
        ENDIF
      ENDDO
      ! CHECK WEIGHTS SUM TO ONE IN CASE NO DEFAULT IS USED
      DO ipoin=1,npoin
        IF (any(send_node(ipoin,:).GT.0)) THEN
          tmp = sum(cpld(cid)%WEIGHTS(ipoin,:))
      ! NOTE THRESHOLD IS HIGH. PROBABLY DUE TO SINGLE PRECISION IN SLF FILE
          IF (abs(tmp-1.0d0).GT.1.0d-4) THEN
            WRITE (lu,*) 'WRONG WEIGHTS IN POINT ', npoin
            WRITE (lu,*) 'WEIGHST ARE: ', cpld(cid)%WEIGHTS(ipoin,:)
            CALL plante(1)
          ELSE ! MAKE SURE THEY ARE ALSO ACCURATE IN DOUBLE PRECISION
            cpld(cid)%WEIGHTS(ipoin,:)=cpld(cid)%WEIGHTS(ipoin,:)/tmp
          ENDIF
        ENDIF
      ENDDO

      RETURN
      END SUBROUTINE
!
!-----------------------------------------------------------------------
!
!                   ***********************
                    SUBROUTINE add_receiver
!                   ***********************
     &         (mesh_recv,file_recv,recv_var,cid)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| MESH_RECV         |-->| MESH (LOCAL) OF THE RECEIVER
!| FILE_RECV         |-->| FILE WITH INFORMATION FROM WHERE TO RECEIVE DATA
!| RECV_VAR          |-->| NAME OF VARIABLE  WHICH CONTAINS THE NODE NR OF THE SENDER
!| CID               |-->| ID OF THE COUPLING (STARTS AT 1)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         ,  INTENT(IN) :: CID
      TYPE(bief_file) ,  INTENT(IN) :: FILE_RECV
      TYPE(bief_mesh),   INTENT(IN) :: MESH_RECV
      CHARACTER(LEN=16), INTENT(IN) :: RECV_VAR
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      INTEGER :: IERR, NPOIN , IPROC, IP_S, IP_R, IP_P
      INTEGER :: MAXP, IPOIN, NP, ITMP, ITAG , IREQ, NPR,NPS, NPP

      INTEGER, ALLOCATABLE, DIMENSION(:) :: NSEND

      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: SEND_NODE
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: TMP_SEND, LOC_SEND_MAP
      LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: DO_SEND
      INTEGER, ALLOCATABLE, DIMENSION(:) :: TMP_KNOLG
      LOGICAL, ALLOCATABLE, DIMENSION(:) :: TMP_LOG

      ! LIST WITH FOR EACH POINT, FROM WHERE TO RECEIVE DATA
      !1.) LOCAL NODE NR
      !2.) PROCES ID
      INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: RECV_LIST
!
!-----------------------------------------------------------------------
!

      npoin = mesh_recv%NPOIN
      maxp  = cpld(cid)%MAXPOIN

      ! READ WEIGHTS AND NODEINDICES FROM SELAFIN FILE
      CALL read_recv(file_recv,recv_var,cid,npoin,send_node,npp)
      cpld(cid)%NPP = npp

      !ALLOCATE DATA
      ALLOCATE(recv_list(npoin,2,npp),stat=ierr)
      CALL check_allocate(ierr,'RECV_LIST')
      recv_list = 0

      ALLOCATE(tmp_knolg(maxp),stat=ierr)
      CALL check_allocate(ierr,'TMP_KNOLG')

!     FIND CORRESPONDING LOCAL NUMBERS ON SENDER
!     (ONLY WORKS IF SENDER IS ALREADY ADDED)

      recv_list = 0
      DO iproc = 0,max(ncsize,1)-1
      ! RECEIVE KNOLG FROM SENDER
      ! BLOCKING RECEIVE TO LIMIT AMOUNT OF DATA
        np = cpld(cid)%SEND_NPOIN(iproc+1)
        IF (iproc.NE.ipid) THEN
          itag = ipid*ncsize+iproc
          !DEFACTO BLOCKING RECEIVE
          CALL  p_ireadi(tmp_knolg,np,
     &         iproc,itag,cpld(cid)%REQLIST(1))
          CALL p_wait_paraco(cpld(cid)%REQLIST(1:1),1)
        ELSE
      ! DATA IS ALREADY ON THE SAME PROCESS
          tmp_knolg(1:np) = cpld(cid)%SEND_KNOLG(1:np)
        ENDIF

        !FIND MATCHING POINTS
        DO  ip_p = 1,npp
          DO ip_r = 1,npoin
            DO ip_s = 1,np
              IF (nint(send_node(ip_r,ip_p)).EQ.tmp_knolg(ip_s)) THEN
                recv_list(ip_r,1,ip_p) = ip_s
                recv_list(ip_r,2,ip_p) = iproc+1
                EXIT
              ENDIF
            ENDDO  !IP_R
          ENDDO  !IP_S
        ENDDO   ! IP_P
      ENDDO ! IPROC

!----------------------------------------------------------------------
! BUILD UP SEND LIST AND RECEIVE LIST

      !COMMUNICATED DATA
      ALLOCATE(cpld(cid)%SEND_MAP(maxp,max(ncsize,1)),stat=ierr)
      CALL check_allocate(ierr,'SEND_MAP')
      ALLOCATE(loc_send_map(maxp,max(ncsize,1)),stat=ierr)
      CALL check_allocate(ierr,'_LOC_SEND_MAP')
      loc_send_map = 0

      ALLOCATE(cpld(cid)%RECV_MAP(max(ncsize,1),npoin,npp),stat=ierr)
      CALL check_allocate(ierr,'RECV_MAP')

      !NUMBER OF POINTS TO COMMUNICATE
      ALLOCATE(cpld(cid)%POIN2SEND(max(ncsize,1)),stat=ierr)
      CALL check_allocate(ierr,'POIN2SEND')
      ALLOCATE(cpld(cid)%POIN2RECV(max(ncsize,1)),stat=ierr)
      CALL check_allocate(ierr,'POIN2RECV')


      ALLOCATE(do_send(max(ncsize,1),maxp),stat=ierr)
      do_send = .false.
      ALLOCATE(tmp_send(max(ncsize,1),maxp),stat=ierr)
      tmp_send =0

      ALLOCATE(nsend(max(ncsize,1)),stat=ierr)
      CALL check_allocate(ierr,'NSEND')
      nsend = 0

      cpld(cid)%POIN2SEND = 0
      cpld(cid)%POIN2RECV = 0
      cpld(cid)%RECV_MAP = 0
      cpld(cid)%SEND_MAP = 0

!----------------------------------------------------------------------
      ! STORE NUMBER OF POINTS
      ! THIS VARIABLE MAPS NUMBER OF POINTS ON SENDER TO NUMBER OF POINTS ON RECEIVER

      !LOOK FOR EACH POINT WHERE IT COMES FROM
      DO ip_p = 1,npp
        DO ip_r=1,npoin
        ! VALUES FROM GLOBAL LIST
          ip_s  = recv_list(ip_r,1,ip_p)
          iproc = recv_list(ip_r,2,ip_p)
          IF (ip_s.GT.0) THEN
            IF (.NOT.do_send(iproc,ip_s)) THEN
          !POINT NOT YET SEND TO THIS PROCESS. ADD TO THE LISTS
              do_send(iproc,ip_s) = .true.
              nsend(iproc)  = nsend(iproc)+1
              cpld(cid)%RECV_MAP(iproc,ip_r,ip_p) = nsend(iproc)
              ! NOTE ORDER OF ARGUMENTS IS OPPOSITE IN LOC_SEND_MAP
              ! THIS IS  NEEDED FOR MPI_SCATTER LATER
              loc_send_map(nsend(iproc),iproc) = ip_s
              tmp_send(iproc,ip_s) = nsend(iproc)
            ELSE
              ! POINT IS ALREADY ON THE LIST. LOOK UP ITS ID
              cpld(cid)%RECV_MAP(iproc,ip_r,ip_p) = tmp_send(iproc,ip_s)
            ENDIF
          ENDIF
        ENDDO !IP_R
      ENDDO !IP_P

      ! STORE NUMBER OF POINTS
      DO iproc = 1,max(ncsize,1)
        cpld(cid)%POIN2RECV(iproc) = nsend(iproc)
      ENDDO

!----------------------------------------------------------------------

      ! COMMUNICATE NUMBER OF POINTS TO OTHER PROCESSES
#if defined (HAVE_MPI)
      DO iproc = 0,max(ncsize,1)-1
        CALL mpi_scatter(cpld(cid)%POIN2RECV,1,mpi_integer,
     &  cpld(cid)%POIN2SEND(iproc+1),1,mpi_integer,iproc,comm,ierr)
        CALL p_sync()
      ENDDO

      !NUMBER OF PROCESSES FOR COMMUNICATION
      cpld(cid)%NP_SEND = count(cpld(cid)%POIN2SEND.GT.0)
      cpld(cid)%NP_RECV = count(cpld(cid)%POIN2RECV.GT.0)
      IF (cpld(cid)%POIN2SEND(ipid+1).GT.0) THEN
        cpld(cid)%NP_SEND_NOLOC = cpld(cid)%NP_SEND - 1
      ELSE
        cpld(cid)%NP_SEND_NOLOC = cpld(cid)%NP_SEND
      ENDIF
      IF (cpld(cid)%POIN2RECV(ipid+1).GT.0) THEN
        cpld(cid)%NP_RECV_NOLOC = cpld(cid)%NP_RECV - 1
      ELSE
        cpld(cid)%NP_RECV_NOLOC = cpld(cid)%NP_RECV
      ENDIF

#else
      cpld(cid)%POIN2SEND(1)  = cpld(cid)%POIN2RECV(1)
      cpld(cid)%NP_SEND = 1
      cpld(cid)%NP_RECV = 1
#endif
      IF (ierr.NE.0) THEN
        WRITE (lu,*) 'ERROR IN MPI_SCATTER POIN2SEND'
        CALL plante(1)
      ENDIF

!----------------------------------------------------------------------
!     COMMUNICATE SEND LIST TO SENDER  (DUMMY VALUES ARE ALSO COPIED)
      ireq = 0
      DO iproc = 0,max(ncsize,1)-1
        np =  cpld(cid)%POIN2RECV(iproc+1)
        IF (np.GT.0) THEN
          IF (iproc.NE.ipid) THEN
            ireq = ireq + 1
            itag = ipid+iproc*ncsize
            CALL  p_iwriti(loc_send_map(1:np,iproc+1),
     &           np,iproc,itag,cpld(cid)%REQLIST(ireq))
          ENDIF
        ENDIF
      ENDDO

      ! RECEIVE LIST AT SENDER
      ireq = 0
      DO iproc = 0,max(ncsize,1)-1
        np =  cpld(cid)%POIN2SEND(iproc+1)
        IF (np.GT.0) THEN
          IF (iproc.NE.ipid) THEN
            ireq = ireq + 1
            itag = iproc+ipid*ncsize
            CALL  p_ireadi(cpld(cid)%SEND_MAP(1:np,iproc+1),np,
     &           iproc,itag,cpld(cid)%REQLIST(ireq))
          ELSE
            ! LOCAL COPY
            DO ipoin = 1,np
              cpld(cid)%SEND_MAP(ipoin,iproc+1) =
     &            loc_send_map(ipoin,iproc+1)
            ENDDO
          ENDIF
        ENDIF
      ENDDO

      IF (ncsize.GT.1) THEN
        IF (ireq.GT.0) THEN
          CALL p_wait_paraco(cpld(cid)%REQLIST(1:ireq),ireq)
        ENDIF
      ENDIF

      ! THIS P_SYNC SOLVES A SEG_FAULT; NO IDEA WHY
      IF (ncsize.GT.1) THEN
        CALL p_sync
      ENDIF

!----------------------------------------------------------------------

      !MAXIMUM NUMBER OF POINTS SEND TO THIS PROCESS

      nps  = cpld(cid)%NP_SEND
      npr  = cpld(cid)%NP_RECV

      !ALLOCATE DATA FOR SENDING AND RECEIVING
      IF (nps.GT.0) THEN
        maxp = maxval(cpld(cid)%POIN2SEND(1:max(ncsize,1)))
        ALLOCATE(cpld(cid)%SEND_BUF(maxvar*maxp,nps),stat=ierr)
        cpld(cid)%SEND_BUF = -99
        CALL check_allocate(ierr,'SEND_DATA')
      ENDIF

      IF (npr.GT.0) THEN
        maxp = maxval(cpld(cid)%POIN2RECV(1:max(ncsize,1)))
        ALLOCATE(cpld(cid)%RECV_BUF(maxvar*maxp,npr),stat=ierr)
        cpld(cid)%RECV_BUF = -999
        CALL check_allocate(ierr,'RECV_BUF')
      ENDIF

      !DEALLOCATE DATA THAT IS NOT NEEDED ANYMORE

      DEALLOCATE (cpld(cid)%SEND_NPOIN,stat=ierr)
      CALL check_allocate(ierr,'DEAL:SEND_NPOIN')

      DEALLOCATE (cpld(cid)%SEND_KNOLG,stat=ierr)
      CALL check_allocate(ierr,'DEAL:SEND_KNOLG')

      DEALLOCATE(tmp_knolg,stat=ierr)
      CALL check_allocate(ierr,'DEAL:TMP_KNOLG')

      DEALLOCATE(send_node,stat=ierr)
      CALL check_allocate(ierr,'DEAL:SEND_NODE')

      DEALLOCATE(recv_list,stat=ierr)
      CALL check_allocate(ierr,'DEAL:RECV_LIST')

      DEALLOCATE(nsend,stat=ierr)
      CALL check_allocate(ierr,'DEAL:NSEND')

      DEALLOCATE(do_send,stat=ierr)
      CALL check_allocate(ierr,'DEAL:DO_SEND')

      DEALLOCATE(tmp_send,stat=ierr)
      CALL check_allocate(ierr,'DEAL:TMP_SEND')

      DEALLOCATE(loc_send_map,stat=ierr)
      CALL check_allocate(ierr,'DEAL:LOC_SEND_MAP')

      !----------------------------------
      ! CHECK THAT ALL POINTS RECEIVE INFORMATION

      ALLOCATE(tmp_log(npoin),stat=ierr)
      CALL check_allocate(ierr,'AL:TMPLOG')
      tmp_log = .false.
      DO ip_p=1,npp
        DO   iproc =1,max(ncsize,1)
          np  = cpld(cid)%POIN2RECV(iproc)
          IF (np.GT.0) THEN
          ! EXTRACT DATA FROM SEND STRUCTURE
            DO ipoin = 1,npoin
              itmp = cpld(cid)%RECV_MAP(iproc,ipoin,ip_p)
              IF (itmp.GT.0) THEN
                tmp_log(ipoin) = .true.
              ENDIF
            ENDDO
          ENDIF
        ENDDO
      ENDDO

      IF (all(tmp_log)) THEN
        WRITE (lu,*) 'ALL NODES RECEIVE INFORMATION'
      ELSE
        WRITE (lu,*) 'WARNING: NO INFORMATION FOR ',
     &      count(.NOT.tmp_log), ' NODES OUT OF ', npoin
      ENDIF
      DEALLOCATE(tmp_log,stat=ierr)
      CALL check_allocate(ierr,'DEAL:TMPLOG')

!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE


!                   **********************
                    SUBROUTINE send_couple
!                   **********************
     &(cid,npoin,nvar,varcouple)
!
!BRIEF    .
!
!
! NOTES: IT IS ASSUMED THAT THE RECEIVE CALL FOLLOWS THE SEND CALL
!         AND THAT TWO CONSECUTIVE SENDS DO NOT OCCUR. HOWEVER, THIS IS
!         NOT CHECKED
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER       , INTENT(IN) :: NPOIN, NVAR,CID
      TYPE(bief_obj), INTENT(IN) :: VARCOUPLE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      INTEGER :: IPOIN, IVAR, I, ITMP, NP, IPROC,IP,IPB, DUMMY
      INTEGER :: IREQ, ITAG , NPV, NOERR, IERR

!
!-----------------------------------------------------------------------

!     CHECK INPUT
      IF (nvar.GT.maxvar) THEN
        WRITE (lu,*) 'TOO MANY VARIABLES'
        CALL plante(1)
      ENDIF
      IF (nvar.EQ.0) THEN
        RETURN
      ENDIF

!     ALLOCATE REQLIST, AS IT IS AUTOMATICALLY DEALLOCATED BY WAIT
      IF (ncsize.GT.1) THEN
          ALLOCATE(cpld(cid)%REQLIST(max(ncsize,1)-1),stat=ierr)
          CALL check_allocate(ierr,'REQLIST')
      ENDIF
!     ALREADY RECEIVE THE DATA TO OPTIMIZE MPI COMMUNICATION
!     THE REAL RECEIVE WILL BE FINISHED BY THE MPI_WAIT CALL
      ip  = 0
      ipb = 0
      DO   iproc =1,max(ncsize,1)
        np  = cpld(cid)%POIN2RECV(iproc)
        npv = np*nvar
        IF (np.GT.0) THEN
          ipb  = ipb + 1
          IF (ipid.NE.iproc-1) THEN
            !RECEIVE DATA   (NON-BLOCKING)
            ip  = ip + 1
            itag = (iproc-1)+ipid*ncsize
            CALL  p_iread(cpld(cid)%RECV_BUF(1:npv,ipb),npv,8,iproc-1,
     &               itag,cpld(cid)%REQLIST(ip))
          ENDIF
        ENDIF
      ENDDO

      IF (ncsize.GT.1) THEN
        CALL p_sync
      ENDIF

      ireq = 0
!     ADD DATA TO SEND STRUCTURE
      ip = 0
      noerr = 0
      DO iproc=1,max(ncsize,1)
        i = 1
        np = cpld(cid)%POIN2SEND(iproc)
        IF (np.GT.0) THEN
          ip = ip +1
          DO ivar =1,nvar
            DO ipoin=1,np
              itmp = cpld(cid)%SEND_MAP(ipoin,iproc)
              IF (itmp.GT.npoin.OR. itmp.LT.1) THEN
                WRITE (lu,*) 'WRONG VALUE FOR SENDMAP: ',
     &             itmp, ' FOR ',ipoin,ivar,ip,i
                noerr = noerr + 1
                cycle
              ENDIF
              cpld(cid)%SEND_BUF(i,ip) = varcouple%ADR(ivar)%P%R(itmp)
              i = i+1
            ENDDO
          ENDDO
        ENDIF
      ENDDO
      IF (ncsize.GT.1) THEN
        CALL p_sync
      ENDIF

!      SEND DATA (NON BLOCKING)
      ip = 0
      DO iproc=1,max(ncsize,1)
        np = cpld(cid)%POIN2SEND(iproc)*nvar
        IF (np.GT.0) THEN
          ip = ip + 1
          IF (iproc-1.NE.ipid) THEN
            itag = (iproc-1)*ncsize+ipid
            ireq = ireq+1
            CALL  p_iwrit(cpld(cid)%SEND_BUF(1:np,ip),
     &           np,8,iproc-1,itag,dummy)
          ELSE
            ! COPY DATA ON THE SAME PROCESSOR
            itmp = count(cpld(cid)%POIN2RECV(1:iproc).GT.0)
            cpld(cid)%RECV_BUF(1:np,itmp) =
     &             cpld(cid)%SEND_BUF(1:np,ip)
          ENDIF
        ENDIF
      ENDDO

      ! SHUT DOWN IF THERE IS AN ERROR
      IF (ncsize.GT.1) THEN
        CALL p_sync
      ENDIF
      IF (noerr.GT.0) THEN
        CALL plante(1)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE
!
!-----------------------------------------------------------------------
!
!                   *************************************
                            SUBROUTINE receive_couple
!                   *************************************
     &               (cid,npoin,nvar,varcouple,default_val)
!
!BRIEF    .
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         ,INTENT(IN   ) :: NPOIN, CID, NVAR
      DOUBLE PRECISION, INTENT(IN), DIMENSION(NVAR), OPTIONAL ::
     &   default_val
      type(bief_obj), INTENT(INOUT) :: varcouple

!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      INTEGER :: IPROC, ITMP , NP, I, IP, IVAR , IPR, NPR
      INTEGER :: NOERR, IPP,NPP
      DOUBLE PRECISION :: WP
      DOUBLE PRECISION, DIMENSION(NVAR) ::  DEF_VAR
!-----------------------------------------------------------------------
!
! TODO: CHECK THAT NUMBER OF VARIABLES SEND AND RECEIVED IS EQUAL

      CALL debug_info('RECVMAP-START')

!     CHECK INPUT
      IF (nvar.GT.maxvar) THEN
        WRITE (lu,*) 'TOO MANY VARIABLES'
        CALL plante(1)
      ENDIF
      IF (nvar.EQ.0) THEN
        RETURN
      ENDIF
 
!     SET DEFAULT VALUES (IN CASE NO VALUE IS SPECIFIED)
      IF (PRESENT(default_val)) THEN
        def_var =  default_val
      ELSE
        def_var =0.0
      ENDIF

      DO ip = 1,npoin
        DO ivar=1,nvar
          varcouple%ADR(ivar)%P%R(ip) = def_var(ivar)
        ENDDO
      ENDDO

      npp = cpld(cid)%NPP

      ! WAIT UNTIL ALL DATA ARE RECEIVED
      npr  = cpld(cid)%NP_RECV_NOLOC
      IF (ncsize.GT.1) THEN
        IF (npr.GT.0) THEN
          CALL p_wait_paraco(cpld(cid)%REQLIST(1:npr),npr)
        ENDIF
      ENDIF

      !UNPACK DATA
      noerr = 0

      DO ipp = 1,npp
        ipr =  0
        DO   iproc =1,max(ncsize,1)
          np = cpld(cid)%POIN2RECV(iproc)
          IF (np.GT.0) THEN
            ipr = ipr + 1
          ! EXTRACT DATA FROM SEND STRUCTURE
            DO ip = 1,npoin
              itmp = cpld(cid)%RECV_MAP(iproc,ip,ipp)
              IF (itmp.GT.0) THEN
                ! ERROR CHECKING
                IF (itmp.GT.np) THEN
                  noerr = noerr + 1
                  cycle
                ENDIF
                wp =  cpld(cid)%WEIGHTS(ip,ipp)
                IF (ipp.EQ.1) THEN
                  DO ivar=1,nvar
                    i = itmp+(ivar-1)*np
                    varcouple%ADR(ivar)%P%R(ip) =
     &                 wp*cpld(cid)%RECV_BUF(i,ipr)
                  ENDDO
                ELSE
                  DO ivar=1,nvar
                    i = itmp+(ivar-1)*np
                    varcouple%ADR(ivar)%P%R(ip) =
     &                 varcouple%ADR(ivar)%P%R(ip) +
     &                 wp*cpld(cid)%RECV_BUF(i,ipr)
                  ENDDO
                ENDIF
              ENDIF
            ENDDO
          ENDIF
        ENDDO
      ENDDO

      ! SHUT DOWN IF THERE IS AN ERROR
      IF (ncsize.GT.1) THEN
          CALL p_sync
      ENDIF
      IF (noerr.GT.0) THEN
        CALL plante(1)
      ENDIF

      CALL debug_info('END UNPACK')
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE

!                   *********************
                    SUBROUTINE end_couple
!                   *********************
!
!BRIEF  DEALLOCATE VARIABLE NEEDED FOR COUPLING
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      INTEGER :: ICPL, IERR
!
!-----------------------------------------------------------------------
!
      !DEACCLOCATE DATA FROM COUPLING STRUCTURES

      DO  icpl =1,nrcouple
        !MAPPING
        DEALLOCATE(cpld(icpl)%POIN2RECV,stat=ierr)
        DEALLOCATE(cpld(icpl)%POIN2SEND,stat=ierr)
        DEALLOCATE(cpld(icpl)%SEND_MAP,stat=ierr)
        DEALLOCATE(cpld(icpl)%RECV_MAP,stat=ierr)

        !DATA BUFFERS
        !TODO: THESE TWO POINTERS MAY NOT BE ALLOCATED
        DEALLOCATE(cpld(icpl)%RECV_BUF,stat=ierr)
        DEALLOCATE(cpld(icpl)%SEND_BUF,stat=ierr)
        DEALLOCATE(cpld(icpl)%REQLIST,stat=ierr)
        DEALLOCATE(cpld(icpl)%WEIGHTS,stat=ierr)
      ENDDO

      DO icpl =1,maxcouple
        DEALLOCATE(cpld(icpl)%SEND_NPOIN,stat=ierr)
      ENDDO
      !
      DEALLOCATE(cpld,stat=ierr)

!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE

!-----------------------------------------------------------------------
      SUBROUTINE debug_info (DEBUG_STRING,DEBUG_INT)

!BRIEF  PRINT INFORMATION ABOUT A VARIABLE.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      IMPLICIT NONE

      CHARACTER (LEN=*), INTENT(IN) :: DEBUG_STRING
      INTEGER, OPTIONAL, INTENT(IN) :: DEBUG_INT
      WRITE (lu,*) ''
      WRITE (lu,*) '***********************************************'
      IF (PRESENT(debug_int)) THEN
        WRITE (lu,*) debug_string, debug_int
      ELSE
        WRITE (lu,*) debug_string
      ENDIF
      WRITE (lu,*) '***********************************************'
      WRITE (lu,*)''
      END SUBROUTINE

      END MODULE couple_mod