quicwind/wind2flux_trans.m

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