CDO

hi Pablo!

this is not about spline vs. linear or something. this is an interpolation from one vertical coordinate system (hybrid sigma pressure) to another one (pure pressure levels). yo can google this term to find out more.

In general the hybrid coordinate use some kind of surface pressure, geopotential higth, some coefficient (a and b) and a formula. ERA-5 should have a proper docu about it's vertical axis. Chere here for an old zoo of vertical coordinate system used on weather and climate.

hth
ralf

Hi Ralf,
thanks a lot for the quick answer! I understand, that the function, depending on the surface pressure, calculates the pressure levels corresponding to the hybrid sigma pressure levels and that this is no interpolation of any type.

Now, I want to have the data on consistent pressure levels, corresponding to the pressure levels of the hybrid sigma levels with a surface pressure of 1013 hPa. If I now pass this array of pressure levels to the ml2pl function it first calculates the pressures at which the model level values are stored (e.g. the pressure levels with a surface pressure of 990 hPa), and then has to interpolate somehow to the array of pressure levels I want the data to be stored in.

My question is, how this last step of interpolation between the pressure levels the model data is stored on that specific point (calculated by e.g. p(i,j,k,n) = = ap(k) + b(k)*ps(n,j,i)) and those I pass to the function. I attached a short Python-Script which shows what I would do with the Numpy function. I hope my question is clearer now.
Thanks again!

Pablo

ah, now I see. sorry for my ignorance. the final step is a linear interpolation. the corresponding code is in src/vertical_interp.cc, line 290 of the latest release (see download area):

template <typename T>
static inline double
vertical_interp_X_kernel(const T *restrict gt, const T *restrict hyb_press, long nl, double pres, long ngp, long nhlev,                                        
                         double missval)
{
  const long nh = nl + ngp;
  return (nl < 0) ? missval
                  : ((nh >= ngp * nhlev) ? gt[nl]
                                         : gt[nl] + (pres - hyb_press[nl]) * (gt[nh] - gt[nl]) / (hyb_press[nh] - hyb_press[nl]));
}

Thanks!