| 1 |
mlosch |
1.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