clear

load_llc2160
tic
    gcmfaces_lines_zonal;
toc    

nx=2160; ny=nx*13; nz=90;
pp='/nobackupp11/dmenemen/DYAMOND/c1440_llc2160/mit_output/STATS/';

%This is very approxiamte in
%1) NO diffusion
%2) U^bar * T^bar instead of (U*T)^bar : waiting UTVT/
%3) T not @ U/V grid ==> FWHT_year0.mat

%year
fin=[pp 'U/U_year.bin'];
disp(fin)
UU=readbin(fin,[nx ny nz]);
fin=[pp 'V/V_year.bin'];
disp(fin)
VV=readbin(fin,[nx ny nz]);
fin=[pp 'THETA/Theta_year.bin'];
disp(fin)
TT=readbin(fin,[nx ny nz]);

%fldU=UVELMASS.*mygrid.mskW; fldV=VVELMASS.*mygrid.mskS;
%replaced by:
fldU=convert2gcmfaces(UU).*mygrid.hFacW; %UVEL==>UVELMASS
fldV=convert2gcmfaces(VV).*mygrid.hFacS;
fldU=fldU.*mygrid.mskW;
fldV=fldV.*mygrid.mskS;
fldT=convert2gcmfaces(TT).*mygrid.mskC;
FLD=exch_T_N(fldT); fldTu=fldT; fldTv=fldT;
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

tic
    gFW=1e-6*calc_MeridionalTransport(fldU,fldV,1);
    gHT=1e-15*4e6*calc_MeridionalTransport(fldU.*fldTu,fldV.*fldTv,1);
toc

tic
    mskC=v4_basin({'atlExt'}); mskC=mk3D(mskC,fldU);
    aFW=1e-6*calc_MeridionalTransport(fldU.*mskC,fldV.*mskC,1);
    aHT=1e-15*4e6*calc_MeridionalTransport(fldU.*fldTu.*mskC,fldV.*fldTv.*mskC,1);
toc

tic
    mskC=v4_basin({'pacExt','indExt'}); mskC=mk3D(mskC,fldU);
    pFW=1e-6*calc_MeridionalTransport(fldU.*mskC,fldV.*mskC,1);
    pHT=1e-15*4e6*calc_MeridionalTransport(fldU.*fldTu.*mskC,fldV.*fldTv.*mskC,1);
toc

save FWHT_year gFW aFW pFW gHT aHT pHT

%next monthly:
gFW=nan(179,  12); aFW=nan(179,  12); pFW=nan(179,  12);
gHT=nan(179,  12); aHT=nan(179,  12); pHT=nan(179,  12);

yr=2020;
for mn=3:14 %
tic
	f_mon=datestr(datenum(yr,mn,1),'yyyymm');
	fin=[pp 'U/U_' f_mon '.bin'];
	disp(fin)
	UU=readbin(fin,[nx ny nz]);
	fin=[pp 'V/V_' f_mon '.bin'];
	disp(fin)
	VV=readbin(fin,[nx ny nz]);
	fin=[pp 'THETA/Theta_' f_mon '.bin'];
	disp(fin)
	TT=readbin(fin,[nx ny nz]);

fldU=convert2gcmfaces(UU).*mygrid.hFacW; %UVEL==>UVELMASS
fldV=convert2gcmfaces(VV).*mygrid.hFacS;
fldU=fldU.*mygrid.mskW;
fldV=fldV.*mygrid.mskS;
fldT=convert2gcmfaces(TT).*mygrid.mskC;
FLD=exch_T_N(fldT); fldTu=fldT; fldTv=fldT;
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

    gFW(:,mn-3+1)=1e-6*calc_MeridionalTransport(fldU,fldV,1);
    gHT(:,mn-3+1)=1e-15*4e6*calc_MeridionalTransport(fldU.*fldTu,fldV.*fldTv,1);

    mskC=v4_basin({'atlExt'}); mskC=mk3D(mskC,fldU);
    aFW(:,mn-3+1)=1e-6*calc_MeridionalTransport(fldU.*mskC,fldV.*mskC,1);
    aHT(:,mn-3+1)=1e-15*4e6*calc_MeridionalTransport(fldU.*fldTu.*mskC,fldV.*fldTv.*mskC,1);

    mskC=v4_basin({'pacExt','indExt'}); mskC=mk3D(mskC,fldU);
    pFW(:,mn-3+1)=1e-6*calc_MeridionalTransport(fldU.*mskC,fldV.*mskC,1);
    pHT(:,mn-3+1)=1e-15*4e6*calc_MeridionalTransport(fldU.*fldTu.*mskC,fldV.*fldTv.*mskC,1);
toc
save FWHT_month gFW aFW pFW gHT aHT pHT
end

%save FWHT_month gFW aFW pFW gHT aHT pHT