NCL 计算整层水汽通量

奚子惠

用NCL计算整层水汽通量报错：



部分脚本如下：
     nlat=73
     nlon=144
     nyear=35   
     nlev=8         
     p       = (/ 1000.,925.,850.,700.,600.,500., \
               400.,300./)
         linlog = 1
     pbot   = 1000.                                             
     ptop   = 300.
         u_file=addfile("uwnd_Dtrend.nc","r")
         v_file=addfile("vwnd_Dtrend.nc","r")
         p_file=addfile("pres_Dtrend.nc","r")
         shum_file=addfile("shum_Dtrend.nc","r")
         uu=u_file->uwnd_Dtrend
         vv=v_file->vwnd_Dtrend
         pres=p_file->pres_Dtrend
         shum=shum_file->shum_Dtrend
       
        ;;; 计算qu = q(比湿)乘以u(u风) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
       
         qu = new((/nyear,nlev,nlat,nlon/),float)
          qv = new((/nyear,nlev,nlat,nlon/),float)
         qu=shum*uu
         qv=shum*vv
         copy_VarCoords(uu,qu)
         copy_VarCoords(vv,qv)
         qu@_FillValue=-9.99E+33  
         qv@_FillValue=-9.99E+33  
        printVarSummary(qu)

         u = vibeta (p,qu(time|:,lat|:,lon|:,level|:),linlog,pres,pbot,ptop)/9.8  ; returns qu(time,lat,lon)
         v = vibeta (p,qv(time|:,lat|:,lon|:,level|:),linlog,pres,pbot,ptop)/9.8  ;;单位是g/s`hPa`cm
       
         copy_VarCoords(uu(:,0,:,:),u)
         copy_VarCoords(vv(:,0,:,:),v)
         u@_FillValue=-9.99E+33  
         v@_FillValue=-9.99E+33  
         
       
         printVarSummary(v)

请教各位大神们知道是什么原因嘛~多谢！

奚子惠

hxyj 发表于 2015-6-21 09:51
请问楼主是否已解决该问题了？能否将计算的完整脚本贴出共享？

是的已经解决了
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
;************************************************
; open file and read in variable
;***********************************************
;levels 1000 925 850 700 600 500 400 300  
; 层数    0   1   2   3   4   5   6   7      
;***********************************************
     nlat=73
     nlon=144
         nyear=35   
     nlev=8         
         
         p       = (/ 1000.,925.,850.,700.,600.,500., \
               400.,300.,250.,200.,150.,100., \
                70.,50.,30.,20.,10. /)
         linlog = 1
     pbot   = 1100.                                             
     ptop   = 300.
         u_file=addfile("uwnd_Dtrend.nc","r")
     v_file=addfile("vwnd_Dtrend.nc","r")
         p_file=addfile("pres_Dtrend.nc","r")
         shum_file=addfile("shum_Dtrend.nc","r")
         xx_file=addfile("uwnd.mon.mean.nc","r")
         
         uu=u_file->uwnd_Dtrend
         vv=v_file->vwnd_Dtrend
         pres=p_file->pres_Dtrend
         shum=shum_file->shum_Dtrend
       
         qu = new((/nyear,nlev,nlat,nlon/),float)
     qv = new((/nyear,nlev,nlat,nlon/),float)
       
         qu=shum*uu
         qv=shum*vv
         copy_VarCoords(uu,qu)
         copy_VarCoords(vv,qv)
         qu@_FillValue=-9.99E+33  
         qv@_FillValue=-9.99E+33  
       
         qu_new = new((/nyear,nlev+9,nlat,nlon/),float)
     qv_new = new((/nyear,nlev+9,nlat,nlon/),float)

         do i=0,7
            qu_new(:,i,:,:)= qu(:,i,:,:)
                qv_new(:,i,:,:)= qv(:,i,:,:)
         end do
         
         do j=8,nlev+9-1
           qu_new(:,j,:,:)= 0.0
           qv_new(:,j,:,:)= 0.0
         end do
         qu_new&level = xx_file->uwnd&level
         qv_new&level = xx_file->uwnd&level
        ; printVarSummary(qu_new)
         
         u = vibeta (p,qu_new(time|:,lat|:,lon|:,level|:),linlog,pres,pbot,ptop)/9.8  ; returns u(time,lat,lon)
         v = vibeta (p,qv_new(time|:,lat|:,lon|:,level|:),linlog,pres,pbot,ptop)/9.8  
       
         copy_VarCoords(uu(:,0,:,:),u)
         copy_VarCoords(vv(:,0,:,:),v)
         u@_FillValue=-9.99E+33  
         v@_FillValue=-9.99E+33  
       
        ; printVarSummary(v)
         asciiwrite ("allqu.txt" , u)
         delete(uu)
         delete(vv)
         delete(pres)
     delete(shum)

end

andrewsoong

没看出这部分代码什么问题~~~请斑竹来说

qingchun

是ntime=35？

风往北吹

看看pres=p_file->pres_Dtrend中是Pa还是hPa。

奚子惠

风往北吹 发表于 2015-5-19 14:51
看看pres=p_file->pres_Dtrend中是Pa还是hPa。

是millibars
           

奚子惠

qingchun 发表于 2015-5-19 13:50
是ntime=35？

这只是变量名而已 应该没关系吧

longlivehj

vibeta采用的是中点公式，需要3个点进行积分。出现这种错误的原因是，地形以上不足3层。可以人为补几层上去。

奚子惠

longlivehj 发表于 2015-5-19 15:50
vibeta采用的是中点公式，需要3个点进行积分。出现这种错误的原因是，地形以上不足3层。可以人为补几层上去 ...

你的意思是将qu，qv再人为的往上补几层吗？那补的那几层是设为缺省还是0比较好？因为我下面还遇到一个问题，就是用完vibeta后，数据出现了部分缺省值，然后我用uv2dv_cfd算完散度之后，想分解成辐散风和旋转风，在用到 dv2uvg(div,ud,vd)这句后 ，算出来的ud全部变成了缺省值，导致图画不出来

奚子惠

longlivehj 发表于 2015-5-19 15:50
vibeta采用的是中点公式，需要3个点进行积分。出现这种错误的原因是，地形以上不足3层。可以人为补几层上去 ...

我查了一下dv2uvg的说明，它的英文意思是“不处理缺省值，原本是缺省值的格点最后那个格点上输出的还是缺省值”，不知道我理解的对不对，但是也不应该我算出来的ud和vd所有格点全部是缺省值诶····

linan404x

新手上路，围观学习