debug/uah_vah.m

   1 % function [uah,vah]=u_v_at_h(z0,u,v,altw,heights)
   2 % [uah,vah]=u_v_at_h(z0,u,v,altw,heights)
   3 % interpolate wind on staggered grid to given
   4 %
   5 % in:
   6 % z0        2D roughness length at cell centers (under w points)
   7 % u         3D wind, staggered in x direction
   8 % v         3D wind, staggered in y direction
   9 % altw      3D altitude, at w points
  10 % heights   1D array of heights above the terrain to interpolate to
  11 %
  12 % out:
  13 % uah       3D uah(:,:,i) is u interpolated to heights(i)
  14 % vah       3D vah(:,:,i) is v interpolated to heights(i)
  15
  16 % staggered dimension x
  17 %                u(1)---z0(1)---u(2)---z0(2)--- ...--u(n)---z0(n)---u(n+1)
  18 % interpolated:                z0u(2)    ...        z0u(n)
  19
  20 % interpolate z0 to under u and v points
  21 z0u = 0.5*(z0(1:end-1,:)+z0(2:end,:));
  22 z0v = 0.5*(z0(:,1:end-1)+z0(:,2:end));
  23
  24 % interpolate altitude to under u and v points
  25 altub=0.5*(altw(1:end-1,:,:)+altw(2:end,:,:));
  26 altvb=0.5*(altw(:,1:end-1,:)+altw(:,2:end,:));
  27
  28 % interpolate altitude vertically to u and v points
  29 hgtu=0.5*(altub(:,:,1:end-1)+altub(:,:,2:end));
  30 hgtv=0.5*(altvb(:,:,1:end-1)+altvb(:,:,2:end));
  31
  32 % subtract the altitude of the ground to get height (above the ground)
  33 for k=1:size(hgtu,3)
  34     hgtu(:,:,k)=hgtu(:,:,k)-altub(:,:,1);
  35     hgtv(:,:,k)=hgtv(:,:,k)-altvb(:,:,1);
  36 end
  37
  38 % log interpolation of the wind at u and v points to height
  39 ua=log_interp_vert(u(2:end-1,:,:),hgtu,z0u,height);
  40 va=log_interp_vert(v(:,2:end-1,:),hgtv,z0v,height);
  41
  42 % extend by extrapolation in staggered dimension
  43 for i=length(height(:)):-1:1
  44     h=height(i);
  45     uah(:,:,i)=[2*ua(1,:,i)-ua(2,:,i);ua(:,:,i);2*ua(end,:,i)-ua(end-1,:,i)];
  46     vah(:,:,i)=[2*va(:,1,i)-va(:,2,i),va(:,:,i),2*va(:,end,i)-va(:,end-1,i)];
  47 end
  48
  49 % end