clc ;clear all;close all;
%===================================================
% perpare : lon, lat, dx, dy, dz
%===================================================
in1='/data3/yiwen/Timcom_exp/EXP_TIMCOM/change_wind_shr_45/GTJRAv2.10/ocn/hist/TOPO.nc'
ncid=netcdf.open(in1,'NC_NOWRITE');
varid = netcdf.inqVarID(ncid,'mask'); mask=netcdf.getVar(ncid,varid,'double'); % size=320x288x45
varid = netcdf.inqVarID(ncid,'lon'); lon=netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'lat'); lat=netcdf.getVar(ncid,varid,'double');
netcdf.close(ncid)

nx = size(mask,1); % 320
ny = size(mask,2); % 288
nz = size(mask,3); % 55

in2='/data2/summer2021/rundata_new.nc';
ncid=netcdf.open(in2,'NC_NOWRITE');
varid = netcdf.inqVarID(ncid,'dy'); dy=netcdf.getVar(ncid,varid,'double'); % size=304x1 (j0_global = 304)
varid = netcdf.inqVarID(ncid,'dx'); dx=netcdf.getVar(ncid,varid,'double'); % size=304x1 (j0_global = 304)
netcdf.close(ncid)


in3='/data2/summer2021/yzgrid_new.nc';
ncid=netcdf.open(in3,'NC_NOWRITE');
varid = netcdf.inqVarID(ncid,'odz'); odz=netcdf.getVar(ncid,varid,'double'); % size=45x1
netcdf.close(ncid)


j0_global = length(dy);   % 304
nghost = 2;
ngy0 = (j0_global-ny)/nghost; % 8
ny_size_0 = j0_global/ngy0;   % 38

n0=0;
for nn = 1:ngy0
 for yy = 2:ny_size_0-1
  n0 = n0+1;
  jj = (nn-1)*ny_size_0+yy;
 DY0(n0) = dy(jj); %size=1x288
  DX0(n0) = dx(jj);
 end
end

for k = 1:nz
 for j = 1:ny
  for i = 1:nx
    DX00(i,j) = DX0(j);
   DY00(i,j) = DY0(j);
  end
 end
end
%===================================================
% AMOC at 26.5N
%===================================================
%========EXP wind stress=======%  
pathia='/data3/yiwen/Timcom_exp/EXP_TIMCOM_renew/test01_control/ocn/'                            %CTAL
% pathia ='/data3/yiwen/Timcom_exp/EXP_TIMCOM_renew/change_wind_shr_45/GTJRAv2.10/ocn/hist/'     %change wind 45
% pathia ='/data3/yiwen/Timcom_exp/EXP_TIMCOM_renew/change_FWF_45/GTJRAv2.10/ocn/hist/'          %change FWF 45
% pathia ='/data3/yiwen/Timcom_exp/EXP_TIMCOM_renew/change_wind_shr_medius/GTJRAv2.10/ocn/hist/' %change wind medius[30-85N]

% pathia ='/data3/yiwen/Timcom_exp/change_wind_shr_NS/GTJRAv2.10/ocn/hist/' %change wind NS
% pathia ='/data3/yiwen/Timcom_exp/change_FWF_NS/GTJRAv2.10/ocn/hist/'
% %change FWF NS\
% pathia='/data3/yiwen/Timcom_exp/EXP_TIMCOM/corrected_data/'

nyra_start = 306; nyra_end = 366;

maxlatIN = find(26.2<lat & lat<26.8);
% maxlatIN = find(45.0<lat & lat<45.5);
maxlat = lat(maxlatIN);

W_lonIN = min(find(lon<360-15 & lon>360-80));
E_lonIN = max(find(lon<360-15 & lon>360-80));

% calculation:
AMOCa = zeros([nz,(nyra_end-nyra_start+1)*12]);

for iyr = nyra_start:nyra_end
for im = 1:12
count = (iyr-nyra_start)*12+im;
FNa = [pathia 'DATA_' num2str(iyr,'%04d') num2str(im,'%02d') '.nc'];
Va = ncread(FNa,'v-velocity');

for iNNZZ = nz:-1:1
  for k = nz:-1:iNNZZ
    for i = W_lonIN:E_lonIN  %Atlantic region
      AMOCa(iNNZZ,count) = AMOCa(iNNZZ,count) - mask(i,maxlatIN,k)*Va(i,maxlatIN,k)*DX00(i,maxlatIN)/odz(k); % unit:cm^3/s
    end % i
  end % k
end

end % im
end % iyr

% output:
AMOC_JRAc = max(AMOCa)/(10^12); % unit: SV=10^6 m^3/s



period = [nyra_start,nyra_end];

for it = 1:length(AMOC_JRAc)/12
SV_yr(it)=mean(AMOC_JRAc((it-1)*12+1:(it-1)*12+12));
end

%% ======================================================
% save to mat
%======================================================
patho='/data3/yiwen/Timcom_exp/matlab_code/timeseries_renew/'

fn21 = [patho 'change_wind_medius'];
save(fn21,'SV_yr','period')



