source: flexpart.git/src/wetdepo.f90 @ 4c64400

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

Updated checks and warning messages for wet deposition parameters given in SPECIES files

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