//========================================================================================
// Athena++ astrophysical MHD code
// Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
// Licensed under the 3-clause BSD License, see LICENSE file for details
//========================================================================================
//! \file ppm_simple.cpp
//  \brief piecewise parabolic reconstruction with modified McCorquodale/Colella limiter
//         for a Cartesian-like coordinate with uniform spacing, Mignone modified original
//         PPM limiter for nonuniform and/or curvilinear coordinate.
//  Operates on the entire nx4 range of a single AthenaArray<Real> input (no MHD).
//  No assumptions of hydrodynamic fluid variable input; no characteristic projection.

// REFERENCES:
// (CW) P. Colella & P. Woodward, "The Piecewise Parabolic Method (PPM) for Gas-Dynamical
// Simulations", JCP, 54, 174 (1984)
//
// (CS) P. Colella & M. Sekora, "A limiter for PPM that preserves accuracy at smooth
// extrema", JCP, 227, 7069 (2008)
//
// (MC) P. McCorquodale & P. Colella,  "A high-order finite-volume method for conservation
// laws on locally refined grids", CAMCoS, 6, 1 (2011)
//
// (CD) P. Colella, M.R. Dorr, J. Hittinger, D. Martin, "High-order, finite-volume methods
// in mapped coordinates", JCP, 230, 2952 (2011)
//
// (Mignone) A. Mignone, "High-order conservative reconstruction schemes for finite volume
// methods in cylindrical and spherical coordinates", JCP, 270, 784 (2014)
//========================================================================================

// C headers

// C++ headers
#include <algorithm>    // max()

// Athena++ headers
#include "../athena.hpp"
#include "../athena_arrays.hpp"
#include "../coordinates/coordinates.hpp"
#include "../eos/eos.hpp"
#include "reconstruction.hpp"

//----------------------------------------------------------------------------------------
//! \fn Reconstruction::PiecewiseParabolicX1()
//  \brief Returns L/R interface values in X1-dir constructed using fourth-order PPM and
//         Colella-Sekora or Mignone limiting over [kl,ku][jl,ju][il,iu]

void Reconstruction::PiecewiseParabolicX1(
    const int k, const int j, const int il, const int iu,
    const AthenaArray<Real> &q,
    AthenaArray<Real> &ql, AthenaArray<Real> &qr) {
  const int nu = q.GetDim4() - 1;

  // CS08 constant used in second derivative limiter, >1 , independent of h
  const Real C2 = 1.25;

  // TODO(felker): renumber scratch array references; not using 2x from ppm.cpp
  // bx (MHD) and wc (characteristic projection)

  // set work arrays used for primitive/characterstic cell-averages to scratch
  AthenaArray<Real> &q_im2 = scr2_ni_, &q_im1 = scr3_ni_,
                     &q_i = scr4_ni_, &q_ip1 = scr5_ni_, &q_ip2 = scr6_ni_,
                &qr_imh = scr7_ni_, &ql_iph = scr8_ni_;

  // set work PPM work arrays to shallow copies of scratch arrays:
  AthenaArray<Real> &dd = scr02_i_, &dd_im1 = scr03_i_, &dd_ip1 = scr04_i_,
                   &dph = scr05_i_, &dph_ip1 = scr06_i_;

  AthenaArray<Real> &d2qc_im1 = scr07_i_, &d2qc = scr08_i_, &d2qc_ip1 = scr09_i_,
                        &d2qf = scr10_i_;

  AthenaArray<Real> &qplus = scr11_i_, &qminus = scr12_i_, &dqf_plus = scr13_i_,
                &dqf_minus = scr14_i_;

  // cache the x1-sliced primitive states for eigensystem calculation
  for (int n=0; n<=nu; ++n) {
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      q_i  (n,i) = q(n,k,j,i  );
      q_im2(n,i) = q(n,k,j,i-2);
      q_im1(n,i) = q(n,k,j,i-1);
      q_ip1(n,i) = q(n,k,j,i+1);
      q_ip2(n,i) = q(n,k,j,i+2);
    }
  }

  //--- Step 1. --------------------------------------------------------------------------
  // Reconstruct interface averages <a>_{i-1/2} and <a>_{i+1/2}
  for (int n=0; n<=nu; ++n) {
    // Compute average slope in i-1, i, i+1 zones
#pragma omp simd simdlen(SIMD_WIDTH)
    for (int i=il; i<=iu; ++i) {
      // nonuniform or uniform Cartesian-like coord reconstruction from volume averages:
      if (!curvilinear[X1DIR]) {
        Real qa = (q_i(n,i) - q_im1(n,i));
        Real qb = (q_ip1(n,i) - q_i(n,i));
        dd_im1(i) = c1i(i-1)*qa + c2i(i-1)*(q_im1(n,i) - q_im2(n,i));
        dd    (i) = c1i(i  )*qb + c2i(i  )*qa;
        dd_ip1(i) = c1i(i+1)*(q_ip2(n,i) - q_ip1(n,i)) + c2i(i+1)*qb;

        // Approximate interface average at i-1/2 and i+1/2 using PPM (CW eq 1.6)
        // KGF: group the biased stencil quantities to preserve FP symmetry
        dph(i) = (c3i(i)*q_im1(n,i) + c4i(i)*q_i(n,i)) +
                 (c5i(i)*dd_im1(i) + c6i(i)*dd(i));
        dph_ip1(i) = (c3i(i+1)*q_i(n,i) + c4i(i+1)*q_ip1(n,i)) +
                     (c5i(i+1)*dd(i) + c6i(i+1)*dd_ip1(i) );
      } else { // radial coordinate
        dph(i) = c1i(i)*q_im2(n,i) + c2i(i)*q_im1(n,i) + c3i(i)*q_i(n,i)
                 + c4i(i)*q_ip1(n,i);
        dph_ip1(i) = c1i(i+1)*q_im1(n,i) + c2i(i+1)*q_i(n,i) + c3i(i+1)*q_ip1(n,i)
                     + c4i(i+1)*q_ip2(n,i);
      }
    }

    //--- Step 2a. -----------------------------------------------------------------------
    // Uniform Cartesian-like coordinate: limit interpolated interface states (CD 4.3.1)
    if (uniform[X1DIR] && !curvilinear[X1DIR]) {
      // approximate second derivative at interfaces for smooth extrema preservation
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu+1; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
        d2qc_im1(i) = q_im2(n,i) + q_i  (n,i) - 2.0*q_im1(n,i);
        d2qc    (i) = q_im1(n,i) + q_ip1(n,i) - 2.0*q_i  (n,i); // (CD eq 85a) (no 1/2)
        d2qc_ip1(i) = q_i  (n,i) + q_ip2(n,i) - 2.0*q_ip1(n,i);
      }

      // i-1/2
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=(iu+1); ++i) {
        Real qa_tmp = dph(i) - q_im1(n,i); // (CD eq 84a)
        Real qb_tmp = q_i(n,i) - dph(i);     // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
        Real qa = 3.0*(q_im1(n,i) + q_i(n,i)  - 2.0*dph(i));  // (CD eq 85b)
        Real qb = d2qc_im1(i);    // (CD eq 85a) (no 1/2)
        Real qc = d2qc(i);   // (CD eq 85c) (no 1/2)
        Real qd = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
          qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
        }
        Real dph_tmp = 0.5*(q_im1(n,i) + q_i(n,i)) - qd/6.0;
        if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at i-1/2 face
          dph(i) = dph_tmp;
        }
      }

      // i+1/2
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=(iu+1); ++i) {
        Real qa_tmp = dph_ip1(i) - q_i(n,i);       // (CD eq 84a)
        Real qb_tmp = q_ip1(n,i) - dph_ip1(i);   // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
        Real qa = 3.0*(q_i(n,i) + q_ip1(n,i) - 2.0*dph_ip1(i));  // (CD eq 85b)
        Real qb = d2qc(i);            // (CD eq 85a) (no 1/2)
        Real qc = d2qc_ip1(i);   // (CD eq 85c) (no 1/2)
        Real qd = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
          qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
        }
        Real dphip1_tmp = 0.5*(q_i(n,i) + q_ip1(n,i)) - qd/6.0;
        if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at i+1/2 face
          dph_ip1(i) = dphip1_tmp;
        }
      }

#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
        d2qf(i) = 6.0*(dph(i) + dph_ip1(i) - 2.0*q_i(n,i)); // a6 coefficient * -2
      }

      //--- Step 2b. ---------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply strict monotonicity constraints
      // (Mignone eq 45)
    } else {
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        dph    (i) = std::min(dph    (i), std::max(q_i(n,i),q_im1(n,i)));
        dph_ip1(i) = std::min(dph_ip1(i), std::max(q_i(n,i),q_ip1(n,i)));

        dph    (i) = std::max(dph    (i), std::min(q_i(n,i),q_im1(n,i)));
        dph_ip1(i) = std::max(dph_ip1(i), std::min(q_i(n,i),q_ip1(n,i)));
      }
    }

    // Cache Riemann states for both non-/uniform limiters
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      qminus(i) = dph(i  );
      qplus(i) =  dph_ip1(i );
    }

    //--- Step 3. ------------------------------------------------------------------------
    // Compute cell-centered difference stencils (MC section 2.4.1)
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      dqf_minus(i) = q_i(n,i) - qminus(i); // (CS eq 25) = -dQ^- in Mignone's notation
      dqf_plus(i)  = qplus(i) - q_i(n,i);
    }

    //--- Step 4a. -----------------------------------------------------------------------
    // For uniform Cartesian-like coordinate: apply CS limiters to parabolic interpolant
    if (uniform[X1DIR] && !curvilinear[X1DIR]) {
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        Real qa_tmp = dqf_minus(i)*dqf_plus(i);
        Real qb_tmp = (q_ip1(n,i) - q_i(n,i))*(q_i(n,i) - q_im1(n,i));

        Real qa = d2qc_im1(i);
        Real qb = d2qc(i);
        Real qc = d2qc_ip1(i);
        Real qd = d2qf(i);
        Real qe = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc) && SIGN(qa) == SIGN(qd)) {
          // Extrema is smooth
          qe = SIGN(qd)* std::min(std::min(C2*std::abs(qa), C2*std::abs(qb)),
                                  std::min(C2*std::abs(qc),
                                           std::abs(qd))); // (CS eq 22)
        }

        // Check if 2nd derivative is close to roundoff error
        qa = std::max(std::abs(q_im1(n,i)), std::abs(q_im2(n,i)));
        qb = std::max(std::max(std::abs(q_i(n,i)), std::abs(q_ip1(n,i))),
                      std::abs(q_ip2(n,i)));

        Real rho = 0.0;
        if (std::abs(qd) > (1.0e-12)*std::max(qa,qb)) {
          // Limiter is not sensitive to roundoff. Use limited ratio (MC eq 27)
          rho = qe/qd;
        }

        Real tmp_m = q_i(n,i) - rho*dqf_minus(i);
        Real tmp_p = q_i(n,i) + rho*dqf_plus(i);
        Real tmp2_m = q_i(n,i) - 2.0*dqf_plus(i);
        Real tmp2_p = q_i(n,i) + 2.0*dqf_minus(i);

        // Check for local extrema
        if ((qa_tmp <= 0.0 || qb_tmp <= 0.0)) {
          // Check if relative change in limited 2nd deriv is > roundoff
          if (rho <= (1.0 - (1.0e-12))) {
            // Limit smooth extrema
            qminus(i) = tmp_m;// (CS eq 23)
            qplus(i) = tmp_p;
          }
          // No extrema detected
        } else {
          // Overshoot i-1/2,R / i,(-) state
          if (std::abs(dqf_minus(i)) >= 2.0*std::abs(dqf_plus(i))) {
            qminus(i) = tmp2_m;
          }
          // Overshoot i+1/2,L / i,(+) state
          if (std::abs(dqf_plus(i)) >= 2.0*std::abs(dqf_minus(i))) {
            qplus(i) = tmp2_p;
          }
        }
      }

      //--- Step 4b. ---------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply Mignone limiters to parabolic
      // interpolant. Note, the Mignone limiter does not check for cell-averaged extrema:
    } else {
      for (int i=il; i<=iu; ++i) {
        Real qa = dqf_minus(i)*dqf_plus(i);
        if (qa <= 0.0) { // Local extrema detected
          qminus(i) = q_i(n,i);
          qplus(i) = q_i(n,i);
        } else { // No extrema detected
          // Overshoot i-1/2,R / i,(-) state
          if (std::abs(dqf_minus(i)) >= hplus_ratio_i(i)*std::abs(dqf_plus(i))) {
            qminus(i) = q_i(n,i) - hplus_ratio_i(i)*dqf_plus(i);
          }
          // Overshoot i+1/2,L / i,(+) state
          if (std::abs(dqf_plus(i)) >= hminus_ratio_i(i)*std::abs(dqf_minus(i))) {
            qplus(i) = q_i(n,i) + hminus_ratio_i(i)*dqf_minus(i);
          }
        }
      }
    }
    //--- Step 5. ------------------------------------------------------------------------
    // Convert limited cell-centered values to interface-centered L/R Riemann states
    // both L/R values defined over [il,iu]
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      ql_iph(n,i ) = qplus(i);
      qr_imh(n,i ) = qminus(i);
    }
  } // end char PPM loop over =nu

  // compute ql_(i+1/2) and qr_(i-1/2)
  for (int n=0; n<=nu; ++n) {
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      ql(n,i+1) = ql_iph(n,i);
      qr(n,i  ) = qr_imh(n,i);
    }
  }
  return;
}



void Reconstruction::PiecewiseParabolicX1(
    const int k, const int j, const int il, const int iu,
    const AthenaArray<Real> &q, const int array_order,
    AthenaArray<Real> &ql, AthenaArray<Real> &qr) {

  if(array_order < 0){

    const int nu = q.GetDim1() - 1;

  // CS08 constant used in second derivative limiter, >1 , independent of h
    const Real C2 = 1.25;

  // TODO(felker): renumber scratch array references; not using 2x from ppm.cpp
  // bx (MHD) and wc (characteristic projection)

  // set work arrays used for primitive/characterstic cell-averages to scratch
    AthenaArray<Real> &q_im2 = scr2_in_, &q_im1 = scr3_in_,
                    &q_i = scr4_in_, &q_ip1 = scr5_in_, &q_ip2 = scr6_in_,
                    &qr_imh = scr7_in_, &ql_iph = scr8_in_;

  // set work PPM work arrays to shallow copies of scratch arrays:
    AthenaArray<Real> &dd = scr02_i_, &dd_im1 = scr03_i_, &dd_ip1 = scr04_i_,
                   &dph = scr05_i_, &dph_ip1 = scr06_i_;

    AthenaArray<Real> &d2qc_im1 = scr07_i_, &d2qc = scr08_i_, &d2qc_ip1 = scr09_i_,
                        &d2qf = scr10_i_;

    AthenaArray<Real> &qplus = scr11_i_, &qminus = scr12_i_, &dqf_plus = scr13_i_,
                &dqf_minus = scr14_i_;

  // cache the x1-sliced primitive states for eigensystem calculation
    for (int i=il; i<=iu; ++i) {
#pragma omp simd
      for (int n=0; n<=nu; ++n) {
        q_i  (n,i) = q(k,j,i  ,n);
        q_im2(n,i) = q(k,j,i-2,n);
        q_im1(n,i) = q(k,j,i-1,n);
        q_ip1(n,i) = q(k,j,i+1,n);
        q_ip2(n,i) = q(k,j,i+2,n);
      }
    }

  //--- Step 1. --------------------------------------------------------------------------
  // Reconstruct interface averages <a>_{i-1/2} and <a>_{i+1/2}
    for (int n=0; n<=nu; ++n) {
    // Compute average slope in i-1, i, i+1 zones
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
      // nonuniform or uniform Cartesian-like coord reconstruction from volume averages:
        if (!curvilinear[X1DIR]) {
          Real qa = (q_i(n,i) - q_im1(n,i));
          Real qb = (q_ip1(n,i) - q_i(n,i));
          dd_im1(i) = c1i(i-1)*qa + c2i(i-1)*(q_im1(n,i) - q_im2(n,i));
          dd    (i) = c1i(i  )*qb + c2i(i  )*qa;
          dd_ip1(i) = c1i(i+1)*(q_ip2(n,i) - q_ip1(n,i)) + c2i(i+1)*qb;

        // Approximate interface average at i-1/2 and i+1/2 using PPM (CW eq 1.6)
        // KGF: group the biased stencil quantities to preserve FP symmetry
          dph(i) = (c3i(i)*q_im1(n,i) + c4i(i)*q_i(n,i)) +
                 (c5i(i)*dd_im1(i) + c6i(i)*dd(i));
          dph_ip1(i) = (c3i(i+1)*q_i(n,i) + c4i(i+1)*q_ip1(n,i)) +
                     (c5i(i+1)*dd(i) + c6i(i+1)*dd_ip1(i) );
        } else { // radial coordinate
          dph(i) = c1i(i)*q_im2(n,i) + c2i(i)*q_im1(n,i) + c3i(i)*q_i(n,i)
                 + c4i(i)*q_ip1(n,i);
          dph_ip1(i) = c1i(i+1)*q_im1(n,i) + c2i(i+1)*q_i(n,i) + c3i(i+1)*q_ip1(n,i)
                     + c4i(i+1)*q_ip2(n,i);
        }
      }

    //--- Step 2a. -----------------------------------------------------------------------
    // Uniform Cartesian-like coordinate: limit interpolated interface states (CD 4.3.1)
      if (uniform[X1DIR] && !curvilinear[X1DIR]) {
      // approximate second derivative at interfaces for smooth extrema preservation
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu+1; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
          d2qc_im1(i) = q_im2(n,i) + q_i  (n,i) - 2.0*q_im1(n,i);
          d2qc    (i) = q_im1(n,i) + q_ip1(n,i) - 2.0*q_i  (n,i); // (CD eq 85a) (no 1/2)
          d2qc_ip1(i) = q_i  (n,i) + q_ip2(n,i) - 2.0*q_ip1(n,i);
        }

      // i-1/2
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=(iu+1); ++i) {
          Real qa_tmp = dph(i) - q_im1(n,i); // (CD eq 84a)
          Real qb_tmp = q_i(n,i) - dph(i);     // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
          Real qa = 3.0*(q_im1(n,i) + q_i(n,i)  - 2.0*dph(i));  // (CD eq 85b)
          Real qb = d2qc_im1(i);    // (CD eq 85a) (no 1/2)
          Real qc = d2qc(i);   // (CD eq 85c) (no 1/2)
          Real qd = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
            qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
          }
          Real dph_tmp = 0.5*(q_im1(n,i) + q_i(n,i)) - qd/6.0;
          if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at i-1/2 face
            dph(i) = dph_tmp;
          }
        }

      // i+1/2
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=(iu+1); ++i) {
          Real qa_tmp = dph_ip1(i) - q_i(n,i);       // (CD eq 84a)
          Real qb_tmp = q_ip1(n,i) - dph_ip1(i);   // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
          Real qa = 3.0*(q_i(n,i) + q_ip1(n,i) - 2.0*dph_ip1(i));  // (CD eq 85b)
          Real qb = d2qc(i);            // (CD eq 85a) (no 1/2)
          Real qc = d2qc_ip1(i);   // (CD eq 85c) (no 1/2)
          Real qd = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
            qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
          }
          Real dphip1_tmp = 0.5*(q_i(n,i) + q_ip1(n,i)) - qd/6.0;
          if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at i+1/2 face
            dph_ip1(i) = dphip1_tmp;
          }
        }

#pragma omp simd
        for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
          d2qf(i) = 6.0*(dph(i) + dph_ip1(i) - 2.0*q_i(n,i)); // a6 coefficient * -2
        }

      //--- Step 2b. ---------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply strict monotonicity constraints
      // (Mignone eq 45)
      } else {
#pragma omp simd
        for (int i=il; i<=iu; ++i) {
          dph    (i) = std::min(dph    (i), std::max(q_i(n,i),q_im1(n,i)));
          dph_ip1(i) = std::min(dph_ip1(i), std::max(q_i(n,i),q_ip1(n,i)));

          dph    (i) = std::max(dph    (i), std::min(q_i(n,i),q_im1(n,i)));
          dph_ip1(i) = std::max(dph_ip1(i), std::min(q_i(n,i),q_ip1(n,i)));
        }
      }

    // Cache Riemann states for both non-/uniform limiters
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        qminus(i) = dph(i  );
        qplus(i) =  dph_ip1(i );
      }

    //--- Step 3. ------------------------------------------------------------------------
    // Compute cell-centered difference stencils (MC section 2.4.1)
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        dqf_minus(i) = q_i(n,i) - qminus(i); // (CS eq 25) = -dQ^- in Mignone's notation
        dqf_plus(i)  = qplus(i) - q_i(n,i);
      }

    //--- Step 4a. -----------------------------------------------------------------------
    // For uniform Cartesian-like coordinate: apply CS limiters to parabolic interpolant
      if (uniform[X1DIR] && !curvilinear[X1DIR]) {
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
          Real qa_tmp = dqf_minus(i)*dqf_plus(i);
          Real qb_tmp = (q_ip1(n,i) - q_i(n,i))*(q_i(n,i) - q_im1(n,i));

          Real qa = d2qc_im1(i);
          Real qb = d2qc(i);
          Real qc = d2qc_ip1(i);
          Real qd = d2qf(i);
          Real qe = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc) && SIGN(qa) == SIGN(qd)) {
          // Extrema is smooth
            qe = SIGN(qd)* std::min(std::min(C2*std::abs(qa), C2*std::abs(qb)),
                                  std::min(C2*std::abs(qc),
                                           std::abs(qd))); // (CS eq 22)
          }

        // Check if 2nd derivative is close to roundoff error
          qa = std::max(std::abs(q_im1(n,i)), std::abs(q_im2(n,i)));
          qb = std::max(std::max(std::abs(q_i(n,i)), std::abs(q_ip1(n,i))),
                      std::abs(q_ip2(n,i)));

          Real rho = 0.0;
          if (std::abs(qd) > (1.0e-12)*std::max(qa,qb)) {
          // Limiter is not sensitive to roundoff. Use limited ratio (MC eq 27)
            rho = qe/qd;
          }

          Real tmp_m = q_i(n,i) - rho*dqf_minus(i);
          Real tmp_p = q_i(n,i) + rho*dqf_plus(i);
          Real tmp2_m = q_i(n,i) - 2.0*dqf_plus(i);
          Real tmp2_p = q_i(n,i) + 2.0*dqf_minus(i);

        // Check for local extrema
          if ((qa_tmp <= 0.0 || qb_tmp <= 0.0)) {
          // Check if relative change in limited 2nd deriv is > roundoff
            if (rho <= (1.0 - (1.0e-12))) {
            // Limit smooth extrema
              qminus(i) = tmp_m;// (CS eq 23)
              qplus(i) = tmp_p;
            }
          // No extrema detected
          } else {
          // Overshoot i-1/2,R / i,(-) state
            if (std::abs(dqf_minus(i)) >= 2.0*std::abs(dqf_plus(i))) {
              qminus(i) = tmp2_m;
            }
          // Overshoot i+1/2,L / i,(+) state
            if (std::abs(dqf_plus(i)) >= 2.0*std::abs(dqf_minus(i))) {
              qplus(i) = tmp2_p;
            }
          }
        }

      //--- Step 4b. ---------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply Mignone limiters to parabolic
      // interpolant. Note, the Mignone limiter does not check for cell-averaged extrema:
      } else {
        for (int i=il; i<=iu; ++i) {
          Real qa = dqf_minus(i)*dqf_plus(i);
          if (qa <= 0.0) { // Local extrema detected
            qminus(i) = q_i(n,i);
            qplus(i) = q_i(n,i);
          } else { // No extrema detected
          // Overshoot i-1/2,R / i,(-) state
            if (std::abs(dqf_minus(i)) >= hplus_ratio_i(i)*std::abs(dqf_plus(i))) {
              qminus(i) = q_i(n,i) - hplus_ratio_i(i)*dqf_plus(i);
            }
          // Overshoot i+1/2,L / i,(+) state
            if (std::abs(dqf_plus(i)) >= hminus_ratio_i(i)*std::abs(dqf_minus(i))) {
              qplus(i) = q_i(n,i) + hminus_ratio_i(i)*dqf_minus(i);
            }
          }
        }
      }
    //--- Step 5. ------------------------------------------------------------------------
    // Convert limited cell-centered values to interface-centered L/R Riemann states
    // both L/R values defined over [il,iu]
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        ql_iph(n,i ) = qplus(i);
        qr_imh(n,i ) = qminus(i);
      }
    } // end char PPM loop over =nu

  // compute ql_(i+1/2) and qr_(i-1/2)
    for (int i=il; i<=iu; ++i) {
#pragma omp simd
      for (int n=0; n<=nu; ++n) {
        ql(i+1,n) = ql_iph(n,i);
        qr(i  ,n) = qr_imh(n,i);
      }
    }
  
  }// End array Ordering
  return;
}
//-------------------------------------------------------------------------------------
//! \fn Reconstruction::PiecewiseParabolicX2()
//  \brief Returns L/R interface values in X2-dir constructed using fourth-order PPM and
//         Colella-Sekora or Mignone limiting over [kl,ku][jl,ju][il,iu]

void Reconstruction::PiecewiseParabolicX2(
    const int k, const int j, const int il, const int iu,
    const AthenaArray<Real> &q,
    AthenaArray<Real> &ql, AthenaArray<Real> &qr) {
  const int nu = q.GetDim4() - 1;
  // CS08 constant used in second derivative limiter, >1 , independent of h
  const Real C2 = 1.25;

  // set work arrays used for primitive/characterstic cell-averages to scratch
  AthenaArray<Real> &q_jm2 = scr2_ni_, &q_jm1 = scr3_ni_,
                     &q_j = scr4_ni_, &q_jp1 = scr5_ni_, &q_jp2 = scr6_ni_,
                &qr_jmh = scr7_ni_, &ql_jph = scr8_ni_;

  // set work PPM work arrays to shallow copies of scratch arrays:
  AthenaArray<Real> &dd = scr02_i_, &dd_jm1 = scr03_i_, &dd_jp1 = scr04_i_,
                   &dph = scr05_i_, &dph_jp1 = scr06_i_;

  AthenaArray<Real> &d2qc_jm1 = scr07_i_, &d2qc = scr08_i_, &d2qc_jp1 = scr09_i_,
                        &d2qf = scr10_i_;

  AthenaArray<Real> &qplus = scr11_i_, &qminus = scr12_i_, &dqf_plus = scr13_i_,
                &dqf_minus = scr14_i_;

  // cache the x1-sliced primitive states for eigensystem calculation
  for (int n=0; n<=nu; ++n) {
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      q_j  (n,i) = q(n,k,j  ,i);
      q_jm2(n,i) = q(n,k,j-2,i);
      q_jm1(n,i) = q(n,k,j-1,i);
      q_jp1(n,i) = q(n,k,j+1,i);
      q_jp2(n,i) = q(n,k,j+2,i);
    }
  }

  //--- Step 1. ------------------------------------------------------------------------
  // Reconstruct interface averages <a>_{j-1/2} and <a>_{j+1/2}
  for (int n=0; n<=nu; ++n) {
    // Compute average slope in j-1, j, j+1 zones
#pragma omp simd simdlen(SIMD_WIDTH)
    for (int i=il; i<=iu; ++i) {
      // nonuniform or uniform Cartesian-like coord reconstruction from volume averages:
      if (!curvilinear[X2DIR]) {
        Real qa = (q_j(n,i) - q_jm1(n,i));
        Real qb = (q_jp1(n,i) - q_j(n,i));
        dd_jm1(i) = c1j(j-1)*qa + c2j(j-1)*(q_jm1(n,i) - q_jm2(n,i));
        dd    (i) = c1j(j  )*qb + c2j(j  )*qa;
        dd_jp1(i) = c1j(j+1)*(q_jp2(n,i) - q_jp1(n,i)) + c2j(j+1)*qb;

        // Approximate interface average at j-1/2 and j+1/2 using PPM (CW eq 1.6)
        // KGF: group the biased stencil quantities to preserve FP symmetry
        dph(i) = (c3j(j)*q_jm1(n,i) + c4j(j)*q_j(n,i)) +
                 (c5j(j)*dd_jm1(i) + c6j(j)*dd(i));
        dph_jp1(i) = (c3j(j+1)*q_j(n,i) + c4j(j+1)*q_jp1(n,i)) +
                     (c5j(j+1)*dd(i) + c6j(j+1)*dd_jp1(i));
      } else { // spherical-polar coordinates meridional x2 direction:
        // TODO(felker): implement weight calculation in reconstruction.cpp
        dph(i) = c1j(j)*q_jm2(n,i) + c2j(j)*q_jm1(n,i) + c3j(j)*q_j(n,i)
                 + c4j(j)*q_jp1(n,i);
        dph_jp1(i) = c1j(j+1)*q_jm1(n,i) + c2j(j+1)*q_j(n,i) + c3j(j+1)*q_jp1(n,i)
                     + c4j(j+1)*q_jp2(n,i);
      }
    }

    //--- Step 2a. ---------------------------------------------------------------------
    // Uniform Cartesian-like coordinate: limit interpolated interface states (CD 4.3.1)
    if (uniform[X2DIR] && !curvilinear[X2DIR]) {
      // approximate second derivative at interfaces for smooth extrema preservation
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
        d2qc_jm1(i) = q_jm2(n,i) + q_j  (n,i) - 2.0*q_jm1(n,i);
        d2qc    (i) = q_jm1(n,i) + q_jp1(n,i) - 2.0*q_j  (n,i); //(CD eq 85a) (no 1/2)
        d2qc_jp1(i) = q_j  (n,i) + q_jp2(n,i) - 2.0*q_jp1(n,i);
      }

      // j-1/2
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        Real qa_tmp = dph(i) - q_jm1(n,i); // (CD eq 84a)
        Real qb_tmp = q_j(n,i) - dph(i);     // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
        Real qa = 3.0*(q_jm1(n,i) + q_j(n,i) - 2.0*dph(i));  // (CD eq 85b)
        Real qb = d2qc_jm1(i);    // (CD eq 85a) (no 1/2)
        Real qc = d2qc(i);   // (CD eq 85c) (no 1/2)
        Real qd = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
          qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
        }
        Real dph_tmp = 0.5*(q_jm1(n,i) + q_j(n,i)) - qd/6.0;
        if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at j-1/2 face
          dph(i) = dph_tmp;
        }
      }
      // j+1/2
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        Real qa_tmp = dph_jp1(i) - q_j(n,i);       // (CD eq 84a)
        Real qb_tmp = q_jp1(n,i) - dph_jp1(i);   // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
        Real qa = 3.0*(q_j(n,i) + q_jp1(n,i)  - 2.0*dph_jp1(i));  // (CD eq 85b)
        Real qb = d2qc(i);            // (CD eq 85a) (no 1/2)
        Real qc = d2qc_jp1(i);   // (CD eq 85c) (no 1/2)
        Real qd = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
          qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
        }
        Real dphjp1_tmp = 0.5*(q_j(n,i) + q_jp1(n,i)) - qd/6.0;
        if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at j+1/2 face
          dph_jp1(i) = dphjp1_tmp;
        }
      }

#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
        d2qf(i) = 6.0*(dph(i) + dph_jp1(i) - 2.0*q_j(n,i)); // a6 coefficient * -2
      }

      //--- Step 2b. -------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply strict monotonicity constraints
      // (Mignone eq 45)
    } else {
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        dph    (i) = std::min(dph    (i), std::max(q_j(n,i),q_jm1(n,i)));
        dph_jp1(i) = std::min(dph_jp1(i), std::max(q_j(n,i),q_jp1(n,i)));

        dph    (i) = std::max(dph    (i), std::min(q_j(n,i),q_jm1(n,i)));
        dph_jp1(i) = std::max(dph_jp1(i), std::min(q_j(n,i),q_jp1(n,i)));
      }
    }

    // Cache Riemann states for both non-/uniform limiters
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      qminus(i) = dph(i  );
      qplus(i) =  dph_jp1(i );
    }

    //--- Step 3. ----------------------------------------------------------------------
    // Compute cell-centered difference stencils (MC section 2.4.1)
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      dqf_minus(i) = q_j(n,i) - qminus(i); // (CS eq 25)
      dqf_plus(i)  = qplus(i) - q_j(n,i);
    }

    //--- Step 4a. ---------------------------------------------------------------------
    // For uniform Cartesian-like coordinate: apply CS limiters to parabolic interpolant
    if (uniform[X2DIR] && !curvilinear[X2DIR]) {
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        Real qa_tmp = dqf_minus(i)*dqf_plus(i);
        Real qb_tmp = (q_jp1(n,i) - q_j(n,i))*(q_j(n,i) - q_jm1(n,i));

        Real qa = d2qc_jm1(i);
        Real qb = d2qc(i);
        Real qc = d2qc_jp1(i);
        Real qd = d2qf(i);
        Real qe = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc) && SIGN(qa) == SIGN(qd)) {
          // Extrema is smooth
          qe = SIGN(qd)* std::min(std::min(C2*std::abs(qa), C2*std::abs(qb)),
                                  std::min(C2*std::abs(qc), std::abs(qd)));// (CS eq 22)
        }

        // Check if 2nd derivative is close to roundoff error
        qa = std::max(std::abs(q_jm1(n,i)), std::abs(q_jm2(n,i)));
        qb = std::max(std::max(std::abs(q_j(n,i)),
                               std::abs(q_jp1(n,i))), std::abs(q_jp2(n,i)));

        Real rho = 0.0;
        if (std::abs(qd) > (1.0e-12)*std::max(qa, qb)) {
          // Limiter is not sensitive to roundoff. Use limited ratio (MC eq 27)
          rho = qe/qd;
        }

        Real tmp_m = q_j(n,i) - rho*dqf_minus(i);
        Real tmp_p = q_j(n,i) + rho*dqf_plus(i);
        Real tmp2_m = q_j(n,i) - 2.0*dqf_plus(i);
        Real tmp2_p = q_j(n,i) + 2.0*dqf_minus(i);

        // Check if relative change in limited 2nd deriv is > roundoff
        // Check for local extrema
        if ((qa_tmp <= 0.0 || qb_tmp <= 0.0)) {
          if (rho <= (1.0 - (1.0e-12))) {
            // Limit smooth extrema
            qminus(i) = tmp_m; // (CS eq 23)
            qplus(i) = tmp_p;
          }
          // No extrema detected
        } else {
          // Overshoot j-1/2,R / j,(-) state
          if (std::abs(dqf_minus(i)) >= 2.0*std::abs(dqf_plus(i))) {
            qminus(i) = tmp2_m;
          }
          // Overshoot j+1/2,L / j,(+) state
          if (std::abs(dqf_plus(i)) >= 2.0*std::abs(dqf_minus(i))) {
            qplus(i) = tmp2_p;
          }
        }
      }

      //--- Step 4b. -------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply Mignone limiters to parabolic
      // interpolant. Note, the Mignone limiter does not check for cell-averaged extrema:
    } else {
      for (int i=il; i<=iu; ++i) {
        Real qa = dqf_minus(i)*dqf_plus(i);
        if (qa <= 0.0) { // Local extrema detected
          qminus(i) = q_j(n,i);
          qplus(i) = q_j(n,i);
        } else { // No extrema detected
          // Overshoot j-1/2,R / j,(-) state
          if (std::abs(dqf_minus(i)) >= hplus_ratio_j(j)*std::abs(dqf_plus(i))) {
            qminus(i) = q_j(n,i) - hplus_ratio_j(j)*dqf_plus(i);
          }
          // Overshoot j+1/2,L / j,(+) state
          if (std::abs(dqf_plus(i)) >= hminus_ratio_j(j)*std::abs(dqf_minus(i))) {
            qplus(i) = q_j(n,i) + hminus_ratio_j(j)*dqf_minus(i);
          }
        }
      }
    }
    //--- Step 5. ----------------------------------------------------------------------
    // Convert limited cell-centered values to interface-centered L/R Riemann states
    // both L/R values defined over [il,iu]
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      ql_jph(n,i ) = qplus(i);
      qr_jmh(n,i ) = qminus(i);
    }
  } // end char PPM loop over =nu


  // compute ql_(j+1/2) and qr_(j-1/2)
  for (int n=0; n<=nu; ++n) {
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      ql(n,i) = ql_jph(n,i);
      qr(n,i) = qr_jmh(n,i);
    }
  }
  return;
}



void Reconstruction::PiecewiseParabolicX2(
    const int k, const int j, const int il, const int iu,
    const AthenaArray<Real> &q, const int array_order,
    AthenaArray<Real> &ql, AthenaArray<Real> &qr) {

  if(array_order < 0){

    const int nu = q.GetDim1() - 1;
  // CS08 constant used in second derivative limiter, >1 , independent of h
    const Real C2 = 1.25;

  // set work arrays used for primitive/characterstic cell-averages to scratch
    AthenaArray<Real> &q_jm2 = scr2_in_, &q_jm1 = scr3_in_,
                    &q_j = scr4_in_, &q_jp1 = scr5_in_, &q_jp2 = scr6_in_,
                    &qr_jmh = scr7_in_, &ql_jph = scr8_in_;

  // set work PPM work arrays to shallow copies of scratch arrays:
    AthenaArray<Real> &dd = scr02_i_, &dd_jm1 = scr03_i_, &dd_jp1 = scr04_i_,
                   &dph = scr05_i_, &dph_jp1 = scr06_i_;

    AthenaArray<Real> &d2qc_jm1 = scr07_i_, &d2qc = scr08_i_, &d2qc_jp1 = scr09_i_,
                        &d2qf = scr10_i_;

    AthenaArray<Real> &qplus = scr11_i_, &qminus = scr12_i_, &dqf_plus = scr13_i_,
                &dqf_minus = scr14_i_;

  // cache the x1-sliced primitive states for eigensystem calculation
    for (int i=il; i<=iu; ++i) {
#pragma omp simd
      for (int n=0; n<=nu; ++n) {
        q_j  (n,i) = q(k,j  ,i,n);
        q_jm2(n,i) = q(k,j-2,i,n);
        q_jm1(n,i) = q(k,j-1,i,n);
        q_jp1(n,i) = q(k,j+1,i,n);
        q_jp2(n,i) = q(k,j+2,i,n);
      }
    }

  //--- Step 1. ------------------------------------------------------------------------
  // Reconstruct interface averages <a>_{j-1/2} and <a>_{j+1/2}
    for (int n=0; n<=nu; ++n) {
    // Compute average slope in j-1, j, j+1 zones
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
      // nonuniform or uniform Cartesian-like coord reconstruction from volume averages:
        if (!curvilinear[X2DIR]) {
          Real qa = (q_j(n,i) - q_jm1(n,i));
          Real qb = (q_jp1(n,i) - q_j(n,i));
          dd_jm1(i) = c1j(j-1)*qa + c2j(j-1)*(q_jm1(n,i) - q_jm2(n,i));
          dd    (i) = c1j(j  )*qb + c2j(j  )*qa;
          dd_jp1(i) = c1j(j+1)*(q_jp2(n,i) - q_jp1(n,i)) + c2j(j+1)*qb;

        // Approximate interface average at j-1/2 and j+1/2 using PPM (CW eq 1.6)
        // KGF: group the biased stencil quantities to preserve FP symmetry
          dph(i) = (c3j(j)*q_jm1(n,i) + c4j(j)*q_j(n,i)) +
                   (c5j(j)*dd_jm1(i) + c6j(j)*dd(i));
          dph_jp1(i) = (c3j(j+1)*q_j(n,i) + c4j(j+1)*q_jp1(n,i)) +
                     (c5j(j+1)*dd(i) + c6j(j+1)*dd_jp1(i));
        } else { // spherical-polar coordinates meridional x2 direction:
        // TODO(felker): implement weight calculation in reconstruction.cpp
          dph(i) = c1j(j)*q_jm2(n,i) + c2j(j)*q_jm1(n,i) + c3j(j)*q_j(n,i)
                   + c4j(j)*q_jp1(n,i);
          dph_jp1(i) = c1j(j+1)*q_jm1(n,i) + c2j(j+1)*q_j(n,i) + c3j(j+1)*q_jp1(n,i)
                     + c4j(j+1)*q_jp2(n,i);
        }
      }

    //--- Step 2a. ---------------------------------------------------------------------
    // Uniform Cartesian-like coordinate: limit interpolated interface states (CD 4.3.1)
      if (uniform[X2DIR] && !curvilinear[X2DIR]) {
      // approximate second derivative at interfaces for smooth extrema preservation
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
          d2qc_jm1(i) = q_jm2(n,i) + q_j  (n,i) - 2.0*q_jm1(n,i);
          d2qc    (i) = q_jm1(n,i) + q_jp1(n,i) - 2.0*q_j  (n,i); //(CD eq 85a) (no 1/2)
          d2qc_jp1(i) = q_j  (n,i) + q_jp2(n,i) - 2.0*q_jp1(n,i);
        }

      // j-1/2
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
          Real qa_tmp = dph(i) - q_jm1(n,i); // (CD eq 84a)
          Real qb_tmp = q_j(n,i) - dph(i);     // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
          Real qa = 3.0*(q_jm1(n,i) + q_j(n,i) - 2.0*dph(i));  // (CD eq 85b)
          Real qb = d2qc_jm1(i);    // (CD eq 85a) (no 1/2)
          Real qc = d2qc(i);   // (CD eq 85c) (no 1/2)
          Real qd = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
            qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
          }
          Real dph_tmp = 0.5*(q_jm1(n,i) + q_j(n,i)) - qd/6.0;
          if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at j-1/2 face
            dph(i) = dph_tmp;
          }
        }
      // j+1/2
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
          Real qa_tmp = dph_jp1(i) - q_j(n,i);       // (CD eq 84a)
          Real qb_tmp = q_jp1(n,i) - dph_jp1(i);   // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
          Real qa = 3.0*(q_j(n,i) + q_jp1(n,i)  - 2.0*dph_jp1(i));  // (CD eq 85b)
          Real qb = d2qc(i);            // (CD eq 85a) (no 1/2)
          Real qc = d2qc_jp1(i);   // (CD eq 85c) (no 1/2)
          Real qd = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
            qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
          }
          Real dphjp1_tmp = 0.5*(q_j(n,i) + q_jp1(n,i)) - qd/6.0;
          if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at j+1/2 face
            dph_jp1(i) = dphjp1_tmp;
          }
        }

#pragma omp simd
        for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
          d2qf(i) = 6.0*(dph(i) + dph_jp1(i) - 2.0*q_j(n,i)); // a6 coefficient * -2
        }

      //--- Step 2b. -------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply strict monotonicity constraints
      // (Mignone eq 45)
      } else {
#pragma omp simd
        for (int i=il; i<=iu; ++i) {
          dph    (i) = std::min(dph    (i), std::max(q_j(n,i),q_jm1(n,i)));
          dph_jp1(i) = std::min(dph_jp1(i), std::max(q_j(n,i),q_jp1(n,i)));

          dph    (i) = std::max(dph    (i), std::min(q_j(n,i),q_jm1(n,i)));
          dph_jp1(i) = std::max(dph_jp1(i), std::min(q_j(n,i),q_jp1(n,i)));
        }
      }

    // Cache Riemann states for both non-/uniform limiters
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        qminus(i) = dph(i  );
        qplus(i) =  dph_jp1(i );
      }

    //--- Step 3. ----------------------------------------------------------------------
    // Compute cell-centered difference stencils (MC section 2.4.1)
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        dqf_minus(i) = q_j(n,i) - qminus(i); // (CS eq 25)
        dqf_plus(i)  = qplus(i) - q_j(n,i);
      }

    //--- Step 4a. ---------------------------------------------------------------------
    // For uniform Cartesian-like coordinate: apply CS limiters to parabolic interpolant
     if (uniform[X2DIR] && !curvilinear[X2DIR]) {
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
          Real qa_tmp = dqf_minus(i)*dqf_plus(i);
          Real qb_tmp = (q_jp1(n,i) - q_j(n,i))*(q_j(n,i) - q_jm1(n,i));

          Real qa = d2qc_jm1(i);
          Real qb = d2qc(i);
          Real qc = d2qc_jp1(i);
          Real qd = d2qf(i);
          Real qe = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc) && SIGN(qa) == SIGN(qd)) {
            // Extrema is smooth
            qe = SIGN(qd)* std::min(std::min(C2*std::abs(qa), C2*std::abs(qb)),
                                  std::min(C2*std::abs(qc), std::abs(qd)));// (CS eq 22)
          }

        // Check if 2nd derivative is close to roundoff error
          qa = std::max(std::abs(q_jm1(n,i)), std::abs(q_jm2(n,i)));
          qb = std::max(std::max(std::abs(q_j(n,i)),
                               std::abs(q_jp1(n,i))), std::abs(q_jp2(n,i)));

          Real rho = 0.0;
          if (std::abs(qd) > (1.0e-12)*std::max(qa, qb)) {
          // Limiter is not sensitive to roundoff. Use limited ratio (MC eq 27)
            rho = qe/qd;
          }

          Real tmp_m = q_j(n,i) - rho*dqf_minus(i);
          Real tmp_p = q_j(n,i) + rho*dqf_plus(i);
          Real tmp2_m = q_j(n,i) - 2.0*dqf_plus(i);
          Real tmp2_p = q_j(n,i) + 2.0*dqf_minus(i);

        // Check if relative change in limited 2nd deriv is > roundoff
        // Check for local extrema
          if ((qa_tmp <= 0.0 || qb_tmp <= 0.0)) {
            if (rho <= (1.0 - (1.0e-12))) {
            // Limit smooth extrema
              qminus(i) = tmp_m; // (CS eq 23)
              qplus(i) = tmp_p;
            }  
          // No extrema detected
          } else {
           // Overshoot j-1/2,R / j,(-) state
            if (std::abs(dqf_minus(i)) >= 2.0*std::abs(dqf_plus(i))) {
              qminus(i) = tmp2_m;
            }
          // Overshoot j+1/2,L / j,(+) state
            if (std::abs(dqf_plus(i)) >= 2.0*std::abs(dqf_minus(i))) {
              qplus(i) = tmp2_p;
            }
          }
        }

      //--- Step 4b. -------------------------------------------------------------------
      // Nonuniform and/or curvilinear coordinate: apply Mignone limiters to parabolic
      // interpolant. Note, the Mignone limiter does not check for cell-averaged extrema:
      } else {
        for (int i=il; i<=iu; ++i) {
          Real qa = dqf_minus(i)*dqf_plus(i);
          if (qa <= 0.0) { // Local extrema detected
            qminus(i) = q_j(n,i);
            qplus(i) = q_j(n,i);
          } else { // No extrema detected
          // Overshoot j-1/2,R / j,(-) state
            if (std::abs(dqf_minus(i)) >= hplus_ratio_j(j)*std::abs(dqf_plus(i))) {
              qminus(i) = q_j(n,i) - hplus_ratio_j(j)*dqf_plus(i);
            }
          // Overshoot j+1/2,L / j,(+) state
            if (std::abs(dqf_plus(i)) >= hminus_ratio_j(j)*std::abs(dqf_minus(i))) {
              qplus(i) = q_j(n,i) + hminus_ratio_j(j)*dqf_minus(i);
            }
          }
        }
      }
    //--- Step 5. ----------------------------------------------------------------------
    // Convert limited cell-centered values to interface-centered L/R Riemann states
    // both L/R values defined over [il,iu]
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        ql_jph(n,i ) = qplus(i);
        qr_jmh(n,i ) = qminus(i);
      }
    } // end char PPM loop over =nu


  // compute ql_(j+1/2) and qr_(j-1/2)
    for (int i=il; i<=iu; ++i) {
#pragma omp simd
      for (int n=0; n<=nu; ++n) {
        ql(i,n) = ql_jph(n,i);
        qr(i,n) = qr_jmh(n,i);
      }
    }

  }// end array ordering
  return;
}

//----------------------------------------------------------------------------------------
//! \fn Reconstruction::PiecewiseParabolicX3()
//  \brief Returns L/R interface values in X3-dir constructed using fourth-order PPM and
//         Colella-Sekora or Mignone limiting over [kl,ku][jl,ju][il,iu]

void Reconstruction::PiecewiseParabolicX3(
    const int k, const int j, const int il, const int iu,
    const AthenaArray<Real> &q,
    AthenaArray<Real> &ql, AthenaArray<Real> &qr) {
  const int nu = q.GetDim4() - 1;
  // CS08 constant used in second derivative limiter, >1 , independent of h
  const Real C2 = 1.25;

  // set work arrays used for primitive/characterstic cell-averages to scratch
  AthenaArray<Real> &q_km2 = scr2_ni_, &q_km1 = scr3_ni_,
                     &q_k = scr4_ni_, &q_kp1 = scr5_ni_, &q_kp2 = scr6_ni_,
                &qr_kmh = scr7_ni_, &ql_kph = scr8_ni_;

  // set work PPM work arrays to shallow copies of scratch arrays:
  AthenaArray<Real> &dd = scr02_i_, &dd_km1 = scr03_i_, &dd_kp1 = scr04_i_,
                   &dph = scr05_i_, &dph_kp1 = scr06_i_;

  AthenaArray<Real> &d2qc_km1 = scr07_i_, &d2qc = scr08_i_,
                    &d2qc_kp1 = scr09_i_, &d2qf = scr10_i_;

  AthenaArray<Real> &qplus = scr11_i_, &qminus = scr12_i_, &dqf_plus = scr13_i_,
                &dqf_minus = scr14_i_;

  // cache the x1-sliced primitive states for eigensystem calculation
  for (int n=0; n<=nu; ++n) {
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      q_k  (n,i) = q(n,k  ,j,i);
      q_km2(n,i) = q(n,k-2,j,i);
      q_km1(n,i) = q(n,k-1,j,i);
      q_kp1(n,i) = q(n,k+1,j,i);
      q_kp2(n,i) = q(n,k+2,j,i);
    }
  }

  //--- Step 1. -------------------------------------------------------------------------
  // Reconstruct interface averages <a>_{k-1/2} and <a>_{k+1/2}
  for (int n=0; n<=nu; ++n) {
    // Compute average slope in k-1, k, k+1 zones
#pragma omp simd simdlen(SIMD_WIDTH)
    for (int i=il; i<=iu; ++i) {
      // nonuniform or uniform Cartesian-like coord reconstruction from volume averages:
      Real qa = (q_k(n,i) - q_km1(n,i));
      Real qb = (q_kp1(n,i) - q_k(n,i));
      dd_km1(i) = c1k(k-1)*qa + c2k(k-1)*(q_km1(n,i) - q_km2(n,i));
      dd    (i) = c1k(k  )*qb + c2k(k  )*qa;
      dd_kp1(i) = c1k(k+1)*(q_kp2(n,i) - q_kp1(n,i)) + c2k(k+1)*qb;

      // Approximate interface average at k-1/2 and k+1/2 using PPM (CW eq 1.6)
      // KGF: group the biased stencil quantities to preserve FP symmetry
      dph(i) = (c3k(k)*q_km1(n,i) + c4k(k)*q_k(n,i)) +
               (c5k(k)*dd_km1(i) + c6k(k)*dd(i));
      dph_kp1(i) = (c3k(k+1)*q_k(n,i) + c4k(k+1)*q_kp1(n,i)) +
                   (c5k(k+1)*dd(i) + c6k(k+1)*dd_kp1(i));
    }

    //--- Step 2a. -----------------------------------------------------------------------
    // Uniform Cartesian-like coordinate: limit interpolated interface states (CD 4.3.1)
    if (uniform[X3DIR]) {
      // approximate second derivative at interfaces for smooth extrema preservation
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
        d2qc_km1(i) = q_km2(n,i) + q_k  (n,i) - 2.0*q_km1(n,i);
        d2qc    (i) = q_km1(n,i) + q_kp1(n,i) - 2.0*q_k  (n,i); // (CD eq 85a) (no 1/2)
        d2qc_kp1(i) = q_k  (n,i) + q_kp2(n,i) - 2.0*q_kp1(n,i);
      }

      // k-1/2
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        Real qa_tmp = dph(i) - q_km1(n,i); // (CD eq 84a)
        Real qb_tmp = q_k(n,i) - dph(i);     // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
        Real qa = 3.0*(q_km1(n,i) + q_k(n,i) - 2.0*dph(i));  // (CD eq 85b)
        Real qb = d2qc_km1(i);    // (CD eq 85a) (no 1/2)
        Real qc = d2qc(i);   // (CD eq 85c) (no 1/2)
        Real qd = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
          qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
        }
        Real dph_tmp = 0.5*(q_km1(n,i) + q_k(n,i)) - qd/6.0;
        if (qa_tmp*qb_tmp < 0.0) {  // Local extrema detected at k-1/2 face
          dph(i) = dph_tmp;
        }
      }
      // k+1/2
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        Real qa_tmp = dph_kp1(i) - q_k(n,i);       // (CD eq 84a)
        Real qb_tmp = q_kp1(n,i) - dph_kp1(i);   // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
        Real qa = 3.0*(q_k(n,i) + q_kp1(n,i) - 2.0*dph_kp1(i));  // (CD eq 85b)
        Real qb = d2qc(i);            // (CD eq 85a) (no 1/2)
        Real qc = d2qc_kp1(i);   // (CD eq 85c) (no 1/2)
        Real qd = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
          qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
        }
        Real dphkp1_tmp = 0.5*(q_k(n,i) + q_kp1(n,i)) - qd/6.0;
        if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at k+1/2 face
          dph_kp1(i) = dphkp1_tmp;
        }
      }

#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
        d2qf(i) = 6.0*(dph(i) + dph_kp1(i) - 2.0*q_k(n,i)); // a6 coefficient * -2
      }
      //--- Step 2b. ---------------------------------------------------------------------
      // Nonuniform coordinate spacing: apply strict monotonicity constraints
      // (Mignone eq 45)
    } else {
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        dph    (i) = std::min(dph    (i), std::max(q_k(n,i),q_km1(n,i)));
        dph_kp1(i) = std::min(dph_kp1(i), std::max(q_k(n,i),q_kp1(n,i)));

        dph    (i) = std::max(dph    (i), std::min(q_k(n,i),q_km1(n,i)));
        dph_kp1(i) = std::max(dph_kp1(i), std::min(q_k(n,i),q_kp1(n,i)));
      }
    }

    // Cache Riemann states for both non-/uniform limiters
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      qminus(i) = dph(i  );
      qplus(i) =  dph_kp1(i );
    }

    //--- Step 3. --------------------------------------------------------------------
    // Compute cell-centered difference stencils (MC section 2.4.1)
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      dqf_minus(i) = q_k(n,i) - qminus(i); // (CS eq 25)
      dqf_plus(i)  = qplus(i) - q_k(n,i);
    }

    //--- Step 4a. -----------------------------------------------------------------------
    // For uniform Cartesian-like coordinate: apply CS limiters to parabolic interpolant
    if (uniform[X3DIR]) {
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
        Real qa_tmp = dqf_minus(i)*dqf_plus(i);
        Real qb_tmp = (q_kp1(n,i) - q_k(n,i))*(q_k(n,i) - q_km1(n,i));

        // Check if extrema is smooth
        Real qa = d2qc_km1(i);
        Real qb = d2qc(i);
        Real qc = d2qc_kp1(i);
        Real qd = d2qf(i);
        Real qe = 0.0;
        if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc) && SIGN(qa) == SIGN(qd)) {
          // Extrema is smooth
          qe = SIGN(qd)* std::min(std::min(C2*std::abs(qa),C2*std::abs(qb)),
                                  std::min(C2*std::abs(qc),std::abs(qd))); // (CS eq 22)
        }

        // Check if 2nd derivative is close to roundoff error
        qa = std::max(std::abs(q_km1(n,i)),std::abs(q_km2(n,i)));
        qb = std::max(std::max(std::abs(q_k(n,i)),
                               std::abs(q_kp1(n,i))), std::abs(q_kp2(n,i)));

        Real rho = 0.0;
        if (std::abs(qd) > (1.0e-12)*std::max(qa,qb)) {
          // Limiter is not sensitive to roundoff. Use limited ratio (MC eq 27)
          rho = qe/qd;
        }

        Real tmp_m = q_k(n,i) - rho*dqf_minus(i);
        Real tmp_p = q_k(n,i) + rho*dqf_plus(i);
        Real tmp2_m = q_k(n,i) - 2.0*dqf_plus(i);
        Real tmp2_p = q_k(n,i) + 2.0*dqf_minus(i);

        // Check for local extrema
        if (qa_tmp <= 0.0 || qb_tmp <= 0.0 ) {
          // Check if relative change in limited 2nd deriv is > roundoff
          if (rho <= (1.0 - (1.0e-12))) {
            // Limit smooth extrema
            qminus(i) = tmp_m; // (CS eq 23)
            qplus(i) = tmp_p;
          }

          // No extrema detected
        } else {
          // Overshoot k-1/2,R / k,(-) state
          if (std::abs(dqf_minus(i)) >= 2.0*std::abs(dqf_plus(i))) {
            qminus(i) = tmp2_m;
          }
          // Overshoot k+1/2,L / k,(+) state
          if (std::abs(dqf_plus(i)) >= 2.0*std::abs(dqf_minus(i))) {
            qplus(i) = tmp2_p;
          }
        }
      }
      //--- Step 4b. ---------------------------------------------------------------------
      // Nonuniform coordinate spacing: apply Mignone limiters to parabolic interpolant
      // Note, the Mignone limiter does not check for cell-averaged extrema:
    } else {
      for (int i=il; i<=iu; ++i) {
        Real qa = dqf_minus(i)*dqf_plus(i);
        if (qa <= 0.0) { // Local extrema detected
          qminus(i) = q_k(n,i);
          qplus(i) = q_k(n,i);
        } else { // No extrema detected
          // TODO(felker): could delete hplus_ratio_k() arrays for curvilinear PPMx3
          // Overshoot k-1/2,R / k,(-) state
          if (std::abs(dqf_minus(i)) >= hplus_ratio_k(k)*std::abs(dqf_plus(i))) {
            qminus(i) = q_k(n,i) - hplus_ratio_k(k)*dqf_plus(i);
          }
          // Overshoot k+1/2,L / k,(+) state
          if (std::abs(dqf_plus(i)) >= hminus_ratio_k(k)*std::abs(dqf_minus(i))) {
            qplus(i) = q_k(n,i) + hminus_ratio_k(k)*dqf_minus(i);
          }
        }
      }
    }
    //--- Step 5. --------------------------------------------------------------------
    // Convert limited cell-centered values to interface-centered L/R Riemann states
    // both L/R values defined over [il,iu]
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      ql_kph(n,i ) = qplus(i);
      qr_kmh(n,i ) = qminus(i);
    }
  } // end char PPM loop over =nu

  // compute ql_(k+1/2) and qr_(k-1/2)
  for (int n=0; n<=nu; ++n) {
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      ql(n,i) = ql_kph(n,i);
      qr(n,i) = qr_kmh(n,i);
    }
  }
  return;
}



void Reconstruction::PiecewiseParabolicX3(
    const int k, const int j, const int il, const int iu,
    const AthenaArray<Real> &q, const int array_order,
    AthenaArray<Real> &ql, AthenaArray<Real> &qr) {

  if(array_order < 0){

    const int nu = q.GetDim1() - 1;
  // CS08 constant used in second derivative limiter, >1 , independent of h
    const Real C2 = 1.25;

  // set work arrays used for primitive/characterstic cell-averages to scratch
    AthenaArray<Real> &q_km2 = scr2_in_, &q_km1 = scr3_in_,
                    &q_k = scr4_in_, &q_kp1 = scr5_in_, &q_kp2 = scr6_in_,
                    &qr_kmh = scr7_in_, &ql_kph = scr8_in_;

  // set work PPM work arrays to shallow copies of scratch arrays:
    AthenaArray<Real> &dd = scr02_i_, &dd_km1 = scr03_i_, &dd_kp1 = scr04_i_,
                   &dph = scr05_i_, &dph_kp1 = scr06_i_;

    AthenaArray<Real> &d2qc_km1 = scr07_i_, &d2qc = scr08_i_,
                    &d2qc_kp1 = scr09_i_, &d2qf = scr10_i_;

    AthenaArray<Real> &qplus = scr11_i_, &qminus = scr12_i_, &dqf_plus = scr13_i_,
                &dqf_minus = scr14_i_;

  // cache the x1-sliced primitive states for eigensystem calculation
    for (int i=il; i<=iu; ++i) {
//#pragma omp simd
      for (int n=0; n<=nu; ++n) {
        q_k  (n,i) = q(k  ,j,i,n);
        q_km2(n,i) = q(k-2,j,i,n);
        q_km1(n,i) = q(k-1,j,i,n);
        q_kp1(n,i) = q(k+1,j,i,n);
        q_kp2(n,i) = q(k+2,j,i,n);
      }
    }

  //--- Step 1. -------------------------------------------------------------------------
  // Reconstruct interface averages <a>_{k-1/2} and <a>_{k+1/2}
    for (int n=0; n<=nu; ++n) {
    // Compute average slope in k-1, k, k+1 zones
#pragma omp simd simdlen(SIMD_WIDTH)
      for (int i=il; i<=iu; ++i) {
      // nonuniform or uniform Cartesian-like coord reconstruction from volume averages:
        Real qa = (q_k(n,i) - q_km1(n,i));
        Real qb = (q_kp1(n,i) - q_k(n,i));
        dd_km1(i) = c1k(k-1)*qa + c2k(k-1)*(q_km1(n,i) - q_km2(n,i));
        dd    (i) = c1k(k  )*qb + c2k(k  )*qa;
        dd_kp1(i) = c1k(k+1)*(q_kp2(n,i) - q_kp1(n,i)) + c2k(k+1)*qb;

      // Approximate interface average at k-1/2 and k+1/2 using PPM (CW eq 1.6)
      // KGF: group the biased stencil quantities to preserve FP symmetry
        dph(i) = (c3k(k)*q_km1(n,i) + c4k(k)*q_k(n,i)) +
               (c5k(k)*dd_km1(i) + c6k(k)*dd(i));
        dph_kp1(i) = (c3k(k+1)*q_k(n,i) + c4k(k+1)*q_kp1(n,i)) +
                   (c5k(k+1)*dd(i) + c6k(k+1)*dd_kp1(i));
      }

    //--- Step 2a. -----------------------------------------------------------------------
    // Uniform Cartesian-like coordinate: limit interpolated interface states (CD 4.3.1)
      if (uniform[X3DIR]) {
      // approximate second derivative at interfaces for smooth extrema preservation
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
          d2qc_km1(i) = q_km2(n,i) + q_k  (n,i) - 2.0*q_km1(n,i);
          d2qc    (i) = q_km1(n,i) + q_kp1(n,i) - 2.0*q_k  (n,i); // (CD eq 85a) (no 1/2)
          d2qc_kp1(i) = q_k  (n,i) + q_kp2(n,i) - 2.0*q_kp1(n,i);
        }

      // k-1/2
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
          Real qa_tmp = dph(i) - q_km1(n,i); // (CD eq 84a)
          Real qb_tmp = q_k(n,i) - dph(i);     // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
          Real qa = 3.0*(q_km1(n,i) + q_k(n,i) - 2.0*dph(i));  // (CD eq 85b)
          Real qb = d2qc_km1(i);    // (CD eq 85a) (no 1/2)
          Real qc = d2qc(i);   // (CD eq 85c) (no 1/2)
          Real qd = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
            qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
          }
          Real dph_tmp = 0.5*(q_km1(n,i) + q_k(n,i)) - qd/6.0;
          if (qa_tmp*qb_tmp < 0.0) {  // Local extrema detected at k-1/2 face
            dph(i) = dph_tmp;
          }
        }
      // k+1/2
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
          Real qa_tmp = dph_kp1(i) - q_k(n,i);       // (CD eq 84a)
          Real qb_tmp = q_kp1(n,i) - dph_kp1(i);   // (CD eq 84b)
        // KGF: add the off-centered quantities first to preserve FP symmetry
          Real qa = 3.0*(q_k(n,i) + q_kp1(n,i) - 2.0*dph_kp1(i));  // (CD eq 85b)
          Real qb = d2qc(i);            // (CD eq 85a) (no 1/2)
          Real qc = d2qc_kp1(i);   // (CD eq 85c) (no 1/2)
          Real qd = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc)) {
            qd = SIGN(qa)* std::min(C2*std::abs(qb),
                                  std::min(C2*std::abs(qc), std::abs(qa)));
          }
          Real dphkp1_tmp = 0.5*(q_k(n,i) + q_kp1(n,i)) - qd/6.0;
          if (qa_tmp*qb_tmp < 0.0) { // Local extrema detected at k+1/2 face
            dph_kp1(i) = dphkp1_tmp;
          }
        }

#pragma omp simd
        for (int i=il; i<=iu; ++i) {
        // KGF: add the off-centered quantities first to preserve FP symmetry
          d2qf(i) = 6.0*(dph(i) + dph_kp1(i) - 2.0*q_k(n,i)); // a6 coefficient * -2
        }
      //--- Step 2b. ---------------------------------------------------------------------
      // Nonuniform coordinate spacing: apply strict monotonicity constraints
      // (Mignone eq 45)
      } else {
#pragma omp simd
        for (int i=il; i<=iu; ++i) {
          dph    (i) = std::min(dph    (i), std::max(q_k(n,i),q_km1(n,i)));
          dph_kp1(i) = std::min(dph_kp1(i), std::max(q_k(n,i),q_kp1(n,i)));

          dph    (i) = std::max(dph    (i), std::min(q_k(n,i),q_km1(n,i)));
          dph_kp1(i) = std::max(dph_kp1(i), std::min(q_k(n,i),q_kp1(n,i)));
        }
      }

    // Cache Riemann states for both non-/uniform limiters
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        qminus(i) = dph(i  );
        qplus(i) =  dph_kp1(i );
      }

    //--- Step 3. --------------------------------------------------------------------
    // Compute cell-centered difference stencils (MC section 2.4.1)
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        dqf_minus(i) = q_k(n,i) - qminus(i); // (CS eq 25)
        dqf_plus(i)  = qplus(i) - q_k(n,i);
      }

    //--- Step 4a. -----------------------------------------------------------------------
    // For uniform Cartesian-like coordinate: apply CS limiters to parabolic interpolant
      if (uniform[X3DIR]) {
#pragma omp simd simdlen(SIMD_WIDTH)
        for (int i=il; i<=iu; ++i) {
          Real qa_tmp = dqf_minus(i)*dqf_plus(i);
          Real qb_tmp = (q_kp1(n,i) - q_k(n,i))*(q_k(n,i) - q_km1(n,i));

        // Check if extrema is smooth
          Real qa = d2qc_km1(i);
          Real qb = d2qc(i);
          Real qc = d2qc_kp1(i);
          Real qd = d2qf(i);
          Real qe = 0.0;
          if (SIGN(qa) == SIGN(qb) && SIGN(qa) == SIGN(qc) && SIGN(qa) == SIGN(qd)) {
          // Extrema is smooth
            qe = SIGN(qd)* std::min(std::min(C2*std::abs(qa),C2*std::abs(qb)),
                                  std::min(C2*std::abs(qc),std::abs(qd))); // (CS eq 22)
          }

        // Check if 2nd derivative is close to roundoff error
          qa = std::max(std::abs(q_km1(n,i)),std::abs(q_km2(n,i)));
          qb = std::max(std::max(std::abs(q_k(n,i)),
                               std::abs(q_kp1(n,i))), std::abs(q_kp2(n,i)));

          Real rho = 0.0;
          if (std::abs(qd) > (1.0e-12)*std::max(qa,qb)) {
          // Limiter is not sensitive to roundoff. Use limited ratio (MC eq 27)
            rho = qe/qd;
          }

          Real tmp_m = q_k(n,i) - rho*dqf_minus(i);
          Real tmp_p = q_k(n,i) + rho*dqf_plus(i);
          Real tmp2_m = q_k(n,i) - 2.0*dqf_plus(i);
          Real tmp2_p = q_k(n,i) + 2.0*dqf_minus(i);

        // Check for local extrema
          if (qa_tmp <= 0.0 || qb_tmp <= 0.0 ) {
          // Check if relative change in limited 2nd deriv is > roundoff
            if (rho <= (1.0 - (1.0e-12))) {
            // Limit smooth extrema
              qminus(i) = tmp_m; // (CS eq 23)
              qplus(i) = tmp_p;
            }

          // No extrema detected
          } else {
          // Overshoot k-1/2,R / k,(-) state
            if (std::abs(dqf_minus(i)) >= 2.0*std::abs(dqf_plus(i))) {
              qminus(i) = tmp2_m;
            }
          // Overshoot k+1/2,L / k,(+) state
            if (std::abs(dqf_plus(i)) >= 2.0*std::abs(dqf_minus(i))) {
              qplus(i) = tmp2_p;
            }
          }
        }
      //--- Step 4b. ---------------------------------------------------------------------
      // Nonuniform coordinate spacing: apply Mignone limiters to parabolic interpolant
      // Note, the Mignone limiter does not check for cell-averaged extrema:
      } else {
        for (int i=il; i<=iu; ++i) {
          Real qa = dqf_minus(i)*dqf_plus(i);
          if (qa <= 0.0) { // Local extrema detected
            qminus(i) = q_k(n,i);
            qplus(i) = q_k(n,i);
          } else { // No extrema detected
          // TODO(felker): could delete hplus_ratio_k() arrays for curvilinear PPMx3
          // Overshoot k-1/2,R / k,(-) state
            if (std::abs(dqf_minus(i)) >= hplus_ratio_k(k)*std::abs(dqf_plus(i))) {
              qminus(i) = q_k(n,i) - hplus_ratio_k(k)*dqf_plus(i);
            }
          // Overshoot k+1/2,L / k,(+) state
            if (std::abs(dqf_plus(i)) >= hminus_ratio_k(k)*std::abs(dqf_minus(i))) {
              qplus(i) = q_k(n,i) + hminus_ratio_k(k)*dqf_minus(i);
            }
          }
        }
      }
    //--- Step 5. --------------------------------------------------------------------
    // Convert limited cell-centered values to interface-centered L/R Riemann states
    // both L/R values defined over [il,iu]
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        ql_kph(n,i ) = qplus(i);
        qr_kmh(n,i ) = qminus(i);
      }
    } // end char PPM loop over =nu

  // compute ql_(k+1/2) and qr_(k-1/2)
    for (int i=il; i<=iu; ++i) {
#pragma omp simd
      for (int n=0; n<=nu; ++n) {
        ql(i,n) = ql_kph(n,i);
        qr(i,n) = qr_kmh(n,i);
      }
    }

  }// end array_order
  return;
}
