load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 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" begin ; barry.h.lynn@gmail.com ; reads WRF ensemble data and plots the wind chill ; note, there is some extraneous code commented out. ; note, you should have just one output file with no sub-directories in each ; wrf output ensemble directory ; then, n loop below should increment by 1 or how many files necessary to reach ; next directory. ; or you need to specify the increment from directory to directory ; ndir_WRF had 180 files, so I looped every 20, other wise ; I read in a time file with dates I specifiy. ; You can read change this to read in a time file anyway you want. ; Or, you can use the times from the WRF files themselves. ; see: http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/Examples/SPECIAL/wrf_panel1.ncl ; pnlres@txString = t@description + " (" + t@units + ")" ; but I use time_var, instead. ; define number of panel plots plot = new(4,graphic) asc_time_file = "wrf_times.file" num_col = 5 num_times = numAsciiRow(asc_time_file) year = new((/num_times/),integer) month = new((/num_times/),integer) day = new((/num_times/),integer) hour = new((/num_times/),integer) minute = new((/num_times/),integer) dims = dimsizes(year) print(dims) ; z1 = asciiread("jan26-28_2015.txt",(/num_obs,num_col/),"float") z1 = asciiread(asc_time_file,(/num_times,num_col/),"integer") do it = 0,num_times-1 year = z1(it,0) end do year = z1(:,0) month = z1(:,1) day = z1(:,2) hour = z1(:,3) minute = z1(:,4) cities = (/" Tel-Aviv"," Jerusalem"," Haifa"," Hermon"," Arad"," Eilat"," Med. Sea"," Amman"," Damascus"," Sefad",\ " Efrat"," Mitzpe Ramon", " Hebron", " Ashquelon"," Ein Gedi"," Netanya"," Beer-7"," El Arish"," Akko", \ " Nahariya"," Tyre"," Zikhron"," Ashdod"/) lat = (/ 32.0670, 31.78300,32.808,33.368,31.2475,29.55,33.75,31.9454,33.4,32.967, \ 31.65,30.617,31.53,31.677,31.450,32.334,31.25, 31.13,32.92,33.006,33.27,32.57, \ 31.80/) lon = (/ 34.767,35.223,35.057,35.783,35.203515,34.95,34.75,35.9284,36.5,35.493, \ 35.133,34.80,35.10,34.583,35.383,34.858,34.80,33.80,35.10,35.095,35.20,34.943, \ 34.655/) ncities = dimsizes(cities) markid = new((/4,ncities/),graphic) textid = new((/4,ncities/),graphic) txres = True ; Set up text resources ; d03 txres@txFontHeightF = 0.008 ; d02 ; txres@txFontHeightF = 0.005 txres@txJust = "CenterLeft" mkres = True ; Set up marker resources mkres@gsMarkerIndex = 16 ; Filled dot mkres@gsMarkerSizeF = .0015 ; Size of dot ; not in use:set bounding box ; bounding box for January 23 2016 ; lats = (/ 39.00, 43.50 /) ; lons = (/ -77.0, -69.50 /) ; bounding box for January 26 2015 ; lats = (/ 37.75, 44.00 /) ; lons = (/ -77.0, -68.00 /) ; bounding box for February 08 2013 ; lats = (/ 38.00, 42.50 /) ; lons = (/ -77.0, -68.75 /) ; set it_beg, it_inc ; 4 km it_beg = 6 it_old = it_beg-1 it_inc = 1 ; my resolution; not in use ; resol = "1.3 km " ; resol = "4.0 km " resol = "4.0 km " ; my_date = "02_08_2013" my_date = "01_22_2016" ; my_date = "01_23_2016" ; my_date = "" ; outfile="t2c_ave_4.0km_" + my_date ; outfile="t2c_ave_1.3km_" + my_date outfile="WRF_temp_4.0km_ens" wks = gsn_open_wks("pdf",outfile) DATADir = "/home/cust100021_vol1/barry/cust100021_vol2/GEFS_4KM" print("DATADir = " + DATADir) ; for when there is no date ; find directories with this date, etc cmd = "find " + DATADir + "/GEFS* -type d -print" dirWRF = systemfunc (cmd) print(dirWRF) ; this is the total number of directories within the listed directory header ; it can be just the number of, for example, "WRF_XX" directories or also all subdirectories within n_dirWRF = dimsizes(dirWRF) ; print("n_dirWR" + n_dirWRF) ; this is the total number of directories within the listed directory header ; it can be just the number of, for example, "WRF_XX" directories or also all subdirectories within n_dirWRF = dimsizes(dirWRF) ; ; this is the name of the file ;;yr = year(0,1,1) ;mt= month(0,1,1) ; dy= day(0,1,1) ; hl= hour(0,1,1) ; wrfname = "wrfout_d02_20" + yr + "-" + mt + "-" dy + "_" + hr + ":00:00" all_files = systemfunc("ls " + dirWRF(0) + "/wrfout_d02*") print(all_files) ; not in use ; wrfname = "wrfout_d03_2016-01-22_00:00:00" ; wrfname = "wrfout_d03_2016-01-23_00:00:00_orig" ; wrfname = "wrfout_d04_2016-01-23_00:00:00" ; wrfname = "wrfout_d02_2016-01-22_00:00:00" ; wrfname = "wrfout_d03_2015-01-26_00:00:00_3dom" ; wrfname = "wrfout_d04_2015-01-26_17:30:00" ; wrfname = "wrfout_d04_2013-02-08_00:00:00" ; wrfname = "wrfout_d04_2013-02-08_00:00:00" ntimes = 30 ; n_loops = 10 n_loops = 9 print("dirWRF(0) = " + dirWRF(0)) filename = all_files(0) print("filename = " + filename) a = addfile(filename,"r") times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) ; number of times in the file ; print ("ntimes = " + ntimes) rainnc_old = wrf_user_getvar(a,"RAINNC",it_old) rain_tot = rainnc_old rain_tot = 0 dims2d=dimsizes(rainnc_old) prob_100=new((/dims2d(0),dims2d(1)/),float) prob_50=new((/dims2d(0),dims2d(1)/),float) prob_25=new((/dims2d(0),dims2d(1)/),float) ; num_times should not be more than ntimes it_val = it_beg ; do it = it_beg,ntimes-1,it_inc do it = it_beg,num_times-1,it_inc ; do it = it_beg,it_beg,it_inc print("it = " + it) print("it_old = " + it_old) i_loop = 0 ; increment by number of sub directories ; do n=0,n_dirWRF -1,20 do n=0,180,20 ; do n=0,0,1 print("i_loop " + i_loop) print("dirWRF(n) = " + dirWRF(n)) all_files = systemfunc("ls " + dirWRF(n) + "/wrfout_d02*") filename = all_files(0) print("filename = " + filename) a = addfile(filename,"r") times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) ; number of times in the file print ("ntimes = " + ntimes) rainnc_old = wrf_user_getvar(a,"RAINNC",it_old) rainc_old = wrf_user_getvar(a,"RAINC",it_old) rain_old = rainnc_old+rainc_old print("it = " + it) print("it_old = " + it_old) rainnc = wrf_user_getvar(a,"RAINNC",it) rainc = wrf_user_getvar(a,"RAINC",it) tc2 = wrf_user_getvar(a,"AFWA_WCHILL",it) ; T2 in Kelvin tc2 = tc2 - 273.16 printMinMax(tc2,False) if (i_loop.eq.0)then t2c_ave=tc2/(n_loops+1) end if if (i_loop.gt.0)then t2c_ave=t2c_ave+tc2/(n_loops+1) end if dims2d=dimsizes(rainnc) if (i_loop.eq.0)then rain_tot=rain_tot do j=1,dims2d(0)-1 do i=1,dims2d(1)-1 if (tc2(j,i).lt.10)then prob_100(j,i) = 100./n_loops end if if (tc2(j,i).gt.10)then prob_100(j,i) = 0. end if if (tc2(j,i).le.5)then prob_50(j,i) = 100./n_loops end if if (tc2(j,i).gt.5)then prob_50(j,i) = 0. end if if (tc2(j,i).le.00)then prob_25(j,i) = 100./n_loops end if if (tc2(j,i).gt.0)then prob_25(j,i) = 0. end if end do end do end if if (i_loop.gt.0)then do j=1,dims2d(0)-1 do i=1,dims2d(1)-1 if (tc2(j,i).le.10)then prob_100(j,i) = prob_100(j,i) + 100./n_loops end if if (tc2(j,i).le.5)then prob_50(j,i) = prob_50(j,i) + 100./n_loops end if if (tc2(j,i).le.0)then prob_25(j,i) = prob_25(j,i) + 100./n_loops end if end do end do end if max_100 = max(prob_100) max_50 = max(prob_50) max_25 = max(prob_25) do j=1,dims2d(0)-1 do i=1,dims2d(1)-1 if (max_100 .eq. 0)then prob_100(j,i) = -1*i end if if (max_50 .eq. 0)then prob_50(j,i) = -1*i end if if (max_25 .eq. 0)then prob_25(j,i) = -1*i end if end do end do i_loop = i_loop + 1 if (i_loop-1.eq.n_loops)then ; t2c_ave = t2c_ave/25.54 xlat2d=wrf_user_getvar(a,"XLAT",it) xlon2d=wrf_user_getvar(a,"XLONG",it) prob_100@lat2d = xlat2d ; this is for plotting prob_100@lon2d = xlon2d ; using gsn_csm routines prob_50@lat2d = xlat2d ; this is for plotting prob_50@lon2d = xlon2d ; using gsn_csm routines prob_25@lat2d = xlat2d ; this is for plotting prob_25@lon2d = xlon2d ; using gsn_csm routines t2c_ave@lat2d = xlat2d ; this is for plotting t2c_ave@lon2d = xlon2d ; using gsn_csm routines rain_tot@lat2d = xlat2d ; this is for plotting rain_tot@lon2d = xlon2d ; using gsn_csm routines ; 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 ; loc = wrf_user_ll_to_ij(a, lons, lats, True) ; x_start = loc(0,0) - 1 ; x_end = loc(0,1) - 1 ; y_start = loc(1,0) - 1 ; y_end = loc(1,1) - 1 ;;;; ;;; create plot (pdf) ;;;; ; copy coordinate variables ; copy_VarCoords(rainnc,t2c_ave) ; copy_VarCoords(rainnc,rain_tot) gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") ; load color table ;;; Font Style res1 = True ; plots modification on res1@txFont = "Helvetica" res1@txFontQuality = "High" ;;; Title ; first calculate time 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) it_val = it_val +1 print("new_time_array = " + new_time_array) ; printVarSummary(new_time_array) format = "" ;; defaults to "%H%M EST %d %c %Y" format = "%H%M IST %d %c %Y" ; 2.16667 = + 2 hours ;+ 10 minutes (for first time on map in EST) ; new_time_array = new_time_array + 2.16667 ; new_time_array = new_time_array -5 new_time_array = new_time_array + 2 time_var=cd_string(new_time_array,format) print("time_var = " + time_var) res1@gsnDraw = False ; do not draw the plot res1@gsnFrame = False ; do not advance the frame ; res1@tiMainString = time_var res1@gsnLeftString = "WRF " + resol + " Mean Wind Chill" res1@gsnRightString = "mm" res1@gsnRightString = "[~S~o~N~C]" ; title string ;;; Map res1 = True res1@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res1@mpDataSetName = "Earth..4" ; high resolution res1@mpDataBaseVersion = "HighRes" ; res1@mpOutlineBoundarySets = "National" ; National borders ; res1@mpOutlineBoundarySets = "AllBoundaries" ; National borders ; res1@mpOutlineBoundarySets = "USStates" ; National borders res1@mpGeophysicalLineColor = "Black" res1@mpNationalLineColor = "Black" res1@mpUSStateLineColor = "Black" res1@mpNationalLineColor = "Black" res1@mpGridLineColor = 0 res1@mpLimbLineColor = "Black" res1@mpPerimLineColor = "Black" res1@mpUSStateLineThicknessF = 2.0 res1@mpGeophysicalLineThicknessF = 1.5 ;;; Temperature ave_t = avg(t2c_ave) res1@cnFillOn = True ; turn on color fill res1@cnLinesOn = False ; turn of contour lines ;;; color fill tresholds and colors res1@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res1@cnLevels = (/ -26.,-24.,-22.,-20.,-18.,-16.,-14.,\ -12.,-10.,-8.,-6.,-4.,-2.,0.,\ 2.,4.,6.,8.,10.,12.,14.,16.,18.,\ 20./) res1@cnFillColors = (/ 19,22,25,31,35,39,43,48,53,\ 61,64,67,71,89,95,99,104,108,111,115,123,\ 130,135,143,148/) ; set the colors to be used ;;; label bars style ;;; if not set standard values are used. ;;; so you normally can delete these lines without any problems res1@cnInfoLabelOn = False ; turn off contour info label res1@lbAutoManage = False ; control label bar res1@pmLabelBarDisplayMode = "Always" ; turns on label bar res1@lbOrientation = "Vertical" res1@pmLabelBarSide = "Right" res1@lbLabelAutoStride = True res1@pmLabelBarWidthF = 0.20 res1@pmLabelBarHeightF = 0.7 res1@lbLabelFontHeightF = .020 res1@lbPerimOn = False res1@lbLabelFont = "Helvetica" ; label font res1@lbTitleOn = True ; turn on title res1@lbTitleString = "Wind Chill [~S~o~N~C]" ; title string res1@lbTitlePosition = "Bottom" ; title position res1@lbTitleFontHeightF= .030 ; make title smaller res1@lbTitleDirection = "Across" ; title direction ;;; Map Projection (min and max Lon and Lat) res1@gsnAddCyclic = True ; data already has cyclic point ; this must also be set for any zoom min_lat = 28.5 max_lat = 34.5 min_lon = 33.0 max_lon = 38.0 ; res1@mpProjection = "LambertConformal" res1@mpMinLatF = min_lat res1@mpMaxLatF = max_lat res1@mpMinLonF = min_lon res1@mpMaxLonF = max_lon res1@mpLimitMode="LatLon" ; to have a common label bar, both plots should be set to the same interval ; b/c the label bar is drawn by default from the interval of the first plot. ; res@cnLevelSelectionMode = "ManualLevels" ; res@cnMinLevelValF = -10. ; res@cnMaxLevelValF = 45. ; res@cnLevelSpacingF = 5. res2 = True res2@gsnDraw = False ; do not draw the plot res2@gsnFrame = False ; do not advance the frame res2@gsnLeftString = "WRF " + resol + "Prob < 0 [~S~o~N~C]" res2@gsnRightString = "%" res2@cnFillOn = True ; turn on color fill res2@cnLinesOn = False ; turn of contour lines ;;; Map res2@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res2@mpDataSetName = "Earth..4" ; high resolution res2@mpDataBaseVersion = "HighRes" ; res1@mpOutlineBoundarySets = "National" ; National borders ; res1@mpOutlineBoundarySets = "AllBoundaries" ; National borders ; res2@mpOutlineBoundarySets = "USStates" ; National borders res2@mpNationalLineColor = "Black" ; National borders color res2@mpGeophysicalLineColor = "Black" res2@mpNationalLineColor = "Black" res2@mpUSStateLineColor = "Black" res2@mpNationalLineColor = "Black" res2@mpGridLineColor = 0 res2@mpLimbLineColor = "Black" res2@mpPerimLineColor = "Black" res2@mpUSStateLineThicknessF = 3.0 res2@mpGeophysicalLineThicknessF = 1.5 ;;; Prec res2@cnFillOn = True ; turn on color fill res2@cnLinesOn = False ; turn of contour lines ;;; color fill tresholds and colors res2@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res2@cnLevels = (/ 1,10., \ 20, 30,40,50,60,80,90/) res2@cnFillColors = (/"White","AntiqueWhite","AntiqueWhite3", \ "chartreuse", \ "chartreuse3","ForestGreen", \ "Yellow","Orange","Red","HotPink3","HotPink1","Violet"/) ;;; label bars style ;;; if not set standard values are used. ;;; so you normally can delete these lines without any problems res2@cnInfoLabelOn = False ; turn off contour info label res2@lbAutoManage = False ; control label bar res2@pmLabelBarDisplayMode = "Always" ; turns on label bar res2@lbOrientation = "Vertical" res2@pmLabelBarSide = "Right" res2@lbLabelAutoStride = True res2@pmLabelBarWidthF = 0.20 res2@pmLabelBarHeightF = 0.7 res2@lbLabelFontHeightF = .020 res2@lbPerimOn = False res2@lbLabelFont = "Helvetica" ; label font res2@lbTitleOn = True ; turn on title res2@lbTitleString = "(%)" ; title string res2@lbTitlePosition = "Bottom" ; title position res2@lbTitleFontHeightF= .030 ; make title smaller res2@lbTitleDirection = "Across" ; title direction ;;; Map Projection (min and max Lon and Lat) res2@gsnAddCyclic = True ; data already has cyclic point ; this must also be set for any zoom ; res2@mpProjection = "LambertConformal" res2@mpMinLatF = min_lat res2@mpMaxLatF = max_lat res2@mpMinLonF = min_lon res2@mpMaxLonF = max_lon res2@mpLimitMode="LatLon" ; res3 = res2 ; res3@gsnLeftString = "GEFS Prob > 0.50 in" res3@gsnLeftString = "WRF " + resol + "Prob < 5 [~S~o~N~C]" res3@lbTitleString = "(%)" ; title string res4 = res2 ; res4@gsnLeftString = "GEFS Prob > 1.00 in" res4@gsnLeftString = "WRF " + resol + "Prob < 10 [~S~o~N~C]" res4@lbTitleString = "(%)" ; title string ; plot(0) = gsn_csm_contour_map(wks,u,res) ; plot(1) = gsn_csm_contour_map(wks,v,res) ;************************************************ ; create panel ;************************************************ resP = True ; modify the panel plot ; resP@gsnPanelMainString = "A plot with a common label bar" ; resP@gsnPanelMainString = time_var ; use this for NCL V6.3.0 and earlier resP@txString = time_var ; resP@gsnPanelLabelBar = True ; add common colorbar ; resP@lbLabelFontHeightF = 0.007 ; make labels smaller plot(0) = gsn_csm_contour_map(wks,t2c_ave,res1) ; create plot plot(1) = gsn_csm_contour_map(wks,prob_25,res2) ; create plot plot(2) = gsn_csm_contour_map(wks,prob_50,res3) ; create plot plot(3) = gsn_csm_contour_map(wks,prob_100,res4) ; create plot ; plot = gsn_csm_contour_map(wks,rain_tot,res1) ; create plot ;For adding markers and text, you don't need to do this in the do loop. You can simply call both procedures once with the array of lat/lon values: do p = 0,3 markid(p,:) = gsn_add_polymarker(wks,plot(p),lon,lat,mkres) textid(p,:) = gsn_add_text(wks,plot(p),cities,lon,lat,txres) end do gsn_panel(wks,plot,(/2,2/),resP) ; now draw as one plot ;Finally, in order to see the markers and the cities, you actually have to draw the map and advance the frame yourself. The gsn_add_xxxx calls just attach the stuff to the map, but nothing gets drawn until you call: ; res@mpOutlineOn = False ; draw(plot) ; frame(wks) end if end do ; end of directory loop it_old = it end do ; end of file read end