%
%  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
%  Modified November 2022 to subtract constant radiative cooling to
%  give zero diurnal mean and to incorporate soil model
%-------------------------------------------------------------------------
%
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=215;% Tropopause temperature
vs=10;     % Surface wind speed (m/s)
Ck=1.5e-3; % Surface sensible and latent heat exchange coefficient
%
%  Soil parameters
%
d1=0.1;    % Normalization surface soil depth (m)
d2=0.2;    % Depth of soil layer (m)
rh2=0.8;   % Initial saturation fraction of deeper soil
rhg=0.8;   % Initial saturation fraction of surface soil
veg=0.3;   % Vegetation cover
stype='Loam';  %  Soil type.  Choices: 'Sand','Loamy sand','Sandy loam',
               % 'Silt loam','Loam','Sandy clay loam','Silty clay loam', 
               % 'Clay loam','Sandy clay','Silty clay','Clay'
%
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)
rhol=1000; % Density of water (kg/m^3)
cp=1005;   % Heat capacity
CV=1e-3;   % Vegetation inverse 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 initial 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;
%
[Ds1,Ms]= SurfaceT2(0.8,D0*(1-D),D0*(1+M),Q1 );  % Diagnose initial surface dry and moist static energies
Ms=max(Ms,Ds1);
Ds1=Ds1/D0-1;  % Normalize output from SurfaceT2
Ms=Ms/D0-1;
%
% Soil properties
%
[wsat, b, Cgsat, p, a, C2ref,C1sat] = soilprop(stype);
%
Cgsat=Cgsat*1e-6;  % I think this corrects error in Table 2 of N & P (1989)
wfl=0.75*wsat;     % Field capacity
w21=rh2*wsat;      % Initial deeper soil volumetric water content
wg1=rhg*wsat;      % Initial surface soil volumetric water content
%
%   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);
D2g=zeros(1,N);
Msg=zeros(1,N);
hg=zeros(1,N);
wgg=zeros(1,N);
w2g=zeros(1,N);
tg=zeros(1,N);
%
h1=1;
h2=1;
emu1=h1*M1;
emu2=emu1;
w22=w21;
wg2=wg1;
Ds2=Ds1;
D21=Ds1;
D22=D21;
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
    %
    % Soil model
    %
    Cg=Cgsat*(wsat/w21).^(b/(2*log(10)));
    CT=1./((1-veg)/Cg+veg/CV);
    wgeq=w21-wsat*a*(w21/wsat)^p*(1-(w21/wsat)^(8*p));
    C1=C1sat*(wsat/wg1)^(0.5*b+1);
    C2=C2ref*(w21/(wsat-w21+0.01));
    alpha=0.5*(1-cos(pi*min(wg1/wfl,1)));
    %
    % Soil forcings
    %
    FDs=cp*CT*rho*h0*(Phi-chi*(Ms-M1))-pi*(Ds1-D21); % Surface temperature
    FD2=0.5*(Ds1-D21);                               % Deeper soil temperature
    Fw2=-(D0*rho*h0*chi/(rhol*Lv*d2))*(Ms-M1-Ds1-D); % Deeper soil moisture
    Fwg=C1*Fw2*d2/d1-0.5*C2*(wg1-wgeq);              % Surface soil moisture

    h3=h1+2.*dt*A*chi*(Ds2+D)/D;         % Step h forward in time
    h2=h3+afilt*(h1+h3-2*h2);            % Apply Asselin filter
    h1=h2;
    %
    D=D1*h2^(-(1+A)/A);                  % Diagnose D
    D=real(D);
    %
    emu3=emu1+2*dt*chi*(Ms-M1);          % Step forward product of PBL depth and MSE
    emu2=emu3+afilt*(emu1+emu3-2*emu2);
    emu1=emu2;
    M1=emu1/h1;                          % Diagnose MSE
    %
    Ds3=Ds1+2*dt*FDs;                    % Step forward soil variables
    D23=D21+2*dt*FD2;
    wg3=wg1+2*dt*Fwg;
    w23=w21+2*dt*Fw2;
    wg3=min(wg3,wsat);
    w23=min(w23,wsat);
    %
    Ds2=Ds3+afilt*(Ds1+Ds3-2*Ds2);       % Apply asselin filter to soil variables
    D22=D23+afilt*(D21+D23-2*D22);
    wg2=wg3+afilt*(wg1+wg3-2*wg2);
    w22=w23+afilt*(w21+w23-2*w22);
    Ds1=Ds2;
    D21=D22;
    wg1=wg2;
    w21=w22;
    %
    Tc=D0*(1+Ds2)/cp-273.15;             % Diagnose saturated surface soil moisture
    es=6.112*exp(17.67*Tc/(243.5+Tc));
    qs=0.622*es/(p0-es);
    Ms=(D0*Ds2+alpha*Lv*qs)/D0;
    %
    %
    % Dimensional vectors for graphics
    %
    hg(n)=h0*h2;
    Dg(n)=D0*(1-D);
    Mg(n)=D0*(1+M1);
    Dsg(n)=D0*(Ds2+1);
    D2g(n)=D0*(D22+1);
    Msg(n)=D0*(1+Ms);
    wgg(n)=wg2;
    w2g(n)=w22;
    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)=[];
D2g(n+1:N)=[];
Mg(n+1:N)=[];
Msg(n+1:N)=[];
wgg(n+1:N)=[];
w2g(n+1:N)=[];
%
Bowen=(Dsg+Dg)./(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','Soil moisture','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,D2g,'--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','Soil 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
        wsatg=zeros(1,nmax)+wsat;
        g6=plot(tg(1:nmax),wgg(1:nmax),'b',tg(1:nmax),w2g(1:nmax),'--b', ...
            tg(1:nmax),wsatg,'g','linewidth',2);
        set(gca,'fontsize',12,'fontweight','bold')
        xlabel('Local Time')
        ylabel('Volumetric Water Content')
        legend('Surface soil moisture','Deeper soil moisture', ...
            'Saturated soil moisture','location','best')
        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 [porosity, b, CGsat, p, a, C2ref,C1sat] = soilprop(stype)
%
% This function returns tabulated soil parameters from Noilhan and Planton
% (1989) and porosity from
% https://hess.copernicus.org/articles/22/4935/2018/hess-22-4935-2018-t04.pdf 
% November, 2022
% 
if strcmpi(stype,'Sand')
    porosity=0.437;
    b=4.05;
    CGsat=3.222;
    p=4;
    a=0.387;
    C2ref=3.9;
    C1sat=0.082;
elseif strcmpi(stype,'Loamy sand')
    porosity=0.437;
    b=4.38;
    CGsat=3.057;
    p=4;
    a=0.404;
    C2ref=3.7;
    C1sat=0.098;
elseif strcmpi(stype,'Sandy loam')
    porosity=0.453;
    b=4.9;
    CGsat=3.56;
    p=4;
    a=0.219;
    C2ref=1.8;
    C1sat=0.132;
elseif strcmpi(stype,'Silt loam')
    porosity=0.501;
    b=5.3;
    CGsat=4.418;
    p=6;
    a=0.102;
    C2ref=0.8;
    C1sat=0.153;    
elseif strcmpi(stype,'Loam')
    porosity=0.463;
    b=5.39;
    CGsat=4.111;
    p=6;
    a=0.148;
    C2ref=0.8;
    C1sat=0.191;
elseif strcmpi(stype,'Sandy clay loam')
    porosity=0.398;
    b=7.12;
    CGsat=3.67;
    p=6;
    a=0.135;
    C2ref=0.8;
    C1sat=0.213;
elseif strcmpi(stype,'Silty clay loam')
    porosity=0.471;
    b=7.75;
    CGsat=3.593;
    p=8;
    a=0.127;
    C2ref=0.4;
    C1sat=0.385;    
elseif strcmpi(stype,'Clay loam')
    porosity=0.464;
    b=8.52;
    CGsat=3.995;
    p=10;
    a=0.084;
    C2ref=0.6;
    C1sat=0.227;
elseif strcmpi(stype,'Sandy clay')
    porosity=0.43;
    b=10.4;
    CGsat=3.058;
    p=8;
    a=0.139;
    C2ref=0.3;
    C1sat=0.421;
elseif strcmpi(stype,'Silty clay')
    porosity=0.479;
    b=10.4;
    CGsat=3.729;
    p=10;
    a=0.075;
    C2ref=0.3;
    C1sat=0.375;    
elseif strcmpi(stype,'Clay')
    porosity=0.475;
    b=11.4;
    CGsat=3.6;
    p=12;
    a=0.083;
    C2ref=0.3;
    C1sat=0.342;    
end
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