How to vertically integrate over a region with different integral bounds at every lat lon point?

I am trying to calculate the vertical integral of the horizontal advection of Moist Static Energy uisng ERA5 data over a particular region. I am using xarray hence I know how to interpolate and integrate:

ds6 = ds5.sel(year=2000)
m = Cp*ds6.t + Lv*ds6.q + ds6.z
hadv = mpcalc.advection(m, u=ds6.u, v=ds6.v)
inter = hadv.interp(level=xs, method="cubic")
inter.integrate('level')

The above integrates from 1000 mb to 100 mb. I have surface pressure data for each coordinate. How do I do the integral at each coordinate only from the surface pressure to 100 mb instead?

Here is what I have tried:

Ps = ds10.sp.sel(time='2021-12-01')
selected = inter.where(inter.level < Ps) #selecting data points above Ps only
selected.dropna(dim='level') #removing Nan so that we can integrate

This does not work because it removes all data points below the maximum Ps of any data point, and spatial variability because of Ps is gone. What do I do?

Solution

Say you have the following dataset ds:

<xarray.Dataset>
Dimensions:  (lon: 10, lat: 20, level: 50)
Coordinates:
  * lon      (lon) int64 0 1 2 3 4 5 6 7 8 9
  * lat      (lat) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
  * level    (level) float64 1.1e+03 1.078e+03 1.056e+03 ... 54.49 32.24 10.0
Data variables:
    sfc_p    (lon, lat) int64 1031 980 1047 1028 988 ... 990 978 1026 1055 1003
    adv      (lon, lat, level) float64 -0.9586 -0.04061 1.627 ... -1.318 0.8371

Then you can apply a function to every individual vertical profile using groupby. Check this out:

all_dims_except_level = [d for d in ds.dims if (d != "level")]
adv_gridpoint = ds.stack(gridpoint=all_dims_except_level)


def integral_from_sfc_p_to_p(data, target_pressure):
    sfc_p_at_gridpoint = data.sfc_p.values.item()

    levels_between_sfc_p_and_upper_limit = data.level.where(
        (data.level > target_pressure) & (data.level < data.sfc_p), drop=True
    ).values

    target_levels = np.concatenate(
        [[sfc_p_at_gridpoint], levels_between_sfc_p_and_upper_limit, [target_pressure]]
    )

    data_interp_to_sfc_p = data.interp(level=target_levels)
    data_integrated = data_interp_to_sfc_p.integrate("level")
    return data_integrated


print(
    adv_gridpoint.groupby("gridpoint")
    .apply(integral_from_sfc_p_to_p, target_pressure=100)
    .unstack()
)

yields:

<xarray.Dataset>
Dimensions:  (lon: 10, lat: 20)
Coordinates:
  * lon      (lon) int64 0 1 2 3 4 5 6 7 8 9
  * lat      (lat) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Data variables:
    sfc_p    (lon, lat) int64 1031 980 1047 1028 988 ... 990 978 1026 1055 1003
    adv      (lon, lat) float64 -172.2 61.96 -109.3 22.05 ... -23.87 182.0 235.1