aboutsummaryrefslogtreecommitdiff
path: root/Source/Initialization/WarpX_parser.F90
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Initialization/WarpX_parser.F90')
-rw-r--r--Source/Initialization/WarpX_parser.F901365
1 files changed, 1365 insertions, 0 deletions
diff --git a/Source/Initialization/WarpX_parser.F90 b/Source/Initialization/WarpX_parser.F90
new file mode 100644
index 000000000..4b95c5ac7
--- /dev/null
+++ b/Source/Initialization/WarpX_parser.F90
@@ -0,0 +1,1365 @@
+! MODULE mod_interpret
+! Parser to read and evaluate math expressions
+
+module mod_interpret
+use iso_c_binding
+use amrex_fort_module, only : amrex_real
+use amrex_error_module, only : amrex_error
+implicit none
+INTEGER, parameter :: plus =1, &
+ minus =2, &
+ multiply =3, &
+ divide =4, &
+ power =5, &
+ greaterthan=6, &
+ lessthan =7, &
+ exponential=8, &
+ logarithm =9, &
+ sine =10, &
+ cosine =11, &
+ tangent =12, &
+ square_root=13, &
+ arccos =14, &
+ arcsin =15, &
+ arctan =16, &
+ sinhyp =17, &
+ coshyp =18, &
+ tanhyp =19, &
+ logten =20
+
+CHARACTER(5), DIMENSION(0:20) :: coper
+INTEGER, parameter :: w_parenthesis = 1, &
+ w_number = 2, &
+ w_operator = 3, &
+ w_intrinsic = 4, &
+ w_variable = 5
+TYPE :: operation_type
+ INTEGER :: a,b,op
+END TYPE operation_type
+
+TYPE :: res_type
+ INTEGER :: nb_op
+ REAL(amrex_real) :: value
+ TYPE(res_type), DIMENSION(:), POINTER :: res
+ LOGICAL, DIMENSION(:), POINTER :: l_res, l_res_do
+ TYPE(operation_type), DIMENSION(:), POINTER :: operation
+ LOGICAL :: a_res,a_l_res,a_l_res_do,a_operation
+END TYPE res_type
+
+TYPE :: res_type_array
+ INTEGER :: nb_op
+ REAL(amrex_real), DIMENSION(:), POINTER :: value
+ TYPE(res_type_array), DIMENSION(:), POINTER :: res
+ LOGICAL, DIMENSION(:), POINTER :: l_res, l_res_do
+ TYPE(operation_type), DIMENSION(:), POINTER :: operation
+ LOGICAL :: a_res,a_l_res,a_l_res_do,a_operation
+END TYPE res_type_array
+
+TYPE(res_type_array) :: res_array
+
+INTEGER :: nb_res=0
+TYPE(res_type), DIMENSION(10) :: table_of_res
+
+contains
+
+RECURSIVE subroutine del_res(resin)
+implicit none
+TYPE(res_type), INTENT(INOUT) :: resin
+INTEGER :: j
+ IF(resin%a_res) then
+ do j = LBOUND(resin%res,1),UBOUND(resin%res,1)
+ call del_res(resin%res(j))
+ end do
+ DEALLOCATE(resin%res)
+ resin%a_res=.FALSE.
+ END if
+ IF(resin%a_l_res) then
+ DEALLOCATE(resin%l_res)
+ resin%a_l_res=.FALSE.
+ END if
+ IF(resin%a_l_res_do) then
+ DEALLOCATE(resin%l_res_do)
+ resin%a_l_res_do=.FALSE.
+ END if
+ IF(resin%a_operation) then
+ DEALLOCATE(resin%operation)
+ resin%a_operation=.FALSE.
+ END if
+return
+end subroutine del_res
+
+RECURSIVE subroutine del_res_array(resin)
+implicit none
+TYPE(res_type_array), INTENT(INOUT) :: resin
+INTEGER :: j
+ IF(resin%a_res) then
+ do j = LBOUND(resin%res,1),UBOUND(resin%res,1)
+ call del_res_array(resin%res(j))
+ end do
+ DEALLOCATE(resin%res)
+ resin%a_res=.FALSE.
+ END if
+ IF(resin%a_l_res) then
+ DEALLOCATE(resin%l_res)
+ resin%a_l_res=.FALSE.
+ END if
+ IF(resin%a_l_res_do) then
+ DEALLOCATE(resin%l_res_do)
+ resin%a_l_res_do=.FALSE.
+ END if
+ IF(resin%a_operation) then
+ DEALLOCATE(resin%operation)
+ resin%a_operation=.FALSE.
+ END if
+ if (associated(resin%value)) then
+ nullify(resin%value)
+ end if
+return
+end subroutine del_res_array
+
+subroutine init_interpret()
+implicit none
+
+coper(0 )='###'
+coper(plus )='+'
+coper(minus )='-'
+coper(multiply )='*'
+coper(divide )='/'
+coper(power )='**'
+coper(square_root)='sqrt'
+coper(exponential)='exp'
+coper(logarithm )='log'
+coper(sine )='sin'
+coper(cosine )='cos'
+coper(tangent )='tan'
+coper(arccos )='acos'
+coper(arcsin )='asin'
+coper(arctan )='atan'
+coper(sinhyp )='sinh'
+coper(coshyp )='cosh'
+coper(tanhyp )='tanh'
+coper(logten )='log10'
+coper(greaterthan)='>'
+coper(lessthan )='<'
+
+return
+end subroutine init_interpret
+
+recursive function calc_res(res,list_var,root) RESULT(calc_res_r)
+implicit none
+REAL(amrex_real) :: calc_res_r
+TYPE(res_type), INTENT(IN OUT) :: res
+REAL(amrex_real), DIMENSION(:), OPTIONAL :: list_var
+LOGICAL, INTENT(IN), OPTIONAL :: root
+
+INTEGER :: j,ia,ib,op
+REAL(amrex_real) :: a,b
+
+do j = 1, res%nb_op
+ ia = res%operation(j)%a
+ ib = res%operation(j)%b
+ op = res%operation(j)%op
+ IF(op<0) then
+ IF(op>-100) res%res(ia)%value=list_var(-op)
+ else
+ IF(res%l_res(ia).and.res%l_res_do(ia)) then
+ IF(PRESENT(list_var)) then
+ a = calc_res(res=res%res(ia),list_var=list_var,root=.false.)
+ else
+ a = calc_res(res=res%res(ia),root=.false.)
+ END if
+ res%l_res_do(ia)=.false.
+ else
+ a = res%res(ia)%value
+ END if
+ IF(res%l_res(ib).and.res%l_res_do(ib)) then
+ IF(PRESENT(list_var)) then
+ b = calc_res(res=res%res(ib),list_var=list_var,root=.false.)
+ else
+ b = calc_res(res=res%res(ib),root=.false.)
+ END if
+ ! This might prevent OMP threading
+ res%l_res_do(ib)=.false.
+ else
+ b = res%res(ib)%value
+ END if
+ ! This might prevent OMP threading
+ res%res(ia)%value=eval(a,b,op)
+ END if
+end do
+IF(res%nb_op==0) then
+ ia=1
+ calc_res_r = res%res(ia)%value
+else
+ calc_res_r = res%res(ia)%value
+END if
+
+IF(.not.PRESENT(root)) call calc_res_init(res)
+
+end function calc_res
+
+recursive subroutine calc_res_init(res)
+implicit none
+TYPE(res_type), INTENT(IN OUT) :: res
+
+INTEGER :: j,ia,ib,op
+
+do j = 1, res%nb_op
+ ia = res%operation(j)%a
+ ib = res%operation(j)%b
+ op = res%operation(j)%op
+ IF(op>=0) then
+ IF(res%l_res(ia).and..not.res%l_res_do(ia)) then
+ call calc_res_init(res%res(ia))
+ res%l_res_do(ia)=.true.
+ END if
+ IF(res%l_res(ib).and..not.res%l_res_do(ib)) then
+ call calc_res_init(res%res(ib))
+ res%l_res_do(ib)=.true.
+ END if
+ END if
+end do
+
+end subroutine calc_res_init
+
+recursive function calc_res_array(res,list_values,root) RESULT(calc_res_array_r)
+implicit none
+REAL(amrex_real), dimension(:), pointer :: calc_res_array_r
+TYPE(res_type_array), INTENT(IN OUT) :: res
+REAL(amrex_real), DIMENSION(:,:), OPTIONAL :: list_values
+LOGICAL, INTENT(IN), OPTIONAL :: root
+
+INTEGER :: j,ia,ib,op,n
+REAL(amrex_real),dimension(:),pointer :: a,b
+
+do j = 1, res%nb_op
+ ia = res%operation(j)%a
+ ib = res%operation(j)%b
+ op = res%operation(j)%op
+ IF(op<0) then
+ IF(op>-100) then
+ n = size(list_values,2)
+ allocate(res%res(ia)%value(n))
+ res%res(ia)%value=list_values(-op,:)
+ end if
+ else
+ IF(res%l_res(ia).and.res%l_res_do(ia)) then
+ IF(PRESENT(list_values)) then
+ a => calc_res_array(res=res%res(ia),list_values=list_values,root=.false.)
+ else
+ a => calc_res_array(res=res%res(ia),root=.false.)
+ END if
+ res%l_res_do(ia)=.false.
+ else
+ a => res%res(ia)%value
+ END if
+ IF(res%l_res(ib).and.res%l_res_do(ib)) then
+ IF(PRESENT(list_values)) then
+ b => calc_res_array(res=res%res(ib),list_values=list_values,root=.false.)
+ else
+ b => calc_res_array(res=res%res(ib),root=.false.)
+ END if
+ res%l_res_do(ib)=.false.
+ else
+ b => res%res(ib)%value
+ END if
+ res%res(ia)%value=>eval_array(a,b,op)
+ END if
+end do
+
+IF(res%nb_op==0) then
+ ia=1
+ calc_res_array_r => res%res(ia)%value
+else
+ calc_res_array_r => res%res(ia)%value
+END if
+
+IF(.not.PRESENT(root)) call calc_res_array_init(res)
+
+end function calc_res_array
+
+recursive subroutine calc_res_array_init(res)
+implicit none
+TYPE(res_type_array), INTENT(IN OUT) :: res
+
+INTEGER :: j,ia,ib,op
+
+do j = 1, res%nb_op
+ ia = res%operation(j)%a
+ ib = res%operation(j)%b
+ op = res%operation(j)%op
+ IF(op>=0) then
+ IF(res%l_res(ia).and..not.res%l_res_do(ia)) then
+ call calc_res_array_init(res%res(ia))
+ res%l_res_do(ia)=.true.
+ END if
+ IF(res%l_res(ib).and..not.res%l_res_do(ib)) then
+ call calc_res_array_init(res%res(ib))
+ res%l_res_do(ib)=.true.
+ END if
+ END if
+end do
+
+end subroutine calc_res_array_init
+
+recursive subroutine eval_res(res,exprin,int_op,list_var,l_root)
+implicit none
+TYPE(res_type) :: res
+CHARACTER(LEN=*), INTENT(IN) :: exprin
+INTEGER, OPTIONAL, INTENT(IN) :: int_op
+CHARACTER(*), DIMENSION(:),INTENT(IN), OPTIONAL :: list_var
+LOGICAL, OPTIONAL :: l_root
+
+CHARACTER(LEN=LEN(exprin)) :: expr_old
+CHARACTER(LEN=LEN(exprin)+2) :: expr
+INTEGER, parameter :: nop_tot=25
+INTEGER :: ln, i, j,ja
+INTEGER :: jtot,jres,jop
+INTEGER, DIMENSION(0:nop_tot) :: next_l,next_r,op,what
+LOGICAL, ALLOCATABLE, DIMENSION(:) :: l_res
+
+IF(.not.PRESENT(l_root)) call del_res(res)
+
+res%a_res=.false.
+res%a_l_res=.false.
+res%a_l_res_do=.false.
+res%a_operation=.false.
+
+op=0
+what=0
+
+expr_old = TRIM(ADJUSTL(C_replace(exprin,' ','')))
+IF(expr_old(1:1)=='-') then
+ expr='0'//expr_old
+else
+ expr='0+'//expr_old
+END if
+ln = LEN(TRIM(ADJUSTL(expr)))
+next_r(0) = 0
+i = 1
+j = 1
+jres = 0
+do WHILE(i<ln+1)
+ IF(PRESENT(list_var)) then
+ call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j),list_var)
+ else
+ call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j))
+ END if
+ next_r(j) = next_r(j)+next_r(j-1)
+ next_l(j) = next_r(j-1)+1+next_l(j)
+ i = next_r(j)+1
+ IF(what(j)/=w_operator) jres=jres+1
+ j = j+1
+end do
+jtot = j-1
+IF(PRESENT(int_op)) then
+ ALLOCATE(res%res(0:jres),res%operation(jres),res%l_res(0:jres),res%l_res_do(0:jres))
+ res%a_res=.true.
+ res%a_l_res=.true.
+ res%a_l_res_do=.true.
+ res%a_operation=.true.
+ res%res(0:jres)%a_res=.false.
+ res%res(0:jres)%a_l_res=.false.
+ res%res(0:jres)%a_l_res_do=.false.
+ res%res(0:jres)%a_operation=.false.
+ res%l_res(0)=.false.
+ res%l_res_do=.true.
+ res%res(0)%value=0.
+ res%nb_op=jres
+ res%operation(jres)%op=int_op
+ res%operation(jres)%a=1
+ res%operation(jres)%b=0
+ res%l_res(jres)=.false.
+ nullify(res%res(0)%res)
+ nullify(res%res(0)%l_res)
+ nullify(res%res(0)%l_res_do)
+ nullify(res%res(0)%operation)
+else
+ ALLOCATE(res%res(jres),res%operation(jres-1),res%l_res(jres),res%l_res_do(jres))
+ res%a_res=.true.
+ res%a_l_res=.true.
+ res%a_l_res_do=.true.
+ res%a_operation=.true.
+ res%res(1:jres)%a_res=.false.
+ res%res(1:jres)%a_l_res=.false.
+ res%res(1:jres)%a_l_res_do=.false.
+ res%res(1:jres)%a_operation=.false.
+ res%nb_op=jres-1
+ res%l_res_do=.true.
+END if
+ALLOCATE(l_res(jres))
+l_res(:) = .true.
+
+res%res(:)%value=0.
+
+jres = 0
+do j=1, jtot
+ select case (what(j))
+ case (w_number)
+ jres = jres+1
+ ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1))
+ res%res(jres)%a_res=.true.
+ res%res(jres)%a_l_res=.true.
+ res%res(jres)%a_l_res_do=.true.
+ res%res(jres)%a_operation=.true.
+ res%res(jres)%res(1)%a_res=.false.
+ res%res(jres)%res(1)%a_l_res=.false.
+ res%res(jres)%res(1)%a_l_res_do=.false.
+ res%res(jres)%res(1)%a_operation=.false.
+ res%res(jres)%nb_op=1
+ res%res(jres)%l_res(1) = .false.
+ res%res(jres)%l_res_do(1) = .true.
+ res%res(jres)%operation(1)%op = -100-op(j)
+ res%res(jres)%operation(1)%a = 1
+ res%res(jres)%operation(1)%b = 0
+ res%l_res(jres)=.true.
+ res%res(jres)%res(1)%value = char2real(expr(next_l(j):next_r(j)))
+ nullify(res%res(jres)%res(1)%res)
+ nullify(res%res(jres)%res(1)%l_res)
+ nullify(res%res(jres)%res(1)%l_res_do)
+ nullify(res%res(jres)%res(1)%operation)
+ case (w_parenthesis)
+ jres = jres+1
+ IF(PRESENT(list_var)) then
+ call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),list_var=list_var,l_root=.false.)
+ else
+ call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),l_root=.false.)
+ endif
+ res%l_res(jres)=.true.
+ case (w_intrinsic)
+ jres = jres+1
+ IF(PRESENT(list_var)) then
+ call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),list_var,l_root=.false.)
+ else
+ call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),l_root=.false.)
+ END if
+ res%l_res(jres)=.true.
+ case (w_variable)
+ jres = jres+1
+ ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1))
+ res%res(jres)%a_res=.true.
+ res%res(jres)%a_l_res=.true.
+ res%res(jres)%a_l_res_do=.true.
+ res%res(jres)%a_operation=.true.
+ res%res(jres)%res(1)%a_res=.false.
+ res%res(jres)%res(1)%a_l_res=.false.
+ res%res(jres)%res(1)%a_l_res_do=.false.
+ res%res(jres)%res(1)%a_operation=.false.
+ res%res(jres)%nb_op=1
+ res%res(jres)%l_res(1) = .false.
+ res%res(jres)%l_res_do(1) = .true.
+ res%res(jres)%operation(1)%op = -op(j)
+ res%res(jres)%operation(1)%a = 1
+ res%res(jres)%operation(1)%b = 0
+ res%l_res(jres)=.true.
+ nullify(res%res(jres)%res(1)%res)
+ nullify(res%res(jres)%res(1)%l_res)
+ nullify(res%res(jres)%res(1)%l_res_do)
+ nullify(res%res(jres)%res(1)%operation)
+ case (w_operator)
+ case default
+ end select
+end do
+jop = 1
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==power) then
+ res%operation(jop)%op = power
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END IF
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==divide) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==multiply) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==plus.or.op(j)==minus) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==greaterthan.or.op(j)==lessthan) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+
+DEALLOCATE(l_res)
+
+return
+end subroutine eval_res
+
+recursive subroutine eval_res_array(res,exprin,int_op,list_var,l_root)
+implicit none
+TYPE(res_type_array) :: res
+CHARACTER(LEN=*), INTENT(IN) :: exprin
+INTEGER, OPTIONAL, INTENT(IN) :: int_op
+CHARACTER(*), DIMENSION(:),INTENT(IN), OPTIONAL :: list_var
+LOGICAL, OPTIONAL :: l_root
+
+CHARACTER(LEN=LEN(exprin)) :: expr_old
+CHARACTER(LEN=LEN(exprin)+2) :: expr
+INTEGER, parameter :: nop_tot=25
+INTEGER :: ln, i, j,ja
+INTEGER :: jtot,jres,jop
+INTEGER, DIMENSION(0:nop_tot) :: next_l,next_r,op,what
+LOGICAL, ALLOCATABLE, DIMENSION(:) :: l_res
+
+IF(.not.PRESENT(l_root)) call del_res_array(res)
+
+res%a_res=.false.
+res%a_l_res=.false.
+res%a_l_res_do=.false.
+res%a_operation=.false.
+
+op=0
+what=0
+
+expr_old = TRIM(ADJUSTL(C_replace(exprin,' ','')))
+IF(expr_old(1:1)=='-') then
+ expr='0'//expr_old
+else
+ expr='0+'//expr_old
+END if
+ln = LEN(TRIM(ADJUSTL(expr)))
+next_r(0) = 0
+i = 1
+j = 1
+jres = 0
+do WHILE(i<ln+1)
+ IF(PRESENT(list_var)) then
+ call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j),list_var)
+ else
+ call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j))
+ END if
+ next_r(j) = next_r(j)+next_r(j-1)
+ next_l(j) = next_r(j-1)+1+next_l(j)
+ i = next_r(j)+1
+ IF(what(j)/=w_operator) jres=jres+1
+ j = j+1
+end do
+jtot = j-1
+IF(PRESENT(int_op)) then
+ ALLOCATE(res%res(0:jres),res%operation(jres),res%l_res(0:jres),res%l_res_do(0:jres))
+ res%a_res=.true.
+ res%a_l_res=.true.
+ res%a_l_res_do=.true.
+ res%a_operation=.true.
+ res%res(0:jres)%a_res=.false.
+ res%res(0:jres)%a_l_res=.false.
+ res%res(0:jres)%a_l_res_do=.false.
+ res%res(0:jres)%a_operation=.false.
+ res%l_res(0)=.false.
+ res%l_res_do=.true.
+ nullify(res%res(0)%value)
+ res%nb_op=jres
+ res%operation(jres)%op=int_op
+ res%operation(jres)%a=1
+ res%operation(jres)%b=0
+ res%l_res(jres)=.false.
+ nullify(res%res(0)%res)
+ nullify(res%res(0)%l_res)
+ nullify(res%res(0)%l_res_do)
+ nullify(res%res(0)%operation)
+else
+ ALLOCATE(res%res(jres),res%operation(jres-1),res%l_res(jres),res%l_res_do(jres))
+ res%a_res=.true.
+ res%a_l_res=.true.
+ res%a_l_res_do=.true.
+ res%a_operation=.true.
+ res%res(1:jres)%a_res=.false.
+ res%res(1:jres)%a_l_res=.false.
+ res%res(1:jres)%a_l_res_do=.false.
+ res%res(1:jres)%a_operation=.false.
+ res%nb_op=jres-1
+ res%l_res_do=.true.
+END if
+ALLOCATE(l_res(jres))
+l_res(:) = .true.
+
+jres = 0
+do j=1, jtot
+ select case (what(j))
+ case (w_number)
+ jres = jres+1
+ ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1))
+ res%res(jres)%a_res=.true.
+ res%res(jres)%a_l_res=.true.
+ res%res(jres)%a_l_res_do=.true.
+ res%res(jres)%a_operation=.true.
+ res%res(jres)%res(1)%a_res=.false.
+ res%res(jres)%res(1)%a_l_res=.false.
+ res%res(jres)%res(1)%a_l_res_do=.false.
+ res%res(jres)%res(1)%a_operation=.false.
+ res%res(jres)%nb_op=1
+ res%res(jres)%l_res(1) = .false.
+ res%res(jres)%l_res_do(1) = .true.
+ res%res(jres)%operation(1)%op = -100-op(j)
+ res%res(jres)%operation(1)%a = 1
+ res%res(jres)%operation(1)%b = 0
+ res%l_res(jres)=.true.
+ allocate(res%res(jres)%res(1)%value(1))
+ res%res(jres)%res(1)%value = char2real(expr(next_l(j):next_r(j)))
+ nullify(res%res(jres)%res(1)%res)
+ nullify(res%res(jres)%res(1)%l_res)
+ nullify(res%res(jres)%res(1)%l_res_do)
+ nullify(res%res(jres)%res(1)%operation)
+ case (w_parenthesis)
+ jres = jres+1
+ IF(PRESENT(list_var)) then
+ call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),list_var=list_var,l_root=.false.)
+ else
+ call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),l_root=.false.)
+ endif
+ res%l_res(jres)=.true.
+ case (w_intrinsic)
+ jres = jres+1
+ IF(PRESENT(list_var)) then
+ call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),list_var,l_root=.false.)
+ else
+ call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),l_root=.false.)
+ END if
+ res%l_res(jres)=.true.
+ case (w_variable)
+ jres = jres+1
+ ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1))
+ res%res(jres)%a_res=.true.
+ res%res(jres)%a_l_res=.true.
+ res%res(jres)%a_l_res_do=.true.
+ res%res(jres)%a_operation=.true.
+ res%res(jres)%res(1)%a_res=.false.
+ res%res(jres)%res(1)%a_l_res=.false.
+ res%res(jres)%res(1)%a_l_res_do=.false.
+ res%res(jres)%res(1)%a_operation=.false.
+ res%res(jres)%nb_op=1
+ res%res(jres)%l_res(1) = .false.
+ res%res(jres)%l_res_do(1) = .true.
+ res%res(jres)%operation(1)%op = -op(j)
+ res%res(jres)%operation(1)%a = 1
+ res%res(jres)%operation(1)%b = 0
+ res%l_res(jres)=.true.
+ nullify(res%res(jres)%res(1)%res)
+ nullify(res%res(jres)%res(1)%l_res)
+ nullify(res%res(jres)%res(1)%l_res_do)
+ nullify(res%res(jres)%res(1)%operation)
+ case (w_operator)
+ case default
+ end select
+end do
+jop = 1
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==power) then
+ res%operation(jop)%op = power
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END IF
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==divide) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==multiply) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==plus.or.op(j)==minus) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+do j=jtot,1,-1
+ IF(what(j)==w_operator) then
+ IF(op(j)==greaterthan.or.op(j)==lessthan) then
+ res%operation(jop)%op = op(j)
+ ja=j/2
+ do WHILE(.not.l_res(ja))
+ ja=ja-1
+ end do
+ res%operation(jop)%a = ja
+ res%operation(jop)%b = j/2+1
+ l_res(j/2+1)=.false.
+ jop = jop+1
+ END if
+ END if
+end do
+
+DEALLOCATE(l_res)
+
+return
+end subroutine eval_res_array
+
+subroutine what_is_next(exprin,what,next_l,next_r,op,list_var)
+implicit none
+CHARACTER(LEN=*), INTENT(IN) :: exprin
+INTEGER, INTENT(OUT) :: what,next_l,next_r,op
+CHARACTER(*), DIMENSION(:),INTENT(IN), optional :: list_var
+
+CHARACTER(LEN=LEN(exprin)) :: expr
+
+INTEGER :: ln, i, asc, istart,ln2,ierror,jv
+
+ln = LEN_TRIM(ADJUSTL(exprin))
+expr = TRIM(ADJUSTL(exprin))
+istart = LEN(exprin)-LEN(ADJUSTL(exprin))
+next_l = istart
+
+i = 1
+asc = IACHAR(expr(1:1))
+select case (asc)
+ ! parenthesis (
+ case (40)
+ what = w_parenthesis
+ next_r = find_close_bracket(expr(i:ln))
+ ! alphabet
+ case (65:90,97:122)
+ ln2 = ln-i+1
+ op = find_intrinsic(expr(i:ln),ln2,ierror)
+ IF(ierror==0) then
+ what = w_intrinsic
+ IF( expr(i+ln2:i+ln2) .ne.'(' ) then
+ WRITE(*,*) 'Error in intrinsic ',exprin,': missing parenthesis.'
+ stop
+ END if
+ next_l = next_l+ln2
+ next_r = ln2+find_close_bracket(expr(i+ln2:ln))
+ ELSEIF(PRESENT(list_var)) then
+ what = w_variable
+ jv = 0
+ op = 0
+ do WHILE(jv+1<=SIZE(list_var))
+ jv=jv+1
+ IF(expr(i:i+ln2-1)==TRIM(list_var(jv)(:))) then
+ op=jv
+ END if
+ end do
+ IF(op==0) then
+ WRITE(*,*) 'Error: variable ',exprin(i:i+ln2-1),' not allowed.'
+ stop
+ END if
+ next_l = next_l+ln2
+ next_r = ln2
+ else
+ WRITE(*,*) 'Error: intrinsic ',exprin(i:i+ln2-1),' does not exists.'
+ stop
+ END if
+ ! number
+ case (48:57)
+ what = w_number
+ next_r = next_l+find_number_end(expr(i:ln))
+ ! operator
+ case (42,43,45,47,60,62)
+ what = w_operator
+ op = find_operator(expr(i:i+1))
+ IF(op==power) THEN
+ next_r = next_l+2
+ else
+ next_r = next_l+1
+ END if
+ case default
+end select
+
+return
+end subroutine what_is_next
+
+function find_close_bracket(a)
+implicit none
+INTEGER :: find_close_bracket
+CHARACTER(*), INTENT(IN) :: a
+INTEGER :: level, i, ln
+level = 0
+ln=LEN(a)
+do i = 1, ln
+ select case (a(i:i))
+ case ('(')
+ level=level+1
+ case (')')
+ level=level-1
+ case default
+ end select
+ IF(level==0) exit
+end do
+IF(i>ln) then
+ WRITE(*,*) 'Error, no closing bracket in expression ',a
+ stop
+END if
+
+find_close_bracket = i
+
+return
+end function find_close_bracket
+
+function find_intrinsic(a,ln,ierror)
+implicit none
+INTEGER :: find_intrinsic
+CHARACTER(*), INTENT(IN) :: a
+CHARACTER(LEN=LEN(a)) :: ac
+INTEGER, INTENT(IN OUT) :: ln
+INTEGER, INTENT(OUT) :: ierror
+INTEGER :: i
+
+ierror=0
+do i = 1, ln
+ select case (IACHAR(a(i:i)))
+ case (65:90,97:122) ! alphabet
+ case default
+ exit
+ end select
+end do
+IF(i>ln+1) then
+ WRITE(*,*) 'Error in expression ',a
+ stop
+END if
+
+ln = i-1
+ac(1:ln) = a(1:ln)
+
+do i = 1, ln
+ IF(IACHAR(a(i:i))<97) ac(i:i)=ACHAR(IACHAR(a(i:i))+32)
+END do
+
+select case (ac(1:ln))
+ case ('sqrt')
+ find_intrinsic = square_root
+ case ('sin')
+ find_intrinsic = sine
+ case ('cos')
+ find_intrinsic = cosine
+ case ('tan')
+ find_intrinsic = tangent
+ case ('exp')
+ find_intrinsic = exponential
+ case ('log','ln')
+ find_intrinsic = logarithm
+ case ('log10')
+ find_intrinsic = logten
+ case ('asin')
+ find_intrinsic = arcsin
+ case ('acos')
+ find_intrinsic = arccos
+ case ('atan')
+ find_intrinsic = arctan
+ case ('sinh')
+ find_intrinsic = sinhyp
+ case ('cosh')
+ find_intrinsic = coshyp
+ case ('tanh')
+ find_intrinsic = tanhyp
+ case default
+ ierror=1
+end select
+return
+end function find_intrinsic
+
+function find_operator(a)
+implicit none
+INTEGER :: find_operator
+CHARACTER(*), INTENT(IN) :: a
+
+ select case (IACHAR(a(1:1)))
+ case (42) ! '*'
+ IF(IACHAR(a(2:2))==42) then
+ find_operator = power
+ else
+ find_operator = multiply
+ END if
+ case (43) ! '+'
+ find_operator = plus
+ case (45) ! '-'
+ find_operator = minus
+ case (47) ! '/'
+ find_operator = divide
+ case (60) ! '<'
+ find_operator = lessthan
+ case (62) ! '>'
+ find_operator = greaterthan
+ case default
+ WRITE(*,*) 'Error, operator ',a(1:1),' does not exist'
+ stop
+ end select
+
+ return
+end function find_operator
+
+function find_number_end(a)
+implicit none
+INTEGER :: find_number_end
+CHARACTER(*), INTENT(IN) :: a
+INTEGER :: i
+LOGICAL :: l_point,l_minus,l_e
+
+l_point=.FALSE.
+l_minus=.FALSE.
+l_e=.FALSE.
+
+do i = 1, len(a)
+ select case (IACHAR(a(i:i)))
+ case (48:57)
+ l_minus = .true.
+ case (46) ! .
+ IF(l_point) exit
+ l_point = .true.
+ l_minus = .true.
+ case (45) ! -
+ IF(l_minus) exit
+ l_minus = .true.
+ case (69,101) ! e
+ IF(l_e) exit
+ l_e = .true.
+ l_minus = .false.
+ case default
+ exit
+ end select
+end do
+
+find_number_end = i-1
+
+return
+end function find_number_end
+
+function char2real(a)
+implicit none
+REAL(amrex_real) :: char2real
+CHARACTER(*), INTENT(IN) :: a
+
+READ(a(:),*) char2real
+
+return
+end function char2real
+
+FUNCTION eval(a,b,op)
+implicit none
+REAL(amrex_real) :: eval
+REAL(amrex_real), INTENT(IN) :: a, b
+INTEGER, INTENT(IN) :: op
+
+select case (op)
+ case (plus)
+ eval=a+b
+ case (minus)
+ eval=a-b
+ case (multiply)
+ eval=a*b
+ case (divide)
+ eval=a/b
+ case (power)
+ eval=a**b
+ case (greaterthan)
+ eval=0.
+ if (a > b) eval=1.
+ case (lessthan)
+ eval=0.
+ if (a < b) eval=1.
+ case (square_root)
+ eval=SQRT(a)
+ case (exponential)
+ eval=EXP(a)
+ case (logarithm)
+ eval=LOG(a)
+ case (sine)
+ eval=SIN(a)
+ case (cosine)
+ eval=COS(a)
+ case (tangent)
+ eval=TAN(a)
+ case (logten)
+ eval=log10(a)
+ case (arcsin)
+ eval=asin(a)
+ case (arccos)
+ eval=acos(a)
+ case (arctan)
+ eval=atan(a)
+ case (sinhyp)
+ eval=sinh(a)
+ case (coshyp)
+ eval=cosh(a)
+ case (tanhyp)
+ eval=tanh(a)
+ case default
+ WRITE(*,*) 'Error in eval, wrong operator ',op
+ stop
+end select
+
+return
+END function eval
+
+
+FUNCTION eval_array(ai,bi,op)
+implicit none
+REAL(amrex_real), dimension(:), pointer :: eval_array
+REAL(amrex_real), INTENT(IN), dimension(:), pointer :: ai, bi
+REAL(amrex_real), dimension(:), pointer :: a, b
+INTEGER, INTENT(IN) :: op
+integer::i1,i2,n
+
+i1 = size(ai)
+i2 = size(bi)
+n = maxval((/i1,i2/))
+
+select case (op)
+ case (plus,minus,multiply,divide,power, greaterthan, lessthan)
+
+ if (i1<n) then
+ allocate(a(n))
+ a(:) = ai(1)
+ else
+ a => ai
+ end if
+ if (i2<n) then
+ allocate(b(n))
+ b(:) = bi(1)
+ else
+ b => bi
+ end if
+
+ case default
+ a => ai
+
+end select
+
+allocate(eval_array(n))
+eval_array=0.
+
+select case (op)
+ case (plus)
+ eval_array=a+b
+ case (minus)
+ eval_array=a-b
+ case (multiply)
+ eval_array=a*b
+ case (divide)
+ eval_array=a/b
+ case (power)
+ eval_array=a**b
+ case (square_root)
+ eval_array=SQRT(a)
+ case (exponential)
+ eval_array=EXP(a)
+ case (logarithm)
+ eval_array=LOG(a)
+ case (sine)
+ eval_array=SIN(a)
+ case (cosine)
+ eval_array=COS(a)
+ case (tangent)
+ eval_array=TAN(a)
+ case (logten)
+ eval_array=log10(a)
+ case (arcsin)
+ eval_array=asin(a)
+ case (arccos)
+ eval_array=acos(a)
+ case (arctan)
+ eval_array=atan(a)
+ case (sinhyp)
+ eval_array=sinh(a)
+ case (coshyp)
+ eval_array=cosh(a)
+ case (tanhyp)
+ eval_array=tanh(a)
+ case default
+ WRITE(*,*) 'Error in eval_array, wrong operator ',op
+ stop
+end select
+
+return
+END function eval_array
+
+function C_up2low(cline)
+CHARACTER(*) :: cline
+CHARACTER(LEN=LEN(cline)) :: C_up2low
+INTEGER :: ln, i, j
+ ln = LEN(cline)
+ do i = 1, ln
+ j = IACHAR(cline(i:i))
+ IF(j<91.and.j>64) then
+ C_up2low(i:i)=ACHAR(j+32)
+ else
+ C_up2low(i:i) = cline(i:i)
+ END if
+ END do
+end function C_up2low
+
+RECURSIVE FUNCTION C_REPLACE(ST,R1,R2) RESULT(C_RES)
+IMPLICIT NONE
+CHARACTER(LEN=*), INTENT(IN) :: ST
+CHARACTER(LEN=LEN(st)) :: c_res
+CHARACTER(LEN=*), INTENT(IN) :: R1
+CHARACTER(LEN=*), INTENT(IN) :: R2
+
+INTEGER :: i,ll
+
+ c_res = ST
+ i=INDEX(TRIM(st),r1)
+ IF(i==0) return
+ ll = LEN(r2)
+ IF(ll>0) c_res(i:i+ll-1) = r2
+ c_res(i+ll:) = c_replace(c_res(i+LEN(r1):),r1,r2)
+
+
+RETURN
+END FUNCTION C_REPLACE
+
+FUNCTION C_nboccur(ST,R1)
+IMPLICIT NONE
+CHARACTER(LEN=*), INTENT(IN) :: ST
+integer :: c_nboccur
+CHARACTER(LEN=*), INTENT(IN) :: R1
+
+INTEGER :: i,is
+
+ C_nboccur=0
+ is=0
+10 i = INDEX(st(is+1:),r1)
+ IF(i/=0) then
+ C_nboccur = C_nboccur+1
+ is=is+i+LEN(r1)
+ GO TO 10
+ END if
+
+RETURN
+END FUNCTION C_nboccur
+
+function stringlist(strin) RESULT(res)
+CHARACTER(*), INTENT(IN) :: strin
+INTEGER :: i1, nwords, i, i2
+
+CHARACTER(LEN(strin)+1) :: str
+CHARACTER(80), DIMENSION(:), pointer :: res
+
+str=ADJUSTL(strin)//' '
+i1=1
+nwords = 0
+do WHILE(TRIM(str(i1:))/='')
+ i2 = 1
+ i2 = INDEX(str(i1:),' ')
+ do WHILE(i2==1)
+ i1 = i1 + 1
+ i2 = INDEX(str(i1:),' ')
+ end do
+ i1=i1+i2
+ nwords = nwords + 1
+end do
+ALLOCATE(res(nwords))
+res(:) = ''
+i1=1
+do i = 1, nwords
+ i2 = INDEX(str(i1:),' ')
+ do WHILE(i2==1)
+ i1 = i1 + 1
+ i2 = INDEX(str(i1:),' ')
+ end do
+ res(i)(1:i2-1) = str(i1:i1+i2-1)
+ i1=i1+i2
+end do
+
+return
+end function stringlist
+
+SUBROUTINE csv2list(strin, strinlen, lofstr)
+CHARACTER(*), INTENT(IN) :: strin
+INTEGER, INTENT(IN) :: strinlen
+INTEGER :: isep, isep_rel, nwords, i, j
+INTEGER, DIMENSION(strinlen) :: word_start, word_end
+CHARACTER(strinlen) :: str
+CHARACTER(len=10), DIMENSION(:), allocatable, INTENT(inout) :: lofstr
+nwords = 1
+str = strin
+isep_rel = INDEX(str, ',')
+isep = isep_rel
+word_start(nwords) = 1
+DO WHILE(isep_rel > 0)
+ word_end(nwords) = isep-1
+ nwords = nwords + 1
+ word_start(nwords) = isep+1
+ str = strin(isep+1:strinlen)
+ isep_rel = INDEX(str, ',')
+ isep = isep_rel + isep
+ENDDO
+word_end(nwords) = strinlen
+
+ALLOCATE( lofstr(nwords) )
+
+DO i=1, nwords
+ lofstr(i) = ""
+ DO j = 1, word_end(i)-word_start(i)+1
+ lofstr(i)(j:j) = strin(word_start(i)+j-1:word_start(i)+j-1)
+ end DO
+ENDDO
+END SUBROUTINE csv2list
+end module mod_interpret
+
+! ---------------------
+! MODULE parser_wrapper
+! Wrappers to use the parser defined in MODULE mod_interpret.
+MODULE parser_wrapper
+USE iso_c_binding
+USE amrex_fort_module, only : amrex_real
+USE amrex_error_module, only : amrex_error
+USE mod_interpret, only : nb_res, table_of_res, calc_res, eval_res, csv2list
+
+IMPLICIT NONE
+
+CONTAINS
+
+! FUNCTION parser_initialize_function RESULT my_index_res
+! Initialize a res_type and assign it a unique identifier.
+! INPUTS:
+!> instr_func : CHAR* for a mathematical expression
+!> e.g. instr_func = "3*cos(x)+y".
+!> instr_var : CHAR* for a comma-separated list of variables in instr_func
+!> e.g. instr_var = "x,y".
+! OUTPUT:
+!> my_index_res : INT parser index. Necessary if this module is used for several
+!> parsers (e.g. one for the electron density, one for the ion
+!> density, one for the laser profile etc.).
+FUNCTION parser_initialize_function(instr_func, length_func, instr_var, length_var) RESULT(my_index_res) &
+ bind(c,name='parser_initialize_function')
+ INTEGER(c_int), INTENT(IN), VALUE :: length_func, length_var
+ CHARACTER(kind=c_char), INTENT(IN) :: instr_func(length_func)
+ CHARACTER(kind=c_char), INTENT(IN) :: instr_var(length_var)
+ INTEGER :: i
+ CHARACTER(len=10), DIMENSION(:), allocatable :: lofstr
+ CHARACTER(len=length_func) :: str_func
+ CHARACTER(len=length_var) :: str_var
+ REAL(amrex_real), DIMENSION(1) :: list_var
+ INTEGER :: my_index_res
+
+ nb_res = nb_res + 1
+ my_index_res = nb_res
+ IF (nb_res>10) THEN
+ call amrex_error('Parser error: cannot have more than 10 parsers.')
+ ENDIF
+
+ str_func = ""
+ DO i=1, length_func
+ str_func(i:i) = instr_func(i)
+ ENDDO
+
+ str_var = ""
+ DO i=1, length_var
+ str_var(i:i) = instr_var(i)
+ ENDDO
+
+ ! Convert variable list from csv string to list of strings.
+ CALL csv2list(str_var, len_trim(str_var), lofstr)
+
+ ! Initialize the res object
+ CALL eval_res(table_of_res(my_index_res), str_func, list_var=lofstr)
+
+END FUNCTION parser_initialize_function
+
+! FUNCTION parser_evaluate_function RESULT out
+! Evaluate parsed function
+! INPUTS:
+!> list_var : REAL* array of values for variables
+!> e.g. list_var = (/3.14_8,2._8/).
+!> nvar : INT number of variables
+!> e.g. nvar = 2.
+!> my_index_res : INT index of the res_type object in table table_of_res.
+! OUTPUT:
+!> out : REAL Result.
+FUNCTION parser_evaluate_function(list_var, nvar, my_index_res) result(out) &
+ bind(c,name='parser_evaluate_function')
+ INTEGER, VALUE, INTENT(IN) :: nvar, my_index_res
+ REAL(amrex_real), INTENT(IN) :: list_var(1:nvar)
+ REAL(amrex_real) :: out
+ ! Evaluate parsed function in table_of_res(my_index_res), of type res_type
+ out = calc_res(table_of_res(my_index_res),list_var)
+
+END FUNCTION parser_evaluate_function
+
+END MODULE parser_wrapper
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-20 14:54:21 +0000