clc;clear all;close all

file='/work/yiwen89419/cesm_POP/OBS_data/SAMBA/MOC_TotalAnomaly_and_constituents.asc'

%======= from 2009/3~2017/4 =======% 
[A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12] = textread(file,'%f%f%f%f%f%f%f%f%f%f%f%f',...
                                           'emptyvalue',NaN,'headerlines',59)
time=[A1,A2,A3,A4];

% Column 5: Total MOC anomaly (relative to record-length average of 14.7 Sv) 
mon_days = [31,28,31,30,31,30,31,31,30,31,30,31];
mon_days2 = [31,29,31,30,31,30,31,31,30,31,30,31];  % 闰年

for year09=1
% ====== (1)2009 ======== %
% MOC093=nanmean(A5(1:13));
% MOC094=nanmean(A5(14:14+mon_days(4)-1));
% MOC095=nanmean(A5(14+mon_days(4):14+mon_days(4)+mon_days(5)-1));
% MOC096=nanmean(A5(14+mon_days(4)+mon_days(5):14+mon_days(4)+mon_days(5)+mon_days(6)-1));
MOC09 = nan(1,10);  
start_idx = 1;
for i = 1:10-3
    if i == 1
        MOC09(i) = nanmean(A5(start_idx:start_idx + 13));
        start_idx = start_idx + 13;
    else
        month = 3 + i;  % month: 4 to 9
        MOC09(i) = nanmean(A5(start_idx:start_idx + mon_days(month) - 1));
        start_idx = start_idx + mon_days(month);
    end
end
end


% for year10=1
% ====== (2)2010~2016 ======== %
% SAMBA_1016=A5(289:2845);
% 
% year_start = 2010;year_end = 2016;
% years = year_start:year_end;
% 
% MOC10_16_mon = [];
% MOC10_16_year = [];
% 
% countm = 0;
% countmdays = 0;
% countydays = 0;
% 
% for iyr = 1:length(years)
%     % 判断是否为闰年
%     if mod(years(iyr), 4) == 0 && (mod(years(iyr), 100) ~= 0 || mod(years(iyr), 400) == 0)
%         days = mon_days2;
%     else
%         days = mon_days;
%     end
%     
%     if years(iyr) >= 2010
%     ms = 1; me = 12;  
%     totaldays = sum(days(ms:me));
%     end
% 
%     for im = ms:me
%         countm = countm + 1;
%         tempA = SAMBA_1016(countmdays+1:countmdays+days(im));
%         
%         MOC10_16_mon(countm) = nanmean(tempA(~isnan(tempA)));
%         countmdays = countmdays + days(im);
%     end
%     
%     tempA = SAMBA_1016(countydays+1:countydays+totaldays);
%     MOC10_16_year(iyr) = nanmean(tempA(~isnan(tempA)));
%     countydays = countydays + totaldays;
% end
% end


for year10=1
SAMBA_1016 = A5(289:2845);

year_start = 2010;
year_end = 2016;
years = year_start:year_end;

MOC10_16_mon = [];
MOC10_16_year = [];

countm = 0;
countmdays = 0;
countydays = 0;

for iyr = 1:length(years)
    % 判断是否为闰年
    if mod(years(iyr), 4) == 0 && (mod(years(iyr), 100) ~= 0 || mod(years(iyr), 400) == 0)
        days = mon_days2;
    else
        days = mon_days;
    end

    % 设置月份范围
    ms = 1; 
    me = 12;
    totaldays = sum(days(ms:me));

    for im = ms:me
        countm = countm + 1;
        start_idx = countmdays + 1;
        end_idx = countmdays + days(im);

        if end_idx > length(SAMBA_1016)
            warning('索引超出数据范围');
            break;
        end

        tempA = SAMBA_1016(start_idx:end_idx);

        % 调试信息
        fprintf('Year: %d, Month: %d, Start Index: %d, End Index: %d\n', ...
                years(iyr), im, start_idx, end_idx);

        MOC10_16_mon(countm) = nanmean(tempA);
        countmdays = countmdays + days(im);
    end

    start_idx = countydays + 1;
    end_idx = countydays + totaldays;

    if end_idx > length(SAMBA_1016)
        warning('索引超出数据范围');
        break;
    end

    tempA = SAMBA_1016(start_idx:end_idx);
    
    % 调试信息
    fprintf('Year: %d, Start Index: %d, End Index: %d\n', ...
            years(iyr), start_idx, end_idx);

    MOC10_16_year(iyr) = nanmean(tempA);
    countydays = countydays + totaldays;
end

end


for year17=1
% ====== (3)2017/4 ======== %
SAMBA_17=A5(2846:end);
MOC17 = nan(1,4);  
start_idx = 1;
for i = 1:4
    if i == 4
        MOC17(i) = nanmean(SAMBA_17(start_idx:start_idx + 28));
        start_idx = start_idx + 28;
    else
%         month = i;  % month: 1 to 3
        MOC17(i) = nanmean(SAMBA_17(start_idx:start_idx + mon_days(i) - 1));
        start_idx = start_idx + mon_days(i);
    end
end

end


%(A) monthly data

SAMBA_mon=[MOC09,MOC10_16_mon,MOC17]
SAMBA_year=[nanmean(MOC09),MOC10_16_year,nanmean(MOC17)]


patho='/work/yiwen89419/cesm_POP/OBS_data/SAMBA/2017/'

readme = {['OBS AMOC at 34.5S'];
          ['year 2009-2017']
          ['-----------------------------------------------'];
          ['OBS_MOC_mon: unit:SV'];
          ['OBS_MOC_year: unit:SV'];
          ['-----------------------------------------------'];
          [''];};

fn21 = [patho 'SAMBA_2017.mat'];
save(fn21,'SAMBA_mon','SAMBA_year');

%%

mon_days = [31,28,31,30,31,30,31,31,30,31,30,31];
mon_days2 = [31,29,31,30,31,30,31,31,30,31,30,31]; % 闰年







