;******************************************************************** ; rcm_3.ncl ;******************************************* ; rcm_3.ncl ; ; Concepts illustrated: ; - Plotting RCM precipitation data ; - Drawing color-filled contours over a lambert conformal map ; - Using a filter to convert Cray COS block binary data to netCDF ;******************************************************************** ; ; These files are loaded by default in NCL V6.2.0 and newer ; 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/contrib/cd_string.ncl" external GET_IJ "./get_ij.so" begin ;plots = new ( 20, graphic ) ;plots = new ( 16, graphic ) plots = new ( 12, graphic ) plots_b = new ( 6, graphic ) ;************************* ; read in precip data (units cm) ;************************* ; a = addfile("prc.1982112500-1983030500.nc","r") lat_beg = 31.0 lat_end = 40.0 lon_beg = -104. ; must use real or you get warning about coerced to the appropriate type lon_end = -95. ; filename_in="../RUN/2KM_ALL_CHARGE/wrfout_d01_2012-04-15_00:00:00" ; filename_in="../RUN/2KM_NSSL_SAU/wrfout_d02_2012-04-15_00:00:00" ; filename_in="../RUN/2KM_NSSL_SAU/wrfout_d02_2012-04-15_00:00:00" ; filename_in="../RUN/wrfout_d02_2012-04-15_00:00:00_saund_2con" ; filename_in="../RUN/wrfout_d02_2012-04-15_00:00:00" ; myFiles_a = systemfunc("ls ../RUN/2KM_CON/wrfout*d02**") myFiles_a = systemfunc("ls ../RUN/2KM_NSSL/wrfout*d02**") myFiles_b = systemfunc("ls ../RUN/wrfout*d02**") myFiles_c = systemfunc("ls /data22/wrf/APRIL_11_13_2023_DUST/RUN/wrfout*d02**") ; myFiles_b = systemfunc("ls /data22/wrf/APRIL_11_13_2023_DUST/RUN/wrfout*d02**") ; myFiles_c = systemfunc("ls ../RUN/1.3KM_2CON/wrfout*d02**") ; myFiles_d = systemfunc("ls ../RUN/1.3KM_NSSL/wrfout*d02**") n_files = dimsizes(myFiles_a) print(myFiles_a) ; print(myFiles_b) print(n_files) filename_in = myFiles_a(0) a_in = addfile(filename_in,"r") ; print(a_in) ; print(a_in) ; print(x) xlat2d=wrf_user_getvar(a_in,"XLAT",0) xlon2d=wrf_user_getvar(a_in,"XLONG",0) dims2d=dimsizes(xlat2d) lat_beg = xlat2d(0,0)+2.25 ; lat_end = xlat2d(dims2d(0)-1,dims2d(1)-1)-0.25 lat_end = 37.5 lon_beg = xlon2d(0,0)+1.25 lon_end = xlon2d(dims2d(0)-1,dims2d(1)-1)-0.25 dims2d = dimsizes(xlat2d) i_loc = 0 j_loc = 0 GET_IJ::get_ij(xlat2d,xlon2d,lat_beg,lon_beg,i_loc,j_loc,dims2d(0),dims2d(1)) i_start = i_loc-1 j_start = j_loc-1 GET_IJ::get_ij(xlat2d,xlon2d,lat_end,lon_end,i_loc,j_loc,dims2d(0),dims2d(1)) i_end = i_loc-1 j_end = j_loc-1 ; wks = gsn_open_wks("pdf","mod_obs_map_2km_April_15_2012_nssl") ; send graphics to PNG file wks = gsn_open_wks("pdf","mod_radar_map_2km_April_11_13_2013") ; send graphics to PNG file ;normalized lightning flash origin density myLFiles = systemfunc("ls ../LIGHT/total*csv") print(myLFiles) no_files = dimsizes(myLFiles) print(no_files) i_file = 0 ii_file = 0 l_file = 0 ; do it = 24,24,4 ; do it = 3,ntimes-1,1 ; do it = 20,ntimes-1,4 ; do it = 72,216,36 ; do it = 72,144,36 ; do it = 36,36,1 do it = 54,54,1 ; do it = 36,360,1 ; do it = 30,30,6 ; do it = 24,24 print("it = "+it) it_beg = it-36 it_end = it a_in = addfile(myFiles_a(it_beg),"r") times = wrf_user_getvar(a_in,"times",-1) ; get times in the file ntimes = dimsizes(times) ; number of times in the file print(times) 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) new_time_units = "hours since 1800-01-01 00:00" new_time_units = new_time_units ; it_val = it_beg it_val = 0 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) ; 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_b=cd_string(new_time_array,format) ; it_val = it it_val = 0 a_in = addfile(myFiles_a(it_end),"r") times = wrf_user_getvar(a_in,"times",-1) ; get times in the file ntimes = dimsizes(times) ; number of times in the file print(times) 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) 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) 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_e=cd_string(new_time_array,format) LAT2D = xlat2d(j_start:j_end,i_start:i_end) LON2D = xlon2d(j_start:j_end,i_start:i_end) ; do i_file = 0,2,1 do i_file = 0,0,1 it_beg = it-36 it_end = it do itt = 0,0,1 if (i_file.le.2)then if (i_file.eq.0)then a_in= addfile(myFiles_a(it-itt),"r") print("myFiles_a = " + myFiles_a(it-itt)) radar = wrf_user_getvar(a_in, "REFL_10CM",0) LICN = dim_max_n(radar,0) end if if (i_file.eq.1)then a_in= addfile(myFiles_b(it-itt),"r") radar = wrf_user_getvar(a_in, "Total_zh",0) LICN = dim_max_n(radar,0) end if if (i_file.eq.2)then a_in= addfile(myFiles_c(it-itt),"r") radar = wrf_user_getvar(a_in, "Total_zh",0) LICN = dim_max_n(radar,0) end if ; if (i_file.eq.3)then ; a_in= addfile(myFiles_c(it-itt),"r") ; end if ; if (i_file.eq.4)then ; a_in= addfile(myFiles_d(it-itt),"r") ; end if ; printVarSummary(LICN) if (itt.eq.0)then total = LICN else total = LICN + total end if end if end do ; if (i_file.eq.4)then if (i_file.eq.3)then total = 0. do itt = 1,6 asci_file = myLFiles(l_file+itt) print("myLFiles = " + myLFiles(l_file+itt)) lines := asciiread(asci_file,-1,"string") n_lines = dimsizes(lines)-1 if (dimsizes(lines).eq.1)then continue else delim = "," str_array := asciiread(asci_file,-1,"string") values_2d := tofloat(str_split_csv(str_array,",",0)) dims_v = dimsizes(values_2d) printVarSummary(values_2d) Latitude := values_2d(1:dims_v(0)-1,6) Longitude := values_2d(1:dims_v(0)-1,7) Amps := 0.001*values_2d(1:dims_v(0)-1,9) type := values_2d(1:dims_v(0)-1,8) i_loc = 1 j_loc = 1 do i_line = 0,n_lines-1 if ((Latitude(i_line).ge.lat_beg.and.Latitude(i_line).le.lat_end) \ .and.(Longitude(i_line).ge.lon_beg.and.Longitude(i_line).le.lon_end))then GET_IJ::get_ij(xlat2d,xlon2d,Latitude(i_line),Longitude(i_line),i_loc,j_loc,dims2d(0),dims2d(1)) total(j_loc-1,i_loc-1) = total(j_loc-1,i_loc-1) + 1 end if end do end if end do l_file = l_file + 6 end if ; if (i_file.eq.0)then printMinMax(total,False) ; print(sum(total)) ; end if xx = total(j_start:j_end,i_start:i_end) ; xx = total(it:it,j_start:j_end,i_start:i_end) ; xx = total(j_start:j_end,i_start:i_end) ; x = dim_sum_n(xx,0) x = xx ; printVarSummary(x) ; don't plot any zeros. since this data doesn't actually contain any missing ; data, we can do the following to avoid plotting zero. x@_FillValue = 0.0 x@_FillValue = -999. ; LAT2D = a->lat2d ; LON2D = a->lon2d ;************************* ; set some parameters and read in date ; removed ;************************* ; index_start = ind(date.eq.date_start) ; get array index of desired times ; index_end = ind(date.eq.date_end ) ; nday = (index_end-index_start)/data_per_day ;*************************************** ; compute monthly average units: mm/day ;*************************************** dims = dimsizes(x) ; get dimension sizes nlat = dims(1) nlon = dims(0) ; y = new((/nlat,nlon/),float) ; allocate memory ; y = (/x(index_end,:,:) - x(index_start,:,:)/)*10./nday ; compute y = x y_bin = y dims_y_bin = dimsizes(y_bin) Levels = ispan(5,55,5) ; set the contour Levels ; Levels = (/1,3,5,7,9,11,13/) ; pnlres@lbLabelStride = 5 ; sets intervals between labels! ; Levels(0) = 0 ; Levels(1) = 1 n_Levels = dimsizes(Levels) num_markers = new((/n_Levels/),integer) num_markers = 0 printMinMax(y_bin,False) do i_level=0,n_Levels-1 ; print("i_level = " + i_level) i_marks = 0 do j = 0,dims_y_bin(0)-1 do i = 0,dims_y_bin(1)-1 check = ismissing(y_bin(j,i)) if (check.eq.False)then if (i_level.ne.n_Levels-1)then if (y_bin(j,i).ge.Levels(i_level).and.y_bin(j,i).lt.Levels(i_level+1))then num_markers(i_level) = num_markers(i_level) +1 i_marks = i_marks + 1 end if else if (y_bin(j,i).ge.Levels(i_level))then num_markers(i_level) = num_markers(i_level) +1 i_marks = i_marks + 1 end if end if end if end do end do ; print("num_markers above = " + i_level + " " + num_markers(i_level)) end do y = where(y.eq.0,y@_Fillvalue,y) y!0 ="lat" ; name dimensions y!1 ="lon" ;******************************************************* ; prepare data for ploting pull out all data except ; last value in x and y (bad points) ;******************************************************* ; precip = y(lat|0:nlat-2,lon|0:nlon-2) precip = y ; precip@long_name = "ENTLN Observations" ; give data long_name and units, if (i_file.le.5)then ; precip@long_name = "Lightning Flash Origin Density" ; give data long_name and units, precip@long_name = "SBM" ; give data long_name and units, precip@long_name = time_var_e ;title end if ; precip@long_name = "" ; precip@units = "(# h~S~-1~N~)" ; so plotted automatically. precip@units = "" ;*********************** ; plot ;*********************** ; wks = gsn_open_wks("png","rcm") ; send graphics to PNG file res = True ; plot mods desired res_colors = (/ \ "SpringGreen1", \ "SpringGreen4", \ "Yellow","Orange","Red1","HotPink","Purple"/) ; ; res@gsnDraw = False ; don't draw ; res@gsnFrame = False ; don't advance frame res@gsnDraw = True res@gsnFrame = True res@lbLabelBarOn = False ; res@cnInfoLabelOn = True ; !!!!! any plot of data that is on a native grid, must use the "corners" ; method of zooming in on map. ; res@cnFillOn = True res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit Levels res@cnLevels = 1.0*tofloat(Levels) ; x_span = fspan(xlon2d(0,0),xlon2d(dims2d(0)-1,dims2d(1)-1),7) ; print(x_span) ; res@mpXBMode = "Explicit" ; dimsX = dimsizes(x_span) ; nx = dimsX ; res@mpXBValues = fspan(0,dims2d(0),nx) ; res@mpXBLabels := sprintf("%.1f",x_span) ; Create labels ; print("res@mpXBValues" + res@mpXBValues) ; print("res@mpXBLabels" + res@mpXBLabels) ; res@mpXBLabelFontHeightF = 0.055 res@pmTickMarkDisplayMode = "Always" ; this insures latitude and longitude lines appear on x and y axis. ; res@mpLimitMode = "Corners" ; choose range of map res@mpLeftCornerLatF = lat_beg res@mpLeftCornerLonF = lon_beg res@mpRightCornerLatF = lat_end res@mpRightCornerLonF = lon_end ; The following 4 pieces of information are REQUIRED to properly display ; data on a native lambert conformal grid. This data should be specified ; somewhere in the model itself. ; WARNING: our local RCM users could not provide us with this information, ; so this is our best guess as to the correct numbers. Use at your own risk. res@mpProjection = "LambertConformal" ; res@mpProjection = "Mercator" res@mpLambertParallel1F = 30 res@mpLambertParallel2F = 58 res@mpLambertMeridianF = 260 ; ; Usually, when data is placed onto a map, it is transformed to the specified ; projection. Since this model is already on a native lambert conformal grid, ; we want to turn OFF the transformation. ; res@tfDoNDCOverlay = True ; do not transform ; res@tfDoNDCOverlay = "NDCViewport" ; NCL V6.5.0 or later res@cnFillOn = True ; color plot desired ; res@cnFillPalette = "gui_default" ; set color map ; res@cnFillPalette = "BlGrYeOrReVi200" ; set color map res@cnLinesOn = False ; no contour lines res@mpGeophysicalLineColor = "red" ; color of continental outlines res@mpPerimOn = True ; draw box around map res@mpGridLineDashPattern = 2 ; lat/lon lines as dashed res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; add state boundaries res@mpUSStateLineColor = "red" ; make them red ; res@tiMainOffsetYF = -0.005 res@tiMainPosition = "Center" res@tiMainString = time_var_b + " to " + time_var_e ;title if (i_file.eq.0)then ; precip@long_name = "Maritime" ; give data long_name and units, res@tiMainString = "NSSL" end if if (i_file.eq.1)then res@tiMainString = "SBM" end if if (i_file.eq.2)then res@tiMainString = "SBM-Dust" end if if (i_file.eq.3)then res@tiMainString = "Observations" end if if (i_file.eq.4)then ; precip@long_name = "Observed Lightning" ; give data long_name and units, res@tiMainString = "Observations" end if plot = gsn_csm_contour_map(wks,precip,res) yes_markers = False if (yes_markers.eq.True) mstring = "y" fontnum = 35 xoffset = 0.0 yoffset = 0.0 ratio = 1.0 size = 0.5 angle = 0.0 new_index = NhlNewMarker(wks, mstring, fontnum, xoffset, yoffset, ratio, size, angle) mkres_l = True ; Set up marker resources ; mkres_l@gsMarkerColor = "black" mkres_l@gsMarkerIndex = new_index mkres_l@gsMarkerSizeF = .003 ; Size of dot mkres_l@gsMarkerOpacityF = 0.15 do i_level=0,n_Levels-1 ;print("in loop " + i_level) ;print("num_markers(i_level) = " + num_markers(i_level)) if (num_markers(i_level).gt.0)then ; print("num_markers in loop + 1 = " +num_markers(i_level)) lat_l := new((/num_markers(i_level)/),float) lon_l := new((/num_markers(i_level)/),float) markid_l := new((/num_markers(i_level)/),graphic) i_mark = 0 do j = 0,dims_y_bin(0)-1 do i = 0,dims_y_bin(1)-1 check = ismissing(y_bin(j,i)) if (check.eq.False)then if (i_level.ne.n_Levels-1)then if (y_bin(j,i).ge.Levels(i_level).and.y_bin(j,i).lt.Levels(i_level+1))then lat_l(i_mark)=LAT2D(j,i) lon_l(i_mark)=LON2D(j,i) ; lat_l(i_mark)=LAT2D(i,j) ; lon_l(i_mark)=LON2D(i,j) i_mark = i_mark+1 end if else if (y_bin(j,i).gt.Levels(i_level))then lat_l(i_mark)=LAT2D(j,i) lon_l(i_mark)=LON2D(j,i) ; lat_l(i_mark)=LAT2D(i,j) ; lon_l(i_mark)=LON2D(i,j) i_mark = i_mark+1 end if end if end if end do end do ; print("i_mark = " + i_mark) mkres_l@gsMarkerColor = res_colors(i_level) ; Size of dot markid_l = gsn_add_polymarker(wks,plot,lon_l,lat_l,mkres_l) ; http://www.ncl.ucar.edu/Document/Graphics/Images/markers.png ; http://www.ncl.ucar.edu/Document/Graphics/marker_styles.shtml tmpstr = unique_string("marker") plot@$tmpstr$ = markid_l ; delete(lat_l) delete(lon_l) delete(markid_l) end if end do end if print("ii_file = " + ii_file) if (it.le.288)then plots(ii_file) = plot else plots_b(ii_file) = plot end if ii_file = ii_file + 1 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 = 1 ; Add white space b/w plots. pnlres@gsnPanelLabelBar = True ; Turn on common labelbar pnlres@lbLabelAutoStride = False ; Spacing of lbar labels. pnlres@lbBoxMinorExtentF = 0.13 pnlres@lbLabelStride = 2 pnlres@lbLabelFontHeightF = .015 pnlres@lbPerimOn = False pnlres@lbLabelFont = "Helvetica" ; label font pnlres@lbTitleOn = True ; turn on title pnlres@lbTitleString = "Maximum Reflectivity (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.55 pnlres@pmLabelBarHeightF = 0.15 pnlres@lbTopMarginF = 0.001 pnlres@pmLabelBarOrthogonalPosF = 0.0 ; position wrt plot ; pnlres@pmLabelBarOrthogonalPosF = 0.01 ; position wrt plot ; pnlres@pmLabelBarOrthogonalPosF = 0.02 ; closer to panels pnlres@pmLabelBarOrthogonalPosF = 0.028 ; closer to panels ; pnlres@pmLabelBarOrthogonalPosF = 0.03 ; closer to panels pnlres@lbTitleOffsetF = -0.10 pnlres@lbTitleOffsetF = -0.20;closer to bar ; pnlres@lbTitleOffsetF = -0.30;closer to bar pnlres@lbTitleOffsetF = -0.35;closer to bar ; pnlres@lbTitleOffsetF = -0.50;closer to bar pnlres@gsnPanelFigureStrings = (/"(a)","(b)","(c)","(d)","(e)","(f)","(g)","(h)","(i)","(j)","(k)","(l)" \ ,"(m)","(n)","(o)","(p)"\ ,"(q)","(r)","(s)","(t)"/) ; pnlres@amJust = "BottomRight" pnlres@amJust = "TopLeft" pnlres@gsnPanelFigureStringsFontHeightF = 0.0100 end do ; end of i_file loop ; if (mod(it,144).eq.0.and.it.le.288)then ; if (mod(it,36).eq.0.and.it.le.288)then if (mod(it,54).eq.0.and.it.le.288)then gsn_panel(wks,plots,(/4,3/),pnlres) ; now draw as one plot ii_file = 0 end if if (mod(it,72).eq.0.and.it.gt.288)then gsn_panel(wks,plots_b,(/2,3/),pnlres) ; now draw as one plot ii_file = 0 end if end do ; end of time loop ; gsn_panel(wks,plots,(/4,5/),pnlres) ; now draw as one plot frame (wks) end