gael/bulkMatlab/ecmwf_load_fluxes.m

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

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

if ycur<2002
rep_in='/net/altix3700/raid4/king/data_1x1_92-03/forcing_ECMWF/';
files_ecmwf=strvcat('ERA40g_ustr_r1_','ERA40g_vstr_r1_','ERA40g_emp_r1_','ERA40g_turbflux_r1_','ERA40g_sw_r1_');
var_ecco=strvcat('ustress','vstress','emp','turbflux','swflux');
else
rep_in='/net/altix3700/raid4/king/data_1x1_92-03/forcing_ECMWF/';
files_ecmwf=strvcat('EAG_ustr_r1_','EAG_vstr_r1_','EAG_emp_r1_','EAG_turbflux_r1_','EAG_sw_r1_');
var_ecco=strvcat('ustress','vstress','emp','turbflux','swflux');
end

ecmwf_grid; jpig=length(x); jpjg=length(y);
lon2D_g=x'*ones(1,jpjg); lat2D_g=ones(jpig,1)*y; 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,:));

for fcur=1:size(files_ecmwf,1)

fid=fopen([rep_in deblank(files_ecmwf(fcur,:)) num2str(ycur)],'r','b');
recl=jpig*jpjg*4; position0=recl*(hcur-1); status=fseek(fid,position0,'bof');
field_mod=fread(fid,[jpig jpjg],'float32'); fclose(fid);
%field_mod=fread(fid,jpig*jpjg*4,'float32'); fclose(fid); field_mod=squeeze(mean(reshape(field_mod,[jpig jpjg 4]),3));

%there are a few NaNs in the files: fix it here by extrapolation
tmp1=find(isnan(field_mod));
while ~isempty(tmp1);
tmp2=field_mod; tmp2(tmp1)=0; tmp3=1*(~isnan(field_mod));
tmp2=tmp2+circshift(tmp2,[0 1])+circshift(tmp2,[0 -1])+circshift(tmp2,[-1 0])+circshift(tmp2,[1 0]);
tmp3=tmp3+circshift(tmp3,[0 1])+circshift(tmp3,[0 -1])+circshift(tmp3,[-1 0])+circshift(tmp3,[1 0]);
tmp1=find(isnan(field_mod)&tmp3>0); 
field_mod(tmp1)=tmp2(tmp1)./tmp3(tmp1); tmp1=find(isnan(field_mod));
end

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

eval([deblank(var_ecco(fcur,:)) '=field_int.*mask;']);

end

rhoConstFresh=999.8;
emp=emp*rhoConstFresh; %m/s->kg/m^2/s
ustress=-ustress.*mask;
vstress=-vstress.*mask;

turbflux=turbflux.*mask;
emp=emp.*mask;
swflux=swflux.*mask;

mystruct_out=struct('ustress',ustress,'vstress',vstress,'emp',emp,'turbflux',turbflux,'swflux',swflux);

1	gforget	1.1	function [mystruct_out]=ecmwf_load_fluxes(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			if ycur<2002
8			rep_in='/net/altix3700/raid4/king/data_1x1_92-03/forcing_ECMWF/';
9			files_ecmwf=strvcat('ERA40g_ustr_r1_','ERA40g_vstr_r1_','ERA40g_emp_r1_','ERA40g_turbflux_r1_','ERA40g_sw_r1_');
10			var_ecco=strvcat('ustress','vstress','emp','turbflux','swflux');
11			else
12			rep_in='/net/altix3700/raid4/king/data_1x1_92-03/forcing_ECMWF/';
13			files_ecmwf=strvcat('EAG_ustr_r1_','EAG_vstr_r1_','EAG_emp_r1_','EAG_turbflux_r1_','EAG_sw_r1_');
14			var_ecco=strvcat('ustress','vstress','emp','turbflux','swflux');
15			end
16
17			ecmwf_grid; jpig=length(x); jpjg=length(y);
18			lon2D_g=x'ones(1,jpjg); lat2D_g=ones(jpig,1)y; mask_g=ones(jpig,jpjg);
19			recl=jpigjpjg4;
20
21			[x1,y1] = meshgrid([lon2D_g(end,1)-360 lon2D_g(:,1)' lon2D_g(1:2,1)'+360],lat2D_g(1,:));
22			[x2,y2] = meshgrid(lon2D_t(:,1),lat2D_t(1,:));
23
24			for fcur=1:size(files_ecmwf,1)
25
26			fid=fopen([rep_in deblank(files_ecmwf(fcur,:)) num2str(ycur)],'r','b');
27			recl=jpigjpjg4; position0=recl*(hcur-1); status=fseek(fid,position0,'bof');
28			field_mod=fread(fid,[jpig jpjg],'float32'); fclose(fid);
29			%field_mod=fread(fid,jpigjpjg4,'float32'); fclose(fid); field_mod=squeeze(mean(reshape(field_mod,[jpig jpjg 4]),3));
30
31			%there are a few NaNs in the files: fix it here by extrapolation
32			tmp1=find(isnan(field_mod));
33			while ~isempty(tmp1);
34			tmp2=field_mod; tmp2(tmp1)=0; tmp3=1*(~isnan(field_mod));
35			tmp2=tmp2+circshift(tmp2,[0 1])+circshift(tmp2,[0 -1])+circshift(tmp2,[-1 0])+circshift(tmp2,[1 0]);
36			tmp3=tmp3+circshift(tmp3,[0 1])+circshift(tmp3,[0 -1])+circshift(tmp3,[-1 0])+circshift(tmp3,[1 0]);
37			tmp1=find(isnan(field_mod)&tmp3>0);
38			field_mod(tmp1)=tmp2(tmp1)./tmp3(tmp1); tmp1=find(isnan(field_mod));
39			end
40
41			tmp1=[field_mod(end,:);field_mod;field_mod(1,:);field_mod(2,:)]';
42			field_int=interp2(x1,y1,tmp1,x2,y2)';
43
44			eval([deblank(var_ecco(fcur,:)) '=field_int.*mask;']);
45
46			end
47
48			rhoConstFresh=999.8;
49			emp=emp*rhoConstFresh; %m/s->kg/m^2/s
50			ustress=-ustress.*mask;
51			vstress=-vstress.*mask;
52
53			turbflux=turbflux.*mask;
54			emp=emp.*mask;
55			swflux=swflux.*mask;
56
57			mystruct_out=struct('ustress',ustress,'vstress',vstress,'emp',emp,'turbflux',turbflux,'swflux',swflux);
58