#pragma once

// [bl-scale-def.h]: boundary-layer scales and dimensionless values definitions
//
// -------------------------------------------------------------------------------------------- //

#include "wstgrid3d.h"


namespace nse
{
	// Time scales
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void time_scale_turbulent(T* _time_scale_turbulent,		// node: [C]

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

	template< typename T >
	void time_scale_svariance(T* _time_scale_svariance,		// node: [C]

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

	// Length scales
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void length_scale_kolmogorov(T* _length_scale_kolmogorov,	// node: [C]

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

	template< typename T >
	void length_scale_mixing(T* _length_scale_mixing,			// node: [W]

		const T* const U2,						// node: [W]
		const T* const V2,						// node: [W]
		const T* const W2,						// node: [W]
		const T* const U,						// node: [C]
		const T* const V,						// node: [C]
		const T* const W,						// node: [W]
		const T* const U_grad,					// node: [W]
		const wstGrid3d< T >& grid);

	template< typename T >
	void length_scale_ellison(T* _length_scale_ellison,			// node: [W]

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

	template< typename T >
	void length_scale_ozmidov(T* _length_scale_ozmidov,			// node: [C]

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

	template< typename T >
	void length_scale_obukhov(T* _length_scale_obukhov,			// node: [W]

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

	// Dimensionless numbers
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void prandtl_turbulent(T* Prandtl_turbulent,		// node: [W]

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

	template< typename T >
	void richardson_gradient(T* Richardson_gradient,	// node: [W]

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

	template< typename T >
	void richardson_flux(T* Richardson_flux,			// node: [W]

		const T* const uw_flux,		// node: [W]
		const T* const U_grad,		// node: [W]
		const T* const Tw_flux,		// node: [W]
		const T c_Richardson, const wstGrid3d< T >& grid);

	template< typename T >
	void reynolds_buoyancy(T* Reynolds_buoyancy,		// node: [C]

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

	template< typename T >
	void froude_horizontal(T* Froude_horizontal,		// node: [C]

		const T* const u_TKE,							// node: [C]
		const T* const v_TKE,							// node: [C]
		const T* const TKE_iso_dissipation,				// node: [C]
		const T* const T_grad,							// node: [W]
		const T c_Richardson, const wstGrid3d< T >& grid);

	template< typename T >
	void mixing_efficiency(T* _mixing_efficiency,		// node: [C]

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

	template< typename T >
	void turbulence_production_ratio(T* turb_production_ratio,	// node: [C]

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


	// pressure-strain models
	// -------------------------------------------------------------------------------------------- //
	template< typename T >
	void Rotta_model(T* u_Rotta, T* v_Rotta, T* w_Rotta,	// node: [C]
		T* uw_Rotta,										// node: [C]

		const T* const TKE, const T* const TKE_iso_dissipation,					// node: [C]
		const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]
		const T* const uw_flux,													// node: [W]

		const T* const u_TKE_exchange,	// node: [C]
		const T* const v_TKE_exchange,	// node: [C] 
		const T* const w_TKE_exchange,	// node: [C]
		const T* const P2Suw_turb_c,	// node: [�]
		const wstGrid3d< T >& grid);

	template< typename T >
	void RDT_model(T* u_RDT, T* v_RDT, T* w_RDT,			// node: [C]
		T* uw_RDT,											// node: [C]

		const T* const TKE_production,		// node: [C]
		const T* const W2_w,				// node: [W]
		const T* const U,					// node: [C]
		const T* const W,					// node: [W]

		const T* const u_TKE_exchange,	// node: [C]
		const T* const v_TKE_exchange,	// node: [C] 
		const T* const w_TKE_exchange,	// node: [C]
		const T* const P2Suw_turb_c,	// node: [�]
		const wstGrid3d< T >& grid);

	template< typename T >
	void Rotta_RDT_model(T* Rotta_RDT_e, T* Rotta_RDT_p,	// node: [C]

		const T* const TKE,					// node: [C]
		const T* const TKE_iso_dissipation,	// node: [C]
		const T* const TKE_production,		// node: [C]
		const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]
		
		const T* const u_TKE_exchange,	// node: [C]
		const T* const v_TKE_exchange,	// node: [C] 
		const T* const w_TKE_exchange,	// node: [C]
		const wstGrid3d< T >& grid);

	template< typename T >
	void Rotta_buoyancy_model(T* Rotta_buoyancy_e, T* Rotta_buoyancy_b,		// node: [C]

		const T* const TKE,					// node: [C]
		const T* const TKE_iso_dissipation,	// node: [C]
		const T* const TKE_heat_flux,		// node: [C]
		const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]

		const T* const u_TKE_exchange,	// node: [C]
		const T* const v_TKE_exchange,	// node: [C] 
		const T* const w_TKE_exchange,	// node: [C]
		const T c_Richardson, const wstGrid3d< T >& grid);

	template< typename T >
	void RDT_buoyancy_model(T* RDT_buoyancy_p, T* RDT_buoyancy_b,			// node: [C]

		const T* const TKE_production,		// node: [C]
		const T* const TKE_heat_flux,		// node: [C]

		const T* const u_TKE_exchange,	// node: [C]
		const T* const v_TKE_exchange,	// node: [C] 
		const T* const w_TKE_exchange,	// node: [C]
		const T c_Richardson, const wstGrid3d< T >& grid);

	template< typename T >
	void Rotta_TPE_model(T* u_Rotta_TPE, T* v_Rotta_TPE, T* w_Rotta_TPE,	// node: [C]

		const T* const TKE,					// node: [C]
		const T* const TPE,					// node: [C]
		const T* const TKE_iso_dissipation,	// node: [C]
		const T* const TPE_iso_dissipation,	// node: [C]
		const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]

		const T* const u_TKE_exchange,	// node: [C]
		const T* const v_TKE_exchange,	// node: [C] 
		const T* const w_TKE_exchange,	// node: [C]
		const wstGrid3d< T >& grid);
	// -------------------------------------------------------------------------------------------- //
}

// Time scales
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::time_scale_turbulent(
	T* _time_scale_turbulent,				// node: [C]

	const T* const TKE,						// node: [C]
	const T* const TKE_iso_dissipation,		// node: [C]
	const wstGrid3d< T >& grid)
	//
	// T = (1/2) * [(u')^2 + (v'2)^2 + (w')^2] / ([dui'/dxj]*[dui'/dxj])
	//
{
	int k;
#pragma omp parallel for private(k) shared(_time_scale_turbulent)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_time_scale_turbulent[k] = TKE[k] / (-TKE_iso_dissipation[k]);
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const C2,						// node: [C]
	const T* const C,						// node: [C]
	const T* const CVA_iso_dissipation,		// node: [C]
	const wstGrid3d< T >& grid)
	//
	// T = (1/2)*[c'^2] / ([dc'/dxj] * [dc'/dxj])
	//
{
	int k;
#pragma omp parallel for private(k) shared(_time_scale_svariance)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_time_scale_svariance[k] = ((T)0.5*(
			C2[k] - C[k] * C[k])) / (-CVA_iso_dissipation[k]);
	}
}
// -------------------------------------------------------------------------------------------- //

// Length scales
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::length_scale_kolmogorov(
	T* _length_scale_kolmogorov,			// node: [C]

	const T* const TKE_iso_dissipation,		// node: [C]
	const T c_kinematic_viscosity, const wstGrid3d< T >& grid)
	//
	// L = [ nu^3 / e(ke) ]^(1/4), e - isotropic dissipation
	//
{
	int k;
#pragma omp parallel for private(k) shared(_length_scale_kolmogorov)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_length_scale_kolmogorov[k] = pow(
			pow(c_kinematic_viscosity, (T)3.0) / (-TKE_iso_dissipation[k]),
			(T)0.25);
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const U2,					// node: [W]
	const T* const V2,					// node: [W]
	const T* const W2,					// node: [W]
	const T* const U,					// node: [C]
	const T* const V,					// node: [C]
	const T* const W,					// node: [W]
	const T* const U_grad,				// node: [W]
	const wstGrid3d< T >& grid)
	//
	// L = sqrt[(u')^2 + (v'2)^2 + (w')^2] / [dU/dz]
	//
	// [U^2, V^2, W^2, W, dU/dz] 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(_length_scale_mixing)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		_length_scale_mixing[k] = sqrt(
			(U2[k] - U[k] * U[k - 1]) +
			(V2[k] - V[k] * V[k - 1]) +
			(W2[k] - W[k] * W[k])
		) / U_grad[k];
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const T2,			// node: [W]
	const T* const Tc,			// node: [C]
	const T* const T_grad,		// node: [W]
	const wstGrid3d< T >& grid)
	//
	// L = sqrt[T'^2] / [dT/dz]
	//
	// [T^2, dT/dz] averages have to be known at all [W] nodes, including walls
	// [T] 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(_length_scale_ellison)
	for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) {	// all [W] nodes
		_length_scale_ellison[k] = sqrt(
			T2[k] - Tc[k] * Tc[k - 1]) / T_grad[k];
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const TKE_iso_dissipation,		// node: [C]
	const T* const T_grad,					// node: [W]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	// L = ( e(TKE) / [Ri(b)*dT/dz]^3/2 )^1/2
	//
	// [dT/dz] average has to be known at all [W] nodes, including walls
{
	int k;
	if (fabs(c_Richardson) > 0) {

#pragma omp parallel for private(k) shared(_length_scale_ozmidov)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			_length_scale_ozmidov[k] = sqrt(
				-TKE_iso_dissipation[k] /
				((T)0.5 * c_Richardson * (T_grad[k] + T_grad[k + 1]) *
				sqrt((T)0.5 * c_Richardson * (T_grad[k] + T_grad[k + 1]))));
		}
	}
	else
	{

#pragma omp parallel for private(k) shared(_length_scale_ozmidov)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			_length_scale_ozmidov[k] = (T) 0.0;
		}
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const uw_flux,		// node: [W]
	const T* const Tw_flux,		// node: [W]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	// L = (u'w' * sqrt(|u'w'|)) / (Ri(b) * T'w')
	//
	// [u'w', T'w'] have to be known at all [W] nodes, excluding walls
{
	// define -k shifts to prevent division by zero //
	const int kbsh = (grid.mpi_com.rank_z == 0) ? 1 : 0;
	const int ktsh = (grid.mpi_com.rank_z == grid.mpi_com.size_z - 1) ? 1 : 0;

	int k;
	if (fabs(c_Richardson) > 0) {
#pragma omp parallel for private(k) shared(_length_scale_obukhov)
		for (k = grid.gcz + kbsh; k <= grid.nz - grid.gcz - ktsh; k++) {	// all [W] nodes, excluding walls
			_length_scale_obukhov[k] = (uw_flux[k] * sqrt(fabs(uw_flux[k]))) / (c_Richardson * Tw_flux[k]);
		}
	}
	else
	{
#pragma omp parallel for private(k) shared(_length_scale_obukhov)
		for (k = grid.gcz + kbsh; k <= grid.nz - grid.gcz - ktsh; k++) {	// all [W] nodes, excluding walls
			_length_scale_obukhov[k] = (T)0;
		}
	}

	if (kbsh)	// bottom b.c.
		_length_scale_obukhov[grid.gcz] = _length_scale_obukhov[grid.gcz + 1];
	if (ktsh)	// top b.c.
		_length_scale_obukhov[grid.nz - grid.gcz] = _length_scale_obukhov[grid.nz - grid.gcz - 1];
}
// -------------------------------------------------------------------------------------------- //

// Dimensionless numbers
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::prandtl_turbulent(
	T* Prandtl_turbulent,		// node: [W]

	const T* const uw_flux,		// node: [W]
	const T* const U_grad,		// node: [W]
	const T* const Tw_flux,		// node: [W]
	const T* const T_grad,		// node: [W]
	const wstGrid3d< T >& grid)
	//
	// Pr-turbulent = (u'w' * dT/dz) / (T'w' * dU/dz)
	//
	// [u'w', dU/dz, T'w', dT/dz] have to be known at all [W] nodes, excluding walls
{
	// define -k shifts to prevent division by zero //
	const int kbsh = (grid.mpi_com.rank_z == 0) ? 1 : 0;
	const int ktsh = (grid.mpi_com.rank_z == grid.mpi_com.size_z - 1) ? 1 : 0;

	int k;
#pragma omp parallel for private(k) shared(Prandtl_turbulent)
	for (k = grid.gcz + kbsh; k <= grid.nz - grid.gcz - ktsh; k++) {	// all [W] nodes, excluding walls
		Prandtl_turbulent[k] = (uw_flux[k] * T_grad[k]) / (Tw_flux[k] * U_grad[k]);
	}

	if (kbsh)	// bottom b.c.
		Prandtl_turbulent[grid.gcz] = Prandtl_turbulent[grid.gcz + 1];
	if (ktsh)	// top b.c.
		Prandtl_turbulent[grid.nz - grid.gcz] = Prandtl_turbulent[grid.nz - grid.gcz - 1];
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const U_grad,	// node: [W]
	const T* const T_grad,	// node: [W]

	const T c_Richardson, 
	const nse_const3d::axisType axis, const wstGrid3d< T >& grid)
	//
	// Ri-gradient = Ri(b) * [(dT/dz) / (dU/dz)^2]
	//
	// [dU/dz, dT/dz] have to be known at all [W] nodes, including walls
{
	if (axis == nse_const3d::axisZ) {

		int k;
#pragma omp parallel for private(k) shared(Richardson_gradient)
		for (k = grid.gcz; k <= grid.nz - grid.gcz; k++) { // all [W] nodes
			Richardson_gradient[k] = c_Richardson * (T_grad[k] / (U_grad[k] * U_grad[k]));
		}
		return;
	}
	if (axis == nse_const3d::axisYZ) {

		int j, k, idx;
#pragma omp parallel for private(j, k, idx) shared(Richardson_gradient)
		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;
				Richardson_gradient[idx] = c_Richardson * (T_grad[idx] / (U_grad[idx] * U_grad[idx]));
			}
		return;
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const uw_flux,	// node: [W]
	const T* const U_grad,	// node: [W]
	const T* const Tw_flux,	// node: [W]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	// Ri-flux = Ri(b) * [T'w' / (u'w' * dU/dz)]
	//
	// [u'w', dU/dz, T'w'] have to be known at all [W] nodes, excluding walls
{
	// define -k shifts to prevent division by zero //
	const int kbsh = (grid.mpi_com.rank_z == 0) ? 1 : 0;
	const int ktsh = (grid.mpi_com.rank_z == grid.mpi_com.size_z - 1) ? 1 : 0;

	int k;
#pragma omp parallel for private(k) shared(Richardson_flux)
	for (k = grid.gcz + kbsh; k <= grid.nz - grid.gcz - ktsh; k++) {	// all [W] node, excluding walls
		Richardson_flux[k] = c_Richardson * (Tw_flux[k] / (uw_flux[k] * U_grad[k]));
	}

	if (kbsh)	// bottom b.c.
		Richardson_flux[grid.gcz] = Richardson_flux[grid.gcz + 1];
	if (ktsh)	// top b.c.
		Richardson_flux[grid.nz - grid.gcz] = Richardson_flux[grid.nz - grid.gcz - 1];
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const TKE_iso_dissipation,	// node: [C]
	const T* const T_grad,				// node: [W]
	const T c_Richardson, const T c_kinematic_viscosity,
	const wstGrid3d< T >& grid)
	//
	// Re-buoyancy = e(TKE) / (nu * N^2) = e(TKE) / (nu*Ri(b)*dT/dz)
	//
	// [dT/dz] has to be known at all [W] nodes, including walls
{
	int k;
	if (fabs(c_Richardson) > 0) {

#pragma omp parallel for private(k) shared(Reynolds_buoyancy)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			Reynolds_buoyancy[k] = (-TKE_iso_dissipation[k]) /
				(c_kinematic_viscosity * c_Richardson * (T)0.5 * (T_grad[k] + T_grad[k + 1]));
		}
	}
	else
	{

#pragma omp parallel for private(k) shared(Reynolds_buoyancy)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			Reynolds_buoyancy[k] = (T) 0.0;
		}
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const u_TKE,				// node: [C]
	const T* const v_TKE,				// node: [C]
	const T* const TKE_iso_dissipation,	// node: [C]
	const T* const T_grad,				// node: [W]
	const T c_Richardson,
	const wstGrid3d< T >& grid)
	//
	// Fr-horizontal = (1/([Ri(b)*dT/dz]^1/2))*[0.5*e(TKE)/(E(u)+E(v))]
	//
	// [dT/dz] has to be known at all [W] nodes, including walls
{
	int k;
	if (fabs(c_Richardson) > 0) {

#pragma omp parallel for private(k) shared(Froude_horizontal)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			Froude_horizontal[k] =
				((T)1.0 / sqrt(c_Richardson * (T)0.5 * (T_grad[k] + T_grad[k + 1]))) *
				((-(T)0.5 * TKE_iso_dissipation[k]) / (u_TKE[k] + v_TKE[k]));
		}
	}
	else
	{

#pragma omp parallel for private(k) shared(Froude_horizontal)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			Froude_horizontal[k] = (T) 0.0;
		}
	}
}
// -------------------------------------------------------------------------------------------- //

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

	const T* const TKE_iso_dissipation,	// node: [C]
	const T* const TVA_iso_dissipation,	// node: [C]
	const T* const T_grad,				// node: [W]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	// e(TPE) / e(TKE) = Ri(b) * [e(TTE) / (e(TKE)*(dT/dz))]
	//
	// [dT/dz] has to be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(_mixing_efficiency)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		_mixing_efficiency[k] = c_Richardson *
			(TVA_iso_dissipation[k] / ((T)0.5 * TKE_iso_dissipation[k] * (T_grad[k] + T_grad[k + 1])));
	}
}
// -------------------------------------------------------------------------------------------- //

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

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

// pressure-strain models
// -------------------------------------------------------------------------------------------- //
template< typename T >
void nse::Rotta_model(
	T* u_Rotta, T* v_Rotta, T* w_Rotta,		// node: [C]
	T* uw_Rotta,							// node: [C]

	const T* const TKE,														// node: [C]
	const T* const TKE_iso_dissipation,										// node: [C]
	const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]
	const T* const uw_flux,													// node: [W]
	
	const T* const u_TKE_exchange,	// node: [C]
	const T* const v_TKE_exchange,	// node: [C] 
	const T* const w_TKE_exchange,	// node: [C]
	const T* const P2Suw_turb_c,	// node: [�]
	const wstGrid3d< T >& grid)
	//
	// Rotta "return-to-isotropy" model constants estimation
	// Qii = (2/3)*(C(r)/t(T))*(TKE - 3*TKE(i))
	//
	// [u'w'] has to be known at all [W] nodes, including walls
{
	int k;
#pragma omp parallel for private(k) shared(u_Rotta, v_Rotta, w_Rotta, uw_Rotta)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		u_Rotta[k] = (T)3.0 * (u_TKE_exchange[k] / (-TKE_iso_dissipation[k])) *
			(TKE[k] / (TKE[k] - (T)3.0 * u_TKE[k]));
		v_Rotta[k] = (T)3.0 * (v_TKE_exchange[k] / (-TKE_iso_dissipation[k])) *
			(TKE[k] / (TKE[k] - (T)3.0 * v_TKE[k]));
		w_Rotta[k] = (T)3.0 * (w_TKE_exchange[k] / (-TKE_iso_dissipation[k])) *
			(TKE[k] / (TKE[k] - (T)3.0 * w_TKE[k]));

		uw_Rotta[k] = (P2Suw_turb_c[k] / (-TKE_iso_dissipation[k])) *
			(TKE[k] / (-(T)0.5 * (uw_flux[k] + uw_flux[k + 1])));
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::RDT_model(T* u_RDT, T* v_RDT, T* w_RDT,	// node: [C] node
	T* uw_RDT,										// node: [C] node

	const T* const TKE_production,											// node: [C]
	const T* const W2_w,													// node: [W]
	const T* const U,														// node: [C]
	const T* const W,														// node: [W]

	const T* const u_TKE_exchange,	// node: [C]
	const T* const v_TKE_exchange,	// node: [C] 
	const T* const w_TKE_exchange,	// node: [C]
	const T* const P2Suw_turb_c,	// node: [�]
	const wstGrid3d< T >& grid)
	//
	// RDT "return-to-isotropy" model constants estimation
	// Qii = - (2/3)*C(p)*((3/2)*Pi - P)
	//
	// [W, W'^2] 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)
{
	T Guw;
	int k;
#pragma omp parallel for private(k, Guw) shared(u_RDT, v_RDT, w_RDT, uw_RDT)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		u_RDT[k] = -(T)3.0 * (u_TKE_exchange[k] / ((T) 2.0 * TKE_production[k]));
		v_RDT[k] = (T)3.0 * (v_TKE_exchange[k] / (TKE_production[k]));
		w_RDT[k] = (T)3.0 * (w_TKE_exchange[k] / (TKE_production[k]));

		Guw =
			(W2_w[k] - W[k] * W[k]) *
			(U[k] - U[k - 1]) * grid.dzmi[k] +

			(W2_w[k + 1] - W[k + 1] * W[k + 1]) *
			(U[k + 1] - U[k]) * grid.dzmi[k + 1];

		uw_RDT[k] = P2Suw_turb_c[k] / Guw;
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::Rotta_RDT_model(T* Rotta_RDT_e, T* Rotta_RDT_p,	// node: [C]

	const T* const TKE,														// node: [C]
	const T* const TKE_iso_dissipation,										// node: [C]
	const T* const TKE_production,											// node: [C]
	const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]
	
	const T* const u_TKE_exchange,	// node: [C]
	const T* const v_TKE_exchange,	// node: [C] 
	const T* const w_TKE_exchange,	// node: [C]
	const wstGrid3d< T >& grid)
	//
	// Rotta-RDT "return-to-isotropy" model constants estimation
	// Qii = (2/3)*(C(r)/t(T))*(TKE - 3*TKE(i)) - (2/3)*C(p)*((3/2)*Pi - P)
	//  (*) 2 equations suffice for constant determination
	//
{
	T Eu, Ev, Gu, Gv, Qu, Qv;
	int k;
#pragma omp parallel for private(k, Eu, Ev, Gu, Gv, Qu, Qv) \
	shared(Rotta_RDT_e, Rotta_RDT_p)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {

		// using Quu, Qvv equations only //
		Eu = (T) 3.0 * u_TKE[k] - TKE[k];
		Ev = (T) 3.0 * v_TKE[k] - TKE[k];

		Gu = (T) 2.0 * TKE_production[k];
		Gv = -TKE_production[k];

		Qu = u_TKE_exchange[k];
		Qv = v_TKE_exchange[k];

		Rotta_RDT_e[k] = (T) 3.0 * (TKE[k] / (-TKE_iso_dissipation[k])) *
			((Gv * Qu - Gu * Qv) / (Ev * Gu - Eu * Gv));
		Rotta_RDT_p[k] = (T) 3.0 *
			((Eu * Qv - Ev * Qu) / (Ev * Gu - Eu * Gv));
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::Rotta_buoyancy_model(
	T* Rotta_buoyancy_e, T* Rotta_buoyancy_b,	// node: [C]

	const T* const TKE,														// node: [C]
	const T* const TKE_iso_dissipation,										// node: [C]
	const T* const TKE_heat_flux,											// node: [C]
	const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]
	
	const T* const u_TKE_exchange,	// node: [C]
	const T* const v_TKE_exchange,	// node: [C] 
	const T* const w_TKE_exchange,	// node: [C]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	// Rotta-buoyancy "return-to-isotropy" model constants estimation
	// Qii = (2/3)*(C(r)/t(T))*(TKE - 3*TKE(i)) - (2/3)*C(b)*((3/2)*Bi - B)
	//  (*) 2 equations suffice for constant determination
	//
{
	T Eu, Ew, Bu, Bw, Qu, Qw;
	int k;

	if (fabs(c_Richardson) > 0) {
#pragma omp parallel for private(k, Eu, Ew, Bu, Bw, Qu, Qw) \
	shared(Rotta_buoyancy_e, Rotta_buoyancy_b)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {

			// using Quu, Qww equations only //
			Eu = (T) 3.0 * u_TKE[k] - TKE[k];
			Ew = (T) 3.0 * w_TKE[k] - TKE[k];

			Bu = -TKE_heat_flux[k];
			Bw = (T) 2.0 * TKE_heat_flux[k];

			Qu = u_TKE_exchange[k];
			Qw = w_TKE_exchange[k];

			Rotta_buoyancy_e[k] = (T) 3.0 * (TKE[k] / (-TKE_iso_dissipation[k])) *
				((Bw * Qu - Bu * Qw) / (Ew * Bu - Eu * Bw));
			Rotta_buoyancy_b[k] = (T) 3.0 *
				((Eu * Qw - Ew * Qu) / (Ew * Bu - Eu * Bw));
		}
	}
	else
	{
#pragma omp parallel for private(k) shared(Rotta_buoyancy_e, Rotta_buoyancy_b)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			Rotta_buoyancy_e[k] = (T) 0;
			Rotta_buoyancy_b[k] = (T) 0;
		}
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::RDT_buoyancy_model(
	T* RDT_buoyancy_p, T* RDT_buoyancy_b,		// node: [C]

	const T* const TKE_production,											// node: [C]
	const T* const TKE_heat_flux,											// node: [C]

	const T* const u_TKE_exchange,	// node: [C]
	const T* const v_TKE_exchange,	// node: [C] 
	const T* const w_TKE_exchange,	// node: [C]
	const T c_Richardson, const wstGrid3d< T >& grid)
	//
	// RDT-buoyancy "return-to-isotropy" model constants estimation
	// Qii = - (2/3)*C(p)*((3/2)*Pi - P) - (2/3)*C(b)*((3/2)*Bi - B)
	//  (*) 2 equations suffice for constant determination
	//
{
	T Gu, Gw, Bu, Bw, Qu, Qw;
	int k;

	if (fabs(c_Richardson) > 0) {
#pragma omp parallel for private(k, Gu, Gw, Bu, Bw, Qu, Qw) \
	shared(RDT_buoyancy_p, RDT_buoyancy_b)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {

			// using Quu, Qww equations only //
			Gu = (T) 2.0 * TKE_production[k];
			Gw = -TKE_production[k];

			Bu = -TKE_heat_flux[k];
			Bw = (T) 2.0 * TKE_heat_flux[k];

			Qu = u_TKE_exchange[k];
			Qw = w_TKE_exchange[k];

			RDT_buoyancy_p[k] = (T) 3.0 *
				((Bw * Qu - Bu * Qw) / (Gw * Bu - Gu * Bw));
			RDT_buoyancy_b[k] = (T) 3.0 *
				((Gu * Qw - Gw * Qu) / (Gw * Bu - Gu * Bw));
		}
	}
	else
	{
#pragma omp parallel for private(k) shared(RDT_buoyancy_p, RDT_buoyancy_b)
		for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
			RDT_buoyancy_p[k] = (T)0;
			RDT_buoyancy_b[k] = (T)0;
		}
	}
}
// -------------------------------------------------------------------------------------------- //

template< typename T >
void nse::Rotta_TPE_model(
	T* u_Rotta_TPE, T* v_Rotta_TPE, T* w_Rotta_TPE,		// node: [C]

	const T* const TKE,						// node: [C]
	const T* const TPE,						// node: [C]
	const T* const TKE_iso_dissipation,		// node: [C]
	const T* const TPE_iso_dissipation,		// node: [C]
	const T* const u_TKE, const T* const v_TKE, const T* const w_TKE,		// node: [C]

	const T* const u_TKE_exchange,	// node: [C]
	const T* const v_TKE_exchange,	// node: [C] 
	const T* const w_TKE_exchange,	// node: [C]
	const wstGrid3d< T >& grid)
	//
	// Rotta-TPE "return-to-isotropy" model constants estimation
	// Qii = (2/3)*(C(r)/t(TKE+TPE))*((TKE+TPE) - 3*(TKE(i) + TPE[i,3]))
	//
{
	int k;
#pragma omp parallel for private(k) shared(u_Rotta_TPE, v_Rotta_TPE, w_Rotta_TPE)
	for (k = grid.gcz; k < grid.nz - grid.gcz; k++) {
		u_Rotta_TPE[k] = (T)3.0 * (u_TKE_exchange[k] / (-TKE_iso_dissipation[k] - TPE_iso_dissipation[k])) *
			((TKE[k] + TPE[k]) / ((TKE[k] + TPE[k]) - (T)3.0 * u_TKE[k]));
		v_Rotta_TPE[k] = (T)3.0 * (v_TKE_exchange[k] / (-TKE_iso_dissipation[k] - TPE_iso_dissipation[k])) *
			((TKE[k] + TPE[k]) / ((TKE[k] + TPE[k]) - (T)3.0 * v_TKE[k]));
		w_Rotta_TPE[k] = (T)3.0 * (w_TKE_exchange[k] / (-TKE_iso_dissipation[k] - TPE_iso_dissipation[k])) *
			((TKE[k] + TPE[k]) / ((TKE[k] + TPE[k]) - (T)3.0 * (w_TKE[k] + TPE[k])));
	}
}
// -------------------------------------------------------------------------------------------- //