88:             double xcos_dxi = 2.0 * xcos(d) / dx(i);

89:             double ycos_dyj = 2.0 * ycos(d) / dy(j);

90:             double zcos_dzk = 2.0 * zcos(d) / dz(k);

91: 

92:             Zone z(zone_layout(*k, *j, *i));

93: 

94:             /* Calculate new zonal flux */

95:             double psi_d_g_z = (rhs(d,g,z)

96:                 + psi_lf(d, g, j, k) * xcos_dxi

97:                 + psi_fr(d, g, i, k) * ycos_dyj

98:                 + psi_bo(d, g, i, j) * zcos_dzk)

99:                 / (xcos_dxi + ycos_dyj + zcos_dzk + sigt(g, z));

100: 

101:             psi(d, g, z) = psi_d_g_z;

102: 

103:             /* Apply diamond-difference relationships */

104:             psi_lf(d, g, j, k) = 2.0 * psi_d_g_z - psi_lf(d, g, j, k);

105:             psi_fr(d, g, i, k) = 2.0 * psi_d_g_z - psi_fr(d, g, i, k);

106:             psi_bo(d, g, i, j) = 2.0 * psi_d_g_z - psi_bo(d, g, i, j);