vis3d/wrfatm2struct.m

   1 function p=wrfatm2struct(filename,timesteps)
   2 % p=wrfatm2struct(filename)
   3 % read WRF NetCDF file atmosphere mesh variables of interest as fields in p
   4 % in:
   5 % filename    string, e.g. 'wrfout_d01_0001-01-01_00:00:00'
   6 % timesteps   vector of integers; use [] for all
   7 % out:
   8 % p           structure with several arrays from WRF and few more
   9 % additional fields in p computed:
  10 % p.alt_at_w  altitude at w-points
  11 % p.altitude  altitide interpolated to theta-points (centers of 3D cells)
  12 % p.uc,vc,wc  wind interpolated to theta-points (centers of 3D cells)
  13 % p.times     reformatted as as cell array of strings
  14 if ~exist('timesteps','var'),
  15     timesteps=[];
  16 end
  17 p=nc2struct(filename,{'U','V','W','PH','PHB','HGT','Z0','T',...
  18     'FZ0','FWH','UF','VF',...
  19     'XLONG','XLAT','GRNHFX','FGRNHFX','FXLONG','FXLAT','Times'},...
  20     {'DX','DY'},timesteps);
  21
  22 % add altitude
  23 p.alt_at_w=(p.ph+p.phb)/9.81; % geopotential altitude at w-points
  24 p.altitude=(p.alt_at_w(:,:,1:end-1,:)+p.alt_at_w(:,:,2:end,:))*0.5; % interpolate to center altitude
  25 % subtract the altitude of the ground to get height (above the ground)
  26 for k=1:size(p.altitude,3)
  27     p.height(:,:,k,:)=p.altitude(:,:,k,:)-p.alt_at_w(:,:,1,:);
  28 end
  29
  30 % refinement ratios
  31 p.sr_x=size(p.fxlong,1)/size(p.xlong,1);
  32 p.sr_y=size(p.fxlong,2)/size(p.xlong,2);
  33
  34 % fire mesh step
  35 p.fdx = p.dx/p.sr_x;
  36 p.fdy = p.dy/p.sr_y;
  37
  38 % add wind at centers (theta points)
  39 p.uc = 0.5*(p.u(1:end-1,:,:,:) + p.u(2:end,:,:,:));
  40 p.vc = 0.5*(p.v(:,1:end-1,:,:) + p.v(:,2:end,:,:));
  41 p.wc = 0.5*(p.w(:,:,1:end-1,:) + p.w(:,:,2:end,:));
  42
  43
  44 for i=1:size(p.times,2), % make Times readable
  45     times{i}=char(p.times(:,1)');
  46 end
  47 p.times=times;
  48 %test
  49 max_rel_err_hgt=big(squeeze(p.alt_at_w(:,:,1,:))-squeeze(p.hgt))/max(big(p.hgt),realmin);
  50 fprintf('relative error of geopotential ground altitude %g\n',max_rel_err_hgt)
  51 end