; 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/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. tkfile = addfile("/global/scratch/umtropea/conus/ctrl/3d/2004/temp/wrf3d_d01_CTRL_TK_20040511.nc","r") zfile = addfile("/global/scratch/umtropea/conus/ctrl/3d/2004/z/wrf3d_d01_CTRL_Z_20040511.nc","r") terfile = addfile("/global/scratch/umtropea/conus/constants/RALconus4km_wrf_constants.nc","r") ufile = addfile("/global/scratch/umtropea/conus/ctrl/3d/2004/wind/wrf3d_d01_CTRL_U_20040511.nc","r") vfile = addfile("/global/scratch/umtropea/conus/ctrl/3d/2004/wind/wrf3d_d01_CTRL_V_20040511.nc","r") wfile = addfile("/global/scratch/umtropea/conus/ctrl/3d/2004/wind/wrf3d_d01_CTRL_W_20040511.nc","r") mergedwindfile = addfile("/global/scratch/umtropea/conus/ctrl/3d/2004/wind/mergedwinds.nc","r") ; Type of plot to generate type = "x11" ; type = "pdf" ; type = "ps" ; type = "ncgm" ; type = "png" wks = gsn_open_wks(type,"vert_cross") ; Set some basic resources res = True ; res@MainTitle = "CTRL May 11, 2004" res@Footer = False pltres = True ter_res = True opts_ter = ter_res opts_ter@gsnYRefLine = 0.0 opts_ter@gsnAboveYRefLineColor = "black" ; opts_ter@gsnAboveYRefLineColor = "darkgreen" opts_ter@gsnDraw = False opts_ter@gsnFrame = False ; gsn_define_colormap (wks,"gui_default") ; gsn_define_colormap (wks,"wgne15") ; gsn_define_colormap (wks,"WhViBlGrYeOrReWh") ; gsn_define_colormap (wks,"ViReOrYeGrAqBl") gsn_define_colormap (wks,"BlAqGrYeOrReVi200") ; gsn_reverse_colormap(wks) ; Reverse the color map. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTime = True FirstTimeMap = True ;;;;;;;;;;;;;;;;;;;;;;;;; Time = wrf_user_getvar(tkfile,"Time",-1) ; get times in the file ntimes = dimsizes(Time) ; number of times in the file mdims = getfilevardimsizes(tkfile,"TK") ; get some dimension sizes for the file nd = dimsizes(mdims) ;Get constants xlat = wrf_user_getvar(terfile,"XLAT",-1) xlon = wrf_user_getvar(terfile,"XLONG",-1) ter = wrf_user_getvar(terfile,"HGT",0) ;--------------------------------------------------------------- ; do it = 0,ntimes-1 ; TIME LOOP do it = 1,1 ; TIME LOOP print("Working on time: " + Time(it) ) res@TimeLabel = Time(it) ; Set Valid time to use on plots ;Get basic variables u = wrf_user_getvar(ufile,"U",it) ; u in m/s v = wrf_user_getvar(vfile,"V",it) ; v in m/s wstag = wrf_user_getvar(wfile,"W",it) ; w in m/s tk = wrf_user_getvar(tkfile,"TK",it) ; T in K zstag = wrf_user_getvar(zfile, "Z",it) ; grid point height ;Convert Kelvin to Celsius and Destagger the Z and W tc = tk - 273.15 z = wrf_user_unstagger(zstag,"Z") w = wrf_user_unstagger(wstag,"Z") ; get height info for labels if ( FirstTime ) then zmin = 0. zmax = 6. ; We are only interested in the first 6km nz = floattoint(zmax + 1) end if opts_ter@trYMaxF = zmax*1000 ;--------------------------------------------------------------- do ip = 1, 3 ;3 plots, center of domain, from A to B plane = new(4,float) plane = (/ 2,2, mdims(nd-1)-2, mdims(nd-2)-2 /) ; start x;y & end x;y point opts = True print(plane) if(ip .eq. 1) then angle = 90. X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" end if if(ip .eq. 2) then angle = 0. X_plane = wrf_user_intrp2d(xlat,plane,angle,opts) X_desc = "latitude" end if if(ip .eq. 3) then angle = 45. X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" end if u_plane = wrf_user_intrp3d(u,z,"v",plane,angle,opts) w_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts) ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) ; print("Max terrain height in plot " + max(ter_plane)) tc_plane2 = tc_plane cross_dims = dimsizes(tc_plane2) rank = dimsizes(cross_dims) ; printVarSummary(tc_plane2) iz_do = 25 do iz = 0,24 iz_do = iz_do-1 do ix = 0,cross_dims(rank-1)-1 if ( ismissing(tc_plane2(iz_do,ix)) ) then tc_plane2(iz_do,ix) = tc_plane2(iz_do+1,ix) end if if ( ismissing(tc_plane2(iz_do,ix)) ) then tc_plane2(iz_do,ix) = tc_plane2(iz_do+1,ix) end if end do end do ; 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) b = ind(zz(:,0) .gt. zmax*1000. ) zmax_pos = b(0) - 1 if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then zspan = b(0) - 1 else zspan = b(0) end if delete(zz) delete(b) FirstTime = False end if ; X-axis labels 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@cnMissingValPerimOn = True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 11 opts_xy@PlotOrientation = tc_plane@Orientation ; X-axis opts_xy@tiXAxisString = X_desc opts_xy@tmXTOn = False opts_xy@tmXBMode = "Explicit" opts_xy@tmXBValues = fspan(0,xspan,nx) opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels opts_xy@tmXBLabelFontHeightF = 0.015 opts_xy@tmXBMajorLengthF = 0.02 opts_xy@tiXAxisFontHeightF = 0.02 ; Y-axis opts_xy@tiYAxisString = "Height (km)" opts_xy@tmYROn = False opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(0,zspan,nz) opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels opts_xy@tmYLLabelFontHeightF = 0.015 opts_xy@tiYAxisFontHeightF = 0.02 opts_xy@tmYLMajorLengthF = 0.02 ; Plotting options for Temperature opts_tc = opts_xy opts_tc@cnInfoLabelZone = 1 opts_tc@cnInfoLabelSide = "Top" opts_tc@cnInfoLabelPerimOn = True opts_tc@cnInfoLabelOrthogonalPosF = -0.00005 opts_tc@ContourParameters = (/ 5. /) opts_tc@cnFillOn = True ; Contour terrain cross section contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) ; Get the contour info for temp contour_tc = wrf_contour(tkfile,wks,tc_plane(0:zmax_pos,:),opts_tc) contour_tc2 = wrf_contour(tkfile,wks,tc_plane2(0:zmax_pos,:),opts_tc) ;------------------------------------------------ ; 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.0075 ; default: 0.05 (LineArrow) vecres@vcLineArrowHeadMaxSizeF = 0.012 ; 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 = "Vert. Circ." ; overwrite field title ; 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(mergedwindfile,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:),vecres) ;--------------------------------------------------------------- ; MAKE PLOTS 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 = True optsM@cnFillOn = True optsM@lbTitleOn = False contour = wrf_contour(terfile,wks,ter,optsM) plot = wrf_map_overlays(terfile,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 ; plot x-section plot = wrf_overlays(tkfile,wks,(/contour_tc,vector/),pltres) ;plot = wrf_overlays(tkfile,wks,(/contour_tc/),pltres) plot = wrf_overlays(tkfile,wks,(/contour_tc,contour_ter/),pltres) ; Delete options and fields, so we don't have carry over delete(opts_xy) delete(opts_tc) delete(u_plane) delete(w_plane) delete(tc_plane) delete(tc_plane2) delete(X_plane) delete(ter_plane) end do ; make next cross section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False end do ; END OF TIME LOOP end