%% Betts2000_EBLmodel.m: % * Script to solve/plot T and q for Equilibrium Boundary Layer % * based on: Betts, 2000, "Idealized Model for Equilibrium Boundary Layer % over Land", Journal of Hydrometeorology, Volume 1, pp. 507-523 % * source referred to hereafter as B00 % % * Coded by Tim Cronin (twc) 1/22/2013-1/25/2013 % % * notes: % * 1) B00 states on p. 511 that "the system of Eqs. (1), (2), (6), (7), % (8), (9), (10), (15a), (18), and (21) can be solved iteratively, ...", % but here I attempt to solve with a mix of prognostic equations and % diagnostic equations, with explicit leapfrog timestepping for the 3 % prognostic equations: % * 1a) Cs Ts_dot = Qs - SHF - LHF % * 1b) Cp P_H/g thetab_dot = SHF + Cp Omega_t/g (thetat - thetab) - Cp P_H/g (Qbrad+Qbevap) % * 1c) P_H/g qb_dot = LHF/Lv + Omega_t/g (qt-qb) + Cp P_H / (g Lv) Qbevap % * Then, the diagnostic equations allow for updating of P_H, thetat, qt, % and Omega_t % * 2) On notation: I use a subscript '_b' or letter 'b', for 'b'oundary or % su'b'cloud layer; B00 uses '_M' or 'M' for 'M'ixed layer % * 3) Reference to 'B02' as a source for an equation indicates: Betts et % al, 2002, "Surface diurnal cycle and boundary layer structure over % Rondonia during the rainy season", JGR, Vol. 107, No. D20, 8065 % % %% Set experiment(s) to run % string parsing on 'simcase' decides what simulations to run, and what % should be plotted % % options for reproduction of Betts, 2000 plots: 'Fig1a', 'Fig1b', 'Fig3-NW', % 'Fig3-W', 'Fig3-SW', 'Fig3-NE', 'Fig3-E', 'Fig3-SE', % -- these all attempt to reproduce the like-numbered (/located) figures from B00 % % options: 'twc-gv', 'twc-Qs', or 'twc-ga' for single-factor experiment across a range % of gv, Qs, or ga values, with comparison to theory from Cronin & % Emanuel 2013, in prep. % %simcase = 'Fig'; simcase = 'Fig1a'; %simcase = 'twc-ga'; %% Initialize Parameters % Reference parameters from Table 1 of B00 % name value units description ga_0 = 0.025; % m/s aerodynamic conductance gv_0 = 1/60; % m/s vegetative conductance to water vapor kent = 0.2; % - entrainment parameter relating ML-top theta-flux to sfc. buoyancy flux Gamma_0 = 6e-4; % K/Pa stability |d \theta/d p| above ML top Qs_0 = 150; % W/m^2 surface net radiation Qbrad0 = 3; % K/day ML radiative cooling rate Qbevap0 = 0; % K/day ML cooling rate by evaporating rain % renormalize Qbrad and Qbevap to K/s Qbrad0 = Qbrad0/86400; Qbevap0 = Qbevap0/86400; % Other reference parameters (from equations/text of B00, and new to this code) % name value units description Cp = 1006; % J/kg/K specific heat capacity of dry air (B13-pc) Rd = 287; % J/kg/K gas constant for dry air RdCp = 0.286; % - Ratio of Rd/Cp Lv = 2.501e6; % J/kg latent heat of vaporization of water at 0 C (B13-pc) g = 9.81; % m/s^2 gravitational constant Th_00 = 303; % K constant potential temperature offset in B00 eqn. (12) PH_00 = 60e2; % Pa reference P_H offset in B00 eqn. (12) P_T0 = 100e2; % Pa reference subsaturation just above ML top (B00 eqn. (13)) eps = 0.622; % - ratio of gas constants / MWs for (water vapor):(dry air) Cs = 1e6; % J/m^2/K surface heat capacity (not in B00) ps = 940e2; % Pa surface pressure p0 = 1e5; % Pa reference pressure for definition of theta nt = 15000; % - number of timesteps dt = 200; % s timestep afilt = 0.1; % - Asselin filter parameter % declare the 'anonymous function' qsat_B00 % saturation mixing ratio from Betts et al, 2002, JGR 107-D20 qstar_B00=@(T_C,p_Pa) (0.622/(0.0016361*p_Pa*exp(-17.67*T_C/(243.5+T_C))-1)); %% Create simulation-specific arrays % (for Gamma, P_T, Qs, Qbrad, Qbevap; gv, Qs, ga) if strcmp(simcase(1:3),'twc') % multiplier factors for gv, Qs, or ga m_facs = (0.7:0.05:1.3); nsims = length(m_facs); m_1s = ones(1,nsims); Gam_arr = [4 5 6 7]*10^(-4); % K/Pa --stability array P_T_arr = [1 1 1 1]*P_T0; % Pa --subsaturation array Qbrad_arr = [1 1 1 1]*Qbrad0; % K/s --radiative cooling array Qbevap_arr = [1 1 1 1]*Qbevap0; % K/s --evaporative cooling array nsens = 4; if strcmp(simcase,'twc-gv') gv_arr = gv_0*m_facs; Qs_arr = Qs_0*m_1s; ga_arr = ga_0*m_1s; elseif strcmp(simcase,'twc-Qs') gv_arr = gv_0*m_1s; Qs_arr = Qs_0*m_facs; ga_arr = ga_0*m_1s; elseif strcmp(simcase,'twc-ga') gv_arr = gv_0*m_1s; Qs_arr = Qs_0*m_1s; ga_arr = ga_0*m_facs; end elseif strcmp(simcase(1:3),'Fig') gv_arr = 1./[60,80,100,150,300,500,700,900]; Qs_arr = Qs_0*ones(1,length(gv_arr)); ga_arr = ga_0*ones(1,length(gv_arr)); if (strcmp(simcase,'Fig1a')||strcmp(simcase,'Fig3-NW')||... strcmp(simcase,'Fig3-W')||strcmp(simcase,'Fig3-SW')) %these plots require Gamma to be varied Gam_arr = [4 5 6 7]*10^(-4);% K/Pa Stability array P_T_arr = [1 1 1 1]*P_T0; Qbrad_arr = [1 1 1 1]*Qbrad0; Qbevap_arr = [1 1 1 1]*Qbevap0; nsens = 4; psyms = {'k-.'; 'k:'; 'k-'; 'k--'}; % plot symbols elseif (strcmp(simcase,'Fig1b')||strcmp(simcase,'Fig3-NE')||... strcmp(simcase,'Fig3-E')||strcmp(simcase,'Fig3-SE')) %these plots require P_T to be varied Gam_arr = [1 1 1]*Gamma_0; P_T_arr = [60 100 140]*1e2; Qbrad_arr = [1 1 1]*Qbrad0; Qbevap_arr = [1 1 1]*Qbevap0; nsens = 3; psyms = {'k:'; 'k-'; 'k--'}; % plot symbols elseif strcmp(simcase,'Fig5-W') % this plot requires Qbrad to be varied Gam_arr = [1 1 1]*Gamma_0; P_T_arr = [1 1 1 ]*P_T0; Qbrad_arr = [1/3 2/3 1]*Qbrad0; Qbevap_arr = [1 1 1]*Qbevap0; nsens = 3; psyms = {'k:'; 'k--'; 'k-'}; % plot symbols elseif strcmp(simcase,'Fig5-SW') % this plot requires Qbevap to be varied Gam_arr = [1 1 1 1]*Gamma_0; P_T_arr = [1 1 1 1]*P_T0; Qbrad_arr = [1 1 1 1]*Qbrad0; Qbevap_arr = [0 1/3 2/3 1]*Qbrad0; nsens = 4; psyms = {'k-'; 'k:'; 'k--'; 'k-.'}; % plot symbols else % do only 1 sensitivity test Gam_arr = Gamma_0; P_T_arr = P_T0; Qbrad_arr = Qbrad0; Qbevap_arr = Qbevap0; nsens = 1; end end nsims = length(gv_arr); %% Loop over values of Gamma, P_T, Qbrad, or Qbevap for isens=1:nsens % set control parameters to their appropriate values Gamma = Gam_arr(isens); P_T = P_T_arr(isens); Qbrad = Qbrad_arr(isens); Qbevap = Qbevap_arr(isens); % initialize plot/output variables to zero for new simulation P_LCL = zeros(length(gv_arr),1); thetab_out = zeros(length(gv_arr),1); qb_out = zeros(length(gv_arr),1); Tb_out = zeros(length(gv_arr),1); Ts_out = zeros(length(gv_arr),1); LHF_out = zeros(length(gv_arr),1); SHF_out = zeros(length(gv_arr),1); delth_out = zeros(length(gv_arr),1); delq_out = zeros(length(gv_arr),1); omega_out = zeros(length(gv_arr),1); %% Loop over single factor values for ivar=1:nsims gv = gv_arr(ivar); Qs = Qs_arr(ivar); ga = ga_arr(ivar); %% Variable initialization % Prognostic variables: 3-element arrays are for Leapfrog/Asselin timestepping % Ts = (zeros(3,1)+Th_00)*(ps/p0)^(RdCp); % surface temperature thetab = zeros(3,1)+Th_00-1; % ML potential temperature Tb = thetab(2)*(ps/p0)^(RdCp); % surface air temperature qstar_b = qstar_B00(Tb-273.16,ps); % svp of surface air RH_b = 0.7; qb = zeros(3,1)+RH_b*qstar_b; % ML mixing ratio -- initialize at 70% RH % prognostic variable tendencies Ts_dot = 0; thetab_dot = 0; qb_dot = 0; % other variables that are used in diagnostic equations only % are declared within the relevant diagnostic equation block %% Integration in time for it=1:nt % Diagnostic equations (D #*) % (D #0) Diagnose surface air density using ideal gas law % (B13-pc) Tb = thetab(2)*(ps/p0)^(RdCp); rho0 = ps/(Rd*Tb*(1+0.608*qb(2))); % ---------- end (D #0) ---------- % (D #1) Diagnose P_H based on thetab and qb using B00 Eqns. (14) and (22) qstar_b = qstar_B00(Tb-273.16,ps); RH_b = qb(2)/qstar_b; % 3 following lines are from B02 Tsat_b = 55+2840/(3.5*log(thetab(2))-log(1000*qb(2)/(0.622+qb(2)))-4.805); psat_b = p0*(Tsat_b/thetab(2))^3.4965; P_H = ps-psat_b; %A = eps*Lv/(2*Cp*Tb); %P_H = ps*(1-RH_b)/(A+(A-1)*RH_b); % N.B. units of P_H are Pascals % check that P_H is +ve and does not exceed 0.9*p0 if P_H<0 P_H = PH_00; elseif P_H>ps P_H = 0.9*ps; end % ---------- end (D #1) ---------- % (D #2) Diagnose thetat based on B00 Eqn. (12) thetat = Th_00 + Gamma*(P_H - PH_00); % check for static stability at ML top if thetat