Draft: adding cpp version of codes

What is the new feature

Creating the C++ version of existing code

The people involved in porting the existing codes to C++ are @b382190, @b382899, @m301159, @k202204, @k202194, @k202174 and @k202170

Mandatory steps before review

• Gitlab CI passes (Hint: use make format for linting)
• New feature is covered by additional unit tests
• Mark the merge request as ready by removing Draft:

Mandatory steps before merge

• Reviewed by a maintainer
• Incorporate review suggestions
• Remember to edit the commit message and select the proper changelog category (feature/bugfix/other)
• Prior to merging, please remove any boilerplate from the MR description, retaining only the Please describe your feature in a couple of words and describe important implementation details of the feature section to maintain

You are not supposed to merge this request by yourself, the maintainers of libiconmath take care of this action!

test/fortran/test_math_utilities.f90

+! 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
+! ---------------------------------------------------------------
+
+MODULE TEST_mo_math_utilities
+
+  USE FORTUTF
+  USE mo_math_types, ONLY: t_cartesian_coordinates, t_geographical_coordinates, &
+    &                      t_line, t_tangent_vectors
+  USE mo_math_constants, ONLY: pi, pi_2, pi_4
+  USE mo_lib_grid_geometry_info
+  ! USE mo_physical_constants, ONLY: earth_radius
+  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+
+  IMPLICIT NONE
+
+  PRIVATE
+
+  PUBLIC :: TEST_cc2gc, TEST_gc2cc, TEST_cc2tv, TEST_gvec2cvec, TEST_cvec2gvec, TEST_tdma_solver_vec
+  REAL(wp), PARAMETER :: earth_radius = 6.371229e6_wp !! [m]    average radius
+
+CONTAINS
+
+  SUBROUTINE TEST_cc2gc
+
+    USE mo_math_utilities, ONLY: cc2gc
+
+    TYPE(t_cartesian_coordinates)    :: coord
+    TYPE(t_grid_geometry_info)       :: geometry_info
+    TYPE(t_geographical_coordinates) :: pos
+
+    REAL(wp) :: lon_ref, lat_ref
+
+    coord%x(1) = 10.0_wp
+    coord%x(2) = 10.0_wp
+    coord%x(3) = 5.0_wp
+
+    geometry_info%geometry_type = sphere_geometry
+
+    pos = cc2gc(coord, geometry_info)
+
+    lon_ref = 0.78539816339744828_wp
+    lat_ref = 0.33983690945412193_wp
+
+    CALL TAG_TEST("TEST_cc2gc_sphere_lon")
+    CALL ASSERT_EQUAL(pos%lon, lon_ref)
+    CALL TAG_TEST("TEST_cc2gc_sphere_lat")
+    CALL ASSERT_EQUAL(pos%lat, lat_ref)
+
+    geometry_info%geometry_type = planar_torus_geometry
+
+    geometry_info%domain_length = 2.0_wp*pi*earth_radius
+    geometry_info%domain_height = 2.0_wp*pi*earth_radius
+
+    coord%x(1) = 1000000.0_wp
+    coord%x(2) = 100.0_wp
+
+    pos = cc2gc(coord, geometry_info)
+
+    lon_ref = 0.15695558894524117_wp
+    lat_ref = -0.17453205322393880_wp
+
+    CALL TAG_TEST("TEST_cc2gc_torus_lon")
+    CALL ASSERT_EQUAL(pos%lon, lon_ref)
+    CALL TAG_TEST("TEST_cc2gc_torus_lat")
+    CALL ASSERT_EQUAL(pos%lat, lat_ref)
+
+  END SUBROUTINE TEST_cc2gc
+
+  SUBROUTINE TEST_gc2cc
+
+    USE mo_math_utilities, ONLY: gc2cc
+
+    TYPE(t_cartesian_coordinates)    :: coord, coord_ref
+    TYPE(t_grid_geometry_info)       :: geometry_info
+    TYPE(t_geographical_coordinates) :: pos
+
+    REAL(wp) :: lon_ref, lat_ref, tol
+    INTEGER :: i
+    CHARACTER(LEN=32) :: tag
+
+    tol = 1d-15
+    pos%lon = 0.78_wp
+    pos%lat = 0.34_wp
+
+    geometry_info%geometry_type = sphere_geometry
+
+    coord = gc2cc(pos, geometry_info)
+
+    coord_ref%x(1) = 0.6702170547483377_wp
+    coord_ref%x(2) = 0.6630199536212522_wp
+    coord_ref%x(3) = 0.3334870921408144_wp
+
+    DO i = 1, SIZE(coord%x)
+!     write(*,"(i4,a,f24.16)") i, ' coord%x(i): ', coord%x(i)
+      WRITE (tag, '(a,i1)') "TEST_gc2cc_sphere_", i
+      CALL TAG_TEST(tag)
+      CALL ASSERT_ALMOST_EQUAL(coord%x(i), coord_ref%x(i), tol)
+    END DO
+
+    geometry_info%geometry_type = planar_torus_geometry
+
+    geometry_info%domain_length = 2.0_wp*pi*earth_radius
+    geometry_info%domain_height = 2.0_wp*pi*earth_radius
+
+    pos%lon = 0.15_wp
+    pos%lat = -0.17_wp
+
+    coord = gc2cc(pos, geometry_info)
+
+    coord_ref%x(1) = 955684.3499999998603016_wp
+    coord_ref%x(2) = 519845.4807382422150113_wp
+    coord_ref%x(3) = 0.0_wp
+
+    DO i = 1, SIZE(coord%x)
+!     write(*,"(i4,a,f24.16)") i, ' coord%x(i): ', coord%x(i)
+      WRITE (tag, '(a,i1)') "TEST_gc2cc_torus_", i
+      CALL TAG_TEST(tag)
+      CALL ASSERT_EQUAL(coord%x(i), coord_ref%x(i))
+    END DO
+
+  END SUBROUTINE TEST_gc2cc
+
+  SUBROUTINE TEST_cc2tv
+    USE mo_math_utilities, ONLY: cc2tv
+
+    TYPE(t_cartesian_coordinates) :: coord
+    TYPE(t_geographical_coordinates) :: pos
+    TYPE(t_tangent_vectors):: tt
+    REAL(wp) :: v1_ref, v2_ref
+
+    pos%lon = pi_2
+    pos%lat = pi_4
+    coord%x(1) = 10.0_wp
+    coord%x(2) = 15.0_wp
+    coord%x(3) = 5.0_wp
+
+    tt = cc2tv(coord, pos)
+
+    v1_ref = -9.9999999999999982_wp
+    v2_ref = -7.0710678118654737_wp
+
+    CALL TAG_TEST("TEST_cc2tv_v1")
+    CALL ASSERT_EQUAL(tt%v1, v1_ref)
+    CALL TAG_TEST("TEST_cc2tv_v2")
+    CALL ASSERT_EQUAL(tt%v2, v2_ref)
+
+  END SUBROUTINE TEST_cc2tv
+
+  SUBROUTINE TEST_gvec2cvec
+    USE mo_math_utilities, ONLY: gvec2cvec
+
+    REAL(wp) :: p_gu, p_gv ! zonal and meridional vec. component
+    REAL(wp) :: p_long, p_lat ! geo. coord. of data point
+    TYPE(t_grid_geometry_info) :: geometry_info
+
+    REAL(wp) :: p_cu, p_cv, p_cw ! Cart. vector
+    REAL(wp) :: p_cu_ref, p_cv_ref, p_cw_ref ! Cart. vector ref
+
+    geometry_info%geometry_type = sphere_geometry
+
+    p_gu = 10.0_wp
+    p_gv = 5.0_wp
+
+    p_long = pi_2
+    p_lat = pi_4
+
+    CALL gvec2cvec(p_gu, p_gv, p_long, p_lat, p_cu, p_cv, p_cw, geometry_info)
+
+    p_cu_ref = -10.0_wp
+    p_cv_ref = -3.5355339059327369_wp
+    p_cw_ref = 3.5355339059327378_wp
+
+    CALL TAG_TEST("TEST_gvec2cvec_sphere_cu")
+    CALL ASSERT_EQUAL(p_cu, p_cu_ref)
+    CALL TAG_TEST("TEST_gvec2cvec_sphere_cv")
+    CALL ASSERT_EQUAL(p_cv, p_cv_ref)
+    CALL TAG_TEST("TEST_gvec2cvec_sphere_cw")
+    CALL ASSERT_EQUAL(p_cw, p_cw_ref)
+
+    geometry_info%geometry_type = planar_torus_geometry
+
+    CALL gvec2cvec(p_gu, p_gv, p_long, p_lat, p_cu, p_cv, p_cw, geometry_info)
+
+    CALL TAG_TEST("TEST_gvec2cvec_torus_cu")
+    CALL ASSERT_EQUAL(p_cu, p_gu)
+    CALL TAG_TEST("TEST_gvec2cvec_torus_cv")
+    CALL ASSERT_EQUAL(p_cv, p_gv)
+    CALL TAG_TEST("TEST_gvec2cvec_torus_cw")
+    CALL ASSERT_EQUAL(p_cw, 0.0_wp)
+
+  END SUBROUTINE TEST_gvec2cvec
+
+  SUBROUTINE TEST_cvec2gvec
+    USE mo_math_utilities, ONLY: cvec2gvec
+
+    REAL(wp) :: p_cu, p_cv, p_cw ! Cart. vector
+    REAL(wp) :: p_long, p_lat ! geo. coord. of data point
+    TYPE(t_grid_geometry_info) :: geometry_info
+
+    REAL(wp) :: p_gu, p_gv ! zonal and meridional vec. comp.
+    REAL(wp) :: p_gu_ref, p_gv_ref
+
+    geometry_info%geometry_type = sphere_geometry
+
+    p_cu = -10.0_wp
+    p_cv = -3.0_wp
+    p_cw = 3.0_wp
+
+    p_long = pi_2
+    p_lat = pi_4
+
+    CALL cvec2gvec(p_cu, p_cv, p_cw, p_long, p_lat, p_gu, p_gv, geometry_info)
+
+!   write(*,"(a,f24.16)") ' p_gu: ', p_gu
+!   write(*,"(a,f24.16)") ' p_gv: ', p_gv
+
+    p_gu_ref = 10.0_wp
+    p_gv_ref = 4.2426406871192857_wp
+
+    CALL TAG_TEST("TEST_cvec2gvec_sphere_gu")
+    CALL ASSERT_EQUAL(p_gu, p_gu_ref)
+    CALL TAG_TEST("TEST_cvec2gvec_sphere_gv")
+    CALL ASSERT_EQUAL(p_gv, p_gv_ref)
+
+    geometry_info%geometry_type = planar_torus_geometry
+
+    CALL cvec2gvec(p_cu, p_cv, p_cw, p_long, p_lat, p_gu, p_gv, geometry_info)
+
+    CALL TAG_TEST("TEST_cvec2gvec_torus_gu")
+    CALL ASSERT_EQUAL(p_gu, p_cu)
+    CALL TAG_TEST("TEST_cvec2gvec_torus_gv")
+    CALL ASSERT_EQUAL(p_gv, p_cv)
+
+  END SUBROUTINE TEST_cvec2gvec
+
+  SUBROUTINE TEST_tdma_solver_vec
+
+    USE mo_math_utilities, ONLY: tdma_solver_vec
+
+    INTEGER, PARAMETER :: n = 10
+    REAL(wp) :: a(n, n), b(n, n), c(n, n), d(n, n), x(n, n)
+    INTEGER :: i, j
+    REAL(wp) :: sum, sum_ref, tol
+    REAL(wp) :: start_time, end_time, elapsed_time
+
+    DO i = 1, n
+      DO j = 1, n
+        a(i, j) = 1.0_wp*(i + j)
+        b(i, j) = 2.0_wp*(i + j)
+        c(i, j) = 1.0_wp*(i + j)
+        d(i, j) = 1.0_wp*(i + j)
+      END DO
+    END DO
+
+    tol = 1d-15
+
+    CALL CPU_TIME(start_time)
+#ifndef __USE_CPP_BINDINGS
+    CALL tdma_solver_vec(a, b, c, d, 1, n, 1, n, x)
+#else
+    CALL tdma_solver_vec(a, b, c, d, 0, n, 0, n, n, n, x, -1)
+#endif
+    CALL CPU_TIME(end_time)
+
+    ! Compute elapsed time
+    elapsed_time = end_time - start_time
+
+    ! Output timing result
+    WRITE (*, *) "Elapsed time for tdma_solver_vec: ", elapsed_time, " seconds"
+
+    sum = 0.0_wp
+    DO i = 1, n
+      DO j = 1, n
+        sum = sum + x(i, j)
+      END DO
+    END DO
+    sum_ref = 27.2727272727272769_wp
+    CALL TAG_TEST("TEST_tdma_solver_vec_full")
+    CALL ASSERT_ALMOST_EQUAL(sum, sum_ref, tol)
+
+    x = 0.0_wp
+#ifndef __USE_CPP_BINDINGS
+    CALL tdma_solver_vec(a, b, c, d, 2, n - 1, 2, n - 1, x)
+#else
+    CALL tdma_solver_vec(a, b, c, d, 1, n - 1, 1, n - 1, n, n, x, -1)
+#endif
+    sum = 0.0_wp
+    DO i = 2, n - 1
+      DO j = 2, n - 1
+        sum = sum + x(i, j)
+      END DO
+    END DO
+    sum_ref = 17.7777777777777679_wp
+
+    CALL TAG_TEST("TEST_tdma_solver_vec_partial")
+    CALL ASSERT_ALMOST_EQUAL(sum, sum_ref, tol)
+
+  END SUBROUTINE TEST_tdma_solver_vec
+
+END MODULE TEST_mo_math_utilities