画ERA5气压层剖面图（含地形）

气象上一般使用气压垂直坐标系，在不同的气压层绘制变量。ERA5再分析数据的最低气压层是1000 hPa。实际上，由于地形起伏，一些位置的地面气压低于1000 hPa，一些位置的地面气压高于1000 hPa。所以，1000 hPa等压面所对应的高度位置，有时候在地面之下，有时候在地面之上。对于那些位于地面之下的数据，没有实际意义的，需要进行剔除。

那么，如何判断等压面的数据是否在地面之上，还是在地面之下呢。

常规的思维，如前所述，就是比较等压面的气压P_level和该位置正下方地面的地面气压P_sp（Surface Pressure），如果等压面气压小于地面气压（即P_level<P_sfc），则该等压面位于地面之上，该气压层的气象要素（比如温度，风速）有实际意义。

但是，实际工作中，一般用的是海平面气压P_sfc，地面气压P_sp这个量一般用不到，尽管ERA5数据有这个量。所以利用气压进行地面高度判断不是很方便。

那么，有什么比较方便的办法呢？使用位势/位势高度进行比较。这是因为观测气压比测量高度更加方便，并且气压一般是随高度递减的，气压和高度一一对应。这也就说明了，气象上通常不画不同高度上的等压线，而画不同等压面的上的等位势/位势高度线。

如何比较呢，可以通过比较不同气压层上某点对应的位势Z和该位置正下方的地面的位势Z_sfc，可以判断该等压面是否在地面之上。若Z>Zsfc，则位于地面上，Z<Zsfc，则位于地面下。

这里以GrADS软件，ERA5再分析数据为例，讲如何绘制带地形覆盖的剖面图（垂直坐标为气压）。

需要用到变量有两个，地面位势，不同等压面的位势。

其中，ERA5的地面位势是"ERA5 hourly data on single levels from 1959 to present"的数据集中，名称为Geopotential的变量（在ncl中变量名为 Z_GDS0_SFC），单位是 m**2 s**-2，表示各格点的地面层所在的位势。这是一个二维变量。

而不同等压面的位势Z是在"ERA5 hourly data on Pressure levels from 1959 to present"数据集中，名称为Geopotential的变量（在ncl中变量名为 Z_GDS0_ISBL_**），单位是 m**2 s**-2，表示各格点的地面层所在的位势。这是一个三维变量。

在GrADS中，二维变量和三维变量不能直接比较（没有类似python的broadcast广播机制），所以需要另外定义一个三维变量，把二维的地形数据spread扩展到三维，再进行比较。

'set z 1 24'

'define z=zsfc'

'd maskout(zprs,zprs-z)'

这样就可以将地形以下的数据mask掉了