clc;clear all;close all

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"

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/' 


file_template = 'G_JRA55_test1.pop.h.%04d-%02d.nc';
nyra_start = 1; nyra_end = 61;  

nx=320;ny=384;nz=60;
n_months = (nyra_end - nyra_start + 1) * 12;
V_vel=zeros([nx,ny,nz,n_months]);
U_vel=zeros([nx,ny,nz,n_months]);

for iyr = nyra_start:nyra_end
    for im = 1:12
    count = (iyr - nyra_start) * 12 + im;
    current_file = sprintf(file_template, iyr, im);
    full_path = [pathi  current_file];
    %Velocity in grid-XY direction%
    uvel = ncread(full_path,'UVEL');  %cm/s
    vvel = ncread(full_path,'VVEL');  %cm/s 
    rot  = ncread(full_path,'ANGLE'); %radians
    %cal true velocity %
    uPop =  uvel.*cos(rot) - vvel.*sin(rot);  
    vPop =  uvel.*sin(rot) + vvel.*cos(rot);
    %choose all data%
    U_vel(:,:,:,count)=uPop; 
    V_vel(:,:,:,count)=vPop;
    end
end
clearvars -except U_vel V_vel
%======= ANNUAL MEAN ========%
for it = 1:size(U_vel,4)/12
Uvel_yr(:,:,:,it)=nanmean(U_vel(:,:,:,(it-1)*12+1:(it-1)*12+12),4);
Vvel_yr(:,:,:,it)=nanmean(V_vel(:,:,:,(it-1)*12+1:(it-1)*12+12),4);
end

for ii=1:61
UVEL_yr{ii} = Uvel_yr(:,:,:,ii);
VVEL_yr{ii} = Vvel_yr(:,:,:,ii);
end


Vvel={Uvel_yr Vvel_yr}
%=====================================
%mean time
%=====================================
%===== every 10-yrs [lon*lat*depth]=====%
for ii=1:2
Vvel_1st{ii}=squeeze(nanmean(Vvel{ii}(:,:,:,1:10),4));  %1958~1967
Vvel_2nd{ii}=squeeze(nanmean(Vvel{ii}(:,:,:,11:20),4)); %1968~1977
Vvel_3rd{ii}=squeeze(nanmean(Vvel{ii}(:,:,:,21:30),4)); %1978~1987
Vvel_4th{ii}=squeeze(nanmean(Vvel{ii}(:,:,:,31:40),4)); %1988~1997
Vvel_5th{ii}=squeeze(nanmean(Vvel{ii}(:,:,:,41:50),4)); %1998~2007
Vvel_10{ii}=squeeze(nanmean(Vvel{ii}(:,:,:,52:61),4));  %last-10yrs
end

Vvel_every10={Vvel_1st Vvel_2nd Vvel_3rd Vvel_4th Vvel_5th Vvel_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 'ISO_Vel.mat'];
save(fn1,'Vvel_every10','UVEL_yr','VVEL_yr','index','-v7.3')
