From cad984553163e5b21a26455885b4d8fa056b4b99 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Fri, 12 Aug 2022 12:08:37 +0200
Subject: [PATCH] Quentin 12/08/2022: Complete cross gradients and shuman as
subroutines with explicit dimensions declarations
---
src/common/aux/gradient_m_phy.F90 | 288 ++++++++++++++++++++++++++++++
src/common/aux/gradient_u_phy.F90 | 142 ++++++++++++++-
src/common/aux/gradient_v_phy.F90 | 137 +++++++++++++-
src/common/aux/gradient_w_phy.F90 | 229 ++++++++++++++++++++++++
src/common/aux/shuman_phy.F90 | 192 +++++++++++++++++++-
src/mesonh/aux/shuman_phy.f90 | 206 +++++++++++++++++++++
6 files changed, 1180 insertions(+), 14 deletions(-)
diff --git a/src/common/aux/gradient_m_phy.F90 b/src/common/aux/gradient_m_phy.F90
index 94dd54aea..ed8d0d2ad 100644
--- a/src/common/aux/gradient_m_phy.F90
+++ b/src/common/aux/gradient_m_phy.F90
@@ -346,4 +346,292 @@ ENDIF
!
IF (LHOOK) CALL DR_HOOK('GY_M_M',1,ZHOOK_HANDLE)
END SUBROUTINE GY_M_M_PHY
+!
+! #######################################################
+ SUBROUTINE GX_M_U_PHY(D,OFLAT,PY,PDXX,PDZZ,PDZX,PGX_M_U)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ##################################################
+!
+!!**** *GX_M_U * - Compute the gradient along x for a variable localized at
+!! a mass point
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute a gradient along x
+! direction for a field PY localized at a mass point. The result PGX_M_U
+! is localized at a x-flux point (u point).
+!
+! ( ____________z )
+! ( ________x )
+! 1 ( dzm(PY) )
+! PGX_M_U = ---- (dxm(PY) - d*zx -------- )
+! d*xx ( d*zz )
+!
+!
+!
+!!** METHOD
+!! ------
+!! We employ the Shuman operators to compute the derivatives and the
+!! averages. The metric coefficients PDXX,PDZX,PDZZ are dummy arguments.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! FUNCTION DXM: compute a finite difference along the x direction for
+!! a variable at a mass localization
+!! FUNCTION DZM: compute a finite difference along the y direction for
+!! a variable at a mass localization
+!! FUNCTION MXM: compute an average in the x direction for a variable
+!! at a mass localization
+!! FUNCTION MZF: compute an average in the z direction for a variable
+!! at a flux side
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (function GX_M_U)
+!!
+!!
+!! AUTHOR
+!! ------
+!! P. Hereil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 05/07/94
+!! Modification 16/03/95 change the order of the arguments
+!! 19/07/00 add the LFLAT switch + inlining(J. Stein)
+!! 20/08/00 optimization (J. Escobar)
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT_TURB
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and result
+! ------------------------------------
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+LOGICAL, INTENT(IN) :: OFLAT
+REAL, DIMENSION(D%NIT*D%NJT*D%NKT), INTENT(IN) :: PY ! variable at the mass point
+REAL, DIMENSION(D%NIT*D%NJT*D%NKT), INTENT(IN) :: PDXX ! metric coefficient dyy
+REAL, DIMENSION(D%NIT*D%NJT*D%NKT), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(D%NIT*D%NJT*D%NKT), INTENT(IN) :: PDZX ! metric coefficient dzy
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGX_M_U ! result at flux
+ ! side
+REAL, DIMENSION(D%NIT*D%NJT*D%NKT) :: ZGX_M_U
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT):: ZY, ZDXX,ZDZZ,ZDZX
+INTEGER IIU,IKU,JI,JK
+!
+INTEGER :: JJK,IJU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+INTEGER :: JI_1JK, JIJK_1, JI_1JK_1, JIJKP1, JI_1JKP1
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE GRADIENT ALONG X
+! -----------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GX_M_U',0,ZHOOK_HANDLE)
+IIU=D%NIT
+IJU=D%NJT
+IKU=D%NKT
+IF (.NOT. OFLAT) THEN
+ JIJKOR = 1 + JPHEXT + IIU*IJU*(JPVEXT_TURB+1 - 1)
+ JIJKEND = IIU*IJU*(IKU-JPVEXT_TURB)
+!CDIR NODEP
+!OCL NOVREC
+ DO JIJK=JIJKOR , JIJKEND
+! indexation
+ JI_1JK = JIJK - 1
+ JIJK_1 = JIJK - IIU*IJU*D%NKL
+ JI_1JK_1 = JIJK - 1 - IIU*IJU*D%NKL
+ JIJKP1 = JIJK + IIU*IJU*D%NKL
+ JI_1JKP1 = JIJK - 1 + IIU*IJU*D%NKL
+!
+ ZGX_M_U(JIJK)= &
+ ( PY(JIJK)-PY(JI_1JK) &
+ -( (PY(JIJK)-PY(JIJK_1)) / PDZZ(JIJK) &
+ +(PY(JI_1JK)-PY(JI_1JK_1)) / PDZZ(JI_1JK) &
+ ) * PDZX(JIJK)* 0.25 &
+ -( (PY(JIJKP1)-PY(JIJK)) / PDZZ(JIJKP1) &
+ +(PY(JI_1JKP1)-PY(JI_1JK)) / PDZZ(JI_1JKP1) &
+ ) * PDZX(JIJKP1)* 0.25 &
+ ) / PDXX(JIJK)
+ END DO
+
+CALL D1D_TO_3D(D,ZGX_M_U,PGX_M_U)
+CALL D1D_TO_3D(D,PDXX,ZDXX)
+CALL D1D_TO_3D(D,PDZZ,ZDZZ)
+CALL D1D_TO_3D(D,PDZX,ZDZX)
+CALL D1D_TO_3D(D,PY,ZY)
+!
+ DO JI=1+JPHEXT,IIU
+ PGX_M_U(JI,:,D%NKU)= ( ZY(JI,:,D%NKU)-ZY(JI-1,:,D%NKU) ) / ZDXX(JI,:,D%NKU)
+ PGX_M_U(JI,:,D%NKA)= -999.
+ END DO
+!
+ PGX_M_U(1,:,:)=PGX_M_U(IIU-2*JPHEXT+1,:,:)
+ELSE
+! PGX_M_U = DXM(PY) / PDXX
+ PGX_M_U(1+JPHEXT:IIU,:,:) = ( ZY(1+JPHEXT:IIU,:,:)-ZY(JPHEXT:IIU-1,:,:) ) &
+ / ZDXX(1+JPHEXT:IIU,:,:)
+!
+ PGX_M_U(1,:,:)=PGX_M_U(IIU-2*JPHEXT+1,:,:)
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GX_M_U',1,ZHOOK_HANDLE)
+END SUBROUTINE GX_M_U_PHY
+!
+ SUBROUTINE GY_M_V_PHY(D,OFLAT,PY,PDYY,PDZZ,PDZY,PGY_M_V)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ##################################################
+!
+!!**** *GY_M_V * - Compute the gradient along y for a variable localized at
+!! a mass point
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute a gradient along y
+! direction for a field PY localized at a mass point. The result PGY_M_V
+! is localized at a y-flux point (v point).
+!
+! ( ____________z )
+! ( ________y )
+! 1 ( dzm(PY) )
+! PGY_M_V = ---- (dym(PY) - d*zy -------- )
+! d*yy ( d*zz )
+!
+!
+!
+!
+!!** METHOD
+!! ------
+!! We employ the Shuman operators to compute the derivatives and the
+!! averages. The metric coefficients PDYY,PDZY,PDZZ are dummy arguments.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! FUNCTION DYM: compute a finite difference along the y direction for
+!! a variable at a mass localization
+!! FUNCTION DZM: compute a finite difference along the y direction for
+!! a variable at a mass localization
+!! FUNCTION MYM: compute an average in the x direction for a variable
+!! at a mass localization
+!! FUNCTION MZF: compute an average in the z direction for a variable
+!! at a flux side
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (function GY_M_V)
+!!
+!!
+!! AUTHOR
+!! ------
+!! P. Hereil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 05/07/94
+!! Modification 16/03/95 change the order of the arguments
+!! 19/07/00 add the LFLAT switch + inlining(J. Stein)
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT_TURB
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+! -------------------------------------
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+LOGICAL, INTENT(IN) :: OFLAT
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDYY !d*yy
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZY !d*zy
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ !d*zz
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PY ! variable at mass
+ ! localization
+
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),INTENT(OUT) :: PGY_M_V ! result at flux
+ ! side
+!REAL, DIMENSION(D%NIT*D%NJT*D%NKT) :: ZGY_M_V
+!REAL, DIMENSION(D%NIT,D%NJT,D%NKT):: ZY, ZDYY,ZDZZ,ZDZY
+INTEGER IJU,IKU,JI,JJ,JK
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE GRADIENT ALONG Y
+! ----------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GY_M_V',0,ZHOOK_HANDLE)
+IJU=D%NJT
+IKU=D%NKT
+IF (.NOT. OFLAT) THEN
+! PGY_M_V = ( DYM(PY) - MZF ( MYM( DZM(PY) /PDZZ ) * PDZY ) )/PDYY
+ DO JK=1+JPVEXT_TURB,IKU-JPVEXT_TURB
+ DO JJ=1+JPHEXT,IJU
+ PGY_M_V(:,JJ,JK)= &
+ ( PY(:,JJ,JK)-PY(:,JJ-1,JK) &
+ -( (PY(:,JJ,JK)-PY(:,JJ,JK-D%NKL)) / PDZZ(:,JJ,JK) &
+ +(PY(:,JJ-1,JK)-PY(:,JJ-D%NKL,JK-D%NKL)) / PDZZ(:,JJ-1,JK) &
+ ) * PDZY(:,JJ,JK)* 0.25 &
+ -( (PY(:,JJ,JK+D%NKL)-PY(:,JJ,JK)) / PDZZ(:,JJ,JK+D%NKL) &
+ +(PY(:,JJ-1,JK+D%NKL)-PY(:,JJ-1,JK)) / PDZZ(:,JJ-1,JK+D%NKL) &
+ ) * PDZY(:,JJ,JK+D%NKL)* 0.25 &
+ ) / PDYY(:,JJ,JK)
+ END DO
+ END DO
+!
+ DO JJ=1+JPHEXT,IJU
+ PGY_M_V(:,JJ,D%NKU)= ( PY(:,JJ,D%NKU)-PY(:,JJ-1,D%NKU) ) / PDYY(:,JJ,D%NKU)
+ PGY_M_V(:,JJ,D%NKA)= -999.
+ END DO
+!
+ PGY_M_V(:,1,:)=PGY_M_V(:,IJU-2*JPHEXT+1,:)
+ELSE
+! PGY_M_V = DYM(PY)/PDYY
+ PGY_M_V(:,1+JPHEXT:IJU,:) = ( PY(:,1+JPHEXT:IJU,:)-PY(:,JPHEXT:IJU-1,:) ) &
+ / PDYY(:,1+JPHEXT:IJU,:)
+!
+ PGY_M_V(:,1,:)=PGY_M_V(:,IJU-2*JPHEXT+1,:)
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GY_M_V',1,ZHOOK_HANDLE)
+END SUBROUTINE GY_M_V_PHY
+!
+SUBROUTINE D1D_TO_3D (D,P1D,P3D)
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+IMPLICIT NONE
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: P1D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: P3D
+
+P3D = P1D
+END SUBROUTINE D1D_TO_3D
END MODULE MODE_GRADIENT_M_PHY
diff --git a/src/common/aux/gradient_u_phy.F90 b/src/common/aux/gradient_u_phy.F90
index 17284c32d..a5937536b 100644
--- a/src/common/aux/gradient_u_phy.F90
+++ b/src/common/aux/gradient_u_phy.F90
@@ -69,7 +69,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! metric coefficient
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
+INTEGER :: JI,JJ,JK, IIB, IIE, IJB, IJE
!
!* 0.2 declaration of local variables
@@ -81,15 +81,147 @@ INTEGER :: JI,JJ,JK
!* 1. DEFINITION of GZ_U_UW
! ---------------------
!
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!
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)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+PGZ_U_UW(IIB:IIE,IJB:IJE,1:D%NKT)= PA_WORK(IIB:IIE,IJB:IJE,1:D%NKT) &
+ / PDZZ_WORK(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!
!----------------------------------------------------------------------------
!
END SUBROUTINE GZ_U_UW_PHY
+!
+ SUBROUTINE GX_U_M_PHY(D,OFLAT,PA,PDXX,PDZZ,PDZX,PGX_U_M)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! #######################################################
+!
+!!**** *GX_U_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian X
+!! direction for a variable placed at the
+!! U point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the X cartesian direction for a field PA placed at the
+! U point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___________x ) )
+! 1 ( (d*zx dzm(PA) ) )
+! PGX_U_M = ---- (dxf(PA) - (------------)) )
+! ___x ( ( ) )
+! d*xx ( ( 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
+!! --------
+!! MXF,MZF : Shuman functions (mean operators)
+!! DXF,DZF : 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 19/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE SHUMAN_PHY, ONLY : DZM_PHY, DXF_PHY, MXF_PHY, MZF_PHY
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+!
+LOGICAL, INTENT(IN) :: OFLAT
+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) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: PGX_U_M ! result mass point
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4
+INTEGER :: IIB,IJB,IIE,IJE
+INTEGER :: JI,JJ,JK
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GX_U_M
+! --------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GX_U_M',0,ZHOOK_HANDLE)
+!
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!
+CALL DXF_PHY(D,PA,ZWORK1)
+CALL MXF_PHY(D,PDXX,ZWORK2)
+
+IF (.NOT. OFLAT) THEN
+ CALL DZM_PHY(D,PA,ZWORK3)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) * PDZX(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,ZWORK3,ZWORK4)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ CALL MZF_PHY(D,ZWORK4,ZWORK3)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ PGX_U_M(IIB:IIE,IJB:IJE,1:D%NKT) = ( ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) - ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)) &
+ / ZWORK2(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)
+ PGX_U_M(IIB:IIE,IJB:IJE,1:D%NKT)= ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GX_U_M',1,ZHOOK_HANDLE)
+END SUBROUTINE GX_U_M_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
index fc8a77c1b..392653606 100644
--- a/src/common/aux/gradient_v_phy.F90
+++ b/src/common/aux/gradient_v_phy.F90
@@ -71,10 +71,14 @@ 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
+INTEGER :: IIB,IJB,IIE,IJE
!
!* 0.2 declaration of local variables
!
-! NONE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
!
!----------------------------------------------------------------------------
!
@@ -84,11 +88,134 @@ INTEGER :: JI,JJ,JK
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)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+PGZ_V_VW(IIB:IIE,IJB:IJE,1:D%NKT)= PA_WORK(IIB:IIE,IJB:IJE,1:D%NKT) &
+ / PDZZ_WORK(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
!----------------------------------------------------------------------------
!
END SUBROUTINE GZ_V_VW_PHY
+ SUBROUTINE GY_V_M_PHY(D,OFLAT,PA,PDYY,PDZZ,PDZY,PGY_V_M)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! #######################################################
+!
+!!**** *GY_V_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Y
+!! direction for a variable placed at the
+!! V point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Y cartesian direction for a field PA placed at the
+! V point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___________y ) )
+! 1 ( (d*zy dzm(PA) ) )
+! PGY_V_M = ---- (dyf(PA) - (------------)) )
+! ___y ( ( ) )
+! d*yy ( ( 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
+!! --------
+!! MYF,MZF : Shuman functions (mean operators)
+!! DYF,DZF : 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 19/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE SHUMAN_PHY, ONLY : DZM_PHY, DYF_PHY, MYF_PHY, MZF_PHY
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+!
+LOGICAL, INTENT(IN) :: OFLAT
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PA ! variable at the V point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: PGY_V_M ! result mass point
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4
+INTEGER :: IIB,IJB,IIE,IJE
+INTEGER :: JI,JJ,JK
+!
+!* 0.2 declaration of local variables
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GY_V_M',0,ZHOOK_HANDLE)
+!
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GY_V_M
+! --------------------
+!
+CALL DYF_PHY(D,PA,ZWORK1)
+CALL MYF_PHY(D,PDYY,ZWORK2)
+!
+IF (.NOT. OFLAT) THEN
+ CALL DZM_PHY(D,PA,ZWORK3)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) * PDZY(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,ZWORK3,ZWORK4)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ CALL MZF_PHY(D,ZWORK4,ZWORK3)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ PGY_V_M(IIB:IIE,IJB:IJE,1:D%NKT) = ( ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) - ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)) &
+ / ZWORK2(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)
+ PGY_V_M(IIB:IIE,IJB:IJE,1:D%NKT)= ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GY_V_M',1,ZHOOK_HANDLE)
+END SUBROUTINE GY_V_M_PHY
+!
END MODULE MODE_GRADIENT_V_PHY
diff --git a/src/common/aux/gradient_w_phy.F90 b/src/common/aux/gradient_w_phy.F90
index f542d3654..e4e3dcf93 100644
--- a/src/common/aux/gradient_w_phy.F90
+++ b/src/common/aux/gradient_w_phy.F90
@@ -1,6 +1,235 @@
MODULE MODE_GRADIENT_W_PHY
IMPLICIT NONE
CONTAINS
+ SUBROUTINE GX_W_UW_PHY(D,OFLAT,PA,PDXX,PDZZ,PDZX,PGX_W_UW)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! #########################################################
+!
+!!**** *GX_W_UW* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian X
+!! direction for a variable placed at the
+!! V 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 X cartesian direction for a field PA placed at the
+! W point. The result is placed at the UW vorticity point.
+!
+!!** 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,MZM,MZF : Shuman functions (mean operators)
+!! DXM,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
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE SHUMAN_PHY, ONLY: MZF_PHY, DZF_PHY, MXM_PHY, DXM_PHY, MZM_PHY, DZM_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 W point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZX ! metric coefficient dzx
+LOGICAL, INTENT(IN) :: OFLAT
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),INTENT(OUT) :: PGX_W_UW ! result UW point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4, ZWORK5
+INTEGER :: IIB,IJB,IIE,IJE
+INTEGER :: JI,JJ,JK
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GX_W_UW
+! ---------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GX_W_UW',0,ZHOOK_HANDLE)
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+CALL MZM_PHY(D,PDXX,ZWORK1)
+CALL DXM_PHY(D,PA,ZWORK2)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK2(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)
+!
+IF (.NOT. OFLAT) THEN
+ CALL MZF_PHY(D,PA,ZWORK2)
+ CALL MXM_PHY(D,ZWORK2,ZWORK4)
+ CALL DZM_PHY(D,ZWORK4,ZWORK5)
+ !
+ CALL MXM_PHY(D,PDZZ,ZWORK2)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ PGX_W_UW(IIB:IIE,IJB:IJE,1:D%NKT)= ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) &
+ - ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)*PDZX(IIB:IIE,IJB:IJE,1:D%NKT) &
+ / (ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT))
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ELSE
+ PGX_W_UW = ZWORK3
+END IF
+!
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GX_W_UW',1,ZHOOK_HANDLE)
+END SUBROUTINE GX_W_UW_PHY
+!
+ SUBROUTINE GY_W_VW_PHY(D,OFLAT,PA,PDYY,PDZZ,PDZY,PGY_W_VW)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! #########################################################
+!
+!!**** *GY_W_VW* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Y
+!! direction for a variable placed at the
+!! W 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 Y cartesian direction for a field PA placed at the
+! W point. The result is placed at the VW vorticity point.
+!
+!!** 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,MZM,MZF : Shuman functions (mean operators)
+!! DYM,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
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE SHUMAN_PHY, ONLY: MZF_PHY, DZF_PHY, MYM_PHY, DYM_PHY, MZM_PHY, DZM_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 W point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDYY ! metric coefficient dxx
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZY ! metric coefficient dzx
+LOGICAL, INTENT(IN) :: OFLAT
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),INTENT(OUT) :: PGY_W_VW ! result UW point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4, ZWORK5
+INTEGER :: IIB,IJB,IIE,IJE
+INTEGER :: JI,JJ,JK
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GY_W_VW
+! ---------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GY_W_VW',0,ZHOOK_HANDLE)
+!IF (.NOT. LFLAT) THEN
+! PGY_W_VW(:,:,:)= DYM(PA(:,:,:))/(MZM(PDYY(:,:,:), KKA, KKU, KL)) &
+! -DZM(MYM(MZF(PA(:,:,:), KKA, KKU, KL)), KKA, KKU, KL)*PDZY(:,:,:) &
+! /( MZM(PDYY(:,:,:), KKA, KKU, KL)*MYM(PDZZ(:,:,:)) )
+!ELSE
+! PGY_W_VW(:,:,:)= DYM(PA(:,:,:))/(MZM(PDYY(:,:,:), KKA, KKU, KL))
+!END IF
+!
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+CALL MZM_PHY(D,PDYY,ZWORK1)
+CALL DYM_PHY(D,PA,ZWORK2)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK2(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)
+!
+IF (.NOT. OFLAT) THEN
+ CALL MZF_PHY(D,PA,ZWORK2)
+ CALL MYM_PHY(D,ZWORK2,ZWORK4)
+ CALL DZM_PHY(D,ZWORK4,ZWORK5)
+ !
+ CALL MYM_PHY(D,PDZZ,ZWORK2)
+ !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ PGY_W_VW(IIB:IIE,IJB:IJE,1:D%NKT)= ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) &
+ - ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)*PDZY(IIB:IIE,IJB:IJE,1:D%NKT) &
+ / (ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT))
+ !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ELSE
+ PGY_W_VW = ZWORK3
+END IF
+
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GY_W_VW',1,ZHOOK_HANDLE)
+END SUBROUTINE GY_W_VW_PHY
+!
SUBROUTINE GZ_W_M_PHY(D,PA,PDZZ,PGZ_W_M)
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
diff --git a/src/common/aux/shuman_phy.F90 b/src/common/aux/shuman_phy.F90
index dbe34cd2c..4f0c0367a 100644
--- a/src/common/aux/shuman_phy.F90
+++ b/src/common/aux/shuman_phy.F90
@@ -156,8 +156,8 @@ IMPLICIT NONE
! ------------------------------------
!
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) :: PMYM ! result at flux localization
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMYM ! result at flux localization
!
!* 0.2 Declarations of local variables
! -------------------------------
@@ -411,8 +411,8 @@ IMPLICIT NONE
! ------------------------------------
!
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) :: PMXM ! result at flux localization
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMXM ! result at flux localization
!
!* 0.2 Declarations of local variables
! -------------------------------
@@ -698,6 +698,190 @@ PDZF(:,:,D%NKU) = -999.
IF (LHOOK) CALL DR_HOOK('DZF',1,ZHOOK_HANDLE)
END SUBROUTINE DZF_PHY
!
+! ###############################
+ SUBROUTINE DYM_PHY(D,PA,PDYM)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *DYM* - Shuman operator : finite difference operator in y direction
+!! for a variable at a mass localization
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a finite difference
+! 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 PDYM(:,j,:) is defined by (PA(:,j,:)-PA(:,j-1,:))
+!! At j=1, PDYM(:,1,:) are replaced by the values of PDYM,
+!! 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 05/07/94
+!! Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS, ONLY: JPHEXT
+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) :: PDYM ! result at flux
+ ! side
+!
+!
+!* 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 DYM
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('DYM',0,ZHOOK_HANDLE)
+IJU=SIZE(PA,2)
+!
+DO JJ=2,IJU
+ PDYM(:,JJ,:) = PA(:,JJ,:) - PA(:,JJ-1,:)
+END DO
+!
+PDYM(:,1,:) = PDYM(:,IJU-2*JPHEXT+1,:)
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('DYM',1,ZHOOK_HANDLE)
+END SUBROUTINE DYM_PHY
+!
+! ###############################
+ SUBROUTINE DXM_PHY(D,PA,PDXM)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###############################
+!
+!!**** *DXM* - Shuman operator : finite difference operator in x direction
+!! for a variable at a mass localization
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a finite difference
+! 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 PDXM(i,:,:) is defined by (PA(i,:,:)-PA(i-1,:,:))
+!! At i=1, PDXM(1,:,:) are replaced by the values of PDXM,
+!! 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 05/07/94
+!! Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS, ONLY: JPHEXT
+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) :: PDXM ! result at flux
+ ! side
+!
+!* 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 DXM
+! ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('DXM',0,ZHOOK_HANDLE)
+IIU = SIZE(PA,1)
+!
+DO JI=2,IIU
+ PDXM(JI,:,:) = PA(JI,:,:) - PA(JI-1,:,:)
+END DO
+!
+PDXM(1,:,:) = PDXM(IIU-2*JPHEXT+1,:,:)
+!
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('DXM',1,ZHOOK_HANDLE)
+END SUBROUTINE DXM_PHY
! ###############################
SUBROUTINE DXF_PHY(D,PA,PDXF)
USE PARKIND1, ONLY : JPRB
diff --git a/src/mesonh/aux/shuman_phy.f90 b/src/mesonh/aux/shuman_phy.f90
index fa3f78966..e2a3a15f0 100644
--- a/src/mesonh/aux/shuman_phy.f90
+++ b/src/mesonh/aux/shuman_phy.f90
@@ -893,6 +893,212 @@ END DO
!
END SUBROUTINE DXF_PHY
!
+! ###############################
+ SUBROUTINE DXM_PHY(D,PA,PDXM)
+! ###############################
+!
+!!**** *DXM* - Shuman operator : finite difference operator in x direction
+!! for a variable at a mass localization
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a finite difference
+! 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 PDXM(i,:,:) is defined by (PA(i,:,:)-PA(i-1,:,:))
+!! At i=1, PDXM(1,:,:) are replaced by the values of PDXM,
+!! 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 05/07/94
+!! Modification to include the periodic case 13/10/94 J.Stein
+!! optimisation 20/08/00 J. Escobar
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS, ONLY: JPHEXT
+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) :: PDXM ! result at flux
+ ! side
+!
+!* 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 DXM
+! ------------------
+!
+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
+ PDXM(JIJK,1,1) = PA(JIJK,1,1) - PA(JIJK-1,1,1)
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JI=1,JPHEXT
+ DO JJK=1,IJU*IKU
+ PDXM(JI,JJK,1) = PDXM(IIU-2*JPHEXT+JI,JJK,1) ! for reprod JPHEXT <> 1
+ END DO
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE DXM_PHY
+!
+! ###############################
+ SUBROUTINE DYM_PHY(D,PA,PDYM)
+! ###############################
+!
+!!**** *DYM* - Shuman operator : finite difference operator in y direction
+!! for a variable at a mass localization
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a finite difference
+! 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 PDYM(:,j,:) is defined by (PA(:,j,:)-PA(:,j-1,:))
+!! At j=1, PDYM(:,1,:) are replaced by the values of PDYM,
+!! 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 05/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, ONLY: JPHEXT
+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) :: PDYM ! result at flux
+ ! side
+!
+!* 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 DYM
+! ------------------
+!
+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
+ PDYM(JIJK,1,1) = PA(JIJK,1,1) - PA(JIJK-IIU,1,1)
+END DO
+!
+DO JJ=1,JPHEXT
+ PDYM(:,JJ,:) = PDYM(:,IJU-2*JPHEXT+JJ,:) ! for reprod JPHEXT <> 1
+END DO
+!
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE DYM_PHY
! ###############################
SUBROUTINE DYF_PHY(D,PA,PDYF)
! ###############################
--
GitLab