From 2e58d3cbacea2657a7db8a952400df3919da6d48 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 20 May 2019 15:28:34 -0500
Subject: [PATCH] Add weighted_sum subroutine, remove remaining sum()
 intrinsics

---
 WENO.F90 | 82 ++++++++++++++++++++++++++++++++++++++------------------
 1 file changed, 56 insertions(+), 26 deletions(-)

diff --git a/WENO.F90 b/WENO.F90
index 9ae69cc..54bb893 100644
--- a/WENO.F90
+++ b/WENO.F90
@@ -394,25 +394,33 @@ subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric
   real, intent(in) :: combined_frozen_metrics
   real, intent(out) :: w(nvars, 3)
 
-  real C(3), IS(3), alpha(3), eps
-  integer p, i
-  real sum_alpha
-  integer j
+  real :: C(3) = (/0.1, 0.6, 0.3/)
+  real :: weights1(0:2) = (/1, -4, 3/)
+  real :: weights2(0:2) = (/-1, 0, 1/)
+  real :: weights3(0:2) = (/-3, 4, 1/)
+  real :: weightsc(0:2) = (/1, -2, 1/)
+
+  real IS(3), alpha(3), eps
+
+  real w1sum, w2sum, w3sum, c1sum, c2sum, c3sum
 
-  C(1) = 0.1
-  C(2) = 0.6
-  C(3) = 0.3
+  integer p, i, j
+  real sum_alpha
 
   eps = 1e-6*combined_frozen_metrics
   p = 2
 
   do i=1,nvars
-    IS(1) = (1.0/4)*(sum([1, -4, 3]*characteristic_fluxes(i,-2:0)))**2 &
-            + (13.0/12)*(sum([1, -2, 1]*characteristic_fluxes(i,-2:0)))**2
-    IS(2) = (1.0/4)*(sum([-1, 0, 1]*characteristic_fluxes(i,-1:1)))**2 &
-            + (13.0/12)*(sum([1, -2, 1]*characteristic_fluxes(i,-1:1)))**2
-    IS(3) = (1.0/4)*(sum([-3, 4, -1]*characteristic_fluxes(i,0:2)))**2 &
-            + (13.0/12)*(sum([1, -2, 1]*characteristic_fluxes(i,0:2)))**2
+    call weighted_sum(3, weights1, characteristic_fluxes(i,-2:0), w1sum)
+    call weighted_sum(3, weightsc, characteristic_fluxes(i,-2:0), c1sum)
+    call weighted_sum(3, weights2, characteristic_fluxes(i,-1:1), w2sum)
+    call weighted_sum(3, weightsc, characteristic_fluxes(i,-1:1), c2sum)
+    call weighted_sum(3, weights3, characteristic_fluxes(i,0:2), w3sum)
+    call weighted_sum(3, weightsc, characteristic_fluxes(i,0:2), c3sum)
+
+    IS(1) = (1.0/4)*w1sum**2 + (13.0/12)*c1sum**2
+    IS(2) = (1.0/4)*w2sum**2 + (13.0/12)*c2sum**2
+    IS(3) = (1.0/4)*w3sum**2 + (13.0/12)*c3sum**2
 
     do j=1,3
       alpha(j) = C(j)/(IS(j) + eps)**p
@@ -489,25 +497,32 @@ subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric
   real, intent(in) :: combined_frozen_metrics
   real, intent(out) :: w(nvars, 3)
 
-  real C(3), IS(3), alpha(3), eps
-  integer p, i
-  real sum_alpha
-  integer j
+  real :: C(3) = (/0.1, 0.6, 0.3/)
+  real :: weights1(0:2) = (/1, -4, 3/)
+  real :: weights2(0:2) = (/-1, 0, 1/)
+  real :: weights3(0:2) = (/-3, 4, 1/)
+  real :: weightsc(0:2) = (/1, -2, 1/)
 
-  C(1) = 0.1
-  C(2) = 0.6
-  C(3) = 0.3
+  real IS(3), alpha(3), eps
+  real w1sum, w2sum, w3sum, c1sum, c2sum, c3sum
+
+  integer p, i, j
+  real sum_alpha
 
   eps = 1e-6*combined_frozen_metrics
   p = 2
 
   do i=1,nvars
-    IS(1) = (1.0/4)*(sum([1, -4, 3]*characteristic_fluxes(i, 3:1:-1)))**2 &
-            + (13.0/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 3:1:-1)))**2
-    IS(2) = (1.0/4)*(sum([-1, 0, 1]*characteristic_fluxes(i, 2:0:-1)))**2 &
-            + (13.0/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 2:0:-1)))**2
-    IS(3) = (1.0/4)*(sum([-3, 4, -1]*characteristic_fluxes(i, 1:-1:-1)))**2 &
-            + (13.0/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 1:-1:-1)))**2
+    call weighted_sum(3, weights1, characteristic_fluxes(i,3:1:-1), w1sum)
+    call weighted_sum(3, weightsc, characteristic_fluxes(i,3:1:-1), c1sum)
+    call weighted_sum(3, weights2, characteristic_fluxes(i,2:0:-1), w2sum)
+    call weighted_sum(3, weightsc, characteristic_fluxes(i,2:0:-1), c2sum)
+    call weighted_sum(3, weights3, characteristic_fluxes(i,1:-1:-1), w3sum)
+    call weighted_sum(3, weightsc, characteristic_fluxes(i,1:-1:-1), c3sum)
+
+    IS(1) = (1.0/4)*w1sum**2 + (13.0/12)*c1sum**2
+    IS(2) = (1.0/4)*w2sum**2 + (13.0/12)*c2sum**2
+    IS(3) = (1.0/4)*w3sum**2 + (13.0/12)*c3sum**2
 
     do j=1,3
       alpha(j) = C(j)/(IS(j) + eps)**p
@@ -572,6 +587,21 @@ subroutine sum_array(n, a, a_sum)
   end do
 end subroutine
 
+subroutine weighted_sum(n, w, a, a_sum)
+  implicit none
+
+  integer, intent(in) :: n
+  real, intent(in) :: a(n), w(n)
+  real, intent(out) :: a_sum
+
+  integer i
+
+  a_sum = 0.0
+  do i=1,n
+    a_sum = a_sum + w(i)*a(i)
+  end do
+end subroutine
+
 !$loopy begin
 !
 ! prg = lp.parse_fortran(lp.c_preprocess(SOURCE), FILENAME)
-- 
GitLab