Contents of /MITgcm_contrib/high_res_cube/matlab-grid-generator/README

code obtained from faulks.lcs.mit.edu:/scratch/adcroft/cubed_sphere

path('/home/dimitri/matlab/adcroft/bin',path);
tic
[dxg,dyg,dxf,dyf,dxc,dyc,dxv,dyu,Ec,Eu,Ev,Ez,latC,lonC,latG,lonG,...
          Q11,Q22,Q12, TUu,TUv,TVu,TVv ]=gengrid_fn(576,4,'conf','c',0,0);
toc

timings                           nireas  faulks
gengrid_fn(576,2,'conf','c',0,0)  47.9    60.0
gengrid_fn(576,4,'conf','c',0,0)  163.9   

current 80S-80N, 1/4-deg isotropic grid:
1440 * 1088 = 1566720 grid points
27.7987 km at Equator, 4.8272 km at 80N

tripolar grid:
2304 grid boxes around Equator
2*pi*6371/2304 = 17.3742 km grid at Equator
switch to bipolar is at 60N where grid is 17.3742*cos(60) =  8.6871 km
at North Pole grid size is 17.3742/pi = 5.5304 km
smallest distance in Canadian Archipelago is approximately 17.3742/6 = 2.8957 km
worst aspect ratio is approximately 3
number of grid points in sperical polar region is 2304*960 = 2211840
number of grid points in bipolar grid is 1152*1152 = 1327104
total number of grid points is 3538944
advantages, grid is perfectly orthogonal
there are no singularities
no modifications to sea-ice thermodynamics

by comparison
the cubed sphere with 2304 grid boxes around Equator
using Alistair's eq 16 has:
576*576*6 = 1990656 grid points (56% of tripolar grid, 127% of 1/4-deg grid)
largest distance is 22.9683 (132% of tripolar grid)
smallest distance is 3.2633 (112% of tripolar grid, 68% of 1/4-deg grid)
worst aspect ratio is 1.6162 (53% of tripolar grid)
but orthogonality is not as good as that of tripolar grid
and sea-ice dynamics need some work
also for the tripolar grid, removing land point will
result in much more dramatic reduction in number of grid
points because the three singularities are over land.


table
first row is min, max, and mean of dxg and dyg
second row is min, max, and mean of dxg./dyg
third row is min, max, and mean of angles, excluding cube corners

conf, nratio=4
    2.1460   18.0293   15.9182
    0.9172    1.0903    1.0000
   89.9978   98.5344   90.0562

conf, nratio=2
    2.1460   18.0293   15.9182
    0.9167    1.0909    1.0000
   89.9978   98.5344   90.0562

q=0, nratio=4
    2.2544   17.5336   15.9487
    0.9172    1.0903    1.0006
   89.9977   98.5344   90.0554

q=0, nratio=2
    2.2544   17.4069   15.9468
    0.9167    1.0909    1.0006
   89.9977   98.5344   90.0554

q=1, nratio=4 (has NaNs)
   13.6379   67.1398   16.1959
    0.2031    4.9226    1.0552
   89.9972   98.5343   90.0497

q=1, nratio=2
   13.6310   67.1393   16.1925
    0.2031    4.9238    1.0551
   89.9972   98.5343   90.0497

q=1/2, nratio=4
    2.8990   21.0171   16.0604
    0.6992    1.4302    1.0137
   89.9974   98.5344   90.0527

q=1/2, nratio=2
    2.8958   21.0001   16.0613
    0.6999    1.4287    1.0136
   89.9974   98.5344   90.0527

q=7/8, nratio=4 (has NaNs)
   13.6379   67.1398   16.1959
    0.2031    4.9226    1.0552
   89.9972   98.5343   90.0497

q=7/8, nratio=2
   13.6310   67.1393   16.1925
    0.2031    4.9238    1.0551
   89.9972   98.5343   90.0497

q=i3, nratio=4
   13.1144   65.1937   16.1921
    0.2092    4.7807    1.0549
   89.9972   98.5343   90.0497

q=i3, nratio=2
   12.4259   62.6094   16.1915
    0.2179    4.5899    1.0546
   89.9973   98.5343   90.0498

tan, nratio=4
    3.2633   22.9683   16.0912
    0.6188    1.6162    1.0203
   89.9974   98.5344   90.0519

tan, nratio=2
    3.2633   22.9683   16.0912
    0.6188    1.6162    1.0203
   89.9974   98.5344   90.0519

tan2, nratio=4
    2.2237   17.7883   15.9327
    0.9174    1.0900    1.0001
   89.9978   98.5344   90.0558

tan2, nratio=2
    2.2237   17.7883   15.9327
    0.9169    1.0906    1.0001
   89.9978   98.5344   90.0558

new, nratio=4
    2.5095   18.8618   15.9191
    0.8889    1.1250    1.0000
   89.9978   98.5344   90.0561

new, nratio=2
    2.8863   20.9484   15.9201
    0.8000    1.2500    1.0002
   89.9978   98.5344   90.0561

 
gengrids                                        % Does it all
gengrids_fn(nx,nratio,projection,orientation,plotfigs,writedata)
[dxg,dyg,dxf,dyf,dxc,dyc,dxv,dyu,Ec,Eu,Ev,Ez,latC,lonC,latG,lonG,...
          Q11,Q22,Q12, TUu,TUv,TVu,TVv ]=gengrid_fn(24,2,'conf','c',1,0);

[X,Y,Z,x,y]=read_scubdat(fn,n);                 % Read SCUB0002.DAT file
[X,Y,Z]=map_xy2xyz(x,y);                        % Map (x,y) -> (X,Y,Z)
[X,Y,Z]=gmap_xy2xyz(x,y);                       % Map (x,y) -> (X,Y,Z)

[lon,lat]=map_xyz2lonlat(X,Y,Z);                % Map (X,Y,Z) -> (lon,lat)
[X,Y,Z]=map_lonlat2xyz(lon,lat);                % Map (lon,lat) -> (X,Y,Z)
[lat,lon]=calc_geocoords(lx,ly,tile)            % Map (x,y) -> (lat,lon)

[X,Y,Z]=rotate_about_xaxis(lX,lY,lZ,angle);     % Rotate about X axis
[X,Y,Z]=rotate_about_yaxis(lX,lY,lZ,angle);     % Rotate about Y axis
[X,Y,Z]=rotate_about_zaxis(lX,lY,lZ,angle);     % Rotate about Z axis

V=coord2vector(XG,YG,ZG);                       % Convert [X,Y,Z] to vector V

angle=angle_between_coords(X1,Y1,Z1,X2,Y2,Z2);  % Angles between coords
angle=angle_between_vectors(V1,V2);             % Angles between vectors V1,V2
excess=excess_of_quad(V1,V2,V3,V4);             % Excess of quadrilateral

q=rescale_coordinate(q,'q=i3');                 % Re-scale single coordinate
[dx,dy,E]=calc_fvgrid(lx,ly);                   % Calculate F.V. grid info
dxg=conf_dx(q);                                 % Calculate conformal dxg
[dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio);         % Sum fine-grid dx -> dxg,...
[Ec,Ez,Ev]=reduce_E(dx,nratio);                 % Sum fine-grid E -> Ec,...

a=bari(b);                                      % Two point average in I dir.
a=barj(b);                                      % Two point average in J dir.

write_tiles('LATC',lat)                         % Write tiled files
displaytiles(lat)                               % Display tiles in comp. layout
plotcube(X,Y,Z,C)                               % Display tiles as 3-D sphere
merccube(lonG,latG,C)                           % Display tiles Mercator proj.
1	code obtained from faulks.lcs.mit.edu:/scratch/adcroft/cubed_sphere
2
3	path('/home/dimitri/matlab/adcroft/bin',path);
4	tic
5	[dxg,dyg,dxf,dyf,dxc,dyc,dxv,dyu,Ec,Eu,Ev,Ez,latC,lonC,latG,lonG,...
6	Q11,Q22,Q12, TUu,TUv,TVu,TVv ]=gengrid_fn(576,4,'conf','c',0,0);
7	toc
8
9	timings nireas faulks
10	gengrid_fn(576,2,'conf','c',0,0) 47.9 60.0
11	gengrid_fn(576,4,'conf','c',0,0) 163.9
12
13	current 80S-80N, 1/4-deg isotropic grid:
14	1440 * 1088 = 1566720 grid points
15	27.7987 km at Equator, 4.8272 km at 80N
16
17	tripolar grid:
18	2304 grid boxes around Equator
19	2pi6371/2304 = 17.3742 km grid at Equator
20	switch to bipolar is at 60N where grid is 17.3742*cos(60) = 8.6871 km
21	at North Pole grid size is 17.3742/pi = 5.5304 km
22	smallest distance in Canadian Archipelago is approximately 17.3742/6 = 2.8957 km
23	worst aspect ratio is approximately 3
24	number of grid points in sperical polar region is 2304*960 = 2211840
25	number of grid points in bipolar grid is 1152*1152 = 1327104
26	total number of grid points is 3538944
27	advantages, grid is perfectly orthogonal
28	there are no singularities
29	no modifications to sea-ice thermodynamics
30
31	by comparison
32	the cubed sphere with 2304 grid boxes around Equator
33	using Alistair's eq 16 has:
34	5765766 = 1990656 grid points (56% of tripolar grid, 127% of 1/4-deg grid)
35	largest distance is 22.9683 (132% of tripolar grid)
36	smallest distance is 3.2633 (112% of tripolar grid, 68% of 1/4-deg grid)
37	worst aspect ratio is 1.6162 (53% of tripolar grid)
38	but orthogonality is not as good as that of tripolar grid
39	and sea-ice dynamics need some work
40	also for the tripolar grid, removing land point will
41	result in much more dramatic reduction in number of grid
42	points because the three singularities are over land.
43
44
45
46	table
47	first row is min, max, and mean of dxg and dyg
48	second row is min, max, and mean of dxg./dyg
49	third row is min, max, and mean of angles, excluding cube corners
50
51	conf, nratio=4
52	2.1460 18.0293 15.9182
53	0.9172 1.0903 1.0000
54	89.9978 98.5344 90.0562
55
56	conf, nratio=2
57	2.1460 18.0293 15.9182
58	0.9167 1.0909 1.0000
59	89.9978 98.5344 90.0562
60
61	q=0, nratio=4
62	2.2544 17.5336 15.9487
63	0.9172 1.0903 1.0006
64	89.9977 98.5344 90.0554
65
66	q=0, nratio=2
67	2.2544 17.4069 15.9468
68	0.9167 1.0909 1.0006
69	89.9977 98.5344 90.0554
70
71	q=1, nratio=4 (has NaNs)
72	13.6379 67.1398 16.1959
73	0.2031 4.9226 1.0552
74	89.9972 98.5343 90.0497
75
76	q=1, nratio=2
77	13.6310 67.1393 16.1925
78	0.2031 4.9238 1.0551
79	89.9972 98.5343 90.0497
80
81	q=1/2, nratio=4
82	2.8990 21.0171 16.0604
83	0.6992 1.4302 1.0137
84	89.9974 98.5344 90.0527
85
86	q=1/2, nratio=2
87	2.8958 21.0001 16.0613
88	0.6999 1.4287 1.0136
89	89.9974 98.5344 90.0527
90
91	q=7/8, nratio=4 (has NaNs)
92	13.6379 67.1398 16.1959
93	0.2031 4.9226 1.0552
94	89.9972 98.5343 90.0497
95
96	q=7/8, nratio=2
97	13.6310 67.1393 16.1925
98	0.2031 4.9238 1.0551
99	89.9972 98.5343 90.0497
100
101	q=i3, nratio=4
102	13.1144 65.1937 16.1921
103	0.2092 4.7807 1.0549
104	89.9972 98.5343 90.0497
105
106	q=i3, nratio=2
107	12.4259 62.6094 16.1915
108	0.2179 4.5899 1.0546
109	89.9973 98.5343 90.0498
110
111	tan, nratio=4
112	3.2633 22.9683 16.0912
113	0.6188 1.6162 1.0203
114	89.9974 98.5344 90.0519
115
116	tan, nratio=2
117	3.2633 22.9683 16.0912
118	0.6188 1.6162 1.0203
119	89.9974 98.5344 90.0519
120
121	tan2, nratio=4
122	2.2237 17.7883 15.9327
123	0.9174 1.0900 1.0001
124	89.9978 98.5344 90.0558
125
126	tan2, nratio=2
127	2.2237 17.7883 15.9327
128	0.9169 1.0906 1.0001
129	89.9978 98.5344 90.0558
130
131	new, nratio=4
132	2.5095 18.8618 15.9191
133	0.8889 1.1250 1.0000
134	89.9978 98.5344 90.0561
135
136	new, nratio=2
137	2.8863 20.9484 15.9201
138	0.8000 1.2500 1.0002
139	89.9978 98.5344 90.0561
140
141
142
143
144
145
146	gengrids % Does it all
147	gengrids_fn(nx,nratio,projection,orientation,plotfigs,writedata)
148	[dxg,dyg,dxf,dyf,dxc,dyc,dxv,dyu,Ec,Eu,Ev,Ez,latC,lonC,latG,lonG,...
149	Q11,Q22,Q12, TUu,TUv,TVu,TVv ]=gengrid_fn(24,2,'conf','c',1,0);
150
151	[X,Y,Z,x,y]=read_scubdat(fn,n); % Read SCUB0002.DAT file
152	[X,Y,Z]=map_xy2xyz(x,y); % Map (x,y) -> (X,Y,Z)
153	[X,Y,Z]=gmap_xy2xyz(x,y); % Map (x,y) -> (X,Y,Z)
154
155	[lon,lat]=map_xyz2lonlat(X,Y,Z); % Map (X,Y,Z) -> (lon,lat)
156	[X,Y,Z]=map_lonlat2xyz(lon,lat); % Map (lon,lat) -> (X,Y,Z)
157	[lat,lon]=calc_geocoords(lx,ly,tile) % Map (x,y) -> (lat,lon)
158
159	[X,Y,Z]=rotate_about_xaxis(lX,lY,lZ,angle); % Rotate about X axis
160	[X,Y,Z]=rotate_about_yaxis(lX,lY,lZ,angle); % Rotate about Y axis
161	[X,Y,Z]=rotate_about_zaxis(lX,lY,lZ,angle); % Rotate about Z axis
162
163	V=coord2vector(XG,YG,ZG); % Convert [X,Y,Z] to vector V
164
165	angle=angle_between_coords(X1,Y1,Z1,X2,Y2,Z2); % Angles between coords
166	angle=angle_between_vectors(V1,V2); % Angles between vectors V1,V2
167	excess=excess_of_quad(V1,V2,V3,V4); % Excess of quadrilateral
168
169	q=rescale_coordinate(q,'q=i3'); % Re-scale single coordinate
170	[dx,dy,E]=calc_fvgrid(lx,ly); % Calculate F.V. grid info
171	dxg=conf_dx(q); % Calculate conformal dxg
172	[dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio); % Sum fine-grid dx -> dxg,...
173	[Ec,Ez,Ev]=reduce_E(dx,nratio); % Sum fine-grid E -> Ec,...
174
175	a=bari(b); % Two point average in I dir.
176	a=barj(b); % Two point average in J dir.
177
178	write_tiles('LATC',lat) % Write tiled files
179	displaytiles(lat) % Display tiles in comp. layout
180	plotcube(X,Y,Z,C) % Display tiles as 3-D sphere
181	merccube(lonG,latG,C) % Display tiles Mercator proj.