clc;clear all;close all;

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"
LAT=ncread(filename,'TLAT');lat=LAT(1,:);
LON=ncread(filename,'TLONG');lon=LON(:,1);
temp=ncread(filename,'TEMP');mask=temp;
mask(~isnan(mask))=1;mask(isnan(mask))=0; %1=ocean & nan=land 

pathi='....~/scratch/cesm2/archive/G_JRA55_test1/ocn/hist/' % pop2 model archive path

file_template = 'G_JRA55_test1.pop.h.%04d-%02d.nc';

nyra_start = 184; nyra_end = 366;

nx=320;ny=384;
n_months = (nyra_end - nyra_start + 1) * 12;
MLD_003=zeros(nx,ny,[(nyra_end-nyra_start+1)*12]);

for iyr = nyra_start:nyra_end
    for im = 1:12
    count = (iyr - nyra_start) * 12 + im;
    current_file = sprintf(file_template, iyr, im);
    full_path = [pathi  current_file];
    hmlx_dr=ncread(full_path,'HMXL_DR'); %cm --> m
    
    MLD_003(:,:,count)= hmlx_dr.*1e-2;
    end
end

%----------- study region ------------%
area=ncread(filename,'TAREA'); %(cm^2)
loc1=find(LON<=360-37.5 & LON>=360-55 & LAT<=57.5 & LAT>=47.5); 

for k=1:60
areamask{k}=area.*mask(:,:,k);    %2D
areamask_IN{k}=areamask{k}(loc1); %1D
sumarea(k)=squeeze(nansum(areamask_IN{k}));
end

for count=1:size(MLD_003,3)
mld_k{count} =MLD_003(:,:,count);
mld_IN{count}=mld_k{count}(loc1);

MLD_03(count)=nansum(nansum(mld_IN{count}.*areamask_IN{1}))/sumarea(1);
end

%------ annual mean ------%
for it = 1:size(MLD_003,3)/12
MLD_yr(it)=nanmean(MLD_03((it-1)*12+1:(it-1)*12+12));
end

%===== DJF ======%
ST=1;END=size(MLD_003,3); 
month12=ST+11:12:END;
month12=[(ST+0:12:END)' (ST+1:12:END)' (ST+11:12:END)']
MLD_DJF=nanmean(MLD_03(month12),2);

%===== March ======%
month03=ST+2:12:END;
MLD_Mar=MLD_03(month03);

patho='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/'
fn2=[patho 'MLD_4_6th_N.mat'];
save(fn2,'MLD_yr','MLD_DJF','MLD_Mar','MLD_003','-v7.3')





