text:
!
!
!  This is a basic agreement governing the use of the fortran function
!  parser developed by Dr Stuart Midgley.
!
!  	AGREEMENT
!
!  	Use of the Function Parser (FP) is subject to the following
!  	conditions
!
!  	This agreement is to remain attached all the source code to the FP
!
!  	You DO NOT have the right to distribute the FP source code
!
!  	You DO NOT have the right to sell the FP source code
!
!  	You DO NOT have the right to sell any source code or compiled
!  	code or binary distribution based on the FP source code
!
!  	The source code is provided AS IS with NO warranty and
!  	NO guarantee that it will work
!
!  	You DO have the right to use and modify the source code
!  	for your own personal NOT-FOR-PROFIT use
!
!  	You CAN send me cash, pizza, post cards, encouragement and
!  	modified code as you see fit ;)
!
!  	END AGREEMENT
!
!  If you don't agree to the above agreement, then don't use it.  Quite
!  simple really.
!
!  If the code proves useful, then great.  Of course, if you want to use it
!  commercially, let me know and we can setup some sort of agreement.
!
!  Payment:  well, pay me what you think the code is worth.  If you want to
!  send cash, please email me for my address.  Otherwise, the
!  occasional email would be good.  I would be really greatful if you have
!  suggestions for improvements.
!
!  Let me know if you have any problems using the code or you find a bug...
!  I am more than happy to help out.
!
!  Dr Stuart Midgley 10 January 2007
!  sdm900@gmail.com
!
!

module USR_FunctionParser
  use PAR_FunctionParser !!((03-B-PAR_FunctionParser.f90))
  use LIB_FunctionParser !!((04-B-LIB_FunctionParser.f90))

  implicit none

  private
  public :: s_createfn, s_errprecision, s_evaluatefn, s_evaluateerr, s_destroyfn

  type cachefn
     real(fdp), dimension(:), pointer :: numbers
     integer, dimension(:), pointer :: commands
  end type cachefn

  character(len=variablelength), dimension(:), allocatable :: variables
  integer, dimension(:), allocatable :: lenvariables
  integer :: numvariables
  real(fdp), dimension(:), allocatable :: varerrors
  character(len=errorlength) :: error
  type(cachefn), dimension(:), pointer, save :: cachedfuncs

contains

  subroutine s_firstrun()
    logical, save :: firstrun = .true.
    integer :: i

    if (firstrun) then
       allocate(cachedfuncs(10))

       do i=1, 10
          nullify(cachedfuncs(i)%numbers, cachedfuncs(i)%commands)
       end do

       firstrun = .false.
    end if
  end subroutine s_firstrun

  subroutine s_seterror(t_error)
    character(len=*), intent(in) :: t_error

    if (error == 'OK' .and. len_trim(t_error) > 0) error = trim(t_error)
  end subroutine s_seterror

  function f_clean(func) result(tmpfunc)
    character(len=*), intent(in) :: func
    character(len=len(func)) :: tmpfunc
    integer :: i, c, l

    tmpfunc = adjustl(func)
    l = len_trim(tmpfunc)
    i = 1

    do while (i <= l)
       if (tmpfunc(i:i)== ' ') then
          tmpfunc(i:) = tmpfunc(i+1:)
          l = l - 1
       else
          c = ichar(tmpfunc(i:i))
          if ((c <= ichar('Z')) .and. (c >= ichar('A'))) tmpfunc(i:i) = achar(c-ichar('A')+ichar('a'))
          i = i + 1
       end if
    end do
  end function f_clean

  function f_numbervars(vars) result(numvars)
    character(len=*), intent(in) :: vars
    integer :: numvars
    character(len=len(vars)) :: tmpvars
    character(len=variablelength) :: tmpvar

    tmpvars = trim(adjustl(vars))
    numvars = 0

    do while (len_trim(tmpvars) > 0)
       read(tmpvars, *) tmpvar
       numvars = numvars + 1
       tmpvars = tmpvars(index(tmpvars, trim(tmpvar))+len_trim(tmpvar):)
    end do
  end function f_numbervars

  subroutine s_errprecision(number)
    real(fdp), intent(inout) :: number
    real(fdp) :: tmp1, tmp2

    !
    ! Format the error to the correct precision
    !

    if (number /= 0.0d0) then
       tmp1 = 10.0d0**(floor(log10(number)))
       tmp2 = anint(number/tmp1)

       if (tmp2 > 1.0d0) then
          number = tmp2*tmp1
       else
          tmp1 = 10.0d0**(floor(log10(number))-1.0d0)
          tmp2 = anint(number/tmp1)
          number = tmp2*tmp1
       end if
    end if
  end subroutine s_errprecision

  subroutine s_readvarerror(errorstring, errorval)
    character(len=*), intent(in) :: errorstring
    real(fdp) :: errorval
    integer :: command, lenerror
    character(len=errorlength) :: tmperror

    errorval = 0.0d0
    lenerror = len_trim(errorstring)
    command = nop_

    if (errorstring(1:1) == '[' .and. errorstring(lenerror:lenerror) == ']') then
       tmperror = errorstring(2:lenerror-1)

       call s_readnumber(tmperror, lenerror, command, errorval)

       if (command /= number_) call s_seterror('Not a number: '//tmperror)

       if (len_trim(tmperror) /= 0) call s_seterror('Not correct error format: '//tmperror)
    else
       call s_seterror('Not a correct error: '//errorstring)
    end if
  end subroutine s_readvarerror

  subroutine s_allocatevariables(vars)
    character(len=*), intent(in) :: vars
    integer :: numvars, i, bracketpos, j, c
    character(len=variablelength+errorlength), dimension(:), allocatable :: tmp

    numvars = f_numbervars(vars)

    allocate(tmp(numvars), variables(numvars), lenvariables(numvars), varerrors(numvars))

    read(vars, *) tmp

    do i = 1, numvars
       do j = 1, len_trim(tmp(i))
          c = ichar(tmp(i)(j:j))
          if ((c <= ichar('Z')) .and. (c >= ichar('A'))) tmp(i)(j:j) = achar(c-ichar('A')+ichar('a'))
       end do

       bracketpos = scan(tmp(i), "[")

       if (bracketpos /= 0) then
          variables(i) = trim(tmp(i)(:bracketpos-1))

          call s_readvarerror(trim(tmp(i)(bracketpos:)), varerrors(i))
       else
          variables(i) = trim(tmp(i))
          varerrors(i) = 0.0d0
       end if
    end do

    lenvariables = len_trim(variables)
    numvariables = numvars

    deallocate(tmp)
  end subroutine s_allocatevariables

  subroutine s_deallocatevariables()
    deallocate(variables, lenvariables, varerrors)
  end subroutine s_deallocatevariables

  subroutine s_funcpush(funcstack, funcstackn, command, number)
    type(cachefn), intent(inout) :: funcstack
    integer, intent(inout) :: funcstackn
    real(fdp), intent(in) :: number
    integer, intent(in) :: command
    type(cachefn) :: tmpstack

    funcstackn = funcstackn + 1

    if (funcstackn >= size(funcstack%numbers)) then
       allocate(tmpstack%numbers(funcstackn+5), tmpstack%commands(funcstackn+5))
       tmpstack%numbers(:funcstackn-1) = funcstack%numbers
       tmpstack%commands(:funcstackn-1) = funcstack%commands
       tmpstack%numbers(funcstackn:) = 0.0d0
       tmpstack%commands(funcstackn:) = end_
       deallocate(funcstack%numbers, funcstack%commands)
       funcstack%numbers => tmpstack%numbers
       funcstack%commands => tmpstack%commands
    end if

    funcstack%commands(funcstackn) = command
    funcstack%numbers(funcstackn) = number
  end subroutine s_funcpush

  subroutine s_checkoperations(func, lenfunc, command, nextcommand)
    integer, intent(in) :: lenfunc
    character(lenfunc), intent(in) :: func
    integer, intent(in) :: command, nextcommand

    select case(command)

    case(number_)
       select case(nextcommand)

       case(number_)
          call s_seterror('Number following number: '//func)

       case(variablestart:variableend)
          call s_seterror('Variable following number: '//func)

       case(bracket_)
          call s_seterror('Bracket following number: '//func)

       case(expressionstart:expressionend)
          call s_seterror('Expression following number: '//func)

       end select

    case(variablestart:variableend)
       select case(nextcommand)

       case(number_)
          call s_seterror('Number following variable: '//func)

       case(variablestart:variableend)
          call s_seterror('Variable following variable: '//func)

       case(bracket_)
          call s_seterror('Bracket following variable: '//func)

       case(expressionstart:expressionend)
          call s_seterror('Expression following variable: '//func)

       end select

    case(conditionalstart:conditionalend)
       select case(nextcommand)

       case(conditionalstart:conditionalend, opsstart:opsend)
          call s_seterror('Operator following operator: '//func)

       case(end_)
          call s_seterror('Premature end of function')

       end select

    case(opsstart:opsend, pmstart:pmend)
       select case(nextcommand)

       case(operatorstart:operatorend)
          call s_seterror('Operator following operator: '//func)

       case(end_)
          call s_seterror('Premature end of function')

       end select
    end select
  end subroutine s_checkoperations

  subroutine s_checkcomponents(components, command)
    integer, intent(in) :: components, command

    if (components /= commandcomps(command)) call s_seterror('Incorrect components for: '//expressions(command))
  end subroutine s_checkcomponents

  subroutine s_readnumber(func, lenfunc, command, number)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    integer, intent(out) :: command
    real(fdp), intent(out) :: number
    integer :: pos
    logical :: dflag, exponent, finish

    exponent = .false.
    dflag = .false.
    finish = .false.
    pos = 0
    number = 0.0d0
    command = nop_

    do while(pos < lenfunc .and. .not. finish)
       pos = pos + 1

       select case(func(pos:pos))

       case ('0':'9')
          if (dflag) then
             exponent = .true.
             dflag = .false.
          end if

       case ('d','e','D','E')
          if (dflag .or. exponent) then
             finish = .true.
          else
             dflag = .true.
          end if

       case('.')
          if (dflag .or. exponent) finish = .true.

       case('+', '-')
          if (dflag) then
             dflag = .false.
             exponent = .true.
          else
             finish = .true.
          end if

       case default
          finish = .true.

       end select
    end do

    if (pos == lenfunc .and. .not. finish) then
       read(func, *) number

       command = number_

       func = ''
       lenfunc = 0
    else
       read(func(:pos-1), *) number

       command = number_

       func = func(pos:)
       lenfunc = lenfunc - pos + 1
    end if

  end subroutine s_readnumber

  subroutine s_readtext(func, lenfunc, command, number)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    integer, intent(out) :: command
    real(fdp), intent(out) :: number
    integer :: i, funccut

    funccut = 0
    command = nop_
    number = 0.0d0

    do i = expressionstart, expressionend
       if (expressionslen(i) < lenfunc) then
          if ((func(:expressionslen(i)) == expressions(i)(:expressionslen(i))) &
               .and. expressionslen(i) > funccut ) then
             command = i
             funccut = expressionslen(i)-1
          end if
       end if
    end do

    if (command == nop_) then
       do i = 1, numvariables
          if (func(:min(lenvariables(i), lenfunc)) == variables(i) .and. &
               &lenvariables(i) > funccut) then
             command = variablestart+i
             number = varerrors(i)
             funccut = lenvariables(i)
          end if
       end do
    end if

    if (command == nop_) then
       call s_seterror('Unknown text: '//func)

       func = ''
       lenfunc = 0
    end if

    if (lenfunc > funccut) then
       func = func(funccut+1:)
       lenfunc = lenfunc - funccut
    else
       func = ''
       lenfunc = 0
    end if
  end subroutine s_readtext

  subroutine s_readoperator(func, lenfunc, command, order)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    integer, intent(out) :: command, order
    integer :: i, funccut

    funccut = 0
    command = nop_
    order = huge(order)

    do i = operatorstart, operatorend
       if (func(:min(operatorslen(i), lenfunc)) == operators(i) .and. &
            &operatorslen(i) > funccut ) then
          command = i
          order = operatorsorder(i)
          funccut = operatorslen(i)
       end if
    end do

    if (command == nop_) then
       call s_seterror('Unknown operator: '//func)

       func = ''
       lenfunc = 0
    end if

    if (lenfunc > funccut) then
       func = func(funccut+1:)
       lenfunc = lenfunc - funccut
    else
       func = ''
       lenfunc = 0
    end if
  end subroutine s_readoperator

  subroutine s_readbracket(func, lenfunc, subfunc, lensubfunc)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    character(len=len(func)), intent(out) :: subfunc
    integer, intent(out) :: lensubfunc
    integer :: pos, brackettype
    integer, dimension(bracketstart:bracketend) :: count

    count = 0
    brackettype = 0
    subfunc = ''

    select case(func(1:1))

    case('(')
       brackettype = roundbracket_

    case('[')
       brackettype = squarebracket_

    case('{')
       brackettype = squigbracket_

    end select

    do pos = 1, lenfunc
       select case(func(pos:pos))

       case('(')
          count(roundbracket_) = count(roundbracket_) + 1

       case(')')
          count(roundbracket_) = count(roundbracket_) - 1

       case('[')
          count(squarebracket_) = count(squarebracket_) + 1

       case(']')
          count(squarebracket_) = count(squarebracket_) - 1

       case('{')
          count(squigbracket_) = count(squigbracket_) + 1

       case('}')
          count(squigbracket_) = count(squigbracket_) - 1

       end select

       if (count(brackettype) == 0) then
          exit
       end if
    end do

    subfunc = func(2:pos-1)
    lensubfunc = pos - 2

    if (sum(abs(count)) /= 0) then
       call s_seterror('Unbalanced brackets: '//func)

       func = ''
       lenfunc = 0
    end if

    if (lensubfunc == 0) then
       call s_seterror('Empty brackets encountered: '//func)

       func = ''
       lenfunc = 0
    end if

    if (lenfunc < pos+1) then
       func = ''
       lenfunc = 0
    else
       func = func(pos+1:)
       lenfunc = lenfunc - pos
    end if
  end subroutine s_readbracket

  subroutine s_readcomma(func, lenfunc, command, order)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    integer, intent(out) :: command, order

    command = comma_
    order = 0

    if (lenfunc < 2) then
       func = ''
       lenfunc = 0
    else
       func = func(2:)
       lenfunc = lenfunc - 1
    end if
  end subroutine s_readcomma

  recursive subroutine s_parseatom(func, lenfunc, funcstack, funcstackn, command, number, order, components)
    type(cachefn), intent(inout) :: funcstack
    integer, intent(inout) :: funcstackn
    integer, intent(inout) :: lenfunc, components
    character(lenfunc), intent(inout) :: func
    integer, intent(in) :: command, order
    real(fdp), intent(in) :: number
    character(len=len(func)) :: tmpfunc
    integer :: lentmpfunc, nextcommand, nextorder, newcomponents
    real(fdp) :: nextnumber

    select case(command)

    case(number_, variablestart:variableend)
       call s_funcpush(funcstack, funcstackn, command, number)

    case(nop_, end_)

    case(comma_)
       call s_parseblock(command, func, lenfunc, funcstack, funcstackn, order, components)

       components = components + 1

    case(bracket_)
       call s_readbracket(func, lenfunc, tmpfunc, lentmpfunc)
       call s_parsefn(tmpfunc, lentmpfunc, funcstack, funcstackn, components)

    case(operatorstart:operatorend)
       call s_parseblock(command, func, lenfunc, funcstack, funcstackn, order, components)
       call s_funcpush(funcstack, funcstackn, command, number)

    case(expressionstart:expressionend)
       call s_getnextatom(func, lenfunc, nextcommand, nextnumber, nextorder)
       call s_checkoperations(func, lenfunc, command, nextcommand)

       newcomponents = 1
       call s_parseatom(func, lenfunc, funcstack, funcstackn, nextcommand, nextnumber, nextorder, newcomponents)
       call s_checkcomponents(newcomponents, command)
       call s_funcpush(funcstack, funcstackn, command, number)

    case default
       call s_seterror('Unknown input: '//func)

       func = ''
       lenfunc = 0

    end select
  end subroutine s_parseatom

  subroutine s_getnextatom(func, lenfunc, command, number, order)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    integer, intent(out) :: command, order
    real(fdp), intent(out) :: number

    command = nop_
    order = huge(order)
    number = 0.0d0

    if (lenfunc > 0) then
       select case(func(1:1))

       case('0':'9', '.')
          call s_readnumber(func, lenfunc, command, number)

       case('+', '-', '/', '*', '^', '>', '!', '<', '=')
          call s_readoperator(func, lenfunc, command, order)

       case('a':'z','A':'Z')
          call s_readtext(func, lenfunc, command, number)

       case('(','[','{')
          command = bracket_

       case(',')
          call s_readcomma(func, lenfunc, command, order)

       case default
          if (lenfunc > 1) call s_seterror('Unknown input: '//func)

          func = ''
          lenfunc = 0
          command = end_

       end select
    else
       command = end_
    end if
  end subroutine s_getnextatom

  subroutine s_peeknextatom(func, lenfunc, command, number, order)
    integer, intent(in) :: lenfunc
    character(lenfunc), intent(in) :: func
    integer, intent(out) :: command, order
    real(fdp), intent(out) :: number
    character(len=lenfunc) :: tmpfunc
    integer :: lentmpfunc

    command = nop_
    order = huge(order)
    tmpfunc = func
    lentmpfunc = lenfunc

    call s_getnextatom(tmpfunc, lentmpfunc, command, number, order)
  end subroutine s_peeknextatom

  recursive subroutine s_parseblock(command, func, lenfunc, funcstack, n, blockorder, components)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    type(cachefn), intent(inout) :: funcstack
    integer, intent(inout) :: n, components
    integer, intent(in) :: command, blockorder
    integer :: tcommand, nextcommand, nextorder
    real(fdp) :: nextnumber

    tcommand = command
    nextcommand = nop_
    nextorder = huge(nextorder)

    call s_peeknextatom(func, lenfunc, nextcommand, nextnumber, nextorder)
    call s_checkoperations(func, lenfunc, tcommand, nextcommand)

    if ((nextcommand == plus_ .or. nextcommand == minus_) .and. &
         count(command == blockseparator) > 0) call s_funcpush(funcstack, n, number_, 0.0d0)

    do while(nextorder > blockorder .and. nextcommand /= end_)
       call s_getnextatom(func, lenfunc, nextcommand, nextnumber, nextorder)
       call s_parseatom(func, lenfunc, funcstack, n, nextcommand, nextnumber, nextorder, components)
       tcommand = nextcommand
       call s_peeknextatom(func, lenfunc, nextcommand, nextnumber, nextorder)
    end do
  end subroutine s_parseblock

  recursive subroutine s_parsefn(func, lenfunc, funcstack, funcstackn, components)
    integer, intent(inout) :: lenfunc
    character(lenfunc), intent(inout) :: func
    type(cachefn), intent(inout) :: funcstack
    integer, intent(inout) :: funcstackn, components
    integer :: order, command
    real(fdp) :: number

    call s_getnextatom(func, lenfunc, command, number, order)

    if (command == plus_ .or. command == minus_) call s_funcpush(funcstack, funcstackn, number_, 0.0d0)

    do while(command /= end_)
       call s_parseatom(func, lenfunc, funcstack, funcstackn, command, number, order, components)
       call s_getnextatom(func, lenfunc, command, number, order)
    end do
  end subroutine s_parsefn

  subroutine s_evaluate(funcstack, funcstackn, varvals, result, errvals)
    type(cachefn), intent(in) :: funcstack
    integer, intent(inout) :: funcstackn
    real(fdp), dimension(:), intent(in) :: varvals
    real(fdp), dimension(:), intent(inout), optional :: errvals
    real(fdp), intent(out) :: result
    real(fdp), dimension(:), allocatable :: tmperrvals
    real(fdp), dimension(:), save, allocatable :: stack
    real(fdp), dimension(maxcommandcomps) :: tmp
    integer :: command, stackn
    real(fdp) :: number

    if (.not. allocated(stack)) then
       allocate(stack(0:size(funcstack%numbers)))
       stack = 0.0d0
    elseif (size(stack) < size(funcstack%numbers)+1) then
       deallocate(stack)
       allocate(stack(0:size(funcstack%numbers)))
       stack = 0.0d0
    end if

    command = nop_
    number = 0.0d0
    result = 0.0d0
    stackn = 0

    if (present(errvals)) then
       allocate(tmperrvals(size(errvals)))
       tmperrvals = errvals
    end if

    do while(command /= end_)
       funcstackn=funcstackn+1;command=funcstack%commands(funcstackn);number=funcstack%numbers(funcstackn)

       select case(command)

       case(number_)
          stackn=stackn+1;stack(stackn)=result
          result = number

       case(variablestart:variableend)
          stackn=stackn+1;stack(stackn)=result
          result = varvals(command-variablestart)
          if (present(errvals)) then
             tmperrvals(command-variablestart) = errvals(command-variablestart)*number
             result = result+tmperrvals(command-variablestart)
          end if

       case(pow_, caret_)
          tmp(1)=stack(stackn);stackn=stackn-1

          if (tmp(1) > 0.0d0) then
             result = tmp(1) ** result
          elseif (tmp(1) < 0.0d0) then
             if (aint(result) == result) then
                result = tmp(1) ** nint(result)
             else
                call s_seterror('x**y, x<0 and y not integer')
             end if
          elseif (result /= 0.0d0) then
             result = 0.0d0
          else
             call s_seterror('Undefined result 0**0')
          end if

       case(neq_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = merge(1.0d0, 0.0d0, tmp(1) /= result)

       case(eqeq_, eq_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = merge(1.0d0, 0.0d0, tmp(1) == result)

       case(gteq_, eqgt_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = merge(1.0d0, 0.0d0, tmp(1) >= result)

       case(lteq_, eqlt_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = merge(1.0d0, 0.0d0, tmp(1) <= result)

       case(gt_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = merge(1.0d0, 0.0d0, tmp(1) > result)

       case(lt_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = merge(1.0d0, 0.0d0, tmp(1) < result)

       case(mult_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = tmp(1) * result

       case(div_)
          tmp(1)=stack(stackn);stackn=stackn-1

          if (result /= 0.0d0) then
             result = tmp(1) / result
          else
             call s_seterror('x/y, y=0')
          end if

       case(plus_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = tmp(1) + result

       case(minus_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = tmp(1) - result

       case(cos_)
          result = cos(result)

       case(sin_)
          result = sin(result)

       case(tan_)
          result = tan(result)

       case(exp_)
          result = exp(result)

       case(log_)
          if (result > 0.0d0) then
             result = log(result)
          else
             call s_seterror('log(x), x<=0')
          end if

       case(log10_)
          if (result > 0.0d0) then
             result = log10(result)
          else
             call s_seterror('log10(x), x<=0')
          end if

       case(sqrt_)
          if (result >= 0.0d0) then
             result = sqrt(result)
          else
             call s_seterror('sqrt(x), x<0')
          end if

       case(cbrt_)
          result = f_cbrt(result)

       case(acos_)
          if (abs(result) <= 1.0d0) then
             result = acos(result)
          else
             call s_seterror('acos(x), |x|>1')
          end if

       case(asin_)
          if (abs(result) <= 1.0d0) then
             result = asin(result)
          else
             call s_seterror('asin(x), |x|>1')
          end if

       case(atan_)
          result = atan(result)

       case(cosh_)
          result = cosh(result)

       case(sinh_)
          result = sinh(result)

       case(tanh_)
          result = tanh(result)

       case(anint_)
          result = anint(result)

       case(aint_)
          result = aint(result)

       case(abs_)
          result = abs(result)

       case(delta_)
          result = merge(1.0d0, 0.0d0, result == 0.0d0)

       case(step_)
          result = merge(0.0d0, 1.0d0, result < 0.0d0)

       case(hat_)
          result = merge(1.0d0, 0.0d0, result <= 0.5d0 .and. result >= -0.5d0)

       case(max_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = max(result, tmp(1))

       case(min_)
          tmp(1)=stack(stackn);stackn=stackn-1
          result = min(result, tmp(1))

       case(besj0_)
          result = f_besj0(result)

       case(besj1_)
          result = f_besj1(result)

       case(besjn_)
          tmp(1)=stack(stackn);stackn=stackn-1
          if (aint(tmp(1)) == tmp(1)) then
             result = f_besjn(tmp(1), result)
          else
             call s_seterror('besjn(n, x), n not integer')
          end if

       case(besy0_)
          result = f_besy0(result)

       case(besy1_)
          result = f_besy1(result)

       case(besyn_)
          tmp(1)=stack(stackn);stackn=stackn-1
          if (aint(tmp(1)) == tmp(1) .and. result > 0.0d0) then
             result = f_besyn(tmp(1), result)
          else
             call s_seterror('besyn(n, x), n not integer or x/>0')
          end if

       case(logn_)
          tmp(1)=stack(stackn);stackn=stackn-1

          if (result > 0.0d0 .and. tmp(1) > 0.0d0) then
             result = log(result)/log(tmp(1))
          else
             call s_seterror('logn(x, y), x/>0 or y/>0')
          end if

       case(erf_)
          result = f_erf(result)

       case(erfc_)
          result = 1.0d0 - f_erf(result)

       case(lgamma_)
          if (result > 0.0d0) then
             result = f_lgamma(result)
          else
             call s_seterror('lgamma(x), x<=0')
          end if

       case(gamma_)
          result = f_gamma(result)

       case(csch_)
          result = sinh(result)

          if (result == 0.0d0) then
             call s_seterror('1/sinh(x), sinh(x) = 0')
          else
             result = 1.0d0/result
          end if

       case(sech_)
          result = 1.0d0/cosh(result)

       case(coth_)
          result = tanh(result)

          if (result == 0.0d0) then
             call s_seterror('1/tanh(x), tanh(x) = 0')
          else
             result = 1.0d0/result
          end if

       case(if_)
          tmp(1)=stack(stackn);stackn=stackn-1
          tmp(2)=stack(stackn);stackn=stackn-1

          result = merge(result, tmp(1), tmp(2) == 0)

       case(gauss_)
          result = exp(-result**2)

       case(sinc_)
          result = merge(1.0d0, sin(result)/result, result == 0.0d0)

       case(besi0_)
          result = f_besi0(result)

       case(besi1_)
          result = f_besi1(result)

       case(besin_)
          tmp(1)=stack(stackn);stackn=stackn-1
          if (aint(tmp(1)) == tmp(1)) then
             result = f_besin(tmp(1), result)
          else
             call s_seterror('besin(n, x), n not integer')
          end if

       case(besk0_)
          result = f_besk0(result)

       case(besk1_)
          result = f_besk1(result)

       case(beskn_)
          tmp(1)=stack(stackn);stackn=stackn-1
          if (aint(tmp(1)) == tmp(1)) then
             result = f_beskn(tmp(1), result)
          else
             call s_seterror('besjk(n, x), n not integer')
          end if

       case(ierfc_)
          if (abs(result-1.0d0) < 1.0d0) then
             result = f_ierfc(result)
          else
             call s_seterror('ierfc(x), 0 < x < 2')
          end if

       case(fresc_)
          result = f_fresc(result)

       case(fress_)
          result = f_fress(result)

       case(expi_)
          result = f_expi(result)

       case(sini_)
          result = f_sini(result)

       case(cosi_)
          result = f_cosi(result)

       case(logi_)
          if (result > 0.0d0) then
             result = f_logi(result)
          else
             call s_seterror('logi(x), x > 0')
          end if

       case(ierf_)
          if (abs(result) < 1.0d0) then
             result = f_ierf(result)
          else
             call s_seterror('ierf(x), -1 < x < 1')
          end if

       case(elle_)
          if (abs(result) < 1.0d0) then
             result = f_elle(result)
          else
             call s_seterror('elle(x), -1 < x < 1')
          end if

       case(ellk_)
          if (abs(result) < 1.0d0) then
             result = f_ellk(result)
          else
             call s_seterror('ellk(x), -1 < x < 1')
          end if

       case(ielle_)
          tmp(1)=stack(stackn);stackn=stackn-1
          if (abs(tmp(1)) < 1.0d0 .and. abs(result) < pi2) then
             result = f_ielle(tmp(1), result)
          else
             call s_seterror('ielle(x, phi), -1 < x < 1, -pi/2 < phi < pi/2')
          end if

       case(iellf_)
          tmp(1)=stack(stackn);stackn=stackn-1
          if (abs(tmp(1)) < 1.0d0 .and. abs(result) < pi2) then
             result = f_iellf(tmp(1), result)
          else
             call s_seterror('iellf(x, phi), -1 < x < 1, -pi/2 < phi < pi/2')
          end if

       end select
    end do

    if (stackn > 1) call s_seterror('Evaluation stack not empty')

    if (present(errvals)) then
       errvals = tmperrvals
       deallocate(tmperrvals)
    end if
  end subroutine s_evaluate

  subroutine s_variableerr(funcstack, number, varvals, errval, i)
    type(cachefn), intent(in) :: funcstack
    real(fdp), intent(inout) :: number
    real(fdp), dimension(:), intent(in) :: varvals
    real(fdp), intent(in) :: errval
    integer, intent(in) :: i
    real(fdp), dimension(size(varvals)) :: errvals
    real(fdp) :: numplus, numminus, rangeplus, rangeminus
    integer :: funcstackn

    number = 0.0d0

    if (error == 'OK' .and. size(varvals) > 0) then
       !
       ! Compute actual error
       !

       errvals = 0.0d0
       errvals(i) = errval
       funcstackn = 0

       call s_evaluate(funcstack, funcstackn, varvals, numplus, errvals)

       rangeplus = errvals(i)
       errvals(i) = -errval
       funcstackn = 0

       call s_evaluate(funcstack, funcstackn, varvals, numminus, errvals)

       rangeminus = errvals(i)

       number = merge(0.5d0*(numplus*rangeminus**2-numminus*rangeplus**2-&
            number*(rangeminus+rangeplus)**2)/rangeplus/rangeminus, &
            0.0d0, (rangeminus /= 0) .and. (rangeplus /= 0))
    end if
  end subroutine s_variableerr

  subroutine s_createfn(func, vars, funcnumber, t_error)
    character(len=*), intent(in) :: func, vars
    integer, intent(out) :: funcnumber
    character(len=*), intent(out), optional :: t_error
    type(cachefn), pointer :: funcstack
    type(cachefn), dimension(:), pointer :: tmpcachedfuncs
    integer :: funcstackn, i, j, n, lentmpfunc, components
    character(len=len(func)) :: tmpfunc

    funcnumber = 0
    error = 'OK'

    call s_firstrun()

    if (trim(vars) /= "") call s_allocatevariables(vars)

    if (error == 'OK') then
       n = size(cachedfuncs)

       do i=1, n
          if (.not. associated(cachedfuncs(i)%numbers)) exit
       end do

       if (i > n) then
          allocate(tmpcachedfuncs(i+10))

          do j = 1, n
             tmpcachedfuncs(j)%commands => cachedfuncs(j)%commands
             tmpcachedfuncs(j)%numbers => cachedfuncs(j)%numbers
          end do

          do j = n+1, size(tmpcachedfuncs)
             nullify(tmpcachedfuncs(j)%commands, tmpcachedfuncs(j)%numbers)
          end do

          deallocate(cachedfuncs)
          cachedfuncs => tmpcachedfuncs
       end if

       funcnumber = i
       funcstack => cachedfuncs(i)
       funcstackn = 0
       tmpfunc = f_clean(func)
       lentmpfunc = len_trim(tmpfunc)

       allocate(funcstack%commands(10), funcstack%numbers(10))

       funcstack%commands = end_
       funcstack%numbers = 0.0d0
       funcstackn = 0
       components = 0

       call s_parsefn(tmpfunc, lentmpfunc, funcstack, funcstackn, components)

!       if (components > 1) call s_seterror('Multiple component function')
    end if

    if (trim(vars) /= "") call s_deallocatevariables()

    if (present(t_error)) t_error = error
  end subroutine s_createfn

  subroutine s_destroyfn(funcnumber)
    integer, intent(in) :: funcnumber

    call s_firstrun()

    if (funcnumber > size(cachedfuncs)) then
       call s_seterror('Function does not exist')
    else if (associated(cachedfuncs(funcnumber)%commands)) then
       deallocate(cachedfuncs(funcnumber)%numbers, cachedfuncs(funcnumber)%commands)
    else
       call s_seterror('Function does not exist')
    end if
  end subroutine s_destroyfn

  subroutine s_evaluatefn(funcnumber, varvals, number, t_error)
    integer, intent(in) :: funcnumber
    real(fdp), dimension(:), intent(in) :: varvals
    real(fdp), intent(out) :: number
    character(len=*), intent(out), optional :: t_error
    type(cachefn) :: funcstack
    integer :: funcstackn

    error = 'OK'

    call s_firstrun()

    number = 0.0d0
    if (present(t_error)) t_error = ''

    if (funcnumber > size(cachedfuncs) .or. funcnumber < 1) then
       call s_seterror('Function does not exist')
    else if (associated(cachedfuncs(funcnumber)%commands)) then
       funcstack = cachedfuncs(funcnumber)
       funcstackn = 0

       call s_evaluate(funcstack, funcstackn, varvals, number)
    else
       call s_seterror('Function does not exist')
    end if

    if (present(t_error)) t_error = error
  end subroutine s_evaluatefn

  subroutine s_evaluateerr(funcnumber, varvals, number, t_error)
    integer, intent(in) :: funcnumber
    real(fdp), dimension(:), intent(in) :: varvals
    real(fdp), intent(out) :: number
    character(len=*), intent(out), optional :: t_error
    type(cachefn) :: funcstack
    real(fdp), dimension(size(varvals)) :: errvals
    real(fdp) :: num, halfnumber, tmp
    integer :: funcstackn, i

    error = 'OK'

    call s_firstrun()

    number = 0.0d0
    if (present(t_error)) t_error = ''
    halfnumber = 0.0d0

    if (funcnumber > size(cachedfuncs) .or. funcnumber < 1) then
       call s_seterror('Function does not exist')
    else if (associated(cachedfuncs(funcnumber)%commands)) then
       funcstack = cachedfuncs(funcnumber)
       errvals = 0.0d0
       funcstackn = 0

       call s_evaluate(funcstack, funcstackn, varvals, num, errvals)

       do i = 1, size(varvals)
          tmp = num

          call s_variableerr(funcstack, tmp, varvals, 1.0d0, i)

          number = number + abs(tmp)

          !
          ! Compute the error with half the range
          !

          tmp = num

          call s_variableerr(funcstack, tmp, varvals, 0.5d0, i)

          halfnumber = halfnumber + abs(tmp)
       end do

       !
       ! Return correct precision error
       !

       halfnumber = 2.0d0*halfnumber

       call s_errprecision(number)
       call s_errprecision(halfnumber)

       !
       ! Verify error value is consistent
       !

       if (halfnumber /= number) call s_seterror('Error value not stable')
    else
       call s_seterror('Function does not exist')
    end if

    if (present(t_error)) t_error = error
  end subroutine s_evaluateerr

end module USR_FunctionParser

! ========================================================
!
!  Include some stubs to allow non-module library linking
!
! ========================================================

!  subroutine createfn(func, vars, funcnumber, t_error)
!    use interpreter
!
!    character(*) :: func, vars
!    integer :: funcnumber
!    character(*), optional :: t_error
!
!    call s_createfn(func, vars, funcnumber, t_error)
!  end subroutine createfn
!
!
!
!  subroutine destroyfn(funcnumber)
!    use interpreter
!
!    integer :: funcnumber
!
!    call s_destroyfn(funcnumber)
!  end subroutine destroyfn
!
!
!
!  subroutine evaluatefn(funcnumber, varvals, number, t_error)
!    use interpreter
!
!    integer :: funcnumber
!    real(8), dimension(:) :: varvals
!    real(8) :: number
!    character(*), optional :: t_error
!
!    call s_evaluatefn(funcnumber, varvals, number, t_error)
!  end subroutine evaluatefn
!
!
!
!  subroutine evaluateerr(funcnumber, varvals, number, t_error)
!    use interpreter
!
!    integer :: funcnumber
!    real(8), dimension(:) :: varvals
!    real(8) :: number
!    character(*),optional :: t_error
!
!    call s_evaluateerr(funcnumber, varvals, number, t_error)
!  end subroutine evaluateerr

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


!---------------------------------------------------------------------!
! OWNER: Ithaca Combustion Enterprise, LLC                            !
! COPYRIGHT: © 2012, Ithaca Combustion Enterprise, LLC                !
! LICENSE: BSD 3-Clause License (The complete text of the license can !
!  be found in the `LICENSE-ICE.txt' file included in the ISAT-CK7    !
!  source directory.)                                                 !
!---------------------------------------------------------------------!

subroutine ell_pt_dist( n, c, gg, p, s )

!  Given the ellipsoid E (centered at c) and the point p, return in s
!  the relative distance from c to the boundary of E along the ray c-p
!  (relative to the distance |p-c|).
!  Let b denote the intersection of the boundary of E and the ray c-p.
!  Then s is:     s = |b-c| / |p-c|.
!
!  E is given by { x | norm(G^T * (x-c) ) <=1 ),  where G is an
!  n x n lower triangular matrix.  The array gg contains
!  the matrix G in packed format.

!  Errors: s = -1 is returned if p = c
!          s = -2 is returned if G is not positive definite

!  S.B. Pope  1/1/06

implicit none

integer, parameter      :: k_dp = kind(1.d0)
integer, intent(in)     :: n
real(k_dp), intent(in)  :: c(n), gg((n*(n+1))/2), p(n)
real(k_dp), intent(out) :: s

real(k_dp) :: y(n), ysq

y = p-c   !   form  y = G^T * [p-c]
call dtpmv( 'L', 'T', 'N', n, gg, y, 1 )

ysq = sum(y*y)

if( ysq > 0.d0 ) then
   s = 1.d0 / sqrt( ysq )   !  s successfully evaluated
elseif( sum( (p-c)*(p-c) ) == 0.d0 ) then
   s = -1.d0                !  p=c
else
   s = -2.d0                !  G is not positive definite
endif

return
end subroutine ell_pt_dist


text:
C
C
CC##       VC01A          25/10/74
CC NAME VC01A(R)                  CHECK
      SUBROUTINE VC01A(X,Y,W,Z,N,A,B,C,G,H,L,M,U)
C###### 25/10/74 LAST LIBRARY UPDATE
C     TAKEN FROM HARWELL SUBROUTINE LIBRARY
C     COMMON U HAS BEEN SUPPRESSED AND ADDED IN CALLING SEQUENCE
C
      DIMENSION X(1),Y(1),W(1),A(1),B(1),C(1),G(1),H(1),L(1)
      DIMENSION U(1),Z(1)
    4 FORMAT(3X,I5,5E16.6,6X,I5)
      LP=0
      EJ=0.0
      FJ=0.0
      GJ=0.0
      DO 1 I=1,N
      IV=I+N
      U(I)=SQRT(W(I))
      Z(I)=Y(I)*U(I)
      EJ=EJ+X(I)*U(I)**2
      FJ=FJ+U(I)*Z(I)
      GJ=GJ+U(I)**2
      U(IV)=0.0
    1 CONTINUE
      A(1)=EJ/GJ
      B(1)=0.0
      C(1)=FJ/GJ
      G(1)=1.0/GJ
      DO 2 J=1,M
      EJ=0.0
      FJ=0.0
      GJ=0.0
      H1=0.0
      L2=0
      DO 3 I=1,N
      Z(I)= Z(I)-C(J) *U(I)
      IV=I+N
      WS=(X(I)-A(J))*U(I)-B(J)*U(IV)
      U(IV)=U(I)
      U(I)=WS
      EJ=EJ+X(I)*U(I)**2
      FJ=FJ+U(I)*Z(I)
      GJ=GJ+U(I)**2
      H1=H1+Z(I)**2
      IF(1-I)13,3,3
   13 L2=(SIGN(1.0,-Z(I)*Z(I-1))+1.0)/2.0+L2
    3 CONTINUE
      A(J+1)=EJ/GJ
      B(J+1)=GJ*G(J)
      C(J+1)=FJ/GJ
      G(J+1)=1.0/GJ
      H(J)=H1
      L(J)=L2
    2 CONTINUE
    5 M1=M+1
      IF(LP.GT.0)WRITE(LP,20)
   20 FORMAT(1H1,8X,4HX(I),12X,4HY(I),12X,4HZ(I),//)
      H1=0.0
      L2=0
      DO 10 I=1,N
      W2=0.0
      YY=C(M+1)
      DO 101 J=1,M
      J1=M-J+1
      W1=W2
      W2=YY
      YY=C(J1)+(X(I)-A(J1))*W2-B(J1+1)*W1
  101 CONTINUE
      Z(I)=YY
      U(I)=(Z(I)-Y(I))
      H1=H1+U(I)*U(I)*W(I)
      IF(1-I)111,110,110
  111 L2=L2+(SIGN(1.0,-U(I)*U(I-1))+1.0)/2.0
  110 IF(LP.GT.0)WRITE(LP,11)X(I),Y(I),Z(I)
   10 CONTINUE
      H(M1)=H1
      L(M1)=L2
      IF(N.EQ.M1) GOTO 22
      H1=H1/(N-M1)
22    DO 21 J=1,M1
      G(J)=G(J)*H1
   21 CONTINUE
      IF(LP.LE.0)GO TO 40
      IF(LP.GT.0)WRITE(LP,7)
    7 FORMAT(1H1,61X,8HVARIANCE,10X,9HRESIDUALS)
      IF(LP.GT.0)WRITE(LP,6)
    6 FORMAT(7X,1HJ,7X,4HA(J),12X,4HB(J),12X,4HC(J),11X,7HOF C(J),6X,
     114HSUM OF SQUARES,3X,12HSIGN CHANGES)
      M1=M+1
      DO 24 I=1,M1
      II=I-1
   24 WRITE(LP,4)II,A(I),B(I),C(I),G(I),H(I),L(I)
   11 FORMAT(3E16.6)
   40 RETURN
      END

ext: f


!> Get the tree for a given halo
program get_tree
  use readtree
  use writetree
  use misc
  use types
  use flap, only : command_line_interface
  implicit none

  integer, parameter :: LEN_LIST = 10000
  !----------------------------------------
  ! cli params
  !----------------------------------------
  type(command_line_interface) :: cli
  integer, allocatable, dimension(:) :: param_halos
  character(len=256) :: param_output, param_tree_file
  integer :: param_verbosity

  !----------------------------------------
  ! loops
  !----------------------------------------
  integer :: i, j
  !----------------------------------------
  ! Tree paramters
  !----------------------------------------
  logical          :: stop_now
  integer(sp)  :: halo_id, BushID, step, level, hosthalo, hostsub
  integer(sp)  :: nbsub, nextsub, nb_of_fathers, nb_of_sons
  real(sp)     :: m, macc, px, py, pz, vx, vy, vz, Lx, Ly, Lz
  real(sp)     :: r, ra, rb, rc, ek, ep, et, spin

  integer(sp) :: prevStep, index

  real(sp), allocatable, dimension(:)    :: m_fathers
  integer(sp), allocatable, dimension(:) :: id_fathers, id_sons

  ! The list of halos to find
  type halo_list
     integer(sp), dimension(LEN_LIST) :: list = 0
  end type halo_list
  type(halo_list), dimension(:), allocatable :: halos ! FIXME: be able to have more
  integer(sp), dimension(LEN_LIST) :: order
  !----------------------------------------
  ! Get the parameters
  !----------------------------------------
  call cli%init(progname='get_tree_of_halo')
  call cli%add(switch='--halo-to-treat', help='Halo to treat', &
       act='store', nargs='*', def='-1')
  call cli%add(switch='--output', switch_ab='-o', help='Name of the output file', &
       act='store', def='data.out')
  call cli%add(switch='--verbose', help='Verbosity', &
       act='store', def='0')
  call cli%add(switch='--tree', help='Tree file', &
       act='store', def='/data33/dubois/H-AGN/MergerTreeHalo/HAGN/tree_rev.dat')

  call cli%get_varying(switch='--halo-to-treat', val=param_halos)
  call cli%get(switch='--output', val=param_output)
  call cli%get(switch='--verbose', val=param_verbosity)
  call cli%get(switch='--tree', val=param_tree_file)

  !----------------------------------------
  ! Read the reversed tree
  !----------------------------------------
  call read_tree(param_tree_file)

  ! allocate data
  allocate(halos(nsteps))
  do i = 1, nsteps-1
     halos(i)%list = 0
  end do
  call append(halos(nsteps), param_halos)

  if (param_verbosity >= 1) then
     write(*, *) 'Looking for halos', param_halos
  end if

  ! Opening output file
  if (param_verbosity >= 1) then
     write(*, *) 'Writing output in ', trim(param_output)
  end if
  call write_tree(param_output, nsteps, nb_of_halos, nb_of_subhalos, aexp, omega_t, age_univ)

  ! iter over all halos
  stop_now = .false.
  prevStep = 0
  do while (.not. stop_now)
     ! Read a part of the tree
     call iter_tree(stop_now, halo_id, BushID, step, level, hosthalo, hostsub, &
          nbsub, nextsub, nb_of_fathers, nb_of_sons, &
          m, macc, px, py, pz, vx, vy, vz, Lx, Ly, Lz, &
          r, ra, rb, rc, ek, ep, et, spin)

     if ((halo_id == param_halos(1)) .and. (step == nsteps)) then
        write(*, *) 'step=', step, 'level=', level, 'm=', m, '#fathers=', nb_of_fathers, '#sons=', nb_of_sons
     end if
     ! at a new step, reorder the halos to find (see below)
     if (step - prevStep /= 0) then
        prevStep = step
        ! sort the halos at the step
        call quick_sort(halos(step)%list, order)

        if (param_verbosity >= 2) then
           write(*, *) ''
           write(*, *)    'step:', step, 'z=', 1/aexp(nsteps-step+1) - 1
           if ((param_verbosity >= 3)) then
              block
                integer :: from, to, s
                s = sizeOf(halos(step))

                from = size(halos(step)%list) - s + 1
                to   = size(halos(step)%list)
                write(*, '(a,1000i8)') 'halos:', halos(step)%list(from:to)
              end block
           end if
        end if

     end if

     ! Try to find the current halo in the halo list
     ! if found, copy it in a new tree
     index = indexOf(halo_id, halos(step)%list)
     if (index > 0) then
        allocate(id_fathers(nb_of_fathers), m_fathers(nb_of_fathers))
        allocate(id_sons(nb_of_sons))

        ! add the parents to the next halos to find
        if (nb_of_fathers > 0) then
           call get_fathers(nb_of_fathers, id_fathers, m_fathers)

           if (nb_of_sons > 0) then
              call get_sons(nb_of_sons, id_sons)
           end if
           if (step-1 >= 1) then
              ! Append the fathers to the halos to treat next time
              call append(halos(step-1), id_fathers)
           end if
        end if

        !----------------------------------------
        ! Write the output
        !----------------------------------------
        call write_iter_tree(.false., halo_id, BushID, step, level, hosthalo, hostsub, &
             nbsub, nextsub, nb_of_fathers, nb_of_sons, &
             m, macc, px, py, pz, vx, vy, vz, Lx, Ly, Lz, &
             r, ra, rb, rc, ek, ep, et, spin, id_fathers, m_fathers, id_sons)

        deallocate(id_fathers, m_fathers, id_sons)
     end if
  end do

  ! Write stop as halo with id == -1
  call write_iter_tree(.true., halo_id, BushID, step, level, hosthalo, hostsub, &
                   nbsub, nextsub, nb_of_fathers, nb_of_sons, &
                   m, macc, px, py, pz, vx, vy, vz, Lx, Ly, Lz, &
                   r, ra, rb, rc, ek, ep, et, spin)

contains
  subroutine append (self, elements)
    type(halo_list), intent(inout) :: self
    integer(sp), intent(in), dimension(:)  :: elements
    integer(sp), dimension(size(elements)) :: tmp_elements, u_tmp_elements, order_els
    integer(sp), dimension(size(self%list)) :: tmp_list, order_list

    integer :: i, from, to, j

    u_tmp_elements = 0
    tmp_elements = elements
    tmp_list     = self%list
    call quick_sort(tmp_elements, order_els)
    call quick_sort(tmp_list, order_list)

    !----------------------------------------
    ! Remove elements already in list
    !----------------------------------------
    j = 0
    do i = 1, size(tmp_elements)
       if (indexOf(tmp_elements(i), tmp_list) <= 0) then
          j = j + 1
          u_tmp_elements(j) = tmp_elements(i)
       else
          ! print*, tmp_elements(i), 'already in list'
       end if
    end do

    ! Reorder array
    call quick_sort(u_tmp_elements, order_els)
    !----------------------------------------
    ! Find first non null index
    !----------------------------------------
    from = 0
    to = size(u_tmp_elements)
    ! Remove 0s from tmp_elements
    do i = 1, to
       if ((u_tmp_elements(i) > 0) .and. from == 0) then
          from = i
       end if
    end do

    ! find first null index in self
    i = sizeOf(self)+1

    if (from > 0) then
       ! TODO: check no overflow
       self%list(i:i+to-from) = u_tmp_elements(from:to)
       ! print*, 'WE:', self%list(1:i+to-from)
    end if
  end subroutine append

  !> Return the index of the last non null element
  function sizeOf(self)
    type(halo_list), intent(in) :: self
    integer :: sizeOf

    integer :: j

    sizeOf = 0

    do j = 1, size(self%list)
       if (self%list(j) > 0) sizeOf = sizeOf + 1
    end do
  end function sizeOf

end program get_tree


text:
module problem_tagging_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  public

contains

  ! This is a template routine for users to set their own tags based on the state.
  ! It will be overwritten by having a copy of this file in the user's problem setup.

  subroutine set_problem_tags(lo, hi, tag,tagl1,tagl2,tagl3,tagh1,tagh2,tagh3, &
                              state,state_l1,state_l2,state_l3,state_h1,state_h2,state_h3,&
                              set,clear,&
                              dx,problo,time,level) &
                              bind(C, name="set_problem_tags")

    use meth_params_module, only: NVAR

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         ,intent(in   ) :: lo(3),hi(3)
    integer         ,intent(in   ) :: state_l1,state_l2,state_l3, &
                                      state_h1,state_h2,state_h3
    integer         ,intent(in   ) :: tagl1,tagl2,tagl3,tagh1,tagh2,tagh3
    real(rt)        ,intent(in   ) :: state(state_l1:state_h1, &
                                      state_l2:state_h2, &
                                      state_l3:state_h3,NVAR)
    integer         ,intent(inout) :: tag(tagl1:tagh1,tagl2:tagh2,tagl3:tagh3)
    real(rt)        ,intent(in   ) :: problo(3),dx(3),time
    integer         ,intent(in   ) :: level,set,clear

  end subroutine set_problem_tags

end module problem_tagging_module

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro