scale_element_operation_tensorprod3D.F90

Go to the documentation of this file.

 
!-------------------------------------------------------------------------------
! Warning: This file was generated from element/scale_element_operation_tensorprod3D.F90.erb.
!          Do not edit this file.
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!> module FElib / Element / Operation with 3D tensor product elements for GPU
!!
!! @par Description
!!           A module for providing mathematical operations assuming a 3D tensor product element with (p+1)^3 DOF for GPU
!!
!! @author Yuta Kawai, Xuanzhengbo Ren, and Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_element_operation_tensorprod3d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
 
  use scale_element_base, only: &
    elementbase3d, &
    elementbase3d_init, elementbase3d_final
 
  use scale_element_operation_base, only: elementoperationbase3d
  use scale_element_modalfilter, only: modalfilter  
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
 
  !> Derived type for elementwise operations with 3D tensor product elements
  type, public, abstract, extends(elementoperationbase3d) :: elementoperationtensorprod3d
  contains
    procedure(elementoperationtensorprod3d_init), deferred, public :: init
    procedure(elementoperationtensorprod3d_final), deferred, public :: final
  end type elementoperationtensorprod3d
 
  interface
    subroutine elementoperationtensorprod3d_init( this, elem3D )
      import elementoperationtensorprod3d
      import elementbase3d
      class(elementoperationtensorprod3d), intent(inout) :: this
      class(elementbase3d), intent(in), target :: elem3D
    end subroutine elementoperationtensorprod3d_init
 
    subroutine elementoperationtensorprod3d_final( this )
      import elementoperationtensorprod3d
      class(elementoperationtensorprod3d), intent(inout) :: this
    end subroutine elementoperationtensorprod3d_final
  end interface
 
  public :: elementoperationtensorprod3d_create
 
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=1
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p1
    real(rp) :: d1d(2,2)
    real(rp) :: d1d_tr(2,2)
    real(rp) :: lift_mat(2,2,2,6)
    real(rp) :: lift_mat_tr(2,2,2,6)
    real(rp) :: intrpmat_vpordm1_tr(2,2)
 
    real(rp) :: mfilter_h1d(2,2)
    real(rp) :: mfilter_h1d_tr(2,2)
    real(rp) :: mfilter_v1d_tr(2,2)
 
    real(rp) :: mfilter_tracer_h1d(2,2)
    real(rp) :: mfilter_tracer_h1d_tr(2,2)
    real(rp) :: mfilter_tracer_v1d_tr(2,2)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p1
    procedure, public :: final => element_operation_tensorprod3d_final_p1
    procedure, public :: dx => element_operation_tensorprod3d_dx_p1
    procedure, public :: dy => element_operation_tensorprod3d_dy_p1
    procedure, public :: dz => element_operation_tensorprod3d_dz_p1
    procedure, public :: lift => element_operation_tensorprod3d_lift_p1
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p1
    procedure, public :: div => element_operation_tensorprod3d_div_p1
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p1
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p1
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p1
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p1
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p1
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p1
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p1
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p1
  end type elementoperationtensorprod3d_p1
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=2
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p2
    real(rp) :: d1d(3,3)
    real(rp) :: d1d_tr(3,3)
    real(rp) :: lift_mat(3,3,3,6)
    real(rp) :: lift_mat_tr(3,3,3,6)
    real(rp) :: intrpmat_vpordm1_tr(3,3)
 
    real(rp) :: mfilter_h1d(3,3)
    real(rp) :: mfilter_h1d_tr(3,3)
    real(rp) :: mfilter_v1d_tr(3,3)
 
    real(rp) :: mfilter_tracer_h1d(3,3)
    real(rp) :: mfilter_tracer_h1d_tr(3,3)
    real(rp) :: mfilter_tracer_v1d_tr(3,3)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p2
    procedure, public :: final => element_operation_tensorprod3d_final_p2
    procedure, public :: dx => element_operation_tensorprod3d_dx_p2
    procedure, public :: dy => element_operation_tensorprod3d_dy_p2
    procedure, public :: dz => element_operation_tensorprod3d_dz_p2
    procedure, public :: lift => element_operation_tensorprod3d_lift_p2
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p2
    procedure, public :: div => element_operation_tensorprod3d_div_p2
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p2
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p2
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p2
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p2
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p2
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p2
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p2
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p2
  end type elementoperationtensorprod3d_p2
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=3
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p3
    real(rp) :: d1d(4,4)
    real(rp) :: d1d_tr(4,4)
    real(rp) :: lift_mat(4,4,4,6)
    real(rp) :: lift_mat_tr(4,4,4,6)
    real(rp) :: intrpmat_vpordm1_tr(4,4)
 
    real(rp) :: mfilter_h1d(4,4)
    real(rp) :: mfilter_h1d_tr(4,4)
    real(rp) :: mfilter_v1d_tr(4,4)
 
    real(rp) :: mfilter_tracer_h1d(4,4)
    real(rp) :: mfilter_tracer_h1d_tr(4,4)
    real(rp) :: mfilter_tracer_v1d_tr(4,4)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p3
    procedure, public :: final => element_operation_tensorprod3d_final_p3
    procedure, public :: dx => element_operation_tensorprod3d_dx_p3
    procedure, public :: dy => element_operation_tensorprod3d_dy_p3
    procedure, public :: dz => element_operation_tensorprod3d_dz_p3
    procedure, public :: lift => element_operation_tensorprod3d_lift_p3
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p3
    procedure, public :: div => element_operation_tensorprod3d_div_p3
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p3
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p3
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p3
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p3
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p3
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p3
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p3
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p3
  end type elementoperationtensorprod3d_p3
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=4
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p4
    real(rp) :: d1d(5,5)
    real(rp) :: d1d_tr(5,5)
    real(rp) :: lift_mat(5,5,5,6)
    real(rp) :: lift_mat_tr(5,5,5,6)
    real(rp) :: intrpmat_vpordm1_tr(5,5)
 
    real(rp) :: mfilter_h1d(5,5)
    real(rp) :: mfilter_h1d_tr(5,5)
    real(rp) :: mfilter_v1d_tr(5,5)
 
    real(rp) :: mfilter_tracer_h1d(5,5)
    real(rp) :: mfilter_tracer_h1d_tr(5,5)
    real(rp) :: mfilter_tracer_v1d_tr(5,5)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p4
    procedure, public :: final => element_operation_tensorprod3d_final_p4
    procedure, public :: dx => element_operation_tensorprod3d_dx_p4
    procedure, public :: dy => element_operation_tensorprod3d_dy_p4
    procedure, public :: dz => element_operation_tensorprod3d_dz_p4
    procedure, public :: lift => element_operation_tensorprod3d_lift_p4
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p4
    procedure, public :: div => element_operation_tensorprod3d_div_p4
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p4
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p4
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p4
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p4
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p4
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p4
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p4
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p4
  end type elementoperationtensorprod3d_p4
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=5
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p5
    real(rp) :: d1d(6,6)
    real(rp) :: d1d_tr(6,6)
    real(rp) :: lift_mat(6,6,6,6)
    real(rp) :: lift_mat_tr(6,6,6,6)
    real(rp) :: intrpmat_vpordm1_tr(6,6)
 
    real(rp) :: mfilter_h1d(6,6)
    real(rp) :: mfilter_h1d_tr(6,6)
    real(rp) :: mfilter_v1d_tr(6,6)
 
    real(rp) :: mfilter_tracer_h1d(6,6)
    real(rp) :: mfilter_tracer_h1d_tr(6,6)
    real(rp) :: mfilter_tracer_v1d_tr(6,6)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p5
    procedure, public :: final => element_operation_tensorprod3d_final_p5
    procedure, public :: dx => element_operation_tensorprod3d_dx_p5
    procedure, public :: dy => element_operation_tensorprod3d_dy_p5
    procedure, public :: dz => element_operation_tensorprod3d_dz_p5
    procedure, public :: lift => element_operation_tensorprod3d_lift_p5
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p5
    procedure, public :: div => element_operation_tensorprod3d_div_p5
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p5
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p5
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p5
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p5
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p5
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p5
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p5
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p5
  end type elementoperationtensorprod3d_p5
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=6
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p6
    real(rp) :: d1d(7,7)
    real(rp) :: d1d_tr(7,7)
    real(rp) :: lift_mat(7,7,7,6)
    real(rp) :: lift_mat_tr(7,7,7,6)
    real(rp) :: intrpmat_vpordm1_tr(7,7)
 
    real(rp) :: mfilter_h1d(7,7)
    real(rp) :: mfilter_h1d_tr(7,7)
    real(rp) :: mfilter_v1d_tr(7,7)
 
    real(rp) :: mfilter_tracer_h1d(7,7)
    real(rp) :: mfilter_tracer_h1d_tr(7,7)
    real(rp) :: mfilter_tracer_v1d_tr(7,7)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p6
    procedure, public :: final => element_operation_tensorprod3d_final_p6
    procedure, public :: dx => element_operation_tensorprod3d_dx_p6
    procedure, public :: dy => element_operation_tensorprod3d_dy_p6
    procedure, public :: dz => element_operation_tensorprod3d_dz_p6
    procedure, public :: lift => element_operation_tensorprod3d_lift_p6
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p6
    procedure, public :: div => element_operation_tensorprod3d_div_p6
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p6
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p6
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p6
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p6
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p6
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p6
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p6
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p6
  end type elementoperationtensorprod3d_p6
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=7
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p7
    real(rp) :: d1d(8,8)
    real(rp) :: d1d_tr(8,8)
    real(rp) :: lift_mat(8,8,8,6)
    real(rp) :: lift_mat_tr(8,8,8,6)
    real(rp) :: intrpmat_vpordm1_tr(8,8)
 
    real(rp) :: mfilter_h1d(8,8)
    real(rp) :: mfilter_h1d_tr(8,8)
    real(rp) :: mfilter_v1d_tr(8,8)
 
    real(rp) :: mfilter_tracer_h1d(8,8)
    real(rp) :: mfilter_tracer_h1d_tr(8,8)
    real(rp) :: mfilter_tracer_v1d_tr(8,8)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p7
    procedure, public :: final => element_operation_tensorprod3d_final_p7
    procedure, public :: dx => element_operation_tensorprod3d_dx_p7
    procedure, public :: dy => element_operation_tensorprod3d_dy_p7
    procedure, public :: dz => element_operation_tensorprod3d_dz_p7
    procedure, public :: lift => element_operation_tensorprod3d_lift_p7
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p7
    procedure, public :: div => element_operation_tensorprod3d_div_p7
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p7
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p7
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p7
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p7
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p7
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p7
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p7
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p7
  end type elementoperationtensorprod3d_p7
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=8
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p8
    real(rp) :: d1d(9,9)
    real(rp) :: d1d_tr(9,9)
    real(rp) :: lift_mat(9,9,9,6)
    real(rp) :: lift_mat_tr(9,9,9,6)
    real(rp) :: intrpmat_vpordm1_tr(9,9)
 
    real(rp) :: mfilter_h1d(9,9)
    real(rp) :: mfilter_h1d_tr(9,9)
    real(rp) :: mfilter_v1d_tr(9,9)
 
    real(rp) :: mfilter_tracer_h1d(9,9)
    real(rp) :: mfilter_tracer_h1d_tr(9,9)
    real(rp) :: mfilter_tracer_v1d_tr(9,9)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p8
    procedure, public :: final => element_operation_tensorprod3d_final_p8
    procedure, public :: dx => element_operation_tensorprod3d_dx_p8
    procedure, public :: dy => element_operation_tensorprod3d_dy_p8
    procedure, public :: dz => element_operation_tensorprod3d_dz_p8
    procedure, public :: lift => element_operation_tensorprod3d_lift_p8
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p8
    procedure, public :: div => element_operation_tensorprod3d_div_p8
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p8
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p8
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p8
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p8
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p8
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p8
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p8
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p8
  end type elementoperationtensorprod3d_p8
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=9
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p9
    real(rp) :: d1d(10,10)
    real(rp) :: d1d_tr(10,10)
    real(rp) :: lift_mat(10,10,10,6)
    real(rp) :: lift_mat_tr(10,10,10,6)
    real(rp) :: intrpmat_vpordm1_tr(10,10)
 
    real(rp) :: mfilter_h1d(10,10)
    real(rp) :: mfilter_h1d_tr(10,10)
    real(rp) :: mfilter_v1d_tr(10,10)
 
    real(rp) :: mfilter_tracer_h1d(10,10)
    real(rp) :: mfilter_tracer_h1d_tr(10,10)
    real(rp) :: mfilter_tracer_v1d_tr(10,10)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p9
    procedure, public :: final => element_operation_tensorprod3d_final_p9
    procedure, public :: dx => element_operation_tensorprod3d_dx_p9
    procedure, public :: dy => element_operation_tensorprod3d_dy_p9
    procedure, public :: dz => element_operation_tensorprod3d_dz_p9
    procedure, public :: lift => element_operation_tensorprod3d_lift_p9
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p9
    procedure, public :: div => element_operation_tensorprod3d_div_p9
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p9
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p9
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p9
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p9
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p9
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p9
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p9
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p9
  end type elementoperationtensorprod3d_p9
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=10
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p10
    real(rp) :: d1d(11,11)
    real(rp) :: d1d_tr(11,11)
    real(rp) :: lift_mat(11,11,11,6)
    real(rp) :: lift_mat_tr(11,11,11,6)
    real(rp) :: intrpmat_vpordm1_tr(11,11)
 
    real(rp) :: mfilter_h1d(11,11)
    real(rp) :: mfilter_h1d_tr(11,11)
    real(rp) :: mfilter_v1d_tr(11,11)
 
    real(rp) :: mfilter_tracer_h1d(11,11)
    real(rp) :: mfilter_tracer_h1d_tr(11,11)
    real(rp) :: mfilter_tracer_v1d_tr(11,11)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p10
    procedure, public :: final => element_operation_tensorprod3d_final_p10
    procedure, public :: dx => element_operation_tensorprod3d_dx_p10
    procedure, public :: dy => element_operation_tensorprod3d_dy_p10
    procedure, public :: dz => element_operation_tensorprod3d_dz_p10
    procedure, public :: lift => element_operation_tensorprod3d_lift_p10
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p10
    procedure, public :: div => element_operation_tensorprod3d_div_p10
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p10
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p10
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p10
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p10
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p10
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p10
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p10
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p10
  end type elementoperationtensorprod3d_p10
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=11
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p11
    real(rp) :: d1d(12,12)
    real(rp) :: d1d_tr(12,12)
    real(rp) :: lift_mat(12,12,12,6)
    real(rp) :: lift_mat_tr(12,12,12,6)
    real(rp) :: intrpmat_vpordm1_tr(12,12)
 
    real(rp) :: mfilter_h1d(12,12)
    real(rp) :: mfilter_h1d_tr(12,12)
    real(rp) :: mfilter_v1d_tr(12,12)
 
    real(rp) :: mfilter_tracer_h1d(12,12)
    real(rp) :: mfilter_tracer_h1d_tr(12,12)
    real(rp) :: mfilter_tracer_v1d_tr(12,12)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p11
    procedure, public :: final => element_operation_tensorprod3d_final_p11
    procedure, public :: dx => element_operation_tensorprod3d_dx_p11
    procedure, public :: dy => element_operation_tensorprod3d_dy_p11
    procedure, public :: dz => element_operation_tensorprod3d_dz_p11
    procedure, public :: lift => element_operation_tensorprod3d_lift_p11
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p11
    procedure, public :: div => element_operation_tensorprod3d_div_p11
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p11
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p11
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p11
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p11
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p11
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p11
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p11
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p11
  end type elementoperationtensorprod3d_p11
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=12
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p12
    real(rp) :: d1d(13,13)
    real(rp) :: d1d_tr(13,13)
    real(rp) :: lift_mat(13,13,13,6)
    real(rp) :: lift_mat_tr(13,13,13,6)
    real(rp) :: intrpmat_vpordm1_tr(13,13)
 
    real(rp) :: mfilter_h1d(13,13)
    real(rp) :: mfilter_h1d_tr(13,13)
    real(rp) :: mfilter_v1d_tr(13,13)
 
    real(rp) :: mfilter_tracer_h1d(13,13)
    real(rp) :: mfilter_tracer_h1d_tr(13,13)
    real(rp) :: mfilter_tracer_v1d_tr(13,13)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p12
    procedure, public :: final => element_operation_tensorprod3d_final_p12
    procedure, public :: dx => element_operation_tensorprod3d_dx_p12
    procedure, public :: dy => element_operation_tensorprod3d_dy_p12
    procedure, public :: dz => element_operation_tensorprod3d_dz_p12
    procedure, public :: lift => element_operation_tensorprod3d_lift_p12
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p12
    procedure, public :: div => element_operation_tensorprod3d_div_p12
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p12
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p12
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p12
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p12
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p12
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p12
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p12
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p12
  end type elementoperationtensorprod3d_p12
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=13
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p13
    real(rp) :: d1d(14,14)
    real(rp) :: d1d_tr(14,14)
    real(rp) :: lift_mat(14,14,14,6)
    real(rp) :: lift_mat_tr(14,14,14,6)
    real(rp) :: intrpmat_vpordm1_tr(14,14)
 
    real(rp) :: mfilter_h1d(14,14)
    real(rp) :: mfilter_h1d_tr(14,14)
    real(rp) :: mfilter_v1d_tr(14,14)
 
    real(rp) :: mfilter_tracer_h1d(14,14)
    real(rp) :: mfilter_tracer_h1d_tr(14,14)
    real(rp) :: mfilter_tracer_v1d_tr(14,14)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p13
    procedure, public :: final => element_operation_tensorprod3d_final_p13
    procedure, public :: dx => element_operation_tensorprod3d_dx_p13
    procedure, public :: dy => element_operation_tensorprod3d_dy_p13
    procedure, public :: dz => element_operation_tensorprod3d_dz_p13
    procedure, public :: lift => element_operation_tensorprod3d_lift_p13
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p13
    procedure, public :: div => element_operation_tensorprod3d_div_p13
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p13
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p13
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p13
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p13
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p13
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p13
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p13
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p13
  end type elementoperationtensorprod3d_p13
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=14
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p14
    real(rp) :: d1d(15,15)
    real(rp) :: d1d_tr(15,15)
    real(rp) :: lift_mat(15,15,15,6)
    real(rp) :: lift_mat_tr(15,15,15,6)
    real(rp) :: intrpmat_vpordm1_tr(15,15)
 
    real(rp) :: mfilter_h1d(15,15)
    real(rp) :: mfilter_h1d_tr(15,15)
    real(rp) :: mfilter_v1d_tr(15,15)
 
    real(rp) :: mfilter_tracer_h1d(15,15)
    real(rp) :: mfilter_tracer_h1d_tr(15,15)
    real(rp) :: mfilter_tracer_v1d_tr(15,15)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p14
    procedure, public :: final => element_operation_tensorprod3d_final_p14
    procedure, public :: dx => element_operation_tensorprod3d_dx_p14
    procedure, public :: dy => element_operation_tensorprod3d_dy_p14
    procedure, public :: dz => element_operation_tensorprod3d_dz_p14
    procedure, public :: lift => element_operation_tensorprod3d_lift_p14
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p14
    procedure, public :: div => element_operation_tensorprod3d_div_p14
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p14
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p14
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p14
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p14
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p14
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p14
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p14
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p14
  end type elementoperationtensorprod3d_p14
  !> Derived type for elementwise operations with 3D tensor product elements with polynomial order p=15
  type, public, extends(elementoperationtensorprod3d) :: elementoperationtensorprod3d_p15
    real(rp) :: d1d(16,16)
    real(rp) :: d1d_tr(16,16)
    real(rp) :: lift_mat(16,16,16,6)
    real(rp) :: lift_mat_tr(16,16,16,6)
    real(rp) :: intrpmat_vpordm1_tr(16,16)
 
    real(rp) :: mfilter_h1d(16,16)
    real(rp) :: mfilter_h1d_tr(16,16)
    real(rp) :: mfilter_v1d_tr(16,16)
 
    real(rp) :: mfilter_tracer_h1d(16,16)
    real(rp) :: mfilter_tracer_h1d_tr(16,16)
    real(rp) :: mfilter_tracer_v1d_tr(16,16)
  contains
    procedure, public :: init => element_operation_tensorprod3d_init_p15
    procedure, public :: final => element_operation_tensorprod3d_final_p15
    procedure, public :: dx => element_operation_tensorprod3d_dx_p15
    procedure, public :: dy => element_operation_tensorprod3d_dy_p15
    procedure, public :: dz => element_operation_tensorprod3d_dz_p15
    procedure, public :: lift => element_operation_tensorprod3d_lift_p15
    procedure, public :: dxdydzlift => element_operation_tensorprod3d_dxdydzlift_p15
    procedure, public :: div => element_operation_tensorprod3d_div_p15
    procedure, public :: div_var5 => element_operation_tensorprod3d_div_var5_p15
    procedure, public :: div_var5_2 => element_operation_tensorprod3d_div_var5_2_p15
    procedure, public :: lift_var5 => element_operation_tensorprod3d_lift_var5_p15
    procedure, public :: vfilterpm1 => element_operation_tensorprod3d_vfilterpm1_p15
    !-
    procedure, public :: setup_modalfilter => element_operation_tensorprod3d_setup_modalfilter_p15
    procedure, public :: setup_modalfilter_tracer => element_operation_tensorprod3d_setup_modalfilter_tracer_p15
    procedure, public :: modalfilter_tracer => element_operation_tensorprod3d_modalfilter_tracer_p15
    procedure, public :: modalfilter_var5 => element_operation_tensorprod3d_modalfilter_var5_p15
  end type elementoperationtensorprod3d_p15
 
contains
  !> Factory subroutine to create an object of ElementOperationTensorProd3D
!OCL SERIAL
  subroutine elementoperationtensorprod3d_create( elem3D, obj )
    use scale_prc, only: prc_abort
    implicit none
    class(elementbase3d), intent(in) :: elem3d                           !< A 3D element object whose polynomial order is used to determine the type of created object
    class(elementoperationtensorprod3d), intent(out), allocatable :: obj !< Created object of ElementOperationTensorProd3D
    !-----------------------------------------------------
 
    if ( elem3d%PolyOrder_h /= elem3d%PolyOrder_v ) then
      log_info('ElemenetOperationTensorprod3D_create',*) 'PolyOrder_h and PolyOrder_v should be same in this module. ElemenetOperationGeneral should be used!'
      call prc_abort
    end if
 
    select case( elem3d%PolyOrder_h )
    case(1)
      allocate( elementoperationtensorprod3d_p1 :: obj )
    case(2)
      allocate( elementoperationtensorprod3d_p2 :: obj )
    case(3)
      allocate( elementoperationtensorprod3d_p3 :: obj )
    case(4)
      allocate( elementoperationtensorprod3d_p4 :: obj )
    case(5)
      allocate( elementoperationtensorprod3d_p5 :: obj )
    case(6)
      allocate( elementoperationtensorprod3d_p6 :: obj )
    case(7)
      allocate( elementoperationtensorprod3d_p7 :: obj )
    case(8)
      allocate( elementoperationtensorprod3d_p8 :: obj )
    case(9)
      allocate( elementoperationtensorprod3d_p9 :: obj )
    case(10)
      allocate( elementoperationtensorprod3d_p10 :: obj )
    case(11)
      allocate( elementoperationtensorprod3d_p11 :: obj )
    case(12)
      allocate( elementoperationtensorprod3d_p12 :: obj )
    case(13)
      allocate( elementoperationtensorprod3d_p13 :: obj )
    case(14)
      allocate( elementoperationtensorprod3d_p14 :: obj )
    case(15)
      allocate( elementoperationtensorprod3d_p15 :: obj )
    case default
      log_info('ElemenetOperationTensorprod3D_create',*) 'The specified polynomial order execceds 15 and is not supported in this module. ElemenetOperationGeneral should be used!'
      call prc_abort
    end select
    call obj%Init( elem3d )
 
    return
  end subroutine elementoperationtensorprod3d_create
  
 
!--- For p=1 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p1( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p1), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      2, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p1
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p1( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p1), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p1 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p1( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p1), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p1 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p1( this )
    implicit none
    class(elementoperationtensorprod3d_p1), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p1
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p1( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p1( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p1
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p1( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p1( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p1
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p1( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p1( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p1
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p1( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p1( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p1
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p1( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p1, &
      element_operation_kernel_matvec_lift_hexahedral_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p1( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p1( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p1
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p1(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p1, &
      element_operation_kernel_matvec_lift_hexahedral_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p1( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p1( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p1
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p1(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p1, &
      element_operation_kernel_matvec_lift_hexahedral_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p1( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p1( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p1
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p1(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p1( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p1
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p1( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p1( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p1
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p1( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p1
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p1( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p1 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p1( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p1
 
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p1( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p1 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p1( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p1
 
    implicit none
    class(elementoperationtensorprod3d_p1), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p1( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p1 
 
 
!--- For p=2 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p2( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p2), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      3, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p2
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p2( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p2), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p2 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p2( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p2), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p2 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p2( this )
    implicit none
    class(elementoperationtensorprod3d_p2), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p2
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p2( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p2( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p2
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p2( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p2( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p2
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p2( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p2( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p2
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p2( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p2( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p2
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p2( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p2, &
      element_operation_kernel_matvec_lift_hexahedral_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p2( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p2( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p2
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p2(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p2, &
      element_operation_kernel_matvec_lift_hexahedral_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p2( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p2( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p2
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p2(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p2, &
      element_operation_kernel_matvec_lift_hexahedral_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p2( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p2( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p2
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p2(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p2( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p2
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p2( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p2( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p2
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p2( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p2
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p2( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p2 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p2( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p2
 
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p2( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p2 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p2( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p2
 
    implicit none
    class(elementoperationtensorprod3d_p2), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p2( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p2 
 
 
!--- For p=3 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p3( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p3), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      4, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p3
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p3( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p3), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p3 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p3( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p3), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p3 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p3( this )
    implicit none
    class(elementoperationtensorprod3d_p3), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p3
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p3( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p3( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p3
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p3( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p3( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p3
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p3( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p3( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p3
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p3( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p3( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p3
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p3( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p3, &
      element_operation_kernel_matvec_lift_hexahedral_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p3( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p3( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p3
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p3(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p3, &
      element_operation_kernel_matvec_lift_hexahedral_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p3( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p3( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p3
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p3(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p3, &
      element_operation_kernel_matvec_lift_hexahedral_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p3( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p3( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p3
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p3(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p3( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p3
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p3( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p3( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p3
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p3( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p3
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p3( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p3 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p3( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p3
 
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p3( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p3 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p3( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p3
 
    implicit none
    class(elementoperationtensorprod3d_p3), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p3( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p3 
 
 
!--- For p=4 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p4( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p4), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      5, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p4
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p4( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p4), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p4 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p4( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p4), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p4 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p4( this )
    implicit none
    class(elementoperationtensorprod3d_p4), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p4
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p4( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p4( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p4
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p4( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p4( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p4
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p4( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p4( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p4
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p4( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p4( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p4
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p4( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p4, &
      element_operation_kernel_matvec_lift_hexahedral_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p4( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p4( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p4
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p4(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p4, &
      element_operation_kernel_matvec_lift_hexahedral_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p4( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p4( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p4
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p4(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p4, &
      element_operation_kernel_matvec_lift_hexahedral_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p4( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p4( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p4
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p4(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p4( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p4
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p4( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p4( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p4
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p4( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p4
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p4( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p4 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p4( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p4
 
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p4( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p4 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p4( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p4
 
    implicit none
    class(elementoperationtensorprod3d_p4), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p4( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p4 
 
 
!--- For p=5 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p5( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p5), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      6, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p5
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p5( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p5), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p5 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p5( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p5), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p5 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p5( this )
    implicit none
    class(elementoperationtensorprod3d_p5), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p5
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p5( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p5( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p5
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p5( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p5( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p5
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p5( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p5( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p5
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p5( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p5( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p5
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p5( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p5, &
      element_operation_kernel_matvec_lift_hexahedral_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p5( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p5( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p5
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p5(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p5, &
      element_operation_kernel_matvec_lift_hexahedral_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p5( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p5( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p5
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p5(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p5, &
      element_operation_kernel_matvec_lift_hexahedral_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p5( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p5( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p5
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p5(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p5( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p5
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p5( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p5( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p5
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p5( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p5
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p5( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p5 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p5( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p5
 
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p5( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p5 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p5( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p5
 
    implicit none
    class(elementoperationtensorprod3d_p5), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p5( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p5 
 
 
!--- For p=6 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p6( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p6), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      7, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p6
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p6( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p6), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p6 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p6( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p6), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p6 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p6( this )
    implicit none
    class(elementoperationtensorprod3d_p6), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p6
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p6( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p6( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p6
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p6( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p6( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p6
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p6( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p6( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p6
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p6( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p6( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p6
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p6( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p6, &
      element_operation_kernel_matvec_lift_hexahedral_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p6( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p6( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p6
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p6(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p6, &
      element_operation_kernel_matvec_lift_hexahedral_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p6( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p6( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p6
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p6(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p6, &
      element_operation_kernel_matvec_lift_hexahedral_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p6( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p6( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p6
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p6(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p6( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p6
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p6( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p6( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p6
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p6( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p6
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p6( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p6 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p6( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p6
 
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p6( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p6 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p6( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p6
 
    implicit none
    class(elementoperationtensorprod3d_p6), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p6( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p6 
 
 
!--- For p=7 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p7( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p7), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      8, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p7
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p7( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p7), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p7 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p7( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p7), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p7 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p7( this )
    implicit none
    class(elementoperationtensorprod3d_p7), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p7
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p7( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p7( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p7
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p7( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p7( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p7
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p7( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p7( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p7
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p7( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p7( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p7
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p7( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p7, &
      element_operation_kernel_matvec_lift_hexahedral_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p7( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p7( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p7
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p7(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p7, &
      element_operation_kernel_matvec_lift_hexahedral_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p7( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p7( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p7
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p7(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p7, &
      element_operation_kernel_matvec_lift_hexahedral_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p7( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p7( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p7
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p7(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p7( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p7
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p7( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p7( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p7
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p7( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p7
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p7( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p7 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p7( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p7
 
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p7( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p7 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p7( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p7
 
    implicit none
    class(elementoperationtensorprod3d_p7), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p7( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p7 
 
 
!--- For p=8 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p8( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p8), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      9, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p8
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p8( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p8), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p8 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p8( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p8), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p8 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p8( this )
    implicit none
    class(elementoperationtensorprod3d_p8), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p8
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p8( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p8( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p8
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p8( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p8( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p8
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p8( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p8( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p8
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p8( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p8( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p8
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p8( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p8, &
      element_operation_kernel_matvec_lift_hexahedral_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p8( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p8( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p8
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p8(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p8, &
      element_operation_kernel_matvec_lift_hexahedral_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p8( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p8( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p8
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p8(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p8, &
      element_operation_kernel_matvec_lift_hexahedral_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p8( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p8( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p8
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p8(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p8( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p8
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p8( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p8( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p8
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p8( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p8
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p8( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p8 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p8( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p8
 
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p8( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p8 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p8( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p8
 
    implicit none
    class(elementoperationtensorprod3d_p8), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p8( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p8 
 
 
!--- For p=9 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p9( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p9), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      10, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p9
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p9( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p9), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p9 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p9( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p9), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p9 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p9( this )
    implicit none
    class(elementoperationtensorprod3d_p9), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p9
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p9( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p9( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p9
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p9( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p9( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p9
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p9( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p9( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p9
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p9( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p9( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p9
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p9( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p9, &
      element_operation_kernel_matvec_lift_hexahedral_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p9( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p9( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p9
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p9(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p9, &
      element_operation_kernel_matvec_lift_hexahedral_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p9( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p9( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p9
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p9(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p9, &
      element_operation_kernel_matvec_lift_hexahedral_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p9( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p9( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p9
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p9(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p9( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p9
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p9( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p9( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p9
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p9( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p9
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p9( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p9 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p9( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p9
 
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p9( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p9 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p9( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p9
 
    implicit none
    class(elementoperationtensorprod3d_p9), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p9( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p9 
 
 
!--- For p=10 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p10( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p10), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      11, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p10
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p10( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p10), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p10 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p10( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p10), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p10 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p10( this )
    implicit none
    class(elementoperationtensorprod3d_p10), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p10
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p10( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p10( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p10
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p10( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p10( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p10
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p10( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p10( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p10
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p10( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p10( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p10
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p10( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p10, &
      element_operation_kernel_matvec_lift_hexahedral_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p10( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p10( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p10
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p10(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p10, &
      element_operation_kernel_matvec_lift_hexahedral_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p10( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p10( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p10
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p10(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p10, &
      element_operation_kernel_matvec_lift_hexahedral_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p10( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p10( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p10
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p10(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p10( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p10
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p10( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p10( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p10
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p10( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p10
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p10( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p10 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p10( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p10
 
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p10( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p10 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p10( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p10
 
    implicit none
    class(elementoperationtensorprod3d_p10), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p10( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p10 
 
 
!--- For p=11 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p11( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p11), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      12, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p11
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p11( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p11), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p11 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p11( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p11), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p11 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p11( this )
    implicit none
    class(elementoperationtensorprod3d_p11), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p11
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p11( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p11( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p11
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p11( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p11( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p11
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p11( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p11( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p11
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p11( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p11( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p11
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p11( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p11, &
      element_operation_kernel_matvec_lift_hexahedral_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p11( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p11( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p11
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p11(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p11, &
      element_operation_kernel_matvec_lift_hexahedral_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p11( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p11( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p11
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p11(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p11, &
      element_operation_kernel_matvec_lift_hexahedral_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p11( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p11( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p11
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p11(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p11( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p11
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p11( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p11( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p11
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p11( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p11
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p11( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p11 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p11( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p11
 
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p11( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p11 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p11( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p11
 
    implicit none
    class(elementoperationtensorprod3d_p11), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p11( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p11 
 
 
!--- For p=12 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p12( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p12), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      13, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p12
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p12( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p12), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p12 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p12( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p12), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p12 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p12( this )
    implicit none
    class(elementoperationtensorprod3d_p12), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p12
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p12( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p12( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p12
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p12( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p12( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p12
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p12( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p12( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p12
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p12( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p12( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p12
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p12( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p12, &
      element_operation_kernel_matvec_lift_hexahedral_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p12( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p12( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p12
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p12(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p12, &
      element_operation_kernel_matvec_lift_hexahedral_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p12( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p12( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p12
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p12(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p12, &
      element_operation_kernel_matvec_lift_hexahedral_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p12( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p12( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p12
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p12(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p12( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p12
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p12( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p12( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p12
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p12( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p12
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p12( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p12 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p12( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p12
 
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p12( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p12 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p12( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p12
 
    implicit none
    class(elementoperationtensorprod3d_p12), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p12( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p12 
 
 
!--- For p=13 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p13( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p13), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      14, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p13
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p13( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p13), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p13 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p13( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p13), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p13 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p13( this )
    implicit none
    class(elementoperationtensorprod3d_p13), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p13
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p13( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p13( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p13
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p13( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p13( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p13
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p13( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p13( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p13
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p13( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p13( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p13
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p13( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p13, &
      element_operation_kernel_matvec_lift_hexahedral_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p13( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p13( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p13
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p13(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p13, &
      element_operation_kernel_matvec_lift_hexahedral_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p13( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p13( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p13
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p13(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p13, &
      element_operation_kernel_matvec_lift_hexahedral_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p13( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p13( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p13
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p13(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p13( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p13
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p13( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p13( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p13
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p13( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p13
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p13( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p13 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p13( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p13
 
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p13( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p13 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p13( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p13
 
    implicit none
    class(elementoperationtensorprod3d_p13), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p13( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p13 
 
 
!--- For p=14 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p14( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p14), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      15, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p14
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p14( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p14), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p14 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p14( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p14), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p14 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p14( this )
    implicit none
    class(elementoperationtensorprod3d_p14), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p14
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p14( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p14( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p14
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p14( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p14( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p14
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p14( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p14( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p14
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p14( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p14( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p14
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p14( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p14, &
      element_operation_kernel_matvec_lift_hexahedral_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p14( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p14( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p14
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p14(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p14, &
      element_operation_kernel_matvec_lift_hexahedral_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p14( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p14( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p14
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p14(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p14, &
      element_operation_kernel_matvec_lift_hexahedral_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p14( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p14( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p14
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p14(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p14( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p14
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p14( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p14( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p14
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p14( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p14
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p14( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p14 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p14( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p14
 
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p14( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p14 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p14( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p14
 
    implicit none
    class(elementoperationtensorprod3d_p14), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p14( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p14 
 
 
!--- For p=15 ------------------------------------
 
!> Initialization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_init_p15( this, elem3D )
    use scale_element_operation_base, only: element_operator_type_tensorprod3d
    implicit none
    class(elementoperationtensorprod3d_p15), intent(inout) :: this
    class(elementbase3d), intent(in), target :: elem3d
    !----------------------------------------------------------
 
    this%elem3D => elem3d
    this%operator_type = element_operator_type_tensorprod3d
    
    call setup_elem_operator( this%D1D, this%D1D_tr, this%Lift_mat, this%Lift_mat_tr, this%IntrpMat_VPOrdM1_tr, &
      16, elem3d )
    !$acc enter data copyin( this )
    !$acc enter data attach( this%elem3D )
    
    return
  end subroutine element_operation_tensorprod3d_init_p15
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_p15( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p15), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_p15 
 
!> Setup modal filter for tracer
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p15( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    use scale_element_hexahedral, only: hexahedralelement
    use scale_element_line, only: lineelement
    implicit none
    class(elementoperationtensorprod3d_p15), intent(inout) :: this
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, &
      mf_etac_h, mf_alpha_h, mf_order_h,               &
      mf_etac_v, mf_alpha_v, mf_order_v,               &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
    !$acc update device( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr )
    return
  end subroutine element_operation_tensorprod3d_setup_modalfilter_tracer_p15 
 
!> Finalization
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_final_p15( this )
    implicit none
    class(elementoperationtensorprod3d_p15), intent(inout) :: this
    !----------------------------------------------------------
    nullify( this%elem3D )
    !$acc exit data delete( this )
    return
  end subroutine element_operation_tensorprod3d_final_p15
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dx_p15( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirx_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirx_p15( this%D1D, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dx_p15
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dy_p15( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_diry_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_diry_p15( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dy_p15
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dz_p15( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_dirz_p15( this%D1D_tr, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_dz_p15
 
!> Calculate the lifting term
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_p15( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_lift_hexahedral_p15( this%Lift_mat, vec_in, vec_out )
    return
  end subroutine element_operation_tensorprod3d_lift_p15
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_dxdydzlift_p15( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_gradlike_dirxyz_p15, &
      element_operation_kernel_matvec_lift_hexahedral_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out_dx(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dy(this%elem3d%np)
    real(rp), intent(out) :: vec_out_dz(this%elem3d%np)
    real(rp), intent(out) :: vec_out_lift(this%elem3d%np)
    !----------------------------------------------------------
    call element_operation_kernel_matvec_gradlike_dirxyz_p15( this%D1D, this%D1D_tr, vec_in, vec_in, &
      vec_out_dx, vec_out_dy, vec_out_dz )
    call element_operation_kernel_matvec_lift_hexahedral_p15( this%Lift_mat, vec_in_lift, &
      vec_out_lift )      
 
    return
  end subroutine element_operation_tensorprod3d_dxdydzlift_p15
 
!> Calculate the 3D divergence
!!
!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_p15(  this, vec_in, vec_in_lift, &
    vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p15, &
      element_operation_kernel_matvec_lift_hexahedral_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot)
    real(rp), intent(out) :: vec_out(this%elem3d%np,4)
    !----------------------------------------------------------
 
    call element_operation_kernel_matvec_lift_hexahedral_p15( this%Lift_mat, vec_in_lift, &
      vec_out(:,4) )      
 
    call element_operation_kernel_matvec_divlike_dirxyz_p15( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
    return
  end subroutine element_operation_tensorprod3d_div_p15
 
!> Calculate the 3D gradient applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_p15(  this, vec_in, vec_in_lift, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p15, &
      element_operation_kernel_matvec_lift_hexahedral_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(in) :: vec_in_lift(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,4,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p15( this%Lift_mat, vec_in_lift(:,iv), &
        vec_out_d(:,4,iv) )      
    end do
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p15( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_p15
 
!> Calculate the 3D divergence applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_div_var5_2_p15(  this, vec_in, &
    vec_out_d )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_divlike_dirxyz_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,3,5)
    real(rp), intent(out) :: vec_out_d(this%elem3d%np,3,5)
 
    integer :: iv
    !----------------------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_divlike_dirxyz_p15( this%D1D, this%D1D_tr, vec_in(:,1,iv), vec_in(:,2,iv), vec_in(:,3,iv), &
        vec_out_d(:,1,iv), vec_out_d(:,2,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_div_var5_2_p15
 
!> Calculate the lifting term applied for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_lift_var5_p15( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_lift_hexahedral_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%nfptot,5)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1, 5
      call element_operation_kernel_matvec_lift_hexahedral_p15( this%Lift_mat, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_tensorprod3d_lift_var5_p15
  
!> Apply a vertical modal filter which a coefficient with the highest vertical mode is set to zero
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_vfilterpm1_p15( this, vec_in, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: element_operation_kernel_matvec_dirz_p15
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------------------------
 
    call element_operation_kernel_matvec_dirz_p15( this%IntrpMat_VPOrdM1_tr, vec_in, &
      vec_out )
    return
  end subroutine element_operation_tensorprod3d_vfilterpm1_p15 
 
!> Apply a modal filter for tracer variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_tracer_p15( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p15
 
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np)
    !---------------------------------------------
 
    call element_operation_kernel_matvec_modalfilter_p15( this%MFilter_tracer_h1D, this%MFilter_tracer_h1D_tr, this%MFilter_tracer_v1D_tr, vec_in(:), vec_work, &
      vec_out(:) )
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_tracer_p15 
  
!> Apply a modal filter for five variables
!!
!OCL SERIAL
  subroutine element_operation_tensorprod3d_modalfilter_var5_p15( this, vec_in, vec_work, vec_out )
    use scale_element_operation_tensorprod3d_kernel, only: &
      element_operation_kernel_matvec_modalfilter_p15
 
    implicit none
    class(elementoperationtensorprod3d_p15), intent(in) :: this
    real(rp), intent(in) :: vec_in(this%elem3d%np,5)
    real(rp), intent(out) :: vec_work(this%elem3d%np)
    real(rp), intent(out) :: vec_out(this%elem3d%np,5)
    integer :: iv
    !---------------------------------------------
 
    do iv=1, 5
      call element_operation_kernel_matvec_modalfilter_p15( this%MFilter_h1D, this%MFilter_h1D_tr, this%MFilter_v1D_tr, vec_in(:,iv), vec_work, &
        vec_out(:,iv) )
    end do
    
    return
  end subroutine element_operation_tensorprod3d_modalfilter_var5_p15 
 
 
!- private -
 
!OCL SERIAL
  subroutine setup_modalfilter( MFilter_h1D, MFilter_h1D_tr, MFilter_v1D_tr, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v, &
    PolyOrder_h, PolyOrder_v )
 
    use scale_element_line, only: lineelement
    implicit none
 
    real(rp), intent(inout) :: mfilter_h1d(:,:)
    real(rp), intent(inout) :: mfilter_h1d_tr(:,:)
    real(rp), intent(inout) :: mfilter_v1d_tr(:,:)
    real(rp), intent(in) :: mf_etac_h
    real(rp), intent(in) :: mf_alpha_h
    integer, intent(in) :: mf_order_h
    real(rp), intent(in) :: mf_etac_v
    real(rp), intent(in) :: mf_alpha_v
    integer, intent(in) :: mf_order_v
    integer, intent(in) :: polyorder_h
    integer, intent(in) :: polyorder_v
 
    type(lineelement) :: elem1d
    type(modalfilter) :: mfilter1d
    !--------------------------------------------------------
 
    call elem1d%Init( polyorder_h, .false. )      
 
    call mfilter1d%Init( elem1d, &
      mf_etac_h, mf_alpha_h, mf_order_h      )
    mfilter_h1d(:,:) = mfilter1d%FilterMat(:,:)
    mfilter_h1d_tr(:,:) = transpose(mfilter_h1d(:,:))
 
    call mfilter1d%Final()
    call elem1d%Final()
 
    !-
    call elem1d%Init( polyorder_v, .false. )      
 
    call mfilter1d%Init( elem1d, &
      mf_etac_v, mf_alpha_v, mf_order_v      )
    mfilter_v1d_tr(:,:) = transpose(mfilter1d%FilterMat(:,:))
 
    call mfilter1d%Final()
    call elem1d%Final()
 
    return
  end subroutine setup_modalfilter
 
!OCL SERIAL
  subroutine setup_elem_operator( D1D, D1D_tr, Lift_mat, Lift_mat_tr, IntrpMat_VPOrdM1_tr, &
    np, elem3D )
    use scale_element_line, only: lineelement
    use scale_sparsemat, only: sparsemat
    implicit none
    integer, intent(in) :: np
    class(elementbase3d), intent(in) :: elem3d
    real(rp), intent(out) :: d1d(np,np)
    real(rp), intent(out) :: d1d_tr(np,np)
    real(rp), intent(out) :: lift_mat(np,np,np,6)
    real(rp), intent(out) :: lift_mat_tr(6,np,np,np)
    real(rp), intent(out) :: intrpmat_vpordm1_tr(np,np)
 
    integer :: p
 
    integer :: p1, p2, p_
    real(rp) :: invv_vpordm1(np,np)
    real(rp) :: intrpmat_vpordm1(np,np)
 
    type(lineelement) :: elem1d
    type(sparsemat) :: lift_sm
    integer :: i, j, k
    !----------------------------------------------------------
 
    call elem1d%Init( np-1, .false. )
    d1d(:,:) = elem1d%Dx1(:,:)
    d1d_tr(:,:) = transpose(d1d)
 
    call lift_sm%Init( elem3d%Lift, storage_format='ELL' )
    do k=1, np
    do j=1, np
    do i=1, np
      p1 = i + (j-1)*np + (k-1)*np**2
      lift_mat_tr(1:6,i,j,k) = &
        (/ elem3d%Lift(p1,i+(k-1)*np), &
           elem3d%Lift(p1,j+(k-1)*np+(2-1)*np**2), &
           elem3d%Lift(p1,i+(k-1)*np+(3-1)*np**2), &
           elem3d%Lift(p1,j+(k-1)*np+(4-1)*np**2), &
           elem3d%Lift(p1,i+(j-1)*np+(5-1)*np**2), &
           elem3d%Lift(p1,i+(j-1)*np+(6-1)*np**2) /)
      lift_mat(i,j,k,1:6) = lift_mat_tr(1:6,i,j,k)
    end do
    end do
    end do
    call lift_sm%Final()
 
    !--
 
    invv_vpordm1(:,:) = elem1d%invV(:,:)
    invv_vpordm1(np,:) = 0.0_rp
    intrpmat_vpordm1(:,:) = matmul(elem1d%V, invv_vpordm1)
    intrpmat_vpordm1_tr(:,:) = transpose(intrpmat_vpordm1)
 
    !-
    call elem1d%Final()
 
    return
  end subroutine setup_elem_operator
 
end module scale_element_operation_tensorprod3d