diff --git a/aux/ini_phyex.f90 b/aux/ini_phyex.f90
deleted file mode 100644
index 95c371288483251224cb7bc880762ee57bac965f..0000000000000000000000000000000000000000
--- a/aux/ini_phyex.f90
+++ /dev/null
@@ -1,296 +0,0 @@
-SUBROUTINE INI_PHYEX(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, KFROM, KTO, &
-                    &PTSTEP, PDZMIN, &
-                    &CMICRO, CSCONV, CTURB, &
-                    &LDCHANGEMODEL, LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT, LDINIT, &
-                    &PHYEX_IN, PHYEX_OUT)
-!
-!IMPORTANT NOTE ON HOW TO DEAL WITH KEYS USED IN THE PHYSICS *AND* IN THE HOST MODEL
-!
-!When this situation occurs:
-!For example, if a physics option (eg logical switch) has an impact on the model setup, this
-!option must be accessible in the host model *and* in the physics.
-!
-!The prefered solution (when it is possible):
-!The option is implemented in a module and in a namelist of PHYEX. And, after the initialisation
-!is done, the option is available in the PHYEX_OUT structure (in addition to the PHYEX module) for use
-!by the host model.
-!
-!When the physics option is needed in the host model before the PHYEX initialisation, or if the physics
-!option must be computed from other options known by the host model, there are several solutions:
-!
-!  If this situation is specific to only one host model:
-!  The option is declared (module) and initialised (namelist) twice: once in the host model, once in PHYEX;
-!  and a consistency check is added after the physics initialisation.
-!  This solution is preferred to overwriting the variable stored in the PHYEX module, because between
-!  initialization and overwriting, the variable may have been used to initialize a parameterization.
-!
-!  If this issue is common to all the host models:
-!  The option is declared and initialised in the host model. Then, the value is given, through INI_PHYEX, to
-!  the parametrisation INIT subroutine (eg TURBN_INIT). In the INIT subroutine the value is assigned
-!  to a module variable specific to PHYEX (in order to be accessible from inside the parametrisations).
-!  In this case the variable *must not* be added in a PHYEX namelist.
-!  An alternative to this solution is not to copy the variable inside a PHYEX module but add it directly
-!  to the parametrisation call arguments.
-!
-USE MODD_PHYEX, ONLY: PHYEX_t
-USE MODD_CST, ONLY: CST, PRINT_CST
-USE MODD_PARAM_ICE_n, ONLY: PARAM_ICE_GOTO_MODEL, PARAM_ICEN_INIT, PARAM_ICEN, &
-                          & CSUBG_AUCV_RC, CSUBG_AUCV_RI
-USE MODD_RAIN_ICE_DESCR_n, ONLY: RAIN_ICE_DESCR_GOTO_MODEL, RAIN_ICE_DESCRN
-USE MODD_RAIN_ICE_PARAM_n, ONLY: RAIN_ICE_PARAM_GOTO_MODEL, RAIN_ICE_PARAMN
-USE MODD_CLOUDPAR_N,     ONLY: CLOUDPAR_GOTO_MODEL, CLOUDPARN
-USE MODD_PARAM_MFSHALL_N,ONLY: PARAM_MFSHALLN, PARAM_MFSHALLN_INIT, PARAM_MFSHALL_GOTO_MODEL
-USE MODD_TURB_N,         ONLY: TURBN, TURBN_INIT, TURB_GOTO_MODEL, LHARAT
-USE MODD_CTURB,          ONLY: CSTURB, CTURB_ASSOCIATE
-USE MODD_NEB_N,          ONLY: NEBN, NEBN_INIT, NEB_GOTO_MODEL, CCONDENS, LSTATNW, LSUBG_COND
-USE MODD_PARAM_LIMA,     ONLY: PARAM_LIMA, PARAM_LIMA_INIT, PARAM_LIMA_ASSOCIATE, LPTSPLIT, LADJ
-USE MODD_PARAM_LIMA_WARM, ONLY: PARAM_LIMA_WARM
-USE MODD_PARAM_LIMA_COLD, ONLY: PARAM_LIMA_COLD
-USE MODD_PARAM_LIMA_MIXED, ONLY: PARAM_LIMA_MIXED
-USE MODD_NSV,            ONLY: TNSV, NSV_ASSOCIATE
-!
-USE MODE_INI_CST, ONLY: INI_CST
-USE MODE_INI_RAIN_ICE, ONLY: INI_RAIN_ICE
-USE MODE_INI_TIWMX, ONLY: INI_TIWMX
-USE MODE_INI_SNOW, ONLY: INI_SNOW
-USE MODE_INI_MFSHALL, ONLY: INI_MFSHALL
-USE MODE_INI_TURB, ONLY: INI_TURB
-USE MODE_LIMA_UPDATE_NSV, ONLY: LIMA_UPDATE_NSV
-USE MODE_INIT_AEROSOL_PROPERTIES, ONLY: INIT_AEROSOL_PROPERTIES
-USE MODE_INI_LIMA, ONLY: INI_LIMA
-!
-USE MODE_MSG, ONLY: PRINT_MSG, NVERB_FATAL
-!
-!!
-!!      *INI_PHYEX* - PHYEX initialisation routine
-!!
-!!    PURPOSE
-!!    -------
-!!      The purpose of this routine is to initialise (default values, namelist, checks, prints)
-!!      the different modules used by the physics
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!     S. Riette
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!    -  Original      Fev 2023
-!!
-!-------------------------------------------------------------------------------
-!
-!**       DECLARATIONS
-!
-IMPLICIT NONE
-
-CHARACTER(LEN=6),  INTENT(IN) :: HPROGRAM     !< Current program
-INTEGER,           INTENT(IN) :: KUNITNML     !< Logical unit to access the namelist
-LOGICAL,           INTENT(IN) :: LDNEEDNAM    !< True to abort if namelist is absent
-INTEGER,           INTENT(IN) :: KLUOUT       !< Logical unit for outputs
-INTEGER,           INTENT(IN) :: KFROM        !< Old model number
-INTEGER,           INTENT(IN) :: KTO          !< New model number
-REAL,              INTENT(IN) :: PTSTEP       !< Timestep
-REAL,              INTENT(IN) :: PDZMIN       !< Minimum thickness
-CHARACTER(LEN=4),  INTENT(IN) :: CMICRO       !< Microphysical scheme to use
-CHARACTER(LEN=4),  INTENT(IN) :: CTURB        !< Turbulence scheme to use
-CHARACTER(LEN=4),  INTENT(IN) :: CSCONV       !< Shallow convection scheme to use
-LOGICAL, OPTIONAL, INTENT(IN) :: LDCHANGEMODEL!< Must we change the active model
-LOGICAL, OPTIONAL, INTENT(IN) :: LDDEFAULTVAL !< Must we initialize variables with default values (defaults to .TRUE.)
-LOGICAL, OPTIONAL, INTENT(IN) :: LDREADNAM    !< Must we read the namelist (defaults to .TRUE.)
-LOGICAL, OPTIONAL, INTENT(IN) :: LDCHECK      !< Must we perform some checks on values (defaults to .TRUE.)
-INTEGER, OPTIONAL, INTENT(IN) :: KPRINT       !< Print level (defaults to 0): 0 for no print, 1 to safely print namelist,
-                                              !! 2 to print informative messages
-LOGICAL, OPTIONAL, INTENT(IN) :: LDINIT       !< Must we call the init routines
-TYPE(PHYEX_t), OPTIONAL, INTENT(IN)    :: PHYEX_IN    !< Structure for constants (IN)
-TYPE(PHYEX_t), OPTIONAL, INTENT(INOUT) :: PHYEX_OUT   !< Structure for constants (OUT)
-
-!IMPORTANT NOTE on PHYEX_OUT arguments.
-!Logically this argument should be declared with INTENT(OUT) but in this case ifort (at least) breaks the
-!execution when the same structure is given for the PHYEX_IN and the PHYEX_OUT argument.
-!When INITENT(INOUT) is used, execution is OK on ifort.
-
-LOGICAL :: LLINIT, LLCHANGEMODEL, LLCHECK
-INTEGER :: IPRINT
-!
-!**       ARGUMENTS
-!
-LLINIT=.TRUE.
-IF(PRESENT(LDINIT)) LLINIT=LDINIT
-LLCHECK=.TRUE.
-IF(PRESENT(LDCHECK)) LLCHECK=LDCHECK
-LLCHANGEMODEL=.TRUE.
-IF(PRESENT(LDCHANGEMODEL)) LLCHANGEMODEL=LDCHANGEMODEL
-IPRINT=0
-IF(PRESENT(KPRINT)) IPRINT=KPRINT
-!
-!**       CST
-!
-IF(LLINIT) THEN
-  IF(IPRINT==2) WRITE(UNIT=KLUOUT,FMT='('' MODD_CST '')')
-  IF(PRESENT(PHYEX_IN)) CST=PHYEX_IN%CST
-  CALL INI_CST()
-  IF(IPRINT==2) CALL PRINT_CST(KLUOUT)
-  IF(PRESENT(PHYEX_OUT)) PHYEX_OUT%CST=CST
-ENDIF
-!
-!**       MICROPHYSICS SCHEME
-!
-IF(CMICRO=='ICE3' .OR. CMICRO=='ICE4' .OR. CMICRO=='LIMA') THEN
-  !The LIMA scheme makes use of the condensation routine of the ICE3/ICE4 scheme
-  !It is why the ICE3/ICE4 initialisation is needed in for the LIMA scheme
-  IF(IPRINT==2) WRITE(UNIT=KLUOUT,FMT='('' MODD_PARAM_ICEN, MODD_RAIN_ICE_DESCRN, MODD_RAIN_ICE_PARAMN, MODD_CLOUDPARN '')')
-  IF(LLCHANGEMODEL) THEN
-    CALL CLOUDPAR_GOTO_MODEL(KFROM, KTO)
-    CALL PARAM_ICE_GOTO_MODEL(KFROM, KTO)
-    CALL RAIN_ICE_DESCR_GOTO_MODEL(KFROM, KTO)
-    CALL RAIN_ICE_PARAM_GOTO_MODEL(KFROM, KTO)
-    IF(CMICRO=='LIMA') THEN
-      CALL PARAM_LIMA_ASSOCIATE()
-    ENDIF
-  ENDIF
-  IF(PRESENT(PHYEX_IN)) THEN
-    PARAM_ICEN=PHYEX_IN%PARAM_ICEN
-    RAIN_ICE_DESCRN=PHYEX_IN%RAIN_ICE_DESCRN
-    RAIN_ICE_PARAMN=PHYEX_IN%RAIN_ICE_PARAMN
-    CLOUDPARN=PHYEX_IN%CLOUDPARN
-    IF(CMICRO=='LIMA') THEN
-      PARAM_LIMA=PHYEX_IN%PARAM_LIMA
-      PARAM_LIMA_WARM=PHYEX_IN%PARAM_LIMA_WARM
-      PARAM_LIMA_COLD=PHYEX_IN%PARAM_LIMA_COLD
-      PARAM_LIMA_MIXED=PHYEX_IN%PARAM_LIMA_MIXED
-    ENDIF
-  ENDIF
-
-  CALL PARAM_ICEN_INIT(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, &
-                      &LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT)
-  IF(CMICRO=='LIMA') THEN
-    CALL PARAM_LIMA_INIT(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, &
-                        &LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT)
-  ENDIF
-
-  IF(LLINIT) THEN
-    CALL INI_RAIN_ICE(KLUOUT, PTSTEP, PDZMIN, CLOUDPARN%NSPLITR, CMICRO)
-    CALL INI_TIWMX
-  
-    IF(RAIN_ICE_PARAMN%XFRMIN(16) > 0.) THEN
-       CALL INI_SNOW(KLUOUT) ! Recalculate snow parameters :  XCCS = XFRMIN(16),XCXS = XFRMIN(17)
-    ENDIF
-
-    IF(CMICRO=='LIMA') THEN
-      CALL INIT_AEROSOL_PROPERTIES
-      CALL INI_LIMA(PTSTEP, PDZMIN, CLOUDPARN%NSPLITR, CLOUDPARN%NSPLITG)
-    ENDIF
-  ENDIF
-
-  IF(PRESENT(PHYEX_OUT)) THEN
-    PHYEX_OUT%PARAM_ICEN=PARAM_ICEN
-    PHYEX_OUT%RAIN_ICE_DESCRN=RAIN_ICE_DESCRN
-    PHYEX_OUT%RAIN_ICE_PARAMN=RAIN_ICE_PARAMN
-    PHYEX_OUT%CLOUDPARN=CLOUDPARN
-    IF(CMICRO=='LIMA') THEN
-      PHYEX_OUT%PARAM_LIMA=PARAM_LIMA
-      PHYEX_OUT%PARAM_LIMA_WARM=PARAM_LIMA_WARM
-      PHYEX_OUT%PARAM_LIMA_COLD=PARAM_LIMA_COLD
-      PHYEX_OUT%PARAM_LIMA_MIXED=PARAM_LIMA_MIXED
-    ENDIF
-  ENDIF
-ENDIF
-!
-!**       NSV init (must be after LIMA init)
-!
-IF(LLINIT) THEN
-  IF(LLCHANGEMODEL) CALL NSV_ASSOCIATE()
-  IF(PRESENT(PHYEX_IN)) TNSV=PHYEX_IN%TNSV
-  TNSV%NSV=0
-  CALL LIMA_UPDATE_NSV(LDINIT=.TRUE., KMI=KTO, KSV=TNSV%NSV, &
-                       &CDCLOUD=CMICRO, LDUPDATE=.TRUE.)
-  TNSV%NSV=TNSV%NSV_LIMA
-  IF(PRESENT(PHYEX_OUT)) PHYEX_OUT%TNSV=TNSV
-ENDIF
-!
-!**       SHALLOW CONVECTION SCHEME
-!
-IF(CSCONV=='EDKF') THEN
-  IF(IPRINT==2) WRITE(UNIT=KLUOUT,FMT='('' MODD_PARAM_MFSHALL_n '')')
-  IF(LLCHANGEMODEL) CALL PARAM_MFSHALL_GOTO_MODEL(KFROM, KTO)
-  IF(PRESENT(PHYEX_IN)) PARAM_MFSHALLN=PHYEX_IN%PARAM_MFSHALLN
-
-  CALL PARAM_MFSHALLN_INIT(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, &
-                          &LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT)
-  IF(LLINIT) THEN
-    CALL INI_MFSHALL()
-  ENDIF
-
-  IF(PRESENT(PHYEX_OUT)) PHYEX_OUT%PARAM_MFSHALLN=PARAM_MFSHALLN
-ENDIF
-!
-!**       CLOUD SCHEME
-!
-IF(.TRUE.) THEN !Placeholder for configuration without cloud scheme or a different one
-  IF(IPRINT==2) WRITE(UNIT=KLUOUT,FMT='('' MODD_NEB_n '')')
-  IF(LLCHANGEMODEL) CALL NEB_GOTO_MODEL(KFROM, KTO)
-  IF(PRESENT(PHYEX_IN)) NEBN=PHYEX_IN%NEBN
-
-  CALL NEBN_INIT(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, &
-                &LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT)
-  IF(LLINIT) THEN
-    !Nothing to do, everything is read from namelist
-  ENDIF
-
-  IF(PRESENT(PHYEX_OUT)) PHYEX_OUT%NEBN=NEBN
-ENDIF
-!
-!**       TURBULENCE SCHEME
-!
-IF(CTURB=='TKEL') THEN
-  IF(IPRINT==2) WRITE(UNIT=KLUOUT,FMT='('' MODD_TURB_n MODD_CTURB '')')
-  IF(LLCHANGEMODEL) THEN
-    CALL TURB_GOTO_MODEL(KFROM, KTO)
-    CALL CTURB_ASSOCIATE()
-  ENDIF
-  IF(PRESENT(PHYEX_IN)) THEN
-    TURBN=PHYEX_IN%TURBN
-    CSTURB=PHYEX_IN%CSTURB
-  ENDIF
-
-  CALL TURBN_INIT(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, &
-                 &LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT)
-  IF(LLINIT) THEN
-    CALL INI_TURB(HPROGRAM)
-  ENDIF
-
-  IF(PRESENT(PHYEX_OUT)) THEN
-    PHYEX_OUT%TURBN=TURBN
-    PHYEX_OUT%CSTURB=CSTURB
-  ENDIF
-ENDIF
-!
-!**       GLOBAL CONSISTENCY TESTS
-!
-IF(LLCHECK) THEN
-  IF((CMICRO=='ICE3' .OR. CMICRO=='ICE4' .OR. CMICRO=='LIMA') .AND. CTURB=='TKEL') THEN
-    IF ((CSUBG_AUCV_RC == 'ADJU' .OR. CSUBG_AUCV_RI == 'ADJU') .AND. CCONDENS /= 'GAUS') THEN
-      CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'INI_PHYEX', &
-                    &"CSUBG_AUCV_RC and/or CSUBG_AUCV_RI cannot be 'ADJU' if CCONDENS is not 'GAUS'")
-    ENDIF
-    IF (.NOT. LHARAT .AND. LSTATNW) THEN
-      CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'INI_PHYEX', &
-                    &'LSTATNW only tested in combination with HARATU and EDMFm!')
-    ENDIF
-  ENDIF
-  !
-  IF(CMICRO=='LIMA') THEN
-   IF (LSUBG_COND .AND. (.NOT. LPTSPLIT)) THEN
-      CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'INI_PHYEX', &                                                                             
-                    &"YOU MUST USE LPTSPLIT=T WITH CMICRO=LIMA AND LSUBG_COND")
-    END IF
-    IF (LSUBG_COND .AND. (.NOT. LADJ)) THEN
-      CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'INI_PHYEX', &
-                    &"YOU MUST USE LADJ=t WITH CMICRO=LIMA AND LSUBG_COND")
-    END IF
-  ENDIF
-ENDIF
-
-END SUBROUTINE INI_PHYEX
diff --git a/ext/ice4_sedimentation_split_momentum.f90 b/ext/ice4_sedimentation_split_momentum.f90
deleted file mode 100644
index 8b137891791fe96927ad78e64b0aad7bded08bdc..0000000000000000000000000000000000000000
--- a/ext/ice4_sedimentation_split_momentum.f90
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/ext/ini_cturb.f90 b/ext/ini_cturb.f90
deleted file mode 100644
index 8b137891791fe96927ad78e64b0aad7bded08bdc..0000000000000000000000000000000000000000
--- a/ext/ini_cturb.f90
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/ext/modd_param_ice.f90 b/ext/modd_param_ice.f90
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/ext/modd_rain_ice_descr.f90 b/ext/modd_rain_ice_descr.f90
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/ext/modd_rain_ice_param.f90 b/ext/modd_rain_ice_param.f90
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/ext/modn_param_ice.f90 b/ext/modn_param_ice.f90
deleted file mode 100644
index 8b137891791fe96927ad78e64b0aad7bded08bdc..0000000000000000000000000000000000000000
--- a/ext/modn_param_ice.f90
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/ext/modn_turb.f90 b/ext/modn_turb.f90
deleted file mode 100644
index 8b137891791fe96927ad78e64b0aad7bded08bdc..0000000000000000000000000000000000000000
--- a/ext/modn_turb.f90
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/ext/modn_turbn.f90 b/ext/modn_turbn.f90
deleted file mode 100644
index 8b137891791fe96927ad78e64b0aad7bded08bdc..0000000000000000000000000000000000000000
--- a/ext/modn_turbn.f90
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/ext/nrcolss.f90 b/ext/nrcolss.f90
deleted file mode 100644
index 4dbeaa1de9d30855774e4761cb6046206225d372..0000000000000000000000000000000000000000
--- a/ext/nrcolss.f90
+++ /dev/null
@@ -1,316 +0,0 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 init 2006/11/23 10:39:56
-!-----------------------------------------------------------------
-!     ###################
-      MODULE MODI_NRCOLSS
-!     ###################
-!
-INTERFACE
-!
-      SUBROUTINE NRCOLSS( KND, PALPHAS, PNUS, PALPHAR, PNUR,                 &
-                         PESR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR,&
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PNRCOLSS, PAG, PBS, PAS                    )
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUS      ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of aggregates collecting rain 
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates (Thompson 2008)
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PNRCOLSS! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-      END SUBROUTINE NRCOLSS
-!
-END INTERFACE
-!
-      END MODULE MODI_NRCOLSS
-!     ########################################################################
-      SUBROUTINE NRCOLSS( KND, PALPHAS, PNUS, PALPHAR, PNUR,                 &
-                         PESR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR,&
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PNRCOLSS, PAG, PBS, PAS                    )
-!     ########################################################################
-!
-!
-!
-!!****  * -  Build up a look-up table containing the scaled fall speed
-!!           difference between size distributed particles of aggregates and Z
-!!
-!!
-!!    PURPOSE
-!!    -------
-!!      The purpose of this routine is to integrate numerically the scaled fall
-!!      speed difference between aggregates and rain  for use in collection
-!!      kernels FOR CONCENTRATIONS. A first integral of the form
-!!
-!!       infty  Dz_max
-!!           / /
-!!           |{|                                                 }
-!!           |{| E_xz (Dx+Dz)^2 |cxDx^dx-czDz^dz| n(Dz) dDz} n(Dx) dDx
-!!           |{|                                                 }
-!!           / /
-!!          0 Dz_min
-!!
-!!      is evaluated and normalised by a second integral of the form
-!!
-!!              infty
-!!             / /
-!!             |{|                          }
-!!             |{| (Dx+Dz)^2 n(Dz) dDz} n(Dx) dDx
-!!             |{|                          }
-!!             / /
-!!              0
-!!
-!!      The result is stored in a two-dimensional array.
-!! 
-!!**  METHOD
-!!    ------
-!!      The free parameters of the size distribution function of aggregates and Z
-!!      (slope parameter LAMBDA) are discretized with a geometrical rate in a
-!!      specific range
-!!            LAMBDA = exp( (Log(LAMBDA_max) - Log(LAMBDA_min))/N_interval )
-!!      The two above integrals are performed using the trapezoidal scheme.
-!!
-!!    EXTERNAL
-!!    --------
-!!      MODI_GENERAL_GAMMA: Generalized gamma distribution law 
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      MODD_CST           : XPI,XRHOLW
-!!      MODD_RAIN_ICE_DESCR: XAS,XAS,XBS
-!!
-!!    REFERENCE
-!!    ---------
-!!      B.S. Ferrier , 1994 : A Double-Moment Multiple-Phase Four-Class
-!!                            Bulk Ice Scheme,JAS,51,249-280.
-!!
-!!    AUTHOR
-!!    ------
-!!      J.-P. Pinty     * Laboratoire d'Aerologie *
-!!
-!!    MODIFICATIONS 
-!!    -------------
-!!      Original    8/11/95
-!!      M. Taufour    03/2022  Adapted from rrcolss for concentration
-!!      J. Wurtz      03/2022  New snow characteristics
-!!
-!-------------------------------------------------------------------------------
-!
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-!
-USE MODI_GENERAL_GAMMA
-!
-USE MODD_CST
-USE MODD_RAIN_ICE_DESCR_n
-!
-IMPLICIT NONE
-!
-!
-!*       0.1   Declarations of dummy arguments 
-!              ------------------------------- 
-!
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUS      ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of aggregates collecting rain 
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates (Thompson 2008)
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PNRCOLSS! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-!
-INTEGER :: JLBDAS  ! Slope index of the size distribution of aggregates
-INTEGER :: JLBDAR  ! Slope index of the size distribution of rain 
-INTEGER :: JDS     ! Diameter index of a particle of aggregates
-INTEGER :: JDR     ! Diameter index of a particle of rain 
-!
-INTEGER :: INR     ! Number of diameter step for the partial integration
-!
-!
-REAL :: ZLBDAS  ! Current slope parameter LAMBDA of aggregates
-REAL :: ZLBDAR  ! Current slope parameter LAMBDA of rain 
-REAL :: ZDLBDAS ! Growth rate of the slope parameter LAMBDA of aggregates
-REAL :: ZDLBDAR ! Growth rate of the slope parameter LAMBDA of rain 
-REAL :: ZDDS    ! Integration step of the diameter of aggregates
-REAL :: ZDDSCALR! Integration step of the diameter of rain  (scaling integral)
-REAL :: ZDDCOLLR! Integration step of the diameter of rain  (fallspe integral)
-REAL :: ZDS     ! Current diameter of the particle aggregates
-REAL :: ZDR     ! Current diameter of the rain 
-REAL :: ZDRMAX  ! Maximal diameter of the raindrops where the integration ends 
-REAL :: ZCOLLR  ! Single integral of the mass weighted fall speed difference 
-                ! over the spectrum of rain 
-REAL :: ZCOLLDRMAX ! Maximum ending point for the partial integral
-REAL :: ZCOLLSR ! Double integral of the mass weighted fall speed difference
-                ! over the spectra of aggregates and rain 
-REAL :: ZSCALR  ! Single integral of the scaling factor over 
-                ! the spectrum of rain 
-REAL :: ZSCALSR ! Double integral of the scaling factor over
-                ! the spectra of aggregates and rain 
-REAL :: ZFUNC   ! Ancillary function
-REAL :: ZCST1
-!
-!
-!-------------------------------------------------------------------------------
-!
-!
-!*       1       COMPUTE THE SCALED VELOCITY DIFFERENCE IN THE MASS
-!*                               COLLECTION KERNEL,
-!                -------------------------------------------------
-!
-!
-!
-!*       1.0     Initialization
-!
-PNRCOLSS(:,:) = 0.0
-ZCST1 = (3.0/XPI)/XRHOLW
-!
-!*       1.1     Compute the growth rate of the slope factors LAMBDA
-!
-ZDLBDAS = EXP( LOG(PLBDASMAX/PLBDASMIN)/REAL(SIZE(PNRCOLSS(:,:),1)-1) )
-ZDLBDAR = EXP( LOG(PLBDARMAX/PLBDARMIN)/REAL(SIZE(PNRCOLSS(:,:),2)-1) )
-!
-!*       1.2     Scan the slope factors LAMBDAX and LAMBDAZ
-!
-DO JLBDAS = 1,SIZE(PNRCOLSS(:,:),1)
-  ZLBDAS = PLBDASMIN * ZDLBDAS ** (JLBDAS-1) 
-!
-!*       1.3     Compute the diameter steps
-!
-  ZDDS   = PDINFTY / (REAL(KND) * ZLBDAS)
-  DO JLBDAR = 1,SIZE(PNRCOLSS(:,:),2)
-    ZLBDAR = PLBDARMIN * ZDLBDAR ** (JLBDAR-1)
-!
-!*       1.4     Initialize the collection integrals
-!
-    ZSCALSR = 0.0
-    ZCOLLSR = 0.0
-!
-!*       1.5     Compute the diameter steps
-!
-    ZDDSCALR = PDINFTY / (REAL(KND) * ZLBDAR)
-!
-!*       1.6     Scan over the diameters DS and DR
-!
-    DO JDS = 1,KND-1
-      ZDS = ZDDS * REAL(JDS)
-      ZSCALR = 0.0
-      ZCOLLR = 0.0
-      DO JDR = 1,KND-1
-        ZDR = ZDDSCALR * REAL(JDR)
-!
-!*       1.7     Compute the normalization factor by integration over the
-!                dimensional spectrum of rain   
-!
-        ZSCALR = ZSCALR + (ZDS+ZDR)**2 * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
-      END DO
-!
-!*       1.8     Compute the scaled fall speed difference by partial
-!                integration over the dimensional spectrum of rain   
-!
-      ZFUNC = PAG - PAS*ZDS**(PBS-3.0) ! approximate limit is Ds=240 microns
-      IF( ZFUNC>0.0 ) THEN
-        ZDRMAX = ZDS*( ZCST1*ZFUNC )**0.3333333
-        ELSE
-        ZDRMAX = PDINFTY / ZLBDAR
-      END IF
-      IF( ZDS>1.0E-4 ) THEN            ! allow computation if Ds>100 microns
-            ! corresponding to a maximal density of the aggregates of XRHOLW
-        IF( ZDRMAX >= 0.5*ZDDSCALR ) THEN
-          INR = CEILING( ZDRMAX/ZDDSCALR )
-          ZDDCOLLR = ZDRMAX / REAL(INR)
-          IF (INR>=KND ) THEN
-            INR = KND
-            ZDDCOLLR = ZDDSCALR
-          END IF
-          DO JDR = 1,INR-1
-            ZDR = ZDDCOLLR * REAL(JDR)
-            ZCOLLR = ZCOLLR + (ZDS+ZDR)**2                                     &
-                       * PESR * ABS(PFALLS*ZDS**PEXFALLS * EXP(-(PFALLEXPS*ZDS)**PALPHAS)-PFALLR*ZDR**PEXFALLR) &
-                                      * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
-          END DO
-          ZCOLLDRMAX = (ZDS+ZDRMAX)**2                                         &
-                    * PESR * ABS(PFALLS*ZDS**PEXFALLS* EXP(-(PFALLEXPS*ZDS)**PALPHAS)-PFALLR*ZDRMAX**PEXFALLR) &
-                                   * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDRMAX)
-          ZCOLLR = (ZCOLLR + 0.5*ZCOLLDRMAX)*(ZDDCOLLR/ZDDSCALR)
-!
-!*       1.9     Compute the normalization factor by integration over the
-!                dimensional spectrum of aggregates
-!
-          ZFUNC   = GENERAL_GAMMA(PALPHAS,PNUS,ZLBDAS,ZDS)
-          ZSCALSR = ZSCALSR + ZSCALR * ZFUNC
-!
-!*       1.10    Compute the scaled fall speed difference by integration over
-!                the dimensional spectrum of aggregates
-!
-          ZCOLLSR = ZCOLLSR + ZCOLLR * ZFUNC
-        END IF
-!
-! Otherwise ZDRMAX = 0.0 so the density of the graupel cannot be reached 
-!                    and so PRRCOLSS(JLBDAS,JLBDAR) = 0.0        !
-!
-      END IF
-    END DO
-!
-!*       1.11    Scale the fall speed difference
-!
-    IF( ZSCALSR>0.0 ) PNRCOLSS(JLBDAS,JLBDAR) = ZCOLLSR / ZSCALSR
-  END DO
-END DO
-!
-END SUBROUTINE NRCOLSS
diff --git a/ext/nscolrg.f90 b/ext/nscolrg.f90
deleted file mode 100644
index 08de78478dc94cb386c188ec3b0d887960c12f43..0000000000000000000000000000000000000000
--- a/ext/nscolrg.f90
+++ /dev/null
@@ -1,317 +0,0 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 init 2006/11/23 10:43:02
-!-----------------------------------------------------------------
-!     ###################
-      MODULE MODI_NSCOLRG
-!     ###################
-!
-INTERFACE
-!
-      SUBROUTINE NSCOLRG( KND, PALPHAS, PZNUS, PALPHAR, PNUR,                &
-                         PESR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR,&
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PNSCOLRG,PAG, PBS, PAS                     )
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PZNUS     ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of the aggregates collecting rain 
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of the aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of the aggregates
-REAL, INTENT(IN) :: PFALLEXPS  ! Fall speed exponential constant of the aggregates
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PNSCOLRG! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-      END SUBROUTINE NSCOLRG
-!
-END INTERFACE
-!
-      END MODULE MODI_NSCOLRG
-!     ########################################################################
-      SUBROUTINE NSCOLRG( KND, PALPHAS, PZNUS, PALPHAR, PNUR,                &
-                         PESR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR,&
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PNSCOLRG,PAG, PBS, PAS                     )
-!     ########################################################################
-!
-!
-!
-!!****  * -  Build up a look-up table containing the scaled fall speed
-!!           difference between size distributed particles of the aggregates and Z
-!!
-!!
-!!    PURPOSE
-!!    -------
-!!      The purpose of this routine is to integrate numerically the scaled fall
-!!      speed difference between aggregates and rain  for use in collection
-!!      kernels. A first integral of the form
-!!
-!!       infty  Dz_max
-!!           / /
-!!           |{|                                                 }
-!!           |{| E_xz (Dx+Dz)^2 |cxDx^dx-czDz^dz| n(Dz) dDz} n(Dx) dDx
-!!           |{|                                                 }
-!!           / /
-!!          0 Dz_min
-!!
-!!      is evaluated and normalised by a second integral of the form
-!!
-!!              infty
-!!             / /
-!!             |{|                          }
-!!             |{| (Dx+Dz)^2 n(Dz) dDz} n(Dx) dDx
-!!             |{|                          }
-!!             / /
-!!              0
-!!
-!!      The result is stored in a two-dimensional array.
-!! 
-!!**  METHOD
-!!    ------
-!!      The free parameters of the size distribution function of the aggregates
-!!      and Z (slope parameter LAMBDA) are discretized with a geometrical rate 
-!!      in a specific range
-!!            LAMBDA = exp( (Log(LAMBDA_max) - Log(LAMBDA_min))/N_interval )
-!!      The two above integrals are performed using the trapezoidal scheme.
-!!
-!!    EXTERNAL
-!!    --------
-!!      MODI_GENERAL_GAMMA: Generalized gamma distribution law 
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      MODD_CST           : XPI,XRHOLW
-!!      MODD_RAIN_ICE_DESCR: XAS,XAS,XBS
-!!
-!!    REFERENCE
-!!    ---------
-!!      B.S. Ferrier , 1994 : A Double-Moment Multiple-Phase Four-Class
-!!                            Bulk Ice Scheme,JAS,51,249-280.
-!!
-!!    AUTHOR
-!!    ------
-!!      J.-P. Pinty     * Laboratoire d'Aerologie *
-!!
-!!    MODIFICATIONS 
-!!    -------------
-!!      Original    8/11/95
-!!      M. Taufour    03/2022  Adapted from rscolrg for concentration
-!!      J. Wurtz      03/2022  New snow characteristics        
-!!
-!-------------------------------------------------------------------------------
-!
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-USE MODI_GENERAL_GAMMA
-!
-USE MODD_CST
-USE MODD_RAIN_ICE_DESCR_n
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of dummy arguments 
-!              ------------------------------- 
-!
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PZNUS     ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of the aggregates collecting rain 
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of the aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of the aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential constant of the aggregates
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PNSCOLRG! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-!
-INTEGER :: JLBDAS  ! Slope index of the size distribution of the aggregates
-INTEGER :: JLBDAR  ! Slope index of the size distribution of rain 
-INTEGER :: JDS     ! Diameter index of a particle of the aggregates
-INTEGER :: JDR     ! Diameter index of a particle of rain 
-!
-INTEGER :: INR     ! Number of diameter step for the partial integration
-!
-REAL :: ZLBDAS  ! Current slope parameter LAMBDA of the aggregates
-REAL :: ZLBDAR  ! Current slope parameter LAMBDA of rain 
-REAL :: ZDLBDAS ! Growth rate of the slope parameter LAMBDA of the aggregates
-REAL :: ZDLBDAR ! Growth rate of the slope parameter LAMBDA of rain 
-REAL :: ZDDS    ! Integration step of the diameter of the aggregates
-REAL :: ZDDSCALR! Integration step of the diameter of rain  (scaling integral)
-REAL :: ZDDCOLLR! Integration step of the diameter of rain  (fallspe integral)
-REAL :: ZDS     ! Current diameter of the particle aggregates
-REAL :: ZDR     ! Current diameter of the raindrops 
-REAL :: ZDRMIN  ! Minimal diameter of the raindrops where the integration starts
-REAL :: ZDRMAX  ! Maximal diameter of the raindrops where the integration ends
-REAL :: ZCOLLR  ! Single integral of the mass weighted fall speed difference 
-                ! over the spectrum of rain 
-REAL :: ZCOLLDRMIN ! Minimum ending point for the partial integral
-REAL :: ZCOLLSR ! Double integral of the mass weighted fall speed difference
-                ! over the spectra of the aggregates and rain 
-REAL :: ZSCALR  ! Single integral of the scaling factor over 
-                ! the spectrum of rain 
-REAL :: ZSCALSR ! Double integral of the scaling factor over
-                ! the spectra of the aggregates and rain 
-REAL :: ZFUNC   ! Ancillary function
-REAL :: ZCST1
-!
-!
-!-------------------------------------------------------------------------------
-!
-!
-!*       1       COMPUTE THE SCALED VELOCITY DIFFERENCE IN THE MASS
-!*                               COLLECTION KERNEL,
-!                -------------------------------------------------
-!
-!
-!*       1.0     Initialization
-!
-PNSCOLRG(:,:) = 0.0
-ZCST1  = (3.0/XPI)/XRHOLW
-!
-!*       1.1     Compute the growth rate of the slope factors LAMBDA
-!
-ZDLBDAR = EXP( LOG(PLBDARMAX/PLBDARMIN)/REAL(SIZE(PNSCOLRG(:,:),1)-1) )
-ZDLBDAS = EXP( LOG(PLBDASMAX/PLBDASMIN)/REAL(SIZE(PNSCOLRG(:,:),2)-1) )
-!
-!*       1.2     Scan the slope factors LAMBDAX and LAMBDAZ
-!
-DO JLBDAR = 1,SIZE(PNSCOLRG(:,:),1)
-  ZLBDAR = PLBDARMIN * ZDLBDAR ** (JLBDAR-1) 
-  ZDRMAX = PDINFTY / ZLBDAR
-!
-!*       1.3     Compute the diameter steps
-!
-  ZDDSCALR = PDINFTY / (REAL(KND) * ZLBDAR)
-  DO JLBDAS = 1,SIZE(PNSCOLRG(:,:),2)
-    ZLBDAS = PLBDASMIN * ZDLBDAS ** (JLBDAS-1)
-!
-!*       1.4     Initialize the collection integrals
-!
-    ZSCALSR = 0.0
-    ZCOLLSR = 0.0
-!
-!*       1.5     Compute the diameter steps
-!
-    ZDDS     = PDINFTY / (REAL(KND) * ZLBDAS)
-!
-!*       1.6     Scan over the diameters DS and DR
-!
-    DO JDS = 1,KND-1
-      ZDS = ZDDS * REAL(JDS)
-      ZSCALR = 0.0
-      ZCOLLR = 0.0
-      DO JDR = 1,KND-1
-        ZDR = ZDDSCALR * REAL(JDR)
-!
-!*       1.7     Compute the normalization factor by integration over the
-!                dimensional spectrum of rain   
-!
-        ZSCALR = ZSCALR + (ZDS+ZDR)**2 * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
-      END DO
-!
-!*       1.8     Compute the scaled fall speed difference by partial
-!                integration over the dimensional spectrum of rain   
-!
-      ZFUNC = PAG - PAS*ZDS**(PBS-3.0) ! approximate limit is Ds=240 microns
-      IF( ZFUNC>0.0 ) THEN
-        ZDRMIN = ZDS*( ZCST1*ZFUNC )**0.3333333
-        ELSE
-        ZDRMIN = 0.0
-      END IF
-      IF( ZDS>1.0E-4 ) THEN            ! allow computation if Ds>100 microns 
-            ! corresponding to a maximal density of the aggregates of XRHOLW
-        IF( (ZDRMAX-ZDRMIN) >= 0.5*ZDDSCALR ) THEN
-          INR = CEILING( (ZDRMAX-ZDRMIN)/ZDDSCALR )
-          ZDDCOLLR = (ZDRMAX-ZDRMIN) / REAL(INR)
-          DO JDR = 1,INR-1
-            ZDR = ZDDCOLLR * REAL(JDR) + ZDRMIN
-            ZCOLLR = ZCOLLR + (ZDS+ZDR)**2                                     &
-                       * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)                &
-                         * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDR**PEXFALLR)
-          END DO
-          IF( ZDRMIN>0.0 ) THEN
-            ZCOLLDRMIN = (ZDS+ZDRMIN)**2                                       &
-                      * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDRMIN)              &
-                      * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDRMIN**PEXFALLR)
-            ELSE
-            ZCOLLDRMIN = 0.0
-          END IF 
-          ZCOLLR = (ZCOLLR + 0.5*ZCOLLDRMIN)*(ZDDCOLLR/ZDDSCALR)
-!
-!*       1.9     Compute the normalization factor by integration over the
-!                dimensional spectrum of the aggregates
-!
-          ZFUNC   = GENERAL_GAMMA(PALPHAS,PZNUS,ZLBDAS,ZDS)  ! MTaufour : !*(ZDS**PEXMASSS)
-          ZSCALSR = ZSCALSR + ZSCALR * ZFUNC
-!
-!*       1.10    Compute the scaled fall speed difference by integration over
-!                the dimensional spectrum of the aggregates
-!
-          ZCOLLSR = ZCOLLSR + ZCOLLR * ZFUNC
-!
-! Otherwise ZDRMIN>ZDRMAX so PRRCOLSS(JLBDAS,JLBDAR) = 0.0        !
-!
-        END IF
-!
-! Otherwise ZDRMAX = 0.0 so the density of the graupel cannot be reached
-!                    and so PRRCOLSS(JLBDAS,JLBDAR) = 0.0        !
-!
-      END IF
-    END DO
-!
-!*       1.10    Scale the fall speed difference
-!
-    IF( ZSCALSR>0.0 ) PNSCOLRG(JLBDAR,JLBDAS) = ZCOLLSR / ZSCALSR
-  END DO
-END DO
-!
-END SUBROUTINE NSCOLRG
diff --git a/ext/rrcolss.f90 b/ext/rrcolss.f90
deleted file mode 100644
index 0dac2fa04b4ba96dba68bb07e5ded522557851ea..0000000000000000000000000000000000000000
--- a/ext/rrcolss.f90
+++ /dev/null
@@ -1,315 +0,0 @@
-!MNH_LIC Copyright 1995-2019 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!     ###################
-      MODULE MODI_RRCOLSS
-!     ###################
-!
-INTERFACE
-!
-      SUBROUTINE RRCOLSS( KND, PALPHAS, PNUS, PALPHAR, PNUR,                 &
-                         PESR, PEXMASSR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PRRCOLSS, PAG, PBS, PAS                    )
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUS      ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of aggregates collecting rain 
-REAL, INTENT(IN) :: PEXMASSR  ! Mass exponent of rain 
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates (Thompson 2008)
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PRRCOLSS! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-      END SUBROUTINE RRCOLSS
-!
-END INTERFACE
-!
-      END MODULE MODI_RRCOLSS
-!     ########################################################################
-      SUBROUTINE RRCOLSS( KND, PALPHAS, PNUS, PALPHAR, PNUR,                 &
-                         PESR, PEXMASSR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PRRCOLSS, PAG, PBS, PAS                    )
-!     ########################################################################
-!
-!
-!
-!!****  * -  Build up a look-up table containing the scaled fall speed
-!!           difference between size distributed particles of aggregates and Z
-!!
-!!
-!!    PURPOSE
-!!    -------
-!!      The purpose of this routine is to integrate numerically the scaled fall
-!!      speed difference between aggregates and rain  for use in collection
-!!      kernels. A first integral of the form
-!!
-!!       infty  Dz_max
-!!           / /
-!!           |{|                                                 }
-!!           |{| E_xz (Dx+Dz)^2 |cxDx^dx-czDz^dz| Dz^bz n(Dz) dDz} n(Dx) dDx
-!!           |{|                                                 }
-!!           / /
-!!          0 Dz_min
-!!
-!!      is evaluated and normalised by a second integral of the form
-!!
-!!              infty
-!!             / /
-!!             |{|                          }
-!!             |{| (Dx+Dz)^2 Dz^bz n(Dz) dDz} n(Dx) dDx
-!!             |{|                          }
-!!             / /
-!!              0
-!!
-!!      The result is stored in a two-dimensional array.
-!! 
-!!**  METHOD
-!!    ------
-!!      The free parameters of the size distribution function of aggregates and Z
-!!      (slope parameter LAMBDA) are discretized with a geometrical rate in a
-!!      specific range
-!!            LAMBDA = exp( (Log(LAMBDA_max) - Log(LAMBDA_min))/N_interval )
-!!      The two above integrals are performed using the trapezoidal scheme.
-!!
-!!    EXTERNAL
-!!    --------
-!!      MODI_GENERAL_GAMMA: Generalized gamma distribution law 
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      MODD_CST           : XPI,XRHOLW
-!!      MODD_RAIN_ICE_DESCR: XAS,XAS,XBS
-!!
-!!    REFERENCE
-!!    ---------
-!!      B.S. Ferrier , 1994 : A Double-Moment Multiple-Phase Four-Class
-!!                            Bulk Ice Scheme,JAS,51,249-280.
-!!
-!!    AUTHOR
-!!    ------
-!!      J.-P. Pinty     * Laboratoire d'Aerologie *
-!!
-!!    MODIFICATIONS 
-!!    -------------
-!!      Original    8/11/95
-!!
-!  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
-!  J. Wurtz       03/2022: new snow characteristics
-!
-!-------------------------------------------------------------------------------
-!
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-!
-USE MODI_GENERAL_GAMMA
-!
-USE MODD_CST
-USE MODD_RAIN_ICE_DESCR_n
-!
-IMPLICIT NONE
-!
-!
-!*       0.1   Declarations of dummy arguments 
-!              ------------------------------- 
-!
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUS      ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of aggregates collecting rain 
-REAL, INTENT(IN) :: PEXMASSR  ! Mass exponent of rain 
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates  (Thompson 2008)
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PRRCOLSS! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-!
-INTEGER :: JLBDAS  ! Slope index of the size distribution of aggregates
-INTEGER :: JLBDAR  ! Slope index of the size distribution of rain 
-INTEGER :: JDS     ! Diameter index of a particle of aggregates
-INTEGER :: JDR     ! Diameter index of a particle of rain 
-!
-INTEGER :: INR     ! Number of diameter step for the partial integration
-!
-!
-REAL :: ZLBDAS  ! Current slope parameter LAMBDA of aggregates
-REAL :: ZLBDAR  ! Current slope parameter LAMBDA of rain 
-REAL :: ZDLBDAS ! Growth rate of the slope parameter LAMBDA of aggregates
-REAL :: ZDLBDAR ! Growth rate of the slope parameter LAMBDA of rain 
-REAL :: ZDDS    ! Integration step of the diameter of aggregates
-REAL :: ZDDSCALR! Integration step of the diameter of rain  (scaling integral)
-REAL :: ZDDCOLLR! Integration step of the diameter of rain  (fallspe integral)
-REAL :: ZDS     ! Current diameter of the particle aggregates
-REAL :: ZDR     ! Current diameter of the rain 
-REAL :: ZDRMAX  ! Maximal diameter of the raindrops where the integration ends 
-REAL :: ZCOLLR  ! Single integral of the mass weighted fall speed difference 
-                ! over the spectrum of rain 
-REAL :: ZCOLLDRMAX ! Maximum ending point for the partial integral
-REAL :: ZCOLLSR ! Double integral of the mass weighted fall speed difference
-                ! over the spectra of aggregates and rain 
-REAL :: ZSCALR  ! Single integral of the scaling factor over 
-                ! the spectrum of rain 
-REAL :: ZSCALSR ! Double integral of the scaling factor over
-                ! the spectra of aggregates and rain 
-REAL :: ZFUNC   ! Ancillary function
-REAL :: ZCST1
-!
-!
-!-------------------------------------------------------------------------------
-!
-!
-!*       1       COMPUTE THE SCALED VELOCITY DIFFERENCE IN THE MASS
-!*                               COLLECTION KERNEL,
-!                -------------------------------------------------
-!
-!
-!
-!*       1.0     Initialization
-!
-PRRCOLSS(:,:) = 0.0
-ZCST1 = (3.0/XPI)/XRHOLW
-!
-!*       1.1     Compute the growth rate of the slope factors LAMBDA
-!
-ZDLBDAS = EXP( LOG(PLBDASMAX/PLBDASMIN)/REAL(SIZE(PRRCOLSS(:,:),1)-1) )
-ZDLBDAR = EXP( LOG(PLBDARMAX/PLBDARMIN)/REAL(SIZE(PRRCOLSS(:,:),2)-1) )
-!
-!*       1.2     Scan the slope factors LAMBDAX and LAMBDAZ
-!
-DO JLBDAS = 1,SIZE(PRRCOLSS(:,:),1)
-  ZLBDAS = PLBDASMIN * ZDLBDAS ** (JLBDAS-1) 
-!
-!*       1.3     Compute the diameter steps
-!
-  ZDDS   = PDINFTY / (REAL(KND) * ZLBDAS)
-  DO JLBDAR = 1,SIZE(PRRCOLSS(:,:),2)
-    ZLBDAR = PLBDARMIN * ZDLBDAR ** (JLBDAR-1)
-!
-!*       1.4     Initialize the collection integrals
-!
-    ZSCALSR = 0.0
-    ZCOLLSR = 0.0
-!
-!*       1.5     Compute the diameter steps
-!
-    ZDDSCALR = PDINFTY / (REAL(KND) * ZLBDAR)
-!
-!*       1.6     Scan over the diameters DS and DR
-!
-    DO JDS = 1,KND-1
-      ZDS = ZDDS * REAL(JDS)
-      ZSCALR = 0.0
-      ZCOLLR = 0.0
-      DO JDR = 1,KND-1
-        ZDR = ZDDSCALR * REAL(JDR)
-!
-!*       1.7     Compute the normalization factor by integration over the
-!                dimensional spectrum of rain   
-!
-        ZSCALR = ZSCALR + (ZDS+ZDR)**2 * ZDR**PEXMASSR                         &
-                                      * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
-      END DO
-!
-!*       1.8     Compute the scaled fall speed difference by partial
-!                integration over the dimensional spectrum of rain   
-!
-      ZFUNC = PAG - PAS*ZDS**(PBS-3.0) ! approximate limit is Ds=240 microns
-      IF( ZFUNC>0.0 ) THEN
-        ZDRMAX = ZDS*( ZCST1*ZFUNC )**0.3333333
-        ELSE
-        ZDRMAX = PDINFTY / ZLBDAR
-      END IF
-      IF( ZDS>1.0E-4 ) THEN            ! allow computation if Ds>100 microns
-            ! corresponding to a maximal density of the aggregates of XRHOLW
-        IF( ZDRMAX >= 0.5*ZDDSCALR ) THEN
-          INR = CEILING( ZDRMAX/ZDDSCALR )
-          ZDDCOLLR = ZDRMAX / REAL(INR)
-          IF (INR>=KND ) THEN
-            INR = KND
-            ZDDCOLLR = ZDDSCALR
-          END IF
-          DO JDR = 1,INR-1
-            ZDR = ZDDCOLLR * REAL(JDR)
-            ZCOLLR = ZCOLLR + (ZDS+ZDR)**2 * ZDR**PEXMASSR                     &
-                       * PESR * ABS(PFALLS*ZDS**PEXFALLS * EXP(-(PFALLEXPS*ZDS)**PALPHAS)-PFALLR*ZDR**PEXFALLR) &
-                                      * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
-          END DO
-          ZCOLLDRMAX = (ZDS+ZDRMAX)**2 * ZDRMAX**PEXMASSR                      &
-                    * PESR * ABS(PFALLS*ZDS**PEXFALLS* EXP(-(PFALLEXPS*ZDS)**PALPHAS)-PFALLR*ZDRMAX**PEXFALLR) &
-                                   * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDRMAX)
-          ZCOLLR = (ZCOLLR + 0.5*ZCOLLDRMAX)*(ZDDCOLLR/ZDDSCALR)
-!
-!*       1.9     Compute the normalization factor by integration over the
-!                dimensional spectrum of aggregates
-!
-          ZFUNC   = GENERAL_GAMMA(PALPHAS,PNUS,ZLBDAS,ZDS)
-          ZSCALSR = ZSCALSR + ZSCALR * ZFUNC
-!
-!*       1.10    Compute the scaled fall speed difference by integration over
-!                the dimensional spectrum of aggregates
-!
-          ZCOLLSR = ZCOLLSR + ZCOLLR * ZFUNC
-        END IF
-!
-! Otherwise ZDRMAX = 0.0 so the density of the graupel cannot be reached 
-!                    and so PRRCOLSS(JLBDAS,JLBDAR) = 0.0        !
-!
-      END IF
-    END DO
-!
-!*       1.11    Scale the fall speed difference
-!
-    IF( ZSCALSR>0.0 ) PRRCOLSS(JLBDAS,JLBDAR) = ZCOLLSR / ZSCALSR
-  END DO
-END DO
-!
-END SUBROUTINE RRCOLSS
diff --git a/ext/rscolrg.f90 b/ext/rscolrg.f90
deleted file mode 100644
index 26969afcd4b3282beafb10e659a0b0d547cfc125..0000000000000000000000000000000000000000
--- a/ext/rscolrg.f90
+++ /dev/null
@@ -1,315 +0,0 @@
-!MNH_LIC Copyright 1995-2019 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-!     ###################
-      MODULE MODI_RSCOLRG
-!     ###################
-!
-INTERFACE
-!
-      SUBROUTINE RSCOLRG( KND, PALPHAS, PZNUS, PALPHAR, PNUR,                &
-                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PRSCOLRG,PAG, PBS, PAS                     )
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PZNUS     ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of the aggregates collecting rain 
-REAL, INTENT(IN) :: PEXMASSS  ! Mass exponent of the aggregates
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of the aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of the aggregates
-REAL, INTENT(IN) :: PFALLEXPS  ! Fall speed exponential constant of the aggregates
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PRSCOLRG! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-      END SUBROUTINE RSCOLRG
-!
-END INTERFACE
-!
-      END MODULE MODI_RSCOLRG
-!     ########################################################################
-      SUBROUTINE RSCOLRG( KND, PALPHAS, PZNUS, PALPHAR, PNUR,                &
-                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
-                         PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
-                         PDINFTY, PRSCOLRG,PAG, PBS, PAS                     )
-!     ########################################################################
-!
-!
-!
-!!****  * -  Build up a look-up table containing the scaled fall speed
-!!           difference between size distributed particles of the aggregates and Z
-!!
-!!
-!!    PURPOSE
-!!    -------
-!!      The purpose of this routine is to integrate numerically the scaled fall
-!!      speed difference between aggregates and rain  for use in collection
-!!      kernels. A first integral of the form
-!!
-!!       infty  Dz_max
-!!           / /
-!!           |{|                                                 }
-!!           |{| E_xz (Dx+Dz)^2 |cxDx^dx-czDz^dz| Dz^bz n(Dz) dDz} n(Dx) dDx
-!!           |{|                                                 }
-!!           / /
-!!          0 Dz_min
-!!
-!!      is evaluated and normalised by a second integral of the form
-!!
-!!              infty
-!!             / /
-!!             |{|                          }
-!!             |{| (Dx+Dz)^2 Dz^bz n(Dz) dDz} n(Dx) dDx
-!!             |{|                          }
-!!             / /
-!!              0
-!!
-!!      The result is stored in a two-dimensional array.
-!! 
-!!**  METHOD
-!!    ------
-!!      The free parameters of the size distribution function of the aggregates
-!!      and Z (slope parameter LAMBDA) are discretized with a geometrical rate 
-!!      in a specific range
-!!            LAMBDA = exp( (Log(LAMBDA_max) - Log(LAMBDA_min))/N_interval )
-!!      The two above integrals are performed using the trapezoidal scheme.
-!!
-!!    EXTERNAL
-!!    --------
-!!      MODI_GENERAL_GAMMA: Generalized gamma distribution law 
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      MODD_CST           : XPI,XRHOLW
-!!      MODD_RAIN_ICE_DESCR: XAS,XAS,XBS
-!!
-!!    REFERENCE
-!!    ---------
-!!      B.S. Ferrier , 1994 : A Double-Moment Multiple-Phase Four-Class
-!!                            Bulk Ice Scheme,JAS,51,249-280.
-!!
-!!    AUTHOR
-!!    ------
-!!      J.-P. Pinty     * Laboratoire d'Aerologie *
-!!
-!!    MODIFICATIONS 
-!!    -------------
-!!      Original    8/11/95
-!!
-!  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
-!  J. Wurtz       03/2022: new snow characteristics
-!
-!-------------------------------------------------------------------------------
-!
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-USE MODI_GENERAL_GAMMA
-!
-USE MODD_CST
-USE MODD_RAIN_ICE_DESCR_n
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of dummy arguments 
-!              ------------------------------- 
-!
-!
-INTEGER, INTENT(IN) :: KND    ! Number of discrete size intervals in DS and DR  
-!
-REAL, INTENT(IN) :: PALPHAS   ! First shape parameter of the aggregates 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PZNUS     ! Second shape parameter of the aggregates
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PALPHAR   ! First shape parameter of the rain  
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PNUR      ! Second shape parameter of the rain 
-                              ! size distribution (generalized gamma law)
-REAL, INTENT(IN) :: PESR      ! Efficiency of the aggregates collecting rain 
-REAL, INTENT(IN) :: PEXMASSS  ! Mass exponent of the aggregates
-REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of the aggregates
-REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of the aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential constant of the aggregates
-REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
-REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
-REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMAX ! Maximun slope of size distribution of rain 
-REAL, INTENT(IN) :: PLBDASMIN ! Minimun slope of size distribution of the aggregates
-REAL, INTENT(IN) :: PLBDARMIN ! Minimun slope of size distribution of rain 
-REAL, INTENT(IN) :: PDINFTY   ! Factor to define the largest diameter up to
-                              ! which the diameter integration is performed
-REAL, INTENT(IN) :: PAG, PBS, PAS
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PRSCOLRG! Scaled fall speed difference in
-                                               ! the mass collection kernel as a
-                                               ! function of LAMBDAX and LAMBDAZ
-!
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-!
-INTEGER :: JLBDAS  ! Slope index of the size distribution of the aggregates
-INTEGER :: JLBDAR  ! Slope index of the size distribution of rain 
-INTEGER :: JDS     ! Diameter index of a particle of the aggregates
-INTEGER :: JDR     ! Diameter index of a particle of rain 
-!
-INTEGER :: INR     ! Number of diameter step for the partial integration
-!
-REAL :: ZLBDAS  ! Current slope parameter LAMBDA of the aggregates
-REAL :: ZLBDAR  ! Current slope parameter LAMBDA of rain 
-REAL :: ZDLBDAS ! Growth rate of the slope parameter LAMBDA of the aggregates
-REAL :: ZDLBDAR ! Growth rate of the slope parameter LAMBDA of rain 
-REAL :: ZDDS    ! Integration step of the diameter of the aggregates
-REAL :: ZDDSCALR! Integration step of the diameter of rain  (scaling integral)
-REAL :: ZDDCOLLR! Integration step of the diameter of rain  (fallspe integral)
-REAL :: ZDS     ! Current diameter of the particle aggregates
-REAL :: ZDR     ! Current diameter of the raindrops 
-REAL :: ZDRMIN  ! Minimal diameter of the raindrops where the integration starts
-REAL :: ZDRMAX  ! Maximal diameter of the raindrops where the integration ends
-REAL :: ZCOLLR  ! Single integral of the mass weighted fall speed difference 
-                ! over the spectrum of rain 
-REAL :: ZCOLLDRMIN ! Minimum ending point for the partial integral
-REAL :: ZCOLLSR ! Double integral of the mass weighted fall speed difference
-                ! over the spectra of the aggregates and rain 
-REAL :: ZSCALR  ! Single integral of the scaling factor over 
-                ! the spectrum of rain 
-REAL :: ZSCALSR ! Double integral of the scaling factor over
-                ! the spectra of the aggregates and rain 
-REAL :: ZFUNC   ! Ancillary function
-REAL :: ZCST1
-!
-!
-!-------------------------------------------------------------------------------
-!
-!
-!*       1       COMPUTE THE SCALED VELOCITY DIFFERENCE IN THE MASS
-!*                               COLLECTION KERNEL,
-!                -------------------------------------------------
-!
-!
-!*       1.0     Initialization
-!
-PRSCOLRG(:,:) = 0.0
-ZCST1  = (3.0/XPI)/XRHOLW
-!
-!*       1.1     Compute the growth rate of the slope factors LAMBDA
-!
-ZDLBDAR = EXP( LOG(PLBDARMAX/PLBDARMIN)/REAL(SIZE(PRSCOLRG(:,:),1)-1) )
-ZDLBDAS = EXP( LOG(PLBDASMAX/PLBDASMIN)/REAL(SIZE(PRSCOLRG(:,:),2)-1) )
-!
-!*       1.2     Scan the slope factors LAMBDAX and LAMBDAZ
-!
-DO JLBDAR = 1,SIZE(PRSCOLRG(:,:),1)
-  ZLBDAR = PLBDARMIN * ZDLBDAR ** (JLBDAR-1) 
-  ZDRMAX = PDINFTY / ZLBDAR
-!
-!*       1.3     Compute the diameter steps
-!
-  ZDDSCALR = PDINFTY / (REAL(KND) * ZLBDAR)
-  DO JLBDAS = 1,SIZE(PRSCOLRG(:,:),2)
-    ZLBDAS = PLBDASMIN * ZDLBDAS ** (JLBDAS-1)
-!
-!*       1.4     Initialize the collection integrals
-!
-    ZSCALSR = 0.0
-    ZCOLLSR = 0.0
-!
-!*       1.5     Compute the diameter steps
-!
-    ZDDS     = PDINFTY / (REAL(KND) * ZLBDAS)
-!
-!*       1.6     Scan over the diameters DS and DR
-!
-    DO JDS = 1,KND-1
-      ZDS = ZDDS * REAL(JDS)
-      ZSCALR = 0.0
-      ZCOLLR = 0.0
-      DO JDR = 1,KND-1
-        ZDR = ZDDSCALR * REAL(JDR)
-!
-!*       1.7     Compute the normalization factor by integration over the
-!                dimensional spectrum of rain   
-!
-        ZSCALR = ZSCALR + (ZDS+ZDR)**2 * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
-      END DO
-!
-!*       1.8     Compute the scaled fall speed difference by partial
-!                integration over the dimensional spectrum of rain   
-!
-      ZFUNC = PAG - PAS*ZDS**(PBS-3.0) ! approximate limit is Ds=240 microns
-      IF( ZFUNC>0.0 ) THEN
-        ZDRMIN = ZDS*( ZCST1*ZFUNC )**0.3333333
-        ELSE
-        ZDRMIN = 0.0
-      END IF
-      IF( ZDS>1.0E-4 ) THEN            ! allow computation if Ds>100 microns 
-            ! corresponding to a maximal density of the aggregates of XRHOLW
-        IF( (ZDRMAX-ZDRMIN) >= 0.5*ZDDSCALR ) THEN
-          INR = CEILING( (ZDRMAX-ZDRMIN)/ZDDSCALR )
-          ZDDCOLLR = (ZDRMAX-ZDRMIN) / REAL(INR)
-          DO JDR = 1,INR-1
-            ZDR = ZDDCOLLR * REAL(JDR) + ZDRMIN
-            ZCOLLR = ZCOLLR + (ZDS+ZDR)**2                                     &
-                       * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)                &
-                         * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDR**PEXFALLR)
-          END DO
-          IF( ZDRMIN>0.0 ) THEN
-            ZCOLLDRMIN = (ZDS+ZDRMIN)**2                                       &
-                      * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDRMIN)              &
-                      * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDRMIN**PEXFALLR)
-            ELSE
-            ZCOLLDRMIN = 0.0
-          END IF 
-          ZCOLLR = (ZCOLLR + 0.5*ZCOLLDRMIN)*(ZDDCOLLR/ZDDSCALR)
-!
-!*       1.9     Compute the normalization factor by integration over the
-!                dimensional spectrum of the aggregates
-!
-          ZFUNC   = (ZDS**PEXMASSS) * GENERAL_GAMMA(PALPHAS,PZNUS,ZLBDAS,ZDS)
-          ZSCALSR = ZSCALSR + ZSCALR * ZFUNC
-!
-!*       1.10    Compute the scaled fall speed difference by integration over
-!                the dimensional spectrum of the aggregates
-!
-          ZCOLLSR = ZCOLLSR + ZCOLLR * ZFUNC
-!
-! Otherwise ZDRMIN>ZDRMAX so PRRCOLSS(JLBDAS,JLBDAR) = 0.0        !
-!
-        END IF
-!
-! Otherwise ZDRMAX = 0.0 so the density of the graupel cannot be reached
-!                    and so PRRCOLSS(JLBDAS,JLBDAR) = 0.0        !
-!
-      END IF
-    END DO
-!
-!*       1.10    Scale the fall speed difference
-!
-    IF( ZSCALSR>0.0 ) PRSCOLRG(JLBDAR,JLBDAS) = ZCOLLSR / ZSCALSR
-  END DO
-END DO
-!
-END SUBROUTINE RSCOLRG