diff --git a/src/gpu/core/GPU/CudaMemoryManager.cpp b/src/gpu/core/GPU/CudaMemoryManager.cpp
index 312e379c476b54e868e97e2d1b5fb88311bb4158..a6adb34c165019372512873d153f45681c01ea11 100644
--- a/src/gpu/core/GPU/CudaMemoryManager.cpp
+++ b/src/gpu/core/GPU/CudaMemoryManager.cpp
@@ -1863,82 +1863,6 @@ void CudaMemoryManager::cudaFreePressX1(int lev)
     checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->QpressX1.RhoBC  ));
     checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->QpressX1.deltaVz));
 }
-//Propeller Velocity
-void CudaMemoryManager::cudaAllocVeloPropeller(int lev)
-{
-    unsigned int mem_size_Propeller_k = sizeof(int)*parameter->getParH(lev)->propellerBC.numberOfBCnodes;
-    unsigned int mem_size_Propeller_q = sizeof(real)*parameter->getParH(lev)->propellerBC.numberOfBCnodes;
-
-    //Host
-    //checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->propellerBC.q27[0]),  parameter->getD3Qxx()*mem_size_Propeller_q ));
-    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->propellerBC.k),                  mem_size_Propeller_k ));
-    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->propellerBC.Vx),                 mem_size_Propeller_q ));
-    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->propellerBC.Vy),                 mem_size_Propeller_q ));
-    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->propellerBC.Vz),                 mem_size_Propeller_q ));
-    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->propellerBC.RhoBC),              mem_size_Propeller_q ));
-
-    //Device
-    //checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->propellerBC.q27[0]),      parameter->getD3Qxx()*mem_size_Propeller_q ));
-    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->propellerBC.k),                      mem_size_Propeller_k ));
-    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->propellerBC.Vx),                     mem_size_Propeller_q ));
-    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->propellerBC.Vy),                     mem_size_Propeller_q ));
-    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->propellerBC.Vz),                     mem_size_Propeller_q ));
-    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->propellerBC.RhoBC),                  mem_size_Propeller_q ));
-
-    //////////////////////////////////////////////////////////////////////////
-    double tmp = (double)mem_size_Propeller_k + 4. * (double)mem_size_Propeller_q;
-    setMemsizeGPU(tmp, false);
-}
-void CudaMemoryManager::cudaCopyVeloPropeller(int lev)
-{
-    unsigned int mem_size_Propeller_k = sizeof(int)*parameter->getParH(lev)->propellerBC.numberOfBCnodes;
-    unsigned int mem_size_Propeller_q = sizeof(real)*parameter->getParH(lev)->propellerBC.numberOfBCnodes;
-
-    //checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->propellerBC.q27[0],  parameter->getParH(lev)->propellerBC.q27[0], parameter->getD3Qxx()* mem_size_Propeller_q,  cudaMemcpyHostToDevice));
-    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->propellerBC.k,       parameter->getParH(lev)->propellerBC.k,                  mem_size_Propeller_k,  cudaMemcpyHostToDevice));
-    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->propellerBC.Vx,      parameter->getParH(lev)->propellerBC.Vx,                 mem_size_Propeller_q,  cudaMemcpyHostToDevice));
-    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->propellerBC.Vy,      parameter->getParH(lev)->propellerBC.Vy,                 mem_size_Propeller_q,  cudaMemcpyHostToDevice));
-    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->propellerBC.Vz,      parameter->getParH(lev)->propellerBC.Vz,                 mem_size_Propeller_q,  cudaMemcpyHostToDevice));
-    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->propellerBC.RhoBC,   parameter->getParH(lev)->propellerBC.RhoBC,              mem_size_Propeller_q,  cudaMemcpyHostToDevice));
-}
-void CudaMemoryManager::cudaFreeVeloPropeller(int lev)
-{
-    //checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->propellerBC.q27[0] ));
-    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->propellerBC.k      ));
-    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->propellerBC.Vx     ));
-    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->propellerBC.Vy     ));
-    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->propellerBC.Vz     ));
-    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->propellerBC.RhoBC  ));
-}
-//Measure Points
-//void CudaMemoryManager::cudaAllocMeasurePoints(int lev, int i)
-//{
-//    //Host
-//    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->MP[i].Vx),                 parameter->getParH(lev)->memSizerealMP ));
-//    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->MP[i].Vy),                 parameter->getParH(lev)->memSizerealMP ));
-//    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->MP[i].Vz),                 parameter->getParH(lev)->memSizerealMP ));
-//    checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->MP[i].Rho),                parameter->getParH(lev)->memSizerealMP ));
-//
-//    //Device
-//    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->MP[i].Vx),                     parameter->getParD(lev)->memSizerealMP ));
-//    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->MP[i].Vy),                     parameter->getParD(lev)->memSizerealMP ));
-//    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->MP[i].Vz),                     parameter->getParD(lev)->memSizerealMP ));
-//    checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->MP[i].Rho),                    parameter->getParD(lev)->memSizerealMP ));
-//}
-//void CudaMemoryManager::cudaCopyMeasurePoints(int lev, int i)
-//{
-//    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->MP[i].Vx,      parameter->getParH(lev)->MP[i].Vx,           parameter->getParH(lev)->memSizerealMP,  cudaMemcpyHostToDevice));
-//    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->MP[i].Vy,      parameter->getParH(lev)->MP[i].Vy,           parameter->getParH(lev)->memSizerealMP,  cudaMemcpyHostToDevice));
-//    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->MP[i].Vz,      parameter->getParH(lev)->MP[i].Vz,           parameter->getParH(lev)->memSizerealMP,  cudaMemcpyHostToDevice));
-//    checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->MP[i].Rho,     parameter->getParH(lev)->MP[i].Rho,          parameter->getParH(lev)->memSizerealMP,  cudaMemcpyHostToDevice));
-//}
-//void CudaMemoryManager::cudaFreeMeasurePoints(int lev, int i)
-//{
-//    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->MP[i].Vx     ));
-//    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->MP[i].Vy     ));
-//    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->MP[i].Vz     ));
-//    checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->MP[i].Rho    ));
-//}
 void CudaMemoryManager::cudaAllocMeasurePointsIndex(int lev)
 {
     //Host
diff --git a/src/gpu/core/GPU/CudaMemoryManager.h b/src/gpu/core/GPU/CudaMemoryManager.h
index dcdfecd6aad6ace700bc5e83e1d6e3f14e1f1bcc..89fed68436e4d40be2bf1abc9959825301a41488 100644
--- a/src/gpu/core/GPU/CudaMemoryManager.h
+++ b/src/gpu/core/GPU/CudaMemoryManager.h
@@ -38,9 +38,6 @@ public:
     void setMemsizeGPU(double admem, bool reset);
     double getMemsizeGPU();
 
-    //void cudaAllocFull(int lev); //DEPRECATED: related to full matrix
-    //void cudaFreeFull(int lev);  //DEPRECATED: related to full matrix
-
     void cudaCopyPrint(int lev);
     void cudaCopyMedianPrint(int lev);
 
@@ -251,14 +248,6 @@ public:
     void cudaCopyPressX1(int lev);
     void cudaFreePressX1(int lev);
 
-    void cudaAllocVeloPropeller(int lev);
-    void cudaCopyVeloPropeller(int lev);
-    void cudaFreeVeloPropeller(int lev);
-
-    void cudaAllocMeasurePoints(int lev, int i);
-    void cudaCopyMeasurePoints(int lev, int i);
-    void cudaFreeMeasurePoints(int lev, int i);
-
     void cudaAllocMeasurePointsIndex(int lev);
     void cudaCopyMeasurePointsIndex(int lev);
     void cudaCopyMeasurePointsToHost(int lev);
diff --git a/src/gpu/core/GPU/GPU_Interface.h b/src/gpu/core/GPU/GPU_Interface.h
index 72e839bcf741b9966a4e38a7b78c1085b20132e2..2f25f5f15cb5f9334c1cef070ba9c32a7c11c68b 100644
--- a/src/gpu/core/GPU/GPU_Interface.h
+++ b/src/gpu/core/GPU/GPU_Interface.h
@@ -933,21 +933,6 @@ void QADPressIncompDev27(  unsigned int numberOfThreads,
                                       unsigned long long numberOfLBnodes, 
                                       bool isEvenTimestep);
 
-void PropVelo(   unsigned int numberOfThreads,
-                            unsigned int* neighborX,
-                            unsigned int* neighborY,
-                            unsigned int* neighborZ,
-                            real* rho,
-                            real* ux,
-                            real* uy,
-                            real* uz,
-                            int* k_Q, 
-                            unsigned int size_Prop,
-                            unsigned long long numberOfLBnodes,
-                            unsigned int* bcMatD,
-                            real* DD,
-                            bool EvenOrOdd);
-
 void ScaleCF27( real* DC, 
                            real* DF, 
                            unsigned int* neighborCX,
diff --git a/src/gpu/core/GPU/GPU_Kernels.cuh b/src/gpu/core/GPU/GPU_Kernels.cuh
index ab8ab0fa3564de523e8fc8e6985f6d98c09fb404..6e007fb6ced0b4301b0f212d4b0199ae4988b728 100644
--- a/src/gpu/core/GPU/GPU_Kernels.cuh
+++ b/src/gpu/core/GPU/GPU_Kernels.cuh
@@ -1264,21 +1264,6 @@ __global__ void QADPressIncomp27(   real* DD,
                                                unsigned long long numberOfLBnodes,
                                                bool isEvenTimestep);
 
-//Propeller BC
-__global__ void PropellerBC(unsigned int* neighborX,
-                                       unsigned int* neighborY,
-                                       unsigned int* neighborZ,
-                                       real* rho,
-                                       real* ux,
-                                       real* uy,
-                                       real* uz,
-                                       int* k_Q,
-                                       unsigned int size_Prop,
-                                       unsigned long long numberOfLBnodes,
-                                       unsigned int* bcMatD,
-                                       real* DD,
-                                       bool EvenOrOdd);
-
 
 
 //coarse to fine
diff --git a/src/gpu/core/GPU/LBMKernel.cu b/src/gpu/core/GPU/LBMKernel.cu
index 403ed748f3d525bb5eaefe35b44ef26391230fcd..a870d43bb36811441afa654a01f4eeb7671ba28e 100644
--- a/src/gpu/core/GPU/LBMKernel.cu
+++ b/src/gpu/core/GPU/LBMKernel.cu
@@ -2688,41 +2688,6 @@ void QPrecursorDevDistributions( LBMSimulationParameter* parameterDevice,
 
 }
 //////////////////////////////////////////////////////////////////////////
-extern "C" void PropVelo(
-    unsigned int numberOfThreads,
-    unsigned int* neighborX,
-    unsigned int* neighborY,
-    unsigned int* neighborZ,
-    real* rho,
-    real* ux,
-    real* uy,
-    real* uz,
-    int* k_Q,
-    unsigned int size_Prop,
-    unsigned long long numberOfLBnodes,
-    unsigned int* bcMatD,
-    real* DD,
-    bool EvenOrOdd)
-{
-    vf::cuda::CudaGrid grid = vf::cuda::CudaGrid(numberOfThreads, size_Prop);
-
-    PropellerBC<<< grid.grid, grid.threads >>>(
-        neighborX,
-        neighborY,
-        neighborZ,
-        rho,
-        ux,
-        uy,
-        uz,
-        k_Q,
-        size_Prop,
-        numberOfLBnodes,
-        bcMatD,
-        DD,
-        EvenOrOdd);
-    getLastCudaError("PropellerBC execution failed");
-}
-//////////////////////////////////////////////////////////////////////////
 void ScaleCF27(
     real* DC,
     real* DF,
diff --git a/src/gpu/core/GPU/VelocityBCs27.cu b/src/gpu/core/GPU/VelocityBCs27.cu
index d1798e84490ffed3e803717e22cc8dbbaccf2222..437c19d81e513c70fbdddbb8e20819933d1c42d0 100644
--- a/src/gpu/core/GPU/VelocityBCs27.cu
+++ b/src/gpu/core/GPU/VelocityBCs27.cu
@@ -5444,401 +5444,3 @@ __global__ void QVelDevice27(
    }
 }
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-////////////////////////////////////////////////////////////////////////////////
-__global__ void PropellerBC(
-    unsigned int* neighborX,
-    unsigned int* neighborY,
-    unsigned int* neighborZ,
-    real* rho,
-    real* ux,
-    real* uy,
-    real* uz,
-    int* k_Q, 
-    unsigned int size_Prop,
-    unsigned long long numberOfLBnodes,
-    unsigned int* bcMatD,
-    real* DD,
-    bool EvenOrOdd)
-{
-   ////////////////////////////////////////////////////////////////////////////////
-   const unsigned  x = threadIdx.x;  // Globaler x-Index 
-   const unsigned  y = blockIdx.x;   // Globaler y-Index 
-   const unsigned  z = blockIdx.y;   // Globaler z-Index 
-
-   const unsigned nx = blockDim.x;
-   const unsigned ny = gridDim.x;
-
-   const unsigned k = nx*(ny*z + y) + x;
-   //////////////////////////////////////////////////////////////////////////
-
-   if(k<size_Prop)
-   {
-    ////////////////////////////////////////////////////////////////////////////////
-        Distributions27 D;
-        if (EvenOrOdd==true)
-        {
-            D.f[dP00] = &DD[dP00 * numberOfLBnodes];
-            D.f[dM00] = &DD[dM00 * numberOfLBnodes];
-            D.f[d0P0] = &DD[d0P0 * numberOfLBnodes];
-            D.f[d0M0] = &DD[d0M0 * numberOfLBnodes];
-            D.f[d00P] = &DD[d00P * numberOfLBnodes];
-            D.f[d00M] = &DD[d00M * numberOfLBnodes];
-            D.f[dPP0] = &DD[dPP0 * numberOfLBnodes];
-            D.f[dMM0] = &DD[dMM0 * numberOfLBnodes];
-            D.f[dPM0] = &DD[dPM0 * numberOfLBnodes];
-            D.f[dMP0] = &DD[dMP0 * numberOfLBnodes];
-            D.f[dP0P] = &DD[dP0P * numberOfLBnodes];
-            D.f[dM0M] = &DD[dM0M * numberOfLBnodes];
-            D.f[dP0M] = &DD[dP0M * numberOfLBnodes];
-            D.f[dM0P] = &DD[dM0P * numberOfLBnodes];
-            D.f[d0PP] = &DD[d0PP * numberOfLBnodes];
-            D.f[d0MM] = &DD[d0MM * numberOfLBnodes];
-            D.f[d0PM] = &DD[d0PM * numberOfLBnodes];
-            D.f[d0MP] = &DD[d0MP * numberOfLBnodes];
-            D.f[d000] = &DD[d000 * numberOfLBnodes];
-            D.f[dPPP] = &DD[dPPP * numberOfLBnodes];
-            D.f[dMMP] = &DD[dMMP * numberOfLBnodes];
-            D.f[dPMP] = &DD[dPMP * numberOfLBnodes];
-            D.f[dMPP] = &DD[dMPP * numberOfLBnodes];
-            D.f[dPPM] = &DD[dPPM * numberOfLBnodes];
-            D.f[dMMM] = &DD[dMMM * numberOfLBnodes];
-            D.f[dPMM] = &DD[dPMM * numberOfLBnodes];
-            D.f[dMPM] = &DD[dMPM * numberOfLBnodes];
-        }
-        else
-        {
-            D.f[dM00] = &DD[dP00 * numberOfLBnodes];
-            D.f[dP00] = &DD[dM00 * numberOfLBnodes];
-            D.f[d0M0] = &DD[d0P0 * numberOfLBnodes];
-            D.f[d0P0] = &DD[d0M0 * numberOfLBnodes];
-            D.f[d00M] = &DD[d00P * numberOfLBnodes];
-            D.f[d00P] = &DD[d00M * numberOfLBnodes];
-            D.f[dMM0] = &DD[dPP0 * numberOfLBnodes];
-            D.f[dPP0] = &DD[dMM0 * numberOfLBnodes];
-            D.f[dMP0] = &DD[dPM0 * numberOfLBnodes];
-            D.f[dPM0] = &DD[dMP0 * numberOfLBnodes];
-            D.f[dM0M] = &DD[dP0P * numberOfLBnodes];
-            D.f[dP0P] = &DD[dM0M * numberOfLBnodes];
-            D.f[dM0P] = &DD[dP0M * numberOfLBnodes];
-            D.f[dP0M] = &DD[dM0P * numberOfLBnodes];
-            D.f[d0MM] = &DD[d0PP * numberOfLBnodes];
-            D.f[d0PP] = &DD[d0MM * numberOfLBnodes];
-            D.f[d0MP] = &DD[d0PM * numberOfLBnodes];
-            D.f[d0PM] = &DD[d0MP * numberOfLBnodes];
-            D.f[d000] = &DD[d000 * numberOfLBnodes];
-            D.f[dMMM] = &DD[dPPP * numberOfLBnodes];
-            D.f[dPPM] = &DD[dMMP * numberOfLBnodes];
-            D.f[dMPM] = &DD[dPMP * numberOfLBnodes];
-            D.f[dPMM] = &DD[dMPP * numberOfLBnodes];
-            D.f[dMMP] = &DD[dPPM * numberOfLBnodes];
-            D.f[dPPP] = &DD[dMMM * numberOfLBnodes];
-            D.f[dMPP] = &DD[dPMM * numberOfLBnodes];
-            D.f[dPMP] = &DD[dMPM * numberOfLBnodes];
-        }
-        //////////////////////////////////////////////////////////////////////////
-        unsigned int KQK = k_Q[k];
-        unsigned int BC  = bcMatD[KQK];
-        if( (BC != GEO_SOLID) && (BC != GEO_VOID))
-        {        
-        //////////////////////////////////////////////////////////////////////////
-        real  vx1 = ux[k];
-        real  vx2 = uy[k];
-        real  vx3 = uz[k];
-        //real  vx1 = -c1o100;
-        //real  vx2 = zero;
-        //real  vx3 = zero;
-        //////////////////////////////////////////////////////////////////////////
-        //index
-        //////////////////////////////////////////////////////////////////////////
-        unsigned int kzero= KQK;
-        unsigned int ke   = KQK;
-        unsigned int kw   = neighborX[KQK];
-        unsigned int kn   = KQK;
-        unsigned int ks   = neighborY[KQK];
-        unsigned int kt   = KQK;
-        unsigned int kb   = neighborZ[KQK];
-        unsigned int ksw  = neighborY[kw];
-        unsigned int kne  = KQK;
-        unsigned int kse  = ks;
-        unsigned int knw  = kw;
-        unsigned int kbw  = neighborZ[kw];
-        unsigned int kte  = KQK;
-        unsigned int kbe  = kb;
-        unsigned int ktw  = kw;
-        unsigned int kbs  = neighborZ[ks];
-        unsigned int ktn  = KQK;
-        unsigned int kbn  = kb;
-        unsigned int kts  = ks;
-        unsigned int ktse = ks;
-        unsigned int kbnw = kbw;
-        unsigned int ktnw = kw;
-        unsigned int kbse = kbs;
-        unsigned int ktsw = ksw;
-        unsigned int kbne = kb;
-        unsigned int ktne = KQK;
-        unsigned int kbsw = neighborZ[ksw];
-        //////////////////////////////////////////////////////////////////////////
-        real f_E,  f_W,  f_N,  f_S,  f_T,  f_B,   f_NE,  f_SW,  f_SE,  f_NW,  f_TE,  f_BW,  f_BE,
-        f_TW, f_TN, f_BS, f_BN, f_TS, f_TNE, f_TSW, f_TSE, f_TNW, f_BNE, f_BSW, f_BSE, f_BNW, f_ZERO;
-
-        f_ZERO= (D.f[d000])[kzero];
-        f_E   = (D.f[dP00])[ke   ];
-        f_W   = (D.f[dM00])[kw   ];
-        f_N   = (D.f[d0P0])[kn   ];
-        f_S   = (D.f[d0M0])[ks   ];
-        f_T   = (D.f[d00P])[kt   ];
-        f_B   = (D.f[d00M])[kb   ];
-        f_NE  = (D.f[dPP0])[kne  ];
-        f_SW  = (D.f[dMM0])[ksw  ];
-        f_SE  = (D.f[dPM0])[kse  ];
-        f_NW  = (D.f[dMP0])[knw  ];
-        f_TE  = (D.f[dP0P])[kte  ];
-        f_BW  = (D.f[dM0M])[kbw  ];
-        f_BE  = (D.f[dP0M])[kbe  ];
-        f_TW  = (D.f[dM0P])[ktw  ];
-        f_TN  = (D.f[d0PP])[ktn  ];
-        f_BS  = (D.f[d0MM])[kbs  ];
-        f_BN  = (D.f[d0PM])[kbn  ];
-        f_TS  = (D.f[d0MP])[kts  ];
-        f_TNE = (D.f[dPPP])[ktne ];
-        f_BSW = (D.f[dMMM])[kbsw ];
-        f_BNE = (D.f[dPPM])[kbne ];
-        f_TSW = (D.f[dMMP])[ktsw ];
-        f_TSE = (D.f[dPMP])[ktse ];
-        f_BNW = (D.f[dMPM])[kbnw ];
-        f_BSE = (D.f[dPMM])[kbse ];
-        f_TNW = (D.f[dMPP])[ktnw ];
-        //f_W    = (D.f[dP00])[ke   ];
-        //f_E    = (D.f[dM00])[kw   ];
-        //f_S    = (D.f[d0P0])[kn   ];
-        //f_N    = (D.f[d0M0])[ks   ];
-        //f_B    = (D.f[d00P])[kt   ];
-        //f_T    = (D.f[d00M])[kb   ];
-        //f_SW   = (D.f[dPP0])[kne  ];
-        //f_NE   = (D.f[dMM0])[ksw  ];
-        //f_NW   = (D.f[dPM0])[kse  ];
-        //f_SE   = (D.f[dMP0])[knw  ];
-        //f_BW   = (D.f[dP0P])[kte  ];
-        //f_TE   = (D.f[dM0M])[kbw  ];
-        //f_TW   = (D.f[dP0M])[kbe  ];
-        //f_BE   = (D.f[dM0P])[ktw  ];
-        //f_BS   = (D.f[d0PP])[ktn  ];
-        //f_TN   = (D.f[d0MM])[kbs  ];
-        //f_TS   = (D.f[d0PM])[kbn  ];
-        //f_BN   = (D.f[d0MP])[kts  ];
-        //f_BSW  = (D.f[dPPP])[ktne ];
-        //f_TNE  = (D.f[dMMM])[kbsw ];
-        //f_TSW  = (D.f[dPPM])[kbne ];
-        //f_BNE  = (D.f[dMMP])[ktsw ];
-        //f_BNW  = (D.f[dPMP])[ktse ];
-        //f_TSE  = (D.f[dMPM])[kbnw ];
-        //f_TNW  = (D.f[dPMM])[kbse ];
-        //f_BSE  = (D.f[dMPP])[ktnw ];
-        //////////////////////////////////////////////////////////////////////////////////
-        real vxo1, vxo2, vxo3, drho;
-        drho   =  /*zero;*/f_TSE + f_TNW + f_TNE + f_TSW + f_BSE + f_BNW + f_BNE + f_BSW +
-                  f_BN + f_TS + f_TN + f_BS + f_BE + f_TW + f_TE + f_BW + f_SE + f_NW + f_NE + f_SW + 
-                  f_T + f_B + f_N + f_S + f_E + f_W + f_ZERO; 
-
-        vxo1   =   (((f_TSE - f_BNW) - (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
-                    ((f_BE - f_TW)   + (f_TE - f_BW))   + ((f_SE - f_NW)   + (f_NE - f_SW)) +
-                    (f_E - f_W) )/ (c1o1 + drho); 
-        
-
-        vxo2   =   ((-(f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
-                    ((f_BN - f_TS)   + (f_TN - f_BS))    + (-(f_SE - f_NW)  + (f_NE - f_SW)) +
-                    (f_N - f_S) )/ (c1o1 + drho); 
-
-        vxo3   =   (((f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) + (f_TSW - f_BNE)) +
-                     (-(f_BN - f_TS)  + (f_TN - f_BS))   + ((f_TE - f_BW)   - (f_BE - f_TW)) +
-                    (f_T - f_B) )/ (c1o1 + drho); 
-
-        real cusq=c3o2*(vxo1*vxo1+vxo2*vxo2+vxo3*vxo3);
-        //vx1 = vx1 * two - vxo1;
-        //vx2 = vx2 * two - vxo2;
-        //vx3 = vx3 * two - vxo3;
-        real cusq2=c3o2*(vx1*vx1+vx2*vx2+vx3*vx3);
-
-         //f_ZERO = ((one+drho) * (   c8over27 *(one+(-cusq2)))) - c8over27;
-         //f_E    = ((one+drho) * (   c2over27 *(one+three*( vx1        )+c9over2*( vx1        )*( vx1        )-cusq2))) - c2over27 ;
-         //f_W    = ((one+drho) * (   c2over27 *(one+three*(-vx1        )+c9over2*(-vx1        )*(-vx1        )-cusq2))) - c2over27 ;
-         //f_N    = ((one+drho) * (   c2over27 *(one+three*(    vx2     )+c9over2*(     vx2    )*(     vx2    )-cusq2))) - c2over27 ;
-         //f_S    = ((one+drho) * (   c2over27 *(one+three*(   -vx2     )+c9over2*(    -vx2    )*(    -vx2    )-cusq2))) - c2over27 ;
-         //f_T    = ((one+drho) * (   c2over27 *(one+three*(         vx3)+c9over2*(         vx3)*(         vx3)-cusq2))) - c2over27 ;
-         //f_B    = ((one+drho) * (   c2over27 *(one+three*(        -vx3)+c9over2*(        -vx3)*(        -vx3)-cusq2))) - c2over27 ;
-         //f_NE   = ((one+drho) * (   c1over54 *(one+three*( vx1+vx2    )+c9over2*( vx1+vx2    )*( vx1+vx2    )-cusq2))) - c1over54 ;
-         //f_SW   = ((one+drho) * (   c1over54 *(one+three*(-vx1-vx2    )+c9over2*(-vx1-vx2    )*(-vx1-vx2    )-cusq2))) - c1over54 ;
-         //f_SE   = ((one+drho) * (   c1over54 *(one+three*( vx1-vx2    )+c9over2*( vx1-vx2    )*( vx1-vx2    )-cusq2))) - c1over54 ;
-         //f_NW   = ((one+drho) * (   c1over54 *(one+three*(-vx1+vx2    )+c9over2*(-vx1+vx2    )*(-vx1+vx2    )-cusq2))) - c1over54 ;
-         //f_TE   = ((one+drho) * (   c1over54 *(one+three*( vx1    +vx3)+c9over2*( vx1    +vx3)*( vx1    +vx3)-cusq2))) - c1over54 ;
-         //f_BW   = ((one+drho) * (   c1over54 *(one+three*(-vx1    -vx3)+c9over2*(-vx1    -vx3)*(-vx1    -vx3)-cusq2))) - c1over54 ;
-         //f_BE   = ((one+drho) * (   c1over54 *(one+three*( vx1    -vx3)+c9over2*( vx1    -vx3)*( vx1    -vx3)-cusq2))) - c1over54 ;
-         //f_TW   = ((one+drho) * (   c1over54 *(one+three*(-vx1    +vx3)+c9over2*(-vx1    +vx3)*(-vx1    +vx3)-cusq2))) - c1over54 ;
-         //f_TN   = ((one+drho) * (   c1over54 *(one+three*(     vx2+vx3)+c9over2*(     vx2+vx3)*(     vx2+vx3)-cusq2))) - c1over54 ;
-         //f_BS   = ((one+drho) * (   c1over54 *(one+three*(    -vx2-vx3)+c9over2*(    -vx2-vx3)*(    -vx2-vx3)-cusq2))) - c1over54 ;
-         //f_BN   = ((one+drho) * (   c1over54 *(one+three*(     vx2-vx3)+c9over2*(     vx2-vx3)*(     vx2-vx3)-cusq2))) - c1over54 ;
-         //f_TS   = ((one+drho) * (   c1over54 *(one+three*(    -vx2+vx3)+c9over2*(    -vx2+vx3)*(    -vx2+vx3)-cusq2))) - c1over54 ;
-         //f_TNE  = ((one+drho) * (   c1over216*(one+three*( vx1+vx2+vx3)+c9over2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cusq2))) - c1over216;
-         //f_BSW  = ((one+drho) * (   c1over216*(one+three*(-vx1-vx2-vx3)+c9over2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cusq2))) - c1over216;
-         //f_BNE  = ((one+drho) * (   c1over216*(one+three*( vx1+vx2-vx3)+c9over2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cusq2))) - c1over216;
-         //f_TSW  = ((one+drho) * (   c1over216*(one+three*(-vx1-vx2+vx3)+c9over2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cusq2))) - c1over216;
-         //f_TSE  = ((one+drho) * (   c1over216*(one+three*( vx1-vx2+vx3)+c9over2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cusq2))) - c1over216;
-         //f_BNW  = ((one+drho) * (   c1over216*(one+three*(-vx1+vx2-vx3)+c9over2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cusq2))) - c1over216;
-         //f_BSE  = ((one+drho) * (   c1over216*(one+three*( vx1-vx2-vx3)+c9over2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cusq2))) - c1over216;
-         //f_TNW  = ((one+drho) * (   c1over216*(one+three*(-vx1+vx2+vx3)+c9over2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cusq2))) - c1over216;
-         f_ZERO = f_ZERO + ((c1o1+drho) * (-  c8o27* (-cusq)                                                                   +   c8o27* (-cusq2)));
-         f_E    = f_E    + ((c1o1+drho) * (-  c2o27* (c3o1*( vxo1          )+c9o2*( vxo1          )*( vxo1          )-cusq) +   c2o27* (c3o1*( vx1        )+c9o2*( vx1        )*( vx1        )-cusq2)));
-         f_W    = f_W    + ((c1o1+drho) * (-  c2o27* (c3o1*(-vxo1          )+c9o2*(-vxo1          )*(-vxo1          )-cusq) +   c2o27* (c3o1*(-vx1        )+c9o2*(-vx1        )*(-vx1        )-cusq2)));
-         f_N    = f_N    + ((c1o1+drho) * (-  c2o27* (c3o1*(      vxo2     )+c9o2*(      vxo2     )*(      vxo2     )-cusq) +   c2o27* (c3o1*(    vx2     )+c9o2*(     vx2    )*(     vx2    )-cusq2)));
-         f_S    = f_S    + ((c1o1+drho) * (-  c2o27* (c3o1*(     -vxo2     )+c9o2*(     -vxo2     )*(     -vxo2     )-cusq) +   c2o27* (c3o1*(   -vx2     )+c9o2*(    -vx2    )*(    -vx2    )-cusq2)));
-         f_T    = f_T    + ((c1o1+drho) * (-  c2o27* (c3o1*(           vxo3)+c9o2*(           vxo3)*(           vxo3)-cusq) +   c2o27* (c3o1*(         vx3)+c9o2*(         vx3)*(         vx3)-cusq2)));
-         f_B    = f_B    + ((c1o1+drho) * (-  c2o27* (c3o1*(          -vxo3)+c9o2*(          -vxo3)*(          -vxo3)-cusq) +   c2o27* (c3o1*(        -vx3)+c9o2*(        -vx3)*(        -vx3)-cusq2)));
-         f_NE   = f_NE   + ((c1o1+drho) * (-  c1o54* (c3o1*( vxo1+vxo2     )+c9o2*( vxo1+vxo2     )*( vxo1+vxo2     )-cusq) +   c1o54* (c3o1*( vx1+vx2    )+c9o2*( vx1+vx2    )*( vx1+vx2    )-cusq2)));
-         f_SW   = f_SW   + ((c1o1+drho) * (-  c1o54* (c3o1*(-vxo1-vxo2     )+c9o2*(-vxo1-vxo2     )*(-vxo1-vxo2     )-cusq) +   c1o54* (c3o1*(-vx1-vx2    )+c9o2*(-vx1-vx2    )*(-vx1-vx2    )-cusq2)));
-         f_SE   = f_SE   + ((c1o1+drho) * (-  c1o54* (c3o1*( vxo1-vxo2     )+c9o2*( vxo1-vxo2     )*( vxo1-vxo2     )-cusq) +   c1o54* (c3o1*( vx1-vx2    )+c9o2*( vx1-vx2    )*( vx1-vx2    )-cusq2)));
-         f_NW   = f_NW   + ((c1o1+drho) * (-  c1o54* (c3o1*(-vxo1+vxo2     )+c9o2*(-vxo1+vxo2     )*(-vxo1+vxo2     )-cusq) +   c1o54* (c3o1*(-vx1+vx2    )+c9o2*(-vx1+vx2    )*(-vx1+vx2    )-cusq2)));
-         f_TE   = f_TE   + ((c1o1+drho) * (-  c1o54* (c3o1*( vxo1     +vxo3)+c9o2*( vxo1     +vxo3)*( vxo1     +vxo3)-cusq) +   c1o54* (c3o1*( vx1    +vx3)+c9o2*( vx1    +vx3)*( vx1    +vx3)-cusq2)));
-         f_BW   = f_BW   + ((c1o1+drho) * (-  c1o54* (c3o1*(-vxo1     -vxo3)+c9o2*(-vxo1     -vxo3)*(-vxo1     -vxo3)-cusq) +   c1o54* (c3o1*(-vx1    -vx3)+c9o2*(-vx1    -vx3)*(-vx1    -vx3)-cusq2)));
-         f_BE   = f_BE   + ((c1o1+drho) * (-  c1o54* (c3o1*( vxo1     -vxo3)+c9o2*( vxo1     -vxo3)*( vxo1     -vxo3)-cusq) +   c1o54* (c3o1*( vx1    -vx3)+c9o2*( vx1    -vx3)*( vx1    -vx3)-cusq2)));
-         f_TW   = f_TW   + ((c1o1+drho) * (-  c1o54* (c3o1*(-vxo1     +vxo3)+c9o2*(-vxo1     +vxo3)*(-vxo1     +vxo3)-cusq) +   c1o54* (c3o1*(-vx1    +vx3)+c9o2*(-vx1    +vx3)*(-vx1    +vx3)-cusq2)));
-         f_TN   = f_TN   + ((c1o1+drho) * (-  c1o54* (c3o1*(      vxo2+vxo3)+c9o2*(      vxo2+vxo3)*(      vxo2+vxo3)-cusq) +   c1o54* (c3o1*(     vx2+vx3)+c9o2*(     vx2+vx3)*(     vx2+vx3)-cusq2)));
-         f_BS   = f_BS   + ((c1o1+drho) * (-  c1o54* (c3o1*(     -vxo2-vxo3)+c9o2*(     -vxo2-vxo3)*(     -vxo2-vxo3)-cusq) +   c1o54* (c3o1*(    -vx2-vx3)+c9o2*(    -vx2-vx3)*(    -vx2-vx3)-cusq2)));
-         f_BN   = f_BN   + ((c1o1+drho) * (-  c1o54* (c3o1*(      vxo2-vxo3)+c9o2*(      vxo2-vxo3)*(      vxo2-vxo3)-cusq) +   c1o54* (c3o1*(     vx2-vx3)+c9o2*(     vx2-vx3)*(     vx2-vx3)-cusq2)));
-         f_TS   = f_TS   + ((c1o1+drho) * (-  c1o54* (c3o1*(     -vxo2+vxo3)+c9o2*(     -vxo2+vxo3)*(     -vxo2+vxo3)-cusq) +   c1o54* (c3o1*(    -vx2+vx3)+c9o2*(    -vx2+vx3)*(    -vx2+vx3)-cusq2)));
-         f_TNE  = f_TNE  + ((c1o1+drho) * (-  c1o216*(c3o1*( vxo1+vxo2+vxo3)+c9o2*( vxo1+vxo2+vxo3)*( vxo1+vxo2+vxo3)-cusq) +   c1o216*(c3o1*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cusq2)));
-         f_BSW  = f_BSW  + ((c1o1+drho) * (-  c1o216*(c3o1*(-vxo1-vxo2-vxo3)+c9o2*(-vxo1-vxo2-vxo3)*(-vxo1-vxo2-vxo3)-cusq) +   c1o216*(c3o1*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cusq2)));
-         f_BNE  = f_BNE  + ((c1o1+drho) * (-  c1o216*(c3o1*( vxo1+vxo2-vxo3)+c9o2*( vxo1+vxo2-vxo3)*( vxo1+vxo2-vxo3)-cusq) +   c1o216*(c3o1*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cusq2)));
-         f_TSW  = f_TSW  + ((c1o1+drho) * (-  c1o216*(c3o1*(-vxo1-vxo2+vxo3)+c9o2*(-vxo1-vxo2+vxo3)*(-vxo1-vxo2+vxo3)-cusq) +   c1o216*(c3o1*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cusq2)));
-         f_TSE  = f_TSE  + ((c1o1+drho) * (-  c1o216*(c3o1*( vxo1-vxo2+vxo3)+c9o2*( vxo1-vxo2+vxo3)*( vxo1-vxo2+vxo3)-cusq) +   c1o216*(c3o1*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cusq2)));
-         f_BNW  = f_BNW  + ((c1o1+drho) * (-  c1o216*(c3o1*(-vxo1+vxo2-vxo3)+c9o2*(-vxo1+vxo2-vxo3)*(-vxo1+vxo2-vxo3)-cusq) +   c1o216*(c3o1*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cusq2)));
-         f_BSE  = f_BSE  + ((c1o1+drho) * (-  c1o216*(c3o1*( vxo1-vxo2-vxo3)+c9o2*( vxo1-vxo2-vxo3)*( vxo1-vxo2-vxo3)-cusq) +   c1o216*(c3o1*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cusq2)));
-         f_TNW  = f_TNW  + ((c1o1+drho) * (-  c1o216*(c3o1*(-vxo1+vxo2+vxo3)+c9o2*(-vxo1+vxo2+vxo3)*(-vxo1+vxo2+vxo3)-cusq) +   c1o216*(c3o1*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cusq2)));
-
-        (D.f[d000])[kzero] =  f_ZERO;
-        (D.f[dP00])[ke   ] =  f_E   ;    // f_W   ;//        
-        (D.f[dM00])[kw   ] =  f_W   ;    // f_E   ;//        
-        (D.f[d0P0])[kn   ] =  f_N   ;    // f_S   ;//        
-        (D.f[d0M0])[ks   ] =  f_S   ;    // f_N   ;//        
-        (D.f[d00P])[kt   ] =  f_T   ;    // f_B   ;//        
-        (D.f[d00M])[kb   ] =  f_B   ;    // f_T   ;//        
-        (D.f[dPP0])[kne  ] =  f_NE  ;    // f_SW  ;//        
-        (D.f[dMM0])[ksw  ] =  f_SW  ;    // f_NE  ;//        
-        (D.f[dPM0])[kse  ] =  f_SE  ;    // f_NW  ;//        
-        (D.f[dMP0])[knw  ] =  f_NW  ;    // f_SE  ;//        
-        (D.f[dP0P])[kte  ] =  f_TE  ;    // f_BW  ;//        
-        (D.f[dM0M])[kbw  ] =  f_BW  ;    // f_TE  ;//        
-        (D.f[dP0M])[kbe  ] =  f_BE  ;    // f_TW  ;//        
-        (D.f[dM0P])[ktw  ] =  f_TW  ;    // f_BE  ;//        
-        (D.f[d0PP])[ktn  ] =  f_TN  ;    // f_BS  ;//        
-        (D.f[d0MM])[kbs  ] =  f_BS  ;    // f_TN  ;//        
-        (D.f[d0PM])[kbn  ] =  f_BN  ;    // f_TS  ;//        
-        (D.f[d0MP])[kts  ] =  f_TS  ;    // f_BN  ;//        
-        (D.f[dPPP])[ktne ] =  f_TNE ;    // f_BSW ;//        
-        (D.f[dMMM])[kbsw ] =  f_BSW ;    // f_BNE ;//        
-        (D.f[dPPM])[kbne ] =  f_BNE ;    // f_BNW ;//        
-        (D.f[dMMP])[ktsw ] =  f_TSW ;    // f_BSE ;//        
-        (D.f[dPMP])[ktse ] =  f_TSE ;    // f_TSW ;//        
-        (D.f[dMPM])[kbnw ] =  f_BNW ;    // f_TNE ;//        
-        (D.f[dPMM])[kbse ] =  f_BSE ;    // f_TNW ;//        
-        (D.f[dMPP])[ktnw ] =  f_TNW ;    // f_TSE ;//        
-
-        //////////////////////////////////////////////////////////////////////////
-        ////(D.f[d000])[kzero] =   c8over27* (drho-cu_sq);
-        //(D.f[dP00])[ke   ] =   three*c2over27* ( vx1        );        //six
-        //(D.f[dM00])[kw   ] =   three*c2over27* (-vx1        );        //six
-        //(D.f[d0P0])[kn   ] =   three*c2over27* (     vx2    );        //six
-        //(D.f[d0M0])[ks   ] =   three*c2over27* (    -vx2    );        //six
-        //(D.f[d00P])[kt   ] =   three*c2over27* (         vx3);        //six
-        //(D.f[d00M])[kb   ] =   three*c2over27* (        -vx3);        //six
-        //(D.f[dPP0])[kne  ] =   three*c1over54* ( vx1+vx2    );        //six
-        //(D.f[dMM0])[ksw  ] =   three*c1over54* (-vx1-vx2    );        //six
-        //(D.f[dPM0])[kse  ] =   three*c1over54* ( vx1-vx2    );        //six
-        //(D.f[dMP0])[knw  ] =   three*c1over54* (-vx1+vx2    );        //six
-        //(D.f[dP0P])[kte  ] =   three*c1over54* ( vx1    +vx3);        //six
-        //(D.f[dM0M])[kbw  ] =   three*c1over54* (-vx1    -vx3);        //six
-        //(D.f[dP0M])[kbe  ] =   three*c1over54* ( vx1    -vx3);        //six
-        //(D.f[dM0P])[ktw  ] =   three*c1over54* (-vx1    +vx3);        //six
-        //(D.f[d0PP])[ktn  ] =   three*c1over54* (     vx2+vx3);        //six
-        //(D.f[d0MM])[kbs  ] =   three*c1over54* (    -vx2-vx3);        //six
-        //(D.f[d0PM])[kbn  ] =   three*c1over54* (     vx2-vx3);        //six
-        //(D.f[d0MP])[kts  ] =   three*c1over54* (    -vx2+vx3);        //six
-        //(D.f[dPPP])[ktne ] =   three*c1over216*( vx1+vx2+vx3);        //six
-        //(D.f[dMMM])[kbsw ] =   three*c1over216*(-vx1-vx2-vx3);        //six
-        //(D.f[dPPM])[kbne ] =   three*c1over216*( vx1+vx2-vx3);        //six
-        //(D.f[dMMP])[ktsw ] =   three*c1over216*(-vx1-vx2+vx3);        //six
-        //(D.f[dPMP])[ktse ] =   three*c1over216*( vx1-vx2+vx3);        //six
-        //(D.f[dMPM])[kbnw ] =   three*c1over216*(-vx1+vx2-vx3);        //six
-        //(D.f[dPMM])[kbse ] =   three*c1over216*( vx1-vx2-vx3);        //six
-        //(D.f[dMPP])[ktnw ] =   three*c1over216*(-vx1+vx2+vx3);        //six
-        //(D.f[d000])[kzero] =   c8over27* (drho-cu_sq);
-        //(D.f[dP00])[ke   ] =   c2over27* (drho+three*( vx1        )+c9over2*( vx1        )*( vx1        )-cu_sq);
-        //(D.f[dM00])[kw   ] =   c2over27* (drho+three*(-vx1        )+c9over2*(-vx1        )*(-vx1        )-cu_sq);
-        //(D.f[d0P0])[kn   ] =   c2over27* (drho+three*(    vx2     )+c9over2*(     vx2    )*(     vx2    )-cu_sq);
-        //(D.f[d0M0])[ks   ] =   c2over27* (drho+three*(   -vx2     )+c9over2*(    -vx2    )*(    -vx2    )-cu_sq);
-        //(D.f[d00P])[kt   ] =   c2over27* (drho+three*(         vx3)+c9over2*(         vx3)*(         vx3)-cu_sq);
-        //(D.f[d00M])[kb   ] =   c2over27* (drho+three*(        -vx3)+c9over2*(        -vx3)*(        -vx3)-cu_sq);
-        //(D.f[dPP0])[kne  ] =   c1over54* (drho+three*( vx1+vx2    )+c9over2*( vx1+vx2    )*( vx1+vx2    )-cu_sq);
-        //(D.f[dMM0])[ksw  ] =   c1over54* (drho+three*(-vx1-vx2    )+c9over2*(-vx1-vx2    )*(-vx1-vx2    )-cu_sq);
-        //(D.f[dPM0])[kse  ] =   c1over54* (drho+three*( vx1-vx2    )+c9over2*( vx1-vx2    )*( vx1-vx2    )-cu_sq);
-        //(D.f[dMP0])[knw  ] =   c1over54* (drho+three*(-vx1+vx2    )+c9over2*(-vx1+vx2    )*(-vx1+vx2    )-cu_sq);
-        //(D.f[dP0P])[kte  ] =   c1over54* (drho+three*( vx1    +vx3)+c9over2*( vx1    +vx3)*( vx1    +vx3)-cu_sq);
-        //(D.f[dM0M])[kbw  ] =   c1over54* (drho+three*(-vx1    -vx3)+c9over2*(-vx1    -vx3)*(-vx1    -vx3)-cu_sq);
-        //(D.f[dP0M])[kbe  ] =   c1over54* (drho+three*( vx1    -vx3)+c9over2*( vx1    -vx3)*( vx1    -vx3)-cu_sq);
-        //(D.f[dM0P])[ktw  ] =   c1over54* (drho+three*(-vx1    +vx3)+c9over2*(-vx1    +vx3)*(-vx1    +vx3)-cu_sq);
-        //(D.f[d0PP])[ktn  ] =   c1over54* (drho+three*(     vx2+vx3)+c9over2*(     vx2+vx3)*(     vx2+vx3)-cu_sq);
-        //(D.f[d0MM])[kbs  ] =   c1over54* (drho+three*(    -vx2-vx3)+c9over2*(    -vx2-vx3)*(    -vx2-vx3)-cu_sq);
-        //(D.f[d0PM])[kbn  ] =   c1over54* (drho+three*(     vx2-vx3)+c9over2*(     vx2-vx3)*(     vx2-vx3)-cu_sq);
-        //(D.f[d0MP])[kts  ] =   c1over54* (drho+three*(    -vx2+vx3)+c9over2*(    -vx2+vx3)*(    -vx2+vx3)-cu_sq);
-        //(D.f[dPPP])[ktne ] =   c1over216*(drho+three*( vx1+vx2+vx3)+c9over2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
-        //(D.f[dMMM])[kbsw ] =   c1over216*(drho+three*(-vx1-vx2-vx3)+c9over2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
-        //(D.f[dPPM])[kbne ] =   c1over216*(drho+three*( vx1+vx2-vx3)+c9over2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
-        //(D.f[dMMP])[ktsw ] =   c1over216*(drho+three*(-vx1-vx2+vx3)+c9over2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
-        //(D.f[dPMP])[ktse ] =   c1over216*(drho+three*( vx1-vx2+vx3)+c9over2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
-        //(D.f[dMPM])[kbnw ] =   c1over216*(drho+three*(-vx1+vx2-vx3)+c9over2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
-        //(D.f[dPMM])[kbse ] =   c1over216*(drho+three*( vx1-vx2-vx3)+c9over2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
-        //(D.f[dMPP])[ktnw ] =   c1over216*(drho+three*(-vx1+vx2+vx3)+c9over2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
-        }
-    }
-}
-//////////////////////////////////////////////////////////////////////////
-
-
diff --git a/src/gpu/core/Parameter/Parameter.cpp b/src/gpu/core/Parameter/Parameter.cpp
index 4177aa2fc9bb5475f025700c8217a36ab4387eca..14aacb94a6da732d581b9820b2d34390b1472aff 100644
--- a/src/gpu/core/Parameter/Parameter.cpp
+++ b/src/gpu/core/Parameter/Parameter.cpp
@@ -397,7 +397,6 @@ void Parameter::initGridPaths(){
     this->setnumberNodes(gridPath + "numberNodes.dat");
     this->setLBMvsSI(gridPath + "LBMvsSI.dat");
     this->setmeasurePoints(gridPath + "measurePoints.dat");
-    this->setpropellerValues(gridPath + "propellerValues.dat");
     this->setcpTop(gridPath + "cpTop.dat");
     this->setcpBottom(gridPath + "cpBottom.dat");
     this->setcpBottom2(gridPath + "cpBottom2.dat");
@@ -957,10 +956,6 @@ void Parameter::setIsOutflowNormal(bool isOutflowNormal)
 {
     this->isOutflowNormal = isOutflowNormal;
 }
-void Parameter::setIsProp(bool isProp)
-{
-    this->isProp = isProp;
-}
 void Parameter::setIsCp(bool isCp)
 {
     this->isCp = isCp;
@@ -1288,18 +1283,6 @@ void Parameter::setperiodicBcValues(std::string periodicBcValues)
 {
     this->periodicBcValues = periodicBcValues;
 }
-void Parameter::setpropellerQs(std::string propellerQs)
-{
-    this->propellerQs = propellerQs;
-}
-void Parameter::setpropellerValues(std::string propellerValues)
-{
-    this->propellerValues = propellerValues;
-}
-void Parameter::setpropellerCylinder(std::string propellerCylinder)
-{
-    this->propellerCylinder = propellerCylinder;
-}
 void Parameter::setmeasurePoints(std::string measurePoints)
 {
     this->measurePoints = measurePoints;
@@ -1356,8 +1339,6 @@ void Parameter::setObj(std::string str, bool isObj)
 {
     if (str == "geo") {
         this->setIsGeo(isObj);
-    } else if (str == "prop") {
-        this->setIsProp(isObj);
     } else if (str == "cp") {
         this->setIsCp(isObj);
     } else if (str == "geoNormal") {
@@ -2246,18 +2227,6 @@ std::string Parameter::getperiodicBcValues()
 {
     return this->periodicBcValues;
 }
-std::string Parameter::getpropellerQs()
-{
-    return this->propellerQs;
-}
-std::string Parameter::getpropellerValues()
-{
-    return this->propellerValues;
-}
-std::string Parameter::getpropellerCylinder()
-{
-    return this->propellerCylinder;
-}
 std::string Parameter::getmeasurePoints()
 {
     return this->measurePoints;
@@ -2409,10 +2378,6 @@ bool Parameter::getCalcHighOrderMoments()
 {
     return this->isHighOrderMoments;
 }
-bool Parameter::getIsProp()
-{
-    return this->isProp;
-}
 bool Parameter::overWritingRestart(uint t)
 {
     return t == getTimeDoRestart();
diff --git a/src/gpu/core/Parameter/Parameter.h b/src/gpu/core/Parameter/Parameter.h
index cc062e2d96b5e049c7a2e692860da026e49ab5b5..35c638878fbd235480132e363df0df0e8e858c08 100644
--- a/src/gpu/core/Parameter/Parameter.h
+++ b/src/gpu/core/Parameter/Parameter.h
@@ -373,7 +373,6 @@ struct LBMSimulationParameter {
     QforBoundaryConditions QInlet, QOutlet, QPeriodic; // DEPRECATED BCs that are not used any more
     unsigned int kInletQread, kOutletQread;            // DEPRECATED
 
-    QforBoundaryConditions propellerBC;                                                 // DEPRECATED
     QforBoundaryConditions geometryBCnormalX, geometryBCnormalY, geometryBCnormalZ;     // DEPRECATED
     QforBoundaryConditions inflowBCnormalX, inflowBCnormalY, inflowBCnormalZ;           // DEPRECATED
     QforBoundaryConditions outflowBCnormalX, outflowBCnormalY, outflowBCnormalZ;        // DEPRECATED
@@ -571,9 +570,6 @@ public:
     void setwallBcValues(std::string wallBcValues);
     void setperiodicBcQs(std::string periodicBcQs);
     void setperiodicBcValues(std::string periodicBcValues);
-    void setpropellerCylinder(std::string propellerCylinder);
-    void setpropellerValues(std::string propellerValues);
-    void setpropellerQs(std::string propellerQs);
     void setmeasurePoints(std::string measurePoints);
     void setnumberNodes(std::string numberNodes);
     void setLBMvsSI(std::string LBMvsSI);
@@ -599,7 +595,6 @@ public:
     void setIsGeoNormal(bool isGeoNormal);
     void setIsInflowNormal(bool isInflowNormal);
     void setIsOutflowNormal(bool isOutflowNormal);
-    void setIsProp(bool isProp);
     void setIsCp(bool isCp);
     void setConcFile(bool concFile);
     void setUseMeasurePoints(bool useMeasurePoints);
@@ -785,9 +780,6 @@ public:
     std::string getwallBcValues();
     std::string getperiodicBcQs();
     std::string getperiodicBcValues();
-    std::string getpropellerQs();
-    std::string getpropellerCylinder();
-    std::string getpropellerValues();
     std::string getmeasurePoints();
     std::string getnumberNodes();
     std::string getLBMvsSI();
@@ -885,7 +877,6 @@ public:
     bool getIsGeoNormal();
     bool getIsInflowNormal();
     bool getIsOutflowNormal();
-    bool getIsProp();
     bool getIsCp();
     bool getIsGeometryValues();
     bool getCalc2ndOrderMoments();
@@ -1003,7 +994,6 @@ private:
     bool doCheckPoint{ false };
     bool readGeo{ false };
     bool isGeo;
-    bool isProp;
     bool isCp;
     bool GeometryValues{ false };
     bool is2ndOrderMoments{ false };
@@ -1060,7 +1050,7 @@ private:
     std::string pressBcPos, pressBcQs, pressBcValue;
     std::string geomBoundaryBcQs, velBcQs;
     std::string geomBoundaryBcValues, velBcValues, pressBcValues, noSlipBcValues;
-    std::string propellerCylinder, propellerValues, propellerQs, measurePoints;
+    std::string measurePoints;
     std::string inletBcQs, inletBcValues;
     std::string outletBcQs, outletBcValues;
     std::string topBcQs, topBcValues;