Changeset 54cbd6c in flexpart.git for src


Ignore:
Timestamp:
Sep 30, 2016, 11:01:54 AM (7 years ago)
Author:
Sabine <sabine.eckhardt@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
9669e1e
Parents:
dced13c
Message:

all f90 files for dry/wet backward mode

Location:
src
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • src/com_mod.f90

    r462f74b r54cbd6c  
    576576  real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
    577577  !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
    579580
    580581  ! numxgrid,numygrid       number of grid points in x,y-direction
     
    595596  ! OHREA                   .true., if OH reaction is switched on
    596597  ! ASSSPEC                 .true., if there are two species asscoiated
    597   ! SCAVDEP                 .true., for bkwd runs, where mass deposited and source regions is calculated
     598  ! DRYBKDEP,WETBKDEP        .true., for bkwd runs, where mass deposited and source regions is calculated - either for dry or for wet deposition
    598599  !                    (i.e. transfer of mass between these two occurs
    599600
  • src/conccalc.f90

    r462f74b r54cbd6c  
    192192        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    193193             (jy.le.numygrid-1)) then
    194           if (SCAVDEP) then
     194          if (DRYBKDEP.or.WETBKDEP) then
    195195             do ks=1,nspec
    196196               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     
    234234          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    235235            w=wx*wy
    236             if (SCAVDEP) then
     236            if (DRYBKDEP.or.WETBKDEP) then
    237237               do ks=1,nspec
    238238                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     
    251251          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    252252            w=wx*(1.-wy)
    253             if (SCAVDEP) then
     253            if (DRYBKDEP.or.WETBKDEP) then
    254254              do ks=1,nspec
    255255                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     
    271271          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    272272            w=(1.-wx)*(1.-wy)
    273             if (SCAVDEP) then
     273            if (DRYBKDEP.or.WETBKDEP) then
    274274               do ks=1,nspec
    275275                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     
    288288          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    289289            w=(1.-wx)*wy
    290             if (SCAVDEP) then
     290            if (DRYBKDEP.or.WETBKDEP) then
    291291               do ks=1,nspec
    292292                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     
    329329          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
    330330               (jy.le.numygridn-1)) then
    331             if (SCAVDEP) then
     331            if (DRYBKDEP.or.WETBKDEP) then
    332332               do ks=1,nspec
    333333                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     
    371371            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    372372              w=wx*wy
    373               if (SCAVDEP) then
     373              if (DRYBKDEP.or.WETBKDEP) then
    374374                 do ks=1,nspec
    375375                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     
    388388            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    389389              w=wx*(1.-wy)
    390               if (SCAVDEP) then
     390              if (DRYBKDEP.or.WETBKDEP) then
    391391                 do ks=1,nspec
    392392                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     
    408408            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    409409              w=(1.-wx)*(1.-wy)
    410               if (SCAVDEP) then
     410              if (DRYBKDEP.or.WETBKDEP) then
    411411                 do ks=1,nspec
    412412                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     
    425425            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    426426              w=(1.-wx)*wy
    427               if (SCAVDEP) then
     427              if (DRYBKDEP.or.WETBKDEP) then
    428428                 do ks=1,nspec
    429429                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
  • src/concoutput.f90

    r462f74b r54cbd6c  
    246246
    247247    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// &
    261252           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// &
    268255           atime//'_'//anspec,form='unformatted')
    269256      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
    271274
    272275    do kp=1,maxpointspec_act
  • src/readcommand.f90

    r462f74b r54cbd6c  
    296296  !Af          1 = mass units
    297297  !Af          2 = mass mixing ratio units
     298  !Af          3 = wet deposition in outputfield
     299  !Af          4 = dry deposition in outputfield
    298300
    299301  if ( ldirect .eq. 1 ) then  ! FWD-Run
     
    318320     endif
    319321  !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
    321326        ind_rel = 1
    322      else ! mass mix
     327     case (2)  !  2 .. mixing ratio at receptor
    323328        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
    325340  endif
    326341
     
    338353    stop
    339354  endif
    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
    355355
    356356  ! Check input dates
     
    401401  if ((iout.lt.1).or.(iout.gt.6)) then
    402402    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        #### '
    404404    write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
    405405    write(*,*) ' #### FOR NETCDF OUTPUT                       #### '
  • src/releaseparticles.f90

    r462f74b r54cbd6c  
    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)
    6059  real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold
    6160  real :: presspart,average_timecorrect
     
    8786  minpart=1
    8887  do i=1,numpoint
    89     do k=1,nspec
    90       rhosum(k)=0
    91     end do
    9288    if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
    9389         (itime.le.ireleaseend(i))) then
     
    186182               xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
    187183                    *timecorrect(k)/average_timecorrect
    188               if (SCAVDEP) then ! if there is no scavenging in wetdepo it will be set to 0
     184              if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
    189185!              if ( henry(k).gt.0 .or. &
    190186!                   crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
    191187!                   ccn_aero(k).gt.0. .or. in_aero(k).gt.0. )  then
    192188                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   
    199189               endif
    200190  ! Assign certain properties to particle
     
    389379              do k=1,nspec
    390380                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
    397381              end do
    398382            endif
    399 
    400383
    401384            numpart=max(numpart,ipart)
     
    40739034      minpart=ipart+1
    408391      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
    418       endif
    419392      endif ! j=1,numrel
    420393  end do
  • src/timemanager.f90

    r462f74b r54cbd6c  
    115115  real :: xold,yold,zold,xmassfract
    116116  real, parameter :: e_inv = 1.0/exp(1.0)
     117  logical :: firstdepocalc
    117118  !double precision xm(maxspec,maxpointspec_act),
    118119  !    +                 xm_depw(maxspec,maxpointspec_act),
     
    171172           write (*,*) 'timemanager> call wetdepo'
    172173        endif     
    173          call wetdepo(itime,lsynctime,loutnext)
     174         call wetdepo(itime,lsynctime,loutnext,.false.)
    174175    endif
    175176
     
    541542        zold=ztra1(j)
    542543
     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
    543587  ! Integrate Lagevin equation for lsynctime seconds
    544588  !*************************************************
    545589
    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))
    549598
    550599  ! Calculate the gross fluxes across layer interfaces
  • src/wetdepo.f90

    r462f74b r54cbd6c  
    2020!**********************************************************************
    2121
    22 subroutine wetdepo(itime,ltsample,loutnext)
    23 !                  i      i        i
     22subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
     23!                  i      i        i            i
    2424!*****************************************************************************
    2525!                                                                            *
     
    3030! This fraction is parameterized from total cloud cover and rates of large   *
    3131! 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                                     *
    3234!                                                                            *
    3335!    Author: A. Stohl                                                        *
     
    9294  integer :: blc_count, inc_count
    9395  real    :: Si_dummy, wetscav_dummy
    94   logical :: readclouds_this_nest
     96  logical :: readclouds_this_nest,forreceptor
    9597
    9698
     
    172174           memtime(1),memtime(2),interp_time,lsp,convp,cc)
    173175    endif
    174 
    175 !  If total precipitation is less than 0.01 mm/h - no scavenging occurs
    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 
    190176
    191177! get the level were the actual particle is in
     
    413399        kp=1
    414400      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
    417404!   depostatistic
    418405!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    419406!   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
    423420!   Correct deposited mass to the last time step when radioactive decay of
    424421!   gridded deposited mass was calculated
     
    426423        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
    427424      endif
    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 
    445425
    446426    end do !all species
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG