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/csm/shea_util.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Script to read in WRF output and make plots. This script can plot one field in color, one field with ; contours, and wind barbs. The colored field is "var1" and is grouped with the wind barbs when ; plotting, and the contoured field is "var2". ; ; THIS SCRIPT MAKES PLOTS OF SIMULATED LOW-LEVEL ROTATION (enter specifics here) OVERLAID WITH MRMS ; OBSERVED LOW-LEVEL ROTATION. ; ; NOTE: This script is based off of the original "create_xyplots_basevars_mem_15kmgrid_varinput.ncl ; but reads in MRMS data to plot observed variables over simulated variables (reflectivity or ; rotation). ; ; NOTE: This script is designed specifically for final analysis of Member 30 (0400 - 0530 UTC) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; begin ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; User Input hhmm_all = new( (/11/) , "string" , -999.0) yymmdd_all = new( (/11/) , "string" , -999.0) start = 0423 step = 0 pi = 3.14159265 ;---------------------------------------------------------; ; This do loop creates the times to plot output. Currently ; set from 0400 - 0530 UTC every 5 mins. ;---------------------------------------------------------; do i = 0,10 ; number of output times for each 1km ensemble member (46 from 0413 - 0459 UTC... every 2 minutes) ; ; add_on = start + step ; hhmm_all(i) = "0" + add_on ; ; yymmdd_all(i) = "150706" ; ; if (step .eq. 55) then ; step = step - 55 ; start = start + 100 ; else ; step = step + 5 ; end if yymmdd_all(i) = "150706" end do hhmm_all = (/"0423","0425","0427","0429","0431","0433","0435","0437","0439","0441","0443"/) mm_all_1km = (/"20","25","25","30","30","35","35","35","40","40","45"/) ; mn = "00" ; valid minute ; commented out because minutes are included in the line above ff = "00" ; forecast in minutes: 00 or 30 expm = "CTRL" ; "CTRL","... clat = 43.531399 ; center of the plot grid clon = -97.875127 ; center of the plot grid domain = 1 ; 1 = 1 km (currently) var1 = "w_levels" ; Shaded contour: "temp","dewp","rh","pmsl","theta" (only 2m),"gph","thetae","vort","cape","refl","rot" lvl1 = "2m" ; "2m" or any specified pressure level (like "850" or "500", for example) var2 = "rot" ; Non-shaded overlay contour...same vars available as above as well as "cin" and "max_uh", and "mrms_refl", "mrms_rot" lvl2 = "2m" plot_width = 0.20 ; degrees longitude (5 /8/30 for 1/3/15 km grid) NOTE: Keep the aspect ratio (width/height) 1.666! plot_height = 0.12 ; degrees latitude (2.5/4/18 for 1/3/15 km grid) tornado_loc_lat = 43.507499 tornado_loc_lon = -97.926625 ; End Input ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; setvalues NhlGetWorkspaceObjectId() "wsMaximumSize": 1000000000 end setvalues days = (/31,28,31,30,31,30,31,31,30,31,30,31/) n_times = dimsizes(hhmm_all) ; number of analysis times to produce do nens = 30,30 ; 1,36 ; work on Member 30 <-- ** THE REST OF THE CODE IS IN THIS DO LOOP ** print("NOW WORKING ON ENSEMBLE MEMBER " + nens + "/36 ... duh") do t = 5,6 ; 0 , n_times-1 ; <-- loop through analysis times (0400 - 0530 UTC) ;------------------------------------------------- ; These are the variables that depend on t (time). ;------------------------------------------------- yymmdd = yymmdd_all(t) hhmm = hhmm_all(t) yymmddc = stringtocharacter(yymmdd) yy = charactertostring(yymmddc(0:1)) mm = charactertostring(yymmddc(2:3)) dd = charactertostring(yymmddc(4:5)) delete(yymmddc) hhmmc = stringtocharacter(hhmm) hh = charactertostring(hhmmc(0:1)) mn = charactertostring(hhmmc(2:3)) delete(hhmmc) imm = stringtoint(mm) idd = stringtoint(dd) ihh = stringtoint(hh) imn = stringtoint(mn) dir = "/scratch/matt.flournoy/WRF/data/333m/fcsts_mem" + nens ; <-- ** THIS IS WHERE YOU SET THE INPUT DIRECTORY ** ; if( ff.eq."00" ) then <-- commented stuff in this paragraph out since these aren't really forecasts wrf_prefix = "wrfout" ; end if ; if( ff.eq."30" ) then ; wrf_prefix = "wrffcst" ; end if ; if ( ff.eq."30" ) then <-- commented stuff in this paragraph out too ; imn=0 ; ihh=ihh+1 ; if ( ihh.ge.24 ) then ; ihh=ihh-24 ; idd=idd+1 ; if( idd.gt.days(imm) ) then ; imm = imm+1 ; idd = 1 ; end if ; end if ; end if mmv = sprinti("%0.2i",imm) ddv = sprinti("%0.2i",idd) hhv = sprinti("%0.2i",ihh) mnv = sprinti("%0.2i",imn) ; Read Ensemble members wrf_name = wrf_prefix + "_d0" + domain + "_20" + yy + "-" + mmv + "-" + ddv + "_" + hhv + ":" + mnv + ":00" ; system("rm -f ncfile.nc") <-- commented the system and "ncfile.nc" stuff out so I can run this script in multiple windows at once ; ncfile = dir+"/"+wrf_name ; print(ncfile) ; system("ln -sf "+ncfile+" ncfile.nc") filename = dir + "/" + wrf_name ens = addfile(filename,"r") ; system("rm -f ncfile.nc") print(wrf_name) loc = wrf_user_latlon_to_ij(ens, clat, clon) if ( ismissing(loc(0)) .or. ismissing(loc(1)) ) print("Center point outside model domain" ) exit end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Extract the variables we need for the plot ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; heights = wrf_user_getvar(ens,"height",0) ; print(heights(:,349,399)) ; printVarSummary(heights) ; sleep(456) p0 = wrf_user_getvar(ens,"pressure",0) p0d = dimsizes(p0) nlvls = p0d(0) ny = p0d(1) nx = p0d(2) delete(p0) delete(p0d) latgrid = new( (/ny,nx/) , "float" ) longrid = new( (/ny,nx/) , "float" ) ; plotvar1 = new( (/ny,nx/) , "float" , -999.0 ) plotvar2 = new( (/ny,nx/) , "float" , -999.0 ) plotvar4 = new( (/ny,nx/) , "float" , -999.0 ) uuplot = new( (/ny,nx/) , "float" , -999.0 ) vvplot = new( (/ny,nx/) , "float" , -999.0 ) hgt_agl = new( (/50,ny,nx/) , "float" , -999.0) lats_mrms = fspan(50.00,40.00,1001) ; lats/lons for MRMS reflectivity data (1.0 degree resolution) lons_mrms = fspan(-102.00,-95.00,701) lats_mrms_rot = fspan(51.00,21.00,6000) ; lats/lons for MRMS rotation data (0.5 degree resolution) lons_mrms_rot = fspan(-127.00,-65.00,12400) ; plotvar2_mrms = new( (/3500,7000/) , "float" , -999.0) ; Special array used for reading in "mrms_refl" and "mrms_rot" for the second variable plotvar1_all = new( (/ny,nx/) , "float" , -999.0 ) ; Special array used for masking refl plotvar1_all_avo = new( (/1,ny,nx/) , "float" , -999.0) ; Special array used for reading in absolute vorticity plotvar1_cape = new( (/4,ny,nx/) , "float" , -999.0 ) ; Special array used for retrieving cape and cin (and lcl/lfc) plotvar2_cape = new( (/4,ny,nx/) , "float" , -999.0 ) ; Special array used for retrieving cape and cin (and lcl/lfc) ; plotvar1_pre = new( (/ny,nx/) , "float" , -999.0 ) ; Special array used for masking cape (see lines 169 & 170) ; plotvar2_pre = new( (/ny,nx/) , "float" , -999.0 ) ; Special array used for masking cin (see lines 169 & 170) latgrid(:,:) = wrf_user_getvar(ens,"lat",0) ; Model latitude longrid(:,:) = wrf_user_getvar(ens,"lon",0) ; Model longriditude ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Extract the first variable ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if ( lvl1.eq."2m" ) then ;Extract the 2 m variable if ( var1.eq."temp" ) then var1title = "temp (c) at 2 m" print("Extracting variables for ensemble member "+nens) plotvar1(:,:) = wrf_user_getvar(ens,"T2",0) plotvar1 = plotvar1 - 273.16 end if if ( var1.eq."dewp" ) then var1title = "dewp (c) at 2 m" print("Extracting variables for ensemble member "+nens) plotvar1(:,:) = wrf_user_getvar(ens,"td2",0) end if if ( var1.eq."rh" ) then var1title = "rh (%) at 2 m" print("Extracting variables for ensemble member "+nens) plotvar1(:,:) = wrf_user_getvar(ens,"rh2",0) end if if ( var1.eq."pmsl" ) then var1title = "pmsl (mb)" print("Extracting variables for ensemble member "+nens) plotvar1(:,:) = wrf_user_getvar(ens,"slp",0) end if if ( var1.eq."theta" ) then var1title = "pot temp (K) at lowest lvl" print("Extracting variables for ensemble member "+nens) theta = new( (/nlvls,ny,nx/) , "float" , -999.0 ) theta = wrf_user_getvar(ens,"theta",0) plotvar1(:,:) = theta(0,:,:) delete(theta) end if if ( var1.eq."cape" ) then var1title = "MLCAPE (J/kg)" print("Extracting variables for ensemble member "+nens) plotvar1_cape(:,:,:) = wrf_user_getvar(ens,"cape_2d",0) plotvar1(:,:) = plotvar1_cape(0,:,:) end if if ( var1.eq."refl" ) then var1title = "Sim. Max Reflectivity (dBZ)" print("Extracting variables for ensemble member "+nens) plotvar1_all(:,:) = wrf_user_getvar(ens,"mdbz",0) plotvar1(:,:) = mask(plotvar1_all,(plotvar1_all.ge.20),True) ; <-- masks all reflectivity below 20 dBZ for plotting purposes end if if ( var1.eq."rot" ) then var1title = "Sim. Low-Level Rotation (10~S~-5~N~ s~S~-1~N~)" plotvar1_all_abs_vort = wrf_user_getvar(ens,"avo",0) ; Absolute Vorticity (10-5 s-1) ... all levels, lats, lons for that time plotvar1_all_avo = plotvar1_all_abs_vort(0,:,:) ; Absolute Vorticity at the lowest model level, all lats/lons for that time ; plotvar1(:,:) = dim_avg_n(plotvar1_all_avo,0) ; Absolute Vorticity average in the first 10 levels at all lats/lons for that time plotvar1(:,:) = plotvar1_all_avo - 2 * 0.000072921159 * sin(43.436648*pi/180) ; convert to relative vorticity ; plotvar1(:,:) = mask(plotvar1,(plotvar1.ge.50),True) ; <-- masks all absolute vorticity below 50*10e-5 for plotting purposes end if if ( var1.eq."rot_horiz" ) then ; <-- use this to calculate 0-1 km AGL horizontal vorticity tilting var1title = "Horizontal Vorticity Tilting (10~S~-5~N~ s~S~-2~N~)" height_w = 500 height_u_1km = 1000 height_v_1km = 1000 print("defining tilting variable...") tilting = new( (/ny,nx/) , "float" , -999.0) ; <-- define tilting array u_horiz_vort = new( (/ny,nx/) , "float" , -999.0) ; <-- define x component of horizontal vorticity v_horiz_vort = new( (/ny,nx/) , "float" , -999.0) ; <-- define y component of horizontal vorticity print("done defining tilting variable, now extracting u, v, and w...") ; plotvar1_all_abs_vort = wrf_user_getvar(ens,"REL_VORT",0) ; Relative Vorticity (10-5 s-1) ... all levels, lats, lons for that time u_alls = wrf_user_getvar(ens,"U",0) ; <-- extract u on mass points from the file v_alls = wrf_user_getvar(ens,"V",0) ; <-- extract v on mass points from the file w_alls = wrf_user_getvar(ens,"W",0) ; <-- extract w on mass points from the file u_all = u_alls(:,:,0:579) ; <-- get rid of one index so they're all on a (50 x 350 x 400) grid v_all = v_alls(:,0:409,:) ; <-- get rid of one index so they're all on a (50 x 350 x 400) grid w_all = w_alls(0:49,:,:) ; <-- get rid of one index so they're all on a (50 x 350 x 400) grid z = wrf_user_getvar(ens,"z",0) ; <-- extract total z on mass points from the file (height AGL + height above sea level) ter_hgt = wrf_user_getvar(ens,"HGT",0) ; <-- extract terrain height at all lat/lon points print("now calculating height AGL") do blah = 0,49 hgt_agl(blah,:,:) = z(blah,:,:) - ter_hgt(:,:) ; <-- height AGL at all lat/lon/level points end do print("done calculating height AGL") ; printVarSummary(z) ; printVarSummary(ter_hgt) ; print(hgt_agl(0:9,200,200)) ; print(hgt_agl(0:9,100,300)) ; plotvar1_all_avo = plotvar1_all_abs_vort(0,:,:) ; Relative Vorticity at the lowest model level, all lats/lons for that time ; plotvar1(:,:) = dim_avg_n(plotvar1_all_avo,0) ; Absolute Vorticity average in the first 10 levels at all lats/lons for that time ; rel_vort_500m = wrf_user_intrp3d(plotvar1_all_abs_vort,z,"h",height_vort,0.,False); <-- interpolate relative vorticity to 500 m, output to rel_vort_500m (1/s) ; plotvar1(:,:) = plotvar1_all_avo - 2 * 0.000072921159 * sin(43.436648*pi/180) ; convert to relative vorticity ; plotvar1(:,:) = mask(plotvar1,(plotvar1.ge.50),True) ; <-- masks all absolute vorticity below 50*10e-5 for plotting purposes w_mid = w_all(1,:,:) ; <-- extract 500 m AGL w v_top = v_all(2,:,:) ; <-- extract 1km AGL v v_bot = v_all(0,:,:) ; <-- extract surface v u_top = u_all(2,:,:) ; <-- extract 1km AGL u u_bot = u_all(0,:,:) ; <-- extract surface u print("done extracting winds, now calculating tilting...") do m = 1,ny-2 ; <-- loop over latitudes (except for the edges) do n = 1,nx-2 ; <-- loop over longitudes (except for the edges) dwdx = ( w_mid(m,n+1) - w_mid(m,n-1) ) / 666 dwdy = ( w_mid(m+1,n) - w_mid(m-1,n) ) / 666 dvdz = ( v_top(m,n) - v_bot(m,n) ) / ( z(2,m,n) - z(0,m,n) ) dudz = ( u_top(m,n) - u_bot(m,n) ) / ( z(2,m,n) - z(0,m,n) ) tilting(m,n) = dwdx * ( dwdy - dvdz ) + dwdy * (dudz - dwdx) u_horiz_vort(m,n) = dwdy - dvdz ; <-- calculate the magnitude of the x component of horizontal vorticity v_horiz_vort(m,n) = dudz - dwdx ; <-- calculate the magnitude of the y component of horizontal vorticity delete(dwdx) delete(dwdy) delete(dvdz) delete(dudz) end do ; <-- end loop over longitudes end do ; <-- end loop over latitudes plotvar1(:,:) = tilting * 100000 ; <-- convert to 10^5 print("0-1 km AGL relative vertical vorticity stretching calculated") delete(z) delete(u_alls) delete(u_all) delete(v_alls) delete(v_all) delete(w_alls) delete(w_all) delete(w_mid) delete(v_top) delete(v_bot) delete(u_top) delete(u_bot) end if if ( var1.eq."p_dyn_NL" ) then print("extracting dynamic, non-linear pressure perturbation for mvS") print("first, define variables...") p_dyn_nl = new( (/ny,nx/) , "float" , -999.0) p_dyn_nl_splat = new( (/ny,nx/) , "float" , -999.0) p_dyn_nl_spin = new( (/ny,nx/) , "float" , -999.0) u_mean = new( (/50/) , "float" , -999.0) v_mean = new( (/50/) , "float" , -999.0) ; w_mean = new( (/50/) , "float" , -999.0) u_prime = new( (/50,ny,nx/) , "float" , -999.0) v_prime = new( (/50,ny,nx/) , "float" , -999.0) w_prime = new( (/50,ny,nx/) , "float" , -999.0) print("next, extract wind fields...") z = wrf_user_getvar(ens,"z",0) ; <-- extract total z on mass points from the file (height AGL + height above sea level) ter_hgt = wrf_user_getvar(ens,"HGT",0) ; <-- extract terrain height at all lat/lon points u_alls = wrf_user_getvar(ens,"U",0) ; <-- extract u on mass points from the file v_alls = wrf_user_getvar(ens,"V",0) ; <-- extract v on mass points from the file w_alls = wrf_user_getvar(ens,"W",0) ; <-- extract w on mass points from the file u_all = wrf_user_unstagger(u_alls,u_alls@stagger) ; <-- unstagger v_all = wrf_user_unstagger(v_alls,v_alls@stagger) ; <-- unstagger w_all = wrf_user_unstagger(w_alls,w_alls@stagger) ; <-- unstagger print("wind fields extracted and unstaggered, now find mean and prime fields...") ;---------------------------------------------------------; ; Now find the indices for the search box on the 1km grid. ;---------------------------------------------------------; SW_lat = 42.143829 SW_lon = -96.071137 NE_lat = 43.429759 NE_lon = -95.774507 print("finding indices for the search box...") filename_1km = "/scratch/matt.flournoy/WRF/data/1km/fcsts_no1kmDA_0300UTC/mem30/wrfout_d01_2015-07-06_04:" + mm_all_1km(t) + ":00" data = addfile(filename_1km,"r") lat2d_1km = data->XLAT(0,:,:) lon2d_1km = data->XLONG(0,:,:) u_1km_alls = wrf_user_getvar(data,"U",0) v_1km_alls = wrf_user_getvar(data,"V",0) w_1km_alls = wrf_user_getvar(data,"W",0) u_1km_all = wrf_user_unstagger(u_1km_alls,u_1km_alls@stagger) ; <-- unstagger v_1km_all = wrf_user_unstagger(v_1km_alls,v_1km_alls@stagger) ; <-- unstagger w_1km_all = wrf_user_unstagger(w_1km_alls,w_1km_alls@stagger) ; <-- unstagger lats_1km = (/SW_lat,NE_lat/) lons_1km = (/SW_lon,NE_lon/) ij_1km = getind_latlon2d(lat2d_1km,lon2d_1km,lats_1km,lons_1km) SW_lat_ind_1km = ij_1km(0,0) SW_lon_ind_1km = ij_1km(0,1) NE_lat_ind_1km = ij_1km(1,0) NE_lon_ind_1km = ij_1km(1,1) print("indices for the search box found, now calculate mean and prime winds...") print(SW_lat_ind_1km) print(SW_lon_ind_1km) print(NE_lat_ind_1km) print(NE_lon_ind_1km) ; do asdf = 0,49 u_mean = dim_avg_n(u_1km_all(:,SW_lat_ind_1km:NE_lat_ind_1km,SW_lon_ind_1km:NE_lon_ind_1km), (/1,2/) ) v_mean = dim_avg_n(v_1km_all(:,SW_lat_ind_1km:NE_lat_ind_1km,SW_lon_ind_1km:NE_lon_ind_1km), (/1,2/) ) ; end do print("mean u and v calculated, now calculate primes...") do fdsa = 0,49 u_prime(fdsa,:,:) = u_all(fdsa,:,:) - u_mean(fdsa) v_prime(fdsa,:,:) = v_all(fdsa,:,:) - v_mean(fdsa) w_prime(fdsa,:,:) = w_all(fdsa,:,:) ; assume mean w is zero end do print("MEAN AND PRIME WIND FIELDS CALCULATED!") print("now calculating the nonlinear, dynamic pressure perturbation... at ~277 m AGL (index 4)") do m = 1,ny-2 ; <-- loop over latitudes (except for the edges) do n = 1,nx-2 ; <-- loop over longitudes (except for the edges) dudx = ( u_prime(4,m,n+1) - u_prime(4,m,n-1) ) / 666 dudy = ( u_prime(4,m+1,n) - u_prime(4,m-1,n) ) / 666 dudz = ( u_prime(5,m,n) - u_prime(3,m,n) ) / ( z(5,m,n) - z(3,m,n) ) dvdx = ( v_prime(4,m,n+1) - v_prime(4,m,n-1) ) / 666 dvdy = ( v_prime(4,m+1,n) - v_prime(4,m-1,n) ) / 666 dvdz = ( v_prime(5,m,n) - v_prime(3,m,n) ) / ( z(5,m,n) - z(3,m,n) ) dwdx = ( w_prime(4,m,n+1) - w_prime(4,m,n-1) ) / 666 dwdy = ( w_prime(4,m+1,n) - w_prime(4,m-1,n) ) / 666 dwdz = ( w_prime(5,m,n) - w_prime(3,m,n) ) / ( z(5,m,n) - z(3,m,n) ) p_dyn_nl_splat(m,n) = (0.25) * ( 4 * ( (dudx)^2 + (dvdy)^2 + (dwdz)^2 ) + 2 * (dudy + dvdx)^2 + 2 * (dudz + dwdx)^2 + 2 * (dwdy + dvdz)^2 ) p_dyn_nl_spin(m,n) = -(0.5) * ( (dwdy-dvdz)^2 + (dudz-dwdx)^2 + (dvdx-dudy)^2 ) p_dyn_nl(m,n) = p_dyn_nl_splat(m,n) + p_dyn_nl_spin(m,n) delete(dudx) delete(dudy) delete(dudz) delete(dvdx) delete(dvdy) delete(dvdz) delete(dwdx) delete(dwdy) delete(dwdz) end do ; <-- end loop over longitudes end do ; <-- end loop over latitudes print("now taking the inverse Laplacian...") ; inds_missing = ind(ismissing(p_dyn_nl)) ; ; print(inds_missing) ; sleep(345) p_dyn_nl_nomissing = p_dyn_nl(1:ny-2,1:nx-2) p_dyn_nl_invlap = p_dyn_nl_nomissing ilapsg(p_dyn_nl_nomissing,0,p_dyn_nl_invlap) print("done taking the inverse Laplacian! Convert to hPa.") plotvar1 = p_dyn_nl_invlap / 100 print("done calculating nonlinear, dynamic pressure perturbation!!") end if if ( var1.eq."helicity" ) then var1title = "0-1 km AGL SR Helicity (m~S~2~N~ s~S~-2~N~)" plotvar1_all(:,:) = wrf_user_getvar(ens,(/"helicity","1000"/),0) plotvar1(:,:) = plotvar1_all(:,:) end if if ( var1.eq."w_levels" ) then var1title = "Updraft/Downdraft at Various Levels" levels = (/1000,2000,3000,4000,5000/) num_levels = dimsizes(levels) w_plot = new( (/num_levels,ny,nx/) , "float" , -999.0 ) print("extracting w at various levels to plot") u_alls = wrf_user_getvar(ens,"U",0) v_alls = wrf_user_getvar(ens,"V",0) w_alls = wrf_user_getvar(ens,"W",0) u_all = wrf_user_unstagger(u_alls,u_alls@stagger) ; <-- unstagger v_all = wrf_user_unstagger(v_alls,v_alls@stagger) ; <-- unstagger w_all = wrf_user_unstagger(w_alls,w_alls@stagger) ; <-- unstagger z = wrf_user_getvar(ens,"z",0) ; <-- extract total z on mass points from the file (height AGL + height above sea level) print("w extracted and ustaggered, now interpolating to desired levels...") do i = 0 , num_levels - 1 w_plot(i,:,:) = wrf_user_intrp3d(w_all,z,"h",levels(i),0.,False) end do print("done interpolating w!") printVarSummary(w_plot) end if if ( var1.eq."thetae" ) then var1title = "epot temp (K) at lowest lvl" print("Extracting variables for ensemble member "+nens) psfc = new( (/1,ny,nx/) , "float" , -999.0 ) t2 = new( (/1,ny,nx/) , "float" , -999.0 ) td2 = new( (/1,ny,nx/) , "float" , -999.0 ) mr = new( (/1,ny,nx/) , "float" , -999.0 ) thetae = new( (/1,ny,nx/) , "float" , -999.0 ) psfc(0,:,:) = wrf_user_getvar(ens,"PSFC",0) t2(0,:,:) = wrf_user_getvar(ens,"T2",0) td2(0,:,:) = wrf_user_getvar(ens,"td2",0) mr(0,:,:) = mixhum_ptd( psfc,td2+273.16,1 ) thetae = wrf_eth(mr,t2,psfc) plotvar1(:,:) = thetae(0,:,:) delete(psfc) delete(t2) delete(td2) delete(mr) delete(thetae) end if ; Extract the surface winds ; u_all = wrf_user_getvar(ens,"U",0) ; <-- extract 3d u field for this time ; v_all = wrf_user_getvar(ens,"V",0) ; <-- extract 3d v field for this time ; uuplot = u_all(0,:,0:579)*1.94386 ; <-- extract surface u field for this time ; vvplot = v_all(0,0:409,:)*1.94386 ; <-- extract surface v field for this time ; delete(u_all) ; delete(v_all) uuplot = u_all(0,:,:)*1.94386 vvplot = v_all(0,:,:)*1.94386 delete(u_all) delete(v_all) else ;Extract the variable for this pressure level plvl = tofloat(lvl1) if ( var1.eq."temp" ) then var1title = "temp (C) at "+lvl1+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) tc = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) tc = wrf_user_getvar(ens,"tc",0) plotvar1(:,:) = wrf_user_intrp3d(tc,p,"h",plvl,0.,False) plotvar1(:,:) = where(psfc.ge.plvl,plotvar1(:,:),plotvar1@_FillValue) delete(tc) delete(p) delete(psfc) end if if ( var1.eq."dewp" ) then var1title = "dewp (C) at "+lvl1+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) td = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) td = wrf_user_getvar(ens,"td",0) plotvar1(:,:) = wrf_user_intrp3d(td,p,"h",plvl,0.,False) plotvar1(:,:) = where(psfc.ge.plvl,plotvar1(:,:),plotvar1@_FillValue) delete(td) delete(p) delete(psfc) end if if ( var1.eq."rh" ) then var1title = "rh (C) at "+lvl1+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) rh = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) rh = wrf_user_getvar(ens,"rh",0) plotvar1(:,:) = wrf_user_intrp3d(rh,p,"h",plvl,0.,False) plotvar1(:,:) = where(psfc.ge.plvl,plotvar1(:,:),plotvar1@_FillValue) delete(rh) delete(p) delete(psfc) end if if ( var1.eq."gph" ) then var1title = "geo height (m) at "+lvl1+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) gph = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) gph = wrf_user_getvar(ens,"geopt",0) ; Full model geopotential (m2 s-2) plotvar1(:,:) = wrf_user_intrp3d(gph,p,"h",plvl,0.,False) plotvar1(:,:) = where(psfc.ge.plvl,plotvar1(:,:),plotvar1@_FillValue) plotvar1(:,:) = (1./9.806)*plotvar1(:,:) delete(gph) delete(p) delete(psfc) end if if ( var1.eq."theta" ) then var1title = "theta (K) at "+lvl1+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) theta = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) theta = wrf_user_getvar(ens,"theta",0) plotvar1(:,:) = wrf_user_intrp3d(theta,p,"h",plvl,0.,False) plotvar1(:,:) = where(psfc.ge.plvl,plotvar1(:,:),plotvar1@_FillValue) delete(theta) delete(p) delete(psfc) end if if ( var1.eq."vort" ) then var1title = "vort (s~S~-1~N~) at "+lvl1+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) vort = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) vort = wrf_user_getvar(ens,"avo",0) ; Absolute Vorticity (10-5 s-1) plotvar1(:,:) = wrf_user_intrp3d(vort,p,"h",plvl,0.,False) plotvar1(:,:) = where(psfc.ge.plvl,plotvar1(:,:),plotvar1@_FillValue) delete(vort) delete(p) delete(psfc) end if if ( var1.eq."thetae" ) then var1title = "epot temp (K) at "+lvl1+" mb" p850 = new( (/1,ny,nx/) , "float" , -999.0 ) p850 = 100.*plvl print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) tc = new( (/nlvls,ny,nx/) , "float" , -999.0 ) td = new( (/nlvls,ny,nx/) , "float" , -999.0 ) tc850 = new( (/1,ny,nx/) , "float" , -999.0 ) td850 = new( (/1,ny,nx/) , "float" , -999.0 ) mr = new( (/1,ny,nx/) , "float" , -999.0 ) thetae = new( (/1,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) tc = wrf_user_getvar(ens,"tc",0) td = wrf_user_getvar(ens,"td",0) p = 100.*wrf_user_getvar(ens,"pressure",0) tc850 = wrf_user_intrp3d(tc,p,"h",plvl,0.,False) td850 = wrf_user_intrp3d(td,p,"h",plvl,0.,False) mr(0,:,:) = mixhum_ptd( p850,td850+273.16,1 ) thetae = wrf_eth(mr,tc850+273.16,p850) plotvar1(:,:) = thetae(0,:,:) plotvar1(:,:) = where(psfc.ge.plvl,plotvar1(:,:),plotvar1@_FillValue) delete(psfc) delete(p) delete(tc) delete(td) delete(tc850) delete(td850) delete(mr) delete(thetae) delete(p850) end if ;Extract the winds for this pressure level print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" ) p = new( (/nlvls,ny,nx/) , "float" ) uvm = new( (/2,nlvls,ny,nx/) , "float" ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) ; Surface Pressure (Pa) p = wrf_user_getvar(ens,"pressure",0) ; Full model pressure (hPa) uvm = wrf_user_getvar(ens,"uvmet",0) ; earth-relative u and v winds uuplot(:,:) = wrf_user_intrp3d(uvm(0,:,:,:),p,"h",plvl,0.,False) vvplot(:,:) = wrf_user_intrp3d(uvm(1,:,:,:),p,"h",plvl,0.,False) uuplot(:,:) = where(psfc.ge.plvl,uuplot(:,:),uuplot@_FillValue) vvplot(:,:) = where(psfc.ge.plvl,vvplot(:,:),vvplot@_FillValue) uuplot(:,:) = uuplot(:,:)*1.94386 vvplot(:,:) = vvplot(:,:)*1.94386 delete(uvm) delete(p) delete(psfc) delete(plvl) end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Extract the second variable ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if ( lvl2.eq."2m" ) then ;Extract the 2 m variable if ( var2.eq."mrms_refl" ) then var2title = "observed reflectivity (dBZ)" print("Extracting observed MRMS reflectivity") mrms_filename = "/scratch/matt.flournoy/MRMS/smoothed_data/July06_0000z-0700z.nc" mrms_file = addfile(mrms_filename,"r") mrms_refl = mrms_file->refl(mrms_index(t),500:1500,2800:3500) ; contains refl data for the analysis time only (2D array) plotvar2_mrms = mrms_refl print("MRMS reflectivity read in") end if if ( var2.eq."mrms_rot" ) then var2title = "observed low-level rotation (s~S~-1~N~)" print("Extracting observed MRMS low-level rotation") mrms_filename = "/scratch/matt.flournoy/MRMS/lowlevelrotation/smoothed_data/20150706-" + hhmm_mrms(t) + ".nc" mrms_file = addfile(mrms_filename,"r") mrms_rot = mrms_file->rot(1000:2000,5000:7000) ; contains rotation data for the analysis time only (2D array) plotvar2_mrms = mrms_rot print("MRMS low-level rotation read in") end if if ( var2.eq."temp" ) then var2title = "temp (c) at 2 m, winds (kt) at 10 m" print("Extracting variables for ensemble member "+nens) plotvar2(:,:) = wrf_user_getvar(ens,"T2",0) plotvar2 = plotvar2 - 273.16 end if if ( var2.eq."dewp" ) then var2title = "dewp (c) at 2 m, winds (kt) at 10 m" print("Extracting variables for ensemble member "+nens) plotvar2(:,:) = wrf_user_getvar(ens,"td2",0) end if if ( var2.eq."rh" ) then var2title = "rh (%) at 2 m, winds (kt) at 10 m" print("Extracting variables for ensemble member "+nens) plotvar2(:,:) = wrf_user_getvar(ens,"rh2",0) end if if ( var2.eq."pmsl" ) then var2title = "pmsl (mb) , winds (kt) at 10 m" print("Extracting variables for ensemble member "+nens) plotvar2(:,:) = wrf_user_getvar(ens,"slp",0) end if if ( var2.eq."theta" ) then var2title = "pot temp (K) at lowest lvl, winds (kt) at 10 m" print("Extracting variables for ensemble member "+nens) theta = new( (/nlvls,ny,nx/) , "float" , -999.0 ) theta = wrf_user_getvar(ens,"theta",0) plotvar2(:,:) = theta(0,:,:) delete(theta) end if if ( var2.eq."cin" ) then var2title = "MLCIN (J/kg), winds (kt) at 10 m" print("Extracting variables for ensemble member "+nens) plotvar2_cape(:,:,:) = wrf_user_getvar(ens,"cape_2d",0) plotvar2(:,:) = plotvar2_cape(1,:,:) end if if ( var2.eq."max_uh" ) then var2title = "Max UH (m~S~2~N~ s~S~-2~N~), winds (kt) at 10 m" print("Exracting second variable for ensemble member "+nens) plotvar2(:,:) = ens->UP_HELI_MAX(0,:,:) end if if ( var2.eq."thetae" ) then var2title = "epot temp (K) at lowest lvl, winds (kt) at 10 m" print("Extracting variables for ensemble member "+nens) psfc = new( (/1,ny,nx/) , "float" , -999.0 ) t2 = new( (/1,ny,nx/) , "float" , -999.0 ) td2 = new( (/1,ny,nx/) , "float" , -999.0 ) mr = new( (/1,ny,nx/) , "float" , -999.0 ) thetae = new( (/1,ny,nx/) , "float" , -999.0 ) psfc(0,:,:) = wrf_user_getvar(ens,"PSFC",0) t2(0,:,:) = wrf_user_getvar(ens,"T2",0) td2(0,:,:) = wrf_user_getvar(ens,"td2",0) mr(0,:,:) = mixhum_ptd( psfc,td2+273.16,1 ) thetae = wrf_eth(mr,t2,psfc) plotvar2(:,:) = thetae(0,:,:) delete(psfc) delete(t2) delete(td2) delete(mr) delete(thetae) end if ; if ( var2.eq."rot" ) then ; var2title = "Sim. Low-Level Rotation (10~S~-5~N~ s~S~-1~N~)" ; plotvar2_all_abs_vort = wrf_user_getvar(ens,"avo",0) ; Absolute Vorticity (10-5 s-1) ... all levels, lats, lons for that time ; plotvar2_all_avo = plotvar2_all_abs_vort(0,:,:) ; Absolute Vorticity at the lowest model level, all lats/lons for that time ; plotvar2(:,:) = dim_avg_n(plotvar2_all_avo,0) ; Absolute Vorticity average in the first 10 levels at all lats/lons for that time ; plotvar2(:,:) = plotvar2_all_avo - 2 * 0.000072921159 * sin(43.436648*pi/180) ; convert to relative vorticity ; plotvar2(:,:) = mask(plotvar2,(plotvar2.ge.50),True) ; <-- masks all absolute vorticity below 50*10e-5 for plotting purposes ; end if if ( var2.eq."rot" ) then var2title = "Sim. 1 km AGL Rotation (10~S~-5~N~ s~S~-1~N~)" z = wrf_user_getvar(ens,"z",0) ; <-- extract z to interpolate relative vorticity to 1 km AGL plotvar2_all_abs_vort = wrf_user_getvar(ens,"REL_VORT",0) ; Relative Vorticity (10-5 s-1) ... all levels, lats, lons for that time plotvar2(:,:) = wrf_user_intrp3d(plotvar2_all_abs_vort,z,"h",1000,0.,False) ; Relative Vorticity at 1 km AGL, all lats/lons for that time plotvar2 = plotvar2 * 100000 ; plotvar2(:,:) = dim_avg_n(plotvar2_all_avo,0) ; Absolute Vorticity average in the first 10 levels at all lats/lons for that time ; plotvar2(:,:) = mask(plotvar2,(plotvar2.ge.50),True) ; <-- masks all absolute vorticity below 50*10e-5 for plotting purposes end if else ;Extract the variable for this pressure level plvl = tofloat(lvl2) if ( var2.eq."temp" ) then var2title = "temp (C) and winds (kt) at "+lvl2+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) tc = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) tc = wrf_user_getvar(ens,"tc",0) plotvar2(:,:) = wrf_user_intrp3d(tc,p,"h",plvl,0.,False) plotvar2(:,:) = where(psfc.ge.plvl,plotvar2(:,:),plotvar2@_FillValue) delete(tc) delete(p) delete(psfc) end if if ( var2.eq."dewp" ) then var2title = "dewp (C) and winds (kt) at "+lvl2+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) td = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) td = wrf_user_getvar(ens,"td",0) plotvar2(:,:) = wrf_user_intrp3d(td,p,"h",plvl,0.,False) plotvar2(:,:) = where(psfc.ge.plvl,plotvar2(:,:),plotvar2@_FillValue) delete(td) delete(p) delete(psfc) end if if ( var2.eq."rh" ) then var2title = "rh (C) and winds (kt) at "+lvl2+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) rh = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) rh = wrf_user_getvar(ens,"rh",0) plotvar2(:,:) = wrf_user_intrp3d(rh,p,"h",plvl,0.,False) plotvar2(:,:) = where(psfc.ge.plvl,plotvar2(:,:),plotvar2@_FillValue) delete(rh) delete(p) delete(psfc) end if if ( var2.eq."gph" ) then var2title = "geo height (m) and winds (kt) at "+lvl2+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) gph = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) gph = wrf_user_getvar(ens,"geopt",0) ; Full model geopotential (m2 s-2) plotvar2(:,:) = wrf_user_intrp3d(gph,p,"h",plvl,0.,False) plotvar2(:,:) = where(psfc.ge.plvl,plotvar2(:,:),plotvar2@_FillValue) plotvar2(:,:) = (1./9.806)*plotvar2(:,:) delete(gph) delete(p) delete(psfc) end if if ( var2.eq."theta" ) then var2title = "theta (K) and winds (kt) at "+lvl2+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) theta = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) theta = wrf_user_getvar(ens,"theta",0) plotvar2(:,:) = wrf_user_intrp3d(theta,p,"h",plvl,0.,False) plotvar2(:,:) = where(psfc.ge.plvl,plotvar2(:,:),plotvar2@_FillValue) delete(theta) delete(p) delete(psfc) end if if ( var2.eq."vort" ) then var2title = "vort (s~S~-1~N~) and winds (kt) at "+lvl2+" mb" print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) vort = new( (/nlvls,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) p = wrf_user_getvar(ens,"pressure",0) vort = wrf_user_getvar(ens,"avo",0) ; Absolute Vorticity (10-5 s-1) plotvar2(:,:) = wrf_user_intrp3d(vort,p,"h",plvl,0.,False) plotvar2(:,:) = where(psfc.ge.plvl,plotvar2(:,:),plotvar2@_FillValue) delete(vort) delete(p) delete(psfc) end if if ( var2.eq."thetae" ) then var2title = "epot temp (K) and winds (kt) at "+lvl2+" mb" p850 = new( (/1,ny,nx/) , "float" , -999.0 ) p850 = 100.*plvl print("Extracting variables for ensemble member "+nens) psfc = new( (/ny,nx/) , "float" , -999.0 ) p = new( (/nlvls,ny,nx/) , "float" , -999.0 ) tc = new( (/nlvls,ny,nx/) , "float" , -999.0 ) td = new( (/nlvls,ny,nx/) , "float" , -999.0 ) tc850 = new( (/1,ny,nx/) , "float" , -999.0 ) td850 = new( (/1,ny,nx/) , "float" , -999.0 ) mr = new( (/1,ny,nx/) , "float" , -999.0 ) thetae = new( (/1,ny,nx/) , "float" , -999.0 ) psfc = 0.01*wrf_user_getvar(ens,"PSFC",0) tc = wrf_user_getvar(ens,"tc",0) td = wrf_user_getvar(ens,"td",0) p = 100.*wrf_user_getvar(ens,"pressure",0) tc850 = wrf_user_intrp3d(tc,p,"h",plvl,0.,False) td850 = wrf_user_intrp3d(td,p,"h",plvl,0.,False) mr(0,:,:) = mixhum_ptd( p850,td850+273.16,1 ) thetae = wrf_eth(mr,tc850+273.16,p850) plotvar2(:,:) = thetae(0,:,:) plotvar2(:,:) = where(psfc.ge.plvl,plotvar2(:,:),plotvar2@_FillValue) delete(psfc) delete(p) delete(tc) delete(td) delete(tc850) delete(td850) delete(mr) delete(thetae) delete(p850) end if end if var1mem = w_plot if (var2 .eq. "mrms_refl" .or. var2 .eq. "mrms_rot") then var2mem = plotvar2_mrms else var2mem = plotvar2 end if uumem = uuplot vvmem = vvplot ; uumem = u_horiz_vort ; vvmem = v_horiz_vort ;;;;;;;;;;;;;;;;;;;;;;;; ; Define color maps ;;;;;;;;;;;;;;;;;;;;;;;; cmap_temp = "cmp_b2r" cmap_dewp = "WhiteGreen" cmap_pmsl = "BlueYellowRed" cmap_rh = "WhiteGreen" cmap_gph = "WhiteBlueGreenYellowRed" cmap_theta = "cmp_b2r" cmap_thetae = "BlGrYeOrReVi200" cmap_vort = "BlWhRe" cmap_refl = "radar" cmap_cape = "BlGrYeOrReVi200" cmap_maxuh = "BlGrYeOrReVi200" cmap_rot = "MPL_OrRd" cmap_rot_tilt = "MPL_RdBu" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Define plot variables that are the same for all plots ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; plotformat = "png" res = True res@cnFillMode = "AreaFill" res@cnMonoFillPattern = True res@cnInfoLabelOn = False ; Contour info label. res@cnLinesOn = False res@cnLineLabelsOn = False res@cnLineLabelInterval = 2 res@cnLineLabelFontHeightF = 0.007 res@cnLineLabelDensityF = 0.8 ; Controls number of line labels on the plot res@cnLineLabelFont = 22 res@cnLineLabelAngleF = -1.0 res@cnLineLabelBackgroundColor = 0 res@cnLineLabelPerimSpaceF = 0.20 res@gsnAddCyclic = False res@gsnDraw = False res@gsnFrame = False res@gsnShape = True res@gsnScale = False res@gsnSpreadColors = True res@lbLabelStride = 2 res@lbLabelAutoStride = True res@lbAutoManage = False res@pmLabelBarSide = "Bottom" res@lbJustification = "TopLeft" res@lbOrientation = "Horizontal" res@pmLabelBarWidthF = 0.40 ; control size of colorbar res@pmLabelBarHeightF = 0.06 res@pmLabelBarOrthogonalPosF = 0.000 ; position wrt plot res@lbLabelFontHeightF = 0.010 ; label bar font size res@lbLabelFont = 22 ; label bar font res@lbPerimOn = False ; no box around label bar res@lbLabelAutoStride = True ; Avoid overlapping labels res@lbLabelOffsetF = 0.05 res@lbLabelPosition = "Top" res@lbTitleFontHeightF = 0.012 ; label bar title font size res@lbTitleFont = 22 ; label bar title font res@lbTitlePosition = "Left" res@lbTitleDirection = "Across" res@lbTitleOffsetF = 0.00 res@lbTitleExtentF = 0.15 res@mpProjection = "Mercator" res@mpLimitMode = "Corners" lat_ll = clat - 0.5*plot_height lat_ur = clat + 0.5*plot_height lon_ll = clon - 0.5*plot_width lon_ur = clon + 0.5*plot_width res@mpLeftCornerLatF = lat_ll res@mpLeftCornerLonF = lon_ll res@mpRightCornerLatF = lat_ur res@mpRightCornerLonF = lon_ur res@mpShapeMode = "FreeAspect" res@mpOutlineOn = True res@mpGridAndLimbOn = False res@mpPerimOn = True res@mpFillOn = False res@mpPerimDrawOrder = "PostDraw" res@mpGeophysicalLineColor = "Black" res@mpNationalLineColor = "Black" res@mpUSStateLineColor = "Black" res@mpGridLineColor = "Black" res@mpLimbLineColor = "Black" res@mpPerimLineColor = "Black" res@mpGeophysicalLineThicknessF = 1.0 res@mpGridLineThicknessF = 1.0 res@mpLimbLineThicknessF = 1.0 res@mpNationalLineThicknessF = 0.5 res@mpUSStateLineThicknessF = 1.0 res@mpDataResolution = "FinestResolution" res@mpOutlineBoundarySets = "AllBoundaries" res@mpUSStateLineThicknessF = 4.0 res@mpDataBaseVersion = "Ncarg4_1" res@mpDataSetName = "Earth..2" res@mpCountyLineThicknessF = 0.5 if ( var1.eq."refl" .and. var2.eq."pmsl" ) then res@cnFillDrawOrder = "PostDraw" ; draw reflectivity overtop of everything else if pressure is the second variable else res@cnFillDrawOrder = "PreDraw" ; draw var1 beneath everything else end if res@tiMainOn = True res@tiMainFontHeightF = 0.025 res@tiMainJust = "TopLeft" res@tiMainPosition = "Left" res@tiMainFontQuality = "High" res@tiMainOffsetYF = 0.01 res@gsnLeftString = "" res@gsnRightString = "" res@vpXF = 0.05 res@vpYF = 0.85 res@vpWidthF = 0.85 res@vpHeightF = 0.70 res@tmXBLabelsOn = False txres = True txres@txBackgroundFillColor = "White" txres@txFont = 22 txres@txPerimSpaceF = 0.3 txres@txPerimThicknessF = 2.0 txres@txPerimColor = "Black" ; polyres = True ; polyres@gsMarkerIndex = 5 ; polyres@gsMarkerSizeF = 0.012 ; polyres@gsMarkerThicknessF = 3.0 ; polyres@gsMarkerColor = (/"Black"/) ; xpt = clon ; ypt = clat ;;;;;;;;;;;;;;;;;;;;;; ; Copy res for res2 ;;;;;;;;;;;;;;;;;;;;;; res2 = True res2@cnFillMode = "AreaFill" res2@cnMonoFillPattern = True res2@cnInfoLabelOn = False ; Contour info label. res2@cnLinesOn = False res2@cnLineLabelsOn = True res2@cnLineLabelInterval = 2 res2@cnLineLabelFontHeightF = 0.007 res2@cnLineLabelDensityF = 0.8 ; Controls number of line labels on the plot res2@cnLineLabelFont = 22 res2@cnLineLabelAngleF = -1.0 res2@cnLineLabelBackgroundColor = 0 res2@cnLineLabelPerimSpaceF = 0.20 res2@gsnAddCyclic = False res2@gsnDraw = False res2@gsnFrame = False res2@gsnShape = True res2@gsnScale = False res2@gsnSpreadColors = True res2@lbLabelStride = 2 res2@lbLabelAutoStride = True res2@lbAutoManage = False res2@pmLabelBarSide = "Bottom" res2@lbJustification = "TopLeft" res2@lbOrientation = "Horizontal" res2@pmLabelBarWidthF = 0.40 ; control size of colorbar res2@pmLabelBarHeightF = 0.06 res2@pmLabelBarOrthogonalPosF = 0.000 ; position wrt plot res2@lbLabelFontHeightF = 0.010 ; label bar font size res2@lbLabelFont = 22 ; label bar font res2@lbPerimOn = False ; no box around label bar res2@lbLabelAutoStride = True ; Avoid overlapping labels res2@lbLabelOffsetF = 0.05 res2@lbLabelPosition = "Top" res2@lbTitleFontHeightF = 0.012 ; label bar title font size res2@lbTitleFont = 22 ; label bar title font res2@lbTitlePosition = "Left" res2@lbTitleDirection = "Across" res2@lbTitleOffsetF = 0.00 res2@lbTitleExtentF = 0.15 res2@mpProjection = "Mercator" res2@mpLimitMode = "Corners" lat_ll = clat - 0.5*plot_height lat_ur = clat + 0.5*plot_height lon_ll = clon - 0.5*plot_width lon_ur = clon + 0.5*plot_width res2@mpLeftCornerLatF = lat_ll res2@mpLeftCornerLonF = lon_ll res2@mpRightCornerLatF = lat_ur res2@mpRightCornerLonF = lon_ur res2@mpShapeMode = "FreeAspect" res2@mpOutlineOn = True res2@mpGridAndLimbOn = False res2@mpPerimOn = True res2@mpFillOn = False res2@mpPerimDrawOrder = "PreDraw" res2@mpOutlineBoundarySets = "AllBoundaries" res2@mpUSStateLineThicknessF = 4.0 res2@tiMainOn = True res2@tiMainFontHeightF = 0.014 res2@tiMainJust = "TopLeft" res2@tiMainPosition = "Left" res2@tiMainFontQuality = "High" res2@tiMainOffsetYF = 0.01 res2@vpXF = 0.05 res2@vpYF = 0.85 res2@vpWidthF = 0.85 res2@vpHeightF = 0.70 res2@tmXBLabelsOn = False ;;;;;;;;;;;;;;;;;;;;;; ; Copy res for res3 ;;;;;;;;;;;;;;;;;;;;;; res3 = True res3@cnFillMode = "AreaFill" res3@cnMonoFillPattern = True res3@cnInfoLabelOn = False ; Contour info label. res3@cnLinesOn = False res3@cnLineLabelsOn = True res3@cnLineLabelInterval = 2 res3@cnLineLabelFontHeightF = 0.007 res3@cnLineLabelDensityF = 0.8 ; Controls number of line labels on the plot res3@cnLineLabelFont = 22 res3@cnLineLabelAngleF = -1.0 res3@cnLineLabelBackgroundColor = 0 res3@cnLineLabelPerimSpaceF = 0.20 res3@gsnAddCyclic = False res3@gsnDraw = False res3@gsnFrame = False res3@gsnShape = True res3@gsnScale = False res3@gsnSpreadColors = True res3@lbLabelStride = 2 res3@lbLabelAutoStride = True res3@lbAutoManage = False res3@pmLabelBarSide = "Bottom" res3@lbJustification = "TopLeft" res3@lbOrientation = "Horizontal" res3@pmLabelBarWidthF = 0.40 ; control size of colorbar res3@pmLabelBarHeightF = 0.06 res3@pmLabelBarOrthogonalPosF = 0.000 ; position wrt plot res3@lbLabelFontHeightF = 0.010 ; label bar font size res3@lbLabelFont = 22 ; label bar font res3@lbPerimOn = False ; no box around label bar res3@lbLabelAutoStride = True ; Avoid overlapping labels res3@lbLabelOffsetF = 0.05 res3@lbLabelPosition = "Top" res3@lbTitleFontHeightF = 0.012 ; label bar title font size res3@lbTitleFont = 22 ; label bar title font res3@lbTitlePosition = "Left" res3@lbTitleDirection = "Across" res3@lbTitleOffsetF = 0.00 res3@lbTitleExtentF = 0.15 res3@mpProjection = "Mercator" res3@mpLimitMode = "Corners" lat_ll = clat - 0.5*plot_height lat_ur = clat + 0.5*plot_height lon_ll = clon - 0.5*plot_width lon_ur = clon + 0.5*plot_width res3@mpLeftCornerLatF = lat_ll res3@mpLeftCornerLonF = lon_ll res3@mpRightCornerLatF = lat_ur res3@mpRightCornerLonF = lon_ur res3@mpShapeMode = "FreeAspect" res3@mpOutlineOn = True res3@mpGridAndLimbOn = False res3@mpPerimOn = True res3@mpFillOn = False res3@mpPerimDrawOrder = "PreDraw" res3@mpOutlineBoundarySets = "AllBoundaries" res3@mpUSStateLineThicknessF = 4.0 res3@tiMainOn = True res3@tiMainFontHeightF = 0.014 res3@tiMainJust = "TopLeft" res3@tiMainPosition = "Left" res3@tiMainFontQuality = "High" res3@tiMainOffsetYF = 0.01 res3@vpXF = 0.05 res3@vpYF = 0.85 res3@vpWidthF = 0.85 res3@vpHeightF = 0.70 res3@tmXBLabelsOn = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;******************************; ; Loop through the plots to make 'em <-- ; MAKE PLOT DIRECTORY HERE ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;******************************; print("VAR1: Min = "+min(var1mem)+" Max = "+max(var1mem)) printVarSummary(var1mem) resolution = "1km" plotdir = "/work/matt.flournoy/plots/FINAL_0300UTC_333m/mem30/w_levels/mvS/" ; <-- used this line for 3 km plots so they are grouped by variable, not by time, ; so they can be viewed in time-loops by ensemble or all members at a certain time. system("mkdir -p "+plotdir) fillvalue = -999.0 var2mem!0 = "lat" var2mem!1 = "lon" var1mem@lat2d = latgrid var1mem@lon2d = longrid var2mem@lat2d = latgrid var2mem@lon2d = longrid ; var2mem&lat = lats_mrms_rot(1000:2000) ; var2mem&lon = lons_mrms_rot(5000:7000) uumem@lat2d = latgrid uumem@lon2d = longrid vvmem@lat2d = latgrid vvmem@lon2d = longrid ; var2mem&lat@units = "degrees north" ; var2mem&lon@units = "degrees east" ; defaults...change below if needed spreadstart = 2 spreadend = -1 reverse_color_map = "False" fillpattern = 0 linethickness1 = 1.0 linethickness2 = 3.0 linethickness3 = 3.0 linethickness4 = 8.0 linedash1 = 0 linedash2 = 0 linedash3 = 0 titlestring = "Updraft/Downdraft ~C~~Z70~Member " + sprinti("%2.0i",nens) + " - " + hh + mn + " UTC" ; titlestring = "300m AGL Non-linear Dynamic p' ~C~~Z70~Member " + sprinti("%2.0i",nens) + " - " + hh + mn + " UTC" if ( var1.eq."temp" ) then colors1 = cmap_temp linecolor1 = "Red" units1 = "C" cnspacing1 = 1.0 if ( lvl1.eq."2m" ) then cnmin1 = 16.0 cnmax1 = 40.0 else plvl = tofloat(lvl1) if ( plvl.ge.900. ) then cnmin1 = 12.0 cnmax1 = 36.0 end if if ( plvl.ge.800. .and. plvl.lt.900. ) then cnmin1 = 6.0 cnmax1 = 30.0 end if if ( plvl.ge.600. .and. plvl.lt.800. ) then cnmin1 = 0.0 cnmax1 = 16.0 end if if ( plvl.ge.400. .and. plvl.lt.600. ) then cnmin1 = -20.0 cnmax1 = -2.0 end if if ( plvl.lt.400. ) then cnmin1 = -60.0 cnmax1 = -40.0 end if end if end if if ( var1.eq."dewp" ) then colors1 = cmap_dewp linecolor1 = "DarkGreen" units1 = "C" cnspacing1 = 1.0 if ( lvl1.eq."2m" ) then cnmin1 = 10.0 cnmax1 = 24.0 else plvl = tofloat(lvl1) if ( plvl.ge.800. ) then cnmin1 = 0.0 cnmax1 = 20.0 end if if ( plvl.ge.600. .and. plvl.lt.800. ) then cnmin1 = -10.0 cnmax1 = 10.0 end if if ( plvl.ge.400. .and. plvl.lt.600. ) then cnmin1 = -40.0 cnmax1 = -10.0 end if if ( plvl.lt.400. ) then cnmin1 = -70.0 cnmax1 = -40.0 end if end if end if if ( var1.eq."pmsl" ) then colors1 = cmap_pmsl reverse_color_map = "True" linecolor1 = "Black" units1 = "%" cnspacing1 = 1.0 cnmin1 = 1000.0 cnmax1 = 1020.0 end if if ( var1.eq."rh" ) then colors1 = cmap_rh linecolor1 = "DarkGreen" units1 = "%" cnspacing1 = 2.5 cnmin1 = 30.0 cnmax1 = 95.0 end if if ( var1.eq."cape" ) then colors1 = cmap_cape linecolor1 = "Black" units1 = "J/kg" cnspacing1 = 200 cnmin1 = 0 cnmax1 = 4000 end if if ( var1.eq."refl" ) then colors1 = cmap_refl linecolor1 = "Black" units1 = "dBZ" cnspacing1 = 5 cnmin1 = 20 cnmax1 = 75 end if if ( var1.eq."gph" ) then colors1 = cmap_gph reverse_color_map = "True" linecolor1 = "Black" units1 = "m" plvl = tofloat(lvl1) if ( plvl.eq.850. ) then cnspacing1 = 10.0 cnmin1 = 1400.0 cnmax1 = 1600.0 end if if ( plvl.eq.700. ) then cnspacing1 = 10.0 cnmin1 = 3000.0 cnmax1 = 3200.0 end if if ( plvl.eq.500. ) then cnspacing1 = 20.0 cnmin1 = 5580.0 cnmax1 = 5940.0 end if if ( plvl.eq.250. ) then cnspacing1 = 40.0 cnmin1 = 10400.0 cnmax1 = 11000.0 end if end if if ( var1.eq."theta" ) then colors1 = cmap_theta linecolor1 = "Red" units1 = "C" cnspacing1 = 1.0 cnmin1 = 290.0 cnmax1 = 320.0 end if if ( var1.eq."vort" ) then colors1 = cmap_vort linecolor1 = "Black" units1 = "10~S~-5~N~ s~S~-1~N~" cnspacing1 = 10 cnmin1 = -70 cnmax1 = 70 end if if ( var1.eq."thetae" ) then colors1 = cmap_thetae linecolor1 = "DarkGreen" units1 = "K" cnspacing1 = 2.0 cnmin1 = 300.0 cnmax1 = 370.0 end if if ( var1.eq."rot" ) then colors1 = cmap_rot linecolor1 = "Black" units1 = "10~S~5~N~ s~S~-1~N~" cnspacing1 = 50.0 cnmin1 = 0 cnmax1 = 1500.0 end if if ( var1.eq."helicity" ) then colors1 = cmap_rot linecolor1 = "Black" units1 = "m~S~2~N~ s~S~-2~N~" cnspacing1 = 50.0 cnmin1 = 0 cnmax1 = 1000.00 end if if ( var1.eq."rot_horiz" ) then ; <-- use these for horizontal vorticity tilting colors1 = cmap_rot_tilt linecolor1 = "Black" units1 = "10~S~-5~N~ s~S~-2~N~" cnspacing1 = 0.5 cnmin1 = -10.0 cnmax1 = 10.0 end if if ( var1.eq."p_dyn_NL") then colors1 = cmap_rot_tilt linecolor1 = "Black" units1 = "hPa" cnspacing1 = 0.10 cnmin1 = -2.0 cnmax1 = 2.0 end if if ( var1.eq."w_levels") then colors1 = cmap_rot_tilt units1 = "m s~S~-1~N~" end if if (isatt(var1mem,"long_name")) then delete(var1mem@long_name) end if if (isatt(var1mem,"units")) then delete(var1mem@units) end if if (isatt(var1mem,"Units")) then delete(var1mem@Units) end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Copy the above crap for var2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if ( var2.eq."temp" ) then colors2 = cmap_temp linecolor2 = "Red" units2 = "C" cnspacing2 = 1.0 if ( lvl2.eq."2m" ) then cnmin2 = 16.0 cnmax2 = 40.0 else plvl = tofloat(lvl2) if ( plvl.ge.900. ) then cnmin2 = 12.0 cnmax2 = 36.0 end if if ( plvl.ge.800. .and. plvl.lt.900. ) then cnmin2 = 6.0 cnmax2 = 30.0 end if if ( plvl.ge.600. .and. plvl.lt.800. ) then cnmin2 = 0.0 cnmax2 = 16.0 end if if ( plvl.ge.400. .and. plvl.lt.600. ) then cnmin2 = -20.0 cnmax2 = -2.0 end if if ( plvl.lt.400. ) then cnmin2 = -60.0 cnmax2 = -40.0 end if end if end if if ( var2.eq."dewp" ) then colors2 = cmap_dewp linecolor2 = "DarkGreen" units2 = "C" cnspacing2 = 1.0 if ( lvl2.eq."2m" ) then cnmin2 = 10.0 cnmax2 = 24.0 else plvl = tofloat(lvl2) if ( plvl.ge.800. ) then cnmin2 = 0.0 cnmax2 = 20.0 end if if ( plvl.ge.600. .and. plvl.lt.800. ) then cnmin2 = -10.0 cnmax2 = 10.0 end if if ( plvl.ge.400. .and. plvl.lt.600. ) then cnmin2 = -40.0 cnmax2 = -10.0 end if if ( plvl.lt.400. ) then cnmin2 = -70.0 cnmax2 = -40.0 end if end if end if if ( var2.eq."pmsl" ) then colors2 = cmap_pmsl linecolor2 = "Black" units2 = "%" cnspacing2 = 1.0 cnmin2 = 1000.0 cnmax2 = 1020.0 end if if ( var2.eq."cin" ) then colors2 = cmap_cape linecolor2 = "Black" units2 = "- J/kg" cnspacing2 = 50 cnmin2 = 0 cnmax2 = 500 end if if ( var2.eq."rh" ) then colors2 = cmap_rh linecolor2 = "DarkGreen" units2 = "%" cnspacing2 = 2.5 cnmin2 = 30.0 cnmax2 = 95.0 end if if ( var2.eq."gph" ) then colors2 = cmap_gph linecolor2 = "Black" units2 = "m" plvl = tofloat(lvl2) if ( plvl.eq.850. ) then cnspacing2 = 10.0 cnmin2 = 1400.0 cnmax2 = 1600.0 end if if ( plvl.eq.700. ) then cnspacing2 = 10.0 cnmin2 = 3000.0 cnmax2 = 3200.0 end if if ( plvl.eq.500. ) then cnspacing2 = 20.0 cnmin2 = 5580.0 cnmax2 = 5940.0 end if if ( plvl.eq.250. ) then cnspacing2 = 40.0 cnmin2 = 10400.0 cnmax2 = 11000.0 end if end if if ( var2.eq."theta" ) then colors2 = cmap_theta linecolor2 = "Red" units2 = "C" cnspacing2 = 1.0 cnmin2 = 290.0 cnmax2 = 320.0 end if if ( var2.eq."vort" ) then colors2 = cmap_vort linecolor2 = "Black" units2 = "10~S~-5~N~ s~S~-1~N~" cnspacing2 = 10 cnmin2 = -70 cnmax2 = 70 end if if ( var2.eq."max_uh" ) then colors2 = cmap_maxuh linecolor2 = "Black" units2 = "m~S~2~N~ s~S~-2~N~" cnspacing2 = 50 cnmin2 = 50 cnmax2 = 400 end if if ( var2.eq."thetae" ) then colors2 = cmap_thetae linecolor2 = "DarkGreen" units2 = "K" cnspacing2 = 2.0 cnmin2 = 300.0 cnmax2 = 370.0 end if if ( var2 .eq. "mrms_refl" ) then ; <-- draw 40 dBZ observed MRMS reflectivity contour colors2 = cmap_refl linecolor2 = "Black" units2 = "dBZ" cnspacing2 = 1.0 cnmin2 = 39.0 cnmax2 = 41.0 end if if ( var2.eq."mrms_rot" ) then ; <-- might need to change this, depending on the values of mrms_rot colors2 = cmap_rot linecolor2 = "Black" units2 = "10~S~5~N~ s~S~-1~N~" ; <-- the data is actually just s^-1, but just label the contours as 1000 and 2000 cnspacing2 = 0.01 cnmin2 = 0.0 cnmax2 = 0.02 end if if ( var2.eq."rot" ) then colors2 = cmap_rot linecolor2 = "maroon1" units2 = "10~S~5~N~ s~S~-1~N~" cnspacing2 = 600.0 cnmin2 = 600 cnmax2 = 4200.0 end if if (isatt(var2mem,"long_name")) then delete(var2mem@long_name) end if if (isatt(var2mem,"units")) then delete(var2mem@units) end if if (isatt(var2mem,"Units")) then delete(var2mem@Units) end if if (isatt(uumem,"long_name")) then delete(uumem@long_name) end if if (isatt(uumem,"units")) then delete(uumem@units) end if if (isatt(uumem,"Units")) then delete(uumem@Units) end if if (isatt(vvmem,"long_name")) then delete(vvmem@long_name) end if if (isatt(vvmem,"units")) then delete(vvmem@units) end if if (isatt(vvmem,"Units")) then delete(vvmem@Units) end if plotfile = plotdir+yymmdd+"_"+hh+mn+"_F"+ff+"_"+expm+"_mem"+nens+"_"+var1+"_"+var2+"_"+lvl1 ; print(plotfile) xwks = gsn_open_wks(plotformat,plotfile) res@cnLevelSelectionMode = "ExplicitLevels" ; Plot the first variable (with shading and wind barbs) ; line colors are next, with each pair having the updraft color and then the downdraft color linecolors = (/"darkgreen","darkorange4","chartreuse4","darkgoldenrod4","chartreuse3","darkgoldenrod3","chartreuse2","darkgoldenrod2","chartreuse1","darkgoldenrod1"/) dashpatterns = (/2,0,2,0,2,0,2,0,2,0/) contour_levels = (/5,-5,5,-5,5,-5,5,-5,5,-5/) var1_index = (/0,0,1,1,2,2,3,3,4,4/) gsn_define_colormap(xwks,colors1) if ( reverse_color_map.eq."True" ) then gsn_reverse_colormap(xwks) end if do i = 0 , (num_levels * 2) - 1 res@cnFillOn = True ; <-- changed this to true to turn on shading for var1 res@cnLinesOn = False ; <-- changed this to false to turn off var1 contours res@cnFillPattern = fillpattern res@cnLineColor = linecolors(i) res@cnLineLabelFontColor = linecolors(i) res@cnLineDashPattern = dashpatterns(i) res@cnLineThicknessF = linethickness1 res@cnLevels = contour_levels ; <-- set +/- 5 m/s as the updraft contour at each level ; res@cnMinLevelValF = cnmin1 ; res@cnMaxLevelValF = cnmax1 ; res@cnLevelSpacingF = cnspacing1 res@lbLabelBarOn = True res@lbTitleString = units1 res@mpUSStateLineColor = "Black" res@pmLabelBarParallelPosF = 0.22 res@tiMainString = titlestring res@vfXCStride = 2 res@vfYCStride = 2 res@vcGlyphStyle = "WindBarb" res@vcWindBarbLineThicknessF = 3.0 res@vcWindBarbColor = "Black" res@vcRefLengthF = 0.020 res@vcRefAnnoOn = False res@vcZeroFLabelOn = False res@vcMinDistanceF = 0.025 ; <-- thins the wind barbs in NDC space res@gsnScalarContour = True res@gsnSpreadColorStart = spreadstart res@gsnSpreadColorEnd = spreadend res@gsnLeftString = "" res@gsnRightString = "" res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False if (i.eq.0) then plot1 = gsn_csm_vector_scalar_map(xwks,uumem,vvmem,var1mem(var1_index(i),:,:),res) else plot2 = gsn_csm_contour(xwks,var1mem(var1_index(i),:,:),res) overlay(plot1,plot2) end if end do ; draw(plot1) ; frame(xwks) ; delete(xwks) <-- commented this out because it caused a segementation fault ... code still seems to work fine delete(colors1) ; delete(plot1) <-- commented this out because it caused a segementation fault ... code still seems to work fine delete(var1mem) ; delete(uumem) ; delete(vvmem) ; Plot the second variable (just contours) res2@cnLevelSelectionMode = "ManualLevels" res2@cnFillOn = False ; <-- changed this to false to turn off shading for var2 res2@cnLinesOn = True ; <-- changed this to True to turn on var2 contours res2@cnFillPattern = fillpattern res2@cnLineColor = linecolor2 res2@cnLineLabelFontColor = linecolor2 res2@cnLineDashPattern = linedash2 res2@cnLineThicknessF = linethickness2 res2@cnLineLabelsOn = False res2@cnMinLevelValF = cnmin2 res2@cnMaxLevelValF = cnmax2 res2@cnLevelSpacingF = cnspacing2 res2@lbLabelBarOn = True res2@lbTitleString = units1 res2@mpUSStateLineColor = "Black" res2@pmLabelBarParallelPosF = 0.22 res2@tiMainString = titlestring res2@vfXCStride = 2 res2@vfYCStride = 2 res2@vcGlyphStyle = "WindBarb" res2@vcWindBarbLineThicknessF = 0.0 ; <-- set this to 0 for var2 so it doesn't plot barbs again (even though it shouldn't since I turned cnFillOn = False) res2@vcWindBarbColor = "Black" res2@vcRefLengthF = 0.020 res2@vcRefAnnoOn = False res2@vcZeroFLabelOn = False res2@vcMinDistanceF = 0.04 ; <-- thins the wind barbs in NDC space (which aren't there for var2 anyway) res2@gsnScalarContour = True res2@gsnSpreadColorStart = spreadstart res2@gsnSpreadColorEnd = spreadend res2@gsnLeftString = "" res2@gsnRightString = "" res2@gsnDraw = False res2@gsnFrame = False res2@cnInfoLabelOn = False plot2 = gsn_csm_contour(xwks,var2mem,res2) ;----------------------------------------------------------------------; ; Insert tornado location with a marker (in NDC coordinates). ;----------------------------------------------------------------------; shpres = True shpres@gsMarkerIndex = 4 ; circle shpres@gsMarkerColor = "Black" shpres@gsMarkerThicknessF = 2.0 shpres@gsMarkerSizeF = 12.0 ; gsn_polymarker_ndc(xwks,0.620,0.49,shpres) overlay(plot1,plot2) ; drawNDCGrid(xwks) ; draw(plot1) ; frame(xwks) ; delete(xwks) <-- commented this out because it caused a segementation fault ... code still seems to work fine ; delete(colors2) ; delete(plot2) <-- commented this out because it caused a segementation fault ... code still seems to work fine delete(var2mem) ; Plot the third variable (just 40 dBZ contour) var3 = "refl" if ( var3.eq."refl" ) then var3title = "Sim. Low-level Reflectivity (dBZ)" print("Extracting variables for ensemble member "+nens) plotvar3_all = wrf_user_getvar(ens,"dbz",0) plotvar3 = plotvar3_all(0,:,:) ; <-- extract lowest model level reflectivity ; plotvar3 = mask(plotvar1_all,(plotvar1_all.ge.20),True) ; <-- masks all reflectivity below 20 dBZ for plotting purposes end if var3mem = plotvar3 var3mem!0 = "lat" var3mem!1 = "lon" var3mem@lat2d = latgrid var3mem@lon2d = longrid if ( var3 .eq. "refl" ) then ; <-- draw 40 dBZ simulated reflectivity contour colors3 = cmap_refl linecolor3 = "gray40" units3 = "dBZ" cnspacing3 = 0.001 cnmin3 = 39.999 cnmax3 = 40.001 end if res3@cnLevelSelectionMode = "ManualLevels" res3@cnFillOn = False ; <-- changed this to false to turn off shading for var2 res3@cnLinesOn = True ; <-- changed this to True to turn on var2 contours res3@cnFillPattern = fillpattern res3@cnLineColor = linecolor3 res3@cnLineLabelFontColor = linecolor3 res3@cnLineDashPattern = linedash3 res3@cnLineThicknessF = linethickness3 res3@cnLineLabelsOn = False res3@cnMinLevelValF = cnmin3 res3@cnMaxLevelValF = cnmax3 res3@cnLevelSpacingF = cnspacing3 res3@lbLabelBarOn = True res3@lbTitleString = units1 res3@mpUSStateLineColor = "Black" res3@pmLabelBarParallelPosF = 0.22 res3@tiMainString = titlestring res3@vfXCStride = 2 res3@vfYCStride = 2 res3@vcGlyphStyle = "WindBarb" res3@vcWindBarbLineThicknessF = 0.0 ; <-- set this to 0 for var2 so it doesn't plot barbs again (even though it shouldn't since I turned cnFillOn = False) res3@vcWindBarbColor = "Black" res3@vcRefLengthF = 0.020 res3@vcRefAnnoOn = False res3@vcZeroFLabelOn = False res3@vcMinDistanceF = 0.04 ; <-- thins the wind barbs in NDC space (which aren't there for var2 anyway) res3@gsnScalarContour = True res3@gsnSpreadColorStart = spreadstart res3@gsnSpreadColorEnd = spreadend res3@gsnLeftString = "" res3@gsnRightString = "" res3@gsnDraw = False res3@gsnFrame = False res3@cnInfoLabelOn = False plot3 = gsn_csm_contour(xwks,var3mem,res3) overlay(plot1,plot3) ; Plot the third variable (updraft) ; var4 = "BLAH" ; ; if ( var4.eq."updraft" ) then ; var4title = "1 km AGL Updraft" ; print("Extracting variables for ensemble member "+nens) ; plotvar4_allw = wrf_user_getvar(ens,"W",0) ; plotvar4_all = plotvar4_allw(0:49,:,:) ; z = wrf_user_getvar(ens,"z",0) ; plotvar4 = wrf_user_intrp3d(plotvar4_all,z,"h",500,0.,False) ; <-- 500 m AGL updraft/downdraft ; ; print("MAX updraft at 500m AGL is: " + max(plotvar4)) ; print("MIN updraft at 500m AGL is: " + min(plotvar4)) ; ; var4mem = plotvar4 ; var4mem!0 = "lat" ; var4mem!1 = "lon" ; var4mem@lat2d = latgrid ; var4mem@lon2d = longrid ; ; res4 = True ; ; res4@cnLevelSelectionMode = "ExplicitLevels" ; res4_levels = (/ -1.0 , 1.0 /) ; ; if ( var4 .eq. "updraft" ) then ; <-- draw 40 dBZ simulated reflectivity contour ; colors4 = cmap_refl ; linecolor4 = (/"burlywood4","forestgreen"/) ; units4 = "m s~S~-1~N~" ; dashpatterns4 = (/ 2 , 0 /) ; end if ; ; res4@cnFillOn = False ; <-- changed this to false to turn off shading for var2 ; res4@cnLinesOn = True ; <-- changed this to True to turn on var2 contours ; res4@cnFillPattern = fillpattern ; res4@cnLineColors = linecolor4 ; res4@cnLineLabelFontColor = linecolor4 ; res4@cnLineDashPatterns = (/ 2 , 0 /) ; res4@cnLineThicknessF = linethickness4 ; res4@cnLineLabelsOn = False ; res4@lbLabelBarOn = True ; res4@lbTitleString = units1 ; res4@mpUSStateLineColor = "Black" ; res4@pmLabelBarParallelPosF = 0.22 ; res4@tiMainString = titlestring ; res4@vfXCStride = 2 ; res4@vfYCStride = 2 ; res4@vcGlyphStyle = "WindBarb" ; res4@vcWindBarbLineThicknessF = 0.0 ; <-- set this to 0 for var2 so it doesn't plot barbs again (even though it shouldn't since I turned cnFillOn = False) ; res4@vcWindBarbColor = "Black" ; res4@vcRefLengthF = 0.020 ; res4@vcRefAnnoOn = False ; res4@vcZeroFLabelOn = False ; res4@vcMinDistanceF = 0.04 ; <-- thins the wind barbs in NDC space (which aren't there for var2 anyway) ; res4@gsnScalarContour = True ; res4@gsnSpreadColorStart = spreadstart ; res4@gsnSpreadColorEnd = spreadend ; res4@gsnLeftString = "" ; res4@gsnRightString = "" ; res4@gsnDraw = False ; res4@gsnFrame = False ; res4@cnInfoLabelOn = False ; ; plot4 = new(2,graphic) ; ; do j = 0,1 ; res4@cnLineColor = linecolor4(j) ; res4@cnLineDashPattern = dashpatterns4(j) ; res4@cnLevels = res4_levels(j) ; ; plot4(j) = gsn_csm_contour(xwks,var4mem,res4) ; overlay(plot1,plot4(j)) ; ; end do ; ; end if print("got to here!!") draw(plot1) frame(xwks) delete(uumem) delete(vvmem) end do ; times (0400 - 0530 UTC) print("done plotting member 30 at " + hhmm) end do ; ensemble members (30) end