CFTN77,Q,T,S
C$CDS ON
C      INCLUDE 'CALPGM.INC'
C   PACKAGE FOR DOUBLET PI : DPI.FOR
C
C   Copyright (C) 1989, California Institute of Technology
C   All rights reserved.  U. S. Government Sponsorship under
C   NASA Contract NAS7-918 is acknowledged
C
C   Herbert M. Pickett  20 Mar 1989
C
      SUBROUTINE HAMX(IBLK,NSIZD,IDPAR,PAR,EGY,T,DEDP,IFDUMP)
C
C .. PACKAGE FOR DOUBLET PI
C
      INTEGER   IBLK,NSIZD,IFDUMP,NSIZE,ITYMAX,IDM,IERR
      REAL*8    IDPAR(0:10),PAR(*),EGY(*),T(NSIZD,*),DEDP(NSIZD,*)
      PARAMETER (ITYMAX=26,IDM=8)
      INTEGER*2 IP(ITYMAX),ISBLK(IDM+1),IWORK(IDM)
      REAL*8    W(ITYMAX,IDM,IDM),WORK(IDM),ZV(2)
      INTEGER   IWHOLE,ISDGN,AII,F,AJI,AJF,I0,I1,I2,L,K1,K3,JQ(IDM),
     &          NPAR,ISIG,I,J,K,IDGN,JS,IQ,IPAR,II,KK,ITMP,KQ,JJ,I4,
     &          JF,OFDIAG,KPAR
      PARAMETER (I0=0,I1=1,I2=2)
      REAL*8  PTY,XJQ,RX,RY,SQJ,X,Z,AA,BB,AB,RT2,AI,FG,
     &          DJ1,DJ2,G,GG,C6JJ,C3JJ,RT,Q,QQ,DDOT
      COMMON /SQN/IWHOLE,ISDGN,AII
C      EMA IDPAR,PAR,EGY,DEDP,T,WORK,W,IWORK
      SAVE NPAR,ISIG,IP,/SQN/
      SAVE W,ZV
	DATA ZV/0.,0./
      NSIZE=NSIZD
      IF(IBLK.EQ.0) THEN
          I=MIN(ISDGN+ISDGN,IDM)
          IF(NSIZE.GT.I) NSIZE=I
          PAR(1)=NSIZE
C         UNSCRAMBLE PARAMETER LIST
          NPAR=IDPAR(0)
          DO 4 I=1,ITYMAX
              DO 5 J=1,NPAR
                  K=ABS(IDPAR(J))
                  IF(K.EQ.I) THEN
                      IP(I)=J
                      GO TO 4
                  ENDIF
5                 CONTINUE
              IP(I)=0
4             CONTINUE
          ISIG=1
          IF(IP(1).NE.0) ISIG=0
          RETURN
      ENDIF
      IF(ISIG.EQ.0) THEN
          K=IP(1)
          ISIG=1
          IF(PAR(K).GT.0) ISIG=-1
      ENDIF
      CALL GETQN(IBLK,I1,JQ,IDGN,I)
      F=IDGN-1
      JS=JQ(1)
      I =JQ(3)
      PTY=I
      IQ=ISIG*I*(1-2*MOD(JS,2))
      DO 1 I=1,NSIZE,2
C ....... SET UP J QUANTUM NUMBER
          JQ(I)=JS*IQ
          JQ(I+1)=JS
          IQ=-IQ
1         JS=JS+1
C .. ZERO HAMILTONIAN AND DERIVATIVE W
      DO 3 I=1,NSIZE
      DO 3 J=I,NSIZE
      DO 3 IPAR=1,NPAR
3         W(IPAR,I,J)=0.
C ... SET UP FINE STRUCTURE TERMS
      DO 6 I=1,NSIZE,2
          II=I+1
          XJQ=JQ(I)+1
          RX=JQ(I)*JQ(I)
          RY=RX-2.
          SQJ=RX-0.25
          X=RX-1.
          Z=SQRT(X)
          DO 7 KK=1,10
              K=IP(KK)
              IF(K.EQ.0) GO TO 7
              GO TO (10,20,30,40,50,60,70,80,90,100),KK
C..A
10            AA=-(1+ISIG)/2
              BB=1+AA
              OFDIAG=0
              GO TO 8
C..AJ
20            AA=-RX
              BB=RY
              OFDIAG=0
              GO TO 8
C..AH
30            AA=RX*RX+X
              BB=-RY*RY-X
              AB=0.
              OFDIAG=1
              GO TO 8
C..B+0.5*Q
40            AA=RX
              BB=RY
              AB=Z
              OFDIAG=1
              GO TO 8
C..D
50            AA=-RX*RX-X
              BB=-RY*RY-X
              AB=-2.*Z*X
              OFDIAG=1
              GO TO 8
C..H
60            AA=(RX*RX+X)*RX+2.*X*X
              BB=(RY*RY+X)*RY+2.*X*X
              AB=((3.*X+1.)*X-1.)*Z
              OFDIAG=1
              GO TO 8
C..P
70            AA=0.5*XJQ
              BB=0.
              AB=0.25*Z
              OFDIAG=1
              GO TO 8
C..Q
80            AA=0.5+JQ(I)
              BB=0.5
              AB=0.5*Z*JQ(I)
              OFDIAG=1
              GO TO 8
C..PD
90            AA=0.5*XJQ*SQJ
              BB=0.
              AB=0.5*Z*SQJ
              OFDIAG=1
              GO TO 8
C..QD
100           AA=0.5*XJQ*XJQ*SQJ
              BB=0.5*X*SQJ
              AB=0.5*Z*XJQ*SQJ
              OFDIAG=1
C
8         IF(OFDIAG.NE.0) W(K,I,II)=AB
          W(K,I,I)=AA
          W(K,II,II)=BB
7         CONTINUE
6         CONTINUE
C
C     CHI PARAMETERS (MAGNETIC INTERACTIONS)
      RT2=SQRT(2.D0)
      AI=0.5*AII
      FG=AI*(AI+1.)*(AI+AI+1.)
      DO 101 I=1,NSIZE,2
          AJI=JQ(I+1)*2 -1
          DJ1=AJI+1
          DO 102 J=I,NSIZE,2
              AJF=JQ(J+1)*2-1
              DJ2=AJF+1
              IF(ABS(AJF-AJI).GT.2) GO TO 102
              G=SQRT(DJ1*DJ2*FG)*C6JJ(AJI,I2,AJF,AII,F,AII)
              ITMP=(AJI+F+AII)/2
              IF(MOD(ITMP,I2).NE.0) G=-G
              GG=G
              ITMP=AJF/2
              IF(MOD(ITMP,I2).NE.0) GG=-GG
              K=IP(11)
              K1=-1
              IF(K.NE.0)W(K,I,J)=GG*2.*C3JJ(AJI,I2,AJF,K1,I0,I1)
              K=IP(12)
              IF(K.NE.0)W(K,I,J)=G*PTY*RT2*C3JJ(AJI,I2,AJF,K1,I2,K1)
              K=IP(13)
              K3=3
              K1=-K3
              IF(K.NE.0)W(K,I+1,J+1)=-C3JJ(AJI,I2,AJF,K1,I0,K3)
     &                               *GG/1.5
              K=IP(14)
              IF(K.EQ.0) GO TO 102
              W(K,I,J+1)=GG*RT2*C3JJ(AJF,I2,AJI,K1,I2,I1)
              IF(J.EQ.I) GO TO 102
              W(K,I+1,J)=GG*RT2*C3JJ(AJI,I2,AJF,K1,I2,I1)
102           CONTINUE
101       CONTINUE
      DO 103 I=1,NSIZE,2
          Z=JQ(I)*JQ(I)-1
          DO 104 KQ=15,19
              K=IP(KQ)
              IF(K.EQ.0) GO TO 104
              IF(KQ.NE.19) THEN
                  AA=Z
                  KK=KQ-4
              ELSE
                  AA=-JQ(I)
                  KK=14
              ENDIF
              J=IP(KK)
              IF(J.EQ.0) GO TO 104
              JJ=KQ-14
              GO TO (106,106,107,108,108),JJ
106               W(K,I,I)=AA*W(J,I,I)
                  IF(JJ.NE.1) GO TO 104
                  J=IP(13)
                  IF(J.EQ.0) GO TO 104
C
107               W(K,I+1,I+1)=AA*W(J,I+1,I+1)
                  GO TO 104
C
108               W(K,I,I+1)=AA*W(J,I,I+1)
104           CONTINUE
103       CONTINUE
C
C     ETA PARAMETERS (QUADRUPOLE INTERACTIONS)
      IF(AII.LT.2)  GO TO 200
      FG=FG*(AI+1.5)/(AI*AI*(AI-0.5))
      RT=PTY/SQRT(24.D0)
      L=4
      DO 111 I=1,NSIZE,2
          AJI=JQ(I+1)*2-1
          DJ1=AJI+1
          DO 112 J=I,NSIZE,2
              AJF=JQ(J+1)*2-1
              DJ2=AJF+1
              IF(ABS(AJF-AJI).GT.L) GO TO 112
              Q=SQRT(DJ1*DJ2*FG)*C6JJ(AJI,L,AJF,AII,F,AII)
              ITMP=(AJI+F+AII)/2
              IF(MOD(ITMP,I2).NE.0) Q=-Q
              QQ=Q
              ITMP=AJF/2
              IF(MOD(ITMP,I2).NE.0) QQ=-QQ
              K=IP(20)
              K1=-1
              IF(K.NE.0)W(K,I,J)=0.25*QQ*C3JJ(AJI,L,AJF,K1,I0,I1)
              K=IP(21)
              K3=3
              K1=-K3
              IF(K.NE.0)W(K,I+1,J+1)=-0.25*QQ*C3JJ(AJI,L,AJF,K1,I0,K3)
              K=IP(22)
              IF(K.EQ.0) GO TO 112
              K3=-1
              W(K,I,J+1)=-RT*Q*C3JJ(AJI,L,AJF,K3,L,K1)
              IF(I.EQ.J) GO TO 112
              I4=4
              W(K,I+1,J)=RT*Q*C3JJ(AJI,I4,AJF,K1,L,K3)
112           CONTINUE
111       CONTINUE
      DO 113 I=1,NSIZE,2
          Z=JQ(I)*JQ(I)-1
          SQJ=Z+0.75
          DO 118 KQ=23,26
              K=IP(KQ)
              IF(K.EQ.0) GO TO 118
              KK=KQ-22
              GO TO (114,115,116,117),KK
114               AA=SQJ
                  BB=SQJ
                  AB=SQJ/JQ(I)
                  OFDIAG=1
                  GO TO 119
115               AA=JQ(I)
                  BB=AA*Z/(Z-6.)
                  OFDIAG=0
                  GO TO 119
116               AA=0.
                  BB=0.
                  AB=-1./JQ(I)
                  OFDIAG=1
                  GO TO 119
117               AA=0.
                  BB=Z/(6.-Z)
                  OFDIAG=0
                  GO TO 119
119           IF(OFDIAG.NE.0) THEN
                KK=IP(22)
                W(K,I,I+1)=AB*W(KK,I,I+1)
              ENDIF
              KK=IP(21)
              W(K,I+1,I+1)=BB*W(KK,I+1,I+1)
              KK=IP(20)
              W(K,I,I)=AA*W(KK,I,I)
118           CONTINUE
113       CONTINUE
C
C  BRING H TOGETHER
C
200   DO 201 I=1,NSIZE
          ISBLK(I)=I
          DO 201 J=I,NSIZE
201           T(J,I)=DDOT(NPAR,W(1,I,J),I1,PAR,I1)
      IF(IFDUMP.NE.0) THEN
C        special to dump Hamiltonian 
          DO 300 I=1,NSIZE
              EGY(I)=T(I,I)
              CALL DCOPY(NPAR,ZV,I0,DEDP,NSIZE)
 300          CONTINUE          
          RETURN
      ENDIF    
C
C DIAGONALIZE AND ORDER EIGENVECTORS
C
      ISBLK(NSIZE+1)=NSIZE+1
      CALL HDIAG(NSIZD,NSIZE,T,EGY,WORK,IERR)
      CALL ORDBLK(NSIZD,NSIZE,T,EGY,NSIZE,ISBLK,WORK,IWORK,I0)
C
C  TRANSFORM W TO DEDP
      KPAR=0
      DO 350 IPAR=1,NPAR
          IF(IDPAR(IPAR).GE.0) THEN
            KPAR=KPAR+1
            KK=IPAR
          ELSE
            AB=PAR(IPAR)/PAR(KK)  
          ENDIF  
          DO 351 I=1,NSIZE
              BB=0.
              DO 352 K=1,NSIZE
                  AA=0.
                  JF=K-1
                  DO 354 J=1,JF
354                   AA=AA+W(IPAR,J,K)*T(J,I)
                  DO 355 J=K,NSIZE
355                   AA=AA+W(IPAR,K,J)*T(J,I)
352               BB=BB+AA*T(K,I)
              IF(IPAR.EQ.KK) THEN
                  DEDP(I,KPAR)=BB 
              ELSE
                  DEDP(I,KPAR)=DEDP(I,KPAR)+BB*AB
              ENDIF
351           CONTINUE                  
350       CONTINUE
      RETURN
      END
C
      FUNCTION INTENS(IBLK,ISIZ,JBLK,JSIZ,IDIP,DIP,S,NSIZD)
C
C ... PACKAGE FOR DOUBLET PI ...ELECTRIC DIPOLE
C
      INTEGER   INTENS,IBLK,ISIZ,JBLK,JSIZ,NSIZD
      INTEGER*4 IDIP(0:*)
      REAL*8 S(NSIZD,*),DIP(*)
C      EMA IDIP,DIP,S
      INTEGER IWHOLE,ISDGN,AII,JJI,JJJ,FI,FJ,I2,
     &        IDGN,JDGN,IX,JX,I,J,K,IM
      REAL*8 DGN,FAC,AJI,OMEG,C6JJ,AJJ
      COMMON /SQN/IWHOLE,ISDGN,AII
      SAVE /SQN/,I2
      DATA I2/ 2/
      INTENS=0
C START OF INTENSITY CALCULATION
      IF(IDIP(0).EQ.0) RETURN
      IF(MOD(IBLK+JBLK,2).EQ.0) RETURN
      IDGN=((IBLK-1)/2)*2+IWHOLE
      JDGN=((JBLK-1)/2)*2+IWHOLE
      IF(ABS(IDGN-JDGN).GT.2) RETURN
      FI=IDGN-1
      FJ=JDGN-1
      DGN=IDGN*JDGN
      IX=ABS(IDGN-ISDGN)
      JX=ABS(JDGN-ISDGN)-IX+2
      FAC=DIP(1)*SQRT(DGN)
      DO 101 I=1,ISIZ
          IM=I-1
          K=MOD(IM,I2)
          JJI=IX+IM-K
          AJI=0.5*JJI
          IF(K.EQ.0) FAC=-FAC
          OMEG=K+0.5
          DO 102 J=1,JSIZ
102           S(I,J)=0.
          J=I-JX
          IF(J.GT.JSIZ) GO TO 101
          IF(J.GT.0) THEN
              JJJ=JJI-2
              S(I,J)=FAC*SQRT((AJI-OMEG)*(AJI+OMEG)/AJI)*
     1               C6JJ(JJI,I2,JJJ,FJ,AII,FI)
          ENDIF
          J=J+2
          IF(J.GT.JSIZ) GO TO 101
          IF(J.GT.0.AND.JJI.GT.K) THEN
              S(I,J)=FAC*OMEG*SQRT((AJI+AJI+1.)/(AJI*AJI+AJI))*
     1               C6JJ(JJI,I2,JJI,FJ,AII,FI)
          ENDIF
          J=J+2
          IF(J.GT.0.AND.J.LE.JSIZ) THEN
              JJJ=JJI+2
              AJJ=AJI+1.
              S(I,J)=-FAC*SQRT((AJJ-OMEG)*(AJJ+OMEG)/AJJ)*
     1               C6JJ(JJI,I2,JJJ,FJ,AII,FI)
          ENDIF
101       CONTINUE
      INTENS=IDGN
      RETURN
      END
C      
      SUBROUTINE SETINT(IFDIAG,NSAV,IDIP)
C         DOES NOTHING TO INITIALIZE
      INTEGER IFDIAG,NSAV,NDIP
      INTEGER*4 IDIP(0:*)
C      EMA IDIP
      IFDIAG=1
      NSAV=1
      NDIP=IDIP(0)
      RETURN
      END
      SUBROUTINE GETQN(IBLK,INDX,IQN,IDGN,IQ)
C
C ..PACKAGE FOR DOUBLET PI
C
      INTEGER   IBLK,INDX,IDGN,IQ,IQN(6)
      INTEGER IWHOLE,ISDGN,AI,KF,MOD
      COMMON /SQN/IWHOLE,ISDGN,AI
      SAVE /SQN/
      KF=(IBLK-1)/2
      IDGN=KF+KF+IWHOLE
      IF(INDX.EQ.0) THEN
          IDGN=2*MIN(ISDGN,IDGN)
          RETURN
      ENDIF
      IQ=1
      IQN(1)=ABS(IDGN-ISDGN)/2+(INDX+1)/2
      IQN(2)=MOD(INDX-1,2)+1
      IQN(3)=1-2*MOD(IBLK,2)
      IQN(4)=KF
      IF(IQN(1).LT.IQN(2)) IDGN=0
      RETURN
      END
C
      FUNCTION SETOPT(LUIN,IQNFMT,ITD,NAMFIL)
C      SET OPTIONS
      INTEGER      SETOPT,LUIN,IQNFMT,ITD
      CHARACTER*80 NAMFIL,CARD
      INTEGER IWHOLE,ISDGN,AII,I,PCARD
      REAL*8 VAL(1)
      COMMON /SQN/IWHOLE,ISDGN,AII
      SAVE /SQN/
      CALL CJJINI
      READ(LUIN,'(A)',END=999) CARD
      VAL(1)=1
      IWHOLE=1
      I=PCARD(CARD,VAL,IWHOLE)
      ISDGN=VAL(1)
      IF(ISDGN.LT.1) ISDGN=1
      NAMFIL='dpi.nam'
      ITD=2
      AII=ISDGN-1
      IWHOLE=MOD(AII,2)
      IQNFMT=804
      IF(IWHOLE.EQ.0) IQNFMT=IQNFMT+10
      IF(ISDGN.EQ.1)  IQNFMT=803
      SETOPT=1
      RETURN
 999  SETOPT=-1
      RETURN
      END
      FUNCTION SETBLK(LU,IDPAR,PAR,NBLKPF,NEGY)
C      SET OPTIONS
      INTEGER SETBLK,LU,NBLKPF,NEGY
      REAL*8 IDPAR(*),PAR(*),DMY(2)
C      EMA IDPAR,DMY
      REAL*4 AI
      INTEGER ITMP,IWHOLE,ISDGN,AII
      COMMON /SQN/IWHOLE,ISDGN,AII
      SAVE /SQN/
      NBLKPF=2
      AI=0.5*AII
      WRITE(LU,900) AI
900   FORMAT(' doublet PI lines for I =',F5.1)
      ITMP=0
      CALL HAMX(ITMP,NEGY,IDPAR,DMY,PAR,PAR,PAR,ITMP)
      NEGY=DMY(1)
      SETBLK=1
      RETURN
      END