diff --git a/src/MNH/aircraft_balloon_evol.f90 b/src/MNH/aircraft_balloon_evol.f90
index 9642bcc4e8de8a731373a2db1fc3b94934abb8ce..43cab9ac901b4752c7356fb9065625063441901a 100644
--- a/src/MNH/aircraft_balloon_evol.f90
+++ b/src/MNH/aircraft_balloon_evol.f90
@@ -144,10 +144,11 @@ USE MODD_NSV,              ONLY : NSV_LIMA_NI,NSV_LIMA_NR,NSV_LIMA_NC
 USE MODD_PARAMETERS
 USE MODD_PARAM_LIMA,       ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                                  XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC
+                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC,LSNOW_T_L=>LSNOW_T
 USE MODD_PARAM_LIMA_COLD,  ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
                                  XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,&
-                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS
+                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,        &
+                                 XLBDAS_MAX,XLBDAS_MIN
 USE MODD_PARAM_LIMA_MIXED, ONLY: XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,&
                                  XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
 USE MODD_PARAM_LIMA_WARM,  ONLY: XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XAR_L=>XAR,&
@@ -318,6 +319,7 @@ REAL                           :: ZAETOT,ZAETMP,ZREFLOC,ZQSCA,ZQBACK,ZQEXT ! tem
 REAL,DIMENSION(:),ALLOCATABLE  :: ZAELOC,ZZMZ ! temporary arrays
 INTEGER                        :: JPTS_GAULAG=7 ! number of points for Gauss-Laguerre quadrature
 REAL                           :: ZLBDA   ! slope distribution parameter
+REAL                           :: ZN   ! number concentration
 REAL                           :: ZFRAC_ICE  ! ice water fraction
 REAL                           :: ZDELTA_EQUIV ! mass-equivalent Gauss-Laguerre point
 REAL                           :: ZFW ! liquid fraction
@@ -1177,7 +1179,17 @@ IF ( TPFLYER%FLY) THEN
                     ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
                   ENDIF
               END SELECT
-              ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+              IF (JLOOP.EQ.5 .AND. CCLOUD=='LIMA' .AND. LSNOW_T_L) THEN
+                 IF (ZTEMPZ(JK)>-10.) THEN
+                    ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+                 ELSE
+                    ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+                 END IF
+                 ZN=XLBS_L*ZRHODREFZ(JK)*ZRZ(JK,JLOOP)*ZLBDA**ZB
+              ELSE
+                 ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+                 ZN=ZCC*ZLBDA**ZCX
+              END IF
               ZREFLOC=0.
               ZAETMP=0.
               DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
@@ -1209,8 +1221,8 @@ IF ( TPFLYER%FLY) THEN
                 ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
                 ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
               END DO
-              ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
-              ZAETMP=ZAETMP  *           XPI    *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
+              ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZN/(4.*GAMMA(ZNU)*.93)
+              ZAETMP=ZAETMP  *           XPI    *ZN/(4.*GAMMA(ZNU))
               TPFLYER%CRARE(IN,JK)=TPFLYER%CRARE(IN,JK)+ZREFLOC
               ZAELOC(JK)=ZAELOC(JK)+ZAETMP
             END IF
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index b8ec5e8e1f18895ebb5815e7510d78b388803191..fdcacc99fae2818fbcc95d09af63b31b6630dff1 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -277,7 +277,7 @@ USE MODD_ALLSTATION_n
 USE MODD_PARAM_LIMA, ONLY : LCOLD, LNUCL, LSEDI, LHHONI, LSNOW, LHAIL, LMEYERS,       &
                             NMOD_IFN, XIFN_CONC, LIFN_HOM, CIFN_SPECIES,              &
                             CINT_MIXING, NMOD_IMM, NIND_SPECIE,                       &
-                            CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,                   &
+                            LSNOW_T, CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,          &
                             XFACTNUC_DEP, XFACTNUC_CON,                               &
                             OWARM=>LWARM, LACTI, ORAIN=>LRAIN, OSEDC=>LSEDC,          &
                             OACTIT=>LACTIT, LBOUND, LSPRO, LADJ, LKHKO,               &
@@ -1026,6 +1026,7 @@ IF (KMI == 1) THEN
   LSEDI  = .TRUE.
   LSNOW  = .TRUE.
   LHAIL  = .FALSE.
+  LSNOW_T = .TRUE.
   CPRISTINE_ICE_LIMA = 'PLAT'
   CHEVRIMED_ICE_LIMA = 'GRAU'
   XFACTNUC_DEP = 1.0  
diff --git a/src/MNH/ini_ice_c1r3.f90 b/src/MNH/ini_ice_c1r3.f90
index b0a35554a88837b36dc5ae44119426abe6a482a6..883e56a1c03a156dcd58f89053a5be119574f0da 100644
--- a/src/MNH/ini_ice_c1r3.f90
+++ b/src/MNH/ini_ice_c1r3.f90
@@ -90,6 +90,7 @@ END MODULE MODI_INI_ICE_C1R3
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -725,15 +726,15 @@ IF( (KACCLBDAS/=NACCLBDAS) .OR. (KACCLBDAR/=NACCLBDAR) .OR. (KND/=IND) .OR. &
     (PACCLBDAS_MIN/=XACCLBDAS_MIN) .OR. (PACCLBDAR_MIN/=XACCLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR,                          &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
   CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR, 0.,                      &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCS                                        )
   CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBS, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBS, XCS, XDS, 0., XCR, XDR,                          &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,XAG, XBS, XAS                         )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -938,7 +939,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAS/=NDRYLBDAS) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAS_MIN/=XDRYLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, XBS, XCG, XDG, XCS, XDS,                              &
+                ZEGS, XBS, XCG, XDG, 0., XCS, XDS, 0.,                      &
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1004,7 +1005,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAR/=NDRYLBDAR) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAR_MIN/=XDRYLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAR, XNUR,                          &
-                ZEGR, XBR, XCG, XDG, XCR, XDR,                              &
+                ZEGR, XBR, XCG, XDG, 0., XCR, XDR, 0.,                      &
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
diff --git a/src/MNH/ini_lima_cold_mixed.f90 b/src/MNH/ini_lima_cold_mixed.f90
index cb427cdb434982b229095adb417eee3d1071b73e..ea9b388c3d4287d04e0e69fe07f77853c6ccc93c 100644
--- a/src/MNH/ini_lima_cold_mixed.f90
+++ b/src/MNH/ini_lima_cold_mixed.f90
@@ -41,6 +41,7 @@ END MODULE MODI_INI_LIMA_COLD_MIXED
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -171,12 +172,27 @@ XF1IS = 0.28
 !
 XAS = 0.02
 XBS = 1.9
+
+IF (LSNOW_T) THEN
+!Cas Gamma generalisee
+XCS = 11.52
+XDS = 0.39
+XFVELOS =0.097
+!Cas MP
+!XCS = 13.2
+!XDS = 0.423       
+!XFVELOS = 25.14
+ELSE
 XCS = 5.
 XDS = 0.27
-!
-XCCS = 5.0
-XCXS = 1.0
-!
+XFVELOS = 0.
+END IF
+
+IF (.NOT. LSNOW_T) THEN
+   XCCS = 5.0
+   XCXS = 1.0
+END IF
+
 XF0S = 0.86
 XF1S = 0.28
 !
@@ -234,8 +250,17 @@ XC1H = 1./2.
 XALPHAI = 3.0  ! Gamma law for the ice crystal volume
 XNUI    = 3.0  ! Gamma law with little dispersion
 !
-XALPHAS = 1.0  ! Exponential law
-XNUS    = 1.0  ! Exponential law
+IF (LSNOW_T) THEN
+!Cas GAMMAGEN
+   XALPHAS = .214   ! Generalized gamma law
+   XNUS    = 43.7   ! Generalized gamma law
+   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
+                            ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )
+ELSE
+   XALPHAS = 1.0  ! Exponential law
+   XNUS    = 1.0  ! Exponential law
+   XTRANS_MP_GAMMAS = 1.
+END IF
 !
 XALPHAG = 1.0  ! Exponential law
 XNUG    = 1.0  ! Exponential law
@@ -248,8 +273,13 @@ XNUH    = 8.0  ! Gamma law with little dispersion
 XLBEXI = 1.0/XBI
 XLBI   = XAI*MOMG(XALPHAI,XNUI,XBI)
 !
-XLBEXS = 1.0/(XCXS-XBS)
-XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS)
+IF (LSNOW_T) THEN
+   XLBEXS = 0. ! Not used
+   XLBS   = 1.0/(XAS*MOMG(XALPHAS,XNUS,XBS))
+ELSE
+   XLBEXS = 1.0/(XCXS-XBS)
+   XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS)
+END IF
 !
 XLBEXG = 1.0/(XCXG-XBG)
 XLBG   = ( XAG*XCCG*MOMG(XALPHAG,XNUG,XBG))**(-XLBEXG)
@@ -266,7 +296,8 @@ IF (GFLAG) THEN
   WRITE(UNIT=ILUOUT0,FMT='(" XLBEXH =",E13.6," XLBH =",E13.6)') XLBEXH,XLBH
 END IF
 !
-XLBDAS_MAX = 500000
+XLBDAS_MAX = 500000. ! LBDAS_MAX doit être compare avec LBDAS avec une forme de Marshall-Palmer
+XLBDAS_MIN = 1000.
 XLBDAG_MAX = 100000.0
 !
 ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc
@@ -312,9 +343,20 @@ XFSEDRI  = XC_I*GAMMA_X0D(XNUI+(XDI+XBI)/XALPHAI)/GAMMA_X0D(XNUI+XBI/XALPHAI)*
 XFSEDCI  = XC_I*GAMMA_X0D(XNUI+XDI/XALPHAI)/GAMMA_X0D(XNUI)*     &
             (ZRHO00)**XCEXVT
 !
-XEXSEDS = (XBS+XDS-XCXS)/(XBS-XCXS)
-XFSEDS  = XCS*XAS*XCCS*MOMG(XALPHAS,XNUS,XBS+XDS)*                         &
+IF (LSNOW_T) THEN
+!HOUZE/HAIC
+   !XEXSEDS = -XDS !(2*XBS+XDS)
+   !XFSEDS  = XCS*MOMG(XALPHAS,XNUS,XBS+XDS)/(MOMG(XALPHAS,XNUS,XBS))    &
+   !            *(ZRHO00)**XCEXVT
+!LH_EXTENDED
+   XEXSEDS = -XDS-XBS
+   XFSEDS  = XCS*MOMG(XALPHAS,XNUS,XBS+XDS)/(MOMG(XALPHAS,XNUS,XBS))    &
+            *(ZRHO00)**XCEXVT
+ELSE
+   XEXSEDS = (XBS+XDS-XCXS)/(XBS-XCXS)
+   XFSEDS  = XCS*XAS*XCCS*MOMG(XALPHAS,XNUS,XBS+XDS)*                         &
             (XAS*XCCS*MOMG(XALPHAS,XNUS,XBS))**(-XEXSEDS)*(ZRHO00)**XCEXVT
+END IF
 !
 XEXSEDG = (XBG+XDG-XCXG)/(XBG-XCXG)
 XFSEDG  = XCG*XAG*XCCG*MOMG(XALPHAG,XNUG,XBG+XDG)*                         &
@@ -324,8 +366,6 @@ XEXSEDH = (XBH+XDH-XCXH)/(XBH-XCXH)
 XFSEDH  = XCH*XAH*XCCH*MOMG(XALPHAH,XNUH,XBH+XDH)*                         &
             (XAH*XCCH*MOMG(XALPHAH,XNUH,XBH))**(-XEXSEDH)*(ZRHO00)**XCEXVT
 !
-!
-!
 XLB(4)    = XLBI
 XLBEX(4)  = XLBEXI
 XD(4)     = XDI
@@ -561,7 +601,7 @@ XR1DEPIS = XC1DEPIS *(XAI*XDICNVS_LIM**XBI)
 !
 ! Harrington parameterization for snow to ice conversion
 !
-XLBDASCNVI_MAX = 6000. ! lbdas max after Field (1999)
+XLBDASCNVI_MAX = 6000.*XTRANS_MP_GAMMAS ! lbdas max after Field (1999)
 !
 XDSCNVI_LIM    = 125.E-6  ! size in microns
 XLBDASCNVI_LIM = (50.0**(1.0/(XALPHAS)))/XDSCNVI_LIM  ! ZLBDAS Limitation
@@ -574,10 +614,10 @@ XR1DEPSI = XC1DEPSI *(XAS*XDSCNVI_LIM**XBS)
 !
 ! Vapor deposition on snow and graupel and hail
 !
-X0DEPS = (4.0*XPI)*XCCS*XC1S*XF0S*MOMG(XALPHAS,XNUS,1.)
-X1DEPS = (4.0*XPI)*XCCS*XC1S*XF1S*SQRT(XCS)*MOMG(XALPHAS,XNUS,0.5*XDS+1.5)
-XEX0DEPS = XCXS-1.0
-XEX1DEPS = XCXS-0.5*(XDS+3.0)
+X0DEPS = XLBS*(4.0*XPI)*XC1S*XF0S*MOMG(XALPHAS,XNUS,1.)
+X1DEPS = XLBS*(4.0*XPI)*XC1S*XF1S*SQRT(XCS)*MOMG(XALPHAS,XNUS,0.5*XDS+1.5)
+XEX0DEPS = XBS-1.0
+XEX1DEPS = -0.5*(XDS+3.0)
 !
 X0DEPG = (4.0*XPI)*XCCG*XC1G*XF0G*MOMG(XALPHAG,XNUG,1.)
 X1DEPG = (4.0*XPI)*XCCG*XC1G*XF1G*SQRT(XCG)*MOMG(XALPHAG,XNUG,0.5*XDG+1.5)
@@ -669,11 +709,12 @@ END IF
 !
 XDCSLIM  = 0.007 ! D_cs^lim = 7 mm as suggested by Farley et al. (1989)
 XCOLCS   = 1.0
-XEXCRIMSS= XCXS-XDS-2.0
-XCRIMSS  = (XPI/4.0)*XCOLCS*XCCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
+XEXCRIMSS = -XDS-2.0
+XCRIMSS  = XLBS * (XPI/4.0)*XCOLCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
+
 XEXCRIMSG= XEXCRIMSS
 XCRIMSG  = XCRIMSS
-XSRIMCG  = XCCS*XAS*MOMG(XALPHAS,XNUS,XBS)
+XSRIMCG  = XAS*MOMG(XALPHAS,XNUS,XBS)
 XEXSRIMCG= XCXS-XBS
 !
 GFLAG = .TRUE.
@@ -684,9 +725,9 @@ IF (GFLAG) THEN
 END IF
 !
 NGAMINC = 80
-XGAMINC_BOUND_MIN = 1.0E-1 ! Minimal value of (Lbda * D_cs^lim)**alpha
-XGAMINC_BOUND_MAX = 1.0E7  ! Maximal value of (Lbda * D_cs^lim)**alpha
-ZRATE = EXP(LOG(XGAMINC_BOUND_MAX/XGAMINC_BOUND_MIN)/REAL(NGAMINC-1))
+XGAMINC_BOUND_MIN = (1000.*XTRANS_MP_GAMMAS*XDCSLIM)**XALPHAS !1.0E-1 ! Minimal value of (Lbda * D_cs^lim)**alpha
+XGAMINC_BOUND_MAX = (50000.*XTRANS_MP_GAMMAS*XDCSLIM)**XALPHAS !1.0E7 ! Maximal value of (Lbda * D_cs^lim)**alpha
+ZRATE = EXP(LOG(XGAMINC_BOUND_MAX/XGAMINC_BOUND_MIN)/FLOAT(NGAMINC-1))
 !
 ALLOCATE( XGAMINC_RIM1(NGAMINC) )
 ALLOCATE( XGAMINC_RIM2(NGAMINC) )
@@ -732,13 +773,13 @@ XHMLINTP2 = 1.0 + XHMLINTP1*LOG( 25.E-6/(XGAMINC_HMC_BOUND_MIN)**(1.0/XALPHAC) )
 !
 !*       7.2    Constants for the accretion of raindrops onto aggregates
 !
-XFRACCSS = ((XPI**2)/24.0)*XCCS*XRHOLW*(ZRHO00**XCEXVT)
+XFRACCSS = XLBS*((XPI**2)/24.0)*XCCR*XRHOLW*(ZRHO00**XCEXVT)
 !
 XLBRACCS1   =    MOMG(XALPHAS,XNUS,2.)*MOMG(XALPHAR,XNUR,3.)
 XLBRACCS2   = 2.*MOMG(XALPHAS,XNUS,1.)*MOMG(XALPHAR,XNUR,4.)
 XLBRACCS3   =                          MOMG(XALPHAR,XNUR,5.)
 !
-XFSACCRG = (XPI/4.0)*XAS*XCCS*(ZRHO00**XCEXVT)
+XFSACCRG = XLBS*(XPI/4.0)*XAS*XCCR*(ZRHO00**XCEXVT)
 !
 XLBSACCR1   =    MOMG(XALPHAR,XNUR,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSACCR2   = 2.*MOMG(XALPHAR,XNUR,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
@@ -750,9 +791,9 @@ XLBSACCR3   =                          MOMG(XALPHAS,XNUS,XBS+2.)
 ! Notice: One magnitude of lambda discretized over 10 points for snow
 !
 NACCLBDAS = 40
-XACCLBDAS_MIN = 5.0E1 ! Minimal value of Lbda_s to tabulate XKER_RACCS
-XACCLBDAS_MAX = 5.0E5 ! Maximal value of Lbda_s to tabulate XKER_RACCS
-ZRATE = LOG(XACCLBDAS_MAX/XACCLBDAS_MIN)/REAL(NACCLBDAS-1)
+XACCLBDAS_MIN = 5.0E2*XTRANS_MP_GAMMAS !5.0E1*XTRANS_MP_GAMMAS ! Minimal value of Lbda_s to tabulate XKER_RACCS
+XACCLBDAS_MAX = 1.0E5*XTRANS_MP_GAMMAS !5.0E5*XTRANS_MP_GAMMAS ! Maximal value of Lbda_s to tabulate XKER_RACCS
+ZRATE = LOG(XACCLBDAS_MAX/XACCLBDAS_MIN)/FLOAT(NACCLBDAS-1)
 XACCINTP1S = 1.0 / ZRATE
 XACCINTP2S = 1.0 - LOG( XACCLBDAS_MIN ) / ZRATE
 NACCLBDAR = 40
@@ -785,15 +826,15 @@ IF( (KACCLBDAS/=NACCLBDAS) .OR. (KACCLBDAR/=NACCLBDAR) .OR. (KND/=IND) .OR. &
     (PACCLBDAS_MIN/=XACCLBDAS_MIN) .OR. (PACCLBDAR_MIN/=XACCLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, XFVELOS, XCR, XDR,                     &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
   CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, XFVELOS, XCR, XDR, 0.,                 &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCS                                        )
   CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBS, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBS, XCS, XDS, XFVELOS, XCR, XDR,                     &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,XAG, XBS, XAS                         )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -936,7 +977,7 @@ XCOLSG   = 0.01 ! Collection efficiency of S+G
 XCOLEXSG = 0.1  ! Temperature factor of the S+G collection efficiency
 WRITE (ILUOUT0, FMT=*) ' NEW Constants for the aggregate collection by the graupeln'
 WRITE (ILUOUT0, FMT=*) ' XCOLSG, XCOLEXSG  = ',XCOLSG,XCOLEXSG
-XFSDRYG = (XPI/4.0)*XCOLSG*XCCG*XCCS*XAS*(ZRHO00**XCEXVT)
+XFSDRYG = XLBS*(XPI/4.0)*XCOLSG*XCCG*XAS*(ZRHO00**XCEXVT)
 !
 XLBSDRYG1   =    MOMG(XALPHAG,XNUG,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSDRYG2   = 2.*MOMG(XALPHAG,XNUG,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
@@ -966,8 +1007,8 @@ ZRATE = LOG(XDRYLBDAR_MAX/XDRYLBDAR_MIN)/REAL(NDRYLBDAR-1)
 XDRYINTP1R = 1.0 / ZRATE
 XDRYINTP2R = 1.0 - LOG( XDRYLBDAR_MIN ) / ZRATE
 NDRYLBDAS = 80
-XDRYLBDAS_MIN = 2.5E1 ! Minimal value of Lbda_s to tabulate XKER_SDRYG
-XDRYLBDAS_MAX = 2.5E9 ! Maximal value of Lbda_s to tabulate XKER_SDRYG
+XDRYLBDAS_MIN = 5.0E2*XTRANS_MP_GAMMAS ! Minimal value of Lbda_s to tabulate XKER_SDRYG
+XDRYLBDAS_MAX = 1.0E5*XTRANS_MP_GAMMAS ! Maximal value of Lbda_s to tabulate XKER_SDRYG
 ZRATE = LOG(XDRYLBDAS_MAX/XDRYLBDAS_MIN)/REAL(NDRYLBDAS-1)
 XDRYINTP1S = 1.0 / ZRATE
 XDRYINTP2S = 1.0 - LOG( XDRYLBDAS_MIN ) / ZRATE
@@ -999,7 +1040,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAS/=NDRYLBDAS) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAS_MIN/=XDRYLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, XBS, XCG, XDG, XCS, XDS,                              &
+                ZEGS, XBS, XCG, XDG, 0., XCS, XDS, XFVELOS,                 &
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1065,7 +1106,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAR/=NDRYLBDAR) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAR_MIN/=XDRYLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAR, XNUR,                          &
-                ZEGR, XBR, XCG, XDG, XCR, XDR,                              &
+                ZEGR, XBR, XCG, XDG,0., XCR, XDR, 0.,                             &
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1126,7 +1167,8 @@ XFWETH = (XPI/4.0)*XCCH*XCH*(ZRHO00**XCEXVT)*MOMG(XALPHAH,XNUH,XDH+2.0)
 !
 !*       9.2.2  Constants for the aggregate collection by the hailstones
 !
-XFSWETH = (XPI/4.0)*XCCH*XCCS*XAS*(ZRHO00**XCEXVT)
+!XFSWETH = (XPI/4.0)*XCCH*XCCS*XAS*(ZRHO00**XCEXVT)
+XFSWETH = XLBS*(XPI/4.0)*XCCH*XAS*(ZRHO00**XCEXVT) ! Wurtz  
 !
 XLBSWETH1   =    MOMG(XALPHAH,XNUH,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSWETH2   = 2.*MOMG(XALPHAH,XNUH,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
@@ -1143,8 +1185,8 @@ XLBGWETH3   =                          MOMG(XALPHAG,XNUG,XBG+2.)
 ! Notice: One magnitude of lambda discretized over 10 points
 !
 NWETLBDAS = 80
-XWETLBDAS_MIN = 2.5E1 ! Minimal value of Lbda_s to tabulate XKER_SWETH
-XWETLBDAS_MAX = 2.5E9 ! Maximal value of Lbda_s to tabulate XKER_SWETH
+XWETLBDAS_MIN = 5.0E2*XTRANS_MP_GAMMAS ! Minimal value of Lbda_s to tabulate XKER_SWETH
+XWETLBDAS_MAX = 1.0E5*XTRANS_MP_GAMMAS ! Maximal value of Lbda_s to tabulate XKER_SWETH
 ZRATE = LOG(XWETLBDAS_MAX/XWETLBDAS_MIN)/REAL(NWETLBDAS-1)
 XWETINTP1S = 1.0 / ZRATE
 XWETINTP2S = 1.0 - LOG( XWETLBDAS_MIN ) / ZRATE
@@ -1182,7 +1224,7 @@ IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAS/=NWETLBDAS) .OR. (KND/=IND) .OR. &
     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAS_MIN/=XWETLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAS, XNUS,                          &
-                ZEHS, XBS, XCH, XDH, XCS, XDS,                              &
+                ZEHS, XBS, XCH, XDH,0., XCS, XDS, XFVELOS,                  &
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1248,7 +1290,7 @@ IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAG/=NWETLBDAG) .OR. (KND/=IND) .OR. &
     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAG_MIN/=XWETLBDAG_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAG, XNUG,                          &
-                ZEHG, XBG, XCH, XDH, XCG, XDG,                              &
+                ZEHG, XBG, XCH, XDH,0., XCG, XDG,   0.,                     &
                 XWETLBDAH_MAX, XWETLBDAG_MAX, XWETLBDAH_MIN, XWETLBDAG_MIN, &
                 ZFDINFTY, XKER_GWETH                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
diff --git a/src/MNH/ini_param_elec.f90 b/src/MNH/ini_param_elec.f90
index bdbd3c6d90293a6d4fde0fd28e4d9187e4712642..02a8b1578cfe2608304c04ccd8d8ebff6fd9025d 100644
--- a/src/MNH/ini_param_elec.f90
+++ b/src/MNH/ini_param_elec.f90
@@ -85,6 +85,7 @@ END MODULE MODI_INI_PARAM_ELEC
 !!        J. Escobar 8/01/2016 bug , missing YDIR='XY' in READ 
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -846,17 +847,17 @@ XLBQSACCRG3 =      MOMG(XALPHAS,XNUS,XFS)    * MOMG(XALPHAR,XNUR,2.)
 ZESR = 1.0
 !
 CALL RRCOLSS (KND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-              ZESR, XFR, XCS, XDS, XCR, XDR,                              &
+              ZESR, XFR, XCS, XDS, 0., XCR, XDR,                          &
               XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
               PFDINFTY, XKER_Q_RACCSS, XAG, XBS, XAS                      )
 !
 CALL RZCOLX  (KND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-              ZESR, XFR, XCS, XDS, XCR, XDR,                              &
+              ZESR, XFR, XCS, XDS, 0., XCR, XDR, 0.,                      &
               XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
               PFDINFTY, XKER_Q_RACCS                                      )
 !
 CALL RSCOLRG (KND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-              ZESR, XFS, XCS, XDS, XCR, XDR,                              &
+              ZESR, XFS, XCS, XDS, 0., XCR, XDR,                          &
               XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
               PFDINFTY, XKER_Q_SACCRG, XAG, XBS, XAS                      )
 !
@@ -878,7 +879,7 @@ XLBQSDRYG3 =      MOMG(XALPHAS,XNUS,XFS)    * MOMG(XALPHAG,XNUG,2.)
 ZEGS = 1. ! also initialized in ini_rain_ice_elec
 !
 CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-             ZEGS, XFS, XCG, XDG, XCS, XDS,                              &
+             ZEGS, XFS, XCG, XDG, 0., XCS, XDS, 0.,                      &
              XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
              PFDINFTY, XKER_Q_SDRYG                                      )
 !
@@ -897,7 +898,7 @@ IF (CNI_CHARGING == 'HELFA') THEN
 !
   IF( .NOT.ALLOCATED(XKER_Q_SDRYGB)) ALLOCATE( XKER_Q_SDRYGB(NDRYLBDAG,NDRYLBDAS) )
   CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-               ZEGS, 0., XCG, XDG, XCS, XDS,                               &
+               ZEGS, 0., XCG, XDG, 0., XCS, XDS, 0.,                       &
                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                PFDINFTY, XKER_Q_SDRYGB                                     )
 ! Delta vqb1_sg
@@ -999,7 +1000,7 @@ IF (CNI_CHARGING == 'TAKAH' .OR. CNI_CHARGING == 'SAP98' .OR. &
   XAUX_LIM3 =      MOMG(XALPHAG,XNUG,2.)
   IF( .NOT.ALLOCATED(XKER_Q_LIMSG)) ALLOCATE( XKER_Q_LIMSG(NDRYLBDAG,NDRYLBDAS) )
   CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-               ZEGS, 0., XCG, XDG, XCS, XDS,                               &
+               ZEGS, 0., XCG, XDG, 0., XCS, XDS, 0.,                       &
                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                PFDINFTY, XKER_Q_LIMSG)
 ENDIF
@@ -1021,7 +1022,7 @@ XLBQRDRYG3 =      MOMG(XALPHAR,XNUR,XFR)    * MOMG(XALPHAG,XNUG,2.)
 ZEGR = 1.0 
 !
 CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAR, XNUR,                            & 
-             ZEGR, XFR, XCG, XDG, XCR, XDR,                                &
+             ZEGR, XFR, XCG, XDG, 0., XCR, XDR, 0.,                        &
              XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN,   &
              PFDINFTY, XKER_Q_RDRYG                                        )
 !
diff --git a/src/MNH/ini_rain_ice.f90 b/src/MNH/ini_rain_ice.f90
index 62cabad5b587f48a6cd6d443c088c8c7cea8c2ca..b89279296f24871305d904846dc9eab06e79b265 100644
--- a/src/MNH/ini_rain_ice.f90
+++ b/src/MNH/ini_rain_ice.f90
@@ -102,6 +102,7 @@ END MODULE MODI_INI_RAIN_ICE
 !!      S. Riette 2016-11: new ICE3/ICE4 options
 !!      P. Wautelet 22/01/2019 bug correction: incorrect write
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -742,15 +743,15 @@ IF( (KACCLBDAS/=NACCLBDAS) .OR. (KACCLBDAR/=NACCLBDAR) .OR. (KND/=IND) .OR. &
     (PACCLBDAS_MIN/=XACCLBDAS_MIN) .OR. (PACCLBDAR_MIN/=XACCLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR,                          &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
   CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR, 0.,                      &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCS                                        )
   CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBS, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBS, XCS, XDS, 0., XCR, XDR,                          &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,  XAG, XBS, XAS                       )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -950,7 +951,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAS/=NDRYLBDAS) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAS_MIN/=XDRYLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, XBS, XCG, XDG, XCS, XDS,                              &
+                ZEGS, XBS, XCG, XDG, 0., XCS, XDS, 0.,                      &
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1016,7 +1017,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAR/=NDRYLBDAR) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAR_MIN/=XDRYLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAR, XNUR,                          &
-                ZEGR, XBR, XCG, XDG, XCR, XDR,                              &
+                ZEGR, XBR, XCG, XDG, 0., XCR, XDR, 0.,                      &
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1154,7 +1155,7 @@ IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAS/=NWETLBDAS) .OR. (KND/=IND) .OR. &
     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAS_MIN/=XWETLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAS, XNUS,                          &
-                ZEHS, XBS, XCH, XDH, XCS, XDS,                              &
+                ZEHS, XBS, XCH, XDH, 0., XCS, XDS, 0.,                      &
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1220,7 +1221,7 @@ IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAG/=NWETLBDAG) .OR. (KND/=IND) .OR. &
     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAG_MIN/=XWETLBDAG_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAG, XNUG,                          &
-                ZEHG, XBG, XCH, XDH, XCG, XDG,                              &
+                ZEHG, XBG, XCH, XDH, 0., XCG, XDG, 0.,                      &
                 XWETLBDAH_MAX, XWETLBDAG_MAX, XWETLBDAH_MIN, XWETLBDAG_MIN, &
                 ZFDINFTY, XKER_GWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1286,7 +1287,7 @@ IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAR/=NWETLBDAR) .OR. (KND/=IND) .OR. &
     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAR_MIN/=XWETLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAR, XNUR,                          &
-                ZEHR, XBR, XCH, XDH, XCR, XDR,                              &
+                ZEHR, XBR, XCH, XDH, 0., XCR, XDR, 0.,                      &
                 XWETLBDAH_MAX, XWETLBDAR_MAX, XWETLBDAH_MIN, XWETLBDAR_MIN, &
                 ZFDINFTY, XKER_RWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
diff --git a/src/MNH/ini_rain_ice_elec.f90 b/src/MNH/ini_rain_ice_elec.f90
index 940caeaeefc96dcda0800b7678b8be775815c116..d6ba287f181b5d1c342e8a64d68d02f00a474400 100644
--- a/src/MNH/ini_rain_ice_elec.f90
+++ b/src/MNH/ini_rain_ice_elec.f90
@@ -87,6 +87,7 @@ END MODULE MODI_INI_RAIN_ICE_ELEC
 !!    Modifications:
 !!      C. Barthe   20/11/09   update to version 4.8.1
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -694,15 +695,15 @@ IF( (KACCLBDAS/=NACCLBDAS) .OR. (KACCLBDAR/=NACCLBDAR) .OR. (KND/=IND) .OR. &
     (PACCLBDAS_MIN/=XACCLBDAS_MIN) .OR. (PACCLBDAR_MIN/=XACCLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR,                          &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
   CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR, 0.,                      &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCS                                        )
   CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBS, XCS, XDS, XCR, XDR,                              &
+                 ZESR, XBS, XCS, XDS, 0., XCR, XDR,                          &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,  XAG, XBS, XAS                       )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -905,7 +906,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAS/=NDRYLBDAS) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAS_MIN/=XDRYLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, XBS, XCG, XDG, XCS, XDS,                              &
+                ZEGS, XBS, XCG, XDG, 0., XCS, XDS, 0.,                      &
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -971,7 +972,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAR/=NDRYLBDAR) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAR_MIN/=XDRYLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAR, XNUR,                          &
-                ZEGR, XBR, XCG, XDG, XCR, XDR,                              &
+                ZEGR, XBR, XCG, XDG, 0., XCR, XDR, 0.,                      &
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1089,7 +1090,7 @@ IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAS/=NWETLBDAS) .OR. (KND/=IND) .OR. &
     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAS_MIN/=XWETLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAS, XNUS,                          &
-                ZEHS, XBS, XCH, XDH, XCS, XDS,                              &
+                ZEHS, XBS, XCH, XDH, 0., XCS, XDS, 0.,                      &
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1155,7 +1156,7 @@ IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAG/=NWETLBDAG) .OR. (KND/=IND) .OR. &
     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAG_MIN/=XWETLBDAG_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAG, XNUG,                          &
-                ZEHG, XBG, XCH, XDH, XCG, XDG,                              &
+                ZEHG, XBG, XCH, XDH, 0., XCG, XDG, 0.,                      &
                 XWETLBDAH_MAX, XWETLBDAG_MAX, XWETLBDAH_MIN, XWETLBDAG_MIN, &
                 ZFDINFTY, XKER_GWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
diff --git a/src/MNH/init_aerosol_concentration.f90 b/src/MNH/init_aerosol_concentration.f90
index e86998c4b18e3a4f712a57c016ea1cbb7e9c14d0..32494739c32a42c280773d7510b518e518d2089b 100644
--- a/src/MNH/init_aerosol_concentration.f90
+++ b/src/MNH/init_aerosol_concentration.f90
@@ -54,7 +54,7 @@ END MODULE MODI_INIT_AEROSOL_CONCENTRATION
 USE MODD_NSV
 USE MODD_PARAM_n,    ONLY : CCLOUD
 USE MODD_PARAM_LIMA, ONLY : LWARM, LACTI, NMOD_CCN, LSCAV, LAERO_MASS,      &
-                            XCCN_CONC, LCCN_HOM,                              &
+                            XCCN_CONC, LCCN_HOM,                            &
                             LCOLD, LNUCL, NMOD_IFN, LMEYERS,                &
                             XIFN_CONC, LIFN_HOM
 USE MODD_PARAMETERS, ONLY : JPVEXT
@@ -79,7 +79,7 @@ INTEGER                                 :: IKB, IKE
 !
 !
 IF ( LWARM .AND. LACTI ) THEN
-  DO JSV = NSV_LIMA_CCN_FREE, NSV_LIMA_CCN_ACTI+NMOD_CCN-1         
+  DO JSV = NSV_LIMA_CCN_FREE, NSV_LIMA_CCN_ACTI+NMOD_CCN-1
     PSVT(:,:,:,JSV) = 0.0
   ENDDO
   IKB = 1+JPVEXT
@@ -112,7 +112,7 @@ END IF ! LWARM AND LACTI
 ! Initialisation des concentrations en IFN
 !
 IF ( LCOLD .AND. LNUCL .AND. (.NOT. LMEYERS) ) THEN
-  DO JSV = NSV_LIMA_IFN_FREE, NSV_LIMA_IFN_NUCL+NMOD_IFN-1         
+  DO JSV = NSV_LIMA_IFN_FREE, NSV_LIMA_IFN_NUCL+NMOD_IFN-1
     PSVT(:,:,:,JSV) = 0.0
   ENDDO
   IKB = 1+JPVEXT
@@ -127,7 +127,7 @@ IF ( LCOLD .AND. LNUCL .AND. (.NOT. LMEYERS) ) THEN
   ELSE
 ! concentration décroissante selon z
     DO JSV = 1, NMOD_IFN
-      WHERE (PZZ(:,:,:) .LE. 1000.) 
+      WHERE (PZZ(:,:,:) .LE. 1000.)
          PSVT(:,:,:,NSV_LIMA_IFN_FREE+JSV-1) = XIFN_CONC(JSV)*1.0E3 / PRHODREF(:,:,:)
       ELSEWHERE (PZZ(:,:,:) .LE. 10000.)
          PSVT(:,:,:,NSV_LIMA_IFN_FREE+JSV-1) = XIFN_CONC(JSV)*1.0E3 &
diff --git a/src/MNH/lidar.f90 b/src/MNH/lidar.f90
index 4a3d987e8873486c7ce58e984cf78dc6b21a49da..1bee4acda2df0fdadc0b7db79d7b24dc5679c0eb 100644
--- a/src/MNH/lidar.f90
+++ b/src/MNH/lidar.f90
@@ -8,7 +8,7 @@
 !      #################
 !
 INTERFACE
-      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PCLDFR,PRT, &
+      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PT,PCLDFR,PRT, &
                        PLIDAROUT,PLIPAROUT,PCT,PDSTC,PDSTD,PDSTS)
 !
 CHARACTER(LEN=*),         INTENT(IN) :: HCLOUD  ! Name of the cloud scheme
@@ -17,6 +17,7 @@ REAL,                     INTENT(IN) :: PALT    ! Altitude of the lidar source
 REAL,                     INTENT(IN) :: PWVL    ! Wavelength of the lidar source
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PZZ     ! Altitude
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PRHO    ! Air density
+REAL, DIMENSION(:,:,:),   INTENT(IN) :: PT      ! Air temperature (C)
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PCLDFR  ! Cloud fraction
 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT     ! Moist variables at t
 REAL, DIMENSION(:,:,:),  INTENT(OUT) :: PLIDAROUT ! Lidar output
@@ -36,7 +37,7 @@ END INTERFACE
 !
 END MODULE MODI_LIDAR
 !     #########################################################
-      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PCLDFR,PRT, &
+      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PT,PCLDFR,PRT, &
                        PLIDAROUT,PLIPAROUT,PCT,PDSTC,PDSTD,PDSTS)
 !     #########################################################
 !
@@ -109,9 +110,11 @@ USE MODD_ICE_C1R3_DESCR,  ONLY : XLBEXI,                      &
 USE MODD_PARAM_LIMA,      ONLY : URTMIN=>XRTMIN, UCTMIN=>XCTMIN, &
                                  UALPHAC=>XALPHAC,UNUC=>XNUC, &
                                  UALPHAR=>XALPHAR,UNUR=>XNUR, &
-                                 UALPHAI=>XALPHAI,UNUI=>XNUI 
+                                 UALPHAI=>XALPHAI,UNUI=>XNUI, &
+                                 USNOW_T=>LSNOW_T
 USE MODD_PARAM_LIMA_COLD, ONLY : UCCS=>XCCS, UCXS=>XCXS, ULBEXS=>XLBEXS, & 
-                                                         ULBS=>XLBS
+                                                         ULBS=>XLBS,     &
+                                 XLBDAS_MAX,XLBDAS_MIN,  UBS=>XBS
 USE MODD_PARAM_LIMA_MIXED,ONLY : UCCG=>XCCG, UCXG=>XCXG, ULBEXG=>XLBEXG, &
                                                          ULBG=>XLBG
 
@@ -130,6 +133,7 @@ REAL,                     INTENT(IN) :: PALT    ! Altitude of the lidar source
 REAL,                     INTENT(IN) :: PWVL    ! Wavelength of the lidar source
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PZZ     ! Altitude
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PRHO    ! Air density
+REAL, DIMENSION(:,:,:),   INTENT(IN) :: PT      ! Air temperature (C)
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PCLDFR  ! Cloud fraction
 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT     ! Moist variables at t
 REAL, DIMENSION(:,:,:),  INTENT(OUT) :: PLIDAROUT ! Lidar output
@@ -523,9 +527,18 @@ SELECT CASE ( HCLOUD )
 !
             YDSD = 'MONOD'
             ZIWC    = PRHO(JI,JJ,JK)*PRT(JI,JJ,JK,5)
-            ZLBDAS  = ZLBS*(ZIWC)**ZLBEXS
+            IF (HCLOUD=='LIMA' .AND. USNOW_T) THEN
+               IF (PT(JI,JJ,JK)>-10.) THEN
+                  ZLBDAS = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(PT(JI,JJ,JK)+273.15))),XLBDAS_MIN)
+               ELSE
+                  ZLBDAS = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(PT(JI,JJ,JK)+273.15))),XLBDAS_MIN)
+               END IF
+               ZCONC=ZLBS*ZIWC*ZLBDAS**UBS
+            ELSE
+               ZLBDAS  = ZLBS*(ZIWC)**ZLBEXS
+               ZCONC   = ZCCS*(ZLBDAS)**ZCXS
+            END IF
             IF (ZLBDAS .GT. 0) THEN
-              ZCONC   = ZCCS*(ZLBDAS)**ZCXS
               ZRADIUS = 0.5*(3.0/ZLBDAS) ! Assume Marshall-Palmer law for Reff
               IANGLE  = 11
               CALL BHMIE_WATER( ZWAVE_LENGTH, ZZREFIND_ICE, YDSD, ZCONC,      &
diff --git a/src/MNH/lima_cold_slow_processes.f90 b/src/MNH/lima_cold_slow_processes.f90
index 9fcacdd5a39b10fc469e99ffba04621f848456a7..a95f91d462b981b9dec4184de9d767cddeecf288 100644
--- a/src/MNH/lima_cold_slow_processes.f90
+++ b/src/MNH/lima_cold_slow_processes.f90
@@ -80,6 +80,7 @@ END MODULE MODI_LIMA_COLD_SLOW_PROCESSES
 !  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -94,8 +95,8 @@ USE MODD_CST,             ONLY: XP00, XRD, XRV, XMV, XMD, XCPD, XCPV,        &
                                 XCL, XCI, XTT, XLSTT, XALPI, XBETAI, XGAMI
 USE MODD_NSV,             ONLY: NSV_LIMA_NI
 USE MODD_PARAMETERS,      ONLY: JPHEXT, JPVEXT
-USE MODD_PARAM_LIMA,      ONLY: LSNOW, XRTMIN, XCTMIN, XALPHAI, XALPHAS,     &
-                                XNUI
+USE MODD_PARAM_LIMA,      ONLY: LSNOW, LSNOW_T, XRTMIN, XCTMIN,              &
+                                XALPHAI, XALPHAS, XNUI, XNUS
 USE MODD_PARAM_LIMA_COLD, ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XBI, XCXS, XCCS, &
                                 XLBDAS_MAX, XDSCNVI_LIM, XLBDASCNVI_MAX,     &
                                 XC0DEPSI, XC1DEPSI, XR0DEPSI, XR1DEPSI,      &
@@ -103,7 +104,8 @@ USE MODD_PARAM_LIMA_COLD, ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XBI, XCXS, XCCS, &
                                 XDICNVS_LIM, XLBDAICNVS_LIM,                 &
                                 XC0DEPIS, XC1DEPIS, XR0DEPIS, XR1DEPIS,      &
                                 XCOLEXIS, XAGGS_CLARGE1, XAGGS_CLARGE2,      &
-                                XAGGS_RLARGE1, XAGGS_RLARGE2
+                                XAGGS_RLARGE1, XAGGS_RLARGE2, XBS,           &
+                                XLBDAS_MIN,XFVELOS,XTRANS_MP_GAMMAS
 
 use mode_budget,          only: Budget_store_init, Budget_store_end
 use mode_tools,           only: Countjv
@@ -316,9 +318,19 @@ IF( IMICRO >= 1 ) THEN
       ZLBDAI(:) = ( XLBI*ZCIT(:) / ZRIT(:) )**XLBEXI
    END WHERE
    ZLBDAS(:)  = 1.E10
-   WHERE (ZRST(:)>XRTMIN(5) )
-      ZLBDAS(:) = XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS
-   END WHERE
+   IF (LSNOW_T) THEN
+      WHERE(ZZT(:)>263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
+      END WHERE
+      WHERE(ZZT(:)<=263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
+      END WHERE
+      ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
+   ELSE
+      WHERE (ZRST(:)>XRTMIN(5) )
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX,XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS),XLBDAS_MIN)
+      END WHERE
+   END IF
 !
    ZKA(:) = 2.38E-2 + 0.0071E-2 * ( ZZT(:) - XTT )          ! k_a
    ZDV(:) = 0.211E-4 * (ZZT(:)/XTT)**1.94 * (XP00/ZPRES(:)) ! D_v
@@ -342,16 +354,11 @@ IF( IMICRO >= 1 ) THEN
           call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'CNVI', pcis(:, :, :) * prhodj(:, :, :) )
       end if
 
-      WHERE ( ZRST(:)>XRTMIN(5) )
-         ZLBDAS(:)  = MIN( XLBDAS_MAX,                                           &
-                           XLBS*( ZRHODREF(:)*MAX( ZRST(:),XRTMIN(5) ) )**XLBEXS )
-      END WHERE
       ZZW(:) = 0.0
       WHERE ( ZLBDAS(:)<XLBDASCNVI_MAX .AND. (ZRST(:)>XRTMIN(5)) &
                                        .AND. (ZSSI(:)<0.0)       )
          ZZW(:) = (ZLBDAS(:)*XDSCNVI_LIM)**(XALPHAS)
-         ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XCCS*ZLBDAS(:)**XCXS)/ZRHODREF(:) * (ZZW(:)**XNUI) &
-                                                               * EXP(-ZZW(:))
+         ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XLBS*ZRST(:)*ZLBDAS(:)**XBS) * (ZZW(:)**XNUI) * EXP(-ZZW(:))
 !
          ZZW(:) = MIN( ( XR0DEPSI+XR1DEPSI*ZCJ(:) )*ZZX(:),ZRSS(:) )
          ZRIS(:) = ZRIS(:) + ZZW(:)
@@ -384,8 +391,11 @@ IF( IMICRO >= 1 ) THEN
 
       ZZW(:) = 0.0
       WHERE ( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>ZRTMIN(5)) )
-         ZZW(:) = ( ZSSI(:)/(ZRHODREF(:)*ZAI(:)) ) *                               &
-                  ( X0DEPS*ZLBDAS(:)**XEX0DEPS + X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS )
+         ZZW(:) = ( ZRST(:)*ZSSI(:)/(ZAI(:)) ) *                                            &
+              ( X0DEPS*ZLBDAS(:)**XEX0DEPS +                                                &
+              (X1DEPS*ZCJ(:)*(1+(XFVELOS/(2.*ZLBDAS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) * &
+              (ZLBDAS(:))**(XEX1DEPS+XBS)))
+
          ZZW(:) =    MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
                    - MIN( ZRSS(:),ABS(ZZW(:)) )*(0.5-SIGN(0.5,ZZW(:)))
          ZRSS(:) = ZRSS(:) + ZZW(:)
@@ -420,8 +430,6 @@ IF( IMICRO >= 1 ) THEN
          ZZW(:) = (ZLBDAI(:)*XDICNVS_LIM)**(XALPHAI)
          ZZX(:) = ( ZSSI(:)/ZAI(:) )*ZCIT(:) * (ZZW(:)**XNUI) *EXP(-ZZW(:))
 !
-! Correction BVIE
-!         ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:)/ZRHODREF(:) &
          ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:) &
                             ,ZRIS(:) ) + ZRTMIN(5), ZRTMIN(5) ) - ZRTMIN(5)
          ZRIS(:) = ZRIS(:) - ZZW(:)
@@ -458,7 +466,7 @@ IF( IMICRO >= 1 ) THEN
       WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>ZRTMIN(4)) &
                                                             .AND. (ZCIS(:)>ZCTMIN(4)) )
          ZZW1(:,3) = (ZLBDAI(:) / ZLBDAS(:))**3
-         ZZW1(:,1) = (ZCIT(:)*(XCCS*ZLBDAS(:)**XCXS)/ZRHODREF(:)*EXP( XCOLEXIS*(ZZT(:)-XTT) )) &
+         ZZW1(:,1) = (ZCIT(:)*(XLBS*ZRST(:)*ZLBDAS(:)**XBS)*EXP(XCOLEXIS*(ZZT(:)-XTT) )) &
                                                     / (ZLBDAI(:)**3)
          ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:,3)),ZCIS(:) )
          ZCIS(:) = ZCIS(:) - ZZW1(:,2)
diff --git a/src/MNH/lima_conversion_melting_snow.f90 b/src/MNH/lima_conversion_melting_snow.f90
index ff5a691461b7a5de36d44febb6c320ce909eee12..307db0255eed1136811b4b159c846dcbfa692c6a 100644
--- a/src/MNH/lima_conversion_melting_snow.f90
+++ b/src/MNH/lima_conversion_melting_snow.f90
@@ -55,15 +55,16 @@ END MODULE MODI_LIMA_CONVERSION_MELTING_SNOW
 !!    -------------
 !!      Original             15/03/2018 
 !!
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT, XMV, XMD, XLVTT, XCPV, XCL, XESTT, XRV
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XNUS, XALPHAS
 USE MODD_PARAM_LIMA_MIXED, ONLY : XFSCVMG
-USE MODD_PARAM_LIMA_COLD,  ONLY : X0DEPS, XEX0DEPS, X1DEPS, XEX1DEPS
+USE MODD_PARAM_LIMA_COLD,  ONLY : X0DEPS, XEX0DEPS, X1DEPS, XEX1DEPS, XBS, XFVELOS
 !
 IMPLICIT NONE
 !
@@ -106,9 +107,10 @@ WHERE( (PRST(:)>XRTMIN(5)) .AND. (PT(:)>XTT) .AND. LDCOMPUTE(:) )
 !
 ! compute RSMLT
 !
-   ZW(:)  = XFSCVMG*MAX( 0.0,( -ZW(:) *             &
-                          ( X0DEPS*           PLBDS(:)**XEX0DEPS +     &
-                            X1DEPS*PCJ(:)*PLBDS(:)**XEX1DEPS ) ))!-    &
+   ZW(:)  = XFSCVMG*MAX( 0.0,( -ZW(:) * PRHODREF(:) * PRST(:) *          &
+                               ( X0DEPS*PLBDS(:)**XEX0DEPS +             &
+                                 X1DEPS*PCJ(:)*PLBDS(:)**(XEX1DEPS+XBS)* &
+                                   (1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) ))
 ! On ne tient pas compte de la collection de pluie et gouttelettes par la neige si T>0 !!!! 
 ! Note that no heat is exchanged because the graupeln produced are still icy!!!
    P_RS_CMEL(:) = - ZW(:)
diff --git a/src/MNH/lima_droplets_riming_snow.f90 b/src/MNH/lima_droplets_riming_snow.f90
index 6bef29df3bfac250078b40d9c8f45d2d00cc4dfa..b1c4a8007bddd182015e9aac542bec4c3e0d1c73 100644
--- a/src/MNH/lima_droplets_riming_snow.f90
+++ b/src/MNH/lima_droplets_riming_snow.f90
@@ -66,6 +66,7 @@ END MODULE MODI_LIMA_DROPLETS_RIMING_SNOW
 !!    -------------
 !!      Original             15/03/2018 
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -73,11 +74,11 @@ END MODULE MODI_LIMA_DROPLETS_RIMING_SNOW
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT, XNUS, XALPHAS
 USE MODD_PARAM_LIMA_MIXED, ONLY : NGAMINC, XRIMINTP1, XRIMINTP2, XGAMINC_RIM1, XGAMINC_RIM2, &
-                                  XCRIMSS, XEXCRIMSS, XSRIMCG, XEXSRIMCG, &
+                                  XCRIMSS, XEXCRIMSS, XSRIMCG, XEXSRIMCG, XSRIMCG2, XSRIMCG3, XEXSRIMCG2, &
                                   XHMLINTP1, XHMLINTP2, XGAMINC_HMC, XHM_FACTS, XHMTMIN, XHMTMAX
-USE MODD_PARAM_LIMA_COLD,  ONLY : XMNU0
+USE MODD_PARAM_LIMA_COLD,  ONLY : XMNU0, XBS, XFVELOS
 !
 IMPLICIT NONE
 !
@@ -171,7 +172,9 @@ WHERE( GRIM )
 !        4.     riming
 !
    ! Cloud droplets collected
-   P_RC_RIM(:) = - XCRIMSS * PRCT(:) * PLBDS(:)**XEXCRIMSS * PRHODREF(:)**(-XCEXVT)
+   P_RC_RIM(:) = - XCRIMSS  * PRCT(:) * PRST(:)*(1+(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
+                                      * PRHODREF(:)**(-XCEXVT+1) &
+                                      * (PLBDS(:)) ** (XEXCRIMSS+XBS)
    P_CC_RIM(:) = P_RC_RIM(:) *(PCCT(:)/PRCT(:)) ! Lambda_c**3
    !
    ! Cloud droplets collected on small aggregates add to snow
@@ -181,7 +184,7 @@ WHERE( GRIM )
    P_RG_RIM(:) = - P_RC_RIM(:) - P_RS_RIM(:) 
    !
    ! Large aggregates collecting droplets add to graupel (instant process ???)
-   ZZW3(:) = XSRIMCG * PLBDS(:)**XEXSRIMCG * (1.0 - ZZW2(:))/(PTSTEP*PRHODREF(:))
+   ZZW3(:) = PRST(:)*(1.0 - ZZW2(:))/PTSTEP
    P_RS_RIM(:) = P_RS_RIM(:) - ZZW3(:) 
    P_RG_RIM(:) = P_RG_RIM(:) + ZZW3(:) 
    !
diff --git a/src/MNH/lima_graupel.f90 b/src/MNH/lima_graupel.f90
index ad114da363f6c1616ce45ee6534f929d64845e2f..578f22bbf93fb04be12fb7e1c7edd9f1957cbe53 100644
--- a/src/MNH/lima_graupel.f90
+++ b/src/MNH/lima_graupel.f90
@@ -129,6 +129,7 @@ END MODULE MODI_LIMA_GRAUPEL
 !!    -------------
 !!      Original             15/03/2018
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -330,12 +331,12 @@ WHERE( GDRY )
        			                                    * (ZVEC1(:) - 1.0)
    ZZW(:) = ZVEC3(:)
 !
-   ZZW3(:) = XFSDRYG * ZZW(:) * EXP( XCOLEXSG*(PT(:)-XTT) )  & ! RSDRYG - rs collected by graupel in dry mode
-                    *( PLBDS(:)**(XCXS-XBS) )*( PLBDG(:)**XCXG )    &
-                    *( PRHODREF(:)**(-XCEXVT-1.) )                      &
-                    *( XLBSDRYG1/( PLBDG(:)**2                 ) + &
-                       XLBSDRYG2/( PLBDG(:)   * PLBDS(:)   ) + &
-                       XLBSDRYG3/(                  PLBDS(:)**2) )
+   ZZW3(:) = XFSDRYG * ZZW(:) * EXP( XCOLEXSG*(PT(:)-XTT) )       & ! RSDRYG - rs collected by graupel in dry mode
+                    *( PRST(:))*( PLBDG(:)**XCXG )                &
+                    *( PRHODREF(:)**(-XCEXVT) )                   &
+                    *( XLBSDRYG1/( PLBDG(:)**2                ) + &
+                       XLBSDRYG2/( PLBDG(:)    * PLBDS(:)     ) + &
+                       XLBSDRYG3/(               PLBDS(:)**2) )
 END WHERE
 !
 !*           1.c  Collection of rr in the dry mode
diff --git a/src/MNH/lima_ice_aggregation_snow.f90 b/src/MNH/lima_ice_aggregation_snow.f90
index 15e01ec84b33a508d8b30285ea944540185b1015..24f74dca1a89816d0ae2f2eced0cdd1080a16531 100644
--- a/src/MNH/lima_ice_aggregation_snow.f90
+++ b/src/MNH/lima_ice_aggregation_snow.f90
@@ -52,6 +52,7 @@ END MODULE MODI_LIMA_ICE_AGGREGATION_SNOW
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             15/03/2018
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 !
@@ -61,7 +62,7 @@ END MODULE MODI_LIMA_ICE_AGGREGATION_SNOW
 USE MODD_CST,             ONLY : XTT
 USE MODD_PARAM_LIMA,      ONLY : XRTMIN, XCTMIN
 USE MODD_PARAM_LIMA_COLD, ONLY : XBI, XCCS, XCXS, XCOLEXIS, XAGGS_CLARGE1, XAGGS_CLARGE2, &
-                                 XAGGS_RLARGE1, XAGGS_RLARGE2
+                                 XAGGS_RLARGE1, XAGGS_RLARGE2, XBS, XLBS
 !
 IMPLICIT NONE
 !
@@ -101,7 +102,7 @@ P_CI_AGGS(:) = 0.
 !
 WHERE ( (PRIT(:)>XRTMIN(4)) .AND. (PRST(:)>XRTMIN(5)) .AND. LDCOMPUTE(:) )
    ZZW1(:) = (PLBDI(:) / PLBDS(:))**3
-   ZZW2(:) = (PCIT(:)*(XCCS*PLBDS(:)**XCXS)/PRHODREF(:)*EXP( XCOLEXIS*(PT(:)-XTT) )) &
+   ZZW2(:) = (PCIT(:)*(XLBS*PRST(:)*PLBDS(:)**XBS)*EXP(XCOLEXIS*(PT(:)-XTT) )) &
         / (PLBDI(:)**3)
    ZZW3(:) = ZZW2(:)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:))
 !
diff --git a/src/MNH/lima_ice_snow_deposition.f90 b/src/MNH/lima_ice_snow_deposition.f90
index 4d92b528ac9aabb0224e61ae9de0c23a5b50f0fb..f31d3175d58108eecb631b515567ca27e669d202 100644
--- a/src/MNH/lima_ice_snow_deposition.f90
+++ b/src/MNH/lima_ice_snow_deposition.f90
@@ -78,21 +78,22 @@ SUBROUTINE LIMA_ICE_SNOW_DEPOSITION (PTSTEP, LDCOMPUTE,                        &
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             15/03/2018
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
 !              ------------
 !
-USE MODD_PARAM_LIMA,      ONLY : XRTMIN, XCTMIN, XALPHAI, XALPHAS, XNUI, XNUS 
-USE MODD_PARAM_LIMA_COLD, ONLY : XCXS, XCCS, &
+USE MODD_PARAM_LIMA,      ONLY : XRTMIN, XCTMIN, XALPHAI, XALPHAS, XNUI, XNUS
+USE MODD_PARAM_LIMA_COLD, ONLY : XCXS, XCCS, XLBS, XBS, &
                                  XLBDAS_MAX, XDSCNVI_LIM, XLBDASCNVI_MAX,     &
                                  XC0DEPSI, XC1DEPSI, XR0DEPSI, XR1DEPSI,      &
                                  XSCFAC, X1DEPS, X0DEPS, XEX1DEPS, XEX0DEPS,  &
                                  XDICNVS_LIM, XLBDAICNVS_LIM,                 &
                                  XC0DEPIS, XC1DEPIS, XR0DEPIS, XR1DEPIS,      &
                                  XCOLEXIS, XAGGS_CLARGE1, XAGGS_CLARGE2,      &
-                                 XAGGS_RLARGE1, XAGGS_RLARGE2  
+                                 XAGGS_RLARGE1, XAGGS_RLARGE2, XFVELOS
 
 !
 IMPLICIT NONE
@@ -166,7 +167,7 @@ WHERE( GMICRO )
    WHERE ( PLBDS(:)<XLBDASCNVI_MAX .AND. (PRST(:)>XRTMIN(5)) &
                                    .AND. (PSSI(:)<0.0)       )
       ZZW(:) = (PLBDS(:)*XDSCNVI_LIM)**(XALPHAS)
-      ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XCCS*PLBDS(:)**XCXS)/PRHODREF(:) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
+      ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XLBS*PRST(:)*PLBDS(:)**XBS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
 !
       ZZW(:) = ( XR0DEPSI+XR1DEPSI*PCJ(:) )*ZZX(:)
 !
@@ -187,8 +188,10 @@ WHERE( GMICRO )
 !
    ZZW(:) = 0.0
    WHERE ( (PRST(:)>XRTMIN(5)) )
-      ZZW(:) = ( PSSI(:)/(PAI(:))/PRHODREF(:) ) * &
-           ( X0DEPS*PLBDS(:)**XEX0DEPS + X1DEPS*PCJ(:)*PLBDS(:)**XEX1DEPS )
+      ZZW(:) =( PRST(:)*PSSI(:)/PAI(:) ) *                               &
+           ( X0DEPS*PLBDS(:)**XEX0DEPS +                &
+           ( X1DEPS*PCJ(:)*(PLBDS(:))**(XBS+XEX1DEPS) * &
+                    (1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)))
       ZZW(:) =    ZZW(:)*(0.5+SIGN(0.5,ZZW(:))) - ABS(ZZW(:))*(0.5-SIGN(0.5,ZZW(:)))
    END WHERE
 !
diff --git a/src/MNH/lima_mixed.f90 b/src/MNH/lima_mixed.f90
index 49024b7b518893f1b9b101175676fb8a5d18f558..4d558eb1da83e2eb5d500c0ee18dbce977d80c42 100644
--- a/src/MNH/lima_mixed.f90
+++ b/src/MNH/lima_mixed.f90
@@ -96,6 +96,7 @@ END MODULE MODI_LIMA_MIXED
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets (no more call to budget in this subroutine)
 !  P. Wautelet 28/05/2020: bugfix: correct array start for PSVT and PSVS
 !  P. Wautelet 02/02/2021: budgets: add missing source terms for SV budgets in LIMA
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -108,9 +109,9 @@ USE MODD_CST,              ONLY: XP00, XRD, XRV, XMV, XMD, XCPD, XCPV,       &
 USE MODD_NSV
 USE MODD_PARAMETERS,       ONLY: JPHEXT, JPVEXT
 USE MODD_PARAM_LIMA,       ONLY: NMOD_IFN, XRTMIN, XCTMIN, LWARM, LCOLD,     &
-                                 NMOD_CCN, NMOD_IMM, LRAIN, LSNOW, LHAIL
+                                 NMOD_CCN, NMOD_IMM, LRAIN, LSNOW, LHAIL, LSNOW_T
 USE MODD_PARAM_LIMA_WARM,  ONLY: XLBC, XLBEXC, XLBR, XLBEXR
-USE MODD_PARAM_LIMA_COLD,  ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XSCFAC
+USE MODD_PARAM_LIMA_COLD,  ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XSCFAC, XLBDAS_MAX, XLBDAS_MIN, XTRANS_MP_GAMMAS
 USE MODD_PARAM_LIMA_MIXED, ONLY: XLBG, XLBEXG, XLBH, XLBEXH
 
 use mode_tools,            only: Countjv
@@ -467,9 +468,19 @@ IF( IMICRO >= 1 ) THEN
       ZLBDAI(:) = ( XLBI*ZCIT(:) / ZRIT(:) )**XLBEXI
    END WHERE
    ZLBDAS(:)  = 1.E10
-   WHERE (ZRST(:)>XRTMIN(5) )
-      ZLBDAS(:) = XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS
-   END WHERE
+   IF (LSNOW_T) THEN
+      WHERE(ZZT(:)>263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
+      END WHERE
+      WHERE(ZZT(:)<=263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
+      END WHERE
+      ZLBDAS(:) = ZLBDAS(:)*XTRANS_MP_GAMMAS
+   ELSE
+      WHERE (ZRST(:)>XRTMIN(5))
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS),XLBDAS_MIN)
+      END WHERE
+   END IF
    ZLBDAG(:)  = 1.E10
    WHERE (ZRGT(:)>XRTMIN(6) )
       ZLBDAG(:) = XLBG*( ZRHODREF(:)*ZRGT(:) )**XLBEXG
diff --git a/src/MNH/lima_mixed_fast_processes.f90 b/src/MNH/lima_mixed_fast_processes.f90
index 09c86c8a20e23fb9cd16adb859604d1aae2c37e5..7cfeffec02fef10d912cece900b69ca2983494bc 100644
--- a/src/MNH/lima_mixed_fast_processes.f90
+++ b/src/MNH/lima_mixed_fast_processes.f90
@@ -142,6 +142,7 @@ END MODULE MODI_LIMA_MIXED_FAST_PROCESSES
 !!      C. Barthe  * LACy *  jan. 2014    add budgets
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -297,9 +298,11 @@ IF( IGRIM>0 ) THEN
 !
    WHERE ( GRIM(:) )
       ZZW1(:,1) = MIN( ZRCS(:),                              &
-  	           XCRIMSS * ZZW(:) * ZRCT(:)                & ! RCRIMSS
-  	                            *   ZLBDAS(:)**XEXCRIMSS &
-    			            * ZRHODREF(:)**(-XCEXVT) )
+             XCRIMSS * ZZW(:) * ZRCT(:)                &
+                                      *  ZRST(:)*(1+(XFVELOS/ZLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
+                                      * ZRHODREF(:)**(-XCEXVT+1) &
+                                      * (ZLBDAS(:)) ** (XEXCRIMSS+XBS) )
+!
       ZRCS(:) = ZRCS(:) - ZZW1(:,1)
       ZRSS(:) = ZRSS(:) + ZZW1(:,1)
       ZTHS(:) = ZTHS(:) + ZZW1(:,1)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RCRIMSS))
@@ -318,14 +321,15 @@ IF( IGRIM>0 ) THEN
 !
 !
    WHERE ( GRIM(:) .AND. (ZRSS(:)>XRTMIN(5)/PTSTEP) )
-      ZZW1(:,2) = MIN( ZRCS(:),                     &
-    	           XCRIMSG * ZRCT(:)                & ! RCRIMSG
-    	                   *  ZLBDAS(:)**XEXCRIMSG  &
-  	                   * ZRHODREF(:)**(-XCEXVT) &
-    		           - ZZW1(:,1)              )
+      ZZW1(:,2) = MIN( ZRCS(:),                     & 
+                   XCRIMSG * ZRCT(:)*  ZRST(:)                & ! RCRIMSG
+                           *(1+(XFVELOS/ZLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSG/XALPHAS)*ZLBDAS(:)**(XBS+XEXCRIMSG)  &
+                           * ZRHODREF(:)**(-XCEXVT+1) &
+                           - ZZW1(:,1)              )
+
       ZZW1(:,3) = MIN( ZRSS(:),                         &
-                       XSRIMCG * ZLBDAS(:)**XEXSRIMCG   & ! RSRIMCG
-   	                       * (1.0 - ZZW(:) )/(PTSTEP*ZRHODREF(:)))
+                       XSRIMCG * XLBS * ZRST(:) * (1.0 - ZZW(:))/PTSTEP )
+
       ZRCS(:) = ZRCS(:) - ZZW1(:,2)
       ZRSS(:) = ZRSS(:) - ZZW1(:,3)
       ZRGS(:) = ZRGS(:) + ZZW1(:,2) + ZZW1(:,3)
@@ -477,7 +481,7 @@ IF( IGACC>0 .AND. LRAIN) THEN
 !
    WHERE ( GACC(:) )
       ZZW1(:,2) = ZCRT(:) *                                           & !! coef of RRACCS
-              XFRACCSS*( ZLBDAS(:)**XCXS )*( ZRHODREF(:)**(-XCEXVT-1.) ) &
+              XFRACCSS*( ZRST(:)*ZLBDAS(:)**XBS )*( ZRHODREF(:)**(-XCEXVT) ) &
          *( XLBRACCS1/((ZLBDAS(:)**2)               ) +                  &
             XLBRACCS2/( ZLBDAS(:)    * ZLBDAR(:)    ) +                  &
             XLBRACCS3/(               (ZLBDAR(:)**2)) )/ZLBDAR(:)**3
@@ -521,7 +525,7 @@ IF( IGACC>0 .AND. LRAIN) THEN
    WHERE ( GACC(:) .AND. (ZRSS(:)>XRTMIN(5)/PTSTEP) )
       ZZW1(:,2) = MAX( MIN( ZRRS(:),ZZW1(:,2)-ZZW1(:,4) ) , 0. )      ! RRACCSG
       ZZW1(:,3) = MIN( ZRSS(:),XFSACCRG*ZZW(:)*                     & ! RSACCRG
-            ( ZLBDAS(:)**(XCXS-XBS) )*( ZRHODREF(:)**(-XCEXVT-1.) ) &
+            ( ZRST(:) )*( ZRHODREF(:)**(-XCEXVT) ) &
            *( XLBSACCR1/((ZLBDAR(:)**2)               ) +           &
               XLBSACCR2/( ZLBDAR(:)    * ZLBDAS(:)    ) +           &
               XLBSACCR3/(               (ZLBDAS(:)**2)) ) )
@@ -573,8 +577,9 @@ WHERE( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>XRTMIN(5)/PTSTEP) .AND. (ZZT(:)>XTT) )
 ! compute RSMLT
 !
    ZZW(:)  = MIN( ZRSS(:), XFSCVMG*MAX( 0.0,( -ZZW(:) *             &
-                          ( X0DEPS*       ZLBDAS(:)**XEX0DEPS +     &
-                            X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS ) -   &
+                          ZRHODREF(:) * ZRST(:)*( X0DEPS*            ZLBDAS(:)**XEX0DEPS +     &
+                          X1DEPS*ZCJ(:)*(1+(XFVELOS/(2.*ZLBDAS(:)))**XALPHAS) &
+                          **(-XNUS+XEX1DEPS/XALPHAS)*(ZLBDAS(:))**(XEX1DEPS+XBS))-    &
                                     ( ZZW1(:,1)+ZZW1(:,4) ) *       &
                              ( ZRHODREF(:)*XCL*(XTT-ZZT(:))) ) /    &
                                             ( ZRHODREF(:)*XLMTT ) ) )
@@ -738,14 +743,14 @@ IF( IGDRY>0 ) THEN
    END DO
    ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
 !
-   WHERE( GDRY(:) )
+   WHERE( GDRY(:) )  
       ZZW1(:,3) = MIN( ZRSS(:),XFSDRYG*ZZW(:)                         & ! RSDRYG
                                       * EXP( XCOLEXSG*(ZZT(:)-XTT) )  &
-                    *( ZLBDAS(:)**(XCXS-XBS) )*( ZLBDAG(:)**XCXG )    &
-                    *( ZRHODREF(:)**(-XCEXVT-1.) )                    &
+                    *( ZRST(:)) )*( ZLBDAG(:)**XCXG )    &
+                    *( ZRHODREF(:)**(-XCEXVT) )                    &
                          *( XLBSDRYG1/( ZLBDAG(:)**2              ) + &
                             XLBSDRYG2/( ZLBDAG(:)   * ZLBDAS(:)   ) + &
-                            XLBSDRYG3/(               ZLBDAS(:)**2) ) )
+                            XLBSDRYG3/(               ZLBDAS(:)**2) )
    END WHERE
    DEALLOCATE(IVEC2)
    DEALLOCATE(IVEC1)
@@ -1185,8 +1190,8 @@ IF( IHAIL>0 ) THEN
 !
       WHERE( GWET(:) )
         ZZW1(:,3) = MIN( ZRSS(:),XFSWETH*ZZW(:)                       & ! RSWETH
-                      *( ZLBDAS(:)**(XCXS-XBS) )*( ZLBDAH(:)**XCXH )  &
-       	                 *( ZRHODREF(:)**(-XCEXVT-1.) )               &
+                      *( ZRST(:))*( ZLBDAH(:)**XCXH )  &
+                        *( ZRHODREF(:)**(-XCEXVT) )               &
                          *( XLBSWETH1/( ZLBDAH(:)**2              ) + &
                             XLBSWETH2/( ZLBDAH(:)   * ZLBDAS(:)   ) + &
                             XLBSWETH3/(               ZLBDAS(:)**2) ) )
diff --git a/src/MNH/lima_rain_accr_snow.f90 b/src/MNH/lima_rain_accr_snow.f90
index 01c31afbe3ff0152065142f33475281538a3c6ac..7c5c9bee51e37bc3c6d453f5a03b6fc07f3fa379 100644
--- a/src/MNH/lima_rain_accr_snow.f90
+++ b/src/MNH/lima_rain_accr_snow.f90
@@ -60,6 +60,7 @@ END MODULE MODI_LIMA_RAIN_ACCR_SNOW
 !!    -------------
 !!      Original             15/03/2018
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -68,7 +69,7 @@ END MODULE MODI_LIMA_RAIN_ACCR_SNOW
 !
 USE MODD_CST,              ONLY : XTT
 USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT
-USE MODD_PARAM_LIMA_COLD,  ONLY : XBS, XCXS
+USE MODD_PARAM_LIMA_COLD,  ONLY : XBS, XCXS, XTRANS_MP_GAMMAS
 USE MODD_PARAM_LIMA_MIXED, ONLY : NACCLBDAS, XACCINTP1S, XACCINTP2S,         &
                                   NACCLBDAR, XACCINTP1R, XACCINTP2R,         &
                                   XKER_RACCSS, XKER_RACCS, XKER_SACCRG,      &
@@ -141,7 +142,7 @@ WHERE( GACC )
 !
 !        1.3.1  select the (ZLBDAS,ZLBDAR) couplet
    !
-   ZVEC1(:) = MAX(MIN(PLBDS(:),5.E5),5.E1)
+   ZVEC1(:) = MAX(MIN(PLBDS(:),5.E5*XTRANS_MP_GAMMAS),5.E1*XTRANS_MP_GAMMAS)
    ZVEC2(:) = PLBDR(:)
 !
 !        1.3.2  find the next lower indice for the ZLBDAS and for the ZLBDAR
@@ -212,42 +213,28 @@ WHERE( GACC )
 !      
 ! BVIE manque PCRT ???????????????????????????????????
 !      ZZW4(:) =                                            & !! coef of RRACCS and RRACCS
-   ZZW4(:) = PCRT(:)                                     & !! coef of RRACCS and RRACCS
-         *  XFRACCSS *( PLBDS(:)**XCXS )*( PRHODREF(:)**(-XCEXVT-1.) ) &
-         *( XLBRACCS1/( PLBDS(:)**2 )                +                     &
-            XLBRACCS2/( PLBDS(:) * PLBDR(:)    ) +                     &
-            XLBRACCS3/(                PLBDR(:)**2 ) ) / PLBDR(:)**3
+   ZZW4(:) = PCRT(:)                                                       & !! coef of RRACCS and RRACCS
+         *  XFRACCSS *( PRST(:)*PLBDS(:)**XBS )*( PRHODREF(:)**(-XCEXVT) ) &
+         *( XLBRACCS1/( PLBDS(:)**2               ) +                      &
+            XLBRACCS2/( PLBDS(:)    * PLBDR(:)    ) +                      &
+            XLBRACCS3/(               PLBDR(:)**2 ) ) / PLBDR(:)**3
 
-!      ZRRS(:) = ZRRS(:) - ZZW1(:,4)
-!      ZRSS(:) = ZRSS(:) + ZZW1(:,4)
-!      ZTHS(:) = ZTHS(:) + ZZW1(:,4)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RRACCSS))
-!
-!      ZCRS(:) = MAX( ZCRS(:)-ZZW1(:,4)*(ZCRT(:)/ZRRT(:)),0.0 ) ! Lambda_r**3 
 !
 !        1.3.6  raindrop accretion-conversion of the large sized aggregates
 !               into graupeln
 !
-   ZZW5(:) = XFSACCRG*ZZW3(:) *                             & ! RSACCRG
-            ( PLBDS(:)**(XCXS-XBS) )*( PRHODREF(:)**(-XCEXVT-1.) ) &
-           *( XLBSACCR1/((PLBDR(:)**2)               ) +           &
-              XLBSACCR2/( PLBDR(:)    * PLBDS(:) ) +           &
-              XLBSACCR3/(                  (PLBDS(:)**2)) )
-      !
-!      P_RR_ACC(:) = - ZZW4(:) * ZZW1(:)           ! RRACCSS
-!      P_CR_ACC(:) = P_RR_ACC(:) * PCRT(:)/PRRT(:) ! Lambda_r**3 
-!      P_RS_ACC(:) = - P_RR_ACC(:) 
-      !
-!      P_RR_ACC(:) = P_RR_ACC(:) - ( ZZW2(:)-P_RS_ACC(:) )
-!      P_CR_ACC(:) = P_CR_ACC(:) - ( ZZW2(:)-P_RS_ACC(:) ) * PCRT(:)/PRRT(:) ! Lambda_r**3
-!      P_RS_ACC(:) = P_RS_ACC(:) - ZZW5(:)
-!      P_RG_ACC(:) = ( ZZW2(:)-P_RS_ACC(:) ) + ZZW5(:)
-      !
+   ZZW5(:) = XFSACCRG*ZZW3(:) *                         & ! RSACCRG
+            ( PRST(:) )*( PRHODREF(:)**(-XCEXVT)    )   &
+           *( XLBSACCR1/( PLBDR(:)**2               ) + &
+              XLBSACCR2/( PLBDR(:)    * PLBDS(:)    ) + &
+              XLBSACCR3/(               PLBDS(:)**2 ) )
+!
    P_RR_ACC(:) = - ZZW4(:) *  ZZW2(:)
    P_CR_ACC(:) = P_RR_ACC(:) * PCRT(:)/PRRT(:)
    P_RS_ACC(:) = ZZW4(:) *  ZZW1(:) - ZZW5(:)
    P_RG_ACC(:) = ZZW4(:) * ( ZZW2(:) - ZZW1(:) ) + ZZW5(:)
    P_TH_ACC(:) = - P_RR_ACC(:) * (PLSFACT(:)-PLVFACT(:))
-   !
+!
 END WHERE
 !
 !
diff --git a/src/MNH/lima_sedimentation.f90 b/src/MNH/lima_sedimentation.f90
index 365ae0f23362e17a84e8f9ab1682d8dc165f38dd..8d48b776d8a46279ab4b6fd8268068faf0411f5e 100644
--- a/src/MNH/lima_sedimentation.f90
+++ b/src/MNH/lima_sedimentation.f90
@@ -66,6 +66,7 @@ END MODULE MODI_LIMA_SEDIMENTATION
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
 !  B. Vie         03/2020: disable temperature change of droplets by air temperature
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -73,10 +74,11 @@ END MODULE MODI_LIMA_SEDIMENTATION
 !
 USE MODD_CST,              ONLY: XRHOLW, XCL, XCI
 USE MODD_PARAMETERS,       ONLY: JPHEXT, JPVEXT
-USE MODD_PARAM_LIMA,       ONLY: XCEXVT, XRTMIN, XCTMIN, NSPLITSED, &
-                                 XLB, XLBEX, XD, XFSEDR, XFSEDC,    &
-                                 XALPHAC, XNUC
-USE MODD_PARAM_LIMA_COLD,  ONLY: XLBEXI, XLBI, XDI
+USE MODD_PARAM_LIMA,       ONLY: XCEXVT, XRTMIN, XCTMIN, NSPLITSED,           &
+                                 XLB, XLBEX, XD, XFSEDR, XFSEDC,              &
+                                 XALPHAC, XNUC, XALPHAS, XNUS, LSNOW_T
+USE MODD_PARAM_LIMA_COLD,  ONLY: XLBEXI, XLBI, XDI, XLBDAS_MAX, XBS, XEXSEDS, &
+                                 XLBDAS_MIN, XTRANS_MP_GAMMAS, XFVELOS
 
 use mode_tools,            only: Countjv
 
@@ -183,11 +185,24 @@ DO JN = 1 ,  NSPLITSED(KID)
          IF (KMOMENTS==2) ZCS(JL) = PCS(I1(JL),I2(JL),I3(JL))
       END DO
 !
-      IF (KMOMENTS==1) ZLBDA(:) = XLB(KID) * ( ZRHODREF(:) * ZRS(:) )**XLBEX(KID)
-      IF (KMOMENTS==2) ZLBDA(:) = ( XLB(KID)*ZCS(:) / ZRS(:) )**XLBEX(KID)
+      IF (KID == 5 .AND. LSNOW_T) THEN
+         ZLBDA(:) = 1.E10
+         WHERE(ZT(:)>263.15 .AND. ZRS(:)>XRTMIN(5))
+            ZLBDA(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(:))),XLBDAS_MIN)
+         END WHERE
+         WHERE(ZT(:)<=263.15 .AND. ZRS(:)>XRTMIN(5))
+            ZLBDA(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(:))),XLBDAS_MIN)
+         END WHERE
+         ZLBDA(:) = ZLBDA(:)*XTRANS_MP_GAMMAS
+         ZZW(:) = XFSEDR(KID) * ZRHODREF(:)**(1.-XCEXVT)*ZRS(:)* &
+              (1 + (XFVELOS/ZLBDA(:))**XALPHAS)**(-XNUS+XEXSEDS/XALPHAS) * ZLBDA(:)**(XBS+XEXSEDS)
+      ELSE
+         IF (KMOMENTS==1) ZLBDA(:) = XLB(KID) * ( ZRHODREF(:) * ZRS(:) )**XLBEX(KID)
+         IF (KMOMENTS==2) ZLBDA(:) = ( XLB(KID)*ZCS(:) / ZRS(:) )**XLBEX(KID)
+         ZZY(:) = ZRHODREF(:)**(-XCEXVT) * ZLBDA(:)**(-XD(KID))
+         ZZW(:) = XFSEDR(KID) * ZRS(:) * ZZY(:) * ZRHODREF(:)
+      END IF ! Wurtz
 !
-      ZZY(:) = ZRHODREF(:)**(-XCEXVT) * ZLBDA(:)**(-XD(KID))
-      ZZW(:) = XFSEDR(KID) * ZRS(:) * ZZY(:) * ZRHODREF(:)
       IF (KMOMENTS==2) ZZX(:) = XFSEDC(KID) * ZCS(:) * ZZY(:) * ZRHODREF(:)
 
       IF (KID==2) THEN
diff --git a/src/MNH/lima_snow_deposition.f90 b/src/MNH/lima_snow_deposition.f90
index d39f714d83639730bbff492e14ae8001b06cdd64..1f7e1dde4d4dee5485ecd441ea0b960cf1bf14ba 100644
--- a/src/MNH/lima_snow_deposition.f90
+++ b/src/MNH/lima_snow_deposition.f90
@@ -63,20 +63,21 @@ SUBROUTINE LIMA_SNOW_DEPOSITION (LDCOMPUTE,                                &
 !!    -------------
 !!      Original             15/03/2018
 !!
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
 !              ------------
 !
-USE MODD_PARAM_LIMA,      ONLY : XRTMIN, XCTMIN, XALPHAI, XALPHAS, XNUI, XNUS 
-USE MODD_PARAM_LIMA_COLD, ONLY : XCXS, XCCS, &
+USE MODD_PARAM_LIMA,      ONLY : XRTMIN, XCTMIN, XALPHAI, XALPHAS, XNUI, XNUS
+USE MODD_PARAM_LIMA_COLD, ONLY : XCXS, XCCS,XLBS,XBS, &
                                  XLBDAS_MAX, XDSCNVI_LIM, XLBDASCNVI_MAX,     &
                                  XC0DEPSI, XC1DEPSI, XR0DEPSI, XR1DEPSI,      &
                                  XSCFAC, X1DEPS, X0DEPS, XEX1DEPS, XEX0DEPS,  &
                                  XDICNVS_LIM, XLBDAICNVS_LIM,                 &
                                  XC0DEPIS, XC1DEPIS, XR0DEPIS, XR1DEPIS,      &
                                  XCOLEXIS, XAGGS_CLARGE1, XAGGS_CLARGE2,      &
-                                 XAGGS_RLARGE1, XAGGS_RLARGE2  
+                                 XAGGS_RLARGE1, XAGGS_RLARGE2, XFVELOS
 
 !
 IMPLICIT NONE
@@ -132,7 +133,7 @@ WHERE( GMICRO )
    WHERE ( PLBDS(:)<XLBDASCNVI_MAX .AND. (PRST(:)>XRTMIN(5)) &
                                    .AND. (PSSI(:)<0.0)       )
       ZZW(:) = (PLBDS(:)*XDSCNVI_LIM)**(XALPHAS)
-      ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XCCS*PLBDS(:)**XCXS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
+      ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XLBS*PRST(:)*PLBDS(:)**XBS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
 !
       ZZW(:) = ( XR0DEPSI+XR1DEPSI*PCJ(:) )*ZZX(:)
 !
@@ -149,8 +150,10 @@ WHERE( GMICRO )
 !
    ZZW(:) = 0.0
    WHERE ( (PRST(:)>XRTMIN(5)) )
-      ZZW(:) = ( PSSI(:)/(PAI(:)*PRHODREF(:)) ) * &
-           ( X0DEPS*PLBDS(:)**XEX0DEPS + X1DEPS*PCJ(:)*PLBDS(:)**XEX1DEPS )
+      ZZW(:) = ( PRST(:)*PSSI(:)/(PAI(:)) ) *           &
+           ( X0DEPS*PLBDS(:)**XEX0DEPS +                &
+           ( X1DEPS*PCJ(:)*(PLBDS(:))**(XBS+XEX1DEPS) * &
+                (1+(XFVELOS/(2.*PLBDS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) )
       ZZW(:) =    ZZW(:)*(0.5+SIGN(0.5,ZZW(:))) - ABS(ZZW(:))*(0.5-SIGN(0.5,ZZW(:)))
    END WHERE
 !
diff --git a/src/MNH/lima_tendencies.f90 b/src/MNH/lima_tendencies.f90
index f9bd9f5856905cc640d6498a6a2e38fbe4054de5..96fda1b108190fb59c1c368a355a8be98336a387 100644
--- a/src/MNH/lima_tendencies.f90
+++ b/src/MNH/lima_tendencies.f90
@@ -233,6 +233,7 @@ SUBROUTINE LIMA_TENDENCIES (PTSTEP, LDCOMPUTE,
 !!      Original             15/03/2018
 !!
 !       Delbeke/Vie     03/2022 : KHKO option
+!       J. Wurtz        03/2022 : new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -240,11 +241,12 @@ SUBROUTINE LIMA_TENDENCIES (PTSTEP, LDCOMPUTE,
 !
 USE MODD_CST,              ONLY : XP00, XRD, XRV, XMD, XMV, XCPD, XCPV, XCL, XCI, XLVTT, XLSTT, XTT, &
                                   XALPW, XBETAW, XGAMW, XALPI, XBETAI, XGAMI
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN,                                                    &
-                                  LCOLD, LNUCL, LSNOW, LHAIL, LWARM, LACTI, LRAIN, LKHKO
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN, XNUS,                                    &
+                                  LCOLD, LNUCL, LSNOW, LHAIL, LWARM, LACTI, LRAIN, LKHKO, LSNOW_T
 USE MODD_PARAM_LIMA_WARM,  ONLY : XLBC, XLBEXC, XLBR, XLBEXR
 USE MODD_PARAM_LIMA_MIXED, ONLY : XLBG, XLBEXG, XLBH, XLBEXH, XLBDAG_MAX
-USE MODD_PARAM_LIMA_COLD,  ONLY : XSCFAC, XLBI, XLBEXI, XLBS, XLBEXS, XLBDAS_MAX
+USE MODD_PARAM_LIMA_COLD,  ONLY : XSCFAC, XLBI, XLBEXI, XLBS, XLBEXS, XLBDAS_MAX, XTRANS_MP_GAMMAS,  &
+                                  XFVELOS, XLBDAS_MIN
 !
 USE MODI_LIMA_DROPLETS_HOM_FREEZING
 USE MODI_LIMA_DROPLETS_SELF_COLLECTION
@@ -505,9 +507,21 @@ WHERE (PRIT(:)>XRTMIN(4) .AND. PCIT(:)>XCTMIN(4) .AND. LDCOMPUTE(:))
    ZLBDI(:) = ( XLBI*PCIT(:) / PRIT(:) )**XLBEXI
 END WHERE
 ZLBDS(:)  = 1.E10
-WHERE (PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:) )
-   ZLBDS(:) = XLBS*( PRHODREF(:)*PRST(:) )**XLBEXS
-END WHERE
+IF (LSNOW_T) THEN
+   WHERE (PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:) )
+      WHERE(ZT(:)>263.15)
+         ZLBDS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(:))),XLBDAS_MIN)
+      END WHERE
+      WHERE(ZT(:)<=263.15)
+         ZLBDS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(:))),XLBDAS_MIN)
+      END WHERE
+   END WHERE
+   ZLBDS(:) =  ZLBDS(:) * XTRANS_MP_GAMMAS
+ELSE
+   WHERE (PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:) )
+      ZLBDS(:) = XLBS*( PRHODREF(:)*PRST(:) )**XLBEXS
+   END WHERE
+END IF
 ZLBDG(:)  = 1.E10
 WHERE (PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:) )
    ZLBDG(:) = XLBG*( PRHODREF(:)*PRGT(:) )**XLBEXG
@@ -644,7 +658,7 @@ END IF
 ! Lambda_s limited for collection processes to prevent too high concentrations
 ! must be changed or removed if C and x modified
 !
-ZLBDS(:) = MIN( XLBDAS_MAX, ZLBDS(:))
+!ZLBDS(:) = MIN( XLBDAS_MAX, ZLBDS(:))
 !
 !
 IF (LCOLD .AND. LSNOW) THEN
diff --git a/src/MNH/modd_param_lima.f90 b/src/MNH/modd_param_lima.f90
index e12d84506bdaf4fbd820e161580996c8efe30087..0bb37d505a5c2da679f0298f9844992a0c990ce6 100644
--- a/src/MNH/modd_param_lima.f90
+++ b/src/MNH/modd_param_lima.f90
@@ -80,6 +80,7 @@ REAL, DIMENSION(:),    SAVE, ALLOCATABLE :: XFRAC_REF       ! AP compostion in P
 !
 ! 1.3 Ice characteristics
 !
+LOGICAL, SAVE :: LSNOW_T                     ! TRUE to enable snow param. after Wurtz 2021
 CHARACTER(LEN=4), SAVE :: CPRISTINE_ICE_LIMA ! Pristine type PLAT, COLU or BURO
 CHARACTER(LEN=4), SAVE :: CHEVRIMED_ICE_LIMA ! Heavily rimed type GRAU or HAIL
 REAL,SAVE              :: XALPHAI,XNUI,    & ! Pristine ice   distribution parameters
diff --git a/src/MNH/modd_param_lima_cold.f90 b/src/MNH/modd_param_lima_cold.f90
index 64494219e13b43678a4600d0f3de3c5a3291241c..65cb7f9a69a6b3e9d2a88e3acfb26ad319f95faa 100644
--- a/src/MNH/modd_param_lima_cold.f90
+++ b/src/MNH/modd_param_lima_cold.f90
@@ -19,6 +19,7 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             ??/??/13 
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 USE MODD_PARAMETERS, ONLY: JPSVNAMELGTMAX
@@ -51,8 +52,9 @@ REAL,SAVE :: XAI,XBI,XC_I,XDI         ,XF0I,XF2I,XC1I ! Cloud ice      charact.
 REAL,SAVE ::                           XF0IS,XF1IS    ! (large Di vent. coef.)
 REAL,SAVE :: XAS,XBS,XCS,XDS,XCCS,XCXS,XF0S,XF1S,XC1S ! Snow/agg.      charact.
 !
-REAL,SAVE :: XLBDAS_MAX               ! Max values allowed for the shape
-                                      ! parameter of snow
+REAL,SAVE :: XLBDAS_MIN, XLBDAS_MAX   ! Max values allowed for the shape parameter of snow
+REAL,SAVE :: XFVELOS                  ! Wurtz - snow fall speed parameterizaed after Thompson 2008
+REAL,SAVE :: XTRANS_MP_GAMMAS         ! Wurtz - change between lambda value for MP and gen. gamma
 !
 CHARACTER(LEN=JPSVNAMELGTMAX),DIMENSION(5),PARAMETER &
                               :: CLIMA_COLD_NAMES=(/'CICE    ','CIFNFREE','CIFNNUCL', &
diff --git a/src/MNH/modd_param_lima_mixed.f90 b/src/MNH/modd_param_lima_mixed.f90
index f13accfc669e88fca83566c2eb72f5c2cc6f4945..57ea4d1a559fcab6f202c4bafb3b6015119732f6 100644
--- a/src/MNH/modd_param_lima_mixed.f90
+++ b/src/MNH/modd_param_lima_mixed.f90
@@ -14,6 +14,7 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             ??/??/13 
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 !
@@ -67,6 +68,7 @@ REAL,SAVE :: XDCSLIM,XCOLCS,                   & ! Constants for the riming of
     	     XEXCRIMSS,XCRIMSS,                & ! the aggregates : RIM
     	     XEXCRIMSG,XCRIMSG,                & !
     	     XEXSRIMCG,XSRIMCG,                & !
+             XSRIMCG2, XSRIMCG3, XEXSRIMCG2,   & ! Murakami 1990
     	     XGAMINC_BOUND_MIN,                & ! Min val. of Lbda_s for RIM
     	     XGAMINC_BOUND_MAX,                & ! Max val. of Lbda_s for RIM
     	     XRIMINTP1,XRIMINTP2                 ! Csts for lin. interpol. of 
diff --git a/src/MNH/modn_param_lima.f90 b/src/MNH/modn_param_lima.f90
index b319bd9f4ed3e4701aeda3e0e9612b210b871b00..7678374f3b782b508b2604c85971ae7aa9f97f06 100644
--- a/src/MNH/modn_param_lima.f90
+++ b/src/MNH/modn_param_lima.f90
@@ -20,7 +20,7 @@ IMPLICIT NONE
 NAMELIST/NAM_PARAM_LIMA/LCOLD, LNUCL, LSEDI, LSNOW, LHAIL, LHHONI, LMEYERS,&
                         NMOD_IFN, XIFN_CONC, LIFN_HOM,                     &
                         CIFN_SPECIES, CINT_MIXING, NMOD_IMM, NIND_SPECIE,  &
-                        CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,            &
+                        LSNOW_T, CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,   &
                         XALPHAI, XNUI, XALPHAS, XNUS, XALPHAG, XNUG,       &
                         XFACTNUC_DEP, XFACTNUC_CON, NPHILLIPS,             &
 !
diff --git a/src/MNH/profilern.f90 b/src/MNH/profilern.f90
index 433a5c6ef02ae8a1630718f497b1f602b536f7a5..e2025c5784c58355b198a1fe13c185e5c1513a1f 100644
--- a/src/MNH/profilern.f90
+++ b/src/MNH/profilern.f90
@@ -112,10 +112,11 @@ USE MODE_FGAU,             ONLY : GAULAG
 USE MODE_FSCATTER,         ONLY: QEPSW,QEPSI,BHMIE,MOMG,MG
 USE MODD_PARAM_LIMA,       ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                                  XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC
+                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC, LSNOW_T_L=>LSNOW_T
 USE MODD_PARAM_LIMA_COLD,  ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
                                  XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,&
-                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS
+                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,        &
+                                 XLBDAS_MIN,XLBDAS_MAX
 USE MODD_PARAM_LIMA_MIXED, ONLY: XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,&
                                  XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
 USE MODD_PARAM_LIMA_WARM,  ONLY: XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XAR_L=>XAR,&
@@ -242,6 +243,7 @@ REAL                           :: ZAETOT,ZAETMP,ZREFLOC,ZQSCA,ZQBACK,ZQEXT ! tem
 REAL,DIMENSION(:),ALLOCATABLE  :: ZAELOC,ZZMZ ! temporary arrays
 INTEGER                        :: JPTS_GAULAG=9 ! number of points for Gauss-Laguerre quadrature
 REAL                           :: ZLBDA   ! slope distribution parameter
+REAL                           :: ZN   ! number cocentration
 REAL                           :: ZFRAC_ICE  ! ice water fraction
 REAL                           :: ZDELTA_EQUIV ! mass-equivalent Gauss-Laguerre point
 REAL                           :: ZFW ! liquid fraction
@@ -724,7 +726,17 @@ IF (GSTORE) THEN
                     ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
                   ENDIF
             END SELECT
-            ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+            IF (JLOOP.EQ.5 .AND. CCLOUD=='LIMA' .AND. LSNOW_T_L) THEN
+               IF (ZTEMPZ(JK)>-10.) THEN
+                  ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+               ELSE
+                  ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+               END IF
+               ZN=XLBS_L*ZRHODREFZ(JK)*ZRZ(JK,JLOOP)*ZLBDA**ZB
+            ELSE
+               ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+               ZN=ZCC*ZLBDA**ZCX
+            END IF
             ZREFLOC=0.
             ZAETMP=0.
             DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
@@ -754,8 +766,8 @@ IF (GSTORE) THEN
                ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
                ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
              END DO
-             ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
-             ZAETMP=ZAETMP  *           XPI    *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
+             ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZN/(4.*GAMMA(ZNU)*.93)
+             ZAETMP=ZAETMP  *           XPI    *ZN/(4.*GAMMA(ZNU))
              TPROFILER%CRARE(IN,JK,I)=TPROFILER%CRARE(IN,JK,I)+ZREFLOC
              ZAELOC(JK)=ZAELOC(JK)+ZAETMP
            END IF
diff --git a/src/MNH/radar_rain_ice.f90 b/src/MNH/radar_rain_ice.f90
index eddac2294d8c62054116b7be77081b610e72aefb..b42848f56435fe5c695e7e015d2498379abe4ac1 100644
--- a/src/MNH/radar_rain_ice.f90
+++ b/src/MNH/radar_rain_ice.f90
@@ -114,7 +114,8 @@ USE MODD_PARAM_LIMA_WARM, ONLY: XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XAR_L=>
                                 XBC_L=>XBC,XAC_L=>XAC,XCR_L=>XCR,XDR_L=>XDR
 USE MODD_PARAM_LIMA_COLD, ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
                                 XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,&
-                                XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,XCS_L=>XCS,XDS_L=>XDS
+                                XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,XCS_L=>XCS,XDS_L=>XDS,&
+                                XLBDAS_MIN,XLBDAS_MAX
 
 USE MODD_PARAM_LIMA_MIXED, ONLY:XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,&
                                 XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG,&
@@ -123,7 +124,7 @@ USE MODD_PARAM_LIMA_MIXED, ONLY:XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=
 
 USE MODD_PARAM_LIMA, ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                            XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                           XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC                      
+                           XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC, LSNOW_T_L=>LSNOW_T                      
 USE MODD_PARAMETERS
 USE MODD_PARAM_n, ONLY : CCLOUD
 USE MODD_LUNIT
@@ -169,6 +170,8 @@ REAL :: ZRHO00                        ! Surface reference air density
 LOGICAL, DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: GRAIN
 REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZLBDA 
                                       ! slope distribution parameter
+REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZN 
+                                      ! number concentration
 REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZW
 REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZREFL_MELT_CONV
 INTEGER                                                       :: JLBDA
@@ -325,7 +328,24 @@ END IF
 !               ---------------
 !
 IF (SIZE(PRT,4) >= 5) THEN
-  IF (CCLOUD=='LIMA') THEN
+  IF (CCLOUD=='LIMA' .AND. LSNOW_T_L) THEN
+    ZDMELT_FACT = ( (6.0*XAS_L)/(XPI*XRHOLW) )**(2.0)
+    ZEXP = 2.0*XBS_L
+    WHERE(PTEMP(:,:,:)>-10. .AND. PRT(:,:,:,5).GT.XRTMIN_L(5))
+       ZLBDA(:,:,:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(PTEMP(:,:,:)+273.15))),XLBDAS_MIN)
+    END WHERE
+    WHERE(PTEMP(:,:,:)<=-10 .AND. PRT(:,:,:,5).GT.XRTMIN_L(5))
+       ZLBDA(:,:,:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(PTEMP(:,:,:)+273.15))),XLBDAS_MIN)
+    END WHERE
+    ZN(:,:,:)=XLBS_L*PRHODREF(:,:,:)*PRT(:,:,:,5)*ZLBDA(:,:,:)**XBS_L
+    WHERE( PRT(:,:,:,5).GT.XRTMIN_L(5) )
+      ZW(:,:,:) = ZEQICE*ZDMELT_FACT                                             &
+                  *1.E18*ZN(:,:,:)*(ZLBDA(:,:,:)**(-ZEXP))*MOMG(XALPHAS_L,XNUS_L,ZEXP)
+      PVDOP(:,:,:) = PVDOP(:,:,:)+ZEQICE*ZDMELT_FACT*MOMG(XALPHAS_L,XNUS_L,ZEXP+XDS_L) &
+                     *1.E18*ZN(:,:,:)*XCS_L*(ZLBDA(:,:,:)**(-ZEXP-XDS_L))
+      PRARE(:,:,:) = PRARE(:,:,:) + ZW(:,:,:)
+    END WHERE
+  ELSEIF (CCLOUD=='LIMA') THEN
     ZDMELT_FACT = ( (6.0*XAS_L)/(XPI*XRHOLW) )**(2.0)
     ZEXP = 2.0*XBS_L
     WHERE( PRT(:,:,:,5).GT.XRTMIN_L(5) )
diff --git a/src/MNH/radar_scattering.f90 b/src/MNH/radar_scattering.f90
index dc3ddc1e5027eed57ce8610829eef632598ac9ac..413bbbb92e5c520ecce803e2b9dfd9e45a4ded9e 100644
--- a/src/MNH/radar_scattering.f90
+++ b/src/MNH/radar_scattering.f90
@@ -119,12 +119,13 @@ USE MODD_RAIN_ICE_DESCR, ONLY: XALPHAR_I=>XALPHAR,XNUR_I=>XNUR,XDR_I=>XDR,XLBEXR
 !!LIMA         
 USE MODD_PARAM_LIMA_WARM, ONLY: XDR_L=>XDR,XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XCR_L=>XCR
 USE MODD_PARAM_LIMA_COLD, ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
-                                XDS_L=>XDS,XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,XCS_L=>XCS
+                                XDS_L=>XDS,XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,XAS_L=>XAS,XBS_L=>XBS,&
+                                XCXS_L=>XCXS,XCS_L=>XCS,XLBDAS_MIN,XLBDAS_MAX
 
 USE MODD_PARAM_LIMA_MIXED, ONLY:XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
 USE MODD_PARAM_LIMA, ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                            XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                           XRTMIN_L=>XRTMIN
+                           XRTMIN_L=>XRTMIN, LSNOW_T_L=>LSNOW_T
 !!LIMA
 USE MODD_RADAR, ONLY:XLAM_RAD,XSTEP_RAD,NBELEV,NDIFF,LATT,NPTS_GAULAG,LQUAD,XVALGROUND,NDGS, &
      LFALL,LWBSCS,LWREFL,XREFLVDOPMIN,XREFLMIN,LSNRT,XSNRMIN
@@ -192,6 +193,7 @@ REAL :: ZDMELT_FACT ! factor used to compute the equivalent melted diameter
 REAL :: ZEQICE=0.224! factor used to convert the ice crystals reflectivity into an equivalent  liquid water reflectivity (from Smith, JCAM 84)
 REAL :: ZEXP        ! anciliary parameter
 REAL :: ZLBDA   ! slope distribution parameter
+REAL :: ZN      ! Number concentration
 REAL :: ZFRAC_ICE,ZD,ZDE ! auxiliary variables
 REAL :: ZQSCA
 REAL,DIMENSION(2) :: ZQEXT
@@ -1282,28 +1284,37 @@ DO JI=1,INBRAD
                     ZDMELT_FACT=6.*ZAS/(XPI*.92*XRHOLW)
                     ZEXP=2.*ZBS !XBS = 1.9 in ini_radar.f90 (bj tab 2.1 p24)
                     !dans ini_rain_ice.f90 :
-                    ZLBDA= ZLBS*(ZM)**ZLBEXS
-
+                    IF (GLIMA .AND. LSNOW_T_L) THEN
+                       IF (PT_RAY(JI,JEL,JAZ,JL,JH,JV)>-10.) THEN
+                          ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(PT_RAY(JI,JEL,JAZ,JL,JH,JV)+273.15))),XLBDAS_MIN)
+                       ELSE
+                          ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(PT_RAY(JI,JEL,JAZ,JL,JH,JV)+273.15))),XLBDAS_MIN)
+                       END IF
+                       ZN=ZLBS*ZM*ZLBDA**ZBS
+                    ELSE
+                       ZLBDA= ZLBS*(ZM)**ZLBEXS
+                       ZN=ZCCS*ZLBDA**ZCXS
+                    END IF
                     ! Rayleigh or Rayleigh-Gans or Rayleigh with 6th order for attenuation
                     IF(NDIFF==0.OR.NDIFF==3.OR.NDIFF==4) THEN
-                      ZREFLOC(1:2)=ZEQICE*.92**2*ZDMELT_FACT**2*1.E18*ZCCS*ZLBDA**(ZCXS-ZEXP)*MOMG(ZALPHAS,ZNUS,ZEXP)
+                      ZREFLOC(1:2)=ZEQICE*.92**2*ZDMELT_FACT**2*1.E18*ZN*ZLBDA**(ZEXP)*MOMG(ZALPHAS,ZNUS,ZEXP)
                       ZREFLOC(3)=0.
                       IF(LWREFL) THEN ! weighting by reflectivities
                         ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
                                -ZCS*SIN(PELEV(JI,JEL,JL,JV))*ZEQICE*.92**2*ZDMELT_FACT**2&
-                               *1.E18*ZCCS*ZLBDA**(ZCXS-ZEXP-ZDS)*MOMG(ZALPHAS,ZNUS,ZEXP+ZDS)
+                               *1.E18*ZN*ZLBDA**(ZEXP-ZDS)*MOMG(ZALPHAS,ZNUS,ZEXP+ZDS)
                       ELSE
-                        ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)+ZCCS*ZLBDA**ZCXS
+                        ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)+ZN
                         ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
                                                             -ZCS*SIN(PELEV(JI,JEL,JL,JV))&
-                                           *ZCCS*ZLBDA**(ZCXS-ZDS)*MOMG(ZALPHAS,ZNUS,ZDS)
+                                           *ZN*ZLBDA**(ZDS)*MOMG(ZALPHAS,ZNUS,ZDS)
                       END IF
                       IF(LATT) THEN
                         IF(NDIFF==0.OR.NDIFF==3) THEN
-                          ZAETMP(:)=ZCCS*ZLBDA**ZCXS*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)&
+                          ZAETMP(:)=ZN*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)&
                                     *MOMG(ZALPHAS,ZNUS,ZBS)/ZLBDA**ZBS)
                         ELSE
-                          ZAETMP(:)=ZCCS*ZLBDA**ZCXS*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)      &
+                          ZAETMP(:)=ZN*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)      &
                                     *MOMG(ZALPHAS,ZNUS,ZBS)/ZLBDA**ZBS                               &
                                     +ZDMELT_FACT**(5./3.)*XPI**4/15./XLAM_RAD(JI)**3                 &
                                     *AIMAG(ZQK**2*(ZQMI**4+27.*ZQMI**2+38.)                             &
@@ -1355,7 +1366,7 @@ DO JI=1,INBRAD
                       ZREFLOC(3)=180.E3/XPI*XLAM_RAD(JI)*ZRE_S22FMS11F
                       ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
                                                         -ZCS*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(1) &
-                                         *1.E18*(XLAM_RAD(JI)/XPI)**4/.93*ZCCS/4./ZLBDA**(3+ZDS)
+                                         *1.E18*(XLAM_RAD(JI)/XPI)**4/.93*(ZN*ZLBDA**(-ZCXS))/4./ZLBDA**(3+ZDS)
                       IF(LATT) THEN
                         ZAETMP(1)=ZIM_S22FT*XLAM_RAD(JI)*2
                         ZAETMP(2)=ZIM_S11FT*XLAM_RAD(JI)*2
@@ -1374,15 +1385,15 @@ DO JI=1,INBRAD
                         ZREFLOC(4)=ZREFLOC(4)+ZQBACK(1)*ZX(JJ)**(ZNUS-1.+2.*ZBS/3./ZALPHAS+ZDS/ZALPHAS)*ZW(JJ)
                         IF(LATT) ZAETMP(:)=ZAETMP(:)+ZQEXT(:)*ZX(JJ)**(ZNUS-1.+2.*ZBS/3./ZALPHAS)*ZW(JJ)
                       END DO ! ****** end loop Gauss-Laguerre quadrature
-                      ZREFLOC(1:2)=1.E18*(XLAM_RAD(JI)/XPI)**4*ZCCS*ZLBDA**(ZCXS-2.*ZBS/3.)/&
+                      ZREFLOC(1:2)=1.E18*(XLAM_RAD(JI)/XPI)**4*ZN*ZLBDA**(-2.*ZBS/3.)/&
                                   (4.*GAMMA(ZNUS)*.93)*ZDMELT_FACT**(2./3.)*ZREFLOC(1:2)               
                       ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
                                             +PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
                                             -ZCS*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(4)   &
-                                            *1.E18*(XLAM_RAD(JI)/XPI)**4*ZCCS          &
-                                            *ZLBDA**(ZCXS-2.*ZBS/3.-ZDS)/              &
+                                            *1.E18*(XLAM_RAD(JI)/XPI)**4*ZN          &
+                                            *ZLBDA**(2.*ZBS/3.-ZDS)/              &
                                             (4.*GAMMA(ZNUS)*.93)*ZDMELT_FACT**(2./3.)                              
-                      IF(LATT) ZAETMP(:)=ZAETMP(:)*XPI*ZCCS*ZLBDA**(ZCXS-2.*ZBS/3.)/(4.*GAMMA(ZNUS))&
+                      IF(LATT) ZAETMP(:)=ZAETMP(:)*XPI*ZN*ZLBDA**(-2.*ZBS/3.)/(4.*GAMMA(ZNUS))&
                                          *ZDMELT_FACT**(2./3.)            
                       ZRE_S22S11_S=0
                       ZIM_S22S11_S=0
diff --git a/src/MNH/rrcolss.f90 b/src/MNH/rrcolss.f90
index 527165111ecf4d225ce5ec0117c09846d2116b9e..00e7b6cc9cb12f8591f36879a97677d8a55fd5f9 100644
--- a/src/MNH/rrcolss.f90
+++ b/src/MNH/rrcolss.f90
@@ -10,7 +10,7 @@
 INTERFACE
 !
       SUBROUTINE RRCOLSS( KND, PALPHAS, PNUS, PALPHAR, PNUR,                 &
-                         PESR, PEXMASSR, PFALLS, PEXFALLS, PFALLR, PEXFALLR, &
+                         PESR, PEXMASSR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, & ! Wurtz
                          PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
                          PDINFTY, PRRCOLSS, PAG, PBS, PAS                    )
 !
@@ -28,6 +28,7 @@ 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
@@ -49,7 +50,7 @@ END INTERFACE
       END MODULE MODI_RRCOLSS
 !     ########################################################################
       SUBROUTINE RRCOLSS( KND, PALPHAS, PNUS, PALPHAR, PNUR,                 &
-                         PESR, PEXMASSR, PFALLS, PEXFALLS, PFALLR, PEXFALLR, &
+                         PESR, PEXMASSR, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
                          PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
                          PDINFTY, PRRCOLSS, PAG, PBS, PAS                    )
 !     ########################################################################
@@ -117,6 +118,7 @@ END INTERFACE
 !!      Original    8/11/95
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -151,6 +153,7 @@ 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
@@ -277,11 +280,11 @@ DO JLBDAS = 1,SIZE(PRRCOLSS(:,:),1)
           DO JDR = 1,INR-1
             ZDR = ZDDCOLLR * REAL(JDR)
             ZCOLLR = ZCOLLR + (ZDS+ZDR)**2 * ZDR**PEXMASSR                     &
-                       * PESR * ABS(PFALLS*ZDS**PEXFALLS-PFALLR*ZDR**PEXFALLR) &
+                       * PESR * ABS(PFALLS*ZDS**PEXFALLS * EXP(-(PFALLEXPS*XDS)**PALPHAS)-PFALLR*ZDR**PEXFALLR) &
                                       * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
           END DO
           ZCOLLDRMAX = (ZDS+ZDRMAX)**2 * ZDRMAX**PEXMASSR                      &
-                    * PESR * ABS(PFALLS*ZDS**PEXFALLS-PFALLR*ZDRMAX**PEXFALLR) &
+                    * PESR * ABS(PFALLS*ZDS**PEXFALLS* EXP(-(PFALLEXPS*XDS)**PALPHAS)-PFALLR*ZDRMAX**PEXFALLR) &
                                    * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDRMAX)
           ZCOLLR = (ZCOLLR + 0.5*ZCOLLDRMAX)*(ZDDCOLLR/ZDDSCALR)
 !
diff --git a/src/MNH/rscolrg.f90 b/src/MNH/rscolrg.f90
index caa868e91d39cbe12010bfa2c265ffe35304dba4..210df03805b8189e969d1712d411b6f7540e07d8 100644
--- a/src/MNH/rscolrg.f90
+++ b/src/MNH/rscolrg.f90
@@ -10,7 +10,7 @@
 INTERFACE
 !
       SUBROUTINE RSCOLRG( KND, PALPHAS, PZNUS, PALPHAR, PNUR,                &
-                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLR, PEXFALLR, &
+                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
                          PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
                          PDINFTY, PRSCOLRG,PAG, PBS, PAS                     )
 !
@@ -28,6 +28,7 @@ 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
@@ -49,7 +50,7 @@ END INTERFACE
       END MODULE MODI_RSCOLRG
 !     ########################################################################
       SUBROUTINE RSCOLRG( KND, PALPHAS, PZNUS, PALPHAR, PNUR,                &
-                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLR, PEXFALLR, &
+                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
                          PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
                          PDINFTY, PRSCOLRG,PAG, PBS, PAS                     )
 !     ########################################################################
@@ -117,6 +118,7 @@ END INTERFACE
 !!      Original    8/11/95
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -149,6 +151,7 @@ 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
@@ -271,12 +274,12 @@ DO JLBDAR = 1,SIZE(PRSCOLRG(:,:),1)
             ZDR = ZDDCOLLR * REAL(JDR) + ZDRMIN
             ZCOLLR = ZCOLLR + (ZDS+ZDR)**2                                     &
                        * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)                &
-                         * PESR * ABS(PFALLS*ZDS**PEXFALLS-PFALLR*ZDR**PEXFALLR)
+                         * 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-PFALLR*ZDRMIN**PEXFALLR)
+                      * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDRMIN**PEXFALLR)
             ELSE
             ZCOLLDRMIN = 0.0
           END IF 
diff --git a/src/MNH/rzcolx.f90 b/src/MNH/rzcolx.f90
index 28658241cf1021a29de694cd5a99b85e9c3340d9..d3977af63a8cccde05450cd09931f1163e82a61e 100644
--- a/src/MNH/rzcolx.f90
+++ b/src/MNH/rzcolx.f90
@@ -10,7 +10,8 @@
 INTERFACE
 !
       SUBROUTINE RZCOLX( KND, PALPHAX, PNUX, PALPHAZ, PNUZ,                  &
-			 PEXZ, PEXMASSZ, PFALLX, PEXFALLX, PFALLZ, PEXFALLZ, &
+                         PEXZ, PEXMASSZ, PFALLX, PEXFALLX, PFALLEXPX,        &
+                         PFALLZ, PEXFALLZ, PFALLEXPZ,                        &
 		         PLBDAXMAX, PLBDAZMAX, PLBDAXMIN, PLBDAZMIN,         &
 		         PDINFTY, PRZCOLX                                    )
 !
@@ -29,8 +30,10 @@ REAL, INTENT(IN) :: PEXZ      ! Efficiency of specy X collecting specy Z
 REAL, INTENT(IN) :: PEXMASSZ  ! Mass exponent of 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
@@ -49,7 +52,8 @@ END INTERFACE
       END MODULE MODI_RZCOLX
 !     ########################################################################
       SUBROUTINE RZCOLX( KND, PALPHAX, PNUX, PALPHAZ, PNUZ,                  &
-			 PEXZ, PEXMASSZ, PFALLX, PEXFALLX, PFALLZ, PEXFALLZ, &
+                         PEXZ, PEXMASSZ, PFALLX, PEXFALLX, PFALLEXPX,        &
+                         PFALLZ, PEXFALLZ, PFALLEXPZ,                        &
 		         PLBDAXMAX, PLBDAZMAX, PLBDAXMIN, PLBDAZMIN,         &
 		         PDINFTY, PRZCOLX                                    )
 !     ########################################################################
@@ -121,6 +125,7 @@ END INTERFACE
 !!      Original    8/11/95
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -152,8 +157,10 @@ REAL, INTENT(IN) :: PEXZ      ! Efficiency of specy X collecting specy Z
 REAL, INTENT(IN) :: PEXMASSZ  ! Mass exponent of 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
@@ -234,20 +241,20 @@ DO JLBDAX = 1,SIZE(PRZCOLX(:,:),1)
       ZSCALZ = 0.0
       ZCOLLZ = 0.0
       DO JDZ = 1,KND-1
-        ZDZ = ZDDZ * REAL(JDZ)
+         ZDZ = ZDDZ * REAL(JDZ)
 !
 !*       1.6     Compute the normalization factor by integration over the
 !                dimensional spectrum of specy Z  
 !
-	ZFUNC  = (ZDX+ZDZ)**2 * ZDZ**PEXMASSZ                            &
-			      * GENERAL_GAMMA(PALPHAZ,PNUZ,ZLBDAZ,ZDZ)
-	ZSCALZ = ZSCALZ + ZFUNC
+         ZFUNC  = (ZDX+ZDZ)**2 * ZDZ**PEXMASSZ                            &
+		               * 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-PFALLZ*ZDZ**PEXFALLZ)
+         ZCOLLZ = ZCOLLZ + ZFUNC * PEXZ * ABS( PFALLX*ZDX**PEXFALLX * EXP(-(ZDX*PFALLEXPX)**PALPHAX) &
+                                             - PFALLZ*ZDZ**PEXFALLZ * EXP(-(ZDZ*PFALLEXPZ)**PALPHAZ)) ! Wurtz
       END DO
 !
 !*       1.8     Compute the normalization factor by integration over the
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index 967b74d5da18fc6a2131f3406a9425d62043a326..7ca40d7a41b47f6939a81a3d9e56f935f12a9457 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -3926,20 +3926,20 @@ IF (LLIDAR) THEN
     ZTMP3(:,:,:,1)=ZSIG_DST(:,:,:,IACCMODE)
     SELECT CASE ( CCLOUD )
     CASE('KESS''ICE3','ICE4')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
                  XRT, ZWORK31, ZWORK32,                                        &
                  PDSTC=ZTMP1,                                                  &
                  PDSTD=ZTMP2,                                                  &
                  PDSTS=ZTMP3)
     CASE('C2R2')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
                  XRT, ZWORK31, ZWORK32,                                        &
                  PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C2R2END),                    &
                  PDSTC=ZTMP1,                                                  &
                  PDSTD=ZTMP2,                                                  &
                  PDSTS=ZTMP3)
     CASE('C3R5')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
                  XRT, ZWORK31, ZWORK32,                                        &
                  PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C1R3END-1),                  &
                  PDSTC=ZTMP1,                                                  &
@@ -3953,7 +3953,7 @@ IF (LLIDAR) THEN
        ZTMP4(:,:,:,3)=XSVT(:,:,:,NSV_LIMA_NR)
        ZTMP4(:,:,:,4)=XSVT(:,:,:,NSV_LIMA_NI)
 !
-       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, MAX(XCLDFR,XICEFR),&
+       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, MAX(XCLDFR,XICEFR),&
             XRT, ZWORK31, ZWORK32,                                  &
             PCT=ZTMP4,                            &
             PDSTC=ZTMP1,                          &
@@ -3964,14 +3964,14 @@ IF (LLIDAR) THEN
   ELSE
     SELECT CASE ( CCLOUD )
     CASE('KESS','ICE3','ICE4')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
            XRT, ZWORK31, ZWORK32)
     CASE('C2R2')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
            XRT, ZWORK31, ZWORK32,                                  &
            PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C2R2END))
     CASE('C3R5')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
            XRT, ZWORK31, ZWORK32,                                  &
            PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C1R3END-1))
     CASE('LIMA')
@@ -3982,7 +3982,7 @@ IF (LLIDAR) THEN
        ZTMP4(:,:,:,3)=XSVT(:,:,:,NSV_LIMA_NR)
        ZTMP4(:,:,:,4)=XSVT(:,:,:,NSV_LIMA_NI)
 !
-       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, MAX(XCLDFR,XICEFR),&
+       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, MAX(XCLDFR,XICEFR),&
             XRT, ZWORK31, ZWORK32,                                  &
             PCT=ZTMP4)
     END SELECT