moist_ape/ 0000755 0013561 0013560 00000000000 13505422437 012061 5 ustar pog ogorman moist_ape/simple_sat_vapor_pressure.m 0000644 0013561 0013560 00000000727 13505422331 017535 0 ustar pog ogorman function simple_sat_vapor_pressure = simple_sat_vapor_pressure(temp) % % Calculates the saturation vapor pressure assuming a constant % latent heat % error(nargchk(1, 2, nargin)) % parameters temp_ref = 273.16; eps = 0.622; e_00 = 610.78; Rd = PARS('gas_constant'); Rv = Rd/eps; l_cond = PARS('l_cond'); % find the saturation vapor pressure simple_sat_vapor_pressure = e_00*exp(l_cond/Rv*(1.0/temp_ref-1.0./temp)); moist_ape/test_moist_mape.m 0000644 0013561 0013560 00000001100 13223212473 015415 0 ustar pog ogorman clear % Test calculation of moist available potential energy for an artifical mean state that varies with latitude and pressure. p = [1000:-50:100]*100; % pressure in Pa lat = [0:5:90]; % latitude in degrees rhum = ones(length(p), length(lat))*0.7; % relative humidity for k=1:length(p) for j=1:length(lat) temp(k,j) = 300-(1e5-p(k))/1e5*70 - lat(j)/90*40; % temperature in K end end do_dry = 0; % include moisture? do_nonconvective = 0; % allow convection? moist_mape_ipcc(p, lat, temp, rhum, ones(length(lat))*1e5, 0, 90, do_dry, do_nonconvective) moist_ape/moist_mape_ipcc.m 0000644 0013561 0013560 00000004242 13223211227 015362 0 ustar pog ogorman function ape = moist_mape_ipcc(pfull, lat, temp_avg, rhum_avg, ps, min_lat, max_lat, do_dry, do_nonconvective) % Example of how to calculate mean ape (MAPE) of a moist atmosphere that varies with pressure and latitude (was applied to IPCC model output) % Result is expressed in J/m^2 error(nargchk(9, 9, nargin)) % pfull: pressure (Pa); first element is highest pressure % lat: latitude (degrees) % temp_avg: mean temperature versus pressure and latitude (K) % rhum_avg: mean relative humidity versus pressure and latitude (fraction rather than percent) % ps: surface pressure versus latitude (Pa) % min_lat: minimum latitude (degrees) % max_lat: maximum latitude (degrees) % do_dry: ignore moisture % do_nonconvective: "non-convective" moist ape discussed in O'Gorman, PNAS, 2010 that does not allow the reordering of the pressure of air parcels originating at a given latitude if do_dry rhum_avg = zeros(size(rhum_avg)); % has to be all zeros (no NaNs) end % interpolate to common grid with even spacing % and unique pressure levels (See Lorenz paper) num_lev = 40; num_lat = 40; pA = 0.05e5; pB = 1.0e5; % exclude subsurface, but try not to lose lowest layer of atmosphere % (important given coarse vertical resolution of IPCC models) dp = pfull(1)-pfull(2); if dp<0 error('unexpected pressure ordering') end for j=1:length(lat) temp_avg(pfull>ps(j)+dp,j)=NaN; end % equal area grid dsin_lat = (sind(max_lat)-sind(min_lat))/(num_lat-1); sin_lat_grid = [sind(min_lat):dsin_lat:sind(max_lat)]; lat_grid = asind(sin_lat_grid); dp = (pB-pA)/num_lev/num_lat; for k=0:num_lev-1 for j=1:num_lat i=k*length(lat_grid)+j; p(i) = pA-dp/2+dp*(k*length(lat_grid) + j); r(i) = interp2(lat, pfull, rhum_avg, lat_grid(j), p(i)); T(i) = interp2(lat, pfull, temp_avg, lat_grid(j), p(i)); lat_parcel(i) = lat_grid(j); end end % exclude NaNs isf = isfinite(T)&isfinite(r); p = p(isf); r = r(isf); T = T(isf); lat_parcel = lat_parcel(isf); disp('divide and conquer') [ape, p_ref] = moist_ape_divide_conquer(T, p, r, do_nonconvective, lat_parcel, lat_grid); ape = ape*(pB-pA)/PARS('gravity'); moist_ape/moist_ape_parcel_move.m 0000644 0013561 0013560 00000004260 13113553272 016572 0 ustar pog ogorman function [Tvf, dh] = parcel_move(T, p, r, e, Tc, pc, wt, dh_to_sat, pf, thermo_type) % move parcel to pressure pf % see moist_ape_thermo for variable definitions error(nargchk(10, 10, nargin)) switch lower(thermo_type) case 'full' cpw = PARS('cp_v'); % specific heat capacity of water vapor case 'simple' % idealized GCM thermo is equivalent to having the same heat capacity % for water vapor and dry air cpw = PARS('cp'); otherwise error(['Unknown type, ', thermo_type, 'of thermodynamic formulation.']) end R = PARS('gas_constant'); Rw = PARS('gas_constant_v'); cp = PARS('cp'); eps = R/Rw; unsat = pf>pc; sat = pf<=pc; % if don't reach saturation kappa = (R+wt.*Rw)./(cp+wt.*cpw); Tf_unsat = T.*(pf./p).^kappa; w = eps*e./(p-e); Tvf_unsat = (1+w./eps).*Tf_unsat./(1+wt); dh_unsat = (cp+wt.*cpw).*(Tf_unsat-T)./(1+wt); % equation 12 Lorenz if all(r==0) % dry mape Tvf = Tvf_unsat; dh = dh_unsat; return end % if reach saturation s_sat_0 = moist_ape_sat_entropy(Tc, pc, wt, thermo_type); % iterate Tf_sat = Tc; % initial guess dT = ones(size(T)); T_perturb = 0.1; % used to evaluated ds_dT while any(abs(dT(sat))>0.01) s_sat = moist_ape_sat_entropy(Tf_sat, pf, wt, thermo_type); Tf_sat_perturb = Tf_sat+T_perturb; s_sat_perturb = moist_ape_sat_entropy(Tf_sat_perturb, pf, wt, thermo_type); ds_sat_dT = (s_sat_perturb-s_sat)/T_perturb; dT = (s_sat_0-s_sat)./ds_sat_dT; Tf_sat = Tf_sat+ dT; end [es_at_Tf, L_at_Tf] = moist_ape_es(Tf_sat, thermo_type); w_at_Tf = eps*es_at_Tf./(pf-es_at_Tf); Tvf_sat = (1+w_at_Tf./eps).*Tf_sat./(1+wt); [es_at_Tc, L_at_Tc] = moist_ape_es(Tc, thermo_type); w_at_Tc = eps*es_at_Tc./(pc-es_at_Tc); % based on equation 16 (doesn't require a c) dh_sat = ((cp+wt*cpw).*(Tf_sat-Tc)-... (wt.*(L_at_Tf-L_at_Tc)-... w_at_Tf.*L_at_Tf+w_at_Tc.*L_at_Tc))./(1+wt)+... dh_to_sat; % combine cases dh = dh_unsat.*unsat + dh_sat.*sat; Tf = Tf_unsat.*unsat + Tf_sat.*sat; Tvf = Tvf_unsat.*unsat + Tvf_sat.*sat; moist_ape/moist_ape_divide_conquer.m 0000644 0013561 0013560 00000011416 13223211017 017265 0 ustar pog ogorman function [moist_ape, p_ref]=moist_ape(T, p, r, do_nonconvective, lat_parcel, lat_grid, thermo_type) % Divide and conquer algorithm for moist available potential energy % The divide and conquer algorithm was introduced in Stansifer, O'Gorman and Holt, QJRMS, 143, 288-292, 2017. The python code used for that paper is available at https://github.com/estansifer/mape/. The version given here differs by using different thermodynamic constants and a different treatment of moisture (either 'full' which includes ice or 'simple' which has constant latent heat of vaporization) % The input profiles are assumed to be 1d, with evenly spaced pressure levels, and no liquid water or ice. An example of how to apply this function to a zonal and time mean state that varies with latitude is given in moist_mape_ipcc.m % The input arrays lat_parcel and lat_grid are only needed if do_nonconvective = 1 (the "non-convective" moist ape discussed in O'Gorman, PNAS, 107, 19176-19180, 2010 that does not allow the reordering of the pressure of air parcels originating at a given latitude) % Variable key: % moist_ape moist available potential energy (or gcape) in J/kg % p_ref reference pressure (Pa) % T temperature (K) % p pressure (Pa) % r generalized relative humidity (ratio of vapor pressures when below one, % but a ratio of mixing ratios when above one) (note fraction rather than percent) % s entropy % h enthalpy % w mixing ratio % L latent heat of evaporation or sublimation % e vapor pressure % wt total mixing ratio % lat_parcel latitudes of parcels in the profile % lat_grid the grid of possible latitudes error(nargchk(3, 7, nargin)) % thermo_type: default is full rather than simple if nargin < 7 thermo_type = 'full'; end if nargin<4 do_nonconvective = 0; % default is full moist ape (allows for convection) end if do_nonconvective&nargin<6 error('must supply latitudes if do_nonconvective=1') end if p(2)