利用python计算假相当位温

yangzi88552

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

1. #新程序

复制代码

1. 之前论坛中提供的python程序：

复制代码

yangzi88552

代码没有显示出来：这是代码部分
新代码:
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

yangzi88552

计算中如果存在问题，请大家多多指出，谢谢~

一大碗年糕

感谢楼主，请问thse=thera*np.exp(((3376./LCL)-2.54)*q*(1.0+0.81*q))这里，Bolton的公式我看是除的抬升凝结温度，metpy返回两个变量是否应该是除lcl_temperature？

口蘑战士

python里面好像就直接有计算假相当位温的公式

wenbinz

口蘑战士 发表于 2023-6-25 17:02
python里面好像就直接有计算假相当位温的公式

那是计算相当位温的公式

wyy713

我用你的方法绘制了假相当位温的剖面图，然后跟metpy计算的相当位温比较，两者的数值差不多，图片也差不多

南南偶卖糕Oo

楼主您好，我的程序run不出来，一直显示：
DatetimeIndex(['2023-06-17 20:00:00'], dtype='datetime64[ns]', freq=None)

请问是什么情况呢

edwardli

wenbinz 发表于 2024-4-25 09:13
那是计算相当位温的公式

metpy的calc，计算的equivalent_potential_temperature()就是中国业务中常用的“假相当位温”