function []=profiles_prep_test_atlas(thresholdVAL);

global rep_out filename_out;
global t_std s_std t_test s_test t_w s_w t_equi s_equi;
global ymd hms pnum lon lat direc dmod ilon ilat imonth;
global fill_value_output initChkFile;

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;


%cost for T :

t_equi=squeeze(T_climatology(imonth,ilon,ilat,:))'; t_equi = interp1q(dmod_ref',t_equi',dmod')';
t_equi2=squeeze(T_climatology2(imonth,ilon,ilat,:))'; t_equi2 = interp1q(dmod_ref',t_equi2',dmod')';
t_equi3=squeeze(T_climatology3(imonth,ilon,ilat,:))'; t_equi3 = interp1q(dmod_ref',t_equi3',dmod')';
t_equi4=squeeze(T_climatology4(imonth,ilon,ilat,:))'; t_equi4 = interp1q(dmod_ref',t_equi4',dmod')';

t_cost=t_w.*((t_std-t_equi).^2); 
t_cost2=t_w.*((t_std-t_equi2).^2);
t_cost3=t_w.*((t_std-t_equi3).^2);
t_cost4=t_w.*((t_std-t_equi4).^2);

%cost for S : 

s_equi=squeeze(S_climatology(imonth,ilon,ilat,:))'; s_equi = interp1q(dmod_ref',s_equi',dmod')';
s_equi2=squeeze(S_climatology2(imonth,ilon,ilat,:))'; s_equi2 = interp1q(dmod_ref',s_equi2',dmod')';
s_equi3=squeeze(S_climatology3(imonth,ilon,ilat,:))'; s_equi3 = interp1q(dmod_ref',s_equi3',dmod')';
s_equi4=squeeze(S_climatology4(imonth,ilon,ilat,:))'; s_equi4 = interp1q(dmod_ref',s_equi4',dmod')';

s_cost=s_w.*((s_std-s_equi).^2);
s_cost2=s_w.*((s_std-s_equi2).^2);
s_cost3=s_w.*((s_std-s_equi3).^2);
s_cost4=s_w.*((s_std-s_equi4).^2);


%remarks: 1) the climatology fields are assumed to be masked with NaN
%	and the next two lines does not do anything for cost=NaN
%	logic: we cannot test the data, this is no reason to reject observations
%	2) the test must fail for all atlases for rejection 
%	logic: an individual atlas may be wrong in some places, this is no reason to reject observations

t_test(find( (t_test==0)&(t_std~=fill_value_output)&(t_cost>thresholdVAL)&(t_cost2>thresholdVAL)&(t_cost3>thresholdVAL)&(t_cost4>thresholdVAL) ))=5;
s_test(find( (s_test==0)&(s_std~=fill_value_output)&(s_cost>thresholdVAL)&(s_cost2>thresholdVAL)&(s_cost3>thresholdVAL)&(s_cost4>thresholdVAL) ))=5;