请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 15995|回复: 14

[经验总结] 关于python计算假相当位温

[复制链接]

新浪微博达人勋

发表于 2022-1-19 21:25:56 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 xpz0212 于 2022-1-19 21:25 编辑

参考各个大佬们对于假相当位温计算(通过温度场,气压值,相对湿度场),(这是链接)
https://mp.weixin.qq.com/s/mLuWEfD4hIt-4CSmphE8FAhttp://bbs.06climate.com/forum.php?mod=viewthread&tid=99672&fromuid=123871


稍微修改了一下,改成计算假相当位温的函数
可以直接调用,放入相对湿度和温度的二维数组,气压值,经向和纬向格点数,返回值为假相当位温值(单位:摄氏度)

下面是代码部分:

def Theta_se(Rh,t,p,nx,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))
   
    for i in range(ny):
        for j in range(nx):      
            
            lew[i,j] = 6.112*(np.e**((17.67*t[i,j])/(t[i,j]+243.5)))
            ee[i,j] = ((lew[i,j]*Rh[i,j])/100)/6.112
            ttd[i,j] = (243.5*(np.log(ee[i,j])))/(17.67-np.log(ee[i,j]))           
            Ltd[i,j] = L0-cl*(ttd[i,j]-t0)                  
            cttd[i,j] = e0*((t0/ttd[i,j])**(cl/Rw))*(np.e**((L0+cl*t0)*(ttd[i,j]-t0)/(Rw*ttd[i,j]*t0)))
            Ftd[i,j] = (0.622*Ltd[i,j]/(cpd*ttd[i,j]))-1            
            tc[i,j] = ttd[i,j]*Ftd[i,j]/(Ftd[i,j]+np.log(ttd[i,j]/t[i,j]))            
            Lc[i,j] = L0-cl*(tc[i,j]-t0)            
            thetad[i,j] = t[i,j]*((1000.0/(p-cttd[i,j]))**(Rd/cpd))
            wc[i,j] = 0.622*cttd[i,j]/(p-cttd[i,j])            
            thetase[i,j] = thetad[i,j]*(np.e**(wc[i,j]*Lc[i,j]/(cpd*tc[i,j])))  
    thetase=thetase-273.15
   
    return thetase



示例图.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2022-1-22 17:38:16 | 显示全部楼层
tulalang 发表于 2022-1-22 14:53
np数组不是可以直接运算么?为何上循环?是有哪方面的考虑呢

是的,可以直接计算,并且计算速度会快很多(是之前直接计算出过错,也不知道问题出在哪)
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2022-1-21 08:42:54 | 显示全部楼层
学习了。。。。
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2022-1-21 17:19:36 | 显示全部楼层
收藏收藏
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2022-1-22 11:39:14 | 显示全部楼层
感谢分享!收藏了,谢谢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-22 14:53:30 | 显示全部楼层
np数组不是可以直接运算么?为何上循环?是有哪方面的考虑呢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-1-22 17:40:45 | 显示全部楼层
xpz0212 发表于 2022-1-22 17:38
是的,可以直接计算,并且计算速度会快很多(是之前直接计算出过错,也不知道问题出在哪)

感谢楼上的建议,使用np数组直接计算更为方便,也更快

def Theta_se(Rh,t,p,nx,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-273.15
   
    return thetase
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-24 09:23:37 | 显示全部楼层
学习了,感谢分享!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-3-14 23:51:20 | 显示全部楼层
本帖最后由 凨ぐ自ピ甴 于 2022-3-15 20:44 编辑

大佬,我按照你的公式计算的ERA5的资料,算出来的结果与micaps的实况资料相差10-20左右,不知道什么原因,求指导啊
QQ截图20220315204227.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-3-14 23:56:58 | 显示全部楼层

发错了

本帖最后由 凨ぐ自ピ甴 于 2022-3-15 00:10 编辑

发错了。。。
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

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

本版积分规则

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

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

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