Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision 0c00f1f8af5469692bfe795d126d70d5205ca6f6)
+++ src/FLEXPART.f90	(revision a9cf4b1e763458eb1d7f20f9d55d1f5bd8a87da6)
@@ -52,7 +52,11 @@
   use com_mod
   use conv_mod
-  use netcdf_output_mod, only: writeheader_netcdf
+
   use random_mod, only: gasdev1
   use class_gribfile
+
+#ifdef USE_NCF
+  use netcdf_output_mod, only: writeheader_netcdf
+#endif
 
   implicit none
@@ -352,4 +356,5 @@
   !******************************************************************
 
+#ifdef USE_NCF
   if (lnetcdfout.eq.1) then 
     call writeheader_netcdf(lnest=.false.)
@@ -365,4 +370,5 @@
     endif
   endif
+#endif
 
   if (verbosity.gt.0) then
Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision 0c00f1f8af5469692bfe795d126d70d5205ca6f6)
+++ src/FLEXPART_MPI.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -53,7 +53,10 @@
   use conv_mod
   use mpi_mod
-  use netcdf_output_mod, only: writeheader_netcdf
   use random_mod, only: gasdev1
   use class_gribfile
+
+#ifdef USE_NCF
+  use netcdf_output_mod, only: writeheader_netcdf
+#endif
 
   implicit none
@@ -64,5 +67,6 @@
   integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
   integer :: detectformat
-
+  integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
+       & tot_i=0
 
 
@@ -203,9 +207,9 @@
 
   if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
-    print *,'ECMWF metdata detected'
+    if (lroot) print *,'ECMWF metdata detected'
   elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
-    print *,'NCEP metdata detected'
+    if (lroot) print *,'NCEP metdata detected'
   else
-    print *,'Unknown metdata format'
+    if (lroot) print *,'Unknown metdata format'
     stop
   endif
@@ -378,21 +382,22 @@
 
   if (mp_measure_time) call mpif_mtime('iotime',0)
+
   if (lroot) then ! MPI: this part root process only
-
-  if (lnetcdfout.eq.1) then 
-    call writeheader_netcdf(lnest=.false.)
-  else 
-    call writeheader
-  end if
-
-  if (nested_output.eq.1) then
-    if (lnetcdfout.eq.1) then
-      call writeheader_netcdf(lnest=.true.)
-    else
-      call writeheader_nest
-    endif
-  endif
-
-!
+#ifdef USE_NCF
+    if (lnetcdfout.eq.1) then 
+      call writeheader_netcdf(lnest=.false.)
+    else 
+      call writeheader
+    end if
+    
+    if (nested_output.eq.1) then
+      if (lnetcdfout.eq.1) then
+        call writeheader_netcdf(lnest=.true.)
+      else
+        call writeheader_nest
+      endif
+    endif
+#endif
+
     if (verbosity.gt.0) then
       print*,'call writeheader'
@@ -402,5 +407,5 @@
 ! FLEXPART 9.2 ticket ?? write header in ASCII format 
     call writeheader_txt
-!if (nested_output.eq.1) call writeheader_nest
+
     if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
     if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
@@ -409,6 +414,4 @@
 
   if (mp_measure_time) call mpif_mtime('iotime',0)
-
-  !open(unitdates,file=path(2)(1:length(2))//'dates')
 
   if (verbosity.gt.0 .and. lroot) then
@@ -481,17 +484,11 @@
 
 ! NIK 16.02.2005 
-  if (lroot) then
-    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
+  if (mp_partgroup_pid.ge.0) then ! Skip for readwind process 
+    call MPI_Reduce(tot_blc_count, tot_b, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
          & mp_comm_used, mp_ierr)
-    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
+    call MPI_Reduce(tot_inc_count, tot_i, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
          & mp_comm_used, mp_ierr)
-  else
-    if (mp_partgroup_pid.ge.0) then ! Skip for readwind process 
-      call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
-           & mp_comm_used, mp_ierr)
-      call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
-           & mp_comm_used, mp_ierr)
-    end if
   end if
+  
 
   if (lroot) then
@@ -500,7 +497,9 @@
       write(*,*) 'Scavenging statistics for species ', species(i), ':'
       write(*,*) 'Total number of occurences of below-cloud scavenging', &
-           & tot_blc_count(i)
+           & tot_b(i)
+!           & tot_blc_count(i)
       write(*,*) 'Total number of occurences of in-cloud    scavenging', &
-           & tot_inc_count(i)
+           & tot_i(i)
+!           & tot_inc_count(i)
       write(*,*) '**********************************************'
     end do
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision 08a38b55d3749abb3d1df461e12c98d18b329a51)
+++ src/concoutput.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -140,4 +140,6 @@
   close(unitdates)  
 
+  ! Overwrite existing dates file on first call, later append to it
+  ! This fixes a bug where the dates file kept growing across multiple runs
   IF (init) THEN
     file_stat='OLD'
@@ -313,5 +315,5 @@
                    drygridsigma(ix,jy)* &
                    sqrt(real(nclassunc))
-125           drygridsigmatotal=drygridsigmatotal+ &
+              drygridsigmatotal=drygridsigmatotal+ &
                    drygridsigma(ix,jy)
             endif
Index: src/concoutput_mpi.f90
===================================================================
--- src/concoutput_mpi.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/concoutput_mpi.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -104,5 +104,5 @@
   real,parameter :: weightair=28.97
   logical :: sp_zer
-  LOGICAL,save :: init=.true.
+  logical,save :: init=.true.
   character :: adate*8,atime*6
   character(len=3) :: anspec
@@ -113,5 +113,4 @@
   integer :: ierr
   character(LEN=100) :: dates_char
-!  character :: dates_char
 
 ! Measure execution time
@@ -129,5 +128,5 @@
 ! This fixes a bug where the dates file kept growing across multiple runs
 
-! If 'dates' file exists, make a backup
+! If 'dates' file exists in output directory, make a backup
   inquire(file=path(2)(1:length(2))//'dates', exist=ldates_file)
   if (ldates_file.and.init) then
@@ -258,21 +257,32 @@
 
     write(anspec,'(i3.3)') ks
-    if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
-      if (ldirect.eq.1) then
-        open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
+
+    if (DRYBKDEP.or.WETBKDEP) then !scavdep output
+      if (DRYBKDEP) & 
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      if (WETBKDEP) & 
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      write(unitoutgrid) itime
+    else
+      if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+        if (ldirect.eq.1) then
+          open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
+               atime//'_'//anspec,form='unformatted')
+        else
+          open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
+               atime//'_'//anspec,form='unformatted')
+        endif
+        write(unitoutgrid) itime
+      endif
+
+      if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
+        open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
              atime//'_'//anspec,form='unformatted')
-      else
-        open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
-             atime//'_'//anspec,form='unformatted')
+
+        write(unitoutgridppt) itime
       endif
-      write(unitoutgrid) itime
-    endif
-
-    if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
-      open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
-           atime//'_'//anspec,form='unformatted')
-
-      write(unitoutgridppt) itime
-    endif
+    endif ! if deposition output
 
     do kp=1,maxpointspec_act
@@ -354,5 +364,5 @@
 ! Concentration output
 !*********************
-        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
 
 ! Wet deposition
@@ -449,8 +459,16 @@
                   endif
                   sp_count_r=sp_count_r+1
+                  if (lparticlecountoutput) then
+                    sparse_dump_r(sp_count_r)= &
+                         sp_fact* &
+                         grid(ix,jy,kz)
+                  else
                   sparse_dump_r(sp_count_r)= &
                        sp_fact* &
                        grid(ix,jy,kz)* &
                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
+                  end if
+
+
 !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
 !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
@@ -638,22 +656,6 @@
 !*************************
 
-  do ks=1,nspec
-    do kp=1,maxpointspec_act
-      do i=1,numreceptor
-        creceptor(i,ks)=0.
-      end do
-      do jy=0,numygrid-1
-        do ix=0,numxgrid-1
-          do l=1,nclassunc
-            do nage=1,nageclass
-              do kz=1,numzgrid
-                gridunc(ix,jy,kz,ks,kp,l,nage)=0.
-              end do
-            end do
-          end do
-        end do
-      end do
-    end do
-  end do
+  creceptor(:,:)=0.
+  gridunc(:,:,:,:,:,:,:)=0.
 
   if (mp_measure_time) call mpif_mtime('rootonly',1)
Index: src/getfields_mpi.f90
===================================================================
--- src/getfields_mpi.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
+++ src/getfields_mpi.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -233,4 +233,8 @@
 40  indmin=indj
 
+   if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
+        call writeprecip(itime,memind(1))
+      endif
+
   else
 
@@ -287,5 +291,8 @@
 60  indmin=indj
 
-    mind3=memstat
+!   if (WETBKDEP.and.lroot) then
+    if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
+      call writeprecip(itime,memind(1))
+    endif
 
   endif
Index: src/makefile
===================================================================
--- src/makefile	(revision ccc9ec26732116bd18049809fcb9b319d3148f23)
+++ src/makefile	(revision a9cf4b1e763458eb1d7f20f9d55d1f5bd8a87da6)
@@ -31,4 +31,8 @@
 #      make [-j] mpi-dbg
 #
+#  NETCDF OUTPUT
+#    To add support for output in netCDF format, append `ncf=yes` to the
+#    `make` command
+# 
 ################################################################################
 
@@ -57,6 +61,5 @@
 	INCPATH2  = ${ROOT_DIR}/include
 	LIBPATH1 = ${ROOT_DIR}/lib
-
-else #ifeq ($(gcc), 5.4)
+else 
 # Compiled libraries under user ~flexpart, gfortran v5.4
 	ROOT_DIR = /homevip/flexpart/
@@ -68,14 +71,14 @@
 	INCPATH2  = /usr/include
 	LIBPATH1 = ${ROOT_DIR}/gcc-5.4.0/lib
-
-#else
-# Default: System libraries at NILU, gfortran v4.6
-#	F90       = /usr/bin/gfortran
-#	MPIF90    = /usr/bin/mpif90.openmpi
-
-#	INCPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/include
-#	INCPATH2 = /usr/include
-#	LIBPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/lib
 endif
+
+
+### Enable netCDF output?
+ifeq ($(ncf), yes)
+	NCOPT = -DUSE_NCF -lnetcdff	
+else
+	NCOPT = -UUSE_NCF
+endif
+
 
 
@@ -89,9 +92,10 @@
 
 ## LIBRARIES
-LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
-
-FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER)  #-Warray-bounds -fcheck=all # -march=native
-
-DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)  #  -ffpe-trap=invalid,overflow,denormal,underflow,zero  -Warray-bounds -fcheck=all
+#LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff
+LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper $(NCOPT)
+
+FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(NCOPT) $(FUSER)  #-Warray-bounds -fcheck=all # -march=native
+
+DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) $(NCOPT) -fbacktrace   -Wall  -fdump-core $(FUSER)  #  -ffpe-trap=invalid,overflow,denormal,underflow,zero  -Warray-bounds -fcheck=all
 
 LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
@@ -139,4 +143,6 @@
 	getfields_mpi.o \
         readwind_ecmwf_mpi.o 
+
+OBJECTS_NCF = netcdf_output_mod.o
 
 OBJECTS = \
@@ -196,5 +202,9 @@
 dynamic_viscosity.o     get_settling.o	\
 initialize_cbl_vel.o	re_initialize_particle.o \
-cbl.o			netcdf_output_mod.o
+cbl.o
+
+ifeq ($(ncf), yes)
+	OBJECTS	:= $(OBJECTS) $(OBJECTS_NCF)
+endif
 
 %.o: %.mod
Index: src/readOHfield.f90
===================================================================
--- src/readOHfield.f90	(revision d7aab4b60376261b6f48658c75153008dfba6bfd)
+++ src/readOHfield.f90	(revision 3f149ccfb522ead3fdd1e69283deccd58b8740fc)
@@ -22,23 +22,25 @@
 subroutine readOHfield
 
-  !*****************************************************************************
-  !                                                                            *
-  ! Reads the OH field into memory                                             *
-  !                                                                            *
-  ! AUTHOR: R.L. Thompson, Nov 2014                                            *
-  !                                                                            *
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  !                                                                            *
-  ! path(numpath)              contains the path names                         *
-  ! lonOH(nxOH)                longitude of OH fields                          *
-  ! latOH(nyOH)                latitude of OH fields                           *
-  ! altOH(nzOH)                altitude of OH fields                           *
-  ! etaOH(nzOH)                eta-levels of OH fields                         *
-  ! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3)                *
-  !                                                                            *
-  !                                                                            *
-  !*****************************************************************************
+!*****************************************************************************
+!                                                                            *
+! Reads the OH field into memory                                             *
+!                                                                            *
+! AUTHOR: R.L. Thompson, Nov 2014                                            *
+!                                                                            *
+! UPDATES:                                                                   *
+!   03/2018 SEC: Converted original netCDF files to binary format            *
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+!                                                                            *
+! path(numpath)              contains the path names                         *
+! lonOH(nxOH)                longitude of OH fields                          *
+! latOH(nyOH)                latitude of OH fields                           *
+! altOH(nzOH)                altitude of OH fields                           *
+! etaOH(nzOH)                eta-levels of OH fields                         *
+! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3)                *
+!                                                                            *
+!                                                                            *
+!*****************************************************************************
 
   use oh_mod
@@ -48,9 +50,5 @@
   implicit none
 
-  include 'netcdf.inc'
-
-  character(len=150) :: thefile
-  character(len=2) :: mm
-  integer :: nid,ierr,xid,yid,zid,vid,m
+  integer :: i,j,k,l,ierr
   real, dimension(:), allocatable :: etaOH
 
@@ -61,142 +59,37 @@
   real, parameter :: scalehgt=7000. ! scale height in metres
 
-  ! Read OH fields and level heights
-  !********************************
 
-  do m=1,12
-  
-    ! open netcdf file
-    write(mm,fmt='(i2.2)') m
-!    thefile=trim(path(1))//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
-    thefile=trim(ohfields_path)//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
-    ierr=nf_open(trim(thefile),NF_NOWRITE,nid)
-    if(ierr.ne.0) then
-      write (*,*) 'The OH field at: '//thefile// ' could not be opened'
-      write (*,*) 'please copy the OH fields there, or change the path in the'
-      write (*,*) 'COMMAND namelist!'
-      write(*,*) nf_strerror(ierr)
-      stop
-    endif
+  open(unitOH,file=trim(ohfields_path) &
+       //'OH_FIELDS/OH_variables.bin',status='old', &
+       form='UNFORMATTED', iostat=ierr, convert='little_endian')
 
-    ! inquire about variables
-    ierr=nf_inq_dimid(nid,'Lon-000',xid)
-    if(ierr.ne.0) then
-      write(*,*) nf_strerror(ierr)
-      stop
-    endif
-    ierr=nf_inq_dimid(nid,'Lat-000',yid)
-    if(ierr.ne.0) then
-      write(*,*) nf_strerror(ierr)
-      stop
-    endif
-    ierr=nf_inq_dimid(nid,'Alt-000',zid)
-    if(ierr.ne.0) then
-      write(*,*) nf_strerror(ierr)
-      stop
-    endif
-
-    if(m.eq.1) then
-
-      ! read dimension sizes
-      ierr=nf_inq_dimlen(nid,xid,nxOH)
-      if(ierr.ne.0) then
-        write(*,*) nf_strerror(ierr)
-        stop
-      endif
-      ierr=nf_inq_dimlen(nid,yid,nyOH)
-      if(ierr.ne.0) then
-        write(*,*) nf_strerror(ierr)
-        stop
-      endif
-      ierr=nf_inq_dimlen(nid,zid,nzOH)
-      if(ierr.ne.0) then
-        write(*,*) nf_strerror(ierr)
-        stop
-      endif  
-
-      ! allocate variables
-      allocate(lonOH(nxOH))
-      allocate(latOH(nyOH))
-      allocate(etaOH(nzOH))
-      allocate(altOH(nzOH))
-      allocate(OH_field(nxOH,nyOH,nzOH,12))
-      allocate(OH_hourly(nxOH,nyOH,nzOH,2))
-
-      ! read dimension variables
-      ierr=nf_inq_varid(nid,'LON',xid)
-      ierr=nf_get_var_real(nid,xid,lonOH)
-      if(ierr.ne.0) then
-        write(*,*) nf_strerror(ierr)
-        stop
-      endif
-      ierr=nf_inq_varid(nid,'LAT',yid)
-      ierr=nf_get_var_real(nid,yid,latOH)
-      if(ierr.ne.0) then
-        write(*,*) nf_strerror(ierr)
-        stop
-      endif
-      ierr=nf_inq_varid(nid,'ETAC',zid)
-      ierr=nf_get_var_real(nid,zid,etaOH)
-      if(ierr.ne.0) then
-        write(*,*) nf_strerror(ierr)
-        stop
-      endif
-
-      ! convert eta-level to altitude (assume surface pressure of 1010 hPa)
-      altOH=log(1010./(etaOH*1010.))*scalehgt
-
-    endif ! m.eq.1
-
-    ! read OH_field
-    ierr=nf_inq_varid(nid,'CHEM-L_S__OH',vid)
-    ierr=nf_get_var_real(nid,vid,OH_field(:,:,:,m))
-    if(ierr.ne.0) then
-      write(*,*) nf_strerror(ierr)
-      stop
-    endif
-
-    ierr=nf_close(nid)
-
-  end do 
- 
-  deallocate(etaOH)
-
-  ! Read J(O1D) photolysis rates
-  !********************************  
-
-  ! open netcdf file
-!  thefile=trim(path(1))//'OH_FIELDS/jrate_average.nc'
-  thefile=trim(ohfields_path)//'OH_FIELDS/jrate_average.nc'
-  ierr=nf_open(trim(thefile),NF_NOWRITE,nid)
   if(ierr.ne.0) then
-    write(*,*) nf_strerror(ierr)
+    write(*,*) 'Cannot read binary OH fields in ',trim(ohfields_path)//'OH_FIELDS/OH_variables.bin'
     stop
   endif
 
-  ! read dimension variables
-  ierr=nf_inq_varid(nid,'longitude',xid)
-  ierr=nf_get_var_real(nid,xid,lonjr)
-  if(ierr.ne.0) then
-    write(*,*) nf_strerror(ierr)
-    stop
-  endif
-  ierr=nf_inq_varid(nid,'latitude',yid)
-  ierr=nf_get_var_real(nid,yid,latjr)
-  if(ierr.ne.0) then
-    write(*,*) nf_strerror(ierr)
-    stop
-  endif
+  read(unitOH) nxOH
+  read(unitOH) nyOH
+  read(unitOH) nzOH
+  write(*,*) nxOH,nyOH,nzOH
 
-  ! read jrate_average
-  ierr=nf_inq_varid(nid,'jrate',vid)
-  ierr=nf_get_var_real(nid,vid,jrate_average)
-  if(ierr.ne.0) then
-    write(*,*) nf_strerror(ierr)
-    stop
-  endif
+! allocate variables
+  allocate(lonOH(nxOH))
+  allocate(latOH(nyOH))
+  allocate(etaOH(nzOH))
+  allocate(altOH(nzOH))
+  allocate(OH_field(nxOH,nyOH,nzOH,12))
+  allocate(OH_hourly(nxOH,nyOH,nzOH,2))
 
-  ierr=nf_close(nid) 
+  read(unitOH) (lonjr(i),i=1,360)
+  read(unitOH) (latjr(i),i=1,180)
+  read(unitOH) (((jrate_average(i,j,k),i=1,360),j=1,180),k=1,12)
+  read(unitOH) (lonOH(i),i=1,nxOH)
+  read(unitOH) (latOH(i),i=1,nyOH)
+  read(unitOH) (lonOH(i),i=1,nxOH)
 
-  return
+  read(unitOH) (altOH(i),i=1,nzOH)
+  read(unitOH) ((((OH_field(i,j,k,l),i=1,nxOH),j=1,nyOH),k=1,nzOH),l=1,12)
+  read(unitOH) ((((OH_hourly(i,j,k,l),i=1,nxOH),j=1,nyOH),k=1,nzOH),l=1,2)
 
 end subroutine readOHfield
Index: src/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision 01e60525c4830ea98f0acaa4ca0894f3845020a9)
+++ src/readcommand.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -332,14 +332,18 @@
      case (3)  ! 3 .. wet deposition in outputfield 
         ind_rel = 3
-         write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
-         write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
-         write(*,*) ' #### Release is performed above ground lev    #### '
+        if (lroot) then
+          write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
+          write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
+          write(*,*) ' #### Release is performed above ground lev    #### '
+        end if
          WETBKDEP=.true.
          allocate(xscav_frac1(maxpart,maxspec))
      case (4)  ! 4 .. dry deposition in outputfield
          ind_rel = 4
-         write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE    #### '
-         write(*,*) ' #### Releaseheight is forced to 0 - 2*href    #### '
-         write(*,*) ' #### Release is performed above ground lev    #### '
+         if (lroot) then
+           write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE    #### '
+           write(*,*) ' #### Releaseheight is forced to 0 - 2*href    #### '
+           write(*,*) ' #### Release is performed above ground lev    #### '
+         end if
          DRYBKDEP=.true.
          allocate(xscav_frac1(maxpart,maxspec))
@@ -392,13 +396,14 @@
   endif
 
-!  check for netcdf output switch (use for non-namelist input only!)
+! Check for netcdf output switch
+!*******************************
   if (iout.ge.8) then
      lnetcdfout = 1
      iout = iout - 8
-! #ifndef NETCDF_OUTPUT
-!      print*,'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
-!      print*,'Please recompile with netcdf library or use standard output format.'
-!      stop
-! #endif
+#ifndef USE_NCF
+     write(*,*) 'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
+     write(*,*) 'Please recompile with netcdf library (`make [...] ncf=yes`) or use standard output format.'
+     stop
+#endif
   endif
 
Index: src/readwind_ecmwf_mpi.f90
===================================================================
--- src/readwind_ecmwf_mpi.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/readwind_ecmwf_mpi.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -103,10 +103,10 @@
 
   integer :: isec1(56),isec2(22+nxmax+nymax)
-  real(kind=4) :: zsec4(jpunp)
-  real(kind=4) :: xaux,yaux
-  real(kind=8) :: xauxin,yauxin
-  real,parameter :: eps=1.e-4
-  real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
-  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
+  real(sp) :: zsec4(jpunp)
+  real(sp) :: xaux,yaux
+  real(dp) :: xauxin,yauxin
+  real(sp),parameter :: eps=1.e-4
+  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
+  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
 
   logical :: hflswitch,strswitch !,readcloud
@@ -123,5 +123,5 @@
 !ZHG test the grib fields that have lcwc without using them
 !  readcloud=.false.
-!hg end
+
   levdiff2=nlev_ec-nwz+1
   iumax=0
Index: src/releaseparticles_mpi.f90
===================================================================
--- src/releaseparticles_mpi.f90	(revision 16b61a5e2d3b1020be1cff080e479ba26f3f7063)
+++ src/releaseparticles_mpi.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -216,5 +216,10 @@
               xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
                    *timecorrect(k)/average_timecorrect
-!             write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
+              if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
+!              if ( henry(k).gt.0 .or. &
+!                   crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
+!                   ccn_aero(k).gt.0. .or. in_aero(k).gt.0. )  then
+                xscav_frac1(ipart,k)=-1.
+               endif
 ! Assign certain properties to particle
 !**************************************
@@ -366,5 +371,5 @@
 !Af ind_rel is defined in readcommand.f
 
-            if (ind_rel .eq. 1) then
+            if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
 
 ! Interpolate the air density
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision fe32dca5c7e5a54238ac82b1dab620991c979f5c)
+++ src/timemanager.f90	(revision a9cf4b1e763458eb1d7f20f9d55d1f5bd8a87da6)
@@ -101,6 +101,8 @@
   use par_mod
   use com_mod
+#ifdef USE_NCF
   use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
        &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
+#endif
 
   implicit none
@@ -389,5 +391,7 @@
           if (surf_only.ne.1) then 
             if (lnetcdfout.eq.1) then 
+#ifdef USE_NCF
               call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+#endif
             else 
               call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
@@ -399,6 +403,8 @@
              write(*,*) 'system clock',count_clock - count_clock0   
             endif
-            if (lnetcdfout.eq.1) then 
+            if (lnetcdfout.eq.1) then
+#ifdef USE_NCF
               call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+#endif
             else
               call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
@@ -419,4 +425,5 @@
               endif
             else
+#ifdef USE_NCF
               if (surf_only.ne.1) then
                 call concoutput_nest_netcdf(itime,outnum)
@@ -424,4 +431,5 @@
                 call concoutput_surf_nest_netcdf(itime,outnum)
               endif
+#endif
             endif
           endif
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision d8eed026ebec594e9a91de6c6fbc45b354a0e4b2)
+++ src/timemanager_mpi.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -102,6 +102,8 @@
   use com_mod
   use mpi_mod
+#ifdef USE_NCF
   use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
        &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
+#endif
 
   implicit none
@@ -113,9 +115,8 @@
   integer :: ip
   integer :: loutnext,loutstart,loutend
-  integer :: ix,jy,ldeltat,itage,nage
+  integer :: ix,jy,ldeltat,itage,nage,idummy
   integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme 
   integer :: numpart_tot_mpi ! for summing particles on all processes
-  real :: outnum,weight,prob(maxspec)
-  real :: decfact
+  real :: outnum,weight,prob(maxspec), prob_rec(maxspec), decfact,wetscav
 
   real(sp) :: gridtotalunc
@@ -123,4 +124,5 @@
        & drydeposit(maxspec)=0_dep_prec
   real :: xold,yold,zold,xmassfract
+  real :: grfraction(3)
   real, parameter :: e_inv = 1.0/exp(1.0)
 
@@ -157,4 +159,6 @@
 !CGZ-lifetime: set lifetime to 0
 
+  if (.not.lusekerneloutput) write(*,*) 'Not using the kernel'
+  if (turboff) write(*,*) 'Turbulence switched off'
 
 
@@ -481,6 +485,8 @@
             if (lroot) then
               if (lnetcdfout.eq.1) then 
+#ifdef USE_NCF
                 call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
                      &drygridtotalunc)
+#endif
               else 
                 call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
@@ -494,6 +500,8 @@
             if (lroot) then
               if (lnetcdfout.eq.1) then
+#ifdef USE_NCF
                 call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
                      &drygridtotalunc)
+#endif
               else
                 call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
@@ -513,5 +521,5 @@
             call mpif_tm_reduce_grid_nest
  
-           if (mp_measure_time) call mpif_mtime('iotime',0)
+            if (mp_measure_time) call mpif_mtime('iotime',0)
 
             if (lnetcdfout.eq.0) then
@@ -526,9 +534,7 @@
               else  ! :TODO: check for zeroing in the netcdf module
                 call concoutput_surf_nest(itime,outnum)
-
               end if
-
             else
-
+#ifdef USE_NCF
               if (surf_only.ne.1) then
                 if (lroot) then              
@@ -544,10 +550,7 @@
                 end if
               endif
-
-
+#endif
             end if
           end if
-          
-
           outnum=0.
         endif
@@ -707,4 +710,40 @@
         zold=ztra1(j)
 
+   
+  ! RECEPTOR: dry/wet depovel
+  !****************************
+  ! Before the particle is moved 
+  ! the calculation of the scavenged mass shall only be done once after release
+  ! xscav_frac1 was initialised with a negative value
+
+      if  (DRYBKDEP) then
+       do ks=1,nspec
+         if  ((xscav_frac1(j,ks).lt.0)) then
+            call get_vdep_prob(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
+            if (DRYDEPSPEC(ks)) then        ! dry deposition
+               xscav_frac1(j,ks)=prob_rec(ks)
+             else
+                xmass1(j,ks)=0.
+                xscav_frac1(j,ks)=0.
+             endif
+         endif
+        enddo
+       endif
+
+       if (WETBKDEP) then 
+       do ks=1,nspec
+         if  ((xscav_frac1(j,ks).lt.0)) then
+            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
+            if (wetscav.gt.0) then
+                xscav_frac1(j,ks)=wetscav* &
+                       (zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
+            else
+                xmass1(j,ks)=0.
+                xscav_frac1(j,ks)=0.
+            endif
+         endif
+        enddo
+       endif
+
 ! Integrate Lagevin equation for lsynctime seconds
 !*************************************************
@@ -808,5 +847,7 @@
                  call initial_cond_calc(itime+lsynctime,j)
             itra1(j)=-999999999
-            !print*, 'terminated particle ',j,'for age'
+            if (verbosity.gt.0) then
+              print*, 'terminated particle ',j,'for age'
+            endif
           endif
         endif
@@ -819,7 +860,7 @@
 
 
-! Added by mc: counter of "unstable" particle velocity during a time scale
-!   of maximumtl=20 minutes (defined in com_mod)
-
+! Counter of "unstable" particle velocity during a time scale
+! of maximumtl=20 minutes (defined in com_mod)
+!************************************************************
     total_nan_intl=0
     i_nan=i_nan+1 ! added by mc to count nan during a time of maxtl (i.e. maximum tl fixed here to 20 minutes, see com_mod)
