Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
+++ src/FLEXPART_MPI.f90	(revision 5184a7c7fc4985dfd336c4fd7fb8acacc26caef4)
@@ -80,5 +80,4 @@
   ! FLEXPART version string
   flexversion_major = '10' ! Major version number, also used for species file names
-!  flexversion='Ver. 10 Beta MPI (2015-05-01)'
   flexversion='Ver. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)'
   verbosity=0
Index: src/advance.f90
===================================================================
--- src/advance.f90	(revision e52967c0a5f99e764feb376b6267ce20286a2850)
+++ src/advance.f90	(revision 5184a7c7fc4985dfd336c4fd7fb8acacc26caef4)
@@ -535,14 +535,18 @@
 
       if (mdomainfill.eq.0) then
-! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
-!              particle number and thus xmass array goes out of bounds
-!        do nsp=1,nspec
-!           if (xmass(nrelpoint,nsp).gt.eps2) goto 887
-!         end do
-! 887     nsp=min(nsp,nspec)
-        if (nspec.eq.1.and.density(1).gt.0.) then
-          call get_settling(itime,real(xt),real(yt),zt,nspec,settling)  !bugfix
+        if (lsettling) then
+          do nsp=1,nspec
+            if (xmass(nrelpoint,nsp).gt.eps2) exit
+          end do
+          if (nsp.gt.nspec) then
+  ! This should never happen          
+            write(*,*) 'advance.f90: ERROR: could not find releasepoint'
+            stop
+          end if
+          if (density(nsp).gt.0.) then
+            call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
+            w=w+settling
+          end if
         end if
-        w=w+settling
       endif
 
@@ -700,20 +704,21 @@
   !*************************************************************************
 
-
-
-    if (mdomainfill.eq.0) then
-! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
-!              particle number and thus xmass array goes out of bounds
-
-!      do nsp=1,nspec
-!         if (xmass(nrelpoint,nsp).gt.eps2) goto 888
-!       end do
-! 888   nsp=min(nsp,nspec)
-!        if (density(nsp).gt.0.) then
-      if (nspec.eq.1.and.density(1).gt.0.) then
-        call get_settling(itime,real(xt),real(yt),zt,nspec,settling)  !bugfix
+  if (mdomainfill.eq.0) then
+    if (lsettling) then
+      do nsp=1,nspec
+        if (xmass(nrelpoint,nsp).gt.eps2) exit
+      end do
+      if (nsp.gt.nspec) then
+  ! This should never happen          
+        write(*,*) 'advance.f90: ERROR: could not find releasepoint'
+        stop
       end if
-      w=w+settling
+      if (density(nsp).gt.0.) then
+        call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
+        w=w+settling
+      end if
     endif
+  end if
+
 
   ! Calculate position at time step itime+lsynctime
@@ -910,16 +915,19 @@
 
   if (mdomainfill.eq.0) then
-! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
-!              particle number and thus xmass array goes out of bounds
-!    do nsp=1,nspec
-!       if (xmass(nrelpoint,nsp).gt.eps2) goto 889
-!     end do
-! 889   nsp=min(nsp,nspec)
-!      if (density(nsp).gt.0.) then
-    if (nspec.eq.1.and.density(1).gt.0.) then
-      call get_settling(itime+ldt,real(xt),real(yt),zt,nspec,settling)  !bugfix
-    end if
-    w=w+settling
-  endif
+    if (lsettling) then
+      do nsp=1,nspec
+        if (xmass(nrelpoint,nsp).gt.eps2) exit
+      end do
+      if (nsp.gt.nspec) then
+  ! This should never happen          
+        write(*,*) 'advance.f90: ERROR: could not find releasepoint'
+        stop
+      end if
+      if (density(nsp).gt.0.) then
+        call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling) !bugfix
+        w=w+settling
+      end if
+    endif
+  end if
 
 
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 2bec33eff0e52c1452ff0940aa3db5c57de9aa74)
+++ src/com_mod.f90	(revision 5184a7c7fc4985dfd336c4fd7fb8acacc26caef4)
@@ -135,4 +135,7 @@
 !ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
   logical :: sumclouds=.false.
+
+!ESO: Disable settling if more than 1 species per release point
+  logical :: lsettling=.true.
 
   logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest
Index: src/readreleases.f90
===================================================================
--- src/readreleases.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
+++ src/readreleases.f90	(revision 5184a7c7fc4985dfd336c4fd7fb8acacc26caef4)
@@ -74,8 +74,9 @@
   implicit none
 
-  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
+  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat,irel,ispc,nsettle
   integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
   real(kind=dp) :: jul1,jul2,julm,juldate
+  real,parameter :: eps2=1.e-9
   character(len=50) :: line
   logical :: old
@@ -93,17 +94,17 @@
   ! declare namelists
   namelist /releases_ctrl/ &
-    nspec, &
-    specnum_rel
+       nspec, &
+       specnum_rel
 
   namelist /release/ &
-    idate1, itime1, &
-    idate2, itime2, &
-    lon1, lon2, &
-    lat1, lat2, &
-    z1, z2, &
-    zkind, &
-    mass, &
-    parts, &
-    comment
+       idate1, itime1, &
+       idate2, itime2, &
+       lon1, lon2, &
+       lat1, lat2, &
+       z1, z2, &
+       zkind, &
+       mass, &
+       parts, &
+       comment
 
   numpoint=0
@@ -132,5 +133,5 @@
   if ((readerror.ne.0).or.(nspec.lt.0)) then
 
-    ! no namelist format, close file and allow reopening in old format
+  ! no namelist format, close file and allow reopening in old format
     close(unitreleases)
     open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
@@ -138,8 +139,8 @@
     readerror=1 ! indicates old format
 
-    ! Check the format of the RELEASES file (either in free format,
-    ! or using a formatted mask)
-    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
-    !**************************************************************************
+  ! Check the format of the RELEASES file (either in free format,
+  ! or using a formatted mask)
+  ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
+  !**************************************************************************
 
     call skplin(12,unitreleases)
@@ -154,6 +155,6 @@
 
 
-    ! Skip first 11 lines (file header)
-    !**********************************
+  ! Skip first 11 lines (file header)
+  !**********************************
 
     call skplin(11,unitreleases)
@@ -193,5 +194,5 @@
       if (old) call skplin(2,unitreleases)
     end do
-    !save compoint only for the first 1000 release points
+  !save compoint only for the first 1000 release points
     read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
     if (old) call skplin(1,unitreleases)
@@ -228,6 +229,6 @@
   specnum_rel2=specnum_rel(1:nspec)
   deallocate(specnum_rel) 
-! eso: BUG, crashes here for nspec=12 and maxspec=6,
-! TODO: catch error and exit
+  ! eso: BUG, crashes here for nspec=12 and maxspec=6,
+  ! TODO: catch error and exit
   allocate(specnum_rel(nspec),stat=stat)
   if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
@@ -280,11 +281,11 @@
 
   if (readerror.ne.0) then
-    ! Skip first 11 lines (file header)
-    !**********************************
+  ! Skip first 11 lines (file header)
+  !**********************************
 
     call skplin(11,unitreleases)
 
-    ! Assign species-specific parameters needed for physical processes
-    !*************************************************************************
+  ! Assign species-specific parameters needed for physical processes
+  !*************************************************************************
 
     read(unitreleases,*,err=998) nspec
@@ -307,17 +308,17 @@
     endif
 
-    ! For backward runs, only 1 species is allowed
-    !*********************************************
-
-    !if ((ldirect.lt.0).and.(nspec.gt.1)) then
-    !write(*,*) '#####################################################'
-    !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
-    !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
-    !write(*,*) '#####################################################'
-    !  stop
-    !endif
-
-    ! Molecular weight
-    !*****************
+  ! For backward runs, only 1 species is allowed
+  !*********************************************
+
+  !if ((ldirect.lt.0).and.(nspec.gt.1)) then
+  !write(*,*) '#####################################################'
+  !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
+  !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
+  !write(*,*) '#####################################################'
+  !  stop
+  !endif
+
+  ! Molecular weight
+  !*****************
 
     if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
@@ -329,17 +330,17 @@
     endif
 
-    ! Radioactive decay
-    !******************
+  ! Radioactive decay
+  !******************
 
     decay(i)=0.693147/decay(i) !conversion half life to decay constant
 
 
-    ! Dry deposition of gases
-    !************************
+  ! Dry deposition of gases
+  !************************
 
     if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
 
-    ! Dry deposition of particles
-    !****************************
+  ! Dry deposition of particles
+  !****************************
 
     vsetaver(i)=0.
@@ -358,6 +359,6 @@
     endif
 
-    ! Dry deposition for constant deposition velocity
-    !************************************************
+  ! Dry deposition for constant deposition velocity
+  !************************************************
 
     dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
@@ -373,11 +374,11 @@
       if (lroot) then
         write (*,*) '  Below-cloud scavenging: ON'
-!  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
+  !  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
       end if
     else
       if (lroot) write (*,*) '  Below-cloud scavenging: OFF'
     endif
-    
-! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
+
+  ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
     if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
       WETDEP=.true.
@@ -385,6 +386,6 @@
       if (lroot) then
         write (*,*) '  In-cloud scavenging: ON'
-!        write (*,*) 'In-cloud scavenging coefficients: ',&
-!           &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
+  !        write (*,*) 'In-cloud scavenging coefficients: ',&
+  !           &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
       end if
     else
@@ -412,5 +413,5 @@
   numpartmax=0
   releaserate=0.
-101   numpoint=numpoint+1
+101 numpoint=numpoint+1
 
   if (readerror.lt.1) then ! reading namelist format
@@ -439,5 +440,5 @@
     compoint(min(1001,numpoint))=comment
 
-    ! namelist output
+  ! namelist output
     if (nmlout.and.lroot) then
       write(unitreleasesout,nml=release)
@@ -472,5 +473,5 @@
       mass(i)=xmass(numpoint,i)
     end do
-    !save compoint only for the first 1000 release points
+  !save compoint only for the first 1000 release points
     if (numpoint.le.1000) then
       read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
@@ -482,5 +483,5 @@
     if (old) call skplin(1,unitreleases)
 
-    ! namelist output
+  ! namelist output
     if (nmlout.and.lroot) then
       idate1=id1
@@ -524,12 +525,12 @@
   !*********************************************************************
   if (WETBKDEP) then
-      zpoint1(numpoint)=0.
-      zpoint2(numpoint)=20000.
-      kindz(numpoint)=1
+    zpoint1(numpoint)=0.
+    zpoint2(numpoint)=20000.
+    kindz(numpoint)=1
   endif
   if (DRYBKDEP) then
-      zpoint1(numpoint)=0.
-      zpoint2(numpoint)=2.*href
-      kindz(numpoint)=1
+    zpoint1(numpoint)=0.
+    zpoint2(numpoint)=2.*href
+    kindz(numpoint)=1
   endif
 
@@ -538,12 +539,12 @@
   !*********************************************************************
 
-   if (xpoint1(numpoint).lt.xlon0) &
-        xpoint1(numpoint)=xpoint1(numpoint)+360.
-   if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
-        xpoint1(numpoint)=xpoint1(numpoint)-360.
-   if (xpoint2(numpoint).lt.xlon0) &
-        xpoint2(numpoint)=xpoint2(numpoint)+360.
-   if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
-        xpoint2(numpoint)=xpoint2(numpoint)-360.
+  if (xpoint1(numpoint).lt.xlon0) &
+       xpoint1(numpoint)=xpoint1(numpoint)+360.
+  if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
+       xpoint1(numpoint)=xpoint1(numpoint)-360.
+  if (xpoint2(numpoint).lt.xlon0) &
+       xpoint2(numpoint)=xpoint2(numpoint)+360.
+  if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
+       xpoint2(numpoint)=xpoint2(numpoint)-360.
 
   ! Determine relative beginning and ending times of particle release
@@ -607,5 +608,5 @@
   goto 101
 
-250   close(unitreleases)
+250 close(unitreleases)
 
   if (nmlout.and.lroot) then
@@ -623,4 +624,27 @@
   endif
 
+  ! Disable settling if more than 1 species at any release point
+  ! or if MQUASILAG and more than one species
+  !*************************************************************
+
+  if (mquasilag.ne.0) then
+    if (nspec.gt.1) lsettling=.false.
+  else
+    do irel=1,numpoint
+      nsettle=0
+      do ispc=1,nspec
+        if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1
+      end do
+      if (nsettle.gt.1) lsettling=.false.
+    end do
+  end if
+
+  if (lroot) then
+    if (.not.lsettling) then 
+      write(*,*) 'WARNING: more than 1 species per release point, settling &
+           &disabled'
+    end if
+  end if
+
   ! Check, whether the total number of particles may exceed totally allowed
   ! number of particles at some time during the simulation
@@ -630,16 +654,16 @@
        0.99*real(maxpart)/real(lage(nageclass))) then
     if (numpartmax.gt.maxpart.and.lroot) then
-  write(*,*) '#####################################################'
-  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
-  write(*,*) '####                                             ####'
-  write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
-  write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
-  write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED   ####'
-  write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
-  write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
-  write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP          ####'
-  write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE  ####'
-  write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT.   ####'
-  write(*,*) '#####################################################'
+      write(*,*) '#####################################################'
+      write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
+      write(*,*) '####                                             ####'
+      write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
+      write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
+      write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED   ####'
+      write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
+      write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
+      write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP          ####'
+      write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE  ####'
+      write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT.   ####'
+      write(*,*) '#####################################################'
       write(*,*) 'Maximum release rate may be: ',releaserate, &
            ' particles per second'
@@ -653,4 +677,5 @@
   endif
 
+
   if (lroot) then
     write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
@@ -659,5 +684,5 @@
   return
 
-994   write(*,*) '#####################################################'
+994 write(*,*) '#####################################################'
   write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
   write(*,*) '####                                             ####'
@@ -667,5 +692,5 @@
   stop
 
-998   write(*,*) '#####################################################'
+998 write(*,*) '#####################################################'
   write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
   write(*,*) '####                                             ####'
@@ -678,5 +703,5 @@
 
 
-999   write(*,*) '#####################################################'
+999 write(*,*) '#####################################################'
   write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
   write(*,*)
