python 基于metpy计算位涡平流项（水平）

使用python 计算位涡水平平流项

使用数据：potential vorticity、水平风场分量：U、V
数据来源：

• ERA5
• hourly ：一天四次
• 空间分辨率：0.5x0.5
• 时间范围：2004-06
  使用的函数：
• metpy.calc.advection(scalar, u=None, v=None, w=None, *, dx=None, dy=None, dz=None, x_dim=- 1, y_dim=- 2, vertical_dim=- 3)
• metpy.calc.advection官方介绍
  与我之前的一个计算水汽散度的帖子的方法差不多：使用metpy计算水汽散度

编程需要注意的点：

• metpy.calc.advection()中传入的变量应该是二维的（lonxlat），而我们下载的位涡数据是四维的（time x level x lat x lon），所以需要做循环计算
• 使用metpy时需要给予各项附上响应的单位，防止计算的量纲出错
• metpy.calc.advection()中的x_dim=- 1, y_dim=- 2，分别对应了你的经度和纬度的维度，默认是经度为最后一维，纬度是倒数第二维
• 虽然metpy.calc.advection()提供了三个方向的计算，但是第三维的计算不太完善，因此本文以下只计算水平项
• 需要计算出dx、dy，也就是网格点的实际距离，这里metpy也给出了对应的函数：mpcalc.lat_lon_grid_deltas

代码如下：

import xarray as xr
import numpy as np
from metpy.units import units
import metpy.constants as constants
import metpy.calc as mpcalc
#########################文件路径##################################################
path=r'D:\PV\download_pv_hourly.nc'

path1=r'D:\PV\uv_0.5x0.5_2004-06.nc'
#########################读取pv############################################################
f_pv = xr.open_dataset(path)
pv   = f_pv.pv.sel(time=slice('2004-06-23','2004-06-30')).sortby('latitude')
pv_units=pv.units

lev  = pv.level

pressure = np.array(lev).reshape((1,37,1,1))*units.hPa
dpdl= np.gradient( pressure,axis=1)
##########################读取水平风场#########################################################
f = xr.open_dataset(path1).sortby('latitude').sel(time=slice('2004-06-23','2004-06-30'))

f = f.transpose('time','level','latitude','longitude')
f_u = f.u
f_v = f.v


lat =  f_u.latitude
lon =  f_u.longitude
time = f_u.time

dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)

pvadv=np.zeros((32,37,361,720))

for i in range(32):
    
    for j in range(37):
        
        print(i,j)  
        pvadv[i,j,:,:] = mpcalc.advection(pv[i,j,:,:],u=f_u[i,j,:,:],v=f_v[i,j,:,:],dx=dx,dy=dy,x_dim=-1,y_dim=-2)


pvadv=pvadv*units['K m**2 kg**-1 s**-1']/units.s
##########################计算平流垂向#####################################
# vertical
p_v=pv*units['K m**2 kg**-1 s**-1']
dpvdl=np.gradient( p_v,axis=1)

pvvdv = f_w*(units.Pa/units.s)*(dpvdl/dpdl)
#########################保存数据#############################################################
pvadv_nc=xr.Dataset({'pvadv':(('time','level','lat','lon'),pvadv)},
                    
                    coords = {
                        'time':time,
                        'level':lev.data,
                        'lat':lat.data,
                        'lon':lon.data,
                        
                        }
                    )
pvadv_nc.attrs['long_name']='potential vorticity advection'
pvadv_nc.to_netcdf('pvadv.nc')
###########################保存数据##########################################################
pvvdv_nc=xr.Dataset({'pvadv':(('time','level','lat','lon'),pvvdv.data)},
                    
                    coords = {
                        'time':time,
                        'level':lev.data,
                        'lat':lat.data,
                        'lon':lon.data,
                        
                        }
                    )
pvvdv_nc.attrs['long_name']='potential vorticity vertical advection'
pvvdv_nc.to_netcdf(r'pvvdv.nc')

Tips：

• 一般的计算结果显示的量级为10的-5次方到10的-8次方

• 虽然ERA5提供了位涡的数据，但是还是自己根据位涡的定义重新计算以下对比验证，如研究中高纬区域可能存在一些问题

    		    				一个努力学习python的ocean er
    	                    	水平有限，欢迎指正！！！
                            	欢迎评论、收藏、点赞、转发、关注。
                            	关注我不后悔，记录学习进步的过程~~