external WRF "./wrf_user_fortran_util_0.so" ; Script to produce standard plots for a WRF real-data run, ; with the ARW coordinate dynamics option. load "WRFOptions.ncl" ; set basic plot options here ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "gsn_code.ncl" ;load "/local/LINUX_PC/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "WRFPlot.ncl" load "WRFUserARW.ncl" load "SkewTFunc.ncl" begin ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. a = addfile("../EXAMPLES/wrfout.nc","r") ; We generate plots, but what kind do we prefer? type = "x11" ; type = "pdf" ; type = "ps" ; type = "ncgm" wks = gsn_open_wks(type,"wrf_real") ; Debug information. debug = False ; debug = True icount = 0 ; Basic Plot Information res@MainTitle = "REAL-TIME WRF" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; What times and how many time steps are in the data set? FirstTime = True times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) ; number of times in the file ; The specific pressure levels that we want the data ; interpolated to. pressure_levels = (/ 850., 700., 500., 300./) ; pressure levels to plot nlevels = dimsizes(pressure_levels) ; number of pressure levels ; This is the big loop over all of the time periods to process. do it = 0,ntimes-1,2 ; do it = 0,2,2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Save some Time Information and create the Map Background time = it if (FirstTime) then times_sav = times(it) end if res@TimeLabel = times(it) map = wrf_map(wks,a,res) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; First get the variables we will need slvl = wrf_user_getvar(a,"slvl",time) ; psl wrf_user_filter2d(slvl, 3) ; smooth psl qv = wrf_user_getvar(a,"QVAPOR",time) ; Qv tc = wrf_user_getvar(a,"tc",time) ; T in C td = wrf_user_getvar(a,"td",time) ; Td in C if ( it .eq. 0 ) then tc2 = tc(0,:,:) ; Use lowest T at time zero else tc2 = wrf_user_getvar(a,"T2",time) ; T2 in Kelvin tc2 = tc2-273.16 ; T2 in C end if if (it .eq. 0) then td2 = td(0,:,:) ; Use lowest Td at time zero else td2 = wrf_user_getvar(a,"td2",time) ; Td2 in C end if u = wrf_user_getvar(a,"ua",time) ; u averaged to mass points v = wrf_user_getvar(a,"va",time) ; v averaged to mass points uvm = wrf_user_getvar(a,"umeta",time) ; umet,vmet averaged to mass points if (it .eq. 0) then u10 = u(0,:,:) ; Use lowest level at time 0 v10 = v(0,:,:) else u10 = wrf_user_getvar(a,"U10",time) ; u at 10 m, mass point v10 = wrf_user_getvar(a,"V10",time) ; v at 10 m, mass point end if w = wrf_user_getvar(a,"wa",time) ; vertical velocity, averaged to half levels w = 100.*w ; vert vel in cm/s p = wrf_user_getvar(a, "p",time) ; pressure is our vertical coordinate z = wrf_user_getvar(a, "Z",time) ; grid point height rh = wrf_user_getvar(a,"rh",time) ; relative humidity ter = wrf_user_getvar(a,"HGT",time) ; terrain height ; Get non-convective, convective and total precipitation ; Calculate tendency values rain_exp = wrf_user_getvar(a,"RAINNC",time) rain_con = wrf_user_getvar(a,"RAINC",time) 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 ; Bookkeeping, just to allow the tendency at the next time step rain_exp_save = rain_exp rain_con_save = rain_con rain_tot_save = rain_tot ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ************************************************************ ; 1) Surface Temp + SLP + wind barbs ; Surface Td + wind barbs ; Plotting options for T tc_plane = 1.8*tc2+32. ; Turn temperature into Fahrenheit - Typical of a surface map opts_tc = res opts_tc@FieldTitle = "Surface Temperature" opts_tc@UnitLabel = "F" opts_tc@ContourParameters = (/ -20., 90., 5./) opts_tc@cnFillOn = True opts_tc@gsnSpreadColorEnd = -3 ; End third from the last color in color map contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) print_opts("opts_tc", opts_tc, debug) ; Plotting options for Td td_plane = 1.8*td2+32. ; Turn temperature into Fahrenheit - Typical of a surface map opts_td = res opts_td@FieldTitle = "Surface Dew Point Temp" opts_td@UnitLabel = "F" opts_td@ContourParameters = (/ -20., 90., 5./) opts_td@cnFillOn = True opts_td@cnLinesOn = True opts_td@cnLineLabelsOn = True opts_td@cnLineLabelBackgroundColor = -1 opts_td@gsnSpreadColorEnd = -3 ; End third from the last color in color map contour_td = wrf_contour(a,wks,td_plane,opts_td) print_opts("opts_td", opts_td, debug) ; Plotting options for SLP opts_psl = res ;opts_psl@FieldTitle = "Sea Level Pressure" ;opts_psl@UnitLabel = "mb" opts_psl@ContourParameters = (/ 900., 1100., 4. /) opts_psl@cnLineColor = "NavyBlue" opts_psl@cnHighLabelsOn = True opts_psl@cnLowLabelsOn = True opts_psl@cnLineLabelBackgroundColor = -1 opts_psl@gsnContourLineThicknessesScale = 2.0 contour_psl = wrf_contour(a,wks,slvl,opts_psl) print_opts("opts_psl", opts_psl, debug) ; Plotting options for Wind Vectors u_plane = u10*1.94386 ; Turn wind into knots - Typical of a surface map v_plane = v10*1.94386 opts_vct = res opts_vct@FieldTitle = "Winds" opts_vct@UnitLabel = "kts" opts_vct@NumVectors = 47 vector = wrf_vector(a,wks,u_plane, v_plane,opts_vct) print_opts("opts_vct", opts_vct, debug) ; MAKE PLOTS wrf_map_overlay(wks,map,(/contour_tc,contour_psl,vector/),False) print_header(icount,debug) wrf_map_overlay(wks,map,(/contour_td,vector/),False) print_header(icount,debug) ; Delete options and fields, so we don't have carry over delete(opts_tc) delete(opts_td) delete(opts_vct) delete(opts_psl) delete(tc_plane) delete(td_plane) delete(u_plane) delete(v_plane) ; ************************************************************ ; 2) 2d precip plots ; Plotting options for Sea Level Pressure opts_psl = res opts_psl@ContourParameters = (/ 900., 1100., 2. /) opts_psl@cnLineColor = "Blue" opts_psl@cnInfoLabelOn = False opts_psl@cnHighLabelsOn = True opts_psl@cnLowLabelsOn = True opts_psl@gsnContourLineThicknessesScale = 1.0 contour_psl = wrf_contour(a,wks,slvl,opts_psl) print_opts("opts_psl",opts_psl, debug) ; Plotting options for Precipitation opts_r = res opts_r@UnitLabel = "mm" opts_r@cnLevelSelectionMode = "ExplicitLevels" opts_r@cnLevels = (/ .1, .2, .4, .8, 1.6, 3.2, 6.4, \ 12.8, 25.6, 51.2, 102.4/) opts_r@cnFillColors = (/"White","White","DarkOliveGreen1", \ "DarkOliveGreen3","Chartreuse", \ "Chartreuse3","Green","ForestGreen", \ "Yellow","Orange","Red","Violet"/) opts_r@cnInfoLabelOn = False opts_r@cnConstFLabelOn = False print_opts("opts_r", opts_r, debug) ; MAKE PLOTS opts_r@cnFillOn = True ; Total Precipitation (color fill) opts_r@FieldTitle = "Total Precipitation" contour_tot = wrf_contour(a,wks, rain_tot, opts_r) ; Precipitation Tendency (color fill) opts_r@FieldTitle = "Precipitation Tendency" opts_r@SubFieldTitle = "from " + times_sav + " to " + times(it) contour_tend = wrf_contour(a,wks, rain_tot_tend,opts_r) ; Non-Convective Precipitation Tendency (color fill) opts_r@FieldTitle = "Explicit Precipitation Tendency" opts_r@SubFieldTitle = "from " + times_sav + " to " + times(it) contour_res = wrf_contour(a,wks,rain_exp_tend,opts_r) ; Convective Precipitation Tendency (contour lines) opts_r@FieldTitle = "Param Precipitation Tendency" opts_r@SubFieldTitle = "from " + times_sav + " to " + times(it) opts_r@cnFillOn = False opts_r@cnLineColor = "Red4" contour_prm = wrf_contour(a,wks,rain_con_tend,opts_r) ; Total Precipitation Tendency + SLP wrf_map_overlay(wks,map,(/contour_tend,contour_psl/),False) print_header(icount,debug) ; Total Precipitation wrf_map_overlay(wks,map,contour_tot,False) print_header(icount,debug) ; Non-Convective and Convective Precipiation Tendencies wrf_map_overlay(wks,map,(/contour_res,contour_prm/),False) print_header(icount,debug) ; Delete options, so we don't have carry over delete(opts_r) delete(opts_psl) ; ************************************************************ ; 3) 3d plots: height, rh, vert vel, temp ; Loop over the selected pressure levels to plot do level = 0,nlevels-1 pressure = pressure_levels(level) tc_plane = wrf_user_intrp3d(tc,p,ter,"h",pressure,0.) z_plane = wrf_user_intrp3d( z,p,ter,"h",pressure,0.) rh_plane = wrf_user_intrp3d(rh,p,ter,"h",pressure,0.) u_plane = wrf_user_intrp3d( u,p,ter,"h",pressure,0.) v_plane = wrf_user_intrp3d( v,p,ter,"h",pressure,0.) spd = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec u_plane = u_plane*1.94386 ; kts v_plane = v_plane*1.94386 ; kts ; Plotting options for T opts_tc = res opts_tc@cnLineColor = "Red" opts_tc@ContourParameters = (/ 5.0 /) opts_tc@cnInfoLabelOrthogonalPosF = 0.06 opts_tc@gsnContourLineThicknessesScale = 2.0 ; Plotting options for Geopotential Heigh opts_z = res opts_z@cnLineColor = "NavyBlue" opts_z@gsnContourLineThicknessesScale = 3.0 ; Plotting options for RH opts_rh = res opts_rh@ContourParameters = (/ 10., 90., 10./) opts_rh@cnFillOn = True opts_rh@cnFillColors = (/"White","White","White", \ "White","Chartreuse","Green",\ "Green3","Green4", \ "ForestGreen","PaleGreen4"/) ; Plotting options for Wind Speed opts_spd = res opts_spd@FieldTitle = "Wind Speed" opts_spd@UnitLabel = "m/s" opts_spd@cnLineColor = "MediumSeaGreen" opts_spd@ContourParameters = (/ 10. /) opts_spd@cnInfoLabelOrthogonalPosF = 0.06 opts_spd@gsnContourLineThicknessesScale = 3.0 ; Plotting options for Wind Vectors opts_vct = res opts_vct@FieldTitle = "Winds" opts_vct@UnitLabel = "kts" opts_vct@NumVectors = 47 ; MAKE PLOTS if ( pressure .eq. 850 ) then ; temp, rh, height, barbs opts_z@ContourParameters = (/ 20.0 /) contour_rh = wrf_contour(a,wks,rh_plane,opts_rh) print_opts("opts_rh", opts_rh, debug) contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) print_opts("opts_tc", opts_tc, debug) contour_height = wrf_contour(a,wks, z_plane,opts_z) print_opts("opts_z", opts_z, debug) vector = wrf_vector(a,wks,u_plane, v_plane,opts_vct) print_opts("opts_vct", opts_vct, debug) wrf_map_overlay(wks,map,(/contour_rh,contour_tc,contour_height, \ vector/),False) print_header(icount,debug) end if if ( pressure .eq. 700 ) then ; temp, height, barbs opts_z@ContourParameters = (/ 30.0 /) contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) print_opts("opts_tc", opts_tc, debug) contour_height = wrf_contour(a,wks, z_plane,opts_z) print_opts("opts_z", opts_z, debug) vector = wrf_vector(a,wks,u_plane, v_plane,opts_vct) print_opts("opts_vct", opts_vct, debug) wrf_map_overlay(wks,map,(/contour_tc,contour_height, \ vector/),False) print_header(icount,debug) end if if ( pressure .eq. 500 ) then ; temp, height, barbs opts_z@ContourParameters = (/ 60.0 /) contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) print_opts("opts_tc", opts_tc, debug) contour_height = wrf_contour(a,wks, z_plane,opts_z) print_opts("opts_z", opts_z, debug) vector = wrf_vector(a,wks,u_plane, v_plane,opts_vct) print_opts("opts_vct", opts_vct, debug) wrf_map_overlay(wks,map,(/contour_tc,contour_height, \ vector/),False) print_header(icount,debug) end if if ( pressure .eq. 300 ) then ; windspeed, height, barbs opts_z@ContourParameters = (/ 60.0 /) contour_spd = wrf_contour(a,wks,spd,opts_spd) print_opts("opts_spd", opts_spd, debug) contour_height = wrf_contour(a,wks, z_plane,opts_z) print_opts("opts_z", opts_z, debug) vector = wrf_vector(a,wks,u_plane, v_plane,opts_vct) print_opts("opts_vct", opts_vct, debug) wrf_map_overlay(wks,map,(/contour_spd,contour_height, \ vector/),False) print_header(icount,debug) end if end do ; END OF 3D LOOP ; Delete options and fields, so we don't have carry over delete(opts_tc) delete(opts_z) delete(opts_rh) delete(opts_spd) delete(opts_vct) delete(z_plane) delete(tc_plane) delete(rh_plane) delete(u_plane) delete(v_plane) delete(spd) ; ************************************************************ ; 4) xsection plots do ip = 1, 2 ; | ; | ; Two different xsections, angle=90 is |, angle=0 is ------ ; | ; | if(ip .eq. 1) then angle = 90. Orientation = "Cross-Section Orientation : West-East" else angle = 0. Orientation = "Cross-Section Orientation : South-North" end if ; Options for XY Plots opts_xy = res opts_xy@tiXAxisString = "Number of Grid Points" opts_xy@tiYAxisString = "Height (km)" opts_xy@PlotType = "XY" opts_xy@AspectRatio = 0.75 opts_xy@PlotOrientation = Orientation opts_xy@cnMissingValPerimOn = True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 11 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = ispan(0,100,10) opts_xy@tmYLLabels = ispan(0,20,2) opts_xy@tiXAxisFontHeightF = 0.020 opts_xy@tiYAxisFontHeightF = 0.020 opts_xy@tmXBMajorLengthF = 0.02 opts_xy@tmYLMajorLengthF = 0.02 ; Build planes of data for temp and rh dimsrh = dimsizes(rh) plane = new(2,float) plane = (/ dimsrh(2)/2, dimsrh(1)/2 /) rh_plane = wrf_user_intrp3d(rh,z,ter,"v",plane,angle) tc_plane = wrf_user_intrp3d(tc,z,ter,"v",plane,angle) ; Plotting options for RH opts_rh = opts_xy opts_rh@ContourParameters = (/ 10., 90., 10. /) opts_rh@cnFillOn = True opts_rh@cnFillColors = (/"White","White","White", \ "White","Chartreuse","Green", \ "Green3","Green4", \ "ForestGreen","PaleGreen4"/) ; Plotting options for Temperature opts_tc = opts_xy opts_tc@ContourParameters = (/ 5. /) ; Get the contour info for the rh and temp contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) print_opts("opts_tc", opts_tc, debug) contour_rh = wrf_contour(a,wks,rh_plane,opts_rh) print_opts("opts_rh", opts_rh, debug) ; Overlay the fields wrf_overlay(wks,(/contour_rh,contour_tc/),False) print_header(icount,debug) ; Delete options and fields, so we don't have carry over delete(opts_tc) delete(opts_rh) delete(tc_plane) delete(rh_plane) end do ; end of xsection plot rotation ; ************************************************************ ; 5) skew-T plots u = uvm(0,:,:,:)*1.94386 v = uvm(1,:,:,:)*1.94386 do ip = 12, 13 ; Define a few skew-T plotting options skewtOpts = True dataOpts = True dataOpts@Parcel = 0 dataOpts@WspdWdir = False ; wind speed and dir [else: u,v] dataOpts@PlotWindH = False ; plot wind barbs at h lvls [pibal; special] dataOpts@HspdHdir = True ; wind speed and dir [else: u,v] ; Choose the lat/lon of the skew-T plots ip_lats = (/ 39.77, 37.77, 41.12, 41.30, 39.07, 35.22, \ 37.23, 35.22, 31.95, 32.50, 43.56, 38.94, \ 40.0 /) ip_lons = (/-104.87, -99.97, -100.67, -95.90, -95.62, -97.45, \ -93.38, -101.72, -102.20, -97.30, -116.24, -77.44, \ -82.90 /) ip_locs = (/"Denver (DEN)", "Dodge City (DDC)", "North Platte (LBF)", \ "Omaha (OMA)", "Topeka (TOP)", "Norman (OUN)", \ "Springfield (SGF)", "Amarillo (AMA)", "Midland (MAF)", \ "Fort Worth (FWD)", "Boise (BOI)", "Washington DC (IAD)", \ "Columbus (CMH)"/) ip_locs = ip_locs + ", Valid at " + times(it) locr = wrf_user_find_ij_lat_long(a, ip_lats(ip-1), ip_lons(ip-1)) loc = floattointeger(locr) ; Get the skew-T background skewtOpts@tiMainString = ip_locs(ip-1) skewt_bkgd = skewT_BackGround (wks, skewtOpts) ; Draw the skew-T plot draw (skewt_bkgd) skewT_data = skewT_PlotData(wks, skewt_bkgd, p(:,loc(0), loc(1)), \ tc(:,loc(0), loc(1)), \ td(:,loc(0), loc(1)), \ z(:,loc(0), loc(1)), \ -u(:,loc(0), loc(1)), \ -v(:,loc(0), loc(1)), \ dataOpts) ; Close the frame frame(wks) print_header(icount,debug) ; Delete the temporary arrays delete(skewtOpts) delete(dataOpts) delete(skewT_data) delete(skewt_bkgd) end do ; ************************************************************ times_sav = times(it) FirstTime = False end do ; end of the time loop end