diff --git a/matlab/ForceBalance.m b/matlab/ForceBalance.m
new file mode 100644
index 0000000..ae2b3cd
--- /dev/null
+++ b/matlab/ForceBalance.m
@@ -0,0 +1,78 @@
+function ForceBalance(M,it)
+%ForceBalance Show the radial force balance
+%   Plot the three radial forces for the given time-step it or averaged
+%   over the range of time steps defined in it
+
+if it==':'
+    it=1:size(M.t2d);
+end
+Eforce=-M.qe*mean(M.Er(:,:,it),3);
+Bforce=-M.qe*mean(M.fluidUTHET(:,:,it),3).*M.Bz';
+N=mean(M.N(:,:,it),3);
+maxdens=max(max(N));
+z=(M.zgrid(1:end-1)+M.zgrid(2:end))/2;
+r=(M.rgrid(1:end-1)+M.rgrid(2:end))/2;
+[R,~]=meshgrid(M.rgrid,M.zgrid);
+Rinv=1./R;
+Rinv(:,1)=0;
+inertforce=M.me*mean(M.fluidUTHET(:,:,it),3).^2.*Rinv';
+
+
+rgridmax=50;
+f=figure();
+
+ax1=subplot(3,1,1);
+surface(ax1,M.zgrid,M.rgrid,Eforce)
+hold on;
+%contour(ax1,z,r,N(1:end-1,1:end-1),[0.8 0.8]*maxdens,'r-','Linewidth',1.3)
+xlabel(ax1,'z [m]')
+ylabel(ax1,'r [m]')
+xlim(ax1,[M.zgrid(1) M.zgrid(end)])
+ylim(ax1,[M.rgrid(1) M.rgrid(rgridmax)])
+colormap(ax1,'jet')
+c = colorbar(ax1);
+c.Label.String= 'F_E [N]';
+title(ax1,'Electric force')
+grid(ax1, 'on')
+view(ax1,2)
+
+
+
+ax2=subplot(3,1,2);
+surface(ax2,M.zgrid,M.rgrid,Bforce)
+xlabel(ax2,'z [m]')
+ylabel(ax2,'r [m]')
+xlim(ax2,[M.zgrid(1) M.zgrid(end)])
+ylim(ax2,[M.rgrid(1) M.rgrid(rgridmax)])
+colormap(ax2,'jet')
+c = colorbar(ax2);
+c.Label.String= 'F_B [N]';
+title(ax2,'Magnetic force')
+grid(ax2, 'on')
+view(ax2,2)
+
+ax3=subplot(3,1,3);
+surface(ax3,M.zgrid,M.rgrid,inertforce)
+xlabel(ax3,'z [m]')
+ylabel(ax3,'r [m]')
+xlim(ax3,[M.zgrid(1) M.zgrid(end)])
+ylim(ax3,[M.rgrid(1) M.rgrid(rgridmax)])
+colormap(ax3,'jet')
+c = colorbar(ax3);
+c.Label.String= 'F_i [N]';
+title(ax3,'Centrifugal force')
+grid(ax3, 'on')
+view(ax3,2)
+
+linkaxes([ax1 ax2 ax3],'xy')
+sgtitle(sprintf('Radial forces t=[%d-%d]s',M.t2d(min(it)),M.t2d(max(it))))
+
+f.PaperOrientation='landscape';
+[~, name, ~] = fileparts(M.file);
+f.PaperUnits='centimeters';
+papsize=[16 14];
+f.PaperSize=papsize;
+%print(f,sprintf('%sforce_balance',name),'-dpdf','-fillpage')
+
+end
+
diff --git a/src/beam_mod.f90 b/src/beam_mod.f90
index f553f2c..4fee857 100644
--- a/src/beam_mod.f90
+++ b/src/beam_mod.f90
@@ -1,1343 +1,1343 @@
 !------------------------------------------------------------------------------
 ! 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 = VZ/50
       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, 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, it0d, it2d, nlend, nr, itparts, nplasma, fluidur, fluiduz, fluiduthet, Presstens, partperiodic, newres
     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)
+    CALL MPI_WAITALL(7, ireducerequest, stat, 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)
+        ekin=ekin+msim*vlight**2*(parts%Gamma(i)-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)
       ! Normalize fluid velocities by density 
       fluiduz=fluiduz/partdensity
       fluidur=fluidur/partdensity
       fluiduthet=fluiduthet/partdensity
     END WHERE
 
     ! Computes the pressure tensor
     Presstens=0
     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
 
     ! 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(modulo(step,it0d).eq. 0 .and. 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. 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 .or. (newres .AND. step .EQ. 1)) THEN
        etot0 = etot
     END IF
     !Calls the communications to send the particles data to the root process for diagnostic file
     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 during 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(parts%Nptot)/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
 
     IF(INT(SUM(Nplocs_all)) .ne. parts%Nptot) THEN
       WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " Nptot:", parts%Nptot
       CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
     END IF
 
     parts%Nptot=Nplocs_all(mpirank)
     WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot
 
   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 lodlinr(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+sqrt(r)*(rb-ra)
   END SUBROUTINE lodlinr
   !________________________________________________________________________________
   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(4*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/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2+sqrt(1-P0*abs(omegac)/H0*deltar2))
         ra(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2-sqrt(1-P0*abs(omegac)/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)
+       Pcomp=P0/rg/me
        Acomp=-SIGN(elchar/me*Athetpos,qsim)
-       VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp)
+       !VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp)
+       VTHET(i)=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/fields_mod.f90 b/src/fields_mod.f90
index 5f8f26e..ac0219a 100644
--- a/src/fields_mod.f90
+++ b/src/fields_mod.f90
@@ -1,324 +1,326 @@
 MODULE fields
 !
 !   Module containing all variables and routines used to obtain |E(r,z) and |B(r,z)
 !
 
 USE constants
 USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, splrz, nlperiod, phinorm, nlPhis, nrank, rhs, mpirank, mpisize
 USE beam, ONLY: parts
 USE bsplines
 USE mumps_bsplines
 USE mpi
 IMPLICIT NONE
 
 DOUBLE PRECISION, allocatable, SAVE :: matcoef(:,:), sol(:), vec1(:), vec2(:)
 TYPE(mumps_mat), SAVE   :: femat                  ! FEM matrix
 
 
 CONTAINS
 !________________________________________________________________________________
 SUBROUTINE fields_init
     USE basic, ONLY: pot, nlperiod, nrank, rhs
     USE bsplines
     USE mumps_bsplines
     INTEGER :: nrz(2)
 
     ! Set up 2d spline splrz
     CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod)
 
     ! Calculate dimension of splines
       nrz(1)=nz
       nrz(2)=nr 
     CALL get_dim(splrz, nrank, nrz, femorder)
 
     ALLOCATE(matcoef(nrank(1),nrank(2)))
     ALLOCATE(pot((nr+1)*(nz+1)))
     ALLOCATE(rhs(nrank(1)*nrank(2)))
     ALLOCATE(sol(nrank(1)*nrank(2))) 
 
     ! Calculate magnetic field components in grid points (Davidson analitical formula employed)
     ALLOCATE(Br((nr+1)*(nz+1)),Bz((nr+1)*(nz+1)))
     ALLOCATE(Athet((nr+1)*(nz+1)))
     ALLOCATE(Er((nr+1)*(nz+1)),Ez((nr+1)*(nz+1)))
     
     CALL magnet
   
   ! Calculate and factorise FEM matrix (depended only on mesh)
     CALL fematrix
     CALL factor(femat)
 
     
 END SUBROUTINE fields_init
 !________________________________________________________________________________
 
 SUBROUTINE rhscon
     USE bsplines
     USE mumps_bsplines
     USE mpi
     USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma
     USE beam, ONLY: parts
     USE mpihelper
     INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch
     INTEGER,DIMENSION(:),ALLOCATABLE::zleft, rleft
+    DOUBLE PRECISION :: chargecoeff
     DOUBLE PRECISION, ALLOCATABLE :: fun(:,:,:), fun2(:,:,:)
     nbunch=100
     ALLOCATE(zleft(nbunch),rleft(nbunch))
     ALLOCATE(fun(1:femorder(1)+1,0:0,nbunch),fun2(1:femorder(2)+1,0:0,nbunch))!Arrays keeping values of b-splines at gauss node
     
     rhs=0 !Initialise rhs
+    chargecoeff=0.5/pi*omegap/omegac*V/nplasma
 ! Assemble rhs
     IF (parts%Nptot .ne. 0) THEN
 !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(zleft, rleft, jw, it, i, iend, irow, jcol, mu, k) PRIVATE(fun, fun2) REDUCTION(+:rhs)
       DO i=1,parts%Nptot,nbunch
          iend=min(i+nbunch-1,parts%Nptot)
          CALL locintv(splrz%sp1, parts%Z(i:iend), zleft)
          CALL locintv(splrz%sp2, parts%R(i:iend), rleft)
          CALL basfun(parts%Z(i:iend), splrz%sp1, fun, zleft+1)
          CALL basfun(parts%R(i:iend), splrz%sp2, fun2, rleft+1)
          
          DO k=1,nbunch
           DO jw=1,(femorder(2)+1)
            DO it=1,(femorder(1)+1)
                irow=zleft(k)+it
                jcol=rleft(k)+jw   
                mu=irow+(jcol-1)*(nrank(1))
-               rhs(mu)=rhs(mu)+fun(it,0,k)*fun2(jw,0,k)/2/pi*omegap/omegac*V/nplasma
+               rhs(mu)=rhs(mu)+fun(it,0,k)*fun2(jw,0,k)*chargecoeff
             END DO
           END DO
          END DO
       END DO
 !$OMP END PARALLEL DO
     END IF
 
     IF(mpisize .gt. 1) THEN
       IF (mpirank.eq.0) THEN
         CALL MPI_REDUCE(MPI_IN_PLACE, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, &
                       & MPI_SUM,  0, MPI_COMM_WORLD,  ierr)
       ELSE
         CALL MPI_REDUCE(rhs, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, &
                       & MPI_SUM,  0, MPI_COMM_WORLD, ierr)
       END IF
     END IF
 ! Dirichlet BC on RHS
     IF(rgrid(0).ne.0.0_db) rhs(1:nrank(1))=potinn
     rhs((nrank(2)-1)*nrank(1)+1:nrank(2)*nrank(1))=potout
     
 
     DEALLOCATE(fun,fun2, zleft, rleft)
   END SUBROUTINE rhscon
 !========================================================================!
 
   SUBROUTINE poisson
     USE bsplines
     USE mumps_bsplines
     USE futils
     INTEGER:: partend, i, ierr, jder(2)
     jder(1)=0
     jder(2)=0
     partend=parts%Nptot
     IF(mpirank.eq.0) CALL bsolve(femat,rhs,sol)
     
     CALL MPI_Bcast(sol(1), nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
 
     matcoef= reshape(sol, (/nrank(1),nrank(2)/))
     CALL gridval(splrz, vec1, vec2, pot, jder, matcoef)
     CALL getgrad(splrz, parts%Z(1:partend), parts%R(1:partend), parts%pot(1:partend), parts%EZ(1:partend), parts%ER(1:partend))
     IF (nlphis) THEN
       Do i=1,partend
         parts%EZ(i)=-parts%Ez(i)
         parts%ER(i)=-parts%ER(i)
       END DO
     ELSE
       Do i=1,partend
         parts%EZ(i)=0
         parts%ER(i)=0
       END DO
     END IF
 
     IF( mpirank .eq. 0) THEN
       CALL gridval(splrz, vec1, vec2, Ez, (/1,0/))
       CALL gridval(splrz, vec1, vec2, Er, (/0,1/))
       Ez=-Ez
       Er=-Er
     END IF
 
   END SUBROUTINE poisson
 !========================================================================!
 
 SUBROUTINE fematrix
 
     USE bsplines
     USE mumps_bsplines
 
     DOUBLE PRECISION, ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:)
     DOUBLE PRECISION, ALLOCATABLE :: f(:,:), aux(:)
     DOUBLE PRECISION, ALLOCATABLE :: coefs(:)
     DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), fun2(:,:)
     DOUBLE PRECISION  :: contrib
     INTEGER, ALLOCATABLE :: idert(:), iderw(:) 
     INTEGER :: i, j, k, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms
     kterms=2 !Number of terms in weak form
 
     ALLOCATE(fun(1:femorder(1)+1,0:1),fun2(1:femorder(2)+1,0:1))!Arrays keeping values of b-splines at gauss node
     ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2)))   !Gaussian nodes and weights arrays 
     ALLOCATE(f((femorder(1)+1)*(femorder(2)+1),2),aux(femorder(1)+1))                                                                  !Auxiliary arrays ordering bsplines
     ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms))                                                                              !Pointers on the order of derivatives 
 
 ! Constuction of auxiliary array ordering bsplines in given interval
     DO i=1,(femorder(1)+1)
        aux(i)=i
     END DO
     DO i=1,(femorder(2)+1)
        f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),1)=aux
        f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),2)=i
     END DO
 
 ! Initialisation of FEM matrix
     CALL init(nrank(1)*nrank(2),2,femat)
    
 ! Assemble FEM matrix
     DO j=1,nr ! Loop on r position
        DO i=1,nz	! Loop on z position
          ! Computation of gauss weight and position in r and z direction for gaussian integration
           CALL get_gauss(splrz%sp1, ngauss(1), i, zg, wzg)
           CALL get_gauss(splrz%sp2, ngauss(2), j, rg, wrg)
           ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration
           DO k=1,ngauss(2)
              xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg
              xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k)
              wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg	
           END DO
           DO igauss=1,ngauss(1)*ngauss(2) ! Loop on gaussian weights and positions
              CALL basfun(xgauss(igauss,1), splrz%sp1, fun, i)
              CALL basfun(xgauss(igauss,2), splrz%sp2, fun2, j)
              CALL coefeq(xgauss(igauss,:), idert, iderw, coefs)
              DO iterm=1,kterms ! Loop on the two integration dimensions
                 DO jt=1,(1+femorder(1))*(femorder(2)+1) 
                    DO iw=1,(1+femorder(1))*(femorder(2)+1)
                       irow=i+f(jt,1)-1; jcol=j+f(jt,2)-1
                       irow2=i+f(iw,1)-1; jcol2=j+f(iw,2)-1
                       mu=irow+(jcol-1)*nrank(1)
                       mu2=irow2+(jcol2-1)*nrank(1)
                       contrib=fun(f(jt,1),idert(iterm))* fun(f(iw,1),idert(iterm))* &
                            &  fun2(f(jt,2),iderw(iterm))*fun2(f(iw,2),iderw(iterm))* &
                            &  wgauss(igauss)*xgauss(igauss,2)
                       CALL updtmat(femat, mu, mu2, contrib)	
                    END DO
                 END DO
              END DO
           END DO
        END DO
     END DO
 
 ! Impose Dirichlet boundary conditions on FEM matrix 
     CALL fe_dirichlet 
 
     DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg)   
     DEALLOCATE(f, aux)                                                                                                                                       
     DEALLOCATE(idert, iderw, coefs, fun,fun2)
   END SUBROUTINE fematrix
 
   SUBROUTINE fe_dirichlet
     DOUBLE PRECISION, ALLOCATABLE :: arr(:)
     INTEGER :: i
     ALLOCATE(arr(nrank(1)*nrank(2)))
     DO i=1,nrank(1)*nrank(2)
       CALL getrow(femat,i,arr)
     END DO
     DO i=1,nrank(1)
        IF(rgrid(0).ne.0.0_db) THEN
          arr=0; arr(i)=1; 
          CALL putrow(femat,i,arr)
        END IF
        arr=0; arr(nrank(1)*nrank(2)+1-i)=1 ; 
        CALL putrow(femat,nrank(1)*nrank(2)+1-i,arr)  
     END DO
     DEALLOCATE(arr)
   END SUBROUTINE fe_dirichlet
 !________________________________________________________________________________
   SUBROUTINE coefeq(x, idt, idw, c)
     DOUBLE PRECISION, INTENT(in) :: x(2)
     INTEGER, INTENT(out) :: idt(:), idw(:)
     DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt))
 
     c(1) = x(2)    
     idt(1) = 0     
     idw(1) = 1     
     c(2) = x(2)     
     idt(2) = 1     
     idw(2) = 0 
   END SUBROUTINE coefeq
 !________________________________________________________________________________
 SUBROUTINE magnet
     USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm
     USE constants, ONLY: Pi
     DOUBLE PRECISION :: rg, zg,halfLz, MirrorRatio
     INTEGER :: i, rindex
     halfLz=(zgrid(nz)+zgrid(0))/2
     MirrorRatio=(Rcurv-1)/(Rcurv+1)
     DO i=1,(nr+1)*(nz+1)
        rindex=(i-1)/(nz+1)
        rg=rgrid(rindex)
        zg=zgrid(i-rindex*(nz+1)-1)-halfLz
        Br(i)=-B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm)/bnorm
        Bz(i)=B0*(1-MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))/bnorm
        Athet(i)=0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm))
     END DO
   END SUBROUTINE magnet
 !________________________________________________________________________________
 !Modified Bessel functions of the first kind of the zero order
 FUNCTION bessi0(x)
     DOUBLE PRECISION :: bessi0,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/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/
     DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2 ,0.916281d-2, &
          &                          -0.2057706d-1, 0.2635537d-1,-0.1647633d-1, 0.392377d-2 /
     if (abs(x).lt.3.75) then
        y=(x/3.75)**2
        bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
     else
        ax=abs(x)
        y=3.75/ax
        bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
     endif
     return
   END FUNCTION bessi0
 !________________________________________________________________________________
 !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 clean_fields
     Use bsplines
     DEALLOCATE(matcoef)
     DEALLOCATE(pot)
     DEALLOCATE(rhs)
     DEALLOCATE(sol)
     DEALLOCATE(Br,Bz)
     DEALLOCATE(Er,Ez)
     DEALLOCATE(vec1,vec2)
     Call DESTROY_SP(splrz)
 
   END SUBROUTINE clean_fields
 
 END MODULE fields
diff --git a/src/resume.f90 b/src/resume.f90
index 0ae911d..3f31670 100644
--- a/src/resume.f90
+++ b/src/resume.f90
@@ -1,80 +1,83 @@
 SUBROUTINE resume
 !
 !   Resume from previous run
 !
   Use beam
   Use basic
   Use fields
   IMPLICIT NONE
  
 !
 !   Local vars and arrays
  INTEGER, DIMENSION(:), ALLOCATABLE:: displs,sendcounts
  INTEGER:: i
 !________________________________________________________________________________
   WRITE(*,'(a)') '   Resume from previous run'
 !________________________________________________________________________________
 ! 
 !   Open and read initial conditions from restart file
-  CALL chkrst(0)  
+  CALL chkrst(0)
+  IF (newres) THEN
+    cstep=0
+  END IF  
   IF(mpisize .gt. 1) THEN
     ALLOCATE(sendcounts(0:mpisize-1), displs(0:mpisize-1)) 
     IF(mpirank.eq.0) THEN
       CALL quicksortparts(1,parts%Nptot)
       CALL distribpartsonprocs
       CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER,MPI_IN_PLACE, 1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       DO i=0,mpisize-1
         displs(i)=i
         sendcounts(i)=2
       END DO
       CALL MPI_Scatterv(Zbounds,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,2,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       displs(0)=0
       sendcounts(0)=Nplocs_all(0)
       DO i=1,mpisize-1
         displs(i)=displs(i-1)+Nplocs_all(i-1)
         sendcounts(i)=Nplocs_all(i)
       END DO
       CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
       CALL MPI_Scatterv(parts%R,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%Z,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%THET,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%UZ,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%UR,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%Rindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%Zindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%partindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
     ELSE
       CALL creat_parts(parts,nplasma)
       CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER, Nplocs_all(mpirank), 1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       sendcounts=2
       CALL MPI_Scatterv(Zbounds,sendcounts,displs,MPI_DOUBLE_PRECISION,Zbounds(mpirank),2,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       parts%Nptot=Nplocs_all(mpirank)
 
       WRITE(*,*) mpirank, "received Nplocs_all : ", Nplocs_all
       WRITE(*,*) mpirank, "received Zbounds : ", Zbounds
 
       CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
       CALL MPI_Scatterv(parts%R,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%R,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%Z,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%Z,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%THET,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%THET,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%UZ,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UZ,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%UR,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UR,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UTHET,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%Rindex,sendcounts, displs,MPI_INTEGER,parts%Rindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%Zindex,sendcounts, displs,MPI_INTEGER,parts%Zindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       CALL MPI_Scatterv(parts%partindex,sendcounts, displs,MPI_INTEGER,parts%partindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
     END IF
     CALL MPI_Bcast(V,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
     CALL MPI_Bcast(qsim,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
     CALL MPI_Bcast(msim,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
   END IF
 ! Compute Self consistent electric field
   CALL bound
   CALL localisation
   ! Init Electric and Magnetic Fields
   CALL fields_init
   CALL rhscon
   CALL poisson
 !
 END SUBROUTINE resume