femwind/read_fmw_3d.m

   1 function p=read_fmw_3d(path,frame)
   2 % p=read_fmw_3d(path,frame)
   3 % read initial and resulting wind at cell centers and compute their coordinates
   4 p=nc2struct(path,{'U0_FMW','V0_FMW','W0_FMW','HT_FMW','ZSF','FWH',...
   5     'U_FMW','W_FMW','V_FMW','UF','VF','FWH','U','V','W','ZS','HGT',...
   6     'PHB','PH','P','PB'},{'DX','DY'},frame);
   7 if ~isempty(p.u0_fmw)
   8     disp('fmw variables, computing derived')
   9 else
  10     return
  11 end
  12 % WRF atm wind at cell centers (theta points)
  13 % see https://wiki.openwfm.org/wiki/How_to_interpret_WRF_variables
  14 % cell centers coordinates
  15 [nx,ny,nz]=size(p.w);
  16 nz=nz-1;
  17 p.xc=zeros(nx,ny,nz);
  18 p.yc=zeros(nx,ny,nz);
  19 p.zc=zeros(nx,ny,nz);
  20 ph = p.phb + p.ph;
  21 p.zw=ph/9.81;
  22 p.zc =0.5*(ph(:,:,1:end-1)+ph(:,:,2:end))/9.81
  23 [p.xc(:,:,1),p.yc(:,:,1)]=ndgrid(p.dx*([1:nx]-0.5),p.dy*([1:ny]-0.5));
  24 for k=2:nz
  25     p.xc(:,:,k) = p.xc(:,:,1);
  26     p.yc(:,:,k) = p.yc(:,:,1);
  27 end
  28 % average the wind from staggered grid to cell centers
  29 p.uc = 0.5*(p.u(1:end-1,:,:) + p.u(2:end,:,:));
  30 p.vc = 0.5*(p.v(:,1:end-1,:) + p.v(:,2:end,:));
  31 p.wc = 0.5*(p.w(:,:,1:end-1) + p.w(:,:,2:end));
  32 p.uvwc=sqrt(p.uc.^2+p.vc.^2+p.wc.^2); % 3D velocity
  33 p.uvc=sqrt(p.uc.^2+p.vc.^2); % horizontal velocity
  34 %*** fmw arrays
  35 % fmw cell center coordinates, same as wind
  36 [nxf,nyf,nzf]=size(p.u0_fmw);
  37 % fire step sizes
  38 p.dxf=p.dx/(nxf/nx);
  39 p.dyf=p.dx/(nxf/nx);
  40 % coordinates of centers
  41 p.xcf=zeros(nxf,nyf,nzf);
  42 p.ycf=zeros(nxf,nyf,nzf);
  43 p.zcf=zeros(nxf,nyf,nzf);
  44 [p.xcf(:,:,1),p.ycf(:,:,1)]=ndgrid(p.dxf*([1:nxf]-0.5),p.dyf*([1:nyf]-0.5));
  45 ht0=[0;p.ht_fmw]; p.htcf=(ht0(2:end)+ht0(1:end-1))/2;
  46 p.zcf(:,:,1) = p.zsf + p.htcf(1);
  47 for k=2:nzf
  48     p.xcf(:,:,k) = p.xcf(:,:,1);
  49     p.ycf(:,:,k) = p.ycf(:,:,1);
  50     p.zcf(:,:,k) = p.zsf + p.htcf(k);
  51 end
  52 %***
  53 % fmw wind velocities
  54 p.uv0=sqrt(p.u0_fmw.^2+p.v0_fmw.^2);
  55 p.uvw0=sqrt(p.u0_fmw.^2+p.v0_fmw.^2+p.w0_fmw.^2);
  56 p.uvf=sqrt(p.u_fmw.^2+p.v_fmw.^2);
  57 p.uvwf=sqrt(p.u_fmw.^2+p.v_fmw.^2+p.w_fmw.^2);
  58 % interpolated to fire wind height fwh
  59 p.ufvf=sqrt(p.uf.^2+p.vf.^2);
  60 end