Warning, /utils/matlab/cs_grid/check_gridfiles.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
0b707f8415 Jean*0001
0002 %- for a symetric cubed-sphere grid (nr=ng=nb), check that:
0003 % 1) all the faces have identical grid-spacing & grid-cell area
0004 % + long & lat match exactly.
b86dbee8a5 Jean*0005 % 2) grid-spacing & grid-cell area (+ long & lat) are symetric
0b707f8415 Jean*0006 % regarding the 3 symetry axes of a square.
0007
b86dbee8a5 Jean*0008 nc=32;
0009 %nc=96;
0b707f8415 Jean*0010
b86dbee8a5 Jean*0011 rDir='./';
0012 inpName='tile.mitgrid';
0013 %rDir='../../cs96_dxC3_dXYa/';
0014 %inpName='cs96_dxC3_dXYa';
0b707f8415 Jean*0015 strict=0; %- set to 1 to only check interior [nc,nc] if cell center var.
0016
0017 np=nc+1; nh=nc/2;
0018
0019 for n=1:6,
0020 %-- read :
b86dbee8a5 Jean*0021 p=strfind(inpName,'.mitgrid');
0b707f8415 Jean*0022 if isempty(p),
0023 namF=sprintf([rDir,inpName,'.face%3.3i.bin'],n);
0024 else
0025 namF=sprintf([rDir,inpName(1:4),'%3.3i',inpName(5:end)],n);
0026 end
0027 fid=fopen(namF,'r','b');
0028 vv1=fread(fid,'real*8');
0029 fclose(fid);
0030 s=size(vv1,1); k1=s/np/np;
0031 fprintf(['read: ',namF,' : size: %i (%ix%ix%i)\n'],s,np,np,k1);
0032 vv1=reshape(vv1,[np np k1]);
b86dbee8a5 Jean*0033 %--- save each face:
0b707f8415 Jean*0034 if n == 1,
b86dbee8a5 Jean*0035 vvf=zeros(np,np,6,k1);
0b707f8415 Jean*0036 end
b86dbee8a5 Jean*0037 vvf(:,:,n,:)=vv1;
0b707f8415 Jean*0038 end
0039
0040 %- XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG:
0041 % 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0042 gpos=[3 3 3 3 3 0 0 0 0 0 1 2 1 2 2 1];
0043 % gpos: 0 = grid-cell corner ; 1 = U point ; 2 = V point ; 3 = Tracer point
0044
b86dbee8a5 Jean*0045 namV={'xC','yC','dxF','dyF','rAc','xG','yG','dxV','dyU','rAz', ...
0046 'dxC','dyC','rAw','rAs','dxG','dyG','AngleCS','AngleSN'};
0047
0048 %-------------------------------------------------------------------------------
0049
0050 %-- Check "mising corner" values (for grid-cell corner variables)
0051 % m1: North-West corner (1,nc+1) on face f1,f3,f5 ;
0052 % m2: South-East corner (nc+1,1) on face f2,f4,f5 ;
0053 % Note: Since these value are not exchanged, they need to be
0054 % identical in the 3 grid-files.
0055 fprintf('-- check (missing) corner values :\n');
0056
0057 for k=1:min(k1,16),
0058 if gpos(k) == 0,
0059 for m=1:2,
0060 if m == 1, tmp=vvf(1,np,[m:2:6],k);
0061 else tmp=vvf(np,1,[m:2:6],k); end
0062 dd=tmp(2:3)-tmp(1);
0063 mxd=max(abs(dd));
0064 if mxd > 0.,
0065 fprintf(' %3s k= %2i : Corner Max diff m%i = %g (f%i,f%i= %g %g)\n', ...
0066 char(namV(k)),k,m,mxd,m+[2 4],dd(1),dd(2))
0067 end
0068 end
0069 end
0070 end
0071
0072 %-------------------------------------------------------------------------------
0073
0074 fprintf('-- check longitude & latitude :\n');
0075
0076 for k=1:min(k1,16),
0b707f8415 Jean*0077 var=vvf(:,:,1,k); avr=mean(var(:));
0078 if k < 3 | k == 6 | k == 7,
0079 %-- check long & latitude:
0080 epsil=3.e-14;
0081 list=[2 4 5];
0082 for n=list,
0083 mxd=0; shift=0;
0084 if k == 1 | k == 6,
0085 shift=90;
0086 if n == 5, shift=-90; end
0087 if n == 4, shift=-180; end
0088 end
0089 if n == 2;
0090 dd=vvf(:,:,n,k)-vvf(:,:,1,k);
0091 if k < 3, dd=dd(1:nc,1:nc); end
0092 else
0093 if gpos(k) == 3,
0094 var=vvf([1:nc],[1:nc],1,k);
0095 dd=vvf([nc:-1:1],[1:nc],n,k);
0096 dd=permute(dd,[2 1])-var;
0097 else
0098 var=vvf(:,:,1,k);
0099 dd=vvf(np:-1:1,:,n,k);
0100 dd=permute(dd,[2 1])-var;
0101 % if n == 4, dd(nh+1,:)=0; end %- skip line of day change
0102 end
0103 end
0104 dd=dd-shift;
0105 mnD=min(dd(:)); MxD=max(dd(:));
0106 ss=360*round(dd/360); dd=dd-ss;
0107 dd=abs(dd);
0108 mxd=max(dd(:));
b86dbee8a5 Jean*0109 if mxd > epsil,
0110 fprintf(' %3s k,n= %2i %2i : Max diff f1 = %g (mn,Mx= %g %g)\n', ...
0111 char(namV(k)),k,n,mxd,mnD,MxD)
0b707f8415 Jean*0112 end
0113 end
0114 n=6;
0115 if gpos(k) == 3,
0116 var=vvf([1:nc],[1:nc],3,k);
0117 tmp=vvf([1:nc],[1:nc],n,k);
0118 else
0119 var=vvf(:,:,3,k);
0120 tmp=vvf(:,:,n,k);
0121 end
0122 dd=tmp(end:-1:1,end:-1:1); dd=permute(dd,[2 1]);
0123 if k == 2 | k == 7,
0124 dd=dd+var;
0125 else
0126 dd=dd-var;
0127 if gpos(k) == 0, dd(nh+1,nh+1)=0; end %- skip the pole longitude
0128 end
0129 mnD=min(dd(:)); MxD=max(dd(:));
0130 ss=360*round(dd/360); dd=dd-ss;
0131 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0132 if mxd > epsil,
0133 fprintf(' %3s k,n= %2i %2i : Max diff f3 = %g (mn,Mx= %g %g)\n', ...
0134 char(namV(k)),k,n,mxd,mnD,MxD)
0b707f8415 Jean*0135 end
0136 %-- check for symetry versus x=Xmid median on face 1
0137 mxd=0;
0138 if gpos(k) == 3,
0139 var=vvf([1:nc],[1:nc],1,k);
0140 sym=var([nc:-1:1],[1:nc]);
0141 else
0142 var=vvf(:,:,1,k);
0143 sym=var([np:-1:1],:);
0144 end
0145 if k == 1 | k == 6,
b86dbee8a5 Jean*0146 dd=sym+var;
0b707f8415 Jean*0147 else
b86dbee8a5 Jean*0148 dd=sym-var;
0b707f8415 Jean*0149 end
0150 mnD=min(dd(:)); MxD=max(dd(:));
0151 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0152 if mxd > epsil,
0153 fprintf(' %3s k= %2i : Sym1.1 Max diff= %g (mn,Mx= %g %g)\n', ...
0154 char(namV(k)),k,mxd,mnD,MxD)
0b707f8415 Jean*0155 end
0156 %-- check for symetry versus y=Ymid median on face 1
0157 mxd=0;
0158 if gpos(k) == 3,
0159 var=vvf([1:nc],[1:nc],1,k);
0160 sym=var([1:nc],[nc:-1:1]);
0161 else
0162 var=vvf(:,:,1,k);
0163 sym=var(:,[np:-1:1]);
0164 end
0165 if k == 1 | k == 6,
b86dbee8a5 Jean*0166 dd=sym-var;
0b707f8415 Jean*0167 else
b86dbee8a5 Jean*0168 dd=sym+var;
0b707f8415 Jean*0169 end
0170 mnD=min(dd(:)); MxD=max(dd(:));
0171 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0172 if mxd > epsil,
0173 fprintf(' %3s k= %2i : Sym1.2 Max diff= %g (mn,Mx= %g %g)\n', ...
0174 char(namV(k)),k,mxd,mnD,MxD)
0b707f8415 Jean*0175 end
0176 %-- check for symetry versus x=Xmid median on face 3
0177 % check only half of the face (enough since check Sym.2 on full face)
0178 mxd=0;
0179 if gpos(k) == 3,
0180 var=vvf([1:nc],[1:nh],3,k);
0181 sym=var([nc:-1:1],[1:nh]);
0182 else
0183 %var=vvf(:,[1:nh+1],3,k);
0184 var=vvf(:,[1:nh],3,k); %- skip the line of day change
0185 sym=var([np:-1:1],:);
0186 end
0187 if k == 1 | k == 6,
0188 dd=sym+var; %- X
0189 dd=dd-180;
0190 else
0191 dd=sym-var; %- Y
0192 end
0193 mnD=min(dd(:)); MxD=max(dd(:));
0194 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0195 if mxd > epsil,
0196 fprintf(' %3s k= %2i : Sym3.1 Max diff= %g (mn,Mx= %g %g)\n', ...
0197 char(namV(k)),k,mxd,mnD,MxD)
0b707f8415 Jean*0198 end
0199 %-- check for symetry versus y=Ymid median on face 3
0200 mxd=0;
0201 if gpos(k) == 3,
0202 var=vvf([1:nc],[1:nc],3,k);
0203 sym=var(:,[nc:-1:1]);
0204 else
0205 var=vvf(:,:,3,k);
0206 sym=var(:,[np:-1:1]);
0207 end
0208 if k == 1 | k == 6,
0209 dd=sym+var; %- X
0210 if gpos(k) == 0, dd(nh:np,nh+1)=0; end
0211 else
0212 dd=sym-var; %- Y
0213 end
0214 mnD=min(dd(:)); MxD=max(dd(:));
0215 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0216 if mxd > epsil,
0217 fprintf(' %3s k= %2i : Sym3.2 Max diff= %g (mn,Mx= %g %g)\n', ...
0218 char(namV(k)),k,mxd,mnD,MxD)
0b707f8415 Jean*0219 end
0220 %-- check for symetry versus y=x diagonal on face 3
0221 % check only 1/4 of the face (enough since already check Sym.1 & Sym.2)
0222 mxd=0;
0223 if gpos(k) == 3,
0224 var=vvf([1:nh],[1:nh],3,k);
0225 sym=permute(var,[2 1]);
0226 else
0227 var=vvf([1:nh+1],[1:nh+1],3,k);
0228 %var=vvf([1:nh],[1:nh],3,k);
0229 sym=permute(var,[2 1]);
0230 end
0231 if k == 1 | k == 6,
0232 dd=sym+var; %- X
0233 dd=dd-90;
0234 if gpos(k) == 0, dd(nh+1,nh+1)=0; end %- skip the pole
0235 else
0236 dd=sym-var; %- Y
0237 end
0238 mnD=min(dd(:)); MxD=max(dd(:));
0239 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0240 if mxd > epsil,
0241 fprintf(' %3s k= %2i : Sym3.3 Max diff= %g (mn,Mx= %g %g)\n', ...
0242 char(namV(k)),k,mxd,mnD,MxD)
0b707f8415 Jean*0243 end
0244 end
b86dbee8a5 Jean*0245 end
0246
0247 %-------------------------------------------------------------------------------
0248
0249 fprintf('-- check grid-spacing & area :\n');
0250
0251 %- do not check the angle (even if in the files, k=17,18)
0252 for k=1:min(k1,16),
0253 var=vvf(:,:,1,k); avr=mean(var(:));
0b707f8415 Jean*0254 %-- check grid-spacing & area:
0255 if k > 2 & k ~= 6 & k ~= 7,
0256 %-- check for identical values across faces
0257 for n=2:6,
0258 dd=vvf(:,:,n,k)-vvf(:,:,1,k);
0259 if strict == 1,
0260 if gpos(k) == 1, dd=dd(:,1:nc); end
0261 if gpos(k) == 2, dd=dd(1:nc,:); end
0262 if gpos(k) == 3, dd=dd(1:nc,1:nc); end
0263 end
0264 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0265 if mxd > 0,
0266 fprintf(' %3s k,n= %2i %2i : Max diff f1 = %g (av= %g)\n', ...
0267 char(namV(k)),k,n,mxd,avr)
0b707f8415 Jean*0268 end
0269 end
0270 %-- check for symetry versus x=Xmid median (only on face 1)
0271 msd=0;
0272 if gpos(k) >= 2,
0273 var=vvf([1:nc],:,1,k);
0274 sym=var([nc:-1:1],:);
0275 else
0276 var=vvf(:,:,1,k);
0277 sym=var([np:-1:1],:);
0278 end
0279 dd=sym-var;
0280 if strict == 1 & rem(gpos(k),2) == 1, dd=dd(:,1:nc); end
0281 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0282 if mxd > 0,
0283 fprintf(' %3s k= %2i : Sym.1 Max diff= %g (av= %g)\n',char(namV(k)),k,mxd,avr)
0b707f8415 Jean*0284 end
0285 %-- check for symetry versus y=Ymid median (only on face 1)
0286 msd=0;
0287 if rem(gpos(k),2) == 1,
0288 var=vvf(:,[1:nc],1,k);
0289 sym=var(:,[nc:-1:1]);
0290 else
0291 var=vvf(:,:,1,k);
0292 sym=var(:,[np:-1:1]);
0293 end
b86dbee8a5 Jean*0294 dd=sym-var;
0b707f8415 Jean*0295 if strict == 1 & gpos(k) >=2, dd=dd(1:nc,:); end
0296 dd=abs(dd); mxd=max(dd(:));
b86dbee8a5 Jean*0297 if mxd > 0,
0298 fprintf(' %3s k= %2i : Sym.2 Max diff= %g (av= %g)\n',char(namV(k)),k,mxd,avr)
0b707f8415 Jean*0299 end
0300 %-- check for symetry versus x=y diagnonal (only on face 1)
0301 msd=0;
0302 if k == 5,
0303 var=vvf(1:nc,1:nc,1,k);
0304 sym=permute(var,[2 1]);
0305 dd=sym-var; dd=abs(dd);
0306 mxd=max(dd(:));
0307 end
0308 if k == 3,
0309 var=vvf(1:nc,1:nc,1,k);
0310 sym=vvf(1:nc,1:nc,1,k+1);
0311 sym=permute(sym,[2 1]);
0312 dd=sym-var; dd=abs(dd);
0313 mxd=max(dd(:));
0314 end
0315 if k == 10,
0316 var=vvf(:,:,1,k);
0317 sym=permute(var,[2 1]);
0318 dd=sym-var; dd=abs(dd);
0319 mxd=max(dd(:));
0320 end
0321 if k == 8,
0322 var=vvf(:,:,1,k);
0323 sym=vvf(:,:,1,k+1);
0324 sym=permute(sym,[2 1]);
0325 dd=sym-var; dd=abs(dd);
0326 mxd=max(dd(:));
0327 end
0328 if k == 11 | k == 13,
0329 var=vvf(:,1:nc,1,k);
0330 sym=vvf(1:nc,:,1,k+1);
0331 sym=permute(sym,[2 1]);
0332 dd=sym-var; dd=abs(dd);
0333 mxd=max(dd(:));
0334 end
0335 if k == 16,
0336 var=vvf(:,1:nc,1,k);
0337 sym=vvf(1:nc,:,1,k-1);
0338 sym=permute(sym,[2 1]);
0339 dd=sym-var; dd=abs(dd);
0340 mxd=max(dd(:));
0341 end
b86dbee8a5 Jean*0342 if mxd > 0,
0343 fprintf(' %3s k= %2i : Sym.3 Max diff= %g (av= %g)\n',char(namV(k)),k,mxd,avr)
0b707f8415 Jean*0344 end
0345 end
0346 end
b86dbee8a5 Jean*0347 %-------------------------------------------------------------------------------
0348
0349 % size of storage array vvf : (np,np,6,k1)
0b707f8415 Jean*0350 %- XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG:
0351 % 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0352 %pos=[0 0 0 0 0 3 3 3 3 3 1 2 1 2 2 1];
0353
0354 return
0355 shift=-1; ccB=[0 0]; cbV=1; AxBx=[-180 180 -90 90]; kEnv=0;
0356 AxBx2=[30 60 25 50];
0357
b86dbee8a5 Jean*0358 n6x=nc*6;
0359 xc0=vvf(1:nc,1:nc,:,1); xc0=reshape(permute(xc0,[1 3 2]),[n6x nc]);
0360 yc0=vvf(1:nc,1:nc,:,2); yc0=reshape(permute(yc0,[1 3 2]),[n6x nc]);
0361 xg0=vvf(1:nc,1:nc,:,6); xg0=reshape(permute(xg0,[1 3 2]),[n6x nc]);
0362 yg0=vvf(1:nc,1:nc,:,7); yg0=reshape(permute(yg0,[1 3 2]),[n6x nc]);
0b707f8415 Jean*0363
0364 figure(2);clf;
0365 xtxt=0;ytxt=-98;
0366 %subplot(311);
b86dbee8a5 Jean*0367 k=5;
0368 var=vvf(1:nc,1:nc,:,k); var=reshape(permute(var,[1 3 2]),[n6x nc]);
0b707f8415 Jean*0369 grph_CS(var,xc0,yc0,xg0,yg0,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
0370 mnV=min(var(:)); MxV=max(var(:));
0371 fprintf('min,Max= %9.5g , %9.5g\n',mnV,MxV);
b86dbee8a5 Jean*0372 text(xtxt,ytxt,sprintf('%s : min,Max= %9.5g , %9.5g',char(namV(k)),mnV,MxV));
0b707f8415 Jean*0373 %title([titexp,'rAc [m^2]'])
