source: trunk/src/get_settling.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 5.8 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine get_settling(itime,xt,yt,zt,nsp,settling)
23  !                          i   i  i  i   i     o
24  !*****************************************************************************
25  !                                                                            *
26  !  This subroutine calculates particle settling velocity.                    *
27  !                                                                            *
28  !    Author: A. Stohl                                                        *
29  !                                                                            *
30  !    May 2010                                                                *
31  !                                                                            *
32  !  Improvement over traditional settling calculation in FLEXPART:            *
33  !  generalize to higher Reynolds numbers and also take into account the      *
34  !  temperature dependence of dynamic viscosity.                              *
35  !                                                                            *
36  !  Based on:                                                                 *
37  !  Naeslund E., and Thaning, L. (1991): On the settling velocity in a        *
38  !  nonstationary atmosphere, Aerosol Science and Technology 14, 247-256.     *
39  !                                                                            *
40  !*****************************************************************************
41  !                                                                            *
42  ! Variables:                                                                 *
43  ! itime [s]          current temporal position                               *
44  ! xt,yt,zt           coordinates position for which wind data shall be cal-  *
45  !                    culated                                                 *
46  !                                                                            *
47  ! Constants:                                                                 *
48  !                                                                            *
49  !*****************************************************************************
50
51  use par_mod
52  use com_mod
53
54  implicit none
55
56  integer :: itime,indz
57  real :: xt,yt,zt
58
59  ! Auxiliary variables needed for interpolation
60  real :: dz1,dz2,dz
61  real :: rho1(2),tt1(2),temperature,airdens,vis_dyn,vis_kin,viscosity
62  real :: settling,settling_old,reynolds,c_d
63  integer :: i,n,nix,njy,indzh,nsp
64
65
66  !*****************************************************************************
67  ! 1. Interpolate temperature and density: nearest neighbor interpolation sufficient
68  !*****************************************************************************
69
70  nix=int(xt)
71  njy=int(yt)
72
73  ! Determine the level below the current position for u,v
74  !*******************************************************
75
76  do i=2,nz
77    if (height(i).gt.zt) then
78      indz=i-1
79      goto 6
80    endif
81  end do
826   continue
83
84
85  ! Vertical distance to the level below and above current position
86  !****************************************************************
87
88  dz=1./(height(indz+1)-height(indz))
89  dz1=(zt-height(indz))*dz
90  dz2=(height(indz+1)-zt)*dz
91
92
93  ! Bilinear horizontal interpolation
94  !**********************************
95
96  ! Loop over 2 levels
97  !*******************
98
99  do n=1,2
100    indzh=indz+n-1
101    rho1(n)=rho(nix,njy,indzh,1)
102    tt1(n)=tt(nix,njy,indzh,1)
103  end do
104
105
106  ! Linear vertical interpolation
107  !******************************
108
109  temperature=dz2*tt1(1)+dz1*tt1(2)
110  airdens=dz2*rho1(1)+dz1*rho1(2)
111
112
113  vis_dyn=viscosity(temperature)
114  vis_kin=vis_dyn/airdens
115
116  reynolds=dquer(nsp)/1.e6*abs(vsetaver(nsp))/vis_kin
117
118
119
120  ! Iteration to determine both Reynolds number and settling velocity
121  !******************************************************************
122
123  settling_old=vsetaver(nsp)    ! initialize iteration with Stokes' law, constant viscosity estimate
124
125  do i=1,20    ! do a few iterations
126
127    if (reynolds.lt.1.917) then
128      c_d=24./reynolds
129    else if (reynolds.lt.500.) then
130      c_d=18.5/(reynolds**0.6)
131    else
132      c_d=0.44
133    endif
134
135    settling=-1.* &
136         sqrt(4*ga*dquer(nsp)/1.e6*density(nsp)*cunningham(nsp)/ &
137         (3.*c_d*airdens))
138
139    if (abs((settling-settling_old)/settling).lt.0.01) goto 11  ! stop iteration
140
141    reynolds=dquer(nsp)/1.e6*abs(settling)/vis_kin
142    settling_old=settling
143  end do
144
14511   continue
146
147end subroutine get_settling
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG