Index: options/COMMAND
===================================================================
--- options/COMMAND	(revision 6d73c4be98f6a308eed735650a968ac230fd19ef)
+++ options/COMMAND	(revision 6d73c4be98f6a308eed735650a968ac230fd19ef)
@@ -16,20 +16,20 @@
  ITSPLIT=        99999999, ! Interval of particle splitting (s) 
  LSYNCTIME=           900, ! All processes are synchronized to this time interval (s)
- CTL=          -5.0000000, ! CTL>1, ABL time step = (Lagrangian timescale (TL))/CTL, uses LSYNCTIME if CTL<0
+ CTL=          -5.0000000, ! Times step smaller than Lagr. time scale, LSYNCTIME if CTL<0, but can be shorter  
  IFINE=                 4, ! Reduction for time step in vertical transport, used only if CTL>1 
  IOUT=                  1, ! Output type: [1]mass 2]pptv 3]1&2 4]plume 5]1&4, +8 for NetCDF output     
  IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end   
  LSUBGRID=              0, ! Increase of ABL heights due to sub-grid scale orographic variations;[0]off 1]on 
- LCONVECTION=          	1, ! Switch for convection parameterization;0]off [1]on    
+ LCONVECTION=          	0, ! Switch for convection parameterization;[0]off 1]on    
  LAGESPECTRA=          	0, ! Switch for calculation of age spectra (needs AGECLASSES);[0]off 1]on  
- IPIN=                  0, ! Warm start from particle dump (needs previous partposit_end file); [0]no 1]yes  
+ IPIN=                  0, ! Warm start from particle dump (needs previous partposit_end file); 0]no [1]yes  
  IOUTPUTFOREACHRELEASE= 1, ! Separate output fields for each location in the RELEASE file; [0]no 1]yes 
  IFLUX=                 0, ! Output of mass fluxes through output grid box boundaries
  MDOMAINFILL=          	0, ! Switch for domain-filling, if limited-area particles generated at boundary
- IND_SOURCE=            1, ! Unit to be used at the source   ;  [1]mass 2]mass mixing ratio 
- IND_RECEPTOR=          1, ! Unit to be used at the receptor; [1]mass 2]mass mixing ratio 3]wet depo. 4]dry depo.
+ IND_SOURCE=            1, ! Unit to be used at the ource   ;  [1]mass 2]mass mixing ratio 
+ IND_RECEPTOR=          1, ! Unit to be used at the receptor;  [1]mass 2]mass mixing ratio 
  MQUASILAG=             0, ! Quasi-Lagrangian mode to track individual numbered particles 
  NESTED_OUTPUT=         0, ! Output also for a nested domain 
- LINIT_COND=            0, ! Output sensitivity to initial conditions (bkw mode only) [0]off 1]conc 2]mmr 
+ LINIT_COND=            0, ! Output sensitivity to initial conditions (backwards mode only)  
  SURF_ONLY=             0, ! Output only for the lowest model layer, used w/ LINIT_COND=1 or 2
  CBLFLAG=               0, ! Skewed, not Gaussian turbulence in the convective ABL, need large CTL and IFINE
Index: options/RELEASES
===================================================================
--- options/RELEASES	(revision 6d73c4be98f6a308eed735650a968ac230fd19ef)
+++ options/RELEASES	(revision 6d73c4be98f6a308eed735650a968ac230fd19ef)
@@ -1,30 +1,30 @@
-***************************************************************************************************************
-*                                                                                                             *
-*                                                                                                             *
-*                                                                                                             *
-*   Input file for the Lagrangian particle dispersion model FLEXPART                                          *
-*                        Please select your options                                                           *
-*                                                                                                             *
-*                                                                                                             *
-*                                                                                                             *
-***************************************************************************************************************
+*************************************************************************
+*                                                                       *
+*                                                                       *
+*                                                                       *
+*   Input file for the Lagrangian particle dispersion model FLEXPART    *
+*                        Please select your options                     *
+*                                                                       *
+*                                                                       *
+*                                                                       *
+*************************************************************************
 &RELEASES_CTRL
- NSPEC      =           1, ! Total number of species
- SPECNUM_REL=          24, ! Species numbers in directory SPECIES
+ NSPEC      =           1,
+ SPECNUM_REL=          24,
  /
-&RELEASE                   ! For each release 
- IDATE1  =       20120101, ! Release start date, YYYYMMDD: YYYY=year, MM=month, DD=day
- ITIME1  =         090000, ! Release start time in UTC HHMISS: HH hours, MI=minutes, SS=seconds
- IDATE2  =       20120101, ! Release end date, same as IDATE1
- ITIME2  =         090000, ! Release end time, same as ITIME1
- LON1    =          0.000, ! Left longitude of release box -180 < LON1 <180
- LON2    =          0.000, ! Right longitude of release box, same as LON1
- LAT1    =         20.000, ! Lower latitude of release box, -90 < LAT1 < 90
- LAT2    =         20.000, ! Upper latitude of release box same format as LAT1 
- Z1      =         50.000, ! Lower height of release box meters/hPa above reference level
- Z2      =         50.000, ! Upper height of release box meters/hPa above reference level
- ZKIND   =              1, ! Reference level 1=above ground, 2=above sea level, 3 for pressure in hPa
- MASS    =       1.0000E8, ! Total mass emitted, only relevant for fwd simulations
- PARTS   =          10000, ! Total number of particles to be released
- COMMENT =    "RELEASE 1", ! Comment, written in the outputfile
+&RELEASE
+ IDATE1  =       20120101,
+ ITIME1  =         090000,
+ IDATE2  =       20120101,
+ ITIME2  =         090000,
+ LON1    =          0.000,
+ LON2    =          0.000,
+ LAT1    =         20.000,
+ LAT2    =         20.000,
+ Z1      =         50.000,
+ Z2      =         50.000,
+ ZKIND   =              1,
+ MASS    =       1.0000E8,
+ PARTS   =          10000, 
+ COMMENT =    "RELEASE 1",
  /
Index: options/SPECIES/SPECIES_008
===================================================================
--- options/SPECIES/SPECIES_008	(revision 79996be9fcfb6fb924aa6f67d27d422820decb2f)
+++ options/SPECIES/SPECIES_008	(revision 8594fb83d316dd4c22c8ec7555c8b029f7f705c6)
@@ -19,7 +19,5 @@
  POHDCONST=-9.9,              ! OH Reaction rate - D [K]
  POHNCONST=2.0,               !  OH Reaction rate -  N (no unit)
- PAREA_HOUR=0.578, 0.491, 0.428, 0.329, 0.384, 0.485, 0.763, 1.103, 1.084, 1.047, 1.096, 1.196, 1.298, 1.357, 1.447, 1.565, 1.636, 1.662, 1.401, 1.168, 1.031, 0.926, 0.816, 0.709,
- PPOINT_HOUR=0.845, 0.806, 0.786, 0.779, 0.793, 0.832, 0.895, 0.977, 1.031, 1.071, 1.105, 1.118, 1.131, 1.136, 1.143, 1.141, 1.133, 1.118, 1.097, 1.091, 1.079, 1.036, 0.966, 0.892,
- PAREA_DOW=1.060, 1.060, 1.060, 1.060, 1.060, 0.900, 0.000,
- PPOINT_DOW=1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000,
+ PAREA_HOUR=0.578, 0.491, 0.428, 0.329, 0.384, 0.485, 0.763, 1.103, 1.084, 1.047, 1.096, 1.196, 1.298, 1.357, 1.447, 1.565, 1.636, 1.662, 1.401, 1.168, 1.031, 0.926, 0.816, 0.709
+ PAREA_DOW=1.060, 1.060, 1.060, 1.060, 1.060, 0.900, 0.000
  /
Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision 50958b85292a7b2ce570403ca6f859592eaa1b5e)
+++ src/FLEXPART.f90	(revision 50958b85292a7b2ce570403ca6f859592eaa1b5e)
@@ -455,19 +455,14 @@
   call timemanager(metdata_format)
 
-  if (verbosity.gt.0) then
 ! NIK 16.02.2005 
-    do i=1,nspec
-      if (tot_inc_count(i).gt.0) then
-         write(*,*) '**********************************************'
-         write(*,*) 'Scavenging statistics for species ', species(i), ':'
-         write(*,*) 'Total number of occurences of below-cloud scavenging', &
-           & tot_blc_count(i)
-         write(*,*) 'Total number of occurences of in-cloud    scavenging', &
-           & tot_inc_count(i)
-         write(*,*) '**********************************************'
-      endif
-    end do
-    write (*,*) 'timemanager> call wetdepo'
-  endif
+  do i=1,nspec
+    write(*,*) '**********************************************'
+    write(*,*) 'Scavenging statistics for species ', species(i), ':'
+    write(*,*) 'Total number of occurences of below-cloud scavenging', &
+         & tot_blc_count(i)
+    write(*,*) 'Total number of occurences of in-cloud    scavenging', &
+         & tot_inc_count(i)
+    write(*,*) '**********************************************'
+  end do
   
   write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision e9e0f0627da3a042c02a5352bf3446bdfbcfa136)
+++ src/com_mod.f90	(revision c2bd55ea39cbd0be88e7e8329b1c9c4e01754957)
@@ -174,6 +174,8 @@
   real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass)
   real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass)
-  real :: rm(maxspec),dryvel(maxspec)
+  real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
   real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
+  ! se  it is possible to associate a species with a second one to make transfer from gas to aerosol
+  integer :: spec_ass(maxspec)
 
   real :: area_hour(maxspec,24),point_hour(maxspec,24)
Index: src/get_wetscav.f90
===================================================================
--- src/get_wetscav.f90	(revision e9e0f0627da3a042c02a5352bf3446bdfbcfa136)
+++ src/get_wetscav.f90	(revision 79e0349a50fd1ee6dfdf4b044ade2b6f1e702855)
@@ -151,7 +151,10 @@
       if (height(il).gt.ztra1(jpart)) then
         hz=il-1
+!        goto 26
         exit
       endif
     end do
+!26  continue
+
 
     if (ngrid.eq.0) then
@@ -200,16 +203,23 @@
 
 
-!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp - 2 and 3 not used removed by SE
+!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp 
 ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
 ! for now they are treated the same
     grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
+    grfraction(2)=max(0.05,cc*(lfr(i)))
+    grfraction(3)=max(0.05,cc*(cfr(j)))
+
 
 ! 2) Computation of precipitation rate in sub-grid cell
 !******************************************************
     prec(1)=(lsp+convp)/grfraction(1)
+    prec(2)=(lsp)/grfraction(2)
+    prec(3)=(convp)/grfraction(3)
+
 
 ! 3) Computation of scavenging coefficients for all species
 !    Computation of wet deposition
 !**********************************************************
+
 
       if (ngrid.gt.0) then
@@ -227,4 +237,5 @@
 !******************************************************************
         if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
+          !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then 
           blc_count(ks)=blc_count(ks)+1
           wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
@@ -260,4 +271,6 @@
           endif
           
+!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
+
         endif ! gas or particle
 !      endif ! positive below-cloud scavenging parameters given in Species file
@@ -272,4 +285,5 @@
         if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then 
           inc_count(ks)=inc_count(ks)+1
+!          write(*,*) 'Incloud: ',inc_count
 ! if negative coefficients (turned off) set to zero for use in equation
           if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
@@ -286,5 +300,5 @@
 ! sec test
 !           cl=1E6*1E-7*prec(1)**0.3 !Sec GFS new
-            cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new, is also suitable for GFS
+            cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new
 !           cl=2E-7*prec(1)**0.36 !Andreas
 !           cl=1.6E-6*prec(1)**0.36 !Henrik
@@ -308,14 +322,25 @@
           if (dquer(ks).gt.0.) then
             S_i= frac_act/cl
+!           write(*,*) 'Si: ',S_i
+
 ! GAS
 !****
           else
+
             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+!REPLACE to switch old/ new scheme 
+          ! S_i=frac_act/cle
             S_i=1/cle
           endif ! gas or particle
 
 ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
-!SEC wetscav fix, the cloud height is no longer needed, it gives wrong results
+!OLD 
+          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
+          else
+!SEC wetscav fix
+             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
+!            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
+          endif
         endif ! positive in-cloud scavenging parameters given in Species file
       endif !incloud
Index: src/netcdf_output_mod.f90
===================================================================
--- src/netcdf_output_mod.f90	(revision 4ad96c54b793af523a98633b88f07aa60c69ea47)
+++ src/netcdf_output_mod.f90	(revision ae43937f34fb948c5aa7c5307dc0a9df9053f861)
@@ -67,5 +67,6 @@
                        ccn_aero,in_aero, & ! wetc_in,wetd_in, &
                        reldiff,henry,f0,density,dquer,dsigma,dryvel,&
-                       weightmolar,ohcconst,ohdconst,vsetaver,&
+!                       weightmolar,ohreact,spec_ass,kao,vsetaver,&
+                       weightmolar,ohcconst,ohdconst,spec_ass,kao,vsetaver,&
                        ! for concoutput_netcdf and concoutput_nest_netcdf
                        nxmin1,nymin1,nz,oro,oron,rho,rhon,&
@@ -113,6 +114,4 @@
   logical, parameter :: write_area = .false.
 
-  ! coordinate transformation from internal to world coord
-  real :: xp1,yp1,xp2,yp2
 contains
 
@@ -513,5 +512,7 @@
         call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
         call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
         call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
 
         if (lnest) then
@@ -534,5 +535,7 @@
         call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
         call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i))) 
         call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
 
         if (lnest) then
@@ -670,12 +673,8 @@
        call nf90_err(nf90_put_var(ncid, relendID, ireleaseend(i), (/i/)))
        call nf90_err(nf90_put_var(ncid, relkindzID, kindz(i), (/i/)))
-       xp1=xpoint1(i)*dx+xlon0
-       yp1=ypoint1(i)*dy+ylat0
-       xp2=xpoint2(i)*dx+xlon0
-       yp2=ypoint2(i)*dy+ylat0
-       call nf90_err(nf90_put_var(ncid, rellng1ID, xp1, (/i/)))
-       call nf90_err(nf90_put_var(ncid, rellng2ID, xp2, (/i/)))
-       call nf90_err(nf90_put_var(ncid, rellat1ID, yp1, (/i/)))
-       call nf90_err(nf90_put_var(ncid, rellat2ID, yp2, (/i/)))
+       call nf90_err(nf90_put_var(ncid, rellng1ID, xpoint1(i), (/i/)))
+       call nf90_err(nf90_put_var(ncid, rellng2ID, xpoint2(i), (/i/)))
+       call nf90_err(nf90_put_var(ncid, rellat1ID, ypoint1(i), (/i/)))
+       call nf90_err(nf90_put_var(ncid, rellat2ID, ypoint2(i), (/i/)))
        call nf90_err(nf90_put_var(ncid, relzz1ID, zpoint1(i), (/i/)))
        call nf90_err(nf90_put_var(ncid, relzz2ID, zpoint2(i), (/i/)))
Index: src/readspecies.f90
===================================================================
--- src/readspecies.f90	(revision d6a245b53fc61a135bafd85f2d253aa94c714293)
+++ src/readspecies.f90	(revision aa8c34a4e29198f0666be476133e8d1bc1d712f8)
@@ -67,8 +67,7 @@
   character(len=16) :: pspecies
   real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
-  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst
+  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pkao
   real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
-  real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
-  integer :: readerror
+  integer :: readerror, pspec_ass
 
 ! declare namelist
@@ -77,6 +76,5 @@
        pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
        preldiff, phenry, pf0, pdensity, pdquer, &
-       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, &
-       parea_dow, parea_hour, ppoint_dow, ppoint_hour
+       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
 
   pspecies="" ! read failure indicator value
@@ -98,16 +96,8 @@
   pohdconst=-9.9E-09
   pohnconst=2.0
+  pspec_ass=-9
+  pkao=-99.99
   pweightmolar=-999.9
 
-  do j=1,24           ! initialize everything to no variation
-    area_hour(pos_spec,j)=1.
-    point_hour(pos_spec,j)=1.
-  end do
-  do j=1,7
-    area_dow(pos_spec,j)=1.
-    point_dow(pos_spec,j)=1.
-  end do
-
-  if (readerror.ne.0) then ! text format input
 ! Open the SPECIES file and read species names and properties
 !************************************************************
@@ -174,16 +164,8 @@
     read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
 !  write(*,*) ohnconst(pos_spec)
-
-! Read in daily and day-of-week variation of emissions, if available
-!*******************************************************************
-
-    read(unitspecies,*,end=22)
-    do j=1,24     ! 24 hours, starting with 0-1 local time
-      read(unitspecies,*) ihour,area_hour(pos_spec,j),point_hour(pos_spec,j)
-    end do
-    read(unitspecies,*)
-    do j=1,7      ! 7 days of the week, starting with Monday
-      read(unitspecies,*) idow,area_dow(pos_spec,j),point_dow(pos_spec,j)
-    end do
+    read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
+!  write(*,*) spec_ass(pos_spec)
+    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
+!       write(*,*) kao(pos_spec)
 
     pspecies=species(pos_spec)
@@ -206,16 +188,8 @@
     pohdconst=ohdconst(pos_spec)
     pohnconst=ohnconst(pos_spec)
-
-
-    do j=1,24     ! 24 hours, starting with 0-1 local time
-      parea_hour(j)=area_hour(pos_spec,j)
-      ppoint_hour(j)=point_hour(pos_spec,j)
-    end do
-    do j=1,7      ! 7 days of the week, starting with Monday
-      parea_dow(j)=area_dow(pos_spec,j)
-      ppoint_dow(j)=point_dow(pos_spec,j)
-    end do
-
-  else ! namelist available
+    pspec_ass=spec_ass(pos_spec)
+    pkao=kao(pos_spec)
+
+  else
 
     species(pos_spec)=pspecies
@@ -238,13 +212,6 @@
     ohdconst(pos_spec)=pohdconst
     ohnconst(pos_spec)=pohnconst
-
-    do j=1,24     ! 24 hours, starting with 0-1 local time
-      area_hour(pos_spec,j)=parea_hour(j)
-      point_hour(pos_spec,j)=ppoint_hour(j)
-    end do
-    do j=1,7      ! 7 days of the week, starting with Monday
-      area_dow(pos_spec,j)=parea_dow(j)
-      point_dow(pos_spec,j)=ppoint_dow(j)
-    end do
+    spec_ass(pos_spec)=pspec_ass
+    kao(pos_spec)=pkao
 
   endif
@@ -336,4 +303,19 @@
   end if
 
+  if (spec_ass(pos_spec).gt.0) then
+    spec_found=.FALSE.
+    do j=1,pos_spec-1
+      if (spec_ass(pos_spec).eq.specnum(j)) then
+        spec_ass(pos_spec)=j
+        spec_found=.TRUE.
+        ASSSPEC=.TRUE.
+      endif
+    end do
+    if (spec_found.eqv..false.) then
+      goto 997
+    endif
+  endif
+
+  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
   if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
 
@@ -348,4 +330,28 @@
 
 
+! Read in daily and day-of-week variation of emissions, if available
+!*******************************************************************
+! HSO: This is not yet implemented as namelist parameters
+
+  do j=1,24           ! initialize everything to no variation
+    area_hour(i,j)=1.
+    point_hour(i,j)=1.
+  end do
+  do j=1,7
+    area_dow(i,j)=1.
+    point_dow(i,j)=1.
+  end do
+
+  if (readerror.ne.0) then ! text format input
+
+    read(unitspecies,*,end=22)
+    do j=1,24     ! 24 hours, starting with 0-1 local time
+      read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j)
+    end do
+    read(unitspecies,*)
+    do j=1,7      ! 7 days of the week, starting with Monday
+      read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j)
+    end do
+
   endif
 
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision c7d1052e9b4a8319033f957f6b77184654bca402)
+++ src/timemanager.f90	(revision a9cf4b1e763458eb1d7f20f9d55d1f5bd8a87da6)
@@ -449,5 +449,5 @@
 
         !write(*,46) float(itime)/3600,itime,numpart
-45      format(i13,' Seconds simulated: ',i13, ' Particles:    Uncertainty: ',3f7.3)
+45      format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES:    Uncertainty: ',3f7.3)
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
         if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
