How to calculate gradient term?

[YH_wcj]

I use wrf-python to extract the variable wind u = wrf.getvar(nc,"ua") and vertically interpolate it to pressure levels. I now wish to compute the horizontal gradients of ua. The horizontal grid spacing configured in my namelist is 27 km.
Is it valid to use dx and dy to perform centered differencing for this calculation?

 

[William.Hatheway]

I use wrf-python to extract the variable wind u = wrf.getvar(nc,"ua") and vertically interpolate it to pressure levels. I now wish to compute the horizontal gradients of ua. The horizontal grid spacing configured in my namelist is 27 km.
Is it valid to use dx and dy to perform centered differencing for this calculation?

Yes, centered differencing is valid, but I would use the grid information already written to the WRF output rather than hard-coding 27 km.

Since
ua = wrf.getvar(nc, "ua")
returns U destaggered to the mass grid, you can use RDX, RDY, MAPFAC_MX, and MAPFAC_MY, assuming those variables are available in your WRF output.

∂u/∂x ≈ 0.5 × RDX × MAPFAC_MX × (u[i+1,j] - u[i-1,j])
∂u/∂y ≈ 0.5 × RDY × MAPFAC_MY × (u[i,j+1] - u[i,j-1])

where
RDX = (Δx)^-1<br>RDY = (Δy)^-1

For a 27 km grid:
RDX = RDY = (27,000 m)^-1

You could also write this as:
∂u/∂x ≈ MAPFAC_MX × (u[i+1,j] - u[i-1,j]) / 54,000 m
∂u/∂y ≈ MAPFAC_MY × (u[i,j+1] - u[i,j-1]) / 54,000 m

One caveat is that this assumes RDX, RDY, MAPFAC_MX, and MAPFAC_MY are present in the WRFOUT file. You can check the file header first with NetCDF commands such as:

ncdump -h wrfout_d01_YYYY-MM-DD_HH:MM:SS | grep -E "RDX|RDY|MAPFAC_MX|MAPFAC_MY"

or simply:

ncdump -h wrfout_d01_YYYY-MM-DD_HH:MM:SS

WRF does not appear to output DUDX or DUDY directly, but if these grid-metric variables are available, they provide what you need to calculate the gradients.

Hope this helps.

 

[YH_wcj]

Thanks a lot, you’re a legend!
Now, I'm using the method you methioned to difference.
And if second derivative , can I use?

mx = MAPFAC_MX ; Because I gain that "mx = my = m" from WRF technote
my = MAPFAC_MY;
m = MAPFAC_M;

 

[William.Hatheway]

Yes, with one important caveat.

You are correct that the WRF Technical Note states that for the isotropic projections, Lambert conformal, polar stereographic, and Mercator,

mx = my = m

So for those projections, `MAPFAC_MX`, `MAPFAC_MY`, and `MAPFAC_M` should represent the same map-scale factor on the mass grid. This does not apply to the latitude-longitude projection.

If the equation in your image is the usual second centered difference multiplied by `m^2`,

∂²u/∂x² ≈ m² × (u[i+1,j] - 2u[i,j] + u[i-1,j]) / (Δx)²

that is valid when `m` can be treated as approximately constant across the finite-difference stencil.

More generally, because the first physical derivative is

∂u/∂x = m × ∂u/∂X

the second derivative should be treated as

∂²u/∂x² = m × ∂/∂X (m × ∂u/∂X)

so if `m` varies appreciably across the domain, you should not simply square the map factor and ignore its spatial variation.

For `ua`, since it is on the mass grid, I would continue using `MAPFAC_MX` for the x direction and `MAPFAC_MY` for the y direction. If you are using Lambert, polar stereographic, or Mercator, then `mx = my = m` is appropriate.

You can also confirm the projection and the available map-factor variables directly from your WRFOUT header, for example:

ncdump -h wrfout_d01_* | grep -E "MAP_PROJ|MAPFAC_M|MAPFAC_MX|MAPFAC_MY"

If `u` is in `m s^-1`, the second horizontal derivative will have units of `m^-1 s^-1`.

Hope this helps.

 

[YH_wcj]

Thanks, I finally get it now.

 

[William.Hatheway]

@kwerner confirm i'm right on this