Changeset c9cf570 in flexpart.git
- Timestamp:
- Nov 18, 2016, 11:29:23 AM (7 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 7062fab
- Parents:
- deaac29
- Location:
- src
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
src/timemanager.f90
rdeaac29 rc9cf570 105 105 ! integer :: ksp 106 106 integer :: loutnext,loutstart,loutend 107 integer :: ix,jy,ldeltat,itage,nage 107 integer :: ix,jy,ldeltat,itage,nage,idummy 108 108 integer :: i_nan=0,ii_nan,total_nan_intl=0 !added by mc to check instability in CBL scheme 109 real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact 109 real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact,wetscav(maxspec) 110 110 ! real :: uap(maxpart),ucp(maxpart),uzp(maxpart) 111 111 ! real :: us(maxpart),vs(maxpart),ws(maxpart) … … 114 114 real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc 115 115 real :: xold,yold,zold,xmassfract 116 real :: grfraction(3) 116 117 real, parameter :: e_inv = 1.0/exp(1.0) 117 118 … … 550 551 551 552 if (DRYBKDEP) then 552 do ks=1,nspec553 do ks=1,nspec 553 554 if ((xscav_frac1(j,ks).lt.0)) then 554 555 call advance_rec(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec) 555 556 if (decay(ks).gt.0.) then ! radioactive decay557 decfact=exp(-real(abs(lsynctime))*decay(ks))558 else559 decfact=1.560 endif561 556 if (DRYDEPSPEC(ks)) then ! dry deposition 562 557 xscav_frac1(j,ks)=prob_rec(ks) … … 565 560 xscav_frac1(j,ks)=0. 566 561 endif 567 ! write (*,*) 'xscav: ',j,ks,xscav_frac1(j,ks)568 562 endif 569 enddo563 enddo 570 564 endif 571 565 572 ! if (WETBKDEP) then 573 ! firstdepocalc=.false. 574 ! do ks=1,nspec 575 ! if ((xscav_frac1(j,ks).lt.0) & 576 ! .and.firstdepocalc.eqv..false.) then 577 ! ! Backward wetdeposition and first timestep after release 578 ! call wetdepo(itime,lsynctime,loutnext,.true.) 579 ! firstdepocalc=.true. 580 ! endif 581 ! enddo 582 ! endif 566 if (WETBKDEP) then 567 do ks=1,nspec 568 if ((xscav_frac1(j,ks).lt.0)) then 569 call get_wetscav(itime,lsynctime,loutnext,j,grfraction,idummy,idummy,wetscav) 570 if (wetscav(ks).gt.0) then 571 xscav_frac1(j,ks)=wetscav(ks)*(zpoint2(j)-zpoint1(j)) 572 else 573 xmass1(j,ks)=0. 574 xscav_frac1(j,ks)=0. 575 endif 576 endif 577 enddo 578 endif 583 579 584 580 ! Integrate Lagevin equation for lsynctime seconds -
src/wetdepo.f90
rdeaac29 rc9cf570 42 42 ! * 43 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 44 ! ix,jy indices of output grid cell for each particle * 48 45 ! itime [s] actual simulation time [s] * 49 46 ! jpart particle index * 50 47 ! 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 48 ! loutnext [s] time for which gridded deposition is next output * 54 49 ! loutstep [s] interval at which gridded deposition is output * 55 ! lsp [mm/h] large scale precipitation rate *56 50 ! ltsample [s] interval over which mass is deposited * 57 ! prec [mm/h] precipitation rate in subgrid, where precipitation occurs*58 51 ! wetdeposit mass that is wet deposited * 59 52 ! wetgrid accumulated deposited mass on output grid * … … 71 64 72 65 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 66 integer :: itage,nage,hz,il,interp_time, n 75 67 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 68 integer :: blc_count, inc_count 69 real :: grfraction(3),wetscav 80 70 real :: wetdeposit(maxspec),restmass 81 71 real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled 82 !save lfr,cfr83 72 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 snow88 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_m91 92 integer :: blc_count, inc_count, firsttimerem93 real :: Si_dummy, wetscav_dummy94 logical :: readclouds_this_nest95 96 firsttimerem=097 73 ! Compute interval since radioactive decay of deposited mass was computed 98 74 !************************************************************************ … … 119 95 endif 120 96 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 126 33 continue 97 ! Determine age class of the particle - nage is used for the kernel 98 !****************************************************************** 99 itage=abs(itra1(jpart)-itramem(jpart)) 100 do nage=1,nageclass 101 if (itage.lt.lage(nage)) goto 33 102 end do 103 33 continue 127 104 128 129 ! Determine which nesting level to be used130 !*****************************************131 132 ngrid=0133 do j=numbnests,1,-1134 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))) then136 ngrid=j137 goto 23138 endif139 end do140 23 continue141 142 143 ! Determine nested grid coordinates144 !**********************************145 readclouds_this_nest=.false.146 147 if (ngrid.gt.0) then148 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 else154 ix=int(xtra1(jpart))155 jy=int(ytra1(jpart))156 endif157 158 159 ! Interpolate large scale precipitation, convective precipitation and160 ! total cloud cover161 ! 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) then166 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 else170 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 endif174 175 ! If total precipitation is less than 0.01 mm/h - no scavenging occurs176 ! sec this is just valid if release is over a point177 if ((lsp.lt.0.01).and.(convp.lt.0.01)) then178 if (WETBKDEP) then179 do ks=1,nspec180 if (xscav_frac1(jpart,ks).lt.0) then ! first timestep no scavenging181 xmass1(jpart,ks)=0.182 xscav_frac1(jpart,ks)=0.183 ! write (*,*) 'paricle removed - no scavenging: ',jpart,ks184 endif185 end do186 endif187 goto 20188 endif189 190 191 ! get the level were the actual particle is in192 do il=2,nz193 if (height(il).gt.ztra1(jpart)) then194 hz=il-1195 ! goto 26196 exit197 endif198 end do199 !26 continue200 201 n=memind(2)202 if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &203 n=memind(1)204 205 if (ngrid.eq.0) then206 clouds_v=clouds(ix,jy,hz,n)207 clouds_h=cloudsh(ix,jy,n)208 else209 clouds_v=cloudsn(ix,jy,hz,n,ngrid)210 clouds_h=cloudshn(ix,jy,n,ngrid)211 endif212 213 ! if there is no precipitation or the particle is above the clouds no214 ! scavenging is done215 216 if (clouds_v.le.1) goto 20217 218 ! 1) Parameterization of the the area fraction of the grid cell where the219 ! precipitation occurs: the absolute limit is the total cloud cover, but220 ! for low precipitation rates, an even smaller fraction of the grid cell221 ! is used. Large scale precipitation occurs over larger areas than222 ! convective precipitation.223 !**************************************************************************224 225 if (lsp.gt.20.) then226 i=5227 else if (lsp.gt.8.) then228 i=4229 else if (lsp.gt.3.) then230 i=3231 else if (lsp.gt.1.) then232 i=2233 else234 i=1235 endif236 237 if (convp.gt.20.) then238 j=5239 else if (convp.gt.8.) then240 j=4241 else if (convp.gt.3.) then242 j=3243 else if (convp.gt.1.) then244 j=2245 else246 j=1247 endif248 249 250 !ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp251 ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms252 ! for now they are treated the same253 grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))254 grfraction(2)=max(0.05,cc*(lfr(i)))255 grfraction(3)=max(0.05,cc*(cfr(j)))256 257 258 ! 2) Computation of precipitation rate in sub-grid cell259 !******************************************************260 prec(1)=(lsp+convp)/grfraction(1)261 prec(2)=(lsp)/grfraction(2)262 prec(3)=(convp)/grfraction(3)263 264 265 ! 3) Computation of scavenging coefficients for all species266 ! Computation of wet deposition267 !**********************************************************268 269 do ks=1,nspec ! loop over species270 wetdeposit(ks)=0.271 wetscav=0.272 273 if (ngrid.gt.0) then274 act_temp=ttn(ix,jy,hz,n,ngrid)275 else276 act_temp=tt(ix,jy,hz,n)277 endif278 279 280 !***********************281 ! BELOW CLOUD SCAVENGING282 !***********************283 if (clouds_v.ge.4) then !below cloud284 285 ! For gas: if positive below-cloud parameters (A or B), and dquer<=0286 !******************************************************************287 if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then288 ! if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then289 blc_count=blc_count+1290 wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)291 292 ! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0293 !*********************************************************************************294 else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then295 blc_count=blc_count+1296 297 !NIK 17.02.2015298 ! For the calculation here particle size needs to be in meter and not um as dquer is299 ! changed in readreleases300 ! For particles larger than 10 um use the largest size defined in the parameterizations (10um)301 dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m302 303 ! Rain:304 if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.) then305 306 ! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,307 ! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file308 wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &309 & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&310 &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))311 312 ! Snow:313 elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.) then314 ! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,315 ! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file316 wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&317 &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&318 &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))319 320 endif321 322 ! write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart323 324 endif ! gas or particle325 ! endif ! positive below-cloud scavenging parameters given in Species file326 endif !end BELOW327 328 !********************329 ! IN CLOUD SCAVENGING330 !********************331 if (clouds_v.lt.4) then ! In-cloud332 ! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are333 ! given in species file, or if gas and positive Henry's constant334 if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then335 inc_count=inc_count+1336 ! if negative coefficients (turned off) set to zero for use in equation337 if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.338 if (in_aero(ks).lt.0.) in_aero(ks)=0.339 340 !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF341 ! nested fields342 if (ngrid.gt.0.and.readclouds_this_nest) then343 cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)344 else if (ngrid.eq.0.and.readclouds) then345 cl=ctwc(ix,jy,n)*(grfraction(1)/cc)346 else !parameterize cloudwater m2/m3347 !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF348 cl=1.6E-6*prec(1)**0.36349 endif350 351 !ZHG: Calculate the partition between liquid and water phase water.352 if (act_temp .le. 253.) then353 liq_frac=0354 else if (act_temp .ge. 273.) then355 liq_frac=1356 else357 liq_frac =((act_temp-273.)/(273.-253.))**2.358 end if359 ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi360 frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)361 362 !ZHG Use the activated fraction and the liqid water to calculate the washout363 364 ! AEROSOL365 !********366 if (dquer(ks).gt.0.) then367 S_i= frac_act/cl368 369 ! GAS370 !****371 else372 373 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl374 !REPLACE to switch old/ new scheme375 ! S_i=frac_act/cle376 S_i=1/cle377 endif ! gas or particle378 379 ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol380 !OLD381 if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then382 wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)383 else384 wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h385 endif386 387 388 endif ! positive in-cloud scavenging parameters given in Species file389 endif !incloud390 !END ZHG TEST391 105 392 106 !************************************************** 393 107 ! CALCULATE DEPOSITION 394 108 !************************************************** 395 ! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!' 396 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v 109 call get_wetscav(itime,ltsample,loutnext,jpart,grfraction,inc_count,blc_count,wetscav) 397 110 398 111 if (wetscav.gt.0.) then 399 112 wetdeposit(ks)=xmass1(jpart,ks)* & 400 113 (1.-exp(-wetscav*abs(ltsample)))*grfraction(1) ! wet deposition 401 !write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &402 ! (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v403 404 405 114 else ! if no scavenging 406 115 wetdeposit(ks)=0. 407 116 endif 408 117 409 118 restmass = xmass1(jpart,ks)-wetdeposit(ks) 410 119 if (ioutputforeachrelease.eq.1) then … … 427 136 endif 428 137 429 if (WETBKDEP) then430 ! the calculation of the scavenged mass shall only be done once after release431 ! xscav_frac1 was initialised with a negative value432 if (xscav_frac1(jpart,ks).lt.0) then433 if (wetdeposit(ks).eq.0) then434 ! terminate particle435 xmass1(jpart,ks)=0.436 xscav_frac1(jpart,ks)=0.437 else438 if (xmass1(jpart,ks).eq.0) then439 firsttimerem=firsttimerem+1440 endif441 !xscav_frac1(jpart,ks)=wetdeposit(ks)/(xmass1(jpart,ks)+wetdeposit(ks))442 xscav_frac1(jpart,ks)=wetscav*(zpoint2(jpart)-zpoint1(jpart))443 ! write (*,*) 'paricle kept: ',jpart,ks,wetdeposit(ks),xscav_frac1(jpart,ks),xmass(jpart,ks)444 endif445 endif446 endif447 448 449 end do !all species450 451 138 ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs 452 139 ! Add the wet deposition to accumulated amount on output grid and nested output grid … … 467 154 tot_inc_count=tot_inc_count+inc_count 468 155 469 ! write (*,*) 'First time removed: ',firsttimerem470 471 156 end subroutine wetdepo
Note: See TracChangeset
for help on using the changeset viewer.