clc;clear all;close all;

load('/data3/yiwen/cesm_final_exp/plot_test/check_diff_cycle/basin_mask.mat')
maskR2(isnan(maskR2))=0;

% in1='/data3/yiwen/cesm_final_exp/GTJRAv2.10c_v2/control_run/hist/TOPO.nc'
ncTOPO ='/data3/yiwen/cesm_final_exp/GTJRAv2.10c_v2/control_run/hist/TOPO.nc';
ncid=netcdf.open(ncTOPO,'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');
zlevel=ncread(ncTOPO,'lev_c');
zlevel=zlevel/100;  % turn cm tSo m
netcdf.close(ncid)

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

ncGRID = '/data2/summer2021/rundata_new.nc';
ncid=netcdf.open(ncGRID,'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)



ncGRIDZ = '/data2/summer2021/yzgrid.nc';
ncid=netcdf.open(ncGRIDZ,'NC_NOWRITE');
varid = netcdf.inqVarID(ncid,'odz'); odz=netcdf.getVar(ncid,varid,'double'); % size=45x1
netcdf.close(ncid);
j0_global = length(dy);   % 304

%%%%%%%%%% OBS data %%%%%%%
%gnid-Cell Area for Ocean Variables
ncOBS = [ '/data3/yiwen/cesm_final_exp/plot_test/OBS_data/areacello_Ofx_TaiESM1-TIMCOM_omip2_r1i1p2f1_gr.nc'];
%Monthly version of HadISST sea ice component
ncSIC =[ '/data3/yiwen/cesm_final_exp/plot_test/OBS_data/HadISST_ice.nc'];

%%
% pathib ='/data3/yiwen/cesm_final_exp/GTJRAv2.10c_v2/control_run/hist/' %control run
% pathib ='/data3/yiwen/cesm_final_exp/GTJRAv2.10c_v2/change_wind_v3/' %*1.3[30-85N]
pathib ='/data3/yiwen/cesm_final_exp/GTJRAv2.10c_v2/change_wind_v4/' %*0.1[30-85N]

% pathJRA = '/data3/yiwen/cesm_final_exp/ocn_GT.JRA.2.10_o3/hist/';


JRA_st = 230;   JRA_en = 244;    lbl='2004to2018';

jraperiod = [JRA_st,JRA_en];
JRAtot = (JRA_en-JRA_st+1)*12;
flistb = strings(1,JRAtot);

tic
count = 1;
for iyr = JRA_st:JRA_en
    for im = 1:12
        fnameb = [pathib,'DATA_' num2str(iyr,'%04d') num2str(im,'%02d') '.nc'];
        disp(fnameb);

        flistb(count) = fnameb;
        count = count +1;
    end
end

%%
%============================================
% CLIMATE FIELD
%============================================
fixZlevel =  1;

%load ([pathMAT,'/GSHHS_COAST220HL1.mat']);
ncid=netcdf.open(ncTOPO,'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');
zlevel=ncread(ncTOPO,'lev_c');
zlevel=zlevel/100  % turn cm to m
sclis = ncread(ncTOPO,'sclis');
tclis = ncread(ncTOPO,'tclis');
netcdf.close(ncid);

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

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

%%
%vinfo = ncinfo(FNa,'salinity'); dim=vinfo.Size; nx=dim(1); ny=dim(2); nz=dim(3);
for fc = 1:JRAtot
        disp(flistb(fc));
        TEMPb = ncread(flistb(fc),'temperature');
        SALTb = ncread(flistb(fc),'salinity');
        for z=1:nz
            %calculate potential temperature
            meanT0b(:,:,z)=Potentmp(SALTb(:,:,z),TEMPb(:,:,z),zlevel(z)*-1.023*0.98,0);
            %calcualte potential density
            rho_b(:,:,z)=rho_eos(squeeze(meanT0b(:,:,z)),squeeze(SALTb(:,:,z)),0);
        end

        %calculate using potential density 0.125
        for i=1:size(rho_b,1)
            for j=1:size(rho_b,2)
                A=min(find(abs(rho_b(i,j,2:end)-rho_b(i,j,1))>0.03));  %  A'=A+1 layer which rho>0.125
                
                if isempty(A)==1;MLDb(i,j,fc)=nan;   % no MLD
                elseif A==1;MLDb(i,j,fc)=zlevel(1);  % MLD is first layer   A'=2
                elseif A==kb(i,j)-1;MLDb(i,j,fc)=zlevel(A+1);else  % all the way to thetopography
                    mlo1=A; mlo2=A+1;depth=[zlevel(mlo1):0.1:zlevel(mlo2)];  %A=2 ,3rd layer which rho >0.125, interpolate between 2 and 3
                    rho_ainterp=interp1(zlevel(mlo1:mlo2),squeeze(rho_b(i,j,mlo1:mlo2)),depth,'linear','extrap');
                    B=min(find(abs(rho_ainterp(:)-rho_b(i,j,1))>0.03));
                    
                    if isempty(B)==1;MLDb(i,j,fc)=zlevel(A+1);else
                        MLDb(i,j,fc)=depth(B);
                    end %if
                end %if
            end%j
        end%i
end%fc

meanMLDb=mean(MLDb(:,:,:),3);   %mean for 230-244(monthly)

%%
%========================================
%save to mat
%========================================
patho='/data3/yiwen/cesm_final_exp/plot_test/check_diff_cycle/WIND_EXP/MLD/';

%%% MLD data %%%
fn21 = [patho, 'MLD_003_v4.mat'];
save(fn21,'lon','lat','mask','MLDb','meanMLDb')











