diff --git a/MNH/contrav.f90 b/MNH/contrav.f90
new file mode 100644
index 0000000000000000000000000000000000000000..1a8dad3f09bae7332b8bf9b1dd796b9234f93917
--- /dev/null
+++ b/MNH/contrav.f90
@@ -0,0 +1,183 @@
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/contrav.f90,v $ $Revision: 1.1.10.1 $
+! MASDEV4_7 operators 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+! ####################
+ MODULE MODI_CONTRAV
+! ####################
+!
+INTERFACE
+!
+ SUBROUTINE CONTRAV(PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PRUCT,PRVCT,PRWCT )
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRUT ! Cartesian comp along x
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Cartesian comp along y
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRWT ! Cartesian comp along z
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRUCT ! Contrav comp along x-bar
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRVCT ! Contrav comp along y-bar
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRWCT ! Contrav comp along z-bar
+!
+END SUBROUTINE CONTRAV
+!
+END INTERFACE
+!
+END MODULE MODI_CONTRAV
+!
+!
+!
+! ##############################################################
+ SUBROUTINE CONTRAV(PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PRUCT,PRVCT,PRWCT )
+! ##############################################################
+!
+!!**** *CONTRAV * - computes the contravariant components from the
+!! cartesian components
+!!
+!! PURPOSE
+!! -------
+! This routine computes the contravariant components of vector
+! defined by its cartesian components (U,V,W) , using the following
+! formulae:
+! UC = U / DXX
+! VC = V / DYY
+! ( ----------x ----------y )
+! ( ---z ---z )
+! 1 ( U V )
+! WC = --- ( W - DZX * --- - DZY * --- )
+! DZZ ( DXX DYY )
+!
+!
+! In the no-topography case, WC = W / DZZ
+!
+!
+!!** METHOD
+!! ------
+!! We employ the Shuman operators to compute the averages. The metric
+!! coefficients PDXX, PDYY, PDZX, PDZY, PDZZ are dummy arguments
+!!
+!!
+!! EXTERNAL
+!! --------
+!! MXF, MYF, MZM : Shuman functions (mean operators)
+!!
+!! Module MODI_SHUMAN : Interface for Shuman functions
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! Module MODD_CONF : contains configuration variable
+!! LFLAT : Logical for topography
+!! = .TRUE. if Zs = 0 (Flat terrain)
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (subroutine CONTRAV)
+!!
+!!
+!! AUTHOR
+!! ------
+!! J.L. Redelsperger * CNRM *
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 27/07/94
+!! Corrections 3/08/94 (by J.P. Lafore)
+!! Corrections 17/10/94 (by J.P. Lafore) WC modified for w-advection
+!----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODD_CONF
+USE MODD_PARAMETERS
+!
+!
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRUT ! Cartesian comp along x
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Cartesian comp along y
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRWT ! Cartesian comp along z
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! Metric coefficients
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRUCT ! Contrav comp along x-bar
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRVCT ! Contrav comp along y-bar
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRWCT ! Contrav comp along z-bar
+!
+!
+!* 0.2 declarations of local variables
+!
+REAL, DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: Z1,Z2
+INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE
+!
+!-----------------------------------------------------------------------
+!
+!* 1. Compute the contravariant components
+! ------------------------------------
+!
+IIB=2
+IJB=2
+IIE=SIZE(PDXX,1)-1
+IJE=SIZE(PDXX,2)-1
+!
+IKB=1+JPVEXT
+IKE=SIZE(PDXX,3) - JPVEXT
+!
+!$acc data region copyin (pdxx,pdyy,pdzz,pdzx,pdzy,prut,prvt,prwt) local (z1,z2)
+!$acc region copyout (PRUCT,PRVCT,prwct)
+PRUCT(:,:,:) = PRUT(:,:,:) / PDXX(:,:,:)
+PRVCT(:,:,:) = PRVT(:,:,:) / PDYY(:,:,:)
+
+IF (LFLAT) THEN
+
+ PRWCT(:,:,:) = PRWT(:,:,:) / PDZZ(:,:,:)
+
+ELSE
+
+ Z1(IIB:IIE,:,IKB:IKE+1)= &
+ (PRUCT(IIB:IIE,:,IKB:IKE+1)+PRUCT(IIB:IIE,:,IKB-1:IKE) ) &
+ *PDZX(IIB:IIE,:,IKB:IKE+1) *0.25 &
+ +(PRUCT(IIB+1:IIE+1,:,IKB:IKE+1)+PRUCT(IIB+1:IIE+1,:,IKB-1:IKE) ) &
+ *PDZX(IIB+1:IIE+1,:,IKB:IKE+1) *0.25
+
+ Z2(:,IJB:IJE,IKB:IKE+1)= &
+ (PRVCT(:,IJB:IJE,IKB:IKE+1)+PRVCT(:,IJB:IJE,IKB-1:IKE) ) &
+ *PDZY(:,IJB:IJE,IKB:IKE+1) *0.25 &
+ +(PRVCT(:,IJB+1:IJE+1,IKB:IKE+1)+PRVCT(:,IJB+1:IJE+1,IKB-1:IKE) ) &
+ *PDZY(:,IJB+1:IJE+1,IKB:IKE+1) *0.25
+
+ PRWCT=0.
+ PRWCT(IIB:IIE,IJB:IJE,IKB:IKE+1) = &
+ ( PRWT(IIB:IIE,IJB:IJE,IKB:IKE+1) &
+ - Z1(IIB:IIE,IJB:IJE,IKB:IKE+1) &
+ - Z2(IIB:IIE,IJB:IJE,IKB:IKE+1) &
+ ) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE+1)
+
+END IF
+!
+
+PRWCT(:,:,1) = - PRWCT(:,:,3) ! Mirror hypothesis
+
+!$acc end region
+!$acc end data region
+!
+!-----------------------------------------------------------------------
+!
+END SUBROUTINE CONTRAV
diff --git a/MNH/get_halo.f90 b/MNH/get_halo.f90
index 2a2ccc43d7a0f7cfb2d1b530d2a7a57863f2c2db..31502d53a8b04e191ca0c79e5f6ddc5ade799b07 100644
--- a/MNH/get_halo.f90
+++ b/MNH/get_halo.f90
@@ -1,7 +1,7 @@
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
-! $Source$ $Revision$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/get_halo.f90,v $ $Revision: 1.1.2.1.2.2 $
! MASDEV4_7 newsrc 2007/03/01 13:18:33
!-----------------------------------------------------------------
! ####################
@@ -87,9 +87,20 @@ INTEGER :: IERROR ! error return code
!
NULLIFY( TZ_PSRC_ll)
!
+! acc update host (PSRC( IIB:IIE , IJB :IJB+JPHEXT , : ))
+! acc update host (PSRC( IIB:IIE , IJE-JPHEXT :IJE , : ))
+! acc update host (PSRC( IIB:IIB+JPHEXT , IJB+JPHEXT+1:IJE-JPHEXT-1 , : ))
+! acc update host (PSRC( IIE-JPHEXT:IIE , IJB+JPHEXT+1:IJE-JPHEXT-1 , : ))
+
CALL ADD3DFIELD_ll(TZ_PSRC_ll,PSRC)
CALL UPDATE_HALO_ll(TZ_PSRC_ll,IERROR, HDIR=HDIR )
CALL CLEANLIST_ll(TZ_PSRC_ll)
+
+! acc update device (PSRC( IIB:IIE , : IJB-1 , : ))
+! acc update device (PSRC( IIB:IIE , IJE+1: , : ))
+! acc update device (PSRC( :IIB-1 , IJB :IJE , : ))
+! acc update device (PSRC( IIE+1: , IJB :IJE , : ))
+
!
END SUBROUTINE GET_HALO
!-----------------------------------------------------------------------
diff --git a/MNH/gradient_m.f90 b/MNH/gradient_m.f90
new file mode 100644
index 0000000000000000000000000000000000000000..59c0b435b1f12b49783eecdda2e206a27266cec9
--- /dev/null
+++ b/MNH/gradient_m.f90
@@ -0,0 +1,743 @@
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/gradient_m.f90,v $ $Revision: 1.1.10.1 $
+! MASDEV4_7 operators 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+! ######################
+ MODULE MODI_GRADIENT_M
+! ######################
+!
+INTERFACE
+!
+!
+FUNCTION GX_M_M(PA,PDXX,PDZZ,PDZX) RESULT(PGX_M_M)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_M_M ! result mass point
+!
+END FUNCTION GX_M_M
+!
+!
+FUNCTION GY_M_M(PA,PDYY,PDZZ,PDZY) RESULT(PGY_M_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_M_M ! result mass point
+!
+END FUNCTION GY_M_M
+!
+!
+FUNCTION GZ_M_M(PA,PDZZ) RESULT(PGZ_M_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_M_M ! result mass point
+!
+END FUNCTION GZ_M_M
+!
+ FUNCTION GX_M_U(PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
+!
+IMPLICIT NONE
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGX_M_U ! result at flux
+ ! side
+END FUNCTION GX_M_U
+!
+!
+ FUNCTION GY_M_V(PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
+!
+IMPLICIT NONE
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY !d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY !d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGY_M_V ! result at flux
+ ! side
+END FUNCTION GY_M_V
+!
+ FUNCTION GZ_M_W(PY,PDZZ) RESULT(PGZ_M_W)
+!
+IMPLICIT NONE
+!
+ ! Metric coefficient:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGZ_M_W ! result at flux
+ ! side
+!
+END FUNCTION GZ_M_W
+!
+END INTERFACE
+!
+END MODULE MODI_GRADIENT_M
+!
+!
+!
+! #######################################################
+ FUNCTION GX_M_M(PA,PDXX,PDZZ,PDZX) RESULT(PGX_M_M)
+! #######################################################
+!
+!!**** *GX_M_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian X
+!! direction for a variable placed at the
+!! mass 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
+! mass point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___x ) )
+! 1 ( _x (d*zx dzm(PA) ) )
+! PGX_M_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
+!! --------
+!! MXM,MXF,MZF : Shuman functions (mean operators)
+!! DXF,DZF : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! 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 18/07/94
+!! 19/07/00 add the LFLAT switch (J. Stein)
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_M_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GX_M_M
+! --------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGX_M_M(:,:,:)= (DXF(MXM(PA(:,:,:))) - &
+ MZF(MXF(PDZX)*DZM(PA(:,:,:))/PDZZ(:,:,:)) ) &
+ /MXF(PDXX(:,:,:))
+ELSE
+ PGX_M_M(:,:,:)=DXF(MXM(PA(:,:,:))) / MXF(PDXX(:,:,:))
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GX_M_M
+!
+!
+! #######################################################
+ FUNCTION GY_M_M(PA,PDYY,PDZZ,PDZY) RESULT(PGY_M_M)
+! #######################################################
+!
+!!**** *GY_M_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Y
+!! direction for a variable placed at the
+!! mass 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
+! mass point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___y ) )
+! 1 ( _y (d*zy dzm(PA) ) )
+! PGY_M_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
+!! --------
+!! MYM,MYF,MZF : Shuman functions (mean operators)
+!! DYF,DZF : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! 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 18/07/94
+!! 19/07/00 add the LFLAT switch (J. Stein)
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODD_CONF
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_M_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GY_M_M
+! --------------------
+!
+!
+IF (.NOT. LFLAT) THEN
+ PGY_M_M(:,:,:)= (DYF(MYM(PA)) - MZF(MYF(PDZY)*DZM(PA)/PDZZ) ) &
+ /MYF(PDYY)
+ELSE
+ PGY_M_M(:,:,:)= DYF(MYM(PA))/MYF(PDYY)
+ENDIF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GY_M_M
+
+!
+!
+!
+! #######################################################
+ FUNCTION GZ_M_M(PA,PDZZ) RESULT(PGZ_M_M)
+! #######################################################
+!
+!!**** *GZ_M_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Z
+!! direction for a variable placed at the
+!! mass point and the result is placed at
+!! the mass 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
+! mass point. The result is placed at the mass point.
+!
+! _________z
+! (dzm(PA))
+! PGZ_M_M = (------ )
+! ( 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
+!! --------
+!! MZF : 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 18/07/94
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_M_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GZ_M_M
+! --------------------
+!
+PGZ_M_M(:,:,:)= MZF( DZM(PA(:,:,:))/PDZZ(:,:,:) )
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GZ_M_M
+!
+!
+! ##################################################
+ FUNCTION GX_M_U(PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
+! ##################################################
+!
+!!**** *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 MODI_SHUMAN
+USE MODD_CONF
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and result
+! ------------------------------------
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGX_M_U ! result at flux
+ ! side
+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
+! -----------------------------
+!
+IIU=SIZE(PY,1)
+IJU=SIZE(PY,2)
+IKU=SIZE(PY,3)
+
+IF (.NOT. LFLAT) THEN
+! PGX_M_U = ( DXM(PY) - MZF ( MXM( DZM(PY) /PDZZ ) * PDZX ) )/PDXX
+
+!$acc data region copyin (pdzz,pdzx,pdxx,py)
+!$acc region
+ DO JK=1+JPVEXT,IKU-JPVEXT
+ DO JI=1+JPHEXT,IIU
+ PGX_M_U(JI,:,JK)= &
+ ( PY(JI,:,JK)-PY(JI-1,:,JK) &
+ -( (PY(JI,:,JK)-PY(JI,:,JK-1)) / PDZZ(JI,:,JK) &
+ +(PY(JI-1,:,JK)-PY(JI-1,:,JK-1)) / PDZZ(JI-1,:,JK) &
+ ) * PDZX(JI,:,JK)* 0.25 &
+ -( (PY(JI,:,JK+1)-PY(JI,:,JK)) / PDZZ(JI,:,JK+1) &
+ +(PY(JI-1,:,JK+1)-PY(JI-1,:,JK)) / PDZZ(JI-1,:,JK+1) &
+ ) * PDZX(JI,:,JK+1)* 0.25 &
+ ) / PDXX(JI,:,JK)
+ END DO
+ END DO
+!$acc end region
+
+!!$ JIJKOR = 1 + JPHEXT + IIU*IJU*(JPVEXT+1 - 1)
+!!$ JIJKEND = IIU*IJU*(IKU-JPVEXT)
+!!$!CDIR NODEP
+!!$!OCL NOVREC
+!!$ DO JIJK=JIJKOR , JIJKEND
+!!$! indexation
+!!$ JI_1JK = JIJK - 1
+!!$ JIJK_1 = JIJK - IIU*IJU
+!!$ JI_1JK_1 = JIJK - 1 - IIU*IJU
+!!$ JIJKP1 = JIJK + IIU*IJU
+!!$ JI_1JKP1 = JIJK - 1 + IIU*IJU
+!!$!
+!!$ PGX_M_U(JIJK,1,1)= &
+!!$ ( PY(JIJK,1,1)-PY(JI_1JK,1,1) &
+!!$ -( (PY(JIJK,1,1)-PY(JIJK_1,1,1)) / PDZZ(JIJK,1,1) &
+!!$ +(PY(JI_1JK,1,1)-PY(JI_1JK_1,1,1)) / PDZZ(JI_1JK,1,1) &
+!!$ ) * PDZX(JIJK,1,1)* 0.25 &
+!!$ -( (PY(JIJKP1,1,1)-PY(JIJK,1,1)) / PDZZ(JIJKP1,1,1) &
+!!$ +(PY(JI_1JKP1,1,1)-PY(JI_1JK,1,1)) / PDZZ(JI_1JKP1,1,1) &
+!!$ ) * PDZX(JIJKP1,1,1)* 0.25 &
+!!$ ) / PDXX(JIJK,1,1)
+!!$ END DO
+
+!
+!$acc region
+ DO JI=1+JPHEXT,IIU
+ PGX_M_U(JI,:,IKU)= ( PY(JI,:,IKU)-PY(JI-1,:,IKU) ) / PDXX(JI,:,IKU)
+ PGX_M_U(JI,:,1)= -999.
+ END DO
+!
+ PGX_M_U(1,:,:)=PGX_M_U(IIU-2*JPHEXT+1,:,:)
+!$acc end region
+
+!$acc end data region
+ELSE
+! PGX_M_U = DXM(PY) / PDXX
+ PGX_M_U(1+JPHEXT:IIU,:,:) = ( PY(1+JPHEXT:IIU,:,:)-PY(JPHEXT:IIU-1,:,:) ) &
+ / PDXX(1+JPHEXT:IIU,:,:)
+!
+ PGX_M_U(1,:,:)=PGX_M_U(IIU-2*JPHEXT+1,:,:)
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION GX_M_U
+!
+!
+! ##################################################
+ FUNCTION GY_M_V(PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
+! ##################################################
+!
+!!**** *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 MODI_SHUMAN
+USE MODD_CONF
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+! -------------------------------------
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY !d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY !d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGY_M_V ! result at flux
+ ! side
+INTEGER IJU,IKU,JJ,JK
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE GRADIENT ALONG Y
+! ----------------------------
+!
+IJU=SIZE(PY,2)
+IKU=SIZE(PY,3)
+IF (.NOT. LFLAT) THEN
+! PGY_M_V = ( DYM(PY) - MZF ( MYM( DZM(PY) /PDZZ ) * PDZY ) )/PDYY
+!$acc region
+ DO JK=1+JPVEXT,IKU-JPVEXT
+ DO JJ=1+JPHEXT,IJU
+ PGY_M_V(:,JJ,JK)= &
+ ( PY(:,JJ,JK)-PY(:,JJ-1,JK) &
+ -( (PY(:,JJ,JK)-PY(:,JJ,JK-1)) / PDZZ(:,JJ,JK) &
+ +(PY(:,JJ-1,JK)-PY(:,JJ-1,JK-1)) / PDZZ(:,JJ-1,JK) &
+ ) * PDZY(:,JJ,JK)* 0.25 &
+ -( (PY(:,JJ,JK+1)-PY(:,JJ,JK)) / PDZZ(:,JJ,JK+1) &
+ +(PY(:,JJ-1,JK+1)-PY(:,JJ-1,JK)) / PDZZ(:,JJ-1,JK+1) &
+ ) * PDZY(:,JJ,JK+1)* 0.25 &
+ ) / PDYY(:,JJ,JK)
+ END DO
+ END DO
+!
+ DO JJ=1+JPHEXT,IJU
+ PGY_M_V(:,JJ,IKU)= ( PY(:,JJ,IKU)-PY(:,JJ-1,IKU) ) / PDYY(:,JJ,IKU)
+ PGY_M_V(:,JJ,1)= -999.
+ END DO
+!
+ PGY_M_V(:,1,:)=PGY_M_V(:,IJU-2*JPHEXT+1,:)
+!$acc end region
+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
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION GY_M_V
+!
+!
+! #########################################
+ FUNCTION GZ_M_W(PY,PDZZ) RESULT(PGZ_M_W)
+! #########################################
+!
+!!**** *GZ_M_W * - Compute the gradient along z direction for a
+!! variable localized at a mass point
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute a gradient along x,y,z
+! directions for a field PY localized at a mass point. The result PGZ_M_W
+! is localized at a z-flux point (w point)
+!
+!
+! dzm(PY)
+! PGZ_M_W = -------
+! d*zz
+!
+!!** METHOD
+!! ------
+!! We employ the Shuman operators to compute the derivatives and the
+!! averages. The metric coefficients PDZZ are dummy arguments.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! FUNCTION DZM : compute a finite difference along the z
+!! direction for a variable at a mass localization
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODI_SHUMAN : interface for the Shuman functions
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (function GZ_M_W)
+!!
+!!
+!! 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 inlining(J. Stein)
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_SHUMAN
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+! -------------------------------------
+!
+ ! Metric coefficient:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGZ_M_W ! result at flux
+ ! side
+!
+INTEGER IKU
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE GRADIENT ALONG Z
+! -----------------------------
+!
+IKU=SIZE(PY,3)
+PGZ_M_W(:,:,1+JPVEXT:IKU) = (PY(:,:,1+JPVEXT:IKU)-PY(:,:,JPVEXT:IKU-1)) &
+ / PDZZ(:,:,1+JPVEXT:IKU)
+PGZ_M_W(:,:,1)=-999.
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION GZ_M_W
+
diff --git a/MNH/metrics.f90 b/MNH/metrics.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b7685197bacefb8b73db3fb62be07126c9197b45
--- /dev/null
+++ b/MNH/metrics.f90
@@ -0,0 +1,174 @@
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/metrics.f90,v $ $Revision: 1.1.8.1.2.1 $ $Date: 2009/04/21 07:42:51 $
+!-----------------------------------------------------------------
+!-----------------------------------------------------------------
+!-----------------------------------------------------------------
+! ###################
+ MODULE MODI_METRICS
+! ###################
+INTERFACE
+!
+SUBROUTINE METRICS(PMAP,PDXHAT,PDYHAT,PZZ, &
+ PDXX,PDYY,PDZX,PDZY,PDZZ)
+REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! Map factor
+REAL, DIMENSION(:), INTENT(IN) :: PDXHAT ! Stretching in x direction
+REAL, DIMENSION(:), INTENT(IN) :: PDYHAT ! Stretching in y direction
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height in z direction
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDZX ! metric coefficient dzx
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDZY ! metric coefficient dzy
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDZZ ! metric coefficient dzz
+!
+END SUBROUTINE METRICS
+!
+END INTERFACE
+!
+END MODULE MODI_METRICS
+!
+!
+!
+! #################################################################
+ SUBROUTINE METRICS(PMAP,PDXHAT,PDYHAT,PZZ, &
+ PDXX,PDYY,PDZX,PDZY,PDZZ)
+! #################################################################
+!
+!!**** *METRICS* - routine to compute metric coefficients
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the metric coefficients
+! dxx,dyy,dzz,dzx,dzy
+!
+!!** METHOD
+!! ------
+!! The horizontal coefficients dxx and dyy (PDXX and PDYY arrays)
+!! are computed according to the thinshell or no thinshell approximation
+!! and to the cartesian or spherical geometry.
+!!
+!! EXTERNAL
+!! --------
+!! MXM,MYM,MZM : Shuman functions (mean operators)
+!! DXM,DYM,DZM : Shuman functions (finite differences operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : contains physical constants
+!!
+!! XRADIUS : earth radius
+!!
+!! Module MODD_CONF : contains configuration variables
+!!
+!! LTHINSHELL : Logical for thinshell approximation
+!! .TRUE. = Thinshell approximation done
+!! LCARTESIAN : Logical for cartesian geometry
+!! .TRUE. = Cartesian geometry used
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (routine METRICS)
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 12/07/94
+!! 14/02/01 (V. Masson and J. Stein) PDZZ initialized below the surface
+!! (influences the 3D turbulence of W) and PDXX,PDYY,PDZZ at the top
+!! 19/03/2008 (J.Escobar) remove spread !!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+USE MODD_CONF
+USE MODD_CST
+!
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! Map factor
+REAL, DIMENSION(:), INTENT(IN) :: PDXHAT ! Stretching in x direction
+REAL, DIMENSION(:), INTENT(IN) :: PDYHAT ! Stretching in y direction
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDZX ! metric coefficient dzx
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDZY ! metric coefficient dzy
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDZZ ! metric coefficient dzz
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: IIU ! Upper dimension in x direction
+INTEGER :: IJU ! Upper dimension in y direction
+INTEGER :: IKU ! Upper dimension in z direction
+REAL :: ZD1 ! DELTA1 (switch 0/1) for thinshell
+ ! approximation
+INTEGER :: JI,JJ,JK
+REAL, DIMENSION(SIZE(PDXHAT),SIZE(PDYHAT),SIZE(PZZ,3)) :: ZDZZ
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE DIMENSIONS OF ARRAYS :
+! ----------------------------
+IIU = SIZE(PDXHAT)
+IJU = SIZE(PDYHAT)
+IKU = SIZE(PZZ,3)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTE PDXX and PDYY :
+! --------------------
+!
+IF (LTHINSHELL) THEN
+ ZD1=0.
+ELSE
+ ZD1=1.
+END IF
+IF (.NOT.LCARTESIAN) THEN
+ ZDZZ(:,:,:) = MZF( 1.+ ZD1*PZZ(:,:,:)/XRADIUS)
+ DO JK=1,IKU ; DO JJ=1,IJU ; DO JI=1,IIU
+ PDXX(JI,JJ,JK) = ZDZZ(JI,JJ,JK) * PDXHAT(JI) /PMAP(JI,JJ)
+ PDYY(JI,JJ,JK) = ZDZZ(JI,JJ,JK) * PDYHAT(JJ) /PMAP(JI,JJ)
+ ENDDO ; ENDDO ; ENDDO
+ PDXX(:,:,:)=MXM(PDXX(:,:,:))
+ PDXX(:,:,IKU)=PDXX(:,:,IKU-1)
+ PDYY(:,:,:)=MYM(PDYY(:,:,:))
+ PDYY(:,:,IKU)=PDYY(:,:,IKU-1)
+ELSE
+ DO JK=1,IKU ; DO JJ=1,IJU ; DO JI=1,IIU
+ PDXX(JI,JJ,JK) = PDXHAT(JI)
+ PDYY(JI,JJ,JK) = PDYHAT(JJ)
+ ENDDO ; ENDDO ; ENDDO
+ PDXX(:,:,:)=MXM(PDXX(:,:,:))
+ PDYY(:,:,:)=MYM(PDYY(:,:,:))
+END IF
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. COMPUTE PDZX AND PDZY :
+! ----------------------
+!
+PDZX(:,:,:) = DXM(PZZ(:,:,:))
+!
+PDZY(:,:,:) = DYM(PZZ(:,:,:))
+!
+!-------------------------------------------------------------------------------
+!
+!* 4. COMPUTE PDZZ :
+! -------------
+!
+PDZZ(:,:,:) = DZM(MZF(PZZ(:,:,:)))
+PDZZ(:,:,IKU) = PZZ(:,:,IKU) - PZZ(:,:,IKU-1) ! same delta z in IKU and IKU -1
+PDZZ(:,:,1) = PDZZ(:,:,2) ! same delta z in 1 and 2
+!-----------------------------------------------------------------------------
+!
+END SUBROUTINE METRICS
diff --git a/MNH/ppm.f90 b/MNH/ppm.f90
index 7761f294fc2c64abba1187303fbb12d7e1ad3ae8..688423f7ffa5100f02bbfcff986ef0b0fd319ae9 100644
--- a/MNH/ppm.f90
+++ b/MNH/ppm.f90
@@ -1,7 +1,7 @@
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
-! $Source$ $Revision$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/ppm.f90,v $ $Revision: 1.1.2.3.2.2 $
! masdev4_7 BUG1 2007/06/15 17:47:18
!-----------------------------------------------------------------
! ###############
@@ -198,6 +198,7 @@ USE MODD_CONF
USE MODD_LUNIT
!END JUAN PPM_LL
USE MODE_MPPDB
+USE MODD_PARAMETERS, ONLY : JPHEXT
!
IMPLICIT NONE
!
@@ -236,6 +237,16 @@ REAL, DIMENSION(SIZE(ZCR,1),SIZE(ZCR,2),SIZE(ZCR,3)) :: ZFPOS, ZFNEG
INTEGER :: ILUOUT,IRESP ! for prints
INTEGER :: IJS,IJN
!END JUAN PPM_LL
+!JUAN ACC
+INTEGER :: I,J,K ,IIU,IJU,IKU
+LOGICAL :: GWEST , GEAST
+! inline shuman with macro
+#define dif2x(DQ,PQ) DQ(IIB:IIE,:,:)=0.5*(PQ(IIB+1:IIE+1,:,:)-PQ(IIB-1:IIE-1,:,:));\
+DQ(IIB-1,:,:)=0.5*(PQ(IIB,:,:)-PQ(IIE-1,:,:));\
+DQ(IIE+1,:,:)=0.5*(PQ(IIB+1,:,:)-PQ(IIE,:,:)) ! DIF2X(DQ,PQ)
+#define dxf(PDXF,PA) PDXF(1:IIU-1,:,:) = PA(2:IIU,:,:) - PA(1:IIU-1,:,:) ; PDXF(IIU,:,:) = PDXF(2*JPHEXT,:,:) ! DXF(PDXF,PA)
+#define mxm(PMXM,PA) PMXM(2:IIU,:,:) = 0.5*( PA(2:IIU,:,:)+PA(1:IIU-1,:,:) ) ; PMXM(1,:,:) = PMXM(IIU-2*JPHEXT+1,:,:) ! MXM(PMXM,PA)
+!JUAN ACC
!-------------------------------------------------------------------------------
!
!* 0.3. COMPUTES THE DOMAIN DIMENSIONS
@@ -244,6 +255,14 @@ INTEGER :: IJS,IJN
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IJS=IJB
IJN=IJE
+
+IIU=size(psrc,1)
+IJU=size(psrc,2)
+IKU=size(psrc,3)
+
+GWEST = LWEST_ll()
+GEAST = LEAST_ll()
+
!
!BEG JUAN PPM_LL
!
@@ -260,167 +279,166 @@ IF(NHALO /= 1) THEN
ENDIF
!
CALL GET_HALO(PSRC,HDIR="01_X")
-!
-PR=PSRC
-ZQL=PSRC
-ZQR=PSRC
-ZDQ=PSRC
-ZQ6=PSRC
-ZDMQ=PSRC
-ZQL0=PSRC
-ZQR0=PSRC
-ZQ60=PSRC
-ZFPOS=PSRC
-ZFNEG=PSRC
-!
-!-------------------------------------------------------------------------------
-!
-SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
-!
-! 1.1 CYCLIC BOUNDARY CONDITIONS IN X DIRECTION
-! -----------------------------------------
-!
-CASE ('CYCL','WALL') ! In that case one must have HLBCX(1) == HLBCX(2)
-!
-! calculate dmq ZDMQ(i) = Fct[ PSRC(i-1),PSRC(i+1) ]
- ZDMQ = DIF2X(PSRC)
-!
-! monotonize the difference followinq eq. 5 in Lin94
-!
-!BEG JUAN PPM_LL01
-!
-! ZDMQ(i) = Fct[ ZDMQ(i),PSRC(i),PSRC(i-1),PSRC(i+1) ]
-!
- ZDMQ(IIB:IIE,IJS:IJN,:) = &
- SIGN( (MIN( ABS(ZDMQ(IIB:IIE,IJS:IJN,:)),2.0*(PSRC(IIB:IIE,IJS:IJN,:) - &
- MIN(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:))), &
- 2.0*(MAX(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:)) - &
- PSRC(IIB:IIE,IJS:IJN,:)) )), ZDMQ(IIB:IIE,IJS:IJN,:) )
-!
-! WEST BOUND
-!
-!!$ ZDMQ(IIB-1,:,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(IIB-1,:,:)), 2.0*(PSRC(IIB-1,:,:) - &
-!!$ MIN(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:))), &
-!!$ 2.0*(MAX(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:)) - &
-!!$ PSRC(IIB-1,:,:)) )), ZDMQ(IIB-1,:,:) )
-!
-! EAST BOUND
-!
-!!$ ZDMQ(IIE+1,:,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(IIE+1,:,:)), 2.0*(PSRC(IIE+1,:,:) - &
-!!$ MIN(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:))), &
-!!$ 2.0*(MAX(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:)) - &
-!!$ PSRC(IIE+1,:,:)) )), ZDMQ(IIE+1,:,:) )
-!
-! update ZDMQ HALO before next/further utilisation
-!
- CALL GET_HALO(ZDMQ,HDIR="01_X") ! PB AVEC HDIR="X1"
- CALL MPPDB_CHECK3DM("PPM::PPM_01_X CYCL ::ZDMQ",PRECISION,ZDMQ)
-!
-! calculate qL and qR with the modified dmq
-!
-! ZQL0(i) = Fct[ PSRC(i),PSRC(i-1),ZDMQ(i),ZDMQ(i-1) ]
-!
- ZQL0(IIB:IIE+1,IJS:IJN,:) = 0.5*(PSRC(IIB:IIE+1,IJS:IJN,:) + PSRC(IIB-1:IIE,IJS:IJN,:)) - &
- (ZDMQ(IIB:IIE+1,IJS:IJN,:) - ZDMQ(IIB-1:IIE,IJS:IJN,:))/6.0
-!
- CALL GET_HALO(ZQL0,HDIR="01_X")
-!
-! WEST BOUND
-!
-!!$ ZQL0(IIB-1,:,:) = ZQL0(IIE,:,:) JUAN PPMLL01
-!
- ZQR0(IIB-1:IIE,IJS:IJN,:) = ZQL0(IIB:IIE+1,IJS:IJN,:)
-!
- CALL GET_HALO(ZQR0,HDIR="01_X")
-!
-! EAST BOUND
-!
-!!$ ZQR0(IIE+1,:,:) = ZQR0(IIB,:,:) JUAN PPMLL01
-!
-! determine initial coefficients of the parabolae
-!
- ZDQ = ZQR0 - ZQL0
- ZQ60 = 6.0*(PSRC - 0.5*(ZQL0 + ZQR0))
-!
-! initialize final parabolae parameters
-!
- ZQL = ZQL0
- ZQR = ZQR0
- ZQ6 = ZQ60
-!
-! eliminate over and undershoots and create qL and qR as in Lin96
-!
- WHERE ( ZDMQ == 0.0 )
- ZQL = PSRC
- ZQR = PSRC
- ZQ6 = 0.0
- ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 )
- ZQ6 = 3.0*(ZQL0 - PSRC)
- ZQR = ZQL0 - ZQ6
- ZQL = ZQL0
- ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
- ZQ6 = 3.0*(ZQR0 - PSRC)
- ZQL = ZQR0 - ZQ6
- ZQR = ZQR0
- END WHERE
-!
-! recalculate coefficients of the parabolae
-!
- ZDQ = ZQR - ZQL
-!
-! and finally calculate fluxes for the advection
-!
-! ZFPOS(i) = Fct[ ZQR(i-1),ZCR(i),ZDQ(i-1),ZQ6(i-1) ]
-!
- ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZQR(IIB-1:IIE,IJS:IJN,:) - 0.5*ZCR(IIB:IIE+1,IJS:IJN,:) * &
- (ZDQ(IIB-1:IIE,IJS:IJN,:) - (1.0 - 2.0*ZCR(IIB:IIE+1,IJS:IJN,:)/3.0) &
- * ZQ6(IIB-1:IIE,IJS:IJN,:))
-!
- CALL GET_HALO(ZFPOS,HDIR="Z1_X")
-!
-! WEST BOUND
-!
-! PPOSX(IIB-1,:,:) is not important for the calc of advection so
-! we set it to 0
-!!$ ZFPOS(IIB-1,:,:) = 0.0 JUANPPMLL01
-!
- ZFNEG(:,IJS:IJN,:) = ZQL(:,IJS:IJN,:) - 0.5*ZCR(:,IJS:IJN,:) * &
- ( ZDQ(:,IJS:IJN,:) + (1.0 + 2.0*ZCR(:,IJS:IJN,:)/3.0) * ZQ6(:,IJS:IJN,:) )
-!
- CALL GET_HALO(ZFNEG,HDIR="Z1_X")
-!
-! advect the actual field in X direction by U*dt
-!
- PR = DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
- CALL GET_HALO(PR,HDIR="01_X")
- CALL MPPDB_CHECK3DM("PPM::PPM_01_X CYCL ::PR",PRECISION,PR)
-!
-!
-!* 1.2 NON-CYCLIC BOUNDARY CONDITIONS IN THE X DIRECTION
-! -------------------------------------------------
-!
-CASE('OPEN')
+
+!$acc data region local (ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG) copyin (psrc,zcr,prho) copyout(pr) !
+!$acc region
+ DO K=1,IKU ; DO J = 1,IJU ; DO I=1,IIU
+PR(I,J,K)=PSRC(I,J,K)
+ZQL(I,J,K)=PSRC(I,J,K)
+ZQR(I,J,K)=PSRC(I,J,K)
+ZDQ(I,J,K)=PSRC(I,J,K)
+ZQ6(I,J,K)=PSRC(I,J,K)
+ZDMQ(I,J,K)=PSRC(I,J,K)
+ZQL0(I,J,K)=PSRC(I,J,K)
+ZQR0(I,J,K)=PSRC(I,J,K)
+ZQ60(I,J,K)=PSRC(I,J,K)
+ZFPOS(I,J,K)=PSRC(I,J,K)
+ZFNEG(I,J,K)=PSRC(I,J,K)
+ENDDO ; ENDDO ; ENDDO
+! acc end region
+!
+!-------------------------------------------------------------------------------
+!
+!!$SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
+!!$!
+!!$! 1.1 CYCLIC BOUNDARY CONDITIONS IN X DIRECTION
+!!$! -----------------------------------------
+!!$!
+!!$CASE ('CYCL','WALL') ! In that case one must have HLBCX(1) == HLBCX(2)
+!!$!
+!!$! calculate dmq ZDMQ(i) = Fct[ PSRC(i-1),PSRC(i+1) ]
+!!$ ZDMQ = DIF2X(PSRC)
+!!$!
+!!$! monotonize the difference followinq eq. 5 in Lin94
+!!$!
+!!$!BEG JUAN PPM_LL01
+!!$!
+!!$! ZDMQ(i) = Fct[ ZDMQ(i),PSRC(i),PSRC(i-1),PSRC(i+1) ]
+!!$!
+!!$ ZDMQ(IIB:IIE,IJS:IJN,:) = &
+!!$ SIGN( (MIN( ABS(ZDMQ(IIB:IIE,IJS:IJN,:)),2.0*(PSRC(IIB:IIE,IJS:IJN,:) - &
+!!$ MIN(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:))), &
+!!$ 2.0*(MAX(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:)) - &
+!!$ PSRC(IIB:IIE,IJS:IJN,:)) )), ZDMQ(IIB:IIE,IJS:IJN,:) )
+!!$!
+!!$! WEST BOUND
+!!$!
+!!$!
+!!$! EAST BOUND
+!!$!
+!!$!
+!!$! update ZDMQ HALO before next/further utilisation
+!!$!
+!!$ CALL GET_HALO(ZDMQ,HDIR="01_X") ! PB AVEC HDIR="X1"
+!!$ CALL MPPDB_CHECK3DM("PPM::PPM_01_X CYCL ::ZDMQ",PRECISION,ZDMQ)
+!!$!
+!!$! calculate qL and qR with the modified dmq
+!!$!
+!!$! ZQL0(i) = Fct[ PSRC(i),PSRC(i-1),ZDMQ(i),ZDMQ(i-1) ]
+!!$!
+!!$ ZQL0(IIB:IIE+1,IJS:IJN,:) = 0.5*(PSRC(IIB:IIE+1,IJS:IJN,:) + PSRC(IIB-1:IIE,IJS:IJN,:)) - &
+!!$ (ZDMQ(IIB:IIE+1,IJS:IJN,:) - ZDMQ(IIB-1:IIE,IJS:IJN,:))/6.0
+!!$!
+!!$ CALL GET_HALO(ZQL0,HDIR="01_X")
+!!$!
+!!$! WEST BOUND
+!!$!
+!!$!
+!!$ ZQR0(IIB-1:IIE,IJS:IJN,:) = ZQL0(IIB:IIE+1,IJS:IJN,:)
+!!$!
+!!$ CALL GET_HALO(ZQR0,HDIR="01_X")
+!!$!
+!!$! EAST BOUND
+!!$!
+!!$!
+!!$! determine initial coefficients of the parabolae
+!!$!
+!!$ ZDQ = ZQR0 - ZQL0
+!!$ ZQ60 = 6.0*(PSRC - 0.5*(ZQL0 + ZQR0))
+!!$!
+!!$! initialize final parabolae parameters
+!!$!
+!!$ ZQL = ZQL0
+!!$ ZQR = ZQR0
+!!$ ZQ6 = ZQ60
+!!$!
+!!$! eliminate over and undershoots and create qL and qR as in Lin96
+!!$!
+!!$ WHERE ( ZDMQ == 0.0 )
+!!$ ZQL = PSRC
+!!$ ZQR = PSRC
+!!$ ZQ6 = 0.0
+!!$ ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 )
+!!$ ZQ6 = 3.0*(ZQL0 - PSRC)
+!!$ ZQR = ZQL0 - ZQ6
+!!$ ZQL = ZQL0
+!!$ ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
+!!$ ZQ6 = 3.0*(ZQR0 - PSRC)
+!!$ ZQL = ZQR0 - ZQ6
+!!$ ZQR = ZQR0
+!!$ END WHERE
+!!$!
+!!$! recalculate coefficients of the parabolae
+!!$!
+!!$ ZDQ = ZQR - ZQL
+!!$!
+!!$! and finally calculate fluxes for the advection
+!!$!
+!!$! ZFPOS(i) = Fct[ ZQR(i-1),ZCR(i),ZDQ(i-1),ZQ6(i-1) ]
+!!$!
+!!$ ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZQR(IIB-1:IIE,IJS:IJN,:) - 0.5*ZCR(IIB:IIE+1,IJS:IJN,:) * &
+!!$ (ZDQ(IIB-1:IIE,IJS:IJN,:) - (1.0 - 2.0*ZCR(IIB:IIE+1,IJS:IJN,:)/3.0) &
+!!$ * ZQ6(IIB-1:IIE,IJS:IJN,:))
+!!$!
+!!$ CALL GET_HALO(ZFPOS,HDIR="Z1_X")
+!!$!
+!!$! WEST BOUND
+!!$!
+!!$! PPOSX(IIB-1,:,:) is not important for the calc of advection so
+!!$! we set it to 0
+!!$!
+!!$ ZFNEG(:,IJS:IJN,:) = ZQL(:,IJS:IJN,:) - 0.5*ZCR(:,IJS:IJN,:) * &
+!!$ ( ZDQ(:,IJS:IJN,:) + (1.0 + 2.0*ZCR(:,IJS:IJN,:)/3.0) * ZQ6(:,IJS:IJN,:) )
+!!$!
+!!$ CALL GET_HALO(ZFNEG,HDIR="Z1_X")
+!!$!
+!!$! advect the actual field in X direction by U*dt
+!!$!
+!!$ PR = DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+!!$ CALL GET_HALO(PR,HDIR="01_X")
+!!$ CALL MPPDB_CHECK3DM("PPM::PPM_01_X CYCL ::PR",PRECISION,PR)
+!!$!
+!!$!
+!!$!* 1.2 NON-CYCLIC BOUNDARY CONDITIONS IN THE X DIRECTION
+!!$! -------------------------------------------------
+!!$!
+!!$CASE('OPEN')
!
! calculate dmq
!
- ZDMQ = DIF2X(PSRC)
+!! ZDMQ = DIF2X(PSRC)
+ ! acc region
+ dif2x(ZDMQ,PSRC)
+ ! acc end region
!
! overwrite the values on the boundary to get second order difference
! for qL and qR at the boundary
!
! WEST BOUND
!
- IF (LWEST_ll()) THEN
+ IF (GWEST) THEN
+ ! acc region
ZDMQ(IIB-1,IJS:IJN,:) = -ZDMQ(IIB,IJS:IJN,:)
+ ! acc end region
ENDIF
!
! EAST BOUND
!
- IF (LEAST_ll()) THEN
+ IF (GEAST) THEN
+ ! acc region
ZDMQ(IIE+1,IJS:IJN,:) = -ZDMQ(IIE,IJS:IJN,:)
+ ! acc end region
ENDIF
!
! monotonize the difference followinq eq. 5 in Lin94
@@ -435,25 +453,15 @@ CASE('OPEN')
!
! WEST BOUND
!
-!!$ ZDMQ(IIB-1,:,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(IIB-1,:,:)), 2.0*(PSRC(IIB-1,:,:) - &
-!!$ MIN(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:))), &
-!!$ 2.0*(MAX(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:)) - &
-!!$ PSRC(IIB-1,:,:)) )), ZDMQ(IIB-1,:,:) )
!
! EAST BOUND
!
-!!$ ZDMQ(IIE+1,:,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(IIE+1,:,:)), 2.0*(PSRC(IIE+1,:,:) - &
-!!$ MIN(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:))), &
-!!$ 2.0*(MAX(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:)) - &
-!!$ PSRC(IIE+1,:,:)) )), ZDMQ(IIE+1,:,:) )
!
!
! update ZDMQ HALO before next/further utilisation
!
- CALL GET_HALO(ZDMQ,HDIR="01_X")
- CALL MPPDB_CHECK3DM("PPM::PPM_01_X OPEN ::ZDMQ",PRECISION,ZDMQ)
+!TEMPO_JUAN CALL GET_HALO(ZDMQ,HDIR="01_X")
+!TEMPO_JUAN CALL MPPDB_CHECK3DM("PPM::PPM_01_X OPEN ::ZDMQ",PRECISION,ZDMQ)
!
! calculate qL and qR
!
@@ -462,162 +470,177 @@ CASE('OPEN')
ZQL0(IIB:IIE+1,IJS:IJN,:) = 0.5*(PSRC(IIB:IIE+1,IJS:IJN,:) + PSRC(IIB-1:IIE,IJS:IJN,:)) - &
(ZDMQ(IIB:IIE+1,IJS:IJN,:) - ZDMQ(IIB-1:IIE,IJS:IJN,:))/6.0
!
- CALL GET_HALO(ZQL0,HDIR="01_X")
+!TEMPO_JUAN CALL GET_HALO(ZQL0,HDIR="01_X")
!
! WEST BOUND
!
- IF (LWEST_ll()) THEN
+ ! acc region
+ IF (GWEST) THEN
+ ! acc region
ZQL0(IIB-1,IJS:IJN,:) = ZQL0(IIB,IJS:IJN,:)
+ ! acc end region
ENDIF
!
ZQR0(IIB-1:IIE,IJS:IJN,:) = ZQL0(IIB:IIE+1,IJS:IJN,:)
!
- CALL GET_HALO(ZQR0,HDIR="01_X")
+!TEMPO_JUAN CALL GET_HALO(ZQR0,HDIR="01_X")
!
! EAST BOUND
!
- IF (LEAST_ll()) THEN
+ IF (GEAST) THEN
+ ! acc region
ZQR0(IIE+1,IJS:IJN,:) = ZQR0(IIE,IJS:IJN,:)
+ ! acc end region
ENDIF
+
+DO K=1,IKU ; DO J = 1,IJU ; DO I=1,IIU
!
! determine initial coefficients of the parabolae
!
- ZDQ = ZQR0 - ZQL0
- ZQ60 = 6.0*(PSRC - 0.5*(ZQL0 + ZQR0))
+ ZDQ (I,J,K)= ZQR0(I,J,K) - ZQL0(I,J,K)
+ ZQ60(I,J,K) = 6.0*(PSRC(I,J,K) - 0.5*(ZQL0(I,J,K) + ZQR0(I,J,K)))
!
! initialize final parabolae parameters
!
- ZQL = ZQL0
- ZQR = ZQR0
- ZQ6 = ZQ60
+ ZQL(I,J,K) = ZQL0(I,J,K)
+ ZQR(I,J,K) = ZQR0(I,J,K)
+ ZQ6(I,J,K) = ZQ60(I,J,K)
!
! eliminate over and undershoots and create qL and qR as in Lin96
!
- WHERE ( ZDMQ == 0.0 )
- ZQL = PSRC
- ZQR = PSRC
- ZQ6 = 0.0
- ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 )
- ZQ6 = 3.0*(ZQL0 - PSRC)
- ZQR = ZQL0 - ZQ6
- ZQL = ZQL0
- ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
- ZQ6 = 3.0*(ZQR0 - PSRC)
- ZQL = ZQR0 - ZQ6
- ZQR = ZQR0
- END WHERE
+ IF ( ZDMQ(I,J,K) == 0.0 ) THEN
+ ZQL(I,J,K) = PSRC(I,J,K)
+ ZQR(I,J,K) = PSRC(I,J,K)
+ ZQ6(I,J,K) = 0.0
+ ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -(ZDQ(I,J,K))**2 ) THEN
+ ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K))
+ ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K)
+ ZQL(I,J,K) = ZQL0(I,J,K)
+ ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) > (ZDQ(I,J,K))**2 ) THEN
+ ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K))
+ ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K)
+ ZQR(I,J,K) = ZQR0(I,J,K)
+ ENDIF
!
! recalculate coefficients of the parabolae
!
- ZDQ = ZQR - ZQL
+ ZDQ(I,J,K) = ZQR(I,J,K) - ZQL(I,J,K)
+ENDDO ; ENDDO ; ENDDO
!
! and finally calculate fluxes for the advection
!
!
! ZFPOS(i) = Fct[ ZQR(i-1),ZCR(i),ZDQ(i-1),ZQ6(i-1) ]
!
-!!$ ZFPOS(IIB+1:IIE+1,:,:) = ZQR(IIB:IIE,:,:) - 0.5*ZCR(IIB+1:IIE+1,:,:) * &
-!!$ (ZDQ(IIB:IIE,:,:) - (1.0 - 2.0*ZCR(IIB+1:IIE+1,:,:)/3.0) &
-!!$ * ZQ6(IIB:IIE,:,:))
ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZQR(IIB-1:IIE,IJS:IJN,:) - 0.5*ZCR(IIB:IIE+1,IJS:IJN,:) * &
(ZDQ(IIB-1:IIE,IJS:IJN,:) - (1.0 - 2.0*ZCR(IIB:IIE+1,IJS:IJN,:)/3.0) &
* ZQ6(IIB-1:IIE,IJS:IJN,:))
!
- CALL GET_HALO(ZFPOS,HDIR="01_X")
+!TEMPO_JUAN CALL GET_HALO(ZFPOSS,HDIR="01_X")
!
!
! WEST BOUND
!
! advection flux at open boundary when u(IIB) > 0
!
- IF (LWEST_ll()) THEN
+ IF (GWEST) THEN
+ ! acc region
ZFPOS(IIB,IJS:IJN,:) = (PSRC(IIB-1,IJS:IJN,:) - ZQR(IIB-1,IJS:IJN,:))*ZCR(IIB,IJS:IJN,:) + &
ZQR(IIB-1,IJS:IJN,:)
+ ! acc end region
! PPOSX(IIB-1,:,:) is not important for the calc of advection so
! we set it to 0
-!!$ ZFPOS(IIB-1,:,:) = 0.0
ENDIF
!
-!!$ ZFNEG(IIB-1:IIE,:,:) = ZQL(IIB-1:IIE,:,:) - 0.5*ZCR(IIB-1:IIE,:,:) * &
-!!$ (ZDQ(IIB-1:IIE,:,:) + (1.0 + 2.0*ZCR(IIB-1:IIE,:,:)/3.0) &
-!!$ * ZQ6(IIB-1:IIE,:,:))
ZFNEG(:,IJS:IJN,:) = ZQL(:,IJS:IJN,:) - 0.5*ZCR(:,IJS:IJN,:) * &
( ZDQ(:,IJS:IJN,:) + (1.0 + 2.0*ZCR(:,IJS:IJN,:)/3.0) * ZQ6(:,IJS:IJN,:) )
!
- CALL GET_HALO(ZFNEG,HDIR="01_X")
+!TEMPO_JUAN CALL GET_HALO(ZFNEG,HDIR="01_X")
!
! EAST BOUND
!
! advection flux at open boundary when u(IIE+1) < 0
- IF (LEAST_ll()) THEN
+ IF (GEAST) THEN
+ ! acc region
ZFNEG(IIE+1,IJS:IJN,:) = (ZQR(IIE,IJS:IJN,:)-PSRC(IIE+1,IJS:IJN,:))*ZCR(IIE+1,IJS:IJN,:) + &
ZQR(IIE,IJS:IJN,:)
+ ! acc end region
ENDIF
!
! advect the actual field in X direction by U*dt
!
- PR = DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+ ! acc region
+ mxm(ZQL,PRHO)
+ ZQR = ZCR* ZQL*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + ZFNEG*(0.5-SIGN(0.5,ZCR)) )
+ dxf(PR,ZQR)
+
+!!$ PR = DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+
+ !$acc end region
+! acc update host (pr)
+!$acc end data region
+
CALL GET_HALO(PR,HDIR="01_X")
!
+!!$END SELECT
!
-END SELECT
-!
-CONTAINS
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-! ########################################################################
- FUNCTION DIF2X(PQ) RESULT(DQ)
-! ########################################################################
-!!
-!!**** DIF2X - leap-frog difference operator in X direction
-!!
-!! Calculates the difference assuming periodic BC (CYCL).
-!!
-!! DQ(I) = 0.5 * (PQ(I+1) - PQ(I-1))
-!!
-!! MODIFICATIONS
-!! -------------
-!!
-!! 18.3.2006. T. Maric - original version
-!! 07/2010 J.Escobar : Correction for reproducility
-!-------------------------------------------------------------------------------
-!
-!
-USE MODE_ll
-!
-IMPLICIT NONE
-!
-!* 0.1 Declarations of dummy arguments :
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
-REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
-!
-!* 0.2 Declarations of local variables :
-!
-INTEGER :: IIB,IJB ! Begining useful area in x,y directions
-INTEGER :: IIE,IJE ! End useful area in x,y directions
-!
-!-------------------------------------------------------------------------------
-!
-!* 1.0. COMPUTE THE DOMAIN DIMENSIONS
-! -----------------------------
-!
-CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
-!
-!-------------------------------------------------------------------------------
-!
-!* 2.0. COMPUTE THE DIFFERENCE
-! ----------------------
-!
-DQ(IIB:IIE,:,:) = PQ(IIB+1:IIE+1,:,:) - PQ(IIB-1:IIE-1,:,:)
-DQ(IIB-1,:,:) = PQ(IIB,:,:) - PQ(IIE-1,:,:)
-DQ(IIE+1,:,:) = PQ(IIB+1,:,:) - PQ(IIE,:,:)
-DQ = 0.5*DQ
-!
-END FUNCTION DIF2X
+!!$CONTAINS
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$! ########################################################################
+!!$ FUNCTION DIF2X(PQ) RESULT(DQ)
+!!$! ########################################################################
+!!$!!
+!!$!!**** DIF2X - leap-frog difference operator in X direction
+!!$!!
+!!$!! Calculates the difference assuming periodic BC (CYCL).
+!!$!!
+!!$!! DQ(I) = 0.5 * (PQ(I+1) - PQ(I-1))
+!!$!!
+!!$!! MODIFICATIONS
+!!$!! -------------
+!!$!!
+!!$!! 18.3.2006. T. Maric - original version
+!!$!! 07/2010 J.Escobar : Correction for reproducility
+!!$!!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!
+!!$USE MODE_ll
+!!$!
+!!$IMPLICIT NONE
+!!$!
+!!$!* 0.1 Declarations of dummy arguments :
+!!$!
+!!$REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
+!!$REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
+!!$!
+!!$!* 0.2 Declarations of local variables :
+!!$!
+!!$INTEGER :: IIB,IJB ! Begining useful area in x,y directions
+!!$INTEGER :: IIE,IJE ! End useful area in x,y directions
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!* 1.0. COMPUTE THE DOMAIN DIMENSIONS
+!!$! -----------------------------
+!!$!
+!!$CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!* 2.0. COMPUTE THE DIFFERENCE
+!!$! ----------------------
+!!$!
+!!$DQ(IIB:IIE,:,:) = PQ(IIB+1:IIE+1,:,:) - PQ(IIB-1:IIE-1,:,:)
+!!$DQ(IIB-1,:,:) = PQ(IIB,:,:) - PQ(IIE-1,:,:)
+!!$DQ(IIE+1,:,:) = PQ(IIB+1,:,:) - PQ(IIE,:,:)
+!!$DQ = 0.5*DQ
+!!$!
+!!$END FUNCTION DIF2X
!
END FUNCTION PPM_01_X
!
@@ -650,6 +673,7 @@ USE MODD_CONF
!BEG JUAN PPM_LL
USE MODD_LUNIT
!END JUAN PPM_LL
+USE MODD_PARAMETERS, ONLY : JPHEXT
USE MODE_MPPDB
!
IMPLICIT NONE
@@ -689,6 +713,18 @@ REAL, DIMENSION(SIZE(ZCR,1),SIZE(ZCR,2),SIZE(ZCR,3)) :: ZFPOS, ZFNEG
INTEGER :: ILUOUT,IRESP ! for prints
INTEGER :: IIW,IIA
!END JUAN PPM_LL
+!
+!JUAN ACC
+INTEGER :: I,J,K ,IIU,IJU,IKU
+LOGICAL :: GSOUTH , GNORTH
+! inline shuman with macro
+#define dif2y(DQ,PQ) DQ(:,IJB:IJE,:) = 0.5*(PQ(:,IJB+1:IJE+1,:) - PQ(:,IJB-1:IJE-1,:)) ; \
+DQ(:,IJB-1,:) = 0.5*(PQ(:,IJB,:) - PQ(:,IJE-1,:)) ; \
+DQ(:,IJE+1,:) = 0.5*(PQ(:,IJB+1,:) - PQ(:,IJE,:)) ! DIF2Y(DQ,PQ)
+#define dyf(PDYF,PA) PDYF(:,1:IJU-1,:) = PA(:,2:IJU,:) - PA(:,1:IJU-1,:); PDYF(:,IJU,:) = PDYF(:,2*JPHEXT,:) ! DYF(PDYF,PA)
+#define mym(PMYM,PA) PMYM(:,2:IJU,:) = 0.5*( PA(:,2:IJU,:)+PA(:,1:IJU-1,:) ) ; PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:) ! MYM(PMYM,PA)
+!JUAN ACC
+!
!-------------------------------------------------------------------------------
!
!* 0.3. COMPUTES THE DOMAIN DIMENSIONS
@@ -697,6 +733,14 @@ INTEGER :: IIW,IIA
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IIW=IIB
IIA=IIE
+
+IIU=size(psrc,1)
+IJU=size(psrc,2)
+IKU=size(psrc,3)
+
+GSOUTH=LSOUTH_ll()
+GNORTH=LNORTH_ll()
+
IF(NHALO /= 1) THEN
CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
WRITE(ILUOUT,*) 'ERROR : PPM ppm_met.f90 --> Juan '
@@ -706,13 +750,14 @@ IF(NHALO /= 1) THEN
CALL ABORT
STOP
ENDIF
-!
CALL GET_HALO(PSRC,HDIR="01_Y")
-!
+
!
!-------------------------------------------------------------------------------
!
!
+!$acc data region local (ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG) copyin (psrc,zcr,prho) copyout(pr) !
+!$acc region ! local (ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG) copyin (psrc,zcr,prho) copyout(pr)
PR=PSRC
ZQL=PSRC
ZQR=PSRC
@@ -724,193 +769,211 @@ ZQR0=PSRC
ZQ60=PSRC
ZFPOS=PSRC
ZFNEG=PSRC
+!#acc end region
!
-SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
-!
-!* 2.1 CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
-! ---------------------------------------------
-!
-CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
-!
-! calculate dmq
- ZDMQ = DIF2Y(PSRC)
-!
-! monotonize the difference followinq eq. 5 in Lin94
-!BEG JUAN PPM_LL01
-!
-! ZDMQ(j) = Fct[ ZDMQ(j),PSRC(j),PSRC(j-1),PSRC(j+1) ]
-!
- ZDMQ(IIW:IIA,IJB:IJE,:) = &
- SIGN( (MIN( ABS(ZDMQ(IIW:IIA,IJB:IJE,:)),2.0*(PSRC(IIW:IIA,IJB:IJE,:) - &
- MIN(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:))), &
- 2.0*(MAX(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:)) - &
- PSRC(IIW:IIA,IJB:IJE,:)) )), ZDMQ(IIW:IIA,IJB:IJE,:) )
-!
-! SOUTH BOUND
-!
-!!$ ZDMQ(:,IJB-1,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(:,IJB-1,:)), 2.0*(PSRC(:,IJB-1,:) - &
-!!$ MIN(PSRC(:,IJE-1,:),PSRC(:,IJB-1,:),PSRC(:,IJB,:))), &
-!!$ 2.0*(MAX(PSRC(:,IJE-1,:),PSRC(:,IJB-1,:),PSRC(:,IJB,:)) - &
-!!$ PSRC(:,IJB-1,:)) )), ZDMQ(:,IJB-1,:) )
-!
-! NORTH BOUND
-!
-!!$ ZDMQ(:,IJE+1,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(:,IJE+1,:)), 2.0*(PSRC(:,IJE+1,:) - &
-!!$ MIN(PSRC(:,IJE,:),PSRC(:,IJE+1,:),PSRC(:,IJB+1,:))), &
-!!$ 2.0*(MAX(PSRC(:,IJE,:),PSRC(:,IJE+1,:),PSRC(:,IJB+1,:)) - &
-!!$ PSRC(:,IJE+1,:)) )), ZDMQ(:,IJE+1,:) )
-!
-! update ZDMQ HALO before next/further utilisation
-!
- CALL GET_HALO(ZDMQ,HDIR="01_Y")
- CALL MPPDB_CHECK3DM("PPM::PPM_01_Y CYCL ::ZDMQ",PRECISION,ZDMQ)
-!
-! calculate qL and qR with the modified dmq
-!
- ZQL0(IIW:IIA,IJB:IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJB:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE,:)) - &
- (ZDMQ(IIW:IIA,IJB:IJE+1,:) - ZDMQ(IIW:IIA,IJB-1:IJE,:))/6.0
-!
- CALL GET_HALO(ZQL0,HDIR="01_Y")
-!
-! SOUTH BOUND
-!
-!!$ ZQL0(:,IJB-1,:) = ZQL0(:,IJE,:) JUAN PPMLL01
-!
- ZQR0(IIW:IIA,IJB-1:IJE,:) = ZQL0(IIW:IIA,IJB:IJE+1,:)
-!
- CALL GET_HALO(ZQR0,HDIR="01_Y")
-!
-! NORTH BOUND
-!
-!!$ ZQR0(:,IJE+1,:) = ZQR0(:,IJB,:) JUAN PPMLL01
-!
-! determine initial coefficients of the parabolae
-!
- ZDQ = ZQR0 - ZQL0
- ZQ60 = 6.0*(PSRC - 0.5*(ZQL0 + ZQR0))
-!
-! initialize final parabolae parameters
-!
- ZQL = ZQL0
- ZQR = ZQR0
- ZQ6 = ZQ60
-!
-! eliminate over and undershoots and create qL and qR as in Lin96
-!
- WHERE ( ZDMQ == 0.0 )
- ZQL = PSRC
- ZQR = PSRC
- ZQ6 = 0.0
- ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 )
- ZQ6 = 3.0*(ZQL0 - PSRC)
- ZQR = ZQL0 - ZQ6
- ZQL = ZQL0
- ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
- ZQ6 = 3.0*(ZQR0 - PSRC)
- ZQL = ZQR0 - ZQ6
- ZQR = ZQR0
- END WHERE
-!
-! recalculate coefficients of the parabolae
-!
- ZDQ = ZQR - ZQL
-!
-! and finally calculate fluxes for the advection
-!
-! ZFPOS(j) = Fct[ ZQR(j-1),ZCR(i),ZDQ(j-1),ZQ6(j-1) ]
-!
- ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZQR(IIW:IIA,IJB-1:IJE,:) - 0.5*ZCR(IIW:IIA,IJB:IJE+1,:) * &
- (ZDQ(IIW:IIA,IJB-1:IJE,:) - (1.0 - 2.0*ZCR(IIW:IIA,IJB:IJE+1,:)/3.0) &
- * ZQ6(IIW:IIA,IJB-1:IJE,:))
-!
- CALL GET_HALO(ZFPOS,HDIR="01_Y")
-!
-! SOUTH BOUND
-!
-! PPOSX(:,IJB-1,:) is not important for the calc of advection so
-! we set it to 0
-!!$ ZFPOS(:,IJB-1,:) = 0.0 JUANPPMLL01
-!
- ZFNEG(IIW:IIA,:,:) = ZQL(IIW:IIA,:,:) - 0.5*ZCR(IIW:IIA,:,:) * &
- ( ZDQ(IIW:IIA,:,:) + (1.0 + 2.0*ZCR(IIW:IIA,:,:)/3.0) * ZQ6(IIW:IIA,:,:) )
-!
- CALL GET_HALO(ZFNEG,HDIR="01_Y")
-!
-! advect the actual field in Y direction by V*dt
-!
- PR = DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
- CALL GET_HALO(PR,HDIR="01_Y")
- CALL MPPDB_CHECK3DM("PPM::PPM_01_Y CYCL ::PR",PRECISION,PR)
-!
-!* 2.2 NON-CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
-! -------------------------------------------------
-!
-CASE('OPEN')
+!!$SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
+!!$!
+!!$!* 2.1 CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
+!!$! ---------------------------------------------
+!!$!
+!!$CASE ('CYCL') ! In that case one must have HLBCY(1) == HLBCY(2)
+!!$!
+!!$! calculate dmq
+!!$!! ZDMQ = DIF2Y(PSRC)
+!!$ dif2y(ZDMQ,PSRC)
+!!$!
+!!$! monotonize the difference followinq eq. 5 in Lin94
+!!$!BEG JUAN PPM_LL01
+!!$!
+!!$! ZDMQ(j) = Fct[ ZDMQ(j),PSRC(j),PSRC(j-1),PSRC(j+1) ]
+!!$!
+!!$ ZDMQ(IIW:IIA,IJB:IJE,:) = &
+!!$ SIGN( (MIN( ABS(ZDMQ(IIW:IIA,IJB:IJE,:)),2.0*(PSRC(IIW:IIA,IJB:IJE,:) - &
+!!$ MIN(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:))), &
+!!$ 2.0*(MAX(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:)) - &
+!!$ PSRC(IIW:IIA,IJB:IJE,:)) )), ZDMQ(IIW:IIA,IJB:IJE,:) )
+!!$
+!!$! update ZDMQ HALO before next/further utilisation
+!!$!
+!!$ CALL GET_HALO(ZDMQ,HDIR="01_Y")
+!!$ CALL MPPDB_CHECK3DM("PPM::PPM_01_Y CYCL ::ZDMQ",PRECISION,ZDMQ)
+!!$!
+!!$! calculate qL and qR with the modified dmq
+!!$!
+!!$ ZQL0(IIW:IIA,IJB:IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJB:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE,:)) - &
+!!$ (ZDMQ(IIW:IIA,IJB:IJE+1,:) - ZDMQ(IIW:IIA,IJB-1:IJE,:))/6.0
+!!$!
+!!$ CALL GET_HALO(ZQL0,HDIR="01_Y")
+!!$!
+!!$! SOUTH BOUND
+!!$!
+!!$ ZQR0(IIW:IIA,IJB-1:IJE,:) = ZQL0(IIW:IIA,IJB:IJE+1,:)
+!!$!
+!!$ CALL GET_HALO(ZQR0,HDIR="01_Y")
+!!$!
+!!$! determine initial coefficients of the parabolae
+!!$!
+!!$ ZDQ = ZQR0 - ZQL0
+!!$ ZQ60 = 6.0*(PSRC - 0.5*(ZQL0 + ZQR0))
+!!$!
+!!$! initialize final parabolae parameters
+!!$!
+!!$ ZQL = ZQL0
+!!$ ZQR = ZQR0
+!!$ ZQ6 = ZQ60
+!!$!
+!!$! eliminate over and undershoots and create qL and qR as in Lin96
+!!$!
+!!$ WHERE ( ZDMQ == 0.0 )
+!!$ ZQL = PSRC
+!!$ ZQR = PSRC
+!!$ ZQ6 = 0.0
+!!$ ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 )
+!!$ ZQ6 = 3.0*(ZQL0 - PSRC)
+!!$ ZQR = ZQL0 - ZQ6
+!!$ ZQL = ZQL0
+!!$ ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
+!!$ ZQ6 = 3.0*(ZQR0 - PSRC)
+!!$ ZQL = ZQR0 - ZQ6
+!!$ ZQR = ZQR0
+!!$ END WHERE
+!!$!
+!!$! recalculate coefficients of the parabolae
+!!$!
+!!$ ZDQ = ZQR - ZQL
+!!$!
+!!$! and finally calculate fluxes for the advection
+!!$!
+!!$! ZFPOS(j) = Fct[ ZQR(j-1),ZCR(i),ZDQ(j-1),ZQ6(j-1) ]
+!!$!
+!!$ ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZQR(IIW:IIA,IJB-1:IJE,:) - 0.5*ZCR(IIW:IIA,IJB:IJE+1,:) * &
+!!$ (ZDQ(IIW:IIA,IJB-1:IJE,:) - (1.0 - 2.0*ZCR(IIW:IIA,IJB:IJE+1,:)/3.0) &
+!!$ * ZQ6(IIW:IIA,IJB-1:IJE,:))
+!!$!
+!!$ CALL GET_HALO(ZFPOS,HDIR="01_Y")
+!!$!
+!!$! SOUTH BOUND
+!!$!
+!!$! PPOSX(:,IJB-1,:) is not important for the calc of advection so
+!!$! we set it to 0
+!!$!
+!!$ ZFNEG(IIW:IIA,:,:) = ZQL(IIW:IIA,:,:) - 0.5*ZCR(IIW:IIA,:,:) * &
+!!$ ( ZDQ(IIW:IIA,:,:) + (1.0 + 2.0*ZCR(IIW:IIA,:,:)/3.0) * ZQ6(IIW:IIA,:,:) )
+!!$!
+!!$ CALL GET_HALO(ZFNEG,HDIR="01_Y")
+!!$!
+!!$! advect the actual field in Y direction by V*dt
+!!$!
+!!$ PR = DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+!!$ CALL GET_HALO(PR,HDIR="01_Y")
+!!$ CALL MPPDB_CHECK3DM("PPM::PPM_01_Y CYCL ::PR",PRECISION,PR)
+
+
+!!$!
+!!$!* 2.2 NON-CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
+!!$! -------------------------------------------------
+!!$!
+!!$CASE('OPEN')
!
! calculate dmq
- ZDMQ = DIF2Y(PSRC)
+!! ZDMQ = DIF2Y(PSRC)
+ !#acc region
+ dif2y(ZDMQ,PSRC)
+ !#acc end region
+
! overwrite the values on the boundary to get second order difference
! for qL and qR at the boundary
!
! SOUTH BOUND
!
- IF (LSOUTH_ll()) THEN
+!!$ IF (LSOUTH_ll()) THEN
+ IF (GSOUTH) THEN
+ !#acc region
ZDMQ(IIW:IIA,IJB-1,:) = -ZDMQ(IIW:IIA,IJB,:)
+ !#acc end region
ENDIF
!
! NORTH BOUND
!
- IF (LNORTH_ll()) THEN
+ IF (GNORTH) THEN
+ !#acc region
ZDMQ(IIW:IIA,IJE+1,:) = -ZDMQ(IIW:IIA,IJE,:)
+ !#acc end region
ENDIF
+
!
! monotonize the difference followinq eq. 5 in Lin94
+ !#acc region
ZDMQ(IIW:IIA,IJB:IJE,:) = &
SIGN( (MIN( ABS(ZDMQ(IIW:IIA,IJB:IJE,:)),2.0*(PSRC(IIW:IIA,IJB:IJE,:) - &
MIN(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:))), &
2.0*(MAX(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:)) - &
PSRC(IIW:IIA,IJB:IJE,:)) )), ZDMQ(IIW:IIA,IJB:IJE,:) )
-!!$ ZDMQ(:,IJB-1,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(:,IJB-1,:)), 2.0*(PSRC(:,IJB-1,:) - &
-!!$ MIN(PSRC(:,IJE-1,:),PSRC(:,IJB-1,:),PSRC(:,IJB,:))), &
-!!$ 2.0*(MAX(PSRC(:,IJE-1,:),PSRC(:,IJB-1,:),PSRC(:,IJB,:)) - &
-!!$ PSRC(:,IJB-1,:)) )), ZDMQ(:,IJB-1,:) )
-!!$ ZDMQ(:,IJE+1,:) = &
-!!$ SIGN( (MIN( ABS(ZDMQ(:,IJE+1,:)), 2.0*(PSRC(:,IJE+1,:) - &
-!!$ MIN(PSRC(:,IJE,:),PSRC(:,IJE+1,:),PSRC(:,IJB+1,:))), &
-!!$ 2.0*(MAX(PSRC(:,IJE,:),PSRC(:,IJE+1,:),PSRC(:,IJB+1,:)) - &
-!!$ PSRC(:,IJE+1,:)) )), ZDMQ(:,IJE+1,:) )
+ !#acc end region
+
+!!$ !#acc region
+!!$ DO K=1,IKU ; DO J=IJB,IJE
+!!$ !#acc do independent
+!!$ DO I=IIW,IIA ; ZDMQ(I,J,K) = SIGN( (MIN( ABS(ZDMQ(I,J,K)),2.0*(PSRC(I,J,K) - &
+!!$ MIN(PSRC(I,J-1,K),PSRC(I,J,K),PSRC(I,J+1,K))), &
+!!$ 2.0*(MAX(PSRC(I,J-1,K),PSRC(I,J,K),PSRC(I,J+1,K)) - &
+!!$ PSRC(I,J,K)) )), ZDMQ(I,J,K) )
+!!$ ENDDO ; ENDDO ; ENDDO
+!!$ !#acc end region
+
!
! update ZDMQ HALO before next/further utilisation
!
- CALL GET_HALO(ZDMQ,HDIR="01_Y")
+ !TEMPO_JUAN CALL GET_HALO(ZDMQ,HDIR="01_Y")
!
! calculate qL and qR with the modified dmq
!
+ !#acc region
ZQL0(IIW:IIA,IJB:IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJB:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE,:)) - &
(ZDMQ(IIW:IIA,IJB:IJE+1,:) - ZDMQ(IIW:IIA,IJB-1:IJE,:))/6.0
-!
- CALL GET_HALO(ZQL0,HDIR="01_Y")
+ !#acc end region
+
+!!$ !#acc region
+!!$ DO K=1,IKU ; DO J=IJB,IJE+1
+!!$ !#acc do independent
+!!$ DO I=IIW,IIA ; ZQL0(I,J,K) = 0.5*(PSRC(I,J,K) + PSRC(I,J-1,K)) - (ZDMQ(I,J,K) - ZDMQ(I,J-1,K))/6.0
+!!$ ENDDO ; ENDDO ; ENDDO
+!!$ !#acc end region
+!
+ !TEMPO_JUAN CALL GET_HALO(ZQL0,HDIR="01_Y")
!
! SOUTH BOUND
!
- IF (LSOUTH_ll()) THEN
+!!$ IF (LSOUTH_ll()) THEN
+ IF (GSOUTH) THEN
+ !#acc region
ZQL0(IIW:IIA,IJB-1,:) = ZQL0(IIW:IIA,IJB,:)
+ !#acc end region
ENDIF
!
+ !#acc region
ZQR0(IIW:IIA,IJB-1:IJE,:) = ZQL0(IIW:IIA,IJB:IJE+1,:)
+ !#acc end region
+
+!!$ !#acc region
+!!$ DO K=1,IKU ; DO J=IJB-1,IJE
+!!$ !#acc do independent
+!!$ DO I=IIW,IIA ; ZQR0(I,J,K) = ZQL0(I,J+1,K)
+!!$ ENDDO ; ENDDO ; ENDDO
+!!$ !#acc end region
!
! NORTH BOUND
!
- IF (LNORTH_ll()) THEN
+ IF (GNORTH) THEN
+ !#acc region
ZQR0(IIW:IIA,IJE+1,:) = ZQR0(IIW:IIA,IJE,:)
+ !#acc end region
ENDIF
!
! determine initial coefficients of the parabolae
!
+ !#acc region
ZDQ = ZQR0 - ZQL0
ZQ60 = 6.0*(PSRC - 0.5*(ZQL0 + ZQR0))
!
@@ -919,135 +982,176 @@ CASE('OPEN')
ZQL = ZQL0
ZQR = ZQR0
ZQ6 = ZQ60
+ !#acc end region
!
! eliminate over and undershoots and create qL and qR as in Lin96
!
+ ! acc region
WHERE ( ZDMQ == 0.0 )
ZQL = PSRC
ZQR = PSRC
ZQ6 = 0.0
- ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 )
+ ENDWHERE
+ WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -(ZDQ)**2 ) )
ZQ6 = 3.0*(ZQL0 - PSRC)
ZQR = ZQL0 - ZQ6
ZQL = ZQL0
- ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
+ ENDWHERE
+ WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > (ZDQ)**2 ) )
ZQ6 = 3.0*(ZQR0 - PSRC)
ZQL = ZQR0 - ZQ6
ZQR = ZQR0
END WHERE
+ ! acc end region
+
!
! recalculate coefficients of the parabolae
!
+ !#acc region
ZDQ = ZQR - ZQL
+ !#acc end region
!
-! and finally calculate fluxes for the advection
-!!$ ZFPOS(:,IJB+1:IJE+1,:) = ZQR(:,IJB:IJE,:) - 0.5*ZCR(:,IJB+1:IJE+1,:) * &
-!!$ (ZDQ(:,IJB:IJE,:) - (1.0 - 2.0*ZCR(:,IJB+1:IJE+1,:)/3.0) &
-!!$ * ZQ6(:,IJB:IJE,:))
+ ! and finally calculate fluxes for the advection
+ !#acc region
ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZQR(IIW:IIA,IJB-1:IJE,:) - 0.5*ZCR(IIW:IIA,IJB:IJE+1,:) * &
(ZDQ(IIW:IIA,IJB-1:IJE,:) - (1.0 - 2.0*ZCR(IIW:IIA,IJB:IJE+1,:)/3.0) &
* ZQ6(IIW:IIA,IJB-1:IJE,:))
+ !#acc end region
+
+!!$ !#acc region
+!!$ DO K=1,IKU ; DO J=IJB,IJE+1
+!!$ !#acc do independent
+!!$ DO I=IIW,IIA ; ZFPOS(I,J,K) = ZQR(I,J-1,K) - 0.5*ZCR(I,J,K) * &
+!!$ (ZDQ(I,J-1,K) - (1.0 - 2.0*ZCR(I,J,K)/3.0) &
+!!$ * ZQ6(I,J-1,K))
+!!$ ENDDO ; ENDDO ; ENDDO
+!!$ !#acc end region
!
- CALL GET_HALO(ZFPOS,HDIR="01_Y")
+ !TEMPO_JUAN CALL GET_HALO(ZFPOS,HDIR="01_Y")
!
!
! advection flux at open boundary when u(IJB) > 0
!
! SOUTH BOUND
!
- IF (LSOUTH_ll()) THEN
+!!$ IF (LSOUTH_ll()) THEN
+ IF (GSOUTH) THEN
+ !#acc region
ZFPOS(IIW:IIA,IJB,:) = (PSRC(IIW:IIA,IJB-1,:) - ZQR(IIW:IIA,IJB-1,:))*ZCR(IIW:IIA,IJB,:) + &
ZQR(IIW:IIA,IJB-1,:)
+ !#acc end region
ENDIF
!
! PPOSX(:,IJB-1,:) is not important for the calc of advection so
! we set it to 0
-!!$ ZFPOS(:,IJB-1,:) = 0.0 ! JUAN PPMLL01
-!
-!!$ ZFNEG(:,IJB-1:IJE,:) = ZQL(:,IJB-1:IJE,:) - 0.5*ZCR(:,IJB-1:IJE,:) * &
-!!$ ( ZDQ(:,IJB-1:IJE,:) + (1.0 + 2.0*ZCR(:,IJB-1:IJE,:)/3.0) * &
-!!$ ZQ6(:,IJB-1:IJE,:) )
+ !#acc region
ZFNEG(IIW:IIA,:,:) = ZQL(IIW:IIA,:,:) - 0.5*ZCR(IIW:IIA,:,:) * &
( ZDQ(IIW:IIA,:,:) + (1.0 + 2.0*ZCR(IIW:IIA,:,:)/3.0) * ZQ6(IIW:IIA,:,:) )
+ !#acc end region
+
+!!$ !#acc region
+!!$ DO K=1,IKU ; DO J=1,IJU
+!!$ !#acc do independent
+!!$ DO I=IIW,IIA ; ZFNEG(I,J,K) = ZQL(I,J,K) - 0.5*ZCR(I,J,K) * &
+!!$ ( ZDQ(I,J,K) + (1.0 + 2.0*ZCR(I,J,K)/3.0) * ZQ6(I,J,K) )
+!!$ ENDDO ; ENDDO ; ENDDO
+!!$ !#acc end region
!
- CALL GET_HALO(ZFNEG,HDIR="01_Y")
-!
-! advection flux at open boundary when u(IJE+1) < 0
-!
-! NORTH BOUND
-!
- IF (LNORTH_ll()) THEN
- ZFNEG(IIW:IIA,IJE+1,:) = (ZQR(IIW:IIA,IJE,:)-PSRC(IIW:IIA,IJE+1,:))*ZCR(IIW:IIA,IJE+1,:) + &
- ZQR(IIW:IIA,IJE,:)
- ENDIF
-!
-! advect the actual field in X direction by U*dt
-!
- PR = DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
-!
- CALL GET_HALO(PR,HDIR="01_Y")
-!
-!
-END SELECT
-!
-CONTAINS
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-! ########################################################################
- FUNCTION DIF2Y(PQ) RESULT(DQ)
-! ########################################################################
-!!
-!!**** DIF2Y - leap-frog difference operator in Y direction
-!!
-!! Calculates the difference assuming periodic BC (CYCL).
-!!
-!! DQ(J) = 0.5 * (PQ(J+1) - PQ(J-1))
-!!
-!!
-!! MODIFICATIONS
-!! -------------
-!!
-!! 18.3.2006. T. Maric - original version, works only for periodic boundary
-!! conditions and on one domain
-!!
-!-------------------------------------------------------------------------------
-!
-!
-USE MODE_ll
-!
-IMPLICIT NONE
-!
-!* 0.1 Declarations of dummy arguments :
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
-REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
-!
-!* 0.2 Declarations of local variables :
-!
-INTEGER :: IIB,IJB ! Begining useful area in x,y directions
-INTEGER :: IIE,IJE ! End useful area in x,y directions
-!
-!-------------------------------------------------------------------------------
-!
-!* 1.0. COMPUTE THE DOMAIN DIMENSIONS
-! -----------------------------
+ !TEMPO_JUAN CALL GET_HALO(ZFNEG,HDIR="01_Y")
!
-CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+! advection flux at open boundary when u(IJE+1) < 0
!
-!-------------------------------------------------------------------------------
+! NORTH BOUND
!
-!* 2.0. COMPUTE THE DIFFERENCE
-! ----------------------
+ IF (GNORTH) THEN
+ !#acc region
+ ZFNEG(IIW:IIA,IJE+1,:) = (ZQR(IIW:IIA,IJE,:)-PSRC(IIW:IIA,IJE+1,:))*ZCR(IIW:IIA,IJE+1,:) + &
+ ZQR(IIW:IIA,IJE,:)
+ !#acc end region
+ ENDIF
+!!$!
+!!$! advect the actual field in X direction by U*dt
+!!$!
+ !#acc region
+ mym(ZQL,PRHO)
+ ZQR = ZCR* ZQL*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + ZFNEG*(0.5-SIGN(0.5,ZCR)) )
+ dyf(PR,ZQR)
+ !#acc end region
+
+!!$ PR = ZDMQ + ZQL0 + ZDQ + ZQ60 + ZQL + ZQR + ZQ6 + ZQR0 + ZFPOS + ZFNEG
+
+!!$ PR = DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+!!$!
+
+!$acc end region
+! acc update host (pr)
+!$acc end data region
+
+ CALL GET_HALO(PR,HDIR="01_Y")
+
!
-DQ(:,IJB:IJE,:) = PQ(:,IJB+1:IJE+1,:) - PQ(:,IJB-1:IJE-1,:)
-DQ(:,IJB-1,:) = PQ(:,IJB,:) - PQ(:,IJE-1,:)
-DQ(:,IJE+1,:) = PQ(:,IJB+1,:) - PQ(:,IJE,:)
-DQ = 0.5 * DQ
+!!$END SELECT
+
+
!
-END FUNCTION DIF2Y
+!!$CONTAINS
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$! ########################################################################
+!!$ FUNCTION DIF2Y(PQ) RESULT(DQ)
+!!$! ########################################################################
+!!$!!
+!!$!!**** DIF2Y - leap-frog difference operator in Y direction
+!!$!!
+!!$!! Calculates the difference assuming periodic BC (CYCL).
+!!$!!
+!!$!! DQ(J) = 0.5 * (PQ(J+1) - PQ(J-1))
+!!$!!
+!!$!!
+!!$!! MODIFICATIONS
+!!$!! -------------
+!!$!!
+!!$!! 18.3.2006. T. Maric - original version, works only for periodic boundary
+!!$!! conditions and on one domain
+!!$!!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!
+!!$USE MODE_ll
+!!$!
+!!$IMPLICIT NONE
+!!$!
+!!$!* 0.1 Declarations of dummy arguments :
+!!$!
+!!$REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
+!!$REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
+!!$!
+!!$!* 0.2 Declarations of local variables :
+!!$!
+!!$INTEGER :: IIB,IJB ! Begining useful area in x,y directions
+!!$INTEGER :: IIE,IJE ! End useful area in x,y directions
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!* 1.0. COMPUTE THE DOMAIN DIMENSIONS
+!!$! -----------------------------
+!!$!
+!!$CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!* 2.0. COMPUTE THE DIFFERENCE
+!!$! ----------------------
+!!$!
+!!$DQ(:,IJB:IJE,:) = PQ(:,IJB+1:IJE+1,:) - PQ(:,IJB-1:IJE-1,:)
+!!$DQ(:,IJB-1,:) = PQ(:,IJB,:) - PQ(:,IJE-1,:)
+!!$DQ(:,IJE+1,:) = PQ(:,IJB+1,:) - PQ(:,IJE,:)
+!!$DQ = 0.5 * DQ
+!!$!
+!!$END FUNCTION DIF2Y
!
END FUNCTION PPM_01_Y
!
@@ -1109,6 +1213,17 @@ REAL, DIMENSION(SIZE(ZCR,1),SIZE(ZCR,2),SIZE(ZCR,3)) :: ZQL0,ZQR0,ZQ60
! advection fluxes
REAL, DIMENSION(SIZE(ZCR,1),SIZE(ZCR,2),SIZE(ZCR,3)) :: ZFPOS, ZFNEG
!
+!JUAN ACC
+INTEGER :: I,J,K ,IIU,IJU,IKU
+LOGICAL :: GWEST , GEAST
+! inline shuman with macro
+#define dif2z(DQ,PQ) DQ(:,:,IKB:IKE) = 0.5*(PQ(:,:,IKB+1:IKE+1) - PQ(:,:,IKB-1:IKE-1)) ; \
+DQ(:,:,IKB-1) = -DQ(:,:,IKB) ;\
+DQ(:,:,IKE+1) = -DQ(:,:,IKE) ! DIF2Z(DQ,PQ)
+#define dzf(PDZF,PA) PDZF(:,:,1:IKU-1) = PA(:,:,2:IKU) - PA(:,:,1:IKU-1) ; PDZF(:,:,IKU) = -999. ! DZF(PDZF,PA)
+#define mzm(PMZM,PA) PMZM(:,:,2:IKU) = 0.5*( PA(:,:,2:IKU)+PA(:,:,1:IKU-1) ) ; PMZM(:,:,1) = -999. ! MZM(PMZM,PA)
+!JUAN ACC
+!
!-------------------------------------------------------------------------------
!
!* 0.3. COMPUTES THE DOMAIN DIMENSIONS
@@ -1116,6 +1231,14 @@ REAL, DIMENSION(SIZE(ZCR,1),SIZE(ZCR,2),SIZE(ZCR,3)) :: ZFPOS, ZFNEG
!
IKB = 1 + JPVEXT
IKE = SIZE(PSRC,3) - JPVEXT
+
+IIU=size(psrc,1)
+IJU=size(psrc,2)
+IKU=size(psrc,3)
+
+
+!$acc data region local (ZDMQ,ZQL0,ZQR0,ZDQ,ZQ60,ZQL,ZQR,ZQ6,ZFPOS,ZFNEG) copyin (psrc,zcr,prho) copyout(pr)
+!$acc region
!
!-------------------------------------------------------------------------------
!
@@ -1123,7 +1246,8 @@ IKE = SIZE(PSRC,3) - JPVEXT
! --------------------------------
!
! calculate dmq
-ZDMQ = DIF2Z(PSRC)
+!! ZDMQ = DIF2Z(PSRC)
+ dif2z(ZDMQ,PSRC)
!
! monotonize the difference followinq eq. 5 in Lin94
! use the periodic BC here, it doesn't matter for vertical (hopefully)
@@ -1166,19 +1290,21 @@ ZQ6 = ZQ60
!
! eliminate over and undershoots and create qL and qR as in Lin96
!
-WHERE ( ZDMQ == 0.0 )
- ZQL = PSRC
- ZQR = PSRC
- ZQ6 = 0.0
-ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 )
- ZQ6 = 3.0*(ZQL0 - PSRC)
- ZQR = ZQL0 - ZQ6
- ZQL = ZQL0
-ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
- ZQ6 = 3.0*(ZQR0 - PSRC)
- ZQL = ZQR0 - ZQ6
- ZQR = ZQR0
-END WHERE
+ WHERE ( ZDMQ == 0.0 )
+ ZQL = PSRC
+ ZQR = PSRC
+ ZQ6 = 0.0
+ END WHERE
+ WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -(ZDQ)**2 ) )
+ ZQ6 = 3.0*(ZQL0 - PSRC)
+ ZQR = ZQL0 - ZQ6
+ ZQL = ZQL0
+ END WHERE
+ WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > (ZDQ)**2 ) )
+ ZQ6 = 3.0*(ZQR0 - PSRC)
+ ZQL = ZQR0 - ZQ6
+ ZQR = ZQR0
+ END WHERE
!
! recalculate coefficients of the parabolae
!
@@ -1208,71 +1334,82 @@ ZFNEG(:,:,IKE+1) = (ZQR(:,:,IKE)-PSRC(:,:,IKE+1))*ZCR(:,:,IKE+1) + &
!
! advect the actual field in Z direction by W*dt
!
-PR = DZF( ZCR*MZM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+
+ mzm(ZQL,PRHO)
+ ZQR = ZCR* ZQL*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + ZFNEG*(0.5-SIGN(0.5,ZCR)) )
+ dzf(PR,ZQR)
+
+!!$PR = DZF( ZCR*MZM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+
+
+!$acc end region
+! acc update host (pr)
+!$acc end data region
+
CALL GET_HALO(PR)
CALL MPPDB_CHECK3DM("PPM::PPM_01_Z ::PR",PRECISION,PR)
!
-CONTAINS
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-! ########################################################################
- FUNCTION DIF2Z(PQ) RESULT(DQ)
-! ########################################################################
-!!
-!!**** DIF2Z - leap-frog difference operator in Z direction
-!!
-!! Calculates the difference assuming periodic BC (CYCL).
-!!
-!! DQ(K) = 0.5 * (PQ(K+1) - PQ(K-1))
-!!
-!!
-!! MODIFICATIONS
-!! -------------
-!!
-!! 18.3.2006. T. Maric - original version
-!!
-!-------------------------------------------------------------------------------
-!
-!
-USE MODE_ll
-USE MODD_CONF
-USE MODD_PARAMETERS
-!
-IMPLICIT NONE
-!
-!* 0.1 Declarations of dummy arguments :
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
-REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
-!
-!* 0.2 Declarations of local variables :
-!
-INTEGER :: IKB ! Begining useful area in z directions
-INTEGER :: IKE ! End useful area in z directions
-!
-!-------------------------------------------------------------------------------
-!
-!* 1.0. COMPUTE THE DOMAIN DIMENSIONS
-! -----------------------------
-!
-!CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
-IKB = 1 + JPVEXT
-IKE = SIZE(PQ,3) - JPVEXT
-!
-!-------------------------------------------------------------------------------
-!
-!* 2.0. COMPUTE THE DIFFERENCE
-! ----------------------
-!
-DQ(:,:,IKB:IKE) = PQ(:,:,IKB+1:IKE+1) - PQ(:,:,IKB-1:IKE-1)
-DQ(:,:,IKB-1) = -DQ(:,:,IKB)
-DQ(:,:,IKE+1) = -DQ(:,:,IKE)
-DQ = 0.5 * DQ
-!
-END FUNCTION DIF2Z
+!!$CONTAINS
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$! ########################################################################
+!!$ FUNCTION DIF2Z(PQ) RESULT(DQ)
+!!$! ########################################################################
+!!$!!
+!!$!!**** DIF2Z - leap-frog difference operator in Z direction
+!!$!!
+!!$!! Calculates the difference assuming periodic BC (CYCL).
+!!$!!
+!!$!! DQ(K) = 0.5 * (PQ(K+1) - PQ(K-1))
+!!$!!
+!!$!!
+!!$!! MODIFICATIONS
+!!$!! -------------
+!!$!!
+!!$!! 18.3.2006. T. Maric - original version
+!!$!!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!
+!!$USE MODE_ll
+!!$USE MODD_CONF
+!!$USE MODD_PARAMETERS
+!!$!
+!!$IMPLICIT NONE
+!!$!
+!!$!* 0.1 Declarations of dummy arguments :
+!!$!
+!!$REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
+!!$REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
+!!$!
+!!$!* 0.2 Declarations of local variables :
+!!$!
+!!$INTEGER :: IKB ! Begining useful area in z directions
+!!$INTEGER :: IKE ! End useful area in z directions
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!* 1.0. COMPUTE THE DOMAIN DIMENSIONS
+!!$! -----------------------------
+!!$!
+!!$!CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+!!$IKB = 1 + JPVEXT
+!!$IKE = SIZE(PQ,3) - JPVEXT
+!!$!
+!!$!-------------------------------------------------------------------------------
+!!$!
+!!$!* 2.0. COMPUTE THE DIFFERENCE
+!!$! ----------------------
+!!$!
+!!$DQ(:,:,IKB:IKE) = PQ(:,:,IKB+1:IKE+1) - PQ(:,:,IKB-1:IKE-1)
+!!$DQ(:,:,IKB-1) = -DQ(:,:,IKB)
+!!$DQ(:,:,IKE+1) = -DQ(:,:,IKE)
+!!$DQ = 0.5 * DQ
+!!$!
+!!$END FUNCTION DIF2Z
!
END FUNCTION PPM_01_Z
!
@@ -1305,6 +1442,7 @@ USE MODD_CONF
USE MODD_LUNIT
USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
!END JUAN PPM_LL
+USE MODD_PARAMETERS, ONLY : JPHEXT
USE MODE_MPPDB
!
IMPLICIT NONE
@@ -1340,6 +1478,14 @@ TYPE(HALO2LIST_ll), POINTER :: TZ_PSRC_HALO2_ll ! halo2 for PSRC
INTEGER :: ILUOUT,IRESP ! for prints
INTEGER :: IJS,IJN
!END JUAN PPM_LL
+!JUAN ACC
+REAL, DIMENSION(SIZE(ZCR,1),SIZE(ZCR,2),SIZE(ZCR,3)) :: ZRHO_MXM , ZCR_MXM , ZCR_DXF
+INTEGER :: I,J,K ,IIU,IJU,IKU
+LOGICAL :: GWEST , GEAST
+! inline shuman with macro
+#define dxf(PDXF,PA) PDXF(1:IIU-1,:,:) = PA(2:IIU,:,:) - PA(1:IIU-1,:,:) ; PDXF(IIU,:,:) = PDXF(2*JPHEXT,:,:) ! DXF(PDXF,PA)
+#define mxm(PMXM,PA) PMXM(2:IIU,:,:) = 0.5*( PA(2:IIU,:,:)+PA(1:IIU-1,:,:) ) ; PMXM(1,:,:) = PMXM(IIU-2*JPHEXT+1,:,:) ! MXM(PMXM,PA)
+!JUAN ACC
!-------------------------------------------------------------------------------
!
!* 0.3. COMPUTES THE DOMAIN DIMENSIONS
@@ -1348,8 +1494,14 @@ INTEGER :: IJS,IJN
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IJS=IJB
IJN=IJE
-!!$IJS=IJB-1
-!!$IJN=IJE+1
+
+IIU=size(psrc,1)
+IJU=size(psrc,2)
+IKU=size(psrc,3)
+
+GWEST = LWEST_ll()
+GEAST = LEAST_ll()
+
!
!BEG JUAN PPM_LL
!
@@ -1366,11 +1518,12 @@ IF(NHALO /= 1) THEN
ENDIF
!
CALL GET_HALO2(PSRC,TZ_PSRC_HALO2_ll)
+!$acc data region local (PR,ZPHAT,ZFPOS,ZFNEG,ZRHO_MXM,ZCR_MXM,ZCR_DXF) copyin (psrc,zcr,prho)
+!$acc region
ZPHAT=PSRC
ZFPOS=PSRC
ZFNEG=PSRC
PR=PSRC
-!!$!
!
!END JUAN PPM_LL
!-------------------------------------------------------------------------------
@@ -1389,132 +1542,98 @@ ZPHAT(IIB+1:IIE,IJS:IJN,:) = ( 7.0 * &
( PSRC(IIB+1:IIE ,IJS:IJN,:) + PSRC(IIB :IIE-1,IJS:IJN,:) ) - &
( PSRC(IIB+2:IIE+1,IJS:IJN,:) + PSRC(IIB-1:IIE-2,IJS:IJN,:) ) ) / 12.0
!
-SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
-CASE ('CYCL','WALL') ! In that case one must have HLBCX(1) == HLBCX(2)
-!
-!!$ ZPHAT(IIB,:,:) = (7.0 * &
-!!$ (PSRC(IIB,:,:) + PSRC(IIB-1,:,:)) - &
-!!$ (PSRC(IIB+1,:,:) + PSRC(IIE-1,:,:))) / 12.0
+!!$SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
+!!$CASE ('CYCL','WALL') ! In that case one must have HLBCX(1) == HLBCX(2)
+!!$!
+!!$!
+!!$! WEST BOUND
+!!$! WEST BOUND
+!!$! i=IIB ( need halo2 psrc(iib-2) )
+!!$! ZPATH(IIB) = Fct[ PSRC(IIB),PSRC(IIB-1),PSRC(IIB+1),PSRC(IIB-2) ]
+!!$!
+!!$!
+!!$ ZPHAT(IIB ,IJS:IJN,:) = ( 7.0 * &
+!!$ ( PSRC(IIB ,IJS:IJN,:) + PSRC(IIB-1,IJS:IJN,:) ) - &
+!!$ ( PSRC(IIB+1,IJS:IJN,:) + TZ_PSRC_HALO2_ll%HALO2%WEST(IJS:IJN,:) ) ) / 12.0
+!!$! <=> WEST BOUND ( PSRC(IIB+1,IJS:IJN,:) + PSRC(IIB-2,IJS:IJN,:) ) ) / 12.0
+!!$!
+!!$! The ZPHAT(IIB-1,:,:) doesn't matter only define an realistic value
+!!$!
+!!$!
+!!$! EAST BOUND
+!!$!
+!!$! The ZPHAT(IIE+1,:,:) doesn't matter only define an realistic value
+!!$!
+!!$!
+!!$! update ZPHAT HALO before next/further utilisation
+!!$!
+!!$! /!\ not needed if ZPHAT(IIB-1) & ZPHAT(IIE+1) doen't matter ???
+!!$!
+!!$CALL GET_HALO(ZPHAT,HDIR="Z0_X")
!!$
+!!$ ZFPOS(IIB:IIE,IJS:IJN,:) = ZPHAT(IIB:IIE,IJS:IJN,:) - &
+!!$ ZCR(IIB:IIE,IJS:IJN,:)*(ZPHAT(IIB:IIE,IJS:IJN,:) - PSRC(IIB-1:IIE-1,IJS:IJN,:)) - &
+!!$ ZCR(IIB:IIE,IJS:IJN,:)*(1.0 - ZCR(IIB:IIE,IJS:IJN,:)) * &
+!!$ (ZPHAT(IIB-1:IIE-1,IJS:IJN,:) - 2.0*PSRC(IIB-1:IIE-1,IJS:IJN,:) + ZPHAT(IIB:IIE,IJS:IJN,:))
!!$!
-!!$ ZPHAT(IIE+1,:,:) = ZPHAT(IIB,:,:)
-!!$ ZPHAT(IIB-1,:,:) = ZPHAT(IIE,:,:)
-!
-! WEST BOUND
-! i=IIB ( need halo2 psrc(iib-2) )
-! ZPATH(IIB) = Fct[ PSRC(IIB),PSRC(IIB-1),PSRC(IIB+1),PSRC(IIB-2) ]
-!
- ZPHAT(IIB ,IJS:IJN,:) = ( 7.0 * &
- ( PSRC(IIB ,IJS:IJN,:) + PSRC(IIB-1,IJS:IJN,:) ) - &
- ( PSRC(IIB+1,IJS:IJN,:) + TZ_PSRC_HALO2_ll%HALO2%WEST(IJS:IJN,:) ) ) / 12.0
-! <=> WEST BOUND ( PSRC(IIB+1,IJS:IJN,:) + PSRC(IIB-2,IJS:IJN,:) ) ) / 12.0
-!
-! The ZPHAT(IIB-1,:,:) doesn't matter only define an realistic value
-!
-!!$ ZPHAT(IIB-1,:,:) = ZPHAT(IIB,:,:) ! JUANTEST1 already done at init by ZPHAT = PSRC
+!!$!
+!!$CALL GET_HALO(ZFPOS,HDIR="Z0_X") ! JUAN
+!!$!
+!!$ ZFNEG(IIB:IIE,IJS:IJN,:) = ZPHAT(IIB:IIE,IJS:IJN,:) + &
+!!$ ZCR(IIB:IIE,IJS:IJN,:)*(ZPHAT(IIB:IIE,IJS:IJN,:) - PSRC(IIB:IIE,IJS:IJN,:)) + &
+!!$ ZCR(IIB:IIE,IJS:IJN,:)*(1.0 + ZCR(IIB:IIE,IJS:IJN,:)) * &
+!!$ (ZPHAT(IIB:IIE,IJS:IJN,:) - 2.0*PSRC(IIB:IIE,IJS:IJN,:) + ZPHAT(IIB+1:IIE+1,IJS:IJN,:))
+!!$!
+!!$! define fluxes for CYCL BC outside physical domain
+!!$!
+!!$CALL GET_HALO(ZFNEG,HDIR="Z0_X") ! JUAN
+!!$!
+!!$! calculate the advection
+!!$!
+!!$ PR = PSRC * PRHO - &
+!!$ DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+!!$ CALL GET_HALO(PR,HDIR="S0_X") ! JUAN
+!!$ CALL MPPDB_CHECK3DM("PPM::PPM_S0_X CYCL ::PR",PRECISION,PR)
+!!$!
+!!$CASE ('OPEN')
!
-! EAST BOUND
!
-! The ZPHAT(IIE+1,:,:) doesn't matter only define an realistic value
+! WEST BOUND
!
-!!$ ZPHAT(IIE+1,:,:) = ZPHAT(IIE,:,:) ! JUANTEST1
+!TEMPO_JUAN IF (.NOT. GWEST) THEN
+!TEMPO_JUAN ZPHAT(IIB ,IJS:IJN,:) = ( 7.0 * &
+!TEMPO_JUAN ( PSRC(IIB ,IJS:IJN,:) + PSRC(IIB-1,IJS:IJN,:) ) - &
+!TEMPO_JUAN ( PSRC(IIB+1,IJS:IJN,:) + TZ_PSRC_HALO2_ll%HALO2%WEST(IJS:IJN,:) ) ) / 12.0
+!TEMPO_JUAN! <=> WEST BOUND ( PSRC(IIB+1,IJS:IJN,:) + PSRC(IIB-2,IJS:IJN,:) ) ) / 12.0
+!TEMPO_JUAN ENDIF
!
! update ZPHAT HALO before next/further utilisation
!
-! /!\ not needed if ZPHAT(IIB-1) & ZPHAT(IIE+1) doen't matter ???
-!
-CALL GET_HALO(ZPHAT,HDIR="Z0_X")
-!
-!!$ ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZPHAT(IIB:IIE+1,IJS:IJN,:) - &
-!!$ ZCR(IIB:IIE+1,IJS:IJN,:)*(ZPHAT(IIB:IIE+1,IJS:IJN,:) - PSRC(IIB-1:IIE,IJS:IJN,:)) - &
-!!$ ZCR(IIB:IIE+1,IJS:IJN,:)*(1.0 - ZCR(IIB:IIE+1,IJS:IJN,:)) * &
-!!$ (ZPHAT(IIB-1:IIE,IJS:IJN,:) - 2.0*PSRC(IIB-1:IIE,IJS:IJN,:) + ZPHAT(IIB:IIE+1,IJS:IJN,:))
-!
- ZFPOS(IIB:IIE,IJS:IJN,:) = ZPHAT(IIB:IIE,IJS:IJN,:) - &
- ZCR(IIB:IIE,IJS:IJN,:)*(ZPHAT(IIB:IIE,IJS:IJN,:) - PSRC(IIB-1:IIE-1,IJS:IJN,:)) - &
- ZCR(IIB:IIE,IJS:IJN,:)*(1.0 - ZCR(IIB:IIE,IJS:IJN,:)) * &
- (ZPHAT(IIB-1:IIE-1,IJS:IJN,:) - 2.0*PSRC(IIB-1:IIE-1,IJS:IJN,:) + ZPHAT(IIB:IIE,IJS:IJN,:))
-!
-!!$ ZFPOS(IIB-1,:,:) = ZFPOS(IIE,:,:) !JUAN
-!
-CALL GET_HALO(ZFPOS,HDIR="Z0_X") ! JUAN
-!
-!!$ ZFNEG(IIB-1:IIE,IJS:IJN,:) = ZPHAT(IIB-1:IIE,IJS:IJN,:) + &
-!!$ ZCR(IIB-1:IIE,IJS:IJN,:)*(ZPHAT(IIB-1:IIE,IJS:IJN,:) - PSRC(IIB-1:IIE,IJS:IJN,:)) + &
-!!$ ZCR(IIB-1:IIE,IJS:IJN,:)*(1.0 + ZCR(IIB-1:IIE,IJS:IJN,:)) * &
-!!$ (ZPHAT(IIB-1:IIE,IJS:IJN,:) - 2.0*PSRC(IIB-1:IIE,IJS:IJN,:) + ZPHAT(IIB:IIE+1,IJS:IJN,:))
-!
- ZFNEG(IIB:IIE,IJS:IJN,:) = ZPHAT(IIB:IIE,IJS:IJN,:) + &
- ZCR(IIB:IIE,IJS:IJN,:)*(ZPHAT(IIB:IIE,IJS:IJN,:) - PSRC(IIB:IIE,IJS:IJN,:)) + &
- ZCR(IIB:IIE,IJS:IJN,:)*(1.0 + ZCR(IIB:IIE,IJS:IJN,:)) * &
- (ZPHAT(IIB:IIE,IJS:IJN,:) - 2.0*PSRC(IIB:IIE,IJS:IJN,:) + ZPHAT(IIB+1:IIE+1,IJS:IJN,:))
-!
-!
-! define fluxes for CYCL BC outside physical domain
-!!$ ZFNEG(IIE+1,:,:) = ZFNEG(IIB,:,:) !JUAN
-!
-CALL GET_HALO(ZFNEG,HDIR="Z0_X") ! JUAN
-!
-! calculate the advection
-!
- PR = PSRC * PRHO - &
- DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
-!
- CALL GET_HALO(PR,HDIR="S0_X")
- CALL MPPDB_CHECK3DM("PPM::PPM_S0_X CYCL ::PR",PRECISION,PR)
-!
-CASE ('OPEN')
-!
-!!$ ZPHAT(IIB,:,:) = 0.5*(PSRC(IIB-1,:,:) + PSRC(IIB,:,:))
-!!$ ZPHAT(IIB-1,:,:) = ZPHAT(IIB,:,:) ! not used
-!!$ ZPHAT(IIE+1,:,:) = 0.5*(PSRC(IIE,:,:) + PSRC(IIE+1,:,:))
-!
-! WEST BOUND
+!TEMPO_JUAN CALL GET_HALO(ZPHATT,HDIR="Z0_X")
!
- IF (.NOT. LWEST_ll()) THEN
- ZPHAT(IIB ,IJS:IJN,:) = ( 7.0 * &
- ( PSRC(IIB ,IJS:IJN,:) + PSRC(IIB-1,IJS:IJN,:) ) - &
- ( PSRC(IIB+1,IJS:IJN,:) + TZ_PSRC_HALO2_ll%HALO2%WEST(IJS:IJN,:) ) ) / 12.0
-! <=> WEST BOUND ( PSRC(IIB+1,IJS:IJN,:) + PSRC(IIB-2,IJS:IJN,:) ) ) / 12.0
- ENDIF
-!
-CALL GET_HALO(ZPHAT,HDIR="Z0_X")
-!
- IF (LWEST_ll()) THEN
+ IF (GWEST) THEN
ZPHAT(IIB ,IJS:IJN,:) = 0.5*(PSRC(IIB-1,IJS:IJN,:) + PSRC(IIB,IJS:IJN,:))
ZPHAT(IIB-1,IJS:IJN,:) = ZPHAT(IIB,IJS:IJN,:)
ENDIF
!
! EAST BOUND
!
- IF (LEAST_ll()) THEN
+ IF (GEAST) THEN
ZPHAT(IIE+1,IJS:IJN,:) = 0.5*(PSRC(IIE,IJS:IJN,:) + PSRC(IIE+1,IJS:IJN,:))
ENDIF
!
-! update ZPHAT HALO before next/further utilisation
-!
-!!$CALL GET_HALO(ZPHAT)
-!
-!!$ ZFPOS(IIB+1:IIE+1,:,:) = ZPHAT(IIB+1:IIE+1,:,:) - &
-!!$ ZCR(IIB+1:IIE+1,:,:)*(ZPHAT(IIB+1:IIE+1,:,:) - PSRC(IIB:IIE,:,:)) - &
-!!$ ZCR(IIB+1:IIE+1,:,:)*(1.0 - ZCR(IIB+1:IIE+1,:,:)) * &
-!!$ (ZPHAT(IIB:IIE,:,:) - 2.0*PSRC(IIB:IIE,:,:) + ZPHAT(IIB+1:IIE+1,:,:))
+! update ZPHAT HALO before next/further utilisation
!
-!!$ ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZPHAT(IIB:IIE+1,IJS:IJN,:) - &
-!!$ ZCR(IIB:IIE+1,IJS:IJN,:)*(ZPHAT(IIB:IIE+1,IJS:IJN,:) - PSRC(IIB-1:IIE,IJS:IJN,:)) - &
-!!$ ZCR(IIB:IIE+1,IJS:IJN,:)*(1.0 - ZCR(IIB:IIE+1,IJS:IJN,:)) * &
-!!$ (ZPHAT(IIB-1:IIE,IJS:IJN,:) - 2.0*PSRC(IIB-1:IIE,IJS:IJN,:) + ZPHAT(IIB:IIE+1,IJS:IJN,:))
-!!$!
ZFPOS(IIB:IIE,IJS:IJN,:) = ZPHAT(IIB:IIE,IJS:IJN,:) - &
ZCR(IIB:IIE,IJS:IJN,:)*(ZPHAT(IIB:IIE,IJS:IJN,:) - PSRC(IIB-1:IIE-1,IJS:IJN,:)) - &
ZCR(IIB:IIE,IJS:IJN,:)*(1.0 - ZCR(IIB:IIE,IJS:IJN,:)) * &
(ZPHAT(IIB-1:IIE-1,IJS:IJN,:) - 2.0*PSRC(IIB-1:IIE-1,IJS:IJN,:) + ZPHAT(IIB:IIE,IJS:IJN,:))
!
-CALL GET_HALO(ZFPOS,HDIR="Z0_X") ! JUAN
+!TEMPO_JUAN CALL GET_HALO(ZFPOS,HDIR="Z0_X") ! JUAN
!
! positive flux on the WEST boundary
- IF (LWEST_ll()) THEN
+ IF (GWEST) THEN
ZFPOS(IIB,IJS:IJN,:) = (PSRC(IIB-1,IJS:IJN,:) - ZPHAT(IIB,IJS:IJN,:))*ZCR(IIB,IJS:IJN,:) + &
ZPHAT(IIB,IJS:IJN,:)
! this is not used
@@ -1522,25 +1641,15 @@ CALL GET_HALO(ZFPOS,HDIR="Z0_X") ! JUAN
ENDIF
!
! negative fluxes
-!!$ ZFNEG(IIB:IIE,:,:) = ZPHAT(IIB:IIE,:,:) + &
-!!$ ZCR(IIB:IIE,:,:)*(ZPHAT(IIB:IIE,:,:) - PSRC(IIB:IIE,:,:)) + &
-!!$ ZCR(IIB:IIE,:,:)*(1.0 + ZCR(IIB:IIE,:,:)) * &
-!!$ (ZPHAT(IIB:IIE,:,:) - 2.0*PSRC(IIB:IIE,:,:) + ZPHAT(IIB+1:IIE+1,:,:))
-!
-!!$ ZFNEG(IIB-1:IIE,IJS:IJN,:) = ZPHAT(IIB-1:IIE,IJS:IJN,:) + &
-!!$ ZCR(IIB-1:IIE,IJS:IJN,:)*(ZPHAT(IIB-1:IIE,IJS:IJN,:) - PSRC(IIB-1:IIE,IJS:IJN,:)) + &
-!!$ ZCR(IIB-1:IIE,IJS:IJN,:)*(1.0 + ZCR(IIB-1:IIE,IJS:IJN,:)) * &
-!!$ (ZPHAT(IIB-1:IIE,IJS:IJN,:) - 2.0*PSRC(IIB-1:IIE,IJS:IJN,:) + ZPHAT(IIB:IIE+1,IJS:IJN,:))
-!
+
ZFNEG(IIB:IIE,IJS:IJN,:) = ZPHAT(IIB:IIE,IJS:IJN,:) + &
ZCR(IIB:IIE,IJS:IJN,:)*(ZPHAT(IIB:IIE,IJS:IJN,:) - PSRC(IIB:IIE,IJS:IJN,:)) + &
ZCR(IIB:IIE,IJS:IJN,:)*(1.0 + ZCR(IIB:IIE,IJS:IJN,:)) * &
(ZPHAT(IIB:IIE,IJS:IJN,:) - 2.0*PSRC(IIB:IIE,IJS:IJN,:) + ZPHAT(IIB+1:IIE+1,IJS:IJN,:))
!
+!TEMPO_JUAN CALL GET_HALO(ZFNEG,HDIR="Z0_X") ! JUAN
!
- CALL GET_HALO(ZFNEG,HDIR="Z0_X") ! JUAN
-!
- IF (LEAST_ll()) THEN
+ IF (GEAST) THEN
!
! in OPEN case ZCR(IIB-1) is not used, so we also set ZFNEG(IIB-1) = 0
!
@@ -1555,13 +1664,18 @@ CALL GET_HALO(ZFPOS,HDIR="Z0_X") ! JUAN
!
! calculate the advection
!
- PR = PSRC * PRHO - &
- DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+ mxm(ZRHO_MXM,PRHO)
+ ZCR_MXM = ZCR* ZRHO_MXM*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + ZFNEG*(0.5-SIGN(0.5,ZCR)) )
+ dxf(ZCR_DXF,ZCR_MXM)
+ PR = PSRC * PRHO - ZCR_DXF
+!!$ PR = PSRC * PRHO - &
+!!$ DXF( ZCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+!
!
! in OPEN case fix boundary conditions
!
- IF (LWEST_ll()) THEN
+ IF (GWEST) THEN
WHERE ( ZCR(IIB,IJS:IJN,:) <= 0. ) ! OUTFLOW condition
PR(IIB-1,IJS:IJN,:) = 2.*PR(IIB,IJS:IJN,:) - PR(IIB+1,IJS:IJN,:)
ELSEWHERE
@@ -1569,18 +1683,22 @@ CALL GET_HALO(ZFPOS,HDIR="Z0_X") ! JUAN
END WHERE
ENDIF
!
- IF (LEAST_ll()) THEN
+ IF (GEAST) THEN
WHERE ( ZCR(IIE,IJS:IJN,:) >= 0. ) ! OUTFLOW condition
PR(IIE+1,IJS:IJN,:) = 2.*PR(IIE,IJS:IJN,:) - PR(IIE-1,IJS:IJN,:)
ELSEWHERE
PR(IIE+1,IJS:IJN,:) = PR(IIE,IJS:IJN,:)
END WHERE
ENDIF
+
+!$acc end region
+!$acc update host (pr)
+!$acc end data region
!
CALL GET_HALO(PR,HDIR="S0_X")
CALL MPPDB_CHECK3DM("PPM::PPM_S0_X OPEN ::PR",PRECISION,PR)
-!
-END SELECT
+!!$!
+!!$END SELECT
!
!-------------------------------------------------------------------------------
CALL DEL_HALO2_ll(TZ_PSRC_HALO2_ll)
@@ -1613,6 +1731,7 @@ USE MODI_GET_HALO
USE MODD_LUNIT
USE MODD_CONF
!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+USE MODD_PARAMETERS, ONLY : JPHEXT
USE MODE_MPPDB
!
IMPLICIT NONE
@@ -1649,6 +1768,14 @@ TYPE(HALO2LIST_ll), POINTER :: TZ_PHAT_HALO2_ll ! halo2 for ZPHAT
INTEGER :: ILUOUT,IRESP ! for prints
INTEGER :: IIW,IIA
!END JUAN PPM_LL
+!JUAN ACC
+REAL, DIMENSION(SIZE(ZCR,1),SIZE(ZCR,2),SIZE(ZCR,3)) :: ZRHO_MYM , ZCR_MYM , ZCR_DYF
+INTEGER :: I,J,K ,IIU,IJU,IKU
+LOGICAL :: GSOUTH , GNORTH
+! inline shuman with macro
+#define dyf(PDYF,PA) PDYF(:,1:IJU-1,:) = PA(:,2:IJU,:) - PA(:,1:IJU-1,:); PDYF(:,IJU,:) = PDYF(:,2*JPHEXT,:) ! DYF(PDYF,PA)
+#define mym(PMYM,PA) PMYM(:,2:IJU,:) = 0.5*( PA(:,2:IJU,:)+PA(:,1:IJU-1,:) ) ; PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:) ! MYM(PMYM,PA)
+!JUAN ACC
!
!-------------------------------------------------------------------------------
!
@@ -1658,8 +1785,14 @@ INTEGER :: IIW,IIA
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IIW=IIB
IIA=IIE
-!!$IIW=IIB-1
-!!$IIA=IIE+1
+
+IIU=size(psrc,1)
+IJU=size(psrc,2)
+IKU=size(psrc,3)
+
+GSOUTH=LSOUTH_ll()
+GNORTH=LNORTH_ll()
+
!
!-------------------------------------------------------------------------------
!
@@ -1672,10 +1805,15 @@ CALL GET_HALO2(PSRC,TZ_PSRC_HALO2_ll)
!
! Initialize with relalistic value all work array
!
+!$acc data region local(PR,ZPHAT,ZFPOS,ZFNEG,ZRHO_MYM,ZCR_MYM,ZCR_DYF) copyin (psrc,zcr,prho)
+!$acc region
+
ZPHAT=PSRC
ZFPOS=PSRC
ZFNEG=PSRC
PR=PSRC
+! acc end region
+! acc end data region
!
!-------------------------------------------------------------------------------
!
@@ -1685,130 +1823,96 @@ ZPHAT(IIW:IIA,IJB+1:IJE,:) = (7.0 * &
(PSRC(IIW:IIA,IJB+1:IJE,:) + PSRC(IIW:IIA,IJB:IJE-1,:)) - &
(PSRC(IIW:IIA,IJB+2:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE-2,:))) / 12.0
!
-SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
-CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
-!
-!!$ ZPHAT(:,IJB,:) = (7.0 * &
-!!$ (PSRC(:,IJB,:) + PSRC(:,IJB-1,:)) - &
-!!$ (PSRC(:,IJB+1,:) + PSRC(:,IJE-1,:))) / 12.0
+!!$SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
+!!$CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
!!$!
-!!$ ZPHAT(:,IJE+1,:) = ZPHAT(:,IJB,:)
-!!$ ZPHAT(:,IJB-1,:) = ZPHAT(:,IJE,:)
-!
-! SOUTH BOUND
-!
- ZPHAT(IIW:IIA,IJB,:) = ( 7.0 * &
- ( PSRC(IIW:IIA,IJB ,:) + PSRC(IIW:IIA,IJB-1,:) ) - &
- ( PSRC(IIW:IIA,IJB+1,:) + TZ_PSRC_HALO2_ll%HALO2%SOUTH(IIW:IIA,:) ) ) / 12.0
-! <=> SOUTH B ( PSRC(IIW:IIA,IJB+1,:) + PSRC(IIW:IIA,IJB-2,:) ) ) / 12.0
-!
-! The ZPHAT(:,IJB-1,:) doesn't matter only define an realistic value
-!
-!!$ ZPHAT(:,IJB-1,:) = ZPHAT(:,IJB,:)
-!
-! NORTH BOUND
-!
-! The ZPHAT(:IJE+1,:) doesn't matter only define an realistic value
-!
-!!$ ZPHAT(:,IJE+1,:) = ZPHAT(:,IJE,:)
-!
-! update ZPHAT HALO before next/further utilisation
-!
-CALL GET_HALO(ZPHAT,HDIR="Z0_Y")
-!
-! calculate the fluxes:
-!
-!!$ ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZPHAT(IIW:IIA,IJB:IJE+1,:) - &
-!!$ ZCR(IIW:IIA,IJB:IJE+1,:)*(ZPHAT(IIW:IIA,IJB:IJE+1,:) - PSRC(IIW:IIA,IJB-1:IJE,:)) - &
-!!$ ZCR(IIW:IIA,IJB:IJE+1,:)*(1.0 - ZCR(IIW:IIA,IJB:IJE+1,:)) * &
-!!$ (ZPHAT(IIW:IIA,IJB-1:IJE,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE,:) + ZPHAT(IIW:IIA,IJB:IJE+1,:))
-!
- ZFPOS(IIW:IIA,IJB:IJE,:) = ZPHAT(IIW:IIA,IJB:IJE,:) - &
- ZCR(IIW:IIA,IJB:IJE,:)*(ZPHAT(IIW:IIA,IJB:IJE,:) - PSRC(IIW:IIA,IJB-1:IJE-1,:)) - &
- ZCR(IIW:IIA,IJB:IJE,:)*(1.0 - ZCR(IIW:IIA,IJB:IJE,:)) * &
- (ZPHAT(IIW:IIA,IJB-1:IJE-1,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE-1,:) + ZPHAT(IIW:IIA,IJB:IJE,:))
-!
-!!$ ZFPOS(:,IJB-1,:) = ZFPOS(:,IJE,:)
-CALL GET_HALO(ZFPOS,HDIR="Z0_Y") ! JUAN
-!
-!!$ ZFNEG(IIW:IIA,IJB-1:IJE,:) = ZPHAT(IIW:IIA,IJB-1:IJE,:) + &
-!!$ ZCR(IIW:IIA,IJB-1:IJE,:)*(ZPHAT(IIW:IIA,IJB-1:IJE,:) - PSRC(IIW:IIA,IJB-1:IJE,:)) + &
-!!$ ZCR(IIW:IIA,IJB-1:IJE,:)*(1.0 + ZCR(IIW:IIA,IJB-1:IJE,:)) * &
-!!$ (ZPHAT(IIW:IIA,IJB-1:IJE,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE,:) +ZPHAT(IIW:IIA,IJB:IJE+1,:))
-!
- ZFNEG(IIW:IIA,IJB:IJE,:) = ZPHAT(IIW:IIA,IJB:IJE,:) + &
- ZCR(IIW:IIA,IJB:IJE,:)*(ZPHAT(IIW:IIA,IJB:IJE,:) - PSRC(IIW:IIA,IJB:IJE,:)) + &
- ZCR(IIW:IIA,IJB:IJE,:)*(1.0 + ZCR(IIW:IIA,IJB:IJE,:)) * &
- (ZPHAT(IIW:IIA,IJB:IJE,:) - 2.0*PSRC(IIW:IIA,IJB:IJE,:) +ZPHAT(IIW:IIA,IJB+1:IJE+1,:))
-!
-!
-! define fluxes for CYCL BC outside physical domain
-!!$ ZFNEG(:,IJE+1,:) = ZFNEG(:,IJB,:)
-CALL GET_HALO(ZFNEG,HDIR="Z0_Y") ! JUAN
-!
-! calculate the advection
-!
- PR = PSRC * PRHO - &
- DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
- CALL GET_HALO(PR,HDIR="S0_Y")
- CALL MPPDB_CHECK3DM("PPM::PPM_S0_Y CYCL ::PR",PRECISION,PR)
-!
-CASE ('OPEN')
-!
-!!$ ZPHAT(:,IJB,:) = 0.5*(PSRC(:,IJB-1,:) + PSRC(:,IJB,:))
-!!$ ZPHAT(:,IJB-1,:) = ZPHAT(:,IJB,:) ! not used
-!!$ ZPHAT(:,IJE+1,:) = 0.5*(PSRC(:,IJE,:) + PSRC(:,IJE+1,:))
+!!$!
+!!$! SOUTH BOUND
+!!$!
+!!$ ZPHAT(IIW:IIA,IJB,:) = ( 7.0 * &
+!!$ ( PSRC(IIW:IIA,IJB ,:) + PSRC(IIW:IIA,IJB-1,:) ) - &
+!!$ ( PSRC(IIW:IIA,IJB+1,:) + TZ_PSRC_HALO2_ll%HALO2%SOUTH(IIW:IIA,:) ) ) / 12.0
+!!$! <=> SOUTH B ( PSRC(IIW:IIA,IJB+1,:) + PSRC(IIW:IIA,IJB-2,:) ) ) / 12.0
+!!$!
+!!$! The ZPHAT(:,IJB-1,:) doesn't matter only define an realistic value
+!!$!
+!!$! NORTH BOUND
+!!$!
+!!$! The ZPHAT(:IJE+1,:) doesn't matter only define an realistic value
+!!$!
+!!$!
+!!$! update ZPHAT HALO before next/further utilisation
+!!$!
+!!$CALL GET_HALO(ZPHAT,HDIR="Z0_Y")
+!!$!
+!!$! calculate the fluxes:
+!!$!
+!!$ ZFPOS(IIW:IIA,IJB:IJE,:) = ZPHAT(IIW:IIA,IJB:IJE,:) - &
+!!$ ZCR(IIW:IIA,IJB:IJE,:)*(ZPHAT(IIW:IIA,IJB:IJE,:) - PSRC(IIW:IIA,IJB-1:IJE-1,:)) - &
+!!$ ZCR(IIW:IIA,IJB:IJE,:)*(1.0 - ZCR(IIW:IIA,IJB:IJE,:)) * &
+!!$ (ZPHAT(IIW:IIA,IJB-1:IJE-1,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE-1,:) + ZPHAT(IIW:IIA,IJB:IJE,:))
+!!$!
+!!$!
+!!$CALL GET_HALO(ZFPOS,HDIR="Z0_Y") ! JUAN
+!!$!
+!!$ ZFPOS(IIW:IIA,IJB:IJE,:) = ZPHAT(IIW:IIA,IJB:IJE,:) - &
+!!$ ZCR(IIW:IIA,IJB:IJE,:)*(ZPHAT(IIW:IIA,IJB:IJE,:) - PSRC(IIW:IIA,IJB-1:IJE-1,:)) - &
+!!$ ZCR(IIW:IIA,IJB:IJE,:)*(1.0 - ZCR(IIW:IIA,IJB:IJE,:)) * &
+!!$ (ZPHAT(IIW:IIA,IJB-1:IJE-1,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE-1,:) + ZPHAT(IIW:IIA,IJB:IJE,:))
+!!$!
+!!$!
+!!$! define fluxes for CYCL BC outside physical domain
+!!$CALL GET_HALO(ZFNEG,HDIR="Z0_Y") ! JUAN
+!!$!
+!!$! calculate the advection
+!!$!
+!!$ PR = PSRC * PRHO - &
+!!$ DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+!!$ CALL GET_HALO(PR,HDIR="S0_Y") ! JUAN
+!!$ CALL MPPDB_CHECK3DM("PPM::PPM_S0_Y CYCL ::PR",PRECISION,PR
+!!$!
+!!$CASE ('OPEN')
!
!
! SOUTH BOUND
!
- IF ( .NOT. LSOUTH_ll()) THEN
- ZPHAT(IIW:IIA,IJB ,:) = (7.0 * &
- (PSRC(IIW:IIA,IJB ,:) + PSRC(IIW:IIA,IJB-1,:)) - &
- (PSRC(IIW:IIA,IJB+1,:) + TZ_PSRC_HALO2_ll%HALO2%SOUTH(IIW:IIA,:) )) / 12.0
-! <=> SOUTH BOUND (PSRC(IIW:IIA,IJB+1,:) + PSRC(IIW:IIA,IJB-2,:) )) / 12.0
- ENDIF
+!TEMPO_JUAN IF ( .NOT. GSOUTH) THEN
+!TEMPO_JUAN ZPHAT(IIW:IIA,IJB ,:) = (7.0 * &
+!TEMPO_JUAN (PSRC(IIW:IIA,IJB ,:) + PSRC(IIW:IIA,IJB-1,:)) - &
+!TEMPO_JUAN (PSRC(IIW:IIA,IJB+1,:) + TZ_PSRC_HALO2_ll%HALO2%SOUTH(IIW:IIA,:) )) / 12.0
+!TEMPO_JUAN! <=> SOUTH BOUND (PSRC(IIW:IIA,IJB+1,:) + PSRC(IIW:IIA,IJB-2,:) )) / 12.0
+!TEMPO_JUAN ENDIF
!
-CALL GET_HALO(ZPHAT,HDIR="Z0_Y")
+!JUAN_TEMP CALL GET_HALO(ZPHAT,HDIR="Z0_Y")
!
- IF (LSOUTH_ll()) THEN
+ IF (GSOUTH) THEN
ZPHAT(IIW:IIA,IJB ,:) = 0.5*(PSRC(IIW:IIA,IJB-1,:) + PSRC(IIW:IIA,IJB,:))
ZPHAT(IIW:IIA,IJB-1,:) = ZPHAT(IIW:IIA,IJB,:)
ENDIF
!
! NORTH BOUND
!
- IF (LNORTH_ll()) THEN
+ IF (GNORTH) THEN
ZPHAT(IIW:IIA,IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJE,:) + PSRC(IIW:IIA,IJE+1,:))
ENDIF
!
!
! update ZPHAT HALO before next/further utilisation
!
-!!$CALL GET_HALO(ZPHAT,HDIR="Z0_Y")
!
! calculate the fluxes:
! positive fluxes
-!!$ ZFPOS(:,IJB+1:IJE+1,:) = ZPHAT(:,IJB+1:IJE+1,:) - &
-!!$ ZCR(:,IJB+1:IJE+1,:)*(ZPHAT(:,IJB+1:IJE+1,:) - PSRC(:,IJB:IJE,:)) - &
-!!$ ZCR(:,IJB+1:IJE+1,:)*(1.0 - ZCR(:,IJB+1:IJE+1,:)) * &
-!!$ (ZPHAT(:,IJB:IJE,:) - 2.0*PSRC(:,IJB:IJE,:) + ZPHAT(:,IJB+1:IJE+1,:))
-!
-!!$ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZPHAT(IIW:IIA,IJB:IJE+1,:) - &
-!!$ ZCR(IIW:IIA,IJB:IJE+1,:)*( ZPHAT(IIW:IIA,IJB:IJE+1,:) - PSRC(IIW:IIA,IJB-1:IJE ,:) ) - &
-!!$ ZCR(IIW:IIA,IJB:IJE+1,:)*( 1.0 - ZCR(IIW:IIA,IJB :IJE+1,:) ) * &
-!!$ (ZPHAT(IIW:IIA,IJB-1:IJE,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE,:) + ZPHAT(IIW:IIA,IJB:IJE+1,:))
!
ZFPOS(IIW:IIA,IJB:IJE,:) = ZPHAT(IIW:IIA,IJB:IJE,:) - &
ZCR(IIW:IIA,IJB:IJE,:)*(ZPHAT(IIW:IIA,IJB:IJE,:) - PSRC(IIW:IIA,IJB-1:IJE-1,:)) - &
ZCR(IIW:IIA,IJB:IJE,:)*(1.0 - ZCR(IIW:IIA,IJB:IJE,:)) * &
(ZPHAT(IIW:IIA,IJB-1:IJE-1,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE-1,:) + ZPHAT(IIW:IIA,IJB:IJE,:))
!
-CALL GET_HALO(ZFPOS,HDIR="Z0_Y") ! JUAN
+!JUAN_TEMP CALL GET_HALO(ZFPOS,HDIR="Z0_Y") ! JUAN
!
! positive flux on the SOUTH boundary
- IF (LSOUTH_ll()) THEN
+ IF (GSOUTH) THEN
ZFPOS(IIW:IIA,IJB,:) = (PSRC(IIW:IIA,IJB-1,:) - ZPHAT(IIW:IIA,IJB,:))*ZCR(IIW:IIA,IJB,:) + &
ZPHAT(IIW:IIA,IJB,:)
!
@@ -1817,24 +1921,15 @@ CALL GET_HALO(ZFPOS,HDIR="Z0_Y") ! JUAN
ENDIF
!
! negative fluxes
-!!$ ZFNEG(:,IJB:IJE,:) = ZPHAT(:,IJB:IJE,:) + &
-!!$ ZCR(:,IJB:IJE,:)*(ZPHAT(:,IJB:IJE,:) - PSRC(:,IJB:IJE,:)) + &
-!!$ ZCR(:,IJB:IJE,:)*(1.0 + ZCR(:,IJB:IJE,:)) * &
-!!$ (ZPHAT(:,IJB:IJE,:) - 2.0*PSRC(:,IJB:IJE,:) +ZPHAT(:,IJB+1:IJE+1,:))
-!
-!!$ ZFNEG(IIW:IIA,IJB-1:IJE,:) = ZPHAT(IIW:IIA,IJB-1:IJE,:) + &
-!!$ ZCR(IIW:IIA,IJB-1:IJE,:)*(ZPHAT(IIW:IIA,IJB-1:IJE,:) - PSRC(IIW:IIA,IJB-1:IJE,:)) + &
-!!$ ZCR(IIW:IIA,IJB-1:IJE,:)*(1.0 + ZCR(IIW:IIA,IJB-1:IJE,:)) * &
-!!$ (ZPHAT(IIW:IIA,IJB-1:IJE,:) - 2.0*PSRC(IIW:IIA,IJB-1:IJE,:) +ZPHAT(IIW:IIA,IJB:IJE+1,:))
!
ZFNEG(IIW:IIA,IJB:IJE,:) = ZPHAT(IIW:IIA,IJB:IJE,:) + &
ZCR(IIW:IIA,IJB:IJE,:)*(ZPHAT(IIW:IIA,IJB:IJE,:) - PSRC(IIW:IIA,IJB:IJE,:)) + &
ZCR(IIW:IIA,IJB:IJE,:)*(1.0 + ZCR(IIW:IIA,IJB:IJE,:)) * &
(ZPHAT(IIW:IIA,IJB:IJE,:) - 2.0*PSRC(IIW:IIA,IJB:IJE,:) +ZPHAT(IIW:IIA,IJB+1:IJE+1,:))
!
- CALL GET_HALO(ZFNEG,HDIR="Z0_Y") ! JUAN
+ !JUAN_TEMP CALL GET_HALO(ZFNEG,HDIR="Z0_Y") ! JUAN
!
- IF (LNORTH_ll()) THEN
+ IF (GNORTH) THEN
! this is not used
ZFNEG(IIW:IIA,IJB-1,:) = 0.0
!
@@ -1845,13 +1940,18 @@ CALL GET_HALO(ZFPOS,HDIR="Z0_Y") ! JUAN
!
! calculate the advection
!
- PR = PSRC * PRHO - &
- DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
- ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+ mym(ZRHO_MYM,PRHO)
+ ZCR_MYM = ZCR* ZRHO_MYM*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + ZFNEG*(0.5-SIGN(0.5,ZCR)) )
+ dyf(ZCR_DYF,ZCR_MYM)
+ PR = PSRC * PRHO - ZCR_DYF
+!!$ PR = PSRC * PRHO - &
+!!$ DYF( ZCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,ZCR)) + &
+!!$ ZFNEG*(0.5-SIGN(0.5,ZCR)) ) )
+!
!
! in OPEN case fix boundary conditions
!
- IF (LSOUTH_ll()) THEN
+ IF (GSOUTH) THEN
WHERE ( ZCR(IIW:IIA,IJB,:) <= 0. ) ! OUTFLOW condition
PR(IIW:IIA,IJB-1,:) = 1.0 * 2.*PR(IIW:IIA,IJB,:) - PR(IIW:IIA,IJB+1,:)
ELSEWHERE
@@ -1859,21 +1959,23 @@ CALL GET_HALO(ZFPOS,HDIR="Z0_Y") ! JUAN
END WHERE
ENDIF
!
- IF (LNORTH_ll()) THEN
+ IF (GNORTH) THEN
WHERE ( ZCR(IIW:IIA,IJE,:) >= 0. ) ! OUTFLOW condition
PR(IIW:IIA,IJE+1,:) = 1.0 * 2.*PR(IIW:IIA,IJE,:) - PR(IIW:IIA,IJE-1,:)
ELSEWHERE
PR(IIW:IIA,IJE+1,:) = PR(IIW:IIA,IJE,:)
END WHERE
ENDIF
+
+!$acc end region
+!$acc update host (pr)
+!$acc end data region
!
CALL GET_HALO(PR,HDIR="S0_Y")
CALL MPPDB_CHECK3DM("PPM::PPM_S0_Y OPEN ::PR",PRECISION,PR)
!
-!
-END SELECT
-!
-!!$CALL GET_HALO(PR)
+
+!!$END SELECT
!
CALL DEL_HALO2_ll(TZ_PSRC_HALO2_ll)
!
diff --git a/MNH/shuman.f90 b/MNH/shuman.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a9e8c285bf4ec46df14e29709c2c8a1106c5712a
--- /dev/null
+++ b/MNH/shuman.f90
@@ -0,0 +1,1249 @@
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/shuman.f90,v $ $Revision: 1.1.10.1 $
+! MASDEV4_7 operators 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+! ##################
+ MODULE MODI_SHUMAN
+! ##################
+!
+INTERFACE
+!
+FUNCTION DXF(PA) RESULT(PDXF)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDXF ! result at mass
+ ! localization
+END FUNCTION DXF
+!
+FUNCTION DXM(PA) RESULT(PDXM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDXM ! result at flux
+ ! side
+END FUNCTION DXM
+!
+FUNCTION DYF(PA) RESULT(PDYF)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDYF ! result at mass
+ ! localization
+END FUNCTION DYF
+!
+FUNCTION DYM(PA) RESULT(PDYM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDYM ! result at flux
+ ! side
+END FUNCTION DYM
+!
+FUNCTION DZF(PA) RESULT(PDZF)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZF ! result at mass
+ ! localization
+END FUNCTION DZF
+!
+FUNCTION DZM(PA) RESULT(PDZM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZM ! result at flux
+ ! side
+END FUNCTION DZM
+!
+FUNCTION MXF(PA) RESULT(PMXF)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMXF ! result at mass
+ ! localization
+END FUNCTION MXF
+!
+FUNCTION MXM(PA) RESULT(PMXM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMXM ! result at flux localization
+END FUNCTION MXM
+!
+FUNCTION MYF(PA) RESULT(PMYF)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMYF ! result at mass
+ ! localization
+END FUNCTION MYF
+!
+FUNCTION MYM(PA) RESULT(PMYM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMYM ! result at flux localization
+END FUNCTION MYM
+!
+FUNCTION MZF(PA) RESULT(PMZF)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZF ! result at mass
+ ! localization
+END FUNCTION MZF
+!
+FUNCTION MZM(PA) RESULT(PMZM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM ! result at flux localization
+END FUNCTION MZM
+!
+END INTERFACE
+!
+END MODULE MODI_SHUMAN
+!
+!
+! ###############################
+ FUNCTION MXF(PA) RESULT(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
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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 + JPHEXT
+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 JJK=1,IJU*IKU
+ PMXF(IIU,JJK,1) = PMXF(2*JPHEXT,JJK,1)
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MXF
+! ###############################
+ FUNCTION MXM(PA) RESULT(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
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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
+!
+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 + JPHEXT
+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 JJK=1,IJU*IKU
+ PMXM(1,JJK,1) = PMXM(IIU-2*JPHEXT+1,JJK,1)
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MXM
+! ###############################
+ FUNCTION MYF(PA) RESULT(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
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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
+!
+PMYF(:,IJU,:) = PMYF(:,2*JPHEXT,:)
+!
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MYF
+! ###############################
+ FUNCTION MYM(PA) RESULT(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
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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
+!
+PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:)
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MYM
+! ###############################
+ FUNCTION MZF(PA) RESULT(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
+! ------------
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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) = -999.
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MZF
+! ###############################
+ FUNCTION MZM(PA) RESULT(PMZM)
+! ###############################
+!
+!!**** *MZM* - Shuman operator : mean operator in z direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute a mean
+! along the z direction (K index) for a field PA localized at a mass
+! point. The result is localized at a z-flux point (w point).
+!
+!!** METHOD
+!! ------
+!! The result PMZM(:,:,k) is defined by 0.5*(PA(:,:,k)+PA(:,:,k-1))
+!! At k=1, PMZM(:,:,1) is defined by -999.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (SHUMAN operators)
+!! Technical specifications Report of The Meso-NH (chapters 3)
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/07/94
+!! optimisation 20/08/00 J. Escobar
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM ! result at flux localization
+!
+!* 0.2 Declarations of local variables
+! -------------------------------
+!
+INTEGER :: JK ! Loop index in z direction
+INTEGER :: IKU ! upper bound in z direction of PA
+!
+INTEGER :: IIU,IJU
+INTEGER :: JIJ,JI,JJ
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. DEFINITION OF MZM
+! ------------------
+!
+!!$IIU = SIZE(PA,1)
+!!$IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!!$!
+!!$JIJKOR = 1 + IIU*IJU
+!!$JIJKEND = IIU*IJU*IKU
+!!$!
+!!$!CDIR NODEP
+!!$!OCL NOVREC
+!!$DO JIJK=JIJKOR , JIJKEND
+!!$ PMZM(JIJK,1,1) = 0.5*( PA(JIJK,1,1)+PA(JIJK-IIU*IJU,1,1) )
+!!$END DO
+!!$
+!!$!
+!!$!CDIR NODEP
+!!$!OCL NOVREC
+!!$DO JIJ=1,IIU*IJU
+!!$ PMZM(JIJ,1,1) = -999.
+!!$END DO
+
+
+!$acc region
+PMZM(:,:,2:IKU) = 0.5*( PA(:,:,2:IKU)+PA(:,:,1:IKU-1) )
+PMZM(:,:,1) = -999.
+!$acc end region
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MZM
+! ###############################
+ FUNCTION DXF(PA) RESULT(PDXF)
+! ###############################
+!
+!!**** *DXF* - Shuman operator : finite difference operator in x direction
+!! for a variable at a flux side
+!!
+!! 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 x-flux
+! point (u point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PDXF(i,:,:) is defined by (PA(i+1,:,:)-PA(i,:,:))
+!! At i=size(PA,1), PDXF(i,:,:) are replaced by the values of PDXF,
+!! 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
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDXF ! 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 DXF
+! ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+JIJKOR = 1 + JPHEXT
+JIJKEND = IIU*IJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+ PDXF(JIJK-1,1,1) = PA(JIJK,1,1) - PA(JIJK-1,1,1)
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JJK=1,IJU*IKU
+ PDXF(IIU,JJK,1) = PDXF(2*JPHEXT,JJK,1)
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION DXF
+! ###############################
+ FUNCTION DXM(PA) RESULT(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
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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
+!
+!
+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 JJK=1,IJU*IKU
+ PDXM(1,JJK,1) = PDXM(IIU-2*JPHEXT+1,JJK,1)
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION DXM
+! ###############################
+ FUNCTION DYF(PA) RESULT(PDYF)
+! ###############################
+!
+!!**** *DYF* - Shuman operator : finite difference operator in y direction
+!! for a variable at a flux side
+!!
+!! 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 y-flux
+! point (v point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PDYF(:,j,:) is defined by (PA(:,j+1,:)-PA(:,j,:))
+!! At j=size(PA,2), PDYF(:,j,:) 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
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDYF ! 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 DYF
+! ------------------
+!
+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
+ PDYF(JIJK-IIU,1,1) = PA(JIJK,1,1) - PA(JIJK-IIU,1,1)
+END DO
+!
+PDYF(:,IJU,:) = PDYF(:,2*JPHEXT,:)
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION DYF
+! ###############################
+ FUNCTION DYM(PA) RESULT(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
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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
+!
+PDYM(:,1,:) = PDYM(:,IJU-2*JPHEXT+1,:)
+!
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION DYM
+! ###############################
+ FUNCTION DZF(PA) RESULT(PDZF)
+! ###############################
+!
+!!**** *DZF* - Shuman operator : finite difference operator in z direction
+!! for a variable at a flux side
+!!
+!! 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 z-flux
+! point (w point). The result is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!! The result PDZF(:,:,k) is defined by (PA(:,:,k+1)-PA(:,:,k))
+!! At k=size(PA,3), PDZF(:,:,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
+! ------------
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZF ! 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 DZF
+! ------------------
+!
+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
+ PDZF(JIJK-IIU*IJU,1,1) = PA(JIJK,1,1)-PA(JIJK-IIU*IJU,1,1)
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJ=1,IIU*IJU
+ PDZF(JIJ,1,IKU) = -999.
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION DZF
+! ###############################
+ FUNCTION DZM(PA) RESULT(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
+! ------------
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of argument and result
+! ------------------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: 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 FUNCTION DZM