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

[update Jan 25th, see comment]

[update Jan 20th see comment]

[update Jan 14th – see comment]

[original version published]

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 m^{2}/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/(m^{2}.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)

on January 10, 2013 at 2:51 pm |Pekka PiriläFor a free clever text editor capable of using colors for Matlab-code I recommend Notepad++. Naming the file as hitran_0_9_21.m helps in that.

on January 10, 2013 at 9:21 pm |Pekka PiriläSoD,

Your code is calling functions

Hitran_read_0_2()

define_atmos_0_5()

continuum1000()

where are they from?

on January 10, 2013 at 9:46 pm |scienceofdoomI updated the article with the 3 functions.

on January 10, 2013 at 9:58 pmPekka PiriläAt least one more function seems to be missing: satvaph2o().

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

on January 10, 2013 at 10:12 pmscienceofdoomI added this function in the main article.

on January 10, 2013 at 10:33 pmPekka PiriläNow I got the code running. After calculation of all layers it turned out that planckmv() is missing.

on January 10, 2013 at 10:48 pmscienceofdoomAdded.

There’s probably some useful Matlab tools to tell me what functions are called..

on January 10, 2013 at 10:55 pmPekka PiriläNow I got to the end and got a few plots. I tested with 4 layers only.

But it’s late here and I continue tomorrow.

on January 10, 2013 at 11:00 pmscienceofdoomGreat news. Look forward to some results.

on January 11, 2013 at 7:46 pmPekka PiriläSoD,

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.

on January 10, 2013 at 10:11 pm |scienceofdoomI’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.

on January 11, 2013 at 6:01 am |gallopingcamelIt has puzzled me that you analyze radiative energy transfer effects in exquisite detail while ignoring simpler theories that explain observations much better.

I entered the name “Nikolov” in your search bar and nothing came back. Surely you must have discussed N&K’s poster if only to refute it.

on January 11, 2013 at 6:45 am |scienceofdoomgallopingcamel,

Actually it’s not really exquisite detail.

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.

Providing clarity is essential. This is the aim of this series.

Do you mean Unified Theory of Climate by Nikolov & Zeller?

Is the link above their paper or is a formal paper available?

A core part of overturning a theory is understanding what the theory is that you think you are overturning.

If this pdf is their ground-breaking theory they could do with reading a textbook to find out what basic atmospheric physics teaches.

Still this is all part of my prediction in New Theory Proves AGW Wrong. Many people continue to entertain, and there will be no shortage of future vaudeville acts on the “greenhouse” effect.

Perhaps I will write an article so their theory can be discussed. Let me know the paper in question.

on January 12, 2013 at 4:00 amRWSoD,

BTW, what’s the objection to L&C 2011?

on January 12, 2013 at 7:02 amscienceofdoomRW,

What have I said about L&C 2011? And I assume that’s Lindzen & Choi 2011?

on January 12, 2013 at 2:25 pmRWI don’t know. I’m just asking. Have you read it?

on January 12, 2013 at 8:14 pmscienceofdoomRW,

I went back through my notes and I read Lindzen & Choi 2011, along with:

Relationships between tropical sea surface temperature and TOA radiation, Trenberth et al 2010Lindzen & Choi 2009

Spencer & Braswell 2010

The Climate sensitivity and its components diagnosed from Earth radiation budget data, Forster & Gregory 2006An observationally based estimate of the climate sensitvity, Gregory et al 2002A new method for diagnosing radiative forcing and climate sensitivity, Gregory et al 2004Spencer & Braswell 2008

Transient climate response estimated from radiative forcing and observed temperature change, Gregory & Forster 2008Estimations of climate sensitivity based on TOA radiation imbalance, Lin et al 2010Can climate sensitivity be estimated from short-term relationships at TOA net radiation and surface temperature, Lin et al 2011Cloud feedbacks in the climate system: A review Article, Stephens 2005On the accuracy of deriving climate feedback parameters from correlations between surface temperature and outgoing radiation, Murphy & Forster 2010.I wrote a few articles, mainly around Spencer & Braswell’s criticism of whether it is possible to measure anything useful.

on January 13, 2013 at 6:52 amscienceofdoomRW,

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.

on January 13, 2013 at 3:02 pmRWSoD,

OK, fair enough, but why is everyone operating under the assumption that the feedback is most likely positive rather than negative?

on January 13, 2013 at 3:09 pmRWThat seems dubious to me to say the very least.

on January 12, 2013 at 4:20 pm |DeWitt Paynegallopingcamel,

The internal search bar isn’t all that useful. Google, OTOH, produced this comment:

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

on January 12, 2013 at 6:55 am |Visualizing Atmospheric Radiation – Part Six – Technical on Line Shapes « The Science of Doom[…] « Visualizing Atmospheric Radiation – Part Five – The Code […]

on January 12, 2013 at 7:03 pm |Pekka PiriläSoD,

You use the constant diffusivity factor of 1.66 to take approximately into account the zenith angle dependence of the pathlength. There’s a recent paper by Zhao and Shi

http://www.sciencedirect.com/science/article/pii/S1350449512000850

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

on January 13, 2013 at 6:36 pm |scienceofdoomPekka,

This is an interesting paper, thanks. I will absorb and look at implementing.

on January 13, 2013 at 6:55 pmPekka Pirilä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.

on January 13, 2013 at 10:39 pmPekka Pirilä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.

on January 13, 2013 at 4:59 am |Visualizing Atmospheric Radiation – Part Seven – CO2 increases « The Science of Doom[…] 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. […]

on January 13, 2013 at 9:59 am |Pekka PiriläSoD,

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

————————–

by the lines

————————–

dt1=1./((vt-v(im(j))).^2+ga^2);

tau(i,:)=tau(i,:)+nummol*S(im(j))*ga*dt1;

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.

on January 13, 2013 at 10:09 am |Pekka Pirilä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.

on January 13, 2013 at 4:42 pm |scienceofdoomPekka,

Very nice changes.

I am updating the code to v0.9.5 (0.9.4 orphan version used for John Eggert’s problem).

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

on January 13, 2013 at 4:57 pmPekka PiriläSoD,

Your change to the continuum part is exactly, what I had in mind.

on January 13, 2013 at 5:06 pmscienceofdoomI also updated the printing section –

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

% Print out various flux values

disp(‘ ‘);

disp([‘Flux up at TOA = ‘ num2str(fluxu,’%4.1f’) ‘ W/m^2, Flux down at surface = ‘ …

num2str(fluxd,’%4.1f’) ‘ W/m^2′]);

disp([‘TOA transmitted flux= ‘ num2str(TOAtr,’%4.1f’) ‘ W/m^2, Surface emission= ‘ …

num2str(sum(rads),’%4.1f’) ‘ W/m^2, Window transmission= ‘ …

num2str(sum(radu1(numz,iw))/sum(rads(iw))*100,’%4.1f’) ‘ %’]);

———————————

on January 13, 2013 at 5:16 pmscienceofdoomPekka,

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.

on January 13, 2013 at 6:03 pmPekka Pirilä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.

on January 14, 2013 at 8:00 am |scienceofdoomI reran the set of results for changing CO2 from 280ppm, 360, 400, 450, 560 ppm at Ts=288K and 300K with the diffusivity factor (1.66) applied to the continuum.

[Pekka, I note your comments of January 13, 2013 at 6:03 pm, need a little time to review and digest].

The results were run at Δv = 1 cm

^{-1}& Δv = 0.1 cm^{-1}. The differences between the two cases are generally small, somewhere around 0.2 – 0.5 W/m^{2}for TOA and 0.1 – 0.3 W/m^{2}for DLR.Then I reran two cases with Δv = 0.01 cm

^{-1}– Ts=300K & CO2 = 280 ppm & 560 ppm.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/m^{2}.on January 14, 2013 at 9:34 am |Pekka PiriläSoD,

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

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

on January 23, 2013 at 11:39 amscienceofdoomPekka,

I’ve taken a look at your code (and the paper). Your code is a step ahead of the paper as you said.

Do you have some results to share comparing the std diffusivity approximation and your more accurate integration over zenith angle?

on January 23, 2013 at 12:54 pmscienceofdoomI’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.

replace diff=1.66; with

c1=1.3515; c2=2.4371; c3=5.7216; a1=.8461; a2=0.4308; % coefficients for the

% approximation of transmission with spherical geometry

…

then:

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

on January 23, 2013 at 3:33 pmPekka PiriläSoD,

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.

on January 23, 2013 at 6:02 pmPekka Pirilä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.

on January 14, 2013 at 12:16 pm |scienceofdoomI updated the code to v0.9.6 – and added a link to a new text file in the main article, and pasted the updated code into the article.

The changes are described briefly in the comments in the code:

a) a few corrections as per Pekka’s comments, the efficiency change in the main loop as per Pekka’s comments, the correction of the diffusivity factor for the continuum.

b) tidying up the output printing to the screen, and the input to the function so that changes in GHGs introduced does not require changes to the code. Examples to follow.

on January 14, 2013 at 12:28 pm |scienceofdoomChanges to the code in the reading of data are:

1. How parameters get passed in:

e.g.

HITRAN_0_9_6(1, 11, [1 2 4 6], [0 280e-6 319e-9 1775e-9], [2 3 2 2], 9.2e4, .8, .4, 1, 288, 6.5, 6.5, 2e4, 5e3, 12*3600, 2, 1);

this says – “[1 2 4 6]” = use H2O, CO2, N2O, CH4

“[0 280e-6 319e-9 1775e-9]” = at mixing ratios of n/a, 280ppm, 319ppb, 1775ppb ” [2 3 2 2]” = the first and second isotopologues of H2O, N2O, CH4 and first 3 isotopologues of CO2.

These 3 vectors should all be the same length of course (1 row x number of GHG molecules being used).

I did a few tests and everything seemed to work ok.

2. How the files get read:

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

on January 14, 2013 at 5:12 pm |Pekka PiriläSod;

A major part from the middle including the main loops seems to be missing from the above comment but not from the downloadable file.

on January 14, 2013 at 6:55 pmscienceofdoomWell picked up. I had another go and the same result. It looks like WordPress has decided there is a maximum lengths for comments (possibly a setting I can change).

Instead I pasted the new code into the article body.

on January 14, 2013 at 11:29 pm |Visualizing Atmospheric Radiation – Part Eight – CO2 Under Pressure « The Science of Doom[…] we have a model (described in Part Five) which incorporates the actual line by line absorption and emission of various GHGs in a realistic […]

on January 20, 2013 at 5:30 am |scienceofdoomI’ve updated this article with the code for v0.10.1.

The article – Part Nine – Reaching Equilibrium – shows some runs from this.

The main changes:

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.

on January 20, 2013 at 7:34 am |Visualizing Atmospheric Radiation – Part Nine – Reaching Equilibrium « The Science of Doom[…] The atmospheric model is described in brief in Part Two and in a comment, then in detail in Part Five – The Code. […]

on January 21, 2013 at 9:17 am |scienceofdoomPekka,

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.

on January 21, 2013 at 10:00 am |Pekka PiriläSoD,

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.

on January 21, 2013 at 10:52 amscienceofdoomPekka,

After reviewing I think the calculation is correct. My notes were sparse, probably it seemed obvious at the time.

Here’s my reasoning for cross-checking when you get the opportunity:

Pierrehumbert notes the equation for the equivalent path in p.260. Then he gives the example of a 1km layer at 300K which he says has an equivalent path of 9.3 kg/m

^{2}.When I calculate (psat/p0).(psat/Rw.T).z

I get 3528

^{2}.1000/(10^{4}.461.4 .300) = 8.99[In this case p0 is the reference pressure where the continuum absorption coefficient was measured – in this case 100 mbar = 10

^{4}Pa, and Rw = gas constant for water vapor = 8.3143/.1802 = 461.4].So this part of the check seems reasonable – because the saturation vapor pressure can be selected as slightly different numbers dependent on the empirical formula used.

The formula in Matlab is:

τ(v) = (psat/p0).(psat/Rw.T) . σ(v) . d/100 . diffusivity coefficient, where d = path length in cm (so d/100 just converts it to meters), σ is in units of m

^{2}/kg.on January 25, 2013 at 9:25 pmPekka PiriläSoD,

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.

on January 25, 2013 at 10:54 pmscienceofdoomPekka,

Thanks, the use of 4 instead of 4.25 was an oversight on my part, which I have now fixed.

on January 23, 2013 at 9:39 am |Visualizing Atmospheric Radiation – Part Ten – “Back Radiation” « The Science of Doom[…] 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. […]

on January 24, 2013 at 7:26 am |scienceofdoomPekka,

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 ====

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

on January 24, 2013 at 10:47 am |Pekka PiriläSoD,

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.

This figure shows the difference of the two approaches (unit mW/m^2).

This figure compares the difference of the two approaches to the simpler model (unit W/m^2).

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.

on January 24, 2013 at 1:00 pmPekka Pirilä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.

on January 24, 2013 at 8:40 am |UliSoD, 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?

on January 24, 2013 at 9:24 am |scienceofdoomUli,

I’m uploading them as I write this comment. They are quite large so it will take a while, they will be in the body of the article.

Because of WordPress file type upload limitations I had to rename them to .doc – so change each file back to .par after saving locally.

If you produce any results and show them elsewhere, please cite,

The HITRAN 2008 molecular spectroscopic database, by L.S. Rothman et al, Journal of Quantitative Spectroscopy & Radiative Transfer (2009).on January 24, 2013 at 10:05 amscienceofdoomFiles for water vapor, CO2, O3, N2O and CH4 are now uploaded in the body of the article.

on January 24, 2013 at 11:04 amUliSoD, thanks.

In between I was able to download the missing data from

http://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.

on January 24, 2013 at 11:11 amscienceofdoomUli,

The wavenumber range of interest is set by vmin & vmax as defined in the program, currently this is 200 – 2500 cm

^{-1}.Separately, as a consequence, any comparison to blackbody emission is within this range.

on January 24, 2013 at 4:08 pm |UliSoD,

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.

on January 24, 2013 at 5:52 pm |Pekka PiriläUli ans Sod,

Increasing the value of convloop seems to be a real improvement that corrects gradually the overshoot of lapse rate in the troposphere. A large value is needed to get that done well.

It seems, however, to be relatively easy to do the lapse rate correction exactly getting rid of the whole loop.

on January 24, 2013 at 8:10 pmscienceofdoomFor 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.

on January 24, 2013 at 8:47 pmPekka PiriläSoD,

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.

I’ll send you the code, when I have got the fluxes in order.

on January 25, 2013 at 11:25 amPekka Pirilä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)

% ==== Calculate adiabatic temperatures and heat capacities ======

AT=zeros(1,numz);

CC=zeros(1,numz);

ACQ=zeros(1,numz);

CC(1)=Cps;

for i=2:numz

AT(i)=-z(i-1)*lapse; % Temperatures relative to surface

% with adiabatic lapse rate

CC(i)=CC(i-1)+Cp(i-1); % Cumulative heat capacity up to layer i-1

ACQ(i)=ACQ(i-1)+AT(i)*Cp(i-1); % Cumulative contribution to heat

% content from AT up to layer i-1

end

% ==== Finished adiabatic temperatures and cumulative heat contents ======

Then the whole convloop section is replaced by the following

% == Calculate convective adjustments based on fixed lapse rate

% calculate the cumulative heat contents after radiative calculation

CQ=zeros(1,numz);

CQ(1)=Ts(h)*Cps;

for i=2:numz

CQ(i)=CQ(i-1)+T(i-1,h)*Cp(i-1);

end

% search top level of adjustment

% calculate surface temperatures that lead to the same

% cumulative heat content up to the level i-1

TsT=zeros(1,numz);

TsT(1)=Ts(h);

for i=2:numz

TsT(i)=(CQ(i)-ACQ(i))/CC(i);

end

% compare temperatures at the top of adjusted range for each

% height of this range. Change in sign of the difference tells

% the correct height for adjustment

adjtop=0;

for i=2:numz

if TsT(i)+AT(i)>T(i-1,h)

adjtop=i-1;

end

end

% calculate changes in heat content and the convective fluxes

% replace temperatures by adjusted values

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.

on January 25, 2013 at 12:14 pmUliPekka,

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.

hitran_0_10_1(1, 11, [1 2 4 6], [0 280e-6 319e-9 1775e-9], [2 3 2 2], 9.2e4, .8, .4, 1, 1 ,288, 6.5, 6.5, 2e4, 5e3, 0.5 ,239 ,12*3600, 2*360, 1);

hitran_0_10_1(1, 11, [1 2 4 6], [0 280e-6 319e-9 1775e-9], [2 3 2 2], 9.2e4, .8, .4, 1, 1 ,288, 6.5, 6.5, 2e4, 5e3, 0.5 ,239/2 ,12*3600, 2*360, 1);

hitran_0_10_1(1, 11, [1 2 4 6], [0 280e-6 319e-9 1775e-9], [2 3 2 2], 9.2e4, .8, .4, 1, 1 ,288, 6.5, 6.5, 2e4, 5e3, 0.5 ,139 ,12*3600, 2*360, 1);

on January 25, 2013 at 12:53 pmPekka PiriläUli,

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.

on January 25, 2013 at 3:31 pmPekka PiriläUli,

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.

on January 25, 2013 at 5:30 pmUliPekka,

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.

on January 25, 2013 at 6:05 pmPekka Pirilä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.

on January 25, 2013 at 11:01 pmscienceofdoomI haven’t tried Pekka’s new code yet but I am not seeing the same problems with getting the lapse rate to match the set value.

I’ll have a closer look soon.

on January 25, 2013 at 11:22 am |scienceofdoomI 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:

on January 25, 2013 at 11:25 am |scienceofdoomFor the example above here is my input statement to the function:

[ vt tau fluxu fluxd ztropo z dz p rho mixh2o zb pb Tb rhob …

Tinit T Ts dE dEs dCE radu radd emitu sheat surfe TOAe TOAf TOAtr alr]…

=HITRAN_0_10_3(1, 21, .9, [1 2], [0 280e-6], [0 560e-6], [1 1 1 1 1], 9.2e4, .8, .4, 6e-6, 1, 1, …

288, 6.5, 6.5, 2e4, 100, 10, 242, 3600*2, 16000, 8000, 1);

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.

on January 25, 2013 at 11:39 am |scienceofdoomOne other point I realized, for more explanation in another post -my understanding was flawed of whether something called “radiative forcing” was actually produced by the program.

The starting point of a given profile is not in equilibrium – as shown in Part Nine – Reaching Equilibrium.

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.

When I did this for 288K, lapse rate=6.5K/km and the above humidity and other conditions the radiative forcing on CO2 change from 280 -> 560ppm was 4 W/m

^{2}.It is too close to the IPCC value not to suggest a fluke at this stage, more evaluation needed.

on January 26, 2013 at 2:04 am |scienceofdoomUli,

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.on January 27, 2013 at 1:08 pm |Visualizing Atmospheric Radiation – Part Eleven – Stratospheric Cooling « The Science of Doom[…] 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. […]

on January 28, 2013 at 10:07 am |scienceofdoomI am interested in doing some comparisons and finally found the detailed data for the 6 std atmospheres often cited for this purpose:

– Tropical

– Mid-latitude summer

– Mid-latitude winter

– Sub-arctic summer

– Sub-arctic winter

– US standard

The profiles are explained (and the values listed in tables) in AFGL atmospheric constituent profiles (0.120 km), by GP Anderson, SA Clough & others, 1986.

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

This file contains:

– height, AFGLz (1×50)

– pressure, AFGLp (6×50 – same for all following)

– temperature, AFGLt

– density, AFGLrho

– water vapor mixing ratio, AFGLh2o

– O3 mixing ratio, AFGLo3

– N2O mixing ratio, AFGLn2o

– CH4 mixing ratio, AFGLch4

The mixing ratios are all in ppv.

If any non-MATLAB people need the data I can supply the .csv files which are 300 values, comma separated. Just ask.

I plan to assess the longwave cooling from each radiatively-active gas using these profiles.

It will need some minor changes to the MATLAB program to allow a defined pressure/temperature/density atmosphere – and to allow the mixing ratio to be defined for each gas as a function of altitude.

on January 28, 2013 at 10:16 am |scienceofdoomAny FORTRAN experts out there??

The LBLRTM code has been very well-verified against detailed measurements.

It uses the Voigt profile as standard, which allows better modeling of the stratospheric line absorption and emission (see Part Six – Technical on Line Shapes).

The bit that claims to do the Voigt profile

And two other (probably) important files for reference:

The Main Bit, I think

Atmospheric ray tracing

I’m interested in implementing the Voigt profile but the standard approaches are very computationally expensive. Maybe this program can shed some light on how to do it..

on January 28, 2013 at 5:33 pm |Pekka PiriläSoD,

Have you had a look on this paper by Zaghloul and Ali

http://dl.acm.org/citation.cfm?id=2049679&picked=formats&CFID=266871991&CFTOKEN=87480039

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.

on January 28, 2013 at 5:40 pm |Pekka PiriläThe Wikipedia page

http://en.wikipedia.org/wiki/Voigt_profile

tells, how the Faddeyeva function is related to the Voigt profile.

on January 29, 2013 at 9:31 am |Pekka Pirilä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.)

on January 29, 2013 at 4:54 pm |Pekka PiriläThe Voigt code seems to be working.

Little effect is expected at pressure levels well above 1 hPa. Therefore I extended the calculation 0.1 hPa (or 10 Pa) and used 50 layers with the uppermost at 46.9-58.8 km.

The temperature profiles deviate in two top layers significantly, otherwise very little changed.

on January 30, 2013 at 8:27 am |Visualizing Atmospheric Radiation – Part Twelve – Heating Rates « The Science of Doom[…] numerous improvements – outlined in Part Five – The Code, I got around to adding some “standard atmospheres” so we can see some comparisons and […]

on January 30, 2013 at 9:38 am |Pekka PiriläSoD,

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

function ans = HITRAN_0_10_3ps(args)

added at the beginning of the code the lines

dv = args.dv;

numz = args.numz;

ratiodp = args.ratiodp;

mol = args.mol;

mix = args.mix;

newmix = args.newmix;

ison = args.ison;

blp = args.blp;

BLH = args.BLH;

FTH = args.FTH;

STaH = args.STaH;

contabs = args.contabs;

linewon = args.linewon;

Tsi = args.Tsi;

lapsekm = args.lapsekm;

lapseikm = args.lapseikm;

ptropo = args.ptropo;

minp = args.minp;

ocd = args.ocd;

src = args.src;

tstep = args.tstep;

nt = args.nt;

ntc = args.ntc;

voigtmaxp = args.voigtmaxp;

ploton = args.ploton;

pin = args.pin;

convloop = args.convloop;

and at the end of the code the lines

ans=struct(‘vt’, vt, ‘tau’, tau, ‘fluxu’, fluxu, ‘fluxd’, fluxd, …

‘ztropo’, ztropo, ‘z’, z, ‘dz’, dz, ‘p’, p, ‘rho’, rho, ‘mixh2o’, mixh2o, …

‘zb’, zb, ‘pb’, pb, ‘Tb’, Tb, ‘rhob’, rhob, ‘Tinit’, Tinit, ‘T’, T, …

‘Ts’, Ts, ‘dE’, dE, ‘dEs’, dEs, ‘dCE’, dCE, ‘radu’, radu, ‘radd’, radd, …

‘emitu’, emitu, ‘sheat’, sheat, ‘surfe’, surfe, ‘TOAe’, TOAe, …

‘TOAf’, TOAf, ‘TOAtr’, TOAtr, ‘alr’, alr);

Then I created the first set of argument structure taking advantage from earlier argument values by the command

args1=struct(‘dv’, 1, ‘numz’, 51, ‘ratiodp’, 0.9, ‘mol’, mol, ‘mix’, mix, …

‘newmix’, newmix, ‘ison’, ison, ‘blp’, blp, ‘BLH’, BLH, ‘FTH’, FTH, …

‘STaH’, STaH, ‘contabs’, contabs, ‘linewon’, 1, ‘Tsi’, Tsi, ‘lapsekm’, 6.5, … ‘lapseikm’, 6.5, ‘ptropo’, ptropo, ‘minp’, 100, ‘ocd’, ocd, ‘src’, src, …

‘tstep’, tstep, ‘nt’, 500, ‘ntc’, 1000, ‘ploton’, 0, ‘pin’, [0], ‘convloop’, 1, …

‘voigtmaxp’, 1000);

and started the run by

ans1=hitran_0_10_3ps(args1);

The name of the argument structure and the answer structure can be defined to describe the nature of the case. Thus it’s easy to figure out later, what they are about.

The structures args1 and ans1 can also be used as arguments for separate plotting functions.

on January 30, 2013 at 10:43 am |scienceofdoomThe latest version of the code posted – v0.10.4 adds 6 std atmospheres as options. Some results in Part Twelve – Heating Rates.

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

on January 30, 2013 at 2:20 pm |Pekka PiriläSoD,

I’m testing my extensions with the new code and have made the required changes in both the basic model code and in my additions.

In doing that I noticed one error in handling mix.

You have the line

mix=repmat(mixc,numz,1);

it should be

mix=repmat(mixc’,1,numz);

on January 30, 2013 at 2:23 pm |Pekka PiriläThere as also another new error. The arrays raddpre, radupre, raddpost, and radupost are undefined, when own profiles are used.

on January 30, 2013 at 10:02 pmscienceofdoomPekka,

I am not sure if this error was here all along (for cases where the GHG concentrations don’t change).

Anyway, I now initialized these variables in the code and updated the file on this page (with the same filename unfortunately).

Note that numz gets defined

duringthe program if a std atmosphere is used, so any vector that needs initialization with ‘numz’ must get initialized after the atmosphere setup.You are correct about mix & mixc, and I have changed that also.

I updated the AFGL file with SI units (meters for z, Pa for pressure) and removed rho because it was not necessary. This file is now in the body of the article.

on January 30, 2013 at 10:07 pm |scienceofdoomI also changed fluxu and fluxd in 0.10.4.

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.

on February 2, 2013 at 1:06 am |Kiehl & Trenberth and the Atmospheric Window « The Science of Doom[…] the basics of a MATLAB line by line calculation of radiative transfer in the atmosphere. And Part Five – The Code gives the specifics, including the […]

on February 2, 2013 at 7:46 am |Visualizing Atmospheric Radiation – Part Four – Water Vapor « The Science of Doom[…] Visualizing Atmospheric Radiation – Part Three – Average Height of Emission Visualizing Atmospheric Radiation – Part Five – The Code […]

on February 2, 2013 at 7:48 am |Visualizing Atmospheric Radiation – Part Three – Average Height of Emission « The Science of Doom[…] Part Five – The Code – code can be downloaded, includes some notes on each release […]

on February 2, 2013 at 7:51 am |Visualizing Atmospheric Radiation – Part One « The Science of Doom[…] Part Five – The Code – code can be downloaded, includes some notes on each release […]

on February 17, 2013 at 11:37 am |Visualizing Atmospheric Radiation – Part Thirteen – Surface Emissivity « The Science of Doom[…] 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 ε = […]

on March 9, 2013 at 2:50 am |scienceofdoomI have an error in the Zhao diffusivity calculation. It’s a Matlab specific thing.

The code reads:

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

It should have a ‘./’ instead of ‘/’ before the (c3

Like this:

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

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.

on March 9, 2013 at 3:21 am |scienceofdoomHere is the equation from Zhao & Shi 2013 (click on the graphics to expand):

And their graph of the diffusivity parameter and % error for their paper vs the exact calculation:

And the result I get using the updated Matlab code:

on March 9, 2013 at 6:29 am |scienceofdoomAfter some more investigation I think that the (revised) algorithm is correctly reproducing Zhao & Shi’s.

For comparison I was calculating the exact integral of flux transmissivity and from that deriving the ‘effective’ value of ‘r’, from r = -ln(t)/τ

This is fine except at τ = 100, transmissivity, t = 10

^{-50}so rounding errors dominate in back-calculating r.And in any case the values are so small that even a large error would have negligible effect on overall radiative transfer calculations.

on March 9, 2013 at 7:23 am |Pekka PiriläSoD,

I state once more my reservations concerning the Zhao & Shi formula.

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.

on March 16, 2014 at 2:50 am |David HodgeDear Sod and Pekka,

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

on March 16, 2014 at 11:18 am |Pekka Pirilä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).

on March 16, 2014 at 1:06 pm |mrdodgePekka,

I would be delighted to help solidify the model. As SoD is traveling I think, we can wait until he is returned and able to comment

on April 4, 2014 at 4:00 am |scienceofdoomDavid,

Sorry for the delay in responding.

I had a check and I produced two minor updates, finishing at v0.10.6.

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.

on June 13, 2014 at 6:15 ammrdodgeah

for some reason i have only just seen this.

Thanks for the reply and please forgive my belated response.

Things have not progressed too much, as i have managed to change countries in the meantime , but I would certainly appreciate the latest copies if you could send them to me.

my challenge is getting all this to work in octave, as I dont have matlab – its nearly there i think.

And I love your work by the way, it one of the best educational experiences I have had in the last few years

on April 4, 2014 at 9:29 am |Pekka Pirilä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.

on June 13, 2014 at 8:05 pm |mrdodgeI will download Pekka’s model later today

Cheers!

on November 23, 2014 at 1:19 am |JuanHi 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.

on November 23, 2014 at 3:53 pm |Pekka PiriläJuan,

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

http://pirila.fi/energy/kuvat/HitranFilter.msi

but using that requires that the full HITRAN-database has been downloaded.

You might find the information that you search for also from SpetralCalc web pages.

on November 23, 2014 at 7:41 pm |DeWitt PayneJuan,

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.