diff --git a/src/MNH/ice4_fast_rg.f90 b/src/MNH/ice4_fast_rg.f90
index c78a154bb43fe2494c7a86baf99a7a0df7676e67..0f7dea2179ad05277f16065f54b35c893dbb4deb 100644
--- a/src/MNH/ice4_fast_rg.f90
+++ b/src/MNH/ice4_fast_rg.f90
@@ -94,6 +94,7 @@ SUBROUTINE ICE4_FAST_RG(KSIZE, LDSOFT, PCOMPUTE, KRR, &
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -335,14 +336,7 @@ ELSE
     END DO
     !
     WHERE(GDRY(:))
-!      PRG_TEND(:, IRSWETG)=XFSDRYG*ZZW(:)                         & ! RSDRYG
-!                                    / XCOLSG &
-!                  *(PLBDAS(:)**(XCXS-XBS))*( PLBDAG(:)**XCXG )    &
-!                  *(PRHODREF(:)**(-XCEXVT-1.))                    &
-!                       *( XLBSDRYG1/( PLBDAG(:)**2              ) + &
-!                          XLBSDRYG2/( PLBDAG(:)   * PLBDAS(:)   ) + &
-!                          XLBSDRYG3/(               PLBDAS(:)**2))
-      PRG_TEND(:, IRSWETG)=XFSDRYG*ZZW(:)                         & ! RSDRYG  ! Modif Wurtz snow diag
+      PRG_TEND(:, IRSWETG)=XFSDRYG*ZZW(:)                         & ! RSDRYG
                                     / XCOLSG &
                   *(PRST(:))*( PLBDAG(:)**XCXG )    &
                   *(PRHODREF(:)**(-XCEXVT))                    &
diff --git a/src/MNH/ice4_fast_rh.f90 b/src/MNH/ice4_fast_rh.f90
index 98106308dc0fd54a23ee2ce45f3e174e36f66f78..a21c85f19308cffede426f7b2510f52ef5cd447a 100644
--- a/src/MNH/ice4_fast_rh.f90
+++ b/src/MNH/ice4_fast_rh.f90
@@ -84,6 +84,7 @@ SUBROUTINE ICE4_FAST_RH(KSIZE, LDSOFT, PCOMPUTE, PWETG, &
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -269,16 +270,10 @@ ELSE
     END DO
     !
     WHERE(GWET(:))
-!      PRH_TEND(:, IRSWETH)=XFSWETH*ZZW(:)                       & ! RSWETH
-!                    *( PLBDAS(:)**(XCXS-XBS) )*( PLBDAH(:)**XCXH )  &
-!                       *( PRHODREF(:)**(-XCEXVT-1.) )               &
-!                       *( XLBSWETH1/( PLBDAH(:)**2              ) + &
-!                          XLBSWETH2/( PLBDAH(:)   * PLBDAS(:)   ) + &
-!                          XLBSWETH3/(               PLBDAS(:)**2) )
-      PRH_TEND(:, IRSWETH)=XFSWETH*ZZW(:)                       & ! RSWETH Modif Wurtz conc snow
+      PRH_TEND(:, IRSWETH)=XFSWETH*ZZW(:)                       &
                     *( PRST(:))*( PLBDAH(:)**XCXH )  &
                        *( PRHODREF(:)**(-XCEXVT) )               &
-                       *( XLBSWETH1/( PLBDAH(:)**2              ) + &	! Il s'agit de moment (?)
+                       *( XLBSWETH1/( PLBDAH(:)**2              ) + &
                           XLBSWETH2/( PLBDAH(:)   * PLBDAS(:)   ) + &
                           XLBSWETH3/(               PLBDAS(:)**2) )
       PRH_TEND(:, IRSDRYH)=PRH_TEND(:, IRSWETH)*(XCOLSH*EXP(XCOLEXSH*(PT(:)-XTT)))
diff --git a/src/MNH/ice4_fast_rs.f90 b/src/MNH/ice4_fast_rs.f90
index b3df38c8c4bdf38c750c3d5a25deb2b83a73d266..50d0cae7d8fc8a15435f8a44f54d300a2ea567f4 100644
--- a/src/MNH/ice4_fast_rs.f90
+++ b/src/MNH/ice4_fast_rs.f90
@@ -77,6 +77,7 @@ SUBROUTINE ICE4_FAST_RS(KSIZE, LDSOFT, PCOMPUTE, &
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -244,13 +245,10 @@ ELSE
     !        5.1.4  riming of the small sized aggregates
     !
     WHERE (GRIM(:))
-!      PRS_TEND(:, IRCRIMSS) = XCRIMSS * ZZW(:) * PRCT(:)                & ! RCRIMSS
-!                                      *   PLBDAS(:)**XEXCRIMSS &
-!                                      * PRHODREF(:)**(-XCEXVT)
-      PRS_TEND(:, IRCRIMSS) = XCRIMSS * ZZW(:) * PRCT(:)                & ! RCRIMSS	 !Wurtz ! Thompson
+      PRS_TEND(:, IRCRIMSS) = XCRIMSS * ZZW(:) * PRCT(:)                & ! RCRIMSS
                                       * PRST(:)*(1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
                                       * PRHODREF(:)**(-XCEXVT+1.) &
-				      * (PLBDAS(:)) ** (XEXCRIMSS+XBS) ! Thompson
+				      * (PLBDAS(:)) ** (XEXCRIMSS+XBS)
     END WHERE
     !
     !        5.1.5  perform the linear interpolation of the normalized
@@ -275,21 +273,17 @@ ELSE
     !
     !
     WHERE(GRIM(:))
-!      PRS_TEND(:, IRCRIMS)=XCRIMSG * PRCT(:)               & ! RCRIMS
-!                                   * PLBDAS(:)**XEXCRIMSG  &
-!                                   * PRHODREF(:)**(-XCEXVT)
       PRS_TEND(:, IRCRIMS) = XCRIMSG * PRCT(:)               & ! RCRIMS
                                    * PRST(:)*(1+(XFVELOS/PLBDAS(:))**(XALPHAS))**(-XNUS+XEXCRIMSG/XALPHAS) &
                                    * PRHODREF(:)**(-XCEXVT+1.) &
-                                   * PLBDAS(:)**(XBS+XEXCRIMSG) ! GAMMAGEN LH_EXTENDED
+                                   * PLBDAS(:)**(XBS+XEXCRIMSG)
       ZZW6(:) = PRS_TEND(:, IRCRIMS) - PRS_TEND(:, IRCRIMSS) ! RCRIMSG
     END WHERE
 
     IF(CSNOWRIMING=='M90 ')THEN
       !Murakami 1990
       WHERE(GRIM(:))
-       ! PRS_TEND(:, IRSRIMCG)=XSRIMCG * PLBDAS(:)**XEXSRIMCG*(1.0-ZZW(:))
-        PRS_TEND(:, IRSRIMCG)=XSRIMCG * PRST(:)*PRHODREF(:)*PLBDAS(:)**(XEXSRIMCG+XBS)*(1.0-ZZW(:)) !Wurtz
+        PRS_TEND(:, IRSRIMCG)=XSRIMCG * PRST(:)*PRHODREF(:)*PLBDAS(:)**(XEXSRIMCG+XBS)*(1.0-ZZW(:))
 
         PRS_TEND(:, IRSRIMCG)=ZZW6(:)*PRS_TEND(:, IRSRIMCG)/ &
                        MAX(1.E-20, &
@@ -394,13 +388,8 @@ ELSE
     !        5.2.4  raindrop accretion on the small sized aggregates
     !
     WHERE(GACC(:))
-!      ZZW6(:) =                                                        & !! coef of RRACCS
-!            XFRACCSS*( PLBDAS(:)**XCXS )*( PRHODREF(:)**(-XCEXVT-1.) ) &
-!       *( XLBRACCS1/((PLBDAS(:)**2)               ) +                  &
-!          XLBRACCS2/( PLBDAS(:)    * PLBDAR(:)    ) +                  &
-!          XLBRACCS3/(               (PLBDAR(:)**2)) )/PLBDAR(:)**4
       ZZW6(:) =                                                        & !! coef of RRACCS
-            XFRACCSS*( PRST(:)*PLBDAS(:)**XBS )*( PRHODREF(:)**(-XCEXVT) ) & ! Wurtz
+            XFRACCSS*( PRST(:)*PLBDAS(:)**XBS )*( PRHODREF(:)**(-XCEXVT) ) &
        *( XLBRACCS1/((PLBDAS(:)**2)               ) +                  &
           XLBRACCS2/( PLBDAS(:)    * PLBDAR(:)    ) +                  &
           XLBRACCS3/(               (PLBDAR(:)**2)) )/PLBDAR(:)**4
@@ -445,12 +434,7 @@ ELSE
     !               into graupeln
     !
     WHERE(GACC(:))
-!      PRS_TEND(:, IRSACCRG) = XFSACCRG*ZZW(:)*                    & ! RSACCRG
-!          ( PLBDAS(:)**(XCXS-XBS) )*( PRHODREF(:)**(-XCEXVT-1.) ) &
-!         *( XLBSACCR1/((PLBDAR(:)**2)               ) +           &
-!            XLBSACCR2/( PLBDAR(:)    * PLBDAS(:)    ) +           &
-!            XLBSACCR3/(               (PLBDAS(:)**2)) )/PLBDAR(:)
-      PRS_TEND(:, IRSACCRG) = XFSACCRG*ZZW(:)*                    & ! RSACCRG	! Modif Wurtz
+      PRS_TEND(:, IRSACCRG) = XFSACCRG*ZZW(:)*                    & ! RSACCRG
           ( PRST(:))*( PRHODREF(:)**(-XCEXVT) ) &
          *( XLBSACCR1/((PLBDAR(:)**2)               ) +           &
             XLBSACCR2/( PLBDAR(:)    * PLBDAS(:)    ) +           &
@@ -516,13 +500,7 @@ ELSE
     !
     ! compute RSMLT
     !
-!    PRSMLTG(:)  = XFSCVMG*MAX( 0.0,( -PRSMLTG(:) *             &
-!                         ( X0DEPS*       PLBDAS(:)**XEX0DEPS +     &
-!                           X1DEPS*PCJ(:)*PLBDAS(:)**XEX1DEPS ) -   &
-!                                   ( PRS_TEND(:, IRCRIMS) + PRS_TEND(:, IRRACCS) ) *       &
-!                            ( PRHODREF(:)*XCL*(XTT-PT(:))) ) /    &
-!                                           ( PRHODREF(:)*XLMTT ) )
-    PRSMLTG(:)  = XFSCVMG*MAX( 0.0,( -PRSMLTG(:) *             & ! Modif GAMMAGEN LH_EXTENDED
+    PRSMLTG(:)  = XFSCVMG*MAX( 0.0,( -PRSMLTG(:) *             &
          PRST(:)*PRHODREF(:) *    &
          ( X0DEPS       *PLBDAS(:)**XEX0DEPS +     &
          X1DEPS*PCJ(:)*(1+(XFVELOS/(2*PLBDAS(:))**XALPHAS))**(XNUS+XEX1DEPS/XALPHAS)*((PLBDAS(:))**(XBS+XEX1DEPS))) -   &
diff --git a/src/MNH/ice4_rsrimcg_old.f90 b/src/MNH/ice4_rsrimcg_old.f90
index 5fc932536a176f4deaf51df66fc320b06f6afa18..d250ddc0e38953cbec90d1615655ba1eb291b313 100644
--- a/src/MNH/ice4_rsrimcg_old.f90
+++ b/src/MNH/ice4_rsrimcg_old.f90
@@ -44,6 +44,7 @@ SUBROUTINE ICE4_RSRIMCG_OLD(KSIZE, ODSOFT, ODCOMPUTE, &
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -131,9 +132,7 @@ IF(.NOT. ODSOFT) THEN
     !
     !
     WHERE(GRIM(:))
-!      PRSRIMCG_MR(:) = XSRIMCG * PLBDAS(:)**XEXSRIMCG   & ! RSRIMCG
-!                               * (1.0 - ZZW(:) )/PRHODREF(:)
-      PRSRIMCG_MR(:) = XSRIMCG * PLBDAS(:)**XEXSRIMCG   & ! RSRIMCG	!Modif Wurtz concentration snow
+      PRSRIMCG_MR(:) = XSRIMCG * PLBDAS(:)**XEXSRIMCG   & ! RSRIMCG
                                * (1.0 - ZZW(:) )*PRST(:)
       PRSRIMCG_MR(:)=MIN(PRST(:), PRSRIMCG_MR(:))
     END WHERE
diff --git a/src/MNH/ice4_sedimentation_split.f90 b/src/MNH/ice4_sedimentation_split.f90
index 9e93d59e577082e5ad1cb1785d7f379899cd5c2f..f9369a8783af23ad117f46ab51489a44afad0436 100644
--- a/src/MNH/ice4_sedimentation_split.f90
+++ b/src/MNH/ice4_sedimentation_split.f90
@@ -6,12 +6,12 @@
 MODULE MODI_ICE4_SEDIMENTATION_SPLIT
 INTERFACE
 SUBROUTINE ICE4_SEDIMENTATION_SPLIT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB, KKTE, KKT, KKL, &
-                                   &PTSTEP, KRR, OSEDIC, ODEPOSC, PVDEPOSC, PDZZ, &
-                                   &PRHODREF, PPABST, PTHT, PRHODJ, &
-				   & PLBDAS, &  ! Modif Wurtz                           
-                                   &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
-                                   &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG, &
-                                   &PSEA, PTOWN,  &
+                                   &PTSTEP, KRR, OSEDIC, ODEPOSC, PVDEPOSC, PDZZ,                 &
+                                   &PRHODREF, PPABST, PTHT, PRHODJ,                               &
+				   & PLBDAS,                                                      &
+                                   &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,   &
+                                   &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG,               &
+                                   &PSEA, PTOWN,                                                  &
                                    &PINPRH, PRHT, PRHS, PFPR)
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB, KKTE, KKT
@@ -26,7 +26,7 @@ REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODREF! Reference den
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PPABST  ! absolute pressure at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PTHT    ! Theta at time t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODJ  ! Dry density * Jacobian
-REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS ! lambda parameter for snow ! Modif Wurtz
+REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS  ! lambda parameter for snow
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRCS    ! Cloud water m.r. source
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRCT    ! Cloud water m.r. at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRRS    ! Rain water m.r. source
@@ -53,12 +53,12 @@ END SUBROUTINE ICE4_SEDIMENTATION_SPLIT
 END INTERFACE
 END MODULE MODI_ICE4_SEDIMENTATION_SPLIT
 SUBROUTINE ICE4_SEDIMENTATION_SPLIT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB, KKTE, KKT, KKL, &
-                                   &PTSTEP, KRR, OSEDIC, ODEPOSC, PVDEPOSC, PDZZ, &
-                                   &PRHODREF, PPABST, PTHT, PRHODJ, &
-				   & PLBDAS, &  ! Modif Wurtz                           
-                                   &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
-                                   &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG, &
-                                   &PSEA, PTOWN,  &
+                                   &PTSTEP, KRR, OSEDIC, ODEPOSC, PVDEPOSC, PDZZ,                 &
+                                   &PRHODREF, PPABST, PTHT, PRHODJ,                               &
+				   & PLBDAS,                                                      &
+                                   &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,   &
+                                   &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG,               &
+                                   &PSEA, PTOWN,                                                  &
                                    &PINPRH, PRHT, PRHS, PFPR)
 !!
 !!**  PURPOSE
@@ -75,6 +75,7 @@ SUBROUTINE ICE4_SEDIMENTATION_SPLIT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB
 !!
 !  P. Wautelet 11/02/2019: dimensions of PINPRC and PINDEP not necessarily KIT,KJT
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -105,7 +106,7 @@ REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODREF! Reference den
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PPABST  ! absolute pressure at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PTHT    ! Theta at time t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODJ  ! Dry density * Jacobian
-REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS ! lambda parameter for snow ! Modif Wurtz
+REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS  ! lambda parameter for snow
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRCS    ! Cloud water m.r. source
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRCT    ! Cloud water m.r. at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRRS    ! Rain water m.r. source
@@ -243,7 +244,7 @@ IF (GSEDIC) THEN
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &2, &
-			  &PLBDAS, & ! Modif Wurtz
+			  &PLBDAS, &
                           &ZRCT, PRCS, PINPRC, ZPRCS, &
                           &ZRAY, ZLBC, ZFSEDC, ZCONC3D, PFPR=PFPR)
 ENDIF
@@ -270,7 +271,7 @@ END IF
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &3, &
-			  &PLBDAS, & ! Modif Wurtz
+			  &PLBDAS, &
                           &ZRRT, PRRS, PINPRR, ZPRRS, &
                           &PFPR=PFPR)
 !
@@ -280,7 +281,7 @@ END IF
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &4, &
-			  &PLBDAS, & ! Modif Wurtz
+			  &PLBDAS, &
                           &ZRIT, PRIS, PINPRI, ZPRIS, &
                           PFPR=PFPR)
 !
@@ -290,7 +291,7 @@ END IF
                         &XSPLIT_MAXCFL, &
                         &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                         &5, &
-                        &PLBDAS, & ! Modif Wurtz
+                        &PLBDAS, &
                         &ZRST, PRSS, PINPRS, ZPRSS, &
                         PFPR=PFPR)
 !
@@ -300,7 +301,7 @@ END IF
                         &XSPLIT_MAXCFL, &
                         &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                         &6, &
-                        &PLBDAS, & ! Modif Wurtz
+                        &PLBDAS, &
                         &ZRGT, PRGS, PINPRG, ZPRGS, &
                         PFPR=PFPR)
 !
@@ -311,7 +312,7 @@ IF (IRR==7) THEN
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &7, &
-			  &PLBDAS, & ! Modif Wurtz
+			  &PLBDAS, &
                           &ZRHT, PRHS, PINPRH, ZPRHS, &
                           PFPR=PFPR)
 ENDIF
@@ -332,9 +333,9 @@ SUBROUTINE INTERNAL_SEDIM_SPLI(KIB,KIE,KIT,KJB,KJE,KJT,KKB,KKTB,KKTE,KKT,KKL,KRR
 !          ------------
 !
 USE MODD_CST,            ONLY: XCPD,XP00,XRD
-USE MODD_RAIN_ICE_DESCR, ONLY: XCC,XCEXVT,XDC,XLBEXC,XRTMIN,XALPHAS,XNUS,XBS
-USE MODD_RAIN_ICE_PARAM, ONLY: XEXCSEDI,XEXSEDG,XEXSEDH,XEXSEDR,XEXSEDS,XFSEDG,XFSEDH,XFSEDI,XFSEDR,XFSEDS,&
-                               XFVELOS
+USE MODD_RAIN_ICE_DESCR, ONLY: XCC,XCEXVT,XDC,XLBEXC,XRTMIN,XALPHAS,XNUS,XBS,XFVELOS
+USE MODD_RAIN_ICE_PARAM, ONLY: XEXCSEDI,XEXSEDG,XEXSEDH,XEXSEDR,XEXSEDS,XFSEDG,XFSEDH,XFSEDI,XFSEDR,XFSEDS
+                               
 !
 IMPLICIT NONE
 !
@@ -456,8 +457,8 @@ DO WHILE (ANY(ZREMAINT>0.))
 
         ZWSED(JI, JJ, JK) = XFSEDS *  &
                               & PRXT(JI,JJ,JK)* &
-                              & PRHODREF(JI,JJ,JK)**(1-XCEXVT) * & !    Modif Wurtz snow
-                              & (1 + (XFVELOS/PLBDAS(JI, JJ, JK))**XALPHAS)** (-XNUS+XEXSEDS/XALPHAS) * & ! GAMMAGEN LH_EXTENDED
+                              & PRHODREF(JI,JJ,JK)**(1-XCEXVT) * &
+                              & (1 + (XFVELOS/PLBDAS(JI, JJ, JK))**XALPHAS)** (-XNUS+XEXSEDS/XALPHAS) * &
 			      & PLBDAS(JI, JJ, JK) ** (XBS+XEXSEDS)
  ! GAMMAGEN_LH_EXTENDED
 
@@ -469,9 +470,6 @@ DO WHILE (ANY(ZREMAINT>0.))
       CASE(3)
         ZFSED=XFSEDR
         ZEXSED=XEXSEDR
-!      CASE(5)
-!        ZFSED=XFSEDS
-!        ZEXSED=XEXSEDS
       CASE(6)
         ZFSED=XFSEDG
         ZEXSED=XEXSEDG
diff --git a/src/MNH/ice4_sedimentation_stat.f90 b/src/MNH/ice4_sedimentation_stat.f90
index 4daa043dab99d656ba93cd941120e3e35c81619d..891223920fd73ea2bbbeed059ce1669aeebd2d51 100644
--- a/src/MNH/ice4_sedimentation_stat.f90
+++ b/src/MNH/ice4_sedimentation_stat.f90
@@ -8,7 +8,7 @@ INTERFACE
 SUBROUTINE ICE4_SEDIMENTATION_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB, KKTE, KKT, KKL, &
                                    &PTSTEP, KRR, OSEDIC, ODEPOSC, PVDEPOSC, PDZZ, &
                                    &PRHODREF, PPABST, PTHT, PRHODJ, &
-				  & PLBDAS, &  ! Modif Wurtz                    
+				  & PLBDAS, &     
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT,&
                                    &PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG, &
@@ -27,7 +27,7 @@ REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODREF! Reference den
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PPABST  ! absolute pressure at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PTHT    ! Theta at time t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODJ  ! Dry density * Jacobian
-REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS ! lambda parameter for snow ! Modif Wurtz
+REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS  ! lambda parameter for snow
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRCS    ! Cloud water m.r. source
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRCT    ! Cloud water m.r. at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRRS    ! Rain water m.r. source
@@ -56,7 +56,7 @@ END MODULE MODI_ICE4_SEDIMENTATION_STAT
 SUBROUTINE ICE4_SEDIMENTATION_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB, KKTE, KKT, KKL, &
                                   &PTSTEP, KRR, OSEDIC, ODEPOSC, PVDEPOSC, PDZZ, &
                                   &PRHODREF, PPABST, PTHT, PRHODJ, &
-				  & PLBDAS, &  ! Modif Wurtz                    
+				  & PLBDAS, & 
                                   &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, &
                                   &PRSS, PRST, PRGS, PRGT,&
                                   &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG, &
@@ -78,6 +78,7 @@ SUBROUTINE ICE4_SEDIMENTATION_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB,
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
 !  P. Wautelet 21/01/2021: initialize untouched part of PFPR
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -104,7 +105,7 @@ REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODREF! Reference den
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PPABST  ! absolute pressure at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PTHT    ! Theta at time t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODJ  ! Dry density * Jacobian
-REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS ! lambda parameter for snow ! Modif Wurtz
+REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS  ! lambda parameter for snow
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRCS    ! Cloud water m.r. source
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRCT    ! Cloud water m.r. at t
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(INOUT)           :: PRRS    ! Rain water m.r. source
@@ -193,7 +194,7 @@ END IF
 CALL INTERNAL_SEDIM_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKTB, KKTE, KKL, &
                         &PRHODREF, PDZZ, ZW, PPABST, PTHT, PTSTEP, &
                         &3, &
-		        &PLBDAS, &!Modif Wurtz
+		        &PLBDAS, &
                         &PRRT, PRRS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -207,7 +208,7 @@ PINPRR(:,:) = ZWSED(:,:,KKB)/XRHOLW                        ! in m/s
 CALL INTERNAL_SEDIM_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKTB, KKTE, KKL, &
                         &PRHODREF, PDZZ, ZW, PPABST, PTHT, PTSTEP, &
                         &4, &
-		        &PLBDAS, &!Modif Wurtz
+		        &PLBDAS, &
                         &PRIT, PRIS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -221,7 +222,7 @@ PINPRI(:,:) = ZWSED(:,:,KKB)/XRHOLW                        ! in m/s
 CALL INTERNAL_SEDIM_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKTB, KKTE, KKL, &
                         &PRHODREF, PDZZ, ZW, PPABST, PTHT, PTSTEP, &
                         &5, &
-		        &PLBDAS, &!Modif Wurtz
+		        &PLBDAS, &
                         &PRST, PRSS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -235,7 +236,7 @@ PINPRS(:,:) = ZWSED(:,:,KKB)/XRHOLW                        ! in m/s
 CALL INTERNAL_SEDIM_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKTB, KKTE, KKL, &
                         &PRHODREF, PDZZ, ZW, PPABST, PTHT, PTSTEP, &
                         &6, &
-		        &PLBDAS, &!Modif Wurtz
+		        &PLBDAS, &
                         &PRGT, PRGS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -250,7 +251,7 @@ IF ( KRR == 7 ) THEN
   CALL INTERNAL_SEDIM_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKTB, KKTE, KKL, &
                           &PRHODREF, PDZZ, ZW, PPABST, PTHT, PTSTEP, &
                           &7, &
-		        &PLBDAS, &!Modif Wurtz
+		        &PLBDAS, &
                           &PRHT, PRHS, ZWSED)
   IF (PRESENT(PFPR)) THEN
     DO JK = KKTB , KKTE
@@ -265,7 +266,7 @@ CONTAINS
   SUBROUTINE INTERNAL_SEDIM_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKTB, KKTE, KKL, &
                                 &PRHODREF, PDZZ, PTSORHODZ, PPABST, PTHT, PTSTEP, &
                                 &KSPE, &
-				&PLBDAS, & !Modif Wurtz
+				&PLBDAS, &
                                 &PRXT, PRXS, PWSED, PSEA, PTOWN)
     !
     !*      0. DECLARATIONS
@@ -284,7 +285,7 @@ CONTAINS
     !
     INTEGER, INTENT(IN) :: KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKT, KKE, KKTB, KKTE, KKL
     REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PRHODREF ! Reference density
-    REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS ! lambda parameter for snow ! Modif Wurtz
+    REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS ! lambda parameter for snow
     REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PDZZ    ! Layer thikness (m)
     REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PTSORHODZ ! TimeStep Over (Rhodref times delta Z)
     REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PPABST
@@ -409,14 +410,14 @@ CONTAINS
 
         IF(PRXT(JI,JJ,JK) > XRTMIN(KSPE)) THEN
             ZWSEDW1(JI,JJ,JK)= XFSEDS *  &
-                              & PRHODREF(JI,JJ,JK)**(-XCEXVT) * & !    Modif Wurtz snow
-                              & (1+(XFVELOS/PLBDAS(JI,JJ,JK))**XALPHAS)**(-XNUS+XEXSEDS/XALPHAS) * & ! GAMMAGEN LH_EXTENDED
+                              & PRHODREF(JI,JJ,JK)**(-XCEXVT) * & 
+                              & (1+(XFVELOS/PLBDAS(JI,JJ,JK))**XALPHAS)**(-XNUS+XEXSEDS/XALPHAS) * & 
 			      & PLBDAS(JI,JJ,JK)**(XBS+XEXSEDS) 
           ENDIF
           IF ( ZQP(JI,JJ) > XRTMIN(KSPE)) THEN
             ZWSEDW2(JI,JJ,JK)= XFSEDS *  &
-                              & PRHODREF(JI,JJ,JK)**(-XCEXVT) * & !    Modif Wurtz snow Thompson
-                              & (1+(XFVELOS/PLBDAS(JI,JJ,JK))**XALPHAS)**(-XNUS+XEXSEDS/XALPHAS) * & ! GAMMAGEN LH_EXTENDED
+                              & PRHODREF(JI,JJ,JK)**(-XCEXVT) * &
+                              & (1+(XFVELOS/PLBDAS(JI,JJ,JK))**XALPHAS)**(-XNUS+XEXSEDS/XALPHAS) * &
                               & PLBDAS(JI,JJ,JK)**(XBS+XEXSEDS) 
 
           ENDIF
@@ -427,9 +428,6 @@ CONTAINS
         IF(KSPE==3) THEN
           ZFSED=XFSEDR
           ZEXSED=XEXSEDR
-      !  ELSEIF(KSPE==5) THEN
-       !   ZFSED=XFSEDS
-        !  ZEXSED=XEXSEDS
         ELSEIF(KSPE==6) THEN
           ZFSED=XFSEDG
           ZEXSED=XEXSEDG
diff --git a/src/MNH/ice4_slow.f90 b/src/MNH/ice4_slow.f90
index df434ef767692cdf4fe4621ffc6e9c12c9e62773..5210967751dbac2c348aecbb013057ae134d62a8 100644
--- a/src/MNH/ice4_slow.f90
+++ b/src/MNH/ice4_slow.f90
@@ -65,6 +65,7 @@ SUBROUTINE ICE4_SLOW(KSIZE, LDSOFT, PCOMPUTE, PRHODREF, PT, &
 !!    MODIFICATIONS
 !!    -------------
 !!
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -173,12 +174,9 @@ IF(LDSOFT) THEN
 ELSE
   PRVDEPS(:) = 0.
   WHERE(ZMASK(:)==1.)
-!    PRVDEPS(:) = ( PSSI(:)/(PRHODREF(:)*PAI(:)) ) *                               &
-!                 ( X0DEPS*PLBDAS(:)**XEX0DEPS + X1DEPS*PCJ(:)*PLBDAS(:)**XEX1DEPS )
-    PRVDEPS(:) = ( PRST(:)*PSSI(:)/PAI(:)) *                               & 	!Modif Wurtz
-                ! ( X0DEPS*PLBDAS(:)**XEX0DEPS + (X1DEPS*PCJ(:)*(PLBDAS(:)+XFVELOS/2.)**(XEX1DEPS)*(PLBDAS(:))**(XNUS+XBS))) ! Thompson
+    PRVDEPS(:) = ( PRST(:)*PSSI(:)/PAI(:)) *                               &
          ( X0DEPS*PLBDAS(:)**XEX0DEPS + (X1DEPS*PCJ(:)*(1+(PLBDAS(:)/(2*XFVELOS)**XALPHAS))**(-XNUS+XEX1DEPS) &
-         *(PLBDAS(:))**(XBS+XEX1DEPS))) ! GAMMAGEN LH_EXTENDED
+         *(PLBDAS(:))**(XBS+XEX1DEPS)))
   END WHERE
 ENDIF
 DO JL=1, KSIZE
@@ -201,15 +199,11 @@ IF(LDSOFT) THEN
 ELSE
   PRIAGGS(:) = 0.
   WHERE(ZMASK(:)==1)
-!    PRIAGGS(:) = XFIAGGS * EXP( XCOLEXIS*(PT(:)-XTT) ) &
-!                         * PRIT(:)                      &
-!                         * PLBDAS(:)**XEXIAGGS          &
-!                         * PRHODREF(:)**(-XCEXVT)
-      PRIAGGS(:) = XFIAGGS * EXP( XCOLEXIS*(PT(:)-XTT) ) &                !Modif Wurtz GAMMAGEN LH_EXTENDED
+      PRIAGGS(:) = XFIAGGS * EXP( XCOLEXIS*(PT(:)-XTT) ) &
                          * PRIT(:)                      &
                          * PRST(:) * (1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXIAGGS/XALPHAS)          &
                          * PRHODREF(:)**(-XCEXVT+1.) &
-                         * ((PLBDAS(:))**(XBS+XEXIAGGS)) ! Thompson
+                         * ((PLBDAS(:))**(XBS+XEXIAGGS))
   END WHERE
 ENDIF
 DO JL=1, KSIZE
diff --git a/src/MNH/ice4_tendencies.f90 b/src/MNH/ice4_tendencies.f90
index 5e5f9000671b65625a576aed5658713cf5b6ff17..f012b218457c4d71c3b38153c6baa048d9d3eaf4 100644
--- a/src/MNH/ice4_tendencies.f90
+++ b/src/MNH/ice4_tendencies.f90
@@ -12,7 +12,6 @@ SUBROUTINE ICE4_TENDENCIES(KSIZE, KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, K
                           &PEXN, PRHODREF, PLVFACT, PLSFACT, K1, K2, K3, &
                           &PPRES, PCF, PSIGMA_RC, &
                           &PCIT, &
-			  &PLBDAS ,& ! Wurtz
                           &PT, PTHT, &
                           &PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PRHT, &
                           &PRVHENI_MR, PRRHONG_MR, PRIMLTC_MR, PRSRIMCG_MR, &
@@ -136,7 +135,6 @@ REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PHLI_LCF
 REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PHLI_HRI
 REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PHLI_LRI
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(OUT)   :: PRAINFR   ! Rain fraction
-REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PLBDAS ! Wurtz
 END SUBROUTINE ICE4_TENDENCIES
 END INTERFACE
 END MODULE MODI_ICE4_TENDENCIES
@@ -147,7 +145,6 @@ SUBROUTINE ICE4_TENDENCIES(KSIZE, KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, K
                           &PEXN, PRHODREF, PLVFACT, PLSFACT, K1, K2, K3, &
                           &PPRES, PCF, PSIGMA_RC, &
                           &PCIT, &
-			  &PLBDAS ,& ! Wurtz
                           &PT, PTHT, &
                           &PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PRHT, &
                           &PRVHENI_MR, PRRHONG_MR, PRIMLTC_MR, PRSRIMCG_MR, &
@@ -177,6 +174,7 @@ SUBROUTINE ICE4_TENDENCIES(KSIZE, KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, K
 !!    -------------
 !
 !  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !
 !
 !*      0. DECLARATIONS
@@ -310,7 +308,6 @@ REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PHLI_LCF
 REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PHLI_HRI
 REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PHLI_LRI
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(OUT)   :: PRAINFR   ! Rain fraction
-REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PLBDAS ! Wurtz
 !
 !*       0.2  declaration of local variables
 !
@@ -321,7 +318,7 @@ REAL, DIMENSION(KSIZE) :: ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT, &
                         & ZRF, &
                         & ZLBDAR, ZLBDAS, ZLBDAG, ZLBDAH, ZLBDAR_RF, &
                         & ZRGSI, ZRGSI_MR
-REAL, DIMENSION(KIT,KJT,KKT) :: ZRRT3D, ZRST3D, ZRGT3D, ZRHT3D
+REAL, DIMENSION(KIT,KJT,KKT) :: ZLBDAS3D, ZRRT3D, ZRST3D, ZRGT3D, ZRHT3D
 INTEGER :: JL
 REAL, DIMENSION(KSIZE) :: ZWETG ! 1. if graupel growths in wet mode, 0. otherwise
 
@@ -404,22 +401,28 @@ ELSE
   !        5.1.6  riming-conversion of the large sized aggregates into graupel (old parametrisation)
   !
   IF(CSNOWRIMING=='OLD ') THEN
-!    ZLBDAS(:)=0.
-!    WHERE(ZRST(:)>0.)
-!      ZLBDAS(:)  = MIN(XLBDAS_MAX, XLBS*(PRHODREF(:)*MAX(ZRST(:), XRTMIN(5)))**XLBEXS)
-!    END WHERE
-    CALL ICE4_RSRIMCG_OLD(KSIZE, ODSOFT, PCOMPUTE==1., &
-                         &PRHODREF, &
-!                         &ZLBDAS, &
-                         &PLBDAS, & ! Wurtz
-                         &ZT, ZRCT, ZRST, &
-                         &PRSRIMCG_MR, PB_RS, PB_RG)
-    DO JL=1, KSIZE
-      ZRST(JL) = ZRST(JL) - PRSRIMCG_MR(JL)
-      ZRGT(JL) = ZRGT(JL) + PRSRIMCG_MR(JL)
-    ENDDO
+     ZLBDAS(:)=0.
+     IF (LSNOW_T) THEN 
+        WHERE (ZRST(:)>XRTMIN(5) .AND. ZT(:)>263.15)
+           ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(:))),XLBDAS_MIN)*XTRANS_MP_GAMMAS
+        END WHERE
+        WHERE (ZRST(:)>XRTMIN(5) .AND. ZT(:)<=263.15)
+           ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(:))),XLBDAS_MIN)*XTRANS_MP_GAMMAS
+        END WHERE
+     ELSE
+        ZLBDAS(:)  = MAX(MIN(XLBDAS_MAX,XLBS*(PRHODREF(:)*MAX(ZRST(:), XRTMIN(5)))**XLBEXS,XLBDAS_MIN)
+     END IF
+     CALL ICE4_RSRIMCG_OLD(KSIZE, ODSOFT, PCOMPUTE==1., &
+                          &PRHODREF, &
+                          &ZLBDAS, &
+                          &ZT, ZRCT, ZRST, &
+                          &PRSRIMCG_MR, PB_RS, PB_RG)
+     DO JL=1, KSIZE
+        ZRST(JL) = ZRST(JL) - PRSRIMCG_MR(JL)
+        ZRGT(JL) = ZRGT(JL) + PRSRIMCG_MR(JL)
+     ENDDO
   ELSE
-    PRSRIMCG_MR(:) = 0.
+     PRSRIMCG_MR(:) = 0.
   ENDIF
 ENDIF
 !
@@ -474,10 +477,17 @@ IF(KSIZE>0) THEN
   !
   !*  compute the slope parameters
   !
-!  ZLBDAS(:)=0.
-!  WHERE(ZRST(:)>0.)
-!    ZLBDAS(:)  = MIN(XLBDAS_MAX, XLBS*(PRHODREF(:)*MAX(ZRST(:), XRTMIN(5)))**XLBEXS)
-!  END WHERE
+  ZLBDAS(:)=0.
+  IF (LSNOW_T) THEN 
+     WHERE (ZRST(:)>XRTMIN(5) .AND. ZT(:)>263.15)
+        ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(:))),XLBDAS_MIN)*XTRANS_MP_GAMMAS
+     END WHERE
+     WHERE (ZRST(:)>XRTMIN(5) .AND. ZT(:)<=263.15)
+        ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(:))),XLBDAS_MIN)*XTRANS_MP_GAMMAS
+     END WHERE
+  ELSE
+     ZLBDAS(:)  = MAX(MIN(XLBDAS_MAX,XLBS*(PRHODREF(:)*MAX(ZRST(:), XRTMIN(5)))**XLBEXS,XLBDAS_MIN)
+  END IF
   ZLBDAG(:)=0.
   WHERE(ZRGT(:)>0.)
     ZLBDAG(:)  = XLBG*(PRHODREF(:)*MAX(ZRGT(:), XRTMIN(6)))**XLBEXG
@@ -508,8 +518,7 @@ ENDIF
 CALL ICE4_SLOW(KSIZE, ODSOFT, PCOMPUTE, PRHODREF, ZT, &
               &PSSI, PLVFACT, PLSFACT, &
               &ZRVT, ZRCT, ZRIT, ZRST, ZRGT, &
-          !    &ZLBDAS, ZLBDAG, &
-              &PLBDAS, ZLBDAG, & ! Wurtz
+              &ZLBDAS, ZLBDAG, &
               &ZAI, ZCJ, PHLI_HCF, PHLI_HRI, &
               &PRCHONI, PRVDEPS, PRIAGGS, PRIAUTS, PRVDEPG, &
               &PA_TH, PA_RV, PA_RC, PA_RI, PA_RS, PA_RG)
@@ -546,8 +555,7 @@ END IF
 CALL ICE4_FAST_RS(KSIZE, ODSOFT, PCOMPUTE, &
                  &PRHODREF, PLVFACT, PLSFACT, PPRES, &
                  &ZDV, ZKA, ZCJ, &
-            !     &ZLBDAR, ZLBDAS, &
-                 &ZLBDAR, PLBDAS, & ! Wurtz
+                 &ZLBDAR, ZLBDAS, &
                  &ZT, ZRVT, ZRCT, ZRRT, ZRST, &
                  &PRIAGGS, &
                  &PRCRIMSS, PRCRIMSG, PRSRIMCG, &
@@ -570,8 +578,7 @@ ENDDO
 CALL ICE4_FAST_RG(KSIZE, ODSOFT, PCOMPUTE, KRR, &
                  &PRHODREF, PLVFACT, PLSFACT, PPRES, &
                  &ZDV, ZKA, ZCJ, PCIT, &
-               !  &ZLBDAR, ZLBDAS, ZLBDAG, &
-                 &ZLBDAR, PLBDAS, ZLBDAG, & ! Wurtz
+                 &ZLBDAR, ZLBDAS, ZLBDAG, &
                  &ZT, ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT, &
                  &ZRGSI, ZRGSI_MR(:), &
                  &ZWETG, &
@@ -590,8 +597,7 @@ IF (KRR==7) THEN
   CALL ICE4_FAST_RH(KSIZE, ODSOFT, PCOMPUTE, ZWETG, &
                    &PRHODREF, PLVFACT, PLSFACT, PPRES, &
                    &ZDV, ZKA, ZCJ, &
-                !   &ZLBDAS, ZLBDAG, ZLBDAR, ZLBDAH, &
-                   &PLBDAS, ZLBDAG, ZLBDAR, ZLBDAH, & ! Wurtz 
+                   &ZLBDAS, ZLBDAG, ZLBDAR, ZLBDAH, &
                    &ZT,  ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT, PRHT, &
                    &PRCWETH, PRIWETH, PRSWETH, PRGWETH, PRRWETH, &
                    &PRCDRYH, PRIDRYH, PRSDRYH, PRRDRYH, PRGDRYH, PRDRYHG, PRHMLTR, &
diff --git a/src/MNH/ini_rain_ice.f90 b/src/MNH/ini_rain_ice.f90
index b2d117da37d7732edd64e794b902bfe437d03da0..5b0b0fc661abb6a7089f893c711ae265e2bf653b 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 with LSNOW_T
 !
 !-------------------------------------------------------------------------------
 !
@@ -297,13 +298,25 @@ XF2I = 0.14
 !
 XAS = 0.02
 XBS = 1.9
-!Values of Thompson (WRF)
-XCS = 12.5 	! Wurtz	GAMMAGEN LH_EXTENDED
-XDS = 0.4	! Wurtz	GAMMAGEN LH_EXTENDED
-XFVELOS = .5 ! Wurtz GAMMAGEN LH_EXTENDED
+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
+XFVELOS = 0.
+END IF
 !
+IF (.NOT. LSNOW_T) THEN
 XCCS = 5.0
 XCXS = 1.0
+END IF
 !
 XF0S = 0.86
 XF1S = 0.28
@@ -379,11 +392,17 @@ XNUR    = 1.0  ! Exponential law
 XALPHAI = 3.0  ! Gamma law for the ice crystal volume
 XNUI    = 3.0  ! Gamma law with little dispersion
 !
-XALPHAS = .214  ! Exponential law
-XNUS    = 43.7  ! Exponential law
-XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
-     ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )    ! Wurtz
-!0.5 * GAMMA(XNUS + 2./XALPHAS) / GAMMA(XNUS+1./XALPHAS) ! Wurtz permet de convertir gamma MP en gamma generalisee
+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
@@ -405,9 +424,13 @@ XLBR   = ( XAR*XCCR*MOMG(XALPHAR,XNUR,XBR) )**(-XLBEXR)
 XLBEXI = 1.0/(-XBI)
 XLBI   = ( XAI*MOMG(XALPHAI,XNUI,XBI) )**(-XLBEXI)
 !
-!XLBEXS = 1.0/(XCXS-XBS) ! Wurtz Not used new snow diagnostic
-!XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS) !Old constant ! OPER
-XLBS   = 1.0/(XAS*MOMG(XALPHAS,XNUS,XBS))  !Wurtz for new diag
+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)
@@ -419,10 +442,9 @@ XLBH   = ( XAH*XCCH*MOMG(XALPHAH,XNUH,XBH) )**(-XLBEXH)
 !
 XLBDAS_MAX = 100000.0
 !
-ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc ! Wurtz verifier par rapport a 1E7
-!XLBDAS_MAX = ( ZCONC_MAX/XCCS )**(1./XCXS) ! OLD CANIAUX
-XLBDAS_MAX = 1.E6 ! Wurtz
-XLBDAS_MIN = 1000. ! Valeurs litterature ! Modif Wurtz ! inutile si on se debrouille bien ?
+ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc
+XLBDAS_MAX = 1.E6
+XLBDAS_MIN = 1000.
 !
 IF (HCLOUD == 'ICE4') THEN
   ALLOCATE( XRTMIN(7) )
@@ -489,15 +511,20 @@ XEXCSEDI =-0.9324*3.0
 WRITE (KLUOUT,FMT=*)' PRISTINE ICE SEDIMENTATION for columns XFSEDI =',XFSEDI
 !
 !
-!XEXSEDS = (XBS+XDS-XCXS)/(XBS-XCXS)	! Wurtz
-!XFSEDS  = XCS*XAS*XCCS*MOMG(XALPHAS,XNUS,XBS+XDS)*                         & ! Wurtz
-!            (XAS*XCCS*MOMG(XALPHAS,XNUS,XBS))**(-XEXSEDS)*(ZRHO00)**XCEXVT  ! Wurtz
-
-!XEXSEDS = -XDS !(2*XBS+XDS) ! Modif Wurtz
-!XEXSEDS = -XDS-XBS-XNUS !(2*XBS+XDS) ! Modif Wurtz Thompson
-XEXSEDS = -XDS-XBS !(2*XBS+XDS) ! Modif Wurtz GAMMAGEN LH_EXTENDED
-XFSEDS  = XCS*MOMG(XALPHAS,XNUS,XBS+XDS)/(MOMG(XALPHAS,XNUS,XBS))    & ! Modif Wurtz
+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)*                         &
@@ -570,17 +597,10 @@ XSCFAC = (0.63**(1./3.))*SQRT((ZRHO00)**XCEXVT) ! One assumes Sc=0.63
 X0DEPI = (4.0*XPI)*XC1I*XF0I*MOMG(XALPHAI,XNUI,1.)
 X2DEPI = (4.0*XPI)*XC1I*XF2I*XC_I*MOMG(XALPHAI,XNUI,XDI+2.0)
 !
-!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)
-
-X0DEPS = XLBS*(4.0*XPI)*XC1S*XF0S*MOMG(XALPHAS,XNUS,1.) ! Wurtz
-X1DEPS = XLBS*(4.0*XPI)*XC1S*XF1S*SQRT(XCS)*MOMG(XALPHAS,XNUS,0.5*XDS+1.5) ! Wurtz
-!XEX0DEPS = XCXS-1.0
-!XEX1DEPS = XCXS-0.5*(XDS+3.0)
-XEX0DEPS = XBS-1.0 ! Wurtz
-!XEX1DEPS = XBS-0.5*(XDS+3.0) ! Wurtz Houze
-!XEX1DEPS = -0.5*(XDS+XNUS+3.0) ! Wurtz Thompson
-XEX1DEPS = -0.5*(XDS+3.0) ! Wurtz GAMMGEN LH_EXTENDED
+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)
@@ -620,11 +640,7 @@ END IF
 !
 XCOLIS   = 0.25 ! Collection efficiency of I+S
 XCOLEXIS = 0.05 ! Temperature factor of the I+S collection efficiency
-!XFIAGGS  = (XPI/4.0)*XCOLIS*XCCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
-XFIAGGS  = XLBS*(XPI/4.0)*XCOLIS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0) ! Wurtz 
-!XEXIAGGS = XCXS-XDS-2.0
-!XEXIAGGS = XBS-XDS-2.0 ! Wurtz Houze
-!XEXIAGGS = -XDS-XNUS-2.0 ! TO MODIFY THOMPSON WURTZ
+XFIAGGS  = XLBS*(XPI/4.0)*XCOLIS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
 XEXIAGGS = -XDS - 2.0 ! GAMMGEN LH_EXTENDED
 !
 GFLAG = .TRUE.
@@ -676,24 +692,15 @@ XEX1EVAR = -1.0-0.5*(XDR+3.0)
 !
 XDCSLIM  = 0.007 ! D_cs^lim = 7 mm as suggested by Farley et al. (1989)
 XCOLCS   = 1.0
-!XEXCRIMSS= XCXS-XDS-2.0 
-!XEXCRIMSS= XBS-XDS-2.0 ! Wurtz Houze
-!XEXCRIMSS= -XDS-XNUS-2.0 ! Modif Wurtz THOMPSON
-XEXCRIMSS= -XDS-2.0 ! Modif Wurtz GAMMAGEN LH_EXTENDED
-!XCRIMSS  = (XPI/4.0)*XCOLCS*XCCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
-XCRIMSS  = XLBS * (XPI/4.0)*XCOLCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0) ! Wurtz
+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  = XLBS*XAS*MOMG(XALPHAS,XNUS,XBS) ! Wurtz
-!XEXSRIMCG= XCXS-XBS
-!XEXSRIMCG= 0. !XBS-XBS Wurtz
-XEXSRIMCG = -XBS ! Wurtz GAMMAGEN LH_EXTENDED
-!XSRIMCG2 = XCCS*XAG*MOMG(XALPHAS,XNUS,XBG)
-XSRIMCG2 = XLBS*XAG*MOMG(XALPHAS,XNUS,XBG) ! Wurtz
+XSRIMCG  = XLBS*XAS*MOMG(XALPHAS,XNUS,XBS)
+XEXSRIMCG = -XBS
+XSRIMCG2 = XLBS*XAG*MOMG(XALPHAS,XNUS,XBG)
 XSRIMCG3 = XFRACM90
-!XEXSRIMCG2=XCXS-XBG
-XEXSRIMCG2=XBS-XBG ! Wurtz
+XEXSRIMCG2=XBS-XBG
 !
 GFLAG = .TRUE.
 IF (GFLAG) THEN
@@ -723,15 +730,13 @@ XRIMINTP2 = 1.0 + XRIMINTP1*LOG( XDCSLIM/(XGAMINC_BOUND_MIN)**(1.0/XALPHAS) )
 !
 !*       7.2    Constants for the accretion of raindrops onto aggregates
 !
-!XFRACCSS = ((XPI**2)/24.0)*XCCS*XCCR*XRHOLW*(ZRHO00**XCEXVT)
-XFRACCSS = XLBS*((XPI**2)/24.0)*XCCR*XRHOLW*(ZRHO00**XCEXVT) ! Wurtz
+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*XCCR*(ZRHO00**XCEXVT)
-XFSACCRG = XLBS*(XPI/4.0)*XAS*XCCR*(ZRHO00**XCEXVT) ! Wurtz
+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.)
@@ -778,15 +783,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, XFVELOS, XCR, XDR,                              & ! Wurtz Thompson
+                 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, XFVELOS, XCR, XDR, 0.,                             & ! Wurtz Thompson
+                 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, XFVELOS, XCR, XDR,                              & ! Wurtz Thompson
+                 ZESR, XBS, XCS, XDS, XFVELOS, XCR, XDR,                     & 
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,  XAG, XBS, XAS                       )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -923,8 +928,7 @@ XCOLSG   = 0.01 ! Collection efficiency of S+G
 XCOLEXSG = 0.1  ! Temperature factor of the S+G collection efficiency
 WRITE (KLUOUT, FMT=*) ' NEW Constants for the aggregate collection by the graupeln'
 WRITE (KLUOUT, 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) ! Wurtz
+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.)
@@ -987,7 +991,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,0., XCS, XDS, XFVELOS,                  & ! Thompson Wurtz
+                ZEGS, XBS, XCG, XDG,0., XCS, XDS, XFVELOS,                  &
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1053,7 +1057,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,0., XCR, XDR, 0.,                              & ! Thompson
+                ZEGR, XBR, XCG, XDG,0., XCR, XDR, 0.,                       &
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1192,7 +1196,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,0., XCS, XDS, XFVELOS,                              & ! Thompson Wurtz
+                ZEHS, XBS, XCH, XDH,0., XCS, XDS, XFVELOS,                  &
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1258,7 +1262,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, 0., XCG, XDG, 0.,                             & ! Thompson
+                ZEHG, XBG, XCH, XDH, 0., XCG, XDG, 0.,                      &
                 XWETLBDAH_MAX, XWETLBDAG_MAX, XWETLBDAH_MIN, XWETLBDAG_MIN, &
                 ZFDINFTY, XKER_GWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1324,7 +1328,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,0., XCR, XDR,0.,                              & ! Thompson
+                ZEHR, XBR, XCH, XDH,0., XCR, XDR,0.,                        &
                 XWETLBDAH_MAX, XWETLBDAR_MAX, XWETLBDAH_MIN, XWETLBDAR_MIN, &
                 ZFDINFTY, XKER_RWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1392,8 +1396,8 @@ IF (GFLAG) THEN
                                                       XAS,XBS
   WRITE(UNIT=KLUOUT,FMT='("                 vitesse: C=",E13.6," D=",E13.6)') &
                                                       XCS,XDS
-  !WRITE(UNIT=KLUOUT,FMT='("           concentration:CC=",E13.6," x=",E13.6)') & ! Modif Wurtz unused
-  !                                                    XCCS,XCXS
+  WRITE(UNIT=KLUOUT,FMT='("           concentration:CC=",E13.6," x=",E13.6)') &
+                                                      XCCS,XCXS
   WRITE(UNIT=KLUOUT,FMT='("            distribution:AL=",E13.6,"NU=",E13.6)') &
                                                       XALPHAS,XNUS
   WRITE(UNIT=KLUOUT,FMT='("            GRAUPEL")')
diff --git a/src/MNH/ini_rain_ice_elec.f90 b/src/MNH/ini_rain_ice_elec.f90
index b99ef0f4c5d47741e5fdde0b44e769ee85a1b281..105c64c68b4c86e70939e5c68506ade9df00c59d 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 with LSNOW_T
 !
 !-------------------------------------------------------------------------------
 !
@@ -271,13 +272,25 @@ XF2I = 0.14
 !
 XAS = 0.02
 XBS = 1.9
-!Values of Thompson (WRF)
-XCS = 12.5 	! Wurtz	GAMMAGEN LH_EXTENDED
-XDS = 0.4	! Wurtz	GAMMAGEN LH_EXTENDED
-XFVELOS = .5 ! Wurtz GAMMAGEN LH_EXTENDED
+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
+XFVELOS = 0.
+END IF
 !
+IF (.NOT. LSNOW_T) THEN
 XCCS = 5.0
 XCXS = 1.0
+END IF
 !
 XF0S = 0.86
 XF1S = 0.28
@@ -344,11 +357,17 @@ XNUR    = 1.0  ! Exponential law
 XALPHAI = 3.0  ! Gamma law for the ice crystal volume
 XNUI    = 3.0  ! Gamma law with little dispersion
 !
-XALPHAS = .214  ! Exponential law
-XNUS    = 43.7  ! Exponential law
-XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
-     ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )    ! Wurtz
-!0.5 * GAMMA(XNUS + 2./XALPHAS) / GAMMA(XNUS+1./XALPHAS) ! Wurtz permet de convertir gamma MP en gamma generalisee
+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
@@ -370,9 +389,13 @@ XLBR   = (XAR * XCCR * MOMG(XALPHAR,XNUR,XBR))**(-XLBEXR)
 XLBEXI = 1.0 / (-XBI)
 XLBI   = (XAI * MOMG(XALPHAI,XNUI,XBI))**(-XLBEXI)
 !
-!XLBEXS = 1.0/(XCXS-XBS) ! Wurtz Not used new snow diagnostic
-!XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS) !Old constant ! OPER
-XLBS   = 1.0/(XAS*MOMG(XALPHAS,XNUS,XBS))  !Wurtz for new diag
+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)
@@ -386,10 +409,9 @@ XLBDAR_MAX = 100000.0
 XLBDAS_MAX = 100000.0
 XLBDAG_MAX = 100000.0
 !
-ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc ! Wurtz verifier par rapport a 1E7
-!XLBDAS_MAX = ( ZCONC_MAX/XCCS )**(1./XCXS) ! OLD CANIAUX
-XLBDAS_MAX = 1.E6 ! Wurtz
-XLBDAS_MIN = 1000. ! Valeurs litterature ! Modif Wurtz ! inutile si on se debrouille bien ?
+ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc
+XLBDAS_MAX = 1.E6
+XLBDAS_MIN = 1000.
 !
 IF (HCLOUD == 'ICE4') THEN
   ALLOCATE( XRTMIN(7) )
@@ -449,17 +471,29 @@ XFSEDI   = 3.89745E11 * MOMG(XALPHAI,XNUI,3.285) *                          &
 XEXCSEDI =-0.9324 * 3.0
 WRITE (KLUOUT,FMT=*)' PRISTINE ICE SEDIMENTATION for columns XFSEDI =',XFSEDI
 !
-XEXSEDS = (XBS + XDS - XCXS) / (XBS - XCXS)
-XFSEDS  = XCS * XAS * XCCS * MOMG(XALPHAS,XNUS,XBS+XDS) *   &
-         (XAS * XCCS * MOMG(XALPHAS,XNUS,XBS))**(-XEXSEDS) * (ZRHO00)**XCEXVT
 !
-XEXSEDG = (XBG + XDG - XCXG) / (XBG - XCXG)
-XFSEDG  = XCG * XAG * XCCG * MOMG(XALPHAG,XNUG,XBG+XDG) *   &
-         (XAG * XCCG * MOMG(XALPHAG,XNUG,XBG))**(-XEXSEDG) * (ZRHO00)**XCEXVT
+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)*                         &
+            (XAG*XCCG*MOMG(XALPHAG,XNUG,XBG))**(-XEXSEDG)*(ZRHO00)**XCEXVT
 !
-XEXSEDH = (XBH + XDH - XCXH) / (XBH - XCXH)
-XFSEDH  = XCH * XAH * XCCH * MOMG(XALPHAH,XNUH,XBH+XDH) *   &
-         (XAH * XCCH * MOMG(XALPHAH,XNUH,XBH))**(-XEXSEDH) * (ZRHO00)**XCEXVT
+XEXSEDH = (XBH+XDH-XCXH)/(XBH-XCXH)
+XFSEDH  = XCH*XAH*XCCH*MOMG(XALPHAH,XNUH,XBH+XDH)*                         &
+            (XAH*XCCH*MOMG(XALPHAH,XNUH,XBH))**(-XEXSEDH)*(ZRHO00)**XCEXVT
 !
 !
 !-------------------------------------------------------------------------------
@@ -524,10 +558,10 @@ XSCFAC = (0.63**(1./3.)) * SQRT((ZRHO00)**XCEXVT) ! One assumes Sc=0.63
 X0DEPI = (4.0 * XPI) * XC1I * XF0I * MOMG(XALPHAI,XNUI,1.)
 X2DEPI = (4.0 * XPI) * XC1I * XF2I * XC_I * MOMG(XALPHAI,XNUI,XDI+2.0)
 !
-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)
@@ -560,9 +594,8 @@ END IF
 !
 XCOLIS   = 0.25 ! Collection efficiency of I+S
 XCOLEXIS = 0.05 ! Temperature factor of the I+S collection efficiency
-XFIAGGS  = (XPI / 4.0) * XCOLIS * XCCS * XCS * (ZRHO00**XCEXVT) * &
-            MOMG(XALPHAS,XNUS,XDS+2.0)
-XEXIAGGS = XCXS - XDS - 2.0
+XFIAGGS  = XLBS*(XPI/4.0)*XCOLIS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
+XEXIAGGS = -XDS - 2.0 ! GAMMGEN LH_EXTENDED
 !
 GFLAG = .TRUE.
 IF (GFLAG) THEN
@@ -612,12 +645,15 @@ XEX1EVAR = -1.0 - 0.5 * (XDR + 3.0)
 !
 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)
-XEXSRIMCG= XCXS - XBS
+XSRIMCG  = XLBS*XAS*MOMG(XALPHAS,XNUS,XBS)
+XEXSRIMCG = -XBS
+XSRIMCG2 = XLBS*XAG*MOMG(XALPHAS,XNUS,XBG)
+XSRIMCG3 = XFRACM90
+XEXSRIMCG2=XBS-XBG
 !
 GFLAG = .TRUE.
 IF (GFLAG) THEN
@@ -648,13 +684,13 @@ XRIMINTP2 = 1.0 + XRIMINTP1 * LOG(XDCSLIM/(XGAMINC_BOUND_MIN)**(1.0/XALPHAS))
 !
 !*       7.2    Constants for the accretion of raindrops onto aggregates
 !
-XFRACCSS = ((XPI**2) / 24.0) * XCCS * XCCR * 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 * XCCR * (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.)
@@ -702,15 +738,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, XFVELOS, XCR, XDR,                              & ! Wurtz Thompson
+                 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, XFVELOS, XCR, XDR, 0.,                             & ! Wurtz Thompson
+                 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, XFVELOS, XCR, XDR,                              & ! Wurtz Thompson
+                 ZESR, XBS, XCS, XDS, XFVELOS, XCR, XDR,                     & 
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,  XAG, XBS, XAS                       )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -828,9 +864,9 @@ XCOLIG   = 0.01 ! Collection efficiency of I+G
 XCOLEXIG = 0.1  ! Temperature factor of the I+G collection efficiency
 WRITE (KLUOUT, FMT=*) ' NEW Constants for the cloud ice collection by the graupeln'
 WRITE (KLUOUT, FMT=*) ' XCOLIG, XCOLEXIG  = ',XCOLIG,XCOLEXIG
-!
-XFIDRYG = (XPI / 4.0) * XCOLIG * XCCG * XCG * (ZRHO00**XCEXVT) * &
-           MOMG(XALPHAG,XNUG,XDG+2.0)
+XFIDRYG = (XPI/4.0)*XCOLIG*XCCG*XCG*(ZRHO00**XCEXVT)*MOMG(XALPHAG,XNUG,XDG+2.0)
+XEXFIDRYG=(XCXG-XDG-2.)/(XCXG-XBG)
+XFIDRYG2=XFIDRYG/XCOLIG*(XAG*XCCG*MOMG(XALPHAG,XNUG,XBG))**(-XEXFIDRYG)
 !
 GFLAG = .TRUE.
 IF (GFLAG) THEN
@@ -847,8 +883,7 @@ XCOLSG   = 0.01 ! Collection efficiency of S+G
 XCOLEXSG = 0.1  ! Temperature factor of the S+G collection efficiency
 WRITE (KLUOUT, FMT=*) ' NEW Constants for the aggregate collection by the graupeln'
 WRITE (KLUOUT, 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.)
@@ -913,7 +948,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,0., XCS, XDS, XFVELOS,                  & ! Thompson Wurtz
+                ZEGS, XBS, XCG, XDG,0., XCS, XDS, XFVELOS,                  &
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -979,7 +1014,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,0., XCR, XDR, 0.,                              & ! Thompson
+                ZEGR, XBR, XCG, XDG,0., XCR, XDR, 0.,                       &
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1041,7 +1076,10 @@ 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)
+XCOLSH   = 0.01 ! Collection efficiency of S+H
+XCOLEXSH = 0.1  ! Temperature factor of the S+H collection efficiency
+!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.)
@@ -1097,7 +1135,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,0., XCS, XDS, XFVELOS,                              & ! Thompson Wurtz
+                ZEHS, XBS, XCH, XDH,0., XCS, XDS, XFVELOS,                  &
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1163,7 +1201,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, 0., XCG, XDG, 0.,                             & ! Thompson
+                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/modd_rain_ice_descr.f90 b/src/MNH/modd_rain_ice_descr.f90
index 96295e4e7841c57028903c03f50731d03b9a9c6c..e754d5a8fac2aad39b74863d924df9df760a27a6 100644
--- a/src/MNH/modd_rain_ice_descr.f90
+++ b/src/MNH/modd_rain_ice_descr.f90
@@ -76,7 +76,11 @@ REAL,SAVE :: XALPHAS,XNUS,XLBEXS,XLBS ! Snow/agg.      distribution parameters
 REAL,SAVE :: XALPHAG,XNUG,XLBEXG,XLBG ! Graupel        distribution parameters 
 REAL,SAVE :: XALPHAH,XNUH,XLBEXH,XLBH ! Hail           distribution parameters
 !
-REAL,SAVE :: XLBDAR_MAX,XLBDAS_MAX,XLBDAG_MAX ! Max values allowed for the shape
+LOGICAL,SAVE :: LSNOW_T         ! Snow parameterization from Wurtz (2021)
+REAL,SAVE :: XFVELOS            ! factor for snow fall speed after Thompson (2008)
+REAL,SAVE :: XTRANS_MP_GAMMAS	! coefficient to convert lambdas for gamma function
+!
+REAL,SAVE :: XLBDAR_MAX,XLBDAS_MIN,XLBDAS_MAX,XLBDAG_MAX ! Max values allowed for the shape
                                               ! parameters (rain,snow,graupeln)
 !
 REAL,DIMENSION(:),SAVE,ALLOCATABLE :: XRTMIN ! Min values allowed for the mixing ratios
diff --git a/src/MNH/modd_rain_ice_param.f90 b/src/MNH/modd_rain_ice_param.f90
index 4394deea2c14a1e1bd0c619220d58ff5c0e57559..434c4bc761d351bc6756fbdd6b5026d1a38141a3 100644
--- a/src/MNH/modd_rain_ice_param.f90
+++ b/src/MNH/modd_rain_ice_param.f90
@@ -39,9 +39,6 @@
 !             ------------
 !
 IMPLICIT NONE
-!certaines constantes a deplacer dans modd_rain_ice_descr.F90 ... peut etre ...
-REAL,SAVE :: XFVELOS,XLBDAS_MIN ! Wurtz
-REAL,SAVE :: XTRANS_MP_GAMMAS			 ! coefficient to convert lambdas for gamma function Wurtz
 REAL,DIMENSION(2),SAVE :: XFSEDC                 ! Constants for sedimentation fluxes of C
 REAL,SAVE :: XFSEDR,XEXSEDR,                   & ! Constants for sedimentation
              XFSEDI,XEXCSEDI,XEXRSEDI,         & ! fluxes of R, I, S and G
diff --git a/src/MNH/rain_ice_red.f90 b/src/MNH/rain_ice_red.f90
index 153617ef6f9135c80f55dede70ccfba54f1a9b11..a44c3af43c419dd7b329d586e7dafffc1b1cb709 100644
--- a/src/MNH/rain_ice_red.f90
+++ b/src/MNH/rain_ice_red.f90
@@ -257,6 +257,7 @@ END MODULE MODI_RAIN_ICE_RED
 !  P. Wautelet 17/01/2020: move Quicksort to tools.f90
 !  P. Wautelet    02/2020: use the new data structures and subroutines for budgets
 !  P. Wautelet 25/02/2020: bugfix: add missing budget: WETH_BU_RRG
+!  J. Wurtz       03/2022: New snow characteristics with LSNOW_T
 !-----------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -273,8 +274,7 @@ USE MODD_CST,            ONLY: XCI,XCL,XCPD,XCPV,XLSTT,XLVTT,XTT
 USE MODD_PARAMETERS,     ONLY: JPVEXT,XUNDEF
 USE MODD_PARAM_ICE,      ONLY: CSUBG_PR_PDF,CSUBG_RC_RR_ACCR,CSUBG_RR_EVAP,LDEPOSC,LFEEDBACKT,LSEDIM_AFTER, &
                                NMAXITER,XMRSTEP,XTSTEP_TS,XVDEPOSC
-USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN,XLBDAS_MAX
-USE MODD_RAIN_ICE_PARAM, ONLY: XLBDAS_MIN,XTRANS_MP_GAMMAS
+USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN,XLBDAS_MIN,XLBDAS_MAX,XTRANS_MP_GAMMAS,LSNOW_T
 USE MODD_VAR_ll,         ONLY: IP
 
 use mode_budget,         only: Budget_store_add, Budget_store_init, Budget_store_end
@@ -357,7 +357,7 @@ REAL, DIMENSION(KIT,KJT,KKT,KRR), OPTIONAL, INTENT(OUT)  :: PFPR ! upper-air pre
 !
 !*       0.2   Declarations of local variables :
 !
-REAL,    DIMENSION(SIZE(PRST,1),SIZE(PRST,2),SIZE(PRST,3)) :: PLBDAS ! Modif  !lbda parameter snow
+REAL,    DIMENSION(SIZE(PRST,1),SIZE(PRST,2),SIZE(PRST,3)) :: ZLBDAS ! Modif  !lbda parameter snow
 
 INTEGER :: IIB           !  Define the domain where is
 INTEGER :: IIE           !  the microphysical sources have to be computed
@@ -402,7 +402,6 @@ REAL, DIMENSION(KSIZE) :: ZRVT,     & ! Water vapor m.r. at t
                         & ZRRT,     & ! Rain water m.r. at t
                         & ZRIT,     & ! Pristine ice m.r. at t
                         & ZRST,     & ! Snow/aggregate m.r. at t
-                        & ZLBDAS,   & ! Snow/aggregate lbdas ! Modif Wurtz
                         & ZRGT,     & ! Graupel m.r. at t
                         & ZRHT,     & ! Hail m.r. at t
                         & ZCIT,     & ! Pristine ice conc. at t
@@ -576,22 +575,25 @@ ELSE
   ENDDO
 ENDIF
 
-!Compute lambda_snow parameter ! Modif Wurtz
+!Compute lambda_snow parameter
 !ZT en KELVIN
-PLBDAS(:,:,:)=1000. !Wurtz
+ZLBDAS(:,:,:)=1000.
 DO JK = 1, KKT
-  DO JJ = 1, KJT
-    DO JI = 1, KIT
-	IF (PRST(JI,JJ,JK)>XRTMIN(5)) THEN
-
-			IF(ZT(JI,JJ,JK)>263.15) THEN
-				PLBDAS(JI,JJ,JK) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(JI,JJ,JK))),XLBDAS_MIN)*XTRANS_MP_GAMMAS ! On le calcul ici car il est utilise dans la sedim en 3D et dans tendencies
-			ELSE
-				PLBDAS(JI,JJ,JK) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(JI,JJ,JK))),XLBDAS_MIN)*XTRANS_MP_GAMMAS
-			END IF
-			END IF
-		END DO
-	END DO
+   DO JJ = 1, KJT
+      DO JI = 1, KIT
+         IF (LSNOW_T) THEN 
+            IF (PRST(JI,JJ,JK)>XRTMIN(5)) THEN
+               IF(ZT(JI,JJ,JK)>263.15) THEN
+                  ZLBDAS(JI,JJ,JK) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(JI,JJ,JK))),XLBDAS_MIN)*XTRANS_MP_GAMMAS
+               ELSE
+                  ZLBDAS(JI,JJ,JK) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(JI,JJ,JK))),XLBDAS_MIN)*XTRANS_MP_GAMMAS
+               END IF
+            END IF
+         ELSE
+            ZLBDAS(JI,JJ,JK)  = MAX(MIN(XLBDAS_MAX,XLBS*(PRHODREF(JI,JK,JL)*MAX(PRST(JI,JK,JL), XRTMIN(5)))**XLBEXS,XLBDAS_MIN)
+         END IF
+      END DO
+   END DO
 END DO
 !
 !-------------------------------------------------------------------------------
@@ -622,7 +624,7 @@ IF(.NOT. LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                   &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, &
                                   &PRHODREF, PPABST, PTHT, PRHODJ, &
-				  &PLBDAS, &  ! Modif Wurtz                                  
+				  &ZLBDAS, &
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -632,7 +634,7 @@ IF(.NOT. LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                   &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, &
                                   &PRHODREF, PPABST, PTHT, PRHODJ, &
-				  &PLBDAS, &  ! Modif Wurtz                                 
+				  &ZLBDAS, &
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -647,7 +649,7 @@ IF(.NOT. LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                    &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, &
                                    &PRHODREF, PPABST, PTHT, PRHODJ, &
-				   &PLBDAS, &  ! Modif Wurtz                                 
+				   &ZLBDAS, &
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &
@@ -656,7 +658,7 @@ IF(.NOT. LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                    &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, &
                                    &PRHODREF, PPABST, PTHT, PRHODJ, &
-				   &PLBDAS, &  ! Modif Wurtz                                 
+				   &ZLBDAS, &
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &
@@ -736,7 +738,6 @@ IF(IMICRO>0) THEN
     ELSE
       ZHLI_LCF(JL)=0.
     ENDIF
-    ZLBDAS(JL) = PLBDAS(I1(JL),I2(JL),I3(JL))
   ENDDO
   IF(GEXT_TEND) THEN
     DO JL=1, IMICRO
@@ -886,7 +887,6 @@ DO WHILE(ANY(ZTIME(:)<PTSTEP)) ! Loop to *really* compute tendencies
                         &ZEXN, ZRHODREF, ZLVFACT, ZLSFACT, I1, I2, I3, &
                         &ZPRES, ZCF, ZSIGMA_RC,&
                         &ZCIT, &
-			&ZLBDAS, & ! Wurtz
                         &ZZT, ZTHT, &
                         &ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT, ZRHT, &
                         &ZRVHENI_MR, ZRRHONG_MR, ZRIMLTC_MR, ZRSRIMCG_MR, &
@@ -1712,7 +1712,7 @@ IF(LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                   &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, &
                                   &PRHODREF, PPABST, PTHT, PRHODJ, &
-                                  &PLBDAS, &   ! Modif Wurtz				   
+                                  &ZLBDAS, &				   
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -1722,7 +1722,7 @@ IF(LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                   &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ,&
                                   &PRHODREF, PPABST, PTHT, PRHODJ, &
-                                  &PLBDAS, &   ! Modif Wurtz				   
+                                  &ZLBDAS, &				   
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -1737,7 +1737,7 @@ IF(LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                    &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, &
                                    &PRHODREF, PPABST, PTHT, PRHODJ, &
-                                   &PLBDAS, &   ! Modif Wurtz				   
+                                   &ZLBDAS, &			   
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &
@@ -1746,7 +1746,7 @@ IF(LSEDIM_AFTER) THEN
       CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, KIT, IJB, IJE, KJT, IKB, IKE, IKTB, IKTE, KKT, KKL, &
                                    &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, &
                                    &PRHODREF, PPABST, PTHT, PRHODJ, &
-                                   &PLBDAS, &   ! Modif Wurtz				   
+                                   &ZLBDAS, &			   
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &