src/horizontal/CMakeLists.txt

@@ -14,7 +14,8 @@ add_library(
   mo_lib_divrot.cpp
   mo_lib_divrot.F90
   mo_lib_laplace.F90
-  mo_lib_gradients.F90)
+  mo_lib_gradients.F90
+  mo_lib_gradients.cpp)
 
 add_library(${PROJECT_NAME}::horizontal ALIAS iconmath-horizontal)
 

src/horizontal/mo_lib_gradients.cpp

+// 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_loopindices.hpp"
+#include "mo_lib_gradients.hpp"
+
+// Computes directional  derivative of a cell centered variable
+// with respect to direction normal to triangle edge.
+//  input: lives on centres of triangles
+//  output:  lives on edges (velocity points)
+template <typename T>
+void grad_fd_norm_lib(const T* psi_c,
+                      const int* edge_cell_idx, const int* edge_cell_blk,
+                      const T* inv_dual_edge_length,
+                      T* grad_norm_psi_e,
+                      // Integer parameters
+                      int i_startblk, int i_endblk,
+                      int i_startidx_in, int i_endidx_in,
+                      int slev, int elev, int nproma,
+                      // Extra array dimensions
+                      int nlev, int nblks_c, int nblks_e,
+                      // OpenACC flag
+                      bool lacc)
+{
+    // Host or device check
+
+    // Wrap raw pointers in unmanaged Kokkos views
+    typedef Kokkos::View<const T***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstT3D;
+    typedef Kokkos::View<const int***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstInt3D;
+    typedef Kokkos::View<const int**, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstT2D;
+    typedef Kokkos::View<T***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedT3D;
+
+    UnmanagedConstT3D psi_c_view(psi_c, nproma, nlev, nblks_c);
+    UnmanagedConstInt3D edge_cell_idx_view(edge_cell_idx, nproma, nblks_e, 2);
+    UnmanagedConstInt3D edge_cell_blk_view(edge_cell_blk, nproma, nblks_e, 2);
+    UnmanagedT3D grad_norm_psi_e_view(grad_norm_psi_e, nproma, nlev, nblks_e);
+
+    for (int jb = i_startblk; jb < i_endblk; ++jb)
+    {
+        int i_startidx, i_endidx;
+        get_indices_e_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, i_endidx});
+        Kokkos::parallel_for("grad_fd_norm_lib_inner", innerPolicy,
+                             KOKKOS_LAMBDA(const int jk, const int je)
+        {
+            // Compute the normal derivative by the finite difference approximation
+            // See Bonaventura and Ringler MWR 2005
+            grad_norm_psi_e_view(je, jk, jb) =
+                psi_c_view(edge_cell_idx_view(je, jb, 1), jk, edge_cell_blk_view(je, jb, 1)) -
+                psi_c_view(edge_cell_idx_view(je, jb, 0), jk, edge_cell_blk_view(je, jb, 0)) *
+                inv_dual_edge_length(je, jb);
+        });
+        // Is a barrier needed here?
+        // Kokkos::fence();
+    }
+}
+

src/horizontal/mo_lib_gradients.hpp

+// 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 <Kokkos_Core.hpp>
+
+template <typename T>
+void grad_fd_norm_lib(const T* psi_c,
+                      const int* edge_cell_idx, const int* edge_cell_blk,
+                      const T* inv_dual_edge_length,
+                      T* grad_norm_psi_e,
+                      // Integer parameters
+                      int i_startblk, int i_endblk,
+                      int i_startidx_in, int i_endidx_in,
+                      int slev, int elev, int nproma,
+                      // Array dimensions
+                      int nlev, int nblks_c, int nblks_e,
+                      // OpenACC flag
+                      bool lacc);
+
+extern "C" void grad_fd_norm_lib();

src/support/mo_lib_loopindices.cpp

@@ -12,18 +12,18 @@
 #include <algorithm> // For std::max
 
 // get_indices_c_lib function
-void get_indices_c_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+void get_indices_c_lib(const int i_startidx_in, const int i_endidx_in, const int nproma,
                         const int i_blk, const int i_startblk, const int i_endblk,
                         int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true) {
-    
+
     //Since code is ported incrementally from Fortran to C++, depending on where the function is called from
     //(either fortran or c++), the first index should be either 0 or 1.
     int first_index;
     if (called_from_cpp)
         first_index = 0;
     else
-        first_index = 1;                   
-    
+        first_index = 1;
+
     if (i_blk == i_startblk) {
         i_startidx_out = std::max(first_index, i_startidx_in);
         i_endidx_out = nproma;
@@ -40,11 +40,11 @@ void get_indices_c_lib(const int i_startidx_in, const int i_endidx_in, const int
 }
 
 // get_indices_e_lib function
-void get_indices_e_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+void get_indices_e_lib(const int i_startidx_in, const int i_endidx_in, const int nproma,
                         const int i_blk, const int i_startblk, const int i_endblk,
                         int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true) {
-    
-    //Since code is ported incrementally from Fortran to C++, depending on where the function is called from, 
+
+    //Since code is ported incrementally from Fortran to C++, depending on where the function is called from,
     //the first index should be either 0 or 1.
     int first_index;
     if (called_from_cpp)
@@ -57,11 +57,11 @@ void get_indices_e_lib(const int i_startidx_in, const int i_endidx_in, const int
 }
 
 // get_indices_v_lib function
-void get_indices_v_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+void get_indices_v_lib(const int i_startidx_in, const int i_endidx_in, const int nproma,
                         const int i_blk, const int i_startblk, const int i_endblk,
                         int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true) {
-    
-    //Since code is ported incrementally from Fortran to C++, depending on where the function is called from, 
+
+    //Since code is ported incrementally from Fortran to C++, depending on where the function is called from,
     //the first index should be either 0 or 1.
     int first_index;
     if (called_from_cpp)
@@ -82,4 +82,4 @@ void get_indices_v_lib(const int i_startidx_in, const int i_endidx_in, const int
         i_startidx_out = first_index;
         i_endidx_out = nproma;
     }
-}
\ No newline at end of file
+}