diff --git a/docs/TODO b/docs/TODO
index 00d21af05f14998f9fe9f1cdfdc1030fedc92379..b7014406c76b03d3f859265e0505d800984996ff 100644
--- a/docs/TODO
+++ b/docs/TODO
@@ -14,6 +14,15 @@ Clé de compilation REPRO48 + REPRO55 ajoutées pour permettre de reproduire le
 - modifient l'organisation de calculs
 Ces clés devront être supprimées
 
+Ecrire doc sur marche à suivre pour intégrer un nouveau développement:
+- dev dans MNH à faire en array-syntax
+- dev dans AROME à faire en boucles do
+- intégration dans PHYEX: en array-syntax avec directives mnh_expand
+- les 3 tests suivants doivent donner les mêmes résultats (au bit près) dans chacun des deux modèles:
+  - compilation directe sans activer mnh_expand
+  - compilation en activant mnh_expand
+  - exécution en changeant le nombre de processeurs
+
 Merge pb:
 - ice4_nucleation_wrapper:
        Tableaux allocatable introduits par Philippe dans meso-nh.
@@ -39,15 +48,23 @@ Etape 2: array syntax -> loop
 
 Pb identifiés à corriger plus tard:
 - deposition devrait être déplacée dans ice4_tendencies
-- non reproduction en changeant le nombre de procs
 - avec les optimisations de Ryad, les tableaux 3D de precip passés à ice4_tendencies
   lorsque HSUBG_RC_RR_ACCR=='PRFR' ne sont  pas utilisables puisque les K1, K2 et K3
   sont relatifs à la boucle IMICRO et que les calculs faits en debut de routine ne
   concernent qu'une partie des points
        => à corriger
-- seules les options oper ont été testées, il manque des test pour sedim_after, nmaxiter, xmrstep, xtstep, autoconv, rainfr
+- seules quelques options sont testées avec les cas test (par exemple, il faudrait tester RMC01 mais
+  l'option n'est pas remontée en namelist)
+- les options CMF_CLOUD='STAT' et LOSIGAMS=.FALSE. semblent cassées en 48 original
 - arome/ini_cmfshall devrait s'appeler ini_param_mfshall
 - th_r_from_thl_rt appelée partout, il faudrait limiter à OTEST
+- la recompilation complète d'AROME n'est pas testée
+- il faudrait inclure un cas test plus conséquent en taille au moins sur belenos
+- doute sur le codage de MODD_PRECISION
+- appel à abort à travers print_msg non testé
+- lignes vides ajoutées après les macros mnh_expand
+- indentation inorrecte dans les blocs mnh_expand
+- sedimentation momentum non branchée
 
 Répertoire arome/ext et mesonh/ext contient les codes non PHYEX qu'il faut modifier dans le pack pour qu'il puisse être compilé.
 Ce répertoire devra être vidé à la fin du phasage, les modifications nécessaires ayadevront avoir été fournies par ailleurs
@@ -61,17 +78,12 @@ Budgets/DDH
 SPP
 - modd_spp_type est pour l'instant dans mpa/micro/externals mais n'est pas de la microphysique
 
-Gradients/shuman:
-- essayer de mettre des abort dans les routines arome (shuman doit suffire)
-
 Nettoyage apl_arome non fait (pb a la compilation) ==> 4 arguments dans aro_turb_mnh supprimés (non utilisés)
 
 turb.F90 : il reste un CALL à SOURCES_NEG_CORRECT à ajouter. Besoin de récupérer CCLOUD dans apl_arome : comment ?
 
 Regarder s'il ne serait pas possible/souhaitable de supprimer modd_lunit de PHYEX. On pourrait se contenter de recevoir le numero d'unité logique
 
-Faire quelque chose de mesonh/micro/modd_blankn.f90: le déplacer dans common ou le supprimer
-
 Nettoyage des répertoires aux nécessaire
 
 Initialiser dans AROME la variable ldiag_in_run de MODD_DIAG_IN_RUN pour pouvoir phaser le modd
diff --git a/src/arome/gmkpack_ignored_files b/src/arome/gmkpack_ignored_files
index 193b04bd4624e14fe92f02563a87d775bf722f2d..44cbe8f4057eaef9834c0d43c528aacc49a4eed3 100644
--- a/src/arome/gmkpack_ignored_files
+++ b/src/arome/gmkpack_ignored_files
@@ -101,6 +101,35 @@ phyex/micro/modi_icecloud.F90
 phyex/micro/tiwmx_tab.F90
 phyex/micro/modi_tiwmx.F90
 phyex/micro/modd_spp_type.F90
+phyex/micro/modd_cst.F90
+phyex/micro/modi_ini_cst.F90
+phyex/micro/ini_cst.F90
+phyex/turb/modi_compute_function_thermo_mf.F90
+phyex/turb/compute_function_thermo_mf.F90
+phyex/turb/modd_cmfshall.F90
+phyex/turb/mf_turb_expl.F90
+phyex/turb/modi_mf_turb_expl.F90
+phyex/turb/mf_turb.F90
+phyex/turb/modi_mf_turb.F90
+phyex/turb/compute_mf_cloud.F90
+phyex/turb/compute_mf_cloud_bigaus.F90
+phyex/turb/compute_mf_cloud_direct.F90
+phyex/turb/compute_mf_cloud_stat.F90
+phyex/turb/modi_compute_mf_cloud.F90
+phyex/turb/modi_compute_mf_cloud_bigaus.F90
+phyex/turb/modi_compute_mf_cloud_direct.F90
+phyex/turb/modi_compute_mf_cloud_stat.F90
+phyex/turb/compute_bl89_ml.F90
+phyex/turb/modi_compute_bl89_ml.F90
+phyex/turb/compute_entr_detr.F90
+phyex/turb/modi_compute_entr_detr.F90
+phyex/turb/compute_updraft_raha.F90
+phyex/turb/modi_compute_updraft_raha.F90
+phyex/turb/modi_compute_updraft_rhcj10.F90
+phyex/turb/compute_updraft_rhcj10.F90
+phyex/turb/modi_compute_updraft.F90
+phyex/turb/compute_updraft.F90
+>>>>>>> 9df244f
 phyex/turb/tke_eps_sources.F90
 phyex/turb/turb_ver.F90
 phyex/turb/prandtl.F90
diff --git a/src/common/micro/ice4_nucleation_elem.func.h b/src/common/micro/ice4_nucleation_elem.func.h
new file mode 100644
index 0000000000000000000000000000000000000000..f9ef73991cae536f9e15786a60d153b783da3c69
--- /dev/null
+++ b/src/common/micro/ice4_nucleation_elem.func.h
@@ -0,0 +1,110 @@
+!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_ELEM(ODCOMPUTE, &
+                           PTHT, PPABST, PRHODREF, PEXN, PLSFACT, PT, &
+                           PRVT, &
+                           PCIT, PRVHENI_MR)
+! ******* TO BE INCLUDED IN THE *CONTAINS* OF A SUBROUTINE, IN ORDER TO EASE AUTOMATIC INLINING ******
+! => Don't use drHook !!!
+!
+
+!!
+!!**  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
+!!     S. Riette Feb 2022: as an include file
+!
+!
+!*      0. DECLARATIONS
+!          ------------
+!
+USE MODD_CST,            ONLY: XALPI, XALPW, XBETAI, XBETAW, XGAMI, XGAMW, XMD, XMV, XTT, XEPSILO
+USE MODD_PARAM_ICE,      ONLY: LFEEDBACKT
+USE MODD_RAIN_ICE_PARAM, ONLY: XALPHA1, XALPHA2, XBETA1, XBETA2, XMNU0, XNU10, XNU20
+USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+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 scalar
+         ZUSW,     & ! Undersaturation over water
+         ZSSI        ! Supersaturation over ice
+!-------------------------------------------------------------------------------
+!
+GNEGT=PT<XTT .AND. PRVT>XRTMIN(1) .AND. ODCOMPUTE
+
+PRVHENI_MR=0.
+IF(GNEGT) THEN
+  ZZW=ALOG(PT)
+  ZUSW=EXP(XALPW - XBETAW/PT - XGAMW*ZZW)          ! es_w
+  ZZW=EXP(XALPI - XBETAI/PT - XGAMI*ZZW)           ! es_i
+
+  ZZW=MIN(PPABST/2., ZZW)             ! safety limitation
+  ZSSI=PRVT*(PPABST-ZZW) / (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
+
+  IF(PT<XTT-5. .AND. ZSSI>0.) THEN
+    ZZW=XNU20*EXP(XALPHA2*ZSSI-XBETA2)
+  ELSEIF(PT<=XTT-2. .AND. PT>=XTT-5. .AND. ZSSI>0.) THEN
+    ZZW=MAX(XNU20*EXP(-XBETA2 ), &
+            XNU10*EXP(-XBETA1*(PT-XTT))*(ZSSI/ZUSW)**XALPHA1)
+  ELSE
+    ZZW=0.
+  ENDIF
+
+  ZZW=ZZW-PCIT
+  ZZW=MIN(ZZW, 50.E3) ! limitation provisoire a 50 l^-1
+  !
+  !*       3.1.2   update the r_i and r_v mixing ratios
+  !
+  PRVHENI_MR=MAX(ZZW, 0.0)*XMNU0/PRHODREF
+  PRVHENI_MR=MIN(PRVT, PRVHENI_MR)
+  !
+  !Limitation due to 0 crossing of temperature
+  !
+  IF(LFEEDBACKT) THEN
+    ZW=MIN(PRVHENI_MR, MAX(0., (XTT/PEXN-PTHT)/PLSFACT)) / &
+       MAX(PRVHENI_MR, 1.E-20)
+    PRVHENI_MR=PRVHENI_MR*ZW
+    ZZW=ZZW*ZW
+  ENDIF
+  !
+  PCIT=MAX(ZZW+PCIT, PCIT)
+ENDIF
+!
+END SUBROUTINE ICE4_NUCLEATION_ELEM
diff --git a/src/common/micro/mode_ice4_compute_pdf.F90 b/src/common/micro/mode_ice4_compute_pdf.F90
index 6fb091d6c2c54691f9cad6371574c0f6417049d3..942c65c25291d9ef50dd052af074aa58f90af95b 100644
--- a/src/common/micro/mode_ice4_compute_pdf.F90
+++ b/src/common/micro/mode_ice4_compute_pdf.F90
@@ -309,7 +309,11 @@ ELSE
   CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'ICE4_COMPUTE_PDF', 'wrong HSUBG_AUCV_RI case' )
 ENDIF
 !
+#ifdef REPRO48
+PRF=PHLC_HCF
+#else
 PRF=MAX(PHLC_HCF,PHLI_HCF)
+#endif
 !
 IF (LHOOK) CALL DR_HOOK('ICE4_COMPUTE_PDF', 1, ZHOOK_HANDLE)
 END SUBROUTINE ICE4_COMPUTE_PDF
diff --git a/src/common/micro/mode_ice4_nucleation.F90 b/src/common/micro/mode_ice4_nucleation.F90
deleted file mode 100644
index 7d54233276052d05f4deec25e9ae7c072c6f4b7a..0000000000000000000000000000000000000000
--- a/src/common/micro/mode_ice4_nucleation.F90
+++ /dev/null
@@ -1,127 +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(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: XALPI, XALPW, XBETAI, XBETAW, XGAMI, XGAMW, XMD, XMV, XTT, XEPSILO
-USE MODD_PARAM_ICE,      ONLY: LFEEDBACKT
-USE MODD_RAIN_ICE_PARAM, ONLY: XALPHA1, XALPHA2, XBETA1, XBETA2, XMNU0, XNU10, XNU20
-USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN
-USE PARKIND1, ONLY : JPRB
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of dummy arguments :
-!
-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=JPRB) :: 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
-!-------------------------------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION', 0, ZHOOK_HANDLE)!
-!
-GNEGT(:)=PT(:)<XTT .AND. PRVT(:)>XRTMIN(1) .AND. ODCOMPUTE(:)
-
-ZUSW(:)=0.
-ZZW(:)=0.
-WHERE(GNEGT(:))
-  ZZW(:)=ALOG(PT(:))
-  ZUSW(:)=EXP(XALPW - XBETAW/PT(:) - XGAMW*ZZW(:))          ! es_w
-  ZZW(:)=EXP(XALPI - XBETAI/PT(:) - XGAMI*ZZW(:))           ! es_i
-END WHERE
-
-ZSSI(:)=0.
-WHERE(GNEGT(:))
-  ZZW(:)=MIN(PPABST(:)/2., ZZW(:))             ! safety limitation
-  ZSSI(:)=PRVT(:)*(PPABST(:)-ZZW(:)) / (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
-
-ZZW(:)=0.
-WHERE(GNEGT(:) .AND. PT(:)<XTT-5.0 .AND. ZSSI(:)>0.0 )
-  ZZW(:)=XNU20*EXP(XALPHA2*ZSSI(:)-XBETA2)
-ELSEWHERE(GNEGT(:) .AND. PT(:)<=XTT-2.0 .AND. PT(:)>=XTT-5.0 .AND. ZSSI(:)>0.0)
-  ZZW(:)=MAX(XNU20*EXP(-XBETA2 ), &
-             XNU10*EXP(-XBETA1*(PT(:)-XTT))*(ZSSI(:)/ZUSW(:))**XALPHA1)
-END WHERE
-WHERE(GNEGT(:))
-  ZZW(:)=ZZW(:)-PCIT(:)
-  ZZW(:)=MIN(ZZW(:), 50.E3) ! limitation provisoire a 50 l^-1
-END WHERE
-
-PRVHENI_MR(:)=0.
-WHERE(GNEGT(:))
-  !
-  !*       3.1.2   update the r_i and r_v mixing ratios
-  !
-  PRVHENI_MR(:)=MAX(ZZW(:), 0.0)*XMNU0/PRHODREF(:)
-  PRVHENI_MR(:)=MIN(PRVT(:), PRVHENI_MR(:))
-END WHERE
-!Limitation due to 0 crossing of temperature
-IF(LFEEDBACKT) THEN
-  ZW(:)=0.
-  WHERE(GNEGT(:))
-    ZW(:)=MIN(PRVHENI_MR(:), &
-              MAX(0., (XTT/PEXN(:)-PTHT(:))/PLSFACT(:))) / &
-              MAX(PRVHENI_MR(:), 1.E-20)
-  END WHERE
-  PRVHENI_MR(:)=PRVHENI_MR(:)*ZW(:)
-  ZZW(:)=ZZW(:)*ZW(:)
-ENDIF
-PCIT(:)=MAX(ZZW(:)+PCIT(:), PCIT(:))
-!
-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_nucleation_wrapper.F90 b/src/common/micro/mode_ice4_nucleation_wrapper.F90
deleted file mode 100644
index f16f0225517060c7d5175a9e4f643174b653d7b4..0000000000000000000000000000000000000000
--- a/src/common/micro/mode_ice4_nucleation_wrapper.F90
+++ /dev/null
@@ -1,139 +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_WRAPPER
-IMPLICIT NONE
-CONTAINS
-SUBROUTINE ICE4_NUCLEATION_WRAPPER(KIT, KJT, KKT, LDMASK, &
-                                   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
-!!    -------------
-!  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
-!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
-!!     R. El Khatib 24-Aug-2021 Optimizations
-!
-!
-!*      0. DECLARATIONS
-!          ------------
-!
-USE MODD_CST,   ONLY: XTT
-USE MODE_TOOLS, ONLY: COUNTJV
-USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
-USE PARKIND1,   ONLY : JPRB
-USE YOMHOOK ,   ONLY : LHOOK, DR_HOOK
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of dummy arguments :
-!
-INTEGER,                        INTENT(IN)    :: KIT, KJT, KKT
-LOGICAL, DIMENSION(KIT,KJT,KKT),INTENT(IN)    :: LDMASK
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(IN)    :: PTHT    ! Theta at t
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(IN)    :: PPABST  ! absolute pressure at t
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(IN)    :: PRHODREF! Reference density
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(IN)    :: PEXN    ! Exner function
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(IN)    :: PLSFACT
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(IN)    :: PT      ! Temperature at time t
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(IN)    :: PRVT    ! Water vapor m.r. at t
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(INOUT) :: PCIT    ! Pristine ice n.c. at t
-REAL, DIMENSION(KIT,KJT,KKT),   INTENT(OUT)   :: PRVHENI_MR ! Mixing ratio change due to the heterogeneous nucleation
-!
-!*       0.2  declaration of local variables
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-INTEGER                            :: JL
-INTEGER                            :: INEGT
-INTEGER, DIMENSION(COUNT(PT<XTT .AND. LDMASK)) :: I1,I2,I3
-LOGICAL, DIMENSION(COUNT(PT<XTT .AND. LDMASK))  :: GLDCOMPUTE ! computation criterium
-LOGICAL, DIMENSION(KIT, KJT, KKT)  :: GNEGT  ! Test where to compute the HEN process
-REAL, DIMENSION(COUNT(PT<XTT .AND. LDMASK))  :: ZZT,      & ! Temperature
-                                      ZPRES,      & ! Pressure
-                                      ZRVT,       & ! Water vapor m.r. at t
-                                      ZCIT,       & ! Pristine ice conc. at t
-                                      ZTHT,       & ! Theta at t
-                                      ZRHODREF,   &
-                                      ZEXN,       &
-                                      ZLSFACT,    &
-                                      ZRVHENI_MR
-!! MNH version INTEGER, DIMENSION(:), ALLOCATABLE :: I1,I2,I3
-!! MNH version LOGICAL, DIMENSION(:), ALLOCATABLE :: GLDCOMPUTE
-!! MNH version LOGICAL, DIMENSION(KIT,KJT,KKT)    :: GNEGT  ! Test where to compute the HEN process
-!! MNH version REAL, DIMENSION(:), ALLOCATABLE    :: ZZT,       & ! Temperature
-!! MNH version                                       ZPRES,      & ! Pressure
-!! MNH version                                       ZRVT,       & ! Water vapor m.r. at t
-!! MNH version                                       ZCIT,       & ! Pristine ice conc. at t
-!! MNH version                                       ZTHT,       & ! Theta at t
-!! MNH version                                       ZRHODREF,   &
-!! MNH version                                       ZEXN,       &
-!! MNH version                                       ZLSFACT,    &
-!! MNH version                                       ZRVHENI_MR
-!
-!-------------------------------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION_WRAPPER', 0, ZHOOK_HANDLE)!
-!
-!
-!  optimization by looking for locations where
-!  the temperature is negative only !!!
-!
-GNEGT(:,:,:)=PT(:,:,:)<XTT .AND. LDMASK
-INEGT = COUNT(GNEGT(:,:,:))
-!
-!! MNH version ALLOCATE(GLDCOMPUTE(INEGT))
-!! MNH version ALLOCATE(I1(INEGT),I2(INEGT),I3(INEGT))
-!! MNH version ALLOCATE(ZZT(INEGT))
-!! MNH version ALLOCATE(ZPRES(INEGT))
-!! MNH version ALLOCATE(ZRVT(INEGT))
-!! MNH version ALLOCATE(ZCIT(INEGT))
-!! MNH version ALLOCATE(ZTHT(INEGT))
-!! MNH version ALLOCATE(ZRHODREF(INEGT))
-!! MNH version ALLOCATE(ZEXN(INEGT))
-!! MNH version ALLOCATE(ZLSFACT(INEGT))
-!! MNH version ALLOCATE(ZRVHENI_MR(INEGT))
-!
-IF(INEGT>0) INEGT=COUNTJV(GNEGT(:,:,:), I1(:), I2(:), I3(:))
-!
-PRVHENI_MR(:,:,:)=0.
-IF(INEGT>0) THEN
-  DO JL=1, INEGT
-    ZRVT(JL)=PRVT(I1(JL), I2(JL), I3(JL))
-    ZCIT(JL)=PCIT(I1(JL), I2(JL), I3(JL))
-    ZPRES(JL)=PPABST(I1(JL), I2(JL), I3(JL))
-    ZTHT(JL)=PTHT(I1(JL), I2(JL), I3(JL))
-    ZRHODREF(JL)=PRHODREF(I1(JL), I2(JL), I3(JL))
-    ZEXN(JL)=PEXN(I1(JL), I2(JL), I3(JL))
-    ZLSFACT(JL)=PLSFACT(I1(JL), I2(JL), I3(JL)) / ZEXN(JL)
-    ZZT(JL)=PT(I1(JL), I2(JL), I3(JL))
-    GLDCOMPUTE(JL)=ZZT(JL)<XTT
-  ENDDO
-  CALL ICE4_NUCLEATION(INEGT, GLDCOMPUTE, &
-                       ZTHT, ZPRES, ZRHODREF, ZEXN, ZLSFACT, ZZT, &
-                       ZRVT, &
-                       ZCIT, ZRVHENI_MR)
-  DO JL=1, INEGT
-    PRVHENI_MR(I1(JL), I2(JL), I3(JL)) = ZRVHENI_MR(JL)
-    PCIT      (I1(JL), I2(JL), I3(JL)) = ZCIT      (JL)
-  END DO
-END IF
-!
-!! MNH versionDEALLOCATE(GLDCOMPUTE)
-!! MNH versionDEALLOCATE(I1,I2,I3)
-!! MNH versionDEALLOCATE(ZZT,ZPRES,ZRVT,ZCIT,ZTHT,ZRHODREF,ZEXN,ZLSFACT,ZRVHENI_MR)
-!
-IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION_WRAPPER', 1, ZHOOK_HANDLE)
-
-END SUBROUTINE ICE4_NUCLEATION_WRAPPER
-END MODULE MODE_ICE4_NUCLEATION_WRAPPER
diff --git a/src/common/micro/mode_ice4_tendencies.F90 b/src/common/micro/mode_ice4_tendencies.F90
index 5257c217c765068c84dc461c88338bc317677cf3..cc53a39aa675b8da205f9e6d69c0333133414068 100644
--- a/src/common/micro/mode_ice4_tendencies.F90
+++ b/src/common/micro/mode_ice4_tendencies.F90
@@ -64,7 +64,6 @@ USE MODD_FIELDS_ADDRESS, ONLY : & ! common fields adress
       & IRG,     & ! Graupel
       & IRH        ! Hail
 !
-USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
 USE MODE_ICE4_RRHONG, ONLY: ICE4_RRHONG
 USE MODE_ICE4_RIMLTC, ONLY: ICE4_RIMLTC
 USE MODE_ICE4_RSRIMCG_OLD, ONLY: ICE4_RSRIMCG_OLD
@@ -207,10 +206,13 @@ ELSE
   !
   !*       2.     COMPUTES THE SLOW COLD PROCESS SOURCES
   !               --------------------------------------
-  CALL ICE4_NUCLEATION(KSIZE, LLCOMPUTE, &
-                       ZVART(:,ITH), PPRES, PRHODREF, PEXN, PLSFACT, ZT, &
-                       ZVART(:,IRV), &
-                       PCIT, PRVHENI_MR)
+!DIR$ VECTOR ALWAYS
+  DO CONCURRENT (JL=1:KSIZE)
+    CALL ICE4_NUCLEATION_ELEM(LLCOMPUTE(JL), &
+                     ZVART(JL,ITH), PPRES(JL), PRHODREF(JL), PEXN(JL), PLSFACT(JL), ZT(JL), &
+                     ZVART(JL,IRV), &
+                     PCIT(JL), PRVHENI_MR(JL))
+  ENDDO
   DO JL=1, KSIZE
     ZVART(JL,ITH)=ZVART(JL,ITH) + PRVHENI_MR(JL)*PLSFACT(JL)
     ZT(JL) = ZVART(JL,ITH) * PEXN(JL)
@@ -304,16 +306,6 @@ CALL ICE4_COMPUTE_PDF(KSIZE, HSUBG_AUCV_RC, HSUBG_AUCV_RI, HSUBG_PR_PDF,&
                       PHLI_HCF, PHLI_LCF, PHLI_HRI, PHLI_LRI, ZRAINFR)
 LLRFR=HSUBG_RC_RR_ACCR=='PRFR' .OR. HSUBG_RR_EVAP=='PRFR'
 IF (LLRFR) THEN
-  CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'MODE_ICE4_TENDENCIES', &
-          'LLRFR case broken by optimisation, see comments in mode_ice4_tendencies to knwon why (and how to reapir)....')
-  !Microphyscs was optimized by introducing chunks of KPROMA size
-  !Thus, in ice4_tendencies, the 1D array represent only a fraction of the points where microphisical species are present
-  !We cannot rebuild the entire 3D arrays here, so we cannot call ice4_rainfr_vert here
-  !A solution would be to suppress optimisation in this case by setting KPROMA=KSIZE in rain_ice
-  !Another solution would be to compute column by column?
-  !Another one would be to cut tendencies in 3 parts: before rainfr_vert, rainfr_vert, after rainfr_vert
-
-
   !Diagnostic of precipitation fraction
   PRAINFR(:,:,:) = 0.
   ZRRT3D (:,:,:) = 0.
@@ -323,8 +315,10 @@ IF (LLRFR) THEN
   DO JL=1,KSIZE
     PRAINFR(K1(JL), K2(JL), K3(JL)) = ZRAINFR(JL)
     ZRRT3D (K1(JL), K2(JL), K3(JL)) = ZVART(JL,IRR)
+#ifndef REPRO48
     ZRST3D (K1(JL), K2(JL), K3(JL)) = ZVART(JL,IRS)
     ZRGT3D (K1(JL), K2(JL), K3(JL)) = ZVART(JL,IRG)
+#endif
   END DO
   IF (KRR==7) THEN
     DO JL=1,KSIZE    
@@ -498,5 +492,7 @@ CALL ICE4_FAST_RI(KSIZE, ODSOFT, PCOMPUTE, &
 !
 IF (LHOOK) CALL DR_HOOK('ICE4_TENDENCIES', 1, ZHOOK_HANDLE)
 !
+CONTAINS
+INCLUDE "ice4_nucleation_elem.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 0fa4b7bd73f38bfa24aba24506b7dbdb5f14ec07..5da4fe37e8978772295b1dcae8f80afe255c519d 100644
--- a/src/common/micro/rain_ice.F90
+++ b/src/common/micro/rain_ice.F90
@@ -205,7 +205,6 @@ USE MODE_ICE4_RAINFR_VERT, ONLY: ICE4_RAINFR_VERT
 USE MODE_ICE4_SEDIMENTATION_STAT, ONLY: ICE4_SEDIMENTATION_STAT
 USE MODE_ICE4_SEDIMENTATION_SPLIT, ONLY: ICE4_SEDIMENTATION_SPLIT
 USE MODE_ICE4_SEDIMENTATION_SPLIT_MOMENTUM, ONLY: ICE4_SEDIMENTATION_SPLIT_MOMENTUM
-USE MODE_ICE4_NUCLEATION_WRAPPER, ONLY: ICE4_NUCLEATION_WRAPPER
 USE MODE_ICE4_TENDENCIES, ONLY: ICE4_TENDENCIES
 !
 IMPLICIT NONE
@@ -318,7 +317,7 @@ REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2)) :: ZINPRI ! Pristine ice instant prec
 LOGICAL :: GEXT_TEND
 LOGICAL :: LSOFT ! Must we really compute tendencies or only adjust them to new T variables
 INTEGER :: INB_ITER_MAX ! Maximum number of iterations (with real tendencies computation)
-REAL :: ZW1D
+REAL :: ZW0D
 REAL :: ZTSTEP ! length of sub-timestep in case of time splitting
 REAL :: ZINV_TSTEP ! Inverse ov PTSTEP
 REAL :: ZTIME_THRESHOLD ! Time to reach threshold
@@ -1111,10 +1110,21 @@ PCIT(:,:,:)=ZCITOUT(:,:,:)
 !*       6.     COMPUTES THE SLOW COLD PROCESS SOURCES OUTSIDE OF ODMICRO POINTS
 !               ----------------------------------------------------------------
 !
-CALL ICE4_NUCLEATION_WRAPPER(KIT, KJT, KKT, .NOT. ODMICRO, &
-                             PTHT, PPABST, PRHODREF, PEXN, ZZ_LSFACT, ZT, &
-                             PRVT, &
-                             PCIT, ZZ_RVHENI_MR)
+DO JK=1, KKT
+  DO JJ=1, KJT
+!DIR$ VECTOR ALWAYS
+    DO CONCURRENT (JI=1:KIT)
+      IF (.NOT. ODMICRO(JI, JJ, JK)) THEN
+        ZW0D=ZZ_LSFACT(JI, JJ, JK)/PEXN(JI, JJ, JK)
+      ENDIF
+      CALL ICE4_NUCLEATION_ELEM(.NOT. ODMICRO(JI, JJ, JK), &
+                                PTHT(JI, JJ, JK), PPABST(JI, JJ, JK), PRHODREF(JI, JJ, JK), &
+                                PEXN(JI, JJ, JK), ZW0D, ZT(JI, JJ, JK), &
+                                PRVT(JI, JJ, JK), &
+                                PCIT(JI, JJ, JK), ZZ_RVHENI_MR(JI, JJ, JK))
+    ENDDO
+  ENDDO
+ENDDO
 !
 !-------------------------------------------------------------------------------
 !
@@ -1766,6 +1776,7 @@ CONTAINS
   IF (LHOOK) CALL DR_HOOK('RAIN_ICE:CORRECT_NEGATIVITIES', 1, ZHOOK_HANDLE)
   !
   END SUBROUTINE CORRECT_NEGATIVITIES
-
+!
+INCLUDE "ice4_nucleation_elem.func.h"
 !
 END SUBROUTINE RAIN_ICE
diff --git a/src/common/turb/mode_compute_entr_detr.F90 b/src/common/turb/mode_compute_entr_detr.F90
index beb77265bf68b4538191a2e7f5f47d6940aa96d8..7d8fae529856f3a50b886d7d64389b5b8f30f7e2 100644
--- a/src/common/turb/mode_compute_entr_detr.F90
+++ b/src/common/turb/mode_compute_entr_detr.F90
@@ -78,7 +78,7 @@ USE MODD_CST
 USE MODD_PARAM_MFSHALL_n
 !
 USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D 
-
+!
 USE MODE_THERMO
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
diff --git a/src/mesonh/micro/modd_blankn.f90 b/src/mesonh/micro/modd_blankn.f90
deleted file mode 100644
index 6428103136f77d7639c070c7032add80316721f5..0000000000000000000000000000000000000000
--- a/src/mesonh/micro/modd_blankn.f90
+++ /dev/null
@@ -1,173 +0,0 @@
-!MNH_LIC Copyright 1996-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 MODD_BLANK_n
-!     #################
-!
-!!****  *MODD_BLANK$n* -  Declarative module for MesoNH developpers namelist
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!      Offer dummy real, integer, logical and character variables for
-!!    test and debugging purposes.
-!!
-!!**  METHOD
-!!    ------
-!!
-!!      Eight dummy real, integer, logical and character*80 variables are
-!!    defined and passed through the namelist read operations. None of the
-!!    MesoNH routines uses any of those variables. When a developper choses
-!!    to introduce temporarily a parameter to some subroutine, he has to
-!!    introduce a USE MODD_BLANK statement into that subroutine. Then he
-!!    can use any of the variables defined here and change them easily via
-!!    the namelist input.
-!!
-!!    REFERENCE
-!!    ---------
-!!      None
-!!
-!!    AUTHOR
-!!    ------
-!!      K. Suhre   *Laboratoire d'Aerologie*
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!    Original 25/04/96
-!!      updated     17/11/00  (P Jabouille) Use dummy array
-!!      updated     26/10/21  (Q.Rodier) Use for n model (grid-nesting)
-!-------------------------------------------------------------------------------
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-USE MODD_PARAMETERS, ONLY : JPDUMMY, JPMODELMAX
-!
-IMPLICIT NONE
-!
-TYPE BLANK_t
-!
-  LOGICAL :: LDUMMY1
-  LOGICAL :: LDUMMY2
-  LOGICAL :: LDUMMY3
-  LOGICAL :: LDUMMY4
-  LOGICAL :: LDUMMY5
-  LOGICAL :: LDUMMY6
-  LOGICAL :: LDUMMY7
-  LOGICAL :: LDUMMY8
-!
-  CHARACTER(len=80) :: CDUMMY1
-  CHARACTER(len=80) :: CDUMMY2
-  CHARACTER(len=80) :: CDUMMY3
-  CHARACTER(len=80) :: CDUMMY4
-  CHARACTER(len=80) :: CDUMMY5
-  CHARACTER(len=80) :: CDUMMY6
-  CHARACTER(len=80) :: CDUMMY7
-  CHARACTER(len=80) :: CDUMMY8
-!
-  INTEGER :: NDUMMY1
-  INTEGER :: NDUMMY2
-  INTEGER :: NDUMMY3
-  INTEGER :: NDUMMY4
-  INTEGER :: NDUMMY5
-  INTEGER :: NDUMMY6
-  INTEGER :: NDUMMY7
-  INTEGER :: NDUMMY8
-!
-  REAL :: XDUMMY1
-  REAL :: XDUMMY2
-  REAL :: XDUMMY3
-  REAL :: XDUMMY4
-  REAL :: XDUMMY5
-  REAL :: XDUMMY6
-  REAL :: XDUMMY7
-  REAL :: XDUMMY8
-!
-END TYPE BLANK_t
-!
-TYPE(BLANK_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: BLANK_MODEL
-!
-LOGICAL, POINTER :: LDUMMY1=>NULL()
-LOGICAL, POINTER :: LDUMMY2=>NULL()
-LOGICAL, POINTER :: LDUMMY3=>NULL()
-LOGICAL, POINTER :: LDUMMY4=>NULL()
-LOGICAL, POINTER :: LDUMMY5=>NULL()
-LOGICAL, POINTER :: LDUMMY6=>NULL()
-LOGICAL, POINTER :: LDUMMY7=>NULL()
-LOGICAL, POINTER :: LDUMMY8=>NULL()
-!
-CHARACTER(len=80), POINTER :: CDUMMY1=>NULL()
-CHARACTER(len=80), POINTER :: CDUMMY2=>NULL()
-CHARACTER(len=80), POINTER :: CDUMMY3=>NULL()
-CHARACTER(len=80), POINTER :: CDUMMY4=>NULL()
-CHARACTER(len=80), POINTER :: CDUMMY5=>NULL()
-CHARACTER(len=80), POINTER :: CDUMMY6=>NULL()
-CHARACTER(len=80), POINTER :: CDUMMY7=>NULL()
-CHARACTER(len=80), POINTER :: CDUMMY8=>NULL()
-!
-INTEGER, POINTER :: NDUMMY1=>NULL()
-INTEGER, POINTER :: NDUMMY2=>NULL()
-INTEGER, POINTER :: NDUMMY3=>NULL()
-INTEGER, POINTER :: NDUMMY4=>NULL()
-INTEGER, POINTER :: NDUMMY5=>NULL()
-INTEGER, POINTER :: NDUMMY6=>NULL()
-INTEGER, POINTER :: NDUMMY7=>NULL()
-INTEGER, POINTER :: NDUMMY8=>NULL()
-!
-REAL, POINTER :: XDUMMY1=>NULL()
-REAL, POINTER :: XDUMMY2=>NULL()
-REAL, POINTER :: XDUMMY3=>NULL()
-REAL, POINTER :: XDUMMY4=>NULL()
-REAL, POINTER :: XDUMMY5=>NULL()
-REAL, POINTER :: XDUMMY6=>NULL()
-REAL, POINTER :: XDUMMY7=>NULL()
-REAL, POINTER :: XDUMMY8=>NULL()
-!
-CONTAINS
-!
-SUBROUTINE BLANK_GOTO_MODEL(KFROM,KTO)
-INTEGER, INTENT(IN) :: KFROM, KTO
-!    
-LDUMMY1=>BLANK_MODEL(KTO)%LDUMMY1
-LDUMMY2=>BLANK_MODEL(KTO)%LDUMMY2
-LDUMMY3=>BLANK_MODEL(KTO)%LDUMMY3
-LDUMMY4=>BLANK_MODEL(KTO)%LDUMMY4
-LDUMMY5=>BLANK_MODEL(KTO)%LDUMMY5
-LDUMMY6=>BLANK_MODEL(KTO)%LDUMMY6
-LDUMMY7=>BLANK_MODEL(KTO)%LDUMMY7
-LDUMMY8=>BLANK_MODEL(KTO)%LDUMMY8
-
-CDUMMY1=>BLANK_MODEL(KTO)%CDUMMY1
-CDUMMY2=>BLANK_MODEL(KTO)%CDUMMY2
-CDUMMY3=>BLANK_MODEL(KTO)%CDUMMY3
-CDUMMY4=>BLANK_MODEL(KTO)%CDUMMY4
-CDUMMY5=>BLANK_MODEL(KTO)%CDUMMY5
-CDUMMY6=>BLANK_MODEL(KTO)%CDUMMY6
-CDUMMY7=>BLANK_MODEL(KTO)%CDUMMY7
-CDUMMY8=>BLANK_MODEL(KTO)%CDUMMY8
-!
-NDUMMY1=>BLANK_MODEL(KTO)%NDUMMY1
-NDUMMY2=>BLANK_MODEL(KTO)%NDUMMY2
-NDUMMY3=>BLANK_MODEL(KTO)%NDUMMY3
-NDUMMY4=>BLANK_MODEL(KTO)%NDUMMY4
-NDUMMY5=>BLANK_MODEL(KTO)%NDUMMY5
-NDUMMY6=>BLANK_MODEL(KTO)%NDUMMY6
-NDUMMY7=>BLANK_MODEL(KTO)%NDUMMY7
-NDUMMY8=>BLANK_MODEL(KTO)%NDUMMY8
-!
-XDUMMY1=>BLANK_MODEL(KTO)%XDUMMY1
-XDUMMY2=>BLANK_MODEL(KTO)%XDUMMY2
-XDUMMY3=>BLANK_MODEL(KTO)%XDUMMY3
-XDUMMY4=>BLANK_MODEL(KTO)%XDUMMY4
-XDUMMY5=>BLANK_MODEL(KTO)%XDUMMY5
-XDUMMY6=>BLANK_MODEL(KTO)%XDUMMY6
-XDUMMY7=>BLANK_MODEL(KTO)%XDUMMY7
-XDUMMY8=>BLANK_MODEL(KTO)%XDUMMY8
-!
-END SUBROUTINE BLANK_GOTO_MODEL
-!
-END MODULE MODD_BLANK_n