Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 07c3e71fd3f385404de511f95c658656679b9e29)
+++ src/com_mod.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -25,5 +25,5 @@
   !****************************************************************
 
-  character :: path(numpath+2*maxnests)*120
+  character :: path(numpath+2*maxnests)*200
   integer :: length(numpath+2*maxnests)
   character(len=256) :: pathfile, flexversion, flexversion_major, arg1, arg2
Index: src/concoutput_mpi.f90
===================================================================
--- src/concoutput_mpi.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/concoutput_mpi.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -59,4 +59,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),
@@ -93,4 +96,5 @@
   character(LEN=8),save :: file_stat='REPLACE'
   logical :: ldates_file
+  logical :: lexist
   integer :: ierr
   character(LEN=100) :: dates_char
@@ -194,4 +198,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
@@ -205,6 +212,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
@@ -259,9 +272,7 @@
         write(unitoutgrid) itime
       endif
-
       if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
         open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
              atime//'_'//anspec,form='unformatted')
-
         write(unitoutgridppt) itime
       endif
@@ -607,4 +618,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/ &
@@ -633,5 +687,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_nest.f90
===================================================================
--- src/concoutput_nest.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/concoutput_nest.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -588,5 +588,4 @@
   close(unitoutfactor)
 
-
   ! Reinitialization of grid
   !*************************
Index: src/concoutput_nest_mpi.f90
===================================================================
--- src/concoutput_nest_mpi.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/concoutput_nest_mpi.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -39,5 +39,4 @@
   !                                                                            *
   !*****************************************************************************
-
 
   use unc_mod
@@ -56,4 +55,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),
@@ -82,4 +84,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
@@ -158,4 +161,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
@@ -169,6 +175,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
@@ -193,4 +205,14 @@
 
   write(anspec,'(i3.3)') ks
+ 
+  if (DRYBKDEP.or.WETBKDEP) then !scavdep output
+      if (DRYBKDEP) &
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_nest_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      if (WETBKDEP) &
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_nest_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      write(unitoutgrid) itime
+  else
   if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     if (ldirect.eq.1) then
@@ -205,4 +227,5 @@
      write(unitoutgrid) itime
    endif
+  endif
 
   if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
@@ -535,28 +558,70 @@
   end do
 
-
-
-  ! 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
+! 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)
+
+  ! 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
+  creceptor(:,:)=0.
+  griduncn(:,:,:,:,:,:,:)=0.
+  
   if (mp_measure_time) call mpif_mtime('iotime',1)
   ! if (mp_measure_time) then
Index: src/concoutput_surf_nest_mpi.f90
===================================================================
--- src/concoutput_surf_nest_mpi.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/concoutput_surf_nest_mpi.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -633,14 +633,15 @@
 
   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.
+    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
@@ -648,5 +649,4 @@
       end do
     end do
-  end do
   end do
 
Index: src/gributils/class_gribfile_mod.f90
===================================================================
--- src/gributils/class_gribfile_mod.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/gributils/class_gribfile_mod.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -78,5 +78,5 @@
         ! module
 
-        USE grib_api
+        use grib_api 
 
         IMPLICIT NONE
Index: src/gridcheck_ecmwf.f90
===================================================================
--- src/gridcheck_ecmwf.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/gridcheck_ecmwf.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -55,5 +55,5 @@
   !**********************************************************************
 
-  use grib_api
+  use eccodes
   use par_mod
   use com_mod
Index: src/gridcheck_gfs.f90
===================================================================
--- src/gridcheck_gfs.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/gridcheck_gfs.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -56,5 +56,5 @@
   !**********************************************************************
 
-  use grib_api
+  use eccodes
   use par_mod
   use com_mod
Index: src/gridcheck_nests.f90
===================================================================
--- src/gridcheck_nests.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/gridcheck_nests.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -19,5 +19,5 @@
   !*****************************************************************************
 
-  use grib_api
+  use eccodes
   use par_mod
   use com_mod
Index: src/initial_cond_output_inversion_mpi.f90
===================================================================
--- src/initial_cond_output_inversion_mpi.f90	(revision 4ab2fbf0313a4a05b598f4d6976c1fb34e8f0f18)
+++ src/initial_cond_output_inversion_mpi.f90	(revision 4ab2fbf0313a4a05b598f4d6976c1fb34e8f0f18)
@@ -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_cond0(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_cond0(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/initial_cond_output_mpi.f90
===================================================================
--- src/initial_cond_output_mpi.f90	(revision 4ab2fbf0313a4a05b598f4d6976c1fb34e8f0f18)
+++ src/initial_cond_output_mpi.f90	(revision 4ab2fbf0313a4a05b598f4d6976c1fb34e8f0f18)
@@ -0,0 +1,151 @@
+!**********************************************************************
+! 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(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
+  real :: sp_fact,fact_recept
+  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+  logical :: sp_zer
+  character(len=3) :: anspec
+
+
+  !*********************************************************************
+  ! 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
+    open(97,file=path(2)(1:length(2))//'grid_initial'// &
+         '_'//anspec,form='unformatted')
+    write(97) itime
+
+    do kp=1,maxpointspec_act
+
+      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_cond0(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_cond0(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)
+
+
+    end do
+
+    close(97)
+
+  end do
+
+
+end subroutine initial_cond_output
Index: src/makefile
===================================================================
--- src/makefile	(revision c56da455adc4af0f5c28c039f48f7b6b73595ced)
+++ src/makefile	(revision 71f21286d3d1d3493b116ea707fc07f5072696df)
@@ -62,13 +62,14 @@
 	LIBPATH1 = ${ROOT_DIR}/lib
 else 
-# Compiled libraries under user ~flexpart, gfortran v5.4
-	ROOT_DIR = /homevip/flexpart/
-
+# Using system installed libraries at njord.nilu.no
 	F90	  = /usr/bin/gfortran
 	MPIF90    = /usr/bin/mpifort
-
-	INCPATH1  = ${ROOT_DIR}/gcc-5.4.0/include	
+	LIBPATH1 = /opt/eccodes/lib
+	INCPATH1  = /opt/eccodes/include
 	INCPATH2  = /usr/include
-	LIBPATH1 = ${ROOT_DIR}/gcc-5.4.0/lib
+
+
+# INCPATH1  = /usr/include	
+# INCPATH2  = /usr/include
 endif
 
@@ -88,10 +89,9 @@
 
 ## OPTIMIZATION LEVEL
-O_LEV = 0 # [0,1,2,3,g,s,fast]
+O_LEV = 2 # [0,1,2,3,g,s,fast]
 O_LEV_DBG = g # [0,g]
 
 ## LIBRARIES
-#LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff
-LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper $(NCOPT)
+LIBS = -leccodes -leccodes_f90 -lm  $(NCOPT)
 
 FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(NCOPT) $(FUSER)  #-Warray-bounds -fcheck=all # -march=native
@@ -128,8 +128,8 @@
 	redist.o		\
 	concoutput_surf.o	concoutput_surf_nest.o	\
-	concoutput_inversion_nest.o 	\
-	concoutput_inversion.o \
 	getfields.o \
-        readwind_ecmwf.o
+        readwind_ecmwf.o \
+	initial_cond_output.o \
+	initial_cond_output_inversion.o
 
 ## For MPI version 
@@ -145,9 +145,12 @@
 	concoutput_surf_mpi.o	concoutput_surf_nest_mpi.o	\
 	getfields_mpi.o \
-        readwind_ecmwf_mpi.o 
+        readwind_ecmwf_mpi.o \
+	initial_cond_output_mpi.o \
+	initial_cond_output_inversion_mpi.o
 
 OBJECTS_NCF = netcdf_output_mod.o
 
 OBJECTS = \
+initial_cond_calc.o \
 advance.o		initialize.o		\
 writeheader.o		writeheader_txt.o  	\
@@ -202,8 +205,10 @@
 drydepokernel_nest.o    zenithangle.o \
 ohreaction.o            getvdep_nests.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 \
-cbl.o
+cbl.o \
+concoutput_inversion_nest.o 	\
+concoutput_inversion.o \
+
 
 ifeq ($(ncf), yes)
@@ -325,6 +330,10 @@
 init_domainfill_mpi.o: com_mod.o mpi_mod.o par_mod.o point_mod.o random_mod.o
 initial_cond_calc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
+initial_cond_calc_mpi.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_mpi.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mpi_mod.o
 initial_cond_output_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o \
+	unc_mod.o
+initial_cond_output_inversion_mpi.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
Index: src/mpi_mod.f90
===================================================================
--- src/mpi_mod.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/mpi_mod.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -2579,4 +2579,56 @@
 601 end subroutine mpif_tm_reduce_grid_nest
 
+
+  subroutine mpif_tm_reduce_initcond
+!***********************************************************************
+! Collect init_cond to PID 0, adding from all processes.
+!
+!
+!***********************************************************************
+    use com_mod
+    use unc_mod
+    use par_mod
+
+    implicit none
+
+    integer :: grid_size
+
+!**********************************************************************
+
+    grid_size=numxgrid*numygrid*numzgrid*maxspec* &
+         & maxpointspec_act
+
+! Time for MPI communications
+    if (mp_measure_time) call mpif_mtime('commtime',0)
+
+
+#ifdef USE_MPIINPLACE
+! Using in-place reduction
+    if (lroot) then
+      call MPI_Reduce(MPI_IN_PLACE, init_cond, grid_size, mp_sp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    else
+      call MPI_Reduce(init_cond, 0, grid_size, mp_sp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    end if
+
+#else
+    call MPI_Reduce(init_cond, init_cond0, grid_size, mp_sp, MPI_SUM, id_root, &
+         & mp_comm_used, mp_ierr)
+    if (mp_ierr /= 0) goto 600
+!    if (lroot) init_cond = init_cond0
+#endif
+
+    if (mp_measure_time) call mpif_mtime('commtime',1)
+
+    goto 601
+
+600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
+    stop
+
+601 end subroutine mpif_tm_reduce_initcond
+    
 
   subroutine mpif_mtime(ident,imode)
Index: src/outgrid_init.f90
===================================================================
--- src/outgrid_init.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/outgrid_init.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -297,4 +297,12 @@
          maxpointspec_act),stat=stat)
     if (stat.ne.0) write(*,*)'ERROR: could not allocate init_cond'
+    if (mpi_mode.gt.0) then
+      if (lroot) then
+        allocate(init_cond0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
+             maxpointspec_act),stat=stat)
+      else
+        allocate(init_cond0(1,1,1,1,1))
+      end if
+    end if
   endif
 
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/par_mod.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -133,5 +133,5 @@
   ! ECMWF
 ! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=359 ! 1.0 degree 92 level
-!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0 ! 1.0 degree 138 level
+  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0 ! 1.0 degree 138 level
 !   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359 ! 1.0 degree 138 level
 ! integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359  ! 0.5 degree 138 level
@@ -139,6 +139,6 @@
 
 ! GFS
-   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
-   integer :: nxshift=0 ! shift not fixed for the executable 
+   ! integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
+   ! integer :: nxshift=0 ! shift not fixed for the executable 
 
 
@@ -205,6 +205,6 @@
   !**************************************************
 
-  integer,parameter :: maxpart=100000
-  integer,parameter :: maxspec=1
+  integer,parameter :: maxpart=10000000
+  integer,parameter :: maxspec=6
 
   real,parameter :: minmass=0.0001
Index: src/readwind_ecmwf.f90
===================================================================
--- src/readwind_ecmwf.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/readwind_ecmwf.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -48,5 +48,5 @@
 !**********************************************************************
 
-  use grib_api
+  use eccodes
   use par_mod
   use com_mod
Index: src/readwind_gfs.f90
===================================================================
--- src/readwind_gfs.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/readwind_gfs.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -47,5 +47,5 @@
   !***********************************************************************
 
-  use grib_api
+  use eccodes
   use par_mod
   use com_mod
Index: src/readwind_nests.f90
===================================================================
--- src/readwind_nests.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/readwind_nests.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
@@ -22,5 +22,5 @@
   !*****************************************************************************
 
-  use grib_api
+  use eccodes
   use par_mod
   use com_mod
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/timemanager.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -90,9 +90,10 @@
   implicit none
 
+  integer(selected_int_kind(16)) :: idummy,idummy2
   integer :: metdata_format
   integer :: j,ks,kp,l,n,itime=0,nstop,nstop1
 ! integer :: ksp
   integer :: loutnext,loutstart,loutend
-  integer :: ix,jy,ldeltat,itage,nage,idummy
+  integer :: ix,jy,ldeltat,itage,nage
   integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme 
   real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact,wetscav
@@ -379,9 +380,9 @@
               call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
             endif
-          else  
+          else ! surf only  
             if (verbosity.eq.1) then
-             print*,'call concoutput_surf '
-             call system_clock(count_clock)
-             write(*,*) 'system clock',count_clock - count_clock0   
+              print*,'call concoutput_surf '
+              call system_clock(count_clock)
+              write(*,*) 'system clock',count_clock - count_clock0   
             endif
             if (lnetcdfout.eq.1) then
@@ -398,5 +399,5 @@
                 endif
               else
-              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+                call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
               endif
               if (verbosity.eq.1) then
@@ -416,6 +417,6 @@
                   call concoutput_inversion_nest(itime,outnum)
                 else 
-                call concoutput_surf_nest(itime,outnum)
-              endif
+                  call concoutput_surf_nest(itime,outnum)
+                endif
               endif
             else
@@ -581,5 +582,5 @@
        do ks=1,nspec
          if  ((xscav_frac1(j,ks).lt.0)) then
-            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
+            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy2,wetscav)
             if (wetscav.gt.0) then
                 xscav_frac1(j,ks)=wetscav* &
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/timemanager_mpi.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -91,4 +91,5 @@
   implicit none
 
+  integer(selected_int_kind(16)) :: idummy,idummy2
   integer :: metdata_format
   logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
@@ -97,5 +98,5 @@
   integer :: ip,irec
   integer :: loutnext,loutstart,loutend
-  integer :: ix,jy,ldeltat,itage,nage,idummy
+  integer :: ix,jy,ldeltat,itage,nage
   integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme 
   integer :: numpart_tot_mpi ! for summing particles on all processes
@@ -159,5 +160,5 @@
 !********************************************************************
 
-    if (mp_dbg_mode) write(*,*) 'myid, itime: ',mp_pid,itime
+!    if (mp_dbg_mode) write(*,*) 'myid, itime: ',mp_pid,itime
     
     if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
@@ -280,5 +281,5 @@
 
     if (lmpreader.and.lmp_use_reader) then
-      if (itime.lt.ideltas*ldirect) then
+      if (abs(itime).lt.ideltas*ldirect) then
         cycle
       else
@@ -480,5 +481,5 @@
               creceptor(:,:)=0.
             end if
-          else 
+          else ! surf only 
             if (lroot) then
               if (lnetcdfout.eq.1) then
@@ -488,5 +489,9 @@
 #endif
               else
-                call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+                if (linversionout.eq.1) then
+                  call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+                else
+                  call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+                end if
               end if
             else
@@ -516,5 +521,17 @@
 
               else
-                call concoutput_surf_nest(itime,outnum)
+                if(linversionout.eq.1) then
+                  if (lroot) then
+                    call concoutput_inversion_nest(itime,outnum)
+                  else
+                    griduncn(:,:,:,:,:,:,:)=0.
+                  end if
+                else
+                  if (lroot) then
+                    call concoutput_surf_nest(itime,outnum)
+                  else
+                    griduncn(:,:,:,:,:,:,:)=0.
+                  end if
+                end if
               end if
             else
@@ -724,5 +741,5 @@
        do ks=1,nspec
          if  ((xscav_frac1(j,ks).lt.0)) then
-            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
+            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy2,wetscav)
             if (wetscav.gt.0) then
                 xscav_frac1(j,ks)=wetscav* &
@@ -877,13 +894,11 @@
 !*****************************************************************************
 
-! eso :TODO: this not implemented yet (transfer particles to PID 0 or rewrite)
-! the tools to do this are already in mpi_mod.f90
-  if (lroot) then 
-    do j=1,numpart
-      if (linit_cond.ge.1) call initial_cond_calc(itime,j)
-    end do
-  end if
-
-
+  do j=1,numpart
+    if (linit_cond.ge.1) call initial_cond_calc(itime,j)
+  end do
+
+! Transfer sum of init_cond field to root process, for output
+  call mpif_tm_reduce_initcond
+    
   if (ipout.eq.2) then
 ! MPI process 0 creates the file, the other processes append to it
@@ -897,6 +912,12 @@
   end if
 
-! eso :TODO: MPI
-  if (linit_cond.ge.1.and.lroot) call initial_cond_output(itime)   ! dump initial cond. field
+  
+  if (linit_cond.ge.1.and.lroot) 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
 
 
Index: src/unc_mod.f90
===================================================================
--- src/unc_mod.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/unc_mod.f90	(revision f3054ea95b97388f57786fce68d9556e51c62ff2)
@@ -35,4 +35,5 @@
 
   real,allocatable, dimension (:,:,:,:,:) :: init_cond
+  real,allocatable, dimension (:,:,:,:,:) :: init_cond0 ! for MPI
 
 end module unc_mod
