source: flexpart.git/src/get_wetscav.f90 @ e9e0f06

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since e9e0f06 was e9e0f06, checked in by Sabine <sabine.eckhardt@…>, 5 years ago

Removed kao and ass_spec and obsolete lines in get_wetscav.f90

  • Property mode set to 100644
File size: 13.8 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        exit
154      endif
155    end do
156
157    if (ngrid.eq.0) then
158      clouds_v=clouds(ix,jy,hz,n)
159      clouds_h=cloudsh(ix,jy,n)
160    else
161      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
162      clouds_h=cloudshn(ix,jy,n,ngrid)
163    endif
164
165! if there is no precipitation or the particle is above the clouds no
166! scavenging is done
167
168    if (clouds_v.le.1) goto 20
169
170! 1) Parameterization of the the area fraction of the grid cell where the
171!    precipitation occurs: the absolute limit is the total cloud cover, but
172!    for low precipitation rates, an even smaller fraction of the grid cell
173!    is used. Large scale precipitation occurs over larger areas than
174!    convective precipitation.
175!**************************************************************************
176
177    if (lsp.gt.20.) then
178      i=5
179    else if (lsp.gt.8.) then
180      i=4
181    else if (lsp.gt.3.) then
182      i=3
183    else if (lsp.gt.1.) then
184      i=2
185    else
186      i=1
187    endif
188
189    if (convp.gt.20.) then
190      j=5
191    else if (convp.gt.8.) then
192      j=4
193    else if (convp.gt.3.) then
194      j=3
195    else if (convp.gt.1.) then
196      j=2
197    else
198      j=1
199    endif
200
201
202!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp - 2 and 3 not used removed by SE
203! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
204! for now they are treated the same
205    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
206
207! 2) Computation of precipitation rate in sub-grid cell
208!******************************************************
209    prec(1)=(lsp+convp)/grfraction(1)
210
211! 3) Computation of scavenging coefficients for all species
212!    Computation of wet deposition
213!**********************************************************
214
215      if (ngrid.gt.0) then
216        act_temp=ttn(ix,jy,hz,n,ngrid)
217      else
218        act_temp=tt(ix,jy,hz,n)
219      endif
220
221!***********************
222! BELOW CLOUD SCAVENGING
223!*********************** 
224      if (clouds_v.ge.4) then !below cloud
225
226! For gas: if positive below-cloud parameters (A or B), and dquer<=0
227!******************************************************************
228        if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
229          blc_count(ks)=blc_count(ks)+1
230          wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
231
232! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
233!*********************************************************************************
234        else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
235          blc_count(ks)=blc_count(ks)+1
236
237!NIK 17.02.2015
238! For the calculation here particle size needs to be in meter and not um as dquer is
239! changed in readreleases
240! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
241          dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
242
243! Rain:
244          if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.)  then
245
246! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
247! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
248            wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
249                 & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
250                 &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
251
252! Snow:
253          elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.)  then
254! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
255! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
256            wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
257                 &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
258                 &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
259
260          endif
261         
262        endif ! gas or particle
263!      endif ! positive below-cloud scavenging parameters given in Species file
264      endif !end BELOW
265
266!********************
267! IN CLOUD SCAVENGING
268!********************
269      if (clouds_v.lt.4) then ! In-cloud
270! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
271! given in species file, or if gas and positive Henry's constant
272        if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
273          inc_count(ks)=inc_count(ks)+1
274! if negative coefficients (turned off) set to zero for use in equation
275          if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
276          if (in_aero(ks).lt.0.) in_aero(ks)=0.
277
278!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
279! nested fields
280          if (ngrid.gt.0.and.readclouds_this_nest) then
281            cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
282          else if (ngrid.eq.0.and.readclouds) then
283            cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
284          else                                  !parameterize cloudwater m2/m3
285!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
286! sec test
287!           cl=1E6*1E-7*prec(1)**0.3 !Sec GFS new
288            cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new, is also suitable for GFS
289!           cl=2E-7*prec(1)**0.36 !Andreas
290!           cl=1.6E-6*prec(1)**0.36 !Henrik
291          endif
292
293!ZHG: Calculate the partition between liquid and water phase water.
294          if (act_temp .le. 253.) then
295            liq_frac=0
296          else if (act_temp .ge. 273.) then
297            liq_frac=1
298          else
299            liq_frac =((act_temp-273.)/(273.-253.))**2.
300          end if
301! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
302          frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
303
304!ZHG Use the activated fraction and the liqid water to calculate the washout
305
306! AEROSOL
307!********
308          if (dquer(ks).gt.0.) then
309            S_i= frac_act/cl
310! GAS
311!****
312          else
313            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
314            S_i=1/cle
315          endif ! gas or particle
316
317! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
318!SEC wetscav fix, the cloud height is no longer needed, it gives wrong results
319            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
320        endif ! positive in-cloud scavenging parameters given in Species file
321      endif !incloud
322
323
32420  continue
325
326end subroutine get_wetscav
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG