      SUBROUTINE BAS9IN(PRTP,IBOUND,IPRINT)
      USE efvs, ONLY: EFV, EFVNAM, EFVUNT, MAPEFV, NEFV
      USE potential, ONLY: IREF, NCONST, NRSQ, VCONST
      USE basis_data, ONLY: ELEVEL, JHALF, JLEVEL, NLEVEL, ROTI
      USE physical_constants, ONLY: Debye_in_SI, Planck_constant_in_SI,
     c                              bohr_magneton, g_e,
     c                              hartree_in_inv_cm,
     c                              inverse_fine_structure_constant,
     c                              kilo_in_SI, nuclear_magneton,
     c                              speed_of_light_in_cm
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE
C-----------------------------------------------------------------------
C  DATE (LAST UPDATE): 11/06/13        STATUS: FINISHED                |
C  AUTHOR: MAYKEL LEONARDO GONZALEZ MARTINEZ                           |
C-----------------------------------------------------------------------
C    THIS SUBROUTINE & ITS ENTRIES SET THE BASIS SET (I.E. QUANTUM     |
C    NUMBERS, COUPLING MATRIX ELEMENTS...) FOR THE INTERACTION TYPES:  |
C                  ATOM(S) + ATOM(L) IN MAGNETIC FIELDS                |
C                          USING THE BASIS SET                         |
C          |i1 ms1> |s1 ms1> |i2 mi2> |(l2 s2)j2 mj2> |L ML>           |
C   HYPERFINE TERMS CAN BE INCLUDED/OMITTED WITH HFFLG = TRUE/FALSE    |
C    * ALL MONOMER QUANTUM NUMBERS & MTOT ARE TWICE THEIR VALUES *     |
C-----------------------------------------------------------------------
C
C VARIOUS UPDATES BY DGG (2015) AND MDF (2019)
C USES
C
      PARAMETER (L2MAX=7)
C EXTERNAL
      EXTERNAL dS1L2d,dS1S2d,dM2S1d
      DIMENSION IXS1L(7),IXSK(6)
C 'EXTRA' VARIABLES:
      CHARACTER(8) PRTP(4),QNAME(10),CL2(L2MAX)
      LOGICAL LEVIN,EIN,LCOUNT
      DIMENSION JSTATE(1),VL(1),IV(1),JSINDX(1),L(1),CENT(1),LAM(1)
      DIMENSION DGVL(*)
      DIMENSION MONQN(NQN2)
C
      CHARACTER A2SP1,B2SP1
      LOGICAL ODD,DIAG2,DIAG3,D1C1R,DHF12,D2C2R,M2S1
c---
c DGG BUGCHKS:
      LOGICAL CALCZEEMAN_S1, CALCZEEMAN_I1
      LOGICAL CALCZEEMAN_L2, CALCZEEMAN_S2
      LOGICAL CALCHF1, CALCSO2, CALCPOT
      LOGICAL CALC_VK0S1,CALC_VK0S2
      LOGICAL CALC_VK2S1,CALC_VK2S2
      LOGICAL WRITEVLMATRICES

c----
      LOGICAL HFFLG,FLGJ2
      DIMENSION AJ(3,3),BJ(3,3)
      DIMENSION IS(10)
C COMMONS
      COMMON /VSO/ A2
C EQUIVALENCES
      EQUIVALENCE (AC2,BF2)
C NAMELISTS
      NAMELIST /BASIS9/IS1,GS1,II1,XN1,BF1,EQQ1,
     &     L2,IS2,GS2,GL2,A2,II2,XN2,AJ,BJ,
     &     BF2,
     &     AC2,AL2,AD2,EQQ2,
     &     PHI02,RL3AV2,RD3AV2,
     &     IS,LMAX,HFFLG,FLGJ2,M2S1,
     &     EZERO,J2ZERO,
     $     CALCZEEMAN_S1, CALCZEEMAN_I1, !DGG BUGCHKS
     $     CALCZEEMAN_L2, CALCZEEMAN_S2,
     $     CALCHF1, CALCSO2, CALCPOT,
     $     CALC_VK0S1,CALC_VK0S2,
     $     CALC_VK2S1,CALC_VK2S2,
     $     WRITEVLMATRICES


! FLGJ2 added by DGG to restrict basis set to J2=J2MX
C USEFUL CONSTANTS (FROM NIST 2006 FUNDAMENTAL PHYSICAL CONSTANTS)
C ELECTRON & ORBITAL G-FACTORS, BOHR & NUCLEAR MAGNETONS (1/(cm G)),
C PLANCK'S CONSTANT (J s), SPEED OF LIGHT IN VACUO (cm/s), DEBYE (C m),
C CONVERSION FACTOR FROM cm^-1 TO kHz,
C INVERSE VALUE OF THE FINE STRUCTURE CONSTANT,
C CONVERSION FACTOR FROM Eh TO cm^-1
C DATA
C CRLS 30/06/11: changed to param so that NIST values can be converted explicitly
C CRLS 30/06/11: values amended to current NIST data
C     parameter (GS=2.00231930436153d0) !(GS=2.00231930436153d0)
C     parameter (alphai=1d0/7.2973525698d-3)
C     DATA GS/2.0023193043622D0/,GLP/1.00107D0/
C     DATA GLP/1.00107D0/ !1.00107D0/
C     DATA BM/0.466864515D-4/,XNM/2.542623616D-8/
C     DATA H/6.6260755D-34/,CCM/2.99792458D10/,DEBYE/3.33564D-30/
C     DATA CMIKHZ/2.99792458D7/,ALPHAI/137.035999679D0/
C     DATA CMIKHZ/2.99792458D7/
C     DATA AEU2CMI/2.194746313705D5/
C
C  20-09-2016: UPDATED TO USE MODULE (physical_constants) THAT CONTAINS
C              CONSISTENT AND UP-TO-DATE VALUES
      PARAMETER (GS=-g_e, ALPHAI=inverse_fine_structure_constant,
     &           BM=bohr_magneton, XNM=nuclear_magneton,
     &           H=Planck_constant_in_SI, CCM=speed_of_light_in_cm,
     &           DEBYE=Debye_in_SI, GLP=1.D0,
     &           CMIKHZ=speed_of_light_in_cm/kilo_in_SI,
     &           AEU2CMI=hartree_in_inv_cm)
      DATA CL2/'S','P','D','F','G','H','I'/
C     BAS9IN IS CALLED ONCE FOR EACH SCATTERING SYSTEM (USUALLY ONCE
C     PER RUN) AND CAN READ IN ANY BASIS SET INFORMATION NOT CONTAINED
C     IN NAMELIST BLOCK &BASIS. IT MUST ALSO HANDLE THE FOLLOWING
C     VARIABLES AND ARRAYS:
C
C     PRTP   SHOULD BE RETURNED AS A CHARACTER STRING DESCRIBING THE
C            INTERACTION TYPE
C     IDENT  CAN BE SET>0 IF AN INTERACTION OF IDENTICAL PARTICLES IS
C            BEING CONSIDERED AND SYMMETRISATION IS REQUIRED.
C            HOWEVER, THIS WOULD REQUIRE EXTRA CODING IN IDPART.
C     IBOUND CAN BE SET>0 IF THE CENTRIFUGAL POTENTIAL IS NOT OF THE
C            FORM L(L+1)/R**2; IF IBOUND>0, THE CENT ARRAY MUST BE
C            RETURNED FROM ENTRY CPL9
C STATEMENT FUNCTIONS:
C (1) TO CHECK WHETHER A NUMBER IS ODD
      ODD(NMBR)=2*(NMBR/2).NE.NMBR
C (2) TO CHECK WHETHER CERTAIN 'DIAGONALITY' CONDITIONS ARE FULFILLED
      DIAG2(J1,K1, J2,K2)=J1.EQ.J2 .AND. K1.EQ.K2
      DIAG3(J1,K1,L1, J2,K2,L2)=J1.EQ.J2 .AND. K1.EQ.K2 .AND. L1.EQ.L2
C
C USES THE 'NEW' MECHANISM IN WAVMAT (08/06) TO HANDLE NON-DIAGONAL HMON
C DEFAULT IS CENTRIFUGAL TERMS FROM CENT ARRAY
      IBOUND=0
      NRSQ=0
C DEFAULT: SPIN-ORBIT + ZEEMAN ARE INCLUDED VIA WAVMAT
      NCONST=2
      MAPEFV=2
CC      EFFLG=ISVEFV.EQ.2 !(THIS VALUE IS CARRIED IN efvs MODULE
C STATES THAT JTOT IS NOT USED FOR WHAT IT IS INTENDED TO (I.E. TOTAL J)
C THIS VERSION USES JTOT FOR |2 x MTOT| (PROJECTION OF TOTAL J)
      JHALF=0
C INITIALISES & READS VARIABLES IN NAMELIST BASIS9
C --- INTERACTION PARTNER 1 (CP1) ---
      IS1=0         ! 2 x ELECTRONIC SPIN
      GS1=GS        ! ELECTRON G-FACTOR
      II1=0         ! 2 x NUCLEAR SPIN
      XN1=0.D0      ! NUCLEAR MAGNETIC MOMENT (IN NUCLEAR MAGNETONS)
      BF1=0.D0      ! HFS: FERMI-CONTACT INTERACTION CONSTANT
      EQQ1=0.D0     ! HFS: NUCLEAR ELECTRIC QUADRUPOLE CONSTANT
C --- INTERACTION PARTNER 2 (CP2) ---
      IS2=0         ! 2 x ELECTRONIC SPIN
      L2=2          ! 2 x ELECTRONIC ORB. ANG. MOMENTUM (P ATOM)
      A2=0.D0       ! SPIN-ORBIT CONSTANT (cm^-1)
      GS2=GS        ! ELECTRON G-FACTOR
      GL2=GLP       ! L G-FACTOR
      II2=0         ! 2 x NUCLEAR SPIN
      XN2=0.D0      ! NUCLEAR MAGNETIC MOMENT (IN NUCLEAR MAGNETONS)
      AJ=0.D0       ! HFS:
      BJ=0.D0       ! HFS:
      AC2=0.D0      ! HFS: FERMI-CONTACT INTERACTION CONSTANT
      AL2=0.D0      ! HFS: ORBITAL COUPLING CONSTANT
      AD2=0.D0      ! HFS: ELECT. SPIN-NUCLEAR SPIN DIPOLAR CONSTANT
      EQQ2=0.D0     ! HFS: NUCLEAR ELECTRIC QUADRUPOLE CONSTANT
      PHI02=0.D0    ! HFS:
      RS3AV=0.D0    ! HFS:
      RL3AV=0.D0    ! HFS:
C --- INTERACTION COMPLEX (CP1 + CP2) ---
      IS(1)=IS1+IS2 ! 2 x TOTAL ELECTRONIC SPIN OF CP1 + CP2
      IS(2)=-1      ! FLAGS MULTIPLE TOTAL SPIN (PESs)
      LMAX=0        ! MAX END-OVER-END ORBITAL ANGULAR MOMENTUM
!DGG BUGCHK
      M2S1=.TRUE. !.FALSE.  ! ONLY .TRUE. IF MULTIPLE TOTAL SPINS ARE CONSIDERED
! END BUGCHK
      HFFLG=.FALSE. ! FLAG ONLY .TRUE. IF HFS TERMS ARE INCLUDED
      FLGJ2=.FALSE. ! IF FLGJ2=1 RESTRICT BASIS TO J2=J2MX, IF FLGJ2=0 THEN NO RESTRICTION
      EZERO=0.D0   ! ZERO OF ENERGY
      J2ZERO=-1     ! IF SET TO .GE.0, EZERO IS SET TO THE ZERO FIELD
                    ! ENERGY OF J2=J2ZERO
C --- DGG BUGCHKS:
      CALCZEEMAN_S1=.TRUE.
      CALCZEEMAN_I1=.TRUE.
      CALCZEEMAN_L2=.TRUE.
      CALCZEEMAN_S2=.TRUE.
      CALCHF1=.TRUE.            ! CALCULATE HYPERFINE ELEMENTS FOR ATOM1 - IF FALSE ELEMENTS WILL BE CALCULATED AS ZERO
      CALCSO2=.TRUE.            ! CALCULATE SPIN-ORBIT ELEMENTS FOR ATOM2 - IF FALSE ELEMENTS WILL BE CALCULATED AS ZERO
      CALCPOT=.TRUE.            ! CALCULATE POTENTIAL MATRIX ELEMENTS. IF FALSE ALL WILL BE ZERO
      CALC_VK0S1=.TRUE.
      CALC_VK0S2=.TRUE.
      CALC_VK2S1=.TRUE.
      CALC_VK2S2=.TRUE.
      WRITEVLMATRICES=.FALSE.
c----

C READS NAMELIST /BASIS9/ & COMPUTES/MODIFIES SEVERAL USEFUL (PRE)FACTORS
      READ(5,BASIS9)

      WRITE(6,*) "HYPERFINE AND J2 FLAGS:", HFFLG, FLGJ2 !DGG BUGCHK
C SETS UP 'INTERACTION TYPE'
      A2SP1=CHAR(IS1+1+48) ! ASSUMES MULTIPLICITY < 10
      B2SP1=CHAR(IS2+1+48) ! ASSUMES MULTIPLICITY < 10
      PRTP(1)=A2SP1//"S + "//B2SP1//CL2(L2/2+1)//" "
      PRTP(2)="IN A MAG"
      PRTP(3)="NETIC FI"
      PRTP(4)="ELD     "
C ERROR: THIS VERSION WORKS FOR L2 <= L2MAX
      IF (L2.GT.2*(L2MAX-1)) THEN
         WRITE(6,107) L2MAX-1
 107     FORMAT(/" *** ERROR *** THIS BASE9 WORKS FOR L2 <= ",I1,
     & ": PLEASE, INCREASE L2MAX AND RECOMPILE")
         STOP
      ENDIF
C WARNING: THIS VERSION DOESN'T WORK WITH L2 != 0, 2 & I2 != 0
      IF (L2.GT.2 .AND. HFFLG .AND. II2.NE.0) THEN
          WRITE(6,108)
 108      FORMAT(/" *** WARNING *** THIS BASE9 ONLY ALLOWS INCLUDING
     & HYPERFINE TERMS FOR S AND P ATOMS: I2 SET TO ZERO")
          II2=0
      ENDIF
C SETS UP MIN/MAX VALUES FOR j2 = l2 + s2
      J2MN=ABS(L2-IS2) ! 2 x MINIMUM j2
      J2MX=L2+IS2      ! 2 x MAXIMUM j2
C CHECKS IF HFS TERMS ARE TO BE INCLUDED
      IF (II1.EQ.0 .AND. II2.EQ.0) HFFLG=.FALSE.
      IF (HFFLG) THEN
         NCONST=NCONST+1
         MAPEFV=MAPEFV+1
         VCONST(2)=1.D0/CMIKHZ
      ELSE
         II1=0
         II2=0
         VCONST(2)=0.D0
      ENDIF
C NUCLEAR AND ELECTRON SPINS, AND ORBITAL ANGULAR MOMENTUM(A)
      XI1=0.5D0*II1
      S1=0.5D0*IS1
      XI2=0.5D0*II2
      S2=0.5D0*IS2
      XL2=0.5D0*L2
C CP1 ELECTRON SPIN-CP2 ELECTRON SPIN PREFACTOR
      IF (IS1.NE.0 .AND. IS2.NE.0) FS1S2=AEU2CMI/ALPHAI**2
C CP1 ELECTRON SPIN-CP2 ORBITAL ANGULAR MOMENTUM PREFACTOR
      IF (IS1.NE.0 .AND. L2.NE.0) FS1L2=(GL2/GS2)*AEU2CMI/ALPHAI**2
C CP2 SPIN-ORBIT CONSTANT
      IF (A2.EQ.0.D0 .AND. ROTI(1).NE.0.D0) A2=ROTI(1)
C CP1 & CP2 ELECTRON ZEEMAN PREFACTORS
C *** NOTE *** - SIGNS ARE BEING INCLUDED IN NUCLEAR ZEEMAN PREFACTORS &
C  1/2 FACTORS TAKE INTO ACCOUNT THAT Q. NUMBERS ARE TWICE THEIR VALUES
C . NUCLEAR SPIN ZEEMAN
      IF (II1.NE.0) THEN
         G1=XN1/XI1
         FZI1=-0.5D0*G1*XNM
      ENDIF
      IF (II2.NE.0) THEN
         G2=XN2/XI2
         FZI2=-0.5D0*G2*XNM
      ENDIF
C . ELECTRON ZEEMAN
      IF (IS1.NE.0) FZS1=0.5D0*GS1*BM
      IF (IS2.NE.0) FZS2=GS2*BM
C . ORBITAL ZEEMAN
      FZL2=GL2*BM
C WRITES OUT SOME GENERAL INFORMATION
      IF (L2.EQ.0) THEN
         WRITE(6,109)
 109     FORMAT(/"  BASIS SET IN GENERAL READS
     & |i1 ms1> |s1 ms1> |i2 mi2> |s2 ms2> |L ML>")
      ELSE
         WRITE(6,110)
 110     FORMAT(/"  BASIS SET IS <(LS)-COUPLED>, WHICH IN GENERAL READS
     & |i1 ms1> |s1 ms1> |i2 mi2> |(l2 s2)j2 mj2> |L ML>")
      ENDIF
      IF (.NOT.HFFLG) THEN
         WRITE(6,111)
 111     FORMAT(/"  HYPERFINE TERMS DONT EXIST OR DELIBERATELY OMITTED")
      ELSE
         WRITE(6,112)
 112     FORMAT(/"  HYPERFINE TERMS INCLUDED")
      ENDIF
      WRITE(6,113)
 113  FORMAT(/"  MSET USED FOR PARITY CASES"
     &       /"  PARITY P = p1p2(-1)^L; p1 = 1, p2 = (-1)^l2")


C VCONST(I = 1-3), IN THE MOST GENERAL CASE:
C I = 1, PREFACTOR SPIN-ORBIT
      VCONST(1)=A2
C I = 2, (IF HFFLG = TRUE) PREFACTOR FOR HYPERFINE TERMS
C     VCONST(2)=CONVERSION FACTOR FROM kHz TO 1/cm
C I = 3, (= 2 IF HFFLG = FALSE), MAGNETIC FIELD IN GAUSS
C     VCONST(MAPEFV)=MAGNETIC (B) FIELD, IT IS SET UP IN DRIVER...
C
C  SET UP ELEMENTS OF efvs MODULE
      NEFV=1
      EFVNAM(1)='MAGNETIC Z FIELD'
      EFVUNT(1)='GAUSS'

      RETURN
C
      ENTRY SET9(LEVIN,EIN,NSTATE,JSTATE,NQN,QNAME,MXPAR,NLABV,IPRINT)
C
C     SET9 IS CALLED ONCE FOR EACH SCATTERING SYSTEM. IT SETS UP:
C     MXPAR, THE NUMBER OF DIFFERENT SYMMETRY BLOCKS
C     NLABV, THE NUMBER OF INDICES NEEDED TO DESCRIBE EACH TERM
C       IN THE POTENTIAL EXPANSION
C     NLEVEL AND JLEVEL, UNLESS LEVIN IS .TRUE.;
C     JSTATE AND NSTATE;
C     ELEVEL, UNLESS EIN IS .TRUE.
C     IF THE LOGICAL VARIABLES ARE .TRUE. ON ENTRY, THE CORRESPONDING
C     QUANTITIES WERE INPUT EXPLICITLY IN NAMELIST BLOCK &BASIS.
C     IF EIN IS .FALSE., THE MOLECULAR CONSTANTS MUST HAVE BEEN SUPPLIED
C     IN THE &BASIS ARRAY ROTI: THE PROGRAMMER MAY USE THESE IN ANY WAY
C     HE LIKES, BUT SHOULD OUTPUT THEM HERE FOR CHECKING.
C     NOTE THAT JLEVEL CONTAINS JUST THE QUANTUM NUMBERS NECESSARY TO
C     SPECIFY THE THRESHOLD ENERGY (AND ELEVEL CONTAINS THE CORRESPONDING
C     ENERGIES) WHEREAS JSTATE CONTAINS ALL THE CHANNEL QUANTUM NUMBERS EXCEPT
C     THE ORBITAL L, WHICH MAY BE A SUPERSET. THE LAST COLUMN OF THE JSTATE
C     ARRAY CONTAINS A POINTER TO THE ENERGY IN THE ELEVEL ARRAY.
C
C *** NOTE *** ALL MONOMER QUANTUM NUMBERS & MTOT ARE TWICE THEIR VALUES
C
      EIN=.FALSE. ! ENERGY LEVELS WILL BE CALCULATED
      MXPAR=2
      NLABV=1
C BASIS SET IS: [|i1 mi1> |s1 ms1> |i2 mi2> |(l2 s2)j2 mj2>] |L ML>
C (OR, IF l2 = 0,: [|i1 mi1> |s1 ms1> |i2 mi2> |s2 ms2>] |L ML>)
C IN PRACTICE, NONE INTERACTION CHANGES e- MOMENTA, ONLY THEIR PROJECTIONS
      NEXT=1
      IF (II1.NE.0) THEN
         QNAME(NEXT)='  2xmi1 '
         NEXT=NEXT+1
      ENDIF
      IF (IS1.NE.0) THEN
         QNAME(NEXT)='  2xms1 '
         NEXT=NEXT+1
      ENDIF
      IF (II2.NE.0) THEN
         QNAME(NEXT)='  2xmi2 '
         NEXT=NEXT+1
      ENDIF
      IF (L2.EQ.0) THEN
         QNAME(NEXT)='  2xms2 '
      ELSE
         QNAME(NEXT)='   2xj2 '
         NEXT=NEXT+1
         QNAME(NEXT)='  2xmj2 '
      ENDIF
C JLEVEL HOLDS ALL POSSIBLE COMBINATIONS OF MONOMER QUANTUM NUMBERS
      NQN=NEXT+1 ! DEFAULT IS THAT ONLY j2 & mj2 (ms2 IF l2=0) ARE INCLUDED
C GETS ACTUAL NUMBER OF (MONOMER) QUANTUM NUMBERS
      NQN1=NQN-1
      NLEVEL=0

      DO 210 MI1=-II1,II1,2
         DO 211 MS1=-IS1,IS1,2
            DO 212 MI2=-II2,II2,2
               DO 213 JJ2=J2MN,J2MX,2
                  IF (FLGJ2 .AND. JJ2.NE.J2MX) THEN
                     CYCLE ! i.e., keep going until correct point in loop
                  ENDIF
                  DO 214 MJ2=-JJ2,JJ2,2
                     ITMP=NQN1*NLEVEL
                     IF (II1.NE.0) THEN
                        JLEVEL(ITMP+1)=MI1
                        ITMP=ITMP+1
                     ENDIF
                     IF (IS1.NE.0) THEN
                        JLEVEL(ITMP+1)=MS1
                        ITMP=ITMP+1
                     ENDIF
                     IF (II2.NE.0) THEN
                        JLEVEL(ITMP+1)=MI2
                        ITMP=ITMP+1
                     ENDIF
                     IF (L2.NE.0) THEN
                        JLEVEL(ITMP+1)=JJ2
                        ITMP=ITMP+1
                     ENDIF
                     JLEVEL(ITMP+1)=MJ2
                     NLEVEL=NLEVEL+1
 214              CONTINUE
 213           CONTINUE
 212        CONTINUE
 211     CONTINUE
 210  CONTINUE
      NSTATE=NLEVEL
C LOOPS OVER LEVELS AGAIN, THIS TIME SETTING UP JSTATE
      II=0
      DO 250 I=1,NLEVEL
         II=II+1
         ITMP=NQN1*(I-1)
         DO 240 KK=0,NQN1-1
            JSTATE(KK*NSTATE+I)=JLEVEL(ITMP+KK+1)
  240    CONTINUE
         JSTATE(NQN1*NSTATE+I)=II
  250 CONTINUE
C PRINTS ATOMIC PARAMETERS
      IF (L2.EQ.0) THEN
         WRITE(6,603)
  603    FORMAT(/"  ENERGY LEVELS CALCULATED FROM e- SPIN ZEEMAN")
      ELSE
         WRITE(6,604) A2
  604    FORMAT(/"  ENERGY LEVELS CALCULATED FROM (IN CM-1):"/
     &           "   A2 (SPIN-ORBIT)         = ",F15.5)
      ENDIF
      IF (HFFLG) WRITE(6,*) " AND HYPERFINE CONSTANTS PROVIDED"
C SETS UP ZERO OF ENERGY
      IF (J2ZERO.GE.0) THEN
        IF (J2ZERO.GE.J2MN.AND.J2ZERO.LE.J2MX) THEN
          EZERO=0.125D0*A2*(J2zero*(J2zero+2)-L2*(L2+2)-IS2*(IS2+2))
          WRITE(6,114) J2ZERO,EZERO
 114  FORMAT(/"  J2ZERO=",I2," SETS ZERO OF ENERGY EZERO=",F15.5)
        ELSE
          WRITE(6,*)"  INVALID VALUE OF J2ZERO"
          STOP
        ENDIF
      ELSEIF (EZERO.NE.0) THEN
        WRITE(6,115) EZERO
      ENDIF
 115  FORMAT(/"  FROM INPUT, ZERO OF ENERGY EZERO=",F15.5)
C SETS UP ELEVEL WITH DIAGONAL ELEMENTS OF HMON
      DO I=1,NLEVEL
         ITMP=I
         IF (II1.NE.0) THEN
c            MI1=JSTATE(ITMP)
            ITMP=ITMP+NSTATE
         ENDIF
         IF (IS1.NE.0) THEN
c            MS1=JSTATE(ITMP)
            ITMP=ITMP+NSTATE
         ENDIF
         IF (II2.NE.0) THEN
c            MI2=JSTATE(ITMP)
            ITMP=ITMP+NSTATE
         ENDIF
         IF (L2.NE.0) THEN
            JJ2=JSTATE(ITMP)
c            ITMP=ITMP+NSTATE
         ENDIF
c         MJ2=JSTATE(ITMP)
         ELEVEL(I)=0.125D0*A2*(JJ2*(JJ2+2)-L2*(L2+2)-IS2*(IS2+2))-ezero
      ENDDO
      RETURN
C
      ENTRY BASE9(LCOUNT,N,JTOT,IBLOCK,JSTATE,NSTATE,NQN,JSINDX,L,
     1            IPRINT)
C
C     BASE9 IS CALLED EITHER TO COUNT THE ACTUAL NUMBER OF CHANNEL BASIS
C     FUNCTIONS OR ACTUALLY TO SET THEM UP (IN THE JSINDX AND L ARRAYS).
C     IT IS CALLED FOR EACH TOTAL J (JTOT) AND SYMMETRY BLOCK (IBLOCK).
C     IF LCOUNT IS .TRUE. ON ENTRY, JUST COUNT THE BASIS FUNCTIONS. OTHERWISE,
C     SET UP JSINDX (POINTER TO JSTATE) AND L (ORBITAL ANGULAR MOMENTUM) FOR
C     EACH CHANNEL.  THIS MUST TAKE INTO ACCOUNT JTOT AND IBLOCK.
C
C     ONE IMPORTANT OVERALL EFFECT IS THAT THE THRESHOLD ENERGY IS IN
C     ELEVEL(JSTATE(NSTATE*(NQN-1)+JSINDX(I)). CHECK THIS!
C
C     THIS VERSION USES JTOT FOR |2 x MTOT| (PROJECTION OF TOTAL J)
C     AND IBLOCK = 1/2 FOR -/+ PARITY
C     SIGN OF MTOT NOT NEEDED BECAUSE FIELD(S) CAN BE SET TO NEGATIVE
C
C *** NOTE *** ALL MONOMER QUANTUM NUMBERS & MTOT ARE TWICE THEIR VALUES
C
      MTOT=JTOT
      N=0
      MI1=0
      MS1=0
      MI2=0
      DO 320 I=1,NSTATE
         ITMP=I
         IF (II1.NE.0) THEN
            MI1=JSTATE(ITMP)
            ITMP=ITMP+NSTATE
         ENDIF
         IF (IS1.NE.0) THEN
            MS1=JSTATE(ITMP)
            ITMP=ITMP+NSTATE
         ENDIF
         IF (II2.NE.0) THEN
            MI2=JSTATE(ITMP)
            ITMP=ITMP+NSTATE
         ENDIF
         IF (L2.NE.0) THEN
c            JJ2=JSTATE(ITMP)
            ITMP=ITMP+NSTATE
         ENDIF
         MJ2=JSTATE(ITMP)
         LMIN=ABS(MTOT-MI1-MS1-MI2-MJ2)
         IF (ODD(LMIN)) THEN
            STOP '*** ERROR-MTOT INCONSISTENT WITH SPIN QUANTUM NUMBERS'
         ELSE
            LMIN=LMIN/2
         ENDIF
         DO 310 LL=LMIN,LMAX
C PARITY P = p1p2(-1)^L: p1 = 1, p2 = (-1)^l2 (JPCA 108, 8941 (2004))
C IBLOCK 1/2 CORRESPONDS TO PARITY -/+
            IF (ODD(IBLOCK+L2/2+LL)) GOTO 310
            N=N+1
            IF (LCOUNT) GOTO 310
            JSINDX(N)=I
            L(N)=LL
 310     CONTINUE
 320  CONTINUE
      RETURN
C
      ENTRY CPL9(N,IBLOCK,NHAM,LAM,MXLAM,NSTATE,JSTATE,JSINDX,L,JTOT,
     1           VL,IV,CENT,DGVL,IBOUND,IEXCH,IPRINT)
C
C     CPL9 IS CALLED AFTER BASE9 FOR EACH JTOT AND IBLOCK, TO SET UP THE
C     POTENTIAL COUPLING COEFFICIENTS VL.
C     IF IBOUND>0, IT ALSO SETS UP THE CENTRIFUGAL COEFFICIENTS CENT.
C     INDICES SPECIFYING THE MXLAM DIFFERENT POTENTIAL SYMMETRIES ARE IN
C     THE FIRST XX*MXLAM ELEMENTS OF LAM; THE STRUCTURE OF THE LAM ARRAY
C     (AND THE VALUE OF XX) IS CHOSEN BY THE PROGRAMMER, AND MUST
C     CORRESPOND WITH THAT USED IN SUBROUTINE POTENL.
C     NHAM IS THE NUMBER OF DIFFERENT POTENTIAL TERMS WHICH CONTRIBUTE TO
C     EACH MATRIX ELEMENT (SEE SUBROUTINE WAVVEC). IT SOMETIMES SAVES
C     A SIGNIFICANT AMOUNT OF SPACE IF IT CAN BE LESS THAN MXLAM.
C
C     THIS VERSION USES JTOT FOR |2 x MTOT| (PROJECTION OF TOTAL J)
C     AND IBLOCK = 1/2 FOR -/+ PARITY
C
C *** NOTE *** ALL MONOMER QUANTUM NUMBERS & MTOT ARE TWICE THEIR VALUES
C
      MTOT=JTOT
      MI1C=0
      MI1R=0
      MS1C=0
      MS1R=0
      MI2C=0
      MI2R=0
C THE VALUES OF MAPEFV AND NCONST ARE FIXED IN BAS9IN
C     MAPEFV=2 ! IF ONLY SPIN-ORBIT & ZEEMAN TERMS ARE INCLUDED
C     MAPEFV=3 ! IF SPIN-ORBIT, HYPERFINE & ZEEMAN TERMS ARE INCLUDED
      IF (IS(2).NE.-1) M2S1=.TRUE.
      IMFLD=MXLAM+MAPEFV
      NHAM=MXLAM+NCONST ! + NRSQ == 0
      WRITE(6,510) NHAM
 510  FORMAT('  NHAM =',I5,' (INCLUDING NCONST TERMS)')
      MXLL=NHAM
      DO 550 LL=1,MXLL
         NNZ=0
         I=LL
C PICKS UP COLUMN QUANTUM NUMBERS
         DO 540 ICOL=1,N
            ITMP=JSINDX(ICOL)
            IF (II1.NE.0) THEN
               MI1C=JSTATE(ITMP)
               ITMP=ITMP+NSTATE
            ENDIF
            IF (IS1.NE.0) THEN
               MS1C=JSTATE(ITMP)
               ITMP=ITMP+NSTATE
            ENDIF
            IF (II2.NE.0) THEN
               MI2C=JSTATE(ITMP)
               ITMP=ITMP+NSTATE
            ENDIF
            IF (L2.NE.0) THEN
               J2C=JSTATE(ITMP)
               ITMP=ITMP+NSTATE
            ELSE
               J2C=IS2
            ENDIF
            MJ2C=JSTATE(ITMP)
            LC=L(ICOL)
            MLC=(MTOT-MI1C-MS1C-MI2C-MJ2C)/2

C PICKS UP ROW QUANTUM NUMBERS
            DO 530 IROW=1,ICOL
               ITMP=JSINDX(IROW)
               IF (II1.NE.0) THEN
                  MI1R=JSTATE(ITMP)
                  ITMP=ITMP+NSTATE
               ENDIF
               IF (IS1.NE.0) THEN
                  MS1R=JSTATE(ITMP)
                  ITMP=ITMP+NSTATE
               ENDIF
               IF (II2.NE.0) THEN
                  MI2R=JSTATE(ITMP)
                  ITMP=ITMP+NSTATE
               ENDIF
               IF (L2.NE.0) THEN
                  J2R=JSTATE(ITMP)
                  ITMP=ITMP+NSTATE
               ELSE
                  J2R=IS2
               ENDIF
               MJ2R=JSTATE(ITMP)
               LR=L(IROW)
               MLR=(MTOT-MI1R-MS1R-MI2R-MJ2R)/2


C SETS UP SOME 'DIAGONALITY SHORTCUTS'
               D1C1R=DIAG2(MI1C,MS1C, MI1R,MS1R)
               D2C2R=DIAG3(MI2C,J2C,MJ2C, MI2R,J2R,MJ2R)
               DHF12=DIAG2(MI1C,MI2C, MI1R,MI2R)
C COMPUTES VL MATRIX ELEMENT
               VL(I)=0.D0
C CONVENTION HERE IS COLUMN <=> "NO PRIME", ROW <=> "PRIME"
C INCLUDES CP1-CP2 INTERACTION TERMS: <a|H12|a'>
      IF (LL.LE.MXLAM) THEN
       LAMBDA=LAM(LL)
C * CP2 SPIN-ORBIT VARIATION WITH CP1 APPROACH: <a|Hso(R)-Hso(inf)|a'>
C *** NOTE *** THIS TERM VARIES WITH R, AND IS SET UP IN VSTAR
       IF (LAMBDA.EQ.-2) THEN
        IF (L2.EQ.0 .OR. IS2.EQ.0) GOTO 520
        IF (DIAG2(LC,MLC, LR,MLR) .AND. D1C1R .AND. D2C2R)
     &   VL(I)=0.125D0*(J2C*(J2C+2)-L2*(L2+2)-IS2*(IS2+2))
C * CP1 MAGNETIC DIPOLE-CP2 MAGNETIC DIPOLE: <a|Vdd|a'>
C   (CALCULATED FROM "DECOUPLED" MATRIX ELEMENTS)
C *** NOTE *** THIS TERM VARIES AS 1/R^3, WHICH IS SET UP IN VSTAR
       ELSEIF (LAMBDA.EQ.-1 .AND. IS1.NE.0) THEN
        IF (.NOT.DHF12) GOTO 520
        IXS1L(1)=IS1
        IXS1L(2)=MS1C
        IXS1L(3)=MS1R
        IXS1L(4)=LC
        IXS1L(5)=MLC
        IXS1L(6)=LR
        IXS1L(7)=MLR
        IF (L2.NE.0) VL(I)
     &   =FS1L2*cHc_from_dHd(0.5D0*L2,0.5D0*IS2,0.5D0*J2C,0.5D0*MJ2C,
     &                       0.5D0*L2,0.5D0*IS2,0.5D0*J2R,0.5D0*MJ2R,
     &                       dS1L2d,IXS1L)
        IF (IS2.NE.0) VL(I)=VL(I)
     &   +FS1S2*cHc_from_dHd(0.5D0*L2,0.5D0*IS2,0.5D0*J2C,0.5D0*MJ2C,
     &                       0.5D0*L2,0.5D0*IS2,0.5D0*J2R,0.5D0*MJ2R,
     &                       dS1S2d,IXS1L)
C * INTERACTION POTENTIAL: <a|V|a'>
       ELSEIF (LAMBDA.GE.0) THEN
        IF (.NOT.DHF12) GOTO 520
        IF (.NOT.CALCPOT) GOTO 520 !DGG ADDED
        MK=MLR-MLC

CC . FOR SEVERAL SPIN MULTIPLICITIES, LAMBDA IS MAPPED AS:
CC   LAMBDA = (I2S1-1)*10 + K, WITH
CC   (1) I2S1 => ~ (I2S1+1)TH SPIN MULTIPLICITY (MAX. 10 DIFFERENT I2S1)
CC   (2) K => 0, 2, ..., 2*L2 TERM IN THE POTENTIAL EXPANSION, GIVEN I2S1
C        I2S1=INT(LAMBDA/10)+1
C        S=0.5D0*IS(I2S1)
C        K=LAMBDA-10*(I2S1-1)

CCC MDF 19/08/19: MAPPING CHANGED TO COPE WITH ATOMS WITH LARGE VALUES
CCC OF L, SUCH AS LANTHANIDES.
C . FOR SEVERAL SPIN MULTIPLICITIES, LAMBDA IS MAPPED AS:
C   LAMBDA = (I2S1-1)*10 + K/2, WITH
C   (1) I2S1 => ~ (I2S1+1)TH SPIN MULTIPLICITY (MAX. 10 DIFFERENT I2S1)
C   (2) K => 0, 2, ..., 2*L2 TERM IN THE POTENTIAL EXPANSION, GIVEN I2S1
        I2S1=INT(LAMBDA/10)+1
        S=0.5D0*IS(I2S1)
        K=(LAMBDA-10*(I2S1-1))*2

C     SWITCHES FOR V_K^S: CALCULATE VL ARRAY OR SET TO ZERO FOR EACH K=0,2; S=3/2,1/2
C        IF (LL.GE.0) THEN        !AND .LE.MXLAM (THIS CONDITION IS ABOVE)
C           IF (LAMBDA.EQ.0 .AND.
C     $          .NOT.CALC_VK0S1) GOTO 520
C           IF (LAMBDA.EQ.2 .AND.
C     $          .NOT.CALC_VK2S1) GOTO 520
C           IF (LAMBDA.EQ.10 .AND.
C     $          .NOT.CALC_VK0S2) GOTO 520
C           IF (LAMBDA.EQ.12 .AND.
C     $          .NOT.CALC_VK2S2) GOTO 520
C       ENDIF

C .. COMPUTES COMMON PRE-FACTOR
C ... EXPLOITS 3-J SYMBOLS PROPERTIES
        IF (ABS(MK).GT.K .OR. ODD(LC+K+LR)) GOTO 520
        FCT1=(L2+1)*SQRT(DBLE((2*LC+1)*(2*LR+1)))
     &       *THREEJ(L2/2,K,L2/2)*THREEJ(LC,K,LR)
     &       *THRJ(DBLE(LC),DBLE(K),DBLE(LR),
     &            -DBLE(MLC),-DBLE(MK),DBLE(MLR))
        IF (M2S1) THEN
C .. INCLUDES MATRIX ELEMENT FOR MULTIPLE SPIN MULTIPLICITIES
C    (CALCULATED FROM "DECOUPLED" MATRIX ELEMENTS)
         IXSK(1)=IS(I2S1)
         IXSK(2)=IS1
         IXSK(3)=MS1C
         IXSK(4)=MS1R
         IXSK(5)=K
         IXSK(6)=MK
         VL(I)=PARSGN(MK-MLC)*FCT1
     &         *cHc_from_dHd(0.5D0*L2,0.5D0*IS2,0.5D0*J2C,0.5D0*MJ2C,
     &                       0.5D0*L2,0.5D0*IS2,0.5D0*J2R,0.5D0*MJ2R,
     &                       dM2S1d,IXSK)

        ELSE !if only single spin multiplicity use the next routine instead
         IF (.NOT.D1C1R) GOTO 520
C .. EQ.(9) IN PHYS. REV. A 68, 013406 (2007)
C ... EXPLOITS 3-J SYMBOLS PROPERTIES
         IF (2*MK.NE.(MJ2C-MJ2R)) GOTO 520
         FCT2=SQRT(DBLE((J2C+1)*(J2R+1)))
         FCT3=THRJ(0.5D0*J2C,DBLE(K),0.5D0*J2R,
     &           -0.5D0*MJ2C,DBLE(MK),0.5D0*MJ2R)
         FCT4=DSIXJ(XL2,0.5D0*J2C,XL2,0.5D0*J2R,S2,DBLE(K))
         VL(I)=PARSGN((IS2+J2C+J2R-MJ2C)/2+K+MK-MLC)*FCT1*FCT2*FCT3*FCT4
        ENDIF ! end if M2S1 condition
       ENDIF ! if lambda .gt. 0

*         write(*,*) LAMBDA, k
!       write(46,*) k,S, VL(I) !DGG bugchk
      ELSE ! if LL .gt. mxlam


C INCLUDES MONOMER TERMS: <a|H1|a'> AND <a|H2|a'>
       IF (DIAG2(LC,MLC, LR,MLR)) THEN
        IF (LL.EQ.MXLAM+1) THEN
C .  CP2
C .. SPIN-ORBIT: <a|Hso|a'>
         IF (D1C1R .AND. L2.NE.0 .AND. IS2.NE.0) THEN
            IF (CALCSO2) THEN !DGG BUGCHK FLAG TO 'TURN ON/OFF' SO ELEMENTS
          IF (D2C2R)
     &     VL(I)=0.125D0*(J2C*(J2C+2)-L2*(L2+2)-IS2*(IS2+2))-EZERO/A2
!     &     VL(I)=0.125D0*(J2C*(J2C+2)-L2*(L2+2)-IS2*(IS2+2))+1.D0
       ENDIF !DGG BUGCHK
         ENDIF
C * HYPERFINE (ONLY IF HFFLG = .TRUE.): <a|Hhf|a'>
C THE kHz TO CM-1 CONVERSION FACTOR IS INCLUDED VIA WAVMAT
        ELSEIF (LL.EQ.MXLAM+2 .AND. HFFLG) THEN
C . CP1
         IF (D2C2R .AND. II1.NE.0) THEN
C .. FERMI(-BREIT)-CONTACT INTERACTION
          IF (IS1.NE.0 .AND. CALCHF1) !DGG BUGCHK CALCHF1
     &     VL(I)=BF1*T1J1oT1J2_D(IS1,MS1C,MS1R,II1,MI1C,MI1R)
         ENDIF
C . CP2
         IF (D1C1R .AND. II2.NE.0) THEN
          IF (L2.EQ.0) THEN
           IF (IS2.NE.0)
     &      VL(I)=VL(I)+BF2*T1J1oT1J2_D(J2C,MJ2C,MJ2R,II2,MI2C,MI2R)
          ELSE
           IF (J2C.EQ.J2R)
     &      VL(I)=VL(I)
     &            +AJ(J2C,J2C)*T1J1oT1J2_D(J2C,MJ2C,MJ2R,II2,MI2C,MI2R)
c           IF (DIAG3(MI2C,J2C,MJ2C, MI2R,J2R+2,MJ2R))
c     &      VL(I)=VL(I)
c     &            +0.25D0*AJ(J2C,J2R)*MI2C*SQRT(DBLE(J2C*J2C-MJ2C*MJ2C))
          ENDIF
         ENDIF
C * ZEEMAN: <a|HZ|a'>
      ELSEIF (LL.EQ.IMFLD) THEN
C . CP1
         IF (D1C1R) THEN
            IF (D2C2R) THEN
               IF (CALCZEEMAN_I1 .AND. II1.NE.0) VL(I)=FZI1*MI1C

               IF (CALCZEEMAN_S1 .AND. IS1.NE.0) VL(I)=VL(I)+FZS1*MS1C
            ENDIF
C . CP2
          IF (II2.NE.0 .AND. D2C2R) VL(I)=VL(I)+FZI2*MI2C
          IF (L2.EQ.0 .AND. D2C2R .AND. CALCZEEMAN_S2) !DGG BUGCHK
     $         VL(I)=VL(I)+0.5D0*FZS2*MJ2C


             IF (L2.NE.0 .AND. DIAG2(MI2C,MJ2C, MI2R,MJ2R)) THEN
                FCT=SQRT(DBLE((J2C+1)*(J2R+1)))
     &               *THRJ(0.5D0*J2C,1.D0,0.5D0*J2R,
     &               -0.5D0*MJ2C,0.D0,0.5D0*MJ2C)

                IF (IS2.NE.0 .AND. CALCZEEMAN_S2) THEN
                   VL(I)=VL(I)+FCT*FZS2*PARSGN(J2C+(L2+IS2-MJ2C)/2+1)
     &                  *SQRT(S2*(S2+1)*(IS2+1))
     &                  *DSIXJ(S2,0.5D0*J2R,S2,0.5D0*J2C,XL2,1.D0)

 321               FORMAT(1x, E12.5, 1x, E12.5)
                ENDIF

                IF (L2.NE.0 .AND. CALCZEEMAN_L2) THEN
                  VL(I)=VL(I)+FCT*FZL2*PARSGN((J2C+J2R+L2+IS2-MJ2C)/2+1)
     &                  *SQRT(XL2*(XL2+1)*(L2+1))
     &                  *DSIXJ(XL2,0.5D0*J2R,XL2,0.5D0*J2C,S2,1.D0)
                ENDIF

             ENDIF
          ENDIF
       ENDIF                   ! LL > MXLAM == MONOMER TERMS
       ENDIF ! DIAG(L, ML)
      ENDIF ! LL
      IF (VL(I).NE.0.D0) NNZ=NNZ+1

  520 I=I+NHAM
      mcountindex=I-NHAM

! === DGG DIAGNOSTICS:
       if (WRITEVLMATRICES .AND. LL.LE.MXLAM) THEN
          IF (LAMBDA.EQ.0) then
             write(66,*) IROW, ICOL, VL(mcountindex)
          ELSEIF (LAMBDA.EQ.2) then
             write(67,*) IROW, ICOL, VL(mcountindex)
          ELSEIF (LAMBDA.EQ.10) then
!             write(*,*) "lambda.eq.10; mcountindex=",mcountindex , !DGG bugchk
!     $            VL(mcountindex)
             write(68,*) IROW, ICOL, VL(mcountindex)
          ELSEIF (LAMBDA.EQ.12) then
             write(69,*) IROW, ICOL, VL(mcountindex)
          ENDIF
       ENDIF

  530 CONTINUE
  540 CONTINUE
      IF (NNZ.EQ.0) WRITE(6,612) JTOT,LL
  612 FORMAT('  * * * NOTE.  FOR JTOT =',I4,',  ALL COUPLING',
     1  ' COEFFICIENTS ARE 0.0 FOR POTENTIAL SYMMETRY',I4)
  550 CONTINUE
      RETURN
C
      ENTRY THRSH9(JREF,MONQN,NQN2,EREF,IPRINT)
C
C  THIS CALCULATES THRESHOLDS FOR A TRIPLET-SIGMA MOLECULE INTERACTING
C  WITH A STRUCTURELESS ATOM (WHICH CONTRIBUTES NOTHING).
C  THE MONOMER QUANUM NUMBERS ARE SPECIFIED IN THE ARRAY MONQN, AND ARE:
C     MONQN(1): f
C     MONQN(2): m_fj
C     MONQN(3): j
C     MONQN(4): m_j
C
      BFIELD=EFV(1)
      IF (JREF.GT.0) THEN
        WRITE(6,*) ' *** ERROR - THRSH9 CALLED WITH POSITIVE IREF'
        STOP
      ENDIF
C
      IF (MONQN(1).EQ.-99999) THEN
        WRITE(6,*) ' *** ERROR - THRSH9 CALLED WITH MONQN UNSET'
        STOP
      ENDIF
C
      IFF1=MONQN(1)
      MF1=MONQN(2)
      JJ2=MONQN(3)
      MJ2=MONQN(4)
C
C  BREIT-RABI FOR ATOM A
C
      IF (IS1.EQ.0.AND.II1.EQ.0) THEN
        E1=0.D0
      ELSEif (is1.eq.1) then
        IF (ABS(IFF1-II1).NE.1) THEN
          WRITE(6,*) ' *** THRSH9: INVALID MONQN(1) =',IFF1
          STOP
        ELSEIF (ABS(MF1).GT.IFF1) THEN
          WRITE(6,*) ' *** THRSH9: MF1 =',MF1,' > F1. STOPPING'
          STOP
        ELSEIF (MOD(MF1+IFF1,2).NE.0) THEN
          WRITE(6,*) ' *** THRSH9: INVALID MONQN(1),MONQN(2) PAIR =',
     1                IFF1,MF1
          STOP
        ENDIF
        HFSPLA=BF1*0.5*(II1+1.0)/CMIKHZ
C
        EA1=-HFSPLA/(2.D0*DBLE(II1+1)) - 0.5D0*G1*XNM*DBLE(MF1)*BFIELD
        BX=BM*BFIELD*(GS1+G1*(XNM/BM))/HFSPLA
        EA2=0.5D0*HFSPLA*SQRT(1.D0+DBLE(MF1+MF1)*BX/DBLE(II1+1)+BX*BX)
C
        IF (ABS(MF1).EQ.II1+1) THEN
          E1=DBLE(II1)*HFSPLA/(2D0*DBLE(II1+1))+SIGN(1.D0,DBLE(MF1))
     1       *BM*BFIELD*(GS1*IS1*0.5D0-G1*(XNM/BM)*II1*0.5D0)
        ELSEIF (IFF1.EQ.II1+1) THEN
          E1=EA1+EA2
        ELSEIF (IFF1.EQ.II1-1) THEN
          E1=EA1-EA2
        ENDIF
      ELSE
        WRITE(6,*) ' *** THRSH9: NOT IMPLEMENTED FOR IS1.GT.1'
      ENDIF
C
      IF (IS2.EQ.0.OR.L2.EQ.0.OR.II2.NE.0) THEN
        WRITE(6,*) ' *** ERROR - THRSH9 NOT IMPLEMENTED FOR THESE SPINS'
        STOP
      ENDIF
C THIS SIMPLY CALCULATES THE LINEAR ZEEMAN EFFECT FOR THE LOW FIELD LIMIT
C FROM THE G FACTOR, AND IGNORES HIGHER ORDER TERMS.
      E2=A2*(0.125D0*(JJ2*(JJ2+2)-L2*(L2+2)-IS2*(IS2+2)))-EZERO
      GJ2=GL2*(JJ2*(JJ2+2)-IS2*(IS2+2)+L2*(L2+2))/(2*JJ2*(JJ2+2))
      GJ2=GJ2+GS2*(JJ2*(JJ2+2)+IS2*(IS2+2)-L2*(L2+2))/(2*JJ2*(JJ2+2))
      E2=E2+0.5D0*MJ2*GJ2*BM*BFIELD
      EREF=E1+E2
      IF (IPRINT.GE.8) THEN
        WRITE(6,*)
        WRITE(6,667) '1',MONQN(1),MONQN(2),E1
        WRITE(6,667) '2',MONQN(3),MONQN(4),E2
        WRITE(6,*) '  WARNING: OFF-DIAGONAL (IE NON-LINEAR) ZEEMAN ',
     &              'TERMS ON ATOM 2 HAVE BEEN NEGLECTED, ENERGIES ',
     &              'MAY BE INACCURATE AT HIGHER FIELDS'
  667   FORMAT('  ATOM ',A1,' WITH DOUBLED QUANTUM NOS',2I3,
     1     ' IS AT ENERGY',F19.12,' CM-1')
      ENDIF

      END
C----------------------------- T1J1oT1J2_D -----------------------------
      FUNCTION T1J1oT1J2_D(J1,M11,M12,J2,M21,M22)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-----------------------------------------------------------------------
C  DATE (LAST UPDATE): 12/07/10        STATUS: FINISHED                |
C  AUTHOR: MAYKEL LEONARDO GONZALEZ MARTINEZ                           |
C         (ADAPTED FROM JMH'S SDOTI2 FUNCTION IN base9.alk.f)          |
C-----------------------------------------------------------------------
C    FUNCTION FOR A MATRIX ELEMENT OF A DOT PRODUCT OF TWO SPHERICAL   |
C       TENSORS OF RANK 1, ACTING DIRECTLY ON TWO ANGULAR MOMENTA      |
C               USING A  TOTALLY DECOUPLED BASIS SET:                  |
C           <j1 m1|<j2 m2|T^1(j1).T^1(j2)|j2 m'2>|j1 m'1>              |
C  (ALL INPUT INTEGERS ARE TWICE THE CORRESPONDING QUANTUM NUMBERS)    |
C-----------------------------------------------------------------------
      T1J1oT1J2_D=0.D0
      IF (M11+M21.EQ.M12+M22) THEN
         IF (M11.EQ.M12) THEN
            T1J1oT1J2_D=0.25D0*(M11*M21)
         ELSEIF (ABS(M11-M12).EQ.2) THEN
            T1J1oT1J2_D=0.125D0*SQRT(DBLE(J1*(J1+2)-M11*M12))
     &       *SQRT(DBLE(J2*(J2+2)-M21*M22))
         ENDIF
      ENDIF
C
      RETURN
      END FUNCTION T1J1oT1J2_D
C----------------------------- T2J1J2oT2C_D ----------------------------
      FUNCTION T2J1J2oT2C_D(ISA,MSAC,MSAR,ISB,MSBC,MSBR,LC,MLC,LR,MLR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-----------------------------------------------------------------------
C  DATE (LAST UPDATE): 10/08/10        STATUS: FINISHED                |
C  AUTHOR: MAYKEL LEONARDO GONZALEZ MARTINEZ                           |
C         (ADAPTED FROM JMH'S SPINSP FUNCTION IN base9.alk2.f)         |
C-----------------------------------------------------------------------
C     FUNCTION FOR A MATRIX ELEMENT OF A PRODUCT OF TWO SPHERICAL      |
C     TENSORS OF RANK 2 (OF THE KIND APPEARING IN SPIN-SPIN TERMS,     |
C    E.G. BROWN & CARRINGTON'S EQ. 8.464), USING A DECOUPLED BASIS     |
C    <L ML|<j1 m1|<j2 m2|T^2(j1,j2).T^2(C)|j2 m'2>|j1 m'1>|L' M'L>     |
C   (ALL INPUT INTEGERS EXCEPT Lx & MLx ARE TWICE THE CORRESPONDING    |
C                            QUANTUM NUMBERS)                          |
C-----------------------------------------------------------------------
      T2J1J2oT2C_D=0.D0
      IF (MSAC+MSBC+2*MLC.NE.MSAR+MSBR+2*MLR) RETURN
      IF (ABS(LC-LR).GT.2 .OR. ABS(MLC-MLR).GT.2) RETURN
C
      XSA=0.5D0*ISA
      XMSAC=0.5D0*MSAC
      XMSAR=0.5D0*MSAR
      XSB=0.5D0*ISB
      XMSBC=0.5D0*MSBC
      XMSBR=0.5D0*MSBR
      XMLC=DBLE(MLC)
      XMLR=DBLE(MLR)
      IQ=MLR-MLC
      FAC1=PARSGN((ISA-MSAC+ISB-MSBC)/2-MLC)
     & *SQRT(XSA*(XSA+1.D0)*(XSA+XSA+1.D0)*XSB*(XSB+1.D0)*(XSB+XSB+1.D0)
     & *(2*LC+1)*(2*LR+1))*THREEJ(LC,2,LR)
      FAC2=THRJ(DBLE(LC),2.D0,DBLE(LR), -XMLC,-DBLE(IQ),XMLR)
      LIQA: DO IQA=-1,1
        IQB=IQ-IQA
        IF (ABS(IQB).GT.1) CYCLE LIQA
        FAC3=THRJ(1.D0,1.D0,2.D0, DBLE(IQA),DBLE(IQB),-DBLE(IQ))
        FAC4=THRJ(XSA,1.D0,XSA, -XMSAC,DBLE(IQA),XMSAR)
        FAC5=THRJ(XSB,1.D0,XSB, -XMSBC,DBLE(IQB),XMSBR)
        T2J1J2oT2C_D=T2J1J2oT2C_D+FAC1*FAC2*FAC3*FAC4*FAC5
      ENDDO LIQA
      T2J1J2oT2C_D=SQRT(30.D0)*T2J1J2oT2C_D
C
      RETURN
      END FUNCTION T2J1J2oT2C_D
C ------------------------------ dS1L2d --------------------------------
      FUNCTION dS1L2d(XJ1C,XM1C,XJ2C,XM2C,XJ1R,XM1R,XJ2R,XM2R,IXS1L)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-----------------------------------------------------------------------
C  DATE (LAST UPDATE): 11/04/13        STATUS: FINISHED                |
C  AUTHOR: MAYKEL LEONARDO GONZALEZ MARTINEZ                           |
C-----------------------------------------------------------------------
      DIMENSION IXS1L(7)
C
      dS1L2d=0.D0
      IF (XJ1C.NE.XJ1R .OR. XJ2C.NE.XJ2R .OR. XM2C.NE.XM2R) RETURN
      L2  =2*XJ1C
      ML2C=2*XM1C
      ML2R=2*XM1R
      IS1 =IXS1L(1)
      MS1C=IXS1L(2)
      MS1R=IXS1L(3)
      LC  =IXS1L(4)
      MLC =IXS1L(5)
      LR  =IXS1L(6)
      MLR =IXS1L(7)
      dS1L2d=-T2J1J2oT2C_D(IS1,MS1C,MS1R,L2,ML2C,ML2R,LC,MLC,LR,MLR)
C
      RETURN
      END FUNCTION dS1L2d
C ------------------------------ dS1S2d --------------------------------
      FUNCTION dS1S2d(XJ1C,XM1C,XJ2C,XM2C,XJ1R,XM1R,XJ2R,XM2R,IXS1L)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-----------------------------------------------------------------------
C  DATE (LAST UPDATE): 11/04/13        STATUS: FINISHED                |
C  AUTHOR: MAYKEL LEONARDO GONZALEZ MARTINEZ                           |
C-----------------------------------------------------------------------
      DIMENSION IXS1L(7)
C
      dS1S2d=0.D0
      IF (XJ1C.NE.XJ1R .OR. XM1C.NE.XM1R .OR. XJ2C.NE.XJ2R) RETURN
      IS2 =2*XJ2C
      MS2C=2*XM2C
      MS2R=2*XM2R
      IS1 =IXS1L(1)
      MS1C=IXS1L(2)
      MS1R=IXS1L(3)
      LC  =IXS1L(4)
      MLC =IXS1L(5)
      LR  =IXS1L(6)
      MLR =IXS1L(7)
      dS1S2d=-T2J1J2oT2C_D(IS1,MS1C,MS1R,IS2,MS2C,MS2R,LC,MLC,LR,MLR)
C
      RETURN
      END FUNCTION dS1S2d
C ------------------------------ dM2S1d --------------------------------
      FUNCTION dM2S1d(XJ1C,XM1C,XJ2C,XM2C,XJ1R,XM1R,XJ2R,XM2R,IXSK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-----------------------------------------------------------------------
C  DATE (LAST UPDATE): 11/06/13        STATUS: FINISHED                |
C  AUTHOR: MAYKEL LEONARDO GONZALEZ MARTINEZ                           |
C-----------------------------------------------------------------------
      DIMENSION IXSK(6)
C
      dM2S1d=0.D0
      IF (XJ1C.NE.XJ1R .OR. XJ2C.NE.XJ2R) RETURN
      S   =0.5D0*IXSK(1)
      S1  =0.5D0*IXSK(2)
      IS1 =IXSK(2)
      MS1C=IXSK(3)
      MS1R=IXSK(4)
      K   =IXSK(5)
      MK  =IXSK(6)
      XL2 =XJ1C
      ML2C=2*XM1C
      ML2R=2*XM1R
      S2  =XJ2C
      IS2 =2*XJ2C
      MS2C=2*XM2C
      MS2R=2*XM2R
      MSC=MS1C+MS2C
      MSR=MS1R+MS2R
      IF (MSC.NE.MSR) RETURN
      FCT1=PARSGN(IS1+IS2+MSC-ML2C/2)*(2*S+1)
      FCT2=THRJ(S1,S2,S,
     &          0.5D0*MS1C,0.5D0*MS2C,-0.5D0*MSC)
      FCT3=THRJ(S1,S2,S,
     &          0.5D0*MS1R,0.5D0*MS2R,-0.5D0*MSC)
      FCT4=THRJ(XL2,DBLE(K),XL2,
     &          -0.5D0*ML2C,DBLE(MK),0.5D0*ML2R)
      dM2S1d=FCT1*FCT2*FCT3*FCT4
C
      RETURN
      END FUNCTION dM2S1d
