quicwind/saddlepoint/test_WRF.m

   1 %ifile = 'wrfinput_d01';
   2 %si = nc2struct(ifile,{'XLONG','XLAT','HGT','U','V','W','PH','PHB'},{'DX','DY'});
   3 %ofile = 'wrfout_d01_2010-07-01_01:00:00';
   4 %so = nc2struct(ofile,{'U','V','W'},{});
   5 s = load('WRF_20.mat');
   6
   7 % transformation matrix from wind at x,y,z coordinates to two constant
   8 % winds at the faces
   9 T=[1,0,0;
  10    1,0,0;
  11    0,1,0;
  12    0,1,0;
  13    0,0,1;
  14    0,0,1];
  15 % intial wind at the faces
  16 v0x = (s.u0(1:end-1,:,:)+s.u0(2:end,:,:))/2;
  17 v0y = (s.v0(:,1:end-1,:)+s.v0(:,2:end,:))/2;
  18 v0z = (s.w0(:,:,1:end-1)+s.w0(:,:,2:end))/2;
  19 v0 = T*[v0x(:),v0y(:),v0z(:)]';
  20 v0 = v0(:);
  21 % final wind at the center from WRF
  22 uWx = (s.u(1:end-1,:,:)+s.u(2:end,:,:))/2;
  23 uWy = (s.v(:,1:end-1,:)+s.v(:,2:end,:))/2;
  24 uWz = (s.w(:,:,1:end-1)+s.w(:,:,2:end))/2;
  25 uW = T*[uWx(:),uWy(:),uWz(:)]';
  26 uW = uW(:);
  27 % geopotential height
  28 gh = (s.ph+s.phb)/9.81;
  29
  30 % dimensions and spacings
  31 n = size(v0x);
  32 h = [s.dx,s.dy];
  33
  34 % create meshes from input file
  35 xm = repmat(s.xlong,[1,1,n(3)]);
  36 ym = repmat(s.xlat,[1,1,n(3)]);
  37 zm = (gh(:,:,1:end-1)+gh(:,:,2:end))/2;
  38 Xm = {xm,ym,zm};
  39 [xx,yy,~] = ndgrid(h(1)*(0:n(1)),h(2)*(0:n(2)),0:n(3));
  40 ghm = (gh(1:end-1,1:end-1,:)+gh(2:end,1:end-1,:)+gh(1:end-1,2:end,:)+gh(2:end,2:end,:))/4;
  41 ghm = [ghm(1,:,:);ghm;ghm(end,:,:)];
  42 ghm = [ghm(:,1,:),ghm,ghm(:,end,:)];
  43 zz = ghm;
  44 X = {xx,yy,zz};
  45
  46 % assembly sparse matrices
  47 [A,D,E,B,C,~] = sparse_assembly(X,h);
  48
  49 % solve
  50 saddle_sparse
  51
  52 % transform resulting fluxes into cartesian winds at the centers
  53 u=E*B*v;
  54
  55 % transform WRF winds at the face to the center
  56 uWc = E*uW;
  57 v0c = E*v0;
  58
  59 % plots
  60 figure(1)
  61 mesh(xx(:,:,1),yy(:,:,1),zz(:,:,1))
  62 hold on
  63 plot_wind_3d(Xm,v0c,'Initial WRF wind',1);
  64 figure(2),
  65 mesh(xx(:,:,1),yy(:,:,1),zz(:,:,1))
  66 hold on
  67 plot_wind_3d(Xm,uWc,'Final WRF wind',1);
  68 figure(3),
  69 mesh(xx(:,:,1),yy(:,:,1),zz(:,:,1))
  70 hold on
  71 plot_wind_3d(Xm,u,'Final mass-consistent wind',1);
  72 figure(4),
  73 mesh(xx(:,:,1),yy(:,:,1),zz(:,:,1))
  74 hold on
  75 plot_wind_3d(Xm,u-uWc,'Mass-consistent vs WRF wind',1);
  76
  77 % differences
  78 %{
  79 ux = reshape(u(1:3:end),n); uy = reshape(u(2:3:end),n); uz = reshape(u(3:3:end),n);
  80 figure(5)
  81 mesh(Xm{1}(:,:,1),Xm{2}(:,:,1),ux(:,:,1)-uWx(:,:,1))
  82 xlabel('x'), ylabel('y'), zlabel('z')
  83 title('Differences in x direction')
  84 figure(6)
  85 mesh(Xm{1}(:,:,1),Xm{2}(:,:,1),uy(:,:,1)-uWy(:,:,1))
  86 xlabel('x'), ylabel('y'), zlabel('z')
  87 title('Differences in y direction')
  88 figure(7)
  89 mesh(Xm{1}(:,:,1),Xm{2}(:,:,1),uz(:,:,1)-uWz(:,:,1))
  90 xlabel('x'), ylabel('y'), zlabel('z')
  91 title('Differences in z direction')
  92 %}
  93
  94 max_diff = max(max(max(u-uWc)))
  95
  96 % save previous results
  97 v0p = v0;
  98 vp = v;
  99 up = u;
 100
 101 % new initial wind
 102 v0 = uW;
 103
 104 % solve new
 105 saddle_sparse
 106
 107 % transform resulting fluxes into cartesian winds at the centers
 108 u=E*B*v;
 109
 110 % plot
 111 figure(5),
 112 mesh(xx(:,:,1),yy(:,:,1),zz(:,:,1))
 113 hold on
 114 plot_wind_3d(Xm,u,'Final mass-consistent wind new initial',1);
 115 figure(6),
 116 mesh(xx(:,:,1),yy(:,:,1),zz(:,:,1))
 117 hold on
 118 plot_wind_3d(Xm,u-uWc,'Mass-consistent new initial vs WRF wind',1);
 119 max_diff_new = max(max(max(u-uWc)))
 120 figure(7),
 121 mesh(xx(:,:,1),yy(:,:,1),zz(:,:,1))
 122 hold on
 123 plot_wind_3d(Xm,up-u,'Mass-consistent vs Mass-consistent new initial',1);
 124 max_diff_mcs = max(max(max(up-u)))
 125
 126 % conditions
 127 plot_conditions(n,B*v0p,C,D,'WRF input')
 128 plot_conditions(n,B*uW,C,D,'WRF final')
 129 plot_conditions(n,vp,C,D,'Mass-consistent')
 130 plot_conditions(n,v,C,D,'Mass-consistent new input')