pw2casino: Source formatting (N. Nemec)

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6348 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
nn245 2010-02-02 16:39:38 +00:00
parent 7a35684787
commit 57ffc146b9
1 changed files with 221 additions and 221 deletions

View File

@ -42,7 +42,7 @@ SUBROUTINE write_casino_pwfn(gather)
LOGICAL :: exst, found
REAL(DP) :: ek, eloc, enl, charge, etotefield
COMPLEX(DP), ALLOCATABLE :: aux(:)
REAL(DP), allocatable :: v_h_new(:,:,:,:), kedtaur_new(:,:)
REAL(DP), ALLOCATABLE :: v_h_new(:,:,:,:), kedtaur_new(:,:)
INTEGER :: ios
INTEGER, EXTERNAL :: atomic_number
REAL (DP), EXTERNAL :: ewald, w1gauss
@ -52,16 +52,16 @@ SUBROUTINE write_casino_pwfn(gather)
REAL(DP), ALLOCATABLE :: g_l(:,:), g_g(:,:), g2(:)
COMPLEX(DP), ALLOCATABLE :: evc_l(:), evc_g(:)
CALL init_us_1
CALL newd
call init_us_1
call newd
ALLOCATE (aux(nrxx))
allocate (aux(nrxx))
call allocate_bec_type ( nkb, nbnd, becp )
! four times npwx should be enough
ALLOCATE (idx (4*npwx) )
ALLOCATE (igtog (4*npwx) )
if (lda_plus_u) allocate(v_h_new(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat))
if (dft_is_meta()) then
allocate (idx (4*npwx) )
allocate (igtog (4*npwx) )
if(lda_plus_u) allocate(v_h_new(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat))
if(dft_is_meta())then
allocate (kedtaur_new(nrxx,nspin))
else
allocate (kedtaur_new(1,nspin))
@ -69,97 +69,97 @@ SUBROUTINE write_casino_pwfn(gather)
idx(:) = 0
igtog(:) = 0
IF( lsda ) THEN
if( lsda )then
nbndup = nbnd
nbnddown = nbnd
nk = nks/2
! nspin = 2
ELSE
else
nbndup = nbnd
nbnddown = 0
nk = nks
! nspin = 1
ENDIF
endif
ek = 0.d0
eloc= 0.d0
enl = 0.d0
demet=0.d0
!
DO ispin = 1, nspin
do ispin = 1, nspin
!
! calculate the local contribution to the total energy
!
! bring rho to G-space
!
aux(:) = CMPLX( rho%of_r(:,ispin), 0.d0,kind=DP)
CALL cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
aux(:) = cmplx( rho%of_r(:,ispin), 0.d0,kind=DP)
call cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
!
DO nt=1,ntyp
DO ig = 1, ngm
do nt=1,ntyp
do ig = 1, ngm
eloc = eloc + vloc(igtongl(ig),nt) * strf(ig,nt) &
* CONJG(aux(nl(ig)))
ENDDO
ENDDO
* conjg(aux(nl(ig)))
enddo
enddo
DO ik = 1, nk
do ik = 1, nk
ikk = ik + nk*(ispin-1)
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio (evc, nwordwfc, iunwfc, ikk, - 1)
CALL init_us_2 (npw, igk, xk (1, ikk), vkb)
CALL calbec ( npw, vkb, evc, becp )
call gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
call davcio (evc, nwordwfc, iunwfc, ikk, - 1)
call init_us_2 (npw, igk, xk (1, ikk), vkb)
call calbec ( npw, vkb, evc, becp )
!
! -TS term for metals (if any)
! -TS term for metals (ifany)
!
IF ( degauss > 0.0_dp) THEN
DO ibnd = 1, nbnd
if( degauss > 0.0_dp)then
do ibnd = 1, nbnd
demet = demet + wk (ik) * &
degauss * w1gauss ( (ef-et(ibnd,ik)) / degauss, ngauss)
END DO
END IF
enddo
endif
!
DO ig =1, npw
IF( igk(ig) > 4*npwx ) &
CALL errore ('pw2casino','increase allocation of index', ig)
do ig =1, npw
if( igk(ig) > 4*npwx ) &
call errore ('pw2casino','increase allocation of index', ig)
idx( igk(ig) ) = 1
ENDDO
enddo
!
! calculate the kinetic energy
!
DO ibnd = 1, nbnd
DO j = 1, npw
ek = ek + CONJG(evc(j,ibnd)) * evc(j,ibnd) * &
do ibnd = 1, nbnd
do j = 1, npw
ek = ek + conjg(evc(j,ibnd)) * evc(j,ibnd) * &
g2kin(j) * wg(ibnd,ikk)
END DO
enddo
!
! Calculate Non-local energy
!
ijkb0 = 0
DO nt = 1, ntyp
DO na = 1, nat
IF (ityp (na) .EQ.nt) THEN
DO ih = 1, nh (nt)
do nt = 1, ntyp
do na = 1, nat
if(ityp (na) .eq. nt)then
do ih = 1, nh (nt)
ikb = ijkb0 + ih
enl=enl+CONJG(becp%k(ikb,ibnd))*becp%k(ikb,ibnd) &
enl=enl+conjg(becp%k(ikb,ibnd))*becp%k(ikb,ibnd) &
*wg(ibnd,ikk)* dvan(ih,ih,nt)
DO jh = ( ih + 1 ), nh(nt)
do jh = ( ih + 1 ), nh(nt)
jkb = ijkb0 + jh
enl=enl + &
(CONJG(becp%k(ikb,ibnd))*becp%k(jkb,ibnd)+&
CONJG(becp%k(jkb,ibnd))*becp%k(ikb,ibnd))&
(conjg(becp%k(ikb,ibnd))*becp%k(jkb,ibnd)+&
conjg(becp%k(jkb,ibnd))*becp%k(ikb,ibnd))&
* wg(ibnd,ikk) * dvan(ih,jh,nt)
END DO
enddo
ENDDO
enddo
ijkb0 = ijkb0 + nh (nt)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
endif
enddo
enddo
enddo
enddo
enddo
#ifdef __PARA
call mp_sum( eloc, intra_pool_comm )
@ -172,12 +172,12 @@ SUBROUTINE write_casino_pwfn(gather)
ek = ek * tpiba2
ngtot = 0
DO ig = 1, 4*npwx
IF( idx(ig) >= 1 ) THEN
do ig = 1, 4*npwx
if( idx(ig) >= 1 )then
ngtot = ngtot + 1
igtog(ngtot) = ig
ENDIF
ENDDO
endif
enddo
!
! compute ewald contribution
@ -187,167 +187,167 @@ SUBROUTINE write_casino_pwfn(gather)
!
! compute hartree and xc contribution
!
CALL v_of_rho( rho, rho_core, rhog_core, &
call v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vnew )
!
etot=(ek + (etxc-etxcc)+ehart+eloc+enl+ewld)+demet
!
IF (ionode.or..not.gather) THEN
if(ionode.or..not.gather)then
io = 77
WRITE (6,'(/,5x,''Writing file pwfn.data for program CASINO'')')
write (6,'(/,5x,''Writing file pwfn.data for program CASINO'')')
CALL seqopn( 77, 'pwfn.data', 'formatted',exst)
call seqopn( 77, 'pwfn.data', 'formatted',exst)
CALL write_header
ENDIF
call write_header
endif
ALLOCATE ( g_l(3,ngtot), evc_l(ngtot) )
DO ig = 1, ngtot
allocate ( g_l(3,ngtot), evc_l(ngtot) )
do ig = 1, ngtot
g_l(:,ig) = g(:,igtog(ig))
ENDDO
enddo
IF(gather)THEN
ALLOCATE ( ngtot_d(nproc_pool), ngtot_cumsum(nproc_pool) )
CALL mp_gather( ngtot, ngtot_d, ionode_id, intra_pool_comm )
CALL mp_bcast( ngtot_d, ionode_id, intra_pool_comm )
if(gather)then
allocate ( ngtot_d(nproc_pool), ngtot_cumsum(nproc_pool) )
call mp_gather( ngtot, ngtot_d, ionode_id, intra_pool_comm )
call mp_bcast( ngtot_d, ionode_id, intra_pool_comm )
id = 0
DO ip = 1,nproc_pool
do ip = 1,nproc_pool
ngtot_cumsum(ip) = id
id = id + ngtot_d(ip)
ENDDO
enddo
ngtot_g = id
ALLOCATE ( g_g(3,ngtot_g), evc_g(ngtot_g) )
CALL mp_gather( g_l, g_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
allocate ( g_g(3,ngtot_g), evc_g(ngtot_g) )
call mp_gather( g_l, g_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
IF (ionode) THEN
ALLOCATE ( indx(ngtot_g) )
CALL create_index2(g_g,indx)
CALL write_gvecs(g_g,indx)
ENDIF
ELSE
ALLOCATE ( indx(ngtot) )
CALL create_index2(g_l,indx)
CALL write_gvecs(g_l,indx)
ENDIF
if(ionode)then
allocate ( indx(ngtot_g) )
call create_index2(g_g,indx)
call write_gvecs(g_g,indx)
endif
else
allocate ( indx(ngtot) )
call create_index2(g_l,indx)
call write_gvecs(g_l,indx)
endif
IF (ionode.or..not.gather)CALL write_wfn_head
if(ionode.or..not.gather)call write_wfn_head
DO ik = 1, nk
IF (ionode.or..not.gather)CALL write_kpt_head
do ik = 1, nk
if(ionode.or..not.gather)call write_kpt_head
DO ispin = 1, nspin
do ispin = 1, nspin
ikk = ik + nk*(ispin-1)
IF( nks > 1 ) THEN
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio(evc,nwordwfc,iunwfc,ikk,-1)
ENDIF
DO ibnd = 1, nbnd
IF (ionode.or..not.gather)CALL write_bnd_head
if( nks > 1 )then
call gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
call davcio(evc,nwordwfc,iunwfc,ikk,-1)
endif
do ibnd = 1, nbnd
if(ionode.or..not.gather)call write_bnd_head
evc_l(:) = (0.d0, 0d0)
DO ig=1, ngtot
do ig=1, ngtot
! now for all G vectors find the PW coefficient for this k-point
find_ig: DO ig7 = 1, npw
IF( igk(ig7) == igtog(ig) )THEN
find_ig: do ig7 = 1, npw
if( igk(ig7) == igtog(ig) )then
evc_l(ig) = evc(ig7,ibnd)
EXIT find_ig
ENDIF
ENDDO find_ig
ENDDO
IF(gather)THEN
CALL mp_gather( evc_l, evc_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
exit find_ig
endif
enddo find_ig
enddo
if(gather)then
call mp_gather( evc_l, evc_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
IF (ionode)CALL write_wfn_data(evc_g,indx)
ELSE
CALL write_wfn_data(evc_l,indx)
ENDIF
ENDDO
ENDDO
ENDDO
IF (ionode.or..not.gather)CLOSE(io)
if(ionode)call write_wfn_data(evc_g,indx)
else
call write_wfn_data(evc_l,indx)
endif
enddo
enddo
enddo
if(ionode.or..not.gather)close(io)
WRITE (stdout,*) 'Kinetic energy ' , ek/e2
WRITE (stdout,*) 'Local energy ', eloc/e2
WRITE (stdout,*) 'Non-Local energy ', enl/e2
WRITE (stdout,*) 'Ewald energy ', ewld/e2
WRITE (stdout,*) 'xc contribution ',(etxc-etxcc)/e2
WRITE (stdout,*) 'hartree energy ', ehart/e2
IF ( degauss > 0.0_dp ) &
WRITE (stdout,*) 'Smearing (-TS) ', demet/e2
WRITE (stdout,*) 'Total energy ', etot/e2
write (stdout,*) 'Kinetic energy ' , ek/e2
write (stdout,*) 'Local energy ', eloc/e2
write (stdout,*) 'Non-Local energy ', enl/e2
write (stdout,*) 'Ewald energy ', ewld/e2
write (stdout,*) 'xc contribution ',(etxc-etxcc)/e2
write (stdout,*) 'hartree energy ', ehart/e2
if( degauss > 0.0_dp ) &
write (stdout,*) 'Smearing (-TS) ', demet/e2
write (stdout,*) 'Total energy ', etot/e2
DEALLOCATE (igtog)
DEALLOCATE (idx)
deallocate (igtog)
deallocate (idx)
call deallocate_bec_type (becp)
DEALLOCATE (aux)
deallocate (aux)
DEALLOCATE ( g_l, evc_l )
IF(gather) DEALLOCATE ( ngtot_d, ngtot_cumsum, g_g, evc_g )
IF(ionode.or..not.gather) DEALLOCATE (indx)
deallocate ( g_l, evc_l )
if(gather) deallocate ( ngtot_d, ngtot_cumsum, g_g, evc_g )
if(ionode.or..not.gather) deallocate (indx)
CONTAINS
SUBROUTINE write_header
WRITE(io,'(a)') title
WRITE(io,'(a)')
WRITE(io,'(a)') ' BASIC INFO'
WRITE(io,'(a)') ' ----------'
WRITE(io,'(a)') ' Generated by:'
WRITE(io,'(a)') ' PWSCF'
WRITE(io,'(a)') ' Method:'
WRITE(io,'(a)') ' DFT'
WRITE(io,'(a)') ' DFT Functional:'
WRITE(io,'(a)') ' unknown'
WRITE(io,'(a)') ' Pseudopotential'
WRITE(io,'(a)') ' unknown'
WRITE(io,'(a)') ' Plane wave cutoff (au)'
WRITE(io,*) ecutwfc/2
WRITE(io,'(a)') ' Spin polarized:'
WRITE(io,*)lsda
IF ( degauss > 0.0_dp ) THEN
WRITE(io,'(a)') ' Total energy (au per primitive cell; includes -TS term)'
WRITE(io,*)etot/e2, demet/e2
ELSE
WRITE(io,'(a)') ' Total energy (au per primitive cell)'
WRITE(io,*)etot/e2
END IF
WRITE(io,'(a)') ' Kinetic energy (au per primitive cell)'
WRITE(io,*)ek/e2
WRITE(io,'(a)') ' Local potential energy (au per primitive cell)'
WRITE(io,*)eloc/e2
WRITE(io,'(a)') ' Non local potential energy(au per primitive cel)'
WRITE(io,*)enl/e2
WRITE(io,'(a)') ' Electron electron energy (au per primitive cell)'
WRITE(io,*)ehart/e2
WRITE(io,'(a)') ' Ion ion energy (au per primitive cell)'
WRITE(io,*)ewld/e2
WRITE(io,'(a)') ' Number of electrons per primitive cell'
WRITE(io,*)NINT(nelec)
! uncomment the following if you want the Fermi energy - KN 2/4/09
! WRITE(io,'(a)') ' Fermi energy (au)'
! WRITE(io,*) ef/e2
WRITE(io,'(a)') ' '
WRITE(io,'(a)') ' GEOMETRY'
WRITE(io,'(a)') ' -------- '
WRITE(io,'(a)') ' Number of atoms per primitive cell '
WRITE(io,*) nat
WRITE(io,'(a)')' Atomic number and position of the atoms(au) '
DO na = 1, nat
write(io,'(a)') title
write(io,'(a)')
write(io,'(a)') ' BASIC INFO'
write(io,'(a)') ' ----------'
write(io,'(a)') ' Generated by:'
write(io,'(a)') ' PWSCF'
write(io,'(a)') ' Method:'
write(io,'(a)') ' DFT'
write(io,'(a)') ' DFT Functional:'
write(io,'(a)') ' unknown'
write(io,'(a)') ' Pseudopotential'
write(io,'(a)') ' unknown'
write(io,'(a)') ' Plane wave cutoff (au)'
write(io,*) ecutwfc/2
write(io,'(a)') ' Spin polarized:'
write(io,*)lsda
if( degauss > 0.0_dp )then
write(io,'(a)') ' Total energy (au per primitive cell; includes -TS term)'
write(io,*)etot/e2, demet/e2
else
write(io,'(a)') ' Total energy (au per primitive cell)'
write(io,*)etot/e2
endif
write(io,'(a)') ' Kinetic energy (au per primitive cell)'
write(io,*)ek/e2
write(io,'(a)') ' Local potential energy (au per primitive cell)'
write(io,*)eloc/e2
write(io,'(a)') ' Non local potential energy(au per primitive cel)'
write(io,*)enl/e2
write(io,'(a)') ' Electron electron energy (au per primitive cell)'
write(io,*)ehart/e2
write(io,'(a)') ' Ion ion energy (au per primitive cell)'
write(io,*)ewld/e2
write(io,'(a)') ' Number of electrons per primitive cell'
write(io,*)nint(nelec)
! uncomment the following ifyou want the Fermi energy - KN 2/4/09
! write(io,'(a)') ' Fermi energy (au)'
! write(io,*) ef/e2
write(io,'(a)') ' '
write(io,'(a)') ' GEOMETRY'
write(io,'(a)') ' -------- '
write(io,'(a)') ' Number of atoms per primitive cell '
write(io,*) nat
write(io,'(a)')' Atomic number and position of the atoms(au) '
do na = 1, nat
nt = ityp(na)
at_num = atomic_number(TRIM(atm(nt)))
WRITE(io,'(i6,3f20.14)') at_num, (alat*tau(j,na),j=1,3)
ENDDO
WRITE(io,'(a)') ' Primitive lattice vectors (au) '
WRITE(io,100) alat*at(1,1), alat*at(2,1), alat*at(3,1)
WRITE(io,100) alat*at(1,2), alat*at(2,2), alat*at(3,2)
WRITE(io,100) alat*at(1,3), alat*at(2,3), alat*at(3,3)
WRITE(io,'(a)') ' '
at_num = atomic_number(trim(atm(nt)))
write(io,'(i6,3f20.14)') at_num, (alat*tau(j,na),j=1,3)
enddo
write(io,'(a)') ' Primitive lattice vectors (au) '
write(io,100) alat*at(1,1), alat*at(2,1), alat*at(3,1)
write(io,100) alat*at(1,2), alat*at(2,2), alat*at(3,2)
write(io,100) alat*at(1,3), alat*at(2,3), alat*at(3,3)
write(io,'(a)') ' '
100 FORMAT (3(1x,f20.15))
100 format (3(1x,f20.15))
END SUBROUTINE
END SUBROUTINE write_header
SUBROUTINE write_gvecs(g,indx)
@ -355,47 +355,47 @@ CONTAINS
INTEGER,INTENT(in) :: indx(:)
INTEGER ig
WRITE(io,'(a)') ' G VECTORS'
WRITE(io,'(a)') ' ---------'
WRITE(io,'(a)') ' Number of G-vectors'
WRITE(io,*) SIZE(g,2)
WRITE(io,'(a)') ' Gx Gy Gz (au)'
DO ig = 1, SIZE(g,2)
WRITE(io,100) tpi/alat*g(1,indx(ig)), tpi/alat*g(2,indx(ig)), &
write(io,'(a)') ' G VECTORS'
write(io,'(a)') ' ---------'
write(io,'(a)') ' Number of G-vectors'
write(io,*) size(g,2)
write(io,'(a)') ' Gx Gy Gz (au)'
do ig = 1, size(g,2)
write(io,100) tpi/alat*g(1,indx(ig)), tpi/alat*g(2,indx(ig)), &
tpi/alat*g(3,indx(ig))
ENDDO
enddo
100 FORMAT (3(1x,f20.15))
100 format (3(1x,f20.15))
WRITE(io,'(a)') ' '
END SUBROUTINE
write(io,'(a)') ' '
END SUBROUTINE write_gvecs
SUBROUTINE write_wfn_head
WRITE(io,'(a)') ' WAVE FUNCTION'
WRITE(io,'(a)') ' -------------'
WRITE(io,'(a)') ' Number of k-points'
WRITE(io,*) nk
END SUBROUTINE
write(io,'(a)') ' WAVE FUNCTION'
write(io,'(a)') ' -------------'
write(io,'(a)') ' Number of k-points'
write(io,*) nk
END SUBROUTINE write_wfn_head
SUBROUTINE write_kpt_head
WRITE(io,'(a)') ' k-point # ; # of bands (up spin/down spin); &
write(io,'(a)') ' k-point # ; # of bands (up spin/down spin); &
& k-point coords (au)'
WRITE(io,'(3i4,3f20.16)') ik, nbndup, nbnddown, &
write(io,'(3i4,3f20.16)') ik, nbndup, nbnddown, &
(tpi/alat*xk(j,ik),j=1,3)
END SUBROUTINE
END SUBROUTINE write_kpt_head
SUBROUTINE write_bnd_head
! KN: if you want to print occupancies, replace these two lines ...
WRITE(io,'(a)') ' Band, spin, eigenvalue (au)'
WRITE(io,*) ibnd, ispin, et(ibnd,ikk)/e2
! KN: ifyou want to print occupancies, replace these two lines ...
write(io,'(a)') ' Band, spin, eigenvalue (au)'
write(io,*) ibnd, ispin, et(ibnd,ikk)/e2
! ...with the following two - KN 2/4/09
! WRITE(io,'(a)') ' Band, spin, eigenvalue (au), occupation number'
! WRITE(io,*) ibnd, ispin, et(ibnd,ikk)/e2, wg(ibnd,ikk)/wk(ikk)
WRITE(io,'(a)') ' Eigenvectors coefficients'
END SUBROUTINE
! write(io,'(a)') ' Band, spin, eigenvalue (au), occupation number'
! write(io,*) ibnd, ispin, et(ibnd,ikk)/e2, wg(ibnd,ikk)/wk(ikk)
write(io,'(a)') ' Eigenvectors coefficients'
END SUBROUTINE write_bnd_head
SUBROUTINE write_wfn_data(evc,indx)
@ -403,22 +403,22 @@ CONTAINS
INTEGER,INTENT(in) :: indx(:)
INTEGER ig
DO ig=1, SIZE(evc,1)
WRITE(io,*)evc(indx(ig))
END DO
END SUBROUTINE
do ig=1, size(evc,1)
write(io,*)evc(indx(ig))
enddo
END SUBROUTINE write_wfn_data
SUBROUTINE create_index2(y,x_index)
DOUBLE PRECISION,INTENT(in) :: y(:,:)
INTEGER,INTENT(out) :: x_index(SIZE(y,2))
DOUBLE PRECISION y2(SIZE(y,2))
INTEGER,INTENT(out) :: x_index(size(y,2))
DOUBLE PRECISION y2(size(y,2))
INTEGER i
DO i = 1,SIZE(y,2)
do i = 1,size(y,2)
y2(i) = sum(y(:,i)**2)
ENDDO
CALL create_index(y2,x_index)
END SUBROUTINE
enddo
call create_index(y2,x_index)
END SUBROUTINE create_index2
SUBROUTINE create_index(y,x_index)
@ -449,7 +449,7 @@ jloop : do j=l+1,ir
if(y(x_index(i))<=yj)then
x_index(i+1)=x_indexj
cycle jloop
endif ! y(x_index(i))<=yj
endif! y(x_index(i))<=yj
x_index(i+1)=x_index(i)
enddo ! i
x_index(l)=x_indexj
@ -493,7 +493,7 @@ jloop : do j=l+1,ir
if(jstack>stacksize)then
write(6,*)'stacksize is too small.'
stop
endif ! jstack>stacksize
endif! jstack>stacksize
if(ir-i+1>=j-l)then
istack(jstack)=ir
istack(jstack-1)=i
@ -502,8 +502,8 @@ jloop : do j=l+1,ir
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif ! ir-i+1>=j-l
endif ! ir-l<ins_sort_thresh
endif! ir-i+1>=j-l
endif! ir-l<ins_sort_thresh
enddo
END SUBROUTINE create_index