 MODULE toolmpi

!!======================================================================
!!                   ***  MODULE toolmpi  ***
!!
!! NPROCS= total no of procs using  =>MYPROCS
!! NPROCJ= number of division on J direction =>MYPROCJ
!! NPROCI= number of division on I direction =>MYPROCI
!! IGROUNDPOINT= BLANK point, i.e. no mpi process is used for that local area
!! NBBLANK=COUNT(PROC==IGROUDPOINT). i.e. total number of non-wet points
!! IWETPOINT= wet point i.e. COUNT(PROC==IWETPOINT)==NPROCS
!! NPROCI * NPROCJ=  NPROCS+NBBLANK
!!
!!======================================================================





























































!! * Modules used
   USE parameters, ONLY : rlg

   USE comvars2d,  ONLY : type_mpi_rsh, type_mpi_rlg
   USE mpi


   IMPLICIT NONE


   PRIVATE

!! * Accessibility
   PUBLIC :: MPI_START1 ,MPI_END1,MPI_START2,MPI_START3
   PUBLIC :: ADD_TAG
   PUBLIC :: MPI_COMM_MARS
   PUBLIC :: MASTER,IPROC,NPROCS
   PUBLIC :: give_ccount,get_ccount,MINcBUDY,MAXcBUDY,cBUDY,MINCPROC,MAXCPROC
   PUBLIC :: give_rcount,get_rcount,MINrBUDY,MAXrBUDY,rBUDY,MINRPROC,MAXRPROC
   PUBLIC :: myiproc ,MYRPROC
   PUBLIC :: myjproc ,MYCPROC
! public:: myi_uproc,myi_dproc,myj_uproc,myj_dproc
   PUBLIC :: LI,LJ,CI,RJ,proc
   PUBLIC :: ciminp2,cimaxm1
   PUBLIC :: rjminp2,rjmaxm1

   PUBLIC :: exchange_rsh ,exchange_i_rsh

   PUBLIC :: ex_j_rsh,ex_j_rlg,ex_j_int
   PUBLIC :: ex_i_rsh,ex_i_rlg,ex_i_int
   PUBLIC :: ex_peri_x_rsh,ex_peri_y_rsh
   PUBLIC :: ex_peri_x_avcv,ex_peri_y_avcv
   PUBLIC :: exchange_vectcpu_int



   PUBLIC :: correct_rsh
   PUBLIC :: sendrecv_rlg,sendrecv_rsh


   PUBLIC :: FINDMAX_MPI_6

   PUBLIC :: BCAST_MPI,FINDMAX_MPI,FINDMAX_MPI_3
   PUBLIC :: ADD_ALL_MPI_RSH,ADD_ALL_MPI_INT
   PUBLIC :: MAX_ALL_MPI_RSH,MAX_ALL_MPI_INT
   PUBLIC :: MIN_ALL_MPI_RSH,MIN_ALL_MPI_INT
   PUBLIC :: ADD_MPI_RSH,ADD_MPI_RLG
   PUBLIC :: CREATE_COMM_ZONAL,ADD_ZONAL_RSH
!PUBLIC :: ADD_MPI_INT
! to check balance between mpi rank (see tooltiming)
   REAL(KIND=rlg), DIMENSION(2), PUBLIC  :: xt_mpi_com ,xt_mpi_wait
   REAL(KIND=rlg), DIMENSION(2), PRIVATE :: ztime1,ztime2,ztime3,ztime4

!! * Interface
   INTERFACE BCAST_MPI
     MODULE PROCEDURE BCAST_MPI_LOGIC
     MODULE PROCEDURE BCAST_MPI_INT
     MODULE PROCEDURE BCAST_MPI_RSH
     MODULE PROCEDURE BCAST_MPI_RLG
     MODULE PROCEDURE BCAST_MPI_CHAR
   END INTERFACE



!! * Private variables
   INTEGER                            :: MPI_COMM_MARS
   INTEGER                            :: NPROCS,IPROC,MAXP,MINP
   INTEGER                            :: MINIPROC,MAXIPROC,MINJPROC,maxjproc,NBBLANK
   INTEGER                            :: MYIPROC,MYJPROC,MYCPROC,MYRPROC
   INTEGER                            :: MINCPROC,MAXCPROC,MINRPROC,MAXRPROC
   INTEGER                            :: MINcBUDY,MAXcBUDY,MINRBUDY,MAXRBUDY
   INTEGER                            :: MINiBUDY,MAXiBUDY,MINjBUDY,MAXjBUDY
   INTEGER                            :: NUMRBUDY,NUMCBUDY
   INTEGER                            :: NUMiBUDY,NUMjBUDY
   INTEGER                            :: ciminp2,cimaxm1 !,cimin,cimax
   INTEGER                            :: rjminp2,rjmaxm1
   INTEGER,PARAMETER                  :: IGROUNDPOINT=-100
   INTEGER,PARAMETER                  :: IWETPOINT=-2
   INTEGER                            :: myi_uproc,myi_dproc,myj_uproc,myj_dproc
   INTEGER                            :: impierror=6
   INTEGER                            :: njbrecv, njhrecv,njbsend,njhsend
   LOGICAL                            :: MASTER=.FALSE.,LPARAMREAD=.FALSE.
   INTEGER,ALLOCATABLE,DIMENSION(:,:) :: LJ,CI,LI,RJ
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: cbudy,give_ccount,get_ccount
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: rbudy,give_rcount,get_rcount
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: ibudy,jbudy
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: proc_i,proc_j,proc_size,proc_displs
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: proc_limin,proc_limax,proc_ljmin,proc_ljmax
   INTEGER,ALLOCATABLE,DIMENSION(:,:) :: PROC

CONTAINS






!!======================================================================

  SUBROUTINE init_mpi_type_size ()
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE  init_mpi_type_size  ***
!&E
!&E ** Purpose : Defines the MPI type for 'rsh' and 'rlg' kind real
!&E              variables.
!&E
!&E       !  2014-12 (M. Honnorat) Refactor MPI routines
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh, rlg

   IMPLICIT NONE

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part
   IF      ( rsh == 4 ) THEN
      type_mpi_rsh = MPI_REAL
   ELSE IF ( rsh == 8 ) THEN
      type_mpi_rsh = MPI_DOUBLE_PRECISION
   ELSE
      WRITE(impierror,*)'MPI exchange works only for rsh=4 and 8, you gave rsh=',rsh
      CALL MPI_FINALIZE(ierr_mpi) ; STOP
   ENDIF

   IF      ( rlg == 4 ) THEN
      type_mpi_rlg = MPI_REAL
   ELSE IF ( rlg == 8 ) THEN
      type_mpi_rlg = MPI_DOUBLE_PRECISION
   ELSE
      WRITE(impierror,*)'MPI exchange works only for rlg=4 and 8, you gave rlg=',rlg
      CALL MPI_FINALIZE(ierr_mpi) ; STOP
   ENDIF

  END SUBROUTINE init_mpi_type_size

!!======================================================================

  SUBROUTINE MPI_START1

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE  MPI_START1  ***
!&E
!&E ** Purpose :
!&E
!&E       !  2013-07 (A. Thevenin - CERFACS) coupling using OASIS3-MCT
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : imax,jmax,l_testcase,num_testcase



   IMPLICIT NONE

!! * Local declarations
   INTEGER :: ierr_mpi
   LOGICAL :: ios, ioshyb
   CHARACTER(len=100) :: filempitxt

!!----------------------------------------------------------------------
!! * Executable part

!Here are the list of trapp



     CALL MPI_INIT(ierr_mpi)

     MPI_COMM_MARS = MPI_COMM_WORLD



   CALL init_mpi_type_size()
   CALL MPI_COMM_SIZE ( MPI_COMM_MARS, NPROCS, ierr_mpi )
   CALL MPI_COMM_RANK ( MPI_COMM_MARS, IPROC, ierr_mpi )
   MINP=0
   MAXP=NPROCS-1
   IF (IPROC.EQ.0) MASTER=.TRUE.
   IF (MASTER) WRITE(*,*) IPROC,"I AM THE MASTER"

   MINIPROC=0
   MINJPROC=0
   NBBLANK=0

   filempitxt='mpi.txt'
   !IF (.NOT. .TRUE.) filempitxt='mpi'//TRIM(Agrif_Cfixed())//'.txt'
   INQUIRE(file     = filempitxt, EXIST = ios)
   filempitxt='hybrid'
   INQUIRE(file     = filempitxt, EXIST = ioshyb)

   IF (l_testcase .AND. .NOT.ioshyb) THEN

     MAXJPROC = -1

     IF (num_testcase==3 .OR. num_testcase==11) THEN
       MAXIPROC=1
       MAXJPROC=1
     ELSE IF (num_testcase==2) THEN
       MAXIPROC=7
       MAXJPROC=0
     ELSE IF (num_testcase==17) THEN
       MAXIPROC=2
       MAXJPROC=2
     ELSE IF (num_testcase==13) THEN   ! casacc
       MAXIPROC=1
       MAXJPROC=1

     ELSE IF (num_testcase==26) THEN
       MAXIPROC=3
       MAXJPROC=0
     ELSE IF (num_testcase==35) THEN
       MAXIPROC=1   ! pour 4 procs
       MAXJPROC=1
!       MAXIPROC=3  ! pour 16 procs
!       MAXJPROC=3
     ELSE IF (num_testcase==37) THEN

       MAXIPROC=1
       MAXJPROC=0

     ELSE
       MAXIPROC = 3
       MAXJPROC = NPROCS/(MAXIPROC+1) - 1
     END IF

     IF(MAXJPROC == -1) THEN
       MAXJPROC = NPROCS/(MAXIPROC+1) - 1
     ENDIF

     ALLOCATE(LI(4,miniproc:maxiproc))
     ALLOCATE(LJ(4,minjproc:maxjproc))
     ALLOCATE(proc(miniproc:maxiproc,minjproc:maxjproc))
     proc=-2

   ELSE IF (.NOT.ios .AND. .NOT.ioshyb) then

     WRITE(*,*)'mpi.txt does not exit'
     WRITE(*,*)'mpi job will go to 1D mode '

!MAXIPROC=3
!MAXJPROC=3
! MAXIPROC=15
! MAXJPROC=0
!IF(nprocs.NE.16)STOP !if(nprocs.ne.25)stop
     IF (jmax.gt.imax) THEN
       MAXIPROC=0
       MAXJPROC=NPROCS-1
     ELSE
       MAXIPROC=NPROCS-1
       MAXJPROC=0
     END IF

     ALLOCATE(LI(4,miniproc:maxiproc))
     ALLOCATE(LJ(4,minjproc:maxjproc))
     ALLOCATE(proc(miniproc:maxiproc,minjproc:maxjproc))
     proc=IWETPOINT

   ELSE IF (.NOT.ios .AND. ioshyb) then

     WRITE(*,*)'openMP and MPI job'

     IF (num_testcase==13) THEN   ! casacc
       MAXIPROC=1
       MAXJPROC=1
     ELSE
       WRITE(*,*)'mpi job will go to 1D mode (along x)'
       MAXIPROC=NPROCS-1
       MAXJPROC=0
     END IF

     ALLOCATE(LI(4,miniproc:maxiproc))
     ALLOCATE(LJ(4,minjproc:maxjproc))
     ALLOCATE(proc(miniproc:maxiproc,minjproc:maxjproc))
     proc=IWETPOINT

   ELSE
     IF (MASTER) THEN
       WRITE(*,*)'mpi.txt  exist'
       WRITE(*,*)'read the mpi.txt'
     END IF
     LPARAMREAD=.TRUE.
     call MPI_PARAMREAD
   END IF



   CALL DEFINE_PROC
   CALL DEFINE_UDPROC



   RETURN

  END SUBROUTINE MPI_START1

!!======================================================================

  SUBROUTINE MPI_END1

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE  MPI_END1  ***
!&E
!&E ** Purpose :
!&E
!&E       !  2013-07 (A. Thevenin - CERFACS) coupling using OASIS3-MCT
!&E
!&E---------------------------------------------------------------------
!! * Modules used

   IMPLICIT NONE

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

   DEALLOCATE(cbudy,LI,LJ,CI,give_ccount,get_ccount)
   DEALLOCATE(rbudy,give_rcount,get_rcount)
   DEALLOCATE(PROC,ibudy,jbudy)
   DEALLOCATE(proc_i,proc_j,proc_size,proc_displs)
   DEALLOCATE(proc_limin,proc_limax,proc_ljmin,proc_ljmax)




   CALL MPI_FINALIZE(ierr_mpi)


   RETURN

  END SUBROUTINE MPI_END1

!!======================================================================

  SUBROUTINE MPI_START2(CIP,LJP)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE MPI_START2  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : imin,imax
   USE parameters, ONLY : limin,limax,liminp1,liminp2,limaxm1,liminm1,limaxp1,liminm2,limaxp2
   USE parameters, ONLY : limaxm1p1,liminp1m1,liminp1m2,ljminp1m2,limaxp3
   USE parameters, ONLY : jmin,jmax
   USE parameters, ONLY : ljmin,ljmax,ljminp1,ljminp2,ljmaxm1,ljmaxm1p1,ljminp1m1
   USE parameters, ONLY : ljminm2,ljmaxp2,ljminm1,ljmaxp1,ljmaxp3
   USE parameters, ONLY : ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north

   IMPLICIT NONE

!! * Arguments
   INTEGER, INTENT(in) :: CIP,LJP

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: i,j,n
   INTEGER :: LIP,RJP,maxfakeproc,countblank
   INTEGER,DIMENSION(imin:imax) :: ISUM
   INTEGER,DIMENSION(jmin:jmax) :: JSUM
   CHARACTER(len=7)  :: tag
   CHARACTER(len=30) :: filename

!!----------------------------------------------------------------------
!! * Executable part

   LIP=LJP
   RJP=CIP

   IF (.not.LPARAMREAD) THEN
! DEFINE LI
     IF (LIP.EQ.0) THEN
       CALL divide_simple1d(imin,imax,miniproc,maxiproc,maxiproc,LI)
     ELSE
       CALL get_isum(isum)
       CALL divide_withsum1d(imin,imax,miniproc,maxiproc,maxiproc,isum,LI)
     END IF
! DEFINE LJ
     IF (LJP.EQ.0) THEN
       CALL divide_simple1d(jmin,jmax,minjproc,maxjproc,maxjproc,LJ)
     ELSE
       CALL get_jsum(jsum)
       CALL divide_withsum1d(jmin,jmax,minjproc,maxjproc,maxjproc,jsum,LJ)
     END IF
   END IF

   limin=LI(1,myiproc)
   liminp1=limin
   liminp2=limin
   liminm2=limin -2
   liminm1=limin -1
   liminp1m1=limin -1
   liminp1m2=limin -2
   limax=LI(2,myiproc)
   limaxm1=limax
   limaxp2=limax +2
   limaxp3=limax +3
   limaxp1=limax +1
   limaxm1p1=limax +1
   IF (myiproc.EQ.miniproc) THEN
     liminp1=imin+1+ishift_obc_west
     liminp2=imin+2+ishift_obc_west
     liminm2=imin
     liminm1=imin
     liminp1m1=imin+1+ishift_obc_west
     liminp1m2=imin+1+ishift_obc_west
   END IF
   IF (myiproc.EQ.maxiproc) THEN
     limaxm1=imax-1-ishift_obc_east
     limaxp2=imax
     limaxp3=imax
     limaxp1=imax
     limaxm1p1=limax-1-ishift_obc_east
   END IF

   ljmin=LJ(1,myjproc)
   ljminp1=ljmin
   ljminp2=ljmin
   ljminm2=ljmin -2
   ljminm1=ljmin -1
   ljminp1m1=ljmin -1
   ljminp1m2=ljmin -2
   ljmax=LJ(2,myjproc)
   ljmaxm1=ljmax
   ljmaxp2=ljmax +2
   ljmaxp3=ljmax +3
   ljmaxp1=ljmax +1
   ljmaxm1p1=ljmax +1
   IF (myjproc.EQ.minjproc) THEN
     ljminp1=jmin+1+jshift_obc_south
     ljminp2=jmin+2+jshift_obc_south
     ljminm2=jmin
     ljminm1=jmin
     ljminp1m1=jmin+1+jshift_obc_south
     ljminp1m2=jmin+1+jshift_obc_south
   END IF
   IF (myjproc.EQ.maxjproc) THEN
     ljmaxm1=jmax-1-jshift_obc_north
     ljmaxp2=jmax
     ljmaxp3=jmax
     ljmaxp1=jmax
     ljmaxm1p1=ljmax-1-jshift_obc_north
   END IF

! Make the shift of the OBC is not too large
! To simplify the exchanges, the OBCs must be inside the cpus
! delimitating the domain
!cval FAUT-IL PRENDRE PLUS DE MARGE ?
! il faut ig+3 dans miniproc (> suffisant)
   IF (myiproc.EQ.miniproc) THEN
     IF ((ishift_obc_west+3)>=limax) THEN
       WRITE(impierror,*) 'ishift_obc_west ,',ishift_obc_west
       WRITE(impierror,*) ', is larger than'
       WRITE(impierror,*) 'limax, ', limax
       WRITE(impierror,*) 'you MUST HAVE ishift_obc_west+3 < limax'
       CALL MPI_FINALIZE(ierr_mpi) ; STOP
     END IF
   END IF
! il faut 3 a gauche de id dans maxiproc (> suffisant)
   IF (myiproc.EQ.maxiproc) THEN
     IF (ishift_obc_east+3>=(limax-limin)) THEN
       WRITE(impierror,*) 'ishift_obc_east ,',ishift_obc_east
       WRITE(impierror,*) ', is larger than the size of the cpu'
       WRITE(impierror,*) 'bounding the eastern boundary , ', limax-limin
       WRITE(impierror,*) 'you MUST HAVE ishift_obc_east < (limax-limin)'
       CALL MPI_FINALIZE(ierr_mpi) ; STOP
     END IF
   END IF
! il faut jb+3 dans minjproc (> suffisant)
   IF (myjproc.EQ.minjproc) THEN
     IF ((jshift_obc_south+3)>=ljmax) THEN
       WRITE(impierror,*) 'jshift_obc_south ,',jshift_obc_south
       WRITE(impierror,*) ', is larger than'
       WRITE(impierror,*) 'ljmax, ', ljmax
       WRITE(impierror,*) 'you MUST HAVE jshift_obc_south+3 < ljmax'
       CALL MPI_FINALIZE(ierr_mpi) ; STOP
     END IF
   END IF
! il faut 3 en bas de jh dans maxjproc (> suffisant)
   IF (myjproc.EQ.maxjproc) THEN
     IF (jshift_obc_north+3>=(ljmax-ljmin)) THEN
       WRITE(impierror,*) 'jshift_obc_north ,',jshift_obc_north
       WRITE(impierror,*) ', is larger than the size of the cpu'
       WRITE(impierror,*) 'bounding the northern boundary , ', ljmax-ljmin
       WRITE(impierror,*) 'you MUST HAVE jshift_obc_north < (ljmax-ljmin)'
       CALL MPI_FINALIZE(ierr_mpi) ; STOP
     END IF
   END IF

! define CI here

   ALLOCATE(CI(4,mincproc:maxcproc))
   CI=0
   COUNTBLANK=COUNT(PROC(MYIPROC,mincproc:maxcproc)==IGROUNDPOINT)
   maxfakeproc=maxcproc-countblank

   IF (CIP.EQ.0) THEN
     CALL divide_simple1d(liminp2,limaxm1,mincproc,maxcproc,maxfakeproc,CI)
   ELSE
     CALL get_isum(isum)
     CALL divide_withsum1d(liminp2,limaxm1,mincproc,maxcproc,maxfakeproc,isum(liminp2:limaxm1),CI)
   END IF
   if(countblank.ne.0)CALL shiftforblank(CI,mincproc,maxcproc,PROC(myiproc,mincproc:maxcproc))


! define RJ here


   ALLOCATE(RJ(4,minrproc:maxrproc))
   RJ=0
   COUNTBLANK=COUNT(PROC(minrproc:maxrproc,myjproc)==IGROUNDPOINT)
   maxfakeproc=maxrproc-countblank
   IF (RJP.EQ.0) THEN
     CALL divide_simple1d(ljminp2,ljmaxm1,minrproc,maxrproc,maxfakeproc,RJ)
   ELSE
     CALL get_jsum(jsum)
     CALL divide_withsum1d(ljminp2,ljmaxm1,minrproc,maxrproc,maxfakeproc,jsum(ljminp2:ljmaxm1),RJ)
   END IF

   if(countblank.ne.0)CALL shiftforblank(RJ,minrproc,maxrproc,PROC(minrproc:maxrproc,myjproc))

   ciminp2=CI(1,mycproc)
   cimaxm1=CI(2,mycproc)
   rjminp2=RJ(1,myrproc)
   rjmaxm1=RJ(2,myrproc)

   IF (LI(2,maxiproc).NE.imax) WRITE(impierror,*) imax,LI(2,maxiproc),"MPI_BUG1 report to tina ASAP" !;stop
   IF (LJ(2,maxjproc).NE.jmax) WRITE(impierror,*) jmax,LJ(2,maxjproc),"MPI_BUG2 report to tina ASAP" !;stop
   IF (CI(2,maxcproc).NE.limaxm1) WRITE(impierror,*) limaxm1,CI(2,maxcproc),"MPI_BUG3 report to tina ASAP "!;stop
   IF (RJ(2,maxrproc).NE.ljmaxm1) WRITE(impierror,*) ljmaxm1,RJ(2,maxrproc),"MPI_BUG4 report to tina ASAP",iproc !;stop

! cval val2tina : any MPI job is stopped with the following lines why ?
!   IF ( LI(2,maxiproc).NE.imax .OR. LJ(2,maxjproc).NE.jmax .OR. CI(2,maxcproc).NE.limaxm1 .OR. RJ(2,maxrproc).NE.ljmaxm1) &
!     CALL MPI_FINALIZE(ierr_mpi) ; STOP

   ALLOCATE(proc_i(minp:maxp),proc_j(minp:maxp))
   ALLOCATE(proc_limin(minp:maxp),proc_ljmin(minp:maxp))
   ALLOCATE(proc_limax(minp:maxp),proc_ljmax(minp:maxp))
   ALLOCATE(proc_size(minp:maxp),proc_displs(minp:maxp+1))

   DO j=minjproc,maxjproc
     DO i=miniproc,maxiproc
       n=proc(i,j)
       IF (n.NE.IGROUNDPOINT) THEN
         proc_i(n)=i
         proc_j(n)=j
         proc_size(n)=LI(3,i)*LJ(3,j)
         proc_limin(n)=LI(1,i)
         proc_limax(n)=LI(2,i)
         proc_ljmin(n)=LJ(1,j)
         proc_ljmax(n)=LJ(2,j)
       END IF
     END DO
   END DO
   proc_displs=0
   DO n=minp+1,maxp+1
     proc_displs(n)=proc_displs(n-1)+proc_size(n-1)
   END DO
   NBBLANK=COUNT(PROC==IGROUNDPOINT)

   filename='mpierror'
   CALL TAG_PROC_IJ(tag)
   filename=TRIM(filename)//'.'//tag
   impierror=5123


   CALL MPI_DEBUG_PRINTPARAMS2

   RETURN

  END SUBROUTINE MPI_START2

!!======================================================================

  SUBROUTINE MPI_START3

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE MPI_START3  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY: ljmin,ljmax
   USE parameters, ONLY: limin,limax
   USE comvars2d , ONLY: jb,jh,ig,id
   IMPLICIT NONE

!! * Local declarations
   INTEGER :: x,y,i,j,k

!!----------------------------------------------------------------------
!! * Executable part

   ALLOCATE(give_ccount(MINcBUDY:MAXcBUDY))
   ALLOCATE(get_ccount(MINcBUDY:MAXcBUDY))

   give_ccount=0
   get_ccount=0

   A: DO x=MINcBUDY,MAXcBUDY
     y=cbudy(x)
     IF (y.EQ.mycproc) CYCLE A
     IF (y.LT.mincproc) CYCLE A
     IF (y.GT.maxcproc) CYCLE A
     k=0
     DO j=ljmin,ljmax
       DO i=CI(1,y),CI(2,y)
         IF (j.LT.jb(i))CYCLE
         IF (j.GT.jh(i))CYCLE
         k=k+1
       END DO
     END DO
     give_ccount(x)=k

     k=0
     DO j=LJ(1,y),LJ(2,y)
       DO i=ciminp2,cimaxm1
         IF (j.LT.jb(i))CYCLE
         IF (j.GT.jh(i))CYCLE
         k=k+1
       END DO
     END DO
     get_ccount(x)=k
   END DO A

   ALLOCATE(give_rcount(MINrBUDY:MAXrBUDY))
   ALLOCATE(get_rcount(MINrBUDY:MAXrBUDY))
   give_rcount=0
   get_rcount=0

   B: DO x=MINrBUDY,MAXrBUDY
     y=rbudy(x)
     IF (y.EQ.myrproc) CYCLE B
     IF (y.LT.minrproc) CYCLE B
     IF (y.GT.maxrproc) CYCLE B
     k=0
     DO j=RJ(1,y),RJ(2,y)
       DO i=limin,limax
         IF (i.LT.ig(j))CYCLE
         IF (i.GT.id(j))CYCLE
         k=k+1
       END DO
     END DO
     give_rcount(x)=k

     k=0
     DO j=rjminp2,rjmaxm1
       DO i=LI(1,y),LI(2,y)
         IF (i.LT.ig(j))CYCLE
         IF (i.GT.id(j))CYCLE
         k=k+1
       END DO
     END DO
     get_rcount(x)=k
   END DO B

   if(master) CALL makeprocmap
   if(master) CALL countwetpoint


   RETURN

  END SUBROUTINE MPI_START3

!!======================================================================

  SUBROUTINE MPI_PARAMREAD

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE MPI_PARAMREAD  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Local declarations
   INTEGER :: ierr_mpi
   LOGICAL :: ios
   INTEGER :: ios2,itoto,itotoo,i,j
   INTEGER :: numiproc,numjproc,numactifs,tmpiproc, tmpjproc, tmpimin,tmpimax, tmpjmin, tmpjmax
   INTEGER :: readnproc
   LOGICAL :: LBLANK
   character(len=100) :: filempitxt

!!----------------------------------------------------------------------
!! * Executable part


   filempitxt='mpi.txt'
   !IF (.NOT. .TRUE.) filempitxt='mpi'//TRIM(Agrif_Cfixed())//'.txt'
   INQUIRE(file     = filempitxt, EXIST = ios)
   IF (.NOT.ios) stop "probleme a l ouverture mpi.txt en ecriture"
   OPEN(unit     = 1000, file     = filempitxt, &
        form     = "formatted",       &
        access   = "sequential",      &
        action   = "read",           &
        position = "rewind" )

   READ(unit = 1000, fmt = '(I12, A12)', iostat = ios2) itoto
   READ(unit = 1000, fmt = '(I12, A12)', iostat = ios2) numiproc
   READ(unit = 1000, fmt = '(I12, A12)', iostat = ios2) numjproc
   READ(unit = 1000, fmt = '(I12, A12)', iostat = ios2) itoto
   READ(unit = 1000, fmt = '(I12, A12)', iostat = ios2) numactifs
   miniproc=0
   minjproc=0
   maxiproc=numiproc-1
   maxjproc=numjproc-1
   ALLOCATE(LI(4,miniproc:maxiproc))
   ALLOCATE(LJ(4,minjproc:maxjproc))
   ALLOCATE(proc(miniproc:maxiproc,minjproc:maxjproc))
   READ(unit = 1000, fmt = '(8A12)', iostat = ios2)
   DO j = 1, numjproc
     DO i = 1, numiproc
       READ(unit = 1000, fmt = '(6I12, L12, I12)', iostat = ios2) &
            tmpiproc, tmpjproc, tmpimin,tmpimax, tmpjmin, tmpjmax, LBLANK
!copy tmpimin imax to mars here??
       IF (.NOT.LBLANK) THEN
         proc(tmpiproc,tmpjproc)=IGROUNDPOINT
       ELSE
         proc(tmpiproc,tmpjproc)=IWETPOINT
         LI(1,tmpiproc)=tmpimin
         LI(2,tmpiproc)=tmpimax
         LJ(1,tmpjproc)=tmpjmin
         LJ(2,tmpjproc)=tmpjmax
       END IF
     END DO
   END DO
   LI(3,:)=LI(2,:)-LI(1,:)+1
   LI(4,:)=LI(1,:)-LI(1,miniproc)
   LJ(3,:)=LJ(2,:)-LJ(1,:)+1
   LJ(4,:)=LJ(1,:)-LJ(1,minjproc)

   CLOSE(unit = 1000)

   IF (numactifs.NE.nprocs) THEN
     WRITE(impierror,*)'numiproc,numjproc'
     WRITE(impierror,*)numiproc,numjproc
     WRITE(impierror,*)'stop : numactifs on mpi.txt does not match nprocs : ',nprocs
     CALL MPI_FINALIZE(ierr_mpi) ; STOP
   END IF
!NBBLANK=6   !read number of BLANK
!MAXIPROC=3  !read max number of I grid
!MAXjPROC=2  !read max number of J grid
!ALLOCATE(proc(miniproc:maxiproc,minjproc:maxjproc))
!proc=-2
!proc(0,0)=IGROUNDPOINT
!proc(0,1)=IGROUNDPOINT
!proc(1,2)=IGROUNDPOINT
!proc(2,2)=IGROUNDPOINT
!proc(3,0)=IGROUNDPOINT
!proc(3,2)=IGROUNDPOINT


   RETURN

  END SUBROUTINE MPI_PARAMREAD

!!======================================================================

  SUBROUTINE DEFINE_PROC

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE DEFINE_PROC  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: n,i,j,countblank
   INTEGER,ALLOCATABLE,DIMENSION(:) :: proctmp !this does not work yet
!to avoid using this i need to get rid of exchange_rsh

!!----------------------------------------------------------------------
!! * Executable part

   n=0
   DO j=minjproc,maxjproc
     DO i=miniproc,maxiproc
       IF (proc(i,j).EQ.IWETPOINT) THEN
         proc(i,j)=n
         n=n+1
       END IF
     END DO
   END DO
   IF (n.NE.nprocs) THEN
     WRITE(impierror,*) 'minjproc,maxjproc,miniproc,maxiproc,n,nprocs'
     WRITE(impierror,*) minjproc,maxjproc,miniproc,maxiproc,n,nprocs
     WRITE(impierror,*)'maxjproc and maxiproc do not correspond to np stop'
     CALL MPI_FINALIZE(ierr_mpi) ; STOP
   END IF

! from here, it defines myiproc,myjproc,mincproc,maxcproc,minrproc,maxcproc
! depending on what is your 'iproc' is
! i.e. each process have differfent value!!!
   DO j=minjproc,maxjproc
     DO i=miniproc,maxiproc
       IF (proc(i,j).EQ.iproc) THEN
         myiproc=i
         myjproc=j
       END IF
     END DO
   END DO

   MYCPROC=MYJPROC
   MYRPROC=MYIPROC
   i=myiproc
   mincproc=minjproc
   A:DO j=minjproc,maxjproc
     IF (PROC(i,j).EQ.IGROUNDPOINT) THEN
       mincproc=j+1
     ELSE
       EXIT A
     END IF
   END DO A
   maxcproc=maxjproc
   B:DO j=maxjproc,minjproc,-1
     IF (PROC(i,j).EQ.IGROUNDPOINT) THEN
       maxcproc=j-1
     ELSE
       EXIT B
     END IF
   END DO B

   j=myjproc
   minrproc=miniproc
   C:DO i=miniproc,maxiproc
     IF (PROC(i,j).EQ.IGROUNDPOINT) THEN
       minrproc=i+1
     ELSE
       EXIT C
     END IF
   END DO C

   maxrproc=maxiproc
   D: DO i=maxiproc,miniproc,-1
     IF (PROC(i,j).EQ.IGROUNDPOINT) THEN
       maxrproc=i-1
     ELSE
       EXIT D
     END IF
   END DO D

! for exchanges in 1D
   MINCBUDY=0
   COUNTBLANK=COUNT(PROC(MYIPROC,mincproc:maxcproc)==IGROUNDPOINT)
   MAXCBUDY=MAXCPROC-MINCPROC-COUNTBLANK
   NUMCBUDY=MAXCBUDY+1
   IF ((numcbudy.NE.1).AND.(MOD(numcbudy,2).EQ.1)) THEN
     MAXCBUDY=MAXCBUDY+1
   END IF
   ALLOCATE(Cbudy(MINcBUDY:MAXcBUDY))
   CALL makebudy(Cbudy,maxcbudy,mincproc,maxcproc,numcbudy,mycproc,PROC(myiproc,mincproc:maxcproc))

   MINRBUDY=0
   COUNTBLANK=COUNT(PROC(minrproc:maxrproc,myjproc)==IGROUNDPOINT)
   MAXRBUDY=MAXRPROC-MINRPROC-COUNTBLANK
   NUMRBUDY=MAXRBUDY+1
   IF ((numRbudy.NE.1).AND.(MOD(numRbudy,2).EQ.1)) THEN
     MAXRBUDY=MAXRBUDY+1
   END IF
   ALLOCATE(rbudy(MINrBUDY:MAXrBUDY))
   CALL makebudy(rbudy,maxrbudy,minrproc,maxrproc,numrbudy,myrproc,PROC(minrproc:maxrproc,myjproc))

! for exchanges in 2D
   MINIBUDY=MINIPROC
   MAXIBUDY=MAXIPROC
   NUMiBUDY=MAXIBUDY+1
   IF ((numibudy.NE.1).AND.(MOD(numibudy,2).EQ.1)) THEN
     MAXiBUDY=MAXiBUDY+1
   END IF
   ALLOCATE(ibudy(MINiBUDY:MAXiBUDY))
   ALLOCATE(proctmp(MINiproc:MAXiproc))
   proctmp=1
   CALL makebudy(ibudy,maxibudy,miniproc,maxiproc,numibudy,myiproc,PROCtmp)
   deallocate(proctmp)

   MINjBUDY=MINjPROC
   MAXjBUDY=MAXjPROC
   NUMjBUDY=MAXjBUDY+1
   IF ((numjbudy.NE.1).AND.(MOD(numjbudy,2).EQ.1)) THEN
     MAXjBUDY=MAXjBUDY+1
   END IF
   ALLOCATE(Jbudy(MINjBUDY:MAXjBUDY))
   ALLOCATE(proctmp(MINjproc:MAXjproc))
   proctmp=1
   CALL makebudy(jbudy,maxjbudy,minjproc,maxjproc,numjbudy,myjproc,PROCtmp)
   deallocate(proctmp)

  END SUBROUTINE DEFINE_PROC

!!======================================================================

  SUBROUTINE DEFINE_UDPROC

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE DEFINE_UDPROC  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!!----------------------------------------------------------------------
!! * Executable part

   IF (myjproc==maxjproc) THEN
     myj_uproc=MPI_PROC_NULL
   ELSE IF (PROC(myiproc,myjproc+1).EQ.IGROUNDPOINT) THEN
     myj_uproc=MPI_PROC_NULL
   ELSE
     myj_uproc=PROC(myiproc,myjproc+1)
   END IF

   IF (myjproc==minjproc) THEN
     myj_dproc=MPI_PROC_NULL
   ELSE IF (PROC(myiproc,myjproc-1).EQ.IGROUNDPOINT) THEN
     myj_dproc=MPI_PROC_NULL
   ELSE
     myj_dproc=PROC(myiproc,myjproc-1)
   END IF

   IF (myiproc==maxiproc) THEN
     myi_uproc=MPI_PROC_NULL
   ELSE IF (PROC(myiproc+1,myjproc).EQ.IGROUNDPOINT) THEN
     myi_uproc=MPI_PROC_NULL
   ELSE
     myi_uproc=PROC(myiproc+1,myjproc)
   END IF

   IF (myiproc==miniproc) THEN
     myi_dproc=MPI_PROC_NULL
   ELSE IF (PROC(myiproc-1,myjproc).EQ.IGROUNDPOINT) THEN
     myi_dproc=MPI_PROC_NULL
   ELSE
     myi_dproc=PROC(myiproc-1,myjproc)
   END IF

   RETURN

  END SUBROUTINE DEFINE_UDPROC

!!======================================================================

  SUBROUTINE MPI_DEBUG_PRINTPARAMS1

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE MPI_DEBUG_PRINTPARAMS1  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Local declarations
   INTEGER           :: i,ii
   CHARACTER(len=7)  :: tag
   CHARACTER(len=30) :: filename

!!----------------------------------------------------------------------
!! * Executable part

   CALL TAG_PROC_IJ(tag)
   filename=TRIM('mpidebugoutput')//'.'//tag
   !IF (.NOT. .TRUE.) filename='mpidebugoutput'//'.'//tag//'_'//TRIM(Agrif_Cfixed())
   i=1000
   OPEN(i,file=filename)
   WRITE(i,*)'mpi_debug_printparams1 start'
   WRITE(i,*)'iproc,nprocs',iproc,nprocs
   WRITE(i,*)'maxiproc,maxjproc',maxiproc,maxjproc
   WRITE(i,*)'myiproc,myjproc',myiproc,myjproc
   WRITE(i,*)'print proc(:,j) here'
   DO ii=minjproc,maxjproc
     WRITE(i,*)'j=',ii
     WRITE(i,*)proc(:,ii)
   END DO

   WRITE(i,*)'mincproc,maxcproc',mincproc,maxcproc
   WRITE(i,*)'mincbudy,maxcbudy',mincbudy,maxcbudy
   WRITE(i,*)'cbudy',cbudy
   WRITE(i,*)'minrproc,maxrprro',minrproc,maxrproc
   WRITE(i,*)'minrbudy,maxrbudy',minrbudy,maxrbudy
   WRITE(i,*)'rbudy',rbudy
   WRITE(i,*)'mpi_debug_printparams1 end'
   WRITE(i,*)

   RETURN

  END SUBROUTINE MPI_DEBUG_PRINTPARAMS1

!!======================================================================

  SUBROUTINE MPI_DEBUG_PRINTPARAMS2

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE MPI_DEBUG_PRINTPARAMS2  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters
   IMPLICIT NONE

!! * Local declarations
   INTEGER :: i,ii,iii,j
   CHARACTER(LEN=lchain) :: file_mpidebugprepare

!!----------------------------------------------------------------------
!! * Executable part

   IF (MASTER) THEN
     iii=1100
     file_mpidebugprepare = 'mpidebugprepare'
     !IF (.NOT. .TRUE.) file_mpidebugprepare='mpidebugprepare_'//TRIM(Agrif_Cfixed())

     OPEN(iii,file=file_mpidebugprepare)
           WRITE(iii,800)nprocs
     DO j=minjproc,maxjproc
       DO i=miniproc,maxiproc
         IF (proc(i,j).GE.0) THEN
           WRITE(iii,800)proc(i,j),i,j,LI(1,i),LI(2,i),LJ(1,j),LJ(2,j)
         END IF
       END DO
     END DO

800  FORMAT(7(i10.10,' '))

     CLOSE(iii)
   END IF



  END SUBROUTINE MPI_DEBUG_PRINTPARAMS2

!!======================================================================

  SUBROUTINE MPI_DEBUG_PRINTPARAMS3

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE MPI_DEBUG_PRINTPARAMS3  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Local declarations
   INTEGER :: i,ii

!!----------------------------------------------------------------------
!! * Executable part

   ii=1000

   WRITE(ii,*)'mpi_debug_printparams3 start'
   WRITE(ii,*)'print give_ccount'
   WRITE(ii,*)give_ccount
   WRITE(ii,*)'print get_ccount'
   WRITE(ii,*)get_ccount
   WRITE(ii,*)'print give_rcount'
   WRITE(ii,*)give_rcount
   WRITE(ii,*)'print get_rcount'
   WRITE(ii,*)get_rcount
   WRITE(ii,*)'mpi_debug_printparams3 end'

   CLOSE(ii)
   RETURN

  END SUBROUTINE MPI_DEBUG_PRINTPARAMS3

!!======================================================================

  SUBROUTINE makebudy(bud,maxbudy,minproc,maxproc,numbudy,myproc,aproc)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE makebudy  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description : ATT: BUDY counts the number of times each procs transpose data;
!&E                  budy will let you know who you should exchange the data each N,i,j
!&E                  N starts always from 0, and ends at maxbudy
!&E                  minproc indicates the first process  which starts the transpose data
!&E                  maxproc indicates the last process which ends the transpose data
!&E
!&E ** Modified variables : bud
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER                     ,INTENT(in)  :: maxbudy,minproc,maxproc,myproc
   INTEGER,DIMENSION(0:maxbudy),INTENT(out) :: bud
   INTEGER,DIMENSION(minproc:maxproc),INTENT(in) :: aproc

!! * Local declarations
   INTEGER                            :: myproc_tmp,numbudy,i,j
   INTEGER,ALLOCATABLE,DIMENSION(:,:) :: budy_all

!!----------------------------------------------------------------------
!! * Executable part
!numbudy=maxproc-minproc+1
   IF (numbudy.EQ.1) THEN
     bud=myproc
     RETURN
   END IF

   ALLOCATE(budy_all(0:maxbudy,0:maxbudy))
   CALL calcbudy(budy_all,maxbudy+1)  !calcbudy calculate maxbudy pairs
!of exchanging process
!myproc_tmp defines the relative numbering of procs
   myproc_tmp=myproc-minproc-count(aproc(minproc:myproc-1)<0)
!   write(*,*)'makebudy',myiproc,myjproc,myproc_tmp
   bud(:)=budy_all(myproc_tmp,0:maxbudy)
   DEALLOCATE(budy_all)
   IF (MOD(numbudy,2).EQ.1) CALL calcbudy_odd(bud,maxbudy)  !if maxbudy is odd
   IF (minproc.NE.0) bud=bud+minproc  !adjust back so that zero becomes minproc
   if(count(aproc(:)<0).gt.0)then
           doi:do i=minproc+1,maxproc
           if(aproc(i).ne.IGROUNDPOINT)cycle doi
           doj:do j=0,maxbudy
           if(bud(j).lt.i)cycle doj
!   write(*,*)'makebudy budmodified',myiproc,i,j
           bud(j)=bud(j)+1
           enddo doj
           enddo doi
   endif
  END SUBROUTINE makebudy

!!======================================================================

  SUBROUTINE calcbudy_odd(bud,n)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE calcbudy_odd  ***
!&E
!&E ** Purpose : exchange the entry with maxbudy to -1
!&E
!&E ** Modified variables : bud
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER               ,INTENT(in)  :: n
   INTEGER,DIMENSION(0:n),INTENT(out) :: bud

!! * Local declarations
   INTEGER :: i,j

!!----------------------------------------------------------------------
!! * Executable part

   DO i=0,n
     j= bud(i)
     IF (j==n) THEN
       bud(i)=-1
     ELSE
       bud(i)=j
     END IF
   END DO

  END SUBROUTINE calcbudy_odd

!!======================================================================

  SUBROUTINE calcbudy(budy,n)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE calcbudy  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference : based on the idea published in
!&E                J. A. M. Schreuder, "Constructing Timetables for Sport Competitions,"
!&E                Mathematical Programming Study 13, pp. 58-67 (1980).
!&E ** History :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER                       ,INTENT(in)  :: n
   INTEGER,DIMENSION(0:n-1,0:n-1),INTENT(out) :: budy

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER                  :: i,j,t,a,b,k,n2
   LOGICAL                  :: test1,test2
   INTEGER,DIMENSION(n,n,n) :: x

!!----------------------------------------------------------------------
!! * Executable part

   IF (MOD(n,2).EQ.1) THEN
     CALL MPI_FINALIZE(ierr_mpi) ; STOP
   END IF

   IF (n==4) THEN
     budy(0,:)=(/ 3,2,1,0 /)
     budy(1,:)=(/ 2,3,0,1 /)
     budy(2,:)=(/ 1,0,3,2 /)
     budy(3,:)=(/ 0,1,2,3 /)
     RETURN
   ELSE IF (n==2) THEN
     budy(0,:)=(/ 1,0 /)
     budy(1,:)=(/ 0,1 /)
     RETURN
   END IF

   test1=.TRUE.
   t=0
   i=0
   x=0
   n2=n/2

10 t=t+1
   i=i+1
   IF (i.EQ.n2) test1=.NOT.test1
   IF (test1) THEN
     x(n,i,t)=1
     test1=.FALSE.
   ELSE
     x(i,n,t)=1
     test1=.TRUE.
   END IF
   test2=.TRUE.
   a=i
   b=i
   DO k=1,n2-1
     a= a-1
     a=MOD(a,n-1)
     b=MOD(b+1,n-1)
!this part is different from original paper
     IF (a.EQ.0)a=n-1
     IF (b.EQ.0)b=n-1
!untill here.

     IF (test2) THEN
       x(a,b,t)=1
       test2=.FALSE.
     ELSE
       x(b,a,t)=1
       test2=.TRUE.
     END IF
   END DO
   IF (t.LT.n-1) GOTO 10

!from now it is adjusted for mars budy.
   budy=0
   DO t=0,n-1
     DO j=1,n
       IF (t.EQ.0) THEN
         budy(j-1,0)=j-1
       ELSE
         DO i=1,n
           IF (x(i,j,t).EQ.1) THEN
             budy(i-1,t)=j-1
             budy(j-1,t)=i-1
           END IF
         END DO
       END IF
     END DO
   END DO

! do i=0,n-1
! write(*,*)"budyi=",i,budy(i,:)
! enddo

  END SUBROUTINE calcbudy

!!======================================================================

  SUBROUTINE TAG_PROC_IJ(tag)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE TAG_PROC_IJ  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : tag
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   CHARACTER(len=7),INTENT(out) :: tag

!!----------------------------------------------------------------------
!! * Executable part


   WRITE(tag, '(i3.3,".",i3.3)') myiproc, myjproc


  END SUBROUTINE TAG_PROC_IJ

!!======================================================================

  SUBROUTINE ADD_TAG(filename, ndim)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ADD_TAG  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : filename
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER, OPTIONAL,  INTENT(in)    :: ndim    ! deprecated, should be removed
   CHARACTER(len=*),   INTENT(inout) :: filename

!! * Local declarations
   CHARACTER(len=7)    :: tag
   INTEGER             :: len_fname, len_var, len_tag

!!----------------------------------------------------------------------
!! * Executable part

   CALL TAG_PROC_IJ(tag)

   len_fname = LEN_TRIM(filename)
   len_var   = LEN(filename)
   len_tag   = LEN_TRIM(tag) + 1

   IF ( (len_fname+len_tag) > len_var) THEN
      WRITE(impierror,*) "WARNING: filename '",TRIM(filename),"' truncated"
      filename = filename(1:len_var-len_tag)//'.'//tag
   ELSE
      filename=TRIM(filename)//'.'//tag
   END IF

  END SUBROUTINE ADD_TAG

!!======================================================================

  SUBROUTINE makeprocmap

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE makeprocmap  ***
!&E
!&E ** Purpose :
!&E           to verify if mpi.txt includes all thepoints between
! ig:id, and jb;jh (skips blank in the middle)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters,ONLY : jmin,jmax,imin,imax,valmanq
   use comvars2d, only : ig,id,jb,jh,h0_g
   IMPLICIT NONE

!! * Arguments
   INTEGER ,DIMENSION(imin:imax,jmin:jmax) :: procmap

!! * Local declarations
   INTEGER :: ip,jp,i,j,ierrr

!!----------------------------------------------------------------------
!! * Executable part

   DO jp=minjproc,maxjproc
   DO ip=miniproc,maxiproc
   procmap(LI(1,ip):LI(2,ip),LJ(1,jp):LJ(2,jp))=proc(ip,jp)

   enddo
   enddo
   ierrr=0
!   write(*,*) 'check if mpi.txt include all the wet point'
   do j=jmin,jmax
   do i=ig(j),id(j)
      if(procmap(i,j).eq.IGROUNDPOINT)then
        if((i.ge.imin+2.and.i.le.imax-2).and.(j.ge.jmin+2.and.j.le.jmax-2))then
! modif pgarreau  pour etre moins mere poule (tina dixit)...
!         if(count(h0(i-2:i+2,j-2:j+2)>-valmanq).ne.0)then
          if(h0_g(i,j) .gt. -valmanq)then
       write(*,*) 'stop: bad mpi.txt from igid i=',i,'j=',j
       ierrr=ierrr+1
       endif
       endif
      endif
   enddo
   enddo
   do i=imin,imax
   do j=jb(i),jh(i)
      if(procmap(i,j).eq.IGROUNDPOINT)then
        if((i.ge.imin+2.and.i.le.imax-2).and.(j.ge.jmin+2.and.j.le.jmax-2))then
! modif pgarreau  pour etre moins mere poule (tina dixit)...
!         if(count(h0(i-2:i+2,j-2:j+2)>-valmanq).ne.0)then
          if(h0_g(i,j) .gt. -valmanq)then
          write(*,*) 'stop: bad mpi.txt from jbjh i=',i,'j=',j
          ierrr=ierrr+1
         endif
        endif
      endif
   enddo
   enddo

   if(ierrr.eq.0)then
!   write(*,*) 'check mpi.txt OK'
   else
   write(*,*) 'check mpi.txt not OK,kill your code and edit your mpi.txt'
   stop
   endif
   return
  END SUBROUTINE makeprocmap

  SUBROUTINE countwetpoint

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE countwetpoint  ***
!&E
!&E ** Purpose :
!&E           to verify if mpi.txt  put almost same no of wet point for each
!nodes
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters,ONLY : jmin,jmax,imin,imax,valmanq
   use comvars2d, only : h0_g
   IMPLICIT NONE

!! * Arguments
   INTEGER ,DIMENSION(imin:imax,jmin:jmax) :: procmap
!   LOGICAL ,DIMENSION(imin-2:imax+2,jmin-2:jmax+2) :: logich0

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: ip,jp,i,j,iwet,ierror

!!----------------------------------------------------------------------
!! * Executable part

! Set to .true. over wet points
!   logich0(:,:)=.FALSE.
!   WHERE(h0 > -valmanq) logich0(imin:imax,jmin:jmax)=.TRUE.
!   write(*,*) 'count wet point for each mpi grid'
   ierror=0

! Loop on domain division
   DO jp=minjproc,maxjproc
     DO ip=miniproc,maxiproc

       iwet=0
! Loop in ljmin,ljmax and limin,limax of each subdomain
       DO j=LJ(1,jp),LJ(2,jp)
         DO i=LI(1,ip),LI(2,ip)
           IF (h0_g(i,j) > -valmanq)iwet=iwet+1
!cval test too large because of +2 -2 if(count( logich0(i-2:i+2,j-2:j+2)).ne.0)iwet=iwet+1
         END DO
       END DO
!       write(*,*) 'ip=',ip,'jp=',jp,'wetpoint=',iwet

! Make sure there is none wet cell (iwet/=0) for sleeping CPU (IGROUNDPOINT=.true.)
       IF (iwet /= 0 .AND. proc(ip,jp) == IGROUNDPOINT) THEN
         WRITE(*,*) 'stop:ip=',ip,'jp=',jp,'You should not have any wet cell in this CPU'
         ierror=ierror+1
       END IF

     END DO
   END DO
   IF (ierror /= 0) THEN
     CALL MPI_FINALIZE(ierr_mpi) ; STOP
   END IF

  END SUBROUTINE countwetpoint

  SUBROUTINE shiftforblank(LA,minaproc,maxaproc,APROC)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE shiftforblank  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : LA
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER           ,INTENT(in)  :: maxaproc,minaproc
   INTEGER ,DIMENSION(4,minaproc:maxaproc) :: LA
   INTEGER ,DIMENSION(minaproc:maxaproc),INTENT(in) ::aproc

!! * Local declarations
   INTEGER :: i

!!----------------------------------------------------------------------
!! * Executable part

   DO i=minaproc,maxaproc
     if(aproc(i).eq.IGROUNDPOINT) then
          LA(:,i+1:maxaproc)=LA(:,i:maxaproc-1)
          LA(:,i)=IGROUNDPOINT
     endif
   enddo

  END SUBROUTINE shiftforblank

!!======================================================================
  SUBROUTINE divide_simple1d(amin,amax,minaproc,maxdimaproc,maxaproc,LA)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE divide_simple1d  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : LA
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                    INTENT(in)  :: amin,amax
   INTEGER,                                    INTENT(in)  :: maxaproc,minaproc,maxdimaproc
   INTEGER, DIMENSION(4,minaproc:maxdimaproc), INTENT(out) :: LA

!! * Local declarations
   INTEGER :: n,m,i,numaproc,ndf

!!----------------------------------------------------------------------
!! * Executable part

   numaproc=maxaproc-minaproc+1
   LA(1,minaproc)=amin
   n=(amax-amin  +1)/ numaproc
   m=MOD((amax-amin+1) ,numaproc )
   ndf=0
   DO i=minaproc,maxaproc
     LA(2,i)=LA(1,i)+n-1
     if(i.eq.maxaproc)then
       if(LA(2,maxaproc).ne.amax) then
           write(*,*)'stop: error maxaproc la2 ne amax'
           write(*,*)maxaproc,amax,iproc
           stop
       else
               cycle
       endif
     endif
   if(m.gt.ndf) LA(2,i)=LA(2,i)+1
     LA(1,i+1)=LA(2,i)+1
     ndf=ndf+1
   enddo

   LA(3,:)=LA(2,:)-LA(1,:)+1
   LA(4,:)=LA(1,:)-LA(1,minaproc)

  END SUBROUTINE divide_simple1d

!!======================================================================

  SUBROUTINE divide_withsum1d(amin,amax,minaproc,maxdimaproc,maxaproc,asum,LA)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE divide_withsum1d  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : LA
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                    INTENT(in)  :: amin,amax
   INTEGER,                                    INTENT(in)  :: maxaproc,minaproc,maxdimaproc
   INTEGER, DIMENSION(amin:amax),              INTENT(in)  :: aSUM
   INTEGER, DIMENSION(4,minaproc:maxdimaproc), INTENT(out) :: LA

!! * Local declarations
   INTEGER ::n,y,x,j,numaproc,iiout

!!----------------------------------------------------------------------
!! * Executable part

   LA(1,minaproc)=amin
   numaproc=maxaproc-minaproc+1
   if(numaproc.eq.1)goto 101

   n=SUM(aSUM(amin:amax))/numaproc
   y=minaproc
   x=0
   DO j=amin,amax
     x=x+aSUM(j)
     IF (x.Ge.n) THEN
       La(2,Y)=j
       y=y+1
       La(1,Y)=j+1
       IF (y.EQ.maxaproc) GOTO 101
       x=0
       n=SUM(aSUM(j+1:amax))/(numaproc-(y-minaproc))
     END IF
   END DO
101 LA(2,maxaproc)=amax
   LA(3,:)=LA(2,:)-LA(1,:)+1
   LA(4,:)=LA(1,:)-LA(1,minaproc)

  END SUBROUTINE divide_withsum1d

!!======================================================================

  SUBROUTINE get_jsum(jsum)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE get_jsum  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : jsum
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : jmin,jmax
   USE comvars2d,  ONLY : ig,id
   IMPLICIT NONE

!! * Arguments
   INTEGER, DIMENSION(jmin:jmax), INTENT(out) :: jsum

!!----------------------------------------------------------------------
!! * Executable part
   jsum(:) = MAX(0, id(:)-ig(:)+1)

  END SUBROUTINE get_jsum

!!======================================================================

  SUBROUTINE get_isum(isum)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE get_isum  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : isum
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : imin,imax
   USE comvars2d,  ONLY : jb,jh
   IMPLICIT NONE

!! * Arguments
   INTEGER, DIMENSION(imin:imax), INTENT(out) :: isum

!!----------------------------------------------------------------------
!! * Executable part
   isum(:) = MAX(0, jh(:)-jb(:)+1)

  END SUBROUTINE get_isum

!!======================================================================

  SUBROUTINE sendrecv_rlg(dest,size_give,give,size_get,get)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sendrecv_rlg  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rlg
   IMPLICIT NONE

!! * Arguments
   INTEGER                            ,INTENT(inout) :: dest,size_give,size_get
   REAL(kind=rlg),DIMENSION(size_give),INTENT(inout) :: GIVE
   REAL(kind=rlg),DIMENSION(size_get) ,INTENT(inout) :: GET

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: send_tag,recv_tag,NDIM_SEND,NDIM_RECV
   INTEGER :: stat(MPI_STATUS_SIZE)
   INTEGER :: TYPE_NDIM_SEND,TYPE_NDIM_RECV

!!----------------------------------------------------------------------
!! * Executable part

   NDIM_SEND=1
   NDIM_RECV=1
   CALL MPI_TYPE_CONTIGUOUS(NDIM_SEND,type_mpi_rlg,TYPE_NDIM_SEND, ierr_mpi )
   CALL MPI_TYPE_CONTIGUOUS(NDIM_RECV,type_mpi_rlg,TYPE_NDIM_RECV, ierr_mpi )
   CALL MPI_TYPE_COMMIT(TYPE_NDIM_SEND, ierr_mpi )
   CALL MPI_TYPE_COMMIT(TYPE_NDIM_RECV, ierr_mpi )

   send_tag=300*iproc+dest
   recv_tag=300*dest+iproc
   CALL MPI_SENDRECV(GIVE,SIZE_GIVE,TYPE_NDIM_SEND,dest,send_tag &
     &,GET,SIZE_GET,TYPE_NDIM_RECV,dest,recv_tag &
     &,MPI_COMM_MARS,stat,ierr_mpi)
   CALL MPI_TYPE_FREE(TYPE_NDIM_SEND, ierr_mpi )
   CALL MPI_TYPE_FREE(TYPE_NDIM_RECV, ierr_mpi )
   RETURN

  END SUBROUTINE sendrecv_rlg

!!======================================================================

  SUBROUTINE sendrecv_one_int(dest, give, get)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sendrecv_one_int  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   INTEGER, INTENT(in)  :: dest
   INTEGER, INTENT(in)  :: give
   INTEGER, INTENT(out) :: get

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: send_tag, recv_tag

!!----------------------------------------------------------------------
!! * Executable part

   send_tag = 300*iproc+dest
   recv_tag = 300*dest+iproc

   CALL MPI_SENDRECV(give, 1, MPI_INTEGER, dest, send_tag, &
                     get,  1, MPI_INTEGER, dest, recv_tag, &
                     MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)

  END SUBROUTINE sendrecv_one_int

!!======================================================================

  SUBROUTINE sendrecv_rsh(dest,size_give,give,size_get,get)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sendrecv_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh
   IMPLICIT NONE

!! * Arguments
   INTEGER                            ,INTENT(inout) :: dest,size_give,size_get
   REAL(kind=rsh),DIMENSION(size_give),INTENT(inout) :: GIVE
   REAL(kind=rsh),DIMENSION(size_get) ,INTENT(inout) :: GET

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: send_tag,recv_tag,NDIM_SEND,NDIM_RECV
   INTEGER :: stat(MPI_STATUS_SIZE)
   INTEGER :: TYPE_NDIM_SEND,TYPE_NDIM_RECV

!!----------------------------------------------------------------------
!! * Executable part

   NDIM_SEND=1
   NDIM_RECV=1
   CALL MPI_TYPE_CONTIGUOUS(NDIM_SEND,type_mpi_rsh,TYPE_NDIM_SEND, ierr_mpi )
   CALL MPI_TYPE_CONTIGUOUS(NDIM_RECV,type_mpi_rsh,TYPE_NDIM_RECV, ierr_mpi )
   CALL MPI_TYPE_COMMIT(TYPE_NDIM_SEND, ierr_mpi )
   CALL MPI_TYPE_COMMIT(TYPE_NDIM_RECV, ierr_mpi )

   send_tag=300*iproc+dest
   recv_tag=300*dest+iproc
   CALL MPI_SENDRECV(GIVE,SIZE_GIVE,TYPE_NDIM_SEND,dest,send_tag &
     &,GET,SIZE_GET,TYPE_NDIM_RECV,dest,recv_tag &
     &,MPI_COMM_MARS,stat,ierr_mpi)
    IF (ierr_mpi.ne.0)  WRITE(impierror,*)'you have error in sendrecv_rsh',dest,ierr_mpi,stat
   CALL MPI_TYPE_FREE(TYPE_NDIM_SEND, ierr_mpi )
   CALL MPI_TYPE_FREE(TYPE_NDIM_RECV, ierr_mpi )
   RETURN

  END SUBROUTINE sendrecv_rsh

!!======================================================================



!!======================================================================

  SUBROUTINE ex_j_rsh(x,y,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_j_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :  x is to exchange down  (for j-1 j-2 ,,,
!&E                   y is to exchange up    (for j+1 j+2 ..
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,ljmin,ljmax
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                     INTENT(in)    :: x, y, ndim
   INTEGER,                                     INTENT(in)    :: i0, i1, j0, j1
   REAL(kind=rsh), DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi, isize
   REAL(kind=rsh), DIMENSION(ndim,i0:i1,-x) :: dr,ds
   REAL(kind=rsh), DIMENSION(ndim,i0:i1, y) :: ur,us

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   
   

!down x(m)   x < 0
   IF ( x /= 0 ) THEN
     isize = size(ds)
     IF (myj_uproc /= MPI_PROC_NULL) ds = A(:,:,ljmax+x+1:ljmax)
     CALL MPI_SENDRECV( &
        ds, isize, type_mpi_rsh, myj_uproc, ljmax,   &
        dr, isize, type_mpi_rsh, myj_dproc, ljmin-1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi )
     IF (myj_dproc /= MPI_PROC_NULL) A(:,:,ljmin+x:ljmin-1) = dr
   END IF

! up y(p)    y > 0
   IF ( y /= 0 ) THEN
     isize = size(us)
     IF (myj_dproc /= MPI_PROC_NULL) us = A(:,:,ljmin:ljmin+y-1)
     CALL MPI_SENDRECV( &
        us, isize, type_mpi_rsh, myj_dproc, ljmin,   &
        ur, isize, type_mpi_rsh, myj_uproc, ljmax+1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)
     IF (myj_uproc /= MPI_PROC_NULL) A(:,:,ljmax+1:ljmax+y) = ur
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_j_rsh

!!======================================================================

  SUBROUTINE ex_i_rsh(x,y,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_i_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :  x is to exchange down  (for i-1 i-2 ,,,
!&E                   y is to exchange up    (for i+1 i+2 ..
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh, limin, limax
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                     INTENT(in)    :: x, y, ndim
   INTEGER,                                     INTENT(in)    :: i0, i1, j0, j1
   REAL(kind=rsh), DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi, jsize
   REAL(kind=rsh), DIMENSION(ndim,-x,j0:j1) :: dr,ds
   REAL(kind=rsh), DIMENSION(ndim, y,j0:j1) :: ur,us

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   
   

!down x(m)   x < 0
   IF ( x /= 0 ) THEN
     jsize = size(ds)
     IF (myi_uproc /= MPI_PROC_NULL) ds = A(:,limax+x+1:limax,:)
     CALL MPI_SENDRECV( &
        ds, jsize, type_mpi_rsh, myi_uproc, limax,   &
        dr, jsize, type_mpi_rsh, myi_dproc, limin-1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi )
     IF (myi_dproc /= MPI_PROC_NULL) A(:,limin+x:limin-1,:) = dr
   END IF

! up y(p)    y > 0
   IF ( y /= 0 ) THEN
     jsize = size(us)
     IF (myi_dproc /= MPI_PROC_NULL) us = A(:,limin:limin+y-1,:)
     CALL MPI_SENDRECV( &
        us, jsize, type_mpi_rsh, myi_dproc, limin,   &
        ur, jsize, type_mpi_rsh, myi_uproc, limax+1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)
     IF (myi_uproc /= MPI_PROC_NULL) A(:,limax+1:limax+y,:) = ur
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_i_rsh

!!======================================================================

  SUBROUTINE ex_j_rlg(x,y,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_j_rlg  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :  x is to exchange down  (for j-1 j-2 ,,,
!&E                   y is to exchange up    (for j+1 j+2 ..
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rlg,ljmin,ljmax
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                     INTENT(in)    :: x, y, ndim
   INTEGER,                                     INTENT(in)    :: i0, i1, j0, j1
   REAL(kind=rlg), DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi, isize
   REAL(kind=rlg), DIMENSION(ndim,i0:i1,-x) :: dr,ds
   REAL(kind=rlg), DIMENSION(ndim,i0:i1, y) :: ur,us

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   
   

!down x(m)   x < 0
   IF ( x /= 0 ) THEN
     isize = size(ds)
     IF (myj_uproc /= MPI_PROC_NULL) ds = A(:,:,ljmax+x+1:ljmax)
     CALL MPI_SENDRECV( &
        ds, isize, type_mpi_rlg, myj_uproc, ljmax,   &
        dr, isize, type_mpi_rlg, myj_dproc, ljmin-1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi )
     IF (myj_dproc /= MPI_PROC_NULL) A(:,:,ljmin+x:ljmin-1) = dr
   END IF

! up y(p)    y > 0
   IF ( y /= 0 ) THEN
     isize = size(us)
     IF (myj_dproc /= MPI_PROC_NULL) us = A(:,:,ljmin:ljmin+y-1)
     CALL MPI_SENDRECV( &
        us, isize, type_mpi_rlg, myj_dproc, ljmin,   &
        ur, isize, type_mpi_rlg, myj_uproc, ljmax+1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)
     IF (myj_uproc /= MPI_PROC_NULL) A(:,:,ljmax+1:ljmax+y) = ur
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_j_rlg

!!======================================================================

  SUBROUTINE ex_i_rlg(x,y,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_i_rlg  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :  x is to exchange down  (for i-1 i-2 ,,,
!&E                   y is to exchange up    (for i+1 i+2 ..
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rlg, limin, limax
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                     INTENT(in)    :: x, y, ndim
   INTEGER,                                     INTENT(in)    :: i0, i1, j0, j1
   REAL(kind=rlg), DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi, jsize
   REAL(kind=rlg), DIMENSION(ndim,-x,j0:j1) :: dr,ds
   REAL(kind=rlg), DIMENSION(ndim, y,j0:j1) :: ur,us

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   
   

!down x(m)   x < 0
   IF ( x /= 0 ) THEN
     jsize = size(ds)
     IF (myi_uproc /= MPI_PROC_NULL) ds = A(:,limax+x+1:limax,:)
     CALL MPI_SENDRECV( &
        ds, jsize, type_mpi_rlg, myi_uproc, limax,   &
        dr, jsize, type_mpi_rlg, myi_dproc, limin-1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi )
     IF (myi_dproc /= MPI_PROC_NULL) A(:,limin+x:limin-1,:) = dr
   END IF

! up y(p)    y > 0
   IF ( y /= 0 ) THEN
     jsize = size(us)
     IF (myi_dproc /= MPI_PROC_NULL) us = A(:,limin:limin+y-1,:)
     CALL MPI_SENDRECV( &
        us, jsize, type_mpi_rlg, myi_dproc, limin,   &
        ur, jsize, type_mpi_rlg, myi_uproc, limax+1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)
     IF (myi_uproc /= MPI_PROC_NULL) A(:,limax+1:limax+y,:) = ur
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_i_rlg

!!======================================================================

  SUBROUTINE ex_j_int(x,y,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_j_int  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :  x is to exchange down  (for j-1 j-2 ,,,
!&E                   y is to exchange up    (for j+1 j+2 ..
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : ljmin,ljmax
   IMPLICIT NONE

!! * Arguments
   INTEGER,                              INTENT(in)    :: x, y, ndim
   INTEGER,                              INTENT(in)    :: i0, i1, j0, j1
   INTEGER, DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi, isize
   INTEGER, DIMENSION(ndim,i0:i1,-x) :: dr,ds
   INTEGER, DIMENSION(ndim,i0:i1, y) :: ur,us

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   
   

!down x(m)   x < 0
   IF ( x /= 0 ) THEN
     isize = size(ds)
     IF (myj_uproc /= MPI_PROC_NULL) ds = A(:,:,ljmax+x+1:ljmax)
     CALL MPI_SENDRECV( &
        ds, isize, MPI_INTEGER, myj_uproc, ljmax,   &
        dr, isize, MPI_INTEGER, myj_dproc, ljmin-1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi )
     IF (myj_dproc /= MPI_PROC_NULL) A(:,:,ljmin+x:ljmin-1) = dr
   END IF

! up y(p)    y > 0
   IF ( y /= 0 ) THEN
     isize = size(us)
     IF (myj_dproc /= MPI_PROC_NULL) us = A(:,:,ljmin:ljmin+y-1)
     CALL MPI_SENDRECV( &
        us, isize, MPI_INTEGER, myj_dproc, ljmin,   &
        ur, isize, MPI_INTEGER, myj_uproc, ljmax+1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)
     IF (myj_uproc /= MPI_PROC_NULL) A(:,:,ljmax+1:ljmax+y) = ur
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_j_int

!!======================================================================

  SUBROUTINE ex_i_int(x,y,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_i_int  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :  x is to exchange down  (for i-1 i-2 ,,,
!&E                   y is to exchange up    (for i+1 i+2 ..
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : limin, limax
   IMPLICIT NONE

!! * Arguments
   INTEGER,                              INTENT(in)    :: x, y, ndim
   INTEGER,                              INTENT(in)    :: i0, i1, j0, j1
   INTEGER, DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi, jsize
   INTEGER, DIMENSION(ndim,-x,j0:j1) :: dr,ds
   INTEGER, DIMENSION(ndim, y,j0:j1) :: ur,us

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   
   

!down x(m)   x < 0
   IF ( x /= 0 ) THEN
     jsize = size(ds)
     IF (myi_uproc /= MPI_PROC_NULL) ds = A(:,limax+x+1:limax,:)
     CALL MPI_SENDRECV( &
        ds, jsize, MPI_INTEGER, myi_uproc, limax,   &
        dr, jsize, MPI_INTEGER, myi_dproc, limin-1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi )
     IF (myi_dproc /= MPI_PROC_NULL) A(:,limin+x:limin-1,:) = dr
   END IF

! up y(p)    y > 0
   IF ( y /= 0 ) THEN
     jsize = size(us)
     IF (myi_dproc /= MPI_PROC_NULL) us = A(:,limin:limin+y-1,:)
     CALL MPI_SENDRECV( &
        us, jsize, MPI_INTEGER, myi_dproc, limin,   &
        ur, jsize, MPI_INTEGER, myi_uproc, limax+1, &
        MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)
     IF (myi_uproc /= MPI_PROC_NULL) A(:,limax+1:limax+y,:) = ur
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_i_int

!!======================================================================

  SUBROUTINE ex_peri_y_rsh(xx,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_peri_y_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description : x is to exchange 1 or 2
!&E **             : j0==ljmin, j1==ljmax
!&E **             : i0 and i1 is equal to the limit where each node holds data.
!&E **             : (take it a bit bigger just in case)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,ljmin,ljmax,jmin,jmax,limin,limax
   USE comvars2d,  ONLY : jb,jh
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                     INTENT(in)    :: xx,ndim
   INTEGER,                                     INTENT(in)    :: i0,i1,j0,j1
   REAL(kind=rsh), DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi
   INTEGER  :: x,y,ds_size,dr_size,stag,rtag,uproc,jn,js
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: dr,ds

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   jn = jh((limax-limin)/2)
   js = jb((limax-limin)/2)

   

! interior of the domain, no OBC
   IF ( (ljmin /= jmin) .AND. (ljmax /= jmax) ) RETURN
! MPI division along vertical stripes
   IF ( (ljmin == jmin) .AND. (ljmax == jmax) ) THEN
! No MPI splicing in j
     IF (xx >= 2) A(:,:,js-1) = A(:,:,jn-2)
     A(:,:,js) = A(:,:,jn-1)
     A(:,:,jn) = A(:,:,js+1)
     IF (xx == 3) THEN
       A(:,:,js-2) = A(:,:,jn-3)
       A(:,:,jn+1) = A(:,:,js+2)
       A(:,:,jn+2) = A(:,:,js+3)
     END IF
     RETURN
   END IF

   IF ( myjproc == minjproc ) THEN
     uproc = PROC(myiproc,maxjproc)
     stag = 0+1000*myiproc
     rtag = 1+1000*myiproc
     x = xx
     y = 1
     IF (xx == 3) y = 3
     ALLOCATE(dr(ndim,i0:i1,x))
     ALLOCATE(ds(ndim,i0:i1,y))
! copy into send buffer
     ds(:,:,1) = A(:,:,js+1)
     IF (xx == 3) ds(:,:,1:3) = A(:,:,js+1:js+3)
   ELSE IF ( myjproc == maxjproc ) THEN
     uproc=PROC(myiproc,minjproc)
     stag = 1+1000*myiproc
     rtag = 0+1000*myiproc
     x = 1
     y = xx
     IF (xx == 3) x=3
     ALLOCATE(dr(ndim,i0:i1,x))
     ALLOCATE(ds(ndim,i0:i1,y))
     IF (xx == 3) THEN
       ds(:,:,1:3) = A(:,:,jn-3:jn-1)
     ELSE IF (xx == 2) THEN
       ds(:,:,1:2) = A(:,:,jn-2:jn-1)
     ELSE IF (xx == 1) THEN
       ds(:,:,1  ) = A(:,:,jn-1)
     END IF
   END IF

   ds_size = size(ds)
   dr_size = size(dr)

   CALL MPI_SENDRECV(                         &
     ds, ds_size, type_mpi_rsh, uproc, stag,  &
     dr, dr_size, type_mpi_rsh, uproc, rtag,  &
     MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)

   IF ( myjproc == minjproc ) THEN
     IF (xx == 1) A(:,:,js     ) = dr(:,:,1)
     IF (xx == 2) A(:,:,js-1:js) = dr(:,:,1:2)
     IF (xx == 3) A(:,:,js-2:js) = dr(:,:,1:3)
   ELSE IF ( myjproc == maxjproc ) THEN
     A(:,:,jn) = dr(:,:,1)
     IF (xx == 3) A(:,:,jn:jn+2) = dr(:,:,1:3)
   END IF
   DEALLOCATE(ds, dr)

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_peri_y_rsh

  SUBROUTINE ex_peri_y_avcv(i0,i1,j0,j1,A,C)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_peri_y_avcv  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description : copy av(jmax-1) into av(jmin+1)
!&E **             : and cv(jmin+1) into cv(jmax-1)
!&E **             : j0==ljmin, j1==ljmax
!&E **             : i0 and i1 is equal to the limit where each node holds data.
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rlg,ljmin,ljmax,jmin,jmax,limin,limax
   USE comvars2d,  ONLY : jb,jh
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                INTENT(in)    :: i0,i1,j0,j1
   REAL(kind=rlg), DIMENSION(i0:i1,j0:j1), INTENT(inout) :: A,C

!! * Local declarations
   INTEGER  :: ierr_mpi
   INTEGER  :: ds_size,dr_size,stag,rtag,uproc,jn,js
   REAL(kind=rlg), DIMENSION(i0:i1) :: dr,ds

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   jn = jh((limax-limin)/2)
   js = jb((limax-limin)/2)

! interior of the domain, no OBC
   IF ( (ljmin /= jmin) .AND. (ljmax /= jmax) ) RETURN
! MPI division along vertical stripes
   IF ( (ljmin == jmin) .AND. (ljmax == jmax) ) THEN
     A(:,js  ) = A(:,jn-1)
     C(:,jn-1) = C(:,js)
     RETURN
   END IF

   IF ( myjproc == minjproc ) THEN
     uproc=PROC(myiproc,maxjproc)
     stag = 0+1000*myiproc
     rtag = 1+1000*myiproc
! copy into send buffer
     ds(:) = C(:,js)
   ELSE IF ( myjproc == maxjproc ) THEN
     uproc=PROC(myiproc,minjproc)
     stag = 1+1000*myiproc
     rtag = 0+1000*myiproc
     ds(:) = A(:,jn-1)
   END IF

   ds_size = size(ds)
   dr_size = size(dr)

   CALL MPI_SENDRECV(                        &
     ds, ds_size, type_mpi_rlg, uproc, stag, &
     dr, dr_size, type_mpi_rlg, uproc, rtag, &
     MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)

   IF ( myjproc == minjproc ) THEN
     A(:,js  ) = dr(:)
   ELSE IF ( myjproc == maxjproc ) THEN
     C(:,jn-1) = dr(:)
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_peri_y_avcv

!!======================================================================

  SUBROUTINE ex_peri_x_rsh(xx,ndim,i0,i1,j0,j1,A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_peri_x_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description : xx is to exchange 1 or 2
!&E **             : i0==limin, i1==limax
!&E **             : j0 and j1 is equal to the limit where each node holds data.
!&E **             : (take it a bit bigger just in case)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,imin,imax,limin,limax,ljmin,ljmax
   USE comvars2d,  ONLY : ig,id
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                     INTENT(in)    :: xx,ndim
   INTEGER,                                     INTENT(in)    :: i0,i1,j0,j1
   REAL(kind=rsh), DIMENSION(ndim,i0:i1,j0:j1), INTENT(inout) :: A

!! * Local declarations
   INTEGER  :: ierr_mpi
   INTEGER  :: x,y,ds_size,dr_size,stag,rtag,uproc,ie,iw
   REAL(kind=rsh), ALLOCATABLE, DIMENSION(:,:,:) :: dr,ds

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   ie = id((ljmax-ljmin)/2)
   iw = ig((ljmax-ljmin)/2)

   

! interior of the domain, no OBC
   IF ( (limin /= imin) .AND. (limax /= imax) ) RETURN
! MPI division along horizontal stripes
   IF ( (limin == imin) .AND. (limax == imax) ) THEN
! No MPI splicing in i
     IF (xx >= 2) A(:,iw-1,:) = A(:,ie-2,:)
     A(:,iw,:) = A(:,ie-1,:)
     A(:,ie,:) = A(:,iw+1,:)
     IF (xx == 3) THEN
       A(:,iw-2,:) = A(:,ie-3,:)
       A(:,ie+1,:) = A(:,iw+2,:)
       A(:,ie+2,:) = A(:,iw+3,:)
     END IF
     RETURN
   END IF

   IF ( myiproc == miniproc ) THEN
     uproc = PROC(maxiproc,myjproc)
     stag = 0+1000*myjproc
     rtag = 1+1000*myjproc
     x = xx
     y = 1
     IF (xx==3) y = 3
     ALLOCATE(dr(ndim,x,j0:j1))
     ALLOCATE(ds(ndim,y,j0:j1))
! copy into send buffer
     ds(:,1,:)=A(:,iw+1,:)
     IF (xx == 3) ds(:,1:3,:) = A(:,iw+1:iw+3,:)
   ELSE IF ( myiproc == maxiproc ) THEN
     uproc=PROC(miniproc,myjproc)
     stag = 1+1000*myjproc
     rtag = 0+1000*myjproc
     x=1
     y=xx
     IF (xx == 3) x=3
     ALLOCATE(dr(ndim,x,j0:j1))
     ALLOCATE(ds(ndim,y,j0:j1))
     IF (xx == 3) THEN
       ds(:,1:3,:) = A(:,ie-3:ie-1,:)
     ELSE IF (xx == 2) THEN
       ds(:,1:2,:) = A(:,ie-2:ie-1,:)
     ELSE IF (xx == 1) THEN
       ds(:,1  ,:) = A(:,ie-1,:)
     END IF
   END IF

   ds_size = size(ds)
   dr_size = size(dr)

! exchange and deallocation of the initial address of send buffer (DS)
   CALL MPI_SENDRECV(                         &
     ds, ds_size, type_mpi_rsh, uproc, stag,  &
     dr, dr_size, type_mpi_rsh, uproc, rtag,  &
     MPI_COMM_MARS , MPI_STATUS_IGNORE, ierr_mpi)

! copy of the initial address of the receive buffer (dr) into the common variable (A)
   IF ( myiproc == miniproc ) THEN
     IF (xx == 1) A(:,iw,     :) = dr(:,1,:)
     IF (xx == 2) A(:,iw-1:iw,:) = dr(:,1:2,:)
     IF (xx == 3) A(:,iw-2:iw,:) = dr(:,1:3,:)
   ELSE IF ( myiproc == maxiproc ) THEN
     A(:,ie,:) = dr(:,1,:)
     IF (xx == 3) A(:,ie:ie+2,:) = dr(:,1:3,:)
   END IF
   DEALLOCATE(ds, dr)

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_peri_x_rsh

!!======================================================================

  SUBROUTINE ex_peri_x_avcv(i0,i1,j0,j1,A,C)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ex_peri_x_avcv  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description : copy av(imax-1) into av(imin+1)
!&E **             : and cv(imin+1) into cv(imax-1)
!&E **             : i0==limin, i1==limax
!&E **             : j0 and j1 is equal to the limit where each node holds data.
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rlg,imin,imax,limin,limax,ljmin,ljmax
   USE comvars2d,  ONLY : ig,id
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                INTENT(in)    :: i0,i1,j0,j1
   REAL(kind=rlg), DIMENSION(i0:i1,j0:j1), INTENT(inout) :: A,C

!! * Local declarations
   INTEGER  :: ierr_mpi
   INTEGER  :: ds_size,dr_size,stag,rtag,uproc,iw,ie
   REAL(kind=rlg), DIMENSION(j0:j1) :: dr,ds

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   ie = id((ljmax-ljmin)/2)
   iw = ig((ljmax-ljmin)/2)

! interior of the domain, no OBC
   IF ( (limin /= imin) .AND. (limax /= imax) ) RETURN
! MPI division along horizontal stripes
   IF ( (limin == imin) .AND. (limax == imax) ) THEN
     A(iw,  :) = A(ie-1,:)
     C(ie-1,:) = C(iw,  :)
     RETURN
   END IF

   IF ( myiproc == miniproc ) THEN
     uproc=PROC(maxiproc,myjproc)
     stag = 0+1000*myjproc
     rtag = 1+1000*myjproc
! copy into send buffer
     ds(:) = C(iw,:)
   ELSE IF ( myiproc == maxiproc ) THEN
     uproc=PROC(miniproc,myjproc)
     stag = 1+1000*myjproc
     rtag = 0+1000*myjproc
     ds(:) = A(ie-1,:)
   END IF

   ds_size = size(ds)
   dr_size = size(dr)

! exchange of the initial address of send buffer (ds)
   CALL MPI_SENDRECV(                        &
     ds, ds_size, type_mpi_rlg, uproc, stag, &
     dr, dr_size, type_mpi_rlg, uproc, rtag, &
     MPI_COMM_MARS, MPI_STATUS_IGNORE, ierr_mpi)

! copy of the initial address of the receive buffer (dr) into the common variables (A and C)
   IF ( myiproc == miniproc ) THEN
     A(iw,  :) = dr(:)
   ELSE IF ( myiproc == maxiproc ) THEN
     C(ie-1,:) = dr(:)
   END IF

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

  END SUBROUTINE ex_peri_x_avcv

!!======================================================================

  SUBROUTINE correct_rsh(A,B,NDIM)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE correct_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : A (, B ?)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,ljmin,ljmax,limin,limax
   USE parameters, ONLY : imin,imax,jmin,jmax
   IMPLICIT NONE

!! * Arguments
   INTEGER                                               ,INTENT(in)    :: ndim
   REAL(kind=rsh),DIMENSION(NDIM,imin:imax,jmin:jmax)    ,INTENT(inout) :: A
   REAL(kind=rsh),DIMENSION(NDIM,limin:limax,ljmin:ljmax),INTENT(in)    :: B
!REAL(kind=rsh),DIMENSION(NDIM,limin:limax,ljmin:ljmax)               :: F

!! * Local declarations
   INTEGER                                     :: ierr_mpi
   INTEGER                                     :: x
   INTEGER                                     :: TYPE_NDIM
   INTEGER                                     :: count_send,count_get,dest,etigo,n
   INTEGER                                     :: stat(MPI_STATUS_SIZE)
   INTEGER,DIMENSION(nprocs)                   :: recvcounts,displs
   INTEGER                                     :: nsize,i,j,nij,sendcount,iroot
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: C,D
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:) :: E
   REAL(kind=rsh)                              :: toto

!!----------------------------------------------------------------------
!! * Executable part


!WRITE(impierror,*)   "i am in correct_rsh"
   IF (nbblank.NE.0) THEN
!  WRITE(impierror,*)   "i am in nbblanc.ne0 correct_rsh"

     CALL MPI_TYPE_CONTIGUOUS(NDIM,type_mpi_rsh,TYPE_NDIM, ierr_mpi )
     CALL MPI_TYPE_COMMIT(TYPE_NDIM, ierr_mpi )
     nsize=SUM(proc_size(minp:maxp))
!tinatestWRITE(impierror,*)   "nsize",nsize
     ALLOCATE(C(NDIM,LI(3,myiproc),LJ(3,myjproc)))

!  DO j=ljmin,ljmax
!     DO i=limin,limax
!      F(1,i,j)=dble(1000*i+j)
!    END DO
!  ENDDO
!C(:,:,:)=F(:,limin:limax,ljmin:ljmax)  !data that i hold, which will be send to everyone
     C(:,:,:)=B(:,limin:limax,ljmin:ljmax)  !data that i hold, which will be send to everyone
     ALLOCATE(E(NDIM,nsize))                ! data box that i keep all the data from everyone
     sendcount=proc_size(iproc)
     recvcounts(1:nprocs)=proc_size(minp:maxp)
     displs(1:nprocs)=proc_displs(minp:maxp)
!WRITE(impierror,*)   "just before gatherv"
!WRITE(impierror,*)   "jusmaster",master,iproc
     iroot=0 !here 0 means master
     CALL MPI_GATHERV(C,sendcount , TYPE_NDIM, &
       & E, recvcounts, displs,     TYPE_NDIM, iroot, &
       &MPI_COMM_MARS, ierr_mpi)
!WRITE(impierror,*)   "just after gatherv",ierr_mpi
     IF (ierr_mpi.NE.0) WRITE(impierror,*)              ierr_mpi,"ierror at exchange"
!WRITE(impierror,*)   "just before mpi_barrier"
!call MPI_BARRIER(MPI_COMM_MARS,ierr_mpi)
!WRITE(impierror,*)   "just after mpi_barrier"
     CALL MPI_TYPE_FREE(TYPE_NDIM, ierr_mpi )
     DEALLOCATE(c)
     nij=0
     DO n=minp,maxp
!IF (n==iproc)CYCLE
!nij=proc_displs(N)
! here i can also use function reshape
       DO j=proc_ljmin(n),proc_ljmax(n)
         DO i=proc_limin(n),proc_limax(n)
           nij=nij+1
           A(:,i,j)=E(:,nij)
!if(A(1,i,j).ne.dble(1000*i+j))then
!  WRITE(impierror,*)   "a.ne.dble",a(1,i,j),i,j
!  if(master)  then
!      WRITE(impierror,*)   "master is a.ne.dble stop"
!   CALL MPI_FINALIZE(ierr_mpi) ; STOP
!  endif
!else
!   WRITE(impierror,*)   "a.eq.dble",a(1,i,j),i,j
!END IF
         END DO
       END DO
       IF (nij.NE.proc_displs(n+1)) THEN
       WRITE(impierror,*)   "nij.ne.proc_displs(n+1) in correct_rsh",nij,proc_displs(n+1)
         CALL MPI_FINALIZE(ierr_mpi) ; STOP
       END IF
     END DO
     DEALLOCATE(e)
!WRITE(impierror,*)   "no error with exchange_rsh"
     RETURN
   ELSE  !this is a subroutin for the case there is only wet points
!correct on row direction(make imin:imax complete on minrproc)
   CALL MPI_TYPE_CONTIGUOUS(1,type_mpi_rsh,TYPE_NDIM, ierr_mpi )
   CALL MPI_TYPE_COMMIT(TYPE_NDIM, ierr_mpi )
   IF (myrproc==minrproc) THEN
     ALLOCATE(C(ndim,imin:imax,ljmin:ljmax))
     C=0.0_rsh
     C(:,limin:limax,ljmin:ljmax)=B(:,limin:limax,ljmin:ljmax)
   END IF

   DO x=minrproc+1,maxrproc

     IF (myrproc.EQ.minrproc) THEN
       count_get=NDIM*LI(3,x)*LJ(3,myjproc)
       ALLOCATE(D(NDIM, LI(3,x),LJ(3,myjproc)))
       dest=PROC(x,myjproc)
       etigo=x*100
       CALL MPI_RECV(D,count_get,TYPE_NDIM,dest,etigo,MPI_COMM_MARS,stat,ierr_mpi)
       C(:,LI(1,x):LI(2,x),ljmin:ljmax)=D
       DEALLOCATE(D)
     ELSE IF (myrproc.EQ.x) THEN
       count_send=NDIM*LI(3,myrproc)*LJ(3,myjproc)
       dest=PROC(minrproc,myjproc)
       etigo=x*100
       CALL MPI_SEND(B,count_send,TYPE_NDIM,dest,etigo,MPI_COMM_MARS,ierr_mpi)
     END IF
   END DO

   IF (iproc==0) A(:,imin:imax,ljmin:ljmax)=C

   IF (myrproc==minrproc) THEN
     DO  x=minjproc+1,maxjproc
       IF (iproc.EQ.0) THEN
         count_get=NDIM*(imax-imin+1)*LJ(3,x)
!must found most left PROC(:,x) which is not -1)
         DO n=miniproc,maxiproc
           dest=PROC(n,x)
           IF (dest.NE.-1)EXIT
         END DO
         etigo=x*1000
         ALLOCATE(D(NDIM, imin:imax,LJ(3,x)))
         CALL MPI_RECV(D(:,:,:),count_get,TYPE_NDIM,dest,etigo,MPI_COMM_MARS,stat,ierr_mpi)
         A(:,:,LJ(1,x):LJ(2,x))=D
         DEALLOCATE(D)
       ELSE IF (myjproc==x) THEN
         count_send=NDIM*(imax-imin+1)*LJ(3,myjproc)
         dest=0
         etigo=x*1000
         CALL MPI_SEND(c(:,:,:),count_send,TYPE_NDIM,dest,etigo,MPI_COMM_MARS,ierr_mpi)
       END IF
     END DO
   END IF

   IF (myrproc==minrproc) THEN
     DEALLOCATE(C)
   END IF

! sendcount=LJ(3,myjproc)
! recvcounts= LJ(3,:)
! displs= LJ(4,:)

! call MPI_GATHERV(B,sendcount , TYPE_NDIM, &
!      & A, recvcounts, displs,     TYPE_NDIM, 0,MPI_COMM_MARS, ierr_mpi)
!      if(ierr_mpi.ne.0) print*,              ierr_mpi,"ierror at exchange"

   CALL MPI_TYPE_FREE(TYPE_NDIM, ierr_mpi )
  RETURN
  ENDIF

  END SUBROUTINE correct_rsh

!!======================================================================

  SUBROUTINE exchange_j_rsh(A,NDIM) !not yet tested!!

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE exchange_j_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : A
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,ljmin,ljmax,limin,limax
   USE parameters, ONLY : jmin,jmax
   IMPLICIT NONE

!! * Arguments
   INTEGER                                             ,INTENT(in)    :: ndim
   REAL(kind=rsh),DIMENSION(NDIM,limin:limax,jmin:jmax),INTENT(inout) :: A

!! * Local declarations
   INTEGER                                                :: ierr_mpi
   INTEGER                                                :: x,y
   INTEGER                                                :: dest,size_give,size_get
   REAL(kind=rsh),DIMENSION(NDIM,limin:limax,ljmin:ljmax) :: B
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:)            :: GIVE,GET

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   B=A(:,:,ljmin:ljmax)
   size_give=ndim*LI(3,myiproc)*LJ(3,myjproc)
   ALLOCATE(GIVE(ndim,LI(3,myiproc),LJ(3,myjproc)))
   doA:DO x=mincbudy,maxcbudy
     y=cbudy(x)
     IF (y.LT.mincproc) CYCLE doA
     IF (y.GT.maxcproc) CYCLE doA
     IF (y.EQ.myjproc) THEN
       A(:,:,ljmin:ljmax)=B
       CYCLE doA
     END IF
     give=B
     size_get=ndim*LI(3,myiproc)*LJ(3,y)
     ALLOCATE(get(ndim,LI(3,myiproc),LJ(3,y)))
     dest=proc(myiproc,y)
!WRITE(impierror,*)   "i am in  exchange_j_rsh"
!WRITE(impierror,*)   "in  exchange_j_rsh just before sendrecv",myiproc,y,dest
     CALL sendrecv_rsh(dest,size_give,give,size_get,get)
!WRITE(impierror,*)   "in  exchange_j_rsh just after sendrecv"
     A(:,:,LJ(1,y):LJ(2,y))=get
     DEALLOCATE(get)
   END DO doA

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

   RETURN

  END SUBROUTINE exchange_j_rsh

!!======================================================================

  SUBROUTINE exchange_i_rsh(A,NDIM) !not yet tested!!

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE exchange_i_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : A
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,ljmin,ljmax,limin,limax
   USE parameters, ONLY : imin,imax
   IMPLICIT NONE

!! * Arguments
   INTEGER                                             ,INTENT(in)    :: ndim
   REAL(kind=rsh),DIMENSION(NDIM,imin:imax,ljmin:ljmax),INTENT(inout) :: A

!! * Local declarations
   INTEGER                                                :: ierr_mpi
   INTEGER                                                :: x,y
   INTEGER                                                :: dest,size_give,size_get
   REAL(kind=rsh),DIMENSION(NDIM,limin:limax,ljmin:ljmax) :: B
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:)            :: GIVE,GET

!!----------------------------------------------------------------------
!! * Executable part


   CALL second_mpi(ztime1)

   B=A(:,limin:limax,:)
   size_give=ndim*LI(3,myiproc)*LJ(3,myjproc)
   ALLOCATE(GIVE(ndim,LI(3,myiproc),LJ(3,myjproc)))
   doA:DO x=minrbudy,maxrbudy
     y=rbudy(x)
     IF (y.LT.minrproc) CYCLE doA
     IF (y.GT.maxrproc) CYCLE doA
     IF (y.EQ.myiproc) CYCLE doA
     give=B
     size_get=ndim*LI(3,y)*LJ(3,myjproc)
     ALLOCATE(get(ndim,LI(3,y),LJ(3,myjproc)))
     dest=proc(y,myjproc)
     CALL sendrecv_rsh(dest,size_give,give,size_get,get)
     A(:,LI(1,y):LI(2,y),:)=get
     DEALLOCATE(get)
   END DO doA

   CALL second_mpi(ztime2)
   xt_mpi_com = xt_mpi_com + ztime2 - ztime1

   RETURN

  END SUBROUTINE exchange_i_rsh

!!======================================================================

  SUBROUTINE exchange_rsh(A,NDIM)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE exchange_rsh  ***
!&E
!&E ** Purpose :
!&E
!&E ** Modified variables : A
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,ljmin,ljmax,limin,limax
   USE parameters, ONLY : imin,imax,jmin,jmax
   IMPLICIT NONE

!! * Arguments
   INTEGER                                           ,INTENT(in)    :: ndim
   REAL(kind=rsh),DIMENSION(NDIM,imin:imax,jmin:jmax),INTENT(inout) :: A

!! * Local declarations
   INTEGER                                     :: ierr_mpi
   INTEGER                                     :: i,j,nsize,n,nij,ni,nj
   INTEGER                                     :: TYPE_NDIM,sendcount
   INTEGER                                     :: size_give,x,y,size_get,dest
   INTEGER,DIMENSION(nprocs)                   :: recvcounts,displs
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: B,C,give,get
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:)   :: D

!!----------------------------------------------------------------------
!! * Executable part

   IF (maxiproc==0.AND.maxjproc==0) THEN
     RETURN
   ELSE IF ((maxiproc==0).AND.(maxjproc.NE.0)) THEN

! a way of avoiding using allgatherv
!        allocate(C(NDIM,imin:imax,ljmin:ljmax))
!         c(:,:,ljmin:ljmax)=A(:,:,ljmin:ljmax)
!          call exchange_j_rsh(C,ndim)
!        A=C
!        deallocate(c)
!        return

     CALL MPI_TYPE_CONTIGUOUS(NDIM*(imax-imin+1),type_mpi_rsh,TYPE_NDIM, ierr_mpi )
     CALL MPI_TYPE_COMMIT(TYPE_NDIM, ierr_mpi )

     sendcount=LJ(3,myjproc)
     recvcounts(1:nprocs)= LJ(3,minjproc:maxjproc)
     displs(1:nprocs)= LJ(4,minjproc:maxjproc)
     ALLOCATE(C(NDIM,1:imax-imin+1,LJ(3,myjproc)))
     C(:,:,:)=A(:,:,ljmin:ljmax)
     CALL MPI_ALLGATHERV(C,sendcount , TYPE_NDIM, &
       & A, recvcounts, displs,     TYPE_NDIM, MPI_COMM_MARS, ierr_mpi)
     IF (ierr_mpi.NE.0) WRITE(impierror,*)              ierr_mpi,"ierror at exchange"
     DEALLOCATE(c)
     CALL MPI_TYPE_FREE(TYPE_NDIM, ierr_mpi )
     RETURN

   ELSE IF (maxjproc==0.AND.maxiproc.NE.0) THEN  !not yet tested
     CALL exchange_i_rsh(A,ndim)
     RETURN

   ELSE IF (nbblank.NE.0) THEN     !tested with nblank.eq.0,only with ljmin:ljmax
!tinatestWRITE(impierror,*)   "i am in nbblanc.ne0 exchange"

     CALL MPI_TYPE_CONTIGUOUS(NDIM,type_mpi_rsh,TYPE_NDIM, ierr_mpi )
     CALL MPI_TYPE_COMMIT(TYPE_NDIM, ierr_mpi )
     nsize=SUM(proc_size(minp:maxp))
!tinatestWRITE(impierror,*)   "nsize",nsize
     ALLOCATE(C(NDIM,LI(3,myiproc),LJ(3,myjproc)))
!test tina
!tinatest DO j=ljmin,ljmax
!tinatest     DO i=limin,limax
!tinatest      A(1,i,j)=dble(1000*i+j)
!tinatest    END DO
!tinatest  ENDDO
     C(:,:,:)=A(:,limin:limax,ljmin:ljmax)  !data that i hold, which will be send to everyone
     ALLOCATE(D(NDIM,nsize))                ! data box that i keep all the data from everyone
     sendcount=proc_size(iproc)
!tinatest WRITE(impierror,*)   "sendcount",sendcount
     recvcounts(1:nprocs)=proc_size(minp:maxp)
     displs(1:nprocs)=proc_displs(minp:maxp)
!tinatest WRITE(impierror,*)   "just before allgatherv"
     CALL MPI_ALLGATHERV(C,sendcount , TYPE_NDIM, &
       & D, recvcounts, displs,     TYPE_NDIM, MPI_COMM_MARS, ierr_mpi)
     IF (ierr_mpi.NE.0) WRITE(impierror,*)              ierr_mpi,"ierror at exchange"
     CALL MPI_TYPE_FREE(TYPE_NDIM, ierr_mpi )
     DEALLOCATE(c)
!tinatest WRITE(impierror,*)   "just after allgatherv"
     nij=0
     DO n=minp,maxp
!IF (n==iproc)CYCLE
!nij=proc_displs(N)
! here i can also use function reshape
       DO j=proc_ljmin(n),proc_ljmax(n)
         DO i=proc_limin(n),proc_limax(n)
           nij=nij+1
           A(:,i,j)=D(:,nij)
!tinatest if(A(1,i,j).ne.dble(1000*i+j))then
!tinatest WRITE(impierror,*)   "a.ne.dble",a(1,i,j),i,j
!tinatest CALL MPI_FINALIZE(ierr_mpi) ; STOP
!tinatest else
!tinatest WRITE(impierror,*)   "a.eq.dble",a(1,i,j),i,j
!tinatest END IF
         END DO
       END DO
       IF (nij.NE.proc_displs(n+1)) THEN
       WRITE(impierror,*)   "nij.ne.proc_displs(n+1)",nij,proc_displs(n+1)
         CALL MPI_FINALIZE(ierr_mpi) ; STOP
       END IF
     END DO
     DEALLOCATE(d)
!tinatest WRITE(impierror,*)   "no error with exchange_rsh"
     RETURN
   ELSE  !this is a subroutin for the case there is only wet points

!WRITE(impierror,*)   "i am  iafter else exchange,nbblank",nbblank
     ALLOCATE(C(NDIM,LI(3,myiproc),jmin:jmax))   !not yet tested
     IF (mincproc==maxcproc) THEN
       C(:,:,ljmin:ljmax)=A(:,limin:limax,ljmin:ljmax)
     ELSE
!WRITE(impierror,*)   "just before exchangej-rsh"
       C(:,:,ljmin:ljmax)=A(:,limin:limax,ljmin:ljmax)
       CALL exchange_j_rsh(C,ndim)
!WRITE(impierror,*)   "just after exchangej-rsh"
     END IF
!WRITE(impierror,*)   "just before mpi_barrier"
     call MPI_BARRIER(MPI_COMM_MARS,ierr_mpi)
!WRITE(impierror,*)   "just after mpi_barrier"
     size_give=ndim*(jmax-jmin+1)*LI(3,myiproc)
     ALLOCATE(GIVE(ndim,LI(3,myiproc),jmin:jmax))
     doA:    DO x=minibudy,maxibudy
       y=ibudy(x)
       IF (y.LT.miniproc) CYCLE doA
       IF (y.GT.maxiproc) CYCLE doA
       IF (y.EQ.myiproc) THEN
         A(:,limin:limax,:)=C
         CYCLE doA
       END IF
       give=C
       size_get=ndim*(jmax-jmin+1)*LI(3,y)
       ALLOCATE(get(ndim,LI(3,y),jmin:jmax))
       dest=proc(y,myjproc)
!WRITE(impierror,*)   "just before sendrecv"
       CALL sendrecv_rsh(dest,size_give,give,size_get,get)
!WRITE(impierror,*)   "just after sendrecv"
       A(:,LI(1,y):LI(2,y),jmin:jmax)=get
       DEALLOCATE(get)
     END DO doA
     DEALLOCATE(give)
     DEALLOCATE(c)


     RETURN
   END IF

  END SUBROUTINE exchange_rsh

!!======================================================================

  SUBROUTINE ADD_MPI_RSH(mine,total)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ADD_MPI_RSH  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh

!! * Arguments
   REAL(kind=rsh),INTENT(in ) ::  mine
   REAL(kind=rsh),INTENT(out) ::  total

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

    CALL MPI_REDUCE(mine,total,1,type_mpi_rsh,MPI_SUM,0,MPI_COMM_MARS,ierr_mpi)

   END SUBROUTINE ADD_MPI_RSH

!!======================================================================

  SUBROUTINE ADD_MPI_RLG(mine,total)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ADD_MPI_RLG  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rlg

!! * Arguments
   REAL(kind=rlg),INTENT(in ) ::  mine
   REAL(kind=rlg),INTENT(out) ::  total

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

    CALL MPI_REDUCE(mine,total,1,type_mpi_rlg,MPI_SUM,0,MPI_COMM_MARS,ierr_mpi)

   END SUBROUTINE ADD_MPI_RLG

!!======================================================================

   SUBROUTINE ADD_MPI_INT(mine,total)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ADD_MPI_INT  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Arguments
   INTEGER,INTENT(in ) ::  mine
   INTEGER,INTENT(out) ::  total

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

   CALL MPI_REDUCE(mine,total,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_MARS,ierr_mpi)

  END SUBROUTINE ADD_MPI_INT

!!======================================================================

  SUBROUTINE BCAST_MPI_LOGIC(val)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE BCAST_MPI_LOGIC  ***
!&E
!&E ** Purpose : The MASTER sends a locical value to every cpus
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   LOGICAL,INTENT(inout) ::  val

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

   CALL MPI_BCAST(val,1,MPI_LOGICAL,0,MPI_COMM_MARS,ierr_mpi)

  END SUBROUTINE  BCAST_MPI_LOGIC

!!======================================================================

  SUBROUTINE BCAST_MPI_INT(val)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE BCAST_MPI_INT  ***
!&E
!&E ** Purpose : The MASTER sends an integer value to every cpus
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh
   IMPLICIT NONE

!! * Arguments
   INTEGER,INTENT(inout) ::  val

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

   CALL MPI_BCAST(val,1,MPI_INTEGER,0,MPI_COMM_MARS,ierr_mpi)

  END SUBROUTINE  BCAST_MPI_INT

!!======================================================================

  SUBROUTINE BCAST_MPI_RSH(val)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE BCAST_MPI_RSH  ***
!&E
!&E ** Purpose : The MASTER sends a RSH (single precision) value to every cpus
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh
   IMPLICIT NONE

!! * Arguments
   REAL(KIND=4),INTENT(inout) ::  val

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

   CALL MPI_BCAST(val,1,MPI_REAL,0,MPI_COMM_MARS,ierr_mpi)

  END SUBROUTINE  BCAST_MPI_RSH

!!======================================================================

  SUBROUTINE BCAST_MPI_RLG(val)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE BCAST_MPI_RLG  ***
!&E
!&E ** Purpose : The MASTER sends a RLG (double precision) value to every cpus
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rlg
   IMPLICIT NONE

!! * Arguments
   REAL(KIND=8),INTENT(inout) ::  val

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

   CALL MPI_BCAST(val,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_MARS,ierr_mpi)

  END SUBROUTINE  BCAST_MPI_RLG

!!======================================================================

  SUBROUTINE BCAST_MPI_CHAR(val)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE BCAST_MPI_CHAR  ***
!&E
!&E ** Purpose : The MASTER sends an integer value to every cpus
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   IMPLICIT NONE

!! * Arguments
   CHARACTER(LEN=*),INTENT(inout) ::  val

!! * Local declarations
   INTEGER :: ierr_mpi

!!----------------------------------------------------------------------
!! * Executable part

   CALL MPI_BCAST(val,1,MPI_CHAR,0,MPI_COMM_MARS,ierr_mpi)

  END SUBROUTINE  BCAST_MPI_CHAR

!!======================================================================

  SUBROUTINE FINDMAX_MPI_3(LV,LLI,LLJ)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE FINDMAX_MPI_3  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh
   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),INTENT(inout) :: LV
   INTEGER       ,INTENT(inout) :: LLI,LLJ

!! * Local declarations
   INTEGER                                   :: ierr_mpi
   REAL(kind=rsh)                            :: LA,LB,LC,BIG,WHYSV
   INTEGER                                   :: k,kbig,WHYSI(2)
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:)   :: V
   INTEGER       ,ALLOCATABLE,DIMENSION(:,:) :: I

!!----------------------------------------------------------------------
!! * Executable part

   WHYSV=LV
   ALLOCATE(V(nprocs))
   CALL MPI_ALLGATHER(WHYSV,1,type_mpi_rsh,V(1), 1,type_mpi_rsh,MPI_COMM_MARS, ierr_mpi)
   WHYSI(1)=LLI
   WHYSI(2)=LLJ
   ALLOCATE(I(2,nprocs))
   CALL MPI_ALLGATHER(WHYSI,2,MPI_INTEGER,I,2,MPI_INTEGER,MPI_COMM_MARS, ierr_mpi)

   BIG=V(1)
   kBIG=1
   DO k=1,NPROCS
     IF (V(k).GT.BIG) THEN
       kBIG=k
       BIG=V(k)
     END IF
   END DO
   LV=V(kBIG)
   LLI=I(1,kBIG)
   LLJ=I(2,kBIG)
   DEALLOCATE(V,I)
   RETURN

  END SUBROUTINE FINDMAX_MPI_3

!!======================================================================



  SUBROUTINE FINDMAX_MPI_6(LV,LLI,LLJ,LLK,LA,LB,LC)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE FINDMAX_MPI_6  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh
   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),INTENT(inout) :: LV,LA,LB,LC
   INTEGER       ,INTENT(inout) :: LLI,LLJ,LLK

!! * Local declarations
   INTEGER                                   :: ierr_mpi
   REAL(kind=rsh)                            :: BIG,WHYSV(4)
   INTEGER                                   :: k,kbig,WHYSI(3)
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:) :: V
   INTEGER       ,ALLOCATABLE,DIMENSION(:,:) :: I

!!----------------------------------------------------------------------
!! * Executable part

   WHYSV(1)=LV
   WHYSV(2)=LA
   WHYSV(3)=LB
   WHYSV(4)=LC
   ALLOCATE(V(4,nprocs))
   CALL MPI_ALLGATHER(WHYSV,4,type_mpi_rsh, V(1,1),4,type_mpi_rsh,MPI_COMM_MARS,ierr_mpi)
   WHYSI(1)=LLI
   WHYSI(2)=LLJ
   WHYSI(3)=LLK
   ALLOCATE(I(3,nprocs))
   CALL MPI_ALLGATHER(WHYSI,3,MPI_INTEGER,I,3,MPI_INTEGER,MPI_COMM_MARS, ierr_mpi)

!use MPI_PACK so that the program only need to send one value? combine u and v
!max search??

   BIG=V(1,1)
   kBIG=1
   DO k=1,NPROCS
     IF (V(1,k).GT.BIG) THEN
       kBIG=k
       BIG=V(1,k)
     END IF
   END DO
   LV=V(1,kBIG)
   LA=V(2,kBIG)
   LB=V(3,kBIG)
   LC=V(4,kBIG)
   LLI=I(1,kBIG)
   LLJ=I(2,kBIG)
   LLK=I(3,kBIG)
   DEALLOCATE(V,I)
   RETURN

  END SUBROUTINE FINDMAX_MPI_6


  SUBROUTINE FINDMAX_MPI(LV,LLI,LLJ,LA,LB,LC)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE FINDMAX_MPI  ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh
   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),INTENT(inout) :: LV,LA,LB,LC
   INTEGER       ,INTENT(inout) :: LLI,LLJ

!! * Local declarations
   INTEGER                                   :: ierr_mpi
   REAL(kind=rsh)                            :: BIG,WHYSV(4)
   INTEGER                                   :: k,kbig,WHYSI(2)
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:) :: V
   INTEGER       ,ALLOCATABLE,DIMENSION(:,:) :: I

!!----------------------------------------------------------------------
!! * Executable part

   WHYSV(1)=LV
   WHYSV(2)=LA
   WHYSV(3)=LB
   WHYSV(4)=LC
   ALLOCATE(V(4,nprocs))
   CALL MPI_ALLGATHER(WHYSV,4,type_mpi_rsh,V(1,1), 4,type_mpi_rsh,MPI_COMM_MARS, ierr_mpi)
   WHYSI(1)=LLI
   WHYSI(2)=LLJ
   ALLOCATE(I(2,nprocs))
   CALL MPI_ALLGATHER(WHYSI,2,MPI_INTEGER,I,2,MPI_INTEGER,MPI_COMM_MARS, ierr_mpi)

!use MPI_PACK so that the program only need to send one value? combine u and v
!max search??

   BIG=V(1,1)
   kBIG=1
   DO k=1,NPROCS
     IF (V(1,k).GT.BIG) THEN
       kBIG=k
       BIG=V(1,k)
     END IF
   END DO
   LV=V(1,kBIG)
   LA=V(2,kBIG)
   LB=V(3,kBIG)
   LC=V(4,kBIG)
   LLI=I(1,kBIG)
   LLJ=I(2,kBIG)
   DEALLOCATE(V,I)
   RETURN

  END SUBROUTINE FINDMAX_MPI

!!======================================================================

  SUBROUTINE CREATE_COMM_ZONAL(COMM_ZONAL,comm2dverif)

   INTEGER, INTENT(INOUT)     :: comm_zonal
   INTEGER, INTENT(INOUT), OPTIONAL     :: comm2dverif

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER                    :: comm1d,comm2d
   INTEGER,PARAMETER          :: ndim2d=2
   INTEGER, DIMENSION(ndim2d) :: dim2d
   LOGICAL, DIMENSION(ndim2d) :: period,subdivision
   LOGICAL                    :: reorder

! create a topology (2D)
   dim2d(2)=MAXIPROC+1
   dim2d(1)=MAXJPROC+1
   period(:)=.FALSE.
   reorder=.FALSE.
   CALL MPI_CART_CREATE(MPI_COMM_MARS,ndim2d,dim2d,period,reorder,comm2d,ierr_mpi)

! moyenne zonale : chaque ligne de la topologie 2D (ou ligne de processeurs) est une topologie cartesienne 1D
! Subdivisions inversees car Tina a range les cpus dans ordre i croisant et non j croissant
   subdivision(1)=.FALSE.   ! Tina a range les cpus dans ordre i croisant et non j croissant
   subdivision(2)=.TRUE.
   CALL MPI_CART_SUB(comm2d,subdivision,comm1d,ierr_mpi)
!print*, 3953,"RANK",iproc," DANS COMM",comm1d ! negatif et identique pour tout proc !!!

   comm_zonal=comm1d

   !PRINT*,'exit CREATE_COMM_ZONAL'

  END SUBROUTINE CREATE_COMM_ZONAL

!!======================================================================

  SUBROUTINE ADD_ZONAL_RSH(mine,mind1,maxd1,mind2,maxd2,comm1d)
   USE parameters, ONLY : rsh

   IMPLICIT NONE

   REAL(KIND=rsh), DIMENSION(mind1:maxd1,mind2:maxd2), INTENT(inout) ::  mine
   INTEGER, INTENT(in) :: mind1,maxd1,mind2,maxd2,comm1d

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER, DIMENSION(1)            :: dim12
   INTEGER                          :: dimmpi
   REAL(KIND=rsh), DIMENSION((maxd1-mind1+1)*(maxd2-mind2+1)) ::  vecmine, vectotal

   !PRINT*,'enter ADD_ZONAL_RSH',mind1,maxd1,mind2,maxd2,comm1d

   dim12=(maxd1-mind1+1)*(maxd2-mind2+1)
   vecmine(:)=RESHAPE(mine,dim12)
   dimmpi = dim12(1)    !cannot be a vector in MPI_ALLREDUCE
   CALL MPI_ALLREDUCE(vecmine,vectotal,dimmpi,type_mpi_rsh,MPI_SUM,comm1d,ierr_mpi)
   mine(:,:)=RESHAPE(vectotal,SHAPE(mine))

   !PRINT*,'exit ADD_ZONAL_RSH'

  END SUBROUTINE ADD_ZONAL_RSH

!!======================================================================

  SUBROUTINE ADD_ALL_MPI_RSH(value)
   USE parameters, ONLY : rsh
   IMPLICIT NONE
   REAL(KIND=rsh), INTENT(inout) :: value

!! * Local declarations
   INTEGER :: ierr_mpi
   REAL(KIND=rsh) :: value_temp

! MPI_IN_PLACE fait planter MPI_ALL_REDUCE avec certaines versions de MPI
!CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,type_mpi_rsh,MPI_SUM,MPI_COMM_MARS,ierr_mpi)
   CALL MPI_ALLREDUCE(value,value_temp,1,type_mpi_rsh,MPI_SUM,MPI_COMM_MARS,ierr_mpi)
   value=value_temp

  END SUBROUTINE ADD_ALL_MPI_RSH

!!======================================================================

  SUBROUTINE ADD_ALL_MPI_INT(value)
   IMPLICIT NONE
   INTEGER, INTENT(inout) :: value

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: value_temp

!CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,MPI_INTEGER,MPI_SUM,MPI_COMM_MARS,ierr_mpi)
   CALL MPI_ALLREDUCE(value,value_temp,1,MPI_INTEGER,MPI_SUM,MPI_COMM_MARS,ierr_mpi)
   value=value_temp

  END SUBROUTINE  ADD_ALL_MPI_INT

!!======================================================================

  SUBROUTINE MAX_ALL_MPI_RSH(value,comm_active)
   USE parameters, ONLY : rsh
   IMPLICIT NONE
   REAL(KIND=rsh), INTENT(inout)        :: value
   INTEGER,        INTENT(in), OPTIONAL :: comm_active

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: communicator
   REAL(KIND=rsh) :: value_temp

   IF (PRESENT(comm_active)) THEN
     communicator=comm_active
   ELSE
     communicator=MPI_COMM_MARS
   END IF
!CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,type_mpi_rsh,MPI_MAX,communicator,ierr_mpi)
   CALL MPI_ALLREDUCE(value,value_temp,1,type_mpi_rsh,MPI_MAX,communicator,ierr_mpi)
   value=value_temp

  END SUBROUTINE  MAX_ALL_MPI_RSH

!!======================================================================

  SUBROUTINE MAX_ALL_MPI_INT(value,comm_active)
   IMPLICIT NONE
   INTEGER, INTENT(inout)        :: value
   INTEGER, INTENT(in), OPTIONAL :: comm_active

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: communicator
   INTEGER :: value_temp

   IF (PRESENT(comm_active)) THEN
     communicator=comm_active
   ELSE
     communicator=MPI_COMM_MARS
   END IF
!CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,MPI_INTEGER,MPI_MAX,communicator,ierr_mpi)
   CALL MPI_ALLREDUCE(value,value_temp,1,MPI_INTEGER,MPI_MAX,communicator,ierr_mpi)
   value=value_temp

  END SUBROUTINE MAX_ALL_MPI_INT

!!======================================================================

  SUBROUTINE MIN_ALL_MPI_RSH(value,comm_active)
   USE parameters, ONLY : rsh
   REAL(KIND=rsh), INTENT(inout)        :: value
   INTEGER,        INTENT(in), OPTIONAL :: comm_active

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: communicator
   REAL(KIND=rsh) :: value_temp

   IF (PRESENT(comm_active)) THEN
     communicator=comm_active
   ELSE
     communicator=MPI_COMM_MARS
   END IF
!CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,type_mpi_rsh,MPI_MIN,communicator,ierr_mpi)
   CALL MPI_ALLREDUCE(value,value_temp,1,type_mpi_rsh,MPI_MIN,communicator,ierr_mpi)
   value=value_temp

  END SUBROUTINE  MIN_ALL_MPI_RSH

!!======================================================================

  SUBROUTINE MIN_ALL_MPI_INT(value,comm_active)
!! * Arguments
   INTEGER, INTENT(inout)        :: value
   INTEGER, INTENT(in), OPTIONAL :: comm_active

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: communicator
   INTEGER :: value_temp

   IF (PRESENT(comm_active)) THEN
     communicator=comm_active
   ELSE
     communicator=MPI_COMM_MARS
   END IF
!CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,MPI_INTEGER,MPI_MIN,communicator,ierr_mpi)
   CALL MPI_ALLREDUCE(value,value_temp,1,MPI_INTEGER,MPI_MIN,communicator,ierr_mpi)
   value=value_temp

  END SUBROUTINE MIN_ALL_MPI_INT

!!======================================================================

  SUBROUTINE exchange_vectcpu_int(A)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE EXCHANGE_VECTCPU_INT ***
!&E
!&E ** Purpose :
!&E
!&E---------------------------------------------------------------------
!! * Arguments
   INTEGER, DIMENSION(0:nprocs-1), INTENT(inout) :: A

!! * Local declarations
   INTEGER :: ierr_mpi
   INTEGER :: V

!!----------------------------------------------------------------------
!! * Executable part

   V=A(iproc)
   CALL MPI_ALLGATHER(V,1,MPI_INTEGER,A,1,MPI_INTEGER,MPI_COMM_MARS,ierr_mpi)

  END SUBROUTINE exchange_vectcpu_int

!!======================================================================

  SUBROUTINE SECOND_MPI(XT)

!! * Arguments
   REAL(KIND=rlg), DIMENSION(2), INTENT(out) :: XT

!!----------------------------------------------------------------------
!! * Executable part

   CALL CPU_TIME(XT(1))
   XT(2) = MPI_Wtime()

  END SUBROUTINE SECOND_MPI

!!======================================================================



 END MODULE toolmpi
