femwind/read_fmw_wrfout.m

   1 function [X,u0]=read_fmw_wrfout(path,timestep)
   2
   3 % if timestep not exist, last timestep
   4 if ~exist('timestep','var'),
   5     timestep=size(nc2struct(path,{'Times'},{}).times,2);
   6 end
   7 % reading wrfout variables
   8 % p=nc2struct(path,{'U0_FMW','V0_FMW','W0_FMW','HT_FMW','ZSF'},{'DX','DY'},timestep);
   9 p=nc2struct(path,{'U0_FMW','V0_FMW','W0_FMW','HT_FMW','ZSF','U_FMW','W_FMW','V_FMW','UF','VF'},{'DX','DY'},timestep);
  10 % looking at sizes
  11 [nx,ny,nz]=size(p.u0_fmw);
  12 % height profile from the ground
  13 htt_fmw=zeros(nz+1,1);
  14 for z=1:nz
  15     htt_fmw(z+1)=2*p.ht_fmw(z)-htt_fmw(z);
  16 end
  17 % compute elevation at corners from midpoints
  18 zsf=midpoints2corners(p.zsf);
  19 % creating grid
  20 [Cx,Cy,Cz]=ndgrid(p.dx*[0:nx],p.dy*[0:ny],htt_fmw);
  21 % adding ZSF
  22 for z=1:nz+1
  23     Cz(:,:,z)=Cz(:,:,z)+zsf;
  24 end
  25 % generate and check mesh
  26 X={Cx,Cy,Cz};
  27 check_mesh(X);
  28 % compute initial wind at the corners from midpoints
  29 %Xu0=midpoints2corners(p.u0_fmw);
  30 %Yu0=midpoints2corners(p.v0_fmw);
  31 %Zu0=midpoints2corners(p.w0_fmw);
  32 % generate initial wind at corners of the cells
  33 %u0={Xu0,Yu0,Zu0};
  34 u0={p.u0_fmw,p.v0_fmw,p.w0_fmw};
  35
  36 end