     int *pfacn, *pfccls;
     double *pidsoln, *pjdsoln, *pdsoln;
     double *pfcnml;
+    double (*pten)[NDIM];
 #if NDIM == 3
     int *pfcnds;
     double *pndcrd, *pfccnd;

            vec[ieq][0] = mat[0][0]*pdsoln[0] + mat[0][1]*pdsoln[1];
            vec[ieq][1] = mat[1][0]*pdsoln[0] + mat[1][1]*pdsoln[1];
#endif
            pdsoln += NDIM;
        };
        // set wall condition in the rotated coordinate;
        vec[0][0] = -vec[0][0];

        mat[1][2] = pndcrd[2] - pfccnd[2];
        len = sqrt(mat[1][0]*mat[1][0] + mat[1][1]*mat[1][1]
                 + mat[1][2]*mat[1][2]);
        mat[1][0] = mvt[0][1] = mat[1][0]/len;
        mat[1][1] = mvt[1][1] = mat[1][1]/len;
        mat[1][2] = mvt[2][1] = mat[1][2]/len;
        mat[2][0] = mvt[0][2] = mat[0][1]*mat[1][2] - mat[0][2]*mat[1][1];
        mat[2][1] = mvt[1][2] = mat[0][2]*mat[1][0] - mat[0][0]*mat[1][2];
        mat[2][2] = mvt[2][2] = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
#else
        mat[1][0] = mvt[0][1] =  pfcnml[1];
        mat[1][1] = mvt[1][1] = -pfcnml[0];
#endif
        // rotate momentum vector.
#if NDIM == 3
