# Calculation of subgrid-scale heat and momentum fluxes

This post was from a previous version of the WRF&MPAS-A Support Forum. New replies have been disabled and if you have follow up questions related to this post, then please start a new thread from the forum home page.

#### SIGSEV

##### New member
Hello,
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!

Hi,
Your assumptions look okay. x and k variables are not coupled. Make sure the theta variable is an uncoupled one. Can you tell me what kinds of values you are seeing for the surface flux and the lowest level?

hi,
To get potential temperature in my diagnostics module I use
Code:
``````IF ( config_flags%use_theta_m == 0 ) THEN
PT = grid%t_2 + 300.0
ELSE
PT(:,:,:) = (grid%t_2(:,:,:) + 300.0)/(1 + rvovrd * grid%moist(:,:,:,P_QV))
END IF``````
The calculation of potential temperature from moist potential temperature is taken from another postition in the WRF code.

In my simulation (dx=40m) I get the following fluxes, at representative grid points during night-time:

surface heat flux HFX = -6.15 W/m² = - 5.2 E-3 Km/s
1st model level SGS heat flux (10m AGL ) = - 1.96 E-5 Km/s
2nd model level SGS heat flux (20m AGL ) = - 1.57 E-6 Km/s
3rd model level SGS heat flux (30m AGL ) = - 1.11 E-9 Km/s

The parameterized heat flux drop on the order of two to four magnitudes between the surface and the 1st model level throughout my domain.

While a positive heat flux should affect many levels by forcing an unstable profile, a negative one like this is
not expected to directly influence more than the surface layer (and perhaps the lowest level due to some background interference). If you're interested in determining the effects at higher levels, then you likely will want to look at an unstable time.

Hello!
I'm am currently testing my SGS heat and momentum flux computation with an convective case, the idealized LES (em_les) and vertical fluxes and flux profiles look ok.
I noticed that the calculation of horizontal gradients (dVAR/dx, dVAR/dy) fails for potential temperature while it delivers valid results for U and V wind components. so far I was unable to find the bug and would be very happy if you could provide some feedback.

To calculate horizontal gradients I use the following subroutine (adapted from dyn_em/module_diffusion_em.F and a similar one for dVARdY based on the the calculation of H2 in module_diffusion_em.F)
Code:
``````  SUBROUTINE calc_dVARdX(dVARdX, mu, config_flags, var, &
msftx, msfty, msfux, msfuy,               &
msfvx, msfvy, rdx, rdy,                   &
fnm, fnp, cf1, cf2, cf3,                  &
zx, zy, rdz, rdzw, dnw, dn,               &
ids, ide, jds, jde, kds, kde,             &
ims, ime, jms, jme, kms, kme,             &
its, ite, jts, jte, kts, kte)

IMPLICIT NONE

TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte
REAL , INTENT(IN   )           ::        cf1, cf2, cf3
REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp
REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::     dn
REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    dnw
REAL , DIMENSION( ims:ime, jms:jme) ,      INTENT(IN   ) ::   msfux
REAL , DIMENSION( ims:ime, jms:jme) ,      INTENT(IN   ) ::   msfuy
REAL , DIMENSION( ims:ime, jms:jme) ,      INTENT(IN   ) ::   msfvx
REAL , DIMENSION( ims:ime, jms:jme) ,      INTENT(IN   ) ::   msfvy
REAL , DIMENSION( ims:ime, jms:jme) ,      INTENT(IN   ) ::   msftx
REAL , DIMENSION( ims:ime, jms:jme) ,      INTENT(IN   ) ::   msfty
REAL , DIMENSION( ims:ime, jms:jme) ,      INTENT(IN   ) ::   mu
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::  rdz, rdzw
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::    var, zx, zy
REAL ,                                        INTENT(IN   ) ::    rdx, rdy
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dVARdX
! Local data
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)    ::     H1avg, H2avg, H1, H2
REAL    :: mrdx, mrdy
REAL    :: tmpzx, tmpzy, tmpzeta_z, rdzu, rdzv
INTEGER :: ktes1,ktes2

ktf=MIN(kte,kde-1)
ktes1=kte-1
ktes2=kte-2

i_start = its
i_end   = MIN(ite,ide-1)
j_start = jts
j_end   = MIN(jte,jde-1)

IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end   = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end   = MIN(jde-2,jte)
IF ( config_flags%periodic_x ) i_start = its
IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)

! H1 = partial var over partial x
DO j = j_start, j_end
DO k = kts+1, ktf
DO i = i_start, i_end + 1
H1avg(i,k,j)=0.5*(fnm(k)*(var(i-1,k  ,j)+var(i,k  ,j))+  &
fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j)))
END DO
END DO
END DO

DO j = j_start, j_end
DO i = i_start, i_end + 1
H1avg(i,kts  ,j)=0.5*(cf1*var(i  ,1,j)+cf2*var(i  ,2,j)+ &
cf3*var(i  ,3,j)+cf1*var(i-1,1,j)+  &
cf2*var(i-1,2,j)+cf3*var(i-1,3,j))

H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
var(i-1,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1))
END DO
END DO

DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end + 1
tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
H1(i,k,j)=msfuy(i,j)*(                      &
rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx*         &
(H1avg(i,k+1,j)-H1avg(i,k,j))*rdzu
END DO
END DO
END DO

DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
dVARdX(i,k,j)= 0.5*(H1(i+1,k,j) + H1(i,k,j))
END DO
END DO
END DO``````

This subroutine is called within my diagnostics module (e.g. for potential temperature) using:
Code:
``````   CALL calc_dVARdX(dPTdX, grid%mu_2, config_flags, PT, .false., &
grid%msftx, grid%msfty, grid%msfux, grid%msfuy, &
grid%msfvx, grid%msfvy, grid%rdx, grid%rdy, &
grid%fnm, grid%fnp, grid%cf1, grid%cf2, grid%cf3, &
grid%zx, grid%zy, grid%rdz, grid%rdzw, grid%dnw, grid%dn, &
ids, ide, jds, jde, kds, kde,          &
ims, ime, jms, jme, kms, kme,          &
its, ite, jts, jte, kts, kte )``````
with
Code:
`` PT = (grid%t_2 + 300.0)/(1 + rvovrd * grid%moist(:,:,:,P_QV))``

When the U and V components are feed to the subroutines, resulting fields are ok but when using potential temperature I get vertical (for dTHETA/dX) and horizontal "stripes" (for dTHETA/dY):
https://ibb.co/w4yLdRC
https://ibb.co/StvsxHR

I assume this is related to some missing/wrong communication between the MPI jobs (WRF is compiled with distributed memory)?
Any advice how to overcome this issue would be greatly appreciated!

Hi,
Can you first try the "use_theta_m = 0" (&dynamics) option to check if the following works:
Try using grid%t_2 directly in the routine, instead of PT. Since t_2 is a dynamic variable, its halo has been declared so that you can access i-1, i+1 (and j-1, j+1) values (PT doesn't have halos). Let me know if this works.

1) Which version of the code you are using
2) Are you interested in the moist theta variable or the dry theta gradient?

Thanks,
Kelly

Hi Kelly,
Thank you for the advice! By using grid%t_2 instead of PT I managed to get rid of the "striped" gradients and thanks to your suggestion I got more insight into the various array dimensions!

I am using WRF V4.0.3 and want to calculate the dry potential temperature gradients in the physical space because I want to calculate the resolved heat fluxes u'theta', v'theta' and w'theta'

Hi,
That is great news! Thank you for letting us know.
Are you content now to move on, or do you still need our help?

Hi Kelly,
does using the dry potential temperature introduce any problesm/issues, especially when doing LES? I'm asking because moist potential temperature is used as the default prognostic value and I'm sure there is some background to that decision?

If I want to use use_theta_m = 1, is it be possible to create a dry potential temperature field with halos from the "memory" arrays of moist potential temperature and specific humidity somehow like this?

Code:
``````! local variable to derive gradients in dry potential temperature
REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: PT

IF ( config_flags%use_theta_m == 0 ) THEN
PT = grid%t_2 + 300.0
ELSE
DO j = jms, jme
DO k = kms, kme
DO i = ims, ime
PT(i,k,j) = ( grid%t_2(i,k,j) + 300.0 ) / ( 1 + rvovrd * grid%moist(i,k,j,P_QV) )
END DO
END DO
END DO
END IF``````

Hi,
Can you tell me from where your new subroutine is called?
Thanks!

Hi,
I created a new diagnostics module (phys/module_diag_LESfluxes.F).
Inside module_diag_LESfluxes.F I have the primary subroutine which is called from module_diagnostics_driver.F
Code:
``````SUBROUTINE diag_les_fluxes(grid, u, v, w, config_flags     &
,ids,ide, jds,jde, kds,kde     &
,ims,ime, jms,jme, kms,kme     &
,its,ite, jts,jte, kts,kte     &
,num_tiles                     )
...
...
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: PT ! full dry potential temperature
...
...
IF ( config_flags%use_theta_m == 0 ) THEN
PT = grid%t_2 + 300.0
ELSE
DO j = jms, jme
DO k = kms, kme
DO i = ims, ime
PT(i,k,j) = ( grid%t_2(i,k,j) + 300.0 ) / ( 1 + rvovrd * grid%moist(i,k,j,P_QV) )
END DO
END DO
END DO
END IF``````

and separate subroutines to determine the gradients in x,y,z direction of
Code:
``````SUBROUTINE calc_dVARdX(....)
...

SUBROUTINE calc_dVARdY(....)
...

SUBROUTINE calc_dVARdZ(....)
...``````

I tested this and do not see "striped" gradients for dry potential temperature. So I think it is okay to get PT in that way and it should include all halos like grid%t_2 does?

thanks!

Hi,
Try passing th_phy to your subroutine and using that instead of PT. This is an array for dry theta, but it is computed before the theta_m is updated. So you still need to do the calculation to get an updated dry theta, but this array should be good for that. You should also add the following call for a halo right after you compute th_phy:

#ifdef DM_PARALLEL
# include "HALO_EM_SBM.inc"
#endif

Since you are borrowing an array for dry theta, it would be better if you move your call to the subroutine later in module_diagnostics_driver.F, so that others may not be affected.

Hi Kelly,
sorry for digging up this old post but I have a follow up question: You wrote

kwerner said:
Try passing th_phy to your subroutine and using that instead of PT. This is an array for dry theta, but it is computed before the theta_m is updated. So you still need to do the calculation to get an updated dry theta, but this array should be good for that.

What do you exactly mean with "do the calculation"?
At the end of the module_diagnostics_driver.F, I can calculate dry theta based on moist theta and humidity and store the result in the array for "th_phy" ?
Somehow like this: grid%th_phy(i,k,j) = ( grid%t_2(i,k,j) + 300.0 ) / ( 1 + rvovrd * grid%moist(i,k,j,P_QV) )

Thanks!

Hi,
Yes that will likely work. I would do a short test to ensure that you are getting updates for dry theta, and the results you expect to make sure.

Hi all,

Is there a standard way of diagonalize subgrid-scale heat and momentum fluxes
in WRF?

Xiang-Yu

Xiang-yu,
I don' think we have a 'standard way' to calculate sub grid heat and momentum fluxes. Such kind of computation must be conducted based on the specific schemes you use to run WRF.

Thank you for clarifying this, Ming.

Best,

Xiang-Yu