fortranparser.f90

Go to the documentation of this file.

!
! Copyright (c) 2000-2008, Roland Schmehl. All rights reserved.
!
! This software is distributable under the BSD license. See the terms of the
! BSD license in the documentation provided with this software.
!
MODULE fortranparser
  !------- -------- --------- --------- --------- --------- --------- --------- -------
  ! Fortran 2008 function parser
  !------- -------- --------- --------- --------- --------- --------- --------- -------
  ! 
  ! This is an OOP Fortran 2008 version of the original fparser by Roland Schmehl. This simple class
  ! wrapping of the original fparser has been developed by Jacopo Chevallard, and it is available on
  ! the GitHub repository https://github.com/jacopo-chevallard/FortranParser.
  ! 
  ! For comments and bug reports, please open an issue on
  ! https://github.com/jacopo-chevallard/FortranParser/issues
  !
  ! This function parser module is intended for applications where a set of mathematical
  ! fortran-style expressions is specified at runtime and is then evaluated for a large 
  ! number of variable values. This is done by compiling the set of function strings 
  ! into byte code, which is interpreted efficiently for the various variable values. 
  !
  ! The source code of the original fparser is available from http://fparser.sourceforge.net
  !
  ! Please send comments, corrections or questions realtive to the original fparser to its author:
  ! Roland Schmehl <roland.schmehl@alumni.uni-karlsruhe.de>
  !
  !------- -------- --------- --------- --------- --------- --------- --------- -------
  ! The function parser concept is based on a C++ class library written by  Juha 
  ! Nieminen <warp@iki.fi> available from http://warp.povusers.org/FunctionParser/
  !------- -------- --------- --------- --------- --------- --------- --------- -------
  USE fortranparser_parameters, ONLY: rn,is               ! Import KIND parameters

  IMPLICIT NONE

  PUBLIC                     :: equationparser

  !------- -------- --------- --------- --------- --------- --------- --------- -------
  PRIVATE

  INTEGER(is),                              PARAMETER :: cimmed   = 1,          &
                                                         cneg     = 2,          &
                                                         cadd     = 3,          & 
                                                         csub     = 4,          & 
                                                         cmul     = 5,          & 
                                                         cdiv     = 6,          & 
                                                         cpow     = 7,          & 
                                                         cabs     = 8,          &
                                                         cexp     = 9,          &
                                                         clog10   = 10,         &
                                                         clog     = 11,         &
                                                         csqrt    = 12,         &
                                                         csinh    = 13,         &
                                                         ccosh    = 14,         &
                                                         ctanh    = 15,         &
                                                         csin     = 16,         &
                                                         ccos     = 17,         &
                                                         ctan     = 18,         &
                                                         casin    = 19,         &
                                                         cacos    = 20,         &
                                                         catan    = 21,         &
                                                         varbegin = 22

  CHARACTER (LEN=1), DIMENSION(cAdd:cPow),  PARAMETER :: ops      = (/ '+',     &
                                                                       '-',     &
                                                                       '*',     &
                                                                       '/',     &
                                                                       '^' /)

  CHARACTER (LEN=5), DIMENSION(cAbs:cAtan), PARAMETER :: funcs    = (/ 'abs  ', &
                                                                       'exp  ', &
                                                                       'log10', &
                                                                       'log  ', &
                                                                       'sqrt ', &
                                                                       'sinh ', &
                                                                       'cosh ', &
                                                                       'tanh ', &
                                                                       'sin  ', &
                                                                       'cos  ', &
                                                                       'tan  ', &
                                                                       'asin ', &
                                                                       'acos ', &
                                                                       'atan ' /)

  INTEGER, parameter  :: max_fun_length = 1024

  TYPE equationparser

    INTEGER(is), POINTER :: bytecode(:) => null()
    INTEGER              :: bytecodesize = 0
    REAL(rn),    POINTER :: immed(:) => null()
    INTEGER              :: immedsize = 0
    REAL(rn),    POINTER :: stack(:) => null()
    INTEGER              :: stacksize = 0
    INTEGER              :: stackptr = 0

    character(len=MAX_FUN_LENGTH) :: funcstring = ''
    character(len=MAX_FUN_LENGTH) :: funcstringorig = ''
    character(len=MAX_FUN_LENGTH), allocatable :: variablenames(:) 
    contains

      private

      procedure, public :: evaluate
      procedure:: parse
      procedure :: compile
      procedure :: addcompiledbyte
      procedure :: compilesubstr
      procedure :: mathitemindex
      procedure :: checksyntax

      final :: finalize

  END TYPE equationparser

  ! Class constructor
  interface equationparser
    procedure constructor
  end interface equationparser

CONTAINS

!*****************************************************************************************
  type (equationparser) function constructor(funcstr, var)

    CHARACTER (LEN=*),               INTENT(in) :: funcstr   ! Function string
    CHARACTER (LEN=*), DIMENSION(:), INTENT(in) :: var       ! Array with variable names

    constructor%ByteCode => null()
    constructor%Immed => null()
    constructor%Stack => null()

    constructor%ByteCodeSize = 0
    constructor%ImmedSize = 0
    constructor%StackSize = 0
    constructor%StackPtr = 0

    constructor%funcString = funcstr
    constructor%funcStringOrig = funcstr

    allocate(constructor%variableNames(size(var)))

    constructor%variableNames(:) = var(:)

    call constructor%parse()

  end function constructor

!*****************************************************************************************
  subroutine finalize(this)

    type(equationparser) :: this

    if (associated(this%ByteCode))  nullify(this%ByteCode)
    if (associated(this%Immed))     nullify(this%Immed)
    if (associated(this%Stack))     nullify(this%Stack)

  end subroutine finalize

!*****************************************************************************************
  SUBROUTINE parse(this)
    ! Parse ith function string FuncStr and compile it into bytecode
    class(equationparser) :: this

    CALL replace('**','^ ', this%funcString)                            ! Exponent into 1-Char. format

    CALL removespaces(this%funcString)                                 ! Condense function string

    CALL this%CheckSyntax()

    CALL this%Compile()                                ! Compile into bytecode

  END SUBROUTINE parse

!*****************************************************************************************
  FUNCTION evaluate(this, Val) RESULT (res)
    ! Evaluate bytecode of ith function for the values passed in array Val(:)
    class(equationparser) :: this
    REAL(rn), DIMENSION(:), INTENT(in) :: val                ! Variable values

    REAL(rn)                           :: res                ! Result
    INTEGER                            :: ip,              & ! Instruction pointer
                                          dp,              & ! Data pointer
                                          sp                 ! Stack pointer
    REAL(rn),                PARAMETER :: zero = 0._rn
    integer       :: evalerrtype

    dp = 1
    sp = 0
    evalerrtype=0

    DO ip=1,this%ByteCodeSize

       SELECT CASE (this%ByteCode(ip))

       CASE (cimmed); sp=sp+1; this%Stack(sp)=this%Immed(dp); dp=dp+1

       CASE   (cneg); this%Stack(sp)=-this%Stack(sp)

       CASE   (cadd); this%Stack(sp-1)=this%Stack(sp-1)+this%Stack(sp); sp=sp-1

       CASE   (csub); this%Stack(sp-1)=this%Stack(sp-1)-this%Stack(sp); sp=sp-1

       CASE   (cmul); this%Stack(sp-1)=this%Stack(sp-1)*this%Stack(sp); sp=sp-1

       CASE   (cdiv)

         IF (this%Stack(sp)==0._rn) THEN 
           evalerrtype=1
           res=zero
           exit
         ENDIF
         this%Stack(sp-1)=this%Stack(sp-1)/this%Stack(sp); sp=sp-1

       CASE   (cpow); this%Stack(sp-1)=this%Stack(sp-1)**this%Stack(sp); sp=sp-1

       CASE   (cabs); this%Stack(sp)=abs(this%Stack(sp))

       CASE   (cexp); this%Stack(sp)=exp(this%Stack(sp))

       CASE (clog10)

         IF (this%Stack(sp)<=0._rn) THEN
           evalerrtype=3
           res=zero
           exit
         ENDIF
         this%Stack(sp)=log10(this%Stack(sp))

       CASE   (clog)

         IF (this%Stack(sp)<=0._rn) THEN
           evalerrtype=3
           res=zero
           exit
         ENDIF
         this%Stack(sp)=log(this%Stack(sp))

       CASE  (csqrt)

         IF (this%Stack(sp)<0._rn) THEN
           evalerrtype=3
           res=zero
           exit
         ENDIF
         this%Stack(sp)=sqrt(this%Stack(sp))

       CASE  (csinh); this%Stack(sp)=sinh(this%Stack(sp))

       CASE  (ccosh); this%Stack(sp)=cosh(this%Stack(sp))

       CASE  (ctanh); this%Stack(sp)=tanh(this%Stack(sp))

       CASE   (csin); this%Stack(sp)=sin(this%Stack(sp))

       CASE   (ccos); this%Stack(sp)=cos(this%Stack(sp))

       CASE   (ctan); this%Stack(sp)=tan(this%Stack(sp))

       CASE  (casin) 

         IF ((this%Stack(sp)<-1._rn) .OR. (this%Stack(sp)>1._rn)) THEN
           evalerrtype=4
           res=zero
           exit
         ENDIF
         this%Stack(sp)=asin(this%Stack(sp))

       CASE  (cacos); 
         IF ((this%Stack(sp)<-1._rn).OR.(this%Stack(sp)>1._rn)) THEN
           evalerrtype=4
           res=zero
           exit
         ENDIF
         this%Stack(sp)=acos(this%Stack(sp))

       CASE  (catan); this%Stack(sp)=atan(this%Stack(sp))

       CASE  default; sp=sp+1; this%Stack(sp)=val(this%ByteCode(ip)-varbegin+1)

       END SELECT

    END DO

    IF (evalerrtype > 0) then
      WRITE(*,*)'*** Error: ',evalerrmsg(evalerrtype)
    else
      res = this%Stack(1)
    endif

  END FUNCTION evaluate

!*****************************************************************************************
  SUBROUTINE checksyntax(this)
    ! Check syntax of function string,  returns 0 if syntax is ok
    class(equationparser) :: this
    INTEGER(is)                                 :: n
    CHARACTER (LEN=1)                           :: c
    REAL(rn)                                    :: r
    LOGICAL                                     :: err
    INTEGER                                     :: parcnt, & ! Parenthesis counter
                                                   j,ib,in,lfunc

    j = 1
    parcnt = 0
    lfunc = len_trim(this%funcString)
    step: DO
       IF (j > lfunc) CALL parseerrmsg(j, this%funcStringOrig)
       c = this%funcString(j:j)
       !-- -------- --------- --------- --------- --------- --------- --------- -------
       ! Check for valid operand (must appear)
       !-- -------- --------- --------- --------- --------- --------- --------- -------
       IF (c == '-' .OR. c == '+') THEN                      ! Check for leading - or +
          j = j+1
          IF (j > lfunc) CALL parseerrmsg(j, this%funcStringOrig, 'Missing operand')
          c = this%funcString(j:j)
          IF (any(c == ops)) CALL parseerrmsg(j, this%funcStringOrig, 'Multiple operators')
       END IF
       n = mathfunctionindex(this%funcString(j:))
       IF (n > 0) THEN                                       ! Check for math function
          j = j+len_trim(funcs(n))
          IF (j > lfunc) CALL parseerrmsg(j, this%funcStringOrig, 'Missing function argument')
          c = this%funcString(j:j)
          IF (c /= '(') CALL parseerrmsg(j, this%funcStringOrig, 'Missing opening parenthesis')
       END IF
       IF (c == '(') THEN                                    ! Check for opening parenthesis
          parcnt = parcnt+1
          j = j+1
          cycle step
       END IF
       IF (scan(c,'0123456789.') > 0) THEN                   ! Check for number
          r = realnum(this%funcString(j:),ib,in,err)
          IF (err) CALL parseerrmsg(j, this%funcStringOrig, 'Invalid number format:  '//this%funcString(j+ib-1:j+in-2))
          j = j+in-1
          IF (j > lfunc) EXIT
          c = this%funcString(j:j)
       ELSE                                                  ! Check for variable
          n = variableindex(this%funcString(j:),this%variableNames,ib,in)
          IF (n == 0) CALL parseerrmsg(j, this%funcStringOrig, 'Invalid element: '//this%funcString(j+ib-1:j+in-2))
          j = j+in-1
          IF (j > lfunc) EXIT
          c = this%funcString(j:j)
       END IF
       DO WHILE (c == ')')                                   ! Check for closing parenthesis
          parcnt = parcnt-1
          IF (parcnt < 0) CALL parseerrmsg(j, this%funcStringOrig, 'Mismatched parenthesis')
          IF (this%funcString(j-1:j-1) == '(') CALL parseerrmsg(j-1, this%funcStringOrig, 'Empty parentheses')
          j = j+1
          IF (j > lfunc) EXIT
          c = this%funcString(j:j)
       END DO
       !-- -------- --------- --------- --------- --------- --------- --------- -------
       ! Now, we have a legal operand: A legal operator or end of string must follow
       !-- -------- --------- --------- --------- --------- --------- --------- -------
       IF (j > lfunc) EXIT
       IF (any(c == ops)) THEN                               ! Check for multiple operators
          IF (j+1 > lfunc) CALL parseerrmsg(j, this%funcStringOrig)
          IF (any(this%funcString(j+1:j+1) == ops)) CALL parseerrmsg(j+1, this%funcStringOrig, 'Multiple operators')
       ELSE                                                  ! Check for next operand
          CALL parseerrmsg(j, this%funcStringOrig, 'Missing operator')
       END IF
       !-- -------- --------- --------- --------- --------- --------- --------- -------
       ! Now, we have an operand and an operator: the next loop will check for another 
       ! operand (must appear)
       !-- -------- --------- --------- --------- --------- --------- --------- -------
       j = j+1
    END DO step
    IF (parcnt > 0) CALL parseerrmsg(j, this%funcStringOrig, 'Missing )')
  END SUBROUTINE checksyntax

!*****************************************************************************************
  FUNCTION evalerrmsg(EvalErrType) RESULT (msg)
    ! Return error message
    integer, intent(in) :: evalerrtype
    CHARACTER (LEN=*), DIMENSION(4), PARAMETER :: m = (/ 'Division by zero                ', &
                                                         'Argument of SQRT negative       ', &
                                                         'Argument of LOG negative        ', &
                                                         'Argument of ASIN or ACOS illegal' /)
    CHARACTER (LEN=LEN(m))                     :: msg
    !----- -------- --------- --------- --------- --------- --------- --------- -------
    IF (evalerrtype < 1 .OR. evalerrtype > SIZE(m)) THEN
       msg = ''
    ELSE
       msg = m(evalerrtype)
    ENDIF

  END FUNCTION evalerrmsg

!*****************************************************************************************
  SUBROUTINE parseerrmsg (j, FuncStr, Msg)
    ! Print error message and terminate program
    INTEGER,                     INTENT(in) :: j
    CHARACTER (LEN=*),           INTENT(in) :: funcstr       ! Original function string
    CHARACTER (LEN=*), OPTIONAL, INTENT(in) :: msg

    INTEGER                                 :: k

    IF (present(msg)) THEN
       WRITE(*,*) '*** Error in syntax of function string: '//msg
    ELSE
       WRITE(*,*) '*** Error in syntax of function string:'
    ENDIF

    WRITE(*,*)
    WRITE(*,'(A)') ' '//funcstr

    WRITE(*,'(A)') '?'
    stop
  END SUBROUTINE parseerrmsg

!*****************************************************************************************
  FUNCTION operatorindex (c) RESULT (n)
    ! Return operator index
    CHARACTER (LEN=1), INTENT(in) :: c
    INTEGER(is)                   :: n,j

    n = 0

    DO j=cadd,cpow
       IF (c == ops(j)) THEN
          n = j
          EXIT
       END IF
    END DO

  END FUNCTION operatorindex

!*****************************************************************************************
  FUNCTION mathfunctionindex (str) RESULT (n)
    ! Return index of math function beginnig at 1st position of string str
    CHARACTER (LEN=*), INTENT(in) :: str

    INTEGER(is)                   :: n,j
    INTEGER                       :: k
    CHARACTER (LEN=LEN(Funcs))    :: fun

    n = 0

    DO j=cabs,catan                                          ! Check all math functions
       k = min(len_trim(funcs(j)), len(str))   
       CALL lowcase(str(1:k), fun)
       IF (fun == funcs(j)) THEN                             ! Compare lower case letters
          n = j                                              ! Found a matching function
          EXIT
       END IF
    END DO

  END FUNCTION mathfunctionindex

!*****************************************************************************************
  FUNCTION variableindex (str, Var, ibegin, inext) RESULT (n)
    !----- -------- --------- --------- --------- --------- --------- --------- -------
    ! Return index of variable at begin of string str (returns 0 if no variable found)
    !----- -------- --------- --------- --------- --------- --------- --------- -------
    IMPLICIT NONE
    CHARACTER (LEN=*),               INTENT(in) :: str       ! String
    CHARACTER (LEN=*), DIMENSION(:), INTENT(in) :: var       ! Array with variable names
    INTEGER(is)                                 :: n         ! Index of variable
    INTEGER, OPTIONAL,              INTENT(out) :: ibegin, & ! Start position of variable name
                                                   inext     ! Position of character after name
    INTEGER                                     :: j,ib,in,lstr
    !----- -------- --------- --------- --------- --------- --------- --------- -------
    n = 0
    lstr = len_trim(str)
    IF (lstr > 0) THEN
       DO ib=1,lstr                                          ! Search for first character in str
          IF (str(ib:ib) /= ' ') EXIT                        ! When lstr>0 at least 1 char in str
       END DO                        
       DO in=ib,lstr                                         ! Search for name terminators
          IF (scan(str(in:in),'+-*/^) ') > 0) EXIT
       END DO
       DO j=1,SIZE(var)
          IF (str(ib:in-1) == var(j)) THEN                     
             n = j                                           ! Variable name found
             EXIT
          END IF
       END DO
    END IF
    IF (present(ibegin)) ibegin = ib
    IF (present(inext))  inext  = in
  END FUNCTION variableindex

!*****************************************************************************************
  SUBROUTINE removespaces (str)
    ! Remove Spaces from string, remember positions of characters in old string
    CHARACTER (LEN=*), INTENT(inout) :: str

    INTEGER                          :: k,lstr

    lstr = len_trim(str)

    k = 1

    DO WHILE (str(k:lstr) /= ' ')                             
       IF (str(k:k) == ' ') THEN
          str(k:lstr)  = str(k+1:lstr)//' '                  ! Move 1 character to left
          k = k-1
       END IF
       k = k+1
    END DO

  END SUBROUTINE removespaces

!*****************************************************************************************
  SUBROUTINE replace (ca,cb,str)
    ! Replace ALL appearances of character set ca in string str by character set cb
    CHARACTER (LEN=*),       INTENT(in) :: ca
    CHARACTER (LEN=LEN(ca)), INTENT(in) :: cb                ! LEN(ca) must be LEN(cb)
    CHARACTER (LEN=*),    INTENT(inout) :: str

    INTEGER                             :: j,lca

    lca = len(ca)

    DO j=1,len_trim(str)-lca+1
       IF (str(j:j+lca-1) == ca) str(j:j+lca-1) = cb
    END DO

  END SUBROUTINE replace

!*****************************************************************************************
  SUBROUTINE compile(this)
    ! Compile i-th function string F into bytecode
    class(equationparser) :: this
    INTEGER                                     :: istat

    IF (ASSOCIATED(this%ByteCode)) DEALLOCATE ( this%ByteCode, &
                                                   this%Immed,    &
                                                   this%Stack     )
    this%ByteCodeSize = 0
    this%ImmedSize    = 0
    this%StackSize    = 0
    this%StackPtr     = 0

    CALL this%CompileSubstr(1,len_trim(this%funcString))               ! Compile string to determine size

    ALLOCATE ( this%ByteCode(this%ByteCodeSize), & 
               this%Immed(this%ImmedSize),       &
               this%Stack(this%StackSize),       &
               stat = istat                            )
    IF (istat /= 0) THEN
       WRITE(*,*) '*** Parser error: Memmory allocation for byte code failed'
       stop
    ELSE
       this%ByteCodeSize = 0
       this%ImmedSize    = 0
       this%StackSize    = 0
       this%StackPtr     = 0
       CALL this%CompileSubstr(1,len_trim(this%funcString))            ! Compile string into bytecode
    END IF

  END SUBROUTINE compile

!*****************************************************************************************
  SUBROUTINE addcompiledbyte(this, b)
    ! Add compiled byte to bytecode
    class(equationparser) :: this
    INTEGER(is), INTENT(in) :: b                             ! Value of byte to be added

    this%ByteCodeSize = this%ByteCodeSize + 1

    IF (ASSOCIATED(this%ByteCode)) then
      this%ByteCode(this%ByteCodeSize) = b
    endif

  END SUBROUTINE addcompiledbyte

!*****************************************************************************************
  FUNCTION mathitemindex(this, b, e) RESULT (n)
    ! Return math item index, if item is real number, enter it into Comp-structure
    class(equationparser) :: this

    INTEGER,           INTENT(in) :: b,e                     ! First and last pos. of substring
    INTEGER(is)                                 :: n         ! Byte value of math item

    n = 0

    IF (scan(this%funcString(b:b),'0123456789.') > 0) THEN                 ! Check for begin of a number
       this%ImmedSize = this%ImmedSize + 1
       IF (ASSOCIATED(this%Immed)) this%Immed(this%ImmedSize) = realnum(this%funcString(b:e))
       n = cimmed
    ELSE                                                     ! Check for a variable
       n = variableindex(this%funcString(b:e), this%variableNames)
       IF (n > 0) n = varbegin+n-1
    END IF

  END FUNCTION mathitemindex

!*****************************************************************************************
  FUNCTION completelyenclosed (F, b, e) RESULT (res)
    ! Check if function substring F(b:e) is completely enclosed by a pair of parenthesis
    CHARACTER (LEN=*), INTENT(in) :: f                       ! Function substring
    INTEGER,           INTENT(in) :: b,e                     ! First and last pos. of substring

    LOGICAL                       :: res
    INTEGER                       :: j,k

    res=.false.

    IF (f(b:b) == '(' .AND. f(e:e) == ')') THEN
       k = 0
       DO j=b+1,e-1
          IF     (f(j:j) == '(') THEN
             k = k+1
          ELSEIF (f(j:j) == ')') THEN
             k = k-1
          END IF
          IF (k < 0) EXIT
       END DO
       IF (k == 0) res=.true.                                ! All opened parenthesis closed
    END IF

  END FUNCTION completelyenclosed

!*****************************************************************************************
  RECURSIVE SUBROUTINE compilesubstr(this, b, e)
    ! Compile i-th function string funcString into bytecode
    class(equationparser) :: this
    INTEGER,                         INTENT(in) :: b,e       ! Begin and end position substring

    INTEGER(is)                                 :: n
    INTEGER                                     :: b2,j,k,io
    CHARACTER (LEN=*),                PARAMETER :: calpha = 'abcdefghijklmnopqrstuvwxyz'// &
                                                            'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
    ! Check for special cases of substring

    IF     (this%funcString(b:b) == '+') THEN                              ! Case 1: funcString(b:e) = '+...'
!      WRITE(*,*)'1. funcString(b:e) = "+..."'
       CALL this%CompileSubstr(b+1, e)
       RETURN
    ELSEIF (completelyenclosed(this%funcString, b, e)) THEN               ! Case 2: funcString(b:e) = '(...)'
!      WRITE(*,*)'2. funcString(b:e) = "(...)"'
       CALL this%CompileSubstr(b+1, e-1)
       RETURN
    ELSEIF (scan(this%funcString(b:b), calpha) > 0) THEN        
       n = mathfunctionindex(this%funcString(b:e))
       IF (n > 0) THEN
          b2 = b+index(this%funcString(b:e),'(')-1
          IF (completelyenclosed(this%funcString, b2, e)) THEN             ! Case 3: funcString(b:e) = 'fcn(...)'
!            WRITE(*,*)'3. funcString(b:e) = "fcn(...)"'
             CALL this%CompileSubstr(b2+1, e-1)
             CALL this%AddCompiledByte(n)
             RETURN
          END IF
       END IF

    ELSEIF (this%funcString(b:b) == '-') THEN
       IF (completelyenclosed(this%funcString, b+1, e)) THEN              ! Case 4: this%funcString(b:e) = '-(...)'
!         WRITE(*,*)'4. this%funcString(b:e) = "-(...)"'
          CALL this%CompileSubstr(b+2, e-1)
          CALL this%AddCompiledByte(cneg)
          RETURN
       ELSEIF (scan(this%funcString(b+1:b+1),calpha) > 0) THEN
          n = mathfunctionindex(this%funcString(b+1:e))
          IF (n > 0) THEN
             b2 = b+index(this%funcString(b+1:e),'(')
             IF (completelyenclosed(this%funcString, b2, e)) THEN          ! Case 5: this%funcString(b:e) = '-fcn(...)'
!               WRITE(*,*)'5. this%funcString(b:e) = "-fcn(...)"'
                CALL this%CompileSubstr(b2+1, e-1);
                CALL this%AddCompiledByte(n)
                CALL this%AddCompiledByte(cneg)
                RETURN
             END IF
          END IF
       ENDIF
    END IF

    ! Check for operator in substring: check only base level (k=0), exclude expr. in ()

    DO io=cadd,cpow                                          ! Increasing priority +-*/^
       k = 0
       DO j=e,b,-1
          IF     (this%funcString(j:j) == ')') THEN
             k = k+1
          ELSEIF (this%funcString(j:j) == '(') THEN
             k = k-1
          END IF
          IF (k == 0 .AND. this%funcString(j:j) == ops(io) .AND. isbinaryop(j, this%funcString)) THEN
             IF (any(this%funcString(j:j) == ops(cmul:cpow)) .AND. this%funcString(b:b) == '-') THEN ! Case 6: this%funcString(b:e) = '-...Op...' with Op > -
!               WRITE(*,*)'6. this%funcString(b:e) = "-...Op..." with Op > -'
                CALL this%CompileSubstr(b+1, e)
                CALL this%AddCompiledByte(cneg)
                RETURN                 
             ELSE                                                        ! Case 7: this%funcString(b:e) = '...BinOp...'
!               WRITE(*,*)'7. Binary operator',this%funcString(j:j)
                CALL this%CompileSubstr(b, j-1)
                CALL this%CompileSubstr(j+1, e)
                CALL this%AddCompiledByte(operatorindex(ops(io)))
                this%StackPtr = this%StackPtr - 1
                RETURN
             END IF
          END IF
       END DO
    END DO

    ! Check for remaining items, i.e. variables or explicit numbers

    b2 = b

    IF (this%funcString(b:b) == '-') b2 = b2+1

    n = this%MathItemIndex(b2, e)

!   WRITE(*,*)'8. AddCompiledByte ',n
    CALL this%AddCompiledByte(n)

    this%StackPtr = this%StackPtr + 1
    IF (this%StackPtr > this%StackSize) this%StackSize = this%StackSize + 1

    IF (b2 > b) CALL this%AddCompiledByte(cneg)

  END SUBROUTINE compilesubstr

!*****************************************************************************************
  FUNCTION isbinaryop (j, F) RESULT (res)
    ! Check if operator F(j:j) in string F is binary operator
    ! Special cases already covered elsewhere:              (that is corrected in v1.1)
    ! - operator character F(j:j) is first character of string (j=1)
    INTEGER,           INTENT(in) :: j                       ! Position of Operator
    CHARACTER (LEN=*), INTENT(in) :: f                       ! String

    LOGICAL                       :: res                     ! Result
    INTEGER                       :: k
    LOGICAL                       :: dflag,pflag

    res=.true.

    IF (f(j:j) == '+' .OR. f(j:j) == '-') THEN               ! Plus or minus sign:
       IF (j == 1) THEN                                      ! - leading unary operator ?
          res = .false.
       ELSEIF (scan(f(j-1:j-1),'+-*/^(') > 0) THEN           ! - other unary operator ?
          res = .false.
       ELSEIF (scan(f(j+1:j+1),'0123456789') > 0 .AND. &     ! - in exponent of real number ?
               scan(f(j-1:j-1),'eEdD')       > 0) THEN
          dflag=.false.; pflag=.false.
          k = j-1
          DO WHILE (k > 1)                                   !   step to the left in mantissa 
             k = k-1
             IF     (scan(f(k:k),'0123456789') > 0) THEN
                dflag=.true.
             ELSEIF (f(k:k) == '.') THEN
                IF (pflag) THEN
                   EXIT                                      !   * EXIT: 2nd appearance of '.'
                ELSE
                   pflag=.true.                              !   * mark 1st appearance of '.'
                ENDIF
             ELSE
                EXIT                                         !   * all other characters
             END IF
          END DO
          IF (dflag .AND. (k == 1 .OR. scan(f(k:k),'+-*/^(') > 0)) res = .false.
       END IF
    END IF
  END FUNCTION isbinaryop

!*****************************************************************************************
  FUNCTION realnum(str, ibegin, inext, error) RESULT (res)
    ! Get real number from string - Format: [blanks][+|-][nnn][.nnn][e|E|d|D[+|-]nnn]
    CHARACTER (LEN=*),  INTENT(in) :: str                    ! String
    REAL(rn)                       :: res                    ! Real number
    INTEGER, OPTIONAL, INTENT(out) :: ibegin,              & ! Start position of real number
                                      inext                  ! 1st character after real number
    LOGICAL, OPTIONAL, INTENT(out) :: error                  ! Error flag

    INTEGER                        :: ib,in,istat
    LOGICAL                        :: bflag,               & ! .T. at begin of number in str
                                      inman,               & ! .T. in mantissa of number
                                      pflag,               & ! .T. after 1st '.' encountered
                                      eflag,               & ! .T. at exponent identifier 'eEdD'
                                      inexp,               & ! .T. in exponent of number
                                      dinman,              & ! .T. if at least 1 digit in mant.
                                      dinexp,              & ! .T. if at least 1 digit in exp.
                                      err                    ! Local error flag
    !----- -------- --------- --------- --------- --------- --------- --------- -------
    bflag=.true.; inman=.false.; pflag=.false.; eflag=.false.; inexp=.false.
    dinman=.false.; dinexp=.false.
    ib   = 1
    in   = 1
    DO WHILE (in <= len_trim(str))
       SELECT CASE (str(in:in))
       CASE (' ')                                            ! Only leading blanks permitted
          ib = ib+1
          IF (inman .OR. eflag .OR. inexp) EXIT
       CASE ('+','-')                                        ! Permitted only
          IF     (bflag) THEN           
             inman=.true.; bflag=.false.                     ! - at beginning of mantissa
          ELSEIF (eflag) THEN               
             inexp=.true.; eflag=.false.                     ! - at beginning of exponent
          ELSE
             EXIT                                            ! - otherwise STOP
          ENDIF
       CASE ('0':'9')                                        ! Mark
          IF     (bflag) THEN           
             inman=.true.; bflag=.false.                     ! - beginning of mantissa
          ELSEIF (eflag) THEN               
             inexp=.true.; eflag=.false.                     ! - beginning of exponent
          ENDIF
          IF (inman) dinman=.true.                           ! Mantissa contains digit
          IF (inexp) dinexp=.true.                           ! Exponent contains digit
       CASE ('.')
          IF     (bflag) THEN
             pflag=.true.                                    ! - mark 1st appearance of '.'
             inman=.true.; bflag=.false.                     !   mark beginning of mantissa
          ELSEIF (inman .AND..NOT.pflag) THEN
             pflag=.true.                                    ! - mark 1st appearance of '.'
          ELSE
             EXIT                                            ! - otherwise STOP
          END IF
       CASE ('e','E','d','D')                                ! Permitted only
          IF (inman) THEN
             eflag=.true.; inman=.false.                     ! - following mantissa
          ELSE
             EXIT                                            ! - otherwise STOP
          ENDIF
       CASE default
          EXIT                                               ! STOP at all other characters
       END SELECT
       in = in+1
    END DO
    err = (ib > in-1) .OR. (.NOT.dinman) .OR. ((eflag.OR.inexp).AND..NOT.dinexp)
    IF (err) THEN
       res = 0.0_rn
    ELSE
       READ(str(ib:in-1),*,iostat=istat) res
       err = istat /= 0
    END IF
    IF (present(ibegin)) ibegin = ib
    IF (present(inext))  inext  = in
    IF (present(error))  error  = err
  END FUNCTION realnum

!*****************************************************************************************
  SUBROUTINE lowcase (str1, str2)
    ! Transform upper case letters in str1 into lower case letters, result is str2
    IMPLICIT NONE
    CHARACTER (LEN=*),  INTENT(in) :: str1
    CHARACTER (LEN=*), INTENT(out) :: str2
    INTEGER                        :: j,k
    CHARACTER (LEN=*),   PARAMETER :: lc = 'abcdefghijklmnopqrstuvwxyz'
    CHARACTER (LEN=*),   PARAMETER :: uc = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'

    str2 = str1

    DO j=1,len_trim(str1)
       k = index(uc,str1(j:j))
       IF (k > 0) str2(j:j) = lc(k:k)
    END DO

  END SUBROUTINE lowcase

  END MODULE fortranparser