#pragma OPENCL EXTENSION cl_khr_fp64: enable
#define PROBLEM_SIZE   64
#define IMAX      PROBLEM_SIZE
#define JMAX      PROBLEM_SIZE
#define KMAX      PROBLEM_SIZE
#define IMAXP     IMAX/2*2
#define JMAXP     JMAX/2*2

//--------------------------------------------------------------------------
// from set_constants()
//--------------------------------------------------------------------------
#define c1        1.4
#define c2        0.4
#define c3        0.1
#define c4        1.0
#define c5        1.4

#define dnxm1     (1.0 / (double)(PROBLEM_SIZE-1))
#define dnym1     (1.0 / (double)(PROBLEM_SIZE-1))
#define dnzm1     (1.0 / (double)(PROBLEM_SIZE-1))

#define c1c2      (c1 * c2)
#define c1c5      (c1 * c5)
#define c3c4      (c3 * c4)
#define c1345     (c1c5 * c3c4)

#define conz1     (1.0-c1c5)

#define tx1       (1.0 / (dnxm1 * dnxm1))
#define tx2       (1.0 / (2.0 * dnxm1))
#define tx3       (1.0 / dnxm1)

#define ty1       (1.0 / (dnym1 * dnym1))
#define ty2       (1.0 / (2.0 * dnym1))
#define ty3       (1.0 / dnym1)

#define tz1       (1.0 / (dnzm1 * dnzm1))
#define tz2       (1.0 / (2.0 * dnzm1))
#define tz3       (1.0 / dnzm1)

#define dx1       0.75
#define dx2       0.75
#define dx3       0.75
#define dx4       0.75
#define dx5       0.75

#define dy1       0.75
#define dy2       0.75
#define dy3       0.75
#define dy4       0.75
#define dy5       0.75

#define dz1       1.0
#define dz2       1.0
#define dz3       1.0
#define dz4       1.0
#define dz5       1.0

#define dxmax     max(dx3, dx4)
#define dymax     max(dy2, dy4)
#define dzmax     max(dz2, dz3)

#define dssp      (0.25 * max(dx1, max(dy1, dz1) ))

#define c4dssp    (4.0 * dssp)
#define c5dssp    (5.0 * dssp)

#define dt        DT_DEFAULT
#define dttx1     (dt*tx1)
#define dttx2     (dt*tx2)
#define dtty1     (dt*ty1)
#define dtty2     (dt*ty2)
#define dttz1     (dt*tz1)
#define dttz2     (dt*tz2)

#define c2dttx1   (2.0*dttx1)
#define c2dtty1   (2.0*dtty1)
#define c2dttz1   (2.0*dttz1)

#define dtdssp    (dt*dssp)

#define comz1     dtdssp
#define comz4     (4.0*dtdssp)
#define comz5     (5.0*dtdssp)
#define comz6     (6.0*dtdssp)

#define c3c4tx3   (c3c4*tx3)
#define c3c4ty3   (c3c4*ty3)
#define c3c4tz3   (c3c4*tz3)

#define dx1tx1    (dx1*tx1)
#define dx2tx1    (dx2*tx1)
#define dx3tx1    (dx3*tx1)
#define dx4tx1    (dx4*tx1)
#define dx5tx1    (dx5*tx1)

#define dy1ty1    (dy1*ty1)
#define dy2ty1    (dy2*ty1)
#define dy3ty1    (dy3*ty1)
#define dy4ty1    (dy4*ty1)
#define dy5ty1    (dy5*ty1)

#define dz1tz1    (dz1*tz1)
#define dz2tz1    (dz2*tz1)
#define dz3tz1    (dz3*tz1)
#define dz4tz1    (dz4*tz1)
#define dz5tz1    (dz5*tz1)

#define c2iv      2.5
#define con43     (4.0/3.0)
#define con16     (1.0/6.0)

#define xxcon1    (c3c4tx3*con43*tx3)
#define xxcon2    (c3c4tx3*tx3)
#define xxcon3    (c3c4tx3*conz1*tx3)
#define xxcon4    (c3c4tx3*con16*tx3)
#define xxcon5    (c3c4tx3*c1c5*tx3)

#define yycon1    (c3c4ty3*con43*ty3)
#define yycon2    (c3c4ty3*ty3)
#define yycon3    (c3c4ty3*conz1*ty3)
#define yycon4    (c3c4ty3*con16*ty3)
#define yycon5    (c3c4ty3*c1c5*ty3)

#define zzcon1    (c3c4tz3*con43*tz3)
#define zzcon2    (c3c4tz3*tz3)
#define zzcon3    (c3c4tz3*conz1*tz3)
#define zzcon4    (c3c4tz3*con16*tz3)
#define zzcon5    (c3c4tz3*c1c5*tz3)


//---------------------------------------------------------------------
// compute eta-direction fluxes
//---------------------------------------------------------------------
__kernel void compute_rhs4(__global const double *g_u,
                           __global const double *g_us,
                           __global const double *g_vs,
                           __global const double *g_ws,
                           __global const double *g_qs,
                           __global const double *g_rho_i,
                           __global const double *g_square,
                           __global double *g_rhs,
                           int gp0,
                           int gp1,
                           int gp2)
{
  int i, j, k, m;
  double p_rhs[5];
  double p_u[5], p_up1[5], p_up2[5], p_um2[5];
  double p_um1[5] = {0.0, 0.0, 0.0, 0.0, 0.0};
  double p_us, p_usp1, p_usm1;
  double p_vs, p_vsp1, p_vsm1;
  double p_ws, p_wsp1, p_wsm1;
  double p_qs, p_qsp1, p_qsm1;
  double p_rho_i, p_rho_ip1, p_rho_im1;
  double p_square, p_squarep1, p_squarem1;

  k = get_global_id(1) + 1;
  i = get_global_id(0) + 1;
  if (k > (gp2-2) || i > (gp0-2)) return;

  __global double (*us)[JMAXP+1][IMAXP+1] =
    (__global double (*)[JMAXP+1][IMAXP+1])g_us;
  __global double (*vs)[JMAXP+1][IMAXP+1] =
    (__global double (*)[JMAXP+1][IMAXP+1])g_vs;
  __global double (*ws)[JMAXP+1][IMAXP+1] =
    (__global double (*)[JMAXP+1][IMAXP+1])g_ws;
  __global double (*qs)[JMAXP+1][IMAXP+1] =
    (__global double (*)[JMAXP+1][IMAXP+1])g_qs;
  __global double (*rho_i)[JMAXP+1][IMAXP+1] =
    (__global double (*)[JMAXP+1][IMAXP+1])g_rho_i;
  __global double (*square)[JMAXP+1][IMAXP+1] =
    (__global double (*)[JMAXP+1][IMAXP+1])g_square;
  __global double (*u)[JMAXP+1][IMAXP+1][5] =
    (__global double (*)[JMAXP+1][IMAXP+1][5])g_u;
  __global double (*rhs)[JMAXP+1][IMAXP+1][5] =
    (__global double (*)[JMAXP+1][IMAXP+1][5])g_rhs;

  p_u[0]   = u[k][0][i][0];
  p_u[1]   = u[k][0][i][1];
  p_u[2]   = u[k][0][i][2];
  p_u[3]   = u[k][0][i][3];
  p_u[4]   = u[k][0][i][4];
  p_up1[0] = u[k][1][i][0];
  p_up1[1] = u[k][1][i][1];
  p_up1[2] = u[k][1][i][2];
  p_up1[3] = u[k][1][i][3];
  p_up1[4] = u[k][1][i][4];
  p_up2[0] = u[k][2][i][0];
  p_up2[1] = u[k][2][i][1];
  p_up2[2] = u[k][2][i][2];
  p_up2[3] = u[k][2][i][3];
  p_up2[4] = u[k][2][i][4];

  p_us   = us[k][0][i];
  p_usp1 = us[k][1][i];
  p_vs   = vs[k][0][i];
  p_vsp1 = vs[k][1][i];
  p_ws   = ws[k][0][i];
  p_wsp1 = ws[k][1][i];
  p_qs   = qs[k][0][i];
  p_qsp1 = qs[k][1][i];
  p_rho_i   = rho_i[k][0][i];
  p_rho_ip1 = rho_i[k][1][i];
  p_square   = square[k][0][i];
  p_squarep1 = square[k][1][i];

#define LOOP_PROLOGUE_FULL         \
  p_um2[0] = p_um1[0];             \
  p_um2[1] = p_um1[1];             \
  p_um2[2] = p_um1[2];             \
  p_um2[3] = p_um1[3];             \
  p_um2[4] = p_um1[4];             \
  p_um1[0] = p_u[0];               \
  p_um1[1] = p_u[1];               \
  p_um1[2] = p_u[2];               \
  p_um1[3] = p_u[3];               \
  p_um1[4] = p_u[4];               \
  p_u[0]   = p_up1[0];             \
  p_u[1]   = p_up1[1];             \
  p_u[2]   = p_up1[2];             \
  p_u[3]   = p_up1[3];             \
  p_u[4]   = p_up1[4];             \
  p_up1[0] = p_up2[0];             \
  p_up1[1] = p_up2[1];             \
  p_up1[2] = p_up2[2];             \
  p_up1[3] = p_up2[3];             \
  p_up1[4] = p_up2[4];             \
  p_up2[0] = u[k][j+2][i][0];      \
  p_up2[1] = u[k][j+2][i][1];      \
  p_up2[2] = u[k][j+2][i][2];      \
  p_up2[3] = u[k][j+2][i][3];      \
  p_up2[4] = u[k][j+2][i][4];      \
                                   \
  p_usm1 = p_us;                   \
  p_us   = p_usp1;                 \
  p_usp1 = us[k][j+1][i];          \
  p_vsm1 = p_vs;                   \
  p_vs   = p_vsp1;                 \
  p_vsp1 = vs[k][j+1][i];          \
  p_wsm1 = p_ws;                   \
  p_ws   = p_wsp1;                 \
  p_wsp1 = ws[k][j+1][i];          \
  p_qsm1 = p_qs;                   \
  p_qs   = p_qsp1;                 \
  p_qsp1 = qs[k][j+1][i];          \
  p_rho_im1 = p_rho_i;             \
  p_rho_i   = p_rho_ip1;           \
  p_rho_ip1 = rho_i[k][j+1][i];    \
  p_squarem1 = p_square;           \
  p_square   = p_squarep1;         \
  p_squarep1 = square[k][j+1][i];

#define LOOP_PROLOGUE_HALF         \
  p_um2[0] = p_um1[0];             \
  p_um2[1] = p_um1[1];             \
  p_um2[2] = p_um1[2];             \
  p_um2[3] = p_um1[3];             \
  p_um2[4] = p_um1[4];             \
  p_um1[0] = p_u[0];               \
  p_um1[1] = p_u[1];               \
  p_um1[2] = p_u[2];               \
  p_um1[3] = p_u[3];               \
  p_um1[4] = p_u[4];               \
  p_u[0]   = p_up1[0];             \
  p_u[1]   = p_up1[1];             \
  p_u[2]   = p_up1[2];             \
  p_u[3]   = p_up1[3];             \
  p_u[4]   = p_up1[4];             \
  p_up1[0] = p_up2[0];             \
  p_up1[1] = p_up2[1];             \
  p_up1[2] = p_up2[2];             \
  p_up1[3] = p_up2[3];             \
  p_up1[4] = p_up2[4];             \
                                   \
  p_usm1 = p_us;                   \
  p_us   = p_usp1;                 \
  p_usp1 = us[k][j+1][i];          \
  p_vsm1 = p_vs;                   \
  p_vs   = p_vsp1;                 \
  p_vsp1 = vs[k][j+1][i];          \
  p_wsm1 = p_ws;                   \
  p_ws   = p_wsp1;                 \
  p_wsp1 = ws[k][j+1][i];          \
  p_qsm1 = p_qs;                   \
  p_qs   = p_qsp1;                 \
  p_qsp1 = qs[k][j+1][i];          \
  p_rho_im1 = p_rho_i;             \
  p_rho_i   = p_rho_ip1;           \
  p_rho_ip1 = rho_i[k][j+1][i];    \
  p_squarem1 = p_square;           \
  p_square   = p_squarep1;         \
  p_squarep1 = square[k][j+1][i];

#define LOOP_BODY                                                      \
   p_rhs[0] = rhs[k][j][i][0] + dy1ty1 *                               \
     (p_up1[0] - 2.0*p_u[0] + p_um1[0]) -                              \
     ty2 * (p_up1[2] - p_um1[2]);                                      \
                                                                       \
   p_rhs[1] = rhs[k][j][i][1] + dy2ty1 *                               \
     (p_up1[1] - 2.0*p_u[1] + p_um1[1]) +                              \
     yycon2 * (p_usp1 - 2.0*p_us + p_usm1) -                           \
     ty2 * (p_up1[1]*p_vsp1 - p_um1[1]*p_vsm1);                        \
                                                                       \
   p_rhs[2] = rhs[k][j][i][2] + dy3ty1 *                               \
     (p_up1[2] - 2.0*p_u[2] + p_um1[2]) +                              \
     yycon2*con43 * (p_vsp1 - 2.0*p_vs + p_vsm1) -                     \
     ty2 * (p_up1[2]*p_vsp1 - p_um1[2]*p_vsm1 +                        \
           (p_up1[4] - p_squarep1 -                                    \
            p_um1[4] + p_squarem1) * c2);                              \
                                                                       \
   p_rhs[3] = rhs[k][j][i][3] + dy4ty1 *                               \
     (p_up1[3] - 2.0*p_u[3] + p_um1[3]) +                              \
     yycon2 * (p_wsp1 - 2.0*p_ws + p_wsm1) -                           \
     ty2 * (p_up1[3]*p_vsp1 - p_um1[3]*p_vsm1);                        \
                                                                       \
   p_rhs[4] = rhs[k][j][i][4] + dy5ty1 *                               \
     (p_up1[4] - 2.0*p_u[4] + p_um1[4]) +                              \
     yycon3 * (p_qsp1 - 2.0*p_qs + p_qsm1) +                           \
     yycon4 * (p_vsp1*p_vsp1       - 2.0*p_vs*p_vs + p_vsm1*p_vsm1) +  \
     yycon5 * (p_up1[4]*p_rho_ip1 -                                    \
             2.0*p_u[4]*p_rho_i +                                      \
               p_um1[4]*p_rho_im1) -                                   \
     ty2 * ((c1*p_up1[4] - c2*p_squarep1) * p_vsp1 -                   \
            (c1*p_um1[4] - c2*p_squarem1) * p_vsm1);

  //---------------------------------------------------------------------
  // add fourth order eta-direction dissipation
  //---------------------------------------------------------------------
  j = 1;
  LOOP_PROLOGUE_FULL
  LOOP_BODY
  for (m = 0; m < 5; m++) {
    rhs[k][j][i][m] = p_rhs[m]- dssp *
      ( 5.0*p_u[m] - 4.0*p_up1[m] + p_up2[m]);
  }

  j = 2;
  LOOP_PROLOGUE_FULL
  LOOP_BODY
  for (m = 0; m < 5; m++) {
    rhs[k][j][i][m] = p_rhs[m] - dssp *
      (-4.0*p_um1[m] + 6.0*p_u[m] -
        4.0*p_up1[m] + p_up2[m]);
  }

  for (j = 3; j <= gp1-4; j++) {
    LOOP_PROLOGUE_FULL
    LOOP_BODY
    for (m = 0; m < 5; m++) {
      rhs[k][j][i][m] = p_rhs[m] - dssp *
        ( p_um2[m] - 4.0*p_um1[m] +
        6.0*p_u[m] - 4.0*p_up1[m] +
          p_up2[m] );
    }
  }

  j = gp1-3;
  LOOP_PROLOGUE_FULL
  LOOP_BODY
  for (m = 0; m < 5; m++) {
    rhs[k][j][i][m] = p_rhs[m] - dssp *
      ( p_um2[m] - 4.0*p_um1[m] +
      6.0*p_u[m] - 4.0*p_up1[m] );
  }

  j = gp1-2;
  LOOP_PROLOGUE_HALF
  LOOP_BODY
  for (m = 0; m < 5; m++) {
    rhs[k][j][i][m] = p_rhs[m] - dssp *
      ( p_um2[m] - 4.0*p_um1[m] + 5.0*p_u[m] );
  }

#undef LOOP_PROLOGUE_FULL
#undef LOOP_PROLOGUE_HALF
#undef LOOP_BODY
}
