!MNH_LIC Copyright 2013-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
!      ##################################
       MODULE MODI_LIMA_PRECIP_SCAVENGING
!      ##################################
!
INTERFACE
   SUBROUTINE LIMA_PRECIP_SCAVENGING (HCLOUD, KLUOUT, KTCOUNT, PTSTEP,        &
                                      PRRT, PRHODREF, PRHODJ, PZZ,            & 
                                      PPABST, PTHT, PSVT, PRSVS, PINPAP )

use modd_nsv, only: nsv_lima_beg

CHARACTER(LEN=4),       INTENT(IN)    :: HCLOUD   ! cloud paramerization
INTEGER,                INTENT(IN)    :: KLUOUT   ! unit for output listing
INTEGER,                INTENT(IN)    :: KTCOUNT  ! iteration count
REAL,                   INTENT(IN)    :: PTSTEP   ! Double timestep except 
                                                  ! for the first time step
!
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRRT     ! Rain mixing ratio at t
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRHODREF ! Air Density [kg/m**3]
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ   ! Dry Density [kg]
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PZZ      ! Altitude
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PPABST   ! Absolute pressure at t
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PTHT     ! Theta at time t 
!
REAL, DIMENSION(:,:,:,NSV_LIMA_BEG:), INTENT(IN)    :: PSVT   ! Particle Concentration [kg-1]
REAL, DIMENSION(:,:,:,NSV_LIMA_BEG:), INTENT(INOUT) :: PRSVS  ! Total Number Scavenging Rate
!
REAL, DIMENSION(:,:),   INTENT(INOUT) :: PINPAP
!
END SUBROUTINE LIMA_PRECIP_SCAVENGING
END INTERFACE
END MODULE MODI_LIMA_PRECIP_SCAVENGING
!
!########################################################################
   SUBROUTINE LIMA_PRECIP_SCAVENGING (HCLOUD, KLUOUT, KTCOUNT, PTSTEP,        &
                                      PRRT, PRHODREF, PRHODJ, PZZ,            &
                                      PPABST, PTHT, PSVT, PRSVS, PINPAP )
!########################################################################x
!
!!    PURPOSE
!!    -------
!!      The purpose of this routine is to compute the total number
!!      below-cloud scavenging rate. 
!!
!!
!!**  METHOD
!!    ------
!!      We assume a generalized gamma distribution law for the raindrop.
!!      The aerosols particles distribution follows a log-normal law.
!!      First, we have to compute the Collision Efficiency, which takes 
!!      account of the three most important wet removal mechanism : 
!!      Brownian diffusion, interception and inertial impaction. 
!!      It is a function of several number (like Reynolds, Schmidt 
!!      or Stokes number for instance). Consequently,
!!      we need first to calculate these numbers. 
!!
!!      Then the scavenging coefficient is deduced from the integration 
!!      of the droplet size distribution, the falling velocity of 
!!      raindrop and aerosol, their diameter, and the collision 
!!      (or collection) efficiency, over the spectrum of droplet
!!      diameters.
!!
!!      The total scavenging rate of aerosol is computed from the 
!!      integration, over all the spectrum of particles aerosols 
!!      diameters, of the scavenging coefficient.
!! 
!!
!!    EXTERNAL
!!    --------
!!      Subroutine SCAV_MASS_SEDIMENTATION
!!
!!      Function COLL_EFFIC    : computes the collision efficiency
!!
!!      Function CONTJV   |
!!      Function GAUHER   |
!!      Function GAULAG   |-> in lima_functions.f90
!!      Function GAMMLN   |
!!     
!!
!!   REFERENCES
!!   ----------
!!   Seinfeld and Pandis
!!   Andronache
!!   
!!   AUTHOR
!!   ------
!!      J.-P. Pinty      * Laboratoire d'Aerologie*
!!      S.    Berthet    * Laboratoire d'Aerologie*
!!      B.    Vié        * Laboratoire d'Aerologie*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original             ??/??/13 
!!
!  P. Wautelet 28/05/2018: corrected truncated integer division (3/2 -> 1.5)
!  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
!  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
!  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
!  P. Wautelet 03/06/2020: bugfix: correct array starts for PSVT and PRSVS
!  P. Wautelet 11/02/2021: bugfix: ZRTMIN was of wrong size (replaced by a scalar)
!-------------------------------------------------------------------------------
!
!*                  0.DECLARATIONS          
!                   --------------
!
use modd_budget,          only: lbudget_sv, NBUDGET_SV1, tbudgets
USE MODD_CST
USE MODD_NSV
USE MODD_PARAMETERS
USE MODD_PARAM_LIMA,      ONLY: NMOD_IFN, NSPECIE, XFRAC,                         &
                                XMDIAM_IFN, XSIGMA_IFN, XRHO_IFN,                 &
                                NMOD_CCN, XR_MEAN_CCN, XLOGSIG_CCN, XRHO_CCN,     &
                                XALPHAR, XNUR,                                    &
                                LAERO_MASS, NDIAMR, NDIAMP, XT0SCAV, XTREF, XNDO, &
                                XMUA0, XT_SUTH_A, XMFPA0, XVISCW, XRHO00,         &
                                XRTMIN, XCTMIN
USE MODD_PARAM_LIMA_WARM, ONLY: XCR, XDR

use mode_budget,          only: Budget_store_init, Budget_store_end
use mode_tools,           only: Countjv

USE MODI_GAMMA
USE MODI_INI_NSV
USE MODI_LIMA_FUNCTIONS

IMPLICIT NONE
!
!*                 0.1 declarations of dummy arguments :
!
CHARACTER(LEN=4),       INTENT(IN)    :: HCLOUD   ! cloud paramerization
INTEGER,                INTENT(IN)    :: KLUOUT   ! unit for output listing
INTEGER,                INTENT(IN)    :: KTCOUNT  ! iteration count
REAL,                   INTENT(IN)    :: PTSTEP   ! Double timestep except 
                                                  ! for the first time step
!
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRRT     ! Rain mixing ratio at t
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRHODREF ! Air Density [kg/m**3]
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ   ! Dry Density [kg]
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PZZ      ! Altitude
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PPABST   ! Absolute pressure at t
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PTHT     ! Theta at time t 
!
REAL, DIMENSION(:,:,:,NSV_LIMA_BEG:), INTENT(IN)    :: PSVT   ! Particle Concentration [/m**3]
REAL, DIMENSION(:,:,:,NSV_LIMA_BEG:), INTENT(INOUT) :: PRSVS  ! Total Number Scavenging Rate
!
REAL, DIMENSION(:,:),   INTENT(INOUT) :: PINPAP
!
!*       0.2   Declarations of local variables :
!
INTEGER :: IIB           !  Define the domain where is 
INTEGER :: IIE           !  the microphysical sources have to be computed
INTEGER :: IJB           ! 
INTEGER :: IJE           !
INTEGER :: IKB           ! 
INTEGER :: IKE           !
!
INTEGER :: JSV               ! CCN or IFN mode 
INTEGER :: J1, J2, IJ, JMOD
!
LOGICAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3)) &
                                 :: GRAIN,  &! Test where rain is present
                                    GSCAV    ! Test where rain is present
INTEGER , DIMENSION(SIZE(GSCAV)) :: I1,I2,I3 ! Used to replace the COUNT
INTEGER                          :: JL       ! and PACK intrinsics
INTEGER                          :: ISCAV
!
REAL                     :: ZDENS_RATIO, & !density ratio 
                            ZNUM,        & !PNU-1.               
                            ZSHAPE_FACTOR
!
REAL,    DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3))  :: ZW     ! work array
REAL,    DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3))  :: PCRT   ! cloud droplet conc.
!
REAL, DIMENSION(:), ALLOCATABLE :: ZLAMBDAR,      &  !slope parameter of the 
                                                     ! generalized Gamma 
                                                     !distribution law for the 
                                                     !raindrop
                                   ZVISC_RATIO,   &  !viscosity ratio 
                                   ZMFPA,         &  !Mean Free Path
                                   ZRHODREF,      &  !Air Density [kg/m**3]
                                   ZVISCA,        &  !Viscosity of Air [kg/(m*s)]
                                   ZT,            &  !Absolute Temperature
                                   ZPABST,        &
                                   ZRRT,          &
                                   ZCONCP,        &
                                   ZCONCR,        &
                                   ZTOT_SCAV_RATE,&
                                   ZTOT_MASS_RATE,&
                                   ZMEAN_SCAV_COEF
!
REAL, DIMENSION(:,:), ALLOCATABLE :: &
                      ZVOLDR,        &  !Mean volumic Raindrop diameter [m]
                      ZBC_SCAV_COEF, &
                      ZCUNSLIP,      &  !CUnningham SLIP correction factor 
                      ZST_STAR,      &  !critical Stokes number for impaction
                      ZSC,           &  !aerosol particle Schmidt number
                      ZRE,           &  !raindrop Reynolds number (for radius)
                      ZFVELR,        &  !Falling VELocity of the Raindrop
                      ZRELT,         &  !RELaxation Time of the particle [s]
                      ZDIFF             !Particle Diffusivity
!
REAL, DIMENSION(NDIAMP) :: ZVOLDP,          & !Mean volumic diameter [m]
                           ZABSCISSP,       & !Aerosol Abscisses
                           ZWEIGHTP           !Aerosol Weights
REAL, DIMENSION(NDIAMR) :: ZABSCISSR,       & !Raindrop Abscisses
                           ZWEIGHTR           !Raindrop Weights
!
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZCOL_EF,     &! Collision efficiency
                                       ZSIZE_RATIO, &! Size Ratio
                                       ZST           ! Stokes number
!
REAL, DIMENSION(SIZE(PRRT,1),SIZE(PRRT,2),SIZE(PRRT,3)) :: ZRRS          
!
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3)) &
                                    :: PMEAN_SCAV_COEF, & !Mean Scavenging 
                                                          ! Coefficient
                                       PTOT_SCAV_RATE,  & !Total Number 
                                                          ! Scavenging Rate
                                       PTOT_MASS_RATE     !Total Mass
                                                          ! Scavenging Rate
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3),NDIAMP) &
                                    ::PBC_SCAV_COEF  !Scavenging Coefficient
REAL, DIMENSION(:), ALLOCATABLE :: ZKNUDSEN ! Knuudsen number
!
! Opt. BVIE
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3))        &
                   :: ZT_3D, ZCONCR_3D, ZVISCA_3D, ZMFPA_3D,               &
                      ZVISC_RATIO_3D, ZLAMBDAR_3D, FACTOR_3D
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3),NDIAMP) &
                   :: ZVOLDR_3D, ZVOLDR_3D_INV, ZVOLDR_3D_POW,             &
                      ZFVELR_3D, ZRE_3D, ZRE_3D_SQRT, ZST_STAR_3D
REAL, DIMENSION(:), ALLOCATABLE   :: FACTOR 
REAL, DIMENSION(:,:), ALLOCATABLE ::     &
                      ZRE_SQRT,          &  ! SQRT of raindrop Reynolds number
                      ZRE_INV,           &  ! INV of raindrop Reynolds number
                      ZSC_INV,           &  ! INV of aerosol particle Schmidt number
                      ZSC_SQRT,          &  ! SQRT of aerosol particle Schmidt number
                      ZSC_3SQRT,         &  ! aerosol particle Schmidt number**(1./3.)
                      ZVOLDR_POW,        &  ! Mean volumic Raindrop diameter [m] **(2+ZDR)
                      ZVOLDR_INV            ! INV of Mean volumic Raindrop diameter [m]
REAL               :: ZDENS_RATIO_SQRT 
INTEGER :: SV_VAR, NM, JM
integer :: idx
REAL :: XMDIAMP 
REAL :: XSIGMAP  
REAL :: XRHOP   
REAL :: XFRACP
!
!
!
!------------------------------------------------------------------------------

if ( lbudget_sv ) then
  do jl = 1, nmod_ccn
    idx = nsv_lima_ccn_free - 1 + jl
    call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + idx), 'SCAV', prsvs(:, :, :, idx) )
  end do
  do jl = 1, nmod_ifn
    idx = nsv_lima_ifn_free - 1 + jl
    call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + idx), 'SCAV', prsvs(:, :, :, idx) )
  end do
end if
!
!*       1.     PRELIMINARY COMPUTATIONS
!   	        ------------------------
!
!
IIB=1+JPHEXT
IIE=SIZE(PRHODREF,1) - JPHEXT
IJB=1+JPHEXT
IJE=SIZE(PRHODREF,2) - JPHEXT
IKB=1+JPVEXT
IKE=SIZE(PRHODREF,3) - JPVEXT
!
! PCRT
PCRT(:,:,:)=PSVT(:,:,:,NSV_LIMA_NR)
!
! Rain mask 
GRAIN(:,:,:) = .FALSE.
GRAIN(IIB:IIE,IJB:IJE,IKB:IKE) = (PRRT(IIB:IIE,IJB:IJE,IKB:IKE)>XRTMIN(3) &
                            .AND. PCRT(IIB:IIE,IJB:IJE,IKB:IKE)>XCTMIN(3) )
!
! Initialize the total mass scavenging rate if LAERO_MASS=T
IF (LAERO_MASS) PTOT_MASS_RATE(:,:,:)  = 0.
!
! Quadrature method: compute absissae and weights
CALL GAUHER(ZABSCISSP,ZWEIGHTP,NDIAMP)
ZNUM = XNUR-1.0E0
CALL GAULAG(ZABSCISSR,ZWEIGHTR,NDIAMR,ZNUM)
!
!
!------------------------------------------------------------------------------
!
!
!*       2.     NUMERICAL OPTIMIZATION
!   	        ----------------------
!
!
! Optimization : compute in advance parameters depending on rain particles and
! environment conditions only, to avoid multiple identical computations in loops
!
!
ZSHAPE_FACTOR = GAMMA_X0D(XNUR+3./XALPHAR)/GAMMA_X0D(XNUR) 
!
WHERE ( GRAIN(:,:,:) )
   !
   ZT_3D(:,:,:)      = PTHT(:,:,:) * ( PPABST(:,:,:)/XP00 )**(XRD/XCPD)
   ZCONCR_3D(:,:,:)  = PCRT(:,:,:) * PRHODREF(:,:,:)                    ![/m3]
   ! Sutherland law for viscosity of air
   ZVISCA_3D(:,:,:)  = XMUA0*(ZT_3D(:,:,:)/XTREF)**1.5*(XTREF+XT_SUTH_A) &
                                                        /(XT_SUTH_A+ZT_3D(:,:,:))
   ! Air mean free path
   ZMFPA_3D(:,:,:)   = XMFPA0*(XP00*ZT_3D(:,:,:))/(PPABST(:,:,:)*XT0SCAV)           
   ! Viscosity ratio
   ZVISC_RATIO_3D(:,:,:) = ZVISCA_3D(:,:,:)/XVISCW   !!!!! inversé par rapport à orig. !
   ! Rain drops parameters
   ZLAMBDAR_3D(:,:,:) = ( ((XPI/6.)*ZSHAPE_FACTOR*XRHOLW*ZCONCR_3D(:,:,:))  &
                          /(PRHODREF(:,:,:)*PRRT(:,:,:)) )**(1./3.)         ![/m]
   FACTOR_3D(:,:,:) = XPI*0.25*ZCONCR_3D(:,:,:)*XCR*(XRHO00/PRHODREF(:,:,:))**(0.4)
   !
END WHERE
!
DO J2=1,NDIAMR
   WHERE ( GRAIN(:,:,:) )
      ! exchange of variables: [m]
      ZVOLDR_3D(:,:,:,J2) = ZABSCISSR(J2)**(1./XALPHAR)/ZLAMBDAR_3D(:,:,:)
      ZVOLDR_3D_INV(:,:,:,J2) = 1./ZVOLDR_3D(:,:,:,J2)
      ZVOLDR_3D_POW(:,:,:,J2) = ZVOLDR_3D(:,:,:,J2)**(2.+XDR)
      ! Raindrop Falling VELocity [m/s]
      ZFVELR_3D(:,:,:,J2) = XCR*(ZVOLDR_3D(:,:,:,J2)**XDR)*(XRHO00/PRHODREF(:,:,:))**(0.4)
      ! Reynolds number
      ZRE_3D(:,:,:,J2) = ZVOLDR_3D(:,:,:,J2)*ZFVELR_3D(:,:,:,J2)          &
                                            *PRHODREF(:,:,:)/(2.0*ZVISCA_3D(:,:,:))   
      ZRE_3D_SQRT(:,:,:,J2) = SQRT( ZRE_3D(:,:,:,J2) )   
      ! Critical Stokes number
      ZST_STAR_3D(:,:,:,J2) = (1.2+(LOG(1.+ZRE_3D(:,:,:,J2)))/12.)        &
                                           /(1.+LOG(1.+ZRE_3D(:,:,:,J2)))
   END WHERE
END DO
!
!
!------------------------------------------------------------------------------
!
!
!*       3.     AEROSOL SCAVENGING
!   	        ------------------
!
!
! Iteration over the aerosol type and mode
!
DO JSV = 1, NMOD_CCN+NMOD_IFN
!
   IF (JSV .LE. NMOD_CCN) THEN
      JMOD = JSV
      SV_VAR = NSV_LIMA_CCN_FREE -1 + JMOD ! Variable number in PSVT
      NM = 1                               ! Number of species (for IFN int. mixing)
   ELSE
      JMOD = JSV - NMOD_CCN
      SV_VAR = NSV_LIMA_IFN_FREE -1 + JMOD
      NM = NSPECIE
   END IF
!
   PBC_SCAV_COEF(:,:,:,:) = 0. 
   PMEAN_SCAV_COEF(:,:,:) = 0.
   PTOT_SCAV_RATE(:,:,:)  = 0.
!
   GSCAV(:,:,:) = .FALSE.
   GSCAV(IIB:IIE,IJB:IJE,IKB:IKE) =GRAIN(IIB:IIE,IJB:IJE,IKB:IKE) .AND.        &
        (PSVT(IIB:IIE,IJB:IJE,IKB:IKE,SV_VAR)>1.0E-2)
   ISCAV = COUNTJV(GSCAV(:,:,:),I1(:),I2(:),I3(:))
!
   IF( ISCAV>=1 ) THEN
      ALLOCATE(ZVISC_RATIO(ISCAV))
      ALLOCATE(ZRHODREF(ISCAV))
      ALLOCATE(ZVISCA(ISCAV))
      ALLOCATE(ZT(ISCAV))
      ALLOCATE(ZRRT(ISCAV))
      ALLOCATE(ZCONCR(ISCAV))
      ALLOCATE(ZLAMBDAR(ISCAV))  
      ALLOCATE(ZCONCP(ISCAV))
      ALLOCATE(ZMFPA(ISCAV))
      ALLOCATE(ZTOT_SCAV_RATE(ISCAV))
      ALLOCATE(ZTOT_MASS_RATE(ISCAV))
      ALLOCATE(ZMEAN_SCAV_COEF(ISCAV))
      ALLOCATE(ZPABST(ISCAV))
      ALLOCATE(ZKNUDSEN(ISCAV))
      ALLOCATE(FACTOR(ISCAV))
!
      ALLOCATE(ZCUNSLIP(ISCAV,NDIAMP))
      ALLOCATE(ZBC_SCAV_COEF(ISCAV,NDIAMP))
      ALLOCATE(ZSC(ISCAV,NDIAMP))
      ALLOCATE(ZSC_INV(ISCAV,NDIAMP))
      ALLOCATE(ZSC_SQRT(ISCAV,NDIAMP))
      ALLOCATE(ZSC_3SQRT(ISCAV,NDIAMP))
      ALLOCATE(ZRELT(ISCAV,NDIAMP))
      ALLOCATE(ZDIFF(ISCAV,NDIAMP))
      ALLOCATE(ZVOLDR(ISCAV,NDIAMR))
      ALLOCATE(ZVOLDR_POW(ISCAV,NDIAMR))
      ALLOCATE(ZVOLDR_INV(ISCAV,NDIAMR))
      ALLOCATE(ZRE(ISCAV,NDIAMR))
      ALLOCATE(ZRE_INV(ISCAV,NDIAMR))
      ALLOCATE(ZRE_SQRT(ISCAV,NDIAMR))
      ALLOCATE(ZST_STAR(ISCAV,NDIAMR))
      ALLOCATE(ZFVELR(ISCAV,NDIAMR)) 
      ALLOCATE(ZST(ISCAV,NDIAMP,NDIAMR))
      ALLOCATE(ZCOL_EF(ISCAV,NDIAMP,NDIAMR))
      ALLOCATE(ZSIZE_RATIO(ISCAV,NDIAMP,NDIAMR))
!
      ZMEAN_SCAV_COEF(:)=0.
      ZTOT_SCAV_RATE(:) =0.
      ZTOT_MASS_RATE(:) =0.
      DO JL=1,ISCAV
         ZRHODREF(JL) =  PRHODREF(I1(JL),I2(JL),I3(JL))
         ZT(JL)       =     ZT_3D(I1(JL),I2(JL),I3(JL))
         ZRRT(JL)     =      PRRT(I1(JL),I2(JL),I3(JL))
         ZPABST(JL)   =    PPABST(I1(JL),I2(JL),I3(JL))
         ZCONCP(JL)   =      PSVT(I1(JL),I2(JL),I3(JL),SV_VAR)*ZRHODREF(JL)![/m3]
         ZCONCR(JL)   = ZCONCR_3D(I1(JL),I2(JL),I3(JL))                       ![/m3]
         ZVISCA(JL)   = ZVISCA_3D(I1(JL),I2(JL),I3(JL))
         ZMFPA(JL)    =  ZMFPA_3D(I1(JL),I2(JL),I3(JL))
         ZVISC_RATIO(JL) = ZVISC_RATIO_3D(I1(JL),I2(JL),I3(JL))
         ZLAMBDAR(JL) = ZLAMBDAR_3D(I1(JL),I2(JL),I3(JL))
         FACTOR(JL)   = FACTOR_3D(I1(JL),I2(JL),I3(JL))
         ZVOLDR(JL,:) = ZVOLDR_3D(I1(JL),I2(JL),I3(JL),:)
         ZVOLDR_POW(JL,:) = ZVOLDR_3D_POW(I1(JL),I2(JL),I3(JL),:)
         ZVOLDR_INV(JL,:) = ZVOLDR_3D_INV(I1(JL),I2(JL),I3(JL),:)
         ZFVELR(JL,:) = ZFVELR_3D(I1(JL),I2(JL),I3(JL),:)
         ZRE(JL,:)    = ZRE_3D(I1(JL),I2(JL),I3(JL),:)
         ZRE_SQRT(JL,:) = ZRE_3D_SQRT(I1(JL),I2(JL),I3(JL),:)
         ZST_STAR(JL,:) = ZST_STAR_3D(I1(JL),I2(JL),I3(JL),:)
      ENDDO
      ZRE_INV(:,:) = 1./ZRE(:,:)

      IF (ANY(ZCONCR .eq. 0.)) print *, 'valeur nulle dans ZLAMBDAR !' 
      IF (ANY(ZLAMBDAR .eq. 0.)) print *, 'valeur nulle dans ZLAMBDAR !' 
!
!------------------------------------------------------------------------------------
!
! Loop over the different species (for IFN int. mixing)
!
      DO JM = 1, NM  ! species (DM1,DM2,BC,O) for IFN
         IF ( JSV .LE. NMOD_CCN ) THEN          ! CCN case
            XRHOP   = XRHO_CCN(JMOD)
            XMDIAMP = 2*XR_MEAN_CCN(JMOD)
            XSIGMAP = EXP(XLOGSIG_CCN(JMOD))
            XFRACP  = 1.0
         ELSE                                   ! IFN case
            XRHOP   = XRHO_IFN(JM)
            XMDIAMP = XMDIAM_IFN(JM)
            XSIGMAP = XSIGMA_IFN(JM)
            XFRACP  = XFRAC(JM,JMOD)
         END IF
      !-----------------------------------------------------------------------------
      ! Loop over the aerosols particles diameters (log normal distribution law) :
      ! 
         DO J1=1,NDIAMP                        
            ! exchange of variables: [m]
            ZVOLDP(J1) = XMDIAMP * EXP(ZABSCISSP(J1)*SQRT(2.)*LOG(XSIGMAP))
            ! Cunningham slip correction factor (1+alpha*Knudsen) 
            ZKNUDSEN(:) = MIN( 20.,ZVOLDP(J1)/ZMFPA(:) )
            ZCUNSLIP(:,J1) = 1.0+2.0/ZKNUDSEN(:)*(1.257+0.4*EXP(-0.55*ZKNUDSEN(:)))
            ! Diffusion coefficient
            ZDIFF(:,J1) = XBOLTZ*ZT(:)*ZCUNSLIP(:,J1)/(3.*XPI*ZVISCA(:)*ZVOLDP(J1))
            ! Schmidt number
            ZSC(:,J1)       = ZVISCA(:)/(ZRHODREF(:)*ZDIFF(:,J1))  
            ZSC_INV(:,J1)   = 1./ZSC(:,J1)
            ZSC_SQRT(:,J1)  = SQRT( ZSC(:,J1) ) 
            ZSC_3SQRT(:,J1) = ZSC(:,J1)**(1./3.)  
            ! Characteristic Time Required for reaching terminal velocity 
            ZRELT(:,J1) = (ZVOLDP(J1)**2)*ZCUNSLIP(:,J1)*XRHOP/(18.*ZVISCA(:))
            ! Density number
            ZDENS_RATIO = XRHOP/XRHOLW
            ZDENS_RATIO_SQRT = SQRT(ZDENS_RATIO)
            ! Initialisation
            ZBC_SCAV_COEF(:,J1)=0.
         !-------------------------------------------------------------------------
         ! Loop over the drops diameters (generalized Gamma distribution) :
         !
            DO J2=1,NDIAMR
               ! Stokes number
               ZST(:,J1,J2) = 2.*ZRELT(:,J1)*(ZFVELR(:,J2)-ZRELT(1,J1)*XG)          &
                                            *ZVOLDR_INV(:,J2)
               ! Size Ratio
               ZSIZE_RATIO(:,J1,J2) = ZVOLDP(J1)*ZVOLDR_INV(:,J2)
               ! Collision Efficiency
               ZCOL_EF(:,J1,J2) = COLL_EFFI(ZRE, ZRE_INV, ZRE_SQRT, ZSC, ZSC_INV,   &
                                       ZSC_SQRT, ZSC_3SQRT, ZST, ZST_STAR,          &
                                       ZSIZE_RATIO, ZVISC_RATIO, ZDENS_RATIO_SQRT) 
               ! Below-Cloud Scavenging Coefficient for a fixed ZVOLDP: [/s]
               ZBC_SCAV_COEF(:,J1) = ZBC_SCAV_COEF(:,J1) +                          &
                                     ZCOL_EF(:,J1,J2) * ZWEIGHTR(J2) * FACTOR(:) * ZVOLDR_POW(:,J2)
            END DO 
         ! End of the loop over the drops diameters
         !--------------------------------------------------------------------------

            ! Total NUMBER Scavenging Rate of aerosol [m**-3.s**-1]
            ZTOT_SCAV_RATE(:) = ZTOT_SCAV_RATE(:) -                                 &
                                ZWEIGHTP(J1)*XFRACP*ZCONCP(:)*ZBC_SCAV_COEF(:,J1)
            ! Total MASS Scavenging Rate of aerosol [kg.m**-3.s**-1]
            ZTOT_MASS_RATE(:) = ZTOT_MASS_RATE(:) +                                 &
                                ZWEIGHTP(J1)*XFRACP*ZCONCP(:)*ZBC_SCAV_COEF(:,J1)   &
                                *XPI/6.*XRHOP*(ZVOLDP(J1)**3)  
         END DO
      ! End of the loop over the drops diameters
      !--------------------------------------------------------------------------

         ! Total NUMBER Scavenging Rate of aerosol [m**-3.s**-1]
         PTOT_SCAV_RATE(:,:,:)=UNPACK(ZTOT_SCAV_RATE(:),MASK=GSCAV(:,:,:),FIELD=0.0)
         ! Free particles (CCN or IFN) [/s]:
         PRSVS(:,:,:,SV_VAR) = max(PRSVS(:,:,:,SV_VAR)+PTOT_SCAV_RATE(:,:,:)  &
                                         * PRHODJ(:,:,:)/PRHODREF(:,:,:) , 0.0 )
         ! Total MASS Scavenging Rate of aerosol which REACH THE FLOOR because of 
         ! rain sedimentation [kg.m**-3.s**-1]
         IF (LAERO_MASS)THEN
            PTOT_MASS_RATE(:,:,:) = PTOT_MASS_RATE(:,:,:) +                         &
                 UNPACK(ZTOT_MASS_RATE(:), MASK=GSCAV(:,:,:), FIELD=0.0)
            CALL SCAV_MASS_SEDIMENTATION( HCLOUD, PTSTEP, KTCOUNT, PZZ, PRHODJ,     &
                                      PRHODREF, PRRT, PSVT(:,:,:,NSV_LIMA_SCAVMASS),&
                                      PRSVS(:,:,:,NSV_LIMA_SCAVMASS), PINPAP        )
            PRSVS(:,:,:,NSV_LIMA_SCAVMASS)=PRSVS(:,:,:,NSV_LIMA_SCAVMASS) +         &
                             PTOT_MASS_RATE(:,:,:)*PRHODJ(:,:,:)/PRHODREF(:,:,:)
         END IF
      ENDDO
! End of the loop over the aerosol species
!--------------------------------------------------------------------------
!
!
!
      DEALLOCATE(FACTOR)
      DEALLOCATE(ZSC_INV)
      DEALLOCATE(ZSC_SQRT)
      DEALLOCATE(ZSC_3SQRT)
      DEALLOCATE(ZRE_INV)
      DEALLOCATE(ZRE_SQRT)
      DEALLOCATE(ZVOLDR_POW)
      DEALLOCATE(ZVOLDR_INV)
!
      DEALLOCATE(ZFVELR)
      DEALLOCATE(ZRE)
      DEALLOCATE(ZST_STAR)
      DEALLOCATE(ZST)
      DEALLOCATE(ZSIZE_RATIO)
      DEALLOCATE(ZCOL_EF)
      DEALLOCATE(ZVOLDR)
      DEALLOCATE(ZDIFF)
      DEALLOCATE(ZRELT)
      DEALLOCATE(ZSC)
      DEALLOCATE(ZCUNSLIP)
      DEALLOCATE(ZBC_SCAV_COEF)
!
      DEALLOCATE(ZTOT_SCAV_RATE)
      DEALLOCATE(ZTOT_MASS_RATE)
      DEALLOCATE(ZMEAN_SCAV_COEF)
!
      DEALLOCATE(ZRRT)
      DEALLOCATE(ZCONCR)
      DEALLOCATE(ZLAMBDAR)
      DEALLOCATE(ZCONCP)
      DEALLOCATE(ZVISC_RATIO)
      DEALLOCATE(ZRHODREF)
      DEALLOCATE(ZVISCA)
      DEALLOCATE(ZPABST)
      DEALLOCATE(ZKNUDSEN)
      DEALLOCATE(ZT)
      DEALLOCATE(ZMFPA)
   ENDIF
ENDDO
!
if ( lbudget_sv ) then
  do jl = 1, nmod_ccn
    idx = nsv_lima_ccn_free - 1 + jl
    call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + idx), 'SCAV', prsvs(:, :, :, idx) )
  end do
  do jl = 1, nmod_ifn
    idx = nsv_lima_ifn_free - 1 + jl
    call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + idx), 'SCAV', prsvs(:, :, :, idx) )
  end do
end if
!------------------------------------------------------------------------------
!
!
!*       3.     SUBROUTINE AND FUNCTION
!   	        -----------------------
!
!
CONTAINS
!
!------------------------------------------------------------------------------
!     ##########################################################################
      SUBROUTINE SCAV_MASS_SEDIMENTATION( HCLOUD, PTSTEP, KTCOUNT, PZZ, PRHODJ,&
                                PRHODREF, PRAIN, PSVT_MASS, PRSVS_MASS, PINPAP )
!     ##########################################################################
!
!!
!!    PURPOSE
!!    -------
!!      The purpose of this routine is to compute the total mass of aerosol 
!!    scavenged by precipitations
!!
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!      None
!!     
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      Module MODD_PARAMETERS
!!          JPHEXT       : Horizontal external points number
!!          JPVEXT       : Vertical external points number
!!      Module MODD_CONF :
!!          CCONF configuration of the model for the first time step
!!
!!    REFERENCE
!!    ---------
!!      Book1 of the documentation ( routine CH_AQUEOUS_SEDIMENTATION )
!!
!!    AUTHOR
!!    ------
!!      J.-P. Pinty      * Laboratoire d'Aerologie*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    22/07/07
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE MODD_PARAMETERS
USE MODD_CONF
!
USE MODD_PARAM_LIMA,      ONLY : XCEXVT, XRTMIN
USE MODD_PARAM_LIMA_WARM, ONLY : XBR, XDR, XFSEDRR
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
!
CHARACTER (LEN=4),        INTENT(IN)    :: HCLOUD  ! Cloud parameterization
REAL,                     INTENT(IN)    :: PTSTEP  ! Time step  
INTEGER,                  INTENT(IN)    :: KTCOUNT ! Current time step number
!
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ     ! Height (z)
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ  ! Dry Density [kg]
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODREF! Reference density
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRAIN   ! Rain water m.r. source
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PSVT_MASS  ! Precip. aerosols at t
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRSVS_MASS ! Precip. aerosols source
!
REAL, DIMENSION(:,:),     INTENT(INOUT) :: PINPAP
!
!*       0.2   Declarations of local variables :
!
INTEGER :: JJ, JK, JN, JRR                ! Loop indexes 
INTEGER :: IIB, IIE, IJB, IJE, IKB, IKE   ! Physical domain
!
REAL    :: ZTSPLITR      ! Small time step for rain sedimentation
REAL    :: ZTSTEP        ! Large time step for rain sedimentation
!
!
LOGICAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) &
                                :: GSEDIM   ! where to compute the SED processes
INTEGER :: ISEDIM 
INTEGER , DIMENSION(SIZE(GSEDIM)) :: I1,I2,I3 ! Used to replace the COUNT
INTEGER                           :: JL       ! and PACK intrinsics
!
!
REAL,    DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3))   &
                                :: ZW,    & ! work array
                                   ZWSED, & ! sedimentation fluxes
                                   ZZS      ! Rain water m.r. source
!
REAL, DIMENSION(:), ALLOCATABLE :: ZRRS,     & ! Rain water m.r. source
                                   ZRHODREF, & ! RHO Dry REFerence
                                   ZZW         ! Work array
!
REAL                            :: ZRTMIN3
!
!
REAL                            :: ZVTRMAX, ZDZMIN, ZT
REAL,    SAVE                   :: ZEXSEDR
LOGICAL, SAVE                   :: GSFIRSTCALL = .TRUE.
INTEGER, SAVE                   :: ISPLITR
!
!-------------------------------------------------------------------------------
!
!*       1.     COMPUTE THE LOOP BOUNDS
!   	        -----------------------
!
IIB=1+JPHEXT
IIE=SIZE(PZZ,1) - JPHEXT
IJB=1+JPHEXT
IJE=SIZE(PZZ,2) - JPHEXT
IKB=1+JPVEXT
IKE=SIZE(PZZ,3) - JPVEXT
!
!-------------------------------------------------------------------------------
!
!*       2.     COMPUTE THE SEDIMENTATION (RS) SOURCE
!	        -------------------------------------
!
!*       2.1    splitting factor for high Courant number C=v_fall*(del_Z/del_T)
!  
firstcall : IF (GSFIRSTCALL) THEN
   GSFIRSTCALL = .FALSE.
   ZVTRMAX = 10.                          
   ZDZMIN = MINVAL(PZZ(IIB:IIE,IJB:IJE,IKB+1:IKE+1)-PZZ(IIB:IIE,IJB:IJE,IKB:IKE))
   ISPLITR = 1
   SPLIT : DO
      ZT = 2.* PTSTEP / REAL(ISPLITR)
      IF ( ZT * ZVTRMAX / ZDZMIN .LT. 1.) EXIT SPLIT
      ISPLITR = ISPLITR + 1
   END DO SPLIT
!
   ZEXSEDR = (XBR+XDR+1.0)/(XBR+1.0) 
!
END IF firstcall
!
!*       2.2    time splitting loop initialization        
!
IF( (KTCOUNT==1) .AND. (CCONF=='START') ) THEN
  ZTSPLITR = PTSTEP / REAL(ISPLITR)       ! Small time step
  ZTSTEP   = PTSTEP                        ! Large time step
  ELSE
  ZTSPLITR= 2. * PTSTEP / REAL(ISPLITR)
  ZTSTEP  = 2. * PTSTEP
END IF
!
!*       2.3    compute the fluxes
! 
!  optimization by looking for locations where
!  the precipitating fields are larger than a minimal value only !!!
!
ZRTMIN3 = XRTMIN(3) / ZTSTEP
ZZS(:,:,:) = PRAIN(:,:,:)
DO JN = 1 , ISPLITR
   GSEDIM(:,:,:) = .FALSE.
   GSEDIM(IIB:IIE,IJB:IJE,IKB:IKE) = ZZS(IIB:IIE,IJB:IJE,IKB:IKE) > ZRTMIN3
! 
   ISEDIM = COUNTJV( GSEDIM(:,:,:),I1(:),I2(:),I3(:))
   IF( ISEDIM >= 1 ) THEN
      IF( JN==1 ) THEN
         ZZS(:,:,:) = ZZS(:,:,:) * ZTSTEP
         DO JK = IKB , IKE-1
            ZW(:,:,JK) =ZTSPLITR*2./(PRHODREF(:,:,JK)*(PZZ(:,:,JK+2)-PZZ(:,:,JK)))
         END DO
         ZW(:,:,IKE)  =ZTSPLITR/(PRHODREF(:,:,IKE)*(PZZ(:,:,IKE+1)-PZZ(:,:,IKE)))
      END IF
      ALLOCATE(ZRRS(ISEDIM)) 
      ALLOCATE(ZRHODREF(ISEDIM))
      DO JL=1,ISEDIM
         ZRRS(JL) = ZZS(I1(JL),I2(JL),I3(JL))
         ZRHODREF(JL) =  PRHODREF(I1(JL),I2(JL),I3(JL))
      ENDDO
      ALLOCATE(ZZW(ISEDIM)) ; ZZW(:) = 0.0
!
!*       2.2.1   for rain
!
      ZZW(:) = XFSEDRR * ZRRS(:)**(ZEXSEDR) * ZRHODREF(:)**(ZEXSEDR-XCEXVT)
      ZWSED(:,:,:) = UNPACK( ZZW(:),MASK=GSEDIM(:,:,:),FIELD=0.0 )
      DO JK = IKB , IKE
         ZZS(:,:,JK) = ZZS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+1)-ZWSED(:,:,JK))
      END DO
      IF( JN==1 ) THEN
         PINPAP(:,:) = ZWSED(:,:,IKB)*                                            &
                             ( PSVT_MASS(:,:,IKB)/MAX(ZRTMIN3,PRRT(:,:,IKB)) )
      END IF
      DEALLOCATE(ZRHODREF)
      DEALLOCATE(ZRRS)
      DEALLOCATE(ZZW)
      IF( JN==ISPLITR ) THEN
         GSEDIM(:,:,:) = .FALSE.
         GSEDIM(IIB:IIE,IJB:IJE,IKB:IKE) = ZZS(IIB:IIE,IJB:IJE,IKB:IKE) > ZRTMIN3
         ZWSED(:,:,:) = 0.0
         WHERE( GSEDIM(:,:,:) ) 
            ZWSED(:,:,:) = 1.0/ZTSTEP - PRAIN(:,:,:)/ZZS(:,:,:)
         END WHERE
      END IF
   END IF
END DO
!
! Apply the rain sedimentation rate to the WR_xxx aqueous species
!
PRSVS_MASS(:,:,:) = PRSVS_MASS(:,:,:) + ZWSED(:,:,:)*PSVT_MASS(:,:,:)
!
END SUBROUTINE SCAV_MASS_SEDIMENTATION
!
!------------------------------------------------------------------------------
!
!###################################################################
  FUNCTION COLL_EFFI (PRE, PRE_INV, PRE_SQRT, PSC, PSC_INV, PSC_SQRT, &
                      PSC_3SQRT, PST, PST_STAR, PSIZE_RATIO,          &
                      PVISC_RATIO, PDENS_RATIO_SQRT) RESULT(PCOL_EF)
!###################################################################
!
!Compute the Raindrop-Aerosol Collision Efficiency
!
!*             0. DECLARATIONS  
!              ---------------
!
  IMPLICIT NONE
!
  INTEGER :: I
!
  REAL, DIMENSION(:,:), INTENT(IN)    :: PRE         
  REAL, DIMENSION(:,:), INTENT(IN)    :: PRE_INV         
  REAL, DIMENSION(:,:), INTENT(IN)    :: PRE_SQRT         
  REAL, DIMENSION(:,:), INTENT(IN)    :: PSC     
  REAL, DIMENSION(:,:), INTENT(IN)    :: PSC_INV    
  REAL, DIMENSION(:,:), INTENT(IN)    :: PSC_SQRT     
  REAL, DIMENSION(:,:), INTENT(IN)    :: PSC_3SQRT     
  REAL, DIMENSION(:,:), INTENT(IN)    :: PST_STAR   
!
  REAL, DIMENSION(:,:,:), INTENT(IN)  :: PST         
  REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSIZE_RATIO 
!
  REAL, DIMENSION(:), INTENT(IN)    :: PVISC_RATIO 
  REAL,               INTENT(IN)    :: PDENS_RATIO_SQRT 
!
  REAL, DIMENSION(SIZE(ZRE,1))      :: PCOL_EF   !result : collision efficiency
!
!-------------------------------------------------------------------------------
!
  PCOL_EF(:) = (4.*PSC_INV(:,J1)*PRE_INV(:,J2)*(1.+0.4*PRE_SQRT(:,J2)              &
                  *PSC_3SQRT(:,J1)+0.16*PRE_SQRT(:,J2)*PSC_SQRT(:,J1)))      &
              +(4.*PSIZE_RATIO(:,J1,J2)*(PVISC_RATIO(:)                    & 
              +(1.+2.*PRE_SQRT(:,J2))*PSIZE_RATIO(:,J1,J2)))  
  DO I=1,ISCAV
    IF (PST(I,J1,J2)>PST_STAR(I,J2)) THEN            
            PCOL_EF(I) = PCOL_EF(I)                                           &
                    +(PDENS_RATIO_SQRT*((PST(I,J1,J2)-PST_STAR(I,J2))        &
                    /(PST(I,J1,J2)-PST_STAR(I,J2)+2./3.))**(3./2.))
    ENDIF
  ENDDO
  END FUNCTION COLL_EFFI
!
!------------------------------------------------------------------------------
!
END SUBROUTINE LIMA_PRECIP_SCAVENGING
