From b993c1dda83271535a38471efac0b83e29c489aa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Riette?= <sebastien.riette@meteo.fr>
Date: Wed, 18 Jan 2023 12:54:17 +0100
Subject: [PATCH] S. Riette 18 Jan 2022: ice4_stepping introduction

---
 docs/TODO                                 |   4 +-
 src/common/micro/mode_ice4_pack.F90       | 355 ++----------------
 src/common/micro/mode_ice4_stepping.F90   | 421 ++++++++++++++++++++++
 src/common/micro/mode_ice4_tendencies.F90 |   2 +-
 src/common/micro/rain_ice.F90             |   2 +-
 5 files changed, 458 insertions(+), 326 deletions(-)
 create mode 100644 src/common/micro/mode_ice4_stepping.F90

diff --git a/docs/TODO b/docs/TODO
index 10508c7a3..419bcb0ab 100644
--- a/docs/TODO
+++ b/docs/TODO
@@ -48,9 +48,7 @@ Pb identifiés à corriger plus tard:
   - th_r_from_thl_rt appelée partout, il faudrait limiter à OTEST
   - La taille du buffer utilisé pour th_r_from_thl_rt doit être mise en module et
     utilisée pour déclarer le buffer dans les routines appelantes et dans th_r_from_thl_rt
-  - rain_ice:
-    - séparer l'avance temporelle du découpage en sous-blocs en créant une couche driver supplémentaire. Cette couche pourrait
-      avoir différentes implémentations (filtre LLMICRO seul, filtre LLMICRO + découpage en sous-blocs, filtre LLMICR
+  - rain_ice: coder d'autres implémentations de ice4_pack (filtre LLMICRO seul et/ou découpage en sous-blocs, rien)
   - shuman, turb, shuman_mf, apl_arome
     - peut-on remettre à jour la partie sous le sol des variables pronostiques avant d'entrer dans EDKFi dans Méso-NH? Si oui, on devrait pouvoir
       récrire les shuman_mf pour qu'ils produisent les mêmes résultats qu'aujourd'hui mais sans utiliser les valeurs sous le sol ou dans "l'espace infini"
diff --git a/src/common/micro/mode_ice4_pack.F90 b/src/common/micro/mode_ice4_pack.F90
index 4a880db26..486cd66ce 100644
--- a/src/common/micro/mode_ice4_pack.F90
+++ b/src/common/micro/mode_ice4_pack.F90
@@ -44,14 +44,11 @@ USE MODD_FIELDS_ADDRESS, ONLY : & ! common fields adress
       & IRG,     & ! Graupel
       & IRH,     & ! Hail
       & IBUNUM,       & ! Number of tendency terms
-      & IBUNUM_MR,    & ! Number of tendency terms expressed as mixing ratio changes
-      & IBUNUM_EXTRA, & ! Number of extra tendency terms
-      & IRREVAV,      & ! Index for the evaporation tendency
-      & IBUEXTRAIND     ! Index indirection
+      & IBUNUM_EXTRA    ! Number of extra tendency terms
 
 USE MODE_MSG,            ONLY: PRINT_MSG, NVERB_FATAL
 
-USE MODE_ICE4_TENDENCIES, ONLY: ICE4_TENDENCIES
+USE MODE_ICE4_STEPPING, ONLY: ICE4_STEPPING
 USE MODE_ICE4_BUDGETS, ONLY: ICE4_BUDGETS
 !
 IMPLICIT NONE
@@ -124,71 +121,39 @@ INTEGER :: IKTB, IKTE, IIJB, IIJE
 INTEGER :: ISTIJ, ISTK
 !
 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 :: ZTSTEP ! length of sub-timestep in case of time splitting
-REAL :: ZINV_TSTEP ! Inverse ov PTSTEP
-REAL :: ZTIME_THRESHOLD ! Time to reach threshold
 !
 !Output packed total mixing ratio change (for budgets only)
 REAL, DIMENSION(KSIZE, IBUNUM-IBUNUM_EXTRA) :: ZBU_PACK
 !
 !For packing
 INTEGER :: IMICRO ! Case r_x>0 locations
-INTEGER :: JL, JV, JJV
-REAL, DIMENSION(KPROMA) :: ZTIME ! Current integration time (starts with 0 and ends with PTSTEP)
+INTEGER :: JL, JV
 REAL, DIMENSION(KPROMA) :: &
-                        & ZMAXTIME, & ! Time on which we can apply the current tendencies
-                        & ZTIME_LASTCALL, &     ! Integration time when last tendecies call has been done
-                        & ZSSI,     &
                         & ZCIT,     & ! Pristine ice conc. at t
                         & ZRHODREF, & ! RHO Dry REFerence
-                        & ZZT,      & ! Temperature
                         & ZPRES,    & ! Pressure
                         & ZEXN,     & ! EXNer Pressure
-                        & ZLSFACT,  & ! L_s/(Pi*C_ph)
-                        & ZLVFACT,  & ! L_v/(Pi*C_ph)
                         & ZSIGMA_RC,& ! Standard deviation of rc at time t
                         & ZCF,      & ! Cloud fraction
                         & ZHLC_HCF, & ! HLCLOUDS : fraction of High Cloud Fraction in grid
-                        & ZHLC_LCF, & ! HLCLOUDS : fraction of Low  Cloud Fraction in grid
-                                      !    note that ZCF = ZHLC_HCF + ZHLC_LCF
                         & ZHLC_HRC, & ! HLCLOUDS : LWC that is High LWC in grid
-                        & ZHLC_LRC, & ! HLCLOUDS : LWC that is Low  LWC in grid
-                                      !    note that ZRC = ZHLC_HRC + ZHLC_LRC
                         & ZHLI_HCF, &
-                        & ZHLI_LCF, &
                         & ZHLI_HRI, &
-                        & ZHLI_LRI
-LOGICAL, DIMENSION(KPROMA) :: LLCOMPUTE ! .TRUE. or points where we must compute tendencies,
+                        & ZRREVAV
 !
 !Output packed tendencies (for budgets only)
 REAL, DIMENSION(KPROMA, IBUNUM-IBUNUM_EXTRA) :: ZBU_SUM
-REAL, DIMENSION(KPROMA, IBUNUM) :: ZBU_INST
 !
 !For mixing-ratio-splitting
-LOGICAL :: LLCPZ0RT
-REAL :: ZTIME_THRESHOLD1D(KPROMA) ! Time to reach threshold
-REAL, DIMENSION(KPROMA, KRR) :: Z0RT ! Mixing-ratios at the beginig of the current loop
-!
 REAL, DIMENSION(KPROMA,0:7) :: &
                         & ZVART, & !Packed variables
-                        & ZEXTPK, & !To take into acount external tendencies inside the splitting
-                        & ZA, ZB
+                        & ZEXTPK   !To take into acount external tendencies inside the splitting
 !
-REAL, DIMENSION(KPROMA, 8) :: ZRS_TEND, ZRG_TEND
-REAL, DIMENSION(KPROMA,10) :: ZRH_TEND
-
-INTEGER, DIMENSION(KPROMA) :: &
-                       & I1,I2, & ! Used to replace the COUNT and PACK intrinsics on variables
-                       & IITER    ! Number of iterations done (with real tendencies computation)
+INTEGER, DIMENSION(KPROMA) :: I1,I2 ! Used to replace the COUNT and PACK intrinsics on variables
 INTEGER, DIMENSION(KSIZE) :: I1TOT, I2TOT ! Used to replace the COUNT and PACK intrinsics
 !
-REAL, DIMENSION(KPROMA) :: ZSUM2, ZMAXB
-REAL :: ZDEVIDE, ZX
-!
 INTEGER :: IC, JMICRO
-LOGICAL :: LLSIGMA_RC, LL_ANY_ITER, LL_AUCV_ADJU
+LOGICAL :: LLSIGMA_RC, LL_AUCV_ADJU
 !
 !-------------------------------------------------------------------------------
 IF (LHOOK) CALL DR_HOOK('ICE4_PACK', 0, ZHOOK_HANDLE)
@@ -201,7 +166,8 @@ IKTE=D%NKTE
 IIJB=D%NIJB
 IIJE=D%NIJE
 GEXT_TEND=.TRUE.
-ZINV_TSTEP=1./PTSTEP
+LLSIGMA_RC=(HSUBG_AUCV_RC=='PDF ' .AND. PARAMI%CSUBG_PR_PDF=='SIGM')
+LL_AUCV_ADJU=(HSUBG_AUCV_RC=='ADJU' .OR. HSUBG_AUCV_RI=='ADJU')
 !
 IF(KPROMA /= KSIZE) THEN
   CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'RAIN_ICE', 'For now, KPROMA must be equal to KSIZE, see code for explanation')
@@ -216,12 +182,11 @@ ENDIF
 IF(BUCONF%LBU_ENABLE) THEN
   DO JV=1, IBUNUM-IBUNUM_EXTRA
     ZBU_PACK(:, JV)=0.
-    ZBU_SUM(:, JV)=0.
   ENDDO
 ENDIF
-!-------------------------------------------------------------------------------
 !
-!***       4.1 Point selection
+!*       2.     POINT SELECTION
+!               ---------------
 !
 !  optimization by looking for locations where
 !  the microphysical fields are larger than a minimal value only !!!
@@ -231,21 +196,9 @@ IF (KSIZE /= COUNT(ODMICRO(IIJB:IIJE,IKTB:IKTE))) THEN
 ENDIF
 
 IF (KSIZE > 0) THEN
-
-  !Maximum number of iterations
-  !We only count real iterations (those for which we *compute* tendencies)
-  INB_ITER_MAX=PARAMI%NMAXITER
-  IF(PARAMI%XTSTEP_TS/=0.)THEN
-    INB_ITER_MAX=MAX(1, INT(PTSTEP/PARAMI%XTSTEP_TS)) !At least the number of iterations needed for the time-splitting
-    ZTSTEP=PTSTEP/INB_ITER_MAX
-    INB_ITER_MAX=MAX(PARAMI%NMAXITER, INB_ITER_MAX) !For the case XMRSTEP/=0. at the same time
-  ENDIF
-  LLSIGMA_RC=(HSUBG_AUCV_RC=='PDF ' .AND. PARAMI%CSUBG_PR_PDF=='SIGM')
-  LL_AUCV_ADJU=(HSUBG_AUCV_RC=='ADJU' .OR. HSUBG_AUCV_RI=='ADJU')
-
-  !-------------------------------------------------------------------------------
   !
-  !***       4.2 Cache-blocking loop
+  !*       3.     CACHE-BLOCKING LOOP
+  !               -------------------
   !
 
   ! starting indexes :
@@ -256,9 +209,9 @@ IF (KSIZE > 0) THEN
   DO JMICRO=1,KSIZE,KPROMA
 
     IMICRO=MIN(KPROMA,KSIZE-JMICRO+1)
-    !-------------------------------------------------------------------------------
     !
-    !***       4.3 Packing
+    !*       4.     PACKING
+    !               -------
     !
 
     ! Setup packing parameters
@@ -333,269 +286,30 @@ IF (KSIZE > 0) THEN
       ISTIJ=IIJB
     ENDDO OUTER_LOOP
 
-    IF (GEXT_TEND) THEN
-      DO JV=0, KRR
-        DO JL=1, IMICRO
-          ZEXTPK(JL, JV)=ZEXTPK(JL, JV)-ZVART(JL, JV)*ZINV_TSTEP
-        ENDDO
-      ENDDO
-    ENDIF
-    IF (LLSIGMA_RC) THEN
-      DO JL=1, IMICRO
-        ZSIGMA_RC(JL)=ZSIGMA_RC(JL)*2.
-      ENDDO 
-    ENDIF
-    IF (LL_AUCV_ADJU) THEN
-      DO JL=1, IMICRO
-        ZHLC_LRC(JL) = ZVART(JL, IRC) - ZHLC_HRC(JL)
-        ZHLI_LRI(JL) = ZVART(JL, IRI) - ZHLI_HRI(JL)
-        IF(ZVART(JL, IRC)>0.) THEN
-          ZHLC_LCF(JL) = ZCF(JL)- ZHLC_HCF(JL)
-        ELSE
-          ZHLC_LCF(JL)=0.
-        ENDIF
-        IF(ZVART(JL, IRI)>0.) THEN
-          ZHLI_LCF(JL) = ZCF(JL)- ZHLI_HCF(JL)
-        ELSE
-          ZHLI_LCF(JL)=0.
-        ENDIF
-      ENDDO
-    ENDIF
-
-    !-------------------------------------------------------------------------------
-    !
-    !***       4.4 temporal loop
     !
+    !*       5.     TENDENCIES COMPUTATION
+    !               ----------------------
     !
-    IITER(1:IMICRO)=0
-    ZTIME(1:IMICRO)=0. ! Current integration time (all points may have a different integration time)
-
-    DO WHILE(ANY(ZTIME(1:IMICRO)<PTSTEP)) ! Loop to *really* compute tendencies
-
-      IF(PARAMI%XTSTEP_TS/=0.) THEN
-        ! In this case we need to remember the time when tendencies were computed
-        ! because when time has evolved more than a limit, we must re-compute tendencies
-        ZTIME_LASTCALL(1:IMICRO)=ZTIME(1:IMICRO)
-      ENDIF
-      DO JL=1, IMICRO
-        IF (ZTIME(JL) < PTSTEP) THEN
-          LLCOMPUTE(JL)=.TRUE. ! Computation (.TRUE.) only for points for which integration time has not reached the timestep
-          IITER(JL)=IITER(JL)+1
-        ELSE
-          LLCOMPUTE(JL)=.FALSE.
-        ENDIF
-      ENDDO
-      LL_ANY_ITER=ANY(IITER(1:IMICRO) < INB_ITER_MAX)
-      LLCPZ0RT=.TRUE.
-      LSOFT=.FALSE. ! We *really* compute the tendencies
-
-      DO WHILE(ANY(LLCOMPUTE(1:IMICRO))) ! Loop to adjust tendencies when we cross the 0°C or when a species disappears
-!$OMP SIMD
-        DO JL=1, IMICRO
-          ZSUM2(JL)=SUM(ZVART(JL,IRI:KRR))
-        ENDDO
-        DO JL=1, IMICRO
-          ZDEVIDE=(CST%XCPD + CST%XCPV*ZVART(JL, IRV) + CST%XCL*(ZVART(JL, IRC)+ZVART(JL, IRR)) + CST%XCI*ZSUM2(JL)) * ZEXN(JL)
-          ZZT(JL) = ZVART(JL, ITH) * ZEXN(JL)
-          ZLSFACT(JL)=(CST%XLSTT+(CST%XCPV-CST%XCI)*(ZZT(JL)-CST%XTT)) / ZDEVIDE
-          ZLVFACT(JL)=(CST%XLVTT+(CST%XCPV-CST%XCL)*(ZZT(JL)-CST%XTT)) / ZDEVIDE
-        ENDDO
-        !-------------------------------------------------------------------------------
-        !
-        !***       4.5 Effective tendencies computation
-        !
-        !
-        ! Tendencies are *really* computed when LSOFT==.FALSE. and only adjusted otherwise
-        CALL ICE4_TENDENCIES(D, CST, PARAMI, ICEP, ICED, BUCONF, &
-                            &KPROMA, IMICRO, &
-                            &KRR, LSOFT, LLCOMPUTE, &
-                            &HSUBG_AUCV_RC, HSUBG_AUCV_RI, &
-                            &ZEXN, ZRHODREF, ZLVFACT, ZLSFACT, I1, I2, &
-                            &ZPRES, ZCF, ZSIGMA_RC, &
-                            &ZCIT, &
-                            &ZZT, ZVART, &
-                            &ZBU_INST, &
-                            &ZRS_TEND, ZRG_TEND, ZRH_TEND, ZSSI, &
-                            &ZA, ZB, &
-                            &ZHLC_HCF, ZHLC_LCF, ZHLC_HRC, ZHLC_LRC, &
-                            &ZHLI_HCF, ZHLI_LCF, ZHLI_HRI, ZHLI_LRI, PRAINFR)
-
-        ! External tendencies
-        IF(GEXT_TEND) THEN
-          DO JV=0, KRR
-            DO JL=1, IMICRO
-              ZA(JL, JV) = ZA(JL, JV) + ZEXTPK(JL, JV)
-            ENDDO
-          ENDDO
-        ENDIF
-        !-------------------------------------------------------------------------------
-        !
-        !***       4.6 Time integration
-        !
-        !
-        ! If we can, we shall use these tendencies until the end of the timestep
-        DO JL=1, IMICRO
-          IF(LLCOMPUTE(JL)) THEN
-            ZMAXTIME(JL)=(PTSTEP-ZTIME(JL)) ! Remaining time until the end of the timestep
-          ELSE
-            ZMAXTIME(JL)=0.
-          ENDIF
-        ENDDO
-
-        !We need to adjust tendencies when temperature reaches 0
-        IF(PARAMI%LFEEDBACKT) THEN
-          DO JL=1, IMICRO
-            !Is ZB(:, ITH) enough to change temperature sign?
-            ZX=CST%XTT/ZEXN(JL)
-            IF ((ZVART(JL, ITH) - ZX) * (ZVART(JL, ITH) + ZB(JL, ITH) - ZX) < 0.) THEN
-              ZMAXTIME(JL)=0.
-            ENDIF
-            !Can ZA(:, ITH) make temperature change of sign?
-            IF (ABS(ZA(JL,ITH)) > 1.E-20 ) THEN
-              ZTIME_THRESHOLD=(ZX - ZB(JL, ITH) - ZVART(JL, ITH))/ZA(JL, ITH)
-              IF (ZTIME_THRESHOLD > 0.) THEN
-                ZMAXTIME(JL)=MIN(ZMAXTIME(JL), ZTIME_THRESHOLD)
-              ENDIF
-            ENDIF
-          ENDDO
-        ENDIF
-
-        !We need to adjust tendencies when a species disappears
-        !When a species is missing, only the external tendencies can be negative (and we must keep track of it)
-        DO JV=1, KRR
-          DO JL=1, IMICRO
-            IF (ZA(JL, JV) < -1.E-20 .AND. ZVART(JL, JV) > ICED%XRTMIN(JV)) THEN
-              ZMAXTIME(JL)=MIN(ZMAXTIME(JL), -(ZB(JL, JV)+ZVART(JL, JV))/ZA(JL, JV))
-            ENDIF
-          ENDDO
-        ENDDO
-
-        !We stop when the end of the timestep is reached
-        DO JL=1, IMICRO
-          IF (ZTIME(JL)+ZMAXTIME(JL) >= PTSTEP) THEN
-            LLCOMPUTE(JL)=.FALSE.
-          ENDIF
-        ENDDO
-        !We must recompute tendencies when the end of the sub-timestep is reached
-        IF (PARAMI%XTSTEP_TS/=0.) THEN
-          DO JL=1, IMICRO
-            IF ((IITER(JL) < INB_ITER_MAX) .AND. (ZTIME(JL)+ZMAXTIME(JL) > ZTIME_LASTCALL(JL)+ZTSTEP)) THEN
-              ZMAXTIME(JL)=ZTIME_LASTCALL(JL)-ZTIME(JL)+ZTSTEP
-              LLCOMPUTE(JL)=.FALSE.
-            ENDIF
-          ENDDO
-        ENDIF
-
-        !We must recompute tendencies when the maximum allowed change is reached
-        !When a species is missing, only the external tendencies can be active and we do not want to recompute
-        !the microphysical tendencies when external tendencies are negative (results won't change because species was already missing)
-        IF (PARAMI%XMRSTEP/=0.) THEN
-          IF (LL_ANY_ITER) THEN
-            ! In this case we need to remember the initial mixing ratios used to compute the tendencies
-            ! because when mixing ratio has evolved more than a threshold, we must re-compute tendencies
-            ! Thus, at first iteration (ie when LLCPZ0RT=.TRUE.) we copy ZVART into Z0RT
-            DO JV=1,KRR
-              IF (LLCPZ0RT) Z0RT(1:IMICRO, JV)=ZVART(1:IMICRO, JV)
-              DO JL=1, IMICRO
-                IF (IITER(JL)<INB_ITER_MAX .AND. ABS(ZA(JL,JV))>1.E-20) THEN
-                  ZTIME_THRESHOLD1D(JL)=(SIGN(1., ZA(JL, JV))*PARAMI%XMRSTEP+ &
-                                        &Z0RT(JL, JV)-ZVART(JL, JV)-ZB(JL, JV))/ZA(JL, JV)
-                ELSE
-                  ZTIME_THRESHOLD1D(JL)=-1.
-                ENDIF
-              ENDDO
-              DO JL=1, IMICRO
-                IF (ZTIME_THRESHOLD1D(JL)>=0 .AND. ZTIME_THRESHOLD1D(JL)<ZMAXTIME(JL) .AND. &
-                   &(ZVART(JL, JV)>ICED%XRTMIN(JV) .OR. ZA(JL, JV)>0.)) THEN
-                  ZMAXTIME(JL)=MIN(ZMAXTIME(JL), ZTIME_THRESHOLD1D(JL))
-                  LLCOMPUTE(JL)=.FALSE.
-                ENDIF
-              ENDDO
-            ENDDO
-            LLCPZ0RT=.FALSE.
-!$OMP SIMD
-            DO JL=1,IMICRO
-              ZMAXB(JL)=MAXVAL(ABS(ZB(JL,1:KRR)))
-            ENDDO
-            DO JL=1, IMICRO
-              IF (IITER(JL)<INB_ITER_MAX .AND. ZMAXB(JL)>PARAMI%XMRSTEP) THEN
-                ZMAXTIME(JL)=0.
-                LLCOMPUTE(JL)=.FALSE.
-              ENDIF
-            ENDDO
-          ENDIF ! LL_ANY_ITER
-        ENDIF ! XMRSTEP/=0.
-        !-------------------------------------------------------------------------------
-        !
-        !***       4.7 New values of variables for next iteration
-        !
-        !
-        DO JV=0, KRR
-          DO JL=1, IMICRO
-            ZVART(JL, JV)=ZVART(JL, JV)+ZA(JL, JV)*ZMAXTIME(JL)+ZB(JL, JV)
-          ENDDO
-        ENDDO
-        DO JL=1, IMICRO
-#ifdef REPRO55
-          ZCIT(JL)=ZCIT(JL) * MAX(0., -SIGN(1., -ZVART(JL,IRI)))
-#else
-          IF (ZVART(JL,IRI)<=0.) ZCIT(JL) = 0.
-#endif
-          ZTIME(JL)=ZTIME(JL)+ZMAXTIME(JL)
-        ENDDO
-        !-------------------------------------------------------------------------------
-        !
-        !***       4.8 Mixing ratio change due to each process
-        !
-        IF(BUCONF%LBU_ENABLE) THEN
-          !Mixing ratio change due to a tendency
-          DO JV=1, IBUNUM-IBUNUM_MR-IBUNUM_EXTRA
-            DO JL=1, IMICRO
-              ZBU_SUM(JL, JV) = ZBU_SUM(JL, JV) + ZBU_INST(JL, JV)*ZMAXTIME(JL)
-            ENDDO
-          ENDDO
-
-          !Mixing ratio change due to a mixing ratio change
-          DO JV=IBUNUM-IBUNUM_MR-IBUNUM_EXTRA+1, IBUNUM-IBUNUM_EXTRA
-            DO JL=1, IMICRO
-              ZBU_SUM(JL, JV) = ZBU_SUM(JL, JV) + ZBU_INST(JL, JV)
-            ENDDO
-          ENDDO
-
-          !Extra contribution as a mixing ratio change
-          DO JV=IBUNUM-IBUNUM_EXTRA+1, IBUNUM
-            JJV=IBUEXTRAIND(JV)
-            DO JL=1, IMICRO
-              ZBU_SUM(JL, JJV) = ZBU_SUM(JL, JJV) + ZBU_INST(JL, JV)
-            ENDDO
-          ENDDO
-        ENDIF
-        !-------------------------------------------------------------------------------
-        !
-        !***       4.9 Next loop
-        !
-        LSOFT=.TRUE. ! We try to adjust tendencies (inner while loop)
-      ENDDO !Iterations on tendency computations (WHILE ANY(LLCOMPUTE))
-    ENDDO !Temporal loop
-
-    IF(GEXT_TEND) THEN
-      !Z..T variables contain the external tendency, we substract it
-      DO JV=0, KRR
-        DO JL=1, IMICRO
-          ZVART(JL, JV) = ZVART(JL, JV) - ZEXTPK(JL, JV) * PTSTEP
-        ENDDO
-      ENDDO
-    ENDIF
-
-    !-------------------------------------------------------------------------------
-    !
-    !***       4.10 Unpacking
+    CALL ICE4_STEPPING(D, CST, PARAMI, ICEP, ICED, BUCONF, &
+                      &LLSIGMA_RC, LL_AUCV_ADJU, GEXT_TEND, &
+                      &KPROMA, IMICRO, PTSTEP, &
+                      &KRR, &
+                      &HSUBG_AUCV_RC, HSUBG_AUCV_RI, &
+                      &ZEXN, ZRHODREF, I1, I2, &
+                      &ZPRES, ZCF, ZSIGMA_RC, &
+                      &ZCIT, &
+                      &ZVART, &
+                      &ZHLC_HCF, ZHLC_HRC, &
+                      &ZHLI_HCF, ZHLI_HRI, PRAINFR, &
+                      &ZEXTPK, ZBU_SUM, ZRREVAV)
     !
+    !*       6.     UNPACKING
+    !               ---------
     !
     DO JL=1, IMICRO
       PCIT  (I1(JL),I2(JL))=ZCIT   (JL)
       IF(PARAMI%LWARM) THEN
-        PEVAP3D(I1(JL),I2(JL))=ZBU_INST(JL, IRREVAV)
+        PEVAP3D(I1(JL),I2(JL))=ZRREVAV(JL)
       ENDIF
       PWR(I1(JL),I2(JL),IRV)=ZVART(JL, IRV)
       PWR(I1(JL),I2(JL),IRC)=ZVART(JL, IRC)
@@ -618,10 +332,9 @@ IF (KSIZE > 0) THEN
 
   ENDDO ! JMICRO
 ENDIF ! KSIZE > 0
-!-------------------------------------------------------------------------------
-!
-!***       4.11 Budgets
 !
+!*       7.     BUDGETS
+!               -------
 !
 IF(BUCONF%LBU_ENABLE) THEN
   !Budgets for the different processes
diff --git a/src/common/micro/mode_ice4_stepping.F90 b/src/common/micro/mode_ice4_stepping.F90
new file mode 100644
index 000000000..fdd598b18
--- /dev/null
+++ b/src/common/micro/mode_ice4_stepping.F90
@@ -0,0 +1,421 @@
+!MNH_LIC Copyright 1995-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_STEPPING
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE ICE4_STEPPING(D, CST, PARAMI, ICEP, ICED, BUCONF, &
+                        &LDSIGMA_RC, LDAUCV_ADJU, LDEXT_TEND, &
+                        &KPROMA, KMICRO, PTSTEP, &
+                        &KRR, &
+                        &HSUBG_AUCV_RC, HSUBG_AUCV_RI, &
+                        &PEXN, PRHODREF, K1, K2, &
+                        &PPRES, PCF, PSIGMA_RC, &
+                        &PCIT, &
+                        &PVART, &
+                        &PHLC_HCF, PHLC_HRC, &
+                        &PHLI_HCF, PHLI_HRI, PRAINFR, &
+                        &PEXTPK, PBU_SUM, PRREVAV)
+
+!  -----------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+
+USE MODD_DIMPHYEX,       ONLY: DIMPHYEX_t
+USE MODD_BUDGET,         ONLY: TBUDGETCONF_t
+USE MODD_CST,            ONLY: CST_t
+USE MODD_PARAM_ICE,      ONLY: PARAM_ICE_t
+USE MODD_RAIN_ICE_DESCR, ONLY: RAIN_ICE_DESCR_t
+USE MODD_RAIN_ICE_PARAM, ONLY: RAIN_ICE_PARAM_t
+USE MODD_FIELDS_ADDRESS, ONLY : & ! common fields adress
+      & ITH,     & ! Potential temperature
+      & IRV,     & ! Water vapor
+      & IRC,     & ! Cloud water
+      & IRR,     & ! Rain water
+      & IRI,     & ! Pristine ice
+      & IBUNUM,       & ! Number of tendency terms
+      & IBUNUM_MR,    & ! Number of tendency terms expressed as mixing ratio changes
+      & IBUNUM_EXTRA, & ! Number of extra tendency terms
+      & IRREVAV,      & ! Index for the evaporation tendency
+      & IBUEXTRAIND     ! Index indirection
+
+USE MODE_MSG,            ONLY: PRINT_MSG, NVERB_FATAL
+
+USE MODE_ICE4_TENDENCIES, ONLY: ICE4_TENDENCIES
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of dummy arguments :
+!
+!
+!
+TYPE(DIMPHYEX_t),         INTENT(IN)    :: D
+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
+TYPE(TBUDGETCONF_t),      INTENT(IN)    :: BUCONF
+LOGICAL,                  INTENT(IN)    :: LDSIGMA_RC
+LOGICAL,                  INTENT(IN)    :: LDAUCV_ADJU
+LOGICAL,                  INTENT(IN)    :: LDEXT_TEND
+INTEGER,                  INTENT(IN)    :: KPROMA ! cache-blocking factor for microphysic loop
+INTEGER,                  INTENT(IN)    :: KMICRO ! Case r_x>0 locations
+REAL,                     INTENT(IN)    :: PTSTEP  ! Double Time step (single if cold start)
+INTEGER,                  INTENT(IN)    :: KRR     ! Number of moist variable
+CHARACTER(LEN=4),         INTENT(IN)    :: HSUBG_AUCV_RC ! Kind of Subgrid autoconversion method
+CHARACTER(LEN=80),        INTENT(IN)    :: HSUBG_AUCV_RI ! Kind of Subgrid autoconversion method
+!
+REAL,    DIMENSION(KPROMA),                     INTENT(IN)    :: PEXN    ! Exner function
+REAL,    DIMENSION(KPROMA),                     INTENT(IN)    :: PRHODREF! Reference density
+INTEGER, DIMENSION(KPROMA),                     INTENT(IN)    :: K1,K2 ! Used to replace the COUNT and PACK intrinsics on variables
+REAL,    DIMENSION(KPROMA),                     INTENT(IN)    :: PPRES
+REAL,    DIMENSION(KPROMA),                     INTENT(IN)    :: PCF ! Cloud fraction
+REAL,    DIMENSION(KPROMA),                     INTENT(INOUT) :: PSIGMA_RC
+REAL,    DIMENSION(KPROMA),                     INTENT(INOUT) :: PCIT
+REAL,    DIMENSION(KPROMA,0:7),                 INTENT(INOUT) :: PVART !Packed variables
+REAL,    DIMENSION(KPROMA),                     INTENT(INOUT) :: PHLC_HRC
+REAL,    DIMENSION(KPROMA),                     INTENT(INOUT) :: PHLC_HCF
+REAL,    DIMENSION(KPROMA),                     INTENT(INOUT) :: PHLI_HRI
+REAL,    DIMENSION(KPROMA),                     INTENT(INOUT) :: PHLI_HCF
+REAL,    DIMENSION(KPROMA),                     INTENT(OUT)   :: PRAINFR
+REAL,    DIMENSION(KPROMA,0:7),                 INTENT(INOUT) :: PEXTPK !To take into acount external tendencies inside the splitting
+REAL,    DIMENSION(KPROMA, IBUNUM-IBUNUM_EXTRA),INTENT(OUT)   :: PBU_SUM
+REAL,    DIMENSION(KPROMA),                     INTENT(OUT)   :: PRREVAV
+!
+!
+!*       0.2   Declarations of local variables :
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+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 :: ZTSTEP ! length of sub-timestep in case of time splitting
+REAL :: ZINV_TSTEP ! Inverse ov PTSTEP
+REAL :: ZTIME_THRESHOLD ! Time to reach threshold
+!
+INTEGER :: JL, JV, JJV
+REAL, DIMENSION(KPROMA) :: &
+                        & ZTIME, & ! Current integration time (starts with 0 and ends with PTSTEP)
+                        & ZMAXTIME, & ! Time on which we can apply the current tendencies
+                        & ZTIME_LASTCALL, &     ! Integration time when last tendecies call has been done
+                        & ZSSI,     &
+                        & ZZT,      & ! Temperature
+                        & ZLSFACT,  & ! L_s/(Pi*C_ph)
+                        & ZLVFACT,  & ! L_v/(Pi*C_ph)
+                        & ZHLC_LCF, & ! HLCLOUDS : fraction of Low  Cloud Fraction in grid
+                                      !    note that PCF = PHLC_HCF + ZHLC_LCF
+                        & ZHLC_LRC, & ! HLCLOUDS : LWC that is Low  LWC in grid
+                                      !    note that ZRC = PHLC_HRC + ZHLC_LRC
+                        & ZHLI_LCF, &
+                        & ZHLI_LRI
+LOGICAL, DIMENSION(KPROMA) :: LLCOMPUTE ! .TRUE. or points where we must compute tendencies,
+!
+!Output packed tendencies (for budgets only)
+REAL, DIMENSION(KPROMA, IBUNUM) :: ZBU_INST
+!
+!For mixing-ratio-splitting
+LOGICAL :: LLCPZ0RT
+REAL :: ZTIME_THRESHOLD1D(KPROMA) ! Time to reach threshold
+REAL, DIMENSION(KPROMA, KRR) :: Z0RT ! Mixing-ratios at the beginig of the current loop
+!
+REAL, DIMENSION(KPROMA,0:7) :: ZA, ZB
+!
+REAL, DIMENSION(KPROMA, 8) :: ZRS_TEND, ZRG_TEND
+REAL, DIMENSION(KPROMA,10) :: ZRH_TEND
+
+INTEGER, DIMENSION(KPROMA) :: IITER    ! Number of iterations done (with real tendencies computation)
+!
+REAL, DIMENSION(KPROMA) :: ZSUM2, ZMAXB
+REAL :: ZDEVIDE, ZX
+!
+LOGICAL :: LL_ANY_ITER
+!
+!-------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('ICE4_STEPPING', 0, ZHOOK_HANDLE)
+!
+!*       1.     GENERALITIES
+!               ------------
+!
+ZINV_TSTEP=1./PTSTEP
+!
+IF(BUCONF%LBU_ENABLE) THEN
+  DO JV=1, IBUNUM-IBUNUM_EXTRA
+    PBU_SUM(:, JV)=0.
+  ENDDO
+ENDIF
+
+!Maximum number of iterations
+!We only count real iterations (those for which we *compute* tendencies)
+INB_ITER_MAX=PARAMI%NMAXITER
+IF(PARAMI%XTSTEP_TS/=0.)THEN
+  INB_ITER_MAX=MAX(1, INT(PTSTEP/PARAMI%XTSTEP_TS)) !At least the number of iterations needed for the time-splitting
+  ZTSTEP=PTSTEP/INB_ITER_MAX
+  INB_ITER_MAX=MAX(PARAMI%NMAXITER, INB_ITER_MAX) !For the case XMRSTEP/=0. at the same time
+ENDIF
+
+IF (LDEXT_TEND) THEN
+  DO JV=0, KRR
+    DO JL=1, KMICRO
+      PEXTPK(JL, JV)=PEXTPK(JL, JV)-PVART(JL, JV)*ZINV_TSTEP
+    ENDDO
+  ENDDO
+ENDIF
+IF (LDSIGMA_RC) THEN
+  DO JL=1, KMICRO
+    PSIGMA_RC(JL)=PSIGMA_RC(JL)*2.
+  ENDDO 
+ENDIF
+IF (LDAUCV_ADJU) THEN
+  DO JL=1, KMICRO
+    ZHLC_LRC(JL) = PVART(JL, IRC) - PHLC_HRC(JL)
+    ZHLI_LRI(JL) = PVART(JL, IRI) - PHLI_HRI(JL)
+    IF(PVART(JL, IRC)>0.) THEN
+      ZHLC_LCF(JL) = PCF(JL)- PHLC_HCF(JL)
+    ELSE
+      ZHLC_LCF(JL)=0.
+    ENDIF
+    IF(PVART(JL, IRI)>0.) THEN
+      ZHLI_LCF(JL) = PCF(JL)- PHLI_HCF(JL)
+    ELSE
+      ZHLI_LCF(JL)=0.
+    ENDIF
+  ENDDO
+ENDIF
+
+!-------------------------------------------------------------------------------
+!
+!***       4.4 temporal loop
+!
+!
+IITER(1:KMICRO)=0
+ZTIME(1:KMICRO)=0. ! Current integration time (all points may have a different integration time)
+
+DO WHILE(ANY(ZTIME(1:KMICRO)<PTSTEP)) ! Loop to *really* compute tendencies
+
+  IF(PARAMI%XTSTEP_TS/=0.) THEN
+    ! In this case we need to remember the time when tendencies were computed
+    ! because when time has evolved more than a limit, we must re-compute tendencies
+    ZTIME_LASTCALL(1:KMICRO)=ZTIME(1:KMICRO)
+  ENDIF
+  DO JL=1, KMICRO
+    IF (ZTIME(JL) < PTSTEP) THEN
+      LLCOMPUTE(JL)=.TRUE. ! Computation (.TRUE.) only for points for which integration time has not reached the timestep
+      IITER(JL)=IITER(JL)+1
+    ELSE
+      LLCOMPUTE(JL)=.FALSE.
+    ENDIF
+  ENDDO
+  LL_ANY_ITER=ANY(IITER(1:KMICRO) < INB_ITER_MAX)
+  LLCPZ0RT=.TRUE.
+  LSOFT=.FALSE. ! We *really* compute the tendencies
+
+  DO WHILE(ANY(LLCOMPUTE(1:KMICRO))) ! Loop to adjust tendencies when we cross the 0°C or when a species disappears
+!$OMP SIMD
+    DO JL=1, KMICRO
+      ZSUM2(JL)=SUM(PVART(JL,IRI:KRR))
+    ENDDO
+    DO JL=1, KMICRO
+      ZDEVIDE=(CST%XCPD + CST%XCPV*PVART(JL, IRV) + CST%XCL*(PVART(JL, IRC)+PVART(JL, IRR)) + CST%XCI*ZSUM2(JL)) * PEXN(JL)
+      ZZT(JL) = PVART(JL, ITH) * PEXN(JL)
+      ZLSFACT(JL)=(CST%XLSTT+(CST%XCPV-CST%XCI)*(ZZT(JL)-CST%XTT)) / ZDEVIDE
+      ZLVFACT(JL)=(CST%XLVTT+(CST%XCPV-CST%XCL)*(ZZT(JL)-CST%XTT)) / ZDEVIDE
+    ENDDO
+    !-------------------------------------------------------------------------------
+    !
+    !***       4.5 Effective tendencies computation
+    !
+    !
+    ! Tendencies are *really* computed when LSOFT==.FALSE. and only adjusted otherwise
+    CALL ICE4_TENDENCIES(D, CST, PARAMI, ICEP, ICED, BUCONF, &
+                        &KPROMA, KMICRO, &
+                        &KRR, LSOFT, LLCOMPUTE, &
+                        &HSUBG_AUCV_RC, HSUBG_AUCV_RI, &
+                        &PEXN, PRHODREF, ZLVFACT, ZLSFACT, K1, K2, &
+                        &PPRES, PCF, PSIGMA_RC, &
+                        &PCIT, &
+                        &ZZT, PVART, &
+                        &ZBU_INST, &
+                        &ZRS_TEND, ZRG_TEND, ZRH_TEND, ZSSI, &
+                        &ZA, ZB, &
+                        &PHLC_HCF, ZHLC_LCF, PHLC_HRC, ZHLC_LRC, &
+                        &PHLI_HCF, ZHLI_LCF, PHLI_HRI, ZHLI_LRI, PRAINFR)
+
+    ! External tendencies
+    IF(LDEXT_TEND) THEN
+      DO JV=0, KRR
+        DO JL=1, KMICRO
+          ZA(JL, JV) = ZA(JL, JV) + PEXTPK(JL, JV)
+        ENDDO
+      ENDDO
+    ENDIF
+    !-------------------------------------------------------------------------------
+    !
+    !***       4.6 Time integration
+    !
+    !
+    ! If we can, we shall use these tendencies until the end of the timestep
+    DO JL=1, KMICRO
+      IF(LLCOMPUTE(JL)) THEN
+        ZMAXTIME(JL)=(PTSTEP-ZTIME(JL)) ! Remaining time until the end of the timestep
+      ELSE
+        ZMAXTIME(JL)=0.
+      ENDIF
+    ENDDO
+
+    !We need to adjust tendencies when temperature reaches 0
+    IF(PARAMI%LFEEDBACKT) THEN
+      DO JL=1, KMICRO
+        !Is ZB(:, ITH) enough to change temperature sign?
+        ZX=CST%XTT/PEXN(JL)
+        IF ((PVART(JL, ITH) - ZX) * (PVART(JL, ITH) + ZB(JL, ITH) - ZX) < 0.) THEN
+          ZMAXTIME(JL)=0.
+        ENDIF
+        !Can ZA(:, ITH) make temperature change of sign?
+        IF (ABS(ZA(JL,ITH)) > 1.E-20 ) THEN
+          ZTIME_THRESHOLD=(ZX - ZB(JL, ITH) - PVART(JL, ITH))/ZA(JL, ITH)
+          IF (ZTIME_THRESHOLD > 0.) THEN
+            ZMAXTIME(JL)=MIN(ZMAXTIME(JL), ZTIME_THRESHOLD)
+          ENDIF
+        ENDIF
+      ENDDO
+    ENDIF
+
+    !We need to adjust tendencies when a species disappears
+    !When a species is missing, only the external tendencies can be negative (and we must keep track of it)
+    DO JV=1, KRR
+      DO JL=1, KMICRO
+        IF (ZA(JL, JV) < -1.E-20 .AND. PVART(JL, JV) > ICED%XRTMIN(JV)) THEN
+          ZMAXTIME(JL)=MIN(ZMAXTIME(JL), -(ZB(JL, JV)+PVART(JL, JV))/ZA(JL, JV))
+        ENDIF
+      ENDDO
+    ENDDO
+
+    !We stop when the end of the timestep is reached
+    DO JL=1, KMICRO
+      IF (ZTIME(JL)+ZMAXTIME(JL) >= PTSTEP) THEN
+        LLCOMPUTE(JL)=.FALSE.
+      ENDIF
+    ENDDO
+    !We must recompute tendencies when the end of the sub-timestep is reached
+    IF (PARAMI%XTSTEP_TS/=0.) THEN
+      DO JL=1, KMICRO
+        IF ((IITER(JL) < INB_ITER_MAX) .AND. (ZTIME(JL)+ZMAXTIME(JL) > ZTIME_LASTCALL(JL)+ZTSTEP)) THEN
+          ZMAXTIME(JL)=ZTIME_LASTCALL(JL)-ZTIME(JL)+ZTSTEP
+          LLCOMPUTE(JL)=.FALSE.
+        ENDIF
+      ENDDO
+    ENDIF
+
+    !We must recompute tendencies when the maximum allowed change is reached
+    !When a species is missing, only the external tendencies can be active and we do not want to recompute
+    !the microphysical tendencies when external tendencies are negative (results won't change because species was already missing)
+    IF (PARAMI%XMRSTEP/=0.) THEN
+      IF (LL_ANY_ITER) THEN
+        ! In this case we need to remember the initial mixing ratios used to compute the tendencies
+        ! because when mixing ratio has evolved more than a threshold, we must re-compute tendencies
+        ! Thus, at first iteration (ie when LLCPZ0RT=.TRUE.) we copy PVART into Z0RT
+        DO JV=1,KRR
+          IF (LLCPZ0RT) Z0RT(1:KMICRO, JV)=PVART(1:KMICRO, JV)
+          DO JL=1, KMICRO
+            IF (IITER(JL)<INB_ITER_MAX .AND. ABS(ZA(JL,JV))>1.E-20) THEN
+              ZTIME_THRESHOLD1D(JL)=(SIGN(1., ZA(JL, JV))*PARAMI%XMRSTEP+ &
+                                    &Z0RT(JL, JV)-PVART(JL, JV)-ZB(JL, JV))/ZA(JL, JV)
+            ELSE
+              ZTIME_THRESHOLD1D(JL)=-1.
+            ENDIF
+          ENDDO
+          DO JL=1, KMICRO
+            IF (ZTIME_THRESHOLD1D(JL)>=0 .AND. ZTIME_THRESHOLD1D(JL)<ZMAXTIME(JL) .AND. &
+               &(PVART(JL, JV)>ICED%XRTMIN(JV) .OR. ZA(JL, JV)>0.)) THEN
+              ZMAXTIME(JL)=MIN(ZMAXTIME(JL), ZTIME_THRESHOLD1D(JL))
+              LLCOMPUTE(JL)=.FALSE.
+            ENDIF
+          ENDDO
+        ENDDO
+        LLCPZ0RT=.FALSE.
+!$OMP SIMD
+        DO JL=1,KMICRO
+          ZMAXB(JL)=MAXVAL(ABS(ZB(JL,1:KRR)))
+        ENDDO
+        DO JL=1, KMICRO
+          IF (IITER(JL)<INB_ITER_MAX .AND. ZMAXB(JL)>PARAMI%XMRSTEP) THEN
+            ZMAXTIME(JL)=0.
+            LLCOMPUTE(JL)=.FALSE.
+          ENDIF
+        ENDDO
+      ENDIF ! LL_ANY_ITER
+    ENDIF ! XMRSTEP/=0.
+    !-------------------------------------------------------------------------------
+    !
+    !***       4.7 New values of variables for next iteration
+    !
+    !
+    DO JV=0, KRR
+      DO JL=1, KMICRO
+        PVART(JL, JV)=PVART(JL, JV)+ZA(JL, JV)*ZMAXTIME(JL)+ZB(JL, JV)
+      ENDDO
+    ENDDO
+    DO JL=1, KMICRO
+#ifdef REPRO55
+      PCIT(JL)=PCIT(JL) * MAX(0., -SIGN(1., -PVART(JL,IRI)))
+#else
+      IF (PVART(JL,IRI)<=0.) PCIT(JL) = 0.
+#endif
+      ZTIME(JL)=ZTIME(JL)+ZMAXTIME(JL)
+    ENDDO
+    !-------------------------------------------------------------------------------
+    !
+    !***       4.8 Mixing ratio change due to each process
+    !
+    IF(BUCONF%LBU_ENABLE) THEN
+      !Mixing ratio change due to a tendency
+      DO JV=1, IBUNUM-IBUNUM_MR-IBUNUM_EXTRA
+        DO JL=1, KMICRO
+          PBU_SUM(JL, JV) = PBU_SUM(JL, JV) + ZBU_INST(JL, JV)*ZMAXTIME(JL)
+        ENDDO
+      ENDDO
+
+      !Mixing ratio change due to a mixing ratio change
+      DO JV=IBUNUM-IBUNUM_MR-IBUNUM_EXTRA+1, IBUNUM-IBUNUM_EXTRA
+        DO JL=1, KMICRO
+          PBU_SUM(JL, JV) = PBU_SUM(JL, JV) + ZBU_INST(JL, JV)
+        ENDDO
+      ENDDO
+
+      !Extra contribution as a mixing ratio change
+      DO JV=IBUNUM-IBUNUM_EXTRA+1, IBUNUM
+        JJV=IBUEXTRAIND(JV)
+        DO JL=1, KMICRO
+          PBU_SUM(JL, JJV) = PBU_SUM(JL, JJV) + ZBU_INST(JL, JV)
+        ENDDO
+      ENDDO
+    ENDIF
+    !-------------------------------------------------------------------------------
+    !
+    !***       4.9 Next loop
+    !
+    LSOFT=.TRUE. ! We try to adjust tendencies (inner while loop)
+  ENDDO !Iterations on tendency computations (WHILE ANY(LLCOMPUTE))
+ENDDO !Temporal loop
+
+IF(LDEXT_TEND) THEN
+  !Z..T variables contain the external tendency, we substract it
+  DO JV=0, KRR
+    DO JL=1, KMICRO
+      PVART(JL, JV) = PVART(JL, JV) - PEXTPK(JL, JV) * PTSTEP
+    ENDDO
+  ENDDO
+ENDIF
+DO JL=1, KMICRO
+  PRREVAV(JL)=ZBU_INST(JL, IRREVAV)
+ENDDO
+!
+IF (LHOOK) CALL DR_HOOK('ICE4_STEPPING', 1, ZHOOK_HANDLE)
+END SUBROUTINE ICE4_STEPPING
+END MODULE MODE_ICE4_STEPPING
diff --git a/src/common/micro/mode_ice4_tendencies.F90 b/src/common/micro/mode_ice4_tendencies.F90
index 6be883b19..b85db3a98 100644
--- a/src/common/micro/mode_ice4_tendencies.F90
+++ b/src/common/micro/mode_ice4_tendencies.F90
@@ -95,7 +95,7 @@ REAL, DIMENSION(KPROMA, IBUNUM),INTENT(INOUT):: PBU_INST
 REAL, DIMENSION(KPROMA, 8),    INTENT(INOUT) :: PRS_TEND
 REAL, DIMENSION(KPROMA, 8),    INTENT(INOUT) :: PRG_TEND
 REAL, DIMENSION(KPROMA, 10),   INTENT(INOUT) :: PRH_TEND
-REAL, DIMENSION(KPROMA),       INTENT(OUT)   :: PSSI
+REAL, DIMENSION(KPROMA),       INTENT(INOUT) :: PSSI
 REAL, DIMENSION(KPROMA,0:7),   INTENT(OUT)   :: PA
 REAL, DIMENSION(KPROMA,0:7),   INTENT(OUT)   :: PB
 REAL, DIMENSION(KPROMA),       INTENT(INOUT)   :: PHLC_HCF
diff --git a/src/common/micro/rain_ice.F90 b/src/common/micro/rain_ice.F90
index 38e184690..904cdedbe 100644
--- a/src/common/micro/rain_ice.F90
+++ b/src/common/micro/rain_ice.F90
@@ -386,7 +386,7 @@ DO JK=IKTB,IKTE
       ZW3D(JIJ, JK)=ZZ_LSFACT(JIJ, JK)/PEXN(JIJ, JK)
 #ifdef REPRO55
 #else
-      PCIT(JIJ,JK)=0.
+      PCIT(JIJ,JK)=0. !ri=0 because where are in the not odmicro case
 #endif
  
     ELSE
-- 
GitLab