; Example series of plotting meteograms with WRF ARW model data ; Add some inrfto to the plot ; Add slp to the plot ; Add a time-Z plot above slp and t2 ; Use wrf_user_ll_to_ij to get the point of interest ;*********************************************** 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" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" ;*********************************************** begin a = addfile("./201202/wrfout_d01_2012-02-17_00:00:00","r") wks = gsn_open_wks ("png","./201202/cross-tke-pbl") pblh = wrf_user_getvar(a,"PBLH",-1) pbl = pblh/1000 tc = wrf_user_getvar(a,"tc",-1) u = wrf_user_getvar(a,"ua",-1) v = wrf_user_getvar(a,"va",-1) tke = wrf_user_getvar(a,"TKE_PBL",-1) dust1 = wrf_user_getvar(a,"DUST_1",-1) dust2 = wrf_user_getvar(a,"DUST_2",-1) dust3 = wrf_user_getvar(a,"DUST_3",-1) dust4 = wrf_user_getvar(a,"DUST_4",-1) dust5 = wrf_user_getvar(a,"DUST_5",-1) z = wrf_user_getvar(a, "height",-1) ;-------Find the ij location for the point if interest-------------- res = True plots = new(4,graphic) llres = True llres@ReturnInt = True ; ip_lon = (/42.23, 42.25, 47.78 ,48.21, 48.74/) ip_lat = (/33.23, 32.05, 30.57, 30.37, 31.34/) ip_locs = (/"Baghdad" , "Nukaib" ,"Basrah", "Abadan" , "Ahvaz"/) do ip = 0 , 4 locij = wrf_user_ll_to_ij(a,ip_lon(ip),ip_lat(ip),llres) locij = locij - 1 ; array pointers in NCL space locX=locij(0) locY=locij(1) locs = ip_locs(ip) taus = ispan(1,25,1) ; integer ; ;;;;;;;;;;;;;;;;;;;;;;;;; get time ;;;;;;;;;;;;;;;;;;;;;;;; ; get time information and strip out the day and hour times_in_file = a->Times dims = dimsizes(times_in_file) times = new(dims(0),string) do i=0,dims(0)-1 times(i) = chartostring(times_in_file(i,2:10)) end do times1 = a->Times Time_3 = wrf_times_c( times1, 3 ) ; integers ddhh = str_get_cols(tostring(Time_3),6,9) ; strings ;----------- extract a time series at a point----------------------------------- pbl_point = pbl(:,locY, locX) tc_point = tc(:,:,locY,locX) tke_point = tke(:,:,locY,locX) u_point = u(:,:,locY, locX) v_point = v(:,:,locY, locX) spd10 = (u10_point*u10_point + v10_point*v10_point)^(0.5) pm1_point = dust1(:,:,locY,locX) pm2_point = dust2(:,:,locY,locX) pm3_point = dust3(:,:,locY,locX) pm4_point = dust4(:,:,locY,locX) pm5_point = dust5(:,:,locY,locX) ;------- Swap the dimensions as we want to plot time on the X axis later tt = tc_point(bottom_top|:,Time|:) ugrid = u_point(bottom_top|:,Time|:) vgrid = v_point(bottom_top|:,Time|:) tke1 = tke_point(bottom_top_stag|:,Time|:) pm1 = pm1_point(bottom_top|:,Time|:) pm3 = pm2_point(bottom_top|:,Time|:) pm2 = pm3_point(bottom_top|:,Time|:) pm4 = pm4_point(bottom_top|:,Time|:) pm5 = pm5_point(bottom_top|:,Time|:) pm = pm1+pm2+pm3+pm4+pm5 spd = (ugrid*ugrid + vgrid*vgrid)^(0.5) ;----------------------------- FirstTime=True if ( FirstTime ) then zmin = 0. zmax = max(z)/1000. ; zmax = 4. nz = floattoint(zmax/2 + 1) FirstTime = False end if ;;;;;;;;;;; ; Find the data span - for use in labels;;;;;;;;;;;;;;;;;;;;;; dim=dimsizes(tt) zspan = dim(0) ;;;;;;;; Find the index where 6km is - only need to do this once ; if ( FirstTime ) then ; zz = wrf_user_intrp3d(z,z,"v",plane,0,opts) ; b = ind(z(:,0) .gt. zmax*1000. ) ; zmax_pos = b(0) - 1 ;if ( abs(z(zmax_pos,0)-zmax*1000.) .lt. abs(z(zmax_pos+1,0)-zmax*1000.) ) then ; zspan = b(0) - 1 ; else ; zspan = b(0) ; end if ; delete(z) ; delete(b) ; FirstTime = False ; end if ;----------------------------------------------------------------------- ; options for XY Plots res2D = True res2D@gsnDraw = False ; Don't draw individual plot. res2D@gsnFrame = False ; Don't advance fram res2D@vpXF = 0.1 ; x location res2D@vpYF = 0.90 ; y location res2D@vpWidthF = 0.75 ; width res2D@vpHeightF = 0.45 ; height ; res2D@gsnMaximize = True res2D@cnInfoLabelOn = False ; res2D@cnFillBackgroundColor = -1 ; res2D@cnConstFLabelBackgroundColor = -1 res2D@tmXTOn = False res2D@tmYLMinorOn = True res2D@tiYAxisString = "Height (km)" res2D@tiYAxisFontHeightF = 0.020 res2D@tmYLMode = "Explicit" res2D@tmYLValues = fspan(0,zspan,nz) ; te tick marks res2D@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels res2D@tmYLMajorLengthF = 0.015 res2D@tmYLLabelFontHeightF = 0.02 res2D@tmYLTickSpacingF = 0.2 res2D@tmYLLabelStride = 0.2 res2D@tiXAxisString = "Day_Time" res2D@tiXAxisFontHeightF = 0.015 res2D@tmXBMode = "Explicit" res2D@tmXBValues = taus res2D@tmXBLabels = ddhh res2D@tmXBLabelJust = "CenterCenter" res2D@tmXBLabelDeltaF = 0.05 ;-- move x-axis labels down res2D@tmXBLabelAngleF = 90. ;-- rotate x-axis labels res2D@tmXBMajorLengthF = 0.01 res2D@tmXBLabelFontHeightF = 0.02 res2D@tmXBTickSpacingF = 1.0 res2D@tmXBMinorOn = False ; no lon minor tickmarks ;;;;;;;;;;;;;;;;;; dust ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; resdust = res2D ; Set basic resources resdust@sfXArray = taus resdust@cnLevelSelectionMode = "ManualLevels" resdust@cnMinLevelValF = 200. resdust@cnMaxLevelValF = 6000. resdust@cnLevelSpacingF = 200. resdust@cnLinesOn = False ;-- turn off contour lines resdust@cnLineLabelsOn = False ;-- turn off line labels resdust@cnLineLabelPlacementMode = "constant" resdust@cnConstFLabelBackgroundColor = 1 resdust@cnFillBackgroundColor = 1 ; resdust@cnConstFLabelBackgroundColor = "Background " resdust@cnFillOn = True ; turns on color fill resdust@cnMissingValFillColor = 0 resdust@gsnSpreadColors = True gsn_define_colormap(wks,"WhBlGrYeRe") resdust@pmLabelBarDisplayMode = "Always" ; Add a label bar resdust@pmLabelBarWidthF = 0.10 resdust@pmLabelBarOrthogonalPosF = 0.02 resdust@pmLabelBarParallelPosF = 0.45 resdust@pmLabelBarSide = "Right" resdust@lbOrientation = "vertical" resdust@lbPerimOn = False resdust@lbBoxLinesOn = False resdust@gsnDraw = False; Don't draw individual plot. resdust@gsnFrame = False ; Don't advance fram ;;;;;;;;;;;; wind ;;;;;;;;;;;;;;;;;;;;;;;;;; uv_res = res2D uv_res@sfXArray = taus ; uv_res@vcRefAnnoOn = False ; turns off the ref vector ; uv_res@vcRefLengthF = 0.020 ; set length of ref vector ; uv_res@vcGlyphStyle = "WindBarb" ; turn on wind barbs uv_res@cnLineLabelPlacementMode = "constant" uv_res@cnLevelSelectionMode = "ManualLevels" uv_res@cnLevelSpacingF = 2 uv_res@cnInfoLabelOn = False uv_res@gsnDraw = False ; Don't draw individual plot. uv_res@gsnFrame = False ; Don't advance fram ;;;;;;;;;;;;;; tke ;;;;;;;;;;;;;;;;;;;;;; restke = res2D restke@sfXArray = taus restke@cnMinLevelValF = 0. restke@cnMaxLevelValF = 3. restke@cnLevelSpacingF = 0.5 restke@cnLineLabelPlacementMode = "constant" restke@cnLineLabelFontColor = "Red" restke@cnConstFLabelFontColor = (/2/) restke@cnLineColor = "Red" restke@cnInfoLabelOn = False restke@cnLineThicknessF = 10.0 restke@gsnDraw = False ; Don't draw individual plot. restke@gsnFrame = False ; Don't advance fram ;;;;;;;;;;;;;;;;;;;;;; plot one dimenstion;;;;;;;;;;;;;;;;; res1D = True res1D@gsnDraw = False ; Don't draw individual plot. res1D@gsnFrame = False ; Don't advance fram res1D@vpXF = 0.1 ; The left side of the box location res1D@vpWidthF = 0.75 ; The Width of the plot box res1D@xyLineColor = 1 ; set line color res1D@xyLineThicknesses = 10 ; increase line thickness ;;;;;;;;;;;;;;;;;;;;;;;; pbl ;;;;;;;;;;;;;;;;;;;;;;;; pbl_res = res1D pbl_res@gsnDraw = False ; Don't draw individual plot. pbl_res@gsnFrame = False ; Don't advance frame pbl_res@sfYArray = fspan(0,zspan,nz) pbl_res@amDataYF = fspan(0,zspan,nz) pbl_res@sfXArray = taus pbl_res@vpYF = 0.3 ; The top side of the plot box pbl_res@vpHeightF = 0.2 pbl_res@tmXBMode = "Explicit" pbl_res@tmXBLabelsOn = True pbl_res@tmXTOn = False pbl_res@tmXBMinorOn = True pbl_res@tmXBMajorLengthF = 0.01 pbl_res@tmXBTickSpacingF = 1.0 pbl_res@tmXBLabelStride = 1 pbl_res@tmXBLabelJust = "CenterCenter" pbl_res@tmXBLabelFontHeightF = 0.02 pbl_res@tiXAxisFontHeightF = 0.02 pbl_res@tiXAxisString = "Day_Time" pbl_res@tiXAxisFontHeightF = 0.015 pbl_res@tmXBValues = taus pbl_res@tmXBLabels = ddhh pbl_res@tmXBLabelDeltaF = 0.11 ;-- move x-axis labels down pbl_res@tmXBLabelAngleF = 90. ;-- rotate x-axis lab pbl_res@tiYAxisString = "Height (km)" pbl_res@tiYAxisFontHeightF = 0.020 pbl_res@tmYLMode = "Explicit" pbl_res@tmYLValues = fspan(0,zspan,nz) ; te tick marks pbl_res@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels pbl_res@tmYLMajorLengthF = 0.015 pbl_res@tmYLLabelFontHeightF = 0.02 pbl_res@tmYLTickSpacingF = 0.2 pbl_res@tmYLLabelStride = 0.2 pbl_res@tmYLMode = "manual" pbl_res@tmYLMajorLengthF = 0.02 pbl_res@tmYLLabelFontHeightF = 0.015 pbl_res@tmYLTickSpacingF = 0.2 pbl_res@tiYAxisString = " PBL Height(km)" pbl_res@tiYAxisFontHeightF = 0.020 pbl_res@tmYLLabelStride = 0.2 ;;;;;;;;;;;;;;;;;;;;;;;;;;; station data ;;;;;;;;;;;;;;;;;; ndcres = True ndcres@txFontHeightF = 0.017 ndcres@txJust = "BottomLeft" gsn_text_ndc(wks, (ip_locs(ip)), 0.28, 0.9, ndcres) ; gsn_text_ndc(wks, ("lon="+(ip_lon(ip))), 0.28, 0.9, ndcres) ; gsn_text_ndc(wks, ("lat="+(ip_lat(ip))), 0.39, 0.9, ndcres) ;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ pmfill = gsn_contour(wks,pm,resdust) windfill = gsn_contour(wks,spd,uv_res) tkefill = gsn_contour(wks,tke1,restke) pbl_plot = gsn_csm_xy(wks,taus,pbl_plot,pbl_res) overlay(pmfill,windfill) overlay(pmfill,tkefill) overlay(pmfill,pbl_plot) draw(pmfill) frame(wks) ; now frame the plot ;----------------------------------------------------------------------- end do end