%%
clear all;
close all;
clc;

load /Volumes/Expansion/case_2007/0.figure_modified/ocean_currents/uvcurrent_winmean_eastnorth.mat;
lon_rho(lon_rho < 0) = lon_rho(lon_rho < 0) + 360;

%将深度m转化为sea pressure
for n=1:32
    p_rho(:,:,n) = gsw_p_from_z(-depth0(:,:,n),lat_rho);
end

%------ 特定transect上对应点的经纬度和计算得到的网格点面积
lon0=lon_rho(473,320:360);
lat0=lat_rho(473,320:360);
lon_sect=repmat(lon0,[32 1]);
lat_sect=repmat(lat0,[32 1]);
pd=squeeze(p_rho(473,320:360,:));%一维的压力值
dis=gsw_distance(lon_sect,lat_sect,pd');%单位是m

depthT=depth0;
depth_sect=squeeze(depthT(473,320:360,:))';
dz=[];
for m=1:31
    dz(m,:)=depth_sect(m,:)-depth_sect(m+1,:);%每层深度层对应的高度值增加值 此时第一层为底部
end
%不同深度层对应的网格的面积，此时用梯形面积公式（上底+下底）*h/2  将dis作为高h 不同点的dz作为上下底边
for n=1:40
    area0(:,n)=(dz(:,n)+dz(:,n+1))*0.5;
end
area=area0.*dis(1:31,:);
area=area';

%------ 特定transect上对应点的uv变量值
uu=squeeze(u(473,320:360,:,:));
vv=squeeze(v(473,320:360,:,:));

%%

for t=1:17
    uut=uu(:,:,t);
    vvt=vv(:,:,t);
    %对垂向网格的数据点进行处理，求网格四个点的平均值作为该网格中心点位置的值，再乘以对应的面积，则得到了export
    for nn=1:40
        %对每一层的两个点求平均值
        uy(nn,:)=(uut(nn,:)+uut(nn+1,:))/2;
        vy(nn,:)=(vvt(nn,:)+vvt(nn+1,:))/2;
    end
    
    for mm=1:31
        %对不同层上下的两个点求平均值,基于每层的平均值来求
        u0(:,mm)=(uy(:,mm)+uy(:,mm+1))/2;
        v0(:,mm)=(vy(:,mm)+vy(:,mm+1))/2;
    end
    
    
    for i=1:size(lon_sect,2)-1
        %将uv在每个网格上进行分解，求出垂直于面的分量 首先求出transect的斜率
        k=(lat_sect(1,1)-lat_sect(1,i+1))/(lon_sect(1,1)-lon_sect(1,i+1));
        alpha=atan(k); %负弧度表示负角 顺时针旋转为负角度
        
        uf=u0(i,:);
        vf=v0(i,:);
        areaf=area(i,:);
        
        up=uf*sin(abs(alpha));
        vp=vf*cos(abs(alpha));
        flux=areaf.*((up-vp));
        fluxt(i,:)=flux;%每一列的flux
    end
    
    index=find(u0(u0<0)); %coastal slope current/westward flow
    flux_total=nansum(fluxt(index));
    
    fw_transport(:,t)=flux_total;
end


save('/Volumes/Expansion/case_2007/0.figure_modified/0. volume_transport/ccvolumetransport_winmean_new.mat','fw_transport');
