vis3d/log_interp_vert.m

   1 function v_levels=log_interp_vert(u,hgtu,z0,levels)
   2 % vertical log interpolation
   3 % u         values given at u poits (half eta levels)
   4 % hgtu      heights at u points
   5 % z0        roughness height
   6 % levels    heights to interpolate to (3rd index)
   7 % Note: the computation runs over all i,j (dimensions 1 and 2 in u)
   8 % and all timesteps (dimensions 4)
   9
  10 % extend u by zeros at the ground
  11 s=size1(u,4);
  12 u0=zeros(s(1),s(2),s(3)+1,s(4));
  13 u0(:,:,2:end,:)=u;
  14
  15 levels=levels(:);
  16 if any(levels<=0),
  17     disp(levels)
  18     error('levels must be positive for log interpolation')
  19 end
  20 log_levels=log(levels); % interpolate to there
  21 n=length(levels);
  22 v_levels=zeros(s(1),s(2),n,s(4));
  23 for t=1:s(4)
  24     for i=1:s(1)
  25         for j=1:s(2)
  26             heights=[z0(i,j,t);squeeze(hgtu(i,j,:,t))];
  27             if heights(2)<=heights(1),
  28                 disp(heights)
  29                 error('first level must be higher than z0')
  30             end
  31             if any(heights<=0),
  32                 disp(heights)
  33                 error('heights must be positive for log interpolation')
  34             end
  35             log_heights=log(heights);
  36             u_ijt=squeeze(u0(i,j,:,t));
  37             v_levels(i,j,:,t)=interp1(log_heights,u_ijt,log_levels);
  38         end
  39     end
  40 end
  41 end