Index: create_tarball.sh
===================================================================
--- create_tarball.sh	(revision a816416efbe37c2603b5f0b62c029f89c8af8ac0)
+++ create_tarball.sh	(revision a2e9de41663ac7ef6e5a44c759e9775b1128a61a)
@@ -2,14 +2,22 @@
 
 #define version number
-version=10.3beta5
+
+githash=$(git rev-parse --short --verify HEAD)
+
+
+version=10.3beta5_$githash
 
 # define tarball name
 targetdir=../flexpart_distribution/
 tarball_tmp=${targetdir}flexpart_v$version 
+#tarball=${targetdir}flexpart_v$version.tar 
+tarball=${tarball_tmp}.tar 
 
 # clean old package
 if [ -d $tarball_tmp ]; then
-  echo $tarball_tmp exists: move to $tarball_tmp.bk   
-  mv $tarball_tmp ${tarball_tmp}.bk 
+  echo $tarball_tmp exists: move to $tarball_tmp.bk and exit  
+  mkdir $tarball_tmp.bk 
+  mv $tarball_tmp ${tarball_tmp}.bk/ 
+  mv $tarball ${tarball_tmp}.bk/ 
   exit 
 fi
@@ -79,4 +87,13 @@
 #cp -r flex_ecmwf_src/* $tarball_tmp/preprocess/flex_extract/
 ## cp -r flex_extract/work/EA* $tarball_tmp/preprocess/flex_extract/work   
+
+echo include flex_extract v7.0.4 b7c1c04a204c91e53759ef590504bf52dfaece64
+flex_extract=../flex_extract_v7.0.4/
+cp $flex_extract/README.md $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/docs $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/grib_templates $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/python $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/src $tarball_tmp/preprocess/flex_extract
+
 
 
@@ -175,5 +192,8 @@
 # cp -r tests/default_cases $tarball_tmp/tests/
 
-
+tar cvf  $tarball  $tarball_tmp  
+
+echo  $tarball complete
+echo exported untarred files in $tarball_tmp 
 exit
 #return
Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/FLEXPART.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -137,5 +137,5 @@
     write(*,*) 'call readpaths'
   endif 
-  call readpaths(pathfile)
+  call readpaths
  
   if (verbosity.gt.1) then !show clock info 
@@ -452,5 +452,7 @@
   endif
 
+  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
   call timemanager(metdata_format)
+ 
 
   if (verbosity.gt.0) then
@@ -467,5 +469,4 @@
       endif
     end do
-    write (*,*) 'timemanager> call wetdepo'
   endif
   
Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision 0c8c7f2f1bd07179b6e1be5a302d9e2d52633f9e)
+++ src/FLEXPART_MPI.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -146,5 +146,5 @@
     write(*,*) 'call readpaths'
   endif 
-  call readpaths(pathfile)
+  call readpaths
  
   if (verbosity.gt.1) then !show clock info 
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/com_mod.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -124,4 +124,7 @@
   ! 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)
 
@@ -176,5 +179,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)
+  real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
   real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
 
@@ -361,5 +364,7 @@
   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)
@@ -383,4 +388,5 @@
   ! 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
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision 20963b15051b22eeb381dcf4e824a9104d92717b)
+++ src/concoutput.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -72,4 +72,7 @@
   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),
@@ -106,4 +109,5 @@
   character(LEN=8),save :: file_stat='REPLACE'
   logical :: ldates_file
+  logical :: lexist
   integer :: ierr
   character(LEN=100) :: dates_char
@@ -203,4 +207,7 @@
         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
@@ -214,6 +221,12 @@
 !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
@@ -353,5 +366,5 @@
 ! Concentration output
 !*********************
-        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
 
 ! Wet deposition
@@ -614,4 +627,47 @@
   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/ &
@@ -640,5 +696,21 @@
   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_inversion.f90
===================================================================
--- src/concoutput_inversion.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/concoutput_inversion.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -0,0 +1,855 @@
+!**********************************************************************
+! 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: src/concoutput_inversion_nest.f90
===================================================================
--- src/concoutput_inversion_nest.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/concoutput_inversion_nest.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -0,0 +1,699 @@
+!**********************************************************************
+! 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 20963b15051b22eeb381dcf4e824a9104d92717b)
+++ src/concoutput_mpi.f90	(revision 6741557b2e77046144879c12ca45ea39587f7099)
@@ -364,5 +364,5 @@
 ! Concentration output
 !*********************
-        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
 
 ! Wet deposition
Index: src/concoutput_nest.f90
===================================================================
--- src/concoutput_nest.f90	(revision 7d02c2f5693a56e1b4e9b3923f9343c0a5e6b1fa)
+++ src/concoutput_nest.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -70,4 +70,7 @@
   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),
@@ -96,4 +99,5 @@
   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
@@ -164,4 +168,7 @@
         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
@@ -175,6 +182,12 @@
     !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
@@ -551,4 +564,45 @@
   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 16b61a5e2d3b1020be1cff080e479ba26f3f7063)
+++ src/concoutput_surf.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -72,4 +72,7 @@
   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),
@@ -101,4 +104,5 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  logical :: lexist
 
 
@@ -180,4 +184,7 @@
         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
@@ -190,6 +197,12 @@
     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
@@ -596,4 +609,47 @@
   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/ &
@@ -622,5 +678,21 @@
   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 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/concoutput_surf_nest.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -70,4 +70,7 @@
   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),
@@ -96,5 +99,5 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
-
+  logical :: lexist
 
   ! Determine current calendar date, needed for the file name
@@ -159,4 +162,7 @@
         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
@@ -169,6 +175,12 @@
       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
@@ -317,6 +329,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
@@ -354,6 +366,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)
 
 
@@ -399,6 +411,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
 
@@ -440,6 +452,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
 
@@ -487,6 +499,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)
 
 
@@ -526,6 +538,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)
 
 
@@ -570,6 +582,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
 
@@ -611,6 +623,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
 
@@ -624,4 +636,46 @@
 
   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 6ecb30a41fb8a07e153abadc712699a09c15f90c)
+++ src/getfields.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -87,4 +87,8 @@
   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
@@ -153,7 +157,7 @@
 40  indmin=indj
 
-   if (WETBKDEP) then
-        call writeprecip(itime,memind(1))
-   endif
+    if (WETBKDEP) then
+      call writeprecip(itime,memind(1))
+    endif
 
   else
@@ -204,9 +208,32 @@
 60  indmin=indj
 
-   if (WETBKDEP) then
-        call writeprecip(itime,memind(1))
-   endif
-
-  endif
+    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
 
   lwindinterv=abs(memtime(2)-memtime(1))
Index: src/initial_cond_output_inversion.f90
===================================================================
--- src/initial_cond_output_inversion.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
+++ src/initial_cond_output_inversion.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -0,0 +1,192 @@
+!**********************************************************************
+! 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 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/makefile	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -128,4 +128,6 @@
 	redist.o		\
 	concoutput_surf.o	concoutput_surf_nest.o	\
+	concoutput_inversion_nest.o 	\
+	concoutput_inversion.o \
 	getfields.o \
         readwind_ecmwf.o
@@ -200,5 +202,5 @@
 drydepokernel_nest.o    zenithangle.o \
 ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
+initial_cond_calc.o     initial_cond_output.o initial_cond_output_inversion.o \
 dynamic_viscosity.o     get_settling.o	\
 initialize_cbl_vel.o	re_initialize_particle.o \
@@ -272,7 +274,9 @@
 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
@@ -321,4 +325,5 @@
 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
Index: src/outg_mod.f90
===================================================================
--- src/outg_mod.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/outg_mod.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -37,4 +37,6 @@
   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 d2a5a833066965c4c4b33df109fd4e461f7275bb)
+++ src/outgrid_init.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -268,4 +268,11 @@
     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 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/par_mod.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -30,7 +30,4 @@
 !        Update 15 August 2013 IP                                              *
 !                                                                              *
-!        ESO 2016:                                                             *
-!          GFS specific parameters moved to gfs_mod.f90                        *
-!          ECMWF specific parameters moved to ecmwf_mod.f90                    *
 !                                                                              *
 !*******************************************************************************
@@ -80,4 +77,7 @@
   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,4 +89,6 @@
   ! 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
@@ -290,4 +292,6 @@
   integer,parameter :: unitboundcond=89
   integer,parameter :: unittmp=101
+! RLT
+  integer,parameter :: unitoutfactor=102
 
 !******************************************************
Index: src/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/readcommand.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -114,4 +114,5 @@
   surf_only, &
   cblflag, &
+  linversionout, &
   ohfields_path
 
@@ -147,4 +148,5 @@
   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/"
 
@@ -414,5 +416,5 @@
   !**********************************************************************
 
-  if ((iout.lt.1).or.(iout.gt.6)) then
+  if ((iout.lt.1).or.(iout.gt.5)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
     write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR        #### '
@@ -474,4 +476,14 @@
       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
Index: src/readpaths.f90
===================================================================
--- src/readpaths.f90	(revision 62e65c73ae7e3d8c4b8b884065ede76e07211ec8)
+++ src/readpaths.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine readpaths !(pathfile)
+subroutine readpaths
 
   !*****************************************************************************
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/timemanager.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -408,5 +408,14 @@
 #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 '
@@ -422,5 +431,9 @@
                 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
@@ -739,5 +752,11 @@
   if (ipout.eq.2) call partoutput(itime)     ! dump particle positions
 
-  if (linit_cond.ge.1) call initial_cond_output(itime)   ! dump initial cond. field
+  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
 
   !close(104)
Index: src/verttransform_ecmwf.f90
===================================================================
--- src/verttransform_ecmwf.f90	(revision 0a94e13d0d122caab219884890c29edcb2893c60)
+++ src/verttransform_ecmwf.f90	(revision 2eefa58ac797f119acc28fade41beff352006040)
@@ -81,4 +81,6 @@
   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
@@ -218,5 +220,4 @@
 !*************************
 
-
   do jy=0,nymin1
     do ix=0,nxmin1
@@ -229,4 +230,6 @@
   wzlev(:,:,1)=0.
   rhoh(:,:,1)=pold/(r_air*tvold)
+  ! RLT add pressure
+  prsh(:,:,1)=ps(:,:,1,n)
 
 
@@ -236,4 +239,6 @@
   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)
@@ -287,4 +292,6 @@
   pv(:,:,1,n)=pvh(:,:,1)
   rho(:,:,1,n)=rhoh(:,:,1)
+! RLT add pressure
+  prs(:,:,1,n)=prsh(:,:,1)
 
   uu(:,:,nz,n)=uuh(:,:,nuvz)
@@ -300,5 +307,6 @@
   pv(:,:,nz,n)=pvh(:,:,nuvz)
   rho(:,:,nz,n)=rhoh(:,:,nuvz)
-
+! RLT
+  prs(:,:,nz,n)=prsh(:,:,nuvz)
 
   kmin=2
@@ -320,4 +328,6 @@
           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
@@ -353,4 +363,6 @@
           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
