Index: options/COMMAND
===================================================================
--- options/COMMAND	(revision 6d73c4be98f6a308eed735650a968ac230fd19ef)
+++ options/COMMAND	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -19,5 +19,5 @@
  IFINE=                 4, ! Reduction for time step in vertical transport, used only if CTL>1 
  IOUT=                  1, ! Output type: [1]mass 2]pptv 3]1&2 4]plume 5]1&4, +8 for NetCDF output     
- IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end   
+ IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end 3]time averaged
  LSUBGRID=              0, ! Increase of ABL heights due to sub-grid scale orographic variations;[0]off 1]on 
  LCONVECTION=          	1, ! Switch for convection parameterization;0]off [1]on    
Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision 50958b85292a7b2ce570403ca6f859592eaa1b5e)
+++ src/FLEXPART.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -68,11 +68,5 @@
   integer :: detectformat
 
-
-
-  ! Initialize arrays in com_mod
-  !*****************************
-  call com_mod_allocate_part(maxpart)
   
-
   ! Generate a large number of random numbers
   !******************************************
@@ -172,4 +166,9 @@
   endif
 
+  ! Initialize arrays in com_mod
+  !*****************************
+  call com_mod_allocate_part(maxpart)
+
+
   ! Read the age classes to be used
   !********************************
Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
+++ src/FLEXPART_MPI.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -77,10 +77,5 @@
   if (mp_measure_time) call mpif_mtime('flexpart',0)
 
-  ! Initialize arrays in com_mod 
-  !*****************************
-
-  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
-
-
+  
   ! Generate a large number of random numbers
   !******************************************
@@ -180,4 +175,9 @@
   endif
 
+  ! Initialize arrays in com_mod 
+  !*****************************
+
+  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
+
 
 ! Read the age classes to be used
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 328fdf90784e3387eb9745d1009d7523f4ed1976)
+++ src/com_mod.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -18,4 +18,6 @@
 
   implicit none
+
+
 
   !****************************************************************
@@ -69,5 +71,5 @@
 
   real :: ctl,fine
-  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
+  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill,ipoutfac
   integer :: mquasilag,nested_output,ind_source,ind_receptor
   integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
@@ -82,4 +84,5 @@
   ! iout     output options: 1 conc. output (ng/m3), 2 mixing ratio (pptv), 3 both
   ! ipout    particle dump options: 0 no, 1 every output interval, 2 only at end
+  ! ipoutfac increase particle dump interval by factor (default 1)
   ! ipin     read in particle positions from dumped file from a previous run
   ! fine     real(ifine)
@@ -128,5 +131,4 @@
 
   logical :: gdomainfill
-
   ! gdomainfill             .T., if domain-filling is global, .F. if not
 
@@ -675,4 +677,12 @@
   real, allocatable, dimension(:,:) :: xscav_frac1
 
+! Variables used for writing out interval averages for partoutput
+!****************************************************************
+
+  integer, allocatable, dimension(:) :: npart_av
+  real, allocatable, dimension(:) :: part_av_cartx,part_av_carty,part_av_cartz,part_av_z,part_av_topo
+  real, allocatable, dimension(:) :: part_av_pv,part_av_qv,part_av_tt,part_av_rho,part_av_tro,part_av_hmix
+  real, allocatable, dimension(:) :: part_av_uu,part_av_vv,part_av_energy
+
   ! eso: Moved from timemanager
   real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
@@ -781,11 +791,19 @@
          & idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
          & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
-         & xmass1(nmpart, maxspec),&
-         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
+         & xmass1(nmpart, maxspec))  ! ,&
+!         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
+
+    if (ipout.eq.3) then
+      allocate(npart_av(nmpart),part_av_cartx(nmpart),part_av_carty(nmpart),&
+           & part_av_cartz(nmpart),part_av_z(nmpart),part_av_topo(nmpart))
+      allocate(part_av_pv(nmpart),part_av_qv(nmpart),part_av_tt(nmpart),&
+           & part_av_rho(nmpart),part_av_tro(nmpart),part_av_hmix(nmpart))
+      allocate(part_av_uu(nmpart),part_av_vv(nmpart),part_av_energy(nmpart))
+    end if
 
 
     allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
          & vs(nmpart),ws(nmpart),cbt(nmpart))
-    
+
   end subroutine com_mod_allocate_part
 
Index: src/init_domainfill.f90
===================================================================
--- src/init_domainfill.f90	(revision 328fdf90784e3387eb9745d1009d7523f4ed1976)
+++ src/init_domainfill.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -86,4 +86,8 @@
     endif
   endif
+
+! Exit here if resuming a run from particle dump
+!***********************************************
+  if (gdomainfill.and.ipin.ne.0) return
 
 ! Do not release particles twice (i.e., not at both in the leftmost and rightmost
Index: src/makefile
===================================================================
--- src/makefile	(revision 7123c700d7191876c1d22402ad01729c688c6e56)
+++ src/makefile	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -118,4 +118,5 @@
 OBJECTS_SERIAL = \
 	releaseparticles.o 	partoutput.o \
+	partoutput_average.o \
 	conccalc.o \
 	init_domainfill.o 	concoutput.o  \
@@ -132,5 +133,5 @@
 ## For MPI version 
 OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
-	conccalc_mpi.o \
+	partoutput_average_mpi.o conccalc_mpi.o \
 	init_domainfill_mpi.o concoutput_mpi.o 	\
 	timemanager_mpi.o FLEXPART_MPI.o	\
@@ -149,5 +150,5 @@
 advance.o		initialize.o		\
 writeheader.o		writeheader_txt.o  	\
-writeprecip.o \
+partpos_average.o	writeprecip.o \
 writeheader_surf.o     	assignland.o\
 part0.o 		gethourlyOH.o\
@@ -348,5 +349,8 @@
 part0.o: par_mod.o
 partdep.o: par_mod.o
+partpos_average.o: com_mod.o par_mod.o
 partoutput.o: com_mod.o par_mod.o
+partoutput_average.o: com_mod.o par_mod.o
+partoutput_average_mpi.o: com_mod.o par_mod.o mpi_mod.o
 partoutput_mpi.o: com_mod.o mpi_mod.o par_mod.o
 partoutput_short.o: com_mod.o par_mod.o
Index: src/mpi_mod.f90
===================================================================
--- src/mpi_mod.f90	(revision 328fdf90784e3387eb9745d1009d7523f4ed1976)
+++ src/mpi_mod.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -88,4 +88,5 @@
 ! Variables for MPI processes in the 'particle' group
   integer, allocatable, dimension(:) :: mp_partgroup_rank
+  integer, allocatable, dimension(:) :: npart_per_process
   integer :: mp_partgroup_comm, mp_partgroup_pid, mp_partgroup_np
 
@@ -125,6 +126,5 @@
 ! mp_time_barrier   Measure MPI barrier time
 ! mp_exact_numpart  Use an extra MPI communication to give the exact number of particles
-!                   to standard output (this does *not* otherwise affect the simulation)
-! mp_rebalance      Attempt to rebalance particle between processes
+!                   to standard output (this does not otherwise affect the simulation)
   logical, parameter :: mp_dbg_mode = .false.
   logical, parameter :: mp_dev_mode = .false.
@@ -133,5 +133,4 @@
   logical, parameter :: mp_measure_time=.false.
   logical, parameter :: mp_exact_numpart=.true.
-  logical, parameter :: mp_rebalance=.true.
 
 ! for measuring CPU/Wall time
@@ -146,6 +145,6 @@
   real(dp),private :: mp_getfields_wtime_beg, mp_getfields_wtime_end, mp_getfields_wtime_total=0.
   real(sp),private :: mp_getfields_time_beg, mp_getfields_time_end, mp_getfields_time_total=0.
-!  real(dp),private :: mp_readwind_wtime_beg, mp_readwind_wtime_end, mp_readwind_wtime_total=0.
-!  real(sp),private :: mp_readwind_time_beg, mp_readwind_time_end, mp_readwind_time_total=0.
+  real(dp),private :: mp_readwind_wtime_beg, mp_readwind_wtime_end, mp_readwind_wtime_total=0.
+  real(sp),private :: mp_readwind_time_beg, mp_readwind_time_end, mp_readwind_time_total=0.
   real(dp),private :: mp_io_wtime_beg, mp_io_wtime_end, mp_io_wtime_total=0.
   real(sp),private :: mp_io_time_beg, mp_io_time_end, mp_io_time_total=0.
@@ -367,4 +366,7 @@
     end if
 
+! Allocate array for number of particles per process    
+    allocate(npart_per_process(0:mp_partgroup_np-1))
+
 ! Write whether MPI_IN_PLACE is used or not
 #ifdef USE_MPIINPLACE
@@ -607,5 +609,5 @@
     integer :: i,jj,nn, num_part=1,m,imin, num_trans
     logical :: first_iter
-    integer,allocatable,dimension(:) :: numparticles_mpi, idx_arr
+    integer,allocatable,dimension(:) :: idx_arr
     real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this
 
@@ -616,8 +618,7 @@
 ! All processes exchange information on number of particles
 !****************************************************************************
-    allocate(numparticles_mpi(0:mp_partgroup_np-1), &
-         &idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
-
-    call MPI_Allgather(numpart, 1, MPI_INTEGER, numparticles_mpi, &
+    allocate( idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
+
+    call MPI_Allgather(numpart, 1, MPI_INTEGER, npart_per_process, &
          & 1, MPI_INTEGER, mp_comm_used, mp_ierr)
 
@@ -625,8 +626,11 @@
 ! Sort from lowest to highest
 ! Initial guess: correct order
-    sorted(:) = numparticles_mpi(:)
+    sorted(:) = npart_per_process(:)
     do i=0, mp_partgroup_np-1
       idx_arr(i) = i
     end do
+
+! Do not rebalance particles for ipout=3    
+    if (ipout.eq.3) return
 
 ! For each successive element in index array, see if a lower value exists
@@ -654,11 +658,11 @@
     m=mp_partgroup_np-1 ! index for last sorted process (most particles)
     do i=0,mp_partgroup_np/2-1
-      num_trans = numparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))
+      num_trans = npart_per_process(idx_arr(m)) - npart_per_process(idx_arr(i))
       if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then
-        if ( numparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.&
-             & real(num_trans)/real(numparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then
+        if ( npart_per_process(idx_arr(m)).gt.mp_min_redist.and.&
+             & real(num_trans)/real(npart_per_process(idx_arr(m))).gt.mp_redist_fract) then
 ! DBG
-          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi', &
-          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi
+          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process', &
+          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process
 ! DBG
           call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)
@@ -668,5 +672,5 @@
     end do
 
-    deallocate(numparticles_mpi, idx_arr, sorted)
+    deallocate(idx_arr, sorted)
 
   end subroutine mpif_calculate_part_redist
@@ -2716,17 +2720,17 @@
       end if
 
-!    case ('readwind')
-!      if (imode.eq.0) then
-!        call cpu_time(mp_readwind_time_beg)
-!        mp_readwind_wtime_beg = mpi_wtime()
-!      else
-!        call cpu_time(mp_readwind_time_end)
-!        mp_readwind_wtime_end = mpi_wtime()
-!
-!        mp_readwind_time_total = mp_readwind_time_total + &
-!             &(mp_readwind_time_end - mp_readwind_time_beg)
-!        mp_readwind_wtime_total = mp_readwind_wtime_total + &
-!             &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
-!      end if
+   case ('readwind')
+     if (imode.eq.0) then
+       call cpu_time(mp_readwind_time_beg)
+       mp_readwind_wtime_beg = mpi_wtime()
+     else
+       call cpu_time(mp_readwind_time_end)
+       mp_readwind_wtime_end = mpi_wtime()
+
+       mp_readwind_time_total = mp_readwind_time_total + &
+            &(mp_readwind_time_end - mp_readwind_time_beg)
+       mp_readwind_wtime_total = mp_readwind_wtime_total + &
+            &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
+     end if
 
     case ('commtime')
Index: src/netcdf_output_mod.f90
===================================================================
--- src/netcdf_output_mod.f90	(revision 4ad96c54b793af523a98633b88f07aa60c69ea47)
+++ src/netcdf_output_mod.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -93,5 +93,5 @@
   character(len=255), parameter :: institution = 'NILU'
 
-  integer            :: tpointer
+  integer            :: tpointer=0
   character(len=255) :: ncfname, ncfnamen
 
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 79e0349a50fd1ee6dfdf4b044ade2b6f1e702855)
+++ src/par_mod.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -280,5 +280,5 @@
 
   integer,parameter :: unitpath=1, unitcommand=1, unitageclasses=1, unitgrid=1
-  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93
+  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93, unitpartout_average=105
   integer,parameter :: unitpartin=93, unitflux=98, unitouttraj=96
   integer,parameter :: unitvert=1, unitoro=1, unitpoin=1, unitreceptor=1
Index: src/partoutput.f90
===================================================================
--- src/partoutput.f90	(revision d2a5a833066965c4c4b33df109fd4e461f7275bb)
+++ src/partoutput.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -71,5 +71,5 @@
   !**************************************
 
-  if (ipout.eq.1) then
+  if (ipout.eq.1.or.ipout.eq.3) then
     open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
          atime,form='unformatted')
Index: src/partoutput_mpi.f90
===================================================================
--- src/partoutput_mpi.f90	(revision d2a5a833066965c4c4b33df109fd4e461f7275bb)
+++ src/partoutput_mpi.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -78,5 +78,5 @@
   !**************************************
 
-  if (ipout.eq.1) then
+  if (ipout.eq.1.or.ipout.eq.3) then
     open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
          atime,form='unformatted',status=file_stat,position='append')
Index: src/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
+++ src/readcommand.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -50,4 +50,5 @@
   ! ipin                 1 continue simulation with dumped particle data, 0 no *
   ! ipout                0 no particle dump, 1 every output time, 3 only at end*
+  ! ipoutfac             increase particle dump interval by factor (default 1) *
   ! itsplit [s]          time constant for particle splitting                  *
   ! loutaver [s]         concentration output is an average over loutaver      *
@@ -97,4 +98,5 @@
   iout, &
   ipout, &
+  ipoutfac, &
   lsubgrid, &
   lconvection, &
@@ -129,4 +131,5 @@
   iout=3
   ipout=0
+  ipoutfac=1
   lsubgrid=1
   lconvection=1
@@ -507,7 +510,7 @@
   !****************************************************************
 
-  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
+  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2).and.(ipout.ne.3)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
-    write(*,*) ' #### IPOUT MUST BE 1, 2 OR 3!                #### '
+    write(*,*) ' #### IPOUT MUST BE 0, 1, 2 OR 3!             #### '
     stop
   endif
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision c7d1052e9b4a8319033f957f6b77184654bca402)
+++ src/timemanager.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -451,5 +451,8 @@
 45      format(i13,' Seconds simulated: ',i13, ' Particles:    Uncertainty: ',3f7.3)
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
-        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
+        if (ipout.ge.1) then
+          if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
+          if (ipout.eq.3) call partoutput_average(itime) ! dump particle positions
+        endif
         loutnext=loutnext+loutstep
         loutstart=loutnext-loutaver/2
@@ -609,4 +612,10 @@
 !        write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
 
+  ! Calculate average position for particle dump output
+  !****************************************************
+
+        if (ipout.eq.3) call partpos_average(itime,j)
+
+
   ! Calculate the gross fluxes across layer interfaces
   !***************************************************
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision 328fdf90784e3387eb9745d1009d7523f4ed1976)
+++ src/timemanager_mpi.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -113,5 +113,5 @@
   integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
 ! integer :: ksp
-  integer :: ip
+  integer :: ip,irec
   integer :: loutnext,loutstart,loutend
   integer :: ix,jy,ldeltat,itage,nage,idummy
@@ -129,4 +129,5 @@
 ! Measure time spent in timemanager
   if (mp_measure_time) call mpif_mtime('timemanager',0)
+
 
 ! First output for time 0
@@ -344,5 +345,5 @@
 ! Check if particles should be redistributed among processes
 !***********************************************************
-    if (mp_rebalance) call mpif_calculate_part_redist(itime)
+    call mpif_calculate_part_redist(itime)
 
 
@@ -593,6 +594,11 @@
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
         if (ipout.ge.1) then
+          irec=0
           do ip=0, mp_partgroup_np-1
-            if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions
+            if (ip.eq.mp_partid) then
+              if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
+              if (ipout.eq.3) call partoutput_average(itime,irec) ! dump particle positions
+            endif
+            if (ipout.eq.3) irec=irec+npart_per_process(ip)
             call mpif_mpi_barrier
           end do
@@ -757,4 +763,9 @@
         if (mp_measure_time) call mpif_mtime('advance',1)
 
+  ! Calculate average position for particle dump output
+  !****************************************************
+
+        if (ipout.eq.3) call partpos_average(itime,j)
+
 
 ! Calculate the gross fluxes across layer interfaces
@@ -895,5 +906,5 @@
     do ip=0, mp_partgroup_np-1
       if (ip.eq.mp_partid) then 
-        !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
+        if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
         call partoutput(itime)    ! dump particle positions
       end if
Index: src/verttransform_ecmwf.f90
===================================================================
--- src/verttransform_ecmwf.f90	(revision 79e0349a50fd1ee6dfdf4b044ade2b6f1e702855)
+++ src/verttransform_ecmwf.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
@@ -73,5 +73,4 @@
   use com_mod
   use cmapf_mod, only: cc2gll
-!  use mpi_mod
 
   implicit none
