             qd[mq] = qold[I_q - ys];
 
             bmdq[mq] = bmdq_trans[I_q];
-            //bpdq[mq] = bpdq_trans[I_q];
         }            
 
         for(int imp = 0; imp < 2; imp++)


        {
            /* ------------------------ Normal solve in Y direction ------------------- */
            double *const qr     = start;                 /* meqn        */
            double *const qd     = qr     + meqn;         /* meqn        */
            double *const auxr   = qd     + meqn;         /* maux        */
            double *const auxd   = auxr   + maux;         /* maux        */
            double *const s      = auxd   + maux;         /* mwaves      */
            double *const wave   = s      + mwaves;       /* meqn*mwaves */
            double *const bmdq   = wave   + meqn*mwaves;  /* meqn        */

            ql[mq] = qold[I_q-1];
            qr[mq] = qold[I_q];

            apdq[mq] = apdq_trans[I_q];
        }            

