function [pr Tr rhor ztropo] = define_atmos_0_5(zr, Ts, ps, ptropo, lapset)
% Calculate pressure, pr; temperature, Tr; density, rhor; from height, zr;

% By assuming a standard temperature profile to a fixed tropopause temp
% v0.3 passes Ttropo from main routine (to allow a more fixed temperature
% structure if required)
% v0.4, Dec 31st 2012, includes prescribed lapse rate in K/m
% v0.5, Jan 1 2013, change to defining ptropo with isothermal above
%       temperature depends on Ts, lapset and ptropo - did some checks
%       and the stable soln is reached in 2 iterations


% Ts=surface temperature
% ps=surface pressure
% zr is a linear vector and starts at the surface
% ptropo = pressure of tropopause 
% lapset = lapse rate in K/m

pr=zeros(1,length(zr)); % allocate space
Tr=zeros(1,length(zr)); % allocate space
rhor=zeros(1,length(zr)); % allocate space

% first calculate temperature profile
% *** htropo=10000; % height of tropopause
Ttropo=215; % *initial* tropopause temperature to allow pressure calculation

num=3; % number of guesses
ztropo=((Ts-Ttropo)/lapset); % estimate height of tropopause in m
for n=1:num
    for i=1:length(zr)
        if zr(i)<ztropo
            Tr(i)=Ts-(zr(i)*lapset);  %  the troposphere
        else
            Tr(i)=Ttropo; % isothermal above tropopause
        end
    end
    
%     subplot(2,2,n), plot(Tr,zr/1000)
%     grid on
    
    % pressure, p = ps * exp(-mg/R*integral(1/T)dz) from 0-z
    mr=28.97e-3; % molar mass of air
    R=8.314; % gas constant
    g=9.81; % gravity constant
    intzt=0; % sum the integral of dz/T in each iteration
    dzr=zr(2)-zr(1); % dzr is linear so delta z is constant
    pr(1)=ps; % surface condition
    
    % calculate pr as a function of z for the initial temperature guess
    for i=2:length(zr)
        intzt=intzt+(dzr/Tr(i));
        pr(i)=ps*exp(-mr*g*intzt/R);
    end
    
    if n<num
        % find the refined height & temperature of the tropopause
        ztropo=zr(find(pr<=ptropo,1)); 
        Ttropo=Ts-ztropo*lapset;
        %disp(['Num= ' num2str(n) ', ztropo= ' num2str(ztropo) ', Ttropo= ' num2str(Ttropo)]);
    end
    
end

% density, rhor = mr.p/RT
rhor=mr.*pr./(R.*Tr);
end

