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;