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)