diff --git a/src/auxval.f90 b/src/auxval.f90
index fd03d9c..990c31e 100644
--- a/src/auxval.f90
+++ b/src/auxval.f90
@@ -1,108 +1,116 @@
 SUBROUTINE auxval
 !
   USE constants
   USE basic
   USE fields
   USE beam
 !
 !   Set auxiliary values
 !
   IMPLICIT NONE
 !
 
   INTEGER:: i
   
 !________________________________________________________________________________
 IF(mpirank .eq. 0) WRITE(*,'(a/)') '=== Set auxiliary values ==='
 !________________________________________________________________________________
 !
 ! Total number of intervals
   nr=sum(nnr)
 
 ! Normalisation constants
   qsim=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)*n0*elchar/nplasma 
   msim=abs(qsim)/elchar*me
   vnorm=vlight
   tnorm=1/sqrt(abs(n0)*elchar*elchar/me/eps_0)
   rnorm=vnorm*tnorm 
   bnorm=B0
   enorm=vlight*bnorm
   phinorm=enorm*rnorm
 
   ! Normalised boundary conditions
   potinn=potinn/phinorm     
   potout=potout/phinorm
 
 ! Characteristic frequencies and normalised volume 
   omegac=qsim/msim*B0
   omegap=sqrt(elchar**2*abs(n0)/me/eps_0)
   V=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)/rnorm/rnorm/rnorm
   IF(mpirank .eq. 0) THEN
     IF(omegap.GT.omegac) THEN
        WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegap/omegac', omegap, omegac, omegap/omegac
     ELSE
        WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegac/omegap', omegap, omegac, omegac/omegap
     END IF
   END IF
   ! Construction of the mesh rgrid in r and zgrid in z and its normalisation
   CALL mesh
   rgrid=rgrid/rnorm
   zgrid=zgrid/rnorm
   dz=dz/rnorm
   dr=dr/rnorm
 
 
   ! Define Zindex bounds for the MPI process in the case of sequential run
   ALLOCATE( Zbounds(0:mpisize), Nplocs_all(0:mpisize-1))
   Zbounds(0) = Zgrid(0)
   Zbounds(mpisize) = Zgrid(nz)
   IF(mpisize.gt.1) THEN
     leftproc=mpirank-1
-    IF(mpirank .eq. 0) leftproc = mpisize-1
+    IF(mpirank .eq. 0 .AND. partperiodic) THEN 
+      leftproc = mpisize-1
+    ELSE IF (mpirank .eq. 0) THEN
+      leftproc=-1
+    END IF
     rightproc = mpirank+1
-    IF(mpirank .eq. mpisize-1) rightproc = 0
+    IF(mpirank .eq. mpisize-1 .AND. partperiodic) THEN 
+      rightproc = 0
+    ELSE IF (mpirank .eq. mpisize-1) THEN 
+      rightproc=-1
+    END IF
   ELSE
     leftproc=-1
     rightproc=-1
   END IF
   WRITE(*,*) "mpirank, left, right", mpirank, leftproc, rightproc
 
 ! Auxiliary vectors
   ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1)))
   DO i=0,nr
      vec1(i*(nz+1)+1:(i+1)*(nz+1))=zgrid!(0:nz)
      vec2(i*(nz+1)+1:(i+1)*(nz+1))=rgrid(i)
   END DO
   ALLOCATE(partdensity(nz*nr))
   ALLOCATE(fluidur(nz*nr))
   ALLOCATE(fluiduz(nz*nr))
   ALLOCATE(fluiduthet(nz*nr))
   ALLOCATE(Presstens(6,nz*nr))
 
 !________________________________________________________________________________
 CONTAINS 
 !________________________________________________________________________________
   SUBROUTINE mesh
     USE basic, ONLY: zgrid, rgrid, nr, nnr, nz, dr, dz, lz, radii
    INTEGER :: j
     ALLOCATE(zgrid(0:nz),rgrid(0:nr))
     dz=(lz(2)-lz(1))/nz
     DO j=0,nz
        zgrid(j)=j*dz+lz(1)
     END DO
     dr(1)=(radii(2)-radii(1))/nnr(1)
     DO j=0,nnr(1)
        rgrid(j)=radii(1)+j*dr(1)
     END DO
     dr(2)=(radii(3)-radii(2))/nnr(2)
     DO j=1,nnr(2)
        rgrid(nnr(1)+j)=radii(2)+j*dr(2)
     END DO
     dr(3)=(radii(4)-radii(3))/nnr(3)
     DO j=1,nnr(3)
        rgrid(nnr(1)+nnr(2)+j)=radii(3)+j*dr(3)
     END DO
   END SUBROUTINE mesh
 !________________________________________________________________________________
 
 END SUBROUTINE auxval
diff --git a/src/basic_mod.f90 b/src/basic_mod.f90
index ed1b91b..3280d35 100644
--- a/src/basic_mod.f90
+++ b/src/basic_mod.f90
@@ -1,306 +1,308 @@
 MODULE basic
 !
   USE hashtable
   USE constants
   USE bsplines
   USE mumps_bsplines
   USE futils
   USE mpihelper
   IMPLICIT NONE
 !
 !   Basic module for time dependent problems
 !
   CHARACTER(len=80) :: label1, label2, label3, label4
 !
 !   BASIC Namelist
 !
   LOGICAL          :: nlres = .FALSE.    ! Restart flag
   LOGICAL          :: nlsave = .TRUE.    ! Checkpoint (save) flag
   LOGICAL          :: newres=.FALSE.   ! New result HDF5 file
   LOGICAL          :: nlxg=.FALSE.      ! Show graphical interface Xgrafix
   INTEGER          :: nrun=1           ! Number of time steps to run
   DOUBLE PRECISION    :: job_time=3600.0  ! Time allocated to this job in seconds
   DOUBLE PRECISION    :: tmax=100000.0    ! Maximum simulation time
   DOUBLE PRECISION    :: extra_time=60.0  ! Extra time allocated
   DOUBLE PRECISION    :: dt=1             ! Time step
   DOUBLE PRECISION    :: time=0           ! Current simulation time (Init from restart file)
 !
 !   Other basic global vars and arrays
 !
   INTEGER :: jobnum                    ! Job number
   INTEGER :: step                      ! Calculation step of this run
   INTEGER :: cstep=0                   ! Current step number (Init from restart file)
   LOGICAL :: nlend                     ! Signal end of run
   INTEGER :: ierr                      ! Integer used for MPI
   INTEGER :: it0d=1, it2d=1, itparts=1 ! Number of iterations between 0d, 2d, particles, values writes to hdf5
   INTEGER :: ittext=1                  ! Number of iterations between text outputs in the console
   INTEGER :: itrestart=10000           ! Number of iterations between save of restart.h5 file
   INTEGER :: itgraph                   ! Number of iterations between graphical interface updates
   INTEGER :: mpirank                   ! MPIrank of the current processus
   INTEGER :: mpisize                   ! Size of the MPI_COMM_WORLD communicator
   INTEGER :: rightproc, leftproc        ! Rank of next and previous processor in the z decomposition
 !
 !  List of logical file units
   INTEGER :: lu_in   = 90              ! File duplicated from STDIN
   INTEGER :: lu_stop = 91              ! stop file, see subroutine TESEND
 !
 !  HDF5 file
   CHARACTER(len=64) :: resfile = "results.h5" ! Main result file
   CHARACTER(len=64) :: rstfile = "restart.h5" ! Restart file
   INTEGER           :: fidres                 ! FID for resfile
   INTEGER           :: fidrst                 ! FID for restart file
   TYPE(BUFFER_TYPE) :: hbuf0                  ! Hashtable for 0d var
 ! 
 ! Plasma parameters
-  LOGICAL           :: nlPhis=.TRUE.          ! Calculate self consistent electric field flag
-  LOGICAL           :: nlclassical=.FALSE.    ! If true, solves the equation of motion classicaly
+  LOGICAL           :: nlPhis= .TRUE.          ! Calculate self consistent electric field flag
+  LOGICAL           :: nlclassical= .FALSE.    ! If true, solves the equation of motion classicaly
   LOGICAL           :: nlperiod(2)=(/.false.,.false./)        ! Set periodic splines on or off
+  LOGICAL           :: partperiodic= .TRUE.   ! Sets if the particles boundary conditions are periodic or open
   INTEGER           :: nplasma                ! Number of superparticles
   INTEGER           :: nblock                 ! Number of intervals in Z for stable distribution initialisation
   DOUBLE PRECISION     :: potinn,potout          ! Potential at inner and outer metalic walls
   DOUBLE PRECISION     :: B0                     ! Max magnitude of magnetic field
   DOUBLE PRECISION, allocatable     :: Bz(:), Br(:)   ! Magnetic field components
   DOUBLE PRECISION, allocatable     :: Athet(:)       ! Magnetic potential
   TYPE(spline2d)    :: splrz                  ! Spline at r and z
   DOUBLE PRECISION, allocatable     :: Ez(:), Er(:)   ! Electric field components
   DOUBLE PRECISION, allocatable     :: pot(:)    ! Electro static potential
   DOUBLE PRECISION     :: radii(4)     ! Inner and outer radius of cylinder and radii of fine mesh region
   DOUBLE PRECISION     :: plasmadim(4)         ! Zmin Zmax Rmin Rmax values for plasma particle loading
   INTEGER           :: distribtype=1               ! Type of distribution function used to load the particles
                                                  !1: gaussian, 2: Stable as defined in 4.95 of Davidson
   DOUBLE PRECISION     :: H0=0, P0=0             ! Initial value of Hamiltonian and canonical angular momentum
                                               ! for distribution 2 
   DOUBLE PRECISION     :: lz(2)                     ! Length of cylinder in z direction
   DOUBLE PRECISION     :: n0                     ! Average density of physical plasma
   DOUBLE PRECISION, ALLOCATABLE:: partdensity(:) ! Number of particles per gridcell computed every it2d
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: fluidur, fluiduz, fluiduthet ! Fluid velocities per gridcell computed every it2d
   DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: Presstens ! Pressure tensor per gridcell computed every it2d
   INTEGER           :: nz, nr, nnr(3)           ! Number of intervals in z and r
   DOUBLE PRECISION     :: dz, dr(3)                ! Cell size in z and r for each region
   DOUBLE PRECISION, allocatable     :: zgrid(:), rgrid(:)   ! Nodes positions
   DOUBLE PRECISION     :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm     ! Normalization constants
   DOUBLE PRECISION     :: qsim                     ! Charge of superparticles
   DOUBLE PRECISION     :: msim                     ! Mass of superparticles
   INTEGER           :: femorder(2)              ! FEM order 
   INTEGER           :: ngauss(2)                ! Number of gauss points 
   LOGICAL           :: nlppform =.TRUE.         ! Argument of set_spline
   INTEGER, SAVE        :: nrank(2)               ! Number of splines in both directions
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: rhs        ! right hand side of the poisson equation solver
   DOUBLE PRECISION     :: omegac                   ! Cyclotronic frequency
   DOUBLE PRECISION     :: omegap                   ! Plasma frequency
   DOUBLE PRECISION     :: V                        ! Normalised volume
   DOUBLE PRECISION     :: temp                     ! Initial temperature of plasma
   DOUBLE PRECISION     :: Rcurv                    ! Magnetic field curvature coefficient
   DOUBLE PRECISION     :: Width                    ! Distance between two magnetic mirrors
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Zbounds              ! Index of bounds for local processus in Z direction
 
   TYPE(BASICDATA) :: bdata
   
 CONTAINS
 !
 !================================================================================
   SUBROUTINE basic_data
 !
 !   Define basic data
 !
     use mpihelper
     IMPLICIT NONE
 !
 !   Local vars and arrays
     CHARACTER(len=256) :: line, inputfilename
 !
     NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg,          &
          &           nplasma, potinn, potout, B0, lz, n0, nz, nnr, femorder, ngauss,             &
          &           nlppform, plasmadim, radii, temp, Rcurv, width, it0d, it2d, itparts, ittext, &
-         &           resfile, rstfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0
+         &           resfile, rstfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0, partperiodic
 !________________________________________________________________________________
 !                   1.   Process Standard Input File
 !
   IF(mpirank .eq. 0) THEN
   IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN
     WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENT REQUIRED, STOPPING'
     STOP
   ENDIF
   CALL GET_COMMAND_ARGUMENT(1,inputfilename)
 
   OPEN(UNIT=lu_in,FILE=trim(inputfilename))
     !DO
     !   READ(*,'(a)', END=110) line
     !   WRITE(lu_in, '(a)') TRIM(line)
     !END DO
 !110 CONTINUE
     !REWIND(lu_in)
 !________________________________________________________________________________
 !                   1.   Label the run
 !
     READ(lu_in,'(a)') label1
     READ(lu_in,'(a)') label2
     READ(lu_in,'(a)') label3
     READ(lu_in,'(a)') label4
 !
     WRITE(*,'(12x,a/)') label1(1:len_trim(label1))
     WRITE(*,'(12x,a/)') label2(1:len_trim(label2))
     WRITE(*,'(12x,a/)') label3(1:len_trim(label3))
     WRITE(*,'(12x,a/)') label4(1:len_trim(label4))
 !________________________________________________________________________________
 !                   2.   Read in basic data specific to run
 !
     READ(lu_in,basic)
     WRITE(*,basic)
     bdata=BASICDATA( nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical, &
       & nplasma, nz, it0d, it2d, itparts, itgraph, distribtype, nblock,     &
       & nrun, job_time, extra_time, tmax, dt, potinn, potout, B0, n0, temp, &
-      & Rcurv, width, H0, P0, femorder, ngauss, nnr, lz, radii, plasmadim, resfile )
+      & Rcurv, width, H0, P0, femorder, ngauss, nnr, lz, radii, plasmadim, resfile, partperiodic)
   END IF
   CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation
   CALL MPI_Bcast(bdata, 1, basicdata_type, 0, MPI_COMM_WORLD, ierr)
   CALL readbdata
 
 !
   END SUBROUTINE basic_data
   !================================================================================
   SUBROUTINE daytim(str)
     !
     !   Print date and time
     !
         IMPLICIT NONE
     !
         CHARACTER(len=*), INTENT(in) :: str
     !
     !   Local vars and arrays
         CHARACTER(len=16) :: d, t, dat, functime
     !________________________________________________________________________________
     !
         CALL DATE_AND_TIME(d,t)
         dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4)
         functime=t(1:2) // ':' // t(3:4) // ':' // t(5:10)
         WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), functime(1:12)
     !
       END SUBROUTINE daytim
 !================================================================================
   SUBROUTINE timera(cntrl, str, eltime)
 !
 !   Timers (cntrl=0/1 to Init/Update)
 !
     IMPLICIT NONE
     INTEGER, INTENT(in)                  :: cntrl
     CHARACTER(len=*), INTENT(in)         :: str
     DOUBLE PRECISION, OPTIONAL, INTENT(out) :: eltime
 !
     INTEGER, PARAMETER                    :: ncmax=128
     INTEGER, SAVE                         :: icall=0, nc=0
     DOUBLE PRECISION, SAVE                   :: startt0=0.0
     CHARACTER(len=16), SAVE               :: which(ncmax)
     INTEGER                               :: lstr, found, i
     DOUBLE PRECISION                         :: seconds
     DOUBLE PRECISION, DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0
 !________________________________________________________________________________
     IF( icall .EQ. 0 ) THEN
        icall = icall+1
        startt0 = seconds()
     END IF
 
     lstr = LEN_TRIM(str)
     IF( lstr .GT. 0 ) found = loc(str)
 !________________________________________________________________________________
 !
     SELECT CASE (cntrl)
 !
     CASE(-1)    !  Current wall time
        IF( PRESENT(eltime) ) THEN
           eltime = seconds() - startt0
        ELSE
           WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', seconds() - startt0
        END IF
 !
     CASE(0)    !  Init Timer
        IF( found .EQ. 0 ) THEN  !  Called for the 1st time for 'str'
           nc = nc+1
           which(nc) = str(1:lstr)
           found = nc
        END IF
        startt(found) = seconds()
 !
     CASE(1)   !  Update timer
        endt(found) = seconds() - startt(found)
        IF( PRESENT(eltime) ) THEN
           eltime = endt(found)
        ELSE
           WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', endt(found)
        END IF
 !
     CASE(2)   !  Update and reset timer
        endt(found) = endt(found) + seconds() - startt(found)
        startt(found) = seconds()
        IF( PRESENT(eltime) ) THEN
           eltime = endt(found)
        END IF
 !
     CASE(9)   !  Display all timers
        IF( nc .GT. 0 ) THEN
           WRITE(*,'(a)') "Timer Summary"
           WRITE(*,'(a)') "============="
           DO i=1,nc
              WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endt(i)
           END DO
        END IF
 !
     END SELECT
 !
   CONTAINS
     INTEGER FUNCTION loc(str)
       CHARACTER(len=*), INTENT(in) :: str
       INTEGER :: i, ind
       loc = 0
       DO i=1,nc
           ind = INDEX(which(i), str(1:lstr))
           IF( ind .GT. 0 .AND. LEN_TRIM(which(i)) .EQ. lstr ) THEN
              loc = i
              EXIT
           END IF
        END DO
     END FUNCTION loc
   END SUBROUTINE timera
 !================================================================================
 SUBROUTINE readbdata
 nlres       = bdata%nlres         
 nlsave      = bdata%nlsave
 newres      = bdata%newres
 nlxg       = bdata%nlxg
 nlppform    = bdata%nlppform
 nlPhis      = bdata%nlPhis
 nlclassical = bdata%nlclassical 
 nplasma     = bdata%nplasma
 nz       = bdata%nz
 it0d        = bdata%it0d 
 it2d        = bdata%it2d
 itparts     = bdata%itparts
 itgraph     = bdata%itgraph
 distribtype = bdata%distribtype
 nblock      = bdata%nblock
 nrun        = bdata%nrun
 job_time    = bdata%job_time
 extra_time  = bdata%extra_time
 tmax        = bdata%tmax
 dt        = bdata%dt
 potinn      = bdata%potinn  
 potout      = bdata%potout
 B0        = bdata%B0
 n0         = bdata%n0
 temp        = bdata%temp
 Rcurv       = bdata%Rcurv
 width       = bdata%width
 H0       = bdata%H0
 P0         = bdata%P0
 femorder    = bdata%femorder    
 ngauss      = bdata%ngauss
 nnr      = bdata%nnr
 lz       = bdata%lz
 radii       = bdata%radii
 plasmadim   = bdata%plasmadim
 resfile     = bdata%resfile
+partperiodic= bdata%partperiodic
 
 END SUBROUTINE readbdata
 
 !================================================================================
 END MODULE basic
diff --git a/src/beam_mod.f90 b/src/beam_mod.f90
index 1162483..65a54d4 100644
--- a/src/beam_mod.f90
+++ b/src/beam_mod.f90
@@ -1,1323 +1,1329 @@
 !------------------------------------------------------------------------------
 ! EPFL/Swiss Plasma Center
 !------------------------------------------------------------------------------
 !
 ! MODULE: beam
 !
 !> @author
 !> Guillaume Le Bars EPFL/SPC
 !
 ! DESCRIPTION: 
 !> Module responsible for loading, advancing and computing the necessary diagnostics for the simulated particles. 
 !
 !------------------------------------------------------------------------------
 
 MODULE beam
 !
   USE constants
   USE mpi
   USE mpihelper
   USE basic, ONLY: mpirank, mpisize
   
   IMPLICIT NONE
 
   !! Stores the particles properties for the run.
   TYPE particles
   INTEGER :: Nptot !! Local number of simulated particles 
   INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex !! Index in the electric potential grid for the R direction
   INTEGER, DIMENSION(:), ALLOCATABLE :: Zindex !! Index in the electric potential grid for the Z direction
   INTEGER, DIMENSION(:), ALLOCATABLE :: partindex !! Index of the particle to be able to follow it when it goes from one MPI host to the other
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
       &   R,Z,THET,   &
       &   WZ, WR, &
       &   pot, Er, Ez
   DOUBLE PRECISION, DIMENSION(:),POINTER:: &
       &   UR, URold, &
       &   UTHET, UTHETold, &
       &   UZ, UZold, &
       &   Gamma, Gammaold
   END TYPE particles
 
 !
   TYPE(particles) :: parts  ! Storage for all the particles
   SAVE :: parts
   TYPE(particle), DIMENSION(:), ALLOCATABLE:: rrecvpartbuff, lrecvpartbuff, rsendpartbuff, lsendpartbuff ! buffer to send and receive particle from left and right processes
   
 !   Diagnostics (scalars)
   DOUBLE PRECISION :: ekin ! Total kinetic energz (J)
   DOUBLE PRECISION :: epot  ! Total potential energy (J)
   DOUBLE PRECISION :: etot, etot0  ! Current and Initial total energy (J)
 
   INTEGER, DIMENSION(7), SAVE :: ireducerequest=MPI_REQUEST_NULL
 !
  INTEGER, DIMENSION(:), ALLOCATABLE :: Nplocs_all !F.B. get local numbers of particles
  LOGICAL:: collected=.false.  !Stores if the particles data have been collected to root process this timestep
 !
 CONTAINS
 
 !---------------------------------------------------------------------------  
    !> @author 
    !> Guillaume Le Bars EPFL/SPC
    !
    ! DESCRIPTION: 
    !> 
    !> @brief
    !> Allocate the memory for the particles variable storing the particles quantities.
    !
    !> @param[inout] p the particles variable needing to be allocated. 
    !> @param[in] nparts the maximum number of particles that will be stored in this variable      
    !---------------------------------------------------------------------------  
    
 
 SUBROUTINE creat_parts(p, nparts)
   TYPE(particles)       :: p
   INTEGER, INTENT(in)   :: nparts
   
   IF (.NOT. ALLOCATED(p%Z) ) THEN
     p%Nptot = nparts
     ALLOCATE(p%Z(nparts))
     ALLOCATE(p%R(nparts))
     ALLOCATE(p%THET(nparts))
     ALLOCATE(p%WZ(nparts))
     ALLOCATE(p%WR(nparts))
     ALLOCATE(p%UR(nparts))
     ALLOCATE(p%UZ(nparts))
     ALLOCATE(p%UTHET(nparts))
     ALLOCATE(p%URold(nparts))
     ALLOCATE(p%UZold(nparts))
     ALLOCATE(p%UTHETold(nparts))
     ALLOCATE(p%Gamma(nparts))
     ALLOCATE(p%Rindex(nparts))
     ALLOCATE(p%Zindex(nparts))
     ALLOCATE(p%partindex(nparts))
     ALLOCATE(p%pot(nparts))
     ALLOCATE(p%Er(nparts))
     ALLOCATE(p%Ez(nparts))
     ALLOCATE(p%GAMMAold(nparts))
     parts%URold=0
     parts%UZold=0
     parts%UTHETold=0
     parts%rindex=0
     parts%zindex=0
     parts%WR=0
     parts%WZ=0
     parts%UR=0
     parts%UZ=0
     parts%UTHET=0
     parts%Z=0
     parts%R=0
     parts%THET=0
     parts%Gamma=1
     parts%Er=0
     parts%Ez=0
     parts%pot=0
     parts%gammaold=1
   END IF
 END SUBROUTINE creat_parts
 !______________________________________________________________________________
 SUBROUTINE load_parts
     ! Loads the particles at the beginning of the simulation and create the parts variable if necessary
     USE basic, ONLY: nplasma, mpirank, ierr, distribtype, mpisize, nlclassical, Rcurv
     USE mpi
 
     INTEGER:: i
     DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
 
     ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma))
 
     ! Select case to define the type of distribution
     SELECT CASE(distribtype)
       CASE(1) ! Gaussian distribution in V and uniform in R
         CALL loaduniformRZ(VR, VZ, VTHET)
       CASE(2) !Stable distribution from Davidson 4.95 p.119
         CALL loadDavidson(VR, VZ, VTHET, lodunir)
       CASE(3) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as 1/R**2
         CALL loadDavidson(VR, VZ, VTHET, lodinvr)
       CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R 
         CALL loadDavidson(VR, VZ, VTHET, lodgausr)
       CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V computed from v_th given by temp
         CALL loadDavidson(VR, VZ, VTHET, lodunir)
       CASE(6) ! Uniform distribution in R and Z and Gaussian distribution in V with Vz<V_perp to satisfy magnetic mirror trapping
         CALL loaduniformRZ(VR, VZ, VTHET)
         VZ = sqrt(2*(Rcurv-1)/(Rcurv+1))*0.1*VZ
       CASE DEFAULT
         IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of distribution:", distribtype
         CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
 
     END SELECT
     Do i=1,nplasma
       parts%partindex(i)=i
     END DO
 
     IF(nlclassical) THEN
       parts%Gamma=1
     ELSE
       parts%Gamma=sqrt(1/(1-VR**2+VZ**2+VTHET**2))
     END IF
     parts%UR=parts%Gamma*VR
     parts%UZ=parts%Gamma*VZ
     parts%UTHET=parts%Gamma*VTHET
     DEALLOCATE(VZ, VR, VTHET)
 
     CALL localisation
     IF(mpisize .gt. 1) THEN
       CALL distribpartsonprocs
     END IF
     
 END SUBROUTINE load_parts
 !________________________________________________________________________________
 SUBROUTINE bound
 
-  USE basic, ONLY: zgrid, nz, Zbounds, cstep, mpirank
+  USE basic, ONLY: zgrid, nz, Zbounds, cstep, mpirank, partperiodic
   INTEGER :: i, rsendnbparts=0, lsendnbparts=0
   INTEGER, DIMENSION(parts%Nptot) :: sendhole
   sendhole=0
   lsendnbparts=0
   rsendnbparts=0
   IF (parts%Nptot .NE. 0) THEN
   ! Boundary condition at z direction
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
     DO i=1,parts%Nptot
       IF(cstep .ne. 0) THEN
         IF (parts%Z(i) .ge. Zbounds(mpirank+1)) THEN
         !$OMP CRITICAL (nbparts)
           rsendnbparts=rsendnbparts+1
           sendhole(lsendnbparts+rsendnbparts)=i
         !$OMP END CRITICAL (nbparts)
         ELSE IF (parts%Z(i) .lt. Zbounds(mpirank)) THEN
         !$OMP CRITICAL (nbparts)
           lsendnbparts=lsendnbparts+1
           sendhole(lsendnbparts+rsendnbparts)=-i
         !$OMP END CRITICAL (nbparts)
         END IF
       END IF
       DO WHILE (parts%Z(i) .GT. zgrid(nz))
          parts%Z(i) = parts%Z(i) - zgrid(nz) + zgrid(0)
       END DO
       DO WHILE (parts%Z(i) .LT. zgrid(0))
          parts%Z(i) = parts%Z(i) + zgrid(nz) - zgrid(0)
       END DO
     END DO
   !$OMP END PARALLEL DO
     IF(mpisize .gt. 1) THEN
       CALL particlescommunication(lsendnbparts, rsendnbparts, sendhole)
+    ELSE IF (.NOT. partperiodic) THEN
+      IF(lsendnbparts+rsendnbparts .gt. 0) THEN
+        DO i=rsendnbparts+lsendnbparts,1,-1
+            CALL move_part(parts%Nptot, abs(sendhole(i)))
+            parts%Nptot=parts%Nptot-1
+        END DO
+      END IF
     END IF
   END IF
 
 END subroutine bound
 !________________________________________________________________________________
 SUBROUTINE localisation
   USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm
   INTEGER :: j,i, nblostparts=0
   INTEGER, DIMENSION(parts%Nptot) :: losthole
   i=1
   losthole=0
 
   IF (parts%Nptot .NE. 0) THEN
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
     DO i=1,parts%Nptot
        parts%zindex(i)=(parts%Z(i)-zgrid(0))/dz
        IF (parts%R(i) .GT. rgrid(0) .AND.  parts%R(i) .LT. rgrid(nnr(1))) THEN
           parts%rindex(i)=(parts%R(i)-rgrid(0))/dr(1)
        ELSE IF(parts%R(i) .GT. rgrid(nnr(1)) .AND.  parts%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN
           parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)))/dr(2)+nnr(1)
        ELSE IF(parts%R(i) .GT. rgrid(nnr(1)+nnr(2)) .AND.  parts%R(i) .LT. rgrid(nr)) THEN
           parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)+nnr(2)))/dr(3)+nnr(1)+nnr(2)
        ELSE
           WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in r:", parts%R(i)*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm 
           !nlend=.TRUE.
           WRITE(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Particle removed  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
           !$OMP CRITICAL (lostparts)
             nblostparts=nblostparts+1
             losthole(nblostparts)=i
           !$OMP END CRITICAL (lostparts)
        END IF
     END DO
   !$OMP END PARALLEL DO
 
     IF(nblostparts.gt.0) THEN
       DO j=1,nblostparts
           CALL move_part(parts%Nptot, losthole(j))
           parts%Nptot=parts%Nptot-1
       END DO
     END IF
 
     !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
     DO i=1,parts%Nptot
       parts%WZ(i)=(parts%Z(i)-zgrid(parts%zindex(i)))/(zgrid(parts%zindex(i)+1)-zgrid(parts%zindex(i)));
       parts%WR(i)=(parts%R(i)-rgrid(parts%rindex(i)))/(rgrid(parts%rindex(i)+1)-rgrid(parts%rindex(i)));
     END DO
     !$OMP END PARALLEL DO
      
   END IF
 END SUBROUTINE localisation
 !________________________________________________________________________________________________________
 SUBROUTINE comp_velocity
 !
 !   Computes the new velocity of the particles due to Lorentz force
 !
   USE basic, ONLY : nlclassical
   ! Store old Velocities
   CALL swappointer(parts%UZold, parts%UZ)  
   CALL swappointer(parts%URold, parts%UR)  
   CALL swappointer(parts%UTHETold, parts%UTHET)  
   CALL swappointer(parts%Gammaold, parts%Gamma)  
 
   IF (nlclassical) THEN
     CALL comp_velocity_fun(gamma_classical)
   ELSE
     CALL comp_velocity_fun(gamma_relativistic)
   END IF 
 
 END SUBROUTINE comp_velocity
 !____________________________________________
 SUBROUTINE comp_velocity_fun(gamma)
 !
 !   Computes the new velocity of the particles due to Lorentz force
 !
   USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz
   interface
     subroutine gamma(gam, UZ, UR, UTHET)
       DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET
       DOUBLE PRECISION, INTENT(OUT):: gam
     end subroutine
   end interface
   DOUBLE PRECISION :: tau
   DOUBLE PRECISION:: BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2
   INTEGER:: J1, J2, J3, J4
   INTEGER:: i
  ! Normalized time increment
   tau=omegac/2/omegap*dt/tnorm
   IF (parts%Nptot .NE. 0) THEN
     !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
     DO i=1,parts%Nptot
        J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1
        J2=J1+1
        J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1
        J4=J3+1
     ! Interpolation for magnetic field
        BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) &
        &   +parts%WZ(i)*(1-parts%WR(i))*Bz(J3)    &
        &   +(1-parts%WZ(i))*parts%WR(i)*Bz(J2)    &
        &   +parts%WZ(i)*parts%WR(i)*Bz(J1)
        BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) &
        &   +parts%WZ(i)*(1-parts%WR(i))*Br(J3)    &
        &   +(1-parts%WZ(i))*parts%WR(i)*Br(J2)    &
        &   +parts%WZ(i)*parts%WR(i)*Br(J1)
     ! First half of electric pulse
        parts%UZ(i)=parts%UZold(i)+parts%Ez(i)*tau
        parts%UR(i)=parts%URold(i)+parts%ER(i)*tau
        CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i))
     ! Rotation along magnetic field
        ZBZ=tau*BRZ/parts%Gamma(i)
        ZBR=tau*BRR/parts%Gamma(i)
        ZPZ=parts%UZ(i)-ZBR*parts%UTHETold(i)                     !u'_{z} 
        ZPR=parts%UR(i)+ZBZ*parts%UTHETold(i)                     !u'_{r}
        ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZ(i)-ZBZ*parts%UR(i))          !u'_{theta}
        SQR=1+ZBZ*ZBZ+ZBR*ZBR
        ZBZ2=2*ZBZ/SQR
        ZBR2=2*ZBR/SQR
        parts%UZ(i)=parts%UZ(i)-ZBR2*ZPTHET                    !u+_{z}
        parts%UR(i)=parts%UR(i)+ZBZ2*ZPTHET                   !u+_{r}
        parts%UTHET(i)=parts%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR)      !u+_{theta}
     
     ! Second half of acceleration
        parts%UZ(i)=parts%UZ(i)+parts%EZ(i)*tau
        parts%UR(i)=parts%UR(i)+parts%ER(i)*tau
 
        CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i)) 
     END DO
     !$OMP END PARALLEL DO
   END IF
   collected=.false.
 END SUBROUTINE comp_velocity_fun
 !______________________________________________________________________________
 SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET)
       DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET
       DOUBLE PRECISION, INTENT(OUT):: gamma
       gamma=1.0
 END SUBROUTINE gamma_classical
 !______________________________________________________________________________
 SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET)
       DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET
       DOUBLE PRECISION, INTENT(OUT):: gamma
       gamma=sqrt(1+UZ**2+UR**2+UTHET**2)
 END SUBROUTINE
 !________________________________________________________________________________
 SUBROUTINE push
     Use basic, ONLY: dt, tnorm
     DOUBLE PRECISION:: XP, YP, COSA, SINA, U1, U2, ALPHA
     INTEGER :: i
 
     IF (parts%Nptot .NE. 0) THEN 
     !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i, XP, YP, COSA, SINA, U1, U2, ALPHA)
     DO i=1,parts%Nptot 
 ! Local Cartesian coordinates
           XP=parts%R(i)+dt/tnorm*parts%UR(i)/parts%Gamma(i)
           YP=dt/tnorm*parts%UTHET(i)/parts%Gamma(i)
 
 ! Conversion to cylindrical coordiantes
           parts%Z(i)=parts%Z(i)+dt/tnorm*parts%UZ(i)/parts%Gamma(i)
           parts%R(i)=sqrt(XP**2+YP**2)
 ! Velocity in rotated reference frame
           IF (parts%R(i) .EQ. 0) THEN 
             COSA=1
             SINA=0
             ALPHA=0
           ELSE
             COSA=XP/parts%R(i)
             SINA=YP/parts%R(i)
             ALPHA=asin(SINA)
           END IF
           parts%THET(i)=parts%THET(i)+ALPHA
 
           U1=COSA*parts%UR(i)+SINA*parts%UTHET(i)
           U2=-SINA*parts%UR(i)+COSA*parts%UTHET(i)
 
           parts%UR(i)=U1
           parts%UTHET(i)=U2 
 
         END DO
         !$OMP END PARALLEL DO
     END IF
     collected=.false.
 END SUBROUTINE push
 !_______________________________________________________________________________
 !_______________________________________________________________________________
 
 SUBROUTINE diagnostics
 !
 !  Compute energies
 !
 !!$    CALL energies
     USE constants, ONLY: vlight
-    USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, partdensity, ierr, step, it2d, nlend, nr, itparts, nplasma, fluidur, fluiduz, fluiduthet, Presstens
+    USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, partdensity, ierr, step, it2d, nlend, nr, itparts, nplasma, fluidur, fluiduz, fluiduthet, Presstens, partperiodic
     DOUBLE PRECISION :: gammamid, vr, vz, vthet
     INTEGER :: i, gridindex
     INTEGER, DIMENSION(:) :: stat(7*MPI_STATUS_SIZE)
     CALL MPI_WAITALL(5, ireducerequest(3:7), stat(2*MPI_STATUS_SIZE+1:7*MPI_STATUS_SIZE), ierr)
     partdensity=0
     fluidur=0
     fluiduz=0
     fluiduthet=0
     ekin=0
     epot=0
     etot=0
 
 
     !! Computation of the kinetic and potential energy as well as  fluid velocities and density
     !$OMP PARALLEL DO REDUCTION(+:epot, ekin, partdensity,fluidur,fluiduz,fluiduthet) DEFAULT(SHARED) PRIVATE(i,gammamid,vr,vz,vthet,gridindex)
     DO i=1,parts%Nptot
-!      
+!      Potential energy
        epot=epot+0.5*qsim*parts%pot(i)*phinorm
+       ! Kinetic energy
        IF(.not. nlclassical) THEN
         gammamid=1/sqrt(1-(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))/(parts%Gammaold(i)*parts%Gamma(i)))
         ekin=ekin+msim*vlight**2*(gammamid-1)
        ELSE
         ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))
        END IF
-
+      
        IF(modulo(step,it2d).eq. 0 .or. nlend) THEN
           gridindex=parts%zindex(i)+1+parts%rindex(i)*nz
+          ! Number of particles per cell
           partdensity(gridindex)=partdensity(gridindex)+1
+          ! Individual particle velocities
           vr=parts%UR(i)/parts%Gamma(i)
           vz=parts%UZ(i)/parts%Gamma(i)
           vthet=parts%UTHET(i)/parts%Gamma(i)
+          ! Fluid velocities
           fluidur(gridindex)=fluidur(gridindex)+vr
           fluiduz(gridindex)=fluiduz(gridindex)+vz
           fluiduthet(gridindex)=fluiduthet(gridindex)+vthet
        END IF
 
     END DO
     !$OMP END PARALLEL DO
 
     IF(mpirank .eq.0 ) THEN
         CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
         & 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
         CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
         & 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
     ELSE 
         CALL MPI_IREDUCE(epot, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
         & 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
         CALL MPI_IREDUCE(ekin, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
         & 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
     END IF
 
-    WHERE (partdensity.gt.0) 
+    WHERE (partdensity.gt.0)
+      ! Normalize fluid velocities by density 
       fluiduz=fluiduz/partdensity
       fluidur=fluidur/partdensity
       fluiduthet=fluiduthet/partdensity
     END WHERE
 
-    ! Stores the particle density on the grid
-    IF(modulo(step,it2d).eq. 0 .or. nlend) THEN
-      IF(mpirank .eq.0 ) THEN
-        CALL MPI_IREDUCE(MPI_IN_PLACE, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
-        CALL MPI_IREDUCE(MPI_IN_PLACE, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(4), ierr)
-        CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(5), ierr)
-        CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduthet, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(6), ierr)
-      ELSE 
-        CALL MPI_IREDUCE(partdensity, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
-        CALL MPI_IREDUCE(fluidur, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(4), ierr)
-        CALL MPI_IREDUCE(fluiduz, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(5), ierr)
-        CALL MPI_IREDUCE(fluiduthet, fluiduthet, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
-        & 0, MPI_COMM_WORLD, ireducerequest(6), ierr)
-      END IF
-    END IF
-
+    ! Computes the pressure tensor
     IF(modulo(step,it2d).eq. 0 .or. nlend) THEN
       !$OMP PARALLEL DO REDUCTION(+: Presstens) DEFAULT(SHARED) PRIVATE(i,vr,vz,vthet,gridindex)
       DO i=1,parts%Nptot
             gridindex=parts%zindex(i)+1+parts%rindex(i)*nz
             vr=parts%UR(i)/parts%Gamma(i)-fluidur(gridindex)
             vz=parts%UZ(i)/parts%Gamma(i)-fluiduz(gridindex)
             vthet=parts%UTHET(i)/parts%Gamma(i)-fluiduthet(gridindex)
             Presstens(1,gridindex)=Presstens(1,gridindex)+vr*vr
             Presstens(2,gridindex)=Presstens(2,gridindex)+vr*vthet
             Presstens(3,gridindex)=Presstens(3,gridindex)+vr*vz
             Presstens(4,gridindex)=Presstens(4,gridindex)+vthet*vthet
             Presstens(5,gridindex)=Presstens(5,gridindex)+vthet*vz
             Presstens(6,gridindex)=Presstens(6,gridindex)+vz*vz
       END DO
       !$OMP END PARALLEL DO
       !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
       DO i=1,6
             WHERE (partdensity.gt.0) Presstens(i,:)=Presstens(i,:)/partdensity
       END DO
       !$OMP END PARALLEL DO
     END IF
-    
-    ! Stores the particle density on the grid
+
+    ! Send the particle density and fluid velocities from the workers to the master
     IF(modulo(step,it2d).eq. 0 .or. nlend) THEN
       IF(mpirank .eq.0 ) THEN
+        CALL MPI_IREDUCE(MPI_IN_PLACE, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
+        CALL MPI_IREDUCE(MPI_IN_PLACE, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(4), ierr)
+        CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(5), ierr)
+        CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduthet, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(6), ierr)
         CALL MPI_IREDUCE(MPI_IN_PLACE, Presstens, 6*nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
         & 0, MPI_COMM_WORLD, ireducerequest(7), ierr)
       ELSE 
+        CALL MPI_IREDUCE(partdensity, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
+        CALL MPI_IREDUCE(fluidur, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(4), ierr)
+        CALL MPI_IREDUCE(fluiduz, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(5), ierr)
+        CALL MPI_IREDUCE(fluiduthet, fluiduthet, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
+        & 0, MPI_COMM_WORLD, ireducerequest(6), ierr)
         CALL MPI_IREDUCE(Presstens, Presstens, 6*nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
         & 0, MPI_COMM_WORLD, ireducerequest(7), ierr)
       END IF
     END IF
 
     CALL MPI_WAIT(ireducerequest(1), stat(1:MPI_STATUS_SIZE), ierr)
     CALL MPI_WAIT(ireducerequest(2), stat(MPI_STATUS_SIZE+1:2*MPI_STATUS_SIZE), ierr)
     etot=epot+ekin
 
     IF(mpisize .gt. 1) THEN
      Nplocs_all(mpirank)=parts%Nptot
      IF(mpirank .eq.0 ) THEN
         CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
         & 0, MPI_COMM_WORLD, ierr)
       ELSE 
         CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
         & 0, MPI_COMM_WORLD, ierr)
       END IF
-      IF(mpirank .eq.0 .and. INT(SUM(Nplocs_all)).ne. nplasma) THEN
+      IF(mpirank .eq. 0 .AND. partperiodic .and. INT(SUM(Nplocs_all)).ne. nplasma) THEN
         WRITE(*,*) "Error particle lost on some processus"
         !CALL MPI_abort(MPI_COMM_WORLD,-1,ierr)
       END IF
     END IF
 !
 !  Shift to Etot at cstep=1 (not valable yet at cstep=0!)
     IF(cstep.EQ.1) THEN
        etot0 = etot
     END IF
     IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend)) THEN
       CALL collectparts
     END IF
   end subroutine diagnostics
 
   SUBROUTINE collectparts
     USE basic, ONLY: mpirank, mpisize, ierr
     INTEGER, DIMENSION(:), ALLOCATABLE :: displs 
     INTEGER:: i
     
     IF(collected) RETURN ! exit subroutine if particles have already been collected this time step
 
     IF(mpirank .ne. 0) THEN
        ALLOCATE(displs(0:mpisize-1))
        displs=0
       ! Send Particles informations to root process
        CALL MPI_Gatherv(parts%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,&    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,&   
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%THET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%THET, Nplocs_all, displs,&   
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,&     
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, &    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, &    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, &    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%Rindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, &    
        & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%Zindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, &    
        & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(parts%partindex, Nplocs_all(mpirank), MPI_INTEGER, parts%partindex, Nplocs_all, displs, &    
        & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
     ELSE
        ALLOCATE(displs(0:mpisize-1))
        displs(0)=0
        Do i=1,mpisize-1
          displs(i)=displs(i-1)+Nplocs_all(i-1)
        END DO
   ! Receive particle information from all processes
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,&    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,&   
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%THET, Nplocs_all, displs,&   
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,&     
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, &    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, &    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, &    
        & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, &    
        & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, &    
        & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
        CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%partindex, Nplocs_all, displs, &    
        & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
     END IF
     collected=.TRUE.
   END SUBROUTINE collectparts
 !_______________________________________________________________________________
   SUBROUTINE adapt_vinit
     !
     !   Computes the velocity at time -dt/2 from velocities computed at time 0
     !
       USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nlclassical, phinorm, nz, distribtype, H0, vnorm
       DOUBLE PRECISION :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, &
       &           SQR, ZBZ2, ZBR2, Vperp, v2
       INTEGER :: J1, J2, J3, J4, i
       DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
       IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3 .OR. distribtype .EQ. 4) THEN
         ALLOCATE(VR(parts%Nptot),VZ(parts%Nptot),VTHET(parts%Nptot))
         CALL loduni(7,VZ)
         VZ=VZ*2*pi
         VTHET=parts%UTHET/parts%Gamma*vnorm
         DO i=1,parts%Nptot
           Vperp=sqrt(MAX(2*H0/me+2*elchar/me*parts%pot(i)*phinorm-VTHET(i)**2,0.0_db))
           VR(i)=Vperp*sin(VZ(i))
           VZ(i)=Vperp*cos(VZ(i))
           IF(nlclassical) THEN
             parts%Gamma(i)=1
           ELSE
             v2=VR(i)**2+VZ(i)**2+VTHET(i)**2
             parts%Gamma(i)=sqrt(1/(1-v2/vnorm**2))
           END IF
           parts%UR(i)=parts%Gamma(i)*VR(i)/vnorm
           parts%UZ(i)=parts%Gamma(i)*VZ(i)/vnorm
           parts%UTHET(i)=parts%Gamma(i)*VTHET(i)/vnorm
         END DO
         DEALLOCATE(VR,VZ,VTHET)
       END IF
 
 
       ! Normalised time increment
       tau=-omegac/2/omegap*dt/tnorm
       ! Store old Velocities
       CALL swappointer(parts%UZold, parts%UZ)  
       CALL swappointer(parts%URold, parts%UR)  
       CALL swappointer(parts%UTHETold, parts%UTHET)  
       CALL swappointer(parts%Gammaold, parts%Gamma)  
 
       IF (parts%Nptot .NE. 0) THEN
       !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
         DO i=1,parts%Nptot
            J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1
            J2=J1+1
            J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1
            J4=J3+1
         
         ! Interpolation for magnetic field
            BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) &
            &   +parts%WZ(i)*(1-parts%WR(i))*Bz(J3)    &
            &   +(1-parts%WZ(i))*parts%WR(i)*Bz(J2)    &
            &   +parts%WZ(i)*parts%WR(i)*Bz(J1)
            BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) &
            &   +parts%WZ(i)*(1-parts%WR(i))*Br(J3)    &
            &   +(1-parts%WZ(i))*parts%WR(i)*Br(J2)    &
            &   +parts%WZ(i)*parts%WR(i)*Br(J1)
 
            ! Half inverse Rotation along magnetic field
            ZBZ=tau*BRZ/parts%Gammaold(i)
            ZBR=tau*BRR/parts%Gammaold(i)
            SQR=1+ZBZ*ZBZ+ZBR*ZBR
            ZPZ=(parts%UZold(i)-ZBR*parts%UTHETold(i))/SQR                     !u-_{z} 
            ZPR=(parts%URold(i)+ZBZ*parts%UTHETold(i))/SQR                     !u-_{r}
            ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZold(i)-ZBZ*parts%URold(i))/SQR          !u-_{theta}
            parts%UZ(i)=ZPZ
            parts%UR(i)=ZPR
            parts%UTHET(i)=ZPTHET
 
            ! half of decceleration
            parts%UZ(i)=parts%UZ(i)+parts%Ez(i)*tau
            parts%UR(i)=parts%UR(i)+parts%Er(i)*tau
 
            IF(.not. nlclassical) THEN
             parts%Gamma(i)=sqrt(1+parts%UZ(i)**2+parts%UR(i)**2+parts%UTHET(i)**2)
            END IF
         
         END DO
         !$OMP END PARALLEL DO
       END IF
       
     END SUBROUTINE adapt_vinit
 !_______________________________________________________________________________
   SUBROUTINE distribpartsonprocs
   ! Computes the start and end indices for the Z boundaries on local processus
   ! Computes the particle indices from initial particle loading vector, that stay in current process
     USE basic, ONLY: nz, nplasma, ierr, Zbounds
 
     INTEGER:: k
     DOUBLE PRECISION:: idealnbpartsperproc
     INTEGER:: totparts ! Total number of particle from zgrid(0) to zgrid(k)
     INTEGER:: partindexstart ! Starting index of local particle in the parts variable used later to move the local particles at the begining of the vector
     INTEGER:: partindexlast   ! Ending index of local particle in the parts variable used later to move the local particles at the begining of the vector
     INTEGER, DIMENSION(0:nz):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1)
     INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction
     INTEGER:: Zperproc
     partspercol=0
     
     
     idealnbpartsperproc = FLOOR(REAL(nplasma)/REAL(mpisize))
     Zmin=MINVAL(parts%Zindex)
     Zmax=MAXVAL(parts%Zindex)
     Zperproc=(Zmax-Zmin)/mpisize
 
     partindexstart=1
     totparts=0
     partindexlast=parts%Nptot
     DO k=0,mpisize-2
       Nplocs_all(k)=idealnbpartsperproc
     END DO
       NPlocs_all(mpisize-1)=idealnbpartsperproc+MODULO(nplasma,mpisize)
     DO k=1, mpisize-1
       Zbounds(k)=parts%Z(INT(k*idealnbpartsperproc))
     END DO
 
     ! Store local end of particle vector
     partindexlast=SUM(Nplocs_all(0:mpirank))
     ! Store local start of particle vector
     IF(mpirank .ne. 0) partindexstart=SUM(Nplocs_all(0:mpirank-1))+1
 
     IF(mpirank .ne. 0) THEN
       DO k= partindexstart, partindexlast
         CALL move_part(k,k-partindexstart+1)
       END DO
     END IF
     parts%Nptot=Nplocs_all(mpirank)
     WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot
 
     IF(INT(SUM(Nplocs_all)) .ne. nplasma) THEN
       WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " nplasma:", nplasma
       CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
     END IF
 
   END SUBROUTINE distribpartsonprocs
 !_______________________________________________________________________________
 
   SUBROUTINE particlescommunication(lsendnbparts, rsendnbparts, sendholes)
     USE mpihelper, ONLY: particle_type
     USE basic, ONLY: rightproc, leftproc, ierr
     INTEGER, INTENT(in)    :: lsendnbparts, rsendnbparts
     INTEGER, INTENT(in) :: sendholes(parts%Nptot)
     INTEGER,  ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0
     INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendrequest, recvrequest
     INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendstatus(2*MPI_STATUS_SIZE), recvstatus(2*MPI_STATUS_SIZE)
     INTEGER :: lsentnbparts, rsentnbparts
     INTEGER :: lreceivednbparts, rreceivednbparts
 
     lsentnbparts=lsendnbparts
     rsentnbparts=rsendnbparts
 
     sendrequest=MPI_REQUEST_NULL
     recvrequest=MPI_REQUEST_NULL
 
     CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, leftproc,  0, MPI_COMM_WORLD, recvrequest(1), ierr)
     CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, recvrequest(2), ierr)
     CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, leftproc,  0, MPI_COMM_WORLD, sendrequest(1), ierr)
     CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, sendrequest(2), ierr)
 
     CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr)
     CALL MPI_Wait(recvrequest(2), recvstatus(2), ierr)
     recvrequest=MPI_REQUEST_NULL
 
     lreceivednbparts=lrecvnbparts
     rreceivednbparts=rrecvnbparts
 
     IF( ALLOCATED(rrecvpartbuff) .AND. (size(rrecvpartbuff) .lt. rrecvnbparts) ) DEALLOCATE(rrecvpartbuff)
     IF(.not. ALLOCATED(rrecvpartbuff) .and. rrecvnbparts .gt. 0) ALLOCATE(rrecvpartbuff(INT(rreceivednbparts*1.3)))
     IF( ALLOCATED(lrecvpartbuff) .AND. (size(lrecvpartbuff) .lt. lrecvnbparts) ) DEALLOCATE(lrecvpartbuff)
     IF((.not. ALLOCATED(lrecvpartbuff) ).and. lrecvnbparts .gt. 0) ALLOCATE(lrecvpartbuff(INT(lreceivednbparts*1.3)))
 
     ! Receive particles from left and right processes
     IF ( lrecvnbparts .gt. 0) THEN
       CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, leftproc,  1, MPI_COMM_WORLD, recvrequest(1), ierr)
     END IF
     IF( rrecvnbparts .gt. 0) THEN
       CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, recvrequest(2), ierr)
     END IF
 
 
     IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN 
       CALL AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes)
     END IF
 
     CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr)
     CALL MPI_Wait(sendrequest(2), sendstatus(MPI_STATUS_SIZE+1), ierr)
 
     IF( lsendnbparts .gt. 0) THEN
         CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, leftproc,  1, MPI_COMM_WORLD, sendrequest(1), ierr)
     END IF
     IF( rsendnbparts .gt. 0) THEN
      CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, sendrequest(2), ierr)
     END IF
     
     IF ( lreceivednbparts .gt. 0) THEN
       CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr)
       IF(ierr .ne. MPI_SUCCESS) THEN
         WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(1)
         CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
       END IF
     END IF
     IF ( rreceivednbparts .gt. 0) THEN 
       CALL MPI_Wait(recvrequest(2), recvstatus(MPI_STATUS_SIZE+1), ierr)
       IF(ierr .ne. MPI_SUCCESS) THEN
         WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(2)
         CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
       END IF
     END IF
 
 
     CALL Addincomingparts(rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, sendholes)
 
     IF( lsendnbparts .gt. 0) THEN 
       CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr)
     END IF
     IF( rsendnbparts .gt. 0) THEN
       CALL MPI_Wait(sendrequest(2), sendstatus(2), ierr)
     END IF
      
 !
 !
   END SUBROUTINE particlescommunication
   !_______________________________________________________________________________
   SUBROUTINE Addincomingparts(rrecvnbparts, lrecvnbparts, sendnbparts, sendholes)
 !
     USE mpihelper
     INTEGER, INTENT(in)    :: rrecvnbparts, lrecvnbparts, sendnbparts
     INTEGER, INTENT(in) :: sendholes(parts%Nptot)
     INTEGER k,partpos, partdiff
 
     ! computes if we received less particles than we sent
     partdiff=sendnbparts-rrecvnbparts-lrecvnbparts
     
 
     IF(rrecvnbparts .gt. 0) THEN
       Do k=1,rrecvnbparts
         IF(k .le. sendnbparts) THEN
           partpos=abs(sendholes(k))
         ELSE
           parts%Nptot=parts%Nptot+1
           partpos=parts%Nptot
         END IF
         CALL Insertincomingpart(rrecvpartbuff, k, partpos)
       END DO
     END IF
 
     IF(lrecvnbparts .gt. 0) THEN
       Do k=1,lrecvnbparts
         IF(k+rrecvnbparts .le. sendnbparts) THEN
           partpos=abs(sendholes(k+rrecvnbparts))
         ELSE
           parts%Nptot=parts%Nptot+1
           partpos=parts%Nptot
         END IF
         CALL Insertincomingpart(lrecvpartbuff, k, partpos)
       END DO
     END IF
 
     IF(partdiff .gt. 0) THEN
       DO k=rrecvnbparts+lrecvnbparts+1, sendnbparts
         CALL move_part(parts%Nptot,abs(sendholes(k)))
         parts%Nptot = parts%Nptot-1
       END DO
     END IF
 
 !
   END SUBROUTINE Addincomingparts
  !_______________________________________________________________________________
   SUBROUTINE Insertincomingpart(buffer, bufferindex, partsindex)
     USE mpihelper
     INTEGER, INTENT(in) :: bufferindex, partsindex
     TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer
       parts%partindex(partsindex) = buffer(bufferindex)%partindex
       parts%Rindex(partsindex) = buffer(bufferindex)%Rindex
       parts%Zindex(partsindex) = buffer(bufferindex)%Zindex
       parts%R(partsindex) = buffer(bufferindex)%R
       parts%Z(partsindex) = buffer(bufferindex)%Z
       parts%THET(partsindex) = buffer(bufferindex)%THET
       parts%UZ(partsindex) = buffer(bufferindex)%UZ
       parts%UR(partsindex) = buffer(bufferindex)%UR
       parts%UTHET(partsindex) = buffer(bufferindex)%UTHET
       parts%Gamma(partsindex) = buffer(bufferindex)%Gamma
 !
 !
   END SUBROUTINE Insertincomingpart
 !_______________________________________________________________________________
   SUBROUTINE Insertsentpart(buffer, bufferindex, partsindex)
     USE mpihelper
     INTEGER, INTENT(in) :: bufferindex, partsindex
     TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer
       buffer(bufferindex)%partindex = parts%partindex(partsindex)
       buffer(bufferindex)%Rindex    = parts%Rindex(partsindex)
       buffer(bufferindex)%Zindex    = parts%Zindex(partsindex)
       buffer(bufferindex)%R         = parts%R(partsindex)
       buffer(bufferindex)%Z         = parts%Z(partsindex)
       buffer(bufferindex)%THET      = parts%THET(partsindex)
       buffer(bufferindex)%UZ        = parts%UZ(partsindex) 
       buffer(bufferindex)%UR        = parts%UR(partsindex) 
       buffer(bufferindex)%UTHET     = parts%UTHET(partsindex)
       buffer(bufferindex)%Gamma     = parts%Gamma(partsindex)
 !
 !
   END SUBROUTINE Insertsentpart
 !_______________________________________________________________________________
   SUBROUTINE AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes)
 !
     USE mpihelper
     INTEGER, INTENT(in)    :: lsendnbparts, rsendnbparts
     INTEGER, INTENT(in) :: sendholes(parts%Nptot)
     INTEGER:: partpos, k
     INTEGER:: lsendpos, rsendpos
     lsendpos=0
     rsendpos=0
     IF( ALLOCATED(rsendpartbuff) .AND. size(rsendpartbuff) .lt. rsendnbparts ) DEALLOCATE(rsendpartbuff)
     IF( ALLOCATED(lsendpartbuff) .AND. size(lsendpartbuff) .lt. lsendnbparts ) DEALLOCATE(lsendpartbuff)
 
     IF(.not. ALLOCATED(rsendpartbuff) .and. rsendnbparts .gt. 0) ALLOCATE(rsendpartbuff(INT(rsendnbparts*1.3)))
     IF(.not. ALLOCATED(lsendpartbuff) .and. lsendnbparts .gt. 0) ALLOCATE(lsendpartbuff(INT(lsendnbparts*1.3)))
 
     Do k=lsendnbparts+rsendnbparts,1,-1
       partpos=abs(sendholes(k))
       IF(sendholes(k) .GT. 0) THEN
         rsendpos=rsendpos+1
         CALL Insertsentpart(rsendpartbuff, rsendpos, partpos)
       ELSE IF(sendholes(k) .LT. 0) THEN 
         lsendpos=lsendpos+1
         CALL Insertsentpart(lsendpartbuff, lsendpos, partpos)
       END IF
     END DO
 !
 !
   END SUBROUTINE AddPartSendBuffers
 !_______________________________________________________________________________
 
   SUBROUTINE exchange_parts(index1, index2)
   INTEGER, INTENT(IN) :: index1, index2
   DOUBLE PRECISION:: R, Z, THET, UR, UZ, UTHET, Gamma
   INTEGER :: Rindex, Zindex, partindex
   !! Exchange particle at index1 with particle at index2
 
   ! Store part at index1 in temporary value
   partindex= parts%partindex(index1)
   Gamma    = parts%Gamma(index1)
   R        = parts%R(index1)
   Z        = parts%Z(index1)
   THET     = parts%THET(index1)
   UR       = parts%UR(index1)
   UTHET    = parts%UTHET(index1)
   UZ       = parts%UZ(index1)
   Rindex   = parts%Rindex(index1)
   Zindex   = parts%Zindex(index1)
  
   ! Move part at index2 in part at index 1
   parts%partindex(index1)= parts%partindex(index2)
   parts%Gamma(index1)    = parts%Gamma(index2)   
   parts%R(index1)        = parts%R(index2)
   parts%Z(index1)        = parts%Z(index2)
   parts%THET(index1)     = parts%THET(index2)
   parts%UR(index1)       = parts%UR(index2)     
   parts%UTHET(index1)    = parts%UTHET(index2)     
   parts%UZ(index1)       = parts%UZ(index2)       
   parts%Rindex(index1)   = parts%Rindex(index2)       
   parts%Zindex(index1)   = parts%Zindex(index2) 
 
   ! Move temporary values from part(index1) to part(index2)
   parts%partindex(index2)= partindex
   parts%Gamma(index2)    = Gamma   
   parts%R(index2)        = R
   parts%Z(index2)        = Z
   parts%THET(index2)     = THET
   parts%UR(index2)       = UR
   parts%UTHET(index2)    = UTHET
   parts%UZ(index2)       = UZ 
   parts%Rindex(index2)   = Rindex
   parts%Zindex(index2)   = Zindex
 
   END SUBROUTINE exchange_parts
 !_______________________________________________________________________________
 
   SUBROUTINE move_part(sourceindex, destindex)
   !! Move particle with index sourceindex to particle with index destindex
   !! This will destroy particle at destindex
   INTEGER, INTENT(IN) :: destindex, sourceindex
 
   ! Move part at sourceindex in part at destindex
   parts%partindex(destindex)= parts%partindex(sourceindex)
   parts%Gamma(destindex)    = parts%Gamma(sourceindex)   
   parts%R(destindex)        = parts%R(sourceindex)
   parts%Z(destindex)        = parts%Z(sourceindex)
   parts%THET(destindex)     = parts%THET(sourceindex)
   parts%UR(destindex)       = parts%UR(sourceindex)     
   parts%UTHET(destindex)    = parts%UTHET(sourceindex)        
   parts%UZ(destindex)       = parts%UZ(sourceindex)                    
   parts%Rindex(destindex)   = parts%Rindex(sourceindex)       
   parts%Zindex(destindex)   = parts%Zindex(sourceindex) 
 
   END SUBROUTINE move_part
 !_______________________________________________________________________________
 !________________________________________________________________________________
 !
   SUBROUTINE loduni(nbase, y)
 !
 !   Load a uniform distribution using the Hammersley's sequence.
 !   (nbase = 0 => Random sampling)
 !
     INTEGER, INTENT(in)        :: nbase
     DOUBLE PRECISION, INTENT(out) :: y(:)
     INTEGER       :: n, i
 !
     n = SIZE(y)
 !
     SELECT CASE (nbase)
     CASE(0)
        CALL RANDOM_NUMBER(y)
     CASE(1)
        DO i=1,n
           y(i) = (i-0.5_db)/n
        END DO
     CASE(2:)
        DO i=1,n
           y(i) = rev(nbase,i)
        END DO
     CASE default
        WRITE(*,'(a,i5)') 'Invalid value of NBASE =', nbase
     END SELECT
 !
   END SUBROUTINE loduni
 !________________________________________________________________________________
   SUBROUTINE lodgaus(nbase, y)
 !
 !   Sample y from the Gauss distributrion
 !
     DOUBLE PRECISION, INTENT(out) :: y(:)
     INTEGER, INTENT(in)        :: nbase
 !
     INTEGER       :: n, i, iflag
     DOUBLE PRECISION :: r(SIZE(y))
     DOUBLE PRECISION :: t
 !
     DOUBLE PRECISION :: c0, c1, c2
     DOUBLE PRECISION :: d1, d2, d3
     DATA c0, c1, c2/2.515517_db, 0.802853_db, 0.010328_db/
     DATA d1, d2, d3/1.432788_db, 0.189269_db, 0.001308_db/
 !
     n = SIZE(y)
     CALL loduni(nbase, r)
     r = 1.0E-5_db + 0.99998_db*r
 !
     DO i=1,n
        iflag = 1
        IF (r(i) .GT. 0.5_db) THEN
           r(i) = 1.0_db - r(i)
           iflag = -1
        END IF
        t = SQRT(LOG(1.0_db/(r(i)*r(i))))
        y(i) = t - (c0+c1*t+c2*t**2) / (1.0_db+d1*t+d2*t**2+d3*t**3)
        y(i) = y(i) * iflag
     END DO
     y = y - SUM(y)/REAL(n,db)
   END SUBROUTINE lodgaus
   !________________________________________________________________________________
   SUBROUTINE lodinvr(nbase, y, ra, rb)
     !
     !   Sample y from the distribution f=1/(r)H(r-ra)H(rb-r) for where H is the heavyside function 
     !
         DOUBLE PRECISION, INTENT(out) :: y(:)
         INTEGER, INTENT(in)        :: nbase
         DOUBLE PRECISION, INTENT(in)  :: rb, ra
     !
         DOUBLE PRECISION :: r(SIZE(y))
     !
         CALL loduni(nbase, r)
         r=r*log(rb/ra)
         y=ra*exp(r)
   END SUBROUTINE lodinvr
   !________________________________________________________________________________
   SUBROUTINE lodunir(nbase, y, ra, rb)
         DOUBLE PRECISION, INTENT(out) :: y(:)
         INTEGER, INTENT(in)        :: nbase
         DOUBLE PRECISION, INTENT(in)  :: rb, ra
     !
         DOUBLE PRECISION :: r(SIZE(y))
     !
         CALL loduni(nbase, r)
         y=ra+(rb-ra)*r
   END SUBROUTINE lodunir
   !________________________________________________________________________________
   SUBROUTINE lodgausr(nbase, y, ra, rb)
         DOUBLE PRECISION, INTENT(out) :: y(:)
         INTEGER, INTENT(in)        :: nbase
         DOUBLE PRECISION, INTENT(in)  :: rb, ra
     !
         DOUBLE PRECISION :: r(SIZE(y))
     !
         CALL lodgaus(nbase, r)
         y=0.5*(rb+ra)+ 0.1*(rb-ra)*r
   END SUBROUTINE lodgausr
 !________________________________________________________________________________
    REAL(db) FUNCTION rev(nbase,i)
 !
 !   Return an element of the Hammersley's sequence
 !
     INTEGER, INTENT(IN) :: nbase   ! Base of the Hammersley's sequence (IN)
     INTEGER, INTENT(IN) :: i       ! Index of the sequence (IN) 
 !
 !   Local vars
     INTEGER :: j1, j2
     REAL(db) :: xs, xsi
 !-----------------------------------------------------------------------
     xs = 0._db
     xsi = 1.0_db
     j2 = i
     DO
        xsi = xsi/nbase
        j1 = j2/nbase
        xs = xs + (j2-nbase*j1)*xsi
        j2 = j1
        IF( j2.LE.0 )  EXIT
     END DO
     rev = xs
   END FUNCTION rev
 !________________________________________________________________________________
   !Modified Bessel functions of the first kind of the first order
   FUNCTION bessi1(x)
     DOUBLE PRECISION ::  bessi1,x
     DOUBLE PRECISION ::  ax
     DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
     SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
     DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
     DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2,  0.163801d-2,-0.1031555d-1, &
          &                           0.2282967d-1, -0.2895312d-1,  0.1787654d-1,-0.420059d-2/
     if (abs(x).lt.3.75) then
        y=(x/3.75)**2
        bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
     else
        ax=abs(x)
        y=3.75/ax
        bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
        if(x.lt.0.)bessi1=-bessi1
     endif
     return
   END FUNCTION bessi1
 !________________________________________________________________________________
 
   SUBROUTINE destroy_parts(p)
     TYPE(particles)       :: p
     p%nptot=0
     IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z) 
     IF(ALLOCATED(p%R)) DEALLOCATE(p%R)
     IF(ALLOCATED(p%THET)) DEALLOCATE(p%THET)
     IF(ALLOCATED(p%WZ)) DEALLOCATE(p%WZ) 
     IF(ALLOCATED(p%WR)) DEALLOCATE(p%WR)
     IF(ASSOCIATED(p%UR)) DEALLOCATE(p%UR) 
     IF(Associated(p%URold)) DEALLOCATE(p%URold)
     IF(Associated(p%UZ)) DEALLOCATE(p%UZ) 
     IF(Associated(p%UZold)) DEALLOCATE(p%UZold)
     IF(Associated(p%UTHET)) DEALLOCATE(p%UTHET) 
     IF(Associated(p%UTHETold)) DEALLOCATE(p%UTHETold)
     IF(Associated(p%Gamma)) DEALLOCATE(p%Gamma) 
     IF(Associated(p%Gammaold)) DEALLOCATE(p%Gammaold) 
     IF(ALLOCATED(p%Rindex)) DEALLOCATE(p%Rindex)
     IF(ALLOCATED(p%Zindex)) DEALLOCATE(p%Zindex)
     IF(ALLOCATED(p%partindex)) DEALLOCATE(p%partindex)
   END SUBROUTINE
 !________________________________________________________________________________
   SUBROUTINE clean_beam
 !
     CALL destroy_parts(parts)
 !
 !
   END SUBROUTINE clean_beam
 !________________________________________________________________________________
   RECURSIVE SUBROUTINE quicksortparts(leftlimit, rightlimit)
   ! Sorts the particle according to their Z position using quicksort algorithm
   INTEGER,INTENT(IN):: leftlimit, rightlimit
   DOUBLE PRECISION:: pivot
   INTEGER::i, cnt, mid
     IF(leftlimit .ge. rightlimit) RETURN
     mid=(leftlimit+rightlimit)/2
     IF(parts%Z(mid).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,mid)
     IF(parts%Z(rightlimit).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,rightlimit)
     IF(parts%Z(mid).lt.parts%Z(rightlimit)) CALL exchange_parts(rightlimit,mid)
 
     ! Store the pivot point for comparison
     pivot=parts%Z(rightlimit)
 
     cnt=leftlimit
     
     ! Move all parts with Z smaller than pivot to the left of pivot
     DO i=leftlimit, rightlimit
       IF(parts%Z(i) .le. pivot) THEN
         CALL exchange_parts(i,cnt)
         cnt=cnt+1
       END IF
     END DO
 
     CALL quicksortparts(leftlimit,cnt-2)
     CALL quicksortparts(cnt,rightlimit)
 
   END SUBROUTINE quicksortparts
 !________________________________________________________________________________
   SUBROUTINE swappointer( pointer1, pointer2)
   DOUBLE PRECISION, DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2 
   DOUBLE PRECISION, DIMENSION(:), POINTER:: temppointer
   temppointer=>pointer1
   pointer1=>pointer2
   pointer2=>temppointer
   END SUBROUTINE swappointer
 !_______________________________________________________________________________
   SUBROUTINE loaduniformRZ(VR,VZ,VTHET)
     USE basic, ONLY: plasmadim, rnorm, temp, vnorm
     USE constants, ONLY: me, kb, elchar
     DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET
     DOUBLE PRECISION:: vth
   ! Initial distribution in z with normalisation
       CALL loduni(1,parts%Z(:))
       parts%Z(:)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(:))/rnorm
     ! Initial distribution in r with normalisation
       CALL loduni(2,parts%R(:))
       parts%R=(plasmadim(3)+parts%R(:)*(plasmadim(4)-plasmadim(3)))/rnorm
   
     ! Initial velocities distribution
       vth=sqrt(3*kb*temp/me)/vnorm        !thermal velocity 
       CALL lodgaus(3,VZ(:))
       CALL lodgaus(5,VR(:))
       CALL lodgaus(7,VTHET(:))
       VZ=VZ*vth
       VR=VR*vth
       VTHET=VTHET*vth
   END SUBROUTINE loaduniformRZ
 !_______________________________________________________________________________
   SUBROUTINE loadDavidson(VR,VZ,VTHET,lodr)
     USE constants, ONLY: me, kb, elchar
     USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, B0, width, qsim, msim, vnorm, &
     &                omegac, zgrid, nz, rnorm, V, n0, nblock, temp
     interface
       subroutine lodr(nbase,y,ra,rb)
         DOUBLE PRECISION, INTENT(out) :: y(:)
         INTEGER, INTENT(in)        :: nbase
         DOUBLE PRECISION, INTENT(in)  :: rb, ra
       end subroutine
     end interface
     DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET
     DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE::ra, rb
     DOUBLE PRECISION :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL, vth
     INTEGER :: i, j, n, blockstart, blockend, addedpart, remainparts
     INTEGER, DIMENSION(:), ALLOCATABLE ::  blocksize
       Allocate(ra(nblock),rb(nblock))
       !r0=(plasmadim(4)+plasmadim(3))/2
       r0=sqrt(2*H0/(me*omegac**2))
       halfLz=(zgrid(nz)+zgrid(0))/2
       MirrorRatio=(Rcurv-1)/(Rcurv+1)
 
       DO n=1,nblock
       ! Compute limits in radius and load radii for each part
         Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1)
         deltar2=1-MirrorRatio*cos(2*pi*Le/width)
         rb(n)=r0*(1+sqrt(1-abs(omegac)*P0*deltar2/H0))/deltar2
         ra(n)=r0*(1-sqrt(1-abs(omegac)*P0*deltar2/H0))/deltar2
       END DO
       
       VOL=SUM(2*pi*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock)
       qsim=VOL*n0*elchar/nplasma 
       msim=abs(qsim)/elchar*me
       V=VOL/rnorm**3
 
       !H0=me*omegac**2*r0**2*0.5
       gamma=1/sqrt(1-omegac**2*r0**2/vlight**2)
       !P0=H0/omegac
       
       blockstart=1
       blockend=0
       ALLOCATE(blocksize(nblock))
       DO n=1,nblock
         blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock
       END DO
       remainparts=parts%Nptot-SUM(blocksize)
       addedpart=1
       n=nblock/2
       j=1
       DO WHILE(remainparts .GT. 0) 
          blocksize(n)=blocksize(n)+addedpart
          remainparts=remainparts-addedpart
          n=n+j
          j=-1*(j+SIGN(1,j))
       END DO
       DO n=1,nblock
         blockstart=blockend+1
         blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot)
     ! Initial distribution in z with normalisation between magnetic mirrors
         CALL loduni(1,parts%Z(blockstart:blockend))
         parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm
         CALL lodr(2,parts%R(blockstart:blockend),ra(n), rb(n))
         parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm
       END DO
       
       
     IF(distribtype .eq. 5) THEN
     ! Initial velocities distribution
       vth=sqrt(3*kb*temp/me)/vnorm        !thermal velocity 
       CALL lodgaus(3,VZ(:))
       CALL lodgaus(5,VR(:))
       CALL lodgaus(7,VTHET(:))
       VZ=VZ*vth/4
       VR=VR*8*vth
       VTHET=VTHET*8*vth
     ELSE
     ! Load velocities theta velocity
     ! Loading of r and z velocity is done in adapt_vinit to have 
     ! access to parts%pot 
       DO i=1,parts%Nptot
     ! Interpolation for Magnetic potential
        rg=parts%R(i)*rnorm
        zg=(parts%Z(i)-halfLz)*rnorm
       
        Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width))
        Pcomp=abs(P0/rg/me)
        Acomp=-SIGN(elchar/me*Athetpos,qsim)
        VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp)
       END DO
       VTHET=VTHET/vnorm
       VZ=0._db
       VR=0._db
     END IF
 
   END SUBROUTINE loadDavidson
 END MODULE beam
diff --git a/src/mpihelper_mod.f90 b/src/mpihelper_mod.f90
index 75f1e9d..c9b752f 100644
--- a/src/mpihelper_mod.f90
+++ b/src/mpihelper_mod.f90
@@ -1,140 +1,145 @@
 MODULE mpihelper
 USE constants
 USE mpi
 IMPLICIT NONE
 
 TYPE BASICDATA
     LOGICAL :: nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical
     INTEGER :: nplasma, nz, it0d, it2d, itparts,  &
             &  itgraph, distribtype, nblock, nrun
     DOUBLE PRECISION :: job_time, extra_time, tmax, dt,   &
                   &  potinn, potout, B0, n0, temp, & 
                   &  Rcurv, width, H0, P0
     INTEGER, DIMENSION(2) :: femorder, ngauss
     INTEGER, DIMENSION(3) :: nnr
     DOUBLE PRECISION, DIMENSION(2):: lz
     DOUBLE PRECISION, DIMENSION(4):: radii, plasmadim
     CHARACTER(len=64) :: resfile= "results.h5"
+    LOGICAL :: partperiodic
   END TYPE BASICDATA
 
   TYPE particle
   INTEGER :: Rindex, Zindex, partindex
   DOUBLE PRECISION :: R, Z, THET,& 
                &  UZ, UR, UTHET, &
                &  Gamma
   END TYPE particle
 
   INTEGER, SAVE :: basicdata_type
   INTEGER, SAVE :: particle_type
   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: rhscol_type
   INTEGER, SAVE:: rhsoverlap_type
   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE:: rhsrequests, rhsoverlaprequests
   INTEGER:: rhscommtag=10, rhsoverlapcommtag = 15
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: rhsoverlap
 
   CONTAINS 
   !=========================================================================
   SUBROUTINE init_mpitypes
 
     CALL init_particlempi
     CALL init_basicdatampi
 
   END SUBROUTINE init_mpitypes
   !=========================================================================
   SUBROUTINE init_particlempi()
       INTEGER :: nblock = 10
       INTEGER:: blocklength(10)
       INTEGER(kind=MPI_ADDRESS_KIND):: displs(10), displ0
       INTEGER, DIMENSION(10) :: types
       TYPE(particle) :: part
       INTEGER:: ierr
        
       CALL mpi_get_address(part%Rindex, displs(1), ierr)
       CALL mpi_get_address(part%Zindex, displs(2), ierr)
       CALL mpi_get_address(part%partindex, displs(3), ierr)
       types(1:3)=MPI_INTEGER
       CALL mpi_get_address(part%R, displs(4), ierr)
       CALL mpi_get_address(part%Z, displs(5), ierr)
       CALL mpi_get_address(part%THET, displs(6), ierr)
       CALL mpi_get_address(part%UZ, displs(7), ierr)
       CALL mpi_get_address(part%UR, displs(8), ierr)
       CALL mpi_get_address(part%UTHET, displs(9), ierr)
       CALL mpi_get_address(part%GAMMA, displs(10), ierr)
       types(4:10)=MPI_DOUBLE_PRECISION
       blocklength(1:10) = 1
       CALL mpi_get_address(part, displ0, ierr)
       displs=displs-displ0
   
 
       CALL MPI_Type_create_struct(nblock, blocklength, displs, types, particle_type, ierr)
       CALL MPI_Type_commit(particle_type,ierr)
 
   END SUBROUTINE init_particlempi
 
 !=========================================================================
   SUBROUTINE init_basicdatampi()
-      INTEGER :: nblock = 36
-      INTEGER:: blocklength(36)
-      INTEGER(kind=MPI_ADDRESS_KIND):: displs(36), displ0
-      INTEGER, DIMENSION(36) :: types
+      INTEGER :: nblock = 37
+      INTEGER:: blocklength(37)
+      INTEGER(kind=MPI_ADDRESS_KIND):: displs(37), displ0
+      INTEGER, DIMENSION(37) :: types
       TYPE(BASICDATA) :: bdata
       INTEGER:: ierr
        
       CALL mpi_get_address(bdata%nlres, displs(1), ierr)
       CALL mpi_get_address(bdata%nlsave, displs(2), ierr)
       CALL mpi_get_address(bdata%newres, displs(3), ierr)
       CALL mpi_get_address(bdata%nlxg, displs(4), ierr)
       CALL mpi_get_address(bdata%nlppform, displs(5), ierr)
       CALL mpi_get_address(bdata%nlPhis, displs(6), ierr)
       CALL mpi_get_address(bdata%nlclassical, displs(7), ierr)
       types(1:7)=MPI_LOGICAL
       CALL mpi_get_address(bdata%nplasma, displs(8), ierr)
       CALL mpi_get_address(bdata%nz, displs(9), ierr)
       CALL mpi_get_address(bdata%it0d, displs(10), ierr)
       CALL mpi_get_address(bdata%it2d, displs(11), ierr)
       CALL mpi_get_address(bdata%itparts, displs(12), ierr)
       CALL mpi_get_address(bdata%itgraph, displs(13), ierr)
       CALL mpi_get_address(bdata%distribtype, displs(14), ierr)
       CALL mpi_get_address(bdata%nblock, displs(15), ierr)
       CALL mpi_get_address(bdata%nrun, displs(16), ierr)
       types(8:16)=MPI_INTEGER
       CALL mpi_get_address(bdata%job_time, displs(17), ierr)
       CALL mpi_get_address(bdata%extra_time, displs(18), ierr)
       CALL mpi_get_address(bdata%tmax, displs(19), ierr)
       CALL mpi_get_address(bdata%dt, displs(20), ierr)
       CALL mpi_get_address(bdata%potinn, displs(21), ierr)
       CALL mpi_get_address(bdata%potout, displs(22), ierr)
       CALL mpi_get_address(bdata%B0, displs(23), ierr)
       CALL mpi_get_address(bdata%n0, displs(24), ierr)
       CALL mpi_get_address(bdata%temp, displs(25), ierr)
       CALL mpi_get_address(bdata%Rcurv, displs(26), ierr)
       CALL mpi_get_address(bdata%width, displs(27), ierr)
       CALL mpi_get_address(bdata%H0, displs(28), ierr)
       CALL mpi_get_address(bdata%P0, displs(29), ierr)
       types(17:29)=MPI_DOUBLE_PRECISION
       blocklength(1:29) = 1
       CALL mpi_get_address(bdata%femorder, displs(30), ierr)
       CALL mpi_get_address(bdata%ngauss, displs(31), ierr)
       CALL mpi_get_address(bdata%nnr, displs(32), ierr)
       types(30:32)=MPI_INTEGER
       blocklength(30:31) = 2
       blocklength(32) = 3
       CALL mpi_get_address(bdata%lz, displs(33), ierr)
       blocklength(33) = 2
       CALL mpi_get_address(bdata%radii, displs(34), ierr)
       CALL mpi_get_address(bdata%plasmadim, displs(35), ierr)
       types(33:35)=MPI_DOUBLE_PRECISION
       blocklength(34:35) = 4
       CALL mpi_get_address(bdata%resfile, displs(36), ierr)
       types(36)=MPI_CHARACTER
       blocklength(36) = 64
+      CALL mpi_get_address(bdata%partperiodic, displs(37), ierr)
+      types(37)=MPI_LOGICAL
+      blocklength(37) = 1
+
       CALL mpi_get_address(bdata, displ0, ierr)
       displs=displs-displ0
 
       CALL MPI_Type_create_struct(nblock, blocklength, displs, types, basicdata_type, ierr)
       CALL MPI_Type_commit(basicdata_type,ierr)
 
   END SUBROUTINE init_basicdatampi
 !=========================================================================
 
 END MODULE mpihelper