DataMuseum.dk

Presents historical artifacts from the history of:

RC4000/8000/9000

This is an automatic "excavation" of a thematic subset of
artifacts from Datamuseum.dk's BitArchive.

See our Wiki for more about RC4000/8000/9000

Excavated with: AutoArchaeologist - Free & Open Source Software.


top - metrics - download

⟦d4fa6f826⟧ TextFile

    Length: 32256 (0x7e00)
    Types: TextFile
    Names: »cpc7«

Derivation

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦this⟧ »cpc7« 

TextFile

c   cpc7
c
c   acrz gf values. oscillator strengths from nomerical mchf radial
c   funktions, with an adaption to calculate electric quadrupole
c   oscillator strengths.  c.f.fischer, k.m.s.saxena, m.godefroid.
c   ref. in comp. phys. commun. 9 (1975) 381 and
c           comp. phys. commun. 15 (1978) 275 and
c           comp. phys. commun. 17 (1979) 426.
c
C
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  *                                                                 *
C  *   JULY 1974                                                     *
C  *                                                                 *
C  *   PROGRAM TO EVALUATE LENGTH AND VELOCITY FORM OF 'GF' VALUES   *
C  *   USING NUMERICAL MCHF RADIAL FUNCTIONS.                        *
C  *                                                                 *
C  *                                                                 *
C  *   BY                                                            *
C  *                                                                 *
C  *                                                                 *
C  *   C. FROESE FISCHER AND K.M.S. SAXENA                           *
C  *                                                                 *
C  *   DEPARTMENT OF APPLIED MATHEMATICS                             *
C  *   UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO  N2L 3G1            *
C  *   C A N A D A.                                                  *
C  *                                                                 *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  *                                                                 *
C  *  DESCRIPTION OF THE COMMON BLOCKS:-                             *
C  *  ================================                               *
C  *                                                                 *
C  *  (1) COMMON BLOCK /PARAM/.                                      *
C  *      --------------------                                       *
C  *                                                                 *
C  *            THIS COMMON BLOCK STORES SEVERAL NUMERICAL CONSTANTS *
C  *  IN DOUBLE PRECISION AND VALUES OF CERTAIN PARAMETERS OF THE    *
C  *  CALCULATIONS.                                                  *
C  *                                                                 *
C  *            THE VARIABLES IN THIS COMMON BLOCK ARE AS FOLLOWS.   *
C  *                                                                 *
C  *     D0 : DOUBLE PRECISION CONSTANT 0                            *
C  *     D1 : DOUBLE PRECISION CONSTANT 1                            *
C  *     D2 : DOUBLE PRECISION CONSTANT 2                            *
C  *     D3 : DOUBLE PRECISION CONSTANT 3                            *
C  *     D4 : DOUBLE PRECISION CONSTANT 4                            *
C  *     D5 : DOUBLE PRECISION CONSTANT 1/2                          *
C  *     D6 : DOUBLE PRECISION CONSTANT 6                            *
C  *     D10: DOUBLE PRECISION CONSTANT 10                           *
c  *     d15: double precision constant 15                           *
c  *     d30: double precision constant 30                           *
C  *     H  : 0.0625,  THE STEP-SIZE DELTA(RHO)                      *
C  *     H1 : H/3                                                    *
C  *     Z  : ATOMIC NUMBER                                          *
C  *     NO : MAXIMUM NUMBER OF POINTS IN THE RANGE OF THE FUNCTIONS *
C  *          =220                                                   *
C  *     ND : NO-2                                                   *
C  *                                                                 *
C  *            MOST OF THESE CONSTANTS ARE SET IN THE BLOCK DATA    *
C  *  SEGMENT.                                                       *
C  *                                                                 *
C  *                                                                 *
C  *  (2) THE BLANK COMMON //.                                       *
C  *      -------------------                                        *
C  *                                                                 *
C  *            THIS COMMON SECTION STORES THE INFORMATION ABOUT     *
C  *  THE NUMERICAL MCHF WAVEFUNCTIONS.                              *
C  *                                                                 *
C  *            THE VARIABLES IN THIS COMMON BLOCK ARE AS FOLLOWS.   *
C  *                                                                 *
C  *     R(220)    : VALUES OF R(J)=EXP(RHO(J))/Z                    *
C  *     RR(220)   : VALUES OF R(J)**2                               *
C  *     R2(220)   : VALUES OF DSQRT(R(J))                           *
C  *     P(50,220) : VALUES OF P(R(J))/R2(J) FOR ORBITALS I=1,50     *
C  *     AZ(50)    : A0=LIMIT OF ( P(R)/(R**(L+1)) ) AS R->0 FOR     *
C  *                 ORBITALS I=1,50                                 *
C  *     L(50)     : L-QUANTUM NUMBERS FOR ORBITALS I=1,50           *
C  *     MAX(50)   : NUMBER OF MESH-POINTS IN THE NUMERICAL          *
C  *                 ORBITALS I=1,50                                 *
C  *                                                                 *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  *                                                                 *
C  *  DESCRIPTION OF THE OTHER IMPORTANT VARIABLES:-                 *
C  *  ============================================                   *
C  *                                                                 *
C  *  ASTER     :  AN INTEGER VARIABLE READ IN FORMAT(A1).           *
C  *  BLANK     :  A BLANK CHARACTER ' ', DEFINED AS AN INTEGER      *
C  *               AND SET BY A DATA STATEMENT IN THE MAIN PROGRAM.  *
C  *  SYMM(5)   :  AN INTEGER VARIABLE STORING THE SPECTROSCOPIC     *
C  *               SYMBOLS 'S','P','D','F', AND 'G' IN FORMAT(A1).   *
C  *               SYMM(I) CORRESPONDS TO THE SYMMETRY WITH L=(I-1). *
C  *  JREAD     :  GENERAL INPUT CHANNEL FOR A SET OF CALCULATIONS.  *
C  *  IREAD1    :  INPUT CHANNEL FOR THE WAVEFUNCTIONS REQUIRED IN   *
C  *               A CASE.                                           *
C  *  IREAD2    :  INPUT CHANNEL FOR THE TRANSITION INTEGRALS DATA   *
C  *               REQUIRED IN A CASE.                               *
C  *  IWRITE    :  OUTPUT CHANNEL REQUIRED TO OUTPUT THE RESULTS OF  *
C  *               A CASE.                                           *
C  *  ATOM      :  ATOMIC SYMBOL/DESIGNATION IN FORMAT(A8).          *
C  *  IZ        :  ATOMIC NUMBER Z.                                  *
C  *  NWF       :  NUMBER OF ORBITALS INVOLVED.                      *
C  *  TRMI      :  TERM SYMBOL FOR THE INITIAL STATE IN FORMAT(A8).  *
C  *  TRMF      :  TERM SYMBOL FOR THE  FINAL  STATE IN FORMAT(A8).  *
C  *  NCFGI     :  NUMBER OF CONFIGURATIONS IN THE INITIAL STATE.    *
C  *  NCFGF     :  NUMBER OF CONFIGURATIONS IN THE  FINAL  STATE.    *
C  *  MULT      :  MULTIPLICITY (2S+1) OF THE TWO STATES.            *
c  *  kctas     :  phase convention                                  *
c  *                  = 0 fano and racah                             *
c  *                  = 1 condon and shortley                        *
C  *  EL(50,3)  :  AN INTEGER VARIABLE STORING THE THREE CHARACTER   *
C  *               ALPHANUMERIC SYMBOL OF THE ORBITALS I=1,50 IN     *
C  *               FORMAT(3A1) FOR EACH OF THEM.                     *
C  *  ETI       :  TOTAL ENERGY (IN A.U.) FOR THE INITIAL STATE.     *
C  *  ETF       :  TOTAL ENERGY (IN A.U.) FOR THE  FINAL  STATE.     *
C  *  CONFAI(30):  CHARACTERS 01-08 (IN A8 FORMAT) OF THE ALPHA-     *
C  *               NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF     *
C  *               THE INITIAL STATE CONFIGURATIONS I=1,30.          *
C  *  CONFBI(30):  CHARACTERS 09-16 (IN A8 FORMAT) OF THE ALPHA-     *
C  *               NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF     *
C  *               THE INITIAL STATE CONFIGURATIONS I=1,30.          *
C  *  CONFCI(30):  CHARACTERS 17-24 (IN A8 FORMAT) OF THE ALPHA-     *
C  *               NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF     *
C  *               THE INITIAL STATE CONFIGURATIONS I=1,30.          *
C  *  CONFAF(30):  CHARACTERS 01-08 (IN A8 FORMAT) OF THE ALPHA-     *
C  *               NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF     *
C  *               THE  FINAL  STATE CONFIGURATIONS I=1,30.          *
C  *  CONFBF(30):  CHARACTERS 09-16 (IN A8 FORMAT) OF THE ALPHA-     *
C  *               NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF     *
C  *               THE  FINAL  STATE CONFIGURATIONS I=1,30.          *
C  *  CONFCF(30):  CHARACTERS 17-24 (IN A8 FORMAT) OF THE ALPHA-     *
C  *               NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF     *
C  *               THE  FINAL  STATE CONFIGURATIONS I=1,30.          *
C  *  WTI(30)   :  WEIGHTS OF THE INITIAL STATE CONFIGURATIONS I=1,30*
C  *  WTF(30)   :  WEIGHTS OF THE  FINAL  STATE CONFIGURATIONS I=1,30*
C  *  COEF      :  MULTILPYING COEFFICIENT IN A TRANSITION INTEGRAL. *
c  *  ktens     :  the order of the tensor                           *
c  *                   = 1 dipole transition                         *
c  *                   = 2 quadropol transition                      *
C  *  KRHO      :  THE NUMBER INDICATING THE ORBITAL EL(KRHO) FROM   *
C  *               THE INITIAL STATE CONFIGURATION APPEARING IN A    *
C  *               TRANSITION INTEGRAL.                              *
C  *  KSIG      :  THE NUMBER INDICATING THE ORBITAL EL(KSIG) FROM   *
C  *               THE  FINAL  STATE CONFIGURATION APPEARING IN A    *
C  *               TRANSITION INTEGRAL.                              *
C  *  JI,JF     :  THE NUMBERS INDICATING THE INITIAL AND FINAL      *
C  *               STATE CONFIGURATIONS FOR A TRANSITION INTEGRAL.   *
C  *  II1,II2,  :  THE NUMBERS INDICATING THE ORBITALS EL(II1),      *
C  *  II3          EL(II2), AND EL(II3) FROM THE INITIAL STATE       *
C  *               CONFIGURATION JI APPEARING IN THE THREE OVERLAP   *
C  *               INTEGRALS IN A TRANSITION INTEGRAL CONTRIBUTION.  *
C  *  IF1,IF2,  :  THE NUMBERS INDICATING THE ORBITALS EL(IF1),      *
C  *  IF3          EL(IF2), AND EL(IF3) FROM THE  FINAL  STATE       *
C  *               CONFIGURATION JF APPEARING IN THE THREE OVERLAP   *
C  *               INTEGRALS IN A TRANSITION INTEGRAL CONTRIBUTION.  *
C  *  IP1,IP2,  :  THE NUMBERS INDICATING THE POWERS OF THE THREE    *
C  *  IP3          OVERLAP INTEGRALS IN A TRANSITION INTEGRAL        *
C  *               CONTRIBUTION.                                     *
C  *                                                                 *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  *                                                                 *
C  *  NON STANDARD FORTRAN:-                                         *
C  *  ====================                                           *
C  *                                                                 *
C  *            THIS PROGRAM IS WRITTEN AS A SYSTEM 360 DOUBLE PRE-  *
C  *  CISION PROGRAM.  HOWEVER, ON COMPUTERS WITH A WORD LENGTH OF   *
C  *  48 BITS OR MORE IT SHOULD BE USED IN SINGLE PRECISION ONLY.    *
C  *  IN ORDER TO CONVERT THIS PROGRAM TO A SINGLE PRECISION PROGRAM *
C  *  THE FOLLOWING SHOULD BE DONE:                                  *
C  *                                                                 *
C  *     (1) REMOVE ALL 'IMPLICIT REAL*8(A-H,O-Z)' CARDS.            *
C  *     (2) REMOVE '*8' FROM THE FUNCTION DEFINITION CARDS, FOR     *
C  *         EXAMPLE, 'REAL FUNCTION GRAD*8(I,J)' WILL BE REPLACED   *
C  *         BY 'REAL FUNCTION GRAD(I,J)'                            *
C  *     (3) CHANGE SYSTEM FUNCTION NAMES 'DSQRT','DABS', ETC. TO    *
C  *         'SQRT','ABS', ETC.                                      *
C  *     (4) CHANGE DOUBLE PRECISION CONSTANTS TO SINGLE PRECISION   *
C  *         CONTANTS.                                               *
C  *     (5) CHANGE D-FORMAT CODES TO E-FORMAT CODES.                *
C  *                                                                 *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     ------------------------------------------------------------------
C *** M A I N
C     ------------------------------------------------------------------
C
      program main
c
c     implicit real*8(a-h,o-z)
      integer BLANK,SYMM(5),ASTER
      long ATOM,TRMI,TRMF
      integer EL(50,3)
      long CONFAI(30),CONFBI(30),CONFCI(30),CONFDI(30),
     1     CONFAF(30),CONFBF(30),CONFCF(30),CONFDF(30)
      real WTI(30),WTF(30),ETI,ETF
      real D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1
      COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1,Z,NO,ND
      COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
      common/inf/writef
C
C *** THE INPUT FORMATS.  IT IS TO BE NOTED HERE THAT THE '102 FORMAT'
C *** USED TO INPUT THE NUMERICAL MCHF RADIAL FUNCTION DATA IS NOT THE
C *** SAME AS ACTUALLY USED IN THE MCHF PROGRAM TO PUNCH THESE RADIAL
C *** FUNCTIONS.  THIS HAS BEEN DONE IN ORDER TO KEEP THE DATA-FIELD
C *** LIMITED TO COLUMNS 1-72 AND USE COLUMNS 73-80 OF THE CARDS FOR
C *** IDENTIFICATION AND SEQUENCE NUMBERS ONLY.  HOWEVER, ONE CAN USE
C *** THE MCHF RADIAL FUNCTION DATA AS IT IS PUNCHED OUT IN THE MCHF
C *** PROGRAM BY CHANGING THE '102 FORMAT' TO:
C102  FORMAT(32X,3A1,I6,22X,D16.8/(7F11.7))
C
C
C *** THE OUTPUT FORMATS.
C
 200  FORMAT(5H1SET ,I3/5H CASE,I3/8H *******//)
 201  FORMAT(23X,37HOSCILLATOR STRENGTHS CALCULATION FOR , A8,4H (Z=,I3,
     12H) ,A8,2H->,A8,12H TRANSITION )
 202  FORMAT(1H ,22X,83(1H*)/)
 203  FORMAT(/9H JREAD  =,I3,38H (GENERAL DATA INPUT UNIT FOR THE SET) /
     19H IREAD1 =,I3,45H (WAVEFUNCTION DATA INPUT UNIT FOR THIS CASE)  /
     29H IREAD2 =,I3,53H (TRANSITION INTEGRALS DATA INPUT UNIT FOR THIS
     2CASE)                                                            /
     39H IWRITE =,I3,36H (PRINTED OUTPUT UNIT FOR THIS CASE)           /
     49H NWF    =,I3,35H (NUMBER OF WAVEFUNCTIONS INVOLVED)            /
     59H NCFGI  =,I3,48H (NUMBER OF CONFIGURATIONS IN THE INITIAL STATE)
     6/9H NCFGF  =,I3,46H (NUMBER OF CONFIGURATIONS IN THE FINAL STATE)/
     79H MULT   =,I3,40H (MULTIPLICITY (2S+1) OF THE TWO STATES),  
     8/9h ICTAS  =,I3,20h (PHASE CONVENTIONS)           )
 204  FORMAT(//5H THE ,I2,28H ORBITALS ARE THE FOLLOWING:/)
 205  FORMAT(2X,9(I1,3X),21(I2,2X))
 206  FORMAT(/2X,30(3A1,1X)/)
 207  FORMAT(2X,30(I2,2X))
 208  FORMAT(//38H CONFIGURATIONS/WEIGHTS OF THE INITIAL,A6,2X,
     124HSTATE WITH TOTAL ENERGY= ,F16.8)
 209  FORMAT(  86H -----------------------------------------------------
     1--------------------------------/)
 210  FORMAT(I7,1H.,4A6,F12.7)
 211  FORMAT(//38H CONFIGURATIONS/WEIGHTS OF THE   FINAL,A6,2X,
     124HSTATE WITH TOTAL ENERGY= ,F16.8)
 212  FORMAT(64H1TERM   INITIAL STATE CONFIGURATION    FINAL STATE CONFI
     1GURATION,11X,12HCONTRIBUTION,18X,6HLENGTH,8X,8HVELOCITY)
 213  FORMAT(128H ----   ---------------------------    ----------------
     1---------     -------------------------         ----------     ---
     2-------//)
 214  FORMAT(//91H ERROR ( L(KRHO)-L(KSIG).NE.+1/-1 )  FOUND IN THE TRAN
     1SITION INTEGRAL PART OF THE FOLLOWING,I4,32H-TH INPUT CARD FROM IR
     2EAD2 FILE./)
 215  FORMAT(//50H ERROR (DIFFERENT L VALUES) FOUND IN THE OVERLAP P,
     120hART OF THE FOLLOWING,I4,32H-TH INPUT CARD FROM IREAD2 FILE./)
 216  FORMAT( 2H (,I2,2H) ,I4,1H.,4A6,4H -> ,I2,1H.,4A6,4H : <,3A1,3H/O/
     1,3A1,1H>,17X,1H:,2X,F13.8,2X,F13.8)
 217  FORMAT(1H+,79X,3H* <,3A1,1H/,3A1,3H>**,I1)
 218  FORMAT(1H ,79X,3H* <,3A1,1H/,3A1,3H>**,I1)
 219  FORMAT(/35X,40HFINAL OSCILLATOR STRENGTHS    (LENGTH) =,F14.8/
     1   63X,12H(VELOCITY) =, F14.8)
 220  FORMAT(/42X,35HENERGY DIFFERENCE OF THE STATES =  ,F16.8,5H CM-1/
     174X,3H=  ,F16.8,10H ANGSTROMS,/74X,3H=  ,F16.8,6H A. U.)
 221  FORMAT(//91H ERROR ( L(KRHO)-L(KSIG).NE.+2/-2/0) FOUND IN THE TRAN
     1SITION INTEGRAL PART OF THE FOLLOWING,I4,32H-1H INPUT CARD FROM IR
     1EAD2 FILE./)
 222  FORMAT(//37H THE ORDER OF THE TENSOR OPERATOR IS ,I4,//)
 300  format(//,1x,' jread skal vaere lig 1 eller 5 ')
 301  format(//,1x,' iwrite skal vaere lig enten 6 eller 7 ')
 302  format(//,1x,' iread1 og iread2 skal vaere enten 1, 2 eller 3. ')
C
C *** SET JREAD EQUAL TO THE CARD INPUT CHANNEL OF YOUR INSTALLATION
C *** & NSET=0, AND START THE DATA READING FOR THE FISRT SET OF THE
C *** CALCULATIONS.
C
      zone readf(128,1,stderror)
      zone readf1(128,1,stderror)
      zone readf2(128,1,stderror)
      zone writef(128,1,stderror)
c
c  set input and output channels
      jread=1
      call zassign(readf,jread)
      call open(readf,4,'inf',0)
      iwrite=7
      call zassign(writef,iwrite)
      call open(writef,4,'outf',0)
      iread1=2
      call zassign(readf1,iread1)
      call open(readf1,4,'ouf',0)
      iread2=3
      call zassign(readf2,iread2)
      call open(readf2,4,'pchf',0)
c
      nset=0
c  ' '=2105376, 'S'=5447712, 'P'=5251104, 'D'=4464672,'F'=4595744
c  'G'=4661280
      aster=2105376
      blank=2105376
      symm(1)=5447712
      symm(2)=5251104
      symm(3)=4464672
      symm(4)=4595744
      symm(5)=4661280
      call bldata
c
C
C *** START THE DATA READING FOR A NEW SET OF CALCULATIONS.
C
1     READ(JREAD,100) ASTER,JREAD
100   format(a1,i4)
      if (readerr.shift.(-8).eq.25) goto 15
      if (jread.eq.1.or.jread.eq.5) goto 900
      write(7,300)
      write(6,300)
      goto 15
900   continue
C
C *** IF 'ASTER' IS NOT A BLANK CHARACTER, ALL SETS OF CALCULATIONS ARE
C *** FINISHED AND THE PROGRAM STOPS.
C *** IF 'ASTER' IS A BLANK CHARACTER, SET NSET=NSET+1 & NCASE=0,
C *** AND START THE DATA READING FOR THE FIRST CASE OF THIS SET.
C
      IF(ASTER.NE.BLANK) GOTO 15
      NSET=NSET+1
      NCASE=0
C
C *** START THE DATA READING FOR A NEW CASE OF THE NSET-TH SET.
C
 2    READ(JREAD,101)  ASTER,IREAD1,IREAD2,IWRITE,ATOM,IZ,NWF,TRMI,TRMF,
     1                NCFGI,NCFGF,MULT,ISO,ICTAS
101   format(a1,3i4,a6,2x,2i4,2(2x,a6),5i4)
      if (readerr.shift.(-8).eq.25) goto 15
      if (iwrite.eq.7.or.iwrite.eq.6) goto 901
      write(7,301)
      write(6,301)
      goto 15
901   continue
      if (iread1.le.2.and.iread2.le.3) goto 902
      write(7,302)
      write(6,302)
      goto 15
902   continue
c
c
C
C *** IF 'ASTER' IS NOT A BLANK CHARACTER, ALL CASES OF CALCULATIONS ARE
C *** FINISHED AND THE PROGRAM PROCEEDS TO TAKE THE NEXT SET OF CASES.
C *** IF 'ASTER' IS A BLANK CHARACTER, SET NCASE=NCASE+1, AND PROCEED
C *** FOR READING ADDITIONAL DATA, AND CALCULATIONS ETC. FOR THIS CASE.
C
      IF (ASTER.NE.BLANK) GOTO 1
      NCASE=NCASE+1
      READ(JREAD,102) ETI
102   format(F9.6)
      READ(JREAD,103) (CONFAI(NC),CONFBI(NC),
     1                  CONFCI(NC),CONFDI(NC),WTI(NC),NC=1,NCFGI)
103   format(3a6,a4,F9.6)
      READ(JREAD,104) ETF
104   format(F9.6)
      READ(JREAD,105) (CONFAF(NC),CONFBF(NC),
     1                  CONFCF(NC),CONFDF(NC),WTF(NC),NC=1,NCFGF)
105   format(3a6,a4,F9.6)
      WRITE(IWRITE,200) NSET,NCASE
      Z=IZ
      WRITE(IWRITE,201) ATOM,IZ,TRMI,TRMF
      WRITE(IWRITE,202)
      WRITE(IWRITE,203) JREAD,IREAD1,IREAD2,IWRITE,NWF,NCFGI,NCFGF
     1  ,MULT,ICTAS
C
C *** CALCULATE THE 'R' VARIABLE-MESH FOR THE NUMERICAL RADIAL ORBITALS
C *** P(R).
C
      RHO=-4.0
      DO 3 J=1,NO
      R(J)=EXP(RHO)/Z
      RR(J)=R(J)*R(J)
      R2(J)=SQRT(R(J))
 3    RHO=RHO+H
C
C *** READ THE RADIAL ORBITALS AND CALCULATE THEIR L-QUANTUM NUMBERS.
C
      DO 5 I=1,NWF
      READ(IREAD1,106) (EL(I,J),J=1,3),M,AZ(I)
      READ(IREAD1,107) (P(I,J), J=1,M)
106   format(3a1,i6,22x,F16.8)
107   format(6F11.7)
      J=3
      IF(EL(I,1).NE.BLANK) J=2
      IF(EL(I,J).EQ.SYMM(1)) L(I)=0
      IF(EL(I,J).EQ.SYMM(2)) L(I)=1
      IF(EL(I,J).EQ.SYMM(3)) L(I)=2
      IF(EL(I,J).EQ.SYMM(4)) L(I)=3
      IF(EL(I,J).EQ.SYMM(5)) L(I)=4
      MAX(I)=M
      IF(M.EQ.NO) GO TO 5
      M=M+1
      DO 4 J=M,NO
 4    P(I,J)=D0
 5    CONTINUE
      read(iread2,108) ktens
108   format(i4)
      write(iwrite,222) ktens
      WRITE(IWRITE,204) NWF
      N30=MIN0(NWF,30)
      WRITE(IWRITE,205) (I,I=1,N30)
      WRITE(IWRITE,206) ((EL(I,J),J=1,3),I=1,N30)
      IF(NWF.LE.30) GO TO 6
      WRITE(IWRITE,207) (I,I=31,NWF)
      WRITE(IWRITE,206) ((EL(I,J),J=1,3),I=31,NWF)
 6    WRITE(IWRITE,208) TRMI,ETI
      WRITE(IWRITE,209)
      WRITE(IWRITE,210) (NC,CONFAI(NC),CONFBI(NC),
     1                      CONFCI(NC),CONFDI(NC),WTI(NC),NC=1,NCFGI)
      WRITE(IWRITE,211) TRMF,ETF
      WRITE(IWRITE,209)
      WRITE(IWRITE,210) (NC,CONFAF(NC),CONFBF(NC),
     1                      CONFCF(NC),CONFDF(NC),WTF(NC),NC=1,NCFGF)
      WRITE(IWRITE,212)
      WRITE(IWRITE,213)
C
C *** INITIALIZE THE TOTAL LENGTH (SL) AND THE TOTAL VELOCITY (SV)
C *** 'GF' VALUES TO ZERO, SET NTERM=0, AND START READING THE TRANSITION
C *** INTEGRALS DATA FOR THE CASE IN QUESTION.
C
      SL = D0
      SV = D0
      NTERM=0
C
C *** READ DATA FOR THE NEXT TRANSITION INTEGRAL TERM.
C
 7    READ(IREAD2,109)  ASTER,COEF,KRHO,JI,KSIG,JF,II1,IF1,IP1,
     1                                             II2,IF2,IP2,
     2                                             II3,IF3,IP3
109   format(a1,F12.8,3x,2i2,1x,2i2,1x,9i4)
      if (readerr.shift.(-8).eq.25) goto 14
C
C *** IF 'ASTER' IS NOT A BLANK CHARACTER, ALL TRANSITION INTEGRAL
C *** TERMS HAVE BEEN READ AND USED. THEN THE PROGRAM PROCEEDS TO DO
C *** THE CALCULATIONS OF THE FINAL RESULTS OF OSCILLATOR STRENGTHS
C *** FOR THE CASE IN QUESTION.
C *** IF 'ASTER' IS A BLANK CHARACTER, SET NTERM=NTERM+1, AND PROCEED
C *** TO EVALUATE THE CONTRIBUTIONS TL AND TV OF THIS TERM TO THE
C *** TOTAL 'GF' VALUES SL AND SV.
C
      IF (ASTER.NE.BLANK) GOTO 14
      NTERM=NTERM+1
C
C *** BEFORE PROCEEDING TO USE THIS TRANSITION INTEGRAL INPUT DATA,
C *** CHECK FOR ANY OBVIOUS ERROR IN IT. IF ERROR IS FOUND, IT IS
C *** PRINTED AND THE PROGRAM STOPS.
C
      go to (25,26),ktens
 25   IF(IABS(L(KRHO)-L(KSIG)).EQ.1) GO TO 8
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
      WRITE(IWRITE,214) NTERM
      GO TO 12
 26   if(iabs(l(krho)-l(ksig)).eq.2.or.(l(krho)-l(ksig)).eq.0) go to 8
      write(iwrite,221) nterm
      go to 12
 8    IF(IP1.EQ.0) GO TO 9
      IF(L(II1).EQ.L(IF1)) GO TO 9
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
      GO TO 11
 9    IF(IP2.EQ.0) GO TO 10
      IF(L(II2).EQ.L(IF2)) GO TO 10
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
      GO TO 11
 10   IF(IP3.EQ.0) GO TO 13
      IF(L(II3).EQ.L(IF3)) GO TO 13
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
 11   WRITE(IWRITE,215) NTERM
 12   WRITE(IWRITE,109) ASTER,COEF,KRHO,JI,KSIG,JF,II1,IF1,IP1,
     1                                             II2,IF2,IP2,
     2                                             II3,IF3,IP3
      CALL EXIT
C
C *** NO ERROR IN THE TRANSITION INTEGRAL DATA IS FOUND AND HENCE
C *** PROCEED TO USE IT FOR EVALUATING TL AND TV VALUES.
C
 13   goto(27,28), ktens
 27   k=2
      l1=l(krho)
      L2 = L(KSIG)
      LL = MAX0(L1,L2)
      if(ictas.eq.1) k=ll-l1
      d=(-1)**k*sqrt(ll*d1)*wti(ji)*wtf(jf)*coef
      goto 29
28    l1=l(krho)
      L2=L(KSIG)
      IF(L1.EQ.L2) GOTO 30
      LL=MAX0(L1,L2)
      XX=D3*LL*(LL-D1)
      XX=XX/(D2*(D2*LL-D1))
      d=sqrt(xx)*wti(ji)*wtf(jf)*coef
      GOTO 29
   30 xx=d1*l1*(l1+d1)*(d2*l1+d1)
      XX=XX/(D2*L1+D3)*(D2*L1-D1)
      k=l1+1
      D=(-1)**K*sqrt(XX)*WTI(JI)*WTF(JF)*COEF
C
c
c***  multiply by the overlaps integral if any
c
 29   IF (IP1 .NE. 0) D = D*QUADR(II1,IF1,0)**IP1
      IF (IP2 .NE. 0) D = D*QUADR(II2,IF2,0)**IP2
      IF (IP3 .NE. 0) D = D*QUADR(II3,IF3,0)**IP3
      goto(31,32),ktens
 31   TL = D*QUADR(KRHO,KSIG,1)
      TV = D*GRAD(KRHO,KSIG)
      goto 33
 32   tl=d*quadr(krho,ksig,2)
      tv=d*grad2(krho,ksig)
C
C *** OUTPUT THE CALCULATED TERM CONTRIBUTIONS TL AND TV.
C
 33   WRITE(IWRITE,216) NTERM,JI,CONFAI(JI),CONFBI(JI),CONFCI(JI),
     1      CONFDI(JI),JF,CONFAF(JF),CONFBF(JF),CONFCF(JF),CONFDF(JF),
     2                  (EL(KRHO,J),J=1,3),(EL(KSIG,J),J=1,3),TL,TV
      IF(IP1.NE.0) WRITE(IWRITE,217) (EL(II1,J),J=1,3),(EL(IF1,J),J=1,3)
     1                              ,IP1
      IF(IP2.NE.0) WRITE(IWRITE,218) (EL(II2,J),J=1,3),(EL(IF2,J),J=1,3)
     1                              ,IP2
      IF(IP3.NE.0) WRITE(IWRITE,218) (EL(II3,J),J=1,3),(EL(IF3,J),J=1,3)
     1                              ,IP3
C
C *** ADD THE EVALUATED TERM CONTRIBUTIONS TL AND TV TO THE TOTAL
C *** 'GF' VALUES SL AND SV.
C
      SL = SL + TL
      SV = SV + TV
C
C *** PROCEED TO TAKE THE NEXT TRANSITION INTEGRAL DATA.
C
      GO TO 7
C
C *** ALL TRANSITION INTEGRALS DATA HAVE BEEN USED AND HENCE PROCEED
C *** TO DO THE FINAL CALCULATIONS AND OUTPUT THEM.
C
 14   D = abs(ETI-ETF)
      goto(34,35),ktens
 34   CL = D2*D*abs(float(MULT))/D3
      CV = D2*abs(float(MULT))/(D*D3)
      goto 36
 35   cl=d1*d**3*abs(float(mult))/d30
      cv=d2*d*abs(float(mult))/d15
 36   gl=cl*sl**2
      GV = CV*SV**2
      if(ktens.ne.2) goto 37
      gl=gl*5.325132272e-05
      gv=gv*5.325132272e-05
 37   write(iwrite,219) gl,gv
      DD = D*D2*109737.3e0
      ANGS = D10**8/DD
      WRITE(IWRITE,220) DD,ANGS,D
C
C *** CALCULATIONS FOR THE PRESENT CASE ARE FINISHED AND THE FINAL
C *** RESULTS HAVE ALSO BEEN OUTPUT.
C *** IF THE NEXT CASE TO BE DONE IS A CASE ISO-ELECTRONIC TO THE ONE
C *** JUST FINISHED, REWIND THE DEVICE 'IREAD2'.
C *** THIS WILL ENABLE US TO REUSE THE TRANSITION INTEGRALS DATA.
C *** PROCEED TO TAKE THE NEXT CASE.
C
      goto 15
      if (iso.ne.1) goto 15
      call close(readf2,.true.)
      call open(readf2,4,'pchf',0)
      goto 2
 15   continue
      call exit
      END
c
c-------------------------------------------------------------------
c                         e x i t
c-------------------------------------------------------------------
c
      subroutine exit
      common/inf/writef
      zone writef(128,1,stderror)
4     format(/,a3)
      eof=1644825
      write(writef,4) eof
      call close(writef,.true.)
      stop
      return
      end

C
C     ------------------------------------------------------------------
C *** B L O C K   D A T A
C     ------------------------------------------------------------------
C
      subroutine bldata
c     IMPLICIT REAL*8(A-H,O-Z)
 
      real D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1
      COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,d15,d30,H,H1,Z,NO,ND
      data d0,d1,d2,d3,d4,d5/0.e0,1.e0,2.e0,3.e0,4.e0,.5e0/
      data d6,d10,d15,d30/6.e0,10.e0,15.e0,30.e0/
      DATA H,H1/.0625e0,.4166666666666667e-01/
      DATA NO,ND/220,218/
      return
      END
C
C     ------------------------------------------------------------------
C *** Q U A D R
C     ------------------------------------------------------------------
C
C *** QUADR INTEGRATES P(I)*P(J)*R**KK BY SIMPSON'S RULE
C
      REAL FUNCTION QUADR(I,J,KK)
c     IMPLICIT REAL*8(A-H,O-Z)
      real D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1
      COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,d15,d30,H,H1,Z,NO,ND
      COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
      K = KK + 2
      LI = L(I)
      LJ = L(J)
      DEN = LI + LJ + 1 + K
      ZR = Z*R(4)
      BI = (P(I,4)/(AZ(I)*R2(4)*R(4)**LI) - D1+ZR/(LI+1) )/ZR**2
      BJ = (P(J,4)/(AZ(J)*R2(4)*R(4)**LJ) - D1+ZR/(LJ+1) )/ZR**2
      ALPHA= (D1/(LI + 1) + D1/(LJ + 1))/(DEN + D1)
      ZR = Z*R(1)
      BETA = (DEN+D1)*ALPHA**2 - D2*(BI+BJ+D1/((LI+1)*(LJ+1)))/(DEN+D2)
      D = P(I,1)*P(J,1)*R(1)**K*(((BETA*ZR+ALPHA)*ZR+D1)/(DEN*H1)+D5)
      M = MIN0(MAX(I),MAX(J)) - 1
      DO 1 JJ = 2,M,2
      JP = JJ + 1
    1 D = D + D2*P(I,JJ)*P(J,JJ)*R(JJ)**K+P(I,JP)*P(J,JP)*R(JP)**K
      QUADR = D*H1
      RETURN
      END
C
C     ------------------------------------------------------------------
C *** G R A D
C     ------------------------------------------------------------------
C
C *** THE GRAD FUNCTION SUBPROGRAM COMPUTES THE FOLLOWING DIRECTLY
C ***        <P(J)!D + L(I)/R !P(I)> WITH L(I) > L(J)
C
      REAL FUNCTION GRAD(I,J)
c     IMPLICIT REAL*8(A-H,O-Z)
      real D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1
      COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,d15,d30,H,H1,Z,NO,ND
      COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
      LL = MAX0(L(I),L(J))
      II = I
      JJ = J
      IF ( L(I) .GT. L(J) ) GO TO 1
      II = J
      JJ = I
1     A1 = (LL+D5)/(LL*(LL+1)*(2*LL+1))
      GRAD = R(1)*P(I,1)*P(J,1)*(D1 + A1*Z*R(1))
      DL = D5*P(I,1)*P(J,1)*R(1)
      MM = MIN0(MAX(I)+1,MAX(J)+1,ND)
      K = 2
      F1 = D5*(P(II,K+1) - P(II,K-1))
      F2 = P(II,K+1) - D2*P(II,K) + P(II,K-1)
      G0 = P(JJ,K)*R(K)
      G1 = D5*(P(JJ,K+1)*R(K+1) - P(JJ,K-1)*R(K-1))
      G2 = P(JJ,K+1)*R(K+1) - D2*P(JJ,K)*R(K) + P(JJ,K-1)*R(K-1)
      GRAD = GRAD + D2*F1*G0 +(D2*F2*G1 + F1*G2)/D3
      DL = DL + D2*P(II,K)*P(JJ,K)*R(K) + P(II,K+1)*P(JJ,K+1)*R(K+1)
      DO 2 K = 4,MM,2
      F1 = D5*(P(II,K+1) - P(II,K-1))
      F2 = P(II,K+1) - D2*P(II,K) + P(II,K-1)
      F3 = D5*(P(II,K+2) - P(II,K-2)) - D2*F1
      F4 = P(II,K+2) + P(II,K-2) - D4*(P(II,K+1) + P(II,K-1))
     1   + D6*P(II,K)
      G0 = P(JJ,K)*R(K)
      G1 = D5*(P(JJ,K+1)*R(K+1) - P(JJ,K-1)*R(K-1))
      G2 = P(JJ,K+1)*R(K+1) - D2*P(JJ,K)*R(K) + P(JJ,K-1)*R(K-1)
      G3 = D5*(P(JJ,K+2)*R(K+2) - P(JJ,K-2)*R(K-2)) -D2*G1
      G4 = P(JJ,K+2)*R(K+2) + P(JJ,K-2)*R(K-2) - D4*(P(JJ,K+1)*R(K+1)
     1   + P(JJ,K-1)*R(K-1)) + D6*P(JJ,K)*R(K)
      GRAD = GRAD + D2*F1*G0 +(D2*F2*G1 + F1*G2)/D3
     1   - (F1*G4-F4*G1 + D4*(F2*G3-F3*G2))/90.E0
 2    DL = DL + D2*P(II,K)*P(JJ,K)*R(K) + P(II,K+1)*P(JJ,K+1)*R(K+1)
      GRAD = GRAD + (LL+D5)*DL*H1
      IF (II .EQ. I) GRAD = - GRAD
      RETURN
      END
C
C     ------------------------------------------------------------------
C *** GRAD2
C     ------------------------------------------------------------------
C
C *** THE GRAD2 FUNCTION SUBPROGRAM COMPUTES THE FOLLOWING DIRECTLY
C ***        <P(I)!R.D + F(DELTA)! P(J)>
C
      REAL FUNCTION GRAD2(I,J)
c     IMPLICIT REAL*8(A-H,O-Z)
      real D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1
      COMMON/PARAM/D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1,Z,NO,ND
      COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
      DIMENSION Q(220)
      JJ=I
      II=J
      do 1 k=1,no
 1    Q(K)=P(JJ,K)*R(K)
      LI=L(I)
      LJ=L(J)
      A1=(LI+LJ+D2)/((LI+D1)*(LJ+D1))
      A2=((LJ+D5)*(LJ+D1)+(LJ+D1+D5)*(LI+D1))/((LI+D1)*(LJ+D1))
      A=A1-A2*(LI+LJ+D3)/((LI+LJ+D4)*(LJ+D5))
      FACT=(LJ+D5)/(LI+LJ+D3)
      G=R(1)**2*P(I,1)*P(J,1)*FACT*(1.+A*Z*R(1))
      MM=MIN0(MAX(I)+1,MAX(J)+1,ND)
      K=2
      F1=D5*(P(II,K+1)-P(II,K-1))
      F2=P(II,K+1)-D2*P(II,K)+P(II,K-1)
      G0=Q(K)*R(K)
      G1=D5*(Q(K+1)*R(K+1)-Q(K-1)*R(K-1))
      G2=Q(K+1)*R(K+1)-D2*Q(K)*R(K)+Q(K-1)*R(K-1)
      G=G+D2*F1*G0+(D2*F2*G1+F1*G2)/D3
      do 2 k=4,mm,2
      F1=D5*(P(II,K+1)-P(II,K-1))
      F2=P(II,K+1)-D2*P(II,K)+P(II,K-1)
      F3=D5*(P(II,K+2)-P(II,K-2))-D2*F1
      F4=P(II,K+2)+P(II,K-2)-D4*(P(II,K+1)+P(II,K-1))
     1   +D6*P(II,K)
      G0=Q(K)*R(K)
      G1=D5*(Q(K+1)*R(K+1)-Q(K-1)*R(K-1))
      G2=Q(K+1)*R(K+1)-D2*Q(K)*R(K)+Q(K-1)*R(K-1)
      G3=D5*(Q(K+2)*R(K+2)-Q(K-2)*R(K-2))-D2*G1
      G4=Q(K+2)*R(K+2)+Q(K-2)*R(K-2)-D4*(Q(K+1)*R(K+1)
     1   +Q(K-1)*R(K-1))+D6*Q(K)*R(K)
      G=G+D2*F1*G0+(D2*F2*G1+F1*G2)/D3
     1  -(F1*G4-F4*G1+D4*(F2*G3-F3*G2))/90.E0
 2    CONTINUE
      U=QUADR(JJ,II,0)
      G=G-D5*U
      DELTA=LJ-LI
      IF(DELTA)10,11,12
 10   GRAD2=G-(LI-D2)*U
      RETURN
 11   GRAD2=G+(D1+D5)*U
      RETURN
 12   GRAD2=G+(LI+D3)*U
      RETURN
      END
▶EOF◀