Index: options/COMMAND
===================================================================
--- options/COMMAND	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ options/COMMAND	(revision 6d73c4be98f6a308eed735650a968ac230fd19ef)
@@ -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 3]time averaged
+ IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end   
  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 2eefa58ac797f119acc28fade41beff352006040)
+++ src/FLEXPART.f90	(revision 50958b85292a7b2ce570403ca6f859592eaa1b5e)
@@ -68,5 +68,11 @@
   integer :: detectformat
 
+
+
+  ! Initialize arrays in com_mod
+  !*****************************
+  call com_mod_allocate_part(maxpart)
   
+
   ! Generate a large number of random numbers
   !******************************************
@@ -137,5 +143,5 @@
     write(*,*) 'call readpaths'
   endif 
-  call readpaths
+  call readpaths(pathfile)
  
   if (verbosity.gt.1) then !show clock info 
@@ -165,9 +171,4 @@
     endif     
   endif
-
-  ! Initialize arrays in com_mod
-  !*****************************
-  call com_mod_allocate_part(maxpart)
-
 
   ! Read the age classes to be used
@@ -452,7 +453,5 @@
   endif
 
-  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
   call timemanager(metdata_format)
- 
 
   if (verbosity.gt.0) then
@@ -469,4 +468,5 @@
       endif
     end do
+    write (*,*) 'timemanager> call wetdepo'
   endif
   
Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/FLEXPART_MPI.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -77,5 +77,10 @@
   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
   !******************************************
@@ -146,5 +151,5 @@
     write(*,*) 'call readpaths'
   endif 
-  call readpaths
+  call readpaths(pathfile)
  
   if (verbosity.gt.1) then !show clock info 
@@ -174,9 +179,4 @@
     endif
   endif
-
-  ! Initialize arrays in com_mod 
-  !*****************************
-
-  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
 
 
@@ -413,5 +413,5 @@
   end if ! (mpif_pid == 0) 
 
-  if (mp_measure_time) call mpif_mtime('iotime',1)
+  if (mp_measure_time) call mpif_mtime('iotime',0)
 
   if (verbosity.gt.0 .and. lroot) then
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/com_mod.f90	(revision e9e0f0627da3a042c02a5352bf3446bdfbcfa136)
@@ -18,6 +18,4 @@
 
   implicit none
-
-
 
   !****************************************************************
@@ -71,5 +69,5 @@
 
   real :: ctl,fine
-  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill,ipoutfac
+  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
   integer :: mquasilag,nested_output,ind_source,ind_receptor
   integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
@@ -84,5 +82,4 @@
   ! 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)
@@ -124,7 +121,4 @@
   ! lnetcdfout   1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
 
-  integer :: linversionout
-  ! linversionout 1 for one grid_time output file for each release containing all timesteps
-
   integer :: nageclass,lage(maxageclass)
 
@@ -134,4 +128,5 @@
 
   logical :: gdomainfill
+
   ! gdomainfill             .T., if domain-filling is global, .F. if not
 
@@ -179,5 +174,5 @@
   real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass)
   real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass)
-  real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
+  real :: rm(maxspec),dryvel(maxspec)
   real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
 
@@ -364,7 +359,5 @@
   real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice      [kg/kg]
   real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0  !combined [m3/m3]
-! RLT add pressure and dry air density
-  real :: prs(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
-  real :: rho_dry(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
+
   real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
@@ -388,5 +381,4 @@
   ! uupol,vvpol [m/s]    wind components in polar stereographic projection
   ! tt [K]               temperature data
-  ! prs                  air pressure
   ! qv                   specific humidity data
   ! pv (pvu)             potential vorticity
@@ -659,5 +651,4 @@
   real :: receptorarea(maxreceptor)
   real :: creceptor(maxreceptor,maxspec)
-  real, allocatable, dimension(:,:) :: creceptor0
   character(len=16) :: receptorname(maxreceptor)
   integer :: numreceptor
@@ -682,12 +673,4 @@
   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
@@ -797,19 +780,11 @@
          & idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
          & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
-         & 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
+         & xmass1(nmpart, maxspec),&
+         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
 
 
     allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
          & vs(nmpart),ws(nmpart),cbt(nmpart))
-
+    
   end subroutine com_mod_allocate_part
 
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/concoutput.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -72,7 +72,4 @@
   real :: sp_fact
   real :: outnum,densityoutrecept(maxreceptor),xl,yl
-! RLT
-  real :: densitydryrecept(maxreceptor)
-  real :: factor_dryrecept(maxreceptor)
 
 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
@@ -109,5 +106,4 @@
   character(LEN=8),save :: file_stat='REPLACE'
   logical :: ldates_file
-  logical :: lexist
   integer :: ierr
   character(LEN=100) :: dates_char
@@ -207,7 +203,4 @@
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
              rho(iix,jjy,kzz-1,mind)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz  
       end do
     end do
@@ -221,12 +214,6 @@
 !densityoutrecept(i)=rho(iix,jjy,1,2)
     densityoutrecept(i)=rho(iix,jjy,1,mind)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
   end do
 
-! RLT
-! conversion factor for output relative to dry air
-  factor_drygrid=densityoutgrid/densitydrygrid
-  factor_dryrecept=densityoutrecept/densitydryrecept
 
 ! Output is different for forward and backward simulations
@@ -366,5 +353,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
@@ -627,47 +614,4 @@
   end do
 
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
-  if (lexist) then
-    ! open and append
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
-            status='old',action='write',access='append')
-  else
-    ! create new
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
-            status='new',action='write')
-  endif
-  sp_count_i=0
-  sp_count_r=0
-  sp_fact=-1.
-  sp_zer=.true.
-  do kz=1,numzgrid
-    do jy=0,numygrid-1
-      do ix=0,numxgrid-1
-        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
-          if (sp_zer.eqv..true.) then ! first value not equal to one
-            sp_count_i=sp_count_i+1
-            sparse_dump_i(sp_count_i)= &
-                  ix+jy*numxgrid+kz*numxgrid*numygrid
-            sp_zer=.false.
-            sp_fact=sp_fact*(-1.)
-          endif
-          sp_count_r=sp_count_r+1
-          sparse_dump_r(sp_count_r)= &
-               sp_fact*factor_drygrid(ix,jy,kz)
-        else ! factor is one
-          sp_zer=.true.
-        endif
-      end do
-    end do
-  end do
-  write(unitoutfactor) sp_count_i
-  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
-  write(unitoutfactor) sp_count_r
-  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
-  close(unitoutfactor)
-
-
   if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
   if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
@@ -696,21 +640,5 @@
   endif
 
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  if (numreceptor.gt.0) then
-    inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
-     if (lexist) then
-     ! open and append
-      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
-              status='old',action='write',access='append')
-    else
-      ! create new
-      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
-              status='new',action='write')
-    endif
-    write(unitoutfactor) itime
-    write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
-    close(unitoutfactor)
-  endif
+
 
 ! Reinitialization of grid
Index: c/concoutput_inversion.f90
===================================================================
--- src/concoutput_inversion.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ 	(revision )
@@ -1,855 +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 concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc, &
-     drygridtotalunc)
-!                        i     i          o             o
-!       o
-!*****************************************************************************
-!                                                                            *
-!     Output of the concentration grid and the receptor concentrations.      *
-!                                                                            *
-!     Author: A. Stohl                                                       *
-!                                                                            *
-!     24 May 1995                                                            *
-!                                                                            *
-!     13 April 1999, Major update: if output size is smaller, dump output    *
-!                    in sparse matrix format; additional output of           *
-!                    uncertainty                                             *
-!                                                                            *
-!     05 April 2000, Major update: output of age classes; output for backward*
-!                    runs is time spent in grid cell times total mass of     *
-!                    species.                                                *
-!                                                                            *
-!     17 February 2002, Appropriate dimensions for backward and forward runs *
-!                       are now specified in file par_mod                    *
-!                                                                            *
-!     June 2006, write grid in sparse matrix with a single write command     *
-!                in order to save disk space                                 *
-!                                                                            *
-!     2008 new sparse matrix format                                          *
-!
-!     January 2017,  Separate files by release but include all timesteps
-!                                                                            *
-!*****************************************************************************
-!                                                                            *
-! Variables:                                                                 *
-! outnum          number of samples                                          *
-! ncells          number of cells with non-zero concentrations               *
-! sparse          .true. if in sparse matrix format, else .false.            *
-! tot_mu          1 for forward, initial mass mixing ration for backw. runs  *
-!                                                                            *
-!*****************************************************************************
-
-  use unc_mod
-  use point_mod
-  use outg_mod
-  use par_mod
-  use com_mod
-  use mean_mod
-
-  implicit none
-
-  real(kind=dp) :: jul
-  integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
-  integer :: sp_count_i,sp_count_r
-  real :: sp_fact
-  real :: outnum,densityoutrecept(maxreceptor),xl,yl
-! RLT
-  real :: densitydryrecept(maxreceptor)
-  real :: factor_dryrecept(maxreceptor)
-
-!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
-!    +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
-!    +    maxageclass)
-!real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
-!    +       maxageclass)
-!real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
-!    +       maxpointspec_act,maxageclass)
-!real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
-!    +       maxpointspec_act,maxageclass),
-!    +     drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
-!    +     maxpointspec_act,maxageclass),
-!    +     wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
-!    +     maxpointspec_act,maxageclass)
-!real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
-!real sparse_dump_r(numxgrid*numygrid*numzgrid)
-!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
-!real sparse_dump_u(numxgrid*numygrid*numzgrid)
-
-  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
-  real,parameter :: weightair=28.97
-  logical :: sp_zer
-  character :: adate*8,atime*6
-  character(len=3) :: anspec
-  logical :: lexist
-  character :: areldate*8,areltime*6
-  logical,save :: lstart=.true.
-  logical,save,allocatable,dimension(:) :: lstartrel
-  integer :: ierr
-  character(LEN=100) :: dates_char
-  integer, parameter :: unitrelnames=654
-
-
-  if(lstart) then
-    allocate(lstartrel(maxpointspec_act))
-    lstartrel(:)=.true.
-  endif
-  print*, 'lstartrel = ',lstartrel
-
-  if (verbosity.eq.1) then
-    print*,'inside concoutput_inversion '
-    CALL SYSTEM_CLOCK(count_clock)
-    WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
-
-! Determine current calendar date
-!**********************************************************
-
-  jul=bdate+real(itime,kind=dp)/86400._dp
-  call caldate(jul,jjjjmmdd,ihmmss)
-  write(adate,'(i8.8)') jjjjmmdd
-  write(atime,'(i6.6)') ihmmss
-!  write(unitdates,'(a)') adate//atime
-
-
-! Prepare output files for dates
-!**********************************************************
-
-  ! Overwrite existing dates file on first call, later append to it
-  ! If 'dates' file exists in output directory, copy to new file dates.old
-  inquire(file=path(2)(1:length(2))//'dates', exist=lexist)
-  if (lexist.and.lstart) then
-    ! copy contents of existing dates file to dates.old
-    print*, 'warning: replacing old dates file'
-    open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
-         &access='sequential', status='old', action='read', iostat=ierr)
-    open(unit=unittmp, file=path(2)(1:length(2))//'dates.old', access='sequential', &
-         &status='replace', action='write', form='formatted', iostat=ierr)
-    do while (.true.)
-      read(unitdates, '(a)', iostat=ierr) dates_char
-      if (ierr.ne.0) exit
-      write(unit=unittmp, fmt='(a)', iostat=ierr, advance='yes') trim(dates_char)
-    end do
-    close(unit=unitdates)
-    close(unit=unittmp)
-    ! create new dates file
-    open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
-         &access='sequential', status='replace', iostat=ierr)
-    close(unit=unitdates)
-  endif
-
-  open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
-  write(unitdates,'(a)') adate//atime
-  close(unitdates)
-
-  !CGZ: Make a filename with names of releases
-  if (lstart) then
-    open(unit=unitrelnames, file=path(2)(1:length(2))//'releases_out',form='formatted', &
-         &access='sequential', status='replace', iostat=ierr)
-    close(unitrelnames)
-  endif
-
-  print*, 'after creating dates files: lstart = ',lstart
-!  print*, 'outnum:',outnum
-!  print*, 'datetime:',adate//atime
-
-
-! For forward simulations, output fields have dimension MAXSPEC,
-! for backward simulations, output fields have dimension MAXPOINT.
-! Thus, make loops either about nspec, or about numpoint
-!*****************************************************************
-
-
-  if (ldirect.eq.1) then
-    do ks=1,nspec
-      do kp=1,maxpointspec_act
-        tot_mu(ks,kp)=1
-      end do
-    end do
-  else
-    do ks=1,nspec
-      do kp=1,maxpointspec_act
-        tot_mu(ks,kp)=xmass(kp,ks)
-      end do
-    end do
-  endif
-
-
-  if (verbosity.eq.1) then
-    print*,'concoutput_inversion 2'
-    CALL SYSTEM_CLOCK(count_clock)
-    WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
-
-!*******************************************************************
-! Compute air density: sufficiently accurate to take it
-! from coarse grid at some time
-! Determine center altitude of output layer, and interpolate density
-! data to that altitude
-!*******************************************************************
-
-  do kz=1,numzgrid
-    if (kz.eq.1) then
-      halfheight=outheight(1)/2.
-    else
-      halfheight=(outheight(kz)+outheight(kz-1))/2.
-    endif
-    do kzz=2,nz
-      if ((height(kzz-1).lt.halfheight).and. &
-           (height(kzz).gt.halfheight)) goto 46
-    end do
-46  kzz=max(min(kzz,nz),2)
-    dz1=halfheight-height(kzz-1)
-    dz2=height(kzz)-halfheight
-    dz=dz1+dz2
-    do jy=0,numygrid-1
-      do ix=0,numxgrid-1
-        xl=outlon0+real(ix)*dxout
-        yl=outlat0+real(jy)*dyout
-        xl=(xl-xlon0)/dx
-        yl=(yl-ylat0)/dy
-        iix=max(min(nint(xl),nxmin1),0)
-        jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
-      end do
-    end do
-  end do
-
-  do i=1,numreceptor
-    xl=xreceptor(i)
-    yl=yreceptor(i)
-    iix=max(min(nint(xl),nxmin1),0)
-    jjy=max(min(nint(yl),nymin1),0)
-    densityoutrecept(i)=rho(iix,jjy,1,2)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
-  end do
-
-! RLT
-! conversion factor for output relative to dry air
-  factor_drygrid=densityoutgrid/densitydrygrid
-  factor_dryrecept=densityoutrecept/densitydryrecept
-
-! Output is different for forward and backward simulations
-  do kz=1,numzgrid
-    do jy=0,numygrid-1
-      do ix=0,numxgrid-1
-        if (ldirect.eq.1) then
-          factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
-        else
-          factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
-        endif
-      end do
-    end do
-  end do
-
-!*********************************************************************
-! Determine the standard deviation of the mean concentration or mixing
-! ratio (uncertainty of the output) and the dry and wet deposition
-!*********************************************************************
-
-  if (verbosity.eq.1) then
-    print*,'concoutput_inversion 3 (sd)'
-    CALL SYSTEM_CLOCK(count_clock)
-    WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
-  gridtotal=0.
-  gridsigmatotal=0.
-  gridtotalunc=0.
-  wetgridtotal=0.
-  wetgridsigmatotal=0.
-  wetgridtotalunc=0.
-  drygridtotal=0.
-  drygridsigmatotal=0.
-  drygridtotalunc=0.
-
-  do ks=1,nspec
-
-    write(anspec,'(i3.3)') ks
-
-    ! loop over releases
-    do kp=1,maxpointspec_act
-
-      print*, 'itime = ',itime
-      !print*, 'lage(1) = ',lage(1)
-      print*, 'ireleasestart(kp) = ',ireleasestart(kp)
-      print*, 'ireleaseend(kp) = ',ireleaseend(kp)
-
-      ! check itime is within release and backward trajectory length
-      if (nageclass.eq.1) then
-        if ((itime.gt.ireleaseend(kp)).or.(itime.lt.(ireleasestart(kp)-lage(1)))) then
-          go to 10
-        endif
-      endif
-
-      ! calculate date of release for filename
-      jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp    ! this is the current day
-      call caldate(jul,jjjjmmdd,ihmmss)     
-      write(areldate,'(i8.8)') jjjjmmdd
-      write(areltime,'(i6.6)') ihmmss
-      print*, 'areldate/areltime = ',areldate//areltime
-
-      ! calculate date of field
-      jul=bdate+real(itime,kind=dp)/86400._dp
-      call caldate(jul,jjjjmmdd,ihmmss)
-      write(adate,'(i8.8)') jjjjmmdd
-      write(atime,'(i6.6)') ihmmss
-!      print*, adate//atime
-
-      if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
-        if (ldirect.eq.1) then
-          ! concentrations
-          inquire(file=path(2)(1:length(2))//'grid_conc_'//areldate// &
-                  areltime//'_'//anspec,exist=lexist)
-          if(lexist.and..not.lstartrel(kp)) then
-            ! open and append to existing file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
-          else
-            ! open new file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='replace',action='write')
-          endif
-        else
-          ! residence times
-          inquire(file=path(2)(1:length(2))//'grid_time_'//areldate// &
-                  areltime//'_'//anspec,exist=lexist)
-          if(lexist.and..not.lstartrel(kp)) then
-            ! open and append to existing file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
-          else
-            ! open new file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='replace',action='write')
-            ! add part of the filename to a file (similar to dates) for easier post-processing
-            open(unit=unitrelnames, file=path(2)(1:length(2))//'releases_out',form='formatted', &
-                 & access='APPEND', iostat=ierr)
-            write(unitrelnames,'(a)') areldate//areltime//'_'//anspec
-            call flush(unit=unitrelnames)
-            close(unitrelnames)
-          endif
-        endif
-        write(unitoutgrid) jjjjmmdd
-        write(unitoutgrid) ihmmss
-      endif
-
-      if ((iout.eq.2).or.(iout.eq.3)) then      
-        ! mixing ratio
-        inquire(file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
-                areltime//'_'//anspec,exist=lexist)
-        if(lexist.and..not.lstartrel(kp)) then
-          ! open and append to existing file
-          open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
-               areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
-        else
-          ! open new file
-          open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
-               areltime//'_'//anspec,form='unformatted',status='replace',action='write')
-        endif   
-        write(unitoutgridppt) jjjjmmdd
-        write(unitoutgridppt) ihmmss
-      endif
-
-      lstartrel(kp)=.false.
-
-      do nage=1,nageclass
-
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-
-!! WET DEPOSITION
-!            if ((WETDEP).and.(ldirect.gt.0)) then
-!              do l=1,nclassunc
-!                auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
-!              end do
-!              call mean(auxgrid,wetgrid(ix,jy), &
-!                   wetgridsigma(ix,jy),nclassunc)
-!! Multiply by number of classes to get total concentration
-!              wetgrid(ix,jy)=wetgrid(ix,jy) &
-!                   *nclassunc
-!              wetgridtotal=wetgridtotal+wetgrid(ix,jy)
-!! Calculate standard deviation of the mean
-!              wetgridsigma(ix,jy)= &
-!                   wetgridsigma(ix,jy)* &
-!                   sqrt(real(nclassunc))
-!              wetgridsigmatotal=wetgridsigmatotal+ &
-!                   wetgridsigma(ix,jy)
-!            endif
-
-!! DRY DEPOSITION
-!            if ((DRYDEP).and.(ldirect.gt.0)) then
-!              do l=1,nclassunc
-!                auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
-!              end do
-!              call mean(auxgrid,drygrid(ix,jy), &
-!                   drygridsigma(ix,jy),nclassunc)
-!! Multiply by number of classes to get total concentration
-!              drygrid(ix,jy)=drygrid(ix,jy)* &
-!                   nclassunc
-!              drygridtotal=drygridtotal+drygrid(ix,jy)
-!! Calculate standard deviation of the mean
-!              drygridsigma(ix,jy)= &
-!                   drygridsigma(ix,jy)* &
-!                   sqrt(real(nclassunc))
-!125           drygridsigmatotal=drygridsigmatotal+ &
-!                   drygridsigma(ix,jy)
-!            endif
-
-! CONCENTRATION OR MIXING RATIO
-            do kz=1,numzgrid
-              do l=1,nclassunc
-                auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
-              end do
-              call mean(auxgrid,grid(ix,jy,kz), &
-                   gridsigma(ix,jy,kz),nclassunc)
-! Multiply by number of classes to get total concentration
-              grid(ix,jy,kz)= &
-                   grid(ix,jy,kz)*nclassunc
-              gridtotal=gridtotal+grid(ix,jy,kz)
-! Calculate standard deviation of the mean
-              gridsigma(ix,jy,kz)= &
-                   gridsigma(ix,jy,kz)* &
-                   sqrt(real(nclassunc))
-              gridsigmatotal=gridsigmatotal+ &
-                   gridsigma(ix,jy,kz)
-            end do
-          end do
-        end do
-
-
-!*******************************************************************
-! Generate output: may be in concentration (ng/m3) or in mixing
-! ratio (ppt) or both
-! Output the position and the values alternated multiplied by
-! 1 or -1, first line is number of values, number of positions
-! For backward simulations, the unit is seconds, stored in grid_time
-!*******************************************************************
-
-        if (verbosity.eq.1) then
-          print*,'concoutput_inversion 4 (output)'
-          CALL SYSTEM_CLOCK(count_clock)
-          WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-        endif
-
-! Concentration output
-!*********************
-
-        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
-
-          if (verbosity.eq.1) then
-            print*,'concoutput_inversion (Wet deposition)'
-            CALL SYSTEM_CLOCK(count_clock)
-            WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-          endif
-
-! Wet deposition
-!          sp_count_i=0
-!          sp_count_r=0
-!          sp_fact=-1.
-!          sp_zer=.true.
-!          if ((ldirect.eq.1).and.(WETDEP)) then
-!            do jy=0,numygrid-1
-!              do ix=0,numxgrid-1
-!! concentration greater zero
-!                if (wetgrid(ix,jy).gt.smallnum) then
-!                  if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)=ix+jy*numxgrid
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1.)
-!                  endif
-!                  sp_count_r=sp_count_r+1
-!                  sparse_dump_r(sp_count_r)= &
-!                       sp_fact*1.e12*wetgrid(ix,jy)/area(ix,jy)
-!                  sparse_dump_u(sp_count_r)= &
-!                       1.e12*wetgridsigma(ix,jy)/area(ix,jy)
-!                else ! concentration is zero
-!                  sp_zer=.true.
-!                endif
-!              end do
-!            end do
-!          else
-!            sp_count_i=0
-!            sp_count_r=0
-!          endif
-!          write(unitoutgrid) sp_count_i
-!          write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
-!          write(unitoutgrid) sp_count_r
-!          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
-!
-!          if (verbosity.eq.1) then
-!            print*,'concoutput_surf (Dry deposition)'
-!            CALL SYSTEM_CLOCK(count_clock)
-!            WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-!          endif
-!! Dry deposition
-!          sp_count_i=0
-!          sp_count_r=0
-!          sp_fact=-1.
-!          sp_zer=.true.
-!          if ((ldirect.eq.1).and.(DRYDEP)) then
-!            do jy=0,numygrid-1
-!              do ix=0,numxgrid-1
-!                if (drygrid(ix,jy).gt.smallnum) then
-!                  if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)=ix+jy*numxgrid
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1.)
-!                  endif
-!                  sp_count_r=sp_count_r+1
-!                  sparse_dump_r(sp_count_r)= &
-!                       sp_fact* &
-!                       1.e12*drygrid(ix,jy)/area(ix,jy)
-!                  sparse_dump_u(sp_count_r)= &
-!                       1.e12*drygridsigma(ix,jy)/area(ix,jy)
-!                else ! concentration is zero
-!                  sp_zer=.true.
-!                endif
-!              end do
-!            end do
-!          else
-!            sp_count_i=0
-!            sp_count_r=0
-!          endif
-!          write(unitoutgrid) sp_count_i
-!          write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
-!          write(unitoutgrid) sp_count_r
-!          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
-
-          if (verbosity.eq.1) then
-            print*,'concoutput_inversion (Concentrations)'
-            CALL SYSTEM_CLOCK(count_clock)
-            WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-          endif
-
-! Concentrations
-
-! surf_only write only 1st layer 
-
-          sp_count_i=0
-          sp_count_r=0
-          sp_fact=-1.
-          sp_zer=.true.
-          do kz=1,1
-            do jy=0,numygrid-1
-              do ix=0,numxgrid-1
-                if (grid(ix,jy,kz).gt.smallnum) then
-                  if (sp_zer.eqv..true.) then ! first non zero value
-                    sp_count_i=sp_count_i+1
-                    sparse_dump_i(sp_count_i)= &
-                         ix+jy*numxgrid+kz*numxgrid*numygrid
-                    sp_zer=.false.
-                    sp_fact=sp_fact*(-1.)
-                  endif
-                  sp_count_r=sp_count_r+1
-                  sparse_dump_r(sp_count_r)= &
-                       sp_fact* &
-                       grid(ix,jy,kz)* &
-                       factor3d(ix,jy,kz)/tot_mu(ks,kp)
-!                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
-!    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
-                  sparse_dump_u(sp_count_r)= &
-                       gridsigma(ix,jy,kz)* &
-                       factor3d(ix,jy,kz)/tot_mu(ks,kp)
-
-                else ! concentration is zero
-                  sp_zer=.true.
-                endif
-              end do
-            end do
-          end do
-          write(unitoutgrid) sp_count_i
-          write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
-          write(unitoutgrid) sp_count_r
-          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!          write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
-
-        endif !  concentration output
-
-! Mixing ratio output
-!********************
-
-        if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
-!
-!! Wet deposition
-!          sp_count_i=0
-!          sp_count_r=0
-!          sp_fact=-1.
-!          sp_zer=.true.
-!          if ((ldirect.eq.1).and.(WETDEP)) then
-!            do jy=0,numygrid-1
-!              do ix=0,numxgrid-1
-!                if (wetgrid(ix,jy).gt.smallnum) then
-!                  if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)= &
-!                         ix+jy*numxgrid
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1.)
-!                  endif
-!                  sp_count_r=sp_count_r+1
-!                  sparse_dump_r(sp_count_r)= &
-!                       sp_fact* &
-!                       1.e12*wetgrid(ix,jy)/area(ix,jy)
-!                  sparse_dump_u(sp_count_r)= &
-!                       1.e12*wetgridsigma(ix,jy)/area(ix,jy)
-!                else ! concentration is zero
-!                  sp_zer=.true.
-!                endif
-!              end do
-!            end do
-!          else
-!            sp_count_i=0
-!            sp_count_r=0
-!          endif
-!          write(unitoutgridppt) sp_count_i
-!          write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
-!          write(unitoutgridppt) sp_count_r
-!          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
-!
-
-! Dry deposition
-!          sp_count_i=0
-!          sp_count_r=0
-!          sp_fact=-1.
-!          sp_zer=.true.
-!          if ((ldirect.eq.1).and.(DRYDEP)) then
-!            do jy=0,numygrid-1
-!              do ix=0,numxgrid-1
-!                if (drygrid(ix,jy).gt.smallnum) then
-!                  if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)= &
-!                         ix+jy*numxgrid
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1)
-!                  endif
-!                  sp_count_r=sp_count_r+1
-!                  sparse_dump_r(sp_count_r)= &
-!                       sp_fact* &
-!                       1.e12*drygrid(ix,jy)/area(ix,jy)
-!                  sparse_dump_u(sp_count_r)= &
-!                       1.e12*drygridsigma(ix,jy)/area(ix,jy)
-!                else ! concentration is zero
-!                  sp_zer=.true.
-!                endif
-!              end do
-!            end do
-!          else
-!            sp_count_i=0
-!            sp_count_r=0
-!          endif
-!          write(unitoutgridppt) sp_count_i
-!          write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
-!          write(unitoutgridppt) sp_count_r
-!          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
-!
-
-! Mixing ratios
-
-! surf_only write only 1st layer 
-
-          sp_count_i=0
-          sp_count_r=0
-          sp_fact=-1.
-          sp_zer=.true.
-          do kz=1,1
-            do jy=0,numygrid-1
-              do ix=0,numxgrid-1
-                if (grid(ix,jy,kz).gt.smallnum) then
-                  if (sp_zer.eqv..true.) then ! first non zero value
-                    sp_count_i=sp_count_i+1
-                    sparse_dump_i(sp_count_i)= &
-                         ix+jy*numxgrid+kz*numxgrid*numygrid
-                    sp_zer=.false.
-                    sp_fact=sp_fact*(-1.)
-                  endif
-                  sp_count_r=sp_count_r+1
-                  sparse_dump_r(sp_count_r)= &
-                       sp_fact* &
-                       1.e12*grid(ix,jy,kz) &
-                       /volume(ix,jy,kz)/outnum* &
-                       weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
-                  sparse_dump_u(sp_count_r)= &
-                       1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ &
-                       outnum*weightair/weightmolar(ks)/ &
-                       densityoutgrid(ix,jy,kz)
-                else ! concentration is zero
-                  sp_zer=.true.
-                endif
-              end do
-            end do
-          end do
-          write(unitoutgridppt) sp_count_i
-          write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
-          write(unitoutgridppt) sp_count_r
-          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!          write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
-
-        endif ! output for ppt
-
-      end do  ! nageclass
-
-      close(unitoutgridppt)
-      close(unitoutgrid)
-
-      ! itime is outside range
-10    continue
-
-    end do  ! maxpointspec_act
-
-  end do  ! nspec
-
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
-  if (lexist.and..not.lstart) then
-    ! open and append
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
-            status='old',action='write',access='append')
-  else
-    ! create new
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
-            status='replace',action='write')
-  endif
-  sp_count_i=0
-  sp_count_r=0
-  sp_fact=-1.
-  sp_zer=.true.
-  do kz=1,1
-    do jy=0,numygrid-1
-      do ix=0,numxgrid-1
-        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
-          if (sp_zer.eqv..true.) then ! first value not equal to one
-            sp_count_i=sp_count_i+1
-            sparse_dump_i(sp_count_i)= &
-                  ix+jy*numxgrid+kz*numxgrid*numygrid
-            sp_zer=.false.
-            sp_fact=sp_fact*(-1.)
-          endif
-          sp_count_r=sp_count_r+1
-          sparse_dump_r(sp_count_r)= &
-               sp_fact*factor_drygrid(ix,jy,kz)
-        else ! factor is one
-          sp_zer=.true.
-        endif
-      end do
-    end do
-  end do
-  write(unitoutfactor) sp_count_i
-  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
-  write(unitoutfactor) sp_count_r
-  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
-  close(unitoutfactor)
-
-
-  if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
-!  if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
-!       wetgridtotal
-!  if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
-!       drygridtotal
-
-! Dump of receptor concentrations
-
-  if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3)  ) then
-    write(unitoutreceptppt) itime
-    do ks=1,nspec
-      write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* &
-           weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
-    end do
-  endif
-
-! Dump of receptor concentrations
-
-  if (numreceptor.gt.0) then
-    write(unitoutrecept) itime
-    do ks=1,nspec
-      write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &
-           i=1,numreceptor)
-    end do
-  endif
-
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  if (numreceptor.gt.0) then
-    inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
-     if (lexist.and..not.lstart) then
-     ! open and append
-      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
-              status='old',action='write',access='append')
-    else
-      ! create new
-      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
-              status='replace',action='write')
-    endif
-    write(unitoutfactor) itime
-    write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
-    close(unitoutfactor)
-  endif
-
-  ! reset lstart
-  if (lstart) then
-    lstart=.false.
-  endif
-  print*, 'after writing output files: lstart = ',lstart
-
-
-! Reinitialization of grid
-!*************************
-
-  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
-
-end subroutine concoutput_inversion
Index: c/concoutput_inversion_nest.f90
===================================================================
--- src/concoutput_inversion_nest.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ 	(revision )
@@ -1,699 +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 concoutput_inversion_nest(itime,outnum)
-  !                        i     i
-  !*****************************************************************************
-  !                                                                            *
-  !     Output of the concentration grid and the receptor concentrations.      *
-  !                                                                            *
-  !     Author: A. Stohl                                                       *
-  !                                                                            *
-  !     24 May 1995                                                            *
-  !                                                                            *
-  !     13 April 1999, Major update: if output size is smaller, dump output    *
-  !                    in sparse matrix format; additional output of           *
-  !                    uncertainty                                             *
-  !                                                                            *
-  !     05 April 2000, Major update: output of age classes; output for backward*
-  !                    runs is time spent in grid cell times total mass of     *
-  !                    species.                                                *
-  !                                                                            *
-  !     17 February 2002, Appropriate dimensions for backward and forward runs *
-  !                       are now specified in file par_mod                    *
-  !                                                                            *
-  !     June 2006, write grid in sparse matrix with a single write command     *
-  !                in order to save disk space                                 *
-  !                                                                            *
-  !     2008 new sparse matrix format                                          *
-  !
-  !     January 2017,  Separate files by release but include all timesteps     *
-  !                                                                            *
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  ! outnum          number of samples                                          *
-  ! ncells          number of cells with non-zero concentrations               *
-  ! sparse          .true. if in sparse matrix format, else .false.            *
-  ! tot_mu          1 for forward, initial mass mixing ration for backw. runs  *
-  !                                                                            *
-  !*****************************************************************************
-
-  use unc_mod
-  use point_mod
-  use outg_mod
-  use par_mod
-  use com_mod
-  use mean_mod
-
-  implicit none
-
-  real(kind=dp) :: jul
-  integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
-  integer :: sp_count_i,sp_count_r
-  real :: sp_fact
-  real :: outnum,densityoutrecept(maxreceptor),xl,yl
-! RLT
-  real :: densitydryrecept(maxreceptor)
-  real :: factor_dryrecept(maxreceptor)
-
-  !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
-  !    +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
-  !    +    maxageclass)
-  !real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
-  !    +       maxageclass)
-  !real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
-  !    +       maxpointspec_act,maxageclass)
-  !real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
-  !    +       maxpointspec_act,maxageclass),
-  !    +     drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
-  !    +     maxpointspec_act,maxageclass),
-  !    +     wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
-  !    +     maxpointspec_act,maxageclass)
-  !real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
-  !real sparse_dump_r(numxgrid*numygrid*numzgrid)
-  !integer sparse_dump_i(numxgrid*numygrid*numzgrid)
-
-  !real sparse_dump_u(numxgrid*numygrid*numzgrid)
-  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
-  real,parameter :: weightair=28.97
-  logical :: sp_zer
-  logical,save :: lnstart=.true.
-  logical,save,allocatable,dimension(:) :: lnstartrel
-  character :: adate*8,atime*6
-  character(len=3) :: anspec
-  logical :: lexist
-  character :: areldate*8,areltime*6
-
-  if(lnstart) then
-    allocate(lnstartrel(maxpointspec_act))
-    lnstartrel(:)=.true.
-  endif
-  print*, 'lnstartrel = ',lnstartrel
-
-  ! 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
-
-  print*, 'outnum:',outnum
-  print*, 'datetime:',adate//atime
-
-  ! For forward simulations, output fields have dimension MAXSPEC,
-  ! for backward simulations, output fields have dimension MAXPOINT.
-  ! Thus, make loops either about nspec, or about numpoint
-  !*****************************************************************
-
-
-    if (ldirect.eq.1) then
-       do ks=1,nspec
-         do kp=1,maxpointspec_act
-           tot_mu(ks,kp)=1
-         end do
-       end do
-   else
-      do ks=1,nspec
-             do kp=1,maxpointspec_act
-               tot_mu(ks,kp)=xmass(kp,ks)
-             end do
-      end do
-    endif
-
-
-  !*******************************************************************
-  ! Compute air density: sufficiently accurate to take it
-  ! from coarse grid at some time
-  ! Determine center altitude of output layer, and interpolate density
-  ! data to that altitude
-  !*******************************************************************
-
-  do kz=1,numzgrid
-    if (kz.eq.1) then
-      halfheight=outheight(1)/2.
-    else
-      halfheight=(outheight(kz)+outheight(kz-1))/2.
-    endif
-    do kzz=2,nz
-      if ((height(kzz-1).lt.halfheight).and. &
-           (height(kzz).gt.halfheight)) goto 46
-    end do
-46   kzz=max(min(kzz,nz),2)
-    dz1=halfheight-height(kzz-1)
-    dz2=height(kzz)-halfheight
-    dz=dz1+dz2
-    do jy=0,numygridn-1
-      do ix=0,numxgridn-1
-        xl=outlon0n+real(ix)*dxoutn
-        yl=outlat0n+real(jy)*dyoutn
-        xl=(xl-xlon0)/dx
-        yl=(yl-ylat0)/dy
-        iix=max(min(nint(xl),nxmin1),0)
-        jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
-      end do
-    end do
-  end do
-
-  do i=1,numreceptor
-    xl=xreceptor(i)
-    yl=yreceptor(i)
-    iix=max(min(nint(xl),nxmin1),0)
-    jjy=max(min(nint(yl),nymin1),0)
-    densityoutrecept(i)=rho(iix,jjy,1,2)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
-  end do
-
-! RLT
-! conversion factor for output relative to dry air
-  factor_drygrid=densityoutgrid/densitydrygrid
-  factor_dryrecept=densityoutrecept/densitydryrecept
-
-  ! Output is different for forward and backward simulations
-    do kz=1,numzgrid
-      do jy=0,numygridn-1
-        do ix=0,numxgridn-1
-          if (ldirect.eq.1) then
-            factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
-          else
-            factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
-          endif
-        end do
-      end do
-    end do
-
-  !*********************************************************************
-  ! Determine the standard deviation of the mean concentration or mixing
-  ! ratio (uncertainty of the output) and the dry and wet deposition
-  !*********************************************************************
-
-  do ks=1,nspec
-
-  write(anspec,'(i3.3)') ks
-
-    do kp=1,maxpointspec_act
-
-      print*, 'itime = ',itime
-      print*, 'lage(1) = ',lage(1)
-      print*, 'ireleasestart(kp) = ',ireleasestart(kp)
-      print*, 'ireleaseend(kp) = ',ireleaseend(kp)
-
-      ! check itime is within release and backward trajectory length
-      if (nageclass.eq.1) then
-        if ((itime.gt.ireleaseend(kp)).or.(itime.lt.(ireleasestart(kp)-lage(1)))) then
-          go to 10
-        endif
-      endif
-
-      ! calculate date of release
-      jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp    ! this is the current day
-      call caldate(jul,jjjjmmdd,ihmmss)
-      write(areldate,'(i8.8)') jjjjmmdd
-      write(areltime,'(i6.6)') ihmmss
-      print*, areldate//areltime
-
-      ! calculate date of field
-      jul=bdate+real(itime,kind=dp)/86400._dp
-      call caldate(jul,jjjjmmdd,ihmmss)
-      write(adate,'(i8.8)') jjjjmmdd
-      write(atime,'(i6.6)') ihmmss
-      print*, adate//atime
-
-      if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
-        if (ldirect.eq.1) then
-          ! concentrations
-          inquire(file=path(2)(1:length(2))//'grid_conc_nest_'//areldate// &
-                  areltime//'_'//anspec,exist=lexist)
-          if(lexist.and..not.lnstartrel(kp)) then
-            ! open and append to existing file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
-          else
-            ! open new file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='replace',action='write')
-          endif
-        else
-          ! residence times
-          inquire(file=path(2)(1:length(2))//'grid_time_nest_'//areldate// &
-                  areltime//'_'//anspec,exist=lexist)
-          if(lexist.and..not.lnstartrel(kp)) then
-            ! open and append to existing file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
-          else
-            ! open new file
-            open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_'//areldate// &
-                 areltime//'_'//anspec,form='unformatted',status='replace',action='write')
-          endif
-        endif
-        write(unitoutgrid) jjjjmmdd
-        write(unitoutgrid) ihmmss
-      endif
-
-      if ((iout.eq.2).or.(iout.eq.3)) then
-        ! mixing ratio
-        inquire(file=path(2)(1:length(2))//'grid_pptv_nest_'//areldate// &
-                areltime//'_'//anspec,exist=lexist)
-        if(lexist.and..not.lnstartrel(kp)) then
-          ! open and append to existing file
-          open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_'//areldate// &
-               areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
-        else
-          ! open new file
-          open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_'//areldate// &
-               areltime//'_'//anspec,form='unformatted',status='replace',action='write')
-        endif
-        write(unitoutgridppt) jjjjmmdd
-        write(unitoutgridppt) ihmmss
-      endif
-
-      lnstartrel(kp)=.false.
-
-      do nage=1,nageclass
-
-        do jy=0,numygridn-1
-          do ix=0,numxgridn-1
-
-!  ! WET DEPOSITION
-!            if ((WETDEP).and.(ldirect.gt.0)) then
-!              do l=1,nclassunc
-!                auxgrid(l)=wetgriduncn(ix,jy,ks,kp,l,nage)
-!              end do
-!              call mean(auxgrid,wetgrid(ix,jy), &
-!                   wetgridsigma(ix,jy),nclassunc)
-!  ! Multiply by number of classes to get total concentration
-!              wetgrid(ix,jy)=wetgrid(ix,jy) &
-!                   *nclassunc
-!  ! Calculate standard deviation of the mean
-!              wetgridsigma(ix,jy)= &
-!                   wetgridsigma(ix,jy)* &
-!                   sqrt(real(nclassunc))
-!            endif
-
-!  ! DRY DEPOSITION
-!            if ((DRYDEP).and.(ldirect.gt.0)) then
-!              do l=1,nclassunc
-!                auxgrid(l)=drygriduncn(ix,jy,ks,kp,l,nage)
-!              end do
-!              call mean(auxgrid,drygrid(ix,jy), &
-!                   drygridsigma(ix,jy),nclassunc)
-!  ! Multiply by number of classes to get total concentration
-!              drygrid(ix,jy)=drygrid(ix,jy)* &
-!                   nclassunc
-!  ! Calculate standard deviation of the mean
-!              drygridsigma(ix,jy)= &
-!                   drygridsigma(ix,jy)* &
-!                   sqrt(real(nclassunc))
-!            endif
-
-  ! CONCENTRATION OR MIXING RATIO
-            do kz=1,numzgrid
-              do l=1,nclassunc
-                auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
-              end do
-              call mean(auxgrid,grid(ix,jy,kz), &
-                     gridsigma(ix,jy,kz),nclassunc)
-  ! Multiply by number of classes to get total concentration
-              grid(ix,jy,kz)= &
-                   grid(ix,jy,kz)*nclassunc
-  ! Calculate standard deviation of the mean
-              gridsigma(ix,jy,kz)= &
-                   gridsigma(ix,jy,kz)* &
-                   sqrt(real(nclassunc))
-            end do
-          end do
-        end do
-
-
-  !*******************************************************************
-  ! Generate output: may be in concentration (ng/m3) or in mixing
-  ! ratio (ppt) or both
-  ! Output the position and the values alternated multiplied by
-  ! 1 or -1, first line is number of values, number of positions
-  ! For backward simulations, the unit is seconds, stored in grid_time
-  !*******************************************************************
-
-  ! Concentration output
-  !*********************
-
-        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
-
-!  ! Wet deposition
-!          sp_count_i=0
-!          sp_count_r=0
-!          sp_fact=-1.
-!          sp_zer=.true.
-!          if ((ldirect.eq.1).and.(WETDEP)) then
-!          do jy=0,numygridn-1
-!            do ix=0,numxgridn-1
-!  ! concentration greater zero
-!              if (wetgrid(ix,jy).gt.smallnum) then
-!                 if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)=ix+jy*numxgridn
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1.)
-!                 endif
-!                 sp_count_r=sp_count_r+1
-!                 sparse_dump_r(sp_count_r)= &
-!                      sp_fact*1.e12*wetgrid(ix,jy)/arean(ix,jy)
-!                 sparse_dump_u(sp_count_r)= &
-!                      1.e12*wetgridsigma(ix,jy)/area(ix,jy)
-!              else ! concentration is zero
-!                  sp_zer=.true.
-!              endif
-!            end do
-!         end do
-!         else
-!            sp_count_i=0
-!            sp_count_r=0
-!         endif
-!         write(unitoutgrid) sp_count_i
-!         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
-
-!  ! Dry deposition
-!         sp_count_i=0
-!         sp_count_r=0
-!         sp_fact=-1.
-!         sp_zer=.true.
-!         if ((ldirect.eq.1).and.(DRYDEP)) then
-!          do jy=0,numygridn-1
-!            do ix=0,numxgridn-1
-!              if (drygrid(ix,jy).gt.smallnum) then
-!                 if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)=ix+jy*numxgridn
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1.)
-!                 endif
-!                 sp_count_r=sp_count_r+1
-!                 sparse_dump_r(sp_count_r)= &
-!                      sp_fact* &
-!                      1.e12*drygrid(ix,jy)/arean(ix,jy)
-!                 sparse_dump_u(sp_count_r)= &
-!                      1.e12*drygridsigma(ix,jy)/area(ix,jy)
-!              else ! concentration is zero
-!                  sp_zer=.true.
-!              endif
-!            end do
-!          end do
-!         else
-!            sp_count_i=0
-!            sp_count_r=0
-!         endif
-!         write(unitoutgrid) sp_count_i
-!         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
-!
-
-  ! Concentrations
-
-  ! surf_only write only 1st layer 
-
-         sp_count_i=0
-         sp_count_r=0
-         sp_fact=-1.
-         sp_zer=.true.
-          do kz=1,1
-            do jy=0,numygridn-1
-              do ix=0,numxgridn-1
-                if (grid(ix,jy,kz).gt.smallnum) then
-                  if (sp_zer.eqv..true.) then ! first non zero value
-                    sp_count_i=sp_count_i+1
-                    sparse_dump_i(sp_count_i)= &
-                         ix+jy*numxgridn+kz*numxgridn*numygridn
-                    sp_zer=.false.
-                    sp_fact=sp_fact*(-1.)
-                   endif
-                   sp_count_r=sp_count_r+1
-                   sparse_dump_r(sp_count_r)= &
-                        sp_fact* &
-                        grid(ix,jy,kz)* &
-                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
-  !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
-  !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
-                   sparse_dump_u(sp_count_r)= &
-                        gridsigma(ix,jy,kz)* &
-                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
-              else ! concentration is zero
-                  sp_zer=.true.
-              endif
-              end do
-            end do
-          end do
-         write(unitoutgrid) sp_count_i
-         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
-         write(unitoutgrid) sp_count_r
-         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
-
-      endif !  concentration output
-
-  ! Mixing ratio output
-  !********************
-
-      if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
-
-!  ! Wet deposition
-!         sp_count_i=0
-!         sp_count_r=0
-!         sp_fact=-1.
-!         sp_zer=.true.
-!         if ((ldirect.eq.1).and.(WETDEP)) then
-!          do jy=0,numygridn-1
-!            do ix=0,numxgridn-1
-!                if (wetgrid(ix,jy).gt.smallnum) then
-!                  if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)= &
-!                         ix+jy*numxgridn
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1.)
-!                 endif
-!                 sp_count_r=sp_count_r+1
-!                 sparse_dump_r(sp_count_r)= &
-!                      sp_fact* &
-!                      1.e12*wetgrid(ix,jy)/arean(ix,jy)
-!                 sparse_dump_u(sp_count_r)= &
-!                      1.e12*wetgridsigma(ix,jy)/area(ix,jy)
-!              else ! concentration is zero
-!                  sp_zer=.true.
-!              endif
-!            end do
-!          end do
-!         else
-!           sp_count_i=0
-!           sp_count_r=0
-!         endif
-!         write(unitoutgridppt) sp_count_i
-!         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
-!
-
-!  ! Dry deposition
-!         sp_count_i=0
-!         sp_count_r=0
-!         sp_fact=-1.
-!         sp_zer=.true.
-!         if ((ldirect.eq.1).and.(DRYDEP)) then
-!          do jy=0,numygridn-1
-!            do ix=0,numxgridn-1
-!                if (drygrid(ix,jy).gt.smallnum) then
-!                  if (sp_zer.eqv..true.) then ! first non zero value
-!                    sp_count_i=sp_count_i+1
-!                    sparse_dump_i(sp_count_i)= &
-!                         ix+jy*numxgridn
-!                    sp_zer=.false.
-!                    sp_fact=sp_fact*(-1)
-!                 endif
-!                 sp_count_r=sp_count_r+1
-!                 sparse_dump_r(sp_count_r)= &
-!                      sp_fact* &
-!                      1.e12*drygrid(ix,jy)/arean(ix,jy)
-!                 sparse_dump_u(sp_count_r)= &
-!                      1.e12*drygridsigma(ix,jy)/area(ix,jy)
-!              else ! concentration is zero
-!                  sp_zer=.true.
-!              endif
-!            end do
-!          end do
-!         else
-!           sp_count_i=0
-!           sp_count_r=0
-!         endif
-!         write(unitoutgridppt) sp_count_i
-!         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
-!
-
-  ! Mixing ratios
-
-    ! surf_only write only 1st layer 
-
-         sp_count_i=0
-         sp_count_r=0
-         sp_fact=-1.
-         sp_zer=.true.
-          do kz=1,1
-            do jy=0,numygridn-1
-              do ix=0,numxgridn-1
-                if (grid(ix,jy,kz).gt.smallnum) then
-                  if (sp_zer.eqv..true.) then ! first non zero value
-                    sp_count_i=sp_count_i+1
-                    sparse_dump_i(sp_count_i)= &
-                         ix+jy*numxgridn+kz*numxgridn*numygridn
-                    sp_zer=.false.
-                    sp_fact=sp_fact*(-1.)
-                 endif
-                 sp_count_r=sp_count_r+1
-                 sparse_dump_r(sp_count_r)= &
-                      sp_fact* &
-                      1.e12*grid(ix,jy,kz) &
-                      /volumen(ix,jy,kz)/outnum* &
-                      weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
-                 sparse_dump_u(sp_count_r)= &
-                      1.e12*gridsigma(ix,jy,kz)/volumen(ix,jy,kz)/ &
-                      outnum*weightair/weightmolar(ks)/ &
-                      densityoutgrid(ix,jy,kz)
-              else ! concentration is zero
-                  sp_zer=.true.
-              endif
-              end do
-            end do
-          end do
-          write(unitoutgridppt) sp_count_i
-          write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
-          write(unitoutgridppt) sp_count_r
-          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!          write(unitoutgridppt) sp_count_r
-!          write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
-
-        endif ! output for ppt
- 
-      end do ! nageclass
-
-      close(unitoutgridppt)
-      close(unitoutgrid)
-
-      ! itime is outside range
-10    continue
-
-    end do ! maxpointspec_act
-
-  end do ! nspec
-
-
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
-  if (lexist.and..not.lnstart) then
-    ! open and append
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
-            status='old',action='write',access='append')
-  else
-    ! create new
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
-            status='replace',action='write')
-  endif
-  sp_count_i=0
-  sp_count_r=0
-  sp_fact=-1.
-  sp_zer=.true.
-  do kz=1,1
-    do jy=0,numygridn-1
-      do ix=0,numxgridn-1
-        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
-          if (sp_zer.eqv..true.) then ! first value not equal to one
-            sp_count_i=sp_count_i+1
-            sparse_dump_i(sp_count_i)= &
-                  ix+jy*numxgridn+kz*numxgridn*numygridn
-            sp_zer=.false.
-            sp_fact=sp_fact*(-1.)
-          endif
-          sp_count_r=sp_count_r+1
-          sparse_dump_r(sp_count_r)= &
-               sp_fact*factor_drygrid(ix,jy,kz)
-        else ! factor is one
-          sp_zer=.true.
-        endif
-      end do
-    end do
-  end do
-  write(unitoutfactor) sp_count_i
-  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
-  write(unitoutfactor) sp_count_r
-  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
-  close(unitoutfactor)
-
-  ! reset lnstart
-  if (lnstart) then
-    lnstart=.false.
-  endif
-
-  ! Reinitialization of grid
-  !*************************
-
-  do ks=1,nspec
-  do kp=1,maxpointspec_act
-    do i=1,numreceptor
-      creceptor(i,ks)=0.
-    end do
-    do jy=0,numygridn-1
-      do ix=0,numxgridn-1
-        do l=1,nclassunc
-          do nage=1,nageclass
-            do kz=1,numzgrid
-              griduncn(ix,jy,kz,ks,kp,l,nage)=0.
-            end do
-          end do
-        end do
-      end do
-    end do
-  end do
-  end do
-
-
-end subroutine concoutput_inversion_nest
-
Index: src/concoutput_mpi.f90
===================================================================
--- src/concoutput_mpi.f90	(revision 6741557b2e77046144879c12ca45ea39587f7099)
+++ src/concoutput_mpi.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -364,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
Index: src/concoutput_nest.f90
===================================================================
--- src/concoutput_nest.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/concoutput_nest.f90	(revision 7d02c2f5693a56e1b4e9b3923f9343c0a5e6b1fa)
@@ -70,7 +70,4 @@
   real :: sp_fact
   real :: outnum,densityoutrecept(maxreceptor),xl,yl
-! RLT
-  real :: densitydryrecept(maxreceptor)
-  real :: factor_dryrecept(maxreceptor)
 
   !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
@@ -99,5 +96,4 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
-  logical :: lexist
   integer :: mind 
 ! mind        eso:added to ensure identical results between 2&3-fields versions
@@ -168,7 +164,4 @@
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
              rho(iix,jjy,kzz-1,mind)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
@@ -182,12 +175,6 @@
     !densityoutrecept(i)=rho(iix,jjy,1,2)
     densityoutrecept(i)=rho(iix,jjy,1,mind)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
   end do
 
-! RLT
-! conversion factor for output relative to dry air
-  factor_drygrid=densityoutgrid/densitydrygrid
-  factor_dryrecept=densityoutrecept/densitydryrecept
 
   ! Output is different for forward and backward simulations
@@ -564,45 +551,4 @@
   end do
 
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
-  if (lexist) then
-    ! open and append
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
-            status='old',action='write',access='append')
-  else
-    ! create new
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
-            status='new',action='write')
-  endif
-  sp_count_i=0
-  sp_count_r=0
-  sp_fact=-1.
-  sp_zer=.true.
-  do kz=1,numzgrid
-    do jy=0,numygridn-1
-      do ix=0,numxgridn-1
-        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
-          if (sp_zer.eqv..true.) then ! first value not equal to one
-            sp_count_i=sp_count_i+1
-            sparse_dump_i(sp_count_i)= &
-                  ix+jy*numxgridn+kz*numxgridn*numygridn
-            sp_zer=.false.
-            sp_fact=sp_fact*(-1.)
-          endif
-          sp_count_r=sp_count_r+1
-          sparse_dump_r(sp_count_r)= &
-               sp_fact*factor_drygrid(ix,jy,kz)
-        else ! factor is one
-          sp_zer=.true.
-        endif
-      end do
-    end do
-  end do
-  write(unitoutfactor) sp_count_i
-  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
-  write(unitoutfactor) sp_count_r
-  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
-  close(unitoutfactor)
 
 
Index: src/concoutput_surf.f90
===================================================================
--- src/concoutput_surf.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/concoutput_surf.f90	(revision 16b61a5e2d3b1020be1cff080e479ba26f3f7063)
@@ -72,7 +72,4 @@
   real :: sp_fact
   real :: outnum,densityoutrecept(maxreceptor),xl,yl
-! RLT
-  real :: densitydryrecept(maxreceptor)
-  real :: factor_dryrecept(maxreceptor)
 
 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
@@ -104,5 +101,4 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
-  logical :: lexist
 
 
@@ -184,7 +180,4 @@
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
              rho(iix,jjy,kzz-1,2)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
       end do
     end do
@@ -197,12 +190,6 @@
     jjy=max(min(nint(yl),nymin1),0)
     densityoutrecept(i)=rho(iix,jjy,1,2)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
   end do
 
-! RLT
-! conversion factor for output relative to dry air
-  factor_drygrid=densityoutgrid/densitydrygrid
-  factor_dryrecept=densityoutrecept/densitydryrecept
 
 ! Output is different for forward and backward simulations
@@ -609,47 +596,4 @@
   end do
 
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
-  if (lexist) then
-    ! open and append
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
-            status='old',action='write',access='append')
-  else
-    ! create new
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
-            status='new',action='write')
-  endif
-  sp_count_i=0
-  sp_count_r=0
-  sp_fact=-1.
-  sp_zer=.true.
-  do kz=1,1
-    do jy=0,numygrid-1
-      do ix=0,numxgrid-1
-        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
-          if (sp_zer.eqv..true.) then ! first value not equal to one
-            sp_count_i=sp_count_i+1
-            sparse_dump_i(sp_count_i)= &
-                  ix+jy*numxgrid+kz*numxgrid*numygrid
-            sp_zer=.false.
-            sp_fact=sp_fact*(-1.)
-          endif
-          sp_count_r=sp_count_r+1
-          sparse_dump_r(sp_count_r)= &
-               sp_fact*factor_drygrid(ix,jy,kz)
-        else ! factor is one
-          sp_zer=.true.
-        endif
-      end do
-    end do
-  end do
-  write(unitoutfactor) sp_count_i
-  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
-  write(unitoutfactor) sp_count_r
-  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
-  close(unitoutfactor)
-
-
   if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
   if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
@@ -678,21 +622,5 @@
   endif
 
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  if (numreceptor.gt.0) then
-    inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
-     if (lexist) then
-     ! open and append
-      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
-              status='old',action='write',access='append')
-    else
-      ! create new
-      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
-              status='new',action='write')
-    endif
-    write(unitoutfactor) itime
-    write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
-    close(unitoutfactor)
-  endif
+
 
 ! Reinitialization of grid
Index: src/concoutput_surf_nest.f90
===================================================================
--- src/concoutput_surf_nest.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/concoutput_surf_nest.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
@@ -70,7 +70,4 @@
   real :: sp_fact
   real :: outnum,densityoutrecept(maxreceptor),xl,yl
-! RLT
-  real :: densitydryrecept(maxreceptor)
-  real :: factor_dryrecept(maxreceptor)
 
   !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
@@ -99,5 +96,5 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
-  logical :: lexist
+
 
   ! Determine current calendar date, needed for the file name
@@ -162,7 +159,4 @@
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
              rho(iix,jjy,kzz-1,2)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
       end do
     end do
@@ -175,12 +169,6 @@
       jjy=max(min(nint(yl),nymin1),0)
       densityoutrecept(i)=rho(iix,jjy,1,2)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
     end do
 
-! RLT
-! conversion factor for output relative to dry air
-  factor_drygrid=densityoutgrid/densitydrygrid
-  factor_dryrecept=densityoutrecept/densitydryrecept
 
   ! Output is different for forward and backward simulations
@@ -329,6 +317,6 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
 
   ! Dry deposition
@@ -366,6 +354,6 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
 
 
@@ -411,6 +399,6 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
          else
 
@@ -452,6 +440,6 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
          endif ! surf_only
 
@@ -499,6 +487,6 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
 
 
@@ -538,6 +526,6 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
 
 
@@ -582,6 +570,6 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
          else
 
@@ -623,6 +611,6 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
          endif ! surf_only
 
@@ -636,46 +624,4 @@
 
   end do
-
-! RLT Aug 2017
-! Write out conversion factor for dry air
-  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
-  if (lexist) then
-    ! open and append
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
-            status='old',action='write',access='append')
-  else
-    ! create new
-    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
-            status='new',action='write')
-  endif
-  sp_count_i=0
-  sp_count_r=0
-  sp_fact=-1.
-  sp_zer=.true.
-  do kz=1,1
-    do jy=0,numygridn-1
-      do ix=0,numxgridn-1
-        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
-          if (sp_zer.eqv..true.) then ! first value not equal to one
-            sp_count_i=sp_count_i+1
-            sparse_dump_i(sp_count_i)= &
-                  ix+jy*numxgridn+kz*numxgridn*numygridn
-            sp_zer=.false.
-            sp_fact=sp_fact*(-1.)
-          endif
-          sp_count_r=sp_count_r+1
-          sparse_dump_r(sp_count_r)= &
-               sp_fact*factor_drygrid(ix,jy,kz)
-        else ! factor is one
-          sp_zer=.true.
-        endif
-      end do
-    end do
-  end do
-  write(unitoutfactor) sp_count_i
-  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
-  write(unitoutfactor) sp_count_r
-  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
-  close(unitoutfactor)
 
 
Index: src/getfields.f90
===================================================================
--- src/getfields.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/getfields.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -87,8 +87,4 @@
   real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
   real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
-! RLT added partial pressure water vapor 
-  real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
-  integer :: kz, ix
-  character(len=100) :: rowfmt
 
   integer :: indmin = 1
@@ -157,7 +153,7 @@
 40  indmin=indj
 
-    if (WETBKDEP) then
-      call writeprecip(itime,memind(1))
-    endif
+   if (WETBKDEP) then
+        call writeprecip(itime,memind(1))
+   endif
 
   else
@@ -208,32 +204,9 @@
 60  indmin=indj
 
-    if (WETBKDEP) then
-      call writeprecip(itime,memind(1))
-    endif
-
-  end if
-
-  ! RLT calculate dry air density
-  pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv)
-  rho_dry=(prs-pwater)/(r_air*tt)
-
-  ! test density
-!  write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))'
-!  if(itime.eq.0) then
-!    open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write')
-!    do kz=1,nzmax
-!      do ix=1,nxmax
-!        write(500,fmt=rowfmt) rho_dry(ix,:,kz,1)
-!      end do
-!    end do
-!    close(500)
-!    open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write')
-!    do kz=1,nzmax
-!      do ix=1,nxmax
-!        write(500,fmt=rowfmt) rho(ix,:,kz,1)
-!      end do
-!    end do
-!    close(500)
-!  endif
+   if (WETBKDEP) then
+        call writeprecip(itime,memind(1))
+   endif
+
+  endif
 
   lwindinterv=abs(memtime(2)-memtime(1))
Index: src/init_domainfill.f90
===================================================================
--- src/init_domainfill.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/init_domainfill.f90	(revision b5127f9034f2a15f9a3ce226fd26edfe845311de)
@@ -86,8 +86,4 @@
     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
@@ -418,5 +414,5 @@
 !***************************************************************************
 
-  if ((ipin.eq.1).and.(.not.gdomainfill)) then
+  if (ipin.eq.1) 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 328fdf90784e3387eb9745d1009d7523f4ed1976)
+++ src/init_domainfill_mpi.f90	(revision b5127f9034f2a15f9a3ce226fd26edfe845311de)
@@ -110,8 +110,4 @@
   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
@@ -217,4 +213,5 @@
         colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)
         colmasstotal=colmasstotal+colmass(ix,jy)
+
       end do
     end do
@@ -469,5 +466,5 @@
 
 ! eso TODO: only needed for root process
-    if ((ipin.eq.1).and.(.not.gdomainfill)) then
+    if (ipin.eq.1) then
       open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
            form='unformatted')
@@ -477,31 +474,25 @@
     endif
 
-    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
+    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))
     
-  end if ! end if(lroot)
-
+  end if ! end if(lroot)  
 
 
 ! Distribute particles to other processes (numpart is 'per-process', not total)
-! 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)
+  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)
 
 ! 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: c/initial_cond_output_inversion.f90
===================================================================
--- src/initial_cond_output_inversion.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ 	(revision )
@@ -1,192 +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 initial_cond_output_inversion(itime)
-  !                                 i
-  !*****************************************************************************
-  !                                                                            *
-  !     Output of the initial condition sensitivity field.                     *
-  !                                                                            *
-  !     Author: A. Stohl                                                       *
-  !                                                                            *
-  !     24 May 1995                                                            *
-  !                                                                            *
-  !     13 April 1999, Major update: if output size is smaller, dump output    *
-  !                    in sparse matrix format; additional output of           *
-  !                    uncertainty                                             *
-  !                                                                            *
-  !     05 April 2000, Major update: output of age classes; output for backward*
-  !                    runs is time spent in grid cell times total mass of     *
-  !                    species.                                                *
-  !                                                                            *
-  !     17 February 2002, Appropriate dimensions for backward and forward runs *
-  !                       are now specified in file par_mod                    *
-  !                                                                            *
-  !     June 2006, write grid in sparse matrix with a single write command     *
-  !                in order to save disk space                                 *
-  !                                                                            *
-  !     2008 new sparse matrix format                                          *
-  !                                                                            *
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  ! ncells          number of cells with non-zero concentrations               *
-  ! sparse          .true. if in sparse matrix format, else .false.            *
-  !                                                                            *
-  !*****************************************************************************
-
-  use unc_mod
-  use point_mod
-  use outg_mod
-  use par_mod
-  use com_mod
-
-  implicit none
-
-  integer :: itime,i,ix,jy,kz,ks,kp,sp_count_i,sp_count_r
-  integer :: jjjjmmdd, ihmmss
-  real(kind=dp) :: jul
-  real :: sp_fact,fact_recept
-  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
-  logical :: sp_zer,lexist
-  logical,save :: listart=.true.
-  logical,save,allocatable,dimension(:) :: listartrel
-  character :: adate*8,atime*6
-  character :: areldate*8,areltime*6
-  character(len=3) :: anspec
-
-  if(listart) then
-    allocate(listartrel(maxpointspec_act))
-    listartrel(:)=.true.
-  endif
-  print*, 'listartrel = ',listartrel
-
-  !*********************************************************************
-  ! Determine the standard deviation of the mean concentration or mixing
-  ! ratio (uncertainty of the output) and the dry and wet deposition
-  !*********************************************************************
-
-  do ks=1,nspec
-
-    write(anspec,'(i3.3)') ks
-
-    do kp=1,maxpointspec_act
-
-      ! calculate date of release
-      jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp    ! this is the current day
-      call caldate(jul,jjjjmmdd,ihmmss)
-      write(areldate,'(i8.8)') jjjjmmdd
-      write(areltime,'(i6.6)') ihmmss
-      print*, areldate//areltime
-
-      ! calculate date of field
-      jul=bdate+real(itime,kind=dp)/86400._dp
-      call caldate(jul,jjjjmmdd,ihmmss)
-      write(adate,'(i8.8)') jjjjmmdd
-      write(atime,'(i6.6)') ihmmss
-      print*, adate//atime
-
-      inquire(file=path(2)(1:length(2))//'grid_initial_'//areldate// &
-         areltime//'_'//anspec,exist=lexist)
-      if(lexist.and..not.listartrel(kp)) then
-        ! open and append to existing file
-        open(97,file=path(2)(1:length(2))//'grid_initial_'//areldate// &
-             areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append') 
-      else
-        ! open new file
-        open(97,file=path(2)(1:length(2))//'grid_initial_'//areldate// &
-             areltime//'_'//anspec,form='unformatted',status='replace',action='write')
-      endif
-      write(97) jjjjmmdd
-      write(97) ihmmss
-
-      listartrel(kp)=.false.
-
-      if (ind_rel.eq.1) then
-        fact_recept=rho_rel(kp)
-      else
-        fact_recept=1.
-      endif
-
-  !*******************************************************************
-  ! Generate output: may be in concentration (ng/m3) or in mixing
-  ! ratio (ppt) or both
-  ! Output the position and the values alternated multiplied by
-  ! 1 or -1, first line is number of values, number of positions
-  ! For backward simulations, the unit is seconds, stored in grid_time
-  !*******************************************************************
-
-  ! Write out dummy "wet and dry deposition" fields, to keep same format
-  ! as for concentration output
-!      sp_count_i=0
-!      sp_count_r=0
-!      write(97) sp_count_i
-!      write(97) (sparse_dump_i(i),i=1,sp_count_i)
-!      write(97) sp_count_r
-!      write(97) (sparse_dump_r(i),i=1,sp_count_r)
-!      write(97) sp_count_i
-!      write(97) (sparse_dump_i(i),i=1,sp_count_i)
-!      write(97) sp_count_r
-!      write(97) (sparse_dump_r(i),i=1,sp_count_r)
-
-
-  ! Write out sensitivity to initial conditions
-      sp_count_i=0
-      sp_count_r=0
-      sp_fact=-1.
-      sp_zer=.true.
-      do kz=1,numzgrid
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-            if (init_cond(ix,jy,kz,ks,kp).gt.smallnum) then
-              if (sp_zer.eqv..true.) then ! first non zero value
-                sp_count_i=sp_count_i+1
-                sparse_dump_i(sp_count_i)= &
-                     ix+jy*numxgrid+kz*numxgrid*numygrid
-                sp_zer=.false.
-                sp_fact=sp_fact*(-1.)
-              endif
-              sp_count_r=sp_count_r+1
-              sparse_dump_r(sp_count_r)=sp_fact* &
-                   init_cond(ix,jy,kz,ks,kp)/xmass(kp,ks)*fact_recept
-            else ! concentration is zero
-              sp_zer=.true.
-            endif
-          end do
-        end do
-      end do
-      write(97) sp_count_i
-      write(97) (sparse_dump_i(i),i=1,sp_count_i)
-      write(97) sp_count_r
-      write(97) (sparse_dump_r(i),i=1,sp_count_r)
-
-      close(97)
-
-    end do
-
-  end do
-
-  ! reset listart
-  if (listart) then
-    listart=.false.
-  endif
-
-end subroutine initial_cond_output_inversion
Index: src/makefile
===================================================================
--- src/makefile	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/makefile	(revision 7123c700d7191876c1d22402ad01729c688c6e56)
@@ -118,5 +118,4 @@
 OBJECTS_SERIAL = \
 	releaseparticles.o 	partoutput.o \
-	partoutput_average.o \
 	conccalc.o \
 	init_domainfill.o 	concoutput.o  \
@@ -128,6 +127,4 @@
 	redist.o		\
 	concoutput_surf.o	concoutput_surf_nest.o	\
-	concoutput_inversion_nest.o 	\
-	concoutput_inversion.o \
 	getfields.o \
         readwind_ecmwf.o
@@ -135,5 +132,5 @@
 ## For MPI version 
 OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
-	partoutput_average_mpi.o conccalc_mpi.o \
+	conccalc_mpi.o \
 	init_domainfill_mpi.o concoutput_mpi.o 	\
 	timemanager_mpi.o FLEXPART_MPI.o	\
@@ -152,5 +149,5 @@
 advance.o		initialize.o		\
 writeheader.o		writeheader_txt.o  	\
-partpos_average.o	writeprecip.o \
+writeprecip.o \
 writeheader_surf.o     	assignland.o\
 part0.o 		gethourlyOH.o\
@@ -202,5 +199,5 @@
 drydepokernel_nest.o    zenithangle.o \
 ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o initial_cond_output_inversion.o \
+initial_cond_calc.o     initial_cond_output.o \
 dynamic_viscosity.o     get_settling.o	\
 initialize_cbl_vel.o	re_initialize_particle.o \
@@ -274,9 +271,7 @@
 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 mean_mod.o
-concoutput_inversion.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 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_inversion_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 mean_mod.o
@@ -325,5 +320,4 @@
 initial_cond_calc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
 initial_cond_output.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
-initial_cond_output_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
 initialize.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o random_mod.o
 initialize_cbl_vel.o: com_mod.o par_mod.o random_mod.o
@@ -354,8 +348,5 @@
 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 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ src/mpi_mod.f90	(revision 0ecc1fed9f7c7c07cef687094373915d51c69489)
@@ -88,5 +88,4 @@
 ! 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
 
@@ -126,5 +125,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.
@@ -191,6 +190,6 @@
 !   mp_np       number of running processes, decided at run-time
 !***********************************************************************
-    use par_mod, only: maxpart, numwfmem, dep_prec, maxreceptor, maxspec
-    use com_mod, only: mpi_mode, verbosity, creceptor0
+    use par_mod, only: maxpart, numwfmem, dep_prec
+    use com_mod, only: mpi_mode, verbosity
 
     implicit none
@@ -338,5 +337,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
@@ -366,14 +365,4 @@
     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
 
@@ -570,5 +559,5 @@
 ! invalid particles at the end of the arrays
 
-601 do i=numpart, 1, -1
+601 do i=num_part, 1, -1
       if (itra1(i).eq.-999999999) then
         numpart=numpart-1
@@ -609,5 +598,5 @@
     integer :: i,jj,nn, num_part=1,m,imin, num_trans
     logical :: first_iter
-    integer,allocatable,dimension(:) :: idx_arr
+    integer,allocatable,dimension(:) :: numparticles_mpi, idx_arr
     real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this
 
@@ -618,7 +607,8 @@
 ! All processes exchange information on number of particles
 !****************************************************************************
-    allocate( idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
-
-    call MPI_Allgather(numpart, 1, MPI_INTEGER, npart_per_process, &
+    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, &
          & 1, MPI_INTEGER, mp_comm_used, mp_ierr)
 
@@ -626,11 +616,8 @@
 ! Sort from lowest to highest
 ! Initial guess: correct order
-    sorted(:) = npart_per_process(:)
+    sorted(:) = numparticles_mpi(:)
     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
@@ -658,11 +645,11 @@
     m=mp_partgroup_np-1 ! index for last sorted process (most particles)
     do i=0,mp_partgroup_np/2-1
-      num_trans = npart_per_process(idx_arr(m)) - npart_per_process(idx_arr(i))
+      num_trans = numparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))
       if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) 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
+        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
 ! DBG
-          ! 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
+          ! 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
 ! DBG
           call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)
@@ -672,5 +659,5 @@
     end do
 
-    deallocate(idx_arr, sorted)
+    deallocate(numparticles_mpi, idx_arr, sorted)
 
   end subroutine mpif_calculate_part_redist
@@ -1974,5 +1961,5 @@
     if (readclouds) then
       j=j+1
-      call MPI_Irecv(ctwc(:,:,mind),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
+      call MPI_Irecv(ctwc(:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
            &MPI_COMM_WORLD,reqs(j),mp_ierr)
       if (mp_ierr /= 0) goto 600
@@ -2339,5 +2326,5 @@
       if (readclouds) then
         j=j+1
-        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
+        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
              &MPI_COMM_WORLD,reqs(j),mp_ierr)
         if (mp_ierr /= 0) goto 600
@@ -2475,5 +2462,25 @@
     end if
 
-  ! Receptor concentrations    
+#else
+
+      call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (lroot) gridunc = gridunc0
+
+#endif
+
+    if ((WETDEP).and.(ldirect.gt.0)) then
+      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    end if
+
+    if ((DRYDEP).and.(ldirect.gt.0)) then
+      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    end if
+
+! Receptor concentrations    
     if (lroot) then
       call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
@@ -2484,31 +2491,4 @@
            & 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
-
-    if ((WETDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
-           & mp_comm_used, mp_ierr)
-      if (mp_ierr /= 0) goto 600
-    end if
-
-    if ((DRYDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
-           & mp_comm_used, mp_ierr)
-      if (mp_ierr /= 0) goto 600
-    end if
-
 
     if (mp_measure_time) call mpif_mtime('commtime',1)
@@ -2720,17 +2700,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')
@@ -2808,8 +2788,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 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/netcdf_output_mod.f90	(revision 4ad96c54b793af523a98633b88f07aa60c69ea47)
@@ -93,5 +93,5 @@
   character(len=255), parameter :: institution = 'NILU'
 
-  integer            :: tpointer=0
+  integer            :: tpointer
   character(len=255) :: ncfname, ncfnamen
 
Index: src/outg_mod.f90
===================================================================
--- src/outg_mod.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/outg_mod.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
@@ -37,6 +37,4 @@
   real,allocatable, dimension (:,:,:) :: areanorth
   real,allocatable, dimension (:,:,:) :: densityoutgrid
-  real,allocatable, dimension (:,:,:) :: densitydrygrid ! added RLT 
-  real,allocatable, dimension (:,:,:) :: factor_drygrid ! added RLT 
   real,allocatable, dimension (:,:,:) :: factor3d
   real,allocatable, dimension (:,:,:) :: grid
Index: src/outgrid_init.f90
===================================================================
--- src/outgrid_init.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/outgrid_init.f90	(revision d2a5a833066965c4c4b33df109fd4e461f7275bb)
@@ -268,11 +268,4 @@
     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
   allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, &
-       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
-! RLT
-  allocate(densitydrygrid(0:max(numxgrid,numxgridn)-1, &
-       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
-  allocate(factor_drygrid(0:max(numxgrid,numxgridn)-1, &
        0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/par_mod.f90	(revision 79e0349a50fd1ee6dfdf4b044ade2b6f1e702855)
@@ -30,4 +30,7 @@
 !        Update 15 August 2013 IP                                              *
 !                                                                              *
+!        ESO 2016:                                                             *
+!          GFS specific parameters moved to gfs_mod.f90                        *
+!          ECMWF specific parameters moved to ecmwf_mod.f90                    *
 !                                                                              *
 !*******************************************************************************
@@ -77,7 +80,4 @@
   real,parameter :: pi=3.14159265, r_earth=6.371e6, r_air=287.05, ga=9.81
   real,parameter :: cpa=1004.6, kappa=0.286, pi180=pi/180., vonkarman=0.4
-  ! additional constants RLT Aug-2017
-  real,parameter :: rgas=8.31447 
-  real,parameter :: r_water=461.495
 
   ! pi                      number "pi"
@@ -89,6 +89,4 @@
   ! kappa                   exponent of formula for potential temperature
   ! vonkarman               von Karman constant
-  ! rgas                    universal gas constant [J/mol/K]
-  ! r_water                 specific gas constant for water vapor [J/kg/K]
 
   real,parameter :: karman=0.40, href=15., convke=2.0
@@ -282,5 +280,5 @@
 
   integer,parameter :: unitpath=1, unitcommand=1, unitageclasses=1, unitgrid=1
-  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93, unitpartout_average=105
+  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93
   integer,parameter :: unitpartin=93, unitflux=98, unitouttraj=96
   integer,parameter :: unitvert=1, unitoro=1, unitpoin=1, unitreceptor=1
@@ -292,6 +290,4 @@
   integer,parameter :: unitboundcond=89
   integer,parameter :: unittmp=101
-! RLT
-  integer,parameter :: unitoutfactor=102
 
 !******************************************************
Index: src/partoutput.f90
===================================================================
--- src/partoutput.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/partoutput.f90	(revision d2a5a833066965c4c4b33df109fd4e461f7275bb)
@@ -71,5 +71,5 @@
   !**************************************
 
-  if (ipout.eq.1.or.ipout.eq.3) then
+  if (ipout.eq.1) then
     open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
          atime,form='unformatted')
Index: c/partoutput_average.f90
===================================================================
--- src/partoutput_average.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ 	(revision )
@@ -1,221 +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 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: c/partoutput_average_mpi.f90
===================================================================
--- src/partoutput_average_mpi.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ 	(revision )
@@ -1,230 +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 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 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/partoutput_mpi.f90	(revision d2a5a833066965c4c4b33df109fd4e461f7275bb)
@@ -78,5 +78,5 @@
   !**************************************
 
-  if (ipout.eq.1.or.ipout.eq.3) then
+  if (ipout.eq.1) then
     open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
          atime,form='unformatted',status=file_stat,position='append')
Index: c/partpos_average.f90
===================================================================
--- src/partpos_average.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ 	(revision )
@@ -1,206 +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 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 2eefa58ac797f119acc28fade41beff352006040)
+++ src/readcommand.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -50,5 +50,4 @@
   ! 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      *
@@ -98,5 +97,4 @@
   iout, &
   ipout, &
-  ipoutfac, &
   lsubgrid, &
   lconvection, &
@@ -114,5 +112,4 @@
   surf_only, &
   cblflag, &
-  linversionout, &
   ohfields_path
 
@@ -132,5 +129,4 @@
   iout=3
   ipout=0
-  ipoutfac=1
   lsubgrid=1
   lconvection=1
@@ -148,5 +144,4 @@
   surf_only=0 
   cblflag=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine
-  linversionout=0
   ohfields_path="../../flexin/"
 
@@ -416,5 +411,5 @@
   !**********************************************************************
 
-  if ((iout.lt.1).or.(iout.gt.5)) then
+  if ((iout.lt.1).or.(iout.gt.6)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
     write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR        #### '
@@ -476,14 +471,4 @@
       write(*,*) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR      ####'
       write(*,*) '#### EACH RELEASE IS FORBIDDEN !             ####'
-      stop
-  endif
-
-  ! Inversion output format only for backward runs
-  !*****************************************************************************
-  
-  if ((linversionout.eq.1).and.(ldirect.eq.1)) then
-      write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND:     ####'
-      write(*,*) '#### INVERSION OUTPUT FORMAT ONLY FOR        ####'
-      write(*,*) '#### BACKWARD RUNS                           ####'
       stop
   endif
@@ -522,7 +507,7 @@
   !****************************************************************
 
-  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2).and.(ipout.ne.3)) then
+  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
-    write(*,*) ' #### IPOUT MUST BE 0, 1, 2 OR 3!             #### '
+    write(*,*) ' #### IPOUT MUST BE 1, 2 OR 3!                #### '
     stop
   endif
Index: src/readpaths.f90
===================================================================
--- src/readpaths.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/readpaths.f90	(revision 62e65c73ae7e3d8c4b8b884065ede76e07211ec8)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine readpaths
+subroutine readpaths !(pathfile)
 
   !*****************************************************************************
Index: src/readwind_gfs.f90
===================================================================
--- src/readwind_gfs.f90	(revision 5d74ed914bb3649af12dc589a95068bd49aa0fdb)
+++ src/readwind_gfs.f90	(revision db91eb78afeb1bb8bafe66068877db8841859062)
@@ -83,5 +83,5 @@
 
   ! NCEP
-  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
+  integer :: numpt,numpu,numpv,numpw,numprh
   real :: help, temp, ew
   real :: elev
@@ -134,5 +134,4 @@
   numpw=0
   numprh=0
-  numpclwch=0
   ifield=0
 10   ifield=ifield+1
@@ -558,20 +557,8 @@
       endif
 ! SEC & IP 12/2018 read GFS clouds
-      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
+      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)
         readclouds=.true.
         sumclouds=.true.
-!        readclouds=.false.
-!       sumclouds=.false.
       endif
 
Index: src/releaseparticles.f90
===================================================================
--- src/releaseparticles.f90	(revision 7873bf796f5ee67aab272b1243107198ceb275b1)
+++ src/releaseparticles.f90	(revision 75a4dedbe36ebfba4e9c67fce8af8641dafcd8eb)
@@ -114,6 +114,5 @@
       average_timecorrect=0.
       do k=1,nspec
-        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
+        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 2eefa58ac797f119acc28fade41beff352006040)
+++ src/timemanager.f90	(revision c7d1052e9b4a8319033f957f6b77184654bca402)
@@ -408,14 +408,5 @@
 #endif
             else
-              if (linversionout.eq.1) then
-                call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-                if (verbosity.eq.1) then
-                  print*,'called concoutput_inversion'
-                  call system_clock(count_clock)
-                  write(*,*) 'system clock',count_clock - count_clock0 
-                endif
-              else
               call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-              endif
               if (verbosity.eq.1) then
                 print*,'called concoutput_surf '
@@ -431,9 +422,5 @@
                 call concoutput_nest(itime,outnum)
               else 
-                if(linversionout.eq.1) then
-                  call concoutput_inversion_nest(itime,outnum)
-                else 
                 call concoutput_surf_nest(itime,outnum)
-              endif
               endif
             else
@@ -464,8 +451,5 @@
 45      format(i13,' Seconds simulated: ',i13, ' Particles:    Uncertainty: ',3f7.3)
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
-        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
+        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
         loutnext=loutnext+loutstep
         loutstart=loutnext-loutaver/2
@@ -625,10 +609,4 @@
 !        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
   !***************************************************
@@ -752,11 +730,5 @@
   if (ipout.eq.2) call partoutput(itime)     ! dump particle positions
 
-  if (linit_cond.ge.1) then
-    if(linversionout.eq.1) then
-      call initial_cond_output_inversion(itime)   ! dump initial cond. field
-    else
-      call initial_cond_output(itime)   ! dump initial cond. fielf
-    endif
-  endif
+  if (linit_cond.ge.1) call initial_cond_output(itime)   ! dump initial cond. field
 
   !close(104)
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ src/timemanager_mpi.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
@@ -113,5 +113,5 @@
   integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
 ! integer :: ksp
-  integer :: ip,irec
+  integer :: ip
   integer :: loutnext,loutstart,loutend
   integer :: ix,jy,ldeltat,itage,nage,idummy
@@ -129,5 +129,4 @@
 ! Measure time spent in timemanager
   if (mp_measure_time) call mpif_mtime('timemanager',0)
-
 
 ! First output for time 0
@@ -533,5 +532,5 @@
                 end if
 
-              else
+              else  ! :TODO: check for zeroing in the netcdf module
                 call concoutput_surf_nest(itime,outnum)
               end if
@@ -594,15 +593,8 @@
 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) 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)
+            if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions
             call mpif_mpi_barrier
           end do
-          if (mp_measure_time) call mpif_mtime('iotime',1)
         end if
 
@@ -765,9 +757,4 @@
         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
@@ -908,5 +895,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 2eefa58ac797f119acc28fade41beff352006040)
+++ src/verttransform_ecmwf.f90	(revision 79e0349a50fd1ee6dfdf4b044ade2b6f1e702855)
@@ -73,4 +73,5 @@
   use com_mod
   use cmapf_mod, only: cc2gll
+!  use mpi_mod
 
   implicit none
@@ -81,6 +82,4 @@
   real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
   real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
-  ! RLT added pressure
-  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: prsh
   real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
   real,dimension(0:nymax-1) :: cosf
@@ -220,4 +219,5 @@
 !*************************
 
+
   do jy=0,nymin1
     do ix=0,nxmin1
@@ -230,6 +230,4 @@
   wzlev(:,:,1)=0.
   rhoh(:,:,1)=pold/(r_air*tvold)
-  ! RLT add pressure
-  prsh(:,:,1)=ps(:,:,1,n)
 
 
@@ -239,6 +237,4 @@
   do kz=2,nuvz
     pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
-    ! RLT add pressure
-    prsh(:,:,kz)=pint
     tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
     rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
@@ -292,6 +288,4 @@
   pv(:,:,1,n)=pvh(:,:,1)
   rho(:,:,1,n)=rhoh(:,:,1)
-! RLT add pressure
-  prs(:,:,1,n)=prsh(:,:,1)
 
   uu(:,:,nz,n)=uuh(:,:,nuvz)
@@ -307,6 +301,5 @@
   pv(:,:,nz,n)=pvh(:,:,nuvz)
   rho(:,:,nz,n)=rhoh(:,:,nuvz)
-! RLT
-  prs(:,:,nz,n)=prsh(:,:,nuvz)
+
 
   kmin=2
@@ -328,6 +321,4 @@
           pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
           rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
-! RLT
-          prs(ix,jy,iz,n)=prs(ix,jy,nz,n)
         else
           innuvz: do kz=idx(ix,jy),nuvz
@@ -363,6 +354,4 @@
           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
           rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
-! RLT add pressure
-          prs(ix,jy,iz,n)=(prsh(ix,jy,kz-1)*dz2+prsh(ix,jy,kz)*dz1)/dz
         endif
       enddo
@@ -665,5 +654,5 @@
         if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
 
-          do kz=nz,2,-1 !go Bottom up!
+          do kz=nz,1,-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 437c545fb9b9ece116fae2e06166d88b979f91a7)
+++ src/verttransform_gfs.f90	(revision db91eb78afeb1bb8bafe66068877db8841859062)
@@ -548,10 +548,10 @@
         if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
 
-          do kz=nz,2,-1 !go Bottom up!
+          do kz=nz,1,-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
