计算整层水汽通量散度的疑问

haloyang

写了3个程序，①参考NCL官网http://www.ncl.ucar.edu/Applications/Scripts/mfc_div_2.ncl，利用函数vibeta垂直积分mfc_adv，mfc_div，然后相加得到整层水汽通量散度；②参考NCL官网http://www.ncl.ucar.edu/Applications/Scripts/mfc_div_1.ncl，计算uq，vq，再利用dim_sum_n(duvq_dpg, 1)垂直积分；③先计算uq，vq，再利用vibeta垂直积分；
结果：①和②图是相近的，而与③的结果是有差异的，请问这是为什么呢？按理说，▽Vq = V▽q+q▽V,①，②，③应该结果相近才是吧？
①的程序：
    ptop = 100              ; 'shum' upper level 300hPa                 
    ptop@units = "hPa"
    g    = 9.80665          ; m/s2

    yrStrt    = 197907
    yrLast    = 201407

    ;---ESRL: CDC data

    diri = "/mnt/f/circulation_data/"
    filq = "hus/JJA_MEAN_GLOBAL/OBS/hus_ERA-interim_JJA_mean_1979-2016.nc"
    filu = "ua/JJA_MEAN_GLOBAL/OBS/ua_ERA-interim_JJA_mean_1979-2016.nc"
    filv = "va/JJA_MEAN_GLOBAL/OBS/va_ERA-interim_JJA_mean_1979-2016.nc"
    filps= "ps/JJA_MEAN_GLOBAL/OBS/ps_ERA-interim_JJA_mean_1979-2016.nc" ;;;;地表气压

    pthu = diri+filu      
    pthv = diri+filv
    pthq = diri+filq     
    pthps= diri+filps   

    fu   = addfile(pthu ,"r")
    fv   = addfile(pthv ,"r")
    fq   = addfile(pthq ,"r")
    fps  = addfile(pthps,"r")

    ;---Time

    yyyymm  = cd_calendar(fu->time, -1)    ; ymd

: human readable
    iStrt   = ind(yyyymm.eq.yrStrt)             ; date for plotting and testing
    iLast   = ind(yyyymm.eq.yrLast)
               
    u    = fu->ua(iStrt:iLast,{1000:ptop},:,:)  ; m/s, (time,level,lat,lon)
    v    = fv->va(iStrt:iLast,{1000:ptop},:,:)
    q    = fq->hus(iStrt:iLast,{1000:ptop},:,:) ; [kg/kg], 1000-300 levels only
    ps   = fps->ps(iStrt:iLast,:,:)                    ; Pa=>[kg/(m-s2)], (time,lat,lon)， surface pressure data

    ;---Vertical levels

    ptop = ptop*100  ;;;30000
    ptop@units = "Pa"

    plev = doubletofloat(q&plev)                      ; hPa  [1000,925,850,700,600,500,400,300]
    plev = plev*100                     ; [100000,...,30000] Pa [kg/(m-s2)]
    plev@units = "Pa"

    ;---Change [kg/kg] to [g/kg]; not necessary: but common units for q

    q    = q*1000            
    q@units = "g/kg"

    ;************************************************
    ; Calculate the MFC_advection term
    ;   MFC_advect = (u*(dq/dx)+v*(dq/dy) )
    ; Internally, gradients are calculated via centered-finite_differences
    ;*************************************************  

    long_name = "MFC_advection"
    units     = "g/(kg-s)"       ; (m/s)*(g/kg)*(1/m) => (m/s)*(g/kg-m) => g/(kg-s)

    gridType  = 1   ; global fixed grid ordered S->N
    opt_adv   = 0   ; return only the advected variable; no gradients
    mfc_adv   = advect_variable_cfd(u, v, q, q&lat, q&lon, False, long_name, units, opt_adv)

    print("--------")

    ;************************************************
    ; Calculate the MFC_divergence term
    ;   MFC_div   = q*((du/dx)+(dv/dy) )      ; con(div)-vergence
    ;*************************************************  
    duv    = uv2dv_cfd(u, v, v&lat, v&lon,2)  

   
    mfc_div   = q*duv
    mfc_div@long_name = "mfc_divergence"
    mfc_div@units     = "g/(kg-s)"  
    copy_VarCoords(u,mfc_div)
    delete(duv)

    printVarSummary(mfc_div)  
    printMinMax(mfc_div, 0)  
    print("--------")

    ;---Integrated moisture flux components
    linlog        = 1
    pbot          = 1100*100                                             
    pbot@units    = "Pa"
    adv           = vibeta(plev,mfc_adv(time|:,lat|:,lon|:,plev|:),linlog,ps,pbot,ptop)/g


    adv@units     = "g/(m2-s)"

    adv          := tofloat(24*3600*adv/1000.) ;;;double
    adv@units     = "mm/day"
    adv@long_name = "Horizontal dry advecrtion"
    copy_VarCoords(v(:,0,:,:), adv)

    div           = vibeta(plev,mfc_div(time|:,lat|:,lon|:,plev|:),linlog,ps,pbot,ptop)/g
    div@units     = "g/(m2-s)"

    div           = 24*3600*div/1000
    div@units     = "mm/day"
    div@long_name = "moisture divergence associated with mass divergence"
    copy_VarCoords(v(:,0,:,:), div)

    printVarSummary(div)  
    printMinMax(div, 0)  
    print("--------")
    printVarSummary(adv)  
    printMinMax(adv, 0)  
    print("--------")

    ;;;;;计算MFC
    iduvq   = adv+div
    iduvq@units = "mm/day"
    iduvq@long_name = "Moisture flux divergence"
    copy_VarCoords(adv,iduvq)


②的程序：
  dp   = dpres_plevel_Wrap(plev, ps, ptop, 0) ; Pa;


  dpg  = dp/g    ;;;;dpg只是一个变量名，代表dp/g
  dpg@long_name = "Layer Mass Weighting"
  dpg@units     = "kg/m2"                     ; dp/g, Pa/(m s-2), reduce to kg m-2

;---Moisture flux components at each pressure level

  uq   = u*q                                  ; [time | 165] x [plev | 7] x [lat | 72] x [lon | 144]
  uq@long_name = "Zonal Moisture Flux [uq]"
  uq@units = "["+u@units+"]["+q@units+"]"     ; [m/s][g/kg]     
  copy_VarCoords(u,uq)                        ; (time,level,lat,lon)

  vq   = v*q                                  ; (:,:,:,:)
  vq@long_name = "Meridional Moisture Flux [vq]"
  vq@units = "["+v@units+"]["+q@units+"]"
  copy_VarCoords(v,vq)                        ; (time,level,lat,lon)

;---Integrated mass weighted moisture flux components
  uq_dpg = uq*dpg               

  iuq    = dim_sum_n(uq_dpg, 1)
  iuq@long_name = "Integrated Zonal UQ [uq*dpg]"
  iuq@LONG_NAME = "Sum: Mass Weighted Integrated Zonal Moisture Flux [uq*dpg]"
  iuq@units     = "[m/s][g/kg]"
  copy_VarCoords(u(:,0,:,:), iuq); (time,lat,lon)
  delete(uq_dpg)

  vq_dpg = vq*dpg              

  ivq    = dim_sum_n(vq_dpg, 1)  
  ivq@long_name = "Integrated Meridional VQ [vq*dpg]"
  ivq@LONG_NAME = "Sum: Mass Weighted Integrated Meridional Moisture Flux [vq*dpg]"
  ivq@units     = "[m/s][g/kg]"
  copy_VarCoords(v(:,0,:,:), ivq)
  delete(vq_dpg)

;;----------计算整层比湿-------------------------------------------------------
  qdp  = q*dp
  iq   = dim_sum_n(qdp, 1)  
  copy_VarCoords(q(:,0,:,:), iq);
   printVarSummary(iq)

;---Divergence of moisture flux: uv2dvF => global 'fixed' rectilinear grid
  duvq  = uv2dv_cfd(uq, vq, uq&lat, uq&lon,2)  
  duvq@long_name = "Divergence of Moisture Flux"
  duvq@units     = "g/(kg-s)"    ; (1/m)*[(m/s)(g/kg)] => [g/(kg-s)]

;---Mass weighted integration [sum] of the divergence of moisture flux

  duvq_dpg = duvq*dpg            ;  [g/(kg-s)][kg/m2] => [g/(m2-s)]

  iduvq    = dim_sum_n(duvq_dpg, 1)
  iduvq@long_name = "Integrated Mass Wgt MFC"
  iduvq@LONG_NAME = "Integrated Mass Weighted Moisture Flux Convergence"
  iduvq@units     = "g/(m2-s)"
  copy_VarCoords(u(:,0,:,:), iduvq)      ; (time,lat,lon)
  delete(duvq_dpg)

  iduvq  = 24*3600*iduvq/1000
  iduvq@units = "mm/day"


③的程序
  uq   = u*q                                 
  uq@long_name = "Zonal Moisture Flux [uq]"
  uq@units = "["+u@units+"]["+q@units+"]"     ; [m/s][g/kg]     
  copy_VarCoords(u,uq)                        ; (time,level,lat,lon)

  vq   = v*q                                  ; (:,:,:,:)
  vq@long_name = "Meridional Moisture Flux [vq]"
  vq@units = "["+v@units+"]["+q@units+"]"
  copy_VarCoords(v,vq)                        ; (time,level,lat,lon)

PRINT_RAW = True
if (PRINT_RAW) then
  printVarSummary(uq)                       
  printMinMax(uq,0)
  print("-----")
  printVarSummary(vq)
  printMinMax(vq,0)
  print("-----")
end if

;---Integrated moisture flux components
  linlog        = 1
  pbot          = 1100*100                                             
  pbot@units    = "Pa"
  iuq           = vibeta(plev,uq(time|:,lat|:,lon|:,plev|:),linlog,ps,pbot,ptop)/g
  iuq@long_name = "Integrated Zonal UQ [uq*dpg]"
  copy_VarCoords(u(:,0,:,:), iuq); (time,lat,lon)

  ivq           = vibeta(plev,vq(time|:,lat|:,lon|:,plev|:),linlog,ps,pbot,ptop)/g
  ivq@long_name = "Integrated Meridional VQ [vq*dpg]"
  ivq@LONG_NAME = "Sum: Mass Weighted Integrated Meridional Moisture Flux [vq*dpg]"
  ivq@units     = "[m/s][g/kg]"
  copy_VarCoords(v(:,0,:,:), ivq); (time,lat,lon)

;;----------计算整层比湿-------------------------------------------------------

;---Divergence of moisture flux: uv2dvF => global 'fixed' rectilinear grid
  iduvq  = uv2dv_cfd(iuq, ivq, ivq&lat, ivq&lon,2)  

  iduvq@long_name = "moisture flux Divergence"
  iduvq@LONG_NAME = "Moisture Flux Divergence"


  iduvq           = 24*3600*iduvq/1000.
  iduvq@units     = "mm/day"
  copy_VarCoords(u(:,0,:,:), iduvq)      ; (time,lat,lon)

WZZ809525579

楼主找到图不一样的原因了吗

凉快了吧

您好，最近也准备用era5的数据分析水汽通量，想请问hus_ERA-interim_JJA_mean_1979-2016.nc这个数据是啥数据啊，还有ps_ERA-interim_JJA_mean_1979-2016.nc这个地表气压是选择的是1000hpa的吗，期待您的回复

haloyang

凉快了吧 发表于 2021-3-18 10:35
您好，最近也准备用era5的数据分析水汽通量，想请问hus_ERA-interim_JJA_mean_1979-2016.nc这个数据是啥数 ...

hus是比湿，ps就是地表气压，不是1000hPa。青藏高原地表气压500多hPa

凉快了吧

haloyang 发表于 2021-3-19 12:32
hus是比湿，ps就是地表气压，不是1000hPa。青藏高原地表气压500多hPa

好的，好的！非常感谢您的回复

凉快了吧

haloyang 发表于 2021-3-19 12:32
hus是比湿，ps就是地表气压，不是1000hPa。青藏高原地表气压500多hPa

好的，好的！非常感谢您的回复

凉快了吧

您好，不好意思打扰您一下，就是在用①的程序运行时，ncl报了这样的错
(0)     advect_variable_cfd: grid is not in South-to-North order.
(0)     advect_variable_cfd fatal error(s) encountered: ier=10000
我使用的是era5的逐月再分析资料 ，其中u,v,q 都是
[time | 111] x [level | 16] x [latitude | 181] x [longitude | 360]
这样的。
想问问您出现过这样的错误吗，谢谢

呆妹小霸王

您好，我可以要一份您使用的数据吗？我最近也在计算水汽通量的散度，但我不知道自己的程序到底对不对(大概率是错的)，我想重复一下看能不能画出和您一样的图。
我不会用NCL，不知道这两种dim_sum_n(duvq_dpg, 1)；vibeta积分函数有什么不一样

凉快了吧

凉快了吧 发表于 2021-3-21 21:16
您好，不好意思打扰您一下，就是在用①的程序运行时，ncl报了这样的错
(0)     advect_variable_cfd: grid ...

u（：，：，：：-1，：）把u的纬度从North to South变成South to North

YQY优优

楼主找到图不一样的原因了吗

来自苹果APP客户端