source: trunk/src/wetdepo.f90

Last change on this file was 24, checked in by igpis, 10 years ago

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

File size: 12.4 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 wetdepo(itime,ltsample,loutnext)
23  !                  i      i        i
24  !*****************************************************************************
25  !                                                                            *
26  ! Calculation of wet deposition using the concept of scavenging coefficients.*
27  ! For lack of detailed information, washout and rainout are jointly treated. *
28  ! It is assumed that precipitation does not occur uniformly within the whole *
29  ! grid cell, but that only a fraction of the grid cell experiences rainfall. *
30  ! This fraction is parameterized from total cloud cover and rates of large   *
31  ! scale and convective precipitation.                                        *
32  !                                                                            *
33  !    Author: A. Stohl                                                        *
34  !                                                                            *
35  !    1 December 1996                                                         *
36  !                                                                            *
37  ! Correction by Petra Seibert, Sept 2002:                                    *
38  ! use centred precipitation data for integration                             *
39  ! Code may not be correct for decay of deposition!                           *
40  !                                                                            *
41  !*****************************************************************************
42  !                                                                            *
43  ! Variables:                                                                 *
44  ! cc [0-1]           total cloud cover                                       *
45  ! convp [mm/h]       convective precipitation rate                           *
46  ! grfraction [0-1]   fraction of grid, for which precipitation occurs        *
47  ! ix,jy              indices of output grid cell for each particle           *
48  ! itime [s]          actual simulation time [s]                              *
49  ! jpart              particle index                                          *
50  ! ldeltat [s]        interval since radioactive decay was computed           *
51  ! lfr, cfr           area fraction covered by precipitation for large scale  *
52  !                    and convective precipitation (dependent on prec. rate)  *
53  ! loutnext [s]       time for which gridded deposition is next output        *
54  ! loutstep [s]       interval at which gridded deposition is output          *
55  ! lsp [mm/h]         large scale precipitation rate                          *
56  ! ltsample [s]       interval over which mass is deposited                   *
57  ! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
58  ! wetdeposit         mass that is wet deposited                              *
59  ! wetgrid            accumulated deposited mass on output grid               *
60  ! wetscav            scavenging coefficient                                  *
61  !                                                                            *
62  ! Constants:                                                                 *
63  !                                                                            *
64  !*****************************************************************************
65
66  use point_mod
67  use par_mod
68  use com_mod
69
70  implicit none
71
72  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
73  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
74  integer :: ks, kp
75  real :: S_i, act_temp, cl, cle ! in cloud scavenging
76  real :: clouds_h ! cloud height for the specific grid point
77  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
78  real :: wetdeposit(maxspec),restmass
79  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
80  save lfr,cfr
81
82
83  real :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
84  real :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
85
86  ! Compute interval since radioactive decay of deposited mass was computed
87  !************************************************************************
88
89  if (itime.le.loutnext) then
90    ldeltat=itime-(loutnext-loutstep)
91  else                                  ! first half of next interval
92    ldeltat=itime-loutnext
93  endif
94
95
96  ! Loop over all particles
97  !************************
98
99  do jpart=1,numpart
100
101    if (itra1(jpart).eq.-999999999) goto 20
102    if(ldirect.eq.1)then
103      if (itra1(jpart).gt.itime) goto 20
104    else
105      if (itra1(jpart).lt.itime) goto 20
106    endif
107  ! Determine age class of the particle
108    itage=abs(itra1(jpart)-itramem(jpart))
109    do nage=1,nageclass
110      if (itage.lt.lage(nage)) goto 33
111    end do
11233  continue
113
114
115  ! Determine which nesting level to be used
116  !*****************************************
117
118    ngrid=0
119    do j=numbnests,1,-1
120      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
121      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
122        ngrid=j
123        goto 23
124      endif
125    end do
12623  continue
127
128
129  ! Determine nested grid coordinates
130  !**********************************
131
132    if (ngrid.gt.0) then
133      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
134      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
135      ix=int(xtn)
136      jy=int(ytn)
137    else
138      ix=int(xtra1(jpart))
139      jy=int(ytra1(jpart))
140    endif
141
142
143  ! Interpolate large scale precipitation, convective precipitation and
144  ! total cloud cover
145  ! Note that interpolated time refers to itime-0.5*ltsample [PS]
146  !********************************************************************
147    interp_time=nint(itime-0.5*ltsample)
148
149    if (ngrid.eq.0) then
150      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
151      1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
152      memtime(1),memtime(2),interp_time,lsp,convp,cc)
153    else
154      call interpol_rain_nests(lsprecn,convprecn,tccn, &
155      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
156      memtime(1),memtime(2),interp_time,lsp,convp,cc)
157    endif
158
159    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
160  ! get the level were the actual particle is in
161    do il=2,nz
162      if (height(il).gt.ztra1(jpart)) then
163        hz=il-1
164        goto 26
165      endif
166    end do
16726  continue
168
169    n=memind(2)
170    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
171    n=memind(1)
172
173  ! if there is no precipitation or the particle is above the clouds no
174  ! scavenging is done
175
176    if (ngrid.eq.0) then
177      clouds_v=clouds(ix,jy,hz,n)
178      clouds_h=cloudsh(ix,jy,n)
179    else
180      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
181      clouds_h=cloudsnh(ix,jy,n,ngrid)
182    endif
183  !write(*,*) 'there is
184  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
185    if (clouds_v.le.1) goto 20
186  !write (*,*) 'there is scavenging'
187
188  ! 1) Parameterization of the the area fraction of the grid cell where the
189  !    precipitation occurs: the absolute limit is the total cloud cover, but
190  !    for low precipitation rates, an even smaller fraction of the grid cell
191  !    is used. Large scale precipitation occurs over larger areas than
192  !    convective precipitation.
193  !**************************************************************************
194
195    if (lsp.gt.20.) then
196      i=5
197    else if (lsp.gt.8.) then
198      i=4
199    else if (lsp.gt.3.) then
200      i=3
201    else if (lsp.gt.1.) then
202      i=2
203    else
204      i=1
205    endif
206
207    if (convp.gt.20.) then
208      j=5
209    else if (convp.gt.8.) then
210      j=4
211    else if (convp.gt.3.) then
212      j=3
213    else if (convp.gt.1.) then
214      j=2
215    else
216      j=1
217    endif
218
219    grfraction=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
220
221  ! 2) Computation of precipitation rate in sub-grid cell
222  !******************************************************
223
224    prec=(lsp+convp)/grfraction
225
226  ! 3) Computation of scavenging coefficients for all species
227  !    Computation of wet deposition
228  !**********************************************************
229
230    do ks=1,nspec      ! loop over species
231      wetdeposit(ks)=0.
232      wetscav=0.   
233
234  ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests
235
236
237  !****i*******************
238  ! BELOW CLOUD SCAVENGING
239  !*********************** 
240      if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime
241 !       for aerosols and not highliy soluble substances weta=5E-6
242        wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient
243      endif !end below-cloud
244
245
246  !********************
247  ! IN CLOUD SCAVENGING
248  !********************
249      if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime
250
251        if (ngrid.gt.0) then
252          act_temp=ttn(ix,jy,hz,n,ngrid)
253        else
254          act_temp=tt(ix,jy,hz,n)
255        endif
256
257! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
258! weta_in=2.0E-07 (default)
259! wetb_in=0.36 (default)
260! wetc_in=0.9 (default)
261! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
262        cl=weta_in(ks)*prec**wetb_in(ks)
263        if (dquer(ks).gt.0) then ! is particle
264          S_i=wetc_in(ks)/cl
265        else ! is gas
266          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
267          S_i=1/cle
268        endif
269        wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
270  !     write(*,*) 'in. wetscav:'
271  !     +  ,wetscav,cle,cl,act_temp,prec,clouds_h
272       
273      endif ! end in-cloud
274
275
276     
277  !**************************************************
278  ! CALCULATE DEPOSITION
279  !**************************************************
280  !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
281  !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
282
283      if (wetscav.gt.0.) then
284        wetdeposit(ks)=xmass1(jpart,ks)* &
285        (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
286      else ! if no scavenging
287        wetdeposit(ks)=0.
288      endif
289
290
291      restmass = xmass1(jpart,ks)-wetdeposit(ks)
292      if (ioutputforeachrelease.eq.1) then
293        kp=npoint(jpart)
294      else
295        kp=1
296      endif
297      if (restmass .gt. smallnum) then
298        xmass1(jpart,ks)=restmass
299  !   depostatistic
300  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
301  !   depostatistic
302      else
303        xmass1(jpart,ks)=0.
304      endif
305  !   Correct deposited mass to the last time step when radioactive decay of
306  !   gridded deposited mass was calculated
307      if (decay(ks).gt.0.) then
308        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
309      endif
310
311
312
313    end do !all species
314
315  ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
316  ! Add the wet deposition to accumulated amount on output grid and nested output grid
317  !*****************************************************************************
318
319    if (ldirect.eq.1) then
320      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
321      real(ytra1(jpart)),nage,kp)
322      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
323        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
324    endif
325
32620  continue
327  end do ! all particles
328
329end subroutine wetdepo
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG