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


Ignore:
Files:
3 deleted
19 edited

Legend:

Unmodified
Added
Removed

  • options/SPECIES/spec_overview

    • Property mode changed from 100755 to 100644

  • src/advance.f90

    r94735e2 ra652cd5  
    117117  real :: usigold,vsigold,wsigold,r,rs
    118118  real :: uold,vold,wold,vdepo(maxspec)
    119   real :: h1(2)
    120119  !real uprof(nzmax),vprof(nzmax),wprof(nzmax)
    121120  !real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
    … …  
    224223
    225224
    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 
    245225  ! Compute maximum mixing height around particle position
    246226  !*******************************************************
    … …  
    249229  if (ngrid.le.0) then
    250230    do k=1,2
    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
     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
    264237    end do
    265238    tropop=tropopause(nix,njy,1,1)
    … …  
    276249  endif
    277250
    278   if (interpolhmix) h=(h1(1)*dt2+h1(2)*dt1)*dtt
    279251  zeta=zt/h
    280252
    … …  
    474446        endif
    475447
    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 
    484448  !****************************************************
    485449  ! Compute turbulent vertical displacement of particle
    … …  
    683647  endif
    684648
    685   if (turboff) then
    686 !sec switch off turbulence
    687     ux=0.0
    688     vy=0.0
    689     wp=0.0
    690   endif
    691649
    692650  ! If particle represents only a single species, add gravitational settling

  • src/com_mod.f90

    r5888298 rc8fc724  
    577577  real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
    578578  !real outheight(maxzgrid),outheighthalf(maxzgrid)
    579 
    580579  logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
    581580       & OHREA,ASSSPEC
    582   logical :: DRYBKDEP,WETBKDEP
    583581
    584582  ! numxgrid,numygrid       number of grid points in x,y-direction
    … …  
    600598  ! OHREA                   .true., if OH reaction is switched on
    601599  ! 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
    603600  !                    (i.e. transfer of mass between these two occurs
    604601
    … …  
    671668  real, allocatable, dimension(:) :: ztra1
    672669  real, allocatable, dimension(:,:) :: xmass1
    673   real, allocatable, dimension(:,:) :: xscav_frac1
    674670
    675671  ! eso: Moved from timemanager
    … …  
    692688  ! xtra1,ytra1,ztra1       spatial positions of the particles
    693689  ! xmass1 [kg]             particle masses
    694   ! xscav_frac1             fraction of particle masse which has been scavenged at receptor
    695  
     690 
    696691
    697692
    … …  
    755750  integer :: mpi_mode=0 ! .gt. 0 if running MPI version
    756751  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  
    762752 
    763753contains

  • src/conccalc.f90

    r9287c01 r4c64400  
    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 !  integer xscav_count
     61
    6262
    6363  ! For forward simulations, make a loop over the number of species;
    … …  
    6565  ! releasepoints
    6666  !***************************************************************************
    67 !  xscav_count=0
     67
    6868  do i=1,numpart
    6969    if (itra1(i).ne.itime) goto 20
    … …  
    767633   continue
    7777
    78 !  if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
    79            
     78
    8079  ! For special runs, interpolate the air density to the particle position
    8180  !************************************************************************
    … …  
    173172      if (yl.lt.0.) jy=jy-1
    174173
     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           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)= &
     188          do ks=1,nspec
     189            gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    197190                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    198191                 xmass1(i,ks)/rhoi*weight
    199              end do
    200           endif
    201         endif
    202 
    203       else                                 ! attribution via uniform kernel
     192          end do
     193        endif
     194
     195      else                                 ! attribution via uniform kernel
    204196
    205197        ddx=xl-real(ix)                   ! distance to left cell border
    … …  
    228220          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    229221            w=wx*wy
    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)= &
     222            do ks=1,nspec
     223              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    239224                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    240225                   xmass1(i,ks)/rhoi*weight*w
    241                end do
    242             endif
     226            end do
    243227          endif
    244228
    245229          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    246230            w=wx*(1.-wy)
    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)= &
     231            do ks=1,nspec
     232              gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    256233                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    257234                   xmass1(i,ks)/rhoi*weight*w
    258                end do
    259              endif
    260           endif
    261         endif !ix ge 0
     235            end do
     236          endif
     237        endif
    262238
    263239
    … …  
    265241          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
    266242            w=(1.-wx)*(1.-wy)
    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)= &
     243            do ks=1,nspec
     244              gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    276245                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    277246                   xmass1(i,ks)/rhoi*weight*w
    278                end do
    279             endif
     247            end do
    280248          endif
    281249
    282250          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
    283251            w=(1.-wx)*wy
    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)= &
     252            do ks=1,nspec
     253              gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    293254                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    294255                   xmass1(i,ks)/rhoi*weight*w
    295                end do
    296             endif
    297           endif
    298         endif !ixp ge 0
    299      endif
     256            end do
     257          endif
     258        endif
     259      endif
     260
     261
    300262
    301263  !************************************
    … …  
    320282        if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
    321283             (xl.gt.real(numxgridn-1)-0.5).or. &
    322              (yl.gt.real(numygridn-1)-0.5).or.(.not.usekernel)) then             ! no kernel, direct attribution to grid cell
     284             (yl.gt.real(numygridn-1)-0.5)) then             ! no kernel, direct attribution to grid cell
    323285          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
    324286               (jy.le.numygridn-1)) then
    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)= &
     287            do ks=1,nspec
     288              griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    334289                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    335290                   xmass1(i,ks)/rhoi*weight
    336                end do
    337             endif
     291            end do
    338292          endif
    339293
    … …  
    365319            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    366320              w=wx*wy
    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)= &
     321              do ks=1,nspec
     322                griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    376323                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    377324                     xmass1(i,ks)/rhoi*weight*w
    378                  end do
    379               endif
     325              end do
    380326            endif
    381327
    382328            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    383329              w=wx*(1.-wy)
    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)= &
     330              do ks=1,nspec
     331                griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    393332                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    394333                     xmass1(i,ks)/rhoi*weight*w
    395                  end do
    396               endif
     334              end do
    397335            endif
    398336          endif
    … …  
    402340            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
    403341              w=(1.-wx)*(1.-wy)
    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)= &
     342              do ks=1,nspec
     343                griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
    413344                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
    414345                     xmass1(i,ks)/rhoi*weight*w
    415                  end do
    416               endif
     346              end do
    417347            endif
    418348
    419349            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
    420350              w=(1.-wx)*wy
    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)= &
     351              do ks=1,nspec
     352                griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
    430353                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
    431354                     xmass1(i,ks)/rhoi*weight*w
    432                  end do
    433               endif
     355              end do
    434356            endif
    435357          endif
    436358        endif
     359
    437360      endif
    438361    endif
    43936220  continue
    440363  end do
    441 !  write(*,*) 'xscav count:',xscav_count
    442364
    443365  !***********************************************************************

  • src/concoutput.f90

    r16b61a5 r16b61a5  
    246246
    247247    write(anspec,'(i3.3)') ks
    248 
    249     if (DRYBKDEP.or.WETBKDEP) then !scavdep output
    250       if (DRYBKDEP) &
    251       open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
     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// &
    252261           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// &
    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
     262
     263      write(unitoutgridppt) itime
     264    endif
    274265
    275266    do kp=1,maxpointspec_act
    … …  
    351342! Concentration output
    352343!*********************
    353         if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
     344        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
    354345
    355346! Wet deposition

  • src/drydepokernel.f90

    r1c0d5e6 r4c64400  
    102102   if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
    103103
    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. &
     104   if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    109105        (jy.le.numygrid-1)) then
    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
     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
    114111
    115     if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
     112  if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
    116113       (jyp.le.numygrid-1)) then
    117114    w=(1.-wx)*(1.-wy)
    118115      drygridunc(ixp,jyp,ks,kp,nunc,nage)= &
    119116           drygridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    120     endif
     117  endif
    121118
    122     if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
     119  if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
    123120       (jy.le.numygrid-1)) then
    124       w=(1.-wx)*wy
     121    w=(1.-wx)*wy
    125122      drygridunc(ixp,jy,ks,kp,nunc,nage)= &
    126123           drygridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
    127     endif
     124  endif
    128125
    129     if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
     126  if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
    130127       (jyp.le.numygrid-1)) then
    131       w=wx*(1.-wy)
     128    w=wx*(1.-wy)
    132129      drygridunc(ix,jyp,ks,kp,nunc,nage)= &
    133130           drygridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    134     endif
     131  endif
    135132
    136     endif ! kernel
    137     endif ! deposit>0
     133  endif
    138134
    139135  end do

  • src/ecmwf_mod.f90

    r842074e r62e65c7  
    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 :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 0.5
    47    integer,parameter :: nxshift=359
    48 !  integer,parameter :: nxshift=718
     45  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
     46  integer,parameter :: nxshift=359
    4947!  integer,parameter :: nxshift=0
    5048

  • src/getfields.f90

    rd1a8707 r8a65cb0  
    13813840  indmin=indj
    139139
    140    if (WETBKDEP) then
    141         call writeprecip(itime,memind(1))
    142    endif
    143 
    144140  else
    145141
    … …  
    17316960  indmin=indj
    174170
    175    if (WETBKDEP) then
    176         call writeprecip(itime,memind(1))
    177    endif
    178 
    179171  endif
    180172

  • src/interpol_rain.f90

    rc7e771d r8a65cb0  
    2121
    2222subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
    23        ny,iwftouse,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
     23       ny,memind,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
    8079
    8180
    … …  
    114113  !***********************
    115114
    116 !  do m=1,2
    117     indexh=iwftouse
     115  do m=1,2
     116    indexh=memind(m)
    118117
    119     y1(1)=p1*yy1(ix ,jy ,level,indexh) &
     118
     119    y1(m)=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(1)=p1*yy2(ix ,jy ,level,indexh) &
     123    y2(m)=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(1)=p1*yy3(ix ,jy ,level,indexh) &
     127    y3(m)=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) - skip to be consistent with clouds
     135  ! 2.) Temporal interpolation (linear)
    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)
    149146
    150147end subroutine interpol_rain

  • src/interpol_vdep.f90

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

  • src/makefile

    rd1a8707 r4c64400  
    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)  #-Warray-bounds -fcheck=all # -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) # -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       \
    146 writeprecip.o \
    147146writeheader_surf.o      assignland.o\
    148147part0.o                 gethourlyOH.o\
    … …  
    153152erf.o                   readavailable.o \
    154153ew.o                    readreleases.o  \
    155 readdepo.o              get_vdep_prob.o   \
    156 get_wetscav.o   \
     154readdepo.o \
    157155psim.o                  outgrid_init.o  \
    158156outgrid_init_nest.o     \
    … …  
    273271
    274272## DEPENDENCIES
    275 get_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
    277273advance.o: cmapf_mod.o com_mod.o hanna_mod.o interpol_mod.o par_mod.o \
    278274        point_mod.o random_mod.o
    … …  
    428424verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
    429425verttransform_nests.o: com_mod.o par_mod.o
    430 get_wetscav.o: com_mod.o par_mod.o point_mod.o
    431426wetdepo.o: com_mod.o par_mod.o point_mod.o
    432427wetdepokernel.o: com_mod.o par_mod.o unc_mod.o
    433428wetdepokernel_nest.o: com_mod.o par_mod.o unc_mod.o
    434429writeheader.o: com_mod.o outg_mod.o par_mod.o point_mod.o
    435 writeprecip.o: com_mod.o par_mod.o point_mod.o
    436430writeheader_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o
    437431writeheader_nest_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o

  • src/par_mod.f90

    rd1a8707 r4c64400  
    194194  integer,parameter :: maxpart=10000000
    195195  integer,parameter :: maxspec=4
    196 
    197196  real,parameter :: minmass=0.0001
    198197
    … …  
    259258  integer,parameter :: unitlsm=1, unitsurfdata=1, unitland=1, unitwesely=1
    260259  integer,parameter :: unitOH=1
    261   integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95, unitprecip=101
     260  integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
    262261  integer,parameter :: unitboundcond=89
    263262  integer,parameter :: unittmp=101

  • src/readcommand.f90

    r8ee24a5 r4c64400  
    298298  !Af          1 = mass units
    299299  !Af          2 = mass mixing ratio units
    300   !            3 = wet deposition in outputfield
    301   !            4 = dry deposition in outputfield
    302300
    303301  if ( ldirect .eq. 1 ) then  ! FWD-Run
    … …  
    322320     endif
    323321  !Af set release-switch
    324      WETBKDEP=.false.
    325      DRYBKDEP=.false.
    326      select case (ind_receptor)
    327      case (1)  !  1 .. concentration at receptor
     322     if (ind_receptor .eq. 1) then !mass
    328323        ind_rel = 1
    329      case (2)  !  2 .. mixing ratio at receptor
     324     else ! mass mix
    330325        ind_rel = 0
    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
     326     endif
    346327  endif
    347328
    … …  
    405386  !**********************************************************************
    406387
    407   if ((iout.lt.1).or.(iout.gt.6)) then
     388  if ((iout.lt.1).or.(iout.gt.5)) then
    408389    write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
    409     write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR        #### '
     390    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR       #### '
    410391    write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
    411392    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 
    537523  ! Check whether x coordinates of release point are within model domain
    538524  !*********************************************************************

  • src/releaseparticles.f90

    r75a4ded r8a65cb0  
    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.
    180178  !*****************************************************************************
    181179            do k=1,nspec
    182180               xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
    183181                    *timecorrect(k)/average_timecorrect
    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
     182  !            write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
    190183  ! Assign certain properties to particle
    191184  !**************************************
    … …  
    209202
    210203            ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
     204
    211205  ! Interpolation of topography and density
    212206  !****************************************
    … …  
    336330  !Af ind_rel is defined in readcommand.f
    337331
    338             if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
     332            if (ind_rel .eq. 1) then
    339333
    340334  ! Interpolate the air density
    … …  
    381375            endif
    382376
     377
    383378            numpart=max(numpart,ipart)
    384379            goto 34      ! Storage space has been found, stop searching
    385380          endif
    386         end do  ! i=1:numpoint
     381        end do
    387382        if (ipart.gt.maxpart) goto 996
    388383
    38938434      minpart=ipart+1
    390       end do ! ipart=minpart,maxpart
    391       endif ! j=1,numrel
     385      end do
     386      endif
    392387  end do
    393388

  • src/timemanager.f90

    rd1a8707 r18adf60  
    105105! integer :: ksp
    106106  integer :: loutnext,loutstart,loutend
    107   integer :: ix,jy,ldeltat,itage,nage,idummy
     107  integer :: ix,jy,ldeltat,itage,nage
    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,wetscav(maxspec)
     109  real :: outnum,weight,prob(maxspec),decfact
    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)
    117116  real, parameter :: e_inv = 1.0/exp(1.0)
    118 
    119117  !double precision xm(maxspec,maxpointspec_act),
    120118  !    +                 xm_depw(maxspec,maxpointspec_act),
    … …  
    147145!CGZ-lifetime: set lifetime to 0
    148146 
    149   if (.not.usekernel) write(*,*) 'Not using the kernel'
    150   if (turboff) write(*,*) 'Turbulence switched off'
    151 
    152147  write(*,46) float(itime)/3600,itime,numpart
    153148
    … …  
    546541        zold=ztra1(j)
    547542
    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 
    584543  ! Integrate Lagevin equation for lsynctime seconds
    585544  !*************************************************
    586545
    587         if (verbosity.gt.0) then
    588            if (j.eq.1) then
    589              write (*,*) 'timemanager> call advance'
    590            endif     
    591         endif
    592      
    593546        call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
    594547             us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
    595548             cbt(j))
    596 !        write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
    597549
    598550  ! Calculate the gross fluxes across layer interfaces

  • src/verttransform.f90

    r6481010 r341f4b7  
    9191
    9292  logical :: init = .true.
    93   logical :: init_w = .false.
    94   logical :: init_r = .false.
    95 
    9693
    9794  !ZHG SEP 2014 tests 
    … …  
    106103  ! CHARACTER(LEN=3)  :: aspec
    107104  ! integer :: virr=0
    108   !real :: tot_cloud_h
    109   !real :: dbg_height(nzmax)
     105  real :: tot_cloud_h
     106  real :: dbg_height(nzmax)
    110107!ZHG
    111108
    … …  
    126123  if (init) then
    127124
    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 
    141125! Search for a point with high surface pressure (i.e. not above significant topography)
    142126! Then, use this point to construct a reference z profile, to be used at all times
    … …  
    177161    end do
    178162
    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
    189163
    190164! Determine highest levels that can be within PBL
    … …  
    204178    init=.false.
    205179
    206 !    dbg_height = height
     180    dbg_height = height
    207181
    208182  endif
    … …  
    624598        convp=convprec(ix,jy,1,n)
    625599        prec=lsp+convp
    626 !        tot_cloud_h=0
     600        tot_cloud_h=0
    627601! Find clouds in the vertical
    628602        do kz=1, nz-1 !go from top to bottom
    … …  
    630604! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
    631605            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
    632 !            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
     606            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
    633607           
    634608!            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        *
    4447! ix,jy              indices of output grid cell for each particle           *
    4548! itime [s]          actual simulation time [s]                              *
    4649! jpart              particle index                                          *
    4750! 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)  *
    4853! loutnext [s]       time for which gridded deposition is next output        *
    4954! loutstep [s]       interval at which gridded deposition is output          *
     55! lsp [mm/h]         large scale precipitation rate                          *
    5056! ltsample [s]       interval over which mass is deposited                   *
     57! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
    5158! wetdeposit         mass that is wet deposited                              *
    5259! wetgrid            accumulated deposited mass on output grid               *
    … …  
    6471
    6572  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    66   integer :: itage,nage,il,interp_time, n
     73  integer :: ngrid,itage,nage,hz,il,interp_time, n
     74  integer(kind=1) :: clouds_v
    6775  integer :: ks, kp
    68   integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count
    69   real :: grfraction(3),wetscav
     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
    7080  real :: wetdeposit(maxspec),restmass
    7181  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
    7296
    7397! Compute interval since radioactive decay of deposited mass was computed
    … …  
    95119    endif
    96120
    97 ! Determine age class of the particle - nage is used for the kernel
     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
     12633  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
     14023  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!**********************************************************
     255
     256    do ks=1,nspec      ! loop over species
     257      wetdeposit(ks)=0.
     258      wetscav=0.
     259
     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
    98276!******************************************************************
    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
    104 
    105     do ks=1,nspec      ! loop over species
    106 
    107       if (WETDEPSPEC(ks).eqv..false.) cycle
     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
    108381
    109382!**************************************************
    110383! CALCULATE DEPOSITION
    111384!**************************************************
    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      
     385!     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
     386!     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
    122387
    123388      if (wetscav.gt.0.) then
    124389        wetdeposit(ks)=xmass1(jpart,ks)* &
    125390             (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
    126395      else ! if no scavenging
    127396        wetdeposit(ks)=0.
    128397      endif
    129  
     398
    130399      restmass = xmass1(jpart,ks)-wetdeposit(ks)
    131400      if (ioutputforeachrelease.eq.1) then
    … …  
    148417      endif
    149418
    150 !    endif ! no deposition
    151     end do ! loop over species
     419
     420    end do !all species
    152421
    153422! Sabine Eckhardt, June 2008 create deposition runs only for forward runs

  • src/wetdepokernel.f90

    r1c0d5e6 r4c64400  
    9595  do ks=1,nspec
    9696
    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. &
     97  if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    10298       (jy.le.numygrid-1)) then
    103       w=wx*wy
     99    w=wx*wy
    104100      wetgridunc(ix,jy,ks,kp,nunc,nage)= &
    105101           wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
    106     endif
     102  endif
    107103
    108     if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
     104  if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
    109105       (jyp.le.numygrid-1)) then
    110       w=(1.-wx)*(1.-wy)
     106    w=(1.-wx)*(1.-wy)
    111107      wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &
    112108           wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    113     endif
     109  endif
    114110
    115     if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
     111  if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
    116112       (jy.le.numygrid-1)) then
    117       w=(1.-wx)*wy
     113    w=(1.-wx)*wy
    118114      wetgridunc(ixp,jy,ks,kp,nunc,nage)= &
    119115           wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
    120     endif
     116  endif
    121117
    122     if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
     118  if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
    123119       (jyp.le.numygrid-1)) then
    124       w=wx*(1.-wy)
     120    w=wx*(1.-wy)
    125121      wetgridunc(ix,jyp,ks,kp,nunc,nage)= &
    126122           wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    127     endif
    128   endif
     123  endif
    129124  end do
    130125  end if
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG