大气视热源Q1与视水汽汇Q2的计算 利用python以ERA5再分析资料为例

python代码： 

import os
import netCDF4 as nc
import numpy as np
import xarray as xr
import pandas as pd
import metpy.constants as constants
from metpy.units import units
from metpy.calc import mixing_ratio_from_specific_humidity, first_derivative, advection


# 加载NetCDF文件
ncfile_path_obs = r"F:\era5\7_pressurelevel.nc"  # 替换为您的文件路径
ncdata_obs = xr.open_dataset(ncfile_path_obs)

#读取数据
T_obs = ncdata_obs['t']
U_obs = ncdata_obs['u']
V_obs = ncdata_obs['v']
W_obs = ncdata_obs['w']
q_obs = ncdata_obs['q']

#计算dT/dt和dq/dt
dTdt = first_derivative(f=T_obs, axis='time')
dqdt = first_derivative(f=q_obs, axis='time')

#计算温度平流和水汽平流，计算结果带符号，为  —V·▽( )
temadvect = advection(scalar=T_obs, u=U_obs, v=V_obs)
qadvect = advection(scalar=q_obs , u=U_obs, v=V_obs)

#将气压reshape，便于之后数组的broadcast计算
pressure = ncdata_obs.level.values.reshape(1,23,1,1)*units.hPa

### 两种计算dT/dp dq/dp的方法，选一种即可
dTdp1 = first_derivative(f=T_obs, axis='level')
dTdp2 = np.gradient(f=T_obs,axis=1)*units.K/np.gradient(pressure,axis=1)
dqdp1 = first_derivative(f=q_obs, axis='level')
dqdp2 = np.gradient(f=q_obs,axis=1)*(units.kg/units.kg)/np.gradient(pressure,axis=1)


#计算平流advection() 偏导first_derivative()函数等 可以直接用T_obs这种Dataarray；
#做其他加减乘除运算时或者比湿dewpoint()函数等，dataarray需要转变成pint.Quantity
#要看具体的函数说明，一般普通气象要素的计算需要pint.Quantity；与x，y，z，t，plevels等维度有关的可以用Dataarray
σ= ((constants.dry_air_gas_constant* T_obs.values*units.K)/constants.dry_air_spec_heat_press/pressure) -  dTdp1
ver_t = W_obs.values/100*(units.hPa/units.s)*σ
ver_q = W_obs.values/100*(units.hPa/units.s)* dqdp1

#依照公式计算Q1和Q2 单位为w/kg
Q1 = constants.dry_air_gas_constant*(dTdt-ver_t-temadvect)
Q2 = constants.water_heat_vaporization*(-dqdt+qadvect-ver_q)

### Q1为 （96*23*281*601）time: 96，level: 23，latitude: 281，longitude: 601 
### ERA数据垂直层次本来就是从高到低 从100hPa到1000hPa  垂直积分不用逆序！！！
### 然后在垂直层次axis=1上积分
total_Q1 = np.trapz(Q1,ncdata_obs.level,axis=1)
total_Q2 = np.trapz(Q2,ncdata_obs.level,axis=1)

#使用numpy垂直积分没有单位，添加单位为w/m**2
total_Q1 = total_Q1*(units.W/units.m/units.m)
total_Q2 = total_Q2*(units.W/units.m/units.m)

大气视热源Q1与视水汽汇Q2的计算公式如下：

 

理论上对任何pressure_level的ERA5再分析资料都可以计算，数据维度要求为 时间*垂直层次*纬度*经度，也就是直接从ERA5网站上下载的nc数据类型。