如何计算WRF台风模拟的假相当位温
如何计算WRF台风模拟的假相当位温
假相当位温
镜像：气象分析3.9 数据：wrf台风模拟数据
🎓 作者：酷炫用户名
📢 版权声明：公益性质转载需联系作者本人获取授权。转载本文时，请务必文字注明“来自：和鲸社区：酷炫用户名”，并附带本项目超链接。
温馨提示
由于可视化代码过长隐藏，可点击如何计算WRF台风模拟的假相当位温运行Fork查看 🔜🔜若没有成功加载可视化图，点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可
前言
为什么写这个
前几日有读者来信想看看假相当位温的计算。我寻思之前写了位温、相当位温，那肯定不能少了假相当位温。凑个位温三兄弟。
故而去查了些资料，准备计算
经典文献中假相当位温定义如下
where T is the temperature (K), p is the pressure (Pa),p0 is the reference pressure (100000 Pa), Tc is the condensation temperature (obtainable from the dewpoint formula) and rv is the water vapor mixing ratio (kg kg-1). When rv = 0, θep = θ, the potential temperature.
Bolton, D. 1980. The computation of equivalent potential temperature. Mon. Wea. Rev.. 108. 1046–1053
当然你在中国气象局的网页看到更简洁的公式
表达式中r为混合比，可见是温度、气压、水汽含量的函数，表示温、压、湿综合的物理量，是预报业务中常用的重要物理量。在同一气压条件下，假相当位温越大空气越暖湿，反之，空气越干冷。850hPa的假相当位温的分布与大小是预报员常关注的重点。暴雨时850hPa的假相当位温值一般在330K以上。反映大气中层结潜在不稳定，暴雨落区一般在0～15K之间。
这次我们就用简短的公式计算试试
关于相当位温和假相当位温的差别
如果你读过Bolton的文献，第一句就说相当位温，又称假相当位温。
也可能看过部分文章认为两者是一个东西。
注意该文献已经是四十年前的，需要用发展的眼光看待定义。
两者的区分可以看这篇微信公众号文章再话“相当位温与假相当位温”
1. 定义函数
In [3]:
import numpy as np
from metpy.calc import vapor_pressure,lcl
from metpy.units import units
def calculate_theta_se(T, p, r,td):
"""
计算假相当位温
参数:
T: 温度（摄氏度）
p: 大气压力（hPa）
r: 混合比
tc : 绝热饱和温度 LCL
返回值:
假相当位温（J/kg）
"""
# 常量
Rd = 287 # J/kg/K
cp_d = 1004 # J/kg/K
L = 2.555*1e6 # J/kg
e = vapor_pressure(p * units.hPa, r * units('g/kg')).magnitude
print('水汽压',e)
tc = lcl(p* units.hPa, (T-273.15) * units.degC, td * units.degC)[1].magnitude
print('lcl',tc)
# 计算假相当位温
theta_se = T * (1000 / (p - e) ) **( Rd / cp_d) * np.exp(L *r / (cp_d * tc))
return theta_se
# 使用示例
T = 230 # 温度
p = 850 # 大气压力
r = 0.005 # 混合比
td = 15
theta_se = calculate_theta_se(T, p, r,td)
print("假相当位温:", theta_se, "J/kg")水汽压 0.0068332158790973254
lcl 31.172628864169724
假相当位温: 362.3898746288438 J/kg2. 实际应用 : WRF后处理提取相关变量计算假相当位温
还是从老伙计wrfout中提取需要的变量：温度 气压 混合比 等等
设置函数
In [4]:
def calculate_theta_se_wrf(T, p, r,td):
"""
计算假相当位温
参数:
T: 温度（摄氏度）
p: 大气压力（hPa）
r: 混合比
tc : 绝热饱和温度 LCL
返回值:
假相当位温（J/kg）
"""
# 常量
Rd = 287 # J/kg/K
cp_d = 1004 # J/kg/K
L = 2.555*1e6 # J/kg
e = vapor_pressure(p * units.hPa, r * units('g/kg')).magnitude
tc = lcl(p* units.hPa, (T-273.15) * units.degC, td * units.degC)[1].magnitude
# 计算假相当位温
theta_se = T * (1000 / (p - e) ) **( Rd / cp_d) * np.exp(L *r / (cp_d * tc))
return theta_se导入库并读数据
In [5]:
from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords,ALL_TIMES
import numpy as np
from netCDF4 import Dataset
import xarray as xr
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cmaps
from glob import glob
import metpy.calc as mpcalc
import metpy.constants as constants
import os
# 定义 WRF 文件夹路径和文件名前缀
wrfout_path = "/home/mw/input/typhoon9537"
filename_prefix = "wrfout_d01_"
wrf_files = sorted([os.path.join(wrfout_path, f) for f in os.listdir(wrfout_path) if f.startswith(filename_prefix)])
wrf_list = [Dataset(f) for f in wrf_files]In [6]:
p_wrf = getvar(wrf_list, 'pressure', timeidx=-1).data
t = getvar(wrf_list, 'tk', timeidx=-1).data
r = (getvar(wrf_list, 'QVAPOR', timeidx=-1)/1000).data
td = getvar(wrf_list, 'td', timeidx=-1).data
lons = getvar(wrf_list, 'lon', timeidx=0)
lats = getvar(wrf_list, 'lat', timeidx=0)In [7]:
result_wrf = calculate_theta_se_wrf(t,p_wrf,r,td)
result_wrf.shapeOut[7]:
(34, 437, 447)插值850hpa
In [8]:
se850 = interplevel(result_wrf, p_wrf, 850.0)绘图部分
3. metpy的相当位温
In [10]:
from metpy.calc import equivalent_potential_temperature
from metpy.units import units
tt = getvar(wrf_list, 'tc', timeidx=0)
result2 = equivalent_potential_temperature(p_wrf*units.hPa, tt *units.degC, td*units.degC)
se2 = interplevel(result2, p_wrf, 850.0)
se2
4. 小结
谢谢观看，还请多多指正 计算不易，还请多多点赞收藏评论