source: flexpart.git/src/wetdepo.f90 @ 462f74b

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 462f74b was 462f74b, checked in by Sabine Eckhardt <sabine@…>, 8 years ago

first version of the backward scavenging

  • Property mode set to 100644
File size: 17.9 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
74  integer(kind=1) :: clouds_v
75  integer :: ks, kp
76!  integer :: n1,n2, icbot,ictop, indcloud !TEST
77  real :: S_i, act_temp, cl, cle ! in cloud scavenging
78  real :: clouds_h ! cloud height for the specific grid point
79  real :: xtn,ytn,lsp,convp,cc,grfraction(3),prec(3),wetscav,totprec
80  real :: wetdeposit(maxspec),restmass
81  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
82!save lfr,cfr
83
84  real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
85  real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
86
87!ZHG aerosol below-cloud scavenging removal polynomial constants for rain and snow
88  real, parameter :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
89  real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
90  real :: frac_act, liq_frac, dquer_m
91
92  integer :: blc_count, inc_count
93  real    :: Si_dummy, wetscav_dummy
94  logical :: readclouds_this_nest
95
96
97! Compute interval since radioactive decay of deposited mass was computed
98!************************************************************************
99
100  if (itime.le.loutnext) then
101    ldeltat=itime-(loutnext-loutstep)
102  else                                  ! first half of next interval
103    ldeltat=itime-loutnext
104  endif
105
106! Loop over all particles
107!************************
108
109  blc_count=0
110  inc_count=0
111
112  do jpart=1,numpart
113
114    if (itra1(jpart).eq.-999999999) goto 20
115    if(ldirect.eq.1)then
116      if (itra1(jpart).gt.itime) goto 20
117    else
118      if (itra1(jpart).lt.itime) goto 20
119    endif
120
121! Determine age class of the particle
122    itage=abs(itra1(jpart)-itramem(jpart))
123    do nage=1,nageclass
124      if (itage.lt.lage(nage)) goto 33
125    end do
12633  continue
127
128
129! Determine which nesting level to be used
130!*****************************************
131
132    ngrid=0
133    do j=numbnests,1,-1
134      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
135           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
136        ngrid=j
137        goto 23
138      endif
139    end do
14023  continue
141
142
143! Determine nested grid coordinates
144!**********************************
145    readclouds_this_nest=.false.
146
147    if (ngrid.gt.0) then
148      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
149      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
150      ix=int(xtn)
151      jy=int(ytn)
152      if (readclouds_nest(ngrid)) readclouds_this_nest=.true.
153    else
154      ix=int(xtra1(jpart))
155      jy=int(ytra1(jpart))
156    endif
157
158
159! Interpolate large scale precipitation, convective precipitation and
160! total cloud cover
161! Note that interpolated time refers to itime-0.5*ltsample [PS]
162!********************************************************************
163    interp_time=nint(itime-0.5*ltsample)
164
165    if (ngrid.eq.0) then
166      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
167           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
168           memtime(1),memtime(2),interp_time,lsp,convp,cc)
169    else
170      call interpol_rain_nests(lsprecn,convprecn,tccn, &
171           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
172           memtime(1),memtime(2),interp_time,lsp,convp,cc)
173    endif
174
175!  If total precipitation is less than 0.01 mm/h - no scavenging occurs
176!  sec this is just valid if release is over a point
177    if ((lsp.lt.0.01).and.(convp.lt.0.01)) then
178          if (SCAVDEP) then
179             do ks=1,nspec
180                if (xscav_frac1(jpart,ks).lt.0) then ! first timestep no scavenging
181                   xmass1(jpart,ks)=0.
182                   xscav_frac1(jpart,ks)=0.
183!                  write (*,*) 'paricle removed - no scavenging: ',jpart,ks
184                endif
185             end do
186          endif
187          goto 20
188    endif
189
190
191! get the level were the actual particle is in
192    do il=2,nz
193      if (height(il).gt.ztra1(jpart)) then
194        hz=il-1
195!        goto 26
196        exit
197      endif
198    end do
199!26  continue
200
201    n=memind(2)
202    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
203         n=memind(1)
204
205    if (ngrid.eq.0) then
206      clouds_v=clouds(ix,jy,hz,n)
207      clouds_h=cloudsh(ix,jy,n)
208    else
209      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
210      clouds_h=cloudshn(ix,jy,n,ngrid)
211    endif
212
213! if there is no precipitation or the particle is above the clouds no
214! scavenging is done
215
216    if (clouds_v.le.1) goto 20
217
218! 1) Parameterization of the the area fraction of the grid cell where the
219!    precipitation occurs: the absolute limit is the total cloud cover, but
220!    for low precipitation rates, an even smaller fraction of the grid cell
221!    is used. Large scale precipitation occurs over larger areas than
222!    convective precipitation.
223!**************************************************************************
224
225    if (lsp.gt.20.) then
226      i=5
227    else if (lsp.gt.8.) then
228      i=4
229    else if (lsp.gt.3.) then
230      i=3
231    else if (lsp.gt.1.) then
232      i=2
233    else
234      i=1
235    endif
236
237    if (convp.gt.20.) then
238      j=5
239    else if (convp.gt.8.) then
240      j=4
241    else if (convp.gt.3.) then
242      j=3
243    else if (convp.gt.1.) then
244      j=2
245    else
246      j=1
247    endif
248
249
250!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
251! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
252! for now they are treated the same
253    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
254    grfraction(2)=max(0.05,cc*(lfr(i)))
255    grfraction(3)=max(0.05,cc*(cfr(j)))
256
257
258! 2) Computation of precipitation rate in sub-grid cell
259!******************************************************
260    prec(1)=(lsp+convp)/grfraction(1)
261    prec(2)=(lsp)/grfraction(2)
262    prec(3)=(convp)/grfraction(3)
263
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      wetscav=0.   
272
273      if (ngrid.gt.0) then
274        act_temp=ttn(ix,jy,hz,n,ngrid)
275      else
276        act_temp=tt(ix,jy,hz,n)
277      endif
278
279
280!***********************
281! BELOW CLOUD SCAVENGING
282!*********************** 
283      if (clouds_v.ge.4) then !below cloud
284
285! For gas: if positive below-cloud parameters (A or B), and dquer<=0
286!******************************************************************
287        if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
288          !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
289          blc_count=blc_count+1
290          wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
291
292! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
293!*********************************************************************************
294        else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
295          blc_count=blc_count+1
296
297!NIK 17.02.2015
298! For the calculation here particle size needs to be in meter and not um as dquer is
299! changed in readreleases
300! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
301          dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
302
303! Rain:
304          if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.)  then
305
306! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
307! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
308            wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
309                 & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
310                 &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
311
312! Snow:
313          elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.)  then
314! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
315! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
316            wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
317                 &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
318                 &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
319
320          endif
321         
322!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
323
324        endif ! gas or particle
325!      endif ! positive below-cloud scavenging parameters given in Species file
326      endif !end BELOW
327
328!********************
329! IN CLOUD SCAVENGING
330!********************
331      if (clouds_v.lt.4) then ! In-cloud
332! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
333! given in species file, or if gas and positive Henry's constant
334        if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
335          inc_count=inc_count+1
336! if negative coefficients (turned off) set to zero for use in equation
337          if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
338          if (in_aero(ks).lt.0.) in_aero(ks)=0.
339
340!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
341! nested fields
342          if (ngrid.gt.0.and.readclouds_this_nest) then
343            cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
344          else if (ngrid.eq.0.and.readclouds) then
345            cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
346          else                                  !parameterize cloudwater m2/m3
347!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
348            cl=1.6E-6*prec(1)**0.36
349          endif
350
351!ZHG: Calculate the partition between liquid and water phase water.
352          if (act_temp .le. 253.) then
353            liq_frac=0
354          else if (act_temp .ge. 273.) then
355            liq_frac=1
356          else
357            liq_frac =((act_temp-273.)/(273.-253.))**2.
358          end if
359! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
360          frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
361
362!ZHG Use the activated fraction and the liqid water to calculate the washout
363
364! AEROSOL
365!********
366          if (dquer(ks).gt.0.) then
367            S_i= frac_act/cl
368
369! GAS
370!****
371          else
372
373            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
374!REPLACE to switch old/ new scheme
375          ! S_i=frac_act/cle
376            S_i=1/cle
377          endif ! gas or particle
378
379! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
380!OLD
381          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
382            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
383          else
384            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
385          endif
386
387
388        endif ! positive in-cloud scavenging parameters given in Species file
389      endif !incloud
390!END ZHG TEST
391
392!**************************************************
393! CALCULATE DEPOSITION
394!**************************************************
395!     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
396!     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
397
398      if (wetscav.gt.0.) then
399        wetdeposit(ks)=xmass1(jpart,ks)* &
400             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
401!write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &
402!             (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v
403
404
405      else ! if no scavenging
406        wetdeposit(ks)=0.
407      endif
408
409      restmass = xmass1(jpart,ks)-wetdeposit(ks)
410      if (ioutputforeachrelease.eq.1) then
411        kp=npoint(jpart)
412      else
413        kp=1
414      endif
415      if (restmass .gt. smallnum) then
416        xmass1(jpart,ks)=restmass
417!   depostatistic
418!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
419!   depostatistic
420      else
421        xmass1(jpart,ks)=0.
422      endif
423!   Correct deposited mass to the last time step when radioactive decay of
424!   gridded deposited mass was calculated
425      if (decay(ks).gt.0.) then
426        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
427      endif
428
429      if (SCAVDEP) then
430! the calculation of the scavenged mass shall only be done once after release
431! xscav_frac1 was initialised with a negative value
432          if (xscav_frac1(jpart,ks).lt.0) then
433             if (wetdeposit(ks).eq.0) then
434! terminate particle
435                xmass1(jpart,ks)=0.
436                xscav_frac1(jpart,ks)=0.
437             else
438                xscav_frac1(jpart,ks)=xscav_frac1(jpart,ks)*(-1.)* &
439                   wetdeposit(ks)/xmass1(jpart,ks)
440!                write (*,*) 'paricle kept: ',jpart,ks,wetdeposit(ks),xscav_frac1(jpart,ks)
441             endif
442          endif
443      endif
444
445
446    end do !all species
447
448! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
449! Add the wet deposition to accumulated amount on output grid and nested output grid
450!*****************************************************************************
451
452    if (ldirect.eq.1) then
453      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
454           real(ytra1(jpart)),nage,kp)
455      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
456           wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
457    endif
458
45920  continue
460  end do ! all particles
461
462! count the total number of below-cloud and in-cloud occurences:
463  tot_blc_count=tot_blc_count+blc_count
464  tot_inc_count=tot_inc_count+inc_count
465
466end subroutine wetdepo
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG