use surface pressure to mask 4D netcdf variable

I've merged a 3D surface pressure field (ERA5, converted from Pa to hPa, function of lat,lon and time) with a 4D variable which is also a function of pressure levels (lat,lon,time,level).

So, my netcdf file has two fields, Temperature which is 4D:

float t(time, level, latitude, longitude)

surface pressure, which is 3d:

float sp(time, latitude, longitude)

The pressure dimension "level" is of course a vector:

int level(level)

What I want to do is make a mask for temperature for all locations where the pressure exceeds the surface pressure.

I know how to use nco to make a mask using a simple threshold:

 ncap2 -s 'mask=(level>800)' t_ps.nc mask.nc

But of course when I try to use the surface pressure

ncap2 -s 'mask=(level>sp)' t_ps.nc mask.nc

I get the error

ncap2: ERROR level and template sp share no dimensions

I think what I need to do is make a new variable like "level3d" which duplicates the pressure "level" to be a function of lat and lon, which I can then use to efficiently make the mask, yes? But I'm not sure how to do this with a dimension (I thought about cdo enlarge but couldn't get it to work).

By the way, instead of posting the data, this is the python api script I used to retrieve it

import cdsapi
c = cdsapi.Client()
c.retrieve(
    'reanalysis-era5-single-levels-monthly-means',
    {
        'format': 'netcdf',
        'product_type': 'monthly_averaged_reanalysis',
        'variable': 'surface_pressure',
        'year': '2020',
        'month': '03',
        'time': '00:00',
    },
    'ps.nc')

c.retrieve(
    'reanalysis-era5-pressure-levels-monthly-means',
    {
        'format': 'netcdf',
        'product_type': 'monthly_averaged_reanalysis',
        'variable': 'temperature',
        'pressure_level': [
            '1', '2', '3',
            '5', '7', '10',
            '20', '30', '50',
            '70', '100', '125',
            '150', '175', '200',
            '225', '250', '300',
            '350', '400', '450',
            '500', '550', '600',
            '650', '700', '750',
            '775', '800', '825',
            '850', '875', '900',
            '925', '950', '975',
            '1000',
        ],
        'year': '2020',
        'month': '03',
        'time': '00:00',
    },
    't.nc')

Solution

Your diagnosis of the NCO behavior is essentially correct. The "broadcast"

ncap2 -s 'mask=(level>sp)' t_ps.nc mask.nc

fails because level and sp are arrays (not scalars) that share no dimensions. The fix would be to create and use a temporary 3D version of level with something like

ncap2 -s 'level_3D[level,latitude,longitude]=level;mask=(level_3D>sp)' t_ps.nc mask.nc