! Copyright (C) 2019 Timothy A. Smith
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.

subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, &
    states, fluxes, metrics, metric_jacobians, flux_derivatives)
  implicit none

  integer, intent(in) :: nvars, ndim, nx, ny, nz, d
  real*8, intent(in) :: states(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
  real*8, intent(in) :: fluxes(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
  real*8, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
  real*8, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
  real*8, intent(out) :: flux_derivatives(nvars, nx, ny, nz)

  real*8 weno_flux_tmp(nvars, 0:nx+1, 0:ny+1, 0:nz+1)

  real*8 metrics_frozen(ndim, ndim)
  real*8 combined_frozen_metrics(ndim)
  real*8 generalized_states_frozen(nvars, -2:3)
  real*8 generalized_fluxes_frozen(nvars, -2:3)

  real*8 R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
  real*8 lambda_pointwise(nvars, -2:3)
  real*8 wavespeeds(nvars)

  real*8 characteristic_fluxes_pos(nvars, -2:3)
  real*8 characteristic_fluxes_neg(nvars, -2:3)

  integer i, j, k, v

  ! for the next three loops, note carefully that the inner loop indices have a different range

  !$loopy begin tagged: flux_x_compute
  do k=1,nz
    do j=1,ny
      do i=0,nx
        call convert_to_generalized_frozen(nvars, ndim, &
                                           states(:, i-2:i+3, j, k), &
                                           fluxes(:, i-2:i+3, j, k), &
                                           metrics(:, :, i:i+1, j, k), &
                                           metric_jacobians(i:i+1, j, k), &
                                           metrics_frozen, &
                                           combined_frozen_metrics, &
                                           generalized_states_frozen, &
                                           generalized_fluxes_frozen)

        call pointwise_eigenvalues(nvars, ndim, 1, states(:, i-2:i+3, j, k), &
                                    metrics(:, :, i-2:i+3, j, k), lambda_pointwise)
        call roe_eigensystem(nvars, ndim, 1, states(:, i:i+1, j, k), &
                              metrics_frozen, R, R_inv, lambda_roe)

        call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)

        call split_characteristic_fluxes(nvars, &
                                         generalized_states_frozen, &
                                         generalized_fluxes_frozen, &
                                         R_inv, &
                                         wavespeeds, &
                                         characteristic_fluxes_pos, &
                                         characteristic_fluxes_neg)

        call weno_flux(nvars, &
                        fluxes(:, i-2:i+3, j, k), &
                        characteristic_fluxes_pos, characteristic_fluxes_neg, &
                        combined_frozen_metrics(1), R, weno_flux_tmp(:, i, j, k))
      end do
    end do
  end do
  !$loopy end tagged: flux_x_compute

  !$loopy begin tagged: flux_x_diff
  do k=1,nz
    do j=1,ny
      do i=1,nx
        do v=1,nvars
          flux_derivatives(v, i, j, k) = &
            (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k))
        end do
      end do
    end do
  end do
  !$loopy end tagged: flux_x_diff

  !$loopy begin tagged: flux_y_compute
  do k=1,nz
    do i=1,nx
      do j=0,ny
        call convert_to_generalized_frozen(nvars, ndim, &
                                           states(:, i, j-2:j+3, k), &
                                           fluxes(:, i, j-2:j+3, k), &
                                           metrics(:, :, i, j:j+1, k), &
                                           metric_jacobians(i, j:j+1, k), &
                                           metrics_frozen, &
                                           combined_frozen_metrics, &
                                           generalized_states_frozen, &
                                           generalized_fluxes_frozen)

        call pointwise_eigenvalues(nvars, ndim, 2, states(:, i, j-2:j+3, k), &
                                    metrics(:, :, i, j-2:j+3, k), lambda_pointwise)
        call roe_eigensystem(nvars, ndim, 2, states(:, i, j:j+1, k), &
                              metrics_frozen, R, R_inv, lambda_roe)

        call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)

        call split_characteristic_fluxes(nvars, &
                                         generalized_states_frozen, &
                                         generalized_fluxes_frozen, &
                                         R_inv, &
                                         wavespeeds, &
                                         characteristic_fluxes_pos, &
                                         characteristic_fluxes_neg)

        call weno_flux(nvars, &
                        fluxes(:, i, j-2:j+3, k), &
                        characteristic_fluxes_pos, characteristic_fluxes_neg, &
                        combined_frozen_metrics(2), R, weno_flux_tmp(:, i, j, k))

      end do
    end do
  end do
  !$loopy end tagged: flux_y_compute

  !$loopy begin tagged: flux_y_diff
  do k=1,nz
    do j=1,ny
      do i=1,nx
        do v=1,nvars
          flux_derivatives(v, i, j, k) = &
            (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k))
        end do
      end do
    end do
  end do
  !$loopy end tagged: flux_y_diff

  !$loopy begin tagged: flux_z_compute
  do j=1,ny
    do i=1,nx
      do k=0,nz
        call convert_to_generalized_frozen(nvars, ndim, &
                                           states(:, i, j, k-2:k+3), &
                                           fluxes(:, i, j, k-2:k+3), &
                                           metrics(:, :, i, j, k:k+1), &
                                           metric_jacobians(i, j, k:k+1), &
                                           metrics_frozen, &
                                           combined_frozen_metrics, &
                                           generalized_states_frozen, &
                                           generalized_fluxes_frozen)

        call pointwise_eigenvalues(nvars, ndim, 3, states(:, i, j, k-2:k+3), &
                                    metrics(:, :, i, j, k-2:k+3), lambda_pointwise)
        call roe_eigensystem(nvars, ndim, 3, states(:, i, j, k:k+1), &
                              metrics_frozen, R, R_inv, lambda_roe)

        call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)

        call split_characteristic_fluxes(nvars, &
                                         generalized_states_frozen, &
                                         generalized_fluxes_frozen, &
                                         R_inv, &
                                         wavespeeds, &
                                         characteristic_fluxes_pos, &
                                         characteristic_fluxes_neg)

        call weno_flux(nvars, &
                        fluxes(:, i, j, k-2:k+3), &
                        characteristic_fluxes_pos, characteristic_fluxes_neg, &
                        combined_frozen_metrics(3), R, weno_flux_tmp(:, i, j, k))
      end do
    end do
  end do
  !$loopy end tagged: flux_z_compute

  !$loopy begin tagged: flux_z_diff
  do k=1,nz
    do j=1,ny
      do i=1,nx
        do v=1,nvars
          flux_derivatives(v, i, j, k) = &
            (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1))
          flux_derivatives(v, i, j, k) = 0.0d0
        end do
      end do
    end do
  end do
  !$loopy end tagged: flux_z_diff
end subroutine

subroutine pointwise_eigenvalues(nvars, ndim, d, states, metrics, lambda_pointwise)
  implicit none

  integer, intent(in) :: nvars
  integer, intent(in) :: ndim
  integer, intent(in) :: d
  real*8, intent(in) :: states(nvars, -2:3)
  real*8, intent(in) :: metrics(ndim, ndim, -2:3)
  real*8, intent(out) :: lambda_pointwise(nvars, -2:3)

  real*8 u(3), c, p, rho, ke, metric_norm
  integer v, k, i

  do k=-2,3
    rho = states(1,k)
    ke = 0.0d0
    metric_norm = 0.0d0
    do i=1,3
      u(i) = states(i+1,k)/rho
      ke = ke + u(i)**2
      metric_norm = metric_norm + metrics(d,i,k)**2
    end do
    metric_norm = sqrt(metric_norm)
    p = (1.4d0 - 1)*(states(nvars,k) - 0.5d0*rho*ke)
    c = sqrt(1.4d0*p/rho)

    do v=1,nvars-2
      lambda_pointwise(v,k) = u(d)
    end do
    lambda_pointwise(nvars-1,k) = u(d) + c*metric_norm
    lambda_pointwise(nvars,k) = u(d) - c*metric_norm
  end do
end subroutine

subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lambda_roe)
  implicit none

  integer, intent(in) :: nvars
  integer, intent(in) :: ndim
  integer, intent(in) :: d
  real*8, intent(in) :: states(nvars, 2)
  real*8, intent(in) :: metrics_frozen(ndim, ndim)
  real*8, intent(out) :: R(nvars, nvars)
  real*8, intent(out) :: R_inv(nvars, nvars)
  real*8, intent(out) :: lambda_roe(nvars)

  integer ik, il, im
  real*8 metric_norm(ndim)
  real*8 p, ke
  real*8 u_orig(ndim, 2), H_orig(2)
  real*8 sqrt_rho(2), sqrt_rho_sum
  real*8 u(ndim), rho, c, H, q
  real*8 u_tilde(ndim), k(ndim), l(ndim), m(ndim)
  integer i, j
  real*8 b1, b2, alpha, beta

  ik = d
  il = d+1
  im = d+2
  if (il > 3) il = il - 3
  if (im > 3) im = im - 3

  do j=1,2
    do i=1,ndim
      u_orig(i,j) = states(i+1,j)/states(1,j)
    end do
  end do

  do j=1,2
    ke = 0.0d0
    do i=1,ndim
      ke = ke + u_orig(i,j)**2
    end do

    p = (1.4d0 - 1)*(states(nvars,j) - 0.5d0*states(1,j)*ke)
    H_orig(j) = (states(nvars,j) + p)/states(1,j)
  end do

  do i=1,ndim
    metric_norm(i) = 0.0d0
    do j=1,ndim
      metric_norm(i) = metric_norm(i) + metrics_frozen(i,j)**2
    end do
    metric_norm(i) = sqrt(metric_norm(i))
  end do

  do i=1,ndim
    k(i) = metrics_frozen(ik,i)/metric_norm(ik)
    l(i) = metrics_frozen(il,i)/metric_norm(il)
    m(i) = metrics_frozen(im,i)/metric_norm(im)
  end do

  sqrt_rho(1) = sqrt(states(1,1))
  sqrt_rho(2) = sqrt(states(1,2))
  sqrt_rho_sum = sqrt_rho(1) + sqrt_rho(2)

  do i=1,ndim
    u(i) = (u_orig(i,1)*sqrt_rho(1) + u_orig(i,2)*sqrt_rho(2))/sqrt_rho_sum
  end do
  rho = sqrt_rho(1)*sqrt_rho(2)
  H = (H_orig(1)*sqrt_rho(1) + H_orig(2)*sqrt_rho(2))/sqrt_rho_sum

  q = 0.0d0
  do i=1,ndim
    q = q + u(i)**2
  end do

  c = sqrt((1.4d0 - 1.0d0)*(H - 0.5d0*q))

  b1 = (1.4d0 - 1.0d0)/(c**2)
  b2 = 1.0d0 + b1*q - b1*H

  u_tilde(1) = 0.0d0
  do i=1,ndim
    u_tilde(1) = u_tilde(1) + k(i)*u(i)
  end do

  u_tilde(2) = 0.0d0
  do i=1,ndim
    u_tilde(2) = u_tilde(2) + l(i)*u(i)
  end do

  u_tilde(3) = 0.0d0
  do i=1,ndim
    u_tilde(3) = u_tilde(3) + m(i)*u(i)
  end do

  alpha = rho/(2.0d0*c)
  beta = 1.0d0/(2.0d0*alpha)

  lambda_roe(1) = u_tilde(1)*metric_norm(ik)
  lambda_roe(2) = u_tilde(1)*metric_norm(ik)
  lambda_roe(3) = u_tilde(1)*metric_norm(ik)
  lambda_roe(4) = u_tilde(1)*metric_norm(ik) + c*metric_norm(ik)
  lambda_roe(5) = u_tilde(1)*metric_norm(ik) - c*metric_norm(ik)

  R(1,1) = 1.0d0
  R(2,1) = u(1)
  R(3,1) = u(2)
  R(4,1) = u(3)
  R(5,1) = H - 1.0d0/b1

  R(1,2) = 0.0d0
  R(2,2) = rho*l(1)
  R(3,2) = rho*l(2)
  R(4,2) = rho*l(3)
  R(5,2) = rho*u_tilde(2)

  R(1,3) = 0.0d0
  R(2,3) = rho*m(1)
  R(3,3) = rho*m(2)
  R(4,3) = rho*m(3)
  R(5,3) = rho*u_tilde(3)

  R(1,4) = alpha
  R(2,4) = alpha*(u(1) + c*k(1))
  R(3,4) = alpha*(u(2) + c*k(2))
  R(4,4) = alpha*(u(3) + c*k(3))
  R(5,4) = alpha*(H + c*u_tilde(1))

  R(1,5) = alpha
  R(2,5) = alpha*(u(1) - c*k(1))
  R(3,5) = alpha*(u(2) - c*k(2))
  R(4,5) = alpha*(u(3) - c*k(3))
  R(5,5) = alpha*(H - c*u_tilde(1))

  R_inv(1,1) = 1.0d0 - b2
  R_inv(1,2) = b1*u(1)
  R_inv(1,3) = b1*u(2)
  R_inv(1,4) = b1*u(3)
  R_inv(1,5) = -b1

  R_inv(2,1) = -u_tilde(2)/rho
  R_inv(2,2) = l(1)/rho
  R_inv(2,3) = l(2)/rho
  R_inv(2,4) = l(3)/rho
  R_inv(2,5) = 0.0d0

  R_inv(3,1) = -u_tilde(3)/rho
  R_inv(3,2) = m(1)/rho
  R_inv(3,3) = m(2)/rho
  R_inv(3,4) = m(3)/rho
  R_inv(3,5) = 0.0d0

  R_inv(4,1) = beta*(b2 - u_tilde(1)/c)
  R_inv(4,2) = -beta*(b1*u(1) - k(1)/c)
  R_inv(4,3) = -beta*(b1*u(2) - k(2)/c)
  R_inv(4,4) = -beta*(b1*u(3) - k(3)/c)
  R_inv(4,5) = beta*b1

  R_inv(5,1) = beta*(b2 + u_tilde(1)/c)
  R_inv(5,2) = -beta*(b1*u(1) + k(1)/c)
  R_inv(5,3) = -beta*(b1*u(2) + k(2)/c)
  R_inv(5,4) = -beta*(b1*u(3) + k(3)/c)
  R_inv(5,5) = beta*b1
end subroutine

subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: lambda_pointwise(nvars, -2:3)
  real*8, intent(in) :: lambda_roe(nvars)
  real*8, intent(out) :: wavespeeds(nvars)

  real*8 kappa
  integer v
  integer k

  kappa = 1.1d0

  do v=1,nvars
    wavespeeds(v) = abs(lambda_roe(v))
    do k=-2,3
      wavespeeds(v) = max(wavespeeds(v), abs(lambda_pointwise(v,k)))
    end do
    wavespeeds(v) = kappa*wavespeeds(v)
  end do
end subroutine

subroutine convert_to_generalized_frozen( &
            nvars, ndim, &
            states, &
            fluxes, &
            metrics, &
            metric_jacobians, &
            metrics_frozen, &
            combined_frozen_metrics, &
            generalized_states_frozen, &
            generalized_fluxes_frozen)
  integer, intent(in) :: nvars, ndim
  real*8, intent(in) :: states(nvars, -2:3)
  real*8, intent(in) :: fluxes(nvars, -2:3)
  real*8, intent(in) :: metrics(ndim, ndim, 2)
  real*8, intent(in) :: metric_jacobians(2)
  real*8, intent(out) :: metrics_frozen(ndim, ndim)
  real*8, intent(out) :: combined_frozen_metrics(ndim)
  real*8, intent(out) :: generalized_states_frozen(nvars, -2:3)
  real*8, intent(out) :: generalized_fluxes_frozen(nvars, -2:3)

  real*8 jacobian_frozen

  integer c, k, v
  integer i, j

  do j=1,ndim
    do i=1,ndim
      metrics_frozen(i,j) = 0.0d0
      do c=1,2
        metrics_frozen(i,j) = metrics_frozen(i,j) + 0.5d0*metrics(i,j,c)/metric_jacobians(c)
      end do
    end do
  end do
  jacobian_frozen = 0.5d0*(metric_jacobians(1) + metric_jacobians(2))

  do i=1,ndim
    combined_frozen_metrics(i) = 0.0d0
    do j=1,ndim
      combined_frozen_metrics(i) = combined_frozen_metrics(i) + metrics_frozen(i,j)**2
    end do
  end do

  do k=-2,3
    do v=1,nvars
      generalized_states_frozen(v,k) = states(v,k)/jacobian_frozen
    end do
  end do

  do k=-2,3
    do v=1,nvars
      ! FIXME: this is a bug
      ! ... we actually need to recompute the fluxes using frozen metrics
      ! ... but this works when the metric matrix is identity
      generalized_fluxes_frozen(v,k) = fluxes(v,k)
    end do
  end do
end subroutine

subroutine split_characteristic_fluxes(nvars, &
                                       generalized_states_frozen, &
                                       generalized_fluxes_frozen, &
                                       R_inv, &
                                       wavespeeds, &
                                       characteristic_fluxes_pos, &
                                       characteristic_fluxes_neg)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: generalized_states_frozen(nvars, -2:3)
  real*8, intent(in) :: generalized_fluxes_frozen(nvars, -2:3)
  real*8, intent(in) :: R_inv(nvars, nvars)
  real*8, intent(in) :: wavespeeds(nvars)
  real*8, intent(out) :: characteristic_fluxes_pos(nvars, -2:3)
  real*8, intent(out) :: characteristic_fluxes_neg(nvars, -2:3)

  integer k, m
  integer l

  do k=-2,3
    do m=1,nvars
      characteristic_fluxes_pos(m,k) = 0.0d0
      characteristic_fluxes_neg(m,k) = 0.0d0
      do l=1,nvars
        characteristic_fluxes_pos(m,k) = characteristic_fluxes_pos(m,k) &
          + 0.5d0*R_inv(m,l) &
          *(generalized_fluxes_frozen(l,k) + wavespeeds(m)*generalized_states_frozen(l,k))
        characteristic_fluxes_neg(m,k) = characteristic_fluxes_neg(m,k) &
          + 0.5d0*R_inv(m,l) &
          *(generalized_fluxes_frozen(l,k) - wavespeeds(m)*generalized_states_frozen(l,k))
      end do
    end do
  end do
end subroutine

subroutine weno_flux(nvars, generalized_fluxes, characteristic_fluxes_pos, &
    characteristic_fluxes_neg, combined_frozen_metric, R, flux)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: generalized_fluxes(nvars, -2:3)
  real*8, intent(in) :: characteristic_fluxes_pos(nvars, -2:3)
  real*8, intent(in) :: characteristic_fluxes_neg(nvars, -2:3)
  real*8, intent(in) :: combined_frozen_metric
  real*8, intent(in) :: R(nvars, nvars)
  real*8, intent(out) :: flux(nvars)

  real*8 consistent(nvars)
  real*8 dissipation_pos(nvars)
  real*8 dissipation_neg(nvars)
  integer v

  call consistent_part(nvars, generalized_fluxes, consistent)
  call dissipation_part_pos(nvars, characteristic_fluxes_pos, &
    combined_frozen_metric, R, dissipation_pos)
  call dissipation_part_neg(nvars, characteristic_fluxes_neg, &
    combined_frozen_metric, R, dissipation_neg)

  do v=1,nvars
    flux(v) = consistent(v) + dissipation_pos(v) + dissipation_neg(v)
  end do
end subroutine

subroutine consistent_part(nvars, generalized_fluxes, consistent)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: generalized_fluxes(nvars, -2:3)
  real*8, intent(out) :: consistent(nvars)

  integer v

  do v=1,nvars
    consistent(v) = generalized_fluxes(v,-2) - 8*generalized_fluxes(v,-1) &
      + 37*generalized_fluxes(v,0) + 37*generalized_fluxes(v,1) &
      - 8*generalized_fluxes(v,2) + generalized_fluxes(v,3)
    consistent(v) = consistent(v)/60.0d0
  end do
end subroutine

subroutine dissipation_part_pos(nvars, characteristic_fluxes, &
  combined_frozen_metric, R, dissipation_pos)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(in) :: combined_frozen_metric
  real*8, intent(in) :: R(nvars, nvars)
  real*8, intent(out) :: dissipation_pos(nvars)

  real*8 combined_fluxes(nvars)

  call weno_combination_pos(nvars, characteristic_fluxes, combined_frozen_metric, combined_fluxes)

  call mult_mat_vec(nvars, nvars, -1.0d0/60, R, combined_fluxes, dissipation_pos)
end subroutine

subroutine weno_combination_pos(nvars, characteristic_fluxes, &
  combined_frozen_metric, combined_fluxes)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(in) :: combined_frozen_metric
  real*8, intent(out) :: combined_fluxes(nvars)

  real*8 w(nvars, 3)
  real*8 flux_differences(nvars, 3)
  integer v

  call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric, w)
  call flux_differences_pos(nvars, characteristic_fluxes, flux_differences)

  do v=1,nvars
    combined_fluxes(v) = (20*w(v,1) - 1)*flux_differences(v,1) &
      - (10*(w(v,1) + w(v,2)) - 5)*flux_differences(v,2) &
      + flux_differences(v,3)
  end do
end subroutine

subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric, w)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(in) :: combined_frozen_metric
  real*8, intent(out) :: w(nvars, 3)

  real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
  real*8 IS(nvars, 3), alpha(3), eps

  integer p, i, j
  real*8 sum_alpha(1)

  eps = 1.0d-6*combined_frozen_metric
  p = 2

  call oscillation_pos(nvars, characteristic_fluxes, IS)

  do i=1,nvars
    sum_alpha(1) = 0.0d0
    do j=1,3
      alpha(j) = C(j)/(IS(i, j) + eps)**p
      sum_alpha(1) = sum_alpha(1) + alpha(j)
    end do
    !call sum_array(3, alpha, sum_alpha)
    do j=1,3
      w(i,j) = alpha(j)/sum_alpha(1)
    end do
  end do
end subroutine

subroutine oscillation_pos(nvars, characteristic_fluxes, oscillation)
  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(out) :: oscillation(nvars, 3)

  real*8 :: weights1(0:2) = (/ 1.0d0, -4.0d0,  3.0d0/)
  real*8 :: weights2(0:2) = (/-1.0d0,  0.0d0,  1.0d0/)
  real*8 :: weights3(0:2) = (/-3.0d0,  4.0d0, -1.0d0/)
  real*8 :: weightsc(0:2) = (/ 1.0d0, -2.0d0,  1.0d0/)

  real*8 w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)

  integer i

  do i=1,nvars
    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)

    oscillation(i, 1) = (1.0d0/4)*w1sum(1)**2 + (13.0d0/12)*c1sum(1)**2
    oscillation(i, 2) = (1.0d0/4)*w2sum(1)**2 + (13.0d0/12)*c2sum(1)**2
    oscillation(i, 3) = (1.0d0/4)*w3sum(1)**2 + (13.0d0/12)*c3sum(1)**2
  end do
end subroutine

subroutine flux_differences_pos(nvars, characteristic_fluxes, flux_differences)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(out) :: flux_differences(nvars, 3)

  integer i, v

  do i=1,3
    do v=1,nvars
      flux_differences(v,i) = -characteristic_fluxes(v,-3+i) + 3*characteristic_fluxes(v,-2+i) &
                                - 3*characteristic_fluxes(v,-1+i) + characteristic_fluxes(v,i)
    end do
  end do
end subroutine

subroutine dissipation_part_neg(nvars, characteristic_fluxes, &
  combined_frozen_metric, R, dissipation_neg)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(in) :: combined_frozen_metric
  real*8, intent(in) :: R(nvars, nvars)
  real*8, intent(out) :: dissipation_neg(nvars)

  real*8 combined_fluxes(nvars)

  call weno_combination_neg(nvars, characteristic_fluxes, combined_frozen_metric, combined_fluxes)

  call mult_mat_vec(nvars, nvars, 1.0d0/60, R, combined_fluxes, dissipation_neg)
end subroutine

subroutine weno_combination_neg(nvars, characteristic_fluxes, &
  combined_frozen_metric, combined_fluxes)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(in) :: combined_frozen_metric
  real*8, intent(out) :: combined_fluxes(nvars)

  real*8 w(nvars, 3)
  real*8 flux_differences(nvars, 3)
  integer v

  call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric, w)
  call flux_differences_neg(nvars, characteristic_fluxes, flux_differences)

  do v=1,nvars
    combined_fluxes(v) = (20*w(v,1) - 1)*flux_differences(v,1) &
      - (10*(w(v,1) + w(v,2)) - 5)*flux_differences(v,2) &
      + flux_differences(v,3)
  end do
end subroutine

subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric, w)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(in) :: combined_frozen_metric
  real*8, intent(out) :: w(nvars, 3)

  real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
  real*8 IS(nvars, 3), alpha(3), eps

  integer p, i, j
  real*8 sum_alpha(1)

  eps = 1.0d-6*combined_frozen_metric
  p = 2

  call oscillation_neg(nvars, characteristic_fluxes, IS)

  do i=1,nvars
    do j=1,3
      alpha(j) = C(j)/(IS(i, j) + eps)**p
    end do
    call sum_array(3, alpha, sum_alpha)
    do j=1,3
      w(i,j) = alpha(j)/sum_alpha(1)
    end do
  end do
end subroutine

subroutine oscillation_neg(nvars, characteristic_fluxes, oscillation)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(out) :: oscillation(nvars, 3)

  real*8 :: weights1(0:2) = (/ 1.0d0, -4.0d0,  3.0d0/)
  real*8 :: weights2(0:2) = (/-1.0d0,  0.0d0,  1.0d0/)
  real*8 :: weights3(0:2) = (/-3.0d0,  4.0d0, -1.0d0/)
  real*8 :: weightsc(0:2) = (/ 1.0d0, -2.0d0,  1.0d0/)

  real*8 w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)

  integer i

  do i=1,nvars
    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)

    oscillation(i, 1) = (1.0d0/4)*w1sum(1)**2 + (13.0d0/12)*c1sum(1)**2
    oscillation(i, 2) = (1.0d0/4)*w2sum(1)**2 + (13.0d0/12)*c2sum(1)**2
    oscillation(i, 3) = (1.0d0/4)*w3sum(1)**2 + (13.0d0/12)*c3sum(1)**2
  end do
end subroutine

subroutine flux_differences_neg(nvars, characteristic_fluxes, flux_differences)
  implicit none

  integer, intent(in) :: nvars
  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real*8, intent(out) :: flux_differences(nvars, 3)

  integer i, v

  do i=1,3
    do v=1,nvars
      flux_differences(v,i) = -characteristic_fluxes(v,1-i) + 3*characteristic_fluxes(v,2-i) &
                                - 3*characteristic_fluxes(v,3-i) + characteristic_fluxes(v,4-i)
    end do
  end do
end subroutine

subroutine mult_mat_vec(m, n, alpha, a, b, c)
  implicit none

  integer, intent(in) :: m, n
  real*8, intent(in) :: alpha
  real*8, intent(in) :: a(m, n)
  real*8, intent(in) :: b(n)
  real*8, intent(out) :: c(m)

  integer i, j
  real*8 accumulator

  do i=1,m
    accumulator = 0.0d0
    do j=1,n
      accumulator = accumulator + alpha*a(i,j)*b(j)
    end do
    c(i) = accumulator
  end do
end subroutine

subroutine sum_array(n, a, a_sum)
  implicit none

  integer, intent(in) :: n
  real*8, intent(in) :: a(n)
  real*8, intent(out) :: a_sum(1)

  integer i

  a_sum(1) = 0.0d0
  do i=1,n
    a_sum(1) = a_sum(1) + a(i)
  end do
end subroutine

subroutine weighted_sum(n, w, a, a_sum)
  implicit none

  integer, intent(in) :: n
  real*8, intent(in) :: a(n), w(n)
  real*8, intent(out) :: a_sum(1)

  integer i

  a_sum(1) = 0.0d0
  do i=1,n
    a_sum(1) = a_sum(1) + w(i)*a(i)
  end do
end subroutine

!$loopy begin
!
! prg = lp.parse_fortran(lp.c_preprocess(SOURCE), FILENAME)
! prg = lp.fix_parameters(prg, ndim=3, nvars=5, _remove=False)
!
! knl = prg["compute_flux_derivatives"]
!
! knl = lp.assume(knl, "nx > 0 and ny > 0 and nz > 0")
!
! knl = lp.set_temporary_scope(knl, "weno_flux_tmp",
!         lp.AddressSpace.GLOBAL)
!
! prg = prg.with_kernel(knl)
!
! RESULT = prg
!
!$loopy end