diff --git a/harmonica/_forward/point.py b/harmonica/_forward/point.py
index 57876d956..2f9bf7bca 100644
--- a/harmonica/_forward/point.py
+++ b/harmonica/_forward/point.py
@@ -9,10 +9,22 @@
 """
 
 import numpy as np
+from choclo.point import (
+    gravity_e,
+    gravity_ee,
+    gravity_en,
+    gravity_eu,
+    gravity_n,
+    gravity_nn,
+    gravity_nu,
+    gravity_pot,
+    gravity_u,
+    gravity_uu,
+)
 from numba import jit, prange
 
 from ..constants import GRAVITATIONAL_CONST
-from .utils import check_coordinate_system, distance_cartesian, distance_spherical_core
+from .utils import check_coordinate_system, distance_spherical_core
 
 
 def point_gravity(
@@ -216,9 +228,11 @@ def point_gravity(
     dispatcher(coordinate_system, parallel)(
         *coordinates, *points, masses, result, kernel
     )
-    result *= GRAVITATIONAL_CONST
+    # Invert sign of gravity_u, gravity_eu, gravity_nu
+    if field in ("g_z", "g_ez", "g_nz"):
+        result *= -1
     # Convert to more convenient units
-    if field in ("g_easting", "g_northing", "g_z"):
+    if field in ("g_e", "g_n", "g_z"):
         result *= 1e5  # SI to mGal
     tensors = ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz", "g_ne", "g_ze", "g_zn")
     if field in tensors:
@@ -249,27 +263,27 @@ def get_kernel(coordinate_system, field):
     """
     kernels = {
         "cartesian": {
-            "potential": kernel_potential_cartesian,
-            "g_z": kernel_g_z_cartesian,
-            "g_northing": kernel_g_northing_cartesian,
-            "g_easting": kernel_g_easting_cartesian,
+            "potential": gravity_pot,
+            "g_z": gravity_u,
+            "g_n": gravity_n,
+            "g_e": gravity_e,
             # diagonal tensor components
-            "g_ee": kernel_g_ee_cartesian,
-            "g_nn": kernel_g_nn_cartesian,
-            "g_zz": kernel_g_zz_cartesian,
+            "g_ee": gravity_ee,
+            "g_nn": gravity_nn,
+            "g_zz": gravity_uu,
             # non-diagonal tensor components
-            "g_en": kernel_g_en_cartesian,
-            "g_ez": kernel_g_ez_cartesian,
-            "g_nz": kernel_g_nz_cartesian,
-            "g_ne": kernel_g_en_cartesian,
-            "g_ze": kernel_g_ez_cartesian,
-            "g_zn": kernel_g_nz_cartesian,
+            "g_en": gravity_en,
+            "g_ez": gravity_eu,
+            "g_nz": gravity_nu,
+            "g_ne": gravity_en,
+            "g_ze": gravity_eu,
+            "g_zn": gravity_nu,
         },
         "spherical": {
             "potential": kernel_potential_spherical,
-            "g_z": kernel_g_z_spherical,
-            "g_northing": None,
-            "g_easting": None,
+            "g_z": kernel_gravity_u_spherical,
+            "g_n": None,
+            "g_e": None,
         },
     }
     if field not in kernels[coordinate_system]:
@@ -281,160 +295,7 @@ def get_kernel(coordinate_system, field):
 
 
 # ------------------------------------------
-# Kernel functions for Cartesian coordinates
-# ------------------------------------------
-
-
-@jit(nopython=True)
-def kernel_potential_cartesian(
-    easting, northing, upward, easting_p, northing_p, upward_p
-):
-    """
-    Kernel function for gravitational potential field in Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    return 1 / distance
-
-
-#  Acceleration components
-#  -----------------------
-
-
-@jit(nopython=True)
-def kernel_g_z_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
-    """
-    Kernel for downward component of gravitational acceleration
-
-    Use Cartesian coords.
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    # Remember that the ``g_z`` field returns the downward component of the
-    # gravitational acceleration. As a consequence, it is multiplied by -1.
-    # Notice that the ``g_z`` does not have the minus signal observed at the
-    # compoents ``g_northing`` and ``g_easting``.
-    return (upward - upward_p) / distance**3
-
-
-@jit(nopython=True)
-def kernel_g_northing_cartesian(
-    easting, northing, upward, easting_p, northing_p, upward_p
-):
-    """
-    Kernel function for northing component of gravitational acceleration
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    return -(northing - northing_p) / distance**3
-
-
-@jit(nopython=True)
-def kernel_g_easting_cartesian(
-    easting, northing, upward, easting_p, northing_p, upward_p
-):
-    """
-    Kernel function for easting component of gravitational acceleration
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    return -(easting - easting_p) / distance**3
-
-
-#  Tensor components
-#  -----------------
-
-
-@jit(nopython=True)
-def kernel_g_ee_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
-    """
-    Kernel function for g_ee component of gravitational tensor
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    return 3 * (easting - easting_p) ** 2 / distance**5 - 1 / distance**3
-
-
-@jit(nopython=True)
-def kernel_g_nn_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
-    """
-    Kernel function for g_nn component of gravitational tensor
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    return 3 * (northing - northing_p) ** 2 / distance**5 - 1 / distance**3
-
-
-@jit(nopython=True)
-def kernel_g_zz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
-    """
-    Kernel function for g_zz component of gravitational tensor
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    return 3 * (upward - upward_p) ** 2 / distance**5 - 1 / distance**3
-
-
-@jit(nopython=True)
-def kernel_g_en_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
-    """
-    Kernel function for g_en component of gravitational tensor
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    return 3 * (easting - easting_p) * (northing - northing_p) / distance**5
-
-
-@jit(nopython=True)
-def kernel_g_ez_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
-    """
-    Kernel function for g_ez component of gravitational tensor
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    # Add a minus sign to account that the z axis points downwards.
-    return -3 * (easting - easting_p) * (upward - upward_p) / distance**5
-
-
-@jit(nopython=True)
-def kernel_g_nz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
-    """
-    Kernel function for g_nz component of gravitational tensor
-
-    Use Cartesian coordinates
-    """
-    distance = distance_cartesian(
-        (easting, northing, upward), (easting_p, northing_p, upward_p)
-    )
-    # Add a minus sign to account that the z axis points downwards.
-    return -3 * (northing - northing_p) * (upward - upward_p) / distance**5
-
-
-# ------------------------------------------
-# Kernel functions for Cartesian coordinates
+# Kernel functions for Spherical coordinates
 # ------------------------------------------
 
 
@@ -448,7 +309,7 @@ def kernel_potential_spherical(
     distance, _, _ = distance_spherical_core(
         longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
     )
-    return 1 / distance
+    return 1 / distance * GRAVITATIONAL_CONST
 
 
 #  Acceleration components
@@ -456,11 +317,11 @@ def kernel_potential_spherical(
 
 
 @jit(nopython=True)
-def kernel_g_z_spherical(
+def kernel_gravity_u_spherical(
     longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
 ):
     """
-    Kernel for downward component of gravitational acceleration
+    Kernel for upward component of gravitational acceleration
 
     Use spherical coordinates
     """
@@ -468,11 +329,19 @@ def kernel_g_z_spherical(
         longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
     )
     delta_z = radius - radius_p * cospsi
-    return delta_z / distance**3
+    return -GRAVITATIONAL_CONST * delta_z / distance**3
 
 
 def point_mass_cartesian(
-    easting, northing, upward, easting_p, northing_p, upward_p, masses, out, kernel
+    easting,
+    northing,
+    upward,
+    easting_p,
+    northing_p,
+    upward_p,
+    masses,
+    out,
+    forward_func,
 ):
     """
     Compute gravitational field of point masses in Cartesian coordinates
@@ -489,19 +358,21 @@ def point_mass_cartesian(
         Array where the gravitational field on each computation point will be
         appended.
         It must have the same size of ``easting``, ``northing`` and ``upward``.
-    kernel : func
-        Kernel function that will be used to compute the gravitational field on
-        the computation points.
+    forward_func : func
+        forward_func function that will be used to compute the gravitational
+        field on the computation points. It could be one of the forward
+        modelling functions in :mod:`choclo.point`.
     """
     for l in prange(easting.size):
         for m in range(easting_p.size):
-            out[l] += masses[m] * kernel(
+            out[l] += forward_func(
                 easting[l],
                 northing[l],
                 upward[l],
                 easting_p[m],
                 northing_p[m],
                 upward_p[m],
+                masses[m],
             )
 
 
diff --git a/harmonica/_forward/tesseroid.py b/harmonica/_forward/tesseroid.py
index 4154f2e29..ddfa175bc 100644
--- a/harmonica/_forward/tesseroid.py
+++ b/harmonica/_forward/tesseroid.py
@@ -10,7 +10,6 @@
 import numpy as np
 from numba import jit, prange
 
-from ..constants import GRAVITATIONAL_CONST
 from ._tesseroid_utils import (
     _adaptive_discretization,
     _check_points_outside_tesseroids,
@@ -23,7 +22,7 @@
     density_based_discretization,
     gauss_legendre_quadrature_variable_density,
 )
-from .point import kernel_g_z_spherical, kernel_potential_spherical
+from .point import kernel_gravity_u_spherical, kernel_potential_spherical
 
 STACK_SIZE = 100
 MAX_DISCRETIZATIONS = 100000
@@ -41,7 +40,7 @@ def tesseroid_gravity(
     dtype=np.float64,
     disable_checks=False,
 ):
-    r"""
+    """
     Compute gravitational field of tesseroids on computation points.
 
     .. warning::
@@ -51,21 +50,15 @@ def tesseroid_gravity(
         It is equivalent to the opposite of the radial component, therefore
         it's positive if the acceleration vector points inside the spheroid.
 
-    .. important::
-
-        - The gravitational potential is returned in
-          :math:`\text{J}/\text{kg}`.
-        - The gravity acceleration components are returned in mgal
-          (:math:`\text{m}/\text{s}^2`).
-
     Parameters
     ----------
-    coordinates : list of arrays
-        List of arrays containing the ``longitude``, ``latitude`` and
-        ``radius`` coordinates of the computation points, defined on
-        a spherical geocentric coordinate system. Both ``longitude`` and
-        ``latitude`` should be in degrees and ``radius`` in meters.
-    tesseroids : list or 2d-array
+    coordinates : list or 1d-array
+        List or array containing ``longitude``, ``latitude`` and ``radius`` of
+        the computation points defined on a spherical geocentric coordinate
+        system.
+        Both ``longitude`` and ``latitude`` should be in degrees and ``radius``
+        in meters.
+    tesseroids : list or 1d-array
         List or array containing the coordinates of the tesseroid:
         ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric
         spherical coordinate system.
@@ -109,8 +102,7 @@ def tesseroid_gravity(
     -------
     result : array
         Gravitational field generated by the tesseroids on the computation
-        points. Gravitational potential is returned in
-        :math:`\text{J}/\text{kg}` and acceleration components in mgal.
+        points.
 
     References
     ----------
@@ -151,7 +143,10 @@ def tesseroid_gravity(
     ... )
 
     """
-    kernels = {"potential": kernel_potential_spherical, "g_z": kernel_g_z_spherical}
+    kernels = {
+        "potential": kernel_potential_spherical,
+        "g_z": kernel_gravity_u_spherical,
+    }
     if field not in kernels:
         raise ValueError("Gravitational field {} not recognized".format(field))
     # Figure out the shape and size of the output array
@@ -193,10 +188,9 @@ def tesseroid_gravity(
         kernels[field],
         dtype,
     )
-    result *= GRAVITATIONAL_CONST
-    # Convert to more convenient units
+    # Convert to more convenient units and Invert sign
     if field == "g_z":
-        result *= 1e5  # SI to mGal
+        result *= -1e5  # SI to mGal
     return result.reshape(cast.shape)
 
 
diff --git a/harmonica/tests/test_point_gravity.py b/harmonica/tests/test_point_gravity.py
index 9f710daea..e7913ecdc 100644
--- a/harmonica/tests/test_point_gravity.py
+++ b/harmonica/tests/test_point_gravity.py
@@ -14,6 +14,18 @@
 import numpy.testing as npt
 import pytest
 import verde as vd
+from choclo.point import (
+    gravity_e,
+    gravity_ee,
+    gravity_en,
+    gravity_eu,
+    gravity_n,
+    gravity_nn,
+    gravity_nu,
+    gravity_pot,
+    gravity_u,
+    gravity_uu,
+)
 
 from .._forward.point import point_gravity
 from .._forward.utils import distance_cartesian
@@ -47,7 +59,7 @@ def test_not_implemented_field():
     point_mass = [0.0, 0.0, 0.0]
     mass = 1.0
     coordinate_system = "spherical"
-    for field in ("g_northing", "g_easting"):
+    for field in ("g_n", "g_e"):
         with pytest.raises(NotImplementedError):
             point_gravity(
                 coordinates,
@@ -159,7 +171,7 @@ def test_potential_symmetry_cartesian():
 
 
 @pytest.mark.use_numba
-@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z"))
+@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z"))
 def test_acceleration_symmetry_cartesian(field):
     """
     Test if the acceleration components verify the expected symmetry
@@ -175,10 +187,10 @@ def test_acceleration_symmetry_cartesian(field):
     easting = point_mass[0] * np.ones(2)
     northing = point_mass[1] * np.ones(2)
     upward = point_mass[2] * np.ones(2)
-    if field == "g_northing":
+    if field == "g_n":
         northing[0] += distance
         northing[1] -= distance
-    elif field == "g_easting":
+    elif field == "g_e":
         easting[0] += distance
         easting[1] -= distance
     elif field == "g_z":
@@ -204,8 +216,8 @@ def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05)
     mass : float
         Mass of the point source.
     field : str
-        Acceleration component that needs to be approximated ("g_easting",
-        "g_northing", "g_z").
+        Acceleration component that needs to be approximated ("g_e",
+        "g_n", "g_z").
     delta : float
         Distance use to compute the finite difference in meters.
 
@@ -219,9 +231,9 @@ def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05)
     # Build a two computation points slightly shifted from the original
     # computation point by a small delta
     coordinates_pair = tuple([coord, coord] for coord in coordinates)
-    if field == "g_easting":
+    if field == "g_e":
         index = 0
-    elif field == "g_northing":
+    elif field == "g_n":
         index = 1
     elif field == "g_z":
         index = 2
@@ -247,7 +259,7 @@ def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05)
 
 
 @pytest.mark.use_numba
-@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z"))
+@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z"))
 @pytest.mark.parametrize(
     "coordinates, point, mass",
     (
@@ -271,7 +283,7 @@ def test_acceleration_finite_diff_cartesian(coordinates, point, mass, field):
 
 
 @pytest.mark.use_numba
-@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z"))
+@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z"))
 def test_acceleration_sign(field):
     """
     Test if acceleration components have the correct sign
@@ -281,9 +293,9 @@ def test_acceleration_sign(field):
     mass = [2670]
     # Define computation points
     coordinates = [np.zeros(3) for i in range(3)]
-    if field == "g_easting":
+    if field == "g_e":
         coordinates[0] = np.array([-150.7, -10, 79])
-    elif field == "g_northing":
+    elif field == "g_n":
         coordinates[1] = np.array([0, 100.2, 210.7])
     elif field == "g_z":
         coordinates[2] = np.array([100.11, -300.7, -400])
@@ -300,7 +312,7 @@ def test_acceleration_sign(field):
     "field",
     # fmt: off
     (
-        "potential", "g_z", "g_northing", "g_easting",
+        "potential", "g_z", "g_n", "g_e",
         "g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"
     ),
     # fmt: on
@@ -764,3 +776,86 @@ def test_point_mass_spherical_parallel():
             parallel=True,
         )
         npt.assert_allclose(result_serial, result_parallel)
+
+
+class TestAgainstChoclo:
+    """
+    Test forward modelling functions against dumb Choclo runs
+    """
+
+    @pytest.fixture()
+    def sample_points(self):
+        """
+        Return three sample prisms
+        """
+        points = np.array(
+            [
+                [-5, 5, -5, 5],
+                [-5, -5, 5, 5],
+                [-10, -10, -10, -5],
+            ],
+            dtype=float,
+        )
+        masses = np.array([100.0, -100.0, 200.0, -300.0], dtype=float)
+        return points, masses
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return four sample observation points
+        """
+        easting = np.array([-5, 10, 0, 15], dtype=float)
+        northing = np.array([14, -4, 11, 0], dtype=float)
+        upward = np.array([9, 6, 6, 12], dtype=float)
+        return (easting, northing, upward)
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize(
+        "field, choclo_func",
+        [
+            ("potential", gravity_pot),
+            ("g_e", gravity_e),
+            ("g_n", gravity_n),
+            ("g_z", gravity_u),
+            ("g_ee", gravity_ee),
+            ("g_nn", gravity_nn),
+            ("g_zz", gravity_uu),
+            ("g_en", gravity_en),
+            ("g_ez", gravity_eu),
+            ("g_nz", gravity_nu),
+        ],
+    )
+    def test_against_choclo(
+        self,
+        field,
+        choclo_func,
+        sample_coordinates,
+        sample_points,
+    ):
+        """
+        Tests forward functions against dumb runs on Choclo
+        """
+        easting, northing, upward = sample_coordinates
+        points, masses = sample_points
+        # Compute expected results with dumb choclo calls
+        expected_result = np.zeros_like(easting)
+        for i in range(easting.size):
+            for j in range(masses.size):
+                expected_result[i] += choclo_func(
+                    easting[i],
+                    northing[i],
+                    upward[i],
+                    points[j, 0],
+                    points[j, 1],
+                    points[j, 2],
+                    masses[j],
+                )
+        if field in ("g_z", "g_ez", "g_nz"):
+            expected_result *= -1  # invert sign
+        if field in ("g_e", "g_n", "g_z"):
+            expected_result *= 1e5  # convert to mGal
+        if field in ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"):
+            expected_result *= 1e9  # convert to Eotvos
+        # Compare with Harmonica results
+        result = point_gravity(sample_coordinates, points, masses, field=field)
+        npt.assert_allclose(result, expected_result)