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

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

interpolate the rain only in space and not in time

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