#rh.cgi by Tim Brice NWS El Paso


#Constants
$Rd = 287.04;
$cp = 1005;

#Sub-routines
sub ctof {
  local($C) = @_;
  $F = (9/5 * $C) + 32;
  return ($F);
}
sub ftoc {
  local($F) = @_;
  $C = 5/9 * ($F - 32);
  return ($C);
}
sub ctok {
  local($C) = @_;
  $K = $C + 273.15;
  return ($K);
}
sub vappress {
  local($ktemp) = @_;
  $expo = 23.832241 - 5.02808*(log($ktemp)/log(10)) - 1.3816*10**-7 *
10**(11.334 - .0303998 * $ktemp);
  $expo = $expo + 8.1328 * 10**-3 * 10**(3.49149 - 1302.8844/$ktemp) -
2949.076/$ktemp;
  $E = 10**($expo);
  return ($E);
}

sub thetae {
  local($ktemp,$press,$mixr1,$klcl) = @_;
  $thetae1 = $ktemp * (1000/$press)**(($Rd/$cp)*(1 - .28 * $mixr1)) *
exp((3376/$klcl - 2.54) * $mixr1 * (1 + .81 * $mixr1));
  return ($thetae1);
}


#change entries to celsius and millibars
if ($rh{airtemp} =~ /-?(\d+)\.?(\d+)|-?(\d+)/)  {
  if ($rh{corf} eq Fahrenheit) {
   &ftoc($rh{airtemp});
   $temp = $C;
} else {
   $temp = $rh{airtemp};
  }
} else {
   print "An illegal character was entered.\n";
   die
}

if ($rh{press} =~ /(\d+)\.?(\d+)|-?(\d+)/)  {
 if ($rh{mborin} eq millibars) {
   $press = $rh{press};
} else {
   $press = $rh{press} * 33.8639;
 }
} else {
 print "An illegal character was entered.\n";
 die
}

if ($rh{rh} =~ /(\d+)\.?(\d+)|-?(\d+)/)  {
  $rh = $rh{rh};
} else {
 print "An illegal character was entered.\n";
 die
}

#Start the math

&ctok($temp);
$ktemp = $K;

&vappress($ktemp);
$etemp = $E;

$edewp = $etemp * ($rh/100);

#start guessing for the Dewpoint
$estdewp = $temp;
$incr = 5;
$prevsign = 1;
$este = 100;

while (abs($edewp - $este) > .01) {
 
  &ctok($estdewp);
  $kestdewp = $K;

  &vappress($kestdewp);
  $este = $E;
  $cursign = $este <=> $edewp;

  if ($prevsign != $cursign) {
    $prevsign = $cursign;
    $incr = $incr * .5;
}
  $estdewp = $estdewp - $incr * $cursign;
 
}
$dewp = $estdewp;
 
#calculate the TLCL and PLCL
$tlcl = $dewp - (.212 + .001571 * $dewp - .000436 * $temp) * ($temp -
$dewp);

&ctok($tlcl);
$klcl = $K;

&ctok($temp);
$ktemp = $K;

$plcl = $press * ($klcl/$ktemp)**(1/($Rd/$cp));

#calculate the vapor press and mixing ratio for the surface parcel
&ctok($dewp);
$kdewp = $K;

&vappress($ktemp);
$esat = $E;

&vappress($kdewp);
$esurf = $E;

$mixr = 621.97 * $esurf/($press - $esurf);

#calculate the potential and equivalent potential temperatures
#for the surface parcel
$theta = $ktemp * (1000/$press)**($Rd/$cp);

$mixr1 = $mixr/1000;
&thetae($ktemp,$press,$mixr1,$klcl);
$thetae = $thetae1;

#let's start the guessing game
$tguess = 0;
$previsign = 1;
$incr = 10;
$thetaediff = 1;
$estthetae = 1;

while (abs($thetaediff) > .05) {
 &ctok($tguess);
 $ktguess = $K;

 $esttheta = $ktguess * (1000/$press)**($Rd/$cp);

 &vappress($ktguess);
 $eguess = $E;

 $mixrguess = 621.97 * $eguess/($press - $eguess);

 $mixrguess1 = $mixrguess/1000;
 &thetae($ktguess,$press,$mixrguess1,$klcl);
 $estthetae = $thetae1;

 $thetaediff = $thetae - $estthetae;
 $cursign = $thetae <=> $estthetae;

  if ($cursign != $prevsign) {
     $prevsign = $cursign;
     $incr = $incr/2;
  }
$tguess = $tguess + $incr * $prevsign;
}



#let's make some pretty output
$~= DP8;
write;

format DP8 =
With an air temperature of @###.# degrees @<<<<<<<<<< a relative
                           $rh{airtemp}   $rh{corf}
humidity of @##.# percent and a station pressure of @###.## @<<<<<<<<<<<<<<<<<:
            $rh                                     $rh{press} $rh{mborin}
.



$~ = DP1;
write;

format DP1 =
You get a wet-bulb temperature of
.

if ($rh{corf} eq Fahrenheit) {
  &ctof($tguess);
  $tguessf = $F;
  $~ = DP2;
  write;
} else {
  $~ = DP3;
  write;
}

format DP2 =
 @###.# degrees Fahrenheit
 $tguessf
.

format DP3 = 
 @###.# degrees Celsius
 $tguess
.

if ($rh{corf} eq Fahrenheit) {
  &ctof($dewp);
  $dewpf = $F;
  $~ = DP4;
  write;
} else {
  $~ = DP5;
  write;
}

format DP4 =
  and a dewpoint temperature of @###.# degrees Fahrenheit.
                                   $dewpf
.

format DP5 = 
  and a dewpoint temperature of @###.# degrees Celsius.
                                   $dewp
.

if ($rh{rh} > 100) {
  print "
The relative humidity is greater than 100 percent, 
  this only rarely occurs in the atmosphere.  You may want to 
  go back and double check your entries.
";
}