diff --git a/src/auxval.f90 b/src/auxval.f90
index 8bd088c..1a16a5a 100644
--- a/src/auxval.f90
+++ b/src/auxval.f90
@@ -1,150 +1,150 @@
 SUBROUTINE auxval
 !
   USE constants
   USE basic
   USE fields
   USE beam
   USE random
   USE omp_lib
 !
 !   Set auxiliary values
 !
   IMPLICIT NONE
 !
 
   INTEGER:: i, nbprocs
   
 !________________________________________________________________________________
 IF(mpirank .eq. 0) WRITE(*,'(a/)') '=== Set auxiliary values ==='
 !________________________________________________________________________________
 !
 ! Total number of intervals
   nr=sum(nnr)
 
 ! Normalisation constants
   if(nplasma .gt. 0) then
     qsim=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)*n0*elchar/nplasma 
   else
     qsim=sign(n0,elchar)
   end if
   msim=abs(qsim)/elchar*partmass
   vnorm=vlight
 
   omegac=sign(elchar,qsim)/partmass*B0
   omegap=sqrt(elchar**2*abs(n0)/partmass/eps_0)
 
   tnorm=min(abs(1/omegac),abs(1/omegap))
   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 
   
   IF(mpirank .eq. 0) THEN
     IF(abs(omegap).GT. abs(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
   invdr=1/dr
   invdz=1/dz
 
 
   ! Define Zindex bounds for the MPI process in the case of sequential run
   ALLOCATE( Zbounds(0:mpisize), Nplocs_all(0:mpisize-1))
   Zbounds=0
   Zbounds(mpisize) = nz
   IF(mpisize.gt.1) THEN
     ! If we use mpi we define the left and right nodes used for communications
     leftproc=mpirank-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 .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
 
 ! Initialize random number generator
 nbprocs = omp_get_max_threads()
 allocate(seed(ran_s,nbprocs), ran_index(nbprocs), ran_array(ran_k,nbprocs))
 
 Do i=1,nbprocs
 ! Generate seed from the default seed-string in random module
   CALL decimal_to_seed(random_seed_str, seed(:,i))
 
 ! Generate a different seed for each processor from the mother seed
 
   CALL next_seed(mpirank*nbprocs+i,seed(:,i))
 
 ! Initialize the random array (first hundred numbers)
 
   CALL random_init(seed(:,i), ran_index(i), ran_array(:,i))
 end do
 !________________________________________________________________________________
 CONTAINS 
 !---------------------------------------------------------------------------
 !> @author
 !> Patryk Kaminski   EPFL/SPC
 !
 ! DESCRIPTION:
 !>
 !> @brief Creates the mesh in r and z direction for calculating the electric and magnetic fields.
 !---------------------------------------------------------------------------
   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
     if (nnr(2).gt.0) then
         dr(2)=(radii(3)-radii(2))/nnr(2)
         DO j=1,nnr(2)
            rgrid(nnr(1)+j)=radii(2)+j*dr(2)
         END DO
     end if
-    if (nnr(2).gt.0) then
+    if (nnr(3).gt.0) then
         dr(3)=(radii(4)-radii(3))/nnr(3)
-        DO j=1,nnr(3)-1
+        DO j=1,nnr(3)
            rgrid(nnr(1)+nnr(2)+j)=radii(3)+j*dr(3)
         END DO
     end if
   END SUBROUTINE mesh
 !________________________________________________________________________________
 
 END SUBROUTINE auxval