Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/FLEXPART.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -49,4 +49,6 @@
   integer :: i,j,ix,jy,inest
   integer :: idummy = -320
+  character(len=256) :: inline_options  !pathfile, flexversion, arg2
+
 
   ! Generate a large number of random numbers
@@ -58,25 +60,106 @@
   call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
 
+  ! 
+  flexversion='Version 9.1.8  (2013-12-08)'
+  !verbosity=0
+  ! Read the pathnames where input/output files are stored
+  !*******************************************************
+
+  inline_options='none'
+  select case (iargc())
+  case (2)
+    call getarg(1,arg1)
+    pathfile=arg1
+    call getarg(2,arg2)
+    inline_options=arg2
+  case (1)
+    call getarg(1,arg1)
+    pathfile=arg1
+    verbosity=0
+    if (arg1(1:1).eq.'-') then
+        write(pathfile,'(a11)') './pathnames'
+        inline_options=arg1 
+    endif
+  case (0)
+    write(pathfile,'(a11)') './pathnames'
+    verbosity=0
+  end select
+  
+    if (inline_options(1:1).eq.'-') then
+      print*, 'inline options=', inline_options       
+      if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
+         print*, 'verbose mode 1: additional information will be displayed'
+         verbosity=1
+      endif
+      if (trim(inline_options).eq.'-v2') then
+         print*, 'verbose mode 2: additional information will be displayed'
+         verbosity=2
+      endif
+      if (trim(inline_options).eq.'-i') then
+         print*, 'info mode: will provide run specific information and stop'
+         verbosity=1
+         info_flag=1
+      endif
+      if (trim(inline_options).eq.'-i2') then
+         print*, 'info mode: will provide run specific information and stop'
+         verbosity=2
+         info_flag=1
+      endif
+    endif
+
+
   ! Print the GPL License statement
   !*******************************************************
-  print*,'Welcome to FLEXPART Version 9.0'
+  print*,'Welcome to FLEXPART', trim(flexversion)
   print*,'FLEXPART is free software released under the GNU Genera'// &
        'l Public License.'
-
-  ! Read the pathnames where input/output files are stored
-  !*******************************************************
-
-  call readpaths
+            
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readpaths'
+  endif 
+  call readpaths(pathfile)
+   
+ 
+  if (verbosity.gt.1) then !show clock info 
+     !print*,'length(4)',length(4)
+     !count=0,count_rate=1000
+     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
+     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
+     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
+     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
+     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
+  endif
+
 
   ! Read the user specifications for the current model run
   !*******************************************************
 
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readcommand'
+  endif
   call readcommand
+  if (verbosity.gt.0) then
+        WRITE(*,*) '    ldirect=', ldirect
+        WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
+        WRITE(*,*) '    iedate,ietime=', iedate,ietime
+        if (verbosity.gt.1) then   
+                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+        endif     
+  endif
 
 
   ! Read the age classes to be used
   !********************************
-
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readageclasses'
+  endif
   call readageclasses
+
+  if (verbosity.gt.1) then   
+                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+  endif     
+ 
 
 
@@ -84,12 +167,28 @@
   !******************************************************************
 
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readavailable'
+  endif  
   call readavailable
-
 
   ! Read the model grid specifications,
   ! both for the mother domain and eventual nests
   !**********************************************
+ 
+  if (verbosity.gt.0) then
+     WRITE(*,*) 'call gridcheck'
+  endif
 
   call gridcheck
+
+  if (verbosity.gt.1) then   
+     CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+     WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+  endif      
+  
+
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call gridcheck_nests'
+  endif  
   call gridcheck_nests
 
@@ -98,13 +197,24 @@
   !************************************
 
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readoutgrid'
+  endif
+
   call readoutgrid
-  if (nested_output.eq.1) call readoutgrid_nest
-
+
+  if (nested_output.eq.1) then
+          call readoutgrid_nest
+    if (verbosity.gt.0) then
+        WRITE(*,*) '# readoutgrid_nest'
+    endif
+  endif
 
   ! Read the receptor points for which extra concentrations are to be calculated
   !*****************************************************************************
 
+  if (verbosity.eq.1) then
+     print*,'call readreceptors'
+  endif
   call readreceptors
-
 
   ! Read the physico-chemical species property table
@@ -117,4 +227,7 @@
   !***************************
 
+  if (verbosity.gt.0) then
+    print*,'call readlanduse'
+  endif
   call readlanduse
 
@@ -123,4 +236,7 @@
   !********************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'call assignland'
+  endif
   call assignland
 
@@ -130,16 +246,25 @@
   !**********************************************
 
+  if (verbosity.gt.0) then
+    print*,'call readreleases'
+  endif
   call readreleases
+
 
   ! Read and compute surface resistances to dry deposition of gases
   !****************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'call readdepo'
+  endif
   call readdepo
-
 
   ! Convert the release point coordinates from geografical to grid coordinates
   !***************************************************************************
 
-  call coordtrafo
+  call coordtrafo  
+  if (verbosity.gt.0) then
+    print*,'call coordtrafo'
+  endif
 
 
@@ -147,4 +272,7 @@
   !*****************************************
 
+  if (verbosity.gt.0) then
+    print*,'Initialize all particles to non-existent'
+  endif
   do j=1,maxpart
     itra1(j)=-999999999
@@ -155,6 +283,12 @@
 
   if (ipin.eq.1) then
+    if (verbosity.gt.0) then
+      print*,'call readpartpositions'
+    endif
     call readpartpositions
   else
+    if (verbosity.gt.0) then
+      print*,'numpart=0, numparticlecount=0'
+    endif    
     numpart=0
     numparticlecount=0
@@ -166,4 +300,8 @@
   !***************************************************************
 
+
+  if (verbosity.gt.0) then
+    print*,'call outgrid_init'
+  endif
   call outgrid_init
   if (nested_output.eq.1) call outgrid_init_nest
@@ -173,6 +311,10 @@
   !******************
 
-  if (OHREA.eqv..TRUE.) &
+  if (OHREA.eqv..TRUE.) then
+    if (verbosity.gt.0) then
+       print*,'call readOHfield'
+    endif
        call readOHfield
+  endif
 
   ! Write basic information on the simulation to a file "header"
@@ -180,7 +322,25 @@
   !******************************************************************
 
+
+  if (verbosity.gt.0) then
+    print*,'call writeheader'
+  endif
+
   call writeheader
-  if (nested_output.eq.1) call writeheader_nest
-  open(unitdates,file=path(2)(1:length(2))//'dates')
+  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
+  call writeheader_txt
+  !if (nested_output.eq.1) call writeheader_nest
+  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
+
+  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
+  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
+
+
+
+  !open(unitdates,file=path(2)(1:length(2))//'dates')
+
+  if (verbosity.gt.0) then
+    print*,'call openreceptors'
+  endif
   call openreceptors
   if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
@@ -190,4 +350,7 @@
   !***************************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'discretize release times'
+  endif
   do i=1,numpoint
     ireleasestart(i)=nint(real(ireleasestart(i))/ &
@@ -200,4 +363,8 @@
   ! Initialize cloud-base mass fluxes for the convection scheme
   !************************************************************
+
+  if (verbosity.gt.0) then
+    print*,'Initialize cloud-base mass fluxes for the convection scheme'
+  endif
 
   do jy=0,nymin1
@@ -218,4 +385,16 @@
   !********************************
 
+  if (verbosity.gt.0) then
+     if (verbosity.gt.1) then   
+       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+       WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+     endif
+     if (info_flag.eq.1) then
+         print*, 'info only mode (stop)'    
+         stop
+     endif
+     print*,'call timemanager'
+  endif
+
   call timemanager
 
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/com_mod.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -7,5 +7,5 @@
 !        June 1996                                                             *
 !                                                                              *
-!        Last update: 9 August 2000                                            *
+!        Last update:15 August 2013 IP                                         *
 !                                                                              *
 !*******************************************************************************
@@ -25,8 +25,11 @@
   character :: path(numpath+2*maxnests)*120
   integer :: length(numpath+2*maxnests)
-
+  character(len=256) :: pathfile, flexversion, arg1, arg2
+  
   ! path                    path names needed for trajectory model
   ! length                  length of path names needed for trajectory model
-
+  ! pathfile                file where pathnames are stored
+  ! flexversion             version of flexpart
+  ! arg                     input arguments from launch at command line
 
   !********************************************************
@@ -65,5 +68,5 @@
   integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
   integer :: mquasilag,nested_output,ind_source,ind_receptor
-  integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond
+  integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
   logical :: turbswitch
 
@@ -91,4 +94,6 @@
   ! mquasilag 0: normal run
   !      1: Particle position output is produced in a condensed format and particles are numbered
+  ! surf_only   switch output in grid_time files for surface only or full vertical resolution
+  !      0=no (full vertical resolution), 1=yes (surface only)
   ! nested_output: 0 no, 1 yes
   ! turbswitch              determines how the Markov chain is formulated
@@ -139,4 +144,6 @@
   real :: decay(maxspec)
   real :: weta(maxspec),wetb(maxspec)
+! NIK: 31.01.2013- parameters for in-cloud scavening
+  real :: weta_in(maxspec), wetb_in(maxspec), wetc_in(maxspec), wetd_in(maxspec)
   real :: reldiff(maxspec),henry(maxspec),f0(maxspec)
   real :: density(maxspec),dquer(maxspec),dsigma(maxspec)
@@ -172,5 +179,7 @@
 
   ! WET DEPOSITION
-  ! weta, wetb              parameters for determining wet scavenging coefficients
+  ! weta, wetb              parameters for determining below-cloud wet scavenging coefficients
+  ! weta_in, wetb_in       parameters for determining in-cloud wet scavenging coefficients
+  ! wetc_in, wetd_in       parameters for determining in-cloud wet scavenging coefficients
 
   ! GAS DEPOSITION
@@ -331,6 +340,10 @@
   real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
   real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
+  !scavenging NIK, PS
   integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
   integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
+      integer icloudbot(0:nxmax-1,0:nymax-1,2)
+      integer icloudthck(0:nxmax-1,0:nymax-1,2)
+
 
   ! uu,vv,ww [m/2]       wind components in x,y and z direction
@@ -346,4 +359,8 @@
   !      rainout  conv/lsp dominated  2/3
   !      washout  conv/lsp dominated  4/5
+! PS 2013
+!c icloudbot (m)        cloud bottom height
+!c icloudthck (m)       cloud thickness     
+
   ! pplev for the GFS version
 
@@ -662,5 +679,4 @@
 
 
-
   !********************
   ! Random number field
@@ -671,4 +687,11 @@
   ! rannumb                 field of normally distributed random numbers
 
+  !********************
+  ! Verbosity, testing flags
+  !********************   
+  integer :: verbosity=0
+  integer :: info_flag=0
+  INTEGER :: count_clock, count_clock0,  count_rate, count_max
+   
 
 end module com_mod
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/concoutput.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -108,6 +108,7 @@
   write(adate,'(i8.8)') jjjjmmdd
   write(atime,'(i6.6)') ihmmss
+  open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
   write(unitdates,'(a)') adate//atime
-
+  close(unitdates)  
 
   ! For forward simulations, output fields have dimension MAXSPEC,
@@ -158,5 +159,5 @@
         yl=outlat0+real(jy)*dyout
         xl=(xl-xlon0)/dx
-        yl=(yl-ylat0)/dx
+        yl=(yl-ylat0)/dy !v9.1.1 
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)
Index: src/concoutput_surf.f90
===================================================================
--- src/concoutput_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/concoutput_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -0,0 +1,744 @@
+!**********************************************************************
+! 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_surf(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                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! 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
+
+  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
+
+  !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 :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
+  real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
+  real :: 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
+
+
+  if (verbosity.eq.1) then
+     print*,'inside concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+
+  ! 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
+  !write(unitdates,'(a)') adate//atime
+
+    open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
+      write(unitdates,'(a)') adate//atime
+    close(unitdates)
+
+  ! 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_surf 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)/dx
+        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
+      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)
+    end do
+
+
+  ! 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_surf 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
+  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+    if (ldirect.eq.1) then
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+    else
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+    endif
+     write(unitoutgrid) itime
+   endif
+
+  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
+   open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
+        atime//'_'//anspec,form='unformatted')
+
+    write(unitoutgridppt) itime
+  endif
+
+  do kp=1,maxpointspec_act
+  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_surf 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_surf (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
+  !oncentraion 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) 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) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+
+  if (verbosity.eq.1) then
+     print*,'concoutput_surf (Concentrations)'
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+  ! Concentrations
+
+  ! if surf_only write only 1st layer 
+
+         if(surf_only.eq.1) then
+         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)
+         else
+
+  ! write full vertical resolution
+  if (verbosity.eq.1) then
+     print*,'concoutput_surf (write full vertical resolution)'
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  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 (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 ! surf_only
+
+    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) 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) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+  ! Mixing ratios
+
+  ! if surf_only write only 1st layer 
+
+         if(surf_only.eq.1) then
+         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)
+         else
+
+  ! write full vertical resolution
+
+         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 (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 ! surf_only
+
+      endif ! output for ppt
+
+  end do
+  end do
+
+    close(unitoutgridppt)
+    close(unitoutgrid)
+
+  end do
+
+  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
+
+
+
+  ! 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_surf
Index: src/concoutput_surf_nest.f90
===================================================================
--- src/concoutput_surf_nest.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/concoutput_surf_nest.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -0,0 +1,652 @@
+!**********************************************************************
+! 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_surf_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                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! 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
+
+  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
+
+  !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 :: 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
+  character :: adate*8,atime*6
+  character(len=3) :: anspec
+
+
+  ! 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
+
+
+  ! 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
+      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)
+    end do
+
+
+  ! 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
+  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+    if (ldirect.eq.1) then
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_' &
+           //adate// &
+           atime//'_'//anspec,form='unformatted')
+    else
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_' &
+           //adate// &
+           atime//'_'//anspec,form='unformatted')
+    endif
+     write(unitoutgrid) itime
+   endif
+
+  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
+   open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_' &
+        //adate// &
+        atime//'_'//anspec,form='unformatted')
+
+    write(unitoutgridppt) itime
+  endif
+
+  do kp=1,maxpointspec_act
+  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
+  !oncentraion 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
+
+  ! if surf_only write only 1st layer 
+
+         if(surf_only.eq.1) then
+         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)
+         else
+
+  ! write full vertical resolution
+
+         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 (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 ! surf_only
+
+
+    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
+
+    ! if surf_only write only 1st layer 
+
+         if(surf_only.eq.1) then
+         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)
+         else
+
+  ! write full vertical resolution
+
+         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 (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 ! surf_only
+
+      endif ! output for ppt
+
+  end do
+  end do
+
+    close(unitoutgridppt)
+    close(unitoutgrid)
+
+  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
+    end do
+  end do
+  end do
+
+
+end subroutine concoutput_surf_nest
+
Index: src/gridcheck.f90
===================================================================
--- src/gridcheck.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/gridcheck.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -98,5 +98,6 @@
 
   integer :: isec1(56),isec2(22+nxmax+nymax)
-  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+  !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+  real(kind=4) :: zsec2(184),zsec4(jpunp)
   character(len=1) :: opt
 
@@ -466,5 +467,4 @@
     k=nlev_ec+1+numskip+i
     akm(nwz-i+1)=zsec2(j)
-  !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
     bkm(nwz-i+1)=zsec2(k)
   end do
Index: src/gridcheck_fnl.f90
===================================================================
--- src/gridcheck_fnl.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/gridcheck_fnl.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -0,0 +1,545 @@
+!**********************************************************************
+! 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 gridcheck
+
+  !**********************************************************************
+  !                                                                     *
+  !             FLEXPART MODEL SUBROUTINE GRIDCHECK                     *
+  !                                                                     *
+  !**********************************************************************
+  !                                                                     *
+  !             AUTHOR:      G. WOTAWA                                  *
+  !             DATE:        1997-08-06                                 *
+  !             LAST UPDATE: 1997-10-10                                 *
+  !                                                                     *
+  !             Update:      1999-02-08, global fields allowed, A. Stohl*
+  !             CHANGE: 17/11/2005, Caroline Forster, GFS data          *
+  !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
+  !                                 ECMWF grib_api                      *
+  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
+  !                                 ECMWF grib_api                      *
+  !                                                                     *
+  !**********************************************************************
+  !                                                                     *
+  ! DESCRIPTION:                                                        *
+  !                                                                     *
+  ! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT      *
+  ! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST-   *
+  ! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE       *
+  ! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES    *
+  ! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT    *
+  ! ANY CALL.                                                           *
+  !                                                                     *
+  ! XLON0                geographical longitude of lower left gridpoint *
+  ! YLAT0                geographical latitude of lower left gridpoint  *
+  ! NX                   number of grid points x-direction              *
+  ! NY                   number of grid points y-direction              *
+  ! DX                   grid distance x-direction                      *
+  ! DY                   grid distance y-direction                      *
+  ! NUVZ                 number of grid points for horizontal wind      *
+  !                      components in z direction                      *
+  ! NWZ                  number of grid points for vertical wind        *
+  !                      component in z direction                       *
+  ! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
+  !                      points of the polar stereographic grid):       *
+  !                      used to check the CFL criterion                *
+  ! VHEIGHT(1)-         heights of gridpoints where u and v are        *
+  ! UVHEIGHT(NUVZ)       given                                          *
+  ! WHEIGHT(1)-          heights of gridpoints where w is given         *
+  ! WHEIGHT(NWZ)                                                        *
+  !                                                                     *
+  !**********************************************************************
+
+  use grib_api
+  use par_mod
+  use com_mod
+  use conv_mod
+  use cmapf_mod, only: stlmbr,stcm2p
+
+  implicit none
+
+  !HSO  parameters for grib_api
+  integer :: ifile
+  integer :: iret
+  integer :: igrib
+  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
+  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
+  !HSO  end
+  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
+  real :: sizesouth,sizenorth,xauxa,pint
+  real :: akm_usort(nwzmax)
+  real,parameter :: eps=0.0001
+
+  ! NCEP GFS
+  real :: pres(nwzmax), help
+
+  integer :: i179,i180,i181
+
+  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
+
+  integer :: isec1(8),isec2(3)
+  real(kind=4) :: zsec4(jpunp)
+  character(len=1) :: opt
+
+  !HSO  grib api error messages
+  character(len=24) :: gribErrorMsg = 'Error reading grib file'
+  character(len=20) :: gribFunction = 'gridcheckwind_gfs'
+  !
+  if (numbnests.ge.1) then
+  write(*,*) ' ###########################################'
+  write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
+  write(*,*) ' NO NESTED WINDFIELDAS ALLOWED FOR GFS!      '
+  write(*,*) ' ###########################################'
+  stop
+  endif
+
+  iumax=0
+  iwmax=0
+
+  if(ideltas.gt.0) then
+    ifn=1
+  else
+    ifn=numbwf
+  endif
+  !
+  ! OPENING OF DATA FILE (GRIB CODE)
+  !
+5   call grib_open_file(ifile,path(3)(1:length(3)) &
+         //trim(wfname(ifn)),'r',iret)
+  if (iret.ne.GRIB_SUCCESS) then
+    goto 999   ! ERROR DETECTED
+  endif
+  !turn on support for multi fields messages
+  call grib_multi_support_on
+
+  ifield=0
+10   ifield=ifield+1
+  !
+  ! GET NEXT FIELDS
+  !
+  call grib_new_from_file(ifile,igrib,iret)
+  if (iret.eq.GRIB_END_OF_FILE )  then
+    goto 30    ! EOF DETECTED
+  elseif (iret.ne.GRIB_SUCCESS) then
+    goto 999   ! ERROR DETECTED
+  endif
+
+  !first see if we read GRIB1 or GRIB2
+  call grib_get_int(igrib,'editionNumber',gribVer,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  if (gribVer.eq.1) then ! GRIB Edition 1
+
+  !read the grib1 identifiers
+  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'level',isec1(8),iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  !get the size and data of the values array
+  call grib_get_real4_array(igrib,'values',zsec4,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  else ! GRIB Edition 2
+
+  !read the grib2 identifiers
+  call grib_get_int(igrib,'discipline',discipl,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'parameterCategory',parCat,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'parameterNumber',parNum,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
+       valSurf,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  !convert to grib1 identifiers
+  isec1(6)=-1
+  isec1(7)=-1
+  isec1(8)=-1
+  if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
+    isec1(6)=33          ! indicatorOfParameter
+    isec1(7)=100         ! indicatorOfTypeOfLevel
+    isec1(8)=valSurf/100 ! level, convert to hPa
+  elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
+    isec1(6)=7           ! indicatorOfParameter
+    isec1(7)=1           ! indicatorOfTypeOfLevel
+    isec1(8)=0
+  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
+       .and.(discipl.eq.2)) then ! LSM
+    isec1(6)=81          ! indicatorOfParameter
+    isec1(7)=1           ! indicatorOfTypeOfLevel
+    isec1(8)=0
+  endif
+
+  if (isec1(6).ne.-1) then
+  !  get the size and data of the values array
+    call grib_get_real4_array(igrib,'values',zsec4,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+  endif
+
+  endif ! gribVer
+
+  if(ifield.eq.1) then
+
+  !get the required fields from section 2
+  !store compatible to gribex input
+  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
+       isec2(2),iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
+       isec2(3),iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+       xaux1in,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
+       xaux2in,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+       yaux1in,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
+       yaux2in,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  xaux1=xaux1in
+!  xaux2=xaux2in
+  xaux2=xaux2in+360.    !!! Transform FNL xauu2in from -1 to 359 (Xuekue Fang, 31 Jan 2013)
+  yaux1=yaux1in
+  yaux2=yaux2in
+  write(*,*) 'xaux1,xaux2,yaux1,yaux2',xaux1,xaux2,yaux1,yaux2
+  nxfield=isec2(2)
+  ny=isec2(3)
+
+  
+
+  if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
+    xaux1=-179.0                             ! 359 DEG EAST ->
+    xaux2=-179.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
+  endif                                      ! TO 180 DEG EAST
+  if (xaux1.gt.180) xaux1=xaux1-360.0
+  if (xaux2.gt.180) xaux2=xaux2-360.0
+  if (xaux1.lt.-180) xaux1=xaux1+360.0
+  if (xaux2.lt.-180) xaux2=xaux2+360.0
+  if (xaux2.lt.xaux1) xaux2=xaux2+360.
+  xlon0=xaux1
+  ylat0=yaux1
+  write(*,*) 'xlon0,ylat0',xlon0,ylat0
+
+  dx=(xaux2-xaux1)/real(nxfield-1)
+  dy=(yaux2-yaux1)/real(ny-1)
+  dxconst=180./(dx*r_earth*pi)
+  dyconst=180./(dy*r_earth*pi)
+  write(*,*) 'dx,dy,dxconst,dyconst',dx,dy,dxconst,dyconst
+  !HSO end edits
+
+
+  ! Check whether fields are global
+  ! If they contain the poles, specify polar stereographic map
+  ! projections using the stlmbr- and stcm2p-calls
+  !***********************************************************
+
+    xauxa=abs(xaux2+dx-360.-xaux1)
+    if (xauxa.lt.0.001) then
+      nx=nxfield+1                 ! field is cyclic
+      xglobal=.true.
+      if (abs(nxshift).ge.nx) &
+           stop 'nxshift in file par_mod is too large'
+      xlon0=xlon0+real(nxshift)*dx
+    else
+      nx=nxfield
+      xglobal=.false.
+      if (nxshift.ne.0) &
+           stop 'nxshift (par_mod) must be zero for non-global domain'
+    endif
+    nxmin1=nx-1
+    nymin1=ny-1
+    if (xlon0.gt.180.) xlon0=xlon0-360.
+    xauxa=abs(yaux1+90.)
+    if (xglobal.and.xauxa.lt.0.001) then
+      sglobal=.true.               ! field contains south pole
+  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
+  ! map scale
+      sizesouth=6.*(switchsouth+90.)/dy
+      call stlmbr(southpolemap,-90.,0.)
+      call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
+           sizesouth,switchsouth,180.)
+      switchsouthg=(switchsouth-ylat0)/dy
+    else
+      sglobal=.false.
+      switchsouthg=999999.
+    endif
+    xauxa=abs(yaux2-90.)
+    if (xglobal.and.xauxa.lt.0.001) then
+      nglobal=.true.               ! field contains north pole
+  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
+  ! map scale
+      sizenorth=6.*(90.-switchnorth)/dy
+      call stlmbr(northpolemap,90.,0.)
+      call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
+           sizenorth,switchnorth,180.)
+      switchnorthg=(switchnorth-ylat0)/dy
+    else
+      nglobal=.false.
+      switchnorthg=999999.
+    endif
+  endif ! ifield.eq.1
+
+  if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
+  if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
+
+  ! NCEP ISOBARIC LEVELS
+  !*********************
+
+  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
+    iumax=iumax+1
+    pres(iumax)=real(isec1(8))*100.0
+  endif
+
+
+  i179=nint(179./dx)
+  if (dx.lt.0.7) then
+    i180=nint(180./dx)+1    ! 0.5 deg data
+  else
+    i180=nint(179./dx)+1    ! 1 deg data
+  endif
+  i181=i180+1
+
+
+  ! NCEP TERRAIN
+  !*************
+
+  if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
+    do jy=0,ny-1
+      do ix=0,nxfield-1
+        help=zsec4(nxfield*(ny-jy-1)+ix+1)
+        if(ix.le.i180) then
+          oro(i179+ix,jy)=help
+          excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+        else
+          oro(ix-i181,jy)=help
+          excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+        endif
+      end do
+    end do
+  endif
+
+  ! NCEP LAND SEA MASK
+  !*******************
+
+  if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
+    do jy=0,ny-1
+      do ix=0,nxfield-1
+        help=zsec4(nxfield*(ny-jy-1)+ix+1)
+        if(ix.le.i180) then
+          lsm(i179+ix,jy)=help
+        else
+          lsm(ix-i181,jy)=help
+        endif
+      end do
+    end do
+  endif
+
+  call grib_release(igrib)
+
+  goto 10                      !! READ NEXT LEVEL OR PARAMETER
+  !
+  ! CLOSING OF INPUT DATA FILE
+  !
+
+  ! HSO
+30   continue
+  call grib_close_file(ifile)
+  ! HSO end edits
+
+  nuvz=iumax
+  nwz =iumax
+  nlev_ec=iumax
+
+  if (nx.gt.nxmax) then
+   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) nx,nxmax
+    stop
+  endif
+
+  if (ny.gt.nymax) then
+   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) ny,nymax
+    stop
+  endif
+
+  if (nuvz.gt.nuvzmax) then
+    write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
+         'direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) nuvz,nuvzmax
+    stop
+  endif
+
+  if (nwz.gt.nwzmax) then
+    write(*,*) 'FLEXPART error: Too many w grid points in z '// &
+         'direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) nwz,nwzmax
+    stop
+  endif
+
+  ! If desired, shift all grids by nxshift grid cells
+  !**************************************************
+
+  if (xglobal) then
+    call shift_field_0(oro,nxfield,ny)
+    call shift_field_0(lsm,nxfield,ny)
+    call shift_field_0(excessoro,nxfield,ny)
+  endif
+
+  ! Output of grid info
+  !********************
+
+  write(*,*)
+  write(*,*)
+  write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
+       nuvz,nwz
+  write(*,*)
+  write(*,'(a)') 'other domain:'
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
+       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
+       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
+  write(*,*)
+
+
+  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
+  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
+
+  numskip=nlev_ec-nuvz  ! number of ecmwf model layers not used
+                        ! by trajectory model
+  do i=1,nwz
+    j=numskip+i
+    k=nlev_ec+1+numskip+i
+    akm_usort(nwz-i+1)=pres(nwz-i+1)
+    bkm(nwz-i+1)=0.0
+  end do
+
+  !******************************
+  ! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
+  !******************************
+      do i=1,nwz
+         if (akm_usort(1).gt.akm_usort(2)) then
+            akm(i)=akm_usort(i)
+         else
+            akm(i)=akm_usort(nwz-i+1)
+         endif
+      end do
+
+  !
+  ! CALCULATION OF AKZ, BKZ
+  ! AKZ,BKZ: model discretization parameters at the center of each model
+  !     layer
+  !
+  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
+  ! i.e. ground level
+  !*****************************************************************************
+
+  do i=1,nuvz
+     akz(i)=akm(i)
+     bkz(i)=bkm(i)
+  end do
+
+  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+  ! upon the transformation to z levels. In order to save computer memory, this is
+  ! not done anymore in the standard version. However, this option can still be
+  ! switched on by replacing the following lines with those below, that are
+  ! currently commented out. For this, similar changes are necessary in
+  ! verttransform.f and verttranform_nests.f
+  !*****************************************************************************
+
+  nz=nuvz
+  if (nz.gt.nzmax) stop 'nzmax too small'
+  do i=1,nuvz
+    aknew(i)=akz(i)
+    bknew(i)=bkz(i)
+  end do
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !nz=nuvz+nwz-1
+  !if (nz.gt.nzmax) stop 'nzmax too small'
+  !do 100 i=1,nwz
+  !  aknew(2*(i-1)+1)=akm(i)
+  !00     bknew(2*(i-1)+1)=bkm(i)
+  !do 110 i=2,nuvz
+  !  aknew(2*(i-1))=akz(i)
+  !10     bknew(2*(i-1))=bkz(i)
+  ! End doubled vertical resolution
+
+
+  ! Determine the uppermost level for which the convection scheme shall be applied
+  ! by assuming that there is no convection above 50 hPa (for standard SLP)
+  !*****************************************************************************
+
+  do i=1,nuvz-2
+    pint=akz(i)+bkz(i)*101325.
+    if (pint.lt.5000.) goto 96
+  end do
+96   nconvlev=i
+  if (nconvlev.gt.nconvlevmax-1) then
+    nconvlev=nconvlevmax-1
+    write(*,*) 'Attention, convection only calculated up to ', &
+         akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
+  endif
+
+  return
+
+999   write(*,*)
+  write(*,*) ' ###########################################'// &
+       '###### '
+  write(*,*) '       TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
+  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
+  write(*,*) ' ###########################################'// &
+       '###### '
+  write(*,*)
+  write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND   !!!'
+  write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE...     !!!'
+  write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
+  write(*,'(a)') '!!! THE "X" KEY...                   !!!'
+  write(*,*)
+  read(*,'(a)') opt
+  if(opt.eq.'X') then
+    stop
+  else
+    goto 5
+  endif
+
+end subroutine gridcheck
Index: src/interpol_rain.f90
===================================================================
--- src/interpol_rain.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/interpol_rain.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -20,8 +20,15 @@
 !**********************************************************************
 
-subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
-       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
-  !                          i   i   i    i    i     i   i
-  !i    i    i  i    i     i      i      i     o     o     o
+!subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
+!       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
+!  !                          i   i   i    i    i     i   i
+!  !i    i    i  i    i     i      i      i     o     o     o
+
+      subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
+     ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
+     intiy1,intiy2,icmv)
+!                               i   i   i    i    i     i   i 
+!     i    i    i  i    i     i      i      i     o     o     o
+
   !****************************************************************************
   !                                                                           *
@@ -40,5 +47,8 @@
   !                                                                           *
   !     30 August 1996                                                        *
-  !                                                                           *
+  !
+  !*  Petra Seibert, 2011/2012:
+  !*  Add interpolation of cloud bottom and cloud thickness
+  !*  for fix to SE's new wet scavenging scheme  *
   !****************************************************************************
   !                                                                           *
@@ -70,10 +80,15 @@
 
   integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
-  integer :: itime,itime1,itime2,level,indexh
+  !integer :: itime,itime1,itime2,level,indexh
+  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 
+  integer :: intiy1,intiy2,ipsum,icmv  
   real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
-  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
-  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
+  integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
+  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
+  !real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
+  real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4  
+  !real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
 
 
@@ -127,4 +142,43 @@
          + p3*yy3(ix ,jyp,level,indexh) &
          + p4*yy3(ixp,jyp,level,indexh)
+
+!CPS clouds:
+        ip1=1
+        ip2=1
+        ip3=1
+        ip4=1
+        if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
+        if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
+        if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
+        if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
+        ipsum= ip1+ip2+ip3+ip4
+        if (ipsum .eq. 0) then
+          yi1(m)=icmv
+        else
+          yi1(m)=(ip1*p1*iy1(ix ,jy ,indexh)    &
+             + ip2*p2*iy1(ixp,jy ,indexh)       &
+             + ip3*p3*iy1(ix ,jyp,indexh)       &
+             + ip4*p4*iy1(ixp,jyp,indexh))/ipsum
+        endif
+        
+        ip1=1
+        ip2=1
+        ip3=1
+        ip4=1
+        if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
+        if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
+        if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
+        if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
+        ipsum= ip1+ip2+ip3+ip4
+        if (ipsum .eq. 0) then
+          yi2(m)=icmv
+        else
+          yi2(m)=(ip1*p1*iy2(ix ,jy ,indexh)  &
+             + ip2*p2*iy2(ixp,jy ,indexh)     &
+             + ip3*p3*iy2(ix ,jyp,indexh)     &
+             + ip4*p4*iy2(ixp,jyp,indexh))/ipsum
+        endif
+!CPS end clouds
+
   end do
 
@@ -134,4 +188,11 @@
   !************************************
 
+      if (abs(itime) .lt. abs(itime1)) then
+        print*,'interpol_rain.f90'
+        print*,itime,itime1,itime2
+        stop 'ITIME PROBLEM'
+      endif
+
+
   dt1=real(itime-itime1)
   dt2=real(itime2-itime)
@@ -143,3 +204,20 @@
 
 
+!PS clouds:
+      intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt
+      if (yi1(1) .eq. float(icmv)) intiy1=yi1(2) 
+      if (yi1(2) .eq. float(icmv)) intiy1=yi1(1) 
+
+      intiy2=(yi2(1)*dt2 + yi2(2)*dt1)/dt
+      if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)  
+      if (yi2(2) .eq. float(icmv)) intiy2=yi2(1) 
+      
+      if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
+        intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
+      else
+        intiy1=icmv
+        intiy2=icmv
+      endif
+!PS end clouds
+
 end subroutine interpol_rain
Index: src/makefile
===================================================================
--- src/makefile	(revision None)
+++ src/makefile	(revision None)
@@ -0,0 +1,119 @@
+SHELL = /bin/bash
+TARGET = local
+WINDS=ecmwf
+#WINDS=gfs
+#WINDS=fnl
+
+FC       = gfortran
+FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+
+ifeq ($(TARGET),dmz)	
+  # options for ganglia
+  INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
+  LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
+  LIBPATH2 = /usr/lib/x86_64-linux-gnu/
+  MAIN = FLEXPART_dmz_
+endif
+ifeq ($(TARGET),local)
+  # local options 	
+  #libs_dir=/.../flexpart/libs/
+  libs_dir=/Users/ignacio/flexpart/libs/
+  INCPATH  = $(libs_dir)/grib_api-1.9.9_dir/include
+  LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
+  LIBPATH2 = $(libs_dir)/jasper_dir/lib
+  MAIN = FLEXPART_local_
+endif
+
+LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
+
+MODOBJS = \
+par_mod.o    com_mod.o \
+conv_mod.o              hanna_mod.o \
+interpol_mod.o          cmapf_mod.o \
+unc_mod.o               oh_mod.o \
+xmass_mod.o             flux_mod.o \
+point_mod.o             outg_mod.o
+
+OBJECTS = \
+writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
+               part0.o \
+caldate.o               partdep.o \
+coordtrafo.o            psih.o \
+raerod.o \
+drydepokernel.o         random.o \
+erf.o                   readavailable.o \
+ew.o                    readcommand.o \
+advance.o               readdepo.o \
+releaseparticles.o      psim.o     \
+FLEXPART.o              readlanduse.o \
+getfields.o             init_domainfill.o\
+interpol_wind.o         readoutgrid.o \
+interpol_all.o          readpaths.o \
+getrb.o                 readreceptors.o \
+getrc.o                 readreleases.o \
+getvdep.o               readspecies.o \
+interpol_misslev.o      \
+conccalc.o              \
+concoutput.o  concoutput_surf.o          scalev.o \
+pbl_profile.o           readOHfield.o\
+juldate.o               timemanager.o \
+interpol_vdep.o         interpol_rain.o \
+partoutput.o \
+hanna.o                 wetdepokernel.o \
+mean.o                  wetdepo.o \
+hanna_short.o           windalign.o \
+hanna1.o                initialize.o \
+                           calcpar_nests.o \
+verttransform_nests.o   interpol_all_nests.o \
+interpol_wind_nests.o   interpol_misslev_nests.o \
+interpol_vdep_nests.o   interpol_rain_nests.o \
+getvdep_nests.o   gridcheck_nests.o \
+readwind_nests.o  \
+readageclasses.o        readpartpositions.o \
+calcfluxes.o            fluxoutput.o \
+qvsat.o                 skplin.o \
+convect43c.o               redist.o \
+sort2.o                 distance.o \
+centerofmass.o          plumetraj.o \
+openouttraj.o           calcpv.o \
+calcpv_nests.o          distance2.o \
+clustering.o            interpol_wind_short.o \
+interpol_wind_short_nests.o shift_field_0.o \
+shift_field.o           outgrid_init.o \
+openreceptors.o         boundcond_domainfill.o\
+partoutput_short.o      readoutgrid_nest.o \
+outgrid_init_nest.o     writeheader_nest.o writeheader_nest_surf.o \
+concoutput_nest.o concoutput_surf_nest.o      wetdepokernel_nest.o \
+drydepokernel_nest.o    zenithangle.o \
+ohreaction.o            getvdep_nests.o \
+initial_cond_calc.o     initial_cond_output.o \
+dynamic_viscosity.o     get_settling.o
+
+
+ifeq ($(WINDS),ecmwf)
+  OBJECTS_WINDS = \
+  calcpar.o          readwind.o \
+  richardson.o       verttransform.o \
+  obukhov.o          gridcheck.o  \
+  convmix.o          calcmatrix.o 
+endif
+
+ifeq ($(WINDS),gfs)
+  OBJECTS_WINDS = \
+  calcpar_gfs.o          readwind_gfs.o \
+  richardson_gfs.o       verttransform_gfs.o \
+  obukhov_gfs.o          gridcheck_gfs.o  \
+  convmix_gfs.o          calcmatrix_gfs.o
+endif
+
+$(MAIN): $(MODOBJS) $(OBJECTS)  $(OBJECTS_WINDS)
+	$(FC) *.o -o $(MAIN)_$(WINDS) $(LDFLAGS)
+
+$(OBJECTS): $(MODOBJS)
+
+%.o: %.f90
+	$(FC) -c $(FFLAGS) $<
+
+clean:
+	rm *.o *.mod
+
Index: src/makefile.ecmwf_absoft
===================================================================
--- src/makefile.ecmwf_absoft	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,92 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART
-INCF	= incl* 
-#
-
-FC       = /opt/absoft10.1/bin/f90
-#FC       = /opt/absoft11.5/bin/f95
-INCPATH  = /nilu2/home/flexpart/lib64/absoft/include
-LIBPATH1 = /nilu2/home/flexpart/lib64/absoft/lib
-LIBPATH2 = /nilu2/home/flexpart/lib64/absoft/lib
-FFLAGS   =   -O2 -I$(INCPATH) -p$(INCPATH) -m64 -mcmodel=medium
-#FFLAGS   =  -g -Rb -Rc -Rs  -s -I$(INCPATH) -p$(INCPATH) -m64 -mcmodel=medium
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o           assignland.o\
-calcpar.o               part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-                        raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
-concoutput.o            scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
-convect43c.o               redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o \
-concoutput_nest.o       wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
-
Index: src/makefile.ecmwf_gfortran
===================================================================
--- src/makefile.ecmwf_gfortran	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,94 +1,0 @@
-SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran
-#
-
-FC       = gfortran
-INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-#LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
-#LIBPATH2 =   /flex_wrk/flexpart/lib64/gfortran/lib/
-#FFLAGS   =   -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o           assignland.o\
-calcpar.o               part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
-concoutput.o            scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
-convect43c.o               redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o \
-concoutput_nest.o       wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
Index: src/makefile.ecmwf_gfortran_jfb
===================================================================
--- src/makefile.ecmwf_gfortran_jfb	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,91 +1,0 @@
-SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran
-#
-
-FC       = gfortran
-INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 = /xnilu_wrk/flex_wrk/lib64/jasper/lib
-#FFLAGS   =   -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o           assignland.o\
-calcpar.o               part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
-concoutput.o            scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
-convect43c.o               redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o \
-concoutput_nest.o       wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
Index: src/makefile.gfs_absoft
===================================================================
--- src/makefile.gfs_absoft	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,93 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART_GFS
-INCF	= incl* 
-#
-
-
-FC       = /opt/absoft10.1/bin/f95
-INCPATH  = /nilu2/home/flexpart/lib64/absoft/include
-LIBPATH1 = /nilu2/home/flexpart/lib64/absoft/lib
-LIBPATH2 = /nilu2/home/flexpart/lib64/absoft/lib
-FFLAGS   =  -s -O2 -I$(INCPATH) -p$(INCPATH) -m64 -mcmodel=medium
-#FFLAGS   =  -g  -Rb -Rc -Rs  -s  -B108  -YEXT_NAMES=LCS -I$(INCPATH)
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o           assignland.o\
-calcpar.o               part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
-concoutput.o            scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
-convect43c.o               redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o \
-concoutput_nest.o       wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-
-%.o: %.f90
-        $(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
Index: src/makefile.gfs_gfortran
===================================================================
--- src/makefile.gfs_gfortran	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,90 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART_GFS_GFORTRAN
-INCF	= incl* 
-#
-
-FC       = gfortran
-INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o           assignland.o\
-calcpar_gfs.o           part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind_gfs.o \
-conccalc.o              richardson_gfs.o \
-concoutput.o            scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform_gfs.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov_gfs.o           gridcheck_gfs.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix_gfs.o               calcmatrix_gfs.o \
-convect43c.o            redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o \
-concoutput_nest.o       wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
Index: src/makefile.henrik
===================================================================
--- src/makefile.henrik	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,100 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART_ICEALOT
-#
-
-FC       = gfortran
-## OUTSIDE TEST PATH
-#INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api_test/include
-#LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api_test/lib
-#LIBPATH2 = /xnilu_wrk/flex_wrk/bin64/jasper/lib
-# OUTSIDE PATH
-INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 = /xnilu_wrk/flex_wrk/bin64/jasper/lib
-# INSIDE PATH
-#INCPATH  = /nilu2/home/flexpart/lib64/gfortran/include
-#LIBPATH1 = /nilu2/home/flexpart/lib64/gfortran/lib
-#LIBPATH2 = /nilu2/home/flexpart/lib64/gfortran/lib
-#FFLAGS   =   -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o           assignland.o\
-calcpar.o               part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
-concoutput.o            scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
-convect43c.o               redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o \
-concoutput_nest.o       wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
Index: src/makefile.test
===================================================================
--- src/makefile.test	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,96 +1,0 @@
-SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran
-#
-
-FC       = gfortran
-#INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-#LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-#LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-INCPATH =   /usr/lib/
-LIBPATH1 =   /usr/lib/
-LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
-#LIBPATH2 =   /flex_wrk/flexpart/lib64/gfortran/lib/
-#FFLAGS   =   -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o           assignland.o\
-calcpar.o               part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
-concoutput.o            scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
-convect43c.o               redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o \
-concoutput_nest.o       wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
Index: src/openouttraj.f90
===================================================================
--- src/openouttraj.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/openouttraj.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -54,7 +54,7 @@
 
   if (ldirect.eq.1) then
-  write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime,'FLEXPART V8.2'
+  write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime, trim(flexversion)
   else
-  write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime,'FLEXPART V8.2'
+  write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime, trim(flexversion)
   endif
   write(unitouttraj,*) method,lsubgrid,lconvection
Index: src/outg_mod.f90
===================================================================
--- src/outg_mod.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/outg_mod.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -43,4 +43,5 @@
   real,allocatable, dimension (:,:) :: wetgridsigma
   real,allocatable, dimension (:) :: sparse_dump_r
+  real,allocatable, dimension (:) :: sparse_dump_u
   integer,allocatable, dimension (:) :: sparse_dump_i
 
Index: src/outgrid_init.f90
===================================================================
--- src/outgrid_init.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/outgrid_init.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -247,4 +247,9 @@
        max(numygrid,numygridn)*numzgrid),stat=stat)
     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
+
+   allocate(sparse_dump_u(max(numxgrid,numxgridn)* &
+       max(numygrid,numygridn)*numzgrid),stat=stat)
+        if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
+
   allocate(sparse_dump_i(max(numxgrid,numxgridn)* &
        max(numygrid,numygridn)*numzgrid),stat=stat)
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/par_mod.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -28,5 +28,5 @@
 !        1997                                                                  *
 !                                                                              *
-!        Last update 10 August 2000                                            *
+!        Last update 15 August 2013 IP                                         *
 !                                                                              *
 !*******************************************************************************
@@ -120,10 +120,13 @@
   ! Maximum dimensions of the input mother grids
   !*********************************************
-
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92
-  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=61,nwzmax=61,nzmax=61
+  
+  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL
+  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
+  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
   !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
-  integer,parameter :: nxshift=359  ! for ECMWF
-  !integer,parameter :: nxshift=0     ! for GFS
+  !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
+
+  integer,parameter :: nxshift=359 ! for ECMWF
+  !integer,parameter :: nxshift=0     ! for GFS or FNL (XF)
 
   integer,parameter :: nconvlevmax = nuvzmax-1
@@ -150,7 +153,7 @@
   !*********************************************
 
-  integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
-  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
-
+  !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
+  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
+  integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
   ! maxnests                maximum number of nested grids
   ! nxmaxn,nymaxn           maximum dimension of nested wind fields in
@@ -194,6 +197,7 @@
   !**************************************************
 
-  integer,parameter :: maxpart=2000000
-  integer,parameter :: maxspec=4
+  !integer,parameter :: maxpart=4000000
+  integer,parameter :: maxpart=2000
+  integer,parameter :: maxspec=6
 
 
@@ -255,6 +259,18 @@
   integer,parameter :: unitlsm=1, unitsurfdata=1, unitland=1, unitwesely=1
   integer,parameter :: unitOH=1
-  integer,parameter :: unitdates=94, unitheader=90, unitshortpart=95
+  integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
   integer,parameter :: unitboundcond=89
 
+!******************************************************
+! integer code for missing values, used in wet scavenging (PS, 2012)
+!******************************************************
+
+    !  integer icmv
+    !  parameter(icmv=-9999)
+      integer,parameter ::  icmv=-9999
+
+! Parameters for testing
+!*******************************************
+!  integer :: verbosity=0
+
 end module par_mod
Index: src/readavailable.f90
===================================================================
--- src/readavailable.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readavailable.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -123,4 +123,6 @@
 
   do k=1,numbnests
+  print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
+  print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
     open(unitavailab,file=path(numpath+2*(k-1)+2) &
          (1:length(numpath+2*(k-1)+2)),status='old',err=998)
@@ -275,5 +277,5 @@
   return
 
-998   write(*,*) ' #### FLEXPART MODEL ERROR! FILE   #### '
+998   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE   #### '
   write(*,'(a)') '     '//path(numpath+2*(k-1)+2) &
        (1:length(numpath+2*(k-1)+2))
@@ -281,5 +283,5 @@
   stop
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
+999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '
   write(*,'(a)') '     '//path(4)(1:length(4))
   write(*,*) ' #### CANNOT BE OPENED           #### '
Index: src/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readcommand.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -80,85 +80,171 @@
   character(len=50) :: line
   logical :: old
-
+  logical :: nmlout=.true. !.false.
+  integer :: readerror
+
+  namelist /command/ &
+    ldirect, &
+    ibdate,ibtime, &
+    iedate,ietime, &
+    loutstep, &
+    loutaver, &
+    loutsample, &
+    itsplit, &
+    lsynctime, &
+    ctl, &
+    ifine, &
+    iout, &
+    ipout, &
+    lsubgrid, &
+    lconvection, &
+    lagespectra, &
+    ipin, &
+    ioutputforeachrelease, &
+    iflux, &
+    mdomainfill, &
+    ind_source, &
+    ind_receptor, &
+    mquasilag, &
+    nested_output, &
+    linit_cond, &
+    surf_only    
+
+  ! Presetting namelist command
+  ldirect=1
+  ibdate=20000101
+  ibtime=0
+  iedate=20000102
+  ietime=0
+  loutstep=10800
+  loutaver=10800
+  loutsample=900
+  itsplit=999999999
+  lsynctime=900
+  ctl=-5.0
+  ifine=4
+  iout=3
+  ipout=0
+  lsubgrid=1
+  lconvection=1
+  lagespectra=0
+  ipin=1
+  ioutputforeachrelease=1
+  iflux=1
+  mdomainfill=0
+  ind_source=1
+  ind_receptor=1
+  mquasilag=0
+  nested_output=0
+  linit_cond=0
+  surf_only=0 
 
   ! Open the command file and read user options
-  !********************************************
-
-
+  ! Namelist input first: try to read as namelist file
+  !**************************************************************************
   open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
-       err=999)
-
-  ! Check the format of the COMMAND file (either in free format,
-  ! or using formatted mask)
-  ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
-  !**************************************************************************
-
-  call skplin(9,unitcommand)
-  read (unitcommand,901) line
-901   format (a)
-  if (index(line,'LDIRECT') .eq. 0) then
-    old = .false.
-  else
-    old = .true.
-  endif
-  rewind(unitcommand)
-
-  ! Read parameters
-  !****************
-
-  call skplin(7,unitcommand)
-  if (old) call skplin(1,unitcommand)
-
-  read(unitcommand,*) ldirect
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ibdate,ibtime
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) iedate,ietime
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) loutstep
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) loutaver
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) loutsample
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) itsplit
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lsynctime
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ctl
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ifine
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) iout
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ipout
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lsubgrid
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lconvection
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lagespectra
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ipin
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ioutputforeachrelease
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) iflux
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) mdomainfill
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ind_source
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ind_receptor
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) mquasilag
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) nested_output
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) linit_cond
+         form='formatted',iostat=readerror)
+  ! If fail, check if file does not exist
+  if (readerror.ne.0) then
+    
+    print*,'***ERROR: file COMMAND not found in ' 
+    print*, path(1)(1:length(1))//'COMMAND'
+    print*, 'Check your pathnames file.'
+    stop
+
+  endif
+
+  read(unitcommand,command,iostat=readerror)
   close(unitcommand)
 
+  ! If error in namelist format, try to open with old input code
+  if (readerror.ne.0) then
+
+    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
+         err=999)
+
+    ! Check the format of the COMMAND file (either in free format,
+    ! or using formatted mask)
+    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
+    !**************************************************************************
+
+    call skplin(9,unitcommand)
+    read (unitcommand,901) line
+  901   format (a)
+    if (index(line,'LDIRECT') .eq. 0) then
+      old = .false.
+    else
+      old = .true.
+    endif
+    rewind(unitcommand)
+
+    ! Read parameters
+    !****************
+
+    call skplin(7,unitcommand)
+    if (old) call skplin(1,unitcommand)
+
+    read(unitcommand,*) ldirect
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ibdate,ibtime
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) iedate,ietime
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) loutstep
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) loutaver
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) loutsample
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) itsplit
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lsynctime
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ctl
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ifine
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) iout
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ipout
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lsubgrid
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lconvection
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lagespectra
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ipin
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ioutputforeachrelease
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) iflux
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) mdomainfill
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ind_source
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ind_receptor
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) mquasilag
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) nested_output
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) linit_cond
+    close(unitcommand)
+
+  endif ! input format
+
+  ! write command file in namelist format to output directory if requested
+  if (nmlout.eqv..true.) then
+    !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
+    open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
+    write(unitcommand,nml=command)
+    close(unitcommand)
+     ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
+     ! write(unitheader,NML=COMMAND)
+     !close(unitheader) 
+  endif
+
   ifine=max(ifine,1)
-
 
   ! Determine how Markov chain is formulated (for w or for w/sigw)
@@ -371,5 +457,5 @@
   endif
 
-  if(lsubgrid.ne.1) then
+  if(lsubgrid.ne.1.and.verbosity.eq.0) then
     write(*,*) '             ----------------               '
     write(*,*) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS'
@@ -505,3 +591,7 @@
   stop
 
+1000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
+       write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+        write(*,'(a)') path(2)(1:length(1))
+        stop
 end subroutine readcommand
Index: src/readoutgrid.f90
===================================================================
--- src/readoutgrid.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readoutgrid.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -91,5 +91,7 @@
   if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
        .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
+    write(*,*) 'outlon0,outlat0:'
     write(*,*) outlon0,outlat0
+    write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
     write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
     write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
Index: src/readpartpositions.f90
===================================================================
--- src/readpartpositions.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readpartpositions.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -77,11 +77,6 @@
 
   read(unitpartin) numpointin
-  if (numpointin.ne.numpoint) then
-     write(*,*) ' ####              WARNING                     #### '
-     write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
-     write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
-     write(*,*) ' #### THIS IS JUST A CHECK TO SEE IF YOU READ  #### '
-     write(*,*) ' #### THE CORRECT PARTICLE DUMP FILE.          #### '
-  endif
+  if (numpointin.ne.numpoint) goto 995
+999 continue 
   do i=1,numpointin
     read(unitpartin)
@@ -152,4 +147,10 @@
   stop
 
+!995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
+995   write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
+  write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
+  write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
+!  stop
+  goto 999 
 
 996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
Index: src/readpaths.f90
===================================================================
--- src/readpaths.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readpaths.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine readpaths
+subroutine readpaths !(pathfile)
 
   !*****************************************************************************
@@ -30,4 +30,7 @@
   !                                                                            *
   !     1 February 1994                                                        *
+  !     last modified                                                          *
+  !     HS, 7.9.2012                                                           *
+  !     option to give pathnames file as command line option                   *
   !                                                                            *
   !*****************************************************************************
@@ -47,15 +50,30 @@
   implicit none
 
-  integer :: i
+  integer   :: i
+  character(256) :: string_test 
+  character(1) :: character_test 
 
   ! Read the pathname information stored in unitpath
   !*************************************************
 
-
-  open(unitpath,file='pathnames',status='old',err=999)
+  open(unitpath,file=trim(pathfile),status='old',err=999)
 
   do i=1,numpath
     read(unitpath,'(a)',err=998) path(i)
     length(i)=index(path(i),' ')-1
+
+    
+    string_test = path(i)
+    character_test = string_test(length(i):length(i))
+    !print*, 'character_test,  string_test ', character_test,  string_test 
+      if ((character_test .NE. '/') .AND. (i .LT. 4))  then
+         print*, 'WARNING: path not ending in /' 
+         print*, path(i)
+         path(i) = string_test(1:length(i)) // '/'
+         length(i)=length(i)+1
+         print*, 'fix: padded with /' 
+         print*, path(i)
+         print*, 'length(i) increased 1' 
+      endif
   end do
 
@@ -70,4 +88,5 @@
     length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
   end do
+  print*,length(5),length(6)
 
 
Index: src/readreleases.f90
===================================================================
--- src/readreleases.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readreleases.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -57,4 +57,6 @@
   ! weta, wetb          parameters to determine the wet scavenging coefficient *
   ! zpoint1,zpoint2     height range, over which release takes place           *
+  ! num_min_discrete    if less, release cannot be randomized and happens at   *
+  !                     time mid-point of release interval                     *
   !                                                                            *
   !*****************************************************************************
@@ -68,6 +70,7 @@
 
   integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
+  integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
-  real(kind=dp) :: jul1,jul2,juldate
+  real(kind=dp) :: jul1,jul2,julm,juldate
   character(len=50) :: line
   logical :: old
@@ -376,4 +379,5 @@
   jul1=juldate(id1,it1)
   jul2=juldate(id2,it2)
+  julm=(jul1+jul2)/2.
   if (jul1.gt.jul2) then
     write(*,*) 'FLEXPART MODEL ERROR'
@@ -391,6 +395,11 @@
         stop
       endif
-      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      if (npart(numpoint).gt.num_min_discrete) then
+        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      else
+        ireleasestart(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+      endif
     else if (ldirect.eq.-1) then
       if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
@@ -401,6 +410,11 @@
         stop
       endif
-      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      if (npart(numpoint).gt.num_min_discrete) then
+        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      else
+        ireleasestart(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+      endif
     endif
   endif
Index: src/readwind.f90
===================================================================
--- src/readwind.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readwind.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -101,4 +101,9 @@
   character(len=24) :: gribErrorMsg = 'Error reading grib file'
   character(len=20) :: gribFunction = 'readwind'
+
+  !HSO conversion of ECMWF etadot to etadot*dp/deta
+  logical :: etacon=.false.
+  real,parameter :: p00=101325.
+  real :: dak,dbk
 
   hflswitch=.false.
@@ -365,4 +370,20 @@
   endif
 
+  ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
+  if(etacon.eqv..true.) then
+    do k=1,nwzmax
+      dak=akm(k+1)-akm(k)
+      dbk=bkm(k+1)-bkm(k)
+      do i=0,nxmin1
+        do j=0,nymin1 
+          wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
+          if (k.gt.1) then 
+            wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
+          endif
+        end do
+      end do
+    end do
+  endif
+
   ! For global fields, assign the leftmost data column also to the rightmost
   ! data column; if required, shift whole grid by nxshift grid points
Index: src/readwind_gfs.f90
===================================================================
--- src/readwind_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/readwind_gfs.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -703,5 +703,4 @@
   endif
 
-
   if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
   if(iumax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/timemanager.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -128,5 +128,13 @@
   !**********************************************************************
 
-  !write (*,*) 'starting simulation'
+
+  if (verbosity.gt.0) then
+    write (*,*) 'timemanager> starting simulation'
+    if (verbosity.gt.1) then
+      CALL SYSTEM_CLOCK(count_clock)
+      WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+    endif     
+  endif
+
   do itime=0,ideltas,lsynctime
 
@@ -142,6 +150,10 @@
   !********************************************************************
 
-    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) &
+    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
+        if (verbosity.gt.0) then
+           write (*,*) 'timemanager> call wetdepo'
+        endif     
          call wetdepo(itime,lsynctime,loutnext)
+    endif
 
     if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
@@ -156,23 +168,56 @@
   !*************************************
 
-      if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) &
-           call convmix(itime)
+   if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
+        if (verbosity.gt.0) then
+           write (*,*) 'timemanager> call convmix -- backward'
+        endif         
+      call convmix(itime)
+        if (verbosity.gt.1) then
+          !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+          CALL SYSTEM_CLOCK(count_clock)
+          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+        endif 
+   endif
 
   ! Get necessary wind fields if not available
   !*******************************************
-
+    if (verbosity.gt.0) then
+           write (*,*) 'timemanager> call getfields'
+    endif 
     call getfields(itime,nstop1)
+        if (verbosity.gt.1) then
+          CALL SYSTEM_CLOCK(count_clock)
+          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+        endif 
     if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
+
   ! Release particles
   !******************
 
+    if (verbosity.gt.0) then
+           write (*,*) 'timemanager>  Release particles'
+    endif 
+
     if (mdomainfill.ge.1) then
       if (itime.eq.0) then
+        if (verbosity.gt.0) then
+          write (*,*) 'timemanager>  call init_domainfill'
+        endif       
         call init_domainfill
       else
+        if (verbosity.gt.0) then
+          write (*,*) 'timemanager>  call boundcond_domainfill'
+        endif   
         call boundcond_domainfill(itime,loutend)
       endif
     else
+      if (verbosity.gt.0) then
+        print*,'call releaseparticles'  
+      endif
       call releaseparticles(itime)
+      if (verbosity.gt.1) then
+        CALL SYSTEM_CLOCK(count_clock)
+        WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+      endif 
     endif
 
@@ -182,7 +227,10 @@
   !**************************************************************
 
-      if ((ldirect.eq.1).and.(lconvection.eq.1)) &
-           call convmix(itime)
-
+   if ((ldirect.eq.1).and.(lconvection.eq.1)) then
+     if (verbosity.gt.0) then
+       write (*,*) 'timemanager> call convmix -- forward'
+     endif    
+     call convmix(itime)
+   endif
 
   ! If middle of averaging period of output fields is reached, accumulated
@@ -299,7 +347,24 @@
       if ((itime.eq.loutend).and.(outnum.gt.0.)) then
         if ((iout.le.3.).or.(iout.eq.5)) then
+          if (surf_only.ne.1) then 
           call concoutput(itime,outnum,gridtotalunc, &
                wetgridtotalunc,drygridtotalunc)
-          if (nested_output.eq.1) call concoutput_nest(itime,outnum)
+          else  
+  if (verbosity.eq.1) then
+     print*,'call concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+          call concoutput_surf(itime,outnum,gridtotalunc, &
+               wetgridtotalunc,drygridtotalunc)
+  if (verbosity.eq.1) then
+     print*,'called concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+          endif
+
+          if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
+          if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
           outnum=0.
         endif
Index: src/verttransform.f90
===================================================================
--- src/verttransform.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/verttransform.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -49,4 +49,6 @@
   ! Sabine Eckhardt, March 2007
   ! added the variable cloud for use with scavenging - descr. in com_mod
+  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
+  ! note that also other subroutines are affected by the fix
   !*****************************************************************************
   !                                                                            *
@@ -70,11 +72,13 @@
 
   integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
-  integer :: rain_cloud_above,kz_inv
+ integer :: rain_cloud_above,kz_inv !SE
+  integer icloudtop !PS
   real :: f_qvsat,pressure
-  real :: rh,lsp,convp
+ !real :: rh,lsp,convp
+  real :: rh,lsp,convp,prec,rhmin  
   real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
   real :: xlon,ylat,xlonr,dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2
+  real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -83,5 +87,7 @@
   real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
   real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+  logical lconvectprec
   real,parameter :: const=r_air/ga
+  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
 
   logical :: init = .true.
@@ -545,38 +551,108 @@
   !   create a cloud and rainout/washout field, clouds occur where rh>80%
   !   total cloudheight is stored at level 0
+
+
+
   do jy=0,nymin1
     do ix=0,nxmin1
-      rain_cloud_above=0
-      lsp=lsprec(ix,jy,1,n)
-      convp=convprec(ix,jy,1,n)
-      cloudsh(ix,jy,n)=0
-      do kz_inv=1,nz-1
-         kz=nz-kz_inv+1
-         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
-         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
-         clouds(ix,jy,kz,n)=0
-         if (rh.gt.0.8) then ! in cloud
-            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-               rain_cloud_above=1
-               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
-                    height(kz)-height(kz-1)
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
-               else
-                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
-               endif
-            else ! no precipitation
-                  clouds(ix,jy,kz,n)=1 ! cloud
+
+
+
+  !    rain_cloud_above=0
+  !    lsp=lsprec(ix,jy,1,n)
+  !    convp=convprec(ix,jy,1,n)
+  !    cloudsh(ix,jy,n)=0
+  !    do kz_inv=1,nz-1
+  !       kz=nz-kz_inv+1
+  !       pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+  !       rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+  !       clouds(ix,jy,kz,n)=0
+  !       if (rh.gt.0.8) then ! in cloud
+  !          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+  !             rain_cloud_above=1
+  !             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+  !                  height(kz)-height(kz-1)
+  !             if (lsp.ge.convp) then
+  !                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+  !             else
+  !                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+  !             endif
+  !          else ! no precipitation
+  !                clouds(ix,jy,kz,n)=1 ! cloud
+  !          endif
+  !       else ! no cloud
+  !          if (rain_cloud_above.eq.1) then ! scavenging
+  !             if (lsp.ge.convp) then
+  !                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+  !             else
+  !                clouds(ix,jy,kz,n)=4 ! convp dominated washout
+  !             endif
+  !          endif
+  !       endif
+  !    end do
+
+
+   ! PS 3012
+
+             lsp=lsprec(ix,jy,1,n)
+          convp=convprec(ix,jy,1,n)
+          prec=lsp+convp
+          if (lsp.gt.convp) then !  prectype='lsp'
+            lconvectprec = .false.
+          else ! prectype='cp '
+            lconvectprec = .true.
+          endif
+          rhmin = 0.90 ! standard condition for presence of clouds
+!PS       note that original by Sabine Eckhart was 80%
+!PS       however, for T<-20 C we consider saturation over ice
+!PS       so I think 90% should be enough          
+          icloudbot(ix,jy,n)=icmv
+          icloudtop=icmv ! this is just a local variable
+98        do kz=1,nz
+            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
+            if (rh .gt. rhmin) then
+              if (icloudbot(ix,jy,n) .eq. icmv) then
+                icloudbot(ix,jy,n)=nint(height(kz))
+              endif
+              icloudtop=nint(height(kz)) ! use int to save memory
             endif
-         else ! no cloud
-            if (rain_cloud_above.eq.1) then ! scavenging
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
-               else
-                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
-               endif
+          enddo
+
+!PS try to get a cloud thicker than 50 m 
+!PS if there is at least .01 mm/h  - changed to 0.002 and put into
+!PS parameter precpmin        
+          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
+              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
+              prec .gt. precmin) then
+            rhmin = rhmin - 0.05
+            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
+          endif
+!PS implement a rough fix for badly represented convection
+!PS is based on looking at a limited set of comparison data
+          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
+             prec .gt. precmin) then 
+
+            if (convp .lt. 0.1) then
+              icloudbot(ix,jy,n) = 500
+              icloudtop =         8000
+            else
+              icloudbot(ix,jy,n) = 0
+              icloudtop =      10000
             endif
-         endif
-      end do
+          endif
+          if (icloudtop .ne. icmv) then
+            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
+          else
+            icloudthck(ix,jy,n) = icmv
+          endif
+!PS  get rid of too thin clouds      
+          if (icloudthck(ix,jy,n) .lt. 50) then
+            icloudbot(ix,jy,n)=icmv
+            icloudthck(ix,jy,n)=icmv
+          endif
+
+
     end do
   end do
@@ -605,3 +681,5 @@
   !104   continue
   ! close(4)
+
+
 end subroutine verttransform
Index: src/wetdepo.f90
===================================================================
--- src/wetdepo.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/wetdepo.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -39,4 +39,7 @@
   ! Code may not be correct for decay of deposition!                           *
   !                                                                            *
+  ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
+  ! scheme, not dated
+  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
   !*****************************************************************************
   !                                                                            *
@@ -71,9 +74,11 @@
 
   integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
-  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
-  integer :: ks, kp
+  integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
+  integer :: kz !PS scheme
+  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
+  integer :: scheme_number ! NIK==1, PS ==2
   real :: S_i, act_temp, cl, cle ! in cloud scavenging
   real :: clouds_h ! cloud height for the specific grid point
-  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
+  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -146,20 +151,38 @@
     interp_time=nint(itime-0.5*ltsample)
 
-    if (ngrid.eq.0) then
-      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
-           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
-           memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    else
-      call interpol_rain_nests(lsprecn,convprecn,tccn, &
-           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
-           memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    endif
-
-    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+!    PS nest case still needs to be implemented!!    
+!    if (ngrid.eq.0) then
+!      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
+!           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+           call interpol_rain(lsprec,convprec,tcc,     &
+            icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
+            memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
+            memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv)  
+!    else
+!      call interpol_rain_nests(lsprecn,convprecn,tccn, &
+!           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+!    endif
+
+
+ !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+ !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
+        prec = lsp+convp
+        precsub = 0.01
+        if (prec .lt. precsub) then
+          goto 20
+        else
+          f = (prec-precsub)/prec
+          lsp = f*lsp
+          convp = f*convp
+        endif
+
 
   ! get the level were the actual particle is in
       do il=2,nz
         if (height(il).gt.ztra1(jpart)) then
-          hz=il-1
+          !hz=il-1
+          kz=il-1
           goto 26
         endif
@@ -174,15 +197,31 @@
   ! scavenging is done
 
-  if (ngrid.eq.0) then
-     clouds_v=clouds(ix,jy,hz,n)
-     clouds_h=cloudsh(ix,jy,n)
-  else
-     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
-     clouds_h=cloudsnh(ix,jy,n,ngrid)
-  endif
-  !write(*,*) 'there is
-  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
-  if (clouds_v.le.1) goto 20
-  !write (*,*) 'there is scavenging'
+!old scheme
+!  if (ngrid.eq.0) then
+!     clouds_v=clouds(ix,jy,hz,n)
+!     clouds_h=cloudsh(ix,jy,n)
+!  else
+!     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
+!     clouds_h=cloudsnh(ix,jy,n,ngrid)
+!  endif
+!  !write(*,*) 'there is
+!  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
+!  if (clouds_v.le.1) goto 20
+!  !write (*,*) 'there is scavenging'
+
+  ! PS: part of 2011/2012 fix 
+
+        if (ztra1(jpart) .le. float(ictop)) then
+          if (ztra1(jpart) .gt. float(icbot)) then
+            indcloud = 2 ! in-cloud
+          else
+            indcloud = 1 ! below-cloud
+          endif
+        elseif (ictop .eq. icmv) then
+          indcloud = 0 ! no cloud found, use old scheme
+        else
+          goto 20 ! above cloud
+        endif
+
 
   ! 1) Parameterization of the the area fraction of the grid cell where the
@@ -228,6 +267,12 @@
   !**********************************************************
 
-    do ks=1,nspec                                  ! loop over species
+    do ks=1,nspec            ! loop over species
       wetdeposit(ks)=0.
+
+    
+    !conflicting changes to the same routine: 1=NIK 2 =PS
+    scheme_number=2   
+    if (scheme_number.eq.1) then !NIK 
+
       if (weta(ks).gt.0.) then
         if (clouds_v.ge.4) then
@@ -236,4 +281,5 @@
           wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
   !        write(*,*) 'bel. wetscav: ',wetscav
+
         else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
   !        IN CLOUD SCAVENGING
@@ -245,12 +291,18 @@
              act_temp=tt(ix,jy,hz,n)
           endif
-          cl=2E-7*prec**0.36
+
+! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
+          cl=weta_in(ks)*prec**wetb_in(ks)
           if (dquer(ks).gt.0) then ! is particle
-            S_i=0.9/cl
+            S_i=wetc_in(ks)/cl
            else ! is gas
             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
             S_i=1/cle
            endif
-           wetscav=S_i*prec/3.6E6/clouds_h
+           wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
   !         write(*,*) 'in. wetscav:'
   !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
@@ -289,4 +341,72 @@
          wetdeposit(ks)=0.
       endif ! weta(k)
+
+     elseif (scheme_number.eq.2) then ! PS
+
+!PS          indcloud=0 ! Use this for FOR TESTING, 
+!PS          will skip the new in/below cloud method !!!             
+
+          if (weta(ks).gt.0.) then
+            if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
+!C               for aerosols and not highliy soluble substances weta=5E-6
+              wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
+!c             write(*,*) 'bel. wetscav: ',wetscav
+            elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
+              if (ngrid.gt.0) then
+                 act_temp=ttn(ix,jy,kz,n,ngrid)
+              else
+                 act_temp=tt(ix,jy,kz,n)
+              endif
+
+! from NIK              
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+
+
+              cl=2E-7*prec**0.36
+              if (dquer(ks).gt.0) then ! is particle
+                S_i=0.9/cl
+              else ! is gas
+                cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+                S_i=1/cle
+              endif
+              wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
+            else ! PS: no cloud diagnosed, old scheme, 
+!CPS          using with fixed a,b for simplicity, one may wish to change!!
+              wetscav = 1.e-4*prec**0.62
+            endif
+
+
+            wetdeposit(ks)=xmass1(jpart,ks)*    &
+            ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
+             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS) 
+            restmass = xmass1(jpart,ks)-wetdeposit(ks)
+            if (ioutputforeachrelease.eq.1) then
+              kp=npoint(jpart)
+            else
+              kp=1
+            endif
+            if (restmass .gt. smallnum) then
+              xmass1(jpart,ks)=restmass
+!cccccccccccccccc depostatistic
+!c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+!cccccccccccccccc depostatistic
+            else
+              xmass1(jpart,ks)=0.
+            endif
+!C Correct deposited mass to the last time step when radioactive decay of 
+!C gridded deposited mass was calculated
+            if (decay(ks).gt.0.) then
+              wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
+            endif
+          else  ! weta(k)<0
+             wetdeposit(ks)=0.
+          endif 
+
+     endif !on scheme
+
+
+
     end do
 
@@ -295,11 +415,21 @@
   !*****************************************************************************
 
-    if (ldirect.eq.1) then
-    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
-         real(ytra1(jpart)),nage,kp)
-    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
-         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
-         nage,kp)
-    endif
+!    if (ldirect.eq.1) then
+!    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
+!         real(ytra1(jpart)),nage,kp)
+!    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
+!         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
+!         nage,kp)
+!    endif
+
+   !PS
+        if (ldirect.eq.1) then
+          call wetdepokernel(nclass(jpart),wetdeposit,     &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
+          if (nested_output.eq.1)                          &
+           call wetdepokernel_nest(nclass(jpart),wetdeposit, &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)  
+        endif
+
 
 20  continue
Index: src/writeheader.f90
===================================================================
--- src/writeheader.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/writeheader.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -67,7 +67,7 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime,'FLEXPART V9.0'
+    write(unitheader) ibdate,ibtime, trim(flexversion)
   else
-    write(unitheader) iedate,ietime,'FLEXPART V9.0'
+    write(unitheader) iedate,ietime, trim(flexversion)
   endif
 
Index: src/writeheader_nest_surf.f90
===================================================================
--- src/writeheader_nest_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/writeheader_nest_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -0,0 +1,156 @@
+!**********************************************************************
+! 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 writeheader_nest_surf
+
+  !*****************************************************************************
+  !                                                                            *
+  !  This routine produces a file header containing basic information on the   *
+  !  settings of the FLEXPART run.                                             *
+  !  The header file is essential and must be read in by any postprocessing    *
+  !  program before reading in the output data.                                *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     7 August 2002                                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  !                                                                            *
+  ! xlon                   longitude                                           *
+  ! xl                     model x coordinate                                  *
+  ! ylat                   latitude                                            *
+  ! yl                     model y coordinate                                  *
+  !                                                                            *
+  !*****************************************************************************
+
+  use point_mod
+  use outg_mod
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  integer :: jjjjmmdd,ihmmss,i,ix,jy,j
+  real :: xp1,yp1,xp2,yp2
+
+
+  !************************
+  ! Open header output file
+  !************************
+
+  open(unitheader,file=path(2)(1:length(2))//'header_nest_grid_time', &
+       form='unformatted',err=998)
+
+
+  ! Write the header information
+  !*****************************
+
+  if (ldirect.eq.1) then
+    write(unitheader) ibdate,ibtime,trim(flexversion)
+  else
+    write(unitheader) iedate,ietime,trim(flexversion)
+  endif
+
+  ! Write info on output interval, averaging time, sampling time
+  !*************************************************************
+
+  write(unitheader) loutstep,loutaver,loutsample
+
+  ! Write information on output grid setup
+  !***************************************
+
+  write(unitheader) outlon0n,outlat0n,numxgridn,numygridn, &
+       dxoutn,dyoutn
+  write(unitheader) 1,(outheight(1),i=1,1)
+
+  call caldate(bdate,jjjjmmdd,ihmmss)
+  write(unitheader) jjjjmmdd,ihmmss
+
+  ! Write number of species, and name for each species (+extra name for depositions)
+  ! Indicate the dimension of the fields (i.e., 1 for deposition fields, numzgrid for
+  ! concentration fields
+  !*****************************************************************************
+
+  write(unitheader) 3*nspec,maxpointspec_act
+  do i=1,nspec
+    write(unitheader) 1,'WD_'//species(i)(1:7)
+    write(unitheader) 1,'DD_'//species(i)(1:7)
+    write(unitheader) 1,species(i)
+  end do
+
+  ! Write information on release points: total number, then for each point:
+  ! start, end, coordinates, # of particles, name, mass
+  !************************************************************************
+
+  write(unitheader) numpoint
+  do i=1,numpoint
+    write(unitheader) ireleasestart(i),ireleaseend(i),kindz(i)
+    xp1=xpoint1(i)*dx+xlon0
+    yp1=ypoint1(i)*dy+ylat0
+    xp2=xpoint2(i)*dx+xlon0
+    yp2=ypoint2(i)*dy+ylat0
+    write(unitheader) xp1,yp1,xp2,yp2,zpoint1(i),zpoint2(i)
+    write(unitheader) npart(i),1
+    if (numpoint.le.1000) then
+      write(unitheader) compoint(i)
+    else
+      write(unitheader) compoint(1001)
+   endif
+    do j=1,nspec
+      write(unitheader) xmass(i,j)
+      write(unitheader) xmass(i,j)
+      write(unitheader) xmass(i,j)
+    end do
+  end do
+
+  ! Write information on some model switches
+  !*****************************************
+
+  write(unitheader) method,lsubgrid,lconvection, &
+       ind_source,ind_receptor
+
+  ! Write age class information
+  !****************************
+
+  write(unitheader) nageclass,(lage(i),i=1,nageclass)
+
+
+  ! Write topography to output file
+  !********************************
+
+  do ix=0,numxgridn-1
+    write(unitheader) (orooutn(ix,jy),jy=0,numygridn-1)
+  end do
+  close(unitheader)
+
+  return
+
+
+998   write(*,*) ' #### FLEXPART MODEL ERROR!   THE FILE         #### '
+  write(*,*) ' #### '//path(2)(1:length(2))//'header'//' #### '
+  write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS    #### '
+  write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
+  write(*,*) ' #### THE PROGRAM AGAIN.                       #### '
+  stop
+
+end subroutine writeheader_nest_surf
Index: src/writeheader_surf.f90
===================================================================
--- src/writeheader_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/writeheader_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -0,0 +1,156 @@
+!**********************************************************************
+! 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 writeheader_surf
+
+  !*****************************************************************************
+  !                                                                            *
+  !  This routine produces a file header containing basic information on the   *
+  !  settings of the FLEXPART run.                                             *
+  !  The header file is essential and must be read in by any postprocessing    *
+  !  program before reading in the output data.                                *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     7 August 2002                                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  !                                                                            *
+  ! xlon                   longitude                                           *
+  ! xl                     model x coordinate                                  *
+  ! ylat                   latitude                                            *
+  ! yl                     model y coordinate                                  *
+  !                                                                            *
+  !*****************************************************************************
+
+  use point_mod
+  use outg_mod
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  integer :: jjjjmmdd,ihmmss,i,ix,jy,j
+  real :: xp1,yp1,xp2,yp2
+
+
+  !************************
+  ! Open header output file
+  !************************
+
+  open(unitheader,file=path(2)(1:length(2))//'header_grid_time', &
+       form='unformatted',err=998)
+
+
+  ! Write the header information
+  !*****************************
+
+  if (ldirect.eq.1) then
+    write(unitheader) ibdate,ibtime, trim(flexversion)
+  else
+    write(unitheader) iedate,ietime, trim(flexversion)
+  endif
+
+  ! Write info on output interval, averaging time, sampling time
+  !*************************************************************
+
+  write(unitheader) loutstep,loutaver,loutsample
+
+  ! Write information on output grid setup
+  !***************************************
+
+  write(unitheader) outlon0,outlat0,numxgrid,numygrid, &
+       dxout,dyout
+  write(unitheader) 1,(outheight(1),i=1,1)
+
+  call caldate(bdate,jjjjmmdd,ihmmss)
+  write(unitheader) jjjjmmdd,ihmmss
+
+  ! Write number of species, and name for each species (+extra name for depositions)
+  ! Indicate the dimension of the fields (i.e., 1 for deposition fields, numzgrid for
+  ! concentration fields
+  !*****************************************************************************
+
+  write(unitheader) 3*nspec,maxpointspec_act
+  do i=1,nspec
+    write(unitheader) 1,'WD_'//species(i)(1:7)
+    write(unitheader) 1,'DD_'//species(i)(1:7)
+    write(unitheader) 1,species(i)
+  end do
+
+  ! Write information on release points: total number, then for each point:
+  ! start, end, coordinates, # of particles, name, mass
+  !************************************************************************
+
+  write(unitheader) numpoint
+  do i=1,numpoint
+    write(unitheader) ireleasestart(i),ireleaseend(i),kindz(i)
+    xp1=xpoint1(i)*dx+xlon0
+    yp1=ypoint1(i)*dy+ylat0
+    xp2=xpoint2(i)*dx+xlon0
+    yp2=ypoint2(i)*dy+ylat0
+    write(unitheader) xp1,yp1,xp2,yp2,zpoint1(i),zpoint2(i)
+    write(unitheader) npart(i),1
+    if (numpoint.le.1000) then
+      write(unitheader) compoint(i)
+    else
+      write(unitheader) compoint(1001)
+    endif
+    do j=1,nspec
+      write(unitheader) xmass(i,j)
+      write(unitheader) xmass(i,j)
+      write(unitheader) xmass(i,j)
+    end do
+  end do
+
+  ! Write information on some model switches
+  !*****************************************
+
+  write(unitheader) method,lsubgrid,lconvection, &
+       ind_source,ind_receptor
+
+  ! Write age class information
+  !****************************
+
+  write(unitheader) nageclass,(lage(i),i=1,nageclass)
+
+
+  ! Write topography to output file
+  !********************************
+
+  do ix=0,numxgrid-1
+    write(unitheader) (oroout(ix,jy),jy=0,numygrid-1)
+  end do
+  close(unitheader)
+
+  return
+
+
+998   write(*,*) ' #### FLEXPART MODEL ERROR!   THE FILE         #### '
+  write(*,*) ' #### '//path(2)(1:length(2))//'header'//' #### '
+  write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS    #### '
+  write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
+  write(*,*) ' #### THE PROGRAM AGAIN.                       #### '
+  stop
+
+end subroutine writeheader_surf
Index: src/writeheader_txt.f90
===================================================================
--- src/writeheader_txt.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/writeheader_txt.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
@@ -0,0 +1,190 @@
+!**********************************************************************
+! 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 writeheader_txt
+
+  !*****************************************************************************
+  !                                                                            *
+  !  This routine produces a file header containing basic information on the   *
+  !  settings of the FLEXPART run.                                             *
+  !  The header file is essential and must be read in by any postprocessing    *
+  !  program before reading in the output data.                                *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     7 August 2002                                                          *
+  !     modified IP 2013 for text output                                       *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  !                                                                            *
+  ! xlon                   longitude                                           *
+  ! xl                     model x coordinate                                  *
+  ! ylat                   latitude                                            *
+  ! yl                     model y coordinate                                  *
+  !                                                                            *
+  !*****************************************************************************
+
+  use point_mod
+  use outg_mod
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  integer :: jjjjmmdd,ihmmss,i,ix,jy,j
+  real :: xp1,yp1,xp2,yp2
+
+
+  !************************
+  ! Open header output file
+  !************************
+
+  open(unitheader,file=path(2)(1:length(2))//'header_txt', &
+       form='formatted',err=998)
+  open(unitheader_txt,file=path(2)(1:length(2))//'header_txt_releases', &
+       form='formatted',err=998)
+
+
+  ! Write the header information
+  !*****************************
+  
+  write(unitheader,*) '# ibdate,ibtime, iedate, ietime, flexversion'
+  write(unitheader,*) ibdate, ibtime, iedate, ietime, trim(flexversion) !  'FLEXPART V9.0'
+  !if (ldirect.eq.1) then
+  !  write(unitheader,*) ibdate,ibtime,trim(flexversion) !  'FLEXPART V9.0'
+  !else
+  !  write(unitheader,*) iedate,ietime,trim(flexversion) ! 'FLEXPART V9.0'
+  !endif
+
+  ! Write info on output interval, averaging time, sampling time
+  !*************************************************************
+  
+  write(unitheader,*) '# interval, averaging time, sampling time'
+  write(unitheader,*) loutstep,loutaver,loutsample
+
+  ! Write information on output grid setup
+  !***************************************
+  
+  write(unitheader,*) '# information on grid setup    '
+  write(unitheader,*) '#outlon0,outlat0,numxgrid,numygrid,dxout,dyout'
+  write(unitheader,*) outlon0,outlat0,numxgrid,numygrid, &
+       dxout,dyout  
+  write(unitheader,*) '# numzgrid, outheight(.) '
+  write(unitheader,*) numzgrid,(outheight(i),i=1,numzgrid)
+
+  write(unitheader,*) '# jjjjmmdd,ihmmss'
+  call caldate(bdate,jjjjmmdd,ihmmss)
+  write(unitheader,*) jjjjmmdd,ihmmss
+
+  ! Write number of species, and name for each species (+extra name for depositions)
+  ! Indicate the vertical dimension of the fields (i.e., 1 for deposition fields, numzgrid for
+  ! concentration fields
+  !*****************************************************************************
+
+  write(unitheader,*) '# information on species'
+  write(unitheader,*) '# 3*nspec,maxpointspec_act'
+  write(unitheader,*) 3*nspec,maxpointspec_act
+  write(unitheader,*) '# for nspec:'
+  write(unitheader,*) '# 1, WD_ '
+  write(unitheader,*) '# 1, DD_ '
+  write(unitheader,*) '# numzgrid,species'
+  do i=1,nspec
+    write(unitheader,*) 1,'WD_'//species(i)(1:7)
+    write(unitheader,*) 1,'DD_'//species(i)(1:7)
+    write(unitheader,*) numzgrid,species(i)
+  end do
+
+  ! Write information on release points: total number, then for each point:
+  ! start, end, coordinates, # of particles, name, mass
+  !************************************************************************
+
+
+  write(unitheader_txt,*) '# information on release points'
+  write(unitheader_txt,*) '# numpoint'
+  write(unitheader_txt,*) numpoint
+  write(unitheader_txt,*) '# for numpoint:'
+  do i=1,numpoint
+    write(unitheader_txt,*) ireleasestart(i),ireleaseend(i),kindz(i)
+    xp1=xpoint1(i)*dx+xlon0
+    yp1=ypoint1(i)*dy+ylat0
+    xp2=xpoint2(i)*dx+xlon0
+    yp2=ypoint2(i)*dy+ylat0
+    write(unitheader_txt,*) xp1,yp1,xp2,yp2,zpoint1(i),zpoint2(i)
+    write(unitheader_txt,*) npart(i),1
+    if (numpoint.le.1000) then
+      write(unitheader_txt,*) compoint(i)
+    else
+      write(unitheader_txt,*) compoint(1001)
+    endif
+    do j=1,nspec
+      write(unitheader_txt,*) xmass(i,j)
+      write(unitheader_txt,*) xmass(i,j)
+      write(unitheader_txt,*) xmass(i,j)
+    end do
+  end do
+
+  ! Write information on model switches
+  !*****************************************
+
+  write(unitheader,*) '# information on model switches'
+  write(unitheader,*) '# method,lsubgrid,lconvection, ind_source,ind_receptor'
+  write(unitheader,*) method,lsubgrid,lconvection, &
+       ind_source,ind_receptor
+
+  ! Write age class information
+  !****************************
+  
+  write(unitheader,*) '# information on age class     '
+  write(unitheader,*) nageclass,(lage(i),i=1,nageclass)
+
+
+  !Do not write topography to text output file. Keep it on the binary one
+  !********************************
+
+  !do ix=0,numxgrid-1
+  !  write(unitheader,*) (oroout(ix,jy),jy=0,numygrid-1)
+  !end do
+
+
+  
+
+
+  close(unitheader)
+  close(unitheader_txt)
+
+
+!  open(unitheader,file=path(2)(1:length(2))//'header_nml', &
+!        form='formatted',err=998)
+!  write(unitheader,NML=COMMAND) 
+!  close(unitheader)
+
+  return
+
+
+998   write(*,*) ' #### FLEXPART MODEL ERROR!   THE FILE         #### '
+  write(*,*) ' #### '//path(2)(1:length(2))//'header_txt'//' #### '
+  write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS    #### '
+  write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
+  write(*,*) ' #### THE PROGRAM AGAIN.                       #### '
+  stop
+
+end subroutine writeheader_txt
