114:           if (node_flux(i, j) < 0.0) {

115:             upwind = i + 2;

116:             donor = i + 1;

117:             downwind = i;

118:             dif = donor;

119:           } else {

120:             upwind = i - 1;

121:             donor = i;

122:             downwind = i + 1;

123:             dif = upwind;

124:           }

125:           sigma = std::fabs(node_flux(i, j)) / (node_mass_pre(donor, j));

126:           width = celldx[i];

127:           vdiffuw = vel1(donor, j) - vel1(upwind, j);

128:           vdiffdw = vel1(downwind, j) - vel1(donor, j);

129:           limiter = 0.0;

130:           if (vdiffuw * vdiffdw > 0.0) {

131:             auw = std::fabs(vdiffuw);

132:             adw = std::fabs(vdiffdw);

133:             wind = 1.0;

134:             if (vdiffdw <= 0.0) wind = -1.0;

135:             limiter =

136:                 wind * std::fmin(std::fmin(width * ((2.0 - sigma) * adw / width + (1.0 + sigma) * auw / celldx[dif]) / 6.0, auw), adw);

137:           }

138:           advec_vel_s = vel1(donor, j) + (1.0 - sigma) * limiter;

139:           mom_flux(i, j) = advec_vel_s * node_flux(i, j);