Added potential_temperature function

This commit is contained in:
Daniele Tricoli 2011-01-15 03:42:38 +01:00
parent bd47fb6f18
commit b3dba6d4e5
3 changed files with 36 additions and 0 deletions

View file

@ -192,3 +192,25 @@ double adiabatic_temperature_gradient(double salinity, double temperature,
temperature + 1.8932e-6) * salinity + ((6.6228e-10 * temperature -
6.836e-8) * temperature + 8.5258e-6) * temperature + 3.5803e-5;
}
double potential_temperature(double salinity, double temperature,
double pressure, double reference_pressure)
{
double h, xk, q;
h = reference_pressure - pressure;
xk = h * adiabatic_temperature_gradient(salinity, temperature, pressure);
temperature += 0.5 * xk;
q = xk;
pressure += 0.5 * h;
xk = h * adiabatic_temperature_gradient(salinity, temperature, pressure);
temperature += 0.29289322 * (xk - q);
q = 0.58578644 * xk + 0.121320344 * q;
xk = h * adiabatic_temperature_gradient(salinity, temperature, pressure);
temperature += 1.707106781 * (xk - q);
q = 3.414213562 * xk - 4.121320344 * q;
pressure += 0.5 * h;
xk = h * adiabatic_temperature_gradient(salinity, temperature, pressure);
return temperature + (xk - 2.0 * q) / 0.6;
}

View file

@ -10,5 +10,7 @@ double freezing_point(double salinity, double pressure);
double specific_heat(double salinity, double temperature, double pressure);
double adiabatic_temperature_gradient(double salinity, double temperature,
double pressure);
double potential_temperature(double salinity, double temperature,
double pressure, double reference_pressure);
#endif

View file

@ -93,6 +93,17 @@ START_TEST(test_adiabatic_temperature_gradient)
}
END_TEST
START_TEST(test_potential_temperature)
{
ck_assert(cmp_double(potential_temperature(25, 0, 0, 0), 0.000000));
ck_assert(cmp_double(potential_temperature(25, 40, 0, 0), 40.000000));
ck_assert(cmp_double(potential_temperature(30, 20, 9000, 0),
18.296537));
ck_assert(cmp_double(potential_temperature(40, 40, 10000, 0.0),
36.998992));
}
END_TEST
Suite *oceanography_suite(void)
{
@ -105,6 +116,7 @@ Suite *oceanography_suite(void)
tcase_add_test(tc_core, test_freezing_point);
tcase_add_test(tc_core, test_specific_heat);
tcase_add_test(tc_core, test_adiabatic_temperature_gradient);
tcase_add_test(tc_core, test_potential_temperature);
suite_add_tcase(s, tc_core);
return s;