source: trunk/LMDZ.TITAN/libf/muphytitan/mm_clouds.f90 @ 4084

! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne
! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
! email of the author : jeremie.burgalat@univ-reims.fr
!
! This software is a computer program whose purpose is to compute
! microphysics processes using a two-moments scheme.
!
! This library is governed by the CeCILL-B license under French law and
! abiding by the rules of distribution of free software.  You can  use,
! modify and/ or redistribute the software under the terms of the CeCILL-B
! license as circulated by CEA, CNRS and INRIA at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and  rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty  and the software's author,  the holder of the
! economic rights,  and the successive licensors  have only  limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading,  using,  modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean  that it is complicated to manipulate,  and  that  also
! therefore means  that it is reserved for developers  and  experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and,  more generally, to use and operate it in the
! same conditions as regards security.
!
! The fact that you are presently reading this means that you have had
! knowledge of the CeCILL-B license and that you accept its terms.

!! file: mm_clouds.f90
!! summary: Clouds microphysics module
!! author: J. Burgalat
!! date: 2013-2015,2017,2022-2023
!! modifications: B. de Batz de Trenquelléon

MODULE MM_CLOUDS
  !! Clouds microphysics module.
  !!
  !! The module contains all definitions of the microphysics processes related to clouds:
  !!
  !! - [nucleation](page/clouds.html#nucleation)
  !! - [condensation](page/clouds.html#condensation)
  !! - [sedimentation](page/clouds.html#sedimentation)
  !!
  !!
  !! The interface methods always use the global variables defined in [[mm_globals(module)]] when values
  !! (temperature, pressure, moments...) over the vertical grid are required.
  !! Consequently, all these functions only deal with output arguments which are most of the time the
  !! tendencies of relevant variables on the atmospheric column.
  !!
  !! @note
  !! Tendencies returned by public methods are always defined from  __TOP__ of the atmosphere to the
  !! __GROUND__.
  USE MM_MPREC
  USE MM_GLOBALS
  USE MM_METHODS
  IMPLICIT NONE

  PRIVATE

  PUBLIC :: mm_cloud_microphysics, mm_cloud_sedimentation, mm_cloud_nucond

CONTAINS

  !============================================================================
  ! CLOUDS MICROPHYSICS INTERFACE SUBROUTINE
  !============================================================================

  SUBROUTINE mm_cloud_microphysics(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs)
    !! Get the evolution of moments tracers through clouds microphysics processes.
    !!
    !! The subroutine is a wrapper to the clouds microphysics methods. It computes the tendencies of moments
    !! tracers for nucleation, condensation and sedimentation processes for the atmospheric column.
    !!
    !! @note
    !! Both __dm3i__ and __dgazs__ are 2D-array with the vertical layers in first dimension and the number
    !! of ice components in the second.
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0a
    !! Tendency of the 0th order moment of the aerosols (fractal mode) (\(m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3a
    !! Tendency of the 3rd order moment of the aerosols distribution (fractal mode) (\(m^{3}.m^{-3}\)) .
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0n
    !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3n
    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dm3i
    !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dgazs
    !! Tendencies of each condensible gaz species (\(mol.mol^{-1}\)).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE   :: zdm0n,zdm3n
    REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm3i
    INTEGER                                    :: i
    dm0a = 0._mm_wp ; dm3a = 0._mm_wp
    dm0n = 0._mm_wp ; dm3n = 0._mm_wp
    dm3i = 0._mm_wp ; dgazs = 0._mm_wp

    IF (mm_w_cloud_nucond) THEN
      ! Calls condensation/nucleation (and update saturation ratio diagnostic)
      ! ADDED : Extraction of nucleation and growth rates
      call mm_cloud_nucond(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs,mm_gazs_sat,mm_nrate,mm_grate)
    ENDIF

    IF (mm_w_cloud_sed) THEN
      ! Calls sedimentation
      ALLOCATE(zdm0n(mm_nla),zdm3n(mm_nla),zdm3i(mm_nla,mm_nesp))
      call mm_cloud_sedimentation(zdm0n,zdm3n,zdm3i)

      ! Computes settling velocity [m.s-1] of clouds (ccn and ices)
      mm_ccn_vsed(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)

      ! Computes flux [kg.m-2.s-1] and precipitation [kg.m-2] of ccn
      mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:))
      mm_ccn_prec = SUM(zdm3n*mm_dzlev*mm_rhoaer)

      ! Computes flux [kg.m-2.s-1] and precipitation [kg.m-2] of ices
      DO i = 1, mm_nesp
        mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho,(3._mm_wp*mm_m3ice(:,i))/(4._mm_wp*mm_pi))
        mm_ice_prec(i) = SUM(zdm3i(:,i)*mm_dzlev*mm_xESPS(i)%rho)
      ENDDO

      ! updates tendencies
      dm0n = dm0n + zdm0n
      dm3n = dm3n + zdm3n
      dm3i = dm3i + zdm3i
    ENDIF

  END SUBROUTINE mm_cloud_microphysics

  !-----------------------------------------------------------------------------
  ! NUCLEATION/CONDENSATION PROCESS RELATED METHODS
  !-----------------------------------------------------------------------------

  SUBROUTINE mm_cloud_nucond(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs,gazsat,nrate,grate)
    !! Get moments tendencies through nucleation/condensation/evaporation.
    !!
    !! The method is a wrapper of [[mm_clouds(module):nc_esp(subroutine)]] which computes the
    !! tendencies of tracers for all the condensible species given in the vector __xESPS__.
    !!
    !! Aerosols and CCN distribution evolution depends on the ice components:
    !!   - For nucleation only creation of CCN can occur.
    !!   - For condensation only loss of CCN can occur.
    !!
    !! We use the simple following rule to compute the variation of CCN and aerosols:
    !! The global variation of CCN (and thus aerosols) is determined from the most intense activity
    !! of the ice components.
    !!
    !! @warning
    !! __xESPS__, __m3i__ and __gazes__ must share the same indexing. For example if __xESPS(IDX)__
    !! corresponds to \(CH_{4}\) properties then __m3i(IDX)__ must be the total volume of solid
    !! \(CH_{4}\) (ice) and  __gazs(IDX)__ its vapor mole fraction.
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0a
    !! Tendency of the 0th order moment of the aerosols (fractal mode) (\(m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3a
    !! Tendency of the 3rd order moment of the aerosols distribution (fractal mode) (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0n
    !! Tendency of the 0th order moment of the aerosols distribution (fractal mode) (\(m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3n
    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dm3i
    !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dgazs
    !! Tendencies of each condensible gaz species (\(mol.mol^{-1}\)) .
    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: gazsat
    !! Saturation ratio of each condensible specie.

    ! ADDED : Extraction of the nucleation and growth rate
    REAL(kind=mm_wp), DIMENSION(:,:),INTENT(out) :: nrate
    !! Nucleation rate values of each condensible species (\(m^{-2}.s^{-1}\)).
    REAL(kind=mm_wp), DIMENSION(:,:),INTENT(out) :: grate
    !! Growth rate values of each condensible species (\(m^{2}.s^{-1}\)).

    INTEGER                                       :: i,idx
    REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm0a,zdm3a,zdm0n,zdm3n

    ALLOCATE(zdm0a(mm_nla,mm_nesp),zdm3a(mm_nla,mm_nesp), &
             zdm0n(mm_nla,mm_nesp),zdm3n(mm_nla,mm_nesp))
    zdm0a(:,:) = 0._mm_wp ; zdm3a(:,:) = 0._mm_wp
    zdm0n(:,:) = 0._mm_wp ; zdm3n(:,:) = 0._mm_wp
    DO i = 1, mm_nesp
      call nc_esp(mm_xESPS(i),mm_gazs(:,i),mm_m3ice(:,i),dgazs(:,i),dm3i(:,i), &
                  zdm0a(:,i),zdm3a(:,i),zdm0n(:,i),zdm3n(:,i),gazsat(:,i),nrate(:,i),grate(:,i))
      ! ADDED : Extraction of nucleation and growth rates
    ENDDO

    DO i=1, mm_nla
      ! retrieve the index of the maximum tendency of CCN where ice variation is not null.
      idx = MAXLOC(zdm0n(i,:),DIM=1,MASK=(dm3i(i,:) /= 0._mm_wp .OR. mm_m3ice(i,:)+dm3i(i,:) >= 0._mm_wp))
      IF (idx == 0) THEN
        dm0n(i) = 0._mm_wp
        dm3n(i) = 0._mm_wp
        dm0a(i) = 0._mm_wp
        dm3a(i) = 0._mm_wp
      ELSE
        IF (mm_debug .AND. ABS(zdm0n(i,idx)) > 1e3) THEN
          WRITE(*,'((a),I2.2,(a),ES10.3,(a))') "Z(",i,") = ",mm_play(i)*1e2, &
            " mbar: Max aer/ccn exchange variation due to specie: "//TRIM(mm_xESPS(idx)%name)
        ENDIF
        dm0n(i) = zdm0n(i,idx)
        dm3n(i) = zdm3n(i,idx)
        dm0a(i) = zdm0a(i,idx)
        dm3a(i) = zdm3a(i,idx)
      ENDIF
    ENDDO

  END SUBROUTINE mm_cloud_nucond

  SUBROUTINE nc_esp(xESP,vapX,m3iX,dvapX,dm3iX,dm0aer,dm3aer,dm0ccn,dm3ccn,Xsat,Xnrate,Xgrate) ! ADDED : arguments nrate, grate
    !! Get moments tendencies through nucleation/condensation/evaporation of a given condensible specie.
    !!
    !! The method computes the global tendencies of the aerosols, ccn and "ice" moments through cloud
    !! microphysics processes (nucleation & condensation).
    !!
    !! @warning
    !! Input quantities __m3iX__,__m3iO__, __m0aer__,__m3aer__, __m0ccn__,__m3ccn__ are assumed to be in
    !! \(X.m^{-3}\) (where X is the unit of the moment that is, a number for M0 and a volume - \(m^3\)
    !! for M3) ; __vapX__ must be expressed in term of molar fraction.
    TYPE(mm_esp), INTENT(in)                   :: xESP
    !! Condensate specie properties.
    REAL(kind=mm_wp),INTENT(in), DIMENSION(:)  :: vapX
    !! Gas specie molar fraction on the vertical grid from __TOP__ to __GROUND__ (\(mol.mol^{-1}\)).
    REAL(kind=mm_wp),INTENT(in), DIMENSION(:)  :: m3iX
    !! 3rd order moment of the ice component (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dvapX
    !! Tendency of gas specie (\(mol.mol^{-1}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3iX
    !! Tendency of the 3rd order moment of the ice component (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm0aer
    !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-3}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3aer
    !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm0ccn
    !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3ccn
    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: Xsat
    !! Saturation ratio values on the vertical grid (--).

    ! ADDED : Extraction of the nucleation and growth rate
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: Xnrate
    !! Nucleation rate values on the vertical grid for the species X (\(m^{-2}.s^{-1}\)).
    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: Xgrate
    !! Growth rate values on the vertical grid for the species X (\(m^{2}.s^{-1}\)).


    INTEGER                                 :: i
    REAL(kind=mm_wp)                        :: bef,aft
    REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: sm0a,sm3a,sm0n,sm3n,sm3iX
    REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: zm0a,zm3a,zm0n,zm3n,zm3iX,zvapX
    REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: pX,sig,qsat,seq,up,down,ctot,newvap,nucr,grate,cm0,cm3,drad

    ! get a copy of drop radius.
    drad(:) = mm_drad(:)

    ! Initialization :
    ! Copy input argument and convert units X.m-3 -> X.kg-1
    ! sXXX is the initial converted value saved
    sm3iX = m3iX/mm_rhoair
    sm0a = mm_m0aer_f/mm_rhoair ; sm3a = mm_m3aer_f/mm_rhoair
    sm0n = mm_m0ccn/mm_rhoair ; sm3n = mm_m3ccn/mm_rhoair
    ! zXXX is our working copy
    zm3iX = sm3iX ; zm0a = sm0a ; zm3a = sm3a ; zm0n = sm0n ; zm3n = sm3n

    ! Molar fraction of X specie is set in mass mixing ratio [kg.kg-1]
    zvapX  = vapX  * xESP%fmol2fmas
    ! Surface tension [N.m-1]
    sig = mm_sigX(mm_temp,xESP)
    ! X specie mass mixing ratio at saturation [kg.kg-1]
    qsat = mm_ysatX(mm_temp,mm_play,xESP) * xESP%fmol2fmas
    ! partial pressure of X specie
    pX = vapX * mm_play
    ! Saturation ratio
    Xsat = zvapX / qsat

    
    ! Gets nucleation rate (ccn radius is the monomer !)
    call nuc_rate((/(mm_rm, i=1,mm_nla)/),mm_temp,xESP,pX,Xsat,nucr)
    ! IMPORTANT: update CCN and aerosols moment from nucleation NOW !
    ! Doing so should prevent a nasty bug that occurs if we want to generate clouds from scratch (i.e. a "dry" atmosphere without any clouds tracers already present).
    !
    ! In such case, we do not produce ice variation on the first call of the method, at most only CCN are produced (i.e. dm3i == 0, dm3n != 0)
    ! But the rules for computing the global tendencies in mm_cloud_nucond state that the global variation for CCN is due to the most active specie exchange.
    !
    ! for nucleation we have the following equations:
    !   dMa(k)/dt = - dMn(k)/dt (conservation of aerosols+ccn)           (1)
    !   dMa(k)/dt = - 4*PI*nucr/rm * Ma(k+3)                             (2)
    !             = - 4*PI*nucr/rm * alpha(k+3)/alpha(k) * rc**3 * Ma(k)
    ! With:
    !   - Ma(k): k-th order moment of aerosol
    !   - Mn(k): k-th order moment of ccn
    !   - nucr : the nucleation rate.
    ! We solve (implicit scheme) :
    !   CST_M(k) = 4*PI*nucr/rm * alpha(k+3)/alpha(k)*rc**3 * dt
    !   Ma(k)[t+dt] = 1/(1+CST_M(k)) * Ma(k)[t]                          (3)
    ! Then, from eq. 2:
    ! Mn(k)[t+dt] = Mn(k)[t] + CST_M(k)/(1+CST_M(k))*Ma(k)[t]            (4)
    
    ! Copies the nucleation rate into an output variable
    Xnrate = nucr
    
    cm0 = 4._mm_wp*mm_pi*nucr/mm_rm*mm_alpha_f(3._mm_wp)*mm_rcf**3*mm_dt
    cm3 = 4._mm_wp*mm_pi*nucr/mm_rm*mm_alpha_f(6._mm_wp)/mm_alpha_f(3._mm_wp)*mm_rcf**3*mm_dt
    zm0a = 1._mm_wp/(1._mm_wp+cm0) * zm0a
    zm3a = 1._mm_wp/(1._mm_wp+cm3) * zm3a
    WHERE (zm0a <= 0._mm_wp .OR. zm3a <= 0._mm_wp)
      zm0a = 0._mm_wp
      zm3a = 0._mm_wp
      zm0n = zm0n + sm0a
      zm3n = zm3n + sm3a
    ELSEWHERE
      zm0n = zm0n + cm0*zm0a
      zm3n = zm3n + cm3*zm3a
    ENDWHERE

    ! update the drop radius (we probably should recompute totally the radius to be in better agreement with the other moments)
    ! We must manage the case where there is no ices and no ccn ==> drop radius is ZERO,
    ! but conditions are met to spawn nucleation process: creation of ccn.
    ! Then we set the drop radius to the monomer radius.
    !
    ! Doing so will prevent a nasty bug to occur later when ice volume is updated !
    WHERE (nucr > 0._mm_wp .AND. drad <= mm_drad_min)
      drad = mm_rm
    ENDWHERE

    ! Equilibrium saturation near the drop
    seq = exp(min(30._mm_wp,2._mm_wp*sig*xESP%masmol/(xESP%rho*mm_rgas*mm_temp*drad)))
    ! Gets growth rate
    call growth_rate(mm_temp,mm_play,pX/Xsat,xESP,seq,drad,grate)
    ctot = zvapx + xESP%rho * zm3iX
    up = zvapx + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad * seq * zm0n
    down = 1._mm_wp + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad / qsat * zm0n
    ! gets new vapor X specie mass mixing ratio : cannot be greater than the
    ! total gas + ice and lower than nothing :)
    newvap = max(min(up/down,ctot),0._mm_wp)
    ! gets "true" growth rate
    grate = grate * (newvap/qsat - seq)
    ! Copies the growth rate into an output variable
    Xgrate = grate

    ! computes tendencies through condensation
    DO i=1,mm_nla
      ! check for the specific case : NO ICE and SUBLIMATION
      IF (zm3iX(i) <= 0._mm_wp .AND. grate(i) <= 0._mm_wp) THEN
        zm3iX(i) = 0._mm_wp
      ELSE
        ! update ice volume ...
        zm3iX(i) = zm3iX(i) + mm_dt*grate(i)*4._mm_wp*mm_pi*drad(i)*zm0n(i)
        ! ... and check if there is ice left in the ccn
        IF (zm3iX(i) <= 0._mm_wp) THEN
          zm3iX(i) = 0._mm_wp
          zm0a(i) = zm0a(i) + zm0n(i) ; zm0n(i) = 0._mm_wp
          zm3a(i) = zm3a(i) + zm3n(i) ; zm3n(i) = 0._mm_wp
        ENDIF
      ENDIF
    ENDDO

    ! Computes balance
    IF (mm_debug) THEN
      DO i=1,mm_nla
        bef = sm0a(i) + sm0n(i)
        aft = zm0a(i)  + zm0n(i)
        IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
          WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
            "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M0 not conserved (z=",i,")",bef," <-> ",aft
        ENDIF
        bef = sm3a(i) + sm3n(i)
        aft = zm3a(i)  + zm3n(i)
        IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
          WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
            "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M3 not conserved (z=",i,")",bef," <-> ",aft
        ENDIF
      ENDDO
    ENDIF

    ! compute tendencies:
    ! all of these tendencies are in X.m-3 !
    dm0aer = (zm0a - sm0a)*mm_rhoair
    dm3aer = (zm3a - sm3a)*mm_rhoair
    dm0ccn = (zm0n - sm0n)*mm_rhoair
    dm3ccn = (zm3n - sm3n)*mm_rhoair

    dm3iX  = (zm3iX - sm3iX)                    ! this one in X.kg-1 (temporary) !
    dvapX  = -xESP%rho * dm3iX / xESP%fmol2fmas ! in order to compute this one in mol.mol-1
    dm3iX  = dm3iX*mm_rhoair                    ! update ice tendencies in X.m-3 !

  END SUBROUTINE nc_esp

  SUBROUTINE nuc_rate(rccn,temp,xESP,pvp,sat,rate)
    !! Get nucleation rate.
    !!
    !! The method computes the heterogeneous nucleation rate for the given specie on a fractal particle
    !! of size __rccn__.
    !! Except __xESP__, all arguments are vectors of the same size (vertical grid).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: rccn !! Radius of the cloud condensation nuclei (m).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: temp !! Temperature (K).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pvp  !! Partial vapor pressure of X specie (Pa).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: sat  !! Saturation ratio of given specie (--).
    TYPE(mm_esp), INTENT(in)                    :: xESP !! X specie properties (--).
    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate !! The nucleation rate (\(m^{-2}.s^{-1}\)).
    INTEGER          :: i
    REAL(kind=mm_wp) :: r,t,s,sig,nX,rstar,gstar,x,zeldov,deltaf,fsh,fstar

    rate(:) = 0._mm_wp
    ! Activation condition
    DO i=1, SIZE(rccn)
      s = sat(i)
      IF (s > 1._mm_wp) THEN
        t = temp(i) ; r = rccn(i)
        sig = mm_sigX(t,xESP)
        nX    = pvp(i)/mm_kboltz/t
        rstar = 2._mm_wp*sig*xESP%vol/(mm_kboltz*t*dlog(s))
        ! curvature radius
        x     = r/rstar
        fsh   = mm_fshape(xESP%mteta,x)
        fstar = (4._mm_wp*mm_pi/3._mm_wp)*sig*(rstar**2.)*fsh
        deltaf=MIN(MAX((2.*mm_fdes-mm_fdif-fstar)/(mm_kboltz*t),-100._mm_wp),100._mm_wp)
        IF (deltaf > -100._mm_wp) THEN
          gstar  = 4._mm_wp*mm_pi*(rstar**3)/(3._mm_wp*xESP%vol)
          zeldov = dsqrt(fstar/(3._mm_wp*mm_pi*mm_kboltz*t*(gstar**2)))
          rate(i)= zeldov*mm_kboltz*t*(nX*rstar)**2._mm_wp*dexp(deltaf)/(fsh*mm_nus*xESP%mas)
        ENDIF
      ENDIF
    ENDDO

    RETURN
  END SUBROUTINE nuc_rate

  SUBROUTINE growth_rate(temp,pres,pXsat,xESP,seq,drad,rate)
    !! Get growth rate through condensation/evaporation process.
    !!
    !! The method computes the growth rate a drop through condensation/evaporation processes:
    !!
    !! $$ r \times \frac{dr}{dt} = g_{rate} \times (S - S_{eq}) $$
    !!
    !! Except __xESP__ which is a scalar, all arguments are vectors of the same size (vertical grid).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: temp  !! Temperature (K).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pres  !! Pressure level (Pa).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pXsat !! Saturation vapor pressure of specie (Pa).
    TYPE(mm_esp), INTENT(in)                    :: xESP  !! Specie properties.
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: seq   !! Equilibrium saturation near the drop.
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: drad  !! Drop radius (m).
    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate  !! Growth rate (\(m^{2}.s^{-1}\)).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: k,knu,slf,rkc,rdc,l,dv
    INTEGER                                     :: n
    n = SIZE(temp)
    ALLOCATE(k(n),knu(n),slf(n),rkc(n),rdc(n),l(n),dv(n))

    ! N2 (air) Thermal conductivity (where does it come from ?)
    k(:) = ( 2.857e-2_mm_wp*temp-0.5428_mm_wp)*4.184e-3_mm_wp
    ! Gas mean free path
    l(:) = mm_lambda_g(temp,pres)
    ! Diffusion coefficient of X gas
    Dv(:) = 1._mm_wp/3._mm_wp*dsqrt(8._mm_wp*mm_rgas*temp(:)/(mm_pi*xESP%masmol))*mm_kboltz*temp(:) / &
      (mm_pi*pres(:)*(mm_air_rad+xESP%ray)**2* dsqrt(1._mm_wp+xESP%fmol2fmas))
    knu(:) = l(:)/drad(:)                                         ! The knudsen number of the drop
    slf(:) = (1.333_mm_wp+0.71_mm_wp/knu(:))/(1._mm_wp+1._mm_wp/knu(:)) ! Slip flow correction
    Dv(:)  = Dv(:)/(1._mm_wp+slf(:)*knu(:))
    ! latent heat resistance coefficient
    rkc(:) = mm_lheatX(temp(:),xESP)**2 * xESP%rho * xESP%masmol / (k(:)*mm_rgas*temp(:)**2)
    ! Diffusion resistance coefficient
    rdc(:) = mm_rgas * temp(:) * xESP%rho / (Dv(:)*pXsat(:)*xESP%masmol)
    ! Growth rate: rdr/dt = rate * (S-Seq) ; rate is returned
    rate(:) = 1._mm_wp / (seq(:)*rkc(:)+rdc(:))
    RETURN
  END SUBROUTINE growth_rate


  !-----------------------------------------------------------------------------
  ! SEDIMENTATION PROCESS RELATED METHODS
  !-----------------------------------------------------------------------------

  SUBROUTINE mm_cloud_sedimentation(dm0n,dm3n,dm3i)
    !! Compute the tendency of _clouds_ related moments through sedimentation process.
    !!
    !! The method computes the tendencies of moments related to cloud microphysics through
    !! sedimentation process. The algorithm used here differs from
    !! [[mm_haze(module):mm_haze_sedimentation(subroutine)]] as all moments settle with the same
    !! terminal velocity which is computed with the average drop radius of the size distribution.
    !! We simply compute an _exchange matrix_ that stores the new positions of each cells through
    !! sedimentation process and then computes the matrix
    !! product with input moments values to get final tendencies.
    REAL(kind=mm_wp), INTENT(out), DIMENSION(:)   :: dm0n
    !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)).
    REAL(kind=mm_wp), INTENT(out), DIMENSION(:)   :: dm3n
    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp), INTENT(out), DIMENSION(:,:) :: dm3i
    !! Tendencies of the 3rd order moment of each ice component of the cloud (\(m^{3}m^{-3}\)).
    INTEGER                                   :: im,nm
    REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: moms, momsf,chg_matrix

    nm = 2 + mm_nesp
    ALLOCATE(moms(mm_nla,nm),momsf(mm_nla,nm),chg_matrix(mm_nla,mm_nla))
    ! Initializes moms
    moms(:,1)  = mm_m0ccn * mm_dzlev
    moms(:,2)  = mm_m3ccn * mm_dzlev
    DO im=1,mm_nesp ; moms(:,2+im) = mm_m3ice(:,im) * mm_dzlev ; ENDDO
    ! Computes exchange matrix
    CALL exchange(mm_drad,mm_drho,mm_dt,chg_matrix)
    ! Computes final moments values
    DO im=1,nm ; momsf(:,im) = MATMUL(chg_matrix,moms(:,im)) ; ENDDO
    ! Computes tendencies (converted in X.m-3)
    dm0n = (momsf(:,1)-moms(:,1))/mm_dzlev
    dm3n = (momsf(:,2)-moms(:,2))/mm_dzlev
    DO im=1,mm_nesp ; dm3i(:,im) = (momsf(:,2+im)-moms(:,2+im))/mm_dzlev ; ENDDO
    RETURN
  END SUBROUTINE mm_cloud_sedimentation

  SUBROUTINE exchange(rad,rhog,dt,matrix)
    !! Compute the exchange matrix.
    !!
    !! The subroutine computes the matrix exchange used by
    !! [[mm_clouds(module):mm_cloud_sedimentation(subroutine)]] to compute moments tendencies
    !! through sedimentation process. Both __rad__ and __rhog__ must be vector with relevant
    !! values over the atmospheric vertical structure. __matrix__ is square 2D-array with same
    !! dimension size than __rad__.
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
    !! Cloud drop radius over the atmospheric vertical structure (m).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rhog
    !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)).
    REAL(kind=mm_wp), INTENT(in)               :: dt
    !! Timestep (s).
    REAL(kind=mm_wp), INTENT(out)              :: matrix(:,:)
    !! The output _exchange matrix_.
    INTEGER                                     :: nz,i,j,jj,jinf,jsup
    REAL(kind=mm_wp)                            :: zni,znip1,xf,xft,xcnt
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: puit
    REAL(kind=mm_wp)                            :: cpte,cpte2
    REAL(kind=mm_wp)                            :: zsurf
    INTEGER, PARAMETER                          :: ichx  = 1
    matrix = 0._mm_wp ; nz = SIZE(rad) ; ALLOCATE(puit(nz))
    zsurf = mm_zlev(nz+1)
    
    ! compute exchange matrix
    DO i=1,nz
      puit(i) = 0._mm_wp
      xcnt = 0._mm_wp
      ! computes drop move (i.e. its new positions)
      CALL getnzs(ichx,i,rad,rhog,dt,zni,znip1)

      ! Peculiar case : Ground level precipitation [znip1 < zsurf && (zni < zsurf || zni > zsurf)]
      ! - complete precipitation [ znip1 <= 0 && zni <= 0 ] :
      IF(zni <= zsurf .and. znip1 <= zsurf) THEN
        xft=0._mm_wp
        xf=1._mm_wp
        xcnt=xcnt+xf
        puit(i)=puit(i)+xf
      ENDIF
      ! - partial precipitation [ znip1 <= zsurf && zni > zsurf ] :
      IF (zni > zsurf .and. znip1 <= zsurf) THEN
        xft=(zni-zsurf)/(zni-znip1)
        xf=(1.-xft)
        xcnt=xcnt+xf
        puit(i)=puit(i)+xf
      ENDIF
      ! General case : no ground precipitation [ znip1 > zsurf && zni > zsurf ]
      IF (zni > zsurf .and. znip1 > zsurf) THEN
        xft = 1._mm_wp ! on a la totalite de la case
        xf  = 0._mm_wp
        xcnt=xcnt+xf
        puit(i)=puit(i)+xf
      ENDIF
      IF (zni < znip1) THEN
        WRITE(*,'("[EXCHANGES] WARNING, missing this case :",2(2X,ES10.3))') zni, znip1
      ENDIF
      
      ! Fix minimum level to the ground
      znip1 = MAX(znip1,zsurf)
      zni   = MAX(zni,zsurf)
      ! Locate new "drop" position in the verical grid
      jsup=nz+1
      jinf=nz+1
      DO j=1,nz
        IF (zni<=mm_zlev(j).and.zni>=mm_zlev(j+1)) jsup=j
        IF (znip1<=mm_zlev(j).and.znip1>=mm_zlev(j+1)) jinf=j
      ENDDO
      ! Volume is out of range: (all drops have touched the ground!)
      ! Note: cannot happen here, it has been treated previously :)
      IF (jsup>=nz+1.and.jinf==jsup) THEN
        WRITE(*,'(a)') "[EXCHANGE] speaking: The impossible happened !"
        call EXIT(666)
      ENDIF
      ! Volume inside a single level
      IF (jsup==jinf.and.jsup<=nz) THEN
        xf=1._mm_wp
        xcnt=xcnt+xft*xf
        matrix(jinf,i)=matrix(jinf,i)+xft*xf
      ENDIF

      ! Volume over 2 levels
      IF (jinf==jsup+1) THEN
        xf=(zni-mm_zlev(jinf))/(zni-znip1)
        xcnt=xcnt+xf*xft
        IF(jsup<=nz) THEN
          matrix(jsup,i)=matrix(jsup,i)+xft*xf
        ENDIF
        xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
        xcnt=xcnt+xf*xft
        IF (jinf<=nz) THEN
          matrix(jinf,i)=matrix(jinf,i)+xft*xf
        ENDIF
      ENDIF

      ! Volume over 3 or more levels
      IF (jinf > jsup+1) THEN
        ! first cell
        xf=(zni-mm_zlev(jsup+1))/(zni-znip1)
        xcnt=xcnt+xf*xft
        matrix(jsup,i)=matrix(jsup,i)+xft*xf
        ! last cell
        xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
        xcnt=xcnt+xf*xft
        matrix(jinf,i)=matrix(jinf,i)+xft*xf
        ! other :)
        DO jj=jsup+1,jinf-1
          xf=(mm_zlev(jj)-mm_zlev(jj+1))/(zni-znip1)
          xcnt=xcnt+xf*xft
          matrix(jj,i)=matrix(jj,i)+xft*xf
        ENDDO
      ENDIF
    ENDDO

    ! checking if everything is alright if debug enabled...
    IF (mm_debug) THEN
      cpte=0._mm_wp ; cpte2=0._mm_wp
      DO j=1,nz
        DO jj=1,nz
          cpte=cpte+matrix(jj,j)
        ENDDO
      ENDDO
      cpte2=cpte+sum(puit)
      IF (abs(cpte2-nz)>1.e-4_mm_wp) THEN
        WRITE(*,'("[EXCHANGE] speaking: tx expl (/nz):",2(2X,ES10.3))') cpte,cpte2
      ENDIF
    ENDIF

    RETURN
  END SUBROUTINE exchange

  SUBROUTINE getnzs(ichx,idx,rad,rho,dt,zni,zns)
    !! Compute displacement of a cell under sedimentation process.
    !!
    !! The method computes the new position of a _drop cell_ through sedimentation process as
    !! descibed in the following scheme:
    !!
    !! ![Cloud sedimentation scheme](|media|/cloud_sed_scheme.svg)
    !!
    !! New positions are returned in __zni__ and __zns__ ouptut arguments.
    !!
    !! @note
    !! The method uses directly [[mm_globals(module):mm_play(variable)]], [[mm_globals(module):mm_plev(variable)]],
    !! [[mm_globals(module):mm_temp(variable)]],[[mm_globals(module):mm_btemp(variable)]],
    !! [[mm_globals(module):mm_zlay(variable)]] and [[mm_globals(module):mm_zlev(variable)]] and uses __idx__ to
    !! get the relevant value to use on the vertical grid.
    INTEGER, INTENT(in)                        :: ichx
    !! Velocity extrapolation control flag (0 for linear, 1 for exponential -preferred -).
    INTEGER, INTENT(in)                        :: idx
    !! Initial position of the drop (subscript of vertical layers vectors).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
    !! Cloud drop radius over the atmospheric vertical structure (m).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rho
    !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)).
    REAL(kind=mm_wp), INTENT(in)               :: dt
    !! Timestep (s).
    REAL(kind=mm_wp), INTENT(out)              :: zni
    !! Final layer upper boundary position (m).
    REAL(kind=mm_wp), INTENT(out)              :: zns
    !! Final layer lower boundary position (m).
    REAL(kind=mm_wp)            :: ws,wi,w,zi,zs
    REAL(kind=mm_wp)            :: alpha,argexp,v0,arg1,arg2
    INTEGER                 :: i,nz
    REAL(kind=mm_wp), PARAMETER :: es = 30._mm_wp

    nz = SIZE(rad)
    ! Linear extrapolation of velocity
    IF (ichx==0) THEN
      ! velocity upper interface
      ws = wsettle(mm_plev(idx),mm_btemp(idx),mm_zlev(idx),rho(idx),rad(idx))
      IF (idx==nz) THEN
        ! veloctity center layer
        wi = wsettle(mm_play(idx),mm_temp(idx),mm_zlay(idx),rho(idx),rad(idx))
      ELSEIF(idx<nz) THEN
        ! velocity lower interface
        wi = wsettle(mm_plev(idx+1),mm_btemp(idx+1),mm_zlev(idx+1), &
          rho(idx+1),rad(idx+1))
      ELSE
        WRITE(*,'(a)') "[getnzs] speaking:"
        WRITE(*,'(a)') "This is the fatal error..."
        WRITE(*,'(a)') "index is higher than number of levels"
        call EXIT(111)
      ENDIF
      w   = (ws+wi)/2._mm_wp
      zni = mm_zlev(idx)-w*dt
      zns = mm_zlev(idx)-mm_dzlev(idx)-w*dt
      RETURN
      ! Exponential extrapolation of velocity

    ELSEIF(ichx==1) THEN
      zs = mm_zlev(idx)
      ws = wsettle(mm_plev(idx),mm_btemp(idx),zs,rho(idx),rad(idx))
      zi=mm_zlay(idx)
      wi = wsettle(mm_play(idx),mm_temp(idx),zi,rho(idx),rad(idx))
      ! ws & wi must be different !
      IF(dabs(wi-ws)/wi <= 1.e-3_mm_wp)  wi=ws/1.001_mm_wp
      IF (wi /= 0._mm_wp) alpha = dlog(ws/wi)/(zs-zi)  ! alpha < 0 if wi > ws
      !   -es < argexp < es
      argexp=MAX(MIN(alpha*zs,es),-es)
      v0 = ws/dexp(argexp)
      arg1=1._mm_wp+v0*alpha*dexp(argexp)*dt
      argexp=MAX(MIN(alpha*(mm_zlev(idx)-mm_dzlev(idx)),es),-es)
      arg2=1._mm_wp+v0*alpha*dexp(argexp)*dt
      IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
        ! correct velocity
        ! divides the velocity argument in arg1 and arg2 :
        ! argX=1+alpha*v0*exp(alpha*z)*dt <==> argX-1=alpha*v0*exp(alpha*z)*dt
        ! ==> argX' = 1 + (argX-1)/2  <==> argX' = (1+argX)/2.
        DO i=1,25
          IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
            IF (mm_debug) &
              WRITE(*,'((a),I2.2,(a))') "[getnzs] must adjust velocity (iter:",i,"/25)"
            arg1=(arg1+1._mm_wp)/2._mm_wp ; arg2=(arg2+1._mm_wp)/2._mm_wp
          ELSE
            EXIT
          ENDIF
        ENDDO
        ! sh** we have to stop
        IF (i>25) THEN
          WRITE(*,'(a)')"[getnzs] speaking:"
          WRITE(*,'(a)') "Cannot adjust velocities"
          call EXIT(111)
        ENDIF
      ENDIF

      zni = mm_zlev(idx)-dlog(arg1)/alpha
      zns = mm_zlev(idx)-mm_dzlev(idx)-dlog(arg2)/alpha

      RETURN
    ENDIF
  END SUBROUTINE getnzs

  ELEMENTAL FUNCTION wsettle(p,t,z,rho,rad) RESULT(w)
    !! Compute the settling velocity of a spherical particle.
    !!
    !! The method computes the effective settling velocity of spherical particle of
    !! radius __rad__. It accounts for the slip-flow transition (no approximation).
    REAL(kind=mm_wp), INTENT(in) :: p   !! The pressure level (Pa).
    REAL(kind=mm_wp), INTENT(in) :: t   !! The temperature (K).
    REAL(kind=mm_wp), INTENT(in) :: z   !! The altitude level (m).
    REAL(kind=mm_wp), INTENT(in) :: rho !! Density of the particle (\(kg.m^{-3}\)).
    REAL(kind=mm_wp), INTENT(in) :: rad !! Radius of the particle (m).
    REAL(kind=mm_wp) :: w               !! Settling velocity (\(m.s^{-1}\)).
    REAL(kind=mm_wp) :: Us, Fc, kn, wtmp, wmax
    REAL(kind=mm_wp), PARAMETER :: ra = 1.75e-10_mm_wp
    
    ! Knudsen number
    kn = (mm_kboltz * t) / (p * 4._mm_wp * sqrt(2._mm_wp) * mm_pi * ra**2) / rad

    ! Computes slip-flow correction : Fc = 1 + (mm_akn * mm_lambda_g(t,p) / rad)
    Fc = (1._mm_wp + 1.2517_mm_wp*kn + 0.4_mm_wp*kn*dexp(-1.1_mm_wp/kn))
    
    ! Computes Stokes settling velocity
    Us = (2._mm_wp * rad**2 * rho * mm_effg(z)) / (9._mm_wp * mm_eta_g(t))
    
    ! Computes settling velocity (correction factor : x3.0)
    w = Us * Fc * 3._mm_wp

    ! Imposes a velocity limit
    wmax = 20._mm_wp ! 20 m/s [Lorenz 1993]
    wtmp = (1._mm_wp / w) + (1._mm_wp / wmax)
    w = 1._mm_wp / wtmp
  END FUNCTION wsettle

  FUNCTION get_mass_flux(rho,m3) RESULT(flx)
    !> Get the mass flux of (clouds related) moment through sedimention.
    !!
    !! @warning
    !! The method is __only__ valid for cloud moments (i.e. ice or ccn). It calls
    !! [[mm_clouds(module):wsettle(function)]] that compute the _mean_ settling velocity of a
    !! cloud drop.
    !!
    !! @note
    !! The computed flux is always positive.
    REAL(kind=mm_wp), INTENT(in)               :: rho
    !! Tracer density (\(kg.m^{-3}\)).
    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3
    !! Vertical profile of the total volume of tracer (i.e. M3) from __TOP__ to __GROUND__ (\(m^{3}.m^{-3}\)).
    REAL(kind=mm_wp), DIMENSION(SIZE(m3))      :: flx
    !! Mass sedimentation fluxes at each layer from __TOP__ to __GROUND__ (\(kg.m^{-2}.s^{-1}\)).
    REAL(kind=mm_wp), SAVE :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp
    flx = rho * pifac * m3 * wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
    RETURN
  END FUNCTION get_mass_flux

END MODULE MM_CLOUDS