data_suport.F90

Go to the documentation of this file.

module data_suport

  use f90_kind
  use call_utils
  use unit_h
  use strings
  use m_mrgrnk

  implicit none

  private

  !data types
  !
  !--- axes_t ----
  type axes_t
     real(rkind)   , allocatable  :: x(:)
     integer(ikind)               :: dim_x
     integer(ikind)               :: lbound
     logical                      :: is_lin=.false.
     logical                      :: is_log=.false.
     logical                      :: is_equi=.false.
     logical                      :: allocated
     integer                      :: last_access
  end type axes_t

  !
  !--- grid_t ----
  ! update this and all nesessary routines below for 5D
  type grid_t
     real(rkind)   , allocatable  :: f1d(:)        ! Values 1D
     real(rkind)   , allocatable  :: f2d(:,:)      ! Values 2D
     real(rkind)   , allocatable  :: f3d(:,:,:)    ! Values 3D
     real(rkind)   , allocatable  :: f4d(:,:,:,:)  ! Values 4D
     type (axes_t) , allocatable  :: axe(:)
     logical                      :: exist_f1d=.false.
     logical                      :: exist_f2d=.false.
     logical                      :: exist_f3d=.false.
     logical                      :: exist_f4d=.false.
     logical                      :: is_lin=.false.
     logical                      :: is_log=.false.
     logical                      :: with_warning=.false.
     character(len=skind)         :: name=""
     integer                      :: interpol_function = 0
     character(len=132)           :: result_label=""
     character(len=132)           :: result_unit=""
     character(len=132)           :: state_label=""
     integer                      :: ndim=0
  end type grid_t

  !
  !--- data_error ----
  type data_error_t
     integer(ikind)               :: ierr
     character (len=128)          :: cerr
  end type data_error_t


  ! update this and all nesessary routines below for 5D

  !! \todo Update this interface for 5D and up
  interface new_grid
     module procedure new_grid_1d,  &
          new_grid_2d,  &
          new_grid_3d,  &
          new_grid_4d
  end interface new_grid

  public :: grid_t, axes_t
  public :: interpol, new_grid, delete
  public :: set_option
  public :: data_error_t

contains


  !-------------------------------------------
  !------- new_grid_1 --------------------------
  !-------------------------------------------

  subroutine new_grid_1d(grid,f,x)
    implicit none
    real(rkind), intent(in) :: f(:)
    real(rkind), intent(in), optional :: x(:)
    type(grid_t),intent(out):: grid
    integer :: optargs

    ! error check
    optargs=0
    call assert(.not.grid%exist_f1d,"Error: grid%f1d is still allocated")
    call assert(.not.grid%exist_f2d,"Error: grid%f2d is still allocated")
    call assert(.not.grid%exist_f3d,"Error: grid%f3d is still allocated")
    call assert(.not.grid%exist_f4d,"Error: grid%f4d is still allocated")
    if(present(x)) optargs=optargs+1
    if(optargs.eq.1) then
       call assert(size(f)==size(x), "Error: size(f)" //size(f)// " /= =size(x) "//size(x))
       call assert(sorted(x(:)),"not sorted x")
       call assert(size(x(:))>=2,"not sorted x")
       allocate(grid%axe(1))
       allocate(grid%axe(1)%x(size(x)))
       grid%axe(1)%x(:) = x(:)
       grid%axe(1)%dim_x = size(x)
       grid%axe(1)%allocated = .true.
    endif
    allocate(grid%f1d(size(f,1)))
    grid%f1d(:)=f(:)
    grid%exist_f1d=.true.
  end subroutine new_grid_1d

  !-------------------------------------------
  !------- new_grid_2d --------------------------
  !-------------------------------------------

  subroutine new_grid_2d(grid,f,x,y)
    implicit none
    real(rkind), intent(in) :: f(:,:)
    real(rkind), intent(in), optional :: x(:)
    real(rkind), intent(in), optional :: y(:)
    type(grid_t),intent(out):: grid
    integer :: optargs

    ! error check
    optargs=0
    call assert(.not.grid%exist_f1d, "Error: grid%f1d is still allocated")
    call assert(.not.grid%exist_f2d, "Error: grid%f2d is still allocated")
    call assert(.not.grid%exist_f3d, "Error: grid%f3d is still allocated")
    call assert(.not.grid%exist_f4d, "Error: grid%f4d is still allocated")
    if(present(x)) optargs=optargs+1
    if(present(y)) optargs=optargs+1
    if(optargs.eq.2) then
       call assert(size(f,1)==size(x), "Error: size(f,1)" //size(f,1)// ' /= =size(x) '//size(x))
       call assert(size(f,2)==size(y), "Error: size(f,2)" //size(f,2)// ' /= =size(y) '//size(y))
       call assert(sorted(x(:)),"not sorted x")
       call assert(sorted(y(:)),"not sorted y")
       call assert(size(x(:))>=2,"not sorted x")
       call assert(size(y(:))>=2,"not sorted y")
       allocate(grid%axe(2))
       allocate(grid%axe(1)%x(size(x)))
       allocate(grid%axe(2)%x(size(y)))
       grid%axe(1)%x(:) = x(:)
       grid%axe(2)%x(:) = y(:)
       grid%axe(1)%dim_x = size(x)
       grid%axe(2)%dim_x = size(y)
       grid%axe(1)%allocated = .true.
       grid%axe(2)%allocated = .true.
    else
       call assert(optargs==0, "Error: wrong number of optionl arguments for new_grid_2d")
    endif
    allocate(grid%f2d(size(f,1),size(f,2)))
    grid%f2d(:,:)=f(:,:)
    grid%exist_f2d=.true.
  end subroutine new_grid_2d

  !-------------------------------------------
  !------- new_grid_3d --------------------------
  !-------------------------------------------

  subroutine new_grid_3d(grid,f,x,y,z)
    implicit none
    real(rkind), intent(in) :: f(:,:,:)
    real(rkind), intent(in), optional :: x(:)
    real(rkind), intent(in), optional :: y(:)
    real(rkind), intent(in), optional :: z(:)
    type(grid_t),intent(out):: grid
    integer :: optargs

    ! error check
    optargs=0
    call assert(.not.grid%exist_f1d, "Error: grid%f1d is still allocated")
    call assert(.not.grid%exist_f2d, "Error: grid%f2d is still allocated")
    call assert(.not.grid%exist_f3d, "Error: grid%f3d is still allocated")
    call assert(.not.grid%exist_f4d, "Error: grid%f4d is still allocated")
    if(present(x)) optargs=optargs+1
    if(present(y)) optargs=optargs+1
    if(present(z)) optargs=optargs+1
    if(optargs.eq.3) then
       call assert(size(f,1)==size(x), "Error: size(fx,1)" //size(f,1)// ' /= =size(x) '//size(x))
       call assert(size(f,2)==size(y), "Error: size(fx,2)" //size(f,2)// ' /= =size(y) '//size(y))
       call assert(size(f,3)==size(z), "Error: size(fx,3)" //size(f,3)// ' /= =size(z) '//size(z))
       call assert(sorted(x(:)),"not sorted x")
       call assert(sorted(y(:)),"not sorted y")
       call assert(sorted(y(:)),"not sorted z")
       call assert(size(x(:))>=2,"not sorted x")
       call assert(size(y(:))>=2,"not sorted y")
       call assert(size(z(:))>=2,"not sorted z")
       allocate(grid%axe(3))
       allocate(grid%axe(1)%x(size(x)))
       allocate(grid%axe(2)%x(size(y)))
       allocate(grid%axe(3)%x(size(z)))
       grid%axe(1)%x(:) = x(:)
       grid%axe(2)%x(:) = y(:)
       grid%axe(3)%x(:) = z(:)
       grid%axe(1)%dim_x = size(x)
       grid%axe(2)%dim_x = size(y)
       grid%axe(3)%dim_x = size(z)
       grid%axe(1)%allocated = .true.
       grid%axe(2)%allocated = .true.
       grid%axe(3)%allocated = .true.
    else
       call assert(optargs==0, "Error: wrong number of optionl arguments for new_grid_3d")
    endif
    allocate(grid%f3d(size(f,1),size(f,2),size(f,3)))
    grid%f3d(:,:,:)=f(:,:,:)
    grid%exist_f3d=.true.
  end subroutine new_grid_3d

  !-------------------------------------------
  !------- new_grid_4d --------------------------
  !-------------------------------------------

  subroutine new_grid_4d(grid,f,w,x,y,z)
    implicit none
    real(rkind), intent(in) :: f(:,:,:,:)
    real(rkind), intent(in), optional :: w(:)
    real(rkind), intent(in), optional :: x(:)
    real(rkind), intent(in), optional :: y(:)
    real(rkind), intent(in), optional :: z(:)
    type(grid_t),intent(out):: grid
    integer :: optargs

    ! error check
    optargs=0
    call assert(.not.grid%exist_f1d, "Error: grid%f1d is still allocated")
    call assert(.not.grid%exist_f2d, "Error: grid%f2d is still allocated")
    call assert(.not.grid%exist_f3d, "Error: grid%f3d is still allocated")
    call assert(.not.grid%exist_f4d, "Error: grid%f4d is still allocated")
    if(present(w)) optargs=optargs+1
    if(present(x)) optargs=optargs+1
    if(present(y)) optargs=optargs+1
    if(present(z)) optargs=optargs+1
    if(optargs.eq.4) then
       call assert(size(f,1)==size(w), "Error: size(fx,1)" //size(f,1)// ' /= =size(w) '//size(w))
       call assert(size(f,2)==size(x), "Error: size(fx,2)" //size(f,2)// ' /= =size(x) '//size(x))
       call assert(size(f,3)==size(y), "Error: size(fx,3)" //size(f,3)// ' /= =size(y) '//size(y))
       call assert(size(f,4)==size(z), "Error: size(fx,4)" //size(f,4)// ' /= =size(z) '//size(z))
       call assert(sorted(w(:)),"not sorted w")
       call assert(sorted(x(:)),"not sorted x")
       call assert(sorted(y(:)),"not sorted y")
       call assert(sorted(z(:)),"not sorted z")
       call assert(size(w(:))>=2,"not sorted w")
       call assert(size(x(:))>=2,"not sorted x")
       call assert(size(y(:))>=2,"not sorted y")
       call assert(size(z(:))>=2,"not sorted z")
       allocate(grid%axe(4))
       allocate(grid%axe(1)%x(size(w)))
       allocate(grid%axe(2)%x(size(x)))
       allocate(grid%axe(3)%x(size(y)))
       allocate(grid%axe(4)%x(size(z)))
       grid%axe(1)%x(:) = w(:)
       grid%axe(2)%x(:) = x(:)
       grid%axe(3)%x(:) = y(:)
       grid%axe(4)%x(:) = z(:)
       grid%axe(1)%dim_x = size(w)
       grid%axe(2)%dim_x = size(x)
       grid%axe(3)%dim_x = size(y)
       grid%axe(4)%dim_x = size(z)
       grid%axe(1)%allocated = .true.
       grid%axe(2)%allocated = .true.
       grid%axe(3)%allocated = .true.
       grid%axe(4)%allocated = .true.
    else
       call assert(optargs==0, "Error: wrong number of optionl arguments for new_grid_4d")
    endif
    allocate(grid%f4d(size(f,1),size(f,2),size(f,3),size(f,4)))
    grid%f4d(:,:,:,:)=f(:,:,:,:)
    grid%exist_f4d=.true.
  end subroutine new_grid_4d


  !-------------------------------------------
  !------- delete ----------------------------
  !-------------------------------------------
  !this routine can also been used when the grid is not set.
  subroutine delete(grid)
    implicit none
    type(grid_t),intent(inout):: grid

    if(allocated(grid%f1d))then
       deallocate(grid%f1d)
       grid%exist_f1d=.false.
    endif
    if(allocated(grid%f2d))then
       deallocate(grid%f2d)
       grid%exist_f2d=.false.
    endif
    if(allocated(grid%f3d))then
       deallocate(grid%f3d)
       grid%exist_f3d=.false.
    endif
    if(allocated(grid%f4d))then
       deallocate(grid%f4d)
       grid%exist_f4d=.false.
    endif
    if(allocated(grid%axe)) then
       deallocate(grid%axe)
    endif
  end subroutine delete

  !---------------------------------
  !----- set_option ----------------
  !---------------------------------
  subroutine set_option(grid,warning)
    implicit none
    type(grid_t)         :: grid
    logical,optional     :: warning

    if(present(warning)) grid%with_warning =warning

  end subroutine set_option


  subroutine interpol(w,x,y,z,fd1,fd2,fd3,fd4,grid,data_error)
    use amns_external_functions
    implicit none
    real (rkind), target           ,intent(in) :: w(:)
    real (rkind), target, optional ,intent(in) :: x(:)
    real (rkind), target, optional ,intent(in) :: y(:)
    real (rkind), target, optional ,intent(in) :: z(:)
    real (rkind),         optional ,intent(out):: fd1(:)
    real (rkind),         optional ,intent(out):: fd2(:,:)
    real (rkind),         optional ,intent(out):: fd3(:,:,:)
    real (rkind),         optional ,intent(out):: fd4(:,:,:,:)
    type (grid_t)                ,intent(inout):: grid
    type (data_error_t),            intent(out):: data_error

    ! interpolation index
    !
    real (rkind) :: dw(size(w))
    integer      :: iw(size(w))
    integer      :: iw_sort(size(w))
    real (rkind), target :: w_s(size(w))
    real (rkind), pointer :: w_p(:)

    real (rkind), allocatable         :: dx(:)  !size(x))
    integer , allocatable             :: ix(:)  !size(x))
    integer , allocatable             :: ix_sort(:)  !size(x))
    real (rkind), allocatable, target :: x_s(:)  !size(x))
    real (rkind), pointer             :: x_p(:)

    real (rkind), allocatable         :: dy(:)  !size(yvec))
    integer , allocatable             :: iy(:)  !size(yvec))
    integer, allocatable              :: iy_sort(:)  !size(yvec))
    real (rkind), allocatable, target :: y_s(:)  !size(yvec))
    real (rkind), pointer             :: y_p(:)

    real (rkind), allocatable         :: dz(:)  !size(zvec))
    integer  , allocatable            :: iz(:)  !size(zvec))
    integer , allocatable             :: iz_sort(:)  !size(zvec))
    real (rkind), allocatable, target :: z_s(:)  !size(zvec))
    real (rkind), pointer             :: z_p(:)

    !local
    logical  :: with_x      ! =present(x)
    logical  :: with_y      ! =present(y)
    logical  :: with_z      ! =present(z)
    logical  :: one_d       ! =present(fd1)
    logical  :: two_d       ! =present(fd2)
    logical  :: three_d     ! =present(fd3)
    logical  :: four_d      ! =present(fd4)
    logical  :: is_sorted_w = .false. ! =present(fgrid3)
    logical  :: is_sorted_x = .false. ! =present(fgrid3)
    logical  :: is_sorted_y = .false. ! =present(fgrid3)
    logical  :: is_sorted_z = .false. ! =present(fgrid3)

    ! help swap
    integer :: iswap
    !loop index
    integer :: i, idpc

    type (fun_err_t) :: fun_err

    data_error%ierr = 0
    data_error%cerr = ''
    fun_err%ierr = 0
    fun_err%cerr = ''

    with_x   =present(x)
    with_y   =present(y)
    with_z   =present(z)
    one_d    =present(fd1)
    two_d    =present(fd2)
    three_d  =present(fd3)
    four_d   =present(fd4)

    !-irregular calls---
    call assert(grid%exist_f1d .or. grid%exist_f2d .or.grid%exist_f3d.or.grid%exist_f4d, &
         "Data not allocated for interpolation")

    if(one_d) then
       call assert(one_d.and. .not.(two_d.or.three_d.or.four_d), &
            "Type mismatch, only one f vector type possible!!")
       if(with_x) call assert(size(x)==size(w), "Size w:"//size(w)// &
            "/= Size x:"//size(x))
       if(with_y) call assert(size(y)==size(w), "Size w:"//size(w)// &
            "/= Size y:"//size(y))
       if(with_z) call assert(size(z)==size(w), "Size w:"//size(x)// &
            "/= Size z:"//size(z))
    endif

    if(two_d)then
       call assert(with_x.and..not.(with_y.or.with_z),"expect x present, y and z not present")
    endif
    if(three_d)then
       call assert(with_y .and. .not. with_z,"expect x and y present and z not present")
    endif
    if(four_d)then
       call assert(with_y .and. with_z,"expect x y and zvec present")
    endif

    select case(grid%interpol_function)

    case (0)
       if(with_x)then
          allocate(dx(size(x)),ix(size(x)),ix_sort(size(x)),x_s(size(x)))
       endif
       if(with_y)then
          allocate(dy(size(y)),iy(size(y)),iy_sort(size(y)),y_s(size(y)))
       endif
       if(with_z)then
          allocate(dz(size(z)),iz(size(z)),iz_sort(size(z)),z_s(size(z)))
       endif

       is_sorted_w=sorted(w(:))
       if(with_x)is_sorted_x=sorted(x(:))
       if(with_y)is_sorted_y=sorted(y(:))
       if(with_z)is_sorted_z=sorted(z(:))

    ! if the vectors are not sorted, we create a sorted vector, and save index map
    ! pointer x_p will show a sorted vector

       if(.not.is_sorted_w) then
          iw_sort(:)= (/ (i,i=1,size(iw_sort(:))   ) /)
          call mrgrnk(w(:),iw_sort(:))
          w_s(:)=w(iw_sort(:))
          w_p => w_s
          if(grid%axe(1)%is_log) w_p = log10(w_p)
       else
          if(grid%axe(1)%is_log) then
             w_s = log10(w)
             w_p => w_s
          else
             w_p => w
          endif
       endif

       if(with_x) then
          if(.not.is_sorted_x) then
             ix_sort(:)= (/ (i,i=1,size(ix_sort(:))   ) /)
             call mrgrnk(x(:),ix_sort(:))
             x_s(:)=x(ix_sort(:))
             x_p => x_s
             if(grid%axe(2)%is_log) x_p = log10(x_p)
          else
             if(grid%axe(2)%is_log) then
                x_s = log10(x)
                x_p => x_s
             else
                x_p => x
             endif
          endif
       endif

       if(with_y) then
          if(.not.is_sorted_y)then
             iy_sort(:)= (/ (i,i=1,size(iy_sort(:))   ) /)
             call mrgrnk(y(:),iy_sort(:))
             y_s(:)=y(iy_sort(:))
             y_p => y_s
             if(grid%axe(3)%is_log) y_p = log10(y_p)
          else
             if(grid%axe(3)%is_log) then
                y_s = log10(y)
                y_p => y_s
             else
                y_p => y
             endif
          endif
       endif

       if(with_z) then
          if(.not.is_sorted_z)then
             iz_sort(:)= (/ (i,i=1,size(iz_sort(:))   ) /)
             call mrgrnk(z(:),iz_sort(:))
             z_s(:)=z(ix_sort(:))
             z_p => z_s
             if(grid%axe(4)%is_log) z_p = log10(z_p)
          else
             z_p => z
             if(grid%axe(4)%is_log) then
                z_s = log10(z)
                z_p => z_s
             else
                z_p => z
             endif
          endif
       endif


    !---get index of w in the grid
       call interpol_index(w_p(:),grid%axe(1)%x(:),iw(:),dw(:),grid)
    ! reset ix to original xvec positions
    ! next comand work, because internaly the f90 compiler works with copy comands

       if(.not. is_sorted_w)then
          iw(iw_sort(:))=iw(:)
          dw(iw_sort(:))=dw(:)
       endif


       if(with_x) then
          idpc=size(x)/2
          call interpol_index(x_p(:),grid%axe(2)%x(:),ix(:),dx(:),grid)
          if(.not. is_sorted_x)then
             ix(ix_sort(:))=ix(:)
             dx(ix_sort(:))=dx(:)
          endif
       endif

       if(with_y) then
          call interpol_index(y_p(:),grid%axe(3)%x(:),iy(:),dy(:),grid)
          if(.not. is_sorted_y)then
             iy(iy_sort(:))=iy(:)
             dy(iy_sort(:))=dy(:)
          endif
       endif

       if(with_z) then
          call interpol_index(z_p(:),grid%axe(4)%x(:),iz(:),dz(:),grid)
          if(.not. is_sorted_z)then
             iz(iz_sort(:))=iz(:)
             dz(iz_sort(:))=dz(:)
          endif
       endif
    !
    !-- now interpolate

       if(with_z)then
          if(one_d)then
             call interpo_4_1(iw(:),ix(:),iy(:),iz(:),dw(:),dx(:),dy(:),dz(:),fd1(:),grid)
          else
             call interpo_4_4(iw(:),ix(:),iy(:),iz(:),dw(:),dx(:),dy(:),dz(:),fd4(:,:,:,:),grid)
          endif
       elseif(with_y)then
          if(one_d)then
             call interpo_3_1(iw(:),ix(:),iy(:),dw(:),dx(:),dy(:),fd1(:),grid)
          else
             call interpo_3_3(iw(:),ix(:),iy(:),dw(:),dx(:),dy(:),fd3(:,:,:),grid)
          endif
       elseif(with_x)then
          if(one_d)then
             call interpo_2_1(iw(:),ix(:),dw(:),dx(:),fd1(:),grid)
          else
             call interpo_2_2(iw(:),ix(:),dw(:),dx(:),fd2(:,:),grid)
          endif
       else
          call interpo_1_1(iw(:),dw(:),fd1(:),grid)
       endif

       if(with_x)then
          deallocate(dx,ix,ix_sort,x_s)
       endif
       if(with_y)then
          deallocate(dy,iy,iy_sort,y_s)
       endif
       if(with_z)then
          deallocate(dz,iz,iz_sort,z_s)
       endif

       if(grid%is_log) then
          if(one_d)   fd1 = 10.0_rkind ** fd1
          if(two_d)   fd2 = 10.0_rkind ** fd2
          if(three_d) fd3 = 10.0_rkind ** fd3
          if(four_d)  fd4 = 10.0_rkind ** fd4
       endif

    case (1001)    ! grid%interpol_function
       call assert(.not.(with_x.or.with_y.or.with_z), "nuclear data '1001' must be a function of 1 parameter")
       call assert(grid%exist_f2d, "nuclear data '1001' requires 2d table of data")
       call assert(one_d, "nuclear data '1001' expected to be passed an effective 1d array")
       call nuclear_data_1001(grid%f2d, w, fd1, grid%with_warning, fun_err)

    case (1002)    ! grid%interpol_function
       call assert(.not.(with_x.or.with_y.or.with_z), "nuclear data '1002' must be a function of 1 parameter")
       call assert(grid%exist_f2d, "nuclear data '1002' requires 2d table of data")
       call assert(one_d, "nuclear data '1002' expected to be passed an effective 1d array")
       call nuclear_data_1002(grid%f2d, w, fd1, grid%with_warning, fun_err)

    case (1003)    ! grid%interpol_function
       call assert(.not.(with_x.or.with_y.or.with_z), "resonant charge data '1003' must be a function of 1 parameter")
       call assert(grid%exist_f1d, "resonant charge transfer data '1003' requires 1d table of data")
       call assert(one_d, "resonant charge transfer '1003' expected to be passed an effective 1d array")
       call rct_data_1003(grid%f1d, w, fd1, grid%with_warning, fun_err)

    case (1004)    ! grid%interpol_function
       call assert(.not.(with_y.or.with_z), "sputter data '1004' must be a function of 2 parameters")
       call assert(grid%exist_f2d, "sputter yield '1004' requires 2d table of data")
       call assert(one_d, "sputter yield '1004' expected to be passed an 2d array")
       call sputter_data_1004(grid%f2d, w, x, fd1, grid%with_warning, fun_err)

    case (1005)    ! grid%interpol_function
       call assert(.not.(with_y.or.with_z), "reflection yield '1005' must be a function of 2 parameters")
       call assert(grid%exist_f2d, "reflection yield '1005' requires 2d table of data")
       call assert(one_d, "reflection yield '1005' expected to be passed an 2d array")
       call reflect_data_1005(grid%f2d, w, x, fd1, grid%with_warning, fun_err)

    case (1006)    ! grid%interpol_function
       call assert(.not.(with_x.or.with_y.or.with_z), "nuclear data '1006' must be a function of 1 parameter")
       call assert(grid%exist_f2d, "nuclear data '1006' requires 2d table of data")
       call assert(one_d, "nuclear data '1006' expected to be passed an effective 1d array")
       call nuclear_data_1006(grid%f2d, w, fd1, grid%with_warning, fun_err)

    case default
       write(*,*) 'Case for grid%interpol_function = ', grid%interpol_function, ' not yet coded'

    end select

    if(fun_err%ierr.ne.0) then
       data_error%ierr = fun_err%ierr
       data_error%cerr = fun_err%cerr
    endif

  end subroutine interpol

  !--------------------------------------
  !--- sorted ---------------------------
  !--------------------------------------
  ! assending order
  function sorted(x) result(result)
    real(rkind), intent (in) :: x(:)
    logical ::result

    !  logical:: order
    integer:: i

    result=.true.
    do i=1,size(x)-1
       if(x(i)>x(i+1))then
          result=.false.
          exit
       endif
    enddo
  end function sorted

  !------------------------------------------
  !------ interpol_index --------------------
  !------------------------------------------
  subroutine interpol_index(x, axe, ix, dx, grid)
    implicit none
    real(rkind)   , intent(in) :: x(:)
    real(rkind)   , intent(out):: axe(:)
    real(rkind)   , intent(out):: dx(:)
    integer(ikind), intent(out):: ix(:)
    type(grid_t)  , intent(in) :: grid


    real(rkind)   :: dx_help(size(dx))
    integer(ikind):: ix_lb,ix_ub

    ! after this call
    !  ix_lb points to the 1st point >= xaxe(1)
    !  ux_ub points to the last point <= xaxe(size(axe))
    !  for these outside points ix is set to 0 or size(axe)+1 and
    !    dx has been set properly
    call bound_check(axe(:), x(:),ix(:),dx(:),ix_lb,ix_ub, grid)

    if(ix_lb>ix_ub) return
    ! find other location of x in respect of bounds
    call locate(axe, x(ix_lb),ix(ix_lb),dx(ix_lb))
    call hunt(axe, x(ix_lb+1:ix_ub),ix(ix_lb+1:ix_ub),dx(ix_lb+1:ix_ub))
    ! else
    !   call locate(axe, xvec(:),ix(:),dx(:))
    ! endif

  end subroutine interpol_index

  !---------------------------------------------
  !--- bound_check ------------------------------
  !---------------------------------------------
  subroutine bound_check(axe, x, ix, dx, lb, ub, grid)
    implicit none
    real(rkind)   , intent(in) :: axe(:)
    real(rkind)   , intent(in) :: x(:)
    real(rkind)   , intent(out):: dx(:)
    integer(ikind), intent(out):: ix(:)
    integer(ikind), intent(out):: lb, ub
    type(grid_t)  , intent(in) :: grid


    integer:: i

    !axe and x are sorted in ascending order
    !left bound
    lb=1
    do i=1, size(x)
       if(x(i)>=axe(1)) then
          exit
       else
          ix(i)=0
          lb=lb+1
          dx(i)=(x(i)-axe(1))/(axe(2)-axe(1))
          if(grid%with_warning) call warning(.false.,"Warning "// x(i) &
               // " < "// axe(1)// 'in grid: ?')
       endif
    enddo

    ! right bound
    ub=size(x)
    do i=size(x),1,-1
       if(x(i)<=axe(size(axe))) then
          exit
       else
          ix(i)=size(axe)+1
          ub=ub-1
          dx(i)=(x(i)-axe(size(axe)-1))/(axe(size(axe))-axe(size(axe)-1))
          if(grid%with_warning) call warning(.false.," Point "// x(i) &
               // " out of bound in grid ")
       endif
    enddo

  end subroutine bound_check

  !------------------------------------------
  !------- interpo_1_1 ----------------------
  !------------------------------------------
  subroutine interpo_1_1(ix,dx,fgrid1,grid)
    implicit none
    integer(ikind), intent(in) :: ix(:)
    real(rkind), intent(out):: fgrid1(:)
    real(rkind), intent(in) :: dx(:)
    type(grid_t),intent(in) :: grid

    integer :: i
    integer :: ii1, ii2
    !lin
    do i=1,size(ix(:))
       ii1=ix(i)
       if(ix(i)<lbound(grid%f1d,1)) ii1=lbound(grid%f1d,1)
       if(ix(i)>ubound(grid%f1d,1)) ii1=ubound(grid%f1d,1)-1
       ii2=ii1+1

       fgrid1(i)=grid%f1d(ii1) &
            +(grid%f1d(ii2) &
            -grid%f1d(ii1))*dx(i)
    enddo
  end subroutine interpo_1_1


  !------------------------------------------
  !------- interpo_2_1 ----------------------
  !------------------------------------------
  subroutine interpo_2_1(ix,iy,dx,dy,fgrid1,grid)
    implicit none
    integer(ikind), intent(in) :: ix(:)
    integer(ikind), intent(in) :: iy(:)
    real(rkind), intent(out):: fgrid1(:)
    real(rkind), intent(in) :: dx(:)
    real(rkind), intent(in) :: dy(:)
    type(grid_t),intent(in) :: grid

    real(rkind)::f1, f2, f3
    integer :: i
    integer :: ii1, ii2, jj1, jj2
    !lin
    do i=1,size(ix(:))
       ii1=ix(i)
       if(ii1<lbound(grid%f2d,1)) ii1=lbound(grid%f2d,1)
       if(ii1>=ubound(grid%f2d,1)) ii1=ubound(grid%f2d,1)-1
       ii2=ii1+1
       jj1=iy(i)
       if(jj1<lbound(grid%f2d,2)) jj1=lbound(grid%f2d,2)
       if(jj1>=ubound(grid%f2d,2)) jj1=ubound(grid%f2d,2)-1
       jj2=jj1+1
!!!
       if(ii1.lt.1.or.ii2.gt.ubound(grid%f2d,1)) then
          write(*,*) 'Error: ii1, ii2 = ', ii1, ii2
       endif
       if(jj1.lt.1.or.jj2.gt.ubound(grid%f2d,2)) then
          write(*,*) 'Error: jj1, jj2 = ', jj1, jj2
       endif
       f1=grid%f2d(ii1,jj1) &
            +(grid%f2d(ii2,jj1) &
            -grid%f2d(ii1,jj1))*dx(i)
       f2=grid%f2d(ii1,jj2) &
            +(grid%f2d(ii2,jj2) &
            -grid%f2d(ii1,jj2))*dx(i)
       fgrid1(i)=f1+(f2-f1)*dy(i)
    enddo
  end subroutine interpo_2_1

  !------------------------------------------
  !------- interpo_2_2 ----------------------
  !------------------------------------------
  subroutine interpo_2_2(ix,iy,dx,dy,fgrid2,grid)
    implicit none
    integer(ikind), intent(in) :: ix(:)
    integer(ikind), intent(in) :: iy(:)
    real(rkind), intent(out):: fgrid2(:,:)
    real(rkind), intent(in) :: dx(:)
    real(rkind), intent(in) :: dy(:)
    type(grid_t),intent(in) :: grid

    real(rkind)::f1, f2
    integer :: i, j
    integer :: ii1, ii2, jj1, jj2
    !lin
    do j=1,size(iy(:))
       jj1=iy(j)
       if(jj1<lbound(grid%f2d,2)) jj1=lbound(grid%f2d,2)
       if(jj1>ubound(grid%f2d,2)) jj1=ubound(grid%f2d,2)-1
       jj2=jj1+1
       do i=1,size(ix(:))
          ii1=ix(i)
          if(ii1<lbound(grid%f2d,1)) ii1=lbound(grid%f2d,1)
          if(ii1>ubound(grid%f2d,1)) ii1=ubound(grid%f2d,1)-1
          ii2=ii1+1
          f1=grid%f2d(ii1,jj1) &
               +(grid%f2d(ii2,jj1) &
               -grid%f2d(ii1,jj1))*dx(i)
          f2=grid%f2d(ii1,jj2) &
               +(grid%f2d(ii2,jj2) &
               -grid%f2d(ii1,jj2))*dx(i)
          fgrid2(i,j)=f1+(f2-f1)*dy(j)
       enddo
    enddo
  end subroutine interpo_2_2

  !------------------------------------------
  !------- interpo_3_1 ----------------------
  !------------------------------------------
  subroutine interpo_3_1(ix,iy,iz,dx,dy,dz,fgrid1,grid)
    implicit none
    integer(ikind), intent(in) :: ix(:)
    integer(ikind), intent(in) :: iy(:)
    integer(ikind), intent(in) :: iz(:)
    real(rkind), intent(out):: fgrid1(:)
    real(rkind), intent(in) :: dx(:)
    real(rkind), intent(in) :: dy(:)
    real(rkind), intent(in) :: dz(:)
    type(grid_t),intent(in) :: grid

    real(rkind)::f1, f2, f3, f4
    integer :: i
    integer :: ii1, ii2, jj1, jj2, kk1, kk2
    !lin
    do i=1,size(ix(:))
       ii1=ix(i)
       if(ii1<lbound(grid%f3d,1)) ii1=lbound(grid%f3d,1)
       if(ii1>ubound(grid%f3d,1)) ii1=ubound(grid%f3d,1)-1
       ii2=ii1+1
       jj1=iy(i)
       if(jj1<lbound(grid%f3d,2)) jj1=lbound(grid%f3d,2)
       if(jj1>ubound(grid%f3d,2)) jj1=ubound(grid%f3d,2)-1
       jj2=jj1+1
       kk1=iz(i)
       if(kk1<lbound(grid%f3d,3)) kk1=lbound(grid%f3d,3)
       if(kk1>ubound(grid%f3d,3)) kk1=ubound(grid%f3d,3)-1
       kk2=kk1+1
       f1=grid%f3d(ii1,jj1,kk1) &
            +(grid%f3d(ii2,jj1,kk1) &
            -grid%f3d(ii1,jj1,kk1))*dx(i)
       f2=grid%f3d(ii1,jj2,kk1) &
            +(grid%f3d(ii2,jj2,kk1) &
            -grid%f3d(ii1,jj2,kk1))*dx(i)
       f3=f1+(f2-f1)*dy(i)
       f1=grid%f3d(ii1,jj1,kk2) &
            +(grid%f3d(ii2,jj1,kk2) &
            -grid%f3d(ii1,jj1,kk2))*dx(i)
       f2=grid%f3d(ii1,jj2,kk2) &
            +(grid%f3d(ii2,jj2,kk2) &
            -grid%f3d(ii1,jj2,kk2))*dx(i)
       f4=f1+(f2-f1)*dy(i)
       fgrid1(i)=f3+(f4-f3)*dz(i)
    enddo
  end subroutine interpo_3_1

  !------------------------------------------
  !------- interpo_3_3 ----------------------
  !------------------------------------------
  subroutine interpo_3_3(ix,iy,iz,dx,dy,dz,fgrid3,grid)
    implicit none
    integer(ikind), intent(in) :: ix(:)
    integer(ikind), intent(in) :: iy(:)
    integer(ikind), intent(in) :: iz(:)
    real(rkind), intent(out):: fgrid3(:,:,:)
    real(rkind), intent(in) :: dx(:)
    real(rkind), intent(in) :: dy(:)
    real(rkind), intent(in) :: dz(:)
    type(grid_t),intent(in) :: grid

    real(rkind)::f1, f2, f3, f4
    integer :: i,j,k
    integer :: ii1, ii2, jj1, jj2, kk1, kk2
    !lin
    do k=1,size(iz(:))
       kk1=iz(k)
       if(kk1<lbound(grid%f3d,3)) kk1=lbound(grid%f3d,3)
       if(kk1>ubound(grid%f3d,3)) kk1=ubound(grid%f3d,3)-1
       kk2=kk1+1
       do j=1,size(iy(:))
          jj1=iy(j)
          if(jj1<lbound(grid%f3d,2)) jj1=lbound(grid%f3d,2)
          if(jj1>ubound(grid%f3d,2)) jj1=ubound(grid%f3d,2)-1
          jj2=jj1+1
          do i=1,size(ix(:))
             ii1=ix(i)
             if(ii1<lbound(grid%f3d,1)) ii1=lbound(grid%f3d,1)
             if(ii1>ubound(grid%f3d,1)) ii1=ubound(grid%f3d,1)-1
             ii2=ii1+1
             f1=grid%f3d(ii1,jj1,kk1)&
                  +(grid%f3d(ii2,jj1,kk1) &
                  -grid%f3d(ii1,jj1,kk1))*dx(i)
             f2=grid%f3d(ii1,jj2,kk1) &
                  +(grid%f3d(ii2,jj2,kk1) &
                  -grid%f3d(ii1,jj2,kk1))*dx(i)
             f3=f1+(f2-f1)*dy(j)
             f1=grid%f3d(ii1,jj1,kk2) &
                  +(grid%f3d(ii2,jj1,kk2) &
                  -grid%f3d(ii1,jj1,kk2))*dx(i)
             f2=grid%f3d(ii1,jj2,kk2) &
                  +(grid%f3d(ii2,jj2,kk2) &
                  -grid%f3d(ii1,jj2,kk2))*dx(i)
             f4=f1+(f2-f1)*dy(j)
             fgrid3(i,j,k)=f3+(f4-f3)*dz(k)
          enddo
       enddo
    enddo
  end subroutine interpo_3_3

  !------------------------------------------
  !------- interpo_4_1 ----------------------
  !------------------------------------------
  subroutine interpo_4_1(iw,ix,iy,iz,dw,dx,dy,dz,fgrid1,grid)
    implicit none
    integer(ikind), intent(in) :: iw(:)
    integer(ikind), intent(in) :: ix(:)
    integer(ikind), intent(in) :: iy(:)
    integer(ikind), intent(in) :: iz(:)
    real(rkind), intent(out):: fgrid1(:)
    real(rkind), intent(in) :: dw(:)
    real(rkind), intent(in) :: dx(:)
    real(rkind), intent(in) :: dy(:)
    real(rkind), intent(in) :: dz(:)
    type(grid_t),intent(in) :: grid

    real(rkind)::f1, f2, f3, f4, f5, f6
    integer :: i
    integer :: ll1,ll2,ii1,ii2,jj1,jj2,kk1,kk2
    !lin
    do i=1,size(iw(:))
       ll1=iw(i)
       if(ll1<lbound(grid%f4d,1)) ll1=lbound(grid%f4d,1)
       if(ll1>ubound(grid%f4d,1)) ll1=ubound(grid%f4d,1)-1
       ll2=ll1+1
       ii1=ix(i)
       if(ii1<lbound(grid%f4d,2)) ii1=lbound(grid%f4d,2)
       if(ii1>ubound(grid%f4d,2)) ii1=ubound(grid%f4d,2)-1
       ii2=ii1+1
       jj1=iy(i)
       if(jj1<lbound(grid%f4d,3)) jj1=lbound(grid%f4d,3)
       if(jj1>ubound(grid%f4d,3)) jj1=ubound(grid%f4d,3)-1
       jj2=jj1+1
       kk1=iz(i)
       if(kk1<lbound(grid%f4d,4)) kk1=lbound(grid%f4d,4)
       if(kk1>ubound(grid%f4d,4)) kk1=ubound(grid%f4d,4)-1
       kk2=kk1+1
       ! building interpolation for first cube
       f1=grid%f4d(ll1,ii1,jj1,kk1) &
            +(grid%f4d(ll2,ii1,jj1,kk1) &
            -grid%f4d(ll1,ii1,jj1,kk1))*dw(i)
       f2=grid%f4d(ll1,ii2,jj1,kk1) &
            +(grid%f4d(ll2,ii2,jj1,kk1) &
            -grid%f4d(ll1,ii2,jj1,kk1))*dw(i)
       f3=f1+(f2-f1)*dx(i)
       f1=grid%f4d(ll1,ii1,jj2,kk1) &
            +(grid%f4d(ll2,ii1,jj2,kk1) &
            -grid%f4d(ll1,ii1,jj2,kk1))*dw(i)
       f2=grid%f4d(ll1,ii2,jj2,kk1) &
            +(grid%f4d(ll2,ii2,jj2,kk1) &
            -grid%f4d(ll1,ii2,jj2,kk1))*dw(i)
       f4=f1+(f2-f1)*dx(i)

       f5=f3+(f4-f3)*dy(i)
       ! done
       ! build second cube
       f1=grid%f4d(ll1,ii1,jj1,kk2) &
            +(grid%f4d(ll2,ii1,jj1,kk2) &
            -grid%f4d(ll1,ii1,jj1,kk2))*dw(i)
       f2=grid%f4d(ll1,ii2,jj1,kk2) &
            +(grid%f4d(ll2,ii2,jj1,kk2) &
            -grid%f4d(ll1,ii2,jj1,kk2))*dw(i)
       f3=f1+(f2-f1)*dx(i)
       f1=grid%f4d(ll1,ii1,jj2,kk2) &
            +(grid%f4d(ll2,ii1,jj2,kk2) &
            -grid%f4d(ll1,ii1,jj2,kk2))*dw(i)
       f2=grid%f4d(ll1,ii2,jj2,kk2) &
            +(grid%f4d(ll2,ii2,jj2,kk2) &
            -grid%f4d(ll1,ii2,jj2,kk2))*dw(i)
       f4=f1+(f2-f1)*dx(i)

       f6=f3+(f4-f3)*dy(i)
       ! done
       fgrid1(i)=f5+(f6-f5)*dz(i)

    enddo

  end subroutine interpo_4_1


  !------------------------------------------
  !------- interpo_4_4 ----------------------
  !------------------------------------------
  subroutine interpo_4_4(iw,ix,iy,iz,dw,dx,dy,dz,fgrid1,grid)
    implicit none
    integer(ikind), intent(in) :: iw(:)
    integer(ikind), intent(in) :: ix(:)
    integer(ikind), intent(in) :: iy(:)
    integer(ikind), intent(in) :: iz(:)
    real(rkind), intent(out):: fgrid1(:,:,:,:)
    real(rkind), intent(in) :: dw(:)
    real(rkind), intent(in) :: dx(:)
    real(rkind), intent(in) :: dy(:)
    real(rkind), intent(in) :: dz(:)
    type(grid_t),intent(in) :: grid

    real(rkind)::f1, f2, f3, f4, f5, f6
    integer :: l,i,j,k
    integer :: ll1,ll2,ii1,ii2,jj1,jj2,kk1,kk2
    !lin
    do k=1,size(iz(:))
       kk1=iz(i)
       if(kk1<lbound(grid%f4d,4)) kk1=lbound(grid%f4d,4)
       if(kk1>ubound(grid%f4d,4)) kk1=ubound(grid%f4d,4)-1
       kk2=kk1+1
       do j=1,size(iy(:))
          jj1=iy(i)
          if(jj1<lbound(grid%f4d,3)) jj1=lbound(grid%f4d,3)
          if(jj1>ubound(grid%f4d,3)) jj1=ubound(grid%f4d,3)-1
          jj2=jj1+1
          do i=1,size(ix(:))
             ii1=ix(i)
             if(ii1<lbound(grid%f4d,2)) ii1=lbound(grid%f4d,2)
             if(ii1>ubound(grid%f4d,2)) ii1=ubound(grid%f4d,2)-1
             ii2=ii1+1
             do l=1,size(iw(:))
                ll1=iw(i)
                if(ll1<lbound(grid%f4d,1)) ll1=lbound(grid%f4d,1)
                if(ll1>ubound(grid%f4d,1)) ll1=ubound(grid%f4d,1)-1
                ll2=ll1+1
  ! building interpolation for first cube
                f1=grid%f4d(ll1,ii1,jj1,kk1) &
                     +(grid%f4d(ll2,ii1,jj1,kk1) &
                     -grid%f4d(ll1,ii1,jj1,kk1))*dw(l)
                f2=grid%f4d(ll1,ii2,jj1,kk1) &
                     +(grid%f4d(ll2,ii2,jj1,kk1) &
                     -grid%f4d(ll1,ii2,jj1,kk1))*dw(l)
                f3=f1+(f2-f1)*dx(i)
                f1=grid%f4d(ll1,ii1,jj2,kk1) &
                     +(grid%f4d(ll2,ii1,jj2,kk1) &
                     -grid%f4d(ll1,ii1,jj2,kk1))*dw(l)
                f2=grid%f4d(ll1,ii2,jj2,kk1) &
                     +(grid%f4d(ll2,ii2,jj2,kk1) &
                     -grid%f4d(ll1,ii2,jj2,kk1))*dw(l)
                f4=f1+(f2-f1)*dx(i)

                f5=f3+(f4-f3)*dy(i)
  ! done
  ! build second cube
                f1=grid%f4d(ll1,ii1,jj1,kk2) &
                     +(grid%f4d(ll2,ii1,jj1,kk2) &
                     -grid%f4d(ll1,ii1,jj1,kk2))*dw(l)
                f2=grid%f4d(ll1,ii2,jj1,kk2) &
                     +(grid%f4d(ll2,ii2,jj1,kk2) &
                     -grid%f4d(ll1,ii2,jj1,kk2))*dw(l)
                f3=f1+(f2-f1)*dx(i)
                f1=grid%f4d(ll1,ii1,jj2,kk2) &
                     +(grid%f4d(ll2,ii1,jj2,kk2) &
                     -grid%f4d(ll1,ii1,jj2,kk2))*dw(l)
                f2=grid%f4d(ll1,ii2,jj2,kk2) &
                     +(grid%f4d(ll2,ii2,jj2,kk2) &
                     -grid%f4d(ll1,ii2,jj2,kk2))*dw(l)
                f4=f1+(f2-f1)*dx(i)

                f6=f3+(f4-f3)*dy(i)
  ! done
                fgrid1(l,i,j,k)=f5+(f6-f5)*dz(i)

             enddo
          enddo
       enddo
    enddo

  end subroutine interpo_4_4

  !------------------------------------------------
  !-------- hunt ----------------------------------
  !------------------------------------------------

  subroutine hunt(axe,x,ix,dx)
    implicit none
    ! be careful, the axe is not necessarily starting with low bound=1
    ! to use the hunt search from num.rec. to get index of x the vector
    ! of axe must be adjusted to start with idx 1. This is done by calling
    ! hunt with the internal x vec. of axe. The calling mechanisme will
    ! addmit then that the pass vector starts with one.
    !
    ! The return indexes are then reset to the correct index after calling hunt_nrv
    !
    ! Additionaly dx is computed for interpolation

    real(rkind)  ,  intent(inout) :: axe(:)
    real(rkind),    intent(in)    :: x(:)
    real(rkind),    intent(out)   :: dx(:)
    integer(ikind), intent(out)   :: ix(:)
    integer(ikind), save :: last_access=1

    integer :: lbnd, ubnd
    integer :: i

    lbnd=lbound(axe(:),1)
    ubnd=ubound(axe(:),1)


    ! call hunt_nrv with internal vector
    call hunt_nrv(axe(:), x(:), ix(:), last_access)
    ! addjust to bound of axe
    ix(:)=ix(:)+lbnd-1

    do i=1,size(x)
       if(ix(i)<ubnd)then
          dx(i)=(x(i)-axe(ix(i)))/(axe(ix(i)+1)-axe(ix(i)))
       else
          dx(i)=0
       endif
    enddo
  end subroutine hunt

  !------------------------------------------------
  !-------- hunt_nrv ------------------------------
  !------------------------------------------------
  ! it is guaranteed by calling that xx and x are indexed from 1 to size(..)
  !
  ! just like hunt in num recip
  ! but x is a vextor and a colecting index vector is pass
  !
  subroutine hunt_nrv(xx,x,ix,jlo)
    implicit none
    !..given an array xx of length n and a value x, this routine returns a value
    !..jlo such that x is between xx(jlo) and xx(jlo+1). the array xx must be
    !..monotonic. j=0 or j=n indicates that x is out of range. on input jlo is a
    !..guess for the table entry, and should be used as input for the next hunt.
    !..one assumes the next table entry is relatively close to the old one.

    real(rkind),    intent(in)    :: xx(:)
    real(rkind),    intent(in)    :: x(:)
    integer(ikind), intent(out)   :: ix(:)
    integer(ikind), intent(inout) :: jlo

    integer :: n,inc,jhi,jm
    logical :: ascnd

    integer :: i

    n=size(xx)

    !..logical function; true if ascending table, false if descending
    ascnd = (xx(n) >= xx(1))

    do i=1,size(x)
       !..initial guess is not useful; go to bisection immediatly
       if (jlo <= 0 .or. jlo > n) then
          jlo=0
          jhi=n+1
       else
          inc=1                       ! set the initial hunting increment

   !..this is the hunt up section
          if (x(i) >= xx(jlo) .eqv. ascnd) then
             do
                jhi=jlo+inc
                if (jhi > n) then             ! done hunting since end of table
                   jhi=n+1
                   exit
                else
                   if (x(i) < xx(jhi) .eqv. ascnd) exit
                   jlo=jhi                     !not done hunting
                   inc=inc+inc                 !so double increment
                end if
             end do                          ! and try again
          else                              ! hunt done
             jhi=jlo
             do
                jlo=jhi-inc
                if (jlo < 1) then             ! hunt done, since end of table
                   jlo=0
                   exit
                else
                   if (x(i) >= xx(jlo) .eqv. ascnd) exit
                   jhi=jlo                     !not done hunting
                   inc=inc+inc                 !so double increment
                end if
             end do                          ! and try again
          end if
       end if                              ! done hunting, value bracketed
       do                                  ! Hunting is done, begin final bisection phase
          if (jhi-jlo <= 1) then
             if (x(i) == xx(n)) jlo=n-1
             if (x(i) == xx(1)) jlo=1
             exit
          else
             jm=(jhi+jlo)/2
             if (x(i) >= xx(jm) .eqv. ascnd) then
                jlo=jm
             else
                jhi=jm
             end if
          end if
       end do
       ix(i)=jlo
    enddo!(i)
  end subroutine hunt_nrv

  !----------------------------------------------------------
  !------------ locate ------------------------------------
  !----------------------------------------------------------
  !
  ! same as hunt search, using midpoint method
  !

  subroutine locate(axe, x, ix, dx)
    implicit none
    !..given an array xx of length n, and a value of x, this routine returns
    !..a value j such that x is between xx(j) and xx(j+1). the array xx must be
    !..monotonic. j=0 or j=n indicates that x is out of range. bisection is used
    !..to find the entry

    real(rkind),    intent(inout) :: axe(:)
    real(rkind),    intent(in)    :: x
    real(rkind),    intent(out)   :: dx
    integer(ikind), intent(out)   :: ix

    integer  :: lbnd, ubnd
    integer  :: i

    lbnd=lbound(axe(:),1)
    ubnd=ubound(axe(:),1)

    ! get index
    ix=locate_nrv(axe(:),x)
    ! addjust to bound of axe
    ix=ix+lbnd-1

    if(ix<ubnd)then
       dx=(x-axe(ix))/(axe(ix+1)-axe(ix))
    else
       dx=0
    endif

    ! last search for hunt search
    ! axe%last_access=ix(size(x))+lbnd-1

  end subroutine locate


  !----------------------------------------------------------
  !------------ locate_nrv ----------------------------------
  !----------------------------------------------------------

  function locate_nrv(xx,x) result (result)
    !..given an array xx of length n, and a value of x, this routine returns
    !..a value j such that x is between xx(j) and xx(j+1). the array xx must be
    !..monotonic. j=0 or j=n indicates that x is out of range. bisection is used
    !..to find the entry

    real(rkind),  intent(in) :: xx(:)
    real(rkind),  intent(in) :: x

    integer      :: result
    integer      :: n,jl,jm,ju
    logical      :: ascnd

    integer  :: i

    n=size(xx)
    ascnd = (xx(n) >= xx(1))
    jl=0
    ju=n+1
    do
       if (ju-jl <= 1) exit
   !..compute a midpoint, and replace either the upper or lower limit
       jm=(ju+jl)/2
       if (ascnd .eqv. (x >= xx(jm))) then
          jl=jm
       else
          ju=jm
       end if
    end do
    if (x == xx(1)) then
       result=1
    else if (x == xx(n)) then
       result=n-1
    else
       result=jl
    end if

  end function locate_nrv

end module data_suport