source: trunk/src/getvdep_nests.f90 @ 28

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