; 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" begin Reanalysis = "NARR" configuration = "2d_UCM_Snudge_WangRec53wave" simname = Reanalysis+configuration simname_short = simname simname_output = Reanalysis+configuration+".2002_2015Mean" gsres = True gsres@gsMarkerIndex = 19 ; circle at first gsres@gsMarkerThicknessF = 1 gsres@gsMarkerSizeF = 0.009105 gsres@gsMarkerColor = "black" ; gsres@gsMarkerColor = "red" Meteo_Lat = 32.7833333 ; Dallas Meteo_Lon = -96.8 idomain = 2 ;1 ; files =systemfunc("ls wrfout_d0"+idomain+"_??_ncea_14yrAugmean_wrfNARR2d_UCM_Snudge_WangRec53wave.nc") ;do ifiles = 7,7 files = systemfunc("ls wrfout_d0"+idomain+"_21_ncea_14yrAugmean_wrfNARR2d_UCM_Snudge_WangRec53wave.nc") do ifiles = 0 , 0; dimsizes(files)-1 a = addfile(files(ifiles)+".nc","r") opt = True ij = wrf_user_ll_to_ij(a,Meteo_Lon, Meteo_Lat,opt) ; print(ij) ; We generate plots, but what kind do we prefer? type = "eps" type@wkOrientation = "portrait" ; type@wkWidth = 810 ; type@wkHeight = 810 figurename = "wrfout_d0"+idomain+"_W_crossSN90_2ndMax_smooth_partial_8";+ifiles wks = gsn_open_wks(type,figurename) gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map ; Set some basic resources res = True ; res@MainTitle = "REAL-TIME WRF" res@Footer = False res@gsnMaximize = False ter_res = True opts_ter = ter_res opts_ter@gsnYRefLine = 0.0 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 = 0,ntimes-1,2 ; TIME LOOP 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) ISLTYP = a->ISLTYP(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) ; height = wrf_user_getvar(a, "height",it) ; 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 wa = (/wa*100./) ; convert m/s to cm/s QRAIN = a->QRAIN(0,:,:,:) QRAIN = (/QRAIN*1e6/) ; convert to g/kg QRAIN@units = "mg kg~S~-1~N~";(/temp*1000/) printVarSummary(wa) printVarSummary(tc) nsmooth = 10 do ismooth = 0, nsmooth wa = smth9(wa, 0.5, 0.25, False) end do ; do ismooth = 0, 0; nsmooth QRAIN = smth9(QRAIN, 0.5, 0.25, False) end do ; ;--------------------------------------------------------------- ; 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 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 printVarSummary(height) print(height(:10,ij(1),ij(0))) print(ter(ij(1),ij(0))) print(xlon(ij(1),ij(0))) print(xlat(ij(1),ij(0))) 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) ua_plane = wrf_user_intrp3d(ua,z,"v",plane,angle,opts) wa_plane = wrf_user_intrp3d(wa,z,"v",plane,angle,opts) ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) QRAIN_plane = wrf_user_intrp3d(QRAIN,z,"v",plane,angle,opts) ; ISLTYP_plane = wrf_user_intrp2d(ISLTYP*1.0,plane,angle,opts) lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts) lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts) printVarSummary(wa_plane) printVarSummary(lon_plane) distance = (lon_plane( :)-Meteo_Lon)^2+(lat_plane(:)-Meteo_Lat)^2 id_d = minind(distance) print("site location = "+id_d) zmin = 0. zmax = max(z)/1000. dim = dimsizes(WSP_plane) ; Find the data span - for use in labels zspan = dim(0) print("zmax="+zmax+" zspan="+zspan) zcordinate = fspan(zmin,zmax,zspan) ; printVarSummary(zcordinate) dim = dimsizes(theta_plane) ; Find the data span - for use in labels zspan = dim(0) res@gsnPaperOrientation = "portrait" ; Options for XY Plots opts_theta = res opts_theta@tiYAxisString = "Height (km)" opts_theta@cnMissingValPerimOn = False ; True opts_theta@cnMissingValFillColor = 0 opts_theta@cnMissingValFillPattern = 11 opts_theta@pmLabelBarOrthogonalPosF = .10 ; move whole thing down opts_theta@tmYLMode = "Explicit" opts_theta@tmYLValues = fspan(0,zspan,20) ; Create tick marks opts_theta@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,20)) ; Create labels opts_theta@tmXBMode = "Explicit" opts_theta@tmXBValues = fspan(0,dim(1)-1,43) ; Create tick marks opts_theta@tmXBLabels = sprintf("%.1f",fspan(lon_plane(0),lon_plane(dim(1)-1),43)) ; Create labels opts_theta@tiYAxisFontHeightF = 0.02820 ; not effective opts_theta@tmYLLabelFontHeightF = 0.0282 opts_theta@tiXAxisFontHeightF = 0.002091020 opts_theta@tmXBLabelFontHeightF = 0.000005109102 opts_theta@lbLabelFontHeightF = 0.0102 opts_theta@lbTitleFontHeightF= .01028 ; make title smaller opts_theta@tmXBMajorLengthF = 0.02 opts_theta@tmYLMajorLengthF = 0.02 opts_theta@PlotOrientation = tc_plane@Orientation n_XMax = 440; id_d+100 n_XMin = 310; id_d-115 opts_theta@trXMaxF = n_XMax opts_theta@trXMinF = n_XMin do ispecies = 0 ,1 ; opts_theta@gsnOrientation = "portrait" opts_theta@gsnDraw = True ; Forces the plot to be drawn opts_theta@gsnFrame = True ; Frame advance opts_theta@tiMainString = "" ;chartostring(a->Times) opts_theta@vpHeightF = 0.35 opts_theta@vpWidthF = 0.5 ; opts_theta@gsnRightString = chartostring(a->Times) ; Plotting options for RH opts_theta@cnFillMode = "RasterFill" ; opts_theta@ContourParameters = (/ -3., 3., .5 /) opts_theta@tmXBFormat = "f" opts_theta@pmLabelBarOrthogonalPosF = -1.0995097 opts_theta@pmLabelBarParallelPosF = 1.069877 opts_theta@cnFillOn = True opts_theta@gsnSpreadColors = True ; use full range of colormap opts_theta@lbLeftMarginF = 0.25 opts_theta@lbRightMarginF = 0.4 opts_theta@lbOrientation = "vertical" opts_theta@lbTitlePosition = "Top" ; title position opts_theta@lbTitleDirection = "Across" ; title direction opts_theta@lbTopMarginF = 0.42 opts_theta@lbBottomMarginF = 0.42 opts_theta@lbTitleOffsetF = 0.012 opts_theta@gsnMaximize = False if(ispecies.eq.0) then opts_theta@ContourParameters = (/ -2., 2., .25 /) n_YMax = 200 opts_theta@trYMaxF = n_YMax opts_theta@lbTitleOn = True opts_theta@lbTitleString = "w~C~cm s~S~-1~N~" opts_theta@lbLabelStride = 4 contour_theta = wrf_contour(a,wks,wa_plane,opts_theta) else opts_theta@ContourParameters = (/ 2., 12.,1. /) n_YMax = 350 opts_theta@trYMaxF = n_YMax opts_theta@lbTitleOn = True opts_theta@lbTitleString = "QRAIN~C~mg kg~S~-1~N~" opts_theta@lbLabelStride = 2 contour_theta = wrf_contour(a,wks,QRAIN_plane,opts_theta) end if 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@trXMinF = 65. resvector@trYMaxF = n_YMax resvector@trXMaxF = n_XMax resvector@trXMinF = n_XMin resvector@vcRefAnnoOn = False ; True resvector@vcRefMagnitudeF = 5. ; 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 vector = wrf_vector(a,wks,ua_plane,wa_plane,resvector) ;Contour terrain cross section opts_ter@trXMaxF = lon_plane(n_XMax) opts_ter@trXMinF = lon_plane(n_XMin) opts_ter@trYMaxF = zcordinate(n_YMax)*1000 opts_ter@vpHeightF = 0.35 opts_ter@vpWidthF = 0.5 opts_ter@xyLineThicknessF = 12.0 ISLTYP_line = ISLTYP(ij(1),:) ; ISLTYP_line = where(ISLTYP_line.ne.1.and.ISLTYP_line.ne.3.and.ISLTYP_line.ne.9.and.ISLTYP_line.ne.12, ISLTYP_line@_FillValue, ISLTYP_line) ; ISLTYP_line(330:345) = ISLTYP_line@_FillValue ISLTYP_line = where(ISLTYP_line.ne.1.and.ISLTYP_line.ne.3.and.ISLTYP_line.ne.9.and.ISLTYP_line.ne.12, 3, ISLTYP_line) ISLTYP_line(330:345) = 3 ISLTYP_line(380:390) = 12 ISLTYP_line(360:365) = 9 ISLTYP_line(371:375) = 3 ISLTYP_line(395:405) = 3 ; contour_ter = gsn_csm_xy(wks,lon_plane,ter_plane,opts_ter) ter_plane_simple = a->HGT(0,ij(1),:) ter_plane_simple12 = where(ISLTYP_line.eq.12.or.ISLTYP_line.eq.9, ter_plane_simple, ter_plane_simple@_FillValue) ; ter_plane_simple9 = where(ISLTYP_line.eq.9, ter_plane_simple, ter_plane_simple@_FillValue) ter_plane_simple3 = where(ISLTYP_line.eq.3.or.ISLTYP_line.eq.1, ter_plane_simple, ter_plane_simple@_FillValue) ; ter_plane_simple1 = where(ISLTYP_line.eq.1, ter_plane_simple, ter_plane_simple@_FillValue) opts_ter@gsnAboveYRefLineColor = "brown";"purple" contour_ter = gsn_csm_xy(wks,xlon(ij(1),:),ter_plane_simple12,opts_ter) ; opts_ter@gsnAboveYRefLineColor = "orange" ; contour_ter9 = gsn_csm_xy(wks,xlon(ij(1),:),ter_plane_simple9,opts_ter) opts_ter@gsnAboveYRefLineColor = "yellow" ; opts_ter@gsnAboveYRefLineColor = "green" contour_ter3 = gsn_csm_xy(wks,xlon(ij(1),:),ter_plane_simple3,opts_ter) ; contour_ter1 = gsn_csm_xy(wks,xlon(ij(1),:),ter_plane_simple1,opts_ter) dummy1 = gsn_add_polymarker(wks,contour_ter,lon_plane(id_d),0.1028,gsres) tres = True tres@txFontColor = "black" tres@txFontHeightF = 0.015005 tres@txBackgroundFillColor = "White" label_xhu = (/"a","b","c","d"/) dummy2 = gsn_add_text(wks,vector,label_xhu(ispecies),313,10,tres) mstrings = (/"u","z","y","q","x","r","Z"/) fontnums = (/34, 35, 35, 35,19, 35,37/) xoffset = 0.0 yoffset = 0.0 ratio = 1.0 sizes = 1.5;(/2.0,2.0,0.75,1.0,2.5,1.5/) angle_marker = 0.0 new_indexes = NhlNewMarker(wks, mstrings, fontnums, xoffset, yoffset, ratio, sizes, angle_marker) opt_soil = opts_ter ; delete(opt_soil@gsnAboveYRefLineColor); = "transparent" opt_soil@gsnAboveYRefLineColor = "transparent" opt_soil@xyLineColors = (/"blue","black","blue"/) opt_soil@xyMarkLineMode = "Markers" ; Markers *and* lines ; opt_soil@xyMarkers = (/19,11,16/) ;19, filled square opt_soil@xyMarkers = (/15,11,16/) ; 3 different markers ; opt_soil@xyMarkerColors := (/"blue","red","green"/) ; 3 different colors 6.3.0 opt_soil@xyMarkerColors = (/"blue","red","green"/) ; 3 different colors 6.3.0 opt_soil@xyLineThicknesses = (/ 3.0, 2.0/) ; make second line thicker opt_soil@tmYRLabelsOn = True opt_soil@tmYROn = True opt_soil@xyDashPatterns = (/0, 16, 12/) opt_soil@trYMaxF = 20. opt_soil@trYMinF = 0. contour_soil = gsn_csm_xy(wks,xlon(ij(1),:),ISLTYP_line,opt_soil) ; MAKE PLOTS pltres = True pltres@NoTitles = True pltres@CommonTitles = True; False pltres@PanelPlot = True pltres@gsnPaperOrientation = "portrait" pltres@gsnMaximize = False pltres@FramePlot = False ; plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter,contour_soil/),pltres) ; plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter,contour_ter9,contour_ter3,contour_ter1,contour_soil/),pltres) ; plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter,contour_ter3,contour_soil/),pltres) plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter,contour_ter3/),pltres) setvalues contour_theta "vpHeightF" : 0.34 "vpWidthF" : 0.65 ; "vpYF" : 0.93-ispecies*0.34084 "vpYF" : 0.9-ispecies*0.3484 ; "vpXF" : 0.15+ispecies*0.355 "tiMainOffsetYF" : -0.01 ; "gsnRightString" : "RString" ; res@TimeLabel "gsnPaperOrientation" : "portrait" "tmXBLabelFontHeightF" : 0.01482 "tiYAxisFontHeightF" : 0.014820 ; effective "tmYLLabelFontHeightF" : 0.01482 "tiXAxisFontHeightF" : 0.01482 "tiYAxisOffsetXF" : 0.0102 "tmYLLabelDeltaF" : -0.68 end setvalues if(ispecies.eq.0) then setvalues contour_theta "tmXBLabelsOn" : False "tiMainString" : "" ; res@TimeLabel + " UTC";"Cross section" ; "lbTitleOn" : True ; "lbTitleString" : "w~C~cm s~S~-1~N~" end setvalues else setvalues contour_theta ; "lbLabelBarOn" : True ; "tmXBLabelsOn" : True "tiMainString" : "" ; res@TimeLabel + " UTC";"Cross section" "tiXAxisString" : "Longitude" end setvalues end if draw(contour_theta) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; END OF TIME LOOP end do ; ispecies delete(ter_plane) delete(wa_plane) delete(tc_plane) delete(theta_plane) delete(va_plane) delete(ua_plane) delete(WSP_plane); delete(zcordinate) frame(wks) ; system("rm "+figurename+".eps") ; system("mv "+figurename+".png /nsftor/xhu/public_html/WRF-UCM/FocusON_OKC/WRFV3.7.1/with_rr_fixed/"+"wrf"+simname_output+"/"+figurename+".png") ; system("mv "+figurename+".eps /nsftor/xhu/public_html/WRF-UCM/FocusON_OKC/WRFV3.7.1/with_rr_fixed/"+"wrf"+simname_output+"/"+figurename+".eps") end do ; ifiles end