source: branches/FLEXPART_9.1.3/src/wetdepo.f90 @ 13

Last change on this file since 13 was 13, checked in by saeck, 11 years ago

update to wetdepo.f90

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