Index: trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90
===================================================================
--- trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90	(revision 2242)
+++ trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90	(revision 3083)
@@ -1,21 +1,21 @@
-! Copyright 2013-2015,2017 Université de Reims Champagne-Ardenne 
+! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
 ! Contributor: J. Burgalat (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, 
+! 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". 
-! 
+! "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. 
-! 
+! liability.
+!
 ! In this respect, the user's attention is drawn to the risks associated
 ! with loading,  using,  modifying and/or developing or reproducing the
@@ -25,8 +25,8 @@
 ! 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. 
-! 
+! 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.
@@ -35,14 +35,14 @@
 !! summary: Parameters and global variables module.
 !! author: J. Burgalat
-!! date: 2013-2015,2017
+!! date: 2013-2015,2017,2022
 
 MODULE MM_GLOBALS
   !! Parameters and global variables module.
-  !!  
+  !!
   !! # Module overview
   !!
   !! The module defines all the parameters and global variables that are common
   !! to all other modules of the library.
-  !! 
+  !!
   !! It is separated in two parts :
   !!
@@ -52,27 +52,27 @@
   !!   method.
   !! - The second part defines a set of vectors that defines the vertical structure of the atmosphere.
-  !!   Each time a new atmospheric column has to be computed (either on a new timestep or on a new couple 
-  !!   of longitude/latitude), these vectors should be intialized with new values by calling 
-  !!   [[mm_globals(module):mm_column_init(function)]] method. 
-  !!   This part is separated in two sets : 
-  !!
-  !!   - The atmospheric structure with temperature, pressure levels and altitude definitions. 
-  !!   - The vertical profiles of tracers with the moments of the two aerosols modes (both \(M_{0}\) 
-  !!     and \(M_{3}\) for a total of 4 vectors), the _clouds_ microphysics moments tracers (i.e. 
+  !!   Each time a new atmospheric column has to be computed (either on a new timestep or on a new couple
+  !!   of longitude/latitude), these vectors should be intialized with new values by calling
+  !!   [[mm_globals(module):mm_column_init(function)]] method.
+  !!   This part is separated in two sets :
+  !!
+  !!   - The atmospheric structure with temperature, pressure levels and altitude definitions.
+  !!   - The vertical profiles of tracers with the moments of the two aerosols modes (both \(M_{0}\)
+  !!     and \(M_{3}\) for a total of 4 vectors), the _clouds_ microphysics moments tracers (i.e.
   !!     \(M_{0}\) and \(M_{3}\) for the ccn and \(M_{3}\) for the ice components).
-  !!     Additionally, the module also stores intermediates variables of interest such as the 
-  !!     characteristic radii of the aerosols modes, the mean drop radius and the drop density, 
+  !!     Additionally, the module also stores intermediates variables of interest such as the
+  !!     characteristic radii of the aerosols modes, the mean drop radius and the drop density,
   !!     the molar fraction of each condensible species (related to ice components) and some
   !!     scalar variables that holds arrays sizes.
   !!
   !! @note
-  !! All the vectors that represent the vertical structure of the atmosphere (altitude, pressure and 
+  !! All the vectors that represent the vertical structure of the atmosphere (altitude, pressure and
   !! temperature...) are oriented from the __TOP__ of the atmosphere to the __GROUND__.
   !!
-  !! @note 
-  !! The module also imports errors module from __FCCP__ library to get definitions of the error object 
+  !! @note
+  !! The module also imports errors module from __SWIFT__ library to get definitions of the error object
   !! everywhere in the library ([[mm_globals(module)]] is always imported, except in [[mm_mprec(module)]]).
   !!
-  !! # Global variables 
+  !! # Global variables
   !!
   !! [[mm_globals(module)]] module contains the declaration of all global/common variable that are shared
@@ -82,6 +82,6 @@
   !! the following sections list all the global variables by category.
   !!
-  !! ## Control flags 
-  !! 
+  !! ## Control flags
+  !!
   !! | Name               | Description
   !! | :----------------- | :-----------------
@@ -93,5 +93,5 @@
   !! | mm_w_clouds_sed    | Enable/Disable clouds microphysics sedimentation
   !! | mm_w_clouds_nucond | Enable/Disable clouds microphysics nucleation/condensation
-  !! | mm_wsed_m0         | Force all aerosols moments to fall at M0 settling velocity 
+  !! | mm_wsed_m0         | Force all aerosols moments to fall at M0 settling velocity
   !! | mm_wsed_m3         | Force all aerosols moments to fall at M3 settling velocity
   !! | mm_no_fiadero_w    | Enable/Disable __Fiadero__ correction
@@ -101,9 +101,9 @@
   !! | Name            | Description
   !! | :-------------- | :-----------------
-  !! | mm_fiadero_min  | Minimum ratio for __Fiadero__'s correction 
+  !! | mm_fiadero_min  | Minimum ratio for __Fiadero__'s correction
   !! | mm_fiadero_max  | Maximum ratio for __Fiadero__'s correction
-  !! | mm_coag_choice  | Coagulation interaction activation flag. It should be a combination of [[mm_globals(module):mm_coag_no(variable)]], [[mm_globals(module):mm_coag_ss(variable)]], [[mm_globals(module):mm_coag_sf(variable)]] and [[mm_globals(module):mm_coag_ff(variable)]]. 
-  !!
-  !! ## Physical constants 
+  !! | mm_coag_choice  | Coagulation interaction activation flag. It should be a combination of [[mm_globals(module):mm_coag_no(variable)]], [[mm_globals(module):mm_coag_ss(variable)]], [[mm_globals(module):mm_coag_sf(variable)]] and [[mm_globals(module):mm_coag_ff(variable)]].
+  !!
+  !! ## Physical constants
   !!
   !! | Name      | Description
@@ -131,5 +131,5 @@
   !! | mm_w_prod   | Angular frequency of the time-dependent production rate.
   !! | mm_ne       | Electric charging of aerosols (\(e^{-}.m^{-1}\)) (unused)
-  !! | mm_rb2ra    | Bulk to apparent radius conversion pre-factor (\(m^X\)) 
+  !! | mm_rb2ra    | Bulk to apparent radius conversion pre-factor (\(m^X\))
   !! | mm_rpla     | Planet radius (m)
   !! | mm_g0       | Planet acceleration due to gravity constant (ground) (\(m.s^{-2}\))
@@ -152,5 +152,5 @@
   ! the following variables are read-only outside this module.
   ! One must call the afferent subroutine to update them.
-    
+
   ! initialization control flags (cannot be updated)
   PROTECTED :: mm_ini,mm_ini_col,mm_ini_aer,mm_ini_cld
@@ -165,7 +165,9 @@
   ! Moments parameters (derived, are updated with moments parameters)
   PROTECTED :: mm_rcs, mm_rcf, mm_drad, mm_drho
-
-  LOGICAL, SAVE :: mm_debug = .true.  !! Enable QnD debug mode (can be used for devel).
-  LOGICAL, SAVE :: mm_log = .false.   !! Enable log mode (for configuration only).
+  ! Thresholds parameters
+  PROTECTED :: mm_m0as_min, mm_m3as_min, mm_rcs_min, mm_m0af_min, mm_m3af_min, mm_rcf_min, mm_m0n_min, mm_m3cld_min
+
+  LOGICAL, SAVE :: mm_debug = .false.  !! Enable QnD debug mode (can be used for devel).
+  LOGICAL, SAVE :: mm_log   = .false.   !! Enable log mode (for configuration only).
 
   LOGICAL, SAVE :: mm_w_haze_prod = .true. !! Enable/Disable haze production.
@@ -182,17 +184,17 @@
   !> Enable/Disable __Fiadero__'s correction.
   !!
-  !! This flag enables/disables the __Fiadero__ correction alogrithm for fractal mode settling velocity 
-  !! computation. 
+  !! This flag enables/disables the __Fiadero__ correction alogrithm for fractal mode settling velocity
+  !! computation.
   !!
   !! @bug
-  !! Currently, the Fiadero correction creates instatibilities on the vertical structure. It seems to be 
+  !! Currently, the Fiadero correction creates instatibilities on the vertical structure. It seems to be
   !! related to the coupling between the two moments. In order to reduce the instabilities, settling
   !! velocity of moments are forced to be the same, see [[mm_globals(module):mm_wsed_m0(variable)]] and
   !! [[mm_globals(module):mm_wsed_m3(variable)]]).
-  LOGICAL, SAVE          :: mm_no_fiadero_w = .false. 
+  LOGICAL, SAVE          :: mm_no_fiadero_w = .false.
 
   !> Minimum ratio for __Fiadero__ correction.
   !!
-  !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the minimum 
+  !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the minimum
   !! value of the moment's ratio between two adjacents vertical cells to be used within the correction.
   REAL(kind=mm_wp), SAVE :: mm_fiadero_min  = 0.1_mm_wp
@@ -200,5 +202,5 @@
   !> Maximum ratio for __Fiadero__ correction.
   !!
-  !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the maximum 
+  !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the maximum
   !! value of the moment's ratio between two adjacents vertical cells to be used within the correction.
   REAL(kind=mm_wp), SAVE :: mm_fiadero_max  = 10._mm_wp
@@ -213,10 +215,10 @@
   INTEGER, PARAMETER :: mm_coag_ff = 4 !! FF mode interaction for coagulation.
   !> Default interactions to activate (all by default).
-  INTEGER, SAVE      :: mm_coag_choice = mm_coag_ss+mm_coag_sf+mm_coag_ff 
+  INTEGER, SAVE      :: mm_coag_choice = mm_coag_ss+mm_coag_sf+mm_coag_ff
 
   !> Pi number.
   REAL(kind=mm_wp), PARAMETER :: mm_pi = 4._mm_wp*atan(1._mm_wp)
   !> Avogadro number.
-  REAL(kind=mm_wp), PARAMETER :: mm_navo = 6.0221367e23_mm_wp 
+  REAL(kind=mm_wp), PARAMETER :: mm_navo = 6.0221367e23_mm_wp
   !> Boltzmann constant (\(J.K^{-1}\)).
   REAL(kind=mm_wp), PARAMETER :: mm_kboltz = 1.3806488e-23_mm_wp
@@ -261,16 +263,48 @@
 
   !> Bulk to apparent radius conversion pre-factor (\(m^{X}\)).
-  !! 
-  !! It is initialized using [[mm_globals(module):mm_rm(variable)]] in 
+  !!
+  !! It is initialized using [[mm_globals(module):mm_rm(variable)]] in
   !! [[mm_globals(module):mm_global_init(interface)]] from the following equation:
   !!
   !! $$ r_{a} = r_{b}^{3/D_{f}}\times r_{m}^{\frac{D_{f}-3}{D_{f}}} $$
   !!
-  !! Where \(r_{a}\) is the apparent radius, \(r_{b}\) the bulk radius and 
+  !! Where \(r_{a}\) is the apparent radius, \(r_{b}\) the bulk radius and
   !! \(rb2ra = r_{m}^{\frac{D_{f}-3}{D_{f}}}\) is the returned pre-factor
-  REAL(kind=mm_wp), SAVE :: mm_rb2ra = 1._mm_wp 
+  REAL(kind=mm_wp), SAVE :: mm_rb2ra = 1._mm_wp
+
+  ! Thresholds !
+
+  !> (min.) Total number of aerosols minimum threshold for the spherical mode.
+  REAL(kind=mm_wp), SAVE :: mm_m0as_min = 1.e-10_mm_wp
+
+  !> (min.) Total volume of aerosols minimum threshold for the spherical mode.
+  REAL(kind=mm_wp), SAVE :: mm_m3as_min = 1.e-40_mm_wp
+
+  !> Characteristic radius minimum threshold for the spherical mode.
+  REAL(kind=mm_wp), SAVE :: mm_rcs_min = 1.e-9_mm_wp
+
+  !> (min.) Total number of aerosols minimum threshold for the fractal mode.
+  REAL(kind=mm_wp), SAVE :: mm_m0af_min = 1.e-10_mm_wp
+
+  !> (min.) Total volume of aerosols minimum threshold for the fractal mode.
+  REAL(kind=mm_wp), SAVE :: mm_m3af_min = 1.e-40_mm_wp
+
+  !> Characteristic radius minimum threshold for the fractal mode.
+  REAL(kind=mm_wp), SAVE :: mm_rcf_min = 1.e-9_mm_wp
+
+  !> (min.) Total number of cloud drop minimum threshold.
+  REAL(kind=mm_wp), SAVE :: mm_m0n_min = 1.e-10_mm_wp
+
+  !> (min.) Total volume of cloud drop minimum threshold.
+  REAL(kind=mm_wp), SAVE :: mm_m3cld_min = 1.e-40_mm_wp
 
   !> Characteristic radius threshold.
   REAL(kind=mm_wp), SAVE :: mm_rc_min = 1.e-200_mm_wp
+
+  !> Minimum cloud drop radius
+  REAL(kind=mm_wp), SAVE :: mm_drad_min = 1.e-9_mm_wp
+
+  !> Maximum cloud drop radius
+  REAL(kind=mm_wp), SAVE :: mm_drad_max = 1.e-3_mm_wp
 
   !> Name of condensible species.
@@ -339,11 +373,11 @@
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_temp
   !>  Air density vertical profile (\(kg.m^{-3}\)).
-  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rhoair 
+  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rhoair
   !> Temperature vertical profil at interfaces (K).
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_btemp
 
   !> Atmospheric levels thickness (m).
-  !! 
-  !! Atmospheric thickness between two adjacent levels (\(m\)) from the 
+  !!
+  !! Atmospheric thickness between two adjacent levels (\(m\)) from the
   !! __TOP__ to the __GROUND__.
   !! @note __mm_dzlev__ is defined on the total number of layers and actually
@@ -352,11 +386,11 @@
 
   !> Atmospheric layers "thickness" (m).
-  !! 
+  !!
   !! Atmospheric thickness between the center of two adjacent layers (\(m\))
   !! from the __TOP__ to the __GROUND__.
-  !! @note 
-  !! __mm_dzlay__ is defined on the total number of layers. The last 
+  !! @note
+  !! __mm_dzlay__ is defined on the total number of layers. The last
   !! value of __mm_dzlay__ is set to twice the altitude of the ground layer.
-  !! @note This value corresponds to the thickness between the center of the 
+  !! @note This value corresponds to the thickness between the center of the
   !! __GROUND__ layer and below the surface. It is arbitrary and not used.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlay
@@ -377,7 +411,7 @@
   !> Ice components 3rd order moments (\(m^{3}.m^{-3}\)).
   !!
-  !! It is a 2D array with the vertical layers in first dimension, and the number of ice 
+  !! It is a 2D array with the vertical layers in first dimension, and the number of ice
   !! components in the second.
-  !! @note 
+  !! @note
   !! Both [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]]
   !! share the same indexing (related to species order).
@@ -387,6 +421,6 @@
   !!
   !! It is a 2D array with the vertical layers in first dimension, and
-  !! the number of condensible species in the second. 
-  !! @note 
+  !! the number of condensible species in the second.
+  !! @note
   !! Both [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]]
   !! share the same indexing (related to species order).
@@ -410,8 +444,8 @@
   !> Spherical mode \(M_{0}\) settling velocity (\(m.s^{-1}\)).
   !!
-  !! It is a vector with the vertical layers that contains the settling velocity for 
+  !! It is a vector with the vertical layers that contains the settling velocity for
   !! the \(0^{th}\) order moment of the spherical mode.
   !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
-  !! @note 
+  !! @note
   !! This variable is always negative.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0as_vsed
@@ -419,8 +453,8 @@
   !> Spherical mode \(M_{3}\) settling velocity (\(m.s^{-1}\)).
   !!
-  !! It is a vector with the vertical layers that contains the settling velocity for the 
+  !! It is a vector with the vertical layers that contains the settling velocity for the
   !! \(3^{rd}\) order moment of the spherical mode.
   !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
-  !! @note 
+  !! @note
   !! This variable is always negative.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3as_vsed
@@ -428,8 +462,8 @@
   !> Fractal mode \(M_{0}\) settling velocity (\(m.s^{-1}\)).
   !!
-  !! It is a vector with the vertical layers that contains the settling velocity for the 
+  !! It is a vector with the vertical layers that contains the settling velocity for the
   !! \(0^{th}\) order moment of the fractal mode.
   !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
-  !! @note 
+  !! @note
   !! This variable is always negative.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0af_vsed
@@ -437,8 +471,8 @@
   !> Fractal mode \(M_{3}\) settling velocity (\(m.s^{-1}\)).
   !!
-  !! It is a vector with the vertical layers that contains the settling velocity for the 
+  !! It is a vector with the vertical layers that contains the settling velocity for the
   !! \(3^{rd}\) order moment of the fractal mode.
   !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
-  !! @note 
+  !! @note
   !! This variable is always negative.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3af_vsed
@@ -448,5 +482,5 @@
   !! It is a vector with the vertical layers that contains the mass fluxes for spherical aerosols.
   !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
-  !! @note 
+  !! @note
   !! This variable is always negative.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_s_flux
@@ -456,5 +490,5 @@
   !! It is a vector with the vertical layers that contains the mass fluxes for fractal aerosols
   !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
-  !! @note 
+  !! @note
   !! This variable is always negative.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_f_flux
@@ -464,16 +498,25 @@
   REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp
 
-  !> CCN mass fluxes (\(kg.m^{-2}.s^{-1}\)).
-  !!
-  !! It is a vector with the vertical layers that contains the 
-  !! mass fluxes for CCN. 
+  !> CCN settling velocity (\(m.s^{-1}\)).
+  !!
+  !! It is a vector with the vertical layers that contains the
+  !! settling velocity for CCN (and ices).
   !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
   !! @note
-  !! This variable is always negative.
+  !! This variable is always positive.
+  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_w
+  
+  !> CCN mass fluxes (\(kg.m^{-2}.s^{-1}\)).
+  !!
+  !! It is a vector with the vertical layers that contains the
+  !! mass fluxes for CCN.
+  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
+  !! @note
+  !! This variable is always positive.
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_flux
 
   !> Ice components precipitations (m).
   !!
-  !! It is a vector of [[mm_globals(module):mm_nesp(variable)]] values which share the same indexing 
+  !! It is a vector of [[mm_globals(module):mm_nesp(variable)]] values which share the same indexing
   !! than [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]].
   !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
@@ -484,5 +527,5 @@
   !> Ice components sedimentation fluxes (\(kg.m^{-2}.s-1\)).
   !!
-  !! It is a 2D-array with the vertical layers in first dimension and the number of ice components 
+  !! It is a 2D-array with the vertical layers in first dimension and the number of ice components
   !! in the second. It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
   !! @note
@@ -492,5 +535,5 @@
   !> Condensible species saturation ratio (--).
   !!
-  !! It is a 2D-array with the vertical layers in first dimension and the number of condensible 
+  !! It is a 2D-array with the vertical layers in first dimension and the number of condensible
   !! species in the second.
   !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
@@ -516,5 +559,5 @@
   INTERFACE mm_cloud_properties
     MODULE PROCEDURE cldprop_sc,cldprop_ve
-  END INTERFACE
+  END INTERFACE mm_cloud_properties
 
   !> Interface to global initialization.
@@ -522,14 +565,13 @@
   !! The method performs the global initialization of the model.
   !! @warning
-  !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it 
+  !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it
   !! initializes global variable that are not thread private.
   !!
-  !! '''
-  !! !$OMP SINGLE
-  !! call mm_global_init(...)
-  !! !$OMP END SINGLE
+  !!   !$OMP SINGLE
+  !!   call mm_global_init(...)
+  !!   !$OMP END SINGLE
   INTERFACE mm_global_init
     MODULE PROCEDURE mm_global_init_0,mm_global_init_1
-  END INTERFACE
+  END INTERFACE mm_global_init
 
   !> Check an option from the configuration system.
@@ -538,11 +580,11 @@
   !! set a default value if the option is not found. This is an overloaded method
   !! that can take in input either a floating point, integer, logical or string
-  !! option value. 
+  !! option value.
   INTERFACE mm_check_opt
     MODULE PROCEDURE check_r1,check_i1,check_l1,check_s1
-  END INTERFACE
+  END INTERFACE mm_check_opt
 
   ! --- OPENMP ---------------
-  ! All variable related to column computations should be private to each thread
+  ! All variables related to column computations should be private to each thread
   !
   !$OMP THREADPRIVATE(mm_ini_col,mm_ini_aer,mm_ini_cld)
@@ -551,7 +593,7 @@
   !$OMP THREADPRIVATE(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_gazs)
   !$OMP THREADPRIVATE(mm_rcs,mm_rcf,mm_drad,mm_drho)
-  !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_flux,mm_ice_prec,mm_ice_fluxes,mm_gazs_sat)
+  !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_w,mm_ccn_flux,mm_ice_prec,mm_ice_fluxes,mm_gazs_sat)
   !$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed)
-
+  !$OMP THREADPRIVATE(mm_m0as_min,mm_m3as_min,mm_rcs_min,mm_m0af_min,mm_m3af_min,mm_rcf_min,mm_m0n_min,mm_m3cld_min)
   !$OMP THREADPRIVATE(mm_nla,mm_nle)
 
@@ -559,88 +601,97 @@
 
 
-  CONTAINS 
+CONTAINS
 
   FUNCTION mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
-                            air_rad,air_mmol,coag_interactions,clouds,spcfile,  &
-                            w_haze_prod,w_haze_sed,w_haze_coag,w_cloud_nucond,  &
-                            w_cloud_sed,force_wsed_to_m0,force_wsed_to_m3,      &
-                            no_fiadero,fiadero_min,fiadero_max) RESULT(err)
+    air_rad,air_mmol,coag_interactions,clouds,spcfile,  &
+    w_haze_prod,w_haze_sed,w_haze_coag,w_cloud_nucond,  &
+    w_cloud_sed,force_wsed_to_m0,force_wsed_to_m3,      &
+    no_fiadero,fiadero_min,fiadero_max,                 &
+    m0as_min,rcs_min,m0af_min,rcf_min,m0n_min,debug) RESULT(err)
     !! Initialize global parameters of the model.
-    !! 
+    !!
     !! The function initializes all the global parameters of the model from direct input.
-    !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their 
-    !! default values are suitable for production runs.  
+    !! Boolean, Fiadero and thresholds parameters are optional as they are rather testing parameters.
+    !! Their default values are suitable for production runs.
     !! @note
     !! If the method fails to initialize parameters (i.e. returned error is not 0). Then the model
     !! should probably be aborted as the global variables of the model will not be correctly setup.
     !! @warning
-    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it 
-    !! initializes global variable that are not thread private.
-    !!
-    !! '''
-    !! !$OMP SINGLE
-    !! call mm_global_init_0(...)
-    !! !$OMP END SINGLE
+    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it
+    !! initializes (only) global variables that are not thread private.
+    !!
+    !!   !$OMP SINGLE
+    !!   call mm_global_init_0(...)
+    !!   !$OMP END SINGLE
     REAL(kind=mm_wp), INTENT(in)           :: dt
-      !! Microphysics timestep in seconds.
+    !! Microphysics timestep in seconds.
     REAL(kind=mm_wp), INTENT(in)           :: df
-      !! Fractal dimension of fractal aerosol.
+    !! Fractal dimension of fractal aerosol.
     REAL(kind=mm_wp), INTENT(in)           :: rm
-      !! Monomer radius in meter.
+    !! Monomer radius in meter.
     REAL(kind=mm_wp), INTENT(in)           :: rho_aer
-      !! Aerosol density in \(kg.m^{-3}\).
+    !! Aerosol density in \(kg.m^{-3}\).
     REAL(kind=mm_wp), INTENT(in)           :: p_prod
-      !!  Aerosol production pressure level in Pa.
+    !!  Aerosol production pressure level in Pa.
     REAL(kind=mm_wp), INTENT(in)           :: tx_prod
-      !! Spherical aerosol mode production rate in \(kg.m^{-2}.s^{-1}\).
+    !! Spherical aerosol mode production rate in \(kg.m^{-2}.s^{-1}\).
     REAL(kind=mm_wp), INTENT(in)           :: rc_prod
-      !! Spherical mode characteristic radius for production in meter.
+    !! Spherical mode characteristic radius for production in meter.
     REAL(kind=mm_wp), INTENT(in)           :: rplanet
-      !! Planet radius in meter
+    !! Planet radius in meter
     REAL(kind=mm_wp), INTENT(in)           :: g0
-      !! Planet gravity acceleration at ground level in \(m.s^{-2}\).
+    !! Planet gravity acceleration at ground level in \(m.s^{-2}\).
     REAL(kind=mm_wp), INTENT(in)           :: air_rad
-      !! Air molecules mean radius in meter.
+    !! Air molecules mean radius in meter.
     REAL(kind=mm_wp), INTENT(in)           :: air_mmol
-      !! Air molecules mean molar mass in \(kg.mol^{-1}\).
+    !! Air molecules mean molar mass in \(kg.mol^{-1}\).
     INTEGER, INTENT(in)                    :: coag_interactions
-      !! Coagulation interactions process control flag.
+    !! Coagulation interactions process control flag.
     LOGICAL, INTENT(in)                    :: clouds
-      !! Clouds microphysics control flag.
+    !! Clouds microphysics control flag.
     CHARACTER(len=*), INTENT(in)           :: spcfile
-      !! Clouds microphysics condensible species properties file.
+    !! Clouds microphysics condensible species properties file.
     REAL(kind=mm_wp), INTENT(in), OPTIONAL :: fiadero_max
-      !! Maximum moment ratio threshold for Fiadero correction (default: 10.) .
+    !! Maximum moment ratio threshold for Fiadero correction (default: 10.) .
     REAL(kind=mm_wp), INTENT(in), OPTIONAL :: fiadero_min
-      !! Minimum moment ratio threshold for Fiadero correction (default: 0.1).
+    !! Minimum moment ratio threshold for Fiadero correction (default: 0.1).
     LOGICAL, INTENT(in), OPTIONAL          :: w_haze_prod
-      !! Haze microphysics production process control flag (default: T).
+    !! Haze microphysics production process control flag (default: T).
     LOGICAL, INTENT(in), OPTIONAL          :: w_haze_sed
-      !! Haze microphysics sedimentation process control flag (default: T).
+    !! Haze microphysics sedimentation process control flag (default: T).
     LOGICAL, INTENT(in), OPTIONAL          :: w_haze_coag
-      !! Haze microphysics coagulation process control flag (default: T).
+    !! Haze microphysics coagulation process control flag (default: T).
     LOGICAL, INTENT(in), OPTIONAL          :: w_cloud_sed
-      !! Cloud microphysics nucleation/conensation process control flag (default: __clouds__ value).
+    !! Cloud microphysics nucleation/conensation process control flag (default: __clouds__ value).
     LOGICAL, INTENT(in), OPTIONAL          :: w_cloud_nucond
-      !! Cloud microphysics production process control flag (default: __clouds__ value).
+    !! Cloud microphysics production process control flag (default: __clouds__ value).
     LOGICAL, INTENT(in), OPTIONAL          :: no_fiadero
-      !! Disable Fiadero correction for haze sedimentation process (default: F).
+    !! Disable Fiadero correction for haze sedimentation process (default: F).
     LOGICAL, INTENT(in), OPTIONAL          :: force_wsed_to_m0
-      !! force __all__ aerosols moments to fall at M0 settling velocity (default: T).
-    LOGICAL, INTENT(in), OPTIONAL          :: force_wsed_to_m3 
-      !! Force __all__ aerosols moments to fall at M3 settling velocity (default: F).
+    !! force __all__ aerosols moments to fall at M0 settling velocity (default: T).
+    LOGICAL, INTENT(in), OPTIONAL          :: force_wsed_to_m3
+    !! Force __all__ aerosols moments to fall at M3 settling velocity (default: F).
+    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: m0as_min
+    !! Minimum threshold for M0 of the spherical mode (default: 1e-10).
+    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: rcs_min
+    !! Minimum threshold for the characteristic radius of the spherical mode in meter (default: 1e-9).
+    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: m0af_min
+    !! Minimum threshold for M0 of the factal mode (default: 1e-10).
+    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: rcf_min
+    !! Minimum threshold for the characteristic radius of the fractal mode in meter (default: _monomer radius_).
+    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: m0n_min
+    !! Minimum threshold for M0 of cloud drop (default: 1e-10).
+    LOGICAL, INTENT(in), OPTIONAL          :: debug
+    !! Debug mode control flag (may print lot of stuff if enabled)
     TYPE(error) :: err
-      !! Error status of the function.
+    !! Error status of the function.
     INTEGER                                           :: i
     TYPE(cfgparser)                                   :: cp
-    CHARACTER(len=st_slen)                            :: spcpath
-    CHARACTER(len=:), ALLOCATABLE                     :: defmsg
     CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
     REAL(kind=mm_wp)                                  :: zfiamin,zfiamax
     LOGICAL                                           :: zwhp,zwhs,zwhc,zwcs,zwcn,znofia, &
                                                          zwstom0,zwstom3
-
     zwhp = .true. ; zwhs = .true. ; zwhc = .true.
-    zwcs = clouds ; zwcn = clouds 
+    zwcs = clouds ; zwcn = clouds
     znofia = .false. ; zfiamin = 0.1_mm_wp ; zfiamax = 10._mm_wp
     zwstom0 = .true. ; zwstom3 = .false.
@@ -652,8 +703,8 @@
 
     ! Store options values in global variables...
-    mm_df          = df 
-    mm_rm          = rm 
+    mm_df          = df
+    mm_rm          = rm
     mm_rb2ra       = mm_rm**((mm_df-3._mm_wp)/mm_df) ! conversion factor for bulk -> fractal radius
-    mm_rhoaer      = rho_aer 
+    mm_rhoaer      = rho_aer
     mm_p_prod      = p_prod
     mm_tx_prod     = tx_prod
@@ -662,5 +713,5 @@
     mm_g0          = g0
     mm_dt          = dt
-    mm_air_rad     = mm_air_rad
+    mm_air_rad     = air_rad
     mm_air_mmol    = air_mmol
     mm_coag_choice = coag_interactions
@@ -670,4 +721,7 @@
       RETURN
     ENDIF
+
+    ! force fractal radius minimum threshold to monomer radius ^^
+    mm_rcf_min = mm_rm
 
     mm_w_clouds = clouds
@@ -680,7 +734,7 @@
         RETURN
       ENDIF
-      ! Reads species properties configuration file  
+      ! Reads species properties configuration file
       err = cfg_read_config(cp,TRIM(spcfile)) ; IF (err /= 0) RETURN
-      err = cfg_get_value(cp,"used_species",species) 
+      err = cfg_get_value(cp,"used_species",species)
       IF (err /= 0) THEN
         err = error("mm_global_init: cannot retrieve 'used_species' values",-1)
@@ -691,5 +745,5 @@
       ALLOCATE(mm_spcname(mm_nesp),mm_xESPS(mm_nesp))
       DO i=1,mm_nesp
-        mm_spcname(i) = to_lower(species(i))
+        mm_spcname(i) = TRIM(species(i))
         IF(.NOT.cfg_has_section(cp,TRIM(mm_spcname(i)))) THEN
           err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1)
@@ -708,66 +762,107 @@
 
     ! optional flags
+    ! debug mode
+    IF (PRESENT(debug)) THEN
+      mm_debug = debug
+    ELSE
+      mm_debug = .false.
+      call printw("mm_debug",to_string(mm_debug))
+    ENDIF
     ! haze control flags
-    IF (PRESENT(w_haze_prod)) THEN 
+    IF (PRESENT(w_haze_prod)) THEN
       mm_w_haze_prod = w_haze_prod
-    ELSE 
-      mm_w_haze_prod = zwhp 
+    ELSE
+      mm_w_haze_prod = zwhp
       call printw("mm_haze_production",to_string(mm_w_haze_prod))
     ENDIF
-    IF (PRESENT(w_haze_sed)) THEN 
+    IF (PRESENT(w_haze_sed)) THEN
       mm_w_haze_sed = w_haze_sed
-    ELSE 
-      mm_w_haze_sed = zwhs 
+    ELSE
+      mm_w_haze_sed = zwhs
       call printw("mm_haze_sedimentation",to_string(mm_w_haze_sed))
     ENDIF
-    IF (PRESENT(w_haze_coag)) THEN 
+    IF (PRESENT(w_haze_coag)) THEN
       mm_w_haze_coag = w_haze_coag
-    ELSE 
+    ELSE
       mm_w_haze_coag = zwhc
       call printw("mm_haze_coagulation",to_string(mm_w_haze_coag))
     ENDIF
-    IF (PRESENT(force_wsed_to_m0)) THEN 
+    IF (PRESENT(force_wsed_to_m0)) THEN
       mm_wsed_m0 = force_wsed_to_m0
-    ELSE 
+    ELSE
       mm_wsed_m0 = zwstom0
       call printw("mm_wsed_m0",to_string(mm_wsed_m0))
     ENDIF
-    IF (PRESENT(force_wsed_to_m3)) THEN 
+    IF (PRESENT(force_wsed_to_m3)) THEN
       mm_wsed_m3 = force_wsed_to_m3
-    ELSE 
+    ELSE
       mm_wsed_m3 = zwstom3
       call printw("mm_wsed_m3",to_string(mm_wsed_m3))
     ENDIF
-    IF (PRESENT(no_fiadero)) THEN 
+    IF (PRESENT(no_fiadero)) THEN
       mm_no_fiadero_w = no_fiadero
-    ELSE 
-      mm_no_fiadero_w = znofia 
+    ELSE
+      mm_no_fiadero_w = znofia
       call printw("mm_no_fiadero",to_string(mm_no_fiadero_w))
     ENDIF
-    IF (PRESENT(fiadero_min)) THEN 
+    IF (PRESENT(fiadero_min)) THEN
       mm_fiadero_min = fiadero_min
-    ELSE 
+    ELSE
       mm_fiadero_min = zfiamin
       call printw("mm_fiadero_min",to_string(mm_fiadero_min))
     ENDIF
-    IF (PRESENT(fiadero_max)) THEN 
+    IF (PRESENT(fiadero_max)) THEN
       mm_fiadero_max = fiadero_max
-    ELSE 
+    ELSE
       mm_fiadero_max = zfiamax
       call printw("mm_fiadero_max",to_string(mm_fiadero_max))
     ENDIF
+
+    ! moments threshold flags
+    IF (PRESENT(m0as_min)) THEN
+      mm_m0as_min = MAX(0._mm_wp,m0as_min)
+    ELSE
+      call printw("mm_m0as_min",to_string(mm_m0as_min))
+    ENDIF
+    IF (PRESENT(rcs_min)) THEN
+      mm_rcs_min = MAX(1.e-9_mm_wp,rcs_min)
+    ELSE
+      call printw("mm_rcs_min",to_string(mm_rcs_min))
+    ENDIF
+    IF (PRESENT(m0af_min)) THEN
+      mm_m0af_min = MAX(0._mm_wp,m0af_min)
+    ELSE
+      call printw("mm_m0af_min",to_string(mm_m0af_min))
+    ENDIF
+    IF (PRESENT(rcf_min)) THEN
+      mm_rcf_min = MAX(rcf_min,mm_rm)
+    ELSE
+      mm_rcf_min = mm_rm
+      call printw("mm_rcf_min",to_string(mm_rcf_min))
+    ENDIF
+    IF (PRESENT(m0n_min)) THEN
+      mm_m0n_min = MAX(0._mm_wp,m0n_min)
+    ELSE
+      call printw("mm_m0n_min",to_string(mm_m0n_min))
+    ENDIF
+
+    ! compute m3 thresholds from user-defined thresholds.
+    mm_m3as_min  =  mm_m0as_min*mm_alpha_s(3._mm_wp) * mm_rcs_min**3._mm_wp
+    mm_m3af_min  =  mm_m0af_min*mm_alpha_f(3._mm_wp) * mm_rcf_min**3._mm_wp
+    mm_m3cld_min =  mm_m0n_min * (4._mm_wp * mm_pi / 3._mm_wp) * mm_drad_min**3._mm_wp
+
     ! clouds control flags
     IF (mm_w_clouds) THEN
-      IF (PRESENT(w_cloud_sed)) THEN 
+      IF (PRESENT(w_cloud_sed)) THEN
         mm_w_cloud_sed = w_cloud_sed
-      ELSE 
-        mm_w_cloud_sed = zwcs 
-        call printw("mm_cloud_sed",to_string(mm_w_cloud_sed)) 
+      ELSE
+        mm_w_cloud_sed = zwcs
+        call printw("mm_cloud_sed",to_string(mm_w_cloud_sed))
       ENDIF
-      IF (PRESENT(w_cloud_nucond)) THEN 
+      IF (PRESENT(w_cloud_nucond)) THEN
         mm_w_cloud_nucond = w_cloud_nucond
-      ELSE 
+      ELSE
         mm_w_cloud_nucond = zwcs
-        call printw("mm_cloud_nucond",to_string(mm_w_cloud_nucond)) 
+        call printw("mm_cloud_nucond",to_string(mm_w_cloud_nucond))
       ENDIF
     ENDIF
@@ -781,5 +876,5 @@
     mm_ini = err == noerror
 
-    CONTAINS
+  CONTAINS
 
     SUBROUTINE printw(string,value)
@@ -788,6 +883,6 @@
       CHARACTER(len=*), INTENT(in) :: value  !! (string) Value of the option.
       IF (mm_log) &
-      WRITE(*,'(a,a,a)') "warning: Parameter "//string//"not given... Using default value: "//value
-    END SUBROUTINE printw 
+        WRITE(*,'(a,a,a)') "warning: Parameter "//string//"not given... Using default value: "//value
+    END SUBROUTINE printw
   END FUNCTION mm_global_init_0
 
@@ -796,14 +891,15 @@
     !!
     !! See [[mm_globals(module):mm_global_init_0(function)]].
-    TYPE(cfgparser), INTENT(in) :: cfg  !! Configuration file path.
-    TYPE(error) :: err                  !! Error status of the function.
+    TYPE(cfgparser), INTENT(in) :: cfg
+    !! Configuration file path.
+    TYPE(error) :: err
+    !! Error status of the function.
     INTEGER                                           :: i
     TYPE(cfgparser)                                   :: spccfg
     CHARACTER(len=st_slen)                            :: spcpath
-    CHARACTER(len=:), ALLOCATABLE                     :: defmsg
     CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
     REAL(kind=mm_wp)                                  :: zfiamin,zfiamax
     LOGICAL                                           :: zwhp,zwhs,zwhc,zwcs,zwcn,znofia, &
-                                                         zwstom0,zwstom3
+      zwstom0,zwstom3
 
     err = noerror
@@ -856,7 +952,7 @@
       ! Gets species property file path
       err = cfg_get_value(cfg,'specie_cfg',spcpath) ; IF (err /= 0) RETURN
-      ! Reads species properties configuration file  
+      ! Reads species properties configuration file
       err = cfg_read_config(spccfg,trim(spcpath)) ; IF (err /= 0) RETURN
-      err = cfg_get_value(spccfg,"used_species",species) 
+      err = cfg_get_value(spccfg,"used_species",species)
       IF (err /= 0) THEN
         err = error("mm_global_init: cannot retrieve 'used_species' values",-1)
@@ -868,5 +964,5 @@
       !mm_spcname(1:mm_nesp) = species(:)
       DO i=1,mm_nesp
-        mm_spcname(i) = to_lower(species(i))
+        mm_spcname(i) = TRIM(species(i))
         IF (.NOT.cfg_has_section(spccfg,TRIM(mm_spcname(i)))) THEN
           err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1)
@@ -890,4 +986,5 @@
 
     ! MP2M Optional parameters
+    err = mm_check_opt(cfg_get_value(cfg,"debug",mm_debug),mm_debug,.false.,wlog=mm_log)
     err = mm_check_opt(cfg_get_value(cfg,"haze_production",mm_w_haze_prod),mm_w_haze_prod,zwhp,wlog=mm_log)
     err = mm_check_opt(cfg_get_value(cfg,"haze_sedimentation",mm_w_haze_sed),mm_w_haze_sed,zwhs,wlog=mm_log)
@@ -901,4 +998,19 @@
     err = mm_check_opt(cfg_get_value(cfg,"fiadero_max_ratio",mm_fiadero_max),mm_fiadero_max,zfiamax,wlog=mm_log)
 
+    err = mm_check_opt(cfg_get_value(cfg,"m0as_min",mm_m0as_min),mm_m0as_min,1e-10_mm_wp,wlog=mm_log)
+    err = mm_check_opt(cfg_get_value(cfg,"rcs_min",mm_rcs_min),mm_rcs_min,1e-9_mm_wp,wlog=mm_log)
+    err = mm_check_opt(cfg_get_value(cfg,"m0af_min",mm_m0af_min),mm_m0af_min,1e-10_mm_wp,wlog=mm_log)
+    err = mm_check_opt(cfg_get_value(cfg,"rcf_min",mm_rcf_min),mm_rcf_min,mm_rm,wlog=mm_log)
+    err = mm_check_opt(cfg_get_value(cfg,"m0n_min",mm_m0n_min),mm_m0n_min,1e-10_mm_wp,wlog=mm_log)
+
+
+    ! force fractal radius minimum threshold to monomer radius ^^
+    mm_rcf_min = MAX(mm_rm,mm_rcf_min)
+
+    ! compute m3 thresholds from user-defined thresholds.
+    mm_m3as_min  =  mm_m0as_min*mm_alpha_s(3._mm_wp) * mm_rcs_min**3._mm_wp
+    mm_m3af_min  =  mm_m0af_min*mm_alpha_f(3._mm_wp) * mm_rcf_min**3._mm_wp
+    mm_m3cld_min =  mm_m0n_min * (4._mm_wp * mm_pi / 3._mm_wp) * mm_drad_min**3._mm_wp
+
     err = noerror
     ! special check for settling velocity
@@ -911,8 +1023,8 @@
   FUNCTION mm_column_init(plev,zlev,play,zlay,temp) RESULT(err)
     !! Initialize vertical atmospheric fields.
-    !! 
+    !!
     !! This subroutine initializes vertical fields needed by the microphysics:
     !!
-    !! 1. Save reversed input field into "local" array 
+    !! 1. Save reversed input field into "local" array
     !! 2. Compute thicknesses layers and levels
     !! 3. Interpolate temperature at levels
@@ -922,5 +1034,5 @@
     !! @attention
     !! All the input vectors should be defined from __GROUND__ to __TOP__ of the atmosphere,
-    !! otherwise nasty things will occur in computations. 
+    !! otherwise nasty things will occur in computations.
     REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: plev !! Pressure levels (Pa).
     REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlev !! Altitude levels (m).
@@ -930,5 +1042,5 @@
     TYPE(error) :: err                                 !! Error status of the function.
     INTEGER :: i
-    mm_ini_col = .false.                          
+    mm_ini_col = .false.
     err = noerror
     IF (.NOT.mm_ini) THEN
@@ -980,6 +1092,6 @@
     ! Hydrostatic equilibrium
     mm_rhoair(1:mm_nla) = (mm_plev(2:mm_nle)-mm_plev(1:mm_nla)) / &
-                          (mm_effg(mm_zlay)*mm_dzlev)
-    mm_ini_col = .true.                          
+      (mm_effg(mm_zlay)*mm_dzlev)
+    mm_ini_col = .true.
     ! write out profiles (only if BOTH debug and log are enabled).
     IF (mm_log.AND.mm_debug) THEN
@@ -1003,15 +1115,15 @@
   FUNCTION mm_aerosols_init(m0aer_s,m3aer_s,m0aer_f,m3aer_f) RESULT(err)
     !! Initialize clouds tracers vertical grid.
-    !! 
-    !! The subroutine initializes aerosols microphysics tracers columns. It allocates variables if 
-    !! required and stores input vectors in reversed order. It also computes the characteristic radii 
-    !! of each mode. 
+    !!
+    !! The subroutine initializes aerosols microphysics tracers columns. It allocates variables if
+    !! required and stores input vectors in reversed order. It also computes the characteristic radii
+    !! of each mode.
     !! @note
-    !! All the input arguments should be defined from ground to top. 
+    !! All the input arguments should be defined from ground to top.
     !!
     !! @attention
     !! [[mm_globals(module):mm_global_init(interface)]] and [[mm_globals(module):mm_column_init(function)]]
     !! must have been called at least once before this method is called. Moreover, this method should be
-    !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the 
+    !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the
     !! vertical atmospheric structure.
     REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_s !! \(0^{th}\) order moment of the spherical mode (\(m^{-2}\)).
@@ -1063,14 +1175,9 @@
     mm_m0aer_f = m0aer_f(mm_nla:1:-1)/mm_dzlev(:)
     mm_m3aer_f = m3aer_f(mm_nla:1:-1)/mm_dzlev(:)
+
+    ! Setup threshold:
+    call mm_set_moments_thresholds()
+
     ! aerosols characteristic radii
-    ! il faudrait peut etre revoir la gestion des zeros ici et la 
-    ! remplacer par une valeur seuil des moments.
-    !
-    !-> JVO 19 : Done. Zero threshold set at espilon from dynamics on the
-    ! input moments in calmufi (safer than here). Might still be some unphysical
-    ! values after the dynamics near the threshold. Could be a could idea to add
-    ! a sanity check filtering too high values of radii.
-    !
-    ! TBD : Add a sanity check for high radii ????
     WHERE(mm_m3aer_s > 0._mm_wp .AND. mm_m0aer_s > 0._mm_wp)
       mm_rcs = mm_get_rcs(mm_m0aer_s,mm_m3aer_s)
@@ -1088,15 +1195,15 @@
   FUNCTION mm_clouds_init(m0ccn,m3ccn,m3ice,gazs) RESULT(err)
     !! Initialize clouds tracers vertical grid.
-    !! 
-    !! The subroutine initializes cloud microphysics tracers columns. It allocates variables if 
-    !! required and stores input vectors in reversed order. It also computes the mean drop radius 
+    !!
+    !! The subroutine initializes cloud microphysics tracers columns. It allocates variables if
+    !! required and stores input vectors in reversed order. It also computes the mean drop radius
     !! and density and allocates diagnostic vectors.
     !! @note
-    !! All the input arguments should be defined from ground to top. 
+    !! All the input arguments should be defined from ground to top.
     !!
     !! @attention
     !! [[mm_globals(module):mm_global_init(interface)]] and [[mm_globals(module):mm_column_init(function)]]
     !! must have been called at least once before this method is called. Moreover, this method should be
-    !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the 
+    !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the
     !! vertical atmospheric structure.
     REAL(kind=mm_wp), DIMENSION(:), INTENT(in)   :: m0ccn !! 0th order moment of the CCN distribution (\(m^{-2}\)).
@@ -1121,5 +1228,5 @@
     !  Actually, mm_nla should always initialized the first time mm_column_init is called, NOT HERE.
     IF (mm_nla < 0)  mm_nla  = SIZE(gazs,DIM=1)
-    ! Note: 
+    ! Note:
     !   here we could check that mm_nesp is the same size of gazs(DIM=2)
     !   Actually, mm_nesp should be always initialized in mm_global_init, NOT HERE.
@@ -1134,4 +1241,7 @@
     IF (.NOT.ALLOCATED(mm_drho))    ALLOCATE(mm_drho(mm_nla))
     ! Allocate memory for diagnostics
+    IF (.NOT.ALLOCATED(mm_ccn_w)) THEN
+      ALLOCATE(mm_ccn_w(mm_nla)) ; mm_ccn_w(:) = 0._mm_wp
+    ENDIF
     IF (.NOT.ALLOCATED(mm_ccn_flux)) THEN
       ALLOCATE(mm_ccn_flux(mm_nla)) ; mm_ccn_flux(:) = 0._mm_wp
@@ -1154,4 +1264,8 @@
       mm_gazs(:,i)  = gazs(mm_nla:1:-1,i)
     ENDDO
+
+    ! Setup threshold :
+    call mm_set_moments_cld_thresholds()
+
     ! drop mean radius
     call mm_cloud_properties(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_drad,mm_drho)
@@ -1162,4 +1276,6 @@
     !! Dump model global parameters on stdout.
     WRITE(*,'(a)')         "========= YAMMS PARAMETERS ============"
+    WRITE(*,'(a,a)')       "mm_fp_precision        : ", mm_wp_s
+    WRITE(*,'(a,L2)')      "mm_debug               : ", mm_debug
     WRITE(*,'(a,L2)')      "mm_w_haze_prod         : ", mm_w_haze_prod
     WRITE(*,'(a,ES14.7)')  "   mm_p_prod           : ", mm_p_prod
@@ -1168,5 +1284,5 @@
     WRITE(*,'(a,L2)')      "mm_w_haze_coag         : ", mm_w_haze_coag
     WRITE(*,'(a,I2.2)')    "   mm_coag_interactions: ", mm_coag_choice
-    WRITE(*,'(a,L2)')      "mm_w_haze_sed          : ", mm_w_haze_sed 
+    WRITE(*,'(a,L2)')      "mm_w_haze_sed          : ", mm_w_haze_sed
     WRITE(*,'(a,L2)')      "   mm_wsed_m0          : ", mm_wsed_m0
     WRITE(*,'(a,L2)')      "   mm_wsed_m3          : ", mm_wsed_m3
@@ -1178,4 +1294,15 @@
     WRITE(*,'(a,L2)')      "   mm_w_cloud_nucond   : ", mm_w_cloud_nucond
     WRITE(*,'(a)')         "---------------------------------------"
+    WRITE(*,'(a)')         "Thresholds spherical mode"
+    WRITE(*,'(a,ES14.7)')  "  mm_m0as_min          : ", mm_m0as_min
+    WRITE(*,'(a,ES14.7)')  "  mm_rcs_min           : ", mm_rcs_min
+    WRITE(*,'(a)')         "Thresholds fractal mode"
+    WRITE(*,'(a,ES14.7)')  "  mm_m0af_min          : ", mm_m0af_min
+    WRITE(*,'(a,ES14.7)')  "  mm_rcf_min           : ", mm_rcf_min
+    WRITE(*,'(a)')         "Thresholds clouds drop"
+    WRITE(*,'(a,ES14.7)')  "  mm_m0n_min           : ", mm_m0n_min
+    WRITE(*,'(a,ES14.7)')  "  mm_drad_min          : ", mm_drad_min
+    WRITE(*,'(a,ES14.7)')  "  mm_drad_max          : ", mm_drad_max
+    WRITE(*,'(a)')         "---------------------------------------"
     WRITE(*,'(a,ES14.7)')  "mm_dt                  : ", mm_dt
     IF (mm_nla > -1) THEN
@@ -1191,15 +1318,70 @@
   END SUBROUTINE mm_dump_parameters
 
+  SUBROUTINE mm_set_moments_thresholds()
+    !! Apply minimum threshold for the aerosols moments.
+    !!
+    !! The method resets moments (for both modes and orders, 0 and 3) values to zero if
+    !! their current value is below the minimum threholds.
+    !!
+    !! See also [[mm_globals(module):mm_m0as_min(variable)]], [[mm_globals(module):mm_rcs_min(variable)]],
+    !! [[mm_globals(module):mm_rcf_min(variable)]] and [[mm_globals(module):mm_m0as_min(variable)]].
+    INTEGER :: i
+    DO i=1,mm_nla
+      IF ((mm_m0aer_s(i) < mm_m0as_min) .OR. (mm_m3aer_s(i) < mm_m3as_min)) THEN
+        mm_m0aer_s(i) = 0._mm_wp ! mm_m0as_min
+        mm_m3aer_s(i) = 0._mm_wp ! mm_m0as_min * mm_rcs_min**3._mm_wp * mm_alpha_s(3._mm_wp)
+      ENDIF
+      IF ((mm_m0aer_f(i) < mm_m0af_min) .OR. (mm_m3aer_f(i) < mm_m3af_min)) THEN
+        mm_m0aer_f(i) = 0._mm_wp ! mm_m0af_min
+        mm_m3aer_f(i) = 0._mm_wp ! mm_m0af_min * mm_rcf_min**3._mm_wp * mm_alpha_f(3._mm_wp)
+      ENDIF
+    ENDDO
+  END SUBROUTINE mm_set_moments_thresholds
+
+  SUBROUTINE mm_set_moments_cld_thresholds()
+    !! Apply minimum threshold for the cloud drop moments.
+    !!
+    !! The method resets moments (for both modes and orders, 0 and 3) values to zero if
+    !! their current value is below the minimum threholds.
+    INTEGER :: i, j
+    REAL(kind=mm_wp) :: m3cld
+
+    DO i = 1, mm_nla
+      m3cld = mm_m3ccn(i)
+      DO j = 1, mm_nesp
+        m3cld = m3cld + mm_m3ice(i,j)
+      ENDDO
+
+      IF ((mm_m0ccn(i) < mm_m0n_min) .OR. (m3cld < mm_m3cld_min)) THEN
+        mm_m0ccn(i) = 0._mm_wp
+        mm_m3ccn(i) = 0._mm_wp
+        DO j = 1, mm_nesp
+          mm_m3ice(i,j) = 0._mm_wp
+        ENDDO
+      ENDIF
+    ENDDO
+  END SUBROUTINE mm_set_moments_cld_thresholds
+
+  ELEMENTAL SUBROUTINE mm_check_tendencies(v,dv)
+    !! Check that tendencies is not greater than value.
+    !!
+    !! the purpose of the subroutine is to update dvalue so that v+dv is not negative.
+    REAL(kind=mm_wp), INTENT(in)    :: v  !! Value to check.
+    REAL(kind=mm_wp), INTENT(inout) :: dv !! Value tendencies to check and update consequently.
+    REAL(kind=mm_wp), PARAMETER :: a = (epsilon(1._mm_wp)-1._mm_wp)
+    IF (v+dv < 0._mm_wp) THEN
+      dv = a*v
+    ENDIF
+  END SUBROUTINE mm_check_tendencies
+
   ELEMENTAL FUNCTION mm_get_rcs(m0,m3) RESULT(res)
     !! Get the characteristic radius for the spherical aerosols size distribution.
-    !! 
+    !!
     !! The method computes the characteristic radius of the size distribution law
     !! of the spherical aerosols mode according to its moments and its inter-moments
     !! relation.
-    REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment 
+    REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment
     REAL(kind=mm_wp), INTENT(in) :: m3 !! \(3^{rd}\) order moment
     REAL(kind=mm_wp) :: res            !! Radius
-    ! arbitrary: if there is no way to compute radius 
-    IF (m3 <= 0._mm_wp .OR. m0 <= 0._mm_wp) res = 1._mm_wp 
     res = (m3/m0/mm_alpha_s(3._mm_wp))**(1._mm_wp/3._mm_wp)
   END FUNCTION mm_get_rcs
@@ -1207,17 +1389,15 @@
   ELEMENTAL FUNCTION mm_get_rcf(m0,m3) RESULT(res)
     !! Get the characteristic radius for the fractal aerosols size distribution.
-    !! 
+    !!
     !! The method computes the characteristic radius of the size distribution law
     !! of the fractal aerosols mode according to its moments and its inter-moments
     !! relation.
-    REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment 
+    REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment
     REAL(kind=mm_wp), INTENT(in) :: m3 !! \(3^{rd}\) order moment
     REAL(kind=mm_wp) :: res            !! Radius
-    ! arbitrary: if there is no way to compute radius 
-    IF (m3 <= 0._mm_wp .OR. m0 <= 0._mm_wp) res = 1._mm_wp 
     res = (m3/m0/mm_alpha_f(3._mm_wp))**(1._mm_wp/3._mm_wp)
   END FUNCTION mm_get_rcf
 
-  ELEMENTAL FUNCTION mm_effg(z) RESULT(effg) 
+  ELEMENTAL FUNCTION mm_effg(z) RESULT(effg)
     !! Compute effective gravitational acceleration.
     REAL(kind=mm_wp), INTENT(in) :: z !! Altitude in meters
@@ -1226,5 +1406,5 @@
     IF (mm_use_effg) effg = effg * (mm_rpla/(mm_rpla+z))**2
     RETURN
-  END FUNCTION mm_effg 
+  END FUNCTION mm_effg
 
   !==================================
@@ -1237,39 +1417,43 @@
     !! The method computes the mean radius and mean density of cloud drops.
     !!
-    !! @bug 
-    !! A possible bug can happen because of threshold snippet. If __drad__ is greater than 
-    !! __drmax__ (== 100 microns) it is automatically set to __drmax__, but computation of 
+    !! @bug
+    !! A possible bug can happen because of threshold snippet. If __drad__ is greater than
+    !! __drmax__ (== 1e3 microns) it is automatically set to __drmax__, but computation of
     !! __drho__ remains unmodified. So __drho__ is not correct in that case.
     !!
-    !! @todo 
-    !! Fix the bug of the subroutine, but it is rather minor, since theoretically we do not 
+    !! @todo
+    !! Fix the bug of the subroutine, but it is rather minor, since theoretically we do not
     !! need the density of the drop.
     !!
-    !! @todo 
-    !! Think about a better implementation of thresholds, and get sure of their consequences in 
-    !! the other parts of the model. 
-    REAL(kind=mm_wp), INTENT(in)               :: m0ccn !! \(0^{th}\) order moment of the ccn 
-    REAL(kind=mm_wp), INTENT(in)               :: m3ccn !! \(3^{rd}\) order moment of the ccn 
+    !! @todo
+    !! Think about a better implementation of thresholds, and get sure of their consequences in
+    !! the other parts of the model.
+    REAL(kind=mm_wp), INTENT(in)               :: m0ccn !! \(0^{th}\) order moment of the ccn
+    REAL(kind=mm_wp), INTENT(in)               :: m3ccn !! \(3^{rd}\) order moment of the ccn
     REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ice !! \(3^{rd}\) order moments of each ice component
-    REAL(kind=mm_wp), INTENT(out)              :: drad  !! Output mean drop radius 
+    REAL(kind=mm_wp), INTENT(out)              :: drad  !! Output mean drop radius
     REAL(kind=mm_wp), INTENT(out), OPTIONAL    :: drho  !! Optional output mean drop density
-    REAL(kind=mm_wp)            :: vtot, wtot, ntot 
-    REAL(kind=mm_wp), PARAMETER :: threshold = 1.e-25_mm_wp,         &  
-                                       drmin = 1.e-10_mm_wp,         &
-                                       drmax = 1.e-4_mm_wp,          &
-                                      athird = 1._mm_wp/3._mm_wp,    &
-                                       pifac = 4._mm_wp/3._mm_wp*mm_pi 
+    REAL(kind=mm_wp)            :: Ntot, Vtot, Wtot
+    REAL(kind=mm_wp), PARAMETER :: athird = 1._mm_wp / 3._mm_wp
+    REAL(kind=mm_wp), PARAMETER :: pifac  = (4._mm_wp * mm_pi) / 3._mm_wp
+    
+    ! Set to zero :
     drad = 0._mm_wp
-    ntot = m0ccn
-    vtot = pifac*m3ccn+SUM(m3ice)
-    wtot = pifac*m3ccn*mm_rhoaer+SUM(m3ice*mm_xESPS(:)%rho)
-    IF (ntot <= threshold .OR. vtot <= 0._mm_wp) THEN
-      drad  = drmin
-      IF (PRESENT(drho)) drho  = mm_rhoaer 
-    ELSE
-      drad = (vtot/ntot/pifac)**athird
-      drad = MAX(MIN(drad,drmax),drmin)
-      IF (PRESENT(drho)) drho = wtot/vtot
-    ENDIF
+    IF (PRESENT(drho)) drho  = 0._mm_wp
+    
+    ! Initialization :
+    Ntot = m0ccn
+    Vtot = pifac * m3ccn + SUM(m3ice)
+    Wtot = pifac * ((m3ccn*mm_rhoaer) + SUM(m3ice*mm_xESPS(:)%rho))
+
+    IF (Ntot <= mm_m0n_min .OR. Vtot <= mm_m3cld_min) THEN
+      drad = mm_drad_min
+      IF (PRESENT(drho)) drho = mm_rhoaer
+    ELSE
+      drad = (Vtot / (pifac*Ntot))**athird
+      drad = MAX(MIN(drad,mm_drad_max),mm_drad_min)
+      IF (PRESENT(drho)) drho = Wtot / Vtot
+    ENDIF
+
     RETURN
   END SUBROUTINE cldprop_sc
@@ -1286,33 +1470,14 @@
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:)           :: drad  !! Output mean drop radius.
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: drho  !! Optional output mean drop density.
-    INTEGER                                     :: i,ns
-    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: vtot,wtot,ntot,rho 
-    REAL(kind=mm_wp), PARAMETER                 :: threshold = 1.e-25_mm_wp,         &  
-                                                       drmin = 1.e-10_mm_wp,         &
-                                                       drmax = 1.e-4_mm_wp,          &
-                                                      athird = 1._mm_wp/3._mm_wp,    &
-                                                       pifac = 4._mm_wp/3._mm_wp*mm_pi
-                                                     
-    ns = SIZE(m0ccn) ; ALLOCATE(vtot(ns),wtot(ns),ntot(ns),rho(ns))
-    drad = 0._mm_wp
-    ntot = m0ccn
-    vtot = pifac*m3ccn+SUM(m3ice,DIM=2)
-    wtot = pifac*m3ccn*mm_rhoaer
-    DO i=1,SIZE(m3ice,DIM=2)
-      wtot = wtot+m3ice(:,i)*mm_xESPS(i)%rho
-    ENDDO
-    WHERE(ntot <= threshold .OR. vtot <= 0._mm_wp) 
-      drad  = drmin
-      rho = mm_rhoaer
-    ELSEWHERE 
-      drad = (vtot/ntot/pifac)**athird
-      drad = MAX(MIN(drad,drmax),drmin)
-      rho = wtot/vtot
-    END WHERE
-    IF (PRESENT(drho)) drho  = rho
+    INTEGER :: i
+    IF (PRESENT(drho)) THEN
+      DO i = 1, SIZE(m0ccn) ; call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i),drho(i)) ; ENDDO
+    ELSE
+      DO i = 1, SIZE(m0ccn) ; call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i)) ; ENDDO
+    ENDIF
     RETURN
   END SUBROUTINE cldprop_ve
 
-  ! For configuration file (requires fccp library).
+  ! For configuration file (requires swift library).
 
   FUNCTION read_esp(parser,sec,pp) RESULT (err)
@@ -1322,23 +1487,23 @@
     TYPE(mm_esp), INTENT(out)     :: pp     !! [[mm_globals(module):mm_esp(type)]] object that stores the parameters.
     TYPE(error)                   :: err    !! Error status of the function.
-    err = cfg_get_value(parser,TRIM(sec)//'name',pp%name)       ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'mas',pp%mas)         ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'vol',pp%vol)         ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'ray',pp%ray)         ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'mas',pp%mas)         ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'vol',pp%vol)         ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'ray',pp%ray)         ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'masmol',pp%masmol)   ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'rho',pp%rho)         ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'tc',pp%tc)           ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'pc',pp%pc)           ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'tb',pp%tb)           ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'w',pp%w)             ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'a_sat',pp%a_sat)     ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'b_sat',pp%b_sat)     ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'c_sat',pp%c_sat)     ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'d_sat',pp%d_sat)     ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'mteta',pp%mteta)     ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'tx_prod',pp%tx_prod) ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/name',pp%name)       ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/mas',pp%mas)         ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/vol',pp%vol)         ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/ray',pp%ray)         ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/mas',pp%mas)         ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/vol',pp%vol)         ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/ray',pp%ray)         ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/masmol',pp%masmol)   ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/rho',pp%rho)         ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/tc',pp%tc)           ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/pc',pp%pc)           ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/tb',pp%tb)           ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/w',pp%w)             ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/a_sat',pp%a_sat)     ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/b_sat',pp%b_sat)     ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/c_sat',pp%c_sat)     ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/d_sat',pp%d_sat)     ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/mteta',pp%mteta)     ; IF (err /= 0) RETURN
+    err = cfg_get_value(parser,TRIM(sec)//'/tx_prod',pp%tx_prod) ; IF (err /= 0) RETURN
     RETURN
   END FUNCTION read_esp
@@ -1346,5 +1511,5 @@
   ! =========================================================================
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !                CONFIGURATION PARSER checking methods 
+  !                CONFIGURATION PARSER checking methods
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ! =========================================================================
@@ -1352,6 +1517,6 @@
   FUNCTION check_r1(err,var,def,wlog) RESULT(ret)
     !! Check an option value (float).
-    !! 
-    !! The method checks an option value and optionally set a default value, __def__ to initialize 
+    !!
+    !! The method checks an option value and optionally set a default value, __def__ to initialize
     !! __var__ on error if given.
     TYPE(error), INTENT(in)                :: err  !! Error object from value getter.
@@ -1360,5 +1525,5 @@
     LOGICAL, INTENT(in), OPTIONAL          :: wlog !! .true. to print warning/error message.
     TYPE(error) :: ret                             !! Input error.
-    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: ' 
+    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
     LOGICAL                     :: zlog
     ret = err
@@ -1376,6 +1541,6 @@
   FUNCTION check_l1(err,var,def,wlog) RESULT(ret)
     !! Check an option value (logical).
-    !! 
-    !! The method checks an option value and optionally set a default value, __def__ to initialize 
+    !!
+    !! The method checks an option value and optionally set a default value, __def__ to initialize
     !! __var__ on error if given.
     TYPE(error), INTENT(in)       :: err  !! Error object from value getter.
@@ -1384,8 +1549,8 @@
     LOGICAL, INTENT(in), OPTIONAL :: wlog !! .true. to print warning/error message.
     TYPE(error) :: ret                    !! Input error.
-    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: ' 
+    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
     LOGICAL                     :: zlog
     ret = err
-     zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
+    zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
     IF (err == 0) RETURN
     IF (PRESENT(def)) THEN
@@ -1400,6 +1565,6 @@
   FUNCTION check_i1(err,var,def,wlog) RESULT(ret)
     !! Check an option value (integer).
-    !! 
-    !! The method checks an option value and optionally set a default value, __def__ to initialize 
+    !!
+    !! The method checks an option value and optionally set a default value, __def__ to initialize
     !! __var__ on error if given.
     TYPE(error), INTENT(in)       :: err  !! Error object from value getter.
@@ -1411,5 +1576,5 @@
     LOGICAL                     :: zlog
     ret = err
-     zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
+    zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
     IF (err == 0) RETURN
     IF (PRESENT(def)) THEN
@@ -1424,6 +1589,6 @@
   FUNCTION check_s1(err,var,def,wlog) RESULT(ret)
     !! Check an option value (string).
-    !! 
-    !! The method checks an option value and optionally set a default value, __def__ to initialize 
+    !!
+    !! The method checks an option value and optionally set a default value, __def__ to initialize
     !! __var__ on error if given.
     TYPE(error), INTENT(in)                      :: err  !! Error object from value getter.
@@ -1432,7 +1597,7 @@
     LOGICAL, INTENT(in), OPTIONAL                :: wlog !! .true. to print warning/error message.
     TYPE(error) :: ret                                   !! Input error.
-    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: ' 
+    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
     LOGICAL                     :: zlog
-    ret = err 
+    ret = err
     zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
     IF (err == 0) RETURN
