From METCALF@crnvma.cern.ch Wed Aug  5 13:45:24 1992
Received: from vm.uni-c.dk by dkuug.dk via EUnet with SMTP (5.64+/8+bit/IDA-1.2.8)
	id AA13558; Wed, 5 Aug 92 13:45:24 +0200
Message-Id: <9208051145.AA13558@dkuug.dk>
Received: from vm.uni-c.dk by vm.uni-c.dk (IBM VM SMTP V2R2) with BSMTP id 3317;
   Wed, 05 Aug 92 13:45:57 DNT
Received: from CERNVM.cern.ch by vm.uni-c.dk (Mailer R2.07) with BSMTP id 4918;
 Wed, 05 Aug 92 13:45:55 DNT
Received: from CERNVM.CERN.CH (METCALF) by CERNVM.cern.ch (Mailer R2.08) with
 BSMTP id 3294; Wed, 05 Aug 92 13:45:58 SET
Date:         Wed, 05 Aug 92 13:40:28 SET
From: Michael Metcalf <METCALF@crnvma.cern.ch>
Subject:      f90 features
To: sc22wg5@dkuug.dk
X-Charset: ASCII
X-Char-Esc: 29

There has always been a lot of discussion about what features of f90
will be of real use to real users. The following module and test
program represent the first example at CERN of a user writing a
package in f90 as a test, from scratch. I think it speaks for itself.

                                    Mike Metcalf

! In this module all global variables are declared.

MODULE PYMODULE

IMPLICIT NONE

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Real numbers should have 64 bit precision throughout the program.
INTEGER, PARAMETER :: PYD = SELECTED_REAL_KIND(12,99)

! Standard constants for real numbers.
REAL(PYD), PARAMETER :: PY0 = 0_PYD, PY1 = 1_PYD, PY2 = 2_PYD, PY3 = 3_PYD, &
& PY4 = 4_PYD, PY5 = 5_PYD, PY6 = 6_PYD, PY7 = 7_PYD, PY8 = 8_PYD, &
& PY9 = 9_PYD, PY10 = 10_PYD, PY11 = 11_PYD, PY12 = 12_PYD, PY13 = 13_PYD, &
& PY14 = 14_PYD, PY15 = 15_PYD, PY16 = 16_PYD, PY17 = 17_PYD, PY18 = 18_PYD, &
& PY19 = 19_PYD, PY20 = 20_PYD, PY21 = 21_PYD, PY22 = 22_PYD, PY23 = 23_PYD
REAL(PYD), PARAMETER :: PYP1 = 0.1_PYD, PYP2 = 0.2_PYD, PYP3 = 0.3_PYD, &
& PYP4 = 0.4_PYD, PYP5 = 0.5_PYD, PYP6 = 0.6_PYD, PYP7 = 0.7_PYD, &
& PYP8 = 0.8_PYD, PYP9 = 0.9_PYD
REAL(PYD), PARAMETER :: PYP25 = 0.25_PYD, PYP75 = 0.75_PYD, PYP95 = 0.95_PYD
REAL(PYD), PARAMETER :: PYE2 = 1E2_PYD, PYE3 = 1E3_PYD, PYE4 = 1E4_PYD, &
& PYE5 = 1E5_PYD, PYE6 = 1E6_PYD, PYE7 = 1E7_PYD, PYE8 = 1E8_PYD, &
& PYE9 = 1E9_PYD, PYE10 = 1E10_PYD, PYE20 = 1E20_PYD
REAL(PYD), PARAMETER :: PYM2 = 1E-2_PYD, PYM3 = 1E-3_PYD, PYM4 = 1E-4_PYD, &
& PYM5 = 1E-5_PYD, PYM6 = 1E-6_PYD, PYM7 = 1E-7_PYD, PYM8 = 1E-8_PYD, &
& PYM9 = 1E-9_PYD, PYM10 = 1E-10_PYD, PYM20 = 1E-20_PYD

! Some basic constants.
REAL(PYD), PARAMETER :: PYPI = 3.141592653589793238_PYD     ! pi

! Machine precision related constants.
REAL(PYD), PARAMETER :: PYTINY = 1E-20_PYD
REAL(PYD), PARAMETER :: PYHUGE = 1E20_PYD

! Information for output.
INTEGER, PARAMETER :: PYUNIT = 6      ! output unit
INTEGER, PARAMETER :: PYNPAG = 55     ! number of lines per page
INTEGER, PRIVATE :: PYII              ! internal loop variable

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Declare PYHIST type to represent a one-dimensional histogram.
TYPE PYHIST
  CHARACTER(60) TITLE
  INTEGER       NCALL, NX
  REAL(PYD)     XMIN, XMAX, DX, WT(-2:100)
END TYPE PYHIST

! Define an empty histogram.
TYPE(PYHIST), PARAMETER :: PYHI0 = PYHIST(' ', 0, 0, PY0, PY0, PY0, &
& (/ (PY0, PYII=1,103) /) )

! Declare array of histograms, initially set empty.
INTEGER, PARAMETER :: PYNHIS = 100
TYPE(PYHIST), DIMENSION(PYNHIS) :: PYHI = &
& (/ (PYHI0, PYII=1,PYNHIS) /)

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Interfaces to all subroutines and functions in the package.
INTERFACE
    SUBROUTINE PYBOOK(ID, TITLE, NX, XMIN, XMAX)
    INTEGER,      INTENT(IN) :: ID, NX
    CHARACTER(*), INTENT(IN) :: TITLE
    REAL(PYD),    INTENT(IN) :: XMIN, XMAX
  END SUBROUTINE PYBOOK

  SUBROUTINE PYFILL(ID, X, WT)
    INTEGER,   INTENT(IN)           :: ID
    REAL(PYD), INTENT(IN)           :: X
    REAL(PYD), INTENT(IN), OPTIONAL :: WT
  END SUBROUTINE PYFILL

  SUBROUTINE PYSCALE(ID, FACTOR)
    INTEGER,   INTENT(IN) :: ID
    REAL(PYD), INTENT(IN) :: FACTOR
  END SUBROUTINE PYSCALE

  SUBROUTINE PYOPERA(ID1, OPER, ID2, ID3, FACTOR1, FACTOR2)
    INTEGER,      INTENT(IN)           :: ID1, ID2, ID3
    CHARACTER(*), INTENT(IN)           :: OPER
    REAL(PYD),    INTENT(IN), OPTIONAL :: FACTOR1, FACTOR2
  END SUBROUTINE PYOPERA

  SUBROUTINE PYCLEAR
  END SUBROUTINE PYCLEAR

  SUBROUTINE PYPRINT(ID)
    INTEGER, INTENT(IN) :: ID
  END SUBROUTINE PYPRINT
END INTERFACE

END MODULE PYMODULE


!*****************************************************************************
!*****************************************************************************


! Subroutine to book one histogram.

SUBROUTINE PYBOOK(ID, TITLE, NX, XMIN, XMAX)

USE PYMODULE, PYDUMMY => PYBOOK
IMPLICIT NONE

INTEGER,      INTENT(IN) :: ID, NX
CHARACTER(*), INTENT(IN) :: TITLE
REAL(PYD),    INTENT(IN) :: XMIN, XMAX

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Check that input sensible.
IF (ID < 1 .OR. ID > PYNHIS) RETURN
IF (NX < 1 .OR. NX > 100 .OR. XMAX-XMIN < NX*PYTINY) RETURN

! Reset histogram.
PYHI(ID) = PYHI0

! Book information on histogram.
PYHI(ID)%TITLE = TITLE
PYHI(ID)%NX    = NX
PYHI(ID)%XMIN  = XMIN
PYHI(ID)%XMAX  = XMAX
PYHI(ID)%DX    = (XMAX-XMIN) /NX

END SUBROUTINE PYBOOK


!*****************************************************************************
!*****************************************************************************


! Subroutine to fill entry into histogram.

SUBROUTINE PYFILL(ID, X, WT)

USE PYMODULE, PYDUMMY => PYFILL
IMPLICIT NONE

INTEGER,   INTENT(IN)           :: ID
REAL(PYD), INTENT(IN)           :: X
REAL(PYD), INTENT(IN), OPTIONAL :: WT

INTEGER   IX
REAL(PYD) WTF

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Check that input sensible. Define weight.
IF (ID < 1 .OR. ID > PYNHIS) RETURN
IF (PYHI(ID)%NX == 0) RETURN
PYHI(ID)%NCALL = PYHI(ID)%NCALL + 1
IF (PRESENT(WT)) THEN; WTF = WT; ELSE; WTF = PY1; ENDIF

! Find bin, store in histogram bin and/or in overflow bin.
IX = 1 + (X-PYHI(ID)%XMIN) /PYHI(ID)%DX
IF (IX < 1) THEN
  PYHI(ID)%WT(-1) = PYHI(ID)%WT(-1) + WTF
ELSEIF (IX > PYHI(ID)%NX) THEN
  PYHI(ID)%WT(-2) = PYHI(ID)%WT(-2) + WTF
ELSE
  PYHI(ID)%WT(0) = PYHI(ID)%WT(0) + WTF
  PYHI(ID)%WT(IX) = PYHI(ID)%WT(IX) + WTF
ENDIF

END SUBROUTINE PYFILL


!*****************************************************************************
!*****************************************************************************


! Subroutine to rescale contents of histogram.

SUBROUTINE PYSCALE(ID, FACTOR)

USE PYMODULE, PYDUMMY => PYSCALE
IMPLICIT NONE

INTEGER,   INTENT(IN) :: ID
REAL(PYD), INTENT(IN) :: FACTOR

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Check that input sensible.
IF (ID < 1 .OR. ID > PYNHIS) RETURN
IF (PYHI(ID)%NX == 0) RETURN

! Rescale contents.
PYHI(ID)%WT = FACTOR * PYHI(ID)%WT

END SUBROUTINE PYSCALE


!*****************************************************************************
!*****************************************************************************


! Subroutine to perform operations on histograms.

SUBROUTINE PYOPERA(ID1, OPER, ID2, ID3, FACTOR1, FACTOR2)

USE PYMODULE, PYDUMMY => PYOPERA
IMPLICIT NONE

INTEGER,      INTENT(IN)           :: ID1, ID2, ID3
CHARACTER(*), INTENT(IN)           :: OPER
REAL(PYD),    INTENT(IN), OPTIONAL :: FACTOR1, FACTOR2

INTEGER   IX
REAL(PYD) FACT1, FACT2, YMIN, WTX

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Check that input sensible.
IF (ID1 < 1 .OR. ID1 > PYNHIS) RETURN
IF (PYHI(ID1)%NX == 0) RETURN
SELECT CASE (OPER)
CASE('+', '-', '*', '/', 'M', 'm')
  IF (ID2 < 1 .OR. ID2 > PYNHIS) RETURN
  IF (PYHI(ID2)%NX == 0) RETURN
END SELECT
IF (ID3 < 1 .OR. ID3 > PYNHIS) RETURN
IF (PYHI(ID3)%NX == 0) RETURN

! Multiplicative factors and number of entries.
IF (PRESENT(FACTOR1)) THEN; FACT1 = FACTOR1; ELSE; FACT1 = PY1; ENDIF
SELECT CASE (OPER)
CASE('+', '-', '*', '/')
  IF (PRESENT(FACTOR2)) THEN; FACT2 = FACTOR2; ELSE; FACT2 = PY1; ENDIF
  PYHI(ID3)%NCALL = PYHI(ID1)%NCALL + PYHI(ID2)%NCALL
CASE('A', 'a', 'S', 's', 'L', 'l')
  IF (PRESENT(FACTOR2)) THEN; FACT2 = FACTOR2; ELSE; FACT2 = PY0; ENDIF
  PYHI(ID3)%NCALL = PYHI(ID1)%NCALL
END SELECT

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Operations on histogram contents.
SELECT CASE (OPER)

CASE('+')
  PYHI(ID3)%WT = FACT1 * PYHI(ID1)%WT + FACT2 * PYHI(ID2)%WT

CASE('-')
  PYHI(ID3)%WT = FACT1 * PYHI(ID1)%WT - FACT2 * PYHI(ID2)%WT

CASE('*')
  PYHI(ID3)%WT = FACT1 * PYHI(ID1)%WT * FACT2 * PYHI(ID2)%WT

CASE('/')
  WHERE ( ABS(FACT2*PYHI(ID2)%WT) > PYTINY)
    PYHI(ID3)%WT = (FACT1 * PYHI(ID1)%WT) / (FACT2 * PYHI(ID2)%WT)
  ELSEWHERE
    PYHI(ID3)%WT = PY0
  ENDWHERE

CASE('A', 'a')
  PYHI(ID3)%WT = FACT1 * PYHI(ID1)%WT + FACT2

CASE('S', 's')
  PYHI(ID3)%WT = FACT1 * SQRT(MAX(PY0, PYHI(ID1)%WT)) + FACT2

CASE('L', 'l')
  YMIN=PYHUGE
  DO IX = 1, PYHI(ID1)%NX
    WTX = PYHI(ID1)%WT(IX)
    IF (WTX < YMIN .AND. WTX > PYTINY) YMIN = PYP8 * WTX
  ENDDO
  PYHI(ID3)%WT = FACT1 * LOG10(MAX(YMIN, PYHI(ID1)%WT)) + FACT2

CASE('M', 'm')
  WHERE (ABS(PYHI(ID1)%WT) > PYTINY)
    PYHI(ID2)%WT = PYHI(ID2)%WT / PYHI(ID1)%WT
    PYHI(ID3)%WT = SQRT(MAX(PY0, PYHI(ID3)%WT/PYHI(ID1)%WT - PYHI(ID2)%WT**2))
  ELSEWHERE
    PYHI(ID2)%WT = PY0
    PYHI(ID3)%WT = PY0
  ENDWHERE
  PYHI(ID1)%WT = FACT1 * PYHI(ID1)%WT
END SELECT

END SUBROUTINE PYOPERA


!*****************************************************************************
!*****************************************************************************


! Subroutine to print and reset histograms.

SUBROUTINE PYCLEAR

USE PYMODULE, PYDUMMY => PYCLEAR
IMPLICIT NONE

INTEGER ID

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Loop over histograms; print and reset contents afterwards.
DO ID = 1, PYNHIS
  IF (PYHI(ID)%NCALL > 0) THEN
    CALL PYPRINT(ID)
    PYHI(ID)%NCALL = 0
    PYHI(ID)%WT    = PY0
  ENDIF
ENDDO

END SUBROUTINE PYCLEAR


!*****************************************************************************
!*****************************************************************************


! Subroutine to print a histogram.

SUBROUTINE PYPRINT(ID)

USE PYMODULE, PYDUMMY => PYPRINT
IMPLICIT NONE

INTEGER, INTENT(IN) :: ID

INTEGER DATI(8), NX, IX, LINES, IPOW, IDY, IROW(100), IFRA(100), &
& IR, IRMIN, IRMAX, NEGWT
REAL(PYD) XMIN, XMAX, DX, XLOW, WT(-2:100), YMIN, YMAX, DY, DYPOW, &
& WTAN, SUMWT, SUMWTX, SUMWTXX, X, XMEAN, XRMS
CHARACTER(100) OUT
CHARACTER(3) FMTNX

! Arrays for step size in histogram and digit characters.
REAL(PYD), DIMENSION(10), PARAMETER :: DYSTEP = (/ .04_PYD, .05_PYD, &
& .06_PYD, .08_PYD, .10_PYD, .12_PYD, .15_PYD, .20_PYD, .25_PYD, .30_PYD /)
CHARACTER(1), DIMENSION(0:10), PARAMETER :: DIGIT = (/ '0', '1', '2', &
& '3', '4', '5', '6', '7', '8', '9', 'X' /)

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Check that input sensible and that histogram not empty.
IF (ID < 1 .OR. ID > PYNHIS) RETURN
IF (PYHI(ID)%NCALL == 0) THEN
  WRITE(PYUNIT, "(/5X,'Histogram no',I4,': no entries')" ) ID
  RETURN
ENDIF

! Print header.
CALL DATE_AND_TIME(VALUES=DATI)
WRITE(PYUNIT, "('1'/5X,'Histogram no',I4,4X,A60,I4,'-',I2,'-',I2,1X,I2, &
& ':',I2/)" ) ID, PYHI(ID)%TITLE, DATI(1), DATI(2), DATI(3), DATI(5), DATI(6)

! Frequently used histogram variables.
NX   = PYHI(ID)%NX
XMIN = PYHI(ID)%XMIN
XMAX = PYHI(ID)%XMAX
DX   = PYHI(ID)%DX
WT   = PYHI(ID)%WT

! Variables for output and formats.
LINES = PYNPAG - 18
WRITE(FMTNX,"(I3.3)") NX

! Find minimum and maximum of histogram contents.
YMIN = WT(1)
YMAX = WT(1)
DO IX = 2, NX
  IF (WT(IX) < YMIN) YMIN = WT(IX)
  IF (WT(IX) > YMAX) YMAX = WT(IX)
ENDDO

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Only print histogram if spread big enough.
IF (YMAX-YMIN > LINES * DYSTEP(1) * PYM9) THEN

! Set power and spacing of histogram scale.
  IF (YMIN > PY0 .AND. YMIN < PYP1*YMAX) YMIN = PY0
  IF (YMAX < PY0 .AND. YMAX > PYP1*YMIN) YMAX = PY0
  IPOW = FLOOR(LOG10(YMAX-YMIN))
  IF (YMAX-YMIN < LINES * DYSTEP(1)  * PY10**IPOW) IPOW = IPOW - 1
  IF (YMAX-YMIN > LINES * DYSTEP(10) * PY10**IPOW) IPOW = IPOW + 1
  DYPOW = DYSTEP(1)
  DO IDY=1,9
    IF (YMAX-YMIN < LINES * DYSTEP(IDY) * PY10**IPOW) EXIT
    DYPOW = DYSTEP(IDY+1)
  ENDDO
  DY = DYPOW * PY10**IPOW

! Find extent of each histogram bin.
  DO IX = 1, NX
    WTAN = ABS(WT(IX)) /DY
    IROW(IX) = SIGN(WTAN + PYP95, WT(IX))
    IFRA(IX) = 1 + PY10*(WTAN + PYP95 - AINT(WTAN+PYP95))
  ENDDO
  IRMIN = SIGN(ABS(YMIN)/DY + PYP95, YMIN)
  IRMAX = SIGN(ABS(YMAX)/DY + PYP95, YMAX)

! Write histogram proper, line by line.
  DO IR = IRMAX, IRMIN, -1
    IF (IR == 0) CYCLE
    OUT = ' '
    DO IX = 1, NX
      IF (IR == IROW(IX)) OUT(IX:IX) = DIGIT(IFRA(IX))
      IF (IR*(IROW(IX)-IR) > 0) OUT(IX:IX) = DIGIT(10)
    ENDDO
    WRITE(PYUNIT, "(3X,F7.2,'*10**',I2,4X,A"//FMTNX//")" ) &
    & IR*DYPOW, IPOW, OUT(1:NX)
  ENDDO

! Write histogram bin contents.
  IPOW = FLOOR( LOG10(MAX(YMAX,-YMIN)) + PYM4)
  OUT=' '
  DO IX = 1, NX
    IF (WT(IX) < -PY10**(IPOW-4)) OUT(IX:IX) = '-'
    IROW(IX) = NINT( PY10**(3-IPOW) * ABS(WT(IX)) )
  ENDDO
  WRITE(PYUNIT, "(/9X,'Contents',4X,A"//FMTNX//")" ) OUT(1:NX)
  DO IR = 4, 1, -1
    DO IX = 1, NX
      OUT(IX:IX) = DIGIT( MOD(IROW(IX),10**IR) /10**(IR-1) )
    ENDDO
    WRITE(PYUNIT, "(10X,'*10**',I2,4X,A"//FMTNX//")" ) &
    & IPOW+IR-4, OUT(1:NX)
  ENDDO

! Write low edge of histogram bins.
  IPOW = FLOOR( LOG10(MAX(-XMIN,XMAX-DX)) + PYM4)
  OUT=' '
  DO IX = 1, NX
    XLOW = XMIN + (IX-1) * DX
    IF (XLOW < -PY10**(IPOW-3)) OUT(IX:IX) = '-'
    IROW(IX) = NINT( PY10**(2-IPOW) * ABS(XLOW) )
  ENDDO
  WRITE(PYUNIT, "(/9X,'Low edge',4X,A"//FMTNX//")" ) OUT(1:NX)
  DO IR = 3, 1, -1
    DO IX = 1, NX
      OUT(IX:IX) = DIGIT( MOD(IROW(IX),10**IR) /10**(IR-1) )
    ENDDO
    WRITE(PYUNIT, "(10X,'*10**',I2,4X,A"//FMTNX//")" ) &
    & IPOW+IR-3, OUT(1:NX)
  ENDDO
ENDIF

! -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

! Calculate histogram statistics.
NEGWT   = 0
SUMWT   = PY0
SUMWTX  = PY0
SUMWTXX = PY0
DO IX = 1, NX
  IF (WT(IX) < -PYTINY) NEGWT = NEGWT + 1
  X = XMIN + (IX-PYP5) * DX
  SUMWT   = SUMWT   + WT(IX)
  SUMWTX  = SUMWTX  + WT(IX) * X
  SUMWTXX = SUMWTXX + WT(IX) * X**2
ENDDO

! Write histogram statistics; mean and rms only if all weights positive.
IF (NEGWT == 0) THEN
  XMEAN = SUMWTX/MAX(PYTINY,SUMWT)
  XRMS  = SQRT(MAX(PY0, SUMWTXX/MAX(PYTINY,SUMWT) - XMEAN**2 ))
  WRITE(PYUNIT, "(/4X,'Entries  = ',I10, 5X,'Mean = ',ES10.3, &
  & 5X,'Underflow = ',ES10.3, 5X,'Low edge  = ',ES10.3 &
  & /4X,'All chan = ',ES10.3, 5X,'RMS  = ',ES10.3, &
  & 5X,'Overflow  = ',ES10.3, 5X,'High edge = ',ES10.3)" ) &
  & PYHI(ID)%NCALL, XMEAN, WT(-1), XMIN, WT(0), XRMS, WT(-2), XMAX
ELSE
  WRITE(PYUNIT, "(/4X,'Entries  = ',I10, 5X,'Mean =  undefined', &
  & 5X,'Underflow = ',ES10.3, 5X,'Low edge  = ',ES10.3 &
  & /4X,'All chan = ',ES10.3, 5X,'RMS  =  undefined', &
  & 5X,'Overflow  = ',ES10.3, 5X,'High edge = ',ES10.3)" ) &
  & PYHI(ID)%NCALL, WT(-1), XMIN, WT(0), WT(-2), XMAX
ENDIF

END SUBROUTINE PYPRINT


!*****************************************************************************
!*****************************************************************************


PROGRAM MAIN

USE PYMODULE
IMPLICIT REAL(PYD) (A-H,O-Z)

CALL PYBOOK(1,'Exponential drop',50,PY0,PY5)
CALL PYBOOK(2,'Sine function',50,PY0,PY5)
CALL PYBOOK(3,'Straight line',50,PY0,PY5)
CALL PYBOOK(4,'Square root of sine',50,PY0,PY5)
CALL PYBOOK(10,'Ratio sine/exponent',50,PY0,PY5)

DO IX = 1, 49
  X  = 0.1*IX-0.05
  F1 = EXP(-X)
  F2 = SIN(X)
  CALL PYFILL(1,X,F1)
  CALL PYFILL(2,X,F2)
  CALL PYFILL(3,X)
  CALL PYFILL(4,X,F2)
ENDDO

CALL PYOPERA(2,'/',1,10)
CALL PYOPERA(4,'S',4,4)
CALL PYCLEAR

END