Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

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.

S

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,
thanks for your quick response!
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.

Please also let me know:
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,
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.
Click to expand...

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.
 
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.
 
You must log in or register to reply here.
Top