| 1 |
function [dxg,dyg,dxf,dyf,dxc,dyc,dxv,dyu,Ec,Eu,Ev,Ez,latC,lonC,latG,lonG,... |
| 2 |
Q11,Q22,Q12, TUu,TUv,TVu,TVv ] ... |
| 3 |
= gengrid_fn(nx,nratio,meth,ornt,plotfigs,writedata) |
| 4 |
|
| 5 |
% Purser algorithm |
| 6 |
% ngen=2; nrat2=1; nx=2*(1+nrat2*(2^ngen))-1; |
| 7 |
% Number nodes along edge of cube |
| 8 |
%nx=33; |
| 9 |
% Ratio of working fine grid to target grid |
| 10 |
%nratio=4; |
| 11 |
% Method |
| 12 |
%meth='conf'; |
| 13 |
%meth='q=1'; |
| 14 |
% Orientation |
| 15 |
%ornt='c'; |
| 16 |
% Number of nodes along edge of fine cube (working grid) |
| 17 |
nxf=nratio*(nx-1)+1; |
| 18 |
|
| 19 |
%plotfigs=0; |
| 20 |
%writedata=1; |
| 21 |
|
| 22 |
% Generate gridded face using pure conformal mapping |
| 23 |
%q=-1:2/(nxf-1):1; % Coordinate -1 < q < +1 spans full face |
| 24 |
q=-1:2/(nxf-1):0; % Coordinate -1 < q < 0 spans first quadrant |
| 25 |
%q=rescale_coordinate(q,'q=0'); |
| 26 |
%q=rescale_coordinate(q,'q=78'); |
| 27 |
%q=rescale_coordinate(q,'q=1'); |
| 28 |
%q=rescale_coordinate(q,'q=i3'); |
| 29 |
%q=rescale_coordinate(q,'tan'); |
| 30 |
%q=rescale_coordinate(q,'conf'); |
| 31 |
q=rescale_coordinate(q,meth); |
| 32 |
|
| 33 |
[lx,ly]=ndgrid(q,q); % Local fine grid curvilinear coordinate (lx,ly) |
| 34 |
|
| 35 |
% Calculate simple finite volume info for fine grid (lx,ly); |
| 36 |
[dx,dy,E] = calc_fvgrid(lx,ly); clear dy; |
| 37 |
|
| 38 |
% Use reflection symmetry of quadrants to fill face |
| 39 |
dx((nxf+1)/2:nxf-1,:)=dx((nxf-1)/2:-1:1,:); |
| 40 |
dx(:,(nxf+1)/2:nxf)=dx(:,(nxf+1)/2:-1:1); |
| 41 |
E((nxf+1)/2:nxf-1,:)=E((nxf-1)/2:-1:1,:); |
| 42 |
E(:,(nxf+1)/2:nxf-1)=E(:,(nxf-1)/2:-1:1); |
| 43 |
|
| 44 |
[dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio); |
| 45 |
[Ec,Ez,Ev]=reduce_E(E,nratio); |
| 46 |
clear dx E % tidy up |
| 47 |
|
| 48 |
% Remaining grid information |
| 49 |
dyg=dxg'; |
| 50 |
dyc=dxc'; |
| 51 |
dyf=dxf'; |
| 52 |
dyu=dxv'; |
| 53 |
Eu=Ev'; |
| 54 |
|
| 55 |
% Calculate geographic coordinates |
| 56 |
switch ornt |
| 57 |
case 'c' |
| 58 |
[LatG,LonG]=calc_geocoords_centerpole(lx,ly,1); |
| 59 |
for n=2:6; [LatG(:,:,n),LonG(:,:,n)]=calc_geocoords_centerpole(lx,ly,n); end; |
| 60 |
LatG=permutetiles(LatG,2); % this is to be consistent with earlier grids |
| 61 |
LonG=permutetiles(LonG,2); % this is to be consistent with earlier grids |
| 62 |
case 'v' |
| 63 |
[LatG,LonG]=calc_geocoords_cornerpole(lx,ly,1); |
| 64 |
for n=2:6; [LatG(:,:,n),LonG(:,:,n)]=calc_geocoords_cornerpole(lx,ly,n); end; |
| 65 |
otherwise |
| 66 |
ornt |
| 67 |
error('Unknown orientation'); |
| 68 |
end |
| 69 |
|
| 70 |
if nratio~=1 |
| 71 |
% Calculate metric tensor |
| 72 |
Q=q;Q(end:nxf)=-q(end:-1:1); |
| 73 |
qx=Q(2+nratio/2:nratio:end)-Q(nratio/2:nratio:end-2); |
| 74 |
[QX,QY]=ndgrid(qx,qx); clear qx Q |
| 75 |
[Xf,Yf,Zf]=map_lonlat2xyz(LonG,LatG); |
| 76 |
dXdx=( Xf(2+nratio/2:nratio:end ,1+nratio/2:nratio:end,:) ... |
| 77 |
-Xf( nratio/2:nratio:end-2,1+nratio/2:nratio:end,:) )./expand(QX); |
| 78 |
dYdx=( Yf(2+nratio/2:nratio:end ,1+nratio/2:nratio:end,:) ... |
| 79 |
-Yf( nratio/2:nratio:end-2,1+nratio/2:nratio:end,:) )./expand(QX); |
| 80 |
dZdx=( Zf(2+nratio/2:nratio:end ,1+nratio/2:nratio:end,:) ... |
| 81 |
-Zf( nratio/2:nratio:end-2,1+nratio/2:nratio:end,:) )./expand(QX); |
| 82 |
dXdy=( Xf(1+nratio/2:nratio:end,2+nratio/2:nratio:end ,:) ... |
| 83 |
-Xf(1+nratio/2:nratio:end, nratio/2:nratio:end-2,:) )./expand(QY); |
| 84 |
dYdy=( Yf(1+nratio/2:nratio:end,2+nratio/2:nratio:end ,:) ... |
| 85 |
-Yf(1+nratio/2:nratio:end, nratio/2:nratio:end-2,:) )./expand(QY); |
| 86 |
dZdy=( Zf(1+nratio/2:nratio:end,2+nratio/2:nratio:end ,:) ... |
| 87 |
-Zf(1+nratio/2:nratio:end, nratio/2:nratio:end-2,:) )./expand(QY); |
| 88 |
Q11=dXdx.*dXdx+dYdx.*dYdx+dZdx.*dZdx; |
| 89 |
Q22=dXdy.*dXdy+dYdy.*dYdy+dZdy.*dZdy; |
| 90 |
Q12=dXdx.*dXdy+dYdx.*dYdy+dZdx.*dZdy; |
| 91 |
clear Xf Yf Zf QX QY Q |
| 92 |
else |
| 93 |
Q11=0; |
| 94 |
Q12=0; |
| 95 |
Q22=0; |
| 96 |
end |
| 97 |
|
| 98 |
% Sub-sample to obtain model coordinates |
| 99 |
latG=LatG(1:nratio:end,1:nratio:end,:); |
| 100 |
lonG=LonG(1:nratio:end,1:nratio:end,:); |
| 101 |
if nratio==1 |
| 102 |
latC=( latG(1:end-1,1:end-1,:) ... |
| 103 |
+latG(2:end ,1:end-1,:) ... |
| 104 |
+latG(1:end-1,2:end ,:) ... |
| 105 |
+latG(2:end ,2:end ,:) )/4; |
| 106 |
lonC=( lonG(1:end-1,1:end-1,:) ... |
| 107 |
+lonG(2:end ,1:end-1,:) ... |
| 108 |
+lonG(1:end-1,2:end ,:) ... |
| 109 |
+lonG(2:end ,2:end ,:) )/4; |
| 110 |
else |
| 111 |
latC=LatG(1+nratio/2:nratio:end,1+nratio/2:nratio:end,:); |
| 112 |
lonC=LonG(1+nratio/2:nratio:end,1+nratio/2:nratio:end,:); |
| 113 |
end |
| 114 |
clear LatG LonG % tidy up |
| 115 |
|
| 116 |
if nratio~=1 |
| 117 |
% Flow rotation tensor |
| 118 |
Xlon=-sin(lonC); Ylon= cos(lonC); Zlon= 0*lonC; |
| 119 |
TUu= (dXdx.*Xlon+dYdx.*Ylon+dZdx.*Zlon); |
| 120 |
TVu=-(dXdy.*Xlon+dYdy.*Ylon+dZdy.*Zlon); |
| 121 |
Xlat=-sin(latC).*cos(lonC); Ylat=-sin(latC).*sin(lonC); Zlat= cos(latC); |
| 122 |
TUv=-(dXdx.*Xlat+dYdx.*Ylat+dZdx.*Zlat); |
| 123 |
TVv= (dXdy.*Xlat+dYdy.*Ylat+dZdy.*Zlat); |
| 124 |
det=sqrt(TUu.*TVv-TUv.*TVu); |
| 125 |
%det=sqrt(TUu.*TUu+TUv.*TUv); |
| 126 |
TUu=TUu./det; |
| 127 |
TUv=TUv./det; |
| 128 |
%det=sqrt(TVv.*TVv+TVu.*TVu); |
| 129 |
TVu=TVu./det; |
| 130 |
TVv=TVv./det; |
| 131 |
clear dXdx dXdy dYdx dYdy dZdx dZdy |
| 132 |
end |
| 133 |
|
| 134 |
% 3D coordinates for plotting |
| 135 |
[XG,YG,ZG]=map_lonlat2xyz(lonG,latG); |
| 136 |
|
| 137 |
% Symmetry? |
| 138 |
stats( dxg-dxg(end:-1:1,:) ) |
| 139 |
stats( dxg(:,1:end)-dxg(:,end:-1:1) ) |
| 140 |
stats( dyg-dyg(:,end:-1:1) ) |
| 141 |
stats( dyg(1:end,:)-dyg(end:-1:1,:) ) |
| 142 |
stats( dxg-dyg' ) |
| 143 |
stats( dxf-dxf(end:-1:1,:) ) |
| 144 |
stats( dxf-dxf(:,end:-1:1) ) |
| 145 |
stats( dyf-dyf(end:-1:1,:) ) |
| 146 |
stats( dxf-dyf' ) |
| 147 |
stats( Ec-Ec(end:-1:1,:) ) |
| 148 |
stats( Ec-Ec(:,end:-1:1) ) |
| 149 |
stats( Ec-Ec' ) |
| 150 |
stats( Eu-Ev' ) |
| 151 |
stats( Eu-Eu(:,end:-1:1) ) |
| 152 |
stats( Ev-Ev(end:-1:1,:) ) |
| 153 |
|
| 154 |
if writedata~=0 |
| 155 |
dir=sprintf('cs-%s-%s-%i-%i/',meth,ornt,nx-1,nxf-1); |
| 156 |
[status,msg]=mkdir(dir); if status==0; disp(msg); end; |
| 157 |
clear status msg |
| 158 |
|
| 159 |
fid=fopen([dir 'grid.info'],'w'); |
| 160 |
fprintf(fid,'nx=%i\n',nx-1); |
| 161 |
fprintf(fid,'nratio=%i\n',nratio); |
| 162 |
fprintf(fid,'nxf=%i\n',(nx-1)*nratio); |
| 163 |
fprintf(fid,'mapping="%s"\n',meth); |
| 164 |
fprintf(fid,'orientation="%s"\n',ornt); |
| 165 |
fclose(fid); |
| 166 |
|
| 167 |
Rsphere=6370e3; |
| 168 |
write_tiles([dir 'LONG'], 180/pi*lonG ) |
| 169 |
write_tiles([dir 'LATG'], 180/pi*latG ) |
| 170 |
write_tiles([dir 'LONC'], 180/pi*lonC ) |
| 171 |
write_tiles([dir 'LATC'], 180/pi*latC ) |
| 172 |
%write_tiles([dir 'DXG'], Rsphere*expand(dxg) ) |
| 173 |
%write_tiles([dir 'DYG'], Rsphere*expand(dyg) ) |
| 174 |
%write_tiles([dir 'DXF'], Rsphere*expand(dxf) ) |
| 175 |
%write_tiles([dir 'DYF'], Rsphere*expand(dyf) ) |
| 176 |
%write_tiles([dir 'DXC'], Rsphere*expand(dxc) ) |
| 177 |
%write_tiles([dir 'DYC'], Rsphere*expand(dyc) ) |
| 178 |
%write_tiles([dir 'DXV'], Rsphere*expand(dxv) ) |
| 179 |
%write_tiles([dir 'DYU'], Rsphere*expand(dyu) ) |
| 180 |
%write_tiles([dir 'RA'], Rsphere^2*expand(Ec) ) |
| 181 |
%write_tiles([dir 'RAW'], Rsphere^2*expand(Eu) ) |
| 182 |
%write_tiles([dir 'RAS'], Rsphere^2*expand(Ev) ) |
| 183 |
%write_tiles([dir 'RAZ'], Rsphere^2*expand(Ez) ) |
| 184 |
write_tile([dir 'DXG'], Rsphere*(dxg) ) |
| 185 |
write_tile([dir 'DYG'], Rsphere*(dyg) ) |
| 186 |
write_tile([dir 'DXF'], Rsphere*(dxf) ) |
| 187 |
write_tile([dir 'DYF'], Rsphere*(dyf) ) |
| 188 |
write_tile([dir 'DXC'], Rsphere*(dxc) ) |
| 189 |
write_tile([dir 'DYC'], Rsphere*(dyc) ) |
| 190 |
write_tile([dir 'DXV'], Rsphere*(dxv) ) |
| 191 |
write_tile([dir 'DYU'], Rsphere*(dyu) ) |
| 192 |
write_tile([dir 'RA'], Rsphere^2*(Ec) ) |
| 193 |
write_tile([dir 'RAW'], Rsphere^2*(Eu) ) |
| 194 |
write_tile([dir 'RAS'], Rsphere^2*(Ev) ) |
| 195 |
write_tile([dir 'RAZ'], Rsphere^2*(Ez) ) |
| 196 |
|
| 197 |
% Save stuff for pre-processing/post-processing |
| 198 |
save([dir 'CUBE_HGRID.mat'],'latG','lonG','latC','lonC','dxg','dyg','dxc','dyc','Ec'); |
| 199 |
save([dir 'CUBE_3DGRID.mat'],'XG','YG','ZG'); |
| 200 |
save([dir 'TUV.mat'],'TUu','TUv','TVu','TVv'); |
| 201 |
end |
| 202 |
|
| 203 |
if plotfigs~=0 |
| 204 |
% Some basic pictures |
| 205 |
figure(1); |
| 206 |
subplot(221);imagesc(dxv');axis xy;title('\Delta x_v');colorbar |
| 207 |
subplot(222);imagesc(dxg');axis xy;title('\Delta x_g');colorbar |
| 208 |
subplot(223);imagesc(dxc');axis xy;title('\Delta x_c');colorbar |
| 209 |
subplot(224);imagesc(dxf');axis xy;title('\Delta x_f');colorbar |
| 210 |
|
| 211 |
figure(2); |
| 212 |
subplot(221);imagesc( ((Ec-dxf.*dxf')./Ec)' );axis xy;title('|E_c-\Delta x_f\Delta y_f|');colorbar |
| 213 |
subplot(222);imagesc(Ec');axis xy;title('E_c');colorbar |
| 214 |
subplot(223);imagesc(Ez');axis xy;title('E_z');colorbar |
| 215 |
subplot(224);imagesc(Ev');axis xy;title('E_v');colorbar |
| 216 |
|
| 217 |
% Plot rotation tensor |
| 218 |
figure(3) |
| 219 |
subplot(221) |
| 220 |
plotcube(XG,YG,ZG,TUu);colorbar;%[az el]=view;view([az+180 el]) |
| 221 |
title('TUu') |
| 222 |
subplot(222) |
| 223 |
plotcube(XG,YG,ZG,TUv);colorbar;%[az el]=view;view([az+180 el]) |
| 224 |
title('TUv') |
| 225 |
subplot(223) |
| 226 |
plotcube(XG,YG,ZG,TVu);colorbar;%[az el]=view;view([az+180 el]) |
| 227 |
title('TVu') |
| 228 |
subplot(224) |
| 229 |
plotcube(XG,YG,ZG,TVv);colorbar;%[az el]=view;view([az+180 el]) |
| 230 |
title('TVv') |
| 231 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph