| 1 |
% t_test, s_test description |
| 2 |
% 0 = valid data |
| 3 |
% 1 = outside +/-79 lat |
| 4 |
% 2 = absurd sal value |
| 5 |
% 3 = doubtful profiler (based on our own evaluations) |
| 6 |
% 4 = doubtful profiler (based on Argo grey list) |
| 7 |
% 5 = high climatology/atlas cost - all four of them |
| 8 |
% 6 = bad Pressure vector |
| 9 |
|
| 10 |
%clear |
| 11 |
|
| 12 |
global rep_out filename_out; |
| 13 |
global t_std s_std t_test s_test t_w s_w t_equi s_equi; |
| 14 |
global ymd hms pnum pnum_txt lon lat direc dmod ilon ilat imonth; |
| 15 |
global fill_value_output initChkFile; |
| 16 |
|
| 17 |
profiles_prep_load_fields; |
| 18 |
|
| 19 |
%depths of standard levels : |
| 20 |
dmod=[[5:10:185] [200:20:500] [550:50:1000] [1100:100:6000]]; |
| 21 |
%maximum level to compute : e.g. 2000m for ARGO profiles |
| 22 |
lev_max=max(find(dmod<=2000)); fprintf(['level max : ' num2str(lev_max) ' \n']); |
| 23 |
dmod=dmod(1:lev_max); |
| 24 |
%fill value for the output files : |
| 25 |
fill_value_output=-9999; |
| 26 |
%do we start from in situ temperature? |
| 27 |
TfromINSITUtoPOT=1; %1 => we will convert to potential temperature |
| 28 |
%do we start from P vertical coordinate? |
| 29 |
PfromPtoZ=1; %1 => we will convert to depth |
| 30 |
|
| 31 |
%loop over the data files : |
| 32 |
for rep_nb=1:3 |
| 33 |
clear fid*; %needed to properly initialize output files |
| 34 |
switch rep_nb |
| 35 |
case 1 |
| 36 |
bassin_cur='INDIAN'; |
| 37 |
case 2 |
| 38 |
bassin_cur='PACIFIC'; |
| 39 |
case 3 |
| 40 |
bassin_cur='ATLANTIC'; |
| 41 |
end |
| 42 |
rep_in=['/net/altix3700/raid4/gforget/ARGO/ifremer/' bassin_cur '/']; |
| 43 |
rep_out='/net/altix3700/raid4/gforget/ARGO/ifremer/ECCOformat/'; filename_out=[bassin_cur '_ARGO_EDW']; |
| 44 |
|
| 45 |
initChkFile=1; |
| 46 |
|
| 47 |
%get the number of files to be treated : |
| 48 |
fid=fopen(['/net/altix3700/raid4/gforget/ARGO/ifremer/NC_list_' bassin_cur],'r'); |
| 49 |
foo=fread(fid); |
| 50 |
nfiles=size(foo); |
| 51 |
nfiles=nfiles(1)/17; |
| 52 |
fclose(fid); |
| 53 |
|
| 54 |
%file with the list of source files : |
| 55 |
fid=fopen(['/net/altix3700/raid4/gforget/ARGO/ifremer/NC_list_' bassin_cur],'r'); |
| 56 |
|
| 57 |
%loop over data files : |
| 58 |
ktQC=0*ones(1,9);ksQC=0*ones(1,9); |
| 59 |
|
| 60 |
for nf=1:nfiles % FILE LOOP |
| 61 |
|
| 62 |
|
| 63 |
if mod(nf,100)==0|nf==1; fprintf([num2str(nf) ' ' num2str(nfiles) '\n']); end; |
| 64 |
|
| 65 |
clear PRES* PSAL* TEMP*; |
| 66 |
TEMP=8888*ones(500,1);PSAL=8888*ones(500,1); |
| 67 |
name=fread(fid,16,'char');name= setstr(name);name=name'; |
| 68 |
junk=fread(fid,1,'char'); |
| 69 |
|
| 70 |
if ( rep_nb~=1|isempty(findstr(name,'20050911_prof')) ) |
| 71 |
|
| 72 |
eval(['ncload ' rep_in name ';']); |
| 73 |
|
| 74 |
%some dimensions inversions exist for _QC variables (!?) |
| 75 |
if size(PRES_QC,1)~=size(PRES,1) |
| 76 |
PRES_QC=PRES_QC'; |
| 77 |
if ~isempty(who('PRES_ADJUSTED_QC'));PRES_ADJUSTED_QC=PRES_ADJUSTED_QC';end; |
| 78 |
if ~isempty(who('PSAL_QC'));PSAL_QC=PSAL_QC';end; |
| 79 |
if ~isempty(who('PSAL_ADJUSTED_QC'));PSAL_ADJUSTED_QC=PSAL_ADJUSTED_QC';end; |
| 80 |
if ~isempty(who('TEMP_QC'));TEMP_QC=TEMP_QC';end; |
| 81 |
if ~isempty(who('TEMP_ADJUSTED_QC'));TEMP_ADJUSTED_QC=TEMP_ADJUSTED_QC';end; |
| 82 |
end |
| 83 |
|
| 84 |
if size(PLATFORM_NUMBER,1)~=size(PRES,1); PLATFORM_NUMBER=PLATFORM_NUMBER'; end; |
| 85 |
|
| 86 |
%get the fillvalues: |
| 87 |
nc = netcdf([rep_in name], 'nowrite'); |
| 88 |
PRES_fillval=fillval(nc{'PRES'}); |
| 89 |
if PSAL(1,1) ~= 8888; PSAL_fillval=fillval(nc{'PSAL'});end; |
| 90 |
if TEMP(1,1) ~= 8888; TEMP_fillval=fillval(nc{'TEMP'});end; |
| 91 |
PLATFORM_NUMBER_fillval=fillval(nc{'PLATFORM_NUMBER'}); |
| 92 |
nc = close(nc); |
| 93 |
|
| 94 |
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); |
| 95 |
VECdateOBS=REFdateOBS+JULD; |
| 96 |
|
| 97 |
xx=size(PRES);mloc=xx(1);mpts=xx(2); |
| 98 |
|
| 99 |
for m=1:mloc % LOCATION LOOP |
| 100 |
|
| 101 |
%date, position, etc |
| 102 |
|
| 103 |
tmp1=greg_0h(VECdateOBS(m)); |
| 104 |
ymd=tmp1(1)*1e4+tmp1(2)*1e2+tmp1(3); |
| 105 |
hms=tmp1(4)*1e4+tmp1(5)*1e2+tmp1(6); |
| 106 |
|
| 107 |
imonth=tmp1(2); |
| 108 |
|
| 109 |
lat=LATITUDE(m); |
| 110 |
lon=LONGITUDE(m); if lon < 0; lon=lon+360;end; |
| 111 |
|
| 112 |
ilon=find(((vec_lon-lon).^2)==min((vec_lon-lon).^2)); ilon=ilon(1); |
| 113 |
ilat=find(((vec_lat-lat).^2)==min((vec_lat-lat).^2)); ilat=ilat(1); |
| 114 |
|
| 115 |
direc=0;if(DIRECTION(m)=='A');direc=1;end;if(DIRECTION(m)=='D');direc=2;end |
| 116 |
|
| 117 |
pnum_txt=deblank(PLATFORM_NUMBER(m,:)); pnum_txt=pnum_txt(pnum_txt~=PLATFORM_NUMBER_fillval); |
| 118 |
pnum=double(pnum_txt); pnum=pnum(find(pnum>=48&pnum<=57)); pnum=str2num(char(pnum)); |
| 119 |
if isempty(pnum); pnum_txt='9999'; pnum=9999; end; |
| 120 |
|
| 121 |
%observations |
| 122 |
|
| 123 |
p=PRES_ADJUSTED(m,:); |
| 124 |
p_QC=PRES_ADJUSTED_QC(m,:); |
| 125 |
if isempty(find(p~=PRES_fillval)); p=PRES(m,:); p_QC=PRES_QC(m,:); end |
| 126 |
|
| 127 |
tmp1=find(isnan(p)); p(tmp1)=PRES_fillval; p_QC(tmp1)='5'; |
| 128 |
for n=1:length(p)-1; |
| 129 |
tmp1=find(p(n+1:end)==p(n)&p(n+1:end)~=PRES_fillval); |
| 130 |
p(n+tmp1)=PRES_fillval;p_QC(n+tmp1)='5'; |
| 131 |
end |
| 132 |
|
| 133 |
bad_P=0; |
| 134 |
tmp1=find(p_QC=='4'); |
| 135 |
if(length(tmp1)<=5); |
| 136 |
%get rid of these few bad points and keep the profile |
| 137 |
p(tmp1)=PRES_fillval;p_QC(tmp1)='5'; |
| 138 |
else; |
| 139 |
%flag the profile (will be masked in the main file) |
| 140 |
%but keep the bad points (to interp and be able to CHECK) |
| 141 |
bad_P=1; |
| 142 |
end; |
| 143 |
|
| 144 |
if TEMP(m,1) ~= 8888 |
| 145 |
t=TEMP_ADJUSTED(m,:); |
| 146 |
t_QC=TEMP_ADJUSTED_QC(m,:); |
| 147 |
t_ERR=TEMP_ADJUSTED_ERROR(m,:);tf=find(t_ERR==TEMP_fillval|isnan(t_ERR));t_ERR(tf)=0; |
| 148 |
if isempty(find(t~=TEMP_fillval)); t=TEMP(m,:); t_QC=TEMP_QC(m,:); end |
| 149 |
end |
| 150 |
|
| 151 |
if PSAL(m,1) ~= 8888 |
| 152 |
s=PSAL_ADJUSTED(m,:); |
| 153 |
s_QC=PSAL_ADJUSTED_QC(m,:); |
| 154 |
s_ERR=PSAL_ADJUSTED_ERROR(m,:);sf=find(s_ERR==PSAL_fillval|isnan(s_ERR));,s_ERR(sf)=0; |
| 155 |
if isempty(find(s~=PSAL_fillval)); s=PSAL(m,:); s_QC=PSAL_QC(m,:); end |
| 156 |
end |
| 157 |
|
| 158 |
%convert pressure to depth (if necessary) |
| 159 |
if PfromPtoZ |
| 160 |
tmp1=find((p~=PRES_fillval)&(~isnan(p))); |
| 161 |
if ~isempty(tmp1); p( tmp1 ) = find_depth(p(tmp1),lat); end; |
| 162 |
end |
| 163 |
|
| 164 |
%interpolation for T : |
| 165 |
z_std=dmod; t_std=NaN*z_std; tE_std=t_std; |
| 166 |
if TEMP(m,1) ~= 8888 |
| 167 |
tmp1=find((t_QC=='1'|t_QC=='2') & (p~=PRES_fillval&~isnan(p))); |
| 168 |
z_in=p(tmp1); t_in=t(tmp1); tE_in=t_ERR(tmp1); |
| 169 |
if length(t_in)>5 %...expected to avoid isolated values |
| 170 |
t_std = interp1(z_in,t_in,z_std); |
| 171 |
t_std=profiles_prep_test_interp(t_std,z_std,z_in,NaN); |
| 172 |
tE_std = interp1(z_in,tE_in,z_std); |
| 173 |
tE_std=profiles_prep_test_interp(tE_std,z_std,z_in,NaN); |
| 174 |
end |
| 175 |
end%z_std=gdept; t_std=NaN*z_std; |
| 176 |
t_std(find(isnan(t_std)))=fill_value_output; |
| 177 |
tE_std(find(isnan(tE_std)))=fill_value_output; |
| 178 |
|
| 179 |
%interpolation for S : |
| 180 |
z_std=dmod; s_std=NaN*z_std; sE_std=s_std; |
| 181 |
if PSAL(m,1) ~= 8888 |
| 182 |
tmp1=find((s_QC=='1'|s_QC=='2') & (p~=PRES_fillval&~isnan(p))); |
| 183 |
z_in=p(tmp1); s_in=s(tmp1); sE_in=s_ERR(tmp1); |
| 184 |
if length(s_in)>5 |
| 185 |
s_std = interp1(z_in,s_in,z_std); |
| 186 |
s_std=profiles_prep_test_interp(s_std,z_std,z_in,NaN); |
| 187 |
sE_std = interp1(z_in,sE_in,z_std); |
| 188 |
sE_std=profiles_prep_test_interp(sE_std,z_std,z_in,NaN); |
| 189 |
end |
| 190 |
end%if PSAL(m,1) ~= 8888 |
| 191 |
s_std(find(isnan(s_std)))=fill_value_output; |
| 192 |
sE_std(find(isnan(sE_std)))=fill_value_output; |
| 193 |
|
| 194 |
%convert T in situ to T potential (if necessary) |
| 195 |
if TfromINSITUtoPOT |
| 196 |
tmpP=0.981*1.027*dmod; tmpS=35*ones(size(t_std)); |
| 197 |
tmpIND=find(t_std~=fill_value_output); |
| 198 |
if ~isempty(tmpIND) |
| 199 |
t_std(tmpIND)=sw_ptmp(tmpS(tmpIND),t_std(tmpIND),tmpP(tmpIND),0); |
| 200 |
end |
| 201 |
end |
| 202 |
|
| 203 |
|
| 204 |
%attribute weights: (modified later in profiles_prep_write_nc) |
| 205 |
% |
| 206 |
% by assumption: the uncertainty field contains non-zero values |
| 207 |
% (which avoid the complication of handling horiz interpolation here) |
| 208 |
% and we do not mask the data (the model will do this, given a mask) |
| 209 |
% |
| 210 |
t_w = squeeze(T_weights(ilon,ilat,:)); t_w = interp1(dmod_ref',t_w',dmod')'.^2; |
| 211 |
tmp1=find(tE_std~=fill_value_output); t_w(tmp1)=t_w(tmp1)+tE_std(tmp1).^2; |
| 212 |
s_w = squeeze(S_weights(ilon,ilat,:)); s_w = interp1(dmod_ref',s_w',dmod')'.^2; |
| 213 |
tmp1=find(sE_std~=fill_value_output); s_w(tmp1)=s_w(tmp1)+sE_std(tmp1).^2; |
| 214 |
t_w=t_w.^-1; s_w=s_w.^-1; |
| 215 |
|
| 216 |
if ~isempty(find(isnan(t_w.*s_w))); fprintf('error in weighting'); return; end; |
| 217 |
|
| 218 |
|
| 219 |
%tests of the data values : |
| 220 |
%-------------------------- |
| 221 |
% t_test>0 or s_test>0 implies the data value is rejected, |
| 222 |
% the value of t_test indicates on what criterium |
| 223 |
t_test=zeros(size(t_std)); s_test=t_test; |
| 224 |
|
| 225 |
%test of the extreme latitudes : |
| 226 |
if abs(lat)>79; t_test(:)=1; s_test(:)=1; end; |
| 227 |
|
| 228 |
%test of "absurd" salinity values : test |
| 229 |
s_test(find( (s_std>42)&(s_std~=fill_value_output) ))=2; |
| 230 |
s_test(find( (s_std<25)&(s_std~=fill_value_output) ))=2; |
| 231 |
|
| 232 |
%tests of doubtful profilers : |
| 233 |
%profiles_prep_test_profiler2; |
| 234 |
%profiles_prep_test_profiler3; |
| 235 |
%profiles_prep_test_profiler4; |
| 236 |
|
| 237 |
%tests of the cost VS the climatology : |
| 238 |
profiles_prep_test_atlas(50); |
| 239 |
|
| 240 |
%bad pressure flag: |
| 241 |
if bad_P==1; t_test(:)=6; s_test(:)=6; end; |
| 242 |
|
| 243 |
%store/write results: |
| 244 |
profiles_prep_write_nc(1); |
| 245 |
|
| 246 |
end % LOCATION LOOP |
| 247 |
end %if ( rep_nb~=1|isempty(findstr(name,'20050911_prof') ) |
| 248 |
end % FILE LOOP |
| 249 |
profiles_prep_write_nc(2); |
| 250 |
end%for rep_nb=1:3 |
| 251 |
|
| 252 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph