; Example script to produce plots for a WRF real-data run, ; with the ARW coordinate dynamics option. ; Plot data on a cross section ; This script will plot data at a set angle through a specified point ; Add some label info to the Y-axis 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/wrf/WRFUserARW.ncl" load "/home/xhu2/tools/nclscripts/wrf/WRFUserARW.ncl" begin gsres = True gsres@gsMarkerIndex = 16 ; circle at first gsres@gsMarkerThicknessF =3 gsres@gsMarkerSizeF = 0.02 gsres@gsMarkerColor = "black" ; gsres@gsMarkerColor = "red" ; gsn_polymarker(wks,plot,-96.8,32.7833333,gsres) ; gsn_text(wks,plot,"D",-96.8,32.7833333+0.1,tres) GapWind_Lat = 30.5319; 32.7833333 GapWind_Lon = 117.1151-0.2-0.35; Anqing°. idomain = 2 ; wrfout_d02_2016-12-06_12:00:00 files = systemfunc("ls wrfout_d0"+idomain+"_2016-12-06_12\:00\:00") a = addfile(files(0)+".nc","r") opt = True ij = wrf_user_ll_to_ij(a,GapWind_Lon, GapWind_Lat,opt) ;do ihour = 0 ,84; dimsizes(files)-1 print(ij) ; Set some basic resources res = True ; res@MainTitle = "REAL-TIME WRF" res@Footer = False ter_res = True opts_ter = ter_res opts_ter@gsnYRefLine = 0.0 opts_ter@gsnAboveYRefLineColor = "transparent" opts_ter@gsnDraw = False opts_ter@gsnFrame = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; times = wrf_user_getvar(a,"times",-1) ; get times in the file ntimes = dimsizes(times) ; number of times in the file FirstTime = True mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file nd = dimsizes(mdims) ;--------------------------------------------------------------- do it = 60,69; ntimes-1 ; TIME LOOP ; We generate plots, but what kind do we prefer? type = "png" type@wkWidth = 800 type@wkHeight = 800 ; figurename = "wrfout_d0"+idomain+"_WSP_crossSNpartial_"+it figurename = "wrfout_d0"+idomain+"_WSP_crossSN16partial_"+it wks = gsn_open_wks(type,figurename) gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map cmap = gsn_retrieve_colormap(wks) do icorlor = 2,4 cmap(icorlor,:) = cmap(0,:) end do gsn_define_colormap(wks,cmap) ; select color map print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots ter = wrf_user_getvar(a,"ter",0) xlon = wrf_user_getvar(a, "XLONG",0) xlat = wrf_user_getvar(a, "XLAT",0) tc = wrf_user_getvar(a,"tc",it) ; T in C theta = wrf_user_getvar(a,"theta",it) ; relative humidity z = wrf_user_getvar(a, "z",it) ; grid point height va = wrf_user_getvar(a,"va",it) ; u on mass points ua = wrf_user_getvar(a,"ua",it) ; u on mass points WSP = sqrt(ua^2+va^2) wa = wrf_user_getvar(a,"wa",it) ; v on mass points pressure = wrf_user_getvar(a,"pressure",it) ; v on mass points if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = max(z)/1000. nz = floattoint(zmax/2 + 1) FirstTime = False ; print(nz) end if z_values = fspan(zmin,zmax,nz) ;printVarSummary(z_values) dimzz = dimsizes(z_values) ;--------------------------------------------------------------- ; Plot a cross session that run south-north through the middle of the plot ; For this we need a pivot point and a angle ; | ; angle=0 is | ; angle = 90 ; E-W angle = 29 ; N-S plane = new(2,float) plane = (/ ij(0), ij(1) /) ; pivot point is center of domain (x,y) opts = False ; start and end points not specified theta_plane = wrf_user_intrp3d(theta,z,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts) WSP_plane = wrf_user_intrp3d(WSP,z,"v",plane,angle,opts) va_plane = wrf_user_intrp3d(va,z,"v",plane,angle,opts) wa_plane = wrf_user_intrp3d(wa,z,"v",plane,angle,opts) pressure_plane = wrf_user_intrp3d(pressure,z,"v",plane,angle,opts) printVarSummary(pressure_plane) print("max="+max(pressure_plane)+" min="+min(pressure_plane)) ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_plane = wrf_user_intrp2d(xlat,plane,angle,opts) distance = (lon_plane( :)-GapWind_Lon)^2+(X_plane(:)-GapWind_Lat)^2 id_d = minind(distance) dim = dimsizes(theta_plane) ; Find the data span - for use in labels zspan = dim(0) zcordinate = fspan(zmin,zmax,zspan) ; print(zspan) ; Options for XY Plots opts_xy = res opts_xy@tiYAxisString = "Height (km)" opts_xy@cnMissingValPerimOn = False ; True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 11 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(0,zspan,nz*2*2) ; Create tick marks opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz*2*2)) ; Create labels opts_xy@tiXAxisFontHeightF = 0.020 opts_xy@tiYAxisFontHeightF = 0.020 opts_xy@tmXBMajorLengthF = 0.02 opts_xy@tmYLMajorLengthF = 0.02 opts_xy@tmYLLabelFontHeightF = 0.015 opts_xy@PlotOrientation = tc_plane@Orientation n_YMax = 50 ; n_XMax = dim(1)-2; 70 n_XMin = 0 n_XMax = 160 n_XMin = 60 opts_xy@trXMaxF = n_XMax opts_xy@trXMinF = n_XMin opts_xy@trYMaxF = n_YMax ; opts_xy@gsnOrientation = "portrait" opts_xy@gsnDraw = True ; Forces the plot to be drawn opts_xy@gsnFrame = True ; Frame advance opts_xy@tiMainString = "" ;chartostring(a->Times) opts_xy@vpHeightF = 0.35 opts_xy@vpWidthF = 0.5 ; opts_xy@gsnRightString = chartostring(a->Times) ; Plotting options for RH opts_theta = opts_xy opts_theta@cnFillMode = "RasterFill" opts_theta@ContourParameters = (/ 4., 14., .5 /) opts_theta@tmXBFormat = "f" opts_theta@pmLabelBarOrthogonalPosF = -0.07 opts_theta@cnFillOn = True opts_theta@gsnSpreadColors = True ; use full range of colormap ; opts_theta@cnFillColors = (/"White","White","White", \ ; "White","Chartreuse","Green", \ ; "Green3","Green4", \ ; "ForestGreen","PaleGreen4"/) ; Plotting options for Temperature opts_tc = opts_xy opts_tc@cnInfoLabelOrthogonalPosF = 0.00 opts_tc@ContourParameters = (/ 5. /) opts_theta@lbTopMarginF = -0.05810482 opts_theta@lbBottomMarginF = 0.21 opts_theta@lbLeftMarginF = 0.551 ; Get the contour info for the theta and temp ; printVarSummary(tc_plane) ; contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) opts_theta@gsnMaximize = False opts_theta@cnLineLabelsOn = False ;True ; print(opts_theta) contour_theta = wrf_contour(a,wks,WSP_plane,opts_theta) opts_theta@ContourParameters = (/ 990., 1026., 1. /) ; delete(opts_theta@ContourParameters ); = (/ .5 /) ; opts_theta@ContourParameters = (/ .5 /) opts_theta@cnFillOn = False ; True opts_theta@cnMonoLineColor = False ; Tells NCL not to draw contour lines in one color ; opts_theta@cnLineColors = span_color_rgba ("NCV_jet",11) ; NCV_jet has 256 colors; span it to get 11 colors opts_theta@cnLineThicknessF = 5.0 ; Make lines thicker opts_theta@cnLineLabelsOn = True opts_theta@cnLevelFlag = "LineAndLabel" ; opts_theta@cnMonoLevelFlag = True opts_theta@cnLineLabelInterval = 1 ; opts_theta@cnLineLabelPlacementMode = "constant" ; opts_theta@cnLineLabelPlacementMode = "computed" opts_theta@cnMonoLineColor = False ; opts_theta@cnLineColors = ispan(2,64,1) ; opts_theta@cnLineColors = (/1,2,3,4/) opts_theta@cnLineColors = ispan(10,250,30) opts_theta@cnLineLabelDensityF = 8.0 if (it.eq.69) then system("rm pressure_plane.nc") fdebug = addfile("pressure_plane.nc","c") fdebug->pressure_plane = pressure_plane(:n_YMax/2,:) end if contour_pressure = wrf_contour(a,wks,pressure_plane,opts_theta) delete(opts_theta@ContourParameters ); = (/ .5 /) resvector = True resvector@vpHeightF = 0.35 resvector@vpWidthF = 0.5 resvector@gsnDraw = True ; Forces the plot to be drawn resvector@gsnFrame = True ; Frame advance resvector@vcGlyphStyle = "CurlyVector" resvector@tmXBFormat = "f" resvector@trYMaxF = n_YMax resvector@trXMaxF = n_XMax resvector@trXMinF = n_XMin resvector@vcRefAnnoOn = True resvector@vcRefMagnitudeF = 8. ; define vector ref mag resvector@vcRefLengthF = 0.025 ; define length of vec ref resvector@vcRefAnnoOrthogonalPosF = -1. ; move ref vector resvector@vcRefAnnoParallelPosF = 0.95 ; move ref vector resvector@vcMinDistanceF = 0.025 ; larger means sparser resvector@vcLineArrowHeadMaxSizeF = 0.0075 ; default: 0.05 (LineArrow), 0.012 (CurlyVector) resvector@gsnMaximize = False resvector@vcRefAnnoString2On = True; False resvector@vcRefAnnoString2 = "m s~S~-1" vector = wrf_vector(a,wks,va_plane,wa_plane*100,resvector) ; dummy = gsn_add_polymarker(wks,vector,ij(1),0.08,gsres) ;Contour terrain cross section ; print(opts_ter) ; opts_ter@trXMaxF = X_plane(n_XMax) ; opts_ter@trXMinF = X_plane(n_XMin) ; opts_ter@trYMaxF = 2600; z_values(30)*1000 opts_ter@trXMaxF = n_XMax opts_ter@trXMinF = n_XMin opts_ter@trYMaxF = zcordinate(n_YMax)*1000 opts_ter@vpHeightF = 0.35 opts_ter@vpWidthF = 0.5 opts_ter@xyLineThicknessF = 4.0 opts_ter@xyLineColors = (/"blue","black","blue"/) opts_ter@xyDashPatterns = (/0, 12, 16/) opts_ter@gsnYRefLine = 0.0 ; contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) contour_ter = gsn_csm_y(wks,ter_plane,opts_ter) ; dummy1 = gsn_add_polymarker(wks,contour_ter,X_plane(ij(1)),0.08,gsres) gsres@gsMarkerIndex = 12 dummy1 = gsn_add_polymarker(wks,contour_ter,id_d,58,gsres) print("id_d="+id_d+" position at cross-section") ; MAKE PLOTS pltres = True pltres@NoTitles = True pltres@CommonTitles = True; False pltres@PanelPlot = True pltres@gsnPaperOrientation = "portrait" pltres@FramePlot = False ; plot = wrf_overlays(a,wks,(/contour_theta,vector/),pltres) plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter, contour_pressure/),pltres) Times_char = a->Times(it,:) CST_hr = mod(stringtoint(chartostring(Times_char(11:12)))+8+24, 24) setvalues contour_theta ; "vpHeightF" : 0.35 ; "vpYF" : 0.9-ifile_loop*0.39 ; "vpWidthF" : 0.35 ; "vpXF" : 0.15+ispecies*0.355 "tiMainString" : Times_char(:12)+ " UTC Jun";"Cross section" "tiMainOffsetYF" : -0.01 ; "gsnRightString" : "RString" ; res@TimeLabel ; "gsnPaperOrientation" : "portrait" end setvalues draw(contour_theta) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; delete(ter_plane) delete(wa_plane) delete(tc_plane) delete(theta_plane) delete(pressure_plane) delete(va_plane) delete(WSP_plane); = wrf_user_intrp3d(WSP,z,"v",plane,angle,opts) delete(zcordinate) frame(wks) ; system("convert -density 300 -rotate -90 -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/wrfNARR3dWSM6_CONUS_UCM_YSU.2011080700/"+figurename+".png") ; system("convert -density 300 -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/wrfNARR3dWSM6_CONUS_UCM_YSU.2011080700/"+figurename+".png") print("finish "+figurename+".eps") end do ; END OF TIME LOOP ;end do ; ihour end