ncl计算整层水汽通量

zhangxiao

ncl绘制整层水汽通量的例子较少，特把写的代码贴出来，请大家帮助完善或者是指出问题，共同进步



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/wrf/WRFUserARW.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"

  begin
  a = addfile("/public3/home/interior/careeri/huser037/WRF/WRFV3S/run/wrfout_d02_2010-06-29_00:00:00","r")
   wks = gsn_open_wks("pdf","test1")

   au = wrf_user_getvar(a,"ua",(/47,263,1/))   ; U component of the wind at 10m
  av = wrf_user_getvar(a,"va",(/47,263,1/))   ; V component of the wind at 10m
  at = wrf_user_getvar(a,"tk",(/47,263,1/))   ; tk component of the wind at 10m
  arh = wrf_user_getvar(a,"rh",(/47,263,1/))   ; rh component of the wind at 10m
  ap  = wrf_user_getvar(a,"pressure",(/47,263,1/))
  p  = wrf_user_getvar(a,"pressure",48)
  sap  = wrf_user_getvar(a,"slp",(/48/))
  linlog=1
  pbot=1001

  printVarSummary(ap)
  printVarSummary(sap)
  sq=mixhum_ptrh(ap,at,arh,-2)
  mflux=sqrt(au^2+av^2)*sq/9.8
  printVarSummary(mflux)
  avg_flux=dim_avg_n_Wrap(mflux,0)
  avg_flux!0="lev"
  avg_flux!1="lat"
  avg_flux!2="lon"
  p!0="lev"
  p!1="lat"
  p!2="lon"
  printVarSummary(avg_flux)
  printVarSummary(p)
  ptop=min(p)
  psfc=sap
  ;csap=conform(p,sap,(/1,2/))
  ;printVarSummary(csap)
  vint = vibeta(p(lat|:,lon|:,lev|:),avg_flux(lat|:,lon|:,lev|:),linlog,psfc,pbot,ptop)
  printVarSummary(vint)
      
      res=True                  
      opts = res                          
      opts@cnFillOn = True  
      opts@pmLabelBarOrthogonalPosF = -0.1
      ;opts@ContourParameters = (/ 10., 90., 10./)
      opts@cnFillColors = (/"White","White","White", \
                              "Chartreuse","Chartreuse2","Green",\
                              "Green3","Green4", \
                              "ForestGreen","PaleGreen4"/)
      contour_psl = wrf_contour(a,wks,vint,opts)
      delete(opts)
  
     pltres = True                              ; Set plot options
     mpres = True                               ; Set map options
     mpres@mpDataBaseVersion="MediumRes"
     mpres@mpDataSetName="Earth..4"
     mpres@mpOutlineSpecifiers=(/"China:states"/)
   
   
    ; MAKE PLOTS                                       
      plot = wrf_map_overlays(a,wks,(/contour_psl/),pltres,mpres)
  


  end    希望大家指出问题，帮忙改正，此程序是绘制整层水汽通量

yangjunhua

谢谢楼主分享

dcx

请问楼主，这个数组是什么意思(/47,263,1/))  。

zhangxiao

随便选择的几个时间

denglin19901030

楼主为啥要做通量的平均值呢？

timevigor

感谢分享学习

timevigor

为什么用的是10m的UV分量

tianlanjiazi

感谢分享，学习学习

暮夜冬雨

{:5_235:}{:5_235:}{:5_235:}{:5_235:}

暮夜冬雨

楼主请问rh是从wrfout读出的吗，不记得wrfout可以输出这个量啊