source: flexpart.git/src/get_wetscav.f90 @ 6985a98

devrelease-10univie
Last change on this file since 6985a98 was 6985a98, checked in by Sabine <sabine.eckhardt@…>, 2 years ago

compiles after merge scavenging into test dev

  • Property mode set to 100644
File size: 14.1 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 get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc_count,wetscav)
23!                          i      i        i     i   i    o           o          o       o
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! lfr, cfr           area fraction covered by precipitation for large scale  *
51!                    and convective precipitation (dependent on prec. rate)  *
52! loutnext [s]       time for which gridded deposition is next output        *
53! loutstep [s]       interval at which gridded deposition is output          *
54! lsp [mm/h]         large scale precipitation rate                          *
55! ltsample [s]       interval over which mass is deposited                   *
56! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
57! wetgrid            accumulated deposited mass on output grid               *
58! wetscav            scavenging coefficient                                  *
59!                                                                            *
60! Constants:                                                                 *
61!                                                                            *
62!*****************************************************************************
63
64  use point_mod
65  use par_mod
66  use com_mod
67
68  implicit none
69
70  integer :: jpart,itime,ltsample,loutnext,i,j,ix,jy
71  integer :: ngrid,hz,il,interp_time, n
72  integer(kind=1) :: clouds_v
73  integer :: ks, kp
74  integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count
75
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 :: 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  real    :: Si_dummy, wetscav_dummy
93  logical :: readclouds_this_nest
94
95
96   wetscav=0.
97
98! Determine which nesting level to be used
99!*****************************************
100    ngrid=0
101    do j=numbnests,1,-1
102      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
103           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
104        ngrid=j
105        goto 23
106      endif
107    end do
10823  continue
109
110
111! Determine nested grid coordinates
112!**********************************
113    readclouds_this_nest=.false.
114
115    if (ngrid.gt.0) then
116      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
117      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
118      ix=int(xtn)
119      jy=int(ytn)
120      if (readclouds_nest(ngrid)) readclouds_this_nest=.true.
121    else
122      ix=int(xtra1(jpart))
123      jy=int(ytra1(jpart))
124    endif
125
126! Interpolate large scale precipitation, convective precipitation and
127! total cloud cover
128! Note that interpolated time refers to itime-0.5*ltsample [PS]
129!********************************************************************
130    interp_time=nint(itime-0.5*ltsample)
131
132    n=memind(2)
133    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
134         n=memind(1)
135
136    if (ngrid.eq.0) then
137      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
138           1,nx,ny,n,real(xtra1(jpart)),real(ytra1(jpart)),1, &
139           memtime(1),memtime(2),interp_time,lsp,convp,cc)
140    else
141      call interpol_rain_nests(lsprecn,convprecn,tccn, &
142           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,n,xtn,ytn,1, &
143           memtime(1),memtime(2),interp_time,lsp,convp,cc)
144    endif
145
146!  If total precipitation is less than 0.01 mm/h - no scavenging occurs
147    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
148
149! get the level were the actual particle is in
150    do il=2,nz
151      if (height(il).gt.ztra1(jpart)) then
152        hz=il-1
153!        goto 26
154        exit
155      endif
156    end do
157!26  continue
158
159
160    if (ngrid.eq.0) then
161      clouds_v=clouds(ix,jy,hz,n)
162      clouds_h=cloudsh(ix,jy,n)
163    else
164      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
165      clouds_h=cloudshn(ix,jy,n,ngrid)
166    endif
167
168! if there is no precipitation or the particle is above the clouds no
169! scavenging is done
170
171    if (clouds_v.le.1) goto 20
172
173! 1) Parameterization of the the area fraction of the grid cell where the
174!    precipitation occurs: the absolute limit is the total cloud cover, but
175!    for low precipitation rates, an even smaller fraction of the grid cell
176!    is used. Large scale precipitation occurs over larger areas than
177!    convective precipitation.
178!**************************************************************************
179
180    if (lsp.gt.20.) then
181      i=5
182    else if (lsp.gt.8.) then
183      i=4
184    else if (lsp.gt.3.) then
185      i=3
186    else if (lsp.gt.1.) then
187      i=2
188    else
189      i=1
190    endif
191
192    if (convp.gt.20.) then
193      j=5
194    else if (convp.gt.8.) then
195      j=4
196    else if (convp.gt.3.) then
197      j=3
198    else if (convp.gt.1.) then
199      j=2
200    else
201      j=1
202    endif
203
204
205!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
206! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
207! for now they are treated the same
208    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
209    grfraction(2)=max(0.05,cc*(lfr(i)))
210    grfraction(3)=max(0.05,cc*(cfr(j)))
211
212
213! 2) Computation of precipitation rate in sub-grid cell
214!******************************************************
215    prec(1)=(lsp+convp)/grfraction(1)
216    prec(2)=(lsp)/grfraction(2)
217    prec(3)=(convp)/grfraction(3)
218
219
220! 3) Computation of scavenging coefficients for all species
221!    Computation of wet deposition
222!**********************************************************
223
224
225      if (ngrid.gt.0) then
226        act_temp=ttn(ix,jy,hz,n,ngrid)
227      else
228        act_temp=tt(ix,jy,hz,n)
229      endif
230
231!***********************
232! BELOW CLOUD SCAVENGING
233!*********************** 
234      if (clouds_v.ge.4) then !below cloud
235
236! For gas: if positive below-cloud parameters (A or B), and dquer<=0
237!******************************************************************
238        if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
239          !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
240          blc_count(ks)=blc_count(ks)+1
241          wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
242
243! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
244!*********************************************************************************
245        else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
246          blc_count(ks)=blc_count(ks)+1
247
248!NIK 17.02.2015
249! For the calculation here particle size needs to be in meter and not um as dquer is
250! changed in readreleases
251! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
252          dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
253
254! Rain:
255          if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.)  then
256
257! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
258! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
259            wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
260                 & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
261                 &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
262
263! Snow:
264          elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.)  then
265! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
266! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
267            wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
268                 &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
269                 &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
270
271          endif
272         
273!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
274
275        endif ! gas or particle
276!      endif ! positive below-cloud scavenging parameters given in Species file
277      endif !end BELOW
278
279!********************
280! IN CLOUD SCAVENGING
281!********************
282      if (clouds_v.lt.4) then ! In-cloud
283! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
284! given in species file, or if gas and positive Henry's constant
285        if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
286          inc_count(ks)=inc_count(ks)+1
287!          write(*,*) 'Incloud: ',inc_count
288! if negative coefficients (turned off) set to zero for use in equation
289          if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
290          if (in_aero(ks).lt.0.) in_aero(ks)=0.
291
292!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
293! nested fields
294          if (ngrid.gt.0.and.readclouds_this_nest) then
295            cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
296          else if (ngrid.eq.0.and.readclouds) then
297            cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
298          else                                  !parameterize cloudwater m2/m3
299!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
300            cl=1.6E-6*prec(1)**0.36
301          endif
302
303!ZHG: Calculate the partition between liquid and water phase water.
304          if (act_temp .le. 253.) then
305            liq_frac=0
306          else if (act_temp .ge. 273.) then
307            liq_frac=1
308          else
309            liq_frac =((act_temp-273.)/(273.-253.))**2.
310          end if
311! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
312          frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
313
314!ZHG Use the activated fraction and the liqid water to calculate the washout
315
316! AEROSOL
317!********
318          if (dquer(ks).gt.0.) then
319            S_i= frac_act/cl
320
321! GAS
322!****
323          else
324
325            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
326!REPLACE to switch old/ new scheme
327          ! S_i=frac_act/cle
328            S_i=1/cle
329          endif ! gas or particle
330
331! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
332!OLD
333          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
334            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
335          else
336            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
337          endif
338        endif ! positive in-cloud scavenging parameters given in Species file
339      endif !incloud
340
341
34220  continue
343
344end subroutine get_wetscav
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG