Changeset db712a8 in flexpart.git for src/wetdepo.f90
- Timestamp:
- Mar 3, 2016, 12:34:56 PM (8 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 38b7917
- Parents:
- b0434e1
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/wetdepo.f90
r41d8574 rdb712a8 21 21 22 22 subroutine wetdepo(itime,ltsample,loutnext) 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 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 65 66 66 use point_mod … … 71 71 72 72 integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy 73 integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v 73 integer :: ngrid,itage,nage,hz,il,interp_time, n 74 integer(kind=1) :: clouds_v 74 75 integer :: ks, kp 75 76 ! integer :: n1,n2, icbot,ictop, indcloud !TEST … … 79 80 real :: wetdeposit(maxspec),restmass 80 81 real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled 81 82 !save lfr,cfr 82 83 83 84 real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/) … … 91 92 integer :: blc_count, inc_count 92 93 real :: Si_dummy, wetscav_dummy 93 94 95 96 97 94 logical :: readclouds_this_nest 95 96 97 ! Compute interval since radioactive decay of deposited mass was computed 98 !************************************************************************ 98 99 99 100 if (itime.le.loutnext) then … … 103 104 endif 104 105 105 106 106 ! Loop over all particles 107 !************************ 107 108 108 109 blc_count=0 … … 118 119 endif 119 120 120 121 ! Determine age class of the particle 121 122 itage=abs(itra1(jpart)-itramem(jpart)) 122 123 do nage=1,nageclass … … 126 127 127 128 128 129 129 ! Determine which nesting level to be used 130 !***************************************** 130 131 131 132 ngrid=0 132 133 do j=numbnests,1,-1 133 134 if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. & 134 (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then135 (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then 135 136 ngrid=j 136 137 goto 23 … … 140 141 141 142 142 143 143 ! Determine nested grid coordinates 144 !********************************** 144 145 145 146 if (ngrid.gt.0) then … … 148 149 ix=int(xtn) 149 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 150 156 else 151 157 ix=int(xtra1(jpart)) … … 154 160 155 161 156 157 158 159 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 !******************************************************************** 160 166 interp_time=nint(itime-0.5*ltsample) 161 167 162 168 if (ngrid.eq.0) then 163 169 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, & 164 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &165 memtime(1),memtime(2),interp_time,lsp,convp,cc)170 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, & 171 memtime(1),memtime(2),interp_time,lsp,convp,cc) 166 172 else 167 173 call interpol_rain_nests(lsprecn,convprecn,tccn, & 168 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &169 memtime(1),memtime(2),interp_time,lsp,convp,cc)174 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 175 memtime(1),memtime(2),interp_time,lsp,convp,cc) 170 176 endif 171 177 … … 173 179 if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 174 180 175 181 ! get the level were the actual particle is in 176 182 do il=2,nz 177 183 if (height(il).gt.ztra1(jpart)) then 178 184 hz=il-1 179 goto 26 185 ! goto 26 186 exit 180 187 endif 181 188 end do 182 26 continue189 !26 continue 183 190 184 191 n=memind(2) … … 190 197 clouds_h=cloudsh(ix,jy,n) 191 198 else 192 ! new removal not implemented for nests yet193 199 clouds_v=cloudsn(ix,jy,hz,n,ngrid) 194 clouds_h=clouds nh(ix,jy,n,ngrid)195 endif 196 197 198 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 199 205 200 206 if (clouds_v.le.1) goto 20 201 202 203 204 205 206 207 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 !************************************************************************** 208 214 209 215 if (lsp.gt.20.) then … … 232 238 233 239 234 235 236 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 237 243 grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp)) 238 244 grfraction(2)=max(0.05,cc*(lfr(i))) … … 240 246 241 247 242 243 248 ! 2) Computation of precipitation rate in sub-grid cell 249 !****************************************************** 244 250 prec(1)=(lsp+convp)/grfraction(1) 245 251 prec(2)=(lsp)/grfraction(2) … … 247 253 248 254 249 250 251 255 ! 3) Computation of scavenging coefficients for all species 256 ! Computation of wet deposition 257 !********************************************************** 252 258 253 259 do ks=1,nspec ! loop over species … … 255 261 wetscav=0. 256 262 257 263 !ZHG test if it nested? 258 264 if (ngrid.gt.0) then 259 265 act_temp=ttn(ix,jy,hz,n,ngrid) … … 261 267 act_temp=tt(ix,jy,hz,n) 262 268 endif 263 264 265 266 267 269 270 271 !*********************** 272 ! BELOW CLOUD SCAVENGING 273 !*********************** 268 274 if (clouds_v.ge.4) then !below cloud 269 275 … … 274 280 if (dquer(ks) .le. 0.) then !is gas 275 281 wetscav=weta(ks)*prec(1)**wetb(ks) 276 282 277 283 !AEROSOL 278 284 else !is particle … … 281 287 ! for particles larger than 10 um use the largest size defined in the parameterizations (10um) 282 288 dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m 283 if (act_temp .ge. 273 .and. weta(ks).gt.0.) then !Rain284 285 286 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 287 293 wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* & 288 294 (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5)) 289 295 290 elseif (act_temp .lt. 273 .and. wetb(ks).gt.0.) then ! Snow291 292 293 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 294 300 wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* & 295 301 (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5)) … … 304 310 305 311 306 307 308 !******************************************************312 !******************** 313 ! IN CLOUD SCAVENGING 314 !******************** 309 315 if (clouds_v.lt.4) then ! In-cloud 310 316 ! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file … … 315 321 if (wetb_in(ks).lt.0.) wetb_in(ks)=0. 316 322 317 !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF 318 if (readclouds) then !icloud_stats(ix,jy,4,n) has units kg/m2 319 ! cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc) 320 ! ESO: stop using icloud_stats 321 cl =clw4(ix,jy,n)*(grfraction(1)/cc) 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) 322 329 else !parameterize cloudwater m2/m3 323 330 !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF 324 331 cl=1.6E-6*prec(1)**0.36 325 332 endif 326 333 327 328 if (act_temp .le. 253) then329 330 else if (act_temp .ge. 273) then331 332 333 liq_frac =((act_temp-273)/(273-253))**2334 endif335 336 337 338 339 340 341 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 !************************************************** 342 349 if (dquer(ks).gt. 0.) then ! is particle 343 350 344 351 S_i= frac_act/cl 345 352 346 347 353 !********************* 354 ! GAS 348 355 else ! is gas 349 356 350 357 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 351 352 358 !REPLACE to switch old/ new scheme 359 ! S_i=frac_act/cle 353 360 S_i=1/cle 354 361 endif ! gas or particle 355 362 356 357 358 if ( readclouds) then363 ! 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 359 366 wetscav=S_i*(prec(1)/3.6E6) 360 367 else … … 362 369 endif 363 370 364 !ZHG 2015 TEST365 ! Si_dummy=frac_act/2E-7*prec(1)**0.36366 ! wetscav_dummy=Si_dummy*(prec(1)/3.6E6)/clouds_h367 ! if (clouds_v.lt.4) then368 ! talltest=talltest+1369 !if(talltest .eq. 1) OPEN(UNIT=199, FILE=utfil,FORM='FORMATTED',STATUS = 'UNKNOWN')370 !if(talltest .lt. 100001) write(199,*) prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90371 !if(talltest .lt. 100001) write(199,*) wetscav, wetscav_dummy,prec(1),ytra1(jpart)-90,clouds_v,cl372 !if(talltest .eq. 100001) CLOSE(199)373 !if(talltest .eq. 100001) STOP374 !375 !write(*,*) 'PREC kg/m2s CLOUD kg/m2', (prec(1)/3.6E6), cl !, '2E-7*prec(1)**0.36', 2E-7*prec(1)**0.36,'2E-7*prec(1)**0.36*clouds_h',2E-7*prec(1)**0.36*clouds_h376 !write(*,*) 'PREC kg/m2s LSP+convp kg/m2', prec(1), convp+lsp377 !write(*,*) wetscav, wetscav_dummy378 !write(*,*) cc, grfraction(1), cc/grfraction(1)379 380 !write(*,*) 'Lmbda_old', (prec(1)/3.6E6)/(clouds_h*2E-7*prec(1)**0.36)381 382 383 !write(*,*) '**************************************************'384 !write(*,*) 'clouds_h', clouds_h, 'clouds_v',clouds_v,'abs(ltsample)', abs(ltsample)385 !write(*,*) 'readclouds', readclouds, 'wetscav',wetscav, 'wetscav_dummy', wetscav_dummy386 !write(*,*) 'S_i', S_i , 'Si_dummy', Si_dummy, 'prec(1)', prec(1)387 388 389 ! write(*,*) 'PRECIPITATION ,cl ECMWF , cl PARAMETIZED, clouds_v, lat' &390 ! ,prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90391 392 !endif393 371 394 372 endif ! positive in-cloud scavenging parameters given in Species file 395 373 endif !incloud 396 374 !END ZHG TEST 397 398 399 400 401 402 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 403 381 404 382 if (wetscav.gt.0.) then … … 421 399 if (restmass .gt. smallnum) then 422 400 xmass1(jpart,ks)=restmass 423 424 425 401 ! depostatistic 402 ! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 403 ! depostatistic 426 404 else 427 405 xmass1(jpart,ks)=0. 428 406 endif 429 430 407 ! Correct deposited mass to the last time step when radioactive decay of 408 ! gridded deposited mass was calculated 431 409 if (decay(ks).gt.0.) then 432 410 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks)) … … 436 414 end do !all species 437 415 438 439 440 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 !***************************************************************************** 441 419 442 420 if (ldirect.eq.1) then … … 450 428 end do ! all particles 451 429 452 430 ! count the total number of below-cloud and in-cloud occurences: 453 431 tot_blc_count=tot_blc_count+blc_count 454 432 tot_inc_count=tot_inc_count+inc_count
Note: See TracChangeset
for help on using the changeset viewer.