I'm working on a diagnostics module to output resolved-scale and subgrid-scale heat and momentum fluxes for a real case WRF-LES. All diagnostics calulations are called fomr the module_diagnostics_driver.F.

I'm using the 3D-TKE based turbulence parameterization (km_opt=2). The subgrid-scale heat and momentum fluxes (in K m/s and m²/s²) are calculated following Deardorff using the eddy viscosity and diffusivity provided by WRF, e.g.

SGS vertical heat flux W_THETA = -xkhv * dTHETA/dZ

or

SGS horizontal heat flux U_THETA = -xkhh * dTHETA/dX

I'm using code taken from the subroutine "horizontal_diffusion_s" from module_diffusion.F to calculate horizontal gradients (H1 respectively H2 in the code) and subroutine "vertical_diffusion_s" to get the vertical gradients (H3 in the code) inside my diagnostics module.

**Is this ok for that purpose? Are there any subroutines readily available to calculate d/dx, d/dy d/dz of three dimensional variables in physical space?**

After calculating gradients for U,V,W and pot. temperature THETA in all three dimension these are multiplied with the appropriate exchange coefficients xk[hm][hv] to get the respective subgrid-scale turbulent fluxes.

However, if a compare the parameterized surface heat flux (converted from W/m² to Km/s) to my derived SGS turbulent heat flux on the first level and further up there is a difference of over three magnitudes with the surface heat flux being larger. In my opinion, it is physical that the surface heat flux is larger than the SGS fluxes. However, I think that the SGS fluxes are way too low on model levels.

**Do I use the exchange coefficients xkmh, xkmv, xkhh, xkhv in the appropriate way or do they needed to be coupled/uncoupled to the dry airmass or something similar because i use them "outside" of solve_em in my diagnsotics module?**

Any hints on how to diagnose the SGS turbulent fluxes will be greatly appreciated!