From 9256b58905f1148e644ff8d49aeb3ea22b16ea7d Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Mon, 6 Mar 2017 11:12:25 -0300
Subject: [PATCH] SOURCES: Plug the GRAVITY module in.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
 src/makefile    |  4 ++--
 src/sources.F90 | 39 ++++++++++++---------------------------
 2 files changed, 14 insertions(+), 29 deletions(-)

diff --git a/src/makefile b/src/makefile
index 708f2c2..0734730 100644
--- a/src/makefile
+++ b/src/makefile
@@ -159,8 +159,8 @@ schemes.o        : schemes.F90 algebra.o coordinates.o equations.o             \
                    interpolations.o timers.o
 shapes.o         : shapes.F90 blocks.o constants.o coordinates.o equations.o   \
                    parameters.o timers.o
-sources.o        : sources.F90 blocks.o coordinates.o equations.o operators.o  \
-                   parameters.o timers.o
+sources.o        : sources.F90 blocks.o coordinates.o equations.o gravity.o    \
+                   operators.o parameters.o timers.o
 timers.o         : timers.F90
 utils.o          : utils.F90 error.o
 
diff --git a/src/sources.F90 b/src/sources.F90
index c04da1f..9c4b64a 100644
--- a/src/sources.F90
+++ b/src/sources.F90
@@ -49,10 +49,6 @@ module sources
 !
   integer     , save :: glm_type    = 0
 
-! gravitational acceleration coefficient
-!
-  real(kind=8), save :: gpoint      = 0.0d+00
-
 ! viscosity coefficient
 !
   real(kind=8), save :: viscosity   = 0.0d+00
@@ -144,10 +140,6 @@ module sources
       glm_type = 0
     end select
 
-! get acceleration coefficient
-!
-      call get_parameter_real("gpoint"     , gpoint   )
-
 ! get viscosity coefficient
 !
       call get_parameter_real("viscosity"  , viscosity)
@@ -161,7 +153,6 @@ module sources
     if (verbose) then
 
       write (*,"(4x,a,1x,a)    ") "glm source terms       =", trim(tglm)
-      write (*,"(4x,a,1x,1e9.2)") "point mass constant    =", gpoint
       write (*,"(4x,a,1x,1e9.2)") "viscosity              =", viscosity
       write (*,"(4x,a,1x,1e9.2)") "resistivity            =", resistivity
 
@@ -242,6 +233,7 @@ module sources
     use equations      , only : nv, inx, iny, inz
     use equations      , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
     use equations      , only : ibx, iby, ibz, ibp
+    use gravity        , only : gravity_enabled, gravitational_acceleration
     use operators      , only : divergence, gradient, laplace, curl
 
 ! local variables are not implicit by default
@@ -257,7 +249,7 @@ module sources
 !
     integer       :: i  , j  , k
     real(kind=8)  :: fc, gc
-    real(kind=8)  :: r2, r3, gx, gy, gz
+    real(kind=8)  :: gx, gy, gz
     real(kind=8)  :: dbx, dby, dbz
     real(kind=8)  :: dvxdx, dvxdy, dvxdz, divv
     real(kind=8)  :: dvydx, dvydy, dvydz
@@ -265,7 +257,7 @@ module sources
 
 ! local arrays
 !
-    real(kind=8), dimension(3)  :: dh
+    real(kind=8), dimension(3)  :: ga, dh
     real(kind=8), dimension(im) :: x
     real(kind=8), dimension(jm) :: y
     real(kind=8), dimension(km) :: z
@@ -280,9 +272,9 @@ module sources
     call start_timer(imu)
 #endif /* PROFILE */
 
-! proceed only if the gravitational acceleration coefficient is not zero
+! proceed only if the gravitational term is enabled
 !
-    if (gpoint /= 0.0d+00) then
+    if (gravity_enabled) then
 
 ! prepare block coordinates
 !
@@ -298,23 +290,16 @@ module sources
         do j = 1, jm
           do i = 1, im
 
-! calculate distance from the origin
+! get gravitational acceleration components
 !
-#if NDIMS == 2
-            r2 = x(i) * x(i) + y(j) * y(j)
-#endif /* NDIMS == 2 */
-#if NDIMS == 3
-            r2 = x(i) * x(i) + y(j) * y(j) + z(k) * z(k)
-#endif /* NDIMS == 3 */
-            r3 = r2 * sqrt(r2)
+            call gravitational_acceleration(x(i), y(j), z(k), ga(1:3))
 
-! calculate gravitational acceleration factors
+! calculate the gravitational source terms
 !
-            gc = gpoint * pdata%q(idn,i,j,k) / max(1.0d-16, r3)
-            gx = gc * x(i)
-            gy = gc * y(j)
+            gx = pdata%q(idn,i,j,k) * ga(1)
+            gy = pdata%q(idn,i,j,k) * ga(2)
 #if NDIMS == 3
-            gz = gc * z(k)
+            gz = pdata%q(idn,i,j,k) * ga(3)
 #endif /* NDIMS == 3 */
 
 ! add source terms to momentum equations
@@ -345,7 +330,7 @@ module sources
         end do ! j = 1, jm
       end do ! k = 1, km
 
-    end if ! gpoint is not zero
+    end if ! gravity enabled
 
 ! proceed only if the viscosity coefficient is not zero
 !