matlab_class/ecco_v4/map_grace.m


xC1=[0.5:359.5]'*ones(1,180);
yC1=ones(360,1)*[-89.5:89.5];
v4_reclen=convert2gcmfaces(mygrid.RAC); v4_reclen=size(v4_reclen); v4_reclen=v4_reclen(1)*v4_reclen(2);

dir_out='./'; %output directory
dir_in='./'; %input directory
file_pref='GRACE_CSR_withland_';

%load error field
[fieldErr_nosmooth,fieldErr]=read_grace_fld([dir_in 'GRACE_CSR_Error.asc']);

%process grace fields
for yy=1992:2012;
for mm=1:12;
tt=sprintf('%4i%02i',yy,mm);

file0=dir([dir_in file_pref tt '.asc']);
if ~isempty(file0);

%read
file0=file0.name;
field0=read_grace_fld([dir_in file0]);

%map
x=[xC1-360;xC1]; y=[yC1;yC1]; z=[field0;field0];
x=[x x x]; y=[y-180 y y+180]; z=[flipdim(z,2) z flipdim(z,2)];

field1=0*mygrid.XC;
for ii=1:5;
xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
zi = interp2(x',y',z',xi,yi);
%zi = griddata(x,y,z,xi,yi);
field1{ii}=zi;
end;

%mask the model poles
tmp1=find(mygrid.RAC<8e8&mygrid.YC>0); field1(tmp1)=NaN;
tmp1=find(mygrid.RAC<2e8&mygrid.YC<0); field1(tmp1)=NaN;

else;
fprintf(['did not find : ' file_pref tt '.asc \n']);
field1=NaN*mygrid.RAC;
end;

%use -999 for mask
field1(isnan(field1))=-999;

%write to file
if mm==1; fid=fopen([dir_out file_pref num2str(yy)],'w','b'); end;
fwrite(fid,convert2gcmfaces(field1),'float32');
if mm==12; fclose(fid); end;

end;
end;

%process error estimate
z=[fieldErr;fieldErr]; z=[flipdim(z,2) z flipdim(z,2)];

field1=0*mygrid.XC;
for ii=1:5;
xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
zi = interp2(x',y',z',xi,yi);
%zi = griddata(x,y,z,xi,yi);
field1{ii}=zi;
end;

%mask the model poles
tmp1=find(mygrid.RAC<8e8&mygrid.YC>0); field1(tmp1)=NaN;
tmp1=find(mygrid.RAC<2e8&mygrid.YC<0); field1(tmp1)=NaN;

%use 0 for mask
field1(isnan(field1))=0;

%write to file
write2file([dir_out file_pref 'err'],convert2gcmfaces(field1));

1	gforget	1.1
2			xC1=[0.5:359.5]'*ones(1,180);
3			yC1=ones(360,1)*[-89.5:89.5];
4			v4_reclen=convert2gcmfaces(mygrid.RAC); v4_reclen=size(v4_reclen); v4_reclen=v4_reclen(1)*v4_reclen(2);
5
6			dir_out='./'; %output directory
7			dir_in='./'; %input directory
8			file_pref='GRACE_CSR_withland_';
9
10			%load error field
11			[fieldErr_nosmooth,fieldErr]=read_grace_fld([dir_in 'GRACE_CSR_Error.asc']);
12
13			%process grace fields
14			for yy=1992:2012;
15			for mm=1:12;
16			tt=sprintf('%4i%02i',yy,mm);
17
18			file0=dir([dir_in file_pref tt '.asc']);
19			if ~isempty(file0);
20
21			%read
22			file0=file0.name;
23			field0=read_grace_fld([dir_in file0]);
24
25			%map
26			x=[xC1-360;xC1]; y=[yC1;yC1]; z=[field0;field0];
27			x=[x x x]; y=[y-180 y y+180]; z=[flipdim(z,2) z flipdim(z,2)];
28
29			field1=0*mygrid.XC;
30			for ii=1:5;
31			xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
32			zi = interp2(x',y',z',xi,yi);
33			%zi = griddata(x,y,z,xi,yi);
34			field1{ii}=zi;
35			end;
36
37			%mask the model poles
38			tmp1=find(mygrid.RAC<8e8&mygrid.YC>0); field1(tmp1)=NaN;
39			tmp1=find(mygrid.RAC<2e8&mygrid.YC<0); field1(tmp1)=NaN;
40
41			else;
42			fprintf(['did not find : ' file_pref tt '.asc \n']);
43			field1=NaN*mygrid.RAC;
44			end;
45
46			%use -999 for mask
47			field1(isnan(field1))=-999;
48
49			%write to file
50			if mm==1; fid=fopen([dir_out file_pref num2str(yy)],'w','b'); end;
51			fwrite(fid,convert2gcmfaces(field1),'float32');
52			if mm==12; fclose(fid); end;
53
54			end;
55			end;
56
57			%process error estimate
58			z=[fieldErr;fieldErr]; z=[flipdim(z,2) z flipdim(z,2)];
59
60			field1=0*mygrid.XC;
61			for ii=1:5;
62			xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
63			zi = interp2(x',y',z',xi,yi);
64			%zi = griddata(x,y,z,xi,yi);
65			field1{ii}=zi;
66			end;
67
68			%mask the model poles
69			tmp1=find(mygrid.RAC<8e8&mygrid.YC>0); field1(tmp1)=NaN;
70			tmp1=find(mygrid.RAC<2e8&mygrid.YC<0); field1(tmp1)=NaN;
71
72			%use 0 for mask
73			field1(isnan(field1))=0;
74
75			%write to file
76			write2file([dir_out file_pref 'err'],convert2gcmfaces(field1));
77