爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3032|回复: 6

[源代码] 利用python计算假相当位温

[复制链接]

新浪微博达人勋

发表于 2023-6-18 20:16:46 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 yangzi88552 于 2023-6-18 21:03 编辑

在参考气象家园相关计算程序时,发现计算结果与实际结果差异过大,因此重新编写了假相当位温的计算程序。
由于假相当位温计算中存在经验公式,尤其是在计算抬升凝结高度时,由于有着不同的方法,所以可能导致结果存在一定差异。
  1. #新程序
复制代码
  1. 之前论坛中提供的python程序:
复制代码


新程序计算结果

新程序计算结果

原始数据结果

原始数据结果

论坛上之前提供的python程序计算结果

论坛上之前提供的python程序计算结果

tmp23061720.024.nc

1.02 MB, 下载次数: 10, 下载积分: 金钱 -5

850hPa温度

thetase23061720.024.nc

1.27 MB, 下载次数: 10, 下载积分: 金钱 -5

850hPa假相当位温

theta_se_forest_n.py

5 KB, 下载次数: 60, 下载积分: 金钱 -5

计算程序

rh23061720.024.nc

1.54 MB, 下载次数: 6, 下载积分: 金钱 -5

850hPa相对湿度

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2023-6-18 20:23:19 | 显示全部楼层
代码没有显示出来:这是代码部分
新代码:
def Theta_se(T0, RH, P0):
    Td = mpcalc.dewpoint_from_relative_humidity(T0, RH)
    thera = mpcalc.equivalent_potential_temperature(P0*units('hPa'),T0,Td)
    # q = mpcalc.mixing_ratio_from_relative_humidity(P0*units('hPa'), T0,RH)
    LCL, lcl_temperature = mpcalc.lcl(P0*units('hPa'),T0, Td)

    P0 = P0*100                        # units: Pa
    T0 = np.array(T0)
    Td = np.array(Td)
    LCL = np.array(LCL)
    tc=T0-273.15
    es = np.where(T0>=273.16,6.1078*np.exp(17.2693882*tc/(T0-35.86)),6.1078*np.exp(21.8745584*tc/(T0-7.66)))
    q = RH*(0.622*es/(P0-es))/100.0
    thse=thera*np.exp(((3376./LCL)-2.54)*q*(1.0+0.81*q))
    return thse

之前论坛中提供的参考python程序:
def Theta_se(Rh,t,p,nx,ny):
    nx,ny = len(nx),len(ny)
    #计算假相当位温
    #Rh:相对湿度(%)|t:温度(K)|p:气压值(hPa)|nx:经向格点数|ny:纬向格点数
    t0 = 273.16
    e0 = 6.1078
    L0 = 2500.79
    cl = 2.3697
    cpd = 1.0048
    Rd = 0.28704
    Rw = 0.4615
   
    #11个变量统一初始化成11个(ny,nx)二维数组
    cttd,Ftd,tc,Lc,thetad,wc,thetase,Ltd,ttd,lew,ee = np.ones((11,ny,nx))
   

    lew = 6.112*(np.e**((17.67*t)/(t+243.5)))
    ee = ((lew*Rh)/100)/6.112
    ttd = (243.5*(np.log(ee)))/(17.67-np.log(ee))           
    Ltd = L0-cl*(ttd-t0)                  
    cttd = e0*((t0/ttd)**(cl/Rw))*(np.e**((L0+cl*t0)*(ttd-t0)/(Rw*ttd*t0)))
    Ftd = (0.622*Ltd/(cpd*ttd))-1            
    tc = ttd*Ftd/(Ftd+np.log(ttd/t))            
    Lc = L0-cl*(tc-t0)            
    thetad = t*((1000.0/(p-cttd))**(Rd/cpd))
    wc = 0.622*cttd/(p-cttd)            
    thetase = thetad*(np.e**(wc*Lc/(cpd*tc)))  
    thetase=thetase
   
    return thetase
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-6-18 20:59:23 | 显示全部楼层
本帖最后由 yangzi88552 于 2023-6-18 21:01 编辑

计算中如果存在问题,请大家多多指出,谢谢~
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-6-22 23:47:33 | 显示全部楼层
感谢楼主,请问thse=thera*np.exp(((3376./LCL)-2.54)*q*(1.0+0.81*q))这里,Bolton的公式我看是除的抬升凝结温度,metpy返回两个变量是否应该是除lcl_temperature?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-6-25 17:02:37 | 显示全部楼层
python里面好像就直接有计算假相当位温的公式
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2024-4-25 09:13:00 | 显示全部楼层
口蘑战士 发表于 2023-6-25 17:02
python里面好像就直接有计算假相当位温的公式

那是计算相当位温的公式
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 3 天前 | 显示全部楼层
我用你的方法绘制了假相当位温的剖面图,然后跟metpy计算的相当位温比较,两者的数值差不多,图片也差不多
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表