爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 18709|回复: 26

[作图] ncl画整层水汽通量求助

[复制链接]

新浪微博达人勋

发表于 2017-4-2 00:58:42 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
参考脚本:
http://bbs.06climate.com/forum.php?mod=viewthread&tid=28890
http://bbs.06climate.com/forum.php?mod=viewthread&tid=35437
资料是FNl的grib2,调了一晚上都没成功,跪求大神帮忙,红色是报错
oad "$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=1   
     nlev=8         
      
grb_file = addfile("/home/fuchao/datafile/FNL/fnl_20140809_00_00.grib2","r")
  tt       = grb_file->TMP_P0_L102_GLL0(:,:,:)
  q    = grb_file->SPFH_P0_2L108_GLL0(:,:)
  u        = grb_file->UGRD_P0_L100_GLL0(:,:,:)
  v        = grb_file->VGRD_P0_L100_GLL0(:,:,:)
slp       = grb_file->PRES_P0_L1_GLL0(:,:)   
qu1 = new((/8,73,144/),float)
        qv1 = new((/8,73,144/),float)
        do iz=0,7
                do i=0,72
                        do j=0,143
                                qu1(iz,i,j) = q(i,j)*u(iz,i,j)
                                qv1(iz,i,j) = q(i,j)*v(iz,i,j)
                        end do
                end do
        end do
   lev =(/1000,925,850,700,600,500,400,300/)
    lat=fspan(-90,90,73)
    lon=fspan(0,357.5,144)
        
        lon!0                 = "lon"
        lon@units         = "degrees_east"
        lat!0                 = "lat"
        lat@units         = "degrees_north"
        lev!0                = "lev"
        lev@units        = "hPa"
        
        qu1!0         = "lev"
        qu1!1         = "lat"
        qu1!2         = "lon"
        qu1&lon = lon
        qu1&lat = lat
        qu1&lev = lev
      
        
        qv1!0         = "lev"
        qv1!1         = "lat"
        qv1!2         = "lon"
        qv1&lon = lon
        qv1&lat = lat
        qv1&lev = lev
;**************************************************************************************
        qu11 = qu1(lat|:, lon|:, lev|:)
        delete(qu1)
        qv11 = qv1(lat|:, lon|:, lev|:)
        delete(qv1)
         p       = (/ 1000.,925.,850.,700.,600.,500., \
               400.,300.,250.,200.,150.,100., \
                70.,50.,30.,20.,10. /)
         linlog = 2
     pbot   = 1100.                                             
     ptop   = 300.
        linlog = 2                                                   
                                                      ; 1为线性积分  2为对数积分
        data11 = (vibeta(p,qu11,linlog,slp,1100,300))/9.8           ; 返回的数值是2维的(lat,lon)
        data12 = (vibeta(p,qv11,linlog,slp,1100,300))/9.8
      
        qu  = new((/73,144/),float)
        qv  = new((/73,144/),float)
        mm  = new((/73,144/),float)
  mm1 = sqrt(data11^2+data12^2)

        do i=0,72
                do j=0,143
                        mm(i,j) = mm1(i,j)
                end do
        end do

        do i=0,72
                do j=0,143
                        qu(i,j) = data11(i,j)
                end do
        end do

        do i=0,72
                do j=0,143
                        qv(i,j) = data12(i,j)
                end do
        end do        

        mm!0 = "lat"
        mm!1 = "lon"
        mm&lat = lat
        mm&lon = lon
        lon!0 = "lon"
        lon@units = "degrees_east"
        lat!0 = "lat"
        lat@units = "degrees_north"        
        qu!0 = "lat"
        qu!1 = "lon"
        qu&lat = lat
        qu&lon = lon
        lon!0 = "lon"
        lon@units = "degrees_east"
        lat!0 = "lat"
        lat@units = "degrees_north"
        qv!0 = "lat"
        qv!1 = "lon"
        qv&lat = lat
        qv&lon = lon
        lon!0 = "lon"
        lon@units = "degrees_east"
        lat!0 = "lat"
        lat@untis = "degrees_north"

wks =gsn_open_wks("pdf","fnl_vapour_flux")
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200")   
  vcres = True
  vcres@gsnAddCyclic  = False
  vcres@gsnDraw      = False                        ; don't draw yet
  vcres@gsnFrame     = False                        ; don't advance frame yet

vcres@mpDataSetName              = "Earth..4"
  vcres@mpDataBaseVersion          = "MediumRes" ; or "Ncarg4_1"
  vcres@mpAreaMaskingOn            = True
  vcres@mpMaskAreaSpecifiers       = (/"China"/)
  vcres@mpOutlineSpecifiers        = (/"China","China:Provinces"/)
   
  vcres@mpLandFillColor            = "white"
  vcres@mpInlandWaterFillColor     = "white"
  vcres@mpOceanFillColor           = "white"
  vcres@mpFillBoundarySets         = "NoBoundaries"
  vcres@mpOutlineBoundarySets      = "NoBoundaries"
  vcres@mpNationalLineColor        = "black"
  vcres@mpProvincialLineColor      = "black"
  vcres@mpGeophysicalLineColor     = "black"
  vcres@mpNationalLineThicknessF   = 2
  vcres@mpProvincialLineThicknessF = 1
  vcres@mpProjection          = "CylindricalEquidistant"
  vcres@tiMainString          =""

  vcres@mpGeophysicalLineThicknessF = 0.35      ; double the thickness of geophysical boundaries
  vcres@mpNationalLineThicknessF    = 0.5      ; double the thickness of national boundaries

  vcres@pmTickMarkDisplayMode = "Always"
          vcres@mpMinLatF         =  0                         ;最小纬度
          vcres@mpMaxLatF         =  50                        ;最大纬度
          vcres@mpMinLonF         =  60                        ;最小经度
          vcres@mpMaxLonF         =  160                       ;最大经度


   vcres@gsnDraw      = False                        ; don't draw yet
  vcres@gsnFrame     = False                        ; don't advance frame yet
                                       

  vcres@lbLabelBarOn = True
  vcres@lbLabelStride       = 2         ; plot every other colar bar label
  vcres@lbOrientation        = "vertical"         ; vertical label bars

  vcres@vcRefMagnitudeF         = 10.0              ; make vectors larger
  vcres@vcRefLengthF            = 0.018            ; ref vec length
  vcres@vcGlyphStyle            = "CurlyVector"    ; turn on curly vectors
  vcres@vcMinDistanceF          = 0.017            ; thin out vectors
  vcres@vcRefAnnoOn = False


  vcres@vcMonoFillArrowFillColor  = True
  vcres@vcLineArrowThicknessF  = 1.0
  vcres@vcRefAnnoOn  =False
  vcres@vcVectorDrawOrder= "PostDraw"
    vcres@gsnLeftString =""
   vcres@gsnRightString =""
  vcres@tiMainString =""
  plot(0)=gsn_csm_vector_map(wks,qu,qv,vcres)  ; create plot
  txres               = True                           
  txres@txFontHeightF = 0.020             ; Set the font height
  txres@txPerimOn     = True
  txres@txBackgroundFillColor =0
  txres@txPerimColor    = 0
;************************************************************************************************
res                 = True                    ; plot mods desired
res@gsnDraw         = False                   ; don't draw yet
res@gsnFrame        = False                   ; don't advance frame yet
res@cnFillOn        = True                    ; turn on color
res@gsnSpreadColors = True                    ; use full colormap
res@gsnSpreadColorStart = 0
res@gsnSpreadColorEnd   = 199
res@cnLinesOn       = False                   ; turn off contour lines
res@cnLineLabelsOn  = False                   ; tuen off line labels
res@cnInfoLabelOn =False

res@cnLevelSelectionMode="ExplicitLevels"
res@cnLevels=(/14,18,22,26,30,34,38,42,46,50,52/)

  res@gsnLeftString =""
  res@gsnRightString =""
  res@tiMainString =""

res@lbLabelBarOn = False
res@lbLabelStride       = 2         ; plot every other colar bar label
res@lbOrientation        = "vertical"         ; vertical label bars

   mm=smth9(mm,0.5,0.25,False)

plot01 = gsn_csm_contour(wks,mm,res)
overlay(plot(0),plot01)
end

2017-04-02 01-03-29 的屏幕截图.png

fnl_20140809_00_00.grib2

12.9 MB, 下载次数: 28, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-3 08:56:08 | 显示全部楼层
slp变量是(lev,lat,,lon)和qu11 = qu1(lat|:, lon|:, lev|:)没有对应上
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-4-3 09:03:40 | 显示全部楼层
牧场物语之主 发表于 2017-4-3 08:56
slp变量是(lev,lat,,lon)和qu11 = qu1(lat|:, lon|:, lev|:)没有对应上

你好,谢谢你的回复,不是那里的问题,是格点数没对上,资料是1度的。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-6 23:57:10 | 显示全部楼层
楼主做出来了吗?可不可以把完整的脚本分享一下?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-9-10 23:11:54 | 显示全部楼层
脚本能分享下吗??谢谢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-9-13 18:45:36 | 显示全部楼层
ccqx2007 发表于 2017-9-10 23:11
脚本能分享下吗??谢谢

   nlat=73
     nlon=144
,这里改成数据里的分辨率就行了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-10-15 06:58:58 | 显示全部楼层
请问你的q    = grb_file->SPFH_P0_2L108_GLL0(:,:)只有两维,我印象中好像是地面的比湿吧?整层全用这个算是不是不太合理啊?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-10-15 08:21:28 | 显示全部楼层
叉一由八 发表于 2017-10-15 06:58
请问你的q    = grb_file->SPFH_P0_2L108_GLL0(:,:)只有两维,我印象中好像是地面的比湿吧?整层全用这个 ...

谢谢提醒,还是用相对湿度去计算要好
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-10-15 11:05:33 | 显示全部楼层
1649518749 发表于 2017-10-15 08:21
谢谢提醒,还是用相对湿度去计算要好

恩恩!另外谢谢你的脚本~我也是用这个资料,你的脚本给我指了好多条明路!谢谢~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-10-15 12:20:56 | 显示全部楼层
谢谢楼主!!!作为NCL小白很有启示
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表