Changes from branches/FLEXPART_9.1.3/src/wetdepo.f90 at r13 to trunk/src/wetdepo.f90 at r4
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/wetdepo.f90
r13 r4 39 39 ! Code may not be correct for decay of deposition! * 40 40 ! * 41 ! Modification by Sabine Eckhart to introduce a new in-/below-cloud42 ! scheme, not dated43 ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification44 41 !***************************************************************************** 45 42 ! * … … 74 71 75 72 integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy 76 integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme 77 integer :: kz !PS scheme 78 integer :: ks, kp, n1,n2, icbot,ictop, indcloud 79 integer :: scheme_number ! NIK==1, PS ==2 73 integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v 74 integer :: ks, kp 80 75 real :: S_i, act_temp, cl, cle ! in cloud scavenging 81 76 real :: clouds_h ! cloud height for the specific grid point 82 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav , precsub,f77 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav 83 78 real :: wetdeposit(maxspec),restmass 84 79 real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled … … 151 146 interp_time=nint(itime-0.5*ltsample) 152 147 153 ! PS nest case still needs to be implemented!! 154 ! if (ngrid.eq.0) then 155 ! call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, & 156 ! 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, & 157 ! memtime(1),memtime(2),interp_time,lsp,convp,cc) 158 call interpol_rain(lsprec,convprec,tcc, & 159 icloudbot,icloudthck,nxmax,nymax,1,nx,ny, & 160 memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), & 161 memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv) 162 ! else 163 ! call interpol_rain_nests(lsprecn,convprecn,tccn, & 164 ! nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 165 ! memtime(1),memtime(2),interp_time,lsp,convp,cc) 166 ! endif 167 168 169 ! if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 170 !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip 171 prec = lsp+convp 172 precsub = 0.01 173 if (prec .lt. precsub) then 174 goto 20 175 else 176 f = (prec-precsub)/prec 177 lsp = f*lsp 178 convp = f*convp 179 endif 180 148 if (ngrid.eq.0) then 149 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, & 150 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, & 151 memtime(1),memtime(2),interp_time,lsp,convp,cc) 152 else 153 call interpol_rain_nests(lsprecn,convprecn,tccn, & 154 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 155 memtime(1),memtime(2),interp_time,lsp,convp,cc) 156 endif 157 158 if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 181 159 182 160 ! get the level were the actual particle is in 183 161 do il=2,nz 184 162 if (height(il).gt.ztra1(jpart)) then 185 !hz=il-1 186 kz=il-1 163 hz=il-1 187 164 goto 26 188 165 endif … … 197 174 ! scavenging is done 198 175 199 !old scheme 200 ! if (ngrid.eq.0) then 201 ! clouds_v=clouds(ix,jy,hz,n) 202 ! clouds_h=cloudsh(ix,jy,n) 203 ! else 204 ! clouds_v=cloudsn(ix,jy,hz,n,ngrid) 205 ! clouds_h=cloudsnh(ix,jy,n,ngrid) 206 ! endif 207 ! !write(*,*) 'there is 208 ! ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz 209 ! if (clouds_v.le.1) goto 20 210 ! !write (*,*) 'there is scavenging' 211 212 ! PS: part of 2011/2012 fix 213 214 if (ztra1(jpart) .le. float(ictop)) then 215 if (ztra1(jpart) .gt. float(icbot)) then 216 indcloud = 2 ! in-cloud 217 else 218 indcloud = 1 ! below-cloud 219 endif 220 elseif (ictop .eq. icmv) then 221 indcloud = 0 ! no cloud found, use old scheme 222 else 223 goto 20 ! above cloud 224 endif 225 176 if (ngrid.eq.0) then 177 clouds_v=clouds(ix,jy,hz,n) 178 clouds_h=cloudsh(ix,jy,n) 179 else 180 clouds_v=cloudsn(ix,jy,hz,n,ngrid) 181 clouds_h=cloudsnh(ix,jy,n,ngrid) 182 endif 183 !write(*,*) 'there is 184 ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz 185 if (clouds_v.le.1) goto 20 186 !write (*,*) 'there is scavenging' 226 187 227 188 ! 1) Parameterization of the the area fraction of the grid cell where the … … 267 228 !********************************************************** 268 229 269 do ks=1,nspec ! loop over species230 do ks=1,nspec ! loop over species 270 231 wetdeposit(ks)=0. 271 272 273 !conflicting changes to the same routine: 1=NIK 2 =PS274 scheme_number=2275 if (scheme_number.eq.1) then !NIK276 277 232 if (weta(ks).gt.0.) then 278 233 if (clouds_v.ge.4) then … … 281 236 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coeff. 282 237 ! write(*,*) 'bel. wetscav: ',wetscav 283 284 238 else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging 285 239 ! IN CLOUD SCAVENGING … … 291 245 act_temp=tt(ix,jy,hz,n) 292 246 endif 293 294 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 295 ! weta_in=2.0E-07 (default) 296 ! wetb_in=0.36 (default) 297 ! wetc_in=0.9 (default) 298 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling) 299 cl=weta_in(ks)*prec**wetb_in(ks) 247 cl=2E-7*prec**0.36 300 248 if (dquer(ks).gt.0) then ! is particle 301 S_i= wetc_in(ks)/cl249 S_i=0.9/cl 302 250 else ! is gas 303 251 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 304 252 S_i=1/cle 305 253 endif 306 wetscav=S_i*prec/3.6E6/clouds_h /wetd_in(ks)254 wetscav=S_i*prec/3.6E6/clouds_h 307 255 ! write(*,*) 'in. wetscav:' 308 256 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h … … 341 289 wetdeposit(ks)=0. 342 290 endif ! weta(k) 343 344 elseif (scheme_number.eq.2) then ! PS345 346 !PS indcloud=0 ! Use this for FOR TESTING,347 !PS will skip the new in/below cloud method !!!348 349 if (weta(ks).gt.0.) then350 if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING351 !C for aerosols and not highliy soluble substances weta=5E-6352 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coeff.353 !c write(*,*) 'bel. wetscav: ',wetscav354 elseif (indcloud .eq. 2) then ! IN CLOUD SCAVENGING355 if (ngrid.gt.0) then356 act_temp=ttn(ix,jy,kz,n,ngrid)357 else358 act_temp=tt(ix,jy,kz,n)359 endif360 361 ! from NIK362 ! weta_in=2.0E-07 (default)363 ! wetb_in=0.36 (default)364 ! wetc_in=0.9 (default)365 366 367 cl=2E-7*prec**0.36368 if (dquer(ks).gt.0) then ! is particle369 S_i=0.9/cl370 else ! is gas371 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl372 S_i=1/cle373 endif374 wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s375 else ! PS: no cloud diagnosed, old scheme,376 !CPS using with fixed a,b for simplicity, one may wish to change!!377 wetscav = 1.e-4*prec**0.62378 endif379 380 381 wetdeposit(ks)=xmass1(jpart,ks)* &382 ! (1.-exp(-wetscav*abs(ltsample)))*fraction ! wet deposition383 (1.-exp(-wetscav*abs(ltsample)))*grfraction ! fraction = grfraction (PS)384 restmass = xmass1(jpart,ks)-wetdeposit(ks)385 if (ioutputforeachrelease.eq.1) then386 kp=npoint(jpart)387 else388 kp=1389 endif390 if (restmass .gt. smallnum) then391 xmass1(jpart,ks)=restmass392 !cccccccccccccccc depostatistic393 !c wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)394 !cccccccccccccccc depostatistic395 else396 xmass1(jpart,ks)=0.397 endif398 !C Correct deposited mass to the last time step when radioactive decay of399 !C gridded deposited mass was calculated400 if (decay(ks).gt.0.) then401 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))402 endif403 else ! weta(k)<0404 wetdeposit(ks)=0.405 endif406 407 endif !on scheme408 409 410 411 291 end do 412 292 … … 415 295 !***************************************************************************** 416 296 417 ! if (ldirect.eq.1) then 418 ! call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), & 419 ! real(ytra1(jpart)),nage,kp) 420 ! if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), & 421 ! wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), & 422 ! nage,kp) 423 ! endif 424 425 !PS 426 if (ldirect.eq.1) then 427 call wetdepokernel(nclass(jpart),wetdeposit, & 428 sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 429 if (nested_output.eq.1) & 430 call wetdepokernel_nest(nclass(jpart),wetdeposit, & 431 sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 432 endif 433 297 if (ldirect.eq.1) then 298 call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), & 299 real(ytra1(jpart)),nage,kp) 300 if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), & 301 wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), & 302 nage,kp) 303 endif 434 304 435 305 20 continue
Note: See TracChangeset
for help on using the changeset viewer.