clc;clear all;close all;
filename='/work1/yiwen89419/scratch/cesm2/archive/KPP_ICE_linear/G_JRA55_test1/ocn/hist/G_JRA55_test1.pop.h.0061-12.nc';
depth=ncread(filename,'moc_z');
lat_depth=ncread(filename,'lat_aux_grid');

%========== POP2 model output ===========%
pathi='....~/scratch/cesm2/archive/G_JRA55_test1/ocn/hist/' % pop2 model archive path
% pathi='....~/scratch/cesm2/archive/KPP_ICE_linear/G_JRA55_test1/ocn/hist/' % pop2 model archive path

nyra_start = 1; nyra_end = 61;

MOC_n1=zeros([395,61,(nyra_end-nyra_start+1)*12]);
count=0;
for iyr = nyra_start:nyra_end
 for im = 1:12
count = count + 1; %0~12
FNa=[pathi 'G_JRA55_test1.pop.h.' num2str(iyr,'%04d') '-' num2str(im,'%02d') '.nc'];
MOC=ncread(FNa,'MOC');

for i=1:length(lat_depth)
 for j=1:length(depth)
   %cal transport(1)Eulerian mean (2)Eddy-induce (3)Submeso
   %Region (1)=GMOC ,(2)=AMOC
   MOC_n1(i,j,count)=MOC_n1(i,j,count)+MOC(i,j,1,2);
 end
end

 end
end

% ---(2)--- %
%==== annual mean =====%
for it = 1:size(MOC_n1,3)/12
 MOC_yr(:,:,it)=nanmean(MOC_n1(:,:,(it-1)*12+1:(it-1)*12+12),3);
end

%===== every 10-yrs =====%
MOC_1st=squeeze(nanmean(MOC_yr(:,:,1:10),3));  %1958~1967
MOC_2nd=squeeze(nanmean(MOC_yr(:,:,11:20),3)); %1968~1977
MOC_3rd=squeeze(nanmean(MOC_yr(:,:,21:30),3)); %1978~1987
MOC_4th=squeeze(nanmean(MOC_yr(:,:,31:40),3)); %1988~1997
MOC_5th=squeeze(nanmean(MOC_yr(:,:,41:50),3)); %1998~2007
MOC_10=squeeze(nanmean(MOC_yr(:,:,52:61),3));  %last-10yrs

MOC_every10={MOC_1st MOC_2nd MOC_3rd MOC_4th MOC_5th MOC_10}
index={'1958_1967','1968_1977','1978_1987','1988_1997','1998_2007','2009_2018'}

patho='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/'
% patho='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/'

fn1= [patho 'MOC_spatial.mat'];
save(fn1,'MOC_every10','index','-v7.3')
