bl-turb-scalar.h

#pragma once

// [bl-turb-scalar.h]: boundary-layer scalar turbulence statistics and budgets
//
// -------------------------------------------------------------------------------------------- //

#include "wstgrid3d.h"

//
// *[Note]: we need gcz >= 2 for computation of production terms:
//         _
//        dC
// - w'w' --
//        dz
//


namespace nse
{
	// Heat eq. balance
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void heat_eq(
		T* heat_balance,				// node: [W]
		T* turbulent_heat_flux,			// node: [W]
		T* heat_stress,					// node: [W]

		const T* const Tw_flux,			// node: [W]
		const T* const T_grad,			// node: [W]
		const T c_diffusivity, const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// SVA production
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void SVA_production(T* _SVA_production,		// node: [C]

		const T* const CW_bottom,				// node: [C (W -- C)]
		const T* const CW_top,					// node: [C (W -- C)]
		const T* const C,						// node: [C]
		const T* const W,						// node: [W]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// SVA transport
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void SVA_transport(T* _SVA_transport,	// node: [C]

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

	// SVA dissipation
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void SVA_dissipation(T* _SVA_dissipation,	// node: [C]

		const T* const C_dissipation,			// node: [C]
		const T* const C,						// node: [C]
		const T c_diffusivity, const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// SVA iso dissipation
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void SVA_iso_dissipation(T* _SVA_iso_dissipation,	// node: [C]

		const T* const C_iso_dissipation,				// node: [C]
		const T* const C,								// node: [C]
		const T c_diffusivity, const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //


	// c'ui' flux budget: production
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void cu_production_shear(T* _cu_production_shear,	// node: [C]

		const T* const CW_bottom_u,		// node: [C (W -- C)]
		const T* const CW_top_u,		// node: [C (W -- C)]
		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 cu_production_gradC(T* _cu_production_gradC,	// node: [C]

		const T* const UW_bottom,		// node: [C (W -- C)]
		const T* const UW_top,			// node: [C (W -- C)]
		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 cv_production_shear(T* _cv_production_shear,	// node: [C]

		const T* const CW_bottom_v,		// node: [C (W -- C)]
		const T* const CW_top_v,		// node: [C (W -- C)]
		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 cv_production_gradC(T* _cv_production_gradC,	// node: [C]

		const T* const VW_bottom,		// node: [C (W -- C)]
		const T* const VW_top,			// node: [C (W -- C)]
		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 cw_production_shear(T* _cw_production_shear,	// 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: [W]
		const T* const W,				// node: [C]
		const wstGrid3d< T >& grid);

	template< typename T >
	void cw_production_gradC(T* _cw_production_gradC,	// node: [W]

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

	// c'ui' flux budget: transport
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void cu_transport(T* _cu_transport,		// node: [C]

		const T* const cuw_flux,			// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void cv_transport(T* _cv_transport,		// node: [C]

		const T* const cvw_flux,			// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void cw_transport(T* _cw_transport,		// node: [W]

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

	// c'ui' flux budget: pressure scrambles
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void cw_pressure_work(T* _cw_pressure_work,		// node: [W]

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


	template< typename T >
	void cu_pressure_gradc(T* _cu_pressure_gradc,	// node: [C]
		// computing as: 
		//   p'*dc'/dx = { d(c'p')/dx = 0 } - c'*dp'/dx
		//											

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

	template< typename T >
	void cv_pressure_gradc(T* _cv_pressure_gradc,	// node: [C]
		// computing as: 
		//   p'*dc'/dy = { d(c'p')/dy = 0 } - c'*dp'/dy
		//											

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

	template< typename T >
	void cw_pressure_gradc(T* _cw_pressure_gradc,	// node: [W]
		// computing as: 
		//   p'*dc'/dz = d(c'p')/dz - c'*dp'/dz
		//

		const T* const _cw_pressure_work,			// node: [W]
		const T* const c_dpdz_turb,					// node: [W]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// c'ui' flux budget: dissipation
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void cu_dissipation(T* _cu_dissipation,			// node: [C]

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

	template< typename T >
	void cv_dissipation(T* _cv_dissipation,			// node: [C]

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

	template< typename T >
	void cw_dissipation(T* _cw_dissipation,			// node: [W]

		const T* const CW_dissipation,				// node: [W]
		const T* const C,							// node: [C]
		const T* const W,							// node: [W]
		const T c_diffusivity, const T c_kinematic_viscosity,
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// c'w' flux budget: buoyancy 
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void cw_buoyancy(T* _cw_buoyancy,		// node: [W]

		const T* const C2_c,				// node: [C]
		const T* const C2_w,				// node: [W]
		const T* const C,					// node: [C]
		const T c_Richardson, const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //


	// TPE structure
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void TPE_structure(T* TPE,							// node: [C]
		T* TPE_heat_flux, T* TPE_diffusion,				// node: [C]
		T* TPE_dissipation, T* TPE_iso_dissipation,		// node: [C]

		const T* const TVA_production,			// node: [C]
		const T* const TVA_diffusion,			// node: [C]
		const T* const TVA_dissipation,			// node: [C]
		const T* const TVA_iso_dissipation,		// node: [C]

		const T* const T2,						// node: [C]
		const T* const Tc,						// node: [C]
		const T* const T_grad,					// node: [W]
		const T c_Richardson, const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //

	// Energy structure
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void energy_structure(T* TKE_share, T* TPE_share,	// node: [C]

		const T* const TKE,		// node: [C]
		const T* const TPE,		// node: [C]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //
}


// Heat eq. balance
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::heat_eq(
	T* heat_balance,				// node: [W]
	T* turbulent_heat_flux,			// node: [W]
	T* heat_stress,					// node: [W]

	const T* const Tw_flux,			// node: [W]
	const T* const T_grad,			// node: [W]
	const T c_diffusivity, const wstGrid3d< T >& grid)
	// integrated [T] averaged heat equation
	//                  _
	//    ____   1  1  dT
	//	- T'w' + -- -- -- = const 
	//			 Re Pr dz
	//    (1)      (2)
	// (1) - turbulent heat flux
	// (2) - heat stress
	// terms are defined at [W] nodes (averaged in [W] nodes)
	//
	// [T'w', dT/dz] have to be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(heat_balance, \
	turbulent_heat_flux, heat_stress)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) { // all [W] nodes
		turbulent_heat_flux[k] = -Tw_flux[k];
		heat_stress[k] = c_diffusivity * T_grad[k];

		heat_balance[k] = turbulent_heat_flux[k] + heat_stress[k];
	}
}
// -------------------------------------------------------------------------------------------- //

// SVA production
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::SVA_production(
	T* _SVA_production,			// node: [C]

	const T* const CW_bottom,	// node: [C (W -- C)]
	const T* const CW_top,		// node: [C (W -- C)]
	const T* const C,			// node: [C]
	const T* const W,			// node: [W]
	const wstGrid3d< T >& grid)
	// production term of scalar variance equation defined at [C] node
	//            _
	//    ____   dC
	//	- c'w' * -- 
	//			 dz
	// [W] average has 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(_SVA_production)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_SVA_production[k] = -(T) 0.5 * (
			// discretization is based on ADV. form //
			// correct for SKEW.form (checked) //
			(CW_bottom[k] - C[k] * W[k]) * (C[k] - C[k - 1]) * grid.dzi[k] +
			(CW_top[k] - C[k] * W[k + 1]) * (C[k + 1] - C[k]) * grid.dzi[k]);
	}
}
// -------------------------------------------------------------------------------------------- //

// SVA transport
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::SVA_transport(
	T* _SVA_transport,			// node: [C]

	const T* const cc_w_flux,	// node: [W]
	const wstGrid3d< T >& grid)
	// transport term of scalar variance equation defined at [C] node
	//       ______________
	//      d[(c')^2 * w')] 
	//	-  ---------------- 
	//			 dz
	// [c'c'w'] flux has to be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(_SVA_transport)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_SVA_transport[k] = -(T) 0.5 *
			(cc_w_flux[k + 1] - cc_w_flux[k]) * grid.dzi[k];
	}
}
// -------------------------------------------------------------------------------------------- //

// SVA dissipation
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::SVA_dissipation(
	T* _SVA_dissipation,			// node: [C]

	const T* const C_dissipation,	// node: [C]
	const T* const C,				// node: [C]
	const T c_diffusivity, const wstGrid3d< T >& grid)
	// dissipation component of scalar variance equation defined at [C] node
	//       ___________
	//          d^2(c')
	//	 k * c' ------- 
	//			dx(j)^2
	// [C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_SVA_dissipation)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_SVA_dissipation[k] = C_dissipation[k] -
			c_diffusivity * (C[k] *
			((C[k + 1] - C[k]) * grid.dzp2i[k] - (C[k] - C[k - 1]) * grid.dzm2i[k]));
	}
}
// -------------------------------------------------------------------------------------------- //

// SVA iso dissipation
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::SVA_iso_dissipation(
	T* _SVA_iso_dissipation,			// node: [C]

	const T* const C_iso_dissipation,	// node: [C]
	const T* const C,					// node: [C]
	const T c_diffusivity, const wstGrid3d< T >& grid)
	// isotropic dissipation component of scalar variance equation defined at [C] node
	//       _____________
	//       d(c')   d(c')
	// - k * ----- * -----  
	//		 dx(j)   dx(j)
	// [C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_SVA_iso_dissipation)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_SVA_iso_dissipation[k] = -(C_iso_dissipation[k] -
			c_diffusivity * (
			(T)0.5 * (C[k + 1] - C[k]) * (C[k + 1] - C[k]) * grid.dzp2i[k] +
			(T)0.5 * (C[k] - C[k - 1]) * (C[k] - C[k - 1]) * grid.dzm2i[k]));
	}
}
// -------------------------------------------------------------------------------------------- //

// c'ui' flux budget: production
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cu_production_shear(
	T* _cu_production_shear,		// node: [C]

	const T* const CW_bottom_u,		// node: [C (W -- C)]
	const T* const CW_top_u,		// node: [C (W -- C)]
	const T* const C,				// node: [C]
	const T* const U,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	// production [shear] term of [c'u'] budget equation defined at [C] node
	//           _
	//   ____   dU
	// - c'w' * --
	//          dz
	//
	// [W] average has to be known at all [W] nodes, including walls
	// [U] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_cu_production_shear)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cu_production_shear[k] = -(T)0.5 * (
			// discretization is based on ADV. form //
			(CW_bottom_u[k] - C[k] * W[k]) * (U[k] - U[k - 1]) * grid.dzi[k] +
			(CW_top_u[k] - C[k] * W[k + 1]) * (U[k + 1] - U[k]) * grid.dzi[k]);
	}
}

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

	const T* const UW_bottom,		// node: [C (W -- C)]
	const T* const UW_top,			// node: [C (W -- C)]
	const T* const C,				// node: [C]
	const T* const U,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	// production [gradC] term of [c'u'] budget equation defined at [C] node
	//           _
	//   ____   dC
	// - u'w' * --
	//          dz
	//
	// [W] average has to be known at all [W] nodes, including walls
	// [C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_cu_production_gradC)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cu_production_gradC[k] = -(T)0.5 * (
			// discretization is based on ADV. form //
			(UW_bottom[k] - U[k] * W[k]) * (C[k] - C[k - 1]) * grid.dzi[k] +
			(UW_top[k] - U[k] * W[k + 1]) * (C[k + 1] - C[k]) * grid.dzi[k]);
	}
}


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

	const T* const CW_bottom_v,		// node: [C (W -- C)]
	const T* const CW_top_v,		// node: [C (W -- C)]
	const T* const C,				// node: [C]
	const T* const V,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	// production [shear] term of [c'v'] budget equation defined at [C] node
	//           _
	//   ____   dV
	// - c'w' * --
	//          dz
	//
	// [W] average has to be known at all [W] nodes, including walls
	// [V] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_cv_production_shear)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cv_production_shear[k] = -(T)0.5 * (
			// discretization is based on ADV. form //
			(CW_bottom_v[k] - C[k] * W[k]) * (V[k] - V[k - 1]) * grid.dzi[k] +
			(CW_top_v[k] - C[k] * W[k + 1]) * (V[k + 1] - V[k]) * grid.dzi[k]);
	}
}

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

	const T* const VW_bottom,		// node: [C (W -- C)]
	const T* const VW_top,			// node: [C (W -- C)]
	const T* const C,				// node: [C]
	const T* const V,				// node: [C]
	const T* const W,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	// production [gradC] term of [c'v'] budget equation defined at [C] node
	//           _
	//   ____   dC
	// - v'w' * --
	//          dz
	//
	// [W] average has to be known at all [W] nodes, including walls
	// [C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_cv_production_gradC)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cv_production_gradC[k] = -(T)0.5 * (
			// discretization is based on ADV. form //
			(VW_bottom[k] - V[k] * W[k]) * (C[k] - C[k - 1]) * grid.dzi[k] +
			(VW_top[k] - V[k] * W[k + 1]) * (C[k + 1] - C[k]) * grid.dzi[k]);
	}
}


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

	const T* const CW_bottom_w,		// node: [W (C -- W)]
	const T* const CW_top_w,		// node: [W (C -- W)]
	const T* C,						// node: [C]
	const T* W,						// node: [W]
	const wstGrid3d< T >& grid)
	//
	// production [shear] term of [c'w'] budget equation defined at [W] node
	//           _
	//   ____   dW
	// - c'w' * --
	//          dz
	//
	// [W] average has to be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(_cw_production_shear)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {	// all [W] nodes, excluding walls
		_cw_production_shear[k] = - (
			// discretization is based on ADV. form //
			(CW_bottom_w[k] - (T)0.25 * (C[k] + C[k - 1]) * (W[k] + W[k - 1])) * (W[k] - W[k - 1]) * grid.dzmi[k] +
			(CW_top_w[k] - (T)0.25 * (C[k] + C[k - 1]) * (W[k] + W[k + 1])) * (W[k + 1] - W[k]) * grid.dzmi[k]);
	}

	w_dirichlet_bc_z(_cw_production_shear, (T)0, (T)0, grid);
}

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

	const T* const W2_w,		// node: [W]
	const T* const W2_c,		// node: [C]
	const T* C,					// node: [C]
	const T* W,					// node: [W]
	const wstGrid3d< T >& grid)
	//
	// production [gradC] term of [c'w'] budget equation defined at [W] node
	//           _
	//   ____   dC
	// - w'w' * --
	//          dz
	//
	// [W^2, W] averages have to be known at all [W] nodes, including walls
	// [C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_cw_production_gradC)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {	// all [W] nodes, excluding walls
		_cw_production_gradC[k] = - (
			// discretization is based on ADV. form //
			(W2_c[k] - W[k] * W[k + 1]) * (C[k + 1] - C[k]) * grid.dziq[k] +
			(W2_c[k - 1] - W[k] * W[k - 1]) * (C[k - 1] - C[k - 2]) * grid.dziq[k - 1] +

			(W2_w[k] - W[k] * W[k]) * (C[k] - C[k - 1]) * (grid.dziq[k] + grid.dziq[k - 1]));
	}

	w_dirichlet_bc_z(_cw_production_gradC, (T)0, (T)0, grid);
}
// -------------------------------------------------------------------------------------------- //

// c'ui' flux budget: transport
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cu_transport(T* _cu_transport,	// node: [C]

	const T* const cuw_flux,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	// transport term of [c'u'] equation defined at [C] node
	//       ______________
	//      d[(c'u') * u'] 
	//	-  ---------------- 
	//			 dz
	// [c'u'w'] flux has be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(_cu_transport)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cu_transport[k] = -(cuw_flux[k + 1] - cuw_flux[k]) * grid.dzi[k];
	}
}

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

	const T* const cvw_flux,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	// transport term of [c'v'] equation defined at [C] node
	//       ______________
	//      d[(c'v') * v'] 
	//	-  ---------------- 
	//			 dz
	// [c'v'w'] flux has be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(_cv_transport)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cv_transport[k] = -(cvw_flux[k + 1] - cvw_flux[k]) * grid.dzi[k];
	}
}

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

	const T* const cww_flux,				// node: [C]
	const wstGrid3d< T >& grid)
	//
	// transport term of [c'w'] budget equation defined at [W] node
	//       ______________
	//      d[(c'w') * w'] 
	//	-  ---------------- 
	//			 dz
	// [c'w'w'] flux has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
#pragma omp parallel for private(k) shared(_cw_transport)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cw_transport[k] = -(cww_flux[k] - cww_flux[k - 1]) * (T)2.0 * grid.dzmi[k];
	}
}
// -------------------------------------------------------------------------------------------- //

// c'ui' flux budget: pressure scrambles
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cw_pressure_work(T* _cw_pressure_work,	// node: [W]

	const T* const cp_flux,							// node: [C]
	const wstGrid3d< T >& grid)
	//
	// pressure work term of [c'w'] budget equation defined at [W]
	//       ______
	//      d[c'p'] 
	//	-  -------- 
	//		  dz
	// [c'p'] flux 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(_cw_pressure_work)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		_cw_pressure_work[k] = -(cp_flux[k] - cp_flux[k - 1]) * (T)2.0 * grid.dzmi[k];
	}
}


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

	const T* const c_dpdx_turb,		// node: [C]
	const wstGrid3d< T >& grid)
	//
	//  p'*dc'/dx = { d(c'p')/dx = 0 } - c'*dp'/dx
	//
{
	int k;
#pragma omp parallel for private(k) shared(_cu_pressure_gradc)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cu_pressure_gradc[k] = - c_dpdx_turb[k];
	}
}

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

	const T* const c_dpdy_turb,			// node: [C]
	const wstGrid3d< T >& grid)
	//
	//  p'*dc'/dy = { d(c'p')/dy = 0 } - c'*dp'/dy
	//
{
	int k;
#pragma omp parallel for private(k) shared(_cv_pressure_gradc)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_cv_pressure_gradc[k] = -c_dpdy_turb[k];
	}
}

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

	const T* const _cw_pressure_work,		// node: [W]
	const T* const c_dpdz_turb,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	//  p'*dc'/dz = d(c'p')/dz - c'*dp'/dz
	//
{
	int k;
#pragma omp parallel for private(k) shared(_cw_pressure_gradc)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		_cw_pressure_gradc[k] = -_cw_pressure_work[k] - c_dpdz_turb[k];
	}
}
// -------------------------------------------------------------------------------------------- //

// c'ui' flux budget: dissipation
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cu_dissipation(
	T* _cu_dissipation,					// node: [C]

	const T* const CU_dissipation,		// node: [C]
	const T* const C,					// node: [C]
	const T* const U,					// node: [C]
	const T c_diffusivity, const T c_kinematic_viscosity,
	const wstGrid3d< T >& grid)
	//
	//       ___________       ___________
	//          d^2(c')           d^2(u')
	//	 k * u' ------- + nu * c' ------- 
	//			dx(j)^2           dx(j)^2
	//
	// [U, C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
	T C_diffusion, U_diffusion;

#pragma omp parallel for private(k, C_diffusion, U_diffusion) shared(_cu_dissipation)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		
		C_diffusion = c_diffusivity *
			((C[k + 1] - C[k]) * grid.dzp2i[k] - (C[k] - C[k - 1]) * grid.dzm2i[k]);
		U_diffusion = c_kinematic_viscosity *
			((U[k + 1] - U[k]) * grid.dzp2i[k] - (U[k] - U[k - 1]) * grid.dzm2i[k]);

		_cu_dissipation[k] = CU_dissipation[k] -
			(C[k] * U_diffusion + C_diffusion * U[k]);
	}
}

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

	const T* const CV_dissipation,		// node: [C]
	const T* const C,					// node: [C]
	const T* const V,					// node: [C]
	const T c_diffusivity, const T c_kinematic_viscosity,
	const wstGrid3d< T >& grid)
	//
	//       ___________       ___________
	//          d^2(c')           d^2(v')
	//	 k * v' ------- + nu * c' ------- 
	//			dx(j)^2           dx(j)^2
	//
	// [V, C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	int k;
	T C_diffusion, V_diffusion;

#pragma omp parallel for private(k, C_diffusion, V_diffusion) shared(_cv_dissipation)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {

		C_diffusion = c_diffusivity *
			((C[k + 1] - C[k]) * grid.dzp2i[k] - (C[k] - C[k - 1]) * grid.dzm2i[k]);
		V_diffusion = c_kinematic_viscosity *
			((V[k + 1] - V[k]) * grid.dzp2i[k] - (V[k] - V[k - 1]) * grid.dzm2i[k]);

		_cv_dissipation[k] = CV_dissipation[k] -
			(C[k] * V_diffusion + C_diffusion * V[k]);
	}
}

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

	const T* const CW_dissipation,		// node: [W]
	const T* const C,					// node: [C]
	const T* const W,					// node: [W]
	const T c_diffusivity, const T c_kinematic_viscosity,
	const wstGrid3d< T >& grid)
	//
	//       ___________       ___________
	//          d^2(c')           d^2(w')
	//	 k * w' ------- + nu * c' ------- 
	//			dx(j)^2           dx(j)^2
	//
	// [CW-dissipation, W] averages have to be known at all [W] nodes, including walls
	// [C] average has be known at all [C] nodes and ghost nodes: (k + 1/2), (k - 1/2)
{
	T *C_diffusion, *W_diffusion;
	int c_buf_id = memStx::get_buf(&C_diffusion, grid.nz);
	int w_buf_id = memStx::get_buf(&W_diffusion, grid.nz);

	null(C_diffusion, grid.nz);
	null(W_diffusion, grid.nz);

	int k;

#pragma omp parallel for private(k) shared(C_diffusion, W_diffusion)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		C_diffusion[k] = c_diffusivity *
			((C[k + 1] - C[k]) * grid.dzp2i[k] - (C[k] - C[k - 1]) * grid.dzm2i[k]);
		W_diffusion[k] = c_kinematic_viscosity *
			((W[k + 1] - W[k]) * grid.dzm2i[k] - (W[k] - W[k - 1]) * grid.dzp2i[k - 1]);
	}

	// assuming W = 0 & ~linear[laplace(W)] = 0 at walls -> no need for laplace(C) b.c.
	//
	w_dirichlet_bc_z(W_diffusion, (T)0, (T)0, grid);

#pragma omp parallel for private(k) shared(_cw_dissipation,\
	C_diffusion, W_diffusion)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		_cw_dissipation[k] = CW_dissipation[k] -
			(T)0.5 * (
			(C[k] + C[k - 1]) * W_diffusion[k] +
			W[k] * (C_diffusion[k] + C_diffusion[k - 1]));
	}

	memStx::free_buf(c_buf_id);
	memStx::free_buf(w_buf_id);
}
// -------------------------------------------------------------------------------------------- //

// c'w' flux budget: buoyancy
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::cw_buoyancy(
	T* _cw_buoyancy,			// node: [W]

	const T* const C2_c,		// node: [C]
	const T* const C2_w,		// node: [W]
	const T* const C,			// node: [C]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	//      ____
	// Ri * c'c'
	//
	// [C^2 at W-node] average has to be known at all [W] nodes, including walls
	// [C, C^2 at C-node] 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(_cw_buoyancy)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		_cw_buoyancy[k] = (T)0.25 * c_Richardson * (
			C2_c[k] - C[k] * C[k] +
			C2_c[k - 1] - C[k - 1] * C[k - 1] + 
			(T)2.0 * (C2_w[k] - C[k] * C[k - 1])
			);
	}
}
// -------------------------------------------------------------------------------------------- //


// TPE structure
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::TPE_structure(
	T* TPE,											// node: [C]
	T* TPE_heat_flux, T* TPE_diffusion,				// node: [C]
	T* TPE_dissipation, T* TPE_iso_dissipation,		// node: [C]

	const T* const TVA_production,		// node: [C]
	const T* const TVA_diffusion,		// node: [C]
	const T* const TVA_dissipation,		// node: [C]
	const T* const TVA_iso_dissipation,	// node: [C]

	const T* const T2,					// node: [C]
	const T* const Tc,					// node: [C]
	const T* const T_grad,				// node: [W]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	// TPE = 1/2*T'T'*Ri(b)/[dT/dz]
	// TPE balance: [temperature variance] * (Ri(b)/[dT/dz])
	//
	// [dT/dz] has to be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(TPE,\
	TPE_heat_flux, TPE_diffusion, TPE_dissipation, TPE_iso_dissipation)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		TPE[k] = (T)0.5 * (T2[k] - Tc[k] * Tc[k]) *
			(c_Richardson / ((T)0.5 * (T_grad[k] + T_grad[k + 1])));

		TPE_heat_flux[k] = TVA_production[k] *
			(c_Richardson / ((T)0.5 * (T_grad[k] + T_grad[k + 1])));
		TPE_diffusion[k] = TVA_diffusion[k] *
			(c_Richardson / ((T)0.5 * (T_grad[k] + T_grad[k + 1])));
		TPE_dissipation[k] = TVA_dissipation[k] *
			(c_Richardson / ((T)0.5 * (T_grad[k] + T_grad[k + 1])));
		TPE_iso_dissipation[k] = TVA_iso_dissipation[k] *
			(c_Richardson / ((T)0.5 * (T_grad[k] + T_grad[k + 1])));
	}
}
// -------------------------------------------------------------------------------------------- //

// Energy structure
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::energy_structure(
	T* TKE_share, T* TPE_share,		// node: [C]

	const T* const TKE,				// node: [C]
	const T* const TPE,				// node: [C]
	const wstGrid3d< T >& grid)
	// 
	// TKE-TPE shares
	//
{
	int k;
#pragma omp parallel for private(k) shared(TKE_share, TPE_share)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		TKE_share[k] = TKE[k] / (TKE[k] + TPE[k]);
		TPE_share[k] = TPE[k] / (TKE[k] + TPE[k]);
	}
}
// -------------------------------------------------------------------------------------------- //