; 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" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl" ;external GET_IJ "./get_ij.so" ; external file required to find center location of cross section undef("set_res1_list") function set_res1_list() begin res = True res@cnFillOn = True res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels Levels = (/0.0,0.1,0.5,1.0,1.5,2.0,3.0,4.0,5.0,6.0,7.0,8.0/) res@cnLevels = Levels n_Levels = dimsizes(Levels) ; cmap := read_colormap_file("NCV_rainbow2") cmap := read_colormap_file("GMT_wysiwygcont") n_cmap = dimsizes(cmap) print(n_cmap) ; cmap := cmap(50:, :) n_int = n_cmap(0)/n_Levels cmap := cmap(0::n_int, :) cmap(0,:) = 1. ; Make first color white ; cmap(1,:) = 1. ; Make first color white res@cnFillPalette = cmap ; set color map ; res@pmLabelBarOrthogonalPosF = -0.1 res@lbOrientation = "vertical" ; vertical label bar ; opts_fft@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@cnLinesOn = False ; this gives a warning (once) that it's not a valid setting, but it is. ; I think it happens because I reuse the opts for different variables. ; res@cnLineLabelsOn = False res@lbLabelBarOn = True res@cnInfoLabelOn = True res@lbOrientation = "vertical" ; vertical label bar return(res) end undef("set_res2_list") function set_res2_list() begin res = True res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbOrientation = "vertical" ; vertical label bar ; res@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; res@cnLevels = ispan(-24,24,2) ; set the contour levels res@cnLevels = (/0.0,0.1,0.5,1.0,1.5,2.0,3.0,4.0,5.0,6.0,7.0,8.0/) ; cmap := read_colormap_file("NCV_rainbow2") ; cmap := read_colormap_file("BlueWhiteOrangeRed") cmap := read_colormap_file("GMT_wysiwygcont") cmap := cmap(50:, :) ; cmap(0,:) = 1. ; Make first color white ; cmap(1,:) = 1. ; Make first color white res@cnFillPalette = cmap ; set color map res@lbLabelBarOn = True res@cnInfoLabelOn = True res@lbOrientation = "vertical" ; vertical label bar res@lbTitleString = "[m s~S~-2~N~]" ; res @lbTitleString = "[x 1.E10]" return(res) end undef("set_res3_list") function set_res3_list() begin res = True res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbOrientation = "vertical" ; vertical label bar ; res@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame ; res@cnFillPalette = "MPL_gist_ncar" res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; res@cnLevels = ispan(0,12,1) ; set the contour levels cmap := read_colormap_file("NCV_rainbow2") cmap := cmap(50:, :) ; cmap(0,:) = 1. ; Make first color white ; cmap(1,:) = 1. ; Make first color white res@cnFillPalette = cmap res@cnLevels = (/-2.0,-1.75,-1.5,-1.25,-1.0,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1.0,1.25,1.50,1.75,2.0,\ 2.5,3.0/) ; res@cnFillColors = (/97,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ ; res@cnFillColors = (/0,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ ; 218,225,232/) res@cnLinesOn = False res@lbLabelBarOn = True res@cnInfoLabelOn = True res@lbOrientation = "vertical" ; vertical label bar ;--- These are required for gsn_csm_contour_map plotting ;res = wrf_map_resources(f,res) res =True res@tfDoNDCOverlay = True res@gsnAddCyclic = False res @lbTitleString = "[x 1.E10]" return(res) end undef("set_res4_list") function set_res4_list() begin res = True res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbOrientation = "vertical" ; vertical label bar ; res@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; res@cnLevels = ispan(-12,12,1) ; set the contour levels res@cnLevels = (/-2.0,-1.75,-1.5,-1.25,-1.0,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1.0,1.25,1.50,1.75,2.0,\ 2.5,3.0/) cmap := read_colormap_file("NCV_rainbow2") cmap := cmap(50:, :) ; cmap(0,:) = 1. ; Make first color white ; cmap(1,:) = 1. ; Make first color white res@lbLabelBarOn = True res@cnInfoLabelOn = True res@lbOrientation = "vertical" ; vertical label bar res @lbTitleString = "[x 1.E10]" return(res) end begin plot_rad1 = new ( 6, graphic ) plot_rad2 = new ( 4, graphic ) plot_rad3 = new ( 4, graphic ) pltres = True pltres@NoTitles = True pltres@PanelPlot = True ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. filename_1 = "wrfout_generic_a" ; weno pos filename_2 = "wrfout_generic_b" filename_3 = "wrfout_generic_c" filename_4 = "wrfout_generic_d" ; mono filename_5 = "wrfout_generic_e" filename_6 = "wrfout_generic_f" ; filename_7 = "wrfout_generic_31" ; sla ; filename_8 = "wrfout_generic_32" ; filename_9 = "wrfout_generic_33" filename_7 = "wrfout_d01_0001-01-01_00:00:00_sla_1km_gamma_0.25" filename_8 = "wrfout_d01_0001-01-01_00:00:00_sla_500m_eta_gamma_0.25" filename_9 = "wrfout_d01_0001-01-01_00:00:00_sla_eta_250m_g_0.25" a = addfile(filename_1,"r") b = addfile(filename_2,"r") c = addfile(filename_3,"r") d = addfile(filename_4,"r") e = addfile(filename_5,"r") f = addfile(filename_6,"r") g = addfile(filename_7,"r") h = addfile(filename_8,"r") i = addfile(filename_9,"r") ; my additional code -- set bounds xlat = wrf_user_getvar(a, "XLAT",0) xlon = wrf_user_getvar(a, "XLONG",0) ter = wrf_user_getvar(a, "HGT",0) dims2d=dimsizes(xlat) printVarSummary(xlat) ; lats = (/ 39.00, 42.50 /) ; lons = (/ -75.0, -69.87 /) ; lats = (/ 39.00, 42.50 /) ; lons = (/ -76.0, -70.0/) ; lats = (/ 31.00, 33.00 /) ; lons = (/ 34.00, 35.75/) ; lats = (/ 31.50, 33.50 /) ; lons = (/ 33.5, 35.5/) ; loc = wrf_user_ll_to_ij(a, lons, lats, True) ; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y) ; subtract one since we want to use it as an index in NCL ; We generate plots, but what kind do we prefer? ; type = "x11" ; type = "png" type = "pdf" type@wkWidth = 2500 ; produces a better looking PNG type@wkHeight = 2500 ; type = "pdf" ; type = "ps" ; type = "ncgm" wks = gsn_open_wks(type,"vert_cross_sections_wa_250m") ; wks = gsn_open_wks(type,"vert_cross_sections_qc_qr_1_gamma_0.25") ; wks = gsn_open_wks(type,"vert_cross_sections_qs_qg_1_015_pos_mon_weno") ; wks = gsn_open_wks(type,"vert_cross_sections_qc_qr_1_015_pos_mon_weno") ; wks = gsn_open_wks(type,"vert_cross_sections_qc_qr_1_015_pos_mon_weno") ; wks = gsn_open_wks(type,"plt_w_elec_cross_sections_00_deg") ; Set some basic resources ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 yr = str_get_cols(times, 0, 3) year = stringtoint(yr) mnth = str_get_cols(times, 5, 6) month = stringtoint(mnth) dy = str_get_cols(times, 8, 9) day = stringtoint(dy) hr = str_get_cols(times, 11, 12) hour = stringtoint(hr) mn = str_get_cols(times, 14, 15) minute = stringtoint(mn) print("hour = " + hour) ; not needed with zoom ; mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file ; nd = dimsizes(mdims) opts_rad = set_res1_list() opts_zdr = set_res1_list() opts_crs = set_res3_list() opts_kdp = set_res4_list() ; ;--------------------------------------------------------------- ; do it = 0,ntimes-1,2 ; TIME LOOP do it = 60,60,1 ; do it = 16,22,1 ; do it = 16,16,1 it_val = it new_time_units = "hours since 1800-01-01 00:00" new_time_units = new_time_units new_time_array = cd_inv_calendar( year(it_val), month(it_val), day(it_val), \ hour(it_val), minute(it_val), 0, new_time_units,0) print("new_time_array = " + new_time_array) ; format = "" ;; defaults to "%H%M IST %d %c %Y" ; format = "%H%M UTC %d %c %Y" format = "%H%M" ; format = "%d %c %Y" new_time_array = new_time_array + 0.01 time_var_1=cd_string(new_time_array,format) print("time_var = " + time_var_1) it_val = it + 5 new_time_array = cd_inv_calendar( year(it_val), month(it_val), day(it_val), \ hour(it_val), minute(it_val), 0, new_time_units,0) print("new_time_array = " + new_time_array) ; format = "" ;; defaults to "%H%M IST %d %c %Y" ; format = "%H%M UTC %d %c %Y" format = "%H%M (5 Min Average)" ; format = "%d %c %Y" new_time_array = new_time_array + 0.01 time_var_2=cd_string(new_time_array,format) print("time_var = " + time_var_2) print("Working on time: " + times(it) ) ; res@TimeLabel = times(it) ; Set Valid time to use on plots dx = a->RDX(0) dx = 1./dx print("dx = " + dx) time = it ff1= wrf_user_getvar(a,"wa",time) print("here 1") ; ff1 = dim_avg_n(ff1_o,0) ff2= wrf_user_getvar(b,"wa",time) print("here 2") ; ff2 = dim_avg_n(ff2_o,0) ff3= wrf_user_getvar(c,"wa",time) ; ff3 = dim_avg_n(ff3_o,0) ; ff4= wrf_user_getvar(d,"wa",time) ; ff4 = dim_avg_n(ff4_o,0) print("here 3") ff5= wrf_user_getvar(e,"wa",time) ; ff5 = dim_avg_n(ff5_o,0) ff6= wrf_user_getvar(f,"wa",time) ; ff6 = dim_avg_n(ff6_o,0) print("here 4") ; ff7_o= g->QSNOW(it:it,:,:,:) ff7 = dim_avg_n(ff7_o,0) ff8_o= h->QSNOW(it:it,:,:,:) ff8 = dim_avg_n(ff8_o,0) ff9_o= i->QSNOW(it:it,:,:,:) ff9 = dim_avg_n(ff9_o,0) print("here 5") ; z = wrf_user_getvar(a, "z",it) ; grid point height tk= wrf_user_getvar(a,"tk",it) p_fl= wrf_user_getvar(a,"p",it) p = p_fl psfc = wrf_user_getvar(a,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho = p/(287.*tk) ; tk:= wrf_user_getvar(b,"tk",it) p_fl:= wrf_user_getvar(b,"p",it) p := p_fl psfc := wrf_user_getvar(b,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) ; print("here 6") tk:= wrf_user_getvar(c,"tk",it) p_fl:= wrf_user_getvar(c,"p",it) p := p_fl psfc := wrf_user_getvar(c,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) print("here 7") ; tk:= wrf_user_getvar(d,"tk",it) p_fl:= wrf_user_getvar(d,"p",it) p := p_fl psfc := wrf_user_getvar(d,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) ; tk:= wrf_user_getvar(e,"tk",it) p_fl:= wrf_user_getvar(e,"p",it) p := p_fl psfc := wrf_user_getvar(e,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) print("here 8") ; tk:= wrf_user_getvar(f,"tk",it) p_fl:= wrf_user_getvar(f,"p",it) p := p_fl psfc := wrf_user_getvar(f,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) print("here 9") ; tk := wrf_user_getvar(g,"tk",it) p_fl := wrf_user_getvar(g,"p",it) p := p_fl psfc := wrf_user_getvar(g,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) ; tk:= wrf_user_getvar(h,"tk",it) p_fl:= wrf_user_getvar(h,"p",it) p := p_fl psfc := wrf_user_getvar(h,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) ; print("here 10") tk:= wrf_user_getvar(i,"tk",it) p_fl:= wrf_user_getvar(i,"p",it) p := p_fl psfc := wrf_user_getvar(i,"PSFC",it) dims3d = dimsizes(tk) p(0,:,:) = 0.5*(p_fl(0,:,:)+psfc(:,:)) do k=0,dims3d(0)-2 p(k+1,:,:) = 0.5*(p_fl(k+1,:,:)+p_fl(k,:,:)) end do rho := p/(287.*tk) ; tc_a = wrf_user_getvar(a,"tc",it) ; T in C tc_b = wrf_user_getvar(b,"tc",it) ; T in C tc_c = wrf_user_getvar(c,"tc",it) ; T in C z_a = wrf_user_getvar(a, "z",it) ; grid point height z_b = wrf_user_getvar(b, "z",it) ; grid point height z_c = wrf_user_getvar(c, "z",it) ; grid point height z_d = wrf_user_getvar(d, "z",it) ; grid point height z_e = wrf_user_getvar(e, "z",it) ; grid point height z_f = wrf_user_getvar(f, "z",it) ; grid point height z_g = wrf_user_getvar(g, "z",it) ; grid point height z_h = wrf_user_getvar(h, "z",it) ; grid point height z_i = wrf_user_getvar(i, "z",it) ; grid point height ; print(x) if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 12. ; see z_ave below for setting height level labels nz = floattoint(zmax + 1) end if ;--------------------------------------------------------------- print("here 5 2") opts = True ; setting start and end times a1D = ndtooned(ff1) dsizes = dimsizes(ff1) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_1 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point a1D := ndtooned(ff4) dsizes = dimsizes(ff4) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_4 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point a1D := ndtooned(ff7) dsizes = dimsizes(ff7) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_7 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point a1D := ndtooned(ff2) dsizes = dimsizes(ff2) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_2 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point a1D := ndtooned(ff5) dsizes = dimsizes(ff5) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_5 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point a1D := ndtooned(ff8) dsizes = dimsizes(ff8) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_8 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point printMinMax(tc_a,False) a1D := ndtooned(ff3) dsizes = dimsizes(ff3) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_3 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point a1D := ndtooned(ff6) dsizes = dimsizes(ff6) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_6 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point a1D := ndtooned(ff9) dsizes = dimsizes(ff9) indices = ind_resolve(maxind(a1D),dsizes) i_max = indices(0,2)-1 j_max = indices(0,1)-1 print("i_max, j_max = " +i_max+ " "+j_max) plane_qc_9 = (/ i_max-61,j_max, i_max+66,j_max /) ; start x;y & end x;y point ; rh_plane = wrf_user_intrp3d(rh,z,"v",plane,0.,opts) ; tc_plane_a = wrf_user_intrp3d(tc_a,z_a,"v",plane_qr_a,0.,opts) ; tc_plane_b = wrf_user_intrp3d(tc_b,z_b,"v",plane_qr_b,0.,opts) ; tc_plane_c = wrf_user_intrp3d(tc_c,z_c,"v",plane_qr_c,0.,opts) dim = dimsizes(ff1) zspan = dim(0) angle = 0. zz = wrf_user_intrp3d(z_a,z,"v",plane_qc_1,angle,opts) b_ind = ind(zz(:,0) .gt. zmax*1000. ) zmax_pos = b_ind(0) - 1 if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then zspan = b_ind(0) - 1 else zspan = b_ind(0) end if print(zspan) ; print(x) ; Options for XY Plots res = True opts_xy = res opts_xy@tiYAxisString = "Height (km)" opts_xy@tiXAxisString = "Distance (km)" opts_xy@AspectRatio = 1.00 opts_xy@cnMissingValPerimOn = True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 11 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@tmYLMode = "Explicit" printVarSummary(zz) dimsz = dimsizes(zz) print(dimsz) ; print(x) zz_span = zz(0:zspan,:) zz_span(0,:)=0. opts_qc = True opts_qc@gsnDraw = False opts_qc@gsnDraw = False opts_qc@FieldTitle = "Snow Mass Content" ; opts_qc = opts_xy opts_qc@TimeLabel = times(it) opts_qc@cnFillOn = True opts_qc@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; opts_qc@cnLevels = ispan(1,7,1) Levels = ispan(-50,50,2) opts_qc@cnLevels = Levels opts_qc@cnLinesOn = False n_Levels = dimsizes(Levels) opts_qc@lbLabelBarOn = False opts_qc@cnInfoLabelOn = False opts_qc@lbOrientation = "vertical" ; vertical label bar opts_qc@pmLabelBarSide = "Right" ; vertical label bar opts_qc@lbTitleString = "[g m~S~-3~N~]" opts_qc@lbTitlePosition = "Bottom" ; title position opts_qc@lbTitleFontHeightF= .015 ; make title smaller ; opts_qc@lbTitleString = "" ; opts_qc@UnitLabel = "g cm~S~-3~N~" ; opts_qc@UnitLabel = "" opts_qc@cnInfoLabelString = "" ; cmap := read_colormap_file("GMT_wysiwygcont") cmap := read_colormap_file("NCV_jaisnd") n_cmap = dimsizes(cmap) print(n_cmap) ; cmap := cmap(50:, :) n_int = n_cmap(0)/n_Levels cmap := cmap(0::n_int, :) cmap(0,:) = 1. ; Make first color white ; cmap(1,:) = 1. ; Make first color white opts_qc@cnFillPalette = cmap ; set color map printVarSummary(zz) dimsz = dimsizes(zz) print(dimsz) ; print(x) z_ave = dim_avg_n(zz,1) z_ave(0) = 0. z_ave = z_ave/1000. print(z_ave) z_dim = dimsizes(z_ave) print(z_dim) found_k = False do k = 1,z_dim-1 if (found_k.eq.False) if (z_ave(k).gt.zmax)then k_max = k-1 found_k = True end if end if end do zz_ave = z_ave(0:k_max) kk_max = dimsizes(zz_ave) print("kk_max = " +kk_max) opts_qc@tmYLMode = "Explicit" opts_qc@tmYLValues = fspan(0,kk_max,6) z_int = kk_max/(6-1) print(opts_qc@tmYLValues) print(dimsz) opts_qc@tmYLLabels = sprintf("%4.0f",zz_ave(0::z_int)) opts_qc@tmXBLabelFontHeightF = 0.05 ; opts_qc@tmYBLabelFontHeightF = 0.05 print("opts_qc@tmYLValues = " + opts_qc@tmYLValues) print("opts_qc@tmYLLabels = " + opts_qc@tmYLLabels) ;getvalues pattern_other ; get plot size for use in creating labelbar ; "vpXF" : vpx ; "vpYF" : vpy ; "vpHeightF" : vph ; "vpWidthF" : vpw ; end getvalues ; opts_qc@vpWidthF = 0.1 ; opts_qc@vpHeightF = 0.2 ; opts_qc@lbBoxMajorExtentF = 0.75 ; puts space between color boxes opts_qc@pmLabelBarWidthF = 0.3 opts_qc@pmLabelBarHeightF = 0.4 ; s@lbMonoFillPattern = False ; lbres@lbFillPatterns = (/0,16,1,2/) opts_qc@lbLabelFontHeightF = 0.025 opts_qc@tiYAxisString = "Height (km)" opts_qc@tiXAxisString = "Distance (km)" opts_qc@tiXAxisFontHeightF= .05 opts_qc@tiYAxisFontHeightF= .05 opts_qc@NoHeaderFooter = True ; opts_tc@NoHeaderFooter = True ; contour_tc_a = wrf_contour(a,wks,tc_plane_a,opts_tc) ; contour_tc_b = wrf_contour(a,wks,tc_plane_b,opts_tc) ; contour_tc_c = wrf_contour(a,wks,tc_plane_c,opts_tc) qc_1_o := wrf_user_intrp3d(ff1,z_a,"v",plane_qc_1,0.,opts) qc_1 = qc_1_o(0:zspan,:) qc_4_o := wrf_user_intrp3d(ff4,z_d,"v",plane_qc_4,0.,opts) qc_4 = qc_4_o(0:zspan,:) qc_7_o := wrf_user_intrp3d(ff7,z_g,"v",plane_qc_7,0.,opts) qc_7 = qc_7_o(0:zspan,:) qc_2_o := wrf_user_intrp3d(ff2,z_b,"v",plane_qc_2,0.,opts) qc_2 = qc_2_o(0:zspan,:) qc_5_o := wrf_user_intrp3d(ff5,z_e,"v",plane_qc_5,0.,opts) qc_5 = qc_5_o(0:zspan,:) qc_8_o := wrf_user_intrp3d(ff8,z_h,"v",plane_qc_8,0.,opts) qc_8 = qc_8_o(0:zspan,:) qc_3_o := wrf_user_intrp3d(ff3,z_c,"v",plane_qc_3,0.,opts) qc_3 = qc_3_o(0:zspan,:) qc_6_o := wrf_user_intrp3d(ff6,z_f,"v",plane_qc_6,0.,opts) qc_6 = qc_6_o(0:zspan,:) qc_9_o := wrf_user_intrp3d(ff9,z_i,"v",plane_qc_9,0.,opts) qc_9 = qc_9_o(0:zspan,:) ; dims_plane_1 = dimsizes(qc_1) print(dims_plane_1) grid_x_1 = new(dims_plane_1(1)+8,float) dx = a->RDX(0) dx = 1./dx do j = 0,dims_plane_1(1) grid_x_1(j) = dx*tofloat(j)/1000. end do dims_plane_2 = dimsizes(qc_2) grid_x_2 = new(dims_plane_2(1)+66,float) dx = b->RDX(0) dx = 1./dx do j = 0,dims_plane_2(1) grid_x_2(j) = dx*tofloat(j)/1000. end do dims_plane_3 = dimsizes(qc_3) grid_x_3 = new(dims_plane_3(1)+66,float) dx = c->RDX(0) dx = 1./dx do j = 0,dims_plane_3(1) grid_x_3(j) = dx*tofloat(j)/1000. end do opts_qc@FieldTitle = "Positive" opts_qc@UnitLabel = "250 m" ; opts_qc@UnitLabel = "" opts_qc@vpWidthF = 0.6 opts_qc@vpHeightF = 0.6 opts_qc@tmXBMode = "Explicit" opts_qc@trXMaxF = 128 opts_qc@tmXBValues = fspan(0,dims_plane_1(1),5) ; Create tick marks opts_qc@tiMainFontHeightF = 0.050 print(grid_x_1(0::8)) opts_qc@tmXBLabels = sprintf("%4.1f",grid_x_1(0::32)) print(opts_qc@tmXBLabels) contour_cloud = wrf_contour(a,wks,qc_1,opts_qc) plot_rad1(0) = wrf_overlays(a,wks,(/contour_cloud/),pltres) opts_qc@FieldTitle = "Monotonic" ; opts_qc@UnitLabel = "500 m" opts_qc@tmXBMode = "Explicit" opts_qc@tmXBValues = fspan(0,dims_plane_2(1),5) ; Create tick marks ; opts_qc@trXMaxF = 64 contour_cloud = wrf_contour(a,wks,qc_2,opts_qc) plot_rad1(1) = wrf_overlays(a,wks,(/contour_cloud/),pltres) opts_qc@FieldTitle = "Weno" ; opts_qc@UnitLabel = "250 m" opts_qc@tmXBMode = "Explicit" opts_qc@tmXBValues = fspan(0,dims_plane_3(1),5) ; Create tick marks ; opts_qc@trXMaxF = 128 contour_cloud = wrf_contour(a,wks,qc_3,opts_qc) plot_rad1(2) = wrf_overlays(a,wks,(/contour_cloud/),pltres) opts_qc@FieldTitle = "Weno Positive" ; opts_qc@UnitLabel = "1 km" ; opts_qc@trXMaxF = 32 opts_qc@tmXBMode = "Explicit" opts_qc@tmXBValues = fspan(0,dims_plane_1(1),5) ; Create tick marks contour_cloud = wrf_contour(a,wks,qc_4,opts_qc) plot_rad1(3) = wrf_overlays(a,wks,(/contour_cloud/),pltres) opts_qc@FieldTitle = "SLA (~F33~g~F~ = 0.25)" ; opts_qc@UnitLabel = "500 m" ; opts_qc@trXMaxF = 64 opts_qc@tmXBMode = "Explicit" opts_qc@tmXBValues = fspan(0,dims_plane_2(1),5) ; Create tick marks contour_cloud = wrf_contour(a,wks,qc_5,opts_qc) plot_rad1(4) = wrf_overlays(a,wks,(/contour_cloud/),pltres) opts_qc@FieldTitle = "SLA (~F33~g~F~ = 0.50)" contour_cloud = wrf_contour(a,wks,qc_6,opts_qc) plot_rad1(5) = wrf_overlays(a,wks,(/contour_cloud/),pltres) ;--------------------------------------------------------------- ; MAKE PLOTS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False ; Panel the WRF plots. Type = "QC + QR" pnlres = True ; pnlres@txString = t@description + " (" + t@units + ")" ; pnlres@txString = Type + " at " + radar_deg + " deg elevation (g/cm^3)" ; pnlres@txString = "Vertical Cross Sections" ; pnlres@txString = "CMC " +time_var_1 + " to " + time_var_2 pnlres@txString = "Vertical Velocity" pnlres@gsnPanelYWhiteSpacePercent = 5 ; Add white space b/w plots. pnlres@gsnPanelLabelBar = True ; Turn on common labelbar pnlres@lbLabelAutoStride = True ; Spacing of lbar labels. pnlres@lbBoxMinorExtentF = 0.13 pnlres@lbLabelFontHeightF = .015 pnlres@lbPerimOn = False pnlres@lbLabelFont = "Helvetica" ; label font pnlres@lbTitleOn = True ; turn on title ; pnlres@lbTitleString = "dBZ" ; title string ; pnlres@lbTitlePosition = "Top" ; title position pnlres@lbTitlePosition = "Bottom" ; title position pnlres@lbTitleFontHeightF= .015 ; make title smaller pnlres@lbTitleDirection = "Across" ; title direction pnlres@lbTitleString = "[g m~S~-3~N~]" pnlres@pmLabelBarWidthF = 0.5 pnlres@pmLabelBarHeightF = 0.15 pnlres@lbTopMarginF = 0.001 pnlres@pmLabelBarOrthogonalPosF = 0.02 ; position wrt plot pnlres@lbTitleOffsetF = -0.4 pnlres@gsnPanelYWhiteSpacePercent =8 ; Add white space b/w plots pnlres@gsnPanelFigureStrings = (/"(a)","(b)","(c)","(d)","(e)","(f)","(g)","(h)","(i)"/) pnlres@amJust = "TopLeft" pnlres@gsnPanelFigureStringsFontHeightF = 0.0175 gsn_panel(wks,(/plot_rad1/),(/3,2/),pnlres) ; gsn_panel(wks,(/plot_rad1/),(/2,2/),pnlres) ; delete(opts_xy) ; delete(opts_tc) ; delete(tc_plane) ; delete(fft_plane) ; delete(ff1_plane) ; delete(ff2_plane) ; delete(ff3_plane) ; delete(X_plane) ; Delete options and fields, so we don't have carry over end do ; END OF TIME LOOP end