; 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 "/usr/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "/usr/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "/usr/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. dir = "~/Desktop/" fili = "wrfout_d04_2015-09-11_2240.nc" filin= dir+fili a = addfile(filin,"r") ; We generate plots, but what kind do we prefer? ; type = "x11" ; type = "pdf" ; type = "ps" ; type = "ncgm" type = "png" ; Set some basic resources res = True res@MainTitle = "REAL-TIME WRF" ;res@Footer = False res@gsnMaximize = True pltres = True ter_res = True opts_ter = ter_res opts_ter@gsnYRefLine = 0. opts_ter@gsnAboveYRefLineColor = "black" opts_ter@gsnDraw = False opts_ter@gsnFrame = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTime = True FirstTimeMap = True times = wrf_user_getvar(a,"times",-1) ; get times in the file ntimes = dimsizes(times) ; number of times in the file 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,ntimes-1,1 ; TIME LOOP it=0 print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots type@wkWidth = 2500 type@wkHeight = 2500 wks = gsn_open_wks(type,"nice"+times(it)) ;pv = wrf_user_getvar(a,"pvo",it) th = wrf_user_getvar(a,"th",it) ; potential temp in K eth = wrf_user_getvar(a,"eth",it) ; equivalent potential temp (K) rh = wrf_user_getvar(a,"rh",it) ; relative humidity z = wrf_user_getvar(a, "z",it) ; grid point height uvm = wrf_user_getvar(a,"ua",it) dbz = wrf_user_getvar(a,"dbz",it) p = wrf_user_getvar(a,"pressure",it) w = wrf_user_getvar(a,"wa", it) ; 3D W at mass points aka upward motion w = w*100. w@units = "cm/s" max_w = max(w) ;----------------------------------------------------------------- if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 10. ; We are only interested in the first 6km nz = floattoint(zmax + 1) end if ;opts_ter@trYMaxF = z_values(dimzz-2)*1000 opts_ter@trYMaxF = zmax*1000 ;--------------------------------------------------------------- minlat = 6.502222 maxlat = 6.433889 minlon = -75.826667 maxlon = -75.650278 opt = True locA = wrf_user_ll_to_ij(a,minlon,minlat,opt) locB = wrf_user_ll_to_ij(a,maxlon,maxlat,opt) lonA = locA(0) - 1 latA = locA(1) - 1 lonB = locB(0) - 1 latB = locB(1) - 1 plane = new(4,float) angle = 0. plane = (/ lonA,latA , lonB,latB /) ;--- start x;y & end x;y point opts = opt X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "Longitude" ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) print("Max terrain height in plot " + max(ter_plane)) opts@extrapolate = True rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(th,z,"v",plane,angle,opts) et_plane = wrf_user_intrp3d(eth,z,"v",plane,angle,opts) ;pv_plane = wrf_user_intrp3d(pv,z,"v",plane,angle,opts) uv_plane = wrf_user_intrp3d(uvm,z,"v",plane,angle,opts) wa_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts) dbzplane = wrf_user_intrp3d(dbz,z,"v",plane,angle,opts) ;ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) ;print("Max terrain height in plot " + max(ter_plane)) rh_plane2 = rh_plane tc_plane2 = tc_plane ;pv_plane2 = pv_plane et_plane2 = et_plane uv_plane2 = uv_plane wa_plane2 = wa_plane dbzplane2 = dbzplane ;epv_plane2= epv_plane cross_dims = dimsizes(rh_plane2) rank = dimsizes(cross_dims) ;printVarSummary(rh_plane2) iz_do = 25 do iz = 0,24 iz_do = iz_do-1 do ix = 0,cross_dims(rank-1)-1 if ( ismissing(rh_plane2(iz_do,ix)) ) then rh_plane2(iz_do,ix) = rh_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 ;if ( ismissing(pv_plane2(iz_do,ix)) ) then ; pv_plane2(iz_do,ix) = pv_plane2(iz_do+1,ix) ;end if if ( ismissing(et_plane2(iz_do,ix)) ) then et_plane2(iz_do,ix) = et_plane2(iz_do+1,ix) end if if ( ismissing(uv_plane2(iz_do,ix)) ) then uv_plane2(iz_do,ix) = uv_plane2(iz_do+1,ix) end if if ( ismissing(wa_plane2(iz_do,ix)) ) then wa_plane2(iz_do,ix) = wa_plane2(iz_do+1,ix) end if if ( ismissing(dbzplane2(iz_do,ix)) ) then dbzplane2(iz_do,ix) = dbzplane2(iz_do+1,ix) end if ;if ( ismissing(epv_plane2(iz_do,ix)) ) then ; epv_plane2(iz_do,ix) = epv_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 lables dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 print(xspan) nx=4 ;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 = tc_plane@Orientation ; Plotting options for RH opts_rh = opts_xy opts_rh@ContourParameters = (/ 70., 98., 2. /) opts_rh@pmLabelBarOrthogonalPosF = -0.1 opts_rh@cnFillOn = True opts_rh@cnFillPalette = "WhiteGreen" ;opts_rh@cnFillColors = (/"White","White","White","White","White","White","White", \ ; "White","White","White","White","White","White","Chartreuse","Green", \ ; "Green3","Green4", \ ; "ForestGreen","PaleGreen4"/) ; 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 = (/ 2.5 /) ;opts_pv = opts_xy ;opts_pv@ContourParameters = (/-10., 10., 1.0/) ;opts_pv@cnFillOn = True ;opts_pv@FieldTitle = "Equivalent Potential Vorticity" opts_uv = opts_xy opts_uv@ContourParameters = (/ -45., 45., 2.5 /) opts_uv@cnFillOn = True opts_uv@cnFillPalette = "NCV_jaisnd" opts_wp = opts_xy opts_wp@ContourParameters = (/-10.,10.,1./) cmap = read_colormap_file("NCV_blue_red") cmap = cmap(::-1,:) ; reverse the color map opts_wp@cnFillOn = True opts_wp@cnFillPalette = cmap opts_r = opts_xy opts_r@cnFillPalette = "WhiteBlueGreenYellowRed" opts_r@ContourParameters = (/ 0., 76., 2. /) opts_r@cnFillOn = True opts_r@gsnSpreadColors = True ;Contour terrain cross section contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) ; Get the contour info for the rh and temp contour_et = wrf_contour(a,wks,et_plane(0:zmax_pos,:),opts_tc) 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) ;contour_pv = wrf_contour(a,wks,pv_plane(0:zmax_pos,:),opts_pv) contour_uv = wrf_contour(a,wks,uv_plane(0:zmax_pos,:),opts_uv) contour_wa = wrf_contour(a,wks,wa_plane(0:zmax_pos,:),opts_wp) contourdbz = wrf_contour(a,wks,dbzplane(0:zmax_pos,:),opts_r) ;contourepv = wrf_contour(a,wks,epv_plane(0:zmax_pos,:),opts_pv) contour_et2 = wrf_contour(a,wks,et_plane2(0:zmax_pos,:),opts_tc) contour_tc2 = wrf_contour(a,wks,tc_plane2(0:zmax_pos,:),opts_tc) contour_rh2 = wrf_contour(a,wks,rh_plane2(0:zmax_pos,:),opts_rh) ;contour_pv2 = wrf_contour(a,wks,pv_plane2(0:zmax_pos,:),opts_pv) contour_uv2 = wrf_contour(a,wks,uv_plane2(0:zmax_pos,:),opts_uv) contour_wa2 = wrf_contour(a,wks,wa_plane2(0:zmax_pos,:),opts_wp) contourdbz2 = wrf_contour(a,wks,dbzplane2(0:zmax_pos,:),opts_r) ;contourepv2 = wrf_contour(a,wks,epv_plane2(0:zmax_pos,:),opts_pv) ;--------------------------------------------------------------- ; 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@cnFillPalette = "OceanLakeLandSnow" optsM@lbTitleOn = False mpres@mpNationalLineColor = "black" mpres@mpUSStateLineColor = "black" mpres@mpUSStateLineThicknessF = 3.0 mpres@mpNationalLineThicknessF = 3.0 mpres@mpGeophysicalLineThicknessF = 3.0 contour = wrf_contour(a,wks,ter,optsM) plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres) lnres = True lnres@gsLineThicknessF = 5.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 = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres) ; plot x-section plot = wrf_overlays(a,wks,(/contour_rh,contour_ter/),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(tc_plane2) delete(rh_plane2) delete(X_plane) delete(ter_plane) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False ;end do ; END OF TIME LOOP