K-epsilon Source term
 

hi every body.
I have trouble in descritizing of Nonlinear terms of source in k and epsilon equations. My cfd code is based on finite volume method. For example in (du/dy)**2 how should I descritize it? It is computed in center of each cell so I face like this:
(du/dy)**2~|du/dy|*|du/dy|
that using Green theory in each cell:
|du/dy|={-[u(i)*dx(i)]}/area of cell
but when I implement, the Fortran code is terminated by math (M6101) error. tanx.