             qr[mq] = qold[I_q];
             qd[mq] = qold[I_q - ys];
 
-            //bmdq[mq] = bmdq_trans[I_q];
             bpdq[mq] = bpdq_trans[I_q];
         }            
 

            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++)

            for(int mq = 0; mq < meqn; mq++)
            {
                int I_q = I + mq*zs;
                qd[mq] = qold[I_q - ys];   /* Down  */  
            }

            for(int m = 0; m < maux; m++)
            {
                int I_aux = I + m*zs;
                auxd[m] = aux[I_aux - ys];
            }               

