C
c
c
      subroutine directionCGM(casco_iter, DomegaDrho, direction, 
     &   nr, ntbc)
c
      INCLUDE 'CASCO.DIM'
c
c  Arrays Passed: 
c     DomegaDrho is given as current gradient; returned unchanged
c     direction  is given as previous direction; it is updated
c
      real DomegaDrho(NRDIM,NTBCDIM), direction(NRDIM,NTBCDIM)
      integer casco_iter, nr, ntbc
c
c  Local Arrays: 
c     DomegaDrho_0 is saved as previous value of the gradient
c
      real DomegaDrho_0(NRDIM,NTBCDIM)
      SAVE DomegaDrho_0

      real gamma, top, bot
      integer i, j
c
c   First time through, use steepest descent
c
      IF(casco_iter.eq.1) THEN
c
         do i=1,nr
         do j=1,ntbc
            direction(i,j)=-DomegaDrho(i,j)
         enddo
         enddo
c
      ELSE
c      
c  Once started, use CGM
C
         top = 0.0
         bot = 0.0
         do i=1,nr
         do j=1,ntbc  
            top=top+DomegaDrho(i,j)*(DomegaDrho(i,j)-DomegaDrho_0(i,j))
            bot=bot+DomegaDrho_0(i,j)*DomegaDrho_0(i,j)         
         enddo
         enddo 
         gamma=top/bot
c  
         do i=1,nr
         do j=1,ntbc  
            direction(i,j)=-DomegaDrho(i,j) + gamma*direction(i,j)
         enddo
         enddo
c
      ENDIF
c
c  Direction is complete;  return
c  
      do i=1,nr
      do j=1,ntbc  
        DomegaDrho_0(i,j) = DomegaDrho(i,j)
      enddo
      enddo
C    
      RETURN
      END
                 

