clear;clc;
filename1='D:\TS_201801_GLB.nc';
filename2='D:\TS_201807_GLB.nc';
ncdisp(filename1); % 在这里改filename1/filename1
lon=ncread(filename1,'LONGITUDE');
lat=ncread(filename1,'LATITUDE');
lat=lat(66:97); % 4.5°N~35.5°N
lon=lon(116:285); % 115.5°E~75.5°W
pre=ncread(filename1,'PRES'); % [db]
% P0=pre(23);
p=[pre(1) pre(6) pre(10) pre(13)];
index=[1 6 10 13];
TOI=ncread(filename1,'TOI');
SOI=ncread(filename1,'SOI');
T=TOI(116:285,66:97,:);
S=SOI(116:285,66:97,:);
for i=1:23
svan(:,:,i)=sw_svan(S(:,:,i),T(:,:,i),pre(i)); % Specific volume anomaly [m3/kg]
end
for i=1:22
svan_ave(:,:,i)=(svan(:,:,i)+svan(:,:,i+1))/2; % δ平均
d_svan(:,:,i)=(pre(i+1)-pre(i))*1e4.*svan_ave(:,:,i); % △δ
end
[a,b,c]=size(d_svan);
[ln,la]=meshgrid(lon,lat);
for i=1:b-1
% for j=1:a
% Lx(i,j)=distance(la(i,j),ln(i,j),la(i+1,j),ln(i,j))/180*pi*6370e3;
x(i,:)=sw_dist(la(i,:),ln(i,:),'km')*1e3; % [m]
% end
end
% for i=1:b
for j=1:a-1
% Ly(i,j)=distance(la(i,j),ln(i,j),la(i,j),ln(i,j+1))/180*pi*6370e3; % [m]
y(:,j)=sw_dist(la(:,j),ln(:,j),'km')*1e3; % [m]
end
% end
f=2*(7.292e-5)*sind(la);
% 先画个第一层试试
for k=1:4
n=index(k);
fai=zeros(a,b);
for j=n:22
fai=fai+d_svan(:,:,j);
end
fai=fai';
for i=1:b-1
for j=1:a-1
V(i,j)=(fai(i,j+1)-fai(i,j))/(f(i,j)*x(i,j));
end
end
for i=1:b-1
for j=1:a-1
U(i,j)=-(fai(i+1,j)-fai(i,j))/(f(i,j)*y(i,j));
end
end
[LON,LAT]=meshgrid([-244:-76],[5:35]);
figure(k);
m_proj('equidistant','lon',[-245 -75],'lat',[4 36]);
m_coast('patch',[.8 .8 .8]);
hold on
m_quiver(LON,LAT,U ,V,'linewi',0.9); % ,'color','b'
m_grid('xtick',8,'tickdir','in','XaxisLocation','bottom','fontsize',11);%添加格网
title({['Flow Field in North Pacific,2018.01(P=',num2str(p(k)),'db)']},'fontsize',20,'fontweight','bold'); % 改标题月份
xlabel('Lontitude(°)','fontsize',13);
ylabel('Latitude (°)','fontsize',13,'position',[-1.56 0.35]);
set(gcf,'position',[50 300 1200 270]);
end
%% 断面流场
delta=d_svan(22:23,:,:); % lon:136.5~137.5°E
[ln,la]=meshgrid([136.5 137.5]',lat);
[d,e,g]=size(delta);
f=2*(7.292e-5)*sind(la);
for j=1:d-1
y(:,j)=sw_dist(la(:,j),ln(:,j),'km')*1e3; % [m]
end
clear U;
for k=1:22
fai=zeros(d,e);
for j=k:22
fai=fai+squeeze(delta(:,:,j));
end
fai=fai';
for i=1:e-1
for j=1:d-1
U(i,j)=-(fai(i+1,j)-fai(i,j))/(f(i,j)*y(i,j));
end
end
UU(:,k)=squeeze(U);
end
UU(:,23)=zeros(31,1);
[cs,h]=contourf(lat(1:31),pre(1:23),UU');
set(gca,'YDir','reverse','ytick',[0:100:1500],'xlim',[5 32],'ylim',[10 1500]);
set(h,'ShowText','on');
caxis([-0.4,0.4]);
colorbar;
xlabel('Latitude(°)','fontsize',16);
ylabel('p(db)','fontsize',16);
title({'Vertical Section in North Pacific 137°E, 2018.01'},'fontsize',25,'fontweight','bold'); % 改标题月份
flow=0;
for i=3:21 % lat:7°N~25°N
for j=1:18 % p:10db~1000db
flow=flow-UU(i,j)*(pre(j+1)-pre(j))*y(i);
end
end
flow./1e6 %[Sv]
评论6