; Example script to produce plots for a WRF real-data run, ; with the ARW coordinate dynamics option. 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/NCL_GETools.ncl" begin year = 2013 month = 9 day = 2 dayst = day ; convert single digit months, days, and hours to double digits if(month .lt. 10) then cmon = "0"+month else cmon = month end if if(day .lt. 10) then cday = "0"+day else cday = day end if if(dayst .lt. 10) then cdayst = "0"+dayst else cdayst = dayst end if ; create output file name fout = "DBZ_nexrad_HCHO_"+cmon+"-"+cday type = "ps" wks = gsn_open_wks(type,fout) setvalues NhlGetWorkspaceObjectId() "wsMaximumSize" : 120000000 end setvalues ; Set some Basic Plot Information ; We would like to switch off the Initial Time and Footer info on the plot ; We also do not want titles on the color label bar latmin = 29.8 latmax = 33.5 lonmin = 360.-93.0 lonmax = 360.-88.0 centerlon = (lonmin + lonmax)/2. llonmin = -93.0 llonmax = -86.0 txres = True txres@txFontHeightF = 0.015 txres@txFont = 21 txres@txJust = "CenterLeft" txres@txPerimOn = False res = True res@gsnFrame = False res@gsnDraw = False ; res@gsnAddCyclic = False res@mpMinLatF = latmin ; min(lat2d) res@mpMaxLatF = latmax ; max(lat2d) res@mpMinLonF = lonmin ; min(lon2d) res@mpMaxLonF = lonmax ; max(lon2d) res@mpCenterLonF = centerlon pltres = True pltres@PanelPlot = True pltres@CommonTitle = True pltres1 = True pltres1@FramePlot = False mpres = res mpres@mpFillOn = False mpres@mpGridAndLimbOn = False mpres@mpOutlineOn = True mpres@mpDataBaseVersion = "MediumRes" mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpGridLineColor = "Blue" mpres@mpLimbLineColor = "Green" mpres@mpPerimLineColor = "Black" mpres@mpUSStateLineThicknessF = 2.0 mpres@mpGridLatSpacingF = 1. mpres@mpGridLonSpacingF = 1. mpres@tmXBLabelFontHeightF = 0.015 mpres@tmYLLabelFontHeightF = 0.015 mpres@tmXBMode = "Explicit" mpres@tmXBValues = (/267,268,269,270,271,272/) mpres@tmXBLabels = (/"93W","92W","91W","90W","89W","88W"/) res@tmXBMinorValues = ispan(255,260,1) mpres@tmYLMode = "Explicit" mpres@tmYLValues = (/30,31,32,33/) mpres@tmYLLabels = (/"30N","31N","32N","33N"/) res@tmYLMinorValues = ispan(29,32,1) plots = new ( 2, graphic ) ; Change color map to something that has a grey scale ;gsn_define_colormap(wks,"wxpEnIR") ;gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ;gsn_define_colormap(wks,"radar") gsn_define_colormap(wks,"radar_1") opts = mpres opts@lbTitleOn = True opts@lbTitleString = "" opts@lbTitleFontHeightF = 0.015 opts@cnFillOn = True ; turn on color opts@cnLinesOn = False ; no contour lines opts@lbLabelFontHeightF = 0.015 opts@lbLabelBarOn = True ;False opts@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels opts@cnMinLevelValF = 5 ; set min contour level opts@cnMaxLevelValF = 75 ; 80 ; set max contour level opts@cnLevelSpacingF = 5 ; contour spacing ; opts@gsnSpreadColorStart = 0 ; with colormap "radar" opts@gsnSpreadColorStart = 9 ; with colormap "radar1" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; The NEXRAD input files after decluttered by Cameron Homeyer's IDL routine ifs = 18 ife = 18 ; 26 minute = "00" do ifl=ifs,ife hr = ifl wrfhr = ifl ims = 0 ime = 50 if(hr .eq. 25) then ime = 30 end if do imn = ims,ime,10 minutes = imn if(minutes .lt. 10) then minute = "0"+minutes else minute = minutes end if if(hr .lt. 10) then hour = "0"+hr else hour = hr end if if(hr .gt. 23) then day = dayst + 1 hr = hr-24 if(hr .lt. 10) then delete(hour) hour = "0"+hr else hour = hr end if end if if(month .lt. 10) then cmon = "0"+month else cmon = month end if if(day .lt. 10) then cday = "0"+day else cday = day end if sdate = year+"-"+cmon+"-"+cday+" "+hour+":"+minute+" UTC" s_start = "0000" s_stop = "0000" if(minute .eq. "00") then ; 10 minutes before / after time hrm = hr-1 hrm1 = "00" hrm1 = hrm if(hrm .lt. 10) then hrm1 = "0"+hrm end if s_start = hrm1+":50" s_stop = hour+":10" else minute_10 = minutes-10 minutep10 = minutes+10 s_start = hour+":"+minute_10 s_stop = hour+":"+minutep10 end if ; nexrad_radar_20120623_0120.nc ; 20130902T1735Z.nc dirname = "/glade/p/mmm/barthm/SEAC4RS/Radar/"+year+cmon+"/" filename = year+cmon+cday+"T"+hour+minute+"Z.nc" fname1 = dirname+filename print(" reading from file "+fname1) a = addfile(fname1, "r") time = 0 lon1d = a->Longitude lat1d = a->Latitude alt1d = a->Altitude refl = a->Reflectivity ; (nz,ny,nx) refl@_FillValue = -9999 if (any(isnan_ieee(refl))) then value = refl@_FillValue replace_ieeenan (refl, value, 0) end if print(" ") print(" refl fill value "+refl@_FillValue) diml = dimsizes(refl) print(" dimsizes refl "+diml+" max refl "+max(refl) ) print(" ") nz = diml(0) ny = diml(1) nx = diml(2) lon2d = new( (/ny,nx/), float) lat2d = new( (/ny,nx/), float) do j=0,ny-1 lon2d(j,:) = lon1d end do do i=0,nx-1 lat2d(:,i) = lat1d end do print(" min/max longitude "+min(lon2d(0,:))+" "+max(lon2d(0,:)) ) print(" min/max latitude "+min(lat2d(0,:))+" "+max(lat2d(0,:)) ) lat2d@units = "degrees_N" lon2d@units = "degrees_E" refl@lon2d = lon2d refl@lat2d = lat2d print(" max refl "+max(refl)) mrefl = dim_max_n( refl, 0) mrefl@description = "d) Maximum Reflectivity at "+sdate mrefl@units = "" ; "(dBZ)" print(" max mrefl "+max(mrefl)) mrefl@lon2d = lon2d mrefl@lat2d = lat2d gsn_define_colormap(wks,"radar_1") plot = gsn_csm_contour_map_ce(wks, mrefl, opts) draw(plot) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Change color map for aircraft data gsn_define_colormap(wks,"BlAqGrYeOrReVi200") colors = new( (/10/), integer) do imkr = 0,9 colors(imkr) = 11 + imkr*20 ; 11, 31, 51, 71, 91, 111, 131, 151, 171, 191 end do colors(0) = 0 colors(1) = 25 ; put red colors in middle of colormap, blue colors at end, green at beginning ; colors(0:7) = colors(2:9) ;; colors(5) = 11 + 4*20 ;; colors(6) = 11 + 3*20 ;; colors(7) = 11 + 2*20 ; colors(8) = 25 ; colors(9) = 11 + 1*20 ; print(" colors "+colors) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Get DC-8 Aircraft Data ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; vars = (/ "co", "o3", "h2o2", "hcho", "hno3", "ch3ooh", "no" /) spcs = (/ "CO", "O3", "H2O2", "HCHO", "HNO3", "CH2O", "NO" /) idc8 = (/ 62, 70, 350, 65, 354, 66, 69 /) ; NO2 = 71 TDLIF, 67 ESRL dc8scl=(/ 1., 1., .001, .001, .001, .001, 1. /) ; scale to ppbv from data file units ilvls =(/ 1, 1, 0, 0, 0, 0, 0 /) ; contour levels selection, see code below ivar = 3 csp = spcs(ivar) var = vars(ivar) print(" variable: "+var+" "+csp+" "+idc8(ivar)+" "+dc8scl(ivar)+" "+ilvls(ivar) ) ;;;;;;;;;;; ; DC8 Data ;;;;;;;;;;; obsdirname = "/glade/p/mmm/barthm/SEAC4RS/Aircraft/" fnamedc8 = "SEAC4RS-mrg01-dc8_merge_2013"+cmon+cdayst+"_R6.ict" print("DC8 file: "+obsdirname+fnamedc8) ictVar_DC8 = read_ict(obsdirname+fnamedc8) ; print(" DC8 varname: "+ictVar_DC8@VariableNames+" "+ictVar_DC8@VariableUnits ) iskip = 0 if(iskip .eq. 0) then dc8_time_s = doubletofloat(ictVar_DC8(:,0)) ; Time, s dc8_galt_s = doubletofloat(ictVar_DC8(:,53)) ; GPS_ALT_MMS dc8_lat_s = doubletofloat(ictVar_DC8(:,51)) ; GPS_LAT_MMS dc8_lon_s = doubletofloat(ictVar_DC8(:,52)) ; GPS_LON_MMS dc8_u_s = doubletofloat(ictVar_DC8(:,46)) ; U_MMS dc8_v_s = doubletofloat(ictVar_DC8(:,47)) ; V_MMS iv = idc8(ivar) dc8_co_s = doubletofloat(ictVar_DC8(:,iv)) if(var .eq. "no") then dc8_no2_s = doubletofloat(ictVar_DC8(:,67)) ;; dc8_no2_s = doubletofloat(ictVar_DC8(:,71)) ; TDLIF data ;; dc8_no2_s = dc8_no2_s * 0.001 ; ppbv dc8_co_s = dc8_co_s + dc8_no2_s ; ppbv end if dc8_var_s = dc8_co_s *dc8scl(ivar) dc8_var_s@_FillValue = -999. dc8_galt_s@_FillValue = -999. dc8_lon_s@_FillValue = -999. dc8_lat_s@_FillValue = -999. ;;;; limit data to an altitude range: dc8_var_s = where(dc8_galt_s .ge. hgt-1, dc8_var_s, dc8_var_s@_FillValue) print(" max DC8 "+var+" "+max(dc8_var_s) ) avgtime = 30. iavgt = 30 dimda = dimsizes(dc8_galt_s) dimdc = dimsizes(dc8_var_s) ndc8_min = floattointeger(dimdc(0)/avgtime) dc8_time = new( (/ndc8_min/), float) dc8_galt = new( (/ndc8_min/), float) dc8_var = new( (/ndc8_min/), float) dc8_lat = new( (/ndc8_min/), float) dc8_lon = new( (/ndc8_min/), float) dc8_u = new( (/ndc8_min/), float) dc8_v = new( (/ndc8_min/), float) dc8_galt@_FillValue = -999. dc8_var@_FillValue = -999. dc8_lat@_FillValue = -999. dc8_lon@_FillValue = -999. dc8_u@_FillValue = -999. dc8_v@_FillValue = -999. do im = 0,ndc8_min-1 i1 = im*iavgt i2 = i1+iavgt-1 i3 = i1+iavgt/2 if(all(ismissing(dc8_var_s(i1:i2)))) then dc8_time(im) = avg(dc8_time_s(i1:i2)) dc8_galt(im) = dc8_galt_s@_FillValue dc8_var(im) = dc8_var_s@_FillValue dc8_lat(im) = dc8_lat_s@_FillValue dc8_lon(im) = dc8_lon_s@_FillValue dc8_u(im) = dc8_u_s@_FillValue dc8_v(im) = dc8_v_s@_FillValue else dc8_time(im) = avg(dc8_time_s(i1:i2)) dc8_galt(im) = avg(dc8_galt_s(i1:i2)) dc8_var(im) = avg(dc8_var_s(i1:i2)) dc8_lat(im) = avg(dc8_lat_s(i1:i2)) dc8_lon(im) = avg(dc8_lon_s(i1:i2)) dc8_u(im) = avg(dc8_u_s(i1:i2)) dc8_v(im) = avg(dc8_v_s(i1:i2)) ; print(" "+im+" "+dc8_time_s(i3)+" "+dc8_galt(im)+" "+dc8_var(im) +" "+dc8_lat(im) +" "+dc8_lon(im) ) end if end do ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Plot var data gsres = True ; find start, stop times for hour of interest wrfsec = wrfhr*3600. + minutes*60. startsec = wrfsec - 600. stopsec = wrfsec + 600. ; startsec = 25.*3600. + 16.*60. ; from outflow times ; stopsec = 25.*3600. + 19.*60. ; from outflow times c = ind(dc8_time .ge. startsec) im1 = c(0) delete(c) c = ind(dc8_time .ge. stopsec) im2 = c(0) delete(c) if(ismissing(im2)) then im2 = ndc8_min-1 end if do im = im1,im2 ; 0,ndc8_min-1 if(ismissing(dc8_var(im))) then else gsres@gsMarkerIndex = 4 ; 8=upward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 4. gsres@gsMarkerColor = "black" gsn_polymarker(wks,plot,dc8_lon(im),dc8_lat(im),gsres) gsres@gsMarkerIndex = 16 ; 7=downward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 2. if(ilvls(ivar) .eq. 0) then contourlvls = (/ .05, .07, .10, .12, .15, .17, .2, .25, .3 /) ; contourlvls = (/ .05, .07, .1, .2, .3, .5, .7, 1., 2. /) contourlvls = (/ .10, .2, .50, .7, 1., 2., 3., 4., 5. /) ; for other species, but adjusted for HCHO else contourlvls = (/ 50., 70., 80., 90., 100., 110., 120., 130., 150./) ; contour levels for CO and O3 end if gsres@gsMarkerColor = colors(0) if(dc8_var(im) .ge. contourlvls(0)) then gsres@gsMarkerColor = colors(1) end if if(dc8_var(im) .ge. contourlvls(1)) then gsres@gsMarkerColor = colors(2) end if if(dc8_var(im) .ge. contourlvls(2)) then gsres@gsMarkerColor = colors(3) end if if(dc8_var(im) .ge. contourlvls(3)) then gsres@gsMarkerColor = colors(4) end if if(dc8_var(im) .ge. contourlvls(4)) then gsres@gsMarkerColor = colors(5) end if if(dc8_var(im) .ge. contourlvls(5)) then gsres@gsMarkerColor = colors(6) end if if(dc8_var(im) .ge. contourlvls(6)) then gsres@gsMarkerColor = colors(7) end if if(dc8_var(im) .ge. contourlvls(7)) then gsres@gsMarkerColor = colors(8) end if if(dc8_var(im) .ge. contourlvls(8)) then gsres@gsMarkerColor = colors(9) end if gsn_polymarker(wks,plot,dc8_lon(im),dc8_lat(im),gsres) ; print(" var "+im+" "+dc8_time(im)+" "+dc8_galt(im)+" "+dc8_lat(im) +" "+dc8_lon(im)+" "+dc8_var(im) +" " +gsres@gsMarkerColor ) end if end do txres@txFontColor = "black" xpt = 0.80 ypt = 0.09 ; 0.148 gsn_text_ndc(wks, "dBZ", xpt, ypt, txres) ; xpt = (dc8_lon(im1) - llonmin)/(llonmax-llonmin) - 0.08 ; ypt = (dc8_lat(im1) - latmin)/(latmax-latmin) - 0.04 ; print(" im1 lon/lat "+dc8_lon(im1)+" "+dc8_lat(im1)+" "+xpt+" "+ypt+" s_start="+s_start ) ; gsn_text_ndc(wks, s_start, xpt, ypt, txres) ; xpt = (dc8_lon(im2) - llonmin)/(llonmax-llonmin) - 0.05 ; ypt = (dc8_lat(im2) - latmin)/(latmax-latmin) - 0.181 ; print(" im2 lon/lat "+dc8_lon(im2)+" "+dc8_lat(im2)+" "+xpt+" "+ypt+" s_stop="+s_stop ) ; gsn_text_ndc(wks, s_stop, xpt, ypt, txres) clvls = 1000.*contourlvls ltxt = flt2string(clvls) do i = 0,9 xpt = 360.-88.85 ypt = 30.05 + 0.235*(i) gsres@gsMarkerIndex = 4 ; 8=upward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 4. gsres@gsMarkerColor = "black" gsn_polymarker(wks,plot,xpt,ypt,gsres) gsres@gsMarkerIndex = 16 ; 7=downward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 2. gsres@gsMarkerColor = colors(i) gsn_polymarker(wks,plot,xpt,ypt,gsres) end do xpt = 0.75 do i = 0,8 ypt = 0.033*i + 0.29 print(" i lon/lat "+dc8_lon(im2)+" "+dc8_lat(im2)+" "+xpt+" "+ypt+" "+ltxt(i) ) gsn_text_ndc(wks, ltxt(i), xpt, ypt, txres) if(i.eq.8) then gsn_text_ndc(wks, "pptv", xpt+0.05, ypt-0.003, txres) end if end do end if ; iskip ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Close the Plot frame(wks) end do ; imn end do ; ifl loop ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end