WRF relative topography... problem with the difference

Question

I currently struggle with the following problem. I want to calculate the relative topography. That means I need the difference between two geopotential heights. In my case H500-H1000 (the number means the pressure level). In case of WRF the pressure levels are topography following... so there is only a pressure level at the corrosponding hight. In that case, you get trouble with mountains.

WRF Python has a methode called interplevl... that works fine if you have the pressure. But how I can solve that with the 1000 hpa level?

My first idea was:

geo_1000 = vinterp(wrfin=nc_meteo_wrf,
                field=geopt,
                vert_coord='pressure',
                interp_levels=[1000.],
                method='cat',
                extrapolate = True,
                squeeze=True)

But the result show also the topography and is independend from the time step...

Also that didn't work:

def poisson_fill(grid, tol=1e-4, max_iter=10000):
    """
    Fill missing values in the grid using Poisson equation via relaxation.
    Parameters:
    grid : 2D numpy array
    Input grid with NaN values indicating missing data.
    tol : float
    Tolerance for convergence.
    max_iter : int
    Maximum number of iterations.
    Returns:
    filled_grid : 2D numpy array
    Grid with missing values filled.
    """
    filled_grid = np.where(np.isnan(grid), 0, grid)
    mask = np.isnan(grid)
    for _ in range(max_iter):
        old_grid = filled_grid.copy()
        laplace_grid = laplace(filled_grid)
        filled_grid[mask] = 0.25 * laplace_grid[mask]
        if np.nanmax(np.abs(filled_grid - old_grid)) < tol:
           break
    return filled_grid

            geo_1000 = vinterp(nc_meteo_wrf, field=geopt, vert_coord="pressure",
                                interp_levels=[1000], extrapolate=True)
            geo_1000 = to_np(geo_1000[0])  # Konvertieren in numpy-Array und die erste Ebene extrahieren
            geopt_1000_filled = poisson_fill(geo_1000)

            geo500_1000_difference = to_np(geo_500) - geopt_1000_filled

Maybe there is only a little error?

Answer 1

Soo.. I found an easy answer.

I calculated the 850 hPa field, after that i used the barometic equation to calculate the the 1000 hPa field. And it works. Here is the code:

# Extrapolation zu 1000 hPa mithilfe der barometrischen Höhenformel def extrapolate_to_surface(z_low, z_high, p_low, p_high, target_pressure):
    # Berechne die Skalenhöhe
    scale_height = (z_high - z_low) / np.log(p_low / p_high)
    # Extrapoliere auf den Ziel-Druck
    z_target = z_low + scale_height * np.log(p_low / target_pressure)
    return z_target

and after that:

        ph = getvar(nc_meteo_wrf, "PH", timeidx=desired_time_index)
        phb = getvar(nc_meteo_wrf, "PHB", timeidx=desired_time_index)
        geopot_phphb = ph + phb
        #z_stag = geopot_phphb/9.81
        z_stag = geopot_phphb
        # Mittelung der gestaggerten Ebenen auf nicht-gestaggerte Ebenen
        z_geo = 0.5 * (z_stag[:-1, :, :] + z_stag[1:, :, :])
        geo_850 = interplevel(z_geo, p, 850)


        target_pressure = 1000
        geo_1000 = extrapolate_to_surface(geo_850, geo_500, 850, 500, target_pressure)
        #geo_1000 = geopt[0]
        geo500_1000_difference = to_np(geo_500) - geo_1000