source: trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_globals.F90 @ 3590

MODULE MP2M_GLOBALS
    !============================================================================
    !
    !     Purpose
    !     -------
    !     Parameters and global variables module.
    !     Defines and initialize all the parameters and global variables that are
    !     common to all other modules of the library.
    !
    !     It is separated in two parts :
    !        - Main parameters and global saved variables. Most of these variables should be initialized
    !          once and should hold the same value during run-time. These variables are completly public and
    !          initialized by [[mm_globals(module):mm_global_init(interface)]] 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.
    !
    !     @note
    !     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__.
    !
    !     Global variables overview:
    !        - Protected variables
    !        - Control flags
    !        - Related free parameters
    !        - Physical constants
    !        - Free parameters
    !        - Miscellaneous parameters
    !        - Vertical structure part
    !
    !     The module also contains thirteen methods:
    !        - mm_global_init_0
    !        - mm_global_init_1
    !        - mm_column_init
    !        - mm_aerosols_init
    !        - mm_alpha_s, mm_alpha_f
    !        - mm_effg
    !        - mm_set_moments_thresholds
    !        - mm_get_rcs, mm_get_rcf
    !        - mm_dump_parameters
    !        - check_r1, check_l1, check_i1, check_s1
    !
    !     Authors
    !     -------
    !     B. de Batz de Trenquelléon, J. Burgalat (11/2024)
    !
    !============================================================================
    USE MP2M_MPREC
    ! From lint
    USE LINT_DATASETS
    ! From swift
    USE SWIFT_CFGPARSE
    USE SWIFT_STRING_OP
    USE SWIFT_ERRORS
    IMPLICIT NONE

    PUBLIC

    PRIVATE :: check_r1,check_i1,check_l1,check_s1

    ! ~~~~~~~~~~~~~~~~~~~
    ! Protected variables
    ! ~~~~~~~~~~~~~~~~~~~
    ! 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
    ! Model parameters (mm_global_init)
    PROTECTED :: mm_dt,mm_rhoaer,mm_df,mm_rm,mm_haze_prod_pCH4,mm_p_prod,mm_rc_prod,mm_tx_prod,mm_rpla,mm_g0,mm_rb2ra
    ! Atmospheric vertical structure (mm_column_init)
    PROTECTED :: mm_nla,mm_nle,mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay
    ! Moments parameters (mm_aerosols_init)
    PROTECTED :: mm_m0aer_s, mm_m3aer_s, mm_m0aer_f, mm_m3aer_f
    ! Moments parameters (derived, are updated with moments parameters)
    PROTECTED :: mm_rcs, mm_rcf
    ! Thresholds parameters
    PROTECTED :: mm_m0as_min, mm_m3as_min, mm_rcs_min, mm_m0af_min, mm_m3af_min, mm_rcf_min

    ! ~~~~~~~~~~~~~
    ! Control flags
    ! ~~~~~~~~~~~~~
    ! Enable/Disable haze production.
    LOGICAL, SAVE :: mm_w_haze_prod = .true.
    ! Enable/Disable haze production from CH4 photolysis.
    LOGICAL, SAVE :: mm_haze_prod_pCH4 = .true.
    ! Enable/Disable haze sedimentation.
    LOGICAL, SAVE :: mm_w_haze_sed = .true.
    ! Enable/Disable haze coagulation.
    LOGICAL, SAVE :: mm_w_haze_coag = .true.
    ! Force all aerosols moments to fall at M0 settling velocity.
    LOGICAL, SAVE :: mm_wsed_m0 = .false.
    ! Force all aerosols moments to fall at M3 settling velocity.
    LOGICAL, SAVE :: mm_wsed_m3 = .false.
    ! Enable/Disable spherical probability transfert.
    LOGICAL, SAVE :: mm_w_ps2s = .true.
    ! Enable/Disable aerosol electric charge correction.
    LOGICAL, SAVE :: mm_w_qe   = .true.

    ! Enable/Disable QnD debug mode (can be used for devel).
    LOGICAL, SAVE :: mm_debug = .false.
    ! Enable/Disable log mode (for configuration only).
    LOGICAL, SAVE :: mm_log   = .false.
    ! Enable/Disable effective G for computations.
    LOGICAL, SAVE :: mm_use_effg = .true.

    ! Initialization control flag [[mm_globals(module):mm_global_init(interface)]].
    LOGICAL, PUBLIC, SAVE :: mm_ini     = .false.
    ! Initialization control flag [[mm_globals(module):mm_column_init(function)]].
    LOGICAL, PUBLIC, SAVE :: mm_ini_col = .false.
    ! Initialization control flag [[mm_globals(module):mm_aerosols_init(function)]].
    LOGICAL, PUBLIC, SAVE :: mm_ini_aer = .false.
  
    ! ~~~~~~~~~~~~~~~~~~~~~~~
    ! Related free parameters
    ! ~~~~~~~~~~~~~~~~~~~~~~~
    ! No mode interaction for coagulation (i.e. no coagulation at all).
    INTEGER, PARAMETER :: mm_coag_no = 0
    ! SS mode interaction for coagulation.
    INTEGER, PARAMETER :: mm_coag_ss = 1
    ! SF mode interaction for coagulation.
    INTEGER, PARAMETER :: mm_coag_sf = 2
    ! FF mode interaction for coagulation.
    INTEGER, PARAMETER :: mm_coag_ff = 4
    ! Default interactions to activate (all by default).
    INTEGER, SAVE      :: mm_coag_choice = mm_coag_ss+mm_coag_sf+mm_coag_ff

    ! ~~~~~~~~~~~~~~~~~~
    ! Physical constants
    ! ~~~~~~~~~~~~~~~~~~
    ! Pi number.
    REAL(kind=mm_wp), PARAMETER :: mm_pi = 4._mm_wp*atan(1._mm_wp)
    ! Avogadro number (mol-1).
    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
    ! Perfect gas constant (J.mol-1.K-1).
    REAL(kind=mm_wp), PARAMETER :: mm_rgas = mm_kboltz * mm_navo
    ! Approximated slip-flow correction coefficient.
    REAL(kind=mm_wp), PARAMETER :: mm_akn = 1.591_mm_wp

    ! ~~~~~~~~~~~~~~~
    ! Free parameters
    ! ~~~~~~~~~~~~~~~
    ! Spherical aerosol production pressure level (Pa).
    REAL(kind=mm_wp), SAVE :: mm_p_prod = 1.e-2_mm_wp
    ! Spherical aerosol production rate (kg.m-2.s-1).
    REAL(kind=mm_wp), SAVE :: mm_tx_prod = 9.8e-14_mm_wp
    ! Spherical aerosol equivalent radius production (m)
    REAL(kind=mm_wp), SAVE :: mm_rc_prod = 1.e-9_mm_wp
    ! Monomer radius (m).
    REAL(kind=mm_wp), SAVE :: mm_rm = 1.e-8_mm_wp  
    ! Fractal dimension of fractal aerosols.
    REAL(kind=mm_wp), SAVE :: mm_df = 2._mm_wp
    ! Aerosol density (kg.m-3).
    REAL(kind=mm_wp), SAVE :: mm_rhoaer = 8.e2_mm_wp

    ! Total number of aerosols minimum threshold for the spherical mode.
    REAL(kind=mm_wp), SAVE :: mm_m0as_min = 1.e-8_mm_wp
    ! Total volume of aerosols minimum threshold for the spherical mode.
    REAL(kind=mm_wp), SAVE :: mm_m3as_min = 1.e-35_mm_wp
    ! Characteristic radius minimum threshold for the spherical mode.
    REAL(kind=mm_wp), SAVE :: mm_rcs_min = 1.e-9_mm_wp
    ! Total number of aerosols minimum threshold for the fractal mode.
    REAL(kind=mm_wp), SAVE :: mm_m0af_min = 1.e-8_mm_wp
    ! Total volume of aerosols minimum threshold for the fractal mode.
    REAL(kind=mm_wp), SAVE :: mm_m3af_min = 1.e-35_mm_wp
    ! Characteristic radius minimum threshold for the fractal mode.
    REAL(kind=mm_wp), SAVE :: mm_rcf_min = 1.e-9_mm_wp
  
    ! Planet radius (m) and gravity acceleration (m.s-2).
    ! WARNING: initialization for Pluto.
    REAL(kind=mm_wp), SAVE :: mm_rpla     = 1187000._mm_wp
    REAL(kind=mm_wp), SAVE :: mm_g0       = 0.617_mm_wp 
    ! Air molecules mean radius (m), molar mass (kg.mol-1), and molecular mass (kg).
    ! WARNING: initialization for N2.
    REAL(kind=mm_wp), SAVE :: mm_air_rad  = 1.75e-10_mm_wp
    REAL(kind=mm_wp), SAVE :: mm_air_mmol = 28.e-3_mm_wp
    REAL(kind=mm_wp), SAVE :: mm_air_mmas = 28.e-3_mm_wp / mm_navo
    ! Microphysical time-step (s).
    REAL(kind=mm_wp), SAVE :: mm_dt       = 180._mm_wp
    
    ! ~~~~~~~~~~~~~~~~~~~~~~~~
    ! Miscellaneous parameters
    ! ~~~~~~~~~~~~~~~~~~~~~~~~
    
    ! Bulk to apparent radius
    ! ~~~~~~~~~~~~~~~~~~~~~~~
    ! Bulk to apparent radius conversion pre-factor (m^X).
    !
    ! With r_a = r_b^(3/Df) . r_m^((Df-3)/(Df))
    ! Then rb2ra = r_m^((Df-3)/(Df))
    REAL(kind=mm_wp), SAVE :: mm_rb2ra = 1._mm_wp

    ! Inter-moment relation
    ! ~~~~~~~~~~~~~~~~~~~~~
    ! Alpha function parameters. 
    ! It stores the parameters of the inter-moments relation functions.
    !
    ! The inter-moments relation function is represented by the sum of exponential quadratic expressions:
    ! alpha(k) = Sum_{i=1}^{n} exp(a_i.k^2 + bi.k^2 + c_i) 
    TYPE, PUBLIC :: aprm
      ! Quadratic coefficients of the quadratic expressions.
      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a
      ! Linear coefficients of the quadratic expressions.
      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: b
      ! Free term of the quadratic expressions.
      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: c
    END TYPE
    ! Inter-moment relation set of parameters for the spherical mode.
    TYPE(aprm), PUBLIC, SAVE :: mm_asp
    ! Inter-moment relation set of parameters for the fractal mode.
    TYPE(aprm), PUBLIC, SAVE :: mm_afp

    ! Transfert probabilities (S --> F)
    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    ! Data set for linear interpolation of transfert probability (M0/CO).
    TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pco0p
    ! Data set for linear interpolation of transfert probability (M3/CO).
    TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pco3p
    ! Data set for linear interpolation of transfert probability (M0/FM).
    TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pfm0p
    ! Data set for linear interpolation of transfert probability (M3/FM). 
    TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pfm3p
  
    ! Mean electric correction
    ! ~~~~~~~~~~~~~~~~~~~~~~~~
    ! Data set for Q_SF(M0).
    TYPE(dset2d), PUBLIC, SAVE, TARGET             :: mm_qbsf0
    ! Extended values of [[mm_gcm(module):mm_qbsf0(variable)]] dataset.
    REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mm_qbsf0_e
    ! Data set for Q_SF(M3).
    TYPE(dset2d), PUBLIC, SAVE, TARGET             :: mm_qbsf3
    ! Extended values of [[mm_gcm(module):mm_qbsf3(variable)]] dataset.
    REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mm_qbsf3_e
    ! Data set for Q_FF(M0).
    TYPE(dset2d), PUBLIC, SAVE, TARGET             :: mm_qbff0
    ! Extended values of [[mm_gcm(module):mm_qbff0(variable)]] dataset.
    REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mm_qbff0_e

    ! btk coefficients
    ! ~~~~~~~~~~~~~~~~
    ! Coefficients for Free-molecular regime kernel approximation (b_0(t)).
    REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(5) :: mm_bt0 = (/1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp/)
    ! Coefficients for Free-molecular regime kernel approximation (b_3(t)).
    REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(5) :: mm_bt3 = (/1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp/)

    ! ~~~~~~~~~~~~~~~~~~~~~~~
    ! Vertical structure part
    ! ~~~~~~~~~~~~~~~~~~~~~~~
    ! Number of vertical layers.
    INTEGER, SAVE :: mm_nla = -1
    ! Number of vertical levels.
    INTEGER, SAVE :: mm_nle = -1
  
    ! Altitude layers (m).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlay
    ! Altitude levels (m).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlev
    !> Pressure layers (Pa).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_play
    ! Pressure levels (Pa).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_plev
    ! Temperature vertical profile (K).
    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
    ! Temperature vertical profil at interfaces (K).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_btemp

    ! Atmospheric levels thickness (m).
    ! @note: __mm_dzlev__ is defined on the total number of layers and actually
    ! corresponds to the thickness of a given layer.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlev
    ! Atmospheric layers thickness (m).
    ! @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.
    ! This value corresponds to the thickness between the center of the
    ! __GROUND__ layer and below the surface.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlay

    ! Spherical mode - 0th order moment (m-3).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_s
    ! Spherical mode - 3rd order moment (m3.m-3).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_s
    ! Spherical mode - characteristic radius (m).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcs
    ! Fractal mode - 0th order moment (m-3).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_f
    ! Fractal mode - 3rd order moment (m3.m-3).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_f
    ! Fractal mode - characteristic radius (m).
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcf

    ! Spherical aerosol precipitation (kg.m-2).
    REAL(kind=mm_wp), SAVE :: mm_aers_prec = 0._mm_wp
    ! Fractal aerosol precipitation (kg.m-2).
    REAL(kind=mm_wp), SAVE :: mm_aerf_prec = 0._mm_wp
  
    ! Spherical mode (M0) settling velocity (m.s-1).
    ! @note: This variable is always negative.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0as_vsed
    ! Spherical mode (M3) settling velocity (m.s-1).
    ! @note: This variable is always negative.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3as_vsed
    ! Fractal mode (M0) settling velocity (m.s-1).
    ! @note: This variable is always negative.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0af_vsed
    ! Fractal mode (M3) settling velocity (m.s-1).
    ! @note: This variable is always negative.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3af_vsed

    ! Spherical aerosol mass fluxes (kg.m-2.s-1).
    ! @note: This variable is always negative.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_s_flux
    ! Fractal aerosol mass fluxes (kg.m-2.s-1).
    ! @note: This variable is always negative.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_f_flux

    ! All variables related to column computations should be private to each thread
    !$OMP THREADPRIVATE(mm_ini_col,mm_ini_aer)
    !$OMP THREADPRIVATE(mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay)
    !$OMP THREADPRIVATE(mm_m0aer_s,mm_m3aer_s,mm_m0aer_f,mm_m3aer_f)
    !$OMP THREADPRIVATE(mm_rcs,mm_rcf)
    !$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed)
    !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux)
    !$OMP THREADPRIVATE(mm_m0as_min,mm_m3as_min,mm_rcs_min,mm_m0af_min,mm_m3af_min,mm_rcf_min)
    !$OMP THREADPRIVATE(mm_nla,mm_nle)

    ! Interface to global initialization.
    ! The method performs the global initialization of the model.
    INTERFACE mm_global_init
      MODULE PROCEDURE mm_global_init_0,mm_global_init_1
    END INTERFACE mm_global_init
  
    ! Check an option from the configuration system.
    ! The method checks for an option in the configuration system and optionally
    ! 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.
    INTERFACE mm_check_opt
      MODULE PROCEDURE check_r1,check_i1,check_l1,check_s1
    END INTERFACE mm_check_opt


    CONTAINS


    !============================================================================
    !                      INITIALIZATION METHODS
    !============================================================================

    FUNCTION mm_global_init_0(dt,df,rm,rho_aer,haze_prod_pCH4,p_prod,tx_prod,rc_prod,&
                              rplanet,g0,air_rad,air_mmol,                           &
                              coag_interactions,w_haze_prod,w_haze_sed,w_haze_coag,  &
                              force_wsed_to_m0,force_wsed_to_m3,                     &
                              m0as_min,rcs_min,m0af_min,rcf_min,debug) RESULT(err)
      !! Initialize global parameters of the model.
      !!
      !! The function initializes all the global parameters of the model from direct input.
      !! Store input values in global variables
      !!
      !! @note
      !! If the method fails to initialize parameters the model should be aborted as the
      !! global variables of the model will not be correctly setup.
      !!
      ! Microphysical timestep (s).
      REAL(kind=mm_wp), INTENT(in) :: dt
      ! Fractal dimension of fractal aerosol (-).
      REAL(kind=mm_wp), INTENT(in) :: df
      ! Monomer radius (m).
      REAL(kind=mm_wp), INTENT(in) :: rm
      ! Aerosol density(kg.m-3).
      REAL(kind=mm_wp), INTENT(in) :: rho_aer
      ! Enable/Disable production from CH4 photolysis.
      LOGICAL, INTENT(in)          :: haze_prod_pCH4
      ! Aerosol production pressure level (Pa).
      REAL(kind=mm_wp), INTENT(in) :: p_prod
      ! Spherical aerosol production rate kg.m-2.s-1).
      REAL(kind=mm_wp), INTENT(in) :: tx_prod
      ! Spherical aerosol characteristic radius of production (m).
      REAL(kind=mm_wp), INTENT(in) :: rc_prod
      ! Planet radius (m).
      REAL(kind=mm_wp), INTENT(in) :: rplanet
      ! Planet gravity acceleration at ground level (m.s-2).
      REAL(kind=mm_wp), INTENT(in) :: g0
      ! Radius of air molecules (m).
      REAL(kind=mm_wp), INTENT(in) :: air_rad
      ! Molar mass of air molecules (kg.mol-1).
      REAL(kind=mm_wp), INTENT(in) :: air_mmol

      ! Coagulation interactions process control flag.
      INTEGER, INTENT(in) :: coag_interactions
      ! Haze production process control flag.
      LOGICAL, INTENT(in) :: w_haze_prod
      ! Haze sedimentation process control flag.
      LOGICAL, INTENT(in) :: w_haze_sed
      ! Haze coagulation process control flag.
      LOGICAL, INTENT(in) :: w_haze_coag
      ! Force __all__ aerosols moments to fall at M0 settling velocity.
      LOGICAL, INTENT(in) :: force_wsed_to_m0
      ! Force __all__ aerosols moments to fall at M3 settling velocity
      LOGICAL, INTENT(in) :: force_wsed_to_m3

      ! Minimum threshold for M0 of the spherical mode (m-3).
      REAL(kind=mm_wp), INTENT(in) :: m0as_min
      ! Minimum threshold for the characteristic radius of the spherical mode (m).
      REAL(kind=mm_wp), INTENT(in) :: rcs_min
      ! Minimum threshold for M0 of the factal mode (m-3).
      REAL(kind=mm_wp), INTENT(in) :: m0af_min
      ! Minimum threshold for the characteristic radius of the fractal mode (m).
      REAL(kind=mm_wp), INTENT(in) :: rcf_min

      ! Debug mode control flag.
      LOGICAL, INTENT(in) :: debug
      
      ! Error status of the function.
      TYPE(error) :: err

      err = noerror

      ! Sanity check:
      IF (mm_ini) THEN
        err = error("mm_global_init: YAMMS global initialization already performed !",-1)
        RETURN
      ENDIF
  
      ! Free parameters:
      mm_dt             = dt
      mm_df             = df
      mm_rm             = rm
      mm_rhoaer         = rho_aer
      mm_haze_prod_pCH4 = haze_prod_pCH4
      mm_p_prod         = p_prod
      mm_tx_prod        = tx_prod
      mm_rc_prod        = rc_prod
      mm_rpla           = rplanet
      mm_g0             = g0
      mm_air_rad        = air_rad
      mm_air_mmol       = air_mmol
      mm_air_mmas       = air_mmol / mm_navo
      
      ! Microphysical processes:
      mm_coag_choice = coag_interactions
      IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN
        err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1)
        RETURN
      ENDIF
      mm_w_haze_prod = w_haze_prod
      mm_w_haze_sed  = w_haze_sed
      mm_w_haze_coag = w_haze_coag
      mm_wsed_m0     = force_wsed_to_m0
      mm_wsed_m3     = force_wsed_to_m3
      
      ! Moment threshold flags:
      mm_m0as_min = MAX(0._mm_wp,m0as_min)
      mm_rcs_min  = MAX(1.e-10_mm_wp,rcs_min)
      mm_m0af_min = MAX(0._mm_wp,m0af_min)
      mm_rcf_min  = MAX(mm_rm,rcf_min)
      
      ! Debug mode:
      mm_debug = debug
      
      ! Computes 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

      ! Computes conversion factor for bulk to apparent radius:
      mm_rb2ra = mm_rm**((mm_df-3._mm_wp)/mm_df)

      ! Sanity check for settling velocity:
      IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN
        err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1)
      ENDIF

      ! End of initialization:
      mm_ini = err == noerror

    END FUNCTION mm_global_init_0


    FUNCTION mm_global_init_1(cfg) RESULT(err)
      !! Set global configuration from a configuration file.
      !!
      !! @note:
      !! See [[mm_globals(module):mm_global_init_0(function)]].
      !!
      
      ! Configuration file path.
      TYPE(cfgparser), INTENT(in) :: cfg
      ! Error status of the function.
      TYPE(error) :: err
      
      err = noerror
      
      ! Sanity check:
      IF (mm_ini) THEN
        err = error("mm_global_init: YAMMS global initialization already performed !",-1)
        RETURN
      ENDIF
  
      ! Free parameters:
      err = mm_check_opt(cfg_get_value(cfg,"timestep",mm_dt),mm_dt,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"df",mm_df),mm_df,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"rm",mm_rm),mm_rm,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"rho_aer",mm_rhoaer),mm_rhoaer,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"call_haze_prod",mm_haze_prod_pCH4),mm_haze_prod_pCH4,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"p_prod",mm_p_prod),mm_p_prod,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"tx_prod",mm_tx_prod),mm_tx_prod,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"rc_prod",mm_rc_prod),mm_rc_prod,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"planet_radius",mm_rpla),mm_rpla,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"g0",mm_g0),mm_g0,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"air_radius",mm_air_rad),mm_air_rad,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"air_molarmass",mm_air_mmol),mm_air_mmol,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"air_molecularmass",mm_air_mmas),mm_air_mmas,wlog=mm_log)
      IF (err/=0) RETURN

      ! Microphysical processes:
      err = mm_check_opt(cfg_get_value(cfg,"haze_coag_interactions",mm_coag_choice),mm_coag_choice,wlog=mm_log)
      IF (err/=0) RETURN
      IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN
        err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1)
        RETURN
      ENDIF
      err = mm_check_opt(cfg_get_value(cfg,"haze_production",mm_w_haze_prod),mm_w_haze_prod,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"haze_sedimentation",mm_w_haze_sed),mm_w_haze_sed,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"haze_coagulation",mm_w_haze_coag),mm_w_haze_coag,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"wsed_m0",mm_wsed_m0),mm_wsed_m0,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"wsed_m3",mm_wsed_m3),mm_wsed_m3,wlog=mm_log)
      IF (err/=0) RETURN
      
      ! Moment threshold flags:
      err = mm_check_opt(cfg_get_value(cfg,"m0as_min",mm_m0as_min),mm_m0as_min,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"rcs_min",mm_rcs_min),mm_rcs_min,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"m0af_min",mm_m0af_min),mm_m0af_min,wlog=mm_log)
      IF (err/=0) RETURN
      err = mm_check_opt(cfg_get_value(cfg,"rcf_min",mm_rcf_min),mm_rcf_min,wlog=mm_log)
      IF (err/=0) RETURN
      
      ! Debug mode:
      err = mm_check_opt(cfg_get_value(cfg,"debug",mm_debug),mm_debug,wlog=mm_log)
      IF (err/=0) RETURN
      
      ! Computes M3 thresholds from user-defined thresholds:
      mm_m0as_min = MAX(0._mm_wp,mm_m0as_min)
      mm_rcs_min  = MAX(1.e-10_mm_wp,mm_rcs_min)
      mm_m0af_min = MAX(0._mm_wp,mm_m0af_min)
      mm_rcf_min  = MAX(mm_rm,mm_rcf_min)
      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

      ! Computes conversion factor for bulk to apparent radius:
      mm_rb2ra = mm_rm**((mm_df-3._mm_wp)/mm_df)

      ! Sanity check for settling velocity:
      IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN
        err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1)
      ENDIF

      ! End of initialization:
      mm_ini = err == noerror

    END FUNCTION mm_global_init_1


    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
      !!    2. Compute thicknesses layers and levels
      !!    3. Interpolate temperature at levels
      !!
      !! @warning
      !! The method should be called whenever the vertical structure of the atmosphere changes.
      !! All the input vectors should be defined from __GROUND__ to __TOP__ of the atmosphere,
      !! 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).
      REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: play ! Pressure layers (Pa).
      REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlay ! Altitude at the center of each layer (m).
      REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: temp ! Temperature at the center of each layer (K).
      TYPE(error) :: err                                 ! Error status of the function.
      INTEGER :: i

      err = noerror
      mm_ini_col = .false.
      
      ! Global initialization must be done before:
      IF (.NOT.mm_ini) THEN
        err = error("mm_column_init: Global initialization not done yet",-1)
        RETURN
      ENDIF

      ! Check number of vertical layers:
      IF (mm_nla < 0) THEN
        mm_nla = SIZE(play)
      ELSE
        IF (mm_nla /= SIZE(play)) THEN
          err = error("mm_column_init: mm_nla cannot be modified dynamically within the run",-1)
          RETURN
        ENDIF
      ENDIF
      
      ! Check number of vertical levels:
      IF (mm_nle < 0) THEN
        mm_nle = SIZE(plev)
      ELSE
        IF (mm_nle /= SIZE(plev)) THEN
          err = error("mm_column_init: mm_nle cannot be modified dynamically within the run",-1)
          RETURN
        ENDIF
      ENDIF

      ! Sanity check:
      IF (mm_nla+1 /= mm_nle) THEN
        err = error("mm_column_init: Inconsistent number of layers/levels",-1)
        RETURN
      ENDIF

      ! Allocates if required:
      IF (.NOT.ALLOCATED(mm_plev))   ALLOCATE(mm_plev(mm_nle))
      IF (.NOT.ALLOCATED(mm_zlev))   ALLOCATE(mm_zlev(mm_nle))
      IF (.NOT.ALLOCATED(mm_play))   ALLOCATE(mm_play(mm_nla))
      IF (.NOT.ALLOCATED(mm_zlay))   ALLOCATE(mm_zlay(mm_nla))
      IF (.NOT.ALLOCATED(mm_temp))   ALLOCATE(mm_temp(mm_nla))
      IF (.NOT.ALLOCATED(mm_btemp))  ALLOCATE(mm_btemp(mm_nle))
      IF (.NOT.ALLOCATED(mm_dzlev))  ALLOCATE(mm_dzlev(mm_nla))
      IF (.NOT.ALLOCATED(mm_dzlay))  ALLOCATE(mm_dzlay(mm_nla))
      IF (.NOT.ALLOCATED(mm_rhoair)) ALLOCATE(mm_rhoair(mm_nla))
      
      ! Saves reversed input vectors:
      mm_plev = plev(mm_nle:1:-1)
      mm_zlev = zlev(mm_nle:1:-1)
      mm_play = play(mm_nla:1:-1)
      mm_zlay = zlay(mm_nla:1:-1)
      mm_temp = temp(mm_nla:1:-1)      
      
      ! Computes temperature vertical profil at interfaces:
      mm_btemp(2:mm_nla)   = (mm_temp(1:mm_nla-1) + mm_temp(2:mm_nla)) / 2._mm_wp
      mm_btemp(1)          = mm_temp(1)
      mm_btemp(mm_nle)     = mm_temp(mm_nla) + (mm_temp(mm_nla) - mm_temp(mm_nla-1)) / 2._mm_wp
      
      ! Computes atmospheric levels thickness:
      mm_dzlev(1:mm_nla)   = mm_zlev(1:mm_nle-1)-mm_zlev(2:mm_nle)
      
      ! Computes atmospheric layers thickness :
      mm_dzlay(1:mm_nla-1) = mm_zlay(1:mm_nla-1)-mm_zlay(2:mm_nla)
      mm_dzlay(mm_nla)     = mm_dzlay(mm_nla-1)
      
      ! Hydrostatic equilibrium:
      mm_rhoair(1:mm_nla) = (mm_plev(2:mm_nle)-mm_plev(1:mm_nla)) / (mm_effg(mm_zlay)*mm_dzlev)
      
      ! Write out profiles for debug and log:
      IF (mm_log.AND.mm_debug) THEN
        WRITE(*,'(a)') '# TEMP             PLAY             ZLAY             DZLAY            RHOAIR'
        DO i=1,mm_nla
          WRITE(*,'(5(ES15.7,2X))') mm_temp(i),mm_play(i),mm_zlay(i),mm_dzlay(i), mm_rhoair(i)
        ENDDO
        WRITE(*,'(a)') '# TEMP             PLEV             ZLEV             DZLEV'
        DO i=1,mm_nle
          IF (i /= mm_nle) THEN
            WRITE(*,'(4(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i),mm_dzlev(i)
          ELSE
            WRITE(*,'(3(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i)
          ENDIF
        ENDDO
      ENDIF

      ! End of initialization:
      mm_ini_col = .true.

      RETURN
    END FUNCTION mm_column_init


    FUNCTION mm_aerosols_init(m0aer_s,m3aer_s,m0aer_f,m3aer_f) RESULT(err)
      !! Initialize aerosol 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.
      !!
      !! @warning
      !! The method should be called after mm_global_init and mm_column_init. Moreover, it should be called
      !! whenever the vertical structure of the atmosphere changes.
      !! All the input vectors should be defined from __GROUND__ to __TOP__ of the atmosphere,
      !! otherwise nasty things will occur in computations.
      !!
      REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_s ! 0th order moment of the spherical mode (m-2).
      REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_s ! 3rd order moment of the spherical mode (m3.m-2).
      REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_f ! 0th order moment of the fractal mode (m-2).
      REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_f ! 3rd order moment of the fractal mode (m3.m-2).
      TYPE(error) :: err                                    ! Error status of the function.
      
      err = noerror
      
      ! Global initialization must be done before:
      IF (.NOT.mm_ini) THEN
        err = error("mm_aerosols_init: Global initialization not done yet",-1) ; RETURN
      ENDIF
      
      ! Column initialization must be done before:
      IF (.NOT.mm_ini_col) THEN
        err = error("mm_aerosols_init: Column initialization not done yet",-1) ; RETURN
      ENDIF

      ! Sanity check:
      IF (SIZE(m0aer_s) /= mm_nla) THEN
        err = error("mm_aerosols_init: Invalid size for input arrays",-1) ; RETURN
      ENDIF
  
      ! Allocate variable if required:
      IF (.NOT.ALLOCATED(mm_m0aer_s)) ALLOCATE(mm_m0aer_s(mm_nla))
      IF (.NOT.ALLOCATED(mm_m3aer_s)) ALLOCATE(mm_m3aer_s(mm_nla))
      IF (.NOT.ALLOCATED(mm_m0aer_f)) ALLOCATE(mm_m0aer_f(mm_nla))
      IF (.NOT.ALLOCATED(mm_m3aer_f)) ALLOCATE(mm_m3aer_f(mm_nla))
      IF (.NOT.ALLOCATED(mm_rcs))     ALLOCATE(mm_rcs(mm_nla))
      IF (.NOT.ALLOCATED(mm_rcf))     ALLOCATE(mm_rcf(mm_nla))
      
      ! Allocate memory for diagnostics if required:
      IF (.NOT.ALLOCATED(mm_m0as_vsed)) THEN
        ALLOCATE(mm_m0as_vsed(mm_nla))
        mm_m0as_vsed(:) = 0._mm_wp
      ENDIF
      IF (.NOT.ALLOCATED(mm_m3as_vsed)) THEN
        ALLOCATE(mm_m3as_vsed(mm_nla))
        mm_m3as_vsed(:) = 0._mm_wp
      ENDIF
      IF (.NOT.ALLOCATED(mm_m0af_vsed)) THEN
        ALLOCATE(mm_m0af_vsed(mm_nla))
        mm_m0af_vsed(:) = 0._mm_wp
      ENDIF
      IF (.NOT.ALLOCATED(mm_m3af_vsed)) THEN
        ALLOCATE(mm_m3af_vsed(mm_nla))
        mm_m3af_vsed(:) = 0._mm_wp
      ENDIF
      IF (.NOT.ALLOCATED(mm_aer_s_flux)) THEN
        ALLOCATE(mm_aer_s_flux(mm_nla))
        mm_aer_s_flux(:) = 0._mm_wp
      ENDIF
      IF (.NOT.ALLOCATED(mm_aer_f_flux)) THEN
        ALLOCATE(mm_aer_f_flux(mm_nla))
        mm_aer_f_flux(:) = 0._mm_wp
      ENDIF
      
      ! Initialization of aerosol tracers:
      ! @note: mm_dzlev is already from top to ground
      mm_m0aer_s = m0aer_s(mm_nla:1:-1)/mm_dzlev(:)
      mm_m3aer_s = m3aer_s(mm_nla:1:-1)/mm_dzlev(:)
      mm_m0aer_f = m0aer_f(mm_nla:1:-1)/mm_dzlev(:)
      mm_m3aer_f = m3aer_f(mm_nla:1:-1)/mm_dzlev(:)
  
      ! Setup threshold (useful for preventing bugs):
      call mm_set_moments_thresholds()
  
      ! Initialization of spherical aerosol characteristic 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)
      ELSEWHERE
        mm_rcs = 0._mm_wp
      ENDWHERE

      ! Initialization of fractal aerosol characteristic radii:
      WHERE(mm_m3aer_f > 0._mm_wp .AND. mm_m0aer_f > 0._mm_wp)
        mm_rcf = mm_get_rcf(mm_m0aer_f,mm_m3aer_f)
      ELSEWHERE
        mm_rcf = 0._mm_wp
      ENDWHERE

      ! End of initialization:
      mm_ini_aer = .true.

    END FUNCTION mm_aerosols_init


    !============================================================================
    !                      INTER-MOMENT RELATION METHODS
    !============================================================================
    
    PURE FUNCTION mm_alpha_s(k) RESULT (res)
      !! Inter-moment relation for spherical aerosols size distribution law.
      !! Mk / M0 = rc^k . alpha(k)
      !!
      REAL(kind=mm_wp), INTENT(in) :: k ! k order of the moment.
      REAL(kind=mm_wp) :: sigma         ! Standard deviation.
      REAL(kind=mm_wp) :: res           ! Alpha value.

      ! Titan's case
      !~~~~~~~~~~~~~
      ! res = SUM(dexp(mm_asp%a*k**2 + mm_asp%b*k + mm_asp%c))

      ! Pluto's case
      !~~~~~~~~~~~~~
      sigma = 0.2_mm_wp
      res = exp(k**2 * sigma**2 / 2._mm_wp)

      RETURN
    END FUNCTION mm_alpha_s


    PURE FUNCTION mm_alpha_f(k) RESULT (res)
      !! Inter-moment relation for fractal aerosols size distribution law.
      !! Mk / M0 = rc^k . alpha(k)
      !!
      REAL(kind=mm_wp), INTENT(in) :: k ! k order of the moment.
      REAL(kind=mm_wp) :: sigma         ! Standard deviation.
      REAL(kind=mm_wp) :: res           ! Alpha value.

      ! Titan's case
      !~~~~~~~~~~~~~
      ! res = SUM(dexp(mm_afp%a*k**2 + mm_afp%b*k + mm_afp%c))
      
      ! Pluto's case
      !~~~~~~~~~~~~~
      sigma = 0.35_mm_wp
      res = exp(k**2 * sigma**2 / 2._mm_wp)
      
      RETURN
    END FUNCTION mm_alpha_f


    !============================================================================
    !                          MISCELLANEOUS METHODS
    !============================================================================

    ELEMENTAL FUNCTION mm_effg(z) RESULT(effg)
      !! Compute effective gravitational acceleration.
      !!
      REAL(kind=mm_wp), INTENT(in) :: z ! Altitude (m)
      REAL(kind=mm_wp) :: effg          ! Effective gravitational acceleration (m.s-2)
      effg = mm_g0
      IF (mm_use_effg) THEN
        effg = effg * (mm_rpla/(mm_rpla+z))**2
      ENDIF
      RETURN
    END FUNCTION mm_effg


    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.
      !!
      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_m3aer_s(i) = 0._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_m3aer_f(i) = 0._mm_wp
        ENDIF
      ENDDO
    END SUBROUTINE mm_set_moments_thresholds


    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 spherical aerosol size distribution
      !! law according to its moments and its inter-moments relation.
      !!
      REAL(kind=mm_wp), INTENT(in) :: m0 ! 0th order moment
      REAL(kind=mm_wp), INTENT(in) :: m3 ! 3rd order moment
      REAL(kind=mm_wp) :: res            ! Radius
      res = (m3 / (m0*mm_alpha_s(3._mm_wp)))**(1._mm_wp/3._mm_wp)
    END FUNCTION mm_get_rcs


    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 fractal aerosol size distribution
      !! law according to its moments and its inter-moments relation.
      !!
      REAL(kind=mm_wp), INTENT(in) :: m0 ! 0th order moment
      REAL(kind=mm_wp), INTENT(in) :: m3 ! 3rd order moment
      REAL(kind=mm_wp) :: res            ! Radius
      res = (m3 / (m0*mm_alpha_f(3._mm_wp)))**(1._mm_wp/3._mm_wp)
    END FUNCTION mm_get_rcf


    SUBROUTINE mm_dump_parameters()
      !! Dump 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)')         "---------------------------------------"
      WRITE(*,'(a)')         "Microphysical control flags"
      
      ! Haze production:
      WRITE(*,'(a,L2)')      "mm_w_haze_prod             : ", mm_w_haze_prod
      WRITE(*,'(a,ES14.7)')  "mm_rc_prod (m)             : ", mm_rc_prod
      WRITE(*,'(a,L2)')      "mm_haze_prod_pCH4          : ", mm_haze_prod_pCH4
      IF (.NOT. mm_haze_prod_pCH4) THEN
        WRITE(*,'(a,ES14.7)')  "   --> mm_p_prod (Pa)         : ", mm_p_prod
        WRITE(*,'(a,ES14.7)')  "   --> mm_tx_prod (kg.m-2.s-1): ", mm_tx_prod
        
      ENDIF
      
      ! Haze coagulation:
      WRITE(*,'(a,L2)')      "mm_w_haze_coag             : ", mm_w_haze_coag
      IF (mm_w_haze_coag) THEN
        WRITE(*,'(a,I2.2)')    "   --> mm_coag_interactions   : ", mm_coag_choice
      ENDIF
      
      ! Haze sedimentation:
      WRITE(*,'(a,L2)')      "mm_w_haze_sed              : ", mm_w_haze_sed
      IF (mm_w_haze_sed) THEN
        WRITE(*,'(a,L2)')      "   --> mm_wsed_m0             : ", mm_wsed_m0
        WRITE(*,'(a,L2)')      "   --> mm_wsed_m3             : ", mm_wsed_m3
      ENDIF
      WRITE(*,'(a)')         "---------------------------------------"
      
      ! Haze threshold:
      WRITE(*,'(a)')         "Spherical aerosol thresholds"
      WRITE(*,'(a,ES14.7)')  "  mm_m0as_min (m-3)        : ", mm_m0as_min
      WRITE(*,'(a,ES14.7)')  "  mm_rcs_min (m)           : ", mm_rcs_min
      WRITE(*,'(a)')         "Fractal aerosol thresholds"
      WRITE(*,'(a,ES14.7)')  "  mm_m0af_min (m-3)        : ", mm_m0af_min
      WRITE(*,'(a,ES14.7)')  "  mm_rcf_min (m)           : ", mm_rcf_min
      WRITE(*,'(a)')         "---------------------------------------"
      
      ! Free parameters:
      WRITE(*,'(a)')         "Free parameters"
      WRITE(*,'(a,ES14.7)')  "mm_rm (m)                  : ", mm_rm
      WRITE(*,'(a,ES14.7)')  "mm_df (-)                  : ", mm_df
      WRITE(*,'(a,ES14.7)')  "mm_rhoaer (kg.m-3)         : ", mm_rhoaer      
      WRITE(*,'(a,ES14.7)')  "mm_rplanete (m)            : ", mm_rpla
      WRITE(*,'(a,ES14.7)')  "mm_g0 (m.s-2)              : ", mm_g0
      WRITE(*,'(a,ES14.7)')  "mm_air_rad (m)             : ", mm_air_rad
      WRITE(*,'(a,ES14.7)')  "mm_air_mmol (kg.mol-1)     : ", mm_air_mmol
      WRITE(*,'(a,ES14.7)')  "mm_air_mmas (kg)           : ", mm_air_mmas
      WRITE(*,'(a,ES14.7)')  "mm_dt (s)                  : ", mm_dt
      IF (mm_nla > -1) THEN
        WRITE(*,'(a,I3.3)')    "mm_nla                   : ", mm_nla
      ELSE
        WRITE(*,'(a)')         "mm_nla                   : not initialized yet"
      ENDIF
      WRITE(*,'(a)')         "======================================="
    END SUBROUTINE mm_dump_parameters


    ! =========================================================================
    ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    !                CONFIGURATION PARSER checking methods
    ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ! =========================================================================
  
    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
      !! __var__ on error if given.
      !!
      TYPE(error), INTENT(in)                :: err  ! Error object from value getter.
      REAL(kind=mm_wp), INTENT(inout)        :: var  ! Input/output option value.
      REAL(kind=mm_wp), INTENT(in), OPTIONAL :: def  ! Default value to set.
      LOGICAL, INTENT(in), OPTIONAL          :: wlog ! .true. to print warning/error message.
      TYPE(error) :: ret                             ! Input error.
      CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
      LOGICAL                     :: zlog
      ret = err
      zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
      IF (err == 0) RETURN
      IF (PRESENT(def)) THEN
        var = def
        IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
        ret = noerror
      ELSE
        IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
      ENDIF
    END FUNCTION check_r1
  
    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
      !! __var__ on error if given.
      !!
      TYPE(error), INTENT(in)       :: err  ! Error object from value getter.
      LOGICAL, INTENT(inout)        :: var  ! Input/output option value.
      LOGICAL, INTENT(in), OPTIONAL :: def  ! Default value to set.
      LOGICAL, INTENT(in), OPTIONAL :: wlog ! .true. to print warning/error message.
      TYPE(error) :: ret                    ! Input error.
      CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
      LOGICAL                     :: zlog
      ret = err
      zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
      IF (err == 0) RETURN
      IF (PRESENT(def)) THEN
        var = def
        IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
        ret = noerror
      ELSE
        IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
      ENDIF
    END FUNCTION check_l1
  
    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
      !! __var__ on error if given.
      !!
      TYPE(error), INTENT(in)       :: err  ! Error object from value getter.
      INTEGER, INTENT(inout)        :: var  ! Input/output option value.
      INTEGER, INTENT(in), OPTIONAL :: def  ! Default value to set.
      LOGICAL, INTENT(in), OPTIONAL :: wlog ! .true. to print warning/error message.
      TYPE(error) :: ret                    ! Input error.
      CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
      LOGICAL                     :: zlog
      ret = err
      zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
      IF (err == 0) RETURN
      IF (PRESENT(def)) THEN
        var = def
        IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
        ret = noerror
      ELSE
        IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
      ENDIF
    END FUNCTION check_i1
  
    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
      !! __var__ on error if given.
      !!
      TYPE(error), INTENT(in)                      :: err  ! Error object from value getter.
      CHARACTER(len=*), INTENT(inout)              :: var  ! Input/output option value.
      CHARACTER(len=*), INTENT(in), OPTIONAL       :: def  ! Default value to set.
      LOGICAL, INTENT(in), OPTIONAL                :: wlog ! .true. to print warning/error message.
      TYPE(error) :: ret                                   ! Input error.
      CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
      LOGICAL                     :: zlog
      ret = err
      zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
      IF (err == 0) RETURN
      IF (PRESENT(def)) THEN
        var = TRIM(def)
        IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,var
        ret = noerror
      ELSE
        IF (zlog) WRITE(*,'(a)') error_to_string(err,'')
      ENDIF
      RETURN
    END FUNCTION check_s1

END MODULE MP2M_GLOBALS