From ae7d5ce97260e0daadffc512cbcc089e79ac4a62 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Mon, 25 Apr 2022 17:54:21 +0200
Subject: [PATCH] Quentin 25/04/2022: Expand in physical points :
tke_eps_sources
---
src/arome/aux/mode_fill_dimphyexn.F90 | 5 +
src/arome/aux/shuman_phy.F90 | 85 ++++++++++++++++
src/common/aux/modd_dimphyexn.F90 | 5 +
src/common/turb/mode_tke_eps_sources.F90 | 124 ++++++++++++-----------
src/common/turb/mode_tridiag_tke.F90 | 99 +++++++++---------
src/common/turb/turb.F90 | 13 ++-
src/mesonh/aux/mode_fill_dimphyexn.F90 | 23 ++++-
src/mesonh/aux/shuman_phy.f90 | 99 ++++++++++++++++++
src/mesonh/ext/phys_paramn.f90 | 2 +-
9 files changed, 342 insertions(+), 113 deletions(-)
create mode 100644 src/arome/aux/shuman_phy.F90
create mode 100644 src/mesonh/aux/shuman_phy.f90
diff --git a/src/arome/aux/mode_fill_dimphyexn.F90 b/src/arome/aux/mode_fill_dimphyexn.F90
index 5df75ef5e..daa5cf480 100644
--- a/src/arome/aux/mode_fill_dimphyexn.F90
+++ b/src/arome/aux/mode_fill_dimphyexn.F90
@@ -73,6 +73,11 @@ YDDIMPHYEX%NKE=1+KVEXT
YDDIMPHYEX%NKTB=1+KVEXT
YDDIMPHYEX%NKTE=KKT-KVEXT
!
+YDDIMPHYEX%NIBC=YDDIMPHYEX%NIB
+YDDIMPHYEX%NJBC=YDDIMPHYEX%NJB
+YDDIMPHYEX%NIEC=YDDIMPHYEX%NIE
+YDDIMPHYEX%NJEC=YDDIMPHYEX%NJE
+!
IF (LHOOK) CALL DR_HOOK('FILL_DIMPHYEX', 1, ZHOOK_HANDLE)
!
END SUBROUTINE FILL_DIMPHYEX
diff --git a/src/arome/aux/shuman_phy.F90 b/src/arome/aux/shuman_phy.F90
new file mode 100644
index 000000000..f7d9c5bd2
--- /dev/null
+++ b/src/arome/aux/shuman_phy.F90
@@ -0,0 +1,85 @@
+MODULE SHUMAN_PHY
+IMPLICIT NONE
+CONTAINS
+! ###############################
+ SUBROUTINE MZM_PHY(D,PA,PMZM)
+! ###############################
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *MZM* - Shuman operator : mean operator in z direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the z direction (K index) for a field PA localized at a mass
+! point. The result is localized at a z-flux point (w point).
+!
+!!** METHOD
+!! ------
+!! The result PMZM(:,:,k) is defined by 0.5*(PA(:,:,k)+PA(:,:,k-1))
+!! At k=1, PMZM(:,:,1) is defined by -999.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (SHUMAN operators)
+!! Technical specifications Report of The Meso-NH (chapters 3)
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/07/94
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMZM ! result at flux localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JK ! Loop index in z direction
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MZM
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MZM',0,ZHOOK_HANDLE)
+DO JK=2,SIZE(PA,3)-1
+ PMZM(:,:,JK) = 0.5*( PA(:,:,JK)+PA(:,:,JK-D%NKL) )
+END DO
+PMZM(:,:,D%NKA) = -999.
+PMZM(:,:,D%NKU) = 0.5*( PA(:,:,D%NKU)+PA(:,:,D%NKU-D%NKL) )
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MZM',1,ZHOOK_HANDLE)
+END SUBROUTINE MZM_PHY
+END MODULE SHUMAN_PHY
diff --git a/src/common/aux/modd_dimphyexn.F90 b/src/common/aux/modd_dimphyexn.F90
index d644ea4d5..d7493616a 100644
--- a/src/common/aux/modd_dimphyexn.F90
+++ b/src/common/aux/modd_dimphyexn.F90
@@ -66,6 +66,11 @@ TYPE DIMPHYEX_t
!* physical levels only from top of atm to ground: DO JK=NKE, NKB, -KKL
!* all (including non physical) following the array ordering: DO JK=1, NKT
!* physical levels only following the array ordering: DO JK=NKTB, NKTE
+ INTEGER :: NIBC ! Computational indices used in DO LOOP
+ INTEGER :: NJBC ! = NIB/NJC/NIE/NJE in all schemes
+ INTEGER :: NIEC ! except in turbulence where external HALO points must be
+ INTEGER :: NJEC ! included so NIBC=NJBC=1 and NIEC/NJEC=NIT/NJT
+!
END TYPE DIMPHYEX_t
!
END MODULE MODD_DIMPHYEX
diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index 526b14d70..5b5d44e76 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -145,6 +145,8 @@ USE MODE_BUDGET, ONLY: BUDGET_STORE_ADD, BUDGET_STORE_END, BUDGET_STORE_INIT
USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
USE MODE_ll
!
+USE SHUMAN_PHY, ONLY : MZM_PHY
+!
USE MODI_GET_HALO
USE MODI_GRADIENT_M
USE MODI_GRADIENT_U
@@ -219,8 +221,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: &
LOGICAL,DIMENSION(D%NIT,D%NJT,D%NKT) :: GTKENEG
! 3D mask .T. if TKE < CSTURB%XTKEMIN
-INTEGER :: IKB ! Index values for the Beginning and End
- ! mass points of the domain
+INTEGER :: IIE,IIB,IJE,IJB,IKB,IKE ! Index value for the mass points of the domain
!
TYPE(LIST_ll), POINTER :: TZFIELDDISS_ll ! list of fields to exchange
INTEGER :: IINFO_ll ! return code of parallel routine
@@ -237,11 +238,16 @@ NULLIFY(TZFIELDDISS_ll)
IF (LHOOK) CALL DR_HOOK('TKE_EPS_SOURCES',0,ZHOOK_HANDLE)
!
IKB=D%NKB
+IKE=D%NKE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
!
! compute the effective diffusion coefficient at the mass point
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ZKEFF(IIB:IIE,IJB:IJE,IKB:IKE) = PLM(IIB:IIE,IJB:IJE,IKB:IKE) * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKB:IKE))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
!----------------------------------------------------------------------------
!
@@ -254,9 +260,9 @@ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
! Complete the sources of TKE with the horizontal turbulent explicit transport
!
IF (HTURBDIM=='3DIM') THEN
- ZTR=PTRH
+ ZTR(IIB:IIE,IJB:IJE,IKB:IKE)=PTRH(IIB:IIE,IJB:IJE,IKB:IKE)
ELSE
- ZTR=0.
+ ZTR(IIB:IIE,IJB:IJE,IKB:IKE)=0.
END IF
!
!
@@ -265,21 +271,21 @@ END IF
! extrapolate the dynamic production with a 1/Z law from its value at the
! W(IKB+1) value stored in PDP(IKB) to the mass localization tke(IKB)
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+D%NKL)/PDZZ(:,:,IKB))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+PDP(IIB:IIE,IJB:IJE,IKB) = PDP(IIB:IIE,IJB:IJE,IKB) * (1. + PDZZ(IIB:IIE,IJB:IJE,IKB+D%NKL)/PDZZ(IIB:IIE,IJB:IJE,IKB))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
! + (Dynamical Production) + (Thermal Production) - (dissipation)
!
-ZMWORK1 = MZM(ZKEFF, D%NKA, D%NKU, D%NKL)
-ZMWORK2 = MZM(PRHODJ, D%NKA, D%NKU, D%NKL)
+CALL MZM_PHY(D,ZKEFF,ZMWORK1)
+CALL MZM_PHY(D,PRHODJ,ZMWORK2)
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-ZFLX(:,:,:) = CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:)
-ZSOURCE(:,:,:) = ( PRTKES(:,:,:) + PRTKEMS(:,:,:) ) / PRHODJ(:,:,:) &
- - PTKEM(:,:,:) / PTSTEP &
- + PDP(:,:,:) + PTP(:,:,:) + ZTR(:,:,:) - PEXPL * ZFLX(:,:,:) * PTKEM(:,:,:)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ZFLX(IIB:IIE,IJB:IJE,IKB:IKE) = CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)) / PLEPS(IIB:IIE,IJB:IJE,IKB:IKE)
+ZSOURCE(IIB:IIE,IJB:IJE,IKB:IKE) = ( PRTKES(IIB:IIE,IJB:IJE,IKB:IKE) + PRTKEMS(IIB:IIE,IJB:IJE,IKB:IKE) ) / PRHODJ(IIB:IIE,IJB:IJE,IKB:IKE) &
+ - PTKEM(IIB:IIE,IJB:IJE,IKB:IKE) / PTSTEP &
+ + PDP(IIB:IIE,IJB:IJE,IKB:IKE) + PTP(IIB:IIE,IJB:IJE,IKB:IKE) + ZTR(IIB:IIE,IJB:IJE,IKB:IKE) - PEXPL * ZFLX(IIB:IIE,IJB:IJE,IKB:IKE) * PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)
!
!* 2.2 implicit vertical TKE transport
!
@@ -287,8 +293,8 @@ ZSOURCE(:,:,:) = ( PRTKES(:,:,:) + PRTKEMS(:,:,:) ) / PRHODJ(:,:,:) &
! Compute the vector giving the elements just under the diagonal for the
! matrix inverted in TRIDIAG
!
-ZA(:,:,:) = - PTSTEP * CSTURB%XCET * ZMWORK1(:,:,:) * ZMWORK2(:,:,:) / PDZZ(:,:,:)**2
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ZA(IIB:IIE,IJB:IJE,IKB:IKE) = - PTSTEP * CSTURB%XCET * ZMWORK1(IIB:IIE,IJB:IJE,IKB:IKE) * ZMWORK2(IIB:IIE,IJB:IJE,IKB:IKE) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE)**2
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
! Compute TKE at time t+deltat: ( stored in ZRES )
!
@@ -298,10 +304,10 @@ CALL GET_HALO(ZRES)
!* diagnose the dissipation
!
IF (LDIAG_IN_RUN) THEN
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- XCURRENT_TKE_DISS(:,:,:) = ZFLX(:,:,:) * PTKEM(:,:,:) &
- *(PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:))
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ XCURRENT_TKE_DISS(IIB:IIE,IJB:IJE,IKB:IKE) = ZFLX(IIB:IIE,IJB:IJE,IKB:IKE) * PTKEM(IIB:IIE,IJB:IJE,IKB:IKE) &
+ *(PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKB:IKE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKB:IKE))
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
CALL ADD3DFIELD_ll( TZFIELDDISS_ll, XCURRENT_TKE_DISS, 'TKE_EPS_SOURCES::XCURRENT_TKE_DISS' )
CALL UPDATE_HALO_ll(TZFIELDDISS_ll,IINFO_ll)
@@ -311,17 +317,17 @@ ENDIF
! TKE must be greater than its minimum value
! CL : Now done at the end of the time step in ADVECTION_METSV for MesoNH
IF(HPROGRAM/='MESONH') THEN
- !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- GTKENEG(:,:,:) = ZRES(:,:,:) <= CSTURB%XTKEMIN
- WHERE ( GTKENEG(:,:,:) )
- ZRES(:,:,:) = CSTURB%XTKEMIN
+ !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ GTKENEG(IIB:IIE,IJB:IJE,IKB:IKE) = ZRES(IIB:IIE,IJB:IJE,IKB:IKE) <= CSTURB%XTKEMIN
+ WHERE ( GTKENEG(IIB:IIE,IJB:IJE,IKB:IKE) )
+ ZRES(IIB:IIE,IJB:IJE,IKB:IKE) = CSTURB%XTKEMIN
END WHERE
- !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
END IF
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PTDISS(:,:,:) = - ZFLX(:,:,:)*(PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+PTDISS(IIB:IIE,IJB:IJE,IKB:IKE) = - ZFLX(IIB:IIE,IJB:IJE,IKB:IKE)*(PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKB:IKE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKB:IKE))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
IF ( OLES_CALL .OR. &
(OTURB_DIAG .AND. TPFILE%LOPENED) ) THEN
@@ -330,9 +336,9 @@ IF ( OLES_CALL .OR. &
!
ZMWORK1 = MZM(ZKEFF, D%NKA, D%NKU, D%NKL)
ZDWORK1 = DZM(PIMPL * ZRES + PEXPL * PTKEM, D%NKA, D%NKU, D%NKL)
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- ZFLX(:,:,:) = - CSTURB%XCET * ZMWORK1(:,:,:) * ZDWORK1(:,:,:) / PDZZ(:,:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ ZFLX(IIB:IIE,IJB:IJE,IKB:IKE) = - CSTURB%XCET * ZMWORK1(IIB:IIE,IJB:IJE,IKB:IKE) * ZDWORK1(IIB:IIE,IJB:IJE,IKB:IKE) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
ZFLX(:,:,IKB) = 0.
ZFLX(:,:,D%NKA) = 0.
@@ -340,9 +346,9 @@ IF ( OLES_CALL .OR. &
! Compute the whole turbulent TRansport of TKE:
!
ZDWORK1 = DZF(MZM(PRHODJ, D%NKA, D%NKU, D%NKL) * ZFLX / PDZZ, D%NKA, D%NKU, D%NKL)
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- ZTR(:,:,:)= ZTR(:,:,:) - ZDWORK1(:,:,:) /PRHODJ(:,:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ ZTR(IIB:IIE,IJB:IJE,IKB:IKE)= ZTR(IIB:IIE,IJB:IJE,IKB:IKE) - ZDWORK1(IIB:IIE,IJB:IJE,IKB:IKE) /PRHODJ(IIB:IIE,IJB:IJE,IKB:IKE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
! Storage in the LES configuration
!
@@ -357,9 +363,9 @@ END IF
!
IF (BUCONF%LBUDGET_TKE) THEN
! Dynamical production
- CALL BUDGET_STORE_ADD( TBUDGETS(NBUDGET_TKE), 'DP', PDP(:, :, :) * PRHODJ(:, :, :) )
+ CALL BUDGET_STORE_ADD( TBUDGETS(NBUDGET_TKE), 'DP', PDP(:,:,:) * PRHODJ(:,:,:) )
! Thermal production
- CALL BUDGET_STORE_ADD( TBUDGETS(NBUDGET_TKE), 'TP', PTP(:, :, :) * PRHODJ(:, :, :) )
+ CALL BUDGET_STORE_ADD( TBUDGETS(NBUDGET_TKE), 'TP', PTP(:,:,:) * PRHODJ(:,:,:) )
! Dissipation
CALL BUDGET_STORE_ADD( TBUDGETS(NBUDGET_TKE), 'DISS',- CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * &
(PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:)) * PRHODJ(:,:,:))
@@ -370,20 +376,20 @@ END IF
!Store the previous source terms in prtkes before initializing the next one
!Should be in IF LBUDGET_TKE only. Was removed out for a correct comput. of PTDIFF in case of LBUDGET_TKE=F in AROME
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PRTKES(:,:,:) = PRTKES(:,:,:) + PRHODJ(:,:,:) * &
- ( PDP(:,:,:) + PTP(:,:,:) &
- - CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * ( PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:) ) )
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+PRTKES(IIB:IIE,IJB:IJE,IKB:IKE) = PRTKES(IIB:IIE,IJB:IJE,IKB:IKE) + PRHODJ(IIB:IIE,IJB:IJE,IKB:IKE) * &
+ ( PDP(IIB:IIE,IJB:IJE,IKB:IKE) + PTP(IIB:IIE,IJB:IJE,IKB:IKE) &
+ - CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)) / PLEPS(IIB:IIE,IJB:IJE,IKB:IKE) * ( PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKB:IKE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKB:IKE) ) )
!
-PTDIFF(:,:,:) = ZRES(:,:,:) / PTSTEP - PRTKES(:,:,:)/PRHODJ(:,:,:) &
- & - PDP(:,:,:)- PTP(:,:,:) - PTDISS(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+PTDIFF(IIB:IIE,IJB:IJE,IKB:IKE) = ZRES(IIB:IIE,IJB:IJE,IKB:IKE) / PTSTEP - PRTKES(IIB:IIE,IJB:IJE,IKB:IKE)/PRHODJ(IIB:IIE,IJB:IJE,IKB:IKE) &
+ & - PDP(IIB:IIE,IJB:IJE,IKB:IKE)- PTP(IIB:IIE,IJB:IJE,IKB:IKE) - PTDISS(IIB:IIE,IJB:IJE,IKB:IKE)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TKE), 'TR', PRTKES(:, :, :) )
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PRTKES(:,:,:) = ZRES(:,:,:) * PRHODJ(:,:,:) / PTSTEP - PRTKEMS(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+PRTKES(IIB:IIE,IJB:IJE,IKB:IKE) = ZRES(IIB:IIE,IJB:IJE,IKB:IKE) * PRHODJ(IIB:IIE,IJB:IJE,IKB:IKE) / PTSTEP - PRTKEMS(IIB:IIE,IJB:IJE,IKB:IKE)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!
! stores the whole turbulent transport
!
@@ -394,10 +400,10 @@ IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TKE), 'TR', PRTK
!* 3. COMPUTE THE DISSIPATIVE HEATING
! -------------------------------
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PRTHLS(:,:,:) = PRTHLS(:,:,:) + CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * &
- (PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:)) * PRHODJ(:,:,:) * PCOEF_DISS(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+PRTHLS(IIB:IIE,IJB:IJE,IKB:IKE) = PRTHLS(IIB:IIE,IJB:IJE,IKB:IKE) + CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)) / PLEPS(IIB:IIE,IJB:IJE,IKB:IKE) * &
+ (PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKB:IKE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKB:IKE)) * PRHODJ(IIB:IIE,IJB:IJE,IKB:IKE) * PCOEF_DISS(IIB:IIE,IJB:IJE,IKB:IKE)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
!----------------------------------------------------------------------------
!
!* 4. STORES SOME DIAGNOSTICS
@@ -405,15 +411,15 @@ PRTHLS(:,:,:) = PRTHLS(:,:,:) + CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:)
!
IF(PRESENT(PTR)) PTR=ZTR
IF(PRESENT(PDISS)) THEN
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- PDISS(:,:,:) = -CSTURB%XCED * (PTKEM(:,:,:)**1.5) / PLEPS(:,:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ PDISS(IIB:IIE,IJB:IJE,IKB:IKE) = -CSTURB%XCED * (PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)**1.5) / PLEPS(IIB:IIE,IJB:IJE,IKB:IKE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
END IF
!
IF(PRESENT(PEDR)) THEN
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- PEDR(:,:,:) = CSTURB%XCED * (PTKEM(:,:,:)**1.5) / PLEPS(:,:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
+ PEDR(IIB:IIE,IJB:IJE,IKB:IKE) = CSTURB%XCED * (PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)**1.5) / PLEPS(IIB:IIE,IJB:IJE,IKB:IKE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKE)
END IF
!
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
diff --git a/src/common/turb/mode_tridiag_tke.F90 b/src/common/turb/mode_tridiag_tke.F90
index e69a64f73..0832e1060 100644
--- a/src/common/turb/mode_tridiag_tke.F90
+++ b/src/common/turb/mode_tridiag_tke.F90
@@ -142,6 +142,7 @@ INTEGER :: JI,JJ,JK ! loop counter
INTEGER :: IKB,IKE ! inner vertical limits
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+INTEGER :: IIE,IIB,IJE,IJB
!
! ---------------------------------------------------------------------------
!
@@ -156,27 +157,31 @@ IKTB=D%NKTB
IKTE=D%NKTE
IKB=D%NKB
IKE=D%NKE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZY(:,:,IKB) = PVARM(:,:,IKB) + PTSTEP*PSOURCE(:,:,IKB) - &
- PEXPL / PRHODJ(:,:,IKB) * PA(:,:,IKB+D%NKL) * (PVARM(:,:,IKB+D%NKL) - PVARM(:,:,IKB))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZY(IIB:IIE,IJB:IJE,IKB) = PVARM(IIB:IIE,IJB:IJE,IKB) + PTSTEP*PSOURCE(IIB:IIE,IJB:IJE,IKB) - &
+ PEXPL / PRHODJ(IIB:IIE,IJB:IJE,IKB) * PA(IIB:IIE,IJB:IJE,IKB+D%NKL) * (PVARM(IIB:IIE,IJB:IJE,IKB+D%NKL) - PVARM(IIB:IIE,IJB:IJE,IKB))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
DO JK=IKTB+1,IKTE-1
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZY(:,:,JK)= PVARM(:,:,JK) + PTSTEP*PSOURCE(:,:,JK) - &
- PEXPL / PRHODJ(:,:,JK) * &
- ( PVARM(:,:,JK-D%NKL)*PA(:,:,JK) &
- -PVARM(:,:,JK)*(PA(:,:,JK)+PA(:,:,JK+D%NKL)) &
- +PVARM(:,:,JK+D%NKL)*PA(:,:,JK+D%NKL) &
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZY(IIB:IIE,IJB:IJE,JK)= PVARM(IIB:IIE,IJB:IJE,JK) + PTSTEP*PSOURCE(IIB:IIE,IJB:IJE,JK) - &
+ PEXPL / PRHODJ(IIB:IIE,IJB:IJE,JK) * &
+ ( PVARM(IIB:IIE,IJB:IJE,JK-D%NKL)*PA(IIB:IIE,IJB:IJE,JK) &
+ -PVARM(IIB:IIE,IJB:IJE,JK)*(PA(IIB:IIE,IJB:IJE,JK)+PA(IIB:IIE,IJB:IJE,JK+D%NKL)) &
+ +PVARM(IIB:IIE,IJB:IJE,JK+D%NKL)*PA(IIB:IIE,IJB:IJE,JK+D%NKL) &
)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END DO
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZY(:,:,IKE)= PVARM(:,:,IKE) + PTSTEP*PSOURCE(:,:,IKE) + &
- PEXPL / PRHODJ(:,:,IKE) * PA(:,:,IKE) * (PVARM(:,:,IKE)-PVARM(:,:,IKE-D%NKL))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZY(IIB:IIE,IJB:IJE,IKE)= PVARM(IIB:IIE,IJB:IJE,IKE) + PTSTEP*PSOURCE(IIB:IIE,IJB:IJE,IKE) + &
+ PEXPL / PRHODJ(IIB:IIE,IJB:IJE,IKE) * PA(IIB:IIE,IJB:IJE,IKE) * (PVARM(IIB:IIE,IJB:IJE,IKE)-PVARM(IIB:IIE,IJB:IJE,IKE-D%NKL))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
!
!* 2. INVERSION OF THE TRIDIAGONAL SYSTEM
@@ -187,55 +192,55 @@ IF ( PIMPL > 1.E-10 ) THEN
!
! going up
!
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZBET(:,:) = 1. + PIMPL * (PDIAG(:,:,IKB)-PA(:,:,IKB+D%NKL) / PRHODJ(:,:,IKB))
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZBET(IIB:IIE,IJB:IJE) = 1. + PIMPL * (PDIAG(IIB:IIE,IJB:IJE,IKB)-PA(IIB:IIE,IJB:IJE,IKB+D%NKL) / PRHODJ(IIB:IIE,IJB:IJE,IKB))
! bet = b(ikb)
- PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ PVARP(IIB:IIE,IJB:IJE,IKB) = ZY(IIB:IIE,IJB:IJE,IKB) / ZBET(IIB:IIE,IJB:IJE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
DO JK = IKB+D%NKL,IKE-D%NKL,D%NKL
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZGAM(:,:,JK) = PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK-D%NKL) / ZBET(:,:)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZGAM(IIB:IIE,IJB:IJE,JK) = PIMPL * PA(IIB:IIE,IJB:IJE,JK) / PRHODJ(IIB:IIE,IJB:IJE,JK-D%NKL) / ZBET(IIB:IIE,IJB:IJE)
! gam(k) = c(k-1) / bet
- ZBET(:,:) = 1. + PIMPL * ( PDIAG(:,:,JK) - &
- ( PA(:,:,JK) * (1. + ZGAM(:,:,JK)) &
- + PA(:,:,JK+D%NKL) &
- ) / PRHODJ(:,:,JK) &
+ ZBET(IIB:IIE,IJB:IJE) = 1. + PIMPL * ( PDIAG(IIB:IIE,IJB:IJE,JK) - &
+ ( PA(IIB:IIE,IJB:IJE,JK) * (1. + ZGAM(IIB:IIE,IJB:IJE,JK)) &
+ + PA(IIB:IIE,IJB:IJE,JK+D%NKL) &
+ ) / PRHODJ(IIB:IIE,IJB:IJE,JK) &
) ! bet = b(k) - a(k)* gam(k)
- PVARP(:,:,JK)= ( ZY(:,:,JK) - PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK) &
- * PVARP(:,:,JK-D%NKL) &
- ) / ZBET(:,:)
+ PVARP(IIB:IIE,IJB:IJE,JK)= ( ZY(IIB:IIE,IJB:IJE,JK) - PIMPL * PA(IIB:IIE,IJB:IJE,JK) / PRHODJ(IIB:IIE,IJB:IJE,JK) &
+ * PVARP(IIB:IIE,IJB:IJE,JK-D%NKL) &
+ ) / ZBET(IIB:IIE,IJB:IJE)
! res(k) = (y(k) -a(k)*res(k-1))/ bet
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END DO
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
! special treatment for the last level
- ZGAM(:,:,IKE) = PIMPL * PA(:,:,IKE) / PRHODJ(:,:,IKE-D%NKL) / ZBET(:,:)
+ ZGAM(IIB:IIE,IJB:IJE,IKE) = PIMPL * PA(IIB:IIE,IJB:IJE,IKE) / PRHODJ(IIB:IIE,IJB:IJE,IKE-D%NKL) / ZBET(IIB:IIE,IJB:IJE)
! gam(k) = c(k-1) / bet
- ZBET(:,:) = 1. + PIMPL * ( PDIAG(:,:,IKE) - &
- ( PA(:,:,IKE) * (1. + ZGAM(:,:,IKE)) ) / PRHODJ(:,:,IKE) &
+ ZBET(IIB:IIE,IJB:IJE) = 1. + PIMPL * ( PDIAG(IIB:IIE,IJB:IJE,IKE) - &
+ ( PA(IIB:IIE,IJB:IJE,IKE) * (1. + ZGAM(IIB:IIE,IJB:IJE,IKE)) ) / PRHODJ(IIB:IIE,IJB:IJE,IKE) &
)
! bet = b(k) - a(k)* gam(k)
- PVARP(:,:,IKE)= ( ZY(:,:,IKE) - PIMPL * PA(:,:,IKE) / PRHODJ(:,:,IKE) &
- * PVARP(:,:,IKE-D%NKL) &
- ) / ZBET(:,:)
+ PVARP(IIB:IIE,IJB:IJE,IKE)= ( ZY(IIB:IIE,IJB:IJE,IKE) - PIMPL * PA(IIB:IIE,IJB:IJE,IKE) / PRHODJ(IIB:IIE,IJB:IJE,IKE) &
+ * PVARP(IIB:IIE,IJB:IJE,IKE-D%NKL) &
+ ) / ZBET(IIB:IIE,IJB:IJE)
! res(k) = (y(k) -a(k)*res(k-1))/ bet
!
! going down
!
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
DO JK = IKE-D%NKL,IKB,-1*D%NKL
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+D%NKL) * PVARP(:,:,JK+D%NKL)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ PVARP(IIB:IIE,IJB:IJE,JK) = PVARP(IIB:IIE,IJB:IJE,JK) - ZGAM(IIB:IIE,IJB:IJE,JK+D%NKL) * PVARP(IIB:IIE,IJB:IJE,JK+D%NKL)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END DO
!
ELSE
!
DO JK=IKTB,IKTE
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- PVARP(:,:,JK) = ZY(:,:,JK)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ PVARP(IIB:IIE,IJB:IJE,JK) = ZY(IIB:IIE,IJB:IJE,JK)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END DO
!
END IF
@@ -244,10 +249,10 @@ END IF
!* 3. FILL THE UPPER AND LOWER EXTERNAL VALUES
! ----------------------------------------
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-PVARP(:,:,D%NKA)=PVARP(:,:,IKB)
-PVARP(:,:,D%NKU)=PVARP(:,:,IKE)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+PVARP(IIB:IIE,IJB:IJE,D%NKA)=PVARP(IIB:IIE,IJB:IJE,IKB)
+PVARP(IIB:IIE,IJB:IJE,D%NKU)=PVARP(IIB:IIE,IJB:IJE,IKE)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
!-------------------------------------------------------------------------------
!
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index 639ea970f..b9e96fd9c 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -459,7 +459,7 @@ REAL, DIMENSION(D%NIT,D%NJT) :: ZTAU11M,ZTAU12M, &
REAL :: ZEXPL ! 1-PIMPL deg of expl.
REAL :: ZRVORD ! RV/RD
!
-INTEGER :: IKB,IKE ! index value for the
+INTEGER :: IIE,IIB,IJE,IJB,IKB,IKE ! index value for the
! Beginning and the End of the physical domain for the mass points
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
@@ -494,6 +494,11 @@ IKTB=D%NKTB
IKTE=D%NKTE
IKB=D%NKB
IKE=D%NKE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!print*,"MINMAX PRTKES = ",MINVAL(PRTKES),MAXVAL(PRTKES)
!
ZEXPL = 1.- PIMPL
ZRVORD= CST%XRV / CST%XRD
@@ -1074,13 +1079,13 @@ END IF
! cloud computation is not statistical
ZWORK1 = MZF(PFLXZTHVMF,D%NKA, D%NKU, D%NKL)
!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PTP(:,:,:) = PTP(:,:,:) + CST%XG / PTHVREF(:,:,:) * ZWORK1(:,:,:)
+PTP(IIB:IIE,IJB:IJE,IKB:IKE) = PTP(IIB:IIE,IJB:IJE,IKB:IKE) + CST%XG / PTHVREF(IIB:IIE,IJB:IJE,IKB:IKE) * ZWORK1(IIB:IIE,IJB:IJE,IKB:IKE)
!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
IF(PRESENT(PTPMF)) THEN
ZWORK1 = MZF(PFLXZTHVMF, D%NKA, D%NKU, D%NKL)
!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- PTPMF(:,:,:)=CST%XG / PTHVREF(:,:,:) * ZWORK1(:,:,:)
+ PTPMF(IIB:IIE,IJB:IJE,IKB:IKE)=CST%XG / PTHVREF(IIB:IIE,IJB:IJE,IKB:IKE) * ZWORK1(IIB:IIE,IJB:IJE,IKB:IKE)
!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
END IF
! 6.2 TKE evolution equation
@@ -1283,7 +1288,7 @@ IF (OLES_CALL) THEN
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
!
-IF(PRESENT(PLEM)) PLEM(:,:,:) = ZLM(:,:,:)
+IF(PRESENT(PLEM)) PLEM(IIB:IIE,IJB:IJE,IKB:IKE) = ZLM(IIB:IIE,IJB:IJE,IKB:IKE)
!----------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
diff --git a/src/mesonh/aux/mode_fill_dimphyexn.F90 b/src/mesonh/aux/mode_fill_dimphyexn.F90
index af7b191d0..351c96d49 100644
--- a/src/mesonh/aux/mode_fill_dimphyexn.F90
+++ b/src/mesonh/aux/mode_fill_dimphyexn.F90
@@ -6,7 +6,7 @@
MODULE MODE_FILL_DIMPHYEX
IMPLICIT NONE
CONTAINS
-SUBROUTINE FILL_DIMPHYEX(YDDIMPHYEX, KIT, KJT, KKT)
+SUBROUTINE FILL_DIMPHYEX(YDDIMPHYEX, KIT, KJT, KKT, LTURB)
! #########################
!
!!
@@ -46,7 +46,9 @@ IMPLICIT NONE
!
TYPE(DIMPHYEX_t), INTENT(OUT) :: YDDIMPHYEX ! Structure to fill in
INTEGER, INTENT(IN) :: KIT, KJT, KKT ! Array dimensions
-
+LOGICAL, INTENT(IN), OPTIONAL :: LTURB ! Flag to replace array dimensions I/JB and I/JE to the full array size
+ ! needed if computation in HALO points (e.g. in turbulence)
+LOGICAL :: YTURB
!
!* 0.2 declaration of local variables
!
@@ -69,6 +71,23 @@ YDDIMPHYEX%NKE=KKT-JPVEXT
YDDIMPHYEX%NKTB=1+JPVEXT
YDDIMPHYEX%NKTE=KKT-JPVEXT
!
+IF(PRESENT(LTURB)) THEN
+ YTURB=LTURB
+ELSE
+ YTURB=.FALSE. ! Default value
+END IF
+!
+IF(YTURB) THEN
+ YDDIMPHYEX%NIBC=1
+ YDDIMPHYEX%NJBC=1
+ YDDIMPHYEX%NIEC=YDDIMPHYEX%NIT
+ YDDIMPHYEX%NJEC=YDDIMPHYEX%NJT
+ELSE
+ YDDIMPHYEX%NIBC=YDDIMPHYEX%NIB
+ YDDIMPHYEX%NJBC=YDDIMPHYEX%NJB
+ YDDIMPHYEX%NIEC=YDDIMPHYEX%NIE
+ YDDIMPHYEX%NJEC=YDDIMPHYEX%NJE
+END IF
IF (LHOOK) CALL DR_HOOK('FILL_DIMPHYEX', 1, ZHOOK_HANDLE)
!
END SUBROUTINE FILL_DIMPHYEX
diff --git a/src/mesonh/aux/shuman_phy.f90 b/src/mesonh/aux/shuman_phy.f90
new file mode 100644
index 000000000..28feda3c7
--- /dev/null
+++ b/src/mesonh/aux/shuman_phy.f90
@@ -0,0 +1,99 @@
+MODULE SHUMAN_PHY
+IMPLICIT NONE
+CONTAINS
+! ###############################
+ SUBROUTINE MZM_PHY(D,PA,PMZM)
+! ###############################
+!
+!!**** *MZM* - Shuman operator : mean operator in z direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the z direction (K index) for a field PA localized at a mass
+! point. The result is localized at a z-flux point (w point).
+!
+!!** METHOD
+!! ------
+!! The result PMZM(:,:,k) is defined by 0.5*(PA(:,:,k)+PA(:,:,k-1))
+!! At k=1, PMZM(:,:,1) is defined by -999.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (SHUMAN operators)
+!! Technical specifications Report of The Meso-NH (chapters 3)
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/07/94
+!! optimisation 20/08/00 J. Escobar
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PMZM ! result at flux localization
+
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JK ! Loop index in z direction
+INTEGER :: IKU ! upper bound in z direction of PA
+!
+INTEGER :: IIU,IJU
+INTEGER :: JIJ,JI,JJ
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MZM
+! ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+JIJKOR = 1 + IIU*IJU
+JIJKEND = IIU*IJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+ PMZM(JIJK,1,1) = 0.5*( PA(JIJK,1,1)+PA(JIJK-IIU*IJU,1,1) )
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJ=1,IIU*IJU
+ PMZM(JIJ,1,1) = -999.
+END DO
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MZM_PHY
+END MODULE SHUMAN_PHY
diff --git a/src/mesonh/ext/phys_paramn.f90 b/src/mesonh/ext/phys_paramn.f90
index 59162c9b7..5e4565959 100644
--- a/src/mesonh/ext/phys_paramn.f90
+++ b/src/mesonh/ext/phys_paramn.f90
@@ -476,7 +476,7 @@ IKB = 1 + JPVEXT
IKE = IKU - JPVEXT
!
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
-CALL FILL_DIMPHYEX(YLDIMPHYEX, SIZE(XTHT,1), SIZE(XTHT,2), SIZE(XTHT,3))
+CALL FILL_DIMPHYEX(YLDIMPHYEX, SIZE(XTHT,1), SIZE(XTHT,2), SIZE(XTHT,3),.TRUE.)
!
ZTIME1 = 0.0_MNHTIME
ZTIME2 = 0.0_MNHTIME
--
GitLab