function [V,w,Mc,S,SS,r,rm,imatch ] = hurrprofile(chi,wc,ep,ro)
% Numerical integration of the outer wind and entropy profiles from simple
% ODEs 
% Inputs: chi, nondimensional subsidence velocity, wc, precipitation
% efficiency, ep, and nondimensional outer radius, ro
%
% Outputs:  Nondimensional values of azimuthal velocity, V, large-scale
%   vertical velocity, w, boundary layer moist static energy,
%   S, saturation moist static energy, SS, radius, r, radius of maximum
%   winds, rm, index (imatch) of matching radius.
%
% Created December, 2021
%----------------------------------------------------------------------
g=0.04;     % Gust factor (fraction of maximum wind speed)
ass=0.2;    % Asselin filter coefficient
dr=0.001;   % Nondimensional radius increment used for outer wind
alpha=1;    % C_k/C_D
%-----------------------------------------------------------------------
V=NaN;w=NaN;Mc=NaN;S=NaN;SS=NaN;r=NaN;rm=NaN;imatch=NaN;  % default values
if chi*g*alpha >= wc
    disp(['chi*g*alpha must be less than wc, stopping'])
    return
end    
rnd=0:dr:ro;
q=max(size(rnd));
%
% Integrate outer wind ODE
%
M=zeros(1,q)+0.5*ro^2;
y=zeros(1,q);
v=zeros(1,q);
w=zeros(1,q)-wc;
h=zeros(1,q);
hs=zeros(1,q);
Res=zeros(1,q);
h0=(alpha*g-wc/chi)/(alpha*g+wc);
h(q)=h0;
%h(q)=0;  % Possible alternative BC if deep convection goes to outer boundary
h(q-1)=h(q);
y(q-1)=0.5*wc*(ro^2-(ro-dr)^2);
%
p=[1,0 g^2-wc^2 2*wc^2*dr -wc^2*dr^2];  %  Checked, Dec 29 2021
rt=roots(p);
for i=1:4
    if abs(imag(rt(i))) < 1e-4 && real(rt(i)) > 0
        v(q-1)=real(rt(i));
    end
end    
%
M(q-1)=rnd(q-1)*v(q-1)+0.5*rnd(q-1)^2;
%
clear imax
ismatch=0;
sreg=0;
for i=q-2:-1:1
    if h(i+1) >= hs(i+1)
        sreg=1;
        ismatch=max(ismatch,i+1);
    end    
    vabs=sqrt(v(i+1)^2+g^2);
    if sreg == 0
        y(i)=0.5*wc*(ro^2-rnd(i)^2);
        M(i)=M(i+2)-2*dr*rnd(i+1)^2*v(i+1)*vabs/y(i+1);
        v(i)=M(i)/rnd(i)-0.5*rnd(i);
        h(i)=h(i+1)+(M(i)-M(i+1))*(wc*(h(i+1)+1/chi)/vabs-alpha*(1-h(i+1))-vabs^2)/(rnd(i+1)*v(i+1));
        h(i+1)=h(i+1)+ass.*(h(i)+h(i+2)-2*h(i+1));
        hs(i)=hs(i+2)+(M(i+2)-M(i))*v(i+1)/rnd(i+1);
        w(i)=-wc;
    else    
        w(i)=(chi*ep*alpha*vabs*(1-hs(i+1)-v(i+1)^2)/(1+chi*hs(i+1))-wc)/(1-ep);
        y(i)=y(i+2)-2*dr*rnd(i+1)*w(i);
        M(i)=M(i+2)-2*dr*rnd(i+1)^2*v(i+1)*vabs/y(i+1);
        v(i)=M(i)/rnd(i)-0.5*rnd(i);
        hs(i)=hs(i+2)+(M(i+2)-M(i))*v(i+1)/rnd(i+1);
        h(i)=hs(i);
    end    
    y(i+1)=y(i+1)+ass.*(y(i)+y(i+2)-2*y(i+1));
    M(i+1)=M(i+1)+ass.*(M(i)+M(i+2)-2*M(i+1));
    v(i+1)=M(i+1)/rnd(i+1)-0.5*rnd(i+1);
    hs(i+1)=hs(i+1)+ass*(hs(i)+hs(i+2)-2*hs(i+1));
    %
    % Calculate difference (Res) between slopes of inner and outer profiles
    %
    vi=1./(v(i)+0.5*rnd(i));
    x=vi+sqrt(vi^2-1+vi*rnd(i));
    dvdr=(2*x/rnd(i)+1)*(1-x^2)/(1+x^2)^2-0.5;
    Res(i)=dvdr-(v(i+1)-v(i))/dr;
    if Res(i+1) <= 0 && Res(i) > 0 
    %if v(i) > 1    
        imax=i;
        imatch=imax;
        rm=rnd(i)/x;
        break
    end
end
%
Mc=(w+wc)/ep;   % Convective mass flux
if min(Mc) < -0.01
    disp(['Warning: Convective mass flux is negative, change outer boundary condition']) %#ok<NBRAK> 
end    
%
if i == 1
    disp(['No matching point found, returning']) %#ok<NBRAK>
    rm=0;
    imatch=1;
    V=zeros(1,q);
    r=zeros(1,q);
    return
end    
%ep1=Res(imax)/(Res(imax)-Res(imax+1));
%rm=rm1(imax+1)*(1-ep1)+ep1*rm1(imax);
%
% Inner profile
%
rt=rnd(1:imax)/rm;
v(1:imax)=(2*rt)./(1+rt.^2)-0.5*rnd(1:imax);
for i=imax-1:-1:1
    M(i)=rnd(i)*v(i)+0.5*rnd(i)^2;
    hs(i)=hs(i+2)+(M(i+2)-M(i))*v(i+1)/rnd(i+1);
    hs(i+1)=hs(i+1)+ass*(hs(i)+hs(i+2)-2*hs(i+1));
    h(i)=hs(i);
    if i > 1
        Mm=rnd(i-1)*v(i-1)+0.5*rnd(i-1)^2;
    else
        Mm=0;
    end    
    w(i)=2*((rnd(i+1)^2*sqrt(g^2+v(i+1)^2)*v(i+1))/(M(i+2)-M(i))-(rnd(i)^2*sqrt(g^2+v(i)^2)* ...
        v(i)/(M(i+1)-Mm)))/(rnd(i)+1e-6);
    Mc(i)=(w(i)+wc)/ep; 
end
%
%w(imax:end)=100*w(imax:end);
%Mc(imax:end)=100*Mc(imax:end);
w(1)=2*w(3)-w(2);
Mc(1)=2*Mc(3)-Mc(2);
w(imax-2)=w(imax-1);
Mc(imax-2)=Mc(imax-1);
r=rnd;
V=v;
S=h;
SS=hs;
%
end