diff --git a/src/common/turb/mode_compute_entr_detr.F90 b/src/common/turb/mode_compute_entr_detr.F90
index 7d8fae529856f3a50b886d7d64389b5b8f30f7e2..77ef02bf56ab1911aa6a9a1dc9948e80473a9c9a 100644
--- a/src/common/turb/mode_compute_entr_detr.F90
+++ b/src/common/turb/mode_compute_entr_detr.F90
@@ -9,7 +9,8 @@
 IMPLICIT NONE
 CONTAINS
 !     ######spl
-          SUBROUTINE COMPUTE_ENTR_DETR(KK,KKB,KKE,KKL,OTEST,OTESTLCL,&
+          SUBROUTINE COMPUTE_ENTR_DETR(D, CST, NEB, PARAMMF,&
+                            KK,KKB,KKE,KKL,OTEST,OTESTLCL,&
                             HFRAC_ICE,PFRAC_ICE,PRHODREF,&
                             PPRE_MINUS_HALF,&
                             PPRE_PLUS_HALF,PZZ,PDZZ,&
@@ -72,10 +73,11 @@ CONTAINS
 !
 !*      0. DECLARATIONS
 !          ------------
-!                         
-USE MODD_CST
 !
-USE MODD_PARAM_MFSHALL_n
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
+USE MODD_CST,             ONLY: CST_t
+USE MODD_NEB,             ONLY: NEB_t
+USE MODD_PARAM_MFSHALL_n, ONLY: PARAM_MFSHALL_t
 !
 USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D 
 !
@@ -89,6 +91,11 @@ IMPLICIT NONE
 !*                    1.1  Declaration of Arguments
 !
 !
+TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
+TYPE(CST_t),            INTENT(IN)   :: CST
+TYPE(NEB_t),            INTENT(IN)   :: NEB
+TYPE(PARAM_MFSHALL_t),  INTENT(IN)   :: PARAMMF
+!
 INTEGER,                INTENT(IN)   :: KK
 INTEGER,                INTENT(IN)   :: KKB          ! near ground physical index
 INTEGER,                INTENT(IN)   :: KKE          ! uppest atmosphere physical index
@@ -173,9 +180,9 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !                ------------------
   IF (LHOOK) CALL DR_HOOK('COMPUTE_ENTR_DETR',0,ZHOOK_HANDLE)
   
-  ZRVORD   = XRV / XRD   !=1.607
-  ZG_O_THVREF(:)=XG/PTHVM(:,KK)
-  ZCOEFFMF_CLOUD=XENTR_MF * XG / XCRAD_MF
+  ZRVORD   = CST%XRV / CST%XRD   !=1.607
+  ZG_O_THVREF(:)=CST%XG/PTHVM(:,KK)
+  ZCOEFFMF_CLOUD=PARAMMF%XENTR_MF * CST%XG / PARAMMF%XCRAD_MF
   
   ZFRAC_ICE(:)=PFRAC_ICE(:) ! to not modify fraction of ice
  
@@ -190,15 +197,15 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
       ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)
     ELSE IF (OTEST(JLOOP) .AND. .NOT. OTESTLCL(JLOOP)) THEN
       !Temperature at flux level KK
-      ZT=PTH_UP(JLOOP)*(PPRE_MINUS_HALF(JLOOP)/XP00) ** (XRD/XCPD)
+      ZT=PTH_UP(JLOOP)*(PPRE_MINUS_HALF(JLOOP)/CST%XP00) ** (CST%XRD/CST%XCPD)
       !Saturating vapor pressure at flux level KK
-      ZFOESW = MIN(EXP( XALPW - XBETAW/ZT - XGAMW*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JLOOP))
-      ZFOESI = MIN(EXP( XALPI - XBETAI/ZT - XGAMI*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JLOOP))
+      ZFOESW = MIN(EXP( CST%XALPW - CST%XBETAW/ZT - CST%XGAMW*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JLOOP))
+      ZFOESI = MIN(EXP( CST%XALPI - CST%XBETAI/ZT - CST%XGAMI*LOG(ZT)  ), 0.99*PPRE_MINUS_HALF(JLOOP))
       !Computation of d.Rsat / dP (partial derivations with respect to P and T
       !and use of T=Theta*(P/P0)**(R/Cp) to transform dT into dP with theta_up
       !constant at the vertical)
-      ZDRSATODP=(XBETAW/ZT-XGAMW)*(1-ZFRAC_ICE(JLOOP))+(XBETAI/ZT-XGAMI)*ZFRAC_ICE(JLOOP)
-      ZDRSATODP=((XRD/XCPD)*ZDRSATODP-1.)*PRSAT_UP(JLOOP)/ &
+      ZDRSATODP=(CST%XBETAW/ZT-CST%XGAMW)*(1-ZFRAC_ICE(JLOOP))+(CST%XBETAI/ZT-CST%XGAMI)*ZFRAC_ICE(JLOOP)
+      ZDRSATODP=((CST%XRD/CST%XCPD)*ZDRSATODP-1.)*PRSAT_UP(JLOOP)/ &
                   &(PPRE_MINUS_HALF(JLOOP)-(ZFOESW*(1-ZFRAC_ICE(JLOOP)) + ZFOESI*ZFRAC_ICE(JLOOP)))
       !Use of d.Rsat / dP and pressure at flux level KK to find pressure (ZPRE)
       !where Rsat is equal to PRT_UP
@@ -246,22 +253,22 @@ DO JLOOP=1,SIZE(OTEST)
 
     !Entr//Detr. computation
     IF (PBUO_INTEG_DRY(JLOOP)>=0.) THEN
-      PENTR(JLOOP) = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
-                 LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(JLOOP,KK))* &
+      PENTR(JLOOP) = 0.5/(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)*&
+                 LOG(1.+ (2.*(PARAMMF%XABUO-PARAMMF%XBENTR*PARAMMF%XENTR_DRY)/PW_UP2(JLOOP,KK))* &
                  PBUO_INTEG_DRY(JLOOP))
       PDETR(JLOOP) = 0.
     ELSE
       PENTR(JLOOP) = 0.
-      PDETR(JLOOP) = 0.5/(XABUO)*&
-                 LOG(1.+ (2.*(XABUO)/PW_UP2(JLOOP,KK))* &
+      PDETR(JLOOP) = 0.5/(PARAMMF%XABUO)*&
+                 LOG(1.+ (2.*(PARAMMF%XABUO)/PW_UP2(JLOOP,KK))* &
                  (-PBUO_INTEG_DRY(JLOOP)))
     ENDIF
-    PENTR(JLOOP) = XENTR_DRY*PENTR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))    
-    PDETR(JLOOP) = XDETR_DRY*PDETR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
+    PENTR(JLOOP) = PARAMMF%XENTR_DRY*PENTR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))    
+    PDETR(JLOOP) = PARAMMF%XDETR_DRY*PDETR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
     !Minimum value of detrainment
     ZWK=PLUP(JLOOP)-0.5*(PZZ(JLOOP,KK)+PZZ(JLOOP,KK+KKL))
     ZWK=SIGN(MAX(1., ABS(ZWK)), ZWK) ! ZWK must not be zero
-    PDETR(JLOOP) = MAX(PPART_DRY(JLOOP)*XDETR_LUP/ZWK, PDETR(JLOOP))
+    PDETR(JLOOP) = MAX(PPART_DRY(JLOOP)*PARAMMF%XDETR_LUP/ZWK, PDETR(JLOOP))
   ELSE
     !No dry part, condensation reached (OTESTLCL)
     PBUO_INTEG_DRY(JLOOP) = 0.
@@ -280,7 +287,7 @@ ENDDO
   !but we are dealing with updraft and not mixture
   ZRCMIX(:)=PRC_UP(:)
   ZRIMIX(:)=PRI_UP(:)
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
+  CALL TH_R_FROM_THL_RT_1D(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
                PPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
                ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
                ZRSATW, ZRSATI,OOCEAN=.FALSE.)
@@ -357,7 +364,7 @@ DO JLOOP=1,SIZE(OTEST)
     ZMIXRT(JLOOP) = 0.1
   ENDIF
 ENDDO
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
+  CALL TH_R_FROM_THL_RT_1D(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
                ZPRE,ZMIXTHL,ZMIXRT,&
                ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
                ZRSATW, ZRSATI,OOCEAN=.FALSE.)
@@ -366,7 +373,7 @@ ENDDO
   !  Compute cons then non cons. var. of mixture at the flux level KK+KKL  with initial ZKIC
   ZMIXTHL(:) = ZKIC_INIT * 0.5*(PTHLM(:,KK)+PTHLM(:,KK+KKL))+(1. - ZKIC_INIT)*PTHL_UP(:)
   ZMIXRT(:)  = ZKIC_INIT * 0.5*(PRTM(:,KK)+PRTM(:,KK+KKL))+(1. - ZKIC_INIT)*PRT_UP(:)
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
+  CALL TH_R_FROM_THL_RT_1D(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
                PPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
                ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
                ZRSATW, ZRSATI,OOCEAN=.FALSE.)
diff --git a/src/common/turb/mode_compute_updraft.F90 b/src/common/turb/mode_compute_updraft.F90
index 134b7471520c6e8d6c3a1d8160d5a374a82fc9be..f53e782e4c65043416f92c5da7db658b11df2aa5 100644
--- a/src/common/turb/mode_compute_updraft.F90
+++ b/src/common/turb/mode_compute_updraft.F90
@@ -296,7 +296,7 @@ IF (OENTR_DETR) THEN
   ! (all or nothing ajustement)
   PRC_UP(:,D%NKB)=0.
   PRI_UP(:,D%NKB)=0.
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
+  CALL TH_R_FROM_THL_RT_1D(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
              PTHL_UP(:,D%NKB),PRT_UP(:,D%NKB),ZTH_UP(:,D%NKB), &
              PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
 
@@ -391,7 +391,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
       ZRC_MIX(:,JK) = ZRC_MIX(:,JK-D%NKL) ! guess of Rc of mixture
       ZRI_MIX(:,JK) = ZRI_MIX(:,JK-D%NKL) ! guess of Ri of mixture
     ENDIF
-    CALL COMPUTE_ENTR_DETR(JK,D%NKB,D%NKE,D%NKL,GTEST,GTESTLCL,HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
+    CALL COMPUTE_ENTR_DETR(D, CST, NEB, PARAMMF, JK,D%NKB,D%NKE,D%NKL,GTEST,GTESTLCL,HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
                            PRHODREF(:,JK),ZPRES_F(:,JK),ZPRES_F(:,JK+D%NKL),&
                            PZZ(:,:),PDZZ(:,:),ZTHVM(:,:),  &
                            PTHLM(:,:),PRTM(:,:),ZW_UP2(:,:),ZTH_UP(:,JK),   &
@@ -488,7 +488,7 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
 ! Compute non cons. var. at level JK+KKL
   ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
   ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
+  CALL TH_R_FROM_THL_RT_1D(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
           PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL),              &
           ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:), OOCEAN=.FALSE.)
   WHERE(GTEST)
diff --git a/src/common/turb/mode_compute_updraft_raha.F90 b/src/common/turb/mode_compute_updraft_raha.F90
index 41414ea157c7fd1d5772da28eca198bb52ce4991..5aedafcbcd151edc3b8261d656539366d4a4fe2f 100644
--- a/src/common/turb/mode_compute_updraft_raha.F90
+++ b/src/common/turb/mode_compute_updraft_raha.F90
@@ -60,6 +60,7 @@ CONTAINS
 !          ------------
 
 USE MODD_CST
+USE MODD_NEB, ONLY: NEB
 USE MODD_PARAM_MFSHALL_n
 
 USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
@@ -301,7 +302,7 @@ GTEST = (ZW_UP2(:,KKB) > ZEPS)
 PRC_UP(:,KKB)=0.
 PRI_UP(:,KKB)=0.
 
-CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,KKB),ZPRES_F(:,KKB), &
+CALL TH_R_FROM_THL_RT_1D(CST,NEB, SIZE(PFRAC_ICE_UP,1), HFRAC_ICE,PFRAC_ICE_UP(:,KKB),ZPRES_F(:,KKB), &
              PTHL_UP(:,KKB),PRT_UP(:,KKB),ZTH_UP(:,KKB), &
              PRV_UP(:,KKB),PRC_UP(:,KKB),PRI_UP(:,KKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
 
@@ -480,7 +481,7 @@ DO JK=KKB,KKE-KKL,KKL
   ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
   ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
   ZRV_UP(:)=PRV_UP(:,JK)
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,JK+KKL),ZPRES_F(:,JK+KKL), &
+  CALL TH_R_FROM_THL_RT_1D(CST,NEB, SIZE(PFRAC_ICE_UP,1), HFRAC_ICE,PFRAC_ICE_UP(:,JK+KKL),ZPRES_F(:,JK+KKL), &
           PTHL_UP(:,JK+KKL),PRT_UP(:,JK+KKL),ZTH_UP(:,JK+KKL),              &
           ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
   WHERE(GTEST)
diff --git a/src/common/turb/mode_compute_updraft_rhcj10.F90 b/src/common/turb/mode_compute_updraft_rhcj10.F90
index 0392012db244776f3720593fba148057f7a72a29..75bc2874f7a0f97048fd08799467aadd74069842 100644
--- a/src/common/turb/mode_compute_updraft_rhcj10.F90
+++ b/src/common/turb/mode_compute_updraft_rhcj10.F90
@@ -61,6 +61,7 @@ SUBROUTINE COMPUTE_UPDRAFT_RHCJ10(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE,       &
 !          ------------
 !
 USE MODD_CST
+USE MODD_NEB, ONLY: NEB
 USE MODD_PARAM_MFSHALL_n
 USE MODD_TURB_n, ONLY : CTURBLEN
 USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
@@ -308,7 +309,7 @@ ZW_UP2(:,KKB) = MAX(0.0001,(2./3.)*ZTKEM_F(:,KKB))
 
 PRC_UP(:,KKB)=0.
 PRI_UP(:,KKB)=0.
-CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,KKB),ZPRES_F(:,KKB), &
+CALL TH_R_FROM_THL_RT_1D(CST,NEB,SIZE(PFRAC_ICE_UP,1),HFRAC_ICE,PFRAC_ICE_UP(:,KKB),ZPRES_F(:,KKB), &
              PTHL_UP(:,KKB),PRT_UP(:,KKB),ZTH_UP(:,KKB), &
              PRV_UP(:,KKB),PRC_UP(:,KKB),PRI_UP(:,KKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
 
@@ -414,7 +415,7 @@ DO JK=KKB,KKE-KKL,KKL
     ZRC_UP(:)   =PRC_UP(:,JK) ! guess
     ZRI_UP(:)   =PRI_UP(:,JK) ! guess 
     ZRV_UP(:)   =PRV_UP(:,JK)
-    CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
+    CALL TH_R_FROM_THL_RT_1D(CST,NEB, SIZE(PFRAC_ICE_UP,1), HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
                PPABSM(:,JK),PTHL_UP(:,JK),PRT_UP(:,JK),&
                ZTH_UP(:,JK),ZRV_UP,ZRC_UP,ZRI_UP,ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)            
     
@@ -514,7 +515,7 @@ DO JK=KKB,KKE-KKL,KKL
   ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
   ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
   ZRV_UP(:)=PRV_UP(:,JK)
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,JK+KKL),ZPRES_F(:,JK+KKL), &
+  CALL TH_R_FROM_THL_RT_1D(CST,NEB, SIZE(PFRAC_ICE_UP,1), HFRAC_ICE,PFRAC_ICE_UP(:,JK+KKL),ZPRES_F(:,JK+KKL), &
           PTHL_UP(:,JK+KKL),PRT_UP(:,JK+KKL),ZTH_UP(:,JK+KKL),              &
           ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
 
diff --git a/src/common/turb/mode_th_r_from_thl_rt_1d.F90 b/src/common/turb/mode_th_r_from_thl_rt_1d.F90
index 1aff6f5943ff5a99cf55a46a60c48fdb72ca3118..1c974d8abab215de0b53bc4e82a2a465f781cc9f 100644
--- a/src/common/turb/mode_th_r_from_thl_rt_1d.F90
+++ b/src/common/turb/mode_th_r_from_thl_rt_1d.F90
@@ -5,8 +5,8 @@
 MODULE MODE_TH_R_FROM_THL_RT_1D
 IMPLICIT NONE
 CONTAINS
-      SUBROUTINE TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,PP,             &
-                                  PTHL, PRT, PTH, PRV, PRL, PRI,         &
+      SUBROUTINE TH_R_FROM_THL_RT_1D(CST, NEB, KT, HFRAC_ICE,PFRAC_ICE,PP, &
+                                  PTHL, PRT, PTH, PRV, PRL, PRI,           &
                                   PRSATW, PRSATI, PRR, PRS, PRG, PRH,OOCEAN)
 !     #################################################################
 !
@@ -50,8 +50,8 @@ CONTAINS
 !
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-USE MODD_CST !, ONLY: XP00, XRD, XCPD, XCPV, XCL, XCI, XLVTT, XTT, XLSTT
-USE MODD_NEB, ONLY: NEB
+USE MODD_CST, ONLY : CST_t
+USE MODD_NEB, ONLY : NEB_t
 USE MODE_THERMO
 !
 IMPLICIT NONE
@@ -59,19 +59,22 @@ IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
 !
+TYPE(CST_t),        INTENT(IN) :: CST
+TYPE(NEB_t),        INTENT(IN) :: NEB
+INTEGER,            INTENT(IN) :: KT
 CHARACTER(LEN=1),   INTENT(IN) :: HFRAC_ICE
 LOGICAL,            INTENT(IN)   ::  OOCEAN       ! switch for Ocean model version
-REAL, DIMENSION(:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:), INTENT(IN) :: PP          ! Pressure
-REAL, DIMENSION(:), INTENT(IN) :: PTHL    ! thetal to transform into th
-REAL, DIMENSION(:),INTENT(IN)  :: PRT    ! Total mixing ratios to transform into rv,rc and ri
-REAL, DIMENSION(:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:), INTENT(OUT):: PTH    ! th
-REAL, DIMENSION(:), INTENT(OUT):: PRV    ! vapor mixing ratio
-REAL, DIMENSION(:), INTENT(INOUT):: PRL    ! vapor mixing ratio
-REAL, DIMENSION(:), INTENT(INOUT):: PRI    ! vapor mixing ratio
-REAL, DIMENSION(:), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
+REAL, DIMENSION(KT), INTENT(INOUT) :: PFRAC_ICE
+REAL, DIMENSION(KT), INTENT(IN) :: PP          ! Pressure
+REAL, DIMENSION(KT), INTENT(IN) :: PTHL    ! thetal to transform into th
+REAL, DIMENSION(KT),INTENT(IN)  :: PRT    ! Total mixing ratios to transform into rv,rc and ri
+REAL, DIMENSION(KT),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
+REAL, DIMENSION(KT), INTENT(OUT):: PTH    ! th
+REAL, DIMENSION(KT), INTENT(OUT):: PRV    ! vapor mixing ratio
+REAL, DIMENSION(KT), INTENT(INOUT):: PRL    ! vapor mixing ratio
+REAL, DIMENSION(KT), INTENT(INOUT):: PRI    ! vapor mixing ratio
+REAL, DIMENSION(KT), INTENT(OUT)  :: PRSATW ! estimated mixing ration at saturation over water
+REAL, DIMENSION(KT), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturation over ice
 !
 !-------------------------------------------------------------------------------
 !
@@ -79,14 +82,14 @@ REAL, DIMENSION(:), INTENT(OUT)  :: PRSATI ! estimated mixing ration at saturati
 INTEGER                       :: II ! Loop control
 INTEGER                       :: JITER ! number of iterations
 INTEGER                       :: J
-REAL, DIMENSION(SIZE(PTHL,1))   :: ZEXN
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZRVSAT,ZCPH,ZRLTEMP,ZCPH2
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZT,ZLVOCPEXN,ZLSOCPEXN
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZDRSATODT,ZDRSATODTW,ZDRSATODTI
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZFOESW, ZFOESI
-REAL, DIMENSION(SIZE(PTHL,1)) :: ZLOGT, Z99PP, Z1PRT
-REAL(KIND=JPRB) :: ZVAR1, ZVAR2, ZTPOW2, ZDELT
-INTEGER, DIMENSION(SIZE(PTHL,1)) :: IERR
+REAL, DIMENSION(KT)   :: ZEXN
+REAL, DIMENSION(KT) :: ZRVSAT,ZCPH,ZRLTEMP,ZCPH2
+REAL, DIMENSION(KT) :: ZT,ZLVOCPEXN,ZLSOCPEXN
+REAL, DIMENSION(KT) :: ZDRSATODT,ZDRSATODTW,ZDRSATODTI
+REAL, DIMENSION(KT) :: ZFOESW, ZFOESI
+REAL, DIMENSION(KT) :: ZLOGT, Z99PP, Z1PRT
+REAL :: ZVAR1, ZVAR2, ZTPOW2, ZDELT
+INTEGER, DIMENSION(KT) :: IERR
 
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !----------------------------------------------------------------------------
@@ -102,24 +105,24 @@ JITER=2
 !
 !Computation of ZCPH2 depending on dummy arguments received
 ZCPH2(:)=0
-IF(PRESENT(PRR)) ZCPH2(:)=ZCPH2(:) + XCL*PRR(:)
-IF(PRESENT(PRS)) ZCPH2(:)=ZCPH2(:) + XCI*PRS(:)
-IF(PRESENT(PRG)) ZCPH2(:)=ZCPH2(:) + XCI*PRG(:)
-IF(PRESENT(PRH)) ZCPH2(:)=ZCPH2(:) + XCI*PRH(:)
+IF(PRESENT(PRR)) ZCPH2(:)=ZCPH2(:) + CST%XCL*PRR(:)
+IF(PRESENT(PRS)) ZCPH2(:)=ZCPH2(:) + CST%XCI*PRS(:)
+IF(PRESENT(PRG)) ZCPH2(:)=ZCPH2(:) + CST%XCI*PRG(:)
+IF(PRESENT(PRH)) ZCPH2(:)=ZCPH2(:) + CST%XCI*PRH(:)
 !
 !Computation of an approximate state thanks to PRL and PRI guess
-ZEXN(:)=(PP(:)/XP00) ** RDSCPD
+ZEXN(:)=(PP(:)/CST%XP00) ** CST%RDSCPD
 
-DO J=1,SIZE(PTHL,1)
-Z99PP(J)=0.99*PP(J)
-PRV(J)=PRT(J)-PRL(J)-PRI(J)
-ZCPH(J)=XCPD+ XCPV * PRV(J)+ XCL * PRL(J) + XCI * PRI(J) + ZCPH2(J)
-ZVAR2=ZCPH(J)*ZEXN(J)
-ZDELT=(PTHL(J)*ZEXN(J))-XTT
-ZLVOCPEXN(J) = (XLVTT + (XCPV-XCL) * ZDELT) /ZVAR2
-ZLSOCPEXN(J) = (XLSTT + (XCPV-XCI) * ZDELT) /ZVAR2 
-PTH(J)=PTHL(J)+ZLVOCPEXN(J)*PRL(J)+ZLSOCPEXN(J)*PRI(J)
-Z1PRT(J)=1+PRT(J)
+DO J=1,KT
+  Z99PP(J)=0.99*PP(J)
+  PRV(J)=PRT(J)-PRL(J)-PRI(J)
+  ZCPH(J)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + ZCPH2(J)
+  ZVAR2=ZCPH(J)*ZEXN(J)
+  ZDELT=(PTHL(J)*ZEXN(J))-CST%XTT
+  ZLVOCPEXN(J) = (CST%XLVTT + (CST%XCPV-CST%XCL) * ZDELT) /ZVAR2
+  ZLSOCPEXN(J) = (CST%XLSTT + (CST%XCPV-CST%XCI) * ZDELT) /ZVAR2 
+  PTH(J)=PTHL(J)+ZLVOCPEXN(J)*PRL(J)+ZLSOCPEXN(J)*PRI(J)
+  Z1PRT(J)=1+PRT(J)
 ENDDO
 !
 !
@@ -134,7 +137,7 @@ DO II=1,JITER
   END IF
   !Computation of liquid/ice fractions
   PFRAC_ICE(:) = 0.
-  DO J=1, SIZE(PFRAC_ICE, 1)
+  DO J=1, KT
     IF(PRL(J)+PRI(J) > 1.E-20) THEN
       PFRAC_ICE(J) = PRI(J) / (PRL(J)+PRI(J))
     ENDIF
@@ -148,55 +151,53 @@ DO II=1,JITER
   ! Log does not vectorize on all compilers:
   ZLOGT(:)=LOG(ZT(:))
 
-  DO J=1,SIZE(PTHL,1)
-
-  ZFOESW(J) = MIN(EXP( XALPW - XBETAW/ZT(J) - XGAMW*ZLOGT(J)  ), Z99PP(J))
-  ZFOESI(J) = MIN(EXP( XALPI - XBETAI/ZT(J) - XGAMI*ZLOGT(J)  ), Z99PP(J))
-  PRSATW(J) = XRD/XRV*ZFOESW(J)/PP(J) / (1.+(XRD/XRV-1.)*ZFOESW(J)/PP(J))
-  PRSATI(J) = XRD/XRV*ZFOESI(J)/PP(J) / (1.+(XRD/XRV-1.)*ZFOESI(J)/PP(J))
-  ZTPOW2=ZT(J)**2
-  ZDRSATODTW(J) = PRSATW(J) / (1.+(XRD/XRV-1.)*ZFOESW(J)/PP(J) ) &
-                   * (XBETAW/ZTPOW2 - XGAMW/ZT(J))*Z1PRT(J)
-  ZDRSATODTI(J) = PRSATI(J) / (1.+(XRD/XRV-1.)*ZFOESI(J)/PP(J) ) &
-                   * (XBETAI/ZTPOW2 - XGAMI/ZT(J))*Z1PRT(J)
-  !PRSATW(J) =  QSAT(ZT(J),PP(J)) !qsatw
-  !PRSATI(J) = QSATI(ZT(J),PP(J)) !qsati
-  !ZDRSATODTW(J) =  DQSAT(ZT(J),PP(J),PRSATW(J))*Z1PRT(J)
-  !ZDRSATODTI(J) = DQSATI(ZT(J),PP(J),PRSATI(J))*Z1PRT(J)
-  PRSATW(J) = PRSATW(J)*Z1PRT(J)
-  PRSATI(J) = PRSATI(J)*Z1PRT(J)
-  ZRVSAT(J) = PRSATW(J)*(1-PFRAC_ICE(J)) + PRSATI(J)*PFRAC_ICE(J)
-  ZDRSATODT(J) = (ZDRSATODTW(J)*(1-PFRAC_ICE(J))+ &
-            & ZDRSATODTI(J)*PFRAC_ICE(J))
+  DO J=1,KT
+    ZFOESW(J) = MIN(EXP( CST%XALPW - CST%XBETAW/ZT(J) - CST%XGAMW*ZLOGT(J)  ), Z99PP(J))
+    ZFOESI(J) = MIN(EXP( CST%XALPI - CST%XBETAI/ZT(J) - CST%XGAMI*ZLOGT(J)  ), Z99PP(J))
+    PRSATW(J) = CST%XRD/CST%XRV*ZFOESW(J)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESW(J)/PP(J))
+    PRSATI(J) = CST%XRD/CST%XRV*ZFOESI(J)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESI(J)/PP(J))
+    ZTPOW2=ZT(J)**2
+    ZDRSATODTW(J) = PRSATW(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESW(J)/PP(J) ) &
+                     * (CST%XBETAW/ZTPOW2 - CST%XGAMW/ZT(J))*Z1PRT(J)
+    ZDRSATODTI(J) = PRSATI(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESI(J)/PP(J) ) &
+                     * (CST%XBETAI/ZTPOW2 - CST%XGAMI/ZT(J))*Z1PRT(J)
+    !PRSATW(J) =  QSAT(ZT(J),PP(J)) !qsatw
+    !PRSATI(J) = QSATI(ZT(J),PP(J)) !qsati
+    !ZDRSATODTW(J) =  DQSAT(ZT(J),PP(J),PRSATW(J))*Z1PRT(J)
+    !ZDRSATODTI(J) = DQSATI(ZT(J),PP(J),PRSATI(J))*Z1PRT(J)
+    PRSATW(J) = PRSATW(J)*Z1PRT(J)
+    PRSATI(J) = PRSATI(J)*Z1PRT(J)
+    ZRVSAT(J) = PRSATW(J)*(1-PFRAC_ICE(J)) + PRSATI(J)*PFRAC_ICE(J)
+    ZDRSATODT(J) = (ZDRSATODTW(J)*(1-PFRAC_ICE(J))+ &
+              & ZDRSATODTI(J)*PFRAC_ICE(J))
 
-  !Computation of new PRL, PRI and PRV
-  !Correction term applied to (PRV(J)-ZRVSAT(J)) is computed assuming that
-  !ZLVOCPEXN, ZLSOCPEXN and ZCPH don't vary to much with T. It takes into account
-  !the variation (estimated linear) of Qsat with T
-  ZRLTEMP(J)=(PRV(J)-ZRVSAT(J))/ &
-                &(1 + ZDRSATODT(J)*ZEXN(J)* &
-                &     (ZLVOCPEXN(J)*(1-PFRAC_ICE(J))+ZLSOCPEXN(J)*PFRAC_ICE(J)))
-  ZRLTEMP(J)=MIN(MAX(-PRL(J)-PRI(J), ZRLTEMP(J)),PRV(J))
-  PRV(J)=PRV(J)-ZRLTEMP(J)
-  PRL(J)=PRL(J)+PRI(J)+ZRLTEMP(J)
-  PRI(J)=PFRAC_ICE(J)     * (PRL(J))
-  PRL(J)=(1-PFRAC_ICE(J)) * (PRT(J) - PRV(J))
+    !Computation of new PRL, PRI and PRV
+    !Correction term applied to (PRV(J)-ZRVSAT(J)) is computed assuming that
+    !ZLVOCPEXN, ZLSOCPEXN and ZCPH don't vary to much with T. It takes into account
+    !the variation (estimated linear) of Qsat with T
+    ZRLTEMP(J)=(PRV(J)-ZRVSAT(J))/ &
+                  &(1 + ZDRSATODT(J)*ZEXN(J)* &
+                  &     (ZLVOCPEXN(J)*(1-PFRAC_ICE(J))+ZLSOCPEXN(J)*PFRAC_ICE(J)))
+    ZRLTEMP(J)=MIN(MAX(-PRL(J)-PRI(J), ZRLTEMP(J)),PRV(J))
+    PRV(J)=PRV(J)-ZRLTEMP(J)
+    PRL(J)=PRL(J)+PRI(J)+ZRLTEMP(J)
+    PRI(J)=PFRAC_ICE(J)     * (PRL(J))
+    PRL(J)=(1-PFRAC_ICE(J)) * (PRT(J) - PRV(J))
 
-  !Computation of Cph (as defined in Meso-NH doc, equation 2.2, to be used with mixing ratios)
-  ZCPH(J)=XCPD+ XCPV * PRV(J)+ XCL * PRL(J) + XCI * PRI(J) + ZCPH2(J)
-
-  !Computation of L/Cph/EXN, then new PTH
-  ZVAR2=ZCPH(J)*ZEXN(J)
-  ZLVOCPEXN(J) = (XLVTT + (XCPV-XCL) * (ZT(J)-XTT)) /ZVAR2
-  ZLSOCPEXN(J) = (XLSTT + (XCPV-XCI) * (ZT(J)-XTT)) /ZVAR2
-  PTH(J)=PTHL(J)+ZLVOCPEXN(J)*PRL(J)+ZLSOCPEXN(J)*PRI(J)
+    !Computation of Cph (as defined in Meso-NH doc, equation 2.2, to be used with mixing ratios)
+    ZCPH(J)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + ZCPH2(J)
 
-  !Computation of estimated mixing ration at saturation
-  !To compute the adjustement a first order development was used
-  ZVAR1=PTH(J)*ZEXN(J)-ZT(J)
-  PRSATW(J)=PRSATW(J) + ZDRSATODTW(J)*ZVAR1
-  PRSATI(J)=PRSATI(J) + ZDRSATODTI(J)*ZVAR1
+    !Computation of L/Cph/EXN, then new PTH
+    ZVAR2=ZCPH(J)*ZEXN(J)
+    ZLVOCPEXN(J) = (CST%XLVTT + (CST%XCPV-CST%XCL) * (ZT(J)-CST%XTT)) /ZVAR2
+    ZLSOCPEXN(J) = (CST%XLSTT + (CST%XCPV-CST%XCI) * (ZT(J)-CST%XTT)) /ZVAR2
+    PTH(J)=PTHL(J)+ZLVOCPEXN(J)*PRL(J)+ZLSOCPEXN(J)*PRI(J)
 
+    !Computation of estimated mixing ration at saturation
+    !To compute the adjustement a first order development was used
+    ZVAR1=PTH(J)*ZEXN(J)-ZT(J)
+    PRSATW(J)=PRSATW(J) + ZDRSATODTW(J)*ZVAR1
+    PRSATI(J)=PRSATI(J) + ZDRSATODTI(J)*ZVAR1
   ENDDO
 ENDDO
 
diff --git a/src/common/turb/mode_th_r_from_thl_rt_2d.F90 b/src/common/turb/mode_th_r_from_thl_rt_2d.F90
index 2ac0c85284102e0ce45d4fa88a97f9fa474e708d..c4911018a8fba842c010932302d1dd557b3c46e9 100644
--- a/src/common/turb/mode_th_r_from_thl_rt_2d.F90
+++ b/src/common/turb/mode_th_r_from_thl_rt_2d.F90
@@ -49,6 +49,8 @@ CONTAINS
 USE MODE_TH_R_FROM_THL_RT_3D, ONLY: TH_R_FROM_THL_RT_3D
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+USE MODD_CST, ONLY: CST
+USE MODD_NEB, ONLY: NEB
 
 IMPLICIT NONE
 !
@@ -98,7 +100,7 @@ IF(PRESENT(PRH)) ZRH(:,:)=PRH(:,:)
 !       2 Call of 1d version
 !         ------------------
 DO JK=1, SIZE(PTHL,2)
-  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE(:,JK),PP(:,JK),             &
+  CALL TH_R_FROM_THL_RT_1D(CST, NEB, SIZE(PTHL,2), HFRAC_ICE,PFRAC_ICE(:,JK),PP(:,JK),             &
                                 PTHL(:,JK), PRT(:,JK), PTH(:,JK),       &
                                 PRV(:,JK), PRL(:,JK), PRI(:,JK),        &
                                 PRSATW(:,JK), PRSATI(:,JK),                &
diff --git a/src/common/turb/mode_th_r_from_thl_rt_3d.F90 b/src/common/turb/mode_th_r_from_thl_rt_3d.F90
index 1e179b139debd6d2c8e7c57146ba47244e7442ac..2c1733ce9dda37de77a20a615c06feb43f40999b 100644
--- a/src/common/turb/mode_th_r_from_thl_rt_3d.F90
+++ b/src/common/turb/mode_th_r_from_thl_rt_3d.F90
@@ -48,6 +48,8 @@ CONTAINS
 USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+USE MODD_CST, ONLY : CST
+USE MODD_NEB, ONLY : NEB
 !
 IMPLICIT NONE
 !
@@ -94,7 +96,7 @@ IF(PRESENT(PRH)) ZRH(:,:,:)=PRH(:,:,:)
 !         ------------------
 DO JK=1, SIZE(PTHL,3)
   DO JJ=1, SIZE(PTHL,2)
-    CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE(:,JJ,JK),PP(:,JJ,JK),             &
+    CALL TH_R_FROM_THL_RT_1D(CST, NEB, SIZE(PTHL,2), HFRAC_ICE,PFRAC_ICE(:,JJ,JK),PP(:,JJ,JK),             &
                                   PTHL(:,JJ,JK), PRT(:,JJ,JK), PTH(:,JJ,JK),       &
                                   PRV(:,JJ,JK), PRL(:,JJ,JK), PRI(:,JJ,JK),        &
                                   PRSATW(:,JJ,JK), PRSATI(:,JJ,JK),                &
diff --git a/src/mesonh/ext/set_rsou.f90 b/src/mesonh/ext/set_rsou.f90
index 9f7cca3e1c76de20267a8a398b309e005b804022..8943cfbcde9ada76f156fb26f2a44b1a176ad123 100644
--- a/src/mesonh/ext/set_rsou.f90
+++ b/src/mesonh/ext/set_rsou.f90
@@ -261,6 +261,7 @@ END MODULE MODI_SET_RSOU
 USE MODD_CONF
 USE MODD_CONF_n
 USE MODD_CST
+USE MODD_NEB, ONLY: NEB
 USE MODD_DYN_n,      ONLY: LOCEAN
 USE MODD_FIELD_n
 USE MODD_GRID
@@ -1592,7 +1593,7 @@ ELSE
   DO JLOOP=1,20 ! loop for pression 
     CALL COMPUTE_EXNER_FROM_GROUND(ZTHVM,ZZMASS_PROFILE(:),ZEXNSURF,ZEXNFLUX,ZEXNMASS)
     ZPRESS(:)=XP00*(ZEXNMASS(:))**(XCPD/XRD)
-    CALL TH_R_FROM_THL_RT_1D('T',ZFRAC_ICE,ZPRESS,ZTHLM,ZMRT,ZTHM,ZMRM,ZMRCM,ZMRIM, &
+    CALL TH_R_FROM_THL_RT_1D(CST,NEB,SIZE(ZPRESS,1),'T',ZFRAC_ICE,ZPRESS,ZTHLM,ZMRT,ZTHM,ZMRM,ZMRCM,ZMRIM, &
                               ZRSATW, ZRSATI,OOCEAN=.FALSE.)
      ZTHVM(:)=ZTHM(:)*(1.+XRV/XRD*ZMRM(:))/(1.+(ZMRM(:)+ZMRIM(:)+ZMRCM(:)))
   ENDDO