%irbroad.m
% Written by:  Shaima Nasiri & S. Ackerman
%              29 March, 1999%
%Broadband flux emissivity code
% Given an atmospheric gas (similar to water), whose broadband 
% emissivity flux is approximated as E = 1.039+0.4246*
% log10(00033 + 0.1155*U), 
% Upward and downward fluxes at every pressure level are calculated.
%
% level(1) is top of atmosphere
% level(numlev) is surface

% Variables:
% by level
%        P = Pressure  (mbar)     
%        T = Temperature Kelvin) 
%        q = specific humidity   
%    Fdown = downward flux
%      Fup = upward flux
% by layer
%       Ta = average temperature
%       Pa = average pressure
%        U = precipitable water (g/cm^2)
%       Us = pressure corrected precipitable water
%        E = broadband flux emissivity
%   deltaE = layer difference of E
% PROFILE% load mid-latitude summer standard atmospheric profile
% mls = Mid-Latitude Summer Atmosphere          
% Height Pressure       Temp   H2O Density     O3 Density
% km       mbar           K     g/m^3            g/m^3
 load mls.txt

 R = 287;                          %gas constant of dry air  {J /(K Kg})

 [maxlev,tmp] = size(mls);         %number of levels in profile
 uselev = size(mls);               %number of levels use
                                   %tropopause around level 12

 P = flipud(mls(1:uselev,2));
 T = flipud(mls(1:uselev,3));
 rho = flipud(mls(1:uselev,4));    %H2O density
 q = rho./(P.*1e2./R./T);
 [numlev temp] = size(P);          % numlev = number of pressure levels 
 numlay = numlev-1;                % numlay = number of atmospheric layers
                                 % P(i+1) - P(i) pressure interval in mbar
 deltaP(1:numlay) = P(2:numlev)-P(1:numlay) ; 


% Constants 
 g = 9.81;                         % gravity in m/s^2 
 sigma = 5.67e-8;                  % W/m^2/K^4 
 conv = 100/100^2;                 % conversion for u to g/cm^2
 m = 0.85;                    % exponent in pressure correction for H2O 

% Calculating values for each layer
 for i=1:numlay    
    Ta(i) = (T(i) + T(i+1))/2;  
    Pa(i) = (P(i) + P(i+1))/2;   
    U(i) = 1/g*(q(i+1)+q(i))/2*deltaP(i)*conv;
    Us(i) = U(i) * (Pa(i)/P(numlev))^m * (T(numlev)/Ta(i))^0.5;  
end
% Initializing flux arrays to zero 
Fdown = zeros(size(P)); 
Fup   = Fdown;
% Calculation of Fdown at pressure level P(count) 
for count=1:numlay
   for i=1:count
      E(i) = 1.039 + 0.4246 * log10(0.0033 + 0.1155*sum(Us(i:count)));
      if E(i)<0;
         E(i)=0.;
         end
   end
   for i=1:count
      if i==count
         deltaE(i) = E(i);
      else
         deltaE(i) = E(i) - E(i+1);
      end
   end
   Fdown(count+1) = sum(sigma*Ta(1:count).^4.*deltaE(1:count));
 end

% Calculation of Fup at pressure level P(count) 
for count=numlay:-1:1
   for i=numlay:-1:count
      E(i) = 1.039 + 0.4246 * log10(0.0033 + 0.1155*sum(Us(i:-1:count)));
            if E(i)<0;
         E(i)=0.;
         end

   end
   for i=numlay:-1:count
      if i==count
         deltaE(i) = E(i);
      else
         deltaE(i) = E(i) - E(i-1);
      end      
   end
    Fup(count) = sigma*T(numlev).^4*(1-sum(deltaE(numlay:-1:count))) + ...
      sum(sigma*Ta(numlay:-1:count).^4.*deltaE(numlay:-1:count));
end
Fup(numlev)=sigma*T(numlev).^4;
plot(Fdown,P,Fup,P,Fdown-Fup,P)
set(gca,'ydir','reverse')
xlabel('Flux (W m^{-2})');
ylabel('Pressure (mb)');
legend ('F\downarrow','F\uparrow','F\downarrow-F\uparrow');

%  class
% calculate heating profile







