restored run_plot_wind_at_h to working again

other/Matlab/vis3d/check_uah_vah.m

@@ -33,7 +33,7 @@
 
 altw=(ph_2+phb)/9.81; % geopotential altitude at w-points
 
-[uahm,vahm]=u_v_at_h(z0,u_2,v_2,altw,height);
+[uahm,vahm]=uah_vah(z0,u_2,v_2,altw,height);
 
 err_uah=big(uah-uahm)
 err_uah=big(vah-vahm)

other/Matlab/vis3d/log_interp_vert.m

@@ -2,7 +2,7 @@
 % vertical log interpolation
 % u         values given at u poits (half eta levels)
 % hgtu      heights at u points 
-% z0        roughtness height 
+% z0        roughness height 
 % levels    heights to interpolate to (3rd index)
 % Note: the computation runs over all i,j (dimensions 1 and 2 in u) 
 % and all timesteps (dimensions 4)

other/Matlab/vis3d/plot_wind_at_h.m

@@ -1,4 +1,4 @@
-function plot_wind_at_h(p,iheights,levels,alpha,timestep,windscale,heightscale,units)
+function plot_wind_at_h(p,iheights,levels,alpha,timestep,windscale,heightscale)
 x=p.xlat;
 y=p.xlong;
 clf
@@ -12,14 +12,14 @@ function plot_wind_at_h(p,iheights,levels,alpha,timestep,windscale,heightscale,u
 % wind at given height
 for i=iheights(:)'
     z=p.heights(i)*ones(size(x));
-    ws=sqrt(p.uch(:,:,i,timestep).^2 + p.vch(:,:,i,timestep).^2)
+    ws=sqrt(p.uch(:,:,i,timestep).^2 + p.vch(:,:,i,timestep).^2);
     print_layer_msg(i,z,ws)
     surf(x,y,z+hgt,ws,'EdgeAlpha',alpha*0,'FaceAlpha',alpha)
     hold on
 end
 for i=levels(:)'
     z=p.altitude(:,:,i,timestep)-hgt;
-    ws=sqrt(p.uc(:,:,i,timestep).^2 + p.vc(:,:,i,timestep).^2)
+    ws=sqrt(p.uc(:,:,i,timestep).^2 + p.vc(:,:,i,timestep).^2);
     print_layer_msg(i,z,ws)
     surf(x,y,z+hgt,ws,'EdgeAlpha',alpha*0,'FaceAlpha',alpha)
     hold on

other/Matlab/vis3d/run_plot_wind_at_h.m

@@ -7,8 +7,8 @@
 % heights       array of heights (m) above the terrain
 % levels        indices of discretization levels
 % alpha         transparency, between 0 and 1
-% windscale     max wind for scale (m/s)
-% heightscale   max height for z-axis scale (m)
+% wscale        max wind speed for color scale (m/s)
+% hscale        max height for z-axis scale (m)
 %
 % out:
 % p          structure with the arrays displayed, and more
@@ -18,6 +18,7 @@
 % run_plot_wind_at_h('wrfout_d01_2006-02-23_12:43:00',35,6.1,[],1.0,10.0,10)
 % run_plot_wind_at_h('wrfout_d01_2006-02-23_12:43:00',35,[6.1],[1:4],0.5,10.0,10)
 % for i=1:35,run_plot_wind_at_h('wrfout_d01_2006-02-23_12:43:00',i,[1,6.1],[],1,10.0,10),end
+% for i=1:35,run_plot_wind_at_h('wrfout_d01_0001-01-01_00:00:00',i,[1,6.1],[],1,10.0,210),end
 
 % run_plot_at_h
 % wrfout='fireflux_small/wrfout_d01_2006-02-23_12:43:00';

other/Matlab/vis3d/wind_uv_at_h.m

@@ -1,5 +1,5 @@
 function p=wind_uv_at_h(p,heights)
-% p=add2struct_wind_at_h(p,heights,levels)
+% p=add2struct_wind_at_h(p,levels)
 % read 
 % in: 
 % p           structure from wrfatm2struct
@@ -13,15 +13,9 @@
 end
 
 % log interpolation of the wind at center points to height
-p.uch=log_interp_vert(p.uc,p.alt_at_w,p.z0,heights);
-p.vch=log_interp_vert(p.vc,p.alt_at_w,p.z0,heights);
+p.uch=log_interp_vert(p.uc,p.height,p.z0,heights);
+p.vch=log_interp_vert(p.vc,p.height,p.z0,heights);
 
 p.heights=heights;
 
 end
-function s=size1(a,ndim)
-% s=size1(a,ndim)
-% return size(a) extended by 1 to ndim dimensions
-s=size(a);
-s=[s,ones(1,ndim-length(s))];
-end 

other/Matlab/vis3d/wrfatm2struct.m

@@ -21,8 +21,13 @@
 % add altitude
 p.alt_at_w=(p.ph+p.phb)/9.81; % geopotential altitude at w-points
 p.altitude=(p.alt_at_w(:,:,1:end-1,:)+p.alt_at_w(:,:,2:end,:))*0.5; % interpolate to center altitude
+% subtract the altitude of the ground to get height (above the ground)
+for k=1:size(p.altitude,3)
+    p.height(:,:,k)=p.altitude(:,:,k)-p.alt_at_w(:,:,1);
+end
+
 
-% add wind
+% add wind at centers (theta points)
 p.uc = 0.5*(p.u(1:end-1,:,:,:) + p.u(2:end,:,:,:));
 p.vc = 0.5*(p.v(:,1:end-1,:,:) + p.v(:,2:end,:,:));
 p.wc = 0.5*(p.w(:,:,1:end-1,:) + p.w(:,:,2:end,:));