fvScalarMatrix hEqn
(
    fvm::ddt(rho,h)
  + fvm::div(mdotf, h) 
  - fvm::laplacian(turbulence->alphaEff(), h)
 + fvm::SuSp(-(fvc::div(mdotf)+fvc::ddt(rho)),h)
 ==
    fvc::div(mdotf/fvc::interpolate(rho)*fvc::interpolate(p))
  - p*fvc::div(mdotf/fvc::interpolate(rho))
);	


hEqn.relax();	 
hEqn.solve();

h.correctBoundaryConditions(); 


//h = dimensionedScalar("h0", dimVelocity*dimVelocity,332310) - magSqr(U)*0.5;

thermo->correct();

gammaGas =  thermo->Cp()/ thermo->Cv();
gammaGas.correctBoundaryConditions();

RGas =  thermo->Cp() -  thermo->Cv();
RGas.correctBoundaryConditions();




