matlab/cs_grid/merccube.m

ca96c02aee Alis*0001 function [] = merccube(XX,YY,C)
062b178e55 Alis*0002 % merccube(x,y,c)
ca96c02aee Alis*0003 %
062b178e55 Alis*0004 % Plots cubed-sphere data in mercator projection. (x,y) are
ca96c02aee Alis*0005 % coordinates, c is cell-centered scalar to be plotted.
062b178e55 Alis*0006 % All arrays (x,y,c) should have dimensions of NxNx6 or 6NxN.
                0007 %
                0008 % The default plotting mode is shading faceted. Using this or
                0009 % shading flat, (x,y) should be the coordinates of grid-corners
                0010 % and can legitimately have dimension (N+1)x(N+1)x6.
                0011 %
                0012 % If using shading interp, then (x,y) must be the coordinates of
                0013 % the cell centers.
ca96c02aee Alis*0014 %
                0015 % e.g.
                0016 %
                0017 % xg=rdmds('XG');
                0018 % yg=rdmds('YG');
                0019 % ps=rdmds('Eta.0000000000');
                0020 % mercube(xg,yg,ps);
                0021 % colorbar;xlabel('Longitude');ylabel('Latitude');
062b178e55 Alis*0022 %
                0023 % xc=rdmds('XC');
                0024 % yc=rdmds('YC');
                0025 % mercube(xc,yc,ps);shading interp
f8d374081f Jean*0026 %
                0027 % Written by adcroft@.mit.edu, 2004.
ca96c02aee Alis*0028 if max(max(max(YY)))-min(min(min(YY))) < 3*pi
                0029  X=XX*180/pi;
                0030  Y=YY*180/pi;
                0031 else
                0032  X=XX;
                0033  Y=YY;
                0034 end
                0035 Q=C;
                0036
22586ff4d2 Alis*0037 if ndims(Q)==2 & size(Q,1)==6*size(Q,2)
ca96c02aee Alis*0038  [nx ny nt]=size(X);
                0039  X=permute( reshape(X,[nx/6 6 ny]),[1 3 2]);
                0040  Y=permute( reshape(Y,[nx/6 6 ny]),[1 3 2]);
                0041  Q=permute( reshape(Q,[nx/6 6 ny]),[1 3 2]);
22586ff4d2 Alis*0042 elseif ndims(Q)==3 & size(Q,2)==6
ca96c02aee Alis*0043  X=permute( X,[1 3 2]);
                0044  Y=permute( Y,[1 3 2]);
                0045  Q=permute( Q,[1 3 2]);
22586ff4d2 Alis*0046 elseif ndims(Q)==3 & size(Q,3)==6
ca96c02aee Alis*0047  [nx ny nt]=size(X);
                0048 else
                0049  size(XX)
                0050  size(YY)
                0051  size(C)
                0052  error('Dimensions should be 2 or 3 dimensions: NxNx6, 6NxN or Nx6xN');
                0053 end
                0054
                0055 if size(X,1)==size(Q,1)
                0056  X(end+1,:,:)=NaN;
                0057  X(:,end+1,:)=NaN;
                0058  X(end,:,[1 3 5])=X(1,:,[2 4 6]);
                0059  X(:,end,[2 4 6])=X(:,1,[3 5 1]);
                0060  X(:,end,[1 3 5])=squeeze(X(1,end:-1:1,[3 5 1]));
                0061  X(end,:,[2 4 6])=squeeze(X(end:-1:1,1,[4 6 2]));
df2ee6ab0f Dani*0062  X(1,end,[1 3 5]) = X(1,1,1);
                0063  X(end,1,[2 4 6]) = X(end,end,2);
ca96c02aee Alis*0064  Y(end+1,:,:)=NaN;
                0065  Y(:,end+1,:)=NaN;
                0066  Y(end,:,[1 3 5])=Y(1,:,[2 4 6]);
                0067  Y(:,end,[2 4 6])=Y(:,1,[3 5 1]);
                0068  Y(:,end,[1 3 5])=squeeze(Y(1,end:-1:1,[3 5 1]));
                0069  Y(end,:,[2 4 6])=squeeze(Y(end:-1:1,1,[4 6 2]));
df2ee6ab0f Dani*0070  Y(1,end,[1 3 5]) = Y(end,end,1);
                0071  Y(end,1,[2 4 6]) = Y(1,1,2);
ca96c02aee Alis*0072 end
                0073 [nx ny nt]=size(X);
                0074
22586ff4d2 Alis*0075 Q(end+1,:,:)=NaN;
                0076 Q(:,end+1,:)=NaN;
ca96c02aee Alis*0077 Q(end,:,[1 3 5])=Q(1,:,[2 4 6]);
                0078 Q(:,end,[2 4 6])=Q(:,1,[3 5 1]);
                0079 Q(:,end,[1 3 5])=squeeze(Q(1,end:-1:1,[3 5 1]));
                0080 Q(end,:,[2 4 6])=squeeze(Q(end:-1:1,1,[4 6 2]));
                0081
                0082 hnx=ceil(nx/2);
                0083 hny=ceil(ny/2);
                0084
                0085 for k=1:6;
                0086  i=1:hnx;
                0087  x=longitude(X(i,:,k));
                0088  pcolor(x,Y(i,:,k),Q(i,:,k))
                0089  axis([-180 180 -90 90])
                0090  hold on
                0091  if max(max(max(x)))>180
                0092   pcolor(x-360,Y(i,:,k),Q(i,:,k))
                0093  end
                0094  i=hnx:nx;
                0095  x=longitude(X(i,:,k));
                0096  pcolor(x,Y(i,:,k),Q(i,:,k))
                0097  if max(max(max(x)))>180
                0098   pcolor(x-360,Y(i,:,k),Q(i,:,k))
                0099  end
                0100 end
                0101 hold off