diff --git a/src/arome/aux/shuman_phy.F90 b/src/arome/aux/shuman_phy.F90
index f7d9c5bd2d7acb3101d1399e01374cc5d541fd57..930a2628011633d5aab2ffc6e305a8f77a9b1fb6 100644
--- a/src/arome/aux/shuman_phy.F90
+++ b/src/arome/aux/shuman_phy.F90
@@ -1,6 +1,194 @@
MODULE SHUMAN_PHY
IMPLICIT NONE
CONTAINS
+! ###############################
+ SUBROUTINE MYF_PHY(D,PA,PMYF)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *MYF* - Shuman operator : mean operator in y direction for a
+!! variable at a flux side
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the y direction (J index) for a field PA localized at a y-flux
+! point (v point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PMYF(i,:,:) is defined by 0.5*(PA(:,j,:)+PA(:,j+1,:))
+!! At j=size(PA,2), PMYF(:,j,:) are replaced by the values of PMYF,
+!! which are the right values in the y-cyclic case
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+!
+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) :: PMYF ! result at flux localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JJ ! Loop index in y direction
+INTEGER :: IJU ! upper bound in y direction of PA
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MYF
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MYF',0,ZHOOK_HANDLE)
+IJU = SIZE(PA,2)
+
+!POUR AROME
+
+PMYF=PA
+
+!
+!DO JJ=1,IJU-1
+! PMYF(:,JJ,:) = 0.5*( PA(:,JJ,:)+PA(:,JJ+1,:) )
+!END DO
+!
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MYF',1,ZHOOK_HANDLE)
+END SUBROUTINE MYF_PHY
+!
+!
+! ###############################
+ SUBROUTINE MYM_PHY(D,PA,PMYM)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *MYM* - Shuman operator : mean operator in y direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the y direction (J index) for a field PA localized at a mass
+! point. The result is localized at a y-flux point (v point).
+!
+!!** METHOD
+!! ------
+!! The result PMYM(:,j,:) is defined by 0.5*(PA(:,j,:)+PA(:,j-1,:))
+!! At j=1, PMYM(:,j,:) are replaced by the values of PMYM,
+!! which are the right values in the y-cyclic case.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+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) :: PMYM ! result at flux localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JJ ! Loop index in y direction
+INTEGER :: IJU ! Size of the array in the y direction
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MYM
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MYM',0,ZHOOK_HANDLE)
+IJU=SIZE(PA,2)
+
+!POUR AROME
+
+PMYM=PA
+
+!
+!DO JJ=2,IJU
+! PMYM(:,JJ,:) = 0.5*( PA(:,JJ,:)+PA(:,JJ-1,:) )
+!END DO
+!
+!PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:)
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MYM',1,ZHOOK_HANDLE)
+END SUBROUTINE MYM_PHY
! ###############################
SUBROUTINE MZM_PHY(D,PA,PMZM)
! ###############################
@@ -82,4 +270,350 @@ 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
+! ###############################
+ SUBROUTINE DZM_PHY(D,PA,PDZM)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *DZM* - Shuman operator : finite difference operator in z direction
+!! for a variable at a mass localization
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a finite difference
+! 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 PDZM(:,j,:) is defined by (PA(:,:,k)-PA(:,:,k-1))
+!! At k=1, PDZM(:,:,k) 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 05/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) :: PDZM ! result at flux
+ ! side
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JK ! Loop index in z direction
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF DZM
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('DZM',0,ZHOOK_HANDLE)
+DO JK=2,SIZE(PA,3)-1
+ PDZM(:,:,JK) = PA(:,:,JK) - PA(:,:,JK-D%NKL)
+END DO
+PDZM(:,:,D%NKA) = -999.
+PDZM(:,:,D%NKU) = PA(:,:,D%NKU) - PA(:,:,D%NKU-D%NKL)
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('DZM',1,ZHOOK_HANDLE)
+END SUBROUTINE DZM_PHY
+
+! ###############################
+ SUBROUTINE MXM_PHY(D,PA,PMXM)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *MXM* - Shuman operator : mean operator in x direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the x direction (I index) for a field PA localized at a mass
+! point. The result is localized at a x-flux point (u point).
+!
+!!** METHOD
+!! ------
+!! The result PMXM(i,:,:) is defined by 0.5*(PA(i,:,:)+PA(i-1,:,:))
+!! At i=1, PMXM(1,:,:) are replaced by the values of PMXM,
+!! which are the right values in the x-cyclic case.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_PARAMETERS
+!
+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) :: PMXM ! result at flux localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JI ! Loop index in x direction
+INTEGER :: IIU ! Size of the array in the x direction
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MXM
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MXM',0,ZHOOK_HANDLE)
+IIU = SIZE(PA,1)
+!
+!POUR AROME
+
+PMXM=PA
+
+!
+!DO JI=2,IIU
+! PMXM(JI,:,:) = 0.5*( PA(JI,:,:)+PA(JI-1,:,:) )
+!END DO
+!
+!PMXM(1,:,:) = PMXM(IIU-2*JPHEXT+1,:,:)
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MXM',1,ZHOOK_HANDLE)
+END SUBROUTINE MXM_PHY
+! ###############################
+ SUBROUTINE MXF_PHY(D,PA,PMXF)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *MXF* - Shuman operator : mean operator in x direction for a
+!! variable at a flux side
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the x direction (I index) for a field PA localized at a x-flux
+! point (u point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PMXF(i,:,:) is defined by 0.5*(PA(i,:,:)+PA(i+1,:,:))
+!! At i=size(PA,1), PMXF(i,:,:) are replaced by the values of PMXF,
+!! which are the right values in the x-cyclic case
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!* 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) :: PMXF ! result at flux localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JI ! Loop index in x direction
+INTEGER :: IIU ! upper bound in x direction of PA
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MXF
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MXF',0,ZHOOK_HANDLE)
+IIU = SIZE(PA,1)
+!
+!POUR AROME
+
+PMXF=PA
+
+!
+!DO JI=1,IIU-1
+! PMXF(JI,:,:) = 0.5*( PA(JI,:,:)+PA(JI+1,:,:) )
+!END DO
+!
+!PMXF(IIU,:,:) = PMXF(2*JPHEXT,:,:)
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MXF',1,ZHOOK_HANDLE)
+END SUBROUTINE MXF_PHY
+! ###############################
+ SUBROUTINE MZF_PHY(D,PA,PMZF)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *MZF* - Shuman operator : mean operator in z direction for a
+!! variable at a flux side
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the z direction (K index) for a field PA localized at a z-flux
+! point (w point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PMZF(:,:,k) is defined by 0.5*(PA(:,:,k)+PA(:,:,k+1))
+!! At k=size(PA,3), PMZF(:,:,k) 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 flux localization
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMZF ! result at mass localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JK ! Loop index in z direction
+INTEGER :: IKT ! upper bound in z direction of PA
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MZF
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MZF',0,ZHOOK_HANDLE)
+IKT = SIZE(PA,3)
+DO JK=2,IKT-1
+ PMZF(:,:,JK) = 0.5*( PA(:,:,JK)+PA(:,:,JK+D%NKL) )
+END DO
+PMZF(:,:,D%NKU) = -999.
+PMZF(:,:,D%NKA) = 0.5*( PA(:,:,D%NKA)+PA(:,:,D%NKA+D%NKL) )
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MZF',1,ZHOOK_HANDLE)
+END SUBROUTINE MZF_PHY
+
END MODULE SHUMAN_PHY
diff --git a/src/common/aux/gradient_u_phy.F90 b/src/common/aux/gradient_u_phy.F90
new file mode 100644
index 0000000000000000000000000000000000000000..7059cbec803bb88dd13e0f7feb5252d42d898bf3
--- /dev/null
+++ b/src/common/aux/gradient_u_phy.F90
@@ -0,0 +1,94 @@
+MODULE MODE_GRADIENT_U_PHY
+IMPLICIT NONE
+CONTAINS
+! #######################################################
+ SUBROUTINE GZ_U_UW_PHY(D,PA,PDZZ,PGZ_U_UW)
+! #######################################################
+!
+!!**** *GZ_U_UW - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Z
+!! direction for a variable placed at the
+!! U point and the result is placed at
+!! the UW vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Z cartesian direction for a field PA placed at the
+! U point. The result is placed at the UW vorticity point.
+!
+! dzm(PA)
+! PGZ_U_UW = ------
+! ____x
+! d*zz
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MXM : Shuman functions (mean operators)
+!! DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE SHUMAN_PHY, ONLY : DZM_PHY, MXM_PHY
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGZ_U_UW ! result UW point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: PA_WORK, PDZZ_WORK
+!
+INTEGER :: JI,JJ,JK
+
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GZ_U_UW
+! ---------------------
+!
+CALL DZM_PHY(D,PA,PA_WORK)
+CALL MXM_PHY(D,PDZZ,PDZZ_WORK)
+!
+!$mnh__expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
+PGZ_U_UW(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)= PA_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT) / PDZZ_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)
+!$mnh_end_expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
+!
+!----------------------------------------------------------------------------
+!
+END SUBROUTINE GZ_U_UW_PHY
+END MODULE MODE_GRADIENT_U_PHY
diff --git a/src/common/aux/gradient_v_phy.F90 b/src/common/aux/gradient_v_phy.F90
new file mode 100644
index 0000000000000000000000000000000000000000..0ed5a5044fa66a1ec9160f601ff768af138c6ffe
--- /dev/null
+++ b/src/common/aux/gradient_v_phy.F90
@@ -0,0 +1,93 @@
+MODULE MODE_GRADIENT_V_PHY
+IMPLICIT NONE
+CONTAINS
+ ! #######################################################
+ SUBROUTINE GZ_V_VW_PHY(D,PA,PDZZ,PGZ_V_VW)
+! #######################################################
+!
+!!**** *GZ_V_VW - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Z
+!! direction for a variable placed at the
+!! V point and the result is placed at
+!! the VW vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Z cartesian direction for a field PA placed at the
+! V point. The result is placed at the VW vorticity point.
+!
+!
+! dzm(PA)
+! PGZ_V_VW = ------
+! ____y
+! d*zz
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MYM : Shuman functions (mean operators)
+!! DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE SHUMAN_PHY, ONLY : DZM_PHY, MYM_PHY
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGZ_V_VW ! result UW point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: PA_WORK, PDZZ_WORK
+!
+INTEGER :: JI,JJ,JK
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GZ_V_VW
+! ---------------------
+!
+CALL DZM_PHY(D,PA,PA_WORK)
+CALL MYM_PHY(D,PDZZ,PDZZ_WORK)
+!
+!$mnh__expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
+PGZ_V_VW(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)= PA_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT) / PDZZ_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)
+!$mnh_end_expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
+!----------------------------------------------------------------------------
+!
+END SUBROUTINE GZ_V_VW_PHY
+END MODULE MODE_GRADIENT_V_PHY
diff --git a/src/common/turb/mode_tridiag_wind.F90 b/src/common/turb/mode_tridiag_wind.F90
index 27faf80c32e35d9d7ac54d09e76536887ac6544a..29c06958b4db42693bce486550481a4b3957b2d8 100644
--- a/src/common/turb/mode_tridiag_wind.F90
+++ b/src/common/turb/mode_tridiag_wind.F90
@@ -126,27 +126,28 @@ IMPLICIT NONE
!* 0.1 declarations of arguments
!
TYPE(DIMPHYEX_t), INTENT(IN) :: D
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM ! variable at t-1
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! upper diag. elements
-REAL, DIMENSION(:,:), INTENT(IN) :: PCOEFS ! implicit coeff for the
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PVARM ! variable at t-1
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PA ! upper diag. elements
+REAL, DIMENSION(D%NIT,D%NJT), INTENT(IN) :: PCOEFS ! implicit coeff for the
! surface flux
REAL, INTENT(IN) :: PTSTEP ! Double time step
REAL, INTENT(IN) :: PEXPL,PIMPL ! weights of the temporal scheme
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJA ! (dry rho)*J averaged
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSOURCE ! source term of PVAR
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRHODJA ! (dry rho)*J averaged
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PSOURCE ! source term of PVAR
!
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PVARP ! variable at t+1
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PVARP ! variable at t+1
!
!* 0.2 declarations of local variables
!
-REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3)) :: ZY ,ZGAM
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZY ,ZGAM
! RHS of the equation, 3D work array
-REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZBET
+REAL, DIMENSION(D%NIT,D%NJT) :: ZBET
! 2D work array
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
!
! ---------------------------------------------------------------------------
!
@@ -161,27 +162,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 / PRHODJA(:,:,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 / PRHODJA(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 / PRHODJA(:,:,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 / PRHODJA(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 / PRHODJA(:,:,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 / PRHODJA(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
@@ -192,53 +197,53 @@ IF ( PIMPL > 1.E-10 ) THEN
!
! going up
!
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZBET(:,:) = 1. - PIMPL * ( PA(:,:,IKB+D%NKL) / PRHODJA(:,:,IKB) &
- + PCOEFS(:,:) * PTSTEP ) ! bet = b(ikb)
- PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZBET(IIB:IIE,IJB:IJE) = 1. - PIMPL * ( PA(IIB:IIE,IJB:IJE,IKB+D%NKL) / PRHODJA(IIB:IIE,IJB:IJE,IKB) &
+ + PCOEFS(IIB:IIE,IJB:IJE) * PTSTEP ) ! bet = b(ikb)
+ 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) / PRHODJA(:,:,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) / PRHODJA(IIB:IIE,IJB:IJE,JK-D%NKL) / ZBET(IIB:IIE,IJB:IJE)
! gam(k) = c(k-1) / bet
- ZBET(:,:) = 1. - PIMPL * ( PA(:,:,JK) * (1. + ZGAM(:,:,JK)) &
- + PA(:,:,JK+D%NKL) &
- ) / PRHODJA(:,:,JK)
+ ZBET(IIB:IIE,IJB:IJE) = 1. - PIMPL * ( PA(IIB:IIE,IJB:IJE,JK) * (1. + ZGAM(IIB:IIE,IJB:IJE,JK)) &
+ + PA(IIB:IIE,IJB:IJE,JK+D%NKL) &
+ ) / PRHODJA(IIB:IIE,IJB:IJE,JK)
! bet = b(k) - a(k)* gam(k)
- PVARP(:,:,JK)= ( ZY(:,:,JK) - PIMPL * PA(:,:,JK) / PRHODJA(:,:,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) / PRHODJA(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) / PRHODJA(:,:,IKE-D%NKL) / ZBET(:,:)
+ ZGAM(IIB:IIE,IJB:IJE,IKE) = PIMPL * PA(IIB:IIE,IJB:IJE,IKE) / PRHODJA(IIB:IIE,IJB:IJE,IKE-D%NKL) / ZBET(IIB:IIE,IJB:IJE)
! gam(k) = c(k-1) / bet
- ZBET(:,:) = 1. - PIMPL * ( PA(:,:,IKE) * (1. + ZGAM(:,:,IKE)) &
- ) / PRHODJA(:,:,IKE)
+ ZBET(IIB:IIE,IJB:IJE) = 1. - PIMPL * ( PA(IIB:IIE,IJB:IJE,IKE) * (1. + ZGAM(IIB:IIE,IJB:IJE,IKE)) &
+ ) / PRHODJA(IIB:IIE,IJB:IJE,IKE)
! bet = b(k) - a(k)* gam(k)
- PVARP(:,:,IKE)= ( ZY(:,:,IKE) - PIMPL * PA(:,:,IKE) / PRHODJA(:,:,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) / PRHODJA(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
@@ -247,10 +252,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/mode_turb_ver_dyn_flux.F90 b/src/common/turb/mode_turb_ver_dyn_flux.F90
index 2dcaf968515fcdcf2d535be6850a53f229010b28..a737845d27b93ffcb1b08086fd72e084e74c76c6 100644
--- a/src/common/turb/mode_turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -216,6 +216,9 @@ USE MODD_OCEANH, ONLY: XSSUFL, XSSUFL_T,XSSVFL
USE MODD_PARAMETERS, ONLY: JPVEXT_TURB,XUNDEF
USE MODD_TURB_n, ONLY: TURB_t
!
+USE SHUMAN_PHY
+USE MODE_GRADIENT_U_PHY, ONLY : GZ_U_UW_PHY
+USE MODE_GRADIENT_V_PHY, ONLY : GZ_V_VW_PHY
!
USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
@@ -318,9 +321,10 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: &
ZSOURCE, & ! source of evolution for the treated variable
ZKEFF, & ! effectif diffusion coeff = LT * SQRT( TKE )
ZWORK1,ZWORK2,&
- ZWORK3,ZWORK4! working var. for shuman operators (array syntax)
+ ZWORK3,ZWORK4,&
+ ZWORK5,ZWORK6! working var. for shuman operators (array syntax)
!
-INTEGER :: IKB,IKE !
+INTEGER :: IIE,IIB,IJE,IJB,IKB,IKE ! index value for the mass points of the domain
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
INTEGER :: JSV,JI,JJ,JK ! scalar loop counter
@@ -339,38 +343,45 @@ TYPE(TFIELDDATA) :: TZFIELD
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',0,ZHOOK_HANDLE)
!
-ZA=XUNDEF
-PDP=XUNDEF
+ZA(:,:,:)=XUNDEF
+PDP(:,:,:)=XUNDEF
!
IKT=D%NKT
IKTB=D%NKTB
IKTE=D%NKTE
IKB=D%NKB
IKE=D%NKE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
!
-ZSOURCE = 0.
-ZFLXZ = 0.
+ZSOURCE(:,:,:) = 0.
+ZFLXZ(:,:,:) = 0.
ZCMFS = CSTURB%XCMFS
IF (OHARAT) ZCMFS=1.
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZDIRSINZW(IIB:IIE,IJB:IJE) = SQRT(1.-PDIRCOSZW(IIB:IIE,IJB:IJE)**2)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
! compute the coefficients for the uncentred gradient computation near the
! ground
!
! With OHARATU length scale and TKE are at half levels so remove MZM
!
IF (OHARAT) THEN
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
- ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-ELSE
- ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), D%NKA, D%NKU, D%NKL)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT) = PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)) + 50*MFMOIST(IIB:IIE,IJB:IJE,1:D%NKT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ELSE
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT))
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ CALL MZM_PHY(D,ZWORK1,ZKEFF)
ENDIF
!
-ZUSLOPEM(:,:,1)=PUSLOPEM(:,:)
-ZVSLOPEM(:,:,1)=PVSLOPEM(:,:)
+ZUSLOPEM(IIB:IIE,IJB:IJE,1)=PUSLOPEM(IIB:IIE,IJB:IJE)
+ZVSLOPEM(IIB:IIE,IJB:IJE,1)=PVSLOPEM(IIB:IIE,IJB:IJE)
!
!----------------------------------------------------------------------------
!
@@ -382,134 +393,136 @@ ZVSLOPEM(:,:,1)=PVSLOPEM(:,:)
!
! Preparation of the arguments for TRIDIAG_WIND
!
-ZWORK1 = MXM( ZKEFF )
-ZWORK2 = MXM( PDZZ )
-ZWORK3 = MXM(MZM(PRHODJ, D%NKA, D%NKU, D%NKL))
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-ZA(:,:,:) = -PTSTEP * ZCMFS * ZWORK1(:,:,:)* ZWORK3(:,:,:) / ZWORK2(:,:,:)**2
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+CALL MXM_PHY(D,ZKEFF,ZWORK1)
+CALL MXM_PHY(D,PDZZ,ZWORK2)
+CALL MZM_PHY(D,PRHODJ,ZWORK3)
+CALL MXM_PHY(D,ZWORK3,ZWORK4)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZA(IIB:IIE,IJB:IJE,1:D%NKT) = -PTSTEP * ZCMFS * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)* ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)**2
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
!
! Compute the source of U wind component
!
! compute the coefficient between the vertical flux and the 2 components of the
! wind following the slope
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) &
- * PCOSSLOPE(:,:)
-ZCOEFFLXV(:,:,1) = PCDUEFF(:,:) * PDIRCOSZW(:,:) * PSINSLOPE(:,:)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZCOEFFLXU(IIB:IIE,IJB:IJE,1) = PCDUEFF(IIB:IIE,IJB:IJE) * (PDIRCOSZW(IIB:IIE,IJB:IJE)**2 - ZDIRSINZW(IIB:IIE,IJB:IJE)**2) &
+ * PCOSSLOPE(IIB:IIE,IJB:IJE)
+ZCOEFFLXV(IIB:IIE,IJB:IJE,1) = PCDUEFF(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE) * PSINSLOPE(IIB:IIE,IJB:IJE)
! prepare the implicit scheme coefficients for the surface flux
-ZCOEFS(:,:,1)= ZCOEFFLXU(:,:,1) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) &
- +ZCOEFFLXV(:,:,1) * PSINSLOPE(:,:)
+ZCOEFS(IIB:IIE,IJB:IJE,1)= ZCOEFFLXU(IIB:IIE,IJB:IJE,1) * PCOSSLOPE(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE) &
+ +ZCOEFFLXV(IIB:IIE,IJB:IJE,1) * PSINSLOPE(IIB:IIE,IJB:IJE)
!
! average this flux to be located at the U,W vorticity point
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZCOEFS(:,:,1:1)=MXM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZWORK1(IIB:IIE,IJB:IJE,1:1)=ZCOEFS(IIB:IIE,IJB:IJE,1:1) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKB)
+CALL MXM_PHY(D,ZWORK1(IIB:IIE,IJB:IJE,1:1),ZCOEFS(IIB:IIE,IJB:IJE,1:1))
!
!
! ZSOURCE= FLUX /DZ
+CALL MXM_PHY(D,PRHODJ,ZWORK1)
IF (OOCEAN) THEN ! OCEAN MODEL ONLY
! Sfx flux assumed to be in SI & at vorticity point
- ZWORK1(:,:,D%NKU:D%NKU) = MXM(PRHODJ(:,:,D%NKU:D%NKU))
- ZWORK2(:,:,IKE:IKE) = MXM(PRHODJ(:,:,IKE:IKE))
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
IF (OCOUPLES) THEN
- ZSOURCE(:,:,IKE) = TURBN%XSSUFL_C(:,:,1)/PDZZ(:,:,IKE) &
- *0.5 * ( 1. + ZWORK1(:,:,D%NKU) / ZWORK2(:,:,IKE))
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = TURBN%XSSUFL_C(IIB:IIE,IJB:IJE,1)/PDZZ(IIB:IIE,IJB:IJE,IKE) &
+ *0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKU) / ZWORK1(IIB:IIE,IJB:IJE,IKE))
ELSE
- ZSOURCE(:,:,IKE) = XSSUFL(:,:)
- ZSOURCE(:,:,IKE) = ZSOURCE(:,:,IKE) /PDZZ(:,:,IKE) &
- *0.5 * ( 1. + ZWORK1(:,:,D%NKU) / ZWORK2(:,:,IKE) )
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = XSSUFL(IIB:IIE,IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = ZSOURCE(IIB:IIE,IJB:IJE,IKE) /PDZZ(IIB:IIE,IJB:IJE,IKE) &
+ *0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKU) / ZWORK1(IIB:IIE,IJB:IJE,IKE) )
ENDIF
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!No flux at the ocean domain bottom
- ZSOURCE(:,:,IKB) = 0.
- ZSOURCE(:,:,IKTB+1:IKTE-1) = 0
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = 0.
+ ZSOURCE(IIB:IIE,IJB:IJE,IKTB+1:IKTE-1) = 0
!
ELSE !ATMOS MODEL ONLY
- ZWORK1(:,:,D%NKA:D%NKA) = MXM(PRHODJ(:,:,D%NKA:D%NKA))
- ZWORK2(:,:,IKB:IKB) = MXM(PRHODJ(:,:,IKB:IKB))
IF (OCOUPLES) THEN
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZSOURCE(:,:,IKB) = TURBN%XSSUFL_C(:,:,1)/PDZZ(:,:,IKB) &
- * 0.5 * ( 1. + ZWORK1(:,:,D%NKA) / ZWORK2(:,:,IKB) )
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = TURBN%XSSUFL_C(IIB:IIE,IJB:IJE,1)/PDZZ(IIB:IIE,IJB:IJE,IKB) &
+ * 0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKA) / ZWORK1(IIB:IIE,IJB:IJE,IKB) )
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
ELSE
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
! compute the explicit tangential flux at the W point
- ZSOURCE(:,:,IKB) = &
- PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
- -PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) &
- -PTAU33M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = &
+ PTAU11M(IIB:IIE,IJB:IJE) * PCOSSLOPE(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE) * ZDIRSINZW(IIB:IIE,IJB:IJE) &
+ -PTAU12M(IIB:IIE,IJB:IJE) * PSINSLOPE(IIB:IIE,IJB:IJE) * ZDIRSINZW(IIB:IIE,IJB:IJE) &
+ -PTAU33M(IIB:IIE,IJB:IJE) * PCOSSLOPE(IIB:IIE,IJB:IJE) * ZDIRSINZW(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
! add the vertical part or the surface flux at the U,W vorticity point
!
- ZWORK3(:,:,IKB:IKB) = MXM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) )
- ZWORK4(:,:,IKB:IKB) = MXM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
- *ZUSLOPEM(:,:,1:1) &
- -ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
- *ZVSLOPEM(:,:,1:1) )
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZSOURCE(:,:,IKB) = &
- ( ZWORK3(:,:,IKB) &
- + ZWORK4(:,:,IKB) &
- - ZCOEFS(:,:,1) * PUM(:,:,IKB) * PIMPL &
- ) * 0.5 * ( 1. + ZWORK1(:,:,D%NKA) / ZWORK2(:,:,IKB) )
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKB)
+ ZWORK3(IIB:IIE,IJB:IJE,IKB:IKB) = ZSOURCE(IIB:IIE,IJB:IJE,IKB:IKB)/PDZZ(IIB:IIE,IJB:IJE,IKB:IKB)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKB)
+ CALL MXM_PHY(D,ZWORK3,ZWORK4)
+ ZWORK5(IIB:IIE,IJB:IJE,IKB:IKB)= ZCOEFFLXU(IIB:IIE,IJB:IJE,1:1) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKB) &
+ *ZUSLOPEM(IIB:IIE,IJB:IJE,1:1) &
+ -ZCOEFFLXV(IIB:IIE,IJB:IJE,1:1) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKB) &
+ *ZVSLOPEM(IIB:IIE,IJB:IJE,1:1)
+ CALL MXM_PHY(D,ZWORK5,ZWORK6)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = &
+ ( ZWORK4(IIB:IIE,IJB:IJE,IKB) &
+ + ZWORK6(IIB:IIE,IJB:IJE,IKB) &
+ - ZCOEFS(IIB:IIE,IJB:IJE,1) * PUM(IIB:IIE,IJB:IJE,IKB) * PIMPL &
+ ) * 0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKA) / ZWORK1(IIB:IIE,IJB:IJE,IKB) )
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
ENDIF
!
- ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
- ZSOURCE(:,:,IKE) = 0.
+ ZSOURCE(IIB:IIE,IJB:IJE,IKTB+1:IKTE-1) = 0.
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = 0.
ENDIF !end ocean or atmosphere cases
!
! Obtention of the split U at t+ deltat
!
CALL TRIDIAG_WIND(D,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
- MXM(PRHODJ),ZSOURCE,ZRES)
+ ZWORK1,ZSOURCE,ZRES)
!
! Compute the equivalent tendency for the U wind component
!
-ZWORK1 = MXM(PRHODJ)
-ZWORK2 = MXM(ZKEFF)
-ZWORK3 = DZM(PIMPL*ZRES + PEXPL*PUM, D%NKA, D%NKU, D%NKL)
-ZWORK4 = MXM(PDZZ)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PRUS(:,:,:)=PRUS(:,:,:)+ZWORK1(:,:,:)*(ZRES(:,:,:)-PUM(:,:,:))/PTSTEP
+CALL MXM_PHY(D,PRHODJ,ZWORK1)
+CALL MXM_PHY(D,ZKEFF,ZWORK2)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)=PIMPL*ZRES(IIB:IIE,IJB:IJE,1:D%NKT) + PEXPL*PUM(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+CALL DZM_PHY(D,ZWORK3,ZWORK4)
+CALL MXM_PHY(D,PDZZ,ZWORK5)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+PRUS(IIB:IIE,IJB:IJE,1:D%NKT)=PRUS(IIB:IIE,IJB:IJE,1:D%NKT)+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*(ZRES(IIB:IIE,IJB:IJE,1:D%NKT)-PUM(IIB:IIE,IJB:IJE,1:D%NKT))/PTSTEP
!
!
!* 5.2 Partial Dynamic Production
!
! vertical flux of the U wind component
!
-ZFLXZ(:,:,:) = -ZCMFS * ZWORK2(:,:,:) * ZWORK3(:,:,:) / ZWORK4(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = -ZCMFS * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
! surface flux
-ZWORK1(:,:,IKB:IKB) = MXM(PDZZ(:,:,IKB:IKB))
-ZWORK2(:,:,D%NKA:D%NKA) = MXM(PRHODJ(:,:,D%NKA:D%NKA))
-ZWORK3(:,:,IKB:IKB) = MXM(PRHODJ(:,:,IKB:IKB))
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZFLXZ(:,:,IKB) = ZWORK1(:,:,IKB) * &
- ( ZSOURCE(:,:,IKB) &
- +ZCOEFS(:,:,1) * ZRES(:,:,IKB) * PIMPL &
- ) / 0.5 / ( 1. + ZWORK2(:,:,D%NKA)/ ZWORK3(:,:,IKB) )
-!
-ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+CALL MXM_PHY(D,PDZZ,ZWORK1)
+CALL MXM_PHY(D,PRHODJ,ZWORK2)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZFLXZ(IIB:IIE,IJB:IJE,IKB) = ZWORK1(IIB:IIE,IJB:IJE,IKB) * &
+ ( ZSOURCE(IIB:IIE,IJB:IJE,IKB) &
+ +ZCOEFS(IIB:IIE,IJB:IJE,1) * ZRES(IIB:IIE,IJB:IJE,IKB) * PIMPL &
+ ) / 0.5 / ( 1. + ZWORK2(IIB:IIE,IJB:IJE,D%NKA)/ ZWORK2(IIB:IIE,IJB:IJE,IKB) )
+!
+ZFLXZ(IIB:IIE,IJB:IJE,D%NKA) = ZFLXZ(IIB:IIE,IJB:IJE,IKB)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
IF (OOCEAN) THEN !ocean model at phys sfc (ocean domain top)
- ZWORK1(:,:,IKE:IKE) = MXM(PDZZ(:,:,IKE:IKE))
- ZWORK2(:,:,D%NKU:D%NKU) = MXM(PRHODJ(:,:,D%NKU:D%NKU))
- ZWORK3(:,:,IKE:IKE) = MXM(PRHODJ(:,:,IKE:IKE))
-!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZFLXZ(:,:,IKE) = ZWORK1(:,:,IKE) * &
- ZSOURCE(:,:,IKE) &
- / 0.5 / ( 1. + ZWORK2(:,:,D%NKU)/ ZWORK3(:,:,IKE) )
- ZFLXZ(:,:,D%NKU) = ZFLXZ(:,:,IKE)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+!
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZFLXZ(IIB:IIE,IJB:IJE,IKE) = ZWORK1(IIB:IIE,IJB:IJE,IKE) * &
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) &
+ / 0.5 / ( 1. + ZWORK2(IIB:IIE,IJB:IJE,D%NKU)/ ZWORK2(IIB:IIE,IJB:IJE,IKE) )
+ ZFLXZ(IIB:IIE,IJB:IJE,D%NKU) = ZFLXZ(IIB:IIE,IJB:IJE,IKE)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END IF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -529,25 +542,38 @@ END IF
!
! first part of total momentum flux
!
-PWU(:,:,:) = ZFLXZ(:,:,:)
+PWU(IIB:IIE,IJB:IJE,1:D%NKT) = ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)
!
! Contribution to the dynamic production of TKE
! compute the dynamic production at the mass point
!
-PDP(:,:,:) = - MZF(MXF(ZFLXZ * GZ_U_UW(PUM,PDZZ, D%NKA, D%NKU, D%NKL)), D%NKA, D%NKU, D%NKL)
+CALL GZ_U_UW_PHY(D,PUM,PDZZ,ZWORK1)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) = ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+CALL MXF_PHY(D,ZWORK2,ZWORK3)
+CALL MZF_PHY(D,ZWORK3,ZWORK4)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+PDP(IIB:IIE,IJB:IJE,1:D%NKT) = -ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
-PDP(:,:,IKB:IKB) = - MXF ( &
- ZFLXZ(:,:,IKB+D%NKL:IKB+D%NKL) * (PUM(:,:,IKB+D%NKL:IKB+D%NKL)-PUM(:,:,IKB:IKB)) &
- / MXM(PDZZ(:,:,IKB+D%NKL:IKB+D%NKL)) &
- )
+CALL MXM_PHY(D,PDZZ,ZWORK1)
+ZWORK2(IIB:IIE,IJB:IJE,IKB) = ZFLXZ(IIB:IIE,IJB:IJE,IKB+D%NKL) * (PUM(IIB:IIE,IJB:IJE,IKB+D%NKL)-PUM(IIB:IIE,IJB:IJE,IKB)) &
+ / ZWORK1(IIB:IIE,IJB:IJE,IKB+D%NKL)
+CALL MXF_PHY(D,ZWORK2,ZWORK3)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+PDP(IIB:IIE,IJB:IJE,IKB) = -ZWORK3(IIB:IIE,IJB:IJE,IKB)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
IF (OOCEAN) THEN
+ ZWORK2(IIB:IIE,IJB:IJE,IKE) = ZFLXZ(IIB:IIE,IJB:IJE,IKE-D%NKL) * (PUM(IIB:IIE,IJB:IJE,IKE)-PUM(IIB:IIE,IJB:IJE,IKE-D%NKL)) &
+ / ZWORK1(IIB:IIE,IJB:IJE,IKE-D%NKL)
+ CALL MXF_PHY(D,ZWORK2,ZWORK3)
! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
- PDP(:,:,IKE:IKE) = - MXF ( &
- ZFLXZ(:,:,IKE-D%NKL:IKE-D%NKL) * (PUM(:,:,IKE:IKE)-PUM(:,:,IKE-D%NKL:IKE-D%NKL)) &
- / MXM(PDZZ(:,:,IKE-D%NKL:IKE-D%NKL)) &
- )
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ PDP(IIB:IIE,IJB:IJE,IKE) = -ZWORK3(IIB:IIE,IJB:IJE,IKE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END IF
!
! Storage in the LES configuration
@@ -661,126 +687,131 @@ END IF
!
! Preparation of the arguments for TRIDIAG_WIND
!!
-ZWORK1 = MYM( ZKEFF )
-ZWORK2 = MYM(MZM(PRHODJ, D%NKA, D%NKU, D%NKL))
-ZWORK3 = MYM( PDZZ )
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-ZA(:,:,:) = - PTSTEP * ZCMFS * ZWORK1(:,:,:) * ZWORK2(:,:,:) / &
- ZWORK3(:,:,:)**2
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+CALL MYM_PHY(D,ZKEFF,ZWORK1)
+CALL MYM_PHY(D,PDZZ,ZWORK2)
+CALL MZM_PHY(D,PRHODJ,ZWORK3)
+CALL MYM_PHY(D,ZWORK3,ZWORK4)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZA(IIB:IIE,IJB:IJE,1:D%NKT) = -PTSTEP * ZCMFS * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)* ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)**2
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
!
!
! Compute the source of V wind component
! compute the coefficient between the vertical flux and the 2 components of the
! wind following the slope
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) &
- * PSINSLOPE(:,:)
-ZCOEFFLXV(:,:,1) = PCDUEFF(:,:) * PDIRCOSZW(:,:) * PCOSSLOPE(:,:)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZCOEFFLXU(IIB:IIE,IJB:IJE,1) = PCDUEFF(IIB:IIE,IJB:IJE) * (PDIRCOSZW(IIB:IIE,IJB:IJE)**2 - ZDIRSINZW(IIB:IIE,IJB:IJE)**2) &
+ * PSINSLOPE(IIB:IIE,IJB:IJE)
+ZCOEFFLXV(IIB:IIE,IJB:IJE,1) = PCDUEFF(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE) * PCOSSLOPE(IIB:IIE,IJB:IJE)
! prepare the implicit scheme coefficients for the surface flux
-ZCOEFS(:,:,1)= ZCOEFFLXU(:,:,1) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) &
- +ZCOEFFLXV(:,:,1) * PCOSSLOPE(:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ZCOEFS(IIB:IIE,IJB:IJE,1)= ZCOEFFLXU(IIB:IIE,IJB:IJE,1) * PSINSLOPE(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE) &
+ +ZCOEFFLXV(IIB:IIE,IJB:IJE,1) * PCOSSLOPE(IIB:IIE,IJB:IJE)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
! average this flux to be located at the V,W vorticity point
-ZCOEFS(:,:,1:1)=MYM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
+ZWORK1(IIB:IIE,IJB:IJE,1:1)=ZCOEFS(IIB:IIE,IJB:IJE,1:1) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKB)
+CALL MYM_PHY(D,ZWORK1(IIB:IIE,IJB:IJE,1:1),ZCOEFS(IIB:IIE,IJB:IJE,1:1))
!
+CALL MYM_PHY(D,PRHODJ,ZWORK1) ! it was MXM(PRHODJ) : bug corrected now ? REPRO55 OCEAN
IF (OOCEAN) THEN ! Ocean case
- ZWORK1(:,:,D%NKU:D%NKU) = MXM(PRHODJ(:,:,D%NKU:D%NKU))
- ZWORK2(:,:,IKE:IKE) = MXM(PRHODJ(:,:,IKE:IKE))
IF (OCOUPLES) THEN
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZSOURCE(:,:,IKE) = TURBN%XSSVFL_C(:,:,1)/PDZZ(:,:,IKE) &
- *0.5 * ( 1. + ZWORK1(:,:,D%NKU) / ZWORK2(:,:,IKE))
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = TURBN%XSSVFL_C(IIB:IIE,IJB:IJE,1)/PDZZ(IIB:IIE,IJB:IJE,IKE) &
+ *0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKU) / ZWORK1(IIB:IIE,IJB:IJE,IKE))
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
ELSE
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZSOURCE(:,:,IKE) = XSSVFL(:,:)
- ZSOURCE(:,:,IKE) = ZSOURCE(:,:,IKE)/PDZZ(:,:,IKE) &
- *0.5 * ( 1. + ZWORK1(:,:,D%NKU) / ZWORK2(:,:,IKE))
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = XSSVFL(IIB:IIE,IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = ZSOURCE(IIB:IIE,IJB:IJE,IKE)/PDZZ(IIB:IIE,IJB:IJE,IKE) &
+ *0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKU) / ZWORK1(IIB:IIE,IJB:IJE,IKE))
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END IF
!No flux at the ocean domain bottom
- ZSOURCE(:,:,IKB) = 0.
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = 0.
ELSE ! Atmos case
- ZWORK1(:,:,D%NKA:D%NKA) = MYM(PRHODJ(:,:,D%NKA:D%NKA))
- ZWORK2(:,:,IKB:IKB) = MYM(PRHODJ(:,:,IKB:IKB))
- ZWORK3(:,:,IKB:IKB) = MYM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
- *ZUSLOPEM(:,:,1:1) &
- +ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
- *ZVSLOPEM(:,:,1:1) )
+ ZWORK3(IIB:IIE,IJB:IJE,IKB:IKB) = ZCOEFFLXU(IIB:IIE,IJB:IJE,1:1) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKB) &
+ *ZUSLOPEM(IIB:IIE,IJB:IJE,1:1) &
+ +ZCOEFFLXV(IIB:IIE,IJB:IJE,1:1) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKB) &
+ *ZVSLOPEM(IIB:IIE,IJB:IJE,1:1)
+ CALL MYM_PHY(D,ZWORK3,ZWORK6)
!
IF (.NOT.OCOUPLES) THEN ! only atmosp without coupling
! compute the explicit tangential flux at the W point
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZSOURCE(:,:,IKB) = &
- PTAU11M(:,:) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
- +PTAU12M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) &
- -PTAU33M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZWORK4(:,:,IKB:IKB) = MYM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) )
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = &
+ PTAU11M(IIB:IIE,IJB:IJE) * PSINSLOPE(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE) * ZDIRSINZW(IIB:IIE,IJB:IJE) &
+ +PTAU12M(IIB:IIE,IJB:IJE) * PCOSSLOPE(IIB:IIE,IJB:IJE) * ZDIRSINZW(IIB:IIE,IJB:IJE) &
+ -PTAU33M(IIB:IIE,IJB:IJE) * PSINSLOPE(IIB:IIE,IJB:IJE) * ZDIRSINZW(IIB:IIE,IJB:IJE) * PDIRCOSZW(IIB:IIE,IJB:IJE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKB)
+ ZWORK3(IIB:IIE,IJB:IJE,IKB:IKB) = ZSOURCE(IIB:IIE,IJB:IJE,IKB:IKB)/PDZZ(IIB:IIE,IJB:IJE,IKB:IKB)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKB:IKB)
+ CALL MYM_PHY(D,ZWORK3,ZWORK5)
!
! add the vertical part or the surface flux at the V,W vorticity point
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZSOURCE(:,:,IKB) = &
- ( ZWORK4(:,:,IKB) &
- + ZWORK3(:,:,IKB) &
- - ZCOEFS(:,:,1) * PVM(:,:,IKB) * PIMPL &
- ) * 0.5 * ( 1. + ZWORK1(:,:,D%NKA) / ZWORK2(:,:,IKB) )
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = &
+ ( ZWORK5(IIB:IIE,IJB:IJE,IKB) &
+ + ZWORK6(IIB:IIE,IJB:IJE,IKB) &
+ - ZCOEFS(IIB:IIE,IJB:IJE,1) * PVM(IIB:IIE,IJB:IJE,IKB) * PIMPL &
+ ) * 0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKA) / ZWORK1(IIB:IIE,IJB:IJE,IKB) )
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
ELSE !atmosphere when coupling
! input flux assumed to be in SI and at vorticity point
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZSOURCE(:,:,IKB) = -TURBN%XSSVFL_C(:,:,1)/(1.*PDZZ(:,:,IKB)) &
- * 0.5 * ( 1. + ZWORK1(:,:,D%NKA) / ZWORK2(:,:,IKB) )
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZSOURCE(IIB:IIE,IJB:IJE,IKB) = -TURBN%XSSVFL_C(IIB:IIE,IJB:IJE,1)/(1.*PDZZ(IIB:IIE,IJB:IJE,IKB)) &
+ * 0.5 * ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKA) / ZWORK1(IIB:IIE,IJB:IJE,IKB) )
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
ENDIF
!No flux at the atmosphere top
- ZSOURCE(:,:,IKE) = 0.
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) = 0.
ENDIF ! End of Ocean or Atmospher Cases
-ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
+ZSOURCE(IIB:IIE,IJB:IJE,IKTB+1:IKTE-1) = 0.
!
! Obtention of the split V at t+ deltat
CALL TRIDIAG_WIND(D,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
- MYM(PRHODJ),ZSOURCE,ZRES)
+ ZWORK1,ZSOURCE,ZRES)
!
! Compute the equivalent tendency for the V wind component
!
-ZWORK1 = MYM(PRHODJ(:,:,:))
-ZWORK2 = MYM(ZKEFF)
-ZWORK3 = DZM(PIMPL*ZRES + PEXPL*PVM, D%NKA, D%NKU, D%NKL)
-ZWORK4 = MYM(PDZZ)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PRVS(:,:,:)=PRVS(:,:,:)+ZWORK1(:,:,:)*(ZRES(:,:,:)-PVM(:,:,:))/PTSTEP
+CALL MYM_PHY(D,PRHODJ,ZWORK1)
+CALL MYM_PHY(D,ZKEFF,ZWORK2)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)=PIMPL*ZRES(IIB:IIE,IJB:IJE,1:D%NKT) + PEXPL*PUM(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+CALL DZM_PHY(D,ZWORK3,ZWORK4)
+CALL MYM_PHY(D,PDZZ,ZWORK5)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+PRVS(IIB:IIE,IJB:IJE,1:D%NKT)=PRVS(IIB:IIE,IJB:IJE,1:D%NKT)+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*(ZRES(IIB:IIE,IJB:IJE,1:D%NKT)-PVM(IIB:IIE,IJB:IJE,1:D%NKT))/PTSTEP
!
!
!* 6.2 Complete 1D dynamic Production
!
! vertical flux of the V wind component
!
-ZFLXZ(:,:,:) = -ZCMFS * ZWORK2(:,:,:) * ZWORK3(:,:,:) / ZWORK4(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = -ZCMFS * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-ZFLXZ(:,:,IKB) = ZWORK4(:,:,IKB) * &
- ( ZSOURCE(:,:,IKB) &
- +ZCOEFS(:,:,1) * ZRES(:,:,IKB) * PIMPL &
- ) / 0.5 / ( 1. + ZWORK1(:,:,D%NKA) / ZWORK1(:,:,IKB) )
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZFLXZ(IIB:IIE,IJB:IJE,IKB) = ZWORK5(IIB:IIE,IJB:IJE,IKB) * &
+ ( ZSOURCE(IIB:IIE,IJB:IJE,IKB) &
+ +ZCOEFS(IIB:IIE,IJB:IJE,1) * ZRES(IIB:IIE,IJB:IJE,IKB) * PIMPL &
+ ) / 0.5 / ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKA) / ZWORK1(IIB:IIE,IJB:IJE,IKB) )
!
!
-ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ZFLXZ(IIB:IIE,IJB:IJE,D%NKA) = ZFLXZ(IIB:IIE,IJB:IJE,IKB)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
IF (OOCEAN) THEN
- !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
- ZFLXZ(:,:,IKE) = ZWORK4(:,:,IKE) * &
- ZSOURCE(:,:,IKE) &
- / 0.5 / ( 1. + ZWORK1(:,:,D%NKU) / ZWORK1(:,:,IKE) )
- ZFLXZ(:,:,D%NKU) = ZFLXZ(:,:,IKE)
- !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZFLXZ(IIB:IIE,IJB:IJE,IKE) = ZWORK5(IIB:IIE,IJB:IJE,IKE) * &
+ ZSOURCE(IIB:IIE,IJB:IJE,IKE) &
+ / 0.5 / ( 1. + ZWORK1(IIB:IIE,IJB:IJE,D%NKU) / ZWORK1(IIB:IIE,IJB:IJE,IKE) )
+ ZFLXZ(IIB:IIE,IJB:IJE,D%NKU) = ZFLXZ(IIB:IIE,IJB:IJE,IKE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END IF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -800,31 +831,43 @@ END IF
!
! second part of total momentum flux
!
-PWV(:,:,:) = ZFLXZ(:,:,:)
+PWV(IIB:IIE,IJB:IJE,1:D%NKT) = ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)
!
! Contribution to the dynamic production of TKE
! compute the dynamic production contribution at the mass point
!
-ZA(:,:,:) = - MZF(MYF(ZFLXZ * GZ_V_VW(PVM,PDZZ, D%NKA, D%NKU, D%NKL)), D%NKA, D%NKU, D%NKL)
+CALL GZ_V_VW_PHY(D,PVM,PDZZ,ZWORK1)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) = ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+CALL MYF_PHY(D,ZWORK2,ZWORK3)
+CALL MZF_PHY(D,ZWORK3,ZWORK4)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZA(IIB:IIE,IJB:IJE,1:D%NKT) = -ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
-ZA(:,:,IKB:IKB) = &
- - MYF ( &
-ZFLXZ(:,:,IKB+D%NKL:IKB+D%NKL) * (PVM(:,:,IKB+D%NKL:IKB+D%NKL)-PVM(:,:,IKB:IKB)) &
- / MYM(PDZZ(:,:,IKB+D%NKL:IKB+D%NKL)) &
- )
+CALL MYM_PHY(D,PDZZ,ZWORK1)
+ZWORK2(IIB:IIE,IJB:IJE,IKB) = ZFLXZ(IIB:IIE,IJB:IJE,IKB+D%NKL) * (PVM(IIB:IIE,IJB:IJE,IKB+D%NKL)-PVM(IIB:IIE,IJB:IJE,IKB)) &
+ / ZWORK1(IIB:IIE,IJB:IJE,IKB+D%NKL)
+CALL MYF_PHY(D,ZWORK2,ZWORK3)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ZA(IIB:IIE,IJB:IJE,IKB) = -ZWORK3(IIB:IIE,IJB:IJE,IKB)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
!
IF (OOCEAN) THEN
! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
- ZA(:,:,IKE:IKE) = - MYF ( &
- ZFLXZ(:,:,IKE-D%NKL:IKE-D%NKL) * (PVM(:,:,IKE:IKE)-PVM(:,:,IKE-D%NKL:IKE-D%NKL)) &
- / MYM(PDZZ(:,:,IKE-D%NKL:IKE-D%NKL)) &
- )
+ ZWORK2(IIB:IIE,IJB:IJE,IKE) = ZFLXZ(IIB:IIE,IJB:IJE,IKE-D%NKL) * (PVM(IIB:IIE,IJB:IJE,IKE)-PVM(IIB:IIE,IJB:IJE,IKE-D%NKL)) &
+ / ZWORK1(IIB:IIE,IJB:IJE,IKE-D%NKL)
+ CALL MYF_PHY(D,ZWORK2,ZWORK3)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+ ZA(IIB:IIE,IJB:IJE,IKE) = -ZWORK3(IIB:IIE,IJB:IJE,IKE)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
END IF
!
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
-!$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=1:D%NKT)
+PDP(IIB:IIE,IJB:IJE,1:D%NKT)=PDP(IIB:IIE,IJB:IJE,1:D%NKT)+ZA(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
! Storage in the LES configuration
!
diff --git a/src/mesonh/aux/shuman_phy.f90 b/src/mesonh/aux/shuman_phy.f90
index 28feda3c7fe85a89e6740c439493eafd4faf54cb..69adba5459df1ae98fc69c8c0e0154339c9094bd 100644
--- a/src/mesonh/aux/shuman_phy.f90
+++ b/src/mesonh/aux/shuman_phy.f90
@@ -1,6 +1,110 @@
MODULE SHUMAN_PHY
IMPLICIT NONE
CONTAINS
+! ###############################
+ SUBROUTINE MXM_PHY(D,PA,PMXM)
+! ###############################
+!
+!!**** *MXM* - Shuman operator : mean operator in x direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the x direction (I index) for a field PA localized at a mass
+! point. The result is localized at a x-flux point (u point).
+!
+!!** METHOD
+!! ------
+!! The result PMXM(i,:,:) is defined by 0.5*(PA(i,:,:)+PA(i-1,:,:))
+!! At i=1, PMXM(1,:,:) are replaced by the values of PMXM,
+!! which are the right values in the x-cyclic case.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!! optimisation 20/08/00 J. Escobar
+!! correction of in halo/pseudo-cyclic calculation for JPHEXT<> 1
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+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) :: PMXM ! result at flux localization
+
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JI,JJ,JK ! Loop index in x direction
+INTEGER :: IIU ! Size of the array in the x direction
+!
+INTEGER :: JJK,IJU,IKU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MXM
+! ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+JIJKOR = 1 + 1
+JIJKEND = IIU*IJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+ PMXM(JIJK,1,1) = 0.5*( PA(JIJK,1,1)+PA(JIJK-1,1,1) )
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JI=1,JPHEXT
+ DO JJK=1,IJU*IKU
+ PMXM(JI,JJK,1) = PMXM(IIU-2*JPHEXT+JI,JJK,1) ! for reprod JPHEXT <> 1
+ END DO
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MXM_PHY
! ###############################
SUBROUTINE MZM_PHY(D,PA,PMZM)
! ###############################
@@ -96,4 +200,497 @@ END DO
!-------------------------------------------------------------------------------
!
END SUBROUTINE MZM_PHY
+! ###############################
+ SUBROUTINE MYM_PHY(D,PA,PMYM)
+! ###############################
+!
+!!**** *MYM* - Shuman operator : mean operator in y direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the y direction (J index) for a field PA localized at a mass
+! point. The result is localized at a y-flux point (v point).
+!
+!!** METHOD
+!! ------
+!! The result PMYM(:,j,:) is defined by 0.5*(PA(:,j,:)+PA(:,j-1,:))
+!! At j=1, PMYM(:,j,:) are replaced by the values of PMYM,
+!! which are the right values in the y-cyclic case.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!! optimisation 20/08/00 J. Escobar
+!! correction of in halo/pseudo-cyclic calculation for JPHEXT<> 1
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+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) :: PMYM ! result at flux localization
+
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JJ ! Loop index in y direction
+INTEGER :: IJU ! Size of the array in the y direction
+!
+!
+INTEGER :: IIU,IKU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MYM
+! ------------------
+!
+IIU=SIZE(PA,1)
+IJU=SIZE(PA,2)
+IKU=SIZE(PA,3)
+!
+JIJKOR = 1 + IIU
+JIJKEND = IIU*IJU*IKU
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+ PMYM(JIJK,1,1) = 0.5*( PA(JIJK,1,1)+PA(JIJK-IIU,1,1) )
+END DO
+!
+DO JJ=1,JPHEXT
+ PMYM(:,JJ,:) = PMYM(:,IJU-2*JPHEXT+JJ,:) ! for reprod JPHEXT <> 1
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MYM_PHY
+! ###############################
+ SUBROUTINE DZM_PHY(D,PA,PDZM)
+! ###############################
+!
+!!**** *DZM* - Shuman operator : finite difference operator in z direction
+!! for a variable at a mass localization
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a finite difference
+! 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 PDZM(:,j,:) is defined by (PA(:,:,k)-PA(:,:,k-1))
+!! At k=1, PDZM(:,:,k) 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 05/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) :: PDZM ! result at flux
+ ! side
+!
+!* 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
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF DZM
+! ------------------
+!
+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
+ PDZM(JIJK,1,1) = PA(JIJK,1,1)-PA(JIJK-IIU*IJU,1,1)
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJ=1,IIU*IJU
+ PDZM(JIJ,1,1) = -999.
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE DZM_PHY
+! ###############################
+ SUBROUTINE MZF_PHY(D,PA,PMZF)
+! ###############################
+!
+!!**** *MZF* - Shuman operator : mean operator in z direction for a
+!! variable at a flux side
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the z direction (K index) for a field PA localized at a z-flux
+! point (w point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PMZF(:,:,k) is defined by 0.5*(PA(:,:,k)+PA(:,:,k+1))
+!! At k=size(PA,3), PMZF(:,:,k) 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 flux localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PMZF ! result at mass 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
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MZF
+! ------------------
+!
+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
+ PMZF(JIJK-IIU*IJU,1,1) = 0.5*( PA(JIJK-IIU*IJU,1,1)+PA(JIJK,1,1) )
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJ=1,IIU*IJU
+ PMZF(JIJ,1,IKU) = PMZF(JIJ,1,IKU-1) !-999.
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MZF_PHY
+! ###############################
+ SUBROUTINE MXF_PHY(D,PA,PMXF)
+! ###############################
+!
+!!**** *MXF* - Shuman operator : mean operator in x direction for a
+!! variable at a flux side
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the x direction (I index) for a field PA localized at a x-flux
+! point (u point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PMXF(i,:,:) is defined by 0.5*(PA(i,:,:)+PA(i+1,:,:))
+!! At i=size(PA,1), PMXF(i,:,:) are replaced by the values of PMXF,
+!! which are the right values in the x-cyclic case
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!! optimisation 20/08/00 J. Escobar
+!! correction of in halo/pseudo-cyclic calculation for JPHEXT<> 1
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+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 flux localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PMXF ! result at mass localization
+
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JI ! Loop index in x direction
+INTEGER :: IIU ! upper bound in x direction of PA
+!
+INTEGER :: JJK,IJU,IKU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MXF
+! ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+JIJKOR = 1 + 1
+JIJKEND = IIU*IJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+ PMXF(JIJK-1,1,1) = 0.5*( PA(JIJK-1,1,1)+PA(JIJK,1,1) )
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JI=1,JPHEXT
+ DO JJK=1,IJU*IKU
+ PMXF(IIU-JPHEXT+JI,JJK,1) = PMXF(JPHEXT+JI,JJK,1) ! for reprod JPHEXT <> 1
+ END DO
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MXF_PHY
+! ###############################
+ SUBROUTINE MYF_PHY(D,PA,PMYF)
+! ###############################
+!
+!!**** *MYF* - Shuman operator : mean operator in y direction for a
+!! variable at a flux side
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the y direction (J index) for a field PA localized at a y-flux
+! point (v point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PMYF(i,:,:) is defined by 0.5*(PA(:,j,:)+PA(:,j+1,:))
+!! At j=size(PA,2), PMYF(:,j,:) are replaced by the values of PMYF,
+!! which are the right values in the y-cyclic case
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT: define the number of marginal points out of the
+!! physical domain along the horizontal directions.
+!!
+!! 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
+!! Modification to include the periodic case 13/10/94 J.Stein
+!! optimisation 20/08/00 J. Escobar
+!! correction of in halo/pseudo-cyclic calculation for JPHEXT<> 1
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+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 flux localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PMYF ! result at mass localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JJ ! Loop index in y direction
+INTEGER :: IJU ! upper bound in y direction of PA
+!
+INTEGER :: IIU,IKU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MYF
+! ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+JIJKOR = 1 + IIU
+JIJKEND = IIU*IJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+ PMYF(JIJK-IIU,1,1) = 0.5*( PA(JIJK-IIU,1,1)+PA(JIJK,1,1) )
+END DO
+!
+DO JJ=1,JPHEXT
+ PMYF(:,IJU-JPHEXT+JJ,:) = PMYF(:,JPHEXT+JJ,:) ! for reprod JPHEXT <> 1
+END DO
+!
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MYF_PHY
END MODULE SHUMAN_PHY