| 1 |
% latlongrid.m creates the model grid and stores the grid information |
| 2 |
% various arrays. Nominally set are: |
| 3 |
% xc, yc, xg, yg, [2d] |
| 4 |
% surface area of grid-cell on sphere (zA) [2d] |
| 5 |
% |
| 6 |
% It also stores the above information in a MATLAB file GRID.mat |
| 7 |
% usage: |
| 8 |
% load GRID.mat (loads everything above) |
| 9 |
% load GRID.mat nxc nyc xc yc (loads just these variables/arrays) |
| 10 |
% |
| 11 |
% Uses data files: n/a |
| 12 |
% Creates data files: GRID.mat |
| 13 |
% Creates arrays: xc,yc,xg,yg,xcf,ycf,zA |
| 14 |
|
| 15 |
load DATAPATH |
| 16 |
|
| 17 |
periodic=1; |
| 18 |
|
| 19 |
% read grid data |
| 20 |
vnms = { 'XC' 'YC' 'XG' 'YG' 'rA'}; |
| 21 |
d2r = pi/180.0; |
| 22 |
ginfo = {}; |
| 23 |
for k = 1:5 |
| 24 |
nc = netcdf(fullfile(gridpath,sprintf(fnamepattern,k)),'nowrite'); |
| 25 |
for j = 1:length(vnms) |
| 26 |
eval(sprintf('ginfo(k).%s = nc{''%s''}(:);',vnms{j},vnms{j})); |
| 27 |
end |
| 28 |
nc = close(nc); |
| 29 |
end |
| 30 |
|
| 31 |
% assemble grid |
| 32 |
xc = [ginfo(1).XC;ginfo(2).XC;ginfo(4).XC(end:-1:1,:)'; ... |
| 33 |
ginfo(5).XC(end:-1:1,:)']; |
| 34 |
yc = [ginfo(1).YC;ginfo(2).YC;ginfo(4).YC(end:-1:1,:)'; ... |
| 35 |
ginfo(5).YC(end:-1:1,:)']; |
| 36 |
xg = [ginfo(1).XG(1:end-1,:);ginfo(2).XG(1:end-1,:); ... |
| 37 |
ginfo(4).XG(end:-1:1,1:end-1)';ginfo(5).XG(end:-1:1,:)']; |
| 38 |
yg = [ginfo(1).YG(1:end-1,:);ginfo(2).YG(1:end-1,:); ... |
| 39 |
ginfo(4).YG(end:-1:1,1:end-1)';ginfo(5).YG(end:-1:1,:)']; |
| 40 |
zA = [ginfo(1).rA;ginfo(2).rA;ginfo(4).rA(end:-1:1,:)'; ... |
| 41 |
ginfo(5).rA(end:-1:1,:)']; |
| 42 |
|
| 43 |
% transform longitudes: |
| 44 |
i1 = find(xc<0); xc(i1) = xc(i1) + 360; |
| 45 |
i1 = find(xg<0); xg(i1) = xg(i1) + 360; |
| 46 |
warning('remember the hack for xg'); |
| 47 |
xg(end,:) = 360; xg(1,:) = 0; |
| 48 |
|
| 49 |
% do the rotation now for the side faces |
| 50 |
xc = xc+s_lon; |
| 51 |
xg = xg+s_lon; |
| 52 |
yc = yc+s_lat; |
| 53 |
yg = yg+s_lat; |
| 54 |
|
| 55 |
disp( sprintf('# of recursive levels = %i',nfinepow) ); |
| 56 |
|
| 57 |
% Create fine grid for interpolation computations |
| 58 |
nfine=2^nfinepow; |
| 59 |
|
| 60 |
% first for the cap (remember the rotation) |
| 61 |
xc3 = ginfo(3).XC+s_lon; |
| 62 |
yc3 = ginfo(3).YC+s_lat; |
| 63 |
xg3 = ginfo(3).XG+s_lon; |
| 64 |
yg3 = ginfo(3).YG+s_lat; |
| 65 |
|
| 66 |
clear xcf3 ycf3 |
| 67 |
xcf = zeros(size(xc3)*nfine); |
| 68 |
ycf = zeros(size(yc3)*nfine); |
| 69 |
xfine = [0.5:1:nfine]/nfine; |
| 70 |
disp('generating fine grid') |
| 71 |
% we need special treatment of the polar caps! |
| 72 |
% $$$ xg1 = xg3; i1 = find(xg1<0); xg1(i1) = xg1(i1) + 360; |
| 73 |
% $$$ i1 = find(xg1>=270); xg1(i1) = NaN; |
| 74 |
% $$$ xg2 = xg3; i2 = find(xg2>=90); xg2(i2) = NaN; |
| 75 |
for ix = 1:size(xc3,1) |
| 76 |
for iy = 1:size(yc3,2) |
| 77 |
xxg1 = xg3(ix:ix+1,iy:iy+1); |
| 78 |
if ~isempty(find(xxg1(:) >= 100)) & ~isempty(find(xxg1(:) < 100)) |
| 79 |
i1 = find(xxg1(:) < 0); xxg1(i1) = xxg1(i1)+360; |
| 80 |
end |
| 81 |
xx = interp2([0 1],[0 1],xxg1',xfine,xfine')'; |
| 82 |
xcf3((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = xx; |
| 83 |
yy = interp2([0 1],[0 1],yg3(ix:ix+1,iy:iy+1)',xfine,xfine')'; |
| 84 |
ycf3((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = yy; |
| 85 |
end |
| 86 |
end |
| 87 |
|
| 88 |
xcf30{1}=xcf3; |
| 89 |
for k=2:4 |
| 90 |
xcf30{k}=xcf30{k-1}; |
| 91 |
ii = find(xcf30{k}< (k-3)*90); |
| 92 |
xcf30{k}(ii) = xcf30{k}(ii) + 360; |
| 93 |
end |
| 94 |
|
| 95 |
xcf3 = [rot90(xcf30{1});xcf30{2};rot90(xcf30{3},3);rot90(xcf30{4},2)]; |
| 96 |
ycf3 = [rot90(ycf3);ycf3;rot90(rot90(rot90(ycf3)));rot90(rot90(ycf3))]; |
| 97 |
%i1 = find(xcf3<0); xcf3(i1) = xcf3(i1) + 360; |
| 98 |
|
| 99 |
|
| 100 |
xcf = zeros(size(xc)*nfine); |
| 101 |
ycf = zeros(size(yc)*nfine); |
| 102 |
xfine = [0.5:1:nfine]/nfine; |
| 103 |
for ix = 1:size(xc,1) |
| 104 |
for iy = 1:size(yc,2) |
| 105 |
xx = interp2([0 1],[0 1],xg(ix:ix+1,iy:iy+1)',xfine,xfine')'; |
| 106 |
xcf((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = xx; |
| 107 |
yy = interp2([0 1],[0 1],yg(ix:ix+1,iy:iy+1)',xfine,xfine')'; |
| 108 |
ycf((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = yy; |
| 109 |
end |
| 110 |
end |
| 111 |
|
| 112 |
xcf = [xcf xcf3]; |
| 113 |
ycf = [ycf ycf3]; |
| 114 |
|
| 115 |
save HFGRID.mat xcf ycf nfinepow |
| 116 |
|
| 117 |
disp( sprintf('Fine grid: %ix%i',size(xcf,1),size(ycf,2) )); |
| 118 |
%disp( sprintf('Fine grid res.: %gx%g',[min(dlonf(:)) min(dlatf(:))]) ); |
| 119 |
|
| 120 |
% now create the "cap" for the coarse grid |
| 121 |
xc3 = ginfo(3).XC; |
| 122 |
xc3 = [rot90(xc3);xc3;rot90(rot90(rot90(xc3)));rot90(rot90(xc3))]; |
| 123 |
yc3 = ginfo(3).YC; |
| 124 |
yc3 = [rot90(yc3);yc3;rot90(rot90(rot90(yc3)));rot90(rot90(yc3))]; |
| 125 |
xg3 = ginfo(3).XG; |
| 126 |
xg3 = [rot90(xg3(:,1:end-1));xg3(1:end-1,:);rot90(rot90(rot90(xg3(:,2:end))));rot90(rot90(xg3))]; |
| 127 |
yg3 = ginfo(3).YG; |
| 128 |
yg3 = [rot90(yg3(:,1:end-1));yg3(1:end-1,:);rot90(rot90(rot90(yg3(:,2:end))));rot90(rot90(yg3))]; |
| 129 |
zA3 = ginfo(3).rA; |
| 130 |
zA3 = [rot90(zA3);zA3;rot90(rot90(rot90(zA3)));rot90(rot90(zA3))]; |
| 131 |
|
| 132 |
i1 = find(xc3<0); xc3(i1) = xc3(i1) + 360; |
| 133 |
i1 = find(xg3<0); xg3(i1) = xg3(i1) + 360; |
| 134 |
|
| 135 |
% assemble |
| 136 |
xc = [xc xc3+s_lon]; |
| 137 |
yc = [yc yc3+s_lat]; |
| 138 |
xg = [xg xg3+s_lon]; |
| 139 |
yg = [yg yg3+s_lat]; |
| 140 |
zA = [zA zA3]; |
| 141 |
|
| 142 |
nxc=size(xc,1); |
| 143 |
disp([sprintf('nxc = %i factors(nxc)=',nxc) sprintf(' %i',factor(nxc))]); |
| 144 |
nyc=size(yc,2); |
| 145 |
disp([sprintf('nyc = %i factors(nyc)=',nyc) sprintf(' %i',factor(nyc))]); |
| 146 |
|
| 147 |
lon_lo=round(min(xg(:))); |
| 148 |
lon_L=round(max(xg(:)))-lon_lo; |
| 149 |
lat_lo=round(min(yg(:))); |
| 150 |
lat_L=round(max(yg(:)))-lat_lo; |
| 151 |
lon_hi=lon_lo+lon_L; |
| 152 |
lat_hi=lat_lo+lat_L; |
| 153 |
|
| 154 |
xu=(xg(:,1:end-1)+xg(:,2:end))/2; |
| 155 |
yu=(yg(:,1:end-1)+yg(:,2:end))/2; |
| 156 |
xv=(xg(1:end-1,:)+xg(2:end,:))/2; |
| 157 |
yv=(yg(1:end-1,:)+yg(2:end,:))/2; |
| 158 |
|
| 159 |
disp( sprintf('[lon_lo lon hi lat_lo lat_hi] = %f %f %f %f',[lon_lo lon_hi lat_lo lat_hi])); |
| 160 |
disp( sprintf('Periodic = %i',periodic) ); |
| 161 |
|
| 162 |
% Structure |
| 163 |
GRID.xc=xc; |
| 164 |
GRID.yc=yc; |
| 165 |
GRID.xg=xg; |
| 166 |
GRID.yg=yg; |
| 167 |
GRID.rac=zA; |
| 168 |
|
| 169 |
% Store data in HGRID.mat |
| 170 |
nt=1; |
| 171 |
save HGRID.mat nxc nyc nt xc yc xu yu xv yv xg yg zA ... |
| 172 |
lon_lo lon_hi lat_lo lat_hi periodic GRID |
| 173 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph