- 积分
- 847
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-5-10
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2014-10-20 17:36:01
|
显示全部楼层
是参考了之前一位同学发布的水汽通量图调用fortran的脚本:
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/wrf/WRFUserARW.ncl"
external MY "/mnt/hgfs/E/term/20110723/vaporflux.so"
begin
;========================
; get list of all files and open as "one big file"
;========================
fall0302 = addfile("wrfout_d01_2011-07-23_00_00_00.nc","r") ; note the "s" of addfile
;========================
; choose how files are combined and read in variable across files
;========================
; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"vaporflux")
;===============================================================
; What times and how many time steps are in the data set?
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
ntimes = dimsizes(times) ; number of times in the file
print(ntimes)
; The specific pressure levels that we want the data interpolated to.
pressure_levels = (/850./) ; pressure levels to plot
nlevels = dimsizes(pressure_levels) ; number of pressure levels
;==========================================================
do it = 0,ntimes-1,1 ; TIME LOOP
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
;=================================================================
; First get the variables we will need
tc = wrf_user_getvar(fall0302,"tc",it) ; T in C
u = wrf_user_getvar(fall0302,"ua",it) ; u averaged to mass points
v = wrf_user_getvar(fall0302,"va",it) ; v averaged to mass points
p = wrf_user_getvar(fall0302, "pressure",it) ; pressure is our vertical coordinate
z = wrf_user_getvar(fall0302, "z",it) ; grid point height
rh = wrf_user_getvar(fall0302,"rh",it) ; relative humidity
;##########################################################################################
radius = 20
g = 9.8
;*******************************************************************************************************
qu =new((/12,41,41/),float)
qv =new((/12,41,41/),float)
mm =new((/12,41,41/),float)
;*******************************************************************************************
xx=41
yy=41
zz=1
tt=12
MY::vaporflux(p,t,rh,u,v,qu,qv,mm,xx,yy,zz,tt,g)
printVarSummary(qu)
copy_VarMeta(u,qu)
copy_VarMeta(v,qv)
copy_VarMeta(rh,mm)
printVarSummary(qu)
;=====================================================================================================
;plot
;=====================================================================================================
res = True
res@gsnAddCyclic = False
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@vcLineArrowThicknessF = 2.0
res@vcRefAnnoOrthogonalPosF = 0 ;ref的位置
res@vcRefAnnoString1On = True
res@vcRefAnnoString2On = False
res@vcRefMagnitudeF = 20 ; make vectors larger
res@vcRefLengthF = 0.1 ; ref vec length
res@vcGlyphStyle = "CurlyVector" ; turn on curly vectors
res@vcMinDistanceF =0.02 ;箭头疏密
res@vcVectorDrawOrder= "PostDraw"
res@gsnLeftString =""
res@gsnRightString =""
res@tiMainString =""
cnres = True ; plot mods desired
cnres@gsnDraw = False ; don't draw yet
cnres@gsnFrame = False ; don't advance frame yet
cnres@cnFillOn = True ; turn on color
cnres@gsnSpreadColors = True ; use full colormap
cnres@gsnSpreadColorStart = 0
cnres@gsnSpreadColorEnd = 201
cnres@cnLinesOn = False ; turn off contour lines
cnres@cnLineLabelsOn = False ; tuen off line labels
cnres@cnInfoLabelOn =False
cnres@cnLevelSelectionMode = "AutomaticLevels"
;cnres@cnMinLevelValF = 4
;cnres@cnMaxLevelValF = 42
;cnres@cnLevelSpacingF =2
cnres@gsnLeftString =""
cnres@gsnRightString =" g s~S~ -1~N~cm~S~ -1~N~hPa~S~ -1"
cnres@tiMainString =""
cnres@lbLabelBarOn = False
mm=smth9(mm,0.5,-0.25,False)
printMinMax (mm*1000, True)
; dir = "/cma/u/Twangdh/hxe/example/"
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
print("====plot..........===========")
plots = new((/2,12/),graphic)
plots(0,5) =wrf_vector(fall0302,wks,qu(5,:,:)*1000,qv(5,:,:)*1000,res) ; create plot
plots(1,5) = wrf_contour(fall0302,wks,mm(5,:,:)*1000,cnres)
overlay(plots(0,5),plots(1,5))
draw(plots)
frame(wks)
;#########################################################################
end |
|