; 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 = "x11" ; 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,"rh_tc_diff_vert_CrossSection") ; ----------------------- ; 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 sst.files/wrfout_d03*") numFiles1 = dimsizes(wrffiles1) print("Number of input files: " + numFiles1 ) do i = 0, numFiles1 -1 wrffiles1(i) = wrffiles1(i)+".nc" end do inpFiles1 = addfiles(wrffiles1,"r") ; ------------------------------------------ ; Second set of input files. wrffiles2 = systemfunc("ls no.sst.files/wrfout_d03*") numFiles2 = dimsizes(wrffiles2) do i = 0, numFiles2 -1 wrffiles2(i) = wrffiles2(i)+".nc" end do inpFiles2 = addfiles(wrffiles2,"r") ; ------------------------------------------ do ifile = 0, numFiles1-1 ; LOOP OVER FILES a = inpFiles1[ifile] b = inpFiles2[ifile] times1 = wrf_user_list_times(a) ; get times in the file ntimes1 = dimsizes(times1) ; number of times in the file ; if the file sizes are different, then use following to set plot times ; for your choice ; ; if(ifile.eq.0) then ; times1 = wrf_user_list_times(a) ; get times in the file ; ntimes1 = dimsizes(times1) ; number of times in the file ; istrt = 1 ; the choice of start ; iend = 114 ; the choice of end ; iend = ntimes1 ; printVarSummary(times1) ; end if ; if(ifile.eq.1) then ; times2 = wrf_user_list_times(a) ; get times in the file ; ntimes2 = dimsizes(times2) ; number of times in the file ; istrt = 1 ; iend = 114 ; iend = ntimes2 ; printVarSummary(times2) ; end if ; irpt = 2 irpt = 1 mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file nd = dimsizes(mdims) xlat = wrf_user_getvar(a, "XLAT",0) xlon = wrf_user_getvar(a, "XLONG",0) ter = wrf_user_getvar(a, "HGT",0) ;--------------------------------------------------------------- ; do it = 0,ntimes1-1,2 ; TIME LOOP ; do it = istrt,iend-1,irpt ; TIME LOOP ; do it = 1,5,2 ; TIME LOOP do it = 30,30,2 ; TIME LOOP print("Working on time: " + times1(it) ) res@TimeLabel = times1(it) ; Set Valid time to use on plots ; if(ifile.eq.0) then ; print("Working on time: " + times1(it) ) ; res@TimeLabel = times1(it) ; Set Valid time to use on plots ; end if ; if(ifile.eq.1) then ; print("Working on time: " + times2(it) ) ; res@TimeLabel = times2(it) ; Set Valid time to use on plots ; end if z = wrf_user_getvar(a, "z",it) ; grid point height u1 = wrf_user_getvar(a,"ua",it) ; T in C v1 = wrf_user_getvar(a,"va",it) ; T in C w1 = wrf_user_getvar(a,"wa",it) ; T in C tc1 = wrf_user_getvar(a,"tc",it) ; T in C rh1 = wrf_user_getvar(a,"rh",it) ; relative humidity u2 = wrf_user_getvar(b,"ua",it) ; T in C v2 = wrf_user_getvar(b,"va",it) ; T in C w2 = wrf_user_getvar(b,"wa",it) ; T in C tc2 = wrf_user_getvar(b,"tc",it) ; T in C rh2 = wrf_user_getvar(b,"rh",it) ; relative humidity ; I just want to focus on the first 6 km if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 6. ; We are only interested in the first 6km nz = floattoint(zmax + 1) end if ;--------------------------------------------------------------- ; do ip = 1, 3 ; three cross sections do ip = 1, 1 ; one specific cross section only ; plane = new(2,float) ; plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y) plane = new(4,float) ; Modify start and end x coordinates to match WRF D03 grid size: ; x1 = 0. ; y1 = 129. ; x2 = 327. ; y2 = 269. x1 = 0. y1 = 119. x2 = 327. y2 = 249. plane = (/ x1,y1,x2,y2 /) opts = False angle = ((y2-y1)/(x2-x1))*180. ; print(angle) if(ip .eq. 1) then ; angle = 90. ; plane = (/ 0,129,327,269 /) X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" end if if(ip .eq. 2) then angle = 0. plane = (/ 130,1, 130,162 /) ; start x;y & end x;y point X_plane = wrf_user_intrp2d(xlat,plane,angle,opts) X_desc = "latitude" end if if(ip .eq. 3) then angle = 45. plane = (/ 49,1, 210,162 /) ; start x;y & end x;y point X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" end if u1_plane = wrf_user_intrp3d(u1,z,"v",plane,angle,opts) w1_plane = wrf_user_intrp3d(w1,z,"v",plane,angle,opts) rh1_plane = wrf_user_intrp3d(rh1,z,"v",plane,angle,opts) tc1_plane = wrf_user_intrp3d(tc1,z,"v",plane,angle,opts) u2_plane = wrf_user_intrp3d(u2,z,"v",plane,angle,opts) w2_plane = wrf_user_intrp3d(w2,z,"v",plane,angle,opts) rh2_plane = wrf_user_intrp3d(rh2,z,"v",plane,angle,opts) tc2_plane = wrf_user_intrp3d(tc2,z,"v",plane,angle,opts) u_plane = u1_plane - u2_plane w_plane = w1_plane - w2_plane rh_plane = rh1_plane - rh2_plane tc_plane = tc1_plane - tc2_plane ; Find the index where 6km is - only need to do this once if ( FirstTime ) then zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts) h = ind(zz(:,0) .gt. zmax*1000. ) zmax_pos = h(0) - 1 if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then zspan = h(0) - 1 else zspan = h(0) end if delete(zz) delete(h) FirstTime = False end if ; X-axis lables dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 1) ;--------------------------------------------------------------- ; Options for XY Plots opts_xy = res opts_xy@tiXAxisString = X_desc opts_xy@tiYAxisString = "Height (km)" opts_xy@cnMissingValPerimOn = True ; opts_xy@cnMissingValFillColor = 0 ; opts_xy@cnMissingValFillPattern = 11 opts_xy@tmXTOn = False opts_xy@tmYROn = False opts_xy@tmXBMode = "Explicit" opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels opts_xy@tmXBLabelFontHeightF = 0.015 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; 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 = rh_plane@Orientation opts_xy@NoHeaderFooter = False opts_xy@FieldTitle = " " ; Plotting options for RH opts_rh = opts_xy ; opts_rh@cnLevelSelectionMode = "AutomaticLevels" opts_rh@cnLevelSelectionMode = "ManualLevels" opts_rh@ContourParameters = (/ -10., 10., 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 = "(sst - no sst) RH Difference " ; 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 = opts_xy ; 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 = "(sst - no sst) Temperature Difference " ; if(ifile.eq.0) then ;; opts_tc@FieldTitle = times1(it)+" - (sst - no sst) T(C) & Vert Circ (m/s)" ; opts_tc@FieldTitle = times1(it)+" - (sst - no sst) Temperature" ; else ;; opts_tc@FieldTitle = times2(it)+" - (sst - no sst) T(C) & Vert Circ (m/s)" ; opts_tc@FieldTitle = times2(it)+" - (sst - no sst) Temperature" ; end if ; Get the contour info for the rh and temp contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc) contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),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 = "(sst - no sst) Vertical Circulation Difference" ; 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,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:),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,wks,ter,optsM) plot = wrf_map_overlays(a,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,wks,(/contour_rh,contour_tc/),pltres) ; plot x-section plot = wrf_overlays(a,wks,(/contour_rh,vector/),pltres) ; plot x-section plot = wrf_overlays(a,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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False end do ; END OF TIME LOOP end do ; END OF FILES LOOP end