function [ vt tau stau fluxu fluxd ztropo z dz p rho mixh2o zb pb Tb rhob ...
    Tinit T radu radd rads radu1 emitu TOAf TOAtr] =...
    HITRAN_0_9_2(dv, numz, mix, blp, BLH, FTH, contabs, Ts, lapsekm, lapseikm, ptropo, minp, tstep, nt, 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

% dv =  resolution of wavenumber for calculations
% numz =  number of boundaries to consider (number of layers = numz-1)
% mix = mixing ratio for each molecule, but first value H2O unused
% blz = boundary layer thickness in m
% RLH =  relative humidity in the boundary layer
% FTH =  relative humidity in the free troposphere (currently layers 2,3..)
% contabs = true or false for including the water vapor continuum absorption
% Ts =  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)
% tstep = timestep in seconds 
% ploton = whether the plot is run or not - lots of other plots are
%         generated via separate programs

disp(['--- New Run v0.9.2 --- ' datestr(now) ' ----']);
disp('  ');

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

% Units in SI unless otherwise noted, apart from 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
diff=1.66; % diffusivity approximation
ems=1.0; % emissivity of surface
cp=1005; % cp (specific heat cap at constant pressure) of air
convadj=true; % if true, convective adjustment to lapse rate = lapse
iso_strato=true; % if true, fixed and isothermal stratospheric temperature
lapse=lapsekm/1000; % lapse rate in K/m
lapsei=lapseikm/1000; % initial lapse rate in K/m 
lapsez=ones(1,numz-1)*lapse; % lapse rate for convective adjustment can be
                 % different for each layer, current all set to 'lapse'
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
if contabs
    conttext='Continuum included';
else
    conttext='No Continuum';
end

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

vmin=200; % 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.. '); 
% define GHGs to use, isotopes to consider, define isotope ratios
mol=[1 2 3 4 6]; % these numbers correspond to molecule numbers in HITRAN list
nmol=length(mol);
mols=['H2O'; 'CO2'; 'O3 '; 'N2O'; 'CH4';]; % names for the molecules in 'mol'
ison=[1 1 1 1 1]; % number of isotopes for each GHG (1 = only the main is used)
isoprop=[0.997 0; 0.984 0; 1.0 0; 0.990 0; 0.988 0]; % isotope proportions, values needed 
            % only to the value of 'ison' for each GHG

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

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

ps=1.013e5; % define surface pressure

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,Ts,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 approximately equal pressure changes
% except for first layer which is the boundary layer at pressure = blp 
% with different humidity (BLH)

% this part of the code was fixed on Jan 9, 2013 for v0.9.2, previously 
% (and up to article 4 in the "Visualization" series) created an arbitrary
% thickness layer above boundary layer due to poor code
dp=(blp-minp)/(numz-2); % finds the pressure change for each height change

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

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

for i=1:numz-1 
    % now calc mixing ratio in molecules for water vapor at prescribed RH
    % = RH * es / p
    % currently satvaph2o() gives es=610.94.*exp(17.625.*Tc./(Tc+243.04)); where Tc in 'C
    if i==1
        mixh2o(i)=BLH*satvaph2o(Tinit(i))/p(i); 
    else
        mixh2o(i)=FTH*satvaph2o(Tinit(i))/p(i);
    end
end
   
if print_layers        % display values at each pressure
    disp(['dp= ' num2str(dp)]);
    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)

% now set optical thickness, [layers x numv]
fprintf('%s', 'Calc optical thickness each layer.. '); 
tau=zeros(numz-1,numv); % optical thickness 
              
% ----------- 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
            mix(m)=mixh2o(i); % precalculated for each layer using prescribed RH
            % disp(['w.v. mixing ratio at ' num2str(p(i)/1e2) ' mbar is ' num2str(mix(m)*1e3)]);
            if contabs % if using continuum absorption
                % method from Pierrehumbert 2010, p.260, uses svp & units
                % in m2/kg
                %disp(['Continuum for layer ' num2str(i)]);
                es=satvaph2o(Tinit(i)); % get sat.vap. pressure in Pa
                %disp(['Sat vap pressure = ' num2str(es)]);
                sig=continuum1000(vt,Tinit(i)); % get capture cross section in m2/kg 
                % ** note non-std units ** and values are only valid for 
                % range defined in function continuum1000                
                if i==1
                    RH=BLH; % use the correct relative humidity
                else
                    RH=FTH; % use the correct relative humidity
                end
                % This value of optical thickness will be most sensitive
                % to temperature changes when time steps take place
                % - so later need to consider having this recalculated after 
                % set dT'C change in temperature in a layer
                tau(i,:)=tau(i,:)+sig.*d/100*(es*RH)^2/(1e4*461.4.*Tinit(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)*na;
            % then calc the actual line width for this temp & pressure
            ga=gama((im(j))).*(p(i)/p0).*(T0/Tinit(i)).^nair((im(j)));
            
            % note that if we adjust temperature considerably over timesteps 
            % then may need an iteration of complete routine because line 
            % widths change with temperature
            
            % -------- Iterate through each wavenumber, k, for this one
            %          line -----------
            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
disp(' '); % to produce a newline

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

rads=ems.*planckmv(vt,Ts); % surface spectral emissive power vs wavenumber, vt
% *** this assumes surface temperature is not changing
T=zeros(nt,numz-1); % define array to store T for each level and time step
T(1,:)=Tinit;   % load temperature profile for start of scenario
    
% 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
    
    
% =========== 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
    %fprintf('%d %s', h, ' '); % display current timestep being calculated
    radu=zeros(numz,numv); % initialize upward intensity at each *boundary* and wavenumber
    radu1=zeros(numz,numv); % initialize different version of radu
    emitu=zeros(numz,numv); % initialize emission upward for each *layer* and wavenumber
    
    radd=zeros(numz,numv); % initialize downward intensity at each *boundary* and wavenumber
    emitd=zeros(numz,numv); % initialize emission downward for each *layer* and wavenumber
    radu(1,:)=rads;  % upward *surface* radiation vs wavenumber
    radu1(1,:)=rads; % evaluating radiation transfer with no emission - for interest
    radd(end,:)=zeros(1,numv);  % downward radiation at TOA vs wavenumber (=zero)
    % units of radu, radd are W/m^2.cm^-1, i.e., spectral emissive power
    
    % Upwards spectral emissive power 
    Eabs=zeros(numz-1); % zero the absorbed energy before we start
    for i=1:numz-1   % each layer
        for j=1:numv   % each wavenumber interval
            
            % calc change in monochromatic wavenumber at next boundary
            % = incoming * transmission + 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(h-1,i))-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 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
            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(h-1,i))-1); 
            % calculate spectral emissive power at next upward boundary
            radd(i,j)=radd(i+1,j)*tran(i,j)+emitd(i,j);
            Eabs(i)=Eabs(i)+(radd(i+1,j)-radd(i,j))*dv;
        end  % each wavenumber interval
        
        dT=Eabs(i)*tstep/(cp*rho(i)*dz(i)); % change in temperature = dQ/heat capacity
        if iso_strato && z(i)>=ztropo  % if we want to keep stratosphere at fixed temperature
            T(h,i)=T(h-1,i);
        else
            T(h,i)=T(h-1,i)+dT; % calculate new temperature
            
            if T(h,i)>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 routine is so fast 
    end  % each layer
    
    % now apply convective adjustment
    if convadj==true  % if convective adjustment chosen..
        % First from surface to layer 1, slightly different code
        if (Ts-T(h,1))/z(1)>lapsez(1)
            T(h,1)=Ts-z(1)*lapsez(1);
        end
        % Now all the atmospheric layers
        for i=2:numz-1            
            if (T(h,i-1)-T(h,i))/(z(i)-z(i-1))>lapsez(i)  % too cold, convection will readjust
                % adjust temperature to lapse rate
                T(h,i)=T(h,i-1)-((z(i)-z(i-1))*lapsez(i));                
            end
        end
    end
    
end   % ---- timestep iterations

% =========== Finish loops to calculate TOA spectrum & flux =====
disp(' '); % to produce a newline

fluxu=sum(radu(end,:))*dv;   % calculate the TOA flux
fluxd=sum(radd(1,:))*dv;   % calculate the DLR total

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

% flux transmitted from surface and in "atmospheric window"
% calculate the surface radiation transmitted to TOA (integral with dv)
TOAtr=dv*rads*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
%TOAf(numz-2,:)=emitu(numz-2,:).*stau(numz-1,:); this gets done inside the
%loop below
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*rads.*exp(-stau(1,:)); % rads is like emitu(0) - the layer below the bottom layer of atmosphere
%sanity check
disp(['Value of TOA flux via curiosity method = ' num2str(sum(sum(TOAf)))]);
figure(2)
zz=[0 z]; % add the surface to the height vector to make it possible to plot
plot(sum(TOAf'),zz/1000,'r*')
xlabel('W/m^2 contribution to TOA flux')
ylabel('Height of layer, km')
grid on
    
disp(['Ts= ' num2str(Ts) ', 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,   ' num2str(nt) ...
    ' timesteps of ' num2str(tstep/3600/24) ' days']);
disp(['CO2= ' num2str(mix(2)*1e6) 'ppm, O3= ' num2str(mix(3)*1e9)...
    'ppbv, N2O= ' num2str(mix(4)*1e9) 'ppbv, CH4= ' num2str(mix(5)*1e9) 'ppbv']);
disp(' ');
disp(['Flux up at TOA = ' num2str(fluxu,'%4.0f') ' W/m^2, Flux down at surface = ' ...
    num2str(fluxd,'%4.0f') ' W/m^2']);
disp(['TOA transmitted flux= ' num2str(TOAtr,'%4.0f') ' W/m^2, Surface emission= ' ...
    num2str(sum(rads),'%4.0f') ' W/m^2, Window transmission= ' ...
    num2str(sum(radu1(numz,iw))/sum(rads(iw))*100,'%4.0f') ' %']);
  
% =========== Plots, if needed ==================
if ploton
    figure(1)
    % first graph - temperatures start and end
    subplot(2,2,1),plot(Tinit,z/1000,T(end,:),z/1000)
    %set(gca,'YDir','reverse')
    grid on
    xlabel('Temperature, K')
    ylabel('Height, km')
    %ylabel('Pressure, hPa')
    title(['Temperature profile, initial Ts=' num2str(Ts),'K, Lapse=' num2str(lapsekm)])
    legend('Initial','Final')
    
    zradu=sum(radu');zradd=sum(radd');
    % second graph - final flux values
    subplot(2,2,2),plot(zradu(2:end),z/1000,zradd(1:end-1),z/1000)
    %set(gca,'YDir','reverse')
    grid on
    xlabel('Flux, W/m^2')
    ylabel('Height, km')
    %ylabel('Pressure, hPa')
    title(['Up & Down Flux, BLH=' num2str(BLH*100) '% FTH=' num2str(FTH*100) ...
        '%, TOA=' num2str(fluxu,'%4.0f') 'W/m^2 DLR=' num2str(fluxd,'%4.0f') 'W/m^2'])
    legend('Up','Down')
    
    % third graph - spectral up
    subplot(2,2,3),plot(vt,radu(end,:))
    xlabel('Wavenumber, cm^-^1')
    ylabel('Spectral intensity')
    title('TOA Upward Spectrum')
    grid on
    
    % fourth graph - spectral down
    subplot(2,2,4),plot(vt,radd(1,:))
    xlabel('Wavenumber, cm^-^1')
    ylabel('Spectral intensity')
    title('Surface Downward Spectrum')
    grid on
end
% ==================================================

disp(['---- Complete End ---- ' datestr(now) ' ----']);
