clear

%gcmface
%if(isdeployed==false)
%p = genpath('/home6/hzhang1/matlab/MITgcm_con2/gcmfaces'); addpath(p);
%end
nx=2160;ny=nx*13;
siz=[nx ny];
dirGrid='/nobackup/hzhang1/pub/llc2160/grid/';
nF=5;fileFormat='compact';
%global mygrid
grid_load(dirGrid,nF,fileFormat,2,1);
%

nx=2160; ny=nx*13; nz=90;
ss=nz/6; %6 segments

pp='/nobackupp11/dmenemen/DYAMOND/c1440_llc2160/mit_output/';

pout='UTVT/';
fldU='U'; fldV='V'; fldT='Theta';
uu=zeros([nx ny ss],'single');
vv=zeros([nx ny ss],'single');
tt=zeros([nx ny ss],'single');

%template:
TT=convert2gcmfaces(tt);
FLD=exch_T_N(TT); fldTu=TT; fldTv=TT;
for iF=1:5
   tmpA=FLD{iF}(2:end-1,2:end-1,:);   
   tmpB=FLD{iF}(1:end-2,2:end-1,:);
   fldTu{iF}=reshape(nanmean([tmpA(:) tmpB(:)],2),size(tmpA));
   tmpB=FLD{iF}(2:end-1,1:end-2,:);
   fldTv{iF}=reshape(nanmean([tmpA(:) tmpB(:)],2),size(tmpA));
end

%output
%yearly: 1 map
ut_yr_mean=zeros([nx ny ss],'single');
vt_yr_mean=zeros([nx ny ss],'single');
%monthly: 12 map
ut_mn_mean=zeros([nx ny ss],'single');
vt_mn_mean=zeros([nx ny ss],'single');
%daily: 365 map
ut_dy_mean=zeros([nx ny ss],'single');
vt_dy_mean=zeros([nx ny ss],'single');

%time
t0 = datenum(2020,1,19,21,0,0);        deltaT = 45;
ts1=(datenum(2020,3,1)-t0)*86400/deltaT+3600/deltaT;
ts2=(datenum(2021,3,1)-t0)*86400/deltaT;
ts=ts1:3600/deltaT:ts2; %length(ts)==365*24

k=1;
seg=(k-1)*ss+(1:ss);

ut_yr_mean=zeros([nx ny ss],'single');
vt_yr_mean=zeros([nx ny ss],'single');

yr=2020;
for mn=3:14 %2020/03 to 2021/02
	ut_mn_mean=zeros([nx ny ss],'single');
	vt_mn_mean=zeros([nx ny ss],'single');

	days=datenum(yr,mn+1,1)-datenum(yr,mn,1);
	for dy=1:days
			ut_dy_mean=zeros([nx ny ss],'single');
			vt_dy_mean=zeros([nx ny ss],'single');
		for hr=1:24
			t=( datenum(yr,mn,dy,hr,0,0)-t0 )*86400/deltaT;

			fn=[pp fldU '/' fldU '.' myint2str(t,10) '.data']; 
			uu = readbin(fn,[nx ny ss],1,'real*4', k -1);
			fn=[pp fldV '/' fldV '.' myint2str(t,10) '.data']; 
			vv = readbin(fn,[nx ny ss],1,'real*4', k -1);
			fn=[pp fldT '/' fldT '.' myint2str(t,10) '.data']; 
			tt = readbin(fn,[nx ny ss],1,'real*4', k -1);

TT=convert2gcmfaces(tt);
FLD=exch_T_N(TT); fldTu=TT; fldTv=TT;
for iF=1:5
   tmpA=FLD{iF}(2:end-1,2:end-1,:);   
   tmpB=FLD{iF}(1:end-2,2:end-1,:);
   fldTu{iF}=reshape(nanmean([tmpA(:) tmpB(:)],2),size(tmpA));
   tmpB=FLD{iF}(2:end-1,1:end-2,:);
   fldTv{iF}=reshape(nanmean([tmpA(:) tmpB(:)],2),size(tmpA));
end
			tt=single(convert2gcmfaces(fldTu));
			ut_dy_mean = ut_dy_mean + uu.*tt;
			ut_yr_mean = ut_yr_mean + uu.*tt;
			tt=single(convert2gcmfaces(fldTv));
			vt_dy_mean = vt_dy_mean + vv.*tt;
			vt_yr_mean = vt_yr_mean + vv.*tt;
		end %hr			
		%daily mean:
			ut_dy_mean = ut_dy_mean/24;
			vt_dy_mean = vt_dy_mean/24;
		f_day=datestr(t0 + t*deltaT/86400-1,'yyyymmdd');
		disp(f_day)
		writebin([pout 'UT_' f_day '_k' myint2str(k,2) '.bin'],ut_dy_mean)
		writebin([pout 'VT_' f_day '_k' myint2str(k,2) '.bin'],vt_dy_mean)

		ut_mn_mean  = ut_mn_mean  + ut_dy_mean;
		vt_mn_mean  = vt_mn_mean  + vt_dy_mean;
	end %dy		
	%monthly mean:
		ut_mn_mean = ut_mn_mean/days;
		vt_mn_mean = vt_mn_mean/days;
	f_mon=datestr(t0 + t*deltaT/86400-1,'yyyymm');
	writebin([pout 'UT_' f_mon '_k' myint2str(k,2) '.bin'],ut_mn_mean)
	writebin([pout 'VT_' f_mon '_k' myint2str(k,2) '.bin'],vt_mn_mean)
end %mn

%yearly mean
ut_yr_mean = ut_yr_mean/length(ts);
vt_yr_mean = vt_yr_mean/length(ts);
%save
writebin([pout 'UT_year_k' myint2str(k,2) '.bin'],ut_yr_mean);
writebin([pout 'VT_year_k' myint2str(k,2) '.bin'],vt_yr_mean);