<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):

1. Compute theta-gradient, using the three model dimensions, converting vertical and horizontal distances to meters.
2. 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.
3. 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.
4. Values of u and v are interpolated - again, using meters (geopotential height) - to the locations given by the tangent vectors.
5. 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.
6. Static stability is computed using values of pressure and theta at current point and vertical neighbour.
7. 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

