- 积分
- 816
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-6-18
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2017-10-14 17:13:34
|
显示全部楼层
本帖最后由 dengxiaohua_121 于 2017-10-14 17:15 编辑
出错的提示就是这样子的。下面是我运行的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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl"
begin
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
ftime=2017011000
start_lat=-15
end_lat=30
start_lon=30
end_lon=100
middle_lon=(start_lon+end_lon)/2.0
path_parent0="/fs01/home/climate/INTERP/data/region_atm/"
path_parent=path_parent0+tostring(ftime)
path_parent1="/fs01/home/climate/INTERP/qixiang/severe_weather"
f_all_file_u=systemfunc("ls "+path_parent0+tostring(ftime)+"/HTTTG*.grb")
nFileList_file_u=dimsizes(f_all_file_u)
f_all_file_v=systemfunc("ls "+path_parent0+tostring(ftime)+"/HTRHG*.grb")
nFileList_file_v=dimsizes(f_all_file_v)
counter_file_tmp=(/nFileList_file_u,nFileList_file_v/)
counter_file=min(counter_file_tmp)
length_str_u=strlen(f_all_file_u(0))
length_str_v=strlen(f_all_file_v(0))
dname_u=(/length_str_u-7,3,4/)
do ii=1,2;counter_file-1;第一个00时刻的不算
output_name_u=str_split_by_length(f_all_file_u(ii),dname_u)
if(stringtodouble(output_name_u(1)) .le. 120.0) then
test_str_input=output_name_u(1)
currentFilePath_u=output_name_u(0)+test_str_input+output_name_u(2)
output_name_v=str_split_by_length(f_all_file_v(0),dname_u)
currentFilePath_v=output_name_v(0)+test_str_input+output_name_v(2)
currentFile_u = addfile(currentFilePath_u, "r")
currentFile_v = addfile(currentFilePath_v, "r")
varNames_u = getfilevarnames (currentFile_u)
varNames_v = getfilevarnames (currentFile_v)
length_str=strlen(currentFilePath_u)
length_str1=strlen(path_parent)
lat=currentFile_u->lat_3
lon=currentFile_u->lon_3
lv_ISBL0=currentFile_u->lv_ISBL0
u_10_tmp1=currentFile_u->TMP_3_ISBL
v_10_tmp1=currentFile_v->R_H_3_ISBL
latitude_new=lat({start_lat:end_lat})
longitude_new=lon({start_lon:end_lon})
counter_lat=dimsizes(latitude_new)
counter_lon=dimsizes(longitude_new)
u_10_tmp2=u_10_tmp1(:,{start_lat:end_lat},{start_lon:end_lon})-273.15
v_10_tmp2=v_10_tmp1(:,{start_lat:end_lat},{start_lon:end_lon})
data_v2=v_10_tmp2
data_v2=(v_10_tmp2-60)*2.5*(u_10_tmp2*(u_10_tmp2+16)/(-64.0))
ic_index_v2 = where(data_v2.ge.5.0 .and. data_v2 .lt. 50.0,1.0, 0.0)
ic_index_v2 = where(data_v2.ge.50.0 .and. data_v2 .lt. 80.0,2.0, ic_index_v2)
ic_index_v2 = where(data_v2.ge.80.0,3.0, ic_index_v2)
;again 6+
do k=0,4;5
if(k .eq. 0)then
kk=24
str_part="975_hPa"
end if
if(k.eq.1)then
kk=20
str_part="850_hPa"
end if
if(k.eq.2)then
kk=17
str_part="700_hPa"
end if
if(k.eq.3)then
kk=13
str_part="500_hPa"
end if
if(k.eq.4)then
kk=9
str_part="300_hPa"
end if
ic_index_v3=reshape(ic_index_v2(kk:kk,:,:),(/counter_lat,counter_lon/))
ic_index_v3!0="lat"
ic_index_v3!1="lon"
ic_index_v3&lat=latitude_new
ic_index_v3&lon=longitude_new
index_aircraft_icing=ic_index_v3
dname=(/length_str1+1,5,10,1,3,4/)
output_name=str_split_by_length(currentFilePath_u,dname)
infor_time=output_name(4)
ti_string=output_name(2)
range_only = ic_index_v3 ; trick to keep cv's and atts
range_only = mask(ic_index_v3,(ic_index_v3.ge.1.0),True)
;绘图
wks = gsn_open_wks ("png", path_parent1+"/Fig/"+"aircraft_icing-fst_"+str_part+"_"+ti_string+"_"+infor_time+"_NI" )
cmap = RGBtoCmap(path_parent1+"/rgb2.txt")
gsn_define_colormap(wks,cmap)
gsn_reverse_colormap(wks) ; Reverse the color map.
res = True ; plot mods desired
res@gsnRightString = str_part+"_"+ti_string+"_"+infor_time ;units
res@gsnLeftString = "aircraft_icing" ;long_name
res@mpDataSetName = "Earth..3" ; This new database contains
res@mpDataBaseVersion = "MediumRes" ; Medium resolution database
res@mpOutlineOn = True ; Turn on map outlines
res@mpProjection = "CylindricalEquidistant"
res@mpGeophysicalLineThicknessF = 2.0 ; double the thickness of geophysical boundaries
res@mpNationalLineThicknessF = 3.0 ; double the thickness of national boundaries
res@cnSmoothingOn=True
res@mpMinLatF = start_lat;-90;-15 ; range to zoom in on
res@mpMaxLatF = end_lat;50
res@mpMinLonF = start_lon;30
res@mpMaxLonF = end_lon;170
res@mpCenterLonF =middle_lon; 120.0
res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res@cnMinLevelValF = 0; ; set min contour level
res@cnMaxLevelValF = 3; ; set max contour level
res@cnLevelSpacingF =1 ; set contour spacing
res@cnFillOn = True ; color plot desired 打开等值线填色
res@cnLinesOn = False ; turn off contour lines 关闭等值线线条
res@lbLabelsOn = True
;res@lbLabelBarOn = False
res@cnLineLabelsOn = False ; turn off contour labels
res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks打开坐标
res@mpLandFillColor = "palegreen4"
res@mpOceanFillColor = "midnightblue";"dodgerblue3"
res@mpFillDrawOrder="PreDraw"
;res@pmLabelBarOrthogonalPosF=0.10
;res@lbLabelStride = 1
res@pmLabelBarWidthF = 0.1
res@pmLabelBarHeightF = 0.45
res@lbOrientation = "Vertical"
res@lbLabelFontHeightF = 0.015
res@gsnAddCyclic = False ;False
res@gsnFrame = False
;res@gsnRightString =ti_string+"_"+infor_time
res@cnSmoothingOn = True
res@cnSmoothingDistanceF = 0.001
res@cnSmoothingTensionF = -0.1
res@cnLineThicknessF =2.0
range_only=smth9(range_only, 0.50, 0, False)
;plot=gsn_csm_contour_map(wks,u_10_tmp2,res)
plot=gsn_csm_contour_map(wks,range_only,res)
draw(plot)
ntxres = True
ntxres@txFontHeightF = 0.018
ntxres@txFontColor = "Black"
gsn_text_ndc(wks,"Aircraft_icing Category: 1-mild; 2-moderate; 3-severe;",0.42,0.22,ntxres)
frame(wks)
delete([/plot,wks,ti_string,output_name,ic_index_v3,dname,range_only/])
end do
end if
end do
end
|
|