source: trunk/src/getvdep.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
  • Property svn:executable set to *
File size: 7.3 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 getvdep(n,ix,jy,ust,temp,pa,L,gr,rh,rr,snow,vdepo)
23  !                   i i  i   i   i   i  i i  i  i    i o
24  !*****************************************************************************
25  !                                                                            *
26  !  This routine calculates the dry deposition velocities.                    *
27  !                                                                            *
28  !     Author: A. Stohl                                                       *
29  !                                                                            *
30  !     20 December 1996                                                       *
31  !     Sabine Eckhardt, Jan 07                                                *
32  !     if the latitude is negative: add half a year to the julian day         *
33  !                                                                            *
34  !*****************************************************************************
35  !                                                                            *
36  ! Variables:                                                                 *
37  ! gr [W/m2]         global radiation                                         *
38  ! L [m]             Obukhov length                                           *
39  ! nyl               kinematic viscosity                                      *
40  ! pa [Pa]           surface air pressure                                     *
41  ! ra [s/m]          aerodynamic resistance                                   *
42  ! raquer [s/m]      average aerodynamic resistance                           *
43  ! rh [0-1]          relative humidity                                        *
44  ! rhoa              density of the air                                       *
45  ! rr [mm/h]         precipitation rate                                       *
46  ! temp [K]          2m temperature                                           *
47  ! tc [C]            2m temperature                                           *
48  ! ust [m/s]         friction velocity                                        *
49  ! snow [m of water equivalent] snow depth                                    *
50  ! xlanduse          fractions of numclasS landuses for each model grid point *
51  !                                                                            *
52  !*****************************************************************************
53
54  use par_mod
55  use com_mod
56
57  implicit none
58
59  integer :: yyyymmdd,hhmmss,yyyy,mmdd,n,lseason,i,j,ix,jy
60  real :: vdepo(maxspec),vd,rb(maxspec),rc(maxspec),raquer,ylat
61  real :: raerod,ra,ust,temp,tc,pa,L,gr,rh,rr,myl,nyl,rhoa,diffh2o,snow
62  real :: slanduse(numclass)
63  real,parameter :: eps=1.e-5
64  real(kind=dp) :: jul
65
66  ! Calculate month and determine the seasonal category
67  !****************************************************
68
69  jul=bdate+real(wftime(n),kind=dp)/86400._dp
70
71  ylat=jy*dy+ylat0
72  if (ylat.lt.0) then
73      jul=jul+365/2
74  endif
75
76
77  call caldate(jul,yyyymmdd,hhmmss)
78  yyyy=yyyymmdd/10000
79  mmdd=yyyymmdd-10000*yyyy
80
81  if ((ylat.gt.-20).and.(ylat.lt.20)) then
82     mmdd=600 ! summer
83  endif
84
85  if ((mmdd.ge.1201).or.(mmdd.le.301)) then
86    lseason=4
87  else if ((mmdd.ge.1101).or.(mmdd.le.331)) then
88    lseason=3
89  else if ((mmdd.ge.401).and.(mmdd.le.515)) then
90    lseason=5
91  else if ((mmdd.ge.516).and.(mmdd.le.915)) then
92    lseason=1
93  else
94    lseason=2
95  endif
96
97  ! Calculate diffusivity of water vapor
98  !************************************
99  diffh2o=2.11e-5*(temp/273.15)**1.94*(101325/pa)
100
101  ! Conversion of temperature from K to C
102  !**************************************
103
104  tc=temp-273.15
105
106  ! Calculate dynamic viscosity
107  !****************************
108
109  if (tc.lt.0) then
110    myl=(1.718+0.0049*tc-1.2e-05*tc**2)*1.e-05
111  else
112    myl=(1.718+0.0049*tc)*1.e-05
113  endif
114
115  ! Calculate kinematic viscosity
116  !******************************
117
118  rhoa=pa/(287.*temp)
119  nyl=myl/rhoa
120
121
122  ! 0. Set all deposition velocities zero
123  !**************************************
124
125  do i=1,nspec
126    vdepo(i)=0.
127  end do
128
129
130  ! 1. Compute surface layer resistances rb
131  !****************************************
132
133  call getrb(nspec,ust,nyl,diffh2o,reldiff,rb)
134
135  ! change for snow
136  do j=1,numclass
137    if (snow.gt.0.001) then ! 10 mm
138       if (j.eq.12) then
139         slanduse(j)=1.
140       else
141         slanduse(j)=0.
142       endif
143    else
144       slanduse(j)=xlanduse(ix,jy,j)
145    endif
146  end do
147
148  raquer=0.
149  do j=1,numclass            ! loop over all landuse classes
150
151    if (slanduse(j).gt.eps)  then
152
153  ! 2. Calculate aerodynamic resistance ra
154  !***************************************
155
156      ra=raerod(L,ust,z0(j))
157      raquer=raquer+ra*slanduse(j)
158
159  ! 3. Calculate surface resistance for gases
160  !******************************************
161
162      call getrc(nspec,lseason,j,tc,gr,rh,rr,rc)
163
164  ! 4. Calculate deposition velocities for gases and ...
165  ! 5. ... sum deposition velocities for all landuse classes
166  !*********************************************************
167
168      do i=1,nspec
169        if (reldiff(i).gt.0.) then
170          if ((ra+rb(i)+rc(i)).gt.0.) then
171            vd=1./(ra+rb(i)+rc(i))
172  ! XXXXXXXXXXXXXXXXXXXXXXXXXX TEST
173  !         vd=1./rc(i)
174  ! XXXXXXXXXXXXXXXXXXXXXXXXXX TEST
175          else
176            vd=9.999
177          endif
178          vdepo(i)=vdepo(i)+vd*slanduse(j)
179        endif
180      end do
181    endif
182  end do
183
184
185  ! 6. Calculate deposition velocities for particles
186  !*************************************************
187
188  call partdep(nspec,density,fract,schmi,vset,raquer,ust,nyl,vdepo)
189
190
191  ! 7. If no detailed parameterization available, take constant deposition
192  !    velocity if that is available
193  !***********************************************************************
194
195  do i=1,nspec
196    if ((reldiff(i).lt.0.).and.(density(i).lt.0.).and. &
197         (dryvel(i).gt.0.)) then
198      vdepo(i)=dryvel(i)
199    endif
200  end do
201
202
203end subroutine getvdep
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG