1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
|
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine wetdepo(itime,ltsample,loutnext)
! i i i
!*****************************************************************************
! *
! Calculation of wet deposition using the concept of scavenging coefficients.*
! For lack of detailed information, washout and rainout are jointly treated. *
! It is assumed that precipitation does not occur uniformly within the whole *
! grid cell, but that only a fraction of the grid cell experiences rainfall. *
! This fraction is parameterized from total cloud cover and rates of large *
! scale and convective precipitation. *
! *
! Author: A. Stohl *
! *
! 1 December 1996 *
! *
! Correction by Petra Seibert, Sept 2002: *
! use centred precipitation data for integration *
! Code may not be correct for decay of deposition! *
! *
!*****************************************************************************
! *
! Variables: *
! cc [0-1] total cloud cover *
! convp [mm/h] convective precipitation rate *
! grfraction [0-1] fraction of grid, for which precipitation occurs *
! ix,jy indices of output grid cell for each particle *
! itime [s] actual simulation time [s] *
! jpart particle index *
! ldeltat [s] interval since radioactive decay was computed *
! lfr, cfr area fraction covered by precipitation for large scale *
! and convective precipitation (dependent on prec. rate) *
! loutnext [s] time for which gridded deposition is next output *
! loutstep [s] interval at which gridded deposition is output *
! lsp [mm/h] large scale precipitation rate *
! ltsample [s] interval over which mass is deposited *
! prec [mm/h] precipitation rate in subgrid, where precipitation occurs*
! wetdeposit mass that is wet deposited *
! wetgrid accumulated deposited mass on output grid *
! wetscav scavenging coefficient *
! *
! Constants: *
! *
!*****************************************************************************
use point_mod
use par_mod
use com_mod
implicit none
integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
integer :: ks, kp
real :: S_i, act_temp, cl, cle ! in cloud scavenging
real :: clouds_h ! cloud height for the specific grid point
real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
real :: wetdeposit(maxspec),restmass
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
save lfr,cfr
real :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
real :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
! Compute interval since radioactive decay of deposited mass was computed
!************************************************************************
if (itime.le.loutnext) then
ldeltat=itime-(loutnext-loutstep)
else ! first half of next interval
ldeltat=itime-loutnext
endif
! Loop over all particles
!************************
do jpart=1,numpart
if (itra1(jpart).eq.-999999999) goto 20
if(ldirect.eq.1)then
if (itra1(jpart).gt.itime) goto 20
else
if (itra1(jpart).lt.itime) goto 20
endif
! Determine age class of the particle
itage=abs(itra1(jpart)-itramem(jpart))
do nage=1,nageclass
if (itage.lt.lage(nage)) goto 33
end do
33 continue
! Determine which nesting level to be used
!*****************************************
ngrid=0
do j=numbnests,1,-1
if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
(ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
ngrid=j
goto 23
endif
end do
23 continue
! Determine nested grid coordinates
!**********************************
if (ngrid.gt.0) then
xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
ix=int(xtn)
jy=int(ytn)
else
ix=int(xtra1(jpart))
jy=int(ytra1(jpart))
endif
! Interpolate large scale precipitation, convective precipitation and
! total cloud cover
! Note that interpolated time refers to itime-0.5*ltsample [PS]
!********************************************************************
interp_time=nint(itime-0.5*ltsample)
if (ngrid.eq.0) then
call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
memtime(1),memtime(2),interp_time,lsp,convp,cc)
else
call interpol_rain_nests(lsprecn,convprecn,tccn, &
nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
memtime(1),memtime(2),interp_time,lsp,convp,cc)
endif
if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
! get the level were the actual particle is in
do il=2,nz
if (height(il).gt.ztra1(jpart)) then
hz=il-1
goto 26
endif
end do
26 continue
n=memind(2)
if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
n=memind(1)
! if there is no precipitation or the particle is above the clouds no
! scavenging is done
if (ngrid.eq.0) then
clouds_v=clouds(ix,jy,hz,n)
clouds_h=cloudsh(ix,jy,n)
else
clouds_v=cloudsn(ix,jy,hz,n,ngrid)
clouds_h=cloudsnh(ix,jy,n,ngrid)
endif
!write(*,*) 'there is
! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
if (clouds_v.le.1) goto 20
!write (*,*) 'there is scavenging'
! 1) Parameterization of the the area fraction of the grid cell where the
! precipitation occurs: the absolute limit is the total cloud cover, but
! for low precipitation rates, an even smaller fraction of the grid cell
! is used. Large scale precipitation occurs over larger areas than
! convective precipitation.
!**************************************************************************
if (lsp.gt.20.) then
i=5
else if (lsp.gt.8.) then
i=4
else if (lsp.gt.3.) then
i=3
else if (lsp.gt.1.) then
i=2
else
i=1
endif
if (convp.gt.20.) then
j=5
else if (convp.gt.8.) then
j=4
else if (convp.gt.3.) then
j=3
else if (convp.gt.1.) then
j=2
else
j=1
endif
grfraction=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
! 2) Computation of precipitation rate in sub-grid cell
!******************************************************
prec=(lsp+convp)/grfraction
! 3) Computation of scavenging coefficients for all species
! Computation of wet deposition
!**********************************************************
do ks=1,nspec ! loop over species
wetdeposit(ks)=0.
if (weta(ks).gt.0.) then
if (clouds_v.ge.4) then
! BELOW CLOUD SCAVENGING
! for aerosols and not highliy soluble substances weta=5E-6
wetscav=weta(ks)*prec**wetb(ks) ! scavenging coeff.
! write(*,*) 'bel. wetscav: ',wetscav
else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
! IN CLOUD SCAVENGING
! BUGFIX tt for nested fields should be ttn
! sec may 2008
if (ngrid.gt.0) then
act_temp=ttn(ix,jy,hz,n,ngrid)
else
act_temp=tt(ix,jy,hz,n)
endif
cl=2E-7*prec**0.36
if (dquer(ks).gt.0) then ! is particle
S_i=0.9/cl
else ! is gas
cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
S_i=1/cle
endif
wetscav=S_i*prec/3.6E6/clouds_h
! write(*,*) 'in. wetscav:'
! + ,wetscav,cle,cl,act_temp,prec,clouds_h
endif
! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
wetdeposit(ks)=xmass1(jpart,ks)* &
(1.-exp(-wetscav*abs(ltsample)))*grfraction ! wet deposition
! new particle mass:
! if (wetdeposit(ks).gt.0) then
! write(*,*) 'wetdepo: ',wetdeposit(ks),ks
! endif
restmass = xmass1(jpart,ks)-wetdeposit(ks)
if (ioutputforeachrelease.eq.1) then
kp=npoint(jpart)
else
kp=1
endif
if (restmass .gt. smallnum) then
xmass1(jpart,ks)=restmass
!ccccccccccccccc depostatistic
! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
!ccccccccccccccc depostatistic
else
xmass1(jpart,ks)=0.
endif
! Correct deposited mass to the last time step when radioactive decay of
! gridded deposited mass was calculated
if (decay(ks).gt.0.) then
wetdeposit(ks)=wetdeposit(ks) &
*exp(abs(ldeltat)*decay(ks))
endif
else ! weta(k)
wetdeposit(ks)=0.
endif ! weta(k)
end do
! Sabine Eckhard, June 2008 create deposition runs only for forward runs
! Add the wet deposition to accumulated amount on output grid and nested output grid
!*****************************************************************************
if (ldirect.eq.1) then
call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
real(ytra1(jpart)),nage,kp)
if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
nage,kp)
endif
20 continue
end do
end subroutine wetdepo
|