      SUBROUTINE POWSAV(HESS, GRAD, XPARAM, PMAT, ILOOP, BMAT, IPOW)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension HESS(MAXPAR,*),GRAD(*),BMAT(MAXPAR,*),IPOW(9),
     1 XPARAM(*), PMAT(*)
**********************************************************************
*
* POWSAV STORES AND RESTORES DATA USED IN THE SIGMA GEOMETRY
*        OPTIMISATION.
*
*  ON INPUT HESS   = HESSIAN MATRIX, PARTIAL OR WHOLE.
*           GRAD   = GRADIENTS.
*           XPARAM = CURRENT STATE OF PARAMETERS.
*           ILOOP  = INDEX OF HESSIAN, OR FLAG OF POINT REACHED SO-FAR.
*           BMAT   = "B" MATRIX!
*           IPOW   = INDICES AND FLAGS.
*           IPOW(9)= 0 FOR RESTORE, 1 FOR DUMP
*
**********************************************************************
      common /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, DUMY(MAXPAR)
      common /ERRFN / ERRFN(MAXPAR), AICORR(MAXPAR)
      common /ELEMTS/ ELEMNT(107)
      common /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),
     1                     LOCDEP(MAXPAR)
      common /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      common /GEOM  / GEO(3,NUMATM)
      common /LOCVAR/ LOCVAR(2,MAXPAR)
      common /KEYWRD/ KEYWRD
      common /VALVAR/ VALVAR(MAXPAR),NUMVAR
      dimension COORD(3,NUMATM)
      character ELEMNT*2, KEYWRD*241, GETNAM*80
      common /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
      common /ALPARM/ ALPARM(3,MAXPAR),X0, X1, X2, JLOOP
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /PATH  / LATOM,LPARAM,REACT(200)
      open(unit=MIOC,file=GETNAM('FOR009'),
     +     status='UNKNOWN',form='UNFORMATTED')
      rewind MIOC
      open(unit=MIOD,file=GETNAM('FOR010'),
     +     status='UNKNOWN',form='UNFORMATTED')
      rewind MIOD
      IR=MIOC
      if (IPOW(9) .NE. 0) then
         if (IPOW(9) .EQ. 1) then
            write(MIOB,
     %        '(//10X,''- - - - - - - TIME UP - - - - - - -'',//)')
            write(MIOB,
     %        '(//10X,'' - THE CALCULATION IS BEING DUMPED TO DISK''/
     %            10X,''   RESTART IT USING THE KEY-WORD "RESTART"'')')
            FUNCT1=SQRT(DOT(GRAD,GRAD,NVAR))
            write(MIOB,'(//10X,''CURRENT VALUE OF GRADIENT NORM =''
     %                        ,F12.6)')  FUNCT1
            do 10 I=1,NVAR
               K=LOC(1,I)
               L=LOC(2,I)
   10       GEO(L,K)=XPARAM(I)
            write(MIOB,'(/10X,''CURRENT VALUE OF GEOMETRY'',/)')
            DEGREE=57.29577951D0
            if (NA(1) .EQ. 99) then
C
C  CONVERT FROM CARTESIAN COORDINATES TO INTERNAL
C
               do 20 I=1,NATOMS
                  do 20 J=1,3
   20          COORD(J,I)=GEO(J,I)
               call XYZINT(COORD,NUMAT,NA,NB,NC,1.D0,GEO)
            end if
            call GEOUT(MIOB)
         end if
         write(IR)IPOW,ILOOP
         write(IR) (XPARAM(I),I=1,NVAR)
         write(IR) (  GRAD(I),I=1,NVAR)
         write(IR) ((HESS(J,I),J=1,NVAR),I=1,NVAR)
         write(IR) ((BMAT(J,I),J=1,NVAR),I=1,NVAR)
         LINEAR=(NVAR*(NVAR+1))/2
         write(IR) (PMAT(I),I=1,LINEAR)
         if (INDEX(KEYWRD,'AIDER').NE.0) write(IR) (AICORR(I),I=1,NVAR)
         LINEAR=(NORBS*(NORBS+1))/2
         write(MIOD) (PA(I),I=1,LINEAR)
         if (NALPHA.NE.0) write(MIOD) (PB(I),I=1,LINEAR)
         if (LATOM .NE. 0) then
            write(IR) ((ALPARM(J,I),J=1,3),I=1,NVAR)
            write(IR)JLOOP,X0, X1, X2
         end if
         close (unit=MIOC)
         close (unit=MIOD)
         return
      else
         write(MIOB,'(//10X,'' RESTORING DATA FROM DISK''/)')
         read(IR)IPOW,ILOOP
         read(IR) (XPARAM(I),I=1,NVAR)
         read(IR) (  GRAD(I),I=1,NVAR)
         read(IR) ((HESS(J,I),J=1,NVAR),I=1,NVAR)
         read(IR) ((BMAT(J,I),J=1,NVAR),I=1,NVAR)
         FUNCT1=SQRT(DOT(GRAD,GRAD,NVAR))
         write(MIOB,'(10X,''FUNCTION ='',F13.6//)') FUNCT1
         LINEAR=(NVAR*(NVAR+1))/2
         read(IR) (PMAT(I),I=1,LINEAR)
         if (INDEX(KEYWRD,'AIDER').NE.0) read(IR) (AICORR(I),I=1,NVAR)
         LINEAR=(NORBS*(NORBS+1))/2
         read(MIOD) (PA(I),I=1,LINEAR)
         if (NALPHA.NE.0) read(MIOD) (PB(I),I=1,LINEAR)
         if (LATOM.NE.0) then
            read(IR) ((ALPARM(J,I),J=1,3),I=1,NVAR)
            read(IR)JLOOP,X0, X1, X2
            ILOOP=ILOOP+1
         end if
         ILOOP=ILOOP+1
         return
      end if
      end
