The dataset has surface pressure and temperature but I'm not sure how to calculate virtual temperature at the surface because the dataset is missing specific humidity at the surface.

It shouldn't be a problem to just use temperature instead. This would introduce a minor error in the surface pressure calculation, but this should be insignificant if pressure is just being used to determine whether reported data are actually under ground.

OK, I'll do that. Thanks

mask_below_surface (3).pdf
These are my notes on the approach I plan to take

from Travis

What you describe in Section 2 is one correct way of accomplishing this, and it is consistent with the way that heights are calculated in models. Section 3 probably isn't a good approach though, since models use an approximation that gravity is constant with height; also, the 'geopotential height' calculated in Section 2 is the actual, physical altitude of a given layer, so there's no need to go further than that. (Some sources refer to geopotential height as g*Z, where Z is the elevation, whereas others refer to it as Z.)

I would add that there are two other ways to accomplish the mask-below-surface operation, depending on what variables are available. In the following list, p_i corresponds to the pressure of a given mesh point (not including the surface pressure in this case), p_s corresponds to the surface pressure (psl in CMIP6/HighResMIP), p_msl corresponds to the mean sea level pressure field (psl in CMIP6/HighResMIP), t_s corresponds to the surface temperature (tas in CMIP6/HighResMIP), and z_s corresponds to the surface elevation

1. If the surface pressure field is available:
   if p_i > p_s: mask the level
2. If the mean sea level pressure field, surface temperature, and surface elevation are available:
   calculate p_s from p_msl, t_s, and z_s; if p_i > p_s: mask the level
3. If neither surface pressure is available, but temperature is:
   calculate geopotential height following your formulation in Section 2; if z_i < z_s: mask the level

Option 1 is preferable since it involves no approximations and is computationally the simplest. Option 2 would be the next most preferable since it involves only a minor approximation to infer surface pressure from the mean sea level pressure. I would consider Option 3 the last resort since it involves a numerical integration and requires more computation than either Option 1 or 2. Also, since the HighResMIP models only output a select set of levels, the numerical integration is not as accurate as the numerical integration that is performed internally in the model to determine mesh heights, so it would be more prone to error.

The attached document from Bill describes how p_s is calculated from p_msl (the end formula is the important one).
Adjusting_psl_to_ps-1.pdf

/global/project/projectdirs/m1949/observationDatasets/SURFACE/USGS_gtopo30_0.23x0.31_remap_c061107.nc

The following details which HighResMIP datasets provide simulation specific orography.

Th models I can find "orog" for are:

CMCC-CM2-HR4

CNRM-CM6-1

CNRM-CM6-1-HR

EC-Earth3

HadGEM3-GC31-HM

HadGEM3-GC31-LL

HadGEM3-GC31-MM

INM-CM5-H

IPSL-CM6A-LR

MPI-ESM1-2-HR

MRI-AGCM3-2-H

MRI-AGCM3-2-S

This is a pretty incomplete list (12 models out of 28 models I downloaded for HighResMIP).