求助NCL计算NCEP逐日6小时水汽通量问题

sunyue19901203

小弟我在用NCEP逐日6小的资料计算水汽通量时遇到点问题，想请教下论坛里的各位大神，我用的是2011年u,v,q,和地面气压p的资料，选取了东亚地区夏季92天，每天第一个时刻，高度为1000,925,850,700,6005个高度的资料进行分析，其中，u和v的单位是m/s，q的单位是kg/kg，p的单位是Pa,下面是我的代码

1. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
2. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   
3. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
   
4. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
   
5. ;============================================================================
   
6. 
   
7. begin
   
8. 
   
9.   u_filepath = "/cma/g4/expert3/sunyue/diag/ncl/model_prep/NCEP/4daily/uwnd/uwnd.2011.nc"
   
10.   v_filepath = "/cma/g4/expert3/sunyue/diag/ncl/model_prep/NCEP/4daily/vwnd/vwnd.2011.nc"
    
11.   q_filepath = "/cma/g4/expert3/sunyue/diag/ncl/model_prep/NCEP/4daily/shum/shum.2011.nc"
    
12.   p_filepath = "/cma/g4/expert3/sunyue/diag/ncl/model_prep/NCEP/4daily/pressure/pres.sfc.2011.nc"
    
13. 
    
14.   u_file     = addfile(u_filepath,"r")
    
15.   v_file     = addfile(v_filepath,"r")
    
16.   q_file     = addfile(q_filepath,"r")
    
17.   p_file     = addfile(p_filepath,"r")
    
18. 
    
19.   u  =  u_file->uwnd(604:971,0:4,{15:50},{70:130})
    
20.   v  =  v_file->vwnd(604:971,0:4,{15:50},{70:130})
    
21.   q  =  q_file->shum(604:971,0:4,{15:50},{70:130})
    
22.   p  =  p_file->pres(604:971,{15:50},{70:130})
    
23. ;===================================================================================================
    
24. 
    
25.   dsize = dimsizes(u)
    
26.   mlat  = dsize(2)
    
27.   nlon  = dsize(3)
    
28.   u_day = new((/92,4,5,mlat,nlon/),"float")
    
29.   v_day = new((/92,4,5,mlat,nlon/),"float")
    
30.   q_day = new((/92,4,5,mlat,nlon/),"float")
    
31.   p_day = new((/92,4,mlat,nlon/),"float")
    
32.   
    
33.   u_day = onedtond(ndtooned(u),(/92,4,5,mlat,nlon/))
    
34.   v_day = onedtond(ndtooned(v),(/92,4,5,mlat,nlon/))
    
35.   q_day = onedtond(ndtooned(q),(/92,4,5,mlat,nlon/))
    
36.   p_day = onedtond(ndtooned(p),(/92,4,mlat,nlon/))
    
37. ;==== time = 00
    
38.   u_day_avg = dim_avg_n(u_day,0)
    
39.   v_day_avg = dim_avg_n(v_day,0)
    
40.   q_day_avg = dim_avg_n(q_day,0)
    
41.   p_day_avg = dim_avg_n(p_day,0)
    
42.   
    
43.   u_00  = u_day_avg(0,:,:,:)
    
44.   v_00  = v_day_avg(0,:,:,:)
    
45.   q_00  = q_day_avg(0,:,:,:)
    
46.   p_00  = p_day_avg(0,:,:)
    
47. 
    
48.   qu1   = q_00 * u_00
    
49.   qv1   = q_00 * v_00
    
50. 
    
51.   printVarSummary(qu1)
    
52.   
    
53.   qu1!0   = "lev"
    
54.   qu1&lev = u&level
    
55.   qu1&lev@units = "hPa"
    
56.   qu1!1   = "lat"
    
57.   qu1&lat = u&lat
    
58.   qu1&lat@units = "degrees_north"
    
59.   qu1!2   = "lon"
    
60.   qu1&lon = u&lon
    
61.   qu1&lon@units = "degrees_east"  
    
62.   
    
63.   qv1!0   = "lev"
    
64.   qv1&lev = v&level
    
65.   qv1&lev@units = "hPa"
    
66.   qv1!1   = "lat"
    
67.   qv1&lat = v&lat
    
68.   qv1&lat@units = "degrees_north"
    
69.   qv1!2   = "lon"
    
70.   qv1&lon = v&lon
    
71.   qv1&lon@units = "degrees_east"
    
72. 
    
73.   qu11 = qu1(lat|:,lon|:,lev|:)
    
74.   qv11 = qv1(lat|:,lon|:,lev|:)
    
75.   delete(qu1)
    
76.   delete(qv1)
    
77.   
    
78.   pres   = (/1000,925,850,700,600/)
    
79.   linlog = 1
    
80. 
    
81.   data11 = (vibeta(pres,qu11,linlog,p_00,1000,600))/9.8
    
82.   data12 = (vibeta(pres,qv11,linlog,p_00,1000,600))/9.8
    
83. 
    
84.   mm = sqrt(data11^2+data12^2)
    
85. 
    
86.   mm!0   = "lat"
    
87.   mm&lat = u&lat
    
88.   mm&lat@units = "degrees_north"
    
89.   mm!1   = "lon"
    
90.   mm&lon = u&lon
    
91.   mm&lon@units = "degrees_east"

复制代码

然后算出来水汽通量的量级最大只有4.4，感觉好小，请问我的程序有没有什么问题;
还有我看有的帖子在积分时将地面气压除了100，换算成hPa，请问需要么？

timevigor

请问UV都是4维的数据，怎么在48.49行处理成3维的？

sunyue19901203

timevigor 发表于 2015-10-7 19:35
请问UV都是4维的数据，怎么在48.49行处理成3维的？

在45,46两行我不是选了时间维是第一个时间么，这不就变成三维的了

gaofengjlsy

谢谢！！！

timevigor

还有再请问一下你最后算出的水汽通量的单位是什么

sunyue19901203

timevigor 发表于 2015-10-14 16:26
还有再请问一下你最后算出的水汽通量的单位是什么

这个我也没搞清楚。。。

luoxincheng

刚开始学NCL，先看看大神们编写的脚本吧