Matlab 画全球风场图并计算散度风和旋度风

clear all;close all;clc
%%
fname = 'wind_uv300.nc';
lat   = ncread(fname, 'lat');
lon   = ncread(fname, 'lon');
uwind = ncread(fname, 'U');
vwind = ncread(fname, 'V');
uwind = uwind(:,:,2);
vwind = vwind(:,:,2);
%
[lat2d, lon2d] = meshgrid(lat,lon);

% plot
figure(1)
set( gcf , 'color' ,'w');
m_proj('miller','lat',[-85, 85],'lon',[-180, 180]);
m_coast('patch',[.7 .7 .7],'edgecolor','none');
m_grid('box','on','linestyle','none');

hold on
% draw data
s0 = 4;
m_quiver( lon2d(1:4:end,1:4:end) , lat2d(1:4:end,1:4:end), ...
          uwind(1:4:end,1:4:end)/s0, vwind(1:4:end,1:4:end)/s0,0,'color','k');%0是拉伸程度
% referenec vector
m_quiver(-60,20, 120/s0,0,0,'color','b');%东经110度,北纬30度的水平箭头,30m/s作为参考,z=0

title( 'Global wind' , 'fontsize' , 18 , 'fontweight' , 'bold' );
%%   计算散度风
[lat2d, lon2d] = meshgrid(double(lat),double(lon));

% grid distance
dd0 = 2*pi*6371000/360;  % distance for one Latitude (unit: meter)
xx = lon2d*dd0.*cosd( lat2d );
yy = lat2d*dd0;
% calculate divergence
div = divergence( xx', yy', uwind', vwind');%要求第二维变化

figure(2)
set( gcf , 'color' ,'w');
m_proj('miller','lat',[-87, 87],'lon',[-180, 180]);
m_contourf(lon2d,lat2d,div',20);
m_coast('patch',[.7 .7 .7]);
m_grid('box','on','linestyle','none');
colormap('jet');colorbar
caxis([-1*10^-4 1*10^-4]);
title('散度风分布图')
%%   计算旋度风
[lat2d, lon2d] = meshgrid(double(lat),double(lon));

% grid distance
dd0 = 2*pi*6371000/360;  % distance for one Latitude (unit: meter)
xx = lon2d*dd0.*cosd( lat2d );
yy = lat2d*dd0;
% calculate rotation
vor = curl(xx', yy', uwind', vwind');

figure(3)
set( gcf , 'color' ,'w');
m_proj('miller','lat',[-87, 87],'lon',[-180, 180]);
m_contourf(lon2d,lat2d,vor',20);
m_coast('patch',[.7 .7 .7]);
m_grid('box','on','linestyle','none');
colormap('jet');colorbar
caxis([-1*10^-4 1*10^-4]);
title('旋度风分布图')