!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!-----------------------------------------------------------------------
! CVS m_Router.F90,v 1.49 2008-05-12 02:09:02 jacob Exp
! CVS MCT_2_6_0
!BOP -------------------------------------------------------------------
!
! !MODULE: m_Router -- Router class 
!
! !DESCRIPTION:
! The Router data type contains all the information needed
! to send an AttrVect between a component on M MPI-processes and a component
! on N MPI-processes.   This module defines the Router datatype and provides
! methods to create and destroy one.
!
! !INTERFACE:

 module m_Router

      use m_realkinds, only : FP

      implicit none

      private   ! except

! !declare a private pointer structure for the real data
      type :: rptr
#ifdef SEQUENCE
        sequence
#endif
        real(FP),dimension(:),pointer :: pr
      end type

! !declare a private pointer structure for the integer data
      type :: iptr
#ifdef SEQUENCE
        sequence
#endif
        integer,dimension(:),pointer :: pi
      end type

! !PUBLIC TYPES:
      public :: Router	        ! The class data structure

      public :: rptr,iptr       ! pointer types used in Router
!\end{verbatim}
!% On return, pe_list is the processor ranks of the other
!% component to receive from/send to.  num_segs is the
!% number of segments out of my local AttrVect which must
!% be sent/received.  (In general, these wont coincide exactly
!% with the segments used to define the GlobalMap)
!% seg_start is the start *in the local AttrVect* of each segment
!%  (start goes from 1 to lsize(GSMap))
!% and seg_lengths is the length.
!\begin{verbatim}

    type Router
#ifdef SEQUENCE
      sequence
#endif
      integer :: comp1id                           ! myid
      integer :: comp2id                           ! id of second component
      integer :: nprocs	                           ! number of procs to talk to
      integer :: maxsize                           ! maximum amount of data going to a processor
      integer :: lAvsize                           ! The local size of AttrVect which can be 
                                                   ! used with this Router in MCT_Send/MCT_Recv
      integer :: numiatt                           ! Number of integer attributes currently in use
      integer :: numratt                           ! Number of real attributes currently in use
      integer,dimension(:),pointer   :: pe_list    ! processor ranks of send/receive in MCT_comm
      integer,dimension(:),pointer   :: num_segs   ! number of segments to send/receive
      integer,dimension(:),pointer   :: locsize    ! total of seg_lengths for a proc
      integer,dimension(:),pointer   :: permarr    ! possible permutation array
      integer,dimension(:,:),pointer :: seg_starts ! starting index
      integer,dimension(:,:),pointer :: seg_lengths! total length
      type(rptr),dimension(:),pointer :: rp1       ! buffer to hold real data
      type(iptr),dimension(:),pointer :: ip1       ! buffer to hold integer data
      integer,dimension(:),pointer   :: ireqs,rreqs  ! buffer for MPI_Requests
      integer,dimension(:,:),pointer :: istatus,rstatus  ! buffer for MPI_Status
    end type Router

! !PUBLIC MEMBER FUNCTIONS:
      public :: init            ! Create a Router
      public :: clean           ! Destroy a Router
      public :: print           ! Print info about a Router


    interface init  ; module procedure  &
        initd_, &       ! initialize a Router between two seperate components
        initp_ 	        ! initialize a Router locally with two GSMaps
    end interface
    interface clean ; module procedure clean_ ; end interface
    interface print ; module procedure print_ ; end interface

! !REVISION HISTORY:
! 15Jan01 - R. Jacob <jacob@mcs.anl.gov> - initial prototype
! 08Feb01 - R. Jacob <jacob@mcs.anl.gov> add locsize and maxsize 
!           to Router type
! 25Sep02 - R. Jacob <jacob@mcs.anl.gov> Remove type string.  Add lAvsize
! 23Jul03 - R. Jacob <jacob@mcs.anl.gov> Add status and reqs arrays used
!           in send/recv to the Router datatype.
! 24Jul03 - R. Jacob <jacob@mcs.anl.gov> Add real and integer buffers
!           for send/recv to the Router datatype.
! 22Jan08 - R. Jacob <jacob@mcs.anl.gov> Add ability to handle an unordered
!           GSMap by creating a new, ordered one and building Router from
!           that.  Save permutation info in Router datatype.
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname='MCT::m_Router'

 contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: initd_ - initialize a Router between two seperate components
!
! !DESCRIPTION:
! The routine {\tt initd\_()} exchanges the {\tt GSMap} with the
! component identified by {\tt othercomp} and then calls {\tt initp\_()}
! to build a Router {\tt Rout} between them.
!
! {\bf N.B.} The  {\tt GSMap} argument must be declared so that the index values
! on a processor are in ascending order.
!
! !INTERFACE:

 subroutine initd_(othercomp,GSMap,mycomm,Rout )
!
! !USES:
!
      use m_GlobalSegMap, only :GlobalSegMap
      use m_ExchangeMaps,only: MCT_ExGSMap => ExchangeMap
      use m_mpif90
      use m_die

      implicit none

! !INPUT PARAMETERS:
!
      integer, intent(in)	       :: othercomp
      integer, intent(in)	       :: mycomm
      type(GlobalSegMap),intent(in)    :: GSMap     ! of the calling comp

! !OUTPUT PARAMETERS:
!
      type(Router), intent(out)        :: Rout

! !REVISION HISTORY:
! 15Jan01 - R. Jacob <jacob@mcs.anl.gov> - initial prototype
! 06Feb01 - R. Jacob <jacob@mcs.anl.gov> - Finish initialization
!           of the Router.  Router now works both ways.
! 25Apr01 - R. Jacob <jacob@mcs.anl.gov> - Eliminate early 
!           custom code to exchange GSMap components and instead
!           the more general purpose routine in m_ExchangeMaps.
!           Use new subroutine OrderedPoints in m_GlobalSegMap
!           to construct the vector of local and remote GSMaps.
!           Clean-up code a little.
! 03May01 - R. Jacob <jacob@mcs.anl.gov> - rename to initd and
!           move most of code to new initp routine
!
!EOP ___________________________________________________________________
!
  character(len=*),parameter :: myname_=myname//'::initd_'

  type(GlobalSegMap)    :: RGSMap  !  the other GSMap
  integer ::		   ier

!--------------------------begin code-----------------------

!!!!!!!!!!!!!!!!!Exchange of global map data 

  call MCT_ExGSMap(GSMap,mycomm,RGSMap,othercomp,ier)
  if(ier /= 0) call die(myname_,'ExGSMap',ier)

!!!!!!!!!!!!!!!!!Begin comparison of globalsegmaps

  call initp_(GSMap,RGSMap, mycomm, Rout)

 end subroutine initd_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: initp_ - initialize a Router from two GlobalSegMaps
!
! !DESCRIPTION:
!
! Given two GlobalSegmentMaps {\tt GSMap} and {\tt RGSMap}, intialize a
! Router {\tt Rout} between them.  Use local communicator {\tt mycomm}.
!
! {\bf N.B.} The two {\tt GSMap} arguments must be declared so that the index values
! on a processor are in ascending order.
!
! !INTERFACE:

 subroutine initp_(inGSMap,inRGSMap,mycomm,Rout )
!
! !USES:
!
      use m_GlobalSegMap,  only :GlobalSegMap
      use m_GlobalSegMap,  only :ProcessStorage
      use m_GlobalSegMap,  only :GSMap_comp_id => comp_id
      use m_GlobalSegMap,  only :GSMap_increasing => increasing
      use m_GlobalSegMap,  only :GlobalSegMap_copy => copy
      use m_GlobalSegMap,  only :GlobalSegMap_init => init
      use m_GlobalSegMap,  only :GlobalSegMap_clean => clean
      use m_GlobalSegMap,  only :GlobalSegMap_OPoints => OrderedPoints
      use m_GlobalSegMap,  only :GlobalSegMap_ngseg => ngseg  ! rml
      use m_GlobalSegMap,  only :GlobalSegMap_nlseg => nlseg  ! rml
      use m_GlobalSegMap,  only :GlobalSegMap_max_nlseg => max_nlseg  ! rml

      use m_GlobalToLocal, only :GlobalToLocalIndex
      use m_MCTWorld,      only :MCTWorld
      use m_MCTWorld,      only :ThisMCTWorld

      use m_Permuter      ,only:Permute
      use m_MergeSorts    ,only:IndexSet
      use m_MergeSorts    ,only:IndexSort

      use m_mpif90
      use m_die

!      use m_zeit


      use m_stdio    ! rml
!      use shr_timer_mod        ! rml timers

      implicit none

! !INPUT PARAMETERS:
!
      type(GlobalSegMap), intent(in)	:: inGSMap
      type(GlobalSegMap), intent(in)	:: inRGSMap
      integer	     ,    intent(in)	:: mycomm

! !OUTPUT PARAMETERS:
!
      type(Router),      intent(out)	:: Rout

! !REVISION HISTORY:
! 03May01 - R.L. Jacob <jacob@mcs.anl.gov> - Initial code brought
!           in from old init routine.
! 31Jul01 - Jace A Mogill <mogill@cray.com>
!           Rewrote to reduce number of loops and temp storage
! 26Apr06 - R. Loy <rloy@mcs.anl.gov> - recode the search through
!           the remote GSMap to improve efficiency
! 05Jan07 - R. Loy <rloy@mcs.anl.gov> - improved bound on size of 
!           tmpsegcount and tmpsegstart
! 15May07 - R. Loy <rloy@mcs.anl.gov> - improved bound on size of
!           rgs_lb and rgs_ub
! 25Jan08 - R. Jacob <jacob@mcs.anl.gov> - Dont die if GSMap is not
!           increasing.  Instead, permute it to increasing and proceed.
!EOP -------------------------------------------------------------------

  character(len=*),parameter :: myname_=myname//'::initp_'
  integer			     :: ier,i,j,k,m,n
  integer			     :: mysize,myPid,othercomp
  integer			     :: lmaxsize,totallength
  integer                            :: maxsegcount,count
  logical, dimension(:), allocatable :: tmppe_list
  integer, dimension(:,:), pointer   :: tmpsegcount,tmpsegstart


  integer :: my_left        ! Left point in local segment (global memory)
  integer :: my_right       ! Right point in local segment (global memory)
  integer :: r_left         ! Left point in remote segment (global memory)
  integer :: r_right        ! Right point in remote segment (global memory)
  integer :: nsegs_overlap  ! Number of segments that overlap between two procs


  integer :: ngseg, nlseg
  integer :: myseg, rseg
  integer :: prev_right         ! Rightmost local point in previous overlapped segment
  integer :: local_left, local_right
  integer,allocatable  :: mygs_lb(:),mygs_ub(:),mygs_len(:),mygs_lstart(:)
  integer :: r_ngseg
  integer :: r_max_nlseg   ! max number of local segments in RGSMap
  integer,allocatable  :: rgs_count(:),rgs_lb(:,:),rgs_ub(:,:)
  integer,allocatable  :: nsegs_overlap_arr(:)

  integer :: overlap_left, overlap_right, overlap_diff

  integer :: proc, nprocs

  integer :: max_rgs_count, max_overlap_segs
  type(GlobalSegMap)	:: GSMap
  type(GlobalSegMap)	:: RGSMap
  integer, dimension(:), pointer   :: gpoints
  integer, dimension(:), pointer   :: permarr
  integer, dimension(:), pointer   :: rpermarr
  integer  :: gmapsize


  integer,save  :: t_initialized=0   ! rml timers
  integer,save  :: t_loop            ! rml timers
  integer,save  :: t_loop2            ! rml timers
  integer,save  :: t_load            ! rml timers

  call MP_comm_rank(mycomm,myPid,ier)
  if(ier/=0) call MP_perr_die(myname_,'MP_comm_rank',ier)

  nullify(Rout%permarr)

  if (.not. GSMap_increasing(inGSMap)) then
    if(myPid == 0) call warn(myname_,'GSMap indices not increasing...Will correct')
    call GlobalSegMap_OPoints(inGSMap,myPid,gpoints)
    gmapsize=ProcessStorage(inGSMap,myPid)
    allocate(permarr(gmapsize), stat=ier)
    if(ier/=0) call die(myname_,'allocate permarr',ier)
    call IndexSet(permarr)
    call IndexSort(permarr,gpoints)
    call Permute(gpoints,permarr,gmapsize)
    call GlobalSegMap_init(GSMap,gpoints,mycomm,inGSMap%comp_id,gsize=inGSMap%gsize)

    allocate(Rout%permarr(gmapsize),stat=ier)
    if(ier/=0) call die(myname_,'allocate Router%permarr',ier)
    Rout%permarr(:)=permarr(:)

    deallocate(gpoints,permarr, stat=ier)
    if(ier/=0) call die(myname_,'deallocate gpoints,permarr',ier)

  else
    call GlobalSegMap_copy(inGSMap,GSMap)
  endif

  if (.not. GSMap_increasing(inRGSMap)) then
    if(myPid == 0) call warn(myname_,'RGSMap indices not increasing...Will correct')
    call GlobalSegMap_OPoints(inRGSMap,myPid,gpoints)
    gmapsize=ProcessStorage(inRGSMap,myPid)
    allocate(rpermarr(gmapsize), stat=ier)
    if(ier/=0) call die(myname_,'allocate rpermarr',ier)
    call IndexSet(rpermarr)
    call IndexSort(rpermarr,gpoints)
    call Permute(gpoints,rpermarr,gmapsize)

    call GlobalSegMap_init(RGSMap,gpoints,mycomm,inRGSMap%comp_id,gsize=inRGSMap%gsize)

    deallocate(gpoints,rpermarr, stat=ier)
    if(ier/=0) call die(myname_,'deallocate gpoints,rpermarr',ier)
  else
    call GlobalSegMap_copy(inRGSMap,RGSMap)
  endif


#if 0
   if (t_initialized == 0) then        ! rml timers
     t_initialized=1                   ! rml timers

     call shr_timer_get(t_loop,"m_Router:initp_ loop")  ! rml timers
     call shr_timer_get(t_loop2,"m_Router:initp_ loop2")  ! rml timers
     call shr_timer_get(t_load,"m_Router:initp_ load")  ! rml timers
   endif
#endif

  mysize = ProcessStorage(GSMap,myPid)
  othercomp = GSMap_comp_id(RGSMap)


!.  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  


!!!!  call zeit_ci('t_loop2')
!  call shr_timer_start(t_loop2)    ! rml timers


!! 
!! determine the global segments on this processor 
!! just once, so the info be used repeatedly below
!! same code was used in m_GlobalToLocal - should make a subroutine...
!!

  ngseg = GlobalSegMap_ngseg(GSMap)
  nlseg = GlobalSegMap_nlseg(GSMap, myPid)

  allocate( mygs_lb(nlseg), mygs_ub(nlseg), mygs_len(nlseg), &
            mygs_lstart(nlseg), stat=ier )
  if(ier/=0) call die(myname_,'allocate mygs',ier)

  n = 0
  do i=1,ngseg
    if (GSMap%pe_loc(i) == myPid ) then
      n=n+1
      mygs_lb(n)=GSMap%start(i)
      mygs_ub(n)=GSMap%start(i) + GSMap%length(i) -1
      mygs_len(n)=GSMap%length(i)
    endif
  enddo

  if (n .ne. nlseg) then
    write(stderr,*) myname_,"mismatch nlseg",n,nlseg
    call die(myname)
  endif

  if (nlseg > 0) mygs_lstart(1)=1
  do i=2,nlseg
    mygs_lstart(i)=mygs_lstart(i-1)+mygs_len(i-1)
  enddo


!!
!! determine the segments in RGSMap that are local to each proc
!!

  nprocs=ThisMCTWorld%nprocspid(othercomp)
  r_ngseg = GlobalSegMap_ngseg(RGSMap)

  !! original size of rgs_lb()/ub() was (r_ngseg,nprocs)
  !! at the cost of looping to compute it (within GlobalSegMap_max_nlseg),
  !!   reduced size to (r_max_nlseg,nprocs)
  !! further reduction could be made by flattening it to one dimension
  !!   of size (r_ngseg) and allocating another array to index into it.
  !!   would not improve overall mem use unless this were also done for
  !!   tmpsegstart()/count() and possibly seg_starts()/lengths (the 
  !!   latter would be a major change).

  r_max_nlseg = GlobalSegMap_max_nlseg(RGSMap)

  allocate( rgs_count(nprocs) , &
            rgs_lb(r_max_nlseg,nprocs), rgs_ub(r_max_nlseg,nprocs), &
            nsegs_overlap_arr(nprocs), stat=ier )
  if(ier/=0) call die(myname_,'allocate rgs, nsegs',ier)


  do proc = 1, nprocs

     rgs_count(proc)=0                !! number of segments in RGSMap local to proc

     do i=1,r_ngseg
       if (RGSMap%pe_loc(i) == (proc-1) ) then

         rgs_count(proc) = rgs_count(proc) +1

         rgs_lb( rgs_count(proc) , proc )=RGSMap%start(i)
         rgs_ub( rgs_count(proc) , proc )=RGSMap%start(i) + RGSMap%length(i) -1

       endif
     enddo

  end do

!!! 
!!! this is purely for error checking

  do proc = 1, nprocs
    if (rgs_count(proc) > r_max_nlseg) then
      write(stderr,*) myname_,"overflow on rgs array",proc,rgs_count(proc)
      call die(myname_,'overflow on rgs',0)
    endif
  enddo

!!!


!!!!!!!!!!!!!!!!!!

! allocate space for searching
!   overlap segments to a given remote proc cannot be more than
!   the max of the local segments and the remote segments

  max_rgs_count=0
  do proc=1,nprocs
    max_rgs_count = max( max_rgs_count, rgs_count(proc) )
  enddo

  max_overlap_segs = max(nlseg,max_rgs_count)

  allocate(tmpsegcount(ThisMCTWorld%nprocspid(othercomp), max_overlap_segs),&
           tmpsegstart(ThisMCTWorld%nprocspid(othercomp), max_overlap_segs),&
 	   tmppe_list(ThisMCTWorld%nprocspid(othercomp)),stat=ier)
  if(ier/=0)  &
    call die( myname_,'allocate tmpsegcount etc. size ', &
              ThisMCTWorld%nprocspid(othercomp), &
              ' by ',max_overlap_segs)


  tmpsegcount=0
  tmpsegstart=0
  count =0
  maxsegcount=0

!!!!!!!!!!!!!!!!!!


  do proc = 1, nprocs
    nsegs_overlap = 0
    tmppe_list(proc) = .FALSE.         ! no overlaps with proc yet

    if ( rgs_count(proc) > 0 ) then
      do myseg = 1, nlseg                  ! loop over local segs on 'myPID'

        my_left = mygs_lb(myseg)
        my_right= mygs_ub(myseg)

        do rseg = 1, rgs_count(proc)       !   loop over remote segs on 'proc'

          r_left  = rgs_lb(rseg,proc)
          r_right = rgs_ub(rseg,proc)

          if (.not. (my_right < r_left   .or.  &     ! overlap
                     my_left  > r_right) ) then

             if (nsegs_overlap == 0) then        ! first overlap w/this proc
               count = count + 1
               tmppe_list(proc) = .TRUE.
               prev_right = -9999
             else
               prev_right = local_right
             endif

             overlap_left=max(my_left, r_left)
             overlap_right=min(my_right, r_right)
             overlap_diff= overlap_right - overlap_left

             local_left  = mygs_lstart(myseg) + (overlap_left - my_left)
             local_right = local_left + overlap_diff

                                                   ! non-contiguous w/prev one
             if (local_left /= (prev_right+1) ) then
               nsegs_overlap = nsegs_overlap + 1
               tmpsegstart(count, nsegs_overlap) = local_left
             endif

             tmpsegcount(count, nsegs_overlap) = &
               tmpsegcount(count, nsegs_overlap) + overlap_diff + 1

           endif
        enddo
      enddo
    endif

    nsegs_overlap_arr(proc)=nsegs_overlap
  enddo

  !! pull this out of the loop to vectorize
  do proc=1,nprocs
    maxsegcount=max(maxsegcount,nsegs_overlap_arr(proc))
  enddo


  if (maxsegcount > max_overlap_segs) &
    call die( myname_,'overran max_overlap_segs =', &
              max_overlap_segs, ' count = ',maxsegcount)

!  write(stderr,*) 'max_overlap_segs =', max_overlap_segs, &
!                  'maxsegcount =',maxsegcount, &
!                  'mysize =',mysize


  deallocate( mygs_lb, mygs_ub, mygs_len, mygs_lstart, &
              rgs_count, rgs_lb, rgs_ub, &
              nsegs_overlap_arr, stat=ier)
  if(ier/=0) call die(myname_,'deallocate mygs,rgs,nsegs',ier)


!  call shr_timer_stop(t_loop2)    ! rml timers
!!!!  call zeit_co('t_loop2')
!!!!  call zeit_flush(6)
! this doesn't seem to work
!  call zeit_allflush(mycomm,0,6)



!.  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  


!!!!!!!!!!!!!!!!!!!!end of search through remote GSMap

! start loading up the Router with data

!  call shr_timer_start(t_load)    ! rml timers

  Rout%comp1id = GSMap_comp_id(GSMap)
  Rout%comp2id = othercomp
  Rout%nprocs = count
  Rout%numiatt = 0
  Rout%numratt = 0

  allocate(Rout%pe_list(count),Rout%num_segs(count), &
    Rout%seg_starts(count,maxsegcount), &
    Rout%seg_lengths(count,maxsegcount), &
    Rout%locsize(count),stat=ier)
  if(ier/=0) call die(myname_,'allocate(Rout..)',ier)

  allocate(Rout%istatus(MP_STATUS_SIZE,count), &
             Rout%rstatus(MP_STATUS_SIZE,count), &
	     Rout%rreqs(count),Rout%ireqs(count),stat=ier)
  if(ier/=0) call die(myname_,'allocate(status,reqs,...)',ier)

! allocate the number of pointers needed
  allocate(Rout%ip1(count),stat=ier)
  if(ier/=0) call die(myname_,'allocate(ip1)',ier)

! allocate the number of pointers needed
  allocate(Rout%rp1(count),stat=ier)
  if(ier/=0) call die(myname_,'allocate(rp1)',ier)


    
  m=0
  do i=1,ThisMCTWorld%nprocspid(othercomp)
      if(tmppe_list(i))then 
      m=m+1
      ! load processor rank in MCT_comm
      Rout%pe_list(m)=ThisMCTWorld%idGprocid(othercomp,i-1)
      endif
    enddo

    lmaxsize=0
    do i=1,count
      totallength=0
      do j=1,maxsegcount
	if(tmpsegcount(i,j) /= 0) then
	 Rout%num_segs(i)=j
 	 Rout%seg_starts(i,j)=tmpsegstart(i,j)
	 Rout%seg_lengths(i,j)=tmpsegcount(i,j)
	 totallength=totallength+Rout%seg_lengths(i,j)
	endif
      enddo
      Rout%locsize(i)=totallength
      lmaxsize=MAX(lmaxsize,totallength)
    enddo

    Rout%maxsize=lmaxsize
    Rout%lAvsize=mysize

      
  deallocate(tmpsegstart,tmpsegcount,tmppe_list,stat=ier)
  if(ier/=0) call die(myname_,'deallocate()',ier)

  call GlobalSegMap_clean(RGSMap)
  call GlobalSegMap_clean(GSMap)


!  call shr_timer_stop(t_load)    ! rml timers

 end subroutine initp_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: clean_ - Destroy a Router
!
! !DESCRIPTION:
! Deallocate Router internal data structures and set integer parts to zero.
!
! !INTERFACE:

    subroutine clean_(Rout,stat)
!
! !USES:
!
      use m_die

      implicit none

!INPUT/OUTPUT PARAMETERS:
      type(Router),      intent(inout) :: Rout

!OUTPUT PARAMETERS:
      integer, optional, intent(out)   :: stat

! !REVISION HISTORY:
! 15Jan01 - R. Jacob <jacob@mcs.anl.gov> - initial prototype
! 08Feb01 - R. Jacob <jacob@mcs.anl.gov> - add code to clean
!           the maxsize and locsize
! 01Mar02 - E.T. Ong <eong@mcs.anl.gov> removed the die to prevent
!           crashes and added stat argument.
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::clean_'
  integer :: ier

  deallocate(Rout%pe_list,Rout%num_segs,Rout%seg_starts, &
  Rout%locsize,Rout%seg_lengths,stat=ier)
  if(present(stat)) then
     stat=ier
  else
     if(ier /= 0) call warn(myname_,'deallocate(Rout%pe_list,...)',ier)
  endif

  deallocate(Rout%rreqs,Rout%ireqs,Rout%rstatus,&
   Rout%istatus,stat=ier)
  if(present(stat)) then
     stat=ier
  else
     if(ier /= 0) call warn(myname_,'deallocate(Rout%rreqs,...)',ier)
  endif

  deallocate(Rout%ip1,Rout%rp1,stat=ier)
  if(present(stat)) then
     stat=ier
  else
     if(ier /= 0) call warn(myname_,'deallocate(Rout%ip1,...)',ier)
  endif

  if(associated(Rout%permarr)) then
     deallocate(Rout%permarr,stat=ier)
     if(present(stat)) then
        stat=ier
     else
        if(ier /= 0) call warn(myname_,'deallocate(Rout%ip1,...)',ier)
     endif
  endif

  Rout%comp1id = 0
  Rout%comp2id = 0
  Rout%nprocs = 0
  Rout%maxsize = 0
  Rout%lAvsize = 0
  Rout%numiatt = 0
  Rout%numratt = 0


 end subroutine clean_


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: print_ - Print router info
!
! !DESCRIPTION:
! Print out communication info about router on unit number 'lun'
! e.g. (source,destination,length)
!
! !INTERFACE:

    subroutine print_(rout,mycomm,lun)
!
! !USES:
!
      use m_die
      use m_mpif90

      implicit none

!INPUT/OUTPUT PARAMETERS:
      type(Router),      intent(in) :: Rout
      integer, intent(in)           :: mycomm
      integer, intent(in)           :: lun

! !REVISION HISTORY:
! 27Jul07 - R. Loy <rloy@mcs.anl.gov>  initial version
!EOP ___________________________________________________________________


    integer iproc
    integer myrank
    integer ier
    character(len=*),parameter :: myname_=myname//'::print_'
  
    call MP_comm_rank(mycomm,myrank,ier)
    if(ier/=0) call MP_perr_die(myname_,'MP_comm_rank',ier)


    do iproc=1,rout%nprocs
      if (rout%num_segs(iproc) > 0) then
        write(lun,*) myrank,rout%pe_list(iproc),rout%locsize(iproc)
      endif
    end do        


  end subroutine print_


 end module m_Router

