假相当位温
镜像：气象分析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')).magnitudetc = 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</code></pre></div></div><div class="rno-markdown-code"><div class="rno-markdown-code-toolbar"><div class="rno-markdown-code-toolbar-info"><div class="rno-markdown-code-toolbar-item is-type"><span class="is-m-hidden">代码语言：</span>javascript</div></div><div class="rno-markdown-code-toolbar-opt"><div class="rno-markdown-code-toolbar-copy"><i class="icon-copy"></i><span class="is-m-hidden">复制</span></div></div></div><div class="developer-code-block"><pre class="prism-token token line-numbers language-javascript"><code class="language-javascript" style="margin-left:0"></code></pre></div></div><h4 id="8nqah" name="%E5%AF%BC%E5%85%A5%E5%BA%93%E5%B9%B6%E8%AF%BB%E6%95%B0%E6%8D%AE">导入库并读数据</h4><p>In [5]:</p><div class="rno-markdown-code"><div class="rno-markdown-code-toolbar"><div class="rno-markdown-code-toolbar-info"><div class="rno-markdown-code-toolbar-item is-type"><span class="is-m-hidden">代码语言：</span>javascript</div></div><div class="rno-markdown-code-toolbar-opt"><div class="rno-markdown-code-toolbar-copy"><i class="icon-copy"></i><span class="is-m-hidden">复制</span></div></div></div><div class="developer-code-block"><pre class="prism-token token line-numbers language-javascript"><code class="language-javascript" style="margin-left:0"></code></pre></div></div><div class="rno-markdown-code"><div class="rno-markdown-code-toolbar"><div class="rno-markdown-code-toolbar-info"><div class="rno-markdown-code-toolbar-item is-type"><span class="is-m-hidden">代码语言：</span>javascript</div></div><div class="rno-markdown-code-toolbar-opt"><div class="rno-markdown-code-toolbar-copy"><i class="icon-copy"></i><span class="is-m-hidden">复制</span></div></div></div><div class="developer-code-block"><pre class="prism-token token line-numbers language-javascript"><code class="language-javascript" style="margin-left:0">from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords,ALL_TIMESimport 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. 小结
谢谢观看，还请多多指正
计算不易，还请多多点赞收藏评论
