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')

If I understand the question properly, you should be able to do this using CDO, using either aexpr or expr along with clev (which should give you access to the pressure level). If one of the variables is a surface variable then CDO should just use it for calculations for all vertical levels when you mix and matched 3D/4D variables. NCO, I guess, is a lot stricter and requires a 100% identical data structure.

The following is likely to work:

cdo -L -aexpr,'mask=(clev(Temperature)>sp)' infile outfile

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

THANKS so much to Robert for pointing me to clev, with this I found the following solution using enlarge too. It is a bit clunky, maybe someone has a more streamlined solution? This works directly now on the two fields downloaded in my example api script.

# make 3d level field
cdo enlarge,t.nc -expr,"l=clev(t)" t.nc level.nc

# merge together
cdo -O  merge t.nc ps.nc level.nc t_ps_p.nc

# make mask, ###note here I have a 100 factor to change Pa->hPa
cdo setctomiss,0 -expr,'mask=(l*100)<sp' t_ps_p.nc mask.nc

# mask the data
cdo mul t.nc mask.nc masked_t.nc