PROGRAM MD
      IMPLICIT NONE
      include 'mpif.h'
      CHARACTER*100 Buffer
      REAL*8,  DIMENSION(:,:,:,:),ALLOCATABLE :: PARTCOORDS
      REAL*8,  DIMENSION(:,:,:,:),ALLOCATABLE :: PARTFORCES
      REAL*8,  DIMENSION(:,:,:,:),ALLOCATABLE :: PARTACCELS
      REAL*8,  DIMENSION(:,:,:,:),ALLOCATABLE :: PARTVELS
      INTEGER, DIMENSION(:,:),    ALLOCATABLE :: CellSizes
      INTEGER, DIMENSION(:,:,:),  ALLOCATABLE :: CellMap 
      INTEGER, DIMENSION(:,:,:),  ALLOCATABLE :: ABSCOORD
      REAL*8,  DIMENSION(:,:,:,:),ALLOCATABLE :: REMOTE_PARTS
      INTEGER, DIMENSION(:,:),    ALLOCATABLE :: NREQR
      INTEGER, DIMENSION(:,:),    ALLOCATABLE :: NREQS
      INTEGER, DIMENSION(:,:,:),  ALLOCATABLE :: REQR
      INTEGER, DIMENSION(:,:,:),  ALLOCATABLE :: REQS
      REAL*8,  DIMENSION(:),      ALLOCATABLE :: PARTBUF
      INTEGER  nCellRows,NCellCols,CellSize,NumIterations
      INTEGER  nParticles,Dimension,nPartPerCell
      REAL*8 ran,centerx,centery
      REAL*8 rx1,ry1,rx,ry,r,rt
      REAL*8 TotPot1,TotPot2,LSUm, TotKin
      REAL*8 TotPotSum,TotKinSum
      REAL*8 ARmass,AReps,ARsigma,CellDim
      REAL*8 timestart,timeend
      INTEGER  nPartsPerCell,MAXPART,MINPART
      INTEGER TotalCells,TotSize
      INTEGER i,j,k,l,cx1,cy1,cx2,cy2,t
      INTEGER maxCellSize
      INTEGER ierror,NCPUS,MY_PE,SQncpus
      INTEGER NLocRows,NLocRowsShort,NLocRowsPE
      INTEGER NLocCols,NLocColsShort,NLocColsPE
      INTEGER icpu,ixend,iyend,ixcount,iycount
      INTEGER ix,iy,ixx,iyy
      INTEGER absx,absy
      PARAMETER(CellDim=12.0)
      PARAMETER(NPartPerCell=10)
      PARAMETER(ARMass=39.948)
 
      CALL MPI_INIT(IERROR)
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NCPUS,ierror)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,MY_PE,ierror)
      
      IF(IARGC() < 2) THEN
         WRITE(6,*)"Incorrect syntax: should be two arguments"
         CALL EXIT(2)
      ENDIF
      
      ! Number of CPUs must be a Perfect Square
      SQncpus = int(sqrt(real(NCPUS)))
      if(MOD(NCPUS,SQncpus).ne.0) then
         if(MY_PE == 0) then
            write(6,*)'The Number of cores must be a perfect square',
     $           '(e.g. 4,16,64,etc...)'
         endif
         CALL MPI_FINALIZE(ierror)
         CALL exit(3)
      endif
      
      CALL GETARG(1,Buffer)
      READ(Buffer,*)nParticles
      CALL GETARG(2,Buffer)
      READ(Buffer,*)NumIterations
      
      ! Lets Allocate the number of cells
      ! Dimension = The total dimensions of the Cells
      Dimension= int(Sqrt(NParticles/real(NPartPerCell)))
      TotalCells = Dimension**2
      ! The values describing the local data are dependent on 
      ! whether the Dimension and the Square root of the number of 
      ! cores is perfectly divisible
      if(MOD(Dimension,SQncpus).ne.0) Then
         NLocRows      = Dimension/SQncpus 
         NLocCols      = Dimension/SQncpus
         NLocColsShort = Dimension-((SQncpus-1)*NLocCols)
         NLocRowsShort = Dimension-((SQncpus-1)*NLocRows)
         if(MOD(MY_PE,SQncpus) == (SQncpus-1)) then
            NLocColsPE = NLocColsShort
         else
            NLocColsPE = NLocCols
         endif
         if(MY_PE >= ((SQncpus*(SQncpus-1)))) then
            NLocRowsPE = NLocRowsShort
         else
            NLocRowsPE = NLocRows
         endif
      else
         NLocRows     = Dimension/SQncpus
         NLocCols     = Dimension/SQncpus
         NLocColsShort= 0
         NLocRowsShort= 0
         NLocRowsPE   = NLocRows
         NLocColsPE   = NLocCols
      endif
      
!! We are going to create the CellMap here 
      ALLOCATE(CellMap(4,Dimension,Dimension))
      icpu = 0
      DO I=1,SQncpus
         ixend = NLocRows
         if((i*NLocRows > Dimension).or.(i.eq.SQncpus)) then
            ixend = NLocRowsShort;
         endif
         DO J=1,SQncpus
            iyend = NLocCols
            if((j*NLocCols > Dimension).or.(j.eq.SQncpus)) then
               iyend = NLocColsShort
            endif
            ixcount = 1
            DO IX = 1,ixend
               iycount = 1
               DO IY = 1,iyend
                  ixx = (i-1)*NLocRows+ix 
                  iyy = (j-1)*NLocCols+iy
                  CellMap(1,iyy,ixx) = icpu
                  CellMap(2,iyy,ixx) = ixcount
                  CellMap(3,iyy,ixx) = iycount
                  CellMap(4,iyy,ixx) = 0
                  iycount = iycount + 1
               ENDDO
               ixcount = ixcount + 1
            ENDDO
            icpu = icpu + 1
         ENDDO
      ENDDO
      
      CALL ComputeAtomsPerCell(CellMap,Dimension,Dimension,
     $                         NParticles,maxCellSize)
      
      if(MY_PE.EQ.0) Then
         WRITE(6,'(A,I6,A,I6,A,I6,A)')' The Total Number of Cells is ',
     $        TotalCells,' With ',maxCellSize,
     $        ' particles per cell, and ',
     $        nParticles,' particles total in system'
         call flush(6)
      endif
    
!OK, now we can allocate the Particle Matrix
!We need some additional arrays here to store values for each cell
!These arrays are analogous to the Cell datastructure in the C++ example
      ALLOCATE(ABSCOORD(2,NLocColsPE,NLocRowsPE))
      ALLOCATE(REMOTE_PARTS((2*maxCellSize)+1,8,NLocColsPE,NLocRowsPE))
      ALLOCATE(NREQR(NLocColsPE,NLocRowsPE))
      ALLOCATE(NREQS(NLocColsPE,NLocRowsPE))
      ALLOCATE(REQR(8,NLocColsPE,NLocRowsPE))
      ALLOCATE(REQS(8,NLocColsPE,NLocRowsPE))
      ALLOCATE(PARTCOORDS(2,maxCellSize,NLocColsPE,NLocRowsPE))
      ALLOCATE(PARTFORCES(2,maxCellSize,NLocColsPE,NLocRowsPE))
      ALLOCATE(PARTACCELS(2,maxCellSize,NLocColsPE,NLocRowsPE))
      ALLOCATE(PARTVELS(2,maxCellSize,NLocColsPE,NLocRowsPE))
!     Zero all of the arrays
      DO i=1,NLocRowsPE
         do j=1,NLocColsPE
            NREQR(j,i) = 0
            NREQS(j,i) = 0
            DO k=1,2
               ABSCOORD(k,j,i) = 0
            ENDDO
            DO k=1,8
               REQR(k,j,i) = 0
               REQS(k,j,i) = 0
               DO l=1,2*maxCellSize+1
                  REMOTE_PARTS(l,k,j,i) = 0.0
               ENDDO
            ENDDO
            
            DO k=1,maxCellSize
               do l=1,2
                  PartCoords(l,k,j,i) = 0.0
                  PartForces(l,k,j,i) = 0.0
                  PartAccels(l,k,j,i) = 0.0
                  PartVels(l,k,j,i)   = 0.0
               enddo
            ENDDO
         enddo
      ENDDO
    
!     WE will be filling the cells, making sure than
!     No atom is less than 4 Angstroms from another

      DO i = 1,NLocRowsPE
         DO j = 1,NLocColsPE
            ABSCOORD(1,j,i) = (INT(MY_PE/SQncpus)*NLocRows)+i
            ABSCOORD(2,j,i) = (MOD(MY_PE,SQncpus)*NLocCols)+j
            
            Centerx = CellDim*real(ABSCOORD(1,j,i)) + 0.5*CellDim
            Centery = CellDim*real(ABSCOORD(2,j,i)) + 0.5*CellDim
            call random_number(ran)
            PartCoords(1,1,j,i) = centerx + ((ran-0.5)*(CellDim-2)) 
            call random_number(ran)
            PartCoords(2,1,j,i) = centery + ((ran-0.5)*(CellDim-2))
            DO k=2,CellMap(4,ABSCOORD(2,j,i),ABSCOORD(1,j,i))
               R=0
               DO WHILE (R < 4.0)
                  R=4.0001
                  call random_number(ran)
                  RX1 = centerx + ((ran-0.5)*(CellDim-2))
                  call random_number(ran)
                  RY1 = centery + ((ran-0.5)*(CellDim-2))
!     check to make sure it is far enough away from the
!     others in the cell
                  DO l=1,k-1
                     RX=RX1-PartCoords(1,l,j,i)
                     RY=RY1-PartCoords(2,l,j,i)
                     RT = RX**2+RY**2
                     IF(RT < R) THEN
                        R=RT
                     ENDIF
                  ENDDO
               ENDDO
               PartCoords(1,k,j,i) = RX1
               PartCoords(2,k,j,i) = RY1
            ENDDO
         ENDDO
      ENDDO

      timestart = MPI_Wtime()
      ! Lets Start Iterating
      WRITE(6,*)'Starting Iterations'
      DO t=1,NumIterations
         ! Zero Energy variables
         TotPot1 = 0.0
         TotPot2 = 0.0
         TotKin = 0.0
         TotKinSum=0.0
         
         
         DO CX1 = 1,NLocRowsPE
            DO CY1 = 1,NLocColsPE 
               absx= ABSCOORD(1,CY1,CX1)
               absy= ABSCOORD(2,CY1,CX1)
               CALL Communicate(CellMap,PartCoords(1,1,CY1,CX1),
     $              CellMap(4,absy,absx),Remote_Parts(1,1,CY1,CX1),
     $              AbsCoord(1,CY1,CX1),Reqr(1,CY1,CX1),
     $              Reqs(1,CY1,CX1),
     $              NReqr(CY1,CX1),NReqs(CY1,CX1),Dimension,
     $              NLocRowsPE,
     $              NLocColsPE,MaxCellSize)
            ENDDO
         ENDDO
         
! Loop over Cells
         DO CX1 = 1, NLocRowsPE
            DO CY1 = 1,NLocColsPE
               absx = ABSCOORD(1,CY1,CX1)
               absy = ABSCOORD(2,CY1,CX1)
!     Cells must be neighbors 
!     Let's do the ones within a cell first
               DO I=1,CellMap(4,absy,absx)
                  DO J=1,CellMap(4,absy,absx)
                     if(i.ne.j)then
                        CALL INTERACT(
     $                       PartCoords(1,j,cy1,cx1),
     $                       PartCoords(2,j,cy1,cx1),
     $                       PartCoords(1,i,cy1,cx1),
     $                       PartCoords(2,i,cy1,cx1),
     $                       PartForces(1,j,cy1,cx1),
     $                       PartForces(2,j,cy1,cx1),
     $                       TotPot1)
                     ENDIF
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
         
         DO CX1=1,NLocRowsPE
            DO CY1=1,NLocColsPE
               CALL POSTWAITS(REQR(1,CY1,CX1),REQS(1,CY1,CX1),
     $              NREQR(CY1,CX1),NREQS(CY1,CX1))
            ENDDO
         ENDDO
         
  ! Lets do the other cells
         DO CX1 = 1,NLocRowsPE
            DO CY1 = 1,NLocColsPE
               absx = ABSCOORD(1,CY1,CX1)
               absy = ABSCOORD(2,CY1,CX1)
               DO CX2 = 1,NREQR(CY1,CX1)
                  DO i = 1,CellMap(4,absy,absx)
                     DO j = 1,int(REMOTE_PARTS(1,CX2,CY1,CX1))
                        CALL INTERACT(
     $                       PartCoords(1,i,cy1,cx1),
     $                       PartCoords(2,i,cy1,cx1),
     $                       REMOTE_PARTS((j-1)*2+2,cx2,cy1,cx1),
     $                       Remote_PARTS((j-1)*2+3,cx2,cy1,cx1),
     $                       PartForces(1,i,cy1,cx1),
     $                       PartForces(2,i,cy1,cx1),
     $                       TotPot2)
                     ENDDO
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
         
                      ! Now apply the forces to integrate to new positions
         DO CX1=1,NLocRowsPE
            DO CY1=1,NLocColsPE
               absx = ABSCOORD(1,CY1,CX1)
               absy = ABSCOORD(2,CY1,CX1)
               DO I=1,CellMap(4,absy,absx)
                  CALL UPDATE(
     $                 PARTFORCES(1,i,cy1,cx1),
     $                 PARTFORCES(2,i,cy1,cx1),
     $                 PARTACCELS(1,i,cy1,cx1),
     $                 PARTACCELS(2,i,cy1,cx1),
     $                 PARTVELS(1,i,cy1,cx1),
     $                 PARTVELS(2,i,cy1,cx1),
     $                 PARTCOORDS(1,i,cy1,cx1),
     $                 PARTCOORDS(2,i,cy1,cx1),
     $                 ARMass,TotKin)
               ENDDO
            ENDDO
         ENDDO
! Collect the energy variables from each processor
         LSum = TotPot2+TotPot1
         CALL MPI_Reduce(LSum,TotPotSum,1,MPI_REAL8,
     $        MPI_SUM,0,MPI_COMM_WORLD,IERROR)
         CALL MPI_Reduce(TotKin,TotKinSum,1,MPI_REAL8,
     $        MPI_SUM,0,MPI_COMM_WORLD,IERROR)
         if(MY_PE==0) Then
            WRITE(6,'(A,I10,A,E20.10,A)')'Iteration ',t,
     $           ' with Total Energy',(TotPotSum+TotKinSum)/NParticles
     $           ,' Per Particle'
         ENDIF
      ENDDO
      timeend = MPI_WTime()
      IF(MY_PE.eq.0) THEN
         WRITE(6,'(A,F20.10)')'The Iteration Time is ',timeend-timestart
      ENDIF
      Call MPI_BARRIER(MPI_COMM_WORLD,IERROR)
      DEALLOCATE(PARTCOORDS)
      DEALLOCATE(PARTFORCES)
      DEALLOCATE(PARTACCELS)
      DEALLOCATE(PARTVELS)
      DEALLOCATE(CellMap)
      DEALLOCATE(ABSCOORD)
      DEALLOCATE(REMOTE_PARTS)
      DEALLOCATE(NREQR)
      DEALLOCATE(NREQS)
      DEALLOCATE(REQR)
      DEALLOCATE(REQS)      
      CALL MPI_FINALIZE(IERROR)
      END PROGRAM
      
      SUBROUTINE INTERACT(x1,y1,x2,y2,fx,fy,TP)
      IMPLICIT NONE
      REAL*8 x1,y1,x2,y2,fx,fy
      REAL*8 rx,ry,r,f,Sig6,Sig12
      REAL*8 AReps,ARsigma,TP
      PARAMETER(AReps=119.8)
      PARAMETER(ARsigma=3.405)
      
      Rx=x1-x2
      Ry=y1-y2
      R=Rx*Rx+Ry*Ry
      IF((R < 0.00001))THEN
         R=0.0
         fx=0.0
         fy=0.0
         RETURN
      ENDIF
      R = SQRT(R)
!Derivative of the potential
      Sig6 = (ARsigma/R)**6
      Sig12 = Sig6**2
      F = ((Sig12 - 0.5*Sig6)*48.0*AReps)/R
      Fx = FX + F*(RX)
      Fy = FY + F*(RY)
      TP = TP + 4.0*AReps*(Sig12-Sig6)
     
      RETURN
      END SUBROUTINE
      
      SUBROUTINE UPDATE(FX,FY,AX,AY,VX,VY,X,Y,MASS,TK)
      REAL*8 FX,FY,AX,AY,VX,VY,X,Y
      REAL*8 TP,TK,TE
      REAL*8 DT
      REAL*8 MASS
      ! We are using a 1.0 fs timestep, this is the conversion factor
      DT=0.000911633
      AX = FX/MASS
      AY = FY/MASS
      !Update Velocities
      VX = VX + 0.5*DT*AX
      VY = VY + 0.5*DT*AY
      !Update Positions
      X  = X + DT*VX
      Y  = Y + DT*VY      
      !Update Energy
      TK= TK + 0.5*MASS*(VX**2+VY**2)
      FX = 0.0
      FY = 0.0
      
      RETURN
      END

      SUBROUTINE ComputeAtomsPerCell(CellMap,NCols,NRows,
     $                               NParts,maxCellSize)
      IMPLICIT NONE
      INTEGER nCols,nRows,maxCellSize,nParts
      INTEGER CellMap(4,NCols,NRows)
      INTEGER nPartsPerCell
      INTEGER molsum,i,j
      PARAMETER(nPartsPerCell=10)
      
      maxCellSize = nPartsPerCell
      DO i=1,NRows
         DO j=1,NCols
            CellMap(4,j,i) =  nPartsPerCell
         ENDDO
      ENDDO
      molsum = NRows*NCols*NPartsPerCell
      !Diivy up the rest of the atoms amongst the cells
      if(molsum < NParts) THEN   
 10      maxCellSize = maxCellSize + 1
         Do i=1,NRows
            DO j=1,NCols
               CellMap(4,j,i) = CellMap(4,j,i) + 1
               molsum = molsum + 1
               if(molsum >=nParts) Then
                  RETURN
               ENDIF
            ENDDO
         ENDDO
         if(molsum < nParts) then
            goto 10
         ENDIF
      ENDIF
      
      RETURN
      END SUBROUTINE
      
      SUBROUTINE COMMUNICATE(CellMap,PartCoords,PSize,RemoteParticles,
     $     AbsCoords,ReqR,ReqS,NReqR,NReqS,
     $     Dims,NRows,NCols,Max)
      IMPLICIT NONE
      include 'mpif.h'
      INTEGER NRows,NCols,Dims,Max,NReqR,NReqS,i,j
      INTEGER CellMap(4,Dims,Dims)
      INTEGER AbsCoords(2),PSize,pcount,ircount,iscount
      INTEGER ReqR(8),ReqS(8),absx,absy
      REAL*8 PartCoords(2,Max)
      REAL*8 RemoteParticles(2*Max+1,8)
      REAL*8, Dimension(:), ALLOCATABLE::PartCs
      INTEGER IERROR,tag
      ALLOCATE(PartCs(max*2 + 1))
      
      PartCs(1) = Psize
      pcount = 2
      Do i=1,PSize
         PartCs(pcount) = PartCoords(1,i)
         pcount = pcount + 1
         PartCs(pcount) = PartCoords(2,i)
         pcount = pcount + 1
      ENDDO
      
      absx = AbsCoords(1)
      absy = AbsCoords(2)
      ircount = 1
!i+1,j
      if(absx.ne.Dims) then
         tag = (absx+1)*Dims+absy
         CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $        MPI_REAL8,CellMap(1,absy,absx+1),tag,
     $        MPI_COMM_WORLD,
     $        ReqR(ircount),IERROR)
         ircount = ircount + 1
      endif
!i-1,j
      if(absx.ne.1) then
         tag = (absx-1)*Dims+absy
         CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $        MPI_REAL8,CellMap(1,absy,absx-1),tag,MPI_COMM_WORLD,
     $        ReqR(ircount),IERROR)
         ircount = ircount + 1
      endif
!     i,j+1
      if(absy.ne.Dims) then
         tag = (absx)*Dims+absy+1
         CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $        MPI_REAL8,CellMap(1,absy+1,absx),tag,MPI_COMM_WORLD,
     $        ReqR(ircount),IERROR)
         ircount = ircount + 1
      endif
!     i,j-1
      if(absy.ne.1) then
         tag = (absx)*Dims+absy-1
         CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $        MPI_REAL8,CellMap(1,absy-1,absx),tag,MPI_COMM_WORLD,
     $        ReqR(ircount),IERROR)
         ircount = ircount + 1
      endif
!     i+1,j+1
      if((absx.ne.Dims).and.(absy.ne.Dims)) then
         tag = (absx+1)*Dims+absy+1
         CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $        MPI_REAL8,CellMap(1,absy+1,absx+1),tag,
     $        MPI_COMM_WORLD,
     $        ReqR(ircount),IERROR)
         ircount = ircount + 1
      endif
!     i+1,j-1
      if((absx.ne.Dims).and.(absy.ne.1)) then
         tag = (absx+1)*Dims+absy-1
         CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $        MPI_REAL8,CellMap(1,absy-1,absx+1),tag,
     $        MPI_COMM_WORLD,
     $        ReqR(ircount),IERROR)
         ircount = ircount + 1
      endif
!     i-1,j+1
      if((absx.ne.1).and.(absy.ne.Dims)) then
         tag = (absx-1)*Dims+absy+1
         CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $        MPI_REAL8,CellMap(1,absy+1,absx-1),tag,
     $        MPI_COMM_WORLD,
     $        ReqR(ircount),IERROR)
         ircount = ircount + 1
      endif
!i-1,j-1
         if((absx.ne.1).and.(absy.ne.1)) then
            tag = (absx-1)*Dims+absy-1
            CALL MPI_IRECV(RemoteParticles(1,ircount),2*Max+1,
     $           MPI_REAL8,CellMap(1,absy-1,absx-1),tag,MPI_COMM_WORLD,
     $           ReqR(ircount),IERROR)
            ircount = ircount + 1
         endif
         
         iscount = 1
      
      tag = absx*Dims+absy
!i+1,j
      if(absx.ne.Dims) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,CellMap(1,absy,absx+1),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif
!i-1,j
      if(absx.ne.1) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,CellMap(1,absy,absx-1),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif
!i,j+1
      if(absy.ne.Dims) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,CellMap(1,absy+1,absx),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif
!i,j-1
      if(absy.ne.1) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,CellMap(1,absy-1,absx),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif
!     i+1,j+1
      if((absx.ne.Dims).and.(absy.ne.Dims)) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,
     $        CellMap(1,absy+1,absx+1),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif 
!i+1,j-1
      if((absx.ne.Dims).and.(absy.ne.1)) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,
     $        CellMap(1,absy-1,absx+1),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif
!     i-1,j+1
      if((absx.ne.1).and.(absy.ne.Dims)) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,
     $        CellMap(1,absy+1,absx-1),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif
!i-1,j-1
      if((absx.ne.1).and.(absy.ne.1)) then
         CALL MPI_ISend(PartCs,2*Max+1,MPI_REAL8,
     $        CellMap(1,absy-1,absx-1),
     $        tag,MPI_COMM_WORLD,ReqS(iscount),IERROR)
         iscount = iscount + 1
      endif
      nreqr = ircount-1
      nreqs = iscount-1
      DEALLOCATE(PartCs)
      RETURN
      END
      
      SUBROUTINE POSTWAITS(REQR,REQS,NREQR,NREQS)
      IMPLICIT NONE
      include 'mpif.h'
      INTEGER REQR(8),REQS(8),NREQR,NREQS
      INTEGER stats(MPI_STATUS_SIZE,8),IERROR
      CALL MPI_WAITALL(NREQS,REQS,STATS,IERROR)
      CALL MPI_WAITALL(NREQR,REQR,STATS,IERROR)
      RETURN
      END