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

cl22

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')

付亚男

你单层的水汽通量就计算错了，水汽通量计算的时候应该是直接用q乘以u、q乘以v，不需要g，论坛上关于水汽通量的帖子，大多数在这个地方都是错的

cl22

付亚男 发表于 2021-10-18 20:17
你单层的水汽通量就计算错了，水汽通量计算的时候应该是直接用q乘以u、q乘以v，不需要g，论坛上关于水汽通 ...

特别感谢呀

puck66

楼主计算正确以后继续分享一下嘛，先谢过

cl22

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')

晚晴山

兄台 我有个疑问 div_total_gc[l,:,:]= np.trapz(div_qv_gc[l,::-1],lev[::-1],axis=0)一行中div_qv_gc是一个四纬向量，那么div_qv_gc[l,::-1]是指什么的？为啥只有两维了？请指教

cl22

晚晴山 发表于 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是三维，一个是天分开算，一个是平均的

nuistdgy

请问一下你这个算出来的整层积分水汽通量散度单位是啥呀？

haodong

楼主 能不能吧文件发一下呀

puck66

请问一下，整层积分qv_u_gc_total=np.trapz(qv_u_gc1[::-1],lev[::-1],axis=0)/9.8*100 这里，为嘛要乘100 ？