ecearth3/sources/tm5mp/proj/cb05/wet_deposition.F90

! 
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
#include "output.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      WET_DEPOSITION
!
! !DESCRIPTION: holds convective precipitation (CP) and large scale
!               precipitation (LSP) budgets variables.
!
!               has all routines to deal with LSP, since it is done here. CP
!               is done in convection.
!
!      **** THIS VERSION DIFFERS FROM THE BASE IN TWO THRESHOLD VALUES ****
!      
!      (1) Scale factor relative to 100% scavenging (of HNO3) for scavenging
!          type = 2 (ie variable factor, using henry solubility) is non-zero and
!          computed if henry coeff > 1 (instead of 10 in base code).
!          
!      (2) Wet removal rates: in case of in-cloud scavenging, test on cloud
!          cover has a threshold of 0.02 instead of 0.05
!\\             
!\\
! !INTERFACE: 
!
MODULE WET_DEPOSITION
  !
  ! !USES:
  !
  use GO,           only : gol, goErr, goPr
  use dims,         only : nregions
  use chem_param,   only : ntracet
  use global_types, only : d3_data
  use tm5_distgrid, only : dgrid, Get_DistGrid, reduce, gather
  use ParTools,     only : isRoot

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public  ::  Wet_Deposition_Init, Wet_Deposition_Done
  public  ::  calc_cvsfac, calculate_lsp_scav, lspscav
  !
  ! !PUBLIC DATA MEMBERS:
  !
  real, public ::  cvsfac(ntracet) = 0.0 ! scavenging efficiencies, used in convection
  real, public ::  cp_scale = 0.5        ! default amount of rain (converted to m/s) with maximum CP removal on 1x1 degree

#ifdef with_budgets
  
  real, dimension(nregions), public ::  sum_wetdep ! global wet dep budget for tracer #1 (includes both LSP and CP; meaningful on root only)

  type, public :: buddep_data
     ! budgets in each column, split vertically...
     real,dimension(:,:,:,:),pointer    :: lsp     ! (i, j, nbud_vg, ntracet) computed here
     real,dimension(:,:,:,:),pointer    :: cp      ! (i, j, nbud_vg, ntracet) computed in convection
  end type buddep_data

  type(buddep_data), dimension(nregions), target, public  :: buddep_dat ! ... for each region
#endif

  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'wet_deposition'
  !
  !                                              Large scale scavenging coefficients [s-1]
  type(d3_data),dimension(nregions) :: rloss1    ! 1: incloud completely soluble gas        
!>>>TvN
  type(d3_data),dimension(nregions) :: rloss1_m7 ! as 1, but with ice as efficient as water 
                                                 ! now needed for M7 aerosols
!<<<TvN
  type(d3_data),dimension(nregions) :: rloss2    ! 2: below cloud completely soluble gas    
  type(d3_data),dimension(nregions) :: rloss3    ! 3: below cloud bulk aerosol (Whitby distribution)
  !
  ! rain-out can not be higher than maximum level of convection
  ! thus maximum level of convection lmax_conv(=>ntot_ed) is used
  !
  !
  ! used from chem_param:
  !
  !   nscav       : selected species for scavenging
  !   nscav_index : index for scavenging:
  !   nscav_type  : type of scavenging:
  !                   0 no scavenging
  !                   1 scavenging 100 % solubility assumed
  !                   2 scavenging henry solubility assumed
  !                   3 scavenging, bulk aerosol removal assumed
  !                   4 scavenging, special case for SO2 with aq phase diss.
  !                5-11 scavenging, M7 aerosol removal
  !
  !----------------------------------------------
  ! acidity needed for explicit calculation of reactive removal of SO2.
  ! Parameterisation calculates reaction of SO2 with H2O2 and O3.
  ! Not yet implemented: for information Frank Dentener    ! see routine wetS
  ! nacid       : selected species for acidity
  ! nacid_index : indices of species for acidity : iso4,imsa,ihno3,inh3,inh4
  ! nacid_eq    : equivalents acid per mole
  ! integer,parameter                    :: nacid=5
  ! integer,parameter,dimension(nacid)   :: &
  !                          nacid_index    = (/iso4,imsa,ihno3,ino3_a,inh3,inh4/)
  ! integer,parameter,dimension(nacid)   :: &
  !                             nacid_eq    = (/   2,   1,    1,  1   ,  -1,  -1/)
  !----------------------------------------------
  !
  ! !REVISION HISTORY:
  !   version March 2003, adapted for TM5 MK from KNMI version
  !   29 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    WET_DEPOSITION_INIT
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE WET_DEPOSITION_INIT( status )
    !
    ! !USES:
    !
    use GO,            only : TrcFile, Init, Done, ReadRc
    use dims,          only : lmax_conv
    use global_data,   only : rcfile
    use meteodata,     only : Set, temper_dat, lwc_dat, iwc_dat, cc_dat, lsp_dat
#ifdef with_budgets
    use budget_global, only : nbud_vg
#endif
    use chem_param,    only : ntrace
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !   29 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Wet_Deposition_Init'

    integer           ::  region, imr,jmr,lmr, i1, i2, j1, j2
    type(TrcFile)     ::  rcF

    ! --- begin ------------------------------------

    ! select meteo to be used
    do region = 1, nregions
       call Set( temper_dat(region), status, used=.true. )
       call Set(    lwc_dat(region), status, used=.true. )
       call Set(    iwc_dat(region), status, used=.true. )
       call Set(     cc_dat(region), status, used=.true. )
       call Set(    lsp_dat(region), status, used=.true. )
    end do

    ! allocate
    do region = 1, nregions
       lmr = lmax_conv
       CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       allocate(rloss1(region)%d3(i1:i2, j1:j2, lmr))
!>>>TvN
       allocate(rloss1_m7(region)%d3(i1:i2, j1:j2, lmr))
!<<<TvN
       allocate(rloss2(region)%d3(i1:i2, j1:j2, lmr))
       allocate(rloss3(region)%d3(i1:i2, j1:j2, lmr))
    end do

    ! read settings from rcfile:
    !  o scaling factor wet removal  (1.-exp(-cp/cp_scale))
    !    cp_scale  : 0.5
    !   (see convection.F90 and wet_deposition.F90) 
    !
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'proces.wet_removal.cp_scale', cp_scale, status, default=0.5 )
    IF_NOTOK_RETURN(status=1)
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    write (gol,*) 'maximum removal CP precip on 1x1 at rain rate (mm/hr) :', cp_scale; call goPr
    cp_scale = cp_scale/(1e3*3600.)  ! to m/s!

    ! Calculate removal rates by convective precipitation
    ! It was commented before : JEW:called too early, KHENRY must be calculated for some species first
    ! Back here, since KHENRY are now calculated before hand with a call to "rates" in sources_sinks_init 
    call calc_cvsfac

    ! budgets
#ifdef with_budgets
    do region = 1, nregions
       CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
       sum_wetdep(region) = 0.0
       allocate( buddep_dat(region)%lsp(i1:i2, j1:j2, nbud_vg, ntracet))
       buddep_dat(region)%lsp = 0.0
       allocate( buddep_dat(region)%cp (i1:i2, j1:j2, nbud_vg, ntracet))
       buddep_dat(region)%cp  = 0.0
    enddo
#endif

    ! ok
    status = 0

  END SUBROUTINE WET_DEPOSITION_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    WET_DEPOSITION_DONE
  !
  ! !DESCRIPTION: deallocate scavenging coeff. & write wet dep and convection
  !               budgets into file.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE WET_DEPOSITION_DONE( status )
    !
    ! !USES:
    !
    use dims,          only : nregions, im, jm, lm
    use chem_param,    only : ntracet
    use partools,      only : par_reduce_element
#ifdef with_budgets
    use budget_global, only : budget_file_global, nbud_vg, budg_dat, nbudg, NHAB
    use budget_global, only : budconvg
#ifdef with_hdf4          
    use file_hdf,      only : THdfFile, TSds
    use file_hdf,      only : Init, Done, WriteAttribute, WriteData, SetDim
#endif
    use Dims,          only : region_name
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !   29 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Wet_Deposition_Done'

#ifdef with_budgets
#ifdef with_hdf4          
    type(THdfFile)                           :: io
    type(TSds)                               :: sds
#endif
    integer                                  :: i1, i2, j1, j2
    integer                                  :: nsend,j,i,n,nzone,nzone_v
    real,dimension(:,:,:,:),allocatable      :: budget4d
    real,dimension(nbudg,nbud_vg,ntracet)    :: budwet_cp,     budwet_lsp     ! for one MPI tile
    real,dimension(nbudg,nbud_vg,ntracet)    :: budconvg_all,  budwet_cp_all, budwet_lsp_all ! 
#endif
    integer                                  :: region,l

    ! --- begin ----------------------------------

    do region = 1, nregions
       deallocate( rloss1(region)%d3 )
       deallocate( rloss1_m7(region)%d3 )
       deallocate( rloss2(region)%d3 )
       deallocate( rloss3(region)%d3 )
    end do

#ifdef with_budgets

#ifdef with_hdf4          
    if(isRoot) then

       call Init(io, budget_file_global, 'write', status)
       IF_NOTOK_RETURN(status=1)

       call WriteAttribute(io,'sum_wetdep',sum_wetdep,status)
       IF_NOTOK_RETURN(status=1)

       call WriteAttribute(io,'cvsfac',cvsfac,status)
       IF_NOTOK_RETURN(status=1)
       
    endif
#endif
    budwet_cp  = 0.0
    budwet_lsp = 0.0


    ! =============================== Aggregate and write column-like Wet Dep budgets
    REG : do region=1,nregions   

       !---- horizontally aggregates LSP and CP budgets
       CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       do n=1, ntracet
       do l=1, nbud_vg
       do j=j1,j2
       do i=i1,i2
          nzone = budg_dat(region)%nzong(i,j)
          budwet_lsp(nzone,l,n) = budwet_lsp(nzone,l,n) + buddep_dat(region)%lsp(i,j,l,n)
          budwet_cp(nzone,l,n)  = budwet_cp(nzone,l,n)  + buddep_dat(region)%cp(i,j,l,n)
       end do
       end do
       end do
       end do

       !-- write Non-Horizontally-Aggregated-Budgets
       if (NHAB) then

          ! global array to gather data
          if (isRoot)then
             allocate(budget4d(im(region), jm(region), nbud_vg, ntracet))
          else
             allocate(budget4d(1,1,1,1))
          end if

          ! gather and write column-like wet dep LSP budgets
          CALL GATHER( dgrid(region), buddep_dat(region)%lsp, budget4d, 0, status)

#ifdef with_hdf4          
          if(isRoot) then
             call Init(Sds,io, 'buddep_dat_lsp_'//region_name(region),(/im(region),jm(region),nbud_vg,ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
             call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
             call SetDim(Sds, 2, 'nbud_vg','budget level', (/(j, j=1,nbud_vg)/), status)
             call SetDim(Sds, 3, 'ntracet','tracer index', (/(j, j=1,ntracet)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget4d, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif
#endif
          ! gather and write column-like wet dep CP budgets
          CALL GATHER( dgrid(region), buddep_dat(region)%cp, budget4d, 0, status)

#ifdef with_hdf4          
          if(isRoot) then
             call Init(Sds,io, 'buddep_dat_cp_'//region_name(region),(/im(region),jm(region),nbud_vg,ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
             call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
             call SetDim(Sds, 2, 'nbud_vg','budget level', (/(j, j=1,nbud_vg)/), status)
             call SetDim(Sds, 3, 'ntracet','tracer index', (/(j, j=1,ntracet)/), status)
             call WriteData(Sds, budget4d, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif
#endif
          deallocate(budget4d)

       endif  ! NHAB
    enddo REG ! regions

    
    !================== Write horizontally aggregated Wet Dep (& convection) budgets

    ! Sum up contribution from each proc into root array
    
    call PAR_REDUCE_ELEMENT( budwet_cp, 'SUM', budwet_cp_all, status )
    IF_NOTOK_RETURN(status=1)
    
    call PAR_REDUCE_ELEMENT( budwet_lsp, 'SUM', budwet_lsp_all, status )
    IF_NOTOK_RETURN(status=1)

    call PAR_REDUCE_ELEMENT( budconvg, 'SUM', budconvg_all, status )
    IF_NOTOK_RETURN(status=1)

    if (isRoot) then

       ! update convection budget
       budconvg_all(:,:,:) = budconvg_all(:,:,:) + budwet_cp(:,:,:)
       
#ifdef with_hdf4          
       call Init(Sds, io, 'budconv_cp', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call SetDim(Sds, 0, 'nbudg', 'horizontal region', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical layer', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'ntracet', 'tracer index', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'comment', 'Convection budget corrected for cp removal', status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, budconvg_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budwet_cp', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call SetDim(Sds, 0, 'nbudg', 'horizontal region', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical layer', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'ntracet', 'tracer index', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, budwet_cp_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budwet_lsp', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call SetDim(Sds, 0, 'nbudg', 'horizontal region', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical layer', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'ntracet', 'tracer index', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, budwet_lsp_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)
       call Done(io, status)
       IF_NOTOK_RETURN(status=1)
#endif
    endif

    do region = 1, nregions
       deallocate( buddep_dat(region)%lsp ) 
       deallocate( buddep_dat(region)%cp ) 
    end do

#endif  /* WITH_BUDGET */

    ! ok
    status = 0

  END SUBROUTINE WET_DEPOSITION_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CALC_CVSFAC
  !
  ! !DESCRIPTION: lookup tables, calculate scale factor for convective scavenging
  !               relative to 100% scavenging (of HNO3),
  !               assuming constant temperature in convective updraft of 273K.
  !               
  ! Factors for different scavenging types are:
  ! 
  ! 0) no/low solubility  cvsfac=0
  ! 1) high solubility    cvsfac=1
  ! 2) henry solubility   cvsfac=variable
  ! 3) bulk aerosol       cvsfac=0.99
  !      For the moment a value of 0.99 is assumed for bulk aerosol.
  !      This is the value for the soluble accumulation mode from Stier et al. (JGR, 2005).
  !      and presents an upper boundary for bulk aerosols.
  !      (A value ~0.9 would probably make more sense, but this would need to be justified 
  !       with some quantitative argument.)
  ! 4) SO2 dissociation   cvsfac=variable depending on T and pH and 
  !                              dissociation of H2SO3<-->HSO3(-)<--->SO3(2-)
  !                              for convective clouds T=273K and pH=5 is chosen
  ! 5-11) M7 modes        cvsfac set equal to the scavenging parameters from Stier et al. (2005)
  !                              for convective in-cloud scavenging
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CALC_CVSFAC
    !
    ! !USES:
    !
    use chem_param, only : nscav, nscav_index, nscav_type
    use chem_param, only : henry, ntlow, ntemp
    use chem_param, only : names
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: iscav,n,k
    real    :: rtl, heff

    cvsfac=0.0

    do iscav=1,nscav
       n=nscav_index(iscav)
       ! skip dummy tracers ..
       if ( n < 0 ) cycle

       ! fill cvsfac given scavenging type:
       select case(nscav_type(iscav))
       case(0)
          cvsfac(n) = 0.0
!>>>TvN
       !case(1,3)
       case(1)
!<<< TvN
          cvsfac(n) = 1.0
       case(2)
          rtl=8.3148e-8*273.   ! phase factor: ratio of aq. phase conc. to total conc.
          ! assuming LWC of 1e-6
          k=nint(273.-float(ntlow))
          k = min(max(1,k),ntemp)

          if ( henry(n,k) > 1. ) then
             cvsfac(n) =  rtl*henry(n,k)/(1.+rtl*henry(n,k))
          else
             cvsfac(n) = 0.0
          end if
!>>>TvN
       case(3)
         cvsfac(n) = 0.99
!<<<TvN
       case(4)
          rtl=8.3148e-8*273.   ! phase factor: ratio of aq. phase conc. to total conc.o
          ! assuming LWC of 1e-6
          k=nint(273.-float(ntlow))
          k = min(max(1,k),ntemp)
          heff = henry(n,k)*3.2e3   ! 3.2e3 factor is dissociation of SO2 and HSO3- at pH=5
          cvsfac(n) =  rtl*heff/(1.+rtl*heff)
!>>>TvN
       case(5)   ! soluble nu
          !cvsfac(n) = 1.0
          cvsfac(n) = 0.2
       case(6)   ! soluble ai
          !cvsfac(n) = 1.0
          cvsfac(n) = 0.6
       case(7)   ! soluble ac
          !cvsfac(n) = 1.0
          cvsfac(n) = 0.99
       case(8)   ! soluble co
          !cvsfac(n) = 1.0
          cvsfac(n) = 0.99
       case(9)   ! insoluble ai
          !cvsfac(n) = 1.0
          cvsfac(n) = 0.2
       case(10)   ! insoluble ac
          !cvsfac(n) = 1.0
          cvsfac(n) = 0.4
       case(11)   ! insoluble co
          !cvsfac(n) = 1.0
          cvsfac(n) = 0.4
!<<<TvN
       case default
          cvsfac(n) = 0.0
       end select
    end do

    ! info ...
    do n = 1, ntracet
       if ( cvsfac(n) > 0.0 ) then
          write (gol,'("  calc_cvsfac: Scavenging factor species ",i3,x,a,"; factor : ",e10.3)') &
               n, names(n), cvsfac(n); call goPr
       end if
    end do

  END SUBROUTINE CALC_CVSFAC
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    LSPSCAV
  !
  ! !DESCRIPTION: Calculation of wet removal by large scale precipitation of
  !               soluble tracers
  !
  ! Remove tracers with previously calculated rainout rate [s-1]
  ! separately for in- and below cloud scavenging and only for the
  ! cloud covered fraction of the gridcell
  !      
  ! Reference:
  !      Langner and Rodhe (1990)
  !      Junge (1959)  
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE LSPSCAV( region )
    !
    ! !USES:
    !
    use binas,         only : rgas
    use dims,          only : im, jm, lm, nchem
    use dims,          only : tref, lmax_conv, dy
    use chem_param,    only : ntrace, ntracet, henry, ntlow, ra, ntemp
    use chem_param,    only : xmhno3 !, ihno3
    use chem_param,    only : nscav, nscav_index, nscav_type
#ifdef with_m7
    use chem_param,    only : inus_n, iais_n, iacs_n, icos_n, iaii_n, iaci_n, icoi_n
#endif
    use meteodata,     only : temper_dat, cc_dat
    use global_data,   only : mass_dat
#ifdef with_budgets
    use budget_global, only : nzon_vg
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)  :: region
    !
    ! !REVISION HISTORY: 
    !  
    ! Programmed by Frank Dentener, August 1995;
    ! Ad Jeuken, KNMI, November 1998
    ! Modifications Bas Henzing, KNMI, 2002
    ! Adapted for TM5, Frank Dentener, JRC, 2002
    ! Paralel, Maarten Krol, Jul 2003
    !   29 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/LSPSCAV'

    real     ::  dt_lagrangian

    !  liquid water content of raining cloud
    !  rgas (8.314 J/mol/K) ---> 0.08314 atm/(mol/l)/K
    !  (thesis Frank Dentener, p. 31)
    !  1e-6 corresponds to 1 g/m3 dimensionless 
    real,parameter     ::  rtl_fac=rgas/1e2*1e-6 

    ! assumed pH of rainwater
    real,parameter     ::  hplus = 1e-5

    ! assumed fraction of in-cloud interstitial aerosol:
    real,parameter     ::  interstitial_fraction = 0.3   

    ! --- local ------------------------------

    real,dimension(:,:,:,:), pointer            :: rm
#ifdef slopes
    real,dimension(:,:,:,:), pointer            :: rxm, rym, rzm
#ifdef secmom    
    real,dimension(:,:,:,:), pointer            :: rxxm, rxym, rxzm, ryym, ryzm, rzzm
#endif
#endif
    real,dimension(:,:,:), pointer              :: t,cc

    real               ::  rtl,f,f1,f2,vol1,vol2,vol3,ahelp1,ahelp2
    real               ::  incloud,belowcloud,newcov,c_old,corr_diff,fnchem
    integer            ::  n,iscav,i,j,k,itemp,nzone_v, nloc
    real               ::  ztr, dkso2, dkhso3,   factor, heff

    ! oldcov: cloud cover in layer above, to calculate below-cloud scaveging.
    real,dimension(:,:),allocatable    :: oldcov
    
    integer            :: i1, i2, j1, j2, status

    ! for wetdep global budget of tracer #1
    real              :: g_sum_wet
    real, allocatable :: l_sum_wet(:,:)
    !
    ! --- begin ------------------------------
    !

!>>> TvN
! Tests have been performed at various resolutions
! using mixing times of 3, 6, 9, 12 and 24 hours.
! Simulations at 1x1 degrees are best reproduced
! by increasing the mixing time with 3 hours at 3x2 degrees
! and with 6 hours at 6x4 degrees.
    dt_lagrangian = 3600.0 * 3.0 ! value at 1x1 degrees or higher resolution 
    if (dy == 2) dt_lagrangian = dt_lagrangian + 3600.0 * 3.0
    if (dy == 4) dt_lagrangian = dt_lagrangian + 3600.0 * 6.0
!<<< TvN

    CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

#ifdef with_budgets    
    allocate(l_sum_wet(i1:i2,j1:j2))
    l_sum_wet = 0.0
#endif    

    rm => mass_dat(region)%rm   ! paralel over tracers
#ifdef slopes
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm
#ifdef secmom
    rxxm => mass_dat(region)%rxxm
    rxym => mass_dat(region)%rxym
    rxzm => mass_dat(region)%rxzm
    ryym => mass_dat(region)%ryym
    ryzm => mass_dat(region)%ryzm
    rzzm => mass_dat(region)%rzzm
#endif
#endif
    t  => temper_dat(region)%data
    cc =>     cc_dat(region)%data

    !$OMP PARALLEL &
    !$OMP  default (none) &
#if defined (with_budgets)
    !$OMP  shared  (nzon_vg, buddep_dat) &
#endif
    !$OMP  shared  (region, ier, isr, jer, jsr, tracer_loc, henry, &
    !$OMP           tracer_active, fnchem, rm, rxm, rym, rzm, t, cc, &
    !$OMP           nchem, rloss1, rloss1_m7, rloss2, rloss3, ra,nscav_type, &
    !$OMP           nscav_index,tref,im ) &
    !$OMP  private (i, j, k, n, nloc, iscav, jss, jee, rtl, itemp, &
    !$OMP           corr_diff, belowcloud, incloud, f, f1, f2, newcov, &
    !$OMP           vol1, vol2, vol3, c_old, nzone_v, ztr, dkso2, dkhso3, &
    !$OMP           factor, heff, oldcov, l_sum_wet)

    fnchem=real(nchem/(2*tref(region)))
    !
    allocate(oldcov(i1:i2,j1:j2))


    do iscav=1,nscav
       !
       n=nscav_index(iscav)                ! tracer number in global count
       nloc = n                            ! tracer number in local count  
       oldcov=0.
       !
       ! assumption no stratiform precipitation above the maximum 
       ! level of convection
       !
       do k=lmax_conv,1,-1   ! top to bottom
          do j=j1,j2
             do i=i1,i2 
                !
                ! calculate, depending on solubility, scale factor relative 
                ! to 100% scavenging (of HNO3)
                !
                ! rtl - composite factor of liquid water content, rgas 
                !       and liquid water content
                rtl=rtl_fac*t(i,j,k)
                ! multiplaction with Henry constant gives phase factor
                itemp=nint(t(i,j,k)-float(ntlow))
                itemp = min(max(1,itemp),ntemp)   ! corrected CMK dec2004 problems on ECMWF
                !
                !corr_diff=sqrt(ra(ihno3)/ra(n)) 
                corr_diff=sqrt(xmhno3/ra(n))   
                ! 
                select case (nscav_type(iscav))
                case(0)
                   incloud    = 0.0
                   belowcloud = 0.0
                case(1)   ! 100% solubility
                   !AJS: note that old code used factor rtl ~ 1,
                   !AJS: computed with huge henry factor ~ 1e7 :
                   !  rtl = rtl*henry(n,itemp) / ( 1.0 + rtl*henry(n,itemp) )   ! near 1.0
                   !  incloud    = rloss1(reion)%d3(i,j,k) * rtl
                   incloud    = rloss1(region)%d3(i,j,k)
                   belowcloud = rloss2(region)%d3(i,j,k)*corr_diff
                case(2)   ! henry solubility assumed
                   rtl = rtl*henry(n,itemp) / ( 1.0 + rtl*henry(n,itemp) )
                   incloud    = rloss1(region)%d3(i,j,k)*rtl
                   belowcloud = rloss2(region)%d3(i,j,k)*rtl*corr_diff
                case(3)   ! bulk aerosol 
                   incloud    = rloss1(region)%d3(i,j,k)*(1.0 - interstitial_fraction)
                   !>>>TvN
                   ! Alternative would be to make the interstitial fraction for bulk aerosols
                   ! consistent the values used for the M7 modes, 
                   ! which are taken from Bourgeois and Bey (JGR, 2011)
                   ! and distinguish between warm, mixed and ice clouds
                   !<<<TvN
                   belowcloud = rloss3(region)%d3(i,j,k)
                case(4)   ! SO2
                   ztr=(1./t(i,j,k)-1./298)
                   dkso2 =1.7e-2*exp(2090.*ztr)        !so2<=>hso3m+hplus
                   dkhso3  = 6.6e-8*exp(1510.*ztr)     !hso3m<=>so3-- + hplus
                   factor = 1.0 + dkso2/hplus + (dkso2*dkhso3)/(hplus**2)
                   heff = factor*henry(n,itemp)
                   rtl = rtl*heff/ ( 1.0 + rtl*heff )
                   incloud    = rloss1(region)%d3(i,j,k)*rtl !
                   belowcloud = rloss2(region)%d3(i,j,k)*rtl*corr_diff
!>>>TvN
! The in-cloud scavenging coefficients are defined as the fraction of the tracer 
! in the cloudy part of the grid box that is embedded in the cloud liquid or ice water,
! i.e. the non-interstitial part.
! We distinguish between liquid, mixed and ice stratiform clouds (Stier et al., 2005),
! depending on the local temperature in the grid cell (Croft et al., ACP, 2010).
! The in-cloud scavenging coefficients depend on size and composition;
! revised values for the M7 modes were provided by Bourgeois and Bey (JGR, 2011).
! For mixed clouds, an alternative method was presented by Zhang et al. (ACP, 2012), 
! which uses a continuous temperature dependency.
! Note that these in-cloud scavenging coefficients account for both nucleation scavenging
! and impaction scavenging (Croft et al., ACP, 2009; 2010).
! Thus, the below-cloud scavenging rates should only account for
! the impaction scavenging by precipitation coming from clouds above the current level.
!
! Estimates for below-cloud scavenging coefficients can be derived 
! from Fig. 2 of Dana and Hales (AE, 1975).
! For estimating these values from the figure, I used aerodynamic radii of 
! 0.007, 0.07, and 0.7 micron as the boundaries of the M7 modes
! (corresponding to a particle density of about 1800 g/cm^3).
! As in Stier et al. (2005), we do not distinguish between soluble and insoluble modes.
! Thus, dry particle radii can be used for estimating the scavenging coefficients from the figure
! (see also the mode boundaries applied in Fig. 2 in Croft et al., 2009).
! I thus arrive at the following rough estimates for below-cloud mass scavenging coefficients
! for the nucleation, aitken, accumulation, and coarse modes: ~0.01, 0.002, 0.01, and 1 mm^-1.
! These numbers are close to the estimates derived earlier from the same figure
! by Elisabetta Vignati, which were previously used: 0.005, 0.002, 0.008, and 1 mm^-1.
!
! However, both sets of estimates based on Dana and Hales are substantially higher
! than the values presented by Croft et al. (2009).
! From the curves presented in their Fig. 2 for the standard Marshall-Palmer rain distribution,
! rough estimates of the mass scavenging coefficients for the four size modes can be obtained.
! My estimates are 0.002, 0.0002, 0.03, and 0.7 mm^-1. 
! Note that especially the value for the accumulation mode is very sensitive to the
! actual mean particle size, and hard to estimate from the figure.
! Since the mean droplet size of the Marshall-Palmer distribution depends on the rain intensity,
! these estimates are only valid for a rain rate of 1 mm/hr.
! For simplicity, we assume that the scavenging coefficients derived from the figure at 1 mm/hr
! can also be applied at other rain intensities.
! 
! In the new implementation particle masses and numbers are scavenged at different rates.
! Rough estimates of the number scavenging coefficients for the four size modes
! can be obtained from the same figure in Croft et al.
! My estimates are 0.02, 0.001, 0.0003, and 0.3 mm^-1.

! Ideally, the below-cloud mass/number scavenging coefficients should be calculated
! using look-up tables to describe the dependence on median radius and precipitation rate, 
! e.g. following the formulation/curves presented by Croft et al.
!
                case(5)   ! soluble nu
                   !belowcloud=0.1*rloss3(region)%d3(i,j,k) ! 0.1*0.05=0.005 mm^-1
#ifdef with_m7
                   if (n /= inus_n) then
                     belowcloud=0.5*rloss3(region)%d3(i,j,k) ! 0.5*0.004 = 0.002 mm^-2
                   else
#endif
                      belowcloud=5.*rloss3(region)%d3(i,j,k) ! 5.*0.004 = 0.02 mm^-2
#ifdef with_m7
                  endif
#endif
                   !incloud=0.0
                   incloud=0.06*rloss1_m7(region)%d3(i,j,k)
                case(6)   ! soluble ai
                   !belowcloud=0.04*rloss3(region)%d3(i,j,k) ! 0.04*0.05=0.002 mm^-1
#ifdef with_m7
                   if (n /= iais_n) then
                     belowcloud=0.05*rloss3(region)%d3(i,j,k) ! 0.05*0.004 = 0.0002 mm^-2
                   else
#endif
                      belowcloud=0.25*rloss3(region)%d3(i,j,k) ! 0.25*0.004 = 0.001 mm^-2
#ifdef with_m7
                  endif
#endif
                   !incloud=0.0
                   if (t(i,j,k).gt.273.15) then 
                      incloud=0.25*rloss1_m7(region)%d3(i,j,k)
                   else
                      incloud=0.06*rloss1_m7(region)%d3(i,j,k)
                   endif
                case(7)   ! soluble ac
                   !belowcloud=0.16*rloss3(region)%d3(i,j,k) ! 0.16*0.05 = 0.008 mm^-1 
#ifdef with_m7
                   if (n /= iacs_n) then
                     belowcloud=7.5*rloss3(region)%d3(i,j,k)  ! 7.5*0.004 = 0.03 mm^-1
                   else  
#endif
                      belowcloud=0.075*rloss3(region)%d3(i,j,k)  ! 0.075*0.004 = 0.0003 mm^-1
#ifdef with_m7
                  endif
#endif
                   !incloud=rloss1(region)%d3(i,j,k)
                   if (t(i,j,k).gt.273.15) then 
                      incloud=0.85*rloss1_m7(region)%d3(i,j,k)
                   else
                      incloud=0.06*rloss1_m7(region)%d3(i,j,k)
                   endif 
                case(8)   ! soluble co
                   !belowcloud=20.*rloss3(region)%d3(i,j,k) ! 20*0.05 = 1 mm^-1 
#ifdef with_m7
                   if (n /= icos_n) then
                     belowcloud=175.*rloss3(region)%d3(i,j,k) ! 175*0.004 = 0.7 mm^-1
                   else
#endif
                      belowcloud=75.*rloss3(region)%d3(i,j,k) ! 75*0.004 = 0.3 mm^-1
#ifdef with_m7
                  endif
#endif
                   !incloud=rloss1(region)%d3(i,j,k)
                   if (t(i,j,k).gt.273.15) then
                      incloud=0.99*rloss1_m7(region)%d3(i,j,k)
                   else if (t(i,j,k).gt.238.15) then
                      incloud=0.75*rloss1_m7(region)%d3(i,j,k)
                   else
                      incloud=0.06*rloss1_m7(region)%d3(i,j,k)
                   endif
                case(9)   ! insoluble ai
                   !belowcloud = 0.04*rloss3(region)%d3(i,j,k) 
#ifdef with_m7
                   if (n /= iaii_n) then
                     belowcloud=0.05*rloss3(region)%d3(i,j,k)
                   else
#endif
                      belowcloud=0.25*rloss3(region)%d3(i,j,k)
#ifdef with_m7
                   endif
#endif
                   !incloud=0.0
                   if (t(i,j,k).gt.273.15) then
                      incloud=0.2*rloss1_m7(region)%d3(i,j,k)
                   else
                      incloud=0.06*rloss1_m7(region)%d3(i,j,k)
                   endif
                   
                case(10)   ! insoluble ac
                   !belowcloud=0.16*rloss3(region)%d3(i,j,k) 
#ifdef with_m7
                   if (n /= iaci_n) then
                     belowcloud=7.5*rloss3(region)%d3(i,j,k)
                   else
#endif
                      belowcloud=0.075*rloss3(region)%d3(i,j,k)
#ifdef with_m7
                  endif
#endif
                   !incloud=0.0
                   if (t(i,j,k).gt.273.15) then
                      incloud=0.4*rloss1_m7(region)%d3(i,j,k)
                   else
                      incloud=0.06*rloss1_m7(region)%d3(i,j,k)
                   endif
                   
                case(11)   ! insoluble co
                   !belowcloud=20.*rloss3(region)%d3(i,j,k) 
#ifdef with_m7
                   if (n /= icoi_n) then
                     belowcloud=175.*rloss3(region)%d3(i,j,k)
                   else
#endif
                      belowcloud=75.*rloss3(region)%d3(i,j,k)
#ifdef with_m7
                  endif
#endif
                   !incloud=0.0
                   if (t(i,j,k).gt.238.15) then
                      incloud=0.4*rloss1_m7(region)%d3(i,j,k)
                   else
                      incloud=0.06*rloss1_m7(region)%d3(i,j,k)
                   endif 
                case default
                   incloud    = 0.0
                   belowcloud = 0.0
                end select

                !if(incloud > 1e-4) then
                !print *, i,j,k,names(n),rtl, rloss1(region)%d3(i,j,k), rtl_fac
                !end if
                !
                ! f1, f2 are the fractions of the tracermass that remain in the 
                ! gridbox after scavenging.
                !
                f1=exp(-dt_lagrangian*incloud)
                f2=exp(-dt_lagrangian*belowcloud)
                !CMK f1=exp(-fnchem*incloud)
                !CMK f2=exp(-fnchem*belowcloud)
                !
                ! A grid box can be divided into three volumes
                ! 1) incloud scavenging       (newcov)
                ! 2) below cloud scavenging   (oldcov-newcov)
                ! 3) no in-cloud scavenging and no below cloud 
                !    scavenging by precipitation (no removal)
                !
                newcov=cc(i,j,k)
                vol1 = newcov
!>>> TvN
! oldcov denotes the area fraction occupied by precipitating clouds above the current level.
! It is determined assuming maximum overlap of the cloudy fractions of the layers above (see below).
! The precipitation rate used for below-cloud scavenging correctly does not include the contribution 
! from the current level.
!<<<TvN
                vol2 = max(0.,oldcov(i,j)-newcov)
                vol3 = 1.-vol1-vol2
                !CMK f=f1*vol1+f2*vol2+vol3
                f=(f1*vol1+f2*vol2+vol3)**(fnchem/dt_lagrangian)
                c_old=rm(i,j,k,nloc)
                rm(i,j,k,nloc)=c_old*f
#ifdef slopes
                rxm(i,j,k,nloc)=rxm(i,j,k,nloc)*f
                rym(i,j,k,nloc)=rym(i,j,k,nloc)*f
                rzm(i,j,k,nloc)=rzm(i,j,k,nloc)*f
#ifdef secmom
                rxxm(i,j,k,nloc)=rxxm(i,j,k,nloc)*f
                rxym(i,j,k,nloc)=rxym(i,j,k,nloc)*f
                rxzm(i,j,k,nloc)=rxzm(i,j,k,nloc)*f
                ryym(i,j,k,nloc)=ryym(i,j,k,nloc)*f
                ryzm(i,j,k,nloc)=ryzm(i,j,k,nloc)*f
                rzzm(i,j,k,nloc)=rzzm(i,j,k,nloc)*f
#endif
#endif
#ifdef with_budgets
                ! _____budget
                nzone_v = nzon_vg(k)
                buddep_dat(region)%lsp(i,j,nzone_v,n)= &
                     buddep_dat(region)%lsp(i,j,nzone_v,n)+ &
                     (c_old-rm(i,j,k,nloc))/ra(n)*1000.   ! in moles
                if ( n == 1 ) l_sum_wet(i,j) = l_sum_wet(i,j) + &
                     (c_old-rm(i,j,k,nloc))   ! in kg
                ! _____budget
#endif

! oldcov is determined assuming maximum overlap of the fractions of precipitating clouds overhead:
                if ( rloss1(region)%d3(i,j,k) > 0.0 ) oldcov(i,j)=max(oldcov(i,j),cc(i,j,k)) 
             end do !i
          end do !j
       end do !k
    end do !iscav

    deallocate(oldcov)

    !$OMP END PARALLEL

#ifdef with_budgets

    call REDUCE( dgrid(region), l_sum_wet, 0, 'SUM', g_sum_wet, status)
    IF_NOTOK_RETURN(status=1)

    if ( isRoot ) sum_wetdep(region) = sum_wetdep(region) + g_sum_wet

    deallocate( l_sum_wet )
#endif

    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#ifdef secmom
    nullify(rxxm)
    nullify(rxym)
    nullify(rxzm)
    nullify(ryym)
    nullify(ryzm)
    nullify(rzzm)
#endif
#endif
    nullify(t)
    nullify(cc)

  END SUBROUTINE LSPSCAV
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CALCULATE_LSP_SCAV
  !
  ! !DESCRIPTION:
  ! 
  ! Calculate wet removal rates rloss1,rloss2,rloss3 (s-1) for 
  ! stratiform precipitation from cloud-top to ground, 
  ! distinguishing between below-cloud and in-cloud scavenging. 
  !
  ! Method: 
  !        adapted from GJ Roelofs and Lelieveld, JGR, October 1995
  !
  ! fills array "rloss" should be called once new precipitation fields 
  ! are available (MK: in trace_after_read)
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CALCULATE_LSP_SCAV( region ) 
    !
    ! !USES:
    !
    use binas,       only : grav, rgas, xmair
    use dims,        only : im,jm,lm,lmax_conv
    use meteodata,   only : temper_dat, lwc_dat, iwc_dat, cc_dat
    use meteodata,   only : lsp_dat
    use global_data, only : emis_data
    use meteodata,   only : phlb_dat
    use partools,    only : isroot
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    !
    ! !REVISION HISTORY: 
    !   Programmed by:  Frank Dentener, IMAU, 1996
    !   Ad Jeuken, KNMI, 1998
    !   Modification, Bas Henzing, KNMI, 2002.
    !   Adapted for TM5, August 2002, Frank Dentener, JRC.
    !   And finally for the new version, MK, IMAU, March 2003
    !   Parallel, MK Jul 2003
    !   25 Jan 2012 - P. Le Sager - adapted for mpi lat-lon domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    integer :: i1, i2, j1, j2
    real,dimension(:,:,:),pointer    :: lsp
    real,dimension(:,:,:),pointer    :: t, lwc, iwc, cc, phlb
    real,parameter                   :: max_lwc=2.e-3        ! kg/m3 
    !
    ! Microphysical parameters:
    !
    !  rdrad2:  square of raindroplet radius (20 microns)
    !  dghno3: in [cm2/s] 
    !  dgair:  in [cm2/s] 
    !  rol:  density of water in [kg/m3]
    !  roi:  density of ice in [kg/m3]
    !
    real,parameter :: rdrad2 = (2E-5)**2
    real,parameter :: dghno3 = 0.136
    real,parameter :: dgair  = 0.133
    real,parameter :: rol    = 1000.
    real,parameter :: roi    = 917.
    !
    ! quantity used in calculation of Sherwood number
    !
    ! bas henzing: bug repair  real,parameter :: znsc=dgair/dghno3**(-3)  
    ! bas henzing: should be: znsc=(dgair/dghno3)**(1./3.); 
    !    znsc is now defined a real
    !
    real,dimension(:),allocatable   :: dzk
    real :: rflx,beta,fac,beta1,beta2,beta3,rlwc,rdrad,rn,ru
!>>>TvN
    real :: fac_m7, beta1_m7
!<<<
    real :: press,aird,dgairx,dghno3x
    real :: znre,znsh,znsc,zkg,totw,sfz,sfz_no
    !
    integer :: nfz,i,j,k,l,no_fz
    integer,dimension(:,:),allocatable :: kcltop
    real,dimension(:,:),allocatable    :: oldcov,fz

    real,dimension(:,:,:),allocatable :: precip ! precipitation per level (kg/m2/s)
    real,dimension(:,:,:),allocatable :: prf ! precipitation formation rate.
    ! evaporation NOT YET AVAILABLE
    !
    ! how much less efficient is tracer scavenged from ice
    ! cloud droplet compared to water cloud droplet. 
    ! This should be tracer dependent. 
    !
    real,parameter  :: ice_eff=0.2
    real            :: inc_rdf
    real,parameter  :: scale_heigth= rgas/grav/xmair*1e3 ! about 29.3*T (m)

    ! ---------------- START ------------------------------------------------
    t    => temper_dat(region)%data
    lwc  =>    lwc_dat(region)%data   !mk: lm, and not lmax_conv
    iwc  =>    iwc_dat(region)%data   !mk: lm, and not lmax_conv
    cc   =>     cc_dat(region)%data     !mk: lm, and not lmax_conv
    phlb =>   phlb_dat(region)%data  !mk: 1:lm+1
    lsp  =>    lsp_dat(region)%data

    CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )    

    allocate( oldcov( i1:i2, j1:j2           ) )
    allocate(     fz( i1:i2, j1:j2           ) )
    allocate( precip( i1:i2, j1:j2, lmax_conv) )
    allocate(    dzk(               lmax_conv) )
    allocate( kcltop( i1:i2, j1:j2           ) )
    allocate(    prf( i1:i2, j1:j2, lmax_conv) )
    !
    ! calculate precipitation formation rate  prf
    !
    call calfk
    !
    ! initialize cloud top
    !
    kcltop(:,:)=lmax_conv
    !
    ! calculate and rescale precip
    !
    sfz=0.
    nfz=0     
    precip(:,:,:)=0.
    !
    do j= j1, j2
       do i= i1, i2
          !
          ! Calculate precipitation intensity at the bottom of each layer
          !
          do k=1,lmax_conv-1
             ! thickness of layer, only correct in troposphere
             dzk(k)=scale_heigth*t(i,j,k)*alog(phlb(i,j,k)/(1.+phlb(i,j,k+1)))
          end do
          !
          do k=lmax_conv-1,1,-1
             precip(i,j,k)=precip(i,j,k+1)+prf(i,j,k)*dzk(k)    !precip: kg/m2/s
          end do
          !
          ! Rescale prf and precip such that these are consistent with ground lsp 
          ! Note that lsp is defined as the total large-scale precipitation (rain+snow)
          ! produced by the cloud scheme.
          ! prf, precip and lsp thus all account for snow/ice 
          !
          no_fz = 0   ! for statistics   !CMK was not initialised!
          sfz_no = 0.0

          fz(i,j)=0.
          !cmk if (lsp(i,j) > 1.e-5) then  old data came in mm/day
          if (lsp(i,j,1)*1e3*3600.*24. > 1.e-5) then   !new data are in m/s
             if (precip(i,j,1) > 0.) then
                fz(i,j)=lsp(i,j,1)*1e3/precip(i,j,1) ! m/s ->kg/m2/s   !new unit...
                !cmk fz(i,j)=lsp(i,j,1)/3600./24./precip(i,j,1) ! mm/day->kg/m2/s
                nfz=nfz+1
                ! precipitation statistics 
                ! (avoid 'strange' statistics by only few values)
                sfz=sfz+min(10.,fz(i,j))
             else
                ! no precipitation calculated, but ECMWF fields contain rain value
                no_fz=no_fz+1
                sfz_no=sfz_no+lsp(i,j,1)*1e3*86400.   !  (in mm/day)
             end if
          else
             precip(i,j,:)=0.
          end if

          do k=1,lmax_conv
             precip(i,j,k)=precip(i,j,k)*fz(i,j)
             prf(i,j,k)=prf(i,j,k)*fz(i,j)
          end do !k

       end do ! i
    end do ! j
    !
    !if(myid == root) then
    !   write(6,*) 'calculate_lsp_scav: region',region
    !   write(6,*) '   rainout: average scale factor for precipitation   = ',sfz/nfz
    !   write(6,*) '   rainout: percentage of columns with precipitation = ', &
    !        100.*nfz/real(im(region)*jm(region)),' %'
    !   if ( no_fz > 0 ) write(6,*) 'rainout: lsp and prf not consistent ', &
    !        no_fz,'average rainfall [mm/day]',sfz_no/no_fz
    !end if
    !
    ! initialise
    !
    oldcov=0.
    ! evap=0.
    rloss1(region)%d3=0.
!>>>TvN
    rloss1_m7(region)%d3=0.
!<<<TvN
    rloss2(region)%d3=0.
    rloss3(region)%d3=0.
    !
    do k=lmax_conv-1,1,-1  ! top-down
       do j=j1,j2
          do i=i1,i2
             !
             ! pressure correction for diffusion coefficient
             !
             press   = (phlb(i,j,k)+phlb(i,j,k+1))/2.
             dgairx  = dgair*1e5/press ! dgair at 1 atmosphere
             dghno3x = dghno3*1e5/press
             beta1=0.
             beta1_m7=0.
             beta2=0.
             beta3=0.
             !
             ! total cloudwater  (kg H2O / kg air)
             !
             totw=lwc(i,j,k)+iwc(i,j,k)
             !
             ! no influx set kcltop to low number  
             !
             if (precip(i,j,k+1)<=1.e-15) kcltop(i,j)=0 
             !
             !==========================================
             !  below-cloud scavenging (beta2 and beta3)
             !==========================================
             !
             ! with evaporation do:
             !    precip(i,j,k+1)=precip(i,j,k+1)-evap(i,j,k))>1E-15
             !
             if( (precip(i,j,k+1)>1e-15) .and. (k<kcltop(i,j)) &
                  .and. (oldcov(i,j)>0.) )then
                !
                ! rflx  in [mm/hr]   !MK?   thought precip was in kg/m2/s ??
                ! rlwc  in [vol mixing ratio]
                ! rdrad in [cm]
                ! ru    in [cm/s] (terminal velocity)
                ! znre = Reynolds number
                ! znsc = (Sherwood number)**(1./3.)
                ! zkg in [cm/s] = dghno3[cm^2/s]/rdrad[cm]
                !

                rflx  = precip(i,j,k+1)/oldcov(i,j)*3600.
                rlwc  = 72.*rflx**0.88*1.e-9
                rdrad = 0.1*0.3659*rflx**0.21
                !*******************************************
                ! correction by Twan van Noije, Bas Henzing:
                !    ru    = 9.58*(1.-exp(-(rdrad*10./0.885)**1.147))
                ! the above equation gives an approximation to the terminal velocity in m/s !!
                ! a conversion factor to cm/s should therefore be included
                !    znre  = 20.*rdrad*ru/dgair
                ! with ru in cm/s, the above expression overestimates the Reynolds number by a factor 10
                ! the combined effect of having ru in m/s and the above expression for the Reynolds number
                ! is to underestimate the Reynolds number by a factor 10, resulting in an underestimation
                ! of the Sherwood number and an overestimation of the below-cloud scavenging
                ! Twan van Noije/Bas Henzing, 24-02-2004
                !*******************************************
                ru    = 100.*9.58*(1.-exp(-(rdrad*10./0.885)**1.147))
                ! see Seinfeld (1986), p. 632
                znre  = 2.*rdrad*ru/dgair
                !WRONG  ru    = 9.58*(1.-exp(-(rdrad*10./0.885)**1.147)) 
                !WRONG  znre  = 20.*rdrad*ru/dgair
                znsc  = (dgair/dghno3)**(1./3.)
                znsh  = 1.+0.3*(znre**0.5)*znsc
                zkg   = dghno3/rdrad*znsh
                beta2 = 3.*zkg*rlwc/rdrad
                !
                ! beta3 for below cloud scavenging of accumulation range aerosol 
                ! (Dana and Hales, Atmos. Env. 1976, pp. 45-50
                ! assuming a Whitby aerosol distrbution r=0.034 sigma=2; 
                ! mass-mean radius r=0.14 microm; 
                ! figure 2 gives a wash-out coefficient of 0.05 mm^-1 (raindepth)
                ! for other sigma and r look-up tables can be generated
!>>>TvN
! A mass washout coefficient of 0.05 mm^-1 in Fig. 2 of Dana and Hales (1975)
! would correspond to a geometric mean/median particle radius of 0.14 micron.
! At a median radius of 0.034 micron and sigma=2, the value becomes ~0.002 mm^-1.
! However, as the aerodynamic radius is used in the plot, these values only apply
! to unit-density particles, with a density equal to that of water (1 g/cm^3).
! The aerodynamic radius equals the actual radius times 
! the square root of ratio of the particle density over that for water.
! For the bulk inorganic aerosols, the particle density is about 1800 g/cm^3,
! so the aerodynamic radius is 1.34 times the actual radius,
! and a median aerodynamic radius of ~0.046 micron should be used.
! This gives a mass washout coefficient of ~0.004 mm^-1.
! Thus, based on the work by Dana and Hales, the value 0.05 mm^-1 seems too low,
! and a value ~0.004 mm^-1 would be more appropriate.
                !
                ! mm-1*[mm hr-1] * [hr/s] => [s-1] 
                !
                !beta3= 0.05*rflx/3600.                          
                beta3= 0.004*rflx/3600.
                !
!<<<TvN
             end if
             !
             !  revaporation ( not implemented yet!, evap put to 0.)
             !
             !  IF ((1.-cc(i,j,k))>1E-20.AND.precip(i,j,k+1)>1E-20) THEN
             !    rev1=max(0.,EVAP(i,j,k)/precip(i,j,k+1))
             !      ! evaporation fraction    
             !    rev(i,j,k)=MIN(1.,rev1)     
             !  END IF
             !
             !==============================
             !  in cloud scavenging (beta1)
             !==============================
             !
             if (totw>1.e-9.and.prf(i,j,k)>0.and.cc(i,j,k)>0.02) then     
                !
                if(k>kcltop(i,j)) kcltop(i,j)=k     !set new cloud top
                !
                ! rlwc: convert mass mixing ratios to volume mixing ratios in cloud
                ! aird: air density in kg air / m^3
                ! max_lwc: in kg H2O / m^3
                ! prf: in kg H2O / m3 / s
                ! beta: in [s-1] = [vol mixing ratio]*[cm2/s]*1e-4/[m2]
                ! fac, beta1: in [s-1]
                !
                aird=press/t(i,j,k)/rgas*xmair*1.e-3
                rlwc=(lwc(i,j,k)/rol+iwc(i,j,k)/roi)*aird/cc(i,j,k)
                rlwc=min(max_lwc/rol,rlwc)
                !
                !bas henzing: bug repair: beta=3.*rlwc*(dghno3*1e-4)/rdrad2
                !bas henzing: should be:  
                !             beta=(3.*rlwc*(dghno3*1e-4)/(2.*rdrad2))*znsh
                !bas henzing: reference:  (Roelofs and Lelieveld, 1995)
                !fd               beta=(3.*rlwc*(dghno3*1e-4)/(2.*rdrad2))*znsh
                ! fd 13/08/2002 use old defenition again (pers. comm Henzing)
                beta=3.*rlwc*(dghno3*1e-4)/rdrad2
                !
                inc_rdf=(iwc(i,j,k)/totw)*ice_eff+lwc(i,j,k)/totw
                fac=prf(i,j,k)*inc_rdf/(totw*aird)
!>>>TvN
! In the calculation of fac (and thus fac_m7) currently the grid-cell average prf is used.
! Question: shouldn't we use the actual value in the cloudy part, i.e. prf/cc ?
! Note that scaling by 1/cc is also applied in the calculation of beta,
! and that for the below-cloud scavenging parameters beta2 and beta3 
! we also use the actual precipitation intensity in the precipitating fraction,
! i.e. precip/oldcov.
!<<<TvN
                !
                !resistance analogy: the slowest determines the overall resistance
                !
                beta1=1./(1./beta+1./fac)
                !  
!>>>TvN
! The scavengy efficiency for in-cloud scavenging of M7 aerosols by ice in stratiform clouds
! is now reduced by applying a lower scavenging coefficient for mixed and ice clouds
! according to Bourgeois and Bey (JGR, 2011).
                fac_m7=prf(i,j,k)/(totw*aird)
                beta1_m7=1./(1./beta+1./fac_m7)
!<<<TvN
                !
                !if no precipitation formation oldcov remains old value
                !
                oldcov(i,j)=max(oldcov(i,j),cc(i,j,k))           
                !
             end if ! in cloud scavenging
             !
             rloss1(region)%d3(i,j,k)=beta1 ! in cloud completely soluble
!>>>TvN
             rloss1_m7(region)%d3(i,j,k)=beta1_m7 ! as rloss1, but with ice as efficient as liquid water
!<<<TvN
             rloss2(region)%d3(i,j,k)=beta2 ! below cloud gas
             rloss3(region)%d3(i,j,k)=beta3 ! below cloud aerosol             
          end do !i
       end do !j
    end do !k 

    !if(myid == root) then
    !   write(*,*) 'calculate_lsp_scav: average rain-out loss rate 1 region', &
    !        region,sum(rloss1(region)%d3)/im(region)/jm(region)/lmax_conv
    !   write(*,*) 'calculate_lsp_scav: average rain-out loss rate 2 region', &
    !        region,sum(rloss2(region)%d3)/im(region)/jm(region)/lmax_conv
    !   write(*,*) 'calculate_lsp_scav: average rain-out loss rate 3 region', &
    !        region,sum(rloss3(region)%d3)/im(region)/jm(region)/lmax_conv
    !end if

    deallocate(oldcov)
    deallocate(fz)
    deallocate(precip)
    deallocate(prf)
    deallocate(dzk)
    deallocate(kcltop)

    nullify(lsp)
    nullify(t)
    nullify(lwc)
    nullify(iwc)
    nullify(cc)
    nullify(phlb)

  contains

    subroutine calfk
      !--------------------------------------------------------------
      ! calculate vertical distribution of large scale precipitation
      !               
      ! OUTPUT:  prf = precipitation formation rate (kg m-3 s-1)
      !
      ! Written by Ad Jeuken, KNMI, May 1998
      ! Adapted for TM5, MK, march 2003
      !--------------------------------------------------------------
      use toolbox, only: lvlpress
      !
      ! microphysical constants
      real,parameter :: cr1=1.e-4    ! s-1
      real,parameter :: cr2=1.0      ! m2 kg-1      
      real,parameter :: qcrs=3.e-4   ! kg/kg
      !bas henzing: replace alfa  real,parameter :: alfa=1.77, beta=0.16
      !bas henzing: new value alfa=3.29 (Heymsfield and Donner, 1990)
      real,parameter :: alfa=3.29, beta=0.16
      !bas henzing: replace b1   real,parameter :: b1=100., b2=0.5, Tberg=268.
      !bas henzing: new value b1=300. (Sunquist et al., 1989)
      real,parameter :: b1=300., b2=0.5, Tberg=268.
      !
      real                 :: plocal,f1,f2,delta_prec,ahelp,rain_frac,c0
      real                 :: pr0,qcr,qcl,qci,dzk,aird,press
      real                 :: qup,qdo,rup,rdo,vtup,vtdo,zcc,ccp,ccm
      integer              :: iclb,iclt,icldtop,i,j,k,l,l1,l2
      real,dimension(:),allocatable   :: zrho,pcl,pci
      !
      allocate(zrho(lmax_conv))
      allocate(pci(lmax_conv))
      allocate(pcl(lmax_conv))

      prf=0.
      do j=j1,j2
         do i=i1,i2
            !
            ! calculate airdensity "zrho" in kg/m3
            !
            do l=1,lmax_conv
               press=(phlb(i,j,l)+phlb(i,j,l+1))/2. 
               zrho(l)=press/t(i,j,l)/rgas*xmair*1.e-3
            end do
            !        
            iclb=0
            iclt=0
            !
            ! do not consider columns with lsp less than 1.e-5 mm/day
            !
            ! if (lsp(i,j,1)>1.e-5) then 
            if ( lsp(i,j,1)*1e3*3600.*24. > 1.e-5 ) then   !new data are in m/s
               ! determine cloud base
               k = 1
               do
                  if ( cc(i,j,k) >= 0.01 ) exit
                  if ( k == lmax_conv ) exit
                  k = k + 1
               end do
               iclb = k
               ! determine cloud top
               k = lmax_conv
               do
                  if ( cc(i,j,k) >= 0.01 ) exit
                  if ( k == 1 ) exit
                  k = k-1
               end do
               iclt = k
               ! check for inconsistency in cloud cover fields
               if ( iclb >= iclt ) iclb=iclt 
               if ( iclb < 1 ) iclb=1
               !mvw: replace fixed iclb/t=14 by 120 hPa level (about 15km)
               !           if (iclb>14) iclb=14  
               !           if (iclt>14) iclt=14
               icldtop=lvlpress(region,12000., 98400.)
               if ( iclb > icldtop ) iclb=icldtop
               if ( iclt > icldtop ) iclt=icldtop
               !
               pr0=0.
               pcl=0.
               pci=0.
               rain_frac=0.00001
               !
               ! start at cloudtop
               do k=iclt,iclb,-1
                  zcc=max(0.01,cc(i,j,k))
                  !
                  ! precipitation formation from ice clouds
                  !
                  ! pci in kg_H2O / (kg_air sec)
                  !
                  if ( ( t(i,j,k) < 258.0 ) .and. ( k > 1 ) ) then
                     ccp=max(0.01,cc(i,j,k+1))
                     ccm=max(0.01,cc(i,j,k-1))
                     qup=(iwc(i,j,k)/zcc+iwc(i,j,k+1)/ccp)*0.5
                     qdo=(iwc(i,j,k)/zcc+iwc(i,j,k-1)/ccm)*0.5
                     rup=(zrho(k)+zrho(k+1))*0.5
                     rdo=(zrho(k)+zrho(k-1))*0.5
                     vtup=alfa*(rup*qup)**beta
                     vtdo=alfa*(rdo*qdo)**beta
                     pci(k)=grav*(vtup*qup*rup-vtdo*qdo*rdo)/ &
                          (phlb(i,j,k+1)-phlb(i,j,k))
                     pci(k)=max(pci(k),0.)
                  end if
                  !
                  ! precipitation formation from liquid clouds
                  ! formulation ECMWF
                  !
                  qcl=lwc(i,j,k)/zcc
                  qcl=min(max_lwc/zrho(k),qcl)
                  qcl=max(0.001e-3/zrho(k),qcl)
                  !
                  ! pcl in kg_H2O / (kg_air sec)
                  !
                  plocal=pr0/rain_frac
                  f1=1.+b1*sqrt(plocal)
                  f2=1.+b2*sqrt(max(0.,Tberg-t(i,j,k)))
                  c0=cr1*f1*f2
                  qcr=qcrs/(f1*f2)
                  pcl(k)=zcc*c0*qcl*(1.-exp(-(qcl/qcr)**2))
                  !
                  ! prec at top next layer in kg m-2 s-1
                  !
                  dzk=29.3*t(i,j,k)*alog(phlb(i,j,k)/(1.+phlb(i,j,k+1)))
                  delta_prec=(pcl(k)+pci(k))*zrho(k)*dzk
                  ahelp=(zcc*delta_prec+rain_frac*pr0)/(delta_prec+pr0) 
                  rain_frac=max(rain_frac,ahelp)
                  pr0=pr0+(pcl(k)+pci(k))*zrho(k)*dzk
                  !
                  ! liquid+ice precipitation formation rates in  kg m-3 s-1
                  !
                  prf(i,j,k)= (pcl(k)+pci(k))*zrho(k)
                  !
               end do ! k
            end if ! lsp gt 1.e-5

         end do ! i
      end do ! j

      deallocate(zrho)
      deallocate(pci)
      deallocate(pcl)
      !
    end subroutine calfk

  END SUBROUTINE CALCULATE_LSP_SCAV
  !EOC

END MODULE WET_DEPOSITION