; 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 ; This script adds lon/lat info along X-axis load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin ; ----------------------- ; Type of plot to generate wks_type = "png" ; wks_type = "ncgm" ;wks_type = "ps" wks_type@wkOrientation = "portrait" wks_type@wkPaperWidthF = 8.5 ; in inches wks_type@wkPaperHeightF = 11.0 ; in inches ; wks_type = "pdf" ; wks_type@wkPaperSize = "A4" ; wks_type = "png" ; wks_type@wkWidth = 1000 ; wks_type@wkHeight = 800 ; ----------------------- wks = gsn_open_wks(wks_type,"uv_rh_tc_diff_850hpa") ; ----------------------- ; gsn_define_colormap (wks,"gui_default") ; gsn_define_colormap (wks,"wgne15") ; gsn_define_colormap (wks,"ViReOrYeGrAqBl") ; gsn_define_colormap (wks,"WhViBlGrYeOrReWh") ; gsn_define_colormap (wks,"BlAqGrYeOrReVi200") ; gsn_reverse_colormap(wks) ; Reverse the color map. ; ----------------------- ; Set some basic resources res = True ; res@MainTitle = "REAL-TIME WRF" res@Footer = False ;pltres = True ;pltres@mpGeophysicalLineColor = "Black" ;pltres@mpGeophysicalLineThicknessF = 2 ;pltres@mpNationalLineThicknessF = 2 ; ;pltres@NationalLineColor = "Black" ; pltres@mpOutLineOn = True ;pltres@mpPerimLineColor = "Black" ;pltres@mpGridSpacingF = 45 ;pltres@mpDataBaseVersion = "Ncarg4_1" ;pltres@mpDataSetName = "Earth..4" res@gsnMaximize = True ; maximize plot in frame ; ----------------------- FirstTime = True FirstTimeMap = True ; ------------------------------------------ ; First set of input files. wrffiles1 = systemfunc("ls -1 /home/bchaimite/lustre/WRFChem/WRFV3/run/wrfout_d02*") inpFiles1 = addfiles(wrffiles1,"r") ListSetType(inpFiles1,"cat") ; ------------------------------------------ ; Second set of input files. wrffiles2 = systemfunc("ls -1 /home/bchaimite/lustre/WRFChem/WRFV3/run/WRFCHEM_BBOF/wrfout_d02*") inpFiles2 = addfiles(wrffiles2,"r") ListSetType(inpFiles2,"cat") a = inpFiles1 ;[ifile] b = inpFiles2 ;[ifile] ; Set some Basic Plot options ;res = True ;res@MainTitle = "REAL-TIME WRF" ;res@Footer = False ;pltres = True ;mpres = pltres res@gsnZonalMean = True ; Draw a Zonal mean res@gsnAddCyclic = False ; Data already has cyclic point res@gsnCenterString = "T63" ; Specify resolution pltres = res ; Set plot options res@PanelPlot = True res@mpGeophysicalLineColor = "Black" res@mpNationalLineColor = "Black" res@mpUSStateLineColor = "Black" res@mpGridLineColor = "Black" res@mpLimbLineColor = "Black" res@mpPerimLineColor = "Black" res@mpGeophysicalLineThicknessF = 2.0 res@mpGridLineThicknessF = 2.0 res@mpLimbLineThicknessF = 2.0 res@mpNationalLineThicknessF = 2.0 res@mpUSStateLineThicknessF = 2.0 res@gsnDraw = False ; don't draw yet res@gsnFrame = False res@mpOutlineDrawOrder = "PreDraw" res@tfDoNDCOverlay = "NDCViewport" ; res@tfDoNDCOverlay = True ; old method ;---Set map resources based on projection on WRF output file ; res = wrf_map_resources(a,res) it=-1 z = wrf_user_getvar(a, "z",it) ; grid point height z_avg = dim_avg_n_Wrap(z,0) printVarSummary(z_avg) u1 = wrf_user_getvar(a,"ua",it) ; T in C u1_avg = dim_avg_n_Wrap(u1,0) printVarSummary(u1_avg) v1 = wrf_user_getvar(a,"va",it) ; T in C v1_avg = dim_avg_n_Wrap(v1,0) printVarSummary(v1_avg) w1 = wrf_user_getvar(a,"wa",it) ; T in C w1_avg = dim_avg_n_Wrap(w1,0) printVarSummary(w1_avg) tc1 = wrf_user_getvar(a,"tc",it) ; T in C tc1_avg = dim_avg_n_Wrap(tc1,0) printVarSummary(tc1_avg) rh1 = wrf_user_getvar(a,"rh",it) ; relative humidity rh1_avg = dim_avg_n_Wrap(rh1,0) printVarSummary(rh1_avg) p1 = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate p1_avg = dim_avg_n_Wrap(p1,0) printVarSummary(p1_avg) printMinMax(p1_avg,False) u2 = wrf_user_getvar(b,"ua",it) ; T in C u2_avg = dim_avg_n_Wrap(u2,0) printVarSummary(u2_avg) v2 = wrf_user_getvar(b,"va",it) ; T in C v2_avg = dim_avg_n_Wrap(v2,0) printVarSummary(v2_avg) w2 = wrf_user_getvar(b,"wa",it) ; T in C w2_avg = dim_avg_n_Wrap(w2,0) printVarSummary(w2_avg) tc2 = wrf_user_getvar(b,"tc",it) ; T in C tc2_avg = dim_avg_n_Wrap(tc2,0) printVarSummary(tc2_avg) rh2 = wrf_user_getvar(b,"rh",it) ; relative humidity rh2_avg = dim_avg_n_Wrap(rh2,0) printVarSummary(rh2_avg) plane=850. angle=0. opts=False u1_plane = wrf_user_intrp3d(u1_avg,p1_avg,"h",plane,angle,opts) v1_plane = wrf_user_intrp3d(v1_avg,p1_avg,"h",plane,angle,opts) rh1_plane = wrf_user_intrp3d(rh1_avg,p1_avg,"h",plane,angle,opts) tc1_plane = wrf_user_intrp3d(tc1_avg,p1_avg,"h",plane,angle,opts) u2_plane = wrf_user_intrp3d(u2_avg,p1_avg,"h",plane,angle,opts) v2_plane = wrf_user_intrp3d(v2_avg,p1_avg,"h",plane,angle,opts) rh2_plane = wrf_user_intrp3d(rh2_avg,p1_avg,"h",plane,angle,opts) tc2_plane = wrf_user_intrp3d(tc2_avg,p1_avg,"h",plane,angle,opts) u_plane = u1_plane - u2_plane u_plane@description = "wind u difference" printVarSummary(u_plane) printMinMax(u_plane,False) v_plane = v1_plane - v2_plane v_plane@description = "wind v difference" printVarSummary(v_plane) printMinMax(v_plane,False) rh_plane = rh1_plane - rh2_plane rh_plane@description = "Relative Humidity Difference" printVarSummary(rh_plane) printMinMax(rh_plane,False) tc_plane = tc1_plane - tc2_plane tc_plane@description="Temperature Differences" printVarSummary(tc_plane) printMinMax(tc_plane,False) ; Plotting options for RH res = wrf_map_resources(a[0],res) opts_rh = res ; opts_rh@cnLevelSelectionMode = "AutomaticLevels" opts_rh@cnLevelSelectionMode = "ManualLevels" opts_rh@ContourParameters = (/ -2., 2., .1 /) ; opts_rh@pmLabelBarOrthogonalPosF = -0.1 opts_rh@cnFillOn = True opts_rh@NoHeaderFooter = False opts_rh@Footer = False ; no footer ; opts_rh@units = " " ; turn off units ; opts_rh@gsnLeftString = "LeftString" ; add the gsn titles ; opts_rh@gsnCenterString = "centerstring" ; opts_rh@gsnRightString = "RightString" opts_rh@FieldTitle = "Relative Humidity Difference at 850hPa" ; overwrite field title ; if(ifile.eq.0) then ;; opts_rh@FieldTitle = times1(it)+" - RH(%) & Vert Circ (m/s)" ;overwrite field title ; opts_rh@FieldTitle = times1(it)+" - RH" ;overwrite field title ; else ;; opts_rh@FieldTitle = times2(it)+" - RH(%) & Vert Circ (m/s)";overwrite field title ; opts_rh@FieldTitle = times2(it)+" - RH";overwrite field title ; end if ; opts_rh@cnFillColors = (/"White","White","White", \ ; "White","Chartreuse","Green", \ ; "Green3","Green4", \ ; "ForestGreen","PaleGreen4"/) ; Plotting options for Temperature opts_tc = res ; opts_tc@cnLevelSelectionMode = "AutomaticLevels" opts_tc@cnLevelSelectionMode = "ManualLevels" opts_tc@ContourParameters = (/ -1., 1., .1 /) ; opts_tc@cnInfoLabelZone = 1 ; opts_tc@cnInfoLabelSide = "Top" ; opts_tc@cnInfoLabelPerimOn = True ; opts_tc@cnInfoLabelOrthogonalPosF = -0.00005 ; opts_tc@units = " " ; turn off units ; opts_tc@gsnLeftString = "LeftString" ; add the gsn titles ; opts_tc@gsnCenterString = "centerstring" ; opts_tc@gsnRightString = "RightString" opts_tc@NoHeaderFooter = False opts_tc@Footer = False ; no footer opts_tc@cnFillOn = True ; opts_tc@ContourParameters = (/ 5. /) opts_tc@FieldTitle = "Temperature Difference at 850hPa " ; Get the contour info for the rh and temp contour_tc = wrf_contour(a[0],wks,tc_plane,opts_tc) contour_rh = wrf_contour(a[0],wks,rh_plane,opts_rh) ;------------------------------------------------ ; curly vector plot ;------------------------------------------------ vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcRefMagnitudeF = 1. ; define vector ref mag vecres@vcRefLengthF = 0.045 ; define length of vec ref vecres@vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector vecres@vcRefAnnoArrowLineColor = "red" ; change ref vector color vecres@vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref vecres@vcMinDistanceF = 0.02 ; larger means sparser vecres@vcLineArrowHeadMaxSizeF = 0.01 ; default: 0.05 (LineArrow) ; 0.012 (CurlyVector) vecres@vcGlyphStyle = "CurlyVector" ; turn on curley vectors vecres@vcLineArrowColor = "black" ; change vector color vecres@vcLineArrowThicknessF = 1.2 ; change vector thickness vecres@vcVectorDrawOrder = "PostDraw" ; draw vectors last vecres@FieldTitle = "Wind Circulation Difference at 850hPa" ; vecres@units = " " ; turn off units ; vecres@gsnLeftString = "LeftString" ; add the gsn titles ; vecres@gsnCenterString = "centerstring" ; vecres@gsnRightString = "RightString" vecres@gsnMaximize = True ; maximize plot in frame vecres@NoHeaderFooter = False ; no model info vecres@Footer = False ; no footer vector = wrf_vector(a[0],wks,u_plane,v_plane,vecres) ;--------------------------------------------------------------- ; MAKE PLOTS ;gsn_define_colormap (wks,"BlAqGrYeOrReVi200") ;if (FirstTimeMap) then ; lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts) ;lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts) ;mpres = True ; pltres = True ; pltres@FramePlot = False ; optsM = res ; optsM@NoHeaderFooter = False ; optsM@cnFillOn = True ; optsM@lbTitleOn = False ; contour = wrf_contour(a[0],wks,ter,optsM) ;plot = wrf_map_overlays(a[0],wks,(/contour/),pltres,mpres) ; lnres = True ; lnres@gsLineThicknessF = 3.0 ;lnres@gsLineColor = "Red" ;do ii = 0,dimsX(0)-2 ; gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) ;end do ;frame(wks) ;delete(lon_plane) ;delete(lat_plane) ;pltres@FramePlot = True ;end if ;gsn_define_colormap (wks,"BlAqGrYeOrReVi200") ; plot = wrf_overlays(a[0],wks,(/contour_rh,contour_tc/),pltres) ; plot x-section plot = wrf_overlays(a[0],wks,(/contour_rh,vector/),pltres) ; plot x-section plot = wrf_overlays(a[0],wks,(/contour_tc,vector/),pltres) ; plot x-section ; plot = wrf_overlays(a,wks,(/contour_rh/),pltres) ; plot x-section ; plot = wrf_overlays(a,wks,(/contour_tc/),pltres) ; plot x-section ; Delete options and fields, so we don't have carry over ;delete(opts_xy) ;delete(opts_tc) ;delete(opts_rh) ;delete(tc_plane) ;delete(rh_plane) ;delete(X_plane) ; end do ; make next cross section ; Panel the WRF plots. pnlres = True pnlres@txString = "" pnlres@gsnPanelYWhiteSpacePercent = 5 ; Add white space b/w plots. pnlres@gsnPanelScalePlotIndex = 1 ;SHAPE FILES: moz_shp_name = "/home/bchaimite/lustre/shapfiles/MOZ_adm3.shp" tza_shp_name = "/home/bchaimite/lustre/shapfiles/TZA_adm0.shp" lnres = True lnres@gsLineColor = "gray25" lnres@gsLineThicknessF = 0.5 moz_id = gsn_add_shapefile_polylines(wks,plot,moz_shp_name,lnres) tza_id = gsn_add_shapefile_polylines(wks,plot,tza_shp_name,lnres) ;gsn_panel(wks,(/plot/),(/4,3/),pnlres) draw(plot) frame(wks) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False ; end do ; END OF TIME LOOP ;end do ; END OF FILES LOOP end