Changes in / [d9f0585:6985a98] in flexpart.git


Ignore:
Files:
3 added
19 edited

Legend:

Unmodified
Added
Removed

  • options/SPECIES/spec_overview

    • Property mode changed from 100644 to 100755

  • src/advance.f90

    ra652cd5 r94735e2  
    117117  real :: usigold,vsigold,wsigold,r,rs
    118118  real :: uold,vold,wold,vdepo(maxspec)
     119  real :: h1(2)
    119120  !real uprof(nzmax),vprof(nzmax),wprof(nzmax)
    120121  !real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
    … …  
    223224
    224225
     226 ! Determine the lower left corner and its distance to the current position
     227  !*************************************************************************
     228
     229  ddx=xt-real(ix)
     230  ddy=yt-real(jy)
     231  rddx=1.-ddx
     232  rddy=1.-ddy
     233  p1=rddx*rddy
     234  p2=ddx*rddy
     235  p3=rddx*ddy
     236  p4=ddx*ddy
     237
     238 ! Calculate variables for time interpolation
     239  !*******************************************
     240
     241  dt1=real(itime-memtime(1))
     242  dt2=real(memtime(2)-itime)
     243  dtt=1./(dt1+dt2)
     244
    225245  ! Compute maximum mixing height around particle position
    226246  !*******************************************************
    … …  
    229249  if (ngrid.le.0) then
    230250    do k=1,2
    231       mind=memind(k) ! eso: compatibility with 3-field version
    232       do j=jy,jyp
    233         do i=ix,ixp
    234           if (hmix(i,j,1,mind).gt.h) h=hmix(i,j,1,mind)
    235         end do
    236       end do
     251       mind=memind(k) ! eso: compatibility with 3-field version
     252       if (interpolhmix) then
     253             h1(k)=p1*hmix(ix ,jy ,1,mind) &
     254                 + p2*hmix(ixp,jy ,1,mind) &
     255                 + p3*hmix(ix ,jyp,1,mind) &
     256                 + p4*hmix(ixp,jyp,1,mind)
     257        else
     258          do j=jy,jyp
     259            do i=ix,ixp
     260               if (hmix(i,j,1,mind).gt.h) h=hmix(i,j,1,mind)
     261            end do
     262          end do
     263        endif
    237264    end do
    238265    tropop=tropopause(nix,njy,1,1)
    … …  
    249276  endif
    250277
     278  if (interpolhmix) h=(h1(1)*dt2+h1(2)*dt1)*dtt
    251279  zeta=zt/h
    252280
    … …  
    446474        endif
    447475
     476        if (turboff) then
     477!sec switch off turbulence
     478          up=0.0
     479          vp=0.0
     480          wp=0.0
     481          delz=0.
     482        endif
     483
    448484  !****************************************************
    449485  ! Compute turbulent vertical displacement of particle
    … …  
    647683  endif
    648684
     685  if (turboff) then
     686!sec switch off turbulence
     687    ux=0.0
     688    vy=0.0
     689    wp=0.0
     690  endif
    649691
    650692  ! If particle represents only a single species, add gravitational settling

  • src/com_mod.f90

    rc8fc724 r5888298  
    577577  real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
    578578  !real outheight(maxzgrid),outheighthalf(maxzgrid)
     579
    579580  logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
    580581       & OHREA,ASSSPEC
     582  logical :: DRYBKDEP,WETBKDEP
    581583
    582584  ! numxgrid,numygrid       number of grid points in x,y-direction
    … …  
    598600  ! OHREA                   .true., if OH reaction is switched on
    599601  ! ASSSPEC                 .true., if there are two species asscoiated
     602  ! DRYBKDEP,WETBKDEP        .true., for bkwd runs, where mass deposited and source regions is calculated - either for dry or for wet deposition
    600603  !                    (i.e. transfer of mass between these two occurs
    601604
    … …  
    668671  real, allocatable, dimension(:) :: ztra1
    669672  real, allocatable, dimension(:,:) :: xmass1
     673  real, allocatable, dimension(:,:) :: xscav_frac1
    670674
    671675  ! eso: Moved from timemanager
    … …  
    688692  ! xtra1,ytra1,ztra1       spatial positions of the particles
    689693  ! xmass1 [kg]             particle masses
    690  
     694  ! xscav_frac1             fraction of particle masse which has been scavenged at receptor
     695 
    691696
    692697
    … …  
    750755  integer :: mpi_mode=0 ! .gt. 0 if running MPI version
    751756  logical :: lroot=.true. ! true if serial version, or if MPI .and. root process
     757 
     758  logical :: usekernel=.false.    ! true if the output kernel shall be switched on
     759  logical :: interpolhmix=.false. ! true if the hmix shall be interpolated
     760  logical :: turboff=.false.       ! true if the turbulence shall be switched off
     761 
    752762 
    753763contains

  • src/conccalc.f90

    r4c64400 r9287c01  
    2121
    2222subroutine conccalc(itime,weight)
    23   !                 i     i
     23  !                      i     i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    5959  real :: xl,yl,wx,wy,w
    6060  real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
    61 
     61!  integer xscav_count
    6262
    6363  ! For forward simulations, make a loop over the number of species;
    … …  
    6565  ! releasepoints
    6666  !***************************************************************************
    67 
     67!  xscav_count=0
    6868  do i=1,numpart
    6969    if (itra1(i).ne.itime) goto 20
    … …  
    767633   continue
    7777
    78 
     78!  if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
     79           
    7980  ! For special runs, interpolate the air density to the particle position
    8081  !************************************************************************
    … …  
    172173      if (yl.lt.0.) jy=jy-1
    173174
    174   ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl
    175175
    176176
    … …  
    183183      if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
    184184           (xl.gt.real(numxgrid-1)-0.5).or. &
    185            (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
     185           (yl.gt.real(numygrid-1)-0.5))) then             ! no kernel, direct attribution to grid cell
    186186        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    187187             (jy.le.numygrid-1)) then
    188           do ks=1,nspec
    189             gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     188          if (DRYBKDEP.or.WETBKDEP) then
     189             do ks=1,nspec
     190               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     191                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
     192                 xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
     193             end do
     194          else
     195             do ks=1,nspec
     196               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    190197                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    191198                 xmass1(i,ks)/rhoi*weight
    192           end do
     199             end do
     200          endif
    193201        endif
    194202
    195       else                                 ! attribution via uniform kernel
     203      else                                 ! attribution via uniform kernel 
    196204
    197205        ddx=xl-real(ix)                   ! distance to left cell border
    … …  
    220228          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    221229            w=wx*wy
    222             do ks=1,nspec
    223               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     230            if (DRYBKDEP.or.WETBKDEP) then
     231               do ks=1,nspec
     232                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     233                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
     234                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
     235               end do
     236            else
     237               do ks=1,nspec
     238                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    224239                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    225240                   xmass1(i,ks)/rhoi*weight*w
    226             end do
     241               end do
     242            endif
    227243          endif
    228244
    229245          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    230246            w=wx*(1.-wy)
    231             do ks=1,nspec
    232               gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     247            if (DRYBKDEP.or.WETBKDEP) then
     248              do ks=1,nspec
     249                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     250                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
     251                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
     252               end do
     253             else
     254              do ks=1,nspec
     255                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    233256                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    234257                   xmass1(i,ks)/rhoi*weight*w
    235             end do
    236           endif
    237         endif
     258               end do
     259             endif
     260          endif
     261        endif !ix ge 0
    238262
    239263
    … …  
    241265          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    242266            w=(1.-wx)*(1.-wy)
    243             do ks=1,nspec
    244               gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     267            if (DRYBKDEP.or.WETBKDEP) then
     268               do ks=1,nspec
     269                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     270                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
     271                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
     272               end do
     273            else
     274               do ks=1,nspec
     275                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    245276                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    246277                   xmass1(i,ks)/rhoi*weight*w
    247             end do
     278               end do
     279            endif
    248280          endif
    249281
    250282          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    251283            w=(1.-wx)*wy
    252             do ks=1,nspec
    253               gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     284            if (DRYBKDEP.or.WETBKDEP) then
     285               do ks=1,nspec
     286                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     287                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
     288                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
     289               end do
     290            else
     291               do ks=1,nspec
     292                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    254293                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    255294                   xmass1(i,ks)/rhoi*weight*w
    256             end do
    257           endif
    258         endif
    259       endif
    260 
    261 
     295               end do
     296            endif
     297          endif
     298        endif !ixp ge 0
     299     endif
    262300
    263301  !************************************
    … …  
    282320        if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
    283321             (xl.gt.real(numxgridn-1)-0.5).or. &
    284              (yl.gt.real(numygridn-1)-0.5)) then             ! no kernel, direct attribution to grid cell
     322             (yl.gt.real(numygridn-1)-0.5).or.(.not.usekernel)) then             ! no kernel, direct attribution to grid cell
    285323          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
    286324               (jy.le.numygridn-1)) then
    287             do ks=1,nspec
    288               griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     325            if (DRYBKDEP.or.WETBKDEP) then
     326               do ks=1,nspec
     327                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     328                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
     329                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
     330               end do
     331            else
     332               do ks=1,nspec
     333                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    289334                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    290335                   xmass1(i,ks)/rhoi*weight
    291             end do
     336               end do
     337            endif
    292338          endif
    293339
    … …  
    319365            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    320366              w=wx*wy
    321               do ks=1,nspec
    322                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     367              if (DRYBKDEP.or.WETBKDEP) then
     368                 do ks=1,nspec
     369                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     370                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
     371                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
     372                 end do
     373              else
     374                do ks=1,nspec
     375                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    323376                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    324377                     xmass1(i,ks)/rhoi*weight*w
    325               end do
     378                 end do
     379              endif
    326380            endif
    327381
    328382            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    329383              w=wx*(1.-wy)
    330               do ks=1,nspec
    331                 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     384              if (DRYBKDEP.or.WETBKDEP) then
     385                 do ks=1,nspec
     386                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     387                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
     388                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
     389                 end do
     390              else
     391                 do ks=1,nspec
     392                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    332393                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    333394                     xmass1(i,ks)/rhoi*weight*w
    334               end do
     395                 end do
     396              endif
    335397            endif
    336398          endif
    … …  
    340402            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    341403              w=(1.-wx)*(1.-wy)
    342               do ks=1,nspec
    343                 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     404              if (DRYBKDEP.or.WETBKDEP) then
     405                 do ks=1,nspec
     406                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
     407                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
     408                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
     409                 end do
     410              else
     411                 do ks=1,nspec
     412                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    344413                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    345414                     xmass1(i,ks)/rhoi*weight*w
    346               end do
     415                 end do
     416              endif
    347417            endif
    348418
    349419            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    350420              w=(1.-wx)*wy
    351               do ks=1,nspec
    352                 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     421              if (DRYBKDEP.or.WETBKDEP) then
     422                 do ks=1,nspec
     423                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
     424                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
     425                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
     426                 end do
     427              else
     428                 do ks=1,nspec
     429                    griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    353430                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    354431                     xmass1(i,ks)/rhoi*weight*w
    355               end do
     432                 end do
     433              endif
    356434            endif
    357435          endif
    358436        endif
    359 
    360437      endif
    361438    endif
    36243920  continue
    363440  end do
     441!  write(*,*) 'xscav count:',xscav_count
    364442
    365443  !***********************************************************************

  • src/concoutput.f90

    r16b61a5 r16b61a5  
    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// &
     248
     249    if (DRYBKDEP.or.WETBKDEP) then !scavdep output
     250      if (DRYBKDEP) &
     251      open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
     252           atime//'_'//anspec,form='unformatted')
     253      if (WETBKDEP) &
     254      open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
     255           atime//'_'//anspec,form='unformatted')
     256      write(unitoutgrid) itime
     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// &
    251261             atime//'_'//anspec,form='unformatted')
    252       else
    253         open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
     262        else
     263          open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
    254264             atime//'_'//anspec,form='unformatted')
     265        endif
     266        write(unitoutgrid) itime
    255267      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// &
     268      if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
     269        open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
    261270           atime//'_'//anspec,form='unformatted')
    262 
    263       write(unitoutgridppt) itime
    264     endif
     271        write(unitoutgridppt) itime
     272      endif
     273    endif ! if deposition output
    265274
    266275    do kp=1,maxpointspec_act
    … …  
    342351! Concentration output
    343352!*********************
    344         if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     353        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
    345354
    346355! Wet deposition

  • src/drydepokernel.f90

    r4c64400 r1c0d5e6  
    102102   if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
    103103
    104    if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
     104    if (.not.usekernel) then
     105       drygridunc(ix,jy,ks,kp,nunc,nage)= &
     106           drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
     107    else
     108      if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    105109        (jy.le.numygrid-1)) then
    106      w=wx*wy
    107      drygridunc(ix,jy,ks,kp,nunc,nage)= &
    108           drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
    109      continue
    110    endif
     110        w=wx*wy
     111        drygridunc(ix,jy,ks,kp,nunc,nage)= &
     112           drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
     113     endif
    111114
    112   if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
     115    if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
    113116       (jyp.le.numygrid-1)) then
    114117    w=(1.-wx)*(1.-wy)
    115118      drygridunc(ixp,jyp,ks,kp,nunc,nage)= &
    116119           drygridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    117   endif
     120    endif
    118121
    119   if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
     122    if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
    120123       (jy.le.numygrid-1)) then
    121     w=(1.-wx)*wy
     124      w=(1.-wx)*wy
    122125      drygridunc(ixp,jy,ks,kp,nunc,nage)= &
    123126           drygridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
    124   endif
     127    endif
    125128
    126   if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
     129    if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
    127130       (jyp.le.numygrid-1)) then
    128     w=wx*(1.-wy)
     131      w=wx*(1.-wy)
    129132      drygridunc(ix,jyp,ks,kp,nunc,nage)= &
    130133           drygridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    131   endif
     134    endif
    132135
    133   endif
     136    endif ! kernel
     137    endif ! deposit>0
    134138
    135139  end do

  • src/ecmwf_mod.f90

    r62e65c7 r842074e  
    4343
    4444!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new
    45   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
    46   integer,parameter :: nxshift=359
     45   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
     46!  integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 0.5
     47   integer,parameter :: nxshift=359
     48!  integer,parameter :: nxshift=718
    4749!  integer,parameter :: nxshift=0
    4850

  • src/getfields.f90

    r8a65cb0 rd1a8707  
    13813840  indmin=indj
    139139
     140   if (WETBKDEP) then
     141        call writeprecip(itime,memind(1))
     142   endif
     143
    140144  else
    141145
    … …  
    16917360  indmin=indj
    170174
     175   if (WETBKDEP) then
     176        call writeprecip(itime,memind(1))
     177   endif
     178
    171179  endif
    172180

  • src/interpol_rain.f90

    r8a65cb0 rc7e771d  
    2121
    2222subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
    23        ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
     23       ny,iwftouse,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
    2424  !                          i   i   i    i    i     i   i
    2525  !i    i    i  i    i     i      i      i     o     o     o
    … …  
    7777  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
    7878  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
     79  integer :: iwftouse
    7980
    8081
    … …  
    113114  !***********************
    114115
    115   do m=1,2
    116     indexh=memind(m)
     116!  do m=1,2
     117    indexh=iwftouse
    117118
    118 
    119     y1(m)=p1*yy1(ix ,jy ,level,indexh) &
     119    y1(1)=p1*yy1(ix ,jy ,level,indexh) &
    120120         + p2*yy1(ixp,jy ,level,indexh) &
    121121         + p3*yy1(ix ,jyp,level,indexh) &
    122122         + p4*yy1(ixp,jyp,level,indexh)
    123     y2(m)=p1*yy2(ix ,jy ,level,indexh) &
     123    y2(1)=p1*yy2(ix ,jy ,level,indexh) &
    124124         + p2*yy2(ixp,jy ,level,indexh) &
    125125         + p3*yy2(ix ,jyp,level,indexh) &
    126126         + p4*yy2(ixp,jyp,level,indexh)
    127     y3(m)=p1*yy3(ix ,jy ,level,indexh) &
     127    y3(1)=p1*yy3(ix ,jy ,level,indexh) &
    128128         + p2*yy3(ixp,jy ,level,indexh) &
    129129         + p3*yy3(ix ,jyp,level,indexh) &
    130130         + p4*yy3(ixp,jyp,level,indexh)
    131   end do
     131!  end do
    132132
    133133
    134134  !************************************
    135   ! 2.) Temporal interpolation (linear)
     135  ! 2.) Temporal interpolation (linear) - skip to be consistent with clouds
    136136  !************************************
    137137
    138   dt1=real(itime-itime1)
    139   dt2=real(itime2-itime)
    140   dt=dt1+dt2
     138!  dt1=real(itime-itime1)
     139!  dt2=real(itime2-itime)
     140!  dt=dt1+dt2
    141141
    142   yint1=(y1(1)*dt2+y1(2)*dt1)/dt
    143   yint2=(y2(1)*dt2+y2(2)*dt1)/dt
    144   yint3=(y3(1)*dt2+y3(2)*dt1)/dt
     142!  yint1=(y1(1)*dt2+y1(2)*dt1)/dt
     143!  yint2=(y2(1)*dt2+y2(2)*dt1)/dt
     144!  yint3=(y3(1)*dt2+y3(2)*dt1)/dt
    145145
     146   yint1=y1(1)
     147   yint2=y2(1)
     148   yint3=y3(1)
    146149
    147150end subroutine interpol_rain

  • src/interpol_vdep.f90

    re200b7a r5844866  
    5454
    5555  ! a) Bilinear horizontal interpolation
    56 
     56! write(*,*) 'interpol: ',dt1,dt2,dtt,lsynctime,ix,jy
    5757  do m=1,2
    5858    indexh=memind(m)

  • src/makefile

    r4c64400 rd1a8707  
    6868LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
    6969
    70 FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) # -march=native
     70FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER)  #-Warray-bounds -fcheck=all # -march=native
    7171
    7272DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)  #  -ffpe-trap=invalid,overflow,denormal,underflow,zero  -Warray-bounds -fcheck=all
    … …  
    144144advance.o               initialize.o            \
    145145writeheader.o           writeheader_txt.o       \
     146writeprecip.o \
    146147writeheader_surf.o      assignland.o\
    147148part0.o                 gethourlyOH.o\
    … …  
    152153erf.o                   readavailable.o \
    153154ew.o                    readreleases.o  \
    154 readdepo.o \
     155readdepo.o              get_vdep_prob.o   \
     156get_wetscav.o   \
    155157psim.o                  outgrid_init.o  \
    156158outgrid_init_nest.o     \
    … …  
    271273
    272274## DEPENDENCIES
     275get_vdep_prob.o: cmapf_mod.o com_mod.o hanna_mod.o interpol_mod.o par_mod.o \
     276        point_mod.o random_mod.o
    273277advance.o: cmapf_mod.o com_mod.o hanna_mod.o interpol_mod.o par_mod.o \
    274278        point_mod.o random_mod.o
    … …  
    424428verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
    425429verttransform_nests.o: com_mod.o par_mod.o
     430get_wetscav.o: com_mod.o par_mod.o point_mod.o
    426431wetdepo.o: com_mod.o par_mod.o point_mod.o
    427432wetdepokernel.o: com_mod.o par_mod.o unc_mod.o
    428433wetdepokernel_nest.o: com_mod.o par_mod.o unc_mod.o
    429434writeheader.o: com_mod.o outg_mod.o par_mod.o point_mod.o
     435writeprecip.o: com_mod.o par_mod.o point_mod.o
    430436writeheader_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o
    431437writeheader_nest_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o

  • src/par_mod.f90

    r4c64400 rd1a8707  
    194194  integer,parameter :: maxpart=10000000
    195195  integer,parameter :: maxspec=4
     196
    196197  real,parameter :: minmass=0.0001
    197198
    … …  
    258259  integer,parameter :: unitlsm=1, unitsurfdata=1, unitland=1, unitwesely=1
    259260  integer,parameter :: unitOH=1
    260   integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
     261  integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95, unitprecip=101
    261262  integer,parameter :: unitboundcond=89
    262263  integer,parameter :: unittmp=101

  • src/readcommand.f90

    r4c64400 r8ee24a5  
    298298  !Af          1 = mass units
    299299  !Af          2 = mass mixing ratio units
     300  !            3 = wet deposition in outputfield
     301  !            4 = dry deposition in outputfield
    300302
    301303  if ( ldirect .eq. 1 ) then  ! FWD-Run
    … …  
    320322     endif
    321323  !Af set release-switch
    322      if (ind_receptor .eq. 1) then !mass
     324     WETBKDEP=.false.
     325     DRYBKDEP=.false.
     326     select case (ind_receptor)
     327     case (1)  !  1 .. concentration at receptor
    323328        ind_rel = 1
    324      else ! mass mix
     329     case (2)  !  2 .. mixing ratio at receptor
    325330        ind_rel = 0
    326      endif
     331     case (3)  ! 3 .. wet deposition in outputfield
     332        ind_rel = 3
     333         write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
     334         write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
     335         write(*,*) ' #### Release is performed above ground lev    #### '
     336         WETBKDEP=.true.
     337         allocate(xscav_frac1(maxpart,maxspec))
     338     case (4)  ! 4 .. dry deposition in outputfield
     339         ind_rel = 4
     340         write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE    #### '
     341         write(*,*) ' #### Releaseheight is forced to 0 - 2*href    #### '
     342         write(*,*) ' #### Release is performed above ground lev    #### '
     343         DRYBKDEP=.true.
     344         allocate(xscav_frac1(maxpart,maxspec))
     345     end select
    327346  endif
    328347
    … …  
    386405  !**********************************************************************
    387406
    388   if ((iout.lt.1).or.(iout.gt.5)) then
     407  if ((iout.lt.1).or.(iout.gt.6)) then
    389408    write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
    390     write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR       #### '
     409    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR        #### '
    391410    write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
    392411    write(*,*) ' #### FOR NETCDF OUTPUT                       #### '

  • src/readreleases.f90

    rc8fc724 rc8fc724  
    521521  endif
    522522
     523  ! If FLEXPART is run for backward deposition force zpoint
     524  !*********************************************************************
     525  if (WETBKDEP) then
     526      zpoint1(numpoint)=0.
     527      zpoint2(numpoint)=20000.
     528      kindz(numpoint)=1
     529  endif
     530  if (DRYBKDEP) then
     531      zpoint1(numpoint)=0.
     532      zpoint2(numpoint)=2.*href
     533      kindz(numpoint)=1
     534  endif
     535
     536
    523537  ! Check whether x coordinates of release point are within model domain
    524538  !*********************************************************************

  • src/releaseparticles.f90

    r8a65cb0 r75a4ded  
    176176  ! scaled. Adjust the mass per particle by the species-dependent time correction factor
    177177  ! divided by the species-average one
     178  ! for the scavenging calculation the mass needs to be multiplied with rho of the particle layer and
     179  ! divided by the sum of rho of all particles.
    178180  !*****************************************************************************
    179181            do k=1,nspec
    180182               xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
    181183                    *timecorrect(k)/average_timecorrect
    182   !            write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
     184              if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
     185!              if ( henry(k).gt.0 .or. &
     186!                   crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
     187!                   ccn_aero(k).gt.0. .or. in_aero(k).gt.0. )  then
     188                xscav_frac1(ipart,k)=-1.
     189               endif
    183190  ! Assign certain properties to particle
    184191  !**************************************
    … …  
    202209
    203210            ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
    204 
    205211  ! Interpolation of topography and density
    206212  !****************************************
    … …  
    330336  !Af ind_rel is defined in readcommand.f
    331337
    332             if (ind_rel .eq. 1) then
     338            if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
    333339
    334340  ! Interpolate the air density
    … …  
    375381            endif
    376382
    377 
    378383            numpart=max(numpart,ipart)
    379384            goto 34      ! Storage space has been found, stop searching
    380385          endif
    381         end do
     386        end do  ! i=1:numpoint
    382387        if (ipart.gt.maxpart) goto 996
    383388
    38438934      minpart=ipart+1
    385       end do
    386       endif
     390      end do ! ipart=minpart,maxpart
     391      endif ! j=1,numrel
    387392  end do
    388393

  • src/timemanager.f90

    r18adf60 rd1a8707  
    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(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)
     118
    117119  !double precision xm(maxspec,maxpointspec_act),
    118120  !    +                 xm_depw(maxspec,maxpointspec_act),
    … …  
    145147!CGZ-lifetime: set lifetime to 0
    146148 
     149  if (.not.usekernel) write(*,*) 'Not using the kernel'
     150  if (turboff) write(*,*) 'Turbulence switched off'
     151
    147152  write(*,46) float(itime)/3600,itime,numpart
    148153
    … …  
    541546        zold=ztra1(j)
    542547
     548   
     549  ! RECEPTOR: dry/wet depovel
     550  !****************************
     551  ! Before the particle is moved
     552  ! the calculation of the scavenged mass shall only be done once after release
     553  ! xscav_frac1 was initialised with a negative value
     554
     555      if  (DRYBKDEP) then
     556       do ks=1,nspec
     557         if  ((xscav_frac1(j,ks).lt.0)) then
     558            call get_vdep_prob(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
     559            if (DRYDEPSPEC(ks)) then        ! dry deposition
     560               xscav_frac1(j,ks)=prob_rec(ks)
     561             else
     562                xmass1(j,ks)=0.
     563                xscav_frac1(j,ks)=0.
     564             endif
     565         endif
     566        enddo
     567       endif
     568
     569       if (WETBKDEP) then
     570       do ks=1,nspec
     571         if  ((xscav_frac1(j,ks).lt.0)) then
     572            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
     573            if (wetscav(ks).gt.0) then
     574                xscav_frac1(j,ks)=wetscav(ks)* &
     575                       (zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
     576            else
     577                xmass1(j,ks)=0.
     578                xscav_frac1(j,ks)=0.
     579            endif
     580         endif
     581        enddo
     582       endif
     583
    543584  ! Integrate Lagevin equation for lsynctime seconds
    544585  !*************************************************
    545586
     587        if (verbosity.gt.0) then
     588           if (j.eq.1) then
     589             write (*,*) 'timemanager> call advance'
     590           endif     
     591        endif
     592     
    546593        call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
    547594             us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
    548595             cbt(j))
     596!        write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
    549597
    550598  ! Calculate the gross fluxes across layer interfaces

  • src/verttransform.f90

    r341f4b7 r6481010  
    9191
    9292  logical :: init = .true.
     93  logical :: init_w = .false.
     94  logical :: init_r = .false.
     95
    9396
    9497  !ZHG SEP 2014 tests 
    … …  
    103106  ! CHARACTER(LEN=3)  :: aspec
    104107  ! integer :: virr=0
    105   real :: tot_cloud_h
    106   real :: dbg_height(nzmax)
     108  !real :: tot_cloud_h
     109  !real :: dbg_height(nzmax)
    107110!ZHG
    108111
    … …  
    123126  if (init) then
    124127
     128
     129    if (init_r) then
     130
     131        open(333,file='heights.txt', &
     132          form='formatted')
     133        do kz=1,nuvz
     134            read(333,*) height(kz)
     135        end do
     136        close(333)
     137        write(*,*) 'height read'
     138    else
     139
     140
    125141! Search for a point with high surface pressure (i.e. not above significant topography)
    126142! Then, use this point to construct a reference z profile, to be used at all times
    … …  
    161177    end do
    162178
     179    if (init_w) then
     180        open(333,file='heights.txt', &
     181          form='formatted')
     182        do kz=1,nuvz
     183              write(333,*) height(kz)
     184        end do
     185        close(333)
     186    endif
     187
     188    endif ! init
    163189
    164190! Determine highest levels that can be within PBL
    … …  
    178204    init=.false.
    179205
    180     dbg_height = height
     206!    dbg_height = height
    181207
    182208  endif
    … …  
    598624        convp=convprec(ix,jy,1,n)
    599625        prec=lsp+convp
    600         tot_cloud_h=0
     626!        tot_cloud_h=0
    601627! Find clouds in the vertical
    602628        do kz=1, nz-1 !go from top to bottom
    … …  
    604630! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
    605631            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
    606             tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
     632!            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
    607633           
    608634!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]

  • src/wetdepo.f90

    rc8fc724 rc8fc724  
    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,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(selected_int_kind(16)), dimension(nspec) :: 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
    83 
    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(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count
    93   real    :: Si_dummy, wetscav_dummy
    94   logical :: readclouds_this_nest
    95 
    9672
    9773! Compute interval since radioactive decay of deposited mass was computed
    … …  
    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
    127 
    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     if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    177 
    178 ! get the level were the actual particle is in
    179     do il=2,nz
    180       if (height(il).gt.ztra1(jpart)) then
    181         hz=il-1
    182 !        goto 26
    183         exit
    184       endif
    185     end do
    186 !26  continue
    187 
    188     n=memind(2)
    189     if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    190          n=memind(1)
    191 
    192     if (ngrid.eq.0) then
    193       clouds_v=clouds(ix,jy,hz,n)
    194       clouds_h=cloudsh(ix,jy,n)
    195     else
    196       clouds_v=cloudsn(ix,jy,hz,n,ngrid)
    197       clouds_h=cloudshn(ix,jy,n,ngrid)
    198     endif
    199 
    200 ! if there is no precipitation or the particle is above the clouds no
    201 ! scavenging is done
    202 
    203     if (clouds_v.le.1) goto 20
    204 
    205 ! 1) Parameterization of the the area fraction of the grid cell where the
    206 !    precipitation occurs: the absolute limit is the total cloud cover, but
    207 !    for low precipitation rates, an even smaller fraction of the grid cell
    208 !    is used. Large scale precipitation occurs over larger areas than
    209 !    convective precipitation.
    210 !**************************************************************************
    211 
    212     if (lsp.gt.20.) then
    213       i=5
    214     else if (lsp.gt.8.) then
    215       i=4
    216     else if (lsp.gt.3.) then
    217       i=3
    218     else if (lsp.gt.1.) then
    219       i=2
    220     else
    221       i=1
    222     endif
    223 
    224     if (convp.gt.20.) then
    225       j=5
    226     else if (convp.gt.8.) then
    227       j=4
    228     else if (convp.gt.3.) then
    229       j=3
    230     else if (convp.gt.1.) then
    231       j=2
    232     else
    233       j=1
    234     endif
    235 
    236 
    237 !ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
    238 ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
    239 ! for now they are treated the same
    240     grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
    241     grfraction(2)=max(0.05,cc*(lfr(i)))
    242     grfraction(3)=max(0.05,cc*(cfr(j)))
    243 
    244 
    245 ! 2) Computation of precipitation rate in sub-grid cell
    246 !******************************************************
    247     prec(1)=(lsp+convp)/grfraction(1)
    248     prec(2)=(lsp)/grfraction(2)
    249     prec(3)=(convp)/grfraction(3)
    250 
    251 
    252 ! 3) Computation of scavenging coefficients for all species
    253 !    Computation of wet deposition
    254 !**********************************************************
     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
    255104
    256105    do ks=1,nspec      ! loop over species
    257       wetdeposit(ks)=0.
    258       wetscav=0.
    259106
    260 ! Cycle loop if wet deposition for the species is off
    261       if (WETDEPSPEC(ks).eqv..false.) cycle
    262 
    263       if (ngrid.gt.0) then
    264         act_temp=ttn(ix,jy,hz,n,ngrid)
    265       else
    266         act_temp=tt(ix,jy,hz,n)
    267       endif
    268 
    269 
    270 !***********************
    271 ! BELOW CLOUD SCAVENGING
    272 !*********************** 
    273       if (clouds_v.ge.4) then !below cloud
    274 
    275 ! For gas: if positive below-cloud parameters (A or B), and dquer<=0
    276 !******************************************************************
    277         if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
    278           !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
    279           blc_count(ks)=blc_count(ks)+1
    280           wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
    281 
    282 ! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
    283 !*********************************************************************************
    284         else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
    285           blc_count(ks)=blc_count(ks)+1
    286 
    287 !NIK 17.02.2015
    288 ! For the calculation here particle size needs to be in meter and not um as dquer is
    289 ! changed in readreleases
    290 ! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
    291           dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
    292 
    293 ! Rain:
    294           if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.)  then
    295 
    296 ! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
    297 ! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
    298             wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
    299                  & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
    300                  &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
    301 
    302 ! Snow:
    303           elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.)  then
    304 ! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
    305 ! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
    306             wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
    307                  &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
    308                  &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
    309 
    310           endif
    311          
    312 !             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
    313 
    314         endif ! gas or particle
    315 !      endif ! positive below-cloud scavenging parameters given in Species file
    316       endif !end BELOW
    317 
    318 !********************
    319 ! IN CLOUD SCAVENGING
    320 !********************
    321       if (clouds_v.lt.4) then ! In-cloud
    322 ! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
    323 ! given in species file, or if gas and positive Henry's constant
    324         if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
    325           inc_count(ks)=inc_count(ks)+1
    326 ! if negative coefficients (turned off) set to zero for use in equation
    327           if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
    328           if (in_aero(ks).lt.0.) in_aero(ks)=0.
    329 
    330 !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
    331 ! nested fields
    332           if (ngrid.gt.0.and.readclouds_this_nest) then
    333             cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
    334           else if (ngrid.eq.0.and.readclouds) then
    335             cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
    336           else                                  !parameterize cloudwater m2/m3
    337 !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
    338             cl=1.6E-6*prec(1)**0.36
    339           endif
    340 
    341 !ZHG: Calculate the partition between liquid and water phase water.
    342           if (act_temp .le. 253.) then
    343             liq_frac=0
    344           else if (act_temp .ge. 273.) then
    345             liq_frac=1
    346           else
    347             liq_frac =((act_temp-273.)/(273.-253.))**2.
    348           end if
    349 ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
    350           frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
    351 
    352 !ZHG Use the activated fraction and the liqid water to calculate the washout
    353 
    354 ! AEROSOL
    355 !********
    356           if (dquer(ks).gt.0.) then
    357             S_i= frac_act/cl
    358 
    359 ! GAS
    360 !****
    361           else
    362 
    363             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    364 !REPLACE to switch old/ new scheme
    365           ! S_i=frac_act/cle
    366             S_i=1/cle
    367           endif ! gas or particle
    368 
    369 ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
    370 !OLD
    371           if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
    372             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
    373           else
    374             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
    375           endif
    376 
    377 
    378         endif ! positive in-cloud scavenging parameters given in Species file
    379       endif !incloud
    380 !END ZHG TEST
     107      if (WETDEPSPEC(ks).eqv..false.) cycle
    381108
    382109!**************************************************
    383110! CALCULATE DEPOSITION
    384111!**************************************************
    385 !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
    386 !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
     112!       wetscav=0.
     113       
     114!        write(*,*) ks,dquer(ks), crain_aero(ks),csnow_aero(ks)
     115!       if (((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) &
     116!          .or. &
     117!          ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.).or. &
     118!            (ccn_aero(ks).gt0) .or. (in_aero(ks).gt.0) .or. (henry(ks).gt.0)))  then
     119
     120      call get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc_count,wetscav)
     121     
    387122
    388123      if (wetscav.gt.0.) then
    389124        wetdeposit(ks)=xmass1(jpart,ks)* &
    390125             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
    391 !write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &
    392 !             (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v
    393 
    394 
    395126      else ! if no scavenging
    396127        wetdeposit(ks)=0.
    397128      endif
    398 
     129 
    399130      restmass = xmass1(jpart,ks)-wetdeposit(ks)
    400131      if (ioutputforeachrelease.eq.1) then
    … …  
    417148      endif
    418149
    419 
    420     end do !all species
     150!    endif ! no deposition
     151    end do ! loop over species
    421152
    422153! Sabine Eckhardt, June 2008 create deposition runs only for forward runs

  • src/wetdepokernel.f90

    r4c64400 r1c0d5e6  
    9595  do ks=1,nspec
    9696
    97   if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
     97  if (.not.usekernel) then
     98      wetgridunc(ix,jy,ks,kp,nunc,nage)= &
     99           wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
     100  else
     101    if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    98102       (jy.le.numygrid-1)) then
    99     w=wx*wy
     103      w=wx*wy
    100104      wetgridunc(ix,jy,ks,kp,nunc,nage)= &
    101105           wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
    102   endif
     106    endif
    103107
    104   if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
     108    if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
    105109       (jyp.le.numygrid-1)) then
    106     w=(1.-wx)*(1.-wy)
     110      w=(1.-wx)*(1.-wy)
    107111      wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &
    108112           wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    109   endif
     113    endif
    110114
    111   if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
     115    if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
    112116       (jy.le.numygrid-1)) then
    113     w=(1.-wx)*wy
     117      w=(1.-wx)*wy
    114118      wetgridunc(ixp,jy,ks,kp,nunc,nage)= &
    115119           wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
    116   endif
     120    endif
    117121
    118   if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
     122    if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
    119123       (jyp.le.numygrid-1)) then
    120     w=wx*(1.-wy)
     124      w=wx*(1.-wy)
    121125      wetgridunc(ix,jyp,ks,kp,nunc,nage)= &
    122126           wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    123   endif
     127    endif
     128  endif
    124129  end do
    125130  end if
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG