/[MITgcm]/MITgcm_contrib/onedegcube/read510.m
ViewVC logotype

Annotation of /MITgcm_contrib/onedegcube/read510.m

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


Revision 1.1 - (hide annotations) (download)
Thu Oct 4 17:56:07 2007 UTC (17 years, 9 months ago) by molod
Branch: MAIN
CVS Tags: HEAD
One degree cube - files used to create, resulting grid files, gcm-generated
grid output files to verify required symmetries

1 molod 1.1 %
2     % Read in CS 510 input fields
3     xc1 = zeros(511,511,6);
4     yc1 = zeros(511,511,6);
5     dxf1 = zeros(511,511,6);
6     dyf1 = zeros(511,511,6);
7     ra1 = zeros(511,511,6);
8     xg1 = zeros(511,511,6);
9     yg1 = zeros(511,511,6);
10     dxv1 = zeros(511,511,6);
11     dyu1 = zeros(511,511,6);
12     raz1 = zeros(511,511,6);
13     dxc1 = zeros(511,511,6);
14     dyc1 = zeros(511,511,6);
15     raw1 = zeros(511,511,6);
16     ras1 = zeros(511,511,6);
17     dxg1 = zeros(511,511,6);
18     dyg1 = zeros(511,511,6);
19     for iface=1:6
20     fid = fopen(['tile00' num2str(iface) '.mitgrid'],'r','b');
21     xc1(:,:,iface) = fread(fid,[511 511],'real*8');
22     yc1(:,:,iface) = fread(fid,[511 511],'real*8');
23     dxf1(:,:,iface) = fread(fid,[511 511],'real*8');
24     dyf1(:,:,iface) = fread(fid,[511 511],'real*8');
25     ra1(:,:,iface) = fread(fid,[511 511],'real*8');
26     xg1(:,:,iface) = fread(fid,[511 511],'real*8');
27     yg1(:,:,iface) = fread(fid,[511 511],'real*8');
28     dxv1(:,:,iface) = fread(fid,[511 511],'real*8');
29     dyu1(:,:,iface) = fread(fid,[511 511],'real*8');
30     raz1(:,:,iface) = fread(fid,[511 511],'real*8');
31     dxc1(:,:,iface) = fread(fid,[511 511],'real*8');
32     dyc1(:,:,iface) = fread(fid,[511 511],'real*8');
33     raw1(:,:,iface) = fread(fid,[511 511],'real*8');
34     ras1(:,:,iface) = fread(fid,[511 511],'real*8');
35     dxg1(:,:,iface) = fread(fid,[511 511],'real*8');
36     dyg1(:,:,iface) = fread(fid,[511 511],'real*8');
37     end
38     %
39     % Now output (approx 1 deg) grid
40     ratio = 5;
41     xc2 = zeros(103,103,6);
42     yc2 = zeros(103,103,6);
43     dxf2 = zeros(103,103,6);
44     dyf2 = zeros(103,103,6);
45     ra2 = zeros(103,103,6);
46     xg2 = zeros(103,103,6);
47     yg2 = zeros(103,103,6);
48     dxv2 = zeros(103,103,6);
49     dyu2 = zeros(103,103,6);
50     raz2 = zeros(103,103,6);
51     dxc2 = zeros(103,103,6);
52     dyc2 = zeros(103,103,6);
53     raw2 = zeros(103,103,6);
54     ras2 = zeros(103,103,6);
55     dxg2 = zeros(103,103,6);
56     anglecos2 = zeros(103,103,6);
57     anglesin2 = zeros(103,103,6);
58     dyg2 = zeros(103,103,6);
59     % Create left, x-center, bottom and y-center indeces into 510 grid
60     ileft=[1:5:510];
61     icent=ileft+2;
62     jbott=[1:5:510];
63     jcent=jbott+2;
64     % First do the interior of each face
65     for iface=1:6
66     % lats and lons
67     xc2(1:102,1:102,iface)=xc1(icent,jcent,iface);
68     yc2(1:102,1:102,iface)=yc1(icent,jcent,iface);
69     xg2(1:102,1:102,iface)=xg1(ileft,jbott,iface);
70     yg2(1:102,1:102,iface)=yg1(ileft,jbott,iface);
71     for ipnt=0:ratio-1
72     dxf2(1:102,1:102,iface)=dxf2(1:102,1:102,iface) + dxf1(ileft+ipnt,jcent,iface);
73     dxg2(1:102,1:102,iface)=dxg2(1:102,1:102,iface) + dxg1(ileft+ipnt,jbott,iface);
74     dxv2(2:102,1:102,iface)=dxv2(2:102,1:102,iface) + dxv1(icent(1:101)+ipnt+1,jbott(1:102),iface);
75     dxc2(2:102,1:102,iface)=dxc2(2:102,1:102,iface) + dxc1(icent(1:101)+ipnt+1,jcent(1:102),iface);
76     end
77     % dx's and dy's
78     for jpnt=0:ratio-1
79     dyf2(1:102,1:102,iface)=dyf2(1:102,1:102,iface) + dyf1(icent,jbott+jpnt,iface);
80     dyg2(1:102,1:102,iface)=dyg2(1:102,1:102,iface) + dyg1(ileft,jbott+jpnt,iface);
81     dyu2(1:102,2:102,iface)=dyu2(1:102,2:102,iface) + dyu1(ileft(1:102),jcent(1:101)+jpnt+1,iface);
82     dyc2(1:102,2:102,iface)=dyc2(1:102,2:102,iface) + dyc1(icent(1:102),jcent(1:101)+jpnt+1,iface);
83     end
84     % Construct values that we don't know yet - use grid symmetries for edge dx's and dy's
85     dxc2(1,1:102,iface)=dxc1(1,jcent(1:102),iface) + 2.*dxc1(2,jcent(1:102),iface) + 2.*dxc1(3,jcent(1:102),iface);
86     dxv2(1,1:102,iface)=dxv1(1,jbott(1:102),iface) + 2.*dxv1(2,jbott(1:102),iface) + 2.*dxv1(3,jbott(1:102),iface);
87     dyc2(1:102,1,iface)=dyc1(icent(1:102),1,iface) + 2.*dyc1(icent(1:102),2,iface) + 2.*dyc1(icent(1:102),3,iface);
88     dyu2(1:102,1,iface)=dyu1(ileft(1:102),1,iface) + 2.*dyu1(ileft(1:102),2,iface) + 2.*dyu1(ileft(1:102),3,iface);
89     % Areas
90     for ipnt=0:ratio-1
91     for jpnt=0:ratio-1
92     ra2(1:102,1:102,iface)=ra2(1:102,1:102,iface) + ra1(ileft+ipnt,jbott+jpnt,iface);
93     ras2(1:102,2:102,iface)=ras2(1:102,2:102,iface) + ras1(ileft(1:102)+ipnt,jcent(1:101)+jpnt+1,iface);
94     raw2(2:102,1:102,iface)=raw2(2:102,1:102,iface) + raw1(icent(1:101)+ipnt+1,jbott(1:102)+jpnt,iface);
95     raz2(2:102,2:102,iface)=raz2(2:102,2:102,iface) + raz1(icent(1:101)+ipnt+1,jcent(1:101)+jpnt+1,iface);
96     end
97     end
98     % Construct values that we don't know yet - use grid symmetries for edge areas
99     for ipnt=0:ratio-1
100     ras2(1:102,1,iface) = ras2(1:102,1,iface) + ...
101     ras1(ileft(1:102)+ipnt,1,iface) + 2.*ras1(ileft(1:102)+ipnt,2,iface) + 2.*ras1(ileft(1:102)+ipnt,3,iface);
102     end
103     for jpnt=0:ratio-1
104     raw2(1,1:102,iface) = raw2(1,1:102,iface) + ...
105     raw1(1,jbott(1:102)+jpnt,iface) + 2.*raw1(2,jbott(1:102)+jpnt,iface) + 2.*raw1(3,jbott(1:102)+jpnt,iface);
106     end
107     % Construct values that we don't know yet - vorticity point area, along face edges than at face corner
108     for jpnt=0:ratio-1
109     raz2(1,2:102,iface) = raz2(1,2:102,iface) + ...
110     raz1(1,jcent(1:101)+jpnt+1,iface) + 2.*raz1(2,jcent(1:101)+jpnt+1,iface) + 2.*raz1(3,jcent(1:101)+jpnt+1,iface);
111     end
112     for ipnt=0:ratio-1
113     raz2(2:102,1,iface) = raz2(2:102,1,iface) + ...
114     raz1(icent(1:101)+ipnt+1,1,iface) + 2.*raz1(icent(1:101)+1+ipnt,2,iface) + 2.*raz1(icent(1:101)+1+ipnt,3,iface);
115     end
116     raz2(1,1,iface) = sum(raz1(1:3,1,iface),1) + 2.*sum(raz1(1,2:3,iface),2) + ...
117     3.*sum(raz1(2:3,3,iface),1) + sum(raz1(2:3,2,iface),1);
118     end
119     % Now the exchanges to fill up the extra column and row
120     % Exchange for lon and lat's at center and corners - no directions (xc is always xc)
121     xc2(1:102,103,[1 3 5])=xc2(1,102:-1:1,[3 5 1]);
122     xc2(103,1:102,[1 3 5])=xc2(1,1:102,[2 4 6]);
123     xc2(1:102,103,[2 4 6])=xc2(1:102,1,[3 5 1]);
124     xc2(103,1:102,[2 4 6])=xc2(102:-1:1,1,[4 6 2]);
125     xc2(103,103,[1 2 3 4 5 6])=xc2(1,1,[3 4 5 6 1 2]);
126     yc2(1:102,103,[1 3 5])=yc2(1,102:-1:1,[3 5 1]);
127     yc2(103,1:102,[1 3 5])=yc2(1,1:102,[2 4 6]);
128     yc2(1:102,103,[2 4 6])=yc2(1:102,1,[3 5 1]);
129     yc2(103,1:102,[2 4 6])=yc2(102:-1:1,1,[4 6 2]);
130     yc2(103,103,[1 2 3 4 5 6])=yc2(1,1,[3 4 5 6 1 2]);
131     xg2(1:102,103,[1 3 5])=xg2(1,102:-1:1,[3 5 1]);
132     xg2(103,1:102,[1 3 5])=xg2(1,1:102,[2 4 6]);
133     xg2(1:102,103,[2 4 6])=xg2(1:102,1,[3 5 1]);
134     xg2(103,1:102,[2 4 6])=xg2(102:-1:1,1,[4 6 2]);
135     xg2(103,103,[1 2 3 4 5 6])=xg2(1,1,[3 4 5 6 1 2]);
136     yg2(1:102,103,[1 3 5])=yg2(1,102:-1:1,[3 5 1]);
137     yg2(103,1:102,[1 3 5])=yg2(1,1:102,[2 4 6]);
138     yg2(1:102,103,[2 4 6])=yg2(1:102,1,[3 5 1]);
139     yg2(103,1:102,[2 4 6])=yg2(102:-1:1,1,[4 6 2]);
140     yg2(103,103,[1 2 3 4 5 6])=yg2(1,1,[3 4 5 6 1 2]);
141     % Exchange for areas at center - no directions
142     ra2(1:102,103,[1 3 5])=ra2(1,102:-1:1,[3 5 1]);
143     ra2(103,1:102,[1 3 5])=ra2(1,1:102,[2 4 6]);
144     ra2(1:102,103,[2 4 6])=ra2(1:102,1,[3 5 1]);
145     ra2(103,1:102,[2 4 6])=ra2(102:-1:1,1,[4 6 2]);
146     ra2(103,103,[1 2 3 4 5 6])=ra2(1,1,[3 4 5 6 1 2]);
147     % Exchange for areas at south and west edges - no directions
148     ras2(1:102,103,[1 3 5])=raw2(1,102:-1:1,[3 5 1]);
149     ras2(103,1:102,[1 3 5])=ras2(1,1:102,[2 4 6]);
150     ras2(1:102,103,[2 4 6])=ras2(1:102,1,[3 5 1]);
151     ras2(103,1:102,[2 4 6])=raw2(102:-1:1,1,[4 6 2]);
152     ras2(103,103,[1 2 3 4 5 6])=ras2(1,1,[3 4 5 6 1 2]);
153     raw2(1:102,103,[1 3 5])=ras2(1,102:-1:1,[3 5 1]);
154     raw2(103,1:102,[1 3 5])=raw2(1,1:102,[2 4 6]);
155     raw2(1:102,103,[2 4 6])=raw2(1:102,1,[3 5 1]);
156     raw2(103,1:102,[2 4 6])=ras2(102:-1:1,1,[4 6 2]);
157     raw2(103,103,[1 2 3 4 5 6])=raw2(1,1,[3 4 5 6 1 2]);
158     % Exchange for areas at vorticity points - no directions
159     % Special case: define upper edge assuming symmetry in y
160     % Done because the edge values are ambiguously defined
161     % for quantities defined on a corner
162     %%%%raz2(1:102,103,[1 3 5])=raz2(1,102:-1:1,[3 5 1]);
163     %%%%raz2(1:102,103,[2 4 6])=raz2(1:102,1,[3 5 1]);
164     %%%%raz2(103,1:102,[1 3 5])=raz2(1,1:102,[2 4 6]);
165     %%%%raz2(103,1:102,[2 4 6])=raz2(102:-1:1,1,[4 6 2]);
166     raz2(1:102,103,:)=raz2(1:102,1,:);
167     raz2(103,1:102,:)=raz2(1,1:102,:);
168     raz2(103,103,[1 2 3 4 5 6])=raz2(1,1,[3 4 5 6 1 2]);
169     %
170     % Exchange for dx's and dy's at center - direction but no sign (dx is sometimes from dy)
171     dxf2(1:102,103,[1 3 5])=dyf2(1,102:-1:1,[3 5 1]);
172     dxf2(103,1:102,[1 3 5])=dxf2(1,1:102,[2 4 6]);
173     dxf2(1:102,103,[2 4 6])=dxf2(1:102,1,[3 5 1]);
174     dxf2(103,1:102,[2 4 6])=dyf2(102:-1:1,1,[4 6 2]);
175     dxf2(103,103,[1 2 3 4 5 6])=dxf2(1,1,[3 4 5 6 1 2]);
176     dyf2(1:102,103,[1 3 5])=dxf2(1,102:-1:1,[3 5 1]);
177     dyf2(103,1:102,[1 3 5])=dyf2(1,1:102,[2 4 6]);
178     dyf2(1:102,103,[2 4 6])=dyf2(1:102,1,[3 5 1]);
179     dyf2(103,1:102,[2 4 6])=dxf2(102:-1:1,1,[4 6 2]);
180     dyf2(103,103,[1 2 3 4 5 6])=dyf2(1,1,[3 4 5 6 1 2]);
181     dxg2(1:102,103,[1 3 5])=dyg2(1,102:-1:1,[3 5 1]);
182     dxg2(103,1:102,[1 3 5])=dxg2(1,1:102,[2 4 6]);
183     dxg2(1:102,103,[2 4 6])=dxg2(1:102,1,[3 5 1]);
184     dxg2(103,1:102,[2 4 6])=dyg2(102:-1:1,1,[4 6 2]);
185     dxg2(103,103,[1 2 3 4 5 6])=dxg2(1,1,[3 4 5 6 1 2]);
186     dyg2(1:102,103,[1 3 5])=dxg2(1,102:-1:1,[3 5 1]);
187     dyg2(103,1:102,[1 3 5])=dyg2(1,1:102,[2 4 6]);
188     dyg2(1:102,103,[2 4 6])=dyg2(1:102,1,[3 5 1]);
189     dyg2(103,1:102,[2 4 6])=dxg2(102:-1:1,1,[4 6 2]);
190     dyg2(103,103,[1 2 3 4 5 6])=dyg2(1,1,[3 4 5 6 1 2]);
191     % Exchange for dx's and dy's at edges - direction but no sign (dx is sometimes from dy)
192     dxc2(1:102,103,[1 3 5])=dyc2(1,102:-1:1,[3 5 1]);
193     dxc2(103,1:102,[1 3 5])=dxc2(1,1:102,[2 4 6]);
194     dxc2(1:102,103,[2 4 6])=dxc2(1:102,1,[3 5 1]);
195     dxc2(103,1:102,[2 4 6])=dyc2(102:-1:1,1,[4 6 2]);
196     dxc2(103,103,[1 2 3 4 5 6])=dxc2(1,1,[3 4 5 6 1 2]);
197     dyc2(1:102,103,[1 3 5])=dxc2(1,102:-1:1,[3 5 1]);
198     dyc2(103,1:102,[1 3 5])=dyc2(1,1:102,[2 4 6]);
199     dyc2(1:102,103,[2 4 6])=dyc2(1:102,1,[3 5 1]);
200     dyc2(103,1:102,[2 4 6])=dxc2(102:-1:1,1,[4 6 2]);
201     dyc2(103,103,[1 2 3 4 5 6])=dyc2(1,1,[3 4 5 6 1 2]);
202     % Special case: define upper edge assuming symmetry in y
203     % Done because the edge values are ambiguously defined
204     % for quantities defined on a corner
205     %%%%dxv2(1:102,103,[1 3 5])=dyu2(1,102:-1:1,[3 5 1]);
206     %%%%dxv2(103,1:102,[1 3 5])=dxv2(1,1:102,[2 4 6]);
207     %%%%dxv2(1:102,103,[2 4 6])=dxv2(1:102,1,[3 5 1]);
208     %%%%dxv2(103,1:102,[2 4 6])=dyu2(102:-1:1,1,[4 6 2]);
209     %%%%dyu2(1:102,103,[1 3 5])=dxv2(1,102:-1:1,[3 5 1]);
210     %%%%dyu2(103,1:102,[1 3 5])=dyu2(1,1:102,[2 4 6]);
211     %%%%dyu2(1:102,103,[2 4 6])=dyu2(1:102,1,[3 5 1]);
212     %%%%dyu2(103,1:102,[2 4 6])=dxv2(102:-1:1,1,[4 6 2]);
213     dxv2(1:102,103,:)=dxv2(1:102,1,:);
214     dxv2(103,1:102,:)=dxv2(1,1:102,:);
215     dxv2(103,103,[1 2 3 4 5 6])=dxv2(1,1,[3 4 5 6 1 2]);
216     dyu2(1:102,103,:)=dyu2(1:102,1,:);
217     dyu2(103,1:102,:)=dyu2(1,1:102,:);
218     dyu2(103,103,[1 2 3 4 5 6])=dyu2(1,1,[3 4 5 6 1 2]);
219     %
220     % All done with grid values - now compute angles
221     ygarg1=reshape(permute(yg2(1:102,1:102,:),[1 3 2]),[612 102]);
222     racarg2=reshape(permute(ra2(1:102,1:102,:),[1 3 2]),[612 102]);
223     dxgarg3=reshape(permute(dxg2(1:102,1:102,:),[1 3 2]),[612 102]);
224     dygarg4=reshape(permute(dyg2(1:102,1:102,:),[1 3 2]),[612 102]);
225     [AngleCS,AngleSN] = cubeCalcAngle(ygarg1,racarg2,dxgarg3,dygarg4);
226     anglecos2(1:102,1:102,:)=permute(reshape(AngleCS,[102 6 102]),[1 3 2]);
227     anglesin2(1:102,1:102,:)=permute(reshape(AngleSN,[102 6 102]),[1 3 2]);
228     %
229     % And now write it all out
230     for iface=1:6
231     fid = fopen(['onedegcube.face00' num2str(iface) '.bin'],'w','b');
232     fwrite(fid,xc2(:,:,iface),'real*8');
233     fwrite(fid,yc2(:,:,iface),'real*8');
234     fwrite(fid,dxf2(:,:,iface),'real*8');
235     fwrite(fid,dyf2(:,:,iface),'real*8');
236     fwrite(fid,ra2(:,:,iface),'real*8');
237     fwrite(fid,xg2(:,:,iface),'real*8');
238     fwrite(fid,yg2(:,:,iface),'real*8');
239     fwrite(fid,dxv2(:,:,iface),'real*8');
240     fwrite(fid,dyu2(:,:,iface),'real*8');
241     fwrite(fid,raz2(:,:,iface),'real*8');
242     fwrite(fid,dxc2(:,:,iface),'real*8');
243     fwrite(fid,dyc2(:,:,iface),'real*8');
244     fwrite(fid,raw2(:,:,iface),'real*8');
245     fwrite(fid,ras2(:,:,iface),'real*8');
246     fwrite(fid,dxg2(:,:,iface),'real*8');
247     fwrite(fid,dyg2(:,:,iface),'real*8');
248     fwrite(fid,anglecos2(:,:,iface),'real*8');
249     fwrite(fid,anglesin2(:,:,iface),'real*8');
250     fclose(fid);
251     end

  ViewVC Help
Powered by ViewVC 1.1.22