请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 2865|回复: 14

matlab写涡度方程各项的代码问题

[复制链接]

新浪微博达人勋

发表于 2023-3-26 23:50:31 | 显示全部楼层 |阅读模式
9金钱

                               
登录/注册后可看大图

因为需要涡度方程,以家园里的分享grads代码想写一个matlab版本的,花了一天时间写了一个M项,但调试了半天不知为什么好几项都是空值,求出的V2项也仅仅是一个值而不是设想的一个数组,报错和代码如下,求大神解答。。。本科生跪谢。。。。。

% 读取数据
nc_path = '202007.nc';  % 数据路径
ncid = netcdf.open(nc_path, 'nowrite');
lon = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'longitude'));
lat = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'latitude'));
level = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'level'));
time = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'time'));
u = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'u'));
v = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'v'));
vo = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'vo'));
netcdf.close(ncid);
% 计算涡度方程各项
[lon1, lat1]=meshgrid(lon,lat);
for x = 1:length(lon1);  % 经度范围
     for y = 1:length(lat1);  % 纬度范围
        z0 = 1;  % 垂直层数
        t0 = 1;  % 时间步长
        fi = lat(y) * pi / 180;  % 纬度转弧度
        beta = 2 * cos(fi) * 7.29e-5;  % beta参数
        % 计算dvo_dx和dvo_dy
        % 计算dvo_dx和dvo_dy
        dx = diff(lon(x) * pi / 180);  % 计算经度间隔
        dy = diff(lat(y) * pi / 180);  % 计算纬度间隔
        dvo_dx = double(diff(vo(x, y, z0, t0), 1, 1)) ./ double(dx) * 6371000 * cos(fi);  % x方向上的vo梯度
        dvo_dy = double(diff(vo(x, y, z0, t0), 1, 2)) ./ double(dy) * 6371000;  % y方向上的vo梯度
        
        V1 = -double(u(x, y, z0, t0)) .* dvo_dx - double(v(x, y, z0, t0)) .* dvo_dy;
        V4 = double(v(x, y, z0, t0)) .* beta;
        V1_all(x, y) = V1;  % 将计算结果存储到数组中
        V4_all(x, y) = V4;
    end
end



                               
登录/注册后可看大图


                               
登录/注册后可看大图



最佳答案

查看完整内容

盲写的,我也不知道能不能解决你的问题,如果有报错就更好了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-3-26 23:50:32 | 显示全部楼层
  1. % 读取数据
  2. nc_path = '202007.nc';  % 数据路径
  3. ncid = netcdf.open(nc_path, 'nowrite');
  4. lon = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'longitude'));
  5. lat = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'latitude'));
  6. level = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'level'));
  7. time = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'time'));
  8. u = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'u'));
  9. v = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'v'));
  10. vo = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'vo'));
  11. netcdf.close(ncid);

  12. % 计算涡度方程各项
  13. [lon1, lat1] = meshgrid(lon, lat);
  14. V1_all = zeros(length(lon), length(lat));
  15. V4_all = zeros(length(lon), length(lat));

  16. for x = 1:length(lon1);  % 经度范围
  17.     for y = 1:length(lat1);  % 纬度范围
  18.         z0 = 1;  % 垂直层数
  19.         t0 = 1;  % 时间步长
  20.         fi = lat(y) * pi / 180;  % 纬度转弧度
  21.         beta = 2 * cos(fi) * 7.29e-5;  % beta参数
  22.         
  23.         % 计算dvo_dx和dvo_dy
  24.         dx = (lon(2) - lon(1)) * pi / 180;  % 计算经度间隔
  25.         dy = (lat(2) - lat(1)) * pi / 180;  % 计算纬度间隔
  26.         dvo_dx = double(diff(vo(x, :, z0, t0), 1, 2)) ./ double(dx) * 6371000 * cos(fi);  % x方向上的vo梯度
  27.         dvo_dy = double(diff(vo(:, y, z0, t0), 1, 1)) ./ double(dy) * 6371000;  % y方向上的vo梯度
  28.         % 补充缺失的最后一个元素
  29.         dvo_dx = [dvo_dx, dvo_dx(end)];
  30.         dvo_dy = [dvo_dy; dvo_dy(end)];

  31.         V1 = -double(u(x, y, z0, t0)) .* dvo_dx - double(v(x, y, z0, t0)) .* dvo_dy;
  32.         V4 = double(v(x, y, z0, t0)) .* beta;
  33.         V1_all(x, y) = V1;  % 将计算结果存储到数组中
  34.         V4_all(x, y) = V4;
  35.     end
  36. end
复制代码

盲写的,我也不知道能不能解决你的问题,如果有报错就更好了
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-3-27 10:05:22 | 显示全部楼层
Sigmas 发表于 2023-3-27 00:51
盲写的,我也不知道能不能解决你的问题,如果有报错就更好了

非常感谢您的回答,我按照您的修改了一下,有v1_all数组了,但依然是空值,不会附图,报错如下:
错误使用  ./
对于此运算,数组的大小不兼容。

出错 V1ANDV4 (第 25 行)
        dvo_dx = double(diff(vo(x, :, z0, t0), 1, 1)) ./ double(dx) * 6371000 * cos(fi);  % x方向上的vo梯度

密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-3-27 10:05:52 | 显示全部楼层
for i=1:length(lon1)
for j=1:length(lat1)
这里是有问题的
应该是[m,n]=size(lon1)
for j=1:m
for j=1:n
xxxxxx
end
end   
这样就能避免你说的空值得问题
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-3-27 10:09:42 | 显示全部楼层
另外 实际上matlab应该不需要这么多层for循环   式子的每一项用一行应该就能搞定 只是那个计算相邻格点距离的时候 可能用循环更好理解一点  用向量化的方法计算格点距离 可能会看起来不好理解
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-3-27 12:03:17 来自手机 | 显示全部楼层
wjy_ecnu 发表于 2023-3-27 10:09
另外 实际上matlab应该不需要这么多层for循环   式子的每一项用一行应该就能搞定 只是那个计算相邻格点距离 ...

不用for循环怎么求那么多格点上的值呢?
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-3-27 12:36:40 | 显示全部楼层
小蛇 发表于 2023-3-27 10:05
非常感谢您的回答,我按照您的修改了一下,有v1_all数组了,但依然是空值,不会附图,报错如下:
错误使 ...
  1. dvo_dx = double(diff(vo(x, :, z0, t0), 1, 2)) * 6371000 * cos(fi) ./ double(dx);  % x方向上的vo梯度
  2. dvo_dy = double(diff(vo(:, y, z0, t0), 1, 1)) * 6371000 ./ double(dy);  % y方向上的vo梯度
复制代码

改一下计算顺序试试呢?
另外以我浅薄的理解也许可以稍微回答一下你对于下面老哥“如何不用for循环”的疑惑,在这个实例中用for循环确实更直观,但改用其他方式运算也许速度可以更快并且更简洁(比如用wjy老哥说到的矩阵向量化)
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-3-27 15:04:05 | 显示全部楼层
散度 divergence  涡度 curl  梯度gradient  梯度注意matlab会自动改变x和y的顺序   但是这些函数不直观
如果研究一下 仍然没有理解 还是用for循环吧  毕竟循环逻辑清楚
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-3-27 17:27:13 | 显示全部楼层
wjy_ecnu 发表于 2023-3-27 15:04
散度 divergence  涡度 curl  梯度gradient  梯度注意matlab会自动改变x和y的顺序   但是这些函数不直观
...

好的,谢谢您的建议。
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-3-27 19:18:10 | 显示全部楼层
本帖最后由 小蛇 于 2023-3-27 19:45 编辑
Sigmas 发表于 2023-3-27 12:36
改一下计算顺序试试呢?
另外以我浅薄的理解也许可以稍微回答一下你对于下面老哥“如何不用for循环” ...

我再试试,非常感谢
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表