Fortran - Re: physical unit/dimension checking and conversion in scientific

Author

Re: physical unit/dimension checking and conversion in scientific

Grant W. Petty

2006-01-24, 7:02 pm

On Sun, 22 Jan 2006, N. Shamsundar wrote:

> Despite the Fortran 77 style of coding used in your example, this solution is
> beyond the reach of standard Fortran 77 compilers (e.g., g77) since the
> solution is based on using program units and USE statements.

That's correct. My only reason for supplying F77 style examples was to
illustrate how easily one can take legacy F77 code and adapt it to use the
PHYSUNITS package. Most of the scientists in my field still use F77. But
yes, you have to compile it with F90.

> about : units that differ not only by a factor but also by an offset, such as
> temperatures, elevations, etc. For example,
>
> type(preal) Tamb
> Tamb=77*u_fahrenheit_degree
> write(*,*)' Atm Temp = ',(Tamb-c_ice_point)/u_kelvin,' C'
>
> does not work for obvious reasons, but no error message is produced.

That's correct. As a meteorologist, I'm painfully aware of the problems
posed by temperature conversions. Unfortunately, there's no simple
solution. The commercial MathCad program, which inspired my foray into
this problem, doesn't have a solution either.

If you understand how things work "under the hood", you can deal with a
temperature conversions in a straightforward way, but if you don't, you can
easily mess up. User beware!

In the example you gave, I would simply write

type(preal) Tamb
Tamb=77*u_fahrenheit_degree + c_ice_point
write(*,*)' Atm Temp = ',(Tamb-c_ice_point)/u_kelvin,' C'

>
> Another common source of error is exemplified by the Rayleigh number:
>
> Ra = g\beta\Delta T L^3/(\alpha\nu)
>
> where
>
> g is acceleration due to gravity
> \beta is thermal expansion coefficient, units of reciprocal abs. temperature.
> \Delta T is temperature difference, same in kelvin or celsius.
> \alph, \nu are diffusivities, units of L^2/time
> L is length
>
> How would you handle this formula in your system?

Not sure where the complication is in this case. Am I overlooking something
subtle here? I assume Ra comes out dimensionless, in which case Ra could
be declared as ordinary REAL.

Just looking at what you wrote, I'd simply assign values with suitable
units to each of the parameters/variables (declared as type PREAL) and then
compute Ra as

Ra = real(g*beta*DeltaT*(L**3)/(alpha*nu))

If the argument to real() computes to something other than a
non-dimensional number, it will generate a run-time error.

- Grant

________________________________________

Prof. Grant W. Petty
Atmospheric and Oceanic Sciences
1225 W. Dayton Street
University of Wisconsin-Madison
Madison, WI 53706

gpetty@aos.wisc.edu
Tel: (608) 263-3265
Fax: (608) 262-0166

N. Shamsundar

2006-01-24, 9:58 pm

Grant W. Petty wrote:
> On Sun, 22 Jan 2006, N. Shamsundar wrote:
>
>
>
> That's correct. My only reason for supplying F77 style examples was to
> illustrate how easily one can take legacy F77 code and adapt it to use
> the PHYSUNITS package. Most of the scientists in my field still use
> F77. But yes, you have to compile it with F90.
>
>
>
> That's correct. As a meteorologist, I'm painfully aware of the problems
> posed by temperature conversions. Unfortunately, there's no simple
> solution. The commercial MathCad program, which inspired my foray into
> this problem, doesn't have a solution either.
>
> If you understand how things work "under the hood", you can deal with a
> temperature conversions in a straightforward way, but if you don't, you
> can easily mess up. User beware!
>
> In the example you gave, I would simply write
>
> type(preal) Tamb
> Tamb=77*u_fahrenheit_degree + c_ice_point
> write(*,*)' Atm Temp = ',(Tamb-c_ice_point)/u_kelvin,' C'
>
>
>
> Not sure where the complication is in this case. Am I overlooking
> something subtle here? I assume Ra comes out dimensionless, in which
> case Ra could be declared as ordinary REAL.
>
It boils down to this subtlety:

0 deg. C = 273.15 K
(when dealing with thermodynamic temperatures)

0 deg. C = 0 K
(when dealing with temperature differences)

There is a delicate situation with \beta = (\partial ln \rho/\partial
T)_p, where T appears in a differential; however, for ideal gases, \beta
evaluates to 1/T_absolute.

> Just looking at what you wrote, I'd simply assign values with suitable
> units to each of the parameters/variables (declared as type PREAL) and
> then compute Ra as
>
> Ra = real(g*beta*DeltaT*(L**3)/(alpha*nu))
>
> If the argument to real() computes to something other than a
> non-dimensional number, it will generate a run-time error.
>
> - Grant
>
> ________________________________________

> Prof. Grant W. Petty
> Atmospheric and Oceanic Sciences
> 1225 W. Dayton Street
> University of Wisconsin-Madison
> Madison, WI 53706
>
> gpetty@aos.wisc.edu
> Tel: (608) 263-3265
> Fax: (608) 262-0166
>
>

Sponsored Links