关于python计算假相当位温

xpz0212

参考各个大佬们对于假相当位温计算（通过温度场，气压值，相对湿度场），（这是链接）
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

xpz0212

tulalang 发表于 2022-1-22 14:53
np数组不是可以直接运算么?为何上循环?是有哪方面的考虑呢

是的，可以直接计算，并且计算速度会快很多（是之前直接计算出过错，也不知道问题出在哪）

mezhe

学习了。。。。

puck66

收藏收藏

凉白开

感谢分享！收藏了，谢谢

tulalang

np数组不是可以直接运算么?为何上循环?是有哪方面的考虑呢

xpz0212

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

凉白开

学习了，感谢分享！

凨ぐ自ピ甴

大佬，我按照你的公式计算的ERA5的资料，算出来的结果与micaps的实况资料相差10-20左右，不知道什么原因，求指导啊

凨ぐ自ピ甴

发错了。。。