diff --git a/src/MNH/ini_lima_cold_mixed.f90 b/src/MNH/ini_lima_cold_mixed.f90
index 41e43f0a6043974e577c3fd39a098fd0209a8cd6..1de5894808114b3bc8026f3a2bb4b4ef4a77394b 100644
--- a/src/MNH/ini_lima_cold_mixed.f90
+++ b/src/MNH/ini_lima_cold_mixed.f90
@@ -293,8 +293,13 @@ else
    XNUG    = 1.0  ! Exponential law
 end if
 !
-XALPHAH = 1.0  ! Gamma law
-XNUH    = 8.0  ! Gamma law with little dispersion
+if (NMOM_H.GE.2) then
+   XALPHAH = 1.0  ! Gamma law
+   XNUH    = 5.0  ! Gamma law with little dispersion
+else
+   XALPHAH = 1.0  ! Gamma law
+   XNUH    = 8.0  ! Gamma law with little dispersion
+end if
 !
 !*       2.2    Constants for shape parameter
 !
@@ -1513,24 +1518,31 @@ END IF
 !*       9.2.1  Constant for the cloud droplet and cloud ice collection
 !               by the hailstones
 !
-XFWETH = (XPI/4.0)*XCCH*XCH*(ZRHO00**XCEXVT)*MOMG(XALPHAH,XNUH,XDH+2.0)
+XFWETH = (XPI/4.0)*XCH*(ZRHO00**XCEXVT)*MOMG(XALPHAH,XNUH,XDH+2.0)
 !
 !*       9.2.2  Constants for the aggregate collection by the hailstones
 !
-!XFSWETH = (XPI/4.0)*XCCH*XCCS*XAS*(ZRHO00**XCEXVT)
 XFSWETH = XNS*(XPI/4.0)*XAS*(ZRHO00**XCEXVT)
+XFNSWETH= XNS*(XPI/4.0)*XAS*(ZRHO00**XCEXVT)
 !
 XLBSWETH1   =    MOMG(XALPHAH,XNUH,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSWETH2   = 2.*MOMG(XALPHAH,XNUH,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
 XLBSWETH3   =                          MOMG(XALPHAS,XNUS,XBS+2.)
+XLBNSWETH1  =    MOMG(XALPHAH,XNUH,2.)                       
+XLBNSWETH2  = 2.*MOMG(XALPHAH,XNUH,1.)*MOMG(XALPHAS,XNUS,1.) 
+XLBNSWETH3  =                          MOMG(XALPHAS,XNUS,2.)
 !
 !*       9.2.3  Constants for the graupel collection by the hailstones
 !
 XFGWETH = XNG*(XPI/4.0)*XAG*(ZRHO00**XCEXVT)
+XFNGWETH= XNG*(XPI/4.0)*XAG*(ZRHO00**XCEXVT)
 !
 XLBGWETH1   =    MOMG(XALPHAH,XNUH,2.)*MOMG(XALPHAG,XNUG,XBG)
 XLBGWETH2   = 2.*MOMG(XALPHAH,XNUH,1.)*MOMG(XALPHAG,XNUG,XBG+1.)
 XLBGWETH3   =                          MOMG(XALPHAG,XNUG,XBG+2.)
+XLBNGWETH1  =    MOMG(XALPHAH,XNUH,2.)                      
+XLBNGWETH2  = 2.*MOMG(XALPHAH,XNUH,1.)*MOMG(XALPHAG,XNUG,1.)
+XLBNGWETH3  =                          MOMG(XALPHAG,XNUG,2.)
 !
 ! Notice: One magnitude of lambda discretized over 10 points
 !
diff --git a/src/MNH/lima_cold.f90 b/src/MNH/lima_cold.f90
index 3e0a6b3e49e007c6b36efdad8d80ab6c33ee0fe9..1dc71a103a347f167fc92949c56f3b2e175a0ef8 100644
--- a/src/MNH/lima_cold.f90
+++ b/src/MNH/lima_cold.f90
@@ -122,6 +122,7 @@ use modd_budget,     only: lbu_enable,
                            tbudgets
 USE MODD_NSV
 USE MODD_PARAM_LIMA
+USE MODD_PARAM_LIMA_COLD, only: XLBS, XLBEXS, XACCS1, XSPONBUDS1, XSPONBUDS2, XSPONBUDS3, XSPONCOEFS2
 
 use mode_budget,          only: Budget_store_init, Budget_store_end
 
@@ -190,10 +191,18 @@ REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3))  &
                                        PCCT,    & ! Cloud water C. at t
                                        PCRT,    & ! Rain water C. at t
                                        PCIT,    & ! Ice crystal C. at t
+                                       PCST,    & ! Snow/aggregates C. at t      
+                                       PCGT,    & ! Graupel C. at t               
+                                       PCHT,    & ! Hail C. at t                  
                                        !
                                        PCCS,    & ! Cloud water C. source
                                        PCRS,    & ! Rain water C. source
-                                       PCIS       ! Ice crystal C. source
+                                       PCIS,    & ! Ice crystal C. source
+                                       PCSS,    & ! Snow/aggregates C. source   
+                                       PCGS,    & ! Graupel C. source           
+                                       PCHS,    & ! Hail C. source               
+                                       !
+                                       ZWLBDS     ! Snow/aggregates lambda       
 !
 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: PNFS     ! CCN C. available source
                                                   !used as Free ice nuclei for
@@ -250,17 +259,29 @@ IF ( KRR .GE. 7 ) PRHS(:,:,:) = PRS(:,:,:,7)
 PCCT(:,:,:) = 0.
 PCRT(:,:,:) = 0.
 PCIT(:,:,:) = 0.
+PCST(:,:,:) = 0.  
+PCGT(:,:,:) = 0. 
+PCHT(:,:,:) = 0. 
 PCCS(:,:,:) = 0.
 PCRS(:,:,:) = 0.
 PCIS(:,:,:) = 0.
+PCSS(:,:,:) = 0. 
+PCGS(:,:,:) = 0.  
+PCHS(:,:,:) = 0.
 !
 IF ( LWARM ) PCCT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NC)
 IF ( LWARM .AND. LRAIN ) PCRT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NR)
 IF ( LCOLD ) PCIT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NI)
+IF ( NMOM_S.GE.2 ) PCST(:,:,:) = PSVT(:,:,:,NSV_LIMA_NS)        
+IF ( NMOM_G.GE.2 ) PCGT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NG)        
+IF ( NMOM_H.GE.2 ) PCHT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NH)        
 !
 IF ( LWARM ) PCCS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NC)
 IF ( LWARM .AND. LRAIN ) PCRS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NR)
 IF ( LCOLD ) PCIS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NI)
+IF ( NMOM_S.GE.2 ) PCSS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NS)       
+IF ( NMOM_G.GE.2 ) PCGS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NG)        
+IF ( NMOM_H.GE.2 ) PCHS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NH)       
 !
 IF ( NMOD_CCN .GE. 1 ) THEN
    ALLOCATE( PNFS(SIZE(PRHODJ,1),SIZE(PRHODJ,2),SIZE(PRHODJ,3),NMOD_CCN) )
@@ -315,15 +336,25 @@ if ( lbu_enable ) then
   if ( lbudget_rh .and. lhail ) call Budget_store_init( tbudgets(NBUDGET_RH), 'SEDI', prhs(:, :, :) * prhodj(:, :, :) )
   if ( lbudget_sv .and. osedi ) &
       call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'SEDI', pcis(:, :, :) * prhodj(:, :, :) )
+  if (NMOM_S.GE.2) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'SEDI', pcss(:, :, :) * prhodj(:, :, :) )
+  if (NMOM_G.GE.2) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'SEDI', pcgs(:, :, :) * prhodj(:, :, :) )
+  if (NMOM_H.GE.2) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'SEDI', pchs(:, :, :) * prhodj(:, :, :) )
 end if
 
-CALL LIMA_COLD_SEDIMENTATION (OSEDI, KSPLITG, PTSTEP, KMI,     &
-                              PZZ, PRHODJ, PRHODREF,           &
-                              PRIT, PCIT,                      &
-                              PRIS, PRSS, PRGS, PRHS, PCIS,    &
-                              PINPRS, PINPRG,&
-                              PINPRH                  )
-
+if (NMOM_S.GE.2 .AND. NMOM_G.GE.2 .AND. NMOM_H.GE.2) then
+   CALL LIMA_COLD_SEDIMENTATION (OSEDI, KSPLITG, PTSTEP, KMI,     &
+                                 PZZ, PRHODJ, PRHODREF,           &
+                                 PRIT, PCIT,                      &
+                                 PRIS, PRSS, PRGS, PRHS, PCIS,    &
+                                 PINPRS, PINPRG, PINPRH,          &
+                                 PCSS=PCSS, PCGS=PCGS, PCHS=PCHS    )
+else 
+   CALL LIMA_COLD_SEDIMENTATION (OSEDI, KSPLITG, PTSTEP, KMI,     &
+                                 PZZ, PRHODJ, PRHODREF,           &
+                                 PRIT, PCIT,                      &
+                                 PRIS, PRSS, PRGS, PRHS, PCIS,    &
+                                 PINPRS, PINPRG, PINPRH)
+end if
 if ( lbu_enable ) then
   if ( lbudget_ri .and. osedi ) call Budget_store_end( tbudgets(NBUDGET_RI), 'SEDI', pris(:, :, :) * prhodj(:, :, :) )
   if ( lbudget_rs .and. lsnow ) call Budget_store_end( tbudgets(NBUDGET_RS), 'SEDI', prss(:, :, :) * prhodj(:, :, :) )
@@ -331,6 +362,9 @@ if ( lbu_enable ) then
   if ( lbudget_rh .and. lhail ) call Budget_store_end( tbudgets(NBUDGET_RH), 'SEDI', prhs(:, :, :) * prhodj(:, :, :) )
   if ( lbudget_sv .and. osedi ) &
       call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'SEDI', pcis(:, :, :) * prhodj(:, :, :) )
+  if (NMOM_S.GE.2) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'SEDI', pcss(:, :, :) * prhodj(:, :, :) )
+  if (NMOM_G.GE.2) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'SEDI', pcgs(:, :, :) * prhodj(:, :, :) )
+  if (NMOM_H.GE.2) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'SEDI', pchs(:, :, :) * prhodj(:, :, :) )
 end if
 !-------------------------------------------------------------------------------
 !
@@ -377,12 +411,42 @@ END IF
 !
 IF (LSNOW) THEN
 !
-   CALL LIMA_COLD_SLOW_PROCESSES(PTSTEP, KMI, PZZ, PRHODJ,                 &
+   IF(NMOM_S.GE.2) THEN
+      CALL LIMA_COLD_SLOW_PROCESSES(PTSTEP, KMI, PZZ, PRHODJ,                 &
                                  PRHODREF, PEXNREF, PPABST,                &
                                  PTHT, PRVT, PRCT, PRRT, PRIT, PRST, PRGT, &
                                  PTHS, PRVS, PRIS, PRSS,                   &
-                                 PCIT, PCIS                                )
+                                 PCIT, PCIS, PCST=PCST, PCSS=PCSS          )  
+   ELSE
+        CALL LIMA_COLD_SLOW_PROCESSES(PTSTEP, KMI, PZZ, PRHODJ,                 &
+                                   PRHODREF, PEXNREF, PPABST,                &
+                                   PTHT, PRVT, PRCT, PRRT, PRIT, PRST, PRGT, &
+                                   PTHS, PRVS, PRIS, PRSS,                   &
+                                   PCIT, PCIS                                )
+   END IF 
+END IF
+!
+IF (NMOM_S.GE.2) THEN
 !
+!        5.    SPONTANEOUS BREAK-UP (NUMERICAL FILTER)
+!              --------------------
+!
+  if ( lbu_enable .and. lbudget_sv) then
+      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'BRKU', pcss(:, :, :) * prhodj(:, :, :) )             
+  end if
+
+   ZWLBDS(:,:,:) = 1.E10
+   WHERE ((PRSS(:,:,:)>XRTMIN(5)/PTSTEP) .AND. (PCSS(:,:,:)>XCTMIN(5)/PTSTEP ))
+      ZWLBDS(:,:,:) = (XLBS * PCSS(:,:,:) / PRSS(:,:,:) )**XLBEXS
+   END WHERE
+   WHERE (ZWLBDS(:,:,:)<(XACCS1/XSPONBUDS1))
+      PCSS(:,:,:) = PCSS(:,:,:)*MAX((1.+XSPONCOEFS2*(XACCS1/ZWLBDS(:,:,:)-XSPONBUDS1)**2),&
+                                                     (XACCS1/ZWLBDS(:,:,:)/XSPONBUDS3)**3)
+   END WHERE
+!
+  if ( lbu_enable .and. lbudget_sv) then
+      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'BRKU', pcss(:, :, :) * prhodj(:, :, :) )             
+  end if
 END IF
 !
 !------------------------------------------------------------------------------
@@ -404,6 +468,9 @@ IF ( KRR .GE. 7 ) PRS(:,:,:,7) = PRHS(:,:,:)
 IF ( LWARM ) PSVS(:,:,:,NSV_LIMA_NC) = PCCS(:,:,:)
 IF ( LWARM .AND. LRAIN ) PSVS(:,:,:,NSV_LIMA_NR) = PCRS(:,:,:)
 IF ( LCOLD ) PSVS(:,:,:,NSV_LIMA_NI) = PCIS(:,:,:)
+IF ( NMOM_S.GE.2 ) PSVS(:,:,:,NSV_LIMA_NS) = PCSS(:,:,:) 
+IF ( NMOM_G.GE.2 ) PSVS(:,:,:,NSV_LIMA_NG) = PCGS(:,:,:) 
+IF ( NMOM_H.GE.2 ) PSVS(:,:,:,NSV_LIMA_NH) = PCHS(:,:,:) 
 !
 IF ( NMOD_CCN .GE. 1 ) THEN
    PSVS(:,:,:,NSV_LIMA_CCN_FREE:NSV_LIMA_CCN_FREE+NMOD_CCN-1) = PNFS(:,:,:,:)
diff --git a/src/MNH/lima_cold_sedimentation.f90 b/src/MNH/lima_cold_sedimentation.f90
index 1cc8f19d16ec957c5cb74c5a3d5b0c4d1add7530..18132ec16944fe1b1111ce3c9cd02c2d354cedc6 100644
--- a/src/MNH/lima_cold_sedimentation.f90
+++ b/src/MNH/lima_cold_sedimentation.f90
@@ -420,7 +420,7 @@ END IF
                 ALLOCATE(ZLBDAH(ISEDIM))                  
                  DO JL = 1,ISEDIM
                     ZRHS(JL) = PRHS(I1(JL),I2(JL),I3(JL))
-                    ZCHS(JL) = PRHS(I1(JL),I2(JL),I3(JL))            
+                    ZCHS(JL) = PCHS(I1(JL),I2(JL),I3(JL))            
                  END DO
                  ZLBDAH(:)  = 1.E10         
                  WHERE( ZRHS(:)>XRTMIN(7) .AND. ZCHS(:)>XCTMIN(7) )
@@ -473,11 +473,15 @@ END IF
          PRSS(:,:,:) = PRSS(:,:,:) / PTSTEP
          PRGS(:,:,:) = PRGS(:,:,:) / PTSTEP
          PRHS(:,:,:) = PRHS(:,:,:) / PTSTEP
+         IF(NMOM_S.GE.2) PCSS(:,:,:) = PCSS(:,:,:) / PTSTEP
+         IF(NMOM_G.GE.2) PCGS(:,:,:) = PCGS(:,:,:) / PTSTEP
+         IF(NMOM_H.GE.2) PCHS(:,:,:) = PCHS(:,:,:) / PTSTEP 
       END IF
    END IF
 END DO
 !++cb++
 DEALLOCATE(ZRTMIN)
+DEALLOCATE(ZCTMIN)
 !--cb--
 !
 END SUBROUTINE LIMA_COLD_SEDIMENTATION
diff --git a/src/MNH/lima_cold_slow_processes.f90 b/src/MNH/lima_cold_slow_processes.f90
index 28203281b11d6fcec53547f12f4e95ff55725e8a..0f74d6562a71c2d67b28f644df40dd449eb50bdd 100644
--- a/src/MNH/lima_cold_slow_processes.f90
+++ b/src/MNH/lima_cold_slow_processes.f90
@@ -255,8 +255,8 @@ IF( IMICRO >= 1 ) THEN
    ALLOCATE(ZTHS(IMICRO))
 !
    ALLOCATE(ZCIS(IMICRO))
-   if (NMOM_S.GE.2) ALLOCATE(ZCST(IMICRO))
-   if (NMOM_S.GE.2) ALLOCATE(ZCSS(IMICRO)) 
+   ALLOCATE(ZCST(IMICRO))
+   ALLOCATE(ZCSS(IMICRO)) 
 ! 
    ALLOCATE(ZRHODREF(IMICRO)) 
    ALLOCATE(ZZT(IMICRO)) 
@@ -331,7 +331,7 @@ IF( IMICRO >= 1 ) THEN
       ZLBDAI(:) = ( XLBI*ZCIT(:) / ZRIT(:) )**XLBEXI
    END WHERE
    ZLBDAS(:)  = 1.E10
-   IF (LSNOW_T) THEN
+   IF (LSNOW_T .AND. NMOM_S.EQ.1) THEN
       WHERE(ZZT(:)>263.15 .AND. ZRST(:)>XRTMIN(5)) 
          ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
       END WHERE
@@ -339,6 +339,8 @@ IF( IMICRO >= 1 ) THEN
          ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
       END WHERE
       ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
+      ZCST(:) = (XNS*ZRST(:)*ZLBDAS(:)**XBS)
+      ZCSS(:) = (XNS*ZRSS(:)*ZLBDAS(:)**XBS)
    ELSE IF (NMOM_S.GE.2) THEN
 	WHERE (ZRST(:)>XRTMIN(5) .AND. ZCST(:)>XCTMIN(5))
                 ZLBDAS(:) = ( XLBS*ZCST(:) / ZRST(:) )**XLBEXS 
@@ -347,6 +349,8 @@ IF( IMICRO >= 1 ) THEN
       WHERE (ZRST(:)>XRTMIN(5) )
          ZLBDAS(:) = MAX(MIN(XLBDAS_MAX,XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS),XLBDAS_MIN)
       END WHERE
+      ZCST(:) = XCCS*ZLBDAS(:)**XCXS / ZRHODREF(:)
+      ZCSS(:) = XCCS*ZLBDAS(:)**XCXS / ZRHODREF(:)
    END IF
 !
    ZKA(:) = 2.38E-2 + 0.0071E-2 * ( ZZT(:) - XTT )          ! k_a
@@ -374,40 +378,19 @@ IF( IMICRO >= 1 ) THEN
       end if
 !
       ZZW(:) = 0.0
-      if(NMOM_S.GE.2) THEN
-        WHERE ( ZLBDAS(:)<XLBDASCNVI_MAX .AND. (ZRST(:)>XRTMIN(5)) .AND.(ZCST(:)>XCTMIN(5))   &
+      WHERE ( ZLBDAS(:)<XLBDASCNVI_MAX .AND. (ZRST(:)>XRTMIN(5)) .AND.(ZCST(:)>XCTMIN(5))   &
                                          .AND. (ZSSI(:)<0.0)       )
-           ZZW(:) = (ZLBDAS(:)*XDSCNVI_LIM)**(XALPHAS)
-           ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (ZCST(:)) * (ZZW(:)**XNUI) & 
-                                                        * EXP(-ZZW(:))
+         ZZW(:) = (ZLBDAS(:)*XDSCNVI_LIM)**(XALPHAS)
+         ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (ZCST(:)) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
 !
-           ZZW(:) = MIN( ( XR0DEPSI+XR1DEPSI*ZCJ(:) )*ZZX(:),ZRSS(:) )
-           ZRIS(:) = ZRIS(:) + ZZW(:)
-           ZRSS(:) = ZRSS(:) - ZZW(:)
+         ZZW(:) = MIN( ( XR0DEPSI+XR1DEPSI*ZCJ(:) )*ZZX(:),ZRSS(:) )
+         ZRIS(:) = ZRIS(:) + ZZW(:)
+         ZRSS(:) = ZRSS(:) - ZZW(:)
 !
-           ZZW(:) = MIN(ZZW(:)*( XC0DEPSI+XC1DEPSI*ZCJ(:) )/( XR0DEPSI+XR1DEPSI*ZCJ(:)),ZCSS(:) )
-           ZCIS(:) = ZCIS(:) + ZZW(:)
-           ZCSS(:) = ZCSS(:) - ZZW(:)
-        END WHERE
-      else
-        WHERE ( ZRST(:)>XRTMIN(5) )
-           ZLBDAS(:)  = MIN( XLBDAS_MAX,                                           &
-                             XLBS*( ZRHODREF(:)*MAX( ZRST(:),XRTMIN(5) ) )**XLBEXS )
-        END WHERE
-        WHERE ( ZLBDAS(:)<XLBDASCNVI_MAX .AND. (ZRST(:)>XRTMIN(5)) &
-                                         .AND. (ZSSI(:)<0.0)       )
-           ZZW(:) = (ZLBDAS(:)*XDSCNVI_LIM)**(XALPHAS)
-           ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XCCS*ZLBDAS(:)**XCXS)/ZRHODREF(:) * (ZZW(:)**XNUI) &
-                                                                       * EXP(-ZZW(:))
-!
-           ZZW(:) = MIN( ( XR0DEPSI+XR1DEPSI*ZCJ(:) )*ZZX(:),ZRSS(:) )
-           ZRIS(:) = ZRIS(:) + ZZW(:)
-           ZRSS(:) = ZRSS(:) - ZZW(:)
-!
-           ZZW(:) = ZZW(:)*( XC0DEPSI+XC1DEPSI*ZCJ(:) )/( XR0DEPSI+XR1DEPSI*ZCJ(:) )
-           ZCIS(:) = ZCIS(:) + ZZW(:)
-        END WHERE
-      end if
+         ZZW(:) = MIN(ZZW(:)*( XC0DEPSI+XC1DEPSI*ZCJ(:) )/( XR0DEPSI+XR1DEPSI*ZCJ(:)),ZCSS(:) )
+         ZCIS(:) = ZCIS(:) + ZZW(:)
+         ZCSS(:) = ZCSS(:) - ZZW(:)
+      END WHERE
 !
 ! Budget storage
       if ( nbumod == kmi .and. lbu_enable ) then
@@ -434,28 +417,17 @@ IF( IMICRO >= 1 ) THEN
       end if
 
       ZZW(:) = 0.0
-      if (NMOM_S.GE.2) then
          WHERE ( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>ZRTMIN(5)) )
-            ZZW(:) = ( ZSSI(:)/(ZRHODREF(:)*ZAI(:)) ) * ZCST *                   &
-                     ( X0DEPS*ZLBDAS(:)**XEX0DEPS + X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS )
-            ZZW(:) =    MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
-                      - MIN( ZRSS(:),ABS(ZZW(:)) )*(0.5-SIGN(0.5,ZZW(:)))
-            ZRSS(:) = ZRSS(:) + ZZW(:)
-            ZRVS(:) = ZRVS(:) - ZZW(:)
-            ZTHS(:) = ZTHS(:) + ZZW(:)*ZLSFACT(:)
-         END WHERE
-      else
-         WHERE ( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>ZRTMIN(5)) )
-            ZZW(:) = ( ZSSI(:)/(ZRHODREF(:)*ZAI(:)) ) *                    &
-                     ( X0DEPS*ZLBDAS(:)**XEX0DEPS + X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS )
-
+            ZZW(:) = ( ZSSI(:)/(ZRHODREF(:)*ZAI(:)) ) * ZCST(:) *  &
+                      ( X0DEPS*ZLBDAS(:)**XEX0DEPS +               &
+                        X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS *        &
+                             (1+0.5*(XFVELOS/ZLBDAS(:))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) )
             ZZW(:) =    MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
                       - MIN( ZRSS(:),ABS(ZZW(:)) )*(0.5-SIGN(0.5,ZZW(:)))
             ZRSS(:) = ZRSS(:) + ZZW(:)
             ZRVS(:) = ZRVS(:) - ZZW(:)
             ZTHS(:) = ZTHS(:) + ZZW(:)*ZLSFACT(:)
          END WHERE
-      end if
 ! Budget storage
       if ( nbumod == kmi .and. lbu_enable ) then
         if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'DEPS', &
@@ -481,42 +453,22 @@ IF( IMICRO >= 1 ) THEN
       end if
 
       ZZW(:) = 0.0
-      if (NMOM_S.GE.2) then
-         WHERE ( (ZLBDAI(:)<XLBDAICNVS_LIM) .AND. (ZCIT(:)>XCTMIN(4)) &
+      WHERE ( (ZLBDAI(:)<XLBDAICNVS_LIM) .AND. (ZCIT(:)>XCTMIN(4)) &
                                             .AND. (ZSSI(:)>0.0)       )
-            ZZW(:) = (ZLBDAI(:)*XDICNVS_LIM)**(XALPHAI)
-            ZZX(:) = ( ZSSI(:)/ZAI(:) )*ZCIT(:) * (ZZW(:)**XNUI) *EXP(-ZZW(:))
+         ZZW(:) = (ZLBDAI(:)*XDICNVS_LIM)**(XALPHAI)
+         ZZX(:) = ( ZSSI(:)/ZAI(:) )*ZCIT(:) * (ZZW(:)**XNUI) *EXP(-ZZW(:))
 !
 ! Correction BVIE
 !         ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:)/ZRHODREF(:) &
-            ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:) &
-                               ,ZRIS(:) ) + ZRTMIN(5), ZRTMIN(5) ) - ZRTMIN(5)
-            ZRIS(:) = ZRIS(:) - ZZW(:)
-            ZRSS(:) = ZRSS(:) + ZZW(:)
+         ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:) , ZRIS(:) ) , 0. )
+         ZRIS(:) = ZRIS(:) - ZZW(:)
+         ZRSS(:) = ZRSS(:) + ZZW(:)
 !
-            ZZW(:) = MIN( ZZW(:)*(( XC0DEPIS+XC1DEPIS*ZCJ(:) )                   &
+         ZZW(:) = MIN( ZZW(:)*(( XC0DEPIS+XC1DEPIS*ZCJ(:) )                   &
                                 /( XR0DEPIS+XR1DEPIS*ZCJ(:) )),ZCIS(:) )
-            ZCIS(:) = ZCIS(:) - ZZW(:)
-            ZCSS(:) = ZCSS(:) + ZZW(:)
+         ZCIS(:) = ZCIS(:) - ZZW(:)
+         ZCSS(:) = ZCSS(:) + ZZW(:)
          END WHERE
-      else
-         WHERE ( (ZLBDAI(:)<XLBDAICNVS_LIM) .AND. (ZCIT(:)>XCTMIN(4)) &
-                                            .AND. (ZSSI(:)>0.0)       )
-            ZZW(:) = (ZLBDAI(:)*XDICNVS_LIM)**(XALPHAI)
-            ZZX(:) = ( ZSSI(:)/ZAI(:) )*ZCIT(:) * (ZZW(:)**XNUI) *EXP(-ZZW(:))
-!
-! Correction BVIE
-!         ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:)/ZRHODREF(:) &
-            ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:) &
-                               ,ZRIS(:) ) + ZRTMIN(5), ZRTMIN(5) ) - ZRTMIN(5)
-            ZRIS(:) = ZRIS(:) - ZZW(:)
-            ZRSS(:) = ZRSS(:) + ZZW(:)
-!
-            ZZW(:) = MIN( ZZW(:)*(( XC0DEPIS+XC1DEPIS*ZCJ(:) )                   &
-                                /( XR0DEPIS+XR1DEPIS*ZCJ(:) )),ZCIS(:) )
-            ZCIS(:) = ZCIS(:) - ZZW(:)
-         END WHERE
-      end if
 !
 ! Budget storage
       if ( nbumod == kmi .and. lbu_enable ) then
@@ -543,35 +495,18 @@ IF( IMICRO >= 1 ) THEN
                                                Unpack( zcis(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
       end if
 !
-      if (NMOM_S.GE.2) then
-         WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>ZRTMIN(4)) &
+      WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>ZRTMIN(4)) &
                                                                .AND. (ZCIS(:)>ZCTMIN(4)) )
-              ZZW1(:,3) = (ZLBDAI(:) / ZLBDAS(:))**3
-              ZZW1(:,1) = (ZCIT(:)*ZCST(:)*EXP( XCOLEXIS*(ZZT(:)-XTT) )) & 
-                                                 / (ZLBDAI(:)**3)
-              ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:,3)),ZCIS(:) )
-              ZCIS(:) = ZCIS(:) - ZZW1(:,2)
-!
-              ZZW1(:,1) = ZZW1(:,1) / ZLBDAI(:)**XBI
-              ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_RLARGE1+XAGGS_RLARGE2*ZZW1(:,3)),ZRIS(:) )
-              ZRIS(:) = ZRIS(:) - ZZW1(:,2)
-              ZRSS(:) = ZRSS(:) + ZZW1(:,2)
-         END WHERE
-      else
-         WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>ZRTMIN(4)) &
-                                                         .AND. (ZCIS(:)>ZCTMIN(4)) )
-              ZZW1(:,3) = (ZLBDAI(:) / ZLBDAS(:))**3
-              ZZW1(:,1) = (ZCIT(:)*(XCCS*ZLBDAS(:)**XCXS)/ZRHODREF(:)*EXP( XCOLEXIS*(ZZT(:)-XTT) )) &
-                                                         / (ZLBDAI(:)**3)
-              ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:,3)),ZCIS(:) )
-              ZCIS(:) = ZCIS(:) - ZZW1(:,2)
-!
-              ZZW1(:,1) = ZZW1(:,1) / ZLBDAI(:)**XBI
-              ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_RLARGE1+XAGGS_RLARGE2*ZZW1(:,3)),ZRIS(:) )
-              ZRIS(:) = ZRIS(:) - ZZW1(:,2)
-              ZRSS(:) = ZRSS(:) + ZZW1(:,2)
-         END WHERE
-      end if
+         ZZW1(:,3) = (ZLBDAI(:) / ZLBDAS(:))**3
+         ZZW1(:,1) = (ZCIT(:)*ZCST(:)*EXP( XCOLEXIS*(ZZT(:)-XTT) )) / (ZLBDAI(:)**3)
+         ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:,3)),ZCIS(:) )
+         ZCIS(:) = ZCIS(:) - ZZW1(:,2)
+!
+         ZZW1(:,1) = ZZW1(:,1) / ZLBDAI(:)**XBI
+         ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_RLARGE1+XAGGS_RLARGE2*ZZW1(:,3)),ZRIS(:) )
+         ZRIS(:) = ZRIS(:) - ZZW1(:,2)
+         ZRSS(:) = ZRSS(:) + ZZW1(:,2)
+      END WHERE
 ! Budget storage
       if ( nbumod == kmi .and. lbu_enable ) then
         if ( lbudget_ri ) call Budget_store_end( tbudgets(NBUDGET_RI), 'AGGS', &
@@ -618,8 +553,8 @@ IF( IMICRO >= 1 ) THEN
   DEALLOCATE(ZRSS)
   DEALLOCATE(ZTHS)
   DEALLOCATE(ZCIS)
-  if (NMOM_S.GE.2)  DEALLOCATE(ZCSS)  
-  if (NMOM_S.GE.2)  DEALLOCATE(ZCST)  
+  DEALLOCATE(ZCSS)  
+  DEALLOCATE(ZCST)  
   DEALLOCATE(ZRHODREF) 
   DEALLOCATE(ZZT) 
   DEALLOCATE(ZPRES) 
diff --git a/src/MNH/lima_mixed.f90 b/src/MNH/lima_mixed.f90
index 5c79c29d892ccce05028b90aa3be65fcd02c5444..2045b158c24f4e476234b284ecd06052ca4cb209 100644
--- a/src/MNH/lima_mixed.f90
+++ b/src/MNH/lima_mixed.f90
@@ -97,6 +97,7 @@ END MODULE MODI_LIMA_MIXED
 !  P. Wautelet 28/05/2020: bugfix: correct array start for PSVT and PSVS
 !  P. Wautelet 02/02/2021: budgets: add missing source terms for SV budgets in LIMA
 !  J. Wurtz       03/2022: new snow characteristics
+!  M. Taufour     07/2022: add concentration for snow, graupel, hail
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -109,10 +110,13 @@ USE MODD_CST,              ONLY: XP00, XRD, XRV, XMV, XMD, XCPD, XCPV,       &
 USE MODD_NSV
 USE MODD_PARAMETERS,       ONLY: JPHEXT, JPVEXT
 USE MODD_PARAM_LIMA,       ONLY: NMOD_IFN, XRTMIN, XCTMIN, LWARM, LCOLD,     &
-                                 NMOD_CCN, NMOD_IMM, LRAIN, LSNOW, LHAIL, LSNOW_T
+                                 NMOD_CCN, NMOD_IMM, LRAIN, LSNOW, LHAIL, LSNOW_T,    &
+                                 NMOM_S, NMOM_G, NMOM_H
 USE MODD_PARAM_LIMA_WARM,  ONLY: XLBC, XLBEXC, XLBR, XLBEXR
-USE MODD_PARAM_LIMA_COLD,  ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XSCFAC, XLBDAS_MAX, XLBDAS_MIN, XTRANS_MP_GAMMAS
-USE MODD_PARAM_LIMA_MIXED, ONLY: XLBG, XLBEXG, XLBH, XLBEXH
+USE MODD_PARAM_LIMA_COLD,  ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XSCFAC, &
+                                 XLBDAS_MAX, XLBDAS_MIN, XTRANS_MP_GAMMAS, &
+                                 XBS, XCCS, XCXS, XNS
+USE MODD_PARAM_LIMA_MIXED, ONLY: XLBG, XLBEXG, XCCG, XCXG, XLBH, XLBEXH, XCCH, XCXH
 
 use mode_tools,            only: Countjv
 
@@ -175,17 +179,16 @@ REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3))  &
                                        PCCT,    & ! Cloud water C. at t
                                        PCRT,    & ! Rain water C. at t
                                        PCIT,    & ! Ice crystal C. at t
-                                       PCST,    & ! Snow/aggregates C. at t       ! MTaufour
-                                       PCGT,    & ! Graupel C. at t               ! MTaufour
-                                       PCHT,    & ! Hail C. at t                  ! MTaufour
+                                       PCST,    & ! Snow/aggregates C. at t   
+                                       PCGT,    & ! Graupel C. at t           
+                                       PCHT,    & ! Hail C. at t 
                                        !
                                        PCCS,    & ! Cloud water C. source
                                        PCRS,    & ! Rain water C. source
                                        PCIS,    & ! Ice crystal C. source
-                                       PCSS,    & ! Snow/aggregates C. source    ! MTaufour
-                                       PCGS,    & ! Graupel C. source            ! MTaufour
-                                       PCHS       ! Hail C. source               ! MTaufour
-!
+                                       PCSS,    & ! Snow/aggregates C. source
+                                       PCGS,    & ! Graupel C. source
+                                       PCHS       ! Hail C. source
 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: PNFS     ! CCN C. available source
                                                   !used as Free ice nuclei for
                                                   !HOMOGENEOUS nucleation of haze
@@ -217,9 +220,9 @@ REAL, DIMENSION(:), ALLOCATABLE   :: ZRHT    ! Hail m.r. at t
 REAL, DIMENSION(:), ALLOCATABLE   :: ZCCT    ! Cloud water conc. at t
 REAL, DIMENSION(:), ALLOCATABLE   :: ZCRT    ! Rain water conc. at t
 REAL, DIMENSION(:), ALLOCATABLE   :: ZCIT    ! Pristine ice conc. at t
-REAL, DIMENSION(:), ALLOCATABLE   :: ZCST    ! Snow/aggregates conc. at t  ! MTaufour
-REAL, DIMENSION(:), ALLOCATABLE   :: ZCGT    ! Graupel conc. at t          ! MTaufour
-REAL, DIMENSION(:), ALLOCATABLE   :: ZCHT    ! Hail conc. at t             ! MTaufour
+REAL, DIMENSION(:), ALLOCATABLE   :: ZCST    ! Snow/aggregates conc. at t 
+REAL, DIMENSION(:), ALLOCATABLE   :: ZCGT    ! Graupel conc. at t
+REAL, DIMENSION(:), ALLOCATABLE   :: ZCHT    ! Hail conc. at t
 !
 REAL, DIMENSION(:), ALLOCATABLE   :: ZRVS    ! Water vapor m.r. source
 REAL, DIMENSION(:), ALLOCATABLE   :: ZRCS    ! Cloud water m.r. source
@@ -234,9 +237,9 @@ REAL, DIMENSION(:), ALLOCATABLE   :: ZTHS    ! Theta source
 REAL, DIMENSION(:),   ALLOCATABLE :: ZCCS    ! Cloud water conc. source
 REAL, DIMENSION(:),   ALLOCATABLE :: ZCRS    ! Rain water conc. source
 REAL, DIMENSION(:),   ALLOCATABLE :: ZCIS    ! Pristine ice conc. source
-REAL, DIMENSION(:),   ALLOCATABLE :: ZCSS    ! Snow/aggregates conc. at t  ! MTaufour
-REAL, DIMENSION(:),   ALLOCATABLE :: ZCGS    ! Graupel conc. at t          ! MTaufour
-REAL, DIMENSION(:),   ALLOCATABLE :: ZCHS    ! Hail conc. at t             ! MTaufour
+REAL, DIMENSION(:),   ALLOCATABLE :: ZCSS    ! Snow/aggregates conc. at t
+REAL, DIMENSION(:),   ALLOCATABLE :: ZCGS    ! Graupel conc. at t
+REAL, DIMENSION(:),   ALLOCATABLE :: ZCHS    ! Hail conc. at t
 !
 REAL, DIMENSION(:,:), ALLOCATABLE :: ZIFS    ! Free Ice nuclei conc. source
 REAL, DIMENSION(:,:), ALLOCATABLE :: ZINS    ! Nucleated Ice nuclei conc. source
@@ -311,29 +314,29 @@ IF ( KRR .GE. 7 ) PRHS(:,:,:) = PRS(:,:,:,7)
 PCCT(:,:,:) = 0.
 PCRT(:,:,:) = 0.
 PCIT(:,:,:) = 0.
-PCST(:,:,:) = 0.  ! MTaufour
-PCGT(:,:,:) = 0.  ! MTaufour
-PCHT(:,:,:) = 0.  ! MTaufour
+PCST(:,:,:) = 0. 
+PCGT(:,:,:) = 0. 
+PCHT(:,:,:) = 0.
 PCCS(:,:,:) = 0.
 PCRS(:,:,:) = 0.
 PCIS(:,:,:) = 0.
-PCSS(:,:,:) = 0.  ! MTaufour
-PCGS(:,:,:) = 0.  ! MTaufour
-PCHS(:,:,:) = 0.  ! MTaufour
+PCSS(:,:,:) = 0. 
+PCGS(:,:,:) = 0. 
+PCHS(:,:,:) = 0. 
 !
 IF ( LWARM ) PCCT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NC) 
 IF ( LWARM .AND. LRAIN ) PCRT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NR)
 IF ( LCOLD ) PCIT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NI)
-IF ( LCOLD .AND. LSNOW ) PCST(:,:,:) = PSVT(:,:,:,NSV_LIMA_NS)          ! MTaufour
-IF ( LCOLD .AND. LSNOW ) PCGT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NG)          ! MTaufour
-IF ( LCOLD .AND. LHAIL ) PCHT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NH)          ! MTaufour
+IF ( NMOM_S.GE.2 ) PCST(:,:,:) = PSVT(:,:,:,NSV_LIMA_NS)
+IF ( NMOM_G.GE.2 ) PCGT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NG)
+IF ( NMOM_H.GE.2 ) PCHT(:,:,:) = PSVT(:,:,:,NSV_LIMA_NH)
 !
 IF ( LWARM ) PCCS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NC)
 IF ( LWARM .AND. LRAIN ) PCRS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NR)
 IF ( LCOLD ) PCIS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NI)
-IF ( LCOLD .AND. LSNOW ) PCSS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NS)          ! MTaufour
-IF ( LCOLD .AND. LSNOW ) PCGS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NG)          ! MTaufour
-IF ( LCOLD .AND. LHAIL ) PCHS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NH)          ! MTaufour
+IF ( NMOM_S.GE.2 ) PCSS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NS)
+IF ( NMOM_G.GE.2 ) PCGS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NG)
+IF ( NMOM_H.GE.2 ) PCHS(:,:,:) = PSVS(:,:,:,NSV_LIMA_NH)
 !
 IF ( NMOD_CCN .GE. 1 ) THEN
    ALLOCATE( PNFS(SIZE(PRHODJ,1),SIZE(PRHODJ,2),SIZE(PRHODJ,3),NMOD_CCN) )
@@ -416,10 +419,10 @@ IF( IMICRO >= 0 ) THEN
    !
    ALLOCATE(ZCCT(IMICRO)) 
    ALLOCATE(ZCRT(IMICRO)) 
-   ALLOCATE(ZCIT(IMICRO)) 
-   ALLOCATE(ZCST(IMICRO)) ! MTaufour 
-   ALLOCATE(ZCGT(IMICRO)) ! MTaufour
-   ALLOCATE(ZCHT(IMICRO)) ! MTaufour   
+   ALLOCATE(ZCIT(IMICRO))
+   ALLOCATE(ZCST(IMICRO)) 
+   ALLOCATE(ZCGT(IMICRO))
+   ALLOCATE(ZCHT(IMICRO))
    !
    ALLOCATE(ZRVS(IMICRO))  
    ALLOCATE(ZRCS(IMICRO)) 
@@ -432,10 +435,10 @@ IF( IMICRO >= 0 ) THEN
    !
    ALLOCATE(ZCCS(IMICRO)) 
    ALLOCATE(ZCRS(IMICRO)) 
-   ALLOCATE(ZCIS(IMICRO)) 
-   ALLOCATE(ZCSS(IMICRO)) ! MTaufour 
-   ALLOCATE(ZCGS(IMICRO)) ! MTaufour
-   ALLOCATE(ZCHS(IMICRO)) ! MTaufour   
+   ALLOCATE(ZCIS(IMICRO))
+   ALLOCATE(ZCSS(IMICRO)) 
+   ALLOCATE(ZCGS(IMICRO)) 
+   ALLOCATE(ZCHS(IMICRO))
    ALLOCATE(ZIFS(IMICRO,NMOD_IFN))
    ALLOCATE(ZINS(IMICRO,NMOD_IFN))
    !
@@ -455,9 +458,9 @@ IF( IMICRO >= 0 ) THEN
       ZCCT(JL) = PCCT(I1(JL),I2(JL),I3(JL))
       ZCRT(JL) = PCRT(I1(JL),I2(JL),I3(JL))
       ZCIT(JL) = PCIT(I1(JL),I2(JL),I3(JL))
-      ZCST(JL) = PCST(I1(JL),I2(JL),I3(JL)) ! MTaufour
-      ZCGT(JL) = PCGT(I1(JL),I2(JL),I3(JL)) ! MTaufour
-      ZCHT(JL) = PCHT(I1(JL),I2(JL),I3(JL)) ! MTaufour      
+      ZCST(JL) = PCST(I1(JL),I2(JL),I3(JL))
+      ZCGT(JL) = PCGT(I1(JL),I2(JL),I3(JL)) 
+      ZCHT(JL) = PCHT(I1(JL),I2(JL),I3(JL))
       !
       ZRVS(JL) = PRVS(I1(JL),I2(JL),I3(JL))
       ZRCS(JL) = PRCS(I1(JL),I2(JL),I3(JL))
@@ -471,9 +474,9 @@ IF( IMICRO >= 0 ) THEN
       ZCCS(JL) = PCCS(I1(JL),I2(JL),I3(JL))
       ZCRS(JL) = PCRS(I1(JL),I2(JL),I3(JL))
       ZCIS(JL) = PCIS(I1(JL),I2(JL),I3(JL))
-      ZCSS(JL) = PCSS(I1(JL),I2(JL),I3(JL)) ! MTaufour
-      ZCGS(JL) = PCGS(I1(JL),I2(JL),I3(JL)) ! MTaufour
-      ZCHS(JL) = PCHS(I1(JL),I2(JL),I3(JL)) ! MTaufour      
+      ZCSS(JL) = PCSS(I1(JL),I2(JL),I3(JL)) 
+      ZCGS(JL) = PCGS(I1(JL),I2(JL),I3(JL)) 
+      ZCHS(JL) = PCHS(I1(JL),I2(JL),I3(JL))       
       DO JMOD_IFN = 1, NMOD_IFN
          ZIFS(JL,JMOD_IFN) = PIFS(I1(JL),I2(JL),I3(JL),JMOD_IFN)
          ZINS(JL,JMOD_IFN) = PINS(I1(JL),I2(JL),I3(JL),JMOD_IFN)
@@ -541,27 +544,51 @@ IF( IMICRO >= 0 ) THEN
       ZLBDAI(:) = ( XLBI*ZCIT(:) / ZRIT(:) )**XLBEXI
    END WHERE
    ZLBDAS(:)  = 1.E10
-   IF (LSNOW_T) THEN
+   IF (LSNOW_T .AND. NMOM_S.EQ.1) THEN
       WHERE(ZZT(:)>263.15 .AND. ZRST(:)>XRTMIN(5)) 
          ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
       END WHERE
       WHERE(ZZT(:)<=263.15 .AND. ZRST(:)>XRTMIN(5)) 
          ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
       END WHERE
-      ZLBDAS(:) = ZLBDAS(:)*XTRANS_MP_GAMMAS
+      ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
+      ZCST(:) = (XNS*ZRST(:)*ZLBDAS(:)**XBS)
+      ZCSS(:) = (XNS*ZRSS(:)*ZLBDAS(:)**XBS)
+   ELSE IF (NMOM_S.GE.2) THEN
+	WHERE (ZRST(:)>XRTMIN(5) .AND. ZCST(:)>XCTMIN(5))
+                ZLBDAS(:) = ( XLBS*ZCST(:) / ZRST(:) )**XLBEXS 
+        END WHERE
    ELSE
-      WHERE (ZRST(:)>XRTMIN(5) .AND. ZCST(:)>XCTMIN(5))
-         ZLBDAS(:) = ( XLBS*ZCST(:) / ZRST(:) )**XLBEXS
+      WHERE (ZRST(:)>XRTMIN(5) )
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX,XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS),XLBDAS_MIN)
       END WHERE
+      ZCST(:) = XCCS*ZLBDAS(:)**XCXS / ZRHODREF(:)
+      ZCSS(:) = XCCS*ZLBDAS(:)**XCXS / ZRHODREF(:)
    END IF
    ZLBDAG(:)  = 1.E10
-   WHERE (ZRGT(:)>XRTMIN(6) .AND. ZCGT(:)>XCTMIN(6) )  ! MTaufour
-      ZLBDAG(:) =  ( XLBG*ZCGT(:) / ZRGT(:) )**XLBEXG  ! MTaufour
-   END WHERE
+   if (NMOM_G.GE.2) then
+       WHERE (ZRGT(:)>XRTMIN(6) .AND. ZCGT(:)>XCTMIN(6) )  
+           ZLBDAG(:) =  ( XLBG*ZCGT(:) / ZRGT(:) )**XLBEXG 
+       END WHERE
+   else
+       WHERE (ZRGT(:)>XRTMIN(6) )
+           ZLBDAG(:) = XLBG*( ZRHODREF(:)*ZRGT(:) )**XLBEXG
+       END WHERE
+      ZCGT(:) = XCCG*ZLBDAG(:)**XCXG / ZRHODREF(:)
+      ZCGS(:) = XCCG*ZLBDAG(:)**XCXG / ZRHODREF(:)
+   end if
    ZLBDAH(:)  = 1.E10
-   WHERE (ZRHT(:)>XRTMIN(7) .AND.  ZCHT(:)>XCTMIN(7))  ! MTaufour
-      ZLBDAH(:) =  ( XLBH*ZCHT(:) / ZRHT(:) )**XLBEXH  ! MTaufour
-   END WHERE
+   if (NMOM_H.GE.2) then
+       WHERE (ZRHT(:)>XRTMIN(7) .AND.  ZCHT(:)>XCTMIN(7))  
+           ZLBDAH(:) =  ( XLBH*ZCHT(:) / ZRHT(:) )**XLBEXH 
+       END WHERE
+   else
+       WHERE (ZRHT(:)>XRTMIN(7) )
+           ZLBDAH(:) = XLBH*( ZRHODREF(:)*ZRHT(:) )**XLBEXH
+       END WHERE
+      ZCHT(:) = XCCH*ZLBDAH(:)**XCXH / ZRHODREF(:)
+      ZCHS(:) = XCCH*ZLBDAH(:)**XCXH / ZRHODREF(:)
+   end if
 ! 
 !-------------------------------------------------------------------------------
 !
@@ -571,12 +598,12 @@ IF( IMICRO >= 0 ) THEN
 !
    CALL LIMA_MIXED_SLOW_PROCESSES(ZRHODREF, ZZT, ZSSI, PTSTEP,        &
                                   ZLSFACT, ZLVFACT, ZAI, ZCJ,         &
-                                  ZRGT, ZRHT, ZCIT, ZCGT, ZCHT,       & ! MTaufour add ZCGT, ZCHT
-                                  ZRVS, ZRCS, ZRIS, ZRGS, ZRHS, ZTHS, & ! MTaufour add ZRHS 
-                                  ZCCS, ZCIS, ZCGS, ZIFS, ZINS,       & ! MTaufour add ZCGS
-                                  ZLBDAI, ZLBDAG, ZLBDAH,             & ! MTaufour add ZLBDAH
+                                  ZRGT, ZRHT, ZCIT, ZCGT, ZCHT,       &
+                                  ZRVS, ZRCS, ZRIS, ZRGS, ZRHS, ZTHS, &
+                                  ZCCS, ZCIS, ZCGS, ZIFS, ZINS,       &
+                                  ZLBDAI, ZLBDAG, ZLBDAH,             &
                                   ZRHODJ, GMICRO, PRHODJ, KMI,        &
-                                  PTHS, PRVS, PRCS, PRIS, PRGS, PRHS, & ! MTaufour add PRHS
+                                  PTHS, PRVS, PRCS, PRIS, PRGS, PRHS, &
                                   PCCS, PCIS                          )
 ! 
 !-------------------------------------------------------------------------------
@@ -586,16 +613,16 @@ IF( IMICRO >= 0 ) THEN
 !   	        ------------------------------------
 !
 IF (LSNOW) THEN
-!   CALL LIMA_MIXED_FAST_PROCESSES(ZRHODREF, ZZT, ZPRES, PTSTEP,             &
-!                                  ZLSFACT, ZLVFACT, ZKA, ZDV, ZCJ,          &
-!                                  ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT,       &
-!                                  ZRHT, ZCCT, ZCRT, ZCIT, ZCST, ZCGT, ZCHT, & ! MTaufour add ZCST, ZCGT, ZCHT
-!                                  ZRCS, ZRRS, ZRIS, ZRSS, ZRGS, ZRHS,       &
-!                                  ZTHS, ZCCS, ZCRS, ZCIS, ZCSS, ZCGS, ZCHS, & ! MTaufour add  ZCSS, ZCGS, ZCHS
-!                                  ZLBDAC, ZLBDAR, ZLBDAI, ZLBDAS, ZLBDAG, ZLBDAH,   & ! MTaufour add ZLBDAI
-!                                  ZRHODJ, GMICRO, PRHODJ, KMI, PTHS,        &
-!                                  PRCS, PRRS, PRIS, PRSS, PRGS, PRHS,       &
-!                                  PCCS, PCRS, PCIS, PCSS, PCGS, PCHS        ) ! MTaufour add PCSS, PCGS, PCHS
+   CALL LIMA_MIXED_FAST_PROCESSES(ZRHODREF, ZZT, ZPRES, PTSTEP,                   &
+                                  ZLSFACT, ZLVFACT, ZKA, ZDV, ZCJ,                &
+                                  ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT,             &
+                                  ZRHT, ZCCT, ZCRT, ZCIT, ZCST, ZCGT, ZCHT,       &
+                                  ZRCS, ZRRS, ZRIS, ZRSS, ZRGS, ZRHS,             &
+                                  ZTHS, ZCCS, ZCRS, ZCIS, ZCSS, ZCGS, ZCHS,       &
+                                  ZLBDAC, ZLBDAR, ZLBDAI, ZLBDAS, ZLBDAG, ZLBDAH, &
+                                  ZRHODJ, GMICRO, PRHODJ, KMI, PTHS,              &
+                                  PRCS, PRRS, PRIS, PRSS, PRGS, PRHS,             &
+                                  PCCS, PCRS, PCIS, PCSS, PCGS, PCHS              )
 END IF
 !
 !-------------------------------------------------------------------------------
@@ -630,12 +657,18 @@ END IF
    PCRS(:,:,:) = UNPACK( ZCRS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
    ZW(:,:,:) = PCIS(:,:,:)
    PCIS(:,:,:) = UNPACK( ZCIS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
-   ZW(:,:,:) = PCSS(:,:,:)                                            ! MTaufour
-   PCSS(:,:,:) = UNPACK( ZCSS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) ) ! MTaufour
-   ZW(:,:,:) = PCGS(:,:,:)                                            ! MTaufour
-   PCGS(:,:,:) = UNPACK( ZCGS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) ) ! MTaufour
-   ZW(:,:,:) = PCHS(:,:,:)                                            ! MTaufour
-   PCHS(:,:,:) = UNPACK( ZCHS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) ) ! MTaufour
+   if (NMOM_S.GE.2) then
+       ZW(:,:,:) = PCSS(:,:,:)
+       PCSS(:,:,:) = UNPACK( ZCSS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) ) 
+   end if
+   if (NMOM_G.GE.2) then
+       ZW(:,:,:) = PCGS(:,:,:)
+       PCGS(:,:,:) = UNPACK( ZCGS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) )
+   end if
+   if (NMOM_H.GE.2) then
+       ZW(:,:,:) = PCHS(:,:,:)      
+       PCHS(:,:,:) = UNPACK( ZCHS(:),MASK=GMICRO(:,:,:),FIELD=ZW(:,:,:) ) 
+   end if
 !
    DO JMOD_IFN = 1, NMOD_IFN
       ZW(:,:,:) = PIFS(:,:,:,JMOD_IFN)
@@ -657,9 +690,9 @@ END IF
    DEALLOCATE(ZCCT) 
    DEALLOCATE(ZCRT) 
    DEALLOCATE(ZCIT)
-   DEALLOCATE(ZCST) ! MTaufour
-   DEALLOCATE(ZCGT) ! MTaufour
-   DEALLOCATE(ZCHT) ! MTaufour   
+   DEALLOCATE(ZCST) 
+   DEALLOCATE(ZCGT)
+   DEALLOCATE(ZCHT) 
 ! 
    DEALLOCATE(ZRVS)  
    DEALLOCATE(ZRCS) 
@@ -673,9 +706,9 @@ END IF
    DEALLOCATE(ZCCS) 
    DEALLOCATE(ZCRS) 
    DEALLOCATE(ZCIS)  
-   DEALLOCATE(ZCSS) ! MTaufour  
-   DEALLOCATE(ZCGS) ! MTaufour     
-   DEALLOCATE(ZCHS) ! MTaufour        
+   DEALLOCATE(ZCSS)  
+   DEALLOCATE(ZCGS)     
+   DEALLOCATE(ZCHS)
    DEALLOCATE(ZIFS)
    DEALLOCATE(ZINS)
 !
@@ -722,9 +755,9 @@ IF ( KRR .GE. 7 ) PRS(:,:,:,7) = PRHS(:,:,:)
 !
 PSVS(:,:,:,NSV_LIMA_NC) = PCCS(:,:,:)
 IF ( LRAIN ) PSVS(:,:,:,NSV_LIMA_NR) = PCRS(:,:,:)
-IF ( LSNOW ) PSVS(:,:,:,NSV_LIMA_NS) = PCSS(:,:,:) ! MTaufour
-IF ( LSNOW ) PSVS(:,:,:,NSV_LIMA_NG) = PCGS(:,:,:) ! MTaufour
-IF ( LHAIL ) PSVS(:,:,:,NSV_LIMA_NH) = PCHS(:,:,:) ! MTaufour
+IF ( NMOM_S.GE.2 ) PSVS(:,:,:,NSV_LIMA_NS) = PCSS(:,:,:)
+IF ( NMOM_G.GE.2 ) PSVS(:,:,:,NSV_LIMA_NG) = PCGS(:,:,:)
+IF ( NMOM_H.GE.2 ) PSVS(:,:,:,NSV_LIMA_NH) = PCHS(:,:,:)
 PSVS(:,:,:,NSV_LIMA_NI) = PCIS(:,:,:)
 !
 IF ( NMOD_CCN .GE. 1 ) THEN
diff --git a/src/MNH/lima_mixed_fast_processes.f90 b/src/MNH/lima_mixed_fast_processes.f90
index 056386015b61db1bb4f7b0292e01bebeacae45cf..f05468c3c6af494ee5e6fc1f318e80a879040336 100644
--- a/src/MNH/lima_mixed_fast_processes.f90
+++ b/src/MNH/lima_mixed_fast_processes.f90
@@ -8,16 +8,16 @@
 !      #####################################
 !
 INTERFACE
-      SUBROUTINE LIMA_MIXED_FAST_PROCESSES (PRHODREF, PZT, PPRES, PTSTEP,           &
-                                            PLSFACT, PLVFACT, PKA, PDV, PCJ,        &
-                                            PRVT1D, PRCT1D, PRRT1D, PRIT1D, PRST1D, &
-                                            PRGT1D, PRTH1D, PCCT1D, PCRT1D, PCIT1D, &
-                                            PRCS1D, PRRS1D, PRIS1D, PRSS1D, PRGS1D, &
-                                            PRHS1D, PTHS1D, PCCS1D, PCRS1D, PCIS1D, &
-                                            PLBDAC, PLBDAR, PLBDAS, PLBDAG, PLBDAH, &
-                                            PRHODJ1D, GMICRO, PRHODJ, KMI, PTHS,    &
-                                            PRCS, PRRS, PRIS, PRSS, PRGS, PRHS,     &
-                                            PCCS, PCRS, PCIS                        )
+      SUBROUTINE LIMA_MIXED_FAST_PROCESSES (PRHODREF, PZT, PPRES, PTSTEP,                   &
+                                            PLSFACT, PLVFACT, PKA, PDV, PCJ,                &
+                                            PRVT1D, PRCT1D, PRRT1D, PRIT1D, PRST1D, PRGT1D, &
+                                            PRHT1D, PCCT1D, PCRT1D, PCIT1D, PCST1D, PCGT1D, PCHT1D, &
+                                            PRCS1D, PRRS1D, PRIS1D, PRSS1D, PRGS1D, PRHS1D, &
+                                            PTHS1D, PCCS1D, PCRS1D, PCIS1D, PCSS1D, PCGS1D, PCHS1D, &
+                                            PLBDAC, PLBDAR, PLBDAI, PLBDAS, PLBDAG, PLBDAH, &
+                                            PRHODJ1D, GMICRO, PRHODJ, KMI, PTHS,            &
+                                            PRCS, PRRS, PRIS, PRSS, PRGS, PRHS,             &
+                                            PCCS, PCRS, PCIS, PCSS, PCGS, PCHS              )
 !
 REAL, DIMENSION(:),   INTENT(IN)    :: PRHODREF  ! RHO Dry REFerence
 REAL, DIMENSION(:),   INTENT(IN)    :: PZT       ! Temperature
@@ -36,11 +36,14 @@ REAL, DIMENSION(:),   INTENT(IN)    :: PRRT1D    ! Rain water m.r. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PRIT1D    ! Pristine ice m.r. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PRST1D    ! Snow/aggregate m.r. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PRGT1D    ! Graupel m.r. at t
-REAL, DIMENSION(:),   INTENT(IN)    :: PRTH1D    ! Hail m.r. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PRHT1D    ! Hail m.r. at t
 !
 REAL, DIMENSION(:),   INTENT(IN)    :: PCCT1D    ! Cloud water conc. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PCRT1D    ! Rain water conc. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PCIT1D    ! Pristine ice conc. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PCST1D    ! Snow/aggregate conc. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PCGT1D    ! Graupel conc. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PCHT1D    ! Hail conc. at t
 !
 REAL, DIMENSION(:),   INTENT(INOUT) :: PRCS1D    ! Cloud water m.r. source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PRRS1D    ! Rain water m.r. source
@@ -54,9 +57,13 @@ REAL, DIMENSION(:),   INTENT(INOUT) :: PTHS1D    ! Theta source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PCCS1D    ! Cloud water conc. source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PCRS1D    ! Rain water conc. source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PCIS1D    ! Pristine ice conc. source
+REAL, DIMENSION(:),   INTENT(INOUT) :: PCSS1D    ! Snow/aggregate conc. source
+REAL, DIMENSION(:),   INTENT(INOUT) :: PCGS1D    ! Graupel conc. source
+REAL, DIMENSION(:),   INTENT(INOUT) :: PCHS1D    ! Hail conc. source
 !
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAC  ! Slope param of the cloud droplet distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAR  ! Slope param of the raindrop  distr
+REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAI  ! Slope param of the ice distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAS  ! Slope param of the aggregate distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAG  ! Slope param of the graupel distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAH  ! Slope param of the hail distr.
@@ -76,22 +83,25 @@ REAL,    DIMENSION(:,:,:), INTENT(IN) :: PRHS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCCS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCRS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCIS
+REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCSS 
+REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCGS 
+REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCHS
 !
 END SUBROUTINE LIMA_MIXED_FAST_PROCESSES
 END INTERFACE
 END MODULE MODI_LIMA_MIXED_FAST_PROCESSES
 !
 !     ###############################################################################
-      SUBROUTINE LIMA_MIXED_FAST_PROCESSES (PRHODREF, PZT, PPRES, PTSTEP,           &
-                                            PLSFACT, PLVFACT, PKA, PDV, PCJ,        &
-                                            PRVT1D, PRCT1D, PRRT1D, PRIT1D, PRST1D, &
-                                            PRGT1D, PRTH1D, PCCT1D, PCRT1D, PCIT1D, &
-                                            PRCS1D, PRRS1D, PRIS1D, PRSS1D, PRGS1D, &
-                                            PRHS1D, PTHS1D, PCCS1D, PCRS1D, PCIS1D, &
-                                            PLBDAC, PLBDAR, PLBDAS, PLBDAG, PLBDAH, &
-                                            PRHODJ1D, GMICRO, PRHODJ, KMI, PTHS,    &
-                                            PRCS, PRRS, PRIS, PRSS, PRGS, PRHS,     &
-                                            PCCS, PCRS, PCIS                        )
+      SUBROUTINE LIMA_MIXED_FAST_PROCESSES (PRHODREF, PZT, PPRES, PTSTEP,                   &
+                                            PLSFACT, PLVFACT, PKA, PDV, PCJ,                &
+                                            PRVT1D, PRCT1D, PRRT1D, PRIT1D, PRST1D, PRGT1D, &
+                                            PRHT1D, PCCT1D, PCRT1D, PCIT1D, PCST1D, PCGT1D, PCHT1D, &
+                                            PRCS1D, PRRS1D, PRIS1D, PRSS1D, PRGS1D, PRHS1D, &
+                                            PTHS1D, PCCS1D, PCRS1D, PCIS1D, PCSS1D, PCGS1D, PCHS1D, &
+                                            PLBDAC, PLBDAR, PLBDAI, PLBDAS, PLBDAG, PLBDAH, &
+                                            PRHODJ1D, GMICRO, PRHODJ, KMI, PTHS,            &
+                                            PRCS, PRRS, PRIS, PRSS, PRGS, PRHS,             &
+                                            PCCS, PCRS, PCIS, PCSS, PCGS, PCHS              )
 !     ###############################################################################
 !
 !!
@@ -146,6 +156,7 @@ END MODULE MODI_LIMA_MIXED_FAST_PROCESSES
 !                          - change the name of some arguments to match the DOCTOR norm
 !                          - change conditions for HMG to occur
 !  J. Wurtz       03/2022: new snow characteristics
+!  M. Taufour     07/2022: add concentration for snow, graupel, hail
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -184,12 +195,15 @@ REAL, DIMENSION(:),   INTENT(IN)    :: PRCT1D    ! Cloud water m.r. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PRRT1D    ! Rain water m.r. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PRIT1D    ! Pristine ice m.r. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PRST1D    ! Snow/aggregate m.r. at t
-REAL, DIMENSION(:),   INTENT(IN)    :: PRGT1D    ! Graupel/hail m.r. at t
-REAL, DIMENSION(:),   INTENT(IN)    :: PRTH1D    ! Hail m.r. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PRGT1D    ! Graupel m.r. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PRHT1D    ! Hail m.r. at t
 !
 REAL, DIMENSION(:),   INTENT(IN)    :: PCCT1D    ! Cloud water conc. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PCRT1D    ! Rain water conc. at t
 REAL, DIMENSION(:),   INTENT(IN)    :: PCIT1D    ! Pristine ice conc. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PCST1D    ! Snow/aggregate conc. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PCGT1D    ! Graupel conc. at t
+REAL, DIMENSION(:),   INTENT(IN)    :: PCHT1D    ! Hail conc. at t
 !
 REAL, DIMENSION(:),   INTENT(INOUT) :: PRCS1D    ! Cloud water m.r. source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PRRS1D    ! Rain water m.r. source
@@ -203,9 +217,13 @@ REAL, DIMENSION(:),   INTENT(INOUT) :: PTHS1D    ! Theta source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PCCS1D    ! Cloud water conc. source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PCRS1D    ! Rain water conc. source
 REAL, DIMENSION(:),   INTENT(INOUT) :: PCIS1D    ! Pristine ice conc. source
+REAL, DIMENSION(:),   INTENT(INOUT) :: PCSS1D    ! Snow/aggregate conc. source
+REAL, DIMENSION(:),   INTENT(INOUT) :: PCGS1D    ! Graupel conc. source
+REAL, DIMENSION(:),   INTENT(INOUT) :: PCHS1D    ! Hail conc. source
 !
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAC  ! Slope param of the cloud droplet distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAR  ! Slope param of the raindrop  distr
+REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAI  ! Slope param of the ice distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAS  ! Slope param of the aggregate distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAG  ! Slope param of the graupel distr.
 REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAH  ! Slope param of the hail distr.
@@ -214,7 +232,7 @@ REAL, DIMENSION(:),   INTENT(IN)    :: PLBDAH  ! Slope param of the hail distr.
 REAL,    DIMENSION(:),     INTENT(IN) :: PRHODJ1D
 LOGICAL, DIMENSION(:,:,:), INTENT(IN) :: GMICRO 
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PRHODJ
-INTEGER,                   INTENT(IN) :: KMI
+INTEGER,                   INTENT(IN) :: KMI 
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PTHS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PRCS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PRRS
@@ -225,7 +243,9 @@ REAL,    DIMENSION(:,:,:), INTENT(IN) :: PRHS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCCS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCRS
 REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCIS
-
+REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCSS 
+REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCGS 
+REAL,    DIMENSION(:,:,:), INTENT(IN) :: PCHS
 !
 !*       0.2   Declarations of local variables :
 !
@@ -234,9 +254,10 @@ INTEGER :: IGRIM, IGACC, IGDRY, IGWET, IHAIL
 INTEGER :: JJ
 INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC1,IVEC2        ! Vectors of indices
 REAL,    DIMENSION(:), ALLOCATABLE :: ZVEC1,ZVEC2, ZVEC3 ! Work vectors
-REAL,    DIMENSION(SIZE(PZT))  :: ZZW, ZZX      
-REAL,    DIMENSION(SIZE(PZT))  :: ZRDRYG, ZRWETG   
-REAL,    DIMENSION(SIZE(PZT),7)  :: ZZW1 
+REAL,    DIMENSION(SIZE(PZT))  :: ZZW, ZZX, ZZNW
+REAL,    DIMENSION(SIZE(PRIT1D)) :: ZZDI, ZZDC
+REAL,    DIMENSION(SIZE(PZT))  :: ZRDRYG, ZRWETG, ZNWETG
+REAL,    DIMENSION(SIZE(PZT),7)  :: ZZW1, ZZNW1 
 REAL :: NHAIL
 REAL :: ZTHRH, ZTHRC
 !
@@ -282,10 +303,13 @@ REAL,    DIMENSION(:), ALLOCATABLE :: ZNI_RDSF,ZRI_RDSF  ! RDSF rates
 REAL,    DIMENSION(:),   ALLOCATABLE :: ZAUX     ! used to distribute
 REAL,    DIMENSION(:,:), ALLOCATABLE :: ZFACT    ! the total concentration in each shape
 REAL,    DIMENSION(:),   ALLOCATABLE :: ZONEOVER_VAR ! for optimization
+LOGICAL :: M2_ICE
 !
 !
 !-------------------------------------------------------------------------------
 !
+M2_ICE = NMOM_S.GE.2 .AND. NMOM_G.GE.2 .AND. NMOM_H.GE.2
+!
 !                         #################
 !                         FAST RS PROCESSES
 !                         #################
@@ -299,6 +323,7 @@ SNOW: IF (LSNOW) THEN
 ZZW1(:,:) = 0.0
 !
 GRIM(:) = (PRCT1D(:)>XRTMIN(2)) .AND. (PRST1D(:)>XRTMIN(5)) .AND. (PRCS1D(:)>XRTMIN(2)/PTSTEP) .AND. (PZT(:)<XTT)
+IF (NMOM_S.GE.2) GRIM(:) = GRIM(:) .AND. (PCST1D(:)>XCTMIN(5))
 IGRIM = COUNT( GRIM(:) )
 !
 IF( IGRIM>0 ) THEN
@@ -313,6 +338,10 @@ IF( IGRIM>0 ) THEN
                                             Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
     if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nc), 'RIM', &
                                             Unpack( pccs1d(:), mask = gmicro(:, :, :), field = pccs(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_S.GE.2 ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'RIM', &
+                                            Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_G.GE.2 ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'RIM', &
+                                            Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
   end if
 !
 !        1.1.0  allocations
@@ -345,16 +374,16 @@ IF( IGRIM>0 ) THEN
 !        1.1.4  riming of the small sized aggregates
 !
   WHERE ( GRIM(:) )
-    ZZW1(:,1) = MIN( PRCS1D(:),                              &
-                 XCRIMSS * ZZW(:) * PRCT1D(:)                & ! RCRIMSS
-                                  * PRST1D(:)*(1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
-                                      * PRHODREF(:)**(-XCEXVT+1) &
-                                      * (PLBDAS(:)) ** (XEXCRIMSS+XBS) )
-    PRCS1D(:) = PRCS1D(:) - ZZW1(:,1)
-    PRSS1D(:) = PRSS1D(:) + ZZW1(:,1)
-    PTHS1D(:) = PTHS1D(:) + ZZW1(:,1) * (PLSFACT(:) - PLVFACT(:)) ! f(L_f*(RCRIMSS))
-!
-    PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,1)*(PCCT1D(:)/PRCT1D(:)),0.0 ) ! Lambda_c**3
+     ZZW1(:,1) = MIN( PRCS1D(:),                          &
+  	           XCRIMSS * ZZW(:) * PRCT1D(:) * PCST1D(:)      & 
+  	                            *   PLBDAS(:)**XEXCRIMSS &
+               * (1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
+    			            * PRHODREF(:)**(-XCEXVT) )
+     PRCS1D(:) = PRCS1D(:) - ZZW1(:,1)
+     PRSS1D(:) = PRSS1D(:) + ZZW1(:,1)
+     PTHS1D(:) = PTHS1D(:) + ZZW1(:,1)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(RCRIMSS))
+!
+     PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,1)*(PCCT1D(:)/PRCT1D(:)),0.0 ) ! Lambda_c**3
   END WHERE
 !
 !        1.1.5  perform the linear interpolation of the normalized
@@ -366,20 +395,23 @@ IF( IGRIM>0 ) THEN
 !
 !        1.1.6  riming-conversion of the large sized aggregates into graupeln
 !
-  WHERE ( GRIM(:) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP) )
-    ZZW1(:,2) = MIN( PRCS1D(:),                     &
-                 XCRIMSG * PRCT1D(:)*  PRST1D(:)                & ! RCRIMSG
-                           *(1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSG/XALPHAS)*PLBDAS(:)**(XBS+XEXCRIMSG)  &
-                           * PRHODREF(:)**(-XCEXVT+1) &
-                           - ZZW1(:,1)              )
-    ZZW1(:,3) = MIN( PRSS1D(:),                         &
-                     XSRIMCG * XNS * PRST1D(:) * (1.0 - ZZW(:))/PTSTEP )
-    PRCS1D(:) = PRCS1D(:) - ZZW1(:,2)
-    PRSS1D(:) = PRSS1D(:) - ZZW1(:,3)
-    PRGS1D(:) = PRGS1D(:) + ZZW1(:,2) + ZZW1(:,3)
-    PTHS1D(:) = PTHS1D(:) + ZZW1(:,2) * (PLSFACT(:) - PLVFACT(:)) ! f(L_f*(RCRIMSG))
-    !
-    PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,2)*(PCCT1D(:)/PRCT1D(:)),0.0 ) ! Lambda_c**3
+  WHERE ( GRIM(:) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP) .AND. (PCSS1D(:)>XCTMIN(5)/PTSTEP))
+     ZZW1(:,2) = MIN( PRCS1D(:),                 &
+    	           XCRIMSG * PRCT1D(:) * PCST1D(:)      & ! RCRIMSG
+    	           *  PLBDAS(:)**XEXCRIMSG*(1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSG/XALPHAS)  &
+  	                   * PRHODREF(:)**(-XCEXVT) &
+    		           - ZZW1(:,1)              )
+     ZZW1(:,3) = MIN( PRSS1D(:),  PCST1D(:) *          &
+                       XSRIMCG * PLBDAS(:)**XEXSRIMCG   & ! RSRIMCG 
+   	                       * (1.0 - ZZW(:) )/(PTSTEP*PRHODREF(:)))
+     PRCS1D(:) = PRCS1D(:) - ZZW1(:,2)
+     PRSS1D(:) = PRSS1D(:) - ZZW1(:,3)
+     PRGS1D(:) = PRGS1D(:) + ZZW1(:,2) + ZZW1(:,3)
+     PTHS1D(:) = PTHS1D(:) + ZZW1(:,2)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(RCRIMSG))
+!
+     PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,2)*(PCCT1D(:)/PRCT1D(:)),0.0 ) ! Lambda_c**3
+     PCSS1D(:) = MAX( PCSS1D(:)-ZZW1(:,3)*(PCST1D(:)/PRST1D(:)),0.0 )
+     PCGS1D(:) = MAX( PCGS1D(:)+ZZW1(:,3)*(PCST1D(:)/PRST1D(:)),0.0 ) !
   END WHERE
   DEALLOCATE(IVEC2)
   DEALLOCATE(IVEC1)
@@ -398,6 +430,10 @@ IF( IGRIM>0 ) THEN
                                            Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
     if ( lbudget_sv ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nc), 'RIM', &
                                            Unpack( pccs1d(:), mask = gmicro(:, :, :), field = pccs(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_S.GE.2) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'RIM', &
+                                           Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_G.GE.2) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'RIM', &
+                                           Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
   end if
 END IF
 !
@@ -709,7 +745,7 @@ IF (ICIBU > 0) THEN
 ! The value of rs removed by CIBU is determined by the mean mass of pristine ice
   WHERE( PRIT1D(:)>XRTMIN(4) .AND. PCIT1D(:)>XCTMIN(4) )
     ZRI_CIBU(:) = MIN( ZRI_CIBU(:), PRSS1D(:), ZNI_CIBU(:)*PRIT1D(:)/PCIT1D(:) )
-  ELSE WHERE
+  ELSEWHERE
     ZRI_CIBU(:) = MIN( ZRI_CIBU(:), PRSS1D(:), MAX( ZNI_CIBU(:)*XMNU0,XRTMIN(4) ) )
   END WHERE
 !
@@ -780,6 +816,10 @@ IF( IGACC>0 .AND. LRAIN) THEN
                                             Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
     if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'ACC', &
                                             Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_S.GE.2) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'ACC', &
+                                            Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_G.GE.2) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'ACC', &
+                                            Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
   end if
 !
 !        1.4.0  allocations
@@ -792,7 +832,11 @@ IF( IGACC>0 .AND. LRAIN) THEN
 !
 !        1.4.1  select the (PLBDAS,PLBDAR) couplet
 !
-  ZVEC1(:) = PACK( PLBDAS(:),MASK=GACC(:) )
+  if (M2_ICE) then
+     ZVEC1(:) = PACK( MAX(MIN(PLBDAS(:),5.E5),5.E1),MASK=GACC(:) )
+  else
+     ZVEC1(:) = PACK( PLBDAS(:),MASK=GACC(:) )
+  end if
   ZVEC2(:) = PACK( PLBDAR(:),MASK=GACC(:) )
 !
 !        1.4.2  find the next lower indice for the PLBDAS and for the PLBDAR
@@ -821,21 +865,44 @@ IF( IGACC>0 .AND. LRAIN) THEN
                                                           * (ZVEC1(JJ) - 1.0)
   END DO
   ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
+  DEALLOCATE(ZVEC3)
+  ALLOCATE(ZVEC3(IGACC))  
+!
+!                                                                         
+!        1.4.3b  perform the bilinear interpolation of the normalized
+!               RACCSS-kernel FOR CONCENTRATION
+!
+  DO JJ = 1,IGACC
+     ZVEC3(JJ) =  (  XKER_N_RACCSS(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                    - XKER_N_RACCSS(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                				 	           * ZVEC1(JJ) &
+                 - (  XKER_N_RACCSS(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                    - XKER_N_RACCSS(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+  	                    			             * (ZVEC1(JJ) - 1.0)
+  END DO
+  ZZNW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
 !
 !        1.4.4  raindrop accretion on the small sized aggregates
 !
   WHERE ( GACC(:) )
-    ZZW1(:,2) = PCRT1D(:) *                                           & !! coef of RRACCS
-              XFRACCSS*( PRST1D(:)*PLBDAS(:)**XBS )*( PRHODREF(:)**(1-XCEXVT) ) &
-         *( XLBRACCS1/((PLBDAS(:)**2)               ) +                  &
-            XLBRACCS2/( PLBDAS(:)    * PLBDAR(:)    ) +                  &
-            XLBRACCS3/(               (PLBDAR(:)**2)) )/PLBDAR(:)**XBR
-    ZZW1(:,4) = MIN( PRRS1D(:),ZZW1(:,2)*ZZW(:) )           ! RRACCSS
-    PRRS1D(:) = PRRS1D(:) - ZZW1(:,4)
-    PRSS1D(:) = PRSS1D(:) + ZZW1(:,4)
-    PTHS1D(:) = PTHS1D(:) + ZZW1(:,4) * (PLSFACT(:) - PLVFACT(:)) ! f(L_f*(RRACCSS))
-!
-    PCRS1D(:) = MAX( PCRS1D(:)-ZZW1(:,4)*(PCRT1D(:)/PRRT1D(:)),0.0 ) ! Lambda_r**3 
+     ZZW1(:,2) = PCRT1D(:) *                                      & !! coef of RRACCS 
+                  XFRACCSS*( PCST1D(:) )*( PRHODREF(:)**(-XCEXVT-1.) ) &
+             *( XLBRACCS1/((PLBDAS(:)**2)               ) +                  &
+                XLBRACCS2/( PLBDAS(:)    * PLBDAR(:)    ) +                  &
+                XLBRACCS3/(               (PLBDAR(:)**2)) )/PLBDAR(:)**3
+!                                                                                
+     ZZNW1(:,2) = PCRT1D(:) *                                           & !! coef of CRACCS
+                  XFNRACCSS*( PCST1D(:) )*( PRHODREF(:)**(-XCEXVT-1.) ) &
+             *( XLBNRACCS1/((PLBDAS(:)**2)               ) +                  &
+               XLBNRACCS2/( PLBDAS(:)    * PLBDAR(:)    ) +                  &
+                XLBNRACCS3/(               (PLBDAR(:)**2)) )
+     ZZW1(:,4) = MIN( PRRS1D(:),ZZW1(:,2)*ZZW(:) )           ! RRACCSS
+     PRRS1D(:) = PRRS1D(:) - ZZW1(:,4)
+     PRSS1D(:) = PRSS1D(:) + ZZW1(:,4)
+     PTHS1D(:) = PTHS1D(:) + ZZW1(:,4)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(RRACCSS))
+!
+     ZZNW1(:,4) = MIN( PCRS1D(:),ZZNW1(:,2)*ZZNW(:) )           ! NRACCSS      
+     PCRS1D(:) = PCRS1D(:)-ZZNW1(:,4)
   END WHERE
 !
 !        1.4.4b perform the bilinear interpolation of the normalized
@@ -850,6 +917,19 @@ IF( IGACC>0 .AND. LRAIN) THEN
                                                          * (ZVEC1(JJ) - 1.0)
   END DO
   ZZW1(:,2) = ZZW1(:,2)*UNPACK( VECTOR=ZVEC3(:),MASK=GACC(:),FIELD=0.0 ) !! RRACCS
+!                                                                                 
+!        1.3.4b2 perform the bilinear interpolation of the normalized
+!               RACCS-kernel
+!
+  DO JJ = 1,IGACC
+     ZVEC3(JJ) =  (   XKER_N_RACCS(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                         -  XKER_N_RACCS(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                                                                         * ZVEC1(JJ) &
+                         - (   XKER_N_RACCS(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                         -  XKER_N_RACCS(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                                                           * (ZVEC1(JJ) - 1.0)
+  END DO
+  ZZNW1(:,2) = ZZNW1(:,2)*UNPACK( VECTOR=ZVEC3(:),MASK=GACC(:),FIELD=0.0 ) !! NRACCS
 !
 !        1.4.5  perform the bilinear interpolation of the normalized
 !               SACCRG-kernel
@@ -864,28 +944,49 @@ IF( IGACC>0 .AND. LRAIN) THEN
   END DO
   ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
 !
+!        1.3.5b  perform the bilinear interpolation of the normalized
+!               SACCRG-kernel
+!
+  DO JJ = 1,IGACC
+     ZVEC3(JJ) =  (  XKER_N_SACCRG(IVEC2(JJ)+1,IVEC1(JJ)+1)* ZVEC1(JJ)          &
+                        - XKER_N_SACCRG(IVEC2(JJ)+1,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
+      			 	                                   * ZVEC2(JJ) &
+                     - (  XKER_N_SACCRG(IVEC2(JJ)  ,IVEC1(JJ)+1)* ZVEC1(JJ)          &
+                        - XKER_N_SACCRG(IVEC2(JJ)  ,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
+			                                     * (ZVEC2(JJ) - 1.0)
+  END DO
+  ZZNW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
+!
 !        1.4.6  raindrop accretion-conversion of the large sized aggregates
 !               into graupeln
 !
-  WHERE ( GACC(:) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP) )
-    ZZW1(:,2) = MAX( MIN( PRRS1D(:),ZZW1(:,2)-ZZW1(:,4) ) , 0. )      ! RRACCSG
-    ZZW1(:,3) = MIN( PRSS1D(:),XFSACCRG*ZZW(:)*                     & ! RSACCRG
-            ( PRST1D(:) )*( PRHODREF(:)**(1-XCEXVT) ) &
-           *( XLBSACCR1/((PLBDAR(:)**2)               ) +           &
-              XLBSACCR2/( PLBDAR(:)    * PLBDAS(:)    ) +           &
-              XLBSACCR3/(               (PLBDAS(:)**2)) ) )
-    PRRS1D(:) = PRRS1D(:) - ZZW1(:,2)
-    PRSS1D(:) = PRSS1D(:) - ZZW1(:,3)
-    PRGS1D(:) = PRGS1D(:) + ZZW1(:,2) + ZZW1(:,3)
-    PTHS1D(:) = PTHS1D(:) + ZZW1(:,2) * (PLSFACT(:) - PLVFACT(:)) ! f(L_f*(RRACCSG))
-!
-    PCRS1D(:) = MAX( PCRS1D(:)-ZZW1(:,2)*(PCRT1D(:)/PRRT1D(:)),0.0 ) ! Lambda_r**3 
+  WHERE ( GACC(:) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP)  .AND. (PCSS1D(:)>XCTMIN(5)/PTSTEP) )
+     ZZW1(:,2) = MAX( MIN( PRRS1D(:),ZZW1(:,2)-ZZW1(:,4) ) , 0. )      ! RRACCSG
+     ZZNW1(:,2) = MAX( MIN( PCRS1D(:),ZZNW1(:,2)-ZZNW1(:,4) ) , 0. )   ! NRACCSG  
+     ZZW1(:,3) = MIN( PRSS1D(:),XFSACCRG*ZZW(:)* PCST1D(:) *           & ! RSACCRG 
+                ( PLBDAS(:)**(-XBS) )*( PRHODREF(:)**(-XCEXVT-1.) )     & 
+                *( XLBSACCR1/((PLBDAR(:)**2)               ) +           &
+                  XLBSACCR2/( PLBDAR(:)    * PLBDAS(:)    ) +           &
+                  XLBSACCR3/(               (PLBDAS(:)**2)) ) )
+     ZZNW1(:,3) = MIN( PCSS1D(:),XFNSACCRG*ZZNW(:)* PCST1D(:) *        & ! NSACCRG 
+                                      ( PRHODREF(:)**(-XCEXVT-1.) )     &            
+               *( XLBNSACCR1/((PLBDAR(:)**2)               ) +          &            
+                  XLBNSACCR2/( PLBDAR(:)    * PLBDAS(:)    ) +          &            
+                  XLBNSACCR3/(               (PLBDAS(:)**2)) ) )           
+     PRRS1D(:) = PRRS1D(:) - ZZW1(:,2)
+     PRSS1D(:) = PRSS1D(:) - ZZW1(:,3)
+     PRGS1D(:) = PRGS1D(:) + ZZW1(:,2)+ZZW1(:,3)
+     PTHS1D(:) = PTHS1D(:) + ZZW1(:,2)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(RRACCSG))
+!
+     PCRS1D(:) = PCRS1D(:)-ZZNW1(:,2) !                                   
+     PCSS1D(:) = PCSS1D(:)-ZZNW1(:,3)
+     PCGS1D(:) = PCGS1D(:)+ZZNW1(:,3)
   END WHERE
-  DEALLOCATE(IVEC2)
-  DEALLOCATE(IVEC1)
-  DEALLOCATE(ZVEC3)
-  DEALLOCATE(ZVEC2)
-  DEALLOCATE(ZVEC1)
+   DEALLOCATE(IVEC2)
+   DEALLOCATE(IVEC1)
+   DEALLOCATE(ZVEC3)
+   DEALLOCATE(ZVEC2)
+   DEALLOCATE(ZVEC1)
   !
   ! Budget storage
   if ( nbumod == kmi .and. lbu_enable ) then
@@ -899,6 +1000,10 @@ IF( IGACC>0 .AND. LRAIN) THEN
                                            Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
     if ( lbudget_sv ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'ACC', &
                                            Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_S.GE.2) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'ACC', &
+                                           Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+    if ( lbudget_sv .and. NMOM_G.GE.2) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'ACC', &
+                                           Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
   end if
 END IF
 !
@@ -911,30 +1016,34 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                           Unpack( prss1d(:), mask = gmicro(:, :, :), field = prss(:, :, :) ) * prhodj(:, :, :) )
   if ( lbudget_rg ) call Budget_store_init( tbudgets(NBUDGET_RG), 'CMEL', &
                                           Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
+  if ( lbudget_sv .and. NMOM_S.GE.2 .and. NMOM_G.GE.2 ) then
+                   call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'CMEL', &
+                                          Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+                  call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'CMEL', &
+                                          Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) ) 
+  end if
 end if
 !
 ZZW(:) = 0.0
 WHERE( (PRST1D(:)>XRTMIN(5)) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP) .AND. (PZT(:)>XTT) )
-  ZZW(:) = PRVT1D(:) * PPRES(:) / ((XMV / XMD) + PRVT1D(:)) ! Vapor pressure
-  ZZW(:) = PKA(:) * (XTT - PZT(:)) +                            &
-         ( PDV(:) * (XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
-                  * (XESTT-ZZW(:))/(XRV*PZT(:)) )
+   ZZW(:) = PRVT1D(:)*PPRES(:)/((XMV/XMD)+PRVT1D(:)) ! Vapor pressure
+   ZZW(:) =  PKA(:)*(XTT-PZT(:)) +                                 &
+              ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
+                          *(XESTT-ZZW(:))/(XRV*PZT(:))             )
 !
 ! compute RSMLT
 !
-  ZZW(:)  = MIN( PRSS1D(:), XFSCVMG*MAX( 0.0,( -ZZW(:) *             &
-                          PRST1D(:)*( X0DEPS*            PLBDAS(:)**XEX0DEPS +     &
-                          X1DEPS*PCJ(:)*(1+0.5*(XFVELOS/PLBDAS(:))**XALPHAS) &
-                          **(-XNUS+XEX1DEPS/XALPHAS)*(PLBDAS(:))**(XEX1DEPS+XBS))-    &
+   ZZW(:)  = MIN( PRSS1D(:), XFSCVMG*MAX( 0.0,( -ZZW(:) * PCST1D(:) *   & 
+                          ( X0DEPS*       PLBDAS(:)**XEX0DEPS +     &
+                            X1DEPS*PCJ(:)*PLBDAS(:)**XEX1DEPS *     &
+                                  (1+0.5*(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) ) -   &
                                     ( ZZW1(:,1)+ZZW1(:,4) ) *       &
                              ( PRHODREF(:)*XCL*(XTT-PZT(:))) ) /    &
                                             ( 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!!!
-!
-  PRSS1D(:) = PRSS1D(:) - ZZW(:)
-  PRGS1D(:) = PRGS1D(:) + ZZW(:)
+   PRSS1D(:) = PRSS1D(:) - ZZW(:)
+   PRGS1D(:) = PRGS1D(:) + ZZW(:)
+   PCSS1D(:) = MAX( PCSS1D(:) - ZZW(:)*(MAX(PCST1D(:),XCTMIN(5))/MAX(PRST1D(:),XRTMIN(5))), 0.0 )
+   PCGS1D(:) = MAX( PCGS1D(:) + ZZW(:)*(MAX(PCST1D(:),XCTMIN(5))/MAX(PRST1D(:),XRTMIN(5))), 0.0 )
 END WHERE
 !
 ! Budget storage
@@ -943,6 +1052,12 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                          Unpack( prss1d(:), mask = gmicro(:, :, :), field = prss(:, :, :) ) * prhodj(:, :, :) )
   if ( lbudget_rg ) call Budget_store_end( tbudgets(NBUDGET_RG), 'CMEL', &
                                          Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
+  if ( lbudget_sv .and. NMOM_S.GE.2 .AND. NMOM_G.GE.2) then
+                   call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'CMEL', &
+                                         Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+                   call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'CMEL', &
+                                         Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) ) 
+  end if
 end if
 
 END IF SNOW
@@ -970,24 +1085,27 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                          Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
   if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'CFRZ', &
                                          Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+  if ( lbudget_sv .and. NMOM_G.GE.2 ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'CFRZ', &
+                                         Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
 end if
 
 ZZW1(:,3:4) = 0.0
 WHERE( (PRIT1D(:)>XRTMIN(4)) .AND. (PRRT1D(:)>XRTMIN(3)) .AND. (PRIS1D(:)>XRTMIN(4)/PTSTEP) .AND. (PRRS1D(:)>XRTMIN(3)/PTSTEP) )
-  ZZW1(:,3) = MIN( PRIS1D(:),XICFRR * PRIT1D(:) * PCRT1D(:)          & ! RICFRRG
-                                  * PLBDAR(:)**XEXICFRR        &
-                                  * PRHODREF(:)**(-XCEXVT-1.0) )
-!
-  ZZW1(:,4) = MIN( PRRS1D(:),XRCFRI * PCIT1D(:) * PCRT1D(:)          & ! RRCFRIG
-                                  * PLBDAR(:)**XEXRCFRI        &
-                                  * PRHODREF(:)**(-XCEXVT-2.0) )
-  PRIS1D(:) = PRIS1D(:) - ZZW1(:,3)
-  PRRS1D(:) = PRRS1D(:) - ZZW1(:,4)
-  PRGS1D(:) = PRGS1D(:) + ZZW1(:,3) + ZZW1(:,4)
-  PTHS1D(:) = PTHS1D(:) + ZZW1(:,4) * (PLSFACT(:) - PLVFACT(:)) ! f(L_f*RRCFRIG)
-!
-  PCIS1D(:) = MAX( PCIS1D(:)-ZZW1(:,3)*(PCIT1D(:)/PRIT1D(:)),0.0 )     ! CICFRRG
-  PCRS1D(:) = MAX( PCRS1D(:)-ZZW1(:,4)*(PCRT1D(:)/PRRT1D(:)),0.0 )     ! CRCFRIG
+   ZZW1(:,3) = MIN( PRIS1D(:),XICFRR * PRIT1D(:) * PCRT1D(:)          & ! RICFRRG
+                                   * PLBDAR(:)**XEXICFRR        &
+                                   * PRHODREF(:)**(-XCEXVT-1.0) )
+!
+   ZZW1(:,4) = MIN( PRRS1D(:),XRCFRI * PCIT1D(:) * PCRT1D(:)          & ! RRCFRIG
+                                   * PLBDAR(:)**XEXRCFRI        &
+                                   * PRHODREF(:)**(-XCEXVT-2.0) )
+   PRIS1D(:) = PRIS1D(:) - ZZW1(:,3)
+   PRRS1D(:) = PRRS1D(:) - ZZW1(:,4)
+   PRGS1D(:) = PRGS1D(:) + ZZW1(:,3)+ZZW1(:,4)
+   PTHS1D(:) = PTHS1D(:) + ZZW1(:,4)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*RRCFRIG)
+!
+   PCIS1D(:) = MAX( PCIS1D(:)-ZZW1(:,3)*(PCIT1D(:)/PRIT1D(:)),0.0 )     ! CICFRRG
+   PCRS1D(:) = MAX( PCRS1D(:)-ZZW1(:,4)*(PCRT1D(:)/PRRT1D(:)),0.0 )     ! CRCFRIG
+   PCGS1D(:) = PCGS1D(:)+ZZW1(:,3)*(PCIT1D(:)/PRIT1D(:)) 
 END WHERE
 !
 if ( nbumod == kmi .and. lbu_enable ) then
@@ -1003,6 +1121,8 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                          Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
   if ( lbudget_sv ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'CFRZ', &
                                          Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+  if ( lbudget_sv .and. NMOM_G.GE.2 ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'CFRZ', &
+                                         Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
 end if
 !
 !
@@ -1115,22 +1235,34 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                           Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
                   call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'WETG', &
                                           Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+       if(NMOM_S.GE.2 .AND. NMOM_G.GE.2 .AND. NMOM_H.GE.2) then
+                  call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'WETG', &
+                                         Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+                  call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'WETG', &
+                                         Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
+                  call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'WETG', &
+                                         Unpack( pchs1d(:), mask = gmicro(:, :, :), field = pchs(:, :, :) ) * prhodj(:, :, :) )                                   
+       end if
   end if
 end if
 !
 ZZW1(:,:) = 0.0
-WHERE( ((PRCT1D(:)>XRTMIN(2)) .AND. (PRGT1D(:)>XRTMIN(6)) .AND. (PRCS1D(:)>XRTMIN(2)/PTSTEP)) .OR. &
-       ((PRIT1D(:)>XRTMIN(4)) .AND. (PRGT1D(:)>XRTMIN(6)) .AND. (PRIS1D(:)>XRTMIN(4)/PTSTEP))      )
-  ZZW(:) = PLBDAG(:)**(XCXG-XDG-2.0) * PRHODREF(:)**(-XCEXVT)
-  ZZW1(:,1) = MIN( PRCS1D(:),XFCDRYG * PRCT1D(:) * ZZW(:) )             ! RCDRYG
-  ZZW1(:,2) = MIN( PRIS1D(:),XFIDRYG * EXP( XCOLEXIG*(PZT(:)-XTT) ) &
-                                   * PRIT1D(:) * ZZW(:) )             ! RIDRYG
+ZZNW1(:,:) = 0.0   
+WHERE( ((PRCT1D(:)>XRTMIN(2)) .AND. (PRGT1D(:)>XRTMIN(6)) .AND. (PRCS1D(:)>XRTMIN(2)/PTSTEP))  .OR. &
+          ((PRIT1D(:)>XRTMIN(4)) .AND. (PRGT1D(:)>XRTMIN(6)) .AND. (PRIS1D(:)>XRTMIN(4)/PTSTEP))      )
+   ZZW(:) = PLBDAG(:)**(-XDG-2.0) * PRHODREF(:)**(-XCEXVT) * PCGT1D(:) 
+   ZZW1(:,1) = MIN( PRCS1D(:),XFCDRYG * PRCT1D(:) * ZZW(:) )             ! RCDRYG
+   ZZNW1(:,1) =  MIN( PCCS1D(:),ZZW1(:,1) * PCCT1D(:) / MAX(PRCT1D(:),XRTMIN(2)) )
+   ZZW1(:,2) = MIN( PRIS1D(:),XFIDRYG * EXP( XCOLEXIG*(PZT(:)-XTT) ) &
+                                    * PRIT1D(:) * ZZW(:) )             ! RIDRYG
+   ZZNW1(:,2) =  MIN( PCIS1D(:),ZZW1(:,2) * PCIT1D(:) / MAX(PRIT1D(:),XRTMIN(4)) )
 END WHERE
 !
 !*       2.3.1  accretion of aggregates on the graupeln
 !        ----------------------------------------------
 !
 GDRY(:) = (PRST1D(:)>XRTMIN(5)) .AND. (PRGT1D(:)>XRTMIN(6)) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP)
+if (M2_ICE) GDRY(:) = GDRY(:) .AND. (PCST1D(:)>XCTMIN(5)) .AND. PCSS1D(:)>XCTMIN(5)/PTSTEP
 IGDRY = COUNT( GDRY(:) )
 !
 IF( IGDRY>0 ) THEN
@@ -1176,13 +1308,38 @@ IF( IGDRY>0 ) THEN
   ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
 !
   WHERE( GDRY(:) )
-    ZZW1(:,3) = MIN( PRSS1D(:),XFSDRYG*ZZW(:)                         & ! RSDRYG
+     ZZW1(:,3) = MIN( PRSS1D(:),XFSDRYG*ZZW(:)                         & ! RSDRYG
                                       * EXP( XCOLEXSG*(PZT(:)-XTT) )  &
-                    *( PRST1D(:)) )*( PLBDAG(:)**XCXG )    &
-                    *( PRHODREF(:)**(-XCEXVT) )                    &
+                    *( PRST1D(:) )*( PCGT1D(:) )      &  
+                    *( PRHODREF(:)**(-XCEXVT+1.) )                    &
                          *( XLBSDRYG1/( PLBDAG(:)**2              ) + &
                             XLBSDRYG2/( PLBDAG(:)   * PLBDAS(:)   ) + &
-                            XLBSDRYG3/(               PLBDAS(:)**2) )
+                            XLBSDRYG3/(               PLBDAS(:)**2) ) )
+  END WHERE
+!
+!                                                                           
+!*       2.2.5b  perform the bilinear interpolation of the normalized
+!               SDRYG-kernel FOR CONCENTRATION
+!
+  DO JJ = 1,IGDRY
+     ZVEC3(JJ) =  (  XKER_N_SDRYG(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                        - XKER_N_SDRYG(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+      	    		 	                                  * ZVEC1(JJ) &
+                     - (  XKER_N_SDRYG(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                        - XKER_N_SDRYG(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                        * (ZVEC1(JJ) - 1.0)
+  END DO
+  ZZNW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
+!
+  WHERE( GDRY(:)   )
+     ZZNW1(:,3) = MIN( PCSS1D(:),XFNSDRYG*ZZNW(:)                       & ! NSDRYG
+                                   * EXP( XCOLEXSG*(PZT(:)-XTT) )  &
+                 *( PCST1D(:)                     )*( PCGT1D(:) )      &  
+                 *( PRHODREF(:)**(-XCEXVT+1.) )                    &
+                      *( XLBNSDRYG1/( PLBDAG(:)**2             ) + &
+                         XLBNSDRYG2/( PLBDAG(:)   * PLBDAS(:)  ) + &
+                         XLBNSDRYG3/(               PLBDAS(:)**2)) )
+!                                                                        
   END WHERE
   DEALLOCATE(IVEC2)
   DEALLOCATE(IVEC1)
@@ -1240,12 +1397,33 @@ IF( IGDRY>0 ) THEN
   ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
 !
   WHERE( GDRY(:) )
-    ZZW1(:,4) = MIN( PRRS1D(:),XFRDRYG*ZZW(:) * PCRT1D(:)          & ! RRDRYG
-                      *( PLBDAR(:)**(-3) )*( PLBDAG(:)**XCXG ) &
-                              *( PRHODREF(:)**(-XCEXVT-1.) )   &
-                  *( XLBRDRYG1/( PLBDAG(:)**2              ) + &
-                     XLBRDRYG2/( PLBDAG(:)   * PLBDAR(:)   ) + &
-                     XLBRDRYG3/(               PLBDAR(:)**2) ) )
+     ZZW1(:,4) = MIN( PRRS1D(:),XFRDRYG*ZZW(:) * PRRT1D(:) * PCGT1D(:) & !
+                                *( PRHODREF(:)**(-XCEXVT+1.) )   &
+                    *( XLBRDRYG1/( PLBDAG(:)**2              ) + &
+                       XLBRDRYG2/( PLBDAG(:)   * PLBDAR(:)   ) + &
+                       XLBRDRYG3/(               PLBDAR(:)**2) ) )    
+  END WHERE
+!
+!                                                                          
+!*       2.2.10b perform the bilinear interpolation of the normalized
+!               RDRYG-kernel FOR CONCENTRATION
+!
+  DO JJ = 1,IGDRY
+     ZVEC3(JJ) =  (  XKER_N_RDRYG(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                        - XKER_N_RDRYG(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                         			 	                  * ZVEC1(JJ) &
+                    - (  XKER_N_RDRYG(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                        - XKER_N_RDRYG(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                                     			     * (ZVEC1(JJ) - 1.0)
+  END DO
+  ZZNW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
+!
+  WHERE( GDRY(:)  )
+     ZZNW1(:,4) = MIN( PCRS1D(:),XFNRDRYG*ZZNW(:) * PCRT1D(:) * PCGT1D(:) & ! NRDRYG
+                                *( PRHODREF(:)**(-XCEXVT+1.) )      &
+                    *( XLBNRDRYG1/( PLBDAG(:)**2              ) +   &
+                       XLBNRDRYG2/( PLBDAG(:)   * PLBDAR(:)   ) +   &
+                       XLBNRDRYG3/(               PLBDAR(:)**2) ) )
   END WHERE
   DEALLOCATE(IVEC2)
   DEALLOCATE(IVEC1)
@@ -1263,24 +1441,30 @@ ZRDRYG(:) = ZZW1(:,1) + ZZW1(:,2) + ZZW1(:,3) + ZZW1(:,4)
 ZZW(:) = 0.0
 ZRWETG(:) = 0.0
 WHERE( PRGT1D(:)>XRTMIN(6) )
-  ZZW1(:,5) = MIN( PRIS1D(:),                                    &
-              ZZW1(:,2) / (XCOLIG*EXP(XCOLEXIG*(PZT(:)-XTT)) ) ) ! RIWETG
-  ZZW1(:,6) = MIN( PRSS1D(:),                                    &
-              ZZW1(:,3) / (XCOLSG*EXP(XCOLEXSG*(PZT(:)-XTT)) ) ) ! RSWETG
-!
-  ZZW(:) = PRVT1D(:)*PPRES(:)/((XMV/XMD)+PRVT1D(:)) ! Vapor pressure
-  ZZW(:) =  PKA(:)*(XTT-PZT(:)) +                                  &
-            ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT ))   &
-                          *(XESTT-ZZW(:))/(XRV*PZT(:))             )
-!
+   ZZW1(:,5) = MIN( PRIS1D(:),                                    &
+               ZZW1(:,2) / (XCOLIG*EXP(XCOLEXIG*(PZT(:)-XTT)) ) ) ! RIWETG
+   ZZNW1(:,5) = MIN( PCIS1D(:),                      &  
+               ZZNW1(:,2) / (XCOLIG*EXP(XCOLEXIG*(PZT(:)-XTT)) ) )  ! NIWETG
+   ZZW1(:,6) = MIN( PRSS1D(:),                                    &
+               ZZW1(:,3) / (XCOLSG*EXP(XCOLEXSG*(PZT(:)-XTT)) ) ) ! RSWETG
+!
+!
+   ZZNW1(:,6) = MIN( PCSS1D(:),                      &       
+               ZZNW1(:,3) / (XCOLSG*EXP(XCOLEXSG*(PZT(:)-XTT)) ) )  ! NSWETG
+!
+   ZZW(:) = PRVT1D(:)*PPRES(:)/((XMV/XMD)+PRVT1D(:)) ! Vapor pressure
+   ZZW(:) =  PKA(:)*(XTT-PZT(:)) +                                  &
+                ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT ))   &
+                           *(XESTT-ZZW(:))/(XRV*PZT(:))             )
+!  
 ! compute RWETG
 !
-  ZRWETG(:)  = MAX( 0.0,                                               &
-                  ( ZZW(:) * ( X0DEPG*       PLBDAG(:)**XEX0DEPG +     &
-                               X1DEPG*PCJ(:)*PLBDAG(:)**XEX1DEPG ) +   &
-                  ( ZZW1(:,5)+ZZW1(:,6) ) *                            &
-                  ( PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PZT(:)))   ) ) / &
-                                  ( PRHODREF(:)*(XLMTT-XCL*(XTT-PZT(:))) )   )
+   ZRWETG(:)  = MAX( 0.0,                                               &
+                   ( ZZW(:) * PCGT1D(:) * ( X0DEPG* PLBDAG(:)**XEX0DEPG + &
+                    X1DEPG*PCJ(:)*PLBDAG(:)**XEX1DEPG ) +   &
+                   ( 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
@@ -1294,37 +1478,44 @@ END WHERE
 ZZW(:) = 0.0
 NHAIL = 0.
 IF (LHAIL) NHAIL = 1. 
-WHERE( PRGT1D(:)>XRTMIN(6) .AND. PZT(:)<XTT .AND. &
-       (ZRDRYG(:)-ZZW1(:,2)-ZZW1(:,3))>=(ZRWETG(:)-ZZW1(:,5)-ZZW1(:,6)) .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
+   ZZW(:) = ZRWETG(:) - ZZW1(:,5) - ZZW1(:,6) ! RCWETG+RRWETG
 !   
 ! limitation of the available rainwater mixing ratio (RRWETH < RRS !)
 !   
-  ZZW1(:,7) = MAX( 0.0,MIN( ZZW(:),PRRS1D(:)+ZZW1(:,1) ) )
-  ZZX(:)    = ZZW1(:,7) / ZZW(:)
-  ZZW1(:,5) = ZZW1(:,5)*ZZX(:)
-  ZZW1(:,6) = ZZW1(:,6)*ZZX(:)
-  ZRWETG(:) = ZZW1(:,7) + ZZW1(:,5) + ZZW1(:,6)
-!   
-  PRCS1D(:) = PRCS1D(:) - ZZW1(:,1)
-  PRIS1D(:) = PRIS1D(:) - ZZW1(:,5)
-  PRSS1D(:) = PRSS1D(:) - ZZW1(:,6)
+   ZZW1(:,7) = MAX( 0.0,MIN( ZZW(:),PRRS1D(:)+ZZW1(:,1) ) )
+   ZZX(:)    = ZZW1(:,7) / ZZW(:)
+   ZZW1(:,5) = ZZW1(:,5)*ZZX(:)
+   ZZNW1(:,5) = ZZNW1(:,5)*ZZX(:) 
+   ZZW1(:,6) = ZZW1(:,6)*ZZX(:)
+   ZZNW1(:,6) = ZZNW1(:,6)*ZZX(:) 
+   ZRWETG(:) = ZZW1(:,7) + ZZW1(:,5) + ZZW1(:,6)
+!     
+   PRCS1D(:) = PRCS1D(:) - ZZW1(:,1)
+   PRIS1D(:) = PRIS1D(:) - ZZW1(:,5)
+   PRSS1D(:) = PRSS1D(:) - ZZW1(:,6)
 !
 ! assume a linear percent of conversion of graupel into hail
 !
-  PRGS1D(:) = PRGS1D(:) + ZRWETG(:)
-  ZZW(:)  = PRGS1D(:)*ZRDRYG(:)*NHAIL/(ZRWETG(:)+ZRDRYG(:)) 
-  PRGS1D(:) = PRGS1D(:) - ZZW(:)                        
-  PRHS1D(:) = PRHS1D(:) + ZZW(:)
-  PRRS1D(:) = MAX( 0.0,PRRS1D(:) - ZZW1(:,7) + ZZW1(:,1) )
-  PTHS1D(:) = PTHS1D(:) + ZZW1(:,7) * (PLSFACT(:) - PLVFACT(:))
-                                                ! f(L_f*(RCWETG+RRWETG))
-!
-  PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,1)*(PCCT1D(:)/MAX(PRCT1D(:),XRTMIN(2))),0.0 )
-  PCIS1D(:) = MAX( PCIS1D(:)-ZZW1(:,5)*(PCIT1D(:)/MAX(PRIT1D(:),XRTMIN(4))),0.0 )
-  PCRS1D(:) = MAX( PCRS1D(:)-MAX( ZZW1(:,7)-ZZW1(:,1),0.0 )                 &
-                        *(PCRT1D(:)/MAX(PRRT1D(:),XRTMIN(3))),0.0 )
+   PRGS1D(:) = PRGS1D(:) + ZRWETG(:)
+   ZZW(:)  = PRGS1D(:)*ZRDRYG(:)*NHAIL/(ZRWETG(:)+ZRDRYG(:)) 
+   PRGS1D(:) = PRGS1D(:) - ZZW(:)                        
+   PRHS1D(:) = PRHS1D(:) + ZZW(:)
+   PRRS1D(:) = MAX( 0.0,PRRS1D(:) - ZZW1(:,7) + ZZW1(:,1) )
+   PTHS1D(:) = PTHS1D(:) + ZZW1(:,7)*(PLSFACT(:)-PLVFACT(:))
+                                                 ! f(L_f*(RCWETG+RRWETG))
+!
+   PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,1)*(PCCT1D(:)/MAX(PRCT1D(:),XRTMIN(2))),0.0 )
+   PCIS1D(:) = MAX( PCIS1D(:)-ZZNW1(:,5),0.0 ) 
+   PCRS1D(:) = MAX( PCRS1D(:)-MAX( ZZW1(:,7)-ZZW1(:,1),0.0 )                 &
+				   *(PCRT1D(:)/MAX(PRRT1D(:),XRTMIN(3))),0.0 )
+   PCSS1D(:) = MAX( PCSS1D(:)-ZZNW1(:,6),0.0 ) 
+   ZZNW(:)  = PCGS1D(:)*ZRDRYG(:)*NHAIL/(ZRWETG(:)+ZRDRYG(:)) 
+   PCGS1D(:) = MAX( PCGS1D(:)-ZZNW(:),0.0 ) 
+   PCHS1D(:) = MAX( PCHS1D(:)+ZZNW(:),0.0 ) 
 END WHERE
 !
 ! Budget storage
@@ -1350,6 +1541,14 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                          Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
                   call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'WETG', &
                                          Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+       if(M2_ICE) then
+                  call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'WETG', &
+                                         Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+                  call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'WETG', &
+                                         Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
+                  call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'WETG', &
+                                         Unpack( pchs1d(:), mask = gmicro(:, :, :), field = pchs(:, :, :) ) * prhodj(:, :, :) )
+       end if
   end if
 end if
 !
@@ -1375,23 +1574,25 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                           Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
                   call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'DRYG', &
                                           Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+    if (M2_ICE)   call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'DRYG', &
+                                          Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
   end if
 end if
 !
-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)
-  PRRS1D(:) = PRRS1D(:) - ZZW1(:,4)
-  PRGS1D(:) = PRGS1D(:) + ZRDRYG(:)
-  PTHS1D(:) = PTHS1D(:) + (ZZW1(:,1)+ZZW1(:,4)) * (PLSFACT(:) - PLVFACT(:)) !
-                                                        ! f(L_f*(RCDRYG+RRDRYG))
-!
-  PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,1)*(PCCT1D(:)/MAX(PRCT1D(:),XRTMIN(2))),0.0 )
-  PCIS1D(:) = MAX( PCIS1D(:)-ZZW1(:,2)*(PCIT1D(:)/MAX(PRIT1D(:),XRTMIN(4))),0.0 )
-  PCRS1D(:) = MAX( PCRS1D(:)-ZZW1(:,4)*(PCRT1D(:)/MAX(PRRT1D(:),XRTMIN(3))),0.0 ) 
-                                                        ! Approximate rates
+WHERE( PRGT1D(:)>XRTMIN(6) .AND. PZT(:)<XTT                              &
+     .AND. (ZRDRYG(:)-ZZW1(:,2)-ZZW1(:,3))<(ZRWETG(:)-ZZW1(:,5)-ZZW1(:,6)) .AND. ZRDRYG(:)>0.0 ) ! case
+   PRCS1D(:) = PRCS1D(:) - ZZW1(:,1)
+   PRIS1D(:) = PRIS1D(:) - ZZW1(:,2)
+   PRSS1D(:) = PRSS1D(:) - ZZW1(:,3)
+   PRRS1D(:) = PRRS1D(:) - ZZW1(:,4)
+   PRGS1D(:) = PRGS1D(:) + ZRDRYG(:)
+   PTHS1D(:) = PTHS1D(:) + (ZZW1(:,1)+ZZW1(:,4))*(PLSFACT(:)-PLVFACT(:)) !
+  						        ! f(L_f*(RCDRYG+RRDRYG))
+!
+   PCCS1D(:) = MAX( PCCS1D(:)-ZZNW1(:,1),0.0 )                                    
+   PCIS1D(:) = MAX( PCIS1D(:)-ZZNW1(:,2),0.0 )                                     
+   PCRS1D(:) = MAX( PCRS1D(:)-ZZNW1(:,4),0.0 )                                    
+   PCSS1D(:) = MAX( PCSS1D(:)-ZZNW1(:,3),0.0 )
 END WHERE
 !
 ! Budget storage
@@ -1415,6 +1616,8 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                          Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
                   call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'DRYG', &
                                          Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+    if (M2_ICE) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'DRYG', &
+                                         Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
   end if
 end if
 !
@@ -1485,31 +1688,33 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                           Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
   if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'GMLT', &
                                           Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
+  if ( lbudget_sv .and. M2_ICE) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'GMLT', &
+                                          Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )              
 end if
 
 ZZW(:) = 0.0
-WHERE( (PRGT1D(:)>XRTMIN(6)) .AND. (PRGS1D(:)>XRTMIN(6)/PTSTEP) .AND. (PZT(:)>XTT) )
-  ZZW(:) = PRVT1D(:)*PPRES(:)/((XMV/XMD)+PRVT1D(:)) ! Vapor pressure
-  ZZW(:) =  PKA(:)*(XTT-PZT(:)) +                                 &
-             ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
-                         *(XESTT-ZZW(:))/(XRV*PZT(:))             )
+   WHERE( (PRGT1D(:)>XRTMIN(6)) .AND. (PRGS1D(:)>XRTMIN(6)/PTSTEP) .AND. (PZT(:)>XTT) )
+      ZZW(:) = PRVT1D(:)*PPRES(:)/((XMV/XMD)+PRVT1D(:)) ! Vapor pressure
+      ZZW(:) =  PKA(:)*(XTT-PZT(:)) +                                 &
+              ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
+                          *(XESTT-ZZW(:))/(XRV*PZT(:))             )
 !
 ! compute RGMLTR
 !
-  ZZW(:)  = MIN( PRGS1D(:), MAX( 0.0,( -ZZW(:) *                     &
-                         ( X0DEPG*       PLBDAG(:)**XEX0DEPG +     &
-                           X1DEPG*PCJ(:)*PLBDAG(:)**XEX1DEPG ) -   &
-                                   ( ZZW1(:,1)+ZZW1(:,4) ) *       &
-                            ( PRHODREF(:)*XCL*(XTT-PZT(:))) ) /    &
-                                           ( PRHODREF(:)*XLMTT ) ) )
-  PRRS1D(:) = PRRS1D(:) + ZZW(:)
-  PRGS1D(:) = PRGS1D(:) - ZZW(:)
-  PTHS1D(:) = PTHS1D(:) - ZZW(:) * (PLSFACT(:) - PLVFACT(:)) ! f(L_f*(-RGMLTR))
-!
-!   PCRS1D(:) = MAX( PCRS1D(:) + ZZW(:)*(XCCG*PLBDAG(:)**XCXG/PRGT1D(:)),0.0 )
-  PCRS1D(:) = PCRS1D(:) + ZZW(:)*5.0E6  ! obtained after averaging
-                                    ! Dshed=1mm and 500 microns
-END WHERE
+      ZZW(:)  = MIN( PRGS1D(:), MAX( 0.0,( -ZZW(:) * PCGT1D(:) *           & 
+                          ( X0DEPG*       PLBDAG(:)**XEX0DEPG +     &
+                            X1DEPG*PCJ(:)*PLBDAG(:)**XEX1DEPG ) -   &
+                                    ( ZZW1(:,1)+ZZW1(:,4) ) *       &
+                             ( PRHODREF(:)*XCL*(XTT-PZT(:))) ) /    &
+                                            ( PRHODREF(:)*XLMTT ) ) )
+      PRRS1D(:) = PRRS1D(:) + ZZW(:)
+      PRGS1D(:) = PRGS1D(:) - ZZW(:)
+      PTHS1D(:) = PTHS1D(:) - ZZW(:)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(-RGMLTR))
+!
+      PCRS1D(:) = PCRS1D(:) + ZZW(:)*5.0E6  ! obtained after averaging
+                                        ! Dshed=1mm and 500 microns
+      PCGS1D(:) = MAX( PCGS1D(:) - ZZW(:)*5.0E6, 0.0)
+   END WHERE    
 !
 if ( nbumod == kmi .and. lbu_enable ) then
   if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'GMLT', &
@@ -1520,6 +1725,8 @@ if ( nbumod == kmi .and. lbu_enable ) then
                                          Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
   if ( lbudget_sv ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'GMLT', &
                                          Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
+  if ( lbudget_sv .and. M2_ICE) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'GMLT', &
+                                          Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )
 end if
 !
 !
@@ -1532,7 +1739,7 @@ end if
 !
 HAIL: IF (LHAIL) THEN
 !
-GHAIL(:) = PRTH1D(:)>XRTMIN(7)
+GHAIL(:) = PRHT1D(:)>XRTMIN(7)
 IHAIL = COUNT(GHAIL(:))
 !
 IF( IHAIL>0 ) THEN
@@ -1562,15 +1769,24 @@ IF( IHAIL>0 ) THEN
                                             Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
                     call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'WETH', &
                                             Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+      if (M2_ICE) then
+                    call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'WETH', &
+                                           Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+                    call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'WETH', &
+                                           Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) )                                     
+      end if
     end if
   end if
 
   ZZW1(:,:) = 0.0
-  WHERE( GHAIL(:) .AND. ( (PRCT1D(:)>XRTMIN(2) .AND. PRCS1D(:)>XRTMIN(2)/PTSTEP) .OR. &
-                          (PRIT1D(:)>XRTMIN(4) .AND. PRIS1D(:)>XRTMIN(4)/PTSTEP) )    )    
-    ZZW(:) = PLBDAH(:)**(XCXH-XDH-2.0) * PRHODREF(:)**(-XCEXVT)
-    ZZW1(:,1) = MIN( PRCS1D(:),XFWETH * PRCT1D(:) * ZZW(:) )             ! RCWETH
-    ZZW1(:,2) = MIN( PRIS1D(:),XFWETH * PRIT1D(:) * ZZW(:) )             ! RIWETH
+  ZZNW1(:,:) = 0.0
+  WHERE( GHAIL(:) .AND. ( (PRCT1D(:)>XRTMIN(2) .AND. PRCS1D(:)>XRTMIN(2)/PTSTEP ) .OR. &
+                              (PRIT1D(:)>XRTMIN(4) .AND. PRIS1D(:)>XRTMIN(4)/PTSTEP)  )    )   
+     ZZW(:) = PCHT1D(:) * PLBDAH(:)**(-XDH-2.0) * PRHODREF(:)**(1-XCEXVT)
+     ZZW1(:,1) = MIN( PRCS1D(:),XFWETH * PRCT1D(:) * ZZW(:) )             ! RCWETH
+     ZZNW1(:,1) = MIN( PCCS1D(:),XFWETH * PCCT1D(:) * ZZW(:) ) !ZZW1(:,1) * ZCCT(:) / MAX(ZRCT(:),XRTMIN(2)) )            ! NCWETH 
+     ZZW1(:,2) = MIN( PRIS1D(:),XFWETH * PRIT1D(:) * ZZW(:) )             ! RIWETH
+     ZZNW1(:,2) = MIN( PCIS1D(:),XFWETH * PCIT1D(:) * ZZW(:) ) !ZZW1(:,2) * ZCIT(:) / MAX(ZRIT(:),XRTMIN(4)) )      ! NIWETH 
   END WHERE
 !
 !*       3.1.1  accretion of aggregates on the hailstones
@@ -1622,12 +1838,33 @@ IF( IHAIL>0 ) THEN
     ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 )
 !
     WHERE( GWET(:) )
-      ZZW1(:,3) = MIN( PRSS1D(:),XFSWETH*ZZW(:)                       & ! RSWETH
-                      *( PRST1D(:))*( PLBDAH(:)**XCXH )  &
-                        *( PRHODREF(:)**(-XCEXVT) )               &
-                         *( XLBSWETH1/( PLBDAH(:)**2              ) + &
-                            XLBSWETH2/( PLBDAH(:)   * PLBDAS(:)   ) + &
-                            XLBSWETH3/(               PLBDAS(:)**2) ) )
+       ZZW1(:,3) = MIN( PRSS1D(:),XFSWETH*ZZW(:) * PRST1D(:) * PCHT1D(:)  & 
+            *( PRHODREF(:)**(-XCEXVT+1.) )               &
+            *( XLBSWETH1/( PLBDAH(:)**2              ) + &
+            XLBSWETH2/( PLBDAH(:)   * PLBDAS(:)   ) + &
+            XLBSWETH3/(               PLBDAS(:)**2) ) )   
+    END WHERE
+!
+!*       3.1.5b  perform the bilinear interpolation of the normalized
+!               SWETH-kernel FOR CONCENTRATION
+!
+    DO JJ = 1,IGWET
+       ZVEC3(JJ) = (  XKER_N_SWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                     - XKER_N_SWETH(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+        		 	                                   * ZVEC1(JJ) &
+                   - ( XKER_N_SWETH(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                     - XKER_N_SWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+         		                                     * (ZVEC1(JJ) - 1.0)
+    END DO
+    ZZNW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 )
+!
+    WHERE( GWET(:) )
+       ZZNW1(:,3) = MIN( PCSS1D(:),XFNSWETH*ZZNW(:) * PCST1D(:)            & ! NSWETH
+                                              *( PCHT1D(:) )              & 
+       	                 *( PRHODREF(:)**(-XCEXVT+1.) )               &
+                         *( XLBNSWETH1/( PLBDAH(:)**2              ) + &
+                            XLBNSWETH2/( PLBDAH(:)   * PLBDAS(:)   ) + &
+                            XLBNSWETH3/(               PLBDAS(:)**2) ) )
     END WHERE
     DEALLOCATE(IVEC2)
     DEALLOCATE(IVEC1)
@@ -1685,13 +1922,34 @@ IF( IHAIL>0 ) THEN
     ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 )
 !
     WHERE( GWET(:) )
-      ZZW1(:,5) = MAX(MIN( PRGS1D(:),XFGWETH*ZZW(:)                       & ! RGWETH
-                    *( PLBDAG(:)**(XCXG-XBG) )*( PLBDAH(:)**XCXH )  &
-                       *( PRHODREF(:)**(-XCEXVT-1.) )               &
-                       *( XLBGWETH1/( PLBDAH(:)**2              ) + &
-                          XLBGWETH2/( PLBDAH(:)   * PLBDAG(:)   ) + &
-                          XLBGWETH3/(               PLBDAG(:)**2) ) ),0. )
+       ZZW1(:,5) = MAX(MIN( PRGS1D(:),XFGWETH*ZZW(:) * PRGT1D(:) * PCHT1D(:) &  
+            *( PRHODREF(:)**(-XCEXVT+1.) )               &
+            *( XLBGWETH1/( PLBDAH(:)**2              ) + &
+            XLBGWETH2/( PLBDAH(:)   * PLBDAG(:)   ) + &
+            XLBGWETH3/(               PLBDAG(:)**2) ) ),0. )
+    END WHERE
+!*       3.1.10b perform the bilinear interpolation of the normalized
+!               GWETH-kernel FOR CONCENTRATION
+!
+    DO JJ = 1,IGWET
+       ZVEC3(JJ) = (  XKER_N_GWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                        - XKER_N_GWETH(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                          			 	                   * ZVEC1(JJ) &
+                     - (  XKER_N_GWETH(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
+                        - XKER_N_GWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
+                                    			     * (ZVEC1(JJ) - 1.0)
+    END DO
+    ZZNW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 )
+!
+    WHERE( GWET(:) )
+       ZZNW1(:,5) = MAX(MIN( PCGS1D(:),XFNGWETH*ZZNW(:) * PCGT1D(:)       & ! NGWETH
+                                                   *( PCHT1D(:) )        &
+                         *( PRHODREF(:)**(-XCEXVT+1.) )               &
+                         *( XLBNGWETH1/( PLBDAH(:)**2              ) + &
+                            XLBNGWETH2/( PLBDAH(:)   * PLBDAG(:)   ) + &
+                            XLBNGWETH3/(               PLBDAG(:)**2) ) ),0. )
     END WHERE
+! 
     DEALLOCATE(IVEC2)
     DEALLOCATE(IVEC1)
     DEALLOCATE(ZVEC3)
@@ -1705,31 +1963,34 @@ IF( IHAIL>0 ) THEN
 !
   ZZW(:) = 0.0
   WHERE( GHAIL(:) .AND. PZT(:)<XTT )
-    ZZW(:) = PRVT1D(:) * PPRES(:) / ((XMV / XMD) + PRVT1D(:)) ! Vapor pressure
-    ZZW(:) = PKA(:) * (XTT - PZT(:)) +                            &
-           ( PDV(:) * (XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
-                    * (XESTT-ZZW(:))/(XRV*PZT(:))             )
+     ZZW(:) = PRVT1D(:)*PPRES(:)/((XMV/XMD)+PRVT1D(:)) ! Vapor pressure
+     ZZW(:) = PKA(:)*(XTT-PZT(:)) +                                 &
+                    ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
+                                *(XESTT-ZZW(:))/(XRV*PZT(:))             )
 !
 ! compute RWETH
 !
-     ZZW(:)  =  MAX(0.,  ( ZZW(:) * ( X0DEPH*       PLBDAH(:)**XEX0DEPH +     &
-                               X1DEPH*PCJ(:)*PLBDAH(:)**XEX1DEPH ) +   &
-                  ( ZZW1(:,2)+ZZW1(:,3)+ZZW1(:,5) ) *                  &
-                  ( PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PZT(:)))   ) ) / &
-                        ( PRHODREF(:)*(XLMTT-XCL*(XTT-PZT(:))) ) )
+     ZZW(:)  =  MAX(0.,  ( ZZW(:) * PCHT1D(:) * ( X0DEPH*       PLBDAH(:)**XEX0DEPH +     &  
+                                    X1DEPH*PCJ(:)*PLBDAH(:)**XEX1DEPH ) +   &
+                       ( ZZW1(:,2)+ZZW1(:,3)+ZZW1(:,5) ) *                  &
+                       ( PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PZT(:)))   ) ) / &
+                             ( PRHODREF(:)*(XLMTT-XCL*(XTT-PZT(:))) ) )            
 !
-     ZZW1(:,6) = MAX( ZZW(:) - ZZW1(:,2) - ZZW1(:,3) - ZZW1(:,5),0.) ! RCWETH+RRWETH
-   END WHERE
+     ZZW1(:,6) = MAX( ZZW(:) - ZZW1(:,2) - ZZW1(:,3) - ZZW1(:,5),0.) !RCWETH+RRWETH
+  END WHERE
    !
-   WHERE ( GHAIL(:) .AND. PZT(:)<XTT  .AND. ZZW1(:,6)/=0.)
+  WHERE ( GHAIL(:) .AND. PZT(:)<XTT  .AND. ZZW1(:,6)/=0.)
 !
 ! limitation of the available rainwater mixing ratio (RRWETH < RRS !)
 !
      ZZW1(:,4) = MAX( 0.0,MIN( ZZW1(:,6),PRRS1D(:)+ZZW1(:,1) ) )
      ZZX(:)    = ZZW1(:,4) / ZZW1(:,6)
      ZZW1(:,2) = ZZW1(:,2) * ZZX(:)
+     ZZNW1(:,2) = ZZNW1(:,2)*ZZX(:) 
      ZZW1(:,3) = ZZW1(:,3) * ZZX(:)
+     ZZNW1(:,3) = ZZNW1(:,3)*ZZX(:)
      ZZW1(:,5) = ZZW1(:,5) * ZZX(:)
+     ZZNW1(:,5) = ZZNW1(:,5)*ZZX(:)  
      ZZW(:)    = ZZW1(:,4) + ZZW1(:,2) + ZZW1(:,3) + ZZW1(:,5)
 !
 !*       3.2.1  integrate the Wet growth of hail
@@ -1743,10 +2004,12 @@ IF( IHAIL>0 ) THEN
      PTHS1D(:) = PTHS1D(:) + ZZW1(:,4) * (PLSFACT(:) - PLVFACT(:)) 
                                                 ! f(L_f*(RCWETH+RRWETH))
 !
-     PCCS1D(:) = MAX( PCCS1D(:)-ZZW1(:,1)*(PCCT1D(:)/MAX(PRCT1D(:),XRTMIN(2))),0.0 )
-     PCIS1D(:) = MAX( PCIS1D(:)-ZZW1(:,2)*(PCIT1D(:)/MAX(PRIT1D(:),XRTMIN(4))),0.0 )
      PCRS1D(:) = MAX( PCRS1D(:)-MAX( ZZW1(:,4)-ZZW1(:,1),0.0 )                 &
-                                     *(PCRT1D(:)/MAX(PRRT1D(:),XRTMIN(3))),0.0 )
+                                       *(PCRT1D(:)/MAX(PRRT1D(:),XRTMIN(3))),0.0 )
+     PCSS1D(:) = MAX( PCSS1D(:)-ZZNW1(:,3),0.0 ) 
+     PCGS1D(:) = MAX( PCGS1D(:)-ZZNW1(:,5),0.0 )   
+     PCCS1D(:) = MAX( PCCS1D(:)-ZZNW1(:,1),0.0 )       
+     PCIS1D(:) = MAX( PCIS1D(:)-ZZNW1(:,2),0.0 )
   END WHERE
 !
   if ( nbumod == kmi .and. lbu_enable ) then
@@ -1771,6 +2034,12 @@ IF( IHAIL>0 ) THEN
                                            Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
                     call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'WETH', &
                                            Unpack( pcis1d(:), mask = gmicro(:, :, :), field = pcis(:, :, :) ) * prhodj(:, :, :) )
+      if (M2_ICE) then
+                    call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'WETH', &
+                                           Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+                    call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'WETH', &
+                                           Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) ) 
+      end if
     end if
   end if
 END IF ! IHAIL>0
@@ -1787,13 +2056,19 @@ IF ( IHAIL>0 ) THEN
                                             Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
     if ( lbudget_rh ) call Budget_store_init( tbudgets(NBUDGET_RH), 'COHG', &
                                             Unpack( prhs1d(:), mask = gmicro(:, :, :), field = prhs(:, :, :) ) * prhodj(:, :, :) )
+    if (lbudget_sv .and. M2_ICE) then
+                      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'COHG', &
+                                           Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) ) 
+                      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'COHG', &
+                                           Unpack( pchs1d(:), mask = gmicro(:, :, :), field = pchs(:, :, :) ) * prhodj(:, :, :) )  
+    end if
   end if
 !
   ZTHRH = 0.01E-3
   ZTHRC = 0.001E-3
   ZZW(:) = 0.0
 !
-  WHERE( PRTH1D(:)<ZTHRH .AND. PRCT1D(:)<ZTHRC .AND. PZT(:)<XTT )
+  WHERE( PRHT1D(:)<ZTHRH .AND. PRCT1D(:)<ZTHRC .AND. PZT(:)<XTT )
     ZZW(:) = MIN( 1.0,MAX( 0.0,1.0-(PRCT1D(:)/ZTHRC) ) )
 !
 ! assume a linear percent conversion rate of hail into graupel
@@ -1802,12 +2077,24 @@ IF ( IHAIL>0 ) THEN
     PRGS1D(:) = PRGS1D(:) + ZZW(:)                      !   partial conversion
     PRHS1D(:) = PRHS1D(:) - ZZW(:)                      ! of hail into graupel
   END WHERE
+    if (M2_ICE) then
+       WHERE( PRHT1D(:)<ZTHRH .AND. PRCT1D(:)<ZTHRC .AND. PZT(:)<XTT )
+          PCGS1D(:) = PCGS1D(:) + ZZW(:)* PCHS1D(:)/MAX(PRHS1D(:),XRTMIN(7)) 
+          PCHS1D(:) = MAX( PCHS1D(:) - ZZW(:)* PCHS1D(:)/MAX(PRHS1D(:),XRTMIN(7)), 0.0 ) 
+       END WHERE
+    end if
 !
   if ( nbumod == kmi .and. lbu_enable ) then
     if ( lbudget_rg ) call Budget_store_end( tbudgets(NBUDGET_RG), 'COHG', &
                                            Unpack( prgs1d(:), mask = gmicro(:, :, :), field = prgs(:, :, :) ) * prhodj(:, :, :) )
     if ( lbudget_rh ) call Budget_store_end( tbudgets(NBUDGET_RH), 'COHG', &
                                            Unpack( prhs1d(:), mask = gmicro(:, :, :), field = prhs(:, :, :) ) * prhodj(:, :, :) )
+    if (lbudget_sv .and. M2_ICE) then
+                      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ng), 'COHG', &
+                                           Unpack( pcgs1d(:), mask = gmicro(:, :, :), field = pcgs(:, :, :) ) * prhodj(:, :, :) ) 
+                      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'COHG', &
+                                           Unpack( pchs1d(:), mask = gmicro(:, :, :), field = pchs(:, :, :) ) * prhodj(:, :, :) )                                   
+    end if
   end if
 END IF
 !
@@ -1815,52 +2102,128 @@ END IF
 !*       3.4    Melting of the hailstones
 !
 IF ( IHAIL>0 ) THEN
-  if ( nbumod == kmi .and. lbu_enable ) then
-    if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'HMLT', &
-                                            Unpack( pths1d(:), mask = gmicro(:, :, :), field = pths(:, :, :) ) * prhodj(:, :, :) )
-    if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR), 'HMLT', &
-                                            Unpack( prrs1d(:), mask = gmicro(:, :, :), field = prrs(:, :, :) ) * prhodj(:, :, :) )
-    if ( lbudget_rh ) call Budget_store_init( tbudgets(NBUDGET_RH), 'HMLT', &
-                                            Unpack( prhs1d(:), mask = gmicro(:, :, :), field = prhs(:, :, :) ) * prhodj(:, :, :) )
-    if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'HMLT', &
-                                            Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
-  end if
-!
-  ZZW(:) = 0.0
-  WHERE( GHAIL(:) .AND. (PRHS1D(:)>XRTMIN(7)/PTSTEP) .AND. (PRTH1D(:)>XRTMIN(7)) .AND. (PZT(:)>XTT) )
-    ZZW(:) = PRVT1D(:) * PPRES(:) / ((XMV / XMD) + PRVT1D(:)) ! Vapor pressure
-    ZZW(:) = PKA(:) * (XTT - PZT(:)) +                              &
-           ( PDV(:) * (XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
-                    * (XESTT - ZZW(:)) / (XRV * PZT(:)) )
+   if ( nbumod == kmi .and. lbu_enable ) then
+      if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'HMLT', &
+           Unpack( pths1d(:), mask = gmicro(:, :, :), field = pths(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR), 'HMLT', &
+           Unpack( prrs1d(:), mask = gmicro(:, :, :), field = prrs(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_rh ) call Budget_store_init( tbudgets(NBUDGET_RH), 'HMLT', &
+           Unpack( prhs1d(:), mask = gmicro(:, :, :), field = prhs(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'HMLT', &
+           Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_sv .and. M2_ICE ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'HMLT', &
+           Unpack( pchs1d(:), mask = gmicro(:, :, :), field = pchs(:, :, :) ) * prhodj(:, :, :) )                                    
+   end if
+!
+   ZZW(:) = 0.0
+   WHERE( GHAIL(:) .AND. (PRHS1D(:)>XRTMIN(7)/PTSTEP) .AND. (PRHT1D(:)>XRTMIN(7)) .AND. (PZT(:)>XTT) )
+      ZZW(:) = PRVT1D(:)*PPRES(:)/((XMV/XMD)+PRVT1D(:)) ! Vapor pressure
+      ZZW(:) = PKA(:)*(XTT-PZT(:)) +                              &
+              ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
+              *(XESTT-ZZW(:))/(XRV*PZT(:))         )
 !
 ! compute RHMLTR
 !
-    ZZW(:)  = MIN( PRHS1D(:), MAX( 0.0,( -ZZW(:) *                     &
-                           ( X0DEPH*       PLBDAH(:)**XEX0DEPH +     &
-                             X1DEPH*PCJ(:)*PLBDAH(:)**XEX1DEPH ) -   &
-                    ZZW1(:,6)*( PRHODREF(:)*XCL*(XTT-PZT(:))) ) /    &
-                                             ( PRHODREF(:)*XLMTT ) ) )
-    PRRS1D(:) = PRRS1D(:) + ZZW(:)
-    PRHS1D(:) = PRHS1D(:) - ZZW(:)
-    PTHS1D(:) = PTHS1D(:) - ZZW(:) * (PLSFACT(:) - PLVFACT(:)) ! f(L_f*(-RHMLTR))
+      ZZW(:)  = MIN( PRHS1D(:), MAX( 0.0,( -ZZW(:) * PCHT1D(:) *           &
+              ( X0DEPH*       PLBDAH(:)**XEX0DEPH +     &
+              X1DEPH*PCJ(:)*PLBDAH(:)**XEX1DEPH ) -   &
+              ZZW1(:,6)*( PRHODREF(:)*XCL*(XTT-PZT(:))) ) /    &
+              ( PRHODREF(:)*XLMTT ) ) )  
+      PRRS1D(:) = PRRS1D(:) + ZZW(:)
+      PRHS1D(:) = PRHS1D(:) - ZZW(:)
+      PTHS1D(:) = PTHS1D(:) - ZZW(:)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(-RHMLTR))
 !
-    PCRS1D(:) = MAX( PCRS1D(:) + ZZW(:)*(XCCH*PLBDAH(:)**XCXH/PRTH1D(:)),0.0 )
-  END WHERE
+      PCRS1D(:) = MAX( PCRS1D(:) + ZZW(:)*(PCHT1D(:)/PRHT1D(:)),0.0 )
+      PCHS1D(:) = MAX( PCHS1D(:) - ZZW(:)*(PCHT1D(:)/PRHT1D(:)),0.0 )
 !
-  if ( nbumod == kmi .and. lbu_enable ) then
-    if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'HMLT', &
-                                           Unpack( pths1d(:), mask = gmicro(:, :, :), field = pths(:, :, :) ) * prhodj(:, :, :) )
-    if ( lbudget_rr ) call Budget_store_end( tbudgets(NBUDGET_RR), 'HMLT', &
-                                           Unpack( prrs1d(:), mask = gmicro(:, :, :), field = prrs(:, :, :) ) * prhodj(:, :, :) )
-    if ( lbudget_rh ) call Budget_store_end( tbudgets(NBUDGET_RH), 'HMLT', &
-                                           Unpack( prhs1d(:), mask = gmicro(:, :, :), field = prhs(:, :, :) ) * prhodj(:, :, :) )
-    if ( lbudget_sv ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'HMLT', &
-                                           Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
-  end if
+   END WHERE
+!
+   if ( nbumod == kmi .and. lbu_enable ) then
+      if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'HMLT', &
+           Unpack( pths1d(:), mask = gmicro(:, :, :), field = pths(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_rr ) call Budget_store_end( tbudgets(NBUDGET_RR), 'HMLT', &
+           Unpack( prrs1d(:), mask = gmicro(:, :, :), field = prrs(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_rh ) call Budget_store_end( tbudgets(NBUDGET_RH), 'HMLT', &
+           Unpack( prhs1d(:), mask = gmicro(:, :, :), field = prhs(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_sv ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'HMLT', &
+           Unpack( pcrs1d(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
+      if ( lbudget_sv .and. M2_ICE ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nh), 'HMLT', &
+           Unpack( pchs1d(:), mask = gmicro(:, :, :), field = pchs(:, :, :) ) * prhodj(:, :, :) )
+   end if
 END IF
 !
 END IF HAIL
 !
+!
+!
+if( M2_ICE) then
+   if ( nbumod == kmi .and. lbu_enable ) then
+      if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'SSC', &
+           Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+   end if
+!*       1.3NEW  Aggregates self-collection  
+!        -------------------------------
+!
+   ZZNW1(:,:) = 0.0
+!
+   GACC(:) = (PCST1D(:)>XCTMIN(5)) .AND. (PZT(:)<XTT)
+   IGACC = COUNT( GACC(:) )
+!
+   IF( IGACC>0 ) THEN
+!
+!        1.3N.0  allocations
+!
+      ALLOCATE(ZVEC1(IGACC))
+      ALLOCATE(IVEC1(IGACC))
+!
+!        1.3N.1  select the (ZLBDAS,ZLBDAS) couplet
+!
+      ZVEC1(:) = PACK( PLBDAS(:),MASK=GACC(:) )
+!
+!        1.3N.2  find the next lower indice for the ZLBDAS and for the ZLBDAS
+!               in the geometrical set of (Lbda_s,Lbda_s) couplet use to
+!               tabulate the SACCS-kernel
+!
+      ZVEC1(1:IGACC) = MAX( 1.0001, MIN( FLOAT(NSCLBDAS)-0.0001,           &
+           XSCINTP1S * LOG( ZVEC1(1:IGACC) ) + XSCINTP2S ) )
+      IVEC1(1:IGACC) = INT( ZVEC1(1:IGACC) )
+      ZVEC1(1:IGACC) = ZVEC1(1:IGACC) - FLOAT( IVEC1(1:IGACC) )
+!
+!        1.3N.3 perform the bilinear interpolation of the normalized
+!               SSCS-kernel
+!
+      ALLOCATE(ZVEC3(IGACC))
+      DO JJ = 1,IGACC
+         ZVEC3(JJ) =  (   XKER_N_SSCS(IVEC1(JJ)+1,IVEC1(JJ)+1)* ZVEC1(JJ)          &
+                    -  XKER_N_SSCS(IVEC1(JJ)+1,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
+                                                                         * ZVEC1(JJ) &
+                 - (   XKER_N_SSCS(IVEC1(JJ)  ,IVEC1(JJ)+1)* ZVEC1(JJ)          &
+                    -  XKER_N_SSCS(IVEC1(JJ)  ,IVEC1(JJ)  )*(ZVEC1(JJ) - 1.0) ) &
+                                                           * (ZVEC1(JJ) - 1.0)
+      END DO
+      ZZNW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC(:),FIELD=0.0 ) !! NSACCS
+      DEALLOCATE(ZVEC3)
+!
+      WHERE( GACC(:) )               
+!          ZZNW1(:,5) = XSSC * (PRST1D(:) * PRHODREF(:))**EXPRS * PCST1D(:)**EXPNS      
+         ZZNW1(:,5) = MIN( PCSS1D(:),XFNSSCS*ZZNW(:)                       & ! NSSCS
+                                      * EXP( XCOLEXSS*(PZT(:)-XTT) )  &
+                    *( PCST1D(:)                     )**( 2 )           &  
+                    *( PRHODREF(:)**(-XCEXVT-1.) )                    &
+                         *( XLBNSSCS1/( PLBDAS(:)**2             ) +  &
+                            XLBNSSCS2/( PLBDAS(:)**2 ) )            )
+!                                                                           
+         PCSS1D(:) = MAX( PCSS1D(:)-ZZNW1(:,5),0.0 )    
+      END WHERE
+      DEALLOCATE(IVEC1)
+      DEALLOCATE(ZVEC1)
+   END IF
+   if ( nbumod == kmi .and. lbu_enable ) then
+      if ( lbudget_sv ) call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ns), 'SSC', &
+           Unpack( pcss1d(:), mask = gmicro(:, :, :), field = pcss(:, :, :) ) * prhodj(:, :, :) )
+   end if
+end if
+!
 !------------------------------------------------------------------------------
 !
 END SUBROUTINE LIMA_MIXED_FAST_PROCESSES
diff --git a/src/MNH/lima_mixed_slow_processes.f90 b/src/MNH/lima_mixed_slow_processes.f90
index 097d575acb9ec376d16d740abf2b96253c71e847..1daf983b673fe3e601bdfe4ccd33a7a5f7f93304 100644
--- a/src/MNH/lima_mixed_slow_processes.f90
+++ b/src/MNH/lima_mixed_slow_processes.f90
@@ -119,6 +119,7 @@ END MODULE MODI_LIMA_MIXED_SLOW_PROCESSES
 !!      C. Barthe  * LACy *  jan. 2014   add budgets
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
 !  P. Wautelet 02/02/2021: budgets: add missing source terms for SV budgets in LIMA
+!  M. Taufour     07/2022: add concentration for snow, graupel, hail
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -203,27 +204,15 @@ INTEGER :: JMOD_IFN
 !
 IF (LSNOW) THEN
    ZZW(:) = 0.0
-   if (NMOM_G.GE.2) then
-      WHERE ( (ZRGT(:)>XRTMIN(6)) .AND. (ZRGS(:)>XRTMIN(6)/PTSTEP) ) 
-          ZZW(:) = ( ZSSI(:)/ZAI(:)/ZRHODREF(:) ) *  ZCGT(:) *               &
-                  ( X0DEPG*ZLBDAG(:)**XEX0DEPG + X1DEPG*ZCJ(:)*ZLBDAG(:)**XEX1DEPG )
-          ZZW(:) =         MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
-                         - MIN( ZRGS(:),ABS(ZZW(:)) )*(0.5-SIGN(0.5,ZZW(:)))
-          ZRGS(:) = ZRGS(:) + ZZW(:)
-          ZRVS(:) = ZRVS(:) - ZZW(:)
-          ZTHS(:) = ZTHS(:) + ZZW(:)*ZLSFACT(:)
-       END WHERE
-   else
-       WHERE ( (ZRGT(:)>XRTMIN(6)) .AND. (ZRGS(:)>XRTMIN(6)/PTSTEP) )
-          ZZW(:) = ( ZSSI(:)/ZAI(:)/ZRHODREF(:) ) * XCCG *                      &
-                   ( X0DEPG*ZLBDAG(:)**(XCXG+XEX0DEPG) + X1DEPG*ZCJ(:)*ZLBDAG(:)**(XCXG+XEX1DEPG) )
-          ZZW(:) =         MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
-                         - MIN( ZRGS(:),ABS(ZZW(:)) )*(0.5-SIGN(0.5,ZZW(:)))
-          ZRGS(:) = ZRGS(:) + ZZW(:)
-          ZRVS(:) = ZRVS(:) - ZZW(:)
-          ZTHS(:) = ZTHS(:) + ZZW(:)*ZLSFACT(:)
-       END WHERE
-   end if
+   WHERE ( (ZRGT(:)>XRTMIN(6)) .AND. (ZRGS(:)>XRTMIN(6)/PTSTEP) ) 
+      ZZW(:) = ( ZSSI(:)/ZAI(:)/ZRHODREF(:) ) *  ZCGT(:) *               &
+           ( X0DEPG*ZLBDAG(:)**XEX0DEPG + X1DEPG*ZCJ(:)*ZLBDAG(:)**XEX1DEPG )
+      ZZW(:) =         MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
+                     - MIN( ZRGS(:),ABS(ZZW(:)) )*(0.5-SIGN(0.5,ZZW(:)))
+      ZRGS(:) = ZRGS(:) + ZZW(:)
+      ZRVS(:) = ZRVS(:) - ZZW(:)
+      ZTHS(:) = ZTHS(:) + ZZW(:)*ZLSFACT(:)
+   END WHERE
 !
 ! Budget storage
   if ( nbumod == kmi .and. lbu_enable ) then