% t_test, s_test description
%   0 = valid data
%   1 = outside +/-79 lat
%   2 = absurd sal value
%   3 = doubtful profiler (based on our own evaluations)
%   4 = doubtful profiler (based on Argo grey list)
%   5 = high climatology/atlas cost - all four of them
%   6 = bad Pressure vector

%clear

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 pnum_txt lon lat direc dmod ilon ilat imonth;
global fill_value_output initChkFile; 

profiles_prep_load_fields;

%depths of standard levels :
dmod=[[5:10:185] [200:20:500] [550:50:1000] [1100:100:6000]];
%maximum level to compute : e.g. 2000m for ARGO profiles
lev_max=max(find(dmod<=2000)); fprintf(['level max : ' num2str(lev_max) ' \n']);
dmod=dmod(1:lev_max);
%fill value for the output files :
fill_value_output=-9999;
%do we start from in situ temperature?
TfromINSITUtoPOT=1; %1 => we will convert to potential temperature
%do we start from P vertical coordinate?
PfromPtoZ=1; %1 => we will convert to depth

%loop over the data files :
for rep_nb=1:3
clear fid*; %needed to properly initialize output files
switch rep_nb
case 1
bassin_cur='INDIAN';
case 2
bassin_cur='PACIFIC';
case 3
bassin_cur='ATLANTIC';
end
rep_in=['/net/altix3700/raid4/gforget/ARGO/ifremer/' bassin_cur '/'];
rep_out='/net/altix3700/raid4/gforget/ARGO/ifremer/ECCOformat/'; filename_out=[bassin_cur '_ARGO_EDW'];

initChkFile=1;

%get the number of files to be treated :
fid=fopen(['/net/altix3700/raid4/gforget/ARGO/ifremer/NC_list_' bassin_cur],'r');
foo=fread(fid);
nfiles=size(foo);
nfiles=nfiles(1)/17;
fclose(fid);

%file with the list of source files :
fid=fopen(['/net/altix3700/raid4/gforget/ARGO/ifremer/NC_list_' bassin_cur],'r');

%loop over data files : 
ktQC=0*ones(1,9);ksQC=0*ones(1,9);

for nf=1:nfiles        % FILE LOOP


if mod(nf,100)==0|nf==1; fprintf([num2str(nf) ' ' num2str(nfiles) '\n']); end;

clear PRES* PSAL* TEMP*;
TEMP=8888*ones(500,1);PSAL=8888*ones(500,1);
name=fread(fid,16,'char');name= setstr(name);name=name';
junk=fread(fid,1,'char');

if ( rep_nb~=1|isempty(findstr(name,'20050911_prof')) )

eval(['ncload ' rep_in name ';']);

%some dimensions inversions exist for _QC variables (!?)
if size(PRES_QC,1)~=size(PRES,1)
PRES_QC=PRES_QC'; 
if ~isempty(who('PRES_ADJUSTED_QC'));PRES_ADJUSTED_QC=PRES_ADJUSTED_QC';end;
if ~isempty(who('PSAL_QC'));PSAL_QC=PSAL_QC';end;
if ~isempty(who('PSAL_ADJUSTED_QC'));PSAL_ADJUSTED_QC=PSAL_ADJUSTED_QC';end;
if ~isempty(who('TEMP_QC'));TEMP_QC=TEMP_QC';end;
if ~isempty(who('TEMP_ADJUSTED_QC'));TEMP_ADJUSTED_QC=TEMP_ADJUSTED_QC';end;
end

if size(PLATFORM_NUMBER,1)~=size(PRES,1); PLATFORM_NUMBER=PLATFORM_NUMBER'; end; 

%get the fillvalues:
nc = netcdf([rep_in name], 'nowrite'); 
PRES_fillval=fillval(nc{'PRES'});
  if PSAL(1,1) ~= 8888; PSAL_fillval=fillval(nc{'PSAL'});end;
  if TEMP(1,1) ~= 8888; TEMP_fillval=fillval(nc{'TEMP'});end;
PLATFORM_NUMBER_fillval=fillval(nc{'PLATFORM_NUMBER'});
nc = close(nc);

REFdateOBS=str2num(strvcat(REFERENCE_DATE_TIME(1:4),REFERENCE_DATE_TIME(5:6),REFERENCE_DATE_TIME(7:8),REFERENCE_DATE_TIME(9:10),REFERENCE_DATE_TIME(11:12),REFERENCE_DATE_TIME(13:14)))'; REFdateOBS=jul_0h(REFdateOBS);
VECdateOBS=REFdateOBS+JULD;

xx=size(PRES);mloc=xx(1);mpts=xx(2);

for m=1:mloc    % LOCATION LOOP

%date, position, etc

tmp1=greg_0h(VECdateOBS(m));
ymd=tmp1(1)*1e4+tmp1(2)*1e2+tmp1(3);
hms=tmp1(4)*1e4+tmp1(5)*1e2+tmp1(6);

imonth=tmp1(2);

lat=LATITUDE(m);
lon=LONGITUDE(m); if lon < 0; lon=lon+360;end;

ilon=find(((vec_lon-lon).^2)==min((vec_lon-lon).^2)); ilon=ilon(1);
ilat=find(((vec_lat-lat).^2)==min((vec_lat-lat).^2)); ilat=ilat(1);

direc=0;if(DIRECTION(m)=='A');direc=1;end;if(DIRECTION(m)=='D');direc=2;end

pnum_txt=deblank(PLATFORM_NUMBER(m,:)); pnum_txt=pnum_txt(pnum_txt~=PLATFORM_NUMBER_fillval);
pnum=double(pnum_txt); pnum=pnum(find(pnum>=48&pnum<=57)); pnum=str2num(char(pnum));
if isempty(pnum); pnum_txt='9999'; pnum=9999; end;

%observations

 p=PRES_ADJUSTED(m,:);
 p_QC=PRES_ADJUSTED_QC(m,:);
 if isempty(find(p~=PRES_fillval)); p=PRES(m,:); p_QC=PRES_QC(m,:); end

        tmp1=find(isnan(p)); p(tmp1)=PRES_fillval; p_QC(tmp1)='5';
        for n=1:length(p)-1;
        tmp1=find(p(n+1:end)==p(n)&p(n+1:end)~=PRES_fillval);
        p(n+tmp1)=PRES_fillval;p_QC(n+tmp1)='5';
        end

	bad_P=0;
	tmp1=find(p_QC=='4');
	if(length(tmp1)<=5);
	%get rid of these few bad points and keep the profile
	p(tmp1)=PRES_fillval;p_QC(tmp1)='5';
	else;
	%flag the profile (will be masked in the main file) 
	%but keep the bad points (to interp and be able to CHECK)
	bad_P=1;
	end;

  if TEMP(m,1) ~= 8888
 t=TEMP_ADJUSTED(m,:);
 t_QC=TEMP_ADJUSTED_QC(m,:);
 t_ERR=TEMP_ADJUSTED_ERROR(m,:);tf=find(t_ERR==TEMP_fillval|isnan(t_ERR));t_ERR(tf)=0;
 if isempty(find(t~=TEMP_fillval)); t=TEMP(m,:); t_QC=TEMP_QC(m,:); end
  end

   if PSAL(m,1) ~= 8888
 s=PSAL_ADJUSTED(m,:);
 s_QC=PSAL_ADJUSTED_QC(m,:);
 s_ERR=PSAL_ADJUSTED_ERROR(m,:);sf=find(s_ERR==PSAL_fillval|isnan(s_ERR));,s_ERR(sf)=0;
 if isempty(find(s~=PSAL_fillval)); s=PSAL(m,:); s_QC=PSAL_QC(m,:); end
  end

%convert pressure to depth (if necessary)
if PfromPtoZ
tmp1=find((p~=PRES_fillval)&(~isnan(p)));
if ~isempty(tmp1); p( tmp1 ) = find_depth(p(tmp1),lat); end;
end

%interpolation for T : 
z_std=dmod; t_std=NaN*z_std; tE_std=t_std;
if TEMP(m,1) ~= 8888
tmp1=find((t_QC=='1'|t_QC=='2') & (p~=PRES_fillval&~isnan(p)));
z_in=p(tmp1); t_in=t(tmp1); tE_in=t_ERR(tmp1);
if length(t_in)>5 %...expected to avoid isolated values
t_std = interp1(z_in,t_in,z_std);
t_std=profiles_prep_test_interp(t_std,z_std,z_in,NaN);
tE_std = interp1(z_in,tE_in,z_std);
tE_std=profiles_prep_test_interp(tE_std,z_std,z_in,NaN);
end
end%z_std=gdept; t_std=NaN*z_std;
t_std(find(isnan(t_std)))=fill_value_output; 
tE_std(find(isnan(tE_std)))=fill_value_output; 

%interpolation for S :
z_std=dmod; s_std=NaN*z_std; sE_std=s_std;
if PSAL(m,1) ~= 8888
tmp1=find((s_QC=='1'|s_QC=='2') & (p~=PRES_fillval&~isnan(p)));
z_in=p(tmp1); s_in=s(tmp1); sE_in=s_ERR(tmp1);
if length(s_in)>5
s_std = interp1(z_in,s_in,z_std);
s_std=profiles_prep_test_interp(s_std,z_std,z_in,NaN);
sE_std = interp1(z_in,sE_in,z_std);
sE_std=profiles_prep_test_interp(sE_std,z_std,z_in,NaN);
end
end%if PSAL(m,1) ~= 8888
s_std(find(isnan(s_std)))=fill_value_output; 
sE_std(find(isnan(sE_std)))=fill_value_output; 

%convert T in situ to T potential (if necessary)
if TfromINSITUtoPOT
tmpP=0.981*1.027*dmod; tmpS=35*ones(size(t_std)); 
tmpIND=find(t_std~=fill_value_output);
if ~isempty(tmpIND)
t_std(tmpIND)=sw_ptmp(tmpS(tmpIND),t_std(tmpIND),tmpP(tmpIND),0);
end
end


%attribute weights: (modified later in profiles_prep_write_nc)
%
%	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)
%
t_w = squeeze(T_weights(ilon,ilat,:)); t_w = interp1(dmod_ref',t_w',dmod')'.^2;
tmp1=find(tE_std~=fill_value_output); t_w(tmp1)=t_w(tmp1)+tE_std(tmp1).^2;
s_w = squeeze(S_weights(ilon,ilat,:)); s_w = interp1(dmod_ref',s_w',dmod')'.^2;
tmp1=find(sE_std~=fill_value_output); s_w(tmp1)=s_w(tmp1)+sE_std(tmp1).^2;
t_w=t_w.^-1; s_w=s_w.^-1;

if ~isempty(find(isnan(t_w.*s_w))); fprintf('error in weighting'); return; end;


%tests of the data values :
%--------------------------
%	t_test>0 or s_test>0 implies the data value is rejected,
% 	the value of t_test indicates on what criterium 
t_test=zeros(size(t_std)); s_test=t_test;

%test of the extreme latitudes : 
if abs(lat)>79; t_test(:)=1; s_test(:)=1; end;

%test of "absurd" salinity values : test
s_test(find( (s_std>42)&(s_std~=fill_value_output) ))=2;
s_test(find( (s_std<25)&(s_std~=fill_value_output) ))=2;

%tests of doubtful profilers :
%profiles_prep_test_profiler2;
%profiles_prep_test_profiler3;
%profiles_prep_test_profiler4;

%tests of the cost VS the climatology :
profiles_prep_test_atlas(50);

%bad pressure flag:
if bad_P==1; t_test(:)=6; s_test(:)=6; end;

%store/write results:
profiles_prep_write_nc(1);

end      % LOCATION LOOP
end %if ( rep_nb~=1|isempty(findstr(name,'20050911_prof') )
end      % FILE LOOP
profiles_prep_write_nc(2);
end%for rep_nb=1:3