C#######################################################################
C#########                                                ##############
C#########        PROGRAM ASYRMO                          ##############
C#########                                                ##############
C
C    VERSION: JULY 1991
C
C    DISTRIBUTED AT THE OAK RIDGE THEORY WORKSHOP, AUG 1991, FOR
C    USE WITH OTHER MODIFIED OSCILLATOR CODES, GAMPN (CALCULATES
C    THE DEFORMED SINGLE PARTICLE ORBITALS, THEIR EIGENVALUES,
C    EIGENVECTORS, AND THE S. P. MATRIX ELEMENTS NEEDED FOR ASYRMO,
C    EG, <J+>, ETC) AND PROBAMO (COMPUTES THE M1/E2 MATRIX ELEMENTS
C    FROM THE PARTICLE + ROTOR EIGENVECTORS GENERATED IN ASYRMO).
C
C
C    ASYRMO DIAGONALIZES THE PARTICLE + TRIAXIAL ROTOR HAMILTONIAN
C    IN THE STRONG COUPLING BASIS, WITH THE S. P. MATRIX ELEMENTS
C    EXPRESSED IN DEFORMED SCHEME. ORIGINAL PROGRAM DESCRIBED IN
C    NUCL. PHYS. A307 (1978) 189. THIS VERSION ACCEPTS INPUT FROM
C    GAMPN (MODIFIED OSCILLATOR POTENTIAL).
C
C    SINGLE-PARTICLE ENERGIES AND VARIOUS MATRIX ELEMENTS AS
C    CALCULATED IN GAMPN ARE READ FROM FILES 16 AND 17. OUPUT DATA
C    FOR M1/E2 MATRIX ELEMENTS, ETC (PROBAMO) WRITTEN ON FILE18
C
C***********************************************************************
C
C    NUMEROUS MODIFICATIONS FROM THE ORIGINAL DESCRIPTION HAVE BEEN
C    MADE, INCLUDING:
C          OPTIONAL TWO-BODY TREATMENT OF THE RECOIL TERMS
C          K-TRUNCTATION
C          EXTENSIONS TO HIGH-SPIN (BUT STILL ONLY 1-QP STATES)
C          VARIOUS ANGULAR MOMENTUM EXPECTATION VALUES CALCULATED
C                  FOR THE ENERGY EIGENVECTORS
C          POSSIBLE VMI (VARIABLE MOMENT OF INERTIA) TREATMENT OF
C                  EVEN-EVEN CORE SPECTRUM (BY AUSTYN GRIFFITHS)
C
C   THE BASIC TWO-BODY TREATMENT OF THE RECOIL TERMS (J+J-, JZ**2, ETC)
C   IS DESCRIBED IN NUCL PHYS A253 (1975) 45. HOWEVER, NOTE THAT THE
C   KRONECKER DELTAS IN THEIR EQN 14 SHOULD BE OMITTED. ALSO, NOTE THAT
C   A SMALL, CONSTANT "CORE CONTRIBUTION" THAT ARISES IN THE TWO-BODY
C   TREATMENT HAS BEEN OMITTED HERE. THIS DOES NOT AFFECT EXCITATION
C   ENERGIES, BUT CAN CAUSE SMALL ERRORS IN THE EXPECTATION VALUES
C   <I.j>, ETC (IN THE SUBROUTINE "EXPECT") ESPECIALLY NEAR BANDHEADS.
C
C   THE VMI TREATMENT FOR THE CORE ENERGY SPECTRUM IS INCORPORATED
C   BY THE METHOD OF NUCL PHYS A253 (1975) 231. THE VMI CORE
C   PARAMETERIZATION IS ALSO GIVEN THERE (INCLUDING TRIAXIAL SHAPES),
C   AND FOR AXIALLY SYMMETRIC SHAPES, THIS IS EQUIVALENT TO THE USUAL
C   METHODS (PHYS REV 178 (1969) 1864)
C
C***********************************************************************
C
C  INPUT DATA (ENERGIES IN MEV):
C
C     1. FILE16, FILE17, FILE18  (FREE FORMAT)
C          FILE16, FILE17 WERE WRITTEN BY GAMPN AND ARE INPUT FILES.
C                  THESE CONTAIN THE S.P. WAVEFUNCTIONS (FILE16, NEEDED
C                  HERE IF THE VMI OPTION IS EXERCISED) AND VARIOUS S.P.
C                 QUANTITIES NEEDED HERE (EG, MATRIX ELEMENTS of j+ ETC)
C          FILE18 : NAME OF THE OUTPUT FILE, NEEDED FOR PROBAMO
C
C     2. IPKT, ISKIP  (FREE FORMAT)
C          IPKT: NUMBER OF DATA SETS (OF GIVEN DEFORMATION AND PARITY)
C                TO BE READ FROM FILE17
C          ISKIP: THE FIRST ISKIP SETS ARE NOT USED
C
C     3. ISTRCH, IREC  (FREE FORMAT)
C          ISTRCH=0: CALCULATION IN SPHERICAL COORDINATES
C                    (DEF PARAM DELTA2, GAMMAD, DELTA4, DELTA6)
C          ISTRCH=1: CALCULATION IN STRETCHED COORDINATES (STANDARD)
C                    (DEF PARAM EPS2, GAMMA, EPS4, EPS6)
C
C          IREC=1:   TREAT THE RECOIL TERMS (J+J-, J+J+, JZ2, ETC) AS
C                    TWO BODY TERMS, IE, CONSTRUCT THEM FROM MATRIX
C                    ELEMENTS OF J+, J-, JZ SUMMED OVER A COMPLETE SET
C                    OF INTERMEDIATE STATES (INCLUDING 1-QP AND 3-QP
C                    STATES), EACH TERM WITH ITS CORRECT PAIRING FACTOR
C          IREC=0:   TREAT THESE TERMS AS ONE-BODY OPERATORS, AS IN
C                    NPA307, 189 WITH SINGLE PARTICLE MATRIX ELEMENTS
C                    MULTIPLIED BY UU'+VV'
C
C    *** RECOMMENDATION: THE IREC=0 OPTION IS SUGGESTED. THE TWO
C    ***                 APPROACHES SEEM TO GIVE VERY SIMILAR
C    ***                 RESULTS, AND THE ONE-BODY TREATMENT IS
C    ***                 SIMPLER. NOTE THAT IF YOU USE IREC=1 HERE,
C    ***                 IREC=1 MUST HAVE BEEN SELECTED IN GAMPN.
C
C     4. VMI,NMIN,NMAX,IARCUT,A00,STIFF  (FREE FORMAT)
C          VMI=1 : CALCULATION OF CORE STATES WITH VMI
C          NMIN,NMAX  : MINIMUM/MAXIMUM N-SHELL IN DEFORMED S.P.
C                       WAVEFUNCTIONS
C          IARCUT: NUMBER OF CORE STATES INCLUDED FOR EACH (CORE) SPIN
C          A00,STIFF : VMI PARAMETERS (MEV AND MEV**3)
C                      A00=1/(2*J0), STIFF=C=1/(2*J1)
C                      WHERE J0 AND C ARE THE VMI PARAMETERS OF
C                      MARISCOTTI ET AL, PHYS REV 178 (1969) 1864.
C                      ALSO, J0 AND J1 ARE THE USUAL HARRIS PARAMETERS,
C                      OBTAINED BY FITTING J(MOM OF INERTIA) VS. hw**2
C                      SOFTNESS=1/(2*STIFF**THETA_00**3)
C
C         STANDARD: VMI=0. FOR MOST APPLICATIONS, THE RIGID ROTOR CORE
C                   IS SUFFICIENT.
C
C     5. Z, AA, IMIN, ISPIN, KMAX, E2PLUS, E2PLUR  (FREE FORMAT)
C          Z: PROTON NUMBER OF THE ODD-MASS NUCLEUS, MAY BE EVEN OR ODD.
C          AA: MASS NUMBER OF THE ODD NUCLEUS
C
C          IMIN:  MINIMUM SPIN CONSIDERED IS IMIN/2
C          ISPIN: MAXIMUM SPIN CONSIDERED IS ISPIN/2
C            ***  ISPIN .LE. 129  ***
C            ***  (ISPIN-IMIN)/2 .LE. 19
C
C          KMAX: MAX K ALLOWED (FOR TRUNCATION;
C                               KMAX = ISPIN => NO TRUNCATION)
C            ***  KMAX  = ...13, 17, 21, ...  KMAX .LE. 29  ***
C
C          E2PLUS: THE MOMENTS OF INERTIA ARE ASSUMED TO VARY ACCORDING
C                  TO THE FORMULA
C                    J(K) = 4B*EPS2**2*(SIN(GAMMA+2K*PI/3))**2
C               A. E2PLUS=0 ==> (HBAR)**2/B = 1225/AA**(7/3)  (MEV)
C                          (GRODZIN'S RELATION)
C               B. E2PLUS .GT. 0 ==> THE B-VALUE IS FIXED TO GIVE
C                          THIS ENERGY FOR THE LOWEST 2+ OF CORE ROTOR
C               C. E2PLUS .LT. 0 ==> (HBAR**2)/B = (-E2PLUS)/AA**2 (MEV)
C
C          E2PLUR: FOR E2PLUR .GT. 0 AND E2PLUS .GT. 0 THE GAMMA ANGLE
C                  IS CALCULATED TO GIVE THE SECOND 2+ ENERGY OF THE
C                  CORE EQUAL TO E2PLUR
C                  (NOTE THAT THE CORE GAMMA THUS BECOMES DIFFERENT
C                  FROM THE GAMMA OF THE S.-P. ENERGIES. THIS OPTION
C                  IS NOT RECOMMENDED.)
C
C     6. GN0, GN1, IPAIR, CHSI, ETA  (FREE FORMAT)
C          THE PAIRING STRENGTH PARAMETER IS
C            G*A = GN0 +/- GN1*(N-Z)/A ;  (PROTONS/NEUTRONS)
C          WITH  SQRT(IPAIR*Z)/SQRT(IPAIR*N)  ORBITALS INCLUDED ABOVE
C          AND BELOW THE FERMI LEVEL IN THE PAIRING CALCULATIONS
C
C            STANDARD:  GN0 = 19.2, GN1 = 7.4, IPAIR = 15
C                                    (NUCL PHYS A131 (1969) 1)
C                       IF (Z .LT. 60 OR N .LT. 60):
C                         GN0=22.0, GN1=8.0, IPAIR=5
C                             (E. G. PHYS SCRIPTA 29 (1984) ...
C
C          IF IPAIR=99 OR GN0=0. THEN DELTA=GN1/SQRT(AA)
C            STANDARD:  GN1 = 12.
C          IF IPAIR=99, THE FERMI LEVEL IS SET EQUAL TO GN0 IN MEV
C          IF GN0=0., THE FERMI LEVEL IS SET ON TOP OF THE
C                     (Z OR N)/2 + 1  LEVEL
C
C      CHSI, ETA: CORIOLIS ATTENUATION PARAMETER
C        ALL (OFF-DIAGONAL) S.P. MATRIX ELEMENTS ARE MULTIPLIED BY CHSI
C        ETA IS AD HOC POWER ON (UU+VV)  (U AND V PAIR FACTORS)
C        (BUT ETA DOES NOT AFFECT THE PAIRING FACTORS IN THE TWO-BODY
C        TERMS IF IREC = 1)
C
C     *** NOTE THAT THE ANGULAR MOMENTUM EXPECTATION VALUES (COMPUTED
C     *** IN THE SUBROUTINE "EXPECT", AFTER THE ENERGY DIAGONALIZATION)
C     *** ARE CALCULATED CORRECTLY EVEN IF CHSI .NE. 1, IE, THE SCALING
C     *** OF THE S.P. MATRIX ELEMENTS (<j+>, ETC) BY CHSI IS REMOVED
C     *** IN "EXPECT". HOWEVER, IF ETA .NE. 1, THE ATTENUATION OF
C     *** THESE MATRIX ELEMENTS IS NOT REMOVED, BUT REMAINS IN THE
C     *** COMPUTED EXPECTATION VALUES <I.j>, ETC.
C
C     *** RECOMMENDATION: AVOID CORIOLIS ATTENUATION IF POSSIBLE,
C     ***                 BUT IF YOU MUST ATTENUATE, USE CHSI.
C     ***                 CHSI APPROX 0.7 IS COMMON.
C
C
C     7. IPAR1, NYIC, (LEVEL2(J),J=1,NYIC)      FORMAT(A1,I2,15I3)
C          IPAR1: PARITY CONSIDERED, "+" OR "-"
C          NYIC: NUMBER OF ORBITALS CONSIDERED IN THE PARTICLE-ROTOR
C                CALCULATION, NU .LE. 15
C          LEVEL2(J): "NUMBERS" ON THE ORBITALS INCLUDED (CALCULATED
C                     FROM THE "BOTTOM" AND FOR THE PARITY "IPAR"
C          CF. INPUT CARD 7 FOR GAMPN THE SET OF ORBITALS IN
C          LEVEL2(J) MUST BE A SUBSET (OFTEN THE FULL SET) OF
C          ORBITALS LEVEL(J) IN GAMPN
C
C     8. NANTJ(I), I=1,20  (FREE FORMAT)
C         NUMBER OF STATES TO BE CALCULATED FOR I=IMIN/2, ..., IMAX/2
C         AND TO BE WRITTEN TO FILE18, AND THUS PASSED ON TO PROBAMO.
C         PROBAMO CURRENTLY ALLOWS UP TO 60 STATES TOTAL.
C         (FOR SMALL VALUES OF ISPIN, THE LAST NUMBERS ARE REDUNDANT)
C        ********** NANT(J) .LE. 10 ***************
C
C     9. NOUTJ(I), I=1,20  (FREE FORMAT)
C         NUMBER OF STATES AT SPIN I=IMIN/2, ..., IMAX/2
C         FOR WHICH THE WAVE-FUNCTIONS ARE PRINTED
C                      (NOUTJ(I) .LE. NANTJ(I))
C
C    10. IPOUT,IPOUT(I) I=1,20  (FREE FORMAT)
C         IPOUT .NE. 0 :CORE STATES PRINTED
C         IPOUT(I)     :NUMBER OF STATES AT SPIN I=2, 3, ...., 20
C         FOR WHICH THE CORE TRIAXIAL ROTOR WAVE-FUNCTIONS ARE PRINTED
C
C ** NOTE: NONE OF THESE CARDS ARE REPEATED, EVEN IF IPKT .GT. 1 **
C
C***********************************************************************
      COMMON/SP2/ S(286),EREN(286),Q20(286),V2(286),UV(286),EQP(286)
      COMMON/INT1/ IMOM,NYIC,J,KLARGE,KSMALL,KMAX
      COMMON/EXPECT1/ PUV(15,15),SQP(15),A1,A2,A3,CHSI
      COMMON/DECOUP/ APLB(15,15),C(15,15),D(15,15),E(15,15),
     $NORM(15,15),JZ(15,15),JZ2(15,15)
CC    COMMON/EIGVEC/ EIV(15,15,10),EIGVEL(225,39)
      COMMON/EIGVEC/ EIV(15,15,10)
      DIMENSION EIGVEL(225,39)
      REAL NORM,JZ,JZ2,NOUSE(15,15)
      REAL JZAV(150),JX2AV(150),JY2AV(150),JZ2AV(150),IDOTJ(150)
      REAL IX2AV(150),IY2AV(150),IZ2AV(150)
      DIMENSION RX2AV(150),RY2AV(150),RZ2AV(150),RDOTJ(150)
      DIMENSION HX(150),HY(150),HZ(150)
      REAL JPLR(161,15),JMNR(161,15),JZR(161,15)
      DIMENSION U1(161), V1(161)
cc      REAL A(25425),VALUE(39)
cc      DIMENSION SPCOEFF(15,11,11,6),SUM(15,39),OVLP(15,15,39)
cc      DIMENSION ROTEN(39,76),ROTEIG(225,39,76)
cc      DIMENSION AOR(76),IPOUT(20)
cc      DIMENSION AOR(76),IPOUT(20)
      REAL A(0:25425),VALUE(51)
      DIMENSION SPCOEFF(15,11,11,6),SUM(15,51),OVLP(15,15,51)
      DIMENSION ROTEN(51,101),ROTEIG(225,51,101)
      DIMENSION AOR(101),IPOUT(20)
      REAL SP(286)
      DIMENSION K3(15),V(225)
      REAL Q22(286),R2(286),RJZ(286)
      DIMENSION MOMENT(5)
      DIMENSION AE(150),AF(150),IAM(150),KAV(150)
      DIMENSION UPAIR(15), VPAIR(15)
      INTEGER LEVEL(15),LEVEL1(15),LEVEL2(15),KAX(15)
      DIMENSION IS(50)
      INTEGER Z,AA,Q
      DIMENSION NANTJ(20),NOUTJ(20)
      INTEGER NYDOM(150),KDOM(150)
      CHARACTER IPAR1*4,KATT*4,KV(286)*4
      CHARACTER FILE16*24,FILE17*24,FILE18*24
      EQUIVALENCE (A,EIV)
C
C     VARIABLES  "KV" ETC NOT USED IN PRESENT PROGRAM
      EQUIVALENCE (KV,Q22,R2,RJZ)
C
      DATA MATDIM/225/
      DATA PI/3.1415926536/
C     DATA IBLANK,ILINE/6H      ,6H------/
CV

      integer count,len,status,iostat
      character(len=132) inputfile
      count = command_argument_count ()
      if(count>0)then
            call get_command_argument (1, inputfile, len, iostat)
      open(unit=5,file=inputfile,status='old',form='formatted',
     1  iostat=iostat)
      endif
      open(unit=6,file='ASYRMO.out',status='unknown',form='formatted')



      CALL HELP
      READ(5,*)FILE16,FILE17,FILE18
      READ(5,*)IPKT,ISKIP
      READ(5,*)ISTRCH,IREC
      READ(5,*)VMI,INMIN,INMAX,IARCUT,AOO,STIFF
CV
      READ(5,*)Z,AA,IMIN,ISPIN,KMAX,E2PLUS,E2PLUR
      READ(5,*)GN0,GN1,IPAIR,CHSI,ETA
      READ(5,1984) IPAR1,NYIC,(LEVEL2(J),J=1,NYIC)
 1984 FORMAT(A1,I2,15I3)
C
C******************************  SHIFT OF GAMMA-BAND? *******
C     READ(5,*) IGAMSH,EGE0,EGO0
C
      READ(5,*) (NANTJ(I),I=1,20)
      READ(5,*) (NOUTJ(I),I=1,20)
      READ(5,*) (IPOUT(I),I=1,20)
C*** APPROPRIATE FOR NORD
CC    OPEN(16,FILE=FILE16,ACCESS='READ')
CC    OPEN(17,FILE=FILE17,ACCESS='READ')
CC    OPEN(18,FILE=FILE18,ACCESS='WRITE',STATUS='UNKNOWN')
C*** APPROPRIATE FOR VAX
      IF (VMI .EQ. 1)
     &   OPEN(16,FILE=FILE16,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(17,FILE=FILE17,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(18,FILE=FILE18,FORM='UNFORMATTED',STATUS='UNKNOWN')
CV
      IF (IPAR1 .EQ. '+') THEN
         INMIN=INMIN/2+1
         INMAX=INMAX/2+1
      ELSE
         INMIN=(INMIN+1)/2
         INMAX=(INMAX+1)/2
      ENDIF
CV
C
C*****LOOP OVER DIFFERENT DEFORMATIONS****
      DO 111 IP=1,IPKT
      READ(17) NU,(LEVEL(J),J=1,NU),IDIMA,(SP(J),J=1,IDIMA),
     $(((NOUSE(MYP,MY),MY=1,NU),MYP=1,NU),I=1,2)
CT    WRITE(1,*) 'NU,LEVEL(J),IDIMA,SP(J)',NU,(LEVEL(J),J=1,NU),
CT   *IDIMA,(SP(J),J=1,IDIMA)
      READ(17) KATT,JREC,NU,(LEVEL(J),J=1,NU),
     $((APLB(MYP,MY),MY=1,NU),MYP=1,NU),((C(MYP,MY),MY=1,NU),MYP=1,NU),
     $((D(MYP,MY),MY=1,NU),MYP=1,NU),((E(MYP,MY),MY=1,NU),MYP=1,NU),
     $((NORM(MYP,MY),MY=1,NU),MYP=1,NU),((JZ(MYP,MY),MY=1,NU),MYP=1,NU
     $),((JZ2(MYP,MY),MY=1,NU),MYP=1,NU),
     $(((NOUSE(MYP,MY),MY=1,NU),MYP=1,NU),I=1,21)
C     IF (ISTRCH .EQ. 1) GO TO 48
      READ (17) IRAP, (NRAP,NRAP,NRAP, I=1,IRAP)
      DO 2 I=1,NU
         READ(17) (RAP,J=1,IRAP)
    2 CONTINUE
      IF (JREC .EQ. 1) READ(17) ((JPLR(NY,MMY),MMY=1,NU),NY=1,IDIMA),
     &                          ((JMNR(NY,MMY),MMY=1,NU),NY=1,IDIMA),
     &                          ((JZR(NY,MMY),MMY=1,NU),NY=1,IDIMA)
   48 CONTINUE
      READ(17) EPSP,GAMMA,EPS4,EPS6,OMOM,Q,(S(K),K=1,Q),(KV(K),K=1,Q),
     $(Q20(K),K=1,Q),(Q22(K),K=1,Q),(R2(K),K=1,Q),(RJZ(K),K=1,Q)
      IF (IREC .EQ. 1 .AND. JREC .NE. 1) THEN
          write(6,40)
   40 FORMAT ('0 IREC=1 BUT JREC=0 => MATRIX ELEMENTS REQUIRED FOR ',
     &     'TWO BODY TREATMENT OF RECOIL TERMS NOT SUPPLIED BY GAMMA')
          STOP
      ENDIF
C
C****MANIPULATIONS TO GET THE S.P. MATRIX ELEMENTS IN THEIR PROPER
C    PLACES  (IF LEVEL2(J) .NE. LEVEL(J) OF GAMMA4)
C
C    15 DIFFERENT PARTICLE STATES CAN BE INCLUDED (NYIC.LE.15)
C
      DO 36 J=1,NYIC
         LEVEL1(J) = LEVEL2(J)
   36 CONTINUE
      DO 31 J=1,NYIC
      DO 32 I=1,NU
         IK=I
         IF (LEVEL1(J) .EQ. LEVEL(I)) GO TO 33
   32 CONTINUE
      write(6,41)LEVEL1(J),IPAR1
   41 FORMAT ('0 NO DECOUPLING PARAMETERS CALCULATED FOR LEVEL:',I3,A1)
      STOP
   33 LEVEL1(J) = IK
   31 CONTINUE
C
C*** ALSO REARRANGE JPLR ETC, BUT ONLY IF NEEDED
C
      IF (IREC .EQ. 1 .AND. VMI .NE. 1) THEN
         DO 144 I=1,NYIC
            IF (LEVEL2(I) .NE. LEVEL(I)) THEN
               II = LEVEL1(I)
               DO 145 J=1,IDIMA
                  JPLR(J,I) = JPLR(J,II)
                  JMNR(J,I) = JMNR(J,II)
                  JZR (J,I) = JZR (J,II)
  145          CONTINUE
            ENDIF
  144    CONTINUE
      ENDIF
C
      DO 34 I=1,NYIC
         MYP = LEVEL1(I)
         DO 35 J=1,NYIC
            MY = LEVEL1(J)
            APLB(I,J) = APLB(MYP,MY)
            C   (I,J) = C   (MYP,MY)
            D   (I,J) = D   (MYP,MY)
            E   (I,J) = E   (MYP,MY)
            NORM(I,J) = NORM(MYP,MY)
            JZ(I,J) = JZ(MYP,MY)
            JZ2(I,J) = JZ2(MYP,MY)
   35    CONTINUE
         LEVEL(I) = LEVEL(MYP)
   34 CONTINUE
C****END MANIPULATIONS
C
CV
CV    READ  S.P EIGENVECTORS
CV
      IF (VMI .EQ. 1) THEN
      READ(16) IDIMM
      DO 820 JJJ=1,NU
      IREAD=0
      DO 823 ILEV=1,NYIC
      IF(LEVEL1(ILEV) .EQ. JJJ) THEN
      READ(16) (IN,ILL,IJJ,IOME,SPCOEFF(ILEV,(IJJ+1)/2,
     $(ABS(IOME)+1)/2,IN),I=1,IDIMM)
      IREAD=1
      ENDIF
  823 CONTINUE
      IF (IREAD .EQ. 0) READ(16) (IN,ILL,IJJ,IOME,DUMMY,I=1,IDIMM)
  820 CONTINUE
      ENDIF
CV
CV
      IF (IP.LE.ISKIP) GO TO 111
C
      write(6,100) Z,AA
  100 FORMAT(1H1,'SPECTRUM OF (Z=',I3,' A=',I3,') IN PARTICLE-ROTOR ',
     *' MODEL WITH MODIFIED OSCILLATOR POTENTIAL',//)
      IF (ISTRCH .EQ. 1) THEN
         write(6,2000) KATT,EPSP,GAMMA,EPS4,EPS6,(LEVEL(J),J=1,NYIC)
      ELSEIF (ISTRCH .EQ. 0) THEN
         write(6,2003) KATT,EPSP,GAMMA,EPS4,EPS6,(LEVEL(J),J=1,NYIC)
      ENDIF
 2000 FORMAT(1X,'PARITY=',A1,3X,'EPS=',F6.3,3X,'GAMMA=',F5.1,3X,
     *'EPS4=',F5.3,'   EPS6=',F5.3,3X,'LEVELS:',15I4)
 2003 FORMAT(1X,'PARITY=',A1,3X,'DELTA=',F6.3,3X,'GAMMAD=',F5.1,3X,
     *'DELTA4=',F5.3,'   DELTA6=',F5.3,3X,'LEVELS:',15I4)
CV
      IF (VMI .EQ. 1) THEN
      write(6,2001) 'PAIRING FACTORS: GN0=',GN0,'  GN1=',GN1,'  IPAIR=',
     $IPAIR,' ***  A00(INPUT)=',AOO,'  CORE SOFTNESS(INPUT)=',4*AOO**3
     $/STIFF
      ELSE
      write(6,2001) 'PAIRING FACTORS: GN0=',GN0,'  GN1=',GN1,'  IPAIR=',
     $IPAIR,' ***  E2PLUS(INPUT)=',E2PLUS,'  E2PLUR(INPUT)=',E2PLUR
      ENDIF
CV
      write(6,2002) 'CORIOLIS ATTENUATION PARAMETERS:  CHSI=',CHSI,
     $'  ETA=',ETA
      IF (ETA .NE. 1) write(6,2004)
 2001 FORMAT(/,' ',2(A,F6.2),A,I3,2(A,F8.3))
 2002 FORMAT(/,' ',A,F5.3,A,F6.2)
 2004 FORMAT(/5X,'FOR ETA .NE. 1, THE ANGULAR MOMENTUM EXPECTATION',
     &           ' VALUES OF THE ENERGY EIGENSTATES ARE CALCULATED',
     &           ' ONLY APPROXIMATELY')
C
C*** DETERMINE ODD NUMBER OF PARTICLES, Z OR N
C
      IF ((Z-1)/2 *2 .EQ. Z-1) THEN
         NODD = Z
         EFAC=(1.-(AA-2.*Z)/FLOAT(AA)/3.)
      ELSE
         NODD = AA-Z
         EFAC=(1.+(AA-2.*Z)/FLOAT(AA)/3.)
      ENDIF
      EFAC=41.0/(FLOAT(AA)**(1./3.))*EFAC*OMOM
C     EFAC=41.0/(FLOAT(AA)**(1./3.))*(1.-(AA-2.*Z)/FLOAT(AA)/3.)*OMOM
      IGAM=0
      IF (GAMMA .EQ. 0. .AND. E2PLUR .LE. 0) IGAM=1
      IF (GAMMA .EQ. 0.) GAMMA=0.1
      IF (GAMMA .EQ. 60.) GAMMA=59.9
      IF (E2PLUS .GT. 0. .AND. E2PLUR .GT. 0.) THEN
C***     ANGLE GAMMA OF CORE
         GR=ASIN(SQRT(1.125*(1.-((E2PLUR/E2PLUS-1.)/
     *   (E2PLUR/E2PLUS+1.))**2)))
         GR=GR/3.
         IF (GAMMA .GT. 30. .AND. GAMMA .LT. 60.) GR = PI/3.-GR
         GRGRAD = GR*180./PI
         write(6,11) GRGRAD
   11    FORMAT(/,'   ********  CORE INERTIAS CORRESPOND TO GAMMA=',
     $   F6.1,'  ********')
      ELSE
         GR=GAMMA*PI/180.
      ENDIF
      XSIN=SIN(3.*GR)**2
C***  MOMENT OF INERTIA PARAMETER B
      H2OB=6.*EPSP**2*E2PLUS*XSIN/(4.5*(3.-SQRT(9.-8.*XSIN)))
      IF (E2PLUS .EQ. 0.) H2OB = 1225./FLOAT(AA)**(7./3.)
      IF (E2PLUS .LT. 0.) H2OB = -E2PLUS/FLOAT(AA*AA)
      BHARM = H2OB*FLOAT(AA)**(7./3.)
      BHARM1 = H2OB*FLOAT(AA*AA)
      IF (E2PLUS .LE. 0.)
     *E2PLUT = 0.75*H2OB*(3.-SQRT(9.-8.*XSIN))/EPSP**2/XSIN
      IF (E2PLUS .GT. 0.) E2PLUT = E2PLUS
CV
      IF (VMI .EQ. 1 ) THEN
      H2OB = 6.*EPSP**2
      CALL INERTIA(2,AOO,STIFF,AO2)
      E2PLUT = 4.5*AO2*(3.-SQRT(9.-8.*XSIN))/XSIN
     $+STIFF/8*(1./AOO-1./AO2)**2
      ENDIF
CV
C
C*****MOMENT OF INERTIAS - GAMMA=0 OR =60 LEADS TO INFINITIES *****
C
      IF (ABS(GAMMA) .LT. 1.0) THEN
         A1=1./3.*H2OB/2./EPSP**2
         A2=A1
         A3=1000.*A2
      ELSEIF (ABS(GAMMA-60.) .LT. 1.0) THEN
         A2=1./3.*H2OB/2./EPSP**2
         A3=A2
         A1=1000.*A2
      ELSE
         A1=1./4./SIN(GR+2.*PI/3.)**2*H2OB/2./EPSP**2
         A2=1./4./SIN(GR-2.*PI/3.)**2*H2OB/2./EPSP**2
         A3=1./4./SIN(GR)**2*H2OB/2./EPSP**2
      ENDIF
CV
      IF (VMI .EQ. 1) THEN
      write(6,151) A1*AOO,A2*AOO,A3*AOO
      ELSE
      write(6, 150) A1,A2,A3
      ENDIF
  151 FORMAT(/,' MOMENT OF INERTIA PARAMETERS:   A1=',
     $F7.3,'*A0R/A00 MEV   A2=',F7.3,'*A0R/A00 MEV
     $  A3=',F9.5,'*A0R/A00 MEV',/)
CV
  150 FORMAT(/,' MOMENT OF INERTIA PARAMETERS:   A1=',
     $F7.3,' MEV   A2=',F7.3,' MEV   A3=',F9.5,' MEV',/)
      DO 10 K=1,Q
         S(K) = S(K)*EFAC
   10 CONTINUE
      AMASS=FLOAT(AA)
C****INITIAL VALUES FOR PAIRING CALCULATIONS, ETC.
C
      L=NODD/2+1
      EF = S(L)
      IEMIN=NODD/2-INT(SQRT(FLOAT(IPAIR*NODD)))
      IF(IEMIN.LT.1) IEMIN=1
      IEMAX=NODD/2+INT(SQRT(FLOAT(IPAIR*NODD)))
      IF (GN0 .GT. 0.01 .AND. IPAIR .NE. 99 ) THEN
         ENE = NODD-2*(IEMIN-1)
      IF (NODD .EQ. Z) THEN
         GM=GN0+GN1*(AA-2.*Z)/AMASS
      ELSE
         GM=GN0-GN1*(AA-2.*Z)/AMASS
      ENDIF
C        GM=GN0+GN1*(AA-2.*Z)/AMASS
         DEL = 12./SQRT(AMASS)
         DEL2=0.0
         TOLER = 0.0001
         GM2=0.
C
C******* SOLUTION OF PAIRING HAMILTONIAN;
C        PARAMETERS USED: IEMIN, IEMAX, AMASS, ENE, GM, DEL, EF, TOLER
C        OTHER PARAMETERS REDUNDANT (OR 0) FOR SOLUTION OF MONOPOLE
C        PAIRING BCS-HAMILTONIAN
C
         CALL PAIRQ2(IEMIN,IEMAX,0,0,0,AMASS,TOLER,
     *   1.,ENE,GM,GM2,ODDE,DEL,DEL2,EF,G0,WE,WWE)
C
      ELSEIF (GN0 .LE. 0.001 .OR. IPAIR .EQ. 99 ) THEN
         DEL=GN1/SQRT(AMASS)
         write(6, 27) 'NO BCS-EQUATION SOLVED'
   27    FORMAT(/,' ',A)
      ENDIF
      IF (IPAIR.EQ.99) EF=GN0
C
      EFHW0 = EF/EFAC
      write(6, 22) DEL,EF,EFAC,EFHW0,E2PLUT,BHARM,BHARM1
   22 FORMAT(/,1X,'DELTA=',F5.3,' MEV',2X,'LAMBDA=',F7.4,' MEV',2X,
     &   ' HW0=',F7.4,' MEV   LAMBDA/HW0=',F7.4,//,' E2ENERGY=',F5.3,
     &' MEV  HBAR**2/B=',F6.0,9H/A**(7/3),1H=,F6.1,5H/A**2,' MEV')

      IF (IREC .EQ. 1 .AND. VMI .NE. 1) write(6,123)
  123 FORMAT(/,1X,'IREC=1  =>  RECOIL TERMS j+j-, j+j+, jzjz, ETC',
     &' TREATED AS TWO-BODY OPERATORS WITH APPROPRIATE PAIRING FACTORS')
      IF (IREC .NE. 1 .AND. VMI .NE. 1) write(6, 124)
  124 FORMAT(/,1X,'IREC=0  =>  RECOIL TERMS j+j-, j+j+, jzjz, ETC',
     &' TREATED AS ONE-BODY OPERATORS WITH UU+VV PAIRING FACTORS')
CC    DO 26 MMY=1,NYIC
CC       MY=LEVEL(MMY)
CC       SQP(MMY)=SQRT((SP(MY)*EFAC-EF)**2+DEL**2)
CC       UPAIR(MMY)=SQRT(0.5*(1.+(SP(MY)*EFAC-EF)/SQP(MMY)))
CC       VPAIR(MMY)=SQRT(0.5*(1.-(SP(MY)*EFAC-EF)/SQP(MMY)))
CC 26 CONTINUE
      SPEMIN = S(IEMIN)/EFAC
      SPEMAX = S(IEMAX)/EFAC
      DO 126 I=1,IDIMA
         IF (SP(I) .LT. SPEMIN) THEN
            V1(I) = 1.0
            U1(I) = 0.0
         ELSEIF (SP(I) .GT. SPEMAX) THEN
            V1(I) = 0.0
            U1(I) = 1.0
         ELSE
            EQP(I)=SQRT((SP(I)*EFAC-EF)**2+DEL**2)
            U1(I)=SQRT(0.5*(1.+(SP(I)*EFAC-EF)/EQP(I)))
            V1(I)=SQRT(0.5*(1.-(SP(I)*EFAC-EF)/EQP(I)))
         ENDIF
  126 CONTINUE
      DO 26 MMY=1,NYIC
         MY=LEVEL(MMY)
         SQP(MMY)=EQP(MY)
         UPAIR(MMY)=U1(MY)
         VPAIR(MMY)=V1(MY)
   26 CONTINUE
      DO 20 MMY=1,NYIC
      DO 20 NNY=1,NYIC
         PUV(MMY,NNY) =
     $   (UPAIR(MMY)*UPAIR(NNY) + VPAIR(MMY)*VPAIR(NNY))**ETA
   20 CONTINUE
C
C*** INCORPORATE UU'+VV' FACTORS DIRECTLY INTO APLB, ETC, TO MAKE
C*** COMPARABLE TO D AND E WHEN IREC=1; NOTE THAT PUV(I,I)=1 SO
C*** THAT THE DIAGONAL ELEMENTS ARE NOT CHANGED
C
      DO 120 I=1,NYIC
      DO 120 J=1,NYIC
         APLB(I,J) = APLB(I,J)*PUV(I,J)
         C   (I,J) = C   (I,J)*PUV(I,J)
         D   (I,J) = D   (I,J)*PUV(I,J)
         E   (I,J) = E   (I,J)*PUV(I,J)
         JZ(I,J) = JZ(I,J)*PUV(I,J)
         JZ2(I,J) = JZ2(I,J)*PUV(I,J)
  120 CONTINUE
C
C*** CONSTRUCT TWO BODY MATRIX ELEMENTS OF JZ2, J+J+, J-J-, ETC
C*** BY SUMMING OVER INTERMEDIATE STATES
C
      IF (IREC .EQ. 1 .AND. VMI .NE. 1) THEN
         DO 130 II=1,NYIC
         DO 130 IF=1,NYIC
            D(II,IF) = 0.0
            E(II,IF) = 0.0
            JZ2(II,IF) = 0.0
  130    CONTINUE
         DO 135 IK=1,IDIMA
         DO 134 II=1,NYIC
            IIY=LEVEL(II)
         DO 133 IF=II,NYIC
            IFY=LEVEL(IF)
            PFAC=((U1(IK)*U1(IIY) + V1(IK)*V1(IIY))*
     &            (U1(IK)*U1(IFY) + V1(IK)*V1(IFY))-
     &            (U1(IK)*V1(IIY) - V1(IK)*U1(IIY))*
     &            (U1(IK)*V1(IFY) - V1(IK)*U1(IFY)))
            D  (IF,II) = D(IF,II) + (JPLR(IK,II)*JPLR(IK,IF) +
     &                   JMNR(IK,II)*JMNR(IK,IF))*PFAC/2.
            E  (IF,II) = E(IF,II) + (JPLR(IK,II)*JMNR(IK,IF) +
     &                   JMNR(IK,II)*JPLR(IK,IF))*PFAC
            JZ2(IF,II) = JZ2(IF,II) + JZR(IK,II)*JZR(IK,IF)*PFAC
  133    CONTINUE
  134    CONTINUE
  135    CONTINUE
         DO 136 II=1,NYIC
         DO 136 IF=II,NYIC
            D(II,IF) = D(IF,II)
            E(II,IF) = E(IF,II)
            JZ2(II,IF) = JZ2(IF,II)
  136    CONTINUE
      ENDIF
C
      DO 29 K=1,IEMAX
         V2(K) = (1.-(S(K)-EF)/SQRT((S(K)-EF)**2+DEL**2))/2.
         IF (K .LT. IEMIN) V2(K)=1.
   29 CONTINUE

      write(6, 23) ETA
   23 FORMAT(///,' S. P. LEVEL   EQP       U     V',TR8,
     $'(UU + VV)**',F4.2,/)
   21 CONTINUE
      DO 25 MY=1,NYIC
         write(6, 24) LEVEL(MY),SQP(MY),UPAIR(MY),VPAIR(MY),
     $   (PUV(MY,NY),NY=1,NYIC)
   24    FORMAT(I7,6X,F6.3,3X,2F6.3,5X,15F6.3)
   25 CONTINUE
CV
CV
CV    CALCULATION OF CORE STATES
CV
CV
      IF (VMI .EQ. 1 .OR. IPOUT(1) .NE. 0) THEN
cc      DO 801 IR=0,21+ISPIN,2
      DO 801 IR=0,17+ISPIN,2
      CALL INERTIA(IR,AOO,STIFF,AOR(IR+1))
      IF (VMI .NE. 1) AOR(IR+1)=1
      IF (VMI .NE. 1) AOO=1
      DO 802 III=0,(IR+2)*(IR+4)/8
      A(III)=0
  802 CONTINUE
      DO 803 KUP=0,IR,2
      DO 803 ICOUNT=0,1
      ANORM=0.25
      IF (KUP .EQ. 0) ANORM=0.35355339
      INDEX=KUP*(IR+3-KUP/2)/4+ICOUNT+1
      IF (ICOUNT .EQ. 0) A(INDEX)=-(AOR(IR+1)*(A1+A2)/2.*(IR*(IR+1)
     $-KUP*KUP)+AOR(IR+1)*A3*KUP*KUP+STIFF/8.*(1./AOO-1./AOR(IR+1))**2)
      IF (ICOUNT .EQ. 1 .AND. KUP .LT. IR)
     $A(INDEX)=-AOR(IR+1)*ANORM*(A1-A2)*(((IR*(IR+1)-KUP*
     $(KUP+1))*(IR*(IR+1)-(KUP+1)*(KUP+2)))**.5)
  803 CONTINUE
      NN=IR/2+1
      IDIMR=IR/2+1
      CALL BIGMAX(A,VALUE,NN,IDIMR,IDIMR,EIGVEL,MATDIM)
      DO 804 L=1,IDIMR
      VALUE(L) = -VALUE(L)
      IF (L .GT. (IARCUT+1)/2) VALUE(L) = 10.*VALUE(L)
      ROTEN(L,IR+1)=VALUE(L)
      DO 804 LL=1,NN
      ROTEIG(LL,L,IR+1)=EIGVEL(LL,L)
  804 CONTINUE
  801 CONTINUE
      DO 805 IR=3,21+ISPIN,2
      CALL INERTIA(IR,AOO,STIFF,AOR(IR+1))
      IF (VMI .NE. 1) AOR(IR+1)=1
      IF (VMI .NE. 1) AOO=1
      DO 806 III=0,(IR-1)*(IR+1)/8
      A(III)=0
  806 CONTINUE
      DO 807 KUP=2,IR,2
      DO 807 ICOUNT=0,1
      ANORM=0.250
      INDEX=(KUP-2)*(2*IR+2-KUP)/8+ICOUNT+1
      IF (ICOUNT .EQ. 0) A(INDEX)=-(AOR(IR+1)*(A1+A2)/2.*(IR*(IR+1)
     $-KUP*KUP)+AOR(IR+1)*A3*KUP*KUP+STIFF/8.*(1./AOO-1./AOR(IR+1))**2)
      IF (ICOUNT .EQ. 1 .AND. KUP .LT. IR-1)
     $A(INDEX)=-AOR(IR+1)*ANORM*(A1-A2)*(((IR*(IR+1)-KUP*
     $(KUP+1))*(IR*(IR+1)-(KUP+1)*(KUP+2)))**.5)
  807 CONTINUE
      NN=(IR-1)/2
      IDIMR=NN
      CALL BIGMAX(A,VALUE,NN,IDIMR,IDIMR,EIGVEL,MATDIM)
      DO 808 L=1,IDIMR
      VALUE(L) = -VALUE(L)
      IF (L .GT. (IARCUT)/2) VALUE(L) = 10.*VALUE(L)
      ROTEN(L,IR+1)=VALUE(L)
      DO 808 LL=1,NN
      ROTEIG(LL+1,L,IR+1)=EIGVEL(LL,L)
  808 CONTINUE
  805 CONTINUE
CV
CV
CV    SHIFT OF GAMMA-BAND?
CV
CV
C     IF (IGAMSH .EQ. 0) GO TO 903
C     DO 901 IR=0,21+ISPIN,2
C        ROTEN(2,IR+1) = ROTEN(1,IR+1) + EGE0
C 901 CONTINUE
C     DO 902 IR=3,21+ISPIN,2
C        ROTEN(1,IR+1) = (ROTEN(1,IR)+ROTEN(1,IR+2))/2. + EGO0
C 902 CONTINUE
C 903 CONTINUE
CV
CV
CV    PRINT OUT CORE STATES
CV
CV
      IF (IPOUT(1 ) .EQ. 0) GOTO 843
      DO 840 IR=2,20
CC    IF (IPOUT(IR) .EQ. 0) CONTINUE
      IF (IPOUT(IR) .EQ. 0) GO TO 840
      IDIM=(IR-1)/2
      IF (IR/2*2 .EQ. IR) IDIM=IR/2+1
      write(6, 841) AOR(IR+1)/AOO,(IK,IK=0,IR,2)
      IPOUT(IR)=MIN0(IPOUT(IR),IDIM)
!      DO 840 II=1,IPOUT(IR)
      DO  II=1,IPOUT(IR)
      write(6, 842) IR,II,ROTEN(II,IR+1),
     $(ROTEIG(III,II,IR+1),III=1,IR/2+1)
      ENDDO
  840 CONTINUE
  841 FORMAT(/,2X,'R CORE   ENERGY   A0R/A00=',F5.3,
     $'    R3=',I2,2H  ,10(I4,2H  ),/,X)
  842 FORMAT(2X,I2,'(',I2,')',2X,F7.3,19X,11F6.3)
  843 CONTINUE
      ENDIF
CV
CV
      IF (IGAM .EQ. 1) THEN
         write(6, 42)
   42 FORMAT(/,10X,'AXIAL SYMMETRY => K RESTRICTED TO OMEGA',//,21X,'K')
         DO 44 J=1,NYIC
            IF (JZ(J,J) .GE. 0) THEN
               KAX(J)=IFIX(2*JZ(J,J)+0.1)
            ELSE
               KAX(J)=IFIX(2*JZ(J,J)-0.1)
            ENDIF
            write(6, 43) LEVEL2(J),JZ(J,J)
   43 FORMAT(5X,'S.P. ORB#:',I3,F6.2)
   44    CONTINUE
      ENDIF
C     HIGHEST SPIN= ISPIN/2       (ISPIN .LE. 129)
C     TEST KMAX
      IF ( (KMAX-1)/4*4 .NE. KMAX-1) THEN
         KMAX = KMAX-2
         write(6, 13) KMAX
      ENDIF
   13 FORMAT(//10X,'WARNING! KMAX RESET TO ',I3)
      write(6, 28)
   28 FORMAT(1H )
      WRITE(18) IMIN,ISPIN,KMAX,NYIC,
     &(LEVEL(I),UPAIR(I),VPAIR(I),I=1,NYIC),IEMAX,(V2(K),K=1,IEMAX)

!      IRK=MSLEFT(XX)
!      call cpu_time(begin_time)

      JKJ = 0
C
C     LOOP OF SPINS, 1/2, 3/2, 5/2, ETC.
C
      DO 4000 IMNT=IMIN,ISPIN,2
         IMOM=IMNT
         INUM=MIN0(IMOM,KMAX)
         IDIM=(INUM+1)/2
         IF (IGAM .EQ. 1) IDIM=1
         IDI=NYIC*IDIM
         II = (IMNT-IMIN+2)/2
         IDIME = NANTJ(II)
         IDIME = MIN0(IDI,IDIME)
         IDIT = IDI*(IDI+1)/2
         DO 711 JJ=1,IDIT
            A(JJ) =0.
  711    CONTINUE
C***********  IMOM REPLACED WITH KMAX HERE
         N=(IMOM+1)/2
         IDIM1=NYIC*IDIM
C        KMAX=IMOM
C        IF((IMOM-1)/4*4.NE.IMOM-1) KMAX=IMOM-2
         LE = MIN0(5,IDIM1)
         DO 15 IJK = 1,LE
            MOMENT(IJK)=IMOM
   15    CONTINUE
         KKMAX=IMOM
         IF((IMOM-1)/4*4.NE.IMOM-1) KKMAX=IMOM-2
         KKMAX=MIN0(KMAX,KKMAX)
CV
CV
CV    CALCULATION OF OVERLAPS
CV
CV
      IF (VMI .EQ. 1) THEN
      ILIMIT1=IDIM
      ILIMIT2=NYIC
      IRSTEP=2
      IF (IGAM .EQ. 1) THEN
      ILIMIT1=NYIC
      ILIMIT2=1
      IRSTEP=4
      ENDIF
      DO 831 IN=INMIN,INMAX
      IJMAX=4*IN-1
      IF (IPAR1 .EQ. '+') IJMAX=4*IN-3
      DO 831 IJ=1,IJMAX,2
      IRMIN=ABS(IMNT-IJ)
      IF (IGAM .EQ. 1 .AND. IRMIN/4*4 .NE. IRMIN) IRMIN=IRMIN+2
      IRMAX=IMNT+IJ
      DO 831 IR=IRMIN,IRMAX,IRSTEP
      SPINFAC=(((FLOAT(IR)+1)/(FLOAT(IMNT)+1))**.5)/2.0
      DO 810 III=1,ILIMIT1
      IKK=KKMAX-4*(III-1)
      IF (IGAM .EQ. 1) IKK=KAX(III)
      IAMAX=(IR/2-1)/2
      IF( IR/4*4 .EQ. IR) IAMAX=IR/4+1
      IF (IGAM .EQ .1) IAMAX=1
      DO 814 ISP=1,NYIC
      DO 814 IA=1,IAMAX
      SUM(ISP,IA)=0
  814 CONTINUE
      NNNMAX=IR/2
      IF (IGAM .EQ .1) NNNMAX=0
      DO 811 NNN=0,NNNMAX,2
      DELTA=1
      IF( NNN .EQ. 0) DELTA=0.707106781
      IF (ABS(IKK-2*NNN) .LE. IJ) THEN
      CGSUM=DELTA*(CLEBI(IR,IJ,IMNT,2*NNN,IKK-2*NNN,IKK)+
     $CLEBI(IR,IJ,IMNT,-2*NNN,2*NNN-IKK,-IKK)*(-1)**((IR+IMNT-IJ)/2))
      DO 812 IISP=1,ILIMIT2
      ISP=IISP
      IF (IGAM .EQ. 1) ISP = III
      SPC=SPCOEFF(ISP,(IJ+1)/2,(ABS(IKK-2*NNN)+1)/2,IN)
      DO 812 IA=1,IAMAX
      SUM(ISP,IA)=SUM(ISP,IA)+ROTEIG(NNN/2+1,IA,IR/2+1)*SPC*CGSUM
  812 CONTINUE
      ENDIF
      IF (ABS(IKK+2*NNN) .LE. IJ) THEN
      CGSUM=DELTA*(CLEBI(IR,IJ,IMNT,-2*NNN,IKK+2*NNN,IKK)*(-1)**(IR/2)+
     $CLEBI(IR,IJ,IMNT,2*NNN,-2*NNN-IKK,-IKK)*(-1)**((IMNT-IJ)/2))
      DO 813 IISP=1,ILIMIT2
      ISP=IISP
      IF (IGAM .EQ. 1) ISP = III
      SPC=SPCOEFF(ISP,(IJ+1)/2,(ABS(IKK+2*NNN)+1)/2,IN)
      DO 813 IA=1,IAMAX
      SUM(ISP,IA)=SUM(ISP,IA)+ROTEIG(NNN/2+1,IA,IR/2+1)*SPC*CGSUM
  813 CONTINUE
      ENDIF
  811 CONTINUE
      DO 810 IISP=1,ILIMIT2
      ISP=IISP
      IF (IGAM .EQ. 1) ISP = III
      DO 810 IA=1,IAMAX
      OVLP((ABS(IKK)+1)/2,ISP,IA)=SUM(ISP,IA)*SPINFAC
  810 CONTINUE
CV
CV    CONSTRUCT PARTIAL MATRIX ELEMENTS
CV
      DO 831 I1=1,NYIC
      DO 831 I2=1,I1
      I3=IDIM*(I1-1)
      I4=IDIM*(I2-1)
      DO 831 III=1,IDIM
      IF (IGAM .EQ. 1) THEN
      KP=KAX(I1)
      ELSE
      KP=KKMAX-4*(III-1)
      ENDIF
      DO 831 JJJ=1,IDIM
      IF (IGAM .EQ. 1) THEN
      K=KAX(I2)
      ELSE
      K=KKMAX-4*(JJJ-1)
      ENDIF
      K3(JJJ)=K
      I=III+I3
      J=JJJ+I4
      INDEX = (2*IDIM1-J)*(J-1)/2 + I
      IF(J.GT.I) GO TO 831
      DO IA=1,IAMAX
      A(INDEX)=A(INDEX)+OVLP((ABS(KP)+1)/2,I1,IA)
     $*ROTEN(IA,IR/2+1)*OVLP((ABS(K )+1)/2,I2,IA)
      ENDDO
  831 CONTINUE
CV
CV    ADD PAIRING AND S.P. FACTORS
CV
      DO 832 I1=1,NYIC
      DO 832 I2=1,I1
      I3=IDIM*(I1-1)
      I4=IDIM*(I2-1)
      DO 832 III=1,IDIM
      IF (IGAM .EQ. 1) THEN
      KP=KAX(I1)
      ELSE
      KP=KKMAX-4*(III-1)
      ENDIF
      DO 832 JJJ=1,IDIM
      IF (IGAM .EQ. 1) THEN
      K=KAX(I2)
      ELSE
      K=KKMAX-4*(JJJ-1)
      ENDIF
      K3(JJJ)=K
      I=III+I3
      J=JJJ+I4
      INDEX = (2*IDIM1-J)*(J-1)/2 + I
      IF(J.GT.I) GO TO 832
C
C     REMOVE SPURIOUS STATES
C
            KTEST=MAX0(IABS(K),IABS(KP))
            IF (KTEST .GT. IMOM .AND. J .EQ. I .AND .I1 .EQ. I2) THEN
            A(INDEX)=50*(A1+A2)*(IMOM*(IMOM+2)+K*K)+SQP(I1)
            A(INDEX)=-A(INDEX)
            ENDIF
            IF (KTEST .GT. IMOM) GO TO 832
C
      A(INDEX) = A(INDEX)*PUV(I1,I2)
      IF (I1 .NE. I2) A(INDEX) = A(INDEX)*CHSI
      IF (I1 .EQ. I2 .AND. I .EQ. J) A(INDEX) = A(INDEX) + SQP(I1)
      A(INDEX) = - A(INDEX)
  832 CONTINUE
CV
CV
      ELSE
C
C**** SETTING UP OF MATRIX TO BE DIAGONALIZED *************
C
         DO 1 I1=1,NYIC
         DO 1 I2=1,I1
            IF(I1.EQ.I2) DELTA=1.
            IF(I1.NE.I2) DELTA=0.0
            I3=IDIM*(I1-1)
            I4=IDIM*(I2-1)
         DO 1 III=1,IDIM
            IF (IGAM .EQ. 1) THEN
            KP=KAX(I1)
            ELSE
            KP=KKMAX-4*(III-1)
            ENDIF
         DO 1 JJJ=1,IDIM
            IF (IGAM .EQ. 1) THEN
            K=KAX(I2)
            ELSE
            K=KKMAX-4*(JJJ-1)
            ENDIF
            K3(JJJ)=K
            RK=FLOAT(K)/2.
            I=III+I3
            J=JJJ+I4
            INDEX = (2*IDIM1-J)*(J-1)/2 + I
            IF(J.GT.I) GO TO 1
C
C     REMOVE SPURIOUS STATES
C
            KTEST=MAX0(IABS(K),IABS(KP))
            IF (KTEST .GT. IMOM .AND. J .EQ. I .AND .I1 .EQ. I2) THEN
            A(INDEX)=50*(A1+A2)*(IMOM*(IMOM+2)+K*K)+SQP(I1)
            A(INDEX)=-A(INDEX)
            ENDIF
            IF (KTEST .GT. IMOM) GO TO 1
C
            IF(K.NE.1) GO TO 3
            IF (KP .EQ. 1) A(INDEX) = (A1+A2)/2.*((IMOM*(IMOM+2.)/4.-
     $      1./4.)*DELTA+APLB(I1,I2)*(-1)**N*N)+
     $      (A1-A2)/2.*C(I1,I2)*(-1)**N*N+
     $      A3*(RK**2*NORM(I1,I2)-2.*RK*JZ(I1,I2)+JZ2(I1,I2))
     $      +(A1+A2)/2.*D(I1,I2)+(A1-A2)/4.*E(I1,I2)
            IF (KP .EQ. -3) A(INDEX) = ((A1-A2)*((IMOM+1.)/4.*DELTA+
     $      APLB(I1,I2)*(-1)**N)+(A1+A2)*C(I1,I2)*(-1)**N)*
     $      SQRT((IMOM-1.)*(IMOM+3.)/4.)/2.
            GO TO 4
    3       CONTINUE
            IF (KP .EQ. K-4) A(INDEX) = (A1-A2)/4.*DELTA*
     $      SQRT((IMOM+K)*(IMOM+K-2.)*(IMOM-K+2)*(IMOM-K+4.)/16.)
            IF (KP .EQ. -K+2)
     *      A(INDEX) = ((A1+A2)*APLB(I1,I2)+(A1-A2)*C(I1,I2))/2.*
     $      (-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K-2.)/4.)
            IF (KP .EQ. -K-2)
     *      A(INDEX) = ((A1-A2)*APLB(I1,I2)+(A1+A2)*C(I1,I2))/2.*
     $      (-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K+2.)/4.)
            IF (KP .EQ. K)
     *      A(INDEX) = (A1+A2)/2.*(IMOM*(IMOM+2.)/4.-K*K/4.)*DELTA+
     $      A3*(RK**2*NORM(I1,I2)-2.*RK*JZ(I1,I2)+JZ2(I1,I2))
     $      +(A1+A2)/2.*D(I1,I2)+(A1-A2)/4.*E(I1,I2)
    4       CONTINUE
CC          A(INDEX) = A(INDEX)*PUV(I1,I2)
            IF (I1 .NE. I2) A(INDEX) = A(INDEX)*CHSI
            IF (I1 .EQ. I2 .AND. I .EQ. J) A(INDEX) = A(INDEX) + SQP(I1)
            A(INDEX) = - A(INDEX)
    1    CONTINUE
      ENDIF
         IF (IGAM .EQ. 1) THEN
         write(6, 19)
         ELSE
         write(6, 17) (K3(L),L=1,IDIM)
         ENDIF
   17    FORMAT(/,12X,' ENERGY     K=',I2,2H/2,14(I4,2H/2))
   19    FORMAT(/,12X,' ENERGY    K= OMEGA')
C
CN   print matrix
C      IPDIM = IDIM1*(IDIM1+1)/2
C      IF (MOD((IMNT+1),24) .EQ. 0) THEN
C         WRITE(6,5219) (A(IPNEW),IPNEW=1,IPDIM)
C      ENDIF
C5219  FORMAT(12F10.4)
CN   END PRINT MATRIX
C
C**** DIAGONALIZATION  *******
C
!      WRITE(6,*) '   IDIM1,IDIME,MATDIM', IDIM1,IDIME,MATDIM

         CALL BIGMAX(A,VALUE,IDIM1,IDIME,IDIME,EIGVEL,MATDIM)
         DO 712  JJ = 1,IDIME
            VALUE(JJ) = -VALUE(JJ)
  712    CONTINUE
C
C**** ADMINISTRATION AND PRINTING *****************
C
      DO 85 L=1,15
      DO 85 NY=1,15
      DO 85 J=1,10
         EIV(L,NY,J) = 0.
   85 CONTINUE
C **** SPIN-PROJ K STORED IN DECREASING VALUES ****
         KLARGE = (KMAX+2-KKMAX)/4+1
CC       KSMALL = KLARGE + IDIM -1
         KSMALL = KLARGE + (INUM -1)/2
C
         DO 5 J=1,IDIME
            IKI = ((KMAX-1)/4*4 + 1 -KKMAX+2)/4
            DO 116 NY=1,NYIC
               IF (IGAM .EQ. 1)THEN
                  III = (KKMAX - KAX(NY))/4 + 1 + IKI
                  EIV(III,NY,J)=EIGVEL(NY,J)
               ELSE
                  DO 16 LLL=1,IDIM
                     LL=LLL+(NY-1)*IDIM
                     L = LLL + IKI
                     EIV(L,NY,J)=EIGVEL(LL,J)
   16             CONTINUE
               ENDIF
  116       CONTINUE
            JKJ = JKJ + 1
            AF(JKJ) = VALUE(J)
            IAM(JKJ) = IMOM
            DO 6 L=1,IDIM1
               V(L)=ABS(EIGVEL(L,J))
    6       CONTINUE
            VMAX=V(1)
            LL=1
            IF(IDIM1.EQ.1) GO TO 8
            DO 9 L=2,IDIM1
               IF(V(L).LT.VMAX) GO TO 9
               VMAX=V(L)
               LL=L
    9       CONTINUE
    8       CONTINUE
            KAV(JKJ)=LL
C
C********** EXPECTATION VALUES OF ANGULAR MOMENTUM OPERATORS  **********
C
            CALL EXPECT(JZAV(JKJ),JX2AV(JKJ),JY2AV(JKJ),JZ2AV(JKJ),
     &       RX2AV(JKJ),RY2AV(JKJ),RZ2AV(JKJ),RDOTJ(JKJ),IX2AV(JKJ),
     &     IY2AV(JKJ),IZ2AV(JKJ),IDOTJ(JKJ),HX(JKJ),HY(JKJ),HZ(JKJ))
C
C********** EXPECTATION VALUES OF CORE ROTATIONAL HAMILTONIAN **********
C
            IF (CHSI .NE. 1. .AND. VMI .NE. 1)
     &          CALL HEXPECT(HX(JKJ),HY(JKJ),HZ(JKJ))
            IF (VMI .EQ. 1) THEN
                HX(JKJ)=0.
                HY(JKJ)=0.
                HZ(JKJ)=0.
            ENDIF
    5    CONTINUE
         LEOUT = MIN0(NOUTJ(II),IDIME)
         LE = MIN0(5,IDIME)
         DO 301 LC=1,LEOUT
            write(6, 28)
            NYA=0
            DO 601 NY=1,NYIC
               EIGMAX = 0.00001
               IF (IGAM .EQ. 1) THEN
                   ABSEIV = ABS(EIGVEL(NY,LC))
                   EIGMAX = MAX(EIGMAX,ABSEIV)
                   IF (EIGMAX .LT. 0.05) GO TO 601
                   NYA = NYA+1
                   IF (NYA .EQ. 1) write(6, 351) IMOM,KATT,LC,VALUE(LC),
     *                             LEVEL(NY),EIGVEL(NY,LC)
                   IF (NYA .GE. 2) write(6, 353) LEVEL(NY),EIGVEL(NY,LC)
               ELSE
                   DO 602 KLL = KLARGE,KSMALL
                      ABSEIV = ABS(EIV(KLL,NY,LC))
                      EIGMAX = MAX(EIGMAX,ABSEIV)
  602              CONTINUE
                   IF (EIGMAX .LT. 0.05) GO TO 601
                       NYA = NYA+1
                       IF (NYA.EQ.1) write(6,351)IMOM,KATT,LC,VALUE(LC),
     *                      LEVEL(NY),(EIV(KLL,NY,LC),KLL=KLARGE,KSMALL)
                       IF (NYA .GE. 2) write(6, 353)
     *                      LEVEL(NY),(EIV(KLL,NY,LC),KLL=KLARGE,KSMALL)
               ENDIF
  601       CONTINUE
  353       FORMAT(22X,I2,1H:,15F6.3)
  351       FORMAT(1X,'I=',I3,'/2',A1,' (',I1,')',F7.3,I4,1H:,15F6.3)
C
C******************************************************************
C
  301    CONTINUE
         WRITE(18) LE
         WRITE(18) (VALUE(JJ),MOMENT(JJ),JJ=1,LE),
     $   (((EIV(KLL,NY,JJ),KLL=1,15),NY=1,NYIC),JJ=1,LE),
     &   ((PUV(MY,NY),MY=1,NYIC),NY=1,NYIC)
 4000 CONTINUE
C
C
!      JRK=MSLEFT(XX)
      TID=(IRK-JRK)*0.001
!      call cpu_time(end_time)
!	elapsed_time=end_time-begin_time
!	TID=elapsed_time

      write(6, 17) (K3(L),L=1,IDIM)
      NI = JKJ
      IDIM=NI
      write(6, 234) NI,TID
  234 FORMAT(1X,I3,' LEVELS    TIME=',F5.1,' SEC')
C
C**** PUTTING THE ENERGY EIGENVALUES IN ORDER
C
CC    CALL SORT6(AF,IAM,KAV,JZAV,JZ2AV,IDOTJ,NI)
      CALL SORT18(AF,IAM,KAV,JZAV,JX2AV,JY2AV,JZ2AV,RX2AV,RY2AV,
     &   RZ2AV,RDOTJ,IX2AV,IY2AV,IZ2AV,IDOTJ,HX,HY,HZ,NI)
      write(6, 350)
CC350 FORMAT(///,' ******************* EXCITATION ENERGIES , DOMINANT ',
CC   *'CONFIGURATION AND EXPECTATION VALUES  *******************',//,
CC   * '  NR  Eexc(keV)  E(keV)    I     K   NY ',
CC   *'    <Jz>  <Jz**2>  <I.j>/<I>',/)
  350 FORMAT(///,' ******************* EXCITATION ENERGIES, DOMINANT ',
     *'CONFIGURATION AND EXPECTATION VALUES (WITH S.P. UU+VV FACTORS) ',
     *'*******************',//,
     * '  NR  Eex(keV)  I   Kdom  NY    <Jz>  <Jx**2> <Jy**2> <Jz**2>',
     *' <I.j>/|I|  <Rx**2> <Ry**2> <Rz**2> <R.j>/|R|',
     *'  <HROTx> <HROTy> <HROTz>',/)
C23456789012345678901234567890123456789012345678901234567890123456789012
CC350 FORMAT(///,' ******************* EXCITATION ENERGIES , DOMINANT ',
CC   *'CONFIGURATION AND EXPECTATION VALUES  *******************',//,
CC   * '  NR  Eexc(keV)  E(keV)    I     K   NY ',
CC   *'    <Jz>  <Jz**2>  <I.j>/<I>',/)
C 350 FORMAT(///,'  EXCITATION ENERGIES AND DOMINANT CONFIGURATION',//,
C    * 3('  NR  Eexc  E (keV)    I     K   NY      '),/)
      EMIN=AF(1)
      DO 7 J=1,NI
         AE(J)=AF(J)*1000.0
         AF(J)=(AF(J)-EMIN)*1000.0
         IF (IGAM .EQ. 1) THEN
            KL=KAV(J)
            NYDOM(J)=LEVEL(KL)
            KDOM(J)=KAX(KL)
         ELSE
            KKMAX=MIN0(IAM(J),KMAX)
            IJL2=(KKMAX+1)/2
            KL=(KAV(J)-1)/IJL2
            KL=KL+1
            KL=LEVEL(KL)
            NYDOM(J) = KL
            KLL=MOD(KAV(J),IJL2)
            IF(KLL.EQ.0) KLL=IJL2
            IF(IAM(J).EQ.1) KLL=1
            IF((IAM(J)-1)/4*4.NE.IAM(J)-1) KKMAX=IAM(J)-2
            KKMAX=MIN0(KKMAX,KMAX)
            KK3=KKMAX-4*(KLL-1)
            KVL=IABS(KK3)
            KDOM(J) = KVL
         ENDIF
    7 CONTINUE
C     LINES = (NI+2)/3
C     DO 107 K=1,LINES
C        KMX = K + 2*LINES
C        IF (KMX .GT. NI) KMX = KMX-LINES
C        PRINT 400,(J,AF(J),AE(J),IAM(J),KATT,KDOM(J),NYDOM(J),
C    $   J=K,KMX,LINES)
C 107 CONTINUE
      DO 107 J=1,NI
CC       PRINT 400,J,AF(J),AE(J),IAM(J),KATT,KDOM(J),NYDOM(J),
CC   $   JZAV(J),JZ2AV(J),IDOTJ(J)
         write(6, 400) J,AF(J),IAM(J),KATT,KDOM(J),NYDOM(J),
     $   JZAV(J),JX2AV(J),JY2AV(J),JZ2AV(J),IDOTJ(J),
     $   RX2AV(J),RY2AV(J),RZ2AV(J),RDOTJ(J),HX(J),HY(J),HZ(J)
      IF (MOD(J,5) .EQ. 0) WRITE(6,*) ' '
  107 CONTINUE
C 400 FORMAT(3(1X,I3,2F7.1,I4,2H/2,A1,I4,2H/2,I4,6X))
CC400 FORMAT(1X,I3,2F9.1,I4,2H/2,A1,I4,2H/2,I4,2X,3F8.2)
  400 FORMAT(1X,I3,F7.1,I4,2H/2,A1,I4,2H/2,I4,2X,F6.2,4(1X,F7.2),
     &            3X,4(1X,F7.2),3X,3(1X,F7.3))
  111 CONTINUE
      STOP
      END
C
C#######################################################################
C
      SUBROUTINE EXPECT(JZAV,JX2AV,JY2AV,JZ2AV,RX2AV,RY2AV,RZ2AV,
     &                  TRDOTJ,IX2AV,IY2AV,IZ2AV,TIDOTJ,HX,HY,HZ)
C
C****** EXPECTATION VALUES OF jz, jx**2, jy**2, jz**2,
C                             Rx**2, Ry**2, Rz**2, R.j/SQRT(R**2)
C                             Ix**2, Iy**2, Iz**2,
C                             and I.j/SQRT(I*(I+1)) ******
C       INCLUDING THEIR SINGLE-PARTICLE PAIRING FACTORS UU' + VV'
C       (NOTE THAT THE PAIRING FACTORS ARE INCLUDED IN APLB, ETC)
C
C NOTE THAT OTHER EXPECTATION VALUES ARE ALSO CALCULATED, SUCH AS
C     Iz, Ixjx, Iyjy, Izjz, I+j-, I-j+, I+j+, I-j-, I+I+, I-I-, I+I-
C     AND EQPAV, THE EXPECTATION VALUE OF THE QUASIPARTICLE HAMILTONIAN
C
C  FOR CHSI=1. AND VMI .NE. 1, THE EXPECTATION VALUES OF THE CORE ROT'L
C       ENERGY ABOUT THE X, Y, AND Z AXES ARE ALSO CALCULATED CORRECTLY
C
      COMMON/INT1/ IMOM,NYIC,LC,KLARGE,KSMALL,KMAX
      COMMON/EXPECT1/ PUV(15,15),SQP(15),A1,A2,A3,CHSI
      COMMON/DECOUP/ APLB(15,15),C(15,15),D(15,15),E(15,15),
     $NORM(15,15),JZ(15,15),JZ2(15,15)
      COMMON/EIGVEC/ EIV(15,15,10)
      REAL NORM,JZ,JZ2
      REAL JZAV,JZ2AV,JX2AV,JY2AV,J2AV,IZAV,IX2AV,IY2AV,IZ2AV
      REAL I2,IXJXAV,IYJYAV,IZJZAV,IPLJMN,IMNJPL,IPLJPL,IMNJMN
      REAL IPL2,IMN2,IPLIMN,IDOTJ
      JZAV=0.
      JZ2AV=0.
      JX2AV=0.
      JY2AV=0.
      IZAV=0.
      IZ2AV=0.
      IZJZAV=0.
      IPLJMN=0.
      IMNJPL=0.
      IPLJPL=0.
      IMNJMN=0.
      IPL2=0.
      IMN2=0.
      IPLIMN=0.
      EQPAV=0.
      I2=FLOAT(IMOM*(IMOM+2))/4.
      N=(IMOM+1)/2
      DO 701 NY=1,NYIC
      DO 702 NYP=1,NYIC
           DELTA=0.0
           IF (NY .EQ. NYP) DELTA=1.0
           EIGS12 = 0.

           DO 703 KLP=KLARGE,KSMALL
              KP = KMAX - (KLP-1)*4
              EIGS12 = EIGS12 + EIV(KLP,NYP,LC)*EIV(KLP,NY,LC)
           DO 704 KL =KLARGE,KSMALL
              K  = KMAX - (KL -1)*4
              EIG12 = EIV(KLP,NYP,LC)*EIV(KL,NY,LC)
CC
CC** SPECIAL TREATMENT FOR K=1/2
CC
CC            IF (K .NE. 1) GOTO 705
CC            IF (KP .EQ. 1) THEN
CC                IMNJPL = IMNJPL - EIG12*
CC   &                     PUV(NY,NYP)*(APLB(NY,NYP)*(-1)**N*N)
CC                IMNJMN = IMNJMN - EIG12*
CC   &                     PUV(NY,NYP)*(C(NY,NYP)*(-1)**N*N)
CC            ELSEIF (KP .EQ.-3) THEN
CC                IPLJMN = IPLJMN - EIG12*PUV(NY,NYP)*
CC   &              C(NY,NYP)*(-1)**N*SQRT((IMOM-1.)*(IMOM+3.)/4.)
CC                IPLJPL = IPLJPL - EIG12*PUV(NY,NYP)*
CC   &              APLB(NY,NYP)*(-1)**N*SQRT((IMOM-1.)*(IMOM+3.)/4.)
CC            ENDIF
CC            GO TO 706
CC705 CONTINUE
C
C** MATRIX ELEMENTS  (-1)**(I-1/2)*<I,KP,NYP|OP|I,-K,NY(CONJUGATE)>
C** INCLUDES THE SPECIAL CASE K=1/2
C
              IF (KP .EQ. -K+2) THEN
                  IMNJPL = IMNJPL - EIG12*APLB(NY,NYP)*
     *                     (-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K-2.)/4.)
                  IMNJMN = IMNJMN - EIG12*C(NY,NYP)*
     *                     (-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K-2.)/4.)

              ELSEIF (KP .EQ. -K-2) THEN
                   IPLJMN = IPLJMN - EIG12*C(NY,NYP)*
     *                      (-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K+2.)/4.)
                   IPLJPL = IPLJPL - EIG12*APLB(NY,NYP)*
     *                      (-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K+2.)/4.)
              ENDIF
  706         CONTINUE
C
C** MATRIX ELEMENTS  <I,KP,NYP|OP|I,K,NY>
C
              IF (KP.EQ.K) THEN
                 IZJZAV = IZJZAV + EIG12*K*JZ(NY,NYP)/2.
                 IZAV = IZAV + EIG12*DELTA*K/2.
                 IZ2AV = IZ2AV + EIG12*DELTA*K*K/4.
                 IPLIMN = IPLIMN + EIG12*DELTA*
     &                    FLOAT(IMOM*(IMOM+2)-K*(K+2))/4.
              ELSEIF (KP .EQ. K-4) THEN
                 IPL2 = IPL2 + EIG12*DELTA*SQRT(FLOAT((IMOM+K)*
     &                     (IMOM-K+2)*(IMOM+K-2)*(IMOM-K+4)))/4.
              ELSEIF (KP .EQ. K+4) THEN
                 IMN2 = IMN2 + EIG12*DELTA*SQRT(FLOAT((IMOM-K)*
     &                     (IMOM+K+2)*(IMOM-K-2)*(IMOM+K+4)))/4.
              ENDIF
  704    CONTINUE
  703    CONTINUE
         JZAV =  JZAV  + EIGS12*JZ(NY,NYP)
         JZ2AV = JZ2AV + EIGS12*JZ2(NY,NYP)
         JX2AV = JX2AV + EIGS12*(D(NY,NYP)/2. + E(NY,NYP)/4.)
         JY2AV = JY2AV + EIGS12*(D(NY,NYP)/2. - E(NY,NYP)/4.)
         EQPAV = EQPAV + EIGS12*DELTA*SQP(NY)
  702 CONTINUE
  701 CONTINUE
      IXJXAV = (IPLJPL + IMNJMN + IPLJMN + IMNJPL)/4.
      IYJYAV = (IPLJMN + IMNJPL - IMNJMN - IPLJPL)/4.
      IX2AV = (IPL2 + IMN2)/4. + (IPLIMN + IZAV)/2.
      IY2AV = (IPLIMN + IZAV)/2. - (IPL2 + IMN2)/4.
      IDOTJ = (IMNJPL + IPLJMN)/2. + IZJZAV
      RX2AV = IX2AV + JX2AV -2.*IXJXAV
      RY2AV = IY2AV + JY2AV -2.*IYJYAV
      RZ2AV = IZ2AV + JZ2AV -2.*IZJZAV

      J2AV = JX2AV + JY2AV + JZ2AV
      RDOTJ = IDOTJ - J2AV
      R2AV = I2 + J2AV - 2.*IDOTJ
C**
C**   Occasionally R2AV is neg when used with the 2-body recoil
C**   treatment of the recoil term because of a "Core contribution"
C**   that has been dropped. Test here
C**
      IF (R2AV .LT. 0.) THEN
         write(6, 999) IMOM
         R2AV = -R2AV
      ENDIF
  999 FORMAT(/10X,'I =',I3,'/2',5X,'NEGATIVE R2AV IN EXPECT')
C***
      TRDOTJ = RDOTJ/SQRT(R2AV)
      TIDOTJ = IDOTJ/SQRT(I2)
C
C** EXPECTATION VALUES OF THE CORE ROTATION HAMILTIONIAN
C**  (NOTE THAT THIS IS CORRECT ONLY IF CHSI = 1.)
C
      HX = A1*RX2AV
      HY = A2*RY2AV
      HZ = A3*RZ2AV

CC    write(6, 721) JZAV,JX2AV,JY2AV,JZ2AV,IZAV,IX2AV,IY2AV,IZ2AV,TIDOTJ,
CC   &          RX2AV,RY2AV,RZ2AV,TRDOTJ,IXJXAV,IYJYAV,IZJZAV,EQPAV
  721 FORMAT(/'  <JZ>=',F7.3,'  <JX**2>=',F7.3,'  <JY**2>=',F7.3,
     &       '  <JZ**2>=',F7.3,'  <IZ>=',F7.3,'  <IX**2>=',F7.3,
     &       '  <IY**2>=',F7.3,/
     &       '  <IZ**2>=',F7.3,'  <I.j>/|I|=',F7.3,'  <RX**2>=',F7.3,
     &       '  <RY**2>=',F7.3,'  <RZ**2>=',F7.3,'  <R.j>/|R|=',F7.3,/
     &       '  <IXjX>=',F7.3,'  <IYjY>=',F7.3,'  <Izjz>=',F7.3,
     &       '  <EQP>=',F7.3)
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE HEXPECT(HX,HY,HZ)
C
C****** EXPECTATION VALUES OF THE CORE ROTATIONAL ENERGY ABOUT THE
C       X, Y, AND Z AXES (<Hx>, <Hy>, <Hz>) IF VMI .NE.1
C       ALSO CALCULATES (BUT NOT PRINTED) THE EXPECTATION VALUE OF
C       THE QUASIPARTICLE HAMILTONIAN (EQPAV)
C
C       NOTE THAT JZ2AV ETC, CONTAIN CHSI AND THUS DIFFER FROM
C                THE EXPECTATION VALUES CALCULATED IN "EXPECT"
C
C**  PUV REMOVED HERE SINCE THESE FACTORS ARE NOW INCLUDED IN APLB, ETC
C
      COMMON/INT1/ IMOM,NYIC,LC,KLARGE,KSMALL,KMAX
      COMMON/EXPECT1/ PUV(15,15),SQP(15),A1,A2,A3,CHSI
      COMMON/DECOUP/ APLB(15,15),C(15,15),D(15,15),E(15,15),
     $NORM(15,15),JZ(15,15),JZ2(15,15)
      COMMON/EIGVEC/ EIV(15,15,10)
      REAL NORM,JZ,JZ2
      REAL JZAV,JZ2AV,JX2AV,JY2AV,J2AV,IZAV,IX2AV,IY2AV,IZ2AV
      REAL I2,IXJXAV,IYJYAV,IZJZAV,IPLJMN,IMNJPL,IPLJPL,IMNJMN
      REAL IPL2,IMN2,IPLIMN,IDOTJ
CC    JZAV=0.
      JZ2AV=0.
      JX2AV=0.
      JY2AV=0.
      IZAV=0.
      IZ2AV=0.
      IZJZAV=0.
      IPLJMN=0.
      IMNJPL=0.
      IPLJPL=0.
      IMNJMN=0.
      IPL2=0.
      IMN2=0.
      IPLIMN=0.
      EQPAV=0.
      I2=FLOAT(IMOM*(IMOM+2))/4.

      N=(IMOM+1)/2
      DO 701 NY=1,NYIC
      DO 702 NYP=1,NYIC
           DELTA=0.0
           IF (NY .EQ. NYP) DELTA=1.0
           CHSI1=CHSI
           IF (NY .EQ. NYP) CHSI1=1.0
           EIGS12 = 0.

           DO 703 KLP=KLARGE,KSMALL
              KP = KMAX - (KLP-1)*4
              EIGS12 = EIGS12 + EIV(KLP,NYP,LC)*EIV(KLP,NY,LC)
           DO 704 KL =KLARGE,KSMALL
              K  = KMAX - (KL -1)*4
              EIG12 = EIV(KLP,NYP,LC)*EIV(KL,NY,LC)
C
C** MATRIX ELEMENTS  (-1)**(I-1/2)*<I,KP,NYP|OP|I,-K,NY(CONJUGATE)>
C** INCLUDES THE SPECIAL CASE K=1/2
C
              IF (KP .EQ. -K+2) THEN
                  IMNJPL = IMNJPL - EIG12*APLB(NY,NYP)*
     *              CHSI1*(-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K-2.)/4.)
                  IMNJMN = IMNJMN - EIG12*C(NY,NYP)*
     *              CHSI1*(-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K-2.)/4.)
              ELSEIF (KP .EQ. -K-2) THEN
                   IPLJMN = IPLJMN - EIG12*C(NY,NYP)*
     *               CHSI1*(-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K+2.)/4.)
                   IPLJPL = IPLJPL - EIG12*APLB(NY,NYP)*
     *               CHSI1*(-1)**N*SQRT(IMOM*(IMOM+2.)/4.-K*(K+2.)/4.)
              ENDIF
  706         CONTINUE
C
C** MATRIX ELEMENTS  <I,KP,NYP|OP|I,K,NY>
C
              IF (KP.EQ.K) THEN
                 IZJZAV = IZJZAV + EIG12*K*JZ(NY,NYP)/2.*CHSI1
                 IZAV = IZAV + EIG12*DELTA*K/2.
                 IZ2AV = IZ2AV + EIG12*DELTA*K*K/4.
                 IPLIMN = IPLIMN + EIG12*DELTA*
     &                    FLOAT(IMOM*(IMOM+2)-K*(K+2))/4.
              ELSEIF (KP .EQ. K-4) THEN
                 IPL2 = IPL2 + EIG12*DELTA*SQRT(FLOAT((IMOM+K)*
     &                     (IMOM-K+2)*(IMOM+K-2)*(IMOM-K+4)))/4.
              ELSEIF (KP .EQ. K+4) THEN
                 IMN2 = IMN2 + EIG12*DELTA*SQRT(FLOAT((IMOM-K)*
     &                     (IMOM+K+2)*(IMOM-K-2)*(IMOM+K+4)))/4.
              ENDIF
  704    CONTINUE
  703    CONTINUE
CC       JZAV =  JZAV  + EIGS12*JZ(NY,NYP)*CHSI1
         JZ2AV = JZ2AV + EIGS12*JZ2(NY,NYP)*CHSI1
         JX2AV = JX2AV + EIGS12*(D(NY,NYP)/2. +E(NY,NYP)/4.)*CHSI1
         JY2AV = JY2AV + EIGS12*(D(NY,NYP)/2. -E(NY,NYP)/4.)*CHSI1
         EQPAV = EQPAV + EIGS12*DELTA*SQP(NY)
  702 CONTINUE
  701 CONTINUE
      IXJXAV = (IPLJPL + IMNJMN + IPLJMN + IMNJPL)/4.
      IYJYAV = (IPLJMN + IMNJPL - IMNJMN - IPLJPL)/4.
      IX2AV = (IPL2 + IMN2)/4. + (IPLIMN + IZAV)/2.
      IY2AV = (IPLIMN + IZAV)/2. - (IPL2 + IMN2)/4.
      RX2AV = IX2AV + JX2AV -2.*IXJXAV
      RY2AV = IY2AV + JY2AV -2.*IYJYAV
      RZ2AV = IZ2AV + JZ2AV -2.*IZJZAV
CC    J2AV = JX2AV + JY2AV + JZ2AV
C
C** EXPECTATION VALUE OF CORE ROTATION HAMILTIONIAN
C**  (NOTE THAT THIS INCLUDES THE UU'+VV' FACTORS AS WELL AS CHSI)
C
CC    H1 = (A1+A2)/2.*(I2 - IZ2AV + J2AV - JZ2AV - IPLJMN - IMNJPL)
CC    H2 = (A1-A2)/4.*(IPL2 + IMN2 + 2.*(JX2AV - JY2AV)
CC   &                             - 2.*(IPLJPL + IMNJMN))
CC    H3 = A3*(IZ2AV + JZ2AV - 2.*IZJZAV)
      HX = A1*RX2AV
      HY = A2*RY2AV
      HZ = A3*RZ2AV
CC    PRINT 721, HX,HY,HZ,EQPAV
CC721 FORMAT(/'  <HX>=',F7.3,' <HY>=',F7.3,' <HZ>=',F7.3,' <EQP>=',F7.3)
      RETURN
      END

C
C#######################################################################
C
C
      SUBROUTINE SORT18(E,KV1,KV2,KV3,KV4,KV5,KV6,KV7,KV8,KV9,
     &                KV10,KV11,KV12,KV13,KV14,KV15,KV16,KV17,NI)
C
C     SORTING OF E, KV1,...., KV17 ACCORDING TO INCREASING VALUES IN E
C       NI: DIMENSION OF E, KV1,..., KV17
C
      DIMENSION E(1),KV1(1),KV2(1),KV3(1),KV4(1),KV5(1),
     &          KV6(1),KV7(1),KV8(1),KV9(1),KV10(1),KV11(1),
     &          KV12(1),KV13(1),KV14(1),KV15(1),KV16(1),KV17(1)
      REAL KV3,KV4,KV5,KV6,KV7,KV8,KV9,KV10,
     &     KV11,KV12,KV13,KV14,KV15,KV16,KV17
      INTEGER Y1,Y2
CT      WRITE(1,*)'UTSKRIFT I SORT15,NI',NI
      N=NI
      M=2**INT(1.4427*ALOG(FLOAT(N-1)))-1
    1 IF(M.LE.0) RETURN
      I=1
    2 J=I
      K=I+M
      Y=E(K)
      Y1=KV1(K)
      Y2=KV2(K)
      Y3=KV3(K)
      Y4=KV4(K)
      Y5=KV5(K)
      Y6=KV6(K)
      Y7=KV7(K)
      Y8=KV8(K)
      Y9=KV9(K)
      Y10=KV10(K)
      Y11=KV11(K)
      Y12=KV12(K)
      Y13=KV13(K)
      Y14=KV14(K)
      Y15=KV15(K)
      Y16=KV16(K)
      Y17=KV17(K)
    3 IF(Y.LT.E(J)) GO TO 5
    4 L=J+M
      E(L)=Y
      KV1(L)=Y1
      KV2(L)=Y2
      KV3(L)=Y3
      KV4(L)=Y4
      KV5(L)=Y5
      KV6(L)=Y6
      KV7(L)=Y7
      KV8(L)=Y8
      KV9(L)=Y9
      KV10(L)=Y10
      KV11(L)=Y11
      KV12(L)=Y12
      KV13(L)=Y13
      KV14(L)=Y14
      KV15(L)=Y15
      KV16(L)=Y16
      KV17(L)=Y17
      I=I+1
      IF(I+M.LE.N) GO TO 2
      M=(M-1)/2
      GO TO 1
    5 L=J+M
      E(L)=E(J)
      KV1(L)=KV1(J)
      KV2(L)=KV2(J)
      KV3(L)=KV3(J)
      KV4(L)=KV4(J)
      KV5(L)=KV5(J)
      KV6(L)=KV6(J)
      KV7(L)=KV7(J)
      KV8(L)=KV8(J)
      KV9(L)=KV9(J)
      KV10(L)=KV10(J)
      KV11(L)=KV11(J)
      KV12(L)=KV12(J)
      KV13(L)=KV13(J)
      KV14(L)=KV14(J)
      KV15(L)=KV15(J)
      KV16(L)=KV16(J)
      KV17(L)=KV17(J)
      J=J-M
      IF(J.GE.1) GO TO 3
      GO TO 4
      END
C
C#######################################################################
C**** SOLUTION OF SYSTEM OF LINEAR EQUATIONS  *******************
C
      SUBROUTINE GJR(A,NC,NR,N,MC,JC,V)
      DIMENSION A(10,10),JC(10),V(2)
C              ---------------------------------------------------------
C               JC IS THE PERMUTATION VECTOR
C              KD IS THE OPTION KEY FOR DETERMINANT EVALUATION
C              KI IS THE OPTION KEY FOR MATRIX INVERSION
C              L IS THE COLUMN CONTROL FOR AX=B
C              M IS THE COLUMN CONTOL FOR MATRIX INVERSION
C              ---------------------------------------------------------
C              INITIALIZATION
C              ---------------------------------------------------------
      IW=V(1)
      M=1
      S=1.
      L=N+(MC-N)*(IW/4)
      KD=2-MOD(IW/2,2)
      IF(KD.EQ.1) V(2)=0.
      KI=2-MOD(IW,2)
      GO TO (5,20),KI
C              ---------------------------------------------------------
C                  INITIALIZE JC FOR INVERSION
C              ---------------------------------------------------------
5     DO  10 I=1,N
10    JC(I)=I
C              ---------------------------------------------------------
C              SEARCH FOR PIVOT ROW
C              ---------------------------------------------------------
20    DO 91  I=1,N
      GO TO (22,21),KI
21    M=I
22    IF (I.EQ.N)  GO TO 60
      X=-1.
      DO 30  J=I,N
      IF (X.GT.ABS(A(J,I)))  GO TO 30
      X=ABS(A(J,I))
      K=J
30    CONTINUE
      IF(K.EQ.I) GO TO 60
      S=-S
      V(1)=-V(1)
      GO TO (35,40),KI
 35   MU=JC(I)
      JC(I)=JC(K)
      JC(K)=MU
C              ---------------------------------------------------------
C              INTERCHANGE ROW I AND ROW K
C              ---------------------------------------------------------
40    DO 50  J=M,L
      X=A(I,J)
      A(I,J)=A(K,J)
50    A(K,J)=X
C              ---------------------------------------------------------
C              TEST FOR SINGULARITY
C              ---------------------------------------------------------
60    IF (ABS(A(I,I)).GT.0.)  GO TO 70
C              ---------------------------------------------------------
C              MATRIX IS SINGULAR
C              ---------------------------------------------------------
      write(6, 211)
  211 FORMAT(//,'  SINGULAR MATRIX FOUND IN GJR')
      STOP
70    GO TO (71,72),KD
C              ---------------------------------------------------------
C              COMPUTE THE DETERMINANT
C              ---------------------------------------------------------
71    IF(A(I,I).LT.0.) S=-S
      V(2)=V(2)+ALOG(ABS(A(I,I)))
72    X=A(I,I)
      A(I,I)=1.
C              ---------------------------------------------------------
C              REDUCTION OF THE I-TH ROW
C              ---------------------------------------------------------
      DO 80 J=M,L
      A(I,J)=A(I,J)/X
C              ---------------------------------------------------------
C              TEST OVERFLOW SWITCH. IF ON
C              RETURN NEGATIVE VALUE OF I IN JC(1)
C              ---------------------------------------------------------
C    OVERFLOW TEST TAKEN AWAY
80    CONTINUE
C              ---------------------------------------------------------
C              REDUCTION OF ALL REMAINING ROWS
C              ---------------------------------------------------------
      DO 91 K=1,N
      IF (K.EQ.I)  GO TO 91
      X=A(K,I)
      A(K,I)=0.
      DO 90  J =M,L
      A(K,J)=A(K,J)-X*A(I,J)
C              -------------------------------------------------------
C              -------------------------------------------------------
C              TEST OVERFLOW SWITCH. IF ON
C              RETURN NEGATIVE VALUE OF I IN JC(1)
C              ---------------------------------------------------------
C    OVERFLOW TEST TAKEN AWAY
90    CONTINUE
91    CONTINUE
C              ---------------------------------------------------------
C              AX=B AND DET.(A) ARE NOW COMPUTED
C              ---------------------------------------------------------
      GO TO (95,140),KI
C              ---------------------------------------------------------
C              PERMUTATION OF THE COLUMNS FOR MATRIX INVERSION
C              ---------------------------------------------------------
95    DO 130  J=1,N
      IF (JC(J).EQ.J)  GO TO 130
      JJ=J+1
      DO 100  I=JJ,N
      IF (JC(I).EQ.J)  GO TO 110
100   CONTINUE
110   JC(I)=JC(J)
      DO 120  K=1,N
      X=A(K,I)
      A(K,I)=A(K,J)
120   A(K,J)=X
130   CONTINUE
140    JC(1)=N
      IF(KD.EQ.1) V(1)=S
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE PAIRQ2(IEMIN,IEMAX,IFIT,ISELF,LFIX,AMAS,TOLER
     *,OMEGA,ENE,GM,G2,ODDE,AD0,AD2,ALA,G0,WE,WWE)
C
C     $$ BCS-CALCULATION FOR MONOPOLE OR MONOPOLE+QUADRUPOLE PAIRING $$
C     $$ ALSO POSSIBILITIES TO DETERMINE G FOR INPUT OF DELTA  $$
C
      COMMON/SP2/  ESP(286),EREN(286), Q(286),V2(286),UV(286),EQP(286)
      DIMENSION X(10,10),FV(2),JC(10)
      NV=IEMAX-IEMIN+1
      IF(NV.GT.0) GO TO 2
    9 FORMAT(1X,'       ERROR IN PAIRQ2      ')
      write(6, 9)
      STOP
    2 CONTINUE
      G0=GM/AMAS
      IF(ABS(G0).LT.1.E-10) G0=1.E-10
      IF(ABS(G2).LT.1.E-10) G2=1.E-10
      NA=0
 220  NA=NA+1
      IF(NA.LT.100) GOTO 201
      write(6, 18)
   18 FORMAT(1X,'       MORE THAN 100 ITERATIONS IN PAIRQ2      ')
      STOP
 201  CONTINUE
      NUM=1
      GM=G0*AMAS
C  ***** DETERMINATION OF DELTA0,DELTA2 AND LAMBDA WHEN G0 AND G2 ARE FI
 250  CONTINUE
      Y0=0.
      YY0=0.
      Y2=0.
      AK1=0.
      A0=0.
      A1=0.
      A2=0.
      B0=0.
      B1=0.
      B2=0.
      B3=0.
      B4=0.
      AK1T=0.
      A0T=0.
      A1T=0.
      A2T=0.
      B0T=0.
      B2T=0.
      DO 260 I=IEMIN,IEMAX
      DELTA=AD0+Q(I)*AD2
      SIG=1.
      IF(DELTA.LT.0.) SIG=-1.
      Q2=Q(I)*Q(I)
      EEP=ESP(I)-ALA
      EPE=EEP
      IX=0
C  *** IF ISELF .NE. 0 THE SELFCONSISTENT CONTRIBUTIONS TO
C ***  THE SINGLE PARTICLE ENERGIES ARE INCLUDED
      IF(ISELF .EQ. 0) GOTO 255
 251  IX=IX+1
      VV=0.5*(1.-EEP/SQRT(EEP*EEP+DELTA*DELTA))
      EEP=ESP(I)-ALA-(G0+G2*Q2)*VV
      IF(ABS(EEP-EPE).LE.0.001) GOTO 255
      EPE=EEP
      IF(IX.LT.11) GOTO 251
      EEP=ESP(I)-ALA
      EEP=EEP-(G0+G2*Q2)*0.5*(1.-EEP/SQRT(EEP*EEP+DELTA*DELTA))
      write(6, 25) I
   25 FORMAT('   BAD RENORMALISATION IN PAIRQ2 FOR I=',I3)
 255  EREN(I)=EEP+ALA
      E2=EEP*EEP+DELTA*DELTA
      E1=SQRT(E2)
      ES1=E1*SIG
      E3=E1*E2*SIG
      Y0=Y0+1./ES1
      YY0=YY0+EEP/E1
      Y2=Y2+Q2/ES1
      AK1=AK1+Q(I)/ES1
      AK1T=AK1T+Q(I)/ES1*SIG
      A0=A0+EEP/E3
      A1=A1+EEP*Q(I)/E3
      A2=A2+EEP*Q2/E3
      B0=B0+1./E3
      B1=B1+Q(I)/E3
      B2=B2+Q2/E3
      B3=B3+Q(I)*Q2/E3
      B4=B4+Q2*Q2/E3
      A0T=A0T+EEP/E3*SIG
      A1T=A1T+EEP*Q(I)/E3*SIG
      A2T=A2T+EEP*Q2/E3*SIG
      B0T=B0T+1./E3*SIG
      B2T=B2T+Q2/E3*SIG
 260  CONTINUE
      AK0=2./G0-Y0
      AK2=2./G2-Y2
      X(1,4)=AK0*AK2-AK1*AK1
      X(2,4)=AD0*AK0-AD2*AK1
      X(3,4)=YY0+ENE-OMEGA*NV
      X(1,1)=2.*A1*AK1+AK2*A0+AK0*A2
      X(1,2)=-2.*AK1*(B1*AD0+AD2*B2)-AK2*(AD0*B0+AD2*B1)-AK0*(AD0*B2+
     1AD2*B3)
      X(1,3)=-2.*AK1*(B2*AD0+AD2*B3)-AK2*(AD0*B1+AD2*B2)-AK0*(AD0*B3+
     1AD2*B4)
      X(2,1)=AD0*A0+AD2*A1
      X(2,2)=-AK0-B0*AD0*AD0-2.*B1*AD0*AD2-B2*AD2*AD2
      X(2,3)=+AK1-B1*AD0*AD0-2.*B2*AD0*AD2-B3*AD2*AD2
      X(3,1)=AD0*AD0*B0+2.*AD0*AD2*B1+AD2*AD2*B2
      X(3,2)=AD0*A0+AD2*A1
      X(3,3)=AD0*A1+AD2*A2
      IF ((NA .EQ. 1) .AND. (NUM .LE. 3)) GO TO 501
      FV(1)=4.
      CALL GJR(X,4,3,3,4,JC,FV)
      IF (JC(1) .NE. 3) write(6, 10)
   10 FORMAT(1X,'       ERROR IN GJR CALLED FROM PAIRQ2       ')
      GO TO 502
  501 CONTINUE
      D=B0*B0*AD0*AD0*AD0+A0*A0*AD0
      DAD0=-(B0*AD0*AD0*AK0-A0*X(3,4))/D
      DALA= (X(3,4)*AD0*B0+AK0*AD0*A0)/D
      X(1,4)=DALA
      X(2,4)=DAD0
      X(3,4)=0.
      IF (NUM .NE. 3) GO TO 502
      AD2=0.
      DO 503 I=IEMIN,IEMAX
      EEP=EREN(I)-ALA-DALA
      E2=EEP*EEP+(AD0+DAD0)**2
      V=0.5*(1.-EEP/SQRT(E2))
      AD2=AD2+G2*Q(I)*SQRT((1.-V)*V)
  503 CONTINUE
  502 CONTINUE
      ALA=ALA+X(1,4)
      AD0=AD0+X(2,4)
      AD2=AD2+X(3,4)
      IF (G2 .LT. 0.000000001) AD2=0.
      IF (ABS(X(1,4)) .LT. 2.) GO TO 274
   13 FORMAT('0            SORRY I FOUND NO SOLUTION            ',
     *//,35X,' PAIRQ2')
      write(6, 13)
      RETURN
  274 CONTINUE
C  ***  NUM COUNTS THE NUMBER OF ITERATIONS WHEN G0 FIXED
      NUM=NUM+1
      IF ((NA .EQ. 1) .AND. (NUM .LE. 4)) GO TO 250
      IF(NUM.LT.100) GO TO 200
      write(6, 18)
      STOP
  200 CONTINUE
      ABSD=ABS(X(1,4))
      BBSD=ABS(X(2,4))
      CBSD=ABS(X(3,4))
      ABCD=AMAX1(BBSD,CBSD)
      ABSD=AMAX1(ABSD,ABCD)
      IF(ABSD.GT.TOLER) GOTO 250
      WE=0.
      SUM=0.
      SUN=0.
      V2S=0.
      Q2V4S=0.
      V4S=0.
      IEMIN1=IEMIN-1
       DO 279 I=1,IEMIN1
      V2(I)=1.
  279 CONTINUE
      DO 280 I=IEMIN,IEMAX
      DELTA=AD0+Q(I)*AD2
      EEP=EREN(I)-ALA
      E2=EEP*EEP+DELTA*DELTA
      EQP(I)=SQRT(E2)
      V2(I)=0.5*(1.-EEP/EQP(I))
      U2=1.-V2(I)
      UV(I)=SQRT(U2*V2(I))
      WE=WE+ ESP(I)*V2(I)
      SUM=SUM+UV(I)
      SUN=SUN+UV(I)*Q(I)
      V2S=V2S+V2(I)
      Q2V4S=Q2V4S+V2(I)*V2(I)*Q(I)*Q(I)
      V4S=V4S+V2(I)*V2(I)
  280 CONTINUE
      DEL1=G0*SUM
      DEL2=G2*SUN
      WE=2.*WE-DEL1*SUM-DEL2*SUN
      WWE=-G0*V4S-G2*Q2V4S
      WE=WE+WWE
C *** WE IS THE TOTAL ENERGY AMONG THE PAIRED PARTICLES   **************
C *** IF IFIT .NE. 0 WE NOW FIT 'DELTA' TO THE ODD-EVEN MASS DIFF.
C *** IF IFIT=0 WE ARE READY
      IF (IFIT .EQ. 0) GO TO 300
      EMIN=100.
      DO 290 I=IEMIN,IEMAX
      DMIN=EMIN
      EMIN=AMIN1(EMIN,EQP(I))
      IF(ABS(DMIN-EMIN).GT.1.D-30) II=I
 290  CONTINUE
      IF(LFIX.GT.0) EMIN=EQP(LFIX)
      IF(LFIX.GT.0) II=LFIX
      DE1=ODDE-EMIN
      IF (LFIX .EQ. -1) DE1=ODDE-AD0
      DE2=ODDE*ODDE-EMIN*EMIN
      IGO=0
      IF(ABS(DE1).LE..0005) IGO=1
      IF(IGO.EQ.1) GOTO 300
      X(1,4)=0.
      X(2,4)=2.*AK2/G0/G0
      X(3,4)=-AD0*Y0-AD2*AK1
      X(1,1)=AD0*AD0*B0T+2.*AD0*AD2*AK1T+AD2*AD2*B2T
      X(1,2)=AD0*A0T+AD2*A1T
      X(1,3)=AD0*A1T+AD2*A2T
      X(2,1)=-A1*AK2-AK0*A2-2.*AK1*A1
      X(2,2)=AD0*AK2*B0+AD2*AK2*B1+AD0*AK0*B2+AD2*AK0*B3+2.*AK1*AD0*B1+
     12.*AD2*AK1*B2
      X(2,3)=AD0*AK2*B1+AD2*AK2*B2+AK0*AD0*B3+AD2*AK0*B4+2.*AK1*AD0*B2+
     12.*AD2*AK1*B3
      X(3,1)=AD0*G0*A0+G0*AD2*A1
      X(3,2)=-G0*AK0-AD0*AD0*G0*B0-2.*G0*AD0*AD2*B1-AD2*AD2*G0*B2
      X(3,3)=-AD0*AD0*G0*B1-2.*AD0*AD2*G0*B2+G0*AK1-G0*AD2*AD2*B3
      FV(1)=4.
      CALL GJR(X,4,3,3,4,JC,FV)
      IF (JC(1) .NE. 3) write(6, 10)
      DERV=2.*(AD0+AD2*Q(II))*(X(2,4)+Q(II)*X(3,4))-2.*(EREN(II)-ALA)*
     1X(1,4)
      IF (LFIX .EQ. -1) DERV=X(2,4)
      IF (LFIX .EQ. -1) DE2=ODDE-AD0
      DG0=DE2/DERV
      IF (NA .GT. 15) DG0=DG0/2.
      G0=G0+DG0
      GOTO 220
 300  CONTINUE
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE BIGMAX(A,VALU,NN,NEIG,NVEC,VECTR,MDIM)
C
C     MATRIX DIAGONALISATION USING HOUSEHOULDER'S METHOD
C
C         A: MATRIX TO BE DIAGONALIZED - TRIANGULAR PART STORED IN
C            ONE DIMENSION; (NN+1)*NN/2
C         VALU: EIGENVALUES ORDERED WITH LARGEST ONE FIRST
C         NN: DIMENSION OF MATRIX; "NN * NN"
C         NEIG: NUMBER OF EIGENVALUES CALCULATED
C         NVEC: NUMBER OF EIGENVECTORS CALCULATED
C         VECTR: EIGENVECTORS
C         MDIM: FIRST DIMENSION AS DECLARED FOR VECTR
C
!      DIMENSION A(1),VALU(1),VALL(225),UPPERD(225),DIAG(225),V(2250),
!     $SQGAMP(1),INTER(1),TMULT(1),T(225,3),SQGAMR(1),UPPER2(450)
!      DIMENSION VECTR(MDIM,1),IFANA(225)
!      EQUIVALENCE (TMULT,INTER,SQGAMP,VALL),(T,SQGAMR),(T(1,2),UPPER2)
!       EQUIVALENCE (TMULTP,IMULTP),(COMPL1,ICMPL1)


      DIMENSION A(0:25425),VALU(1),
	1  VALL(225),UPPERD(225),DIAG(225),V(0:2250),
     1  SQGAMP(225),INTER(225),TMULT(225),T(225,3)
     1 ,SQGAMR(0:255),UPPER2(0:450)
      DIMENSION VECTR(MDIM,1),IFANA(225)
      EQUIVALENCE (TMULT,INTER,SQGAMP,VALL)
	equivalence (T,SQGAMR(1)),(T(1,2),UPPER2(1))
       EQUIVALENCE (TMULTP,IMULTP),(COMPL1,ICMPL1)
	v(0)=0.
	sqgamr(0)=0.
	upper2(0)=0.
	A(0)=0.


      NZ=0
      N=NN
      IF(N.LE.1) GO TO 41
      NP1=N+1
      NM1=N-1
      NM2=N-2
      NT2P1=N*2+1
      IX=0
      DO 10 I=1,NM2
      SIGMA2=0.
      IP1=I+1
      DO 1 J=IP1,N
      IJ=IX+J
    1 SIGMA2=SIGMA2+A(IJ)**2
      SIGMA=SQRT(SIGMA2)
      II=IX+I
      DIAG(I)=A(II)
      IIP1=IX+I+1
      UPPERD(I)=-SIGN(SIGMA,A(IIP1))
      UPPER2(I)=SIGMA2
      IF(ABS(SIGMA).GT.ABS(A(IIP1)))GO TO 2
      UPPERD(I)=A(IIP1)
      A(IIP1)=0.
      GO TO 10
    2 A(IIP1)=SQRT(1.+ABS(A(IIP1))/SIGMA)
      SQTGAM=-SIGN(SIGMA*A(IIP1),UPPERD(I))
      IP2=I+2
      DO 3 J=IP2,N
      IJ=IX+J
    3 A(IJ)=A(IJ)/SQTGAM
      JK1=I*(2*N-I-1)/2
      JX=JK1
      IIX=JK1
      DO 5 J=IP1,N
      SQGAMP(J)=0.
      JK=JK1+J
      DO 4 K=IP1,J
      IK=IX+K
      SQGAMP(J)=SQGAMP(J)+A(JK)*A(IK)
    4 JK=JK+N-K
      IF(J.EQ.N)GO TO 6
      CALL LOOP1(J+2,NP1,SQGAMP(J),A(JX),A(IX))
    5 JX=JX+N-J
    6 DELGAM=0.
      DO 7 J=IP1,N
      IJ=IX+J
    7 DELGAM=DELGAM+A(IJ)*SQGAMP(J)
      DGO2=.5*DELGAM
      DO 8 J=IP1,N
      IJ=IX+J
    8 SQGAMR(J)=SQGAMP(J)-DGO2*A(IJ)
      DO 9 II=IP1,N
      III=IX+II
      CALL LOOP2(A(IIX),A(IX),SQGAMR(NZ),SQGAMR(II),A(III),II+1,NP1)
    9 IIX=IIX+N-II
   10 IX=IX+N-I
      M=N*(N+1)/2
      UPPERD(NM1)=A(M-1)
      UPPER2(NM1)=UPPERD(NM1)**2
      DIAG(NM1)=A(M-2)
      DIAG(N)=A(M)
      ENORM=AMAX1(ABS(DIAG(1))+ABS(UPPERD(1)),ABS(DIAG(N))+ABS(UPPERD(NM
     11)))
      DO 11 I=2,NM1
      ENRTMP=ABS(DIAG(I))+ABS(UPPERD(I))+ABS(UPPERD(I-1))
   11 IF(ENRTMP.GT.ENORM)ENORM=ENRTMP
      DO 12 I=1,NEIG
      VALU(I)=ENORM
   12 VALL(I)=-ENORM
      DO 24 I=1,NEIG
   13 ROOT=.5*(VALU(I)+VALL(I))
      IF(ROOT.EQ.VALL(I).OR.ROOT.EQ.VALU(I))GO TO 24
      NAGREE=0
      PM2=0.
      PM1=1.
      DO 21 J=1,N
      IF(PM2.NE.0.)GO TO 15
   14 PM1=SIGN(1.,PM1)
      GO TO 17
   15 IF(PM1.NE.0.)GO TO 17
   16 P=-SIGN(1.,PM2)
      PM2=0.
      IF(UPPER2(J-1))18,14,18
   17 P=DIAG(J)-ROOT-UPPER2(J-1)*PM2/PM1
      PM2=1.
   18 IF(P)21,19,20
   19 PM2=PM1
      IF(PM2)21,20,20
   20 NAGREE=NAGREE+1
   21 PM1=P
      DO 23 J=I,NEIG
      IF(J.LE.NAGREE)GO TO 22
      IF(VALU(J).LE.ROOT)GO TO 13
      VALU(J)=ROOT
      GO TO 23
   22 VALL(J)=ROOT
   23 CONTINUE
      GO TO 13
   24 CONTINUE
      IF(NVEC.EQ.0)GO TO 39
      EPSLON=ENORM*1.E-6
       ICMPL1=NOT(1)
      DO 38 I=1,NVEC
      DO 25 J=1,N
      V(J)=1.
      T(J,2)=DIAG(J)-VALU(I)
      IF(J.EQ.N)GO TO 26
      T(J,3)=UPPERD(J)
   25 T(J+1,1)=UPPERD(J)
   26 T(N,3)=0.
      DO 29 J=1,N
      IF(ABS(T(J,2)).LT.1.E-15)T(J,2)=EPSLON
      T(J,1)=T(J,2)
      T(J,2)=T(J,3)
      T(J,3)=0.
      IF(J.EQ.N)GO TO 30
      INTER(J)=0
       IFANA(J)=0
      JP1=J+1
      IF(ABS(T(JP1,1)).LE.ABS(T(J,1)))GO TO 28
CO    INTER(J)=1
       IFANA(J)=1
      DO 27 K=1,3
      TEMP=T(J,K)
      T(J,K)=T(JP1,K)
   27 T(JP1,K)=TEMP
   28 TMULTP=T(JP1,1)/T(J,1)
CO    TMULT(J)=IOR(INTER(J),IAND(IMULTP,ICMPL1))
       TMULT(J)=TMULTP
      T(JP1,2)=T(JP1,2)-TMULTP*T(J,2)
   29 T(JP1,3)=T(JP1,3)-TMULTP*T(J,3)
   30 ITER=1
   31 DO 32 J=1,N
      L=N+1-J
      V(L)=(V(L)-T(L,2)*V(L+1)-T(L,3)*V(L+2))/T(L,1)
  32   CONTINUE
      VNORM=0.
      DO 33 L=1,N
   33 VNORM=VNORM+V(L)**2
      VNORM=SQRT(VNORM)
      DO 34 J=1,N
   34 V(J)=V(J)/VNORM
      IF(ITER.EQ.2)GO TO 36
      ITER=2
      DO 35 L=2,N
      LM1=L-1
CO    IF(IAND(INTER(LM1),1).EQ.0.)GO TO 35
       IF(IFANA(LM1).EQ.0)GO TO 35
      VTEMP=V(LM1)
      V(LM1)=V(L)
      V(L)=VTEMP
   35 V(L)=V(L)-TMULT(LM1)*V(LM1)
      GO TO 31
   36 IF(VNORM.EQ.0.)V(I)=1.
      IIX=(N*N-N-6)/2
      DO 37 KK=1,NM2
      IIP1=N-KK
      UTV=0.
      CALL LOOP3(UTV,A(IIX),V(NZ),IIP1+1,NP1)
      CALL LOOP4(A(IIX),V(NZ),NP1,IIP1+1,UTV)
   37 IIX=IIX+IIP1-N-2
      DO 100 J=1,N
 100  VECTR(J,I)=V(J)
   38 CONTINUE
      GO TO 39
   41 VECTR(1,1) = 1.0
      VALU(1) = A(1)
   39 RETURN
       END
C
C#######################################################################
C***** SUBROUTINES FOR BIGMAX: LOOP1, LOOP2, LOOP3 AND LOOP4  **********
C
      SUBROUTINE LOOP1(JP2,NP1,SGAMPJ,AJX,AIX)
      DIMENSION AJX(1), AIX(1)
      DO 1 L=JP2,NP1
    1 SGAMPJ=SGAMPJ+AJX(L)*AIX(L)
      RETURN
       END
      SUBROUTINE LOOP2(AIIX,AIX,S,SI,AIII,IP1,NP1)
      DIMENSION AIIX(1),AIX(1),S(1)
      DO 2 JJ=IP1,NP1
    2 AIIX(JJ)=AIIX(JJ)-AIII*S(JJ)-SI*AIX(JJ)
      RETURN
       END
      SUBROUTINE LOOP3(UTV,AIIX,V,IIP2,NP1)
      DIMENSION AIIX(1), V(1)
      DO 3 J=IIP2,NP1
    3 UTV=UTV+AIIX(J)*V(J)
      RETURN
       END
      SUBROUTINE LOOP4(AIIX,V,NP1,IIP2,UTV)
      DIMENSION AIIX(1),V(1)
      DO 4 K=IIP2,NP1
    4 V(K)=V(K)-AIIX(K)*UTV
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE HELP
C
C     NOTE THAT F HAS BEEN ENLARGED TO ALLOW HIGH SPIN CALCNS
C     (IN CLEBI). TESTED BY PBS
C
      COMMON /LU/ LU
      COMMON /PTR/ F(135)
C*****COMPUTATION OF F(N+1)=LN(N^)
      F(1) = 0.
      DO 101 I=1,134
      AII = FLOAT(I)
      F(I+1) = F(I) + ALOG(AII)
  101 CONTINUE
      RETURN
      END
C
C#######################################################################
C
      FUNCTION CLEBI(I1,I2,I3,N1,N2,N3)
C
CALCULATION OF CLEBSCH-GORDAN COEFF <I1/2 N1/2 I2/2 N2/2 ^ I3/2 N3/2>
C
C     NOTE THAT FCT HAS BEEN INCREASED TO ALLOW CALCULATIONS TO HIGH
C     SPIN. TESTED BY PBS FOR I1.LE.121, I2=2 OR 4, N1.LE.13

      COMMON /PTR/ FCT(135)
      INTEGER Z,ZMIN,ZMAX
      J1=I1
      J2=I2
      J =I3
      N=135
      M1=N1
      M2=N2
      M=-N3
      CC=0.0
      JSUM=J1+J2+J
      JM1 =J1-IABS(M1)
      JM2 =J2-IABS(M2)
      JM3 =J -IABS(M )
      IF((MOD(JSUM,2).NE.0).OR.(MOD(JM1,2).NE.0).OR.(MOD(JM2,2).NE.0)
     1.OR.(MOD(JM3,2).NE.0).OR.(JM1.LT.0).OR.(JM2.LT.0).OR.(JM3.LT.0))
     2GO TO 1
      IF((M1+M2+M.NE.0).OR.(J.GT.J1+J2).OR.(J.LT.IABS(J1-J2))) GO TO 1
      ZMIN=0
      IF(J-J2+M1.LT.0) ZMIN=-J+J2-M1
      IF(J-J1-M2+ZMIN.LT.0) ZMIN=-J+J1+M2
      ZMAX=J1+J2-J
      IF(J2+M2-ZMAX.LT.0) ZMAX=J2+M2
      IF(J1-M1-ZMAX.LT.0) ZMAX=J1-M1
      JA=(J1+M1)/2+1
      JB=JA-M1
      JC=(J2+M2)/2+1
      JD=JC-M2
      JE=(J +M )/2+1
      JF=JE-M
      JG=(J1+J2-J)/2+1
      JH=JA+JB-JG
      JI=JC+JD-JG
      JJ=JE+JF+JG-1
      IF(JJ.GT.N) THEN
        write(6, 12) JJ
        GO TO 5
      ENDIF
   12 FORMAT(/5X,' ERROR IN CLEBI: JJ= ',I4)
      IA=ZMIN/2
      IB=JG-IA+1
      IC=JB-IA+1
      ID=JC-IA+1
      IE=JA-JG+IA
      IF=JD-JG+IA
      FASE=1.0
      IF(MOD(IA,2).EQ.0) FASE=-FASE
      Z =ZMIN
      ARII=FLOAT(J+1)
      ARI=(FCT(JA)+FCT(JB)+FCT(JC)+FCT(JD)+FCT(JE)+FCT(JF)+FCT(JG)
     1 +FCT(JH)+FCT(JI)-FCT(JJ)+ALOG(ARII))/2.0
    2 IA=IA+1
      IB=IB-1
      IC=IC-1
      ID=ID-1
      IE=IE+1
      IF=IF+1
      FASE=-FASE
      ARP=-FCT(IA)-FCT(IB)-FCT(IC)-FCT(ID)-FCT(IE)-FCT(IF)
      ARS=ARI+ARP
      IF(ARS.GT.80.0) THEN
       write(6, 14) ARS
       GO TO 5
      ENDIF
   14 FORMAT(/5X,' ERROR IN CLEBI: ARS= ', F9.4)
      IF(ARS.LT.-80.0) GO TO 10
      CC=CC+FASE*EXP(ARS)
   10 CONTINUE
      Z=Z+2
      IF(Z.LE.ZMAX) GO TO 2
    1 CLEBI=CC
      RETURN
C   3 WRITE(6,4)
C   4 FORMAT(26H0 ERROR IN ARGUMENT OF VCC)
C     GO TO 7
    5 write(6, 6)
      STOP
    6 FORMAT(29H0 ERROR - FACTORIALS EXCEEDED)
C   7 WRITE(61,8)I1,I2,I3,N1,N2,N3
C   8 FORMAT(10I10)
C     RETURN
      END
CV
CV######################################################################
CV
      SUBROUTINE INERTIA(IR,AOO,STIFF,AOR)
      R=IR*(IR+1)/4./STIFF+1/216./AOO**3
      T=R+1./216./AOO**3
      T=T*IR*(IR+1)/4./STIFF
      T=T**.5
c	write(60,*)ir,r+t,r-t,aoo
c	can get rounding error and (R-T) accidentally going negative
c	aes 30.3.93
c
	R_T=R-T
	if((R_T.lt.0.0).and.(R_T.gt.-0.01))R_T=0.0
      AOR=(R+T)**(1./3.)+(R_T)**(1./3.)+1./6./AOO
      AOR=1./2./AOR
      RETURN
      END
CV
CV######################################################################
CV
