; ; Script to produce standard plots for a WRF real-data run, ; with the MASS 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/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" undef("get_wrfoutfname") function get_wrfoutfname(domain:string, yr:integer, month:integer, day:integer, hr:integer, minutes:integer) begin 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 do while(hr .gt. 23) day = day + 1 if(day .lt. 10) then cday = "0"+day else cday = day end if hr = hr-24 print(" minutes "+minutes+" hr "+hr+" day "+day) if(hr .lt. 10) then delete(hour) hour = "0"+hr else hour = hr end if end do if(yr .lt. 10) then year = "000"+yr else if(yr .lt. 100) then year = "00"+yr else if(yr .lt. 1000) then year = "0"+yr else year = yr end if end if end if if(month .lt. 10) then mon = "0"+month else mon = month end if if(day .lt. 10) then cday = "0"+day else cday = day end if filename = "wrfout_"+domain+"_"+year+"-"+mon+"-"+cday+"_"+hour+":"+minute+":00.nc" print(" reading from file "+filename) delete(year) delete(mon) delete(cday) delete(hour) delete(minute) return(filename) end ;------------------------------------------------------------------------------------------- ;****************************************************************************************************; begin minutes = 0 if(minutes .lt. 10) then cmin = "0"+minutes else cmin = minutes end if if(hr .lt. 10) then chr = "0"+hr else chr = hr end if if(day .lt. 10) then cday = "0"+day else cday = day end if if(month .lt. 10) then mon = "0"+month else mon = month end if cdate = year+"_"+mon+"-"+cday+"_"+chr+cmin dirname1 = "./wrfout_linked/" fout = "Frappe_XSEC_loc_"+cdate wks = gsn_open_wks("png" ,fout) setvalues NhlGetWorkspaceObjectId() "wsMaximumSize" : 90000000 end setvalues ;******************************************************** ; Set up options for plotting ;******************************************************** res = True ;res@InitTime = False ;res@Footer = False res@NoHeaderFooter = True res@vpWidthF = 0.8 ;0.22 res@vpHeightF = 0.6 ;0.22 res@vpXF = 0.05 ; 0.1, 0.37, 0.64 res@vpYF = 0.9 res@UnitLabel = "ppbv" res@cnFillMode = "CellFill" res@cnLinesOn = False ; no contour lines res@cnInfoLabelOn = False res@cnLineLabelsOn = False res@cnLevels = (/ 10., 15., 20., 25., 30., 35., 40., 45., 50., 55. /) res@cnFillColors = (/"white","dodgerblue","cyan","forestgreen","green","greenyellow","yellow","gold","orange","orangered","red" /) pltres = True pltres1 = True pltres1@FramePlot = False mpres = True mpres1 = True mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpCountyLineColor = "Black" mpres@mpUSStateLineThicknessF = 3 mpres@mpGridLineColor = "Gray" mpres@mpLimbLineColor = "Gray" mpres@mpPerimLineColor = "Black" mpres@mpDataBaseVersion = "Ncarg4_1" mpres@mpDataSetName = "Earth..2" mpres@mpOutlineBoundarySets = "AllBoundaries" ;******************************************************** ; Get data ;******************************************************** ; dirname1 = "/glade/p/work/stacy/WRF_F/" ;dirname1 = "/glade/scratch/mizzi/FRAPPE_FORECASTING/realtime/gfs/2014061612/wrf_rundir/" fname1 = get_wrfoutfname("d02", year, month, day, hr, minutes) print(" reading from file "+dirname1+fname1) a1 = addfile (dirname1 + fname1, "r") times = wrf_user_list_times(a1) ; get times in the file ntimes = dimsizes(times) ; number of times in the file it = 0 time = it vegfrac = wrf_user_getvar(a1, "VEGFRA", time) lat2d = wrf_user_getvar(a1, "XLAT", time) lon2d = wrf_user_getvar(a1, "XLONG", time) z = wrf_user_getvar(a1,"z",time) ht = wrf_user_getvar(a1,"HGT",time) dimh = dimsizes(ht) ny = dimh(0) nx = dimh(1) mgrpl = wrf_user_getvar(a1,"GRPL_MAX",time) print(" max mgrpl "+max(mgrpl)) ;************************************************ ; create the plot ;************************************************ ; gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") gsn_define_colormap(wks,"WhiteBlueGreenYellowRed") ; 254 colors ; gsn_define_colormap(wks,"topo_15lev") ; i = NhlNewColor(wks,0.8,0.8,0.8) ; add gray to color map ; res@gsnSpreadColorStart = 5 opts = True opts@InitTime = False opts@Footer = False opts@vpWidthF = 0.8 opts@vpHeightF = 0.6 opts@vpXF = 0.05 opts@vpYF = 0.9 opts@UnitLabel = "meters" opts@cnFillMode = "RasterFill" opts@cnLinesOn = False opts@cnInfoLabelOn = False opts@cnLineLabelsOn = False opts@cnFillOn = True ; turn on color opts@lbLabelBarOn = True ; True opts@lbTitleOn = False opts@lbTitleString = "" opts@lbTitleFontHeightF = 0.015 opts@lbLabelFontHeightF = .015 opts@lbLabelStride = 2 opts@MainTitle = "" ; times(time) opts@cnLevelSelectionMode = "AutomaticLevels" ; set manual contour levels ; opts@cnLevels = fspan(0,3000.,13) ; Create tick marks ; opts@cnLevels = fspan(0,70.,15) ; Create tick marks ; opts@cnLevels = (/ 0.1, 1., 2., 5., 7., 10., 15., 20., 25., 30./) ; opts@cnFillColors = (/"white","dodgerblue","cyan","springgreen","green","greenyellow","yellow","gold","orange","orangered","red" /) ; opts@cnLevels = (/ 0.001 /) print("ter.hgt max "+max(ht)) print("grpl.max max "+max(mgrpl)) contour0 = wrf_contour(a1,wks,ht,opts) resb = True ; polyline mods desired resb@gsLineThicknessF = 5.0 ; thickness of lines resb@gsLineColor = "deeppink3" resb@gsLineDashPattern = 0 three_lons = (/-109., -105.7, -105.2, -104.7/) three_lats = (/ 38.5, 39.5, 40., 40.5/) sitex = three_lons sitey = three_lats nsites = dimsizes(three_lons) ;; do iloc = 0,dimsizes(three_lats)-1 nm = getind_latlon2d (lat2d,lon2d, three_lats(0), three_lons(0)) n = nm(0,0) m = nm(0,1) plane1 = (/ m, n /) nm = getind_latlon2d (lat2d,lon2d, three_lats(1), three_lons(1)) n = nm(0,0) m = nm(0,1) plane2 = (/ m, n /) nm = getind_latlon2d (lat2d,lon2d, three_lats(2), three_lons(2)) n = nm(0,0) m = nm(0,1) plane3 = (/ m, n /) nm = getind_latlon2d (lat2d,lon2d, three_lats(3), three_lons(3)) n = nm(0,0) m = nm(0,1) plane4 = (/ m, n /) angles = (/ 0., 90. /) optx = False ; start and end points not specified ixsec=0 lon1_vplane = wrf_user_intrp2d(lon2d,plane1,angles(ixsec),optx) lat1_vplane = wrf_user_intrp2d(lat2d,plane1,angles(ixsec),optx) lon2_vplane = wrf_user_intrp2d(lon2d,plane2,angles(ixsec),optx) lat2_vplane = wrf_user_intrp2d(lat2d,plane2,angles(ixsec),optx) lon3_vplane = wrf_user_intrp2d(lon2d,plane3,angles(ixsec),optx) lat3_vplane = wrf_user_intrp2d(lat2d,plane3,angles(ixsec),optx) lon4_vplane = wrf_user_intrp2d(lon2d,plane4,angles(ixsec),optx) lat4_vplane = wrf_user_intrp2d(lat2d,plane4,angles(ixsec),optx) diml = dimsizes(lat1_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat1_vplane(0), lon1_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat1_vplane(nl), lon1_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc101 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon1_vplane, lat1_vplane/]) diml = dimsizes(lat2_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat2_vplane(0), lon2_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat2_vplane(nl), lon2_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc102 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon2_vplane, lat2_vplane/]) diml= dimsizes(lat3_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat3_vplane(0), lon3_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat3_vplane(nl), lon3_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc103 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon3_vplane, lat3_vplane/]) diml= dimsizes(lat4_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat4_vplane(0), lon4_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat4_vplane(nl), lon4_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc104 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon4_vplane, lat4_vplane/]) ixsec=1 lon_vplane = wrf_user_intrp2d(lon2d,plane1,angles(ixsec),optx) lat_vplane = wrf_user_intrp2d(lat2d,plane1,angles(ixsec),optx) diml = dimsizes(lat_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat_vplane(0), lon_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat_vplane(nl), lon_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc111 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon_vplane, lat_vplane/]) lon_vplane = wrf_user_intrp2d(lon2d,plane2,angles(ixsec),optx) lat_vplane = wrf_user_intrp2d(lat2d,plane2,angles(ixsec),optx) diml = dimsizes(lat_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat_vplane(0), lon_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat_vplane(nl), lon_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc112 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon_vplane, lat_vplane/]) lon_vplane = wrf_user_intrp2d(lon2d,plane3,angles(ixsec),optx) lat_vplane = wrf_user_intrp2d(lat2d,plane3,angles(ixsec),optx) diml = dimsizes(lat_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat_vplane(0), lon_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat_vplane(nl), lon_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc113 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon_vplane, lat_vplane/]) lon_vplane = wrf_user_intrp2d(lon2d,plane4,angles(ixsec),optx) lat_vplane = wrf_user_intrp2d(lat2d,plane4,angles(ixsec),optx) diml = dimsizes(lat_vplane) nl = diml-1 nm = getind_latlon2d (lat2d,lon2d, lat_vplane(0), lon_vplane(0)) yb1 = nm(0,0) xb1 = nm(0,1) nm = getind_latlon2d (lat2d,lon2d, lat_vplane(nl), lon_vplane(nl)) yb2 = nm(0,0) xb2 = nm(0,1) vxsc114 = gsn_add_polyline(wks,contour0,(/xb1, xb2/), (/yb1, yb2/), resb) delete([/lon_vplane, lat_vplane/]) plot = wrf_map_overlays(a1,wks,(/contour0/),pltres1,mpres) dum = new(2,graphic) colors=(/"black","white"/) types = (/16,16/) sizes = (/0.009 ,0.007/) do idum=0,1 pmres = True pmres@gsMarkerColor = colors(idum) pmres@gsMarkerIndex = types(idum) pmres@gsMarkerSizeF = sizes(idum) dum(idum) = gsn_add_polymarker(wks, plot, sitex, sitey, pmres) delete(pmres) end do draw(plot) frame(wks) ;;;end do delete(opts) delete(optx) delete(contour0) delete(lat2d) delete(lon2d) delete(z) ; it loop delete(times) delete(fname1) delete(a1) print(" "+fout+".png") end