% Phys189L_humidity_Solution (from Dr. Minschwaner)
% Template matlab script file for Phys189, 2017
% This script reads in two sounding data files and plots RH for each
%
% Here is the name of the ascii text file that contains Socorro sounding
%fname='Airsonde_sounding_23_Oct_2017.txt' ;
fname='Airsonde_sounding_23_Oct_2017-noheader.txt' ;
% This command reads the data into the array "a". The data delimiter is four blank spaces
%a = dlmread(fname,' ',1,0); 
a = dlmread(fname) ; 
% Get the number of rows and columns in the array "a", in case we need it
nrow=size(a,1) ; 
ncol=size(a,2) ;
% These commands split the columns of "a" into two vectors, "p1" and "t1"
% where "p1"=pressure in mb, "t1"=temperature in deg C, "wb"=wet-bulb temp
  p1 = a(:,4) ; 
t1 = a(:,2) ;
wb = a(:,3) ;
% Compute altitudes using barometric law, assume H=scale height = 7.5 km
h=7.5 ;
% Also assume surface is 1.42 km above mean sea level
zsurface=1.42 ;
z1=zsurface + h*log(p1(1)./p1) ;
% Relative humidity based on temperature, wet-bulb temp, and pressure
% involves a number of steps
% First, we calculate the SATURATION water vapor pressure at each level
ps=6.108*exp((17.27.*t1)./(237.0+t1)) ;
% Next, get the SATURATION water vapor pressure at the wet bulb temperature
pswb=6.108*exp((17.27.*wb)./(237.0+wb)) ;
% Now we calculate the ENVIRONMENTAL water vapor pressure at each level
pse=pswb + (-0.00066 - 7.59e-7.*wb).*p1.*(t1-wb) ;
% Finally, relative humidity is the ration of ENVIRONMENTAL to SATURATION vapor pressure
rh1=pse./ps ;
% Multiply by 100 to convert to percent
rh1=100.*rh1 ;
% Here is the name of the ascii text file that contains albuquerque sounding
%fname='KABQ_12Z_23_Oct_2017.txt' ;
fname='KABQ_12Z_23_Oct_2017_noheaderfooter.txt' ;
% This command reads the data into the array "b". The data delimiter is two blank spaces
%b = dlmread(fname,' ',[8 0 146 10]) ;
b = dlmread(fname) ; 
% Get the number of rows and columns in the array "b", in case we need it
nrow=size(b,1) ;
ncol=size(b,2) ;
% These commands split the columns of "b" into three vectors, "z2" and "t2" and "rh2"
% where "z2"=altitude in meters, "t2"=temperature in deg C, "rh2"=relative humidity in %
z2 = b(:,2) ;
t2 = b(:,3) ;
rh2 = b(:,5) ;
% Convert meters to km
z2=1.e-3*z2 ;
% Plot the socorro and albuquerque relative humidity data
plot(rh1,z1,'b-',rh2,z2,'r:') ;
axis([0 100 1 11]) ;
xlabel('Relative Humidity (%)') ;
ylabel('Altitude (km)') ;
title('Humidity Soundings Oct 23, 2017') ;
legend('Socorro 21Z','Albuquerque 12Z') ;
%
