C#######################################################################
C##########                                           ##################
C##########       PROGRAM GAMPN                       ##################
C##########                                           ##################
C##########       VERSION: JULY, 1991                 ##################
C##########                                           ##################
C##########     RECEIVED AT ANU FROM PAUL SEMMES      ##################
C##########       4-MAR-1993   A.E. STUCHBERY         ##################
C##########                                           ##################
C##########        PC version AES May 2007            ##################
C######      I have replaced time with elapsed time       ##############
C##########                                           ##################
C##########                                           ##################
C#######################################################################
C***    ORIGINAL PROGRAM DESCRIBED IN NUCL PHYS A307 (1978) 189     ****
C                        AND PHYSICA SCRIPTA 8 (1973) 17               *
C                                                                      *
C  PLEASE REPORT ERRORS, INCONSISTENCIES OR PROBLEMS TO:               *
C                                                                      *
C      INGEMAR RAGNARSSON                     PAUL SEMMES              *
C      DEPT OF MATH PHYSICS        OR         PHYSICS DEPT, BOX 5051   *
C      LTH, BOX 118                           TENESSEE TECH. UNIV.     *
C      S-22100 LUND, SWEDEN                   COOKEVILLE, TN 38505 USA *
C                                                                      *
C                                                                      *
C***********************************************************************
C
C     VERSION: JULY, 1991
C
C     DISTRIBUTED AT THE OAK RIDGE THEORY WORKSHOP, AUG 1991, FOR USE
C     WITH OTHER MODIFIED-OSCILLATOR CODES ASYRMO (DIAGONALIZES THE
C     PARTICLE + TRIAXIAL ROTOR HAMILTONIAN) AND PROBAMO (CALCULATES
C     M1/E2 MATRIX ELEMENTS FROM THE ENERGY EIGENVECTORS OF ASYRMO)
C
C
C     GAMPN CALCULATES OF SINGLE-PARTICLE ENERGIES AND VARIOUS MATRIX
C     ELEMENTS IN A MODIFIED OSCILLATOR (NILSSON) POTENTIAL.
C     PROGRAM MAINLY DESIGNED TO PROVIDE INPUT FOR SUBSEQUENT ODD-A
C     AND ODD-ODD PARTICLE-ROTOR CALCULATIONS
C
C     SINGLE-PARTICLE TRANSFORMATION:  /N l j omega> basis   --->
C      /N l lambda sigma> basis WRITTEN ON FILE2 (NEEDED FOR ODD-ODD
C      CALCULATIONS WITH A PROTON-NEUTRON RESIDUAL INTERACTION)
C
C     DATA FOR VMI CALCULATIONS (IN ASYRMO/ASYRPN) WRITTEN ON FILE16
C     DATA FOR PART-ROTOR PROGRAMS (ASYRPN AND PROBA) WRITTEN ON FILE17
C
C     ALSO POSSIBLE TO CALCULATE ENERGIES IN A CRANKED POTENTIAL
C                                              (OMROT .NE. 0)
C
C***********************************************************************
C
C  DESCRIPTION OF THE INPUT DATA (FREE FORMAT UNLESS SPECIFIED):
C
C  CARD 1. FILE2, FILE16, FILE17
C           OUTPUT FILES TO BE USED FOR SUBSEQUENT APPLICATIONS
C            (RECALL THAT ON A VAX FREE FORMAT CHARACTER
C             DATA MUST BE ENCLOSED WITH APOSTROPHES)
C           FILE2 CONTAINS THE S.P. WAVE FUNCTIONS IN THE
C                 |N,L,LAMBDA,SIGMA> BASIS NEEDED FOR Vpn CALCULATIONS
C           FILE16 CONTAINS THE S.P. WAVEFUNCTIONS NEEDED FOR THE VMI
C                 OPTION IN ASYR
C           FILE17 CONTAINS THE USUAL OUTPUT NEEDED FOR ASYR, IE,
C                  SINGLE PARTICLE MATRIX ELEMENTS LIKE <j+>, <j->, etc
C
C       2. ISTRCH, ICORR, IREC
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                        (NOTE THAT THE SHAPES GENERATED ARE DIFFE-
C                         RENT IN THE TWO CASES; FURTHERMORE THE
C                         L*S AND L*L TERMS ARE SOMEWHAT DIFFERENT
C                         WITH L DEF IN STRETCHED OR SPHER COORDINATES)
C
C           ICORR=0:  NORMAL CALCULATION OF J+/J-/JZ/J+J_, ETC, IE,
C                     APPROXIMATE WITH J+/J-/JZ (STRETCHED)
C           ICORR=1:  INCLUDE CORRECTION TERMS FOR J+/J-/JZ FROM THE
C                     STRETCHED BASIS (SEE APPENDIX, NUCL PHYS A307,189)
C                     AND CONSTRUCT CORRECTED J+J-, JZ2, ETC FROM THESE
C
C           IREC = 0: USUAL TREATMENT OF JX2, JY2, JZ2 IN ASYR, IE, AS
C                     ONE-BODY OPERATORS WITH UU' + VV' PAIRING FACTORS
C           IREC = 1: "RECOIL" TERMS JX2, JY2, JZ2 ARE TREATED AS TWO-
C                     BODY OPERATORS IN ASYRMO => MORE j+, j_, jz
C                     MATRIX ELEMENTS MUST BE CALCULATED HERE (FOR THE
C                     SUMS OVER INTERMEDIATE STATES IN ASYRMO)
C
C  ** RECOMMENDATIONS: FOR NORMAL DEFORMATIONS, THE DIFFERENCES BETWEEN
C  ***                 THE STRETCHED AND PHYSICAL j+, j_, jz OPERATORS
C  ***                 ARE NEGLIGIBLE ==> ICORR=0 IS STANDARD.
C  ***                    SIMILARLY, THE SIMPLER ONE-BODY TREATMENT OF
C  ***                 THE RECOIL TERM SEEMS QUITE ADEQUATE ==> IREC=0
C  ***                 IS STANDARD. (NOTE THAT IF IREC=1 IS SELECTED
C  ***                 HERE, IREC=0 CAN STILL BE CHOSEN IN ASYRMO.)
C
C       3. NKAMY, NNEUPR
C           NKAMY .LE. 0: SAME KAPPA AND MY FOR ALL N-SHELLS
C           NKAMY .GT. 1: DIFFERENT KAPPA AND MY FOR DIFFERENT N-SHELLS
C                                          (SEE CARD 3 BELOW)
C           FOR NKAMY .GT. 1:(NKAMY-1) .GE. MAX(NPROT,NNEUTR); (CARD 9)
C                             (NKAMY-1) .LE. 10
C           NNEUPR=1:  CALCULATION FOR PROTONS
C           NNEUPR=-1: CALCULATION FOR NEUTRONS
C           NNEUPR=2:  CALCULATION FOR PROTONS AND NEUTRONS (FOR ODD-ODD)
C
C       4. KAPPAP(I), MYP(I), KAPPAN(I), MYN(I)
C           OMITTED IF NKAMY .LE. 1; OTHERWISE NKAMY CARDS NEEDED
C           KAPPAP AND MYP (KAPPAN AND MYN) VALUES OF L.S AND L**2
C           COUPLING PARAMETERS FOR PROTONS (NEUTRONS) FOR THE
C           DIFFERENT SHELLS, N=0,1,2,3,..., (NKAMY-1)
C
C          4A. NKAMYL
C                NUMBER OF N-SHELLS WITH KAPPA AND MY L-DEPENDENT
C
C            CARDS 3B-3C REPEATED NKAMYL TIMES (OMIT FOR NKAMYL=0):
C              4B. NSHELL
C                    N QUANTUM NUMBER FOR SHELLS WITH L-DEPENDENCE
C              4C. KAPPAP(N,L), MYP(N,L), KAPPAN(N,L), MYN(N,L)
C                    L=0,2,.. OR L=1,3,..        (N=NSHELL)
C
C  *** FOR "STANDARD" KAPPA'S & MU 'S, SEE NUCL PHYS A436 (1985) 14 ***
C
C       5. KAPPAP, MYP, KAPPAN, MYN
C           THIS CARD MAY NEVER BE OMITTED BUT IS DUMMY FOR NKAMY .GT. 1
C           FOR NKAMY .LE. 1 USED AS CARDS 3 BUT REFERRING TO ALL SHELLS
C
C       6. NUU, IPKT, NOYES, ITRANS
C           NUU: NUMBER OF OSCILLATOR SHELLS COUPLED, NUU=4 STAND VALUE
C                  (REDUNDANT FOR EPS4=EPS6=OMROT=0 AND ISTRCH=1 BECAUSE
C                   THE N-SHELLS ARE PURE IN THIS CASE)
C           IPKT: NUMBER OF DEFORMATIONS, CF. CARD 9 BELOW
C           NOYES=1 IN PART-ROT CALCULATIONS
C                  (FOR NOYES=0 ONLY SINGLE PARTICLE ENERGIES BUT NO
C                   WAVE-FUNCTIONS OR MATRIX ELEMENTS ARE CALCULATED)
C           ITRANS: CONTROLS OPTION OF TRANFORMING THE LEVELP(J) AND
C                   LEVELN(K) PROTON AND NEUTRON ORBITALS TO THE
C                   /N,L,LAM,SIGMA> BASIS (NEEDED IF SUBSEQUENT Vpn
C                   CALCULATIONS ARE DESIRED, OUTPUT TO FILE2)
C           ITRANS=0: NO SUCH TRANSFOMATION
C           ITRANS=1: WAVEFUNCTIONS ARE TRANSFORMED, BUT NOT PRINTED
C           ITRANS=2: W.F.'S ARE TRANSFORMED AND PRINTED OUT
C
C       7. EMIN, EMAX
C           FOR ENERGIES BETWEEN EMIN AND EMAX (IN OSCILLATOR UNITS) THE
C           WAVE-FUNCTIONS ARE PRINTED
C
C       8a. IPARP, NORBP, (LEVELP(J),J=1,NORBP)      FORMAT(A1,I2,15I3)
C           IPARP: PARITY CONSIDERED FOR PROTONS, "+" OR "-"
C           NORBP: NUMBER OF PROTON ORBITALS CONSIDERED IN THE
C               PARTICLE-ROTOR CALCULATION, NORBP .LE. 15
C           LEVELP(J): "NUMBERS" ON THE ORBITALS INCLUDED
C             (CALCULATED  FROM THE "BOTTOM" AND FOR THE PARITY "IPARP")
C
C       8b. IPARN,NORBN,(LEVELN(J),J=1,NORBN)
C                      ANALOGOUS TO 7a BUT FOR NEUTRONS
C
C       9. Z,A
C           Z: PROTON NUMBER;   A: MASS NUMBER
C           FOR Z AND/OR N=A-Z PARTICLES (Z AND A EVEN), THE MICROSCOPIC
C           QUADRUPOLE MOMENTS Q0 AND Q2 ARE CALCULATED (WITH A SHARP
C           FERMI SURFACE). AN ASYMMETRY PARAMETER GAMMA CORRESPONDING
C           TO THE MATTER DISTRIBUTION IS THEN DEDUCED
C
C      10. IPKT CARDS:
C         EPS2, GAMMA, EPS4, EPS6, OMROT, NPROT, NNEUTR, NSHELP, NSHELN
C                        (ISTRCH = 1)
C                              OR
C         DELTA2,GAMMAD,DELTA4,DELTA6,OMROT,NPROT,NNEUTR,NSHELP,NSHELN
C                        (ISTRCH = 0)
C             FOUR DEFORMATION PARAMETERS
C                (EPS6 (DELTA6) IS NOT DEFINED TO HAVE THE CORRECT
C                TRANSFORMATION PROPERTIES FOR GAMMA .NE. 0 AND SHOULD
C                ONLY BE USED FOR SMALL GAMMA, SAY GAMMA .LE. 15)
C             OMROT: ROTATIONAL FREQUENCY (=0 IN PART-ROTOR CALC)
C             NPROT: NUMBER OF OSCILLATOR SHELLS INCLUDED FOR PROTONS
C             NNEUTR: NUMBER OF OSCILLATOR SHELLS INCLUDED FOR NEUTRONS
C                MAX NUMBER OF SHELLS IN PRESENT VERSION: 10
C             NSHELP,NSHELN:USED IN THE TRANSFORMATION TO /N,L,LAM,SIG>
C                BASIS IF THE N-SHELLS ARE PURE. NSHELP IS THEN THE
C                N-SHELL FOR THE LEVELP(J) PROTON ORBITALS, NSHELN
C                FOR THE LEVELN(K) NEUTRON ORBITALS
C
C***********************************************************************
C
      COMMON /COUPL/ NUU
      COMMON /LU/ LU
      COMMON /PI/ PI
      COMMON /CMN1/ KAPPA,OMOM,MY,L2F(18)
      COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22
      COMMON /CMN7/ EPSP,GAMMA
      COMMON /QSTR/ Q
      COMMON/E/ EMIN,EMAX
      COMMON/LEVEL/ NU,LEVEL(15)
      COMMON /NUCL/ Z,A
      COMMON /PAR/ IPAR
      COMMON/OMROT/ OMROT
      COMMON/EQFAC/ EFAC,QFAC,ISO
      COMMON/ITRANS/ ITRANS,NSHELL
      DIMENSION LEVELP(15),LEVELN(15)
      COMMON/BAS/ ISTRCH,ICORR,IREC
      COMMON/XKAMY/ XKAPPA(13,7),XMY(13,7),AVKAMY(13)
      INTEGER Z,A
      REAL KAPPA,MY,KAPPAP,MYP,KAPPAN,MYN
      REAL XKAP(13),XMYP(13),XKAN(13),XMYN(13),
     $     XKAPL(13,7),XMYPL(13,7),XKANL(13,7),XMYNL(13,7),
     $     AVKMYP(13),AVKMYN(13)
      CHARACTER IPAR*4,IPARP*4,IPARN*4,FILE2*24,FILE16*24,FILE17*24

      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='GAMPN.out',status='unknown',form='formatted')


      LU=6
      PI=3.1415926536

      READ(5,*)FILE2,FILE16,FILE17
C
C  OPEN FILES  (APPROPRIATE FOR VAX)
C
      OPEN(2,FILE=FILE2,  FORM='FORMATTED',  STATUS='unknown')
      OPEN(16,FILE=FILE16,FORM='UNFORMATTED',STATUS='unknown')
      OPEN(17,FILE=FILE17,FORM='UNFORMATTED',STATUS='unknown')

!      OPEN(2,FILE=FILE2,  FORM='FORMATTED',  STATUS='NEW')
!      OPEN(16,FILE=FILE16,FORM='UNFORMATTED',STATUS='NEW')
!      OPEN(17,FILE=FILE17,FORM='UNFORMATTED',STATUS='NEW')
C
C
CALCULATIN OF FACTORIALS, GAMMA FUNCTION AND
C     AVERAGE WITHIN SHELLS OF L*L
C
      CALL HELP
C
      EPS1=0.0
      EPS3=0.0
      EPS5=0.0
      READ(5,*) ISTRCH,ICORR,IREC
      READ(5,*) NKAMY,NNEUPR
      IF (NKAMY .LE. 1) GO TO 51
      write(6,16)
   16 FORMAT(///,'         KAPPAP      MYP     KAPPAN      MYN')
      DO 52 IN=1,NKAMY
      READ(5,*)XKAP(IN),XMYP(IN),XKAN(IN),XMYN(IN)
      NSHELL=IN-1
      write(6,15) NSHELL,XKAP(IN),XMYP(IN),XKAN(IN),XMYN(IN)
   15 FORMAT(I5,4F10.4)
      NL = (IN+1)/2
      DO 57 IL = 1,NL
         XKAPL(IN,IL) = XKAP(IN)
         XKANL(IN,IL) = XKAN(IN)
         XMYPL(IN,IL) = XMYP(IN)
         XMYNL(IN,IL) = XMYN(IN)
         AVKMYP(IN) = XMYP(IN)*XKAP(IN)
         AVKMYN(IN) = XMYN(IN)*XKAN(IN)
   57 CONTINUE
   52 CONTINUE
      READ(5,*) NKAMYL
      IF (NKAMYL .EQ. 0) GO TO 61
      write(6,18) '   SHELLS WITH KAPPA L-DEPENDENT',
     $' N     L    KAPPAP    MYP       KAPPAN       MYN'
   18 FORMAT(/,A,/,TR19,A,/)
      DO 63 I = 1,NKAMYL
         READ(5,*) NSHELL
         IN = NSHELL+1
         NL = NSHELL/2 + 1
         SUMMYP = 0.
         SUMMYN = 0.
         SUMORB = 0.
         DO 62 IL = 1,NL
            READ(5,*)XKAPL(IN,IL),XMYPL(IN,IL),XKANL(IN,IL),XMYNL(IN,IL)
            IF (MOD(NSHELL,2) .EQ. 0) THEN
               L = 2*(IL-1)
            ELSE
               L = 2*IL-1
            ENDIF
            SUMMYP = SUMMYP + L*(L+1)*(2*L+1)*XMYPL(IN,IL)*XKAPL(IN,IL)
            SUMMYN = SUMMYN + L*(L+1)*(2*L+1)*XMYNL(IN,IL)*XKANL(IN,IL)
            SUMORB = SUMORB + L*(L+1)*(2*L+1)
            write(6,17) NSHELL,L,
     $      XKAPL(IN,IL),XMYPL(IN,IL),XKANL(IN,IL),XMYNL(IN,IL)
   17       FORMAT(15X,2I6,4F10.4)
   62    CONTINUE
         AVKMYP(IN) = SUMMYP/SUMORB
         AVKMYN(IN) = SUMMYN/SUMORB
   63 CONTINUE
   61 CONTINUE
   51 CONTINUE
      READ(5,*)KAPPAP,MYP,KAPPAN,MYN
      READ(5,*)NUU,IPKT,NOYES,ITRANS
      READ(5,*)EMIN,EMAX
      READ(5,1952) IPARP,NORBP,(LEVELP(J),J=1,NORBP)
      READ(5,1952) IPARN,NORBN,(LEVELN(J),J=1,NORBN)
 1952 FORMAT(A1,I2,15I3)
      READ(5,*)Z,A
      N=A-Z
      DO 1 I=1,IPKT
!      IRK=MSLEFT(XX)




      READ(5,*)EPSP,GAMMA,EPS4,EPS6,OMROT,NPROT,NNEUTR,NSHELP,NSHELN
      IF(NPROT.LT.3.OR.NNEUTR.LT.3) GO TO 111
      IF(NPROT.GT.12.OR.NNEUTR.GT.12) GO TO 111
      GAMMAR=GAMMA*PI/180.0
      EPS2=EPSP*COS(GAMMAR)
      EPS22=EPSP*SIN(GAMMAR)/SQRT(0.75)
C
C     CALCULATION OF VOLUME CONSERVATION FACTOR, OMOM
C
      CALL OMO(OMOM)
C
C*************SINGLE PARTICLE ENERGY LEVELS CALCULATION
      IF (I .NE. 1) write(6,12)
   12 FORMAT(1H1)
      IF (NNEUPR .LT. 0) GO TO 71
      IF (ISTRCH .EQ. 1) THEN
         write(6,3)
      ELSEIF (ISTRCH .EQ. 0) THEN
         write(6,13)
      ENDIF
    3 FORMAT(//,'     KAPPA    MY     EPS   GAMMA    EPS4   ',
     *'  EPS6     W0/W00   NMAX  COUPL     OMROT      EFAC      QFAC')
   13 FORMAT(//,'     KAPPA    MY    DELTA  GAMMAD  DELTA4   ',
     *'DELTA6    W0/W00   NMAX  COUPL     OMROT      EFAC      QFAC')
      ISO=1
      IPAR = IPARP
      KAPPA=KAPPAP
      MY=MYP
      NPROT1 = NPROT + 1
      DO 53 IN=1,NPROT1
      NL = (IN-1)/2+1
      DO 58 IL=1,NL
      IF (NKAMY .GE. 2) THEN
         XKAPPA(IN,IL) = XKAPL(IN,IL)
         XMY(IN,IL) = XMYPL(IN,IL)
         AVKAMY(IN) = AVKMYP(IN)
      ELSE
         XKAPPA(IN,IL) = KAPPAP
         XMY(IN,IL) = MYP
         AVKAMY(IN) = MYP*KAPPAP
      ENDIF
   58 CONTINUE
   53 CONTINUE
C********FACTORS FOR TRANSFORMATION FROM OSC UNITS TO PHYSICAL UNITS
         EFAC=41.0/FLOAT(A)**(1./3.)*(1.-FLOAT(N-Z)/A/3.)*OMOM
         QFAC=41.5/EFAC
      NU = NORBP
      DO 56 IP=1,NU
         LEVEL(IP) = LEVELP(IP)
   56 CONTINUE
      NSHELL=NSHELP
      IF (NKAMY .LE. 0) write(6,100) KAPPA,MY,
     *EPSP,GAMMA,EPS4,EPS6, OMOM,NPROT,NUU,OMROT,EFAC,QFAC
  100 FORMAT(2X,F7.4,2F8.3,F7.1,F8.3,F9.3,F11.5,2I6,F12.5,2F10.5,/)
      IF (NKAMY .GT. 1) write(6,101) '    SEE TABLE  ',
     *EPSP,GAMMA,EPS4,EPS6, OMOM,NPROT,NUU,OMROT,EFAC,QFAC
  101 FORMAT(2X,A,F8.3,F7.1,F8.3,F9.3,F11.5,2I6,F12.5,2F10.5,/)
      IF (ISTRCH .EQ. 0) write(6,201)
      IF (ISTRCH .EQ. 1) write(6,202)
  201 FORMAT(/,20X,' **** PROTONS ****',20X,
     &             ' **** SPHERICAL COORDINATES ****')
  202 FORMAT(/,20X,' **** PROTONS ****',20X,
     &             ' **** STRETCHED COORDINATES ****')

C
CALCULATION OF ENERGY LEVELS, MATRIX ELEMENTS ETC + OUTPUT
C
      CALL ENERLJ(OMOM,NPROT,NOYES)
C
!      JRK=MSLEFT(XX)
      TID=(IRK-JRK)*0.001
      write(6,10) TID
   10 FORMAT(1X,16H     TOTAL TIME:,F6.1,4H SEC)
      IF (NNEUPR .EQ. 1) GO TO 1
   71 CONTINUE
C
CALCULATIONS FOR NEUTRONS
C
!      IRK=MSLEFT(XX)
      ISO=2
      IPAR = IPARN
      KAPPA = KAPPAN
      MY = MYN
      NNEUT1 = NNEUTR + 1
      DO 55 IN = 1,NNEUT1
      NL = (IN-1)/2+1
      DO 59 IL=1,NL
      IF (NKAMY .GT. 2) THEN
         XKAPPA(IN,IL) = XKANL(IN,IL)
         XMY(IN,IL) = XMYNL(IN,IL)
         AVKAMY(IN) = AVKMYN(IN)
      ELSE
         XKAPPA(IN,IL) = KAPPAN
         XMY(IN,IL) = MYN
         AVKAMY(IN) = MYN*KAPPAN
      ENDIF
   59 CONTINUE
   55 CONTINUE
      NU = NORBN
      DO 54 IN=1,NU
         LEVEL(IN) = LEVELN(IN)
   54 CONTINUE
      NSHELL=NSHELN
C********FACTORS FOR TRANSFORMATION FROM OSC UNITS TO PHYSICAL UNITS
         EFAC=41.0/FLOAT(A)**(1./3.)*(1.+FLOAT(N-Z)/A/3.)*OMOM
         QFAC=41.5/EFAC
      write(6,12)
      IF (ISTRCH .EQ. 1) THEN
         write(6,3)
      ELSEIF (ISTRCH .EQ. 0) THEN
         write(6,13)
      ENDIF
      IF (NKAMY .LE. 0) write(6,100) KAPPA,MY,
     *EPSP,GAMMA,EPS4,EPS6, OMOM,NNEUTR,NUU,OMROT,EFAC,QFAC
      IF (NKAMY .GT. 1) write(6,101) '    SEE TABLE  ',
     *EPSP,GAMMA,EPS4,EPS6, OMOM,NNEUTR,NUU,OMROT,EFAC,QFAC
      IF (ISTRCH .EQ. 0) write(6,203)
      IF (ISTRCH .EQ. 1) write(6,204)
  203 FORMAT(/,20X,' **** NEUTRONS ****',19X,
     &             ' **** SPHERICAL COORDINATES ****')
  204 FORMAT(/,20X,' **** NEUTRONS ****',19X,
     &             ' **** STRETCHED COORDINATES ****')
      CALL ENERLJ(OMOM,NNEUTR,NOYES)
C
c      JRK=MSLEFT(XX)
      TID=(IRK-JRK)*0.001
      write(6,10) TID
    1 CONTINUE
      GO TO 112
  111 CONTINUE
  113 FORMAT(1X,30HINPUT CARD INCORRECTLY PUNCHED)
      write(6,113)
      STOP
  112 CONTINUE
 1946 FORMAT(  )
      STOP
      END
C
C#######################################################################
C
      SUBROUTINE QNUMLJ(BB,NP,L,JJ,JOMEGA,SN)
C
C     SET UP OF BASIS FUNCTIONS (ONE FOR EACH CALL):
C     ON OUTPUT: NP=N, L=L, JJ=2*J, JOMEGA=2*OMEGA
C        INPUT: BB=NUMBERING OF WAVE-FUNCTIONS
C               SN=MAX NUMBER OF SHELLS (FOR GIVEN PARITY)
C
      INTEGER B,SN,ORDN,F,D,BB
      F=0
      B=BB
      OMEGA=SN+0.5
      IF((SN/2)*2.NE.SN) OMEGA=SN-0.5
      ORDN=SN-OMEGA+1.6
      NP=SN
      I1=ORDN
      I2=B-1
      IF(I1.GT.I2) GO TO 2
      D=I1
    1 CONTINUE
      F=F+ORDN
      OMEGAM=-NP+0.5
      IF((SN/2)*2.NE.SN) OMEGAM=-NP-0.5
      IF(ABS(OMEGA-OMEGAM).GT.0.1) GO TO 10
      NP=NP-2
      OMEGA=NP+0.5
      IF((NP/2)*2.NE.NP) OMEGA=NP-0.5
      GO TO 11
   10 CONTINUE
      OMEGA=OMEGA-2.0
   11 CONTINUE
      ORDN=NP-ABS(OMEGA)+1.6
      D=D+ORDN
      IF(D.LE.I2) GO TO 1
    2 CONTINUE
      B=B-F
      RJ=FLOAT(NP)+0.5
      IF((B/2)*2.EQ.B) GO TO 3
      S=-0.5
      GO TO 4
    3 CONTINUE
      S=0.5
    4 CONTINUE
      RJ=RJ+1.0-FLOAT(B)
      L=RJ+S+0.1
      JJ=2.*RJ+0.1
      JOMEGA=2.*OMEGA+0.1*SIGN(1.0,OMEGA)
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE QTME
C
CALCULATION OF QUADRUPOLE TRANSITION MATRIX ELEMENTS WITHIN THOSE
C     ORBITALS CONSIDERED IN PART-ROT CALCULATION
C***** QUADRUPOLE OPERATORS DEFINED AS:
C        Q20=SQRT(16*PI/5)*R2*Y20           (=2Z2-X2-Y2)
C        Q22=SQRT(4*PI/5)*R2*(Y22+Y2-2)     (=SQRT(3/2)*(X2-Y2))
C        Q20, Q22, R2 CALCULATED IN SPHERICAL BASIS (PHYSICAL VALUES)
C        Q20S, Q22S AND RHO2 SAME MATRIX ELEMENTS IN STRETCHED BASIS
C                                               (ISTRCH=1)
C        RJZ: MATRIX ELEMENTS OF JZ
C        KV1=N, KV2=L, KV3=2*J, KV4=2*OMEGA, KATT=PARITY ('+' OR '-')
C
      COMMON/CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22
      COMMON/DIM/ IDIM,KATT
      CHARACTER KATT*4
      COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161)
      COMMON/EIG/ EIGVEC(161,161),VALU(161)
      COMMON/BAS/ ISTRCH,ICORR,IREC
      COMMON/FACTOR/ GAC,HAC
      COMMON/QUADR/ Q20S(15,15),Q22S(15,15),Q20(15,15),Q22(15,15),
     *RHO2(15,15),R2(15,15),RJZ(15,15)
      COMMON/LEVEL/ NU,LEVEL(15)
      INTEGER RNP,RLP,RJP,OMEGAP,RN,RL,RJ,OMEGA
      DATA PI/3.1415926536/
      FAC=(1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)*(1.-2.*EPS2/3.)
      IF (NU .LE. 15) GO TO 4
      write(6,10)
   10 FORMAT(1H1,'QTME STOP')
      STOP
    4 CONTINUE
      DO 3 MYP=1,NU
      DO 3 MY =1,NU
      Q20S(MYP,MY)=0.0
      Q22S(MYP,MY)=0.0
      Q20(MYP,MY)=0.0
      Q22(MYP,MY)=0.0
      RHO2(MYP,MY)=0.0
      R2(MYP,MY)=0.0
      RJZ(MYP,MY)=0.0
    3 CONTINUE
      DO 2 IP=1,IDIM
      RNP=KV1(IP)
      RLP=KV2(IP)
      RJP=KV3(IP)
      OMEGAP=KV4(IP)
      DO 2 I=1,IDIM
      RN=KV1(I)
      RL=KV2(I)
      RJ=KV3(I)
      OMEGA=KV4(I)
C***************** IS THIS COMBINATION OF BASIS VECTORS SIGNIFICANT? ***
C     EIGMAX = 1.E-8
C     DO 25 MMY=1,NU
C     DO 24 NNY=MMY,NU
C        MY = LEVEL(MMY)
C        NY = LEVEL(NNY)
C        EIGABS = ABS(EIGVEC(IP,MY)*EIGVEC(I,NY))
C        EIGMAX = MAX(EIGABS,EIGMAX)
C  24 CONTINUE
C  25 CONTINUE
C     IF (EIGMAX .LT. 1.E-5) GO TO 2
C
C     IF ((ABS(EPS4)+ABS(EPS6)) .LT. 0.001 .AND. RNP .NE. RN .AND.
C    *  (ISTRCH .EQ. 1)) GO TO 2
      IF(IABS(RNP-RN).GT.2.OR.IABS(RJP-RJ).GT.4.OR.IABS(RLP-RL).GT.2.
     $OR.IABS(OMEGAP-OMEGA).GT.4) GO TO 6
      CGC=0.0
      CGC1=0.0
      CGC2=0.0
      RME=0.0
      HLZ=0.0
      HJZ=0.0
C
CALCULATION OF RADIAL MATRIX ELEMENTS, R2
C
      RME=RMATEL(2*RNP,2*RLP,2*RN,2*RL,2)
C
      IF(OMEGAP.EQ.OMEGA) CGC=CLEBI(RJ,4,RJP,-1,0,-1)*
     $CLEBI(RJ,4,RJP,OMEGA,0,OMEGA)
      IF(OMEGAP.EQ.OMEGA+4) CGC1=CLEBI(RJ,4,RJP,-1,0,-1)*
     $CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)
      IF(OMEGAP.EQ.OMEGA-4) CGC2=CLEBI(RJ,4,RJP,-1,0,-1)*
     $CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP)
      IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA)
     $HJZ=FLOAT(OMEGA)/2.
      IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.OMEGAP.EQ.OMEGA) HLZ=
     $(OMEGA-1.)/2.*CLEBI(2*RL,1,RJP,OMEGA-1,1,OMEGA)*
     $CLEBI(2*RL,1,RJ,OMEGA-1,1,OMEGA)+
     $(OMEGA+1.)/2.*CLEBI(2*RL,1,RJP,OMEGA+1,-1,OMEGA)*
     $CLEBI(2*RL,1,RJ,OMEGA+1,-1,OMEGA)
      DO 5 MYP=1,NU
      DO 5 MY=1,NU
      NYP=LEVEL(MYP)
      NY=LEVEL(MY)
      Q20S(MYP,MY)=Q20S(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)*
     $2.0*RME*SQRT((RJ+1.)/(RJP+1.))*CGC

      Q22S(MYP,MY)=Q22S(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)*
     $RME*SQRT((RJ+1.)/(RJP+1.))*(CGC1+CGC2)

      IF(RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA)
     $RHO2(MYP,MY)=RHO2(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)*RME

      RJZ(MYP,MY)=RJZ(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)*
     $(HJZ+GAC*HLZ-HAC/SQRT(1.5)*(CGC1-CGC2)*SQRT((RJ+1.)/(RJP+1.))*RME)

    1 CONTINUE
    5 CONTINUE
    6 CONTINUE
    2 CONTINUE
      IF (ISTRCH .EQ. 0) GO TO 21
      DO 7 MYP=1,NU
      DO 7 MY=1,NU
      Q20(MYP,MY)=((1.+EPS2/3.-EPS22*EPS22/6.)*Q20S(MYP,MY)+
     $EPS22/SQRT(6.)*(1.-2.*EPS2/3.)*Q22S(MYP,MY)+
     $(2.*EPS2/3.*(1.+EPS2/3.)-EPS22*EPS22/6.)*RHO2(MYP,MY))/FAC
      Q22(MYP,MY)=((1.+EPS2/3.)*Q22S(MYP,MY)+
     $EPS22/SQRT(24.)*Q20S(MYP,MY)-
     $EPS22/SQRT(6.)*RHO2(MYP,MY))/FAC*(1.-2.*EPS2/3.)
      R2(MYP,MY)=(((1.+EPS2/3.)*(1.-EPS2/3.)-1./3.*(EPS22/2.)**2)*
     $RHO2(MYP,MY)+((1.+EPS2/3.)*EPS2/3.-1./3.*(EPS22/2.)**2)*
     $Q20S(MYP,MY)-EPS22/2.*(1.-2.*EPS2/3.)*SQRT(2./3.)*Q22S(MYP,MY))/
     $FAC
    7 CONTINUE
      GO TO 22
   21 CONTINUE
      DO 23 MYP=1,NU
      DO 23 MY=1,NU
      Q20(MYP,MY) = Q20S(MYP,MY)
      Q22(MYP,MY) = Q22S(MYP,MY)
       R2(MYP,MY) = RHO2(MYP,MY)
   23 CONTINUE
   22 CONTINUE
      WRITE(17) NU,(LEVEL(J),J=1,NU),IDIM,(VALU(J),J=1,IDIM),
     $((Q20(MYP,MY),MY=1,NU),MYP=1,NU),((Q22(MYP,MY),MY=1,NU),MYP=1,NU)
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE Q20Q22
C
C      SAME AS QTME BUT ONLY DIAG MATR ELEM FOR ALL LOW-ENERGY ORBITALS
C***** QUADRUPOLE MOMENTS IN UNITS OF 'HBAR/MW0'
C        Q20=SQRT(16*PI/5)*R2*Y20   (=2Z2-X2-Y2)
C        Q22=SQRT(4*PI/5)*R2*(Y22+Y2-2)   (=SQRT(3/2)*(X2-Y2))
C
      COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22
      COMMON /DIM/ IDIM,KATT
      CHARACTER KATT*4
      COMMON /KVANT/ KV1(161),KV2(161),KV3(161),KV4(161)
      COMMON/EIG/ EIGVEC(161,161),VALU(161)
      COMMON /QMOM/ Q20S(70),Q22S(70),Q20(70),Q22(70),RHO2(70)
     $,R2(70),RJZ(70)
      COMMON /NVEC/ NVEC,NQVEC
      COMMON/FACTOR/ GAC,HAC
      COMMON/BAS/ ISTRCH,ICORR,IREC
      COMMON/OMROT/ OMROT

      INTEGER RNP,RLP,RJP,OMEGAP,RN,RL,RJ,OMEGA
      DATA PI/3.1415926536/
      FAC=(1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)*(1.-2.*EPS2/3.)
      DO 3 NY=1,NQVEC
      Q20S(NY)=0.0
      Q22S(NY)=0.0
      Q20(NY)=0.0
      Q22(NY)=0.0
      RHO2(NY)=0.0
      R2(NY)=0.0
      RJZ(NY)=0.0
    3 CONTINUE
      DO 2 IP=1,IDIM
      RNP=KV1(IP)
      RLP=KV2(IP)
      RJP=KV3(IP)
      OMEGAP=KV4(IP)
      DO 2 I=1,IDIM
      RN=KV1(I)
      RL=KV2(I)
      RJ=KV3(I)
      OMEGA=KV4(I)
C***************** IS THIS COMBINATION OF BASIS VECTORS SIGNIFICANT? ***
C     EIGMAX = 1.E-8
C     DO 25 NY=1,NVEC
C        EIGABS = ABS(EIGVEC(IP,NY)*EIGVEC(I,NY))
C        EIGMAX = MAX(EIGABS,EIGMAX)
C  25 CONTINUE
C     IF (EIGMAX .LT. 1.E-5) GO TO 2
C
C     ABCOUP = ABS(EPS4) + ABS(EPS6) + ABS(OMROT)
C     IF ((ABCOUP) .LT. 0.001 .AND. RNP .NE. RN .AND. ISTRCH .EQ. 1)
C    *GO TO 2
      IF(IABS(RNP-RN).GT.2.OR.IABS(RJP-RJ).GT.4.OR.IABS(RLP-RL).GT.2.
     $OR.IABS(OMEGAP-OMEGA).GT.4) GO TO 6
      CGC=0.0
      CGC1=0.0
      CGC2=0.0
      RME=0.0
      HLZ=0.0
      HJZ=0.0
      RME=RMATEL(2*RNP,2*RLP,2*RN,2*RL,2)
      IF(OMEGAP.EQ.OMEGA) CGC=CLEBI(RJ,4,RJP,-1,0,-1)*
     $CLEBI(RJ,4,RJP,OMEGA,0,OMEGA)
      IF(OMEGAP.EQ.OMEGA+4) CGC1=CLEBI(RJ,4,RJP,-1,0,-1)*
     $CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)
      IF(OMEGAP.EQ.OMEGA-4) CGC2=CLEBI(RJ,4,RJP,-1,0,-1)*
     $CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP)
      IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA)
     $HJZ=FLOAT(OMEGA)/2.
      IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.OMEGAP.EQ.OMEGA) HLZ=
     $(OMEGA-1.)/2.*CLEBI(2*RL,1,RJP,OMEGA-1,1,OMEGA)*
     $CLEBI(2*RL,1,RJ,OMEGA-1,1,OMEGA)+
     $(OMEGA+1.)/2.*CLEBI(2*RL,1,RJP,OMEGA+1,-1,OMEGA)*
     $CLEBI(2*RL,1,RJ,OMEGA+1,-1,OMEGA)
      DO 5 NY=1,NQVEC
      Q20S(NY)=Q20S(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)*
     $2.0*RME*SQRT((RJ+1.)/(RJP+1.))*CGC
      Q22S(NY)=Q22S(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)*
     $RME*SQRT((RJ+1.)/(RJP+1.))*(CGC1+CGC2)
      IF(RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA)
     $RHO2(NY)=RHO2(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)*RME
      RJZ(NY)=RJZ(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)*
     $(HJZ+GAC*HLZ-HAC/SQRT(1.5)*(CGC1-CGC2)*SQRT((RJ+1.)/(RJP+1.))*RME
     &*(RNP-RN)/2.)
    5 CONTINUE
    6 CONTINUE
    2 CONTINUE
      IF (ISTRCH .EQ. 0) GO TO 21
      DO 7 NY=1,NQVEC
      Q20(NY)=((1.+EPS2/3.-EPS22*EPS22/6.)*Q20S(NY)+
     $EPS22/SQRT(6.)*(1.-2.*EPS2/3.)*Q22S(NY)+
     $(2.*EPS2/3.*(1.+EPS2/3.)-EPS22*EPS22/6.)*RHO2(NY))/FAC
      Q22(NY)=((1.+EPS2/3.)*Q22S(NY)+
     $EPS22/SQRT(24.)*Q20S(NY)-
     $EPS22/SQRT(6.)*RHO2(NY))/FAC*(1.-2.*EPS2/3.)
      R2(NY)=(((1.+EPS2/3.)*(1.-EPS2/3.)-1./3.*(EPS22/2.)**2)*
     $RHO2(NY)+((1.+EPS2/3.)*EPS2/3.-1./3.*(EPS22/2.)**2)*Q20S(NY)-
     $EPS22/2.*(1.-2.*EPS2/3.)*SQRT(2./3.)*Q22S(NY))/FAC
    7 CONTINUE
      GO TO 22
   21 CONTINUE
      DO 23 NY=1,NQVEC
      Q20(NY) = Q20S(NY)
      Q22(NY) = Q22S(NY)
       R2(NY) = RHO2(NY)
   23 CONTINUE
   22 CONTINUE
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE DECOUP
C
C***** CALCULATION OF GENERALIZED DECOUPLING FACTORS ETC.
C      IN THE CASE OF STRETCHED COORDINATES (ISTRCH=1) THE DECOUPLING
C      FACTORS ARE APPROXIMATED BY THEIR "STRETCHED VALUES"
C               (CF. APPENDIX, NUCL PHYS 307 (1978) 189
C      IF ICORR=1, <J+>/<J->/<JZ> MATRIX ELEMENTS ARE CORRECTED AS
C               DESCRIBED IN THE APPENDIX, NPA 307,189 (1978)
C               AND CORRECTED D, E, AND JZ2 ARE CONTRUCTED FROM
C               J+, J-, JZ BY SUMMING OVER INTERMEDIATE STATES
C
C      GENERALISED DECOUPLING FACTORS: APLB, C, D, E
C      APLB = A+B = <J+>,  C = <J->, D = <J**2-JZ**2>, E = <J+**2+J-**2>J-**2>
C      NORM: NORMALIZATION INTEGRAL
C      JZ = <JZ>, JZ2 = <JZ**2>
C      SNOLL = <SZ>, SPLUS = <S+>, SMINUS = <S->
C      R2SZ = <(R**2)*SZ>,.... R2SMIN = <(R**2)*S->,... ETC
C      D2 = -(SQRT(10)/2)(S*C(2))**1  -  TENSOR OF RANK 1
C      *** SPHERICAL COMPONENTS CALC FOR D2, S+,S-,J+,J- FOR S AND J **
C      C(2M) = SQRT(4*PI/5)*Y(2M) WHERE Y(2M) SPHERICAL HARMONIC
C      Q20 ETC. TRANSFERRED FROM QTME FOR PRINTING
C
C      IF IREC=1, THE 'RECOIL TERMS' (J+J-, J+J+, JZ2, ETC) ARE TREATED AS
C      TWO BODY TERMS IN ASYR, WHICH REQUIRES THAT THE MATRIX ELEMENTS OF
C      J+, J- AND JZ BE CALCULATED FOR A COMPLETE SET OF INTERMEDIATE STATES.
C
      COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161)
      COMMON/EIG/ EIGVEC(161,161),VALU(161)
      COMMON/DIM/ IDIM,KATT
      CHARACTER KATT*4
      COMMON/QUADR/ Q20S(15,15),Q22S(15,15),Q20(15,15),Q22(15,15),
     *RHO2(15,15),R2(15,15),RJZ(15,15)
      COMMON/MATREL/ APLB(15,15),C(15,15),D(15,15),E(15,15),
     *NORM(15,15),JZ(15,15),JZ2(15,15),
     *SNOLL(15,15),SPLUS(15,15),SMINUS(15,15),
     *R2SZ(15,15),R2SPLU(15,15),R2SMIN(15,15),
     *R4SZ(15,15),R4SPLU(15,15),R4SMIN(15,15),
     *R2JZ(15,15),R2JPLU(15,15),R2JMIN(15,15),
     *R4JZ(15,15),R4JPLU(15,15),R4JMIN(15,15),
     *R2DZ(15,15),R2DPLU(15,15),R2DMIN(15,15),
     *R4DZ(15,15),R4DPLU(15,15),R4DMIN(15,15)
      REAL JZ,JZ2,NORM
C
C**** EQUIVA IS INTRODUCED TO REFER TO ALL MATRIX ELEMENTS
C**** SIMULTANEOUSLY: INITIALIZE TO ZERO, MAKE SYMMETRIC AND OUTPUT
      REAL EQUIVA(15,15,28)
      EQUIVALENCE (APLB,EQUIVA)
C
C***  JZR, JPLR, AND JMNR ARE NEEDED TO CALCULATE MATRIX ELEMENTS OF
C***  THE RECOIL TERM IN ASYR IF IT IS TREATED AS A TWO BODY TERM.
C***  SZR, ...,FZR ARE NEEDED ALSO IF THE CORRECTION DUE TO STRETCHING
C***  IS DESIRED
C
C****      FMNR, FPLR, FZR AS DEFINED IN APP., NUCL PHYS A307, 189
C**** NOTE THAT FPLR IS PROPORTIONAL TO Y2,-1
C****           FMNR                    Y2,+1
C****           FZR                     Y2,+2 - Y2,-2
C
      DIMENSION JPLR(161,15),JMNR(161,15),JZR(161,15),
     &          SPLR(161,15),SMNR(161,15),SZR(161,15),
     &          FPLR(161,15),FMNR(161,15),FZR(161,15)
      REAL JZR,JPLR,JMNR
      REAL JPLRS
C
CQS** Q2TEST = Q20S, SZ=SNOLL, SPLU=-SPLUS/SQRT(2), SMIN=SMINUS/SQRT(2)
CQS**                                          INTRODUCED AS TEST CASES
CQS   REAL Q2TEST(15,15),
CQS  *SZ(15,15),SPLU(15,15),SMIN(15,15)

      INTEGER RNP,RLP,RJP,OMEGAP,RN,RL,RJ,OMEGA
      COMMON/LEVEL/ NU,LEVEL(15)
      COMMON/BAS/ ISTRCH,ICORR,IREC
      COMMON/CMN7/ EPSP,GAMMA
      DATA PI/3.1415926536/
      ICORRS=ICORR*ISTRCH
      IF (ICORRS .EQ. 1) write(6,101)
  101 FORMAT(1H ,20(1H*),' NOTE: ICORR=1 => J+, J-, JZ, J+J-, J+J-,'
     & ' JZ2, CORRECTIONS  DUE TO STRETCHED BASIS APPLIED ',20(1H*),///)
      write(6,100)
  100 FORMAT(1H ,46(1H*),'  SINGLE PARTICLE MATRIX ELEMENTS  ',47(1H*),
     *///'   ENY   EMY     A+B     C      D      E    JZ  ',
     *'  S0    S+    S-  <R2S> <R4S>  R2DZ  R2D+  R2D- ',
CC   *' JZ2    S+    S-  <R2S> <R4S>  R2DZ  R2D+  R2D- ',
     *' Q20S  Q22S   Q20   Q22  RHO2    R2',/)
C 100 FORMAT(1H ,42(1H*),'  SINGLE PARTICLE MATRIX ELEMENTS  ',42(1H*),
C    *///'      ENY    EMY     A+B    C      D      E  ',
C    *'  <JZ>   <JZ2>  <S0>   <S+>  <S-> ',
C    *'<R2SZ> <R2S+> <R2S-> <R4SZ> <R4S+> <R4S->',/,
C    *'  <R2JZ> <R2J+> <R2J-> <R4JZ> <R4J+> <R4J->',
C    *' <R2DZ> <R2D+> <R2D-> <R4DZ> <R4D+> <R4D->',
C    *' Q20S   Q22S    Q20    Q22    RHO2     R2',/)
      IF (NU .LE. 15) GO TO 10
      write(6,80) NU
   80 FORMAT(1H1,'NYIC=',I3)
      STOP
   10 CONTINUE
      DO 1001 I=1,NU
      DO 1001 J=1,NU
      DO 1001 K=1,28
         EQUIVA(I,J,K) = 0.
 1001 CONTINUE
      IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN
         DO 1002 I=1,IDIM
         DO 1002 J=1,NU
            JZR (I,J) = 0.0
            JPLR(I,J) = 0.0
            JMNR(I,J) = 0.0
            SZR (I,J) = 0.0
            SPLR(I,J) = 0.0
            SMNR(I,J) = 0.0
            FZR (I,J) = 0.0
            FPLR(I,J) = 0.0
            FMNR(I,J) = 0.0
 1002    CONTINUE
      ENDIF
CQS   DO 1000 NN=1,NU
CQS   DO 1000 MM=1,NU
CQS      SZ(NN,MM)=0.0
CQS      SMIN(NN,MM)=0.0
CQS      SPLU(NN,MM)=0.0
CQS      Q2TEST(NN,MM)=0.0
C1000 CONTINUE
C
C***** COMPUTE MATRIX ELEMENTS <N,L,J,OMEGA / OP / NP,LP,JP,OMEGAP>
C*****                     AND <N,L,J,OMEGA / OP / NP,LP,JP,-OMEGAP>
C
C***** FIRST LOOP OVER BASIS FUNCTIONS (INITIAL STATE)
C
      DO 11 IP = 1,IDIM
      RNP = KV1(IP)
      RLP = KV2(IP)
      RJP = KV3(IP)
      OMEGAP = KV4(IP)
C
C***** SECOND LOOP OVER BASIS FUNCTIONS (FINAL STATE)
C
      DO 12 I=1,IDIM
      RN = KV1(I)
      RL = KV2(I)
      RJ = KV3(I)
      OMEGA = KV4(I)
C***************** IS THIS COMBINATION OF BASIS VECTORS SIGNIFICANT? ***
C     EIGMAX = 1.E-8
C     DO 21 MMY=1,NU
C     DO 22 NNY=MMY,NU
C        MY = LEVEL(MMY)
C        NY = LEVEL(NNY)
C        EIGABS = ABS(EIGVEC(IP,MY)*EIGVEC(I,NY))
C        EIGMAX = MAX(EIGABS,EIGMAX)
C  22 CONTINUE
C  21 CONTINUE
C     IF (EIGMAX .LT. 1.E-5) GO TO 12
C
CQ    IF (ABS(RL-RLP) .GT. 2) GO TO 12
CC    IF (ABS(RJ-RJP) .GT. 2) GO TO 12
C
C*********************************************RADIAL MATRIX ELEMENT*****
      IF (RN .NE. RNP .OR. RL .NE. RLP) THEN
         RAD0 = 0.
      ELSE
         RAD0 = 1.
      ENDIF
      RAD2 = RMATEL(2*RNP,2*RLP,2*RN,2*RL,2)
      CALL MRHOT(RNP,RLP,0,RN,RL,0,4,RAD4)
C
C**** REDUCED MATRIX ELEMENT FOR TENSOR OPERATOR D *** SEE EDMONDS FOR *
C                                                  *** DEFINITIONS *****
      REDMD = REDMD1(2*RL,1,RJ,2*RLP,1,RJP)
C
C**** REDUCED MATRIX ELEMENT OF S USED AS A TEST CASE ***************
CS
CS    IF (RL .EQ. RLP) REDMS = 0.5*SQRT(3.)*(-1.)**((RJP-OMEGAP)/2)*
CS   *(CLEBI(2*RL,1,RJ,0,1,1)*CLEBI(2*RLP,1,RJP,0,1,1)-
CS   *CLEBI(2*RL,1,RJ,2,-1,1)*CLEBI(2*RLP,1,RJP,2,-1,1))/
CS   *CLEBI(RJ,RJP,2,1,-1,0)
C
      IF (OMEGAP .NE. OMEGA) GO TO 41
C
C**** MATRIX ELEM OF SCALAR OPERATORS AND ZERO COMP OF TENSOR OPERATORS
C
      X = (RJ*(RJ+2.) - OMEGA*OMEGA)/4.
      XA = OMEGA/2.
      XB = FLOAT(OMEGA*OMEGA)/4.
      XD = (-1.)**((RJP-OMEGAP)/2)/SQRT(3.)*
     *CLEBI(RJ,RJP,2,OMEGA,-OMEGAP,0)
      CG = CLEBI(2*RL,1,RJP,OMEGAP-1, 1,OMEGAP)*
     *     CLEBI(2*RL,1,RJ ,OMEGA -1, 1,OMEGA )-
     *     CLEBI(2*RL,1,RJP,OMEGAP+1,-1,OMEGAP)*
     *     CLEBI(2*RL,1,RJ ,OMEGA +1,-1,OMEGA )
CQ    MJ=OMEGA
CQ    XQ=2.*SQRT((2*RLP+1.)/(2*RL+1.))*
CQ   *(CLEBI(2*RL,1,RJ,MJ-1,1,MJ)*CLEBI(2*RLP,1,RJP,MJ-1,1,MJ)*
CQ   *CLEBI(4,2*RLP,2*RL,0,MJ-1,MJ-1) +
CQ   *CLEBI(2*RL,1,RJ,MJ+1,-1,MJ)*CLEBI(2*RLP,1,RJP,MJ+1,-1,MJ)*
CQ   *CLEBI(4,2*RLP,2*RL,0,MJ+1,MJ+1))*
CQ   *CLEBI(4,2*RLP,2*RL,0,0,0)

C******************************* START LOOPS OF NILSSON ORBITALS *******
C
      DO 31 MMY=1,NU
      DO 32 NNY=MMY,NU
      MY = LEVEL(MMY)
      NY = LEVEL(NNY)
      EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
      IF (RNP .NE. RN .OR. RLP .NE. RL .OR. RJP .NE. RJ) GO TO 51
         NORM(NNY,MMY) = NORM(NNY,MMY) + EIGV12
         D   (NNY,MMY) = D   (NNY,MMY) + EIGV12*X
         JZ2(NNY,MMY) = JZ2(NNY,MMY) + EIGV12*XB
   51 CONTINUE
      IF (RLP .NE. RL) GO TO 52
         SNOLL(NNY,MMY) = SNOLL(NNY,MMY) + EIGV12*CG*0.5*RAD0
         R2SZ (NNY,MMY) = R2SZ (NNY,MMY) + EIGV12*CG*0.5*RAD2
         R4SZ (NNY,MMY) = R4SZ (NNY,MMY) + EIGV12*CG*0.5*RAD4
CS       SZ(NNY,MMY) =   SZ(NNY,MMY) + EIGV12*XD*REDMS*RAD0
   52 CONTINUE
      IF (RJP .NE. RJ) GO TO 53
         JZ  (NNY,MMY) = JZ  (NNY,MMY) + EIGV12*XA*RAD0
         R2JZ(NNY,MMY) = R2JZ(NNY,MMY) + EIGV12*XA*RAD2
         R4JZ(NNY,MMY) = R4JZ(NNY,MMY) + EIGV12*XA*RAD4
   53 CONTINUE
         R2DZ(NNY,MMY) = R2DZ(NNY,MMY) + EIGV12*XD*REDMD*RAD2
         R4DZ(NNY,MMY) = R4DZ(NNY,MMY) + EIGV12*XD*REDMD*RAD4
CQ       Q2TEST(NNY,MMY) = Q2TEST(NNY,MMY) + EIGV12*XQ*RAD2
   32 CONTINUE
   31 CONTINUE
C
      IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN
         DO 131 MMY=1,NU
            MY = LEVEL(MMY)
            IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 132
caes            DO 132 NY=1,IDIM
            DO NY=1,IDIM ! AES
               EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
               IF (RLP .EQ. RL)
     &            SZR(NY,MMY) = SZR(NY,MMY) + EIGV12*CG*0.5*RAD0
               IF (RJP .EQ. RJ)
     &            JZR(NY,MMY) = JZR(NY,MMY) + EIGV12*XA*RAD0
	      ENDDO   ! AES
  132       CONTINUE
  131    CONTINUE
      ENDIF
C
C********************************** END LOOPS OVER NILSSON ORBITALS ****
   41 CONTINUE
      IF (OMEGAP .NE. -OMEGA+2) GO TO 42
C
C**** (+1) COMPONENT OF TENSOR OPERATORS ******************************
C
      X = (-1.)**((RJ-OMEGA)/2)*0.5*SQRT((RJ+OMEGA)*(RJ-OMEGA+2.))
      XD = (-1.)/SQRT(3.)*CLEBI(RJ,RJP,2,OMEGA,OMEGAP,2)
      CGX=CLEBI(2*RL,1,RJP,OMEGAP-1,1,OMEGAP)*
     *CLEBI(2*RL,1,RJ,-(OMEGA-1),-1,-OMEGA)*(-1.)**((RJ-OMEGA)/2)
C
      IF (IABS(RNP-RN) .EQ. 2 .AND. IABS(RLP-RL) .LE. 2) THEN
         CGF = (-1)**((RJP-OMEGAP)/2)*FLOAT(RNP-RN)/SQRT(1.5)*
     &SQRT((RJP+1.0)/(RJ+1.0))*CLEBI(RJP,4,RJ,-OMEGAP,2,OMEGA)*
     &CLEBI(RJP,4,RJ,-1,0,-1)
      ELSE
         CGF = 0.0
      ENDIF
C****************************** START LOOPS OVER NILSSON ORBITALS ******
C
      DO 33 MMY=1,NU
      DO 34 NNY=MMY,NU
      MY = LEVEL(MMY)
      NY = LEVEL(NNY)
      EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
      IF (RLP .NE. RL) GO TO 62
         SPLUS(NNY,MMY) = SPLUS(NNY,MMY) + EIGV12*CGX*RAD0
         R2SPLU(NNY,MMY) = R2SPLU(NNY,MMY) + EIGV12*CGX*RAD2
         R4SPLU(NNY,MMY) = R4SPLU(NNY,MMY) + EIGV12*CGX*RAD4
CS       SPLU(NNY,MMY) = SPLU(NNY,MMY) + EIGV12*XD*REDMS*RAD0
   62 CONTINUE
      IF (RJP .NE. RJ) GO TO 63
          APLB(NNY,MMY) = APLB(NNY,MMY) + EIGV12*X*RAD0
          R2JPLU(NNY,MMY) = R2JPLU(NNY,MMY) + EIGV12*X*RAD2
          R4JPLU(NNY,MMY) = R4JPLU(NNY,MMY) + EIGV12*X*RAD4
   63 CONTINUE
          R2DPLU(NNY,MMY) = R2DPLU(NNY,MMY) + EIGV12*XD*REDMD*RAD2
          R4DPLU(NNY,MMY) = R4DPLU(NNY,MMY) + EIGV12*XD*REDMD*RAD4
   34 CONTINUE
   33 CONTINUE
C
      IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN
         DO 133 MMY=1,NU
            MY = LEVEL(MMY)
            IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 134
            DO NY=1,IDIM
               EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
               IF (RLP .EQ. RL)
     &            SPLR(NY,MMY) = SPLR(NY,MMY) + EIGV12*CGX*RAD0
               IF (RJP .EQ. RJ)
     &            JPLR(NY,MMY) = JPLR(NY,MMY) + EIGV12*X*RAD0
               FMNR(NY,MMY) = FMNR(NY,MMY) + EIGV12*CGF*RAD2
	      ENDDO
  134       CONTINUE
  133    CONTINUE
      ENDIF
C
C***************************** END LOOPS OVER NILSSON ORBITALS *********
   42 CONTINUE
      IF (OMEGAP .NE. -OMEGA-2) GO TO 43
C
C***** (-1) COMPONENT OF VECTOR OPERATORS ***********************
C
      X = (-1.)**((RJ-OMEGA)/2)*0.5*SQRT((RJ-OMEGA)*(RJ+OMEGA+2.))
      XD = (-1.)/SQRT(3.)*CLEBI(RJ,RJP,2,OMEGA,OMEGAP,-2)
      CGX=CLEBI(2*RL,1,RJP,OMEGAP+1,-1,OMEGAP)*
     *CLEBI(2*RL,1,RJ,-(OMEGA+1),1,-OMEGA)*(-1.)**((RJ-OMEGA)/2)
      IF (IABS(RNP-RN) .EQ. 2 .AND. IABS(RLP-RL) .LE. 2) THEN
         CGF = (-1)**((RJP-OMEGAP)/2)*FLOAT(RNP-RN)/SQRT(1.5)*
     &SQRT((RJP+1.0)/(RJ+1.0))*CLEBI(RJP,4,RJ,-OMEGAP,-2,OMEGA)*
     &CLEBI(RJP,4,RJ,-1,0,-1)
      ELSE
         CGF = 0.0
      ENDIF
C******************************* START LOOPS OVER NILSSON ORBITALS *****
C
      DO 35 MMY=1,NU
      DO 36 NNY=MMY,NU
      MY = LEVEL(MMY)
      NY = LEVEL(NNY)
      EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
      IF (RLP .NE. RL) GO TO 72
         SMINUS(NNY,MMY) = SMINUS(NNY,MMY) + EIGV12*CGX*RAD0
         R2SMIN(NNY,MMY) = R2SMIN(NNY,MMY) + EIGV12*CGX*RAD2
         R4SMIN(NNY,MMY) = R4SMIN(NNY,MMY) + EIGV12*CGX*RAD4
CS       SMIN(NNY,MMY) = SMIN(NNY,MMY) + EIGV12*XD*REDMS*RAD0
   72 CONTINUE
      IF (RJP .NE. RJ) GO TO 73
         C   (NNY,MMY) = C   (NNY,MMY) + EIGV12*X*RAD0
         R2JMIN(NNY,MMY) = R2JMIN(NNY,MMY) + EIGV12*X*RAD2
         R4JMIN(NNY,MMY) = R4JMIN(NNY,MMY) + EIGV12*X*RAD4
   73 CONTINUE
         R2DMIN(NNY,MMY) = R2DMIN(NNY,MMY) + EIGV12*XD*REDMD*RAD2
         R4DMIN(NNY,MMY) = R4DMIN(NNY,MMY) + EIGV12*XD*REDMD*RAD4
   36 CONTINUE
   35 CONTINUE
C
      IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN
         DO 135 MMY=1,NU
            MY = LEVEL(MMY)
            IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 136
            DO NY=1,IDIM
               EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
               IF (RLP .EQ. RL)
     &            SMNR(NY,MMY) = SMNR(NY,MMY) + EIGV12*CGX*RAD0
               IF (RJP .EQ. RJ)
     &            JMNR(NY,MMY) = JMNR(NY,MMY) + EIGV12*X*RAD0
               FPLR(NY,MMY) = FPLR (NY,MMY) + EIGV12*CGF*RAD2
	      ENDDO
  136       CONTINUE
  135    CONTINUE
      ENDIF
C
C***************************** END LOOPS OVER NILSSON ORBITALS *********
   43 CONTINUE
      IF (OMEGAP .NE.  OMEGA-4 .OR. ICORRS .EQ. 1) GO TO 44
C
C***** -2 AND +2 COMPONENTS OF TENSOR OPERATORS ******************
C
      IF (RJ .NE. RJP .OR. RN .NE. RNP .OR. RL .NE. RLP) GO TO 44
      IF ((RJ+OMEGA-2) .LT. 0) GO TO 44
      X = SQRT((RJ+OMEGA)*(RJ-OMEGA+2.)*(RJ+OMEGA-2)*(RJ-OMEGA+4.))/4.
C******************************** START LOOPS OVER NILSSON ORBITALS ****
C
      DO 38 MMY=1,NU
      DO 37 NNY=MMY,NU
      MY = LEVEL(MMY)
      NY = LEVEL(NNY)
      EIGV21 = EIGVEC(IP,NY)*EIGVEC(I,MY)
      EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
      E(NNY,MMY) = E(NNY,MMY) + (EIGV12+EIGV21)*X
   37 CONTINUE
   38 CONTINUE
C
C********************************** END LOOPS OVER NILSSON ORBITALS ****
   44 CONTINUE
C
C***** CORRECTION TO JZ DUE TO STRETCHING OF BASIS, NONAXIALITY
C
      IF (ICORRS .NE. 1 .OR. IABS(RNP-RN) .NE. 2) GO TO 45
      MU = OMEGA - OMEGAP
      IF (IABS(MU) .NE. 4 .OR. IABS(RLP-RL) .GT. 2) GO TO 45
      CGF = FLOAT(RNP-RN)/SQRT(6.)*SQRT((RJP+1.0)/(RJ+1.0))*
     *      CLEBI(RJP,4,RJ,-1,0,-1)*CLEBI(RJP,4,RJ,OMEGAP,MU,OMEGA)
      IF (MU .EQ. 4) CGF = -CGF
C******************************** START LOOPS OVER NILSSON ORBITALS ****
C
      DO 139 MMY=1,NU
         MY = LEVEL(MMY)
         IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 140
         DO NY=1,IDIM
            EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY)
            FZR(NY,MMY) = FZR(NY,MMY) + EIGV12*CGF*RAD2
	   ENDDO
  140    CONTINUE
  139 CONTINUE
C
C********************************** END LOOPS OVER NILSSON ORBITALS ****
   45 CONTINUE
C
   12 CONTINUE
   11 CONTINUE
C
C***************************** END LOOPS OVER BASIS FUNCTIONS *********

C**** CORRECT J+/J-/JZ DUE TO STRETCHED BASIS IF ICORR=1
C
      IF (ICORRS .EQ. 1) THEN
         GAMMAR = GAMMA*PI/180.0
         OMX = 1.-2./3.*EPSP*COS(GAMMAR+2.*PI/3.)
         OMY = 1.-2./3.*EPSP*COS(GAMMAR-2.*PI/3.)
         OMZ = 1.-2./3.*EPSP*COS(GAMMAR)
         F1 = SQRT(OMZ/OMY)
         F2 = SQRT(OMZ/OMX)
C
         FAC1 = (F1 + 1./F1 + 1./F2 + F2)/4.
         FAC2 = (F1 + 1./F1 - 1./F2 - F2)/4.
         FAC3 = (F1 - 1./F1 + 1./F2 - F2)/4.
         FAC4 = (F1 - 1./F1 - 1./F2 + F2)/4.
         FAC5 = (F2/F1 + F1/F2)/2.
         FAC6 = (F2/F1 - F1/F2)/2.
C
         DO 125 MMY=1,NU
         DO 125 NY=1,IDIM
            JPLRS = JPLR(NY,MMY)
            JPLR(NY,MMY) = FAC1*JPLRS        + (1.-FAC1)*SPLR(NY,MMY)
     &                    +FAC2*JMNR(NY,MMY) -     FAC2 *SMNR(NY,MMY)
     &                    +FAC3*FPLR(NY,MMY) +     FAC4 *FMNR(NY,MMY)
            JMNR(NY,MMY) = FAC1*JMNR(NY,MMY) + (1.-FAC1)*SMNR(NY,MMY)
     &                    +FAC2*JPLRS        -     FAC2 *SPLR(NY,MMY)
     &                    +FAC3*FMNR(NY,MMY) +     FAC4 *FPLR(NY,MMY)
            JZR (NY,MMY) = FAC5*JZR(NY,MMY)  + (1.-FAC5)* SZR(NY,MMY)
     &                    +FAC6*FZR(NY,MMY)
  125    CONTINUE
         DO 126 MMY=1,NU
         DO 126 NNY=1,NU
            D  (NNY,MMY) = 0.0
            E  (NNY,MMY) = 0.0
            JZ2(NNY,MMY) = 0.0
  126    CONTINUE
         DO 129 MMY=1,NU
            MY = LEVEL(MMY)
         DO 128 NNY=MMY,NU
            NY = LEVEL(NNY)
            APLB(NNY,MMY) = JPLR(NY,MMY)
            C   (NNY,MMY) = JMNR(NY,MMY)
            JZ  (NNY,MMY) = JZR (NY,MMY)
         DO 127 IK=1,IDIM
            D  (NNY,MMY) = D(NNY,MMY) + (JPLR(IK,MMY)*JPLR(IK,NNY)
     &                    + JMNR(IK,MMY)*JMNR(IK,NNY))/2.
            E  (NNY,MMY) = E(NNY,MMY) + JPLR(IK,MMY)*JMNR(IK,NNY)
     &                    + JMNR(IK,MMY)*JPLR(IK,NNY)
            JZ2(NNY,MMY) = JZ2(NNY,MMY)  + JZR(IK,MMY)*JZR(IK,NNY)
  127    CONTINUE
  128    CONTINUE
  129    CONTINUE
      ENDIF
C************** START LOOPS OVER NILSSON ORBITALS FOR PRINTING ********
C
      DO  1 MMY=1,NU
      MY = LEVEL(MMY)
      DO  1 NNY=MMY,NU
      AVR2S=
     $(ABS(R2SZ(NNY,MMY))+ABS(R2SPLU(NNY,MMY))+ABS(R2SMIN(NNY,MMY)))/
     $(ABS(SNOLL(NNY,MMY))+ABS(SPLUS(NNY,MMY))+ABS(SMINUS(NNY,MMY)))
      AVR4S=
     $(ABS(R4SZ(NNY,MMY))+ABS(R4SPLU(NNY,MMY))+ABS(R4SMIN(NNY,MMY)))/
     $(ABS(SNOLL(NNY,MMY))+ABS(SPLUS(NNY,MMY))+ABS(SMINUS(NNY,MMY)))
      NY = LEVEL(NNY)
      IF(MY.EQ.NY) write(6,200) VALU(NY),KATT,APLB(NNY,MMY),
     $C(NNY,MMY),D(NNY,MMY),E(NNY,MMY),JZ(NNY,MMY),
CC   $JZ2(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S,
     $SNOLL(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S,
     $R2DZ(NNY,MMY),R2DPLU(NNY,MMY),R2DMIN(NNY,MMY),
     $Q20S(NNY,MMY),Q22S(NNY,MMY),Q20(NNY,MMY),Q22(NNY,MMY),
     $RHO2(NNY,MMY),R2(NNY,MMY)
      IF(NY.GT.MY) write(6,300) VALU(NY),KATT,VALU(MY),KATT,
     $APLB(NNY,MMY),C(NNY,MMY),D(NNY,MMY),E(NNY,MMY),JZ(NNY,MMY),
CC   $JZ2(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S,
     $SNOLL(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S,
     $R2DZ(NNY,MMY),R2DPLU(NNY,MMY),R2DMIN(NNY,MMY),
     $Q20S(NNY,MMY),Q22S(NNY,MMY),Q20(NNY,MMY),Q22(NNY,MMY),
     $RHO2(NNY,MMY),R2(NNY,MMY)
CC    PRINT 423,APLB(NNY,MMY),                                          ,
CC   $C(NNY,MMY),D(NNY,MMY),E(NNY,MMY),JZ(NNY,MMY),JZ2(NNY,MMY)
  423 FORMAT(' ',10F9.4)
  200 FORMAT(/,1X,F5.2,A1,' MY=NY',2F7.3,2F7.2,15F6.2)
  300 FORMAT(1X,2(F5.2,A1),2F7.3,2F7.2,15F6.2)
      DO 3 K=1,28
      EQUIVA(MMY,NNY,K) = EQUIVA(NNY,MMY,K)
    3 CONTINUE
    1 CONTINUE
C
C************* END LOOPS OVER NILSSON ORBITALS FOR PRINTING ***********
C
      WRITE(17) KATT,IREC,NU,(LEVEL(J),J=1,NU),
     $(((EQUIVA(MYP,MY,I),MY=1,NU),MYP=1,NU),I=1,28)
C     IF (ISTRCH .EQ. 1) GO TO 9
      WRITE(17) IDIM,(KV1(I),KV3(I),KV4(I),I=1,IDIM)
      DO 2 MMY=1,NU
      MY=LEVEL(MMY)
      WRITE(17) (EIGVEC(I,MY),I=1,IDIM)
    2 CONTINUE
    9 CONTINUE
C
      IF (IREC .EQ. 1) WRITE(17) ((JPLR(NY,MMY),MMY=1,NU),NY=1,IDIM),
     &                           ((JMNR(NY,MMY),MMY=1,NU),NY=1,IDIM),
     &                           ((JZR(NY,MMY),MMY=1,NU),NY=1,IDIM)
      RETURN
      END
C
C#######################################################################
C
      FUNCTION REDMD1(L,S,J,LP,SP,JP)
C
C**** REDUCED MATRIX ELEMENT FOR TENSOR OPERATOR D1 OF RANK ONE FORMED *
C**** FROM THE COUPLING OF C2 AND S; D1 = -SQRT(10)/2*(C2*S)1          *
C****   DEFINITIONS  FROM EDMOMDS: ANGULAR MOMENTUM IN QUANT MECHANICS *
C
C**** L S,J,LP,SP,JP HAVE TWICE THEIR PHYSICAL VALUES, I.E. S=SP=1

      INTEGER L,S,J,LP,SP,JP
      REAL NINEJ
      IF   ((LP+L) .GE. 4 .AND. ABS(LP-L) .LE. 4
     * .AND.(JP+J) .GE. 2 .AND. ABS(JP-J) .LE. 2) GO TO 1
      REDMD1 = 0.
      RETURN
    1 CONTINUE
C
C**** SIX-J SYMBOLS FROM CLOSED FORMULAS IN EDMONDS, (6.3.3) AND (6.3.4)
C
      IF (JP .EQ. (LP-1)) THEN
         XT = (J+LP+5.)*(LP+3.-J)/4.
         XN = LP*(LP+1.)*12.
         SIXJA = (-1.)**((J+LP+3)/2)*SQRT(XT/XN)
      ELSEIF (JP .EQ. (LP+1)) THEN
         XT = (J+3.-LP)*(J+LP-1.)/4.
         XN = (LP+1.)*(LP+2.)*12.
         SIXJA = (-1.)**((J+LP+3)/2)*SQRT(XT/XN)
      ELSE
         STOP 6
      ENDIF
C
      IF (J .EQ. (L-1)) THEN
         XT = (LP+L+6.)*(L+4.-LP)/4.
         XN = L*(L+1.)*20.
         SIXJB = (-1.)**((L+LP)/2)*SQRT(XT/XN)
      ELSEIF (J .EQ. (L+1)) THEN
         XT = (LP+4.-L)*(LP+L-2.)/4.
         XN = (L+1.)*(L+2.)*20.
         SIXJB = (-1.)**((L+LP)/2)*SQRT(XT/XN)
      ELSE
         STOP 66
      ENDIF
C
      SIXJC = SQRT(3.)/6
C
C******************** NINE-J SYMBOL CALCULATED FROM EDMONDS (6.4.17)

      NINEJ = -SIXJA*SIXJB/SIXJC/3.
C
C*** REDUCED MATRIX ELEMENTS OF C2 =SQRT(4*PI/5)*Y2 AND S
C
      REDMC2 = SQRT(L+1.)*CLEBI(L,4,LP,0,0,0)
      REDMS1 = SQRT(1.5)
C
C************************ REDUCED MATRIX ELEMENT OF PRODUCT OPERATOR D1,
C                        ********************* SEE EDMOMDS (7.1.5) *****
      C = -SQRT(10.)/2.
      REDMD1 = C*REDMC2*REDMS1*SQRT((JP+1.)*(J+1.)*3.)*NINEJ
C
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE ENERLJ(OMOM,SNI,NOYES)
C
C***** CALCULATION OF SINGLE PARTICLE ENERGIES, WAVE FUNCTIONS
C      AND MATRIX ELEMENTS
C   ******  MAX 10 N-SHELLS  *******************************************
C
CV
      COMMON/LEVEL/NU,LEVEL(15)
CV
      COMMON /LU/ LU
      COMMON /COUPL/ NUU
      COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22
      COMMON/CMN7/ EPSP,GAMMA
      COMMON /QSTR/ Q
      COMMON/CMN1/ KAPPA,DUMMY,MY,L2F(18)
      COMMON/E/ EMIN,EMAX
      COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161)
C     COMMON/EIG/ EIGVEC(203,70),VALU(203)
      COMMON/EIG/ EIGVEC(161,161),VALU(161)
      COMMON/DIM/ IDIM,KATT
      COMMON /NUCL/ Z,A
      COMMON /QMOM/ Q20S(70),Q22S(70),Q20(70),Q22(70),RHO2(70)
     $,R2(70),RJZ(70)
      COMMON/EQFAC/ EFAC,QFAC,ISO
      COMMON /PAR/ IPAR
      COMMON /NVEC/ NVEC,NQVEC
      COMMON/OMROT/ OMROT
      COMMON/FACTOR/ FAC,HAC
      COMMON/XKAMY/ XKAPPA(13,7),XMY(13,7),AVKAMY(13)
      COMMON/ITRANS/ ITRANS,NSHELL
      COMMON/BAS/ ISTRCH,ICORR,IREC
!      DIMENSION AM(13041)
      DIMENSION AM(0:25425)
      DIMENSION E(572),EIGVAL(70)
      DIMENSION KVA(161),LEVNUM(280)
      DIMENSION EB(280)
      INTEGER SNI,SN,Q,RN,RL,RJ,RNP,RLP,RJP,OMEGA,OMEGAP
      DIMENSION DIM(2),TIME(2),TID(2)
      INTEGER DIM,DIME
      INTEGER Z,A,Z2
      INTEGER QQ
      REAL KAPPA,MY
      COMMON/Q2022/ QQ20(280),QQ22(280)
     $,RR2(280),RRJZ(280)
      CHARACTER IPLUS*4,IMINUS*4,KATT*4,IPAR*4,KVANT(2)*4,KVB(280)*4
      CHARACTER KV(728)*4
      DATA MATDIM/161/
      DATA IPLUS,IMINUS/'+','-'/
C     DATA DELTA/0.0316/
      DATA DELTA/0.000316/
      DATA PI/3.1415926536/
      SN=SNI
      IF (SNI .LE. 10) GO TO 200
  199 FORMAT(1X,' N IS GREATER THAN MAX VALUE 10  ...  ENERLJ ')
      write(6,199)
      STOP
  200 CONTINUE
  100 FORMAT(/4X,'OMROT=',F5.3)
      GAMMAR=GAMMA*PI/180.0
      OMX=1.-2./3.*EPSP*COS(GAMMAR+2.*PI/3.)
      OMY=1.-2./3.*EPSP*COS(GAMMAR-2.*PI/3.)
      OMZ=1.-2./3.*EPSP*COS(GAMMAR)
      FAC=0.5*(SQRT(OMY/OMX)+SQRT(OMX/OMY))-1.
      HAC=0.5*(SQRT(OMY/OMX)-SQRT(OMX/OMY))
      OMROT=OMROT/OMOM
      Q=0
      QQ=0
      IR=1
      IF(ABS(OMROT).GT.0.0001) IR=2
C
C        LOOP IN SIGNATURES (OMROT .GT. 0)
C
      DO 201 ISYM=1,IR
         IS=ISYM
         IB=0
C
C        LOOP IN PARITIES
C
      DO 201 NUM=1,2
 7132    CONTINUE
         MATNUM=IABS(NUM-3*IB)
         SN=SNI+1-MATNUM
         IF((SN/2)*2.EQ.SN) GO TO 7133
         IDIM=(SN+1)*(2*SN+7)*(SN+3)/24
         KATT=IMINUS
         GO TO 7135
 7133    CONTINUE
         IDIM=(SN+2)*(2*SN+3)*(SN+4)/24
         KATT=IPLUS
 7134    CONTINUE
 7135    CONTINUE
         IF(NUM.EQ.1.AND.KATT.EQ.IPAR.AND.MATNUM.EQ.1) THEN
            IB=1
            GO TO 7132
         ENDIF
         IF(NOYES.EQ.0) GO TO 101
         IF (NUM .EQ. 2 .OR. IS .EQ. 2) write(6,397)
  397    FORMAT(//)
  101    CONTINUE
         DIM(MATNUM)=IDIM
         KVANT(MATNUM)=KATT
         KLAFS=IDIM*(IDIM+1)/2
         DO 2 KK=1,KLAFS
            AM(KK)=0.0
    2    CONTINUE
         DO 1 I=1,IDIM
            CALL QNUMLJ(I,RN,RL,RJ,OMEGA,SN)
            KV1(I)=RN
            KV2(I)=RL
            KV3(I)=RJ
            KV4(I)=OMEGA*(3-2*IS)
    1    CONTINUE
!         IRK=MSLEFT(XX)
         KK=0
C
C           LOOP OVER BASIS FUNCTIONS, CREATION OF HAM MATR ELEMENTS
C
         DO 210 I=1,IDIM
            RNP=KV1(I)
            RLP=KV2(I)
            RJP=KV3(I)
            OMEGAP=KV4(I)
C
C           SECOND LOOP OVER BASIS FUNCTIONS (J .GE. I)
C
         DO 210 J=I,IDIM
            RN=KV1(J)
            RL=KV2(J)
            RJ=KV3(J)
            OMEGA=KV4(J)
            KK=KK+1
            IF(RNP.GT.RN+NUU.OR.OMEGAP.GT.OMEGA+8.OR.
     $      OMEGAP.LT.OMEGA-8) GO TO 303
            E20=0.0
            E22=0.0
            E40=0.0
            E42=0.0
            E44=0.0
            E60=0.0
            ELZ=0.0
            END=0.0
            IF(OMEGAP.EQ.OMEGA+8.OR.OMEGAP.EQ.OMEGA-8) GO TO 305
            IF(OMEGAP.EQ.OMEGA+4.OR.OMEGAP.EQ.OMEGA-4) GO TO 306
            IF(RNP.EQ.RN) E20=-2.0*EPS2/3.0*CLEBI(RJ,4,RJP,-1,0,-1)*
     $      CLEBI(RJ,4,RJP,OMEGA,0,OMEGAP)
            IF (ISTRCH .EQ. 1) GO TO 901
            IF(RNP.EQ.RN+2) E20=-2.0*EPS2/3.0*CLEBI(RJ,4,RJP,-1,0,-1)*
     $      CLEBI(RJ,4,RJP,OMEGA,0,OMEGAP)
  901       CONTINUE
C           E40=EPS4*(COS(1.5*GAMMAR)**2+0.375*SIN(1.5*GAMMAR)**2)*
C**** EPS4 FOR GAMMA .NE. 0 OR 60 AS DEFINED E.G. IN NUCL PHYS
C     A436 (1985) 14   ###### E40, E42 AND E44 ####
C                             SEE ALSO FUNCTION FVCON  ###########
C
C**** ALSO NOTE TYPO IN THAT PAPER : V4 TERM OMITS THE FACTOR
C     2*SQRT(9/4*PI), REQUIRED TO GIVE 2*EPS4*P4 FOR GAMMA=0
C     AS IS CONVENTIONAL, E.G. NUCL. PHYS. A129 (1969) 445
C
            E40=EPS4*((5.*COS(GAMMAR)**2+1.)/6.)*
     $      CLEBI(RJ,8,RJP,-1,0,-1)*CLEBI(RJ,8,RJP,OMEGA,0,OMEGAP)
            E60=EPS6*CLEBI(RJ,12,RJP,-1,0,-1)*
     $      CLEBI(RJ,12,RJP,OMEGA,0,OMEGAP)
            IF(RNP.EQ.RN.AND.RLP.EQ.RL) ELZ=-OMROT*FAC*
     $      ((OMEGA-1.)/2.*CLEBI(2*RL,1,RJP,OMEGA-1,1,OMEGA)*
     $      CLEBI(2*RL,1,RJ,OMEGA-1,1,OMEGA)+
     $      (OMEGA+1.)/2.*CLEBI(2*RL,1,RJP,OMEGA+1,-1,OMEGA)*
     $      CLEBI(2*RL,1,RJ,OMEGA+1,-1,OMEGA))
            GO TO 304
  306       CONTINUE
            IF(RNP.EQ.RN) E22=EPS22/SQRT(6.0)*CLEBI(RJ,4,RJP,-1,0,-1)*
     $      (CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)+
     $      CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP))
            IF (ISTRCH .EQ. 1) GO TO 902
            IF(RNP.EQ.RN+2) E22=EPS22/SQRT(6.0)*CLEBI(RJ,4,RJP,-1,0,-1)*
     $      (CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)+
     $      CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP))
  902       CONTINUE
            IF(RNP.EQ.RN+2)
     $      END= OMROT/SQRT(1.5)*HAC*CLEBI(RJ,4,RJP,-1,0,-1)*
     $      (CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)-
     $      CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP))
C           E42=-EPS4*SQRT(10.)/8.*SIN(1.5*GAMMAR)**2*
            E42=-EPS4*SQRT(30.)/12.*SIN(2.*GAMMAR)*
     $      CLEBI(RJ,8,RJP,-1,0,-1)*
     $      (CLEBI(RJ,8,RJP,OMEGA,4,OMEGAP)+
     $      CLEBI(RJ,8,RJP,OMEGA,-4,OMEGAP))
            GO TO 304
  305       CONTINUE
C           E44=EPS4*SQRT(70.)/16.*SIN(1.5*GAMMAR)**2*
            E44=EPS4*SQRT(70.)/12.*SIN(GAMMAR)**2*
     $      CLEBI(RJ,8,RJP,-1,0,-1)*
     $      (CLEBI(RJ,8,RJP,OMEGA,8,OMEGAP)+
     $      CLEBI(RJ,8,RJP,OMEGA,-8,OMEGAP))
  304       CONTINUE
            CALL MRHOT(RNP,RLP,0,RN,RL,0,2,RME)
            ARJ=FLOAT(RJ)
            ARJP = FLOAT(RJP)
            ARL = FLOAT(RL)
            AM(KK) = RME*SQRT((ARJ+1.)/(ARJP+1.))*
     *      (E20+E22+E40+E42+E44+E60+END)+ELZ
            IF(I.EQ.J) AM(KK)=AM(KK)+FLOAT(RN)+1.5-XKAPPA(RN+1,RL/2+1)/
     *      OMOM*(ARJ*(ARJ+2.)/4.-ARL*(ARL+1.)-0.75)-
     *      XKAPPA(RN+1,RL/2+1)/OMOM*XMY(RN+1,RL/2+1)*ARL*(ARL+1.)+
     *      AVKAMY(RN+1)/OMOM*FLOAT(L2F(RN+1))-
     *      OMROT*FLOAT(OMEGA)/2.
C
C           AM IS SYMMETRICAL MATRIX TO BE DIAGONALIZED
C                  (ONLY STORED FOR (J .GE. I)
C
            AM(KK)=-AM(KK)
  303       CONTINUE
  210       CONTINUE
C
C
c         JRK=MSLEFT(XX)
         TIME(MATNUM)=(IRK-JRK)*0.001

         IDI=IDIM
         NQVEC=MIN0(IDI,70)
         IF (IREC .EQ. 1) THEN
            NVEC=IDIM
         ELSE
            NVEC=NQVEC
         ENDIF
         IF(NOYES.EQ.0) GO TO 11
C*****    PRINT AV BASFUNKTIONER
         DIME=DIM(MATNUM)
         IMIN=(DIME+4)/5
         KVR1=3-2*ISYM
         write(6,9) KVR1
    9    FORMAT(4X,'SYMMETRY R1=',I2,/)
         write(6,7) KVANT(MATNUM)
    7    FORMAT(/,'   BASIS FUNCTIONS /N,L,J,OMEGA>,   PARITY ',A1,//)
         DO 113 K=1,IMIN
            IMAX = K + 4*IMIN
            IF (IMAX .GT. DIME) IMAX = IMAX - IMIN
            write(6,16)(J,KV1(J),KV2(J),KV3(J),KV4(J),J=K,IMAX,IMIN)
   16       FORMAT(1X,5(I3,1X,1H/,2I2,I3,2H/2,I3,2H/2,1H>,3X))
  113    CONTINUE
   11    CONTINUE
         IF(NOYES.EQ.0) NVEC=0
!         IRK=MSLEFT(XX)
C
C        DIAGONALIZATION
C
         CALL BIGMAX(AM,VALU,IDIM,IDI,NVEC,EIGVEC,MATDIM)
C
c         JRK=MSLEFT(XX)
         TID(MATNUM)=(IRK-JRK)*0.001


         DO 405 J=1,IDI
            VALU(J)=-VALU(J)
  405    CONTINUE
         IF(NOYES.EQ.0) GO TO 12
C*****    PRINT AV VAVE-FUNCTIONS
         write(6,398) KATT,EMIN,EMAX
  398    FORMAT(//,'    SINGLE PARTICLE (PARITY',A1,') WAVE FUNCTIONS ',
     *   'IN THE ENERGY INTERVAL',F4.1,3H<E<,F3.1)
         K=0
         DO 404 J=1,NVEC
            IF(VALU(J).LT.EMIN.OR.VALU(J).GT.EMAX) GO TO 403
            K=K+1
            EIGVAL(K)=VALU(J)
            KK=0
            DO 402 I=1,IDIM
               IF(ABS(EIGVEC(I,J)).LT.DELTA) GO TO 401
               KK=KK+1
               AM(KK)=EIGVEC(I,J)
               KVA(KK)=I
  401          CONTINUE
  402       CONTINUE
            KKMAX=KK
            write(6,600) EIGVAL(K),(AM(KK),KVA(KK),KK=1,KKMAX)
C 600       FORMAT(/1X,2HE=,F6.4,2X,5(9(F6.3,1H/,I3,1H>)/11X))
  600       FORMAT(/1X,2HE=,F6.4,2X,5(7(F9.6,1H/,I3,1H>)/11X))
  403       CONTINUE
  404    CONTINUE
CV
CV    PRINT OF SINGLE-PARTICLE WAVEFUNCTIONS FOR VMI
CV
      IF (KATT .NE. IPAR) GOTO 801
      WRITE(16) IDIM
      DO 800 II=1,NU
      IIII=LEVEL(II)
      IF (IPAR .EQ. IPLUS) WRITE(16) (KV1(III)/2+1,KV2(III),
     $KV3(III),KV4(III),EIGVEC(III,IIII),III=1,IDIM)
      IF (IPAR .EQ. IMINUS) WRITE(16) ((KV1(III)+1)/2,KV2(III),
     $KV3(III),KV4(III),EIGVEC(III,IIII),III=1,IDIM)
  800 CONTINUE
  801 CONTINUE
CV
CV
C*** TRANSFORM TO /N,L,LAMBDA,SIGMA> BASIS (NECESSARY FOR Vpn CALCN)
CV
         ITEST=NUU*(ABS(EPS4)+ABS(EPS6))
         IF (ITEST .EQ. 0) THEN
            NMAX=NSHELL
            IDIM1=((NSHELL+1)*(NSHELL+2))/2
         ELSE
            IDIM1=IDIM
            NMAX=SN
         ENDIF
         IF (KATT .EQ. IPAR .AND. ITRANS .NE. 0)
     &   CALL TRANSF(NMAX,IDIM1)
C
C*****   CALCULATION OF QUADRUPOLE MATRIX ELEMENTS
C           SUBROUTINES Q20Q22 AND QTME
C
!         IRKA=MSLEFT(XX)
         CALL Q20Q22
!         IRK = MSLEFT(XX)
         IF(KATT .EQ. IPAR .AND. OMROT .LT. 0.0001) THEN
            CALL QTME
         ENDIF
c         JRKA=MSLEFT(XX)
         SEC = (IRKA-IRK)*0.001
         SEC1 = (IRK-JRKA)*0.001
         write(6,55)SEC,SEC1
   55    FORMAT(/,1X,8H Q20Q22:,F6.2,13H SEC    QTME:,F6.2,4H SEC)
CT       WRITE(1,*)'NQVEC,QQ',NQVEC,QQ
         DO 204 I=1,NQVEC
            IQ=I+QQ
            EB(IQ)=VALU(I)
            KVB(IQ)=KATT
            QQ20(IQ)=Q20(I)
            QQ22(IQ)=Q22(I)
            RR2(IQ)=R2(I)
            RRJZ(IQ)=RJZ(I)
  204    CONTINUE
         QQ=QQ+NQVEC
         IF(NUM.NE.2.) GO TO 406
         IF(IS.NE.IR) GO TO 406
CT       WRITE(1,*)'NQVEC,QQ,NUM,IS,IR,Z',NQVEC,QQ,NUM,IS,IR,Z
CT       WRITE(1,*)'EB,KVB',(EB(I),I=1,QQ),(KVB(I),I=1,QQ)
C
C        SORTING ACCORDING TO SINGLE-PARTICLE ENERGIES
C
         CALL SORT4(EB,KVB,QQ20,QQ22,RR2,RRJZ,QQ)
         LEVPOS=0
         LEVNEG=0
         DO 309 I=1,QQ
         IF (KVB(I).EQ.IPLUS) THEN
            LEVPOS=LEVPOS+1
            LEVNUM(I)=LEVPOS
         ELSE
            LEVNEG=LEVNEG+1
            LEVNUM(I)=LEVNEG
         ENDIF
  309    CONTINUE
C
CT       WRITE(1,*)'EB,KVB',(EB(I),I=1,QQ),(KVB(I),I=1,QQ)
CT       WRITE(1,*)'A',A
         N=A-Z
C********QUADRUPOLE MOMENTS OF NUCLEUS
         Q0=0.0
         Q2=0.0
         IF (ISO .EQ. 1) THEN
            Z2=Z/2
         ELSEIF (ISO .EQ. 2) THEN
            Z2=N/2
         ENDIF
         IZ=Z2*IS
         DO 205 I=1,IZ
            Q0=Q0+QQ20(I)*2./FLOAT(IS)
            Q2=Q2+QQ22(I)*2./FLOAT(IS)
  205    CONTINUE
         GAM=180./PI*ATAN(-SQRT(2.)*Q2/Q0)
         IF (ISO .EQ. 1) THEN
            write(6,35) Q0,Q2,GAM,Z,A
         ELSEIF (ISO .EQ. 2) THEN
            write(6,35) Q0,Q2,GAM,N,A
         ENDIF
   35    FORMAT(/1X,3HQ0=,F8.3,' HBAR/MW0',5X,3HQ2=,F8.3,' HBAR/MW0',
     *   5X,6HGAMMA=,F5.1,' DEG   FOR',I3,' PARTICLES;  A=',I3,/)
  406    CONTINUE
         IF(KATT.NE.IPAR) GO TO 202
C
C*****   CALC OF GENERALIZED DECOUPLING FACTORS, SUBROUTINE DECOUP
!         IRKA = MSLEFT(XX)
         IF(KATT.EQ.IPAR.AND.OMROT.LT.0.0001) THEN
             write(6,399)
             CALL DECOUP
         ENDIF
!         JRKA=MSLEFT(XX)
         SEC=(IRKA-JRKA)*0.001
         write(6,65) SEC
   65    FORMAT(/,1X,'DECOUP:',F6.2,' SEC')
  399    FORMAT(1H1)
   19    CONTINUE
  202    CONTINUE
   12    CONTINUE
         DO 203 I=1,IDI
            IQ=I+Q
            E(IQ)=VALU(I)
            KV(IQ)=KATT
  203    CONTINUE
         Q=Q+IDI
  201 CONTINUE
      IF(NOYES.EQ.0) GO TO 17
      write(6,8)
    8 FORMAT(///,1X,2('  #   ENERGY +/-(#)    <Q20>    <Q22>     <R2>',
     *'     <JZ>',4X),/)
      IMIN = (QQ+1)/2
      IF (IMIN .GT. 40) IMIN = 40
      DO 111 K=1,IMIN
         IMAX = K + IMIN
         IF (IMAX .GT. QQ) IMAX = IMIN
         write(6,14) (J,EB(J),KVB(J),LEVNUM(J),QQ20(J),QQ22(J),RR2(J),
     *   RRJZ(J),J = K,IMAX,IMIN)
C        IF (MOD(K,5) .EQ. 0) PRINT 24
  24     FORMAT(/)
  14     FORMAT(1X,2(I3,F9.4,2X,A1,'(#',I2,')',4F9.4,3X))
  111 CONTINUE
      WRITE(17) EPSP,GAMMA,EPS4,EPS6,OMOM,QQ,
     $(EB(I),I=1,QQ),(KVB(I),I=1,QQ),
     $(QQ20(I),I=1,QQ),(QQ22(I),I=1,QQ),(RR2(I),I=1,QQ),(RRJZ(I),I=1,QQ)
   17 CONTINUE
C
C***** PRINT OF S P ENERGIES AND QUANTUM NUMBERS  -  NOYES = 0
      IF(NOYES.NE.0) GO TO 114
      CALL SORT(E,KV,Q)
      IF(NOYES.NE.0) write(6,399)
      write(6,6)
    6 FORMAT(1X,'    LEVEL,ENERGY,PARITY',/)
      IMIN=MIN0(44,Q)
      DO 112 K=1,IMIN
         IMAX=K+176
         IMAX=MIN0(IMAX,Q)/44*44+K
         IF(IMAX.GT.Q) IMAX=IMAX-44
         write(6,15) (J,E(J),KV(J),J=K,IMAX,44)
   15    FORMAT(1X,5(I7,F9.4,4X,A1,3X))
  112 CONTINUE
  114 CONTINUE
C***** PRINT AV TIDER.
      write(6,40)
   40 FORMAT(//)
      DO 10 I=1,2
         write(6,30) KVANT(I),DIM(I),TIME(I),TID(I)
   30    FORMAT(4X,9H PARITY: ,A1,15H     DIMENSION:,I4,12H     MATRIX:
     $   ,F6.1,4H SEC,21H     DIAGONALISATION:,F6.1,4H SEC)
   10 CONTINUE
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE QNUMLL(NUM,NP,LP,LAMP,ISIGP,NMAX)
C
C     SETS UP /N,L,LAMBDA,SIGMA> BASIS FUNCTIONS
C     (ONE FOR EACH CALL) BUT WITH THE RESTRICTION
C     OMEGA = LAMBDA + SIGMA =..., -3/2, +1/2, +5/2, ...
C     ON OUTPUT: NP=N, LP=L, LAMP=LAMBDA, ISIGP=2*SIGMA
C        INPUT:  NUM=NUMBERING ON BASIS FUNCTION
C                NMAX=MAXIMUM N-SHELL (PARITY=(-1)**NMAX)
C
      NUMP=NUM
      NTRY=NMAX
    1 CONTINUE
      NTEST=(NTRY+1)*(NTRY+2)/2
      IF(NUMP .LE. NTEST) GO TO 2
      NUMP=NUMP-NTEST
      NTRY=NTRY-2
      GO TO 1
    2 CONTINUE
      NP=NTRY
C
      LTRY=NP
    3 CONTINUE
      LTEST=2*LTRY +1
      IF(NUMP .LE. LTEST) GO TO 4
      NUMP=NUMP-LTEST
      LTRY=LTRY-2
      GO TO 3
    4 CONTINUE
      LP=LTRY
      LAMP=LP+1-NUMP
      ISIGP=-1
      IF((LAMP/2)*2 .EQ. LAMP) ISIGP=1
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE TRANSF(NMAX,IDIM1)
C
C     TRANSFORMS S.P. STATES FROM THE /N,L,J,OMEGA> BASIS
C     TO THE /N,L,LAMBDA,SIGMA> BASIS (CALLED FROM ENERLJ)
C
      COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161)
      COMMON/EIG/ EIGVEC(161,161),VALU(161)
      COMMON/DIM/ IDIM,KATT
      COMMON/LEVEL/ NU,LEVEL(15)
      COMMON/ITRANS/ ITRANS,NSHELL
      DIMENSION KVL1(161),KVL2(161),KVL3(161),KVL4(161),KDOM(161)
      DIMENSION T1(161,15),AV(15),AM(161)
      INTEGER OMEGA,OMEGAP
      CHARACTER KATT*4
      DATA DELTA/0.00316/
      DO 1 I=1,IDIM1
         CALL QNUMLL(I,N,L,LAM,ISIG,NMAX)
         KVL1(I)=N
         KVL2(I)=L
         KVL3(I)=LAM
         KVL4(I)=ISIG
    1 CONTINUE
      DO 3 I=1,IDIM1
      DO 3 JJ=1,NU
         T1(I,JJ)=0.0
    3 CONTINUE
      DO 150 IJ=1,NU
         JJ=LEVEL(IJ)
         AV(IJ)=VALU(JJ)
         DO 140 I=1,IDIM
            IF (ABS(EIGVEC(I,JJ)) .LT. 0.00001) GO TO 139
            N = KV1(I)
            L = KV2(I)
            J = KV3(I)
            OMEGA = KV4(I)
            DO 130 K=1,IDIM1
               NP = KVL1(K)
               LP = KVL2(K)
               LAM = KVL3(K)
               ISIG = KVL4(K)
               OMEGAP = 2*LAM + ISIG
               IF (NP.NE.N .OR. LP.NE.L .OR. OMEGAP.NE.OMEGA)
     &            GO TO 129
               CG1 = CLEBI(2*L,1,J,2*LAM,ISIG,OMEGA)
               T1(K,IJ) = T1(K,IJ) + CG1*EIGVEC(I,JJ)
  129       CONTINUE
  130       CONTINUE
  139    CONTINUE
  140    CONTINUE
  149 CONTINUE
  150 CONTINUE
         IF (ITRANS .EQ. 1) GO TO 406
C*****    PRINT AV BASFUNKTIONER
         write(6,398)
  398 FORMAT(//5X,'SELECTED SINGLE-PARTICLE STATES EXPANDED IN ',
     &       '/N,L,LAMBDA,SIGMA> BASIS ',//)
         IMIN=(IDIM1+4)/5
         write(6,7)
    7    FORMAT(/,'   BASIS FUNCTIONS /N,L,LAMBDA,SIGMA> ',//)
         DO 113 K=1,IMIN
            IMAX = K + 4*IMIN
            IF (IMAX .GT. IDIM1) IMAX = IMAX - IMIN
            write(6,16)(J,KVL1(J),KVL2(J),KVL3(J),KVL4(J),J=K,IMAX,
     &                IMIN)
   16       FORMAT(1X,5(I3,1X,1H/,4I3,2H/2,1H>,3X))
  113    CONTINUE
C*****    PRINT OUT WAVE-FUNCTIONS
         K=0
         DO 404 J=1,NU
            KK=0
            DO 402 K=1,IDIM1
               IF(ABS(T1(K,J)) .LT. DELTA) GO TO 401
               KK=KK+1
               AM(KK)=T1(K,J)
               KDOM(KK)=K
  401          CONTINUE
  402       CONTINUE
            KKMAX=KK
            write(6,600) AV(J), (AM(KK),KDOM(KK),KK=1,KKMAX)
  600       FORMAT(/2X,'E=',F6.4,1X,5(7(F9.6,1H/,I3,1H>)/11X))
  404    CONTINUE
  406    CONTINUE
C***  WRITE TO FILE2
         WRITE(2,*) IDIM1
         DO 410 J=1,IDIM1
            WRITE(2,*) KVL1(J),KVL2(J),KVL3(J),KVL4(J)
  410    CONTINUE
         WRITE(2,*) NU
         WRITE(2,*) (LEVEL(I),I=1,NU)
         DO 420 K=1,NU
            WRITE(2,*) (T1(J,K),J=1,IDIM1)
  420    CONTINUE
      RETURN
      END
C
C#######################################################################
C
      FUNCTION RMATEL(N,L,NP,LP,LAMBDA)
C
C**** N, L, NP, LP ARE INTEGERS WITH TWICE THEIR PHYSICAL VALUE.
C     RADIAL MATRIX ELEMENTS OF R**LAMBDA (R DIMENSIONLESS)
C     (LAMBDA = 2,3; CLOSED FORMULAS USED)
C     N/2 ,NP/2 - NUMBER OF OSCILLATOR QUANTA
C     L/2, LP/2 - ORBITAL ANGULAR MOMENTUM
C
      RMATEL=0.
      IF (LAMBDA.EQ.2) GO TO 200
      IF (LAMBDA.EQ.3) GO TO 300
      write(6,100)LAMBDA
  100 FORMAT(/,'RADIAL MATRIX ELEMENTS OF R',2H**,I2,
     *' CAN BE SUPPLIED BY MRHOT')
      STOP
  200 IF (NP-N) 210,220,230
  210 IF (NP.NE.N-4) GO TO 9
      IF (LP-L) 211,212,213
  211 IF (LP.NE.L-4) GO TO 9
      RMATEL=SQRT(FLOAT((N+L+2)*(N+L-2)))/4.
      GO TO 9
  212 RMATEL=SQRT(FLOAT((N-L)*(N+L+2)))/4.
      GO TO 9
  213 IF (LP.NE.L+4) GO TO 9
      RMATEL=SQRT(FLOAT((N-L)*(N-L-4)))/4.
      GO TO 9
  220 IF (LP-L) 221,222,223
  221 IF (LP.NE.L-4) GO TO 9
      RMATEL=SQRT(FLOAT((N-L+4)*(N+L+2)))/2.
      GO TO 9
  222 RMATEL=FLOAT(N+3)/2.
      GO TO 9
  223 IF (LP.NE.L +4) GO TO 9
      RMATEL=SQRT(FLOAT((N-L)*(N+L+6)))/2.
      GO TO 9
  230 IF (NP.NE.N+4) GO TO 9
      IF (LP-L) 231,232,233
  231 IF (LP.NE.L-4) GO TO 9
      RMATEL=SQRT(FLOAT((N-L+4)*(N-L+8)))/4.
      GO TO 9
  232 RMATEL=SQRT(FLOAT((N-L+4)*(N+L+6)))/4.
      GO TO 9
  233 IF (LP.NE.L+4) GO TO 9
      RMATEL=SQRT(FLOAT((N+L+6)*(N+L+10)))/4.
      GO TO 9
  300 ND=NP-N
      LD=LP-L
      NA=IABS(ND)
      LA=IABS(LD)
      IF (NA.EQ.2) GO TO 301
      IF (NA.EQ.6) GO TO 303
      GO TO 9
  301 IF (LA.EQ.2) GO TO 310
      IF (LA.EQ.6) GO TO 320
      GO TO 9
  310 IF (LD*ND) 311,9,312
  311 RMATEL=FLOAT(3*N+L+10+ND)*SQRT(FLOAT(N-L+2+ND))/8.
      GO TO 9
  312 RMATEL= FLOAT(3*N-L+8+ND)*SQRT(FLOAT(N+L+4+ND))/8.
      GO TO 9
  320 IF (LD*ND) 321,9,322
  321 RMATEL= SQRT(FLOAT((N+L+4-ND)*(N-L+2*ND)*(N-L+4+2*ND)))*3./8.
      GO TO 9
  322 RMATEL=SQRT(FLOAT((N+L+2+2*ND)*(N+L+6+2*ND)*(N-L+2-ND)))*3./8.
      GO TO 9
  303 IF (LA.EQ.2) GO TO 330
      IF (LA.EQ.6) GO TO 340
      GO TO 9
  330 ND=ND/3
      IF (LD*ND) 331,9,332
  331 RMATEL= SQRT(FLOAT((N+L+4+ND)*(N-L+2*ND)*(N-L+4+2*ND)))/8.
      GO TO 9
  332 RMATEL=SQRT(FLOAT((N+L+2+2*ND)*(N+L+6+2*ND)*(N-L+2+ND)))/8.
      GO TO 9
  340 ND=ND/3
      IF (LD*ND) 341,9,342
  341 RMATEL=SQRT(FLOAT((N-L+4*ND)*(N-L+4+4*ND)*(N-L+2+ND)))/8.
      GO TO 9
  342 RMATEL= SQRT(FLOAT((N+L+2+4*ND)*(N+L+6+4*ND)*(N+L+4+ND)))/8.
    9 RETURN
      END
C
C#######################################################################
C
      SUBROUTINE SORT(E,KV,NI)
C
C     SORTING OF E(1:NI), KV(1:NI) ACCORDING TO INCREASING VALUES IN E
C
      DIMENSION KV(1),E(1)
      INTEGER D,I,J
      REAL Y
	character KV*4,Y1*4
      N=NI
      D=2**INT(1.4427*ALOG(FLOAT(N-1)))-1
  341 IF (D .LE. 0) RETURN
      I=1
  342 J=I
      IPD=I+D
      Y=E(IPD)
      Y1=KV(IPD)
  343 IF (Y.LT.  E(J)) GO TO 344
  345 JPD=J+D
      E(JPD)=Y
      KV(JPD)=Y1
      I=I+1
      IF (I+D .LE. N) GO TO 342
      D=(D-1)/2
      GO TO 341
  344 JPD=J+D
      E(JPD)=E(J)
      KV(JPD)=KV(J)
      J=J-D
      IF (J .GE. 1) GO TO 343
      GO TO 345
      END
C
C#######################################################################
C
      SUBROUTINE SORT4(E,KV1,KV2,KV3,KV4,KV5,NI)
C
C     SORTING OF E, KV1,...., KV5 ACCORDING TO INCREASING VALUES IN E
C       NI: DIMENSION OF E, KV1,..., KV5
C
      DIMENSION E(1),KV1(1),KV2(1),KV3(1),KV4(1),KV5(1)
      REAL KV2,KV3,KV4,KV5
      CHARACTER Y1*4,KV1*4
CT      WRITE(1,*)'UTSKRIFT I SORT4,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)
      Y4=KV4(K)
      Y5=KV5(K)
      Y3=KV3(K)
    3 IF(Y.LT.E(J)) GO TO 5
    4 L=J+M
      E(L)=Y
      KV1(L)=Y1
      KV2(L)=Y2
      KV4(L)=Y4
      KV5(L)=Y5
      KV3(L)=Y3
      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)
      KV4(L)=KV4(J)
      KV5(L)=KV5(J)
      KV3(L)=KV3(J)
      J=J-M
      IF(J.GE.1) GO TO 3
      GO TO 4
      END
C
C#######################################################################
C
      FUNCTION BVCON(FVCON,XA,XB,YA,YB,NX,NY,X,WX,Y,WY,IX,IY)
C
C     DOUBLE INTEGRATION SUBROUTINE FOR VOLUME CONSERVATION CONDITION.
C     FVCON   IS THE EXTERNAL FUNCTION TO BE INTEGRATED.
C     X,WX    ARE THE ARRAYS OF NORMALIZED POINTS AND WEIGHTS FOR X-VARI
C     Y,WY    ARE THE ARRAYS OF NORMALIZED POINTS AND WEIGHTS FOR Y-VARI
C
      DIMENSION X(1),WX(1),Y(1),WY(1)
C     CALCULATION OF THE DENORMALIZING X AND Y SLOPES
      HX=(XB-XA)/(2.*IX)
      HY=(YB-YA)/(2.*IY)
      BVCON=0.0
      DO 10 L=1,IY
C     CALCULATION OF THE DENORMALIZING Y INTERCEPT
      BY=YA+(2.*L-1.)*HY
      DO 10 K=1,IX
C     CALCULATION OF THE DENORMALIZING X INTERCEPT
      BX=XA+(2.*K-1.)*HX
      DO 10 J=1,NY
C     CALCULATION OF THE DENORMALIZED POINT Y
      A=HY*Y(J)+BY
      DO 10 I=1,NX
C     CALCULATION OF THE DENORMALIZED POINT X AND SUMMATION OF THE
C     APPROXIMATING SUM
   10 BVCON=BVCON+WY(J)*WX(I)*FVCON(HX*X(I)+BX,A)
C     ADJUSTMENT OF THE CALCULATED SUM.
      BVCON=HX*HY*BVCON
      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
      COMMON /PTR/ FCT(57)
      INTEGER Z,ZMIN,ZMAX
      J1=I1
      J2=I2
      J =I3
      N=57
      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) GO TO 5
      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) GO TO 5
      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(61,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
C
C#######################################################################
C
      FUNCTION FVCON(XX,YY)
C
C     FUNCTION TO BE INTEGRATED IN VOLUME CALC (STRETCHED COORD)
C
      COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22
      COMMON/CMN7/ EPSP,GAMMA
      X=XX
      Y=YY
      PI=3.1415926536
      GAMMAR=GAMMA*PI/180.0
C     ANM=1.0-EPS2/3.0*(3.*X*X-1.)+0.5*EPS22*(1.-X*X)*COS(2.*Y)+
C    $EPS4/4.0*(COS(1.5*GAMMAR)**2+0.375*SIN(1.5*GAMMAR)**2)*
C    $(35.*X*X*X*X-30.*X*X+3.)-5.*EPS4/8.*SIN(1.5*GAMMAR)**2*(7.*X*X-1.)
C    $*(1.-X*X)*COS(2.*Y)+35.*EPS4/32.*SIN(1.5*GAMMAR)**2*(1.-X*X)*
C    $(1.-X*X)*COS(4.*Y)
      ANM=1.0-EPS2/3.0*(3.*X*X-1.)+0.5*EPS22*(1.-X*X)*COS(2.*Y)+
     $EPS4/4.0*((5.*COS(GAMMAR)**2+1.)/6.)*
     $(35.*X*X*X*X-30.*X*X+3.)
     $-5.*EPS4/12.*SQRT(3.)*SIN(2.*GAMMAR)*(7.*X*X-1.)
     $*(1.-X*X)*COS(2.*Y)+35.*EPS4/24.*SIN(GAMMAR)**2*(1.-X*X)*
     $(1.-X*X)*COS(4.*Y)
      FVCON=1./(4.*PI)/SQRT((1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)*
     $(1.-2.*EPS2/3.))/(ANM**(1.5))
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE HELP
      COMMON /LU/ LU
      COMMON /PTR/ F(57)
      COMMON /GM/ GM(150)
      COMMON /CMN1/ KAPPA,OMOM,MY,L2F(18)
      LU = 6
C*****COMPUTATION OF F(N+1)=LN(N^)
      F(1) = 0.
      DO 101 I=1,56
      AII = FLOAT(I)
      F(I+1) = F(I) + ALOG(AII)
  101 CONTINUE
C*****COMPUTATION OF (L2)AV = L2F(N+1)
      L2F(1) = 0
      DO 106 I = 1,17
      L2F(I+1) = L2F(I) + I + 1
  106 CONTINUE
C*****COMPUTATION OF GAMMA-FUNCTION IN INTEGER AND HALF-INTEGER POINTS.
      GM(52)=1.0
      GM(51)=1.77245385
      DO 34 I=1,28
C     KOLLA OM 28  R EN TILLR CKLIGT H\G GR NS^
      GM(2*I+52)=I*GM(2*I+50)
   34 GM(2*I+51)=(I-0.5)*GM(2*I+49)
      GM(49)=-2.0*GM(51)
      DO 35 I=1,20
      IR=-2*I+49
      IIR=IR+2
   35 GM(IR)=GM(IIR)/(-I-0.5)
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE LGAUSS(X,W,N)
C
C     SUBROUTINE LGAUSS COMPUTES THE ZEROES OF THE LEGENDRE POLYNOMIAL
C     AND THEIR ASSOCIATED WEIGHTS FOR A GAUSSIAN QUADRATURE
C
      DIMENSION X(1),W(1)
      IF(N-1) 1,2,3
C     REQUEST FOR A ZERO POUNT FORMULA IS MEANINGLESS SO RETURN
    1 write(6,11)
   11 FORMAT(//,'REQUEST FOR ZERO POINT FORMULA IN LGAUSS')
      STOP
  2   X(1)=0.0
      W(1)=2.0
      RETURN
C     FOR A ONE POINT FORMULA SEND BACK RESULTS WITHOUT COMPUTING
  3   R=N
      G=-1.
C     THE INITIAL GUESS FOR THE SMALLEST ROOT OF P(N) IS TAKEN AS -1.
      DO 147 I=1,N
      TEST=-2.
      IC=N+1-I
C     WHENEVER WE FIND A ROOT OF THE POLYNOMIAL, ITS NEGATIVE IS ALSO A
C     THE INDEX IC TELLS WHERE TO STORE THE OTHER ROOT
      IF(IC.LT.I) GO TO 150
  4   S=G
      T=1.
      U=1.
      V=0.
C     EVALUATION OF THE N-TH LEGENDRE POLYNOMIAL AND ITS FIRST DERIVATIV
C     WHERE   U=DS/DX
C             V=DT/DX
C             DP=DP/DX
      DO 50 K=2,N
      A=K
      P=((2.0*A-1.0)*S*G-(A-1.0)*T)/A
      DP=((2.0*A-1.0)*(S+G*U)-(A-1.0)*V)/A
      V=U
      U=DP
      T=S
  50  S=P
      IF (ABS(TEST+G) .LT. 0.00000005) GO TO 100
      IF (ABS((TEST-G)/(TEST+G)).LT.0.0000005) GO TO 100
      SUM=0.
      IF(I.EQ.1) GO TO 52
C     THE FOLLOWING COMPUTES THE REDUCED LEGENDRE POLYNOMIAL AND ITS DER
      DO 51 K=2,I
  51  SUM=SUM+1./(G-X(K-1))
  52  TEST=G
      G=G-P/(DP-P*SUM)
      GO TO 4
  100 X(IC)=-G
      X(I)=G
      W(I)=2./(R*T*DP)
      W(IC)=W(I)
  147 G=G-R*T/((R+2.)*G*DP+R*V-2.*R*T*SUM)
  150 RETURN
      END
C
C#######################################################################
C
      SUBROUTINE MRHOT(NNP,LLP,LAMP,NN,LL,LAM,T,ME)
C
C     RADIAL MATRIX ELEMENT, ME, OF R**T (R DIMENSIONLESS)
C     N,NP - NUMBER OF OSCILLATOR QUANTA
C     L,LP - ORBITAL ANGULAR MOMENTUM
C     LAM, LAMP - PROJECTION OF L AND LP (DUMMY IN THIS CASE)
C
      COMMON /LU/ LU
      REAL ME
    1 FORMAT(12I6)
      INTEGER T,TA,E,C,D,SLUT
      COMMON /GM/ GM(150)
  300 FORMAT(7H DANGER)
  301 FORMAT(31H TERMINATION SOMETHING IS WRONG)
  302 FORMAT(33H IS NOT THIS A WONDERFUL COMPILER)
      NANP=NNP
      LALP=LLP
      LAAMP=LAMP
      NAN=NN
      LAL=LL
      LAAM=LAM
      TA=T
      IF(LAAM-LAAMP) 231,201,231
  201 CONTINUE
      NP=(NANP-LALP)/2+1
      N=(NAN-LAL)/2+1
      IF(NP-N)202,203,203
  202 NPNP=NP
      NP=N
      N=NPNP
      L=LALP
      LALP=LAL
      LAL=L
  203 CONTINUE
      C=(LAL-LALP-TA)/2
      D=(LALP-LAL-TA)/2
      E=D+NP-N-1
      IF(2*C.EQ.(LAL-LALP-TA)) GO TO 401
      ME=0.0
      DO 3 K=1,N
      K1=2*(N-K+1)+50
      K2=LAL-LALP-TA+2*(K-1)+50
      K3=LALP-LAL-TA+2*(K+NP-N-1)+50
      K4=2*(NP-N+K)+50
      K5=2*(N-K+1)+LAL+LALP+1+TA+50
      K6=LAL-LALP-TA+50
      K7=LALP-LAL+50-TA
    3 ME=ME+GM(2*N+50)*GM(2*NP+50)/GM(K1)*GM(K2)/GM(2*K+50)*GM(K3)
     */GM(K4)*GM(K5)/GM(K6)/GM(K7)
      K1=2*N+2*LAL+51
      K2=2*NP+2*LALP+51
      ME=ME/SQRT(GM(2*N+50)*GM(K1)*GM(2*NP+50)*GM(K2))*(-1)**(NP-N)
      GO TO 230
  401 CONTINUE
      IF(C.LE.0.AND.D.LE.0)  GO TO 204
      GO TO 208
  204 CONTINUE
      IF(E.GT.-1) GO TO 206
      GO TO 209
  206 SLUT=0
      GO TO  220
  209 CONTINUE
      IF(C.LE.E+1) GO TO 207
      GO TO 221
  207 SLUT=-E
      GO TO 220
  221 SLUT=1-C
      GO TO 210
  208 CONTINUE
      IF(C.LE.0.AND.D.GT.0) GO TO  211
      GO TO 212
  211 SLUT=1-C
      GO TO 213
  212 CONTINUE
      IF(C.GT.0.AND.D.LE.0) GO TO 214
      GO TO 215
  214 CONTINUE
      IF(E.GT.-1) GO TO 216
      GO TO 217
  216 SLUT=0
      GO TO 218
  217 SLUT=-E
  218 CONTINUE
      GO TO 219
  215 CONTINUE
      WRITE(LU,300)
      WRITE(LU,301)
      WRITE(LU,302)
      STOP
  219 CONTINUE
  210 CONTINUE
  220 CONTINUE
  213 CONTINUE
      IF(SLUT.EQ.0) GO TO 231
      IF(SLUT.LT.0) GO TO 222
      GO TO 223
  222 WRITE(LU,300)
      WRITE(LU,301)
      WRITE(LU,302)
      STOP
  223 ME=0.0
      IF(SLUT.GT.N)  SLUT=N
      DO 234 K=1,SLUT
      KK1=2*(N-K+1)+50
      KK2=2*(N-K+1)+LAL+LALP+1+TA+50
      KK3=2*(NP-N+K)+50
      KK4=2-2*C+50
      KK5=2*(2-K-C)+50
      KK6=2*(C+K-1)+50
      KK7=2*C+50
      IF(C.LE.0) GO TO 224
      GO TO 225
  224 GLI=(-1)**(K+1)*GM(KK4)/GM(KK5)
      GO TO 226
  225 GLI=GM(KK6)/GM(KK7)
  226 CONTINUE
      KK8=2-2*D+50
      KK9=2*(2-D-K-NP+N)+50
      KK10=2*(D+K+NP-N-1)+50
      KK11=2*D+50
      IF(D.LE.0.) GO TO 227
      GO TO 228
  227 GLA=GM(KK8)/GM(KK9)*(-1)**(1+K+NP-N)
      GO TO 229
  228 GLA=GM(KK10)/GM(KK11)
  229 CONTINUE
      ME=ME+GM(2*N+50)/GM(2*K+50)*GM(2*NP+50)/GM(KK1)*GM(KK2)/GM(KK3)*
     *GLI*GLA
  234 CONTINUE
      KK12=2*(N+LAL)+51
      KK13=2*(NP+LALP)+51
      ME=ME*(-1)**(N-NP)/SQRT(GM(2*N+50)*GM(2*NP+50)*GM(KK12)*GM(KK13))
      GO TO 230
  231 ME=0.0
  230 CONTINUE
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE OMO(OMOM)
C
C     VOLUME CONSERVATION CALCULATION
C
      COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22
      DIMENSION Z(9),WZ(9)
      COMMON/BAS/ ISTRCH,ICORR,IREC
      EXTERNAL FVCON
C
C     VOLUME CONSERVATION FACTOR CALCULATED WITH DELTA4, DELTA6
C     = 0 IN THE CASE OF SPHERICAL COORDINATES
C
      IF (ISTRCH .EQ. 0) GO TO 6
      IF(ABS(EPS4).LT.0.0001) GO TO 4
      PI=3.1415926536
      CALL LGAUSS(Z,WZ,9)
      NX=1
      NY=1
      A=BVCON(FVCON,-1.0,0.0,0.0,PI/2.,9,9,Z,WZ,Z,WZ,NX,NY)
      OMOM=(8.*A)**(1.0/3.0)
      RETURN
    4 OMOM=((1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)*(1.-2.*EPS2/3.))
     $**(-1.0/3.0)
      RETURN
    6 CONTINUE
      FAC = 2.*EPS2/3.
      OMOM=((1.+FAC+EPS22)*(1.+FAC-EPS22)*(1.-2.*FAC))**(-1./6.)
      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(225),
!     $SQGAMP(225),INTER(1),TMULT(1),T(225,3),SQGAMR(225),UPPER2(225) 
!      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
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      APPROPRIATE FOR NORD COMPUTER
C
CC     FUNCTION MSLEFT(XX)
C      DOUBLE INTEGER IUSE,TUSED
C      IUSE=TUSED(GOLEG)
C      MSLEFT=-IUSE*20
C      RETURN
C      END
C
C***********************************************************************
C     APPROPRIATE FOR VAX (FROM TORD)
C
!      FUNCTION MSLEFT(XX)
! !     LOGICAL STATUS
!      DATA I0,IPID,JPI_CPU/0,0,1031/
!      IF(IPID.EQ.0)THEN
!          STATUS = LIB$GETJPI(I0,IPID,,ISEX,,)
!          IF(.NOT.STATUS)THEN
!              WRITE(*,*)' something wrong with lib$getjpi!'
!          ENDIF
!      ENDIF
!      STATUS = LIB$GETJPI(JPI_CPU,IPID,,ISEX,,)
!      IF(.NOT.STATUS)THEN
!          WRITE(*,*)' something wrong with lib$getjpi!'
!      ENDIF
!      MSLEFT=-ISEX*10
!      RETURN
!      END
C
C***********************************************************************
C     APPROPRIATE FOR VAX (FROM TORD)
C
c      FUNCTION SECOND(XX)
c      LOGICAL STATUS
c      DATA I0,IPID,JPI_CPU/0,0,1031/
c      IF(IPID.EQ.0)THEN
c          STATUS = LIB$GETJPI(I0,IPID,,ISEX,,)
c          IF(.NOT.STATUS)THEN
c              WRITE(*,*)' something wrong with lib$getjpi!'
c          ENDIF
c      ENDIF
c      STATUS = LIB$GETJPI(JPI_CPU,IPID,,ISEX,,)
c      IF(.NOT.STATUS)THEN
c          WRITE(*,*)' something wrong with lib$getjpi!'
c      ENDIF
c      SECOND=FLOAT(ISEX)/100.
c      RETURN
c      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
