Index: src/advance.f90
===================================================================
--- src/advance.f90	(revision 18adf6049d9d47356934a37eb7784506f795e220)
+++ src/advance.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
@@ -493,13 +493,12 @@
 
       if (mdomainfill.eq.0) then
-! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
-!              particle number and thus xmass array goes out of bounds
-!        do nsp=1,nspec
-!           if (xmass(nrelpoint,nsp).gt.eps2) goto 887
-!         end do
-! 887     nsp=min(nsp,nspec)
-        if (nspec.eq.1.and.density(1).gt.0.) then
-          call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
-        end if
+        do nsp=1,nspec
+          if (xmass(nrelpoint,nsp).gt.eps2) goto 887
+        end do
+887     nsp=min(nsp,nspec)
+!!$        if (density(nsp).gt.0.) &
+!!$             call get_settling(itime,xts,yts,zt,nsp,settling)    !old
+        if (density(nsp).gt.0.) &
+             call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
         w=w+settling
       endif
@@ -655,15 +654,12 @@
 
     if (mdomainfill.eq.0) then
-! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
-!              particle number and thus xmass array goes out of bounds
-
-!      do nsp=1,nspec
-!         if (xmass(nrelpoint,nsp).gt.eps2) goto 888
-!       end do
-! 888   nsp=min(nsp,nspec)
-!        if (density(nsp).gt.0.) then
-      if (nspec.eq.1.and.density(1).gt.0.) then
-        call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
-      end if
+      do nsp=1,nspec
+        if (xmass(nrelpoint,nsp).gt.eps2) goto 888
+      end do
+888   nsp=min(nsp,nspec)
+!!$      if (density(nsp).gt.0.) &
+!!$           call get_settling(itime,xts,yts,zt,nsp,settling)    !old
+      if (density(nsp).gt.0.) &
+           call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
       w=w+settling
     endif
@@ -862,14 +858,12 @@
 
   if (mdomainfill.eq.0) then
-! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
-!              particle number and thus xmass array goes out of bounds
-!    do nsp=1,nspec
-!       if (xmass(nrelpoint,nsp).gt.eps2) goto 889
-!     end do
-! 889   nsp=min(nsp,nspec)
-!      if (density(nsp).gt.0.) then
-    if (nspec.eq.1.and.density(1).gt.0.) then
-      call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling)  !bugfix
-    end if
+    do nsp=1,nspec
+      if (xmass(nrelpoint,nsp).gt.eps2) goto 889
+    end do
+889   nsp=min(nsp,nspec)
+!!$    if (density(nsp).gt.0.) &
+!!$         call get_settling(itime+ldt,xts,yts,zt,nsp,settling)    !old
+    if (density(nsp).gt.0.) &
+         call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling)  !bugfix
     w=w+settling
   endif
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 18adf6049d9d47356934a37eb7784506f795e220)
+++ src/par_mod.f90	(revision 62e65c73ae7e3d8c4b8b884065ede76e07211ec8)
@@ -186,5 +186,5 @@
   !**************************************************
 
-  integer,parameter :: maxpart=1000000
+  integer,parameter :: maxpart=20000000
   integer,parameter :: maxspec=6
   real,parameter :: minmass=0.0001
Index: src/partoutput_short.f90
===================================================================
--- src/partoutput_short.f90	(revision 18adf6049d9d47356934a37eb7784506f795e220)
+++ src/partoutput_short.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
@@ -146,5 +146,5 @@
 
 
-  !write(*,*) numshortout,numshortall
+  write(*,*) numshortout,numshortall
 
   close(unitshortpart)
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision 18adf6049d9d47356934a37eb7784506f795e220)
+++ src/timemanager.f90	(revision 05cf28de9ff073ff7e3cdc064f0ed59c61a6dd09)
@@ -589,6 +589,6 @@
             endif
 
-! Skip check on mass fraction when npoint represents particle number
-            if (mdomainfill.eq.0.and.mquasilag.eq.0) then
+
+            if (mdomainfill.eq.0) then
               if (xmass(npoint(j),ks).gt.0.) &
                    xmassfract=max(xmassfract,real(npart(npoint(j)))* &
@@ -605,6 +605,6 @@
 !ZHG 2015
             else
-              xmassfract=1.0
-            end if
+              xmassfract=1.
+            endif
           end do
 
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision 18adf6049d9d47356934a37eb7784506f795e220)
+++ src/timemanager_mpi.f90	(revision 05cf28de9ff073ff7e3cdc064f0ed59c61a6dd09)
@@ -546,12 +546,7 @@
 ! Decide whether to write an estimate of the number of particles released, 
 ! or exact number (require MPI reduce operation)
-        if (mp_dev_mode) then
-          numpart_tot_mpi = numpart
-        else
-          numpart_tot_mpi = numpart*mp_partgroup_np
-        end if
-
-        if (mp_exact_numpart.and..not.(lmpreader.and.lmp_use_reader).and.&
-             &.not.mp_dev_mode) then
+        numpart_tot_mpi = numpart*mp_partgroup_np
+
+        if (mp_exact_numpart.and..not.(lmpreader.and.lmp_use_reader)) then
           call MPI_Reduce(numpart, numpart_tot_mpi, 1, MPI_INTEGER, MPI_SUM, id_root, &
                & mp_comm_used, mp_ierr)
@@ -743,6 +738,6 @@
             endif
 
-! Skip check on mass fraction when npoint represents particle number
-            if (mdomainfill.eq.0.and.mquasilag.eq.0) then
+
+            if (mdomainfill.eq.0) then
               if (xmass(npoint(j),ks).gt.0.)then 
                    xmassfract=max(xmassfract,real(npart(npoint(j)))* &
