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

mask(mask==0)=nan;
for k = 1:nz %
 for j = 1:ny
  for i = 1:nx
    DX00(i,j) = DX0(j);
    DY00(i,j) = DY0(j);
    VOLUME(i,j,k) = DX00(i,j)*DY00(i,j)*mask(i,j,k)*(1/odz(k));
  end
 end
end
%choose north Atlantic region%
W_lonIN = min(find(lon<360-15 & lon>360-80));
E_lonIN = max(find(lon<360-15 & lon>360-80));
lon_loc=W_lonIN:E_lonIN;
% lon_loc=1:320;
lat_loc=235:263;         %45N-65N
VOLUME_At=VOLUME(lon_loc,lat_loc,:);
SURFACE_VOLUME1=nansum(nansum(VOLUME_At(:,:,1)));


pathib =['/data3/yiwen/Timcom_exp/EXP_TIMCOM_renew/test01_control/ocn/']
nyrb_start = 346; nyrb_end = 366;
% calculation:
surfQDOTb = zeros([nx,ny,(nyrb_end-nyrb_start+1)*12]);
surfNWaterb= zeros([nx,ny,(nyrb_end-nyrb_start+1)*12]);
surfSUMSb= zeros([nx,ny,(nyrb_end-nyrb_start+1)*12]);
surfSUMTb= zeros([nx,ny,(nyrb_end-nyrb_start+1)*12]);
count = 0;
for iyr = nyrb_start:nyrb_end
for im = 1:12
count = count +1
FNb = [pathib 'DATA_' num2str(iyr,'%04d') num2str(im,'%02d') '.nc'];
qdot = ncread(FNb,'qdot');qdot2 =ncread(FNb,'qdot2');
Evap=ncread(FNb,'evapo');Precip=ncread(FNb,'prec');
TEMP = ncread(FNb,'temperature');SALT = ncread(FNb,'salinity');
    for i =1:nx
       for j = 1:ny
        surfQDOTb(i,j,count) = surfQDOTb(i,j,count) + VOLUME(i,j,1)*(qdot(i,j)+qdot2(i,j))*mask(i,j,1);
        surfNWaterb(i,j,count)= surfNWaterb(i,j,count) + VOLUME(i,j,1)*(-Evap(i,j)-Precip(i,j))*mask(i,j,1);
        
        surfSUMSb(i,j,count) = surfSUMSb(i,j,count) + VOLUME(i,j,1)*SALT(i,j,1)*mask(i,j,1);
        surfSUMTb(i,j,count) = surfSUMTb(i,j,count) + VOLUME(i,j,1)*TEMP(i,j,1)*mask(i,j,1);
       end % j
    end % i
end % im
end % iyr


% output:
SurfQDOTb = squeeze(nansum(nansum(surfQDOTb(lon_loc,lat_loc,:))))./SURFACE_VOLUME1;
SurfNWaterb= squeeze(nansum(nansum(surfNWaterb(lon_loc,lat_loc,:))))./SURFACE_VOLUME1;
SurfSUMSb= squeeze(nansum(nansum(surfSUMSb(lon_loc,lat_loc,:))))./SURFACE_VOLUME1;
SurfSUMTb= squeeze(nansum(nansum(surfSUMTb(lon_loc,lat_loc,:))))./SURFACE_VOLUME1;


rho0=1025   %kg/m3
Cp=4000     %J/kg
for count=1:size(surfQDOTb,2)
thermal(count)= SurfQDOTb(count)/Cp/rho0*(0.049762-0.0030063*SurfSUMSb(count)+...
                2*SurfSUMTb(count)*(-0.0075911+0.000037*SurfSUMSb(count))+3*SurfSUMTb(count)^2*0.000035);           
           
haline(count)=(0.8056-0.0030063*SurfSUMTb(count)+...
               0.000037*SurfSUMTb(count)^2)*SurfSUMSb(count)/(1-SurfSUMSb(count))*SurfNWaterb(count);           
end
D=thermal+haline

for iyr = 1:length(D)/12
 DD (iyr)= mean(D((iyr-1)*12+1:(iyr-1)*12+12));
 DT (iyr)= mean(thermal((iyr-1)*12+1:(iyr-1)*12+12));
 DH (iyr)= mean(haline((iyr-1)*12+1:(iyr-1)*12+12));
end

% plot(DD)
%%
set(0,'defaultfigurecolor','w')
plot([1998:2018],DD,'color',[0/255,0/255,205/255],...
    'linestyle','-','marker','none','MarkerSize',35,'MarkerFaceColor',...
     [0/255,0/255,205/255],'Linewidth',2);hold on
plot([1998:2018],DT); 
plot([1998:2018],DH,'linestyle','--','marker','none','MarkerSize',35); 
legend('Densiy total','thermal','haline')
set(gca,'fontsize',12);
set(gca,'XMinorTick','on','YMinorTick','on')
set(gca,'TickDir','out')
ytickformat('%.2f')
% ylabel({'qdot +qdot2'},'color','k','Fontsize',12)
set(gca, {'XColor', 'YColor'}, {'k','k'});
grid on;
axis([1998 2018 -5*1e-6 20*1e-6])




