Changeset 462f74b in flexpart.git for src


Ignore:
Timestamp:
Jul 12, 2016, 11:02:42 AM (8 years ago)
Author:
Sabine Eckhardt <sabine@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
842074e
Parents:
f28aa0a
Message:

first version of the backward scavenging

Location:
src
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • src/com_mod.f90

    r62e65c7 r462f74b  
    576576  real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
    577577  !real outheight(maxzgrid),outheighthalf(maxzgrid)
    578   logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC
     578  logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC,SCAVDEP
    579579
    580580  ! numxgrid,numygrid       number of grid points in x,y-direction
     
    595595  ! OHREA                   .true., if OH reaction is switched on
    596596  ! ASSSPEC                 .true., if there are two species asscoiated
     597  ! SCAVDEP                 .true., for bkwd runs, where mass deposited and source regions is calculated
    597598  !                    (i.e. transfer of mass between these two occurs
    598599
     
    665666  real, allocatable, dimension(:) :: ztra1
    666667  real, allocatable, dimension(:,:) :: xmass1
     668  real, allocatable, dimension(:,:) :: xscav_frac1
    667669
    668670  ! eso: Moved from timemanager
     
    685687  ! xtra1,ytra1,ztra1       spatial positions of the particles
    686688  ! xmass1 [kg]             particle masses
    687  
     689  ! xscav_frac1             fraction of particle masse which has been scavenged at receptor
     690 
    688691
    689692
  • src/conccalc.f90

    rfdc0f03 r462f74b  
    2121
    2222subroutine conccalc(itime,weight)
    23   !                 i     i
     23  !                      i     i
    2424  !*****************************************************************************
    2525  !                                                                            *
     
    6161
    6262
     63  integer :: usekernel
     64
     65  usekernel=0
     66  if (usekernel.ne.1) then
     67     write (*,*) 'NOT USING THE KERNEL!'
     68  endif
    6369  ! For forward simulations, make a loop over the number of species;
    6470  ! for backward simulations, make an additional loop over the
     
    181187  !*****************************************************************************
    182188
    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. &
    184190           (xl.gt.real(numxgrid-1)-0.5).or. &
    185            (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
     191           (yl.gt.real(numygrid-1)-0.5)).or.(usekernel.eq.0)) then             ! no kernel, direct attribution to grid cell
    186192        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    187193             (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)= &
    190203                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    191204                 xmass1(i,ks)/rhoi*weight
    192           end do
     205             end do
     206          endif
    193207        endif
    194208
     
    220234          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    221235            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)= &
    224245                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    225246                   xmass1(i,ks)/rhoi*weight*w
    226             end do
     247               end do
     248            endif
    227249          endif
    228250
    229251          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    230252            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)= &
    233262                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    234263                   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
    238268
    239269
     
    241271          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    242272            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)= &
    245282                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    246283                   xmass1(i,ks)/rhoi*weight*w
    247             end do
     284               end do
     285            endif
    248286          endif
    249287
    250288          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    251289            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)= &
    254299                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    255300                   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
    262306
    263307  !************************************
     
    285329          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
    286330               (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)= &
    289340                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    290341                   xmass1(i,ks)/rhoi*weight
    291             end do
     342               end do
     343            endif
    292344          endif
    293345
     
    319371            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    320372              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)= &
    323382                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    324383                     xmass1(i,ks)/rhoi*weight*w
    325               end do
     384                 end do
     385              endif
    326386            endif
    327387
    328388            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    329389              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)= &
    332399                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    333400                     xmass1(i,ks)/rhoi*weight*w
    334               end do
     401                 end do
     402              endif
    335403            endif
    336404          endif
     
    340408            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    341409              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)= &
    344419                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    345420                     xmass1(i,ks)/rhoi*weight*w
    346               end do
     421                 end do
     422              endif
    347423            endif
    348424
    349425            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    350426              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)= &
    353436                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    354437                     xmass1(i,ks)/rhoi*weight*w
    355               end do
     438                 end do
     439              endif
    356440            endif
    357441          endif
    358442        endif
    359 
    360443      endif
    361444    endif
  • src/concoutput.f90

    rdb712a8 r462f74b  
    264264    endif
    265265
     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
    266272    do kp=1,maxpointspec_act
    267273      do nage=1,nageclass
     
    342348! Concentration output
    343349!*********************
    344         if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     350        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
    345351
    346352! Wet deposition
  • src/readcommand.f90

    rc04b739 r462f74b  
    339339  endif
    340340
     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
    341356  ! Check input dates
    342357  !******************
     
    384399  !**********************************************************************
    385400
    386   if ((iout.lt.1).or.(iout.gt.5)) then
     401  if ((iout.lt.1).or.(iout.gt.6)) then
    387402    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     #### '
    389404    write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
    390405    write(*,*) ' #### FOR NETCDF OUTPUT                       #### '
  • src/releaseparticles.f90

    r8a65cb0 r462f74b  
    5757  real :: xaux,yaux,zaux,rfraction
    5858  real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4
     59  real :: rhosum(nspec)
    5960  real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold
    6061  real :: presspart,average_timecorrect
     
    8687  minpart=1
    8788  do i=1,numpoint
     89    do k=1,nspec
     90      rhosum(k)=0
     91    end do
    8892    if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
    8993         (itime.le.ireleaseend(i))) then
     
    176180  ! scaled. Adjust the mass per particle by the species-dependent time correction factor
    177181  ! 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.
    178184  !*****************************************************************************
    179185            do k=1,nspec
    180186               xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
    181187                    *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
    183200  ! Assign certain properties to particle
    184201  !**************************************
     
    372389              do k=1,nspec
    373390                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
    374397              end do
    375398            endif
     
    379402            goto 34      ! Storage space has been found, stop searching
    380403          endif
    381         end do
     404        end do  ! i=1:numpoint
    382405        if (ipart.gt.maxpart) goto 996
    383406
    38440734      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
    386418      endif
     419      endif ! j=1,numrel
    387420  end do
    388421
  • src/timemanager.f90

    r18adf60 r462f74b  
    376376
    377377      if ((itime.eq.loutend).and.(outnum.gt.0.)) then
    378         if ((iout.le.3.).or.(iout.eq.5)) then
     378        if ((iout.le.3.).or.(iout.eq.5).or.(iout.eq.6)) then
    379379          if (surf_only.ne.1) then
    380380            if (lnetcdfout.eq.1) then
  • src/wetdepo.f90

    r05cf28d r462f74b  
    174174
    175175!  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
    177190
    178191! get the level were the actual particle is in
     
    414427      endif
    415428
     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
    416445
    417446    end do !all species
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG