source: flexpart.git/src/wetdepo.f90 @ d9f0585

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

Bugfix: using namelist species files would sometimes set incorrect below-cloud scavenging parameters when using more than 1 species. Update: per-species scavenging statistics are written out at end of simulation.

  • Property mode set to 100644
File size: 17.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 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(selected_int_kind(16)), dimension(nspec) :: 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! Cycle loop if wet deposition for the species is off
261      if (WETDEPSPEC(ks).eqv..false.) cycle
262
263      if (ngrid.gt.0) then
264        act_temp=ttn(ix,jy,hz,n,ngrid)
265      else
266        act_temp=tt(ix,jy,hz,n)
267      endif
268
269
270!***********************
271! BELOW CLOUD SCAVENGING
272!*********************** 
273      if (clouds_v.ge.4) then !below cloud
274
275! For gas: if positive below-cloud parameters (A or B), and dquer<=0
276!******************************************************************
277        if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
278          !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
279          blc_count(ks)=blc_count(ks)+1
280          wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
281
282! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
283!*********************************************************************************
284        else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
285          blc_count(ks)=blc_count(ks)+1
286
287!NIK 17.02.2015
288! For the calculation here particle size needs to be in meter and not um as dquer is
289! changed in readreleases
290! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
291          dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
292
293! Rain:
294          if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.)  then
295
296! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
297! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
298            wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
299                 & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
300                 &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
301
302! Snow:
303          elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.)  then
304! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
305! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
306            wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
307                 &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
308                 &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
309
310          endif
311         
312!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
313
314        endif ! gas or particle
315!      endif ! positive below-cloud scavenging parameters given in Species file
316      endif !end BELOW
317
318!********************
319! IN CLOUD SCAVENGING
320!********************
321      if (clouds_v.lt.4) then ! In-cloud
322! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
323! given in species file, or if gas and positive Henry's constant
324        if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
325          inc_count(ks)=inc_count(ks)+1
326! if negative coefficients (turned off) set to zero for use in equation
327          if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
328          if (in_aero(ks).lt.0.) in_aero(ks)=0.
329
330!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
331! nested fields
332          if (ngrid.gt.0.and.readclouds_this_nest) then
333            cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
334          else if (ngrid.eq.0.and.readclouds) then
335            cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
336          else                                  !parameterize cloudwater m2/m3
337!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
338            cl=1.6E-6*prec(1)**0.36
339          endif
340
341!ZHG: Calculate the partition between liquid and water phase water.
342          if (act_temp .le. 253.) then
343            liq_frac=0
344          else if (act_temp .ge. 273.) then
345            liq_frac=1
346          else
347            liq_frac =((act_temp-273.)/(273.-253.))**2.
348          end if
349! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
350          frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
351
352!ZHG Use the activated fraction and the liqid water to calculate the washout
353
354! AEROSOL
355!********
356          if (dquer(ks).gt.0.) then
357            S_i= frac_act/cl
358
359! GAS
360!****
361          else
362
363            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
364!REPLACE to switch old/ new scheme
365          ! S_i=frac_act/cle
366            S_i=1/cle
367          endif ! gas or particle
368
369! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
370!OLD
371          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
372            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
373          else
374            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
375          endif
376
377
378        endif ! positive in-cloud scavenging parameters given in Species file
379      endif !incloud
380!END ZHG TEST
381
382!**************************************************
383! CALCULATE DEPOSITION
384!**************************************************
385!     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
386!     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
387
388      if (wetscav.gt.0.) then
389        wetdeposit(ks)=xmass1(jpart,ks)* &
390             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
391!write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &
392!             (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v
393
394
395      else ! if no scavenging
396        wetdeposit(ks)=0.
397      endif
398
399      restmass = xmass1(jpart,ks)-wetdeposit(ks)
400      if (ioutputforeachrelease.eq.1) then
401        kp=npoint(jpart)
402      else
403        kp=1
404      endif
405      if (restmass .gt. smallnum) then
406        xmass1(jpart,ks)=restmass
407!   depostatistic
408!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
409!   depostatistic
410      else
411        xmass1(jpart,ks)=0.
412      endif
413!   Correct deposited mass to the last time step when radioactive decay of
414!   gridded deposited mass was calculated
415      if (decay(ks).gt.0.) then
416        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
417      endif
418
419
420    end do !all species
421
422! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
423! Add the wet deposition to accumulated amount on output grid and nested output grid
424!*****************************************************************************
425
426    if (ldirect.eq.1) then
427      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
428           real(ytra1(jpart)),nage,kp)
429      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
430           wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
431    endif
432
43320  continue
434  end do ! all particles
435
436! count the total number of below-cloud and in-cloud occurences:
437  tot_blc_count(1:nspec)=tot_blc_count(1:nspec)+blc_count(1:nspec)
438  tot_inc_count(1:nspec)=tot_inc_count(1:nspec)+inc_count(1:nspec)
439
440end subroutine wetdepo
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG