用matlab求水汽通量散度

wmh951022

在搜程序的时候看见有人这么写来求水汽通量散度，但是程序运行到这一步的时候显示未定义函数或变量。请问，这是我的matlab里面缺少这个函数，还是什么问题呢？初学者，劳烦大家解答我的疑惑了。谢谢！

wmh951022

wmh951022 发表于 2020-11-11 21:49
clear all
clc

我也不是大气专业的，只是在网上找到了一个版本，然后根据我自己的需要尝试弄得，不是特别懂。找到的程序，贴在这里，供大家参考

二爷名声在外

是

gyq1991

你的水汽通量散度弄出来了吗，能分享吗？

wmh951022

gyq1991 发表于 2020-8-26 11:02
你的水汽通量散度弄出来了吗，能分享吗？

可以，你加我qq1057110485，我发给你

bibibol

希望楼主可以将程序分享出来~

Chwithin

希望楼主可以将程序分享出来~

wmh951022

clear all
clc

% q       g/kg
% u       m/s
% v       m/s
% p      millibar


q=ncread('shum.mon.mean.nc','shum');
u=ncread('uwnd.mon.mean.nc','uwnd');
v=ncread('vwnd.mon.mean.nc','vwnd');

% ncdisp('shum.mon.mean.nc','shum');
% ncdisp('uwnd.mon.mean.nc','uwnd');
% ncdisp('vwnd.mon.mean.nc','vwnd');
% ncdisp('pres.sfc.mon.mean.nc','pres');

lev=ncread('shum.mon.mean.nc','level');
ps0=ncread('pres.sfc.mon.mean.nc','pres');

Ofnm = 'E:pres.sfc.mon.mean.nc';
glon=ncread(Ofnm,'lon');
glat=ncread(Ofnm,'lat');

[Wlat,Wlon]=meshgrid(glat, glon); % ？
ind1= [157:168,229:240,373:384,481:492,565:576]; % low:1961、1967、1979、1988、1995(数据从1948年1月开始)


% % % % %
q1=q(:,:,1:8,ind1); % 1:8是1000hpa-300hpa
u1=u(:,:,1:8,ind1);
v1=v(:,:,1:8,ind1);
ps1=ps0(:,:,ind1);

% % % % %
qs1 = zeros(144,73,8,15);% 5年乘以3个月（6、7、8）
us1 = zeros(144,73,8,15);
vs1 = zeros(144,73,8,15);
ps11 = zeros(144,73,15);

% % % % %
i1 = 0;
mi1 = 0;
for yeari1 = 1:5
    for monthi1 = 1:12
        mi1 = mi1 + 1;
        if (monthi1==6)||(monthi1==7)||(monthi1==8)
            i1 = i1 +1;
            qs1(:,:,:,i1) = q1(:,:,:,mi1);
            us1(:,:,:,i1) = u1(:,:,:,mi1);
            vs1(:,:,:,i1) = v1(:,:,:,mi1);
            ps11(:,:,i1)  = ps1(:,:,mi1);

        end

    end
end

% % % % %
qu1 = zeros(144,73,8,15);
qv1 = zeros(144,73,8,15);

% % % % %
for i1=1:8
    for j1=1:15
        for y1 = 1:73
            for x1 = 1:144
                if (qs1(x1,y1,i1,j1)>=0) &&  (-70<us1(x1,y1,i1,j1)<125) && (-70<vs1(x1,y1,i1,j1)<70)
                    qu1(x1,y1,i1,j1)=qs1(x1,y1,i1,j1)*us1(x1,y1,i1,j1);
                    qv1(x1,y1,i1,j1)=qs1(x1,y1,i1,j1)*vs1(x1,y1,i1,j1);
                else
                    qu1(x1,y1,i1,j1)=NaN;
                    qv1(x1,y1,i1,j1)=NaN;
                end
            end
        end
    end
end

% 各网格点建立气压层数据，并用地表气压订正
% lev=ncread('shum.mon.mean.nc','level');
% ps=ncread('pres.sfc.mon.mean.nc','pres');
lps1=size(ps11);
for k=1:lps1(1)
    for kk=1:lps1(2)
        for kkk=1:lps1(3)
            pss1=ps11(k,kk,kkk);
            if pss1>=0   
                lll=[pss1;lev];
                dl=lll(1:end-1)-lll(2:end);
                dp1(k,kk,:,kkk)=dl;
            else
                lll=[1000;lev];
                dl=lll(1:end-1)-lll(2:end);
                dp1(k,kk,:,kkk)=dl;
            end
        end
    end
end
dp1(dp1<0)=0;

% 水汽积分是离散求和计算而得。因此这一步在于计算每一层层高。层高是下层减上层之间的差值。
% 由于高海拔地区如青藏高原地表气压可能小于地面气压数据，因此在用地表气压订正的时候，该地区可能出现负值层高。
% 这行代码用来置零。

% % % % %
tqup_vint1=(100/9.8*sum(qu1.*dp1,3)/1000);                                                               
tqvp_vint1=(100/9.8*sum(qv1.*dp1,3)/1000);   % unit: kg/(ms)


tqup1=squeeze(tqup_vint1);
tqvp1=squeeze(tqvp_vint1);
%  tqup111 = sum(tqup1,3)/5;


tqup_ind1 = mean(tqup1,3); % 3？

tqvp_ind1 = mean(tqvp1,3);

puck66

谢谢楼主分享好物！run一下试试

duruyi

太厉害了，一直想看matlab程序