- Timestamp:
- Sep 30, 2016, 11:01:54 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:
- 9669e1e
- Parents:
- dced13c
- Location:
- src
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
src/com_mod.f90
r462f74b r54cbd6c 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,SCAVDEP 578 logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC 579 logical :: DRYBKDEP,WETBKDEP 579 580 580 581 ! numxgrid,numygrid number of grid points in x,y-direction … … 595 596 ! OHREA .true., if OH reaction is switched on 596 597 ! ASSSPEC .true., if there are two species asscoiated 597 ! SCAVDEP .true., for bkwd runs, where mass deposited and source regions is calculated598 ! DRYBKDEP,WETBKDEP .true., for bkwd runs, where mass deposited and source regions is calculated - either for dry or for wet deposition 598 599 ! (i.e. transfer of mass between these two occurs 599 600 -
src/conccalc.f90
r462f74b r54cbd6c 192 192 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & 193 193 (jy.le.numygrid-1)) then 194 if ( SCAVDEP) then194 if (DRYBKDEP.or.WETBKDEP) then 195 195 do ks=1,nspec 196 196 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & … … 234 234 if ((jy.ge.0).and.(jy.le.numygrid-1)) then 235 235 w=wx*wy 236 if ( SCAVDEP) then236 if (DRYBKDEP.or.WETBKDEP) then 237 237 do ks=1,nspec 238 238 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & … … 251 251 if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then 252 252 w=wx*(1.-wy) 253 if ( SCAVDEP) then253 if (DRYBKDEP.or.WETBKDEP) then 254 254 do ks=1,nspec 255 255 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & … … 271 271 if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then 272 272 w=(1.-wx)*(1.-wy) 273 if ( SCAVDEP) then273 if (DRYBKDEP.or.WETBKDEP) then 274 274 do ks=1,nspec 275 275 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & … … 288 288 if ((jy.ge.0).and.(jy.le.numygrid-1)) then 289 289 w=(1.-wx)*wy 290 if ( SCAVDEP) then290 if (DRYBKDEP.or.WETBKDEP) then 291 291 do ks=1,nspec 292 292 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & … … 329 329 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & 330 330 (jy.le.numygridn-1)) then 331 if ( SCAVDEP) then331 if (DRYBKDEP.or.WETBKDEP) then 332 332 do ks=1,nspec 333 333 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & … … 371 371 if ((jy.ge.0).and.(jy.le.numygridn-1)) then 372 372 w=wx*wy 373 if ( SCAVDEP) then373 if (DRYBKDEP.or.WETBKDEP) then 374 374 do ks=1,nspec 375 375 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & … … 388 388 if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then 389 389 w=wx*(1.-wy) 390 if ( SCAVDEP) then390 if (DRYBKDEP.or.WETBKDEP) then 391 391 do ks=1,nspec 392 392 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & … … 408 408 if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then 409 409 w=(1.-wx)*(1.-wy) 410 if ( SCAVDEP) then410 if (DRYBKDEP.or.WETBKDEP) then 411 411 do ks=1,nspec 412 412 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & … … 425 425 if ((jy.ge.0).and.(jy.le.numygridn-1)) then 426 426 w=(1.-wx)*wy 427 if ( SCAVDEP) then427 if (DRYBKDEP.or.WETBKDEP) then 428 428 do ks=1,nspec 429 429 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & -
src/concoutput.f90
r462f74b r54cbd6c 246 246 247 247 write(anspec,'(i3.3)') ks 248 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 249 if (ldirect.eq.1) then 250 open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// & 251 atime//'_'//anspec,form='unformatted') 252 else 253 open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// & 254 atime//'_'//anspec,form='unformatted') 255 endif 256 write(unitoutgrid) itime 257 endif 258 259 if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio 260 open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// & 248 249 if (DRYBKDEP.or.WETBKDEP) then !scavdep output 250 if (DRYBKDEP) & 251 open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// & 261 252 atime//'_'//anspec,form='unformatted') 262 263 write(unitoutgridppt) itime 264 endif 265 266 if (iout.eq.6) then !scavdep output 267 open(unitoutgrid,file=path(2)(1:length(2))//'grid_scav_'//adate// & 253 if (WETBKDEP) & 254 open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// & 268 255 atime//'_'//anspec,form='unformatted') 269 256 write(unitoutgrid) itime 270 endif 257 else 258 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 259 if (ldirect.eq.1) then 260 open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// & 261 atime//'_'//anspec,form='unformatted') 262 else 263 open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// & 264 atime//'_'//anspec,form='unformatted') 265 endif 266 write(unitoutgrid) itime 267 endif 268 if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio 269 open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// & 270 atime//'_'//anspec,form='unformatted') 271 write(unitoutgridppt) itime 272 endif 273 endif ! if deposition output 271 274 272 275 do kp=1,maxpointspec_act -
src/readcommand.f90
r462f74b r54cbd6c 296 296 !Af 1 = mass units 297 297 !Af 2 = mass mixing ratio units 298 !Af 3 = wet deposition in outputfield 299 !Af 4 = dry deposition in outputfield 298 300 299 301 if ( ldirect .eq. 1 ) then ! FWD-Run … … 318 320 endif 319 321 !Af set release-switch 320 if (ind_receptor .eq. 1) then !mass 322 WETBKDEP=.false. 323 DRYBKDEP=.false. 324 select case (ind_receptor) 325 case (1) ! 1 .. concentration at receptor 321 326 ind_rel = 1 322 else ! mass mix327 case (2) ! 2 .. mixing ratio at receptor 323 328 ind_rel = 0 324 endif 329 case (3) ! 3 .. wet deposition in outputfield 330 ind_rel = 3 331 write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE #### ' 332 WETBKDEP=.true. 333 allocate(xscav_frac1(maxpart,maxspec)) 334 case (4) ! 4 .. dry deposition in outputfield 335 ind_rel = 4 336 write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE #### ' 337 DRYBKDEP=.true. 338 allocate(xscav_frac1(maxpart,maxspec)) 339 end select 325 340 endif 326 341 … … 338 353 stop 339 354 endif 340 341 if ((ldirect.eq.-1).and.(iout.eq.6)) then342 if ((ind_receptor .eq. 1) .and. (ind_source .eq. 1)) then343 write(*,*) ' #### FLEXPART SCAVENGING DEPOSIT BACKWARD MODE #### '344 SCAVDEP=.true.345 allocate(xscav_frac1(maxpart,maxspec))346 else347 write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'348 write(*,*) '#### FOR SCAVDEP MODE ind_source and ####'349 write(*,*) '#### ind_receptor have to be 1 ! ####'350 stop351 endif352 else353 SCAVDEP=.false.354 endif355 355 356 356 ! Check input dates … … 401 401 if ((iout.lt.1).or.(iout.gt.6)) then 402 402 write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' 403 write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 , 5 OR 6 FOR#### '403 write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR #### ' 404 404 write(*,*) ' #### STANDARD FLEXPART OUTPUT OR 9 - 13 #### ' 405 405 write(*,*) ' #### FOR NETCDF OUTPUT #### ' -
src/releaseparticles.f90
r462f74b r54cbd6c 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)60 59 real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold 61 60 real :: presspart,average_timecorrect … … 87 86 minpart=1 88 87 do i=1,numpoint 89 do k=1,nspec90 rhosum(k)=091 end do92 88 if ((itime.ge.ireleasestart(i)).and. &! are we within release interval? 93 89 (itime.le.ireleaseend(i))) then … … 186 182 xmass1(ipart,k)=xmass(i,k)/real(npart(i)) & 187 183 *timecorrect(k)/average_timecorrect 188 if ( SCAVDEP) then ! if there is no scavenging in wetdepo it will be set to 0184 if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0 189 185 ! if ( henry(k).gt.0 .or. & 190 186 ! crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. & 191 187 ! ccn_aero(k).gt.0. .or. in_aero(k).gt.0. ) then 192 188 xscav_frac1(ipart,k)=-1. 193 ! write(*,*) '190: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i194 ! xscav_frac1(ipart,k)=(-1.)/real(npart(i)) &195 ! *timecorrect(k)/average_timecorrect196 ! else197 ! xscav_frac1(ipart,k)=0198 ! endif199 189 endif 200 190 ! Assign certain properties to particle … … 389 379 do k=1,nspec 390 380 xmass1(ipart,k)=xmass1(ipart,k)*rhoout 391 if (SCAVDEP) then392 xscav_frac1(ipart,k)=xscav_frac1(ipart,k)393 !mctest xscav_frac1(ipart,k)=xscav_frac1(ipart,k)*rhoout394 rhosum(k)=rhosum(k)+rhoout395 ! write(*,*) '391: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i396 endif397 381 end do 398 382 endif 399 400 383 401 384 numpart=max(numpart,ipart) … … 407 390 34 minpart=ipart+1 408 391 end do ! ipart=minpart,maxpart 409 if (SCAVDEP) then410 do ipart=minpart,maxpart411 do k=1,nspec412 if (xscav_frac1(ipart,k).lt.0) then413 !mctest xscav_frac1(ipart,k)=xscav_frac1(ipart,k)/rhosum(k)414 ! write(*,*) '409: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i415 endif416 end do417 end do418 endif419 392 endif ! j=1,numrel 420 393 end do -
src/timemanager.f90
r462f74b r54cbd6c 115 115 real :: xold,yold,zold,xmassfract 116 116 real, parameter :: e_inv = 1.0/exp(1.0) 117 logical :: firstdepocalc 117 118 !double precision xm(maxspec,maxpointspec_act), 118 119 ! + xm_depw(maxspec,maxpointspec_act), … … 171 172 write (*,*) 'timemanager> call wetdepo' 172 173 endif 173 call wetdepo(itime,lsynctime,loutnext )174 call wetdepo(itime,lsynctime,loutnext,.false.) 174 175 endif 175 176 … … 541 542 zold=ztra1(j) 542 543 544 545 ! RECEPTOR: dry/wet depovel 546 !**************************** 547 ! Before the particle is moved 548 ! the calculation of the scavenged mass shall only be done once after release 549 ! xscav_frac1 was initialised with a negative value 550 do ks=1,nspec 551 if (DRYBKDEP.and.(xscav_frac1(j,ks).lt.0)) then 552 if (ks.eq.1) then 553 call advance_rec(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), & 554 us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, & 555 cbt(j)) 556 endif 557 if (decay(ks).gt.0.) then ! radioactive decay 558 decfact=exp(-real(abs(lsynctime))*decay(ks)) 559 else 560 decfact=1. 561 endif 562 if (DRYDEPSPEC(ks)) then ! dry deposition 563 drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact 564 xscav_frac1(j,ks)=xscav_frac1(j,ks)*(-1.)* & 565 drydeposit(ks)/xmass1(j,ks) 566 if (decay(ks).gt.0.) then ! correct for decay (see wetdepo) 567 drydeposit(ks)=drydeposit(ks)* & 568 exp(real(abs(ldeltat))*decay(ks)) 569 endif 570 else 571 xmass1(j,ks)=0 572 xscav_frac1(j,ks)=0. 573 endif 574 endif 575 enddo 576 577 firstdepocalc=.false. 578 do ks=1,nspec 579 if ((WETBKDEP).and.(xscav_frac1(j,ks).lt.0) & 580 .and.firstdepocalc.eqv..false.) then 581 ! Backward wetdeposition and first timestep after release 582 call wetdepo(itime,lsynctime,loutnext,.true.) 583 firstdepocalc=.true. 584 endif 585 enddo 586 543 587 ! Integrate Lagevin equation for lsynctime seconds 544 588 !************************************************* 545 589 546 call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), & 547 us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, & 548 cbt(j)) 590 if (verbosity.gt.0) then 591 if (j.eq.1) then 592 write (*,*) 'timemanager> call advance' 593 endif 594 endif 595 call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), & 596 us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, & 597 cbt(j)) 549 598 550 599 ! Calculate the gross fluxes across layer interfaces -
src/wetdepo.f90
r462f74b r54cbd6c 20 20 !********************************************************************** 21 21 22 subroutine wetdepo(itime,ltsample,loutnext )23 ! i i i 22 subroutine wetdepo(itime,ltsample,loutnext,forreceptor) 23 ! i i i i 24 24 !***************************************************************************** 25 25 ! * … … 30 30 ! This fraction is parameterized from total cloud cover and rates of large * 31 31 ! scale and convective precipitation. * 32 ! SEC: if forrecptor is true then the wetdeposition fraction is only applied * 33 ! on the xscav_frac and not on the xmass * 32 34 ! * 33 35 ! Author: A. Stohl * … … 92 94 integer :: blc_count, inc_count 93 95 real :: Si_dummy, wetscav_dummy 94 logical :: readclouds_this_nest 96 logical :: readclouds_this_nest,forreceptor 95 97 96 98 … … 172 174 memtime(1),memtime(2),interp_time,lsp,convp,cc) 173 175 endif 174 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 (SCAVDEP) 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 176 191 177 ! get the level were the actual particle is in … … 413 399 kp=1 414 400 endif 415 if (restmass .gt. smallnum) then 416 xmass1(jpart,ks)=restmass 401 if (forreceptor .eqv. .false.) then 402 if (restmass .gt. smallnum) then 403 xmass1(jpart,ks)=restmass 417 404 ! depostatistic 418 405 ! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 419 406 ! depostatistic 420 else 421 xmass1(jpart,ks)=0. 422 endif 407 else 408 xmass1(jpart,ks)=0. 409 endif 410 else ! for the backward deposition calculation 411 if (wetdeposit(ks).gt.0) then ! deposition occured 412 xscav_frac1(jpart,ks)=xscav_frac1(jpart,ks)*(-1.)* & 413 wetdeposit(ks)/xmass1(jpart,ks) 414 ! write (*,*) 'paricle kept: ',jpart,ks,wetdeposit(ks),xscav_frac1(jpart,ks) 415 else 416 xmass1(jpart,ks)=0. 417 xscav_frac1(jpart,ks)=0. 418 endif 419 endif 423 420 ! Correct deposited mass to the last time step when radioactive decay of 424 421 ! gridded deposited mass was calculated … … 426 423 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks)) 427 424 endif 428 429 if (SCAVDEP) 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 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 endif442 endif443 endif444 445 425 446 426 end do !all species
Note: See TracChangeset
for help on using the changeset viewer.