Skip to content
54 changes: 46 additions & 8 deletions matlab_codes/calculate_piecewisefit.m
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur

%
% Cecile Cabanes, June 2013 : calculate off-diagonal terms for error estimate: add horizontal covariance to track changes: see "change config 129"
% Delphine Dobler (DD), September 2024 :
% 4 - save computed theta levels in a file

%pn_float_dir='testfloats/';
%pn_float_name='robbins4900178';
Expand Down Expand Up @@ -156,10 +158,15 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
return
end

% DD (2024/09-4): initialisation of arrays to be saved in a file
Theta_nseq=NaN*zeros(10,n_seq);
Index=NaN*zeros(10,n);
Plevels=NaN*zeros(10,n_seq);

% loop through sequences of calseries, ie if time series is broken into segments --
for i=1:n_seq
for iseq=1:n_seq

calindex = find(calseries==unique_cal(i));
calindex = find(calseries==unique_cal(iseq));
k = length(calindex);

% choose 10 float theta levels to use in the piecewise linear fit --------
Expand All @@ -180,6 +187,31 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur

[Theta, P, index, var_s_th, th] =...
find_10thetas( unique_SAL, unique_PTMP, unique_PRES, unique_la_ptmp, use_theta_gt, use_theta_lt, use_pres_gt, use_pres_lt, use_percent_gt);

% DD (2024/09-4): Save Theta-related information for plotting functions
Index(:,calindex)=index;
Plevels(:,iseq)=P;
[mvth,~]=size(var_s_th);
Theta_nseq(:,iseq)=Theta;
if iseq==1
Var_s_Thetas=NaN.*ones(mvth,n_seq);
Thetas=NaN.*ones(mvth,n_seq);
Var_s_Thetas(:,iseq)=var_s_th;
Thetas(:,iseq)=th;
else
[mVth,~]=size(Var_s_Thetas);
tmp1=Var_s_Thetas;
tmp2=Thetas;
mm=max(mVth,mvth);
Var_s_Thetas=NaN.*ones(mm,n_seq);
Thetas=NaN.*ones(mm,n_seq);
Var_s_Thetas(1:mVth,:)=tmp1;
Var_s_Thetas(1:mvth,iseq)=var_s_th;
Thetas(1:mVth,:)=tmp2;
Thetas(1:mvth,iseq)=th;
end
%


pp=find(isnan(index)==0);
if(isempty(pp)==0) % only proceed when there are valid levels ----
Expand Down Expand Up @@ -229,18 +261,18 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
if isempty(breaks)
[xfit(calindex), pcond_factor(calindex), pcond_factor_err(calindex), time_deriv(calindex), ...
time_deriv_err(calindex), sta_mean(calindex), sta_rms(calindex), NDF, fitcoef, fitbreaks] = ...
fit_cond(x, y, err, covariance, 'max_no_breaks', max_breaks(i));
fit_cond(x, y, err, covariance, 'max_no_breaks', max_breaks(iseq));
else
breaks_in = breaks(i,:);
breaks_in = breaks(iseq,:);
breaks_in = breaks_in(find(isfinite(breaks_in)));
if isempty(max_breaks(i))
if isempty(max_breaks(iseq))
[xfit(calindex), pcond_factor(calindex), pcond_factor_err(calindex), time_deriv(calindex), ...
time_deriv_err(calindex), sta_mean(calindex), sta_rms(calindex), NDF, fitcoef, fitbreaks] = ...
fit_cond(x, y, err, covariance, 'breaks', breaks_in);
else
[xfit(calindex), pcond_factor(calindex), pcond_factor_err(calindex), time_deriv(calindex), ...
time_deriv_err(calindex), sta_mean(calindex), sta_rms(calindex), NDF, fitcoef, fitbreaks] = ...
fit_cond(x, y, err, covariance, 'breaks', breaks_in, 'max_no_breaks', max_breaks(i));
fit_cond(x, y, err, covariance, 'breaks', breaks_in, 'max_no_breaks', max_breaks(iseq));
end
end

Expand All @@ -266,16 +298,22 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
sta_SAL1(:,calindex) = sw_salt( (sta_COND(:,calindex)+sta_COND_err(:,calindex))/sw_c3515, unique_PTMP, 0 );
sta_SAL_err(:,calindex) = abs(sta_SAL(:,calindex)-sta_SAL1(:,calindex));

fcoef(i,1:length(fitcoef)) = fitcoef;
fcoef(iseq,1:length(fitcoef)) = fitcoef;
if ~isempty(fitbreaks)
fbreaks(i,1:length(fitbreaks)) = fitbreaks;
fbreaks(iseq,1:length(fitbreaks)) = fitbreaks;
end

end

end %if there are valid levels
end %for each unique_cal

% DD (2024/09-4) : save theta levels information inside a file.
Theta=Theta_nseq;
ls_theta = strcat( po_system_configuration.FLOAT_CALIB_DIRECTORY, ...
pn_float_dir, 'selected_theta_', pn_float_name,'.mat');
save(ls_theta, 'Theta','Index','Plevels','Var_s_Thetas','Thetas')
% end of DD (2024/09-4)

% save calibration data --------------------------------

Expand Down
111 changes: 111 additions & 0 deletions matlab_codes/create_la_wmo_boxes_file.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
function create_la_wmo_boxes_file( po_system_configuration )

% function create_la_wmo_boxes_file( po_system_configuration )
%
% Delphine Dobler, August 2024

clim_ctd_=strcat(po_system_configuration.HISTORICAL_DIRECTORY,po_system_configuration.HISTORICAL_CTD_PREFIX);
clim_argo_=strcat(po_system_configuration.HISTORICAL_DIRECTORY,po_system_configuration.HISTORICAL_ARGO_PREFIX);
clim_bottle_=strcat(po_system_configuration.HISTORICAL_DIRECTORY,po_system_configuration.HISTORICAL_BOTTLE_PREFIX);
la_wmo_boxes_file=strcat( po_system_configuration.CONFIG_DIRECTORY, po_system_configuration.CONFIG_WMO_BOXES);

la_wmo_boxes=zeros(648,4);
la_wmo_boxes(:,1)=[1800;1700;1600;1500;1400;1300;1200;1100;1000;3000;3100;...
3200;3300;3400;3500;3600;3700;3800;1801;1701;1601;1501;...
1401;1301;1201;1101;1001;3001;3101;3201;3301;3401;3501;...
3601;3701;3801;1802;1702;1602;1502;1402;1302;1202;1102;...
1002;3002;3102;3202;3302;3402;3502;3602;3702;3802;1803;...
1703;1603;1503;1403;1303;1203;1103;1003;3003;3103;3203;...
3303;3403;3503;3603;3703;3803;1804;1704;1604;1504;1404;...
1304;1204;1104;1004;3004;3104;3204;3304;3404;3504;3604;...
3704;3804;1805;1705;1605;1505;1405;1305;1205;1105;1005;...
3005;3105;3205;3305;3405;3505;3605;3705;3805;1806;1706;...
1606;1506;1406;1306;1206;1106;1006;3006;3106;3206;3306;...
3406;3506;3606;3706;3806;1807;1707;1607;1507;1407;1307;...
1207;1107;1007;3007;3107;3207;3307;3407;3507;3607;3707;...
3807;1808;1708;1608;1508;1408;1308;1208;1108;1008;3008;...
3108;3208;3308;3408;3508;3608;3708;3808;1809;1709;1609;...
1509;1409;1309;1209;1109;1009;3009;3109;3209;3309;3409;...
3509;3609;3709;3809;1810;1710;1610;1510;1410;1310;1210;...
1110;1010;3010;3110;3210;3310;3410;3510;3610;3710;3810;...
1811;1711;1611;1511;1411;1311;1211;1111;1011;3011;3111;...
3211;3311;3411;3511;3611;3711;3811;1812;1712;1612;1512;...
1412;1312;1212;1112;1012;3012;3112;3212;3312;3412;3512;...
3612;3712;3812;1813;1713;1613;1513;1413;1313;1213;1113;...
1013;3013;3113;3213;3313;3413;3513;3613;3713;3813;1814;...
1714;1614;1514;1414;1314;1214;1114;1014;3014;3114;3214;...
3314;3414;3514;3614;3714;3814;1815;1715;1615;1515;1415;...
1315;1215;1115;1015;3015;3115;3215;3315;3415;3515;3615;...
3715;3815;1816;1716;1616;1516;1416;1316;1216;1116;1016;...
3016;3116;3216;3316;3416;3516;3616;3716;3816;1817;1717;...
1617;1517;1417;1317;1217;1117;1017;3017;3117;3217;3317;...
3417;3517;3617;3717;3817;7817;7717;7617;7517;7417;7317;...
7217;7117;7017;5017;5117;5217;5317;5417;5517;5617;5717;...
5817;7816;7716;7616;7516;7416;7316;7216;7116;7016;5016;...
5116;5216;5316;5416;5516;5616;5716;5816;7815;7715;7615;...
7515;7415;7315;7215;7115;7015;5015;5115;5215;5315;5415;...
5515;5615;5715;5815;7814;7714;7614;7514;7414;7314;7214;...
7114;7014;5014;5114;5214;5314;5414;5514;5614;5714;5814;...
7813;7713;7613;7513;7413;7313;7213;7113;7013;5013;5113;...
5213;5313;5413;5513;5613;5713;5813;7812;7712;7612;7512;...
7412;7312;7212;7112;7012;5012;5112;5212;5312;5412;5512;...
5612;5712;5812;7811;7711;7611;7511;7411;7311;7211;7111;...
7011;5011;5111;5211;5311;5411;5511;5611;5711;5811;7810;...
7710;7610;7510;7410;7310;7210;7110;7010;5010;5110;5210;...
5310;5410;5510;5610;5710;5810;7809;7709;7609;7509;7409;...
7309;7209;7109;7009;5009;5109;5209;5309;5409;5509;5609;...
5709;5809;7808;7708;7608;7508;7408;7308;7208;7108;7008;...
5008;5108;5208;5308;5408;5508;5608;5708;5808;7807;7707;...
7607;7507;7407;7307;7207;7107;7007;5007;5107;5207;5307;...
5407;5507;5607;5707;5807;7806;7706;7606;7506;7406;7306;...
7206;7106;7006;5006;5106;5206;5306;5406;5506;5606;5706;...
5806;7805;7705;7605;7505;7405;7305;7205;7105;7005;5005;...
5105;5205;5305;5405;5505;5605;5705;5805;7804;7704;7604;...
7504;7404;7304;7204;7104;7004;5004;5104;5204;5304;5404;...
5504;5604;5704;5804;7803;7703;7603;7503;7403;7303;7203;...
7103;7003;5003;5103;5203;5303;5403;5503;5603;5703;5803;...
7802;7702;7602;7502;7402;7302;7202;7102;7002;5002;5102;...
5202;5302;5402;5502;5602;5702;5802;7801;7701;7601;7501;...
7401;7301;7201;7101;7001;5001;5101;5201;5301;5401;5501;...
5601;5701;5801;7800;7700;7600;7500;7400;7300;7200;7100;...
7000;5000;5100;5200;5300;5400;5500;5600;5700;5800];

for i_hist_type=2:4
if i_hist_type==2
clim_dir=clim_ctd_;
clim_type="ctd_";
prefix=po_system_configuration.HISTORICAL_CTD_PREFIX;
end
if i_hist_type==4
clim_dir=clim_argo_;
clim_type="argo_";
prefix=po_system_configuration.HISTORICAL_ARGO_PREFIX;
end
if i_hist_type==3
clim_dir=clim_bottle_;
clim_type="bot_";
prefix=po_system_configuration.HISTORICAL_BOTTLE_PREFIX;
end

prefix_str=split(prefix,"/");
clim_dir=strcat(po_system_configuration.HISTORICAL_DIRECTORY,strjoin(prefix_str(1:end-1),"/"));

listing=dir(clim_dir);
n_files=size(listing,1);
for ifile=1:n_files
l_name=listing(ifile).name;
if strlength(l_name) < strlength(clim_type)
continue
end
if l_name(1:strlength(clim_type))== clim_type
box_no=strrep(l_name,clim_type,'');
box_no=strrep(box_no,'.mat','');
box_no=str2double(box_no);
i_box=(la_wmo_boxes(:,1)==box_no);
la_wmo_boxes(i_box,i_hist_type)=1;
end

end
end

save(la_wmo_boxes_file,'la_wmo_boxes')
30 changes: 18 additions & 12 deletions matlab_codes/get_region_ow.m
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@

function [ pa_grid_lat, pa_grid_long, pa_grid_dates ] = get_region_ow(pa_wmo_numbers, pa_float_name, po_config_data) ;
function [ pa_grid_lat, pa_grid_long, pa_grid_dates, pa_grid_box_id ] = get_region_ow(pa_wmo_numbers, pa_float_name, po_config_data)

% function [ pa_grid_lat, pa_grid_long, pa_grid_dates ] = get_region_ow( pa_wmo_numbers, pa_float_name, po_config_data ) ;
%
Expand All @@ -18,12 +18,16 @@
% Annie Wong, December 2007
% Breck Owens, December 2006
% C.Cabanes (2015) to save time: load only long, lat, dates and source data
% D.Dobler (DD), August 2024: to save more time later in retr_region_ow: also output
% corresponding boxes ids. Additional minor modifications: comment date_hist as it is unused, and
% correct indentation twice.

pa_grid_lat = [ ] ;
pa_grid_long = [ ] ;
pa_grid_dates = [ ] ;
pa_grid_lat = [ ] ;
pa_grid_long = [ ] ;
pa_grid_dates = [ ] ;
pa_grid_box_id = [ ] ;

[m,n]=size(pa_wmo_numbers);
[m,~]=size(pa_wmo_numbers);



Expand All @@ -40,7 +44,7 @@
%toc

not_use=[];
date_hist = changedates(lo_box_data.dates);
%date_hist = changedates(lo_box_data.dates);
lo_box_data.lat(not_use)=[];
lo_box_data.long(not_use)=[];
lo_box_data.dates(not_use)=[];
Expand All @@ -55,22 +59,24 @@
not_use=[];
for i=1:length(lo_box_data.lat)
profile=lo_box_data.source{i};
jj=findstr(profile,'_');
jj=findstr(profile,'_');
ref_float=profile(1:jj-1);
kk=findstr(pa_float_name, ref_float);
if(isempty(kk)==0)
not_use=[not_use,i];
not_use=[not_use,i];
end
end
lo_box_data.lat(not_use)=[];
lo_box_data.long(not_use)=[];
lo_box_data.dates(not_use)=[];
%-----------------------------------------------------------------------
end

pa_grid_lat = [ pa_grid_lat, lo_box_data.lat ] ;
pa_grid_long = [ pa_grid_long, lo_box_data.long ] ;
pa_grid_dates = [ pa_grid_dates, lo_box_data.dates ] ;

n_observations = size(lo_box_data.lat,2);
pa_grid_box_id = [ pa_grid_box_id, pa_wmo_numbers(ln_index,1)*ones(1,n_observations)] ;
pa_grid_lat = [ pa_grid_lat, lo_box_data.lat ] ;
pa_grid_long = [ pa_grid_long, lo_box_data.long ] ;
pa_grid_dates = [ pa_grid_dates, lo_box_data.dates ] ;
end
end
fclose('all');
Expand Down
11 changes: 9 additions & 2 deletions matlab_codes/map_data_grid.m
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,10 @@
% Bretherton, etal, 1975) is removed. The error estimate includes
% the contributions from both the mapping and the mean
% estimatation.
%
% Delphine Dobler (DD), August 2024:
% 3.3 - Performance: Inside map_data_grid, comment calculation
% of vdataerror as this is never used afterwards.
% **********************************************************************

[m,n]=size(v);
Expand Down Expand Up @@ -72,8 +76,11 @@
Cmd = s*covarxyt_pv(posdata, posdata, gs, ts, q, phi, pv);
vdata = Cmd*wght + vm;
% include error in mean (eqn 24 of Bretherton, etal)
vdataerror = repmat(sqrt(s-diag(Cmd*Cdd*Cmd') ...
+(1- sum(Cmd*Cdd,2)).^2/sum_Cdd),1,1);
%vdataerror = repmat(sqrt(s-diag(Cmd*Cdd*Cmd') ...
% +(1- sum(Cmd*Cdd,2)).^2/sum_Cdd),1,1);
% DD (2024/08-3.3): never used afterwards, save some time (Uncomment to
% compute)
vdataerror = [];

% map to posgrid -----

Expand Down
Loading