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

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

Completed handling of nested wind fields with cloud water (for wet deposition).

  • Property mode set to 100644
File size: 16.5 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
146    if (ngrid.gt.0) then
147      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
148      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
149      ix=int(xtn)
150      jy=int(ytn)
151      if (readclouds_nest(ngrid)) then
152        readclouds_this_nest=.true.
153      else
154        readclouds_this_nest=.false.
155      end if
156    else
157      ix=int(xtra1(jpart))
158      jy=int(ytra1(jpart))
159    endif
160
161
162! Interpolate large scale precipitation, convective precipitation and
163! total cloud cover
164! Note that interpolated time refers to itime-0.5*ltsample [PS]
165!********************************************************************
166    interp_time=nint(itime-0.5*ltsample)
167
168    if (ngrid.eq.0) then
169      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
170           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
171           memtime(1),memtime(2),interp_time,lsp,convp,cc)
172    else
173      call interpol_rain_nests(lsprecn,convprecn,tccn, &
174           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
175           memtime(1),memtime(2),interp_time,lsp,convp,cc)
176    endif
177
178!  If total precipitation is less than 0.01 mm/h - no scavenging occurs
179    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
180
181! get the level were the actual particle is in
182    do il=2,nz
183      if (height(il).gt.ztra1(jpart)) then
184        hz=il-1
185!        goto 26
186        exit
187      endif
188    end do
189!26  continue
190
191    n=memind(2)
192    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
193         n=memind(1)
194
195    if (ngrid.eq.0) then
196      clouds_v=clouds(ix,jy,hz,n)
197      clouds_h=cloudsh(ix,jy,n)
198    else
199      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
200      clouds_h=cloudshn(ix,jy,n,ngrid)
201    endif
202
203! if there is no precipitation or the particle is above the clouds no
204! scavenging is done
205
206    if (clouds_v.le.1) goto 20
207
208! 1) Parameterization of the the area fraction of the grid cell where the
209!    precipitation occurs: the absolute limit is the total cloud cover, but
210!    for low precipitation rates, an even smaller fraction of the grid cell
211!    is used. Large scale precipitation occurs over larger areas than
212!    convective precipitation.
213!**************************************************************************
214
215    if (lsp.gt.20.) then
216      i=5
217    else if (lsp.gt.8.) then
218      i=4
219    else if (lsp.gt.3.) then
220      i=3
221    else if (lsp.gt.1.) then
222      i=2
223    else
224      i=1
225    endif
226
227    if (convp.gt.20.) then
228      j=5
229    else if (convp.gt.8.) then
230      j=4
231    else if (convp.gt.3.) then
232      j=3
233    else if (convp.gt.1.) then
234      j=2
235    else
236      j=1
237    endif
238
239
240!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
241! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
242! for now they are treated the same
243    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
244    grfraction(2)=max(0.05,cc*(lfr(i)))
245    grfraction(3)=max(0.05,cc*(cfr(j)))
246
247
248! 2) Computation of precipitation rate in sub-grid cell
249!******************************************************
250    prec(1)=(lsp+convp)/grfraction(1)
251    prec(2)=(lsp)/grfraction(2)
252    prec(3)=(convp)/grfraction(3)
253
254
255! 3) Computation of scavenging coefficients for all species
256!    Computation of wet deposition
257!**********************************************************
258
259    do ks=1,nspec      ! loop over species
260      wetdeposit(ks)=0.
261      wetscav=0.   
262
263!ZHG test if it nested?
264      if (ngrid.gt.0) then
265        act_temp=ttn(ix,jy,hz,n,ngrid)
266      else
267        act_temp=tt(ix,jy,hz,n)
268      endif
269
270
271!***********************
272! BELOW CLOUD SCAVENGING
273!*********************** 
274      if (clouds_v.ge.4) then !below cloud
275
276        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then !if positive below-cloud parameters given in SPECIES file (either A or B)
277          blc_count=blc_count+1
278
279!GAS
280          if (dquer(ks) .le. 0.) then  !is gas
281            wetscav=weta(ks)*prec(1)**wetb(ks)
282
283!AEROSOL
284          else !is particle
285!NIK 17.02.2015
286! For the calculation here particle size needs to be in meter and not um as dquer is changed to 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            if (act_temp .ge. 273. .and. weta(ks).gt.0.)  then !Rain
290! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
291! the below-cloud scavenging (rain efficienty)
292! parameter A (=weta) from SPECIES file
293              wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* &
294                   (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5))
295
296            elseif (act_temp .lt. 273. .and. wetb(ks).gt.0.)  then ! Snow
297! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
298! the below-cloud scavenging (Snow efficiency)
299! parameter B (=wetb) from SPECIES file
300              wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* &
301                   (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
302
303            endif
304
305!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
306
307          endif !gas or particle
308        endif ! positive below-cloud scavenging parameters given in Species file
309      endif !end BELOW
310
311
312!********************
313! IN CLOUD SCAVENGING
314!********************
315      if (clouds_v.lt.4) then ! In-cloud
316! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file
317        if (weta_in(ks).gt.0. .or. wetb_in(ks).gt.0.) then
318          inc_count=inc_count+1
319! if negative coefficients (turned off) set to zero for use in equation
320          if (weta_in(ks).lt.0.) weta_in(ks)=0.
321          if (wetb_in(ks).lt.0.) wetb_in(ks)=0.
322
323!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
324! nested fields
325          if (ngrid.gt.0.and.readclouds_this_nest) then
326            cl=clw4n(ix,jy,n,ngrid)*(grfraction(1)/cc)
327          else if (ngrid.eq.0.and.readclouds) then
328            cl=clw4(ix,jy,n)*(grfraction(1)/cc)
329          else                                  !parameterize cloudwater m2/m3
330!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
331            cl=1.6E-6*prec(1)**0.36
332          endif
333
334!ZHG: Calculate the partition between liquid and water phase water.
335          if (act_temp .le. 253.) then
336            liq_frac=0
337          else if (act_temp .ge. 273.) then
338            liq_frac=1
339          else
340            liq_frac =((act_temp-273.)/(273.-253.))**2.
341          end if
342! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
343          frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
344
345!ZHG Use the activated fraction and the liqid water to calculate the washout
346
347! AEROSOL
348!**************************************************
349          if (dquer(ks).gt. 0.) then ! is particle
350
351            S_i= frac_act/cl
352
353!*********************
354! GAS
355          else ! is gas
356
357            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
358!REPLACE to switch old/ new scheme
359! S_i=frac_act/cle
360            S_i=1/cle
361          endif ! gas or particle
362
363! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
364!OLD
365          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
366            wetscav=S_i*(prec(1)/3.6E6)
367          else
368            wetscav=S_i*(prec(1)/3.6E6)/clouds_h
369          endif
370
371
372        endif ! positive in-cloud scavenging parameters given in Species file
373      endif !incloud
374!END ZHG TEST
375
376!**************************************************
377! CALCULATE DEPOSITION
378!**************************************************
379!     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
380!     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
381
382      if (wetscav.gt.0.) then
383        wetdeposit(ks)=xmass1(jpart,ks)* &
384             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
385!write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &
386!             (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v
387
388
389      else ! if no scavenging
390        wetdeposit(ks)=0.
391      endif
392
393      restmass = xmass1(jpart,ks)-wetdeposit(ks)
394      if (ioutputforeachrelease.eq.1) then
395        kp=npoint(jpart)
396      else
397        kp=1
398      endif
399      if (restmass .gt. smallnum) then
400        xmass1(jpart,ks)=restmass
401!   depostatistic
402!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
403!   depostatistic
404      else
405        xmass1(jpart,ks)=0.
406      endif
407!   Correct deposited mass to the last time step when radioactive decay of
408!   gridded deposited mass was calculated
409      if (decay(ks).gt.0.) then
410        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
411      endif
412
413
414    end do !all species
415
416! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
417! Add the wet deposition to accumulated amount on output grid and nested output grid
418!*****************************************************************************
419
420    if (ldirect.eq.1) then
421      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
422           real(ytra1(jpart)),nage,kp)
423      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
424           wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
425    endif
426
42720  continue
428  end do ! all particles
429
430! count the total number of below-cloud and in-cloud occurences:
431  tot_blc_count=tot_blc_count+blc_count
432  tot_inc_count=tot_inc_count+inc_count
433
434end subroutine wetdepo
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG