;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This is to plot 4 figures at one page, i.e., panel plot. ; ; mchen, 2012/10 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ; load "./WRFUserARW.ncl" begin ; define experiment. exp_num=1 plot the first exp. exp_num=2 plot the two exp. If three more experiments needs to ; plot here, line_modes, line_markers, line_colors, legend_labels are needed to change to match exp_num. titlenum=4 titlename=(/"Total Precip","Surface T and WIND","300mb Field","500mb Fields"/) plot=new(4,graphic) wks = gsn_open_wks("x11","2010061200") gsn_define_colormap(wks,"wh-bl-gr-ye-re") ; Basic Plot Information resB = True resB@MainTitle = " " resB@Footer = False pltres = True ; pltres@PanelPlot = True ; pltres@gsnDraw = False pltres@gsnFrame = False pltres@gsnMaximize = True mpres = True DATADir = "./" FILES = systemfunc (" ls -1 " + DATADir + "wrfout_d01* ") numFILES = dimsizes(FILES) print("numFILES = " + numFILES) print(FILES) pressure_levels = (/ 500., 300./) ; pressure levels to plot nlevels = dimsizes(pressure_levels) ; number of pressure levels FirstTime = True do ifil = 0, numFILES-1 ; BIG FILES LOOP a = addfile(FILES(ifil)+".nc","r") times = wrf_user_getvar(a,"times",-1) ; get all times in the file ntimes = dimsizes(times) ; number of times in the file print("Working on time: " + times(0) ) resB@TimeLabel = times(0) ; Set Valid time to use on plots if (FirstTime) then times_sav = times(0) end if res = resB res@cnFillOn = False it=0 slp = wrf_user_getvar(a,"slp",it) ; psl wrf_smooth_2d(slp, 3) ; smooth psl td2 = wrf_user_getvar(a,"td2",it) ; Td2 in C tc2 = wrf_user_getvar(a,"T2",it) ; T2 in Kelvin tc2 = tc2-273.16 ; T2 in C u10 = wrf_user_getvar(a,"U10",it) ; u at 10 m, mass point v10 = wrf_user_getvar(a,"V10",it) ; v at 10 m, mass point tf2 = 1.8*tc2+32. ; Turn temperature into Fahrenheit tf2@description = "Surface Temperature" tf2@units = "F" td_f = 1.8*td2+32. ; Turn temperature into Fahrenheit td_f@description = "Surface Dew Point Temp" td_f@units = "F" u10 = u10*1.94386 ; Turn wind into knots v10 = v10*1.94386 u10@units = "kts" v10@units = "kts" tc = wrf_user_getvar(a,"tc",it) ; 3D tc u = wrf_user_getvar(a,"ua",it) ; 3D U at mass points v = wrf_user_getvar(a,"va",it) ; 3D V at mass points p = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate z = wrf_user_getvar(a, "z",it) ; grid point height rh = wrf_user_getvar(a,"rh",it) ; relative humidity rain_exp = wrf_user_getvar(a,"RAINNC",it) rain_con = wrf_user_getvar(a,"RAINC",it) snow_exp = wrf_user_getvar(a,"SNOWNC",it) rain_tot = rain_exp + rain_con if( FirstTime ) then rain_exp_save = rain_exp rain_con_save = rain_con rain_tot_save = rain_tot end if rain_exp_tend = rain_exp - rain_exp_save rain_con_tend = rain_con - rain_con_save rain_tot_tend = rain_tot - rain_tot_save rain_exp_save = rain_exp rain_con_save = rain_con rain_tot_save = rain_tot do level = 0,nlevels-1 ; LOOP OVER LEVELS pressure = pressure_levels(level) tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False) z_plane = wrf_user_intrp3d( z,p,"h",pressure,0.,False) rh_plane = wrf_user_intrp3d(rh,p,"h",pressure,0.,False) u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False) v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False) spd = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec spd@description = "Wind Speed" spd@units = "m/s" u_plane = u_plane*1.94386 ; kts v_plane = v_plane*1.94386 ; kts u_plane@units = "kts" v_plane@units = "kts" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Plotting options for T opts = res opts@cnFillOn = True opts@ContourParameters = (/ -65., 5., 5./) opts@gsnSpreadColorEnd = -3 ; End third from the last color in color map contour_tc = wrf_contour(a,wks,tc_plane,opts) delete(opts) ; Plotting options for RH opts = res opts@cnFillOn = True opts@pmLabelBarOrthogonalPosF = -0.1 opts@ContourParameters = (/ 10., 90., 10./) opts@cnFillColors = (/"White","White","White", \ "White","Chartreuse","Green",\ "Green3","Green4", \ "ForestGreen","PaleGreen4"/) contour_rh = wrf_contour(a,wks,rh_plane,opts) delete(opts) ; Plotting options for Wind Speed opts = res opts@cnLineColor = "MediumSeaGreen" opts@ContourParameters = (/ 10. /) opts@cnInfoLabelOrthogonalPosF = 0.07 ; offset second label information opts@gsnContourLineThicknessesScale = 3.0 contour_spd = wrf_contour(a,wks,spd,opts) delete(opts) ; Plotting options for Wind Vectors opts = res opts@FieldTitle = "Wind" ; overwrite Field Title opts@NumVectors = 47 ; wind barb density opts@vcRefAnnoArrowLineColor = "white" ; change ref vector color opts@vcLineArrowColor = "white" ; change vector color vector = wrf_vector(a,wks,u_plane,v_plane,opts) delete(opts) ; Plotting options for Geopotential Heigh opts_z = res opts_z@cnLineColor = "Red" opts_z@gsnContourLineThicknessesScale = 3.0 ; MAKE PLOTS if ( pressure .eq. 850 ) then ; plot temp, rh, height, wind barbs opts_z@ContourParameters = (/ 20.0 /) contour_height = wrf_contour(a,wks,z_plane,opts_z) plot(level) = wrf_map_overlays(a,wks,(/contour_rh,contour_tc,contour_height, \ vector/),pltres,mpres) end if if ( pressure .eq. 700 ) then ; plot temp, height, wind barbs opts_z@ContourParameters = (/ 30.0 /) contour_height = wrf_contour(a,wks, z_plane,opts_z) plot(level) = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \ vector/),pltres,mpres) end if if ( pressure .eq. 500 ) then ; plot temp, height, wind barbs opts_z@ContourParameters = (/ 60.0 /) contour_height = wrf_contour(a,wks, z_plane,opts_z) plot(level) = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \ vector/),pltres,mpres) end if if ( pressure .eq. 300 ) then ; plot windspeed, height, wind barbs opts_z@ContourParameters = (/ 60.0 /) contour_height = wrf_contour(a,wks, z_plane,opts_z) plot(level) = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \ vector/),pltres,mpres) end if delete(opts_z) end do ; END OF LEVEL LOOP ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Plotting options for Sea Level Pressure opts = res opts@ContourParameters = (/ 900., 1100., 10. /) opts@cnLineColor = "Blue" opts@cnInfoLabelOn = False ;opts@cnHighLabelsOn = True ;opts@cnLowLabelsOn = True opts@gsnContourLineThicknessesScale = 1.0 contour_psl = wrf_contour(a,wks,slp,opts) delete(opts) ; Plotting options for Precipitation opts = res opts@cnLevelSelectionMode = "ExplicitLevels" opts@cnLevels = (/ .1, .2, .4, .8, 1.6, 3.2, 6.4, \ 12.8, 25.6, 51.2, 102.4/) opts@cnFillColors = (/"White","White","DarkOliveGreen1", \ "DarkOliveGreen3","Chartreuse", \ "Chartreuse3","Green","ForestGreen", \ "Yellow","Orange","Red","Violet"/) opts@cnInfoLabelOn = False opts@cnConstFLabelOn = False ; Precipitation (color fill) opts@cnFillOn = True opts@FieldTitle = "Total Precipitation" contour_tot = wrf_contour(a,wks, rain_tot, opts) delete(opts) ; MAKE PLOTS plot(2) = wrf_map_overlays(a,wks,contour_tot,pltres,mpres) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Plotting options for T opts = res opts@cnFillOn = True opts@ContourParameters = (/ -20., 90., 5./) opts@gsnSpreadColorEnd = -3 ; End third from the last color in color map contour_tc = wrf_contour(a,wks,tf2,opts) delete(opts) ; Plotting options for SLP opts = res opts@cnLineColor = "Blue" opts@cnHighLabelsOn = True opts@cnLowLabelsOn = True opts@ContourParameters = (/ 900., 1100., 4. /) opts@cnLineLabelBackgroundColor = -1 opts@gsnContourLineThicknessesScale = 2.0 contour_psl = wrf_contour(a,wks,slp,opts) delete(opts) ; Plotting options for Wind Vectors opts = res opts@FieldTitle = "Wind" ; overwrite Field Title opts@NumVectors = 47 ; density of wind barbs vector = wrf_vector(a,wks,u10,v10,opts) delete(opts) ; MAKE PLOTS plot(3) = wrf_map_overlays(a,wks,(/contour_tc,contour_psl,vector/),pltres,mpres) end do ; END IF BIG FILE LOOP ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; panel plot pnlres = True ; plot mods desired pnlres@gsnPanelYWhiteSpacePercent = 5 pnlres@gsnPanelXWhiteSpacePercent = 5 gsn_panel(wks,(/plot/),(/2,2/),pnlres) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end