source: flexpart.git/src/wetdepo.f90 @ 1c3c778

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 1c3c778 was 1c3c778, checked in by Espen Sollum ATMOS <eso@…>, 8 years ago

Added 'incloud_ratio' for wet deposition

  • Property mode set to 100644
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!*****************************************************************************
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    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
177
178! get the level were the actual particle is in
179    do il=2,nz
180      if (height(il).gt.ztra1(jpart)) then
181        hz=il-1
182!        goto 26
183        exit
184      endif
185    end do
186!26  continue
187
188    n=memind(2)
189    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
190         n=memind(1)
191
192    if (ngrid.eq.0) then
193      clouds_v=clouds(ix,jy,hz,n)
194      clouds_h=cloudsh(ix,jy,n)
195    else
196      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
197      clouds_h=cloudshn(ix,jy,n,ngrid)
198    endif
199
200! if there is no precipitation or the particle is above the clouds no
201! scavenging is done
202
203    if (clouds_v.le.1) goto 20
204
205! 1) Parameterization of the the area fraction of the grid cell where the
206!    precipitation occurs: the absolute limit is the total cloud cover, but
207!    for low precipitation rates, an even smaller fraction of the grid cell
208!    is used. Large scale precipitation occurs over larger areas than
209!    convective precipitation.
210!**************************************************************************
211
212    if (lsp.gt.20.) then
213      i=5
214    else if (lsp.gt.8.) then
215      i=4
216    else if (lsp.gt.3.) then
217      i=3
218    else if (lsp.gt.1.) then
219      i=2
220    else
221      i=1
222    endif
223
224    if (convp.gt.20.) then
225      j=5
226    else if (convp.gt.8.) then
227      j=4
228    else if (convp.gt.3.) then
229      j=3
230    else if (convp.gt.1.) then
231      j=2
232    else
233      j=1
234    endif
235
236
237!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
238! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
239! for now they are treated the same
240    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
241    grfraction(2)=max(0.05,cc*(lfr(i)))
242    grfraction(3)=max(0.05,cc*(cfr(j)))
243
244
245! 2) Computation of precipitation rate in sub-grid cell
246!******************************************************
247    prec(1)=(lsp+convp)/grfraction(1)
248    prec(2)=(lsp)/grfraction(2)
249    prec(3)=(convp)/grfraction(3)
250
251
252! 3) Computation of scavenging coefficients for all species
253!    Computation of wet deposition
254!**********************************************************
255
256    do ks=1,nspec      ! loop over species
257      wetdeposit(ks)=0.
258      wetscav=0.   
259
260!ZHG test if it nested?
261      if (ngrid.gt.0) then
262        act_temp=ttn(ix,jy,hz,n,ngrid)
263      else
264        act_temp=tt(ix,jy,hz,n)
265      endif
266
267
268!***********************
269! BELOW CLOUD SCAVENGING
270!*********************** 
271      if (clouds_v.ge.4) then !below cloud
272
273        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then !if positive below-cloud parameters given in SPECIES file (either A or B)
274          blc_count=blc_count+1
275
276!GAS
277          if (dquer(ks) .le. 0.) then  !is gas
278            wetscav=weta(ks)*prec(1)**wetb(ks)
279
280!AEROSOL
281          else !is particle
282!NIK 17.02.2015
283! For the calculation here particle size needs to be in meter and not um as dquer is changed to in readreleases
284! for particles larger than 10 um use the largest size defined in the parameterizations (10um)
285            dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
286            if (act_temp .ge. 273. .and. weta(ks).gt.0.)  then !Rain
287! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
288! the below-cloud scavenging (rain efficienty)
289! parameter A (=weta) from SPECIES file
290              wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* &
291                   (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5))
292
293            elseif (act_temp .lt. 273. .and. wetb(ks).gt.0.)  then ! Snow
294! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
295! the below-cloud scavenging (Snow efficiency)
296! parameter B (=wetb) from SPECIES file
297              wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* &
298                   (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
299
300            endif
301
302!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
303
304          endif !gas or particle
305        endif ! positive below-cloud scavenging parameters given in Species file
306      endif !end BELOW
307
308
309!********************
310! IN CLOUD SCAVENGING
311!********************
312      if (clouds_v.lt.4) then ! In-cloud
313! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file
314        if (weta_in(ks).gt.0. .or. wetb_in(ks).gt.0.) then
315          inc_count=inc_count+1
316! if negative coefficients (turned off) set to zero for use in equation
317          if (weta_in(ks).lt.0.) weta_in(ks)=0.
318          if (wetb_in(ks).lt.0.) wetb_in(ks)=0.
319
320!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
321! nested fields
322          if (ngrid.gt.0.and.readclouds_this_nest) then
323            cl=clw4n(ix,jy,n,ngrid)*(grfraction(1)/cc)
324          else if (ngrid.eq.0.and.readclouds) then
325            cl=clw4(ix,jy,n)*(grfraction(1)/cc)
326          else                                  !parameterize cloudwater m2/m3
327!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
328            cl=1.6E-6*prec(1)**0.36
329          endif
330
331!ZHG: Calculate the partition between liquid and water phase water.
332          if (act_temp .le. 253.) then
333            liq_frac=0
334          else if (act_temp .ge. 273.) then
335            liq_frac=1
336          else
337            liq_frac =((act_temp-273.)/(273.-253.))**2.
338          end if
339! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
340          frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
341
342!ZHG Use the activated fraction and the liqid water to calculate the washout
343
344! AEROSOL
345!**************************************************
346          if (dquer(ks).gt. 0.) then ! is particle
347
348            S_i= frac_act/cl
349
350!*********************
351! GAS
352          else ! is gas
353
354            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
355!REPLACE to switch old/ new scheme
356! S_i=frac_act/cle
357            S_i=1/cle
358          endif ! gas or particle
359
360! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
361!OLD
362          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
363            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
364          else
365            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
366          endif
367
368
369        endif ! positive in-cloud scavenging parameters given in Species file
370      endif !incloud
371!END ZHG TEST
372
373!**************************************************
374! CALCULATE DEPOSITION
375!**************************************************
376!     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
377!     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
378
379      if (wetscav.gt.0.) then
380        wetdeposit(ks)=xmass1(jpart,ks)* &
381             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
382!write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &
383!             (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v
384
385
386      else ! if no scavenging
387        wetdeposit(ks)=0.
388      endif
389
390      restmass = xmass1(jpart,ks)-wetdeposit(ks)
391      if (ioutputforeachrelease.eq.1) then
392        kp=npoint(jpart)
393      else
394        kp=1
395      endif
396      if (restmass .gt. smallnum) then
397        xmass1(jpart,ks)=restmass
398!   depostatistic
399!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
400!   depostatistic
401      else
402        xmass1(jpart,ks)=0.
403      endif
404!   Correct deposited mass to the last time step when radioactive decay of
405!   gridded deposited mass was calculated
406      if (decay(ks).gt.0.) then
407        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
408      endif
409
410
411    end do !all species
412
413! Sabine Eckhardt, 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)),nage,kp)
422    endif
423
42420  continue
425  end do ! all particles
426
427! count the total number of below-cloud and in-cloud occurences:
428  tot_blc_count=tot_blc_count+blc_count
429  tot_inc_count=tot_inc_count+inc_count
430
431end subroutine wetdepo
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG