Changeset c9cf570 in flexpart.git


Ignore:
Timestamp:
Nov 18, 2016, 11:29:23 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:
7062fab
Parents:
deaac29
Message:

first version with seperated get_wetscav

Location:
src
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • src/timemanager.f90

    rdeaac29 rc9cf570  
    105105! integer :: ksp
    106106  integer :: loutnext,loutstart,loutend
    107   integer :: ix,jy,ldeltat,itage,nage
     107  integer :: ix,jy,ldeltat,itage,nage,idummy
    108108  integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme
    109   real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact
     109  real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact,wetscav(maxspec)
    110110  ! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
    111111  ! real :: us(maxpart),vs(maxpart),ws(maxpart)
     
    114114  real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc
    115115  real :: xold,yold,zold,xmassfract
     116  real :: grfraction(3)
    116117  real, parameter :: e_inv = 1.0/exp(1.0)
    117118
     
    550551
    551552      if  (DRYBKDEP) then
    552       do ks=1,nspec
     553       do ks=1,nspec
    553554         if  ((xscav_frac1(j,ks).lt.0)) then
    554555            call advance_rec(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
    555 
    556             if (decay(ks).gt.0.) then             ! radioactive decay
    557                 decfact=exp(-real(abs(lsynctime))*decay(ks))
    558             else
    559                 decfact=1.
    560             endif
    561556            if (DRYDEPSPEC(ks)) then        ! dry deposition
    562557               xscav_frac1(j,ks)=prob_rec(ks)
     
    565560                xscav_frac1(j,ks)=0.
    566561             endif
    567 !         write (*,*) 'xscav: ',j,ks,xscav_frac1(j,ks)
    568562         endif
    569        enddo
     563        enddo
    570564       endif
    571565
    572 !      if (WETBKDEP) then
    573 !      firstdepocalc=.false.
    574 !      do ks=1,nspec
    575 !         if ((xscav_frac1(j,ks).lt.0) &
    576 !                .and.firstdepocalc.eqv..false.) then
    577 !            ! Backward wetdeposition and first timestep after release
    578 !            call wetdepo(itime,lsynctime,loutnext,.true.)
    579 !            firstdepocalc=.true.
    580 !         endif
    581 !      enddo
    582 !      endif
     566       if (WETBKDEP) then
     567       do ks=1,nspec
     568         if  ((xscav_frac1(j,ks).lt.0)) then
     569            call get_wetscav(itime,lsynctime,loutnext,j,grfraction,idummy,idummy,wetscav)
     570            if (wetscav(ks).gt.0) then
     571                xscav_frac1(j,ks)=wetscav(ks)*(zpoint2(j)-zpoint1(j))
     572            else
     573                xmass1(j,ks)=0.
     574                xscav_frac1(j,ks)=0.
     575            endif
     576         endif
     577        enddo
     578       endif
    583579
    584580  ! Integrate Lagevin equation for lsynctime seconds
  • src/wetdepo.f90

    rdeaac29 rc9cf570  
    4242!                                                                            *
    4343! 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        *
    4744! ix,jy              indices of output grid cell for each particle           *
    4845! itime [s]          actual simulation time [s]                              *
    4946! jpart              particle index                                          *
    5047! 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)  *
    5348! loutnext [s]       time for which gridded deposition is next output        *
    5449! loutstep [s]       interval at which gridded deposition is output          *
    55 ! lsp [mm/h]         large scale precipitation rate                          *
    5650! ltsample [s]       interval over which mass is deposited                   *
    57 ! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
    5851! wetdeposit         mass that is wet deposited                              *
    5952! wetgrid            accumulated deposited mass on output grid               *
     
    7164
    7265  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    73   integer :: ngrid,itage,nage,hz,il,interp_time, n
    74   integer(kind=1) :: clouds_v
     66  integer :: itage,nage,hz,il,interp_time, n
    7567  integer :: ks, kp
    76 !  integer :: n1,n2, icbot,ictop, indcloud !TEST
    77   real :: S_i, act_temp, cl, cle ! in cloud scavenging
    78   real :: clouds_h ! cloud height for the specific grid point
    79   real :: xtn,ytn,lsp,convp,cc,grfraction(3),prec(3),wetscav,totprec
     68  integer :: blc_count, inc_count
     69  real :: grfraction(3),wetscav
    8070  real :: wetdeposit(maxspec),restmass
    8171  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    82 !save lfr,cfr
    8372
    84   real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
    85   real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
    86 
    87 !ZHG aerosol below-cloud scavenging removal polynomial constants for rain and snow
    88   real, parameter :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
    89   real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
    90   real :: frac_act, liq_frac, dquer_m
    91 
    92   integer :: blc_count, inc_count, firsttimerem
    93   real    :: Si_dummy, wetscav_dummy
    94   logical :: readclouds_this_nest
    95 
    96   firsttimerem=0
    9773! Compute interval since radioactive decay of deposited mass was computed
    9874!************************************************************************
     
    11995    endif
    12096
    121 ! Determine age class of the particle
    122     itage=abs(itra1(jpart)-itramem(jpart))
    123     do nage=1,nageclass
    124       if (itage.lt.lage(nage)) goto 33
    125     end do
    126 33  continue
     97! Determine age class of the particle - nage is used for the kernel
     98!******************************************************************
     99     itage=abs(itra1(jpart)-itramem(jpart))
     100     do nage=1,nageclass
     101       if (itage.lt.lage(nage)) goto 33
     102     end do
     103 33  continue
    127104
    128 
    129 ! Determine which nesting level to be used
    130 !*****************************************
    131 
    132     ngrid=0
    133     do j=numbnests,1,-1
    134       if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
    135            (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
    136         ngrid=j
    137         goto 23
    138       endif
    139     end do
    140 23  continue
    141 
    142 
    143 ! Determine nested grid coordinates
    144 !**********************************
    145     readclouds_this_nest=.false.
    146 
    147     if (ngrid.gt.0) then
    148       xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
    149       ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
    150       ix=int(xtn)
    151       jy=int(ytn)
    152       if (readclouds_nest(ngrid)) readclouds_this_nest=.true.
    153     else
    154       ix=int(xtra1(jpart))
    155       jy=int(ytra1(jpart))
    156     endif
    157 
    158 
    159 ! Interpolate large scale precipitation, convective precipitation and
    160 ! total cloud cover
    161 ! Note that interpolated time refers to itime-0.5*ltsample [PS]
    162 !********************************************************************
    163     interp_time=nint(itime-0.5*ltsample)
    164 
    165     if (ngrid.eq.0) then
    166       call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
    167            1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
    168            memtime(1),memtime(2),interp_time,lsp,convp,cc)
    169     else
    170       call interpol_rain_nests(lsprecn,convprecn,tccn, &
    171            nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
    172            memtime(1),memtime(2),interp_time,lsp,convp,cc)
    173     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 (WETBKDEP) 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 
    190 
    191 ! get the level were the actual particle is in
    192     do il=2,nz
    193       if (height(il).gt.ztra1(jpart)) then
    194         hz=il-1
    195 !        goto 26
    196         exit
    197       endif
    198     end do
    199 !26  continue
    200 
    201     n=memind(2)
    202     if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    203          n=memind(1)
    204 
    205     if (ngrid.eq.0) then
    206       clouds_v=clouds(ix,jy,hz,n)
    207       clouds_h=cloudsh(ix,jy,n)
    208     else
    209       clouds_v=cloudsn(ix,jy,hz,n,ngrid)
    210       clouds_h=cloudshn(ix,jy,n,ngrid)
    211     endif
    212 
    213 ! if there is no precipitation or the particle is above the clouds no
    214 ! scavenging is done
    215 
    216     if (clouds_v.le.1) goto 20
    217 
    218 ! 1) Parameterization of the the area fraction of the grid cell where the
    219 !    precipitation occurs: the absolute limit is the total cloud cover, but
    220 !    for low precipitation rates, an even smaller fraction of the grid cell
    221 !    is used. Large scale precipitation occurs over larger areas than
    222 !    convective precipitation.
    223 !**************************************************************************
    224 
    225     if (lsp.gt.20.) then
    226       i=5
    227     else if (lsp.gt.8.) then
    228       i=4
    229     else if (lsp.gt.3.) then
    230       i=3
    231     else if (lsp.gt.1.) then
    232       i=2
    233     else
    234       i=1
    235     endif
    236 
    237     if (convp.gt.20.) then
    238       j=5
    239     else if (convp.gt.8.) then
    240       j=4
    241     else if (convp.gt.3.) then
    242       j=3
    243     else if (convp.gt.1.) then
    244       j=2
    245     else
    246       j=1
    247     endif
    248 
    249 
    250 !ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
    251 ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
    252 ! for now they are treated the same
    253     grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
    254     grfraction(2)=max(0.05,cc*(lfr(i)))
    255     grfraction(3)=max(0.05,cc*(cfr(j)))
    256 
    257 
    258 ! 2) Computation of precipitation rate in sub-grid cell
    259 !******************************************************
    260     prec(1)=(lsp+convp)/grfraction(1)
    261     prec(2)=(lsp)/grfraction(2)
    262     prec(3)=(convp)/grfraction(3)
    263 
    264 
    265 ! 3) Computation of scavenging coefficients for all species
    266 !    Computation of wet deposition
    267 !**********************************************************
    268 
    269     do ks=1,nspec      ! loop over species
    270       wetdeposit(ks)=0.
    271       wetscav=0.   
    272 
    273       if (ngrid.gt.0) then
    274         act_temp=ttn(ix,jy,hz,n,ngrid)
    275       else
    276         act_temp=tt(ix,jy,hz,n)
    277       endif
    278 
    279 
    280 !***********************
    281 ! BELOW CLOUD SCAVENGING
    282 !*********************** 
    283       if (clouds_v.ge.4) then !below cloud
    284 
    285 ! For gas: if positive below-cloud parameters (A or B), and dquer<=0
    286 !******************************************************************
    287         if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
    288           !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
    289           blc_count=blc_count+1
    290           wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
    291 
    292 ! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
    293 !*********************************************************************************
    294         else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
    295           blc_count=blc_count+1
    296 
    297 !NIK 17.02.2015
    298 ! For the calculation here particle size needs to be in meter and not um as dquer is
    299 ! changed in readreleases
    300 ! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
    301           dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
    302 
    303 ! Rain:
    304           if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.)  then
    305 
    306 ! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
    307 ! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
    308             wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
    309                  & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
    310                  &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
    311 
    312 ! Snow:
    313           elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.)  then
    314 ! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
    315 ! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
    316             wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
    317                  &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
    318                  &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
    319 
    320           endif
    321          
    322 !             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
    323 
    324         endif ! gas or particle
    325 !      endif ! positive below-cloud scavenging parameters given in Species file
    326       endif !end BELOW
    327 
    328 !********************
    329 ! IN CLOUD SCAVENGING
    330 !********************
    331       if (clouds_v.lt.4) then ! In-cloud
    332 ! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
    333 ! given in species file, or if gas and positive Henry's constant
    334         if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
    335           inc_count=inc_count+1
    336 ! if negative coefficients (turned off) set to zero for use in equation
    337           if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
    338           if (in_aero(ks).lt.0.) in_aero(ks)=0.
    339 
    340 !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
    341 ! nested fields
    342           if (ngrid.gt.0.and.readclouds_this_nest) then
    343             cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
    344           else if (ngrid.eq.0.and.readclouds) then
    345             cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
    346           else                                  !parameterize cloudwater m2/m3
    347 !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
    348             cl=1.6E-6*prec(1)**0.36
    349           endif
    350 
    351 !ZHG: Calculate the partition between liquid and water phase water.
    352           if (act_temp .le. 253.) then
    353             liq_frac=0
    354           else if (act_temp .ge. 273.) then
    355             liq_frac=1
    356           else
    357             liq_frac =((act_temp-273.)/(273.-253.))**2.
    358           end if
    359 ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
    360           frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
    361 
    362 !ZHG Use the activated fraction and the liqid water to calculate the washout
    363 
    364 ! AEROSOL
    365 !********
    366           if (dquer(ks).gt.0.) then
    367             S_i= frac_act/cl
    368 
    369 ! GAS
    370 !****
    371           else
    372 
    373             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    374 !REPLACE to switch old/ new scheme
    375           ! S_i=frac_act/cle
    376             S_i=1/cle
    377           endif ! gas or particle
    378 
    379 ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
    380 !OLD
    381           if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
    382             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
    383           else
    384             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
    385           endif
    386 
    387 
    388         endif ! positive in-cloud scavenging parameters given in Species file
    389       endif !incloud
    390 !END ZHG TEST
    391105
    392106!**************************************************
    393107! CALCULATE DEPOSITION
    394108!**************************************************
    395 !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
    396 !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
     109      call get_wetscav(itime,ltsample,loutnext,jpart,grfraction,inc_count,blc_count,wetscav)
    397110
    398111      if (wetscav.gt.0.) then
    399112        wetdeposit(ks)=xmass1(jpart,ks)* &
    400113             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
    401 !write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &
    402 !             (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v
    403 
    404 
    405114      else ! if no scavenging
    406115        wetdeposit(ks)=0.
    407116      endif
    408 
     117 
    409118      restmass = xmass1(jpart,ks)-wetdeposit(ks)
    410119      if (ioutputforeachrelease.eq.1) then
     
    427136      endif
    428137
    429       if (WETBKDEP) 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                  if (xmass1(jpart,ks).eq.0) then
    439                      firsttimerem=firsttimerem+1
    440                 endif
    441                  !xscav_frac1(jpart,ks)=wetdeposit(ks)/(xmass1(jpart,ks)+wetdeposit(ks))
    442                  xscav_frac1(jpart,ks)=wetscav*(zpoint2(jpart)-zpoint1(jpart))
    443 !                write (*,*) 'paricle kept: ',jpart,ks,wetdeposit(ks),xscav_frac1(jpart,ks),xmass(jpart,ks)
    444              endif
    445           endif
    446       endif
    447 
    448 
    449     end do !all species
    450 
    451138! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
    452139! Add the wet deposition to accumulated amount on output grid and nested output grid
     
    467154  tot_inc_count=tot_inc_count+inc_count
    468155
    469 !  write (*,*) 'First time removed: ',firsttimerem
    470 
    471156end subroutine wetdepo
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG