From db35d8614da80a00d931339dd687e9ba3793e1c8 Mon Sep 17 00:00:00 2001
From: ESCOBAR Juan <escj@nuwa>
Date: Wed, 26 Oct 2022 18:13:44 +0200
Subject: [PATCH] Juan 26/10/2022:add orig
advec_4th_order_aux.f90/advecuvw_4th.f90 in ZSOLVER for GPU port of CEN4TH
---
src/ZSOLVER/advec_4th_order_aux.f90 | 735 ++++++++++++++++++++++++++++
src/ZSOLVER/advecuvw_4th.f90 | 382 +++++++++++++++
2 files changed, 1117 insertions(+)
create mode 100644 src/ZSOLVER/advec_4th_order_aux.f90
create mode 100644 src/ZSOLVER/advecuvw_4th.f90
diff --git a/src/ZSOLVER/advec_4th_order_aux.f90 b/src/ZSOLVER/advec_4th_order_aux.f90
new file mode 100644
index 000000000..041ac4517
--- /dev/null
+++ b/src/ZSOLVER/advec_4th_order_aux.f90
@@ -0,0 +1,735 @@
+!MNH_LIC Copyright 2005-2022 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! ###############################
+ MODULE MODI_ADVEC_4TH_ORDER_AUX
+! ###############################
+!
+INTERFACE
+!
+ SUBROUTINE ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PFIELDT, KGRID, &
+ PMEANX, PMEANY,TPHALO2 )
+!
+USE MODD_ARGSLIST_ll, ONLY : HALO2_ll
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMEANX, PMEANY ! fluxes
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFIELDT ! variable at t
+INTEGER, INTENT(IN) :: KGRID ! C grid localisation
+!
+TYPE(HALO2_ll), POINTER :: TPHALO2 ! halo2 for the field at t
+!
+END SUBROUTINE ADVEC_4TH_ORDER_ALGO
+!
+!-------------------------------------------------------------------------------
+!
+ FUNCTION MZF4(PA) RESULT(PMZF4)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
+ ! side
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZF4 ! result at mass
+ ! localization
+!
+END FUNCTION MZF4
+!
+!-------------------------------------------------------------------------------
+!
+ FUNCTION MZM4(PA) RESULT(PMZM4)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM4 ! result at flux
+ ! localization
+END FUNCTION MZM4
+!
+!-------------------------------------------------------------------------------
+!
+END INTERFACE
+!
+END MODULE MODI_ADVEC_4TH_ORDER_AUX
+!
+!-------------------------------------------------------------------------------
+!
+! ########################################################################
+ SUBROUTINE ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PFIELDT, KGRID, &
+ PMEANX, PMEANY,TPHALO2 )
+! ########################################################################
+!!
+!!**** *ADVEC_4TH_ORDER_ALGO * - routine used to compute 4th order horizontal
+!! advection fluxes of 3D prognostic variables
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this routine is to compute 2sd or 4th order horizontal
+!! advection fluxes of a prognostic variable.
+!!
+!!** METHOD
+!! ------
+!! In case of cyclic LBCs, the routine returns the scalar component of the
+!! advection fluxes by applying a 4th order horizontal averaging operator to
+!! the prognostic variable on each grid level. In the case of open LBCs, the
+!! averaging operator degenerates to a 2nd order one on the first ring
+!! inside the computationnal domain.
+!! The "halo2" (or the second layer of the halo) of the prognostic
+!! variable is passed as argument.
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! MODULE MODD_ARGSLIST
+!! HALO2LIST_ll : type for a list of "HALO2_lls"
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation ( routine ADVEC_4TH_ORDER )
+!! User Interface for the MesoNH Parallel Package
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/10/05
+! J. Escobar 21/03/2013: for HALOK comment all NHALO=1 test
+! P. Wautelet 21/11/2019: TPHALO2 dummy argument is no longer optional
+!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_ARGSLIST_ll, ONLY: HALO2_ll
+USE MODD_CONF
+!
+#ifdef MNH_OPENACC
+USE MODE_DEVICE
+#endif
+use mode_ll, only: GET_INDICE_ll, LWEST_ll, LEAST_ll, LNORTH_ll, LSOUTH_ll
+#ifdef MNH_OPENACC
+USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
+#endif
+use mode_mppdb
+#ifdef MNH_OPENACC
+use mode_msg
+#endif
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMEANX, PMEANY ! fluxes
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFIELDT ! variable at t
+INTEGER, INTENT(IN) :: KGRID ! C grid localisation
+!
+TYPE(HALO2_ll), POINTER :: TPHALO2 ! halo2 for the field at t
+!
+!* 0.2 Declarations of local variables :
+!
+INTEGER:: IW,IE,IS,IN,IT,IB,IWF,IEF,ISF,INF ! Coordinate of forth order diffusion area
+!
+INTEGER:: IIB,IJB ! Begining useful area in x,y directions
+INTEGER:: IIE,IJE ! End useful area in x,y directions
+!
+INTEGER:: ILUOUT,IRESP ! for prints
+!
+! JUAN ACC
+LOGICAL :: GWEST , GEAST
+LOGICAL :: GSOUTH , GNORTH
+#ifndef MNH_OPENACC
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZHALO2_WEST, ZHALO2_EAST
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZHALO2_SOUTH, ZHALO2_NORTH
+#else
+REAL, DIMENSION(:,:), pointer, contiguous :: ZHALO2_WEST, ZHALO2_EAST
+REAL, DIMENSION(:,:), pointer, contiguous :: ZHALO2_SOUTH, ZHALO2_NORTH
+#endif
+!
+
+!$acc data present( PMEANX, PMEANY, PFIELDT )
+
+IF (MPPDB_INITIALIZED) THEN
+ !Check all IN arrays
+ CALL MPPDB_CHECK(PFIELDT,"ADVEC_4TH_ORDER_ALGO beg:PFIELDT")
+END IF
+
+#ifndef MNH_OPENACC
+allocate( zhalo2_west ( size( pfieldt, 2 ), size( pfieldt, 3 ) ) )
+allocate( zhalo2_east ( size( pfieldt, 2 ), size( pfieldt, 3 ) ) )
+allocate( zhalo2_south( size( pfieldt, 2 ), size( pfieldt, 3 ) ) )
+allocate( zhalo2_north( size( pfieldt, 2 ), size( pfieldt, 3 ) ) )
+#else
+!Pin positions in the pools of MNH memory
+CALL MNH_MEM_POSITION_PIN( 'ADVEC_4TH_ORDER_ALGO' )
+
+CALL MNH_MEM_GET( zhalo2_west, size( pfieldt, 2 ), size( pfieldt, 3 ) )
+CALL MNH_MEM_GET( zhalo2_east, size( pfieldt, 2 ), size( pfieldt, 3 ) )
+CALL MNH_MEM_GET( zhalo2_south, size( pfieldt, 2 ), size( pfieldt, 3 ) )
+CALL MNH_MEM_GET( zhalo2_north, size( pfieldt, 2 ), size( pfieldt, 3 ) )
+
+!$acc data present ( zhalo2_west, zhalo2_east, zhalo2_south, zhalo2_north )
+#endif
+
+!-------------------------------------------------------------------------------
+!
+!* 0.3. COMPUTES THE DOMAIN DIMENSIONS
+! ------------------------------
+!
+#ifdef MNH_OPENACC
+CALL INIT_ON_HOST_AND_DEVICE(ZHALO2_WEST,-1e99,'ADVEC_4TH_ORDER_ALGO::ZHALO2_WEST')
+CALL INIT_ON_HOST_AND_DEVICE(ZHALO2_EAST,-1e99,'ADVEC_4TH_ORDER_ALGO::ZHALO2_EAST')
+CALL INIT_ON_HOST_AND_DEVICE(ZHALO2_SOUTH,-1e99,'ADVEC_4TH_ORDER_ALGO::ZHALO2_SOUTH')
+CALL INIT_ON_HOST_AND_DEVICE(ZHALO2_NORTH,-1e99,'ADVEC_4TH_ORDER_ALGO::ZHALO2_NORTH')
+#endif
+!
+CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+!
+GWEST = LWEST_ll()
+GEAST = LEAST_ll()
+GSOUTH = LSOUTH_ll()
+GNORTH = LNORTH_ll()
+!
+!-------------------------------------------------------------------------------
+!
+!* 0.4. INITIALIZE THE FIELDS
+! ---------------------
+!
+!$acc kernels present(PMEANX,PMEANY)
+PMEANX(:,:,:) = 0.0
+PMEANY(:,:,:) = 0.0
+!$acc end kernels
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. CALCULATE THE NUMERICAL MEAN IN THE X DIRECTION
+! -----------------------------------------------
+!
+SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
+!
+!* 1.1 CYCLIC CASE IN THE X DIRECTION:
+!
+CASE ('CYCL') ! In that case one must have HLBCX(1) == HLBCX(2)
+!
+#ifdef MNH_OPENACC
+call Print_msg( NVERB_WARNING, 'GEN', 'ADVEC_4TH_ORDER_ALGO', 'OpenACC: HLBCX(1) AND CYCL not yet tested' )
+#endif
+ZHALO2_WEST(:,:) = TPHALO2%WEST(:,:)
+ZHALO2_EAST(:,:) = TPHALO2%EAST(:,:)
+!$acc update device (ZHALO2_WEST,ZHALO2_EAST)
+!
+!$acc kernels present(PMEANX)
+ IW=IIB+1
+ IE=IIE
+!
+ IF(KGRID == 2) THEN
+ IWF=IW-1
+ IEF=IE-1
+ ELSE
+ IWF=IW
+ IEF=IE
+ END IF
+!
+!* lateral boundary conditions
+ PMEANX(IWF-1,:,:) = (7.0*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) ) - &
+ ( PFIELDT(IW,:,:)+ZHALO2_WEST(:,:) ) )/12.0
+!
+ PMEANX(IEF+1,:,:) = (7.0*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) ) - &
+ ( ZHALO2_EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
+!
+!* inner domain
+ PMEANX(IWF:IEF,:,:) = (7.0*( PFIELDT(IW:IE,:,:)+PFIELDT(IW-1:IE-1,:,:) ) - &
+ ( PFIELDT(IW+1:IE+1,:,:)+PFIELDT(IW-2:IE-2,:,:) ) )/12.0
+!$acc end kernels
+!
+!!$!
+!!$
+!!$ IF(NHALO == 1) THEN
+!!$ PMEANX(IWF-1,:,:) = (7.0*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) ) - &
+!!$ ( PFIELDT(IW,:,:)+ZPHALO2_WEST(:,:) ) )/12.0
+!!$!
+!!$ PMEANX(IEF+1,:,:) = (7.0*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) ) - &
+!!$ ( ZPHALO2_EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
+!!$ ENDIF
+!!$!
+!!$ PMEANX(IWF:IEF,:,:) = (7.0*( PFIELDT(IW:IE,:,:)+PFIELDT(IW-1:IE-1,:,:) ) - &
+!!$ ( PFIELDT(IW+1:IE+1,:,:)+PFIELDT(IW-2:IE-2,:,:) ) )/12.0
+!!$!
+!* 1.2 NON CYCLIC CASE IN THE X DIRECTION
+!
+CASE ('OPEN','WALL','NEST')
+!
+ZHALO2_WEST(:,:) = TPHALO2%WEST(:,:)
+ZHALO2_EAST(:,:) = TPHALO2%EAST(:,:)
+!$acc update device (ZHALO2_WEST,ZHALO2_EAST)
+!
+!$acc kernels present(PMEANX)
+ IF (GWEST) THEN
+ IF(KGRID == 2) THEN
+ IW=IIB+2 ! special case of C grid
+ ELSE
+ IW=IIB+1
+ END IF
+ ELSE
+!!$ IF(NHALO == 1) THEN
+ IW=IIB+1
+!!$ ELSE
+!!$ IW=IIB
+!!$ ENDIF
+ ENDIF
+!!$ IF (GEAST .OR. NHALO == 1) THEN
+ IF (GEAST) THEN
+! T. Maric
+! IE=IIE-1 ! original
+ IE=IIE
+ ELSE
+ IE=IIE
+ END IF
+!
+ IF(KGRID == 2) THEN
+ IWF=IW-1
+ IEF=IE-1
+ ELSE
+ IWF=IW
+ IEF=IE
+ END IF
+!
+! T. Maric. 16.1.2006.
+! write(*,*)' IW, IE, IWF, IEF = ',IW, IE, IWF, IEF
+! stop 'Stopping in advec_4th_order_aux.f90'
+!
+!* Use a second order scheme at the physical border
+!
+ IF (GWEST) THEN
+ PMEANX(IWF-1,:,:) = 0.5*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) )
+ ! T. Maric
+ ! PMEANX(1,:,:) = PMEANX(IWF-1,:,:)
+ ! extrapolate
+ !PMEANX(1,:,:) = 0.5*(3.0*PFIELDT(1,:,:) - PFIELDT(2,:,:))
+!!$ ELSE IF (NHALO == 1) THEN
+ ELSE
+ PMEANX(IWF-1,:,:) = (7.0*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) ) - &
+ ( PFIELDT(IW,:,:)+ZHALO2_WEST(:,:) ) )/12.0
+ ENDIF
+!
+ IF (GEAST) THEN
+ PMEANX(IEF+1,:,:) = 0.5*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) )
+!!$ ELSEIF (NHALO == 1) THEN
+ ELSE
+ PMEANX(IEF+1,:,:) = (7.0*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) ) - &
+ ( ZHALO2_EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
+ ENDIF
+!
+!* Use a fourth order scheme elsewhere
+!
+ PMEANX(IWF:IEF,:,:) = (7.0*( PFIELDT(IW:IE,:,:)+PFIELDT(IW-1:IE-1,:,:) ) - &
+ ( PFIELDT(IW+1:IE+1,:,:)+PFIELDT(IW-2:IE-2,:,:) ) )/12.0
+!$acc end kernels
+END SELECT
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTES THE 4TH ORDER MEAN IN THE Y DIRECTION
+! ----------------------------------------------
+!
+IF ( .NOT. L2D ) THEN
+ SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
+!
+!* 2.1 CYCLIC CASE IN THE Y DIRECTION:
+!
+ CASE ('CYCL') ! In that case one must have HLBCY(1) == HLBCY(2)
+!
+#ifdef MNH_OPENACC
+call Print_msg( NVERB_WARNING, 'GEN', 'ADVEC_4TH_ORDER_ALGO', 'OpenACC: HLBCX(2) AND CYCL not yet tested' )
+#endif
+ZHALO2_SOUTH(:,:) = TPHALO2%SOUTH(:,:)
+ZHALO2_NORTH(:,:) = TPHALO2%NORTH(:,:)
+!$acc update device (ZHALO2_SOUTH,ZHALO2_NORTH)
+!
+!$acc kernels present(PMEANY)
+!
+!
+ IS=IJB+1
+ IN=IJE
+!
+ IF(KGRID == 3) THEN
+ ISF=IS-1
+ INF=IN-1
+ ELSE
+ ISF=IS
+ INF=IN
+ END IF
+!
+!* lateral boundary conditions
+ PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS-1,:)+PFIELDT(:,IS-2,:) ) - &
+ ( PFIELDT(:,IS,:)+ZHALO2_SOUTH(:,:) ) )/12.0
+!
+ PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:) ) - &
+ ( ZHALO2_NORTH(:,:)+PFIELDT(:,IN-1,:) ) )/12.0
+!
+!* inner domain
+ PMEANY(:,ISF:INF,:) = (7.0*( PFIELDT(:,IS:IN,:)+PFIELDT(:,IS-1:IN-1,:)) - &
+ ( PFIELDT(:,IS+1:IN+1,:)+PFIELDT(:,IS-2:IN-2,:) ))/12.0
+!$acc end kernels
+!!$!
+!!$ IF(NHALO == 1) THEN
+!!$ PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS,:)+PFIELDT(:,IS-1,:) ) - &
+!!$ ( PFIELDT(:,IS+1,:)+ZPHALO2_SOUTH(:,:) ) )/12.0
+!!$!
+!!$ PMEANY(:,ISF+1,:) = (7.0*( PFIELDT(:,IS,:)+PFIELDT(:,IS-1,:) ) - &
+!!$ ( ZPHALO2_NORTH(:,:)+PFIELDT(:,IS-2,:) ) )/12.0
+!!$ ENDIF
+!!$!
+!!$ PMEANY(:,ISF:INF,:) = (7.0*( PFIELDT(:,IS:IN,:)+PFIELDT(:,IS-1:IN-1,:)) - &
+!!$ ( PFIELDT(:,IS+1:IN+1,:)+PFIELDT(:,IS-2:IN-2,:) ))/12.0
+!!$!
+!* 2.2 NON CYCLIC CASE IN THE Y DIRECTION
+!
+ CASE ('OPEN','WALL','NEST')
+!
+ZHALO2_SOUTH(:,:) = TPHALO2%SOUTH(:,:)
+ZHALO2_NORTH(:,:) = TPHALO2%NORTH(:,:)
+!$acc update device (ZHALO2_SOUTH,ZHALO2_NORTH)
+!
+!$acc kernels present(PMEANY)
+ IF (GSOUTH) THEN
+ IF(KGRID == 3) THEN
+ IS=IJB+2 ! special case of C grid
+ ELSE
+ IS=IJB+1
+ END IF
+ ELSE
+!!$ IF(NHALO == 1) THEN
+ IS=IJB+1
+!!$ ELSE
+!!$ IS=IJB
+!!$ ENDIF
+ ENDIF
+!!$ IF (GNORTH .OR. NHALO == 1) THEN
+ IF (GNORTH) THEN
+! T. Maric
+! IN=IJE-1 ! original
+ IN=IJE
+ ELSE
+ IN=IJE
+ END IF
+!
+ IF(KGRID == 3) THEN
+ ISF=IS-1
+ INF=IN-1
+ ELSE
+ ISF=IS
+ INF=IN
+ END IF
+!
+!* Use a second order scheme at the physical border
+!
+ IF (GSOUTH) THEN
+ PMEANY(:,ISF-1,:) = 0.5*( PFIELDT(:,IS-1,:)+PFIELDT(:,IS-2,:) )
+ ! T. Maric
+ ! PMEANY(:,1,:) = PMEANY(:,ISF-1,:)
+ ! extrapolate
+ !PMEANY(:,1,:) = 0.5*(3.0*PFIELDT(:,1,:) - PFIELDT(:,2,:))
+!!$ ELSEIF (NHALO == 1) THEN
+ ELSE
+!!$ PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS,:)+PFIELDT(:,IS-1,:)) - &
+!!$ ( PFIELDT(:,IS+1,:)+TPHALO2%SOUTH(:,:) ))/12.0
+ PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS-1,:)+PFIELDT(:,IS-2,:)) - &
+ ( PFIELDT(:,IS,:)+ZHALO2_SOUTH(:,:) ))/12.0
+ ENDIF
+!
+ IF (GNORTH) THEN
+ PMEANY(:,INF+1,:) = 0.5*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:) )
+!!$ ELSEIF (NHALO == 1) THEN
+ ELSE
+!!$ PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN,:)+PFIELDT(:,IN-1,:)) - &
+!!$ ( TPHALO2%NORTH(:,:)+PFIELDT(:,IN-2,:) ))/12.0
+ PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:)) - &
+ ( ZHALO2_NORTH(:,:)+PFIELDT(:,IN-1,:) ))/12.0
+ ENDIF
+!
+!* Use a fourth order scheme elsewhere
+!
+ PMEANY(:,ISF:INF,:) = (7.0*( PFIELDT(:,IS:IN,:)+PFIELDT(:,IS-1:IN-1,:)) - &
+ ( PFIELDT(:,IS+1:IN+1,:)+PFIELDT(:,IS-2:IN-2,:) ))/12.0
+!$acc end kernels
+!
+ END SELECT
+ELSE
+!$acc kernels present(PMEANY)
+ PMEANY(:,:,:) = 0.0
+!$acc end kernels
+ENDIF
+!
+IF (MPPDB_INITIALIZED) THEN
+ !Check all OUT arrays
+ CALL MPPDB_CHECK(PMEANX,"ADVEC_4TH_ORDER_ALGO end:PMEANX")
+ CALL MPPDB_CHECK(PMEANY,"ADVEC_4TH_ORDER_ALGO end:PMEANY")
+END IF
+
+!$acc end data
+
+#ifdef MNH_OPENACC
+!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
+CALL MNH_MEM_RELEASE( 'ADVEC_4TH_ORDER_ALGO' )
+#endif
+
+!$acc end data
+
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE ADVEC_4TH_ORDER_ALGO
+!
+!-------------------------------------------------------------------------------
+!
+! ################################
+ FUNCTION MZF4(PA) RESULT(PMZF4)
+! ################################
+!
+!!**** *MZF4* - 4th order Shuman operator : mean operator in z direction for a
+!! variable at a flux side
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this function is to compute a 4th order mean value
+!! 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 PMZF4(:,:,k) is defined by
+!! PMZF4(:,:,k)=0.5*(PA(:,:,k)+PA(:,:,k+1)) at k=1 and size(PA,3)-1
+!! PMZF4(:,:,k)=-999. at k=size(PA,3)
+!! PMZF4(:,:,k)=7/12*(PA(:,:,k)+PA(:,:,k+1))
+!! -1/12*(PA(:,:,k-1)+PA(:,:,k+2)) elsewhere
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (SHUMAN operators)
+!! Technical specifications Report of The Meso-NH (chapters 3)
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Lab Aerologie *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/10/05
+!!
+!-------------------------------------------------------------------------------
+!
+!* 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)) :: PMZF4 ! 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,IIJU ! upper bounds in the x and y directions of PA
+INTEGER :: JIJ,JIJK ! running loop indexes after linearisation
+INTEGER :: JIJKOR1,JIJKEND1 ! loop boundaries
+INTEGER :: JIJKOR2,JIJKEND2 ! loop boundaries
+INTEGER :: JIJKOR3,JIJKEND3 ! loop boundaries
+!
+!-------------------------------------------------------------------------------
+
+!$acc data present( PA, PMZF4 )
+!
+!* 1. DEFINITION OF MZF4
+! ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+IIJU = IIU*IJU
+!
+JIJKOR1 = 1 + IIJU
+JIJKEND1 = 2*IIJU
+!
+!$acc kernels
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR1 , JIJKEND1
+ PMZF4(JIJK-IIJU,1,1) = 0.5*( PA(JIJK-IIJU,1,1)+PA(JIJK,1,1) )
+END DO
+!
+JIJKOR2 = 1 + JIJKEND1
+JIJKEND2 = IIJU*IKU - IIJU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR2 , JIJKEND2
+ PMZF4(JIJK-IIJU,1,1) = (7.0*( PA(JIJK,1,1)+PA(JIJK-IIJU,1,1) ) - &
+ ( PA(JIJK+IIJU,1,1)+PA(JIJK-2*IIJU,1,1) ) )/12.0
+END DO
+!
+JIJKOR3 = 1 + JIJKEND2
+JIJKEND3 = IIJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR3 , JIJKEND3
+ PMZF4(JIJK-IIJU,1,1) = 0.5*( PA(JIJK-IIJU,1,1)+PA(JIJK,1,1) )
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJ=1,IIJU
+ PMZF4(JIJ,1,IKU) = -999.
+END DO
+!$acc end kernels
+
+!$acc end data
+
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MZF4
+!
+!-------------------------------------------------------------------------------
+!
+! ################################
+ FUNCTION MZM4(PA) RESULT(PMZM4)
+! ################################
+!
+!!**** *MZM4* - 4th order Shuman operator : mean operator in z direction for a
+!! mass variable
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this function is to compute a 4th order mean value
+!! 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 PMZM4(:,:,k) is defined by
+!! PMZM4(:,:,k)=0.5*(PA(:,:,k)+PA(:,:,k+1)) at k=2 and size(PA,3)
+!! PMZM4(:,:,k)=-999. at k=1
+!! PMZM4(:,:,k)=7/12*(PA(:,:,k)+PA(:,:,k+1))
+!! -1/12*(PA(:,:,k-1)+PA(:,:,k+2)) elsewhere
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS(PMEANX,PMEANY)
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (SHUMAN operators)
+!! Technical specifications Report of The Meso-NH (chapters 3)
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Lab Aerologie *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/10/05
+!!
+!-------------------------------------------------------------------------------
+!
+!* 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)) :: PMZM4 ! 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,IIJU ! upper bounds in the x and y directions of PA
+INTEGER :: JIJ,JIJK ! running loop indexes after linearisation
+INTEGER :: JIJKOR1,JIJKEND1 ! loop boundaries
+INTEGER :: JIJKOR2,JIJKEND2 ! loop boundaries
+!
+!-------------------------------------------------------------------------------
+
+!$acc data present( PA, PMZM4 )
+!
+!* 1. DEFINITION OF MZM4
+! ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+IIJU = IIU*IJU
+!
+JIJKOR1 = 1 + IIJU
+JIJKEND1 = JIJKOR1 + IIJU
+!
+!$acc kernels
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR1 , JIJKEND1
+ PMZM4(JIJK,1,1) = 0.5*( PA(JIJK,1,1)+PA(JIJK-IIJU,1,1) )
+END DO
+!
+JIJKOR2 = 1 + JIJKEND1
+JIJKEND2 = IIJU*IKU - IIJU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR2 , JIJKEND2
+ PMZM4(JIJK,1,1) = (7.0*( PA(JIJK,1,1)+PA(JIJK-IIJU,1,1) ) - &
+ ( PA(JIJK+IIJU,1,1)+PA(JIJK-2*IIJU,1,1) ) )/12.0
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJ=1,IIJU
+ PMZM4(JIJ,1,IKU) = 0.5*( PA(JIJ,1,IKU)+PA(JIJ-IIJU,1,IKU) )
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJ=1,IIJU
+ PMZM4(JIJ,1,1) = -999.
+END DO
+!$acc end kernels
+
+!$acc end data
+
+!-------------------------------------------------------------------------------
+!
+END FUNCTION MZM4
diff --git a/src/ZSOLVER/advecuvw_4th.f90 b/src/ZSOLVER/advecuvw_4th.f90
new file mode 100644
index 000000000..d523b56f5
--- /dev/null
+++ b/src/ZSOLVER/advecuvw_4th.f90
@@ -0,0 +1,382 @@
+!MNH_LIC Copyright 2005-2022 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! ###########################
+ MODULE MODI_ADVECUVW_4TH
+! ###########################
+!
+INTERFACE
+!
+ SUBROUTINE ADVECUVW_4TH ( HLBCX, HLBCY, PRUCT, PRVCT, PRWCT, &
+ PUT, PVT, PWT, PRUS, PRVS, PRWS, TPHALO2LIST )
+!
+USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRUCT ! contravariant
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVCT ! components
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRWCT ! of momentum
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! U,V,W at t
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Source terms
+!
+TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST ! list for diffusion
+!
+END SUBROUTINE ADVECUVW_4TH
+!
+END INTERFACE
+!
+END MODULE MODI_ADVECUVW_4TH
+!
+!
+! ######################################################################
+ SUBROUTINE ADVECUVW_4TH ( HLBCX, HLBCY, PRUCT, PRVCT, PRWCT, &
+ PUT, PVT, PWT, PRUS, PRVS, PRWS, TPHALO2LIST )
+! ######################################################################
+!
+!!**** *ADVECUVW_4TH * - routine to compute the 4th order centered
+!! advection tendency of momentum (U,V,W)
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this routine is to call the ADVEC_4TH_ORDER_ALGO
+!! routine for the horizontal advection and the MZM4 and MZF4 functions for
+!! the vertical advection of momentum. The code is
+!! parallelized and works for various boundary conditions.
+!!
+!!** METHOD
+!! ------
+!! For each wind component the ADVECUVW_4TH routine calls
+!! the ADVEC_4TH_ORDER_ALGO routine which computes the numerical advection
+!! of any 3D field.
+!! The following variables are passed as argument to ADVEC_4TH_ORDER_ALGO :
+!!
+!! -- The variable at t
+!! -- The second layer of the halo of the field at t
+!! -- The horizontal advection fluxes
+!! -- The localisation on the model grid :
+!!
+!! IGRID = 1 for mass grid point
+!! IGRID = 2 for U grid point
+!! IGRID = 3 for V grid point
+!! IGRID = 4 for W grid point
+!!
+!! EXTERNAL
+!! --------
+!! BUDGET : Stores the different budget components
+!! (not used in current version)
+!! ADVEC_4TH_ORDER_ALGO : computes the horizontal advection fluxes
+!! MZF4 and MZM4 : computes the vertical advection fluxes
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODULE MODD_BUDGET:
+!! NBUMOD : model in which budget is calculated
+!! CBUTYPE : type of desired budget
+!! 'CART' for cartesian box configuration
+!! 'MASK' for budget zone defined by a mask
+!! 'NONE' ' for no budget
+!! NBUPROCCTR : process counter used for each budget variable
+!! Switches for budgets activations:
+!!
+!! MODULE MODD_ARGSLIST
+!! HALO2LIST_ll : type for a list of "HALO2_lls"
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation ( routine ADVECUVW_4TH )
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/10/05
+!! Modif
+!! J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+USE MODD_CONF
+USE MODD_GRID_n
+USE MODD_PARAMETERS
+
+USE MODE_ll
+#ifdef MNH_OPENACC
+USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
+#endif
+use mode_mppdb
+
+USE MODI_ADVEC_4TH_ORDER_AUX
+#ifndef MNH_OPENACC
+USE MODI_SHUMAN
+#else
+USE MODI_SHUMAN_DEVICE
+#endif
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRUCT ! contravariant
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVCT ! components
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRWCT ! of momentum
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! Variables at t
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Source terms
+!
+TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST ! list for diffusion
+!
+!* 0.2 Declarations of local variables :
+!
+TYPE(HALO2LIST_ll), POINTER :: TZHALO2LIST
+!
+INTEGER :: IGRID ! localisation on the model grid
+#ifndef MNH_OPENACC
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZMEANX, ZMEANY ! fluxes
+#else
+INTEGER :: IIU, IJU, IKU
+REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZMEANX, ZMEANY ! fluxes
+!
+REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZTEMP1, ZTEMP2, ZTEMP3, ZTEMP4
+#endif
+!
+#if 0
+#define dxm(PDXM,PA) PDXM(2:IIU,:,:) = PA(2:IIU,:,:) - PA(1:IIU-1,:,:) ; PDXM(1,:,:) = PDXM(IIU-2*JPHEXT+1,:,:) ! DXM(PDXM,PA)
+#define mxf(PMXF,PA) PMXF(1:IIU-1,:,:) = 0.5*( PA(2:IIU,:,:)+PA(1:IIU-1,:,:) ) ; PMXF(IIU,:,:) = PMXF(2*JPHEXT,:,:) ! MXF(PMXF,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)
+#define dyf(PDYF,PA) PDYF(:,1:IJU-1,:) = PA(:,2:IJU,:) - PA(:,1:IJU-1,:) ; PDYF(:,IJU,:) = PDYF(:,2*JPHEXT,:) ! DYF(PDYF,PA)
+#define dzf(PDZF,PA) PDZF(:,:,1:IKU-1) = PA(:,:,2:IKU) - PA(:,:,1:IKU-1) ; PDZF(:,:,IKU) = -999. ! DZF(PDZF,PA)
+#define mzm4(PMZM4,PA) PMZM4(:,:,3:IKU-1) = (7.0*( PA(:,:,3:IKU-1)+PA(:,:,2:IKU-2) ) - (PA(:,:,4:IKU)+PA(:,:,1:IKU-3) ) )/12.0 ; \
+ PMZM4(:,:,2) = 0.5*( PA(:,:,2)+PA(:,:,1) ) ; PMZM4(:,:,IKU) = 0.5*( PA(:,:,IKU)+PA(:,:,IKU-1) ) ; PMZM4(:,:,1) = -999.
+#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)
+#define dxf(PDXF,PA) PDXF(1:IIU-1,:,:) = PA(2:IIU,:,:) - PA(1:IIU-1,:,:) ; PDXF(IIU,:,:) = PDXF(2*JPHEXT,:,:) ! DXF(PDXF,PA)
+#define myf(PMYF,PA) PMYF(:,1:IJU-1,:) = 0.5*( PA(:,1:IJU-1,:)+PA(:,2:IJU,:) ) ; PMYF(:,IJU,:) = PMYF(:,2*JPHEXT,:) ! MYF(PMYF,PA)
+#define dym(PDYM,PA) PDYM(:,2:IJU,:) = PA(:,2:IJU,:) - PA(:,1:IJU-1,:) ; PDYM(:,1,:) = PDYM(:,IJU-2*JPHEXT+1,:) ! DYM(PDYM,PA)
+#define mzm(PMZM,PA) PMZM(:,:,2:IKU) = 0.5*( PA(:,:,2:IKU)+PA(:,:,1:IKU-1) ) ; PMZM(:,:,1) = -999. ! MZM(PMZM,PA)
+#define mzf(PMZF,PA) PMZF(:,:,1:IKU-1) = 0.5*( PA(:,:,1:IKU-1)+PA(:,:,2:IKU) ) ; PMZF(:,:,IKU) = -999. ! MZF(PMZF,PA)
+#define dzm(PDZM,PA) PDZM(:,:,2:IKU) = PA(:,:,2:IKU) - PA(:,:,1:IKU-1) ; PDZM(:,:,1) = -999. ! DZM(PDZM,PA)
+#define mzf4(PMZF4,PA) PMZF4(:,:,2:IKU-2) = (7.0*( PA(:,:,3:IKU-1)+PA(:,:,2:IKU-2) ) - (PA(:,:,4:IKU)+PA(:,:,1:IKU-3) ) )/12.0 ; \
+ PMZF4(:,:,1) = 0.5*( PA(:,:,2)+PA(:,:,1) ) ; PMZF4(:,:,IKU-1) = 0.5*( PA(:,:,IKU)+PA(:,:,IKU-1) ) ; PMZF4(:,:,IKU) = -999.
+#endif
+!
+IF (MPPDB_INITIALIZED) THEN
+ !Check all IN arrays
+ CALL MPPDB_CHECK(PRUCT,"ADVECUVW_4TH beg:PRUCT")
+ CALL MPPDB_CHECK(PRVCT,"ADVECUVW_4TH beg:PRVCT")
+ CALL MPPDB_CHECK(PRWCT,"ADVECUVW_4TH beg:PRWCT")
+ CALL MPPDB_CHECK(PUT,"ADVECUVW_4TH beg:PUT")
+ CALL MPPDB_CHECK(PVT,"ADVECUVW_4TH beg:PVT")
+ CALL MPPDB_CHECK(PWT,"ADVECUVW_4TH beg:PWT")
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PRUS,"ADVECUVW_4TH beg:PRUS")
+ CALL MPPDB_CHECK(PRVS,"ADVECUVW_4TH beg:PRVS")
+ CALL MPPDB_CHECK(PRWS,"ADVECUVW_4TH beg:PRWS")
+END IF
+
+#ifdef MNH_OPENACC
+IIU = SIZE( PUT, 1 )
+IJU = SIZE( PUT, 2 )
+IKU = SIZE( PUT, 3 )
+
+!Pin positions in the pools of MNH memory
+CALL MNH_MEM_POSITION_PIN()
+
+CALL MNH_MEM_GET( ZMEANX, IIU, IJU, IKU )
+CALL MNH_MEM_GET( ZMEANY, IIU, IJU, IKU )
+
+CALL MNH_MEM_GET( ZTEMP1, IIU, IJU, IKU )
+CALL MNH_MEM_GET( ZTEMP2, IIU, IJU, IKU )
+CALL MNH_MEM_GET( ZTEMP3, IIU, IJU, IKU )
+CALL MNH_MEM_GET( ZTEMP4, IIU, IJU, IKU )
+#endif
+
+!$acc data present( PRUCT, PRVCT, PRWCT, PUT, PVT, PWT, PRUS, PRVS, PRWS, ZMEANX, ZMEANY, ZTEMP1, ZTEMP2, ZTEMP3, ZTEMP4 )
+
+!-------------------------------------------------------------------------------
+!
+!* 2. CALL THE ADVEC_4TH_ORDER_ALGO ROUTINE FOR MOMENTUM
+! --------------------------------------------------
+!
+IGRID = 2
+!!$IF(NHALO == 1) THEN
+ TZHALO2LIST => TPHALO2LIST
+ CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PUT, IGRID, ZMEANX, ZMEANY, &
+ TZHALO2LIST%HALO2 )
+!!$ELSE
+!!$ CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PUT, IGRID, ZMEANX, ZMEANY)
+!!$ENDIF
+!
+#ifndef MNH_OPENACC
+PRUS(:,:,:) = PRUS(:,:,:) &
+ -DXM( MXF(PRUCT(:,:,:))*ZMEANX(:,:,:) )
+!
+PRUS(:,:,:) = PRUS(:,:,:) &
+ -DYF( MXM(PRVCT(:,:,:))*ZMEANY(:,:,:) )
+!
+PRUS(:,:,:) = PRUS(:,:,:) &
+ -DZF( MXM(PRWCT(:,:,:))*MZM4(PUT(:,:,:)) )
+#else
+call mxf_device(PRUCT,ZTEMP1)
+!$acc kernels
+ZTEMP2 = ZTEMP1 * ZMEANX
+!$acc end kernels
+call dxm_device(ZTEMP2,ZTEMP3)
+!$acc kernels
+PRUS(:,:,:) = PRUS(:,:,:) - ZTEMP3
+!$acc end kernels
+!
+call mxm_device(PRVCT,ZTEMP1)
+!$acc kernels
+ZTEMP2 = ZTEMP1 * ZMEANY
+!$acc end kernels
+call dyf_device(ZTEMP2,ZTEMP3)
+!$acc kernels
+PRUS(:,:,:) = PRUS(:,:,:) - ZTEMP3
+!$acc end kernels
+!
+ZTEMP1 = MZM4( PUT )
+call mxm_device(PRWCT,ZTEMP2)
+!$acc kernels
+ZTEMP3 = ZTEMP1 * ZTEMP2
+!$acc end kernels
+call dzf_device( ZTEMP3, ZTEMP4 )
+!$acc kernels
+PRUS(:,:,:) = PRUS(:,:,:) - ZTEMP4
+!$acc end kernels
+#endif
+!
+!
+IGRID = 3
+!!$IF(NHALO == 1) THEN
+ TZHALO2LIST => TZHALO2LIST%NEXT
+ CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PVT, IGRID, ZMEANX, ZMEANY, &
+ TZHALO2LIST%HALO2 )
+!!$ELSE
+!!$ CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PVT, IGRID, ZMEANX, ZMEANY)
+!!$ENDIF
+!
+#ifndef MNH_OPENACC
+PRVS(:,:,:) = PRVS(:,:,:) &
+ -DXF( MYM(PRUCT(:,:,:))*ZMEANX(:,:,:) )
+!
+PRVS(:,:,:) = PRVS(:,:,:) &
+ -DYM( MYF(PRVCT(:,:,:))*ZMEANY(:,:,:) )
+!
+PRVS(:,:,:) = PRVS(:,:,:) &
+ -DZF( MYM(PRWCT(:,:,:))*MZM4(PVT(:,:,:)) )
+#else
+call mym_device(PRUCT,ZTEMP1)
+!$acc kernels
+ZTEMP2 = ZTEMP1 * ZMEANX
+!$acc end kernels
+call dxf_device(ZTEMP2,ZTEMP3)
+!$acc kernels
+PRVS(:,:,:) = PRVS(:,:,:) - ZTEMP3
+!$acc end kernels
+!
+call myf_device(PRVCT,ZTEMP1)
+!$acc kernels
+ZTEMP2 = ZTEMP1 * ZMEANY
+!$acc end kernels
+call dym_device(ZTEMP2,ZTEMP3)
+!$acc kernels
+PRVS(:,:,:) = PRVS(:,:,:) - ZTEMP3
+!$acc end kernels
+!
+call mym_device(PRWCT,ZTEMP1)
+ZTEMP2 = MZM4( PVT )
+!$acc kernels
+ZTEMP3 = ZTEMP1 * ZTEMP2
+!$acc end kernels
+call dzf_device( ZTEMP3, ZTEMP4 )
+!$acc kernels
+PRVS(:,:,:) = PRVS(:,:,:) - ZTEMP4
+!$acc end kernels
+#endif
+CALL MPPDB_CHECK(PRUCT,"ADVECUVW_4TH 02: PRUCT")
+!
+!
+IGRID = 4
+!
+!!$IF(NHALO == 1) THEN
+ TZHALO2LIST => TZHALO2LIST%NEXT
+ CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PWT, IGRID, ZMEANX, ZMEANY, &
+ TZHALO2LIST%HALO2 )
+!!$ELSE
+!!$ CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PWT, IGRID, ZMEANX, ZMEANY)
+!!$ENDIF
+!
+#ifndef MNH_OPENACC
+PRWS(:,:,:) = PRWS(:,:,:) &
+ -DXF( MZM(PRUCT(:,:,:))*ZMEANX(:,:,:) )
+!
+PRWS(:,:,:) = PRWS(:,:,:) &
+ -DYF( MZM(PRVCT(:,:,:))*ZMEANY(:,:,:) )
+!
+PRWS(:,:,:) = PRWS(:,:,:) &
+ -DZM( MZF(PRWCT(:,:,:))*MZF4(PWT(:,:,:)) )
+#else
+call mzm_device(PRUCT,ZTEMP1)
+!$acc kernels
+ZTEMP2 = ZTEMP1 * ZMEANX
+!$acc end kernels
+call dxf_device(ZTEMP2,ZTEMP3)
+!$acc kernels
+PRWS(:,:,:) = PRWS(:,:,:) - ZTEMP3
+!$acc end kernels
+!
+call mzm_device(PRVCT,ZTEMP1)
+!$acc kernels
+ZTEMP2 = ZTEMP1 * ZMEANY
+!$acc end kernels
+call dyf_device(ZTEMP2,ZTEMP3)
+!$acc kernels
+PRWS(:,:,:) = PRWS(:,:,:) - ZTEMP3
+!$acc end kernels
+!
+call mzf_device( PRWCT, ZTEMP1 )
+ZTEMP2 = MZF4( PWT )
+!$acc kernels
+ZTEMP1 = ZTEMP1 * ZTEMP2
+!$acc end kernels
+call dzm_device( ZTEMP1, ZTEMP4 )
+!$acc kernels
+PRWS(:,:,:) = PRWS(:,:,:) - ZTEMP4
+!$acc end kernels
+#endif
+
+!$acc end data
+
+#ifdef MNH_OPENACC
+!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
+CALL MNH_MEM_RELEASE()
+#endif
+
+IF (MPPDB_INITIALIZED) THEN
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PRUS,"ADVECUVW_4TH end:PRUS")
+ CALL MPPDB_CHECK(PRVS,"ADVECUVW_4TH end:PRVS")
+ CALL MPPDB_CHECK(PRWS,"ADVECUVW_4TH end:PRWS")
+END IF
+
+!-------------------------------------------------------------------------------
+!
+!
+END SUBROUTINE ADVECUVW_4TH
--
GitLab