Index: src/conccalc.f90
===================================================================
--- src/conccalc.f90	(revision e52967c0a5f99e764feb376b6267ce20286a2850)
+++ src/conccalc.f90	(revision 08a38b55d3749abb3d1df461e12c98d18b329a51)
@@ -193,15 +193,22 @@
              (jy.le.numygrid-1)) then
           if (DRYBKDEP.or.WETBKDEP) then
-             do ks=1,nspec
-               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
-                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
-                 xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
-             end do
+            do ks=1,nspec
+              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
+            end do
           else
-             do ks=1,nspec
-               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
-                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
-                 xmass1(i,ks)/rhoi*weight
-             end do
+            if (lparticlecountoutput) then
+              do ks=1,nspec
+                gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
+              end do
+            else
+              do ks=1,nspec
+                gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight
+              end do
+            end if
           endif
         endif
@@ -336,12 +343,19 @@
                end do
             else
-               do ks=1,nspec
-                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
-                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
-                   xmass1(i,ks)/rhoi*weight
-               end do
-            endif
-          endif
-
+              if (lparticlecountoutput) then
+                do ks=1,nspec
+                  griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                       griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
+                end do
+              else
+                do ks=1,nspec
+                  griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                       griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                       xmass1(i,ks)/rhoi*weight
+                end do
+              endif
+            endif
+          endif
+          
         else                                 ! attribution via uniform kernel
 
Index: src/conccalc_mpi.f90
===================================================================
--- src/conccalc_mpi.f90	(revision e52967c0a5f99e764feb376b6267ce20286a2850)
+++ src/conccalc_mpi.f90	(revision 08a38b55d3749abb3d1df461e12c98d18b329a51)
@@ -205,9 +205,24 @@
         if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
              (jy.le.numygrid-1)) then
-          do ks=1,nspec
-            gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
-                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
-                 xmass1(i,ks)/rhoi*weight
-          end do
+          if (DRYBKDEP.or.WETBKDEP) then
+            do ks=1,nspec
+              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
+            end do
+          else
+            if (lparticlecountoutput) then
+              do ks=1,nspec
+                gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
+              end do
+            else
+              do ks=1,nspec
+                gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight
+              end do
+            end if
+          endif
         endif
 
@@ -239,18 +254,34 @@
           if ((jy.ge.0).and.(jy.le.numygrid-1)) then
             w=wx*wy
-            do ks=1,nspec
-              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+               do ks=1,nspec
+                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+            endif
           endif
 
           if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
             w=wx*(1.-wy)
-            do ks=1,nspec
-              gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+              do ks=1,nspec
+                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+               end do
+             else
+              do ks=1,nspec
+                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+             endif
           endif
         endif
@@ -260,18 +291,34 @@
           if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
             w=(1.-wx)*(1.-wy)
-            do ks=1,nspec
-              gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+               do ks=1,nspec
+                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+            endif
           endif
 
           if ((jy.ge.0).and.(jy.le.numygrid-1)) then
             w=(1.-wx)*wy
-            do ks=1,nspec
-              gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+               do ks=1,nspec
+                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+            endif
           endif
         endif
@@ -304,9 +351,24 @@
           if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
                (jy.le.numygridn-1)) then
-            do ks=1,nspec
-              griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                    griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
-                   xmass1(i,ks)/rhoi*weight
-            end do
+                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+              if (lparticlecountoutput) then
+                do ks=1,nspec
+                  griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                       griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
+                end do
+              else
+                do ks=1,nspec
+                  griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                       griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                       xmass1(i,ks)/rhoi*weight
+                end do
+              endif
+            endif
           endif
 
@@ -338,18 +400,34 @@
             if ((jy.ge.0).and.(jy.le.numygridn-1)) then
               w=wx*wy
-              do ks=1,nspec
-                griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                do ks=1,nspec
+                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
 
             if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
               w=wx*(1.-wy)
-              do ks=1,nspec
-                griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                 do ks=1,nspec
+                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
           endif
@@ -359,18 +437,34 @@
             if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
               w=(1.-wx)*(1.-wy)
-              do ks=1,nspec
-                griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                 do ks=1,nspec
+                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
 
             if ((jy.ge.0).and.(jy.le.numygridn-1)) then
               w=(1.-wx)*wy
-              do ks=1,nspec
-                griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                 do ks=1,nspec
+                    griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
           endif
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision 16b61a5e2d3b1020be1cff080e479ba26f3f7063)
+++ src/concoutput.f90	(revision 08a38b55d3749abb3d1df461e12c98d18b329a51)
@@ -446,8 +446,16 @@
                   endif
                   sp_count_r=sp_count_r+1
-                  sparse_dump_r(sp_count_r)= &
-                       sp_fact* &
-                       grid(ix,jy,kz)* &
-                       factor3d(ix,jy,kz)/tot_mu(ks,kp)
+                  if (lparticlecountoutput) then
+                    sparse_dump_r(sp_count_r)= &
+                         sp_fact* &
+                         grid(ix,jy,kz)
+                  else
+                    sparse_dump_r(sp_count_r)= &
+                         sp_fact* &
+                         grid(ix,jy,kz)* &
+                         factor3d(ix,jy,kz)/tot_mu(ks,kp)
+                  end if
+
+
 !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
 !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision d8eed026ebec594e9a91de6c6fbc45b354a0e4b2)
+++ src/par_mod.f90	(revision 08a38b55d3749abb3d1df461e12c98d18b329a51)
@@ -60,6 +60,11 @@
   logical, parameter :: lnokernel=.false.
 
+  !*********************************************************************
+  ! Set to T to change output units to number of particles per grid cell
+  !*********************************************************************
+  logical, parameter :: lparticlecountoutput=.false.
+
   !***********************************************************
-  ! Number of directories/files used for FLEXPART input/output
+  ! number of directories/files used for FLEXPART input/output
   !***********************************************************
 
@@ -141,13 +146,13 @@
   !*********************************************
   
-!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new 
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
-
-!  integer,parameter :: nxmax=181,nymax=91,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
-
-!  INTEGER,PARAMETER :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !NCEP data
-
-!  integer,parameter :: nxshift=359  ! for ECMWF
-  integer,parameter :: nxshift=0     ! for GFS
+  ! ECMWF
+! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=359 ! 1.0 degree 92 level
+   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359 ! 1.0 degree 138 level
+! integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359  ! 0.5 degree 138 level
+!  integer,parameter :: nxmax=181,nymax=91,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=0  ! CERA 2.0 degree 92 level
+
+! GFS
+!   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0
+
 
   !*********************************************
@@ -155,5 +160,5 @@
   !*********************************************
 
-  integer,parameter :: maxnests=1,nxmaxn=451,nymaxn=226
+  integer,parameter :: maxnests=0,nxmaxn=451,nymaxn=226
 
   ! nxmax,nymax        maximum dimension of wind fields in x and y
@@ -213,5 +218,5 @@
   !**************************************************
 
-  integer,parameter :: maxpart=10000000
+  integer,parameter :: maxpart=3000000
   integer,parameter :: maxspec=4
 
