PROGRAM realtype
	IMPLICIT NONE

	! This program discusses the representation of the type REAL
	! on computers and how code can be made portable.
	! Real numbers on a computer are represented in the form:
	! 		+-.d1 d2 d3 ....dP * B**r
	! where 
	!	B is the base (10 for the user)
	!	P is the precision (number of decimal places)
	! 	r lies in the range -R <= r <= +R where R is the 
	!	maximum exponent
	!	dn are decimal digits such that 
	!  	d1 is one of 1,2,..,9  and  dn is one of 0,1,...,9

	! The two parameters P and R need to be specified.

	! The computer uses the binary system so that B=2, d1=1 and 
	! dn= 0 or 1. Note that conversion between binary and decimal 
	! is not usually exact.
!--------------------------------------------------------------------------
! Declarations

	! First, consider real numbers which assume default values of P,R

	REAL :: real0

	! Second, consider reals where the precision and exponent range 
	! are user-prescribed. 
	! Use the intrinsic function SELECTED_REAL_KIND(P,R)
	! This ensures portability since any computer will guarantee 
	! the above minimum precision and range if it is supported.

	! Define the following KIND numbers, known as KIND TYPE PARAMETERS:
	INTEGER, PARAMETER :: real_6_30=SELECTED_REAL_KIND(P=6,R=30)
	INTEGER, PARAMETER :: real_10_100=SELECTED_REAL_KIND(P=10,R=100)
	INTEGER, PARAMETER :: real_15_300=SELECTED_REAL_KIND(P=15,R=300)

	! If the precision P cannot be supported then -1 is returned; 
	! if the exponent range cannot be supported -2 is returned; 
	! if neither is supported then -3 is returned.
	! The computer will support at least two levels of precision and 
	! exponent range - historically referred to as single and double 
	! precision. On some machines, there may be a third, called 
	! 'extended precision'. These will have an associated KIND number 
	! (which may be different on different machines). 
	! This number can be discovered by using the the intrinsic function 
	! KIND(x) which returns the integer kind number of the variable x. 

	! Now declare example variables of each type.
	REAL(KIND=4) :: real1, pi_4, x4, square4
	REAL(KIND=8) :: real2, real3, pi_8, pi_wrong, x8, square8, square4_to_8
	REAL(KIND=real_6_30) :: real4
	REAL(KIND=real_10_100) :: real5
	REAL(KIND=real_15_300) :: real6

!--------------------------------------------------------------------------
! real, defaults
	real0 = 0.12345e31

	! Also note how the largest (HUGE(real0)) and smallest (TINY(real0)) 
	! numbers of same type as real0 can be obtained.
	! The machine epsilon (EPSILON(real0)) is the smallest number which 
	! when added to one returns a number different from one.

	PRINT*
	PRINT*, 'Default real numbers:'
	PRINT*
	PRINT*, 'real0 = ',real0
	PRINT*, 'kind = ',KIND(real0)
	PRINT*, 'huge = ',HUGE(real0)
	PRINT*, 'tiny = ',TINY(real0)
	PRINT*, 'machine epsilon = ',EPSILON(real0)
	PRINT*

!--------------------------------------------------------------------------
! real, kind = 4; note that this is the default kind number.
	real1 = 0.12345678987654321e31

	PRINT*
	PRINT*, 'Kind=4 real numbers (same as default):'
	PRINT*
	PRINT*, 'real1 = ',real1
	PRINT*, 'kind = ',KIND(real1)
	PRINT*, 'huge = ',HUGE(real1)
	PRINT*, 'tiny = ',TINY(real1)
	PRINT*, 'machine epsilon = ',EPSILON(real1)
	PRINT*

!--------------------------------------------------------------------------
! real, kind = 8 (double precision)
	! Note how the kind number (here: _8)  must be attached to 
	! the end of constants. 

	real2 = 0.12345678987654321e-20_8

	! Observe what happens if this is omitted (treated as kind=4 default,
	! is truncated accordingly, and when assigned to real3 (of kind=8) 
	! has lower order decimal places inserted which correspond to a 
	! different number.

	real3 = 0.12345678987654321e-20

	PRINT*
	PRINT*, 'Kind=8 real numbers (double precision):'
	PRINT*
	PRINT*, 'real2 = ',real2
	PRINT*, 'real3 = ',real3
	PRINT*, 'kind = ',KIND(real2)
	PRINT*, 'huge = ',HUGE(real2)
	PRINT*, 'tiny = ',TINY(real2)
	PRINT*, 'machine epsilon = ',EPSILON(real2)
	PRINT*

!--------------------------------------------------------------------------
! real kind=real_6_30
	! Suppose we wish to ensure a minimum precision and exponent 
	! range for our code which is to be used on a variety of platforms.
	! The above kind number will ensure AT LEAST 6 decimal places of
	! precision and an exponent range of AT LEAST +-30.
	! This will translate to kind=4 on Tower.

	real4 = 1.234e12_real_6_30

	PRINT*
	PRINT*, 'Kind=real_6_30:'
	PRINT*
	PRINT*, 'real4 = ',real4
	PRINT*, 'kind = ',KIND(real4)
	PRINT*, 'huge = ',HUGE(real4)
	PRINT*, 'tiny = ',TINY(real4)
	PRINT*, 'machine epsilon = ',EPSILON(real4)
	PRINT*

!--------------------------------------------------------------------------
! real kind=real_10_100
	! The above kind number will ensure AT LEAST 10 decimal places of
	! precision and an exponent range of AT LEAST +-100.
	! This will translate to kind=8 on Tower.

	real5 = 0.1234567898e-95_real_10_100

	PRINT*
	PRINT*, 'Kind=real_10_100:'
	PRINT*
	PRINT*, 'real5 = ',real5
	PRINT*, 'kind = ',KIND(real5)
	PRINT*, 'huge = ',HUGE(real5)
	PRINT*, 'tiny = ',TINY(real5)
	PRINT*, 'machine epsilon = ',EPSILON(real5)
	PRINT*

!--------------------------------------------------------------------------
! real kind=real_15_300
	! Ofcourse it is sensible to exploit the full precision and 
	! exponent range offered by a machine so that a choice of P,R 
	! approaching those corresponding to kind=8 is appropriate.
	! The above kind number will ensure AT LEAST 15 decimal places of
	! precision and an exponent range of AT LEAST +-300.
	! This will translate to kind=8 on Tower.

	real6 = 0.123456789876543e-95_real_10_100

	PRINT*
	PRINT*, 'Kind=real_15_300:'
	PRINT*
	PRINT*, 'real6 = ',real6
	PRINT*, 'kind = ',KIND(real6)
	PRINT*, 'huge = ',HUGE(real6)
	PRINT*, 'tiny = ',TINY(real6)
	PRINT*, 'machine epsilon = ',EPSILON(real6)
	PRINT*

!--------------------------------------------------------------------------
! NOTE: I stress again that when working to higher precision make sure 
! that you represent constants to that higher precision by appending the 
! the kind number. 

	!Observe the following: why is pi_wrong wrong?
	pi_4 = 4.0_4*ATAN(1.0)
	pi_8 = 4.0_8*ATAN(1.0_8)
	pi_wrong = 4.0*ATAN(1.0)

	PRINT*
	PRINT*, pi_4, pi_8, pi_wrong
	PRINT*

	! Another example:
	x4 = SQRT(0.2_4)
	x8 = SQRT(0.2_8)
	square4 = x4*x4
	square8 = x8*x8
	square4_to_8 = x4*x4

	PRINT*
	PRINT*, ' The expected answer is 0.2 but the actual answers are '
	PRINT*
	PRINT*, square4, square8, square4_to_8
	PRINT*

END PROGRAM realtype
!--------------------------------------------------------------------------
