diff --git a/src/MNH/lima_hail.f90 b/src/MNH/lima_hail.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9db77ce1fb2fb8a715a265b79732fd2f32b11c8f
--- /dev/null
+++ b/src/MNH/lima_hail.f90
@@ -0,0 +1,574 @@
+!MNH_LIC Copyright 2018-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_LIMA_HAIL
+! #################################
+!
+INTERFACE
+ SUBROUTINE LIMA_HAIL (PTSTEP, LDCOMPUTE, &
+ PRHODREF, PPRES, PT, PKA, PDV, PCJ, &
+ PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PRHT, &
+ PCCT, PCRT, PCIT, PCST, PCGT, PCHT, &
+ PLBDC, PLBDR, PLBDS, PLBDG, PLBDH, &
+ PLVFACT, PLSFACT, &
+ P_TH_WETH, P_RC_WETH, P_CC_WETH, P_RR_WETH, P_CR_WETH, &
+ P_RI_WETH, P_CI_WETH, P_RS_WETH, P_CS_WETH, P_RG_WETH, P_CG_WETH, P_RH_WETH, &
+ P_RG_COHG, P_CG_COHG, &
+ P_TH_HMLT, P_RR_HMLT, P_CR_HMLT, P_CH_HMLT, &
+ PA_TH, PA_RC, PA_CC, PA_RR, PA_CR, &
+ PA_RI, PA_CI, PA_RS, PA_CS, PA_RG, PA_CG, PA_RH, PA_CH )
+!
+REAL, INTENT(IN) :: PTSTEP
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF !
+REAL, DIMENSION(:), INTENT(IN) :: PPRES !
+REAL, DIMENSION(:), INTENT(IN) :: PT !
+REAL, DIMENSION(:), INTENT(IN) :: PKA !
+REAL, DIMENSION(:), INTENT(IN) :: PDV !
+REAL, DIMENSION(:), INTENT(IN) :: PCJ !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRVT !
+REAL, DIMENSION(:), INTENT(IN) :: PRCT !
+REAL, DIMENSION(:), INTENT(IN) :: PRRT !
+REAL, DIMENSION(:), INTENT(IN) :: PRIT !
+REAL, DIMENSION(:), INTENT(IN) :: PRST !
+REAL, DIMENSION(:), INTENT(IN) :: PRGT !
+REAL, DIMENSION(:), INTENT(IN) :: PRHT !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PCCT !
+REAL, DIMENSION(:), INTENT(IN) :: PCRT !
+REAL, DIMENSION(:), INTENT(IN) :: PCIT !
+REAL, DIMENSION(:), INTENT(IN) :: PCST !
+REAL, DIMENSION(:), INTENT(IN) :: PCGT !
+REAL, DIMENSION(:), INTENT(IN) :: PCHT !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PLBDC !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDR !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDS !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDG !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDH !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PLVFACT !
+REAL, DIMENSION(:), INTENT(IN) :: PLSFACT !
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RC_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CC_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RR_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CR_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RS_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CS_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RG_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CG_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RH_WETH
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RG_COHG
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CG_COHG
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_HMLT
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RR_HMLT
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CR_HMLT
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CH_HMLT
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_TH
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RC
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CC
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RR
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CR
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RI
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CI
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RS
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CS
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RG
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CG
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RH
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CH
+!
+END SUBROUTINE LIMA_HAIL
+END INTERFACE
+END MODULE MODI_LIMA_HAIL
+!
+! #################################################################################
+ SUBROUTINE LIMA_HAIL (PTSTEP, LDCOMPUTE, &
+ PRHODREF, PPRES, PT, PKA, PDV, PCJ, &
+ PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PRHT, &
+ PCCT, PCRT, PCIT, PCST, PCGT, PCHT, &
+ PLBDC, PLBDR, PLBDS, PLBDG, PLBDH, &
+ PLVFACT, PLSFACT, &
+ P_TH_WETH, P_RC_WETH, P_CC_WETH, P_RR_WETH, P_CR_WETH, &
+ P_RI_WETH, P_CI_WETH, P_RS_WETH, P_CS_WETH, P_RG_WETH, P_CG_WETH, P_RH_WETH, &
+ P_RG_COHG, P_CG_COHG, &
+ P_TH_HMLT, P_RR_HMLT, P_CR_HMLT, P_CH_HMLT, &
+ PA_TH, PA_RC, PA_CC, PA_RR, PA_CR, &
+ PA_RI, PA_CI, PA_RS, PA_CS, PA_RG, PA_CG, PA_RH, PA_CH )
+! #################################################################################
+!
+!! PURPOSE
+!! -------
+!! Compute the wet/dry growth of graupel, associated Hallett-Mossop ice production,
+!! and graupel melting rates
+!!
+!! AUTHOR
+!! ------
+!! J.-M. Cohard * Laboratoire d'Aerologie*
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!! S. Berthet * Laboratoire d'Aerologie*
+!! B. Vié * CNRM *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 15/03/2018
+! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+! J. Wurtz 03/2022: new snow characteristics
+!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST, ONLY : XTT, XMD, XMV, XRD, XRV, XLVTT, XLMTT, XESTT, XCL, XCI, XCPV
+USE MODD_PARAM_LIMA, ONLY : XRTMIN, XCTMIN, XCEXVT, LHAIL
+USE MODD_PARAM_LIMA_MIXED, ONLY : NWETLBDAG, XWETINTP1G, XWETINTP2G, &
+ NWETLBDAH, X0DEPH, X1DEPH, XDH, XEX0DEPH, XEX1DEPH, &
+ XFWETH, XWETINTP1H, XWETINTP2H, &
+ NWETLBDAS, XWETINTP1S, XWETINTP2S, &
+ XKER_N_GWETH, XKER_GWETH, XKER_N_SWETH, XKER_SWETH, &
+ XFSWETH, XLBSWETH1, XLBSWETH2, XLBSWETH3, &
+ XFNSWETH, XLBNSWETH1, XLBNSWETH2, XLBNSWETH3, &
+ XFGWETH, XLBGWETH1, XLBGWETH2, XLBGWETH3, &
+ XFNGWETH, XLBNGWETH1, XLBNGWETH2, XLBNGWETH3
+
+USE MODD_PARAM_LIMA_COLD, ONLY : XMNU0, XCXS, XBS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+REAL, INTENT(IN) :: PTSTEP
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF !
+REAL, DIMENSION(:), INTENT(IN) :: PPRES !
+REAL, DIMENSION(:), INTENT(IN) :: PT !
+REAL, DIMENSION(:), INTENT(IN) :: PKA !
+REAL, DIMENSION(:), INTENT(IN) :: PDV !
+REAL, DIMENSION(:), INTENT(IN) :: PCJ !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRVT !
+REAL, DIMENSION(:), INTENT(IN) :: PRCT !
+REAL, DIMENSION(:), INTENT(IN) :: PRRT !
+REAL, DIMENSION(:), INTENT(IN) :: PRIT !
+REAL, DIMENSION(:), INTENT(IN) :: PRST !
+REAL, DIMENSION(:), INTENT(IN) :: PRGT !
+REAL, DIMENSION(:), INTENT(IN) :: PRHT !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PCCT !
+REAL, DIMENSION(:), INTENT(IN) :: PCRT !
+REAL, DIMENSION(:), INTENT(IN) :: PCIT !
+REAL, DIMENSION(:), INTENT(IN) :: PCST !
+REAL, DIMENSION(:), INTENT(IN) :: PCGT !
+REAL, DIMENSION(:), INTENT(IN) :: PCHT !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PLBDC !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDR !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDS !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDG !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDH !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PLVFACT !
+REAL, DIMENSION(:), INTENT(IN) :: PLSFACT !
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RC_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CC_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RR_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CR_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RS_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CS_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RG_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CG_WETH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RH_WETH
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RG_COHG
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CG_COHG
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_HMLT
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RR_HMLT
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CR_HMLT
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CH_HMLT
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_TH
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RC
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CC
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RR
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CR
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RI
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CI
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RS
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CS
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RG
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CG
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_RH
+REAL, DIMENSION(:), INTENT(INOUT) :: PA_CH
+!
+!* 0.2 Declarations of local variables :
+!
+LOGICAL, DIMENSION(SIZE(PRCT)) :: GWET
+INTEGER :: JJ
+!
+REAL, DIMENSION(SIZE(PRCT)) :: Z1, Z2, Z3, Z4
+REAL, DIMENSION(SIZE(PRCT)) :: ZZX, ZZW, ZZW1, ZZW2, ZZW3, ZZW4, ZZW5, ZZW6
+REAL, DIMENSION(SIZE(PRCT)) :: ZZW3N, ZZW4N, ZZW6N
+REAL, DIMENSION(SIZE(PRCT)) :: ZRWETH
+!
+INTEGER, DIMENSION(SIZE(PRCT)) :: IVEC1,IVEC2 ! Vectors of indices
+REAL, DIMENSION(SIZE(PRCT)) :: ZVEC1,ZVEC2, ZVEC3 ! Work vectors
+!
+INTEGER :: NHAIL
+REAL :: ZTHRH, ZTHRC
+!
+!-------------------------------------------------------------------------------
+!
+!
+P_TH_WETH(:) = 0.
+P_RC_WETH(:) = 0.
+P_CC_WETH(:) = 0.
+P_RR_WETH(:) = 0.
+P_CR_WETH(:) = 0.
+P_RI_WETH(:) = 0.
+P_CI_WETH(:) = 0.
+P_RS_WETH(:) = 0.
+P_CS_WETH(:) = 0.
+P_RG_WETH(:) = 0.
+P_CG_WETH(:) = 0.
+P_RH_WETH(:) = 0.
+!
+P_RG_COHG(:) = 0.
+P_CG_COHG(:) = 0.
+!
+P_TH_HMLT(:) = 0.
+P_RR_HMLT(:) = 0.
+P_CR_HMLT(:) = 0.
+P_CH_HMLT(:) = 0.
+!
+ZZW1(:) = 0. ! RCWETH
+ZZW2(:) = 0. ! RIWETH
+ZZW3(:) = 0. ! RSWETH
+ZZW3N(:) = 0.! NSWETH
+ZZW4(:) = 0. ! RGWETH
+ZZW4N(:) = 0.! NGWETH
+ZZW5(:) = 0. ! RCWETH+RRWETH
+ZZW6(:) = 0. ! RSWETH
+!
+ZRWETH(:) = 0.
+!
+!
+!* 1. Hail growth by collection (wet case only ?)
+! -----------------------------------------------
+!
+! 1.a Collection of rc and ri
+! ---------------------------
+!
+WHERE( PRHT(:)>XRTMIN(7) .AND. PCHT(:)>XCTMIN(7) .AND. LDCOMPUTE(:) )
+ ZZW(:) = PCHT(:) * PLBDH(:)**(-XDH-2.0) * PRHODREF(:)**(1-XCEXVT)
+ ZZW1(:) = XFWETH * PRCT(:) * ZZW(:) ! RCWETH
+ ZZW2(:) = XFWETH * PRIT(:) * ZZW(:) ! RIWETH
+END WHERE
+!
+!* 1.b Collection of rs
+! --------------------
+!
+GWET(:) = PRST(:)>XRTMIN(5) .AND. PRHT(:)>XRTMIN(7) .AND. LDCOMPUTE(:) .AND. &
+ PCST(:)>XCTMIN(5) .AND. PCHT(:)>XCTMIN(7)
+!
+WHERE( GWET )
+!
+!* Select the (ZLBDAG,ZLBDAS) couplet
+!
+ ZVEC1(:) = PLBDH(:)
+ ZVEC2(:) = PLBDS(:)
+!
+!* find the next lower indice for the ZLBDAG and for the ZLBDAS
+! in the geometrical set of (Lbda_g,Lbda_s) couplet use to
+! tabulate the SDRYG-kernel
+!
+ ZVEC1(:) = MAX( 1.0001, MIN( REAL(NWETLBDAH)-0.0001, &
+ XWETINTP1H * LOG( ZVEC1(:) ) + XWETINTP2H ) )
+ IVEC1(:) = INT( ZVEC1(:) )
+ ZVEC1(:) = ZVEC1(:) - REAL( IVEC1(:) )
+!
+ ZVEC2(:) = MAX( 1.0001, MIN( REAL(NWETLBDAS)-0.0001, &
+ XWETINTP1S * LOG( ZVEC2(:) ) + XWETINTP2S ) )
+ IVEC2(:) = INT( ZVEC2(:) )
+ ZVEC2(:) = ZVEC2(:) - REAL( IVEC2(:) )
+!
+!* perform the bilinear interpolation of the normalized
+! SDRYG-kernel
+ !
+ Z1(:) = GET_XKER_SWETH(IVEC1(:)+1,IVEC2(:)+1)
+ Z2(:) = GET_XKER_SWETH(IVEC1(:)+1,IVEC2(:) )
+ Z3(:) = GET_XKER_SWETH(IVEC1(:) ,IVEC2(:)+1)
+ Z4(:) = GET_XKER_SWETH(IVEC1(:) ,IVEC2(:) )
+ ZVEC3(:) = ( Z1(:)* ZVEC2(:) &
+ - Z2(:)*(ZVEC2(:) - 1.0) ) &
+ * ZVEC1(:) &
+ - ( Z3(:)* ZVEC2(:) &
+ - Z4(:)*(ZVEC2(:) - 1.0) ) &
+ * (ZVEC1(:) - 1.0)
+ ZZW(:) = ZVEC3(:)
+ ZZW3(:) = XFSWETH * ZZW(:) & ! RSWETH
+ * PRST(:) * PCHT(:) &
+ * PRHODREF(:)**(1-XCEXVT) &
+ *( XLBSWETH1/( PLBDH(:)**2 ) + &
+ XLBSWETH2/( PLBDH(:) * PLBDS(:) ) + &
+ XLBSWETH3/( PLBDS(:)**2) )
+!
+ Z1(:) = GET_XKER_N_SWETH(IVEC1(:)+1,IVEC2(:)+1)
+ Z2(:) = GET_XKER_N_SWETH(IVEC1(:)+1,IVEC2(:) )
+ Z3(:) = GET_XKER_N_SWETH(IVEC1(:) ,IVEC2(:)+1)
+ Z4(:) = GET_XKER_N_SWETH(IVEC1(:) ,IVEC2(:) )
+ ZVEC3(:) = ( Z1(:)* ZVEC2(:) &
+ - Z2(:)*(ZVEC2(:) - 1.0) ) &
+ * ZVEC1(:) &
+ - ( Z3(:)* ZVEC2(:) &
+ - Z4(:)*(ZVEC2(:) - 1.0) ) &
+ * (ZVEC1(:) - 1.0)
+ ZZW(:) = ZVEC3(:)
+ ZZW3N(:) = XFNSWETH * ZZW(:) & ! NSWETH
+ * PCST(:) * PCHT(:) &
+ * PRHODREF(:)**(1-XCEXVT) &
+ *( XLBNSWETH1/( PLBDH(:)**2 ) + &
+ XLBNSWETH2/( PLBDH(:) * PLBDS(:) ) + &
+ XLBNSWETH3/( PLBDS(:)**2) )
+END WHERE
+!
+!* 1.c Collection of rg
+! ---------------------
+!
+GWET(:) = PRGT(:)>XRTMIN(6) .AND. PRHT(:)>XRTMIN(7) .AND. LDCOMPUTE(:) .AND. &
+ PCGT(:)>XCTMIN(6) .AND. PCHT(:)>XCTMIN(7)
+!
+WHERE( GWET )
+!
+!* Select the (ZLBDAG,ZLBDAR) couplet
+!
+ ZVEC1(:) = PLBDH(:)
+ ZVEC2(:) = PLBDG(:)
+!
+!* Find the next lower indice for the ZLBDAG and for the ZLBDAR
+! in the geometrical set of (Lbda_g,Lbda_r) couplet use to
+! tabulate the RDRYG-kernel
+!
+ ZVEC1(:) = MAX( 1.0001, MIN( REAL(NWETLBDAH)-0.0001, &
+ XWETINTP1H * LOG( ZVEC1(:) ) + XWETINTP2H ) )
+ IVEC1(:) = INT( ZVEC1(:) )
+ ZVEC1(:) = ZVEC1(:) - REAL( IVEC1(:) )
+!
+ ZVEC2(:) = MAX( 1.0001, MIN( REAL(NWETLBDAG)-0.0001, &
+ XWETINTP1G * LOG( ZVEC2(:) ) + XWETINTP2G ) )
+ IVEC2(:) = INT( ZVEC2(:) )
+ ZVEC2(:) = ZVEC2(:) - REAL( IVEC2(:) )
+!
+!* Perform the bilinear interpolation of the normalized
+! RDRYG-kernel
+!
+ Z1(:) = GET_XKER_GWETH(IVEC1(:)+1,IVEC2(:)+1)
+ Z2(:) = GET_XKER_GWETH(IVEC1(:)+1,IVEC2(:) )
+ Z3(:) = GET_XKER_GWETH(IVEC1(:) ,IVEC2(:)+1)
+ Z4(:) = GET_XKER_GWETH(IVEC1(:) ,IVEC2(:) )
+ ZVEC3(:) = ( Z1(:)* ZVEC2(:) &
+ - Z2(:)*(ZVEC2(:) - 1.0) ) &
+ * ZVEC1(:) &
+ - ( Z3(:)* ZVEC2(:) &
+ - Z4(:)*(ZVEC2(:) - 1.0) ) &
+ * (ZVEC1(:) - 1.0)
+ ZZW(:) = ZVEC3(:)
+ ZZW4(:) = XFGWETH * ZZW(:) & ! RGWETH
+ * PRGT(:) * PCHT(:) &
+ * PRHODREF(:)**(1-XCEXVT) &
+ *( XLBGWETH1/( PLBDH(:)**2 ) + &
+ XLBGWETH2/( PLBDH(:) * PLBDG(:) ) + &
+ XLBGWETH3/( PLBDG(:)**2) )
+!
+ Z1(:) = GET_XKER_N_GWETH(IVEC1(:)+1,IVEC2(:)+1)
+ Z2(:) = GET_XKER_N_GWETH(IVEC1(:)+1,IVEC2(:) )
+ Z3(:) = GET_XKER_N_GWETH(IVEC1(:) ,IVEC2(:)+1)
+ Z4(:) = GET_XKER_N_GWETH(IVEC1(:) ,IVEC2(:) )
+ ZVEC3(:) = ( Z1(:)* ZVEC2(:) &
+ - Z2(:)*(ZVEC2(:) - 1.0) ) &
+ * ZVEC1(:) &
+ - ( Z3(:)* ZVEC2(:) &
+ - Z4(:)*(ZVEC2(:) - 1.0) ) &
+ * (ZVEC1(:) - 1.0)
+ ZZW(:) = ZVEC3(:)
+ ZZW4N(:) = XFNGWETH * ZZW(:) & ! NGWETH
+ * PCGT(:) * PCHT(:) &
+ * PRHODREF(:)**(1-XCEXVT) &
+ *( XLBNGWETH1/( PLBDH(:)**2 ) + &
+ XLBNGWETH2/( PLBDH(:) * PLBDG(:) ) + &
+ XLBNGWETH3/( PLBDG(:)**2) )
+
+END WHERE
+!
+! 1.e Wet growth of hail
+! ----------------------
+!
+ZZW(:) = 0.0
+WHERE( PRGT(:)>XRTMIN(6) .AND. PCGT(:) >XCTMIN(6) .AND. LDCOMPUTE(:) )
+ ZZW(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
+ ZZW(:) = PKA(:)*(XTT-PT(:)) + &
+ ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PT(:) - XTT )) &
+ *(XESTT-ZZW(:))/(XRV*PT(:)) )
+!
+! Total mass gained by hail in wet mode
+ ZRWETH(:) = MAX( 0.0, &
+ ( ZZW(:) * PCHT(:) * ( X0DEPH* PLBDH(:)**XEX0DEPH + &
+ X1DEPH*PCJ(:)*PLBDH(:)**XEX1DEPH ) + &
+ ( ZZW2(:)+ZZW3(:)+ZZW4(:) ) * ( XLMTT + (XCI-XCL)*(XTT-PT(:)) ) ) &
+ / (XLMTT-XCL*(XTT-PT(:))) )
+ ! We must agregate, at least, the cold species
+ ZRWETH(:)=MAX(ZRWETH(:), ZZW2(:)+ZZW3(:)+ZZW4(:))
+ ! Mass of rain and cloud droplets frozen by hail (RCWETH + RRWETH)
+ ZZW5(:) = ZRWETH(:) - ZZW2(:) - ZZW3(:) - ZZW4(:)
+END WHERE
+!
+ZZW(:) = 0.0
+WHERE( LDCOMPUTE(:) .AND. PT(:)<XTT .AND. ZZW5(:)>0.0 )
+ P_RC_WETH(:) = - ZZW1(:)
+ P_CC_WETH(:) = P_RC_WETH(:) * PCCT(:)/MAX(PRCT(:),XRTMIN(2))
+ P_RR_WETH(:) = - ZZW5(:) + ZZW1(:)
+ P_CR_WETH(:) = P_RR_WETH(:) * PCRT(:)/MAX(PRRT(:),XRTMIN(3))
+ P_RI_WETH(:) = - ZZW2(:)
+ P_CI_WETH(:) = P_RI_WETH(:) * PCIT(:)/MAX(PRIT(:),XRTMIN(4))
+ P_RS_WETH(:) = - ZZW3(:)
+ P_CS_WETH(:) = - ZZW3N(:)
+ P_RG_WETH(:) = - ZZW4(:)
+ P_CG_WETH(:) = - ZZW4N(:)
+ P_RH_WETH(:) = ZRWETH(:)
+ !
+ P_TH_WETH(:) = ZZW5(:) * (PLSFACT(:)-PLVFACT(:))
+END WHERE
+!
+!
+!* 2. Hail -> graupel conversion (COHG)
+! -------------------------------------
+!
+ZTHRH=0.01E-3
+ZTHRC=0.001E-3
+ZZW(:) = 0.0
+GWET(:) = PRHT(:)<ZTHRH .AND. PRCT(:)<ZTHRC .AND. PT(:)<XTT
+WHERE( GWET(:) )
+ P_RG_COHG(:) = PRHT * MIN( 1.0,MAX( 0.0,1.0-(PRCT(:)/ZTHRC) ) )
+ P_CG_COHG(:) = P_RG_COHG(:) * PCHT(:)/MAX(PRHT(:),XRTMIN(7))
+END WHERE
+!
+!
+!* 3. Hail Melting (HMLT)
+! -----------------------
+!
+ZZW(:) = 0.0
+WHERE( PRHT(:)>XRTMIN(6) .AND. PCHT(:)>XCTMIN(7) .AND. PT(:)>XTT .AND. LDCOMPUTE(:) )
+ ZZW(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
+ ZZW(:) = PKA(:)*(XTT-PT(:)) + &
+ ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PT(:) - XTT )) &
+ *(XESTT-ZZW(:))/(XRV*PT(:)) )
+!
+! compute RHMLTR
+!
+ ZZW(:) = MAX( 0.0,( -ZZW(:) * PCHT(:) * &
+ ( X0DEPH* PLBDH(:)**XEX0DEPH + &
+ X1DEPH*PCJ(:)*PLBDH(:)**XEX1DEPH ) - &
+ ( ZZW5(:) ) * ( XCL*(XTT-PT(:))) ) &
+ / XLMTT )
+ P_RR_HMLT(:) = ZZW(:)
+ P_CR_HMLT(:) = ZZW(:) * 5.0E6 ! obtained after averaging, Dshed=1mm and 500 microns
+ P_CH_HMLT(:) = - ZZW(:) * PCHT(:)/PRHT(:)
+ !
+ P_TH_HMLT(:) = - P_RR_HMLT(:) * (PLSFACT(:)-PLVFACT(:))
+END WHERE
+!
+!
+!
+!
+PA_RC(:) = PA_RC(:) + P_RC_WETH(:)
+PA_CC(:) = PA_CC(:) + P_CC_WETH(:)
+PA_RR(:) = PA_RR(:) + P_RR_WETH(:) + P_RR_HMLT(:)
+PA_CR(:) = PA_CR(:) + P_CR_WETH(:) + P_CR_HMLT(:)
+PA_RI(:) = PA_RI(:) + P_RI_WETH(:)
+PA_CI(:) = PA_CI(:) + P_CI_WETH(:)
+PA_RS(:) = PA_RS(:) + P_RS_WETH(:)
+PA_CS(:) = PA_CS(:) + P_CS_WETH(:)
+PA_RG(:) = PA_RG(:) + P_RG_WETH(:) + P_RG_COHG(:)
+PA_CG(:) = PA_CG(:) + P_CG_WETH(:) + P_CG_COHG(:)
+PA_RH(:) = PA_RH(:) + P_RH_WETH(:) - P_RG_COHG(:) - P_RR_HMLT(:)
+PA_CH(:) = PA_CH(:) - P_CG_COHG(:) + P_CH_HMLT(:)
+PA_TH(:) = PA_TH(:) + P_TH_WETH(:) + P_TH_HMLT(:)
+!
+!-------------------------------------------------------------------------------
+!
+CONTAINS
+ FUNCTION GET_XKER_SWETH(GRAUPEL,SNOW) RESULT(RET)
+ INTEGER, DIMENSION(:) :: GRAUPEL
+ INTEGER, DIMENSION(:) :: SNOW
+ REAL, DIMENSION(SIZE(SNOW)) :: RET
+ !
+ INTEGER I
+ !
+ DO I=1,SIZE(GRAUPEL)
+ RET(I) = XKER_SWETH(MAX(MIN(GRAUPEL(I),SIZE(XKER_SWETH,1)),1),MAX(MIN(SNOW(I),SIZE(XKER_SWETH,2)),1))
+ END DO
+ END FUNCTION GET_XKER_SWETH
+!
+!-------------------------------------------------------------------------------
+!
+ FUNCTION GET_XKER_N_SWETH(GRAUPEL,SNOW) RESULT(RET)
+ INTEGER, DIMENSION(:) :: GRAUPEL
+ INTEGER, DIMENSION(:) :: SNOW
+ REAL, DIMENSION(SIZE(SNOW)) :: RET
+ !
+ INTEGER I
+ !
+ DO I=1,SIZE(GRAUPEL)
+ RET(I) = XKER_N_SWETH(MAX(MIN(GRAUPEL(I),SIZE(XKER_N_SWETH,1)),1),MAX(MIN(SNOW(I),SIZE(XKER_N_SWETH,2)),1))
+ END DO
+ END FUNCTION GET_XKER_N_SWETH
+!
+!-------------------------------------------------------------------------------
+!
+ FUNCTION GET_XKER_GWETH(GRAUPEL,SNOW) RESULT(RET)
+ INTEGER, DIMENSION(:) :: GRAUPEL
+ INTEGER, DIMENSION(:) :: SNOW
+ REAL, DIMENSION(SIZE(SNOW)) :: RET
+ !
+ INTEGER I
+ !
+ DO I=1,SIZE(GRAUPEL)
+ RET(I) = XKER_GWETH(MAX(MIN(GRAUPEL(I),SIZE(XKER_GWETH,1)),1),MAX(MIN(SNOW(I),SIZE(XKER_GWETH,2)),1))
+ END DO
+ END FUNCTION GET_XKER_GWETH
+!
+!-------------------------------------------------------------------------------
+!
+ FUNCTION GET_XKER_N_GWETH(GRAUPEL,SNOW) RESULT(RET)
+ INTEGER, DIMENSION(:) :: GRAUPEL
+ INTEGER, DIMENSION(:) :: SNOW
+ REAL, DIMENSION(SIZE(SNOW)) :: RET
+ !
+ INTEGER I
+ !
+ DO I=1,SIZE(GRAUPEL)
+ RET(I) = XKER_N_GWETH(MAX(MIN(GRAUPEL(I),SIZE(XKER_N_GWETH,1)),1),MAX(MIN(SNOW(I),SIZE(XKER_N_GWETH,2)),1))
+ END DO
+ END FUNCTION GET_XKER_N_GWETH
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE LIMA_HAIL
diff --git a/src/MNH/lima_hail_deposition.f90 b/src/MNH/lima_hail_deposition.f90
new file mode 100644
index 0000000000000000000000000000000000000000..1b411138fdfc344c7ab8cdc4043bc37e737baa2e
--- /dev/null
+++ b/src/MNH/lima_hail_deposition.f90
@@ -0,0 +1,100 @@
+!MNH_LIC Copyright 2018-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-------------------------------------------------------------------------------
+! #################################
+ MODULE MODI_LIMA_HAIL_DEPOSITION
+! #################################
+!
+INTERFACE
+ SUBROUTINE LIMA_HAIL_DEPOSITION (LDCOMPUTE, PRHODREF, &
+ PRHT, PCHT, PSSI, PLBDH, PAI, PCJ, PLSFACT, &
+ P_TH_DEPH, P_RH_DEPH )
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHT ! hail mr
+REAL, DIMENSION(:), INTENT(IN) :: PCHT ! hail conc
+REAL, DIMENSION(:), INTENT(IN) :: PSSI !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDH !
+REAL, DIMENSION(:), INTENT(IN) :: PAI !
+REAL, DIMENSION(:), INTENT(IN) :: PCJ !
+REAL, DIMENSION(:), INTENT(IN) :: PLSFACT !
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_DEPH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RH_DEPH
+!!
+END SUBROUTINE LIMA_HAIL_DEPOSITION
+END INTERFACE
+END MODULE MODI_LIMA_HAIL_DEPOSITION
+!
+! ###########################################################################
+ SUBROUTINE LIMA_HAIL_DEPOSITION (LDCOMPUTE, PRHODREF, &
+ PRHT, PCHT, PSSI, PLBDH, PAI, PCJ, PLSFACT, &
+ P_TH_DEPH, P_RH_DEPH )
+! ###########################################################################
+!
+!! PURPOSE
+!! -------
+!! Deposition of water vapour on graupel
+!!
+!!
+!! AUTHOR
+!! ------
+!! J.-M. Cohard * Laboratoire d'Aerologie*
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!! S. Berthet * Laboratoire d'Aerologie*
+!! B. Vié * CNRM *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 15/03/2018
+!!
+! M. Taufour 07/2022 add concentration for snow, graupel, hail
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAM_LIMA, ONLY : XRTMIN, XCTMIN
+USE MODD_PARAM_LIMA_MIXED, ONLY : X0DEPH, XEX0DEPH, X1DEPH, XEX1DEPH
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF !
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHT ! hail mr
+REAL, DIMENSION(:), INTENT(IN) :: PCHT ! hail conc
+REAL, DIMENSION(:), INTENT(IN) :: PSSI !
+REAL, DIMENSION(:), INTENT(IN) :: PLBDH !
+REAL, DIMENSION(:), INTENT(IN) :: PAI !
+REAL, DIMENSION(:), INTENT(IN) :: PCJ !
+REAL, DIMENSION(:), INTENT(IN) :: PLSFACT !
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_DEPH
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RH_DEPH
+!
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. Deposition of vapour on graupel
+! -------------------------------
+!
+P_TH_DEPH(:) = 0.0
+P_RH_DEPH(:) = 0.0
+WHERE ( PRHT(:)>XRTMIN(7) .AND. PCHT(:)>XCTMIN(7) .AND. LDCOMPUTE(:) )
+ P_RH_DEPH(:) = PSSI(:) / PAI(:) * PCHT(:) * &
+ ( X0DEPH*PLBDH(:)**XEX0DEPH + X1DEPH*PCJ(:)*PLBDH(:)**XEX1DEPH )
+ P_TH_DEPH(:) = P_RH_DEPH(:)*PLSFACT(:)
+END WHERE
+!
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE LIMA_HAIL_DEPOSITION
diff --git a/src/MNH/lima_snow_self_collection.f90 b/src/MNH/lima_snow_self_collection.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9c499d66648bb390997c5918487a4d2d46ed69fc
--- /dev/null
+++ b/src/MNH/lima_snow_self_collection.f90
@@ -0,0 +1,147 @@
+!MNH_LIC Copyright 2018-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-------------------------------------------------------------------------------
+! #################################
+ MODULE MODI_LIMA_SNOW_SELF_COLLECTION
+! #################################
+!
+INTERFACE
+ SUBROUTINE LIMA_SNOW_SELF_COLLECTION (LDCOMPUTE, &
+ PRHODREF, PT, &
+ PCST, PLBDS, &
+ P_CS_SSC )
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF ! Reference Exner function
+REAL, DIMENSION(:), INTENT(IN) :: PT ! Temperature
+!
+REAL, DIMENSION(:), INTENT(IN) :: PCST ! Snow C. at t
+REAL, DIMENSION(:), INTENT(IN) :: PLBDS !
+!
+REAL, DIMENSION(:), INTENT(OUT) :: P_CS_SSC
+!
+END SUBROUTINE LIMA_SNOW_SELF_COLLECTION
+END INTERFACE
+END MODULE MODI_LIMA_SNOW_SELF_COLLECTION
+!
+! #############################################################
+ SUBROUTINE LIMA_SNOW_SELF_COLLECTION (LDCOMPUTE, &
+ PRHODREF, PT, &
+ PCST, PLBDS, &
+ P_CS_SSC )
+! #############################################################
+!
+!! PURPOSE
+!! -------
+!! Compute the self-collection and physical break-up of snow
+!!
+!!
+!! AUTHOR
+!! ------
+!! J.-M. Cohard * Laboratoire d'Aerologie*
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!! S. Berthet * Laboratoire d'Aerologie*
+!! B. Vié * CNRM *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 15/03/2018
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST, ONLY : XTT
+USE MODD_PARAM_LIMA, ONLY : XCTMIN, XCEXVT
+USE MODD_PARAM_LIMA_COLD, ONLY : NSCLBDAS, XSCINTP1S, XSCINTP2S, XKER_N_SSCS, XFNSSCS, XCOLEXSS, &
+ XLBNSSCS1, XLBNSSCS2
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF ! Reference Exner function
+REAL, DIMENSION(:), INTENT(IN) :: PT ! Temperature
+!
+REAL, DIMENSION(:), INTENT(IN) :: PCST ! Snow C. at t
+REAL, DIMENSION(:), INTENT(IN) :: PLBDS !
+!
+REAL, DIMENSION(:), INTENT(OUT) :: P_CS_SSC
+!
+!* 0.2 Declarations of local variables :
+!
+REAL, DIMENSION(SIZE(PCST)) :: &
+ ZW1, & ! work arrays
+ ZW2
+LOGICAL, DIMENSION(SIZE(PCST)) :: GSSC
+INTEGER :: IGSSC, JJ
+INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC1,IVEC2 ! Vectors of indices
+REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1,ZVEC2, ZVEC3 ! Work vectors
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. Snow self-collection and break-up
+! ---------------------------------
+!
+!
+P_CS_SSC(:)=0.
+!
+ZW1(:) =0.
+ZW2(:) =0.
+!
+GSSC(:) = PCST(:)>XCTMIN(5)
+IGSSC = COUNT(GSSC(:))
+!
+IF( IGSSC>0 ) THEN
+!
+! 1.3N.0 allocations
+!
+ ALLOCATE(ZVEC1(IGSSC))
+ ALLOCATE(IVEC1(IGSSC))
+!
+! 1.3N.1 select the (ZLBDAS,ZLBDAS) couplet
+!
+ ZVEC1(:) = PACK( PLBDS(:),MASK=GSSC(:) )
+!
+! 1.3N.2 find the next lower indice for the ZLBDAS and for the ZLBDAS
+! in the geometrical set of (Lbda_s,Lbda_s) couplet use to
+! tabulate the SACCS-kernel
+!
+ ZVEC1(1:IGSSC) = MAX( 1.0001, MIN( FLOAT(NSCLBDAS)-0.0001, &
+ XSCINTP1S * LOG( ZVEC1(1:IGSSC) ) + XSCINTP2S ) )
+ IVEC1(1:IGSSC) = INT( ZVEC1(1:IGSSC) )
+ ZVEC1(1:IGSSC) = ZVEC1(1:IGSSC) - FLOAT( IVEC1(1:IGSSC) )
+!
+! 1.3N.3 perform the bilinear interpolation of the normalized
+! SSCS-kernel
+!
+ ALLOCATE(ZVEC3(IGSSC))
+ DO JJ = 1,IGSSC
+ ZVEC3(JJ) = ( XKER_N_SSCS(IVEC1(JJ)+1,IVEC1(JJ)+1)* ZVEC1(JJ) &
+ - XKER_N_SSCS(IVEC1(JJ)+1,IVEC1(JJ) )*(ZVEC1(JJ) - 1.0) ) &
+ * ZVEC1(JJ) &
+ - ( XKER_N_SSCS(IVEC1(JJ) ,IVEC1(JJ)+1)* ZVEC1(JJ) &
+ - XKER_N_SSCS(IVEC1(JJ) ,IVEC1(JJ) )*(ZVEC1(JJ) - 1.0) ) &
+ * (ZVEC1(JJ) - 1.0)
+ END DO
+ ZW1(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GSSC(:),FIELD=0.0 ) !! NSACCS
+ DEALLOCATE(ZVEC3)
+!
+ WHERE( GSSC(:) )
+ P_CS_SSC(:) = - XFNSSCS * ZW1(:) * EXP( XCOLEXSS*(PT(:)-XTT) ) * PCST(:)**2 &
+ * PRHODREF(:)**(-XCEXVT-1.) * (XLBNSSCS1+XLBNSSCS2) / PLBDS(:)**2
+ END WHERE
+ DEALLOCATE(IVEC1)
+ DEALLOCATE(ZVEC1)
+END IF
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE LIMA_SNOW_SELF_COLLECTION
diff --git a/src/MNH/nrcolss.f90 b/src/MNH/nrcolss.f90
new file mode 100644
index 0000000000000000000000000000000000000000..d21be2bd3b12000a5c5871c47e9ba5faa22d84f5
--- /dev/null
+++ b/src/MNH/nrcolss.f90
@@ -0,0 +1,316 @@
+!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
+!
+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/src/MNH/nscolrg.f90 b/src/MNH/nscolrg.f90
new file mode 100644
index 0000000000000000000000000000000000000000..790a01f76a32fef7603ff99d7ad1339a657fbb80
--- /dev/null
+++ b/src/MNH/nscolrg.f90
@@ -0,0 +1,317 @@
+!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
+!
+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/src/MNH/nzcolx.f90 b/src/MNH/nzcolx.f90
new file mode 100644
index 0000000000000000000000000000000000000000..e4493c2e60f617e81998e2466d534e887891feaf
--- /dev/null
+++ b/src/MNH/nzcolx.f90
@@ -0,0 +1,278 @@
+!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/05/18 13:07:25
+!-----------------------------------------------------------------
+! ##################
+ MODULE MODI_NZCOLX
+! ##################
+!
+INTERFACE
+!
+ SUBROUTINE NZCOLX( KND, PALPHAX, PNUX, PALPHAZ, PNUZ, &
+ PEXZ, PFALLX, PEXFALLX, PFALLEXPX, &
+ PFALLZ, PEXFALLZ, PFALLEXPZ, &
+ PLBDAXMAX, PLBDAZMAX, PLBDAXMIN, PLBDAZMIN, &
+ PDINFTY, PNZCOLX )
+!
+INTEGER, INTENT(IN) :: KND ! Number of discrete size intervals in DX and DZ
+!
+!
+REAL, INTENT(IN) :: PALPHAX ! First shape parameter of the specy X
+ ! size distribution (generalized gamma law)
+REAL, INTENT(IN) :: PNUX ! Second shape parameter of the specy X
+ ! size distribution (generalized gamma law)
+REAL, INTENT(IN) :: PALPHAZ ! First shape parameter of the specy Z
+ ! size distribution (generalized gamma law)
+REAL, INTENT(IN) :: PNUZ ! Second shape parameter of the specy Z
+REAL, INTENT(IN) :: PEXZ ! Efficiency of specy X collecting specy Z
+REAL, INTENT(IN) :: PFALLX ! Fall speed constant of specy X
+REAL, INTENT(IN) :: PEXFALLX ! Fall speed exponent of specy X
+REAL, INTENT(IN) :: PFALLEXPX ! Fall speed exponential constant of specy X
+REAL, INTENT(IN) :: PFALLZ ! Fall speed constant of specy Z
+REAL, INTENT(IN) :: PEXFALLZ ! Fall speed exponent of specy Z
+REAL, INTENT(IN) :: PFALLEXPZ ! Fall speed exponential constant of specy Z
+REAL, INTENT(IN) :: PLBDAXMAX ! Maximun slope of size distribution of specy X
+REAL, INTENT(IN) :: PLBDAZMAX ! Maximun slope of size distribution of specy Z
+REAL, INTENT(IN) :: PLBDAXMIN ! Minimun slope of size distribution of specy X
+REAL, INTENT(IN) :: PLBDAZMIN ! Minimun slope of size distribution of specy Z
+REAL, INTENT(IN) :: PDINFTY ! Factor to define the largest diameter up to
+ ! which the diameter integration is performed
+!
+REAL, DIMENSION(:,:), INTENT(INOUT) :: PNZCOLX ! Scaled fall speed difference in
+ ! the mass collection kernel as a
+ ! function of LAMBDAX and LAMBDAZ
+!
+ END SUBROUTINE NZCOLX
+!
+END INTERFACE
+!
+ END MODULE MODI_NZCOLX
+! ################################################################
+ SUBROUTINE NZCOLX( KND, PALPHAX, PNUX, PALPHAZ, PNUZ, &
+ PEXZ, PFALLX, PEXFALLX, PFALLEXPX, &
+ PFALLZ, PEXFALLZ, PFALLEXPZ, &
+ PLBDAXMAX, PLBDAZMAX, PLBDAXMIN, PLBDAZMIN, &
+ PDINFTY, PNZCOLX )
+! ################################################################
+!
+!
+!
+!!**** * - Build up a look-up table containing the scaled fall speed
+!! difference between size distributed particles of specy X and Z
+!!
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this routine is to integrate numerically the scaled fall
+!! speed difference between specy X and specy Z for use in collection
+!! kernels. A first integral of the form
+!!
+!! infty
+!! / /
+!! |{| }
+!! |{| E_xz (Dx+Dz)^2 |cxDx^dx-czDz^dz| g(Dz) dDz} g(Dx) dDx
+!! |{| }
+!! / /
+!! 0
+!!
+!! is evaluated and normalised by a second integral of the form
+!!
+!! infty
+!! / /
+!! |{| }
+!! |{| (Dx+Dz)^2 g(Dz) dDz} g(Dx) dDx
+!! |{| }
+!! / /
+!! 0
+!!
+!! where E_xz is a collection efficiency, g(D) is the generalized Gamma
+!! distribution law. The 'infty' diameter is defined according to the
+!! current value of the Lbda that is D_x=PDINFTY/Lbda_x or
+!! D_z=PINFTY/Lbda_z.
+!! The result is stored in a two-dimensional array.
+!!
+!!** METHOD
+!! ------
+!! The free parameters of the size distribution function of specy X 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 and
+!! the [0,infty] interval is discretized over KND values of D_x or D_z.
+!!
+!! EXTERNAL
+!! --------
+!! MODI_GENERAL_GAMMA: Generalized gamma distribution law
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! None
+!!
+!! 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 rzcolx for concentration
+!! J. Wurtz 03/2022: new snow characteristics
+!!
+!-------------------------------------------------------------------------------
+!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_GENERAL_GAMMA
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 Declarations of dummy arguments
+! -------------------------------
+!
+!
+INTEGER, INTENT(IN) :: KND ! Number of discrete size intervals in DX and DZ
+!
+!
+REAL, INTENT(IN) :: PALPHAX ! First shape parameter of the specy X
+ ! size distribution (generalized gamma law)
+REAL, INTENT(IN) :: PNUX ! Second shape parameter of the specy X
+ ! size distribution (generalized gamma law)
+REAL, INTENT(IN) :: PALPHAZ ! First shape parameter of the specy Z
+ ! size distribution (generalized gamma law)
+REAL, INTENT(IN) :: PNUZ ! Second shape parameter of the specy Z
+ ! size distribution (generalized gamma law)
+REAL, INTENT(IN) :: PEXZ ! Efficiency of specy X collecting specy Z
+REAL, INTENT(IN) :: PFALLX ! Fall speed constant of specy X
+REAL, INTENT(IN) :: PEXFALLX ! Fall speed exponent of specy X
+REAL, INTENT(IN) :: PFALLEXPX ! Fall speed exponential constant of specy X
+REAL, INTENT(IN) :: PFALLZ ! Fall speed constant of specy Z
+REAL, INTENT(IN) :: PEXFALLZ ! Fall speed exponent of specy Z
+REAL, INTENT(IN) :: PFALLEXPZ ! Fall speed exponential constant of specy Z
+REAL, INTENT(IN) :: PLBDAXMAX ! Maximun slope of size distribution of specy X
+REAL, INTENT(IN) :: PLBDAZMAX ! Maximun slope of size distribution of specy Z
+REAL, INTENT(IN) :: PLBDAXMIN ! Minimun slope of size distribution of specy X
+REAL, INTENT(IN) :: PLBDAZMIN ! Minimun slope of size distribution of specy Z
+REAL, INTENT(IN) :: PDINFTY ! Factor to define the largest diameter up to
+ ! which the diameter integration is performed
+!
+REAL, DIMENSION(:,:), INTENT(INOUT) :: PNZCOLX ! Scaled fall speed difference in
+ ! the mass collection kernel as a
+ ! function of LAMBDAX and LAMBDAZ
+!
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+!
+INTEGER :: JLBDAX ! Slope index of the size distribution of specy X
+INTEGER :: JLBDAZ ! Slope index of the size distribution of specy Z
+INTEGER :: JDX ! Diameter index of a particle of specy X
+INTEGER :: JDZ ! Diameter index of a particle of specy Z
+!
+!
+REAL :: ZLBDAX ! Current slope parameter LAMBDA of specy X
+REAL :: ZLBDAZ ! Current slope parameter LAMBDA of specy Z
+REAL :: ZDLBDAX ! Growth rate of the slope parameter LAMBDA of specy X
+REAL :: ZDLBDAZ ! Growth rate of the slope parameter LAMBDA of specy Z
+REAL :: ZDDX ! Integration step of the diameter of specy X
+REAL :: ZDDZ ! Integration step of the diameter of specy Z
+REAL :: ZDX ! Current diameter of the particle specy X
+REAL :: ZDZ ! Current diameter of the particle specy Z
+REAL :: ZCOLLZ ! Single integral of the mass weighted fall speed difference
+ ! over the spectrum of specy Z
+REAL :: ZCOLLXZ ! Double integral of the mass weighted fall speed difference
+ ! over the spectra of specy X and specy Z
+REAL :: ZSCALZ ! Single integral of the scaling factor over
+ ! the spectrum of specy Z
+REAL :: ZSCALXZ ! Double integral of the scaling factor over
+ ! the spectra of specy X and specy Z
+REAL :: ZFUNC ! Ancillary function
+!
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1 COMPUTE THE SCALED VELOCITZ DIFFERENCE IN THE MASS
+!* COLLECTION KERNEL,
+! -------------------------------------------------
+!
+!
+!
+!* 1.1 Compute the growth rate of the slope factors LAMBDA
+!
+ZDLBDAX = EXP( LOG(PLBDAXMAX/PLBDAXMIN)/REAL(SIZE(PNZCOLX(:,:),1)-1) )
+ZDLBDAZ = EXP( LOG(PLBDAZMAX/PLBDAZMIN)/REAL(SIZE(PNZCOLX(:,:),2)-1) )
+!
+!* 1.2 Scan the slope factors LAMBDAX and LAMBDAZ
+!
+DO JLBDAX = 1,SIZE(PNZCOLX(:,:),1)
+ ZLBDAX = PLBDAXMIN * ZDLBDAX ** (JLBDAX-1)
+ DO JLBDAZ = 1,SIZE(PNZCOLX(:,:),2)
+ ZLBDAZ = PLBDAZMIN * ZDLBDAZ ** (JLBDAZ-1)
+!
+!* 1.3 Initialize the collection integrals
+!
+ ZSCALXZ = 0.0
+ ZCOLLXZ = 0.0
+!
+!* 1.4 Compute the diameter steps
+!
+ ZDDX = PDINFTY / (REAL(KND) * ZLBDAX)
+ ZDDZ = PDINFTY / (REAL(KND) * ZLBDAZ)
+!
+!* 1.5 Scan over the diameters DX and DZ
+!
+ DO JDX = 1,KND-1
+ ZDX = ZDDX * REAL(JDX)
+!
+ ZSCALZ = 0.0
+ ZCOLLZ = 0.0
+ DO JDZ = 1,KND-1
+ ZDZ = ZDDZ * REAL(JDZ)
+!
+!* 1.6 Compute the normalization factor by integration over the
+! dimensional spectrum of specy Z
+!
+ ZFUNC = (ZDX+ZDZ)**2 * GENERAL_GAMMA(PALPHAZ,PNUZ,ZLBDAZ,ZDZ)
+ ZSCALZ = ZSCALZ + ZFUNC
+!
+!* 1.7 Compute the scaled fall speed difference by integration over
+! the dimensional spectrum of specy Z
+!
+ ZCOLLZ = ZCOLLZ + ZFUNC * PEXZ * ABS( PFALLX*ZDX**PEXFALLX * EXP(-(ZDX*PFALLEXPX)**PALPHAX) &
+ - PFALLZ*ZDZ**PEXFALLZ * EXP(-(ZDZ*PFALLEXPZ)**PALPHAZ))
+ END DO
+!
+!* 1.8 Compute the normalization factor by integration over the
+! dimensional spectrum of specy X
+!
+ ZFUNC = GENERAL_GAMMA(PALPHAX,PNUX,ZLBDAX,ZDX)
+ ZSCALXZ = ZSCALXZ + ZSCALZ * ZFUNC
+!
+!* 1.9 Compute the scaled fall speed difference by integration over
+! the dimensional spectrum of specy X
+!
+ ZCOLLXZ = ZCOLLXZ + ZCOLLZ * ZFUNC
+ END DO
+!
+!* 1.10 Scale the fall speed difference
+!
+ PNZCOLX(JLBDAX,JLBDAZ) = ZCOLLXZ / ZSCALXZ
+ END DO
+END DO
+!
+END SUBROUTINE NZCOLX