diff --git a/src/common/micro/ice4_nucleation.func.h b/src/common/micro/ice4_nucleation.func.h
new file mode 100644
index 0000000000000000000000000000000000000000..758c515754d4594686a3fc0d78c39eefdb0ddb7a
--- /dev/null
+++ b/src/common/micro/ice4_nucleation.func.h
@@ -0,0 +1,123 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+ELEMENTAL SUBROUTINE ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, ODCOMPUTE, &
+                           PTHT, PPABST, PRHODREF, PEXN, PLSFACT, PT, &
+                           PRVT, &
+                           PCIT, PRVHENI_MR)
+!!
+!!**  PURPOSE
+!!    -------
+!!      Computes the nucleation
+!!
+!!    AUTHOR
+!!    ------
+!!      S. Riette from the splitting of rain_ice source code (nov. 2014)
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!
+!!     R. El Khatib 24-Aug-2021 Optimizations
+!
+!
+!*      0. DECLARATIONS
+!          ------------
+!
+USE MODD_CST,            ONLY: CST_t
+USE MODD_PARAM_ICE_n,      ONLY: PARAM_ICE_t
+USE MODD_RAIN_ICE_DESCR_n, ONLY: RAIN_ICE_DESCR_t
+USE MODD_RAIN_ICE_PARAM_n, ONLY: RAIN_ICE_PARAM_t
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+TYPE(CST_t),              INTENT(IN)    :: CST
+TYPE(PARAM_ICE_t),        INTENT(IN)    :: PARAMI
+TYPE(RAIN_ICE_PARAM_t),   INTENT(IN)    :: ICEP
+TYPE(RAIN_ICE_DESCR_t),   INTENT(IN)    :: ICED
+LOGICAL, INTENT(IN)    :: ODCOMPUTE
+REAL,    INTENT(IN)    :: PTHT    ! Theta at t
+REAL,    INTENT(IN)    :: PPABST  ! absolute pressure at t
+REAL,    INTENT(IN)    :: PRHODREF! Reference density
+REAL,    INTENT(IN)    :: PEXN    ! Exner function
+REAL,    INTENT(IN)    :: PLSFACT
+REAL,    INTENT(IN)    :: PT      ! Temperature at time t
+REAL,    INTENT(IN)    :: PRVT    ! Water vapor m.r. at t
+REAL,    INTENT(INOUT) :: PCIT    ! Pristine ice n.c. at t
+REAL,    INTENT(OUT)   :: PRVHENI_MR ! Mixing ratio change due to the heterogeneous nucleation
+!
+!*       0.2  declaration of local variables
+!
+REAL  :: ZW ! work array
+LOGICAL  :: GNEGT  ! Test where to compute the HEN process
+REAL   :: ZZW,      & ! Work array
+                           ZUSW,     & ! Undersaturation over water
+                           ZSSI        ! Supersaturation over ice
+!-------------------------------------------------------------------------------
+!
+!
+  IF (ODCOMPUTE) THEN
+GNEGT=PT<CST%XTT .AND. PRVT>ICED%XRTMIN(1)
+ELSE
+GNEGT=.FALSE.
+END IF
+
+ZUSW=0.
+ZZW=0.
+
+IF (GNEGT) THEN
+ ZZW=ALOG(PT)
+ ZUSW=EXP(CST%XALPW - CST%XBETAW/PT - CST%XGAMW*ZZW)          ! es_w
+ ZZW=EXP(CST%XALPI - CST%XBETAI/PT - CST%XGAMI*ZZW)           ! es_i
+END IF
+  
+ZSSI=0.
+IF (GNEGT) THEN
+  ZZW=MIN(PPABST/2., ZZW)             ! safety limitation
+  ZSSI=PRVT*(PPABST-ZZW) / (CST%XEPSILO*ZZW) - 1.0
+                                               ! Supersaturation over ice
+  ZUSW=MIN(PPABST/2., ZUSW)            ! safety limitation
+  ZUSW=(ZUSW/ZZW)*((PPABST-ZZW)/(PPABST-ZUSW)) - 1.0
+                             ! Supersaturation of saturated water vapor over ice
+  ZSSI=MIN(ZSSI, ZUSW) ! limitation of SSi according to SSw=0
+END IF
+  
+ZZW=0.
+
+IF(GNEGT) THEN
+  IF(PT<CST%XTT-5.0 .AND. ZSSI>0.0) THEN
+    ZZW=ICEP%XNU20*EXP(ICEP%XALPHA2*ZSSI-ICEP%XBETA2)
+  ELSEIF(PT<=CST%XTT-2.0 .AND. PT>=CST%XTT-5.0 .AND. ZSSI>0.0) THEN
+    ZZW=MAX(ICEP%XNU20*EXP(-ICEP%XBETA2 ), &                                                                                       
+                ICEP%XNU10*EXP(-ICEP%XBETA1*(PT-CST%XTT))*(ZSSI/ZUSW)**ICEP%XALPHA1)
+  ENDIF
+ENDIF
+IF (GNEGT) THEN
+  ZZW=ZZW-PCIT
+  ZZW=MIN(ZZW, 50.E3) ! limitation provisoire a 50 l^-1
+END IF
+
+PRVHENI_MR=0.
+
+IF (GNEGT) THEN
+  PRVHENI_MR=MAX(ZZW, 0.0)*ICEP%XMNU0/PRHODREF
+  PRVHENI_MR=MIN(PRVT, PRVHENI_MR)
+END IF
+!Limitation due to 0 crossing of temperature
+IF(PARAMI%LFEEDBACKT) THEN
+  ZW=0.
+  IF (GNEGT) THEN
+    ZW=MIN(PRVHENI_MR, &
+              MAX(0., (CST%XTT/PEXN-PTHT)/PLSFACT)) / &
+              MAX(PRVHENI_MR, 1.E-20)
+  END IF
+  PRVHENI_MR=PRVHENI_MR*ZW
+  ZZW=ZZW*ZW
+ENDIF
+IF (GNEGT) THEN
+  PCIT=MAX(ZZW+PCIT, PCIT)
+END IF
+END SUBROUTINE ICE4_NUCLEATION
diff --git a/src/common/micro/mode_ice4_nucleation.F90 b/src/common/micro/mode_ice4_nucleation.F90
deleted file mode 100644
index 6453550f45de02343cc4775c7821bfd951347508..0000000000000000000000000000000000000000
--- a/src/common/micro/mode_ice4_nucleation.F90
+++ /dev/null
@@ -1,155 +0,0 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!-----------------------------------------------------------------
-MODULE MODE_ICE4_NUCLEATION
-IMPLICIT NONE
-CONTAINS
-SUBROUTINE ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, KSIZE, ODCOMPUTE, &
-                           PTHT, PPABST, PRHODREF, PEXN, PLSFACT, PT, &
-                           PRVT, &
-                           PCIT, PRVHENI_MR)
-!!
-!!**  PURPOSE
-!!    -------
-!!      Computes the nucleation
-!!
-!!    AUTHOR
-!!    ------
-!!      S. Riette from the splitting of rain_ice source code (nov. 2014)
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!     R. El Khatib 24-Aug-2021 Optimizations
-!
-!
-!*      0. DECLARATIONS
-!          ------------
-!
-USE MODD_CST,            ONLY: CST_t
-USE MODD_PARAM_ICE_n,      ONLY: PARAM_ICE_t
-USE MODD_RAIN_ICE_DESCR_n, ONLY: RAIN_ICE_DESCR_t
-USE MODD_RAIN_ICE_PARAM_n, ONLY: RAIN_ICE_PARAM_t
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of dummy arguments :
-!
-TYPE(CST_t),              INTENT(IN)    :: CST
-TYPE(PARAM_ICE_t),        INTENT(IN)    :: PARAMI
-TYPE(RAIN_ICE_PARAM_t),   INTENT(IN)    :: ICEP
-TYPE(RAIN_ICE_DESCR_t),   INTENT(IN)    :: ICED
-INTEGER,                  INTENT(IN)    :: KSIZE
-LOGICAL, DIMENSION(KSIZE),INTENT(IN)    :: ODCOMPUTE
-REAL, DIMENSION(KSIZE),   INTENT(IN)    :: PTHT    ! Theta at t
-REAL, DIMENSION(KSIZE),   INTENT(IN)    :: PPABST  ! absolute pressure at t
-REAL, DIMENSION(KSIZE),   INTENT(IN)    :: PRHODREF! Reference density
-REAL, DIMENSION(KSIZE),   INTENT(IN)    :: PEXN    ! Exner function
-REAL, DIMENSION(KSIZE),   INTENT(IN)    :: PLSFACT
-REAL, DIMENSION(KSIZE),   INTENT(IN)    :: PT      ! Temperature at time t
-REAL, DIMENSION(KSIZE),   INTENT(IN)    :: PRVT    ! Water vapor m.r. at t
-REAL, DIMENSION(KSIZE),   INTENT(INOUT) :: PCIT    ! Pristine ice n.c. at t
-REAL, DIMENSION(KSIZE),   INTENT(OUT)   :: PRVHENI_MR ! Mixing ratio change due to the heterogeneous nucleation
-!
-!*       0.2  declaration of local variables
-!
-REAL, DIMENSION(KSIZE) :: ZW ! work array
-REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
-LOGICAL, DIMENSION(KSIZE) :: GNEGT  ! Test where to compute the HEN process
-REAL, DIMENSION(KSIZE)  :: ZZW,      & ! Work array
-                           ZUSW,     & ! Undersaturation over water
-                           ZSSI        ! Supersaturation over ice
-INTEGER :: JL
-!-------------------------------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION', 0, ZHOOK_HANDLE)!
-!
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(ODCOMPUTE(:))
-  GNEGT(:)=PT(:)<CST%XTT .AND. PRVT(:)>ICED%XRTMIN(1)
-ELSEWHERE
-  GNEGT(:)=.FALSE.
-ENDWHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
-
-ZUSW(:)=0.
-ZZW(:)=0.
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  ZZW(:)=ALOG(PT(:))
-  ZUSW(:)=EXP(CST%XALPW - CST%XBETAW/PT(:) - CST%XGAMW*ZZW(:))          ! es_w
-  ZZW(:)=EXP(CST%XALPI - CST%XBETAI/PT(:) - CST%XGAMI*ZZW(:))           ! es_i
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
-
-ZSSI(:)=0.
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  ZZW(:)=MIN(PPABST(:)/2., ZZW(:))             ! safety limitation
-  ZSSI(:)=PRVT(:)*(PPABST(:)-ZZW(:)) / (CST%XEPSILO*ZZW(:)) - 1.0
-                                               ! Supersaturation over ice
-  ZUSW(:)=MIN(PPABST(:)/2., ZUSW(:))            ! safety limitation
-  ZUSW(:)=(ZUSW(:)/ZZW(:))*((PPABST(:)-ZZW(:))/(PPABST(:)-ZUSW(:))) - 1.0
-                             ! Supersaturation of saturated water vapor over ice
-  !
-  !*       3.1     compute the heterogeneous nucleation source RVHENI
-  !
-  !*       3.1.1   compute the cloud ice concentration
-  !
-  ZSSI(:)=MIN(ZSSI(:), ZUSW(:)) ! limitation of SSi according to SSw=0
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
-
-ZZW(:)=0.
-DO JL=1,KSIZE
-  IF(GNEGT(JL)) THEN
-    IF(PT(JL)<CST%XTT-5.0 .AND. ZSSI(JL)>0.0) THEN
-      ZZW(JL)=ICEP%XNU20*EXP(ICEP%XALPHA2*ZSSI(JL)-ICEP%XBETA2)
-    ELSEIF(PT(JL)<=CST%XTT-2.0 .AND. PT(JL)>=CST%XTT-5.0 .AND. ZSSI(JL)>0.0) THEN
-      ZZW(JL)=MAX(ICEP%XNU20*EXP(-ICEP%XBETA2 ), &                                                                                       
-                  ICEP%XNU10*EXP(-ICEP%XBETA1*(PT(JL)-CST%XTT))*(ZSSI(JL)/ZUSW(JL))**ICEP%XALPHA1)
-    ENDIF
-  ENDIF
-ENDDO
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  ZZW(:)=ZZW(:)-PCIT(:)
-  ZZW(:)=MIN(ZZW(:), 50.E3) ! limitation provisoire a 50 l^-1
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
-
-PRVHENI_MR(:)=0.
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  !
-  !*       3.1.2   update the r_i and r_v mixing ratios
-  !
-  PRVHENI_MR(:)=MAX(ZZW(:), 0.0)*ICEP%XMNU0/PRHODREF(:)
-  PRVHENI_MR(:)=MIN(PRVT(:), PRVHENI_MR(:))
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
-!Limitation due to 0 crossing of temperature
-IF(PARAMI%LFEEDBACKT) THEN
-  ZW(:)=0.
-  !$mnh_expand_where(JL=1:KSIZE)
-  WHERE(GNEGT(:))
-    ZW(:)=MIN(PRVHENI_MR(:), &
-              MAX(0., (CST%XTT/PEXN(:)-PTHT(:))/PLSFACT(:))) / &
-              MAX(PRVHENI_MR(:), 1.E-20)
-  END WHERE
-  PRVHENI_MR(:)=PRVHENI_MR(:)*ZW(:)
-  ZZW(:)=ZZW(:)*ZW(:)
-  !$mnh_end_expand_where(JL=1:KSIZE)
-ENDIF
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  PCIT(:)=MAX(ZZW(:)+PCIT(:), PCIT(:))
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
-!
-IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION', 1, ZHOOK_HANDLE)
-END SUBROUTINE ICE4_NUCLEATION
-END MODULE MODE_ICE4_NUCLEATION
diff --git a/src/common/micro/mode_ice4_tendencies.F90 b/src/common/micro/mode_ice4_tendencies.F90
index 9395bf08fd0c026ad77a9a7776bfb41764af0d33..816170d7f6a3017deda3a76977b6126f8e898770 100644
--- a/src/common/micro/mode_ice4_tendencies.F90
+++ b/src/common/micro/mode_ice4_tendencies.F90
@@ -54,7 +54,6 @@ USE MODE_ICE4_FAST_RS, ONLY: ICE4_FAST_RS
 USE MODE_ICE4_FAST_RG, ONLY: ICE4_FAST_RG
 USE MODE_ICE4_FAST_RH, ONLY: ICE4_FAST_RH
 USE MODE_ICE4_FAST_RI, ONLY: ICE4_FAST_RI
-USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
 !
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
 !
@@ -135,10 +134,10 @@ ELSE
   !
   !*       2.     COMPUTES THE SLOW COLD PROCESS SOURCES
   !               --------------------------------------
-  CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, KSIZE, LDCOMPUTE(:), &
-                   ZVART(:,ITH), PPRES(:), PRHODREF(:), PEXN(:), PLSFACT(:), ZT(:), &
-                   ZVART(:,IRV), &
-                   PCIT(:), PBU_INST(:, IRVHENI_MR))
+  CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, LDCOMPUTE(1:KSIZE), &
+                   ZVART(1:KSIZE,ITH), PPRES(1:KSIZE), PRHODREF(1:KSIZE), PEXN(1:KSIZE), PLSFACT(1:KSIZE), ZT(1:KSIZE), &
+                   ZVART(1:KSIZE,IRV), &
+                   PCIT(1:KSIZE), PBU_INST(1:KSIZE, IRVHENI_MR))
   DO JL=1, KSIZE
     ZVART(JL,ITH)=ZVART(JL,ITH) + PBU_INST(JL, IRVHENI_MR)*PLSFACT(JL)
     ZT(JL) = ZVART(JL,ITH) * PEXN(JL)
@@ -549,5 +548,7 @@ ENDDO
 !
 IF (LHOOK) CALL DR_HOOK('ICE4_TENDENCIES', 1, ZHOOK_HANDLE)
 !
+CONTAINS
+INCLUDE "ice4_nucleation.func.h"
 END SUBROUTINE ICE4_TENDENCIES
 END MODULE MODE_ICE4_TENDENCIES
diff --git a/src/common/micro/rain_ice.F90 b/src/common/micro/rain_ice.F90
index eae16c20362caf547d0088b7923291d5864743a8..48a652e514305cc6c94ba9a2e0a6ee3a2bf9875b 100644
--- a/src/common/micro/rain_ice.F90
+++ b/src/common/micro/rain_ice.F90
@@ -203,7 +203,6 @@ USE MODE_ICE4_RAINFR_VERT,   ONLY: ICE4_RAINFR_VERT
 USE MODE_ICE4_COMPUTE_PDF,   ONLY: ICE4_COMPUTE_PDF
 USE MODE_ICE4_SEDIMENTATION, ONLY: ICE4_SEDIMENTATION
 USE MODE_ICE4_PACK, ONLY: ICE4_PACK
-USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
 USE MODE_ICE4_CORRECT_NEGATIVITIES, ONLY: ICE4_CORRECT_NEGATIVITIES
 !
 IMPLICIT NONE
@@ -418,7 +417,7 @@ DO JK=IKTB,IKTE
     ENDIF
   ENDDO
 ENDDO
-CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, D%NIJT*D%NKT, LLW3D(:,:), &
+CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, LLW3D(:,:), &
                      PTHT(:, :), PPABST(:, :), PRHODREF(:, :), &                                       
                      PEXN(:, :), ZW3D(:, :), ZT(:, :), &                                                           
                      PRVT(:, :), &                                                                                 
@@ -669,4 +668,6 @@ ENDIF
 
 IF (LHOOK) CALL DR_HOOK('RAIN_ICE', 1, ZHOOK_HANDLE)
 !
+CONTAINS
+INCLUDE "ice4_nucleation.func.h"
 END SUBROUTINE RAIN_ICE