continuum1000 – Continuum absorption coefficient of water vapor (self-broadened) around 1000 cm-1, from Pierrehumbert 2010, p260, calculates capture cross section, kh2o in m2/kg for each wavenumber in v between minv, maxv, at each temperature in T value of kh2o will be zero outside of minv – maxv

define_atmos_0_5 – creates a high resolution atmosphere using the surface temperature, lapse rate, ideal gas laws (pressure with height also depends on the temperature) and some imposed assumptions about the stratosphere – the pressure of the tropopause and the temperature profile above it (isothermal currently)

Update Jan 24th of the HITRAN data files – these were .par files, but WordPress only allows certain file types to be uploaded. To use them with the Matlab program, save and rename each to .par. Each file is the line data for one molecule, with the number corresponding to its HITRAN number (see the papers in the references).

The code for the main function is pasted inline at the end of the article for casual perusal. But WordPress doesn’t format it particularly well, so load it up and look at it in a text editor if you want to understand it.

In MATLAB it’s easy to see the code vs comments (%..) because the color is changed. Not so easy in a text editor unless you have a clever one.

The variables on the right – dv, numz, mix, blp, BLH, FTH, contabs, Ts, lapsekm, minp, tstep, nt, ploton – are the inputs to the function. Some are vectors (MATLAB is fine with x being a number or a 1d vector or a 2d matrix etc).

The function is called either from the command line or from another program adjusting parameters each time and re-running.

‘…’ at the end of line links to the next line (to keep readability).

Vectors can be multiplied together to get the dot product or to get a1b1, a2b2, a3b3 as a new vector. It’s the hardest thing to get used to in Matlab when you have been used to indexing everything, e.g., for i=1 to n; c(i)=a(i)*b(i); end – ask if you want to know what a specific line does.

ii) Using surface pressure, minp (pressure of TOA) and number of layers (numz-1) – create a set of layers of roughly equal pressure using data from the high res atmosphere – get p, T, rho for the midpoint and the boundaries

Calculate transmission through each layer plus emission from each layer to get spectrum at each boundary in both directions.

Transmissivity = exp(-tau)

Absorptivity= 1-exp(-tau); Beer Lambert law

Emissivity = Absorptivity; from Kirchhoff

The above values are for each wavenumber.

6. If more than tstep>2 the change in temperature will be calculated at each layer from change in energy. dT = (change in upward flux through layer + change in downward flux through layer) x timestep / (rho x cp x thickness of layer). It’s just the first law of thermodynamics.

Then convective adjustment checks whether the lapse rate is too great, which is unstable, and if so, adjusts the temperature back to the lapse rate. This critical lapse rate can be different at different heights, currently it is set to one value.

Code

==== v 0.10.4 ====== added to article Jan 30th ===

function [ vt tau fluxu fluxd neth z dz p rho mix zb pb Tb …
Tinit T Ts dE dEs dCE radu radd emitu sheat surfe TOAe rf radupre raddpre …
radupost raddpost TOAf TOAtr alr] =…
HITRAN_0_10_4(dv, std_atm, numz, ratiodp, mol, mixc, newmix, ison, blp, BLH, FTH, STaH, contabs,…
linewon, Tsi, lapsekm, lapseikm, ptropo, minp, ocd, src, tstep, nt, ntc, ploton)
% Radiative transfer through atmosphere, up and down
% Uses the HITRANS database for absorption and emission lines with
% Lorenzian broadening, not yet using Voigt profile; Pierrehumbert 2010 for
% continuum absorption.
% Generally creates a standard atmosphere of pressure, temperature to
% calculate up and down spectral fluxes
% Can work from starting point to recalculate new temperature profile via
% “radiative convective” model using the defined lapse rate
%
%
% v0.1 – Basic calculations to get out of the starting blocks
% v0.2 – Tidy up and allow a “slab” to be at any pressure and temperature
% with line width determined by the pressure & temperature
% v0.3 – Allow multiple layers
% v0.4 – Orphan branch to evaluate wings of lines vs weak lines
% v0.5 – More GHGs, select which ones to pull data in
% v0.6 – Orphan branch not used here, adds CFC-11 & 12 which was difficult
% because the minor GHGs are in a different database
% v0.7 – Dec 31st 2012, significant update to calc emission and absorption
% through atmosphere (rather than just calc transmissivity) as a
% radiative-convective model
% v0.7.1- Added the energy absorbed/emitted by layer/wavenumber from
% rte_0_5_3
% v0.8 – Set boundary layer top pressure so changes to no of layers doesn’t
% alter calcs (BLH can be different from FTH), fix stratosphere temp
% in lieu of solar heating (via iso_strato parameter)
% v0.9 – Introduce continuum absorption of water vapor
% v0.9.1- Adding some interim calcs like effective emission height, % of
% emitted flux from each height that makes up TOA; likewise for DLR;
% v0.9.2- Clean up code comments for publishing, fix the calculation of
% layers just above the boundary layer (as noted by Frank), change
% surface pressure to 1013hPa, fixed convective adjustment to include
% the surface to first atmospheric layer, add different lapse rate
% for each layer
% v0.9.3- Allow stratosphere temp to move, but stay isothermal
% v0.9.4- Orphan version
% v0.9.5- Minor improvements and code efficiency suggested by Pekka
% v0.9.6- Change to input parameters to tidy up code – “mol” is a vector of
% the HITRAN molecule numbers, to match vector “mix” which is the
% matching mixing ratio. This avoids editing code when different
% GHGs are considered.
% v0.9.7- Added option to turn off line width formula, primarily to show
% what happens when no narrowing of lines due to pressure (or
% temperature changes)
% v0.10- Better tracking of heat movement when model run to equilibrium
% Removed a few unnecessary variables passed back from function,
% added comments for all variables passed back, note – Ts => Tsi in
% this version
% v0.10.1- Added very basic parameterized solar heating by level, otherwise
% other tradeoffs have to be made (e.g.,lose isothermal statrosphere)
% v0.10.2- Now create new GHG conditions, newmix, at time, ntc. Use (via
% separate function “optical_2″ better diffusivity calculation)
% v0.10.3- New layer model, previously the pressure difference for each
% layer was equal, now allow this to change for extending better
% into the stratosphere
% v0.10.4- Further layer model changes – can use one of 6 std atmospheres,
% make GHG concentration now a function of level

% how many molecules are being considered
nmol=length(mol);
% select the molecule names being used from the list in ‘mol’
mols=Hmols(mol,:);
% select isotopologues for each molecule, each molecule, m, will only be
% considered to the depth of ison(m) defined in the function parameters
isoprop=Hisoprop(mol,1:max(ison));

% We don’t know the final size of the vectors and some GHGs have very large
% lists, while others are very small so space considerations mean creating
% a [1 x total lines] array rather than [nmol x max of largest]
for i=1:nmol
if i==1 % first value established the vector

if std_atm>0 % if we want to use one of the 6 AFGL std atmospheres
load(‘AFGL’) % this has height, pressure, temperature and concentrations
% of water vapor, O3, NO2, CH4 – not CO2, set from mixc
numz=find(AFGLp(std_atm,:)<minp,1); % number of boundaries
zb=zeros(1,numz); pb=zeros(1,numz); Tb=zeros(1,numz); rhob=zeros(1,numz);
zb=AFGLz(1:numz); % height of boundaries in meters
pb=AFGLp(std_atm,1:numz); % pressure of boundaries in Pa
Tb=AFGLt(std_atm,1:numz); % temperature of boundaries
Tsi=AFGLt(std_atm,1); % set ocean starting temperature

% calculate pressure, height, temperature, density in midpoint of layer
% [could have used the specified points as the middle of layers but need
% some consistency with the existing calculation of atmospheres]
dz=zeros(1,numz-1);z=zeros(1,numz-1);Tinit=zeros(1,numz-1);
p=zeros(1,numz-1);rho=zeros(1,numz-1);
mix=zeros(nmol,numz-1); % matrix of mixing ratios for molecule, height

p=(pb(2:numz)+pb(1:numz-1))/2; % pressure in midpoint of layer
Tinit=(Tb(2:numz)+Tb(1:numz-1))/2; % initial temperature in midpoint of layer
z=(zb(2:numz)+zb(1:numz-1))/2; % height in midpoint of layer
rho=p./(Tinit*287); % density in midpoint of layer in kg/m^3 from ideal gas law
dz=zb(2:numz)-zb(1:numz-1);

% Create matrix ‘mix'(each molecule,each height)
loch2o=find(mol==1, 1); % location in mol vector of H2O
if ~(isempty(loch2o)) % if water vapor used
% water vapor in mid-point of boundary
mix(loch2o,:)=(AFGLh2o(std_atm,1:numz-1)+AFGLh2o(std_atm,2:numz))/2;
end
locco2=find(mol==2, 1); % location in mol vector of CO2
if ~(isempty(locco2)) % if CO2 is used
mix(locco2,:)=mixc(locco2); % set the concentration of CO2 at all heights the same
end
loco3=find(mol==3, 1); % location in mol vector of O3
if ~(isempty(loco3)) % if O3 is used
% set the concentration of O3 to that atmosphere
mix(loco3,:)=(AFGLo3(std_atm,1:numz-1)+AFGLo3(std_atm,2:numz))/2;
end
locn2o=find(mol==4, 1); % location in mol vector of N2O
if ~(isempty(locn2o)) % if N2O is used
% set the concentration of N2O to that atmosphere
mix(locn2o,:)=(AFGLn2o(std_atm,1:numz-1)+AFGLn2o(std_atm,2:numz))/2;
end
locch4=find(mol==6, 1); % location in mol vector of CH4
if ~(isempty(locch4)) % if CH4 is used
% set the concentration of CH44 to that atmosphere
mix(locch4,:)=(AFGLch4(std_atm,1:numz-1)+AFGLch4(std_atm,2:numz))/2;
end
% are needed ptropo & iztropo needed – not defined with std atmospheres?

% We divide the atmosphere into either approximately equal pressure changes
% or reduce dp for each successive layer by factor ratiodp
% Note: first layer is the boundary layer at pressure = blp
% with different humidity (BLH)

if ratiodp==1 % if equal pressure change for each layer
dp=dptot/(numz-2); % finds the pressure change for each height change

for i=3:numz % locate each value
zi(i)=find(pr<=(pr(zi(2))-(i-2)*dp), 1); % gets the location in the vector where
% pressure is that value
end
elseif ratiodp<1 && ratiodp>0 % if progressively reduced pressure change
dp=(dptot*(1-ratiodp))/(1-ratiodp^(numz-2)); % calculate value of dp to
% meet the total pressure change with this geometric series & ratiodp
ptemp=zeros(1,numz);ptemp(1)=ps;ptemp(2)=blp;
for i=3:numz
ptemp(i)=ptemp(i-1)-dp*ratiodp^(i-3);
zi(i)=find(pr<=ptemp(i),1);
end

% now copy the vector mixc to each height
mix=repmat(mixc’,1,numz);
% now calc mixing ratio in molecules for water vapor at prescribed RH
for i=1:numz-1
% = RH * es / p
% currently satvaph2o() gives es=610.94.*exp(17.625.*Tc./(Tc+243.04)); where Tc in ‘C
if i==1
mix(1,i)=BLH*satvaph2o(Tinit(i))/p(i);
elseif i<iztropo
mix(1,i)=FTH*satvaph2o(Tinit(i))/p(i);
else
mix(1,i)=STaH;
end
end
end

lapsez=ones(1,numz-1)*lapse; % lapse rate for convective adjustment can be
% different for each layer, current all set to ‘lapse’
radupre=zeros(numz,numv);radupost=zeros(numz,numv); % create matrix for some spectral
raddpre=zeros(numz,numv);raddpost=zeros(numz,numv); % values, that only get assigned
% if a jump in GHGs

% ———— Solar Heating of the Atmosphere

sheat=zeros(1,numz-1); % W/m^2 solar heating of each layer of the atmosphere
sheath2o=zeros(1,numz-1); % W/m^2 solar heating by water vapor of each layer of the atmosphere

if solaratm % if we are using solar heating of the atmosphere
% ozone and water vapor, first pass, using solar heating rates from
% Petty (2006), p.313
% first ozone, use 0’C/day at 10km, straight line through 2’C/day at
% 30km and convert to W/m^2
sheat=(z>=10000).*((z./1000-10)/7).*Cp/86400;
%sheat=(z>=10000).*((z./1000-10)/10).*Cp/86400;
% now water vapor – this will be the most badly modeled due to the
% variability of water vapor – convert to W/m^2
sheath2o=(((z<5000).*(z./1000*.6/5+.7))+(z>=5000 & z<9000).*((9000-z)./1000*1.1/4+.3)…
+(z>=9000 & z<15000).*((15000-z)./1000*.3/5)).*Cp/86400;
sheat=sheat+sheath2o;
end

% =========== Finished defining standard atmosphere against height ======
% =========== Absorption and Emission – each layer of the atmosphere, i =====
% note: work in cm for volumes and distances (absorption coefficients are in cm)

T(:,1)=Tinit; % load temperature profile for start of scenario
Ts(1)=Tsi; % ocean temperature for first time step (this could have been the bottom
% layer of the atmospheric model, but feature added after code developed

% first, we timestep
% second, we work through each layer
% third, through each wavenumber
% fourth, through each absorbing molecule
% currently calculating surface radiation absorption up to TOA AND
% downward radiation from TOA (at TOA = 0)

% Change in radiated energy = dI(v) * dv (per second)
% accumulate through each wavenumber
% if the upwards radiation entering the layer is more than
% the upwards radiation leaving the layer, then a heating
dE(i,h)=dE(i,h)+(radu(i,j)-radu(i+1,j))*dv*tstep;
% pre 0.10; Eabs(i)=Eabs(i)+(radu(i,j)-radu(i+1,j))*dv;
end % each wavenumber interval
end % each layer

% Downwards spectral emissive power
for i=numz-1:-1:1 % each layer from the top down
for j=1:numv % each wavenumber interval
% calculate emission downward (same as upward but in case we
% play around later)
emitd(i,j)=abso(i,j)*3.7418e-8.*vt(j)^3/(exp(vt(j)*1.4388/T(i,h-1))-1);
% calculate spectral emissive power at next upward boundary
radd(i,j)=radd(i+1,j)*tran(i,j)+emitd(i,j);
dE(i,h)=dE(i,h)+(radd(i+1,j)-radd(i,j))*dv*tstep;
end % each wavenumber interval
dE(i,h)=dE(i,h)+sheat(i)*tstep; % add in solar absorbed for this layer

if iso_strato && i>iztropo % if we want to keep stratosphere at fixed temperature
% but allow whole stratosphere temp to move in this version,
% defined by the value of temperature in the tropopause layer
T(i,h)=T(iztropo,h-1); % layers above fixed to last iteration of tropopause layer
% calculating from top down so can’t use current time period value

% **** need to account for movement of heat in this part of the
% atmosphere under this scenario – not yet done
else
T(i,h)=T(i,h-1)+dT; % calculate new temperature

if T(i,h)>500 % picks up finite element problem
disp([‘Terminated at h= ‘ num2str(h) ‘, z(i)= ‘ num2str(z(i)) ‘, i = ‘ num2str(i)]);
disp([‘time = ‘ num2str(h*tstep/3600) ‘ hrs; = ‘ num2str(h*tstep/3600/24) ‘ days’]);
disp(datestr(now));
return
end
end
% might need a step to see how close to an equilibrium we are
% getting – not yet implemented because this part of routine is so fast
end % each layer

for j=1:convloop % we need to iterate to a solution
tc=zeros(1,numz-1); % temporary value of convective energy

% Have already worked out the new temperatures for timestep,h,
% based on radiated energy from previous time step, so use this
% timestep and adjust for convection
% Problem is 3 simultaneous equations with 3 unknowns because both
% the temperature above and below change and the lapse rate must be
% valid for the new temperatures

% First from surface to layer 1, slightly different code
if (Ts(h)-T(1,h))/z(1)>lapsez(1) % if lapse rate is too high

end
% Now for all the atmospheric layers
for i=2:numz-1
if (T(i-1,h)-T(i,h))/(z(i)-z(i-1))>lapsez(i) % too cold, convection will readjust
% calculate the convective heat transfer that satisfies the
% simulataneous equation
tc(i)=(T(i-1,h)-T(i,h)-(z(i)-z(i-1))*lapsez(i))/(1/Cp(i)+1/Cp(i-1));

% flux proportions transmitted from each layer & from the surface, by
% wavenumber
stau=zeros(numz-1,numv); % optical thickness by wavenumber summed from top down
TOAf=zeros(numz,numv); % TOA flux proportion, has to include surface so is number of layers + 1
TOAf(numz,:)=dv*emitu(numz-1,:); % get the top layer, no attenuation
stau(numz-1,:)=tau(numz-1,:); % set the first tau value for the top layer

for i=numz-2:-1:1
% flux from layer i+1 transmitted to TOA is emitu(i)*stau(i+1)
% needed to shift everything up 1 layer so that the surface can be
% included
TOAf(i+1,:)=dv*emitu(i,:).*exp(-stau(i+1,:));
% sum tau going downwards
stau(i,:)=stau(i+1,:)+tau(i,:); %
end
TOAf(1,:)=dv*radu(1,:).*exp(-stau(1,:)); % radu(1,:)= surface intensity (below the bottom layer of atmosphere)

% Calculate the lapse rates for each layer at each timestep
alr=zeros(numz-1,nt);
alr(1,:)=1000*(Ts-T(1,:))./z(1); % ocean to layer 1
for i=2:numz-1
alr(i,:)=1000*(T(i-1,:)-T(i,:))./(z(i)-z(i-1)); % value in K/km
end

% Print out inputs to simulation
disp(‘ ‘); % to produce a newline
if std_atm==0 % if user designed profile

na=(Na*rho(i)/mair)/1e6; % number of air molecules per cm^3
d=100*dz(i); % path length in cm

% and so for this slab of depth, d:
% ———– Iterate through each GHG for this slab, m ————–
for m=1:nmol
% fprintf(‘%d %s’, m, ‘ ‘); % display molecule being calculated
if mol(m)==1 % if water vapor, special case
if contabs % if using continuum absorption
% method from Pierrehumbert 2010, p.260, uses svp & units
% in m2/kg

% mixh2o = partial pressure of w.v. / pressure, so
% the partial pressure of water vapor = mixh2o(i) x p(i)
tau(i,:)=tau(i,:)+sig.*pi/100*(mix(m,i)*p(i))^2/(1e4*461.4.*Tlayer(i));
end
end
im=find(molx==mol(m)); % the index for all values of this GHG
immax=length(im); % the number of lines for this GHG
%disp([‘Layer ‘ num2str(i) ‘. Mixing ratio for molecule ‘ num2str(m) ‘ = ‘ num2str(mix(m))]);
% ——- Iterate through each absorption line for this GHG, j —-

for j=1:immax % each absorption line
% im(j) is the index for v, S, etc
% v(im(j))=line center, S(im(j))=line strength, gama(im(j))=half width
% now a code inefficient method – calculate the profile for each
% across the entire wavenumber range

% for this one line we need to find iso # and therefore the relevant
% proportion then x mixing ratio x no of air molecules
nummol=isoprop(m,iso(im(j)))*mix(m,i)*na;
% then calc the actual line width for this temp & pressure
if linewon % normal physics
ga=gama((im(j))).*(p(i)/p0).*(T0/Tlayer(i)).^nair((im(j)));
else % non-real physics for comparison
ga=gama(im(j));
end

% Pekka’s code improvement for v0.9.5
dt1=1./((vt-v(im(j))).^2+ga^2); % line shape across all wavenumbers
tau(i,:)=tau(i,:)+nummol*S(im(j))*ga*dt1; % change to tau for
% this layer and line for all wavenumbers
end % end of each absorption line, j
end % end of each GHG, m

(I decided to test the code with more restricted data where weak lines are left out. As I don’t have the full Hitran database, I replaced Hitran_read by reading the vectors directly from the more limited data set.)

I had to try a few more times before the results were right. I didn’t first realize that the input values for pressure must are given in Pa. That’s mentioned in the list of variables for minp, but not for the ptropo and I missed that for a while. Furthermore your list mentions blz, while the function requires blp, i.e the pressure at the top of the boundary layer rather than the thickness.of boundary layer.

When I got those details right the results are close enough to what you show in your posts.

I may use parts of the code to look at more restricted issues rather than those related to the full atmosphere. Such calculations might be useful in explaining how physics fundamentals determine the properties of the atmosphere.

I’ve updated the code slightly (not published) to allow the stratospheric temperature to change, but still remain isothermal. So the temperature of the layer that contains pressure = pstrato can be adjusted by radiative-convective means, and then the layers above are moved in the next time step to this temperature.

This was changed to check that if we started with a more isothermal troposphere, the model would change the tropospheric temperature over time back to the lapse rate. (It does).

I’m not happy with the treatment of the stratosphere, especially as it has an impact on the troposphere, and one that may be quite significant when CO2 is changed. I recall a paper by Lindzen that put some value on it. I have 18 papers by Lindzen on my PC so not sure about which one, or even if my memory is correct.

Many people argue against the basics because the conceptual model in their heads differs from physics. Or they have grasped 1 point out of 100 in atmospheric physics and don’t understand where it applies and where it doesn’t and what the other 99 say.

Many people claim CO2 “cannot have that effect because” – with theories like water vapor overwhelms CO2, CO2 is already saturated, climate science ignores convection, the 2nd law of thermodynamics… it’s a long list.

Dealing with confusion requires a) demonstrating via theory and experiment and b) trying to help with mental models.

The articles I wrote can be seen under the subheading “Feedback” in the Roadmap section.

In essence I mostly agree with Spencer & Braswell that you can’t really measure climate sensitivity, although their second paper on the subject I am not sure about.

I recommend reading as many of the above papers as you can. A lot of people think that Lindzen first worked out how to do it and then people tried to shoot him down. The reality is that we’ve had over 20 years of different researchers trying to calculate the climate response from decent quality measurements (ERBE onwards).

Not much point reading only one paper and deciding you know the answer (not saying that is you), essential to understand why and how different people got to different answers with the same data sets.

However, and this is just an unfortunate anecdotal history thing, I do find that the “next argument” against “something to do with global warming” is 95% likely to be from someone who hasn’t read a textbook, and has no familiarity with the last 30 years of climate science. In the most recent case, Nikolov & Zeller have been lauded by the internet fraternity for their free-thinking. But honestly, I really don’t want to get past page 5 of their paper, it’s so bad.

Wasn’t the rather thorough trashing of N&Z at WUWT (ignoring comments from the clueless), good enough for you?

that proposes another better approximation. This approximation is based on a rational + logarithm parametrization of the theoretically derived function. They report a 10% increase in computational time, but that applies naturally only to their case.

The paper is short (4 pp.) and contains the derivation of the formulas.

The authors summarize

In summary, the diffusivity factor problem has been reinvestigated by using the bridging function method. A new formula is proposed based on the exact results in the transparent and opaque limits. It is shown that the new scheme is much more accurate than other previously proposed methods. By applied into a one dimensional radiative model, it is found that the new diffusivity factor can produce very accurate results in flux and cooling rate, in the main spectral ranges for CO2, O3, and H2O. Also it is shown that the new proposed diffusivity factor is efficient in computing.

As I wrote in another comment their formula would be an improvement but not the full solution. Its limitations become clear when it’s realized that the outcome depends on the number of layers.

For radiation from the surface to space the formula should be applied to the optical depth of the whole atmosphere. Similarly it should always be applied to the optical depth between altitude of emission and altitude of absorption. This is obviously not possible in the standard approach. The error is most significant at wavelengths with high transmissivity through single layers but much lower for the full atmosphere.

The more correct approach would require having the azimuthal angle as an additional variable. Radiation intensities and transmissivities should be calculated separately for each azimuthal angle. Tau would be calculated for a three-dimensional grid rather than two-dimensional, altitude, frequency, and now in addition azimuthal angle. That’s simple in principle but a lot more to compute.

[…] I created some simulations of different CO2 concentrations using the atmospheric radiation model described (briefly) in Part Two and in detail in Visualizing Atmospheric Radiation – Part Five – The Code. […]

It’s easy to make the code hugely more efficient by taking advantage of the optimized vector handling of Matlab. I replaced the lines

——————–

% ——– Iterate through each wavenumber, k, for this one
% line ———–
dt1=1./((vt-v(im(j))).^2+ga^2);
tau(i,:)=tau(i,:)+nummol*S(im(j))*ga*dt1;
for k=1:numv % each wavenumber
% vt(k) is the wavenumber under consideration
dtau=diff*nummol*d*S(im(j))*ga/pi/((vt(k)-v(im(j)))^2+ga^2);
% note that prior to 0.7.1 the diffusivity factor was not
% used
tau(i,k)=tau(i,k)+dtau;
end

end % end of each absorption line, j
end % end of each GHG, m
end % end of each layer, i

end % end of each absorption line, j
end % end of each GHG, m
tau(i,:)=d/pi*tau(i,:);
tau(i,:)=diff*tau(i,:);
end % end of each layer, i

————————–

In that I got rid of the inner explicit loop that took almost all the time and replaced it by vector operations that Matlab does much more efficiently.

I moved also multiplications by some factors that are constant for each layer to the end to get them out from the innermost loop. That requires a small change in the continuum section (replacing d by pi). The multiplication by diff is now the final operation that applies also to the continuum part. Having it at the end makes it possible to replace the multiplication by a nonlinear correction using the formula of Zhao and Shi. This change would not increase much the computational time as it’s done once for each frequency at every level.

Based on my own experience I add a few comments for those who would like to try the model themselves.

My own experience comes from using an up-to-date version of Matlab that I have based on an university license (while retired I’m still working a little at the university on part-time basis). There’s a free alternative, Octave, that might be usable. I have no experience on that. Even if it works it may be far less efficient, as Matlab does most of its work in highly optimized libraries as long as vectors and matrices are used properly in the code and explicit indexing avoided as far as possible as the example of my above message explains.

So far I have not used directly the HITRAN-database. I got it recently but so far I use data from textfiles that I downloaded from SpectralCalc Line List Browser. I skipped over that section of SoD’s code that reads data from database and replaced it by a few lines of code that read the data directly and without any filtering from the textfile.

In obtaining data from SpectralCalc the spectral range should be set to 200-2500 1/cm, It’s also useful to put a lower limit for line intensity. Reasonable values are H20; 1e-26, CO2: 1e-24, N2O: 1e-22 and CH4: 1e-22. Those four are enough for most. There’s a limit on the amount of data that can be downloaded in one day from one IP-address without subscription. Three downloads are allowed, but I think that the limit is really on the number of molecules rather than downloads. Thus one of those four must be left to a second time.

After getting the data, I opened them in Excel and removed unnecessary columns leaving molecule, isotope, frequency, intensity, air_halfwidth, and t_exponent. Then I joined all molecules in one file and left the column titles only at the top row.

After that I opened the file in Matlab, asked Matlab to create the code to extract the data, and inserted the line

[molx,iso,v,S,gama,nair] = importfile(‘HiTrandata.txt’,2, 93817);

where importfile is the Matlab-generated code and 93817 is the number of lines in my datafile (my above advice leads to less lines). This line replaces the for-loop where the original code reads the database.

SoD might improve the list of variables a little correcting the error that the list contains blz while blp is needed and BLH is misspelled as RLH. It should also be clear that all input pressures are given in Pa.

Input values must be given to the parameters of the main function. Most are just single numbers, but mix is a vector of five elements. I would prefer that SoD tells what are his base values for the parameters as some of mine are not well justified (they are close enough to give essentially the same results).

Then the function can be called listing also the output variables as they are listed in the function declaration.

That’s all when Matlab is available. Octave requires probably some changes in the way the program is invoked.

When I reviewed the changes it looked to me like I made a basic error in 0.9.3 – the continuum absorption does not have diffusivity applied and therefore – using this approximation – the impact on tau for the continuum is understated by a factor of 1.66.

Is this your understanding?

Old code for continuum:
tau(i,:)=tau(i,:)+sig.*d/100*(es*RH)^2/(1e4*461.4.*Tinit(i))

New code for continuum in v0.9.5:
tau(i,:)=tau(i,:)+sig.*pi/100*(es*RH)^2/(1e4*461.4.*Tinit(i));

I also fixed up the minor points noted above (like BLH, blz and explicit on Pa needed in values).

i) when changing which GHG’s are present I can easily cause an error message at the end if I hadn’t carefully updated the printing list – so just had a 5 hr run with dv=0.01 give an error message at the end and produce no values..

ii) updated the flux values to 1 decimal place to more easily compare different runs

————————————-

% Print out GHG concentrations
for i=1:nmol
if mol(i)>1 % i.e. is not water vapor, which is printed via BLH, FTH
fprintf([mols(i,:) ‘= ‘ num2str(mix(i)*1e6) ‘ ppm, ‘],’%s’);
end
end
disp(‘ ‘);

Great code improvement!
Now a 10-layer model with dv=1cm-1 runs in 1:40 instead of about 4+ minutes.

I see that especially my DLR has increased as a result of the continuum correction. This may also be the main reason why the cooling rates in the model were understated in the lower troposphere.
I will rerun some past results and see the impact.

Actually the diffusion factor is 2.0 for the continuum as it’s always when the optical depth is very small. It’s relatively easy to understand both limits. The value 2 for very small optical depth follows from the Lambert’s cosine law in the limit of small optical depth, and the value 1 for very large optical depth from the fact that in that limit the shortest path dominates. The value 1.66 is just a compromise that works reasonably overall.

To get this more correct the Zhao Shi formula could be used. A much simpler parametrization that agrees roughly with their formula would be essentially as good, but the whole approach of using a multiplicative factor can never be made accurate. Calculating correctly the transmission of radiation that passes trough several layers up to the whole atmosphere is requires a different approach.

A calculation where all simplifying assumptions are dropped is computationally very heavy. Fortunately the more efficient approaches are good enough for most purposes.

The runs took 4.5 hours each. That’s 230,000 individual calculations of optical thickness for each of 4 molecules (water vapor, CO2, NO2, CH4) at each of 10 levels. About 10 million calculations in total for optical thickness / transmissivity / emissivity.

The change in TOA and DLR from the Δv = 0.1 cm-1 results in both cases was 0.0 W/m2.

The short comment of 10:39 pm goes a little further on that point. Thinking a little more on that I realized that the extra computational load would perhaps not be too heavy as it’s not necessary to repeat that step that presently takes most of the time. The factor diff would depend on the azimuthal angle but not the step where the line-by-line summation is done for a given layer and wavenumber.

The next step of calculation of tran and abso should be done separately for each azimuthal angle and the multiplication by diff should be done at this stage. tau would remain a two-dimensional matrix while tran and abso would be three-dimensional. All later steps would need to handle the extra dimension. Whenever a sum over wavenumber occurs the sum should be done also over azimuthal angle.

The main point in the above is that that the ability of a photon to pass trough a layer depends both on the wavenumber and on the azimuthal angle. In wave-parallel model both these factors enter the later steps of the calculation in an identical way.

I’ve implemented the Zhao & Shi 2013 method on a new version of the code.

The time impact is zero on models with 10 layers, Δv = 1cm-1 and water vapor, CO2, N2O & CH4.

After trying a few results I see larger differences with high DLR – to be expected as the optical thickness is larger and this is where the errors are highest using the simple diffusivity approximation.

There is no reason not to keep the improvement as it is just 2 more lines of code.

Adding the azimuthal angles has a similar influence on computational time as increasing the number of frequencies included, i.e. using 10 azimuthal angles has the same influence as decreasing dv to one tenth. Therefore I haven’t done many comparisons and I haven’t archived well those runs that I have done.

The main conclusion was that the effect on the final results is not very large, not large enough to change any of the conclusions that you have done using the single diffusivity factor.

A decided to do a new comparison combining my earlier modifications with your more recent code. I have a test running right now, but there may still be some problems in the code. I’ll know only, when the whole code has executed successfully.

Now I have got partial results from one comparison (partial due to a failure when creating plots). The case was close to your part 9, but keeping surface temperature fixed at 288 K. I did extend the calculation up to the 1 hPa level to see, whether changes are larger in stratosphere.

Flux up at TOA changed from 237.0 to 237.7 W/m^2 and down at surface from 282.3 to 281.8 W/m^2. The temperature rose by 0.2 K at tropopause and above. The surface temperature was fixed and the difference grew smoothly going up in the troposphere. I don’t know, why the model allows the lapse rate to be a little higher than 6.5 K/km. The excess was a little less in the model with ten azimuthal segments than with the single diff=1.66.

The differences are small enough to be ignored in most calculations made for educational purposes but it’s certainly nice to know that they are so small.

% how many molecules are being considered
nmol=length(mol);
% select the molecule names being used from the list in ‘mol’
mols=Hmols(mol,:);
% select isotopologues for each molecule, each molecule, m, will only be
% considered to the depth of ison(m) defined in the function parameters
isoprop=Hisoprop(mol,1:max(ison));

1. The atmosphere now includes a very basic parameterized solar heating of the atmosphere. In the model, unlike the real atmosphere, the heating rate does NOT depend on the GHG concentrations. This means it is unrealistic. However, it is more realistic than zero solar atmospheric heating.

2. The isothermal stratosphere can still be turned on, but I’ve turned it off for this and future model runs. Because the model now includes TOA energy accounting, if it is turned on then the TOA energy will never show as balancing because secret solar heating of the stratosphere will be occurring.

3. There is an ocean (isothermal ocean) of depth ‘ocd’, and the model now changes the temperature of the ocean as energy is gained or lost. (So the surface temperature is not constant any longer).

4. Energy changes are tracked in the ocean (‘dEs’) and each atmospheric level (‘dE’) for each time step. Convective energy changes (‘dCE’) between layers are separately logged.

5. TOA and surface flux accounting for each time step are plotted and returned as ‘surfe’ and ‘TOAe’.

6. The code has been better documented at the start with all of the returned variables by the function noted and array sizes documented.

Do you think I have correctly coded the water vapor continuum, as per Pierrehumbert 2010?

I am now questioning my maths.

I had to rethink it when creating the next version of code – where the optical thickness, τ, is created as a separate function, and water vapor mixing ratio is supplied into the function. This will allow us to change optical thickness during the time steps as CO2, for example, is changed and decide whether to keep RH or mixing ratio constant.

So far I have not tried to verify all parts of your code. I have checked with some care those parts that I have modified in some way, but not evrything else. Gradually I’ll go trough the rest as well, and that includes the continuum code.

I had some problems in getting JavaHAWKS to run for filtering parts of HITRAN database in smaller files (on every attempt it stopped having done only fraction of the task). I ended up in coding a simple Windows application that reads in the whole database (all 2.5 million lines) keeping it in memory and allowing filtering based on molecules, number of isotopes frequency and strength of the line. The application tells with unnoticeable delay the number of lines included in selection and the sum of their strengths. Then one can write the selected lines in a file including any number of molecules in the same file. That way it seemed to be easy to create a file of all significant lines keeping the file size reasonable.

I could make the application freely available on my site, if anyone is interested.

More or less the same can be done at SpetralCalc site, but there are some limitations in that. At least it’s much faster to try alternative filters with my application.

I went trough your continuum code and it seems to agree with Pierrehumbert’s book in other ways, but Pierrehumbert tells that the temperature dependence is (296/T)^4.25 while you have the exponent 4.0. The difference is not large as long as the temperature is reasonably close to 296 K. Furthermore the continuum is not important at much lower temperatures, but why not follow Pierrehumbert on this detail as well.

[…] Here is the DLR for 4 different surface temperatures. In each case there is a lapse rate of 6.5 K/km, the boundary layer humidity (BLH) = 100%, the free tropospheric humidity (FTH) = 40% and there were 10 atmospheric layers in the model with a top of atmosphere at 50 hPa. More about the model in Part Two and Part Five – The Code. […]

Re one of your earlier comments, I also realised when trying to handle the stratosphere that the constant Δp was a handicap.

I have updated the code, as yet unpublished, with a variant shown below (not very elegant but gets the job done). This requires the parameter ‘ratiodp’ as an input to the function, which, if =1, will be like the old code, but if less than 1 will reduce the pressure in each subsequent layer by facto ‘ratiodp’. Just a geometric series.

This seems to create a reasonable stratosphere.

My stratosphere, even up to 35km, declines with altitude when I let it relax towards equilibrium. This is why I needed to create a decent one to play around with.

I see two related issues:
1. The stratosphere should increase with temperature
2. The radiative forcing from doubling CO2 is currently too high by about 50%, but we don’t have the correct stratospheric adjustment so can’t determine if that is the reason.

===== new code for creating the pressure profile of the atmosphere ====

if ratiodp==1 % if equal pressure change for each layer
dp=dptot/(numz-2); % finds the pressure change for each height change

for i=3:numz % locate each value
zi(i)=find(pr<=(pr(zi(2))-(i-2)*dp), 1); % gets the location in the vector where
% pressure is that value
end
elseif ratiodp0 % if progressively reduced pressure change
dp=(dptot*(1-ratiodp))/(1-ratiodp^(numz-2)); % calculate value of dp to
% meet the total pressure change with this geometric series & ratiodp
ptemp=zeros(1,numz);ptemp(1)=ps;ptemp(2)=blp;
for i=3:numz
ptemp(i)=ptemp(i-1)-dp*ratiodp^(i-3);
zi(i)=find(pr0 && <=1');
return
end

I have done my calculations of the stratosphere in a similar way using the input value of tropopause pressure as the dividing line between equal pressure difference and equal pressure ratio. Then I had another way of fixing the the share of the layers in stratosphere.

Do you have a value for what the radiative forcing should be for a clear sky atmosphere of a specific surface temperature and simple moisture profile? I have no idea of how much it should deviate from the Myhre calculation that is a global average and takes clouds into account in some way.

I did look a little more at the effect that using azimuthal segments has for the results. As I told already the most important results change little. One place where a difference can be seen is in the calculation of emission from the surface directly to space.

I would expect that the Zhao and Shi formula gives very similar results to the method based on 10 segments in theta. For that it should be used with the optical depth of the atmosphere as input.

We can see that there’s a systematic difference in one direction in the atmospheric window and deviations in the other direction at many discrete lines. The systematic difference in the window region could probably be removed by changing diff to a little larger value than 1.66. In the present case the sum of all frequencies changes from 76.38 to 76.97 W/m^2, when the opposite deviations cancel to a significant extent.

I made one further model run changing the parameter diff from 1.66 to 1.7. I wanted to see, whether that makes the results of the simplified model to agree better with the multi-theta model on the atmospheric transmissivity in the region of atmospheric window. The agreement got, indeed, better in that frequency range but in general the deviation of the simplified model from the more detailed model got larger.

The more important summary quantities differ roughly twice as much from the better model with diff=1.7 as with diff=1.66. Extrapolating from the observations it seems that the best agreement would be obtained for diff in the range 1.62-1.64. We see that the results are rather sensitive to the value of diff. A more extensive study might tell that 1.66 is the optimal compromise or that a slightly lower value would be better. Whatever the answer to that point is, it’s clear that the simple model gives good results only when diff is close to it’s recommended value 1.66.

SoD, great that you provide the code for the radiation calculation. I’ve tried to run the code but failed due missing files of the HITRAN data. e.g. 01_hit08.par and so on. Where and how can I get these data files?

SoD, thanks.
In between I was able to download the missing data fromhttp://www.spectralcalc.com/spectral_browser/db_data.php
I used the range 0 to 10000/cm. I hope this is sufficient.
I have not yet data for all the molecules because the download is limited, but I could run the program successful.

SoD,
so far I have two suggestions.
First increasing the iterations for the convection perhaps to convloop=30; or even higher. This takes not much extra time but is more accurate.

The free tropospheric humidity is also applied to the stratosphere. This results in a too wet and therefore too cool stratosphere, because of too high emissivity. Maybe better would be a constant mixing ratio as at the tropopause.

For the convective adjustment, I expect it to be a simultaneous equation similar to the one given for the 2 layer convective adjustment where the application is for n adjacent layers identified.

Pekka, have you already solved this?

Uli has a good point on the stratospheric water vapor. Without checking it is typically around 6ppm, I will dig out some measurements.

I have a new version of code which produces better graphs when 20 or so layers, includes the diffusivity work of Zhao & Shi, divides up the atmosphere more sensibly considering the stratosphere and allows for changes in GHGs at some point during the time run. I’ll publish that shortly.

I have a solution based on assumption that the convection always starts from the surface. The resulting equation is really simple: some summation and one division. It must be tried for all numbers of layers starting from the surface to figure out the top of convection, but that’s easy again. The code converges with much smaller number of time steps. What’s is still not working is the calculation of energy fluxes within the convective region.

The basic idea is to calculate cumulative sums of heat contents starting from the ocean (inclusive). Then a formula is used to calculate surface temperatures that lead to the similar sums applying convective lapse rate. From the difference of these values the change in surface temperature is calculated that makes the sums equal. That’s repeated for each number of layers calculating also the temperature at the top. This temperature is compared with the radiative solution. Where the calculated temperature and the radiative temperature cross is the tropopause.

The approach might fail it the lapse rate would be exceeded in the middle, but not all the way from surface, but I don’t see any need to calculate such cases.

There may still be some difference in the way I handle the layers, but that would be a minor difference in the way discretization is done and could later be corrected.

It turns out that your original iteration requires hundreds repetitions to get correct results. It’s so fast that putting convloop=500 does not take much time and gives reasonably accurate results. As I tried to explain earlier, the whole iteration can be replaced by a direct calculation of the energy balances and the requirement of exactly right lapse rate. Comparison of convloop=500 with my new code confirms that the two approaches give the same results at the level expected with convloop=500.

I added following initializations before the main timestep iteration loop starts (I start indexing from the surface)

if adjtop>0
dEs(h)=dEs(h)+Cps*(TsT(adjtop+1)-Ts(h));
dESum=-Cps*(TsT(adjtop+1)-Ts(h));
Ts(h)=TsT(adjtop+1);
for i=1:adjtop
dCE(i,h)=dESum;
dESum=dESum-Cp(i)*(Ts(h)+AT(i+1)-T(i,h));
dE(i,h)=dE(i,h)+Cp(i)*(Ts(h)+AT(i+1)-T(i,h));
T(i,h)=Ts(h)+AT(i+1);
end
end

% == end of convective adjustment

I left the original version in making the calculation use my version, when convloop=1 and the old version with higher values. That makes comparisons easier and allows for use of the old version (with a large value of convloop) if the case is pathological in such a way that my assumptions are not valid.

for the lapse rate calculations, I find it important that cases with different convection regions can be done.
Below are my test cases. First the normal case. The convection starts from ground. The other cases are for polar conditions. In the second case the convection is only in the middle. But in the last case there are convection in the middle and at the ground. I think the lapse rate correction should be able to calculate all these cases without wrong results.

One possibility is to just use the original code with a high number of iterations. 500 is not excessive at all.

It would certainly be possible to improve my approach to cover the more complex cases, but it gets significantly more complex. First should the code search for regions of excess in the lapse rate and then apply appropriate code allowing both the lower and the upper level to vary.

I tried what happens when the initial ocean temperature is changed in your third example to be cold already initially. The results of the original model change totally. Thus it seems that the model cannot calculate that case leaving open, whether any case that the model can describe violates the rule that the convective range starts from the surface. I would consider it probable that it cannot.

It’s certainly possible to calculate the radiative fluxes even for complex temperature profiles but the convective corrections may become meaningless in other cases. The model is not a dynamic model of the atmosphere, it’s dynamic behavior does not correspond to that of the atmosphere.

Where the dynamics works better is in calculating warming of the surface, when the rate of change of atmospheric temperatures is so much faster that it can be considered instantaneous. Using a high value for convloop and using my code are two ways of implementing this assumption.

Pekka,
I agree, the final state in the third case depends on the initial state. But I’m not sure if it is an inability of the code, of the model, or it is a real world bistability.

The wavenumber range from vmin=200 to vmax=2500 seems to small. A substantial amount of energy is radiated outside especially at low wavenumbers. I suggest to set vmin=1.
vmax may increased to 3000 but this seems not so important, also because it would slow the code further.

The present model gets its energy from solar radiation. The absorption of solar radiation in the atmosphere is not sufficient for maintaining the adiabatic lapse rate, only the surface can do that as long as we are looking at the Earth atmosphere and not some other strange planet.

There are important physical situations where this is not true, but they are related to subsiding flows and driven by pressure gradients as part of large scale circulation. Modeling a single cell might be a good major next step beyond this model. That would require more from the model in various ways.

Within the range of applicability of the present model it’s probably safe to assume that convection starts always from surface. Inappropriate parameterizations might violate that assumption, and an initial state that’s very far from the stationary solution might prevent reaching the solution. What happens in such cases may depend on the approach to convective adjustments.

I made an update to the code to v0.10.3 – pasted inline and as a link.

1. In this version of the function I have extracted the calculation of optical thickness as optical_2, because I wanted to enable increases in some GHGs at particular times during the run. So as a separate function this can be done at the start and at some particular time (and at multiple times with minor code changes).
This optical_2 code is also linked and appended inline after v0.10.3

2. optical_2 now includes the Zhao & Shi 2013 improvement to the diffusivity approximation for spherical geometry under the plane parallel assumption (more on that in another article).

3. The function now has [newmix] and [ntc] as inputs.
[newmix] is the new version of [mix] to be applied at timestep=ntc.
So for example if mol=[1 2] i.e., H2O & CO2 ; mix=[0 280e-6] and newmix=[0 560e-6] then at the start CO2 is set to 280ppm and at timestep=ntc CO2 is increased to 560ppm. Water vapor of course is set by other parameters (BLH, FTH and a new one..).

4. The function has another new input ‘ratiodp’, which I’ve set to 0.9 in recent runs with 20 layers (numz=21).
This sets a geometric progression to the change in pressure for each layer. So with 0.9 each layer has 0.9 of the pressure drop of the previous layer. This works well for creating a better resolution stratosphere.

At ratiodp=0.8 I found that the convective adjustment wasn’t working well, which needs more investigation.

5. I’ve improved the graphs a little, with a selection of layers to be plotted mostly set by the vector ‘iz’ in the plotting routine.

6. Parameter ‘STaH’ is introduced- stratospheric absolute humidity, as a mixing ratio. This was suggested by Uli and an excellent solution to the taxing problem of why the stratosphere was declining in temperature.
With STaH = 6e-6 (i.e., 6ppmv) the stratosphere increases in temperature:

This is 1300 days at 2 hour timesteps with a change to CO2 from 280ppm to 560ppm at 650 days.

The starting temperature is 288K with a 6.5K lapse rate, a tropopause at 200 hPa, TOA at 1 hPa, 20 layers with each layer only 0.9 of the pressure drop of the previous layer, boundary layer relative humidity at 80%, free tropospheric relative humidity at 40% and stratospheric humidity at a realistic 6ppm.

So, for example, if we run 280ppm CO2 with 288K and conditions…. then we do the same run with 560pmm CO2 with 288k and same conditions… then the difference in TOA flux is NOT radiative forcing.

This is even apart from allowing stratospheric adjustment and considerations of whether cloudy sky RF is significantly different from clear sky RF..

The calculation of RF should be by first finding the equilibrium condition for 280ppm and assessing the TOA outgoing radiation – and THEN changing CO2 to 560ppm and capturing the instantaneous change in TOA radiation.

The wavenumber range from vmin=200 to vmax=2500 seems to small. A substantial amount of energy is radiated outside especially at low wavenumbers. I suggest to set vmin=1.
vmax may increased to 3000 but this seems not so important, also because it would slow the code further.

This is another good call.
I changed vmin to 1 cm-1 and for a run with 20 layers over 400 days in 2 hr steps (4800 time steps) the time increased from 5 mins 25 secs to about 5:50, so I will use this from now on.

[…] Understanding atmospheric radiation is not so simple. But now we have a line by line model of absorption and emission of radiation in the atmosphere we can do some “experiments”. See Part Two and Part Five – The Code. […]

The data can be found in the LBLRTM (publicly available) which is a FORTRAN program. I can’t decipher the useful stuff because I never used FORTRAN.

I did find the data files for these 6 std atmospheres, which, with a bit of text editing and MATLAB reading I turned into a set of MATLAB arrays – see main article (file ‘AFGL’ – “Save as” this file but change from a .doc to a .m file).

They have published a Matlab function for calculating the Voigt lineshape. (Actually their function calculates a complex function of complex argument. The real part of the function is the Voigt function, imaginary part of the argument tells the ratio of the widths.)

It may be worthwhile to limit the use of this function to the high altitudes and to cut off the tails at a reasonable distance from the line.

I have started to implement the Zaghloul and Ali code. At the moment I get badly erroneous results, but it seems already clear that their code can be implemented keeping the model efficient enough when unnecessary computations are avoided. (In this way the outcome may be even faster code than the present one as the present amount of unnecessary computations is very large.)

It has been a bit cumbersome to keep several runs in order and available for comparison. Using structures to include all arguments and all results seems to make that much easier. Thus I changed the function declaration to

The function that calculates τ has a slight update because instead of mixh2o holding the water vapor concentration, mix(molecule, layer) now does that for all GHGs. This allows other GHGs, like ozone, to vary with height (which ozone does in real life).

So now the code has the parameter ‘std_atm’ passed in, if this is ≠0 then a std atmosphere is loaded and various parameters are not used (as noted in the comments to the code). The only GHG which then uses ‘mixc’ is CO2, the other GHGs use the values from one of the six AFGL atmospheres.

CO2 is also prescribed (concentration vs height) in those standard atmospheres but the concentration is mostly constant with height and this is something I want to vary even when starting with a pre-defined atmosphere.

The value of a standard atmosphere is not to calculate changes with GHGs – this will only be relevant if the solar radiation absorption through the atmosphere & the average convection actually makes the temperature profile an equilibrium condition. The value is for comparisons with other calculations, specifically on heating rates, TOA / DLR flux, and spectrum at TOA and surface.

This version doesn’t yet use Pekka’s Voigt (line-shape) improvement or the parameter passing via a data-type “structure”.

These were (in prior versions) just a single value – TOA OLR and surface DLR. Now they are the flux upward and downward at each boundary on the last time step. This allows the plotting of the upward and downward flux through the atmosphere and they are returned by the function. Only a tiny code change for this.

‘Neth’ was also added, which is the net heating rate per layer in ‘C/day.

[…] the Matlab model already created and explained – in brief in Part Two and with the code in Part Five – The Code (note 2). The change in surface emissivity is assumed to be wavelength independent (so if ε = […]

I found it because I went back to review Pekka’s comment about this topic and ran some tests to compare with the Zhao paper.

The code even like this has an error but it’s much closer to their published results. The outstanding error shows up in higher values of τ – the published % error on the diffusivity parameter around τ = 100 is around 10-5%, whereas my calculation shows about 1% error.

This is all still better than the usual diffusivity = 1.66 approximation, but I’d like to resolve it, a necessary step to deal with the question Pekka raised some time ago.

The formula is derived for transmission trough a single layer. In the limit of thin layer it gives the value of 2 for all wavelengths. Therefore it gets worse and worse for the calculation of the whole atmosphere when the layers are made thinner and thinner. That’s a very bad property as we should expect that the accuracy of the calculation is improved when the number of layers is increased.

Being derived for single layer the formula is essentially perfect for the calculation of radiation directly from the surface to the space. In that use no diffusivity factor should be used for the optical depth before the single step where this formula is applied to the full optical depth of the atmosphere.

Spending a lot of computational time, it would be possible to compare both the Zhao & Shi formula and the fixed diffusivity factor to a calculation done with azimuthal segments. That should be done with a variable number of layers. I would not be surprised at all, if the outcome is that the fixed diffusivity factor turns out to be better for everything else except the calculation of transmission trough the whole atmosphere.

I wonder if the code above is the latest version, or have there been enhancements? I am currently going through Pierrehumberts book and this model touches on many aspects of chapter 4, hence my interest. By the way, there is a site called github where code files can be hosted, shared and changes easily tracked. I would be happy to create a repository if there is interest

David,
For my part I can tell that I have made several changes to that code. Some of them are changes or additions to the model, others are made to make the interface better suited for the kind of use I have for the model. The changes are not well finalized or elegant.

Most important additions to the code that I have made are:

1) Addition of the Voigt lineshape that affects the results in the upper stratosphere, but not at lower altitudes

2) Possibility of handling azimuthal angle explicitly. This increases the time of computation quite a lot. It turned out that most results of the model change very little from the one-dimensional approximation used in the original model

3) A different way of handling the lapse rate in calculations, where the temperature profile is allowed to change during computation. As the model is not really capable of determining the profile either with or without my modifications, i do not consider this an essential change although my modifications make it both much faster and more realistic in my view.

4) Putting both arguments and results in structures. This helps greatly saving the results of several different runs for further comparisons.

5) Introducing a cut-off to the lineshape. This is very important, when the CO2 concentration is very high like tens of times higher than in the real Earth atmosphere, but not for any of the foreseen future CO2 concentrations.

It might be worthwhile to work together with SoD to make a better finalized update. As it stands my version is perhaps not in a form suitable for open distribution. I can, however, make the present code available to anyone interested on a more private basis (in a zip-file on web at a address to be supplied by email).

v0.10.5 was fixing something to do with the ‘mix’ variable (no recollection on this without comparing both versions in detail)

v0.10.6 was adding non-black body emissivity to the earth’s surface and reflected DLR – to see what the effect was

Neither produced anything of note, of course the non-blackbody emissivity is (should be) interesting for all the people who are convinced a 0.95 surface emissivity is a “game changer” for understanding the surface and TOA energy balance..

Certainly Pekka has produced more comprehensive versions than me, advanced a few ideas, and has often fixed up my cludgy or erroneous use of Matlab.

I would definitely recommend taking his changes on board. The only reason I haven’t done it is lack of time and then moving onto other topics, like ice ages.

And currently I don’t have a lot of time even for ice ages.

So, in summary I am happy to email you v0.10.5 and v0.10.6 and try and answer any questions on them, but suggest instead taking Pekka’s more advanced versions and building from them.

I decided to put my present model version with all data needed to run it in a zip file. Workspace_1.mat contains model set of parameters.

After opening the workspace the following command starts the calculation of the tropical standard atmosphere

anstropic=hitran_0_10_4p2(argstropic);

Argument sets args1 and args2 lead to the calculation of a non-standard atmosphere based on radiative energy balances and convective corrections calculated by the model. Convection is calculated enforcing a maximum lapse rate given as parameter (the value is 6.5 in these parameter sets). args3 uses the Voigt profile at low pressures, other parameter sets do not use it.

As I noted already, I haven’t documented well all changes that I have done. I try to answer any specific questions people have in trying to use this code.

Hi to everyone, I’m Juan and I’m working with absorption spectroscopy in combustion environments. Anyway, I would like to know if someone could send me a code to read the Hitran database to plot line strength and absorbance against wavenumber for CO,CO2 and H2O.

The subroutine optical_2 of SoD’s radiative transfer model calculates the optical thickness of a layer of air. Emission line strength involves also the Planck emission spectrum calculated by the subroutine planckmv. (In my version of the model the optical_2 is replaced by optical_3p, but the additional features of my version are not of interest in your application.) My zip-file (link in a comment above) contains also files 01_hit08.par, 02_hit08.par, and 05_hit08.par which contain part of the lines H2O, CO2, and CO -lines of the full HITRANo8 database. Very weak lines have been excluded and so have high wavenumbers (over 4000 1/cm), which might be important in combustion environment.

I have coded also a small application that can be used to filter out from the full HITRAN08 database those lines that are of interest. That’s available from

Why not use Hottel/Leckner emissivities? I thought that was the standard method for radiative transfer calculations in a combustion environment. If you still want to work with individual lines, you need the HITEMP and HITEMP supplement database, not the HITRAN database.