Index: src/get_settling.f90
===================================================================
--- src/get_settling.f90	(revision 67b13e9201ae6e3332b2fa2fbd119f57f68f525e)
+++ src/get_settling.f90	(revision 67b13e9201ae6e3332b2fa2fbd119f57f68f525e)
@@ -0,0 +1,147 @@
+!**********************************************************************
+! 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 get_settling(itime,xt,yt,zt,nsp,settling)
+  !                          i   i  i  i   i     o
+  !*****************************************************************************
+  !                                                                            *
+  !  This subroutine calculates particle settling velocity.                    *
+  !                                                                            *
+  !    Author: A. Stohl                                                        *
+  !                                                                            *
+  !    May 2010                                                                *
+  !                                                                            *
+  !  Improvement over traditional settling calculation in FLEXPART:            *
+  !  generalize to higher Reynolds numbers and also take into account the      *
+  !  temperature dependence of dynamic viscosity.                              *
+  !                                                                            *
+  !  Based on:                                                                 *
+  !  Naeslund E., and Thaning, L. (1991): On the settling velocity in a        *
+  !  nonstationary atmosphere, Aerosol Science and Technology 14, 247-256.     *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! itime [s]          current temporal position                               *
+  ! xt,yt,zt           coordinates position for which wind data shall be cal-  *
+  !                    culated                                                 *
+  !                                                                            *
+  ! Constants:                                                                 *
+  !                                                                            *
+  !*****************************************************************************
+
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  integer :: itime,indz
+  real :: xt,yt,zt
+
+  ! Auxiliary variables needed for interpolation
+  real :: dz1,dz2,dz
+  real :: rho1(2),tt1(2),temperature,airdens,vis_dyn,vis_kin,viscosity
+  real :: settling,settling_old,reynolds,c_d
+  integer :: i,n,nix,njy,indzh,nsp
+
+
+  !*****************************************************************************
+  ! 1. Interpolate temperature and density: nearest neighbor interpolation sufficient
+  !*****************************************************************************
+
+  nix=int(xt)
+  njy=int(yt)
+
+  ! Determine the level below the current position for u,v
+  !*******************************************************
+
+  do i=2,nz
+    if (height(i).gt.zt) then
+      indz=i-1
+      goto 6
+    endif
+  end do
+6   continue
+
+
+  ! Vertical distance to the level below and above current position
+  !****************************************************************
+
+  dz=1./(height(indz+1)-height(indz))
+  dz1=(zt-height(indz))*dz
+  dz2=(height(indz+1)-zt)*dz
+
+
+  ! Bilinear horizontal interpolation
+  !**********************************
+
+  ! Loop over 2 levels
+  !*******************
+
+  do n=1,2
+    indzh=indz+n-1
+    rho1(n)=rho(nix,njy,indzh,1)
+    tt1(n)=tt(nix,njy,indzh,1)
+  end do
+
+
+  ! Linear vertical interpolation
+  !******************************
+
+  temperature=dz2*tt1(1)+dz1*tt1(2)
+  airdens=dz2*rho1(1)+dz1*rho1(2)
+
+
+  vis_dyn=viscosity(temperature)
+  vis_kin=vis_dyn/airdens
+
+  reynolds=dquer(nsp)/1.e6*abs(vsetaver(nsp))/vis_kin
+
+
+
+  ! Iteration to determine both Reynolds number and settling velocity
+  !******************************************************************
+
+  settling_old=vsetaver(nsp)    ! initialize iteration with Stokes' law, constant viscosity estimate
+
+  do i=1,20    ! do a few iterations
+
+    if (reynolds.lt.1.917) then
+      c_d=24./reynolds
+    else if (reynolds.lt.500.) then
+      c_d=18.5/(reynolds**0.6)
+    else
+      c_d=0.44
+    endif
+
+    settling=-1.* &
+         sqrt(4*ga*dquer(nsp)/1.e6*density(nsp)*cunningham(nsp)/ &
+         (3.*c_d*airdens))
+
+    if (abs((settling-settling_old)/settling).lt.0.01) goto 11  ! stop iteration
+
+    reynolds=dquer(nsp)/1.e6*abs(settling)/vis_kin
+    settling_old=settling
+  end do
+
+11   continue
+
+end subroutine get_settling
