source: trunk/src/wetdepo.f90 @ 20

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

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

File size: 16.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  ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
42  ! scheme, not dated
43  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
44  !*****************************************************************************
45  !                                                                            *
46  ! Variables:                                                                 *
47  ! cc [0-1]           total cloud cover                                       *
48  ! convp [mm/h]       convective precipitation rate                           *
49  ! grfraction [0-1]   fraction of grid, for which precipitation occurs        *
50  ! ix,jy              indices of output grid cell for each particle           *
51  ! itime [s]          actual simulation time [s]                              *
52  ! jpart              particle index                                          *
53  ! ldeltat [s]        interval since radioactive decay was computed           *
54  ! lfr, cfr           area fraction covered by precipitation for large scale  *
55  !                    and convective precipitation (dependent on prec. rate)  *
56  ! loutnext [s]       time for which gridded deposition is next output        *
57  ! loutstep [s]       interval at which gridded deposition is output          *
58  ! lsp [mm/h]         large scale precipitation rate                          *
59  ! ltsample [s]       interval over which mass is deposited                   *
60  ! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
61  ! wetdeposit         mass that is wet deposited                              *
62  ! wetgrid            accumulated deposited mass on output grid               *
63  ! wetscav            scavenging coefficient                                  *
64  !                                                                            *
65  ! Constants:                                                                 *
66  !                                                                            *
67  !*****************************************************************************
68
69  use point_mod
70  use par_mod
71  use com_mod
72
73  implicit none
74
75  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
76  integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
77  integer :: kz !PS scheme
78  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
79  integer :: scheme_number ! NIK==1, PS ==2
80  real :: S_i, act_temp, cl, cle ! in cloud scavenging
81  real :: clouds_h ! cloud height for the specific grid point
82  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
83  real :: wetdeposit(maxspec),restmass
84  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
85  save lfr,cfr
86
87
88  real :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
89  real :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
90
91  ! Compute interval since radioactive decay of deposited mass was computed
92  !************************************************************************
93
94  if (itime.le.loutnext) then
95    ldeltat=itime-(loutnext-loutstep)
96  else                                  ! first half of next interval
97    ldeltat=itime-loutnext
98  endif
99
100
101  ! Loop over all particles
102  !************************
103
104  do jpart=1,numpart
105    if (itra1(jpart).eq.-999999999) goto 20
106    if(ldirect.eq.1)then
107      if (itra1(jpart).gt.itime) goto 20
108    else
109      if (itra1(jpart).lt.itime) goto 20
110    endif
111  ! Determine age class of the particle
112    itage=abs(itra1(jpart)-itramem(jpart))
113    do nage=1,nageclass
114      if (itage.lt.lage(nage)) goto 33
115    end do
11633   continue
117
118
119  ! Determine which nesting level to be used
120  !*****************************************
121
122    ngrid=0
123    do j=numbnests,1,-1
124      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
125           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
126        ngrid=j
127        goto 23
128      endif
129    end do
13023   continue
131
132
133  ! Determine nested grid coordinates
134  !**********************************
135
136    if (ngrid.gt.0) then
137      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
138      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
139      ix=int(xtn)
140      jy=int(ytn)
141    else
142      ix=int(xtra1(jpart))
143      jy=int(ytra1(jpart))
144    endif
145
146
147  ! Interpolate large scale precipitation, convective precipitation and
148  ! total cloud cover
149  ! Note that interpolated time refers to itime-0.5*ltsample [PS]
150  !********************************************************************
151    interp_time=nint(itime-0.5*ltsample)
152
153!    PS nest case still needs to be implemented!!   
154!    if (ngrid.eq.0) then
155!      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
156!           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
157!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
158           call interpol_rain(lsprec,convprec,tcc,     &
159            icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
160            memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
161            memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv) 
162!    else
163!      call interpol_rain_nests(lsprecn,convprecn,tccn, &
164!           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
165!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
166!    endif
167
168
169 !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
170 !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
171        prec = lsp+convp
172        precsub = 0.01
173        if (prec .lt. precsub) then
174          goto 20
175        else
176          f = (prec-precsub)/prec
177          lsp = f*lsp
178          convp = f*convp
179        endif
180
181
182  ! get the level were the actual particle is in
183      do il=2,nz
184        if (height(il).gt.ztra1(jpart)) then
185          !hz=il-1
186          kz=il-1
187          goto 26
188        endif
189      end do
19026     continue
191
192  n=memind(2)
193  if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
194       n=memind(1)
195
196  ! if there is no precipitation or the particle is above the clouds no
197  ! scavenging is done
198
199!old scheme
200!  if (ngrid.eq.0) then
201!     clouds_v=clouds(ix,jy,hz,n)
202!     clouds_h=cloudsh(ix,jy,n)
203!  else
204!     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
205!     clouds_h=cloudsnh(ix,jy,n,ngrid)
206!  endif
207!  !write(*,*) 'there is
208!  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
209!  if (clouds_v.le.1) goto 20
210!  !write (*,*) 'there is scavenging'
211
212  ! PS: part of 2011/2012 fix
213
214        if (ztra1(jpart) .le. float(ictop)) then
215          if (ztra1(jpart) .gt. float(icbot)) then
216            indcloud = 2 ! in-cloud
217          else
218            indcloud = 1 ! below-cloud
219          endif
220        elseif (ictop .eq. icmv) then
221          indcloud = 0 ! no cloud found, use old scheme
222        else
223          goto 20 ! above cloud
224        endif
225
226
227  ! 1) Parameterization of the the area fraction of the grid cell where the
228  !    precipitation occurs: the absolute limit is the total cloud cover, but
229  !    for low precipitation rates, an even smaller fraction of the grid cell
230  !    is used. Large scale precipitation occurs over larger areas than
231  !    convective precipitation.
232  !**************************************************************************
233
234    if (lsp.gt.20.) then
235      i=5
236    else if (lsp.gt.8.) then
237      i=4
238    else if (lsp.gt.3.) then
239      i=3
240    else if (lsp.gt.1.) then
241      i=2
242    else
243      i=1
244    endif
245
246    if (convp.gt.20.) then
247      j=5
248    else if (convp.gt.8.) then
249      j=4
250    else if (convp.gt.3.) then
251      j=3
252    else if (convp.gt.1.) then
253      j=2
254    else
255      j=1
256    endif
257
258    grfraction=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
259
260  ! 2) Computation of precipitation rate in sub-grid cell
261  !******************************************************
262
263    prec=(lsp+convp)/grfraction
264
265  ! 3) Computation of scavenging coefficients for all species
266  !    Computation of wet deposition
267  !**********************************************************
268
269    do ks=1,nspec            ! loop over species
270      wetdeposit(ks)=0.
271
272   
273    !conflicting changes to the same routine: 1=NIK 2 =PS
274    scheme_number=2   
275    if (scheme_number.eq.1) then !NIK
276
277      if (weta(ks).gt.0.) then
278        if (clouds_v.ge.4) then
279  !          BELOW CLOUD SCAVENGING
280  !          for aerosols and not highliy soluble substances weta=5E-6
281          wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
282  !        write(*,*) 'bel. wetscav: ',wetscav
283
284        else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
285  !        IN CLOUD SCAVENGING
286  ! BUGFIX tt for nested fields should be ttn
287  ! sec may 2008
288          if (ngrid.gt.0) then
289             act_temp=ttn(ix,jy,hz,n,ngrid)
290          else
291             act_temp=tt(ix,jy,hz,n)
292          endif
293
294! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
295! weta_in=2.0E-07 (default)
296! wetb_in=0.36 (default)
297! wetc_in=0.9 (default)
298! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
299          cl=weta_in(ks)*prec**wetb_in(ks)
300          if (dquer(ks).gt.0) then ! is particle
301            S_i=wetc_in(ks)/cl
302           else ! is gas
303            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
304            S_i=1/cle
305           endif
306           wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
307  !         write(*,*) 'in. wetscav:'
308  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
309        endif
310
311
312  !      if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
313  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
314        wetdeposit(ks)=xmass1(jpart,ks)* &
315             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
316  !      new particle mass:
317  !      if (wetdeposit(ks).gt.0) then
318  !         write(*,*) 'wetdepo: ',wetdeposit(ks),ks
319  !      endif
320        restmass = xmass1(jpart,ks)-wetdeposit(ks)
321         if (ioutputforeachrelease.eq.1) then
322           kp=npoint(jpart)
323         else
324           kp=1
325         endif
326        if (restmass .gt. smallnum) then
327          xmass1(jpart,ks)=restmass
328  !ccccccccccccccc depostatistic
329  !       wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
330  !ccccccccccccccc depostatistic
331        else
332          xmass1(jpart,ks)=0.
333        endif
334  ! Correct deposited mass to the last time step when radioactive decay of
335  ! gridded deposited mass was calculated
336       if (decay(ks).gt.0.) then
337          wetdeposit(ks)=wetdeposit(ks) &
338               *exp(abs(ldeltat)*decay(ks))
339       endif
340      else  ! weta(k)
341         wetdeposit(ks)=0.
342      endif ! weta(k)
343
344     elseif (scheme_number.eq.2) then ! PS
345
346!PS          indcloud=0 ! Use this for FOR TESTING,
347!PS          will skip the new in/below cloud method !!!             
348
349          if (weta(ks).gt.0.) then
350            if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
351!C               for aerosols and not highliy soluble substances weta=5E-6
352              wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
353!c             write(*,*) 'bel. wetscav: ',wetscav
354            elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
355              if (ngrid.gt.0) then
356                 act_temp=ttn(ix,jy,kz,n,ngrid)
357              else
358                 act_temp=tt(ix,jy,kz,n)
359              endif
360
361! from NIK             
362! weta_in=2.0E-07 (default)
363! wetb_in=0.36 (default)
364! wetc_in=0.9 (default)
365
366
367              cl=2E-7*prec**0.36
368              if (dquer(ks).gt.0) then ! is particle
369                S_i=0.9/cl
370              else ! is gas
371                cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
372                S_i=1/cle
373              endif
374              wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
375            else ! PS: no cloud diagnosed, old scheme,
376!CPS          using with fixed a,b for simplicity, one may wish to change!!
377              wetscav = 1.e-4*prec**0.62
378            endif
379
380
381            wetdeposit(ks)=xmass1(jpart,ks)*    &
382            ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
383             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS)
384            restmass = xmass1(jpart,ks)-wetdeposit(ks)
385            if (ioutputforeachrelease.eq.1) then
386              kp=npoint(jpart)
387            else
388              kp=1
389            endif
390            if (restmass .gt. smallnum) then
391              xmass1(jpart,ks)=restmass
392!cccccccccccccccc depostatistic
393!c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
394!cccccccccccccccc depostatistic
395            else
396              xmass1(jpart,ks)=0.
397            endif
398!C Correct deposited mass to the last time step when radioactive decay of
399!C gridded deposited mass was calculated
400            if (decay(ks).gt.0.) then
401              wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
402            endif
403          else  ! weta(k)<0
404             wetdeposit(ks)=0.
405          endif
406
407     endif !on scheme
408
409
410
411    end do
412
413  ! Sabine Eckhard, June 2008 create deposition runs only for forward runs
414  ! Add the wet deposition to accumulated amount on output grid and nested output grid
415  !*****************************************************************************
416
417!    if (ldirect.eq.1) then
418!    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
419!         real(ytra1(jpart)),nage,kp)
420!    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
421!         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
422!         nage,kp)
423!    endif
424
425   !PS
426        if (ldirect.eq.1) then
427          call wetdepokernel(nclass(jpart),wetdeposit,     &
428             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
429          if (nested_output.eq.1)                          &
430           call wetdepokernel_nest(nclass(jpart),wetdeposit, &
431             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 
432        endif
433
434
43520  continue
436  end do
437
438end subroutine wetdepo
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG