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