diff --git a/src/MNH/ice4_fast_rg.f90 b/src/MNH/ice4_fast_rg.f90
index b84dda857e7a679561f062eabbe0570f7fc22408..c78a154bb43fe2494c7a86baf99a7a0df7676e67 100644
--- a/src/MNH/ice4_fast_rg.f90
+++ b/src/MNH/ice4_fast_rg.f90
@@ -335,12 +335,19 @@ ELSE
     END DO
     !
     WHERE(GDRY(:))
-      PRG_TEND(:, IRSWETG)=XFSDRYG*ZZW(:)                         & ! RSDRYG
+!      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
                                     / XCOLSG &
-                  *(PLBDAS(:)**(XCXS-XBS))*( PLBDAG(:)**XCXG )    &
-                  *(PRHODREF(:)**(-XCEXVT-1.))                    &
-                       *( XLBSDRYG1/( PLBDAG(:)**2              ) + &
-                          XLBSDRYG2/( PLBDAG(:)   * PLBDAS(:)   ) + &
+                  *(PRST(:))*( PLBDAG(:)**XCXG )    &
+                  *(PRHODREF(:)**(-XCEXVT))                    &
+                       *( XLBSDRYG1/( PLBDAG(:)**2              ) + &	! Il s'agit de moments (?)
+                          XLBSDRYG2/( PLBDAG(:)   *  PLBDAS(:)   ) + &
                           XLBSDRYG3/(               PLBDAS(:)**2))
       PRG_TEND(:, IRSDRYG)=PRG_TEND(:, IRSWETG)*XCOLSG*EXP(XCOLEXSG*(PT(:)-XTT))
     END WHERE
diff --git a/src/MNH/ice4_fast_rh.f90 b/src/MNH/ice4_fast_rh.f90
index fcac937485414ba29fd691cb0774a32cb3ea4a3c..98106308dc0fd54a23ee2ce45f3e174e36f66f78 100644
--- a/src/MNH/ice4_fast_rh.f90
+++ b/src/MNH/ice4_fast_rh.f90
@@ -269,10 +269,16 @@ ELSE
     END DO
     !
     WHERE(GWET(:))
-      PRH_TEND(:, IRSWETH)=XFSWETH*ZZW(:)                       & ! RSWETH
-                    *( PLBDAS(:)**(XCXS-XBS) )*( PLBDAH(:)**XCXH )  &
-                       *( PRHODREF(:)**(-XCEXVT-1.) )               &
-                       *( XLBSWETH1/( PLBDAH(:)**2              ) + &
+!      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
+                    *( PRST(:))*( PLBDAH(:)**XCXH )  &
+                       *( PRHODREF(:)**(-XCEXVT) )               &
+                       *( XLBSWETH1/( PLBDAH(:)**2              ) + &	! Il s'agit de moment (?)
                           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 6d71c7b61b8188969aa488a3b22d65ad15d7cc26..15084073ac3aa6aeb1427dad2e1d3d0adf756a3f 100644
--- a/src/MNH/ice4_fast_rs.f90
+++ b/src/MNH/ice4_fast_rs.f90
@@ -170,8 +170,8 @@ ELSE
     PRS_TEND(:, IFREEZ1)=PKA(:)*(XTT-PT(:)) +                              &
              (PDV(:)*(XLVTT+(XCPV-XCL)*(PT(:)-XTT)) &
                            *(XESTT-PRS_TEND(:, IFREEZ1))/(XRV*PT(:))           )
-    PRS_TEND(:, IFREEZ1)=PRS_TEND(:, IFREEZ1)* ( X0DEPS*       PLBDAS(:)**XEX0DEPS +     &
-                           X1DEPS*PCJ(:)*PLBDAS(:)**XEX1DEPS )/ &
+    PRS_TEND(:, IFREEZ1)=PRS_TEND(:, IFREEZ1)* PRST(:) * ( X0DEPS*       PLBDAS(:)**XEX0DEPS +     &
+                           X1DEPS*PCJ(:)*PLBDAS(:) **(XBS+XEX1DEPS)*(1+(XFVELOS/(2*PLBDAS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS))/ &
                           ( PRHODREF(:)*(XLMTT-XCL*(XTT-PT(:))) )
     PRS_TEND(:, IFREEZ2)=(PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PT(:)))   ) / &
                           ( PRHODREF(:)*(XLMTT-XCL*(XTT-PT(:))) )
@@ -218,7 +218,7 @@ ELSE
     !        5.1.1  select the PLBDAS
     !
     DO JJ = 1, IGRIM
-      ZVEC1(JJ) = PLBDAS(I1(JJ))
+      ZVEC1(JJ) = PLBDAS(I1(JJ))**XALPHAS + XFVELOS**XALPHAS)**(1./XALPHAS)
     END DO
     !
     !        5.1.2  find the next lower indice for the PLBDAS in the geometrical
@@ -243,9 +243,13 @@ 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
+!                                      *   PLBDAS(:)**XEXCRIMSS &
+!                                      * PRHODREF(:)**(-XCEXVT)
+      PRS_TEND(:, IRCRIMSS) = XCRIMSS * ZZW(:) * PRCT(:)                & ! RCRIMSS	 !Wurtz ! Thompson
+                                      * PRST(:)*(1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
+                                      * PRHODREF(:)**(-XCEXVT+1.) &
+				      * (PLBDAS(:)) ** (XEXCRIMSS+XBS) ! Thompson
     END WHERE
     !
     !        5.1.5  perform the linear interpolation of the normalized
@@ -270,19 +274,25 @@ ELSE
     !
     !
     WHERE(GRIM(:))
-      PRS_TEND(:, IRCRIMS)=XCRIMSG * PRCT(:)               & ! RCRIMS
-                                   * PLBDAS(:)**XEXCRIMSG  &
-                                   * PRHODREF(:)**(-XCEXVT)
+!      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
       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 * PLBDAS(:)**XEXSRIMCG*(1.0-ZZW(:))
+        PRS_TEND(:, IRSRIMCG)=XSRIMCG * PRST(:)*PRHODREF(:)*PLBDAS(:)**(XEXSRIMCG+XBS)*(1.0-ZZW(:)) !Wurtz
+
         PRS_TEND(:, IRSRIMCG)=ZZW6(:)*PRS_TEND(:, IRSRIMCG)/ &
                        MAX(1.E-20, &
-                           XSRIMCG3*XSRIMCG2*PLBDAS(:)**XEXSRIMCG2*(1.-ZZW2(:)) - &
+                           XSRIMCG3*XSRIMCG2*PRST(:)*PRHODREF(:)*PLBDAS(:)**XEXSRIMCG2*(1.-ZZW2(:)) - &
                            XSRIMCG3*PRS_TEND(:, IRSRIMCG))
       END WHERE
     ELSE
@@ -383,8 +393,13 @@ 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*( PLBDAS(:)**XCXS )*( PRHODREF(:)**(-XCEXVT-1.) ) &
+            XFRACCSS*( PRST(:)*PLBDAS(:)**XBS )*( PRHODREF(:)**(-XCEXVT) ) & ! Wurtz
        *( XLBRACCS1/((PLBDAS(:)**2)               ) +                  &
           XLBRACCS2/( PLBDAS(:)    * PLBDAR(:)    ) +                  &
           XLBRACCS3/(               (PLBDAR(:)**2)) )/PLBDAR(:)**4
@@ -429,8 +444,13 @@ ELSE
     !               into graupeln
     !
     WHERE(GACC(:))
-      PRS_TEND(:, IRSACCRG) = XFSACCRG*ZZW(:)*                    & ! RSACCRG
-          ( PLBDAS(:)**(XCXS-XBS) )*( PRHODREF(:)**(-XCEXVT-1.) ) &
+!      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
+          ( PRST(:))*( PRHODREF(:)**(-XCEXVT) ) &
          *( XLBSACCR1/((PLBDAR(:)**2)               ) +           &
             XLBSACCR2/( PLBDAR(:)    * PLBDAS(:)    ) +           &
             XLBSACCR3/(               (PLBDAS(:)**2)) )/PLBDAR(:)
@@ -495,12 +515,23 @@ 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(:) *             &
+!                         ( 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
+         PRST(:)*PRHODREF(:) *    &
+         ( X0DEPS       *PLBDAS(:)**XEX0DEPS +     &
+         X1DEPS*PCJ(:)*(1+(XFVELOS/(2*PLBDAS(:))**XALPHAS))**(XNUS+XEX1DEPS/XALPHAS)*((PLBDAS(:))**(XBS+XEX1DEPS))) -   &
+         ( PRS_TEND(:, IRCRIMS) + PRS_TEND(:, IRRACCS)) *       &
+         ( PRHODREF(:)*XCL*(XTT-PT(:))) ) /    &
+         ( PRHODREF(:)*XLMTT ) )
+    !
+    ! note that RSCVMG = RSMLT*XFSCVMG but no heat is exchanged (at the rate RSMLT)
+    ! because the graupeln produced by this process are still icy!!!
+    !
     ! When T < XTT, rc is collected by snow (riming) to produce snow and graupel
     ! When T > XTT, if riming was still enabled, rc would produce snow and graupel with snow becomming graupel (conversion/melting) and graupel becomming rain (melting)
     ! To insure consistency when crossing T=XTT, rc collected with T>XTT must be transformed in rain.
diff --git a/src/MNH/ice4_rsrimcg_old.f90 b/src/MNH/ice4_rsrimcg_old.f90
index cf88792b1b42827bdce381c4e2ad5644b3fed376..5fc932536a176f4deaf51df66fc320b06f6afa18 100644
--- a/src/MNH/ice4_rsrimcg_old.f90
+++ b/src/MNH/ice4_rsrimcg_old.f90
@@ -131,8 +131,10 @@ IF(.NOT. ODSOFT) THEN
     !
     !
     WHERE(GRIM(:))
-      PRSRIMCG_MR(:) = XSRIMCG * PLBDAS(:)**XEXSRIMCG   & ! RSRIMCG
-                               * (1.0 - ZZW(:) )/PRHODREF(:)
+!      PRSRIMCG_MR(:) = XSRIMCG * PLBDAS(:)**XEXSRIMCG   & ! RSRIMCG
+!                               * (1.0 - ZZW(:) )/PRHODREF(:)
+      PRSRIMCG_MR(:) = XSRIMCG * PLBDAS(:)**XEXSRIMCG   & ! RSRIMCG	!Modif Wurtz concentration snow
+                               * (1.0 - ZZW(:) )*PRST(:)
       PRSRIMCG_MR(:)=MIN(PRST(:), PRSRIMCG_MR(:))
     END WHERE
   END IF
diff --git a/src/MNH/ice4_sedimentation_split.f90 b/src/MNH/ice4_sedimentation_split.f90
index cb0a147d070b865e8e9391bacf7ef143c3727311..ae003910312484142c3a255be712cc3fa1f2f1a2 100644
--- a/src/MNH/ice4_sedimentation_split.f90
+++ b/src/MNH/ice4_sedimentation_split.f90
@@ -8,6 +8,7 @@ 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,  &
@@ -25,6 +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(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,6 +55,7 @@ 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,  &
@@ -102,6 +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(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
@@ -239,6 +243,7 @@ IF (GSEDIC) THEN
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &2, &
+			  &PLBDAS, & ! Modif Wurtz
                           &ZRCT, PRCS, PINPRC, ZPRCS, &
                           &ZRAY, ZLBC, ZFSEDC, ZCONC3D, PFPR=PFPR)
 ENDIF
@@ -265,8 +270,9 @@ END IF
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &3, &
+			  &PLBDAS, & ! Modif Wurtz
                           &ZRRT, PRRS, PINPRR, ZPRRS, &
-                          PFPR=PFPR)
+                          &PFPR=PFPR)
 !
 !*       2.3   for pristine ice
 !
@@ -274,6 +280,7 @@ END IF
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &4, &
+			  &PLBDAS, & ! Modif Wurtz
                           &ZRIT, PRIS, PINPRI, ZPRIS, &
                           PFPR=PFPR)
 !
@@ -283,6 +290,7 @@ END IF
                         &XSPLIT_MAXCFL, &
                         &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                         &5, &
+                        &PLBDAS, & ! Modif Wurtz
                         &ZRST, PRSS, PINPRS, ZPRSS, &
                         PFPR=PFPR)
 !
@@ -292,6 +300,7 @@ END IF
                         &XSPLIT_MAXCFL, &
                         &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                         &6, &
+                        &PLBDAS, & ! Modif Wurtz
                         &ZRGT, PRGS, PINPRG, ZPRGS, &
                         PFPR=PFPR)
 !
@@ -302,6 +311,7 @@ IF (IRR==7) THEN
                           &XSPLIT_MAXCFL, &
                           &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
                           &7, &
+			  &PLBDAS, & ! Modif Wurtz
                           &ZRHT, PRHS, PINPRH, ZPRHS, &
                           PFPR=PFPR)
 ENDIF
@@ -315,7 +325,7 @@ CONTAINS
 !
 SUBROUTINE INTERNAL_SEDIM_SPLI(KIB,KIE,KIT,KJB,KJE,KJT,KKB,KKTB,KKTE,KKT,KKL,KRR, &
                               &PMAXCFL,PRHODREF,POORHODZ,PDZZ,PPABST,PTHT,PTSTEP, &
-                              &KSPE,PRXT,PRXS,PINPRX,PPRXS,                       &
+                              &KSPE,PLBDAS,PRXT,PRXS,PINPRX,PPRXS,                &
                               &PRAY,PLBC,PFSEDC,PCONC3D,PFPR)
 !
 !*      0. DECLARATIONS
@@ -344,6 +354,7 @@ REAL, DIMENSION(KIT,KJT),     INTENT(OUT)             :: PINPRX ! instant precip
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PPRXS ! external tendencie
 REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN), OPTIONAL    :: PRAY, PLBC, PFSEDC, PCONC3D
 REAL, DIMENSION(KIT,KJT,KKT,KRR), INTENT(INOUT), OPTIONAL :: PFPR    ! upper-air precipitation fluxes
+    REAL, DIMENSION(KIT,KJT,KKT), INTENT(IN)              :: PLBDAS ! lambda parameter for snow ! Modif Wurtz
 !
 !*       0.2  declaration of local variables
 !
@@ -433,15 +444,33 @@ DO WHILE (ANY(ZREMAINT>0.))
                             &      ALOG(PRHODREF(JI,JJ,JK)*PRXT(JI,JJ,JK)) )**XEXCSEDI
       ENDIF
     ENDDO
+    ELSEIF(KSPE==5) THEN
+      ! ******* for snow
+      PWSED(:,:,:) = 0.
+      DO JL=1, KSEDIM
+        JI=I1(JL)
+        JJ=I2(JL)
+        JK=I3(JL)
+        IF(PRXT(JI,JJ,JK)> XRTMIN(KSPE)) THEN
+
+        PWSED(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
+			      & PLBDAS(JI, JJ, JK) ** (XBS+XEXSEDS)
+ ! GAMMAGEN_LH_EXTENDED
+
+        ENDIF
+      ENDDO
   ELSE
     ! ******* for other species
     SELECT CASE(KSPE)
       CASE(3)
         ZFSED=XFSEDR
         ZEXSED=XEXSEDR
-      CASE(5)
-        ZFSED=XFSEDS
-        ZEXSED=XEXSEDS
+!      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 3cbb31493eac8295e718f2e1438e4e3a269520e7..4daa043dab99d656ba93cd941120e3e35c81619d 100644
--- a/src/MNH/ice4_sedimentation_stat.f90
+++ b/src/MNH/ice4_sedimentation_stat.f90
@@ -8,6 +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                    
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT,&
                                    &PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG, &
@@ -26,6 +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(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
@@ -54,6 +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                    
                                   &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, &
                                   &PRSS, PRST, PRGS, PRGT,&
                                   &PINPRC, PINDEP, PINPRR, PINPRI, PINPRS, PINPRG, &
@@ -83,7 +86,8 @@ SUBROUTINE ICE4_SEDIMENTATION_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKTB,
 USE MODD_CST
 
 USE MODE_MSG
-
+USE MODD_RAIN_ICE_DESCR
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -100,6 +104,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(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
@@ -160,6 +165,7 @@ IF (OSEDIC) THEN
   CALL INTERNAL_SEDIM_STAT(KIB, KIE, KIT, KJB, KJE, KJT, KKB, KKE, KKT, KKTB, KKTE, KKL, &
                           &PRHODREF, PDZZ, ZW, PPABST, PTHT, PTSTEP, &
                           &2, &
+			  &PLBDAS, &!Modif Wurtz
                           &PRCT, PRCS, ZWSED, PSEA, PTOWN)
   IF (PRESENT(PFPR)) THEN
     DO JK = KKTB , KKTE
@@ -187,6 +193,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
                         &PRRT, PRRS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -200,6 +207,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
                         &PRIT, PRIS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -213,6 +221,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
                         &PRST, PRSS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -226,6 +235,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
                         &PRGT, PRGS, ZWSED)
 IF (PRESENT(PFPR)) THEN
   DO JK = KKTB , KKTE
@@ -240,6 +250,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
                           &PRHT, PRHS, ZWSED)
   IF (PRESENT(PFPR)) THEN
     DO JK = KKTB , KKTE
@@ -254,6 +265,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
                                 &PRXT, PRXS, PWSED, PSEA, PTOWN)
     !
     !*      0. DECLARATIONS
@@ -272,6 +284,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)              :: 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
@@ -386,14 +399,37 @@ CONTAINS
                               &      ALOG(PRHODREF(JI,JJ,JK)*ZQP(JI,JJ)) )**XEXCSEDI
           ENDIF
         ENDDO
+
+      ELSEIF(KSPE==5) THEN
+       ! ******* for snow
+        DO JL=1, JCOUNT
+          JI=I1(JL)
+          JJ=I2(JL)
+          !calculation of w
+
+        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
+			      & 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
+                              & PLBDAS(JI,JJ,JK)**(XBS+XEXSEDS) 
+
+          ENDIF
+        ENDDO
+
       ELSE
         ! ******* for other species
         IF(KSPE==3) THEN
           ZFSED=XFSEDR
           ZEXSED=XEXSEDR
-        ELSEIF(KSPE==5) THEN
-          ZFSED=XFSEDS
-          ZEXSED=XEXSEDS
+      !  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 15d0cd78e495cb255015fb6ed29fbfdb5361c748..6396ede4915ce9ee5dc5befe44d37364b2680a7a 100644
--- a/src/MNH/ice4_slow.f90
+++ b/src/MNH/ice4_slow.f90
@@ -173,8 +173,12 @@ IF(LDSOFT) THEN
 ELSE
   PRVDEPS(:) = 0.
   WHERE(ZMASK(:)==1.)
-    PRVDEPS(:) = ( PSSI(:)/(PRHODREF(:)*PAI(:)) ) *                               &
-                 ( X0DEPS*PLBDAS(:)**XEX0DEPS + X1DEPS*PCJ(:)*PLBDAS(:)**XEX1DEPS )
+!    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
+         ( X0DEPS*PLBDAS(:)**XEX0DEPS + (X1DEPS*PCJ(:)*(1+(PLBDAS(:)/(2*XFVELOS)**XALPHAS))**(-XNUS+XEX1DEPS) &
+         *(PLBDAS(:))**(XBS+XEX1DEPS))) ! GAMMAGEN LH_EXTENDED
   END WHERE
 ENDIF
 DO JL=1, KSIZE
@@ -197,10 +201,15 @@ IF(LDSOFT) THEN
 ELSE
   PRIAGGS(:) = 0.
   WHERE(ZMASK(:)==1)
-    PRIAGGS(:) = XFIAGGS * EXP( XCOLEXIS*(PT(:)-XTT) ) &
+!    PRIAGGS(:) = XFIAGGS * EXP( XCOLEXIS*(PT(:)-XTT) ) &
+!                         * PRIT(:)                      &
+!                         * PLBDAS(:)**XEXIAGGS          &
+!                         * PRHODREF(:)**(-XCEXVT)
+      PRIAGGS(:) = XFIAGGS * EXP( XCOLEXIS*(PT(:)-XTT) ) &                !Modif Wurtz GAMMAGEN LH_EXTENDED
                          * PRIT(:)                      &
-                         * PLBDAS(:)**XEXIAGGS          &
-                         * PRHODREF(:)**(-XCEXVT)
+                         * PRST(:) * (1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXIAGGS/XALPHAS)          &
+                         * PRHODREF(:)**(-XCEXVT+1.) &
+                         * ((PLBDAS(:))**(XBS+XEXIAGGS)) ! Thompson
   END WHERE
 ENDIF
 DO JL=1, KSIZE
diff --git a/src/MNH/ice4_tendencies.f90 b/src/MNH/ice4_tendencies.f90
index 49cd599235d8a67c664bd2710e7864adc0a82123..5e5f9000671b65625a576aed5658713cf5b6ff17 100644
--- a/src/MNH/ice4_tendencies.f90
+++ b/src/MNH/ice4_tendencies.f90
@@ -12,6 +12,7 @@ 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, &
@@ -135,6 +136,7 @@ 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
@@ -145,6 +147,7 @@ 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, &
@@ -307,6 +310,7 @@ 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
 !
@@ -400,13 +404,14 @@ 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
+!    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, &
+!                         &ZLBDAS, &
+                         &PLBDAS, & ! Wurtz
                          &ZT, ZRCT, ZRST, &
                          &PRSRIMCG_MR, PB_RS, PB_RG)
     DO JL=1, KSIZE
@@ -469,10 +474,10 @@ 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.
+!  WHERE(ZRST(:)>0.)
+!    ZLBDAS(:)  = MIN(XLBDAS_MAX, XLBS*(PRHODREF(:)*MAX(ZRST(:), XRTMIN(5)))**XLBEXS)
+!  END WHERE
   ZLBDAG(:)=0.
   WHERE(ZRGT(:)>0.)
     ZLBDAG(:)  = XLBG*(PRHODREF(:)*MAX(ZRGT(:), XRTMIN(6)))**XLBEXG
@@ -503,7 +508,8 @@ ENDIF
 CALL ICE4_SLOW(KSIZE, ODSOFT, PCOMPUTE, PRHODREF, ZT, &
               &PSSI, PLVFACT, PLSFACT, &
               &ZRVT, ZRCT, ZRIT, ZRST, ZRGT, &
-              &ZLBDAS, ZLBDAG, &
+          !    &ZLBDAS, ZLBDAG, &
+              &PLBDAS, ZLBDAG, & ! Wurtz
               &ZAI, ZCJ, PHLI_HCF, PHLI_HRI, &
               &PRCHONI, PRVDEPS, PRIAGGS, PRIAUTS, PRVDEPG, &
               &PA_TH, PA_RV, PA_RC, PA_RI, PA_RS, PA_RG)
@@ -540,7 +546,8 @@ END IF
 CALL ICE4_FAST_RS(KSIZE, ODSOFT, PCOMPUTE, &
                  &PRHODREF, PLVFACT, PLSFACT, PPRES, &
                  &ZDV, ZKA, ZCJ, &
-                 &ZLBDAR, ZLBDAS, &
+            !     &ZLBDAR, ZLBDAS, &
+                 &ZLBDAR, PLBDAS, & ! Wurtz
                  &ZT, ZRVT, ZRCT, ZRRT, ZRST, &
                  &PRIAGGS, &
                  &PRCRIMSS, PRCRIMSG, PRSRIMCG, &
@@ -563,7 +570,8 @@ ENDDO
 CALL ICE4_FAST_RG(KSIZE, ODSOFT, PCOMPUTE, KRR, &
                  &PRHODREF, PLVFACT, PLSFACT, PPRES, &
                  &ZDV, ZKA, ZCJ, PCIT, &
-                 &ZLBDAR, ZLBDAS, ZLBDAG, &
+               !  &ZLBDAR, ZLBDAS, ZLBDAG, &
+                 &ZLBDAR, PLBDAS, ZLBDAG, & ! Wurtz
                  &ZT, ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT, &
                  &ZRGSI, ZRGSI_MR(:), &
                  &ZWETG, &
@@ -582,7 +590,8 @@ IF (KRR==7) THEN
   CALL ICE4_FAST_RH(KSIZE, ODSOFT, PCOMPUTE, ZWETG, &
                    &PRHODREF, PLVFACT, PLSFACT, PPRES, &
                    &ZDV, ZKA, ZCJ, &
-                   &ZLBDAS, ZLBDAG, ZLBDAR, ZLBDAH, &
+                !   &ZLBDAS, ZLBDAG, ZLBDAR, ZLBDAH, &
+                   &PLBDAS, ZLBDAG, ZLBDAR, ZLBDAH, & ! Wurtz 
                    &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_ice_c1r3.f90 b/src/MNH/ini_ice_c1r3.f90
index b0a35554a88837b36dc5ae44119426abe6a482a6..a500485c8d88c257931f3a60712039ac597295c9 100644
--- a/src/MNH/ini_ice_c1r3.f90
+++ b/src/MNH/ini_ice_c1r3.f90
@@ -109,6 +109,8 @@ USE MODD_REF
 !
 use mode_msg
 !
+USE MODD_RAIN_ICE_PARAM, ONLY : XFVELOS
+!
 USE MODI_GAMMA
 USE MODI_GAMMA_INC
 USE MODI_READ_XKER_RACCS
@@ -725,15 +727,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='("*****************************************")')
@@ -938,7 +940,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='("*****************************************")')
@@ -1004,7 +1006,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..e095021226a58c6f52ac287c0b0b2316e06a6bd6 100644
--- a/src/MNH/ini_lima_cold_mixed.f90
+++ b/src/MNH/ini_lima_cold_mixed.f90
@@ -58,6 +58,8 @@ USE MODD_REF
 !
 use mode_msg
 !
+USE MODD_RAIN_ICE_PARAM, ONLY : XFVELOS
+
 USE MODI_LIMA_FUNCTIONS
 USE MODI_GAMMA
 USE MODI_GAMMA_INC
@@ -171,8 +173,11 @@ XF1IS = 0.28
 !
 XAS = 0.02
 XBS = 1.9
-XCS = 5.
-XDS = 0.27
+!XCS = 5.
+!XDS = 0.27
+XCS = 40 	! Wurtz	Thompson
+XDS = 0.55	! Wurtz	Thompson	
+XFVELOS = 125 ! Wurtz Thompson
 !
 XCCS = 5.0
 XCXS = 1.0
@@ -785,15 +790,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,                              & ! Thompson
                  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.,                             & ! Wurtz Thompson
                  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,                              & ! Wurtz Thompson
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,XAG, XBS, XAS                         )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -999,7 +1004,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,                              & ! Thompson
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1065,7 +1070,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.,                              & ! Thompson
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1182,7 +1187,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,                              & ! Thompson Wurtz
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1248,7 +1253,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.,                             & ! Thompson
                 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..b2ae256b43dab046bac7dea508baaac426b23316 100644
--- a/src/MNH/ini_param_elec.f90
+++ b/src/MNH/ini_param_elec.f90
@@ -845,20 +845,20 @@ 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,                              &
-              XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
-              PFDINFTY, XKER_Q_RACCSS, XAG, XBS, XAS                      )
+!CALL RRCOLSS (KND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
+!              ZESR, XFR, XCS, XDS, 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,                              &
-              XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
-              PFDINFTY, XKER_Q_RACCS                                      )
+!CALL RZCOLX  (KND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
+!              ZESR, XFR, XCS, XDS, XCR, XDR,                              &
+!              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,                              &
-              XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
-              PFDINFTY, XKER_Q_SACCRG, XAG, XBS, XAS                      )
+!CALL RSCOLRG (KND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
+!              ZESR, XFS, XCS, XDS, XCR, XDR,                              &
+!              XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
+!              PFDINFTY, XKER_Q_SACCRG, XAG, XBS, XAS                      )
 !
 !-------------------------------------------------------------------------------
 !
@@ -877,10 +877,10 @@ 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,                              &
-             XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-             PFDINFTY, XKER_Q_SDRYG                                      )
+!CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+!             ZEGS, XFS, XCG, XDG, XCS, XDS,                              &
+!             XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+!             PFDINFTY, XKER_Q_SDRYG                                      )
 !
 !
 !*      11.2    NI process: Heldson et Farley (1987) parameterization
@@ -896,10 +896,10 @@ IF (CNI_CHARGING == 'HELFA') THEN
   XLBQSDRYGB6H =      MOMG(XALPHAG,XNUG,2.) 
 !
   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,                               &
-               XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-               PFDINFTY, XKER_Q_SDRYGB                                     )
+!  CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+!               ZEGS, 0., XCG, XDG, XCS, XDS,                               &
+!               XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+!               PFDINFTY, XKER_Q_SDRYGB                                     )
 ! Delta vqb1_sg
 ENDIF
 !
@@ -918,10 +918,10 @@ IF (CNI_CHARGING == 'GARDI') THEN
   XLBQSDRYGB6G =      MOMG(XALPHAS,XNUS,6.)    
 !
   IF( .NOT.ALLOCATED(XKER_Q_SDRYGB)) ALLOCATE( XKER_Q_SDRYGB(NDRYLBDAG,NDRYLBDAS) )
-  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, 4., XCG, XDG, XCS, XDS, 4.,                           &
-                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-                PFDINFTY, XKER_Q_SDRYGB                                     )
+!  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+!                ZEGS, 4., XCG, XDG, XCS, XDS, 4.,                           &
+!                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+!                PFDINFTY, XKER_Q_SDRYGB                                     )
 END IF
 !
 !
@@ -941,16 +941,16 @@ IF (CNI_CHARGING == 'SAUN1' .OR. CNI_CHARGING == 'SAUN2' .OR. &
   IF( .NOT.ALLOCATED(XKER_Q_SDRYGB2)) ALLOCATE( XKER_Q_SDRYGB2(NDRYLBDAG,NDRYLBDAS) )
 !
 ! Positive charging region
-  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, XSMP, XCG, XDG, XCS, XDS, (1.+XSNP),                  &
-                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-                PFDINFTY, XKER_Q_SDRYGB1                                    )
+!  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+!                ZEGS, XSMP, XCG, XDG, XCS, XDS, (1.+XSNP),                  &
+!                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+!                PFDINFTY, XKER_Q_SDRYGB1                                    )
 !
 ! Negative charging region
-  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, XSMN, XCG, XDG, XCS, XDS, (1.+XSNN),                  &
-                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-                PFDINFTY, XKER_Q_SDRYGB2                                    )
+!  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+!                ZEGS, XSMN, XCG, XDG, XCS, XDS, (1.+XSNN),                  &
+!                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+!                PFDINFTY, XKER_Q_SDRYGB2                                    )
 ENDIF
 !
 !
@@ -979,10 +979,10 @@ IF (CNI_CHARGING == 'TAKAH') THEN
   XFQSDRYGBT11 = 2. * MOMG(XALPHAG,XNUG,1.) * MOMG(XALPHAS,XNUS,3.)
 !
   IF( .NOT.ALLOCATED(XKER_Q_SDRYGB)) ALLOCATE( XKER_Q_SDRYGB(NDRYLBDAG,NDRYLBDAS) )
-  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, 2., XCG, XDG, XCS, XDS, 2.,                           &
-                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-                PFDINFTY, XKER_Q_SDRYGB                                     )
+!  CALL VQZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+!                ZEGS, 2., XCG, XDG, XCS, XDS, 2.,                           &
+!                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+!                PFDINFTY, XKER_Q_SDRYGB                                     )
 END IF
 !
 !
@@ -998,10 +998,10 @@ IF (CNI_CHARGING == 'TAKAH' .OR. CNI_CHARGING == 'SAP98' .OR. &
   XAUX_LIM2 = 2. * MOMG(XALPHAS,XNUS,1.) * MOMG(XALPHAG,XNUG,1.) 
   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,                               &
-               XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-               PFDINFTY, XKER_Q_LIMSG)
+!  CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+!               ZEGS, 0., XCG, XDG, XCS, XDS,                               &
+!               XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+!               PFDINFTY, XKER_Q_LIMSG)
 ENDIF
 !
 !
@@ -1020,10 +1020,10 @@ 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,                                &
-             XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN,   &
-             PFDINFTY, XKER_Q_RDRYG                                        )
+!CALL RZCOLX (KND, XALPHAG, XNUG, XALPHAR, XNUR,                            & 
+!             ZEGR, XFR, XCG, XDG, XCR, XDR,                                &
+!             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..b2d117da37d7732edd64e794b902bfe437d03da0 100644
--- a/src/MNH/ini_rain_ice.f90
+++ b/src/MNH/ini_rain_ice.f90
@@ -297,8 +297,10 @@ XF2I = 0.14
 !
 XAS = 0.02
 XBS = 1.9
-XCS = 5.1
-XDS = 0.27
+!Values of Thompson (WRF)
+XCS = 12.5 	! Wurtz	GAMMAGEN LH_EXTENDED
+XDS = 0.4	! Wurtz	GAMMAGEN LH_EXTENDED
+XFVELOS = .5 ! Wurtz GAMMAGEN LH_EXTENDED
 !
 XCCS = 5.0
 XCXS = 1.0
@@ -377,8 +379,11 @@ XNUR    = 1.0  ! Exponential law
 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
+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
 !
 XALPHAG = 1.0  ! Exponential law
 XNUG    = 1.0  ! Exponential law
@@ -400,8 +405,9 @@ XLBR   = ( XAR*XCCR*MOMG(XALPHAR,XNUR,XBR) )**(-XLBEXR)
 XLBEXI = 1.0/(-XBI)
 XLBI   = ( XAI*MOMG(XALPHAI,XNUI,XBI) )**(-XLBEXI)
 !
-XLBEXS = 1.0/(XCXS-XBS)
-XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS)
+!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
 !
 XLBEXG = 1.0/(XCXG-XBG)
 XLBG   = ( XAG*XCCG*MOMG(XALPHAG,XNUG,XBG) )**(-XLBEXG)
@@ -413,8 +419,10 @@ 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
-XLBDAS_MAX = ( ZCONC_MAX/XCCS )**(1./XCXS)
+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 ?
 !
 IF (HCLOUD == 'ICE4') THEN
   ALLOCATE( XRTMIN(7) )
@@ -481,9 +489,15 @@ 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
+!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
+            *(ZRHO00)**XCEXVT
 !
 XEXSEDG = (XBG+XDG-XCXG)/(XBG-XCXG)
 XFSEDG  = XCG*XAG*XCCG*MOMG(XALPHAG,XNUG,XBG+XDG)*                         &
@@ -556,10 +570,17 @@ 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 = (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
 !
 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)
@@ -599,8 +620,12 @@ 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  = (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
+XEXIAGGS = -XDS - 2.0 ! GAMMGEN LH_EXTENDED
 !
 GFLAG = .TRUE.
 IF (GFLAG) THEN
@@ -651,15 +676,24 @@ 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= 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
 XEXCRIMSG= XEXCRIMSS
 XCRIMSG  = XCRIMSS
-XSRIMCG  = XCCS*XAS*MOMG(XALPHAS,XNUS,XBS)
-XEXSRIMCG= XCXS-XBS
-XSRIMCG2 = XCCS*XAG*MOMG(XALPHAS,XNUS,XBG)
+!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
 XSRIMCG3 = XFRACM90
-XEXSRIMCG2=XCXS-XBG
+!XEXSRIMCG2=XCXS-XBG
+XEXSRIMCG2=XBS-XBG ! Wurtz
 !
 GFLAG = .TRUE.
 IF (GFLAG) THEN
@@ -689,13 +723,15 @@ 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 = ((XPI**2)/24.0)*XCCS*XCCR*XRHOLW*(ZRHO00**XCEXVT)
+XFRACCSS = XLBS*((XPI**2)/24.0)*XCCR*XRHOLW*(ZRHO00**XCEXVT) ! Wurtz
 !
 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 = (XPI/4.0)*XAS*XCCS*XCCR*(ZRHO00**XCEXVT)
+XFSACCRG = XLBS*(XPI/4.0)*XAS*XCCR*(ZRHO00**XCEXVT) ! Wurtz
 !
 XLBSACCR1   =    MOMG(XALPHAR,XNUR,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSACCR2   = 2.*MOMG(XALPHAR,XNUR,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
@@ -742,15 +778,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,                              & ! Wurtz Thompson
                  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.,                             & ! Wurtz Thompson
                  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,                              & ! Wurtz Thompson
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,  XAG, XBS, XAS                       )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -887,7 +923,8 @@ 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 = (XPI/4.0)*XCOLSG*XCCG*XCCS*XAS*(ZRHO00**XCEXVT)
+XFSDRYG = XLBS*(XPI/4.0)*XCOLSG*XCCG*XAS*(ZRHO00**XCEXVT) ! Wurtz
 !
 XLBSDRYG1   =    MOMG(XALPHAG,XNUG,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSDRYG2   = 2.*MOMG(XALPHAG,XNUG,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
@@ -950,7 +987,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,                  & ! Thompson Wurtz
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1016,7 +1053,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.,                              & ! Thompson
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1082,7 +1119,8 @@ XFWETH = (XPI/4.0)*XCCH*XCH*(ZRHO00**XCEXVT)*MOMG(XALPHAH,XNUH,XDH+2.0)
 !
 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 = (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.)
@@ -1154,7 +1192,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,                              & ! Thompson Wurtz
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1220,7 +1258,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.,                             & ! Thompson
                 XWETLBDAH_MAX, XWETLBDAG_MAX, XWETLBDAH_MIN, XWETLBDAG_MIN, &
                 ZFDINFTY, XKER_GWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1286,7 +1324,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.,                              & ! Thompson
                 XWETLBDAH_MAX, XWETLBDAR_MAX, XWETLBDAH_MIN, XWETLBDAR_MIN, &
                 ZFDINFTY, XKER_RWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1354,8 +1392,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)') &
-                                                      XCCS,XCXS
+  !WRITE(UNIT=KLUOUT,FMT='("           concentration:CC=",E13.6," x=",E13.6)') & ! Modif Wurtz unused
+  !                                                    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 940caeaeefc96dcda0800b7678b8be775815c116..b99ef0f4c5d47741e5fdde0b44e769ee85a1b281 100644
--- a/src/MNH/ini_rain_ice_elec.f90
+++ b/src/MNH/ini_rain_ice_elec.f90
@@ -271,8 +271,10 @@ XF2I = 0.14
 !
 XAS = 0.02
 XBS = 1.9
-XCS = 5.1
-XDS = 0.27
+!Values of Thompson (WRF)
+XCS = 12.5 	! Wurtz	GAMMAGEN LH_EXTENDED
+XDS = 0.4	! Wurtz	GAMMAGEN LH_EXTENDED
+XFVELOS = .5 ! Wurtz GAMMAGEN LH_EXTENDED
 !
 XCCS = 5.0
 XCXS = 1.0
@@ -342,8 +344,11 @@ XNUR    = 1.0  ! Exponential law
 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
+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
 !
 XALPHAG = 1.0  ! Exponential law
 XNUG    = 1.0  ! Exponential law
@@ -365,8 +370,9 @@ XLBR   = (XAR * XCCR * MOMG(XALPHAR,XNUR,XBR))**(-XLBEXR)
 XLBEXI = 1.0 / (-XBI)
 XLBI   = (XAI * MOMG(XALPHAI,XNUI,XBI))**(-XLBEXI)
 !
-XLBEXS = 1.0 / (XCXS - XBS)
-XLBS   = (XAS * XCCS * MOMG(XALPHAS,XNUS,XBS))**(-XLBEXS)
+!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
 !
 XLBEXG = 1.0 / (XCXG - XBG)
 XLBG   = (XAG * XCCG * MOMG(XALPHAG,XNUG,XBG))**(-XLBEXG)
@@ -380,8 +386,10 @@ 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
-XLBDAS_MAX = (ZCONC_MAX / XCCS)**(1./XCXS)
+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 ?
 !
 IF (HCLOUD == 'ICE4') THEN
   ALLOCATE( XRTMIN(7) )
@@ -694,15 +702,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,                              & ! Wurtz Thompson
                  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.,                             & ! Wurtz Thompson
                  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,                              & ! Wurtz Thompson
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,  XAG, XBS, XAS                       )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -905,7 +913,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,                  & ! Thompson Wurtz
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -971,7 +979,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.,                              & ! Thompson
                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
                 ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1089,7 +1097,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,                              & ! Thompson Wurtz
                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
                 ZFDINFTY, XKER_SWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
@@ -1155,7 +1163,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.,                             & ! Thompson
                 XWETLBDAH_MAX, XWETLBDAG_MAX, XWETLBDAH_MIN, XWETLBDAG_MIN, &
                 ZFDINFTY, XKER_GWETH                                        )
   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
diff --git a/src/MNH/modd_rain_ice_param.f90 b/src/MNH/modd_rain_ice_param.f90
index 7568e2a68efc91b9081f88ca55d426724038c3fc..4394deea2c14a1e1bd0c619220d58ff5c0e57559 100644
--- a/src/MNH/modd_rain_ice_param.f90
+++ b/src/MNH/modd_rain_ice_param.f90
@@ -39,7 +39,9 @@
 !             ------------
 !
 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 2416326653944702f6e731bfda3d47c2f2f7d9d3..d0445e41e67890d1c0bac3b8001049b0d4de93e8 100644
--- a/src/MNH/rain_ice_red.f90
+++ b/src/MNH/rain_ice_red.f90
@@ -356,6 +356,8 @@ 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
+
 INTEGER :: IIB           !  Define the domain where is
 INTEGER :: IIE           !  the microphysical sources have to be computed
 INTEGER :: IJB           !
@@ -368,7 +370,7 @@ INTEGER :: JI, JJ, JK
 !For packing
 INTEGER :: IMICRO ! Case r_x>0 locations
 INTEGER, DIMENSION(KSIZE) :: I1,I2,I3 ! Used to replace the COUNT
-INTEGER                   :: JL       ! and PACK intrinsics
+INTEGER                   :: JI, JJ, JK, JL       ! and PACK intrinsics
 !
 !Arrays for nucleation call outisde of LDMICRO points
 REAL,    DIMENSION(KIT, KJT, KKT) :: ZW ! work array
@@ -399,6 +401,7 @@ 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
@@ -571,6 +574,24 @@ ELSE
     ENDDO
   ENDDO
 ENDIF
+
+!Compute lambda_snow parameter ! Modif Wurtz
+!ZT en KELVIN
+PLBDAS(:,:,:)=1000. !Wurtz
+DO JK = 1, IKT
+  DO JJ = 1, IJT
+    DO JI = 1, IIT
+	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
+END DO
 !
 !-------------------------------------------------------------------------------
 !
@@ -600,6 +621,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                                  
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -609,6 +631,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                                 
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -623,6 +646,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                                 
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &
@@ -631,6 +655,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                                 
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &
@@ -710,6 +735,7 @@ 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
@@ -859,6 +885,7 @@ 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, &
@@ -1684,6 +1711,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				   
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -1693,6 +1721,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				   
                                   &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,&
                                   &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,&
                                   &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
@@ -1707,6 +1736,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				   
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &
@@ -1715,6 +1745,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				   
                                    &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
                                    &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, &
                                    &PSEA=PSEA, PTOWN=PTOWN, &
diff --git a/src/MNH/rrcolss.f90 b/src/MNH/rrcolss.f90
index 527165111ecf4d225ce5ec0117c09846d2116b9e..141c644275c398f10e69f719336a42f867b5b3fb 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, &
                          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)
 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                    )
 !     ########################################################################
@@ -151,6 +152,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)
 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 +279,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*ZDS)**PALPHAS)-PFALLR*ZDR**PEXFALLR) & ! GAMMAGEN LH_EXTENDED
                                       * 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*ZDS)**PALPHAS)-PFALLR*ZDRMAX**PEXFALLR) & ! GAMMAGEN LH_EXTENDED
                                    * 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..ac07d23b22365a41e8c1ea8f6fa216f65f2ae42f 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 (Thompson)
 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                     )
 !     ########################################################################
@@ -149,6 +150,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 (Thompson)
 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 +273,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) ! GAMMAGEN LH_EXTENDED
           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) ! GAMMAGEN LH_EXTENDED
             ELSE
             ZCOLLDRMIN = 0.0
           END IF 
diff --git a/src/MNH/rzcolx.f90 b/src/MNH/rzcolx.f90
index 28658241cf1021a29de694cd5a99b85e9c3340d9..c47dd0a3e36a9cab3f89f585af68871fba3415a7 100644
--- a/src/MNH/rzcolx.f90
+++ b/src/MNH/rzcolx.f90
@@ -10,7 +10,7 @@
 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 +29,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 (Thompson)
 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 exponent of specy Z (Thompson)
 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 +51,7 @@ 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                                    )
 !     ########################################################################
@@ -152,8 +154,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 (Thompson)
 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 exponent of specy Z (Thompson)
 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
@@ -246,8 +250,9 @@ DO JLBDAX = 1,SIZE(PRZCOLX(:,:),1)
 !*       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)) ! GAMMAGEN LH_EXTENDED
       END DO
 !
 !*       1.8     Compute the normalization factor by integration over the