Code:
#include <stdio.h>
#include <math.h>
int main()
{
float P, h, T, DpT; // P=pressure, h=height, T=temperature, DpT=dewpoint temp
int i; // counter for loop
int j; // counter for loop
float K, km; // K=Temp in Kelvin, km=height in kilometers, LR=Lapse Rate
float esw; // saturation vapor pressure over water
float esi; // saturation vapor pressure over ice
float height[500]; // array of heights in kilometers
float tempC[500]; // array of temp in celcius
float tempK[500]; // array of temp in kelvin
float pressure[500]; // array of pressures
float dewpt[500]; // array of dewpoint temps
float BP = 373; // Boil point at 1 atm
float TP = 273.16; // Triple point
float estw = 1013.25; // saturation vapor pressure at the steam-point pressure
float esi0 = 6.1173; // saturation vapor pressure at the ice-point pressure
float part1, part2, part3, part4, part5, part6; // parts of Goff-Gratch formula
const int header_length = 0;
char buff[256];
j = 0;
FILE *input;
input = fopen("KLM_ama.txt","r"); //input data file
if(input == NULL) //input check
{
printf("check your file name!\n");
return(0);
}
for(i=0; i<header_length; i++)
{
fgets(buff,265,input);
}
while(fgets(buff,256,input))
{
fscanf(input,"%f %f %f %f", &P,&h,&T,&DpT);
tempC[j] = T;
pressure[j] = P;
dewpt[j] = DpT;
K = tempC[j] + 273; // Convert to Kelvin
km = h/1000; // Convert height to kilometer
tempK[j] = K;
height[j] = km;
// printf("%5.0f %5.2f %f \n",pressure[j],height[j],tempK[j]);
j++;
}
for(j=0; j<200; j++)
{
if(tempK[j]>223 && tempK[j]<375) // calculation over water
{
part1 = -7.90298*((BP/tempK[j])-1)+5.02808*log10(BP/tempK[j]);
part2 = 0.00000013816*((pow(10,11.344*(1-(tempK[j]/BP))))-1);
part3 = 0.0081328*((pow(10,-3.49149*((BP/tempK[j])-1)))-1);
esw = (pow(10,part1-part2+part3+log10(estw)));
// printf("Saturation vapor pressure over water is %5.2f hPa at %5.0f mb \n", esw, pressure[j]);
}
if(tempK[j]<273 && tempK[j]>173) // calculation over ice
{
part4 = -9.09718*((TP/tempK[j])-1);
part5 = 3.56654*log10(TP/tempK[j]);
part6 = .876793*(1-(tempK[j]/TP));
esi = (pow(10,part4-part5+part6+log10(esi0)));
// printf("Saturation vapor pressure over ice is %5.3f hPa at %5.0f mb \n", esi, pressure[j]);
}
}
for(j=0; j<1; j++) // calculation of LCL
{
float LCL;
LCL = ((tempC[j]-dewpt[j])/2.5)*304.8;
printf("LCL hieght is %5.2f meters\n", LCL);
}
for(j=0; j<200; j++) // calculation of freezing levels
{
if(tempC[j]<0)
{
printf("freezing levels at: %5.2f hpa\n", pressure[j]);
}
}
}