global dmod_ref vec_lon vec_lat tmask3D;
global T_climatology S_climatology T_climatology2 S_climatology2;
global T_climatology3 S_climatology3 T_climatology4 S_climatology4;
global T_weights S_weights;

if isempty(T_climatology)

%	a) ECCOv2 depth levels
dmod_ref=1.0e+03 * [  0.0050   0.0150   0.0275   0.0450   0.0650   0.0875   0.1175  ...
0.1600   0.2225   0.3100   0.4350   0.6100   0.8475   1.1600   1.5425   1.9750  ...
2.4500   2.9500   3.4500   3.9500   4.4500   4.9500   5.4500 ];

%       b) ECCOv2 grid point coordinates :
vec_lon=[0.5:359.5]; vec_lat=[-79.5:79.5];

%       c) load (ECCOv2-grid) monthly atlases :
rep_clim='/bay/data1/king/GOURETSKI_climatology/gridded/outliers/', rep_clim2='/bay/data1/king/NCEP_model/LEVITUS/',
rep_clim3='/raid0/gforget/DATAbin/HYDROatlases/woa2005/',rep_clim4='/net/ross/raid2/gforget/gaelmodel0406clim/',
T_climatology=zeros(12,360,160,23); S_climatology=T_climatology;
T_climatology2=zeros(12,360,160,23); S_climatology2=T_climatology2;
T_climatology3=zeros(12,360,160,23); S_climatology3=T_climatology3;
T_climatology4=zeros(12,360,160,23); S_climatology4=T_climatology4;
fid0=fopen([rep_clim 'LEV_wghc_ptmp_r4'],'r','b');
fid1=fopen([rep_clim 'LEV_wghc_salt_r4'],'r','b');
fid02=fopen([rep_clim2 'LEVITUS_1x1_ptmp_NEW_corK'],'r','b');
fid12=fopen([rep_clim2 'LEVITUS_1x1_salt_NEW_corK'],'r','b');
fid03=fopen([rep_clim3 'wod05_ptmp23'],'r','b');
fid13=fopen([rep_clim3 'wod05_salt23'],'r','b');
fid04=fopen([rep_clim4 'gm0406_ptmp23'],'r','b');
fid14=fopen([rep_clim4 'gm0406_salt23'],'r','b');
for tcur=1:12
field_tmp=fread(fid0,[360 160*23],'float32');
T_climatology(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid1,[360 160*23],'float32');
S_climatology(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid02,[360 160*23],'float32');
T_climatology2(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid12,[360 160*23],'float32');
S_climatology2(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid03,[360 160*23],'float32');
T_climatology3(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid13,[360 160*23],'float32');
S_climatology3(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid04,[360 160*23],'float32');
T_climatology4(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid14,[360 160*23],'float32');
S_climatology4(tcur,:,:,:)=reshape(field_tmp,[360 160 23]);
end%for tcur=1:12
fclose(fid0);fclose(fid1); fclose(fid02);fclose(fid12);
fclose(fid03);fclose(fid13); fclose(fid04);fclose(fid14);

test1=find((T_climatology.*S_climatology)==0);
T_climatology(test1)=NaN; S_climatology(test1)=NaN;
test1=find((T_climatology2.*S_climatology2)==0);
T_climatology2(test1)=NaN; S_climatology2(test1)=NaN;
test1=find((T_climatology3.*S_climatology3)==0);
T_climatology3(test1)=NaN; S_climatology3(test1)=NaN;
test1=find((T_climatology4.*S_climatology4)==0);
T_climatology4(test1)=NaN; S_climatology4(test1)=NaN;

%       d) load (ECCOv2-grid) monthly atlases :
%
%       by assumption: the uncertainty field contains non-zero values
%       (which avoid the complication of handling horiz interpolation here)
%       and we do not mask the data (the model will do this, given a mask)
rep_weight='/net/ross/raid0/gforget/MYweights/';
T_weights=zeros(360,160,23); S_weights=T_weights;
fid0=fopen([rep_weight 'wvar_Tr3.bin'],'r','b');
fid1=fopen([rep_weight 'wvar_Sr3.bin'],'r','b');
field_tmp=fread(fid0,[360 160*23],'float32');
T_weights=reshape(field_tmp,[360 160 23]);
field_tmp=fread(fid1,[360 160*23],'float32');
S_weights=reshape(field_tmp,[360 160 23]);

%       e) compute mask for tests:
tmask3D=squeeze(mean(S_climatology,1)).*squeeze(mean(S_climatology2,1)).*...
squeeze(mean(S_climatology3,1)).*squeeze(mean(S_climatology4,1)).*...
squeeze(mean(T_climatology,1)).*squeeze(mean(T_climatology2,1)).*...
squeeze(mean(T_climatology3,1)).*squeeze(mean(T_climatology4,1));
tmask3D(find(~isnan(tmask3D)))=1;

for tcur=1:12
T_climatology(tcur,:,:,:)=squeeze(T_climatology(tcur,:,:,:)).*tmask3D; S_climatology(tcur,:,:,:)=squeeze(S_climatology(tcur,:,:,:)).*tmask3D;
T_climatology2(tcur,:,:,:)=squeeze(T_climatology2(tcur,:,:,:)).*tmask3D; S_climatology2(tcur,:,:,:)=squeeze(S_climatology2(tcur,:,:,:)).*tmask3D;
T_climatology3(tcur,:,:,:)=squeeze(T_climatology3(tcur,:,:,:)).*tmask3D; S_climatology3(tcur,:,:,:)=squeeze(S_climatology3(tcur,:,:,:)).*tmask3D;
T_climatology4(tcur,:,:,:)=squeeze(T_climatology4(tcur,:,:,:)).*tmask3D; S_climatology4(tcur,:,:,:)=squeeze(S_climatology4(tcur,:,:,:)).*tmask3D;
end

end