; 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 CALC_LPI "./calclpi3d.so" 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 ; res@cnLevels = ispan(tointeger(1E0),tointeger(1E2),tointeger(2E1)) res@cnLevels = ispan(10,610,60) res@cnInfoLabelOn = False 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 ; 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 res@lbTitleOn = True ; turn on title res@lbTitleString = "J/kg" 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(tointeger(1E0),tointeger(1E2),tointeger(2E1)) ; cmap := read_colormap_file("NCV_rainbow2") cmap := read_colormap_file("BlueWhiteOrangeRed") 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 ; res1@lbTitleString = "[~S~o~N~C]" res @lbTitleString = "[x 1.E-6]" 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@tfDoNDCOverlay = True res@gsnAddCyclic = False res @lbTitleString = "[x 1.E-6]" 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 ) ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. filename_1 = "wrfout_generic_ccn" filename_2 = "wrfout_generic_dust" filename_3 = "wrfout_generic_dust_urb" a = addfile(filename_1,"r") b = addfile(filename_2,"r") c = addfile(filename_3,"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.75, 33.25 /) ; lons = (/ 34.50, 35.25/) lats = (/ 31.75, 33.25 /) lons = (/ 34.50, 35.25/) 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 x_start = loc(1,0) - 1 x_end = loc(1,1) - 1 y_start = loc(0,0) - 1 y_end = loc(0,1) - 1 ; note: x is y here print("x_start = " + x_start) print("x_end = " + x_end) print("y_start = " + y_start) print("y_end = " + y_end) ; latitude longitude of Tel-Aviv xlat_zoom = xlat(x_start:x_end,y_start:y_end) xlon_zoom = xlon(x_start:x_end,y_start:y_end) ter_zoom = ter(x_start:x_end,y_start:y_end) dims_zoom = dimsizes(xlat_zoom) ; new point further west ; Zichron Yaakov ; Tel-Aviv ; Herzilya ; Netanya (average of Netanya and Tel-Aviv, and over the sea) ; lat_beg = 32.08 lat_beg = 32.15 lon_beg = 34.60 i_loc = 1 j_loc = 1 GET_IJ::get_ij(xlat_zoom,xlon_zoom,lat_beg,lon_beg,i_loc,j_loc,dims_zoom(0),dims_zoom(1)) i_point = i_loc-1 j_point = j_loc-1 print("i_point = " + i_point) print("j_point = " + j_point) print(xlat(j_point,i_point)) print(xlon(j_point,i_point)) ; We generate plots, but what kind do we prefer? ; type = "x11" ; type = "png" ; wtype = "pdf" type = "pdf" type@wkWidth = 2500 ; produces a better looking PNG type@wkHeight = 2500 ; type = "ps" ; type = "ncgm" ; wks = gsn_open_wks(type,"elec_vert_cross_sections_45_deg_01_05_2018_mar_exp_decay") ; wks = gsn_open_wks(type,"elec_vert_cross_sections_dusty") ; wks = gsn_open_wks(type,"radar_lpi_vert_cross_sections_1000") ; wks = gsn_open_wks(type,"radar_lpi_vert_cross_sections_0700") ; wks = gsn_open_wks(type,"radar_lpi_vert_cross_sections_0800") wks = gsn_open_wks(type,"radar_lpi_vert_cross_sections_1000") ; wks = gsn_open_wks(type,"radar_lpi_vert_cross_sections_0530") ; wks = gsn_open_wks(type,"radar_lpi_vert_cross_sections_1030") ; wks = gsn_open_wks(type,"radar_lpi_vert_cross_sections_0700") ; 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_res1_list() opts_kdp = set_res1_list() ; ;--------------------------------------------------------------- ; do it = 0,ntimes-1,2 ; TIME LOOP ; do it = 219,219,1 ; do it = 220,220,1 ; do it = 232,232,1 ; do it = 220,220,1 ; do it = 224,224,1 do it = 232,232,1 ; do it = 214,214,1 ; do it = 234,234,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 UTC" ; format = "%d %c %Y" new_time_array = new_time_array + 0.01 time_var=cd_string(new_time_array,format) print("time_var = " + time_var) 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 UTC" ; format = "%d %c %Y" new_time_array = new_time_array + 0.01 time_var_1=cd_string(new_time_array,format) print("time_var_1 = " + time_var_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 UTC" ; format = "%d %c %Y" new_time_array = new_time_array + 0.01 time_var_2=cd_string(new_time_array,format) print("time_var_2 = " + time_var_2) print("Working on time: " + times(it) ) ; res@TimeLabel = times(it) ; Set Valid time to use on plots ; rh = wrf_user_getvar(a,"rh",it) ; relative humidity ; fft = wrf_user_getvar(a,"LIGHTNING_NEU",it) ; total reflectivity ; printVarSummary(fft) ; ff1 = wrf_user_getvar(a,"LIGHTNING_NEU",it+2) ; liquid ; ff2 = wrf_user_getvar(b,"LIGHTNING_NEU",it) ; total reflectivity ; ff3 = wrf_user_getvar(b,"LIGHTNING_NEU",it+2) ; liquid ; ff2 = wrf_user_getvar(a,"INDUC",it) ; plates ; ff3 = wrf_user_getvar(a,"NONINDUC",it) ; plates ; fft = fft*1.E-6 ; ff1 = ff1*1.E-6 ; ff2 = ff2*1.E-6 ; ff3 = ff3*1.E-6 time = it rad_h= wrf_user_getvar(a,"Total_zh",time) lpi_1 = rad_h lpi := wrf_user_getvar(a,"QGRAUP",time) ph := wrf_user_getvar(a,"PH",time) phb := wrf_user_getvar(a,"PHB",time) p := wrf_user_getvar(a,"P",time) pb := wrf_user_getvar(a,"PB",time) qc := wrf_user_getvar(a,"QCLOUD",time) qr := wrf_user_getvar(a,"QRAIN",time) qi := wrf_user_getvar(a,"QICE",time) qs := wrf_user_getvar(a,"QSNOW",time) qg := wrf_user_getvar(a,"QGRAUP",time) w := wrf_user_getvar(a,"W",time) t := wrf_user_getvar(a,"T",time) u_a := wrf_user_getvar(a,"ua",time) v_a := wrf_user_getvar(a,"va",time) dims3d := dimsizes(qg) dims3dw := dimsizes(w) dims3dph := dimsizes(phb) ni = dims3d(2) nj = dims3d(1) nk = dims3d(0) nkk = dims3dph(0) CALC_LPI::calclpi(qc, qr, qi, qs, qg,w,ph,phb,p,pb,t,lpi,ni,nj,nk,nkk) lpi_2 = lpi time = it rad_h= wrf_user_getvar(b,"Total_zh",time) lpi_3 = rad_h ph := wrf_user_getvar(b,"PH",time) phb := wrf_user_getvar(b,"PHB",time) p := wrf_user_getvar(b,"P",time) pb := wrf_user_getvar(b,"PB",time) qc := wrf_user_getvar(b,"QCLOUD",time) qr := wrf_user_getvar(b,"QRAIN",time) qi := wrf_user_getvar(b,"QICE",time) qs := wrf_user_getvar(b,"QSNOW",time) qg := wrf_user_getvar(b,"QGRAUP",time) w := wrf_user_getvar(b,"W",time) t := wrf_user_getvar(b,"T",time) u_b := wrf_user_getvar(b,"ua",time) v_b := wrf_user_getvar(b,"va",time) dims3d := dimsizes(qg) dims3d := dimsizes(qg) dims3dw := dimsizes(w) dims3dph := dimsizes(phb) ni = dims3d(2) nj = dims3d(1) nk = dims3d(0) nkk = dims3dph(0) CALC_LPI::calclpi(qc, qr, qi, qs, qg,w,ph,phb,p,pb,t,lpi,ni,nj,nk,nkk) lpi_4 = lpi rad_h= wrf_user_getvar(c,"Total_zh",time) lpi_5 = rad_h ph := wrf_user_getvar(c,"PH",time) phb := wrf_user_getvar(c,"PHB",time) p := wrf_user_getvar(c,"P",time) pb := wrf_user_getvar(c,"PB",time) qc := wrf_user_getvar(c,"QCLOUD",time) qr := wrf_user_getvar(c,"QRAIN",time) qi := wrf_user_getvar(c,"QICE",time) qs := wrf_user_getvar(c,"QSNOW",time) qg := wrf_user_getvar(c,"QGRAUP",time) w := wrf_user_getvar(c,"W",time) t := wrf_user_getvar(c,"T",time) u_c := wrf_user_getvar(c,"ua",time) v_c := wrf_user_getvar(c,"va",time) dims3d := dimsizes(qg) dims3dw := dimsizes(w) dims3dph := dimsizes(phb) ni = dims3d(2) nj = dims3d(1) nk = dims3d(0) nkk = dims3dph(0) CALC_LPI::calclpi(qc, qr, qi, qs, qg,w,ph,phb,p,pb,t,lpi,ni,nj,nk,nkk) lpi_6 = lpi tc_a = wrf_user_getvar(a,"tc",it) ; T in C z_a = wrf_user_getvar(a, "z",it) ; grid point height tc_b = wrf_user_getvar(b,"tc",it) ; T in C z_b = wrf_user_getvar(b, "z",it) ; grid point height tc_c = wrf_user_getvar(c,"tc",it) ; T in C z_c = wrf_user_getvar(c, "z",it) ; grid point height fft = lpi_1 ff1 = lpi_2 ff2 = lpi_3 ff3 = lpi_4 ff4 = lpi_5 ff5 = lpi_6 printMinMax(fft,False) printMinMax(ff1,False) printMinMax(ff2,False) printMinMax(ff3,False) ; print(x) tc_zoom_a = tc_a(:,x_start:x_end,y_start:y_end) tc_zoom_b = tc_b(:,x_start:x_end,y_start:y_end) tc_zoom_c = tc_c(:,x_start:x_end,y_start:y_end) u_zoom_a = u_a(:,x_start:x_end,y_start:y_end) u_zoom_b = u_b(:,x_start:x_end,y_start:y_end) u_zoom_c = u_c(:,x_start:x_end,y_start:y_end) v_zoom_a = v_a(:,x_start:x_end,y_start:y_end) v_zoom_b = v_b(:,x_start:x_end,y_start:y_end) v_zoom_c = v_c(:,x_start:x_end,y_start:y_end) fft_zoom = fft(:,x_start:x_end,y_start:y_end) ff1_zoom = ff1(:,x_start:x_end,y_start:y_end) ff2_zoom = ff2(:,x_start:x_end,y_start:y_end) ff3_zoom = ff3(:,x_start:x_end,y_start:y_end) ff4_zoom = ff4(:,x_start:x_end,y_start:y_end) ff5_zoom = ff5(:,x_start:x_end,y_start:y_end) z_zoom_a = z_a(:,x_start:x_end,y_start:y_end) z_zoom_b = z_b(:,x_start:x_end,y_start:y_end) z_zoom_c = z_c(:,x_start:x_end,y_start:y_end) printMinMax(fft_zoom,False) printMinMax(ff1_zoom,False) printMinMax(ff2_zoom,False) printMinMax(ff3_zoom,False) ;---Convert to 1D a1D = ndtooned(fft_zoom) dsizes_a = dimsizes(fft_zoom) ;---Resolve the 1D indices back to their original 3D array. indices = ind_resolve(maxind(a1D),dsizes_a) print(indices) indices_t = indices print(fft_zoom(indices(0,0),indices(0,1),indices(0,2))) ;---Convert to 1D a1D = ndtooned(ff1_zoom) dsizes_a = dimsizes(ff1_zoom) ;---Resolve the 1D indices back to their original 3D array. indices = ind_resolve(maxind(a1D),dsizes_a) print(indices) print(ff1_zoom(indices(0,0),indices(0,1),indices(0,2))) indices_1 = indices ;---Convert to 1D a1D = ndtooned(ff2_zoom) dsizes_a = dimsizes(ff2_zoom) ;---Resolve the 1D indices back to their original 3D array. indices = ind_resolve(maxind(a1D),dsizes_a) print(indices) indices_2 = indices print(ff2_zoom(indices(0,0),indices(0,1),indices(0,2))) ;---Convert to 1D a1D = ndtooned(ff3_zoom) dsizes_a = dimsizes(ff3_zoom) ;---Resolve the 1D indices back to their original 3D array. indices = ind_resolve(maxind(a1D),dsizes_a) print(indices) print(ff3_zoom(indices(0,0),indices(0,1),indices(0,2))) indices_3 = indices ;---Convert to 1D a1D = ndtooned(ff4_zoom) dsizes_a = dimsizes(ff4_zoom) ;---Resolve the 1D indices back to their original 3D array. indices = ind_resolve(maxind(a1D),dsizes_a) print(indices) print(ff4_zoom(indices(0,0),indices(0,1),indices(0,2))) indices_4 = indices ;---Convert to 1D a1D = ndtooned(ff5_zoom) dsizes_a = dimsizes(ff5_zoom) ;---Resolve the 1D indices back to their original 3D array. indices = ind_resolve(maxind(a1D),dsizes_a) print(indices) print(ff5_zoom(indices(0,0),indices(0,1),indices(0,2))) indices_5 = indices ;---Convert to 1D if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 12. ; We are only interested in the first 6km nz = floattoint(zmax + 1) end if ;--------------------------------------------------------------- ; do ip = 1, 3 ; we are doing 3 plots do ip = 3, 3 ; we are doing 3 plots ; do ip = 1, 1 ; we are doing 3 plots ; all with the pivot point (plane) in the center of the domain ; at angles 0, 45 and 90 ; ; | ; angle=0 is | ; | ; plane = new(2,float) ; note this part of the code assumes that the arrays have been created (x,y) ; plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y) plane = (/ i_point, j_point /) ; pivot point is Upton opts = False if(ip .eq. 1) then angle = 90. X_plane = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) X_desc = "Longitude" end if if(ip .eq. 2) then angle = 0. X_plane = wrf_user_intrp2d(xlat_zoom,plane,angle,opts) X_desc = "latitude" end if if(ip .eq. 3) then ; angle = 45. ; angle = 10. angle = 90. ; West to East! ; angle = 135. X_desc = "Longitude" end if plane = (/ indices_t(0,2), indices_t(0,1) /) ; pivot point is Upton X_plane_t = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) fft_plane = wrf_user_intrp3d(fft_zoom,z_zoom_a,"v",plane,angle,opts) tc_plane_t = wrf_user_intrp3d(tc_zoom_a,z_zoom_a,"v",plane,angle,opts) u_plane_a = wrf_user_intrp3d(u_zoom_a,z_zoom_a,"v",plane,angle,opts) v_plane_a = wrf_user_intrp3d(v_zoom_a,z_zoom_a,"v",plane,angle,opts) plane = (/ indices_1(0,2), indices_t(0,1) /) ; pivot point is Upton ; plane = (/ indices_1(0,2), indices_1(0,1) /) ; pivot point is Upton X_plane_1 = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) ff1_plane = wrf_user_intrp3d(ff1_zoom,z_zoom_a,"v",plane,angle,opts) tc_plane_1 = wrf_user_intrp3d(tc_zoom_a,z_zoom_a,"v",plane,angle,opts) plane = (/ indices_2(0,2), indices_2(0,1) /) ; pivot point is Upton X_plane_2 = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) ff2_plane = wrf_user_intrp3d(ff2_zoom,z_zoom_b,"v",plane,angle,opts) tc_plane_2 = wrf_user_intrp3d(tc_zoom_b,z_zoom_b,"v",plane,angle,opts) u_plane_b = wrf_user_intrp3d(u_zoom_b,z_zoom_b,"v",plane,angle,opts) v_plane_b = wrf_user_intrp3d(v_zoom_b,z_zoom_b,"v",plane,angle,opts) ; plane = (/ indices_3(0,2), indices_3(0,1) /) ; pivot point is Upton plane = (/ indices_3(0,2), indices_3(0,1) /) ; pivot point is Upton X_plane_3 = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) ff3_plane = wrf_user_intrp3d(ff3_zoom,z_zoom_b,"v",plane,angle,opts) tc_plane_3 = wrf_user_intrp3d(tc_zoom_b,z_zoom_b,"v",plane,angle,opts) printVarSummary(fft_plane) plane = (/ indices_4(0,2), indices_4(0,1) /) ; pivot point is Upton X_plane_4 = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) ff4_plane = wrf_user_intrp3d(ff4_zoom,z_zoom_b,"v",plane,angle,opts) tc_plane_4 = wrf_user_intrp3d(tc_zoom_c,z_zoom_c,"v",plane,angle,opts) u_plane_c = wrf_user_intrp3d(u_zoom_c,z_zoom_c,"v",plane,angle,opts) v_plane_c = wrf_user_intrp3d(v_zoom_c,z_zoom_c,"v",plane,angle,opts) printVarSummary(ff4_plane) plane = (/ indices_5(0,2), indices_5(0,1) /) ; pivot point is Upton X_plane_5 = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) ff5_plane = wrf_user_intrp3d(ff5_zoom,z_zoom_c,"v",plane,angle,opts) tc_plane_5 = wrf_user_intrp3d(tc_zoom_c,z_zoom_c,"v",plane,angle,opts) printVarSummary(ff5_plane) ; Find the index where 6km is - only need to do this once if ( FirstTime ) then zz_a = wrf_user_intrp3d(z_zoom_a,z_a,"v",plane,angle,opts) zz_b = wrf_user_intrp3d(z_zoom_b,z_b,"v",plane,angle,opts) b_a = ind(zz_a(:,0) .gt. zmax*1000. ) b_b = ind(zz_b(:,0) .gt. zmax*1000. ) zmax_pos = b_a(0) - 1 zmax_pos = b_a(0) - 1 if ( abs(zz_a(zmax_pos,0)-zmax*1000.) .lt. abs(zz_a(zmax_pos+1,0)-zmax*1000.) ) then zspan = b_a(0) - 1 else zspan = b_a(0) end if delete(zz_a) delete(b_a) FirstTime = False end if ; X-axis lables X_plane =X_plane_t print("x_plane = " + X_plane) dimsX = dimsizes(X_plane) print("dimsX = " + dimsX) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 3) xspan_t = xspan nx_t = nx xmax_t = xmax xmin_t = xmin print("xmax_t = " + xmax_t) print("xmin_t = " + xmin_t) ; print(x) X_plane :=X_plane_1 dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 3) xspan_1 = xspan nx_1 = nx xmax_1 = xmax xmin_1 = xmin X_plane :=X_plane_2 dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 3) xspan_2 = xspan nx_2 = nx xmax_2 = xmax xmin_2 = xmin X_plane :=X_plane_3 dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 3) xspan_3 = xspan nx_3 = nx xmax_3 = xmax xmin_3 = xmin X_plane :=X_plane_4 dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 3) xspan_4 = xspan nx_4 = nx xmax_4 = xmax xmin_4 = xmin X_plane :=X_plane_5 dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 3) xspan_5 = xspan nx_5 = nx xmax_5 = xmax xmin_5 = xmin ; ; xmin_new = nx-xmin ; xmax_new = nx+xmin ; print("xmin_new = " + xmin_new) ; print("xmax_new = " + xmax_new) ;--------------------------------------------------------------- ; Options for XY Plots opts_xy = True 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" print("opts_xy@tmXBLabels" + opts_xy@tmXBLabels) 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_t = opts_xy opts_xy_1 = opts_xy opts_xy_2 = opts_xy opts_xy_3 = opts_xy opts_xy_4 = opts_xy opts_xy_5 = opts_xy opts_xy_t@tmXBValues = fspan(0,xspan_t,nx_t) ; Create tick marks opts_xy_t@tmXBLabels = sprintf("%.2f",fspan(xmin_t,xmax_t,nx_t)) ; Create labels opts_xy_1@tmXBValues = fspan(0,xspan_1,nx_1) ; Create tick marks opts_xy_1@tmXBLabels = sprintf("%.2f",fspan(xmin_t,xmax_1,nx_1)) ; Create labels opts_xy_2@tmXBValues = fspan(0,xspan_2,nx_2) ; Create tick marks opts_xy_2@tmXBLabels = sprintf("%.2f",fspan(xmin_t,xmax_2,nx_2)) ; Create labels opts_xy_3@tmXBValues = fspan(0,xspan_3,nx_3) ; Create tick marks opts_xy_3@tmXBLabels = sprintf("%.2f",fspan(xmin_t,xmax_3,nx_3)) ; Create labels opts_xy_4@tmXBValues = fspan(0,xspan_4,nx_4) ; Create tick marks opts_xy_4@tmXBLabels = sprintf("%.2f",fspan(xmin_4,xmax_4,nx_4)) ; Create labels opts_xy_5@tmXBValues = fspan(0,xspan_5,nx_5) ; Create tick marks opts_xy_5@tmXBLabels = sprintf("%.2f",fspan(xmin_5,xmax_5,nx_5)) ; Create labels ; opts_xy@PlotOrientation = tc_plane@Orientation ; opts_xy@PlotOrientation = tc_plane@Orientation ; opts_xy@PlotOrientation = tc_plane@Orientation ; Plotting options for RH ;;; opts_xy opts_rad@gsnLeftStringFontHeightF = 0.03 opts_rad@gsnRightStringFontHeightF = 0.03 opts_rad@tiXAxisString = X_desc opts_rad@tiYAxisString = "Height (km)" opts_rad@cnMissingValPerimOn = True opts_rad@cnMissingValFillColor = 0 opts_rad@cnMissingValFillPattern = 11 opts_rad@tmXTOn = False opts_rad@tmYROn = False opts_rad@tmXBMode = "Explicit" print("opts_rad@tmXBLabels " + opts_rad@tmXBLabels) opts_rad@tmXBLabelFontHeightF = 0.015 opts_rad@tmYLMode = "Explicit" opts_rad@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_rad@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels opts_rad@tiXAxisFontHeightF = 0.040 opts_rad@tiYAxisFontHeightF = 0.040 opts_rad@tmXBMajorLengthF = 0.02 opts_rad@tmYLMajorLengthF = 0.02 opts_rad@tmYLLabelFontHeightF = 0.035 opts_rad@tmXBLabelFontHeightF = 0.035 opts_rad@cnInfoLabelOn = False ; don't need four, just 2 of these for opts_rad opts_rad_t = opts_rad opts_rad_1 = opts_rad opts_rad_2 = opts_rad opts_rad_3 = opts_rad opts_rad_4 = opts_rad opts_rad_5 = opts_rad opts_rad_t@tmXBValues = fspan(0,xspan_t,nx_t) ; Create tick marks opts_rad_t@tmXBLabels = sprintf("%.1f",fspan(xmin_t,xmax_t,nx_t)) ; Create labels opts_rad_1@tmXBValues = fspan(0,xspan_1,nx_1) ; Create tick marks opts_rad_1@tmXBLabels = sprintf("%.1f",fspan(xmin_1,xmax_1,nx_1)) ; Create labels opts_rad_2@tmXBValues = fspan(0,xspan_2,nx_2) ; Create tick marks opts_rad_2@tmXBLabels = sprintf("%.1f",fspan(xmin_2,xmax_2,nx_2)) ; Create labels opts_rad_3@tmXBValues = fspan(0,xspan_3,nx_3) ; Create tick marks opts_rad_3@tmXBLabels = sprintf("%.1f",fspan(xmin_3,xmax_3,nx_3)) ; Create labels opts_rad_4@tmXBValues = fspan(0,xspan_4,nx_4) ; Create tick marks opts_rad_4@tmXBLabels = sprintf("%.1f",fspan(xmin_4,xmax_4,nx_4)) ; Create labels opts_rad_5@tmXBValues = fspan(0,xspan_5,nx_5) ; Create tick marks opts_rad_5@tmXBLabels = sprintf("%.1f",fspan(xmin_5,xmax_5,nx_5)) ; Create labels opts_zdr@gsnLeftStringFontHeightF = 0.03 opts_zdr@gsnRightStringFontHeightF = 0.03 opts_zdr@tiXAxisString = X_desc opts_zdr@tiYAxisString = "Height (km)" opts_zdr@cnMissingValPerimOn = True opts_zdr@cnMissingValFillColor = 0 opts_zdr@cnMissingValFillPattern = 11 opts_zdr@tmXTOn = False opts_zdr@tmYROn = False opts_zdr@tmXBMode = "Explicit" opts_zdr@cnInfoLabelOn = False ; opts_zdr@tmXBValues = fspan(0,xspan,nx) ; Create tick marks ; opts_zdr@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels print("opts_zdr@tmXBLabels" + opts_zdr@tmXBLabels) opts_zdr@tmXBLabelFontHeightF = 0.015 opts_zdr@tmYLMode = "Explicit" opts_zdr@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_zdr@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels opts_zdr@tiXAxisFontHeightF = 0.040 opts_zdr@tiYAxisFontHeightF = 0.040 opts_zdr@tmXBMajorLengthF = 0.02 opts_zdr@tmYLMajorLengthF = 0.02 opts_zdr@tmYLLabelFontHeightF = 0.035 opts_zdr@tmXBLabelFontHeightF = 0.035 ; opts_zdr@PlotOrientation = tc_plane@Orientation opts_zdr_t = opts_zdr opts_zdr_1 = opts_zdr opts_zdr_2 = opts_zdr opts_zdr_3 = opts_zdr opts_zdr_4 = opts_zdr opts_zdr_5 = opts_zdr opts_zdr_t@tmXBValues = fspan(0,xspan_t,nx_t) ; Create tick marks opts_zdr_t@tmXBLabels = sprintf("%.1f",fspan(xmin_t,xmax_t,nx_t)) ; Create labels opts_zdr_1@tmXBValues = fspan(0,xspan_1,nx_1) ; Create tick marks opts_zdr_1@tmXBLabels = sprintf("%.1f",fspan(xmin_1,xmax_1,nx_1)) ; Create labels opts_zdr_1 = opts_zdr_t opts_zdr_2@tmXBValues = fspan(0,xspan_2,nx_2) ; Create tick marks opts_zdr_2@tmXBLabels = sprintf("%.1f",fspan(xmin_2,xmax_2,nx_2)) ; Create labels opts_zdr_3@tmXBValues = fspan(0,xspan_3,nx_3) ; Create tick marks opts_zdr_3@tmXBLabels = sprintf("%.1f",fspan(xmin_3,xmax_3,nx_3)) ; Create labels opts_zdr_4@tmXBValues = fspan(0,xspan_4,nx_4) ; Create tick marks opts_zdr_4@tmXBLabels = sprintf("%.1f",fspan(xmin_4,xmax_4,nx_4)) ; Create labels opts_zdr_5@tmXBValues = fspan(0,xspan_5,nx_5) ; Create tick marks opts_zdr_5@tmXBLabels = sprintf("%.1f",fspan(xmin_5,xmax_5,nx_5)) ; Create labels opts_crs@tiXAxisString = X_desc opts_crs@tiYAxisString = "Height (km)" opts_crs@cnMissingValPerimOn = True opts_crs@cnMissingValFillColor = 0 opts_crs@cnMissingValFillPattern = 11 opts_crs@tmXTOn = False opts_crs@tmYROn = False opts_crs@tmXBMode = "Explicit" opts_crs@tmXBValues = fspan(0,xspan,nx) ; Create tick marks opts_crs@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels print("opts_crs@tmXBLabels" + opts_crs@tmXBLabels) opts_crs@tmXBLabelFontHeightF = 0.015 opts_crs@tmYLMode = "Explicit" opts_crs@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_crs@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels opts_crs@tiXAxisFontHeightF = 0.020 opts_crs@tiYAxisFontHeightF = 0.020 opts_crs@tmXBMajorLengthF = 0.02 opts_crs@tmYLMajorLengthF = 0.02 opts_crs@tmYLLabelFontHeightF = 0.015 ; opts_crs@PlotOrientation = tc_plane@Orientation opts_kdp@tiXAxisString = X_desc opts_kdp@tiYAxisString = "Height (km)" opts_kdp@cnMissingValPerimOn = True opts_kdp@cnMissingValFillColor = 0 opts_kdp@cnMissingValFillPattern = 11 opts_kdp@tmXTOn = False opts_kdp@tmYROn = False opts_kdp@tmXBMode = "Explicit" opts_kdp@tmXBValues = fspan(0,xspan,nx) ; Create tick marks opts_kdp@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels print("opts_kdp@tmXBLabels" + opts_kdp@tmXBLabels) opts_kdp@tmXBLabelFontHeightF = 0.015 opts_kdp@tmYLMode = "Explicit" opts_kdp@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_kdp@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels opts_kdp@tiXAxisFontHeightF = 0.020 opts_kdp@tiYAxisFontHeightF = 0.020 opts_kdp@tmXBMajorLengthF = 0.02 opts_kdp@tmYLMajorLengthF = 0.02 opts_kdp@tmYLLabelFontHeightF = 0.015 ; opts_kdp@PlotOrientation = tc_plane@Orientation ; Plotting options for Temperature opts_tc = opts_xy_1 opts_tc@cnInfoLabelZone = 1 opts_tc@cnInfoLabelSide = "Top" opts_tc@cnInfoLabelPerimOn = True opts_tc@cnInfoLabelOrthogonalPosF = -0.00005 ; opts_tc@ContourParameters = (/ 5. /) opts_tc@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels opts_tc@cnLevels = ispan(-60,40,4) ; set the contour levels opts_tc@gsnDraw = False ; don't draw opts_tc@gsnFrame = False ; don't advance frame opts_tc@cnInfoLabelOn = False opts_tc@gsnContourZeroLineThicknessF = 3.5 vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = 20 ; define vector ref mag vecres@vcRefLengthF = 0.045 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " ; turn off axis label vecres@vcRefAnnoOrthogonalPosF = -.535 ; move ref vector into plot ; Get the contour info for the rh and temp ; plot_tc1 = gsn_contour(wks,vect_plane_a(0:zmax_pos,:),opts_tc) opts_tc = opts_xy_1 plot_tc2 = gsn_contour(wks,tc_plane_1(0:zmax_pos,:),opts_tc) opts_tc = opts_xy_2 ; plot_tc3 = gsn_contour(wks,vect_plane_b(0:zmax_pos,:),opts_tc) opts_tc = opts_xy_3 plot_tc4 = gsn_contour(wks,tc_plane_3(0:zmax_pos,:),opts_tc) opts_tc = opts_xy_4 ; plot_tc5 = gsn_contour(wks,vect_plane_c(0:zmax_pos,:),opts_tc) opts_tc = opts_xy_5 plot_tc6 = gsn_contour(wks,tc_plane_5(0:zmax_pos,:),opts_tc) plot_a = gsn_vector(wks,u_plane_a(0:zmax_pos,::4),v_plane_a(0:zmax_pos,::4),vecres) plot_b = gsn_vector(wks,u_plane_b(0:zmax_pos,::4),v_plane_b(0:zmax_pos,::4),vecres) plot_c = gsn_vector(wks,u_plane_c(0:zmax_pos,::4),v_plane_c(0:zmax_pos,::4),vecres) ; plot_c = gsn_vector(wks,u_plane_c(::4,::4),v_plane_c(::4,::4),vecres) opts_rad = opts_rad_t opts_rad@gsnRightString = time_var_1 opts_rad@gsnLeftString = "Control" opts_rad@cnLevels := ispan(15,65,5) opts_rad@lbTitleString = "dBZ" opts_zdr@lbTitleString = "J/kg" plot_rad1(0) = gsn_csm_contour(wks,fft_plane(0:zmax_pos,:),opts_rad) opts_zdr = opts_zdr_1 opts_zdr@cnLevels = ispan(10,610,60) opts_zdr@gsnRightString = time_var_2 opts_zdr@gsnLeftString = "Control" opts_zdr@cnInfoLabelOn = False plot_rad1(1) = gsn_csm_contour(wks,ff1_plane(0:zmax_pos,:),opts_zdr) opts_rad = opts_rad_2 opts_rad@lbTitleString = "dBZ" opts_rad@cnLevels := ispan(15,65,5) opts_rad@gsnRightString = time_var_1 opts_rad@gsnLeftString = "Control + Desert Dust" plot_rad1(2) = gsn_csm_contour(wks,ff2_plane(0:zmax_pos,:),opts_rad) opts_zdr = opts_zdr_3 opts_zdr@lbTitleString = "J/kg" opts_zdr@gsnRightString = time_var_2 opts_zdr@gsnLeftString = "Control + Desert Dust" opts_zdr@cnLevels = ispan(10,610,60) plot_rad1(3) = gsn_csm_contour(wks,ff3_plane(0:zmax_pos,:),opts_zdr) opts_rad = opts_rad_2 opts_rad@lbTitleString = "dBZ" opts_rad@cnLevels := ispan(15,65,5) opts_rad@gsnRightString = time_var_1 opts_rad@gsnLeftString = "Control +" plot_rad1(4) = gsn_csm_contour(wks,ff4_plane(0:zmax_pos,:),opts_rad) opts_zdr = opts_zdr_3 opts_zdr@lbTitleString = "J/kg" opts_zdr@gsnRightString = time_var_2 opts_zdr@gsnLeftString = "Control +" opts_zdr@cnLevels = ispan(10,610,60) plot_rad1(5) = gsn_csm_contour(wks,ff5_plane(0:zmax_pos,:),opts_zdr) ;--------------------------------------------------------------- ; 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 = time_var pnlres@gsnPanelYWhiteSpacePercent = 5 ; Add white space b/w plots. pnlres@gsnPanelLabelBar = False ; Turn on common labelbar pnlres@lbLabelAutoStride = False ; 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@pmLabelBarWidthF = 0.5 pnlres@pmLabelBarHeightF = 0.15 pnlres@lbTopMarginF = 0.001 pnlres@pmLabelBarOrthogonalPosF = 0.02 ; position wrt plot pnlres@lbTitleOffsetF = -0.4 pnlres@amJust = "TopLeft" pnlres@gsnPanelFigureStringsFontHeightF = 0.0175 pnlres@gsnPanelFigureStrings = (/"(a)","(b)","(c)","(d)","(e)","(f)"/) overlay(plot_rad1(0),plot_a) overlay(plot_rad1(1),plot_tc2) overlay(plot_rad1(2),plot_b) overlay(plot_rad1(3),plot_tc4) overlay(plot_rad1(4),plot_c) overlay(plot_rad1(5),plot_tc6) gsn_panel(wks,(/plot_rad1/),(/3,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) end do ; make next cross section ; Delete options and fields, so we don't have carry over end do ; END OF TIME LOOP end