#pragma once

// [bl-flux-def.h]: boundary-layer turbulent fluxes definitions
//
// -------------------------------------------------------------------------------------------- //

#include "wstgrid3d.h"



namespace nse
{
	// U,V,W,P 2nd order moments
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void uv_flux(T* flux,		// node: [C || V]

		const T* const UV,		// node: [C || V]
		const T* const U,		// node: [C || C]
		const T* const V,		// node: [C || V]
		const nse_const3d::axisType axis, const wstGrid3d< T >& grid);	// [axisZ || axisYZ]

	template< typename T >
	void uw_flux(T* flux,		// node: [W]

		const T* const UW,		// node: [W]
		const T* const U,		// node: [C]
		const T* const W,		// node: [W]
		const nse_const3d::axisType axis, const wstGrid3d< T >& grid);	// [axisZ || axisYZ]

	template< typename T >
	void vw_flux(T* flux,		// node: [W || VW]

		const T* const VW,		// node: [W || VW]
		const T* const V,		// node: [C || V]
		const T* const W,		// node: [W || W]
		const nse_const3d::axisType axis, const wstGrid3d< T >& grid);	// [axisZ || axisYZ]


	template< typename T >
	void pu_flux(T* flux,		// node: [C]

		const T* const PU,		// node: [C]
		const T* const P,		// node: [C]
		const T* const U,		// node: [C]
		const wstGrid3d< T >& grid);

	template< typename T >
	void pv_flux(T* flux,		// node: [C]

		const T* const PV,		// node: [C]
		const T* const P,		// node: [C]
		const T* const V,		// node: [C]
		const wstGrid3d< T >& grid);

	template< typename T >
	void pw_flux(T* flux,		// node: [W]

		const T* const PW,		// node: [W]
		const T* const P,		// node: [C]
		const T* const W,		// node: [W]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- /

	// U,V,W,C 2nd order moments
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void cu_flux(T* flux,		// node: [C]

		const T* const CU,		// node: [C]
		const T* const C,		// node: [C]
		const T* const U,		// node: [C]
		const wstGrid3d< T >& grid);

	template< typename T >
	void cv_flux(T* flux,		// node: [C]

		const T* const CV,		// node: [C]
		const T* const C,		// node: [C]
		const T* const V,		// node: [C]
		const wstGrid3d< T >& grid);

	template< typename T >
	void cw_flux(T* flux,		// node: [W]

		const T* const CW,		// node: [W]
		const T* const C,		// node: [C]
		const T* const W,		// node: [W]
		const nse_const3d::axisType axis, const wstGrid3d< T >& grid);	// [axisZ || axisYZ]

	template< typename T >
	void cc_flux(T* flux,		// node: [C]

		const T* const CC,		// node: [C]
		const T* const Ca,		// node: [C]
		const T* const Cb,		// node: [C]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// [U^2,V^2,W^2] * W - fluxes
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void u2w_flux(T* flux,		// node: [W]

		const T* const U2W,		// node: [W]
		const T* const U2,		// node: [W]
		const T* const UW,		// node: [W]
		const T* const UW_adv,	// node: [W]
		const T* const U,	 	// node: [C]
		const T* const W,		// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void v2w_flux(T* flux,		// node: [W]

		const T* const V2W,		// node: [W]
		const T* const V2,		// node: [W]
		const T* const VW,		// node: [W]
		const T* const VW_adv,	// node: [W]
		const T* const V,		// node: [C]
		const T* const W,		// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void w2w_flux(T* flux,		// node: [C]

		const T* const W3,		// node: [C]
		const T* const W2_c,	// node: [C]
		const T* const W2_w,	// node: [W]
		const T* const W,		// node: [W]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// [C^2] * W - fluxes
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void c2w_flux(T* flux,		// node: [W]

		const T* const C2W,		// node: [W]
		const T* const C2,		// node: [W]
		const T* const CW,		// node: [W]
		const T* const CW_adv,	// node: [W]
		const T* const C,		// node: [C]
		const T* const W,		// node: [W]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// [UV, UW, VW] * W - fluxes
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void uvw_flux(T* flux,				// node: [W]

		const T* const UVW,				// node: [W]
		const T* const UW,				// node: [W]
		const T* const VW,				// node: [W]
		const T* const UV_uvw,			// node: [W]
		const T* const UW_uvw,			// node: [W]
		const T* const VW_uvw,			// node: [W]
		const T* const U,				// node: [C]
		const T* const V,				// node: [C]
		const T* const W,				// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void uww_flux(T* flux,				// node: [C]
		
		const T* const UWW,				// node: [C]
		const T* const W2_w,			// node: [W]
		const T* const W2_c,			// node: [C]
		const T* const W2_u,			// node: [C]
		const T* const W2_uw,			// node: [W]
		const T* const UW,				// node: [W]
		const T* const UW_bottom_uw,	// node: [W (C -- W)]
		const T* const UW_top_uw,		// node: [W (C -- W)]
		const T* const U,				// node: [C]
		const T* const W,				// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void vww_flux(T* flux,				// node: [C]

		const T* const VWW,				// node: [C]
		const T* const W2_w,			// node: [W]
		const T* const W2_c,			// node: [C]
		const T* const W2_v,			// node: [C]
		const T* const W2_vw,			// node: [W]
		const T* const VW,				// node: [W]
		const T* const VW_bottom_vw,	// node: [W (C -- W)]
		const T* const VW_top_vw,		// node: [W (C -- W)]
		const T* const V,				// node: [C]
		const T* const W,				// node: [W]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// [CU, CV, CW] * W - fluxes
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void cuw_flux(T* flux,				// node: [W]

		const T* const CUW,				// node: [W]
		const T* const UW,				// node: [W]
		const T* const CU_uw,			// node: [W]
		const T* const CW,				// node: [W]
		const T* const CW_uw,			// node: [W]
		const T* const C,				// node: [C]
		const T* const U,				// node: [C]
		const T* const W,				// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void cvw_flux(T* flux,				// node: [W]

		const T* const CVW,				// node: [W]
		const T* const VW,				// node: [W]
		const T* const CV_vw,			// node: [W]
		const T* const CW,				// node: [W]
		const T* const CW_vw,			// node: [W]
		const T* const C,				// node: [C]
		const T* const V,				// node: [C]
		const T* const W,				// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void cww_flux(T* flux,				// node: [C]

		const T* const CWW,				// node: [C]
		const T* const W2_w,			// node: [W]
		const T* const W2_c,			// node: [C]
		const T* const CW,				// node: [W]
		const T* const CW_bottom_w,		// node: [W (C -- W)]
		const T* const CW_top_w,		// node: [W (C -- W)]
		const T* const C,				// node: [C]
		const T* const W,				// node: [W]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// c - pressure gradient covariances
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void c_w_pressure_gradient_turb(T* c_dpdz_turb,		// node: [W]

		const T* const C_dPdz,							// node: [W]
		const T* const C,								// node: [C]
		const T* const Pressure,						// node: [C]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //
}


// U,V,W,P 2nd order moments
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::uv_flux(
	T* flux,					// node: [C || V]

	const T* const UV,			// node: [C || V]
	const T* const U,			// node: [C || C]
	const T* const V,			// node: [C || V]
	const nse_const3d::axisType axis, const wstGrid3d< T >& grid)	// [axisZ || axisYZ]
	// ____   __   _   _ 
	// u'v' = UV - U * V
{
	if (axis == nse_const3d::axisZ) {

		int k;
#pragma omp parallel for private(k) shared(flux)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			flux[k] = UV[k] - U[k] * V[k];
		}
		return;
	}
	if (axis == nse_const3d::axisYZ) {
		// [UV, V] averages have to be known at all [V] nodes, including walls
		// [U] average has to be known at all [C] nodes and ghost nodes: (j + 1/2), (j - 1/2)

		int j, k, idx;
#pragma omp parallel for private(j, k, idx) shared(flux)
		for (j = grid.gcy; j <= grid.ny - grid.gcy; j++)	// all [V] nodes
			for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
				idx = j * grid.nz + k;
				flux[idx] = UV[idx] - (T) 0.5 * (U[idx] + U[idx - grid.nz]) * V[idx];
			}
		return;
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::uw_flux(
	T* flux,					// node: [W]

	const T* const UW,			// node: [W]
	const T* const U,			// node: [C]
	const T* const W,			// node: [W]
	const nse_const3d::axisType axis, const wstGrid3d< T >& grid)	// [axisZ || axisYZ]
	// ____   __   _   _ 
	// u'w' = UW - U * W
	// [UW, W] averages have to be known at all [W] nodes, including walls
	// [U] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	if (axis == nse_const3d::axisZ) {

		int k;
#pragma omp parallel for private(k) shared(flux)
		for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
			flux[k] = UW[k] - (T) 0.5 * (U[k] + U[k - 1]) * W[k];
		}
		return;
	}
	if (axis == nse_const3d::axisYZ) {

		int j, k, idx;
#pragma omp parallel for private(j, k, idx) shared(flux)
		for (j = grid.gcy; j < grid.ny - grid.gcy; j++)
			for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
				idx = j * grid.nz + k;
				flux[idx] = UW[idx] - (T) 0.5 * (U[idx] + U[idx - 1]) * W[idx];
			}
		return;
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::vw_flux(
	T* flux,					// node: [W || VW]

	const T* const VW,			// node: [W || VW]
	const T* const V,			// node: [C || V]
	const T* const W,			// node: [W || W]

	const nse_const3d::axisType axis, const wstGrid3d< T >& grid)	// [axisZ || axisYZ]
	// ____   __   _   _ 
	// v'w' = VW - V * W
{
	if (axis == nse_const3d::axisZ) {
		// [VW, W] averages have to be known at all [W] nodes, including walls
		// [V] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)

		int k;
#pragma omp parallel for private(k) shared(flux)
		for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
			flux[k] = VW[k] - (T) 0.5 * (V[k] + V[k - 1]) * W[k];
		}
		return;
	}
	if (axis == nse_const3d::axisYZ) {
		// [VW] average has to be known at all [VW] nodes, including walls
		// [W] average has to be known at all [W] nodes and ghost nodes: (j + 1/2), (j - 1/2)
		// [V] average has to be known at all [V] nodes and ghost nodes: (k + 1/2), (k - 1/2)

		int j, k, idx;
#pragma omp parallel for private(j, k, idx) shared(flux)
		for (j = grid.gcy; j <= grid.ny - grid.gcy; j++)	// all [V] nodes
			for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
				idx = j * grid.nz + k;
				flux[idx] = VW[idx] - (T) 0.25 *
					(V[idx] + V[idx - 1]) * (W[idx] + W[idx - grid.nz]);
			}
		return;
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::pu_flux(T* flux,	// node: [C]

	const T* const PU,		// node: [C]
	const T* const P,		// node: [C]
	const T* const U,		// node: [C]
	const wstGrid3d< T >& grid)
	// ____   __   _   _ 
	// p'u' = PU - P * U
	//
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		flux[k] = PU[k] - P[k] * U[k];
	}
}

template< typename T >
void nse::pv_flux(T* flux,	// node: [C]

	const T* const PV,		// node: [C]
	const T* const P,		// node: [C]
	const T* const V,		// node: [C]
	const wstGrid3d< T >& grid)
	// ____   __   _   _ 
	// p'v' = PV - P * V
	//
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		flux[k] = PV[k] - P[k] * V[k];
	}
}

template< typename T >
void nse::pw_flux(
	T* flux,						// node: [W]

	const T* const PW,				// node: [W]
	const T* const P,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	// ____   __   _   _ 
	// p'w' = PW - P * W
	// ____     __ 
	// p'w' and PW - are defined at [W] node for the following reason:
	// for staggered grid interpolation to [C] node in TKE equation results in: 
	//	______ 1z        ___ 1z
	//	    dp     d(w *  p )        dw
	//	w * --  =  --------    - p * --
	//	    dz        dz             dz
	//
	// [PW, W] averages have to be known at all [W] nodes, including walls
	// [P] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		flux[k] = PW[k] - (T) 0.5 * (P[k] + P[k - 1]) * W[k];
	}
}
// -------------------------------------------------------------------------------------------- //

// U,V,W,C 2nd order moments
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cu_flux(
	T* flux,					// node: [C]

	const T* const CU,			// node: [C]
	const T* const C,			// node: [C]
	const T* const U,			// node: [C]
	const wstGrid3d< T >& grid)
	// ____   __   _   _ 
	// c'u' = CU - C * U
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		flux[k] = CU[k] - C[k] * U[k];
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::cv_flux(
	T* flux,					// -z [C] node

	const T* const CV,			// -z [C] node
	const T* const C,			// -z [C] node
	const T* const V,			// -z [C] node
	const wstGrid3d< T >& grid)
	// ____   __   _   _ 
	// c'v' = CV - C * V
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		flux[k] = CV[k] - C[k] * V[k];
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::cw_flux(
	T* flux,					// node: [W]

	const T* const CW,			// node: [W]
	const T* const C,			// node: [C]
	const T* const W,			// node: [W]
	const nse_const3d::axisType axis, const wstGrid3d< T >& grid)	// [axisZ || axisYZ]
	// ____   __   _   _ 
	// c'w' = CW - C * W
	// [CW, W] averages have to be known at all [W] nodes, including walls
	// [C] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	if (axis == nse_const3d::axisZ) {

		int k;
#pragma omp parallel for private(k) shared(flux)
		for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
			flux[k] = CW[k] - (T) 0.5 * (C[k] + C[k - 1]) * W[k];
		}
		return;
	}
	if (axis == nse_const3d::axisYZ) {

		int j, k, idx;
#pragma omp parallel for private(j, k, idx) shared(flux)
		for (j = grid.gcy; j < grid.ny - grid.gcy; j++)
			for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
				idx = j * grid.nz + k;
				flux[idx] = CW[idx] - (T) 0.5 * (C[idx] + C[idx - 1]) * W[idx];
			}
		return;
	}
}
// -------------------------------------------------------------------------------------------- //

// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cc_flux(T* flux,		// node: [C]

	const T* const CC,		// node: [C]
	const T* const Ca,		// node: [C]
	const T* const Cb,		// node: [C]
	const wstGrid3d< T >& grid)
	// ______   __   __   __ 
	// ca'cb' = CC - Cb * Ca
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		flux[k] = CC[k] - Ca[k] * Cb[k];
	}
}
// -------------------------------------------------------------------------------------------- //

// [U^2,V^2,W^2] * W - fluxes
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::u2w_flux(
	T* flux,					// node: [W]

	const T* const U2W,			// node: [W]
	const T* const U2,			// node: [W]
	const T* const UW,			// node: [W]
	const T* const UW_adv,		// node: [W]
	const T* const U,	 		// node: [C]
	const T* const W,			// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _       _   __       _   _   _
	// u'u'w' = UUW - UU * W - 2 * U * UW + 2 * U * U * W
	//                ___     _   _   _
	// calculation of UUW and U * U * W product:
	//	~~~~~1z   __1x
	//  (U * U) * W
	//
	// [U2W, U2, UW, UW_adv, W] averages have to be known at all [W] nodes, including walls
	// [U] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		flux[k] =
			U2W[k] - U2[k] * W[k] - (U[k] + U[k - 1]) * UW[k] +
			(T) 2.0 * (U[k] * U[k - 1]) * W[k] +

			// correction term:
			UW_adv[k] * (U[k] - U[k - 1]) * grid.dzh[k];
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::v2w_flux(
	T* flux,				// node: [W]

	const T* const V2W,		// node: [W]
	const T* const V2,		// node: [W]
	const T* const VW,		// node: [W]
	const T* const VW_adv,	// node: [W]
	const T* const V,		// node: [C]
	const T* const W,		// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _       _   __       _   _   _
	// v'v'w' = VVW - VV * W - 2 * V * VW + 2 * V * V * W
	//                ___			  _   _   _
	// calculation of VVW product and V * V * W:
	//	~~~~~1z   __1y
	//  (V * V) * W
	//
	// [V2W, V2, VW, VW_adv, W] averages have to be known at all [W] nodes, including walls
	// [V] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		flux[k] =
			V2W[k] - V2[k] * W[k] - (V[k] + V[k - 1]) * VW[k] +
			(T) 2.0 * (V[k] * V[k - 1]) * W[k] +

			// correction term:
			VW_adv[k] * (V[k] - V[k - 1]) * grid.dzh[k];
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::w2w_flux(
	T* flux,				// node: [C]

	const T* const W3,		// node: [C]
	const T* const W2_c,	// node: [C]
	const T* const W2_w,	// node: [W]
	const T* const W,		// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___       _   __       _   _   _
	// w'w'w' = WWW - 3 * W * WW + 2 * W * W * W
	//                ___     _   _   _
	// calculation of WWW and W * W * W product:
	//	~~~~~1z   __1z
	//  (W * W) * W
	//
	// [W2_w, W] averages have to be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		flux[k] =
			W3[k] - (W[k] + W[k + 1]) * W2_c[k] -
			(T) 0.5 * (W[k] * W2_w[k + 1] + W[k + 1] * W2_w[k]) +

			W[k] * W[k + 1] * (W[k] + W[k + 1]);
	}
}
// -------------------------------------------------------------------------------------------- //


// [C^2] * W - fluxes
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::c2w_flux(
	T* flux,				// node: [W]

	const T* const C2W,		// node: [W]
	const T* const C2,		// node: [W]
	const T* const CW,		// node: [W]
	const T* const CW_adv,	// node: [W]
	const T* const C,		// node: [C]
	const T* const W,		// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _       _   __       _   _   _
	// c'c'w' = CCW - CC * W - 2 * C * CW + 2 * C * C * W
	//                ___			  _   _   _
	// calculation of CCW product and C * C * W:
	//	~~~~~1z   
	//  (C * C) * W
	//
	// [C2W, C2, CW, CW_adv, W] averages have to be known at all [W] nodes, including walls
	// [C] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		flux[k] =
			C2W[k] - C2[k] * W[k] - (C[k] + C[k - 1]) * CW[k] +
			(T) 2.0 * (C[k] * C[k - 1]) * W[k] +

			// correction term:
			CW_adv[k] * (C[k] - C[k - 1]) * grid.dzh[k];
	}
}
// -------------------------------------------------------------------------------------------- //


// [UV, UW, VW] * W - fluxes
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::uvw_flux(
	T* flux,				// node: [W]

	const T* const UVW,		// node: [W]
	const T* const UW,		// node: [W]
	const T* const VW,		// node: [W]
	const T* const UV_uvw,	// node: [W]
	const T* const UW_uvw,	// node: [W]
	const T* const VW_uvw,	// node: [W]
	const T* const U,		// node: [C]
	const T* const V,		// node: [C]
	const T* const W,		// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _   __   _   __   _       _   _   _
	// u'v'w' = UVW - UW * V - VW * U - UV * W + 2 * U * V * W 
	//
	// [UW, VW, UV, W] averages have to be known at all [W] nodes, including walls
	// [U, V] averages have to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes

		flux[k] = UVW[k] -
			(T)0.25 * (UW[k] + UW_uvw[k]) * (V[k] + V[k - 1]) -
			(T)0.25 * (VW[k] + VW_uvw[k]) * (U[k] + U[k - 1]) -
			UV_uvw[k] * W[k] +

			(T)0.5 * W[k] * (U[k] + U[k - 1]) * (V[k] + V[k - 1]);
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::uww_flux(
	T* flux,						// node: [C]

	const T* const UWW,				// node: [C]
	const T* const W2_w,			// node: [W]
	const T* const W2_c,			// node: [C]
	const T* const W2_u,			// node: [C]
	const T* const W2_uw,			// node: [W]
	const T* const UW,				// node: [W]
	const T* const UW_bottom_uw,	// node: [W (C -- W)]
	const T* const UW_top_uw,		// node: [W (C -- W)]
	const T* const U,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _   _   __   __   _       _   _   _
	// u'w'w' = UWW - WW * U - W * UW - UW * W + 2 * U * W * W
	//
	// [W^2 (at -UW and -W nodes), W, UW] averages have to be known at all [W] nodes, including walls
	// [U] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
	// [UW-top] average has to be known at (gcz, nz - gcz - 1) nodes
	// [UW-bottom] average has to be known at (gcz + 1, nz - gcz) nodes
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {

		flux[k] = UWW[k] -

			// __   _
			// WW * U
			(T)0.5 * (

			(T)0.25 * (U[k + 1] + (T)2.0 * U[k] + U[k - 1]) *
			(T)0.25 * (W2_w[k] + (T)2.0 * W2_c[k] + W2_w[k + 1])
			+
			(T)0.5 *
			(T)0.25 * (W2_uw[k + 1] + W2_u[k]) * (U[k] + U[k + 1])
			+
			(T)0.5 *
			(T)0.25 * (W2_uw[k] + W2_u[k]) * (U[k] + U[k - 1])

			)

			-

			//   _   __   __   _
			// - W * UW - UW * W
			(
			(T)0.25 * (W[k] + W[k + 1]) * (UW_bottom_uw[k + 1] + UW_top_uw[k]) +
			(T)0.25 * (W[k + 1] * UW_bottom_uw[k + 1] + W[k] * UW_top_uw[k]) +
			(T)0.125 * (W[k] + W[k + 1]) * (UW[k] + UW[k + 1])
			)

			+

			//     _   _   _
			// 2 * W * W * U
			(T)0.5 * (W[k] + W[k + 1]) *
			(
			(T)0.125 * (W[k] + W[k + 1]) * (U[k + 1] + (T)2.0 * U[k] + U[k - 1]) +
			(T)0.25 * W[k + 1] * (U[k] + U[k + 1]) + (T)0.25 * W[k] * (U[k] + U[k - 1])
			);
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::vww_flux(
	T* flux,						// node: [C]

	const T* const VWW,				// node: [C]
	const T* const W2_w,			// node: [W]
	const T* const W2_c,			// node: [C]
	const T* const W2_v,			// node: [C]
	const T* const W2_vw,			// node: [W]
	const T* const VW,				// node: [W]
	const T* const VW_bottom_vw,	// node: [W (C -- W)]
	const T* const VW_top_vw,		// node: [W (C -- W)]
	const T* const V,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _   _   __   __   _       _   _   _
	// v'w'w' = VWW - WW * V - W * VW - VW * W + 2 * V * W * W
	//
	// [W^2 (at -VW and -W nodes), W, VW] averages have to be known at all [W] nodes, including walls
	// [V] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
	// [VW-top] average has to be known at (gcz, nz - gcz - 1) nodes
	// [VW-bottom] average has to be known at (gcz + 1, nz - gcz) nodes
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {

		flux[k] = VWW[k] -

			// __   _
			// WW * V
			(T)0.5 * (

			(T)0.25 * (V[k + 1] + (T)2.0 * V[k] + V[k - 1]) *
			(T)0.25 * (W2_w[k] + (T)2.0 * W2_c[k] + W2_w[k + 1])
			+
			(T)0.5 *
			(T)0.25 * (W2_vw[k + 1] + W2_v[k]) * (V[k] + V[k + 1])
			+
			(T)0.5 *
			(T)0.25 * (W2_vw[k] + W2_v[k]) * (V[k] + V[k - 1])

			)

			-

			//   _   __   __   _
			// - W * VW - VW * W
			(
			(T)0.25 * (W[k] + W[k + 1]) * (VW_bottom_vw[k + 1] + VW_top_vw[k]) +
			(T)0.25 * (W[k + 1] * VW_bottom_vw[k + 1] + W[k] * VW_top_vw[k]) +
			(T)0.125 * (W[k] + W[k + 1]) * (VW[k] + VW[k + 1])
			)

			+

			//     _   _   _
			// 2 * W * W * V
			(T)0.5 * (W[k] + W[k + 1]) *
			(
			(T)0.125 * (W[k] + W[k + 1]) * (V[k + 1] + (T)2.0 * V[k] + V[k - 1]) +
			(T)0.25 * W[k + 1] * (V[k] + V[k + 1]) + (T)0.25 * W[k] * (V[k] + V[k - 1])
			);
	}
}
// -------------------------------------------------------------------------------------------- //


// [CU, CV, CW] * W - fluxes
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cuw_flux(
	T* flux,						// node: [W]

	const T* const CUW,				// node: [W]
	const T* const UW,				// node: [W]
	const T* const CU_uw,			// node: [W]
	const T* const CW,				// node: [W]
	const T* const CW_uw,			// node: [W]
	const T* const C,				// node: [C]
	const T* const U,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _   _   __   _   __       _   _   _
	// c'u'w' = CUW - CU * W - C * UW - U * CW + 2 * C * U * W
	//
	// [CUW, UW, CW, CU, W] averages have to be known at all [W] nodes, including walls
	// [U, C] averages have to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {

		flux[k] = CUW[k] - (T)0.5 * UW[k] * (C[k] + C[k - 1])
			- (T)0.25 * (CW_uw[k] + CW[k]) * (U[k] + U[k - 1])
			- CU_uw[k] * W[k] +

			(T)2.0 * (T)0.25 * (U[k] + U[k - 1]) * (C[k] + C[k - 1]) * W[k];
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::cvw_flux(
	T* flux,						// node: [W]

	const T* const CVW,				// node: [W]
	const T* const VW,				// node: [W]
	const T* const CV_vw,			// node: [W]
	const T* const CW,				// node: [W]
	const T* const CW_vw,			// node: [W]
	const T* const C,				// node: [C]
	const T* const V,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _   _   __   _   __       _   _   _
	// c'v'w' = CVW - CV * W - C * VW - V * CW + 2 * C * V * W
	//
	// [CVW, VW, CW, CV, W] averages have to be known at all [W] nodes, including walls
	// [V, C] averages have to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {

		flux[k] = CVW[k] - (T)0.5 * VW[k] * (C[k] + C[k - 1])
			- (T)0.25 * (CW_vw[k] + CW[k]) * (V[k] + V[k - 1])
			- CV_vw[k] * W[k] +

			(T)2.0 * (T)0.25 * (V[k] + V[k - 1]) * (C[k] + C[k - 1]) * W[k];
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::cww_flux(
	T* flux,						// node: [C]

	const T* const CWW,				// node: [C]
	const T* const W2_w,			// node: [W]
	const T* const W2_c,			// node: [C]
	const T* const CW,				// node: [W]
	const T* const CW_bottom_w,		// node: [W (C -- W)]
	const T* const CW_top_w,		// node: [W (C -- W)]
	const T* const C,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	// ______   ___   __   _   _   __   __   _       _   _   _
	// c'w'w' = CWW - WW * C - W * CW - WC * W + 2 * C * W * W
	//
	// [W^2, CW, W] averages have to be known at all [W] nodes, including walls
	// [C] average has to be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
	// [CW-top] average has to be known at (gcz, nz - gcz - 1) nodes
	// [CW-bottom] average has to be known at (gcz + 1, nz - gcz) nodes
{
	int k;
#pragma omp parallel for private(k) shared(flux)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		flux[k] = CWW[k] 
			
			// __   _
			// WW * C
			- (T)0.5 * (
			(T)0.25 * (T)0.25 * (W2_w[k] + W2_w[k + 1] + (T)2.0 * W2_c[k]) * (C[k + 1] + (T)2.0 * C[k] + C[k - 1]) +
			(T)0.25 * (T)0.5 * (W2_w[k] + W2_c[k]) * (C[k] + C[k - 1]) +
			(T)0.25 * (T)0.5 * (W2_w[k + 1] + W2_c[k]) * (C[k] + C[k + 1]))

			// _   __
			// W * CW
			- (T)0.125 * (W[k] + W[k + 1]) * (CW_bottom_w[k + 1] + CW_top_w[k] + CW[k] + CW[k + 1])
			- (T)0.25 * (
			(T)0.5 * (W[k] + W[k + 1]) * (CW_bottom_w[k + 1] + CW_top_w[k]) + 
			W[k] * CW_top_w[k] + W[k + 1] * CW_bottom_w[k + 1])
			
			//     _   _   _
			// 2 * W * W * C
			+ (T)0.5 * (W[k] + W[k + 1]) * (
			(T)0.5 * (T)0.25 * (W[k] + W[k + 1]) * (C[k + 1] + (T)2.0 * C[k] + C[k - 1]) +
			(T)0.25 * W[k] * (C[k] + C[k - 1]) + (T)0.25 * W[k + 1] * (C[k] + C[k + 1]));
	}
}
// -------------------------------------------------------------------------------------------- //


// c - pressure gradient covariances
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::c_w_pressure_gradient_turb(
	T* c_dpdz_turb,				// node: [W]

	const T* const C_dPdz,		// node: [W]
	const T* const C,			// node: [C]
	const T* const Pressure,	// node: [C]
	const wstGrid3d< T >& grid)
	// C-pressure gradient covariance
	// __________
	//       dp'
	// c' * ----
	//       dz
	//
	// [C, P] averages have be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(c_dpdz_turb)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		c_dpdz_turb[k] = C_dPdz[k] -
			(C[k] + C[k - 1]) * (Pressure[k] - Pressure[k - 1]) * grid.dzmi[k];
	}
}
// -------------------------------------------------------------------------------------------- //