quicwind/wind2flux.m

   1 function F=wind2flux(U,X)
   2 % F=wind2flux(U,X)
   3 % compute fluxes at midpoints of side in a hexa grid
   4 % assuming all sides in the y direction are straight vertical
   5 % while horizontal sides may be slanted
   6 % in:
   7 %     U{l}(i,j,k) flow component l at lower side midpoint of cell (i,j,k)
   8 %                 such as output from grad3z
   9 %     X{l}(i,j,k) coordinate l of node (i,j,k), if l=1,2 does not depend on k
  10 % out:
  11 %     F{l}(i,j,k) flux in normal direction through lower side midpoint of cell (i,j,k)
  12 %                 flux through the bottom is zero
  13 %
  14 % Jan Mandel, August 2019
  15
  16 u = U{1}; v = U{2}; w = U{3};
  17 x = X{1}; y = X{2}; z = X{3};
  18 [nx1,ny1,nz1] = size(x);
  19 nx = nx1-1; ny = ny1-1; nz = nz1-1;
  20
  21 % test inputs
  22 if ndims(u)~=3|ndims(v)~=3|ndims(w)~=3|ndims(x)~=3|ndims(y)~=3|ndims(z)~=3,
  23     error('wind2flux: all input arrays must be 3D')
  24 end
  25 if any(size(u)~=[nx1,ny,nz])|any(size(v)~=[nx,ny1,nz])|any(size(w)~=[nx,ny,nz1])
  26     error('wind2flux: arrays u v w must be staggered dimensioned with x y z')
  27 end
  28
  29 check_mesh(X)
  30
  31 err=false;
  32 for k=2:nz+1
  33     if any(x(:,:,k)~=x(:,:,1))|any(y(:,:,k)~=y(:,:,1))
  34         x(:,:,k)=x(:,:,1);
  35         y(:,:,k)=y(:,:,1);
  36         err=true;
  37     end
  38 end
  39 if err,
  40     error('wind2flux: arrays x and y must be constant in 3rd dimension')
  41 end
  42
  43
  44 %                                             ^ z,w,k
  45 %     (i,j+1,k+1)---------(i+1,j+1,k+1        |
  46 %     /  |               / |                  |    /
  47 %    /   |              /  |                  |   / y,v,j
  48 %   /    |             /   |                  |  /
  49 %(i,j,k+1)----------(i+1,j,k+1)               | /
  50 %  |    /|            |    |                  |--------> x,u,i
  51 %  |  dy |            |    |
  52 %  | /   |            |    |
  53 %  |/    |            |    |
  54 %  |   (i,j+1,k)---------(i+1,j+1,k)
  55 %  |   /              |  /
  56 %  |  /               | /
  57 %  | /                |/
  58 %(i,j,k)----------(i+1,j,k)
  59
  60
  61 s = cell_sizes(X);
  62
  63 % flux through vertical sides left to rigth (u,x,i direction)
  64 f_u = u .* s.area_u;
  65
  66 % flux through vertical sides front to back (v,y,j direction)
  67 f_v = v .* s.area_v;
  68
  69 % slope in x direction
  70 dzdx = zeros(nx,ny+1,nz+1);
  71 for k=1:nz+1
  72     for j=1:ny+1
  73         for i=1:nx
  74             dzdx(i,j,k) = ((z(i+1,j,k)-z(i,j,k))/(x(i+1,j,k)-x(i,j,k)));
  75         end
  76     end
  77 end
  78
  79 % slope in y direction
  80 dzdy = zeros(nx+1,ny,nz+1);
  81 for k=1:nz+1
  82     for j=1:ny
  83         for i=1:nx+1
  84             dzdy(i,j,k) = ((z(i,j+1,k)-z(i,j,k))/(y(i,j+1,k)-y(i,j,k)));
  85         end
  86     end
  87 end
  88 % normal flux through slanted horizontal bottom to top j direction
  89 f_w = zeros(nx,ny,nz+1);
  90 % continue one layer up
  91 u(:,:,nz+1)=u(:,:,nz);
  92 v(:,:,nz+1)=v(:,:,nz);
  93 s.dz_at_u(:,:,nz+1)=s.dz_at_u(:,:,nz);
  94 s.dz_at_v(:,:,nz+1)=s.dz_at_v(:,:,nz);
  95 for k=2:nz+1   % zero normal flux on the ground, k=1
  96     for j=1:ny
  97         for i=1:nx
  98             % mesh horizontal mesh cell sizes of layer k
  99             % u, v averaged from above and below
 100             u_at_k = 0.5*( (u(i,j,k)   + u(i+1,j,k))  *s.dz_at_u(i,j,k) ...
 101                          + (u(i,j,k-1) + u(i+1,j,k-1))*s.dz_at_u(i,j,k-1) ...
 102                          )/ (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k-1));
 103             v_at_k = 0.5*( (v(i,j,k)   +v(i,j+1,k))   *s.dz_at_v(i,j,k) ...
 104                          + (v(i,j,k-1) +v(i,j+1,k-1)) *s.dz_at_v(i,j,k-1) ...
 105                          )/ (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k-1));
 106             % average slope at midpoints from two sides
 107             dzdx_at_k = 0.5*(dzdx(i,j,k) + dzdx(i,j+1,k));
 108             dzdy_at_k = 0.5*(dzdy(i,j,k) + dzdy(i+1,j,k));
 109             f_w(i,j,k)=s.area_w(i,j,k)*(w(i,j,k) - ...
 110                 dzdx_at_k * u_at_k - dzdy_at_k * v_at_k);
 111         end
 112     end
 113 end
 114 F = {f_u, f_v, f_w};
 115 end
 116
 117
 118
 119
 120
 121