!!  Copyright (C) Stichting Deltares, 2005-2024.
!!
!!  This file is part of iMOD.
!!
!!  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 3 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: imod.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands.

module pks_iarmwp
implicit none

logical :: liarmwp = .false. 

logical :: lfirst = .true.

type xpType
   integer :: xprnk = -1
   integer :: xpbuf = 0
   integer :: nid   = 0
   integer, pointer, dimension(:)   :: idx => null()
   integer, pointer, dimension(:)   :: id  => null()
   integer, pointer, dimension(:,:) :: lrc => null()    
end type xpType

! serial
type xpaType
   integer                             :: nrxp = 0 
   integer, pointer, dimension(:)      :: xpflg => null() 
   type(xpType), pointer, dimension(:) :: xp => null()
end type xpaType
type(xpaType), pointer, dimension(:) :: xpa => null()

! parallel
integer :: nrxp = 0
type(xpType), pointer, dimension(:) :: xp => null()

! buffer
integer :: lenb
real, dimension(:), allocatable :: rb, sb

save

end module pks_iarmwp
    
module pks_imod_utl
use imod_idf
implicit none

contains
 
!###====================================================================
 subroutine pks_imod_utl_idfmerge_init()
!###====================================================================
 ! locals
 logical :: pks7mpimasterwrite, lpks
 !......................................................................
 call pks7mpiactive(lpks)
 if (lpks) then
    if (pks7mpimasterwrite()) midfout = 1 ! .true.
 end if   
 
! write(*,*) 'lidfout ini =',midfout
 
 end subroutine pks_imod_utl_idfmerge_init
 
!###====================================================================
 subroutine pks_imod_utl_write_idfmergefile(wd,fname)
!###====================================================================
 ! modules
 use pksmpi_mod
 ! arguments
 character(len=*), intent(in) :: wd
 character(len=*), intent(out) :: fname
 ! locals
 logical :: lpks
 integer :: i, j, iu
 character(len=1024) :: s
 character(len=1024), dimension(10) :: sa
 logical :: pks7mpimasterwrite
 !......................................................................
 
 call pks7mpiactive(lpks)
 if (.not.lpks) return
 !## only master process write the pksidfout.txt file
 if(.not.pks7mpimasterwrite())return

 write(fname,'(2a)') trim(wd), '\pksidfout.txt' 
 call osd_s_filename(fname)
 
 iu = 0
 call imod_utl_openasc(iu,fname,'w')

 write(s,*) nrproc
 write(iu,'(a)') trim(adjustl(s))
 
 do i = 1, nidfout
    s = idfout(i) 
    call imod_utl_swapslash(s)
    write(iu,'(a)') trim(s)
 end do
 close(iu)
 
 if(associated(idfout)) deallocate(idfout)
 
 return
 end subroutine pks_imod_utl_write_idfmergefile 
  
!###====================================================================
 subroutine pks_imod_utl_idfmerge(fname,ipksdelete)
!###====================================================================
 ! arguments
 character(len=*), intent(in) :: fname
 integer,intent(in) :: ipksdelete
 ! functions
 logical :: pks7mpimasterwrite
 ! locals
 integer :: iu
 logical :: lpks, lex
 type(idfobj) :: gidf, lidf
 integer :: ios, gncol, gnrow, nrproc, iproc, ic1, ic2, ir1, ir2, ip0, i, j, n, m, iirow, iicol, icol, irow
! real :: xmin, xmax, ymin, ymax, dx, dy
 character(len=1024) :: s, fp0, f
 real :: nodata
 character(len=1) :: slash 
 REAL(KIND=8),POINTER,DIMENSION(:) :: XCOL=>null(),YROW=>null(),XCOL_BU=>null(),YROW_BU=>null()
 REAL(KIND=8) :: DX1,DX2,DY1,DY2,XC,YC
 
 !......................................................................

 call pks7mpibarrier()
 if (pks7mpimasterwrite()) then
    call imod_utl_printtext('Start merging PKS output IDF files',0)
    
    inquire(file=fname,exist=lex)
    if(.not.lex) call imod_utl_printtext('File '//trim(fname)//' does not exist',2)
    
    iu = getunit(); open(unit=iu,file=fname,status='old')
    read(iu,*,iostat=ios) nrproc
    if(ios.ne.0) call imod_utl_printtext('Error reading '//trim(fname),2)

    CALL IDFNULLIFY(GIDF); CALL IDFNULLIFY(LIDF)
        
    !## determine common network
    do
    
       read(iu,'(a1024)',iostat=ios) s; if (ios.ne.0) exit
       if(len_trim(s).eq.0) cycle; if(s(1:1).eq.'#') cycle
       if(index(s,'"').gt.0)then
        read(s,*) fp0  != s
       else
        fp0=s
       endif
       
!       !## skip all others than heads
!       if(index(fp0,'HEAD_').le.0)cycle
       
       !## start point to insert process-number
       ip0=len_trim(fp0)-9
       
       N=0; M=0; IIROW=0; IICOL=0
       
       do iproc = 1, nrproc
          f = fp0(1:ip0)
          write(s,'(a,i3.3)') '_p',iproc-1
          f=trim(f)//trim(s)//'.idf'
          ! check is file exist
          inquire(file=f,exist=lex)
          if (.not.lex)cycle

          ! read the IDF header info only
          if (.not.idfread(lidf,f,0)) call imod_utl_printtext('Error reading '//trim(f),2); close(lidf%iu)

          N=N+LIDF%NCOL+1
          IF(ASSOCIATED(XCOL))THEN
           IF(N.GT.SIZE(XCOL))THEN
            ALLOCATE(XCOL_BU(N)); DO ICOL=1,SIZE(XCOL); XCOL_BU(ICOL)=XCOL(ICOL); ENDDO; DEALLOCATE(XCOL); XCOL=>XCOL_BU
           ENDIF
          ELSE
           IF(.NOT.ASSOCIATED(XCOL))ALLOCATE(XCOL(N))
          ENDIF
          M=M+LIDF%NROW+1
          IF(ASSOCIATED(YROW))THEN
           IF(M.GT.SIZE(YROW))THEN
            ALLOCATE(YROW_BU(M)); DO ICOL=1,SIZE(YROW); YROW_BU(ICOL)=YROW(ICOL); ENDDO; DEALLOCATE(YROW); YROW=>YROW_BU
           ENDIF
          ELSE
           IF(.NOT.ASSOCIATED(YROW))ALLOCATE(YROW(M))
          ENDIF
          
          IF(.NOT.IDFFILLSXSY(LIDF))call imod_utl_printtext('Error allocating sx/sy vectors '//trim(f),2)
          !## store x and y's
          DO IROW=0,LIDF%NROW
           IIROW=IIROW+1; YROW(IIROW)=LIDF%SY(IROW)
          ENDDO
          DO ICOL=0,LIDF%NCOL
           IICOL=IICOL+1; XCOL(IICOL)=LIDF%SX(ICOL)
          ENDDO
          CALL IDFDEALLOCATEX(LIDF)
          
       end do
    
       IF(.NOT.ASSOCIATED(XCOL))THEN
        STOP 'NOTHING TO MERGE, PROBABLY MERGE-FILE IS EMPTY !!!!'
       ENDIF
    
       !## get network - sort x and y-coordinates
       CALL IMOD_UTL_QKSORT(N,N,XCOL)
       CALL IMOD_UTL_QKSORT(M,M,YROW)
       !## remove duplicates
       IICOL=1; DO ICOL=2,SIZE(XCOL)
        IF(XCOL(ICOL).GT.XCOL(IICOL))THEN; IICOL=IICOL+1; XCOL(IICOL)=XCOL(ICOL); ENDIF
       ENDDO
       IIROW=1; DO IROW=2,SIZE(YROW)
        IF(YROW(IROW).GT.YROW(IIROW))THEN; IIROW=IIROW+1; YROW(IIROW)=YROW(IROW); ENDIF
       ENDDO
       LIDF%NCOL=IICOL-1; LIDF%NROW=IIROW-1
       IF(.NOT.IDFALLOCATESXY(LIDF))call imod_utl_printtext('Error allocating new network sx/sy vector',2)
       DO ICOL=1,IICOL; LIDF%SX(ICOL-1)=XCOL(ICOL); ENDDO; LIDF%XMIN=LIDF%SX(0); LIDF%XMAX=LIDF%SX(LIDF%NCOL)
       IR1=IIROW+1; DO IROW=1,IIROW
        IR1=IR1-1; LIDF%SY(IROW-1)=YROW(IR1)
       ENDDO; LIDF%YMAX=LIDF%SY(0); LIDF%YMIN=LIDF%SY(LIDF%NROW)
       !## determine whether equi- or non-equi
       DX1=10.0D10; DX2=0.0D0
       DO ICOL=1,LIDF%NCOL
        DX1=MIN(DX1,LIDF%SX(ICOL)-LIDF%SX(ICOL-1))
        DX2=MAX(DX2,LIDF%SX(ICOL)-LIDF%SX(ICOL-1))
       ENDDO
       DY1=10.0D10; DY2=0.0D0
       DO IROW=1,LIDF%NROW
        DY1=MIN(DY1,LIDF%SY(IROW-1)-LIDF%SY(IROW))
        DY2=MAX(DY2,LIDF%SY(IROW-1)-LIDF%SY(IROW))
       ENDDO
       !## equidistantial network
       IF(DX1.EQ.DX2.AND.DY1.EQ.DY2)THEN
        LIDF%IEQ=0; LIDF%DX=DX1; LIDF%DY=DY1
       !## non-equidistantial network
       ELSE
        LIDF%IEQ=1; LIDF%DX=0.0D0; LIDF%DY=0.0D0
       ENDIF
       !## copy settings to the global IDF-file
       CALL IDFDEALLOCATEX(GIDF); CALL IDFCOPY(LIDF,GIDF)
       IF(.NOT.IDFALLOCATEX(GIDF))call imod_utl_printtext('Error allocating new network x array',2)
       
       i=0; do iproc = 1, nrproc
          f = fp0(1:ip0)
          write(s,'(a,i3.3)') '_p',iproc-1
          f=trim(f)//trim(s)//'.idf'
       
          ! check is file exist
          inquire(file=f,exist=lex)
          if (.not.lex)cycle

          i=i+1
          
          !## read the IDF
          if (.not.idfread(lidf,f,1))call imod_utl_printtext('Error reading '//trim(f),2)   
          !## get indices of model to fit in
          CALL UTL_IDFGETLOC(LIDF,1        ,1        ,XC,YC); CALL UTL_IDFIROWICOL(GIDF,IR1,IC1,XC,YC)
          CALL UTL_IDFGETLOC(LIDF,LIDF%NROW,LIDF%NCOL,XC,YC); CALL UTL_IDFIROWICOL(GIDF,IR2,IC2,XC,YC)
          GIDF%NODATA =LIDF%NODATA
          GIDF%TOP=LIDF%TOP
          GIDF%BOT=LIDF%BOT
          !## set all to nodata for the first time
          IF(i.eq.1)GIDF%X=GIDF%NODATA

          GIDF%COMMENT=LIDF%COMMENT 

          GIDF%X(IC1:IC2,IR1:IR2) = LIDF%X
          IF(IPKSDELETE.EQ.1)THEN
           LIDF%IU=IMOD_UTL_GETUNIT(); OPEN(UNIT=LIDF%IU,FILE=F); CLOSE(LIDF%IU,STATUS='DELETE')
          ENDIF
          CALL IDFDEALLOCATEX(LIDF)
       END DO
       
       !## write the global IDF
       f = fp0(1:ip0)//'.idf' !; F=IMOD_UTL_CAPF(F,'U')
       call imod_utl_printtext('Writing '//trim(f)//'...',0)
       lidfout=2
       if (.not.idfwrite(gidf,f,0)) call imod_utl_printtext('Could not write '//trim(f),2)

       IF(ASSOCIATED(XCOL))DEALLOCATE(XCOL)
       IF(ASSOCIATED(YROW))DEALLOCATE(YROW)
       
    end do
    
    !## close file
    call idfdeallocatex(gidf)
    close(iu) !,status='delete')
        
    call imod_utl_printtext('Done merging PKS output IDF files',0)
 end if
 call pks7mpibarrier() 

 return
 end subroutine pks_imod_utl_idfmerge 
 
!###====================================================================
 subroutine pks_imod_utl_iarmwp_xch_init(iarmwp)
!###====================================================================
! modules
 use pks_iarmwp
! arguments
 integer, intent(in) :: iarmwp
! locals
 logical :: lpks
!......................................................................
 call pks7mpiactive(lpks)
 if (.not.lpks) return
 
 if (iarmwp.eq.1) then
    liarmwp = .true.
 end if
 
 return   
 end subroutine pks_imod_utl_iarmwp_xch_init
 
 !###====================================================================
 subroutine pks_imod_utl_iarmwp_xch_disable()
!###====================================================================
! modules
 use pks_iarmwp, only: liarmwp
!......................................................................
 
 liarmwp = .false.
 
 return   
 end subroutine pks_imod_utl_iarmwp_xch_disable
 
!###====================================================================
 subroutine pks_imod_utl_iarmwp_xch_store(msir,msic,mfil,mfir,mfic,dxcid,ncol,nrow,nlay,iact)
!###====================================================================
 ! modules
 use rf2mf_module, only: pks
 use pks_iarmwp
 use pksmpi_mod, only: nrproc  
 ! arguments
 integer, intent(in) :: msir, msic, mfil, mfir, mfic, ncol, nrow, nlay, iact
 integer, dimension(ncol,nrow,nlay), intent(inout) :: dxcid 
 ! locals
 logical :: lpks, pmsovl(nrproc), pmfovl(nrproc)
 integer :: iproc, jproc, pmsnovl, pmfnovl, i, n, mfid
 integer :: ic1, ic2, ir1, ir2 ! non-overlapping
 integer :: jc1, jc2, jr1, jr2 ! overlapping
 !......................................................................
 call pks7mpiactive(lpks)
 if (.not.lpks) return
 
 if (.not.associated(xpa)) then
     allocate(xpa(nrproc)) !A1
     do iproc = 1, nrproc
        allocate(xpa(iproc)%xpflg(nrproc)) !A2
        xpa(iproc)%xpflg = 0
        allocate(xpa(iproc)%xp(nrproc)) ! A3
     end do
 end if
 
 if (iact.eq.2.and.lfirst) then
    lfirst = .false.
     do iproc = 1, nrproc
        do jproc = 1, nrproc 
           n = xpa(iproc)%xp(jproc)%nid
           if(n.gt.0) then
              allocate(xpa(iproc)%xp(jproc)%id(n))
           end if
           xpa(iproc)%xp(jproc)%nid = 0
        end do
        xpa(iproc)%xpflg = 0
     end do
 end if    
 
 pmsovl = .false.; pmfovl = .false.
 pmsnovl = -1; pmfnovl = -1
 do iproc = 1, nrproc
    ic1 = pks%partminmax(iproc,1); ic2 = pks%partminmax(iproc,2) ! non overlapping 
    ir1 = pks%partminmax(iproc,3); ir2 = pks%partminmax(iproc,4) ! non overlapping 
    jc1 = max(1,ic1-1); jc2 = min(ncol,ic2+1)                    ! overlapping
    jr1 = max(1,ir1-1); jr2 = min(nrow,ir2+1)                    ! overlapping  
     
    if (msic.ge.jc1 .and. msic.le.jc2 .and. &                    ! overlapping
        msir.ge.jr1 .and. msir.le.jr2) pmsovl(iproc) = .true.

    if (mfic.ge.jc1 .and. mfic.le.jc2 .and. &                    ! overlapping
        mfir.ge.jr1 .and. mfir.le.jr2) pmfovl(iproc) = .true.
    
    if (msic.ge.ic1 .and. msic.le.ic2 .and. &                    ! non overlapping
        msir.ge.ir1 .and. msir.le.ir2) pmsnovl = iproc

    if (mfic.ge.ic1 .and. mfic.le.ic2 .and. &                    ! non overlapping
        mfir.ge.ir1 .and. mfir.le.ir2) pmfnovl = iproc
 end do
 if(pmfnovl.le.0.and.pmsnovl.gt.0) then
    write(*,*) 'Warning IARMWP=1, MODFLOW cell does not belong to any partition!',mfic,mfir
    return
 end if
 if(pmfnovl.gt.0.and.pmsnovl.le.0.) then
    write(*,*) 'Warning IARMWP=1, MetaSWAP svat does not belong to any partition!',msic,msir
 end if
 if(pmfnovl.le.0.or.pmsnovl.le.0.) then
    return
 end if

 mfid = abs(dxcid(mfic,mfir,mfil))
 do iproc = 1, nrproc
    jproc = -1 
    if (     pmsovl(iproc).and.     pmfovl(iproc)) cycle
    if (.not.pmsovl(iproc).and..not.pmfovl(iproc)) cycle
    ! svat not in overlapping parition & cell     in overlapping parition
    if (.not.pmsovl(iproc).and.     pmfovl(iproc)) then
       jproc = pmsnovl         
    end if    
    ! svat     in overlapping parition & cell not in overlapping parition
    if (     pmsovl(iproc).and..not.pmfovl(iproc)) then
       jproc = pmfnovl
    end if
    if (jproc.le.0) then
       write(*,*) 'Program error, pks_imod_utl_iarmwp_xch_store',pmsnovl,pmfnovl,msic,msir,mfic,mfir    
       stop 
    end if
    if(iproc.eq.jproc) cycle
    xpa(iproc)%xpflg(jproc) = 1
    xpa(iproc)%xp(jproc)%nid = xpa(iproc)%xp(jproc)%nid + 1 
    if (iact.eq.2) then
       i = xpa(iproc)%xp(jproc)%nid  
       xpa(iproc)%xp(jproc)%id(i) = mfid    
    end if
       
    xpa(jproc)%xpflg(iproc) = 1
    xpa(jproc)%xp(iproc)%nid = xpa(jproc)%xp(iproc)%nid + 1 
    if (iact.eq.2) then
       i = xpa(jproc)%xp(iproc)%nid  
       xpa(jproc)%xp(iproc)%id(i) = mfid    
    end if
    
!    if(iact.eq.2) then
!       dxcid(mfic,mfir,mfil) = -abs(dxcid(mfic,mfir,mfil))    
!    end if
    
 end do
 
 return   
 end subroutine pks_imod_utl_iarmwp_xch_store

!###====================================================================
 subroutine pks_imod_utl_iarmwp_xch_write(dxcid,ncol,nrow,nlay,ndxc,modwd)
!###====================================================================
 ! modules
 use pks_iarmwp
 use pksmpi_mod, only: nrproc
 use imod_utl, only: imod_utl_openasc
 ! arguments
 integer, intent(in) :: ncol, nrow, nlay, ndxc
 integer, dimension(ncol,nrow,nlay), intent(in) :: dxcid
 character(len=*), intent(in) :: modwd
 ! locals
 logical :: lpks
 integer :: iproc, jproc, id, n, i, j, lun, ilay, irow, icol
 character(len=1000) :: fname, s
 character(len=100), dimension(4) :: sa
 integer, dimension(:,:), allocatable :: iwrk
 ! functions
 integer :: cfn_unique_i
 !......................................................................
 call pks7mpiactive(lpks)
 if (.not.lpks) return

 allocate(iwrk(3,ndxc))
 do ilay = 1, nlay
    do irow = 1, nrow
       do icol = 1, ncol
          id = abs(dxcid(icol,irow,ilay))
          if (id.ne.0) then
             if (id.gt.ndxc) then
                write(*,*) 'Program error, stopping'
                stop
             end if 
             iwrk(1,id) = ilay 
             iwrk(2,id) = irow 
             iwrk(3,id) = icol 
          end if
       end do
    end do
 end do
 
 do iproc = 1, nrproc
    !write(*,'(a,i3.3,a)') '========= iproc = ',iproc,':' 
    do jproc = 1, nrproc
       if (xpa(iproc)%xpflg(jproc).eq.1) then
          xpa(iproc)%nrxp = xpa(iproc)%nrxp + 1    
       end if 
    end do
    !write(*,*) 'nxp = ',xp(iproc)%nrxp 
    do jproc = 1, nrproc
       if (xpa(iproc)%xpflg(jproc).eq.1) then
          n = cfn_unique_i(xpa(iproc)%xp(jproc)%id,xpa(iproc)%xp(jproc)%nid,0) 
          xpa(iproc)%xp(jproc)%nid = n 
          !write(*,*) jproc,':', xp(iproc)%xp(jproc)%nid
          !do i = 1, xp(iproc)%xp(jproc)%nid
          !   write(*,*) '-->',xp(iproc)%xp(jproc)%id(i)
          !end do
       end if 
    end do
    
    lun = 0
    write(fname,'(2a,i3.3)') trim(modwd),'pks_xch_iarmwp.p', iproc-1 
    call osd_s_filename(fname)
    write(*,'(a,1x,2a)') 'Writing',trim(fname),'...'
    call imod_utl_openasc(lun,fname,'w')
    write(s,*) xpa(iproc)%nrxp
    write(lun,'(a)') trim(adjustl(s))
    do jproc = 1, nrproc
       if (xpa(iproc)%xpflg(jproc).eq.1) then
          write(sa(1),*) jproc-1
          write(sa(2),*) xpa(iproc)%xp(jproc)%nid
          write(lun,'(1(a,1x),a)')(trim(adjustl(sa(j))),j=1,2)
          do i = 1, xpa(iproc)%xp(jproc)%nid
             id = xpa(iproc)%xp(jproc)%id(i)
             write(sa(1),*) id
             write(sa(2),*) iwrk(1,id)
             write(sa(3),*) iwrk(2,id)
             write(sa(4),*) iwrk(3,id)
             write(lun,'(3(a,1x),a)')(trim(adjustl(sa(j))),j=1,4)
          end do
       end if
    end do
    close(lun)
 end do

 ! cleanup
 do iproc = 1, nrproc
    if (associated(xpa(iproc)%xpflg)) deallocate(xpa(iproc)%xpflg)
    do jproc = 1, nrproc
       if (associated(xpa(iproc)%xp(jproc)%id)) &
          deallocate(xpa(iproc)%xp(jproc)%id)    
    end do
 end do
 deallocate(xpa)
 deallocate(iwrk)
 
 return   
 end subroutine pks_imod_utl_iarmwp_xch_write
  
!###====================================================================
 subroutine pks_imod_utl_iarmwp_xch_read(modwd)
!###====================================================================
 ! modules
 use pks_iarmwp
 use pksmpi_mod, only: myrank
 ! arguments
 character(len=*), intent(in) :: modwd
 ! locals
 character(len=1000) :: fname
 logical :: lpks
 integer :: lun, ixp, i, xprnk, nid, id, il, ir, ic, n
!......................................................................
 call pks7mpiactive(lpks)
 if (.not.lpks) return
 if (.not.liarmwp) return
 
 lun = 0
 write(fname,'(2a,i3.3)') trim(modwd),'\pks_xch_iarmwp.p', myrank 
 call osd_s_filename(fname)
 write(*,'(a,1x,2a)') 'Reading',trim(fname),'...'
 call imod_utl_openasc(lun,fname,'r')
 
 read(lun,*) nrxp
 allocate(xp(max(nrxp,1)))
 
 n = 0
 do ixp = 1, nrxp     
    read(lun,*) xprnk, nid  
    n = n + nid
    xp(ixp)%xprnk = xprnk; xp(ixp)%nid = nid
    allocate(xp(ixp)%idx(nid))
    allocate(xp(ixp)%id(nid))
    allocate(xp(ixp)%lrc(3,nid))
    do i = 1, nid 
       read(lun,*) id, il, ir, ic
       xp(ixp)%id(i)    = id
       xp(ixp)%lrc(1,i) = il 
       xp(ixp)%lrc(2,i) = ir 
       xp(ixp)%lrc(3,i) = ic 
    end do 
 end do
 close(lun)
 
 lenb = max(n,1)
 allocate(sb(lenb),rb(lenb))
 
 return   
 end subroutine pks_imod_utl_iarmwp_xch_read 
 
!###====================================================================
subroutine pks_imod_utl_iarmwp_xch(x,xcht)
!###====================================================================
! modules
 use pks_iarmwp
 use pksmpi_mod, only: myrank
 use pksmpiwrp_mod

! arguments
 real, dimension(*), intent(inout) :: x
 character(len=1), intent(in) :: xcht
! locals
 logical :: lpks
 integer :: ixp, bufptr, buflen, i, j, nrxp2
!......................................................................
 
 call pks7mpiactive(lpks)
 if (.not.lpks) return

 bufptr = 1
 nrxp2 = 0
 do ixp = 1, nrxp ! loop over exchange partners
    nrxp2 = nrxp2 + 1 
    xp(ixp)%xpbuf = bufptr
    buflen = xp(ixp)%nid
    call pks7mpiwrpirecvr( rb(bufptr),&
                           buflen,&
                           xp(ixp)%xprnk,&
                           0,&
                           pks7mpiwrpcomm_world,&
                           rreq(nrxp2) )
    bufptr = bufptr + buflen
 end do ! ixp

 ! check
 bufptr = bufptr - 1
 if (bufptr.gt.lenb) then
    write(*,*) 'Program error pks_imod_utl_iarmwp_xch: lenb!'
    call pksstop(' ')
 end if

 ! pack
 bufptr = 1
 do ixp = 1, nrxp
    do i = 1, xp(ixp)%nid
       j = xp(ixp)%idx(i) 
       if (j.eq.0) then
          write(*,*) 'Program error pks_imod_utl_iarmwp_xch: packing!',&
                      myrank, ixp, nrxp, j
          write(*,*) 'index ',j
          call pksstop(' ')
       end if
       sb(bufptr) = x(abs(j))
       if (xcht.eq.'q') then
          ! write(*,*) '@@@ PCK Q:', myrank, j, xp(ixp)%id(i), sb(bufptr)    
       end if    
      bufptr = bufptr + 1
    end do
 end do ! ixp
 
 ! check
 bufptr = bufptr - 1
 if (bufptr.gt.lenb) then
    write(*,*) 'Program error pks_imod_utl_iarmwp_xch: lenb!'
    call pksstop(' ')
 end if
      
 bufptr = 1

 nrxp2 = 0
 do ixp = 1, nrxp
    nrxp2 = nrxp2 + 1
    buflen = xp(ixp)%nid
    call pks7mpiwrpisendr( sb(bufptr),&
                           buflen,&
                           xp(ixp)%xprnk,&
                           0,&
                           pks7mpiwrpcomm_world,&
                           sreq(nrxp2) )
    bufptr = bufptr + buflen
 end do ! ixp

 ! check
 bufptr = bufptr - 1
 if (bufptr.gt.lenb) then
    write(*,*) 'Program error pks_imod_utl_iarmwp_xch: lenb!'
    call pksstop(' ')
 end if

 call pks7mpiwrpwaitall( nrxp2, rreq )

 ! unpacking
 do ixp = 1, nrxp
    bufptr = xp(ixp)%xpbuf
    do i = 1, xp(ixp)%nid
       j = xp(ixp)%idx(i)
       if (j.eq.0) then
          write(*,*) 'Program error pks_imod_utl_iarmwp_xch: unpacking!',&
                      myrank, ixp, nrxp, j
          call pksstop(' ')
       end if
       if (xcht.eq.'q') then
          x(abs(j)) = x(abs(j)) + rb(bufptr)
          ! write(*,*) '@@@ UPCK Q:', myrank, j, xp(ixp)%id(i), rb(bufptr), x(abs(j))
       end if   
       if (xcht.eq.'h') then
          if (j.gt.0) then ! I do not have this cell
             x(j) = rb(bufptr)   
             ! write(*,*) '@@@ UPCK H:', myrank, j, xp(ixp)%id(i), rb(bufptr)
          end if  
       end if      
       bufptr = bufptr + 1
    end do
 end do ! ixp
      
 ! check
 bufptr = bufptr - 1
 if (bufptr.gt.lenb) then
    write(*,*) 'Program error pks_imod_utl_iarmwp_xch: lenb!'
    call pksstop(' ')
 end if

 call pks7mpiwrpwaitall( nrxp2, sreq ) 
 
 return
 end subroutine pks_imod_utl_iarmwp_xch
 
 end module pks_imod_utl