C#######################################################################
C########                                                    ###########
C########           PROGRAM PROBAMO                          ###########
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 LATER, EG,
C      <j+>, ETC), AND ASYRMO (DIAGONALIZES THE PARTICLE + TRIAXIAL
C      ROTOR HAMILTONIAN).
C
C########                                                    ###########
C#######################################################################
C
C***** PROBAMO CALCULATES M1/E2 MATRIX ELEMENTS IN THE PARTICLE-ROTOR
C      FORMAILISM FROM THE ENERGY EIGENVECTORS SUPPLIED BY ASYRMO.
C      BOTH DIAGONAL AND OFF-DIAGONAL MATRIX ELEMENTS ARE CALCULATED
C      (IE, STATIC MOMENTS, TRANSITION RATES, MIXING RATIOS, ETC).
C
C      SPECTROSCOPIC FACTORS FOR SINGLE NUCLEON TRANSFER REACTIONS
C      CAN ALSO BE CALCULATED, AS WELL AS THE HYPERFINE ANOMALIES
C      (DUE TO THE DISTRIBUTION OF MAGNETIZATION INSIDE THE NUCLEUS;
C      SEE PHYS REV 123 (1961) 1326)
C
C***********************************************************************
C
C     INPUT OF SINGLE-PARTICLE QUANTITIES FROM GAMPN ON FILE17
C     INPUT OF PARTICLE-ROTOR  QUANTITIES FROM ASYRMO ON FILE18
C
C   INPUT DATA (PARTLY OVERLAPPING WITH ASYRMO):
C
C  CARD 1. FILE17, FILE18  (FREE FORMAT)
C           FILE17: OUTPUT FILE FROM GAMPN,
C                   CONTAINS S.P. MATRIX ELEMENTS
C           FILE18: OUTPUT FILE FROM ASYRMO,
C                   CONTAINS ODD-A WAVE FUNCTIONS
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 (ON FILE17) ARE NOT USED
C
C       3. ISTRCH      (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       4. Z, AA       (FREE FORMAT)
C           Z: THE PROTON NUMBER OF THE ODD-MASS NUCLEUS
C           AA: MASS NUMBER OF THE ODD NUCLEUS
C
C       5. ISPEC, CUTOFF, IQ, GSFAC, GR   (FREE FORMAT)
C           ISPEC: CONTROLS CALCULATION OF SPECTROSCOPIC FACTORS;
C           ISPEC=0 - NO SPEC. FACTORS CALCULATED
C           ISPEC=1 - CALCULATE SPEC. FACTOR IF ISTRCH=0 ALSO
C
C           CUTOFF: ENERGY CUTOFF IN keV FOR TRANS. PROB. AND MOMENTS
C
C           IQ=0 - MACROSCOPIC CALC OF QUADRUPOLE MOMENTS FOR CORE
C                    (FROM A HOMOGENOUS SHARP SURFACE BODY)
C           IQ=1 - THE CORE QUADR MOMENTS ARE CALCULATED FROM THE SINGLE
C                  PARTICLE MOMENTS ASSUMING A SHARP FERMI SURFACE
C           FOR ODD-NEUTRON NUCLEUS,  IQ=0 MUST BE USED
C
C           GSFAC: SCALING FACTOR OF THE INTERNAL SPIN G-FACTOR
C             GS(EFF) = GSFAC*GS(FREE);
C                     (STANDARD: (0.6 .LT. GSFAC .LT. 1.0))
C           GR : G-FACTOR FOR CORE;
C                IF GR<0 ON INPUT THEN STANDARD VALUE, GR= Z/AA, CHOSEN
C
C           NOTE THAT IMIN, ISPIN & KMAX ARE ALL READ FROM FILE18
C           IMIN : MINIMUM SPIN CONSIDERED IS IMIN/2
C           ISPIN: MAXIMUM SPIN CONSIDERED IS ISPIN/2
C           KMAX : MAXIMIUM K CONSIDERED (K TRUNCATION)
C
C       6. BS2(S-STATE), BS4(S-STATE), BS2(P-STATE), BS4(P-STATE)
C           (FREE FORMAT)
C           THE FACTORS OF R**2 AND R**4 FOR THE SPIN-OPERATOR
C           WHEN CALCULATING THE HYPERFINE ANOMALY.
C                 GIVEN FOR ELECTRONIC S- AND P-STATE
C           IF (BS2(I) .EQ. 0) NO CALCULATION OF HYPERFINE ANOMALIES
C
C***********************************************************************
C
      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
C
C**** EQUIVA IS INTRODUCED TO REFER TO ALL MATRIX ELEMENTS
C**** SIMULTANEOUSLY: INPUT, ETC
      REAL EQUIVA(15,15,28)
      EQUIVALENCE (APLB,EQUIVA)
C
      DIMENSION SUM1(15,15),SUM2(15,15),SUM3(15,15)
      DIMENSION SUM11(15,15),SUM12(15,15),SUM13(15,15),SUM14(15,15)
      DIMENSION Q0T(15,15),Q2T(15,15),VALU(286)
      DIMENSION S(286),Q20(286),Q22(286),R2(286),RJZ(286),V2(286)
      DIMENSION K3(15),F(60),IM(60),EIV(15,15,60)
      DIMENSION LEVEL(15),LEVEL1(15)
C     DIMENSION KV1(161), KV2(161), KV3(161), KV4(161), EIGVEC(161,15)
      DIMENSION KV1(161), KV3(161), KV4(161), EIGVEC(161,15)
      DIMENSION SPIKUP(6), SSTRIP(6), UPAIR(15), VPAIR(15)
      REAL C1(20,5,15),C2(20,5,15),C3(20,5,15),
     *     C4(20,5,15),C5(20,5,15),C6(20,5,15)
      REAL PUV(15,15)
      REAL BS2(2), BS4(2), BL2(2), BL4(2)
      REAL SUMI(60),SUMSR0(60),SUMSR2(60),SUMSR4(60),SUMDR2(60),
     $ SUMDR4(60),SUMJR0(60),SUMJR2(60),SUMJR4(60)
      REAL MAGMR,MAGMS1,MAGMS2,MAGML,MAGM1,MAGM2
      REAL NOUSE(161,15)
      INTEGER G,H,Z,AA,Q
      CHARACTER KATT*4,KV(286)*4
      CHARACTER FILE17*24,FILE18*24
      EQUIVALENCE (S,KV,R2,RJZ)
      DATA PI/3.1415928536/


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



C
      READ(5,*)FILE17,FILE18
C
CC    OPEN FOR VAX
C
      OPEN(17,FILE=FILE17,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(18,FILE=FILE18,FORM='UNFORMATTED',STATUS='OLD')
CC
CC** APPROPRIATE FOR NORD COMPUTER
CC
CC    OPEN(17,FILE=FILE17,ACCESS='READ')
CC    OPEN(18,FILE=FILE18,ACCESS='READ')
      CALL HELP
C**************************************** INPUT ********
      READ(5,*)IPKT,ISKIP
      READ(5,*)ISTRCH
      READ(5,*)Z,AA
      READ(5,*)ISPEC,CUTOFF,IQ,GSFAC,GR
      READ(5,*) (BS2(I), BS4(I), I=1,2)
C*************************** SINGLE PARTICLE UNITS OF B(E2) AND B(M1) **
      BE2W = 0.5184/PI*0.36*AA**(4./3.)*0.0001
      BM1MW = 0.01/PI*0.5625
C********************* RADIUS USED FOR CALC OF HYP ANAMOLIES B-FACT ****
      R2B = (1.07*FLOAT(AA)**(1./3.) + 1.5)**2
C*************************** COLLECTIVE FACTORS FOR HYP ANOMALIES ******
      ROTNY1 = 0.7
      ROTNY2 = 0.8*ROTNY1
CC      IF (Z .LE. AA/2) THEN
C
C*** TEST FOR ODD Z OR ODD N
C
      IF ((Z-1)/2 *2 .EQ. Z-1) THEN
         GSFREE = 5.59
         GS = GSFAC*GSFREE
         GL=1.0
         EP=1.0
         EFAC1=(1.-FLOAT(AA-2*Z)/FLOAT(AA)/3.)
         NODD = Z
      ELSE
         GSFREE = -3.83
         GS = GSFAC*GSFREE
         GL=0.
         EP=0.
         EFAC1=(1.+FLOAT(AA-2*Z)/FLOAT(AA)/3.)
         NODD = AA-Z
      ENDIF
         IF (GR .LT. 0.) GR=FLOAT(Z)/FLOAT(AA)
C
C*********************** LOOP OVER DIFFERENT DEFORMATIONS *******
      DO 111 IP=1,IPKT
!      call cpu_time(begin_time)
!      IRK=MSLEFT(XX)
      READ(17) NU,(LEVEL(J),J=1,NU),IDIM,(VALU(J),J=1,IDIM),
     $((Q0T(MYP,MY),MY=1,NU),MYP=1,NU),((Q2T(MYP,MY),MY=1,NU),MYP=1,NU)
CC    READ(17) KATT,NU,(LEVEL(J),J=1,NU),
      READ(17) KATT,JREC,NU,(LEVEL(J),J=1,NU),
     $(((EQUIVA(MYP,MY,I),MY=1,NU),MYP=1,NU),I=1,28)
C    $((APLB(MYP,MY),MY=1,NU),MYP=1,NU),((C(MYP,MY),MY=1,NU),MYP=1,NU),
C    $((D(MYP,MY),MY=1,NU),MYP=1,NU),((E(MYP,MY),MY=1,NU),MYP=1,NU),
C    $((NORM(MYP,MY),MY=1,NU),MYP=1,NU),((JZ(MYP,MY),MY=1,NU),MYP=1,NU
C    $),((JZ2(MYP,MY),MY=1,NU),MYP=1,NU),((SNOLL(MYP,MY),MY=1,NU),MYP=1
C    $,NU),((SPLUS(MYP,MY),MY=1,NU),MYP=1,NU),((SMINUS(MYP,MY),MY=1,NU),
C    $MYP=1,NU)
C*******  CHANGE BY PBS
C     IF (ISTRCH .EQ. 1) GO TO 3
C     READ(17) IDIM, (KV1(K), KV2(K), KV3(K), KV4(K), K=1,IDIM)
      READ(17) IDIM, (KV1(K), KV3(K), KV4(K), K=1,IDIM)
      DO 2 MY = 1,NU
         READ(17) (EIGVEC(K,MY), K=1,IDIM)
    2 CONTINUE
      IF (JREC .EQ. 1) READ(17) (((NOUSE(I,J),J=1,NU),I=1,IDIM),K=1,3)
    3 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 (IP .LE. ISKIP) GO TO 111
      READ(18) IMIN, ISPIN, KMAX, NYIC,
     &(LEVEL1(I),UPAIR(I),VPAIR(I),I=1,NYIC),IEMAX,(V2(K),K=1,IEMAX)
C
C****MANIPULATIONS TO GET THE S.P. MATRIX ELEMENTS IN THEIR PROPER
C    PLACES  (IF LEVEL1(J) .NE. LEVEL(J) OF GAMMA4)
      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),KATT
   41    FORMAT ('0 NO DECOUPLING PARAMETERS CALCULATED FOR LEVEL:',
     $   I3,A1)
         STOP
   33    LEVEL1(J) = IK
   31 CONTINUE
      DO 34 I=1,NYIC
         MYP = LEVEL1(I)
         DO 35 J=1,NYIC
            MY = LEVEL1(J)
            Q0T(I,J) = Q0T(MYP,MY)
            Q2T(I,J) = Q2T(MYP,MY)
            DO 35 K = 1,28
               EQUIVA(I,J,K) = EQUIVA(MYP,MY,K)
   35    CONTINUE
         DO 36 K=1,IDIM
            EIGVEC(K,I) = EIGVEC(K,MYP)
   36    CONTINUE
         LEVEL(I) = LEVEL(MYP)
   34 CONTINUE
      NU = NYIC
C******************** END MANIPULATIONS ****************
      IDIM1 = (KMAX+1)/2
      DO 15 L=1,IDIM1
         K3(L) = KMAX - (L-1)*4
   15 CONTINUE
      IA=0
      IB=0
C********************* LOOP OVER STATES OF NUCLEUS ******
      DO 17 IJK=IMIN,ISPIN,2
         READ(18) LE
         IA=IB+1
         IB=IB+LE
         READ(18) (F(J),IM(J),J=IA,IB),(((EIV(L,NY,J),L=1,15),NY=1,NU),
     $   J=IA,IB), ((PUV(MY,NY), MY=1,NU), NY=1,NU)

C OUTPUT

CC    DO 1777 J=IA,IB
CC    WRITE(6,2777) IA,IB,LE,NU,IMIN,ISPIN
C2777 FORMAT('0  IA,IB,LE,NU,IMIN,IMAX=',6I5,/)
CC       DO 1777 NY=1,NU
CC          WRITE(6,2788) (EIV(L,NY,J),L=1,15)
C2788 FORMAT(1X,15F7.4)
C1777 CONTINUE
C
C OUTPUT
C
   17 CONTINUE
      IF (ISTRCH .EQ. 1) THEN
         write(6,2000) KATT,EPSP,GAMMA,EPS4,EPS6,(LEVEL(J),J=1,NU)
      ELSEIF (ISTRCH .EQ. 0) THEN
         write(6,2004) KATT,EPSP,GAMMA,EPS4,EPS6,(LEVEL(J),J=1,NU)
      ENDIF
 2000 FORMAT(' PARITY=',A1,3X,'EPS=',F6.3,3X,'GAMMA=',F5.1,3X,
     *'EPS4=',F5.3,'   EPS6=',F5.3,3X,'LEVELS:',15I4)
 2004 FORMAT(' PARITY=',A1,3X,'DELTA=',F6.3,3X,'GAMMAD=',F5.1,3X,
     *'DELTA4=',F5.3,'   DELTA6=',F5.3,3X,'LEVELS:',15I4)
      write(6,2001) Z,AA,IQ,GSFAC,GR,IMIN,ISPIN,KMAX
 2001 FORMAT(6H0   Z=,I3,5H   A=,I3,6H   IQ=,I1,9H   GSFAC=,F5.2,
     & 3X,'GR= ',F5.2,3X,'MIN SPIN=',I3,'/2',3X,'MAX SPIN=',I3,'/2',
     & 3X,'KMAX=',I3,'/2')
      NYIC=NU
      EMIN = F(1)
      DO 25 J=2,IB
         IF (EMIN .GT. F(J)) EMIN = F(J)
   25 CONTINUE
      DO 60 J=1,IB
         F(J)=(F(J)-EMIN)*1000.
   60 CONTINUE
      NI=IB
      IF (ISTRCH .EQ. 0 .AND. ISPEC .EQ. 1) THEN
C
C      CALCULATION OF SPECTROSCOPIC FACTORS
C
         I1=MOD(KV1(1),2)
         write(6,400) (K, K=I1,10,2)
         write(6,402)
  400    FORMAT(/,'SPECTROSCOPIC FACTOR = SUM**2/(I+1/2) FOR STRIPPING',
     *   /,24X,'2*SUM**2',5X,'FOR PICKUP',//,
     *   ' SUM(NL J=I):STRIPPING N=',6I15)
  402    FORMAT(12X,'PICKUP',/)
         DO 50 J=1,NI
            II = IM(J)
            DO 51 N=1,6
               SPIKUP(N) = 0.
               SSTRIP(N) = 0.
   51       CONTINUE
            DO 52 K=1,IDIM
               IF (KV3(K) .NE. II) GO TO 52
               N = KV1(K)/2 + 1
               IOM = KV4(K)
               IF (MOD((KMAX-IOM)/2,2) .EQ. 1) IOM = -IOM
               L = (KMAX-IOM)/4 + 1
               DO 53 NY = 1,NU
                  X = EIV(L,NY,J)*EIGVEC(K,NY)
                  SSTRIP(N) = SSTRIP(N) + UPAIR(NY)*X
                  SPIKUP(N) = SPIKUP(N) + VPAIR(NY)*X
   53          CONTINUE
   52       CONTINUE
            write(6,401)IM(J),F(J),(SSTRIP(N),N=1,6), (SPIKUP(N),N=1,6)
   50    CONTINUE
  401    FORMAT(1X,'I=',I2,'/2 E=',F6.0,' KEV     ',6F15.8,/,
     *   25X,6F15.8,/)
      ENDIF
C
      GAMR=GAMMA*PI/180.
      IF (ISTRCH .EQ. 1) THEN
         A1=1./(1.+EPSP*COS(GAMR)/3.+EPSP*SIN(GAMR)/SQRT(3.))
         A2=1./(1.+EPSP*COS(GAMR)/3.-EPSP*SIN(GAMR)/SQRT(3.))
         A3=1./(1.-EPSP*COS(GAMR)/1.5)
      ELSEIF (ISTRCH .EQ. 0) THEN
         A1=1./SQRT(1.+2.*EPSP*COS(GAMR)/3.+2.*EPSP*SIN(GAMR)/SQRT(3.))
         A2=1./SQRT(1.+2.*EPSP*COS(GAMR)/3.-2.*EPSP*SIN(GAMR)/SQRT(3.))
         A3=1./SQRT(1.-4.*EPSP*COS(GAMR)/3.)
      ENDIF
C
CC    EFAC=41./(FLOAT(AA)**(1./3.))*(1.-FLOAT(AA-2*Z)/FLOAT(AA)/3.)*OMOM
      EFAC=41./(FLOAT(AA)**(1./3.))*EFAC1*OMOM
      R2FAC=41.5/EFAC
      Q00=1.44*(FLOAT(AA)/(A1*A2*A3))**(2./3.)/5.
C      IF (2*Z .LT. AA) THEN
C         Q00=Q00*FLOAT(Z)
C         QFAC = R2FAC
C      ELSE
C **** B(E2) AND Q FROM PROTON DISTRIB ALSO FOR ODD NUETRON ****
C         Q00=Q00*FLOAT(AA-Z)
C         QFAC = 41.5*(FLOAT(AA)**(1./3.))/
C     $   (41.* (1.+FLOAT(AA-2*Z)/FLOAT(AA)/3.)*OMOM)
C      ENDIF
         Q00=Q00*FLOAT(Z)
         QFAC = R2FAC
C     PRINT*, 'R2B=',R2B,'R2FAC=',R2FAC,'QFAC=',QFAC
   79 CONTINUE
C******************* MACROSCOPIC QUADRUPOLE MOMENTS ******
      Q0=Q00*(2.*A3**2-A1**2-A2**2)
      Q2= Q00*SQRT(1.5)*(A1**2-A2**2)
      IF (IQ.EQ.0) GO TO 226
C******************* MICROSCOPIC QUADRUPOLE MOMENTS ******
      Q0=0.0
      Q2=0.0
      ANTNUC=0.
      DO 225 IJ=1,IEMAX
         Q0=Q0+2.0*Q20(IJ)*V2(IJ)*QFAC
         Q2=Q2+2.0*Q22(IJ)*V2(IJ)*QFAC
         ANTNUC=ANTNUC+2.0*V2(IJ)
  225 CONTINUE
C**** CORRECTION IF PARTICLE NUMBER NOT CORRECT, I.E. IF DELTA IS  *****
C**** INPUT AND BCS-EQUATION NOT SOLVED                            *****
      Q0 = Q0*FLOAT(NODD)/ANTNUC
      Q2 = Q2*FLOAT(NODD)/ANTNUC
  226 CONTINUE
      IF (IQ .EQ. 0) write(6,2002) Q0, Q2
      IF (IQ .EQ. 1) write(6,2003) Q0, Q2
 2002 FORMAT(/,' MACR QUADR MOMENTS:  Q0=',F8.3,'  Q2=',F8.3,' FM**2')
 2003 FORMAT(/,' MICR QUADR MOMENTS:  Q0=',F8.3,'  Q2=',F8.3,' FM**2')

      write(6,1000)
 1000 FORMAT(//'   B(E2) AND Q(E2) IN E2B2 AND EB,   B(M1) AND M(M1) IN
     *NM2 AND NM,   T(E2) AND T(M1) IN 1/PS',/)
      write(6,450)
  450 FORMAT(1X,'      E1     II      EF     IF     B(E2)',13X,
     *'B(M1)',14X,2HGD,8X,1HD,7X,'EREL    T(E2)      T(M1)',/,
     *43X,'  WU               MWU',/)
C
C**** STORING OF CLEBSCH-GORDAN COEFFICIENTS
      IIMIN=(IMIN+1)/2
      IIMAX=(ISPIN+1)/2
      DO 81 I=IIMIN,IIMAX
         II=2*I-1
         IIB=I-(IMIN-1)/2
         DO 82 J=1,5
            JF = II- 6 + 2*J
            IF (JF .LT. 0) GO TO 82
            IF (JF .GT. ISPIN) GO TO 82
            DO 83 L=1,IDIM1
               K= K3(L)
               IF (IABS(K) .GT. II) GO TO 83
               C1(IIB,J,L) = CLEBI(II,2,JF,K,0,K)
               C2(IIB,J,L) = CLEBI(II,2,JF,K,-2,K-2)
               C3(IIB,J,L) = CLEBI(II,2,JF,K,2,K+2)
               C4(IIB,J,L) = CLEBI(II,4,JF,K,0,K)
               C5(IIB,J,L) = CLEBI(II,4,JF,K,-4,K-4)
               C6(IIB,J,L) = CLEBI(II,4,JF,K,4,K+4)
   83       CONTINUE
   82    CONTINUE
   81 CONTINUE
      IJTOT = 0
C
C********************  LOOP OF ENERGY LEVELS OF ODD NUCLEUS
      DO 13 I=1,NI
         II=IM(I)
         IIA = (II+1)/2
         IIB = IIA-(IMIN-1)/2
         SUMI(I) = 0.
         SUMSR0(I) = 0.
         SUMSR2(I) = 0.
         SUMSR4(I) = 0.
         SUMDR2(I) = 0.
         SUMDR4(I) = 0.
         SUMJR0(I) = 0.
         SUMJR2(I) = 0.
         SUMJR4(I) = 0.
C********************  SECOND LOOP OF ENERGY LEVELS
      DO 13 J=1,NI
         IF=IM(J)
         IFA = (IF+6-II)/2
         FAS1 = (-1.)**((IF+1)/2)
         DO 21 G=1,NYIC
         DO 21 H=1,NYIC
            SUM1(H,G)=0.0
            SUM2(H,G)=0.0
            SUM3(H,G)=0.0
            SUM11(H,G)=0.0
            SUM12(H,G)=0.0
            SUM13(H,G)=0.0
            SUM14(H,G)=0.0
   21    CONTINUE
         IF(IABS(II-IF).GT.4) GO TO 13
C        CUTOFF IN B(E2) AND B(M1) PRINT-OUT.
         IF(F(I).GT.CUTOFF.OR.F(J).GT.CUTOFF.OR.F(I).LT.F(J)-0.0001)
     &   GO TO 13
         IJTOT = IJTOT+1
         RM1=0.0
         RE2=0.0
C************************  TWO LOOPS OVER S.P. ORBITALS
         DO 12 I2=1,NYIC
         DO 12 I1=1,NYIC
            UUVV=UPAIR(I1)*UPAIR(I2) - VPAIR(I1)*VPAIR(I2)
C************************  TWO LOOPS OVER SPIN PROJECTION, K
            DO 114 L=1,IDIM1
               K=K3(L)
               IF (IABS(K) .GT. II) GO TO 114
            DO 14 LP=1,IDIM1
               KP=K3(LP)
               IF (IABS(KP) .GT. IF) GO TO 14
               IF(IABS(II-IF).GT.2) GO TO 16
               IF(KP.EQ.K) SUM11(I1,I2)=SUM11(I1,I2)+
     *         EIV(L,I1,I)*EIV(LP,I2,J)
               IF(KP.EQ.K) SUM12(I1,I2)=SUM12(I1,I2)+
     *         EIV(L,I1,I)*EIV(LP,I2,J)*C1(IIB,IFA,L)
               IF(KP.EQ.-K+2) SUM13(I1,I2)=SUM13(I1,I2)+
     *         EIV(L,I1,I)*EIV(LP,I2,J)*C2(IIB,IFA,L)
               IF(KP.EQ.-K-2) SUM14(I1,I2)=SUM14(I1,I2)+
     *         EIV(L,I1,I)*EIV(LP,I2,J)*C3(IIB,IFA,L)
   16          CONTINUE
               IF(KP.EQ.K) SUM1(I1,I2)=SUM1(I1,I2)+
     *         EIV(L,I1,I)*EIV(LP,I2,J)*C4(IIB,IFA,L)
               IF(KP.EQ.K-4) SUM2(I1,I2)=SUM2(I1,I2)+
     *         EIV(L,I1,I)*EIV(LP,I2,J)*C5(IIB,IFA,L)
               IF(KP.EQ.K+4) SUM3(I1,I2)=SUM3(I1,I2)+
     *         EIV(L,I1,I)*EIV(LP,I2,J)*C6(IIB,IFA,L)
   14       CONTINUE
  114       CONTINUE
            IF(II.EQ.IF) DELTA=1.0
            IF(II.NE.IF) DELTA=0.0
            IF(I1.EQ.I2) DELTB=1.0
            IF(I1.NE.I2) DELTB=0.0

            DRM1 = SQRT(0.75/PI)*
     $      (GR*SQRT(II*(II+2.)/4.)*DELTA*DELTB*SUM11(I1,I2)+
     $      SUM12(I1,I2)*((GS-GL)*SNOLL(I1,I2)+(GL-GR)*JZ(I1,I2))+
     $      SUM13(I1,I2)*(-1.0)**((IF-1)/2)/SQRT(2.)*((GS-GL)*
     $      SPLUS(I1,I2)+(GL-GR)*APLB(I1,I2))+
     $      SUM14(I1,I2)*(-1.0)**((IF+1)/2)/SQRT(2.)*
     $      ((GS-GL)*SMINUS(I1,I2)+(GL-GR)*C(I1,I2)))
            RM1 = RM1 + DRM1*PUV(I1,I2)

            RE2=RE2 + SQRT(5./16./PI)*
     $      ((Q0*DELTB+EP*UUVV*Q0T(I1,I2)*QFAC)*SUM1(I1,I2)+
     $      (Q2*DELTB+EP*UUVV*Q2T(I1,I2)*QFAC)*
     $      (SUM2(I1,I2)+SUM3(I1,I2)))
            IF (I .EQ. J .AND. (BS2(1)+BS2(2)) .GT. 0.000001) THEN
C
C******************* CONTRIBUTIONS TO HYPERFINE ANOMALY ***************
C
            SUMI(I)=SUMI(I)+
     $      SQRT(II*(II+2.)/4.)*DELTA*DELTB*SUM11(I1,I2)*PUV(I1,I2)
            SUMSR0(I) = SUMSR0(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*SNOLL(I1,I2) + FAS1/SQRT(2.)*
     $      (-SUM13(I1,I2)*SPLUS(I1,I2)+SUM14(I1,I2)*SMINUS(I1,I2)))
            SUMSR2(I) = SUMSR2(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*R2SZ(I1,I2) + FAS1/SQRT(2.)*
     $      (-SUM13(I1,I2)*R2SPLU(I1,I2)+SUM14(I1,I2)*R2SMIN(I1,I2)))
            SUMSR4(I) = SUMSR4(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*R4SZ(I1,I2) + FAS1/SQRT(2.)*
     $      (-SUM13(I1,I2)*R4SPLU(I1,I2)+SUM14(I1,I2)*R4SMIN(I1,I2)))
            SUMDR2(I) = SUMDR2(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*R2DZ(I1,I2) + FAS1*
     $      (SUM13(I1,I2)*R2DPLU(I1,I2)+SUM14(I1,I2)*R2DMIN(I1,I2)))
            SUMDR4(I) = SUMDR4(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*R4DZ(I1,I2) + FAS1*
     $      (SUM13(I1,I2)*R4DPLU(I1,I2)+SUM14(I1,I2)*R4DMIN(I1,I2)))
            SUMJR0(I) = SUMJR0(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*JZ(I1,I2) + FAS1/SQRT(2.)*
     $      (-SUM13(I1,I2)*APLB(I1,I2)+SUM14(I1,I2)*C(I1,I2)))
            SUMJR2(I) = SUMJR2(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*R2JZ(I1,I2) + FAS1/SQRT(2.)*
     $      (-SUM13(I1,I2)*R2JPLU(I1,I2)+SUM14(I1,I2)*R2JMIN(I1,I2)))
            SUMJR4(I) = SUMJR4(I) +
     $      PUV(I1,I2)*(SUM12(I1,I2)*R4JZ(I1,I2) + FAS1/SQRT(2.)*
     $      (-SUM13(I1,I2)*R4JPLU(I1,I2)+SUM14(I1,I2)*R4JMIN(I1,I2)))
            ENDIF
C
C**********************************************************************
C
   12    CONTINUE
C************************** TRANSITION PROBABILITIES ******************
         IF(IABS(II-IF).GT.2) RM1=0.0
         BM1=RM1**2
         BE2=RE2**2
         BE2=BE2/10000.0
         EST=F(I)-F(J)
         IF(F(I)-F(J))101,102,101
  101    CONTINUE
         TE2=12.2*(EST/1000.0)**5*BE2
         TM1=17.6*(EST/1000.0)**3*BM1
         IF (TE2 .NE. 0. .OR. TM1 .NE. 0.) GO TO 107
         X1 =1.0
         GO TO 102
  107    CONTINUE
         X1=TE2/(TE2+TM1)
         IF(ABS(X1-1.0).LE.0.000001) GO TO 102
         PE1=(SIGN(1.,RM1*RE2))*SQRT(X1/(1.-X1))
         GO TO 105
  102    PE1=0.000
C****************** STATIC MOMENTS ************
C
  105    IF (I.EQ.J) RN1=RM1*SQRT(4.*PI/3.)*CLEBI(II,2,II,II,0,II)
         IF(I.EQ.J) RQS=SQRT(16.*PI/5.)*RE2*CLEBI(II,4,II,II,0,II)/100.0
         BE2WU = BE2/BE2W
         BM1MWU = BM1/BM1MW

         IF (I .NE. J) write(6,501) F(I),IM(I),F(J),IM(J),BE2,BE2WU,BM1,
     *   BM1MWU,X1,PE1,EST,TE2,TM1
  501    FORMAT(1X,F9.1,I4,2H/2,2H -,F7.1,I4,2H/2,F10.4,F8.2,F10.4,F8.2,
     *   2F9.3,F9.1,2E11.3)
         IF(I.EQ.J) write(6,500) F(I),IM(I),F(J),IM(J),RQS,RN1
  500    FORMAT(1X,F9.1,I4,2H/2,2H -,F7.1,I4,2H/2,F10.4,8X,F10.4)
         IF (MOD(IJTOT,5) .EQ. 0) write(6,502)
  502    FORMAT(/)
   13 CONTINUE
C************************* HYPERFINE ANOMALY (AND MAGN MOM) **********
      DO 19 IE=1,2
         IF (BS2(IE) .LT. 0.000001) GO TO 19
         IF (IE .EQ. 1) write(6,703) BS2(IE), BS4(IE)
         IF (IE .EQ. 2) write(6,704) BS2(IE), BS4(IE)
  703 FORMAT(///,' ELECTRONIC S-STATE, B-FACTORS FOR HYPERFINE ANOMALY:'
     $'  ,B2=',F7.4,'  B4=',F7.4,
     $   '   UPPER LINE: GS=GSFREE,  LOWER LINE: GS=0.7*GSFREE',/)
  704 FORMAT(///,' ELECTRONIC P-STATE, B-FACTORS FOR HYPERFINE ANOMALY:'
     $'  ,B2=',F7.4,'  B4=',F7.4,
     $   '   UPPER LINE: GS=GSFREE,  LOWER LINE: GS=0.7*GSFREE',/)
         write(6,705)
  705    FORMAT(18X,'** CONTRIBUTION TO MAGN MOM **   ',10(1H*),
     $   ' 100*(CONTRIBUTION TO HYPERFINE ANOMALY) ',10(1H*),/)
         write(6,702)
  702    FORMAT('    ENERGY    I       R       S       L    TOTAL',
     $   '      R     SR2     SR4     LR2     LR4    ',
     $   ' DR2     DR4    TOTAL    EPS(%)',/)
         BL2(IE) = 0.6*BS2(IE)
         BL4(IE) = (3./7.)*BS4(IE)
         BSR2 = BS2(IE)/(R2B/R2FAC)
         BSR4 = BS4(IE)/(R2B/R2FAC)**2
         BLR2 = BL2(IE)/(R2B/R2FAC)
         BLR4 = BL4(IE)/(R2B/R2FAC)**2
         BDR2 = 0.4*BS2(IE)/(R2B/R2FAC)
         BDR4 = (4./7.)*BS4(IE)/(R2B/R2FAC)**2
         DO 18 I=1,NI
            IF (F(I) .GT. CUTOFF) GO TO 18
            II = IM(I)
            CGCOFF = CLEBI(II,2,II,II,0,II)
C
            MAGMR = (SUMI(I)-SUMJR0(I))*GR*CGCOFF
            MAGMS1 = SUMSR0(I)*GSFREE*CGCOFF
            MAGMS2 = 0.7*MAGMS1
            MAGML  = (SUMJR0(I)-SUMSR0(I))*GL*CGCOFF
            MAGM1 = MAGMR + MAGMS1 + MAGML
            MAGM2 = MAGMR + MAGMS2 + MAGML
C
            EPSR = (SUMI(I)-SUMJR0(I))
     $      *GR*(BL2(IE)*ROTNY1+BL4(IE)*ROTNY2)*CGCOFF
            EP1SR2 = SUMSR2(I)*BSR2*GSFREE*CGCOFF
            EP1SR4 = SUMSR4(I)*BSR4*GSFREE*CGCOFF
            EP2SR2 = 0.7*EP1SR2
            EP2SR4 = 0.7*EP1SR4
            EPSLR2 = (SUMJR2(I)-SUMSR2(I))*GL*BLR2*CGCOFF
            EPSLR4 = (SUMJR4(I)-SUMSR4(I))*GL*BLR4*CGCOFF
            EP1DR2 = -SUMDR2(I)*GSFREE*BDR2*CGCOFF
            EP1DR4 = -SUMDR4(I)*GSFREE*BDR4*CGCOFF
            EP2DR2 = 0.7*EP1DR2
            EP2DR4 = 0.7*EP1DR4
            EPS1=EPSR+EP1SR2+EP1SR4+EPSLR2+EPSLR4+EP1DR2+EP1DR4
            EPS2=EPSR+EP2SR2+EP2SR4+EPSLR2+EPSLR4+EP2DR2+EP2DR4
            HYPAN1 = -EPS1/MAGM1
            HYPAN2 = -EPS2/MAGM2
C
            write(6,700) F(I),IM(I),MAGMR,MAGMS1,MAGML,MAGM1,
     $      EPSR,EP1SR2,EP1SR4,EPSLR2,EPSLR4,EP1DR2,EP1DR4,EPS1,HYPAN1
            write(6,701) MAGMS2,MAGM2,EP2SR2,EP2SR4,EP2DR2,EP2DR4,
     $      EPS2,HYPAN2
  700       FORMAT(F10.2,I4,'/2',12F8.3,F10.3)
  701       FORMAT(24X,F8.3,8X,F8.3,8X,2F8.3,16X,3F8.3,F10.3)
   18    CONTINUE
   19    CONTINUE
C
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,600) TID
  600 FORMAT(/,1X,'TIME=',F5.1,' SEC')
  111 CONTINUE
      STOP
      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)
      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,134
      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
      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
C
C#######################################################################
C   APPROPRIATE FOR NORD COMPUTER
C
CC      FUNCTION MSLEFT(XX)
CC      DOUBLE INTEGER IUSE,TUSED
CC      MSLEFT=-IUSE*20
CC      RETURN
CC      END
C
C***********************************************************************
C    APPROPRIATE FOR VAX (WRITTEN BY 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
