| 1 |
% |
| 2 |
% %Id: $ |
| 3 |
% |
| 4 |
% Ed Hill |
| 5 |
% |
| 6 |
% Create a netCDF file defining a CS32 grid for SCRIP. |
| 7 |
|
| 8 |
clear all |
| 9 |
close all |
| 10 |
|
| 11 |
debug_lev = 0; |
| 12 |
|
| 13 |
%================================================================== |
| 14 |
% Read the tile00?.mitgrid files |
| 15 |
gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ... |
| 16 |
'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' }; |
| 17 |
ne = 32; |
| 18 |
nep1 = ne + 1; |
| 19 |
iface = 1; |
| 20 |
for iface = 1:6 |
| 21 |
fname = sprintf('CS32/tile%03d.mitgrid', iface); |
| 22 |
gid = fopen(fname, 'r', 'ieee-be'); |
| 23 |
tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]); |
| 24 |
fclose(gid); |
| 25 |
% surf(tmp(:,:,1)), view(2), shading interp |
| 26 |
% for jj = 1:length(gvars) |
| 27 |
for jj = 1:7 |
| 28 |
comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ... |
| 29 |
[gvars{jj}], iface, jj); |
| 30 |
eval(comm); |
| 31 |
end |
| 32 |
end |
| 33 |
% surf(XC(:,:,1)), view(2), shading interp |
| 34 |
% subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2) |
| 35 |
% subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2) |
| 36 |
% surf(YC(:,:,1)), view(2), shading interp |
| 37 |
% surf(XG(:,:,1)), view(2), shading interp |
| 38 |
% surf(YG(:,:,1)), view(2), shading interp |
| 39 |
%================================================================== |
| 40 |
|
| 41 |
|
| 42 |
grid_size = ne * ne * 6; |
| 43 |
grid_rank = 2; |
| 44 |
grid_corners = 4; |
| 45 |
|
| 46 |
! rm -f remap_grid_CS32.cdf |
| 47 |
|
| 48 |
nc = netcdf(['remap_grid_CS32.cdf'], 'clobber'); |
| 49 |
nc.title = 'Global CS32 grid'; |
| 50 |
nc.author = 'Ed Hill'; |
| 51 |
nc('grid_size') = grid_size; |
| 52 |
nc('grid_rank') = grid_rank; |
| 53 |
nc('grid_corners') = grid_corners; |
| 54 |
|
| 55 |
nc{'grid_dims'} = ncint('grid_rank'); |
| 56 |
nc{'grid_dims'}(:) = [ ne 6*ne ]; |
| 57 |
|
| 58 |
nc{'grid_imask'} = ncint('grid_size'); |
| 59 |
nc{'grid_imask'}(:) = ones(length(grid_size),1); |
| 60 |
|
| 61 |
nc{'grid_center_lat'} = ncdouble('grid_size'); |
| 62 |
nc{'grid_center_lon'} = ncdouble('grid_size'); |
| 63 |
nc{'grid_center_lat'}.units = 'degrees'; |
| 64 |
nc{'grid_center_lon'}.units = 'degrees'; |
| 65 |
nc{'grid_corner_lat'} = ncdouble('grid_size', 'grid_corners'); |
| 66 |
nc{'grid_corner_lon'} = ncdouble('grid_size', 'grid_corners'); |
| 67 |
nc{'grid_corner_lat'}.units = 'degrees'; |
| 68 |
nc{'grid_corner_lon'}.units = 'degrees'; |
| 69 |
|
| 70 |
cen_lat = zeros(grid_size,1); |
| 71 |
cen_lon = zeros(grid_size,1); |
| 72 |
cor_lat = zeros(grid_size,4); |
| 73 |
cor_lon = zeros(grid_size,4); |
| 74 |
j = 0; |
| 75 |
for k = 1:6 |
| 76 |
for iy = 1:ne |
| 77 |
for ix = 1:ne |
| 78 |
j = j + 1; |
| 79 |
cen_lat(j) = YC(ix,iy,k); |
| 80 |
cen_lon(j) = XC(ix,iy,k); |
| 81 |
% Note: grid points *MUST* be ordered CCW and here |
| 82 |
% we start from the lower-left corner |
| 83 |
a = [ YG(ix,iy,k) YG(ix+1,iy,k) YG(ix+1,iy+1,k) YG(ix,iy+1,k) ]; |
| 84 |
b = [ XG(ix,iy,k) XG(ix+1,iy,k) XG(ix+1,iy+1,k) XG(ix,iy+1,k) ]; |
| 85 |
% fixup the poles |
| 86 |
for ii = 1:4 |
| 87 |
if (abs(b(ii)-XC(ix,iy,k)) > 150) |
| 88 |
b(ii) = b(ii) + sign(XC(ix,iy,k))*360; |
| 89 |
end |
| 90 |
if (90 - abs(a(ii))) < 0.01 |
| 91 |
b(ii) = XC(ix,iy,k); |
| 92 |
end |
| 93 |
end |
| 94 |
cor_lat(j,:) = a; |
| 95 |
cor_lon(j,:) = b; |
| 96 |
% if debug_lev > 1 |
| 97 |
% pause(1) |
| 98 |
% hold on |
| 99 |
% plot([cor_lon(j,:) cor_lon(j,1)],... |
| 100 |
% [cor_lat(j,:) cor_lat(j,1)],'o-b') |
| 101 |
% plot(cen_lon(j),cen_lat(j),'xr') |
| 102 |
% hold off |
| 103 |
% end |
| 104 |
end |
| 105 |
end |
| 106 |
end |
| 107 |
nc{'grid_center_lat'}(:) = cen_lat; |
| 108 |
nc{'grid_center_lon'}(:) = cen_lon; |
| 109 |
nc{'grid_corner_lat'}(:) = cor_lat; |
| 110 |
nc{'grid_corner_lon'}(:) = cor_lon; |
| 111 |
|
| 112 |
close(nc) |
| 113 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph