%  TC PBL Model
%  Integrates nondimensional PBL equations in time until steady state 
%  achieved, given a gradient wind profile. Plots all three velocity
%  components at discrete intervals in time. 
%  K. Emanuel, July 2012, last modified April, 2013
%----------------------------------------
vm=70;      % Maximum gradient wind speed (m/s)
rm=50;      % Radius of maximum winds (km)
f=5e-5;     % Coriolis parameter (1/s)
h=1000;     % PBL depth (m)
lh=1500;    % Horizontal mixing length (m)
CD=1.5e-3;  % Surface drag coefficient
ckcd=2;     % Ratio of C_k to C_D (used in the specification of the gradient wind profile)
%-----------------------------------------
%
%  The following are calculated from the dimensional parameters specified
%  above. Or, their values can be specified here directly. 
%
pnu=f^3*lh^2*h/(CD*vm^3);    % pnu parameter
chi=(f*h/(CD*vm))^2;         % chi parameter
rm=1000*rm*f/vm;             % non-dimensional radius of maximum winds
%----------------------------------------
% Numerical and other parameters
%
dr=0.004;       % radial grid spacing (nondimensional)
dt=0.001;       % time step (nondimensional)
endtime=2;      % nondimensional end time
plot_dt=0.2;    % nondimensional time interval used in plots
eyesolid='n';   % Impose solid body gradient wind inside rm?
af=0.2;         % Asselin filter value for time marching
%--------------------------------------------------------------------------
r=0:dr:1;
t=0:dt:endtime;
%
qr=max(size(r));
qt=max(size(t));
dtp=floor(qt.*plot_dt/endtime);
dmgdr=zeros(1,qr);
vg=zeros(1,qr);
v=zeros(qt,qr);
u=zeros(qt,qr);
m=zeros(qt,qr);
wdisp=zeros(qt,qr);
dri=1./dr;
%
%  Specify gradient angular momentum profile from Emanuel & Rotunno, 2011
%
rn=r./rm;
if ckcd == 2  % (Special case; limit of full formula as ckcd goes to 2)
    mg=rm.*exp(-(1-rn.^2)./(2.*rn.^2));
elseif ckcd < 2    
    mg=rm.*(2.*rn.^2./(2-ckcd+ckcd.*rn.^2)).^(1./(2-ckcd));
else
    mg=rm.*max(0.5.*(2-ckcd+ckcd.*rn.^2)./rn.^2,0).^(1./(ckcd-2));
end    
%
% Set to solid body rotation inside rm if desired
%
if strcmp(eyesolid,'y');
    for i=1:qr,
        if r(i) <= rm
            mg(i)=r(i).^2./rm;
        end
    end    
end
%
mg=max(mg,0.5.*r.^2);  %  Don't let gradient wind be negative
%
m(1,1:qr)=mg;          % Initial angular momentum
vg(2:end)=mg(2:end)./r(2:end)-0.5.*r(2:end);  % Calculate gradient wind from m
v(1,1:qr)=vg(1:qr);   %  Set initial PBL v to vg
dmgdr(2:qr-1)=(mg(3:qr)-mg(1:qr-2))./(2.*dr); % Radial gradient of m
u(1,2:qr-1)=-r(2:qr-1).*v(1,2:qr-1).^2./dmgdr(2:qr-1); % Initial estimate of u
ug=u;
%
% Time march
%
%  First time step forward
%
Sp=(dri.*(v(1,3:qr)-v(1,2:qr-1))-(v(1,3:qr)+v(1,2:qr-1))./(r(3:qr)+r(2:qr-1))).^2 + ...
    (2./chi).*((dri.*(u(1,3:qr)-u(1,2:qr-1))).^2+((u(1,3:qr)+u(1,2:qr-1))./ ...
    (r(3:qr)+r(2:qr-1))).^2);
Sm=(dri.*(v(1,2:qr-1)-v(1,1:qr-2))-(v(1,2:qr-1)+v(1,1:qr-2))./(r(2:qr-1)+r(1:qr-2))).^2 + ...
    (2./chi).*((dri.*(u(1,2:qr-1)-u(1,1:qr-2))).^2+((u(1,2:qr-1)+u(1,1:qr-2))./ ...
    (r(2:qr-1)+r(1:qr-2))).^2);
Sp=sqrt(Sp);
Sm=sqrt(Sm);
S=0.5.*(Sm+Sp);
we=-0.5.*dri.*(r(3:qr).*u(1,3:qr)-r(1:qr-2).*u(1,1:qr-2))./r(2:qr-1);
we=min(we,0);
u(2,2:qr-1)=u(1,2:qr-1)+dt.*(chi.*(m(1,2:qr-1).^2-mg(2:qr-1).^2)./r(2:qr-1).^3- ...
    (v(1,2:qr-1).^2+u(1,2:qr-1).^2./chi).^0.5.*u(1,2:qr-1) - ...
    (u(1,3:qr).^2-u(1,1:qr-2).^2).*dri + ...
    2.*pnu.*(dri.^2.*(Sp.*(r(2:qr-1)+0.5.*dr).*(u(1,3:qr)-u(1,2:qr-1))-Sm.* ...
    (r(2:qr-1)-0.5.*dr).*(u(1,2:qr-1)-u(1,1:qr-2)))-S.*u(1,2:qr-1))./r(2:qr-1)+ ...
    we.*u(1,2:qr-1));
m(2,2:qr-1)=m(1,2:qr-1)+dt.*(-r(2:qr-1).*v(1,2:qr-1).*(v(1,2:qr-1).^2+ ...
    u(1,2:qr-1).^2./chi).^0.5-u(1,2:qr-1).*0.5.*dri.*(m(1,3:qr)-m(1,1:qr-2)) + ...
    pnu.*dri.^2.*(Sp.*(r(2:qr-1)+0.5.*dr).^3.*(v(1,3:qr)./r(3:qr)-v(1,2:qr-1)./r(2:qr-1)) - ...
    Sm.*(r(2:qr-1)-0.5.*dr).^3.*(v(1,2:qr-1)./r(2:qr-1)-v(1,1:qr-2)./max(r(1:qr-2),0.4*dr)))./r(2:qr-1)- ...
    we.*(mg(2:qr-1)-m(1,2:qr-1))); 
v(2,2:qr-1)=m(2,2:qr-1)./r(2:qr-1)-0.5.*r(2:qr-1);
%
disp('Calculating... Wait for Menu to Appear...')
%
%  March forward
%
for i=2:qt-1,
    Sp=(dri.*(v(i-1,3:qr)-v(i-1,2:qr-1))-(v(i-1,3:qr)+v(i-1,2:qr-1))./(r(3:qr)+r(2:qr-1))).^2 + ...
        (2./chi).*((dri.*(u(i-1,3:qr)-u(i-1,2:qr-1))).^2+((u(i-1,3:qr)+u(i-1,2:qr-1))./ ...
        (r(3:qr)+r(2:qr-1))).^2);
    Sm=(dri.*(v(i-1,2:qr-1)-v(i-1,1:qr-2))-(v(i-1,2:qr-1)+v(i-1,1:qr-2))./(r(2:qr-1)+r(1:qr-2))).^2 + ...
        (2./chi).*((dri.*(u(i-1,2:qr-1)-u(i-1,1:qr-2))).^2+((u(i-1,2:qr-1)+u(i-1,1:qr-2))./ ...
        (r(2:qr-1)+r(1:qr-2))).^2);
    Sp=sqrt(Sp);
    Sm=sqrt(Sm);
    S=0.5.*(Sm+Sp);
    %
    w=-0.5.*dri.*(r(3:qr).*u(i,3:qr)-r(1:qr-2).*u(i,1:qr-2))./r(2:qr-1);
    wdisp(i,2:qr-1)=w;
    we=min(w,0);
    %
    %uadv=(u(i,3:qr).^2-u(i,1:qr-2).^2).*dri;  % Centered differencing
    %
    % Upstream differencing of advection term
    %
    uadv=0.5.*dri.*(min(u(i,3:qr)+u(i,2:qr-1),0).*(u(i,3:qr)-u(i,2:qr-1))+ ...
        max(u(i,2:qr-1)+u(i,1:qr-2),0).*(u(i,2:qr-1)-u(i,1:qr-2)));
    %
    u(i+1,2:qr-1)=u(i-1,2:qr-1)+2.*dt.*(chi.*(m(i,2:qr-1).^2-mg(2:qr-1).^2)./r(2:qr-1).^3- ...
        (v(i-1,2:qr-1).^2+u(i-1,2:qr-1).^2./chi).^0.5.*u(i-1,2:qr-1) - uadv + ...
        2.*pnu.*(dri.^2.*(Sp.*(r(2:qr-1)+0.5.*dr).*(u(i-1,3:qr)-u(i-1,2:qr-1))-Sm.* ...
        (r(2:qr-1)-0.5.*dr).*(u(i-1,2:qr-1)-u(i-1,1:qr-2)))-S.*u(i-1,2:qr-1))./r(2:qr-1)+ ...
        we.*u(i-1,2:qr-1));
    %
    madv1=u(i,2:qr-1).*0.5.*dri.*(m(i,3:qr)-m(i,1:qr-2)); % Centered differencing
    %
    % Upstream differencing of advection term
    %
    madv2=0.5.*dri.*(min(u(i,3:qr)+u(i,2:qr-1),0).*(m(i,3:qr)-m(i,2:qr-1))+ ...
        max(u(i,2:qr-1)+u(i,1:qr-2),0).*(m(i,2:qr-1)-m(i,1:qr-2)));
    %
    madv=0.5.*(madv1+madv2);  % Average of differencing methods seems to work best here
    %
    m(i+1,2:qr-1)=m(i-1,2:qr-1)+2.*dt.*(-r(2:qr-1).*v(i-1,2:qr-1).*(v(i-1,2:qr-1).^2+ ...
        u(i-1,2:qr-1).^2./chi).^0.5 - madv + ...
        pnu.*dri.^2.*(Sp.*(r(2:qr-1)+0.5.*dr).^3.*(v(i-1,3:qr)./r(3:qr)-v(i-1,2:qr-1)./r(2:qr-1)) - ...
        Sm.*(r(2:qr-1)-0.5.*dr).^3.*(v(i-1,2:qr-1)./r(2:qr-1)-v(i-1,1:qr-2)./max(r(1:qr-2),0.4*dr)))./r(2:qr-1)- ...
        we.*(mg(2:qr-1)-m(i-1,2:qr-1)));
    m(i+1,2:qr-1)=max(m(i+1,2:qr-1),0.5.*r(2:qr-1).^2);
    v(i+1,2:qr-1)=m(i+1,2:qr-1)./r(2:qr-1)-0.5.*r(2:qr-1);
    %
    %  Apply Asselin filter
    %
    u(i,2:qr-1)=u(i,2:qr-1)+af.*(u(i+1,2:qr-1)+u(i-1,2:qr-1)-2.*u(i,2:qr-1));
    m(i,2:qr-1)=m(i,2:qr-1)+af.*(m(i+1,2:qr-1)+m(i-1,2:qr-1)-2.*m(i,2:qr-1));
    v(i,2:qr-1)=m(i,2:qr-1)./r(2:qr-1)-0.5.*r(2:qr-1);
    %v(i,2:qr-1)=max(v(i,2:qr-1),0);
end    
%
% Display plots
%
t1=floor(0.8.*qt);  % Times to display variables at
t2=floor(0.9.*qt);
t3=qt;
%
[~,imax]=max(v(end,:));
[~,izero]=min(max(v(end,imax:end),0));
%
choice = 1;
while choice ~= 4
choice=menu('Plot Menu','Azimuthal Velocity','Radial Velocity','Vertical Velocity','Exit');
%
if choice == 1
    h=plot(r,v(1:dtp:end,:));
    set(h,'linewidth',1.5)
    hold on
      h1=plot(r,vg,'k',r,v(end,:),'b');
      set(h1,'linewidth',2);
      legend(h1,'Gradient wind','Azimuthal Wind at Final Time','location','best')
    hold off
    set(gca,'xlim',[0 r(izero)])
    set(gca,'fontweight','bold')
    xlabel('Nondimensional Radius')
    ylabel('Nondimensional v')
    title('Non-dimensional Boundary Layer Azimuthal Velocity')
    %legend('80% of End Time','90% of End Time','End Time','Gradient Wind','location','best');
    %
elseif choice == 2    
    h=plot(r,u(1:dtp:end,:));
    set(h,'linewidth',1.5)
    hold on
      h1=plot(r,u(end,:),'k');
      set(h1,'linewidth',2);
      legend(h1,'Radial Wind at Final Time','location','best')
    hold off
    set(gca,'xlim',[0 r(izero)])
    set(gca,'fontweight','bold')
    xlabel('Nondimensional Radius')
    ylabel('Nondimensional u')
    title('Non-dimensional Boundary Layer Radial Velocity')
    %legend('80% of End Time','90% of End Time','End Time','location','best');
elseif choice == 3    
    h=plot(r,wdisp(dtp:dtp:end-dtp,:));
    set(h,'linewidth',1.5)
    hold on
      h1=plot(r,wdisp(end-dtp,:),'k');
      set(h1,'linewidth',2);
      legend(h1,'Vertical Velocity at Final Time','location','best')
    hold off    
    set(gca,'xlim',[0 r(izero)])
    set(gca,'fontweight','bold')
    xlabel('Nondimensional Radius')
    ylabel('Nondimensional w')
    title('Non-dimensional Boundary Layer Vertical Velocity')    
end    
end 