;***********Load the ncl internal libraries*********; load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" begin ;***data*** pre_data = addfile("SRF_IST_2009090100-2009091200_tpr_daily.nc","r") topo_data = addfile("DOMAIN_IST.nc","r") yagis_all = asciiread("yagis_station.txt",(/28,7/),"float") ;**ERA lat lon*** latt1 = 28 latt2= 76 lonn1 = 50 lonn2 = 94 lat2d = topo_data->xlat(latt1:latt2,lonn1:lonn2) lon2d = topo_data->xlon(latt1:latt2,lonn1:lonn2) ;***pre sum model*** pre = pre_data->tpr(6,latt1:latt2,lonn1:lonn2) pre1 = pre_data->tpr(7,latt1:latt2,lonn1:lonn2) pre2 = pre_data->tpr(8,latt1:latt2,lonn1:lonn2) pre_sum = pre+pre1+pre2 ;***sta sum *** yagis_07 = yagis_all(:,3) yagis_08 = yagis_all(:,4) yagis_09 = yagis_all(:,5) yagis_10 = yagis_all(:,6) yagis_sum = yagis_07+yagis_08+yagis_09 lon = yagis_all(:,0) lat = yagis_all(:,1) ;***get dimensions*** dimvar= dimsizes(pre) jlat = dimvar(0) ilon = dimvar(1) ;***markers*** arr= (/5,10,15,20,30,40,50,60,80,100,120/);precip_11lev colors = (/2,3,4,5,6,7,8,9,10,11,12,13/);precip_11lev dum = new(dimsizes(colors),graphic) labels = new(dimsizes(arr)+1,string) num_distinct_markers = dimsizes(arr)+1 lat_new = new((/num_distinct_markers,28/),float,-999) lon_new = new((/num_distinct_markers,28/),float,-999) staval=yagis_08 do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(staval.lt.arr(0)) labels(i) = "x < " + arr(0) end if if (i.eq.num_distinct_markers-1) then indexes = ind(staval.ge.max(arr)) labels(i) = "x >= " + max(arr) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(staval.ge.arr(i-1).and.staval.lt.arr(i)) labels(i) = arr(i-1) + " <= x < " + arr(i) end if if (.not.any(ismissing(indexes))) then npts_range = dimsizes(indexes) lat_new(i,0:npts_range-1) = lat(indexes) lon_new(i,0:npts_range-1) = lon(indexes) end if delete(indexes) end do wks = gsn_open_wks("ps","stations_pre_7_8_9_colormap") gsn_define_colormap (wks,"precip_11lev") ;***colormap resources*** resmap = True resmap@cnFillOn = True resmap@cnLinesOn = False resmap@cnLineLabelsOn = False resmap@cnInfoLabelOn = False resmap@lbLabelBarOn = False resmap@mpFillOn = False resmap@mpLimitMode = "Corners" resmap@mpLeftCornerLatF = lat2d(0,0) resmap@mpLeftCornerLonF = lon2d(0,0) resmap@mpRightCornerLatF = lat2d(jlat-1,ilon-1) resmap@mpRightCornerLonF = lon2d(jlat-1,ilon-1) resmap@pmTickMarkDisplayMode = "Always" resmap@tmXBLabelFontHeightF = 0.012 resmap@tmYLLabelFontHeightF = 0.012 resmap@gsnCenterStringFontHeightF = 0.028 resmap@gsnMaximize = True resmap@tmXTOn = False resmap@tmYROn = False resmap@tiMainString = "" resmap@tiMainFontHeightF = 0.018 resmap@gsnDraw = False resmap@gsnFrame = False resmap@mpProjection = "LambertConformal" resmap@mpLambertParallel1F = 30. resmap@mpLambertParallel2F = 60. resmap@mpLambertMeridianF =33 resmap@cnLevelSelectionMode = "ExplicitLevels" resmap@cnLevels = (/5,10,15,20,30,40,50,60,80,100,120/) resmap@gsnScale = True resmap@gsnAddCyclic = False resmap@tfDoNDCOverlay = True resmap@mpDataBaseVersion = "HighRes" plot = gsn_csm_contour_map(wks,pre_sum,resmap) ;***marker resources*** polyres = True polyres@gsMarkerIndex = 16 polyres@gsMarkerSizeF = 5. do i = 0, num_distinct_markers-1 if (.not.ismissing(lat_new(i,0))) polyres@gsMarkerColor = colors(i) polyres@gsMarkerThicknessF = 2.8*(i+1) dum(i)=gsn_add_polymarker(wks,plot,lon_new(i,:),lat_new(i,:),polyres) end if end do ;***panel resources*** resP = True resP@gsnPanelLabelBar = True resP@gsnFrame = False resP@lbOrientation = "Horizontal" resP@lbLabelAutoStride = True resP@pmLabelBarWidthF = 0.5 resP@pmLabelBarHeightF = 0.05 resP@lbTitleOn = True resP@lbTitleString = "mm" resP@lbTitlePosition = "Right" resP@lbTitleFontHeightF= .012 resP@lbTitleDirection = "Across" resP@lbLabelStride = 2 resP@pmLabelBarParallelPosF = 0.02 resP@lbLabelFontHeightF = 0.014 resP@pmLabelBarOrthogonalPosF = -0.05 resP@gsnStringsFontHeightF = .002 gsn_panel(wks,plot,(/1,1/),resP) frame(wks) end