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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3711|回复: 14

[求助] 请问有人用python算过程整层水汽通量散度吗?算了一下感觉不太对

[复制链接]

新浪微博达人勋

发表于 2021-10-18 19:59:16 | 显示全部楼层 |阅读模式

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

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

x
u=xr.open_dataset('E:/all_data_1/data/uwnd/uwnd.2018.nc')
v=xr.open_dataset('E:/all_data_1/data/vwnd/vwnd.2018.nc')
sh=xr.open_dataset('E:/all_data_1/data/shum/shum.2018.nc')
lat=sh.lat
lon=sh.lon
lev=sh.level[0:8]
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
div_qv_gc=np.zeros((3,8,73,144))
qv_u_gc=np.zeros((3,8,73,144))
qv_v_gc=np.zeros((3,8,73,144))


qv_u_gc1=np.zeros((8,73,144))
qv_v_gc1=np.zeros((8,73,144))
# for j in pd.date_range(data1.loc[i]['start_day'],data1.loc[i]['end_day'],freq='D'):  
fu=u['uwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:]
fv=v['vwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:]

fu1=u['uwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:].mean(axis=0)
fv1=v['vwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:].mean(axis=0)


fq=sh['shum'].loc['2018-1-18':'2018-1-20',1000:300,:,:]
fq1=sh['shum'].loc['2018-1-18':'2018-1-20',1000:300,:,:].mean(axis=0)

qv_u_gc[:,:,:,:]=fu*fq/constants.g*100
qv_v_gc[:,:,:,:]=fv*fq/constants.g*100

qv_u_gc1[:,:,:]=fu1*fq1/constants.g*100
qv_v_gc1[:,:,:]=fv1*fq1/constants.g*100
        # 计算每次过程的每层水汽通量散度
for j in range(3):
    for k in range(lev.shape[0]):
        div_qv_gc[j,k,:,:] = mpcalc.divergence(u = qv_u_gc[j,k,:,:],v =qv_v_gc[j,k,:,:],dx = dx ,dy = dy)*100000
#计算整层水汽通量散度
div_total_gc=np.zeros((3,73,144))
for l in range(3):
    div_total_gc[l,:,:]= np.trapz(div_qv_gc[l,::-1],lev[::-1],axis=0)
div_mean=div_total_gc.mean(axis=0)
# div_all[i,:,:]=div_mean
# div1=div_all.mean(axis=0)         
div_ano1=xr.DataArray(div_mean , coords=[lat,lon], dims=['lat','lon'],name='div')
div_ano1.to_netcdf('E:/all_data_1/out_data/nc/qiang_year/vaper_div_2018.nc')     
        #计算整层散度通量
qv_u_gc_total=np.zeros((3,73,144))
qv_v_gc_total=np.zeros((3,73,144))
qv_u_gc_total=np.trapz(qv_u_gc1[::-1],lev[::-1],axis=0)
qv_v_gc_total=np.trapz(qv_v_gc1[::-1],lev[::-1],axis=0)

qu_ano=xr.DataArray(qv_u_gc_total, coords=[lat,lon], dims=['lat','lon'],name='qu')
qv_ano=xr.DataArray(qv_v_gc_total, coords=[lat,lon], dims=['lat','lon'],name='qv')
qu_ano.to_netcdf('E:/all_data_1/out_data/nc/qiang_year/qu_2018.nc')  
qv_ano.to_netcdf('E:/all_data_1/out_data/nc/qiang_year/qv_2018.nc')

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

新浪微博达人勋

发表于 2021-10-18 20:17:47 | 显示全部楼层
你单层的水汽通量就计算错了,水汽通量计算的时候应该是直接用q乘以u、q乘以v,不需要g,论坛上关于水汽通量的帖子,大多数在这个地方都是错的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-10-19 09:31:46 | 显示全部楼层
付亚男 发表于 2021-10-18 20:17
你单层的水汽通量就计算错了,水汽通量计算的时候应该是直接用q乘以u、q乘以v,不需要g,论坛上关于水汽通 ...

特别感谢呀
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 22:57:46 | 显示全部楼层
楼主计算正确以后继续分享一下嘛,先谢过
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-10-20 10:28:18 | 显示全部楼层
u=xr.open_dataset('F:/Rpython/lp36/data/water/uwnd.2018.nc')
v=xr.open_dataset('F:/Rpython/lp36/data/water/vwnd.2018.nc')
sh=xr.open_dataset('F:/Rpython/lp36/data/water/shum.2018.nc')
lat=sh.lat
lon=sh.lon
lev=sh.level[0:8]
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
div_qv_gc=np.zeros((3,8,73,144))
qv_u_gc=np.zeros((3,8,73,144))
qv_v_gc=np.zeros((3,8,73,144))
qv_u_gc1=np.zeros((8,73,144))
qv_v_gc1=np.zeros((8,73,144))
# for j in pd.date_range(data1.loc[i]['start_day'],data1.loc[i]['end_day'],freq='D'):  
fu=u['uwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:]
fv=v['vwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:]
fu1=u['uwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:].mean(axis=0)
fv1=v['vwnd'].loc['2018-1-18':'2018-1-20',1000:300,:,:].mean(axis=0)
fq=sh['shum'].loc['2018-1-18':'2018-1-20',1000:300,:,:]
fq1=sh['shum'].loc['2018-1-18':'2018-1-20',1000:300,:,:].mean(axis=0)
qv_u_gc[:,:,:,:]=fu*fq
qv_v_gc[:,:,:,:]=fv*fq
qv_u_gc1[:,:,:]=fu1*fq1
qv_v_gc1[:,:,:]=fv1*fq1
# 计算每次过程的每层水汽通量散度
for j in range(3):
    for k in range(lev.shape[0]):
        div_qv_gc[j,k,:,:] = mpcalc.divergence(u = qv_u_gc[j,k,:,:],v=qv_v_gc[j,k,:,:],dx = dx ,dy = dy)*100000
#计算整层水汽通量散度
div_total_gc=np.zeros((3,73,144))
for l in range(3):
    div_total_gc[l,:,:]= np.trapz(div_qv_gc[l,::-1],lev[::-1],axis=0)
# div_mean=div_total_gc.mean(axis=0)/constants.g*100
div_mean=div_total_gc.mean(axis=0)/9.8*100     
div_ano1=xr.DataArray(div_mean , coords=[lat,lon], dims=['lat','lon'],name='div')
div_ano1.to_netcdf('F:/Rpython/lp36/data/water/water2/vaper_div_2018.nc')   
#计算整层散度通量
qv_u_gc_total=np.zeros((3,73,144))
qv_v_gc_total=np.zeros((3,73,144))
# qv_u_gc_total=np.trapz(qv_u_gc1[::-1],lev[::-1],axis=0)/constants.g*100
# qv_v_gc_total=np.trapz(qv_v_gc1[::-1],lev[::-1],axis=0)/constants.g*100
qv_u_gc_total=np.trapz(qv_u_gc1[::-1],lev[::-1],axis=0)/9.8*100
qv_v_gc_total=np.trapz(qv_v_gc1[::-1],lev[::-1],axis=0)/9.8*100
qu_ano=xr.DataArray(qv_u_gc_total, coords=[lat,lon], dims=['lat','lon'],name='qu')
qv_ano=xr.DataArray(qv_v_gc_total, coords=[lat,lon], dims=['lat','lon'],name='qv')
qu_ano.to_netcdf('F:/Rpython/lp36/data/water/water2/qu_2018.nc')  
qv_ano.to_netcdf('F:/Rpython/lp36/data/water/water2/qv_2018.nc')
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-11-30 21:54:29 | 显示全部楼层
兄台 我有个疑问 div_total_gc[l,:,:]= np.trapz(div_qv_gc[l,::-1],lev[::-1],axis=0)一行中div_qv_gc是一个四纬向量,那么div_qv_gc[l,::-1]是指什么的?为啥只有两维了?请指教
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-12-1 17:04:38 | 显示全部楼层
晚晴山 发表于 2021-11-30 21:54
兄台 我有个疑问 div_total_gc[l,:,:]= np.trapz(div_qv_gc[l,::-1],lev[::-1],axis=0)一行中div_qv_gc是一 ...

div_total_gcd是三维,一个是天分开算,一个是平均的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-12-8 10:28:14 | 显示全部楼层
请问一下你这个算出来的整层积分水汽通量散度单位是啥呀?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-2-7 11:12:39 | 显示全部楼层
楼主 能不能吧文件发一下呀
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-7-16 00:25:46 | 显示全部楼层
请问一下,整层积分qv_u_gc_total=np.trapz(qv_u_gc1[::-1],lev[::-1],axis=0)/9.8*100 这里,为嘛要乘100 ?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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