Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision 6b22af900a57ac3fc1a984db76f4ef2cb0436af0)
+++ src/concoutput.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -63,4 +63,5 @@
   use par_mod
   use com_mod
+  use mean_mod
 
   implicit none
@@ -90,7 +91,8 @@
 
 !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
-  real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
-  real :: drygridtotal,drygridsigmatotal,drygridtotalunc
+  real(dep_prec) :: auxgrid(nclassunc)
+  real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
+  real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
+  real(dep_prec) :: drygridtotal,drygridsigmatotal,drygridtotalunc
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/concoutput_mpi.f90
===================================================================
--- src/concoutput_mpi.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/concoutput_mpi.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -68,4 +68,5 @@
   use com_mod
   use mpi_mod
+  use mean_mod
 
   implicit none
@@ -95,7 +96,8 @@
 
 !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
-  real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
-  real :: drygridtotal,drygridsigmatotal,drygridtotalunc
+  real(dep_prec) :: auxgrid(nclassunc)
+  real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
+  real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
+  real(dep_prec) :: drygridtotal,drygridsigmatotal,drygridtotalunc
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/concoutput_nest.f90
===================================================================
--- src/concoutput_nest.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/concoutput_nest.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -61,4 +61,5 @@
   use par_mod
   use com_mod
+  use mean_mod
 
   implicit none
@@ -88,5 +89,5 @@
 
   !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc)
+  real(dep_prec) :: auxgrid(nclassunc)
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/concoutput_nest_mpi.f90
===================================================================
--- src/concoutput_nest_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/concoutput_nest_mpi.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -65,4 +65,5 @@
   use com_mod
   use mpi_mod
+  use mean_mod
 
   implicit none
@@ -92,5 +93,5 @@
 
   !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc)
+  real(dep_prec) :: auxgrid(nclassunc)
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/concoutput_surf.f90
===================================================================
--- src/concoutput_surf.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/concoutput_surf.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -63,4 +63,5 @@
   use par_mod
   use com_mod
+  use mean_mod
 
   implicit none
@@ -90,7 +91,8 @@
 
   !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
-  real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
-  real :: drygridtotal,drygridsigmatotal,drygridtotalunc
+  real(dep_prec) :: auxgrid(nclassunc)
+  real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
+  real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
+  real(dep_prec) :: drygridtotal,drygridsigmatotal,drygridtotalunc
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/concoutput_surf_mpi.f90
===================================================================
--- src/concoutput_surf_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/concoutput_surf_mpi.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -69,4 +69,5 @@
   use com_mod
   use mpi_mod
+  use mean_mod
 
   implicit none
@@ -96,7 +97,8 @@
 
   !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
-  real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
-  real :: drygridtotal,drygridsigmatotal,drygridtotalunc
+  real(dep_prec) :: auxgrid(nclassunc)
+  real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
+  real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
+  real(dep_prec) :: drygridtotal,drygridsigmatotal,drygridtotalunc
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/concoutput_surf_nest.f90
===================================================================
--- src/concoutput_surf_nest.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/concoutput_surf_nest.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -61,4 +61,5 @@
   use par_mod
   use com_mod
+  use mean_mod
 
   implicit none
@@ -88,5 +89,5 @@
 
   !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc)
+  real(dep_prec) :: auxgrid(nclassunc)
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/concoutput_surf_nest_mpi.f90
===================================================================
--- src/concoutput_surf_nest_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/concoutput_surf_nest_mpi.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -64,4 +64,5 @@
   use com_mod
   use mpi_mod
+  use mean_mod
 
   implicit none
@@ -91,5 +92,5 @@
 
   !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  real :: auxgrid(nclassunc)
+  real(dep_prec) :: auxgrid(nclassunc)
   real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Index: src/makefile
===================================================================
--- src/makefile	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/makefile	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -82,5 +82,5 @@
 xmass_mod.o             flux_mod.o \
 point_mod.o             outg_mod.o \
-	 		random_mod.o
+mean_mod.o	 	random_mod.o
 
 MPI_MODOBJS = \
@@ -167,5 +167,5 @@
 interpol_vdep.o         interpol_rain.o \
 hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
+                  wetdepo.o \
 hanna_short.o           windalign.o \
 hanna1.o                	\
@@ -293,16 +293,16 @@
 conccalc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
 conccalc_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o unc_mod.o
-concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
+concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
 concoutput_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
-	unc_mod.o
-concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
+	unc_mod.o mean_mod.o
+concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
 concoutput_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
-	unc_mod.o
-concoutput_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
+	unc_mod.o mean_mod.o
+concoutput_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
 concoutput_surf_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
-	unc_mod.o
-concoutput_surf_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
+	unc_mod.o mean_mod.o
+concoutput_surf_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
 concoutput_surf_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o \
-	point_mod.o unc_mod.o
+	point_mod.o unc_mod.o mean_mod.o
 conv_mod.o: par_mod.o
 convect43c.o: conv_mod.o par_mod.o
@@ -359,6 +359,7 @@
 interpol_wind_short_nests.o: com_mod.o interpol_mod.o par_mod.o
 juldate.o: par_mod.o
+mean_mod.o: par_mod.o
 mpi_mod.o: com_mod.o par_mod.o unc_mod.o
-netcdf_output_mod.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
+netcdf_output_mod.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
 obukhov.o: par_mod.o
 obukhov_gfs.o: par_mod.o
@@ -366,4 +367,5 @@
 openouttraj.o: com_mod.o par_mod.o point_mod.o
 openreceptors.o: com_mod.o par_mod.o
+outg_mod.o: par_mod.o
 outgrid_init.o: com_mod.o flux_mod.o oh_mod.o outg_mod.o par_mod.o unc_mod.o
 outgrid_init_nest.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
@@ -376,5 +378,5 @@
 partoutput_short_mpi.o: com_mod.o mpi_mod.o par_mod.o
 pbl_profile.o: par_mod.o
-plumetraj.o: com_mod.o par_mod.o point_mod.o
+plumetraj.o: com_mod.o par_mod.o point_mod.o mean_mod.o
 psih.o: par_mod.o
 psim.o: par_mod.o
Index: c/mean.f90
===================================================================
--- src/mean.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,66 +1,0 @@
-!**********************************************************************
-! 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 mean(x,xm,xs,number)
-
-  !*****************************************************************************
-  !                                                                            *
-  !  This subroutine calculates mean and standard deviation of a given element.*
-  !                                                                            *
-  !      AUTHOR: Andreas Stohl, 25 January 1994                                *
-  !                                                                            *
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  ! x(number)           field of input data                                    *
-  ! xm                  mean                                                   *
-  ! xs                  standard deviation                                     *
-  ! number              number of elements of field x                          *
-  !                                                                            *
-  ! Constants:                                                                 *
-  ! eps                 tiny number                                            *
-  !                                                                            *
-  !*****************************************************************************
-
-  implicit none
-
-  integer :: number,i
-  real :: x(number),xm,xs,xl,xq,xaux
-  real,parameter :: eps=1.0e-30
-
-  xl=0.
-  xq=0.
-  do i=1,number
-    xl=xl+x(i)
-    xq=xq+x(i)*x(i)
-  end do
-
-  xm=xl/real(number)
-
-  xaux=xq-xl*xl/real(number)
-
-  if (xaux.lt.eps) then
-    xs=0.
-  else
-    xs=sqrt(xaux/real(number-1))
-  endif
-
-end subroutine mean
Index: src/mean_mod.f90
===================================================================
--- src/mean_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/mean_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -0,0 +1,195 @@
+!**********************************************************************
+! 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/>.   *
+!**********************************************************************
+
+module mean_mod
+  public
+
+! Interface to select single or double precision version of the 'mean'
+! function from type of input arguments ("function overloading")
+!********************************************************************
+  interface mean
+    module procedure mean_sp
+    module procedure mean_dp
+    module procedure mean_mixed_prec
+  end interface mean
+
+contains
+
+  subroutine mean_sp(x_sp,xm,xs,number)
+
+!*****************************************************************************
+!                                                                            *
+!  This subroutine calculates mean and standard deviation of a given element.*
+!                                                                            *
+!      AUTHOR: Andreas Stohl, 25 January 1994                                *
+!                                                                            *
+!      Single precision version ESO 2016                                     *
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+! x_sp(number)        field of input data                                    *
+! xm                  mean                                                   *
+! xs                  standard deviation                                     *
+! number              number of elements of field x_sp                       *
+!                                                                            *
+! Constants:                                                                 *
+! eps                 tiny number                                            *
+!                                                                            *
+!*****************************************************************************
+
+    use par_mod, only: sp
+
+    implicit none
+
+    ! integer :: number,i
+    ! real(sp) :: x_sp(number),xm,xs,xl,xq,xaux
+    ! real(sp),parameter :: eps=1.0e-30
+
+    real(sp), intent(in) :: x_sp(number)
+    real(sp), intent(out) ::xm,xs
+    integer,intent(in) :: number
+    real(sp) :: xl,xq,xaux
+    real(sp),parameter :: eps=1.0e-30
+    integer :: i
+
+    xl=0.
+    xq=0.
+    do i=1,number
+      xl=xl+x_sp(i)
+      xq=xq+x_sp(i)*x_sp(i)
+    end do
+
+    xm=xl/real(number,kind=sp)
+
+    xaux=xq-xl*xl/real(number,kind=sp)
+
+    if (xaux.lt.eps) then
+      xs=0.
+    else
+      xs=sqrt(xaux/real(number-1,kind=sp))
+    endif
+
+  end subroutine mean_sp
+
+  subroutine mean_dp(x_dp,xm,xs,number)
+
+!*****************************************************************************
+!                                                                            *
+!  This subroutine calculates mean and standard deviation of a given element.*
+!                                                                            *
+!      AUTHOR: Andreas Stohl, 25 January 1994                                *
+!                                                                            *
+!      Double precision version ESO 2016                                     *
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+! x_dp(number)        field of input data                                    *
+! xm                  mean                                                   *
+! xs                  standard deviation                                     *
+! number              number of elements of field x_dp                       *
+!                                                                            *
+! Constants:                                                                 *
+! eps                 tiny number                                            *
+!                                                                            *
+!*****************************************************************************
+
+    use par_mod, only: dp
+
+    implicit none
+
+    real(dp), intent(in) :: x_dp(number)
+    real(dp), intent(out) ::xm,xs
+    integer,intent(in) :: number
+    real(dp) :: xl,xq,xaux
+    real(dp),parameter :: eps=1.0e-30
+    integer :: i
+
+    xl=0._dp
+    xq=0._dp
+    do i=1,number
+      xl=xl+x_dp(i)
+      xq=xq+x_dp(i)*x_dp(i)
+    end do
+
+    xm=xl/real(number,kind=dp)
+
+    xaux=xq-xl*xl/real(number,kind=dp)
+
+    if (xaux.lt.eps) then
+      xs=0._dp
+    else
+      xs=sqrt(xaux/real(number-1,kind=dp))
+    endif
+
+  end subroutine mean_dp
+
+  subroutine mean_mixed_prec(x_dp,xm,xs,number)
+
+!*****************************************************************************
+!                                                                            *
+!  This subroutine calculates mean and standard deviation of a given element.*
+!                                                                            *
+!      AUTHOR: Andreas Stohl, 25 January 1994                                *
+!                                                                            *
+!      Mixed precision version ESO 2016 (dp input, sp output)                *
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+! x_dp(number)        field of input data                                    *
+! xm                  mean                                                   *
+! xs                  standard deviation                                     *
+! number              number of elements of field x_dp                       *
+!                                                                            *
+! Constants:                                                                 *
+! eps                 tiny number                                            *
+!                                                                            *
+!*****************************************************************************
+
+    use par_mod, only: sp,dp
+
+    implicit none
+
+    real(dp), intent(in) :: x_dp(number)
+    real(sp), intent(out) ::xm,xs
+    integer,intent(in) :: number
+    real(sp) :: xl,xq,xaux
+    real(sp),parameter :: eps=1.0e-30
+    integer :: i
+
+    xl=0._sp
+    xq=0._sp
+    do i=1,number
+      xl=xl+x_dp(i)
+      xq=xq+x_dp(i)*x_dp(i)
+    end do
+
+    xm=xl/real(number,kind=sp)
+
+    xaux=xq-xl*xl/real(number,kind=sp)
+
+    if (xaux.lt.eps) then
+      xs=0._sp
+    else
+      xs=sqrt(xaux/real(number-1,kind=sp))
+    endif
+
+  end subroutine mean_mixed_prec
+end module mean_mod
Index: src/mpi_mod.f90
===================================================================
--- src/mpi_mod.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/mpi_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -53,4 +53,5 @@
 !                         unless a dedicated process runs getfields/readwind *
 ! lmp_sync                If .false., use asynchronous MPI                   *
+! mp_cp                   Real precision to use for deposition fields        *
 !                                                                            *
 !                                                                            *
@@ -84,5 +85,5 @@
   integer :: mp_seed=0
   integer, parameter :: mp_sp=MPI_REAL4, mp_dp=MPI_REAL8
-  integer, parameter :: mp_pp=mp_sp
+  integer :: mp_cp
   integer, parameter :: id_root=0 ! master process
 
@@ -163,5 +164,5 @@
 !   mp_np       number of running processes, decided at run-time
 !***********************************************************************
-    use par_mod, only: maxpart, numwfmem
+    use par_mod, only: maxpart, numwfmem, dep_prec
     use com_mod, only: mpi_mode
 
@@ -170,5 +171,5 @@
     integer :: i,j,s,addmaxpart=0
 
-! each process gets an ID (mp_pid) in the range 0,..,mp_np-1
+! Each process gets an ID (mp_pid) in the range 0,..,mp_np-1
     call MPI_INIT(mp_ierr)
     if (mp_ierr /= 0) goto 100
@@ -179,5 +180,5 @@
 
 
-! this variable is used to handle subroutines common to parallel/serial version
+! Variable mpi_mode is used to handle subroutines common to parallel/serial version
     if (lmp_sync) then
       mpi_mode=2 ! hold 2 windfields in memory
@@ -188,4 +189,17 @@
     if (mp_pid.ne.0) then
       lroot = .false.
+    end if
+
+! Set MPI precision to use for transferring deposition fields
+!************************************************************
+    if (dep_prec==dp) then
+      mp_cp = MPI_REAL8
+      if (lroot) write(*,*) 'Using double precision for deposition fields'
+    else if (dep_prec==sp) then
+      mp_cp = MPI_REAL4
+      if (lroot) write(*,*) 'Using single precision for deposition fields'
+    else
+      write(*,*) 'ERROR: something went wrong setting MPI real precision'
+      stop
     end if
 
@@ -1671,5 +1685,5 @@
 
     if ((WETDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_sp, MPI_SUM, id_root, &
+      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
            & mp_comm_used, mp_ierr)
       if (mp_ierr /= 0) goto 600
@@ -1677,5 +1691,5 @@
 
     if ((DRYDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_sp, MPI_SUM, id_root, &
+      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
            & mp_comm_used, mp_ierr)
       if (mp_ierr /= 0) goto 600
@@ -1747,5 +1761,5 @@
 
     if ((WETDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(wetgriduncn, wetgriduncn0, grid_size2d, mp_sp, MPI_SUM, id_root, &
+      call MPI_Reduce(wetgriduncn, wetgriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
            & mp_comm_used, mp_ierr)
       if (mp_ierr /= 0) goto 600
@@ -1753,5 +1767,5 @@
 
     if ((DRYDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(drygriduncn, drygriduncn0, grid_size2d, mp_sp, MPI_SUM, id_root, &
+      call MPI_Reduce(drygriduncn, drygriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
            & mp_comm_used, mp_ierr)
       if (mp_ierr /= 0) goto 600
Index: src/netcdf_output_mod.f90
===================================================================
--- src/netcdf_output_mod.f90	(revision d8107c2d5e0d3760432a003d5f2a37585183732d)
+++ src/netcdf_output_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -36,4 +36,8 @@
   !    switch write_releases
   !  - additional updates for FLEXPART 9.x
+  ! 
+  ! ESO 2016
+  !  - Deposition fields can be calculated in double precision, see variable
+  !    'dep_prec' in par_mod
   !*****************************************************************************
 
@@ -48,5 +52,5 @@
                        wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
                        area,arean,volumen, orooutn
-  use par_mod,   only: dp, maxspec, maxreceptor, nclassunc,&
+  use par_mod,   only: dep_prec, sp, dp, maxspec, maxreceptor, nclassunc,&
                        unitoutrecept,unitoutreceptppt, nxmax,unittmp
   use com_mod,   only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
@@ -70,4 +74,6 @@
                        nested_output, ipout, surf_only, linit_cond, &
                        flexversion,mpi_mode
+
+  use mean_mod
 
   implicit none
@@ -723,14 +729,19 @@
   integer, intent(in) :: itime
   real, intent(in)    :: outnum
-  real, intent(out)   :: gridtotalunc,wetgridtotalunc,drygridtotalunc
+  real(dep_prec),intent(out):: wetgridtotalunc,drygridtotalunc
+  real, intent(out)   :: gridtotalunc
   real                :: densityoutrecept(maxreceptor)
   integer             :: ncid,kp,ks,kz,ix,jy,iix,jjy,kzz,kzzm1,ngrid
-  integer             :: nage,i,l, jj
+  integer             :: nage,i,l,jj
   real                :: tot_mu(maxspec,maxpointspec_act)
   real                :: halfheight,dz,dz1,dz2
   real                :: xl,yl,xlrot,ylrot,zagnd,zagndprev
-  real                :: auxgrid(nclassunc),gridtotal,gridsigmatotal
-  real                :: wetgridtotal,wetgridsigmatotal
-  real                :: drygridtotal,drygridsigmatotal
+  real(dep_prec)      :: auxgrid(nclassunc)
+  real(dep_prec)      :: gridtotal,gridsigmatotal
+  real(dep_prec)      :: wetgridtotal,wetgridsigmatotal
+  real(dep_prec)      :: drygridtotal,drygridsigmatotal
+  ! real(sp)            :: gridtotal,gridsigmatotal
+  ! real(sp)            :: wetgridtotal,wetgridsigmatotal
+  ! real(sp)            :: drygridtotal,drygridsigmatotal
 
   real, parameter     :: weightair=28.97
@@ -889,10 +900,10 @@
                    wetgridsigma(ix,jy),nclassunc)
               ! Multiply by number of classes to get total concentration
-              wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc)
+              wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc,kind=dep_prec)
               wetgridtotal=wetgridtotal+wetgrid(ix,jy)
               ! Calculate standard deviation of the mean
               wetgridsigma(ix,jy)= &
                    wetgridsigma(ix,jy)* &
-                   sqrt(real(nclassunc))
+                   sqrt(real(nclassunc,kind=dep_prec))
               wetgridsigmatotal=wetgridsigmatotal+ &
                    wetgridsigma(ix,jy)
@@ -1061,5 +1072,6 @@
   integer, intent(in) :: itime
   real, intent(in)    :: outnum
-  real, intent(out)   :: gridtotalunc,wetgridtotalunc,drygridtotalunc
+  real(sp), intent(out)   :: gridtotalunc
+  real(dep_prec), intent(out)   :: wetgridtotalunc,drygridtotalunc
 
   print*,'Netcdf output for surface only not yet implemented'
@@ -1117,5 +1129,6 @@
   real                :: halfheight,dz,dz1,dz2
   real                :: xl,yl,xlrot,ylrot,zagnd,zagndprev
-  real                :: auxgrid(nclassunc),gridtotal
+  real(dep_prec)      :: auxgrid(nclassunc)
+  real                :: gridtotal
   real, parameter     :: weightair=28.97
 
Index: src/outg_mod.f90
===================================================================
--- src/outg_mod.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/outg_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -22,4 +22,6 @@
 module outg_mod
 
+  use par_mod, only: dep_prec
+
   implicit none
 
@@ -37,9 +39,9 @@
   real,allocatable, dimension (:,:,:) :: factor3d
   real,allocatable, dimension (:,:,:) :: grid
-  real,allocatable, dimension (:,:) :: wetgrid
-  real,allocatable, dimension (:,:) :: drygrid
+  real(dep_prec),allocatable, dimension (:,:) :: wetgrid
+  real(dep_prec),allocatable, dimension (:,:) :: drygrid
   real,allocatable, dimension (:,:,:) :: gridsigma
-  real,allocatable, dimension (:,:) :: drygridsigma
-  real,allocatable, dimension (:,:) :: wetgridsigma
+  real(dep_prec),allocatable, dimension (:,:) :: drygridsigma
+  real(dep_prec),allocatable, dimension (:,:) :: wetgridsigma
   real,allocatable, dimension (:) :: sparse_dump_r
   real,allocatable, dimension (:) :: sparse_dump_u
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/par_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -41,13 +41,16 @@
 
   !****************************************************************
-  ! Parameter defining KIND parameter for double/single precision
-  ! 'cp' sets the precision for some internal calculations (sp or dp),
-  ! sp is default but dp can be used for increased precision.
-  ! WARNING: for now, do not use 'dp'
+  ! Parameters defining KIND parameter for double/single precision
   !****************************************************************
 
   integer,parameter :: dp=selected_real_kind(P=15)
   integer,parameter :: sp=selected_real_kind(6)
-  integer,parameter :: cp=sp ! sp is default
+
+  !****************************************************************
+  ! dep_prec sets the precision for deposition calculations (sp or 
+  ! dp). sp is default, dp can be used for increased precision.
+  !****************************************************************
+
+  integer,parameter :: dep_prec=dp
 
   !***********************************************************
@@ -79,5 +82,5 @@
   real,parameter :: hmixmin=100., hmixmax=4500., turbmesoscale=0.16
   real,parameter :: d_trop=50., d_strat=0.1
-  real,parameter :: rho_water=1000 !ZHG 2015 [kg/m3]
+  real,parameter :: rho_water=1000. !ZHG 2015 [kg/m3]
   ! karman                  Karman's constant
   ! href [m]                Reference height for dry deposition
Index: src/plumetraj.f90
===================================================================
--- src/plumetraj.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/plumetraj.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -55,4 +55,5 @@
   use par_mod
   use com_mod
+  use mean_mod
 
   implicit none
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/timemanager.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -111,6 +111,7 @@
   ! real :: us(maxpart),vs(maxpart),ws(maxpart)
   ! integer(kind=2) :: cbt(maxpart)
-  real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
-  real :: drygridtotalunc,xold,yold,zold,xmassfract
+  real(sp) :: gridtotalunc
+  real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc
+  real :: xold,yold,zold,xmassfract
   real, parameter :: e_inv = 1.0/exp(1.0)
   !double precision xm(maxspec,maxpointspec_act),
@@ -595,5 +596,5 @@
 !ZHG 2015
                   !CGZ-lifetime: Check mass fraction left/save lifetime
-                   ! if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.inv_e.and.checklifetime(j,ks).eq.0.)then
+                   ! if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.e_inv.and.checklifetime(j,ks).eq.0.)then
                        !Mass below 1% of initial >register lifetime
                        ! checklifetime(j,ks)=abs(itra1(j)-itramem(j))
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/timemanager_mpi.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -114,6 +114,7 @@
   real :: decfact
 
-  real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
-  real :: drygridtotalunc,xold,yold,zold,xmassfract
+  real(sp) :: gridtotalunc
+  real(dep_prec) :: drygridtotalunc,wetgridtotalunc,drydeposit(maxspec)
+  real :: xold,yold,zold,xmassfract
   real, parameter :: e_inv = 1.0/exp(1.0)
 !double precision xm(maxspec,maxpointspec_act),
@@ -261,5 +262,5 @@
 ! TODO: z0(7) changes with time, so should be dimension (numclass,2) to
 ! allow transfer of the future value in the background
-        call MPI_Bcast(z0,numclass,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
+        call MPI_Bcast(z0,numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
         call mpif_gf_request
       end if
@@ -756,5 +757,5 @@
                    
                    !CGZ-lifetime: Check mass fraction left/save lifetime
-                   ! if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.inv_e.and.checklifetime(j,ks).eq.0.)then
+                   ! if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.e_inv.and.checklifetime(j,ks).eq.0.)then
                        !Mass below 1% of initial >register lifetime
                    !     checklifetime(j,ks)=abs(itra1(j)-itramem(j))
Index: src/unc_mod.f90
===================================================================
--- src/unc_mod.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/unc_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -22,5 +22,5 @@
 module unc_mod
 
-  use par_mod, only:cp
+  use par_mod, only:dep_prec
 
   implicit none
@@ -28,14 +28,14 @@
   real,allocatable ,dimension (:,:,:,:,:,:,:) :: gridunc
   real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygridunc
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgridunc
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: wetgridunc
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn
 
 ! For sum of individual contributions, used for the MPI version
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygridunc0
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn0
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgridunc0
-  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn0
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc0
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn0
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: wetgridunc0
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn0
 
   real,allocatable, dimension (:,:,:,:,:) :: init_cond
Index: src/verttransform.f90
===================================================================
--- src/verttransform.f90	(revision 41d8574bea1b0f91e61dffd5dd1d17882407f826)
+++ src/verttransform.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -772,5 +772,5 @@
 !          if (lsp.gt.convp) then !  prectype='lsp'
 !            lconvectprec = .false.
-!          else ! prectype='cp '
+!          else ! prectype='cp'
 !            lconvectprec = .true.
 !           endif
