tic,close all, clc, clear,format long

%============================================
% CLIMATE FIELD
%============================================
fixZlevel =  1;
in1='/work/yhtseng00/scratch/taiesm_timcom/archive/BT1850_adjustcloud_v2.10_shallow2_narrow_medit/ocn/hist/TOPO.nc'
%in1='/data2/giant/OCEAN_OUTPUT/CLOGIN/yht_GTIAF_QVOLa/TOPO.nc';
sclis = ncread(in1,'sclis');
tclis = ncread(in1,'tclis');
kb=ncread(in1,'kb')
dim=size(sclis)
 nx=dim(1); ny=dim(2); nz=dim(3); nm=dim(4);

mask_2 = double(ncread(in1,'mask'));
lat = ncread(in1,'lat');
lon = ncread(in1,'lon');
[lon_m,lat_m]=meshgrid(lon,lat);

CLIMATE_S=zeros(nx,ny);
CLIMATE_T=zeros(nx,ny);

for im=1:nm
 CLIMATE_S = CLIMATE_S + (sclis(:,:,fixZlevel,im).*mask(:,:,fixZlevel))./12;
 CLIMATE_T = CLIMATE_T + (tclis(:,:,fixZlevel,im).*mask(:,:,fixZlevel))./12;
end

%pathib='/data2/giant/OCEAN_OUTPUT/CLOGIN/yht_GTJRA_QVOLa/';
casename='BT1850_adjustcloud_v2.10_shallow2_narrow_medit'
pathib =['/work/yhtseng00/scratch/taiesm_timcom/archive/' casename '/ocn/hist/']
nyrb_start = 1; nyrb_end = 10;
nyear=nyrb_end-nyrb_start+1;
TEMP3Dd = zeros([nx,ny,nz]);
SALT3Dd = zeros([nx,ny,nz]);
count = 0;
for iyr = nyrb_start:nyrb_end
for im = 1:12
FNb = [pathib 'DATA_' num2str(iyr,'%04d') num2str(im,'%02d') '.nc'];
%vinfo = ncinfo(FN,'salinity'); dim=vinfo.Size; nx=dim(1); ny=dim(2); nz=dim(3);
TEMPd = ncread(FNb,'temperature');
SALTd = ncread(FNb,'salinity');
%SSTb = SSTb +TEMPb(:,:,fixZlevel).*mask(:,:,fixZlevel)./((nyrb_end-nyrb_start+1)*12);
%SSSb = SSSb +SALTb(:,:,fixZlevel).*mask(:,:,fixZlevel)./((nyrb_end-nyrb_start+1)*12);
TEMP3Dd=TEMP3Dd+TEMPd;
SALT3Dd=SALT3Dd+SALTd;
end
end
meanTEMPd=TEMP3Dd/(nyear*12);
meanSALTd=SALT3Dd/(nyear*12);