/[MITgcm]/MITgcm_contrib/enderton/Diagnostics/DiagUtility/merccube_mod.m
ViewVC logotype

Annotation of /MITgcm_contrib/enderton/Diagnostics/DiagUtility/merccube_mod.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Thu Aug 11 14:38:51 2005 UTC (19 years, 11 months ago) by enderton
Branch: MAIN
CVS Tags: HEAD
Same as merccube.m, but can handle small grid irregularities such as values near 0 and 180 and -180, but not the exact values.

1 enderton 1.1 function [] = merccube(XX,YY,C)
2     % merccube(x,y,c)
3     %
4     % Plots cubed-sphere data in mercator projection. (x,y) are
5     % coordinates, c is cell-centered scalar to be plotted.
6     % All arrays (x,y,c) should have dimensions of NxNx6 or 6NxN.
7     %
8     % The default plotting mode is shading faceted. Using this or
9     % shading flat, (x,y) should be the coordinates of grid-corners
10     % and can legitimately have dimension (N+1)x(N+1)x6.
11     %
12     % If using shading interp, then (x,y) must be the coordinates of
13     % the cell centers.
14     %
15     % e.g.
16     %
17     % xg=rdmds('XG');
18     % yg=rdmds('YG');
19     % ps=rdmds('Eta.0000000000');
20     % mercube(xg,yg,ps);
21     % colorbar;xlabel('Longitude');ylabel('Latitude');
22     %
23     % xc=rdmds('XC');
24     % yc=rdmds('YC');
25     % mercube(xc,yc,ps);shading interp
26    
27     XXdim = size(XX);
28     XX = reshape(XX,[prod(XXdim),1]);
29     XX(find(abs(XX) < 0.01)) = 0;
30     XX(find(abs(XX) > 179.99)) = 180;
31     XX = reshape(XX,XXdim);
32    
33     if max(max(max(YY)))-min(min(min(YY))) < 3*pi
34     X=XX*180/pi;
35     Y=YY*180/pi;
36     else
37     X=XX;
38     Y=YY;
39     end
40     Q=C;
41    
42     if ndims(Q)==2 & size(Q,1)==6*size(Q,2)
43     [nx ny nt]=size(X);
44     X=permute( reshape(X,[nx/6 6 ny]),[1 3 2]);
45     Y=permute( reshape(Y,[nx/6 6 ny]),[1 3 2]);
46     Q=permute( reshape(Q,[nx/6 6 ny]),[1 3 2]);
47     elseif ndims(Q)==3 & size(Q,2)==6
48     X=permute( X,[1 3 2]);
49     Y=permute( Y,[1 3 2]);
50     Q=permute( Q,[1 3 2]);
51     elseif ndims(Q)==3 & size(Q,3)==6
52     [nx ny nt]=size(X);
53     else
54     size(XX)
55     size(YY)
56     size(C)
57     error('Dimensions should be 2 or 3 dimensions: NxNx6, 6NxN or Nx6xN');
58     end
59    
60     if size(X,1)==size(Q,1)
61     X(end+1,:,:)=NaN;
62     X(:,end+1,:)=NaN;
63     X(end,:,[1 3 5])=X(1,:,[2 4 6]);
64     X(:,end,[2 4 6])=X(:,1,[3 5 1]);
65     X(:,end,[1 3 5])=squeeze(X(1,end:-1:1,[3 5 1]));
66     X(end,:,[2 4 6])=squeeze(X(end:-1:1,1,[4 6 2]));
67     Y(end+1,:,:)=NaN;
68     Y(:,end+1,:)=NaN;
69     Y(end,:,[1 3 5])=Y(1,:,[2 4 6]);
70     Y(:,end,[2 4 6])=Y(:,1,[3 5 1]);
71     Y(:,end,[1 3 5])=squeeze(Y(1,end:-1:1,[3 5 1]));
72     Y(end,:,[2 4 6])=squeeze(Y(end:-1:1,1,[4 6 2]));
73     end
74     [nx ny nt]=size(X);
75    
76     Q(end+1,:,:)=NaN;
77     Q(:,end+1,:)=NaN;
78     Q(end,:,[1 3 5])=Q(1,:,[2 4 6]);
79     Q(:,end,[2 4 6])=Q(:,1,[3 5 1]);
80     Q(:,end,[1 3 5])=squeeze(Q(1,end:-1:1,[3 5 1]));
81     Q(end,:,[2 4 6])=squeeze(Q(end:-1:1,1,[4 6 2]));
82    
83     hnx=ceil(nx/2);
84     hny=ceil(ny/2);
85    
86     for k=1:6;
87     i=1:hnx;
88     x=longitude(X(i,:,k));
89     pcolor(x,Y(i,:,k),Q(i,:,k))
90     axis([-180 180 -90 90])
91     hold on
92     if max(max(max(x)))>180
93     pcolor(x-360,Y(i,:,k),Q(i,:,k))
94     end
95     i=hnx:nx;
96     x=longitude(X(i,:,k));
97     pcolor(x,Y(i,:,k),Q(i,:,k))
98     if max(max(max(x)))>180
99     pcolor(x-360,Y(i,:,k),Q(i,:,k))
100     end
101     end
102     hold off

  ViewVC Help
Powered by ViewVC 1.1.22