Rossby radius calculation

gistfile1.matlab

%-----calculate the pressure at the tropopause
% do lat=ncread, pfull=ncread first!
pfull_range=20; % levels of pressure in modelling
lat_point=49; % 
T_p=zonal_temp(:,lat_point);
d_T=T_p(1:pfull_range-2)-T_p(3:pfull_range);
d_pfull=pfull(1:pfull_range-2)-pfull(3:pfull_range);
dT_dz=-9.8/287*pfull(2:pfull_range-1).*d_T./d_pfull./T_p(2:pfull_range-1);
subplot(1,2,1)
plot(pfull(2:pfull_range-1),1000*dT_dz,'linewidth',2)
xlabel('pressure','fontsize',16)
ylabel('lapse rate (K/km)','fontsize',16)
set(gca,'fontsize',16)

n_lower=2;
n_reverse=find(dT_dz==min(dT_dz));
P_trop=interp1(1000*dT_dz(n_lower:n_reverse),pfull(1+n_lower:(1+n_reverse)),-2,'pchip');
hold on
plot(interp1(1000*dT_dz(n_lower:n_reverse),pfull(1+n_lower:(1+n_reverse)),...
    linspace(1000*dT_dz(1),min(1000*dT_dz),10),'pchip'),linspace(1000*dT_dz(1),min(1000*dT_dz),10),...
    '--r','linewidth',2);

set(gca,'fontsize',16)
title(strcat('topopause=',...
    num2str(P_trop,3),', lat=',num2str(lat(lat_point),2)),...
    'fontsize',16)

%---potential temperature profile, calculate the Rossby radius
poten_temp=T_p.*(1000./pfull).^(2/7);
subplot(1,2,2)
plot(pfull,poten_temp,'linewidth',2)
set(gca,'fontsize',16)
xlabel('pressure','fontsize',16)
ylabel('potential temperature','fontsize',16)


y_p_temp=poten_temp(round(pfull_range/2):pfull_range-2);
x_pfull=pfull(round(pfull_range/2):pfull_range-2);

% linear fit to get d_theta/d_p
[xData, yData] = prepareCurveData( x_pfull, y_p_temp );
ft = fittype( 'poly1' );
opts = fitoptions( ft );
opts.Lower = [-Inf -Inf];
opts.Upper = [Inf Inf];

[fitresult, gof] = fit( xData, yData, ft, opts );

Ps=700; % unit mb
f=2*pi/(24*3600)*sin(pi/180*lat(lat_point)); % changes with latitude
Np_2=-287/Ps*(Ps/1000)^(2/7)*fitresult.p1;
LR=sqrt(Np_2)*(Ps-P_trop)/f/1000 % unit km
title(strcat('L_R=',num2str(LR,4),', Np=',num2str(sqrt(Np_2),3)),...
    'fontsize',16)

%Rossby_wavenumber=2*pi*sin(pi/4)*6400/LR

%% calculate the Rossby radius for a range of latitude
exp_num=['1','2','3','4'];
filename_prefix='Nov28_exp';
pfull_range=20; % levels of pressure in modelling
%lat_range=43:1:56; % from 29N to 66N
lat_range_set=[34,60;
    34,60;
    34,60;
    34,60]; % use different range for different experiments
color_scheme={'r','b','g','k','m','c','--r','--b','--g','--k','--m','--c',':r',':b',':g',':k',':m',':c'...
    '--sr','--sb','--sg','--sk','--sm','--sc'};

for loop=1:length(exp_num)
    load(strcat(filename_prefix,num2str(exp_num(loop)),'_zonal_temp.mat'),'zonal_temp');
    lat_range=lat_range_set(loop,1):1:lat_range_set(loop,2);
    Rossby_radius_set=zeros(length(lat_range),1);
    
    for i=1:length(lat_range)
        lat_point=lat_range(i);
        T_p=zonal_temp(:,lat_point);
        d_T=T_p(1:pfull_range-2)-T_p(3:pfull_range);
        d_pfull=pfull(1:pfull_range-2)-pfull(3:pfull_range);
        dT_dz=-9.8/287*pfull(2:pfull_range-1).*d_T./d_pfull./T_p(2:pfull_range-1);

        n_lower=1;
        n_reverse=find(dT_dz==min(dT_dz));
        P_trop=interp1(1000*dT_dz(n_lower:n_reverse),pfull(1+n_lower:(1+n_reverse)),-2,'pchip');


        %---potential temperature profile, calculate the Rossby radius
        poten_temp=T_p.*(1000./pfull).^(2/7);
        y_p_temp=poten_temp(round(pfull_range/2):pfull_range-2);
        x_pfull=pfull(round(pfull_range/2):pfull_range-2);

        % linear fit to get d_theta/d_p
        [xData, yData] = prepareCurveData( x_pfull, y_p_temp );
        ft = fittype( 'poly1' );
        opts = fitoptions( ft );
        opts.Lower = [-Inf -Inf];
        opts.Upper = [Inf Inf];

        [fitresult, gof] = fit( xData, yData, ft, opts );

        Ps=700; % unit mb
        f=2*pi/(24*3600)*sin(pi/180*lat(lat_point));
        Np_2=-287/Ps*(Ps/1000)^(2/7)*fitresult.p1;
        LR=sqrt(Np_2)*(Ps-P_trop)/f/1000; % unit km
        Rossby_radius_set(i)=LR;
    end

    plot(lat(lat_range),Rossby_radius_set,char(color_scheme(loop)),'linewidth',2)
    hold on
end

xlabel('latitude','fontsize',16)
ylabel('Rossby radius (km)','fontsize',16)
set(gca,'fontsize',16)