clc;clear all;close all;
%===================================================
% OBS AMOC at 26.5N
%===================================================
%input:
% pathi ='/data3/yiwen/cesm_final_exp/plot_test/OBS_data/moc_transports.nc';
% pathi='/work/yiwen89419/cesm_POP/OBS_data/RAPID/2018/moc_transports.nc'
pathi='/work/yiwen89419/cesm_POP/OBS_data/RAPID/2022/moc_transports.nc'
time=ncread(pathi,'time');
baseDate = datetime(2004, 4, 1,00,00,00);
finalDate = baseDate + days(time);  % 將天數加到基準日期上
disp(finalDate);


%<overturning transport> size= 10535x1, unit='sv'
moc_mar_hc10 = ncread(pathi,'moc_mar_hc10'); 
year_start=2004; month_start=4; day_start=1;time_end=0;
% year_end=2018;month_end=9;day_end=30; time_end=0;
year_end=2022;month_end=1;day_end=15; time_end=0;

% calculation:
OBS_MOC_mon = [];
OBS_MOC_year = [];
mon_days = [31,28,31,30,31,30,31,31,30,31,30,31]; %from 1-12 month
mon_days2= [31,29,31,30,31,30,31,31,30,31,30,31];
years = [year_start:1:year_end];

countm =0;
countmdays=0;
countydays=0;
for iyr=1:length(years)
%-------------------------
if (mod(years(iyr),4))==0
days = mon_days2;
else
days = mon_days;
end
%-------------------------
if (years(iyr) == 2004)
ms=4;me=12;
totaldays = sum(days(ms:me));
elseif (years(iyr) == 2022) %end-year
% ms=1;me=8; %2004-2018/9
ms=1;me=1;   %2004-2022/2
totaldays = sum(days(ms:me));
else
ms=1;me=12;
totaldays = sum(days);
end
%------------------------
for im = ms:me
countm = countm+1;
%------------------------
%[years(iyr) im countmdays+1 countmdays+days(im)*2 days(im)*2]
tempA=[];tempA = moc_mar_hc10(countmdays+1:countmdays+days(im)*2);

OBS_MOC_mon(countm) = mean(tempA(~isnan(tempA)));
countmdays = countmdays + days(im)*2;
%------------------------
end
[years(iyr) countydays+1 countydays+totaldays*2 totaldays*2]
tempA=[];tempA = moc_mar_hc10(countydays+1:countydays+totaldays*2);
OBS_MOC_year(iyr) = mean(tempA(~isnan(tempA)));
countydays = countydays + totaldays*2;
end

%%
% output:
%======================================================
% save to mat
%======================================================
% patho ='/data3/yiwen/cesm_final_exp/plot_test/timeseries_result/';
patho='/work/yiwen89419/cesm_POP/OBS_data/2022/'

readme = {['OBS AMOC at 26.5N'];
          ['year 2004-2018']
          ['-----------------------------------------------'];
          ['OBS_MOC_mon: unit:SV'];
          ['OBS_MOC_year: unit:SV'];
          ['-----------------------------------------------'];
          [''];};

fn21 = [patho 'MOC_OBS_time.mat'];
save(fn21,'OBS_MOC_mon','OBS_MOC_year');

