clc;clear all;close all

pathi1='/data3/yiwen/Timcom_exp/test01_control/ice/';                      %control run
pathi2='/data3/yiwen/Timcom_exp/change_wind_shr_45/GTJRAv2.10/ice/hist/';  %change wind 45 
pathi3='/data3/yiwen/Timcom_exp/change_wind_shr_NS/GTJRAv2.10/ice/hist/';  %change wind NS 
pathi4='/data3/yiwen/Timcom_exp/change_FWF45/GTJRAv2.10/ice/hist/';        %change FWF45 

% pathi4='/data3/yiwen/Timcom_exp/test01_control/ocn/';                    %control run
% pathi5='/data3/yiwen/Timcom_exp/change_wind_shr_45/GTJRAv2.10/ocn/hist/';%change wind 45 
% pathi6='/data3/yiwen/Timcom_exp/change_wind_shr_NS/GTJRAv2.10/ocn/hist/';%change wind NS 


nyra_start = 346; nyra_end = 366;  
nx=320;ny=288;

V_vel1=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);U_vel1=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);
V_vel2=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);U_vel2=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);
V_vel3=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);U_vel3=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);
V_vel4=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);U_vel4=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);

AREA_TD1=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);
AREA_TD2=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);
AREA_TD3=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);
AREA_TD4=zeros([nx,ny,(nyra_end-nyra_start+1)*12]);


count=0
for iyr = nyra_start:nyra_end
for im = 1:12
count = count +1;
FNb1 = [pathi1 'GTJRAv2.10.cice.h.' num2str(iyr,'%04d') '-' num2str(im,'%02d') '.nc']; %CTAL
FNb2 = [pathi2 'GTJRAv2.10.cice.h.' num2str(iyr,'%04d') '-' num2str(im,'%02d') '.nc']; %change 45
FNb3 = [pathi3 'GTJRAv2.10.cice.h.' num2str(iyr,'%04d') '-' num2str(im,'%02d') '.nc']; %change NS
FNb3 = [pathi3 'GTJRAv2.10.cice.h.' num2str(iyr,'%04d') '-' num2str(im,'%02d') '.nc']; %change NS
FNb3 = [pathi4 'GTJRAv2.10.cice.h.' num2str(iyr,'%04d') '-' num2str(im,'%02d') '.nc']; %change NS
% FNb4 = [pathi4 'DATA_' num2str(iyr,'%04d') num2str(im,'%02d') '.nc'];
% FNb5 = [pathi5 'DATA_' num2str(iyr,'%04d') num2str(im,'%02d') '.nc'];
% FNb6 = [pathi6 'DATA_' num2str(iyr,'%04d') num2str(im,'%02d') '.nc'];

uvel1 = ncread(FNb1,'uvel'); uvel2 = ncread(FNb2,'uvel'); 
uvel3 = ncread(FNb3,'uvel'); uvel4 = ncread(FNb4,'uvel');  %(cm/s)
vvel1 = ncread(FNb1,'vvel'); vvel2 = ncread(FNb2,'vvel'); 
vvel3 = ncread(FNb3,'vvel'); vvel4 = ncread(FNb4,'vvel');  %(cm/s)

area_td1= ncread(FNb1,'daidtd'); %area tendency dynamics [%/day]
area_td2= ncread(FNb2,'daidtd');
area_td3= ncread(FNb3,'daidtd');
area_td4= ncread(FNb4,'daidtd');

for i=1:nx
 for j=1:ny
   U_vel1(i,j,count)=U_vel1(i,j,count)+uvel1(i,j);
   V_vel1(i,j,count)=V_vel1(i,j,count)+vvel1(i,j);
   U_vel2(i,j,count)=U_vel2(i,j,count)+uvel2(i,j);
   V_vel2(i,j,count)=V_vel2(i,j,count)+vvel2(i,j);
   U_vel3(i,j,count)=U_vel3(i,j,count)+uvel3(i,j);
   V_vel3(i,j,count)=V_vel3(i,j,count)+vvel3(i,j);
   U_vel4(i,j,count)=U_vel4(i,j,count)+uvel4(i,j);
   V_vel4(i,j,count)=V_vel4(i,j,count)+vvel4(i,j);

   
   AREA_TD1(i,j,count)=AREA_TD1(i,j,count)+area_td1(i,j);
   AREA_TD2(i,j,count)=AREA_TD2(i,j,count)+area_td2(i,j);
   AREA_TD3(i,j,count)=AREA_TD3(i,j,count)+area_td3(i,j);
   AREA_TD4(i,j,count)=AREA_TD4(i,j,count)+area_td4(i,j);
 
 end
end


end
end

%=====================================
%mean time -->
%=====================================
Uvel_C=nanmean(U_vel1,3);Vvel_C=nanmean(V_vel1,3);
Uvel_45=nanmean(U_vel2,3);Vvel_45=nanmean(V_vel2,3);
Uvel_NS=nanmean(U_vel3,3);Vvel_NS=nanmean(V_vel3,3);
Uvel_FWF45=nanmean(U_vel4,3);Vvel_FWF45=nanmean(V_vel4,3);

%(%day)
AREATD_C=nanmean(AREA_TD1,3);
AREATD_45=nanmean(AREA_TD2,3);
AREATD_NS=nanmean(AREA_TD3,3);
AREATD_FWF45=nanmean(AREA_TD4,3);

%%
period = [nyra_start,nyra_end];
patho='/data3/yiwen/Timcom_exp/matlab_code/SIC/';
fn2 = [patho 'SIC_V&F'];
save(fn2,'Uvel_C','Vvel_C','Uvel_45','Vvel_45','Uvel_NS','Vvel_NS',...
         'Uvel_FWF45','Vvel_FWF45',...
         'AREATD_C','AREATD_45','AREATD_NS','AREATD_FWF45','period')


