matlab计算水汽通量积分

AnneChen

全新手。
这段时间为了计算整层水汽通量，一会儿摸索Grads，一会儿用matlab。
Grads自带的vint用来做整层水汽通量积分是最好，但是后续我需要做距平，grads就有些不好用了。
由于技术有限，我尝试将Grads计算出来的水汽通量dat文件导入matlab中，可是维数不对，也不知道dat的数据是如何存储的。
因此，我又想用matlab写出类似vint的积分方程，从而利用它来计算水汽通量积分。
以下是我写的代码，不知道对不对，还请指教。

数据来自NCEP再分析数据
单位分别为：
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');

%计算各层qu，qv
q1=q(:,:,1:8,1:819);
u1=u(:,:,1:8,1:819);
v1=v(:,:,1:8,1:819);
        for i=1:8
            for j=1:819
                if q1(:,:,i,j)>=0 &  -70<u1(:,:,i,j)<125 & -70<v1(:,:,i,j)<70
                    qu(:,:,i,j)=q1(:,:,i,j).*u1(:,:,i,j);
                    qv(:,:,i,j)=q1(:,:,i,j).*v1(:,:,i,j);
                else
                    qu=NaN;
                    qv=NaN;
                end
            end
        end


%各网格点建立气压层数据，并用地表气压订正。
lev=ncread('shum.mon.mean.nc','level');
ps=ncread('pres.sfc.mon.mean.nc','pres');

lps=size(ps);
for k=1:lps(1)
    for kk=1:lps(2)
        for kkk=1:lps(3)
            pss=ps(k,kk,kkk);
            if pss>=0   
                lll=[pss;lev];
                dl=lll(1:end-1)-lll(2:end);
                dp(k,kk,:,kkk)=dl; %#ok<SAGROW>
            else
                lll=[1000;lev];
                dl=lll(1:end-1)-lll(2:end);
                dp(k,kk,:,kkk)=dl; %#ok<SAGROW>
            end
        end
    end
end
dp(dp<0)=0;  % 水汽积分是离散求和计算而得。因此这一步在于计算每一层层高。层高是下层减上层之间的差值。
                       % 由于高海拔地区如青藏高原地表气压可能小于地面气压数据，因此在用地表气压订正的时候，该地区可能出现负值层高。
                       % 这行代码用来置零。

%对整层水汽积分
tqup_vint=(100/9.8*sum(qu.*dp,3)/1000); %不知道为何vint里需要乘以100作为scale factor，而且vint函数对其他因子的单位标准是什么？还请指导。
                                                                     %乘以100的目的应该是将milibar 转换为100hpa。(感谢网友热心回答)
tqvp_vint=(100/9.8*sum(qv.*dp,3)/1000);   % unit: kg/(ms)
tqup=squeeze(tqup_vint);
tqvp=squeeze(tqvp_vint);

微醉的阳光

楼主，你最后选择sum的方式结果可能不太好，由于层数太少，导致累加与积分结果可能有些误差，建议采用trapz

qixianger

你好，我觉得你程序的积分不是从地面开始的，你每次的level除去最前面的地面气压其他的都还是从1000hPa开始的，所以有问题，我修改了一下，认为是这个样子的lev1=ncread('shum.mon.mean.nc','level');
ps=ncread('pres.mon.mean.nc','pres');

lps=size(ps);
for k=1:lps(1)
    for kk=1:lps(2)
        for kkk=1:lps(3)
            pss=ps(k,kk,kkk);
            lev = lev1;
            if pss>=0
                lev(lev>pss)=pss;
                lll=[pss;lev];
                dl=lll(1:end-1)-lll(2:end);
                dp(k,kk,:,kkk)=dl; %#ok<SAGROW>
            else
                lll=[1000;lev];
                dl=lll(1:end-1)-lll(2:end);
                dp(k,kk,:,kkk)=dl; %#ok<SAGROW>
            end
        end
    end
end

你可以算一下青藏高原的水汽通量结合文献，看一下你的。如果有什么不对的地方，欢迎指正。

zhchlue

可以使用积分trapz

AnneChen

zhchlue 发表于 2016-5-10 18:14
可以使用积分trapz

trapz（）函数是用梯形公式计算数值积分，对矩阵使用trapz（）函数时，相当于将行或列相邻元素相加除以2（即相邻元素的平均值）再求和。设A为m行n列矩阵，trapz（A，1）对列进行处理，输出为n个分量的行向量；trapz(A,2)对行进行处理，输出为m个分量的列向量。 所以这个是面积分函数？

kongfeng0824

水汽通量散度可能更有意义

kongfeng0824

这是第一次看到用MATLAB计算水汽通量的帖子。鼓励！

AnneChen

kongfeng0824 发表于 2016-5-12 13:25
水汽通量散度可能更有意义

一个是量，一个是变化量。两个一起更好。
不知道这样积分对不对呢，还请指教。

zhchlue

AnneChen 发表于 2016-5-10 18:58
trapz（）函数是用梯形公式计算数值积分，对矩阵使用trapz（）函数时，相当于将行或列相邻元素相加除以2 ...

是的，你理解的很对

kongfeng0824

请问用matlab解决了水汽通量及其散度的算法了吗？可以写一篇帖子分享一下吗？

玄德来了

话说再分析资料grib2文件里面没有各层的Q，还需要自己计算，NC里面有的？

AnneChen

kongfeng0824 发表于 2016-7-31 09:29
请问用matlab解决了水汽通量及其散度的算法了吗？可以写一篇帖子分享一下吗？

以上就是我计算的整层水汽通量算法。不过还没有写散度的问题。