%
%  This script solves the non-dimensional equations for the evolution of
%  the boundary layer into an initially dry, dry adiabatic atmosphere under
%  radiative surface heating. K. Emanuel, July, 2016
%  This version allows time dependence of surface radiation. December, 2021
%  Last modified November 2022 to subtract constant raidiative cooling to
%  give zero diurnal mean
%-------------------------------------------------------------------------
%
T0=305;    % Surface temperature of initial dry-adiabatic atmosphere
delT=15;   % Initial depression of surface air temperature (K) in new PBL
RHinit=70; % Initial relative humidity (%) of surface air at new PBL top
Q=700;     % Peak (noon) surface radiative flux (W/m^2)
h0=300;    % Initial depth (m) of the PBL 
hD=4000;   % Height (m) above which dry static energy increases to tropopause 
           % Note: Level of free convection assumed to be below this level
Tanvil=200;% Tropopause temperature
vs=10;      % Surface wind speed (m/s)
alpha=0.8; % Coefficient multiplying surface evaporation (soil moisture parameter)
Ck=1.5e-3; % Surface sensible and latent heat exchange coefficient
%
A=0.2;     % Lilly coefficient governing convective PBL top entrainment
dt=0.002;  % Nondimensional time step
endtime=1; % Nondimensional end time (integration stops before then, though, when D=1)
%-------------------------------------------------------------------------
rho=1.1;   % Surface density (kg/m^3)
cp=1005;   % Heat capacity
R=287;     % Gas constant for dry air
Rv=461;    % Gas constant for water vapor
Lv=2.5e6;  % Latent heat of vaporization
g=9.8;     % Gravity
p0=1005;   % Surface pressure (hPa)
afilt=0.1; % Asselin coefficient
%-------------------------------------------------------------------------
%
D0=cp*T0;  % Dry static energy of initial dry adiabatic environment
%
T=12*3600; % Length of day (12 hours)
%
% Nondimensional numbers
%
chi=Ck*vs*T/h0;
Phi0=Q*T/(rho*D0*h0);
%
% Calculate initial value of D assuming wet bulb temperature at top of PBL
%
Th=T0-g*h0/cp;   %  Initial temperature at altitude of PBL top
%
p=p0*(cp*Th/D0).^(cp/287); % Pressure at PBL top
%
Q1=Q/(rho*Ck*vs);  % Normalized heating used to find surface temperature
%
D=cp*(T0-delT);    % Initial dimensional dry static energy of new PBL 
Tc=T0-delT-g*h0/cp-273.15;
es=6.112*exp(17.67*Tc/(243.5+Tc));
qinit=0.01*RHinit*0.622*es/(p-es);
%
M=(D+Lv*qinit-D0)/D0;    % Initial normalized moist static energy of new PBL
%
%   Normalize initial dry static energy
%
D=1-D/D0;
%
D1=D;  % Initial values of normalized static energies
M1=M;
%
%   Step forward in time using leap-frog scheme with Asselin filter
%
N=1+ceil(endtime/dt); % Upper bound on number of points in time
Mg=zeros(1,N);
Dg=zeros(1,N);
Dsg=zeros(1,N);
Msg=zeros(1,N);
hg=zeros(1,N);
tg=zeros(1,N);
%
h1=1;
h2=1;
t=0;
n=0;
%
while t <= endtime && D > 1e-4 && h0*h2 < hD && h1 > 0
    n=n+1;
    t=t+dt;
    Phi=Phi0*sin(pi*t)-Phi0/pi; % Diurnally varying radiative forcing, mean =0 over 24 hours
    %
    [Ds,Ms]= SurfaceT2(alpha,D0*(1-D),D0*(1+M),Q1 );  % Diagnose surface dry and moist static energies
    Ds=Ds/D0-1;  % Normalize output from SurfaceT2
    Ms=Ms/D0-1;
    %
    h3=h1+2.*dt*A*chi*(Ds+D)/D;    % Step h forward in time
    h2=h3+afilt*(h1+h3-2*h2);      % Apply Asselin filter
    h1=h2;
    %
    M=(M1+(Phi0/pi)*(1-cos(pi*t)))/h2;     % Diagnose M
    D=D1*h2^(-(1+A)/A);                    % Diagnose D
    D=real(D);
    %
    % Dimensional vectors for graphics
    %
    hg(n)=h0*h2;
    Dg(n)=D0*(1-D);
    Mg(n)=D0*(1+M);
    Dsg(n)=D0*(1+Ds);
    Msg(n)=D0*(1+Ms);
    tg(n)=6+t*12;
    %
end    
%
% Truncate graphics vectors
%
hg(n+1:N)=[];
tg(n+1:N)=[];
Dg(n+1:N)=[];
Dsg(n+1:N)=[];
Mg(n+1:N)=[];
Msg(n+1:N)=[];
%
Bowen=(Dsg-Dg)./((alpha+1e-6)*(Msg-Mg-Dsg+Dg));   % Bowen Ratio
%
% Calculate CIN and CAPE
%
z=0:5:1.5e4;  % Altitude matrix for CIN and CAPE calculations
%
[~,jmin]=min(abs(z-hD)); % Level above which dry static energy increases
                         % linearly with falling temperature
%
qz=max(size(z));
hlfc=zeros(1,n);
hpbl=zeros(1,n);
CIN=zeros(1,n);
CAPE=zeros(1,n);
%
for i=1:n   % Cycle through time steps
    %
    Tmp=(D0-g*z)/cp;            % Temperature of free troposphere below hD
    p=p0*(cp*Tmp/D0).^(cp/287); % Pressure
    Tc=Tmp-273.15;
    es=6.112*exp(17.67*Tc./(243.5+Tc));
    qs=0.622*es./(p-es);
    %
    Ms2=D0+Lv*qs;               % Saturation moist entropy below hD
    %
    % Ratio of moist to dry adiabatic lapse rates
    %
    Gam=(1+Lv*qs./(R*Tmp))./(1+Lv^2.*qs./(cp*Rv*Tmp.^2));
    %
    %  Find level of free convection
    %
    jlfc=qz;
    for j=qz-1:-1:1
        if Mg(i) <= Ms2(j) && Mg(i) > Ms2(j+1)
            jlfc=j;
            break
        end
    end    
    jlfc=min(jlfc,qz-1);
    hlfc(i)=z(jlfc);
    Tlfc=Tmp(jlfc);
    %
    % Find nearest point in vector z to boundary layer top
    %
    [~,jpbl]=min(abs(z-hg(i)));
    %
    % Now calculate convective inhibition (CIN)
    %
    % First calculate LCL of surface air
    %
    Tsfc=Dg(i)/cp;      
    TCsfc=Tsfc-273.15;
    esfc=6.112*exp(17.67*TCsfc/(243.5+TCsfc));
    qssfc=0.62*esfc/(p0-esfc);
    qsfc=(Mg(i)-Dg(i))/Lv;
    rh=qsfc/qssfc;
    Tstar=55+1/(1/(Tsfc-55)-log(rh)/2840);  % Bolton formula for T at LCL
    zlcl=(Tsfc-Tstar)*cp/g;                 % Altitude of LCL
    [~,jlcl]=min(abs(z-zlcl));              % LCL index
    %
    % Calculate contribution to CIN from saturated ascent above LCL
    % 
    if jlcl < jlfc
      CIN(i)=CIN(i)+sum(Gam(jlcl:jlfc).*(Ms2(jlcl:jlfc)-Mg(i)).* ...
          log(Tmp(jlcl:jlfc)./Tmp(jlcl+1:jlfc+1)));
    end  
    % 
    % Calculate contribution to CIN from dy ascent between PBL top and LCL
    %
    if jpbl < jlcl
        jtop=min(jlcl,jlfc);
        CIN(i)=CIN(i)+sum(((D0-Dg(i))).*log(Tmp(jpbl:jtop)./Tmp(jpbl+1:jtop+1)));
    end    
    %
    % Calculate CAPE
    %
    for j=jlfc:qz
        if Tmp(j) > Tanvil
            Dvar=D0;
            if j > jmin
                Dvar=D0+max((Mg(i)-D0)*(Tmp(jmin)-Tmp(j))/(Tmp(jmin)-Tanvil),0);
            end    
            CAPE(i)=CAPE(i)+Gam(j)*(Mg(i)-Dvar-Lv*qs(j))*(Tmp(j-1)-Tmp(j))/Tmp(j);
        end    
    end    
end 
CIN=max(CIN,0);
CAPE=CAPE-CIN;
CAPE=max(CAPE,0);
CIN(CAPE==0)=NaN;   % Don't plot CIN where it cannot be defined
%
% Graph solutions
%
nmax=n;
titext={['Ts = ',num2str(T0),' K,  F = ',num2str(Q),' W m^{-2}']; ...
    ['\alpha = ',num2str(alpha),',  h_0 = ',num2str(h0),' m,  V_{sfc}= ', num2str(vs), ' m s^{-1}']};
%
ng=3;
while ng > 0
    ng=menu('Graphs','Stop','PBL Depth','Static Energies','CAPE and CIN','Bowen Ratio', ...
        'Bowen Ratio vs. PBL Depth','Soundings');
    if ng == 1
        return
    elseif ng == 2
        g2=plot(tg(1:nmax),hg(1:nmax),'b','linewidth',2);
        set(gca,'fontsize',12,'fontweight','bold')
        set(gca,'xlim',[5.5 18.5])
        xlabel('Local Time')
        ylabel('PBL Depth (m)')
        title(titext)
    elseif ng == 3
        Dtemp=zeros(1,max(size(tg)))+D0;
        g3=plot(tg,Dsg,'r',tg,Dg,'b',tg,Dtemp,'--b',tg, Mg,'g','linewidth',2);
        set(gca,'fontsize',12,'fontweight','bold')
        set(gca,'xlim',[5.5 18.5])
        xlabel('Local Time')
        ylabel('Static Energies (J/Kg)')
        legend('Surface D','D','D_0','M','location','best')
        title(titext)
    elseif ng == 4
        g4=plot(tg(1:nmax),CAPE(1:nmax),'b',tg(1:nmax),CIN(1:nmax),'r','linewidth',2);
        set(gca,'fontsize',12,'fontweight','bold')
        set(gca,'xlim',[5.5 18.5])
        ylim=get(gca,'ylim');
        ylim(1)=0;
        set(gca,'ylim',ylim)
        xlabel('Local Time')
        ylabel('CAPE and CIN (J/Kg)')
        legend('CAPE','CIN','location','best')
        title(titext)
    elseif ng == 5
        g5=plot(tg(1:nmax),Bowen(1:nmax),'b','linewidth',2);
        set(gca,'fontsize',12,'fontweight','bold')
        set(gca,'xlim',[5.5 18.5])
        xlabel('Local Time')
        ylabel('Bowen Ratio')
        title(titext)
    elseif ng == 6
        g6=plot(hg(1:nmax),Bowen(1:nmax),'b','linewidth',2);
        set(gca,'fontsize',12,'fontweight','bold')
        xlabel('PBL Depth')
        ylabel('Bowen Ratio')
        title(titext)
    elseif ng == 7
        hour=1;
        while hour > 0 && ~isempty(hour)
            answer = inputdlg('Hour to plot sounding at (zero to stop)','Sounding Hour'); 
            hour=str2double(answer);
            if isempty(hour)
                hour=0;
            end    
            if hour == 0
                break
            end    
            i=1+round(n*(hour-6)/(tg(n)-6));
            i=max(i,1);
            i=min(i,n);
            [~,jpbl]=min(abs(z-hg(i)));
            Temp=zeros(1,qz);qmix=zeros(1,qz)+0.15*qs;
            Temp(1:jpbl)=(Dg(i)-g*z(1:jpbl))/cp;
            Temp(jpbl+1:jmin)=(D0-g*z(jpbl+1:jmin))/cp;
            Dvar=D0+max((Mg(i)-D0)*(z(jmin+1:qz)-z(jmin))/(16000-z(jmin+1)),0);
            Temp(jmin+1:qz)=(Dvar-g*z(jmin+1:qz))/cp;
            %
            qmix(1:jpbl)=(Mg(i)-Dg(i))/Lv;
            %
            success=skewT(p,Temp-273.15,qmix);
            %set(gca,'yscale','linear')
            xlim=get(gca,'xlim');
            xlim(1)=-45;
            set(gca,'xlim',xlim)
            tit7=['Hour ',num2str(hour),',  CAPE= ',num2str(round(CAPE(i))),' J/Kg,  CIN= ', ...
                num2str(round(CIN(i))),' J/Kg'];
            title(tit7)
        end   
    end
end    
%
function [ Ds, Ms ] = SurfaceT2(alpha,D,M,Q )
%
%   Finds dry static energy of a zero-heat-capacity surface given
%   the net radiative input (Q) normalized by surface wind, Ck, and
%   density, moist coefficient (alpha) of surface (0-1), and dry and moist  
%   static energies of the boundary layer (D and M). Newton-Raphson method.
%   July, 2016
%--------------------------------------------------------------------
thresh=1e-2;  % Convergence threshold
p0=1005;      % Surface pressure (hPa)
cp=1005;      % Heat capacity
Lv=2.5e6;     % Latent heat
Rv=461;       % Gas constant for water vapor
%
Ds1=0.99*D;   % Initial guess of surface dry static
%
Ds2=2*Ds1;
Ds1old=Ds1;
%
coef=alpha*Lv^2*cp/Rv;
%
while abs(Ds2-Ds1old) >= thresh
    Ts=Ds1/cp;
    Tsc=Ts-273.15;
    es=6.112*exp(17.67*Tsc./(243.5+Tsc));
    qs=0.622*es./(p0-es);
    Ms=Ds1+Lv.*qs;
    %
    Ds2=Ds1-((1-alpha)*(Ds1-D)+alpha*(Ms-M)-Q)/(1+coef/Ds1^2);
    Ds1old=Ds1;
    Ds1=Ds2;
    %disp(num2str(Ds2-Ds1old))
end
Ds=Ds2;
Ms=real(Ms);
Ds=real(Ds);
end
%
function success=skewT(pz,tz,qz)
%
% Calculates and displays a sounding on a skew-T chart. Pseudo-adiabats
% from entropy definition by George Bryan (Mon. Wea. Rev., 2008). 
%
% pz = pressure (hPa), tz = temperature (C)
% qz=relative humidity (%) if lowest level value is in range 0.5-100; 
%   otherwise qz is assumed to be specific humidity (g/g).
%
% Modified December, 2021
%
%------------------------------------------------------------------------
%
%   First check to see if qz is relative or specific humidity
%
qsize=max(size(qz));
imin=100;
for i=1:qsize
    if ~isnan(qz(i)) && qz(i) ~= 0
        imin=min(i,imin);
    end
end  
if imin == 100
    disp(' ')
    disp('No humidity data, quitting')
    return
end    
if qz(imin) > 0.5
    ez=6.112.*exp(17.67.*tz./(243.5+tz));
    qz=0.01*qz.*0.622.*ez./(pz-ez);
end
%
chi=log(pz.*qz./(6.112.*(0.622+qz)));
tdz=243.5.*chi./(17.67-chi);
%
p=1050:-25:50;
pplot=transpose(p);
t0=-48:2:50;
tmin=(273.15+min(t0))*(0.001*min(p)).^0.287-273.15;
v=10*floor(tmin/10):10:10*ceil(max(t0)/10);
[ps1,ps2]=size(p);
ps=max(ps1,ps2);
[ts1,ts2]=size(t0);
ts=max(ts1,ts2);
Lvcp=2.555e3./1005;
Rdcp=287/1005;
tvem=zeros(ts,ps);
thet=zeros(ts,ps);
thetaea=zeros(ts,ps);
q=zeros(ts,ps);
for i=1:ts
   for j=1:ps
      tv=t0(i)+30.*log(0.001.*p(j));
      TK=(273.15+tv);
      Tc=TK-273.15;
      es=6.112.*exp(17.67.*Tc./(243.5+Tc));
      qg=0.622.*es./(p(j)-es);
      q(i,j)=1000*qg;
      thet(i,j)=TK.*(1000./p(j)).^Rdcp;
      thetaea(i,j)=thet(i,j).*exp(Lvcp.*q(i,j)./TK);
      tvem(i,j)=tv;
   end
end
t0=transpose(t0);
temp=transpose(tvem);
theta=transpose(thet);
thetae=transpose(thetaea);
qs=transpose(sqrt(q));
[~]=contour(t0,pplot,temp,v,'k');
set(gca,'ylim',[50,max(pz)+10])
hold on
set(gca,'ytick',1000:100:100)
set(gca,'yscale','log','ydir','reverse')
set(gca,'fontweight','bold')
set(gca,'ytick',100:100:1000)
set(gca,'ygrid','on')
hold on
[~]=contour(t0,pplot,theta,24,'b');
[~]=contour(t0,pplot,qs,24,'g');
[~]=contour(t0,pplot,thetae,24,'r');
tvz=(tz+273.15)-273.15;
tzm=tvz-30.*log(0.001.*pz);
tdzm=tdz-30.*log(0.001.*pz);
h=plot(tzm,pz,'r',tdzm,pz,'g');
set(h,'linewidth',2)
hold off
xlabel('Density Temperature (C)','fontweight','bold')
ylabel('Pressure (mb)','fontweight','bold')
success='yes';
end