Atmospheric Vortex Engine

Mathcad Calculations - Atmospheric Thermodynamics


Atmospheric Thermodynamic (AT_MCD) is Mathcad program for calculating:

  • The thermodynamic properties of air containing water in any phase,
  • The heat and work required to change the state of humid air masses,
  • The work produced when the position of air masses changes,
  • And for other related calculations.

The AT_MCD requires Mathcad version 5.0 or higher.

Mathcad is the leading programming language for mid-size technical calculations. Mathcad programs are easy to document because the equations appear in almost normal algebraic notation; and because results can be displayed or plotted right after they are calculated. Functions including solver functions can be defined at the beginning of the program and re-used whenever required. Mathcad has a powerful solver which is well suited for handling thermodynamic flashes where the number of component involved is small. Standard thermodynamic symbols have been used whenever possible to make the AT_MCD easier to follow.

More information on Mathcad is available at:


This section explains the basis for the calculations.

The first section of the program consists of standard functions which can be invoked later as required.

The standard function section includes:

  • Thermodynamic constants
  • Thermodynamic functions
  • Solver Blocks

Specific procedures usually follow the standard function section.

Calculating any thermodynamic function such as entropy, enthalpy, or virtual temperature is simply a matter of invoking the function.

The calculations are based on Thermodynamique de l'Atmosphère by L. Dufour and J. Van Meighem. The equations used to calculate thermodynamic properties are described in Appendix 1 of Michaud (1995).

Thermodynamic property are usually per unit mass of pure air. Three properties must be known to calculate any other thermodynamic properties of humid air. The program uses pressure (P) in kPa, temperature (T) in K, and total mixing ratio (M) in kg/kg, called the PKM arguments as standard arguments. Other humidity parameters such as: PCD, PCW, PCU, where pressure is in kPa, C is the temperature , D is the dew point, W is the wet bulb, all in degree C, and U is the relative humidity in fraction can be converted to the standard PKM format.

Units are usually in the base SI unit or an engineering multiple thereof. Pressures are in kPa. Mixing ratio is in kg/kg, relative humidity is in fraction to avoid complicating equations with multiples of 10. Displayed mixing ratio are multiplied by 1000 and displayed relative humidity are multiplied by 100 for readability. Mathcad is case sensitive, constant and variable names are usually in upper case. The functions defined at the beginning of the program start with a lower case "f".

The functions have been written so that any thermodynamic property can be calculated from these same three arguments always in the PKM order. The number 3 at the end of the function name is used as a reminder that the three standards properties are required as arguments, for example fST3(P,T,M) is the function to calculate the total entropy of 1 kilogram of air including the entropy of its water content in any phase. fHT3 is the enthalpy, fTV3 is the virtual temperature, and fU3 is the relative humidity. fMV3, fML3, and fMI3 are the moisture content in the vapour, liquid and solid phase. The letters T and M are used to designate property per unit mass of air plus its water content, and property per unit mass of substance respectively. For example fST3 is the entropy per unit mass of air, fSM3 is the entropy per unit mass of substance.

The mixing ratio (M) is the total water content in kg-water/kg-air, where the water can be in any phase. The calculations are based on the various phases being in equilibrium. If the mixing ratio is less than the saturation mixing ratio, the water is taken to be all in the vapour phase. If the mixing ratio is more than the saturation mixing ratio, the amount of water in the vapour phase is the saturation amount and the program calculates how much water there is in each phase.

The program is valid for 3 kinds of air:

1. Pure air where the moisture content is zero,

2. Moist air where the moisture is less than the saturated amount,

3. Saturated air where the water content is equal or greater than the saturation amount.

The three standard properties (PKM) are used for the three kinds of air. When the air is pure the third property (M) is zero. The phase rule requires that three properties be specified to describe humid air but limits the number of properties that can be specified to two at saturation. When the air is saturated, the third property is used to determine the quantity of condensed water. The program checks if M is beyond the saturation amount, if so the water beyond the saturated amount is taken to be in one or both of the condensed phases. Calculated properties include the contribution of the water in any phase.

The freezing temperature (TF) and the freezing band (FB) are defined in the global definitions. Any condensed water is taken to be in the liquid phase if the temperature is above the freezing temperature. If the temperature is below the bottom of the freezing band, the condensed water is taken to be all in the ice phase. The freezing point (TF) is the top of the freezing band. Within the freezing band the quantity of condensed water in each phase is taken to be proportional to the position in the band, at the middle of the freezing band the condensed water is half liquid and half solid. Using a freezing band eliminates a singularity when the condensed water freezes suddenly and corresponds to what happens in the atmosphere.

The effect of using a freezing band is equivalent to the transition phase suggested by Ooyama. A freezing point of -10 C with a freezing band of 20 C are used for most calculations. The effect of changing the freezing temperature and the freezing band can be checked by changing TF and FB. When there are three phases present the water in the gas phase is in equilibrium with the ice phase. The water in the liquid phase is not strictly in equilibrium but the error is negligible. Donner used a transition zone where the condensed water freezes linearly between 258 K and 248 K. The freezing temperature (TF) and the freezing band (FB) were placed at the beginning of the standard function section so that they can be changed if desired; nothing else in the standard function section should be changed. Freezing can be eliminated by setting TFC to -150 C. Eliminating the singularity when the condensed water freezes, requires that the freezing band (FB) be at least 5 K.

Duhem's rule states that for a closed system of known composition knowing any two thermodynamic properties is sufficient to calculate all others. The closed system considered consist of 1 kilogram of air and M kilograms of water, knowing the mass of water and any two properties all others thermodynamic properties can be calculated.

The program makes extensive use of the Mathcad solver. Mathcad requires that solver functions be defined before they are invoked. For this reason the standard are placed at the beginning of the program. Once defined with a GIVEN statement the solver functions can be used in the same way as global function.

The solver is used to calculate temperature during adiabatic expansion. fTSOL(S,P,M) is used to calculate isentropic expansion temperature or potential temperature. Temperature is calculated from entropy (S), pressure (P), and water content (M). Once temperature is known any other thermodynamic property can be calculated from the standard PKM arguments.

"Pseudo-Adiabatic" expansion where the water separates (open system) as it condenses is calculated by expanding in ~2 kPa steps and separating the condensed water after each step, and recalculating the entropy of the gas phase before the next step.

The solver approach is powerful.

1. There is no need to derive equations for entropy during adiabatic expansion. Entropy is always calculated from its fundamental definition. There is no need to calculate the latent heat for specific temperatures.

2. There is no need to use approximations when trying to re-arrange equations to isolate variables. The calculation process can always be reversed yielding back the initial conditions. The functions are sufficiently linear that the solver rarely has difficulty finding the solution. The initial guesses (PG, TG, and MG) can be the same for all cases.

3. There is no need to use different equations when going from moist air, to saturated air containing liquid water, to saturated air containing ice. There no need to make entropy adjustments when passing from the liquid to the ice phase.

4. Atmospheric processes tend to be of the type where some properties are conserved some are varied. For "true-adiabatic" expansion, the solver solves for temperature given entropy, pressure, and total mixing ratio. For mixing processes, the solver solves for temperature given enthalpy, pressure, and total water content.

5. The calculation of potential temperature is identical to the calculation of isentropic expansion temperature except that the final pressure is 100 kPa. Potential temperature thus includes the effect of water in any phase. There is no need for functions such as liquid water potential temperature, the potential temperature of a mass of air containing condensed water is the liquid water potential temperature provided the initial condition included the total water content. The use of potential temperatures is avoided because of ambiguity in their definitions make them more subject to misuse than standard thermodynamic properties. The program includes functions to calculate isentropic desiccation temperatures and equivalent temperatures and several other properties.

6. Thermodynamic properties calculated using the three PKM arguments are unambiguous. It is rarely clear whether potential temperature calculations include the effect of water in the vapour or liquid phase. The need to provide three arguments forces the user to carefully define the process.

The reference state, the zero entropy and the zero enthalpy state, for air is 0 C and 100 kPa. The reference state for water is liquid water at 0 C. The reference pressure for water vapour is the saturation vapour pressure of liquid water at 0 C, 0.61070 kPa.

AT1_MCD can be used to solve a wide variety of problems. The calculations in the articles by L. Michaud were done with Mathcad and checked with a similar HP48SX program. AT1_MCD is designed to be similar to the earlier Hewlett-Packard HP48SX calculator Atmospheric Thermodynamics program (AT1_HP48). AT1_HP48 and AT1_MCD use similar parameter and function names. The Mathcad version is faster and easier to document, but the HP calculator version has more powerful looping capability and the calculator version is more flexible and more convenient for testing new ideas. Mathcad has trouble with some arrays, experimentation was necessary to establish the Mathcad limitations.


Dufour, L., et J. Van Mieghem, 1975: Thermodynamique de l'Atmosphère.

Institut Royal Météorologique de Belgique, Bruxelles.

Ooyama, K.V., 1990: A thermodynamic foundation for modelling the moist atmosphere. J. Atmos. Sci., 47, 2580-2593.

Donner, L.J., 1993: A cumulus parameterization including mass fluxes, vertical momentum dynamics, and mesoscale effects. J. Atmos. Sci., 50, 889-906.

Michaud, L.M., 1995: Heat to work conversion during upward heat-convection. Part I: Carnot engine method. Atmos. Res., 39, 157-178.

Michaud, L.M., 1999: Vortex process for capturing mechanical energy during upward heat-convection in the atmosphere. Applied Energy, 62(4), 241-251.


Standard Mathcad functions

Standard functions with a numerical example of the function immediately following the function.

The standard function are used as the starting point in many of the atmospheric thermodynamic programs.
The Mathcad programs used in the Michaud articles will be added when time permits.



Vortex process for capturing mechanical energy during upward heat-convection in the atmosphere. Applied Energy, 1999, 62(4), 241-251.
Program used to calculate the work produced when air is raised reversibly from the surface to the 20 kPa level in tropical-oceanic areas. This program assumed a standard 20 kPa geopotential height and was used to generate Table 1 of the article. The calculations are based on condensed water not separating from the air and being allowed to freeze (true adiabatic with freezing).


Total Energy Equation Method for Calculating Hurricane Intensity.
Program used to calculate the minimum surface pressure when air is raised from the surface to the 20 or 10 kPa levels in tropical-oceanic areas. The minimum surface pressure is calculated for air approaching equilibrium at normal surface, and for air approaching equilibrium at the reduced eyewall pressure. The condensed water is not seperated from the rising air and is allowed to freeze. This program was used to generate Table 1 and Table 2 of the above paper.
Similar to the above except that a sounding is used in lieu of standard geo-potential heights.
Programs "MPI_S" can be used with any soundings. Holland (1997) used a Willis island average January sounding, but did not give the sounding. The above note used standard tropical geo-potential heights because the average sounding was not available. The Willis Island January 17, 1999, 0000Z sounding which is typical of January is used to illustrate the use of the program with a sounding. This program was used to generate Fig. 2 of the above paper.

Willis Island sounding WILLIS.prn
Sounding must be placed in the same directory as the Mathcad files.

Unedited complete WILLIS Sounding in Word format.

The sounding can be obtained at the University of Wyoming Atmospheric Sounding site.
The sample sounding is the 17 January, 1999, 0000Z Willis island sounding. Willis Island is North of East of Australia, station number 94299.


Michaud, L.M., 1998: Entrainment and detrainment required to explain updraft properties and work dissipation. Tellus A., 50A, 283-301.
Tellus Web Site. Abstract

Program to calculate Table 1 of the Article.Program:

Universal program to calculate the work done when air from the base of a specific sounding is raised with a specified entrainment, specified detrainment, and specified condensed water separtion is raised reversibly. Program:

Universal program to calculate the work done when air from the base of a specific sounding is raised with a specified entrainment, specified detrainment, and specified condensed water separtion is raised irreversibly. Program:

Sample Sounding, KAV92121800.prn. Kavieng 1992, 18 December, 0000Z sounding obtained from CODIAC and edited to 2 kPA interval .

Note: The following two recent articles showed that CAPE in tropical oceanic areas stays within the narrow 500 to 1800 J/kg range, and that the primary factor affecting the intensity of convection and precipitation is humidity aloft. These articles used a two dimensional model.
Program permits varying humidity aloft and can be used to produced result similar to those in the second of the two articles from a one dimensional model.

Lucas, C., and E. Zipser, 2000: Environmental variability during TOGA COARE. J. Atmos. Sci., 57, 2333-2350.
Lucas, C., E. Zipser, and B. Ferrier, 2000: Sensitivity of tropical west Pacific squall lines to tropospheric wind and moisture profiles. J. Atmos. Sci., 57, 2351-2373.

Mean TOGA COARE cluster 10 sounding. t10c0_2.prn provided by Christopher Lucas and edited to 2 kPa interval.
Program shows that for the above sounding, the mass flux peaks at 3 times the mass flux at base of the sounding at the 50 kPa level.
Increasing the humidity between between 96 and 20 kPa by 12%, increases the peak mass flux to 10 times the mass flux at the base of the sounding and increases the level to 42 kPa.
Increasing the relative humidity aloft by 12% increases the the mass flux by a factor of 3, and also increased cloud top height.