********************************************************************************
*                                                                              *
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### ortho ###
c
c     This is a orthogonalisation routine. 
c
      subroutine ortho (iorb1,psi,f4,wgt2,wk0,iprt)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension psi(*),f4(*),wgt2(*),wk0(*)
      dimension ovla(60),istp(60)

      jor=0
      ibeg1 = i1b (iorb1)
c ! orbitals are stored in reverse order
      do iorb2=1,norb
         iorb3=norb-iorb2+1
         if (iorb3.le.iorb1) go to 101
         
         ngrid= min (i1si(iorb1),i1si(iorb3))	  	    
         if (mgx(6,iorb3).ne.mgx(6,iorb1)) go to 101
         if (ige(iorb3).ne.ige(iorb1)) go to 101
         
         jor=jor+1
         ibeg3= i1b (iorb3 )
         istp(jor)=ibeg3
         
         call prod2 (ngrid,psi(ibeg1),psi(ibeg3),wk0)
         call prod  (ngrid,f4,wk0)
         
         ovla(jor)=ddot (ngrid,wgt2,1,wk0,1)

c        write final overlap
         
         if (iprt.ne.0) then
c           write(*,60) iorn(iorb1),bond(iorb1),gut(iorb1),
c           &                   iorn(iorb3),bond(iorb3),gut(iorb3),ovla(jor)
c           
            write(*,60) iorn(iorb1),bond(iorb1),gut(iorb1),
     &           iorn(iorb3),bond(iorb3),gut(iorb3),ovla(jor)
 60         format(4x,' <',i2,1x,a5,1x,a1,1x,'|',
     &           1x,i2,1x,a5,1x,a1,1x,'> = ',d10.2)     
         endif
 101     continue
      enddo

      if (iprt.eq.2) return

c     start orthogonalisation

      if (jor.eq.0) return

      ano=0.d0
      do ira=1,jor
         ngrid= min (i1si(iorb1),i1si(ira))	  	    	
         ano=ano+ovla(ira)*ovla(ira)
         call daxpy (ngrid,-ovla(ira),psi(istp(ira)),1,psi(ibeg1),1)
      enddo

c     normalize the orthogonalized orbital isa  but be sure the
c     unorthogonalized isa was normalized otherwise use <isa/isa>
c     instead of 1. in ano

      ano=sqrt(1.d0+ano)/(1.d0-ano*ano)
      ngrid = i1si(iorb1)
      call dscal (ngrid,ano,psi(ibeg1),1)

      return
      end