c PSTGICF.f :: Parallel version of stgicf. 23 OCT 19 c D.M. Mitnik, D.C. Griffin, N.R. Badnell c v2.5/4.5 c........multiple-file opening (no global files) c........Now IMODE=-1, only IMODE=1 not present (little used). c C PROGRAM STGICF V4.6 DCG & NRB 13 APR 22 C C PROGRAM TO READ UNPHYSICAL K-MATRICES IN LS COUPLING GENERATED BY C STGF, DETERMINE THE UNPHYSICAL K-MATRICES IN PURE JK COUPLING OR C INTERMEDIATE COUPLING AND FINALLY USE THESE TO DETERMINE AND PRINT C THE COLLISION STRENGTHS BETWEEN INDIVIDUAL LEVELS. THIS PROGRAM C WILL READ K-MATRICES FROM EITHER THE REGULAR R-MATRIX CODES OR THE C NO-EXCHANGE CODES. C PROGRAM MAIN IMPLICIT REAL*8 (A-H,O-Z) C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C C PARAMETER DIMENSIONS C NEMXC--MAXIMUM NUMBER OF ENERGIES IN COARSE MESH READ FROM C STGF RUN C NEMXF--MAXIMUM NUMBER OF ENERGIES IN FINE MESH (IMODE.NE.0) C NMTM --MAXIMUM NUMBER OF TARGET TERMS C NMLT --MAXIMUM NUMBER OF JK OR IC LEVELS C NMLV --MAXIMUM NUMBER OF LEVELS FOR A GIVEN TARGET TERM C NMPL --MAXIMUM NUMBER OF LSPI PARTIAL WAVES C NMPJ --MAXIMUM NUMBER OF JPI PARTIAL WAVES C NMCL --MAXIMUM NUMBER OF CHANNEL IN LSPI K-MATRICES C NMCJ --MAXIMUM NUMBER OF CHANNELS IN JPI K-MATRICES C NMIL --MAXIMUM NUMBER OF CHANNELS FOR A GIVEN TERM FOR A GIVEN C LSPI PARTIAL WAVE C NMIJ --MAXIMUM NUMBER OF CHANNELS FOR A GIVEN LEVEL FOR A GIVEN C JPI PARTIAL WAVE C NMLS --MAXIMUM NUMBER OF LSPI PARTIAL WAVES FOR WHICH THE LS C MATRIX ELEMENTS CAN CONTRIBUTE TO THE JK K-MATRIX ELEMENTS C FOR A GIVEN JPI C NMLSE--MAXIMUM NUMBER OF ENERGIES FOR WHICH THE LS K-MATRICES ARE C HELD IN MEMORY, TO REDUCE I/O. C INJM --MAXIMUM VALUE OF 2J+1 FOR ANY LEVEL C INLM --MAXIMUM VALUE OF L+1 FOR THE CONTINUUM ELECTRON C NMKV --MAXIMUM NUMBER OF NUMBER OF K VALUES FOR A GIVEN PARTIAL C WAVE, A GIVEN TERM, A GIVEN J-VALUE, AND A GIVEN VALUE OF C L OF THE CONTINUUM ELECTRON - THIS SHOULD ALWAYS BE 2 C NMFT --MAXIMUM NUMBER OF DIFFERENT FACTORIALS USED IN COUPLING C COEFFICIENT CALCULATIONS - 800 SHOULD BE ENOUGH C C C I/O UNITS (ALL INPUT UNLESS SPECIFIED AS OUTPUT. HOWEVER, APART FROM C ********* UNIT5, THEY SHOULD ALL HAVE BEEN CREATED AUTOMATICALLY.) C C 1 SCRATCH (DIRECT ACCESS) C 4 RESERVED C 5 PROGRAM CONTROL. FILE dstgicf (THIS OPEN MAYBE COMMENTED-OUT) C 6 PRINTED OUTPUT. FILE routicf C 7 TERM COUPLING COEFFICIENTS & LEVEL LIST. FILE TCCDW.DAT C 8 LS DIPOLE LINE STRENGTHS. FILE strength.dat C 9 BASIC LS DATA. FILE jbinls C 10 DETAILED LS K-MATRIX OUTPUT. FILE kmtls.prnt C 11 DETAILED JK K-MATRIX OUTPUT. FILE kmtjk.prnt C 12 DETAILED IC K-MATRIX OUTPUT. FILE kmtic.prnt C 14 DETAILED JK COLLISION STRENGTH OUTPUT. FILE omgjk.prnt C 19 DETAILED IC COLLISION STRENGTH OUTPUT. FILE omgic.prnt C 13 TERM LIST. FILE term.dat C 15/17 BASIC JK/IC DATA jbinjk/jbinic C 16/18 JK/IC K-MATRICES (UNPHYSICAL). FILE kmtjk/ic.dat COLD 21 LS K-MATRICES (UNPHYSICAL). FILE kmtls.dat CNEW 101,102,... kmtls.001, kmtls.002,... etc by symmetry. C 20 TOTAL COLLISION STRENGTH OUTPUT, JK/IC. FILE OMEGA C 23 TOTAL COLLISION STRENGTH WITHOUT TOP-UP, JK/IC. FILE OMEGANOTOP C 22 ADAS ADF04 TEMPLATE OUTPUT. FILE adasexj.in.form C 40 SCRATCH (SEQUENTIAL) C 41 SCRATCH (SEQUENTIAL) C 42 SCRATCH (SEQUENTIAL) C C LOGICAL BFORM,EX CHARACTER TITLE*80,ELAS*3,PRINT*6,FILNAM*6 C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK3/GAMMA(NMFT) COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK8/ITM(NMPL,NMCL),LC(NMPL,NMCL) COMMON/BLOK9/NCHTNE(NMPL), 1 ISTOTNE(NMPL),LTTNE(NMPL),LPTTNE(NMPL),ITMNE(NMPL,NMCL), 2 LCNE(NMPL,NMCL),INWNE(NMPL),LSPI(NMPL),INW(NMPL,2), 3 ITMLC(NMPL,NMCL),KV(NMPL,NMCL),INC(NMPL,NMCL),INWX(NMPL),NPWNE COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK11/NVEC(NMLT),IVEC(NMTM,NMLT),CVEC(NMTM,NMLT),TOLTCC COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/HIGHPW/INOEXCH,JMNTWO COMMON/BLCHR/TITLE COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT) COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/HOLDLS/RLSHLD(NMLSE*MXRLS+1) COMMON/INTER/RJ1(NOMCJ),RJ2(NOMCJ),IEK0 COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/NRBJKT/ITL(NMLT),FJTL(NMLT),ILLV(NMLT) COMMON/NRBORD/LORDER(NMTM),JORDER(NMLT),NORDER(NMLT),IORDER X,ICORR(NMTM) COMMON/NRBRND/ELEVR(NMLT) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/TOPINF/RMAX(NOMLT),TOPR1,TOPR2,ETOP,ITRAN(NOMLT) X,IEMAX(NOMLT),IOMSET COMMON/NRBWRK/WORK(NMCJ) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX COMMON/WORK1/OMST(NOMLT) C EQUIVALENCE (IDIP,IBIGE) C NAMELIST/STGIC/ITCC,IMODE,IDIP,INOEXCH,JMNTWO,ITOP,IPRTLS,IPRTJK, X IPRTIC,IPTRM,IPTCC,IPTOMJK,IPTOMIC,IBUG,IPRTOP,LRGLAM,LCBE,IBIGE, X IUP,NOMWRT,JMXTWO,TOPR1,TOPR2,ETOP,ELAS,PRINT,IRD0,IJBIN,IPRINT, X TOLTCC, X ieq NAMELIST/MESH1/MXE,E0,EINCR,ieq C c **** parallel **** common/parainfo/iam,nproc include 'mpif.h' common/blockstop/istop character*1 filec,filed,fileu c **** parallel **** C NPOS(I,J,IONE)=((J-1)*(J-2*IONE))/2+I C c **** parallel **** call mpi_init(ierr) call mpi_comm_rank(mpi_comm_world,iam,ierr) call mpi_comm_size(mpi_comm_world,nproc,ierr) istop=0 c ich0 = ichar('0') ic = iam/100 id = (iam - 100*ic)/10 iu = (iam - 100*ic - 10*id) ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) c **** parallel **** C c write(0,*)' BEGIN proc=:',iam OPEN(UNIT=5,FILE="dstgicf",FORM="FORMATTED",STATUS='OLD') c **** parallel **** if (iam.eq.0) then OPEN(UNIT=6,FILE="routicf",FORM="FORMATTED",STATUS='UNKNOWN') OPEN(UNIT=22,FILE="adasexj.in.form",FORM="FORMATTED", 1 STATUS='UNKNOWN') else c **** routicfg etc. for all but processor zero !! OPEN(UNIT=6,FILE="routicfg",FORM="FORMATTED", 1 STATUS='UNKNOWN') OPEN(UNIT=22,FILE="adasexj.in.formg",FORM="FORMATTED", 1 STATUS='UNKNOWN') endif c **** parallel **** C CALL CPU_TIME(TIMEI) C C CALCULATE THE LOGARITHMS OF FACTORIALS AND STORE THEM IN THE C ARRAY GAMMA - USED IN CALCULATION OF 6-J SYMBOLS C GAMMA(1)=TZERO GAMMA(2)=TZERO DO K=3,NMFT J=K-1 AA=J GAMMA(K)=GAMMA(J)+LOG(AA) ENDDO C C GENERAL INPUT C C TITLE(10)--TITLE FOR RUN - 80 CHARACTERS C WRITE(6,199) 199 FORMAT(/34X,'PSTGICF v2.5'/34X,12('*')//) READ(5,100) TITLE 100 FORMAT(A80) WRITE(6,200) TITLE 200 FORMAT(/80('*')//A80//80('*')//) C C ITCC -- IF EQUAL TO ZERO, THE CODE CALCULATES THE COLLISION C STRENGTHS IN PURE JK COUPLING. IF NOT EQUAL TO ZERO, C IT CALCULATES INTERMEDIATE COUPLING COLLISION C STRENGTHS. FOR THIS OPTION THE TERM-COUPLING C COEFFICIENTS MUST BE PRESENT IN THE FILE TCCDW.DAT. C (DEFAULT: 1) C IMODE -- IF EQUAL TO ZERO, THE CODE READS THE LS UNPHYSICAL C K-MATRICES OUTPUT BY STGF, DETERMINES THE UNPHYSICAL C K-MATRICES IN PURE JK COUPLING, AND IF ITCC.NE.0, IN C IC COUPLING, USING THE TERM-COUPLING COEFFICIENTS. C ***THEN IF IJBIN.NE.0 (DEFAULT =0) C IT WRITES THE ASSOCIATED DATA TO FILE JBINJK C (IF ITCC.EQ.0) OR JBINIC (IF ITCC.NE.0) AND THE C K-MATRICES TO THE DIRECT ACCESS FILE KMTJK.DAT C (IF ITCC.EQ.0) OR KMTIC.DAT (IF ITCC.NE.0).*** C CALCULATES THE COLLISION STRENGTHS IN EITHER JK OR IC C COUPLING ON THE ORIGINAL STGF ENERGY MESH AND WRITES C THEM TO THE FILE OMEGA. c c **** parallel **** C NOTE: IMODE>0 IS NOT PRESENTLY IMPLEMENTED. c **** parallel **** c C -- IF LESS THAN ZERO, THE CODES READS THE LS UNPHYSICAL C K-MATRICES OUTPUT BY STGF, DETERMINES THE UNPHYSICAL C K-MATRICES IN PURE JK COUPLING, AND IT ITCC.NE.0 IN IC C COUPLING, USING THE TERM-COUPLING COEFFICIENTS. THEN C IT READS THE ENERGY MESH PARAMETERS IN THE SAME FORMAT C AS WITH THE IMESH=1 MODE USED IN STGF, AND INTERPOLATES C THE UNPHYICAL K-MATRICES IN EITHER JK OR IC COUPLING ON C TO THE FINE ENERGY MESH. FINALLY, IT GENERATES THE C COLLISION STRENGTHS ON THE FINE ENERGY MESH AND WRITES C THEM TO THE OMEGA FILE. THE FINE ENERGY MESH SHOULD C BE MUCH FINER THAN THE COARSE ENERGY MESH USED IN STGF, C THE MINIMUM ENERGY SHOULD BE GREATER THAN OR EQUAL TO C THE MINIMUM ENERGY USED IN THE COARSE ENERGY MESH, AND C THE MAXIMUM ENERGY SHOULD BE LESS THAN OR EQUAL TO C THE MAXIMUM ENERGY USED IN THE COARSE MESH. THIS IS C THE PREFERRED METHOD FOR GENERATING COLLISION STRENGTHS C ON THE FINE ENERGY MESH SINCE IT AVOIDS THE WRITING OF C THE ABSOLUTE FILES KMTJK.DAT OR KMTIC.DAT, WHICH CAN C BECOME VERY LARGE FOR LARGE PROBLEMS. C -- IF GREATER THAN ZERO, THE CODE ASSUMES THAT THERE HAS C BEEN A PREVIOUS IMODE=0 RUN AND THAT JBINJK AND C KMTJK.DAT EXIST IF ITTC.EQ.0 OR JBINIC AND KMTIC.DAT C EXIST, IF ITCC.NE.0. IT THEN READS THE ASSOCIATED DATA C FROM THE FILE JBINJK AND THE UNPHYSICAL K-MATRICES FROM C THE FILE KMTJK.DAT (ITCC.EQ.0) OR THE ASSOCIATED DATA C FROM THE FILE JBINIC AND THE UNPHYSICAL K-MATRICES FROM C THE FILE KMTIC.DAT (ITCC.NE.0) AND CALCULATES THE C COLLISION STRENGTHS ON AN ENERGY MESH TO BE READ IN THE C SAME FORMAT AS WITH THE IMESH=1 MODE USED IN STGF. C A LINEAR INTERPOLATION OF THE UNPHYSICAL K-MATRICES C READ FROM KMTJK.DAT OR KMTIC.DAT IS THEN PERFORMED AND C THE COLLISION STRENGTHS ARE CALCULATED ON THE FINE C ENERGY MESH AND WRITTEN TO THE FILE OMEGA. THE ENERGY C MESH READ IN SHOULD BE MUCH FINER THAN THE COARSE C ENERGY MESH USED IN THE FILE KMTJK.DAT OR THE FILE C KMTIC.DAT, AND THE MINIMUM ENERGY SHOULD BE GREATER C THAN OR EQUAL TO THE MINIMUM ENERGY USED IN THE COARSE C ENERGY MESH AND THE MAXIMUM ENERGY SHOULD BE LESS THAN C OR EQUAL TO THE MAXIMUM ENERGY USED IN THE COARSE MESH. C (DEFAULT: 0) c *** parallel c ieq -- No. of steps to be inserted between coarse stgf mesh c energies, case imode=-1. Default 1 (fine=coarse mesh). c *** parallel C IRD0 -- IF .GE. 50 ASSUMES K-MATRIX FILES SPLIT BY SYMMETRY C kmtls.001, kmtls.002, ... C AND OPENS THEM ON UNIT NUMBERS IRD0+1, IRD0+2 ... C IF .LT. 50 ASSUMES ALL SYMMETRIES IN SAME FILE C (kmtls.dat) AND OPENS UNIT 21. C (Serial DEFAULT: 50) I.E. SPLIT BY SYMMETRY, c parallel=21 i.e. not split by symmetry, only by proc. C IDIP/IBIGE -- IF NOT EQUAL TO ZERO, READ LS DIPOLE LINE STRENGTHS C WITH SIGNS EQUAL TO THE SIGNS OF THE DIPOLE TRANSITION C MATRIX ELEMENTS FROM THE FILE STRENGTH.DAT FOR C THE CALCULATION OF RADIATIVE TRANSITION DATA. C AND OUTPUTS INFINITE ENERGY COLLISION STRENGTH TO C FILE OMEGANOTOP. IN ADDITION, IF .GT. 0 OUTPUTS TO C FILE OMEGA AS WELL. THE LATTER IS NOT RECOMMENDED AS C THIS IS FOR DIPOLE TRANSITIONS ONLY AND BORN DATA IS C DESIRED AS WELL FOR NON-DIPOLE TRANSITIONS. C (DEFAULT: -1) C INOEXCH -- IF NOT EQUAL TO ZERO, THEN READ LS K-MATRICES FROM THE C NO-EXCHANGE CODES; IN THIS CASE, ONE ALSO NEEDS THE C MINIMUM VALUE OF J FOR THE (N+1)-ELECTRON SYSTEM (SEE C THE INPUT OF JMNTWO BELOW) AND THE CODE WILL CHECK TO C SEE THAT JMNTWO IS GREATER THAN OR EQUAL TO ZERO IF C INOEXCH IS NOT EQUAL TO ZERO. C (DEFAULT: 0) C JMNTWO -- THIS IS TWO TIMES THE MINIMUM VALUE OF J FOR THE C (N+1)-ELECTRON SYSTEM, TO BE INCLUDED IN THE STGIC RUN. C IT IS USED FOR A SEPARATE RUN ON THE HIGH PARTIAL-WAVES C WITH EITHER THE NO-EXCHANGE OR EXCHANGE CODES. HERE C THE MINIMUM VALUE OF J SHOULD BE EQUAL TO JMAX+1, WHERE C JMAX IS THE MAXIMUM VALUE OF J FROM THE CORRESPONDING C LOW PARTIAL-WAVE RUN. IF JMNTWO IS GREATER THAN OR C EQUAL TO ZERO, THE CODE WILL AUTOMATICALLY CHECK C TO BE SURE THAT ALL THE LS COUPLING K-MATRICES FROM THE C LS R-MATRIX RUN THAT ARE NEEDED FOR THIS VALUE OF C JMNTWO ARE THERE. THIS REQUIRES THE USER TO DETERMINE C THE LOWEST VALUE OF L FOR THE (N+1)-ELECTRON SYTEM THAT C CAN LEAD TO THIS VALUE OF J; IT WILL MEAN THAT SOME OF C THE SAME VALUES OF L WILL BE NEEDED FOR BOTH THE LOW C AND HIGH PARTIAL-WAVE LS R-MATRIX RUNS. IF JMNTWO C IS GREATER THAN OR EQUAL TO ZERO, THEN THE CODE C WILL NOT WRITE THE FILES (JBINJK AND KMTJK.DAT, IF C ITCC.EQ.0; OR JBINIC AND KMTIC.DAT, IF ITCC.EQ.1) FOR A C LATER IMODE=1 RUN, SINCE THERE IS NO POINT IN C LATER GENERATING THE OMEGA FILE ON A FINE ENERGY MESH C FOR THE HIGH PARTIAL WAVES. NOTE: THE TWO OMEGA FILES C FROM THE LOW AND HIGH PARTIAL-WAVE RUNS CAN BE MERGED C INTO ONE OMEGA FILE USING THE PROGRAM OMADD. C (DEFAULT: -1) C JMXTWO -- MAX VALUE OF J TIMES 2 C (DEFAULT: 999, I.E. ALL) C ITOP -- IF NOT EQUAL TO ZERO USE A SIMPLE GEOMETRIC SERIES C TOP UP PROCEDURE TO TOP UP THE COLLISION STRENGTHS. C IF .GT. 0 THE RATIO IS THAT OF SUCCESSIVE OMEGAS C IF .LT. 0 THE RATIO IS E-FINAL/E-INITIAL. C (SWITCHES TO DEGENERATE ENERGY CASE AS E-RATIO -> 1) C (DEFAULT: 0) HOWEVER, IF IN ADDITION, OR INSTEAD, C LRGLAM -- IF .GE. ZERO THEN BURGESS DIPOLE TOP-UP IS APPLIED AT C 2*JTOT=LRGLAM, AND ITOP IS RESET=-2 IF ZERO. IT IS C REQUIRED THAT THE LS LINE STRENGTHS BE PRESENT IN THE C FILE strength.dat BUT IDIP/IBIGE CAN BE EQUAL TO ZERO. C (DEFAULT: -1). SET LRGLAM LARGE (E.G. 999) TO LET CODE C RESET IT CONSISTENT WITH THE MAX JTOT AVAILABLE. C HOWEVER, IF IN ADDITION C LCBE -- IF .GT. ZERO THEN COULOMB BETHE PARTIALS ARE CALCULATED C FOR 2*JTOT.GE.LCBE, AND THE OMEGA RATIOS CBE/CC PRINTED C IF IPRTOP.NE.0, AND THE TOP-UP IS DONE WITH THE CBE C OMEGAS. IF IPRTOP.NE.2 ONLY LCBE=LRGLAM (OR -1) MAKES C ANY SENSE. SET LCBE.LT.0 TO USE CC OMEGAS INSTEAD. C (DEFAULT: 0 AND LCBE IS RESET TO LRGLAM) C TOPR1,2 ARE THE RATIOS OF (OMEGA+TOP-UP)/OMEGA ABOVE C WHICH DETAILS ARE LOGGED, TOPR1 IS APPLIED TO DIPOLE C TRANSITIONS AND TOPR2 IS APPLIED TO NON-DIPOLE. C ALSO TESTS FOR NEGATIVE DIPOLE TOP-UP. C (DEFAULT: -1.3,-2.0 OFF, NEGATIVE DIPOLE TOP-UP TEST ON) C IOMSET -- IF .EQ. 1 THEN OMEGA(TOT) IS RESET TO SATISFY THE C TOPR1,2 TEST. C (DEFAULT: 0, NO RESETTING OF OMEGAS) C ETOP -- MINIMUM ENERGY AT WHICH THE TOP-UP TESTS ARE APPLIED. C IN AN NX RUN THE TOP-UP MAY BE A LARGE FRACTION OF THE C NX TOTAL BUT "SMALL" WHEN THE EXCHANGE CONTRIBUTION IS C ADDED-IN, ESPECIALLY AT LOW ENERGIES, SO IT MAY BE C PREFERABLE TO SKIP THE TEST THEN. C (DEFAULT: -1.) C IUP -- IN JK-COUPLING (ITCC=0) LEVELS WITHIN A TERM ARE C ORDERED ASCENDING (IUP=1, DEFAULT) OR DESCENDING IN C J (IUP=-1). C NOMWRT -- IF .EQ. 0 NO OMEGAS ARE WRITTEN C -- IF .GT. 0, NOMWRT OMEGAS ARE WRITTEN ROW-WISE, INC ZEROES C -- EDIT SR.OMEG TO DROP ZEROES (DEFAULT, ALL) C -- IF .LT. 0, -NOMWRT OMEGAS BETWEEN OPEN LEVELS ARE WRITTEN, C COLUMN-WISE. C ELAS -- IF .EQ. 'NO', EXCLUDE ELASTIC TRANSITIONS (DEFAULT). C IF .EQ. 'YES' THEN INCLUDE THEM. C IPRINT -- IF .GE. 0 WRITE PROGRESS INFO TO SCREEN (DEFAULT 0) C IPRTLS -- IF NOT EQUAL TO ZERO PRINT K-MATRICES IN LS COUPLING TO C FILE KMTLS.PRNT, ONLY APPLIES WHEN IMODE.EQ.0 C (DEFAULT: 0) C IPRTJK -- IF NOT EQUAL TO ZERO PRINT K-MATRICES IN JK COUPLING TO C FILE KMTJK.PRNT, ONLY APPLIES WHEN IMODE.EQ.0 C (DEFAULT 0). C IPRTIC -- IF NOT EQUAL TO ZERO PRINT K-MATRICES IN IC COUPLING TO C FILE KMTIC.PRNT, ONLY APPLIES WHEN IMODE.EQ.0. C (DEFAULT: 0) C IPTRM -- IF NOT EQUAL TO ZERO PRINT TERMS AND ENERGIES C (DEFAULT: 0) C IPTCC -- IF NOT EQUAL TO ZERO AND ITCC NOT EQUAL TO ZERO, THEN C PRINT LEVELS, ENERGIES, AND TERM COUPLING COEFFICIENTS C ADDITIONALLY, IF SOME OF THE ELEMENTS OF THE TCC'S ARE C ELIMINATED, THE TCC'S WILL BE RE-ORTHONORMALIZED AND C ONE CAN PRINT THE EIGENVECTOR OVERLAPS BEFORE AND AFTER C RE-ORTHONORMALIZATION BY SETTING IPCC NEGATIVE. C ONLY APPLIES WHEN IMODE.EQ.0. C (DEFAULT: 0) C IPTOMJK -- IF EQUAL TO 1 PRINT TOTAL OMEGAS IN JK COUPLING C IF EQUAL TO 2 PRINT PARTIAL AND TOTAL OMEGAS IN JK C COUPLING TO FILE OMGJK.PRNT, ONLY APPLIES WHEN ITCC.EQ.0. C (DEFAULT: 0) C IPTOMIC -- IF EQUAL TO 1 PRINT TOTAL OMEGAS IN IC COUPLING C IF EQUAL TO 2 PRINT PARTIAL AND TOTAL OMEGAS IN IC C COUPLING TO FILE OMGIC.PRNT, ONLY APPLIES WHEN IMODE.EQ.0. C (DEFAULT: 0) C IPRTOP -- IF EQUAL TO 0 PRINT SUMMARY DETAILS OF TOP-UP C IF EQUAL TO 1 PRINT DETAILS OF DIPOLE TOP-UP (SR.TOP1) C IF EQUAL TO 2 PRINT DETAILS OF CC/CBe OMEGA (SR.TOP2) C IF EQUAL TO 3 PRINT DETAILS OF NON-DIPOLE (SR.TOP3) C (DEFAULT: 0) C IBUG -- IF NOT EQUAL TO 0, PRINT J VALUES FOR EACH TERM AND C THE NUMBER OF CHANNELS IN EACH JPI INTERMEDIATE- C COUPLING K-MATRIX. ONLY APPLIES WHEN IMODE.EQ.0. C (DEFAULT: 0) C PRINT -- .EQ. 'FORM' THEN THE OMEGA FILE IS FORMATTED (DEFAULT). C 'UNFORM' THEN AN UNFORMATTED OMEGAU IS CREATED INSTEAD. C C****OPTIONS BELOW ARE FOR EXPERIENCED USERS ONLY: C C TOLTCC -- IF GREATER THAN 0 THEN DROP TERM COUPLING COEFFICIENTS C WHOSE ABS VALUE IS LESS THAN TOLCC, AND REORTHONORMALIZE C THE SET. THIS IS IN ADDITION TO THE ELIMINATION OF C THOSE PRESENT FOR TERMS WHICH WERE INCLUDED IN THE LS CI C EXPANSION ONLY (AND SO MUST BE DROPPED, AUTOMATICALLY); C I.E. IT OPERATES ON THE CC SET. C ITCC=1 IMODE=0 IDIP=-1 INOEXCH=0 JMNTWO=-1 ITOP=0 IPRTLS=0 IPRTJK=0 IPRTIC=0 IPTRM=0 IPTCC=0 IPTOMJK=0 IPTOMIC=0 IBUG=0 IPRTOP=0 LRGLAM=-1 LCBE=0 IUP=1 NOMWRT=9999999 JMXTWO=999 TOPR1=-1.3 TOPR2=-2.0 IOMSET=0 ETOP=-1. ELAS='NO' PRINT='FORM' IPRINT=0 IRD0=21 !Parallel default, do not split by symmetry IJBIN=0 TOLTCC=-ONE ieq=1 !parallel, or read in mesh1 C READ(5,STGIC) c c **** parallel **** if(imode.gt.0) stop 'Only imode.le.0 implemented in parallel code' c **** parallel **** c C WRITE(6,210) IMODE,ITCC,INOEXCH,JMNTWO,JMXTWO,IDIP,IPRTLS,IPRTJK, 1 IPRTIC,IPTRM,IPTCC,IPTOMJK,IPTOMIC,IPRTOP,ITOP,LRGLAM,LCBE,TOPR1 2,TOPR2,ETOP,IOMSET C 210 FORMAT(/' IMODE =',I2,2X,'ITCC =',I2,2X,'INOEXCH =',I2,2X X,'JMNTWO =',I3,2X,'JMXTWO =',I4,2X,'IBIGE =',I2//' IPRTLS =',I2,2X X,'IPRTJK =',I2,2X,'IPRTIC =',I2,2X,'IPTRM =',I2,2X,'IPTCC =' X,I2,2X,'IPTOMJK =',I2,2X,'IPTOMIC =',I2,2X,'IPRTOP =',I2// X' ITOP =',I2,2X,'LRGLAM =',I4,2X,'LCBE =',I4,2X,'TOPR1 =',F4.1 X,2X,'TOPR2 =',F4.1,2X,'ETOP =',F5.2,2X,'IOMSET =',I2///) C IF(IRD0.LT.50)IRD0=21 !SET kmtls.dat UNIT C BFORM=PRINT.EQ.'FORM'.OR.PRINT.EQ.'form' C IF(INOEXCH.NE.0.AND.JMNTWO.LT.0) THEN WRITE(6,205) 205 FORMAT(/'*** ERROR: FOR A NO-EXCHANGE RUN JMNTWO MUST BE SET') STOP 'SET JMNTWO FOR THIS NO-EXCHANGE RUN' END IF C IF(TOPR1*TOPR2.LT.0.0)THEN !NEED BOTH ON OR OFF WRITE(6,*)'*** INPUT ERROR, TOPR1 & TOPR2 MUST BE THE SAME SIGN' STOP '*** INPUT ERROR, TOPR1 & TOPR2 MUST BE THE SAME SIGN' ENDIF IF(IOMSET.EQ.1)THEN IF(TOPR1.LT.0.0)THEN WRITE(6,*)'*** ATTENTION, IOMSET RESET TO ZERO AS TOPR1,2' X ,' ARE NEGATIVE ***' IOMSET=0 ENDIF WRITE(6,195)TOPR1,TOPR2 195 FORMAT(/'*** ATTENTION: IOMSET=1, TOP-UP RATIO LIMITED BY: ' X ,' TOPR1=',F4.1,3X,'TOPR2=',F4.1,' ***'/) ENDIF C IONE=1 IF(ELAS.EQ.'YES')IONE=0 C IF(IDIP.GT.0)IBIGE=1 IF(IDIP.LT.0)IBIGE=-1 IF(LRGLAM.GE.0) THEN IDIPOLE=1 IF(ITOP.EQ.0)ITOP=-2 IF(LCBE.EQ.0)LCBE=LRGLAM IF(LCBE.LT.0)WRITE(6,198) 198 FORMAT(/10X,62('*')//10X,'STRONG WARNING: CC OMEGAS ARE BEING' X ,' USED FOR THE DIPOLE TOP-UP'/25X, X ' RECOMMEND SETTING LCBE=LRGLAM'//10X,62('*')) IF(ITOP.GT.0)WRITE(6,197) 197 FORMAT(/8X,66('*')//8X,'STRONG WARNING: CC OMEGAS ARE BEING' X ,' USED FOR THE NON-DIPOLE TOP-UP'/27X, X ' RECOMMEND SETTING ITOP=-2'//8X,66('*')) ELSE IF(ITOP.NE.0)THEN !DO NOT ALLOW GEOMETRIC SERIES LRGLAM=999 !SWITCH-ON BURGESS DIPOLE TOP-UP IDIPOLE=1 IF(LCBE.EQ.0)LCBE=LRGLAM IF(LCBE.LT.0)WRITE(6,198) IF(ITOP.GT.0)WRITE(6,197) ELSE IF(IDIP.LT.0)IDIP=0 IDIPOLE=IDIP ENDIF END IF IF(LRGLAM.GE.0.AND.LCBE.GT.LRGLAM)THEN WRITE(6,208) 208 FORMAT(/'*** ATTENTION: LCBE HAS BEEN INPUT GREATER THAN LRGLAM' X,' AND SO HAS BEEN RESET TO LRGLAM'/) WRITE(6,*)'**RESET** LCBE NEGATIVE TO SWITCH-OFF USE OF CBE' WRITE(6,*)' ' LCBE=LRGLAM ENDIF C IF(IPRTLS.NE.0) X OPEN(UNIT=10,FILE="kmtls.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") IF(IPRTJK.NE.0) X OPEN(UNIT=11,FILE="kmtjk.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") IF(IPRTIC.NE.0) X OPEN(UNIT=12,FILE="kmtic.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") IF(IPTOMJK.NE.0) X OPEN(UNIT=19,FILE="omgjk.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") IF(IPTOMIC.NE.0) X OPEN(UNIT=14,FILE="omgic.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") C IF(IMODE.LE.0) THEN C C IF IMODE.LE.0: C FIRST OPEN THE FILE JBINLS CONTAINING CHANNEL INFO IN LS COUPLING C THEN CHECK THAT kmtls.001 (and hence kmtls.002 etc) EXISTS. C c **** parallel **** OPEN(UNIT=9,FILE="jbinls"//filec//filed//fileu, + FORM="UNFORMATTED",STATUS='OLD') FILNAM="kmtls." IF(IRD0.GE.50)THEN INQUIRE(FILE=FILNAM//"001."//filec//filed//fileu,EXIST=EX) IF(.NOT.EX)THEN WRITE(6,*)'***ERROR: CANNOT FIND FILE kmtls.001' STOP'***ERROR: CANNOT FIND FILE kmtls.001' ENDIF ELSE INQUIRE(FILE=FILNAM//"dat."//filec//filed//fileu,EXIST=EX) IF(.NOT.EX)THEN WRITE(6,*)'***ERROR: CANNOT FIND FILE kmtls.dat' STOP'***ERROR: CANNOT FIND FILE kmtls.dat' ENDIF OPEN(UNIT=IRD0,FILE=FILNAM//"dat."//filec//filed//fileu, X FORM="UNFORMATTED",STATUS='OLD') ENDIF c **** parallel **** C C IF ITCC.EQ.0: C OPEN THE FILE CONTAINING THE TERMS, LEVELS, AND TERM- C COUPLING COEFFICIENTS IF ITCC.NE.0 OR THE TERM DATA C IF(ITCC.NE.0) THEN OPEN(UNIT=7,FILE="TCCDW.DAT",FORM='FORMATTED', 1 STATUS='UNKNOWN') END IF OPEN(UNIT=13,FILE="term.dat",FORM='FORMATTED', 1 STATUS='UNKNOWN') ENDIF C EINCR=TZERO C C IF IMODE.NE.0: C FIRST READ THE NEW ENERGY MESH IN THE STANDARD IMESH=1 STGF C FORMAT. C IF(IMODE.NE.0) THEN !.gt.0 not get possible in parallel mxe=0 !parallel READ(5,MESH1) IF(MXE.GT.NEMXF) THEN WRITE(6,216) NEMXF,MXE 216 FORMAT(/'**** THE PARAMETER NEMXF MUST BE INCREASED FROM ', 1 I4,' TO ',I4,' ****') STOP 'STOP TO INCREASE THE PARAMETER NEMXF' END IF ENDIF IF(IMODE.GT.0)THEN C C IF IMODE.GT.0: C THEN OPEN THE APPROPRIATE JBIN FILE C IF(ITCC.EQ.0) THEN OPEN(UNIT=17,FILE='jbinjk',FORM="UNFORMATTED",STATUS='OLD') ELSE OPEN(UNIT=17,FILE='jbinic',FORM="UNFORMATTED",STATUS='OLD') END IF END IF C C IF ITOP .GT. 0 OPEN TWO SCRATCH FILES TO STORE PARTIAL OMEGAS C FOR JMAX-1, ODD AND EVEN AND, IF ITOP .NE. 0, ONE FOR SUM TO JMAX-1. C IF(ITOP.GT.0)THEN OPEN(UNIT=40,STATUS='SCRATCH',FORM="UNFORMATTED") OPEN(UNIT=41,STATUS='SCRATCH',FORM="UNFORMATTED") ENDIF IF(ITOP.NE.0)OPEN(UNIT=42,STATUS='SCRATCH',FORM="UNFORMATTED") C C RADIATIVE LINE STRENGTHS C ********* C IF IDIPOLE.NE.0 THEN DETERMINE THE LS DIPOLE TRANSITION C MATRIX ELEMENTS FOR USE IN CALCULATING THE LINE STRENGTHS C IF(IDIPOLE.NE.0) THEN OPEN(UNIT=8,FILE="strength.dat",FORM="FORMATTED" 1 ,STATUS='UNKNOWN') READ(8,102,END=217) NTERMS,NOMTLS 102 FORMAT(2I5) GO TO 218 217 WRITE(6,*)'***ERROR: LINE STRENGTH FILE (strength.dat) MISSING!' STOP '*** ERROR: LINE STRENGTH FILE (strength.dat) MISSING!' C 218 NOMTM0=(NTERMS*(NTERMS-1))/2 IF(NOMTLS.EQ.0)NOMTLS=NOMTM0 IF(ABS(NOMTLS).EQ.NOMTM0)THEN !NO ELASTIC IONELS=1 ELSE NOMTM0=(NTERMS*(NTERMS+1))/2 IF(NOMTLS.NE.ABS(NOMTM0))THEN !ERROR WRITE(6,219)NTERMS, NOMTLS 219 FORMAT('***ERROR: UNABLE TO DECODE LINE STRENGTH FILE' X ,' (strength.dat) FOR NTERMS, NOMTLS= ',2I6) STOP '***ERROR: UNABLE TO DECODE LINE STRENGTH FILE' ENDIF IONELS=0 !ELASTIC IF(IONE.EQ.1)THEN WRITE(6,220) 220 FORMAT('***ERROR: UNABLE TO DECODE LINE STRENGTH FILE' X ,' (strength.dat)'/' RE-RUN STGF WITH' X ,' ELASTIC TRANSITIONS OFF OR SWITCH-ON ELASTIC HERE') STOP '***ERROR: UNABLE TO DECODE LINE STRENGTH FILE' ENDIF ENDIF NOMTM0=(NTERMS*(NTERMS+1-2*IONE))/2 DO N=1,NOMTM0 SSLS(N)=TZERO ENDDO IF(NOMTLS.GT.0)READ(8,103) X ((SSLS(NPOS(I,J,IONE)),J=I+IONELS,NTERMS),I=1,NTERMS-IONELS) IF(NOMTLS.LT.0)READ(8,103) X ((SSLS(NPOS(I,J,IONE)),I=1,J-IONELS),J=1+IONELS,NTERMS) 103 FORMAT(6E12.5) C IF(ITCC.NE.0) THEN DO N=1,NOMTM0 SSIGN=ONE IF(SSLS(N).LT.TZERO) THEN SSIGN=-ONE SSLS(N)=-SSLS(N) END IF SSLS(N)=SSIGN*SQRT(SSLS(N)) END DO ELSE C C REMOVE THE SIGN ON SS C DO N=1,NOMTM0 SSLS(N)=ABS(SSLS(N)) END DO ENDIF ENDIF C C *** IF IMODE.GT.0, READ BASIC DATA FROM JBIN FILE *** C IF(IMODE.GT.0) THEN IF(ITCC.EQ.0) THEN CALL RDJK0(E0,EINCR) GO TO 95 ELSE CALL RDIC0(E0,EINCR) GO TO 95 ENDIF ENDIF C C *** SET-UP TARGET TERMS AND LEVELS *** C CALL TARGET(IUP,IPTRM) C C *** READ PRELIMINARY INFORMATION FROM jbinls *** C C NUMK -- NUMBER OF INCIDENT ENERGIES FOR WHICH LS K-MATRICES WILL C BE READ IN C EMESH(I) -- ENERGY MESH C READ(9) NUMK,NZ,NELC IF(NUMK.GT.NEMXC) THEN WRITE(6,206) NEMXC,NUMK 206 FORMAT(/'**** THE PARAMETER NEMXC MUST BE INCREASED FROM ', 1 I4,' TO ',I4,' ****') STOP 'STOP TO INCREASE THE PARAMETER NEMXC' END IF AZ=MAX(NZ-NELC,1) AZSQ=AZ*AZ READ(9) (EMESH(I),I=1,NUMK) C C SET UP THE NEW ENERGY MESH C IF(IMODE.LT.0)THEN CSER EMIN=EMESH(1) CSER EMAX=EMESH(NUMK) CSER IEE=0 CSER EM=E0-EINCR CSER EINCRH=EINCR/TWO CSER DO IE=1,MXE CSER EM=EM+EINCR CSER IF(EM.GT.EMIN-EINCRH.AND.EM.LT.EMAX+EINCRH) THEN CSER IEE=IEE+1 CSER EMSH(IEE)=EM CSER END IF CSER END DO CSER MXE=IEE c c***parallel c if(mxe.gt.0.and.iam.eq.0)then write(*,*)'***ATTENTION: user input mxe (& e0, eincr) are' x ,' ignored by interpolation - only ieq is used, along with' x ,' the PSTGF mesh' write(6,*)'***ATTENTION: user input mxe (& e0, eincr) are' x ,' ignored by interpolation - only ieq is used, along with' x ,' the PSTGF mesh' endif c ieq=iabs(ieq) !incase user sets <0 as stgf eincc=emesh(2)-emesh(1) ncint=1 do i=2,numk !=nproce einct=emesh(i)-emesh(i-1) if(einct.gt.eincc*1.1)then ncint=ncint+1 endif enddo if(ncint.eq.1.and.ieq.gt.1.and.iam.eq.0)then write(*,*)'****WARNING: interpolation specified but ncint=1' x ,' was nseq and/or ncint used specify a suitable stgf mesh?' write(6,*)'****WARNING: interpolation specified but ncint=1' x ,' was nseq and/or ncint used specify a suitable stgf mesh?' endif c nseq=numk/ncint if(nseq*ncint.ne.numk)stop 'mis-match on energy mesh' eincr=eincc/ieq nseq=nseq-1 nseqf=ieq*nseq mxe=nseqf*ncint if(mxe.gt.nemxf) then write(6,216) nemxf,mxe stop 'stop to increase the parameter nemxf' end if c c write(70+iam,*)ncint,eincc,eincr,nseqf c iee=0 e0=emesh(1) eint=nproc*nseq*eincc do i=1,ncint do j=1,nseq ec=e0+(j-1)*eincc do k=1,ieq iee=iee+1 em=ec+(k-1)*eincr emsh(iee)=em enddo enddo e0=e0+eint enddo if(iam+1.eq.nproc)mxe=mxe-ieq+1 c c write(60+iam,*)iam,mxe,iee c do i=1,mxe c write(60+iam,*)iam,i,emsh(i) c enddo c***parallel c ENDIF C C WRITE PRELIMINARY INFORMATION TO jbinjk/ic C IF(IJBIN.NE.0) THEN C C IF ITCC.EQ.0 OPEN FILE JBINJK AND WRITE PRELIMINARY DATA TO IT C IF(ITCC.EQ.0) THEN OPEN(UNIT=15,FILE='jbinjk',FORM='UNFORMATTED', 1 STATUS='UNKNOWN') WRITE(15) NTLS,NUMK,IBIGE,NZ,NELC WRITE(15) (ENT(I),NJV(I),I=1,NTLS) DO I=1,NTLS WRITE(15) (FJT(I,IJ),IJ=1,NJV(I)) END DO WRITE(15) (EMESH(I),I=1,NUMK) ELSE C C ELSE OPEN FILE JBINIC AND WRITE PRELIMINARY DATA TO IT C OPEN(UNIT=15,FILE='jbinic',FORM='UNFORMATTED', 1 STATUS='UNKNOWN') WRITE(15) NTLS WRITE(15) (ENT(I),I=1,NTLS) WRITE(15) NLEV,NUMK,IBIGE,NZ,NELC WRITE(15) (ELEV(I),JLEV(I),JPLEV(I),I=1,NLEV) WRITE(15) (EMESH(I),I=1,NUMK) END IF END IF C C IF FINE MESH EQUALS COARSE MESH, THEN TRANSFER C IF(IMODE.EQ.0)THEN MXE=NUMK DO I=1,NUMK EMSH(I)=EMESH(I) ENDDO ncint=1 nseq=mxe nseqf=nseq ENDIF C C WRITE PRELIMINARY DATA TO OMEGA FILE C 95 IF(NOMWRT.NE.0)THEN NOMT=NLEV*(NLEV+1-2*IONE)/2 NOMT=MIN(NOMT,ABS(NOMWRT)) IF(NOMWRT.LT.0)NOMT=-NOMT NOMWRT=NOMT DO I=1,NLEV EINL(I)=ELEV(I)/AZSQ NCHJX(I)=(JLEV(I)+1)*(1-2*JPLEV(I)) ENDDO C MXE0=MXE IF(IBIGE.GT.0)MXE0=MXE0+1 MXE1=MXE+IABS(IBIGE) C c **** parallel **** IF(BFORM)THEN OPEN(UNIT=20,FILE="OMEGA"//filec//filed//fileu, + FORM="FORMATTED",STATUS='UNKNOWN') WRITE(20,*) NZ,NELC WRITE(20,*) NLEV,MXE0,NOMWRT WRITE(20,*) (' 0',NCHJX(I),I=1,NLEV) WRITE(20,270) (EINL(I),I=1,NLEV) 270 FORMAT(1P5E16.6) IF(ITOP.NE.0)THEN OPEN(UNIT=23,FILE="OMEGANOTOP"//filec//filed//fileu, + FORM="FORMATTED",STATUS='UNKNOWN') WRITE(23,*) NZ,NELC WRITE(23,*) NLEV,MXE1,NOMWRT WRITE(23,*) (' 0',NCHJX(I),I=1,NLEV) WRITE(23,270) (EINL(I),I=1,NLEV) ENDIF ELSE OPEN(UNIT=20,FILE="OMEGAU"//filec//filed//fileu, + FORM="UNFORMATTED",STATUS='UNKNOWN') IZERO=0 WRITE(20) NZ,NELC WRITE(20) NLEV,MXE0,NOMWRT WRITE(20) (IZERO,NCHJX(I),I=1,NLEV) WRITE(20) (EINL(I),I=1,NLEV) IF(ITOP.NE.0)THEN OPEN(UNIT=23,FILE="OMEGANOTOPU"//filec//filed//fileu, + FORM="UNFORMATTED",STATUS='UNKNOWN') WRITE(23) NZ,NELC WRITE(23) NLEV,MXE1,NOMWRT WRITE(23) (IZERO,NCHJX(I),I=1,NLEV) WRITE(23) (EINL(I),I=1,NLEV) ENDIF ENDIF c **** parallel **** C INITIALIZE OMEGA CALL OMEG(0,EMSH,MXE,EINL,NLEV,OMEGA,NOMWRT,IONE) ENDIF C C *** IF IMODE.LE.0, READ BASIC DATA FROM JBINLS FILE *** C IF(IMODE.LE.0)THEN IF(INOEXCH.EQ.0)CALL RDLS IF(INOEXCH.NE.0)CALL RDLSNE CALL INIT(JMXTWO) ENDIF C IF(IJBIN.NE.0) WRITE(15)FJTMAX IF(IMODE.GT.0) READ(17)FJTMAX JTMAX2=2*FJTMAX+0.1 C C RESET LRGLAM AND LCBE IF NECESSARY C IF(LRGLAM.GT.JTMAX2-2)THEN LRGLAM=JTMAX2-2 WRITE(6,207)LRGLAM 207 FORMAT(/'***ATTENTION: LRGLAM HAS BEEN INPUT GREATER/EQUAL TO ' X,'(TWICE) THE MAX JTOT FOR THIS RUN AND SO HAS BEEN RESET TO ', XI3/) WRITE(6,*)'**RESET** LRGLAM NEGATIVE TO SWITCH-OFF DIPOLE TOP-UP' WRITE(6,*)' ' ENDIF IF(LCBE.GT.JTMAX2-2)THEN LCBE=JTMAX2-2 WRITE(6,204)LCBE 204 FORMAT(/'***ATTENTION: LCBE HAS BEEN INPUT GREATER/EQUAL TO ' X,'(TWICE) THE MAX JTOT FOR THIS RUN AND SO HAS BEEN RESET TO ', XI3/) WRITE(6,*)'**RESET** LCBE NEGATIVE TO SWITCH-OFF USE OF CBE' WRITE(6,*)' ' ENDIF C C FINDING STARTING POINT ON COARSE MESH (serial only) C (WE ASSUME THAT THE FINE MESH IS NOT COARSER THAN THE COARSE!) C IE0=1 IEK0=0 IF(IMODE.LT.0)then IE0=0 iek0=1 endif EMSH(0)=999999.0 EMESH(0)=-1. EINCRH=EINCR/TWO C CSER IF(IMODE.NE.0)THEN CSER DO I=1,NUMK-1 CSER IF(EMSH(1).GT.EMESH(I)-EINCRH.AND. CSER X EMSH(1).LT.EMESH(I+1)+EINCRH)THEN CSER IEK0=I CSER GO TO 89 CSER ENDIF CSER ENDDO CSER STOP "MAIN: CAN'T FIND REQUESTED ENERGY..." CSER ENDIF CSER89 continue C C *** LOOP OVER JP SYMMETRIES *** C DO IWJ=1,NPWJ IF(IMODE.LE.0)THEN CALL INITJK(IWJ) IF(ITCC.NE.0)CALL INITIC(IWJ) ELSE IF(ITCC.EQ.0)CALL READJK(IWJ) IF(ITCC.NE.0)CALL READIC(IWJ) ENDIF c write(70+iam,8889)iwj c8889 format(2x,'symmetry=: ',i4) C C *** LOOP OVER INCIDENT ELECTRON ENERGIES *** C IEE=0 do nc=1,ncint IEK=IEK0+(nc-1)*(nseq+1) IF(IMODE.LT.0)IEK=IEK-1 eincrh=abs(eincrh) DO ie=IE0,nseqf IEE=IEE+1 c write(70+iam,*) iam,ie,iee,emsh(iee) c8888 format(2x,'iee=: ',i5,' energy=: ',1pg16.4,' proc=: ',i4) ISW=0 IF(EMSH(IEE).GT.EMESH(IEK)-EINCRH)THEN IEK=IEK+1 c write(70+iam,*)iee,iek,emesh(iek) IF(IMODE.LE.0)THEN ISW=1 IF(INOEXCH.EQ.0)CALL RDKMTX(IEK,IWJ,FILNAM,IRD0) IF(INOEXCH.NE.0)CALL RDKMNE(IEK,IWJ,FILNAM,IRD0) ENDIF ENDIF IF(ISW.EQ.1)THEN CALL KMTRJK(IEK,IWJ) IF(ITCC.NE.0)CALL KMTRIC(IEK,IWJ) ENDIF IF(ie.GT.0)THEN IF(IMODE.NE.0)THEN IF(ITCC.EQ.0)CALL INTERJK(IEE,EINCR,IWJ,IEK) IF(ITCC.NE.0)CALL INTERIC(IEE,EINCR,IWJ,IEK) ENDIF IF(ITCC.EQ.0)CALL OMGJK(IEE,IWJ,BFORM) IF(ITCC.NE.0)CALL OMGIC(IEE,IWJ,BFORM) ENDIF IF(ie.GT.0)EINCRH=-ABS(EINCRH) if(ie.eq.0)iee=iee-1 if(iee.eq.mxe)go to 90 ENDDO enddo C C *** END LOOP OVER INCIDENT ELECTRON ENERGIES *** C 90 EINCRH=ABS(EINCRH) C IF(ITCC.EQ.0)NNN=NCHPJ(IWJ) IF(ITCC.NE.0)NNN=NCHPJX(IWJ) C WRITE(6,276) IWJ,NNN,FJTT(IWJ),JPTT(IWJ) IF(iam.eq.0.and.IPRINT.GE.0) X WRITE(0,276) IWJ,NNN,FJTT(IWJ),JPTT(IWJ) 276 FORMAT('COMPLETED SYMMETRY NO.',I3,' NCHAN =',I4, X' WITH JPI =',F4.1,I3) CALL FLUSH(6) C ENDDO C C *** END LOOP OVER JP SYMMETRIES *** C C WRITE(6,277) JTMAX2 277 FORMAT(/5X,'THE MAXIMUM TOTAL J (X2) FROM THIS RUN IS:' X ,I3/48X,'**'/) C CALL CPU_TIME(TIMEF) TIME=TIMEF-TIMEI C TIME=TIME/60.0 WRITE(6,290) TIME,nproc 290 FORMAT(//1X,'CPU TIME=',F9.3,' MIN -- processors=:',i4) CALL FLUSH(6) CALL FLUSH(20) !FLUSH OMEGXXX FOR TIMESTAMP c c **** parallel **** call mpi_finalize(ierr) c **** parallel **** STOP END C C ****************************************************************** C REAL*8 FUNCTION ARGAM(L,A) C IMPLICIT REAL*8 (A-H,O-Z) C C CALCULATES ARGGAMMA(L+1+I*A) C WHERE L IS AN INTEGER NOT LESS THAN ZERO C B=ABS(A) B=250.0D0*B**0.25D0-A*A J0=L+1 C=J0 D=C*C Z=0.0D0 IF(D -B)1,6,6 1 B=SQRT (B) J1=B DO 5 J=J0,J1 D=J D=A/D G1=ABS(D) IF(G1-0.1D0)2,3,3 2 G1=D*D G2=-35.0D0*G1+45.0D0 G2=-G1*G2+63.0D0 G2=-G1*G2+105.0D0 G1=D -D*G1*G2/315.0D0 GO TO 4 3 G1=ATAN (D) 4 Z=Z+G1 5 CONTINUE J0=J1+1 6 D=J0 G0=D*D U=A*A G1=1.0D0/(G0+U) G2=G1*G1 G3=10.0D0*G0*G0-20.0D0*G0*U+2.0D0*U*U G3=G3*G2-21.0D0*G0+7.0D0*U G3=G3*G2+210.0D0 G1=A*G3*G1/2520.0D0 ARGAM=-Z+0.5D0*A*LOG(G0+U)+(D -0.5D0)*ATAN(A/D)-A-G1 RETURN END C C ****************************************************************** C SUBROUTINE DIPIC IMPLICIT REAL*8 (A-H,O-Z) C C CALCULATE THE DIPOLE LINE STRENGTHS BETWEEN LEVELS IN IC C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) PARAMETER (EINF=1.0D6) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK11/NVEC(NMLT),IVEC(NMTM,NMLT),CVEC(NMTM,NMLT),TOLTCC COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX C NPOS(I,J,IONE)=((J-1)*(J-2*IONE))/2+I IPOS(I,J,IONE,N)=J+N*(I-1)-(I*(I-1+2*IONE))/2 C NOMLT0=(NLEV*(NLEV+1-2*IONE))/2 DO N=1,NOMLT0 IDIPX(N)=0 SS(N)=TZERO ENDDO C IF(ITOP.NE.0)THEN !DETERMINE TRANSITION MULTIPOLE N=0 DO ILF=1+IONE,NLEV IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 IF(JLEV(ILI)+JLEV(ILF).NE.0)THEN JDIF=IABS(JLEV(ILI)-JLEV(ILF))/2 JSUM=(JLEV(ILI)+JLEV(ILF))/2 IF(JPLEV(ILI).EQ.JPLEV(ILF))THEN !EVEN IF(JDIF.EQ.0)JDIF=2 IF((-1)**JDIF.LT.0)JDIF=JDIF+1 ELSE !ODD IF(JDIF.LE.1)JDIF=3 !CASE SS(N)=0. IF((-1)**JDIF.GT.0)JDIF=JDIF+1 ENDIF IF(JDIF.GT.JSUM)JDIF=0 IDIPX(N)=JDIF ENDIF ENDDO ENDDO ENDIF C IF(IDIPOLE.EQ.0)RETURN C LRGLMN=1000 INDM=0 N=0 DO ILF=1+IONE,NLEV FJF=JLEV(ILF) FJF=FJF/TWO IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 SS(N)=TZERO IF(JPLEV(ILF).EQ.JPLEV(ILI)) GO TO 412 IF(JLEV(ILF)+JLEV(ILI).EQ.0.OR.ABS(JLEV(ILF)-JLEV(ILI)).GT.2) 1 GO TO 412 C FJI=JLEV(ILI) FJI=FJI/TWO SUM=TZERO DO IVKI=1,NVEC(ILI) ITI=IVEC(IVKI,ILI) FLTI=LT(ITI) FSPI=LST(ITI) FSPI=(FSPI-ONE)/TWO DO IVKF=1,NVEC(ILF) ITF=IVEC(IVKF,ILF) FLTF=LT(ITF) FSPF=LST(ITF) FSPF=(FSPF-ONE)/TWO IF(FSPF.EQ.FSPI) THEN PHASE=ONE IF(ENT(ITI).GT.ENT(ITF)) THEN IF(MOD(FLTI+FSPI+FJF+ONE,TWO).NE.TZERO) PHASE=-ONE ELSE IF(MOD(FLTF+FSPI+FJF,TWO).NE.TZERO) PHASE=-ONE END IF ITXF=MAX(ITI,ITF) ITXI=MIN(ITI,ITF) NT=NPOS(ITXI,ITXF,IONE) DLINE=PHASE*SQRT((TWO*FJI+ONE)*(TWO*FJF+ONE))* 1 S6J(FLTI,FSPI,FJI,FJF,ONE,FLTF)/SQRT(TWO*FSPI+ONE) SUM=SUM+CVEC(IVKI,ILI)*DLINE*SSLS(NT)*CVEC(IVKF,ILF) ENDIF ENDDO ENDDO SS(N)=SUM**2 C C IF(AZSQ*SS(N).GT.1.D-7)THEN IF(SS(N).GT.1.D-10)THEN IF(LRGLAM.GE.0) THEN INDM=INDM+1 NTOP(INDM,1)=ILI NTOP(INDM,2)=ILF LLM=MIN(JLEV(ILI),JLEV(ILF)) LRGLMN=MIN(LRGLMN,LRGLAM-2*LLM) LLM=-LLM+1 LITLAM(INDM)=(LRGLAM+LLM)/2 ENDIF ELSE IF(IPRTOP.EQ.1)WRITE(6,*)'NOTE, LINE STRENGTH SET TO ZERO' X ,' FOR: ILI=',ILI,' ILF=',ILF,' SS=',SS(N) SS(N)=TZERO ENDIF C 412 IF(SS(N).NE.TZERO)IDIPX(N)=1 !OVERWRITE DIPOLE IF(NOMWRT.LT.0)IPT=NPOS(ILI,ILF,IONE) IF(NOMWRT.GT.0)IPT=IPOS(ILI,ILF,IONE,NLEV) OMGINF(IPT)=(FOUR/THREE)*SS(N)*LOG(EINF*AZSQ) ENDDO ENDDO C RETURN END C C ****************************************************************** C SUBROUTINE DIPJK IMPLICIT REAL*8 (A-H,O-Z) C C CALCULATE THE DIPOLE LINE STRENGTHS BETWEEN LEVELS IN PURE JK C COUPLING C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) PARAMETER (EINF=1.0D6) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX COMMON/NRBJKT/ITL(NMLT),FJTL(NMLT),ILLV(NMLT) C NPOS(I,J,IONE)=((J-1)*(J-2*IONE))/2+I IPOS(I,J,IONE,N)=J+N*(I-1)-(I*(I-1+2*IONE))/2 C NOMLT0=(NLEV*(NLEV+1-2*IONE))/2 DO N=1,NOMLT0 IDIPX(N)=0 SS(N)=TZERO ENDDO C IF(ITOP.NE.0)THEN !DETERMINE TRANSITION MULTIPOLE N=0 DO ILF=1+IONE,NLEV ITF=ITL(ILF) IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 ITI=ITL(ILI) IF(LT(ITI)+LT(ITF).NE.0.AND.LST(ITI).EQ.LST(ITF))THEN LDIF=IABS(LT(ITI)-LT(ITF)) LSUM=LT(ITI)+LT(ITF) IF(LPT(ITI).EQ.LPT(ITF))THEN !EVEN IF(LDIF.EQ.0)LDIF=2 IF((-1)**LDIF.LT.0)LDIF=LDIF+1 ELSE !ODD IF(LDIF.LE.1)LDIF=3 !CASE SS(N)=0. IF((-1)**LDIF.GT.0)LDIF=LDIF+1 ENDIF IF(LDIF.GT.LSUM)LDIF=0 IDIPX(N)=LDIF ENDIF ENDDO ENDDO ENDIF C IF(IDIPOLE.EQ.0)RETURN C LRGLMN=1000 INDM=0 N=0 DO ILF=1+IONE,NLEV ITF=ITL(ILF) FLTF=LT(ITF) FSPF=LST(ITF) FSPF=(FSPF-ONE)/TWO IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 ITI=ITL(ILI) FLTI=LT(ITI) IF(LST(ITI).NE.LST(ITF).OR.LPT(ITI).EQ.LPT(ITF)) THEN SS(N)=TZERO ELSE NT=NPOS(ITI,ITF,IONE) SS(N)=(TWO*FJTL(ILI)+ONE)*(TWO*FJTL(ILF)+ONE)* 1 (S6J(FLTI,FSPF,FJTL(ILI),FJTL(ILF),ONE,FLTF)**2) 2 *SSLS(NT)/(TWO*FSPF+ONE) END IF C C IF(AZSQ*SS(N).GT.1.D-7)THEN IF(SS(N).GT.1.D-10)THEN IF(LRGLAM.GE.0)THEN INDM=INDM+1 NTOP(INDM,1)=ILI NTOP(INDM,2)=ILF LLM=MIN(JLEV(ILI),JLEV(ILF)) LRGLMN=MIN(LRGLMN,LRGLAM-2*LLM) LLM=-LLM+1 LITLAM(INDM)=(LRGLAM+LLM)/2 ENDIF ELSE IF(IPRTOP.EQ.1)WRITE(6,*)'NOTE, LINE STRENGTH SET TO TZERO' X ,' FOR: ILI=',ILI,' ILF=',ILF,' SS=',SS(N) SS(N)=TZERO ENDIF C IF(SS(N).NE.TZERO)IDIPX(N)=1 !OVERWRITE DIPOLE IF(NOMWRT.LT.0)IPT=NPOS(ILI,ILF,IONE) IF(NOMWRT.GT.0)IPT=IPOS(ILI,ILF,IONE,NLEV) OMGINF(IPT)=(FOUR/THREE)*SS(N)*LOG(EINF*AZSQ) ENDDO ENDDO C RETURN END C C*********************************************************************** C REAL*8 FUNCTION F21(A,B,C,D,EPS,IFAIL) IMPLICIT REAL*8 (A-H,O-Z) C IPRINT=IFAIL IFAIL=0 T=(A*B*D)/C DD=1.0D0/(1.0D0-D) SUM=1.0D0+T TN1=0.0D0 I=1 3 AI=I T=T*(A+AI)*(B+AI)*D/((C+AI)*(1.0D0+AI)) TN2=T*DD F21=SUM+TN2 SUM=SUM+T AT=ABS(T+TN2-TN1) AS=ABS(F21)*EPS IF(AS-AT)1,2,2 1 TN1=TN2 I=I+1 IF(I-300)3,3,4 4 IF(IPRINT.GT.0)WRITE(6,100) IFAIL=3 100 FORMAT(' FAILED TO CONVERGE IN F21') 2 RETURN END C C********************************************************** C REAL*8 FUNCTION FDIP(EK1,L1,EK2,L2,IFAIL) C IMPLICIT REAL*8 (A-H,O-Z) C C ALAN BURGESS DEPT. OF APPLIED MATHS. AND THEORETICAL PHYSICS,CAMBRIDGE C CALCULATES THE FUNCTION I(KAPPA1,L1,KAPPA2,L2,1) DEFINED IN PHIL. C TRANS. ROY. SOC. A226,255,1970, WHERE EK1=KAPPA1**2 AND EK2=KAPPA2**2. C IT IS SUITABLE FOR USE IN EQUATIONS (8),(9),(10) OR (11) OF C J. PHYS. B. 7,L364,1974. C NRB - IFAIL DATA EPS/1.D-4/ C IPRINT=IFAIL IF(EK1+EK2-1.0D-40) 11,11,12 11 FDIP=0.0D0 IFAIL=1 IF(IPRINT.EQ.-2)WRITE(6,100)IFAIL RETURN 12 IF(EK1-EK2) 1,1,2 1 EMIN=EK1 EMAX=EK2 GO TO 3 2 EMIN=EK2 EMAX=EK1 3 T=EMIN/EMAX IF(T-0.02944D0) 4,4,5 4 FDIP=FDIP1(EK1,L1,EK2,L2) GO TO 9 5 IF(T-0.16667D0) 7,6,6 6 FDIP=FDIP2(EK1,L1,EK2,L2) GO TO 9 7 FDIP=FDIP1(EK1,L1,EK2,L2) IF(FDIP*FDIP-1.0D-40) 6,6,8 8 IF(FDIP.LT.0.0.OR.FDIP.GT.1.)THEN IFAIL=3 IF(IPRINT.EQ.-2)WRITE(6,100)IFAIL FDIP=0.0 RETURN ENDIF FA=FDIPA(EK1,L1,EK2,L2) IFAIL=0 IF(FA.EQ.0.0)THEN FA=FDIP0(EK1,L1,EK2,L2,EPS,IFAIL) IFAIL=-IFAIL IF(FA.EQ.0.0)RETURN ENDIF RAT=FDIP/FA IF(RAT.GT.10.)THEN IFAIL=4 IF(IPRINT.EQ.-2)WRITE(6,100)IFAIL FDIP=0.0 ENDIF RETURN 9 IF(FDIP*FDIP-1.0D-40) 10,10,8 10 IFAIL=2 IF(IPRINT.EQ.-2)WRITE(6,100)IFAIL RETURN 100 FORMAT('***FDIP FAILURE: IFAIL=',I2) END C C*********************************************************************** C REAL*8 FUNCTION FDIP0(EK1,L1,EK2,L2,EPS,IFAIL) C IMPLICIT REAL*8 (A-H,O-Z) C C ALAN BURGESS,DEPT OF APPLIED MATHS. AND THEORETICAL PHYSICS,CAMBRIDGE C CALCULATES THE FUNCTION I0(K1,L1,K2,L2,1) DEFINED IN PHIL. TRANS. C ROY. SOC. A266,255,1970, WHERE EK1=K1*K1, EK2=K2*K2, AND THE RELATIVE C ACCURACY IS APPROXIMATELY EPS. C IT IS SUITABLE FOR USE IN EQUATIONS (13) ETC. OF J.PHYS.B. 7,L364,1974 C NRB - IFAIL C IPRINT=IFAIL IFAIL=0 IF(L1-L2)1,2,4 1 L=L1 GO TO 5 2 IF(IPRINT.EQ.-2)WRITE(6,100)L1 IFAIL=1 FDIP0=0.0D0 3 RETURN 4 L=L2 5 EL=L FDIP0=0.5D0/(EL+1.0D0) IF(EK1-EK2)6,3,7 6 E=EK1/EK2 P=L1-L GO TO 8 7 E=EK2/EK1 P=L2-L 8 FDIP0=FDIP0*E**((EL+P+0.5D0)*0.5D0) C TO OBTAIN THE FUNCTION EK1 OF M.J. SEATON, PROC. PHYS. SOC. A68,457, C 1955, REMOVE THE 'C' ON THE NEXT LINE. C FDIP0=1.0D0 IF(E -0.5D0)21,20,20 20 P1=P-0.5D0 T=P1*(EL+1.0D0)*(E -1.0D0) I0=L+1 H0=0.0D0 DO 9 I=1,I0 TI=I H0=H0+1.0D0/TI 9 CONTINUE X=1.0D0-E H=1.0D0-(P+P+H0+LOG(0.25D0*X)) S=1.0D0+T*H A=EL+1.0D0 B=P1 C=1.0D0 D=0.0D0 10 A=A+1.0D0 B=B+1.0D0 C=C+1.0D0 D=D+1.0D0 T=T*A*B*X/(C*D) H=H+P1/(D*B)+EL/(C*A) T1=T*H S=S+T1 IF(ABS(T1)-EPS*ABS(S))13,11,11 11 IF(C-300.0D0)10,12,12 12 IF(IPRINT.EQ.-2)WRITE(6,101) IFAIL=2 13 FDIP0=FDIP0*S RETURN 21 A=EL+1.0D0 B=P-0.5D0 C=EL+P+1.5D0 F=F21(A,B,C,E,EPS,IFAIL) L=L+1 EL=L IF(P-0.5D0)23,23,24 23 C1=EL+EL+1.0D0 GO TO 25 24 C1=1.0D0 25 DO 22 I=1,L AI=I AII=AI+AI C1=C1*AI*AI*4.0D0/(AII*(AII+1.0D0)) 22 CONTINUE FDIP0=FDIP0*F*C1 RETURN 100 FORMAT(' FAILED IN FDIP0, L1=L2=',I5) 101 FORMAT(' FAILED TO CONVERGE IN FDIP0') END C C ****************************************************************** C REAL*8 FUNCTION FDIP1(EK1,L1,EK2,L2) C IMPLICIT REAL*8 (A-H,O-Z) C IF(L1-L2)1,2,3 1 L=L1 A1=EK1 A2=EK2 GO TO 4 2 FDIP1=0.0D0 RETURN 3 L=L2 A1=EK2 A2=EK1 4 LP=L+1 ELP=LP B1=SQRT(1.0D0+ELP*ELP*A2)*FMON1(EK1,EK2,L) B2=SQRT(1.0D0+ELP*ELP*A1)*FMON1(EK1,EK2,LP) IF(B1*B2-1.0D-40)5,5,6 5 FDIP1=0.0D0 RETURN 6 FDIP1=(B1-B2)/ELP RETURN END C C ****************************************************************** C REAL*8 FUNCTION FDIP2(EK1,L1,EK2,L2) C IMPLICIT REAL*8 (A-H,O-Z) C BIG=1.D250 WMAX=200.0D0 ETA1=1.0D0/SQRT(EK1) ETA2=1.0D0/SQRT(EK2) W1=ETA2-ETA1 PI=3.141592653589793D0 A=ABS(W1) B=PI*A IF(B-0.01D0)1,1,2 1 C=3.0D0/(3.0D0-B*(3.0D0-B*(2.0D0-B))) C=SQRT(C) GO TO 5 2 IF(B-14.0D0)4,3,3 3 C=SQRT(B+B) GO TO 5 4 B=B+B C1=1.0D0-EXP(-B) C=SQRT(B/C1) 5 C=0.5D0*C/SQRT(ETA1*ETA2) C2=ETA1+ETA2 C1=4.0D0*ETA1*ETA2/(C2*C2) L=L1 IF(L2-L1)6,6,7 6 L=L2 T1=ETA1 ETA1=ETA2 ETA2=T1 W1=-W1 7 C=C*C1**(L+1) U0=L+1 U1=ETA1 V0=U0 V1=-ETA2 W0=1.0D0 X0=W1/(C2*C2) Y2=-ETA2-ETA2 Y0=-U0*W1+Y2 Y1=ETA2*W1 T1=X0/(1.0D0+W1*W1) Z0=U0*T1 Z1=U1*T1 T=Z0-Z1*W1 Z1=Z0*W1+Z1 Z0=T Q0=-1.0D0+Z0*Y0-Z1*Y1 Q1=Z0*Y1+Z1*Y0 X=W1*X0 8 U0=U0+1.0D0 V0=V0+1.0D0 W0=W0+1.0D0 IF(W0-WMAX)21,21,20 20 FDIP2=0.0D0 RETURN 21 IF(T1-BIG)22,22,20 22 Y0=Y0+Y2 T=Z0*U0-Z1*U1 Z1=Z0*U1+Z1*U0 Z0=T T=Z0*V0-Z1*V1 Z1=Z0*V1+Z1*V0 Z0=T T=Z0*W0-Z1*W1 Z1=Z0*W1+Z1*W0 Z0=T X0=X/(W0*(W0*W0+W1*W1)) Z0=Z0*X0 Z1=Z1*X0 T0=Z0*Y0-Z1*Y1 T1=Z0*Y1+Z1*Y0 Q0=Q0+T0 Q1=Q1+T1 T1=T0*T0+T1*T1 T0=Q0*Q0+Q1*Q1 IF(T0-1.0D+24*T1)8,8,9 9 J1=0 J2=L+1 P=ARGAM(J1,W1)+ARGAM(L,ETA1)-ARGAM(J2,ETA2) IW0=W0 IF(A-1.0D-40)11,11,10 10 P=P+W1*LOG(C2/A) 11 P0=COS(P) P1=SIN(P) T=P0*Q0-P1*Q1 Q1=P0*Q1+P1*Q0 Q0=T FDIP2=C*Q1 RETURN END C C*********************************************************************** C REAL*8 FUNCTION FDIPA(EK1,L1,EK2,L2) IMPLICIT REAL*8 (A-H,O-Z) C C N.R. BADNELL C ASYMPTOTIC EXPRESSION FOR I(KAPPA1,L1,KAPPA2,L2,1) BASED ON A40,1 OF BHT C IF(EK1*EK2.GT.1.0D-50)THEN X1=1.0D0/SQRT(EK1) X2=1.0D0/SQRT(EK2) XP=ABS(X1-X2) IF(XP.GT.1.D2)GO TO 9 PI=ACOS(-1.0D0) XP=EXP(0.5D0*PI*XP) IF(EK1-EK2)1,1,2 1 E=EK1/EK2 IF(L1-L2)3,3,4 3 L=L1 GO TO 7 4 L=L2 GO TO 8 2 E=EK2/EK1 IF(L1-L2)5,5,6 5 L=L1 GO TO 8 6 L=L2 GO TO 7 C A40 7 TL=L T0=1.0D0-E IF(TL*T0.LT.E)GO TO 9 T=PI*TL EE=SQRT(E) F0=SQRT(T*T0*EE)*EE**L TL=L+L+1 FDIPA=F0*XP/TL RETURN C A41 8 T0=1.0D0-E TL=L IF(TL*T0.LT.E)GO TO 9 T0=1.0D0/T0 T=TL*PI EE=SQRT(E) F0=SQRT(T*T0*EE)*EE**(L+1) TL=L+L+1 TL2=L+L+3 FDIPA=F0*XP/(TL*TL2) RETURN ENDIF 9 FDIPA=0.0D0 RETURN END C C*********************************************************************** C REAL*8 FUNCTION FMON1(EK1,EK2,L) C IMPLICIT REAL*8 (A-H,O-Z) C IF(EK1+EK2-1.0D-40)28,28,29 28 FMON1=1.0D+50 RETURN 29 CONTINUE VMAX=200.0D0 X1=SQRT(EK1) X2=SQRT(EK2) X3=X1+X2 X4=X3*X3 X5=X1*X2 X6=X2-X1 X7=4.0D0/X4 PI=3.141592653589793D0 IF(EK1-EK2)1,1,2 1 ETA=1.0D0/X2 GO TO 3 2 ETA=1.0D0/X1 3 G=0.5D0*PI*EXP(-PI*ETA) IF(G.EQ.0.0D0)GO TO 20 !NRB A1=1.0D0 A2=1.0D0 MG=0 MA1=0 MA2=0 M=-1 4 M=M+1 EM=M T=EM+EM+1.0D0 G=G*X7/(T*(T+1.0D0)) EMM=EM*EM A1=A1*(1.0D0+EMM*EK1) A2=A2*(1.0D0+EMM*EK2) 30 IF(G-0.015625D0) 31,32,32 31 G=64.0D0*G MG=MG-1 GO TO 30 32 IF(G-64.0D0) 34,34,33 33 G=0.015625D0*G MG=MG+1 GO TO 32 34 IF(A1-64.0D0) 36,36,35 35 A1=0.015625D0*A1 MA1=MA1+1 GO TO 34 36 IF(A2-64.0D0) 38,38,37 37 A2=0.015625D0*A2 MA2=MA2+1 GO TO 36 38 CONTINUE IF(M-L)4,5,5 5 G=G*(T+1.0D0) IF(X1-300.0D0)7,6,6 6 B=PI/X1 A1=1.5D0*A1/(B*(3.0D0-B*(3.0D0-B*(2.0D0-B)))) GO TO 9 7 IF(X1-0.2D0)9,9,8 8 B=-PI/X1 A1=A1/(1.0D0-EXP(B+B)) 9 IF(X2-300.0D0)11,10,10 10 B=PI/X2 A2=1.5D0*A2/(B*(3.0D0-B*(3.0D0-B*(2.0D0-B)))) GO TO 13 11 IF(X2-0.2D0)13,13,12 12 B=-PI/X2 A2=A2/(1.0D0-EXP(B+B)) 13 G=G*SQRT(A1*A2)*(8.0D0)**(MG+MG+MA1+MA2) S0=1.0D0 S1=0.0D0 U=L V=0.0D0 W=U+U+1.0D0 T0=1.0D0 T1=0.0D0 14 U=U+1.0D0 V=V+1.0D0 W=W+1.0D0 IF(V-VMAX)21,21,20 20 FMON1=0.0D0 RETURN 21 CONTINUE U0=U*U*X5+1.0D0 U1=U*X6 T=T0*U0-T1*U1 T1=T0*U1+T1*U0 T0=T T=X7/(V*W) T0=T*T0 T1=T*T1 S0=S0+T0 S1=S1+T1 S=S0*S0+S1*S1 T=T0*T0+T1*T1 SM=1.0D0/S TM=1.0D0/T IF(SM*TM.EQ.0.0D0)GO TO 20 !NRB IF(S-1.0D+24*T)14,15,15 15 FMON1=G*SQRT(S) IV=V RETURN END C C ****************************************************************** C SUBROUTINE INIT(JMXTWO) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE THE POSSIBLE JP SYMMETRIES FROM INPUT LSP. C C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (HALF=0.5) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK11/NVEC(NMLT),IVEC(NMTM,NMLT),CVEC(NMTM,NMLT),TOLTCC COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/HIGHPW/INOEXCH,JMNTWO COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX C DIMENSION FJJMN(NMPL,2),FJJMX(NMPL,2),NWP(2),FK(2) C C FIRST DETERMINE THE VALUES OF J FOR EACH JPI PARTIAL WAVE C IE=0 IO=0 DO 10 IW=1,NPW IF(LPTT(IW).EQ.0) THEN IE=IE+1 FJJMN(IE,1)=FJTTMN(IW) FJJMX(IE,1)=FJTTMX(IW) ELSE IO=IO+1 FJJMN(IO,2)=FJTTMN(IW) FJJMX(IO,2)=FJTTMX(IW) END IF 10 CONTINUE NWP(1)=IE NWP(2)=IO IF(JMNTWO.GE.0) THEN if(mod(nint(2*fjttmn(1)),2).ne.mod(jmntwo,2))then write(6,*) x 'Warning: Input JMNTWO=',jmntwo,' inconsistent with spins,' x ,'increasing it to JMNTWO=',jmntwo+1 write(0,*) x 'Warning: Input JMNTWO=',jmntwo,' inconsistent with spins,' x ,'increasing it to JMNTWO=',jmntwo+1 jmntwo=jmntwo+1 endif FJTMIN=JMNTWO FJTMIN=FJTMIN/TWO DO 14 IP=1,2 DO 12 I=1,NWP(IP) IF(FJJMN(I,IP).LT.FJTMIN) FJJMN(I,IP)=FJTMIN 12 CONTINUE 14 CONTINUE END IF C C BUBBLE SORT TO PUT FJJMN(I,IP) AND FJJMX(I,IP) IN ORDER OF C INCREASEING FJJMN(I,IP) C DO 40 IP=1,2 KMX=NWP(IP)-1 DO 30 K=1,KMX IMX=NWP(IP)-K DO 20 I=1,IMX IF(FJJMN(I,IP).GT.FJJMN(I+1,IP)) THEN S=FJJMN(I,IP) T=FJJMX(I,IP) FJJMN(I,IP)=FJJMN(I+1,IP) FJJMX(I,IP)=FJJMX(I+1,IP) FJJMN(I+1,IP)=S FJJMX(I+1,IP)=T END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE C C TEST ON MAX TOTAL J C IF(JMXTWO.GE.0)THEN if(mod(nint(2*fjttmn(1)),2).ne.mod(jmxtwo,2))then write(6,*) x 'Warning: Input JMXTWO=',jmxtwo,' inconsistent with spins,' x ,'reducing it to JMXTWO=',jmxtwo-1 write(0,*) x 'Warning: Input JMXTWO=',jmxtwo,' inconsistent with spins,' x ,'reducing it to JMXTWO=',jmxtwo-1 jmxtwo=jmxtwo-1 endif FJMXIN=JMXTWO FJMXIN=FJMXIN*0.5 FJTMAX=MIN(FJTMAX,FJMXIN) ENDIF C C STORE TOTAL J-VALUES AND PARITIES FOR EACH JPI PARTIAL WAVE C K=0 DO 60 IP=1,2 K=K+1 FJ=FJJMN(1,IP) FJTT(K)=FJ JPTT(K)=IP-1 IMX=NWP(IP) DO 50 I=1,IMX IF(FJ.GE.FJJMX(I,IP).OR.FJ.GE.FJTMAX) GO TO 50 IF(FJ.LT.FJJMN(I,IP)) FJ=FJJMN(I,IP)-1.0 45 FJ=FJ+1.0 IF(FJ.GT.FJTMAX) GO TO 50 K=K+1 IF(K.LT.NMPJ) THEN FJTT(K)=FJ JPTT(K)=IP-1 END IF IF(FJ.LT.FJJMX(I,IP)) GO TO 45 50 CONTINUE 60 CONTINUE NPWJ=K !NPWJ NOW DEFINED IF(NPWJ.GT.NMPJ) THEN WRITE(6,550) NMPJ,NPWJ 550 FORMAT(/'**** THE PARAMETER NMPJ MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') STOP 'STOP TO INCREASE THE PARAMETER NMPJ' END IF C C DO BUBBLE SORT TO PUT TOTAL J-VALUES AND PARITIES IN ORDER C OF INCREASING J AND THEN PARITY VALUE C IWMX=NPWJ-1 DO 80 IW=1,IWMX KWMX=NPWJ-IW DO 70 KW=1,KWMX IPJ1=2*FJTT(KW)+JPTT(KW)+0.1 IPJ2=2*FJTT(KW+1)+JPTT(KW+1)+0.1 IF(IPJ1.GT.IPJ2) THEN S=FJTT(KW) IS=JPTT(KW) FJTT(KW)=FJTT(KW+1) JPTT(KW)=JPTT(KW+1) FJTT(KW+1)=S JPTT(KW+1)=IS END IF 70 CONTINUE 80 CONTINUE C C FOR EACH JPI PARTIAL WAVE DETERMINE THE NO OF JK CHANNELS C FLCMN=LCMN FLCMX=LCMX ICMX=0 C C CORRECTION FOR THE NPWJ VARIABLE C NPWJCORR DETERMINES THE NUMBER OF PARTIAL WAVES THAT MUST BE C DELETED IN NPWJ FOR HIGH--J CASES C NPWJCORR=0 DO 150 IWJ=1,NPWJ IF(FJTT(IWJ).GE.HALF) THEN NK=2 FK(1)=FJTT(IWJ)-HALF FK(2)=FJTT(IWJ)+HALF ELSE NK=1 FK(1)=FJTT(IWJ)+HALF END IF IC=0 DO 140 IT=1,NTLS IPD=IABS(JPTT(IWJ)-LPT(IT)) JVMX=NJV(IT) DO 130 JV=1,JVMX DO 120 KV=1,NK FLMN=ABS(FK(KV)-FJT(IT,JV)) IF(FLMN.LT.FLCMN) FLMN=FLCMN FLMX=FK(KV)+FJT(IT,JV) C C IF FLMX EXCEEDS THE MAXIMUM VALUE ALLOWED BY THE INPUT DATA C DELETE THIS PARTIAL WAVE FROM THE FINAL LIST C IF(FLMX.GT.FLCMX) THEN NPWJCORR = NPWJCORR + 1 GO TO 150 ENDIF FL=FLMN-ONE 100 FL=FL+ONE LPD=0 IF(NINT(MOD(FL,TWO)).NE.0) LPD=1 IF(LPD.NE.IPD) GO TO 110 IC=IC+1 110 IF(FL.LT.FLMX) GO TO 100 120 CONTINUE 130 CONTINUE 140 CONTINUE C NCHPJ(IWJ) - - TOTAL NUMBER OF CHANNELS FOR THE IWJ PARTIAL WAVE NCHPJ(IWJ)=IC ICMX=MAX(ICMX,IC) IF(IBUG.NE.0) WRITE(6,599) IWJ,IC 599 FORMAT(/'FOR JK PARTIAL WAVE NO.',I4,' THERE ARE ',I4,' CHANNELS') 150 CONTINUE C IF(ICMX.GT.NMCJ) THEN WRITE(6,551) NMCJ,ICMX 551 FORMAT(/'**** THE PARAMETER NMCJ MUST BE INCREASED FROM ', 1 I4,' TO ',I4,' ****') STOP 'STOP TO INCREASE THE PARAMETER NMCJ' END IF C C CORRECTION OF THE NUMBER OF PARTIAL WAVES C NPWJ = NPWJ - NPWJCORR C C FOR EACH JPI PARTIAL WAVE DETERMINE THE NO OF IC CHANNELS C IF(ITCC.NE.0)THEN FLCMN=LCMN FLCMX=LCMX ICMX=0 C DO 250 IWJ=1,NPWJ IF(FJTT(IWJ).GE.HALF) THEN NK=2 FK(1)=FJTT(IWJ)-HALF FK(2)=FJTT(IWJ)+HALF ELSE NK=1 FK(1)=FJTT(IWJ)+HALF END IF IC=0 DO 240 IL=1,NLEV IPD=IABS(JPTT(IWJ)-JPLEV(IL)) DO 220 KV=1,NK FLEV=JLEV(IL) FLEV=JLEV(IL)/TWO FLMN=ABS(FK(KV)-FLEV) IF(FLMN.LT.FLCMN) FLMN=FLCMN FLMX=FK(KV)+FLEV IF(FLMX.GT.FLCMX) FLMX=FLCMX FL=FLMN-ONE 200 FL=FL+ONE LPD=0 IF(NINT(MOD(FL,TWO)).NE.0) LPD=1 IF(LPD.NE.IPD) GO TO 210 IC=IC+1 210 IF(FL.LT.FLMX) GO TO 200 220 CONTINUE 240 CONTINUE NCHPJX(IWJ)=IC ICMX=MAX(ICMX,IC) IF(IBUG.NE.0) WRITE(6,699) IWJ,IC 699 FORMAT(/'FOR IC PARTIAL WAVE NO.',I4,' THERE ARE ',I4,' CHANNELS') 250 CONTINUE C IF(ICMX.GT.NMCJ) THEN WRITE(6,551) NMCJ,ICMX STOP 'STOP TO INCREASE THE PARAMETER NMCJ' END IF ENDIF C C WRITE DATA TO JBINJK/IC C IF(IJBIN.NE.0) THEN IF(ITCC.EQ.0)THEN C C OPEN FILE KMTJK.DAT AND WRITE HEADER DATA TO FILE JBINJK C OPEN(UNIT=16,FILE='kmtjk.dat',STATUS='UNKNOWN', X FORM="UNFORMATTED") C DO IT=1,NTLS WRITE(15) LT(IT),LST(IT),LPT(IT) END DO WRITE(15) NPWJ WRITE(15) (NCHPJ(IWJ),FJTT(IWJ),JPTT(IWJ),IWJ=1,NPWJ) ELSE C C OPEN FILE KMTIC.DAT AND WRITE HEADER DATA TO FILE JBINIC C OPEN(UNIT=16,FILE='kmtic.dat',STATUS='UNKNOWN', X FORM="UNFORMATTED") C DO IT=1,NTLS WRITE(15) LT(IT),LST(IT),LPT(IT) END DO DO IL=1,NLEV WRITE(15) NVEC(IL) WRITE(15) (IVEC(I,IL),CVEC(I,IL),I=1,NVEC(IL)) END DO WRITE(15) NPWJ WRITE(15) (NCHPJ(IWJ),NCHPJX(IWJ),FJTT(IWJ) X ,JPTT(IWJ),IWJ=1,NPWJ) ENDIF END IF C RETURN END C C ****************************************************************** C SUBROUTINE INITIC(IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C RECOUPLE FROM JK COUPLING TO IC COUPLING C C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (HALF=0.5) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/HIGHPW/INOEXCH,JMNTWO DIMENSION FK(2) C C DETERMINE THE CHANNEL QUANTUM NUMBERS FOR THIS JPI PARTIAL WAVE C FLCMN=LCMN FLCMX=LCMX ILVMX=0 IWJ=IWJ0 C IF(FJTT(IWJ).GE.HALF) THEN NK=2 FK(1)=FJTT(IWJ)-HALF FK(2)=FJTT(IWJ)+HALF ELSE NK=1 FK(1)=FJTT(IWJ)+HALF END IF IC=0 DO 140 IL=1,NLEV IPD=IABS(JPTT(IWJ)-JPLEV(IL)) ILV=0 DO 120 KV=1,NK FLEV=JLEV(IL) FLEV=JLEV(IL)/TWO FLMN=ABS(FK(KV)-FLEV) IF(FLMN.LT.FLCMN) FLMN=FLCMN FLMX=FK(KV)+FLEV IF(FLMX.GT.FLCMX) FLMX=FLCMX FL=FLMN-ONE 100 FL=FL+ONE LPD=0 IF(NINT(MOD(FL,TWO)).NE.0) LPD=1 IF(LPD.NE.IPD) GO TO 110 IC=IC+1 ILVJX(IC)=IL FJVJX(IC)=FLEV FLCJX(IC)=FL FKVJX(IC)=FK(KV) KVSV(IC)=KV ILV=ILV+1 IF(ILV.GT.NMIJ) THEN IF(ILVMX.LT.ILV) ILVMX=ILV ELSE ICHJX(IL,ILV)=IC END IF 110 IF(FL.LT.FLMX) GO TO 100 120 CONTINUE NCHJX(IL)=ILV 140 CONTINUE IF(ILVMX.NE.0) GO TO 200 C IF(IJBIN.NE.0) THEN C C WRITE SOME MORE DATA TO JBINIC FOR THIS PARTIAL WAVE C WRITE(15) (ILVJX(IC),FLCJX(IC),FKVJX(IC) X ,IC=1,NCHPJX(IWJ)) WRITE(15) (NCHJX(IL),IL=1,NLEV) DO IL=1,NLEV IJMX=NCHJX(IL) IF(IJMX.GT.0) THEN WRITE(15) (ICHJX(IL,IJ),IJ=1,IJMX) END IF END DO END IF 200 CONTINUE C IF(ILVMX.NE.0) THEN WRITE(6,312) NMIJ,ILVMX 312 FORMAT(/'**** THE PARAMETER NMIJ MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') STOP 'STOP TO INCREASE THE PARAMETER NMIJ' END IF C RETURN END C C ****************************************************************** C SUBROUTINE INITJK(IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C RECOUPLE TO JK COUPLING FROM LS COUPLING C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (HALF=0.5) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/HIGHPW/INOEXCH,JMNTWO DIMENSION FK(2) C SAVE MXNCT DATA MXNCT/0/ C C DETERMINE THE CHANNEL QUANTUM NUMBERS FOR EACH JPI PARTIAL WAVE IWJ C IWJ=IWJ0 FLCMN=LCMN FLCMX=LCMX ICMX=0 ILVMX=0 C IF(FJTT(IWJ).GE.HALF) THEN NK=2 FK(1)=FJTT(IWJ)-HALF FK(2)=FJTT(IWJ)+HALF ELSE NK=1 FK(1)=FJTT(IWJ)+HALF END IF IC=0 ILL=0 DO 140 IT=1,NTLS IPD=IABS(JPTT(IWJ)-LPT(IT)) JVMX=NJV(IT) DO 130 JV=1,JVMX ILV=0 ILL=ILL+1 DO 120 KV=1,NK FLMN=ABS(FK(KV)-FJT(IT,JV)) IF(FLMN.LT.FLCMN) FLMN=FLCMN FLMX=FK(KV)+FJT(IT,JV) FL=FLMN-ONE 100 FL=FL+ONE LPD=0 IF(NINT(MOD(FL,TWO)).NE.0) LPD=1 IF(LPD.NE.IPD) GO TO 110 IC=IC+1 C C STORE REQUIRED INFORMATION ABOUT THE K-MATRIX IN JK COUPLING: C C ITMJ(IC) - - TERM NUMBER FOR THE PARTIAL WAVE IWJ AND C CHANNEL IC C ILVJ(IC) - - LEVEL NUMBER FOR THE PARTIAL WAVE IWJ AND C CHANNEL IC C FJVJ(IC) - - J VALUE OF LEVEL FOR THE PARTIAL WAVE IWJ AND C CHANNEL IC C FLCJ(IC) - - L VALUE OF THE CONTINUUM ELECTRON FOR THE C PARTIAL WAVE IWJ AND THE CHANNEL IC C FKVJ(IC) - - K VALUE OF THE PARTIAL WAVE IWJ AND CHANNEL IC. C ICHN(IT,JV,INL,KV) -- CHANNEL NUMBER FOR THE IWJ PARTIAL WAVE, C FOR A GIVEN LEVEL JV OF TERM IT, FOR A C GIVEN VALUE OF INL (L+1 FOR THE C CONTINUUM ELECTRON) AND A GIVEN C K-VALUE. C ICHJ(IT,JV,ILV) - - CHANNEL NUMBER FOR THE IWJ PARTIAL WAVE, C FOR TERM NUMBER IT, FOR LEVEL JV OF TERM C IT; ILV IS JUST AN INDEX FOR EACH CHANNEL C WITH THIS SET OF IWJ, IT, AND JV, AND C NCHJ(IT,JV) - - IS THE TOTAL NUMBER OF CHANNELS WITH THESE C VALUES OF IWJ, IT, AND JV. C ITMJ(IC)=IT ILVJ(IC)=ILL FJVJ(IC)=FJT(IT,JV) FLCJ(IC)=FL FKVJ(IC)=FK(KV) ILV=ILV+1 IF (ILV.GT.NMIJ)THEN IF(ILVMX.LT.ILV)ILVMX=ILV ELSE ICHJ(IT,JV,ILV)=IC ENDIF INL=FL+1.1 ICHN(IT,JV,INL,KV)=IC 110 IF(FL.LT.FLMX) GO TO 100 120 CONTINUE NCHJ(IT,JV)=ILV 130 CONTINUE 140 CONTINUE C IF(ILVMX.GT.NMIJ)THEN WRITE(6,249)NMIJ,ILVMX 249 FORMAT(/'**** THE PARAMETER NMIJ MUST BE INCREASED FROM', 1 I3,' TO ',I3,' ****') STOP 'STOP TO INCREASE THE PARAMETER NMIJ' ENDIF C IF(ITCC.EQ.0.AND.IJBIN.NE.0) THEN C C WRITE DATA FOR THIS PARTIAL WAVE TO FILE JBINJK C WRITE(15) (ITMJ(ICJ),FLCJ(ICJ),FKVJ(ICJ) X ,ILVJ(ICJ),ICJ=1,NCHPJ(IWJ)) DO IT=1,NTLS WRITE(15) (NCHJ(IT,IJ),IJ=1,NJV(IT)) DO IJ=1,NJV(IT) JIMX=NCHJ(IT,IJ) IF(JIMX.NE.0) THEN WRITE(15) (ICHJ(IT,IJ,JI),JI=1,JIMX) END IF END DO END DO END IF C C NOW GENERATE JK COUPLING COEFFICIENT FOR THIS JPI PARTIAL WAVE C NPWFMX=0 ICJMX=NCHPJ(IWJ) NPWF=0 DO IW=1,NPW IF(FJTT(IWJ).GE.FJTTMN(IW).AND.FJTT(IWJ).LE.FJTTMX(IW).AND. 1 JPTT(IWJ).EQ.LPTT(IW)) THEN NPWF=NPWF+1 IF(NPWF.GT.NMLS) THEN IF(NPWFMX.LT.NPWF) NPWFMX=NPWF ELSE IWF(NPWF)=IW END IF END IF ENDDO IF (NPWFMX.GT.0) GO TO 280 NPWFSV(IWJ)=NPWF C NCT=0 DO IL=1,NPWF IW=IWF(IL) NCT=NCT+(NCHT(IW)*(NCHT(IW)+1))/2 DO ICJ=1,ICJMX IT=ITMJ(ICJ) ICFJ(ICJ,IL)=0 CLS(ICJ,IL)=TZERO IIMX=NCH(IW,IT) DO II=1,IIMX IC=ICH(IW,IT,II) IF(NINT(ABS(FLCJ(ICJ)-FLC(IW,IC))).EQ.0) THEN COEFF=ONE IF(NINT(MOD((FJTT(IWJ)+HALF-FLCJ(ICJ)- 1 FJVJ(ICJ)),TWO)).NE.0)COEFF=-ONE COEFF=COEFF*SQRT((TWO*FLTT(IW)+ONE)*(TWO*FSTT(IW)+ONE)* 1 (TWO*FJVJ(ICJ)+ONE)*(TWO*FKVJ(ICJ)+ONE)) S6J1=S6J(FLTT(IW),FST(IT),FKVJ(ICJ),HALF,FJTT(IWJ), 1 FSTT(IW)) S6J2=S6J(FLCJ(ICJ),FLT(IT),FLTT(IW),FST(IT),FKVJ(ICJ), 1 FJVJ(ICJ)) CLS(ICJ,IL)=COEFF*S6J1*S6J2 ICFJ(ICJ,IL)=IC END IF ENDDO ENDDO ENDDO C MXNCT=MAX(NCT,MXNCT) IF(IWJ.EQ.NPWJ)THEN WRITE(6,254)MXRLS,MXNCT 254 FORMAT(/' MXRLS=',I8,' NEEDED =',I8/) ENDIF C IF(NCT.GT.MXRLS)THEN WRITE(6,253) MXRLS,NCT 253 FORMAT(/'**** THE PARAMETER MXRLS MUST BE INCREASED FROM ', 1 I6,' TO ',I6,' ****') STOP 'STOP TO INCREASE THE PARAMETER MXRLS' ENDIF C 280 IF(NPWFMX.GT.0) THEN WRITE(6,252) NMLS,NPWFMX 252 FORMAT(/'**** THE PARAMETER NMLS MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') STOP 'STOP TO INCREASE THE PARAMETER NMLS' END IF C RETURN END C C ****************************************************************** C SUBROUTINE INTERIC(IEE,EINCR,IWJ0,IEK) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO CALCULATE THE UNPHYSICAL K-MATRICES ON A NEW ENERGY C MESH FROM LINEAR INTERPOLATION OF THE UNPHYSICAL K-MATRICES. C C INCLUDE 'PARAM' C PARAMETER (TWO=2.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/INTER/RJ1(NOMCJ),RJ2(NOMCJ),IEK0 SAVE KLAST C C DETERMINE THE NUMBER OF OPEN CHANNELS C IWJ=IWJ0 NCHOJ(IWJ)=0 DO ICJ=1,NCHPJX(IWJ) IL=ILVJX(ICJ) IF(AZSQ*EMSH(IEE).GE.ELEV(IL)) NCHOJ(IWJ)=NCHOJ(IWJ)+1 END DO ICJMX=(NCHPJX(IWJ)*(NCHPJX(IWJ)+1))/2 C C IMODE.GT.0 SEE IF WE NEED A NEW ENERGY YET C IF(IMODE.GT.0)THEN IF(IEE.EQ.1)THEN !ADVANCE TO FIRST ENERGY DO N=1,IEK0-1 READ(18) ENDDO KLAST=0 READ(18) (RJ1(N),N=1,ICJMX) ENDIF IF(IEK.GT.KLAST)THEN !READ NEW ENERGY KLAST=IEK IF(MOD(IEK-IEK0+1,2).EQ.0)THEN READ(18) (RJ2(N),N=1,ICJMX) ELSE READ(18) (RJ1(N),N=1,ICJMX) ENDIF ENDIF IF(IEE.EQ.MXE)THEN !ADVANCE TO START OF NEW SYMMETRY BLOCK DO N=1,NUMK-IEK READ(18) ENDDO ENDIF ENDIF C C LINEAR INTERPOLATION TO DETERMINE THE UNPHYSICAL K-MATRICES C IF(MOD(IEK-IEK0+1,2).EQ.0)THEN DO N=1,ICJMX SLOPE=(RJ2(N)-RJ1(N))/(EMESH(IEK)-EMESH(IEK-1)) RJ(N)=RJ1(N)+SLOPE*(EMSH(IEE)-EMESH(IEK-1)) END DO ELSE DO N=1,ICJMX SLOPE=(RJ1(N)-RJ2(N))/(EMESH(IEK)-EMESH(IEK-1)) RJ(N)=RJ2(N)+SLOPE*(EMSH(IEE)-EMESH(IEK-1)) END DO ENDIF C RETURN END C C ****************************************************************** C SUBROUTINE INTERJK(IEE,EINCR,IWJ0,IEK) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO CALCULATE THE UNPHYSICAL K-MATRICES ON A NEW ENERGY C MESH FROM LINEAR INTERPOLATION OF THE UNPHYSICAL K-MATRICES. C C INCLUDE 'PARAM' C PARAMETER (TWO=2.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/INTER/RJ1(NOMCJ),RJ2(NOMCJ),IEK0 SAVE KLAST C C DETERMINE THE NUMBER OF OPEN CHANNELS C IWJ=IWJ0 NCHOJ(IWJ)=0 DO ICJ=1,NCHPJ(IWJ) ITJ=ITMJ(ICJ) IF(AZSQ*EMSH(IEE).GE.ENT(ITJ)) NCHOJ(IWJ)=NCHOJ(IWJ)+1 END DO ICJMX=(NCHPJ(IWJ)*(NCHPJ(IWJ)+1))/2 C C IMODE.GT.0 SEE IF WE NEED A NEW ENERGY YET C IF(IMODE.GT.0)THEN IF(IEE.EQ.1)THEN !ADVANCE TO FIRST ENERGY DO N=1,IEK0-1 READ(18) ENDDO KLAST=0 READ(18) (RJ1(N),N=1,ICJMX) ENDIF IF(IEK.GT.KLAST)THEN !READ NEW ENERGY KLAST=IEK IF(MOD(IEK-IEK0+1,2).EQ.0)THEN READ(18) (RJ2(N),N=1,ICJMX) ELSE READ(18) (RJ1(N),N=1,ICJMX) ENDIF ENDIF IF(IEE.EQ.MXE)THEN !ADVANCE TO START OF NEW SYMMETRY BLOCK DO N=1,NUMK-IEK READ(18) ENDDO ENDIF ENDIF C C LINEAR INTERPOLATION TO DETERMINE THE UNPHYSICAL K-MATRICES C IF(MOD(IEK-IEK0+1,2).EQ.0)THEN DO N=1,ICJMX SLOPE=(RJ2(N)-RJ1(N))/(EMESH(IEK)-EMESH(IEK-1)) RJ(N)=RJ1(N)+SLOPE*(EMSH(IEE)-EMESH(IEK-1)) END DO ELSE DO N=1,ICJMX SLOPE=(RJ1(N)-RJ2(N))/(EMESH(IEK)-EMESH(IEK-1)) RJ(N)=RJ2(N)+SLOPE*(EMSH(IEE)-EMESH(IEK-1)) END DO ENDIF C RETURN END C C ****************************************************************** C SUBROUTINE KMTRIC(IEK,IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE THE K-MATRICES IN IC COUPLING FROM THE K-MATRICES IN JK C COUPLING C C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (HALF=0.5) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK11/NVEC(NMLT),IVEC(NMTM,NMLT),CVEC(NMTM,NMLT),TOLTCC COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/INTER/RJ1(NOMCJ),RJ2(NOMCJ),IEK0 COMMON/HIGHPW/INOEXCH,JMNTWO COMMON/NRBWRK/WORK(NMCJ) C IPOS(I,J,N)=J+N*(I-1)-(I*(I-1))/2 C C GENERATE K-MATRIX IN IC COUPLING FOR JPI PARTIAL WAVE IWJ0 C N=0 IWJ=IWJ0 NCHOJ(IWJ)=0 ICJMX=NCHPJX(IWJ) DO IC=1,ICJMX IL=ILVJX(IC) IF(AZSQ*EMESH(IEK).GE.ELEV(IL)) NCHOJ(IWJ)=NCHOJ(IWJ)+1 INJ=NINT(TWO*FJVJX(IC)+ONE) INL=NINT(FLCJX(IC)+ONE) KV=KVSV(IC) ICJM=NCHPJ(IWJ) C DO ICCP=1,ICJM WORK(ICCP)=TZERO ENDDO DO IVK=1,NVEC(IL) IT=IVEC(IVK,IL) JV=JVSV(IT,INJ) IF(JV.GT.0)THEN ICC=ICHN(IT,JV,INL,KV) CVP=CVEC(IVK,IL) NN=(ICC*(ICC-1))/2 DO ICCP=1,ICC NN=NN+1 WORK(ICCP)=WORK(ICCP)+CVP*A(NN) ENDDO C NB ICCP CAN BE .GT. ICC NN=ICJM*(ICC-1)-NN+2*ICC DO ICCP=ICC+1,ICJM NN=NN+1 WORK(ICCP)=WORK(ICCP)+CVP*RJ(NN) ENDDO ENDIF ENDDO C DO ICP=IC,ICJMX !LOWER HALF N=N+1 ILP=ILVJX(ICP) INJP=NINT(TWO*FJVJX(ICP)+ONE) INLP=NINT(FLCJX(ICP)+ONE) KVP=KVSV(ICP) TM(N)=TZERO DO IVKP=1,NVEC(ILP) ITP=IVEC(IVKP,ILP) JVP=JVSV(ITP,INJP) IF(JVP.GT.0) THEN ICCP=ICHN(ITP,JVP,INLP,KVP) TM(N)=TM(N)+WORK(ICCP)*CVEC(IVKP,ILP) ENDIF ENDDO ENDDO ENDDO C C NOW STORE THE INTERMEDIATE COUPLING UNPHYSICAL K-MATRIX IN C RJ FOR THIS PARTIAL WAVE. C NMAX=N DO N=1,NMAX RJ(N)=TM(N) END DO C C WRITE THE UNPHYSICAL K-MATRIX IN IC COUPLING TO KMTIC C IF(IPRTIC.NE.0) THEN WRITE(12,300) FJTT(IWJ),JPTT(IWJ),ICJMX,NCHOJ(IWJ) 300 FORMAT(//48X,'K-MATRIX FOR J =',F5.1,' PARITY =',I3/ 1 43X,I4,' TOTAL CHANNELS AND ',I4,' OPEN CHANNELS'/ 2 48X,35('*')) DO 310 IC=1,ICJMX WRITE(12,311) IC,ILVJX(IC),FJVJX(IC),FLCJX(IC),FKVJX(IC) 1 ,(RJ(IPOS(MIN(IC,ICP),MAX(IC,ICP),ICJMX)),ICP=1,ICJMX) 311 FORMAT(1X,I2,I3,3F6.1,1PE12.4,6E13.4/(24X,E12.4,6E13.4)) 310 CONTINUE END IF C C WRITE THE UNPHYSICAL K-MATRIX TO FILE KMTIC.DAT C CASE IJBIN=IMODE.EQ.0.AND.JMNTWO.LT.0 C IF(IJBIN.NE.0)WRITE(16)(RJ(N),N=1,NMAX) C C STORE IF IMODE.LT.0. C IF(IMODE.LT.0)THEN IF(MOD(IEK-IEK0+1,2).EQ.1)THEN DO N=1,NMAX RJ1(N)=RJ(N) ENDDO ELSE N=0 DO N=1,NMAX RJ2(N)=RJ(N) ENDDO ENDIF ENDIF C RETURN END C C ****************************************************************** C SUBROUTINE KMTRJK(IEK,IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE THE UNPHYSICAL K-MATRICES IN JK COUPLING FROM THE C THE UNPHYSICAL K-MATRICES IN LS COUPLING C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (HALF=0.5) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/HIGHPW/INOEXCH,JMNTWO COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/INTER/RJ1(NOMCJ),RJ2(NOMCJ),IEK0 C NPOS(ICMN,ICMX)=(ICMX*(ICMX-1))/2+ICMN IPOS(I,J,N)=J+N*(I-1)-(I*(I-1))/2 C C GENERATE K-MATRIX IN JK COUPLING FOR JPI PARTIAL WAVE IWJ0 C IWJ=IWJ0 ICJM=NCHPJ(IWJ) C C DETERMINE THE NUMBER OF OPEN CHANNELS IN JK COUPLING FOR THIS C PARTIAL WAVE C NCHOJ(IWJ)=0 DO ICJ=1,ICJM ITJ=ITMJ(ICJ) IF(AZSQ*EMESH(IEK).GE.ENT(ITJ)) NCHOJ(IWJ)=NCHOJ(IWJ)+1 END DO C NMAX=(ICJM*(ICJM+1))/2 DO N=1,NMAX RJ(N)=TZERO ENDDO IF(ITCC.NE.0)THEN DO N=1,NMAX A(N)=TZERO ENDDO ENDIF C DO IL=1,NPWFSV(IWJ) N=0 DO ICJ=1,ICJM IC=ICFJ(ICJ,IL) DO ICJP=ICJ,ICJM !LOWER HALF N=N+1 ICP=ICFJ(ICJP,IL) IF(IC.NE.0.AND.ICP.NE.0) THEN C NB ICP CAN BE .GT. IC, SO ICMN=MIN(IC,ICP) ICMX=MAX(IC,ICP) NLS=NPOS(ICMN,ICMX) !LS UPPER HALF FROM STGF... NLS=NLS+LSPOS(IL-1) RJ(N)=RJ(N)+CLS(ICJ,IL)*CLS(ICJP,IL)*RLS(NLS) END IF ENDDO ENDDO ENDDO C IF(ITCC.NE.0)THEN N=0 DO J=1,ICJM DO I=J,ICJM N=N+1 NN=NPOS(J,I) A(NN)=RJ(N) ENDDO ENDDO ENDIF C C WRITE THE UNPHYSICAL K-MATRIX IN JK COUPLING TO KMTJK C IF(IPRTJK.NE.0) THEN WRITE(11,300) FJTT(IWJ),JPTT(IWJ),ICJM,NCHOJ(IWJ) 300 FORMAT(//48X,'K-MATRIX FOR J =',F5.1,' PARITY =',I3,/ 1 43X,I4,' TOTAL CHANNELS AND ',I4,' OPEN CHANNELS'/ 2 43X,44('*')) DO 240 IC=1,ICJM WRITE(11,310) IC,ITMJ(IC),FJVJ(IC),FLCJ(IC),FKVJ(IC), 1 (RJ(IPOS(MIN(IC,ICP),MAX(IC,ICP),ICJM)),ICP=1,ICJM) 310 FORMAT(1X,I2,I3,3F6.1,1PE12.4,6E13.4/(24X,E12.4,6E13.4)) 240 CONTINUE END IF C C IF ITCC.EQ.0 WRITE THE UNPHYSICAL K-MATRIX TO FILE KMTJK.DAT C CASE IJBIN=IMODE.EQ.0.AND.JMNTWO.LT.0 C IF(ITCC.EQ.0.AND.IJBIN.NE.0)WRITE(16)(RJ(N),N=1,NMAX) C C STORE IF ITCC.EQ.0 AND IMODE.LT.0. C IF(ITCC.EQ.0.AND.IMODE.LT.0)THEN IF(MOD(IEK-IEK0+1,2).EQ.1)THEN DO N=1,NMAX RJ1(N)=RJ(N) ENDDO ELSE N=0 DO N=1,NMAX RJ2(N)=RJ(N) ENDDO ENDIF ENDIF C RETURN END C****************************************************************** C SUBROUTINE LUBS(A,B,LB,NB,IERR) C ________________________________________________________ C | | C | SOLVE AN LU FACTORED SYMMETRIC SYSTEM WITHOUT PIVOTING | C | | C | INPUT: | C | | C | A --LUS'S OUTPUT | C | | C | B --RIGHT SIDE (DESTROYED) | C | | C | LB --LEADING (ROW) DIMENSION OF ARRAY B | C | | C | NB --DIMENSION OF MATRIX STORED IN B | C | | C | OUTPUT: | C | | C | B --SOLUTION | C |________________________________________________________| C REAL*8 A(*),B(LB,*),T,TZERO,ONE INTEGER I,J,K,L,N,LB,NB,IB,I1,IERR,KB,NB0 C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C IERR=0 NB0=ABS(NB) C I1 = NINT(A(1)) IF ( ABS(I1) .NE. 1233 ) THEN IERR=1 !ERROR, MUST FACTOR BEFORE SOLVING RETURN ENDIF C ----------------------------- C |*** FORWARD ELIMINATION ***| C ----------------------------- DO 20 IB = 1,NB0 KB=1 IF(NB.LT.0)KB=IB N = NINT(A(2)) L = 3 K=1 IF ( I1 .LT. 0 ) GOTO 80 30 IF ( K .EQ. N ) GOTO 50 T = B(K,IB)/A(K+L) J = L L = L + N - K K = K + 1 IF ( T .EQ. TZERO ) GOTO 30 DO 40 I = K,N 40 B(I,IB) = B(I,IB) - T*A(I+J) GOTO 30 C ----------------------------------- C |*** BACK SUBSTITUTION BY ROWS ***| C ----------------------------------- 50 B(N,IB) = B(N,IB)/A(K+L) 60 IF ( K .EQ. KB ) GO TO 20 J = K K = K - 1 L = L + K - N T = B(K,IB) DO 70 I = J,N 70 T = T - B(I,IB)*A(I+L) B(K,IB) = T/A(K+L) GOTO 60 C ----------------------------- C |*** COMPUTE NULL VECTOR ***| C ----------------------------- 80 IF ( A(K+L) .EQ. TZERO ) GOTO 90 L = L + N - K K = K + 1 GOTO 80 90 DO 100 I = 1,N 100 B(I,IB) = TZERO B(K,IB) = ONE GOTO 60 20 CONTINUE IF(NB.LT.0)THEN DO J=1,IB DO I=J,IB B(J,I)=B(I,J) ENDDO ENDDO ENDIF RETURN END C****************************************************************** C SUBROUTINE LUS(A,LA,N,W,IERR) C C ________________________________________________________ C | | C | LU FACTOR A SYMMETRIC MATRIX WITHOUT PIVOTING | C | | C | INPUT: | C | | C | A --ARRAY CONTAINING MATRIX | C | (ONLY THE LOWER HALF NEED BE DEFINED) | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | SET .LE. 0 IF A ALREADY PACKED | C | | C | N --MATRIX DIMENSION | C | | C | W --WORK ARRAY WITH LENGTH AT LEAST N | C | | C | OUTPUT: | C | | C | A --INVERSE (IN LOWER HALF ONLY) | C |________________________________________________________| C REAL*8 A(*),W(*),R,S,T INTEGER G,H,I,J,K,L,M,N REAL*8 TZERO C PARAMETER (TZERO=0.0) C IERR=0 C IF(LA.LE.0)GO TO 4 !ALREADY PACKED C ---------------- CNRB |*** PACK A ***| C ---------------- H=LA-N I=0 M=0 L=N G=(N*(N+1))/2 2 IF(L.EQ.G)GO TO 4 K=L+1 M=M+1 L=L+N-M I=I+H+M DO J=K,L A(J)=A(I+J) ENDDO GO TO 2 C C ------------------------ C |*** COMPUTE 1-NORM ***| C ------------------------ 4 DO 10 I = 1,N 10 W(I) = TZERO I = -N K = 0 R = TZERO S = TZERO 20 I = I + N - K K = K + 1 J = K S = ABS(A(I+J)) 30 IF ( J .EQ. N ) GOTO 40 J = J + 1 T = ABS(A(I+J)) S = S + T W(J) = W(J) + T GOTO 30 40 S = S + W(K) IF ( R .LT. S ) R = S IF ( K .LT. N ) GOTO 20 J = 3 + (N+N*N)/2 C ----------------------------------- C |*** SHIFT MATRIX DOWN 3 SLOTS ***| C ----------------------------------- 50 A(J) = A(J-3) J = J - 1 IF ( J .GT. 3 ) GOTO 50 A(1) = 1233 A(2) = N A(3) = R H = N K = 4 60 IF ( H .EQ. 1 ) GOTO 90 C -------------------------- C |*** SAVE PIVOT ENTRY ***| C -------------------------- S = A(K) K = K + H G = K H = H - 1 M = H IF ( S .EQ. TZERO ) GOTO 100 J = 0 70 J = J - M M = M - 1 L = G + M T = A(G+J)/S C --------------------------- C |*** ELIMINATE BY ROWS ***| C --------------------------- DO 80 I = G,L 80 A(I) = A(I) - T*A(I+J) G = L + 1 IF ( M .GT. 0 ) GOTO 70 GOTO 60 90 IF ( A(K) .NE. TZERO ) RETURN A(1) = -1233 RETURN 100 A(1) = -1233 GOTO 60 END C C ****************************************************************** C SUBROUTINE MQDTK(IEE,IWJ) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE THE PHYSICAL K MATRIX FROM THE UNPHYSICAL K-MATRIX C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (THOU=1000.) PARAMETER (FNUMAX=THOU*THOU) parameter (fnumin=0.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/NRBWRK/WORK(NMCJ) C CBL ALLOCATABLE :: TEMP(:) !BLAS (FOR UNSLICED DGEMM) C NCHOP=NCHOJ(IWJ) C PI=ACOS(-ONE) IF(ITCC.EQ.0)NCC=NCHPJ(IWJ)-NCHOP IF(ITCC.NE.0)NCC=NCHPJX(IWJ)-NCHOP C IF(NCC.LT.1) RETURN C C STORE KCC FIRST C NCMAX=(NCC*(NCC+1))/2 NOMAX=(NCHOP*(NCHOP+1))/2 NN=NOMAX+NCC*NCHOP DO N=1,NCMAX NN=NN+1 A(N)=RJ(NN) !PACK LOWER HALF END DO C C ADD TAN(NU*PI) DOWN THE DIAGONAL C N=1 ICT=NCHOP DO ICC=NCC,1,-1 ICT=ICT+1 if(fnu(ict).gt.fnumin)then A(N)=A(N)+TAN(FNU(ICT)*PI) else do j=n+1,n+icc-1 a(j)=tzero enddo a(n)=one endif N=N+ICC END DO C C INITIALIZE KCO C N=0 NN=0 DO IC=NCHOP,1,-1 NN=NN+IC DO ICC=1,NCC N=N+1 NN=NN+1 ICT=ICC+NCHOP IF(FNU(ICT).LE.FNUMAX.and.fnu(ict).gt.fnumin) THEN TM(N)=RJ(NN) ELSE TM(N)=TZERO ENDIF ENDDO ENDDO C CSTRTNBL CALL LUS(A,-NMCJ,NCC,WORK,IERR) IF (IERR.NE.0) THEN WRITE(6,600) STOP 'ERROR IN LUS' END IF CALL LUBS(A,TM,NCC,NCHOP,IERR) IF (IERR.NE.0) THEN WRITE(6,601) STOP 'ERROR IN LUBS' END IF CENDNBL C CSTRTBL CBL CALL DSPTRF('L',NCC,A,IPIVOT,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,602) INFO CBL STOP 'FAILURE IN BLAS ROUTINE DSPTRF' CBL ENDIF CBL CALL DSPTRS('L',NCC,NCHOP,A,IPIVOT,TM,NCC,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,603) INFO CBL STOP 'FAILURE IN BLAS ROUTINE DSPTRS' CBL ENDIF CENDBL C C C INITIALIZE KOC^T=KCO C N=0 NN=0 DO IC=NCHOP,1,-1 NN=NN+IC DO ICC=1,NCC N=N+1 NN=NN+1 ICT=ICC+NCHOP IF(FNU(ICT).LE.FNUMAX.and.fnu(ict).gt.fnumin) THEN A(N)=RJ(NN) ELSE A(N)=TZERO ENDIF ENDDO ENDDO C C PACK KOO C N=NCHOP NN=NCHOP DO ICP=2,NCHOP NN=NN+NCC DO IC=ICP,NCHOP N=N+1 NN=NN+1 RJ(N)=RJ(NN) ENDDO ENDDO C C NOW FINISH UP THE CALCULATION OF THE PHYSICAL K-MATRIX C CSTRTNBL N=0 NN=0 DO ICP=1,NCHOP NI=NN DO IC=ICP,NCHOP N=N+1 NJ=NN DO ICC=1,NCC NI=NI+1 NJ=NJ+1 RJ(N)=RJ(N)-A(NI)*TM(NJ) ENDDO ENDDO NN=NN+NCC ENDDO CENDNBL C CSTRTBL CBL ALLOCATE (TEMP(NCHOP*NCHOP),stat=ierr) CBL if(ierr.ne.0)stop 'SR.MQDTK: Failure to allocate temp' CBL CALL DGEMM('T','N',NCHOP,NCHOP,NCC,ONE,A,NCC,TM,NCC CBL X ,TZERO,TEMP,NCHOP) CBL N=0 CBL NN=0 CBL DO ICP=1,NCHOP CBL DO IC=ICP,NCHOP CBL N=N+1 CBL NN=NN+1 CBL RJ(N)=RJ(N)-TEMP(NN) CBL ENDDO CBL NN=NN+ICP CBL ENDDO CBL DEALLOCATE (TEMP,stat=ierr) CBL if(ierr.ne.0)stop 'SR.MQDTK: Failure to deallocate temp' CENDBL C RETURN C 600 FORMAT(' SR.MQDTK: LUS RETURNED WITH INFO =',I6) 601 FORMAT(' SR.MQDTK: LUBS RETURNED WITH INFO =',I6) CB602 FORMAT(//10X,10('*'),' SR.MQDTK: DSPTRF RETURNED WITH INFO =',I6) CB603 FORMAT(//10X,10('*'),' SR.MQDTK: DSPTRS RETURNED WITH INFO =',I6) END C****************************************************************** C SUBROUTINE MULTS(Y,X,A,N) C ________________________________________________________ C | | C | SYMMETRIC MATRIX MULTIPLICATION: Y=AX | C | | C |NRB: ADAPTED FROM SYMMETRIC MATRIX TIMES VECTOR: MULTSV | C | | C | INPUT: | C | | C | A,X --ARRAYS PACKED WITH ELEMENTS CONTAINED | C | IN EACH ROW, ON DIAGONAL AND TO RIGHT, | C | OF COEFFICIENT MATRIX | C | | C | N --MATRIX DIMENSION | C | | C | OUTPUT: | C | | C | Y --PRODUCT BETWEEN A AND X (PACKED) | C |________________________________________________________| C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(*),X(*),Y(*) C PARAMETER (TZERO=0.0) C NN=N*(N+1)/2 DO I=1,NN Y(I)=TZERO ENDDO C LM=0 DO M=1,N C K = 1 L = 0 KM=1 IF ( N .EQ. 1 ) GOTO 40 DO J = 2,M T = X(M+L) IM=LM DO I = M,N IM=IM+1 Y(IM) = Y(IM) + T*A(I+L) ENDDO L = L + N - K K = J ENDDO LM=LM+N-M+1 KM=K+L DO J = M+1,N T = X(KM) S = A(K+L)*T IM=KM DO I = J,N IM=IM+1 R = A(I+L) S = S + R*X(IM) Y(IM) = Y(IM) + R*T ENDDO Y(KM) = Y(KM) + S L = L + N - K K = J KM=KM+1 ENDDO 40 Y(KM) = Y(KM) + A(K+L)*X(KM) ENDDO C RETURN END C*************************************************************** C SUBROUTINE OMEG(IE,EMESH,MXE,ENAT,NAST,OMEGA,NOMWRT,IONE) C C USE MEMORY FOR SUMMING OMEGA, WITH OVERFLOW ONTO DISK C (ADAPTED BY NRB, FROM KAB, FOR ELASTIC TRANSITIONS C AND USE OF NOMT=ABS(NOMWRT) TO REDUCE I/O) C IE = 0 TO INITIALIZE C IE .GT.0 TO RETRIEVE OMEGA FOR ENERGY IE C IF NOMWRT .LT. 0 THEN RECOVER -NOMWRT NON-ZERO OMEGAS FOR C FINAL WRITE RUNNING DOWN THE COLUMNS C IF NOMWRT .GT. 0 THEN RECOVER NOMWRT OMEGAS (INC ZEROES) FOR C FINAL WRITE RUNNING ALONG THE ROWS C (EXC ZEROES CAN BE IMPLEMENTED BY SWITCHING JJJJ BELOW) C EMESH(MXE) = ENERGY MESH C ENAT(NAST) = EXCITATION THRESHOLDS: NEEDED TO MINIMISE STORAGE C OMEGA(NOMT)= RETURNS OMEGA MATRIX AS AN ARRAY C /MEMORY/ STORAGE IN OMEM(MXOM) C MPOS(IE) = LAST POSITION OCCUPIED IN OMEM AT ENERGY IE, OR C IF MPOS.LT.0 = -RECORD NUMBER ON DA SCRATCH FILE C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) C COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX C DIMENSION EMESH(0:MXE),ENAT(NAST),OMEGA(NOMLT) COMMON/WORK1/OMST(NOMLT) C NOMT=ABS(NOMWRT) N=0 IF(NOMWRT.GT.0)THEN DO I=1,NAST DO J=I+IONE,NAST N=N+1 OMEGA(N)=TZERO ENDDO IF(N.GE.NOMT)GO TO 20 ENDDO I=NAST ENDIF IF(NOMWRT.LT.0)THEN DO J=1+IONE,NAST DO I=1,J-IONE N=N+1 OMEGA(N)=TZERO ENDDO IF(N.GE.NOMT)GO TO 20 ENDDO J=NAST ENDIF C C INITIALISE C 20 IF(IE.EQ.0) THEN IF(NOMWRT.GT.0)THEN ITMAX=I JTMAX=NAST ELSE ITMAX=NAST JTMAX=J ENDIF MPOS(0)=0 NPOS=0 NREC=0 NST=1 DO I=1,MXE 4 IF(NST.LE.NAST)THEN IF(EMESH(I).GE.ENAT(NST)) THEN NP=MIN(NST,ITMAX,JTMAX) NTRAN=(NP*(NP+1-2*IONE))/2 IF(NST.GT.ITMAX)NTRAN=NTRAN+(NST-ITMAX)*ITMAX NST=NST+1 GO TO 4 ENDIF ENDIF NPOS=NPOS+NTRAN IF(NPOS.LE.MXOM) THEN MPOS(I)=NPOS ELSE NREC=NREC+1 MPOS(I)=-NREC ENDIF ENDDO C MMM=MIN(NPOS,MXOM) DO N=1,MMM OMEM(N)=TZERO ENDDO C IF(NREC.GT.0) THEN MMEG=NPOS/1000000 MMEG=MMEG+1 NPOS=MXOM WRITE(6,100)MMEG NLEN=MIN(NOMT,(NAST*(NAST+1-2*IONE))/2) OPEN(1,STATUS='SCRATCH',ACCESS='DIRECT',RECL=MZREC*NLEN, X FORM='UNFORMATTED') DO MREC=1,NREC CALL OMWRIT(OMEGA,NLEN,MREC) ENDDO ENDIF C WRITE(6,101)MXOM,NPOS C ELSE C C FINALLY, RETRIEVE AND RETURN OMEGA C IF(MPOS(IE).EQ.0) THEN DO N=1,NOMT OMEGA(N)=TZERO ENDDO GO TO 25 ENDIF C JJJJ=NAST !INC ZEROES NOMWRT.GT.0 C JJJJ=JTMAX !EXC ZEROES NOMWRT.GT.0 C IF(MPOS(IE).GT.0) THEN M=MPOS(IE-1) K=0 DO JT=1+IONE,JJJJ IP=MIN(JT-IONE,ITMAX) DO IT=1,IP M=M+1 IF(NOMWRT.GT.0)THEN K=JT+1-IONE+(JJJJ+1-IONE)*(IT-1)-(IT*(IT+1))/2 ELSE K=K+1 ENDIF OMEGA(K)=OMEM(M) ENDDO IF(M.GE.MPOS(IE)) GO TO 25 ENDDO ELSE MREC=-MPOS(IE) CALL OMREAD(OMEGA,NOMT,MREC) IF(NOMWRT.GT.0)THEN DO I=1,NOMT OMST(I)=OMEGA(I) ENDDO M=0 DO JT=1+IONE,JJJJ IP=MIN(JT-IONE,ITMAX) DO IT=1,IP M=M+1 K=JT+1-IONE+(JJJJ+1-IONE)*(IT-1)-(IT*(IT+1))/2 OMEGA(K)=OMST(M) ENDDO IF(M.GE.NOMT) GO TO 25 ENDDO ENDIF ENDIF ENDIF C C 25 RETURN C 100 FORMAT(//' ****OPENING SCRATCH FILE, COULD AVOID' X ,' BY INCREASING MZMEG TO:',I4//) 101 FORMAT(//' MXOM =', I9,' USED =',I9//) END C C ****************************************************************** C SUBROUTINE OMGIC(IEE,IWJ0,BFORM) IMPLICIT REAL*8 (A-H,O-Z) C C CALCULATE THE COLLISION STRENGTHS FOR TRANSITIONS BETWEEN LEVELS C IN IC COUPLING C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) PARAMETER (CNVEV=13.6058) PARAMETER (EINF=1.0D6) C LOGICAL BFORM C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT) COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/TOPINF/RMAX(NOMLT),TOPR1,TOPR2,ETOP,ITRAN(NOMLT) X,IEMAX(NOMLT),IOMSET COMMON/NRBRND/ELEVR(NMLT) COMMON/WORK1/OMST(NOMLT) C SAVE ICALL C DATA ICALL/1/ C NPOS(I,J,IONE)=((J-1)*(J-2*IONE))/2+I IPOS(I,J,N)=J+N*(I-1)-(I*(I-1))/2 C IWJ=IWJ0 LRGL2=TWO*FJTT(IWJ)+0.1 !JTOT=LRGL2/2 COEFF=(TWO*FJTT(IWJ)+ONE)/TWO C ENGMSH=EMSH(IEE) EC=AZSQ*ENGMSH C NLVOP=0 DO I=1,NLEV IF(EC.GT.ELEV(I)) THEN EINL(I)=EC+ELEV(1)-ELEV(I) NLVOP=NLVOP+1 ENDIF ENDDO C NOMT=(NLVOP*(NLVOP+1-2*IONE))/2 JTMAX=NLVOP IF(NOMWRT.LT.0.AND.-NOMWRT.LT.NOMT)THEN NOMT=-NOMWRT K=0 DO J=1+IONE,NLVOP DO I=1,J-IONE K=K+1 ENDDO IF(K.GE.NOMT)GO TO 11 ENDDO 11 JTMAX=MIN(J,NLVOP) ENDIF IF(NOMWRT.GT.0.AND.NOMWRT.LT.(NLEV*(NLEV+1-2*IONE))/2)THEN K=0 DO J=1+IONE,NLVOP IP=MIN(ITMAX,J-IONE) DO I=1,IP K=K+1 ENDDO ENDDO NOMT=K ENDIF C DO I=NCHOJ(IWJ)+1,NCHPJX(IWJ) ILV=ILVJX(I) EE=ELEV(ILV)/AZSQ-EMSH(IEE) FNU(I)=ONE/SQRT(EE) ENDDO C C INITIALIZE PARTIAL OMEGAS C DO I=1,NOMT OMGPW(I)=TZERO END DO C C INITIALIZE LINE STRENGTH AND BURGESS TOP-UP (ONCE ONLY) C IF(ICALL.EQ.1)CALL DIPIC C IF(IEE.EQ.1)THEN IF(LRGLAM.EQ.LRGL2)CALL TOP1IC(IWJ) ENDIF C C REWIND ANY SCRATCH FILES IF APPROPRIATE C IF(IEE.EQ.1)THEN IF(ITOP.GT.0)THEN IF(IWJ.EQ.NPWJ-1)REWIND(40) IF(IWJ.EQ.NPWJ)REWIND(41) ENDIF IF(ITOP.NE.0.AND.IWJ.EQ.NPWJ)REWIND(42) ENDIF C ICALL=0 IF(NLVOP.LE.IONE)GO TO 270 C IF(NCHOJ(IWJ).LE.IONE) GO TO 270 C C DETERMINE THE PHYSICAL K MATRIX C CALL MQDTK(IEE,IWJ) C C CALCULATE THE T-MATRIX ELEMENTS SQUARED C CALL TMTRIX(NCHOJ(IWJ)) C C BURGESS DIPOLE TOP-UP (CBE OR CC) C CALL TOP2IC(IWJ,NLVOP,LRGL2) C C CALCULATE JPI PARTIAL COLLISION STRENGTHS C DO 90 ILF=1+IONE,JTMAX JFMX=NCHJX(ILF) IF(JFMX.EQ.0) GO TO 90 IP=MIN(ILF-IONE,ITMAX) DO 80 ILI=1,IP JIMX=NCHJX(ILI) IF(JIMX.EQ.0) GO TO 80 N=NPOS(ILI,ILF,IONE) DO JF=1,JFMX !LOOP OVER FINAL CHANNELS ICF=ICHJX(ILF,JF) DO JI=1,JIMX !LOOP OVER INITIAL CHANNELS ICI=ICHJX(ILI,JI) I=IPOS(ICI,ICF,NCHOJ(IWJ)) OMGPW(N)=OMGPW(N)+TM(I) ENDDO ENDDO OMGPW(N)=COEFF*OMGPW(N) 80 ENDDO 90 ENDDO C C WRITE PARTIAL OMEGAS TO SCRATCH FOR NON-DIPOLE TOP C IF(ITOP.GT.0)THEN IF(IWJ.EQ.NPWJ-3)WRITE(40)((OMGPW(NPOS(ILI,ILF,IONE)) X ,ILI=1,MIN(ILF-IONE,ITMAX)),ILF=1+IONE,JTMAX) IF(IWJ.EQ.NPWJ-2)WRITE(41)((OMGPW(NPOS(ILI,ILF,IONE)) X ,ILI=1,MIN(ILF-IONE,ITMAX)),ILF=1+IONE,JTMAX) ENDIF C C ADD IN GEOMETRIC SERIES TOPUP FOR NON-DIPOLE TRANSITIONS C IF(ITOP.NE.0.AND.IWJ.GE.NPWJ-1)CALL TOP3(IEE,IWJ,NOMT) C C STORE PARTIAL COLLISION STRENGTHS C IF(MPOS(IEE).GT.0)THEN K0=MPOS(IEE-1) K=K0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 OMEM(K)=OMEM(K)+OMGPW(K-K0) if(k.gt.mpos(iee))stop 'mis-match of storage in omem' ENDDO ENDDO ELSEIF(MPOS(IEE).LT.0)THEN NREC=-MPOS(IEE) CALL OMREAD(OMEGA,NOMT,NREC) K=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 OMEGA(K)=OMEGA(K)+OMGPW(K) ENDDO ENDDO CALL OMWRIT(OMEGA,NOMT,NREC) ENDIF C C WRITE SUM TO LRGLAM/2-1 TO SCRATCH I.E. WITHOUT TOP-UP C IF(ITOP.NE.0.AND.LRGL2.EQ.LRGLAM-2.AND.JPTT(IWJ).EQ.1)THEN IF(MPOS(IEE).GT.0)THEN K=MPOS(IEE-1) WRITE(42)(OMEM(K+N),N=1,NOMT) ELSEIF(MPOS(IEE).LT.0)THEN WRITE(42)(OMEGA(N),N=1,NOMT) ENDIF ENDIF C C PRINT COLLISION STRENGTHS FOR EACH PARTIAL WAVE C (N.B. LRGL2.EQ.LRGLAM WILL INLCUDE TOP-UP) C IF(IPTOMIC.EQ.2) THEN WRITE(19,320) IEE,IWJ 320 FORMAT(//80('*')/5X,'COLLISION STRENGTHS BETWEEN LEVELS ', 1 'AT ENERGY MESH POINT NUMBER ',I5,5X, 2 'FOR IC PARTIAL WAVE',I3/) N=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 EINEV=CNVEV*EINL(ILI) ETH=CNVEV*(ELEV(ILF)-ELEV(ILI)) WRITE(14,330) ILI,0.5*JLEV(ILI),ILF,0.5*JLEV(ILF),EINEV,ETH 330 FORMAT(//5X,'TRANSITION FROM LEVEL',I3,' J =',F5.1, 1 ' TO LEVEL',I3,' J =',F5.1, 2 /5X,'ELECTRON ENERGY = ',1PE15.8,' EV', 3 /5X,'THRESHOLD ENERGY =',1PE15.8,' EV') WRITE(14,340) 340 FORMAT(/4X,'J',3X,'PI',10X,'OMEGA (IC)') WRITE(14,350) 350 FORMAT(4X,'-',3X,'--',10X,'----------') WRITE(14,360) FJTT(IWJ),JPTT(IWJ),OMGPW(N) 360 FORMAT(1X,F5.1,1X,I2,7X,1PE15.6) ENDDO ENDDO END IF C 270 IF(IWJ.LT.NPWJ)RETURN IF(NLVOP.LE.IONE)GO TO 275 C C DETERMINE TOP-UP FRACTIONS AND OPTIONALLY RESET OMEGAS. C IF(ITOP.NE.0)CALL TOP4(IEE,IWJ,NOMT) C C PRINT TOTAL COLLISION STRENGTHS ONLY C IF(IPTOMIC.EQ.1.OR.IPTOMIC.EQ.2) THEN WRITE(6,300) IEE 300 FORMAT(///80('*')/5X,'TOTAL COLLISION STRENGTHS BETWEEN LEVELS' 1 ,' (IC)'/5X,'AT ENERGY MESH POINT NUMBER ',I5) 305 FORMAT(5X,'FINAL ELECTRON ENERGY =',1PE15.6,' EV') 310 FORMAT(1X,'LEVEL',I3,2X,'JI =',F5.1,4X,'LEVEL',I4, 1 2X,'JF =',F5.1,4X,'OMEGA =',1PE15.6) C IF(MPOS(IEE).GT.0)THEN K=MPOS(IEE-1) DO ILF=1+IONE,JTMAX EFNEV=CNVEV*EINL(ILF) WRITE(6,305) EFNEV IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 WRITE(6,310) ILI,0.5*JLEV(ILI),ILF,0.5*JLEV(ILF) 1 ,OMEM(K) ENDDO ENDDO ELSEIF(MPOS(IEE).LT.0)THEN K=0 DO ILF=1+IONE,JTMAX EFNEV=CNVEV*EINL(ILF) WRITE(6,305) EFNEV IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 WRITE(6,310) ILI,0.5*JLEV(ILI),ILF,0.5*JLEV(ILF) 1 ,OMEGA(K) ENDDO ENDDO ENDIF END IF C C WRITE TOTAL COLLISION STRENGTHS TO FILE OMEGA, AND OMEGANOTOP (ITOP.NE.0) C 275 IF(NOMWRT.EQ.0)RETURN C IF(NOMWRT.LT.0)THEN IF(BFORM)THEN IF(IEE.EQ.1)THEN !GET ROUNDED ENERGIES NREC=1+(NLEV-1)/5 DO N=1,NREC BACKSPACE(20) ENDDO READ(20,370)(ELEVR(I),I=1,NLEV) ENDIF WRITE(20,380)ENGMSH BACKSPACE(20) READ(20,380)E BACKSPACE(20) DO IT=1,NLEV IF(E.LT.ELEVR(IT))GO TO 375 ENDDO IT=NLEV+1 375 NLVOPR=IT-1 IF(NLVOPR.NE.NLVOP)THEN NOMTR=(NLVOPR*(NLVOPR+1-2*IONE))/2 NOMTR=MIN(NOMTR,-NOMWRT) ELSE NOMTR=NOMT ENDIF ENDIF IF(ITOP.NE.0)THEN IF(BFORM)THEN IF(NOMTR.GT.NOMT)THEN DO N=NOMT+1,NOMTR OMST(N)=TZERO ENDDO ENDIF WRITE(23,380) ENGMSH,(OMST(N),N=1,NOMTR) ENDIF IF(.NOT.BFORM)WRITE(23) ENGMSH,(OMST(N),N=1,NOMT) ENDIF CALL OMEG(IEE,EMSH,MXE,EINL,NLEV,OMEGA,-NOMT,IONE) IF(BFORM)THEN IF(NOMTR.GT.NOMT)THEN DO N=NOMT+1,NOMTR OMEGA(N)=TZERO ENDDO ENDIF WRITE(20,380) ENGMSH,(OMEGA(N),N=1,NOMTR) ENDIF IF(.NOT.BFORM)WRITE(20) ENGMSH,(OMEGA(N),N=1,NOMT) ENDIF C IF(NOMWRT.GT.0)THEN IF(ITOP.NE.0)THEN !NEED TO MAP INTERNAL ORDER TO EXTERNAL JJJJ=NLEV !INC ZEROES NOMWRT.GT.0 C JJJJ=JTMAX !EXC ZEROES NOMWRT.GT.0 M=0 DO JT=1+IONE,NLEV IP=MIN(JT-IONE,ITMAX) DO IT=1,IP M=M+1 K=JT+JJJJ*(IT-1)-(IT*(IT-1+2*IONE))/2 OMEGA(K)=OMST(M) ENDDO ENDDO IF(BFORM)WRITE(23,380) ENGMSH,(OMEGA(N),N=1,NOMWRT) IF(.NOT.BFORM)WRITE(23) ENGMSH,(OMEGA(N),N=1,NOMWRT) ENDIF CALL OMEG(IEE,EMSH,MXE,EINL,NLEV,OMEGA,NOMWRT,IONE) IF(BFORM)WRITE(20,380) ENGMSH,(OMEGA(N),N=1,NOMWRT) IF(.NOT.BFORM)WRITE(20) ENGMSH,(OMEGA(N),N=1,NOMWRT) ENDIF 370 FORMAT(5E16.6) 380 FORMAT(1PE14.8,6(1PE11.3)/(14X,6(E11.3))) C IF(IDIPOLE.NE.0.AND.IEE.EQ.MXE) THEN C C PRINT DIPOLE LINE STRENGTHS C IF(IPTOMIC.EQ.1) THEN WRITE(6,385) 385 FORMAT(///80('*')/5X,'DIPOLE LINE STRENGTHS BETWEEN LEVELS' 1 ,' (IC)') N=0 DO ILF=1+IONE,NLEV IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 WRITE(6,390) ILI,0.5*JLEV(ILI),ILF,0.5*JLEV(ILF),SS(N) 390 FORMAT(1X,'LEVEL',I3,2X,'JI =',F5.1,4X,'LEVEL',I4, 1 2X,'JF =',F5.1,4X,'SLINE =',1PE15.6) ENDDO ENDDO END IF C C PRINT OMEGA AT INFINITE ENERGY IN FILE OMEGA/OMEGANOTOP C NOMLT0=(NLEV*(NLEV+1-2*IONE))/2 NOMLT0=MIN(NOMLT0,ABS(NOMWRT)) IF(BFORM)THEN IF(IDIP.GT.0)WRITE(20,380) EINF,(OMGINF(I),I=1,NOMLT0) IF(IDIP.NE.0)WRITE(23,380) EINF,(OMGINF(I),I=1,NOMLT0) ELSE IF(IDIP.GT.0)WRITE(20) EINF,(OMGINF(I),I=1,NOMLT0) IF(IDIP.NE.0)WRITE(23) EINF,(OMGINF(I),I=1,NOMLT0) ENDIF END IF C RETURN END C C ****************************************************************** C SUBROUTINE OMGJK(IEE,IWJ0,BFORM) IMPLICIT REAL*8 (A-H,O-Z) C C CALCULATE THE COLLISION STRENGTHS FOR TRANSITIONS BETWEEN LEVELS C IN JK COUPLING C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) PARAMETER (CNVEV=13.6058) PARAMETER (EINF=1.0D6) C LOGICAL BFORM C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK3/GAMMA(NMFT) COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT) COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/TOPINF/RMAX(NOMLT),TOPR1,TOPR2,ETOP,ITRAN(NOMLT) X,IEMAX(NOMLT),IOMSET COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX COMMON/NRBJKT/ITL(NMLT),FJTL(NMLT),ILLV(NMLT) COMMON/NRBRND/ELEVR(NMLT) COMMON/WORK1/OMST(NOMLT) C SAVE ICALL C DATA ICALL/1/ C NPOS(I,J,IONE)=((J-1)*(J-2*IONE))/2+I IPOS(I,J,N)=J+N*(I-1)-(I*(I-1))/2 C IWJ=IWJ0 LRGL2=TWO*FJTT(IWJ)+0.1 !JTOT=LRGL2/2 COEFF=(TWO*FJTT(IWJ)+ONE)/TWO C ENGMSH=EMSH(IEE) EC=AZSQ*ENGMSH C C STORE INFORMATION FOR EACH JK LEVEL C IL=0 NLVOP=0 DO I=1,NTLS EINC=EC+ENT(1)-ENT(I) IJMX=NJV(I) DO IJ=1,IJMX IF(EC.GT.ENT(I)) NLVOP=NLVOP+1 IL=IL+1 EINL(IL)=EINC ITL(IL)=I ILLV(IL)=IJ FJTL(IL)=FJT(I,IJ) ENDDO ENDDO IF(IL.NE.NLEV)STOP 'TARGET LEVEL MIS-MATCH' C NOMT=(NLVOP*(NLVOP+1-2*IONE))/2 JTMAX=NLVOP IF(NOMWRT.LT.0.AND.-NOMWRT.LT.NOMT)THEN NOMT=-NOMWRT K=0 DO J=1+IONE,NLVOP DO I=1,J-IONE K=K+1 ENDDO IF(K.GE.NOMT)GO TO 11 ENDDO 11 JTMAX=MIN(J,NLVOP) ENDIF IF(NOMWRT.GT.0.AND.NOMWRT.LT.(NLEV*(NLEV+1-2*IONE))/2)THEN K=0 DO J=1+IONE,NLVOP IP=MIN(ITMAX,J-IONE) DO I=1,IP K=K+1 ENDDO ENDDO NOMT=K ENDIF C DO I=NCHOJ(IWJ)+1,NCHPJ(IWJ) ILV=ILVJ(I) EE=ELEV(ILV)/AZSQ-EMSH(IEE) FNU(I)=ONE/SQRT(EE) ENDDO C C INITIALIZE PARTIAL OMEGAS C DO I=1,NOMT OMGPW(I)=TZERO END DO C C INITIALIZE LINE STRENGTH AND BURGESS TOP-UP (ONCE ONLY) C IF(ICALL.EQ.1)CALL DIPJK C IF(IEE.EQ.1)THEN IF(LRGLAM.EQ.LRGL2)CALL TOP1JK(IWJ) ENDIF C C REWIND ANY SCRATCH FILES IF APPROPRIATE C IF(IEE.EQ.1)THEN IF(ITOP.GT.0)THEN IF(IWJ.EQ.NPWJ-1)REWIND(40) IF(IWJ.EQ.NPWJ)REWIND(41) ENDIF IF(ITOP.NE.0.AND.IWJ.EQ.NPWJ)REWIND(42) ENDIF C ICALL=0 IF(NLVOP.LE.IONE) GO TO 270 C IF(NCHOJ(IWJ).LE.IONE) GO TO 270 C C DETERMINE THE PHYSICAL K MATRIX C CALL MQDTK(IEE,IWJ) C C CALCULATE THE T-MATRIX ELEMENTS SQUARED C CALL TMTRIX(NCHOJ(IWJ)) C C BURGESS DIPOLE TOP-UP (CBE OR CC) C CALL TOP2JK(IWJ,NLVOP,LRGL2) C C CALCULATE JPI PARTIAL COLLISION STRENGTHS C ILI=0 DO 90 I=1,NTLS IF(EC.LT.ENT(I)) GO TO 100 NLVI=NJV(I) DO 85 II=1,NLVI ILI=ILI+1 JIMX=NCHJ(I,II) IF(JIMX.EQ.0) GO TO 85 ILF=0 DO 80 IP=1,NTLS IF(EC.LT.ENT(IP)) GO TO 85 NLVF=NJV(IP) DO 75 IF=1,NLVF ILF=ILF+1 IF(ILF.LT.ILI+IONE) GO TO 75 JFMX=NCHJ(IP,IF) IF(JFMX.EQ.0) GO TO 75 N=NPOS(ILI,ILF,IONE) DO JF=1,JFMX !LOOP OVER FINAL CHANNELS ICF=ICHJ(IP,IF,JF) DO JI=1,JIMX !LOOP OVER INITIAL CHANNELS ICI=ICHJ(I,II,JI) IN=IPOS(ICI,ICF,NCHOJ(IWJ)) OMGPW(N)=OMGPW(N)+TM(IN) ENDDO ENDDO OMGPW(N)=COEFF*OMGPW(N) 75 ENDDO 80 ENDDO 85 ENDDO 90 ENDDO C C WRITE PARTIAL OMEGAS TO SCRATCH FOR NON-DIPOLE TOP C 100 IF(ITOP.GT.0)THEN IF(IWJ.EQ.NPWJ-3)WRITE(40)((OMGPW(NPOS(ILI,ILF,IONE)) X ,ILI=1,MIN(ILF-IONE,ITMAX)),ILF=1+IONE,JTMAX) IF(IWJ.EQ.NPWJ-2)WRITE(41)((OMGPW(NPOS(ILI,ILF,IONE)) X ,ILI=1,MIN(ILF-IONE,ITMAX)),ILF=1+IONE,JTMAX) ENDIF C C ADD IN GEOMETRIC SERIES TOPUP FOR NON-DIPOLE TRANSITIONS C IF(ITOP.NE.0.AND.IWJ.GE.NPWJ-1)CALL TOP3(IEE,IWJ,NOMT) C C STORE PARTIAL COLLISION STRENGTHS C IF(MPOS(IEE).GT.0)THEN K0=MPOS(IEE-1) K=K0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 OMEM(K)=OMEM(K)+OMGPW(K-K0) ENDDO ENDDO ELSEIF(MPOS(IEE).LT.0)THEN NREC=-MPOS(IEE) CALL OMREAD(OMEGA,NOMT,NREC) K=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 OMEGA(K)=OMEGA(K)+OMGPW(K) ENDDO ENDDO CALL OMWRIT(OMEGA,NOMT,NREC) ENDIF C C WRITE SUM TO LRGLAM/2-1 TO SCRATCH I.E. WITHOUT TOP-UP C IF(ITOP.NE.0.AND.LRGL2.EQ.LRGLAM-2.AND.JPTT(IWJ).EQ.1)THEN IF(MPOS(IEE).GT.0)THEN K=MPOS(IEE-1) WRITE(42)(OMEM(K+N),N=1,NOMT) ELSEIF(MPOS(IEE).LT.0)THEN WRITE(42)(OMEGA(N),N=1,NOMT) ENDIF ENDIF C C PRINT COLLISION STRENGTHS FOR EACH PARTIAL WAVE C (N.B. LRGL2.EQ.LRGLAM WILL INLCUDE TOP-UP) C IF(IPTOMJK.EQ.2) THEN WRITE(19,320) IEE,IWJ 320 FORMAT(//80('*')/5X,'COLLISION STRENGTHS BETWEEN LEVELS ', 1 'AT ENERGY MESH POINT NUMBER ',I5,5X, 2 'FOR JK PARTIAL WAVE',I3/) N=0 DO ILF=1+IONE,JTMAX ITF=ITL(ILF) IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 EINEV=CNVEV*EINL(ILI) ITI=ITL(ILI) ETH=CNVEV*(ENT(ITF)-ENT(ITI)) WRITE(19,330) ITI,FJTL(ILI),ITF,FJTL(ILF),EINEV,ETH 330 FORMAT(//5X,'TRANSITION FROM TERM NO. ',I2,' J =',F5.1, 1 ' TO TERM NO. ',I2,' J =', 2 F5.1/5X,'ELECTRON ENERGY = ',1PE15.8,'EV'/ 3 5X,'THRESHOLD ENERGY =',1PE15.8,'EV') WRITE(19,340) 340 FORMAT(/4X,'J',3X,'PI',10X,'OMEGA (JK)') WRITE(19,350) 350 FORMAT(4X,'-',3X,'--',10X,'----------') WRITE(19,360) FJTT(IWJ),JPTT(IWJ),OMGPW(N) 360 FORMAT(1X,F5.1,1X,I2,7X,1PE15.6) ENDDO ENDDO END IF C 270 IF(IWJ.LT.NPWJ)RETURN IF(NLVOP.LE.IONE)GO TO 275 C C DETERMINE TOP-UP FRACTIONS AND OPTIONALLY RESET OMEGAS. C IF(ITOP.NE.0)CALL TOP4(IEE,IWJ,NOMT) C C PRINT TOTAL COLLISION STRENGTHS ONLY C IF(IPTOMJK.EQ.1.OR.IPTOMJK.EQ.2) THEN WRITE(6,300) IEE 300 FORMAT(///80('*')/5X,'TOTAL COLLISION STRENGTHS BETWEEN LEVELS' 1 ,' (JK)'/5X,'AT ENERGY MESH POINT NUMBER',I5) 305 FORMAT(5X,'FINAL ELECTRON ENERGY =',1PE15.6,' EV') 310 FORMAT(1X,'TERM ',I2,'.',2X,'JI =',F5.1,4X,'TERM ',I2,'.', 1 2X,'JF =',F5.1,4X,'OMEGA =',1PE15.6) C IF(MPOS(IEE).GT.0)THEN K=MPOS(IEE-1) DO ILF=1+IONE,JTMAX EFNEV=CNVEV*EINL(ILF) WRITE(6,305) EFNEV IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 WRITE(6,310) ITL(ILI),FJTL(ILI),ITL(ILF),FJTL(ILF) 1 ,OMEM(K) ENDDO ENDDO ELSEIF(MPOS(IEE).LT.0)THEN K=0 DO ILF=1+IONE,JTMAX EFNEV=CNVEV*EINL(ILF) WRITE(6,305) EFNEV IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP K=K+1 WRITE(6,310) ITL(ILI),FJTL(ILI),ITL(ILF),FJTL(ILF) 1 ,OMEGA(K) ENDDO ENDDO ENDIF END IF C C WRITE TOTAL COLLISION STRENGTHS TO FILE OMEGA, AND OMEGANOTOP (ITOP.NE.0) C 275 IF(NOMWRT.EQ.0)RETURN C IF(NOMWRT.LT.0)THEN IF(BFORM)THEN IF(IEE.EQ.1)THEN !GET ROUNDED ENERGIES NREC=1+(NLEV-1)/5 DO N=1,NREC BACKSPACE(20) ENDDO READ(20,370)(ELEVR(I),I=1,NLEV) ENDIF WRITE(20,380)ENGMSH BACKSPACE(20) READ(20,380)E BACKSPACE(20) DO IT=1,NLEV IF(E.LT.ELEVR(IT))GO TO 375 ENDDO IT=NLEV+1 375 NLVOPR=IT-1 IF(NLVOPR.NE.NLVOP)THEN NOMTR=(NLVOPR*(NLVOPR+1-2*IONE))/2 NOMTR=MIN(NOMTR,-NOMWRT) ELSE NOMTR=NOMT ENDIF ENDIF IF(ITOP.NE.0)THEN IF(BFORM)THEN IF(NOMTR.GT.NOMT)THEN DO N=NOMT+1,NOMTR OMST(N)=TZERO ENDDO ENDIF WRITE(23,380) ENGMSH,(OMST(N),N=1,NOMTR) ENDIF IF(.NOT.BFORM)WRITE(23) ENGMSH,(OMST(N),N=1,NOMT) ENDIF CALL OMEG(IEE,EMSH,MXE,EINL,NLEV,OMEGA,-NOMT,IONE) IF(BFORM)THEN IF(NOMTR.GT.NOMT)THEN DO N=NOMT+1,NOMTR OMEGA(N)=TZERO ENDDO ENDIF WRITE(20,380) ENGMSH,(OMEGA(N),N=1,NOMTR) ENDIF IF(.NOT.BFORM)WRITE(20) ENGMSH,(OMEGA(N),N=1,NOMT) ENDIF C IF(NOMWRT.GT.0)THEN IF(ITOP.NE.0)THEN !NEED TO MAP INTERNAL ORDER TO EXTERNAL JJJJ=NLEV !INC ZEROES NOMWRT.GT.0 C JJJJ=JTMAX !EXC ZEROES NOMWRT.GT.0 M=0 DO JT=1+IONE,NLEV IP=MIN(JT-IONE,ITMAX) DO IT=1,IP M=M+1 K=JT+JJJJ*(IT-1)-(IT*(IT-1+2*IONE))/2 OMEGA(K)=OMST(M) ENDDO ENDDO IF(BFORM)WRITE(23,380) ENGMSH,(OMEGA(N),N=1,NOMWRT) IF(.NOT.BFORM)WRITE(23) ENGMSH,(OMEGA(N),N=1,NOMWRT) ENDIF CALL OMEG(IEE,EMSH,MXE,EINL,NLEV,OMEGA,NOMWRT,IONE) IF(BFORM)WRITE(20,380) ENGMSH,(OMEGA(N),N=1,NOMWRT) IF(.NOT.BFORM)WRITE(20) ENGMSH,(OMEGA(N),N=1,NOMWRT) ENDIF 370 FORMAT(5E16.6) 380 FORMAT(1PE14.8,6(1PE11.3)/(14X,6(E11.3))) C IF(IDIPOLE.NE.0.AND.IEE.EQ.MXE) THEN C C PRINT DIPOLE LINE STRENGTHS C IF(IPTOMJK.EQ.1) THEN WRITE(6,385) 385 FORMAT(///80('*')/5X,'DIPOLE LINE STRENGTHS BETWEEN LEVELS' 1 ,' (JK)') N=0 DO ILF=1+IONE,NLEV IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 WRITE(6,390) ITL(ILI),FJTL(ILI),ITL(ILF),FJTL(ILF),SS(N) 390 FORMAT(1X,'TERM ',I2,'.',2X,'JI =',F5.1,4X,'TERM ',I2,'.', 1 2X,'JF =',F5.1,4X,'SLINE =',1PE15.6) ENDDO ENDDO END IF C C PRINT OMEGA AT INFINITE ENERGY IN FILE OMEGA/OMEGANOTOP C NOMLT0=(NLEV*(NLEV+1-2*IONE))/2 NOMLT0=MIN(NOMLT0,ABS(NOMWRT)) IF(BFORM)THEN IF(IDIP.GT.0)WRITE(20,380) EINF,(OMGINF(I),I=1,NOMLT0) IF(IDIP.NE.0)WRITE(23,380) EINF,(OMGINF(I),I=1,NOMLT0) ELSE IF(IDIP.GT.0)WRITE(20) EINF,(OMGINF(I),I=1,NOMLT0) IF(IDIP.NE.0)WRITE(23) EINF,(OMGINF(I),I=1,NOMLT0) ENDIF END IF C RETURN END C C************************************************************************ SUBROUTINE OMREAD(OMEGA,NOMT,IE) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION OMEGA(NOMT) C READ(1,REC=IE)OMEGA C RETURN END C************************************************************************ SUBROUTINE OMWRIT(OMEGA,NOMT,IE) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION OMEGA(NOMT) C WRITE(1,REC=IE)OMEGA C RETURN END C C ****************************************************************** C SUBROUTINE ORDER(EN,NORDER,NDIM,IUP) IMPLICIT REAL*8 (A-H,O-Z) C C C ------------------------------------------------------------------ C C RETURNS NORDER(I)=POINTER TO I-TH ENERGY IN EN ARRAY, C IUP=1 FOR ASCENDING ENERGIES, IUP=-1 FOR DESCENDING ENERGIES C C ------------------------------------------------------------------ DIMENSION EN(NDIM),NORDER(NDIM) C ------------------------------------------------------------------ DO 40 K = 1,NDIM J = K J1 = J - 1 IF (J1.EQ.0) GOTO 30 E = EN(J) + IUP*1.0D-7 DO 20 I = 1,J1 IF (J.LT.K) GOTO 10 N = NORDER(I) IF (IUP.GT.0 .AND. E.GT.EN(N)) GOTO 20 IF (IUP.LT.0 .AND. E.LT.EN(N)) GOTO 20 10 CONTINUE NORDER(J) = NORDER(J-1) J = J - 1 20 CONTINUE 30 NORDER(J) = K 40 CONTINUE C END C C ****************************************************************** C SUBROUTINE RDIC0(E0,EINCR) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO READ BASIC IC DATA FROM JBINIC C C INCLUDE 'PARAM' C PARAMETER (TWO=2.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK11/NVEC(NMLT),IVEC(NMTM,NMLT),CVEC(NMTM,NMLT),TOLTCC COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) C C READ BASIC DATA FROM JBINIC C READ(17) NTLS READ(17) (ENT(IT),IT=1,NTLS) READ(17) NLEV,NUMK,IBIGE,NZ,NELC AZ=MAX(NZ-NELC,1) AZSQ=AZ*AZ READ(17) (ELEV(IL),JLEV(IL),JPLEV(IL),IL=1,NLEV) READ(17) (EMESH(IE),IE=1,NUMK) DO IT=1,NTLS READ(17) LT(IT),LST(IT),LPT(IT) END DO DO IL=1,NLEV READ(17) NVEC(IL) READ(17) (IVEC(I,IL),CVEC(I,IL),I=1,NVEC(IL)) END DO READ(17) NPWJ READ(17) (NCHPJ(IWJ),NCHPJX(IWJ),FJTT(IWJ),JPTT(IWJ),IWJ=1,NPWJ) C C OPEN FILE KMTIC.DAT C OPEN(UNIT=18,FILE='kmtic.dat',STATUS='OLD',FORM="UNFORMATTED") C C SET UP THE NEW ENERGY MESH C EMIN=EMESH(1) EMAX=EMESH(NUMK) IEE=0 EM=E0-EINCR EINCRH=EINCR/TWO DO IE=1,MXE EM=EM+EINCR IF(EM.GT.EMIN-EINCRH.AND.EM.LT.EMAX+EINCRH) THEN IEE=IEE+1 EMSH(IEE)=EM END IF END DO MXE=IEE C RETURN END C C ****************************************************************** C SUBROUTINE RDJK0(E0,EINCR) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO READ BASIC JK DATA FROM JBINJK C C INCLUDE 'PARAM' C PARAMETER (TWO=2.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) C C READ BASIC DATA FROM JBINJK C READ(17) NTLS,NUMK,IBIGE,NZ,NELC AZ=MAX(NZ-NELC,1) AZSQ=AZ*AZ READ(17) (ENT(IT),NJV(IT),IT=1,NTLS) DO IT=1,NTLS READ(17) (FJT(IT,IJ),IJ=1,NJV(IT)) END DO READ(17) (EMESH(IE),IE=1,NUMK) DO IT=1,NTLS READ(17) LT(IT),LST(IT),LPT(IT) END DO READ(17) NPWJ READ(17) (NCHPJ(IWJ),FJTT(IWJ),JPTT(IWJ),IWJ=1,NPWJ) C NLEV=0 DO I=1,NTLS DO J=1,NJV(I) NLEV=NLEV+1 JLEV(NLEV)=TWO*FJT(I,J)+0.1 JPLEV(NLEV)=LPT(I) ELEV(NLEV)=ENT(I) END DO END DO C C OPEN FILE KMTJK.DAT C OPEN(UNIT=18,FILE='kmtjk.dat',STATUS='OLD',FORM="UNFORMATTED") C C SET UP THE NEW ENERGY MESH C EMIN=EMESH(1) EMAX=EMESH(NUMK) IEE=0 EM=E0-EINCR EINCRH=EINCR/TWO DO IE=1,MXE EM=EM+EINCR IF(EM.GT.EMIN-EINCRH.AND.EM.LT.EMAX+EINCRH) THEN IEE=IEE+1 EMSH(IEE)=EM END IF END DO MXE=IEE C RETURN END C C ****************************************************************** C SUBROUTINE RDKMNE(IEK,IWJ0,FILNAM,IRD0) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO READ THE K-MATRIX FILE FROM kmtls.001 etc FROM THE C NO-EXCHANGE RUN WITH PARTIAL WAVES BY L AND PI OF THE (N+1)- C ELECTRON SYSTEM AND S OF THE N-ELECTRON SYSTEM AND GENERATE C AND STORE THE K-MATRICES BY L, PI, AND S FOR THE (N+1)-ELECTRON C SYSTEM IN C C RLS(N) -- UNPHYSICAL K-MATRIX ELEMENTS, N IS VECTOR STORE OF C (I,IP,IW) C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CHARACTER FILNAM*6 C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK8/ITM(NMPL,NMCL),LC(NMPL,NMCL) COMMON/BLOK9/NCHTNE(NMPL), 1 ISTOTNE(NMPL),LTTNE(NMPL),LPTTNE(NMPL),ITMNE(NMPL,NMCL), 2 LCNE(NMPL,NMCL),INWNE(NMPL),LSPI(NMPL),INW(NMPL,2), 3 ITMLC(NMPL,NMCL),KV(NMPL,NMCL),INC(NMPL,NMCL),INWX(NMPL),NPWNE COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/HIGHPW/INOEXCH,JMNTWO C c **** parallel **** character*1 filec,filed,fileu common/parainfo/iam,nproc include 'mpif.h' common/blockstop/istop c **** parallel **** c SAVE IWJ,IRD DATA IWJ/0/,IRD/0/ C NPOS(I,J)=(J*(J-1))/2+I C C READ THE UNPHYSICAL K-MATRICES FOR THE LSPI PARTIAL WAVES C REQUIRED FOR THE INPUT JP SYMMETRY. C DO N=1,MXRLS RLS(N)=TZERO END DO C IF(IRD0.GE.50)THEN !MULTIPLE LS SYMMETRY FILES IF(IWJ0.GT.IWJ)THEN !NEW SYMMETRY c **** parallel **** ich0 = ichar('0') ic = iam/100 id = (iam - 100*ic)/10 iu = (iam - 100*ic - 10*id) ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) c **** parallel **** IF(IWJ.GT.0)THEN !CLOSE-OFF OLD FILES DO IL=IRD0+1,IRD CLOSE(UNIT=IL,STATUS='KEEP') ENDDO ENDIF IWJ=IWJ0 NREC=IEK-1 IC0=ICHAR('0') IRD=IRD0 IW0=0 DO IL=1,NPWFSV(IWJ) !OPEN NEW FILES IW=IWF(IL) !IW IS MONOTONICALLY INCREASING DO IWMX=1,INWX(IW) IWOO=INW(IW,IWMX) !IWOO SHOULD BE MONOTONIC IF(IWOO.NE.0) THEN !DCG 7/10/00 (INW(IW,1) CAN BE ZERO) if(iwoo.lt.iw0)stop'rdkmne.f: mis-read on rlsne' IF(IWOO.GT.IW0)THEN IC1=IWOO/100 IC2=(IWOO-100*IC1)/10 IC3=IWOO-100*IC1-10*IC2 IC1=IC1+IC0 IC2=IC2+IC0 IC3=IC3+IC0 IRD=IRD+1 OPEN(UNIT=IRD,FILE=FILNAM//CHAR(IC1)//CHAR(IC2)//CHAR(IC3) X //'.'//filec//filed//fileu,FORM='UNFORMATTED',STATUS='OLD') DO I=1,NREC !MOVE TO CORRECT ENERGY READ(IRD) ENDDO IW0=IWOO ENDIF ENDIF ENDDO ENDDO ENDIF C ELSE IWJ=IWJ0 REWIND(IRD0) NREC=IEK-1 ENDIF C IRD=IRD0 IW0=0 LSPOS(0)=0 DO IL=1,NPWFSV(IWJ) IW=IWF(IL) LSPOS(IL)=LSPOS(IL-1)+(NCHT(IW)*(NCHT(IW)+1))/2 DO IWMX=1,INWX(IW) IWOO=INW(IW,IWMX) IF(IWOO.NE.0) THEN if(iwoo.lt.iw0)stop'rdkmne.f: mis-read on rlsne' IF(IWOO.GT.IW0)THEN IF(IRD0.GE.50)THEN IRD=IRD+1 ELSE NREC=NREC+(IWOO-IW0-1)*NUMK !SKIP BLOCKS NOT NEEDED DO I=1,NREC READ(IRD) ENDDO NREC=NUMK-1 !MOVE TO CORRECT ENERGY ON NEXT SYMMETRY ENDIF NRD=(NCHTNE(IWOO)*(NCHTNE(IWOO)+1))/2 READ(IRD)(A(IC),IC=1,NRD) !N.B. UPPER HALF LS FROM STGF... IW0=IWOO ENDIF DO IC=1,NCHT(IW) K=KV(IW,IC) IWO=INW(IW,K) IF(IWO.EQ.IWOO)THEN ICO=INC(IW,IC) DO ICP=1,IC KP=KV(IW,ICP) IWOP=INW(IW,KP) IF(IWOP.EQ.IWO) THEN ICOP=INC(IW,ICP) ICMN=MIN(ICO,ICOP) ICMX=MAX(ICO,ICOP) NN=NPOS(ICMN,ICMX) N=NPOS(ICP,IC) N=N+LSPOS(IL-1) RLS(N)=A(NN) END IF END DO ENDIF END DO END IF ENDDO ENDDO C C PRINT THE K-MATRIX TO FILE KMTLS.PRNT C IF(IPRTLS.NE.0) THEN DO IL=1,NPWFSV(IWJ) IW=IWF(IL) LSTT=TWO*FSTT(IW)+ONE LTT=FLTT(IW) WRITE(10,300) LTT,LSTT,LPTT(IW) 300 FORMAT(//44X,'K-MATRIX FOR L =',I3,' 2S+1 =',I3,' PARITY =', 1 I3/44X,43('*')) DO I=1,NCHT(IW) WRITE(10,310) I,ITM(IW,I),LC(IW,I), 1 (RLS(LSPOS(IL-1)+NPOS(MIN(I,IP),MAX(I,IP))),IP=1,NCHT(IW)) 310 FORMAT(1X,I2,2I3,1PE12.4,8E13.4/(9X,E12.4,8E13.4)) END DO ENDDO END IF C RETURN END C C ****************************************************************** C SUBROUTINE RDKMTX(IEK,IWJ0,FILNAM,IRD0) IMPLICIT REAL*8 (A-H,O-Z) C C READ THE K-MATRICES IN LS COUPLING FROM kmtls.001,002... FOR THE C IEKTH ENERGY NECESSARY TO GENERATE JP SYMMETRY IWJ0 AND STORE THEM IN C C RLS(N) -- UNPHYSICAL K-MATRIX ELEMENTS, N IS VECTOR STORE OF C (I,IP,IW) C C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CHARACTER FILNAM*6 C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK8/ITM(NMPL,NMCL),LC(NMPL,NMCL) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/HOLDLS/RLSHLD(NMLSE*MXRLS+1) C c **** parallel **** character*1 filec,filed,fileu common/parainfo/iam,nproc include 'mpif.h' common/blockstop/istop c **** parallel **** c SAVE IEK1,IEK2,IWJ DATA IWJ/0/ C NPOS(I,J)=(J*(J-1))/2+I C IF(IRD0.GE.50)THEN !MULTIPLE LS SYMMETRY FILES IF(IWJ0.GT.IWJ)THEN !NEW SYMMETRY c **** parallel **** ich0 = ichar('0') ic = iam/100 id = (iam - 100*ic)/10 iu = (iam - 100*ic - 10*id) ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) c **** parallel **** IF(IWJ.GT.0)THEN !CLOSE-OFF OLD FILES DO IL=1,NPWFSV(IWJ) CLOSE(UNIT=IRD0+IL,STATUS='KEEP') ENDDO ENDIF IWJ=IWJ0 NREC=IEK-1 IC0=ICHAR('0') DO IL=1,NPWFSV(IWJ) !OPEN NEW FILES IW=IWF(IL) !IW IS MONOTONICALLY INCREASING NRD=(NCHT(IW)*(NCHT(IW)+1))/2 LSPOS(IL)=LSPOS(IL-1)+NRD IC1=IW/100 IC2=(IW-100*IC1)/10 IC3=IW-100*IC1-10*IC2 IC1=IC1+IC0 IC2=IC2+IC0 IC3=IC3+IC0 IRD=IRD0+IL OPEN(UNIT=IRD,FILE=FILNAM//CHAR(IC1)//CHAR(IC2)//CHAR(IC3) X //'.'//filec//filed//fileu,FORM='UNFORMATTED',STATUS='OLD') DO I=1,NREC !MOVE TO CORRECT ENERGY READ(IRD) ENDDO ENDDO ENDIF C ELSE !SINGLE FILE OF LS SYMMETRIES IF(IWJ0.GT.IWJ)IEK2=0 IWJ=IWJ0 IF(IEK.LE.IEK2)THEN !HOLDING K-MATRIX M=IEK-IEK1+1 DO IL=1,NPWFSV(IWJ) N1=LSPOS(IL-1)+1 N2=LSPOS(IL) NN=NMLSE*LSPOS(IL-1)+(M-1)*(N2-N1+1) DO N=N1,N2 NN=NN+1 RLS(N)=RLSHLD(NN) ENDDO ENDDO GO TO 100 ELSE !NEED TO READ NEW ENERGY REWIND(IRD0) NREC=IEK-1 !CORRECT ENERGY OF FIRST SYMMETRY IW0=0 IEK1=IEK+1 MM=MIN(NMLSE,NUMK-IEK) IEK2=IEK+MM ENDIF ENDIF C IRD=IRD0 LSPOS(0)=0 DO IL=1,NPWFSV(IWJ) IW=IWF(IL) !IW IS MONOTONICALLY INCREASING IF(IRD0.GE.50)THEN IRD=IRD+1 ELSE NREC=NREC+(IW-IW0-1)*NUMK !SKIP BLOCKS NOT NEEDED DO I=1,NREC READ(IRD) ENDDO NREC=NUMK-MM-1 !MOVE TO CORRECT ENERGY ON NEXT SYMMETRY IW0=IW ENDIF NRD=(NCHT(IW)*(NCHT(IW)+1))/2 LSPOS(IL)=LSPOS(IL-1)+NRD N1=LSPOS(IL-1)+1 N2=LSPOS(IL) READ(IRD)(RLS(N),N=N1,N2) !N.B. UPPER HALF LS FROM STGF... IF(IRD0.LT.50)THEN !READ EXTRA ENERGIES TO SAVE ON I/O N2=NMLSE*LSPOS(IL-1) DO M=1,MM N1=N2+1 N2=N1+NRD-1 READ(IRD)(RLSHLD(N),N=N1,N2) ENDDO ENDIF ENDDO C C PRINT THE K-MATRIX TO FILE KMTLS.PRNT C 100 IF(IPRTLS.NE.0) THEN DO IL=1,NPWFSV(IWJ) IW=IWF(IL) LSTT=TWO*FSTT(IW)+ONE LTT=FLTT(IW) WRITE(10,300) LTT,LSTT,LPTT(IW) 300 FORMAT(//44X,'K-MATRIX FOR L =',I3,' 2S+1 =',I3,' PARITY =', 1 I3/44X,43('*')) DO I=1,NCHT(IW) WRITE(10,310) I,ITM(IW,I),LC(IW,I), 1 (RLS(LSPOS(IL-1)+NPOS(MIN(I,IP),MAX(I,IP))),IP=1,NCHT(IW)) 310 FORMAT(1X,I2,2I3,1PE12.4,8E13.4/(9X,E12.4,8E13.4)) ENDDO ENDDO END IF C RETURN END C C ****************************************************************** C SUBROUTINE RDLS IMPLICIT REAL*8 (A-H,O-Z) C C READ THE LS SYMMETRIES FROM JBINLS AND GENERATE JP SYMMETRIES C C STORE THE FOLLOWING INFORMATION ON LS K-MATRIX ELEMENTS FOR LATER C USE: C C NCHT(IW) -- TOTAL NUMBER OF CHANNELS, OPEN AND CLOSED, FOR C THE PARTIAL WAVE IW C FLTT(IW) -- TOTAL ANGULAR MOMENTUM FOR PARTIAL WAVE NUMBER IW C FSTT(IW) -- TOTAL SPIN FOR PARTIAL WAVE NUMBER IW C LPTT(IW) -- TOTAL PARITY FOR THE PARTIAL WAVE NUMBER IW C FLC(IW,I) -- CONTINUUM ELECTRON ANGULAR MOMENTUM OF THE ITH C CHANNEL OF PARTIAL WAVE IW C ICH(IW,IT,I) -- CHANNEL NUMBER FOR TERM NUMBER IT IN PARTIAL WAVE C IW. I IS JUST AN INDEX FOR EACH CHANNEL C ASSOCIATED WITH TERM NUMBER IT C NCH(IW,IT) -- THE NUMBER OF CHANNELS ASSOCIATED WITH TERM C NUMBER IT IN PARTIAL WAVE IW. C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THOU=1000.) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK8/ITM(NMPL,NMCL),LC(NMPL,NMCL) COMMON/HIGHPW/INOEXCH,JMNTWO C FLTMAX=TZERO FSTMAX=TZERO FLTMIN=THOU LCMN=1000 LCMX=0 NCHPMX=0 IERROR=0 KMX=0 ITSTLC=0 C C LOOP OVER THE LSP PARTIAL WAVES C NPW=0 DO IW=1,1000 C READ(9,end=31) NCHP,ISTOT,LTT,NPTY C IF(NCHP.LT.0) GO TO 31 !TERMINATOR IF(ISTOT.LT.0)THEN !NX DATA FOUND WRITE(6,*)'*** ERROR: INOEXCH.EQ.0 BUT NX DATA FOUND ON' X ,' FILE KMTLS.DAT ***' STOP 'INOEXCH.EQ.0 BUT NX DATA FOUND ON KMTLS.DAT' ENDIF NPW=NPW+1 IF(NPW.GT.NMPL) THEN WRITE(6,205) NMPL,NPW STOP 'STOP TO INCREASE THE PARAMETER NMPL' END IF IF(NCHP.GT.NMCL.AND.NCHP.GT.NCHPMX)THEN IERROR=1 NCHPMX=NCHP END IF DO IC=1,NCHP C READ(9,end=333) ITMRD,LCRD C IF(IC.LE.NMCL) THEN ITM(IW,IC)=ITMRD LC(IW,IC)=LCRD IF(LC(IW,IC).LT.LCMN) LCMN=LC(IW,IC) IF(LC(IW,IC).GT.LCMX) LCMX=LC(IW,IC) END IF ENDDO IF((LCMX+1).GT.INLM) THEN ITSTLC=1 IERROR=1 END IF C C DETERMINE THE CHANNEL NUMBERS FOR EACH TERM FOR THIS LSPI C PARTIAL WAVE C DO IT=1,NTLS K=0 DO IC=1,NCHP IF(IT.EQ.ITM(IW,IC)) THEN K=K+1 IF(K.GT.NMIL) THEN IF(KMX.LT.K) KMX=K IERROR=1 ELSE ICH(IW,IT,K)=IC END IF ENDIF ENDDO NCH(IW,IT)=K ENDDO C IF(IERROR.EQ.0) THEN C NCHT(IW)=NCHP FLTT(IW)=LTT LSTT=IABS(ISTOT) FSTOT=LSTT FSTT(IW)=(FSTOT-ONE)/TWO LPTT(IW)=NPTY C C DETERMINE AND STORE THE RANGE OF TOTAL J VALUES FOR EACH LSPI C PARTIAL WAVE C FJTTMN(IW)=ABS(FLTT(IW)-FSTT(IW)) FJTTMX(IW)=FLTT(IW)+FSTT(IW) IF(FLTT(IW).LT.FLTMIN) FLTMIN=FLTT(IW) IF(FLTT(IW).GT.FLTMAX) FLTMAX=FLTT(IW) IF(FSTT(IW).GT.FSTMAX) FSTMAX=FSTT(IW) C C STORE THE VALUES OF THE CONTINUUM ELECTRON ORBITAL ANGULAR C MOMENTUM FOR EACH CHANNEL AND EACH PARTIAL WAVE C DO I=1,NCHP FLC(IW,I)=LC(IW,I) ENDDO C ENDIF C ENDDO !END PW LOOP C 31 IF(JMNTWO.GE.0) THEN FJTMIN=JMNTWO FJTMIN=FJTMIN/TWO FLTMNV=ABS(FJTMIN-FSTMAX) IF(FLTMIN.GT.FLTMNV) THEN WRITE(6,201) FJTMIN,FLTMNV,FLTMIN,FLTMNV STOP 'STOP IN RDKMTX - PROBLEMS IN KMTLS.DAT' END IF END IF FJTMAX=FLTMAX-FSTMAX C IF(IERROR.NE.0) THEN IF(NCHPMX.GT.0) THEN WRITE(6,312) NMCL,NCHPMX END IF IF(ITSTLC.NE.0) THEN LCMX1=LCMX+1 WRITE(6,313) INLM,LCMX1 END IF IF(KMX.NE.0) THEN WRITE(6,314) NMIL,KMX END IF STOP 'STOP TO INCREASE THE PARAMETER NMIL' END IF C RETURN C 333 NPW=NPW-1 GO TO 31 C 201 FORMAT(/80('*')/ 1 ' IN ORDER TO GENERATE K-MATRICES FOR J = ',F4.1,' AND', 2 ' ABOVE, YOU MUST INCLUDE'/'LS K-MATRICES FOR L = ',F4.1, 3 ' AND ABOVE. THE INPUT LS K-MATRICES HAD A MINIMUM L OF ', 4 F4.1/' REDO YOUR HIGH-L RUN WITH A MINIMUM VALUE', 5 ' OF L = ',F4.1/80('*')/) 205 FORMAT(/'**** THE PARAMETER NMPL MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') 312 FORMAT(/'**** THE PARAMETER NMCL MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') 313 FORMAT(/'**** THE PARAMETER INLM MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') 314 FORMAT(/'**** THE PARAMETER NMIL MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') C END C C ****************************************************************** C SUBROUTINE RDLSNE IMPLICIT REAL*8 (A-H,O-Z) C C READ THE NO-EXCHANGE LS SYMMETRIES FROM JBINLS AND GENERATE JP C C STORE THE FOLLOWING INFORMATION ON LS K-MATRIX ELEMENTS FOR LATER C USE: C C NCHT(IW) -- TOTAL NUMBER OF CHANNELS, OPEN AND CLOSED, FOR C THE PARTIAL WAVE IW C FLTT(IW) -- TOTAL ANGULAR MOMENTUM FOR PARTIAL WAVE NUMBER IW C FSTT(IW) -- TOTAL SPIN FOR PARTIAL WAVE NUMBER IW C LPTT(IW) -- TOTAL PARITY FOR THE PARTIAL WAVE NUMBER IW C FLC(IW,I) -- CONTINUUM ELECTRON ANGULAR MOMENTUM OF THE ITH C CHANNEL OF PARTIAL WAVE IW C ICH(IW,IT,I) -- CHANNEL NUMBER FOR TERM NUMBER IT IN PARTIAL WAVE C IW. I IS JUST AN INDEX FOR EACH CHANNEL C ASSOCIATED WITH TERM NUMBER IT C NCH(IW,IT) -- THE NUMBER OF CHANNELS ASSOCIATED WITH TERM C NUMBER IT IN PARTIAL WAVE IW. C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THOU=1000.) PARAMETER (HALF=0.5) PARAMETER (THREEHF=1.5) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK4/NCHT(NMPL),FLTT(NMPL),FSTT(NMPL), 1 LPTT(NMPL),FJTTMN(NMPL),FJTTMX(NMPL),FLC(NMPL,NMCL), 2 NCH(NMPL,NMTM),ICH(NMPL,NMTM,NMIL) COMMON/BLOK8/ITM(NMPL,NMCL),LC(NMPL,NMCL) COMMON/BLOK9/NCHTNE(NMPL), 1 ISTOTNE(NMPL),LTTNE(NMPL),LPTTNE(NMPL),ITMNE(NMPL,NMCL), 2 LCNE(NMPL,NMCL),INWNE(NMPL),LSPI(NMPL),INW(NMPL,2), 3 ITMLC(NMPL,NMCL),KV(NMPL,NMCL),INC(NMPL,NMCL),INWX(NMPL),NPWNE COMMON/HIGHPW/INOEXCH,JMNTWO C FLTMAX=TZERO FSTMAX=TZERO FLTMNNE=THOU FSTMXNE=TZERO LCMN=1000 LCMX=0 NCHMX=0 NMX=0 C C READ AND STORE THE K-MATRIX ELEMENTS FROM THE NO-EXCHANGE RUN C FOR THIS ENERGY. C C LOOP OVER PARTIAL WAVES C NPWNE=0 DO IWNE=1,1000 C READ(9,end=31) NCHP,ISTOT,LTT,NPTY C IF(NCHP.LT.0)GO TO 31 !TERMINATOR IF(ISTOT.GT.0)THEN !EXCHANGE DATA FOUND WRITE(6,*)'*** ERROR: INOEXCH.NE.0 BUT EXCHANGE DATA' X ,' FOUND ON FILE KMTLS.DAT ***' STOP 'INOEXCH.NE.0 BUT EXCHANGE DATA FOUND ON KMTLS.DAT' ENDIF NPWNE=NPWNE+1 IF(NCHP.GT.NMCL) THEN WRITE(6,200) NMCL,NCHP STOP 'STOP TO INCREASE THE PARAMETER NMCL' END IF IF(NPWNE.LE.NMPL)THEN DO IC=1,NCHP C READ(9,end=333) ITMRD,LCRD C ITMNE(IWNE,IC)=ITMRD LCNE(IWNE,IC)=LCRD IF(LCRD.LT.LCMN) LCMN=LCRD IF(LCRD.GT.LCMX) LCMX=LCRD END DO FSTOT=IABS(ISTOT) FSTOT=(FSTOT-ONE)/TWO IF(FSTOT.GT.FSTMXNE) FSTMXNE=FSTOT FLTOT=LTT IF(FLTOT.LT.FLTMNNE) FLTMNNE=FLTOT NCHTNE(IWNE)=NCHP ISTOTNE(IWNE)=IABS(ISTOT) LTTNE(IWNE)=LTT LPTTNE(IWNE)=NPTY INWNE(IWNE)=IWNE LSPI(IWNE)=10000*LPTTNE(IWNE)+100*LTTNE(IWNE)+ISTOTNE(IWNE) ELSE DO IC=1,NCHP READ(9,end=333) ITMRD,LCRD ENDDO ENDIF END DO C 31 IF(NPWNE.GT.NMPL) THEN WRITE(6,205) NMPL,NPWNE STOP 'STOP TO INCREASE THE PARAMETER NMPL' END IF C FJTMIN=JMNTWO FJTMIN=FJTMIN/TWO FLTMNV=ABS(FJTMIN-FSTMXNE-HALF) IF(FLTMNNE.GT.FLTMNV) THEN WRITE(6,201) FJTMIN,FLTMNV,FLTMNNE,FLTMNV STOP 'STOP IN RDKMNE - PROBLEMS IN KMTLS.DAT' END IF IF(LCMX+1.GT.INLM) THEN LCMX1=LCMX+1 WRITE(6,202) INLM,LCMX1 STOP 'STOP TO INCREASE THE PARAMETER INLM' END IF C C NOW GENERATE THE LS K-MATRICES BY L, PI, AND S OF THE (N+1)- C ELECTRON SYSTEM RATHER THAN BY L AND PI OF THE (N+1)-ELECTRON C SYSTEM AND S FOR THE N-ELECTRON SYSTEM. C DO J=1,2 DO I=1,NMPL INW(I,J)=0 ENDDO ENDDO DO J=1,NMCL DO I=1,NMPL INC(I,J)=0 END DO END DO DO I=1,NMPL INWX(I)=0 ENDDO C IWMX=NPWNE-1 DO IW=1,IWMX KWMX=NPWNE-IW DO KW=1,KWMX IF(LSPI(KW).GT.LSPI(KW+1)) THEN INWSV=INWNE(KW) LSPISV=LSPI(KW) INWNE(KW)=INWNE(KW+1) LSPI(KW)=LSPI(KW+1) INWNE(KW+1)=INWSV LSPI(KW+1)=LSPISV END IF END DO END DO C IW=0 DO IPWNE=1,NPWNE LLP=INT(LSPI(IPWNE)/10000) LLT=INT((LSPI(IPWNE)-10000*LLP)/100) IST=INT((LSPI(IPWNE)-10000*LLP)-100*LLT) IW=IW+1 FLTT(IW)=LLT LPTT(IW)=LLP FSTOT=IST FSTOT=(FSTOT-ONE)/TWO IF(FSTOT.EQ.TZERO) THEN FSTT(IW)=HALF INW(IW,1)=INWNE(IPWNE) INWX(IW)=1 ELSE IF(IW.GT.1)THEN !BREAK IF INTO 2 PARTS FOR STUPID COMPILERS IF(FLTT(IW).EQ.FLTT(IW-1).AND.LPTT(IW).EQ.LPTT(IW-1)) THEN INW(IW-1,2)=INWNE(IPWNE) INWX(IW-1)=2 FSTT(IW)=FSTOT+HALF INW(IW,1)=INWNE(IPWNE) INWX(IW)=1 GO TO 50 ENDIF ENDIF FSTT(IW)=FSTOT-HALF INW(IW,2)=INWNE(IPWNE) INWX(IW)=2 IW=IW+1 FLTT(IW)=LLT LPTT(IW)=LLP FSTT(IW)=FSTOT+HALF INW(IW,1)=INWNE(IPWNE) INWX(IW)=1 50 CONTINUE END IF END DO C NPW=IW IF(NPW.GT.NMPL) THEN WRITE(6,203) NMPL,NPW STOP 'STOP TO INCREASE THE PARAMTER NMPL' END IF C DO IW=1,NPW IC=0 NCHT(IW)=0 DO K=1,2 IF(INW(IW,K).NE.0) THEN IWO=INW(IW,K) NCHT(IW)=NCHT(IW)+NCHTNE(IWO) NCHO=NCHTNE(IWO) DO ICO=1,NCHO IC=IC+1 ITMLC(IW,IC)=100*ITMNE(IWO,ICO)+LCNE(IWO,ICO) KV(IW,IC)=K INC(IW,IC)=ICO END DO END IF END DO IF(NCHMX.LT.NCHT(IW)) NCHMX=NCHT(IW) END DO C DO IW=1,NPW ICMX=NCHT(IW)-1 DO IC=1,ICMX KCMX=NCHT(IW)-IC DO KC=1,KCMX IF(ITMLC(IW,KC).GT.ITMLC(IW,KC+1)) THEN ITMSV=ITMLC(IW,KC) KVSV=KV(IW,KC) INCSV=INC(IW,KC) ITMLC(IW,KC)=ITMLC(IW,KC+1) KV(IW,KC)=KV(IW,KC+1) INC(IW,KC)=INC(IW,KC+1) ITMLC(IW,KC+1)=ITMSV KV(IW,KC+1)=KVSV INC(IW,KC+1)=INCSV END IF END DO END DO END DO C IF(NCHMX.GT.NMCL) THEN WRITE(6,204) NMCL,NCHMX STOP 'STOP TO INCREASE THE PARAMETER NMCL' END IF C DO IW=1,NPW DO IC=1,NCHT(IW) ITM(IW,IC)=INT(ITMLC(IW,IC)/100) LC(IW,IC)=INT(ITMLC(IW,IC)-100*ITM(IW,IC)) FLC(IW,IC)=LC(IW,IC) END DO DO IT=1,NTLS N=0 DO IC=1,NCHT(IW) IF(IT.EQ.ITM(IW,IC)) THEN N=N+1 IF(NMX.LT.N) NMX=N ICH(IW,IT,N)=IC END IF END DO NCH(IW,IT)=N END DO C C DETERMINE AND STORE THE RANGE OF TOTAL J VALUES FOR EACH LSPI C PARTIAL WAVE C FJTTMN(IW)=ABS(FLTT(IW)-FSTT(IW)) FJTTMX(IW)=FLTT(IW)+FSTT(IW) IF(FLTT(IW).GT.FLTMAX) FLTMAX=FLTT(IW) IF(FSTT(IW).GT.FSTMAX) FSTMAX=FSTT(IW) END DO FJTMAX=FLTMAX-FSTMAX C IF(NMX.GT.NMIL) THEN WRITE(6,206) NMIL,NMX STOP 'STOP TO INCREASE THE PARAMETER NMIL' END IF C RETURN C 333 NPWNE=NPWNE-1 GO TO 31 C 200 FORMAT('NMCL EXCEEDED IN READING THE NO-EXCHANGE K-', 1 'MATRICES'/'INCREASE NMCL ABOVE ITS PRESENT VALUE OF ', 2 I4,' TO AT LEAST ',I4) 201 FORMAT(/80('*')/ 1 ' IN ORDER TO GENERATE K-MATRICES FOR J = ',F4.1,' AND', 2 ' ABOVE, YOU MUST INCLUDE'/'LS K-MATRICES FOR L = ',F4.1, 3 ' AND ABOVE. THE INPUT LS K-MATRICES HAD A MINIMUM L OF ', 4 F4.1/' REDO YOUR NO-EXCHANGE RUN WITH A MINIMUM VALUE' 5,' OF L = ',F4.1/80('*')/) 202 FORMAT(/'**** THE PARAMETER INLM MUST BE INCREASED FROM',I3, 1 'TO ', I3) 203 FORMAT('THE NUMBER OF PARTIAL WAVES GENERATED IN RDKMNE IS' 1,' TOO LARGE, INCREASE NMPL FROM ',I3,' TO ',I3) 204 FORMAT('NMCL EXCEEDED IN RDKMNE - INCREASE NMCL FROM ITS', 1 'PRESENT VALUE OF ',I4,' TO AT LEAST ',I4) 205 FORMAT(/'**** THE PARAMETER NMPL MUST BE INCREASED FROM ', 1I3,' TO ',I3,' ****') 206 FORMAT(/'**** THE PARAMETER NMIL MUST BE INCREASED FROM ', 1 I3,' TO ',I3,' ****') C END C C ****************************************************************** C SUBROUTINE RDVEC(NTLSRD) IMPLICIT REAL*8 (A-H,O-Z) C C READ THE TERM-COUPLING CONSTANTS AND STORE THEM C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (ADMN=1.0D-14) C LOGICAL BORT C COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK11/NVEC(NMLT),IVEC(NMTM,NMLT),CVEC(NMTM,NMLT),TOLTCC COMMON/NRBORD/LORDER(NMTM),JORDER(NMLT),NORDER(NMLT),IORDER X,ICORR(NMTM) C DIMENSION IVECRD(NMTM),CVECRD(NMTM),OVLPSV(NMLT) C BORT=.TRUE. !ALLOW ANY ORTHONORMALIZATION C IF(TOLTCC.GE.TZERO)WRITE(6,215)TOLTCC C DO I=1,4 READ(7,*) ENDDO C ICVMAX=0 JCVMAX=0 CVMAX=TZERO JFMIN=0 FNMIN=ONE C C READ TERM COUPLING COEFFICIENTS C DO J=1,NLEV C READ(7,*) READ(7,*)JX,NVECRD,(IVECRD(I),CVECRD(I),I=1,MIN(4,NVECRD)) C 130 FORMAT(2I5,4(I4,F11.7)) IF(NVECRD.GT.NMTM)THEN WRITE(6,207)NMTM,NVECRD STOP 'STOP TO INCREASE THE PARAMETER NMTM' ENDIF IF(NVECRD.GT.4)THEN IMIN=1 19 IMIN=IMIN+4 IMAX=IMIN+3 IF(IMAX.GT.NVECRD)IMAX=NVECRD READ(7,*)(IVECRD(I),CVECRD(I),I=IMIN,IMAX) C 140 FORMAT(10X,4(I4,F11.7)) IF(IMAX.LT.NVECRD)GO TO 19 ENDIF C C TEST TO DETERMINE IF ANY COMPONENTS OF THE TERM COUPLING C COEFFICIENT FOR THIS LEVEL ARISE FROM TERMS ABOVE NTLS. C IF NOT, SIMPLY STORE THE TERM-COUPLING COEFFIECIENTS; IF SO, C ELIMINATE THOSE COMPONENTS AND LATER RE-ORTHONORMALIZE THE TERM- C COUPLING COEFFICIENTS. C JJ=JORDER(J) C IF(NTLS.EQ.NTLSRD.AND.TOLTCC.LT.TZERO)THEN NVEC(JJ)=NVECRD DO I=1,NVECRD IVEC(I,JJ)=LORDER(IVECRD(I)) CVEC(I,JJ)=CVECRD(I) ENDDO ELSE IF(NTLSRD.LT.NTLS)THEN WRITE(6,180) !NOT A NORMAL ROUTE STOP 'STOP IN RDVEC BECAUSE OF PROBLEMS IN TCCDW.DAT' ELSE NVEC(JJ)=0 FNORM=TZERO DO I=1,NVECRD IV=IVECRD(I) IF(ICORR(LORDER(IV)).GT.0.AND.IV.LE.NTLSRD X .AND.ABS(CVECRD(I)).GT.TOLTCC)THEN NVEC(JJ)=NVEC(JJ)+1 IVEC(NVEC(JJ),JJ)=ICORR(LORDER(IV)) FNORM=FNORM+CVECRD(I)**2 CVEC(NVEC(JJ),JJ)=CVECRD(I) ELSE IF(ABS(CVECRD(I)).GT.CVMAX)THEN CVMAX=ABS(CVECRD(I)) JCVMAX=JJ ICVMAX=IVECRD(I) ENDIF ENDIF ENDDO FNORM=SQRT(FNORM) IF(FNORM.LT.FNMIN)THEN FNMIN=FNORM JFMIN=JJ ENDIF IF(BORT)THEN DO I=1,NVEC(JJ) CVEC(I,JJ)=CVEC(I,JJ)/FNORM ENDDO ENDIF ENDIF ENDIF C ENDDO C IF(BORT.AND.(NTLS.NE.NTLSRD.OR.TOLTCC.GE.TZERO))THEN C C NOW PRINT THE OVERLAPS BEFORE ORTHOGONALIZATION C IF(IPTCC.LT.0)THEN WRITE(6,260) DO J=2,NLEV DO K=1,J-1 IF(JLEV(K).EQ.JLEV(J).AND.JPLEV(K).EQ.JPLEV(J))THEN OVLP=TZERO DO I=1,NVEC(J) DO L=1,NVEC(K) IF(IVEC(L,K).EQ.IVEC(I,J))THEN OVLP=OVLP+CVEC(I,J)*CVEC(L,K) GO TO 25 ENDIF ENDDO 25 ENDDO WRITE(6,280)J,K,OVLP ENDIF ENDDO ENDDO ENDIF C C COMPONENTS OF THE TCC'S HAVE BEEN ELIMINATED; THEREFORE, C RE-ORTHONORMALIZE THE TCC'S C DO J=2,NLEV C NVSV=NVEC(J) C DO K=1,J-1 IF(JLEV(K).EQ.JLEV(J).AND.JPLEV(K).EQ.JPLEV(J))THEN OVLPSV(K)=TZERO DO I=1,NVEC(J) DO L=1,NVEC(K) IF(IVEC(L,K).EQ.IVEC(I,J))THEN OVLPSV(K)=OVLPSV(K)+CVEC(I,J)*CVEC(L,K) GO TO 40 ENDIF ENDDO 40 ENDDO DO L=1,NVEC(K) DO I=1,NVEC(J) IF(IVEC(I,J).EQ.IVEC(L,K))GO TO 44 ENDDO ADVC=ABS(OVLPSV(K)*CVEC(L,K)) IF(ADVC.GT.ADMN)THEN NVEC(J)=NVEC(J)+1 IVEC(NVEC(J),J)=IVEC(L,K) CVEC(NVEC(J),J)=TZERO ENDIF 44 ENDDO ENDIF ENDDO C C IF MORE COMPONENTS HAVE BEEN ADDED TO THE JTH EIGENVECTOR, C SORT THEM C IF(NVEC(J).GT.NVSV)THEN IMX=NVEC(J)-1 DO I=1,IMX LMX=NVEC(J)-I DO L=1,LMX IF(IVEC(L,J).GT.IVEC(L+1,J))THEN IVTMP=IVEC(L,J) CVTMP=CVEC(L,J) IVEC(L,J)=IVEC(L+1,J) CVEC(L,J)=CVEC(L+1,J) IVEC(L+1,J)=IVTMP CVEC(L+1,J)=CVTMP ENDIF ENDDO ENDDO ENDIF C C RE-ORTHOGONALIZE C DO K=1,J-1 IF(JLEV(K).EQ.JLEV(J).AND.JPLEV(K).EQ.JPLEV(J))THEN DO I=1,NVEC(J) DO L=1,NVEC(K) IF(IVEC(L,K).EQ.IVEC(I,J))THEN CVEC(I,J)=CVEC(I,J)-OVLPSV(K)*CVEC(L,K) GO TO 48 ENDIF ENDDO 48 ENDDO ENDIF ENDDO C C RE-NORMALIZE C FNORM=TZERO DO I=1,NVEC(J) FNORM=FNORM+CVEC(I,J)**2 ENDDO FNORM=SQRT(FNORM) DO I=1,NVEC(J) CVEC(I,J)=CVEC(I,J)/FNORM ENDDO C ENDDO C C NOW PRINT THE OVERLAPS AFTER RE-ORTHOGONALIZATION C IF(IPTCC.LT.0)THEN WRITE(6,265) DO J=2,NLEV DO K=1,J-1 IF(JLEV(K).EQ.JLEV(J).AND.JPLEV(K).EQ.JPLEV(J))THEN OVLP=TZERO DO I=1,NVEC(J) DO L=1,NVEC(K) IF(IVEC(L,K).EQ.IVEC(I,J))THEN OVLP=OVLP+CVEC(I,J)*CVEC(L,K) GO TO 80 ENDIF ENDDO 80 ENDDO WRITE(6,290)J,K,OVLP ENDIF ENDDO ENDDO ENDIF C ENDIF C C PRINT THE LEVEL ENERGIES AND TERM COUPLING COEFFICIENTS C IF(CVMAX.GT.0.1)WRITE(6,222)CVMAX,JCVMAX,ICVMAX IF(CVMAX.GT.0.3)STOP 'TCC COMPONENT ZEROED-OUT TOO LARGE' IF(IPTCC.NE.0)THEN IF(JFMIN.NE.0)WRITE(6,225)JFMIN,ONE/FNMIN IF(JCVMAX.NE.0)WRITE(6,220)CVMAX,JCVMAX,ICVMAX WRITE(6,200)NLEV WRITE(6,210) DO J=1,NLEV WRITE(6,230)J,JLEV(J),JPLEV(J),ELEV(J) WRITE(6,240)(IVEC(I,J),CVEC(I,J),I=1,NVEC(J)) ENDDO ENDIF C RETURN C 180 FORMAT('THE NUMBER OF TERMS IN TCCDW.DAT IS LESS THAN THE'/ 1 'THE NUMBER OF TERMS IN THE STGF RUN - MUST STOP') 200 FORMAT(//5X,'NLEV=',I5) 207 FORMAT(/'**** THE PARAMETER NMTM MUST BE INCREASED FROM ', 1 I4,' TO ',I4,' ****') 210 FORMAT(/'THE LEVELS ARE:') 215 FORMAT(/'MINIMUM TERM COUPLING COEFFICIENT RETAINED: TOLTCC=' X,1PD8.1) 220 FORMAT(/'THE ABS LARGEST TCC COMPONENT ZEROED-OUT IS:',F10.6/ X'THIS OCCURRED FOR LEVEL',I5,' FOR ORIGINAL TERM NO. CPT',I4/) 222 FORMAT(/'*** WARNING: TCC COMPONENT ZEROED-OUT:',F10.6/ X'THIS OCCURRED FOR LEVEL',I5,' FOR ORIGINAL TERM NO. CPT',I4/ X'THIS TERM SHOULD BE INCLUDED IN THE CC EXPANSION OR ', X'THE LEVEL DROPPED ***') 225 FORMAT(/'THE LARGEST RENORMALIZATION OF TCCS OCCURRED FOR' X,' LEVEL',I5,' THE FACTOR BEING:',F10.6) 230 FORMAT(/I5,5X,'2J=',I5,5X,'PIE=',I5,5X,'ELEV=',F16.7,' RY') 240 FORMAT(4(I5,F11.7)) 260 FORMAT(/'OVERLAPS BEFORE ORTHOGONALIZATION') 265 FORMAT(/'OVERLAPS AFTER ORTHOGONALIZATION') 280 FORMAT('THE OVERLAP BETWEEN EIGENVECTORS ',I3,' AND ',I3,' = ', X 1PE15.5) 290 FORMAT('THE OVERLAP BETWEEN EIGENVECTORS ',I3,' AND ',I3,' = ', X 1PE15.5) C END C C ****************************************************************** C SUBROUTINE READIC(IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO READ PARTIAL WAVE IC DATA FROM JBINIC C C INCLUDE 'PARAM' C C COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) C IWJ=IWJ0 C READ(17) (ILVJX(IC),FLCJX(IC),FKVJX(IC) X ,IC=1,NCHPJX(IWJ)) READ(17) (NCHJX(IL),IL=1,NLEV) DO IL=1,NLEV IJMX=NCHJX(IL) IF(IJMX.GT.0) THEN READ(17) (ICHJX(IL,IJ),IJ=1,IJMX) END IF END DO C RETURN END C C ****************************************************************** C SUBROUTINE READJK(IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO READ PARTIAL WAVE JK DATA FROM JBINJK C C INCLUDE 'PARAM' C C COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) C IWJ=IWJ0 C READ(17) (ITMJ(ICJ),FLCJ(ICJ),FKVJ(ICJ) X ,ILVJ(ICJ),ICJ=1,NCHPJ(IWJ)) DO IT=1,NTLS READ(17) (NCHJ(IT,IJ),IJ=1,NJV(IT)) DO IJ=1,NJV(IT) JIMX=NCHJ(IT,IJ) IF(JIMX.NE.0) THEN READ(17) (ICHJ(IT,IJ,JI),JI=1,JIMX) END IF END DO END DO C RETURN END C C ****************************************************************** C REAL*8 FUNCTION S6J(FJ1,FJ2,FJ3,FL1,FL2,FL3) IMPLICIT REAL*8 (A-H,O-Z) C C CALCULATES 6-J SYMBOL(FJ1,FJ2,FJ3,FL1,FL2,FL3) USING C A 'NESTED' ALGORITHM C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (HALF=0.5) C COMMON /BLOK3/GAMMA(NMFT) C S6J=TZERO I1=NINT(FJ1+FJ2+FJ3+TWO) I2=NINT(FJ1+FJ2-FJ3+ONE) I3=NINT(FJ1+FJ3-FJ2+ONE) I4=NINT(FJ3+FJ2-FJ1+ONE) IF(I2.LT.1.OR.I3.LT.1.OR.I4.LT.1) GO TO 10 D1=GAMMA(I2)+GAMMA(I3)+GAMMA(I4)-GAMMA(I1) I1=NINT(FJ1+FL2+FL3+TWO) I2=NINT(FJ1+FL2-FL3+ONE) I3=NINT(FJ1+FL3-FL2+ONE) I4=NINT(FL3+FL2-FJ1+ONE) IF(I2.LT.1.OR.I3.LT.1.OR.I4.LT.1) GO TO 10 D2=GAMMA(I2)+GAMMA(I3)+GAMMA(I4)-GAMMA(I1) I1=NINT(FJ2+FL1+FL3+TWO) I2=NINT(FJ2+FL1-FL3+ONE) I3=NINT(FJ2+FL3-FL1+ONE) I4=NINT(FL3+FL1-FJ2+ONE) IF(I2.LT.1.OR.I3.LT.1.OR.I4.LT.1) GO TO 10 D3=GAMMA(I2)+GAMMA(I3)+GAMMA(I4)-GAMMA(I1) I1=NINT(FL2+FL1+FJ3+TWO) I2=NINT(FL2+FL1-FJ3+ONE) I3=NINT(FL2+FJ3-FL1+ONE) I4=NINT(FJ3+FL1-FL2+ONE) IF(I2.LT.1.OR.I3.LT.1.OR.I4.LT.1) GO TO 10 D4=GAMMA(I2)+GAMMA(I3)+GAMMA(I4)-GAMMA(I1) DELT=(D1+D2+D3+D4)*HALF K1=NINT(FJ1+FJ2+FJ3) K2=NINT(FJ1+FL2+FL3) K3=NINT(FJ2+FL1+FL3) K4=NINT(FL2+FL1+FJ3) K5=NINT(FJ1+FJ2+FL2+FL1) K6=NINT(FJ2+FL2+FJ3+FL3) K7=NINT(FJ1+FL1+FJ3+FL3) LMIN=MAX0(0,K1,K2,K3,K4) LMAX=MIN0(K5,K6,K7) L=LMIN+1 LM=LMIN-1 M1=L+1 M2=L-K1 M3=L-K2 M4=L-K3 M5=L-K4 M6=K5-LM M7=K6-LM M8=K7-LM DELT=DELT+GAMMA(M1)-GAMMA(M2)-GAMMA(M3)-GAMMA(M4)-GAMMA(M5) X -GAMMA(M6)-GAMMA(M7)-GAMMA(M8) A=ONE IF(LMAX.NE.LMIN)THEN JJ=LMAX DO J=L,LMAX JK=JJ-1 C=(JJ-K1)*(JJ-K2) A=ONE-A*(JJ+1)*(K5-JK)*(K6-JK)*(K7-JK)/(C*(JJ-K3)*(JJ-K4)) JJ=JK ENDDO ENDIF ISH=IABS(MOD(LM+1,2)) ISH=1-ISH-ISH S6J=A* EXP(DELT)*ISH C 10 RETURN END C C ****************************************************************** C SUBROUTINE TARGET(IUP,IPTRM) IMPLICIT REAL*8 (A-H,O-Z) C C SET-UP TARGET TERMS AND LEVELS. C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (CONVCM=1.097373D+05) C CHARACTER LBCD(10)*1 COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK11/NVEC(NMLT),IVEC(NMTM,NMLT),CVEC(NMTM,NMLT),TOLTCC COMMON/NRBORD/LORDER(NMTM),JORDER(NMLT),NORDER(NMLT),IORDER X,ICORR(NMTM) C DATA (LBCD(I),I=1,10)/"S","P","D","F","G","H","I","K","L","M"/ C C READ THE NUMBER OF LS TERMS (NTLS) FROM TERM.DAT AND, IF ITCC.NE.0, C READ THE NUMBER OF LS TERMS (NTLSTCC) INCLUDED IN THE TCCDW.DAT FILE. C IF NTLSTCC IS GREATER THAN NTLS THEN, AFTER THE TCCS ARE C READ, ELIMINATE THOSE COMPONENTS ARISING FROM TERMS NOT C INCLUDED IN TERM.DAT AND RENORMALIZE THE TERM-COUPLING C COEFFICIENTS - THIS IS THE DEFAULT OPERATION, NTLSRD0.EQ.0. C THIS ALLOWS FOR THE CASE OF CI TERMS INTERSPERSED WITH CC TERMS. C C THE FOLLOWING OPTIONS FOR NTLSRD0 ARE SUPPORTED (THE USER MUST ADD C THIS TO LINE 1 OF THE TCCDW.DAT FILE FORMAT I5 AFTER 'LS TERMS'). C NTLSRD=IABS(NTLSRD0) TERM ENERGIES ETC ARE READ FROM TCCDW.DAT. C NTLSRD MUST BE LARGE ENOUGH TO COVER ALL OF THE CC TERMS. C IF NTLSRD0 .GE. 0 THESE ENERGIES ARE COMPARED (AFTER BEING PUT C INTO ENERGY ORDER) WITH THOSE IN TERM.DAT. THOSE THAT DO NOT C MATCH ARE ASSUMED TO BE CORRELATION AND THE TERM NUMBER IS TAGGED C NEGATIVE SO THAT THESE TCC COMPONENTS CAN BE ELIMINATED. C IF NTLSRD0 .LT. 0 THEN IT IS ASSUMED THAT THE USER HAS ALREADY C TAGGED THE (COLUMN 1) TERM NUMBERS NEGATIVE FOR THE CORRELATION C TERMS WITHIN THE FIRST NTLSRD TERMS. C FOR NON-ZERO NTLSRD0 IT IS ASSUMED THAT ALL TERMS ABOVE NTLSRD C ARE CORRELATION AND THERE IS NO NEED TO TAG THEM SO. C C IF TERM.DAT IS NOT PRESENT (CASE ITCC.NE.0 ONLY) THEN USER MUST C EITHER SET NTLSRD=NTLS, AND IF NTLSTCC.GT.NTLS IT IS ASSUMED THE C FIRST NTLS TERMS ARE CC, OR SET NTLSRD<0 AND FLAGGED CORRELATION. C C A NOTE ON USING OBSERVED ENERGIES: C ITCC.EQ.0: IF OBSERVED TERM ENERGIES ARE USED IN STG3 THEN C NOTHING FURTHER IS NEEDED. IF NOT, OBSERVED TERM ENERGIES CAN C REPLACE THE CALCULATED ONES IN FILE TERM.DAT BUT **NOT** IF IT C CHANGES THE ENERGY ORDER OF A TERM. C ITCC.NE.0: IF OBSERVED TERM ENERGIES ARE USED IN STG3 THEN C THEY **MUST** BE ENTERED INTO THE TCCDW.DAT FILE IF THE ENERGY C ORDER OF A TERM CHANGES SINCE TCCDW.DAT WILL STILL CONTAIN THE C CALCULATED TERM ENERGIES. C OBSERVED LEVEL ENERGIES CAN ALSO BE ENTERED INTO TCCDW.DAT TO C REPLACE THE CALCULATED ONES, THE ENERGY ORDER OF A LEVEL CAN C CHANGE. C C IF(ITCC.NE.0)THEN IORDER=0 DO I=1,NMTM LORDER(I)=I ENDDO DO I=1,NMLT JORDER(I)=I ENDDO ELSE IORDER=1 ENDIF C NTLS=0 C READ(13,*,END=1,ERR=1) NTLS0 C NTLS=NTLS0 C 1 IF(NTLS.GT.NMTM) THEN WRITE(6,207) NMTM,NTLS STOP 'STOP TO INCREASE THE PARAMETER NMTM' END IF C IF(ITCC.NE.0) THEN !TERM COUPLING READ(7,113) NTLSTCC,NTLSRD0 NTLSRD=IABS(NTLSRD0) IF(NTLS.EQ.0)THEN !NO term.dat ASSUME NON-CC ALREADY FLAGGED IF(NTLSRD0.GT.0)THEN !ASSUME ALL CC BELOW ALL CORR NTLS=NTLSRD NTLSRD0=0 IF(NTLSTCC.GT.NTLS)NTLSRD=NTLS+1 !FOR RDVEC ELSEIF(NTLSRD0.LT.0)THEN !USER HAS FLAGGED C NTLS=NTLSTCC !AS ONLY,.ELSE TAKE WHAT WE FIND ELSE WRITE(6,*)' ERROR: REQUIRE term.dat FILE IF NTLSRD.EQ.0' STOP ' ERROR: REQUIRE term.dat FILE IF NTLSRD.EQ.0' ENDIF ELSEIF(NTLSRD.LE.NTLS)THEN !MATCH CC TO term.dat IF(NTLSRD0.EQ.0)THEN !CHECK ALL NTLSRD=NTLSTCC NTLSRD0=NTLSRD ELSE NTLSRD=NTLS IF(NTLSTCC.GT.NTLS)NTLSRD=NTLS+1 !FOR RDVEC ENDIF ENDIF IF(NTLSRD.GT.NMTM) THEN WRITE(6,207) NMTM,NTLSRD STOP 'STOP TO INCREASE THE PARAMETER NMTM' END IF ELSE !ALGEBRAIC RECOUPLING IF(NTLS.EQ.0)THEN WRITE(6,*)' ERROR: REQUIRE term.dat FILE IF ITCC.EQ.0' STOP ' ERROR: REQUIRE term.dat FILE IF ITCC.EQ.0' ENDIF NTLSRD=NTLS END IF C C READ INFORMATION ON TARGET LS TERMS & J LEVELS FROM UNIT 7 C C LST(I)--2S+1 FOR THE ITH TARGET TERM C LT(I)--TOTAL ANGULAR MOMENTUM OF ITH TARGET TERM C LPT(I)--PARITY OF ITH TARGET TERM (0-EVEN, 1-ODD) C ENT(I)--ENERGY OF ITH TARGET TERM IN RYDBERGS C C JLEV--2J FOR LEVEL C JPLEV--PARITY OF LEVEL (0-EVEN,1-ODD) C ELEV--ENERGY OF LEVEL IN RYDBERGS C IF(ITCC.NE.0) THEN 202 READ(7,*) READ(7,*) DO I=1,NTLSRD READ(7,*) ICORR(LORDER(I)),LST(LORDER(I)),LT(LORDER(I)) X ,LPT(LORDER(I)),ENT(LORDER(I)) END DO DO I=NTLSRD+1,NTLSTCC READ(7,*) ENDDO READ(7,*) READ(7,*) NLEV IF(NLEV.GT.NMLT) THEN WRITE(6,205) NMLT,NLEV STOP 'STOP TO INCREASE THE PARAMETER NMLT' END IF READ(7,*) READ(7,*) DO I=1,NLEV READ(7,*) IX,JLEV(JORDER(I)),JPLEV(JORDER(I)) X ,ELEV(JORDER(I)) ENDDO IF(IORDER.EQ.0)THEN CALL ORDER(ENT,NORDER,NTLSRD,3) !FOR RYD VS A.U DO I=1,NTLSRD LORDER(NORDER(I))=I !REVERSE INDEX ENDDO CALL ORDER(ELEV,NORDER,NLEV,3) !FOR RYD VS A.U DO I=1,NLEV JORDER(NORDER(I))=I !REVERSE INDEX ENDDO REWIND(7) READ(7,*) IORDER=1 GO TO 202 ENDIF IF(NTLSRD0.GT.0)THEN !MATCH TCCDW.DAT & TERM.DAT TOLE=1.5E-6*SQRT(AZSQ) I=0 J=0 10 READ(13,*)LSTD,LTD,LPTD,ENTD 11 I=I+1 IF(I.GT.NTLSRD)THEN WRITE(6,240)NTLS,NTLSRD STOP' UNABLE TO DETERMINE CC TERMS' ENDIF IF(ABS(ENT(I)-ENTD).GT.TOLE.OR. !CORRELATION X LST(I).NE.LSTD.OR.LT(I).NE.LTD.OR.LPT(I).NE.LPTD)THEN ICORR(I)=-IABS(ICORR(I)) !TAG GO TO 11 ENDIF J=J+1 ICORR(I)=J !POSITION IN CC LIST IF(J.LT.NTLS)GO TO 10 REWIND(13) READ(13,*) ELSEIF(NTLSRD0.LT.0)THEN !LOOK FOR ANY TAGGED J=0 DO I=1,-NTLSRD0 IF(ICORR(I).GT.0)THEN J=J+1 ICORR(I)=J ENDIF ENDDO IF(J*NTLS.NE.NTLS*NTLS)THEN WRITE(6,241)NTLS,-NTLSRD0,J STOP' UNABLE TO DETERMINE CC TERMS' ENDIF ELSE DO I=1,NTLS ICORR(I)=I ENDDO IF(NTLSTCC.GT.NTLS)ICORR(NTLSRD)=-1 !FOR RDVEC ENDIF ENDIF C C PUT CC TERMS IN FIRST NTLS POSITIONS (ALL ITCC) C (TCCDW DOES NOT "NEED" OBSERVED ENERGIES IF NO REORDERING) C DO I=1,NTLS c write(6,*) i,lorder(i),lst(i),lt(i),lpt(i),ent(i) READ(13,*) LST(I),LT(I),LPT(I),ENT(I) END DO DO I=2,NTLS IF(ENT(I).LT.ENT(I-1)-2.5D-7)THEN !FOR SR.ORDER WRITE(6,219)I STOP 'ERROR: TERMS IN TERM.DAT MUST BE IN ENERGY ORDER' ELSEIF(ENT(I).LT.ENT(I-1)+0.5D-7)THEN !RE-SET DEGENERATE ENT(I)=ENT(I-1) ENDIF ENDDO C IF(IPTRM.NE.0) THEN WRITE(6,220) DO I=1,NTLS WRITE(6,230) I,LST(I),LT(I),LPT(I),ENT(I) ENDDO END IF C C DETERMINE AND STORE THE J-VALUES FOR EACH TARGET TERM C DO INJ=1,INJM DO I=1,NMTM JVSV(I,INJ)=0 ENDDO ENDDO IF(IBUG.NE.0) WRITE(6,250) INJMX=0 NJVMX=0 NLEV0=0 DO I=1,NTLS FLSPT=LST(I) FST(I)=(FLSPT-ONE)/TWO FLT(I)=LT(I) FM=MIN(FLT(I),FST(I)) NJV(I)=TWO*FM+ONE+0.1 IF(NJV(I).GT.NJVMX) NJVMX=NJV(I) FJ=ABS(FLT(I)-IUP*FST(I))-IUP*ONE DO J=1,NJV(I) NLEV0=NLEV0+1 FJ=FJ+IUP*ONE FJT(I,J)=FJ INJV=TWO*FJ+ONE+0.1 JVSV(I,INJV)=J IF(INJV.GT.INJMX) INJMX=INJV ENDDO IF(IBUG.NE.0) THEN WRITE(6,260) I,(FJT(I,J),J=1,NJV(I)) ENDIF ENDDO IF(NJVMX.GT.NMLV) THEN WRITE(6,208) NMLV,NJVMX STOP 'STOP TO INCREASE THE PARAMETER NMLV' END IF IF(INJMX.GT.INJM) THEN WRITE(6,209) INJM,INJMX STOP 'STOP TO INCREASE THE PARAMETER INJM' END IF C C WRITE BASIC INPUT FILE FOR ADASEX AND READ TCCS (ITCC.NE.0) C IF(ITCC.NE.0) THEN C CALL RDVEC(NTLSRD) !READ TCCS C IF(NLEV0.NE.NLEV)THEN WRITE(6,216)NLEV0,NLEV WRITE(0,216)NLEV0,NLEV ENDIF C WRITE(22,300)NLEV DO J=1,NLEV CVECMX=TZERO DO I=1,NVEC(J) IF(ABS(CVEC(I,J)).GT.CVECMX) THEN CVECMX=ABS(CVEC(I,J)) ITMNO=IVEC(I,J) END IF END DO LVAL=LT(ITMNO)+1 FJVAL=JLEV(J) FJVAL=FJVAL/TWO EVAL=ELEV(J)*CONVCM WRITE(22,310) LST(ITMNO),LBCD(LVAL),FJVAL,EVAL,JORDER(J) END DO ELSE NLEV=0 WRITE(22,300)NLEV0 DO I=1,NTLS LVAL=LT(I)+1 EVAL=ENT(I)*CONVCM DO J=1,NJV(I) NLEV=NLEV+1 JLEV(NLEV)=TWO*FJT(I,J)+0.1 ELEV(NLEV)=ENT(I) WRITE(22,310) LST(I),LBCD(LVAL),FJT(I,J),EVAL,LORDER(NLEV) END DO END DO END IF WRITE(22,311) CLOSE(UNIT=22) C RETURN C 113 FORMAT(I5,9X,I5) 205 FORMAT(/'**** THE PARAMETER NMLT MUST BE INCREASED FROM', 1 I5,' TO',I5,' ****') 207 FORMAT(/'**** THE PARAMETER NMTM MUST BE INCREASED FROM', 1 I5,' TO',I5,' ****') 208 FORMAT(/'**** THE PARAMETER NMLV MUST BE INCREASED FROM', 1 I5,' TO',I5,' ****') 209 FORMAT(/'**** THE PARAMETER INJM MUST BE INCREASED FROM', 1 I5,' TO',I5,' ****') 216 FORMAT(/'****WARNING***: THE NUMBER OF TARGET LEVELS HAS BEEN', X ' REDUCED FROM',I5,' TO',I5/' THE FORMER IS THE', X ' NUMBER OF LEVELS ASSOCIATED WITH THE CC TERMS WHILE'/ X ' THE LATTER IS THE NUMBER PRESENT IN THE TCCDW.DAT FILE') 219 FORMAT(/'*** ERROR: TERM ENERGIES IN TERM.DAT MUST BE' X ,' IN ENERGY ORDER'/' TERM',I5,' IS OUT OF ORDER') 220 FORMAT(12X,'LS TERMS INCLUDED IN K-MATRICES') 230 FORMAT(3X,I2,'.',3X,'2S+1=',I2,3X,'L=',I2,3X,'PARITY=', 1 I2,3X,'ENERGY(RY)=',1PE15.6) 240 FORMAT('***SR.TARGET: UNABLE TO DETERMINE CC TERMS: ', X 'NTLS=',I5,' NTLSRD=',I5) 241 FORMAT('***SR.TARGET: UNABLE TO DETERMINE CC TERMS: ', X 'NTLS=',I5,' NTLSRD=',I5,' NUMBER OF CC FOUND=',I5) 250 FORMAT(///5X,'J-VALUES FOR EACH TERM') 260 FORMAT(5X,'TERM ',I5,'.',7(5X,F5.1)) 300 FORMAT(' &ADASEX NLEVS=',I5,' &END') 310 FORMAT(6X,'(',I1,A1,')',5X,F5.1,F11.0,I5) 311 FORMAT('NAME:'/'DATE:'/'.') C END C C ****************************************************************** C SUBROUTINE TMTRIX(NCHOP) IMPLICIT REAL*8 (A-H,O-Z) C C SUBROUTINE TO CALCULATE THE T-MATRIX FROM THE K-MATRIX C NON-BLAS VERSION. ALSO CONTAINS (COMMENTED-OUT) A BLAS VERSION C WHICH USES NO ADDITIONAL MEMORY, UNLIKE THE EXPLICIT BLAS VERSION, C BUT IT IS SLOWER. C C T = -2I*R/(1-IR) C TR = 2*R*(1+R**2)^(-1)*R C TI= -2*R*(1+R**2)^(-1) C C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (FOUR=4.0) C COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK5/RLS(MXRLS),RJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/NRBWRK/WORK(NMCJ) C NOMAX=(NCHOP*(NCHOP+1))/2 C CALL MULTS(A,RJ,RJ,NCHOP) C N=1 DO J=NCHOP,1,-1 A(N)=A(N)+ONE N=N+J END DO C CSTRTNBL CALL VERTS(A,-NMCJ,NCHOP,WORK,IERR) IF (IERR.NE.0) THEN WRITE(6,100) STOP 'STOP BECAUSE NO INVERSE FOUND IN TMTRIX' END IF CENDNBL C CSTRTBL CBL CALL DSPTRF('L',NCHOP,A,IPIVOT,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,602) INFO CBL STOP 'FAILURE IN BLAS ROUTINE DSPTRF' CBL ENDIF CBL CALL DSPTRI('L',NCHOP,A,IPIVOT,WORK,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,603) INFO CBL STOP 'FAILURE IN BLAS ROUTINE DSPTRI' CBL ENDIF CENDBL C CALL MULTS(TM,RJ,A,NCHOP) C N=1 DO J=NCHOP,1,-1 A(N)=A(N)-ONE N=N+J END DO C DO N=1,NOMAX TM(N)=FOUR*(TM(N)*TM(N)+A(N)*A(N)) ENDDO C RETURN C 100 FORMAT(' THE MATRIX: 1 + K^2 HAS NO INVERSE IN SUBROUTINE ' 1 ,'TMTRIX - MUST STOP') CB602 FORMAT(//10X,10('*'),' SR.TMTRIX: DSPTRF RETURNED WITH INFO =',I6) CB603 FORMAT(//10X,10('*'),' SR.TMTRIX: DSPTRI RETURNED WITH INFO =',I6) END C C ****************************************************************** C SUBROUTINE TOP1IC(IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE ENERGY-INDEPENDENT PART OF BURGESS TOP-UP C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) PARAMETER (CNVEV=13.6058) PARAMETER (TOLE=1.E-10) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) C SAVE ICALL C DATA ICALL/1/ C IWJ=IWJ0 C IF(ICALL.EQ.1.AND.IPRTOP.EQ.1)THEN WRITE(6,6000)LRGLAM 6000 FORMAT(//10X,'TOP-UP AT (LARGE L).EQ.',I2, + ', (SMALL L).GT.LITLAM'// + 15X,'ALLOWED TRANSITIONS ARE'// + 15X,'INDEX TARGET STATES LITLAM') C DO I=1,INDM WRITE(6,6100)I,NTOP(I,1),NTOP(I,2),LITLAM(I) 6100 FORMAT(I18,I9,I5,I12) ENDDO ENDIF C ICALL=0 C DO I=1,INDM TOPA(I)=TZERO TOPB(I)=TZERO ENDDO C IF(IPRTOP.EQ.1)THEN WRITE(6,6700)0,LRGLAM,JPTT(IWJ) 6700 FORMAT(//12X,'(2S+1) L/2J P =',3I3/12X,24('*')//) WRITE(6,6710) WRITE(6,6720)(I,ILVJX(I),FLCJX(I) X ,2.*FKVJX(I),I=1,NCHPJX(IWJ)) 6710 FORMAT(12X,'CHANNEL',2X,'TARGET',4X,'SMALL', 1 /12X,'INDEX ',3X,'INDEX ',5X,'L',5X,'2K'/) 6720 FORMAT(7X,I8,I9,F11.0,F7.0) C WRITE(6,6200) 6200 FORMAT(/15X,'TOP-UP IN SMALL L FOR'/ + 15X,'INDEX CHANNELS') ENDIF C DO 1000 I=1,INDM LAM=LITLAM(I) ILI=NTOP(I,1) ILF=NTOP(I,2) JIMX=NCHJX(ILI) JFMX=NCHJX(ILF) DO 1520 JF=JFMX,1,-1 !LOOP OVER FINAL CHANNELS ICF=ICHJX(ILF,JF) DO 1510 JI=JIMX,1,-1 !LOOP OVER INITIAL CHANNELS ICI=ICHJX(ILI,JI) C KI=KF IF(ABS(FKVJX(ICI)-FKVJX(ICF)).GT.0.1) GO TO 1510 KJ=2*FKVJX(ICI)+0.1 LLI=FLCJX(ICI)+0.1 !LI CONTINUUM LLF=FLCJX(ICF)+0.1 !LF CONTINUUM DELE=(ELEV(ILF)-ELEV(ILI))/AZSQ DELE=MAX(DELE,TOLE) IF(LLI.EQ.LAM.AND.LLF.EQ.LLI-1.AND.TOPA(I).EQ.TZERO)THEN ITOPA(I)=JPTT(IWJ) NTOPA(I,1)=ICI NTOPA(I,2)=ICF W=DBLE(LRGLAM+1)/DBLE(2*KJ+2) TOPA(I)=ONE/(DBLE(LAM*LAM)*DELE*W* X WSQ2(JLEV(ILI),JLEV(ILF),2*LLI,2*LLF,KJ,ISIGN)) IF(IPRTOP.EQ.1)WRITE(6,6300)I,ICI,ICF,TOPA(I) 6300 FORMAT(I18,I9,I5,10X,'TOPA = ',1PE12.4) GO TO 1000 ENDIF DELE=-DELE IF(LLF.EQ.LAM.AND.LLI.EQ.LLF-1.AND.TOPB(I).EQ.TZERO)THEN ITOPB(I)=JPTT(IWJ) NTOPB(I,1)=ICI NTOPB(I,2)=ICF W=DBLE(LRGLAM+1)/DBLE(2*KJ+2) TOPB(I)=ONE/(DBLE(LAM*LAM)*DELE*W* X WSQ2(JLEV(ILI),JLEV(ILF),2*LLI,2*LLF,KJ,ISIGN)) IF(IPRTOP.EQ.1)WRITE(6,6310)I,ICI,ICF,TOPB(I) 6310 FORMAT(I18,I9,I5,10X,'TOPB = ',1PE12.4) GO TO 1000 ENDIF 1510 ENDDO 1520 ENDDO 1000 ENDDO C RETURN END C C ****************************************************************** C SUBROUTINE TOP1JK(IWJ0) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE ENERGY-INDEPENDENT PART OF BURGESS TOP-UP C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) PARAMETER (CNVEV=13.6058) PARAMETER (TOLE=1.E-10) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/NRBJKT/ITL(NMLT),FJTL(NMLT),ILLV(NMLT) C SAVE ICALL C DATA ICALL/1/ C IWJ=IWJ0 C IF(ICALL.EQ.1.AND.IPRTOP.EQ.1)THEN WRITE(6,6000)LRGLAM 6000 FORMAT(//10X,'TOP-UP AT (LARGE L).EQ.',I2, + ', (SMALL L).GT.LITLAM'// + 15X,'ALLOWED TRANSITIONS ARE'// + 15X,'INDEX TARGET STATES LITLAM') C DO I=1,INDM WRITE(6,6100)I,NTOP(I,1),NTOP(I,2),LITLAM(I) 6100 FORMAT(I18,I9,I5,I12) ENDDO ENDIF C ICALL=0 C DO I=1,INDM TOPA(I)=TZERO TOPB(I)=TZERO ENDDO C IF(IPRTOP.EQ.1)THEN WRITE(6,6700)0,LRGLAM,JPTT(IWJ) 6700 FORMAT(//12X,'(2S+1) L/2J P =',3I3/12X,24('*')//) WRITE(6,6710) WRITE(6,6720)(I,ITMJ(I),ILVJ(I),FLCJ(I) X ,2.*FKVJ(I),I=1,NCHPJ(IWJ)) 6710 FORMAT(12X,'CHANNEL',2X,'TARGET',4X,'SMALL', 1 /12X,'INDEX ',3X,'INDEX ',5X,'L',5X,'2K'/) 6720 FORMAT(7X,I8,I8,I4,F8.0,F7.0) C WRITE(6,6200) 6200 FORMAT(/15X,'TOP-UP IN SMALL L FOR'/ + 15X,'INDEX CHANNELS') ENDIF C DO 1000 I=1,INDM LAM=LITLAM(I) ILI=NTOP(I,1) ILF=NTOP(I,2) ITI=ITL(ILI) ITF=ITL(ILF) ILLI=ILLV(ILI) ILLF=ILLV(ILF) JIMX=NCHJ(ITI,ILLI) JFMX=NCHJ(ITF,ILLF) DO 1520 JF=JFMX,1,-1 !LOOP OVER FINAL CHANNELS ICF=ICHJ(ITF,ILLF,JF) DO 1510 JI=JIMX,1,-1 !LOOP OVER INITIAL CHANNELS ICI=ICHJ(ITI,ILLI,JI) C KI=KF IF(ABS(FKVJ(ICI)-FKVJ(ICF)).GT.0.1)GO TO 1510 KJ=2*FKVJ(ICI)+0.1 LLI=FLCJ(ICI)+0.1 !LI CONTINUUM LLF=FLCJ(ICF)+0.1 !LF CONTINUUM DELE=(ENT(ITF)-ENT(ITI))/AZSQ DELE=MAX(DELE,TOLE) IF(LLI.EQ.LAM.AND.LLF.EQ.LLI-1.AND.TOPA(I).EQ.TZERO)THEN ITOPA(I)=JPTT(IWJ) NTOPA(I,1)=ICI NTOPA(I,2)=ICF W=DBLE(LRGLAM+1)/DBLE(2*KJ+2) TOPA(I)=ONE/(DBLE(LAM*LAM)*DELE*W* X WSQ2(JLEV(ILI),JLEV(ILF),2*LLI,2*LLF,KJ,ISIGN)) IF(IPRTOP.EQ.1)WRITE(6,6300)I,ICI,ICF,TOPA(I) 6300 FORMAT(I18,I9,I5,10X,'TOPA = ',1PE12.4) GO TO 1000 ENDIF DELE=-DELE IF(LLF.EQ.LAM.AND.LLI.EQ.LLF-1.AND.TOPB(I).EQ.TZERO)THEN ITOPB(I)=JPTT(IWJ) NTOPB(I,1)=ICI NTOPB(I,2)=ICF W=DBLE(LRGLAM+1)/DBLE(2*KJ+2) TOPB(I)=ONE/(DBLE(LAM*LAM)*DELE*W* X WSQ2(JLEV(ILI),JLEV(ILF),2*LLI,2*LLF,KJ,ISIGN)) IF(IPRTOP.EQ.1)WRITE(6,6310)I,ICI,ICF,TOPB(I) 6310 FORMAT(I18,I9,I5,10X,'TOPB = ',1PE12.4) GO TO 1000 ENDIF 1510 ENDDO 1520 ENDDO 1000 ENDDO C RETURN END C C ****************************************************************** C SUBROUTINE TOP2IC(IWJ0,NLVOP,LRGL2) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE COULOMB-BETHE OMEGAS AND BURGESS DIPOLE TOP-UP C IN IC COUPLING C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/BLOK12/NCHJX(NMLT),ICHJX(NMLT,NMIJ), 1 ILVJX(NMCJ),FLCJX(NMCJ),FJVJX(NMCJ), 2 FKVJX(NMCJ),KVSV(NMCJ),NCHPJX(NMPJ) COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) C NPOS(I,J,IONE)=((J-1)*(J-2*IONE))/2+I IPOS(I,J,N)=J+N*(I-1)-(I*(I-1))/2 C IWJ=IWJ0 COEFF=(TWO*FJTT(IWJ)+ONE)/TWO C IF(LCBE.GE.0.AND.LRGL2.GE.LCBE.AND.(IPRTOP.EQ.2.OR. X LRGL2.GE.LRGLAM.AND.LRGLAM.GE.0))THEN C C DETERMINE COULOMB-BETHE OMEGAS C IF(IPRTOP.EQ.2)THEN WRITE(6,6700)0,LRGL2,JPTT(IWJ) 6700 FORMAT(//12X,'(2S+1) L/2J P =',3I3/12X,24('*')//) WRITE(6,7200) 7200 FORMAT(3X,'I',3X,'J',4X,' IT',1X,' JT',7X X ,'OMEGA-CC',2X,'OMEGA-CBE',6X,'CBE/CC') ENDIF C DO ICF=1+IONE,NCHOJ(IWJ) ILF=ILVJX(ICF) LLF=FLCJX(ICF)+0.1 KF=2*FKVJX(ICF)+0.1 DO ICI=1,ICF-IONE ILI=ILVJX(ICI) LLI=FLCJX(ICI)+0.1 KI=2*FKVJX(ICI)+0.1 N=NPOS(ILI,ILF,IONE) IF(IDIPX(N).EQ.1.AND.KI.EQ.KF.AND.ABS(LLI-LLF).EQ.1)THEN NT=IPOS(ICI,ICF,NCHOJ(IWJ)) EI=EINL(ILI)/AZSQ EF=EINL(ILF)/AZSQ IFAIL=IPRTOP F=FDIP(EI,LLI,EF,LLF,IFAIL) IF(IFAIL.NE.0.AND.IPRTOP.EQ.2) X WRITE(6,7205)IFAIL,ICI,ICF,EI,LLI,EF,LLF 7205 FORMAT(' FDIP FAILURE: IFAIL=',I2,' FOR I,J=',2I4,' E,L=' X ,2(1PE13.5,I3)) W=FOUR*MAX(LLI,LLF)/(DBLE(KI+1)*THREE) CF2=SS(N)*W*WSQ2(JLEV(ILI),JLEV(ILF),2*LLI,2*LLF,KI,ISIGN) TCBE=FOUR*CF2*F*F IF(IPRTOP.EQ.2)WRITE(6,7210)ICI,ICF,ILI,ILF,COEFF* X TM(NT),COEFF*TCBE,TCBE/TM(NT) 7210 FORMAT(2I4,3X,2I4,4X,2(1PE11.3),0PF12.5) C IF(LRGL2.GE.LRGLAM.AND.LRGLAM.GE.0)THEN C C REPLACE OMEGA-CC BY OMEGA-CBE FOR TOP-UP C TM(NT)=ABS(TCBE) ENDIF ENDIF ENDDO ENDDO ENDIF C IF(LRGL2.GE.LRGLMN.AND.LRGLAM.GE.0)THEN C C ZERO-OUT PARTIALS THAT ARE ALREADY INCLUDED IN BURGESS DIPOLE C TOP-UP C DO I=1,INDM IF(NTOP(I,2).GT.NLVOP)GO TO 100 ILI=NTOP(I,1) ILF=NTOP(I,2) JIMX=NCHJX(ILI) JFMX=NCHJX(ILF) DO JF=JFMX,1,-1 !LOOP OVER FINAL CHANNELS ICF=ICHJX(ILF,JF) DO JI=JIMX,1,-1 !LOOP OVER INITIAL CHANNELS ICI=ICHJX(ILI,JI) IF(MAX(FLCJX(ICI),FLCJX(ICF))-.1.GT.LITLAM(I))THEN NT=IPOS(ICI,ICF,NCHOJ(IWJ)) TM(NT)=TZERO ENDIF ENDDO ENDDO ENDDO ENDIF C 100 IF(LRGL2.EQ.LRGLAM)THEN C C BURGESS DIPOLE TOP-UP C DO I=1,INDM IF(NTOP(I,2).GT.NLVOP)GO TO 200 IF(TOPA(I).NE.TZERO)THEN ILF=NTOP(I,2) ICI=NTOPA(I,1) ICF=NTOPA(I,2) NT=IPOS(ICI,ICF,NCHOJ(IWJ)) TM(NT)=TM(NT)*(ONE+LITLAM(I)**2*EINL(ILF)/AZSQ)*TOPA(I) ENDIF IF(TOPB(I).NE.TZERO)THEN ILI=NTOP(I,1) ICI=NTOPB(I,1) ICF=NTOPB(I,2) NT=IPOS(ICI,ICF,NCHOJ(IWJ)) TM(NT)=TM(NT)*(ONE+LITLAM(I)**2*EINL(ILI)/AZSQ)*TOPB(I) ENDIF ENDDO ENDIF C 200 RETURN END C C ****************************************************************** C SUBROUTINE TOP2JK(IWJ0,NLVOP,LRGL2) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE COULOMB-BETHE OMEGAS AND BURGESS DIPOLE TOP-UP C IN JK COUPLING C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK7/TM(NOMCJ),A(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT) COMMON/BLOKTR/IWF(NMLS),ICFJ(NMCJ,NMLS),FLCJ(NMCJ), 1 CLS(NMCJ,NMLS),ITMJ(NMCJ),NPWFSV(NMPJ),FKVJ(NMCJ) X,FJVJ(NMCJ),ILVJ(NMCJ) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/NRBJKT/ITL(NMLT),FJTL(NMLT),ILLV(NMLT) C NPOS(I,J,IONE)=((J-1)*(J-2*IONE))/2+I IPOS(I,J,N)=J+N*(I-1)-(I*(I-1))/2 C IWJ=IWJ0 COEFF=(TWO*FJTT(IWJ)+ONE)/TWO C IF(LCBE.GE.0.AND.LRGL2.GE.LCBE.AND.(IPRTOP.EQ.2.OR. X LRGL2.GE.LRGLAM.AND.LRGLAM.GE.0))THEN C C DETERMINE COULOMB-BETHE OMEGAS C IF(IPRTOP.EQ.2)THEN WRITE(6,6700)0,LRGL2,JPTT(IWJ) 6700 FORMAT(//12X,'(2S+1) L/2J P =',3I3/12X,24('*')//) WRITE(6,7200) 7200 FORMAT(3X,'I',3X,'J',4X,' IT',1X,' JT',7X X ,'OMEGA-CC',2X,'OMEGA-CBE',6X,'CBE/CC') ENDIF C DO ICF=1+IONE,NCHOJ(IWJ) ILF=ILVJ(ICF) ITF=ITL(ILF) LLF=FLCJ(ICF)+0.1 KF=2*FKVJ(ICF)+0.1 DO ICI=1,ICF-IONE ILI=ILVJ(ICI) ITI=ITL(ILI) LLI=FLCJ(ICI)+0.1 KI=2*FKVJ(ICI)+0.1 N=NPOS(ILI,ILF,IONE) IF(IDIPX(N).EQ.1.AND.KI.EQ.KF.AND.ABS(LLI-LLF).EQ.1)THEN NT=IPOS(ICI,ICF,NCHOJ(IWJ)) EI=EINL(ILI)/AZSQ EF=EINL(ILF)/AZSQ IFAIL=IPRTOP F=FDIP(EI,LLI,EF,LLF,IFAIL) IF(IFAIL.NE.0.AND.IPRTOP.EQ.2) X WRITE(6,7205)IFAIL,ICI,ICF,EI,LLI,EF,LLF 7205 FORMAT(' FDIP FAILURE: IFAIL=',I2,' FOR I,J=',2I4,' E,L=' X ,2(1PE13.5,I3)) W=FOUR*MAX(LLI,LLF)/(DBLE(KI+1)*THREE) CF2=SS(N)*W*WSQ2(JLEV(ILI),JLEV(ILF),2*LLI,2*LLF,KI,ISIGN) TCBE=FOUR*CF2*F*F IF(IPRTOP.EQ.2)WRITE(6,7210)ICI,ICF,ILI,ILF,COEFF* X TM(NT),COEFF*TCBE,TCBE/TM(NT) 7210 FORMAT(2I4,3X,2I4,4X,2(1PE11.3),0PF12.5) C IF(LRGL2.GE.LRGLAM.AND.LRGLAM.GE.0)THEN C C REPLACE OMEGA-CC BY OMEGA-CBE FOR TOP-UP C TM(NT)=ABS(TCBE) ENDIF ENDIF ENDDO ENDDO ENDIF C IF(LRGL2.GE.LRGLMN.AND.LRGLAM.GE.0)THEN C C ZERO-OUT PARTIALS THAT ARE ALREADY INCLUDED IN BURGESS DIPOLE C TOP-UP C DO I=1,INDM IF(NTOP(I,2).GT.NLVOP)GO TO 100 ILI=NTOP(I,1) ITI=ITL(ILI) ILLI=ILLV(ILI) ILF=NTOP(I,2) ITF=ITL(ILF) ILLF=ILLV(ILF) JIMX=NCHJ(ITI,ILLI) JFMX=NCHJ(ITF,ILLF) DO JF=JFMX,1,-1 !LOOP OVER FINAL CHANNELS ICF=ICHJ(ITF,ILLF,JF) DO JI=JIMX,1,-1 !LOOP OVER INITIAL CHANNELS ICI=ICHJ(ITI,ILLI,JI) IF(MAX(FLCJ(ICI),FLCJ(ICF))-.1.GT.LITLAM(I))THEN NT=IPOS(ICI,ICF,NCHOJ(IWJ)) TM(NT)=TZERO ENDIF ENDDO ENDDO ENDDO ENDIF C 100 IF(LRGL2.EQ.LRGLAM)THEN C C BURGESS DIPOLE TOP-UP C DO I=1,INDM IF(NTOP(I,2).GT.NLVOP)GO TO 200 IF(TOPA(I).NE.TZERO)THEN ILF=NTOP(I,2) ICI=NTOPA(I,1) ICF=NTOPA(I,2) NT=IPOS(ICI,ICF,NCHOJ(IWJ)) TM(NT)=TM(NT)*(ONE+LITLAM(I)**2*EINL(ILF)/AZSQ)*TOPA(I) ENDIF IF(TOPB(I).NE.TZERO)THEN ILI=NTOP(I,1) ICI=NTOPB(I,1) ICF=NTOPB(I,2) NT=IPOS(ICI,ICF,NCHOJ(IWJ)) TM(NT)=TM(NT)*(ONE+LITLAM(I)**2*EINL(ILI)/AZSQ)*TOPB(I) ENDIF ENDDO ENDIF C 200 RETURN END C C ****************************************************************** C SUBROUTINE TOP3(IEE,IWJ,NOMT) IMPLICIT REAL*8 (A-H,O-Z) C C SERIES TOPUP FOR NON-DIPOLE TRANSITIONS C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK6/FJTT(NMPJ),JPTT(NMPJ),NCHPJ(NMPJ),NCHOJ(NMPJ), 1 JVSV(NMTM,INJM),ICHN(NMTM,NMLV,INLM,NKVM), 2 NCHJ(NMTM,NMLV),ICHJ(NMTM,NMLV,NMIJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX COMMON/WORK1/OMST(NOMLT) C DIMENSION ITFLG(NOMLT) SAVE ITFLG,IFAST,ICALL C DATA ICALL/999999/ C IF(IEE.LT.ICALL)THEN ! INITIALIZE IF(LRGLAM.LT.0)STOP 'DO NOT USE SERIES FOR DIPOLE TOP-UP' II=0 IF(ITOP.LT.0)II=-1 IFAST=0 NOMLT0=(NLEV*(NLEV+1-2*IONE))/2 DO N=1,NOMLT0 ITFLG(N)=II ENDDO ENDIF C ICALL=IEE C IF(ITOP.GT.0)THEN ! USE OMEGA RATIO IF(IWJ.EQ.NPWJ-1)READ(40)(OMST(N),N=1,NOMT) IF(IWJ.EQ.NPWJ)READ(41)(OMST(N),N=1,NOMT) N=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 IF(ITFLG(N).EQ.0.AND.IDIPX(N).GT.1)THEN AQ=OMGPW(N)/OMST(N) BQ=FJTT(IWJ)/(FJTT(IWJ)+1) BQ=BQ**(2*IDIPX(N)-1) IF(AQ.LT.MIN(BQ,EINL(ILF)/EINL(ILI)))THEN O1=ONE/(ONE-AQ) OMGPW(N)=OMGPW(N)*O1 ELSE ITFLG(N)=IEE IFAST=1 ENDIF ENDIF END DO END DO IF(IFAST.EQ.0)RETURN !QUICK RETURN ENDIF C IF(IPRTOP.EQ.3)WRITE(6,801)EMSH(IEE) 801 FORMAT(//'2**LAM-POLE TOP-UP',3X,'I1',3X,'I2',6X,'AQ',8X,'BQ' X ,8X,'O1',8X,'O2',8X,'O3'/1PE14.8) C ! USE ENERGY RATIO N=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 IF(IDIPX(N).GT.1)THEN C C INTERPOLATE BETWEEN DEGENERATE AND NON-DEGENERATE LIMITS WHEN C L.LT.2*EINL(ILF)/(EINL(ILI)-EINL(ILF)) C IF(ITOP.EQ.-1)THEN AQ=EINL(ILF)/EINL(ILI) O1=ONE+FJTT(IWJ)/(IDIPX(N)-ONE) O1=O1/TWO IF(AQ.GT.0.99)THEN OMGPW(N)=OMGPW(N)*O1 IF(IPRTOP.EQ.3)WRITE(6,803)ILI,ILF,AQ,O1 803 FORMAT(18X,2I5,F10.3,10X,F10.3) ELSE O2=ONE/(ONE-AQ) BQ=AQ*O2 IF(FJTT(IWJ).GT.TWO*BQ)THEN OMGPW(N)=OMGPW(N)*O2 IF(IPRTOP.EQ.3)WRITE(6,804)ILI,ILF,BQ,O2 804 FORMAT(18X,2I5,10X,F10.3,10X,F10.3) ELSE T=FJTT(IWJ)/(BQ*TWO) O3=O2*T+O1*(ONE-T) OMGPW(N)=OMGPW(N)*O3 IF(IPRTOP.EQ.3)WRITE(6,802)ILI,ILF,AQ,BQ,O1,O2,O3 802 FORMAT(18X,2I5,5F10.3) ENDIF ENDIF C C INTERPOLATE BETWEEN DEGENERATE AND NON-DEGENERATE LIMITS WHEN C ENERGY-RATIO EXCEEDS J-RATIO C ELSEIF(ITFLG(N).NE.0)THEN AQ=EINL(ILF)/EINL(ILI) BQ=FJTT(IWJ)/(FJTT(IWJ)+1) BQ=BQ**(2*IDIPX(N)-1) IF(AQ.LT.BQ)THEN O1=ONE/(ONE-AQ) OMGPW(N)=OMGPW(N)*O1 IF(IPRTOP.EQ.3)WRITE(6,803)ILI,ILF,AQ,O1 ELSE O2=ONE+FJTT(IWJ)/(IDIPX(N)-ONE) O2=O2/TWO IF(AQ.LT.ONE)THEN O1=ONE/(ONE-AQ) O3=O1*((ONE-AQ)/(ONE-BQ))**2 X +O2*(AQ-BQ)*(TWO-AQ-BQ)/(ONE-BQ)**2 ELSE O3=O2 ENDIF OMGPW(N)=OMGPW(N)*O3 IF(IPRTOP.EQ.3)WRITE(6,802)ILI,ILF,AQ,BQ,O1,O2,O3 ENDIF ENDIF ENDIF END DO END DO C C WRITE INFO ON SWITCH FROM OMEGA RATIO TO ENERGY/J RATIO C IF(IPRTOP.EQ.0.AND.ITOP.GT.0.AND.IEE.EQ.MXE.AND.IFAST.NE.0)THEN WRITE(6,805)FJTT(IWJ),JPTT(IWJ) 805 FORMAT(//'DETAILS OF ENERGY AT WHICH NON-DIPOLE TRANSITIONS', X' SWITCHED TO USING ENERGY AND/OR J RATIO INSTEAD OF OMEGA RATIO'/ X/'CASE J =',F5.1,' PARITY =',I2//' ILI ILF IE ENERGY'/) N=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 IF(ITFLG(N).GT.0)WRITE(6,806)ILI,ILF,ITFLG(N),EMSH(ITFLG(N)) 806 FORMAT(2I5,I7,F12.7) END DO END DO ENDIF RETURN C END C C ****************************************************************** C SUBROUTINE TOP4(IEE,IWJ,NOMT) IMPLICIT REAL*8 (A-H,O-Z) C C DETERMINE TOP-UP FRACTIONS AND OPTIONALLY RESET OMEGAS. C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C COMMON/BLOK0/NZ,NELC,AZSQ,NUMK,MXE COMMON/BLOK1/FJTMAX,IDIP,IPRTLS,IPRTJK,IPRTIC,LLMN,LLMX,NPW,NTLS, 1 NPWJ,INJMX,LCMN,LCMX,IPTCC,IPTOMJK,IPTOMIC,IBUG,IMODE,ITOP,IDUM, 2 ITCC,IJBIN,IDIPOLE,IONE COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX COMMON/NRBTOP/TOPA(NOMLT),TOPB(NOMLT),OMGINF(NOMLT),NTOPA(NOMLT,2) X,NTOPB(NOMLT,2),ITOPA(NOMLT),ITOPB(NOMLT),NTOP(NOMLT,2) X,LITLAM(NOMLT),LRGLAM,LRGLMN,LCBE,IPRTOP,INDM,NOMWRT,IDIPX(NOMLT) COMMON/TOPINF/RMAX(NOMLT),TOPR1,TOPR2,ETOP,ITRAN(NOMLT) X,IEMAX(NOMLT),IOMSET COMMON/WORK1/OMST(NOMLT) C READ(42)(OMST(N),N=1,NOMT) C IF(TOPR1.EQ.TZERO.AND.TOPR2.LE.TZERO)RETURN !NO TESTS C IF(IEE.EQ.1)THEN !INITIALIZE NOMLT0=(NLEV*(NLEV+1-2*IONE))/2 DO N=1,NOMLT0 ITRAN(N)=0 RMAX(N)=TZERO IEMAX(N)=0 ENDDO ENDIF C IF(EMSH(IEE).GT.ETOP)THEN IF(MPOS(IEE).GT.0)THEN K=MPOS(IEE-1) DO N=1,NOMT IF(OMST(N).GT.TZERO)THEN R=OMEM(K+N)/OMST(N) TOPR=TOPR2 IF(IDIPX(N).EQ.1)THEN TOPR=TOPR1 IF(R.LT.0.99)THEN WRITE(6,777)IEE,N,OMST(N),OMEM(K+N) 777 FORMAT('***ATTENTION: NEGATIVE DIPOLE TOP-UP:' X ,2I7,1PE12.4,E12.4) OMEM(K+N)=OMST(N) ENDIF ENDIF IF(TOPR.GT.TZERO.AND.R.GT.TOPR)THEN IF(IOMSET.EQ.1)OMEM(K+N)=OMST(N)*TOPR IF(ITRAN(N).EQ.0)ITRAN(N)=IEE IF(R.GT.RMAX(N))THEN RMAX(N)=R IEMAX(N)=IEE ENDIF ENDIF ENDIF ENDDO ELSEIF(MPOS(IEE).LT.0)THEN DO N=1,NOMT IF(OMST(N).GT.TZERO)THEN R=OMEGA(N)/OMST(N) TOPR=TOPR2 IF(IDIPX(N).EQ.1)THEN TOPR=TOPR1 IF(R.LT.0.99)THEN WRITE(6,777)IEE,N,OMST(N),OMEGA(N) OMEGA(N)=OMST(N) ENDIF ENDIF IF(TOPR.GT.TZERO.AND.R.GT.TOPR)THEN IF(IOMSET.EQ.1)OMEGA(N)=OMST(N)*TOPR IF(ITRAN(N).EQ.0)ITRAN(N)=IEE IF(R.GT.RMAX(N))THEN RMAX(N)=R IEMAX(N)=IEE ENDIF ENDIF ENDIF ENDDO NREC=-MPOS(IEE) CALL OMWRIT(OMEGA,NOMT,NREC) ENDIF ENDIF C IF(TOPR1.LE.TZERO.AND.TOPR2.LE.TZERO)RETURN C C PRINT DETAILS OF TOP-UP RATIOS C IF(IEE.EQ.MXE)THEN JTR=9999999 NMIN=0 N=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 IF(ITRAN(N).NE.0.AND.ITRAN(N).LT.JTR)THEN NMIN=N IT=ILI JT=ILF JTR=ITRAN(N) ENDIF ENDDO ENDDO IF(NMIN.GT.0)THEN WRITE(6,706)IT,JT,JTR,EMSH(JTR) N=0 DO ILF=1+IONE,JTMAX IP=MIN(ILF-IONE,ITMAX) DO ILI=1,IP N=N+1 IF(ITRAN(N).NE.0)THEN IE=ITRAN(N) II=1 IF(IDIPX(N).EQ.1)II=-1 WRITE(6,707)IE,II*ILI,II*ILF,IEMAX(N),RMAX(N) ENDIF ENDDO ENDDO ENDIF ENDIF C 706 FORMAT(//'THE FIRST TRANSITION AND ENERGY TO FAIL TOP-UP TEST IS:' X' ILI=',I4,2X,'ILF=',I4,2X,'AT',2X,'IE=',I5 X,' ENERGY=',F8.5//'FOR EACH TRANSITION: ENERGY AT WHICH IT ' X,'FIRST FAILED TOP-UP AND ENERGY AND VALUE OF WORST FAILURE' X//' IE',' ILI',' ILF',' IEMAX',' RMAX') 707 FORMAT(I7,2I5,I7,F8.3) C RETURN END C C ****************************************************************** C SUBROUTINE VERTS(V,LV,N,W,IERR) C C ________________________________________________________ C | | C | INVERT A SYMMETRIC MATRIX WITHOUT PIVOTING | C | | C | INPUT: | C | | C | V --ARRAY CONTAINING MATRIX | C | (ONLY THE LOWER HALF NEED BE DEFINED) | C | | C | LV --LEADING (ROW) DIMENSION OF ARRAY V | C | SET .LE. 0 IF V ALREADY PACKED | C | | C | N --MATRIX DIMENSION | C | | C | W --WORK ARRAY WITH LENGTH AT LEAST N | C | | C | OUTPUT: | C | | C | V --INVERSE (IN LOWER HALF ONLY) | C |________________________________________________________| C REAL*8 V(*),W(*),S,T INTEGER G,H,I,J,K,L,M,N,LV,IERR REAL*8 TZERO,ONE C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C IERR=0 C IF(LV.LE.0)GO TO 4 !ALREADY PACKED C ---------------- CNRB |*** PACK V ***| C ---------------- H=LV-N I=0 M=0 L=N G=(N*(N+1))/2 2 IF(L.EQ.G)GO TO 4 K=L+1 M=M+1 L=L+N-M I=I+H+M DO J=K,L V(J)=V(I+J) ENDDO GO TO 2 C 4 H = N K = 1 10 IF ( H .EQ. 1 ) GOTO 40 C -------------------------- C |*** SAVE PIVOT ENTRY ***| C -------------------------- S = V(K) K = K + H G = K H = H - 1 M = H IF ( S .EQ. TZERO ) GOTO 50 J = 0 20 J = J - M M = M - 1 L = G + M T = V(G+J)/S C --------------------------- C |*** ELIMINATE BY ROWS ***| C --------------------------- DO 30 I = G,L 30 V(I) = V(I) - T*V(I+J) G = L + 1 IF ( M .GT. 0 ) GOTO 20 GOTO 10 40 IF ( V(K) .NE. TZERO ) GOTO 60 IERR=2 RETURN 50 IERR=1 RETURN C ------------------------------------------ C |*** SOLVE FOR ROWS OF INVERSE MATRIX ***| C ------------------------------------------ 60 G = N + N DO 150 M = 1,N L = ((G-M)*(M-1))/2 H = L K = M DO 70 I = M,N 70 W(I) = TZERO W(M) = ONE 80 IF ( K .EQ. N ) GOTO 100 T = W(K)/V(K+L) J = L L = L + N - K K = K + 1 IF ( T .EQ. TZERO ) GOTO 80 DO 90 I = K,N 90 W(I) = W(I) - T*V(I+J) GOTO 80 C ----------------------------------- C |*** BACK SUBSTITUTION BY ROWS ***| C ----------------------------------- 100 W(N) = W(N)/V(K+L) 110 IF ( K .EQ. M ) GOTO 130 J = K K = K - 1 L = L + K - N T = W(K) DO 120 I = J,N 120 T = T - W(I)*V(I+L) W(K) = T/V(K+L) GOTO 110 130 DO 140 I = M,N 140 V(I+H) = W(I) 150 CONTINUE C IF(LV.LE.0)RETURN !LEAVE PACKED C ------------------ CNRB |*** UNPACK V ***| C ------------------ H=LV-N I=(N-1)*H+(N*(N-1))/2 M=N-1 L=(N*(N+1))/2 K=L 200 IF(I.EQ.0)RETURN DO J=L,K,-1 V(I+J)=V(J) ENDDO I=I-H-M L=K-1 M=M-1 K=K-N+M GO TO 200 C END C C ****************************************************************** C REAL*8 FUNCTION WSQ2(A,B,C,D,F,ISIGN) C C CALCULATES 3*(2F+1)*W(A,B,C,D,1,F)**2 WHERE W IS RACAH C COEFFICIENT. C NRB: C ISIGN IS THE SIGN OF W BEFORE SQUARING, C *** AND *** A,B,C,D,E,F ARE INPUT TWICE THEIR ACTUAL VALUE. C IMPLICIT INTEGER(A-N) IMPLICIT REAL*8 (W) C C ISIGN=(-1)**((F-A-C)/2) C IF(B-A)10,20,30 C 10 IF(D-C)11,12,13 20 IF(D-C)21,22,23 30 IF(D-C)31,32,33 C 11 WSQ2=.75*DBLE((F+1)*(B+D+F+4)*(B+D+F+6)*(B+D-F+2)*(B+D-F+4))/ + DBLE((B+2)*(B+1)*(B+3)*(D+2)*(D+1)*(D+3))/4. RETURN 12 WSQ2=.75*DBLE((F+1)*(B+D+F+4)*(-B+D+F)*(B-D+F+2)*(B+D-F+2))/ + DBLE((B+2)*(B+1)*(B+3)*D*(D+2)*(D+1))/2. ISIGN=-ISIGN RETURN 13 WSQ2=.75*DBLE((F+1)*(-B+D+F-2)*(-B+D+F)*(B-D+F+2)*(B-D+F+4))/ + DBLE((B+2)*(B+1)*(B+3)*D*(D-1)*(D+1))/4. RETURN C 21 WSQ2=.75*DBLE((F+1)*(B+D+F+4)*(-B+D+F+2)*(B-D+F)*(B+D-F+2))/ + DBLE(B*(B+2)*(B+1)*(D+2)*(D+1)*(D+3))/2. ISIGN=-ISIGN RETURN 22 G=-(B*(B+2)+D*(D+2)-F*(F+2)) IF(G.LT.0)ISIGN=-ISIGN WSQ2=.75*DBLE((F+1)*G*G)/ !(F+1) IS NOT SQUARED-NRB + DBLE(B*(B+2)*(B+1)*D*(D+2)*(D+1)) RETURN 23 WSQ2=.75*DBLE((F+1)*(B+D+F+2)*(-B+D+F)*(B-D+F+2)*(B+D-F))/ + DBLE(B*(B+2)*(B+1)*D*(D-1)*(D+1))/2. RETURN C 31 WSQ2=.75*DBLE((F+1)*(-B+D+F+2)*(-B+D+F+4)*(B-D+F-2)*(B-D+F))/ + DBLE(B*(B-1)*(B+1)*(D+2)*(D+1)*(D+3))/4 RETURN 32 WSQ2=.75*DBLE((F+1)*(B+D+F+2)*(-B+D+F+2)*(B-D+F)*(B+D-F))/ + DBLE(B*(B-1)*(B+1)*D*(D+2)*(D+1))/2. RETURN 33 WSQ2=.75*DBLE((F+1)*(B+D+F)*(B+D+F+2)*(B+D-F-2)*(B+D-F))/ + DBLE(B*(B-1)*(B+1)*D*(D-1)*(D+1))/4. RETURN C END