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

  integer, intent(in) :: nvars, ndim, nx, ny, nz
  real, intent(in) :: states(nvars, nx, ny, nz)
  real, intent(in) :: fluxes(nvars, ndim, nx, ny, nz)
  real, intent(in) :: metrics(ndim, ndim, nx, ny, nz)
  real, intent(in) :: metric_jacobians(nx, ny, nz)
  real, intent(out) :: flux_derivatives(nvars, ndim, nx, ny, nz)

  real flux_derivatives_generalized(nvars, ndim, nx, ny, nz)
  real generalized_fluxes(nvars, ndim, nx, ny, nz)
  real generalized_states_frozen(nvars, -2:3)
  real generalized_fluxes_frozen(nvars, ndim, -2:3)
  real R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
  real lambda_pointwise(nvars, -2:3)
  real lambda(nvars)
  integer i, j, k
  real delta_xi, delta_eta, delta_zeta
  real flux(nvars)
  real characteristic_fluxes_pos(nvars, -2:3)
  real characteristic_fluxes_neg(nvars, -2:3)
  real combined_frozen_metrics(ndim)
  integer v

  ! grid spacing in generalized coordinates
  delta_xi = 1.0
  delta_eta = 1.0
  delta_zeta = 1.0

  flux_derivatives_generalized = 0.0

  call convert_to_generalized(nvars, ndim, nx, ny, nz, &
    fluxes, metrics, metric_jacobians, generalized_fluxes)

  do k=1,nz
    do j=1,ny
      do i=0,nx
        call pointwise_eigenvalues(nvars, 1, states(:, i-2:i+3, j, k), lambda_pointwise)
        call roe_eigensystem(nvars, 1, states(:, i:i+1, j, k), R, R_inv, lambda_roe)

        call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, lambda)

        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), &
                                           combined_frozen_metrics, &
                                           generalized_states_frozen, &
                                           generalized_fluxes_frozen)
        call split_characteristic_fluxes(nvars, &
                                         generalized_states_frozen, &
                                         generalized_fluxes_frozen(:, 1, :), &
                                         R_inv, &
                                         lambda, &
                                         characteristic_fluxes_pos, &
                                         characteristic_fluxes_neg)

        call weno_flux(nvars, &
                        generalized_fluxes(:, 1, i-2:i+3, j, k), &
                        characteristic_fluxes_pos, characteristic_fluxes_neg, &
                        combined_frozen_metrics(1), R, flux)

        do v=1,nvars
          flux_derivatives_generalized(v, 1, i, j, k) = &
            flux_derivatives_generalized(v, 1, i, j, k) + flux(v) / delta_xi
          flux_derivatives_generalized(v, 1, i+1, j, k) = &
            flux_derivatives_generalized(v, 1, i+1, j, k) - flux(v) / delta_xi
        end do
      end do
    end do
  end do

  ! two more loops:
  !   j is inner index, filling flux_derivatives_generalized(:,2,:,:,:), using delta_eta
  !   k is inner index, filling flux_derivatives_generalized(:,3,:,:,:), using delta_zeta

  call convert_from_generalized(nvars, ndim, nx, ny, nz, &
                                flux_derivatives_generalized, &
                                metrics, &
                                metric_jacobians, &
                                flux_derivatives)
end subroutine

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

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

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

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

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

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

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

  real u(ndim), rho, c, H, q
  real u_tilde(ndim), k(ndim), l(ndim), m(ndim)
  integer i
  real b1, b2, alpha, beta

  k = 
  l = 
  m = 

  u = 
  rho = 
  H = 
  c = 

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

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

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

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

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

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

  lambda_roe(1) = u(1)
  lambda_roe(2) = u(1)
  lambda_roe(3) = u(1)
  lambda_roe(4) = u(1) + c
  lambda_roe(5) = u(1) - c

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

  R(1,2) = 0.0
  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.0
  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.0 - 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.0

  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.0

  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, lambda)
  implicit none

  integer, intent(in) :: nvars
  real, intent(in) :: lambda_pointwise(nvars, -2:3)
  real, intent(in) :: lambda_roe(nvars)
  real, intent(out) :: lambda(nvars)

  real kappa
  integer v
  integer k

  kappa = 1.1

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

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

  real metrics_frozen(ndim, ndim)
  real jacobian_frozen

  integer c, k, v
  integer i, j

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

  do i=1,ndim
    combined_frozen_metrics(i) = 0.0
    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
      call mult_mat_mat(ndim, 1, ndim, 1.0, metrics_frozen, &
        fluxes(v,:,k), generalized_fluxes_frozen(v,:,k))
    end do
  end do
end subroutine

subroutine convert_to_generalized(nvars, ndim, nx, ny, nz, &
                                  fluxes, &
                                  metrics, &
                                  metric_jacobians, &
                                  generalized_fluxes)
  implicit none

  integer, intent(in) :: nvars, ndim, nx, ny, nz
  real, intent(in) :: fluxes(nvars, ndim, nx, ny, nz)
  real, intent(in) :: metrics(ndim, ndim, nx, ny, nz)
  real, intent(in) :: metric_jacobians(nx, ny, nz)
  real, intent(out) :: generalized_fluxes(nvars, ndim, nx, ny, nz)

  integer i, j, k, v

  do k=1,nz
    do j=1,ny
      do i=1,nx
        do v=1,nvars
          call mult_mat_mat(ndim, 1, ndim, 1.0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
            fluxes(v,:,i,j,k), generalized_fluxes(v,:,i,j,k))
        end do
      end do
    end do
  end do
end subroutine

subroutine convert_from_generalized(nvars, ndim, nx, ny, nz, &
                                    flux_derivatives_generalized, &
                                    metrics, &
                                    metric_jacobians, &
                                    flux_derivatives)
  implicit none

  integer, intent(in) :: nvars, ndim, nx, ny, nz
  real, intent(in) :: flux_derivatives_generalized(nvars, ndim, nx, ny, nz)
  real, intent(in) :: metrics(ndim, ndim, nx, ny, nz)
  real, intent(in) :: metric_jacobians(nx, ny, nz)
  real, intent(out) :: flux_derivatives(nvars, ndim, nx, ny, nz)

  integer i, j, k, v
  integer d

  do k=1,nz
    do j=1,ny
      do i=1,nx
        do v=1,nvars
          call solve3x3(metrics(:,:,i,j,k), flux_derivatives_generalized(v,:,i,j,k), &
            flux_derivatives(v,:,i,j,k))
          do d=1,ndim
            flux_derivatives(v,d,i,j,k) = flux_derivatives(v,d,i,j,k)*metric_jacobians(i,j,k)
          end do
        end do
      end do
    end do
  end do
end subroutine

subroutine solve3x3(a, b, x)
  implicit none

  real, intent(in) :: a(3,3), b(3)
  real, intent(out) :: x(3)

  real temp(3,3)
  real det_a, det_temp
  integer k
  integer i,j

  call determinant3x3(a, det_a)
  do k=1,3
    do j=1,3
      do i=1,3
        temp(i,j) = a(i,j)
        temp(i,k) = b(k)
      end do
    end do
    call determinant3x3(temp, det_temp)
    x(k) = det_temp/det_a
  end do
end subroutine

subroutine determinant3x3(a, det)
  implicit none

  real, intent(in) :: a(3,3)
  real, intent(out) :: det

  det = a(1,1)*(a(2,2)*a(3,3) - a(2,3)*a(3,2)) &
        + a(1,2)*(a(2,3)*a(3,1) - a(2,1)*a(3,3)) &
        + a(1,3)*(a(2,1)*a(3,2) - a(2,2)*a(3,1))
end subroutine

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

  integer, intent(in) :: nvars
  real, intent(in) :: generalized_states_frozen(nvars, -2:3)
  real, intent(in) :: generalized_fluxes_frozen(nvars, -2:3)
  real, intent(in) :: R_inv(nvars, nvars)
  real, intent(in) :: lambda(nvars)
  real, intent(out) :: characteristic_fluxes_pos(nvars, -2:3)
  real, 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.0
      characteristic_fluxes_neg(m,k) = 0.0
      do l=1,nvars
        characteristic_fluxes_pos(m,k) = characteristic_fluxes_pos(m,k) &
          + 0.5*R_inv(m,l) &
          *(generalized_fluxes_frozen(l,k) + lambda(m)*generalized_states_frozen(l,k))
        characteristic_fluxes_neg(m,k) = characteristic_fluxes_neg(m,k) &
          + 0.5*R_inv(m,l) &
          *(generalized_fluxes_frozen(l,k) - lambda(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_metrics, R, flux)
  implicit none

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

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

  call consistent_part(nvars, generalized_fluxes, consistent)
  call dissipation_part_pos(nvars, characteristic_fluxes_pos, &
    combined_frozen_metrics, R, dissipation_pos)
  call dissipation_part_neg(nvars, characteristic_fluxes_neg, &
    combined_frozen_metrics, 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, intent(in) :: generalized_fluxes(nvars, -2:3)
  real, 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)
  end do
end subroutine

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

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

  real combined_fluxes(nvars)

  call weno_combination_pos(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes)

  call mult_mat_mat(nvars, 1, nvars, -1.0/60, R, combined_fluxes, dissipation_pos)
end subroutine

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

  integer, intent(in) :: nvars
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, intent(in) :: combined_frozen_metrics
  real, intent(out) :: combined_fluxes(nvars)

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

  call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metrics, 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_metrics, w)
  implicit none

  integer, intent(in) :: nvars
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, intent(in) :: combined_frozen_metrics
  real, intent(out) :: w(nvars, 3)

  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

  integer p, i, j
  real sum_alpha

  eps = 1e-6*combined_frozen_metrics
  p = 2

  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)

    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
    end do
    call sum_array(3, alpha, sum_alpha)
    do j=1,3
      w(i,j) = alpha(j)/sum_alpha
    end do
  end do
end subroutine

subroutine flux_differences_pos(nvars, characteristic_fluxes, flux_differences)
  implicit none

  integer, intent(in) :: nvars
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, 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_metrics, R, dissipation_neg)
  implicit none

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

  real combined_fluxes(nvars)

  call weno_combination_neg(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes)

  call mult_mat_mat(nvars, 1, nvars, 1.0/60, R, combined_fluxes, dissipation_neg)
end subroutine

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

  integer, intent(in) :: nvars
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, intent(in) :: combined_frozen_metrics
  real, intent(out) :: combined_fluxes(nvars)

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

  call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metrics, 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_metrics, w)
  implicit none

  integer, intent(in) :: nvars
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, intent(in) :: combined_frozen_metrics
  real, intent(out) :: w(nvars, 3)

  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

  integer p, i, j
  real sum_alpha

  eps = 1e-6*combined_frozen_metrics
  p = 2

  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)

    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
    end do
    call sum_array(3, alpha, sum_alpha)
    do j=1,3
      w(i,j) = alpha(j)/sum_alpha
    end do
  end do
end subroutine

subroutine flux_differences_neg(nvars, characteristic_fluxes, flux_differences)
  implicit none

  integer, intent(in) :: nvars
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, 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_mat(m, n, l, alpha, a, b, c)
  implicit none

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

  integer i, j, k

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

subroutine sum_array(n, a, a_sum)
  implicit none

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

  integer i

  a_sum = 0.0
  do i=1,n
    a_sum = a_sum + a(i)
  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)
! RESULT = prg
!
!$loopy end