From a5ae0ed68b6105d2975e08883cfe28623d94a1c4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Riette?= <sebastien.riette@meteo.fr>
Date: Fri, 22 Apr 2022 10:48:55 +0200
Subject: [PATCH] S. Riette 22/4/2022 shuman_mf functions tranformed into
 subroutine

---
 src/common/turb/mode_compute_bl89_ml.F90      |   4 +-
 src/common/turb/mode_compute_entr_detr.F90    |   1 -
 .../turb/mode_compute_mf_cloud_bigaus.F90     |  17 +-
 .../turb/mode_compute_mf_cloud_stat.F90       |  17 +-
 src/common/turb/mode_compute_updraft.F90      |  29 +--
 src/common/turb/mode_compute_updraft_raha.F90 |  18 +-
 .../turb/mode_compute_updraft_rhcj10.F90      |  27 +--
 src/common/turb/mode_mf_turb.F90              |  33 ++--
 src/common/turb/mode_mf_turb_expl.F90         |  19 +-
 src/common/turb/mode_tridiag_massflux.F90     |   4 +-
 src/common/turb/shuman_mf.F90                 | 173 +++++++++---------
 11 files changed, 187 insertions(+), 155 deletions(-)

diff --git a/src/common/turb/mode_compute_bl89_ml.F90 b/src/common/turb/mode_compute_bl89_ml.F90
index 231af2f26..8a8449608 100644
--- a/src/common/turb/mode_compute_bl89_ml.F90
+++ b/src/common/turb/mode_compute_bl89_ml.F90
@@ -87,13 +87,13 @@ REAL    :: ZTEST,ZTEST0,ZTESTM  !test for vectorization
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('COMPUTE_BL89_ML',0,ZHOOK_HANDLE)
 !
-ZDELTVPT(:,:)=DZM_MF(PVPT(:,:), D%NKA, D%NKU, D%NKL)
+CALL DZM_MF(D, PVPT(:,:), ZDELTVPT(:,:))
 ZDELTVPT(:,D%NKA)=0.
 WHERE (ABS(ZDELTVPT(:,:))<CSTURB%XLINF)
   ZDELTVPT(:,:)=CSTURB%XLINF
 END WHERE
 !
-ZHLVPT(:,:)=MZM_MF(PVPT(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PVPT(:,:), ZHLVPT(:,:))
 !
 !We consider that gradient between mass levels KKB and KKB+KKL is the same as
 !the gradient between flux level KKB and mass level KKB
diff --git a/src/common/turb/mode_compute_entr_detr.F90 b/src/common/turb/mode_compute_entr_detr.F90
index 9f3d89ccb..cdcad823a 100644
--- a/src/common/turb/mode_compute_entr_detr.F90
+++ b/src/common/turb/mode_compute_entr_detr.F90
@@ -81,7 +81,6 @@ USE MODD_PARAM_MFSHALL_n, ONLY: PARAM_MFSHALL_t
 !
 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/common/turb/mode_compute_mf_cloud_bigaus.F90 b/src/common/turb/mode_compute_mf_cloud_bigaus.F90
index 0e1e0aa08..bfe4f3a91 100644
--- a/src/common/turb/mode_compute_mf_cloud_bigaus.F90
+++ b/src/common/turb/mode_compute_mf_cloud_bigaus.F90
@@ -62,9 +62,6 @@ USE MODD_CST,             ONLY: CST_t
 USE MODD_PARAM_MFSHALL_n, ONLY: PARAM_MFSHALL_t
 !
 USE MODI_SHUMAN_MF, ONLY: MZF_MF, GZ_M_W_MF
-!USE MODI_GAMMA_INC
-!
-!USE MODE_THERMO
 !
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -112,15 +109,15 @@ IF (LHOOK) CALL DR_HOOK('COMPUTE_MF_CLOUD_BIGAUS',0,ZHOOK_HANDLE)
 !
 !
 !Vertical gradient of RT, result on mass points
-ZW1(:,:)=GZ_M_W_MF(PRTM(:,:), PDZZ(:,:), D%NKA, D%NKU, D%NKL)
-ZGRAD_Z_RT(:,:)=MZF_MF(ZW1(:,:), D%NKA, D%NKU, D%NKL)
+CALL GZ_M_W_MF(D, PRTM(:,:), PDZZ(:,:), ZW1(:,:))
+CALL MZF_MF(D, ZW1(:,:), ZGRAD_Z_RT(:,:))
 
 !Interpolation on mass points
-ZTHV_UP_M(:,:) = MZF_MF(PTHV_UP(:,:), D%NKA, D%NKU, D%NKL)
-ZRSAT_UP_M(:,:)= MZF_MF(PRSAT_UP(:,:), D%NKA, D%NKU, D%NKL)
-ZRT_UP_M(:,:)  = MZF_MF(PRT_UP(:,:), D%NKA, D%NKU, D%NKL)
-ZEMF_M(:,:)    = MZF_MF(PEMF(:,:), D%NKA, D%NKU, D%NKL)
-ZFRAC_ICE_UP_M(:,:) = MZF_MF(PFRAC_ICE_UP(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZF_MF(D, PTHV_UP(:,:), ZTHV_UP_M(:,:))
+CALL MZF_MF(D, PRSAT_UP(:,:), ZRSAT_UP_M(:,:))
+CALL MZF_MF(D, PRT_UP(:,:), ZRT_UP_M(:,:))
+CALL MZF_MF(D, PEMF(:,:), ZEMF_M(:,:))
+CALL MZF_MF(D, PFRAC_ICE_UP(:,:), ZFRAC_ICE_UP_M(:,:))
 
 !computation of omega star up
 ZOMEGA_UP_M(:)=0.
diff --git a/src/common/turb/mode_compute_mf_cloud_stat.F90 b/src/common/turb/mode_compute_mf_cloud_stat.F90
index 5a75cbd63..c15549060 100644
--- a/src/common/turb/mode_compute_mf_cloud_stat.F90
+++ b/src/common/turb/mode_compute_mf_cloud_stat.F90
@@ -89,6 +89,7 @@ REAL, DIMENSION(D%NIT,D%NKT),   INTENT(OUT)  :: PSIGMF                  ! SQRT(v
 REAL, DIMENSION(D%NIT,D%NKT) :: ZFLXZ
 REAL, DIMENSION(D%NIT,D%NKT) :: ZT
 REAL, DIMENSION(D%NIT,D%NKT) :: ZAMOIST, ZATHETA
+REAL, DIMENSION(D%NIT,D%NKT) :: ZWK
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
 !*                    0.2 initialisation
@@ -112,14 +113,16 @@ IF (KRRL > 0)  THEN
 !
 
 !
-    ZFLXZ(:,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(PTHLM(:,:), D%NKA, D%NKU, D%NKL)) * &
-                      GZ_M_W_MF(PTHLM(:,:),PDZZ(:,:), D%NKA, D%NKU, D%NKL)
+    CALL MZM_MF(D, PTHLM(:,:), ZFLXZ(:,:))
+    CALL GZ_M_W_MF(D, PTHLM(:,:), PDZZ(:,:), ZWK(:,:))
+    ZFLXZ(:,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(:,:)*(PTHL_UP(:,:)-ZFLXZ(:,:)) * ZWK(:,:)
 !
 !   Avoid negative values
     ZFLXZ(:,:) = MAX(0.,ZFLXZ(:,:))
 
 
-    PSIGMF(:,:) = MZF_MF(ZFLXZ(:,:), D%NKA, D%NKU, D%NKL) * ZATHETA(:,:)**2
+    CALL MZF_MF(D, ZFLXZ(:,:), PSIGMF(:,:))
+    PSIGMF(:,:) = PSIGMF(:,:) * ZATHETA(:,:)**2
 
 !
 !
@@ -128,14 +131,16 @@ IF (KRRL > 0)  THEN
 !
 !
 !
-    ZFLXZ(:,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(:,:)*(PRT_UP(:,:)-MZM_MF(PRTM(:,:), D%NKA, D%NKU, D%NKL)) * &
-                      GZ_M_W_MF(PRTM(:,:),PDZZ(:,:), D%NKA, D%NKU, D%NKL)
+    CALL MZM_MF(D, PRTM(:,:), ZFLXZ(:,:))
+    CALL GZ_M_W_MF(D, PRTM(:,:), PDZZ(:,:), ZWK(:,:))
+    ZFLXZ(:,:) = -2 * PARAMMF%XTAUSIGMF * PEMF(:,:)*(PRT_UP(:,:)-ZFLXZ(:,:)) * ZWK(:,:)
 !
 !   Avoid negative values
     ZFLXZ(:,:) = MAX(0.,ZFLXZ(:,:))
 !
 
-    PSIGMF(:,:) = PSIGMF(:,:) + ZAMOIST(:,:) **2 * MZF_MF(ZFLXZ(:,:), D%NKA, D%NKU, D%NKL)
+    CALL MZF_MF(D, ZFLXZ(:,:), ZWK(:,:))
+    PSIGMF(:,:) = PSIGMF(:,:) + ZAMOIST(:,:) **2 * ZWK(:,:)
 !
 !        1.3  Vertical part of Sigma_s
 !
diff --git a/src/common/turb/mode_compute_updraft.F90 b/src/common/turb/mode_compute_updraft.F90
index 36e5e153f..76c8df61c 100644
--- a/src/common/turb/mode_compute_updraft.F90
+++ b/src/common/turb/mode_compute_updraft.F90
@@ -192,6 +192,9 @@ REAL  :: ZTMAX,ZRMAX  ! control value
 
 REAL, DIMENSION(D%NIT) :: ZSURF
 REAL, DIMENSION(D%NIT,D%NKT) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
+!
+REAL, DIMENSION(D%NIT,D%NKT) :: ZWK
+!
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
 IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT',0,ZHOOK_HANDLE)
@@ -249,15 +252,15 @@ END IF
 
 ! Initialisation of environment variables at t-dt
 ! variables at flux level
-ZTHLM_F(:,:) = MZM_MF(PTHLM(:,:), D%NKA, D%NKU, D%NKL)
-ZRTM_F (:,:) = MZM_MF(PRTM(:,:), D%NKA, D%NKU, D%NKL)
-ZUM_F  (:,:) = MZM_MF(PUM(:,:), D%NKA, D%NKU, D%NKL)
-ZVM_F  (:,:) = MZM_MF(PVM(:,:), D%NKA, D%NKU, D%NKL)
-ZTKEM_F(:,:) = MZM_MF(PTKEM(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PTHLM(:,:), ZTHLM_F(:,:))
+CALL MZM_MF(D, PRTM(:,:), ZRTM_F (:,:))
+CALL MZM_MF(D, PUM(:,:), ZUM_F  (:,:))
+CALL MZM_MF(D, PVM(:,:), ZVM_F  (:,:))
+CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 
 DO JSV=1,KSV
   IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
-  ZSVM_F(:,:,JSV) = MZM_MF(PSVM(:,:,JSV), D%NKA, D%NKU, D%NKL)
+ CALL MZM_MF(D, PSVM(:,:,JSV), ZSVM_F(:,:,JSV))
 END DO
 !                     
 !          Initialisation of updraft characteristics 
@@ -276,10 +279,10 @@ PRT_UP(:,D%NKB) = ZRTM_F(:,D%NKB)+MAX(0.,MIN(ZRMAX,(PSFRV(:)/SQRT(ZTKEM_F(:,D%NK
 
 
 IF (OENTR_DETR) THEN
-  ZTHM_F (:,:) = MZM_MF(PTHM (:,:), D%NKA, D%NKU, D%NKL)
-  ZPRES_F(:,:) = MZM_MF(PPABSM(:,:), D%NKA, D%NKU, D%NKL)
-  ZRHO_F (:,:) = MZM_MF(PRHODREF(:,:), D%NKA, D%NKU, D%NKL)
-  ZRVM_F (:,:) = MZM_MF(PRVM(:,:), D%NKA, D%NKU, D%NKL)
+  CALL MZM_MF(D, PTHM (:,:), ZTHM_F (:,:))
+  CALL MZM_MF(D, PPABSM(:,:), ZPRES_F(:,:))
+  CALL MZM_MF(D, PRHODREF(:,:), ZRHO_F (:,:))
+  CALL MZM_MF(D, PRVM(:,:), ZRVM_F (:,:))
 
   ! thetav at mass and flux levels
   ZTHVM_F(:,:)=ZTHM_F(:,:)*((1.+ZRVORD*ZRVM_F(:,:))/(1.+ZRTM_F(:,:)))
@@ -314,8 +317,10 @@ IF (OENTR_DETR) THEN
   ZTKEM_F(:,D%NKB)=0.
   !
   IF(TURB%CTURBLEN=='RM17') THEN
-    ZDUDZ = MZF_MF(GZ_M_W_MF(PUM,PDZZ, D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL)
-    ZDVDZ = MZF_MF(GZ_M_W_MF(PVM,PDZZ, D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL)
+    CALL GZ_M_W_MF(D, PUM, PDZZ, ZWK)
+    CALL MZF_MF(D, ZWK, ZDUDZ)
+    CALL GZ_M_W_MF(D, PVM, PDZZ, ZWK)
+    CALL MZF_MF(D, ZWK, ZDVDZ)
     ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
   ELSE
     ZSHEAR = 0. !no shear in bl89 mixing length
diff --git a/src/common/turb/mode_compute_updraft_raha.F90 b/src/common/turb/mode_compute_updraft_raha.F90
index 3cbacbc0c..c2a5cadab 100644
--- a/src/common/turb/mode_compute_updraft_raha.F90
+++ b/src/common/turb/mode_compute_updraft_raha.F90
@@ -241,11 +241,11 @@ PRSAT_UP(:,:)=PRVM(:,:) ! should be initialised correctly but is (normaly) not u
 ! Initialisation of environment variables at t-dt
 
 ! variables at flux level
-ZTHLM_F(:,:) = MZM_MF(PTHLM(:,:), D%NKA, D%NKU, D%NKL)
-ZRTM_F (:,:) = MZM_MF(PRTM(:,:), D%NKA, D%NKU, D%NKL)
-ZUM_F  (:,:) = MZM_MF(PUM(:,:), D%NKA, D%NKU, D%NKL)
-ZVM_F  (:,:) = MZM_MF(PVM(:,:), D%NKA, D%NKU, D%NKL)
-ZTKEM_F(:,:) = MZM_MF(PTKEM(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PTHLM(:,:), ZTHLM_F(:,:))
+CALL MZM_MF(D, PRTM(:,:), ZRTM_F(:,:))
+CALL MZM_MF(D, PUM(:,:), ZUM_F(:,:))
+CALL MZM_MF(D, PVM(:,:), ZVM_F(:,:))
+CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 
 !DO JSV=1,ISV 
 ! IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
@@ -272,10 +272,10 @@ PRT_UP(:,D%NKB) = ZRTM_F(:,D%NKB)+MAX(0.,MIN(ZRMAX,(PSFRV(:)/SQRT(ZTKEM_F(:,D%NK
 ZQT_UP(:) = PRT_UP(:,D%NKB)/(1.+PRT_UP(:,D%NKB))
 ZTHS_UP(:,D%NKB)=PTHL_UP(:,D%NKB)*(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(:))
 
-ZTHM_F (:,:) = MZM_MF(PTHM (:,:), D%NKA, D%NKU, D%NKL)
-ZPRES_F(:,:) = MZM_MF(PPABSM(:,:), D%NKA, D%NKU, D%NKL)
-ZRHO_F (:,:) = MZM_MF(PRHODREF(:,:), D%NKA, D%NKU, D%NKL)
-ZRVM_F (:,:) = MZM_MF(PRVM(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PTHM (:,:), ZTHM_F(:,:))
+CALL MZM_MF(D, PPABSM(:,:), ZPRES_F(:,:))
+CALL MZM_MF(D, PRHODREF(:,:), ZRHO_F(:,:))
+CALL MZM_MF(D, PRVM(:,:), ZRVM_F(:,:))
 
 ! thetav at mass and flux levels 
 ZTHVM_F(:,:)=ZTHM_F(:,:)*((1.+ZRVORD*ZRVM_F(:,:))/(1.+ZRTM_F(:,:)))
diff --git a/src/common/turb/mode_compute_updraft_rhcj10.F90 b/src/common/turb/mode_compute_updraft_rhcj10.F90
index 6de2c738f..57b31c12d 100644
--- a/src/common/turb/mode_compute_updraft_rhcj10.F90
+++ b/src/common/turb/mode_compute_updraft_rhcj10.F90
@@ -190,6 +190,9 @@ REAL  :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
 REAL  :: ZTMAX,ZRMAX, ZEPS  ! control value
 
 REAL, DIMENSION(D%NIT,D%NKT) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
+!
+REAL, DIMENSION(D%NIT,D%NKT) :: ZWK
+!
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT_RHCJ10',0,ZHOOK_HANDLE)
 
@@ -244,11 +247,11 @@ PRSAT_UP(:,:)=PRVM(:,:) ! should be initialised correctly but is (normaly) not u
 ! Initialisation of environment variables at t-dt
 
 ! variables at flux level
-ZTHLM_F(:,:) = MZM_MF(PTHLM(:,:), D%NKA, D%NKU, D%NKL)
-ZRTM_F (:,:) = MZM_MF(PRTM(:,:), D%NKA, D%NKU, D%NKL)
-ZUM_F  (:,:) = MZM_MF(PUM(:,:), D%NKA, D%NKU, D%NKL)
-ZVM_F  (:,:) = MZM_MF(PVM(:,:), D%NKA, D%NKU, D%NKL)
-ZTKEM_F(:,:) = MZM_MF(PTKEM(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PTHLM(:,:), ZTHLM_F(:,:))
+CALL MZM_MF(D, PRTM(:,:), ZRTM_F(:,:))
+CALL MZM_MF(D, PUM(:,:), ZUM_F(:,:))
+CALL MZM_MF(D, PVM(:,:), ZVM_F(:,:))
+CALL MZM_MF(D, PTKEM(:,:), ZTKEM_F(:,:))
 
 ! This updraft is not yet ready to use scalar variables
 !DO JSV=1,ISV
@@ -283,10 +286,10 @@ DO JI=1,D%NIT
   !ZTHS_UP(JI,KKB)=PTHL_UP(JI,KKB)*(1.+XLAMBDA_MF*ZQT_UP(JI))
 ENDDO
 
-ZTHM_F (:,:) = MZM_MF(PTHM (:,:), D%NKA, D%NKU, D%NKL)
-ZPRES_F(:,:) = MZM_MF(PPABSM(:,:), D%NKA, D%NKU, D%NKL)
-ZRHO_F (:,:) = MZM_MF(PRHODREF(:,:), D%NKA, D%NKU, D%NKL)
-ZRVM_F (:,:) = MZM_MF(PRVM(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PTHM (:,:), ZTHM_F(:,:))
+CALL MZM_MF(D, PPABSM(:,:), ZPRES_F(:,:))
+CALL MZM_MF(D, PRHODREF(:,:), ZRHO_F(:,:))
+CALL MZM_MF(D, PRVM(:,:), ZRVM_F(:,:))
 
 ! thetav at mass and flux levels 
 DO JK=1,D%NKT
@@ -336,8 +339,10 @@ GLMIX=.TRUE.
 ZTKEM_F(:,D%NKB)=0.
 !
 IF(TURB%CTURBLEN=='RM17') THEN
-  ZDUDZ = MZF_MF(GZ_M_W_MF(PUM,PDZZ, D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL)
-  ZDVDZ = MZF_MF(GZ_M_W_MF(PVM,PDZZ, D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL)
+  CALL GZ_M_W_MF(D, PUM, PDZZ, ZWK)
+  CALL MZF_MF(D, ZWK, ZDUDZ)
+  CALL GZ_M_W_MF(D, PVM, PDZZ, ZWK)
+  CALL MZF_MF(D, ZWK, ZDVDZ)
   ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
 ELSE
   ZSHEAR = 0. !no shear in bl89 mixing length
diff --git a/src/common/turb/mode_mf_turb.F90 b/src/common/turb/mode_mf_turb.F90
index 4720c6753..60eb017ee 100644
--- a/src/common/turb/mode_mf_turb.F90
+++ b/src/common/turb/mode_mf_turb.F90
@@ -150,15 +150,20 @@ PSVDT = 0.
 !   ( Resulting fluxes are in flux level (w-point) as PEMF and PTHL_UP )
 !
 
-PFLXZTHMF(:,:) = PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(PTHLM(:,:), D%NKA, D%NKU, D%NKL))
+CALL MZM_MF(D, PTHLM(:,:), PFLXZTHMF(:,:))
+PFLXZTHMF(:,:) = PEMF(:,:)*(PTHL_UP(:,:)-PFLXZTHMF(:,:))
 
-PFLXZRMF(:,:) =  PEMF(:,:)*(PRT_UP(:,:)-MZM_MF(PRTM(:,:), D%NKA, D%NKU, D%NKL))
+CALL MZM_MF(D, PRTM(:,:), PFLXZRMF(:,:))
+PFLXZRMF(:,:) =  PEMF(:,:)*(PRT_UP(:,:)-PFLXZRMF(:,:))
 
-PFLXZTHVMF(:,:) = PEMF(:,:)*(PTHV_UP(:,:)-MZM_MF(PTHVM(:,:), D%NKA, D%NKU, D%NKL))
+CALL MZM_MF(D, PTHVM(:,:), PFLXZTHVMF(:,:))
+PFLXZTHVMF(:,:) = PEMF(:,:)*(PTHV_UP(:,:)-PFLXZTHVMF(:,:))
 
 IF (OMIXUV) THEN
-  PFLXZUMF(:,:) =  PEMF(:,:)*(PU_UP(:,:)-MZM_MF(PUM(:,:), D%NKA, D%NKU, D%NKL))
-  PFLXZVMF(:,:) =  PEMF(:,:)*(PV_UP(:,:)-MZM_MF(PVM(:,:), D%NKA, D%NKU, D%NKL))
+  CALL MZM_MF(D, PUM(:,:), PFLXZUMF(:,:))
+  PFLXZUMF(:,:) =  PEMF(:,:)*(PU_UP(:,:)-PFLXZUMF(:,:))
+  CALL MZM_MF(D, PVM(:,:), PFLXZVMF(:,:))
+  PFLXZVMF(:,:) =  PEMF(:,:)*(PV_UP(:,:)-PFLXZVMF(:,:))
 ELSE
   PFLXZUMF(:,:) = 0.
   PFLXZVMF(:,:) = 0.
@@ -180,7 +185,8 @@ ENDIF
 CALL TRIDIAG_MASSFLUX(D,PTHLM,PFLXZTHMF,-PEMF,PTSTEP,PIMPL,  &
                       PDZZ,PRHODJ,ZVARS )
 ! compute new flux
-PFLXZTHMF(:,:) = PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(ZVARS(:,:), D%NKA, D%NKU, D%NKL))
+CALL MZM_MF(D, ZVARS(:,:), PFLXZTHMF(:,:))
+PFLXZTHMF(:,:) = PEMF(:,:)*(PTHL_UP(:,:)-PFLXZTHMF(:,:))
 
 !!! compute THL tendency
 !
@@ -192,7 +198,8 @@ PTHLDT(:,:)= (ZVARS(:,:)-PTHLM(:,:))/PTSTEP
 CALL TRIDIAG_MASSFLUX(D,PRTM(:,:),PFLXZRMF,-PEMF,PTSTEP,PIMPL,  &
                                  PDZZ,PRHODJ,ZVARS )
 ! compute new flux
-PFLXZRMF(:,:) =  PEMF(:,:)*(PRT_UP(:,:)-MZM_MF(ZVARS(:,:), D%NKA, D%NKU, D%NKL))
+CALL MZM_MF(D, ZVARS(:,:), PFLXZRMF(:,:))
+PFLXZRMF(:,:) =  PEMF(:,:)*(PRT_UP(:,:)-PFLXZRMF(:,:))
 
 !!! compute RT tendency
 PRTDT(:,:) = (ZVARS(:,:)-PRTM(:,:))/PTSTEP
@@ -207,7 +214,8 @@ IF (OMIXUV) THEN
   CALL TRIDIAG_MASSFLUX(D,PUM,PFLXZUMF,-PEMF,PTSTEP,PIMPL,  &
                                  PDZZ,PRHODJ,ZVARS )
   ! compute new flux
-  PFLXZUMF(:,:) = PEMF(:,:)*(PU_UP(:,:)-MZM_MF(ZVARS(:,:), D%NKA, D%NKU, D%NKL))
+  CALL MZM_MF(D, ZVARS(:,:), PFLXZUMF(:,:))
+  PFLXZUMF(:,:) = PEMF(:,:)*(PU_UP(:,:)-PFLXZUMF(:,:))
 
   ! compute U tendency
   PUDT(:,:)= (ZVARS(:,:)-PUM(:,:))/PTSTEP
@@ -221,7 +229,8 @@ IF (OMIXUV) THEN
   CALL TRIDIAG_MASSFLUX(D,PVM,PFLXZVMF,-PEMF,PTSTEP,PIMPL,  &
                                  PDZZ,PRHODJ,ZVARS )
   ! compute new flux
-  PFLXZVMF(:,:) = PEMF(:,:)*(PV_UP(:,:)-MZM_MF(ZVARS(:,:), D%NKA, D%NKU, D%NKL))
+  CALL MZM_MF(D, ZVARS(:,:), PFLXZVMF(:,:))
+  PFLXZVMF(:,:) = PEMF(:,:)*(PV_UP(:,:)-PFLXZVMF(:,:))
 
   ! compute V tendency
   PVDT(:,:)= (ZVARS(:,:)-PVM(:,:))/PTSTEP
@@ -237,7 +246,8 @@ DO JSV=1,KSV
   !*     compute mean flux of scalar variables at time t-dt
   !   ( Resulting fluxes are in flux level (w-point) as PEMF and PTHL_UP )
 
-  PFLXZSVMF(:,:,JSV) = PEMF(:,:)*(PSV_UP(:,:,JSV)-MZM_MF(PSVM(:,:,JSV), D%NKA, D%NKU, D%NKL))
+  CALL MZM_MF(D, PSVM(:,:,JSV), PFLXZSVMF(:,:,JSV))
+  PFLXZSVMF(:,:,JSV) = PEMF(:,:)*(PSV_UP(:,:,JSV)-PFLXZSVMF(:,:,JSV))
   
   !
   ! 3.5 Compute the tendency for scalar variables
@@ -246,7 +256,8 @@ DO JSV=1,KSV
   CALL TRIDIAG_MASSFLUX(D,PSVM(:,:,JSV),PFLXZSVMF(:,:,JSV),&
                         -PEMF,PTSTEP,PIMPL,PDZZ,PRHODJ,ZVARS )
   ! compute new flux
-  PFLXZSVMF(:,:,JSV) = PEMF(:,:)*(PSV_UP(:,:,JSV)-MZM_MF(ZVARS, D%NKA, D%NKU, D%NKL))
+  CALL MZM_MF(D, ZVARS, PFLXZSVMF(:,:,JSV))
+  PFLXZSVMF(:,:,JSV) = PEMF(:,:)*(PSV_UP(:,:,JSV)-PFLXZSVMF(:,:,JSV))
 
   ! compute Sv tendency
   PSVDT(:,:,JSV)= (ZVARS(:,:)-PSVM(:,:,JSV))/PTSTEP
diff --git a/src/common/turb/mode_mf_turb_expl.F90 b/src/common/turb/mode_mf_turb_expl.F90
index b865f5925..d0f349a11 100644
--- a/src/common/turb/mode_mf_turb_expl.F90
+++ b/src/common/turb/mode_mf_turb_expl.F90
@@ -128,22 +128,27 @@ PVDT   = 0.
 !          -----------------------------------------------
 !   ( Resulting fluxes are in flux level (w-point) as PEMF and PTHL_UP )
 
-ZRTM_F (:,:) = MZM_MF(PRTM (:,:), D%NKA, D%NKU, D%NKL)
-ZTHLM_F(:,:) = MZM_MF(PTHLM(:,:), D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PRTM (:,:), ZRTM_F(:,:))
+CALL MZM_MF(D, PTHLM(:,:), ZTHLM_F(:,:))
 ZQTM   (:,:) = ZRTM_F (:,:)/(1.+ZRTM_F (:,:))
 ZQT_UP (:,:) = PRT_UP (:,:)/(1.+PRT_UP (:,:))
 ZTHS_UP(:,:) = PTHL_UP(:,:)*(1.+PARAMMF%XLAMBDA_MF*ZQT_UP(:,:))
 ZTHSM  (:,:) = ZTHLM_F(:,:)*(1.+PARAMMF%XLAMBDA_MF*ZQTM(:,:))
 
-PFLXZTHLMF(:,:)  = PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(PTHLM(:,:), D%NKA, D%NKU, D%NKL))  ! ThetaL
-PFLXZRMF(:,:)    = PEMF(:,:)*(PRT_UP (:,:)-MZM_MF(PRTM (:,:), D%NKA, D%NKU, D%NKL))  ! Rt
-PFLXZTHVMF(:,:)  = PEMF(:,:)*(PTHV_UP(:,:)-MZM_MF(PTHVM(:,:), D%NKA, D%NKU, D%NKL))  ! ThetaV
+CALL MZM_MF(D, PTHLM(:,:), PFLXZTHLMF(:,:))
+PFLXZTHLMF(:,:)  = PEMF(:,:)*(PTHL_UP(:,:)-PFLXZTHLMF(:,:))  ! ThetaL
+CALL MZM_MF(D, PRTM (:,:), PFLXZRMF(:,:))
+PFLXZRMF(:,:)    = PEMF(:,:)*(PRT_UP (:,:)-PFLXZRMF(:,:))  ! Rt
+CALL MZM_MF(D, PTHVM(:,:), PFLXZTHVMF(:,:))
+PFLXZTHVMF(:,:)  = PEMF(:,:)*(PTHV_UP(:,:)-PFLXZTHVMF(:,:))  ! ThetaV
 
 ZFLXZTHSMF(:,:)  = PEMF(:,:)*(ZTHS_UP(:,:)-ZTHSM(:,:))    ! Theta S flux
 
 IF (OMIXUV) THEN
-  PFLXZUMF(:,:) =  PEMF(:,:)*(PU_UP(:,:)-MZM_MF(PUM(:,:), D%NKA, D%NKU, D%NKL))  ! U
-  PFLXZVMF(:,:) =  PEMF(:,:)*(PV_UP(:,:)-MZM_MF(PVM(:,:), D%NKA, D%NKU, D%NKL))  ! V
+  CALL MZM_MF(D, PUM(:,:), PFLXZUMF(:,:))
+  PFLXZUMF(:,:) =  PEMF(:,:)*(PU_UP(:,:)-PFLXZUMF(:,:))  ! U
+  CALL MZM_MF(D, PVM(:,:), PFLXZVMF(:,:))
+  PFLXZVMF(:,:) =  PEMF(:,:)*(PV_UP(:,:)-PFLXZVMF(:,:))  ! V
 ELSE
   PFLXZUMF(:,:) = 0.
   PFLXZVMF(:,:) = 0.
diff --git a/src/common/turb/mode_tridiag_massflux.F90 b/src/common/turb/mode_tridiag_massflux.F90
index 230caa4c1..88fd02999 100644
--- a/src/common/turb/mode_tridiag_massflux.F90
+++ b/src/common/turb/mode_tridiag_massflux.F90
@@ -123,7 +123,7 @@ SUBROUTINE TRIDIAG_MASSFLUX(D,PVARM,PF,PDFDT,PTSTEP,PIMPL,  &
 !
 USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
 !
-USE MODI_SHUMAN_MF
+USE MODI_SHUMAN_MF, ONLY: MZM_MF
 !
 IMPLICIT NONE
 !
@@ -160,7 +160,7 @@ INTEGER                              :: JK            ! loop counter
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('TRIDIAG_MASSFLUX',0,ZHOOK_HANDLE)
-ZMZM_RHODJ  = MZM_MF(PRHODJ, D%NKA, D%NKU, D%NKL)
+CALL MZM_MF(D, PRHODJ, ZMZM_RHODJ)
 ZRHODJ_DFDT_O_DZ = ZMZM_RHODJ*PDFDT/PDZZ
 !
 ZA=0.
diff --git a/src/common/turb/shuman_mf.F90 b/src/common/turb/shuman_mf.F90
index 8c1c1ade0..2618e4052 100644
--- a/src/common/turb/shuman_mf.F90
+++ b/src/common/turb/shuman_mf.F90
@@ -8,51 +8,56 @@
 !
 INTERFACE
 !
-FUNCTION DZF_MF(PA, KKA, KKU, KKL)  RESULT(PDZF)
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux side
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZF   ! result at mass
+SUBROUTINE DZF_MF(D, PA, PDZF)
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
+IMPLICIT NONE
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at flux side
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PDZF   ! result at mass
                                                  ! localization
-END FUNCTION DZF_MF
+END SUBROUTINE DZF_MF
 !
-FUNCTION DZM_MF(PA, KKA, KKU, KKL)  RESULT(PDZM)
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass localization
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZM   ! result at flux
+SUBROUTINE DZM_MF(D, PA, PDZM)
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
+IMPLICIT NONE
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PDZM   ! result at flux
                                                  ! side
-END FUNCTION DZM_MF
+END SUBROUTINE DZM_MF
 !
-FUNCTION MZF_MF(PA, KKA, KKU, KKL)  RESULT(PMZF)
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux side
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZF   ! result at mass
+SUBROUTINE MZF_MF(D, PA, PMZF)
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
+IMPLICIT NONE
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at flux side
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PMZF   ! result at mass
                                                  ! localization
-END FUNCTION MZF_MF
-!
-FUNCTION MZM_MF(PA, KKA, KKU, KKL)  RESULT(PMZM)
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass localization
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZM   ! result at flux localization
-END FUNCTION MZM_MF
-!
-FUNCTION GZ_M_W_MF(PY,PDZZ, KKA, KKU, KKL) RESULT(PGZ_M_W)
-REAL, DIMENSION(:,:), INTENT(IN)  :: PDZZ ! Metric coefficient d*zz
-REAL, DIMENSION(:,:), INTENT(IN)  :: PY   ! variable at mass localization
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)  :: KKL  ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2)) :: PGZ_M_W  ! result at flux side
-END FUNCTION GZ_M_W_MF
+END SUBROUTINE MZF_MF
+!
+SUBROUTINE MZM_MF(D, PA, PMZM)
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
+IMPLICIT NONE
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PMZM   ! result at flux localization
+END SUBROUTINE MZM_MF
+!
+SUBROUTINE GZ_M_W_MF(D, PY, PDZZ, PGZ_M_W)
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
+IMPLICIT NONE
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PDZZ ! Metric coefficient d*zz
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PY   ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PGZ_M_W  ! result at flux side
+END SUBROUTINE GZ_M_W_MF
 !
 END INTERFACE
 !
 END MODULE MODI_SHUMAN_MF
 !
 !     ###############################
-      FUNCTION MZF_MF(PA, KKA, KKU, KKL)  RESULT(PMZF)
+      SUBROUTINE MZF_MF(D, PA, PMZF)
 !     ###############################
 !
 !!****  *MZF* -  SHUMAN_MF operator : mean operator in z direction for a
@@ -98,15 +103,15 @@ END MODULE MODI_SHUMAN_MF
 !*       0.    DECLARATIONS
 !              ------------
 !
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
 IMPLICIT NONE
 !
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux side
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZF   ! result at mass
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at flux side
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PMZF   ! result at mass
                                                  ! localization
 !
 !*       0.2   Declarations of local variables
@@ -120,17 +125,17 @@ INTEGER :: JK             ! Loop index in z direction
 !*       1.    DEFINITION OF MZF
 !              ------------------
 !
-DO JK=2,SIZE(PA,2)-1
-  PMZF(:,JK) = 0.5*( PA(:,JK)+PA(:,JK+KKL) )
+DO JK=2,D%NKT-1
+  PMZF(:,JK) = 0.5*( PA(:,JK)+PA(:,JK+D%NKL) )
 END DO
-PMZF(:,KKA) = 0.5*( PA(:,KKA)+PA(:,KKA+KKL) )
-PMZF(:,KKU) = PA(:,KKU)
+PMZF(:,D%NKA) = 0.5*( PA(:,D%NKA)+PA(:,D%NKA+D%NKL) )
+PMZF(:,D%NKU) = PA(:,D%NKU)
 !
 !-------------------------------------------------------------------------------
 !
-END FUNCTION MZF_MF
+END SUBROUTINE MZF_MF
 !     ###############################
-      FUNCTION MZM_MF(PA, KKA, KKU, KKL)  RESULT(PMZM)
+      SUBROUTINE MZM_MF(D, PA, PMZM)
 !     ###############################
 !
 !!****  *MZM* -  SHUMAN_MF operator : mean operator in z direction for a
@@ -176,15 +181,15 @@ END FUNCTION MZF_MF
 !*       0.    DECLARATIONS
 !              ------------
 !
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
 IMPLICIT NONE
 !
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass localization
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZM   ! result at flux localization
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PMZM   ! result at flux localization
 !
 !*       0.2   Declarations of local variables
 !              -------------------------------
@@ -197,17 +202,17 @@ INTEGER :: JK             ! Loop index in z direction
 !*       1.    DEFINITION OF MZM
 !              ------------------
 !
-DO JK=2,SIZE(PA,2)-1
-  PMZM(:,JK) = 0.5*( PA(:,JK)+PA(:,JK-KKL) )
+DO JK=2,D%NKT-1
+  PMZM(:,JK) = 0.5*( PA(:,JK)+PA(:,JK-D%NKL) )
 END DO
-PMZM(:,KKA) = PA(:,KKA)
-PMZM(:,KKU) = 0.5*( PA(:,KKU)+PA(:,KKU-KKL) )
+PMZM(:,D%NKA) = PA(:,D%NKA)
+PMZM(:,D%NKU) = 0.5*( PA(:,D%NKU)+PA(:,D%NKU-D%NKL) )
 !
 !-------------------------------------------------------------------------------
 !
-END FUNCTION MZM_MF
+END SUBROUTINE MZM_MF
 !     ###############################
-      FUNCTION DZF_MF(PA, KKA, KKU, KKL)  RESULT(PDZF)
+      SUBROUTINE DZF_MF(D, PA, PDZF)
 !     ###############################
 !
 !!****  *DZF* -  SHUMAN_MF operator : finite difference operator in z direction
@@ -253,15 +258,15 @@ END FUNCTION MZM_MF
 !*       0.    DECLARATIONS
 !              ------------
 !
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
 IMPLICIT NONE
 !
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux side
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZF   ! result at mass
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at flux side
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PDZF   ! result at mass
                                                  ! localization
 !
 !*       0.2   Declarations of local variables
@@ -274,17 +279,17 @@ INTEGER :: JK           ! Loop index in z direction
 !*       1.    DEFINITION OF DZF
 !              ------------------
 !
-DO JK=2,SIZE(PA,2)-1
-  PDZF(:,JK) = PA(:,JK+KKL) - PA(:,JK)
+DO JK=2,D%NKT-1
+  PDZF(:,JK) = PA(:,JK+D%NKL) - PA(:,JK)
 END DO
-PDZF(:,KKA) = PA(:,KKA+KKL) - PA(:,KKA)
-PDZF(:,KKU) = 0.
+PDZF(:,D%NKA) = PA(:,D%NKA+D%NKL) - PA(:,D%NKA)
+PDZF(:,D%NKU) = 0.
 !
 !-------------------------------------------------------------------------------
 !
-END FUNCTION DZF_MF
+END SUBROUTINE DZF_MF
 !     ###############################
-      FUNCTION DZM_MF(PA, KKA, KKU, KKL)  RESULT(PDZM)
+      SUBROUTINE DZM_MF(D, PA, PDZM)
 !     ###############################
 !
 !!****  *DZM* -  SHUMAN_MF operator : finite difference operator in z direction
@@ -330,15 +335,15 @@ END FUNCTION DZF_MF
 !*       0.    DECLARATIONS
 !              ------------
 !
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
 IMPLICIT NONE
 !
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass localization
-INTEGER,              INTENT(IN)       :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)       :: KKL    ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZM   ! result at flux
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PA     ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PDZM   ! result at flux
                                                  ! side
 !
 !*       0.2   Declarations of local variables
@@ -351,18 +356,18 @@ INTEGER :: JK            ! Loop index in z direction
 !*       1.    DEFINITION OF DZM
 !              ------------------
 !
-DO JK=2,SIZE(PA,2)-1
-  PDZM(:,JK) = PA(:,JK) - PA(:,JK-KKL)
+DO JK=2,D%NKT-1
+  PDZM(:,JK) = PA(:,JK) - PA(:,JK-D%NKL)
 END DO
-PDZM(:,KKA) = 0.
-PDZM(:,KKU) = PA(:,KKU) - PA(:,KKU-KKL)
+PDZM(:,D%NKA) = 0.
+PDZM(:,D%NKU) = PA(:,D%NKU) - PA(:,D%NKU-D%NKL)
 !
 !-------------------------------------------------------------------------------
 !
-END FUNCTION DZM_MF
+END SUBROUTINE DZM_MF
 
 !     ###############################
-      FUNCTION GZ_M_W_MF(PY,PDZZ, KKA, KKU, KKL) RESULT(PGZ_M_W)
+      SUBROUTINE GZ_M_W_MF(D, PY, PDZZ, PGZ_M_W)
 !     ###############################
 !
 !!****  *GZ_M_W * - Compute the gradient along z direction for a
@@ -408,16 +413,16 @@ END FUNCTION DZM_MF
 !
 !-------------------------------------------------------------------------------
 !
+USE MODD_DIMPHYEX,        ONLY: DIMPHYEX_t
 IMPLICIT NONE
 !
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
-REAL, DIMENSION(:,:), INTENT(IN)  :: PDZZ ! Metric coefficient d*zz
-REAL, DIMENSION(:,:), INTENT(IN)  :: PY   ! variable at mass localization
-INTEGER,              INTENT(IN)  :: KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,              INTENT(IN)  :: KKL  ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2)) :: PGZ_M_W  ! result at flux side
+TYPE(DIMPHYEX_t),             INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PDZZ ! Metric coefficient d*zz
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN)  :: PY   ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PGZ_M_W  ! result at flux side
 !
 !*       0.2   Declarations of local variables
 !              -------------------------------
@@ -428,12 +433,12 @@ INTEGER  JK
 !*       1.    COMPUTE THE GRADIENT ALONG Z
 !              -----------------------------
 !
-DO JK=2,SIZE(PY,2)-1
-  PGZ_M_W(:,JK) = (PY(:,JK) - PY(:,JK-KKL)) / PDZZ(:,JK)
+DO JK=2,D%NKT-1
+  PGZ_M_W(:,JK) = (PY(:,JK) - PY(:,JK-D%NKL)) / PDZZ(:,JK)
 END DO
-PGZ_M_W(:,KKA) = 0.
-PGZ_M_W(:,KKU) = (PY(:,KKU) - PY(:,KKU-KKL)) / PDZZ(:,KKU)
+PGZ_M_W(:,D%NKA) = 0.
+PGZ_M_W(:,D%NKU) = (PY(:,D%NKU) - PY(:,D%NKU-D%NKL)) / PDZZ(:,D%NKU)
 !
 !-------------------------------------------------------------------------------
 !
-END FUNCTION GZ_M_W_MF
+END SUBROUTINE GZ_M_W_MF
-- 
GitLab