subroutine compute_flux_derivatives(states, fluxes, flux_derivatives)
  real, intent(in) :: states(nvars, nx, ny, nz)
  real, intent(in) :: fluxes(nvars, 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(ndim, nx, ny, nz)

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

  flux_derivatives_generalized = 0.0

  do k=1,nz
    do j=1,ny
      do i=0,nx
        call weno_flux(states(:, i-2:i+3, j, k),
                       fluxes(:, i-2:i+3, j, k),
                       metrics(:, :, i-2:i+3, j, k),
                       metric_jacobians(i-2:i+3, j, k),
                       flux)

        flux_derivatives_generalized(1, i, j, k) += flux / delta_xi
        flux_derivatives_generalized(1, i+1, j, k) -= flux / delta_xi
      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(flux_derivatives_generalized,
                                metrics,
                                metric_jacobians,
                                flux_derivatives)
end subroutine

subroutine weno_flux(states, fluxes, metrics, metric_jacobians, flux)
  real, intent(in) :: states(nvars, -2:3)
  real, intent(in) :: fluxes(nvars, -2:3)
  real, intent(in) :: metrics(ndim, ndim, -2:3)
  real, intent(in) :: metric_jacobians(-2:3)
  real, intent(out) :: flux

  ! compute frozen metrics
  call roe_eigenvectors(..., R)
  ! compute characteristic fluxes
  call split_characteristic_fluxes(..., characteristic_fluxes_pos, characteristic_fluxes_neg)

  call consistent_part(fluxes, metrics, metric_jacobians, consistent)
  call dissipation_part_pos(characteristic_fluxes_pos, combined_frozen_metrics, R, dissipation_pos)
  call dissipation_part_neg(characteristic_fluxes_neg, combined_frozen_metrics, R, dissipation_neg)

  flux = consistent + dissipation_pos + dissipation_neg
end subroutine

subroutine roe_eigenvectors(..., R)
  ! returns eigenvectors of Roe matrix
end subroutine

subroutine consistent_part(fluxes, metrics, metric_jacobians, consistent)
  call convert_to_generalized(fluxes, metrics, metric_jacobians, generalized_fluxes)

  consistent = sum([1, -8, 37, 37, -8, 1]*generalized_fluxes)/60
end subroutine

subroutine convert_to_generalized(fluxes,
                                  metrics,
                                  metric_jacobians,
                                  generalized_fluxes)
  real, intent(in) :: fluxes(nvars, -2:3)
  real, intent(in) :: metrics(ndim, ndim, -2:3)
  real, intent(in) :: metric_jacobians(-2:3)
  real, intent(out) :: generalized_fluxes(-2:3)

  ! transforms fluxes from physical coordinates to generalized coordinates
  ! uses metrics at each point
end subroutine

subroutine dissipation_part_pos(states, fluxes, metrics, metric_jacobians, dissipation_pos)
  real, intent(in) :: characteristic_fluxes(-2:3)
  real, intent(in) :: combined_frozen_metrics
  real, intent(in) :: R(nvars, nvars)

  call weno_combination_pos(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)

  dissipation_pos = -matmul(R, combined_fluxes)/60
end subroutine

subroutine weno_combination_pos(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
  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)

  call weno_weights_pos(characteristic_fluxes, combined_frozen_metrics, w)
  call flux_differences_pos(characteristic_fluxes, flux_differences)

  combined_fluxes = (20*w(:,1) - 1)*flux_differences(:,1)
                      - (10*(w(:,1) + w(:,2)) - 5)*flux_differences(:,2)
                      + flux_differences(:,3)
end subroutine

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

  C = [0.1, 0.6, 0.3]
  eps = 1e-6*combined_frozen_metrics
  p = 2

  do i=1,nvars
    IS(1) = (1/4)*(sum([1, -4, 3]*characteristic_fluxes(i,-2:0)))**2
            + (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i,-2:0)))**2
    IS(2) = (1/4)*(sum([-1, 0, 1]*characteristic_fluxes(i,-1:1)))**2
            + (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i,-1:1)))**2
    IS(3) = (1/4)*(sum([-3, 4, -1]*characteristic_fluxes(i,0:2)))**2
            + (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i,0:2)))**2

    alpha = C/(IS + eps)**p
    w(i,:) = alpha/sum(alpha)
  end do
end subroutine

subroutine flux_differences_pos(characteristic_fluxes, flux_differences)
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, intent(out) :: flux_differences(nvars, 3)

  do i=1,3
    flux_differences(:,i) = sum([-1, 3, -3, 1]*characteristic_fluxes(:,-3+i:i))
  end do
end subroutine

subroutine dissipation_part_neg(states, fluxes, metrics, metric_jacobians, dissipation_neg)
  real, intent(in) :: characteristic_fluxes(-2:3)
  real, intent(in) :: combined_frozen_metrics
  real, intent(in) :: R(nvars, nvars)

  call weno_combination_neg(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)

  dissipation_neg = matmul(R, combined_fluxes)/60
end subroutine

subroutine weno_combination_neg(characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
  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)

  call weno_weights_neg(characteristic_fluxes, combined_frozen_metrics, w)
  call flux_differences_neg(characteristic_fluxes, flux_differences)

  combined_fluxes = (20*w(:,1) - 1)*flux_differences(:,1)
                      - (10*(w(:,1) + w(:,2)) - 5)*flux_differences(:,2)
                      + flux_differences(:,3)
end subroutine

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

  C = [0.1, 0.6, 0.3]
  eps = 1e-6*combined_frozen_metrics
  p = 2

  do i=1,nvars
    IS(1) = (1/4)*(sum([1, -4, 3]*characteristic_fluxes(i, 3:1)))**2
            + (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 3:1)))**2
    IS(2) = (1/4)*(sum([-1, 0, 1]*characteristic_fluxes(i, 2:0)))**2
            + (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 2:0)))**2
    IS(3) = (1/4)*(sum([-3, 4, -1]*characteristic_fluxes(i, 1:-1)))**2
            + (13/12)*(sum([1, -2, 1]*characteristic_fluxes(i, 1:-1)))**2

    alpha = C/(IS + eps)**p
    w(i,:) = alpha/sum(alpha)
  end do
end subroutine

subroutine flux_differences_neg(characteristic_fluxes, flux_differences)
  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
  real, intent(out) :: flux_differences(nvars, 3)

  do i=1,3
    flux_differences(:,i) = sum([-1, 3, -3, 1]*characteristic_fluxes(:,1-i:4-i))
  end do
end subroutine

subroutine convert_from_generalized(flux_derivatives_generalized,
                                    metrics,
                                    metric_jacobians,
                                    flux_derivatives)
  ! transforms flux derivatives from generalized coordinates to physical coordinates
end subroutine