混合层深度

蓬蒿人

参考的文献是：

[1] Huang P Q, Lu Y Z, Zhou S Q. An objective method for determining ocean mixed layer depth with applications to WOCE data[J]. Journal of Atmospheric and Oceanic Technology, 2018, 35(3): 441-458.

仅供参考，中文也是我翻译的，不一定说就准确哦。

发帖子目的也只是因为我差金币了。

蓬蒿人

ncdisp("E:\Data\soda3.4.2_mn_ocean_reg_2017.nc");

%读取盐度：实用盐度

salt=ncread("E:\Data\soda3.4.2_mn_ocean_reg_2017.nc","salt");

%读取经纬度、深度及时间

lon=ncread("E:\Data\soda3.4.2_mn_ocean_reg_2017.nc","xt_ocean");

lat=ncread("E:\Data\soda3.4.2_mn_ocean_reg_2017.nc","yt_ocean");

depth=ncread("E:\Data\soda3.4.2_mn_ocean_reg_2017.nc","st_ocean");

time=ncread("E:\Data\soda3.4.2_mn_ocean_reg_2017.nc","time");

%读取温度，此时温度的单位为位温

temp=ncread("E:\Data\soda3.4.2_mn_ocean_reg_2017.nc","temp");

%截取南海附近的经纬度及1月份的数据；

lon_need=find(lon>=100&lon<=120);

lat_need=find(lat>=4&lat<=21);

pt_need=temp(lon_need,lat_need,:,1);

salt_need=salt(lon_need,lat_need,:,1);

p=zeros(size(salt_need));

SA=zeros(size(salt_need));

CT=zeros(size(salt_need));

temp_need=zeros(size(salt_need));

for i=1:length(lon_need);

    for j=1:length(lat_need);

        for k=1:length(depth);

    p(i,j,k)=gsw_p_from_z(-depth(k),lat(lat_need(j)));

    SA(i,j,k)=gsw_SA_from_SP(salt_need(i,j,k),p(i,j,k),lon(i),lat(j));

    CT(i,j,k)=gsw_CT_from_pt(SA(i,j,k),pt_need(i,j,k));

    temp_need(i,j,k)=gsw_t_from_CT(SA(i,j,k),CT(i,j,k),p(i,j,k));

        end

    end

end

P=zeros(size(temp_need));

R=zeros(size(temp_need));

X=zeros(size(temp_need));

%求解标准差，最大偏差及相对偏差

for i=1:size(temp_need,1)

    for j=1:size(temp_need,2)

        for k=1:size(temp_need,3)

            %P是标准差

        P(i,j,k)=nanstd(temp_need(i,j,1:k));

        %R是最大偏差

        R(i,j,k)=max(temp_need(i,j,1:k))-min(temp_need(i,j,1:k));

        %X是相对偏差

        X(i,j,k)=P(i,j,k)/R(i,j,k);

        end

    end

end

%求解zn1

zn1=zeros(size(temp_need,1),size(temp_need,2));

deep_need=zeros(size(temp_need,1),size(temp_need,2));

for i=1:size(temp_need,1);

    for j=1:size(temp_need,2);

        for k=1:(size(temp_need,3)-1);

            if X(i,j,k)-X(i,j,k+1)<10^(-2);

                zn1(i,j)=zn1(i,j);

            else

                zn1(i,j)=depth(k+1);

            end

            if zn1(i,j)>5000;

                zn1(i,j)=nan;

            end

            %将zn1转化为深度

            if isnan(zn1(i,j));

            deep_need(i,j)=nan;

            else

                deep_need(i,j)=find(depth(:,1)==zn1(i,j));

            end

        end

    end

end

%计算zn2;

zn2=zeros(size(temp_need,1),size(temp_need,2));

for i=1:size(temp_need,1);

    for j=1:size(temp_need,2);

        for k=1:deep_need(i,j);

            if isnan(deep_need(i,j))

                zn2(i,j)=nan;

            else

                if ((P(i,j,k+1)-P(i,j,k))/(P(i,j,(deep_need(i,j)+1))-P(i,j,deep_need(i,j))))>0.5;

                    zn2(i,j)=zn2(i,j);

                else zn2(i,j)=depth(k);

                end

                if zn2(i,j)==0;

                    zn2(i,j)=nan;

                end

            end

        end

    end

end

figure(1);

hold on;

m_proj('Equidistant Cylindrical','latitudes',[4 21],'longtitudes',[100 120])

%海岸线

m_coast('patch',[0.90 0.90 0.90] ,'edgecolor',[0.38 0.38 0.38],'linewidth',1.0);

%网格

m_grid( 'box ' , 'on' , 'linest' , 'none ' , 'linewidth' ,0.1 , 'tickdir' , 'in' , ...

    'backcolor',[0.98 0.98 0.98],'ytick',[0:5:21], 'xtick',[100:5:120], 'fontsize',25);

[lat1,lon1]=meshgrid(lat(lat_need),lon(lon_need));

[c,h]=m_contourf(lon1,lat1,zn2(:,:),50,'linestyle','nFone');

clabel(c,h,[0:20:160]);

hBar=colorbar('location','eastoutside','fontsize',20);

title("南海及周边地区混合层的深度",'fontsize',20);