Computing Isentropic Potential Vorticity from Hybrid Sigma Pressure levels

Question

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);
    }
    }

}

Answer 1

Based on an answer from ECMWF support which OP can confirm to my best understanding relative vorticity in the ECMWF model is not calculated using grid points and finite differences (centered and forward and backward). Instead it is calculated using spectral approaches in meteorology. There is a package in Fortran called spherepack and a python wrapper as well pyspharm that one can use to convert grid points to spectral approach but the data has to be global. Output of regional mesoscale models will not work.

So if you are using data from a GCM that is using grid points I would suggest you convert to spectral approach first, calculate the vorticity(HINT : use windspharm for this purpose). As Jareth Holt mentions in the comments you will need to smoothen the relative vorticity and here is a link that shows you how to do that - Smoothen Relative Vorticity. If values of absolute vorticity appear correct(and by that I mean by and large positive in NH and negative in SH) then proceed to calculate the isentropic PV by interpolating from model levels to isentropic levels using Newton Raphson interpolation.

If you are going to stick with grid points and finite differences the formulae for relative vorticity must include map scale factors assuming data is available as a function of latitude and longitude and you need some formula for calculating vorticity at the poles. Your code as it stands does not include map scale factors in the relative vorticity calculation.

There are few other issues that you will need to take into consideration. When these isentropic surfaces intersect the earth's surface the values for temperature and surface pressure will need to be extrapolated. The necessary ECMWF documentation for that is available on a google search for the same. Finally all super adiabatic layers must be eliminated prior to isentropic interpolation.

Instead of showing us the entire code please verify if the several "little" calculations that go into calculating PV(relative vorticity, absolute vorticity, calculation of static stability, isentropic interpolation)are correct and then the calculations sorts itself out.