Feeds:
Posts
Comments

Archive for January, 2013

The MATLAB program for the previous articles in this series is attached as a .doc file.

It’s not a .doc, it is a text file (but WordPress doesn’t give the option of .txt uploads), so save as, rename to .txt and open in a text editor like Notepad++

[update Jan 30th – see comment]

HITRAN_0_10_4 – replaced

AFGL

optical_3

[update Jan 25th, see comment]

HITRAN_0_10_3

optical_2

[update Jan 20th see comment]

HITRAN_0_10_1

[update Jan 14th – see comment]

HITRAN_0_9_6

[original version published]

HITRAN_0_9_2

calling functions:

Hitran_read_0_2 – reads HITRAN data

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

satvaph2o – Saturation vapor pressure in Pa at temperature T in K, from formula in °C: es=610.94.*exp(17.625.*Tc./(Tc+243.04));

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)

planckmv – planck function for matrix of wavenumber, v, cm-1, and temp, t (K), result in spectral emissive power (pi * intensity), in W/(m2.cm-1)

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).

Water vapor – 01_hit08

CO2 – 02_hit08

O3 – 03_hit08

N2O – 04_hit08

CH4 – 06_hit08

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.

Summary of Code Operation

Matlab specifics

1. The variables on the left of the function start – vt tau stau fluxu fluxd ztropo z dz p rho mixh2o zb pb Tb rhob Tinit T radu radd rads radu1 emitu TOAf TOAtr – are variables being brought back. Some are important, some are for checking, some are curiosity value.

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.

Code function

2. Read HITRAN database for the relevant wavenumbers and molecules – to get line strength, position, air-broadening and temperature exponent (line width changes with pressure & temperature).

3. a) i) Create a high res atmosphere from surface temperature, lapse rate and pressure of tropopause – with isothermal stratosphere. Pressure is calculated from ideal gas law.

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

OR

c)  ii) Use one of six standard atmospheres (AFGL) – more info in Part Twelve – Heating Rates

4. Calculate optical thickness for each layer at each wavenumber. dv defines the spacing for the calculation.

For each layer, get each absorption line of each GHG:

– then calculate the actual strength from the line shape equation at each wavenumber of interest across the whole band

– calculate delta tau from that line at that wavenumber and add it to tau(layer,wavenumber)

5. Initialize surface upwards radiation from temperature of surface & emissivity, downward from TOA to be zero.

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

% ====== Input conditions (passed to function) ====================
% dv = resolution of wavenumber for calculations in cm^-1
% std_atm = 1-6 to use tropical; mid-latitude summer/winter; subarctic
% summer/winter; US std 1976; *OR* =0 if we define our own from params:
% numz, ratiodp, blp, BLH, FTH, STaH, Tsi, lapseikm, ptropo
% i.e., the above list is NOT used if std_atm>0
% numz = number of boundaries to consider (number of layers = numz-1)
% ratiodp = scaling for the dp across each layer, <=1 (1=each layer the
% same, 0.9 means each layer 10% less pressure than the one before)
% mol = vector of molecule # to be considered from HITRAN numbers: 1=H2O,
% 2=CO2, 3=O3, 4=N2O, 5=CO, 6=CH4, must match mixing ration from
% vector ‘mix’, e.g., mol=[1 2 4 6]; when mix=[0 360e-6 319e-9 1775e-9]
% which means H2O (ignore because set by BLH/FTH), CO2 = 360ppm, N2O
% = 319 ppbv, CH4 = 1775 ppbv
% mixc = mixing ratio for each molecule as described for ‘mol’, but note
% if first value is for H2O it is unused (BLH & FTH are used)
% newmix = new mixing ratio to be introduced (as step function for now) at
% timestep ntc
% ison = number of isotopes for each GHG, must match ‘mol’: 1 = only the
% main is used), e.g., [1 1 1 1] means first isotopologue for the
% 4 GHGs being; if mol=[1 2 3 4] & ison=[1 3 1 1] then CO2 will be
% considered to the 3rd isotopologue but the others only to first
% blp = boundary layer thickness in Pa
% BLH = relative humidity in the boundary layer
% FTH = relative humidity in the free troposphere (currently layers 2,3..)
% STaH = mixing ratio of water vapor in stratosphere
% contabs = true or false for including the water vapor continuum absorption
% linewon = whether to keep linewidth adjustment on or off (for real
% physics = true, for fake physics = false)
% Tsi = initial surface temperature
% lapsekm = lapse rate in troposphere, K/km
% lapseikm = initial lapse rate in troposphere, K/km, to setup initial temp
% profile
% ptropo = defined tropopause pressure
% minp = top of atmosphere to consider for radiative calcs, in pressure (Pa)
% ocd = ocean depth in m – for running dynamic simulations
% src = average solar absorbed radiation (over surface area of earth)
% tstep = timestep in seconds
% nt = number of timesteps
% ntc = number of timesteps where ‘newmix’ of GHGs is introduced, so ntc < nt
% or, if newmix is not wanted, make ntc>nt and newmix will not be
% invoked
% ploton = whether the plot is run or not – lots of other plots are
% generated via separate programs

% ====== Output variables (returned from function) ====================
% vt(numv) = wavenumber vector
% tau(numz-1,numv) = optical thickness by layer and by wavenumber
% fluxu = upward flux at each boundary
% fluxd = downward flux at each boundary
% neth = net radiative heating rate ‘C/day, each layer
% z(numz-1) = height of each layer
% dz(numz-1) = thickness of each layer
% p(numz-1) = pressure of each layer
% rho(numz-1) = density of each layer
% mix(nmol,numz-1) = mixing ratio of each GHG at each layer
% zb(numz) = height at each boundary
% pb(numz) = pressure at each boundary
% Tb(numz) = temperature at each boundary
% Tinit(numz-1) = starting temperature of each layer
% T(numz-1,nt) = temperature of each atmosphere layer at each time step, T(1,:)=Tinit
% Ts(nt) = temperature of ocean
% dE(numz-1,nt) = change in energy each time step in each layer
% dEs(nt) = change in energy each time step in ocean
% dCE(numz-1,nt) = change in convective heat each time step upward to each
% layer (1=from ocean to atmosphere layer 1)
% radu(numz,numv) = upward flux at each boundary, boundary 1 = surface
% radd(numz-1,numv) = downward flux at each boundary
% emitu(numz-1,numv) = emitted radiation from within each layer, layer 1= first
% atmospheric layer
% sheat(numz-1) = absorbed solar in each layer in W/m^2
% surfe(nt) = net flux at surface each timestep (+ve downward)
% TOAe(nt) = net flux at TOA each timestep (+ve downward)
% rf = radiative forcing – from the TOA flux pre and post change in GHG, =0
% if this condition doesn’t take place (GHG change)
% radupre, raddpre, radupost, raddpost – spectrum up and down pre and post
% GHG change
% TOAf(numz,numv) = TOA flux contributed from each boundary
% TOAtr(numv) = surface radiation transmitted to TOA
% alr(numz-1,nt) = actual lapse rate between each layer at each time step, in K/km

disp(‘ ‘);
disp([‘— New Run v0.10.4 — ‘ datestr(now) ‘ —-‘]);
disp(‘ ‘);

% =========== Basic definitions =====================

% Units in SI unless otherwise noted including all pressure values in Pa
% except for SI, frequency in cm^-1
p0=1.013e5; % std pressure for the HITRANS database in Pa
T0=296; % std temperature for HITRANS database in K
Na=6.02e23; % Avogradro’s number
mair=28.96e-3; % molar mass of air
% mh2o=18.02e-3; % molar mass of water vapor
ems=1.0; % emissivity of surface
cp=1005; % cp (specific heat cap at constant pressure) of air
cpw=4180; % cp of air, J/kg.K
rhow=1000; % density of water, kg/m^3
solaratm=true; % solar heating of the atmosphere, currently very parameterized
convadj=true; % if true, convective adjustment to lapse rate = lapse
convloop=10; % no of times to repeat convective loop to get lapse rate correct
iso_strato=false; % if true, fixed and isothermal stratospheric temperature
% if solaratm is true, this should be set to false
lapse=lapsekm/1000; % lapse rate in K/m
lapsei=lapseikm/1000; % initial lapse rate in K/m
rf=0; %radiative forcing if GHG change
print_layers=false; % if true, prints out z, p, T, rho for each layer
windmax=10000/8; windmin=10000/12; % ‘atmospheric window’ region in wavenumber
% for “curiosity value calculation
% set up a few print strings
if contabs
conttext=’Continuum included’;
else
conttext=’No Continuum’;
end
if convadj
convtext=’With convective adjustment’;
else
convtext=’No convective adjustment’;
end
if linewon
linetext=’ Linewidth depends on p,T; ‘;
else
linetext=’ Constant linewidth; ‘;
end
dbh=false; % debug heat stuff added in 0.9.8, if true print a lot out

% =========== Define wavenumber ranges ==============

vmin=10; % min wavenumber to consider (overall)
vmax=2500; % max wavenumber to consider (overall)
% because absorption lines below the lower boundary (and above the upper
% boundary) can still affect the wavenumbers under consideration we need to
% extend the boundaries of which absorption lines to consider
% rvmin=1; % ratio below lower boundary to consider – not useful for a wide range
% rvmax=1; % ratio above upper boundary to consider

vt=vmin:dv:vmax; % vt = wavenumbers we will calculate optical depth for
numv=length(vt); % number of wavenumbers

% =========== Read in GHG data ======================

fprintf(‘%s’, ‘Loading HITRAN data.. ‘);

% names for the HITRAN molecules up to molecule 6, need to be 3 characters
% each because of how Matlab handles strings
Hmols=[‘H2O’; ‘CO2’; ‘O3 ‘; ‘N2O’; ‘CO ‘; ‘CH4’;];

% atmospheric proportions of HITRAN isotopes up to molecule 6, data from
% Rothman et al (2005) table 6, vector ‘ison’ extracts the values needed
Hisoprop=[0.9973 0.0020 0.0004; 0.9842 0.0111 0.0039; 0.9929 0.0040 0.0020; …
0.9903 0.0036 0.0036; 0.9865 0.0111 0.0020; 0.9883 0.0111 0.0006];

% select data on GHGs from vector “mol”, isotopes to consider, define isotope ratios

% 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));

% Now read in the line data via ‘Hitran_read’ function

% v=line position, S=line strength, iso=isotopologue, gama = air-broadened
% half-width, nair=temperature-dependence exponent for gama

% 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

[ v S iso gama nair ] = Hitran_read_0_2(mol(i), vmin, vmax, ison(i));
% now set an identifier
molx=ones(1,length(v))*mol(i); % identifies all these values as this molecule

else % now subsequent values get appended
[ v1 S1 iso1 gama1 nair1 ] = Hitran_read_0_2(mol(i), vmin, vmax, ison(i));

molx1=ones(1,length(v1))*mol(i); % identifies all these values as this molecule
v=[v v1]; iso=[iso iso1]; gama=[gama gama1]; S=[S S1]; nair=[nair nair1];
molx=[molx molx1];
end
end

fprintf(‘%s’, ‘finished.. ‘);
% =========== Finished reading GHG data ======================

% =========== Define standard atmosphere against height =============

ps=1.013e5; % define surface pressure

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?

else % we want to create our own atmosphere based on the surface temperature,
% Tsi, the lapse rate, lapsekm, boundary layer pressure, blp

% —- First create a “high resolution” atmosphere
% zr = height, pr = pressure, Tr = temperature, rhor = density

maxzr=50e3; % height of high resolution atmosphere
numzr=5001; % number of points used to define high res. atmosphere
zr=linspace(0,maxzr,numzr); % height vector from sea level to maxzr

% use ‘define_atmos’ function to determine pr, Tr, rhor & ztropo (height of
% tropopause)
[pr Tr rhor ztropo] = define_atmos_0_5(zr,Tsi,ps,ptropo,lapsei);
% for consideration – if the radiative-convective model changes temperature
% significantly, should the pressure be re-adjusted?

% —- Second, create “coarser resolution” atmosphere – this reduces
% computation requirements for absorption & emission of radiation
% – this atmosphere has numz boundaries and numz-1 layers

% z, p, Tinit, rho; subset of values used for RTE calcs

% 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)

zi=zeros(1,numz); % zi = lookup vector to “select” heights, pressures etc
zi(1)=1; % make first value the surface value
zi(2)=find(pr<=blp, 1); % gets the location of boundary layer top edge
dptot=blp-minp;

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

else % bad parameter
disp(‘Error, ratiodp must be >0 && <=1′);
return
end

% now create the vectors of coarser resolution atmosphere
% zb is height of boundaries, zb(1) = surface; zb(numz) = TOA
% T, p, rho are all the midpoint (pressure) between the boundaries

zb=zr(zi); % height at boundaries
pb=pr(zi); % pressure at boundaries
Tb=Tr(zi); % starting temperature at boundaries
rhob=rhor(zi); % density at boundaries

% now calculate density, pressure, temperature and water vapor saturation
% mixing ratio at midpoint of each layer – midpoint by pressure
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);
im=zeros(1,numz-1); % index for getting midpoints at midpoint pressure

% calculate pressure in midpoint of layer
for i=1:numz-1
p(i)=(pb(i+1)+pb(i))/2; % pressure in midpoint of layer
dz(i)=zb(i+1)-zb(i); % thickness of each layer
end

% lookup midpoint values for z, T, rho at the midpoint pressure
for i=1:numz-1
im(i)=find(pr<=p(i), 1); % gets the location in the vector of midpoint pressure
end

z=zr(im); % calculate height of midpoint pressure
Tinit=Tr(im); % temperature of midpoint pressure
rho=rhor(im); % density of midpoint pressure
iztropo=find(pb<ptropo,1)-1; % location in the p vector of the first layer above the tropopause

% 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

% calculate Heat capacity for each layer, Cp = rho x dz x cp in J/K
Cp=zeros(1,numz-1);
Cp=rho.*dz.*cp; % of each atmospheric layer
Cps=rhow*ocd*cpw; % of ocean

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

if print_layers % display values at each pressure
disp(‘ ‘);
if std_atm==0
disp([‘dp= ‘ num2str(dp)]);
disp([‘Pressure of tropopause layer = ‘ num2str(z(iztropo)) ‘, at ‘ num2str(iztropo)]);
end
for i=1:numz-1
disp([‘p= ‘ num2str(p(i)/100,’%4.0f’) ‘, z= ‘ num2str(z(i),’%4.0f’)…
‘, T= ‘ num2str(Tinit(i)) ‘, rho= ‘ num2str(rho(i),’%4.2f’)]);
end
% for i=1:numz
% disp([‘pb= ‘ num2str(pb(i)/100,’%4.0f’) ‘, z= ‘ num2str(zb(i),’%4.0f’)]);
% end
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)

% call optical thickness for each layer of the atmosphere
% this function can be repeated if GHGs change (or temperature changes)
% have to pass the atmospheric layers (p, T, rho, mix, numz-1), the
% HITRAN data (v, S, iso, gama, nair) and the GHG concentrations (mix,
% ison)

tau=optical_3(vt, v, S, iso, gama, nair, molx, numz, numv, p, p0, Tinit, …
T0, Na, rho, mair, dz, mol, nmol, mix, isoprop, contabs, linewon);
% now calculate the transmissivity and absorptivity for each layer and
% wavenumber
% i = layer, k = wavenumber
tran=zeros(numz-1,numv); abso=zeros(numz-1,numv); % set array
for i=1:numz-1 % each layer
for k=1:numv % each wavenumber interval
tran(i,k)=exp(-tau(i,k)); % transmissivity = 0 – 1
abso(i,k)=(1-tran(i,k)); % absorptivity = emissivity = 1-transmissivity
end
end

% =========== Calculate absorption, emission & temperature changes =========

Ts=zeros(1,nt); % define array to store ocean temperature
T=zeros(numz-1,nt); % define array to store T for each atmospheric level and time step
dEs=zeros(1,nt); % energy change for ocean for each time step
dE=zeros(numz-1,nt); % energy change for each atmospheric level
dCE=zeros(numz-1,nt); % convective energy tracking

sr=ones(1,nt)*src; % solar flux at TOA vs time (later can vary with time)
sr0=sr-sum(sheat); % define solar flux absorbed at surface (=solar at
% TOA – atmospheric absorbed)

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

TOAe=zeros(1,nt); % TOA instantaneous flux balance vs time
surfe=zeros(1,nt); % Surface instantaneous flux balance vs time

% =========== Main loops to calculate TOA spectrum & flux =====

% 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)

% h = timestep, i = layer, j = wavenumber
for h=2:nt % main iterations to achieve equilibrium
if dbh
disp([‘Timestep ‘ num2str(h)]);
%fprintf(‘%d %s’, h, ‘ ‘); % display current timestep being calculated
end

if h==ntc % change to GHG concentrations, currently keep water vapor
% concentration the same and use same temperature as
% starting conditions

tau=optical_3(vt, v, S, iso, gama, nair, molx, numz, numv, p, p0, Tinit,…
T0, Na, rho, mair, dz, mol, nmol, newmix, isoprop, contabs, linewon);

% now calculate the transmissivity and absorptivity for each layer and
% wavenumber
% i = layer, k = wavenumber
tran=zeros(numz-1,numv); abso=zeros(numz-1,numv); % set array
for i=1:numz-1 % each layer
for k=1:numv % each wavenumber interval
tran(i,k)=exp(-tau(i,k)); % transmissivity = 0 – 1
abso(i,k)=(1-tran(i,k)); % absorptivity = emissivity = 1-transmissivity
end
end
end

radu=zeros(numz,numv); % initialize upward intensity at each *boundary* and wavenumber
radu1=zeros(numz,numv); % initialize curiosity version of radu (no emission)
emitu=zeros(numz-1,numv); % initialize emission upward for each *layer* and wavenumber

radd=zeros(numz,numv); % initialize downward intensity at each *boundary* and wavenumber
emitd=zeros(numz-1,numv); % initialize emission downward for each *layer* and wavenumber
radu(1,:)=ems.*planckmv(vt,Ts(h-1)); % upward *surface* radiation vs wavenumber
radu1(1,:)=radu(1,:); % evaluating radiation transfer with no emission – for interest
radd(end,:)=zeros(1,numv); % downward radiation at TOA vs wavenumber (=zero)
% units of radu, radd, emitu, emitd are W/m^2.cm^-1, i.e., spectral emissive power

% Upwards spectral emissive power
for i=1:numz-1 % each layer
for j=1:numv % each wavenumber interval

% calc intensity (each wavenumber) at next boundary
% = incoming * transmissivity + emission
% emission upward = Planck formula x emissivity (=absorptivity)
emitu(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
radu(i+1,j)=radu(i,j)*tran(i,j)+emitu(i,j);
radu1(i+1,j)=radu1(i,j)*tran(i,j); % transmitted only, for curiosity

% 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

dT=dE(i,h)/Cp(i); % change in temperature = dQ/heat capacity
if dbh
disp([‘Layer ‘ num2str(i) ‘, dE = ‘ num2str(dE(i,h)) ‘, dT= ‘ num2str(dT)]);
end

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

% now ocean heat changes (via radiation) = time x (Solar + Rdown -Rup)
dEs(h)=(sr0(h)+(sum(radd(1,:))-sum(radu(1,:)))*dv)*tstep;
% Temperature change = Energy change / Cp
Ts(h)=Ts(h-1)+ dEs(h)/Cps;
if dbh
disp([‘Ocean dE = ‘ num2str(dEs(h)) ‘, new Ts= ‘ num2str(Ts(h))]);
end

% now apply convective adjustment
if convadj==true % if convective adjustment chosen..

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

tc(1)=(Ts(h)-T(1,h)-z(1)*lapsez(1))/(1/Cps+1/Cp(1));
if tc(1)<0 % should always be positive ***
disp([‘Error in convective adj. layer 1, h= ‘ num2str(h)…
‘, i= ‘ num2str(i) ‘, dCE= ‘ num2str(tc(1))]);
end

if dbh
disp([‘Lapse to layer 1= ‘ num2str((Ts(h)-T(1,h))/z(1)*1000) ‘K/km, newT= ‘…
num2str(newT) ‘, old atm T= ‘ num2str(T(1,h)) ‘, T(ocean)= ‘ num2str(Ts(h))]);
end

dEs(h)=dEs(h)-tc(1); % reduce heat tracking stored in ocean
dE(1,h)=dE(1,h)+tc(1); % increase heat tracking stored in atmos layer 1
T(1,h)=T(1,h)+tc(1)/Cp(1); % calc revised atmospheric temperature
Ts(h)=Ts(h)-tc(1)/Cps; % calc revised ocean temperature

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));

if tc(i)<0 % should always be positive ***
disp([‘Error in atm convective adj, h= ‘ num2str(h)…
‘, i= ‘ num2str(i) ‘, dCE= ‘ num2str(tc(i))]);
end
if dbh
disp([‘Lapse from layer ‘ num2str(i) ‘ to ‘ num2str(i-1)…
‘ = ‘ num2str((T(i-1,h)-T(i,h))/(z(i)-z(i-1))*1000) ‘ K/km, newT= ‘ …
num2str(newT) ‘, old atm T= ‘ num2str(T(i,h))]);
end

dE(i-1,h)=dE(i-1,h)-tc(i); % reduce heat tracking stored in layer below
dE(i,h)=dE(i,h)+tc(i); % increase heat tracking stored in this layer
T(i,h)=T(i,h)+tc(i)/Cp(i); % set revised atmospheric temperature
% also adjust atmospheric layer below (has lost heat)
T(i-1,h)=T(i-1,h)-tc(i)/Cp(i-1);

end
end
% now accumulate the convective adjustments
dCE(:,h)=dCE(:,h)+tc’;
end
end

% instantaneous imbalance at TOA (+ve downward) in W/m^2 = (solar rate – olr rate)
TOAe(h)=sr(h)-sum(radu(numz,:))*dv;
% instantaneous imbalance at surface in W/m^2 +ve downward = solar + DLR – Surface
% emitted – convected (per second)
surfe(h)=sr0(h)+sum((radd(1,:))-radu(1,:))*dv-dCE(1,h)/tstep;

% 1st law of thermodynamics check; heat gained = Es +sum(E)
dqgain=dEs(h)+sum(dE(:,h)); % net energy absorbed in the time step
if abs(TOAe(h)*tstep-dqgain)>1 % if more than 1W discrepancy
disp([‘Time step ‘ num2str(h) ‘, Net Energy missing = ‘ num2str(TOAe(h)*tstep-dqgain)]);
end

% if 1st or middle timestep, or just before & after change in GHGs..
if h==2
% print out initial results of flux
disp([‘ Time= ‘ num2str(h*tstep/86400,’%4.2f’) ‘ days, Flux up at TOA = ‘…
num2str(sum(radu(end,:))*dv,’%4.1f’) ‘ W/m^2, Flux down at surface = ‘ …
num2str(sum(radd(1,:))*dv,’%4.1f’) ‘ W/m^2’]);
elseif h==ntc-1
% display values
disp([‘ Just before GHG change: time= ‘ num2str(h*tstep/86400,’%4.2f’)…
‘ days, Flux up at TOA = ‘ num2str(sum(radu(end,:))*dv,’%4.1f’)…
‘ W/m^2, Flux down at surface = ‘ num2str(sum(radd(1,:))*dv,’%4.1f’) ‘ W/m^2’]);
% capture upwards and downwards spectrum
raddpre=radd;radupre=radu;
elseif h==ntc
disp([‘ Just after GHG change: time= ‘ num2str(h*tstep/86400,’%4.2f’)…
‘ days, Flux up at TOA = ‘ num2str(sum(radu(end,:))*dv,’%4.1f’)…
‘ W/m^2, Flux down at surface = ‘ num2str(sum(radd(1,:))*dv,’%4.1f’) ‘ W/m^2’]);
rf=TOAe(ntc-1)-TOAe(h); % radiative forcing from the change
disp([‘Pseudo radiative forcing= ‘ num2str(rf,’%4.1f’) ‘ W/m^2’]);
disp(‘ ‘);
% capture upwards and downwards spectrum
raddpost=radd;radupost=radu;
end
end % —- timestep iterations

% =========== Finish loops to calculate TOA spectrum, flux =====

% Calculate TOA and DLR flux
fluxu=sum(radu’)*dv; % calculate the final up flux at each boundary
fluxd=sum(radd’)*dv; % calculate the final down flux at each boundary

% Calculate heating rates
neth=86400*(diff(fluxd)-diff(fluxu))./Cp; % +ve number = heating, ‘C/day

% =========== Calculate some curiosity values ===============

% flux transmitted from surface and in “atmospheric window”
% calculate the surface radiation transmitted to TOA (integral with dv)
TOAtr=dv*radu(1,:)*exp(-sum(tau))’;
iw=find(vt>=windmin & vt<=windmax); % find the window region in the vector

% 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

disp([‘Tsi= ‘ num2str(Tsi) ‘, Lapse rate= ‘ num2str(lapsekm) ‘, ‘ conttext ‘, BLH= ‘ num2str(BLH*100)…
‘%, FTH= ‘ num2str(FTH*100) ‘%, blp = ‘ num2str(blp/100) ‘hPa, TOA = ‘ num2str(minp/100) ‘hPa’]);
disp([num2str(numz-1) ‘ layers, dv= ‘ num2str(dv) ‘cm-1, ‘ linetext num2str(nt) …
‘ timesteps of ‘ num2str(tstep/3600/24,’%4.2f’) ‘ days’]);
disp([‘Solar absorbed= ‘ num2str(src) ‘ W/m^2, Absorbed in atmosphere= ‘ …
num2str(sum(sheat),’%4.2f’) ‘ W/m^2’]);
disp(convtext);

else % using a std atmosphere
disp([‘Using AFGL ‘ AFGLname(std_atm,:)]);
disp([‘Surface temperature= ‘ num2str(Tsi) ‘ K, TOA = ‘ num2str(minp/100) ‘hPa’]);
disp([num2str(numz-1) ‘ layers, dv= ‘ num2str(dv) ‘cm-1, ‘ linetext num2str(nt) …
‘ timesteps of ‘ num2str(tstep/3600/24,’%4.2f’) ‘ days’]);
disp([‘Solar absorbed= ‘ num2str(src) ‘ W/m^2, Absorbed in atmosphere= ‘ …
num2str(sum(sheat),’%4.2f’) ‘ W/m^2’]);
disp([convtext ‘, ‘ conttext]);
end

% 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(‘ ‘);

% Print out various flux values
disp(‘ ‘);
disp([‘ Final Time= ‘ num2str(h*tstep/86400,’%4.2f’) ‘ days’]);
disp([‘Flux up at TOA = ‘ num2str(fluxu(numz),’%4.1f’) ‘ W/m^2, Flux down at surface = ‘ …
num2str(fluxd(1),’%4.1f’) ‘ W/m^2’]);
disp([‘TOA transmitted surface flux= ‘ num2str(TOAtr,’%4.1f’) ‘ W/m^2, Surface emission= ‘ …
num2str(sum(radu(1,:)),’%4.1f’) ‘ W/m^2, Window transmission= ‘ …
num2str(sum(radu1(numz,iw))/sum(radu(1,iw))*100,’%4.1f’) ‘ %’]);
disp([‘Starting flux imbalance TOA= ‘ num2str(TOAe(2),’%4.2f’) ‘ W/m^2, Final imbalance TOA= ‘ …
num2str(TOAe(end),’%4.2f’) ‘ W/m^2’]);

% Print out temperature values & lapse rate
disp(‘ ‘);
disp(‘Temperature change, start to finish ‘);
for i=numz-1:-1:1
disp([‘Layer ‘ num2str(i) ‘ = ‘ num2str(T(i,nt)-T(i,1),’%4.2f’)]);
end
disp([‘Temp change ocean = ‘ num2str(Ts(nt)-Ts(1),’%4.2f’)]);
disp(‘Final lapse rate ‘);
for i=numz-1:-1:2
disp([‘Lapse rate ‘ num2str(i-1) ‘-‘ num2str(i) ‘ = ‘ …
num2str(alr(i,end),’%4.2f’) ‘ K/km’]);
end
disp([‘Lapse rate ocean – 1 = ‘ num2str(alr(1,end),’%4.2f’) ‘K/km’]);

% =========== Plots, if needed ==================
if ploton

% ====== Graphs vs time – Temperature different layers, TOA & surface
% forcing energy change each layer, convective energy
figure(1)
tv=2:nt; % time vector, a lot of values are first calculated on timestep 2

iz=[1 4 8 12 16 20]; % atmospheric layers to consider (to reduce graph clutter)
nump=length(iz); % number of items plotted

% Create relevant graph legends for the selected layers
textleg(1)={‘Ocean’};
for i=1:nump
textleg(i+1)={[num2str(iz(i))]};
end
nump=length(iz); % number of items plotted
textleg1(1)={[‘Ocean->’ num2str(iz(1))]};
for i=2:nump
textleg1(i)={[num2str(iz(i)-1) ‘-‘ num2str(iz(i))]};
end

% 1. Temperature anomaly vs time
subplot(2,2,1),plot((1:nt)*tstep/86400, Ts-Ts(1), …
(1:nt)*tstep/86400, T(iz,:)-(repmat(T(iz,1),1,nt)))
ylabel(‘Temperature Anomaly, K’)
xlabel(‘Time’)
title([‘Temperature anomaly, initial Ts=’ num2str(Tsi),’K, Lapse=’ num2str(lapsekm)])
grid on
%legend(‘Ocean’,’1′,’2′,’3′,’4′,’5′,’6′,’7′,’8′,’9′, ’10’)
legend(textleg);

% 2. TOA & surface imbalance vs time
subplot(2,2,2),plot(tv*tstep/86400,TOAe(tv),tv*tstep/86400,surfe(tv))
xlabel(‘Time, days’)
ylabel(‘Flux imbalance, W/m^2’)
grid on
legend(‘TOA’,’Surface’)

% 3. Energy change each layer vs time (converted to instantaneous)
subplot(2,2,3),plot(tv*tstep/86400,dEs(tv)/tstep,tv*tstep/86400,dE(iz,tv)/tstep)
xlabel(‘Time, days’)
ylabel(‘Instantaneous Energy Gain, W/m^2’)
title([‘Each time step= ‘ num2str(tstep/86400,’%3.2f’) ‘ days’])
grid on
legend(textleg);

% 4. Convective flux between layers vs time
subplot(2,2,4),plot(tv*tstep/86400,dCE(iz,tv)/tstep)
xlabel(‘Time, days’)
ylabel(‘W/m^2’)
grid on
title([‘Convective flux between Layers, Each time step= ‘ num2str(tstep/86400,’%3.2f’) ‘ days’])
legend(textleg1);

% ====== Cumulative transmitted flux to TOA & Emitted flux each layer
figure(2)
% Cumulative transmitted flux to TOA
zz=[0 z]; % add the surface to the height vector to make it possible to plot
zTOA=sum(TOAf’).*dv; % get the integrated value across wavelengths, each boundary
accTOA=zeros(1,numz); % accumulated values
accTOA(1)=zTOA(1);
for i=2:numz
accTOA(i)=accTOA(i-1)+zTOA(i); % calculated accumulated values (as per Part 3 blog comment)
end
subplot(1,2,1),plot(accTOA,zz/1000,’r*’)
xlabel(‘Cumulative flux reaching TOA, W/m^2’)
ylabel(‘Height of layer, km’)
title(‘Cumulative transmitted flux from each layer to TOA’)
grid on
% Emitted flux from each atmospheric layer & surface
% ’emitu’ is from each atmos. layer, we need to add surface
emitted=[sum(radu(1,:))*dv sum(emitu’)*dv];
subplot(1,2,2),plot(emitted,zz/1000,’b*’)
xlabel(‘Emitted flux, W/m^2’)
ylabel(‘Height of layer, km’)
title(‘Emitted flux from each layer’)
grid on

% ====== Plot lapse rate as a check, and temperature profile vs height
figure(3)
% Actual lapse rate – for checking
% Not so useful when past 10 layers, so just the first 9
subplot(1,2,1),plot((1:nt)*tstep/86400,alr(1:10,:))
xlabel(‘Time, days’)
ylabel(‘Lapse rate, K/km’)
title(‘Lapse rate change with time for bottom set of layers’)
grid on
legend(‘Ocean->Layer 1’, ‘1-2′,’2-3′,’3-4′,’4-5′,’5-6’, ‘6-7′,’7-8′,’8-9′,’9-10’)

% Vertical temperature profile at 3 times
subplot(1,2,2),plot(T(:,1),z/1000,T(:,ceil(nt/2)),z/1000,T(:,end),z/1000)
legend(‘Start’,’Mid’,’Finish’)
ylabel(‘Height, km’)
xlabel(‘Temperature’)
title(‘Temperature Profile’)
grid on

% ======== Plot flux and cooling to space curve

figure(4)
subplot(1,2,1),plot(fluxu,zz/1000,fluxd,zz/1000,’Linewidth’,2)
xlabel(‘Flux, W/m^2’)
ylabel(‘Height, km’)
title(‘Up and Down Flux vs Height’)
legend(‘Up’,’Down’)
grid on
subplot(1,2,2),plot(neth,z/1000,’Linewidth’,2,’Color’,’Red’)
xlabel(‘Heating Rate, \circC/day’)
ylabel(‘Height, km’)
title(‘Longwave Radiative Heating Rate’)
grid on

% ======== Plot standard atmosphere data if applicable
if std_atm>0

for i=1:nmol % create a legend for the GHGs on 3rd subplot
textleg2(i)={Hmols(mol(i),:)};
end

figure(5)
subplot(2,2,1),plot(Tinit,z/1000,’Linewidth’,2,’Color’,’Red’)
xlabel(‘Temperature, K’)
ylabel(‘Height, km’)
grid on
title([strtrim(AFGLname(std_atm,:)) ‘, Temperature vs Height’])

subplot(2,2,2),semilogy(Tinit,p/100,’Linewidth’,2,’Color’,’Red’)
set(gca,’YDir’,’reverse’)
xlabel(‘Temperature, K’)
ylabel(‘Pressure, hPa’)
grid on
title([strtrim(AFGLname(std_atm,:)) ‘, Temperature vs Pressure’])

subplot(2,2,3),semilogx(mix*1e6,z/1000,’Linewidth’,2)
xlabel(‘ppmv’)
grid on
title([strtrim(AFGLname(std_atm,:)) ‘, GHG concentration vs Height’])
legend(textleg2)

subplot(2,2,4),semilogx(rho,z/1000,’Linewidth’,2)
xlabel(‘Density, kg/m^3’)
grid on
title([strtrim(AFGLname(std_atm,:)) ‘, Atmospheric Density vs Height’])
end
end
% ==================================================

disp([‘—- End of Run —- ‘ datestr(now) ‘ —-‘]);

===== end of v 0.10.4======

===== Optical_3 ======

function [tau]=optical_3(vt, v, S, iso, gama, nair, molx, numz, numv, p, p0, Tlayer, T0,…
Na, rho, mair, dz, mol, nmol, mix, isoprop, contabs, linewon)

% Optical calculates tau as a function of wavenumber, v and layer,
% i=1-numz-1. Uses data read from the HITRAN database; p, T, rho from each
% layer; mol, mix, ison as the molecules with the mixing ratio and prescribed
% isotopes to consider

% v3 – minor change: mix is now a matrix mix(nmol,numz-1), i.e., mixing ratio for each
% molecule in each layer, & water vapor is now mix(1,1:numz-1), not mixh2o

c1=1.3515; c2=2.4371; c3=5.7216; a1=.8461; a2=0.4308; % coefficients for the
% approximation of transmission with spherical geometry
tau=zeros(numz-1,numv); % optical thickness
fprintf(‘%s’, ‘Calc optical thickness each layer.. ‘);
tau=zeros(numz-1,numv); % optical thickness matrix

% ———– Each layer ————————
for i=1:numz-1
fprintf(‘%d %s’, i, ‘ ‘); % display layer being calculated (this section takes the longest)

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

sig=continuum1000(vt,Tlayer(i)); % get capture cross section in m2/kg
% ** note non-std units ** and values are only valid for
% range defined in function continuum1000

% 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

tau(i,:)=d/pi*tau(i,:); % now include the 3 constants
% Zhao & Shi 2013 approximation (for spherical geometry)
r=1+1./tau(i,:).*log(1+tau(i,:)/2+tau(i,:)./(2*tau(i,:)+2+4.*(c1*tau(i,:).^a1+c2*tau(i,:))/(c3*tau(i,:).^a2+c2*tau(i,:)+1)));
tau(i,:)=r.*tau(i,:);
end % end of each layer, i
disp(‘ ‘);

end

===== end of Optical_3 ======

Related Articles

Part One – some background and basics

Part Two – some early results from a model with absorption and emission from basic physics and the HITRAN database

Part Three – Average Height of Emission – the complex subject of where the TOA radiation originated from, what is the “Average Height of Emission” and other questions

Part Four – Water Vapor – results of surface (downward) radiation and upward radiation at TOA as water vapor is changed

Part Six – Technical on Line Shapes – absorption lines get thineer as we move up through the atmosphere..

Part Seven – CO2 increases – changes to TOA in flux and spectrum as CO2 concentration is increased

Part Eight – CO2 Under Pressure – how the line width reduces (as we go up through the atmosphere) and what impact that has on CO2 increases

Part Nine – Reaching Equilibrium – when we start from some arbitrary point, how the climate model brings us back to equilibrium (for that case), and how the energy moves through the system

Part Ten – “Back Radiation” – calculations and expectations for surface radiation as CO2 is increased

Part Eleven – Stratospheric Cooling – why the stratosphere is expected to cool as CO2 increases

Part Twelve – Heating Rates – heating rate (‘C/day) for various levels in the atmosphere – especially useful for comparisons with other models.

References

The data used to create these graphs comes from the HITRAN database.

The HITRAN 2008 molecular spectroscopic database, by L.S. Rothman et al, Journal of Quantitative Spectroscopy & Radiative Transfer (2009)

The HITRAN 2004 molecular spectroscopic database, by L.S. Rothman et al., Journal of Quantitative Spectroscopy & Radiative Transfer (2005)

Advertisements

Read Full Post »

We’ve looked, via the model, at how radiation travels through, and interacts with, the atmosphere. But this has been for one set of atmospheric conditions which are listed in Part Two.

Water vapor is the most important atmospheric “greenhouse” gas.  But its effect on the surface radiation and TOA radiation are quite complex.

  • Water vapor has absorption lines throughout most of the wavelengths of interest for terrestrial radiation (4-50 μm = 2500-200 cm-1 )
  • Whereas the CO2 concentration is “well-mixed”, i.e. broadly speaking the same mixing ratio or ppmv everywhere in the atmosphere, water vapor is concentrated much more in the lower atmosphere, especially in the planetary boundary layer – so the mixing ratio of water vapor might be 1000 times higher at the surface than at the top of the troposphere
  • The water vapor continuum absorbs as a function of the square of the number of water vapor molecules – other gases like CO2 absorb as a linear function of the number of CO2 molecules

I ran the model described in Part Two many times, changing the boundary layer humidity (BLH), the free tropospheric humidity (FTH) and the surface temperature.

Boundary Layer Humidity

Here is how the TOA flux changes as the boundary layer humidity changes from 20-100% at free tropospheric humidity =60%:

Atmospheric-radiation-6a-Flux-TOA-vs-BLH

Figure 1

Is this surprising?

Now the downward longwave radiation (DLR) at the surface with the same conditions:

Atmospheric-radiation-6b-Flux-DLR-vs-BLH

Figure 2

I could plot some more comparisons with different FTH values, but regardless of the FTH value the graphs of TOA flux vs BLH have the same characteristics.

Why, if water vapor is such a strong GHG, does the top of atmosphere radiation stay almost the same for a boundary layer almost dry through to completely saturated?

The reason is simple. The surface is emitting a blackbody radiation spectrum for that surface temperature (see note 2 in Part Two). The boundary layer is at almost the same temperature as the surface (2-3K difference in this case) so at a given wavelength radiation either gets transmitted, or gets absorbed and re-emitted (usually a combination). But if re-emitted, it is at almost the same temperature as the surface.

The radiative transfer equations show us that the change in flux as radiation travels through a body is caused by the difference between the temperature of the source of the radiation and the temperature of the body in question. This is explained a little more in Part Two and is an essential point to grasp.

So the upward radiation at TOA is the almost the same whether the boundary layer is saturated or dry.

But the surface DLR experiences a big difference as the boundary layer moisture changes. The reason why should be obvious by now, but this subject is quite difficult to take in when you are new to it.

What downward radiation is incident on the boundary layer from above? The answer is – nothing like as much as is incident on the boundary layer from below. The downward atmospheric radiation on this boundary layer is just the emission from the various GHGs (water vapor, CO2, O3, etc) and the water vapor concentration above is quite a lot lower than the boundary layer.

When the boundary layer is saturated with water vapor it emits very strongly across many bands. So the humidity of the boundary layer has a big impact on the DLR, but very little on the TOA outgoing radiation.

Free Tropospheric Humidity

Now let’s see the TOA changes as free tropospheric humidity, FTH, is changed – with BLH fixed at 100%:

Atmospheric-radiation-6c-Flux-TOA-vs-TFH

Figure 3

So now as the atmosphere above the boundary layer gets more water vapor it absorbs (and re-emits) more strongly. But the re-emission is from colder layers of the atmosphere so as this GHG increases in concentration up through the atmosphere, the TOA radiation reduces significantly. And if we reduce the outgoing radiation from the planet then, all other things being equal, the planet warms.

Let’s see the reasons a little more clearly, for the extreme case when we change FTH from 100% to 0%, with Ts=300K. The left side has the spectra for selected heights for both cases, and the right side has the difference between the two cases (at the same heights):

Atmospheric-radiation-6d-Spectrum-TOA-vs-TFH

Figure 4 – Click to enlarge

[Note the title on the top right is slightly incorrect. It is not ΔTOA, it is the Δ (difference) at a few different heights including TOA.]

We can see that in the center of the CO2 band (600-700 cm-1) there is zero change between the two cases at all heights – as expected.

The differences between the two cases occur:

  • strongly around 200-550 cm-1 (50-18 μm) where water vapor absorption is quite strong
  • somewhat around the “atmospheric window” – which still absorbs due to the continuum, especially at the high water vapor saturation pressure when the surface temperature is near 300K
  • strongly around the 1500 cm-1 (7μm) region where there are lots of water vapor absorption lines – however the surface emission is not so high in the first place so the overall effect is reduced

Now let’s look at DLR at the surface:

Atmospheric-radiation-6f-DLR-flux-vs-FTH

Figure 5

Even though the boundary layer is saturated, changes in water vapor above this layer still have a significant impact on the surface DLR.

Let’s see why by comparing the spectra of the 100% and 0% cases at 300K:

Atmospheric-radiation-6e-Spectrum-DLR-vs-FTH

Figure 6 – Click to enlarge

The dark blue curve is the one we are measuring at the surface. The top right is the difference between the two cases at four different heights. The bottom right shows only the difference in spectra at the surface (it is hard to see it in the top right graph).

What is clear is that the difference is caused by the “mid-strength” absorbing/emitting “atmospheric window around 800-1200 cm-1. The “background” DLR from the very top of atmosphere is of course zero (see figure 6 in Part Two). So anything that adds to the DLR on the way down assists the surface spectra.

This is not the case for very strongly emitting regions – see the region around 500 cm-1. The boundary layer emits and absorbs due to water vapor so strongly in this wavenumber region that incident radiation from above is irrelevant to the surface measurement.

With and Without Water Vapor

Now let’s look at the upward spectra at various heights with (saturated) and without (dry) water vapor:

Atmospheric-radiation-6f-Spectrum-TOA-saturated-dry

Figure 7 – Click to enlarge

The bottom right graph is just the top layer in the top right graph shown separately for clarity.

Conclusion

Most people’s untrained intuition about how different radiatively-active gases (=”greenhouse” gases) in different concentrations at different locations change the interaction of radiation with the atmosphere are wrong. Intuition needs to be informed by measurement and theory. It’s just not an intuitive subject.

We’ve seen that water vapor can have very different effects on TOA radiation and the surface. And we’ve seen that water vapor in different places has very different effects. Also we’ve seen the difference between a dry and saturated atmosphere.

Three important points:

1. On a technical note – this model has at least one important flaw, which is to do with how absorption lines change near the top of the troposphere – the Voigt profile vs the Lorenzian profile.

The Voigt profile is not yet implemented because when I last looked at it it made my head hurt trying to implement it in a useful manner (my attempt at the Voigt profile turned a surprisingly fast model given the 287,000 lines calculated in each of 10 layers into a bucket of sludge).

I am going to have another crack at this headache-inducing puzzle before trying to do lots of “what happens when CO2 concentrations are changed by small amounts” scenarios.

Actually, if there are any maths whizzes out there who would like to do the heavy lifting – or even just explain what seems simple to someone who has forgotten almost all maths ever learnt – please let me know here or via email at scienceofdoom – the usual bit – gmail.com. Probably you will get your name on the top of one of the graphs or something, no promises on the font size yet.

2. These scenarios we have seen are all from a “snapshot” of the climate in 1D without running it to a new equilibrium. Obviously if you change the water vapor from 0% to 100% the surface temperature won’t stay the same. Everything will change. When we have been comparing scenarios we have had mostly the exact same surface and atmospheric temperature. This is for a good reason. Small steps first. Actually, grasping atmospheric radiative transfer is a big step.

3. Changes in water vapor affect not just radiative transfer but latent heat and convection. The deep convection that “cranks the engine” on the important tropical circulation is from solar heating over the warmest oceans (see Clouds & Water Vapor – Part Five – Back of the envelope calcs from Pierrehumbert). Radiative transfer is just one piece of the puzzle.

Related Articles

Part One – some background and basics

Part Two – some early results from a model with absorption and emission from basic physics and the HITRAN database

Part Three – Average Height of Emission – the complex subject of where the TOA radiation originated from, what is the “Average Height of Emission” and other questions

Part Five – The Code – code can be downloaded, includes some notes on each release

Part Six – Technical on Line Shapes – absorption lines get thineer as we move up through the atmosphere..

Part Seven – CO2 increases – changes to TOA in flux and spectrum as CO2 concentration is increased

Part Eight – CO2 Under Pressure – how the line width reduces (as we go up through the atmosphere) and what impact that has on CO2 increases

Part Nine – Reaching Equilibrium – when we start from some arbitrary point, how the climate model brings us back to equilibrium (for that case), and how the energy moves through the system

Part Ten – “Back Radiation” – calculations and expectations for surface radiation as CO2 is increased

Part Eleven – Stratospheric Cooling – why the stratosphere is expected to cool as CO2 increases

Part Twelve – Heating Rates – heating rate (‘C/day) for various levels in the atmosphere – especially useful for comparisons with other models.

References

The data used to create these graphs comes from the HITRAN database.

The HITRAN 2008 molecular spectroscopic database, by L.S. Rothman et al, Journal of Quantitative Spectroscopy & Radiative Transfer (2009)

The HITRAN 2004 molecular spectroscopic database, by L.S. Rothman et al., Journal of Quantitative Spectroscopy & Radiative Transfer (2005)

Read Full Post »

If you’ve just stumbled across this article without reading the earlier posts, please take a few minutes to review:

Most people find the actual results of radiative transfer in the atmosphere non-intuitive. Intuition is not a good guide for this topic. So a lot of misconceptions arise because the results of atmospheric physics disagree with the mental models in people’s heads. Obviously the physics must be wrong or probably climate scientists haven’t understood the basics.. Shaking of heads.

For people interested in reality, read on.

We are still looking at how radiation travels and interacts with the atmosphere before anything changes.

There is a lot of fascination in the subject of the “average height of emission” of terrestrial radiation to space. If we take a very simple view, as the atmosphere gets more opaque to radiation (with more “greenhouse” gases) the emission to space must take place from a higher altitude. And higher altitudes are colder, so the magnitude of radiation emitted will be a lesser value. And so the earth emits less radiation and so warms up.

This “average height of emission” is often supplied as a mental model and it’s a good initial starting point.

Here is the result of the atmospheric model created with a surface temperature of 288K (15°C), 80% humidity in the boundary layer and 40% humidity above that (the “free troposphere). This is a cloud-free sample – clouds are very common, but really make life complicated and we are trying to provide a small level of enlightenment. Simple stuff first.

The model is the same as in Part Two – but with 20 layers instead of 10. More layers just means better resolution plus a little bit more accuracy. Each layer contains roughly the same number of molecules (same pressure differential between each layer), so each higher layer is progressively thicker.

The graph shows how much radiation (“flux”) makes it from the surface and from each atmospheric layer in the model to the top of the atmosphere (TOA) – [update Jan 9th, see revised graph in comments].

Atmospheric-radiation-5-Flux-contribution-to-TOA-inc-surface

Figure 1

And here we’ve zoomed in by expanding the x-axis:

Atmospheric-radiation-5-Flux-contribution-to-TOA-excl-surface

Figure 2

The TOA flux = 239.5 W/m², so what is the level where half of this value comes from below and half from above?

If we include the surface and the first 5 layers we don’t have quite half (48%), and if we go to 6 layers we get just over half (51%).  Layer 5 is centered at 1.9km with the top of this layer at 2.1km. Layer 6 is centered at 2.4km.

So let’s say the “average” height of emission to space is just over 2 km (in this example).

There’s probably a better mathematical way of expressing it (this is more like the “median height”) but in fact this “average emission height” is really a curiosity value number anyway. In the words of guru commenter Pekka Pirilä (on another topic):

Any number that is not observable and that’s not used as an input or intermediate value in any calculation that aims to produce observable results is of curiosity value only by definition.

So it’s interesting but you don’t find it a key subject of any climate science papers. Still, being as so many people find it fascinating we will see how it changes as “greenhouse” gases vary in concentration and temperature profiles change.

While we are looking at this, let’s see what wavenumbers from what levels make the largest contribution to the TOA flux. That is, let’s look at the spectral distribution vs height.

First the TOA spectra for these conditions (Ts=288K, Boundary layer humidity=80%, Free tropospheric humidity=40%):

Atmospheric-radiation-5a-Flux-Basic-reference-TOA

Figure 3

Now to see where this all originated from we divide up the wavenumbers into bands of 100 cm-1, and we see the contribution to the TOA flux by band and height in the atmosphere (note that height in km is now ‘lying on the side’ to the left and wavenumber to the right, lost the axes labels somewhere along the way):

Atmospheric-radiation-5a-Flux-3dcontribution-to-TOA

Figure 4

Zooming in a little:

Atmospheric-radiation-5c-Flux-3dcontribution-to-TOA-zoom

Figure 5

We see that in the “atmospheric window” between 800 cm-1 to 1200 cm-1 the surface transmits almost “straight through” (62% of surface flux makes it straight through to the top of atmosphere in this wavenumber range). A small component comes from around the center of the CO2 band (667 cm-1) from the top layer. The rest mostly comes from the “wings” of the CO2 band and where the water vapor absorption is not so strong, around 400 cm-1.

Conclusion

Hopefully seeing the actual data in these different ways helps to see that “average height of emission” is not a real concept or a particularly useful concept. Perhaps it’s a bit like averaging the kg of food consumed per day per person in the entire world. You get a value but the components that made it up are so wide ranging the average has lost anything useful. It’s not like average height of male 20-year olds in Latvia.

Transmission and emission of atmospheric radiation is extremely wavelength dependent.

Related Articles

Part One – some background and basics

Part Two – some early results from a model with absorption and emission from basic physics and the HITRAN database

Part Four – Water Vapor – results of surface (downward) radiation and upward radiation at TOA as water vapor is changed

Part Five – The Code – code can be downloaded, includes some notes on each release

Part Six – Technical on Line Shapes – absorption lines get thineer as we move up through the atmosphere..

Part Seven – CO2 increases – changes to TOA in flux and spectrum as CO2 concentration is increased

Part Eight – CO2 Under Pressure – how the line width reduces (as we go up through the atmosphere) and what impact that has on CO2 increases

Part Nine – Reaching Equilibrium – when we start from some arbitrary point, how the climate model brings us back to equilibrium (for that case), and how the energy moves through the system

Part Ten – “Back Radiation” – calculations and expectations for surface radiation as CO2 is increased

Part Eleven – Stratospheric Cooling – why the stratosphere is expected to cool as CO2 increases

Part Twelve – Heating Rates – heating rate (‘C/day) for various levels in the atmosphere – especially useful for comparisons with other models.

References

The data used to create these graphs comes from the HITRAN database.

The HITRAN 2008 molecular spectroscopic database, by L.S. Rothman et al, Journal of Quantitative Spectroscopy & Radiative Transfer (2009)

The HITRAN 2004 molecular spectroscopic database, by L.S. Rothman et al., Journal of Quantitative Spectroscopy & Radiative Transfer (2005)

Read Full Post »

We reviewed some simple concepts in Part One.

I’ve created a MATLAB model which can do a reasonable job of calculating radiative transfer through the atmosphere. More details about the model to follow, but first let’s look at an actual result and the implications.

There are a whole set of starting conditions, some of which are:

  • 10 layers (of roughly equal pressure change)
  • surface temperature = 288K (15ºC)
  • boundary layer humidity (BLH) = 80%, and boundary layer top of 920hPa
  • free tropospheric humidity (FTH) = 40%
  • lapse rate (the temperature profile in the atmosphere) = 6.5 K/km
  • tropopause at 11.7km, isothermal atmosphere above at 212K and TOA at 50hPa

Atmospheric-radiation-3a

Figure 1

There’s a little too much information when we see all layers, so here are just four (see note 1):

Atmospheric-radiation-3

Figure 2 – Selected layers

What do we see?

Start at the top – the blue line – this is the emission of radiation upwards from the surface. In this case, for simplicity, the surface emissivity = 1.0 (see note 2) so this is the Planck function at 288K. The next curve down is at 2800m up where the temperature has dropped to 270K. The red curve is at 6740m & 244K, and the bottom curve is at 23km & 212K, well into the stratosphere.

Let’s zoom in on one region of wavenumbers/wavelengths:

Atmospheric-radiation-3b-expanded-550-850-499px

Figure 3 – Expanded view

First, the region 640-700 cm-1 (14.3-15.6μm). The upward radiation at each higher altitude (which corresponds to each lower curve in the figure) is at the Planck blackbody function for the temperature of that layer.

The reason is that the incident radiation gets completely absorbed. Nothing gets out the other side. Transmissivity = 0, absorptivity = 1. It is “saturated”.

But we don’t see zero radiation. Why not?

The atmosphere is a strong absorber at these wavelengths, and therefore a strong emitter at these wavelengths. So each layer emits as a blackbody (in this region of wavelengths). We can easily see the temperature of the atmosphere from the Planck function if we are able to measure the radiation from these highly absorbing/emitting wavelengths.

Second, the region near 850 cm-1 (below 12μm). See that the upward radiation at each altitude is almost at the surface radiation value. This is in the “atmospheric window” where the absorption is very low. The atmosphere is almost transparent at these wavelengths. So the absorption is low and the emission is low. But the starting point, if we can use that term, is the emission at the surface temperature of 288K. And so, in this wavelength region, at any point in the atmosphere the upwards radiation is close to the Planck curve of 288K. Basically, the intensity of radiation stays the same as it travels upward through the atmosphere because there is little absorption.

Transmitted Radiation Only

Just to make the subject of emission even clearer, here is a calculation where the atmosphere magically does not emit any radiation – compare this with figure 2:

Atmospheric-radiation-3c-no-emission

Figure 4 – No emission by the atmosphere

Even if we just look at the first 3 layers of the model (1.8km) we get pretty much the same view – i.e., most of the surface radiation is absorbed before we get very far through the atmosphere, but of course it is very wavelength dependent:

Atmospheric-radiation-3c-no-emission-3layers

Figure 5 – No emission by the atmosphere

My calculation says that of 376 W/m² of surface emitted radiation between 200 cm-1 and 2500 cm-1, 75 W/m² (20%) gets transmitted to the top of atmosphere (note 4).

This is not all through the “atmospheric window” – you can see the wavelength dependence in figure 4. I calculate 61 W/m² through the atmospheric window (8-12μm), which means in that wavelength range 62% of surface radiation is being transmitted.

Up and Down Flux

Let’s look at the total (longwave) flux up and down through the atmosphere (note 3):

Atmospheric-radiation-4a-Flux

Figure 6

Notice that the downward flux is zero at the top of atmosphere. This is a boundary condition – there is no (significant) source of longwave radiation coming from outside the atmosphere. As we go down through the atmosphere it gets warmer and so the atmosphere emits more and more. Also as we go down the atmosphere there is much more water vapor, meaning the emissivity of the atmosphere increases significantly. So the atmosphere emits ever more radiation the closer we get to the surface.

We have already considered the upward transmission of radiation. Here the blue line on the graph is simply the sum (the “integral”) of the spectral components we saw in earlier graphs.

Why does the flux reduce with height? Because the absorption of upward radiation is greater than the emission of radiation upwards at each height.

If this point is not clear, please reread this article and Part One – if you are confused over this fundamental point it will be impossible to make good progress in understanding atmospheric radiation.

The absorptivity (the ability of the atmosphere to absorb radiation) is equal to the emissivity (the ability of the atmosphere to emit radiation) at any given wavelength. So why isn’t emission = absorption?

Because the incident upward radiation on a given layer comes from a higher temperature source:

  • Absorption =incident radiation x absorptivity
  • Emission = Planck function (blackbody radiation value) at the temperature of the gas x emissivity

Please ask if this is not crystal clear.

Net Flux & Heating or Cooling

If we want to do any heat transfer calculations we need to look at how the flux changes through the atmosphere. How much radiation enters and how much leaves (see note 5). Anything different from zero for a given layer means there must be heating or cooling by radiation. (This could be balanced by convection – and by absorbed solar radiation).

Let’s see the flux changes for each layer:

Atmospheric-radiation-4b-Net-Flux

Figure 7

What this is showing is the calculation of (radiation in – radiation out) for each layer. As should be obvious from the previous figure, the upward path of longwave radiation is heating the atmosphere (more is absorbed than is emitted), whereas the downward path of longwave radiation is cooling the atmosphere (more is emitted than is absorbed).

When we sum both up we find that the atmosphere is cooling via radiation. “Greenhouse” gases are cooling the atmosphere! If only climate science considered the basics!

Each of the layers in the model contains a similar number of molecules – this is because I divided the atmosphere up into approximately equal pressure sections. This means that 10 W/m² cooling in any layer should equate to similar temperature changes in each layer, but let’s do that calculation anyway (heating rate per unit area/[specific heat capacity x density x depth of layer]):

Atmospheric-radiation-4c-Heating

Figure 8 – Heating (cooling) from longwave radiation

The atmosphere is not actually transparent to solar radiation and you can find similar graphs of Net Shortwave Heating per Day in many climate science textbooks and papers.  The humid lower atmosphere gets a strong solar heating via water vapor. See Atmospheric Radiation and the “Greenhouse” Effect – Part Eleven – Heating Rates.

My graph doesn’t actually reproduce the magnitude of the cooling rates seen for standard atmospheres – typically around 2°C/day in the lower atmosphere but I’m pleased with getting the profile quite similar – remember that the “divergence” is the difference between two values. In this case, the up and down fluxes are in the 200-400 W/m² range, while the net is around 10-20 W/m².

To see what actual difference there is from a more complete model we would need to plug in one of the “standard atmospheres” and compare. The exact profile of water vapor concentration and atmospheric temperature have a big effect – something we will be looking at in detail anyway in later articles.

Convection

The calculation of radiative transfer in the atmosphere can be done for a given profile without knowing anything about convection. That is, if we know where we are right now – without knowing how we got here – we can still do an accurate calculation of how energy moves through the atmosphere by radiation.

If we want to predict the result of how the radiative heating/cooling changes the surface and atmospheric temperature then of course we need a model of convection – and atmospheric circulation.

Conclusion

What the model has done so far is taken:

  • a given temperature profile
  • a given concentration of GHGs including water vapor
  • the large spectroscopic database of absorption lines complied by professionals over decades

– and used basic theory well-known and proven for many decades to calculate the upward and downward path of radiation through the atmosphere.

There’s lots to consider further. But the points and subjects in this article are all fundamental to understanding atmospheric radiation. So if anything is not clear, please ask questions.

Related Articles

Part One – some background and basics

Part Two – some early results from a model with absorption and emission from basic physics and the HITRAN database

Part Three – Average Height of Emission – the complex subject of where the TOA radiation originated from, what is the “Average Height of Emission” and other questions

Part Four – Water Vapor – results of surface (downward) radiation and upward radiation at TOA as water vapor is changed

Part Five – The Code – code can be downloaded, includes some notes on each release

Part Six – Technical on Line Shapes – absorption lines get thineer as we move up through the atmosphere..

Part Seven – CO2 increases – changes to TOA in flux and spectrum as CO2 concentration is increased

Part Eight – CO2 Under Pressure – how the line width reduces (as we go up through the atmosphere) and what impact that has on CO2 increases

Part Nine – Reaching Equilibrium – when we start from some arbitrary point, how the climate model brings us back to equilibrium (for that case), and how the energy moves through the system

Part Ten – “Back Radiation” – calculations and expectations for surface radiation as CO2 is increased

Part Eleven – Stratospheric Cooling – why the stratosphere is expected to cool as CO2 increases

Part Twelve – Heating Rates – heating rate (‘C/day) for various levels in the atmosphere – especially useful for comparisons with other models.

 

References

The data used to create these graphs comes from the HITRAN database.

The HITRAN 2008 molecular spectroscopic database, by L.S. Rothman et al, Journal of Quantitative Spectroscopy & Radiative Transfer (2009)

The HITRAN 2004 molecular spectroscopic database, by L.S. Rothman et al., Journal of Quantitative Spectroscopy & Radiative Transfer (2005)

Notes

Note 1: There is a slight inconsistency in the data presentation. There are 11 boundaries and therefore 10 layers. The spectra are calculated at the boundaries. The water vapor mixing ratio is calculated in the middle of the layer.

The calculation of emission of radiation is based on the temperature and the concentration of each GHG, including water vapor, in the mid-layer (the mid-pressure point in each layer) .

Note 2: The surface emissivity for the ocean, for example, is about 0.96 – see Emissivity of the Ocean. In some parts of the blogworld assuming an emissivity of 1.0 is a heresy that demonstrates what these inappropriately-named “skeptics” have known all along, climate science assumes an “unphysical blackbody model of the world”! And therefore cannot be taken seriously. More exclamation marks and so on.

I could have set an emissivity of 0.96 for the surface and this would have reduced the emitted upward radiation from the surface by 4%. But then for a radiative transfer calculation I would need to reflect 4% of the downward atmospheric radiation upwards (what is not absorbed or transmitted must be reflected). So in fact the upward radiation difference for the two cases (emissivity of 1.0 and 0.96) is quite small, less than 1% and not particularly useful for this calculation.

Note 3: Climate science uses the conventions of shortwave and longwave radiation. Shortwave is wavelengths less than 4μm (wavenumbers greater than 2500cm-1), while longwave is greater than 4μm.

99% of solar radiation is shortwave, while 99% of all terrestrial radiation is longwave. This makes it easy to separate the two. See The Sun and Max Planck Agree – Part Two.

Note 4: The emission of thermal radiation by a surface at 288K with an emissivity of 1.0 is 390 W/m². This is across all wavelengths. The model looks at the range of wavenumbers that equates to 4-50μm to ease up the calculation effort required. Almost all of the “missing spectrum” is in the far infra-red (longer wavelengths/lower wavenumbers), and is subject to relatively high absorption from water vapor.

Note 5: If you hear the technical term flux divergence it is essentially the same thing. Flux divergence is per unit volume so it isn’t such a useful value. Instead the most common term is heating rate which divides the gain (loss) in radiation energy by the heat capacity to calculate the radiative heating (cooling) rate per unit of time (typically per day).

Read Full Post »

This could be considered as a continuation of the earlier series – Atmospheric Radiation and the “Greenhouse” Effect – but I’ve elected to start a new series.

It’s clear that many people have conceptual problems with the subject of what is, in technical terms called radiative transfer. That is, how radiation travels through the atmosphere and is affected by the atmosphere.

Radiation and Gas in a Box

First, let’s consider what happens as we shine an intense beam of infrared radiation at a narrow range of wavelengths (let’s say somewhere in the region of 15μm) through a box of CO2 gas at room temperature:

Atmospheric-radiation-2

Figure 1

The red arrow on the left is the incident radiation. The graph indicates the spectrum. The spectrum on the right is made up of two main parts:

  • the transmitted radiation – the incident radiation attenuated by the absorbing gas
  • the emitted radiation due to the temperature and emissivity of the gas at these wavelengths

The yellow spectrum shows what we would measure from one of the sides. Note that the transmitted radiation that goes from left to right has no effect on this yellow spectrum (except in so far as absorption of the incident radiation affects the temperature of the gas).

If we increase the length of the box (left to right) – and keep the density the same – the transmitted radiation from the right side would decrease in intensity. If we reduce the length of the box (again, same density) the transmitted radiation from the right would increase in intensity.

But the emitted radiation from the top is only dependent on the temperature of the gas and its emission/absorption lines.

And the temperature of the gas is of course affected by the balance between absorbed and emitted radiation as well as any heat transfer from the surroundings via convection and conduction.

Hopefully, this is clear. If anyone thinks this simple picture is wrong, now is the time to make a comment. Confusion over this part means that you can’t make any progress in understanding atmospheric radiation.

Scattering is insignificant for longwave radiation (4μm and up). Stimulated emission is insignificant for intensities seen in the atmosphere.

Radiation in the Atmosphere

How does radiation travel through the atmosphere?

Atmospheric-radiation-1

Figure 1

The idea shown here is a spectrum of radiation at different wavelengths incident on a “layer” of the atmosphere (see note 1). The atmosphere has lots of absorption lines of many different strengths. As a result the transmitted radiation making it out of the other side is some proportion of the incident radiation. The proportion varies with the wavelength.

The atmosphere also emits radiation, and the emission lines are the same as the absorption lines. More about that in Planck, Stefan-Boltzmann, Kirchhoff and LTE.

However, the emission depends on the temperature of the gas in the layer (as well as the absorption/emission lines). But the absorption depends on the intensity of incident radiation (as well as the absorption/emission lines), which in turn depends on the temperature of the source of the radiation.

So in almost every case, the sum of transmitted plus emitted radiation is not equal to the incident radiation. By the way, the spectrum at the top is just a raggedy freehand drawing to signify that the outgoing spectrum is not like the incoming spectrum. It’s not meant to be representative of actual intensity vs wavelength.

And – it’s a two way street. I only showed one half of the story in figure 1. The same physics affects downward radiation in exactly the same way.

Considering One Wavelength at a Time

To calculate the actual transmission of radiation through the layer we simply work out the transmissivity, tλ, of the layer at each wavelength, λ (tλ simply indicates that t will vary for each value of λ we consider). We do that by looking up values calculated by spectroscopic professionals. These values are per molecule, or per kg of particular molecules so we need to find out how much of each absorbing gas is present.

1. The incident radiation making it through the layer = Iλ x tλ – for example, it could be 90% making it through, or 20%.

2. The “new” radiation emitted from each side of the layer equals the “Planck blackbody function at the temperature of the layer and the wavelength of interest” x Emissivity of the gas at that wavelength.

In case people are interested this can be written as Bλ(T).ελ, where ελ = emissivity at that wavelength, and Bλ(T) is the “Planck function” at that temperature and wavelength. Well, emissivity = absorptivity (at the same wavelength) and absorptivity = 1-transmissivity, so the same equation can be written as Bλ(T).(1-tλ).

Digression

Perhaps (the Planck function showing up in an equation) this is where many blogs (Parady blogs?) get the idea, and promote and endorse the idea, that climate science depends on the assumption that the atmosphere emits as a blackbody. There are some cases where the atmospheric emission is not far from the “blackbody assumption” (e.g., in clouds), but that is due to reality not assumption. There is no “blackbody assumption for the atmosphere” in climate science. But there is a movement of people who believe it to be true.

Misleadingly, they like to be known as “skeptics”.

End of digression..

Doing the Calculation

So it’s not really that hard to understand how radiation travels through the atmosphere. It is difficult to calculate it, mostly due to having to read a million absorption lines, figure out the correct units, get a model of the atmosphere (temperature profile + concentration of different “greenhouse” gases at each height), write a finite element program and work out a solution.

But that is just tedious details, it’s not as hard as having to understand general relativity or potential vorticity.

One important point – it is not possible to do this calculation in your head. If you think you have done it in your head, even to a close approximation, please go back and read this section again, then look up all of the absorption lines.

Still convinced – post your answer in a comment here.

In the next article I’ll explain radiative-convective models and show some results from the atmospheric model I built in MATLAB which uses the HITRAN database. Now that we have a model which calculates realistic values for emission, absorption and transmission we can slice and dice the results any way we want.

Does water vapor mask out the effects of CO2? What proportion of radiation is transmitted through the atmospheric window? What is the average emission height to space?

Related Articles

Part Two – some early results from a model with absorption and emission from basic physics and the HITRAN database

Part Three – Average Height of Emission – the complex subject of where the TOA radiation originated from, what is the “Average Height of Emission” and other questions

Part Four – Water Vapor – results of surface (downward) radiation and upward radiation at TOA as water vapor is changed

Part Five – The Code – code can be downloaded, includes some notes on each release

Part Six – Technical on Line Shapes – absorption lines get thineer as we move up through the atmosphere..

Part Seven – CO2 increases – changes to TOA in flux and spectrum as CO2 concentration is increased

Part Eight – CO2 Under Pressure – how the line width reduces (as we go up through the atmosphere) and what impact that has on CO2 increases

Part Nine – Reaching Equilibrium – when we start from some arbitrary point, how the climate model brings us back to equilibrium (for that case), and how the energy moves through the system

Part Ten – “Back Radiation” – calculations and expectations for surface radiation as CO2 is increased

Part Eleven – Stratospheric Cooling – why the stratosphere is expected to cool as CO2 increases

Part Twelve – Heating Rates – heating rate (‘C/day) for various levels in the atmosphere – especially useful for comparisons with other models.

Notes

Note 1 – What is a layer of atmosphere? Isn’t the thickness of this layer somewhat arbitrary? What if we change the thickness? And doesn’t radiation go in all directions, not just up?

These are all good questions.

In typical physics terms the actual equation of “radiative transfer” is a differential equation, which expresses continual change. In practical terms, solving a differential equation in most real world cases requires a numerical solution which has finite thicknesses for each layer.

People trying to solve these kind of problems usually check what happens to the solution as they go for more of thinner layers vs less of thicker layers. There is a trade-off between accuracy and speed.

Radiation does go in all directions. The plane parallel assumption has very strong justification and – in simple terms – mathematically resolves to a vertical solution with a correction factor. You can see the plane parallel assumption and the derivation of the equations of radiative transfer in Understanding Atmospheric Radiation and the “Greenhouse” Effect – Part Six – The Equations.

Read Full Post »

« Newer Posts