diff --git a/src/interpolation/CMakeLists.txt b/src/interpolation/CMakeLists.txt
index 37c3ad0997aadfb480d01f544e94a7655f8482e1..10515165cbe3a16354226aa9f11756200bfa33fd 100644
--- a/src/interpolation/CMakeLists.txt
+++ b/src/interpolation/CMakeLists.txt
@@ -12,11 +12,13 @@
add_library(
iconmath-interpolation
mo_lib_interpolation_scalar.F90
+ mo_lib_interpolation_scalar.cpp
mo_lib_interpolation_vector.F90
mo_lib_interpolation_vector.cpp
mo_lib_intp_rbf.F90
mo_lib_intp_rbf-rbf_vec_interpol_vertex_lib.cpp
mo_lib_intp_rbf-rbf_vec_interpol_vertex_lib_bindings.cpp
+ interpolation_bindings.cpp
)
add_library(${PROJECT_NAME}::interpolation ALIAS iconmath-interpolation)
diff --git a/src/interpolation/interpolation_bindings.cpp b/src/interpolation/interpolation_bindings.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..628f4115978bc8981084b3418478d8e1844b2232
--- /dev/null
+++ b/src/interpolation/interpolation_bindings.cpp
@@ -0,0 +1,328 @@
+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#include "interpolation_bindings.h"
+#include "mo_lib_interpolation_scalar.hpp"
+#include "mo_lib_interpolation_vector.hpp"
+
+// This is the binding for mo_interpolation_vector::edges2cells_vector_lib
+// (wp=dp)
+void edges2cells_vector_lib_dp(const double *p_vn_in, const double *p_vt_in,
+ const int *cell_edge_idx,
+ const int *cell_edge_blk,
+ const double *e_bln_c_u, const double *e_bln_c_v,
+ double *p_u_out, double *p_v_out, int i_startblk,
+ int i_endblk, int i_startidx_in, int i_endidx_in,
+ int slev, int elev, int nproma, int nlev,
+ int nblks_e, int nblks_c) {
+
+ edges2cells_vector_lib<double>(
+ p_vn_in, p_vt_in, cell_edge_idx, cell_edge_blk, e_bln_c_u, e_bln_c_v,
+ p_u_out, p_v_out, i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev,
+ elev, nproma, nlev, nblks_e, nblks_c);
+}
+
+// This is the binding for mo_interpolation_vector::edges2cells_vector_lib
+// (wp=sp)
+void edges2cells_vector_lib_sp(const float *p_vn_in, const float *p_vt_in,
+ const int *cell_edge_idx,
+ const int *cell_edge_blk, const float *e_bln_c_u,
+ const float *e_bln_c_v, float *p_u_out,
+ float *p_v_out, int i_startblk, int i_endblk,
+ int i_startidx_in, int i_endidx_in, int slev,
+ int elev, int nproma, int nlev, int nblks_e,
+ int nblks_c) {
+
+ edges2cells_vector_lib<float>(
+ p_vn_in, p_vt_in, cell_edge_idx, cell_edge_blk, e_bln_c_u, e_bln_c_v,
+ p_u_out, p_v_out, i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev,
+ elev, nproma, nlev, nblks_e, nblks_c);
+}
+
+// This is the binding for mo_interpolation_scalar::verts2edges_scalar_lib
+// (wp=dp)
+void verts2edges_scalar_lib_dp(
+ const double *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const double *coeff_int, double *p_edge_out,
+ const int i_startblk, const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_v, const int nblks_e, const bool lacc) {
+
+ verts2edges_scalar_lib<double>(p_vertex_in, edge_vertex_idx, edge_vertex_blk,
+ coeff_int, p_edge_out, i_startblk, i_endblk,
+ i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_v, nblks_e, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::verts2edges_scalar_lib
+// (wp=sp)
+void verts2edges_scalar_lib_sp(
+ const float *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const float *coeff_int, float *p_edge_out,
+ const int i_startblk, const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_v, const int nblks_e, const bool lacc) {
+
+ verts2edges_scalar_lib<float>(p_vertex_in, edge_vertex_idx, edge_vertex_blk,
+ coeff_int, p_edge_out, i_startblk, i_endblk,
+ i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_v, nblks_e, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2edges_scalar_dp_lib
+void cells2edges_scalar_lib_dp(
+ const double *p_cell_in, const int *edge_cell_idx, const int *edge_cell_blk,
+ const double *coeff_int, double *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in, const int *i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e, const int patch_id,
+ const bool l_limited_area, const bool lfill_latbc, const bool lacc) {
+
+ cells2edges_scalar_lib<double, double>(
+ p_cell_in, edge_cell_idx, edge_cell_blk, coeff_int, p_edge_out,
+ i_startblk_in, i_endblk_in, i_startidx_in, i_endidx_in, slev, elev,
+ nproma, nlev, nblk_c, nblks_e, patch_id, l_limited_area, lfill_latbc,
+ lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2edges_scalar_sp_lib
+void cells2edges_scalar_lib_sp(const float *p_cell_in, const int *edge_cell_idx,
+ const int *edge_cell_blk, const float *coeff_int,
+ float *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in,
+ const int *i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e,
+ const int patch_id, const bool l_limited_area,
+ const bool lfill_latbc, const bool lacc) {
+
+ cells2edges_scalar_lib<float, float>(
+ p_cell_in, edge_cell_idx, edge_cell_blk, coeff_int, p_edge_out,
+ i_startblk_in, i_endblk_in, i_startidx_in, i_endidx_in, slev, elev,
+ nproma, nlev, nblk_c, nblks_e, patch_id, l_limited_area, lfill_latbc,
+ lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2edges_scalar_sp2dp_lib
+void cells2edges_scalar_lib_sp2dp(
+ const float *p_cell_in, const int *edge_cell_idx, const int *edge_cell_blk,
+ const double *coeff_int, double *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in, const int *i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e, const int patch_id,
+ const bool l_limited_area, const bool lfill_latbc, const bool lacc) {
+
+ cells2edges_scalar_lib<float, double>(
+ p_cell_in, edge_cell_idx, edge_cell_blk, coeff_int, p_edge_out,
+ i_startblk_in, i_endblk_in, i_startidx_in, i_endidx_in, slev, elev,
+ nproma, nlev, nblk_c, nblks_e, patch_id, l_limited_area, lfill_latbc,
+ lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::edges2verts_scalar_lib
+// (wp=dp)
+void edges2verts_scalar_lib_dp(
+ const double *p_edge_in, const int *vert_edge_idx, const int *vert_edge_blk,
+ const double *v_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v, const bool lacc) {
+
+ edges2verts_scalar_lib<double>(p_edge_in, vert_edge_idx, vert_edge_blk, v_int,
+ p_vert_out, i_startblk, i_endblk,
+ i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_e, nblks_v, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::edges2verts_scalar_lib
+// (wp=sp)
+void edges2verts_scalar_lib_sp(const float *p_edge_in, const int *vert_edge_idx,
+ const int *vert_edge_blk, const float *v_int,
+ float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v,
+ const bool lacc) {
+
+ edges2verts_scalar_lib<float>(p_edge_in, vert_edge_idx, vert_edge_blk, v_int,
+ p_vert_out, i_startblk, i_endblk, i_startidx_in,
+ i_endidx_in, slev, elev, nproma, nlev, nblks_e,
+ nblks_v, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::edges2cells_scalar_dp_lib
+void edges2cells_scalar_lib_dp(const double *p_edge_in, const int *edge_idx,
+ const int *edge_blk, const double *coeff_int,
+ double *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c,
+ const bool lacc) {
+ edges2cells_scalar_lib<double>(p_edge_in, edge_idx, edge_blk, coeff_int,
+ p_cell_out, i_startblk, i_endblk,
+ i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_e, nblks_c, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::edges2cells_scalar_sp_lib
+void edges2cells_scalar_lib_sp(const float *p_edge_in, const int *edge_idx,
+ const int *edge_blk, const float *coeff_int,
+ float *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c,
+ const bool lacc) {
+
+ edges2cells_scalar_lib<float>(p_edge_in, edge_idx, edge_blk, coeff_int,
+ p_cell_out, i_startblk, i_endblk, i_startidx_in,
+ i_endidx_in, slev, elev, nproma, nlev, nblks_e,
+ nblks_c, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2verts_scalar_dp_lib
+void cells2verts_scalar_lib_dp(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async) {
+ cells2verts_scalar_lib<double, double>(
+ p_cell_in, vert_cell_idx, vert_cell_blk, coeff_int, p_vert_out,
+ i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_c, nblks_v, lacc, acc_async);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2verts_scalar_dp2sp_lib
+void cells2verts_scalar_lib_dp2sp(
+ const float *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async) {
+ cells2verts_scalar_lib<float, double>(
+ p_cell_in, vert_cell_idx, vert_cell_blk, coeff_int, p_vert_out,
+ i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_c, nblks_v, lacc, acc_async);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2verts_scalar_sp_lib
+void cells2verts_scalar_lib_sp(const float *p_cell_in, const int *vert_cell_idx,
+ const int *vert_cell_blk, const float *coeff_int,
+ float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v,
+ const bool lacc, const bool acc_async) {
+ cells2verts_scalar_lib<float, float>(
+ p_cell_in, vert_cell_idx, vert_cell_blk, coeff_int, p_vert_out,
+ i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_c, nblks_v, lacc, acc_async);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2verts_scalar_ri_lib
+// (wp=dp, vp=dp)
+void cells2verts_scalar_ri_lib_dp(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async) {
+ cells2verts_scalar_ri_lib<double, double>(
+ p_cell_in, vert_cell_idx, vert_cell_blk, coeff_int, p_vert_out,
+ i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_c, nblks_v, lacc, acc_async);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2verts_scalar_ri_lib
+// (wp=dp, vp=sp)
+void cells2verts_scalar_ri_lib_dp2sp(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async) {
+ cells2verts_scalar_ri_lib<double, float>(
+ p_cell_in, vert_cell_idx, vert_cell_blk, coeff_int, p_vert_out,
+ i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_c, nblks_v, lacc, acc_async);
+}
+
+// This is the binding for mo_interpolation_scalar::cells2verts_scalar_ri_lib
+// (wp=sp, vp=sp)
+void cells2verts_scalar_ri_lib_sp(
+ const float *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const float *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async) {
+ cells2verts_scalar_ri_lib<float, float>(
+ p_cell_in, vert_cell_idx, vert_cell_blk, coeff_int, p_vert_out,
+ i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev, elev, nproma,
+ nlev, nblks_c, nblks_v, lacc, acc_async);
+}
+
+// This is the binding for mo_interpolation_scalar::verts2cells_scalar_lib
+// (wp=dp)
+void verts2cells_scalar_lib_dp(
+ const double *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const double *coeff_int, double *p_cell_out,
+ const int nblks_c, const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v, const bool lacc) {
+ verts2cells_scalar_lib<double>(p_vert_in, cell_index_idx, cell_vertex_blk,
+ coeff_int, p_cell_out, nblks_c, npromz_c, slev,
+ elev, nproma, nlev, nblks_v, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::verts2cells_scalar_lib
+// (wp=sp)
+void verts2cells_scalar_lib_sp(
+ const float *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const float *coeff_int, float *p_cell_out,
+ const int nblks_c, const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v, const bool lacc) {
+ verts2cells_scalar_lib<float>(p_vert_in, cell_index_idx, cell_vertex_blk,
+ coeff_int, p_cell_out, nblks_c, npromz_c, slev,
+ elev, nproma, nlev, nblks_v, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::cell_avg_lib (wp=dp)
+void cell_avg_lib_dp(const double *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const double *avg_coeff,
+ double *avg_psi_c, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_c,
+ const bool lacc) {
+ cell_avg_lib<double>(psi_c, cell_neighbor_idx, cell_neighbor_blk, avg_coeff,
+ avg_psi_c, i_startblk, i_endblk, i_startidx_in,
+ i_endidx_in, slev, elev, nproma, nlev, nblks_c, lacc);
+}
+
+// This is the binding for mo_interpolation_scalar::cell_avg_lib (wp=sp)
+void cell_avg_lib_sp(const float *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const float *avg_coeff,
+ float *avg_psi_c, const int i_startblk, const int i_endblk,
+ const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_c, const bool lacc) {
+ cell_avg_lib<float>(psi_c, cell_neighbor_idx, cell_neighbor_blk, avg_coeff,
+ avg_psi_c, i_startblk, i_endblk, i_startidx_in,
+ i_endidx_in, slev, elev, nproma, nlev, nblks_c, lacc);
+}
diff --git a/src/interpolation/interpolation_bindings.h b/src/interpolation/interpolation_bindings.h
new file mode 100644
index 0000000000000000000000000000000000000000..7cb873d8ac2f8ba1bffa0553d28847b2b9b11c58
--- /dev/null
+++ b/src/interpolation/interpolation_bindings.h
@@ -0,0 +1,188 @@
+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+#pragma once
+
+extern "C" {
+
+// mo_lib_interpolation_vector.F90
+void edges2cells_vector_lib_dp(const double *p_vn_in, const double *p_vt_in,
+ const int *cell_edge_idx,
+ const int *cell_edge_blk,
+ const double *e_bln_c_u, const double *e_bln_c_v,
+ double *p_u_out, double *p_v_out, int i_startblk,
+ int i_endblk, int i_startidx_in, int i_endidx_in,
+ int slev, int elev, int nproma, int nlev,
+ int nblks_e, int nblks_c);
+
+void edges2cells_vector_lib_sp(const float *p_vn_in, const float *p_vt_in,
+ const int *cell_edge_idx,
+ const int *cell_edge_blk, const float *e_bln_c_u,
+ const float *e_bln_c_v, float *p_u_out,
+ float *p_v_out, int i_startblk, int i_endblk,
+ int i_startidx_in, int i_endidx_in, int slev,
+ int elev, int nproma, int nlev, int nblks_e,
+ int nblks_c);
+
+// mo_lib_interpolation_scalar.F90
+void verts2edges_scalar_lib_dp(
+ const double *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const double *coeff_int, double *p_edge_out,
+ const int i_startblk, const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_v, const int nblks_e, const bool lacc);
+
+void verts2edges_scalar_lib_sp(
+ const float *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const float *coeff_int, float *p_edge_out,
+ const int i_startblk, const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_v, const int nblks_e, const bool lacc);
+
+void cells2edges_scalar_lib_dp(
+ const double *p_cell_in, const int *edge_cell_idx, const int *edge_cell_blk,
+ const double *coeff_int, double *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in, const int *i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e, const int patch_id,
+ const bool l_limited_area, const bool lfill_latbc, const bool lacc);
+
+void cells2edges_scalar_lib_sp(const float *p_cell_in, const int *edge_cell_idx,
+ const int *edge_cell_blk, const float *coeff_int,
+ float *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in,
+ const int *i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e,
+ const int patch_id, const bool l_limited_area,
+ const bool lfill_latbc, const bool lacc);
+
+void cells2edges_scalar_lib_sp2dp(
+ const float *p_cell_in, const int *edge_cell_idx, const int *edge_cell_blk,
+ const double *coeff_int, double *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in, const int *i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e, const int patch_id,
+ const bool l_limited_area, const bool lfill_latbc, const bool lacc);
+
+void edges2verts_scalar_lib_dp(
+ const double *p_edge_in, const int *vert_edge_idx, const int *vert_edge_blk,
+ const double *v_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v, const bool lacc);
+
+void edges2verts_scalar_lib_sp(const float *p_edge_in, const int *vert_edge_idx,
+ const int *vert_edge_blk, const float *v_int,
+ float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v,
+ const bool lacc);
+
+void edges2cells_scalar_lib_dp(const double *p_edge_in, const int *edge_idx,
+ const int *edge_blk, const double *coeff_int,
+ double *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c,
+ const bool lacc);
+void edges2cells_scalar_lib_sp(const float *p_edge_in, const int *edge_idx,
+ const int *edge_blk, const float *coeff_int,
+ float *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c,
+ const bool lacc);
+
+/////////////////////////////////////////////
+
+void cells2verts_scalar_lib_dp(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+void cells2verts_scalar_lib_dp2sp(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const float *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+void cells2verts_scalar_lib_sp(const float *p_cell_in, const int *vert_cell_idx,
+ const int *vert_cell_blk, const float *coeff_int,
+ float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v,
+ const bool lacc, const bool acc_async);
+
+/////////////////////////////////////////////
+
+void cells2verts_scalar_ri_lib_dp(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+void cells2verts_scalar_ri_lib_dp2sp(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+void cells2verts_scalar_ri_lib_sp(
+ const float *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const float *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+/////////////////////////////////////////////
+
+void verts2cells_scalar_lib_dp(
+ const double *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const double *coeff_int, double *p_cell_out,
+ const int nblks_c, const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v, const bool lacc);
+
+void verts2cells_scalar_lib_sp(
+ const float *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const float *coeff_int, float *p_cell_out,
+ const int nblks_c, const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v, const bool lacc);
+
+/////////////////////////////////////////////
+
+void cell_avg_lib_dp(const double *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const double *avg_coeff,
+ double *avg_psi_c, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_c,
+ const bool lacc);
+void cell_avg_lib_sp(const float *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const float *avg_coeff,
+ float *avg_psi_c, const int i_startblk, const int i_endblk,
+ const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_c, const bool lacc);
+}
diff --git a/src/interpolation/mo_lib_interpolation_scalar.cpp b/src/interpolation/mo_lib_interpolation_scalar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9e4e6c5ab4a3a531cb0cad254cc46b049c5ee58f
--- /dev/null
+++ b/src/interpolation/mo_lib_interpolation_scalar.cpp
@@ -0,0 +1,753 @@
+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#include "mo_lib_interpolation_scalar.hpp"
+#include "mo_lib_loopindices.hpp"
+#include <Kokkos_Core.hpp>
+#include <iostream>
+
+//-----------------------------------------------------------------------
+//
+// ! averaging and interpolation routines and
+// ! routines needed to compute the coefficients therein
+//
+//-----------------------------------------------------------------------
+
+//-----------------------------------------------------------------------
+//>
+/// Performs average of scalar fields from vertices to velocity points.
+///
+/// The coefficients are given by coeff_int.
+///
+template <typename T>
+void verts2edges_scalar_lib(const T *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const T *coeff_int,
+ T *p_edge_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_v, const int nblks_e,
+ const bool lacc) {
+
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedT3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ UnmanagedConstT3D p_vertex_in_view(p_vertex_in, nproma, nlev, nblks_v);
+ UnmanagedConstInt3D iidx_view(edge_vertex_idx, nproma, nblks_e, 4);
+ UnmanagedConstInt3D iblk_view(edge_vertex_blk, nproma, nblks_e, 4);
+ UnmanagedConstT3D coeff_int_view(coeff_int, nproma, 2, nblks_e);
+ UnmanagedT3D p_edge_out_view(p_edge_out, nproma, nlev, nblks_e);
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+
+ int i_startidx, i_endidx;
+ get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, jb, i_startblk,
+ i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "verts2edges_scalar", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int je) {
+ p_edge_out_view(je, jk, jb) =
+ coeff_int_view(je, 0, jb) *
+ p_vertex_in_view(iidx_view(je, jb, 0), jk,
+ iblk_view(je, jb, 0)) +
+ coeff_int_view(je, 1, jb) *
+ p_vertex_in_view(iidx_view(je, jb, 1), jk,
+ iblk_view(je, jb, 1));
+ });
+ Kokkos::fence();
+ }
+}
+
+//------------------------------------------------------------------------
+//>
+/// Computes average of scalar fields from centers of triangular faces to.
+///
+/// Computes average of scalar fields from centers of triangular faces to
+/// velocity points.
+///
+template <typename T, typename S>
+void cells2edges_scalar_lib(const T *p_cell_in, const int *edge_cell_idx,
+ const int *edge_cell_blk, const S *coeff_int,
+ S *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in,
+ const int *i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e,
+ const int patch_id, const bool l_limited_area,
+ const bool lfill_latbc, const bool lacc) {
+
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<const S ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstS3D;
+ typedef Kokkos::View<S ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedS3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ UnmanagedConstT3D p_cell_in_view(p_cell_in, nproma, nlev, nblk_c);
+ UnmanagedConstInt3D iidx_view(edge_cell_idx, nproma, nblks_e, 2);
+ UnmanagedConstInt3D iblk_view(edge_cell_blk, nproma, nblks_e, 2);
+ UnmanagedConstS3D coeff_int_view(coeff_int, nproma, 2, nblks_e);
+ UnmanagedS3D p_edge_out_view(p_edge_out, nproma, nlev, nblks_e);
+
+ // Fill outermost nest boundary
+ int i_startblk, i_endblk;
+ if ((l_limited_area || patch_id > 0) && (lfill_latbc)) {
+ i_startblk = i_startblk_in[0];
+ i_endblk = i_endblk_in[0];
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+
+ int i_startidx, i_endidx;
+ get_indices_e_lib(i_startidx_in[0], i_endidx_in[0], nproma, jb,
+ i_startblk, i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "cells2edges_scalar", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int je) {
+ if (iidx_view(je, jb, 0) >= 0 && iblk_view(je, jb, 0) >= 0) {
+ p_edge_out_view(je, jk, jb) = p_cell_in_view(
+ iidx_view(je, jb, 0), jk, iblk_view(je, jb, 0));
+ } else if (iidx_view(je, jb, 1) >= 0 && iblk_view(je, jb, 1) >= 0) {
+ p_edge_out_view(je, jk, jb) = p_cell_in_view(
+ iidx_view(je, jb, 1), jk, iblk_view(je, jb, 1));
+ } else {
+ std::cerr << "mo_interpolation:cells2edges_scalar_lib: error in "
+ "lateral boundary filling"
+ << std::endl;
+ std::exit(EXIT_FAILURE);
+ }
+ });
+ Kokkos::fence();
+ }
+ } else {
+ // Process the remaining grid points for which a real interpolation is
+ // possible
+ i_startblk = i_startblk_in[1];
+ i_endblk = i_endblk_in[1];
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+
+ int i_startidx, i_endidx;
+ get_indices_e_lib(i_startidx_in[1], i_endidx_in[1], nproma, jb,
+ i_startblk, i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "cells2edges_scalar", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int je) {
+ p_edge_out_view(je, jk, jb) =
+ coeff_int_view(je, 0, jb) *
+ p_cell_in_view(iidx_view(je, jb, 0), jk,
+ iblk_view(je, jb, 0)) +
+ coeff_int_view(je, 1, jb) *
+ p_cell_in_view(iidx_view(je, jb, 1), jk,
+ iblk_view(je, jb, 1));
+ });
+ Kokkos::fence();
+ }
+ }
+}
+
+//------------------------------------------------------------------------
+//>
+/// Computes average of scalar fields from velocity points to.
+///
+/// Computes average of scalar fields from velocity points to
+/// centers of dual faces.
+///
+template <typename T>
+void edges2verts_scalar_lib(const T *p_edge_in, const int *vert_edge_idx,
+ const int *vert_edge_blk, const T *v_int,
+ T *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v,
+ const bool lacc) {
+
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedT3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ UnmanagedConstT3D p_edge_in_view(p_edge_in, nproma, nlev, nblks_e);
+ UnmanagedConstInt3D iidx_view(vert_edge_idx, nproma, nblks_v, 6);
+ UnmanagedConstInt3D iblk_view(vert_edge_blk, nproma, nblks_v, 6);
+ UnmanagedConstT3D v_int_view(v_int, nproma, 6, nblks_v);
+ UnmanagedT3D p_vert_out_view(p_vert_out, nproma, nlev, nblks_v);
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+
+ int i_startidx, i_endidx;
+ get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, jb, i_startblk,
+ i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "edges2verts_scalar", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int jv) {
+ p_vert_out_view(jv, jk, jb) =
+ v_int_view(jv, 0, jb) * p_edge_in_view(iidx_view(jv, jb, 0), jk,
+ iblk_view(jv, jb, 0)) +
+ v_int_view(jv, 1, jb) * p_edge_in_view(iidx_view(jv, jb, 1), jk,
+ iblk_view(jv, jb, 1)) +
+ v_int_view(jv, 2, jb) * p_edge_in_view(iidx_view(jv, jb, 2), jk,
+ iblk_view(jv, jb, 2)) +
+ v_int_view(jv, 3, jb) * p_edge_in_view(iidx_view(jv, jb, 3), jk,
+ iblk_view(jv, jb, 3)) +
+ v_int_view(jv, 4, jb) * p_edge_in_view(iidx_view(jv, jb, 4), jk,
+ iblk_view(jv, jb, 4)) +
+ v_int_view(jv, 5, jb) * p_edge_in_view(iidx_view(jv, jb, 5), jk,
+ iblk_view(jv, jb, 5));
+ });
+ Kokkos::fence();
+ }
+}
+
+//------------------------------------------------------------------------
+//>
+/// Computes interpolation from edges to cells
+///
+/// Computes interpolation of scalar fields from velocity points to
+/// cell centers via given interpolation weights
+///
+template <typename T>
+void edges2cells_scalar_lib(const T *p_edge_in, const int *edge_idx,
+ const int *edge_blk, const T *coeff_int,
+ T *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c,
+ const bool lacc) {
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedT3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ // edge based scalar input field, dim: (nproma,nlev,nblks_e)
+ UnmanagedConstT3D p_edge_in_view(p_edge_in, nproma, nlev, nblks_e);
+
+ // line indices of edges of triangles, dim: (nproma,nblks_c, 3)
+ UnmanagedConstInt3D iidx_view(edge_idx, nproma, nblks_c, 3); // edge_idx_view
+
+ // block indices of edges of triangles, dim: (nproma,nblks_c, 3)
+ UnmanagedConstInt3D iblk_view(edge_blk, nproma, nblks_c, 3); // edge_blk_view
+
+ // coefficients for (area weighted) interpolation, dim:
+ // (nproma,3-cell_type,nblks_c)
+ UnmanagedConstT3D coeff_int_view(coeff_int, nproma, 3, nblks_c);
+
+ // cell based scalar output field, dim: (nproma,nlev,nblks_c)
+ UnmanagedT3D p_cell_out_view(p_cell_out, nproma, nlev, nblks_c);
+
+ int i_startidx, i_endidx;
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+ get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, jb, i_startblk,
+ i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "edges2cells_scalar_lib_inner", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int jc) {
+ p_cell_out_view(jc, jk, jb) =
+ coeff_int_view(jc, 0, jb) * p_edge_in_view(iidx_view(jc, jb, 0),
+ jk,
+ iblk_view(jc, jb, 0)) +
+ coeff_int_view(jc, 1, jb) * p_edge_in_view(iidx_view(jc, jb, 1),
+ jk,
+ iblk_view(jc, jb, 1)) +
+ coeff_int_view(jc, 2, jb) * p_edge_in_view(iidx_view(jc, jb, 2),
+ jk,
+ iblk_view(jc, jb, 2));
+ });
+ Kokkos::fence();
+ }
+}
+
+//------------------------------------------------------------------------
+//>
+/// Computes average of scalar fields from centers of cells to vertices.
+///
+template <typename T, typename S>
+void cells2verts_scalar_lib(const T *p_cell_in, const int *vert_cell_idx,
+ const int *vert_cell_blk, const S *coeff_int,
+ S *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v,
+ const bool lacc, const bool acc_async) {
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<const S ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstS3D;
+ typedef Kokkos::View<S ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedS3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ // cell based scalar input field, dim: (nproma,nlev,nblks_c)
+ UnmanagedConstT3D p_cell_in_view(p_cell_in, nproma, nlev, nblks_c);
+
+ // line indices of cells around each vertex, dim: (nproma,nblks_v, 6)
+ UnmanagedConstInt3D iidx_view(vert_cell_idx, nproma, nblks_v,
+ 6); // vert_cell_idx_view
+
+ // block indices of cells around each vertex, dim: (nproma,nblks_v, 6)
+ UnmanagedConstInt3D iblk_view(vert_cell_blk, nproma, nblks_v,
+ 6); // vert_cell_blk_view
+
+ // coefficients for interpolation, dim: (nproma,9-cell_type,nblks_v)
+ UnmanagedConstS3D coeff_int_view(coeff_int, nproma, 6, nblks_v);
+
+ // vertex based scalar output field, dim: (nproma,nlev,nblks_c)
+ UnmanagedS3D p_vert_out_view(p_vert_out, nproma, nlev, nblks_c);
+
+ int i_startidx, i_endidx;
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+
+ get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, jb, i_startblk,
+ i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "cells2verts_scalar_lib", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int jv) {
+ p_vert_out_view(jv, jk, jb) =
+ coeff_int_view(jv, 0, jb) * p_cell_in_view(iidx_view(jv, jb, 0),
+ jk,
+ iblk_view(jv, jb, 0)) +
+ coeff_int_view(jv, 1, jb) * p_cell_in_view(iidx_view(jv, jb, 1),
+ jk,
+ iblk_view(jv, jb, 1)) +
+ coeff_int_view(jv, 2, jb) * p_cell_in_view(iidx_view(jv, jb, 2),
+ jk,
+ iblk_view(jv, jb, 2)) +
+ coeff_int_view(jv, 3, jb) * p_cell_in_view(iidx_view(jv, jb, 3),
+ jk,
+ iblk_view(jv, jb, 3)) +
+ coeff_int_view(jv, 4, jb) * p_cell_in_view(iidx_view(jv, jb, 4),
+ jk,
+ iblk_view(jv, jb, 4)) +
+ coeff_int_view(jv, 5, jb) * p_cell_in_view(iidx_view(jv, jb, 5),
+ jk,
+ iblk_view(jv, jb, 5));
+ });
+ Kokkos::fence();
+ }
+}
+
+//-------------------------------------------------------------------------
+//>
+/// Same as above, but provides output optionally in single precision and
+/// assumes reversed index order of the output field in loop exchange mode
+///
+template <typename T, typename S>
+void cells2verts_scalar_ri_lib(const T *p_cell_in, const int *vert_cell_idx,
+ const int *vert_cell_blk, const T *coeff_int,
+ S *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v,
+ const bool lacc, const bool acc_async) {
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<S ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedS3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ // cell based scalar input field, dim: (nproma,nlev,nblks_c)
+ UnmanagedConstT3D p_cell_in_view(p_cell_in, nproma, nlev, nblks_c);
+
+ // line indices of cells around each vertex, dim: (nproma,nblks_v, 6)
+ UnmanagedConstInt3D iidx_view(vert_cell_idx, nproma, nblks_v,
+ 6); // vert_cell_idx_view
+
+ // block indices of cells around each vertex, dim: (nproma,nblks_v, 6)
+ UnmanagedConstInt3D iblk_view(vert_cell_blk, nproma, nblks_v,
+ 6); // vert_cell_blk_view
+
+ // coefficients for interpolation, dim: (nproma,9-cell_type,nblks_v)
+ UnmanagedConstT3D coeff_int_view(coeff_int, nproma, 6, nblks_v);
+
+ // vertex based scalar output field, dim: (nproma,nlev,nblks_c)
+#ifdef __LOOP_EXCHANGE
+ UnmanagedS3D p_vert_out_view(p_vert_out, nproma, nlev, nblks_c);
+#else
+ UnmanagedS3D p_vert_out_view(p_vert_out, nlev, nproma, nblks_c);
+#endif
+
+ int i_startidx, i_endidx;
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+
+ get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, jb, i_startblk,
+ i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "cells2verts_scalar_ri_lib", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int jv) {
+
+#ifdef __LOOP_EXCHANGE
+ p_vert_out_view(jv, jk, jb) =
+#else
+ p_vert_out_view(jk, jv, jb) =
+#endif
+ coeff_int_view(jv, 0, jb) * p_cell_in_view(iidx_view(jv, jb, 0),
+ jk,
+ iblk_view(jv, jb, 0)) +
+ coeff_int_view(jv, 1, jb) * p_cell_in_view(iidx_view(jv, jb, 1),
+ jk,
+ iblk_view(jv, jb, 1)) +
+ coeff_int_view(jv, 2, jb) * p_cell_in_view(iidx_view(jv, jb, 2),
+ jk,
+ iblk_view(jv, jb, 2)) +
+ coeff_int_view(jv, 3, jb) * p_cell_in_view(iidx_view(jv, jb, 3),
+ jk,
+ iblk_view(jv, jb, 3)) +
+ coeff_int_view(jv, 4, jb) * p_cell_in_view(iidx_view(jv, jb, 4),
+ jk,
+ iblk_view(jv, jb, 4)) +
+ coeff_int_view(jv, 5, jb) * p_cell_in_view(iidx_view(jv, jb, 5),
+ jk,
+ iblk_view(jv, jb, 5));
+ });
+ Kokkos::fence();
+ }
+}
+
+//-------------------------------------------------------------------------
+//>
+/// Computes average of scalar fields from vertices to centers of cells.
+///
+template <typename T>
+void verts2cells_scalar_lib(const T *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const T *coeff_int,
+ T *p_cell_out, const int nblks_c,
+ const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v,
+ const bool lacc) {
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedT3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ // cell based scalar input field, dim: (nproma,nlev,nblks_v)
+ UnmanagedConstT3D p_vert_in_view(p_vert_in, nproma, nlev, nblks_v);
+
+ // line indices of vertices of triangles, dim: (nproma,nblks_c, 3)
+ UnmanagedConstInt3D iidx_view(cell_index_idx, nproma, nblks_c,
+ 3); // cell_vertex_idx
+
+ // block indices of vertices of triangles, dim: (nproma,nblks_c, 3)
+ UnmanagedConstInt3D iblk_view(cell_vertex_blk, nproma, nblks_c,
+ 3); // cell_vertex_blk
+
+ // coefficients for interpolation, dim: (nproma, 3, nblks_c)
+ UnmanagedConstT3D coeff_int_view(coeff_int, nproma, 3, nblks_c);
+
+ // vertex based scalar output field, dim: (nproma,nlev,nblks_c)
+ UnmanagedT3D p_cell_out_view(p_cell_out, nproma, nlev, nblks_c);
+
+ for (int jb = 0; jb < nblks_c; ++jb) {
+
+ int nlen;
+ if (jb != nblks_c) {
+ nlen = nproma;
+ } else {
+ nlen = npromz_c;
+ }
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy({slev, 0},
+ {elev + 1, nlen});
+
+ Kokkos::parallel_for(
+ "cell_avg_lib_inner", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int jc) {
+ p_cell_out_view(jc, jk, jb) =
+ coeff_int_view(jc, 0, jb) * p_vert_in_view(iidx_view(jc, jb, 0),
+ jk,
+ iblk_view(jc, jb, 0)) +
+ coeff_int_view(jc, 1, jb) * p_vert_in_view(iidx_view(jc, jb, 1),
+ jk,
+ iblk_view(jc, jb, 1)) +
+ coeff_int_view(jc, 2, jb) * p_vert_in_view(iidx_view(jc, jb, 2),
+ jk,
+ iblk_view(jc, jb, 2));
+ });
+ Kokkos::fence();
+ }
+}
+
+//-------------------------------------------------------------------------
+//>
+/// Computes the average of a cell-based variable.
+///
+/// Computes the average of a cell-based variable
+/// over its original location and the neighboring triangles.
+/// Version with variable weighting coefficients, computed such that
+/// linear horizontal gradients are not aliased into a checkerboard noise
+/// input: lives on centers of triangles
+/// output: lives on centers of triangles
+///
+template <typename T>
+void cell_avg_lib(const T *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const T *avg_coeff,
+ T *avg_psi_c, const int i_startblk, const int i_endblk,
+ const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_c, const bool lacc) {
+ // Wrap raw pointers in unmanaged Kokkos Views.
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedT3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
+
+ // cell based variable before averaging, dim: (nproma,nlev,nblks_c)
+ UnmanagedConstT3D psi_c_view(psi_c, nproma, nlev, nblks_c);
+ // line indices of triangles next to each cell, dim: (nproma,nblks_c, 3)
+ UnmanagedConstInt3D iidx_view(cell_neighbor_idx, nproma, nblks_c,
+ 3); // cell_neighbour_idx
+ // block indices of triangles next to each cell, dim: (nproma,nblks_c, 3)
+ UnmanagedConstInt3D iblk_view(cell_neighbor_blk, nproma, nblks_c,
+ 3); // cell_neighbour_blk
+ // averaging coefficients, dim: (nproma,nlev,nblks_c)
+ UnmanagedConstT3D avg_coeff_view(avg_coeff, nproma, nlev, nblks_c);
+
+ // cell based variable after averaging, dim: (nproma,nlev,nblks_c)
+ UnmanagedT3D avg_psi_c_view(avg_psi_c, nproma, nlev, nblks_c);
+
+ int i_startidx, i_endidx;
+
+ for (int jb = i_startblk; jb < i_endblk + 1; ++jb) {
+ get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, jb, i_startblk,
+ i_endblk, i_startidx, i_endidx);
+
+ Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
+ {slev, i_startidx}, {elev + 1, i_endidx + 1});
+
+ Kokkos::parallel_for(
+ "cell_avg_lib_inner", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int jc) {
+ // calculate the weighted average
+
+ avg_psi_c_view(jc, jk, jb) =
+ psi_c_view(jc, jk, jb) * avg_coeff_view(jc, 0, jb) +
+ psi_c_view(iidx_view(jc, jb, 0), jk, iblk_view(jc, jb, 0)) *
+ avg_coeff_view(jc, 1, jb) +
+ psi_c_view(iidx_view(jc, jb, 1), jk, iblk_view(jc, jb, 1)) *
+ avg_coeff_view(jc, 2, jb) +
+ psi_c_view(iidx_view(jc, jb, 2), jk, iblk_view(jc, jb, 2)) *
+ avg_coeff_view(jc, 3, jb);
+ });
+ Kokkos::fence();
+ }
+}
+
+//-----------------------------------------------------------------------
+//
+// Explicit Instantiations
+//
+//-----------------------------------------------------------------------
+
+template void verts2edges_scalar_lib<double>(
+ const double *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const double *coeff_int, double *p_edge_out,
+ const int i_startblk, const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_v, const int nblks_e, const bool lacc);
+
+template void verts2edges_scalar_lib<float>(
+ const float *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const float *coeff_int, float *p_edge_out,
+ const int i_startblk, const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_v, const int nblks_e, const bool lacc);
+
+template void cells2edges_scalar_lib<double, double>(
+ const double *p_cell_in, const int *edge_cell_idx, const int *edge_cell_blk,
+ const double *coeff_int, double *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in, const int *i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e, const int patch_id,
+ const bool l_limited_area, const bool lfill_latbc, const bool lacc);
+
+template void cells2edges_scalar_lib<float, float>(
+ const float *p_cell_in, const int *edge_cell_idx, const int *edge_cell_blk,
+ const float *coeff_int, float *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in, const int *i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e, const int patch_id,
+ const bool l_limited_area, const bool lfill_latbc, const bool lacc);
+
+// sp2dp
+template void cells2edges_scalar_lib<float, double>(
+ const float *p_cell_in, const int *edge_cell_idx, const int *edge_cell_blk,
+ const double *coeff_int, double *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in, const int *i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e, const int patch_id,
+ const bool l_limited_area, const bool lfill_latbc, const bool lacc);
+
+template void edges2verts_scalar_lib<double>(
+ const double *p_edge_in, const int *vert_edge_idx, const int *vert_edge_blk,
+ const double *v_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v, const bool lacc);
+
+template void edges2verts_scalar_lib<float>(
+ const float *p_edge_in, const int *vert_edge_idx, const int *vert_edge_blk,
+ const float *v_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v, const bool lacc);
+
+template void edges2cells_scalar_lib<double>(
+ const double *p_edge_in, const int *edge_idx, const int *edge_blk,
+ const double *coeff_int, double *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c, const bool lacc);
+
+template void edges2cells_scalar_lib<float>(
+ const float *p_edge_in, const int *edge_idx, const int *edge_blk,
+ const float *coeff_int, float *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c, const bool lacc);
+
+template void cells2verts_scalar_lib<double, double>(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+template void cells2verts_scalar_lib<float, double>(
+ const float *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+template void cells2verts_scalar_lib<float, float>(
+ const float *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const float *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+template void cells2verts_scalar_ri_lib<double, double>(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, double *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+template void cells2verts_scalar_ri_lib<double, float>(
+ const double *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const double *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+template void cells2verts_scalar_ri_lib<float, float>(
+ const float *p_cell_in, const int *vert_cell_idx, const int *vert_cell_blk,
+ const float *coeff_int, float *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v, const bool lacc,
+ const bool acc_async);
+
+template void verts2cells_scalar_lib<double>(
+ const double *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const double *coeff_int, double *p_cell_out,
+ const int nblks_c, const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v, const bool lacc);
+
+template void verts2cells_scalar_lib<float>(
+ const float *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const float *coeff_int, float *p_cell_out,
+ const int nblks_c, const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v, const bool lacc);
+
+template void cell_avg_lib<double>(
+ const double *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const double *avg_coeff, double *avg_psi_c,
+ const int i_startblk, const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_c, const bool lacc);
+
+template void
+cell_avg_lib<float>(const float *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const float *avg_coeff,
+ float *avg_psi_c, const int i_startblk, const int i_endblk,
+ const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_c, const bool lacc);
diff --git a/src/interpolation/mo_lib_interpolation_scalar.hpp b/src/interpolation/mo_lib_interpolation_scalar.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8c8d2deaefcab69e6c6db0b3da051e0a21a070f8
--- /dev/null
+++ b/src/interpolation/mo_lib_interpolation_scalar.hpp
@@ -0,0 +1,90 @@
+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#pragma once
+
+template <typename T>
+void verts2edges_scalar_lib(const T *p_vertex_in, const int *edge_vertex_idx,
+ const int *edge_vertex_blk, const T *coeff_int,
+ T *p_edge_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_v, const int nblks_e,
+ const bool lacc);
+;
+
+template <typename T, typename S>
+void cells2edges_scalar_lib(const T *p_cell_in, const int *edge_cell_idx,
+ const int *edge_cell_blk, const S *coeff_int,
+ S *p_edge_out, const int *i_startblk_in,
+ const int *i_endblk_in, const int *i_startidx_in,
+ const int *i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblk_c, const int nblks_e,
+ const int patch_id, const bool l_limited_area,
+ const bool lfill_latbc, const bool lacc);
+
+template <typename T>
+void edges2verts_scalar_lib(const T *p_edge_in, const int *vert_edge_idx,
+ const int *vert_edge_blk, const T *v_int,
+ T *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_v,
+ const bool lacc);
+
+template <typename T>
+void edges2cells_scalar_lib(const T *p_edge_in, const int *edge_idx,
+ const int *edge_blk, const T *coeff_int,
+ T *p_cell_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_e, const int nblks_c,
+ const bool lacc);
+
+template <typename T, typename S>
+void cells2verts_scalar_lib(const T *p_cell_in, const int *vert_cell_idx,
+ const int *vert_cell_blk, const S *coeff_int,
+ S *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v,
+ const bool lacc, const bool acc_async);
+
+template <typename T, typename S>
+void cells2verts_scalar_ri_lib(const T *p_cell_in, const int *vert_cell_idx,
+ const int *vert_cell_blk, const T *coeff_int,
+ S *p_vert_out, const int i_startblk,
+ const int i_endblk, const int i_startidx_in,
+ const int i_endidx_in, const int slev,
+ const int elev, const int nproma, const int nlev,
+ const int nblks_c, const int nblks_v,
+ const bool lacc, const bool acc_async);
+
+template <typename T>
+void verts2cells_scalar_lib(const T *p_vert_in, const int *cell_index_idx,
+ const int *cell_vertex_blk, const T *coeff_int,
+ T *p_cell_out, const int nblks_c,
+ const int npromz_c, const int slev, const int elev,
+ const int nproma, const int nlev, const int nblks_v,
+ const bool lacc);
+
+template <typename T>
+void cell_avg_lib(const T *psi_c, const int *cell_neighbor_idx,
+ const int *cell_neighbor_blk, const T *avg_coeff,
+ T *avg_psi_c, const int i_startblk, const int i_endblk,
+ const int i_startidx_in, const int i_endidx_in,
+ const int slev, const int elev, const int nproma,
+ const int nlev, const int nblks_c, const bool lacc);
diff --git a/src/interpolation/mo_lib_interpolation_vector.cpp b/src/interpolation/mo_lib_interpolation_vector.cpp
index 00a914a54866703339b53f0f079556fc2a672293..8e6a28e997808f84dfa34feacf81c1ce66a3e1ac 100644
--- a/src/interpolation/mo_lib_interpolation_vector.cpp
+++ b/src/interpolation/mo_lib_interpolation_vector.cpp
@@ -9,34 +9,27 @@
// SPDX-License-Identifier: BSD-3-Clause
// ---------------------------------------------------------------
-#include "mo_lib_loopindices.hpp"
#include "mo_lib_interpolation_vector.hpp"
-// The templated C++ function using Kokkos.
-// Raw pointer arguments are wrapped into unmanaged Kokkos::Views.
-// Note: The dimensions below must match the Fortran arrays.
-// - p_vn_in and p_vt_in: dimensions [nproma, nlev, nblks_e]
-// - cell_edge_idx and cell_edge_blk: dimensions [nproma, nblks_c, 3]
-// - e_bln_c_u and e_bln_c_v: dimensions [nproma, 6, nblks_c]
-// - p_u_out and p_v_out: dimensions [nproma, nlev, nblks_c]
template <typename T>
-void edges2cells_vector_lib(
- const T* p_vn_in, const T* p_vt_in,
- const int* cell_edge_idx, const int* cell_edge_blk,
- const T* e_bln_c_u, const T* e_bln_c_v,
- T* p_u_out, T* p_v_out,
- // Additional integer parameters.
- int i_startblk, int i_endblk,
- int i_startidx_in, int i_endidx_in,
- int slev, int elev,
- int nproma,
- // Dimensions for the arrays.
- int nlev, int nblks_e, int nblks_c)
-{
+void edges2cells_vector_lib(const T *p_vn_in, const T *p_vt_in,
+ const int *cell_edge_idx, const int *cell_edge_blk,
+ const T *e_bln_c_u, const T *e_bln_c_v, T *p_u_out,
+ T *p_v_out,
+ // Additional integer parameters.
+ int i_startblk, int i_endblk, int i_startidx_in,
+ int i_endidx_in, int slev, int elev, int nproma,
+ // Dimensions for the arrays.
+ int nlev, int nblks_e, int nblks_c) {
+
// Wrap raw pointers in unmanaged Kokkos Views.
- typedef Kokkos::View<const T***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstT3D;
- typedef Kokkos::View<T***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedT3D;
- typedef Kokkos::View<const int***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstInt3D;
+ typedef Kokkos::View<const T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedConstT3D;
+ typedef Kokkos::View<T ***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged>
+ UnmanagedT3D;
+ typedef Kokkos::View<const int ***, Kokkos::LayoutLeft,
+ Kokkos::MemoryUnmanaged>
+ UnmanagedConstInt3D;
UnmanagedConstT3D p_vn_in_view(p_vn_in, nproma, nlev, nblks_e);
UnmanagedConstT3D p_vt_in_view(p_vt_in, nproma, nlev, nblks_e);
@@ -54,88 +47,76 @@ void edges2cells_vector_lib(
for (int jb = i_startblk; jb <= i_endblk; ++jb) {
// Call get_indices_c_lib to get inner loop indices for block jb.
int i_startidx, i_endidx;
- get_indices_c_lib(i_startidx_in, i_endidx_in, nproma,
- jb, i_startblk, i_endblk,
- i_startidx, i_endidx);
+ get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, jb, i_startblk,
+ i_endblk, i_startidx, i_endidx);
Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy(
{slev, i_startidx}, {elev + 1, i_endidx + 1});
- Kokkos::parallel_for("edges2cells_inner", innerPolicy,
- KOKKOS_LAMBDA(const int jk, const int jc) {
- // Compute the bilinear interpolation for cell (jc, jk, jb).
- p_u_out_view(jc, jk, jb) =
- e_bln_c_u_view(jc, 0, jb) *
- p_vn_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk, cell_edge_blk_view(jc, jb, 0) - 1) +
- e_bln_c_u_view(jc, 1, jb) *
- p_vt_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk, cell_edge_blk_view(jc, jb, 0) - 1) +
- e_bln_c_u_view(jc, 2, jb) *
- p_vn_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk, cell_edge_blk_view(jc, jb, 1) - 1) +
- e_bln_c_u_view(jc, 3, jb) *
- p_vt_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk, cell_edge_blk_view(jc, jb, 1) - 1) +
- e_bln_c_u_view(jc, 4, jb) *
- p_vn_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk, cell_edge_blk_view(jc, jb, 2) - 1) +
- e_bln_c_u_view(jc, 5, jb) *
- p_vt_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk, cell_edge_blk_view(jc, jb, 2) - 1);
+ Kokkos::parallel_for(
+ "edges2cells_inner", innerPolicy,
+ KOKKOS_LAMBDA(const int jk, const int jc) {
+ // Compute the bilinear interpolation for cell (jc, jk, jb).
+ p_u_out_view(jc, jk, jb) =
+ e_bln_c_u_view(jc, 0, jb) *
+ p_vn_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk,
+ cell_edge_blk_view(jc, jb, 0) - 1) +
+ e_bln_c_u_view(jc, 1, jb) *
+ p_vt_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk,
+ cell_edge_blk_view(jc, jb, 0) - 1) +
+ e_bln_c_u_view(jc, 2, jb) *
+ p_vn_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk,
+ cell_edge_blk_view(jc, jb, 1) - 1) +
+ e_bln_c_u_view(jc, 3, jb) *
+ p_vt_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk,
+ cell_edge_blk_view(jc, jb, 1) - 1) +
+ e_bln_c_u_view(jc, 4, jb) *
+ p_vn_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk,
+ cell_edge_blk_view(jc, jb, 2) - 1) +
+ e_bln_c_u_view(jc, 5, jb) *
+ p_vt_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk,
+ cell_edge_blk_view(jc, jb, 2) - 1);
- p_v_out_view(jc, jk, jb) =
- e_bln_c_v_view(jc, 0, jb) *
- p_vn_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk, cell_edge_blk_view(jc, jb, 0) - 1) +
- e_bln_c_v_view(jc, 1, jb) *
- p_vt_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk, cell_edge_blk_view(jc, jb, 0) - 1) +
- e_bln_c_v_view(jc, 2, jb) *
- p_vn_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk, cell_edge_blk_view(jc, jb, 1) - 1) +
- e_bln_c_v_view(jc, 3, jb) *
- p_vt_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk, cell_edge_blk_view(jc, jb, 1) - 1) +
- e_bln_c_v_view(jc, 4, jb) *
- p_vn_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk, cell_edge_blk_view(jc, jb, 2) - 1) +
- e_bln_c_v_view(jc, 5, jb) *
- p_vt_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk, cell_edge_blk_view(jc, jb, 2) - 1);
- });
+ p_v_out_view(jc, jk, jb) =
+ e_bln_c_v_view(jc, 0, jb) *
+ p_vn_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk,
+ cell_edge_blk_view(jc, jb, 0) - 1) +
+ e_bln_c_v_view(jc, 1, jb) *
+ p_vt_in_view(cell_edge_idx_view(jc, jb, 0) - 1, jk,
+ cell_edge_blk_view(jc, jb, 0) - 1) +
+ e_bln_c_v_view(jc, 2, jb) *
+ p_vn_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk,
+ cell_edge_blk_view(jc, jb, 1) - 1) +
+ e_bln_c_v_view(jc, 3, jb) *
+ p_vt_in_view(cell_edge_idx_view(jc, jb, 1) - 1, jk,
+ cell_edge_blk_view(jc, jb, 1) - 1) +
+ e_bln_c_v_view(jc, 4, jb) *
+ p_vn_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk,
+ cell_edge_blk_view(jc, jb, 2) - 1) +
+ e_bln_c_v_view(jc, 5, jb) *
+ p_vt_in_view(cell_edge_idx_view(jc, jb, 2) - 1, jk,
+ cell_edge_blk_view(jc, jb, 2) - 1);
+ });
// Optionally fence after each block if required.
Kokkos::fence();
}
}
-extern "C" void edges2cells_vector_lib_dp(
- const double* p_vn_in, const double* p_vt_in,
- const int* cell_edge_idx, const int* cell_edge_blk,
- const double* e_bln_c_u, const double* e_bln_c_v,
- double* p_u_out, double* p_v_out,
- int i_startblk, int i_endblk,
- int i_startidx_in, int i_endidx_in,
- int slev, int elev,
- int nproma,
- int nlev, int nblks_e, int nblks_c)
-{
- edges2cells_vector_lib<double>(p_vn_in, p_vt_in,
- cell_edge_idx, cell_edge_blk,
- e_bln_c_u, e_bln_c_v,
- p_u_out, p_v_out,
- i_startblk, i_endblk,
- i_startidx_in, i_endidx_in,
- slev, elev,
- nproma,
- nlev, nblks_e, nblks_c);
-}
+template void edges2cells_vector_lib<double>(
+ const double *p_vn_in, const double *p_vt_in, const int *cell_edge_idx,
+ const int *cell_edge_blk, const double *e_bln_c_u, const double *e_bln_c_v,
+ double *p_u_out, double *p_v_out,
+ // Additional integer parameters.
+ int i_startblk, int i_endblk, int i_startidx_in, int i_endidx_in, int slev,
+ int elev, int nproma,
+ // Dimensions for the arrays.
+ int nlev, int nblks_e, int nblks_c);
-extern "C" void edges2cells_vector_lib_sp(
- const float* p_vn_in, const float* p_vt_in,
- const int* cell_edge_idx, const int* cell_edge_blk,
- const float* e_bln_c_u, const float* e_bln_c_v,
- float* p_u_out, float* p_v_out,
- int i_startblk, int i_endblk,
- int i_startidx_in, int i_endidx_in,
- int slev, int elev,
- int nproma,
- int nlev, int nblks_e, int nblks_c)
-{
- edges2cells_vector_lib<float>(p_vn_in, p_vt_in,
- cell_edge_idx, cell_edge_blk,
- e_bln_c_u, e_bln_c_v,
- p_u_out, p_v_out,
- i_startblk, i_endblk,
- i_startidx_in, i_endidx_in,
- slev, elev,
- nproma,
- nlev, nblks_e, nblks_c);
-}
+template void edges2cells_vector_lib<float>(
+ const float *p_vn_in, const float *p_vt_in, const int *cell_edge_idx,
+ const int *cell_edge_blk, const float *e_bln_c_u, const float *e_bln_c_v,
+ float *p_u_out, float *p_v_out,
+ // Additional integer parameters.
+ int i_startblk, int i_endblk, int i_startidx_in, int i_endidx_in, int slev,
+ int elev, int nproma,
+ // Dimensions for the arrays.
+ int nlev, int nblks_e, int nblks_c);
\ No newline at end of file
diff --git a/src/interpolation/mo_lib_interpolation_vector.hpp b/src/interpolation/mo_lib_interpolation_vector.hpp
index 0d19b243bfec356c06c4898cf04e4b3e722dd139..91869979969f779d720c3b0e1ad58dfe2fabda36 100644
--- a/src/interpolation/mo_lib_interpolation_vector.hpp
+++ b/src/interpolation/mo_lib_interpolation_vector.hpp
@@ -8,42 +8,26 @@
// See LICENSES/ for license information
// SPDX-License-Identifier: BSD-3-Clause
// ---------------------------------------------------------------
+#pragma once
+#include "mo_lib_loopindices.hpp"
#include <Kokkos_Core.hpp>
#include <vector>
+// The templated C++ function using Kokkos.
+// Raw pointer arguments are wrapped into unmanaged Kokkos::Views.
+// Note: The dimensions below must match the Fortran arrays.
+// - p_vn_in and p_vt_in: dimensions [nproma, nlev, nblks_e]
+// - cell_edge_idx and cell_edge_blk: dimensions [nproma, nblks_c, 3]
+// - e_bln_c_u and e_bln_c_v: dimensions [nproma, 6, nblks_c]
+// - p_u_out and p_v_out: dimensions [nproma, nlev, nblks_c]
template <typename T>
-void edges2cells_vector_lib(
- const T* p_vn_in, const T* p_vt_in,
- const int* cell_edge_idx, const int* cell_edge_blk,
- const T* e_bln_c_u, const T* e_bln_c_v,
- T* p_u_out, T* p_v_out,
- // Additional integer parameters.
- int i_startblk, int i_endblk,
- int i_startidx_in, int i_endidx_in,
- int slev, int elev,
- int nproma,
- // Dimensions for the arrays.
- int nlev, int nblks_e, int nblks_c);
-
-extern "C" void edges2cells_vector_lib_dp(
- const double* p_vn_in, const double* p_vt_in,
- const int* cell_edge_idx, const int* cell_edge_blk,
- const double* e_bln_c_u, const double* e_bln_c_v,
- double* p_u_out, double* p_v_out,
- int i_startblk, int i_endblk,
- int i_startidx_in, int i_endidx_in,
- int slev, int elev,
- int nproma,
- int nlev, int nblks_e, int nblks_c);
-
-extern "C" void edges2cells_vector_lib_sp(
- const float* p_vn_in, const float* p_vt_in,
- const int* cell_edge_idx, const int* cell_edge_blk,
- const float* e_bln_c_u, const float* e_bln_c_v,
- float* p_u_out, float* p_v_out,
- int i_startblk, int i_endblk,
- int i_startidx_in, int i_endidx_in,
- int slev, int elev,
- int nproma,
- int nlev, int nblks_e, int nblks_c);
+void edges2cells_vector_lib(const T *p_vn_in, const T *p_vt_in,
+ const int *cell_edge_idx, const int *cell_edge_blk,
+ const T *e_bln_c_u, const T *e_bln_c_v, T *p_u_out,
+ T *p_v_out,
+ // Additional integer parameters.
+ int i_startblk, int i_endblk, int i_startidx_in,
+ int i_endidx_in, int slev, int elev, int nproma,
+ // Dimensions for the arrays.
+ int nlev, int nblks_e, int nblks_c);
\ No newline at end of file
diff --git a/test/c/CMakeLists.txt b/test/c/CMakeLists.txt
index 13c5dfeca1423c055a5a0f5f74cc1b2723e7c4da..c9320cb7eec1402eb107dee72f0b68d7ba3a9908 100644
--- a/test/c/CMakeLists.txt
+++ b/test/c/CMakeLists.txt
@@ -21,11 +21,16 @@ FetchContent_MakeAvailable(googletest)
# Find Kokkos (or use your existing Kokkos installation)
# find_package(Kokkos REQUIRED)
+if(IM_ENABLE_LOOP_EXCHANGE)
+ target_compile_definitions(iconmath-interpolation PRIVATE __LOOP_EXCHANGE)
+endif()
+
set(SOURCES
main.cpp
test_tdma_solver.cpp
test_interpolation_vector.cpp
test_intp_rbf.cpp
+ test_interpolation_scalar.cpp
)
# Create the test executable from your test files, including main.cpp.
add_executable(iconmath_test_c ${SOURCES})
diff --git a/test/c/test_interpolation_scalar.cpp b/test/c/test_interpolation_scalar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0ee7fa37cad59ded759911ddf10b4e986adbd35a
--- /dev/null
+++ b/test/c/test_interpolation_scalar.cpp
@@ -0,0 +1,532 @@
+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#include "mo_lib_interpolation_scalar.hpp"
+#include <Kokkos_Core.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+// Free-function helpers for 3D and 4D array sizes (assumed column-major)
+template <typename T> size_t num_elements_3d(int d1, int d2, int d3) {
+ return static_cast<size_t>(d1) * d2 * d3;
+}
+
+template <typename T> size_t num_elements_4d(int d1, int d2, int d3, int d4) {
+ return static_cast<size_t>(d1) * d2 * d3 * d4;
+}
+
+// Define a helper struct that holds the two types.
+template <typename InT, typename OutT> struct MixedPrecision {
+ using in_type = InT;
+ using out_type = OutT;
+};
+
+// Define the list of type pairs we want to test.
+typedef ::testing::Types<MixedPrecision<double, double>,
+ MixedPrecision<double, float>,
+ MixedPrecision<float, float>>
+ MixedTypes;
+
+typedef ::testing::Types<MixedPrecision<double, double>,
+ MixedPrecision<float, double>,
+ MixedPrecision<float, float>>
+ MixedTypesSP2DP;
+
+// Shared dimensions for all routines and classes
+class interp_dimensions {
+public:
+ // Constant dimensions.
+ static constexpr int nproma = 16; // inner loop length
+ static constexpr int nlev = 7; // number of vertical levels
+ static constexpr int nblks_c = 2; // number of cell blocks
+ static constexpr int nblks_e = 2; // number of edge blocks (for p_e_in)
+ static constexpr int nblks_v =
+ 2; // number of vertex blocks (for rbf arrays and outputs)
+
+ // Parameter values.
+ const int i_startblk = 0;
+ const int i_endblk = 1; // Test blocks [0, 1]
+ const int i_startidx = 2;
+ const int i_endidx = nproma - 3; // Partial range: 2 .. nproma-3
+ const int slev = 1;
+ const int elev = nlev - 1; // Partial vertical range (1 .. nlev-1)
+ const bool lacc = false; // Not using ACC-specific behavior.
+ const bool acc_async = false; // No asynchronous execution.
+};
+
+template <typename T>
+class InterpolationScalarTypedTestFixture : public ::testing::Test,
+ public interp_dimensions {
+public:
+ // Arrays used for verts2edges
+ std::vector<T> p_vertex_in; // Dimensions: (nproma, nlev, nblks_v)
+ std::vector<int> edge_vertex_idx; // Dimensions: (nproma, nblks_e, 4)
+ std::vector<int> edge_vertex_blk; // Dimensions: (nproma, nblks_e, 4)
+ std::vector<T> coeff_int_edges; // Dimensions: (nproma, 2, nblks_e)
+ std::vector<T> p_edge_out; // Dimensions: (nproma, nlev, nblks_e)
+
+ // Arrays used for edges2verts
+ std::vector<T> p_edge_in; // Dimensions: (nproma, nlev, nblks_e)
+ std::vector<int> edge_vert_idx; // Dimensions: (nproma, nblks_e, 6)
+ std::vector<int> edge_vert_blk; // Dimensions: (nproma, nblks_e, 6)
+ std::vector<T> v_int; // Dimensions: (nproma, 6, nblks_v)
+ std::vector<T> p_vert_out; // Dimensions: (nproma, nlev, nblks_v)
+
+ // Arrays used for edges2cells
+ // std::vector<T> p_edge_in; // Dimensions: (nproma, nlev, nblks_e)
+ std::vector<int> edge_idx; // Dimensions: (nproma, nblks_c, 3)
+ std::vector<int> edge_blk; // Dimensions: (nproma, nblks_c, 3)
+ std::vector<T> coeff_int_cells; // Dimensions: (nproma, 3, nblks_c)
+ std::vector<T> p_cell_out; // Dimensions: (nproma, nlev, nblks_c)
+
+ // Arrays used for verts2cells
+ std::vector<T> p_vert_in; // Dimensions: (nproma, nlev, nblks_v)
+ std::vector<int> cell_index_idx; // Dimensions: (nproma, nblks_c, 3)
+ std::vector<int> cell_index_blk; // Dimensions: (nproma, nblks_c, 3)
+
+ // Arrays used for avg_lib
+ std::vector<T> psi_c; // Dimensions: (nproma, nlev, nblks_c)
+ std::vector<int> cell_neighbor_idx; // Dimensions: (nproma, nblks_c, 3)
+ std::vector<int> cell_neighbor_blk; // Dimensions: (nproma, nblks_c, 3)
+ std::vector<T> avg_coeff; // Dimensions: (nproma, nlev, nblks_c)
+ std::vector<T> avg_psi_c; // Dimensions: (nproma, nlev, nblks_c)
+
+ const int cell_type = 6;
+ const int npromz_c = 32;
+
+ InterpolationScalarTypedTestFixture() {
+ // Allocate and initialize arrays needed for verts2edges
+ p_vertex_in.resize(num_elements_3d<T>(nproma, nlev, nblks_v),
+ static_cast<T>(1));
+ edge_vertex_idx.resize(num_elements_3d<int>(nproma, nblks_e, 4), 1);
+ edge_vertex_blk.resize(num_elements_3d<int>(nproma, nblks_e, 4), 0);
+ coeff_int_edges.resize(num_elements_3d<T>(nproma, 2, nblks_e),
+ static_cast<T>(1));
+
+ p_edge_out.resize(num_elements_3d<T>(nproma, nlev, nblks_e),
+ static_cast<T>(0));
+
+ // Allocate & Initialize arrays needed for edges2verts
+ p_edge_in.resize(num_elements_3d<T>(nproma, nlev, nblks_e),
+ static_cast<T>(1));
+ edge_vert_idx.resize(num_elements_3d<int>(nproma, nblks_e, 6), 1);
+ edge_vert_blk.resize(num_elements_3d<int>(nproma, nblks_e, 6), 0);
+ v_int.resize(num_elements_3d<T>(nproma, 6, nblks_v), static_cast<T>(1));
+
+ p_vert_out.resize(num_elements_3d<T>(nproma, nlev, nblks_v),
+ static_cast<T>(0));
+
+ // Allocate & Initialize arrays needed for edges2cells
+ edge_idx.resize(num_elements_3d<int>(nproma, nblks_c, 3), 1);
+ edge_blk.resize(num_elements_3d<int>(nproma, nblks_c, 3), 0);
+ coeff_int_cells.resize(num_elements_3d<T>(nproma, 3, nblks_c),
+ static_cast<T>(1));
+
+ p_cell_out.resize(num_elements_3d<T>(nproma, nlev, nblks_c),
+ static_cast<T>(0));
+
+ // Allocate and initialize arrays needed for verts2cells
+ p_vert_in.resize(num_elements_3d<T>(nproma, nlev, nblks_v),
+ static_cast<T>(1));
+ cell_index_idx.resize(num_elements_3d<int>(nproma, nblks_c, 3), 1);
+ cell_index_blk.resize(num_elements_3d<int>(nproma, nblks_c, 3), 0);
+
+ // Allocate and initialize arrays needed for avg_lib
+ psi_c.resize(num_elements_3d<T>(nproma, nlev, nblks_c), static_cast<T>(1));
+ cell_neighbor_idx.resize(num_elements_3d<int>(nproma, nblks_c, 3), 1);
+ cell_neighbor_blk.resize(num_elements_3d<int>(nproma, nblks_c, 3), 0);
+ avg_coeff.resize(num_elements_3d<T>(nproma, nlev, nblks_c),
+ static_cast<T>(1));
+
+ // Allocate output arrays and initialize to zero.
+ avg_psi_c.resize(num_elements_3d<T>(nproma, nlev, nblks_c),
+ static_cast<T>(0));
+ }
+};
+
+typedef ::testing::Types<float, double> SingleType;
+
+TYPED_TEST_SUITE(InterpolationScalarTypedTestFixture, SingleType);
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// ! verts2edges
+//
+////////////////////////////////////////////////////////////////////////////////
+
+TYPED_TEST(InterpolationScalarTypedTestFixture, Verts2Edges) {
+
+ verts2edges_scalar_lib<TypeParam>(
+ this->p_vertex_in.data(), this->edge_vertex_idx.data(),
+ this->edge_vertex_blk.data(), this->coeff_int_edges.data(),
+ this->p_edge_out.data(), this->i_startblk, this->i_endblk,
+ this->i_startidx, this->i_endidx, this->slev, this->elev, this->nproma,
+ this->nlev, this->nblks_v, this->nblks_e, this->lacc);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx, i_endidx] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+ // Since every contribution is 1 and there are 2 stencil points,
+ // expect 2.
+ EXPECT_NEAR(this->p_edge_out[idx], static_cast<TypeParam>(2),
+ static_cast<TypeParam>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// ! edges2verts
+//
+////////////////////////////////////////////////////////////////////////////////
+
+TYPED_TEST(InterpolationScalarTypedTestFixture, Edges2Verts) {
+
+ edges2verts_scalar_lib<TypeParam>(
+ this->p_edge_in.data(), this->edge_vert_idx.data(),
+ this->edge_vert_blk.data(), this->v_int.data(), this->p_vert_out.data(),
+ this->i_startblk, this->i_endblk, this->i_startidx, this->i_endidx,
+ this->slev, this->elev, this->nproma, this->nlev, this->nblks_e,
+ this->nblks_v, this->lacc);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx_in, i_endidx_in] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+ // Since every contribution is 1 and there are 6 stencil points,
+ // expect 6.
+ EXPECT_NEAR(this->p_vert_out[idx], static_cast<TypeParam>(6),
+ static_cast<TypeParam>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// ! edges2cells
+//
+////////////////////////////////////////////////////////////////////////////////
+
+TYPED_TEST(InterpolationScalarTypedTestFixture, Edges2Cells) {
+
+ edges2cells_scalar_lib<TypeParam>(
+ this->p_edge_in.data(), this->edge_idx.data(), this->edge_blk.data(),
+ this->coeff_int_cells.data(), this->p_cell_out.data(), this->i_startblk,
+ this->i_endblk, this->i_startidx, this->i_endidx, this->slev, this->elev,
+ this->nproma, this->nlev, this->nblks_e, this->nblks_c, this->lacc);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx_in, i_endidx_in] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+ // Since every contribution is 1 and there are 3 stencil points,
+ // expect 3.
+ EXPECT_NEAR(this->p_cell_out[idx], static_cast<TypeParam>(3),
+ static_cast<TypeParam>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
+
+TYPED_TEST(InterpolationScalarTypedTestFixture, Verts2Cells) {
+
+ verts2cells_scalar_lib<TypeParam>(
+ this->p_vert_in.data(), this->cell_index_idx.data(),
+ this->cell_index_blk.data(), this->coeff_int_cells.data(),
+ this->p_cell_out.data(), this->nblks_c, this->npromz_c, this->slev,
+ this->elev, this->nproma, this->nlev, this->nblks_v, this->lacc);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx_in, i_endidx_in] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+ // Since every contribution is 1 and there are 3 stencil points,
+ // expect 3.
+ EXPECT_NEAR(this->p_cell_out[idx], static_cast<TypeParam>(3),
+ static_cast<TypeParam>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// ! cell_avg
+//
+////////////////////////////////////////////////////////////////////////////////
+
+TYPED_TEST(InterpolationScalarTypedTestFixture, AvgLib) {
+
+ // Call the function
+ cell_avg_lib<TypeParam>(this->psi_c.data(), this->cell_neighbor_idx.data(),
+ this->cell_neighbor_blk.data(),
+ this->avg_coeff.data(), this->avg_psi_c.data(),
+ this->i_startblk, this->i_endblk, this->i_startidx,
+ this->i_endidx, this->slev, this->elev, this->nproma,
+ this->nlev, this->nblks_c, this->lacc);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx_in, i_endidx_in] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+ // Since every contribution is 1 and there are 4 stencil points,
+ // expect 4.
+ EXPECT_NEAR(this->avg_psi_c[idx], static_cast<TypeParam>(4),
+ static_cast<TypeParam>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
+
+template <typename TypePair>
+class InterpolationScalarMixedTestFixture : public ::testing::Test,
+ public interp_dimensions {
+public:
+ using InType = typename TypePair::in_type;
+ using OutType = typename TypePair::out_type;
+
+ // Arrays used for cells2edges
+ std::vector<InType> p_cell_in; // Dimensions: (nproma, nlev, nblks_c)
+ std::vector<int> edge_cell_idx; // Dimensions: (nproma, nblks_e, 2)
+ std::vector<int> edge_cell_blk; // Dimensions: (nproma, nblks_e, 2)
+ std::vector<OutType> coeff_int_edges; // Dimensions: (nproma, 2, nblks_e)
+ std::vector<OutType> p_edge_out; // Dimensions: (nproma, nlev, nblks_e)
+
+ // Further parameters for cells2edges
+ const int patch_id = 0;
+ const bool l_limited_area = false;
+ const bool lfill_latbc = false;
+ std::vector<int> i_startblk_in; // Dimensions: (2)
+ std::vector<int> i_endblk_in; // Dimensions: (2)
+ std::vector<int> i_startidx_in; // Dimensions: (2)
+ std::vector<int> i_endidx_in; // Dimensions: (2)
+
+ // Arrays used for cells2verts
+ std::vector<int> vert_cell_idx; // Dimensions: (nproma, nblks_v, 6)
+ std::vector<int> vert_cell_blk; // Dimensions: (nproma, nblks_v, 6)
+ std::vector<OutType> coeff_int_verts; // Dimensions: (nproma, 6, nblks_v)
+ std::vector<OutType> p_vert_out; // Dimensions: (nproma, nlev, nblks_v)
+
+ InterpolationScalarMixedTestFixture() {
+ // Allocate and initialize arrays needed for cells2edges
+ p_cell_in.resize(num_elements_3d<InType>(nproma, nlev, nblks_c),
+ static_cast<InType>(1));
+ edge_cell_idx.resize(num_elements_3d<int>(nproma, nblks_e, 2), 1);
+ edge_cell_blk.resize(num_elements_3d<int>(nproma, nblks_e, 2), 0);
+ coeff_int_edges.resize(num_elements_3d<InType>(nproma, 2, nblks_e),
+ static_cast<OutType>(1));
+
+ p_edge_out.resize(num_elements_3d<OutType>(nproma, nlev, nblks_e),
+ static_cast<OutType>(0));
+
+ // Allocate neighbour indexes for cells2edges
+ i_startblk_in.resize(2, i_startblk);
+ i_endblk_in.resize(2, i_endblk);
+ i_startidx_in.resize(2, i_startidx);
+ i_endidx_in.resize(2, i_endidx);
+
+ // Allocate & Initialize arrays needed for cells2verts
+ vert_cell_idx.resize(num_elements_3d<int>(nproma, nblks_v, 6), 1);
+ vert_cell_blk.resize(num_elements_3d<int>(nproma, nblks_v, 6), 0);
+ coeff_int_verts.resize(num_elements_3d<InType>(nproma, 6, nblks_v),
+ static_cast<OutType>(1));
+
+ p_vert_out.resize(num_elements_3d<OutType>(nproma, nlev, nblks_v),
+ static_cast<OutType>(0));
+ }
+};
+
+TYPED_TEST_SUITE(InterpolationScalarMixedTestFixture, MixedTypesSP2DP);
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// ! cells2edges
+//
+////////////////////////////////////////////////////////////////////////////////
+
+TYPED_TEST(InterpolationScalarMixedTestFixture, cells2edges) {
+ using InType = typename TestFixture::InType;
+ using OutType = typename TestFixture::OutType;
+
+ // Call the function
+ cells2edges_scalar_lib<InType, OutType>(
+ this->p_cell_in.data(), this->edge_cell_idx.data(),
+ this->edge_cell_blk.data(), this->coeff_int_edges.data(),
+ this->p_edge_out.data(), this->i_startblk_in.data(),
+ this->i_endblk_in.data(), this->i_startidx_in.data(),
+ this->i_endidx_in.data(), this->slev, this->elev, this->nproma,
+ this->nlev, this->nblks_c, this->nblks_e, this->patch_id,
+ this->l_limited_area, this->lfill_latbc, this->lacc);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx_in, i_endidx_in] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+ // Since every contribution is 1 and there are 2 stencil points,
+ // expect 2.
+ EXPECT_NEAR(this->p_edge_out[idx], static_cast<OutType>(2),
+ static_cast<OutType>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// ! cells2verts
+//
+////////////////////////////////////////////////////////////////////////////////
+
+TYPED_TEST(InterpolationScalarMixedTestFixture, cells2verts) {
+ using InType = typename TestFixture::InType;
+ using OutType = typename TestFixture::OutType;
+
+ cells2verts_scalar_lib<InType, OutType>(
+ this->p_cell_in.data(), this->vert_cell_idx.data(),
+ this->vert_cell_blk.data(), this->coeff_int_verts.data(),
+ this->p_vert_out.data(), this->i_startblk, this->i_endblk,
+ this->i_startidx, this->i_endidx, this->slev, this->elev, this->nproma,
+ this->nlev, this->nblks_c, this->nblks_v, this->lacc, this->acc_async);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx_in, i_endidx_in] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+ // Since every contribution is 1 and there are 6 stencil points,
+ // expect 6.
+ EXPECT_NEAR(this->p_vert_out[idx], static_cast<OutType>(6),
+ static_cast<OutType>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// ! cells2verts ri
+//
+////////////////////////////////////////////////////////////////////////////////
+
+// The test for cells2verts_ri is similar to cells2verts, but is done here
+// separtely to avoid as a differebt template instantiation is needed for the
+// function call
+template <typename Types>
+class Cells2vertsriScalarLibTestFixture : public testing::Test,
+ public interp_dimensions {
+public:
+ using InType = typename Types::in_type;
+ using OutType = typename Types::out_type;
+
+ // Arrays stored in std::vector.
+ std::vector<InType> p_cell_in; // Dimensions: (nproma, nlev, nblks_c)
+ std::vector<int> vert_cell_idx; // Dimensions: (nproma, nblks_v, 6)
+ std::vector<int> vert_cell_blk; // Dimensions: (nproma, nblks_v, 6)
+ std::vector<InType> coeff_int; // Dimensions: (nproma, 6, nblks_v)
+ std::vector<OutType> p_vert_out; // Dimensions: (nproma, nlev, nblks_v)
+
+ Cells2vertsriScalarLibTestFixture() {
+ // Allocate and initialize inputs.
+ p_cell_in.resize(num_elements_3d<InType>(nproma, nlev, nblks_c),
+ static_cast<InType>(1));
+ vert_cell_idx.resize(num_elements_3d<int>(nproma, nblks_v, 6), 1);
+ vert_cell_blk.resize(num_elements_3d<int>(nproma, nblks_v, 6), 0);
+ coeff_int.resize(num_elements_3d<InType>(nproma, 6, nblks_v),
+ static_cast<InType>(1));
+
+ // Allocate output arrays and initialize to zero.
+ p_vert_out.resize(num_elements_3d<OutType>(nproma, nlev, nblks_v),
+ static_cast<OutType>(0));
+ }
+};
+
+// Add test suite
+TYPED_TEST_SUITE(Cells2vertsriScalarLibTestFixture, MixedTypes);
+
+// Add test
+TYPED_TEST(Cells2vertsriScalarLibTestFixture, cells2verts_ri) {
+ using InType = typename TestFixture::InType;
+ using OutType = typename TestFixture::OutType;
+
+ // Call the function
+ cells2verts_scalar_ri_lib<InType, OutType>(
+ this->p_cell_in.data(), this->vert_cell_idx.data(),
+ this->vert_cell_blk.data(), this->coeff_int.data(),
+ this->p_vert_out.data(), this->i_startblk, this->i_endblk,
+ this->i_startidx, this->i_endidx, this->slev, this->elev, this->nproma,
+ this->nlev, this->nblks_c, this->nblks_v, this->lacc, this->acc_async);
+
+ // Check the outputs only for blocks in the range
+ // { [i_startblk, i_endblk], [slev,elev], [i_startidx_in, i_endidx_in] }
+ for (int block = this->i_startblk; block <= this->i_endblk; ++block) {
+ for (int level = this->slev; level < this->elev; ++level) {
+ for (int i = this->i_startidx; i < this->i_endidx; ++i) {
+ // Compute the linear index for a 3D array in column-major order:
+#ifdef __LOOP_EXCHANGE
+ size_t idx =
+ i + level * this->nproma + block * this->nproma * this->nlev;
+#else
+ size_t idx = level + i * this->nlev + block * this->nproma * this->nlev;
+#endif
+ // Since every contribution is 1 and there are 6 stencil points,
+ // expect 6.
+ EXPECT_NEAR(this->p_vert_out[idx], static_cast<OutType>(6),
+ static_cast<OutType>(1e-5))
+ << "Failure at block " << block << ", level " << level << ", index "
+ << i;
+ }
+ }
+ }
+}
diff --git a/test/c/test_interpolation_vector.cpp b/test/c/test_interpolation_vector.cpp
index 0eb5a8dc3ed9340ab34eab49da7c84d1b5d97207..680fb6e5ac669549b7e96f3fd5c94ba7a69edd3e 100644
--- a/test/c/test_interpolation_vector.cpp
+++ b/test/c/test_interpolation_vector.cpp
@@ -9,29 +9,31 @@
// SPDX-License-Identifier: BSD-3-Clause
// ---------------------------------------------------------------
-#include <gtest/gtest.h>
#include <Kokkos_Core.hpp>
+#include <gtest/gtest.h>
#include <vector>
+
#include "mo_lib_interpolation_vector.hpp"
// Dimensions for the test (small, trivial test).
-// We assume Fortran ordering: column-major, but our C wrappers will wrap raw pointers into Kokkos::Views with LayoutLeft.
+// We assume Fortran ordering: column-major, but our C wrappers will wrap raw
+// pointers into Kokkos::Views with LayoutLeft.
constexpr int nproma = 2;
-constexpr int nlev = 3;
-constexpr int nblks_e = 2; // For the edge arrays (p_vn_in, p_vt_in)
-constexpr int nblks_c = 2; // For the cell arrays and interpolation coefficients
+constexpr int nlev = 3;
+constexpr int nblks_e = 2; // For the edge arrays (p_vn_in, p_vt_in)
+constexpr int nblks_c = 2; // For the cell arrays and interpolation coefficients
// For the get_indices_c_lib inputs.
-constexpr int i_startblk = 0;
-constexpr int i_endblk = 1; // two blocks: indices 0 and 1
+constexpr int i_startblk = 0;
+constexpr int i_endblk = 1; // two blocks: indices 0 and 1
constexpr int i_startidx_in = 0;
-constexpr int i_endidx_in = nproma - 1; // 0 and 1
-constexpr int slev = 0;
-constexpr int elev = nlev - 1; // 0 .. 2
+constexpr int i_endidx_in = nproma - 1; // 0 and 1
+constexpr int slev = 0;
+constexpr int elev = nlev - 1; // 0 .. 2
-// Helper to compute total number of elements for a 3D array stored in column-major order.
-template<typename T>
-size_t num_elements(int dim1, int dim2, int dim3) {
+// Helper to compute total number of elements for a 3D array stored in
+// column-major order.
+template <typename T> size_t num_elements(int dim1, int dim2, int dim3) {
return static_cast<size_t>(dim1) * dim2 * dim3;
}
@@ -46,12 +48,13 @@ TEST(Edges2CellsTest, DPTest) {
// Here we set cell_edge_idx to 1, 2, 1 for every triple.
for (int i = 0; i < num_elements<int>(nproma, nblks_c, 3); i += 3) {
- cell_edge_idx[i] = 1;
- cell_edge_idx[i+1] = 2;
- cell_edge_idx[i+2] = 1;
+ cell_edge_idx[i] = 1;
+ cell_edge_idx[i + 1] = 2;
+ cell_edge_idx[i + 2] = 1;
}
- // Similarly, set cell_edge_blk to all ones (valid since nblks_e=2, so index 1 means block 0 after subtracting 1).
- // e_bln_c_u and e_bln_c_v: dimensions [nproma, 6, nblks_c]
+ // Similarly, set cell_edge_blk to all ones (valid since nblks_e=2, so index 1
+ // means block 0 after subtracting 1). e_bln_c_u and e_bln_c_v: dimensions
+ // [nproma, 6, nblks_c]
std::vector<double> e_bln_c_u(num_elements<double>(nproma, 6, nblks_c), 1.0);
std::vector<double> e_bln_c_v(num_elements<double>(nproma, 6, nblks_c), 1.0);
// Output arrays: dimensions [nproma, nlev, nblks_c]
@@ -62,16 +65,11 @@ TEST(Edges2CellsTest, DPTest) {
std::vector<double> p_v_ref(num_elements<double>(nproma, nlev, nblks_c), 6.0);
// Call the dp (double precision) version.
- edges2cells_vector_lib_dp(
- p_vn_in.data(), p_vt_in.data(),
- cell_edge_idx.data(), cell_edge_blk.data(),
- e_bln_c_u.data(), e_bln_c_v.data(),
- p_u_out.data(), p_v_out.data(),
- i_startblk, i_endblk,
- i_startidx_in, i_endidx_in,
- slev, elev,
- nproma,
- nlev, nblks_e, nblks_c);
+ edges2cells_vector_lib<double>(
+ p_vn_in.data(), p_vt_in.data(), cell_edge_idx.data(),
+ cell_edge_blk.data(), e_bln_c_u.data(), e_bln_c_v.data(), p_u_out.data(),
+ p_v_out.data(), i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev,
+ elev, nproma, nlev, nblks_e, nblks_c);
// Check that for each computed cell in p_u_out and p_v_out, the value is 6.
// This is because for each cell, the kernel adds 6 terms of 1*1.
@@ -90,9 +88,9 @@ TEST(Edges2CellsTest, SPTest) {
std::vector<int> cell_edge_blk(num_elements<int>(nproma, nblks_c, 3), 1);
// Set cell_edge_idx values to 1, 2, 1.
for (int i = 0; i < num_elements<int>(nproma, nblks_c, 3); i += 3) {
- cell_edge_idx[i] = 1;
- cell_edge_idx[i+1] = 2;
- cell_edge_idx[i+2] = 1;
+ cell_edge_idx[i] = 1;
+ cell_edge_idx[i + 1] = 2;
+ cell_edge_idx[i + 2] = 1;
}
std::vector<float> e_bln_c_u(num_elements<float>(nproma, 6, nblks_c), 1.0f);
std::vector<float> e_bln_c_v(num_elements<float>(nproma, 6, nblks_c), 1.0f);
@@ -103,16 +101,11 @@ TEST(Edges2CellsTest, SPTest) {
std::vector<float> p_v_ref(num_elements<float>(nproma, nlev, nblks_c), 6.0f);
// Call the sp (float precision) version.
- edges2cells_vector_lib_sp(
- p_vn_in.data(), p_vt_in.data(),
- cell_edge_idx.data(), cell_edge_blk.data(),
- e_bln_c_u.data(), e_bln_c_v.data(),
- p_u_out.data(), p_v_out.data(),
- i_startblk, i_endblk,
- i_startidx_in, i_endidx_in,
- slev, elev,
- nproma,
- nlev, nblks_e, nblks_c);
+ edges2cells_vector_lib<float>(
+ p_vn_in.data(), p_vt_in.data(), cell_edge_idx.data(),
+ cell_edge_blk.data(), e_bln_c_u.data(), e_bln_c_v.data(), p_u_out.data(),
+ p_v_out.data(), i_startblk, i_endblk, i_startidx_in, i_endidx_in, slev,
+ elev, nproma, nlev, nblks_e, nblks_c);
// Verify that every computed output equals 6.
for (size_t idx = 0; idx < p_u_out.size(); ++idx) {