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

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

Minor cosmetic edits

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