gael/bulkMatlab/gpcp_load_field.m

function [field_out]=gpcp_load_fields(ycur,hcur);
%loads the fields and does the interpolation to the 1x1 ECCO grid
domaine_global_def;

mask=squeeze(tmask3D(:,:,1));

rep_in='/raid0/gforget/DATAbin/forcing/GPCP/';

jpig=144; jpjg=72;
lon2D_g=[1.25:2.5:358.75]'*ones(1,jpjg); lat2D_g=ones(jpig,1)*[-88.75:2.5:88.75];
mask_g=ones(jpig,jpjg); recl=jpig*jpjg*4;

[x1,y1] = meshgrid([lon2D_g(end,1)-360  lon2D_g(:,1)' lon2D_g(1:2,1)'+360],lat2D_g(1,:));
[x2,y2] = meshgrid(lon2D_t(:,1),lat2D_t(1,:));


%interpolate from the monthly field:
%       coeffs etc:
dcur=ceil((hcur)/4); [mcur,mcurW]=coeff_MonthlyAtlasInterp(dcur);
ycur2=[ycur ycur];
if ~isempty(find(mcur==1))&dcur>300; ycur2(find(mcur==1))=ycur+1; end;
if ~isempty(find(mcur==12))&dcur<100; ycur2(find(mcur==12))=ycur-1; end;
tmp1=find(ycur2==1991); ycur2(tmp1)=1992; mcur(tmp1)=1;
tmp1=find(ycur2==2005); ycur2(tmp1)=2004; mcur(tmp1)=12;
%       load the two fields:
fid=fopen([rep_in 'gpcp_v2_psg.' num2str(ycur2(1))],'r','b');
status=fseek(fid,576,'bof'); %skip header
position0=recl*(mcur(1)-1); status=fseek(fid,position0,'bof');
field1=flipdim(fread(fid,[jpig jpjg],'float32'),2); fclose(fid);
fid=fopen([rep_in 'gpcp_v2_psg.' num2str(ycur2(2))],'r','b');
status=fseek(fid,576,'bof'); %skip header
position0=recl*(mcur(2)-1); status=fseek(fid,position0,'bof');
field2=flipdim(fread(fid,[jpig jpjg],'float32'),2); fclose(fid);
field_cur=(1-mcurW)*field1+mcurW*field2;

%interpolate:
tmp1=[field_cur(end,:);field_cur;field_cur(1,:);field_cur(2,:)]';
field_int=interp2(x1,y1,tmp1,x2,y2)';

field_out=field_int.*mask/86400;


1	gforget	1.1	function [field_out]=gpcp_load_fields(ycur,hcur);
2			%loads the fields and does the interpolation to the 1x1 ECCO grid
3			domaine_global_def;
4
5			mask=squeeze(tmask3D(:,:,1));
6
7			rep_in='/raid0/gforget/DATAbin/forcing/GPCP/';
8
9			jpig=144; jpjg=72;
10			lon2D_g=[1.25:2.5:358.75]'ones(1,jpjg); lat2D_g=ones(jpig,1)[-88.75:2.5:88.75];
11			mask_g=ones(jpig,jpjg); recl=jpigjpjg4;
12
13			[x1,y1] = meshgrid([lon2D_g(end,1)-360 lon2D_g(:,1)' lon2D_g(1:2,1)'+360],lat2D_g(1,:));
14			[x2,y2] = meshgrid(lon2D_t(:,1),lat2D_t(1,:));
15
16
17			%interpolate from the monthly field:
18			% coeffs etc:
19			dcur=ceil((hcur)/4); [mcur,mcurW]=coeff_MonthlyAtlasInterp(dcur);
20			ycur2=[ycur ycur];
21			if ~isempty(find(mcur==1))&dcur>300; ycur2(find(mcur==1))=ycur+1; end;
22			if ~isempty(find(mcur==12))&dcur<100; ycur2(find(mcur==12))=ycur-1; end;
23			tmp1=find(ycur2==1991); ycur2(tmp1)=1992; mcur(tmp1)=1;
24			tmp1=find(ycur2==2005); ycur2(tmp1)=2004; mcur(tmp1)=12;
25			% load the two fields:
26			fid=fopen([rep_in 'gpcp_v2_psg.' num2str(ycur2(1))],'r','b');
27			status=fseek(fid,576,'bof'); %skip header
28			position0=recl*(mcur(1)-1); status=fseek(fid,position0,'bof');
29			field1=flipdim(fread(fid,[jpig jpjg],'float32'),2); fclose(fid);
30			fid=fopen([rep_in 'gpcp_v2_psg.' num2str(ycur2(2))],'r','b');
31			status=fseek(fid,576,'bof'); %skip header
32			position0=recl*(mcur(2)-1); status=fseek(fid,position0,'bof');
33			field2=flipdim(fread(fid,[jpig jpjg],'float32'),2); fclose(fid);
34			field_cur=(1-mcurW)field1+mcurWfield2;
35
36			%interpolate:
37			tmp1=[field_cur(end,:);field_cur;field_cur(1,:);field_cur(2,:)]';
38			field_int=interp2(x1,y1,tmp1,x2,y2)';
39
40			field_out=field_int.*mask/86400;
41
42