用grads和ncl画水汽通量图出现的差异

wow哈哈

相同的计算公式，相同的地图区域，为什么二者画出来的图差异很大，grads应该是对的，那么求教高手为什么ncl出问题了。（还有小问题，参考的别的高手的脚本，不明白#是什么意思，和*一样吗）
Grads脚本：
'reinit'
'open d:\study\data\201607\fnl_20160704_00_00.ctl'
'enable print d:\study\data\201607\fnl_20160704_00_00.gmf'
'set lat 20 55 '
'set lon 85 135'
'set lev 850'
'set grads off'
'set grid off'
'set mpdset cnworld'
#'set map  15 5 1'
#'draw map'
#'set cthick 9'
#'set xlopts 1 4 0.12'
#'set ylopts 1 4 0.12'
#'set cint 0.5'

'define prs=lev'
'define es=6.112*exp(17.67*(TMPprs-273.15)/(TMPprs-29.65))'
'define qs=0.622*es/(lev-0.278*es)'
'define q=RHprs*qs/100'
'define a=mag(UGRDprs,VGRDprs)'
'define g=9.8'
'set gxout shaded'
'set cmin 10'
'd a*q*1000/9.8'
'run cbarn.gs'
'set gxout contour'
'set ccolor 0'
'set cmin 6'
'd a*q*1000/9.8'
'print'
pull dummy
*t=t+1
*endwhile
'disable print'

Ncl脚本：
begin
files=systemfunc("ls *.grib2")
  a=addfile(files(3),"r")
  time=03
  U=a->UGRD_P0_L100_GLL0
  V=a->VGRD_P0_L100_GLL0
  RH=a->RH_P0_L100_GLL0
  TK=a->TMP_P0_L100_GLL0
  T=TK-273.15
  ;lt=a->lv_ISBL0
  ;lr=a->lv_ISBL4
  ;print(lt)
  ;print(lr)
  ;lat=a->lat_0
  ;lon=a->lon_0
  es=6.112*exp((17.67*T)/(T+243.5))
  qs=0.622*es(25,:,:)/(850-0.278*es(25,:,:))
  q=RH(25,:,:)*qs/100
  v=sqrt(U(25,:,:)^2+V(25,:,:)^2)
  Q=v*q*1000/9.8
  Q@units="g/kg"
  lat=fspan(90,-90,181)
  lon=fspan(0,359,360)
  lat@units="degrees_north"
  lon@units="degrees_east"
  RH!0="lat"
  ;q&lat=lat
  RH!1="lon"
  ;q&lon=lon
  wks=gsn_open_wks("x11","Qvapor_flux" + time)
  ;plot=new(6,graphic)
  gsn_define_colormap(wks,"WhiteBlue")
  ;gsn_define_colormap(wks,"BlueRedGray")
  res=True
  res@gsnFrame=False
  res@gsnDraw=False
  ;res@PanelPlot=True
res@cnLinesOn=False
  res@cnLineLabelsOn=False
  res@cnFillOn=True
  res@cnLevelSelectionMode="ExplicitLevels"
  ;res@cnFillColors=(/18,34,50,66,82,98,114,130,146,162,178,194,210,226/)
  res@cnLevels=(/6,8,10,12,14,16,18,20,22,24,26/)
  
  res@gsnAddCyclic=False
  res@tfDoNDCOverlay= True
  res@mpAreaMaskingOn = True
  ;res@mpMaskAreaSpecifiers = (/"States"/)
  res@mpDataBaseVersion = "MediumRes"
  res@mpDataSetName = "Earth..4"
  res@mpOutlineOn = True
  res@mpOutlineSpecifiers = (/"china","jiangsu","anhui","shandong","henan","zhejiang",\
                              "jiangxi","hubei","fujian","shanxi","guangxi",\
                              "hunan","shanxi","guizhou","sichuan","gansu"/)
  res@mpProvincialLineThicknessF = 1
  res@mpProvincialLineColor="grey39"
  res@mpGeophysicalLineThicknessF  = 1
  res@mpGeophysicalLineColor="grey39"
  res@mpLimitMode= "Corners"
  res@mpLeftCornerLatF  = lat(70)
  res@mpLeftCornerLonF  = lon(85)
  res@mpRightCornerLatF = lat(35)
  res@mpRightCornerLonF = lon(135)
  
  res@lbLabelBarOn=True
  res@tmYROn=False
  res@tmXTOn=False
  res@tmXBOn=True
  res@tmXBLabelsOn=True
  res@tmXBMode="Explicit"
  res@tmXBValues=(/85,90,95,100,105,110,115,120,125,130,135/)
  res@tmXBLabels=(/"85E","90E","95E","100E","105E","110E","115E","120E","125E","130E","135E"/)
  res@tmYLOn=True
  res@tmYLLabelsOn=True
  res@tmYLMode="Explicit"
  res@tmYLValues=(/20,25,30,35,40,45,50,55/)
  res@tmYLLabels=(/"20N","25N","30N","35N","40N","45N","50N","55N"/)

  plot=gsn_csm_contour_map(wks,Q,res)

draw(wks)
  frame(wks)
end
附图

yl8809357

作为一个小白  完全没看懂这些代码是什么意思

Lancelot

RH!0="lat"
  ;q&lat=lat
  RH!1="lon"
  ;q&lon=lon
因为这里被注释掉了

wow哈哈

谢谢你的回复，不过好像不是这个问题，这个注释去掉就报错了

Lancelot

wow哈哈 发表于 2017-5-22 11:23
谢谢你的回复，不过好像不是这个问题，这个注释去掉就报错了

不，就是这里的问题，改程序的时候你搞错了没有设置Q的经纬属性

Eudora99

qs=0.622*es(25,:,:)/(850-0.278*es(25,:,:))这一步应该为 qs=0.622*es(25,:,:)/(850-0.378*es(25,:,:))吧?是0.378，而不是0.278.

柳絮

请问楼主问题解决了吗？

wow哈哈

柳絮 发表于 2018-4-26 23:40
请问楼主问题解决了吗？

解决了，就是plot里面没有给Q的范围，地图和数据大小未匹配

柳絮

wow哈哈 发表于 2018-4-27 09:05
解决了，就是plot里面没有给Q的范围，地图和数据大小未匹配

哦 好的 谢谢 能把脚本分享一下吗？

wow哈哈

柳絮 发表于 2018-4-27 10:08
哦 好的 谢谢 能把脚本分享一下吗？

太久之前的问题了，现在脚本已经改的面目全非了

柳絮

wow哈哈 发表于 2018-5-2 21:20
太久之前的问题了，现在脚本已经改的面目全非了

好的 谢谢 我已画出来了