<p>My task is to compute isentropic potential vorticity on a grid.
All data (ECMWF ens) are given on a hybrid sigma pressure grid.
Since I really cannot find reference, I assume (!) that relative vorticity (vo, id=138) is evaluated on model levels and therefore is useless for me.</p>

<p>My approach (for each grid point):</p>

<ol>
<li><p>Compute theta-gradient, using the three model dimensions, converting vertical and horizontal distances to meters.</p></li>
<li><p>Define two vectors using three constraints each:</p>

<ul>
<li>Scalar product with theta gradient vanishes.</li>
<li>First vector has no meridional, second one no zonal component.</li>
<li>Distance to adjacent horizontal grid point defines magnitude of 
meridional/zonal component.</li>
</ul></li>
<li>These two vectors now span a plane tangential to the isentropic surface at the current grid point. They point to positions only vertically displaced from the horizontal neighbours of the current grid point.</li>
<li>Values of u and v are interpolated - again, using meters (geopotential height) - to the locations given by the tangent vectors.</li>
<li>The geometrical length of these tangent vectors is determined and used together with the interpolated values of u and v to compute isentropic relative vorticity.</li>
<li>Static stability is computed using values of pressure and theta at current point and vertical neighbour.</li>
<li>Using gravitational acceleration and the conversion factor of 10^-6, I use the equation given in the famous book by Wallace &amp; Hobbs for Potential Vorticity: $g(f + \eta_{\theta})\frac{\partial \theta}{\partial p}$</li>
</ol>

<p>My problem:</p>

<p>The result does not look like what I expect.
Values range from around -30 to +50 pvu throughout the atmosphere, while the 1.5pvu surface clearly does not depict the tropopause. It looks more like an accident. I am using "omp" to increase computation speed, but things look the same if I don't use it, too, so it cannot be index confusion.</p>

<h3>My question:</h3>

<p>What could be wrong?
Can you spot systematic errors? Maybe I got blind with time.
Do you know reference for how relative vorticity is evaluated in the ECMWF model? (I tried, really...)</p>

<p>To be more specific:
My method is to be implemented in C++ into an already existing environment to visualize NWP data in 3D in realtime.
The basic datafields - given in hybrid sigma coordinates - are read, processed, written into a derived data field, interpolated to and rendered on pressure levels as triggered by user interaction.
The program should possess the processing methods on its own, which is why I cannot use Python packages, for example.</p>

<p>Here is the source code of my approach.  If something is unclear, please ask!
Maybe somebody will spot an error.</p>

<pre><code> // ISENTROPIC POTENTIAL VORTICITY
// Input grids : 0 t, 1 q, 2 u, 3 v, 4 sfcPhi, 5 sfc t_d, 6 sfc t, 7 sfc p, 8 sfc u, 9 sfc v
else if ( (derivedVarName == "Isentropic PV (an)")
     || (derivedVarName == "Isentropic PV (fc)")
     || (derivedVarName == "Isentropic PV (ens)") )
{
    unsigned int k, j, i;
    int n, l;
#pragma omp parallel
        {
        #pragma omp for private(k, j, i, n, l)
    for (k = 1; k &lt; derivedGrid-&gt;getNumLevels() - 1; k++)
        for (j = 0; j &lt; derivedGrid-&gt;getNumLats() - 1; j++)
            for (i = 0; i &lt; derivedGrid-&gt;getNumLons() - 1; i++)
    {

        // Current grid point.
        double T000 = inputGrids.at(0)-&gt;getValue(k, j, i);
        double p000 = inputGrids.at(0)-&gt;getPressure(k, j, i);
        double q000 = inputGrids.at(1)-&gt;getValue(k, j, i);
        double theta000 = potTemp(T000, p000);

        // Upper neighbour.
        double T100 = inputGrids.at(0)-&gt;getValue(k-1, j, i);
        double p100 = inputGrids.at(0)-&gt;getPressure(k-1, j, i);
        double q100 = inputGrids.at(1)-&gt;getValue(k-1, j, i);
        double theta100 = potTemp(T100, p100);

        // x-neighbour.
        double theta001 = potTemp(inputGrids.at(0)-&gt;getValue(k, j, i+1),
                                  inputGrids.at(0)-&gt;getPressure(k, j, i+1));

        // y-neighbour.
        double theta010 = potTemp(inputGrids.at(0)-&gt;getValue(k, j+1, i),
                                  inputGrids.at(0)-&gt;getPressure(k, j+1, i));

        // (1) Compute theta gradient at current point, using geopotential thickness dZ [m] for the vertical,
        // great circle distance [m] for the horizontal distances.

            // Compute geopotential thickness between current point and upper neighbour.
            double T_v = (virtualTempFromSpecHum(T000, q000) + virtualTempFromSpecHum(T100, q100)) / 2.;
            double dPhi = geopotThickness(T_v, p100, p000);
            double dZ = geopotToZ(dPhi);

            // Compute gradient.
            double* grad = gradient(theta000,
                                    inputGrids.at(0)-&gt;getEastInterfaceLon(i),
                                    inputGrids.at(0)-&gt;getNorthInterfaceLat(j), // Current grid point.

                                    theta001,
                                    inputGrids.at(0)-&gt;getEastInterfaceLon(i+1), // xlon-neighbour.

                                    theta010,
                                    inputGrids.at(0)-&gt;getNorthInterfaceLat(j+1), // ylat-neighbour.

                                    theta100,
                                    dZ // Z-neighbour.
                                    );

            // Gradient components.
            double xgrad = grad[2];
            double ygrad = grad[1];
            double Zgrad = grad[0];

        // (2) Define tangent vectors and compute their Z-components.
        // Definition:
        // 1) Scalar product of tangent vectors with gradient vanishes.
        // 2) One tangent vector in y-z-plane, the other in x-z-plane.
        // 3) Length of horizontal components equals spacing between grid points.
        double radius_m = MetConstants::EARTH_RADIUS_km * 1000.;

            // Tangent vector in y-z-plane.
            double dx = gcDistance_deg(inputGrids.at(0)-&gt;getEastInterfaceLon(i),
                                       inputGrids.at(0)-&gt;getNorthInterfaceLat(j),
                                       inputGrids.at(0)-&gt;getEastInterfaceLon(i+1),
                                       inputGrids.at(0)-&gt;getNorthInterfaceLat(j),
                                       radius_m);

            // Z-component.
            double dZ_x = -dx * xgrad / Zgrad;

            // Tangent vector in x-z-plane.
            double dy = gcDistance_deg(inputGrids.at(0)-&gt;getEastInterfaceLon(i),
                                       inputGrids.at(0)-&gt;getNorthInterfaceLat(j),
                                       inputGrids.at(0)-&gt;getEastInterfaceLon(i),
                                       inputGrids.at(0)-&gt;getNorthInterfaceLat(j+1),
                                       radius_m);

            // Z-component.
            double dZ_y = -dy * ygrad / Zgrad;

        // Tangent vectors now point from current grid point to locations only
        // vertically displaced from x- or y-neighbour, respectively.

        // (3) Find levels whose Z values enclose geopotential height values and
        // interpolate u and v to these values.
        // The geopotential height of the current grid point is computed in the
        // x-coordinate section and reused in the y-coordinate section.

                // x-coordinate
                double vdx; // Interpolated value of v.
                double Z_x; // Interpolation value of Z.

                // (3a) Compute geopotential height of current grid point and add Z component
                // of x-Z-plane tangent vector.
                    // PV : 0 t, 1 q, 2 u, 3 v, 4 sfcPhi, 5 sfc t_d, 6 sfc t, 7 sfc p

                    // Surface geopotential.
                    double geopotSurf_J = inputGrids.at(4)-&gt;getValue(0, j, i);

                    // Assuming temperature ranging from -45 to +60 degrees Celsius.

                    // Raphido: Only sfc T_v is used. Use mean value instead, also in geopotential height method.
                    // Surface virtual temperature.
                    double pSurf_hPa = inputGrids.at(7)-&gt;getValue(0, j, i) / 100.;
                    double virtualTempSurf_K = virtualTempFromDewPoint(inputGrids.at(5)-&gt;getValue(0, j, i),
                                                                       pSurf_hPa,
                                                                       inputGrids.at(6)-&gt;getValue(0, j, i));

                    // Surface layer geopotential thickness.
                    double dPhiSurfLayer = geopotThickness(virtualTempSurf_K,
                                                           inputGrids.at(0)-&gt;getPressure(inputGrids.at(0)-&gt;getNumLevels() - 1, j, i),
                                                           pSurf_hPa);

                    double geopotential000_J = 0.;

                    double virtualTemp_av_K = 0.;
                    double dPhi_J = 0.;

                    // For each point, loop from surface pressure level to current pressure level (k).
                    for (l = derivedGrid-&gt;getNumLevels() - 2; l &gt; k-1; l--)
                    {

                        // Compute average virtual temperature between current pressure level
                        // and lower level.
                        virtualTemp_av_K = ( virtualTempFromSpecHum(inputGrids.at(0)-&gt;getValue(l, j, i),
                                                                    inputGrids.at(1)-&gt;getValue(l, j, i))
                                           + virtualTempFromSpecHum(inputGrids.at(0)-&gt;getValue(l+1, j, i),
                                                                    inputGrids.at(1)-&gt;getValue(l+1, j, i)) ) / 2.;

                        // Compute the current thickness.
                        dPhi_J = geopotThickness(virtualTemp_av_K,
                                                 inputGrids.at(0)-&gt;getPressure(l, j, i),
                                                 inputGrids.at(0)-&gt;getPressure(l+1, j, i));

                        // Sum up the differentials.
                        geopotential000_J += dPhi_J;
                    }

                    // Add surface layer geopotential thickness
                    geopotential000_J += dPhiSurfLayer;

                    // Add surface geopotential
                    geopotential000_J += geopotSurf_J;

                    // Convert to geopotential height by dividing through g_0.
                    Z_x = geopotToZ(geopotential000_J);

                    // Add gradient contribution.
                    Z_x += dZ_x;

                // (3b) Loop through model levels to find adjacent grid points enclosing Z_x.
                // Starting from the lowermost level, compute geopotential height of adjacent grid points.
                for(n = derivedGrid-&gt;getNumLevels() - 1; n &gt; 1; n--)
                {

                        // Surface geopotential.
                        double geopotSurf_J = inputGrids.at(4)-&gt;getValue(0, j, i+1);

                        // Assuming temperature ranging from -45 to +60 degrees Celsius

                        // Surface virtual temperature.
                        double pSurf_hPa = inputGrids.at(7)-&gt;getValue(0, j, i+1) / 100.;
                        double virtualTempSurf_K = virtualTempFromDewPoint(inputGrids.at(5)-&gt;getValue(0, j, i+1),
                                                                          pSurf_hPa,
                                                                          inputGrids.at(6)-&gt;getValue(0, j, i+1));

                        // Surface layer geopotential thickness
                        double dPhiSurfLayer = geopotThickness(virtualTempSurf_K,
                                                              inputGrids.at(0)-&gt;getPressure(inputGrids.at(0)-&gt;getNumLevels() - 1, j, i+1),
                                                              pSurf_hPa);

                        // Lower neighbour.
                        double geopotential_J_lower = geopotSurf_J;
                        double Zx_bot = geopotToZ(geopotential_J_lower);

                        // Upper neighbour.
                        double geopotential_J_upper = geopotential_J_lower + dPhiSurfLayer;
                        double Zx_top = geopotToZ(geopotential_J_upper);

                        // Check whether adjacent grid points enclose Z_x.
                        // If they do, check whether surface is involved or not, then interpolate.
                        if((Zx_bot &lt; Z_x) &amp;&amp; (Z_x &lt; Zx_top))
                        {

                            // Surface level involved.
                            if(n == derivedGrid-&gt;getNumLevels() - 1)
                            {
                                double v_top = inputGrids.at(3)-&gt;getValue(inputGrids.at(0)-&gt;getNumLevels() - 1, j, i+1);
                                double v_bot = inputGrids.at(9)-&gt;getValue(0, j, i+1);

                                // Interpolated value of v.
                                vdx = v_bot + (v_top - v_bot) * (Z_x - Zx_bot) / (Zx_top - Zx_bot);

                                break;
                            }

                            // Surface level not involved.
                            else
                            {
                                double v_top = inputGrids.at(3)-&gt;getValue(n-1, j, i+1);
                                double v_bot = inputGrids.at(3)-&gt;getValue(n, j, i+1);

                                // Interpolated value of v.
                                vdx = v_bot + (v_top - v_bot) * (Z_x - Zx_bot) / (Zx_top - Zx_bot);

                                break;
                            }

                        }
                        // Z_x not enclosed, go one level up.
                        else
                        {
                            geopotential_J_lower = geopotential_J_upper;
                            Zx_bot = geopotToZ(geopotential_J_lower);

                            double virtualTemp_av_K = 0.;
                            double dPhi_J = 0.;



                            // Compute average virtual temperature between current pressure level
                            // and lower level.
                            virtualTemp_av_K = ( virtualTempFromSpecHum(inputGrids.at(0)-&gt;getValue(n, j, i+1),
                                                                        inputGrids.at(1)-&gt;getValue(n, j, i+1))
                                               + virtualTempFromSpecHum(inputGrids.at(0)-&gt;getValue(n-1, j, i+1),
                                                                        inputGrids.at(1)-&gt;getValue(n-1, j, i+1)) ) / 2.;

                            // Compute the current thickness.
                            dPhi_J = geopotThickness(virtualTemp_av_K,
                                                     inputGrids.at(0)-&gt;getPressure(n-1, j, i+1),
                                                     inputGrids.at(0)-&gt;getPressure(n, j, i+1));


                            geopotential_J_upper += dPhi_J;
                            Zx_top = geopotToZ(geopotential_J_upper);
                        }

                }

                // y-coordinate
                double udy; // Interpolated value of u.
                double Z_y; // Interpolation value of Z.

                // (3c) Compute geopotential height of current grid point and add Z component
                // of y-Z-plane tangent vector

                    // Geopotential at current grid point already known.
                    // Convert to geopotential height.
                    Z_y = geopotToZ(geopotential000_J);

                    // Add gradient contribution.
                    Z_y += dZ_y;

                // (3d) Loop through model levels to find adjacent grid points enclosing Z_y.
                // Starting from the lowermost level, compute geopotential height of adjacent grid points.
                for(n = derivedGrid-&gt;getNumLevels() - 1; n &gt; 1; n--)
                {

                        // Surface geopotential.
                        double geopotSurf_J = inputGrids.at(4)-&gt;getValue(0, j+1, i);

                        // Assuming temperature ranging from -45 to +60 degrees Celsius.

                        // Surface virtual temperature.
                        double pSurf_hPa = inputGrids.at(7)-&gt;getValue(0, j+1, i) / 100.;
                        double virtualTempSurf_K = virtualTempFromDewPoint(inputGrids.at(5)-&gt;getValue(0, j+1, i),
                                                                           pSurf_hPa,
                                                                           inputGrids.at(6)-&gt;getValue(0, j+1, i));

                        // Surface layer geopotential thickness
                        double dPhiSurfLayer = geopotThickness(virtualTempSurf_K,
                                                               inputGrids.at(0)-&gt;getPressure(inputGrids.at(0)-&gt;getNumLevels() - 1, j+1, i),
                                                               pSurf_hPa);

                        // Lower neighbour.
                        double geopotential_J_lower = geopotSurf_J;
                        double Zy_bot = geopotToZ(geopotential_J_lower);

                        // Upper neighbour.
                        double geopotential_J_upper = geopotential_J_lower + dPhiSurfLayer;
                        double Zy_top = geopotToZ(geopotential_J_upper);

                        // Check whether adjacent grid points enclose Z_x.
                        // If they do, check whether surface is involved or not, then interpolate.
                        if((Zy_bot &lt; Z_y) &amp;&amp; (Z_y &lt; Zy_top))
                        {

                            // Surface level involved.
                            if(n == derivedGrid-&gt;getNumLevels() - 1)
                            {
                                double u_top = inputGrids.at(2)-&gt;getValue(inputGrids.at(0)-&gt;getNumLevels() - 1, j+1, i);
                                double u_bot = inputGrids.at(8)-&gt;getValue(0, j+1, i);

                                // Interpolated value of v.
                                udy = u_bot + (u_top - u_bot) * (Z_y - Zy_bot) / (Zy_top - Zy_bot);

                                break;
                            }

                            // Surface level not involved.
                            else
                            {
                                double u_top = inputGrids.at(2)-&gt;getValue(n-1, j+1, i);
                                double u_bot = inputGrids.at(2)-&gt;getValue(n, j+1, i);

                                // Interpolated value of u.
                                udy = u_bot + (u_top - u_bot) * (Z_y - Zy_bot) / (Zy_top - Zy_bot);

                                break;
                            }

                        }
                        // Z_y not enclosed, go one level up.
                        else
                        {
                            geopotential_J_lower = geopotential_J_upper;
                            Zy_bot = geopotToZ(geopotential_J_lower);

                            double virtualTemp_av_K = 0.;
                            double dPhi_J = 0.;



                            // Compute average virtual temperature between current pressure level
                            // and lower level.
                            virtualTemp_av_K = ( virtualTempFromSpecHum(inputGrids.at(0)-&gt;getValue(n, j+1, i),
                                                                        inputGrids.at(1)-&gt;getValue(n, j+1, i))
                                               + virtualTempFromSpecHum(inputGrids.at(0)-&gt;getValue(n-1, j+1, i),
                                                                        inputGrids.at(1)-&gt;getValue(n-1, j+1, i)) ) / 2.;

                            // Compute the current thickness.
                            dPhi_J = geopotThickness(virtualTemp_av_K,
                                                     inputGrids.at(0)-&gt;getPressure(n-1, j+1, i),
                                                     inputGrids.at(0)-&gt;getPressure(n, j+1, i));


                            geopotential_J_upper += dPhi_J;
                            Zy_top = geopotToZ(geopotential_J_upper);
                        }

                }

        // OLD: Compute geometrical length of tangent vectors.
        // NEW: Write new method using the already computed dZ's.
        double dx_Z = vecLengthZ(dx, 0, dZ_x);

        double dy_Z = vecLengthZ(0, dy, dZ_y);

        // Compute wind differences
        double du = udy - inputGrids.at(2)-&gt;getValue(k, j, i);
        double dv = vdx - inputGrids.at(3)-&gt;getValue(k, j, i);

        // Compute isentropic relative vorticity.
        double vo_rel = (du / dy_Z) - (dv / dx_Z);

        // Compute static stability.
        double dTheta = theta100 - theta000;
        double dp = 100. * (p100 - p000);

        // Compute planetary vorticity.
        // vo_plan = 2. * w_0 * sin(2. * M_PI * lat / 360.);
        double w_0 = MetConstants::EARTH_ROTATION_RATE;
        double lat000 = inputGrids.at(0)-&gt;getNorthInterfaceLat(j);
        double vo_plan = 2. * w_0 * sin(2. * M_PI * lat000 / 360.);

        // Compute PV[pvu]
        double factor = MetConstants::PVU_FACTOR;
        double g = MetConstants::GRAVITY_ACCELERATION;
        double PV_pvu = factor * g * (vo_rel + vo_plan) * (dTheta / dp);

        derivedGrid-&gt;setValue(k, j, i, PV_pvu);
    }
    }

}
</code></pre>

CC BY-SA 3.0
2017-01-03T09:11:36.930Z
@gansub I edited my post and added the formula, without writing out the conversion factor! I also had a look at this Python package: as far as I can see, the method to compute isentropic PV is not explicitly shown and data are already given on isentropic levels! But thanks, I will try and see if this package can help me!
RaphiK
CC BY-SA 3.0
2017-01-03T10:33:14.620Z
@gansub thanks for your effort! I added some more specific information to my post in order to give a better description of my problem!
RaphiK
CC BY-SA 3.0
2017-01-04T05:39:37.160Z
Expanding on the link from @gansub: They mention both that the fields have to be interpolated to isentropic surfaces and that they "spectrally" compute relative vorticity. Relative vorticity is super-sensitive to small-scale changes; you'll probably have to smooth the fields to get reasonable PV values. idk if your environment already has smoothing routines, but I'm sure there are a ton in C/C++.
Jareth Holt
<meteorology><numerical-modelling><data-analysis><vorticity>
Computing Isentropic Potential Vorticity from Hybrid Sigma Pressure levels
RaphiK
# Computing Isentropic Potential Vorticity from Hybrid Sigma Pressure levels

(Asked by: username_0 on 2017-01-02)

My task is to compute isentropic potential vorticity on a grid.
All data (ECMWF ens) are given on a hybrid sigma pressure grid.
Since I really cannot find reference, I assume (!) that relative vorticity (vo, id=138) is evaluated on model levels and therefore is useless for me.
My approach (for each grid point):

Compute theta-gradient, using the three model dimensions, converting vertical and horizontal distances to meters.
Define two vectors using three constraints each:

Scalar product with theta gradient vanishes.
First vector has no meridional, second one no zonal component.
Distance to adjacent horizontal grid point defines magnitude of 
meridional/zonal component.

These two vectors now span a plane tangential to the isentropic surface at the current grid point. They point to positions only vertically displaced from the horizontal neighbours of the current grid point.
Values of u and v are interpolated - again, using meters (geopotential height) - to the locations given by the tangent vectors.
The geometrical length of these tangent vectors is determined and used together with the interpolated values of u and v to compute isentropic relative vorticity.
Static stability is computed using values of pressure and theta at current point and vertical neighbour.
Using gravitational acceleration and the conversion factor of 10^-6, I use the equation given in the famous book by Wallace & Hobbs for Potential Vorticity: $g(f + \eta_{\theta})\frac{\partial \theta}{\partial p}$

My problem:
The result does not look like what I expect.
Values range from around -30 to +50 pvu throughout the atmosphere, while the 1.5pvu surface clearly does not depict the tropopause. It looks more like an accident. I am using "omp" to increase computation speed, but things look the same if I don't use it, too, so it cannot be index confusion.
My question:
What could be wrong?
Can you spot systematic errors? Maybe I got blind with time.
Do you know reference for how relative vorticity is evaluated in the ECMWF model? (I tried, really...)
To be more specific:
My method is to be implemented in C++ into an already existing environment to visualize NWP data in 3D in realtime.
The basic datafields - given in hybrid sigma coordinates - are read, processed, written into a derived data field, interpolated to and rendered on pressure levels as triggered by user interaction.
The program should possess the processing methods on its own, which is why I cannot use Python packages, for example.
Here is the source code of my approach.  If something is unclear, please ask!
Maybe somebody will spot an error.
 // ISENTROPIC POTENTIAL VORTICITY
// Input grids : 0 t, 1 q, 2 u, 3 v, 4 sfcPhi, 5 sfc t_d, 6 sfc t, 7 sfc p, 8 sfc u, 9 sfc v
else if ( (derivedVarName == "Isentropic PV (an)")
     || (derivedVarName == "Isentropic PV (fc)")
     || (derivedVarName == "Isentropic PV (ens)") )
{
    unsigned int k, j, i;
    int n, l;
#pragma omp parallel
        {
        #pragma omp for private(k, j, i, n, l)
    for (k = 1; k < derivedGrid->getNumLevels() - 1; k++)
        for (j = 0; j < derivedGrid->getNumLats() - 1; j++)
            for (i = 0; i < derivedGrid->getNumLons() - 1; i++)
    {

        // Current grid point.
        double T000 = inputGrids.at(0)->getValue(k, j, i);
        double p000 = inputGrids.at(0)->getPressure(k, j, i);
        double q000 = inputGrids.at(1)->getValue(k, j, i);
        double theta000 = potTemp(T000, p000);

        // Upper neighbour.
        double T100 = inputGrids.at(0)->getValue(k-1, j, i);
        double p100 = inputGrids.at(0)->getPressure(k-1, j, i);
        double q100 = inputGrids.at(1)->getValue(k-1, j, i);
        double theta100 = potTemp(T100, p100);

        // x-neighbour.
        double theta001 = potTemp(inputGrids.at(0)->getValue(k, j, i+1),
                                  inputGrids.at(0)->getPressure(k, j, i+1));

        // y-neighbour.
        double theta010 = potTemp(inputGrids.at(0)->getValue(k, j+1, i),
                                  inputGrids.at(0)->getPressure(k, j+1, i));

        // (1) Compute theta gradient at current point, using geopotential thickness dZ [m] for the vertical,
        // great circle distance [m] for the horizontal distances.

            // Compute geopotential thickness between current point and upper neighbour.
            double T_v = (virtualTempFromSpecHum(T000, q000) + virtualTempFromSpecHum(T100, q100)) / 2.;
            double dPhi = geopotThickness(T_v, p100, p000);
            double dZ = geopotToZ(dPhi);

            // Compute gradient.
            double* grad = gradient(theta000,
                                    inputGrids.at(0)->getEastInterfaceLon(i),
                                    inputGrids.at(0)->getNorthInterfaceLat(j), // Current grid point.

                                    theta001,
                                    inputGrids.at(0)->getEastInterfaceLon(i+1), // xlon-neighbour.

                                    theta010,
                                    inputGrids.at(0)->getNorthInterfaceLat(j+1), // ylat-neighbour.

                                    theta100,
                                    dZ // Z-neighbour.
                                    );

            // Gradient components.
            double xgrad = grad[2];
            double ygrad = grad[1];
            double Zgrad = grad[0];

        // (2) Define tangent vectors and compute their Z-components.
        // Definition:
        // 1) Scalar product of tangent vectors with gradient vanishes.
        // 2) One tangent vector in y-z-plane, the other in x-z-plane.
        // 3) Length of horizontal components equals spacing between grid points.
        double radius_m = MetConstants::EARTH_RADIUS_km * 1000.;

            // Tangent vector in y-z-plane.
            double dx = gcDistance_deg(inputGrids.at(0)->getEastInterfaceLon(i),
                                       inputGrids.at(0)->getNorthInterfaceLat(j),
                                       inputGrids.at(0)->getEastInterfaceLon(i+1),
                                       inputGrids.at(0)->getNorthInterfaceLat(j),
                                       radius_m);

            // Z-component.
            double dZ_x = -dx * xgrad / Zgrad;

            // Tangent vector in x-z-plane.
            double dy = gcDistance_deg(inputGrids.at(0)->getEastInterfaceLon(i),
                                       inputGrids.at(0)->getNorthInterfaceLat(j),
                                       inputGrids.at(0)->getEastInterfaceLon(i),
                                       inputGrids.at(0)->getNorthInterfaceLat(j+1),
                                       radius_m);

            // Z-component.
            double dZ_y = -dy * ygrad / Zgrad;

        // Tangent vectors now point from current grid point to locations only
        // vertically displaced from x- or y-neighbour, respectively.

        // (3) Find levels whose Z values enclose geopotential height values and
        // interpolate u and v to these values.
        // The geopotential height of the current grid point is computed in the
        // x-coordinate section and reused in the y-coordinate section.

                // x-coordinate
                double vdx; // Interpolated value of v.
                double Z_x; // Interpolation value of Z.

                // (3a) Compute geopotential height of current grid point and add Z component
                // of x-Z-plane tangent vector.
                    // PV : 0 t, 1 q, 2 u, 3 v, 4 sfcPhi, 5 sfc t_d, 6 sfc t, 7 sfc p

                    // Surface geopotential.
                    double geopotSurf_J = inputGrids.at(4)->getValue(0, j, i);

                    // Assuming temperature ranging from -45 to +60 degrees Celsius.

                    // Raphido: Only sfc T_v is used. Use mean value instead, also in geopotential height method.
                    // Surface virtual temperature.
                    double pSurf_hPa = inputGrids.at(7)->getValue(0, j, i) / 100.;
                    double virtualTempSurf_K = virtualTempFromDewPoint(inputGrids.at(5)->getValue(0, j, i),
                                                                       pSurf_hPa,
                                                                       inputGrids.at(6)->getValue(0, j, i));

                    // Surface layer geopotential thickness.
                    double dPhiSurfLayer = geopotThickness(virtualTempSurf_K,
                                                           inputGrids.at(0)->getPressure(inputGrids.at(0)->getNumLevels() - 1, j, i),
                                                           pSurf_hPa);

                    double geopotential000_J = 0.;

                    double virtualTemp_av_K = 0.;
                    double dPhi_J = 0.;

                    // For each point, loop from surface pressure level to current pressure level (k).
                    for (l = derivedGrid->getNumLevels() - 2; l > k-1; l--)
                    {

                        // Compute average virtual temperature between current pressure level
                        // and lower level.
                        virtualTemp_av_K = ( virtualTempFromSpecHum(inputGrids.at(0)->getValue(l, j, i),
                                                                    inputGrids.at(1)->getValue(l, j, i))
                                           + virtualTempFromSpecHum(inputGrids.at(0)->getValue(l+1, j, i),
                                                                    inputGrids.at(1)->getValue(l+1, j, i)) ) / 2.;

                        // Compute the current thickness.
                        dPhi_J = geopotThickness(virtualTemp_av_K,
                                                 inputGrids.at(0)->getPressure(l, j, i),
                                                 inputGrids.at(0)->getPressure(l+1, j, i));

                        // Sum up the differentials.
                        geopotential000_J += dPhi_J;
                    }

                    // Add surface layer geopotential thickness
                    geopotential000_J += dPhiSurfLayer;

                    // Add surface geopotential
                    geopotential000_J += geopotSurf_J;

                    // Convert to geopotential height by dividing through g_0.
                    Z_x = geopotToZ(geopotential000_J);

                    // Add gradient contribution.
                    Z_x += dZ_x;

                // (3b) Loop through model levels to find adjacent grid points enclosing Z_x.
                // Starting from the lowermost level, compute geopotential height of adjacent grid points.
                for(n = derivedGrid->getNumLevels() - 1; n > 1; n--)
                {

                        // Surface geopotential.
                        double geopotSurf_J = inputGrids.at(4)->getValue(0, j, i+1);

                        // Assuming temperature ranging from -45 to +60 degrees Celsius

                        // Surface virtual temperature.
                        double pSurf_hPa = inputGrids.at(7)->getValue(0, j, i+1) / 100.;
                        double virtualTempSurf_K = virtualTempFromDewPoint(inputGrids.at(5)->getValue(0, j, i+1),
                                                                          pSurf_hPa,
                                                                          inputGrids.at(6)->getValue(0, j, i+1));

                        // Surface layer geopotential thickness
                        double dPhiSurfLayer = geopotThickness(virtualTempSurf_K,
                                                              inputGrids.at(0)->getPressure(inputGrids.at(0)->getNumLevels() - 1, j, i+1),
                                                              pSurf_hPa);

                        // Lower neighbour.
                        double geopotential_J_lower = geopotSurf_J;
                        double Zx_bot = geopotToZ(geopotential_J_lower);

                        // Upper neighbour.
                        double geopotential_J_upper = geopotential_J_lower + dPhiSurfLayer;
                        double Zx_top = geopotToZ(geopotential_J_upper);

                        // Check whether adjacent grid points enclose Z_x.
                        // If they do, check whether surface is involved or not, then interpolate.
                        if((Zx_bot < Z_x) && (Z_x < Zx_top))
                        {

                            // Surface level involved.
                            if(n == derivedGrid->getNumLevels() - 1)
                            {
                                double v_top = inputGrids.at(3)->getValue(inputGrids.at(0)->getNumLevels() - 1, j, i+1);
                                double v_bot = inputGrids.at(9)->getValue(0, j, i+1);

                                // Interpolated value of v.
                                vdx = v_bot + (v_top - v_bot) * (Z_x - Zx_bot) / (Zx_top - Zx_bot);

                                break;
                            }

                            // Surface level not involved.
                            else
                            {
                                double v_top = inputGrids.at(3)->getValue(n-1, j, i+1);
                                double v_bot = inputGrids.at(3)->getValue(n, j, i+1);

                                // Interpolated value of v.
                                vdx = v_bot + (v_top - v_bot) * (Z_x - Zx_bot) / (Zx_top - Zx_bot);

                                break;
                            }

                        }
                        // Z_x not enclosed, go one level up.
                        else
                        {
                            geopotential_J_lower = geopotential_J_upper;
                            Zx_bot = geopotToZ(geopotential_J_lower);

                            double virtualTemp_av_K = 0.;
                            double dPhi_J = 0.;



                            // Compute average virtual temperature between current pressure level
                            // and lower level.
                            virtualTemp_av_K = ( virtualTempFromSpecHum(inputGrids.at(0)->getValue(n, j, i+1),
                                                                        inputGrids.at(1)->getValue(n, j, i+1))
                                               + virtualTempFromSpecHum(inputGrids.at(0)->getValue(n-1, j, i+1),
                                                                        inputGrids.at(1)->getValue(n-1, j, i+1)) ) / 2.;

                            // Compute the current thickness.
                            dPhi_J = geopotThickness(virtualTemp_av_K,
                                                     inputGrids.at(0)->getPressure(n-1, j, i+1),
                                                     inputGrids.at(0)->getPressure(n, j, i+1));


                            geopotential_J_upper += dPhi_J;
                            Zx_top = geopotToZ(geopotential_J_upper);
                        }

                }

                // y-coordinate
                double udy; // Interpolated value of u.
                double Z_y; // Interpolation value of Z.

                // (3c) Compute geopotential height of current grid point and add Z component
                // of y-Z-plane tangent vector

                    // Geopotential at current grid point already known.
                    // Convert to geopotential height.
                    Z_y = geopotToZ(geopotential000_J);

                    // Add gradient contribution.
                    Z_y += dZ_y;

                // (3d) Loop through model levels to find adjacent grid points enclosing Z_y.
                // Starting from the lowermost level, compute geopotential height of adjacent grid points.
                for(n = derivedGrid->getNumLevels() - 1; n > 1; n--)
                {

                        // Surface geopotential.
                        double geopotSurf_J = inputGrids.at(4)->getValue(0, j+1, i);

                        // Assuming temperature ranging from -45 to +60 degrees Celsius.

                        // Surface virtual temperature.
                        double pSurf_hPa = inputGrids.at(7)->getValue(0, j+1, i) / 100.;
                        double virtualTempSurf_K = virtualTempFromDewPoint(inputGrids.at(5)->getValue(0, j+1, i),
                                                                           pSurf_hPa,
                                                                           inputGrids.at(6)->getValue(0, j+1, i));

                        // Surface layer geopotential thickness
                        double dPhiSurfLayer = geopotThickness(virtualTempSurf_K,
                                                               inputGrids.at(0)->getPressure(inputGrids.at(0)->getNumLevels() - 1, j+1, i),
                                                               pSurf_hPa);

                        // Lower neighbour.
                        double geopotential_J_lower = geopotSurf_J;
                        double Zy_bot = geopotToZ(geopotential_J_lower);

                        // Upper neighbour.
                        double geopotential_J_upper = geopotential_J_lower + dPhiSurfLayer;
                        double Zy_top = geopotToZ(geopotential_J_upper);

                        // Check whether adjacent grid points enclose Z_x.
                        // If they do, check whether surface is involved or not, then interpolate.
                        if((Zy_bot < Z_y) && (Z_y < Zy_top))
                        {

                            // Surface level involved.
                            if(n == derivedGrid->getNumLevels() - 1)
                            {
                                double u_top = inputGrids.at(2)->getValue(inputGrids.at(0)->getNumLevels() - 1, j+1, i);
                                double u_bot = inputGrids.at(8)->getValue(0, j+1, i);

                                // Interpolated value of v.
                                udy = u_bot + (u_top - u_bot) * (Z_y - Zy_bot) / (Zy_top - Zy_bot);

                                break;
                            }

                            // Surface level not involved.
                            else
                            {
                                double u_top = inputGrids.at(2)->getValue(n-1, j+1, i);
                                double u_bot = inputGrids.at(2)->getValue(n, j+1, i);

                                // Interpolated value of u.
                                udy = u_bot + (u_top - u_bot) * (Z_y - Zy_bot) / (Zy_top - Zy_bot);

                                break;
                            }

                        }
                        // Z_y not enclosed, go one level up.
                        else
                        {
                            geopotential_J_lower = geopotential_J_upper;
                            Zy_bot = geopotToZ(geopotential_J_lower);

                            double virtualTemp_av_K = 0.;
                            double dPhi_J = 0.;



                            // Compute average virtual temperature between current pressure level
                            // and lower level.
                            virtualTemp_av_K = ( virtualTempFromSpecHum(inputGrids.at(0)->getValue(n, j+1, i),
                                                                        inputGrids.at(1)->getValue(n, j+1, i))
                                               + virtualTempFromSpecHum(inputGrids.at(0)->getValue(n-1, j+1, i),
                                                                        inputGrids.at(1)->getValue(n-1, j+1, i)) ) / 2.;

                            // Compute the current thickness.
                            dPhi_J = geopotThickness(virtualTemp_av_K,
                                                     inputGrids.at(0)->getPressure(n-1, j+1, i),
                                                     inputGrids.at(0)->getPressure(n, j+1, i));


                            geopotential_J_upper += dPhi_J;
                            Zy_top = geopotToZ(geopotential_J_upper);
                        }

                }

        // OLD: Compute geometrical length of tangent vectors.
        // NEW: Write new method using the already computed dZ's.
        double dx_Z = vecLengthZ(dx, 0, dZ_x);

        double dy_Z = vecLengthZ(0, dy, dZ_y);

        // Compute wind differences
        double du = udy - inputGrids.at(2)->getValue(k, j, i);
        double dv = vdx - inputGrids.at(3)->getValue(k, j, i);

        // Compute isentropic relative vorticity.
        double vo_rel = (du / dy_Z) - (dv / dx_Z);

        // Compute static stability.
        double dTheta = theta100 - theta000;
        double dp = 100. * (p100 - p000);

        // Compute planetary vorticity.
        // vo_plan = 2. * w_0 * sin(2. * M_PI * lat / 360.);
        double w_0 = MetConstants::EARTH_ROTATION_RATE;
        double lat000 = inputGrids.at(0)->getNorthInterfaceLat(j);
        double vo_plan = 2. * w_0 * sin(2. * M_PI * lat000 / 360.);

        // Compute PV[pvu]
        double factor = MetConstants::PVU_FACTOR;
        double g = MetConstants::GRAVITY_ACCELERATION;
        double PV_pvu = factor * g * (vo_rel + vo_plan) * (dTheta / dp);

        derivedGrid->setValue(k, j, i, PV_pvu);
    }
    }

}

**Tags:** meteorology, numerical-modelling, data-analysis, vorticity

> **Comments:**
> **username_0**: @gansub I edited my post and added the formula, without writing out the conversion factor! I also had a look at this Python package: as far as I can see, the method to compute isentropic PV is not explicitly shown and data are already given on isentropic levels! But thanks, I will try and see if this package can help me!  
> **username_0**: @gansub thanks for your effort! I added some more specific information to my post in order to give a better description of my problem!  
> **username_1**: Expanding on the link from @gansub: They mention both that the fields have to be interpolated to isentropic surfaces and that they "spectrally" compute relative vorticity. Relative vorticity is super-sensitive to small-scale changes; you'll probably have to smooth the fields to get reasonable PV values. idk if your environment already has smoothing routines, but I'm sure there are a ton in C/C++.  

---

## 0 answers

