等压面线性插值

1. 背景与目的

最近疯狂处理ICON模式数据，包括垂直层的数据，对于他的模式设置和输出数据格式有了更深的了解，下面简单记录一下。

在大气模式输出（如AMIP）或再分析资料中，动力与热力变量通常定义在 模式垂直坐标（hybrid σ–p、η 或 model level）上，而非固定的等压面。  这里特别的，以ICON风暴解析的AMIP输出为例，请看下图：

• 低层密，高层稀疏
• 从上到下设置垂直层数编号

在 ICON 中，选择的是基于高度的坐标系，该坐标系遵循地形，因此顶部和底部三角形面相对于球体上的切平面倾斜。由于随着距下边界距离的增加，模型水平逐渐变为恒定高度的水平，因此网格盒的顶部和底部三角形面也彼此稍微倾斜。每个网格框的确切高度取决于地球上的地理位置。顶面和底面称为垂直网格的半层，盒子的中心称为垂直网格的全层，参见图 3.2 的说明。请注意，全级和半级的编号是自上而下的，从上半级和全级的 k = 1 开始。垂直方向采用洛伦兹型交错，这意味着水平速度、位温和密度定义为全水平，而垂直速度定义为半水平。而且，模式默认输出的坐标为模式层数

然而，大多数诊断分析（例如散度、波动结构、垂直剖面比较）均基于标准气压层进行。通常两种做法是：

• 1、将模式层转换为具体的高度，这在只分析海洋数据时非常方便。因为海洋上，每个点的高度是完全相同的，直接将zg分配为新的垂直维度即可。

mse_cntl_new = mse_cntl.assign_coords(level = zg_profile.values)

• 2、根据气压数据，将其插值到对应的压力层上: 插值到

下面讲怎么实现第二种方法

2. 插值基本假设

采用 垂直线性插值（linear interpolation in pressure coordinate），并基于以下假设：

1. 压力在垂直方向单调变化（通常自上而下递增）；
2. 在相邻两个模式层之间，变量随压力变化近似线性；
3. 插值仅沿垂直方向进行

3. 数学原理

设目标等压面为，位于两个相邻模式层之间：

对应的变量值为：

则目标层上的插值结果为：

其中线性权重  定义为：

该表达式等价于压力坐标下的一阶插值近似：根据目标点在两个已知点之间的“距离比例”，对两个端点的变量值进行加权平均。

4. 代码实现

4.1 垂直维度自动识别

通过匹配维度名称中的关键字（如 "level"），自动识别变量与压力场的垂直维度：

4.2 确定插值区间

对每个时间与水平格点，搜索满足：

的第一个模式层索引 ( k )，并取：

• 上层：
• 下层：

这两个层即构成插值所需的垂直区间。

4.3 构造插值权重

利用上下层压力值计算线性权重：

该权重在空间和时间上完全向量化计算

4.4 变量插值

对物理变量本身进行线性组合：

4.5 构造标准等压面维度

插值结果新增 plev 维度，并附加 CF 规范元数据：

• standard_name = air_pressure
• units = hPa
• axis = Z
• positive = down

4.6 多压力层插值策略

针对多个目标气压层：

1. 对每个压力层独立构建插值计算图；
2. 使用 xr.concat 一次性合并为多层等压面结果；
3. 在内存受限环境下，可选择逐层计算并拼接。

defhacky_plvl_interpolation_plus(data_var, data_pfull, plvl_target, unit='hPa', 
                             label='pressure level', height_name='level'):
"""
    Interpolate variable to a specific pressure level.

    Parameters
    ----------
    data_var : xr.DataArray
        Variable to interpolate
    data_pfull : xr.DataArray
        Pressure field (in same units as plvl_target)
    plvl_target : float
        Target pressure level
    unit : str
        Unit of pressure level
    label : str
        Long name for pressure level
    height_name : str
        Substring to identify vertical dimension

    Returns
    -------
    xr.DataArray
        Interpolated variable at target pressure level
    """

    var_height_dim = next((dim for dim in data_var.dims if height_name in dim), None)
    pfull_height_dim = next((dim for dim in data_pfull.dims if height_name in dim), None)

if var_height_dim isNoneor pfull_height_dim isNone:
raise ValueError(f"Could not find dimension containing '{height_name}'")

# Find levels for interpolation 
    level_above = (data_pfull > plvl_target).argmax(dim=pfull_height_dim).compute()
    level_below = level_above - 1

# Get values at surrounding levels
    value_above = data_pfull.isel({pfull_height_dim: level_above})
    value_below = data_pfull.isel({pfull_height_dim: level_below})

# Linear interpolation weight
    f = (plvl_target - value_below) / (value_above - value_below)

# Interpolate variable
    data_interpolated = (
        (1 - f) * data_var.isel({var_height_dim: level_below}) + 
        f * data_var.isel({var_height_dim: level_above})
    )

# Clean up and add dimension
    data_interpolated = (
        data_interpolated
        .drop_vars(pfull_height_dim, errors='ignore')
        .expand_dims(dim={"plev": [plvl_target]}, axis=0)
    )

# Add metadata
    data_interpolated['plev'].attrs = {
'standard_name': 'air_pressure',
'long_name': label,
'units': unit,
'axis': 'Z',
'positive': 'down'
    }

return data_interpolated


definterpolate_to_pressure_levels(data_var, data_pfull, plvl_list, 
                                   pressure_unit_conversion=100, 
                                   height_name='level'):
"""
    Interpolate variable to multiple pressure levels efficiently.

    Parameters
    ----------
    data_var : xr.DataArray
        Variable to interpolate
    data_pfull : xr.DataArray
        Pressure field
    plvl_list : list of float
        Target pressure levels in hPa
    pressure_unit_conversion : float
        Factor to convert data_pfull to hPa (default: 100 for Pa->hPa)
    height_name : str
        Substring to identify vertical dimension

    Returns
    -------
    xr.DataArray
        Interpolated variable at all target pressure levels
    """
# Convert pressure to hPa if needed
if pressure_unit_conversion != 1:
        data_pfull = data_pfull / pressure_unit_conversion

# Build list of interpolations (lazy operations)
    interpolated_levels = [
        hacky_plvl_interpolation_plus(data_var, data_pfull, plvl, height_name=height_name)
for plvl in plvl_list
    ]

# Concatenate all at once (much faster than iterative concat)
    result = xr.concat(interpolated_levels, dim='plev')

return result

调用函数：

plvl_list = [5e1, 1e2, 2e2, 2.5e2, 3e2, 3.5e2, 4e2, 5e2, 6e2, 7e2, 8e2, 8.5e2, 9e2, 1e3]# in hPa ... pfull is in Pa
plvl_list
results = []
level_results = []
for i, plvl in enumerate(plvl_list):
    print(f"  Level {i+1}/{len(plvl_list)}: {plvl} hPa", flush=True)
    level_data = hacky_plvl_interpolation_plus(
        div_data['div'], 
        pre_data['phalf'] / 100, 
        plvl, 
        height_name='level'
    )
    level_results.append(level_data)
result = xr.concat(level_results, dim='plev')

一个例子，使用wa来进行插值，维度都是levelxtimexlatxlon,别问为什么用的是pfull，他输出的对应的就是pfull....

p_list = [ 5e2,8.5e2]# in hPa ... pfull is in Pa
p_list
results = []
level_results = []
for i, plvl in enumerate(p_list):
    print(f"  Level {i+1}/{len(p_list)}: {plvl} hPa", flush=True)
    level_data = hacky_plvl_interpolation_plus(
        wa_cntl['wa'], 
        p_full['pfull'] / 100, 
        plvl, 
        height_name='level'
    )
    level_results.append(level_data)
result = xr.concat(level_results, dim='plev')