utils.F90

Go to the documentation of this file.

! This file is part of the 'atomes' software.
!
! 'atomes' is free software: you can redistribute it and/or modify it under the terms
! of the GNU Affero General Public License as published by the Free Software Foundation,
! either version 3 of the License, or (at your option) any later version.
!
! 'atomes' is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
! See the GNU General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License along with 'atomes'.
! If not, see <https://www.gnu.org/licenses/>
!
! Copyright (C) 2022-2024 by CNRS and University of Strasbourg
!
 
!
! This file contains various tools
!
 
!********************************************************************
!
! Transformer un entier en chaine de charactères
! Transform an integer into character
!
 
SUBROUTINE charint(Word, Num)
 
IMPLICIT NONE
 
INTEGER :: lc
 
INTEGER, INTENT(IN) :: Num
 
CHARACTER (LEN=16) :: Pseudo
 
CHARACTER (LEN=15), INTENT(INOUT) :: Word
 
 
word=' '
write (word, *) num
lc=1
do while (word(1:lc) .eq. ' ')
 
  lc= lc+1
 
enddo
 
pseudo=word
word=' '
 
write (word, *) pseudo(lc:len_trim(pseudo))
 
END SUBROUTINE
 
SUBROUTINE creat_ring (RING_INIT, ELEM_CR)
 
USE parameters
 
IMPLICIT NONE
 
TYPE (RING), INTENT(INOUT) :: RING_INIT
INTEGER, INTENT(IN) :: ELEM_CR
 
ring_init%SPEC = 0
ring_init%ATOM = elem_cr
nullify(ring_init%PAST)
nullify(ring_init%NEXT)
 
END SUBROUTINE
 
SUBROUTINE do_ring (THE_RING, ELEM_DO)
 
USE parameters
 
IMPLICIT NONE
 
TYPE (RING), POINTER, INTENT(INOUT) :: THE_RING
!INTEGER, DIMENSION(NA), INTENT(INOUT) :: RPT
INTEGER, INTENT(IN) :: ELEM_DO
TYPE (RING), POINTER :: NEW
allocate(new, stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Subroutine: DO_RING"//char(0), "Pointer: NEW"//char(0))
  alc=.true.
  goto 001
endif
nullify(new%NEXT)
nullify(new%PAST)
the_ring%NEXT => new
new%PAST => the_ring
new%ATOM = elem_do
new%SPEC = the_ring%SPEC+1
the_ring => new
!if (RPT(THE_RING%ATOM) .eq. 0) then
!  RPT(THE_RING%ATOM) = 1
!endif
 
001 continue
 
END SUBROUTINE
 
!********************************************************************
!
! Calculer une distance
! Compute a distance including PBC
!
 
SUBROUTINE calcrij(AT1, AT2, STEP_1, STEP_2, SID)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: AT1, AT2, STEP_1, STEP_2, SID
INTEGER :: R1
DOUBLE PRECISION, DIMENSION(3) :: COORDA, COORDB
DOUBLE PRECISION, DIMENSION(3) :: Aij, Bij, Nij
 
if (step_1 .eq. -2) then
 
  do r1=1, 3
    coorda(r1) = nfullpos(at1,r1,step_2)
    coordb(r1) = nfullpos(at2,r1,step_2)
  enddo
 
elseif (step_1 .eq. -1) then
 
! Prepos for MSD
  do r1=1, 3
    coorda(r1) = pob(at2,r1)
    coordb(r1) = poa(at1,r1)
  enddo
 
else
! Standard inter-atomic distance
  do r1=1, 3
    coorda(r1) = fullpos(at1,r1,step_1)
    coordb(r1) = fullpos(at2,r1,step_2)
  enddo
 
endif
 
do r1=1, 3
  rij(r1) = coorda(r1) - coordb(r1)
enddo
 
if (.not.pbc) goto 001
 
! Test if the simulation box is cubic-like symmetric
! For that we check the lattice angles
if (the_box(sid)%GLASS) then
 
  do r1=1, 3
    rij(r1)=rij(r1) - anint(rij(r1)/(the_box(sid)%modv(r1)))*(the_box(sid)%modv(r1))
  enddo
 
else
 
  aij = matmul(coorda,the_box(sid)%carttofrac)
  bij = matmul(coordb,the_box(sid)%carttofrac)
  do r1=1,3
    nij(r1) = aij(r1) - bij(r1)
    nij(r1) = nij(r1) - anint(nij(r1))
  enddo
  rij = matmul(nij,the_box(sid)%fractocart)
 
endif
 
001 continue
 
dij=0.0d0
do r1=1, 3
  dij=dij+rij(r1)**2
enddo
 
END SUBROUTINE
 
DOUBLE PRECISION FUNCTION  calcdij(R12, AT1, AT2, STEP_1, STEP_2, SID)
 
USE parameters
 
IMPLICIT NONE
 
DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: r12
INTEGER, INTENT(IN) :: at1, at2, step_1, step_2, sid
INTEGER :: r1
DOUBLE PRECISION, DIMENSION(3) :: coorda, coordb
DOUBLE PRECISION, DIMENSION(3) :: aij, bij, nij
 
if (step_1 .eq. -2) then
 
  do r1=1, 3
    coorda(r1) = nfullpos(at1,r1,step_2)
    coordb(r1) = nfullpos(at2,r1,step_2)
  enddo
 
elseif (step_1 .eq. -1) then
 
! Prepos for MSD
  do r1=1, 3
    coorda(r1) = pob(at2,r1)
    coordb(r1) = poa(at1,r1)
  enddo
 
else
! Standard inter-atomic distance
  do r1=1, 3
    coorda(r1) = fullpos(at1,r1,step_1)
    coordb(r1) = fullpos(at2,r1,step_2)
  enddo
 
endif
 
do r1=1, 3
  r12(r1) = coorda(r1) - coordb(r1)
enddo
 
if (.not.pbc) goto 001
 
! Test if the simulation box is cubic-like symetric
! For that we check the lattice angles
if (the_box(sid)%GLASS) then
 
  do r1=1, 3
    r12(r1)=r12(r1) - anint(r12(r1)/(the_box(sid)%modv(r1)))*(the_box(sid)%modv(r1))
  enddo
 
else
 
  aij = matmul(coorda,the_box(sid)%carttofrac)
  bij = matmul(coordb,the_box(sid)%carttofrac)
  do r1=1,3
    nij(r1) = aij(r1) - bij(r1)
    nij(r1) = nij(r1) - anint(nij(r1))
  enddo
  r12 = matmul(nij,the_box(sid)%fractocart)
 
endif
 
001 continue
 
calcdij = 0.0d0
do r1=1, 3
  calcdij = calcdij + r12(r1)**2
enddo
 
END FUNCTION
 
!********************************************************************
!
! Home made arcosine calculation
!
 
DOUBLE PRECISION FUNCTION sacos (ANG)
 
USE parameters
 
DOUBLE PRECISION, INTENT(IN) :: ang
DOUBLE PRECISION :: scos
 
if (ang .lt. -1.0d0) then
  scos = -2.0d0 - ang
else if (ang .gt. 1.0d0) then
  scos = 2.0d0 - ang
else
  scos = ang
endif
sacos = dacos(scos)*180.d0/pi
 
END FUNCTION
 
!********************************************************************
!
! Calculer un angle direct
! Compute a direct angle
!
 
DOUBLE PRECISION FUNCTION angijk(ATG1, ATG2, ATG3, ASTEP)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: atg1, atg2, atg3, astep
INTEGER :: ag1, ag2
DOUBLE PRECISION :: daa, dbb, vang
DOUBLE PRECISION, DIMENSION(3) :: raa, rbb
 
INTERFACE
  DOUBLE PRECISION FUNCTION  calcdij(R12, AT1, AT2, STEP_1, STEP_2, SID)
    DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: r12
    INTEGER, INTENT(IN) :: at1, at2, step_1, step_2, sid
  END FUNCTION
  DOUBLE PRECISION FUNCTION sacos (ANG)
    DOUBLE PRECISION, INTENT(IN) :: ang
  END FUNCTION
END INTERFACE
 
if (ncells .gt. 1) then
  ag2 = astep
else
  ag2 = 1
endif
 
daa = calcdij(raa, atg2,atg1,astep,astep,ag2)
daa=sqrt(daa)
dbb = calcdij(rbb, atg2,atg3,astep,astep,ag2)
vang=0.0d0
do ag1=1, 3
  vang=vang+raa(ag1)*rbb(ag1)
enddo
dbb=sqrt(dbb)
 
angijk = sacos(vang/(daa*dbb))
 
END FUNCTION
 
!********************************************************************
!
! Calculer un angle dièdre
! Compute a dihedral angle
!
 
DOUBLE PRECISION FUNCTION diedre(DG1, DG2, DG3, DG4, DSTEP)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: dg1, dg2, dg3, dg4, dstep
INTEGER :: ag, ah
DOUBLE PRECISION :: dh1, dh2, dh3, vdi
DOUBLE PRECISION, DIMENSION(3) :: d12, d23, d34
DOUBLE PRECISION, DIMENSION(3) :: vh1, vh2
 
INTERFACE
  DOUBLE PRECISION FUNCTION  calcdij(R12, AT1, AT2, STEP_1, STEP_2, SID)
    DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: r12
    INTEGER, INTENT(IN) :: at1, at2, step_1, step_2, sid
  END FUNCTION
  DOUBLE PRECISION FUNCTION sacos (ANG)
    DOUBLE PRECISION, INTENT(IN) :: ang
  END FUNCTION
END INTERFACE
 
if (ncells .gt. 1) then
  ah = dstep
else
  ah = 1
endif
 
dh1 = calcdij(d12,dg1,dg2,dstep,dstep,ah)
dh2 = calcdij(d23,dg2,dg3,dstep,dstep,ah)
dh3 = calcdij(d34,dg3,dg4,dstep,dstep,ah)
 
vh1(1) = d12(2)*d23(3) - d12(3)*d23(2)
vh1(2) = d12(3)*d23(1) - d12(1)*d23(3)
vh1(3) = d12(1)*d23(2) - d12(2)*d23(1)
 
vh2(1) = d23(2)*d34(3) - d23(3)*d34(2)
vh2(2) = d23(3)*d34(1) - d23(1)*d34(3)
vh2(3) = d23(1)*d34(2) - d23(2)*d34(1)
 
dh1=0.0d0
dh2=0.0d0
do ag=1, 3
  dh1 = dh1 + vh1(ag)*vh1(ag)
  dh2 = dh2 + vh2(ag)*vh2(ag)
enddo
dh1=sqrt(dh1)
dh2=sqrt(dh2)
if (dh1.eq.0.0d0 .or. dh2.eq.0.0d0) then
  diedre = 0.0d0
else
  vdi = 0.0d0
  do ag=1, 3
    vdi=vdi+vh1(ag)*vh2(ag)
  enddo
  vdi = vdi/(dh1*dh2)
  diedre = sacos(vdi)
endif
 
END FUNCTION
 
! Cas particulier lors des statistiques d'anneaux
 
SUBROUTINE ect_type_rings(MOYENNE, TABLEAU, LONGTE, LREPRES, EC_TYPE)
 
INTEGER :: pose
INTEGER, INTENT(IN) :: LONGTE, LREPRES
 
DOUBLE PRECISION, INTENT(IN) :: MOYENNE
DOUBLE PRECISION, INTENT(IN), DIMENSION(LONGTE) :: TABLEAU
DOUBLE PRECISION, INTENT(INOUT) :: EC_TYPE
 
ec_type=0.0d0
 
do pose=1, longte
 
  if (tableau(pose) .ne. 0.0) ec_type= ec_type + (tableau(pose) - moyenne)**2
 
enddo
 
ec_type=sqrt(ec_type/dble(lrepres-1))
 
END SUBROUTINE
 
!********************************************************************
!
! Générateur de nombre aléatoire
! Numerical recipes random number generator 3
!
 
DOUBLE PRECISION FUNCTION ran3 (idnum)
 
INTEGER, INTENT(IN) :: idnum
INTEGER :: idum
INTEGER :: mbig,mseed,mz
INTEGER ::  iii,iff,ii,inext,inextp,k
INTEGER :: mj,mk,ma(55)
 
DOUBLE PRECISION :: fac
 
parameter(mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig)
SAVE  iff,inext,inextp,ma
DATA  iff /0/
 
idum = idnum
if(idum.lt.0 .or. iff.eq.0)then
 
  iff=1
  mj=mseed-iabs(idum)
  mj=mod(mj,mbig)
  ma(55)=mj
  mk=1
  do 11 iii=1,54
 
    ii=mod(21*iii,55)
    ma(ii)=mk
    mk=mj-mk
    if(mk.lt.mz)mk=mk+mbig
    mj=ma(ii)
 
11  continue
 
  do 13 k=1,4
 
    do 12 iii=1,55
 
      ma(iii)=ma(iii)-ma(1+mod(iii+30,55))
      if(ma(iii).lt.mz)ma(iii)=ma(iii)+mbig
 
12    continue
 
13  continue
  inext=0
  inextp=31
  idum=1
 
endif
 
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran3=mj*fac
 
END FUNCTION
 
REAL (KIND=c_double) FUNCTION random3(seed) bind (C,NAME='random3_')
 
use, INTRINSIC :: iso_c_binding
 
INTEGER (KIND=c_int), INTENT(IN) :: seed
 
INTERFACE
  DOUBLE PRECISION FUNCTION ran3 (idnum)
    INTEGER, INTENT(IN) :: idnum
  END FUNCTION
END INTERFACE
 
random3 = ran3(seed)
 
END FUNCTION
 
!********************************************************************
!
! Calculer la somme d'une série de valeurs
! Compute sum of a list of values
!
 
SUBROUTINE somme(TABS, LONGTS, SOMTAB)
 
INTEGER :: posm
INTEGER, INTENT(IN) :: LONGTS
 
DOUBLE PRECISION, INTENT(INOUT) :: SOMTAB
DOUBLE PRECISION, INTENT(IN), DIMENSION(LONGTS) :: TABS
 
somtab=0.0d0
 
do posm=1, longts
 
somtab=somtab+tabs(posm)
 
enddo
 
END SUBROUTINE
 
!********************************************************************
!
! Calculer la moyenne d'une série de valeurs
! Compute the average of list of values
!
 
SUBROUTINE moyenne(TABLEAU, LONGTM, MOYTAB)
 
INTEGER :: posm
INTEGER, INTENT(IN) :: LONGTM
 
DOUBLE PRECISION, INTENT(INOUT) :: MOYTAB
DOUBLE PRECISION, INTENT(IN), DIMENSION(LONGTM) :: TABLEAU
 
moytab=0.0d0
 
do posm=1, longtm
 
moytab=moytab+tableau(posm)
 
enddo
 
moytab=moytab/dble(longtm)
 
END SUBROUTINE
 
!********************************************************************
!
! Calculer l'écart type sur une série de valeurs
! Compute the standard error on a list of values
!
 
SUBROUTINE ect_type(MOYENNE, TABLEAU, LONGTE, EC_TYPE)
 
INTEGER :: pose
INTEGER, INTENT(IN) :: LONGTE
 
DOUBLE PRECISION, INTENT(IN) :: MOYENNE
DOUBLE PRECISION, INTENT(IN), DIMENSION(LONGTE) :: TABLEAU
DOUBLE PRECISION, INTENT(INOUT) :: EC_TYPE
 
ec_type=0.0d0
 
do pose=1, longte
 
ec_type= ec_type + (tableau(pose) - moyenne)**2
 
enddo
 
ec_type=sqrt(ec_type/dble(longte-1))
 
END SUBROUTINE
 
!********************************************************************
!
! Trier une liste d'entiers
! Sort a list of integers
! This is the most simple sort
! For an improved routine (huge data to sort) see SORT3 in this file
!
 
SUBROUTINE tri(TAB, DIMTAB)
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: DIMTAB
INTEGER, DIMENSION(DIMTAB), INTENT(INOUT) :: TAB
INTEGER :: SORTX, SORTY, VALUE
 
do sortx=2, dimtab
 
  VALUE=tab(sortx)
 
  do sorty=sortx-1, 1, -1
 
    if (tab(sorty) .le. VALUE) exit
 
    tab(sorty+1)=tab(sorty)
 
  enddo
 
  tab(sorty+1)=VALUE
 
enddo
 
END SUBROUTINE
 
!********************************************************************
!
! Lisser une courbe / To smooth a curve
!
 
LOGICAL FUNCTION smooth (TABTOLISS, GTOLISS, DIMTOLISS, SIGMALISS)
 
! Lissage selon une gaussienne
! Gaussian smoothing - zero based curve
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: dimtoliss
DOUBLE PRECISION, INTENT(IN) :: sigmaliss
DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMTOLISS) :: gtoliss
DOUBLE PRECISION, INTENT(INOUT), DIMENSION(DIMTOLISS) :: tabtoliss
 
INTEGER :: inda, indb, err
DOUBLE PRECISION, PARAMETER :: pi=acos(-1.0)
DOUBLE PRECISION :: dq, dq2
DOUBLE PRECISION :: factliss
 
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: newtab
 
allocate(newtab(dimtoliss), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: SMOOTH"//char(0), "Table: NEWTAB"//char(0))
  smooth=.false.
  goto 001
endif
 
do inda=1, dimtoliss
  newtab(inda)=0.0d0
enddo
 
factliss=1.0d0/(sigmaliss*sqrt(2.0d0*pi))
 
do inda=1, dimtoliss
  do indb=1, dimtoliss
 
    if(indb .eq. 1)then
      dq= gtoliss(2)-gtoliss(1)
    elseif (indb .eq. dimtoliss) then
      dq = gtoliss(dimtoliss)- gtoliss(dimtoliss-1)
    else
      dq = (gtoliss(indb+1)- gtoliss(indb-1))*0.5
    endif
 
    dq2=(gtoliss(inda)-gtoliss(indb))*(gtoliss(inda)-gtoliss(indb))
    newtab(inda) = newtab(inda) + exp(-dq2/(2.0d0*sigmaliss*sigmaliss))*tabtoliss(indb)*dq
 
  enddo
  newtab(inda)=factliss*newtab(inda)
enddo
 
do inda=1, dimtoliss
  tabtoliss(inda)=newtab(inda)
enddo
 
smooth=.true.
 
001 continue
 
if (allocated(newtab)) deallocate(newtab)
 
END FUNCTION
 
!********************************************************************
!
! Sort routine adapted from:
! Numerical sort array 3 - Quicksort INTEGER
!
 
SUBROUTINE sort3(arr, ndum)
 
INTEGER, INTENT(IN) :: ndum
INTEGER, DIMENSION(ndum) :: index
 
INTEGER, DIMENSION(ndum), INTENT(INOUT) :: arr
 
call indexx_i(arr, index, ndum)
 
arr=arr(index)
 
END SUBROUTINE
 
SUBROUTINE indexx_i(arr, index, n)
 
  IMPLICIT NONE
  INTEGER, PARAMETER :: NPAR_ARTH=16,npar2_arth=8
  INTEGER, DIMENSION(n), INTENT(INOUT) :: index
  INTEGER, PARAMETER :: NN=15, nstack=50
  INTEGER, INTENT(IN) :: n
  INTEGER :: k,i,j,indext,jstack,l,r
  INTEGER, DIMENSION(NSTACK) :: istack
  INTEGER :: a
  INTEGER, DIMENSION(n), INTENT(IN) :: arr
 
  call arth(1,1,n,index)
  jstack=0
  l=1
  r=n
  do
    if (r-l < nn) then
      do j=l+1,r
        indext=index(j)
        a=arr(indext)
        do i=j-1,l,-1
          if (arr(index(i)) <= a) exit
          index(i+1)=index(i)
        end do
        index(i+1)=indext
      end do
      if (jstack .eq. 0) RETURN
      r=istack(jstack)
      l=istack(jstack-1)
      jstack=jstack-2
    else
      k=(l+r)/2
      call swap(index(k),index(l+1))
      call icomp_xchg(index(l),index(r))
      call icomp_xchg(index(l+1),index(r))
      call icomp_xchg(index(l),index(l+1))
      i=l+1
      j=r
      indext=index(l+1)
      a=arr(indext)
      do
        do
          i=i+1
          if (arr(index(i)) >= a) exit
        end do
        do
          j=j-1
          if (arr(index(j)) <= a) exit
        end do
        if (j < i) exit
        call swap(index(i),index(j))
      end do
      index(l+1)=index(j)
      index(j)=indext
      jstack=jstack+2
      if (jstack > nstack) write (6, *) 'Sort3 error: indexx: NSTACK too small'
      if (r-i+1 >= j-l) then
        istack(jstack)=r
        istack(jstack-1)=i
        r=j-1
      else
        istack(jstack)=j-1
        istack(jstack-1)=l
        l=i
      end if
    end if
  end do
  CONTAINS
!BL
  SUBROUTINE icomp_xchg(i, j)
  INTEGER, INTENT(INOUT) :: i,j
  INTEGER :: swp
  if (arr(j) < arr(i)) then
    swp=i
    i=j
    j=swp
  end if
  END SUBROUTINE icomp_xchg
 
  SUBROUTINE arth(first, increment, m, arth_i)
  INTEGER, INTENT(IN) :: first,increment,m
  INTEGER, DIMENSION(m) :: arth_i
  INTEGER :: k,k2,temp
 
  if (m > 0) arth_i(1)=first
  if (m <= npar_arth) then
    do k=2,m
      arth_i(k)=arth_i(k-1)+increment
    end do
  else
    do k=2,npar2_arth
      arth_i(k)=arth_i(k-1)+increment
    end do
    temp=increment*npar2_arth
    k=npar2_arth
    do
      if (k >= m) exit
      k2=k+k
      arth_i(k+1:min(k2,m))=temp+arth_i(1:min(k,m-k))
      temp=temp+temp
      k=k2
    end do
  end if
  END SUBROUTINE
 
  SUBROUTINE swap(a,b)
  INTEGER, INTENT(INOUT) :: a,b
  INTEGER :: dum
  dum=a
  a=b
  b=dum
  END SUBROUTINE swap
 
END SUBROUTINE indexx_i
!********************************************************************
 
!********************************************************************
!
! Sort routine adapted from:
! Numerical sort array 3 - Quicksort DOUBLE PRECISION
!
 
SUBROUTINE sort3_dp (arr, ndum)
 
INTEGER, INTENT(IN) :: ndum
INTEGER, DIMENSION(ndum) :: index
 
DOUBLE PRECISION, DIMENSION(ndum), INTENT(INOUT) :: arr
 
call indexx_dp (arr, index, ndum)
 
arr=arr(index)
 
END SUBROUTINE
 
SUBROUTINE indexx_dp (arr, index, n)
 
  IMPLICIT NONE
  INTEGER, PARAMETER :: NPAR_ARTH=16,npar2_arth=8
  INTEGER, DIMENSION(n), INTENT(INOUT) :: index
  INTEGER, PARAMETER :: NN=15, nstack=50
  INTEGER, INTENT(IN) :: n
  INTEGER :: k,i,j,indext,jstack,l,r
  INTEGER, DIMENSION(NSTACK) :: istack
  DOUBLE PRECISION :: a
  DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: arr
 
  call arth(1,1,n,index)
  jstack=0
  l=1
  r=n
  do
    if (r-l < nn) then
      do j=l+1,r
        indext=index(j)
        a=arr(indext)
        do i=j-1,l,-1
          if (arr(index(i)) <= a) exit
          index(i+1)=index(i)
        end do
        index(i+1)=indext
      end do
      if (jstack .eq. 0) RETURN
      r=istack(jstack)
      l=istack(jstack-1)
      jstack=jstack-2
    else
      k=(l+r)/2
      call swap(index(k),index(l+1))
      call icomp_xchg(index(l),index(r))
      call icomp_xchg(index(l+1),index(r))
      call icomp_xchg(index(l),index(l+1))
      i=l+1
      j=r
      indext=index(l+1)
      a=arr(indext)
      do
        do
          i=i+1
          if (arr(index(i)) >= a) exit
        end do
        do
          j=j-1
          if (arr(index(j)) <= a) exit
        end do
        if (j < i) exit
        call swap(index(i),index(j))
      end do
      index(l+1)=index(j)
      index(j)=indext
      jstack=jstack+2
      if (jstack > nstack) write (6, *) 'Sort3 error: indexx: NSTACK too small'
      if (r-i+1 >= j-l) then
        istack(jstack)=r
        istack(jstack-1)=i
        r=j-1
      else
        istack(jstack)=j-1
        istack(jstack-1)=l
        l=i
      end if
    end if
  end do
  CONTAINS
!BL
  SUBROUTINE icomp_xchg(i, j)
  INTEGER, INTENT(INOUT) :: i,j
  INTEGER :: swp
  if (arr(j) < arr(i)) then
    swp=i
    i=j
    j=swp
  end if
  END SUBROUTINE icomp_xchg
 
  SUBROUTINE arth(first, increment, m, arth_i)
  INTEGER, INTENT(IN) :: first,increment,m
  INTEGER, DIMENSION(m) :: arth_i
  INTEGER :: k,k2,temp
 
  if (m > 0) arth_i(1)=first
  if (m <= npar_arth) then
    do k=2,m
      arth_i(k)=arth_i(k-1)+increment
    end do
  else
    do k=2,npar2_arth
      arth_i(k)=arth_i(k-1)+increment
    end do
    temp=increment*npar2_arth
    k=npar2_arth
    do
      if (k >= m) exit
      k2=k+k
      arth_i(k+1:min(k2,m))=temp+arth_i(1:min(k,m-k))
      temp=temp+temp
      k=k2
    end do
  end if
  END SUBROUTINE
 
  SUBROUTINE swap(a,b)
  INTEGER, INTENT(INOUT) :: a,b
  INTEGER :: dum
  dum=a
  a=b
  b=dum
  END SUBROUTINE swap
 
END SUBROUTINE indexx_dp
!********************************************************************