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','PHB','PH'},{'DX','DY'},frame);
   6
   7 % WRF atm wind at cell centers (theta points)
   8 % see https://wiki.openwfm.org/wiki/How_to_interpret_WRF_variables
   9 % cell centers coordinates
  10 [nx,ny,nz]=size(p.w);
  11 nz=nz-1;
  12 p.xc=zeros(nx,ny,nz);
  13 p.yc=zeros(nx,ny,nz);
  14 p.zc=zeros(nx,ny,nz);
  15 ph = p.phb + p.ph;
  16 p.zc =0.5*(ph(:,:,1:end-1)+ph(:,:,2:end))/9.81
  17 [p.xc(:,:,1),p.yc(:,:,1)]=ndgrid(p.dx*([1:nx]-0.5),p.dy*([1:ny]-0.5));
  18 for k=2:nz
  19     p.xc(:,:,k) = p.xc(:,:,1);
  20     p.yc(:,:,k) = p.yc(:,:,1);
  21 end
  22 % average the wind from staggered grid to cell centers
  23 p.uc = 0.5*(p.u(1:end-1,:,:) + p.u(2:end,:,:));
  24 p.vc = 0.5*(p.v(:,1:end-1,:) + p.v(:,2:end,:));
  25 p.wc = 0.5*(p.w(:,:,1:end-1) + p.w(:,:,2:end));
  26 p.uvwc=sqrt(p.uc.^2+p.vc.^2+p.wc.^2); % 3D velocity
  27 p.uvc=sqrt(p.uc.^2+p.vc.^2); % horizontal velocity
  28 %*** fmw arrays
  29 % fmw cell center coordinates, same as wind
  30 [nxf,nyf,nzf]=size(p.u0_fmw);
  31 % fire step sizes
  32 p.dxf=p.dx/(nxf/nx);
  33 p.dyf=p.dx/(nxf/nx);
  34 % coordinates of centers
  35 p.xcf=zeros(nxf,nyf,nzf);
  36 p.ycf=zeros(nxf,nyf,nzf);
  37 p.zcf=zeros(nxf,nyf,nzf);
  38 [p.xcf(:,:,1),p.ycf(:,:,1)]=ndgrid(p.dxf*([1:nxf]-0.5),p.dyf*([1:nyf]-0.5));
  39 ht0=[0;p.ht_fmw]; p.htcf=(ht0(2:end)+ht0(1:end-1))/2;
  40 p.zcf(:,:,1) = p.zsf + p.htcf(1);
  41 for k=2:nzf
  42     p.xcf(:,:,k) = p.xcf(:,:,1);
  43     p.ycf(:,:,k) = p.ycf(:,:,1);
  44     p.zcf(:,:,k) = p.zsf + p.htcf(k);
  45 end
  46 %***
  47 % fmw wind velocities
  48 p.uv0=sqrt(p.u0_fmw.^2+p.v0_fmw.^2);
  49 p.uvw0=sqrt(p.u0_fmw.^2+p.v0_fmw.^2+p.w0_fmw.^2);
  50 p.uv0=sqrt(p.u_fmw.^2+p.v_fmw.^2);
  51 p.uvw0f=sqrt(p.u_fmw.^2+p.v_fmw.^2+p.w_fmw.^2);
  52 % interpolated to fire wind height fwh
  53 p.uvf=sqrt(p.uf.^2+p.vf.^2);
  54 end