diff --git a/src/gpu/VirtualFluids_GPU/GPU/PressBCs27.cu b/src/gpu/VirtualFluids_GPU/GPU/PressBCs27.cu
index e105a1a08219f2c6e0009e994dc1f83af5f57b33..6a14cebd465a1e79acf14ea019abb34e66d9c85f 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/PressBCs27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/PressBCs27.cu
@@ -3057,124 +3057,86 @@ __global__ void QPressZeroRhoOutflowDevice27(  real* rhoBC,
    if(k>=numberOfBCnodes) return;
    ////////////////////////////////////////////////////////////////////////////////
    //index
-   unsigned int KQK  = k_Q[k];
-   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];
+
+   uint k_000 = k_Q[k];
+   uint k_M00 = neighborX[k_000];
+   uint k_0M0 = neighborY[k_000];
+   uint k_00M = neighborZ[k_000];
+   uint k_MM0 = neighborY[k_M00];
+   uint k_M0M = neighborZ[k_M00];
+   uint k_0MM = neighborZ[k_0M0];
+   uint k_MMM = neighborZ[k_MM0];
+
    ////////////////////////////////////////////////////////////////////////////////
-   //index1
-   unsigned int K1QK  = k_N[k];
-   // unsigned int k1zero= K1QK;
-   unsigned int k1e   = K1QK;
-   unsigned int k1w   = neighborX[K1QK];
-   unsigned int k1n   = K1QK;
-   unsigned int k1s   = neighborY[K1QK];
-   unsigned int k1t   = K1QK;
-   unsigned int k1b   = neighborZ[K1QK];
-   unsigned int k1sw  = neighborY[k1w];
-   unsigned int k1ne  = K1QK;
-   unsigned int k1se  = k1s;
-   unsigned int k1nw  = k1w;
-   unsigned int k1bw  = neighborZ[k1w];
-   unsigned int k1te  = K1QK;
-   unsigned int k1be  = k1b;
-   unsigned int k1tw  = k1w;
-   unsigned int k1bs  = neighborZ[k1s];
-   unsigned int k1tn  = K1QK;
-   unsigned int k1bn  = k1b;
-   unsigned int k1ts  = k1s;
-   unsigned int k1tse = k1s;
-   unsigned int k1bnw = k1bw;
-   unsigned int k1tnw = k1w;
-   unsigned int k1bse = k1bs;
-   unsigned int k1tsw = k1sw;
-   unsigned int k1bne = k1b;
-   unsigned int k1tne = K1QK;
-   unsigned int k1bsw = neighborZ[k1sw];
+   //index of neighbor
+   uint kN_000 = k_N[k];
+   uint kN_M00 = neighborX[k_000];
+   uint kN_0M0 = neighborY[k_000];
+   uint kN_00M = neighborZ[k_000];
+   uint kN_MM0 = neighborY[k_M00];
+   uint kN_M0M = neighborZ[k_M00];
+   uint kN_0MM = neighborZ[k_0M0];
+   uint kN_MMM = neighborZ[k_MM0];
    ////////////////////////////////////////////////////////////////////////////////
    Distributions27 dist;
    getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);   
-   real f1[27], f[27];   
+   real f[27], fN[27];   
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-   f1[DIR_P00] = (dist.f[DIR_P00])[k1e   ];
-   f1[DIR_M00] = (dist.f[DIR_M00])[k1w   ];
-   f1[DIR_0P0] = (dist.f[DIR_0P0])[k1n   ];
-   f1[DIR_0M0] = (dist.f[DIR_0M0])[k1s   ];
-   f1[DIR_00P] = (dist.f[DIR_00P])[k1t   ];
-   f1[DIR_00M] = (dist.f[DIR_00M])[k1b   ];
-   f1[DIR_PP0] = (dist.f[DIR_PP0])[k1ne  ];
-   f1[DIR_MM0] = (dist.f[DIR_MM0])[k1sw  ];
-   f1[DIR_PM0] = (dist.f[DIR_PM0])[k1se  ];
-   f1[DIR_MP0] = (dist.f[DIR_MP0])[k1nw  ];
-   f1[DIR_P0P] = (dist.f[DIR_P0P])[k1te  ];
-   f1[DIR_M0M] = (dist.f[DIR_M0M])[k1bw  ];
-   f1[DIR_P0M] = (dist.f[DIR_P0M])[k1be  ];
-   f1[DIR_M0P] = (dist.f[DIR_M0P])[k1tw  ];
-   f1[DIR_0PP] = (dist.f[DIR_0PP])[k1tn  ];
-   f1[DIR_0MM] = (dist.f[DIR_0MM])[k1bs  ];
-   f1[DIR_0PM] = (dist.f[DIR_0PM])[k1bn  ];
-   f1[DIR_0MP] = (dist.f[DIR_0MP])[k1ts  ];
-   // f1[DIR_000] = (dist.f[DIR_000])[k1zero];
-   f1[DIR_PPP] = (dist.f[DIR_PPP])[k1tne ];
-   f1[DIR_MMP] = (dist.f[DIR_MMP])[k1tsw ];
-   f1[DIR_PMP] = (dist.f[DIR_PMP])[k1tse ];
-   f1[DIR_MPP] = (dist.f[DIR_MPP])[k1tnw ];
-   f1[DIR_PPM] = (dist.f[DIR_PPM])[k1bne ];
-   f1[DIR_MMM] = (dist.f[DIR_MMM])[k1bsw ];
-   f1[DIR_PMM] = (dist.f[DIR_PMM])[k1bse ];
-   f1[DIR_MPM] = (dist.f[DIR_MPM])[k1bnw ];
+   f[DIR_000] = (dist.f[DIR_000])[k_000];
+   f[DIR_P00] = (dist.f[DIR_P00])[k_000];
+   f[DIR_M00] = (dist.f[DIR_M00])[k_M00];
+   f[DIR_0P0] = (dist.f[DIR_0P0])[k_000];
+   f[DIR_0M0] = (dist.f[DIR_0M0])[k_0M0];
+   f[DIR_00P] = (dist.f[DIR_00P])[k_000];
+   f[DIR_00M] = (dist.f[DIR_00M])[k_00M];
+   f[DIR_PP0] = (dist.f[DIR_PP0])[k_000];
+   f[DIR_MM0] = (dist.f[DIR_MM0])[k_MM0];
+   f[DIR_PM0] = (dist.f[DIR_PM0])[k_0M0];
+   f[DIR_MP0] = (dist.f[DIR_MP0])[k_M00];
+   f[DIR_P0P] = (dist.f[DIR_P0P])[k_000];
+   f[DIR_M0M] = (dist.f[DIR_M0M])[k_M0M];
+   f[DIR_P0M] = (dist.f[DIR_P0M])[k_00M];
+   f[DIR_M0P] = (dist.f[DIR_M0P])[k_M00];
+   f[DIR_0PP] = (dist.f[DIR_0PP])[k_000];
+   f[DIR_0MM] = (dist.f[DIR_0MM])[k_0MM];
+   f[DIR_0PM] = (dist.f[DIR_0PM])[k_00M];
+   f[DIR_0MP] = (dist.f[DIR_0MP])[k_0M0];
+   f[DIR_PPP] = (dist.f[DIR_PPP])[k_000];
+   f[DIR_MPP] = (dist.f[DIR_MPP])[k_M00];
+   f[DIR_PMP] = (dist.f[DIR_PMP])[k_0M0];
+   f[DIR_MMP] = (dist.f[DIR_MMP])[k_MM0];
+   f[DIR_PPM] = (dist.f[DIR_PPM])[k_00M];
+   f[DIR_MPM] = (dist.f[DIR_MPM])[k_M0M];
+   f[DIR_PMM] = (dist.f[DIR_PMM])[k_0MM];
+   f[DIR_MMM] = (dist.f[DIR_MMM])[k_MMM];
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-   f[DIR_P00] = (dist.f[DIR_P00])[ke   ];
-   f[DIR_M00] = (dist.f[DIR_M00])[kw   ];
-   f[DIR_0P0] = (dist.f[DIR_0P0])[kn   ];
-   f[DIR_0M0] = (dist.f[DIR_0M0])[ks   ];
-   f[DIR_00P] = (dist.f[DIR_00P])[kt   ];
-   f[DIR_00M] = (dist.f[DIR_00M])[kb   ];
-   f[DIR_PP0] = (dist.f[DIR_PP0])[kne  ];
-   f[DIR_MM0] = (dist.f[DIR_MM0])[ksw  ];
-   f[DIR_PM0] = (dist.f[DIR_PM0])[kse  ];
-   f[DIR_MP0] = (dist.f[DIR_MP0])[knw  ];
-   f[DIR_P0P] = (dist.f[DIR_P0P])[kte  ];
-   f[DIR_M0M] = (dist.f[DIR_M0M])[kbw  ];
-   f[DIR_P0M] = (dist.f[DIR_P0M])[kbe  ];
-   f[DIR_M0P] = (dist.f[DIR_M0P])[ktw  ];
-   f[DIR_0PP] = (dist.f[DIR_0PP])[ktn  ];
-   f[DIR_0MM] = (dist.f[DIR_0MM])[kbs  ];
-   f[DIR_0PM] = (dist.f[DIR_0PM])[kbn  ];
-   f[DIR_0MP] = (dist.f[DIR_0MP])[kts  ];
-   f[DIR_000] = (dist.f[DIR_000])[kzero];
-   f[DIR_PPP] = (dist.f[DIR_PPP])[ktne ];
-   f[DIR_MMP] = (dist.f[DIR_MMP])[ktsw ];
-   f[DIR_PMP] = (dist.f[DIR_PMP])[ktse ];
-   f[DIR_MPP] = (dist.f[DIR_MPP])[ktnw ];
-   f[DIR_PPM] = (dist.f[DIR_PPM])[kbne ];
-   f[DIR_MMM] = (dist.f[DIR_MMM])[kbsw ];
-   f[DIR_PMM] = (dist.f[DIR_PMM])[kbse ];
-   f[DIR_MPM] = (dist.f[DIR_MPM])[kbnw ];
+   fN[DIR_000] = (dist.f[DIR_000])[kN_000];
+   fN[DIR_P00] = (dist.f[DIR_P00])[kN_000];
+   fN[DIR_M00] = (dist.f[DIR_M00])[kN_M00];
+   fN[DIR_0P0] = (dist.f[DIR_0P0])[kN_000];
+   fN[DIR_0M0] = (dist.f[DIR_0M0])[kN_0M0];
+   fN[DIR_00P] = (dist.f[DIR_00P])[kN_000];
+   fN[DIR_00M] = (dist.f[DIR_00M])[kN_00M];
+   fN[DIR_PP0] = (dist.f[DIR_PP0])[kN_000];
+   fN[DIR_MM0] = (dist.f[DIR_MM0])[kN_MM0];
+   fN[DIR_PM0] = (dist.f[DIR_PM0])[kN_0M0];
+   fN[DIR_MP0] = (dist.f[DIR_MP0])[kN_M00];
+   fN[DIR_P0P] = (dist.f[DIR_P0P])[kN_000];
+   fN[DIR_M0M] = (dist.f[DIR_M0M])[kN_M0M];
+   fN[DIR_P0M] = (dist.f[DIR_P0M])[kN_00M];
+   fN[DIR_M0P] = (dist.f[DIR_M0P])[kN_M00];
+   fN[DIR_0PP] = (dist.f[DIR_0PP])[kN_000];
+   fN[DIR_0MM] = (dist.f[DIR_0MM])[kN_0MM];
+   fN[DIR_0PM] = (dist.f[DIR_0PM])[kN_00M];
+   fN[DIR_0MP] = (dist.f[DIR_0MP])[kN_0M0];
+   fN[DIR_PPP] = (dist.f[DIR_PPP])[kN_000];
+   fN[DIR_MPP] = (dist.f[DIR_MPP])[kN_M00];
+   fN[DIR_PMP] = (dist.f[DIR_PMP])[kN_0M0];
+   fN[DIR_MMP] = (dist.f[DIR_MMP])[kN_MM0];
+   fN[DIR_PPM] = (dist.f[DIR_PPM])[kN_00M];
+   fN[DIR_MPM] = (dist.f[DIR_MPM])[kN_M0M];
+   fN[DIR_PMM] = (dist.f[DIR_PMM])[kN_0MM];
+   fN[DIR_MMM] = (dist.f[DIR_MMM])[kN_MMM];
    //////////////////////////////////////////////////////////////////////////
    real drho = vf::lbm::getDensity(f);
    
@@ -3187,75 +3149,75 @@ __global__ void QPressZeroRhoOutflowDevice27(  real* rhoBC,
    switch(direction)
    {
       case MZZ:
-         (dist.f[DIR_P00])[ke   ] = computeOutflowDistribution(f, f1, DIR_P00  , rhoCorrection, cs, c2o27);
-         (dist.f[DIR_PM0])[kse  ] = computeOutflowDistribution(f, f1, DIR_PM0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_PP0])[kne  ] = computeOutflowDistribution(f, f1, DIR_PP0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_P0M])[kbe  ] = computeOutflowDistribution(f, f1, DIR_P0M, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_P0P])[kte  ] = computeOutflowDistribution(f, f1, DIR_P0P, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_PMP])[ktse ] = computeOutflowDistribution(f, f1, DIR_PMP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_PPP])[ktne ] = computeOutflowDistribution(f, f1, DIR_PPP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_PMM])[kbse ] = computeOutflowDistribution(f, f1, DIR_PMM, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_PPM])[kbne ] = computeOutflowDistribution(f, f1, DIR_PPM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_P00])[k_000] = computeOutflowDistribution(f, fN, DIR_P00  , rhoCorrection, cs, c2o27);
+         (dist.f[DIR_PM0])[k_0M0] = computeOutflowDistribution(f, fN, DIR_PM0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_PP0])[k_000] = computeOutflowDistribution(f, fN, DIR_PP0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_P0M])[k_00M] = computeOutflowDistribution(f, fN, DIR_P0M, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_P0P])[k_000] = computeOutflowDistribution(f, fN, DIR_P0P, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_PMP])[k_0M0] = computeOutflowDistribution(f, fN, DIR_PMP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_PPP])[k_000] = computeOutflowDistribution(f, fN, DIR_PPP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_PMM])[k_0MM] = computeOutflowDistribution(f, fN, DIR_PMM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_PPM])[k_00M] = computeOutflowDistribution(f, fN, DIR_PPM, rhoCorrection, cs, c1o216);
          break;
 
       case PZZ:
-         (dist.f[DIR_M00])[kw   ] = computeOutflowDistribution(f, f1, DIR_M00, rhoCorrection, cs, c2o27);
-         (dist.f[DIR_MM0])[ksw  ] = computeOutflowDistribution(f, f1, DIR_MM0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_MP0])[knw  ] = computeOutflowDistribution(f, f1, DIR_MP0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_M0M])[kbw  ] = computeOutflowDistribution(f, f1, DIR_M0M, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_M0P])[ktw  ] = computeOutflowDistribution(f, f1, DIR_M0P, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_MMP])[ktsw ] = computeOutflowDistribution(f, f1, DIR_MMP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MPP])[ktnw ] = computeOutflowDistribution(f, f1, DIR_MPP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MMM])[kbsw ] = computeOutflowDistribution(f, f1, DIR_MMM, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MPM])[kbnw ] = computeOutflowDistribution(f, f1, DIR_MPM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_M00])[k_M00] = computeOutflowDistribution(f, fN, DIR_M00, rhoCorrection, cs, c2o27);
+         (dist.f[DIR_MM0])[k_MM0] = computeOutflowDistribution(f, fN, DIR_MM0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_MP0])[k_M00] = computeOutflowDistribution(f, fN, DIR_MP0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_M0M])[k_M0M] = computeOutflowDistribution(f, fN, DIR_M0M, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_M0P])[k_M00] = computeOutflowDistribution(f, fN, DIR_M0P, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_MMP])[k_MM0] = computeOutflowDistribution(f, fN, DIR_MMP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MPP])[k_M00] = computeOutflowDistribution(f, fN, DIR_MPP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MMM])[k_MMM] = computeOutflowDistribution(f, fN, DIR_MMM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MPM])[k_M0M] = computeOutflowDistribution(f, fN, DIR_MPM, rhoCorrection, cs, c1o216);
          break;
 
       case ZMZ:
-         (dist.f[DIR_0P0])[kn   ] = computeOutflowDistribution(f, f1, DIR_0P0, rhoCorrection, cs, c2o27);
-         (dist.f[DIR_PP0])[kne  ] = computeOutflowDistribution(f, f1, DIR_PP0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_MP0])[knw  ] = computeOutflowDistribution(f, f1, DIR_MP0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0PP])[ktn  ] = computeOutflowDistribution(f, f1, DIR_0PP, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0PM])[kbn  ] = computeOutflowDistribution(f, f1, DIR_0PM, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_PPP])[ktne ] = computeOutflowDistribution(f, f1, DIR_PPP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MPP])[ktnw ] = computeOutflowDistribution(f, f1, DIR_MPP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_PPM])[kbne ] = computeOutflowDistribution(f, f1, DIR_PPM, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MPM])[kbnw ] = computeOutflowDistribution(f, f1, DIR_MPM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_0P0])[k_000] = computeOutflowDistribution(f, fN, DIR_0P0, rhoCorrection, cs, c2o27);
+         (dist.f[DIR_PP0])[k_000] = computeOutflowDistribution(f, fN, DIR_PP0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_MP0])[k_M00] = computeOutflowDistribution(f, fN, DIR_MP0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0PP])[k_000] = computeOutflowDistribution(f, fN, DIR_0PP, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0PM])[k_00M] = computeOutflowDistribution(f, fN, DIR_0PM, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_PPP])[k_000] = computeOutflowDistribution(f, fN, DIR_PPP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MPP])[k_M00] = computeOutflowDistribution(f, fN, DIR_MPP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_PPM])[k_00M] = computeOutflowDistribution(f, fN, DIR_PPM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MPM])[k_M0M] = computeOutflowDistribution(f, fN, DIR_MPM, rhoCorrection, cs, c1o216);
          break;  
 
       case ZPZ:   
-         (dist.f[DIR_0M0])[ks   ] =computeOutflowDistribution(f, f1, DIR_0M0, rhoCorrection, cs, c2o27);
-         (dist.f[DIR_PM0])[kse  ] =computeOutflowDistribution(f, f1, DIR_PM0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_MM0])[ksw  ] =computeOutflowDistribution(f, f1, DIR_MM0, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0MP])[kts  ] =computeOutflowDistribution(f, f1, DIR_0MP, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0MM])[kbs  ] =computeOutflowDistribution(f, f1, DIR_0MM, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_PMP])[ktse ] =computeOutflowDistribution(f, f1, DIR_PMP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MMP])[ktsw ] =computeOutflowDistribution(f, f1, DIR_MMP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_PMM])[kbse ] =computeOutflowDistribution(f, f1, DIR_PMM, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MMM])[kbsw ] =computeOutflowDistribution(f, f1, DIR_MMM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_0M0])[k_0M0] =computeOutflowDistribution(f, fN, DIR_0M0, rhoCorrection, cs, c2o27);
+         (dist.f[DIR_PM0])[k_0M0] =computeOutflowDistribution(f, fN, DIR_PM0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_MM0])[k_MM0] =computeOutflowDistribution(f, fN, DIR_MM0, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0MP])[k_0M0] =computeOutflowDistribution(f, fN, DIR_0MP, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0MM])[k_0MM] =computeOutflowDistribution(f, fN, DIR_0MM, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_PMP])[k_0M0] =computeOutflowDistribution(f, fN, DIR_PMP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MMP])[k_MM0] =computeOutflowDistribution(f, fN, DIR_MMP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_PMM])[k_0MM] =computeOutflowDistribution(f, fN, DIR_PMM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MMM])[k_MMM] =computeOutflowDistribution(f, fN, DIR_MMM, rhoCorrection, cs, c1o216);
          break;
 
       case ZZM:
-         (dist.f[DIR_00P])[kt   ] = computeOutflowDistribution(f, f1, DIR_00P, rhoCorrection, cs, c2o27);
-         (dist.f[DIR_P0P])[kte  ] = computeOutflowDistribution(f, f1, DIR_P0P, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_M0P])[ktw  ] = computeOutflowDistribution(f, f1, DIR_M0P, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0PP])[ktn  ] = computeOutflowDistribution(f, f1, DIR_0PP, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0MP])[kts  ] = computeOutflowDistribution(f, f1, DIR_0MP, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_PPP])[ktne ] = computeOutflowDistribution(f, f1, DIR_PPP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MPP])[ktnw ] = computeOutflowDistribution(f, f1, DIR_MPP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_PMP])[ktse ] = computeOutflowDistribution(f, f1, DIR_PMP, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MMP])[ktsw ] = computeOutflowDistribution(f, f1, DIR_MMP, rhoCorrection, cs, c1o216); 
+         (dist.f[DIR_00P])[k_000] = computeOutflowDistribution(f, fN, DIR_00P, rhoCorrection, cs, c2o27);
+         (dist.f[DIR_P0P])[k_000] = computeOutflowDistribution(f, fN, DIR_P0P, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_M0P])[k_M00] = computeOutflowDistribution(f, fN, DIR_M0P, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0PP])[k_000] = computeOutflowDistribution(f, fN, DIR_0PP, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0MP])[k_0M0] = computeOutflowDistribution(f, fN, DIR_0MP, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_PPP])[k_000] = computeOutflowDistribution(f, fN, DIR_PPP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MPP])[k_M00] = computeOutflowDistribution(f, fN, DIR_MPP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_PMP])[k_0M0] = computeOutflowDistribution(f, fN, DIR_PMP, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MMP])[k_MM0] = computeOutflowDistribution(f, fN, DIR_MMP, rhoCorrection, cs, c1o216); 
          break;
 
       case ZZP:
-         (dist.f[DIR_00M])[kb   ] = computeOutflowDistribution(f, f1, DIR_00M, rhoCorrection, cs, c2o27);
-         (dist.f[DIR_P0M])[kbe  ] = computeOutflowDistribution(f, f1, DIR_P0M, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_M0M])[kbw  ] = computeOutflowDistribution(f, f1, DIR_M0M, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0PM])[kbn  ] = computeOutflowDistribution(f, f1, DIR_0PM, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_0MM])[kbs  ] = computeOutflowDistribution(f, f1, DIR_0MM, rhoCorrection, cs, c1o54);
-         (dist.f[DIR_PPM])[kbne ] = computeOutflowDistribution(f, f1, DIR_PPM, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MPM])[kbnw ] = computeOutflowDistribution(f, f1, DIR_MPM, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_PMM])[kbse ] = computeOutflowDistribution(f, f1, DIR_PMM, rhoCorrection, cs, c1o216);
-         (dist.f[DIR_MMM])[kbsw ] = computeOutflowDistribution(f, f1, DIR_MMM, rhoCorrection, cs, c1o216);     
+         (dist.f[DIR_00M])[k_00M] = computeOutflowDistribution(f, fN, DIR_00M, rhoCorrection, cs, c2o27);
+         (dist.f[DIR_P0M])[k_00M] = computeOutflowDistribution(f, fN, DIR_P0M, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_M0M])[k_M0M] = computeOutflowDistribution(f, fN, DIR_M0M, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0PM])[k_00M] = computeOutflowDistribution(f, fN, DIR_0PM, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_0MM])[k_0MM] = computeOutflowDistribution(f, fN, DIR_0MM, rhoCorrection, cs, c1o54);
+         (dist.f[DIR_PPM])[k_00M] = computeOutflowDistribution(f, fN, DIR_PPM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MPM])[k_M0M] = computeOutflowDistribution(f, fN, DIR_MPM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_PMM])[k_0MM] = computeOutflowDistribution(f, fN, DIR_PMM, rhoCorrection, cs, c1o216);
+         (dist.f[DIR_MMM])[k_MMM] = computeOutflowDistribution(f, fN, DIR_MMM, rhoCorrection, cs, c1o216);     
          break;
       default:
          break;