这里只给部分主程序和调用的两个子程序

%% NCEP Wind Data Processing

clear,clc

uwndFileName = 'uwnd.mon.mean.nc';

vwndFileName = 'vwnd.mon.mean.nc';

lon = double(ncread(vwndFileName,'lon'));

lat = double(ncread(vwndFileName,'lat'));

level = ncread(vwndFileName,'level');

time = ncread(vwndFileName,'time');

vwnd = ncread(vwndFileName,'vwnd');

uwnd = ncread(uwndFileName,'uwnd');

vwnd = vwnd(:,:,:,1:length(time)-11);% 1948.01-2018.12

uwnd = uwnd(:,:,:,1:length(time)-11);% 1948.01-2018.12

time = (time/24 + datenum(1800,1,1,0,0,0));

time = time(1:length(time)-11);% 1948.01-2018.12

time_vec = datevec(time);

[Lo,La,z,t] = size(vwnd);

vwnd_clim = squeeze(nanmean(reshape(vwnd,[Lo,La,z,12,t/12]),5));% Climatology

uwnd_clim = squeeze(nanmean(reshape(uwnd,[Lo,La,z,12,t/12]),5));% Climatology

% monthly wind vector for year n and month m.

sigmamn = NaN(Lo,La,t,z);

refLevel = 850;% 850 hpa

v1 = squeeze(vwnd_clim(:,:,level == refLevel,1)); % Jan v

v7 = squeeze(vwnd_clim(:,:,level == refLevel,7)); % Jul v

u1 = squeeze(uwnd_clim(:,:,level == refLevel,1)); % Jan u

u7 = squeeze(uwnd_clim(:,:,level == refLevel,7)); % Jul u

sigma = Cal_sigma(lat,u1,v1,u7,v7);

for j = 1:z

refLevel = level(j);% 850;% 850 hpa

v1 = squeeze(vwnd_clim(:,:,level == refLevel,1)); % Jan v

v7 = squeeze(vwnd_clim(:,:,level == refLevel,7)); % Jul v

u1 = squeeze(uwnd_clim(:,:,level == refLevel,1)); % Jan u

u7 = squeeze(uwnd_clim(:,:,level == refLevel,7)); % Jul u

for i = 1:t

    v_temp = squeeze(vwnd(:,:,level == refLevel,i)); % mn v

    u_temp = squeeze(uwnd(:,:,level == refLevel,i)); % mn u

    [~,sigmamn(:,:,i,j)] = Cal_sigma(lat,u1,v1,u7,v7,u_temp,v_temp);

end

disp(['Level = ',num2str(level(j)),'hpa finished ...'])

end

% sigma = sigma* + 2

sigma = sigma -2; sigmamn = sigmamn - 2;

%% Index

year = (1948:1:2018)';

% South Asia Summer Moonsoon Index (SASMI) Jun-Sep 5-22.5N 35-97.5E 850hPa

SASMI = Cal_index(sigmamn,[35 97.5],[5 22.5],850,lon,lat,level,6,9);

% East Asia Summer Moonsson Index (EASMI) Jun-Aug 10-40N 110-140E 850hPa

EASMI = Cal_index(sigmamn,[110 140],[10 40],850,lon,lat,level,6,8);

function [sigma,sigma_mn] =  Cal_sigma(lat,u1,v1,u7,v7,varargin)

% Input: 

%    lat: wind field latitude grid    size: lat*1

%    u1 v1 u7 v7 January and July climatological wind vector  size: lon*lat

% 

%    varargin(1) = u_mn varargin(2) = v_mn   size: lon*lat

%    monthly wind vector for year n and month m.

% Output: 

%    Dynamical Normalized Seasonality (DNS)

% written by 2020.01.03 xiaoxin

% ref: Li and Zeng 2002 Geophysical Research Letters

% eg. sigma = Cal_sigma(lat,u1,v1,u7,v7);

% eg. sigma in month j level i    sigmamn size:lon*lat*zLevel*t

% [~,sigmamn(:,:,i,j)] = Cal_sigma(lat,u1,v1,u7,v7,u_temp,v_temp);

% South Asia Summer Moonsoon Index (SASMI) Jun-Sep 5-22.5N 35-97.5E 850hPa

% SASMI = Cal_index(sigmamn,[35 97.5],[5 22.5],850,lon,lat,level,6,9);

v_minus = v1 - v7; u_minus = u1 - u7;

if length(varargin) == 2

    v_temp = varargin{2}; u_temp = varargin{1};

    vmn_minus = v1 - v_temp; umn_minus = u1 - u_temp;

end

v_mean = (v1 + v7)./2; u_mean = (u1 + u7)./2;

% top: ||V1_bar - V7_bar||; Caution radian in angle

[Lo,La] = size(u1);

top = NaN(Lo,La); top_mn = NaN(Lo,La); bottom = NaN(Lo,La);

for i = 1:Lo

    for j = 2:La-1

        if i == 1 % i-1 = Lo

        top(i,j) = (v_minus(Lo,j).^2+u_minus(Lo,j).^2 + ...

            4*(v_minus(i,j).^2+u_minus(i,j).^2) + ...

            v_minus(i+1,j).^2+u_minus(i+1,j).^2).*cosd(lat(j)) + ...

            (v_minus(i,j-1).^2+u_minus(i,j-1).^2).*cosd(lat(j-1)) + ...

            (v_minus(i,j+1).^2+u_minus(i,j+1).^2).*cosd(lat(j+1));

        top(i,j) = sqrt(top(i,j));

        bottom(i,j) = (v_mean(Lo,j).^2+u_mean(Lo,j).^2 + ...

            4*(v_mean(i,j).^2+u_mean(i,j).^2) + ...

            v_mean(i+1,j).^2+u_mean(i+1,j).^2).*cosd(lat(j)) + ...

            (v_mean(i,j-1).^2+u_mean(i,j-1).^2).*cosd(lat(j-1)) + ...

            (v_mean(i,j+1).^2+u_mean(i,j+1).^2).*cosd(lat(j+1));

        bottom(i,j) = sqrt(bottom(i,j));

        elseif i == Lo

        top(i,j) = (v_minus(i-1,j).^2+u_minus(i-1,j).^2 + ...

            4*(v_minus(i,j).^2+u_minus(i,j).^2) + ...

            v_minus(1,j).^2+u_minus(1,j).^2).*cosd(lat(j)) + ...

            (v_minus(i,j-1).^2+u_minus(i,j-1).^2).*cosd(lat(j-1)) + ...

            (v_minus(i,j+1).^2+u_minus(i,j+1).^2).*cosd(lat(j+1));

        top(i,j) = sqrt(top(i,j));

        bottom(i,j) = (v_mean(i-1,j).^2+u_mean(i-1,j).^2 + ...

            4*(v_mean(i,j).^2+u_mean(i,j).^2) + ...

            v_mean(1,j).^2+u_mean(1,j).^2).*cosd(lat(j)) + ...

            (v_mean(i,j-1).^2+u_mean(i,j-1).^2).*cosd(lat(j-1)) + ...

            (v_mean(i,j+1).^2+u_mean(i,j+1).^2).*cosd(lat(j+1));

        bottom(i,j) = sqrt(bottom(i,j));

        else

        top(i,j) = (v_minus(i-1,j).^2+u_minus(i-1,j).^2 + ...

            4*(v_minus(i,j).^2+u_minus(i,j).^2) + ...

            v_minus(i+1,j).^2+u_minus(i+1,j).^2).*cosd(lat(j)) + ...

            (v_minus(i,j-1).^2+u_minus(i,j-1).^2).*cosd(lat(j-1)) + ...

            (v_minus(i,j+1).^2+u_minus(i,j+1).^2).*cosd(lat(j+1));

        top(i,j) = sqrt(top(i,j));

        bottom(i,j) = (v_mean(i-1,j).^2+u_mean(i-1,j).^2 + ...

            4*(v_mean(i,j).^2+u_mean(i,j).^2) + ...

            v_mean(i+1,j).^2+u_mean(i+1,j).^2).*cosd(lat(j)) + ...

            (v_mean(i,j-1).^2+u_mean(i,j-1).^2).*cosd(lat(j-1)) + ...

            (v_mean(i,j+1).^2+u_mean(i,j+1).^2).*cosd(lat(j+1));

        bottom(i,j) = sqrt(bottom(i,j));

        end

        if length(varargin) == 2

            if i == 1

            top_mn(i,j) = (vmn_minus(Lo,j).^2+umn_minus(Lo,j).^2 + ...

                4*(vmn_minus(i,j).^2+umn_minus(i,j).^2) + ...

                vmn_minus(i+1,j).^2+umn_minus(i+1,j).^2).*cosd(lat(j)) + ...

                (vmn_minus(i,j-1).^2+umn_minus(i,j-1).^2).*cosd(lat(j-1)) + ...

                (vmn_minus(i,j+1).^2+umn_minus(i,j+1).^2).*cosd(lat(j+1));

            top_mn(i,j) = sqrt(top_mn(i,j));

            elseif i == Lo

            top_mn(i,j) = (vmn_minus(i-1,j).^2+umn_minus(i-1,j).^2 + ...

                4*(vmn_minus(i,j).^2+umn_minus(i,j).^2) + ...

                vmn_minus(1,j).^2+umn_minus(1,j).^2).*cosd(lat(j)) + ...

                (vmn_minus(i,j-1).^2+umn_minus(i,j-1).^2).*cosd(lat(j-1)) + ...

                (vmn_minus(i,j+1).^2+umn_minus(i,j+1).^2).*cosd(lat(j+1));

            top_mn(i,j) = sqrt(top_mn(i,j));

            else

            top_mn(i,j) = (vmn_minus(i-1,j).^2+umn_minus(i-1,j).^2 + ...

                4*(vmn_minus(i,j).^2+umn_minus(i,j).^2) + ...

                vmn_minus(i+1,j).^2+umn_minus(i+1,j).^2).*cosd(lat(j)) + ...

                (vmn_minus(i,j-1).^2+umn_minus(i,j-1).^2).*cosd(lat(j-1)) + ...

                (vmn_minus(i,j+1).^2+umn_minus(i,j+1).^2).*cosd(lat(j+1));

            top_mn(i,j) = sqrt(top_mn(i,j));

            end

        end

    end

end

sigma = top./bottom;

if length(varargin) == 2

    sigma_mn = top_mn./bottom;

end

end

function Index = Cal_index(sigmamn,lonlim,latlim,refLevel,lon,lat,level,m1,m2)

% Calculate normalized Index at different region

% Input

%    sigmamn size:lon*lat*zLevel*t

%    area: lon1-lon2 lat1-lat2 reflevel=zlevel i.e. 850hpa

%    lon lat: Longitude and latitude matrix n*1 m*1

%    level: z*1 [1000;925;850;700;600;500;400;300;250;200;150;100;70;50;30;20;10]

%    m1 m2: monthlim

% Output

%    Index: Normalized index t*1

% written by 2020.01.03 xiaoxin

% ref: Li and Zeng 2002 Geophysical Research Letters

% eg. sigma in month j level i    sigmamn size:lon*lat*zLevel*t

% [~,sigmamn(:,:,i,j)] = Cal_sigma(lat,u1,v1,u7,v7,u_temp,v_temp);

% South Asia Summer Moonsoon Index (SASMI) Jun-Sep 5-22.5N 35-97.5E 850hPa

% SASMI = Cal_index(sigmamn,[35 97.5],[5 22.5],850,lon,lat,level,6,9);

lon1 = lonlim(1); lon2 = lonlim(2);

lat1 = latlim(1); lat2 = latlim(2);

if lon1<=lon2

    Index = sigmamn(lon>=lon1 & lon<=lon2,lat>=lat1 & lat<=lat2,:,level == refLevel);

    Index = reshape(squeeze(nanmean(nanmean(Index,1),2)),[12 length(Index)/12])';

    Index = nanmean(Index(:,m1:m2),2); % summer:JJA mean

%     Index = nanmean(Index(:,1:12),2); % All-year mean

    Index = zscore(Index);% Normalized

else

    Index =  cat(1,sigmamn(lon >= lon1,lat>=lat1 & lat<=lat2,:,level == refLevel),...

     sigmamn(lon <= lon2,lat>=lat1 & lat<=lat2,:,level == refLevel));

    Index = reshape(squeeze(nanmean(nanmean(Index,1),2)),[12 length(Index)/12])';

    Index = nanmean(Index(:,m1:m2),2); % summer:JJA mean

%     Index = nanmean(Index(:,1:12),2); % All-year mean

    Index = zscore(Index);% Normalized

end