diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index f11672ace69e60ffc7f4acb39b5655b6081e5a47..51b902548ce2544ef0f8ee84e7409faab45b4fd0 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -277,7 +277,7 @@ USE MODD_ALLSTATION_n
 !
 USE MODD_PARAM_LIMA, ONLY : LCOLD, LNUCL, LSEDI, LHHONI, LSNOW, LHAIL, LMEYERS,       &
                             NMOD_IFN, NMOM_I, XIFN_CONC, LIFN_HOM, CIFN_SPECIES,      &
-                            CINT_MIXING, NMOD_IMM, NIND_SPECIE,                       &
+                            CINT_MIXING, NMOD_IMM, NIND_SPECIE, LMURAKAMI,            &
                             YSNOW_T=>LSNOW_T, CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA, &
                             XFACTNUC_DEP, XFACTNUC_CON,                               &
                             OWARM=>LWARM, LACTI, ORAIN=>LRAIN, OSEDC=>LSEDC,          &
@@ -1030,6 +1030,7 @@ IF (KMI == 1) THEN
   LSNOW  = .TRUE.
   LHAIL  = .FALSE.
   YSNOW_T = .TRUE.
+  LMURAKAMI = .TRUE.
   CPRISTINE_ICE_LIMA = 'PLAT'
   CHEVRIMED_ICE_LIMA = 'GRAU'
   XFACTNUC_DEP = 1.0  
diff --git a/src/MNH/ini_lima_cold_mixed.f90 b/src/MNH/ini_lima_cold_mixed.f90
index 2a59420e0d0029f87a8837da855080c88da52b10..8aaa67fb33df1ede4e7acbb8228c7fade4bb052d 100644
--- a/src/MNH/ini_lima_cold_mixed.f90
+++ b/src/MNH/ini_lima_cold_mixed.f90
@@ -741,9 +741,13 @@ XCRIMSS  = XNS * (XPI/4.0)*XCOLCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0
 
 XEXCRIMSG= XEXCRIMSS
 XCRIMSG  = XCRIMSS
-XSRIMCG  = XAS*MOMG(XALPHAS,XNUS,XBS)
-XEXSRIMCG= XCXS-XBS
+XSRIMCG  = XNS*XAS*MOMG(XALPHAS,XNUS,XBS)
+XEXSRIMCG= -XBS
 !
+!Pour Murakami 1990
+XSRIMCG2 = XNS*XAG*MOMG(XALPHAS,XNUS,XBG)
+XSRIMCG3 = 0.1
+XEXSRIMCG2=XBS-XBG
 GFLAG = .TRUE.
 IF (GFLAG) THEN
   WRITE(UNIT=ILUOUT0,FMT='("      riming of the aggregates")')
@@ -758,11 +762,13 @@ ZRATE = EXP(LOG(XGAMINC_BOUND_MAX/XGAMINC_BOUND_MIN)/FLOAT(NGAMINC-1))
 !
 ALLOCATE( XGAMINC_RIM1(NGAMINC) )
 ALLOCATE( XGAMINC_RIM2(NGAMINC) )
+ALLOCATE( XGAMINC_RIM4(NGAMINC) )
 !
 DO J1=1,NGAMINC
   ZBOUND = XGAMINC_BOUND_MIN*ZRATE**(J1-1)
   XGAMINC_RIM1(J1) = GAMMA_INC(XNUS+(2.0+XDS)/XALPHAS,ZBOUND)
   XGAMINC_RIM2(J1) = GAMMA_INC(XNUS+XBS/XALPHAS      ,ZBOUND)
+  XGAMINC_RIM4(J1) = GAMMA_INC(XNUS+XBG/XALPHAS      ,ZBOUND) ! Pour Murakami 1990
 END DO
 !
 XRIMINTP1 = XALPHAS / LOG(ZRATE)
diff --git a/src/MNH/lima_cold_slow_processes.f90 b/src/MNH/lima_cold_slow_processes.f90
index 0c9fd34fb6a3e3cd94ff36edfd03abcc79eb09a4..a74b94686bd4113686802ef74b2d3efa9d0133c1 100644
--- a/src/MNH/lima_cold_slow_processes.f90
+++ b/src/MNH/lima_cold_slow_processes.f90
@@ -358,7 +358,7 @@ IF( IMICRO >= 1 ) THEN
       WHERE ( ZLBDAS(:)<XLBDASCNVI_MAX .AND. (ZRST(:)>XRTMIN(5)) &
                                        .AND. (ZSSI(:)<0.0)       )
          ZZW(:) = (ZLBDAS(:)*XDSCNVI_LIM)**(XALPHAS)
-         ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XNS*ZRST(:)*ZLBDAS(:)**XBS) * (ZZW(:)**XNUI) * EXP(-ZZW(:))
+         ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XNS*ZRST(:)*ZLBDAS(:)**XBS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
 !
          ZZW(:) = MIN( ( XR0DEPSI+XR1DEPSI*ZCJ(:) )*ZZX(:),ZRSS(:) )
          ZRIS(:) = ZRIS(:) + ZZW(:)
@@ -466,7 +466,7 @@ IF( IMICRO >= 1 ) THEN
       WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>ZRTMIN(4)) &
                                                             .AND. (ZCIS(:)>ZCTMIN(4)) )
          ZZW1(:,3) = (ZLBDAI(:) / ZLBDAS(:))**3
-         ZZW1(:,1) = (ZCIT(:)*(XNS*ZRST(:)*ZLBDAS(:)**XBS)*EXP(XCOLEXIS*(ZZT(:)-XTT) )) &
+         ZZW1(:,1) = (ZCIT(:)*(XNS*ZRST(:)*ZLBDAS(:)**XBS)*EXP(XCOLEXIS*(ZZT(:)-XTT) ))*ZRHODREF(:) &
                                                     / (ZLBDAI(:)**3)
          ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:,3)),ZCIS(:) )
          ZCIS(:) = ZCIS(:) - ZZW1(:,2)
diff --git a/src/MNH/lima_droplets_riming_snow.f90 b/src/MNH/lima_droplets_riming_snow.f90
index b1c4a8007bddd182015e9aac542bec4c3e0d1c73..d5d32ad8f7a5cb3b01f7a46d415a5fcc91f6f9e4 100644
--- a/src/MNH/lima_droplets_riming_snow.f90
+++ b/src/MNH/lima_droplets_riming_snow.f90
@@ -74,8 +74,8 @@ END MODULE MODI_LIMA_DROPLETS_RIMING_SNOW
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT, XNUS, XALPHAS
-USE MODD_PARAM_LIMA_MIXED, ONLY : NGAMINC, XRIMINTP1, XRIMINTP2, XGAMINC_RIM1, XGAMINC_RIM2, &
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT, XNUS, XALPHAS, LMURAKAMI
+USE MODD_PARAM_LIMA_MIXED, ONLY : NGAMINC, XRIMINTP1, XRIMINTP2, XGAMINC_RIM1, XGAMINC_RIM2, XGAMINC_RIM4, &
                                   XCRIMSS, XEXCRIMSS, XSRIMCG, XEXSRIMCG, XSRIMCG2, XSRIMCG3, XEXSRIMCG2, &
                                   XHMLINTP1, XHMLINTP2, XGAMINC_HMC, XHM_FACTS, XHMTMIN, XHMTMAX
 USE MODD_PARAM_LIMA_COLD,  ONLY : XMNU0, XBS, XFVELOS
@@ -110,119 +110,120 @@ REAL, DIMENSION(:),   INTENT(OUT)   :: P_RS_HMS
 !
 !*       0.2   Declarations of local variables :
 !
-LOGICAL, DIMENSION(SIZE(PRCT))  :: GRIM
-!
-REAL,    DIMENSION(SIZE(PRCT))  :: ZZW1, ZZW2, ZZW3, ZZW4
+REAL,    DIMENSION(SIZE(PRCT))  :: ZZW1, ZZW2, ZZW3, ZZW4, ZZW5
 !
 INTEGER, DIMENSION(SIZE(PRCT))  :: IVEC1,IVEC2      ! Vectors of indices
 REAL,    DIMENSION(SIZE(PRCT))  :: ZVEC1,ZVEC2      ! Work vectors
+INTEGER                         :: JI
 !
 !-------------------------------------------------------------------------------
 !
 !
-P_TH_RIM(:) = 0.
-P_RC_RIM(:) = 0.
-P_CC_RIM(:) = 0.
-P_RS_RIM(:) = 0.
-P_RG_RIM(:) = 0.
-!
-P_RI_HMS(:) = 0.
-P_CI_HMS(:) = 0.
-P_RS_HMS(:) = 0.
 !
-ZZW1(:) = 0.
-ZZW2(:) = 0.
-ZZW3(:) = 0.
-ZZW4(:) = 0.
+DO JI = 1, SIZE(PRCT)
 !
 !*       Cloud droplet riming of the aggregates  
 !        --------------------------------------
 !
+   IF ( PRCT(JI)>XRTMIN(2) .AND. PRST(JI)>XRTMIN(5) .AND. PT(JI)<XTT .AND. LDCOMPUTE(JI) ) THEN
 !
-GRIM(:) = .False.
-GRIM(:) = (PRCT(:)>XRTMIN(2)) .AND. (PRST(:)>XRTMIN(5)) .AND. (PT(:)<XTT) .AND. LDCOMPUTE(:)
-!
-WHERE( GRIM )
-!
-   ZVEC1(:) = PLBDS(:)
+      ZVEC1(JI) = PLBDS(JI)
 !
 !        1.     find the next lower indice for the ZLBDAS in the geometrical
 !               set of Lbda_s used to tabulate some moments of the incomplete 
 !               gamma function
 !
-   ZVEC2(:) = MAX( 1.0001, MIN( REAL(NGAMINC)-0.0001,           &
-                       XRIMINTP1 * LOG( ZVEC1(:) ) + XRIMINTP2 ) )
-   IVEC2(:) = INT( ZVEC2(:) )
-   ZVEC2(:) = ZVEC2(:) - REAL( IVEC2(:) )
+      ZVEC2(JI) = MAX( 1.0001, MIN( REAL(NGAMINC)-0.0001,           &
+                       XRIMINTP1 * LOG( ZVEC1(JI) ) + XRIMINTP2 ) )
+      IVEC2(JI) = INT( ZVEC2(JI) )
+      ZVEC2(JI) = ZVEC2(JI) - REAL( IVEC2(JI) )
 !
 !        2.     perform the linear interpolation of the normalized
 !               "2+XDS"-moment of the incomplete gamma function
 !
-   ZVEC1(:) =   XGAMINC_RIM1( IVEC2(:)+1 )* ZVEC2(:)      &
-                    - XGAMINC_RIM1( IVEC2(:)   )*(ZVEC2(:) - 1.0)
-   ZZW1(:) = ZVEC1(:)
+      ZVEC1(JI) =   XGAMINC_RIM1( IVEC2(JI)+1 )* ZVEC2(JI)      &
+                 - XGAMINC_RIM1( IVEC2(JI)   )*(ZVEC2(JI) - 1.0)
+      ZZW1(JI) = ZVEC1(JI)
 !
 !        3.     perform the linear interpolation of the normalized
 !               "XBS"-moment of the incomplete gamma function
 !
-   ZVEC1(:) =  XGAMINC_RIM2( IVEC2(:)+1 )* ZVEC2(:)      &
-                   - XGAMINC_RIM2( IVEC2(:)   )*(ZVEC2(:) - 1.0)
-   ZZW2(:) = ZVEC1(:)
+      ZVEC1(JI) =  XGAMINC_RIM2( IVEC2(JI)+1 )* ZVEC2(JI)      &
+                - XGAMINC_RIM2( IVEC2(JI)   )*(ZVEC2(JI) - 1.0)
+      ZZW2(JI) = ZVEC1(JI)
 !
 !        4.     riming
 !
    ! Cloud droplets collected
-   P_RC_RIM(:) = - XCRIMSS  * PRCT(:) * PRST(:)*(1+(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
-                                      * PRHODREF(:)**(-XCEXVT+1) &
-                                      * (PLBDS(:)) ** (XEXCRIMSS+XBS)
-   P_CC_RIM(:) = P_RC_RIM(:) *(PCCT(:)/PRCT(:)) ! Lambda_c**3
+      P_RC_RIM(JI) = - XCRIMSS  * PRCT(JI) * PRST(JI)*(1+(XFVELOS/PLBDS(JI))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
+                                         * PRHODREF(JI)**(-XCEXVT+1) &
+                                         * (PLBDS(JI)) ** (XEXCRIMSS+XBS)
+      P_CC_RIM(JI) = P_RC_RIM(JI) *(PCCT(JI)/PRCT(JI)) ! Lambda_c**3
    !
    ! Cloud droplets collected on small aggregates add to snow
-   P_RS_RIM(:) = - P_RC_RIM(:) * ZZW1(:)
+      P_RS_RIM(JI) = - P_RC_RIM(JI) * ZZW1(JI)
    !
    ! Cloud droplets collected on large aggregates add to graupel
-   P_RG_RIM(:) = - P_RC_RIM(:) - P_RS_RIM(:) 
+      P_RG_RIM(JI) = - P_RC_RIM(JI) - P_RS_RIM(JI) 
    !
-   ! Large aggregates collecting droplets add to graupel (instant process ???)
-   ZZW3(:) = PRST(:)*(1.0 - ZZW2(:))/PTSTEP
-   P_RS_RIM(:) = P_RS_RIM(:) - ZZW3(:) 
-   P_RG_RIM(:) = P_RG_RIM(:) + ZZW3(:) 
+      IF (LMURAKAMI) THEN
+   ! Graupel formation based on Murakami
+         ZVEC1(JI) =  XGAMINC_RIM4( IVEC2(JI)+1 )* ZVEC2(JI)      &
+                   - XGAMINC_RIM4( IVEC2(JI)   )*(ZVEC2(JI) - 1.0)
+         ZZW5(JI) = ZVEC1(JI)
+         ZZW3(JI) = XSRIMCG *PRHODREF(JI)*PRST(JI)* PLBDS(JI)**(XEXSRIMCG+XBS) * (1.0 - ZZW2(JI))!/(PTSTEP*PRHODREF(JI))
+         ZZW3(JI) = (P_RG_RIM(JI))*ZZW3(JI)/ &
+                       MAX(1.E-10, & !-20
+                           XSRIMCG3*XSRIMCG2*PRST(JI)*PRHODREF(JI)*PLBDS(JI)**(XEXSRIMCG2)*(1.-ZZW5(JI))- &
+                           XSRIMCG3*ZZW3(JI))
+      ELSE
+      ! Large aggregates collecting droplets add to graupel (instant process ???)
+         ZZW3(JI) = PRST(JI)*(1.0 - ZZW2(JI))/PTSTEP
+      END IF
+      P_RS_RIM(JI) = P_RS_RIM(JI) - ZZW3(JI) 
+      P_RG_RIM(JI) = P_RG_RIM(JI) + ZZW3(JI) 
    !
-   P_TH_RIM(:) = - P_RC_RIM(:)*(PLSFACT(:)-PLVFACT(:))
-END WHERE
-!
+      P_TH_RIM(JI) = - P_RC_RIM(JI)*(PLSFACT(JI)-PLVFACT(JI))
+   ELSE
+      P_TH_RIM(:) = 0.
+      P_RC_RIM(:) = 0.
+      P_CC_RIM(:) = 0.
+      P_RS_RIM(:) = 0.
+      P_RG_RIM(:) = 0.
+   END IF
 !
 !*       Hallett-Mossop ice production (HMS)  
 !        -----------------------------------
 !
-!
-GRIM(:) = .False.
-GRIM(:) = (PT(:)<XHMTMAX)     .AND. (PT(:)>XHMTMIN)     .AND. &
-              (PRST(:)>XRTMIN(5)) .AND. (PRCT(:)>XRTMIN(2)) .AND. &
-               LDCOMPUTE(:)
-!
-WHERE ( GRIM )
-!
-   ZVEC1(:) = PLBDC(:)
-   ZVEC2(:) = MAX( 1.0001, MIN( REAL(NGAMINC)-0.0001,           &
-                         XHMLINTP1 * LOG( ZVEC1(:) ) + XHMLINTP2 ) )
-   IVEC2(:) = INT( ZVEC2(:) )
-   ZVEC2(:) = ZVEC2(:) - REAL( IVEC2(:) )
-   ZVEC1(:) =   XGAMINC_HMC( IVEC2(:)+1 )* ZVEC2(:)      &
-                    - XGAMINC_HMC( IVEC2(:)   )*(ZVEC2(:) - 1.0)
-   ZZW4(:) = ZVEC1(:) ! Large droplets
-!
-   WHERE ( ZZW4(:)<0.99 )
-      P_CI_HMS(:) = - P_RC_RIM(:) * (PCCT(:)/PRCT(:)) * (1.0-ZZW4(:)) * XHM_FACTS * &
-                          MAX( 0.0, MIN( (PT(:)-XHMTMIN)/3.0,(XHMTMAX-PT(:))/2.0 ) ) ! CCHMSI
-!
-      P_RI_HMS(:) = P_CI_HMS(:) * XMNU0                                     ! RCHMSI
-      P_RS_HMS(:) = - P_RI_HMS(:) 
-   END WHERE
-
-END WHERE
-!
-!
+   IF ( PT(JI)<XHMTMAX .AND. PT(JI)>XHMTMIN .AND.                      &
+        PRST(JI)>XRTMIN(5) .AND. PRCT(JI)>XRTMIN(2) .AND. LDCOMPUTE(JI) ) THEN
+!
+      ZVEC1(JI) = PLBDC(JI)
+      ZVEC2(JI) = MAX( 1.0001, MIN( REAL(NGAMINC)-0.0001,           &
+                         XHMLINTP1 * LOG( ZVEC1(JI) ) + XHMLINTP2 ) )
+      IVEC2(JI) = INT( ZVEC2(JI) )
+      ZVEC2(JI) = ZVEC2(JI) - REAL( IVEC2(JI) )
+      ZVEC1(JI) =   XGAMINC_HMC( IVEC2(JI)+1 )* ZVEC2(JI)      &
+                 - XGAMINC_HMC( IVEC2(JI)   )*(ZVEC2(JI) - 1.0)
+      ZZW4(JI) = ZVEC1(JI) ! Large droplets
+!
+      IF ( ZZW4(JI)<0.99 ) THEN
+         P_CI_HMS(JI) = - P_RC_RIM(JI) * (PCCT(JI)/PRCT(JI)) * (1.0-ZZW4(JI)) * XHM_FACTS * &
+                          MAX( 0.0, MIN( (PT(JI)-XHMTMIN)/3.0,(XHMTMAX-PT(JI))/2.0 ) ) ! CCHMSI
+!
+         P_RI_HMS(JI) = P_CI_HMS(JI) * XMNU0                                     ! RCHMSI
+         P_RS_HMS(JI) = - P_RI_HMS(JI)
+      ELSE
+         P_RI_HMS(:) = 0.
+         P_CI_HMS(:) = 0.
+         P_RS_HMS(:) = 0.
+      END IF
+   ELSE
+      P_RI_HMS(:) = 0.
+      P_CI_HMS(:) = 0.
+      P_RS_HMS(:) = 0.
+   END IF
+END DO
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/lima_graupel.f90 b/src/MNH/lima_graupel.f90
index 578f22bbf93fb04be12fb7e1c7edd9f1957cbe53..1883891a1a537d02e912df06c08e33cb07c15406 100644
--- a/src/MNH/lima_graupel.f90
+++ b/src/MNH/lima_graupel.f90
@@ -414,6 +414,8 @@ WHERE( PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:) )
                    ( ZZW5(:)+ZZW6(:) ) *                            &
                    ( PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PT(:)))   ) ) / &
                                    ( PRHODREF(:)*(XLMTT-XCL*(XTT-PT(:))) )   )
+  !We must agregate, at least, the cold species
+   ZRWETG(:)=MAX(ZRWETG(:), ZZW5(:)+ZZW6(:))
 END WHERE
 !
 !            1.f Wet mode and partial conversion to hail
@@ -422,8 +424,8 @@ END WHERE
 ZZW(:) = 0.0
 NHAIL = 0.
 IF (LHAIL) NHAIL = 1. 
-WHERE( LDCOMPUTE(:)      .AND. PRGT(:)>XRTMIN(6)        .AND. PT(:)<XTT       &
-                             .AND. ZRDRYG(:)>=ZRWETG(:) .AND. ZRWETG(:)>0.0 ) 
+WHERE( LDCOMPUTE(:) .AND. PRGT(:)>XRTMIN(6) .AND. PT(:)<XTT .AND. &
+       (ZRDRYG(:)-ZZW2(:)-ZZW3(:))>=(ZRWETG(:)-ZZW5(:)-ZZW6(:)) .AND. ZRWETG(:)>0.0 ) 
 !
 ! Mass of rain and cloud droplets frozen by graupel in wet mode : RCWETG + RRWETG = RWETG - RIWETG - RSWETG
    ZZW7(:) = ZRWETG(:) - ZZW5(:) - ZZW6(:)
@@ -449,8 +451,8 @@ END WHERE
 !            1.g Dry mode
 !            ------------
 !
-WHERE( LDCOMPUTE(:)      .AND. PRGT(:)>XRTMIN(6)       .AND. PT(:)<XTT       &
-                             .AND. ZRDRYG(:)<ZRWETG(:) .AND. ZRDRYG(:)>0.0 )
+WHERE( LDCOMPUTE(:) .AND. PRGT(:)>XRTMIN(6) .AND. PT(:)<XTT .AND.                  &
+       (ZRDRYG(:)-ZZW2(:)-ZZW3(:))<(ZRWETG(:)-ZZW5(:)-ZZW6(:)) .AND. ZRDRYG(:)>0.0 )
    !
    P_RC_DRYG(:) = - ZZW1(:)
    P_CC_DRYG(:) = P_RC_DRYG(:) * PCCT(:)/MAX(PRCT(:),XRTMIN(2))
@@ -470,8 +472,8 @@ END WHERE
 !
 ! BVIE test ZRDRYG<ZZW ?????????????????????????
 !GDRY(:) = (PT(:)<XHMTMAX) .AND. (PT(:)>XHMTMIN)    .AND. (ZRDRYG(:)<ZZW(:))&
-GDRY(:) = (PT(:)<XHMTMAX) .AND. (PT(:)>XHMTMIN)     .AND. (ZRDRYG(:)<ZRWETG(:))&
-               .AND. (PRGT(:)>XRTMIN(6)) .AND. (PRCT(:)>XRTMIN(2)) .AND. LDCOMPUTE(:)
+GDRY(:) = PT(:)<XHMTMAX .AND. PT(:)>XHMTMIN .AND. PRGT(:)>XRTMIN(6) .AND. PRCT(:)>XRTMIN(2) .AND. LDCOMPUTE(:) .AND. &
+          (ZRDRYG(:)-ZZW2(:)-ZZW3(:))<(ZRWETG(:)-ZZW5(:)-ZZW6(:))
 
 ZZX(:)=9999.
 ZVEC1(:)=0.
diff --git a/src/MNH/lima_ice_aggregation_snow.f90 b/src/MNH/lima_ice_aggregation_snow.f90
index 74f1594b561566b8b30373c78a3b82206c50322b..af815a498702f7fab7fe73ea159f7bb2c272b2f6 100644
--- a/src/MNH/lima_ice_aggregation_snow.f90
+++ b/src/MNH/lima_ice_aggregation_snow.f90
@@ -114,7 +114,7 @@ IF (NMOM_I.EQ.1) THEN
 ELSE
    WHERE ( (PRIT(:)>XRTMIN(4)) .AND. (PRST(:)>XRTMIN(5)) .AND. LDCOMPUTE(:) )
       ZZW1(:) = (PLBDI(:) / PLBDS(:))**3
-      ZZW2(:) = (PCIT(:)*(XNS*PRST(:)*PLBDS(:)**XBS)*EXP(XCOLEXIS*(PT(:)-XTT) )) &
+      ZZW2(:) = (PCIT(:)*(XNS*PRST(:)*PLBDS(:)**XBS)*EXP(XCOLEXIS*(PT(:)-XTT) ))*PRHODREF(:) &
            / (PLBDI(:)**3)
       ZZW3(:) = ZZW2(:)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:))
 !
diff --git a/src/MNH/lima_mixed_fast_processes.f90 b/src/MNH/lima_mixed_fast_processes.f90
index 4ce1a80932b529416248070daee6e26af33c1e86..94b622ef3ed710bb825964bb9aeadab84b7025a8 100644
--- a/src/MNH/lima_mixed_fast_processes.f90
+++ b/src/MNH/lima_mixed_fast_processes.f90
@@ -1281,6 +1281,8 @@ WHERE( PRGT1D(:)>XRTMIN(6) )
                   ( ZZW1(:,5)+ZZW1(:,6) ) *                            &
                   ( PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PZT(:)))   ) ) / &
                                   ( PRHODREF(:)*(XLMTT-XCL*(XTT-PZT(:))) )   )
+  !We must agregate, at least, the cold species
+   ZRWETG(:)=MAX(ZRWETG(:), ZZW1(:,5)+ZZW1(:,6))
 END WHERE
 !
 !
@@ -1292,8 +1294,8 @@ END WHERE
 ZZW(:) = 0.0
 NHAIL = 0.
 IF (LHAIL) NHAIL = 1. 
-WHERE( PRGT1D(:)>XRTMIN(6) .AND. PZT(:)<XTT                               &
-                         .AND. ZRDRYG(:)>=ZRWETG(:) .AND. ZRWETG(:)>0.0 ) 
+WHERE( PRGT1D(:)>XRTMIN(6) .AND. PZT(:)<XTT .AND. &
+       (ZRDRYG(:)-ZZW1(:,2)-ZZW1(:,3))>=(ZRWETG(:)-ZZW1(:,5)-ZZW1(:,6)) .AND. ZRWETG(:)>0.0 ) 
 !   
   ZZW(:) = ZRWETG(:) - ZZW1(:,5) - ZZW1(:,6) ! RCWETG+RRWETG
 !   
@@ -1376,8 +1378,8 @@ if ( nbumod == kmi .and. lbu_enable ) then
   end if
 end if
 !
-WHERE( PRGT1D(:)>XRTMIN(6) .AND. PZT(:)<XTT                              &
-                         .AND. ZRDRYG(:)<ZRWETG(:) .AND. ZRDRYG(:)>0.0 ) ! case
+WHERE( PRGT1D(:)>XRTMIN(6) .AND. PZT(:)<XTT .AND. &
+       (ZRDRYG(:)-ZZW1(:,2)-ZZW1(:,3))<(ZRWETG(:)-ZZW1(:,5)-ZZW1(:,6)) .AND. ZRDRYG(:)>0.0 )
   PRCS1D(:) = PRCS1D(:) - ZZW1(:,1)
   PRIS1D(:) = PRIS1D(:) - ZZW1(:,2)
   PRSS1D(:) = PRSS1D(:) - ZZW1(:,3)
@@ -1429,8 +1431,9 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                           Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
 end if
 
-GDRY(:) = (PZT(:)<XHMTMAX) .AND. (PZT(:)>XHMTMIN)      .AND. (ZRDRYG(:)<ZRWETG(:))&
-                           .AND. (PRGT1D(:)>XRTMIN(6)) .AND. (PRCT1D(:)>XRTMIN(2))
+GDRY(:) = (PZT(:)<XHMTMAX) .AND. (PZT(:)>XHMTMIN)            &
+     .AND. (PRGT1D(:)>XRTMIN(6)) .AND. (PRCT1D(:)>XRTMIN(2)) &
+     .AND. (ZRDRYG(:)-ZZW1(:,2)-ZZW1(:,3))<(ZRWETG(:)-ZZW1(:,5)-ZZW1(:,6))
 IGDRY = COUNT( GDRY(:) )
 IF( IGDRY>0 ) THEN
   ALLOCATE(ZVEC1(IGDRY))
diff --git a/src/MNH/lima_snow_deposition.f90 b/src/MNH/lima_snow_deposition.f90
index 89fa334ae85e15bf4ac6ac6633333e2a5be18838..fd28d8c6d889fa51b7ae00afe61a9e709772ebb9 100644
--- a/src/MNH/lima_snow_deposition.f90
+++ b/src/MNH/lima_snow_deposition.f90
@@ -149,7 +149,7 @@ ELSE
 !
          ZZW(:) = ( XR0DEPSI+XR1DEPSI*PCJ(:) )*ZZX(:)
 !
-         ZZW2(:) = ZZW(:)*( XC0DEPSI+XC1DEPSI*PCJ(:) )/( XR0DEPSI+XR1DEPSI*PCJ(:) )
+         ZZW2(:)= ( XC0DEPSI+XC1DEPSI*PCJ(:) )*ZZX(:)
       END WHERE
 !
       P_RI_CNVI(:) = ZZW(:)
diff --git a/src/MNH/modd_param_lima.f90 b/src/MNH/modd_param_lima.f90
index 3d6b4421e8b2dc1225896579b8b085c7925815bf..cb12356ec2fc80a1ba6268b923f486edc1020b6b 100644
--- a/src/MNH/modd_param_lima.f90
+++ b/src/MNH/modd_param_lima.f90
@@ -85,6 +85,7 @@ REAL, DIMENSION(:),    SAVE, ALLOCATABLE :: XFRAC_REF       ! AP compostion in P
 ! 1.3 Ice characteristics
 !
 LOGICAL, SAVE :: LSNOW_T                     ! TRUE to enable snow param. after Wurtz 2021
+LOGICAL, SAVE :: LMURAKAMI                   ! snow + liq -> graupel after Murakami (as in RAIN_ICE_RED)
 CHARACTER(LEN=4), SAVE :: CPRISTINE_ICE_LIMA ! Pristine type PLAT, COLU or BURO
 CHARACTER(LEN=4), SAVE :: CHEVRIMED_ICE_LIMA ! Heavily rimed type GRAU or HAIL
 REAL,SAVE              :: XALPHAI,XNUI,    & ! Pristine ice   distribution parameters
diff --git a/src/MNH/modd_param_lima_mixed.f90 b/src/MNH/modd_param_lima_mixed.f90
index b92e8603ad712a757f8e6c8a320583ba5bb4c39e..6074d09fbc9df88464f18c6cf44873fa0c239477 100644
--- a/src/MNH/modd_param_lima_mixed.f90
+++ b/src/MNH/modd_param_lima_mixed.f90
@@ -115,6 +115,7 @@ INTEGER,SAVE :: NGAMINC                          ! Number of tab. Lbda_s
 REAL, DIMENSION(:), SAVE, ALLOCATABLE          &
                        :: XGAMINC_RIM1,        & ! Tab. incomplete Gamma funct.
                           XGAMINC_RIM2,        & ! for XDS+2 and for XBS
+                          XGAMINC_RIM4,        & ! Murakami
                           XGAMINC_HMC            ! and for the HM process
 !
 REAL,SAVE :: XFRACCSS,                         & ! Constants for the accretion 
diff --git a/src/MNH/modn_param_lima.f90 b/src/MNH/modn_param_lima.f90
index 84a2d83e23282a1bc7c2dcc84f5f53d1a87a658d..e516ad088c7c2e6af4260abccd5bba28c0abe424 100644
--- a/src/MNH/modn_param_lima.f90
+++ b/src/MNH/modn_param_lima.f90
@@ -23,7 +23,7 @@ NAMELIST/NAM_PARAM_LIMA/LCOLD, LNUCL, LSEDI, LSNOW, LHAIL, LHHONI, LMEYERS,&
                         LSNOW_T, CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,   &
                         XALPHAI, XNUI, XALPHAS, XNUS, XALPHAG, XNUG,       &
                         XFACTNUC_DEP, XFACTNUC_CON, NPHILLIPS,             &
-                        LCIBU, XNDEBRIS_CIBU, LRDSF,                       &
+                        LCIBU, XNDEBRIS_CIBU, LRDSF, LMURAKAMI,            &
                         LWARM, LACTI, LRAIN, LSEDC, LACTIT, LBOUND, LSPRO, &
                         LADJ, LKHKO,                                       &
                         NMOD_CCN, XCCN_CONC,                               &