vis3d/vprofile.m

   1 function vprofile(p,x1,x2,tstep)
   2 % vprofile(filename,x1,x2,tstep)
   3 % compute and plot vertical profiles at given location and time
   4 % input
   5 % p             string, WRF file name, or structure p where the file was read using the command
   6 %               p=wrfatm2struct(filename,tstep);
   7 % x1            distance in m from lower left end of domain in coordinate 1
   8 % x2            distance in m from lower left end of domain in coordinate 2
   9 % tstep         the index in input files to compute this at, or empty
  10 %
  11 % examples:
  12 %             vprofile('wrfinput_d01',100,100,[])
  13 %             vprofile('wrfout_d01_0001-01-01_00:00:00',100,100,1)
  14
  15 % the data for this time step
  16 if ~isstruct(p),
  17     p=wrfatm2struct(filename,tstep);
  18 end
  19 % horizontal coordinates in terms of mesh index, cell centered
  20 i1=x1/p.dx-0.5;
  21 i2=x2/p.dy-0.5;
  22 fi1=x1/p.fdx-0.5;
  23 fi2=x2/p.fdy-0.5;
  24 % interpolate the vertical profiles
  25 altitude=interp_12(p.altitude,i1,i2);
  26
  27 disp(['Time is ',p.times{1}])
  28
  29 % wind profile
  30 u=interp_12(p.u,i1+0.5,i2);  % staggered in x1
  31 v=interp_12(p.v,i1,i2+0.5);  % staggered in x2
  32 [direction,speed]=pol(u,v); % wind from north to south is direction zero
  33
  34 layers=5;
  35 fprintf('Wind profile at (%gm %gm) from lower left corner of domain\n',x1,x2)
  36 z0 = interp2(p.z0, i1, i2);
  37 dz0=z0;  % for display
  38 fprintf('Rougness height from LANDUSE %5.3fm\n',z0);
  39 fz0 = interp2(p.fz0, fi1, fi2);
  40 if fz0 > 0,
  41 fprintf('Roughness height on fire mesh %5.3fm\n',fz0);
  42   dz0=fz0;
  43 end
  44 fwh = interp2(p.fwh, fi1, fi2);
  45 if fwh > 0,
  46   fprintf('Wind height fire mesh        %5.3fm\n',fwh);
  47 end
  48 uf = interp2(p.uf, fi1, fi2);
  49 vf = interp2(p.vf, fi1, fi2);
  50 [df,sf]=pol(uf,vf);
  51
  52 if sf<=0,
  53   fprintf('No fuel at this location, wind not interpolated.\n')
  54 end
  55
  56 fprintf('layer altitude     U       V   speed  direction (degrees)\n')
  57 fprintf('z0   %7.3f %7.3f %7.3f %7.3f\n',dz0,0,0,0)  % at roughness height
  58
  59 if fwh>0,
  60      fprintf('fwh  %7.3f %7.3f %7.3f %7.3f %7.3f\n',fwh,uf,vf,sf,df)
  61 elseif sf>0,
  62      fprintf('reduced      %7.3f %7.3f %7.3f %7.3f\n',uf,vf,sf,df)
  63 end
  64 for i=1:layers,
  65      fprintf('%3i  %7.3f %7.3f %7.3f %7.3f %7.3f\n',i,altitude(i),u(i),v(i),speed(i),direction(i))
  66 end
  67
  68 mm=[layers,size(p.u,3)];
  69 for i=1:2
  70 m=mm(i);
  71 speedfile=sprintf('speed_%i.png',m);
  72 directionfile=sprintf('direction_%i.png',m);
  73 windufile=sprintf('windu_%i.png',m);
  74 windvfile=sprintf('windv_%i.png',m);
  75 figure(1)
  76 plot(speed(1:m),altitude(1:m))
  77 if fwh>0,
  78     hold on
  79     plot(sf,fwh,'*')
  80     hold off
  81 end
  82 xlabel('speed m/s')
  83 ylabel('altitude m')
  84 title('Wind speed');
  85 print('-dpng',speedfile)
  86 figure(2)
  87 plot(direction(1:m),altitude(1:m))
  88 if fwh>0,
  89     hold on
  90     plot(df,fwh,'*')
  91     hold off
  92 end
  93 xlabel('direction degrees')
  94 ylabel('altitude m')
  95 title('Wind direction');
  96 print('-dpng',directionfile)
  97 figure(3)
  98 plot(u(1:m),altitude(1:m))
  99 if fwh>0,
 100     hold on
 101     plot(uf,fwh,'*')
 102     hold off
 103 end
 104 xlabel('speed m/s')
 105 ylabel('alitude m')
 106 title('Wind U');
 107 print('-dpng',windufile)
 108 figure(4)
 109 plot(v(1:m),altitude(1:m))
 110 if fwh>0,
 111     hold on
 112     plot(vf,fwh,'*')
 113     hold off
 114 end
 115 xlabel('speed m/s')
 116 ylabel('altitude m')
 117 title('Wind V');
 118 print('-dpng',windvfile)
 119 end
 120
 121 disp('PNG files with figures created')
 122
 123 % no temperature yet
 124 % t=interp2(t,i1,i2);
 125 % q=interp2(qvapor,i1,i2);
 126
 127 end
 128
 129 function [d,s]=pol(u,v)
 130      [d,s]=cart2pol(-v,u); % wind from north to south is direction zero
 131      d=180*d/pi;  % convert to degrees
 132 end
 133
 134