Real dRho1 = rho[i][t] - rho[i-1][t];
Real dU1 = u[i][t] - u[i-1][t];
Real dP1 = p[i][t] - p[i-1][t];
Real dRho2 = rho[i+1][t] - rho[i][t];
Real dU2 = u[i+1][t] - u[i][t];
Real dP2 = p[i+1][t] - p[i][t];
dRho[i] = max( 0.0, min( parameter*dRho2, max( dRho1, min( dRho2, parameter*dRho1 ))))
+min( 0.0, max( parameter*dRho2, min( dRho1, max( dRho2, parameter*dRho1 ))));
dU[i] = max( 0.0, min( parameter*dU2, max( dU1, min( dU2, parameter*dU1 ))))
+min( 0.0, max( parameter*dU2, min( dU1, max( dU2, parameter*dU1 ))));
dP[i] = max( 0.0, min( parameter*dP2, max( dP1, min( dP2, parameter*dP1 ))))
+min( 0.0, max( parameter*dP2, min( dP1, max( dP2, parameter*dP1 ))));