- Timestamp:
- Jul 12, 2016, 11:02:42 AM (8 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 842074e
- Parents:
- f28aa0a
- Location:
- src
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
src/com_mod.f90
r62e65c7 r462f74b 576 576 real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn 577 577 !real outheight(maxzgrid),outheighthalf(maxzgrid) 578 logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC 578 logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC,SCAVDEP 579 579 580 580 ! numxgrid,numygrid number of grid points in x,y-direction … … 595 595 ! OHREA .true., if OH reaction is switched on 596 596 ! ASSSPEC .true., if there are two species asscoiated 597 ! SCAVDEP .true., for bkwd runs, where mass deposited and source regions is calculated 597 598 ! (i.e. transfer of mass between these two occurs 598 599 … … 665 666 real, allocatable, dimension(:) :: ztra1 666 667 real, allocatable, dimension(:,:) :: xmass1 668 real, allocatable, dimension(:,:) :: xscav_frac1 667 669 668 670 ! eso: Moved from timemanager … … 685 687 ! xtra1,ytra1,ztra1 spatial positions of the particles 686 688 ! xmass1 [kg] particle masses 687 689 ! xscav_frac1 fraction of particle masse which has been scavenged at receptor 690 688 691 689 692 -
src/conccalc.f90
rfdc0f03 r462f74b 21 21 22 22 subroutine conccalc(itime,weight) 23 ! i i23 ! i i 24 24 !***************************************************************************** 25 25 ! * … … 61 61 62 62 63 integer :: usekernel 64 65 usekernel=0 66 if (usekernel.ne.1) then 67 write (*,*) 'NOT USING THE KERNEL!' 68 endif 63 69 ! For forward simulations, make a loop over the number of species; 64 70 ! for backward simulations, make an additional loop over the … … 181 187 !***************************************************************************** 182 188 183 if (( itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &189 if (((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. & 184 190 (xl.gt.real(numxgrid-1)-0.5).or. & 185 (yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell191 (yl.gt.real(numygrid-1)-0.5)).or.(usekernel.eq.0)) then ! no kernel, direct attribution to grid cell 186 192 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & 187 193 (jy.le.numygrid-1)) then 188 do ks=1,nspec 189 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 194 if (SCAVDEP) then 195 do ks=1,nspec 196 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 197 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 198 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*xscav_frac1(i,ks) 199 end do 200 else 201 do ks=1,nspec 202 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 190 203 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 191 204 xmass1(i,ks)/rhoi*weight 192 end do 205 end do 206 endif 193 207 endif 194 208 … … 220 234 if ((jy.ge.0).and.(jy.le.numygrid-1)) then 221 235 w=wx*wy 222 do ks=1,nspec 223 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 236 if (SCAVDEP) then 237 do ks=1,nspec 238 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 239 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 240 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*w*weight*xscav_frac1(i,ks) 241 end do 242 else 243 do ks=1,nspec 244 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 224 245 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 225 246 xmass1(i,ks)/rhoi*weight*w 226 end do 247 end do 248 endif 227 249 endif 228 250 229 251 if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then 230 252 w=wx*(1.-wy) 231 do ks=1,nspec 232 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 253 if (SCAVDEP) then 254 do ks=1,nspec 255 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 256 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 257 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks) 258 end do 259 else 260 do ks=1,nspec 261 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 233 262 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 234 263 xmass1(i,ks)/rhoi*weight*w 235 end do 236 endif 237 endif 264 end do 265 endif 266 endif 267 endif !ix ge 0 238 268 239 269 … … 241 271 if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then 242 272 w=(1.-wx)*(1.-wy) 243 do ks=1,nspec 244 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 273 if (SCAVDEP) then 274 do ks=1,nspec 275 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 276 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 277 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*w*weight*xscav_frac1(i,ks) 278 end do 279 else 280 do ks=1,nspec 281 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 245 282 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 246 283 xmass1(i,ks)/rhoi*weight*w 247 end do 284 end do 285 endif 248 286 endif 249 287 250 288 if ((jy.ge.0).and.(jy.le.numygrid-1)) then 251 289 w=(1.-wx)*wy 252 do ks=1,nspec 253 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & 290 if (SCAVDEP) then 291 do ks=1,nspec 292 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & 293 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 294 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks) 295 end do 296 else 297 do ks=1,nspec 298 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & 254 299 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 255 300 xmass1(i,ks)/rhoi*weight*w 256 end do 257 endif 258 endif 259 endif 260 261 301 end do 302 endif 303 endif 304 endif !ixp ge 0 305 endif 262 306 263 307 !************************************ … … 285 329 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & 286 330 (jy.le.numygridn-1)) then 287 do ks=1,nspec 288 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 331 if (SCAVDEP) then 332 do ks=1,nspec 333 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 334 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 335 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*xscav_frac1(i,ks) 336 end do 337 else 338 do ks=1,nspec 339 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 289 340 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 290 341 xmass1(i,ks)/rhoi*weight 291 end do 342 end do 343 endif 292 344 endif 293 345 … … 319 371 if ((jy.ge.0).and.(jy.le.numygridn-1)) then 320 372 w=wx*wy 321 do ks=1,nspec 322 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 373 if (SCAVDEP) then 374 do ks=1,nspec 375 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 376 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 377 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks) 378 end do 379 else 380 do ks=1,nspec 381 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & 323 382 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 324 383 xmass1(i,ks)/rhoi*weight*w 325 end do 384 end do 385 endif 326 386 endif 327 387 328 388 if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then 329 389 w=wx*(1.-wy) 330 do ks=1,nspec 331 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 390 if (SCAVDEP) then 391 do ks=1,nspec 392 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 393 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 394 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks) 395 end do 396 else 397 do ks=1,nspec 398 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 332 399 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 333 400 xmass1(i,ks)/rhoi*weight*w 334 end do 401 end do 402 endif 335 403 endif 336 404 endif … … 340 408 if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then 341 409 w=(1.-wx)*(1.-wy) 342 do ks=1,nspec 343 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 410 if (SCAVDEP) then 411 do ks=1,nspec 412 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 413 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 414 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks) 415 end do 416 else 417 do ks=1,nspec 418 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & 344 419 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & 345 420 xmass1(i,ks)/rhoi*weight*w 346 end do 421 end do 422 endif 347 423 endif 348 424 349 425 if ((jy.ge.0).and.(jy.le.numygridn-1)) then 350 426 w=(1.-wx)*wy 351 do ks=1,nspec 352 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & 427 if (SCAVDEP) then 428 do ks=1,nspec 429 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & 430 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 431 xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks) 432 end do 433 else 434 do ks=1,nspec 435 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & 353 436 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ & 354 437 xmass1(i,ks)/rhoi*weight*w 355 end do 438 end do 439 endif 356 440 endif 357 441 endif 358 442 endif 359 360 443 endif 361 444 endif -
src/concoutput.f90
rdb712a8 r462f74b 264 264 endif 265 265 266 if (iout.eq.6) then !scavdep output 267 open(unitoutgrid,file=path(2)(1:length(2))//'grid_scav_'//adate// & 268 atime//'_'//anspec,form='unformatted') 269 write(unitoutgrid) itime 270 endif 271 266 272 do kp=1,maxpointspec_act 267 273 do nage=1,nageclass … … 342 348 ! Concentration output 343 349 !********************* 344 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5) ) then350 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then 345 351 346 352 ! Wet deposition -
src/readcommand.f90
rc04b739 r462f74b 339 339 endif 340 340 341 if ((ldirect.eq.-1).and.(iout.eq.6)) then 342 if ((ind_receptor .eq. 1) .and. (ind_source .eq. 1)) then 343 write(*,*) ' #### FLEXPART SCAVENGING DEPOSIT BACKWARD MODE #### ' 344 SCAVDEP=.true. 345 allocate(xscav_frac1(maxpart,maxspec)) 346 else 347 write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' 348 write(*,*) '#### FOR SCAVDEP MODE ind_source and ####' 349 write(*,*) '#### ind_receptor have to be 1 ! ####' 350 stop 351 endif 352 else 353 SCAVDEP=.false. 354 endif 355 341 356 ! Check input dates 342 357 !****************** … … 384 399 !********************************************************************** 385 400 386 if ((iout.lt.1).or.(iout.gt. 5)) then401 if ((iout.lt.1).or.(iout.gt.6)) then 387 402 write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' 388 write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR#### '403 write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, 5 OR 6 FOR #### ' 389 404 write(*,*) ' #### STANDARD FLEXPART OUTPUT OR 9 - 13 #### ' 390 405 write(*,*) ' #### FOR NETCDF OUTPUT #### ' -
src/releaseparticles.f90
r8a65cb0 r462f74b 57 57 real :: xaux,yaux,zaux,rfraction 58 58 real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4 59 real :: rhosum(nspec) 59 60 real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold 60 61 real :: presspart,average_timecorrect … … 86 87 minpart=1 87 88 do i=1,numpoint 89 do k=1,nspec 90 rhosum(k)=0 91 end do 88 92 if ((itime.ge.ireleasestart(i)).and. &! are we within release interval? 89 93 (itime.le.ireleaseend(i))) then … … 176 180 ! scaled. Adjust the mass per particle by the species-dependent time correction factor 177 181 ! divided by the species-average one 182 ! for the scavenging calculation the mass needs to be multiplied with rho of the particle layer and 183 ! divided by the sum of rho of all particles. 178 184 !***************************************************************************** 179 185 do k=1,nspec 180 186 xmass1(ipart,k)=xmass(i,k)/real(npart(i)) & 181 187 *timecorrect(k)/average_timecorrect 182 ! write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k 188 if (SCAVDEP) then ! if there is no scavenging in wetdepo it will be set to 0 189 ! if ( henry(k).gt.0 .or. & 190 ! crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. & 191 ! ccn_aero(k).gt.0. .or. in_aero(k).gt.0. ) then 192 xscav_frac1(ipart,k)=-1. 193 ! write(*,*) '190: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i 194 ! xscav_frac1(ipart,k)=(-1.)/real(npart(i)) & 195 ! *timecorrect(k)/average_timecorrect 196 ! else 197 ! xscav_frac1(ipart,k)=0 198 ! endif 199 endif 183 200 ! Assign certain properties to particle 184 201 !************************************** … … 372 389 do k=1,nspec 373 390 xmass1(ipart,k)=xmass1(ipart,k)*rhoout 391 if (SCAVDEP) then 392 xscav_frac1(ipart,k)=xscav_frac1(ipart,k) 393 !mctest xscav_frac1(ipart,k)=xscav_frac1(ipart,k)*rhoout 394 rhosum(k)=rhosum(k)+rhoout 395 ! write(*,*) '391: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i 396 endif 374 397 end do 375 398 endif … … 379 402 goto 34 ! Storage space has been found, stop searching 380 403 endif 381 end do 404 end do ! i=1:numpoint 382 405 if (ipart.gt.maxpart) goto 996 383 406 384 407 34 minpart=ipart+1 385 end do 408 end do ! ipart=minpart,maxpart 409 if (SCAVDEP) then 410 do ipart=minpart,maxpart 411 do k=1,nspec 412 if (xscav_frac1(ipart,k).lt.0) then 413 !mctest xscav_frac1(ipart,k)=xscav_frac1(ipart,k)/rhosum(k) 414 ! write(*,*) '409: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i 415 endif 416 end do 417 end do 386 418 endif 419 endif ! j=1,numrel 387 420 end do 388 421 -
src/timemanager.f90
r18adf60 r462f74b 376 376 377 377 if ((itime.eq.loutend).and.(outnum.gt.0.)) then 378 if ((iout.le.3.).or.(iout.eq.5) ) then378 if ((iout.le.3.).or.(iout.eq.5).or.(iout.eq.6)) then 379 379 if (surf_only.ne.1) then 380 380 if (lnetcdfout.eq.1) then -
src/wetdepo.f90
r05cf28d r462f74b 174 174 175 175 ! If total precipitation is less than 0.01 mm/h - no scavenging occurs 176 if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 176 ! sec this is just valid if release is over a point 177 if ((lsp.lt.0.01).and.(convp.lt.0.01)) then 178 if (SCAVDEP) then 179 do ks=1,nspec 180 if (xscav_frac1(jpart,ks).lt.0) then ! first timestep no scavenging 181 xmass1(jpart,ks)=0. 182 xscav_frac1(jpart,ks)=0. 183 ! write (*,*) 'paricle removed - no scavenging: ',jpart,ks 184 endif 185 end do 186 endif 187 goto 20 188 endif 189 177 190 178 191 ! get the level were the actual particle is in … … 414 427 endif 415 428 429 if (SCAVDEP) then 430 ! the calculation of the scavenged mass shall only be done once after release 431 ! xscav_frac1 was initialised with a negative value 432 if (xscav_frac1(jpart,ks).lt.0) then 433 if (wetdeposit(ks).eq.0) then 434 ! terminate particle 435 xmass1(jpart,ks)=0. 436 xscav_frac1(jpart,ks)=0. 437 else 438 xscav_frac1(jpart,ks)=xscav_frac1(jpart,ks)*(-1.)* & 439 wetdeposit(ks)/xmass1(jpart,ks) 440 ! write (*,*) 'paricle kept: ',jpart,ks,wetdeposit(ks),xscav_frac1(jpart,ks) 441 endif 442 endif 443 endif 444 416 445 417 446 end do !all species
Note: See TracChangeset
for help on using the changeset viewer.