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 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
@@ -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
@@ -413,5 +413,5 @@
   end if ! (mpif_pid == 0) 
 
-  if (mp_measure_time) call mpif_mtime('iotime',0)
+  if (mp_measure_time) call mpif_mtime('iotime',1)
 
   if (verbosity.gt.0 .and. lroot) then
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision e9e0f0627da3a042c02a5352bf3446bdfbcfa136)
+++ 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
 
@@ -651,4 +653,5 @@
   real :: receptorarea(maxreceptor)
   real :: creceptor(maxreceptor,maxspec)
+  real, allocatable, dimension(:,:) :: creceptor0
   character(len=16) :: receptorname(maxreceptor)
   integer :: numreceptor
@@ -673,4 +676,12 @@
   real, allocatable, dimension(:,:) :: xmass1
   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
@@ -780,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 b5127f9034f2a15f9a3ce226fd26edfe845311de)
+++ 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
@@ -414,5 +418,5 @@
 !***************************************************************************
 
-  if (ipin.eq.1) then
+  if ((ipin.eq.1).and.(.not.gdomainfill)) then
     open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
          form='unformatted')
Index: src/init_domainfill_mpi.f90
===================================================================
--- src/init_domainfill_mpi.f90	(revision b5127f9034f2a15f9a3ce226fd26edfe845311de)
+++ src/init_domainfill_mpi.f90	(revision 328fdf90784e3387eb9745d1009d7523f4ed1976)
@@ -110,4 +110,8 @@
   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
 ! grid cell) for a global domain
@@ -213,5 +217,4 @@
         colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)
         colmasstotal=colmasstotal+colmass(ix,jy)
-
       end do
     end do
@@ -466,5 +469,5 @@
 
 ! eso TODO: only needed for root process
-    if (ipin.eq.1) then
+    if ((ipin.eq.1).and.(.not.gdomainfill)) then
       open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
            form='unformatted')
@@ -474,25 +477,31 @@
     endif
 
-    numpart = numpart/mp_partgroup_np
-    if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1
-
-  else ! Allocate dummy arrays for receiving processes 
-    allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),&
-         & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),&
-         & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),&
-         & xmass1_tmp(nullsize, nullsize))
+    if (ipin.eq.0) then    
+      numpart = numpart/mp_partgroup_np
+      if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1
+    end if
+
+  else ! Allocate dummy arrays for receiving processes
+    if (ipin.eq.0) then    
+      allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),&
+           & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),&
+           & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),&
+           & xmass1_tmp(nullsize, nullsize))
+    end if
     
-  end if ! end if(lroot)  
+  end if ! end if(lroot)
+
 
 
 ! Distribute particles to other processes (numpart is 'per-process', not total)
-  call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
-! eso TODO: xmassperparticle: not necessary to send
-  call MPI_Bcast(xmassperparticle, 1, mp_sp, id_root, mp_comm_used, mp_ierr)
-  call mpif_send_part_properties(numpart)
+! Only if not restarting from previous run
+  if (ipin.eq.0) then
+    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
+    call mpif_send_part_properties(npart(1)/mp_partgroup_np)
 
 ! Deallocate the temporary arrays used for all particles
-  deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
+    deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
          & itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)
+  end if
 
 
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 0ecc1fed9f7c7c07cef687094373915d51c69489)
+++ src/mpi_mod.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
@@ -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,5 +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) 
+!                   to standard output (this does not otherwise affect the simulation)
   logical, parameter :: mp_dbg_mode = .false.
   logical, parameter :: mp_dev_mode = .false.
@@ -190,6 +191,6 @@
 !   mp_np       number of running processes, decided at run-time
 !***********************************************************************
-    use par_mod, only: maxpart, numwfmem, dep_prec
-    use com_mod, only: mpi_mode, verbosity
+    use par_mod, only: maxpart, numwfmem, dep_prec, maxreceptor, maxspec
+    use com_mod, only: mpi_mode, verbosity, creceptor0
 
     implicit none
@@ -337,5 +338,5 @@
 
 ! Set maxpart per process
-! eso 08/2016: Increase maxpart per process, in case of unbalanced distribution
+! ESO 08/2016: Increase maxpart per process, in case of unbalanced distribution
     maxpart_mpi=int(mp_maxpart_factor*real(maxpart)/real(mp_partgroup_np))
     if (mp_np == 1) maxpart_mpi = maxpart
@@ -365,4 +366,14 @@
     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
+    if (lroot) write(*,*) 'Using MPI_IN_PLACE operations'
+#else
+    if (lroot) allocate(creceptor0(maxreceptor,maxspec))
+    if (lroot) write(*,*) 'Not using MPI_IN_PLACE operations'
+#endif
     goto 101
 
@@ -559,5 +570,5 @@
 ! invalid particles at the end of the arrays
 
-601 do i=num_part, 1, -1
+601 do i=numpart, 1, -1
       if (itra1(i).eq.-999999999) then
         numpart=numpart-1
@@ -598,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
 
@@ -607,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)
 
@@ -616,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
@@ -645,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)
@@ -659,5 +672,5 @@
     end do
 
-    deallocate(numparticles_mpi, idx_arr, sorted)
+    deallocate(idx_arr, sorted)
 
   end subroutine mpif_calculate_part_redist
@@ -1961,5 +1974,5 @@
     if (readclouds) then
       j=j+1
-      call MPI_Irecv(ctwc(:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
+      call MPI_Irecv(ctwc(:,:,mind),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
            &MPI_COMM_WORLD,reqs(j),mp_ierr)
       if (mp_ierr /= 0) goto 600
@@ -2326,5 +2339,5 @@
       if (readclouds) then
         j=j+1
-        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
+        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
              &MPI_COMM_WORLD,reqs(j),mp_ierr)
         if (mp_ierr /= 0) goto 600
@@ -2462,9 +2475,25 @@
     end if
 
+  ! Receptor concentrations    
+    if (lroot) then
+      call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
+           & mp_comm_used,mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    else
+      call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, &
+           & mp_comm_used,mp_ierr)
+    end if
+
 #else
 
       call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
            & mp_comm_used, mp_ierr)
+      if (mp_ierr /= 0) goto 600
       if (lroot) gridunc = gridunc0
+
+      call MPI_Reduce(creceptor, creceptor0,rcpt_size,mp_sp,MPI_SUM,id_root, &
+           & mp_comm_used,mp_ierr)
+      if (mp_ierr /= 0) goto 600
+      if (lroot) creceptor = creceptor0
 
 #endif
@@ -2482,13 +2511,4 @@
     end if
 
-! Receptor concentrations    
-    if (lroot) then
-      call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
-           & mp_comm_used,mp_ierr)
-      if (mp_ierr /= 0) goto 600
-    else
-      call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, &
-           & mp_comm_used,mp_ierr)
-    end if
 
     if (mp_measure_time) call mpif_mtime('commtime',1)
@@ -2700,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')
@@ -2788,8 +2808,8 @@
           write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR GETFIELDS:',&
                & mp_getfields_time_total
-          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
-               & mp_readwind_wtime_total
-          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
-               & mp_readwind_time_total
+!          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
+!               & mp_readwind_wtime_total
+!          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
+!               & mp_readwind_time_total
           write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FILE IO:',&
                & mp_io_wtime_total
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_average.f90
===================================================================
--- src/partoutput_average.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ src/partoutput_average.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
@@ -0,0 +1,221 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine partoutput_average(itime)
+  !                             i
+  !*****************************************************************************
+  !                                                                            *
+  !     Dump all particle positions                                            *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     12 March 1999                                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  !                                                                            *
+  !*****************************************************************************
+
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,i,j,jjjjmmdd,ihmmss
+  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
+  real :: xlon,ylat
+  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
+  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
+  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
+  real :: tr(2),tri,zlim
+  character :: adate*8,atime*6
+
+  integer(kind=2) :: ishort_xlon(maxpart),ishort_ylat(maxpart),ishort_z(maxpart)
+  integer(kind=2) :: ishort_topo(maxpart),ishort_tro(maxpart),ishort_hmix(maxpart)
+  integer(kind=2) :: ishort_pv(maxpart),ishort_rho(maxpart),ishort_qv(maxpart)
+  integer(kind=2) :: ishort_tt(maxpart),ishort_uu(maxpart),ishort_vv(maxpart)
+  integer(kind=2) :: ishort_energy(maxpart)
+
+
+  ! Determine current calendar date, needed for the file name
+  !**********************************************************
+
+  jul=bdate+real(itime,kind=dp)/86400._dp
+  call caldate(jul,jjjjmmdd,ihmmss)
+  write(adate,'(i8.8)') jjjjmmdd
+  write(atime,'(i6.6)') ihmmss
+
+
+  ! Some variables needed for temporal interpolation
+  !*************************************************
+
+  dt1=real(itime-memtime(1))
+  dt2=real(memtime(2)-itime)
+  dtt=1./(dt1+dt2)
+
+  ! Open output file and write the output
+  !**************************************
+
+    open(unitpartout_average,file=path(2)(1:length(2))//'partposit_average_'//adate// &
+         atime,form='unformatted',access='direct',status='replace',recl=24)
+
+
+  ! Write current time to file
+  !***************************
+
+!  write(unitpartout_average) itime,numpart
+  do i=1,numpart
+
+  ! Take only valid particles
+  !**************************
+
+    if (itra1(i).eq.itime) then
+      part_av_topo(i)=part_av_topo(i)/float(npart_av(i))
+      part_av_z(i)=part_av_z(i)/float(npart_av(i))
+      part_av_pv(i)=part_av_pv(i)/float(npart_av(i))
+      part_av_qv(i)=part_av_qv(i)/float(npart_av(i))
+      part_av_tt(i)=part_av_tt(i)/float(npart_av(i))
+      part_av_uu(i)=part_av_uu(i)/float(npart_av(i))
+      part_av_vv(i)=part_av_vv(i)/float(npart_av(i))
+      part_av_rho(i)=part_av_rho(i)/float(npart_av(i))
+      part_av_tro(i)=part_av_tro(i)/float(npart_av(i))
+      part_av_hmix(i)=part_av_hmix(i)/float(npart_av(i))
+      part_av_energy(i)=part_av_energy(i)/float(npart_av(i))
+
+      part_av_cartx(i)=part_av_cartx(i)/float(npart_av(i))
+      part_av_carty(i)=part_av_carty(i)/float(npart_av(i))
+      part_av_cartz(i)=part_av_cartz(i)/float(npart_av(i))
+
+! Project Cartesian coordinates back onto Earth's surface
+! #######################################################
+      xlon=atan2(part_av_cartx(i),-1.*part_av_carty(i))
+      ylat=atan2(part_av_cartz(i),sqrt(part_av_cartx(i)*part_av_cartx(i)+ &
+           part_av_carty(i)*part_av_carty(i)))
+      xlon=xlon/pi180
+      ylat=ylat/pi180
+
+
+  ! Convert all data to integer*2 variables (from -32768 to 32767) for compressed output
+  !*************************************************************************************
+
+      if (xlon.gt.180.) xlon=xlon-360.
+      if (xlon.lt.-180.) xlon=xlon+360.
+      ishort_xlon(i)=nint(xlon*180.)
+      ishort_ylat(i)=nint(ylat*360.)
+
+      zlim=(part_av_z(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_z(i)=nint(zlim)
+
+      zlim=(part_av_topo(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_topo(i)=nint(zlim)
+
+      zlim=(part_av_tro(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tro(i)=nint(zlim)
+
+      zlim=(part_av_hmix(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_hmix(i)=nint(zlim)
+
+      zlim=(part_av_rho(i)*20000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_rho(i)=nint(zlim)
+
+      zlim=(part_av_qv(i)*1000000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_qv(i)=nint(zlim)
+
+      zlim=(part_av_pv(i)*100.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_pv(i)=nint(zlim)
+
+      zlim=((part_av_tt(i)-273.15))*300.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tt(i)=nint(zlim)
+
+      zlim=(part_av_uu(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_uu(i)=nint(zlim)
+
+      zlim=(part_av_vv(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_vv(i)=nint(zlim)
+
+      zlim=(part_av_energy(i)-300000.)/30.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_energy(i)=nint(zlim)
+
+  ! Turn on for readable test output
+  !*********************************
+
+!        write(119,*) itime,i,xlon,ylat,part_av_z(i),part_av_topo(i),part_av_tro(i), &
+!        part_av_hmix(i),part_av_rho(i),part_av_qv(i),part_av_pv(i),part_av_tt(i), &
+!        ishort_uu(i),ishort_vv(i)
+    endif
+
+! Re-initialize averages, and set number of steps to zero
+    npart_av(i)=0
+    part_av_topo(i)=0.
+    part_av_z(i)=0.
+    part_av_cartx(i)=0.
+    part_av_carty(i)=0.
+    part_av_cartz(i)=0.
+    part_av_pv(i)=0.
+    part_av_qv(i)=0.
+    part_av_tt(i)=0.
+    part_av_uu(i)=0.
+    part_av_vv(i)=0.
+    part_av_rho(i)=0.
+    part_av_tro(i)=0.
+    part_av_hmix(i)=0.
+    part_av_energy(i)=0.
+    
+  end do
+
+  ! Write the output
+  !*****************
+  do i=1,numpart
+    if (itra1(i).eq.itime) then
+      write(unitpartout_average,rec=i) ishort_xlon(i),ishort_ylat(i),ishort_z(i), &
+           ishort_topo(i),ishort_tro(i),ishort_hmix(i),ishort_rho(i),ishort_qv(i),ishort_pv(i), &
+           ishort_tt(i),ishort_uu(i),ishort_vv(i)  ! ,ishort_energy(i)
+    endif
+  enddo
+
+
+  close(unitpartout_average)
+
+end subroutine partoutput_average
Index: src/partoutput_average_mpi.f90
===================================================================
--- src/partoutput_average_mpi.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ src/partoutput_average_mpi.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
@@ -0,0 +1,230 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine partoutput_average(itime,irec)
+  !                             i
+  !*****************************************************************************
+  !                                                                            *
+  !     Dump all particle positions                                            *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     12 March 1999                                                          *
+  !     03/2019 AST: Version for MPI                                           *
+  !                  processes sequentially access and append data to file     *
+  !                                                                            *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  !                                                                            *
+  !*****************************************************************************
+
+  use par_mod
+  use com_mod
+  use mpi_mod
+
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,i,j,jjjjmmdd,ihmmss,irec
+  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
+  real :: xlon,ylat
+  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
+  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
+  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
+  real :: tr(2),tri,zlim
+  character :: adate*8,atime*6
+  character(LEN=8) :: file_stat='OLD'
+
+  integer(kind=2) :: ishort_xlon(maxpart_mpi),ishort_ylat(maxpart_mpi),ishort_z(maxpart_mpi)
+  integer(kind=2) :: ishort_topo(maxpart_mpi),ishort_tro(maxpart_mpi),ishort_hmix(maxpart_mpi)
+  integer(kind=2) :: ishort_pv(maxpart_mpi),ishort_rho(maxpart_mpi),ishort_qv(maxpart_mpi)
+  integer(kind=2) :: ishort_tt(maxpart_mpi),ishort_uu(maxpart_mpi),ishort_vv(maxpart_mpi)
+  integer(kind=2) :: ishort_energy(maxpart_mpi)
+
+  ! MPI root process creates/overwrites the file, other processes append to it
+  if (lroot) file_stat='REPLACE'
+
+  ! Determine current calendar date, needed for the file name
+  !**********************************************************
+
+  jul=bdate+real(itime,kind=dp)/86400._dp
+  call caldate(jul,jjjjmmdd,ihmmss)
+  write(adate,'(i8.8)') jjjjmmdd
+  write(atime,'(i6.6)') ihmmss
+
+
+  ! Some variables needed for temporal interpolation
+  !*************************************************
+
+  dt1=real(itime-memtime(1))
+  dt2=real(memtime(2)-itime)
+  dtt=1./(dt1+dt2)
+
+  ! Open output file and write the output
+  !**************************************
+
+  open(unitpartout_average,file=path(2)(1:length(2))//'partposit_average_'//adate// &
+       atime,form='unformatted',access='direct',status=file_stat,recl=24)
+
+
+  ! Write current time to file
+  !***************************
+
+!  write(unitpartout_average) itime,numpart
+  do i=1,numpart
+
+  ! Take only valid particles
+  !**************************
+
+    if (itra1(i).eq.itime) then
+      part_av_topo(i)=part_av_topo(i)/float(npart_av(i))
+      part_av_z(i)=part_av_z(i)/float(npart_av(i))
+      part_av_pv(i)=part_av_pv(i)/float(npart_av(i))
+      part_av_qv(i)=part_av_qv(i)/float(npart_av(i))
+      part_av_tt(i)=part_av_tt(i)/float(npart_av(i))
+      part_av_uu(i)=part_av_uu(i)/float(npart_av(i))
+      part_av_vv(i)=part_av_vv(i)/float(npart_av(i))
+      part_av_rho(i)=part_av_rho(i)/float(npart_av(i))
+      part_av_tro(i)=part_av_tro(i)/float(npart_av(i))
+      part_av_hmix(i)=part_av_hmix(i)/float(npart_av(i))
+      part_av_energy(i)=part_av_energy(i)/float(npart_av(i))
+
+      part_av_cartx(i)=part_av_cartx(i)/float(npart_av(i))
+      part_av_carty(i)=part_av_carty(i)/float(npart_av(i))
+      part_av_cartz(i)=part_av_cartz(i)/float(npart_av(i))
+
+! Project Cartesian coordinates back onto Earth's surface
+! #######################################################
+      xlon=atan2(part_av_cartx(i),-1.*part_av_carty(i))
+      ylat=atan2(part_av_cartz(i),sqrt(part_av_cartx(i)*part_av_cartx(i)+ &
+           part_av_carty(i)*part_av_carty(i)))
+      xlon=xlon/pi180
+      ylat=ylat/pi180
+
+
+  ! Convert all data to integer*2 variables (from -32768 to 32767) for compressed output
+  !*************************************************************************************
+
+      if (xlon.gt.180.) xlon=xlon-360.
+      if (xlon.lt.-180.) xlon=xlon+360.
+      ishort_xlon(i)=nint(xlon*180.)
+      ishort_ylat(i)=nint(ylat*360.)
+
+      zlim=(part_av_z(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_z(i)=nint(zlim)
+
+      zlim=(part_av_topo(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_topo(i)=nint(zlim)
+
+      zlim=(part_av_tro(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tro(i)=nint(zlim)
+
+      zlim=(part_av_hmix(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_hmix(i)=nint(zlim)
+
+      zlim=(part_av_rho(i)*20000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_rho(i)=nint(zlim)
+
+      zlim=(part_av_qv(i)*1000000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_qv(i)=nint(zlim)
+
+      zlim=(part_av_pv(i)*100.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_pv(i)=nint(zlim)
+
+      zlim=((part_av_tt(i)-273.15))*300.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tt(i)=nint(zlim)
+
+      zlim=(part_av_uu(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_uu(i)=nint(zlim)
+
+      zlim=(part_av_vv(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_vv(i)=nint(zlim)
+
+      zlim=(part_av_energy(i)-300000.)/30.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_energy(i)=nint(zlim)
+
+
+  ! Turn on for readable test output
+  !*********************************
+
+!        write(119,*) itime,i,xlon,ylat,part_av_z(i),part_av_topo(i),part_av_tro(i), &
+!        part_av_hmix(i),part_av_rho(i),part_av_qv(i),part_av_pv(i),part_av_tt(i), &
+!        ishort_uu(i),ishort_vv(i)
+    endif
+
+! Re-initialize averages, and set number of steps to zero
+    npart_av(i)=0
+    part_av_topo(i)=0.
+    part_av_z(i)=0.
+    part_av_cartx(i)=0.
+    part_av_carty(i)=0.
+    part_av_cartz(i)=0.
+    part_av_pv(i)=0.
+    part_av_qv(i)=0.
+    part_av_tt(i)=0.
+    part_av_uu(i)=0.
+    part_av_vv(i)=0.
+    part_av_rho(i)=0.
+    part_av_tro(i)=0.
+    part_av_hmix(i)=0.
+    part_av_energy(i)=0.
+    
+  end do
+
+  ! Write the output
+  !*****************
+
+  do i=1,numpart
+    if (itra1(i).eq.itime) then
+      write(unitpartout_average,rec=irec+i) ishort_xlon(i),ishort_ylat(i),ishort_z(i), &
+           ishort_topo(i),ishort_tro(i),ishort_hmix(i),ishort_rho(i),ishort_qv(i),ishort_pv(i), &
+           ishort_tt(i),ishort_uu(i),ishort_vv(i)  ! ,ishort_energy(i)
+    endif
+  enddo
+
+  close(unitpartout_average)
+
+end subroutine partoutput_average
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/partpos_average.f90
===================================================================
--- src/partpos_average.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ src/partpos_average.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
@@ -0,0 +1,206 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+
+subroutine partpos_average(itime,j)
+
+
+!**********************************************************************
+! This subroutine averages particle quantities, to be used for particle
+! dump (in partoutput.f90). Averaging is done over output interval.
+!**********************************************************************
+
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  integer :: itime,j,ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
+  real :: xlon,ylat,x,y,z
+  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
+  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
+  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
+  real :: uu1(2),uuprof(2),uui,vv1(2),vvprof(2),vvi
+  real :: tr(2),tri,energy
+
+
+
+ ! Some variables needed for temporal interpolation
+  !*************************************************
+
+  dt1=real(itime-memtime(1))
+  dt2=real(memtime(2)-itime)
+  dtt=1./(dt1+dt2)
+
+  xlon=xlon0+xtra1(j)*dx
+  ylat=ylat0+ytra1(j)*dy
+
+  !*****************************************************************************
+  ! Interpolate several variables (PV, specific humidity, etc.) to particle position
+  !*****************************************************************************
+
+  ix=xtra1(j)
+  jy=ytra1(j)
+  ixp=ix+1
+  jyp=jy+1
+  ddx=xtra1(j)-real(ix)
+  ddy=ytra1(j)-real(jy)
+  rddx=1.-ddx
+  rddy=1.-ddy
+  p1=rddx*rddy
+  p2=ddx*rddy
+  p3=rddx*ddy
+  p4=ddx*ddy
+
+! eso: Temporary fix for particle exactly at north pole
+  if (jyp >= nymax) then
+  !  write(*,*) 'WARNING: conccalc.f90 jyp >= nymax'
+    jyp=jyp-1
+  end if
+
+  ! Topography
+  !***********
+
+  topo=p1*oro(ix,jy)+p2*oro(ixp,jy)+p3*oro(ix,jyp)+p4*oro(ixp,jyp)
+
+ ! Potential vorticity, specific humidity, temperature, and density
+  !*****************************************************************
+
+  do il=2,nz
+    if (height(il).gt.ztra1(j)) then
+      indz=il-1
+      indzp=il
+      goto 6
+    endif
+  end do
+6 continue
+
+  dz1=ztra1(j)-height(indz)
+  dz2=height(indzp)-ztra1(j)
+  dz=1./(dz1+dz2)
+
+
+  do ind=indz,indzp
+    do m=1,2
+      indexh=memind(m)
+
+  ! Potential vorticity
+      pv1(m)=p1*pv(ix ,jy ,ind,indexh) &
+            +p2*pv(ixp,jy ,ind,indexh) &
+            +p3*pv(ix ,jyp,ind,indexh) &
+            +p4*pv(ixp,jyp,ind,indexh)
+  ! Specific humidity
+      qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
+            +p2*qv(ixp,jy ,ind,indexh) &
+            +p3*qv(ix ,jyp,ind,indexh) &
+            +p4*qv(ixp,jyp,ind,indexh)
+  ! Temperature
+      tt1(m)=p1*tt(ix ,jy ,ind,indexh) &
+            +p2*tt(ixp,jy ,ind,indexh) &
+            +p3*tt(ix ,jyp,ind,indexh) &
+            +p4*tt(ixp,jyp,ind,indexh)
+  ! U wind
+      uu1(m)=p1*uu(ix ,jy ,ind,indexh) &
+            +p2*uu(ixp,jy ,ind,indexh) &
+            +p3*uu(ix ,jyp,ind,indexh) &
+            +p4*uu(ixp,jyp,ind,indexh)
+  ! V wind
+      vv1(m)=p1*vv(ix ,jy ,ind,indexh) &
+            +p2*vv(ixp,jy ,ind,indexh) &
+            +p3*vv(ix ,jyp,ind,indexh) &
+            +p4*vv(ixp,jyp,ind,indexh)
+  ! Density
+      rho1(m)=p1*rho(ix ,jy ,ind,indexh) &
+             +p2*rho(ixp,jy ,ind,indexh) &
+             +p3*rho(ix ,jyp,ind,indexh) &
+             +p4*rho(ixp,jyp,ind,indexh)
+    end do
+    pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt
+    qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
+    ttprof(ind-indz+1)=(tt1(1)*dt2+tt1(2)*dt1)*dtt
+    uuprof(ind-indz+1)=(uu1(1)*dt2+uu1(2)*dt1)*dtt
+    vvprof(ind-indz+1)=(vv1(1)*dt2+vv1(2)*dt1)*dtt
+    rhoprof(ind-indz+1)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
+  end do
+  pvi=(dz1*pvprof(2)+dz2*pvprof(1))*dz
+  qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
+  tti=(dz1*ttprof(2)+dz2*ttprof(1))*dz
+  uui=(dz1*uuprof(2)+dz2*uuprof(1))*dz
+  vvi=(dz1*vvprof(2)+dz2*vvprof(1))*dz
+  rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
+
+  ! Tropopause and PBL height
+  !**************************
+
+  do m=1,2
+    indexh=memind(m)
+
+  ! Tropopause
+    tr(m)=p1*tropopause(ix ,jy ,1,indexh) &
+        + p2*tropopause(ixp,jy ,1,indexh) &
+        + p3*tropopause(ix ,jyp,1,indexh) &
+        + p4*tropopause(ixp,jyp,1,indexh)
+
+  ! PBL height
+    hm(m)=p1*hmix(ix ,jy ,1,indexh) &
+        + p2*hmix(ixp,jy ,1,indexh) &
+        + p3*hmix(ix ,jyp,1,indexh) &
+        + p4*hmix(ixp,jyp,1,indexh)
+  end do
+
+  hmixi=(hm(1)*dt2+hm(2)*dt1)*dtt
+  tri=(tr(1)*dt2+tr(2)*dt1)*dtt
+
+
+  energy=tti*cpa+(ztra1(j)+topo)*9.81+qvi*2501000.+(uui**2+vvi**2)/2.
+
+  ! Add new values to sum and increase counter by one
+  !**************************************************
+
+  npart_av(j)=npart_av(j)+1
+
+  ! Calculate Cartesian 3D coordinates suitable for averaging
+  !**********************************************************
+
+  xlon=xlon*pi180
+  ylat=ylat*pi180
+  x = cos(ylat)*sin(xlon)
+  y = -1.*cos(ylat)*cos(xlon)
+  z = sin(ylat)
+
+  part_av_cartx(j)=part_av_cartx(j)+x
+  part_av_carty(j)=part_av_carty(j)+y
+  part_av_cartz(j)=part_av_cartz(j)+z
+  part_av_z(j)=part_av_z(j)+ztra1(j)
+  part_av_topo(j)=part_av_topo(j)+topo
+  part_av_pv(j)=part_av_pv(j)+pvi
+  part_av_qv(j)=part_av_qv(j)+qvi
+  part_av_tt(j)=part_av_tt(j)+tti
+  part_av_uu(j)=part_av_uu(j)+uui
+  part_av_vv(j)=part_av_vv(j)+vvi
+  part_av_rho(j)=part_av_rho(j)+rhoi
+  part_av_tro(j)=part_av_tro(j)+tri
+  part_av_hmix(j)=part_av_hmix(j)+hmixi
+  part_av_energy(j)=part_av_energy(j)+energy
+
+
+return
+end subroutine partpos_average 
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/readwind_gfs.f90
===================================================================
--- src/readwind_gfs.f90	(revision db91eb78afeb1bb8bafe66068877db8841859062)
+++ src/readwind_gfs.f90	(revision 5d74ed914bb3649af12dc589a95068bd49aa0fdb)
@@ -83,5 +83,5 @@
 
   ! NCEP
-  integer :: numpt,numpu,numpv,numpw,numprh
+  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
   real :: help, temp, ew
   real :: elev
@@ -134,4 +134,5 @@
   numpw=0
   numprh=0
+  numpclwch=0
   ifield=0
 10   ifield=ifield+1
@@ -557,8 +558,20 @@
       endif
 ! SEC & IP 12/2018 read GFS clouds
-      if(isec1(6).eq.153) then  !! CLWCR  Cloud liquid water content [kg/kg]
-        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+      if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg] 
+         if((i.eq.0).and.(j.eq.0)) then
+            do ii=1,nuvz
+              if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
+            end do
+        endif
+        help=zsec4(nxfield*(ny-j-1)+i+1)
+        if(i.le.i180) then
+          clwch(i179+i,j,numpclwch,n)=help
+        else
+          clwch(i-i181,j,numpclwch,n)=help
+        endif
         readclouds=.true.
         sumclouds=.true.
+!        readclouds=.false.
+!       sumclouds=.false.
       endif
 
Index: src/releaseparticles.f90
===================================================================
--- src/releaseparticles.f90	(revision 75a4dedbe36ebfba4e9c67fce8af8641dafcd8eb)
+++ src/releaseparticles.f90	(revision 7873bf796f5ee67aab272b1243107198ceb275b1)
@@ -114,5 +114,6 @@
       average_timecorrect=0.
       do k=1,nspec
-        if (zpoint1(i).gt.0.5) then      ! point source
+        if(abs(xpoint2(i)-xpoint1(i)).lt.1.E-4.and.abs(ypoint2(i)-ypoint1(i)).lt.1.E-4) then
+!        if (zpoint1(i).gt.0.5) then      ! point source
           timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek)
         else                             ! area source
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 20963b15051b22eeb381dcf4e824a9104d92717b)
+++ src/timemanager_mpi.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
@@ -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
@@ -532,5 +533,5 @@
                 end if
 
-              else  ! :TODO: check for zeroing in the netcdf module
+              else
                 call concoutput_surf_nest(itime,outnum)
               end if
@@ -593,8 +594,15 @@
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
         if (ipout.ge.1) then
+          if (mp_measure_time) call mpif_mtime('iotime',0)
+          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
+          if (mp_measure_time) call mpif_mtime('iotime',1)
         end if
 
@@ -757,4 +765,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 +908,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
@@ -654,5 +653,5 @@
         if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
 
-          do kz=nz,1,-1 !go Bottom up!
+          do kz=nz,2,-1 !go Bottom up!
             if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 
Index: src/verttransform_gfs.f90
===================================================================
--- src/verttransform_gfs.f90	(revision db91eb78afeb1bb8bafe66068877db8841859062)
+++ src/verttransform_gfs.f90	(revision 437c545fb9b9ece116fae2e06166d88b979f91a7)
@@ -548,10 +548,10 @@
         if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
 
-          do kz=nz,1,-1 !go Bottom up!
+          do kz=nz,2,-1 !go Bottom up!
             if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
               clouds(ix,jy,kz,n)=1                               ! is a cloud
               if (lsp.ge.convp) then
-                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
+                clouds(ix,jy,kz,n)=3                             ! lsp in-cloud
               else
                 clouds(ix,jy,kz,n)=2                             ! convp in-cloud
