C PROGRAM PSTGICFDAMP V2.3 NRB 23 OCT 19 C PROGRAM STGICFDAMP V4.6 NRB & DCG 13 APR 22 C C PROGRAM TO READ DAMPED (I.E. COMPLEX) UNPHYSICAL K/S-MATRICES IN LS C COUPLING GENERATED BY STGFDAMP, DETERMINE THE UNPHYSICAL K/S-MATRICES C IN PURE JK COUPLING OR INTERMEDIATE COUPLING AND FINALLY USE THESE TO C DETERMINE AND PRINT THE COLLISION STRENGTHS BETWEEN INDIVIDUAL LEVELS. C THIS PROGRAM WILL READ K/S-MATRICES FROM EITHER THE REGULAR R-MATRIX C CODES OR THE NON-EXCHANGE CODES. C PROGRAM MAIN IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (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 FOR AUGER (CORE) RATES. FILE AUGER. C 5 PROGRAM CONTROL. FILE dstgicfdamp (THIS OPEN MAYBE COMMENTED-OUT) C 6 PRINTED OUTPUT. FILE routicfdamp 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 k/smtls.prnt C 11 DETAILED JK K-MATRIX OUTPUT. FILE k/smtjk.prnt C 12 DETAILED IC K-MATRIX OUTPUT. FILE k/smtic.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-/S-MATRICES (UNPHYSICAL). FILE zk/smtjk/ic.dat COLD 21 LS K-/S-MATRICES (UNPHYSICAL). FILE z/kmtls.dat, smtls.dat CNEW 101,102,... kmtls.001, kmtls.002,... etc by symmetry. C 22 ADAS ADF04 TEMPLATE OUTPUT. FILE adasexj.in.form C 20 TOTAL COLLISION STRENGTH OUTPUT, JK/IC. FILE OMEGA C 23 TOTAL COLLISION STRENGTH WITHOUT TOP-UP, JK/IC. FILE OMEGANOTOP C 24 DR COLLISION STRENGTHS. FILE OMEGDR C 40 SCRATCH (SEQUENTIAL) C 41 SCRATCH (SEQUENTIAL) C 42 SCRATCH (SEQUENTIAL) C C CHARACTER TITLE*80,ELAS*3,PRINT*6,FILNAM*6,KAR*4 LOGICAL EX,BFORM C COMMON/AUGER/AAUGER(NMLT),IAUGER 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,IQDT 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/ZRLS(MXRLS),ZRJ(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/ZTM(NOMCJ),ZA(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,IREADK 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),ARDEC(NMLT),NEWAR,NTYP1 COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/HOLDLS/ZRLSHLD(NMLSE*MXRLS+1) COMMON/INTER/ZRJ1(NOMCJ),ZRJ2(NOMCJ),IEK0 COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/NRBDR/PDR(NMCJ),OMEGDR(NMET,NEMXF),NDRMET COMMON/NRBJKT/ITL(NMLT),FJTL(NMLT),ILLV(NMLT) COMMON/NRBORD/LORDER(NMTM),JORDER(NMLT),NORDER(NMLT),IORDER X,ICORR(NMTM) COMMON/NRBQDT/QNMAX,FNUMAX 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/ZWORK(NMCJ),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 IREADK,NTYP1,NEWAR,IAUGER,NDRMET,IQDT,FNUMAX X,ieq NAMELIST/MESH1/MXE,E0,EINCR,QNMAX,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="dstgicfdamp",FORM="FORMATTED",STATUS='OLD') c **** parallel **** if (iam.eq.0) then OPEN(UNIT=6,FILE="routicfdamp",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="routicfdampg",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,'PSTGICFDAMP v2.3'/34X,16('*')//) READ(5,100) TITLE 100 FORMAT(A80) WRITE(6,200) TITLE 200 FORMAT(/80('*')//A80//80('*')//) C C C READ DATA FROM UNIT 5 C ********************* C C DAMPING VARIABLES: C C NTYP1 = 1 TYPE-I DAMPING SWITCHED-ON (DEFAULT) C = 0 SWITCH-OFF TYPE-I DAMPING. C IAUGER= 0 (DEFAULT) C = 1 READ CORE RATES FROM FILE AUGER: AAUNITS, THEN ITARG, C AA(AAUNITS) WHERE Z-SCALED AA(RYD)=AA(AAUNITS)/AAUNITS C I.E. AAUNITS/Z^2= 1 RYD, 0.5 AU, 2.06706E16 /SEC (1/HBAR). C NDRMET -- THE NUMBER OF INITIAL METASTABLE STATES FOR WHICH DR IS C CALCULATED. (DEFAULT=0). C IQDT = 0 LOOK FOR A K-MX FILE FIRST AND THEN AN S-MX FILE (DEFAULT). C = 1 LOOK FOR UNPHYSICAL S-MX FILE smtls.001,smtls.002,.. C AND STOP IF NOT FOUND. C = 2 LOOK FOR UNPHYSICAL K-MX FILE z/kmtls.001,z/kmtls.002,.. C AND STOP IF NOT FOUND. C QNMAX -- IF .GT. 0 THEN APPLY GAILITIS AVERAGE WHEN FNU .GT. QNMAX C (DEFAULT=-1., NO GAILITIS AVERAGE) C FNUMAX -- CUT-OFF DR AT THIS FNU (IQDT=2 ONLY). DEFAULT=1000000. C C END DAMPING SPECIFIC VARIABLES 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 -- 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/S-MATRIX FILES SPLIT BY SYMMETRY C zk/smtls.001, zk/smtls.002, ... C AND OPENS THEM ON UNIT NUMBERS IRD0+1, IRD0+2 ... C IF .LT. 50 ASSUMES ALL SYMMETRIES IN SAME FILE C (zk/smtls.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 STGICFDAMP: C IN MANY/MOST CASES, THE ONLY DAMPING AT THE HIGH-L OF C AN NX RUN WILL BE NTYP1 I.E. THE K-MATRIX IS REAL SO, C SET INOEXCH.GT.0 TO READ DAMPED K-MATRICES ON zkmls.001, C zkmls.002,...(PRODUCED BY A STGFDAMP RUN.) C SET INOEXCH.LT.0 TO READ UNDAMPED K-MATRICES ON kmtls.001, C kmtls.002,...(PRODUCED BY A STGF RUN.) C (DEFAULT: 0) C IREADK -- AN EXCHANGE RUN WITH ONLY AUGER DAMPING ONLY NEEDS TO C READ REAL K-MATRIX. IREADK .GT./.LT.0 ACTS AS INOEXCH THEN. C DEFAULT: READS DAMPED K-MATRICES. TO AVOID CONFLICT IN NX C CASE, IF JREADK IS SET THEN IT MUST BE SAME SIGN AS INOEXCH. 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 S/K-MATRICES IN LS COUPLING TO C FILE S/KMTLS.PRNT, ONLY APPLIES WHEN IMODE.EQ.0 C (DEFAULT: 0) C IPRTJK -- IF NOT EQUAL TO ZERO PRINT S/K-MATRICES IN JK COUPLING TO C FILE S/KMTJK.PRNT, ONLY APPLIES WHEN IMODE.EQ.0 C (DEFAULT 0). C IPRTIC -- IF NOT EQUAL TO ZERO PRINT S/K-MATRICES IN IC COUPLING TO C FILE S/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 NEWAR = 1 FORCES RECALCULATION OF ARDEC AT EVERY ENERGY KEEPING C ONLY FINAL STATES BELOW THE IONIZATION LIMIT (DEFAULT). C = 0 DOES NOT, AND SUMS OVER ALL ARAD. C THIS IS ***EXTREMELY DANGEROUS*** FOR DR C SINCE RADIATIVE TRANSITIONS TO AUTOIONIZING STATES C ARE THEN COUNTED AS RECOMBINED. THIS CAN LEAD TO A VERY C LARGE OVER-ESTIMATE OF THE DR CROSS SECTION. C ***HOWEVER***, IF YOU ARE INTERESTED IN (DAMPED) C EXCITATION, THEN NEWAR=0 MIGHT BE MORE APPROPRIATE. 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 IREADK=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' NEWAR=1 NTYP1=1 IAUGER=0 NDRMET=0 IQDT=0 FNUMAX=1000000. QNMAX=-1.0 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 IF(INOEXCH*IREADK.LT.0)THEN WRITE(6,*)' *** CONFLICT BETWEEN INOEXCH AND IREADK SETTINGS:' X ,INOEXCH,IREADK STOP ' *** CONFLICT BETWEEN INOEXCH AND IREADK SETTINGS:' ELSE IF(INOEXCH.NE.0)IREADK=INOEXCH ENDIF C IF(IRD0.LT.50)IRD0=21 !SET zk/smtls.dat UNIT C IQDT=MOD(IQDT,100) !JUST IN CASE... 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,IQDT,NTYP1,IAUGER,NEWAR 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// X'IQDT =',I2,2X,'NTYP1 =',I2,2X,'IAUGER =',I2,2X,'NEWAR =',I2///) C BFORM=PRINT.EQ.'FORM'.OR.PRINT.EQ.'form' C IF(NDRMET.GT.NMET)THEN WRITE(6,*)' ***DIMENSION ERROR: NDRMET= ',NDRMET,' .GT. NMET= ' X ,NMET STOP ' ***DIMENSION ERROR: NDRMET .GT. NMET' ENDIF IF(NDRMET.GT.0.AND.ELAS.NE.'YES')THEN WRITE(6,209) 209 FORMAT(/'*** ATTENTION: OMEGA FILE CONTAINS ELASTIC ' X ,'TRANSITIONS BECAUSE DR IS SWITCHED-ON') ELAS='YES' !ENSURE ELASTIC PRESENT FOR DR ENDIF C IF(NTYP1.LT.0)NTYP1=0 IF(NTYP1.EQ.0)WRITE(6,212) 212 FORMAT(/'*** WARNING: NTYP1.EQ. 0, TYPE-I DAMPING SWITCHED-OFF'/) C IF(NEWAR.LT.0)NEWAR=0 IF(NEWAR.EQ.0)WRITE(6,211) 211 FORMAT(/'*** WARNING: NEWAR.EQ.0, DAMPING TO ALL FINAL STATES,' X,' BOUND AND AUTOIONIZING'/' *** DR LIKELY TO BE OVERESTIMATED'/) 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(NTYP1.GT.0)IDIPOLE=1 !DAMPING ON STGICFDAMP C IF(IQDT.EQ.1)THEN IF(IPRTLS.NE.0) X OPEN(UNIT=10,FILE="smtls.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") IF(IPRTJK.NE.0) X OPEN(UNIT=11,FILE="smtjk.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") IF(IPRTIC.NE.0) X OPEN(UNIT=12,FILE="smtic.prnt"//filec//filed//fileu, + FORM="FORMATTED",STATUS="UNKNOWN") ENDIF IF(IQDT.EQ.2)THEN 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") ENDIF 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 WHICH FILES TO USE FOR UNPHYSICAL MATRICES AND IQDT C CONSISTENTCY. C OPEN(UNIT=9,FILE="jbinls"//filec//filed//fileu, !par + FORM="UNFORMATTED",STATUS='OLD') C IF(IRD0.GE.50)KAR='001.' IF(IRD0.LT.50)KAR='dat.' IF(IREADK.GE.0)THEN IF(IQDT.NE.1)THEN FILNAM="zkmls." INQUIRE(FILE=FILNAM//KAR//filec//filed//fileu,EXIST=EX) !par IF(EX)THEN IQDT=2 ELSE IF(IQDT.EQ.2)THEN WRITE(6,*)'*** ERROR: IQDT=2 BUT K-MX FILE zkmls.xxx' X ,' NOT FOUND' STOP'ERROR: IQDT=2 BUT K-MX FILE zkmls.xxx NOT FOUND' ENDIF ENDIF ENDIF IF(IQDT.NE.2)THEN FILNAM="smtls." INQUIRE(FILE=FILNAM//KAR//filec//filed//fileu,EXIST=EX) !par IF(EX)THEN IQDT=1 ELSE IF(IQDT.EQ.1)THEN WRITE(6,*)'*** ERROR: IQDT=1 BUT S-MX FILE smtls.xxx' X ,' NOT FOUND' STOP'ERROR: IQDT=1 BUT S-MX FILE smtls.xxx NOT FOUND' ELSE WRITE(6,*)'*** ERROR: IQDT=',IQDT,' BUT NO LS K-/S- ' X ,'MATRIX DATA FOUND...' STOP 'CANNOT FIND AN LS K-/S- MATRIX FILE' ENDIF ENDIF ENDIF ELSE FILNAM="kmtls." INQUIRE(FILE=FILNAM//KAR//filec//filed//fileu,EXIST=EX) !par IF(EX)THEN IF(IQDT.EQ.1)THEN WRITE(6,*)'*** ERROR: USE OF IREADK.LT.0 (kmtls.xxx) AND' X ,' IQDT=1 (smtls.xxx) IS INCONSISTENT.' WRITE(6,*)'RESET ONE OR OTHER VARIABLE DEPENDING ON' X ,' WHICH FILE YOU MEAN TO USE.' STOP 'INCONSISTENT INPUT: IREADK.LT.0 AND IQDT=1' ENDIF IQDT=2 ELSE WRITE(6,*)'***ERROR: REAL K-MX FILE kmtls.xxx' X ,' NOT FOUND' STOP'***ERROR: REAL K-MX FILE kmtls.xxx NOT FOUND' ENDIF ENDIF C IF(IRD0.LT.50) X OPEN(UNIT=IRD0,FILE=FILNAM//KAR//filec//filed//fileu, !par X FORM="UNFORMATTED",STATUS='OLD') 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 AUGER WIDTHS; INPUT AS AAUNITS THEN CONVERTED TO A.U. C ************ C IF(IAUGER.GT.0)THEN OPEN(4,FILE='AUGER',STATUS='UNKNOWN') C AAUNITS/Z^2 = 1 RYD, 0.5 AU, 2.06706E16 /SEC I.E. 1/HBAR DO I=1,NMLT AAUGER(I)=TZERO ENDDO READ(4,*)AAUNITS DO I=1,NMLT READ(4,*,END=39)II,AAA AAUGER(II)=AAA/(TWO*AAUNITS) ENDDO 39 CONTINUE CLOSE(4) 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/S-MATRICES C WILL 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(70+iam,*)iam,mxe,iee c do i=1,mxe c write(70+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,IQDT 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,IQDT 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 DO I=1,NLEV EINL(I)=ELEV(I)/AZSQ NCHJX(I)=(JLEV(I)+1)*(1-2*JPLEV(I)) ENDDO C 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 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, X 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, X 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, X 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, X 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 WRITE PRELIMINARY DATA TO OMEGDR FILE C IF(NDRMET.GT.0)THEN OPEN(UNIT=24,FILE="OMEGDR"//filec//filed//fileu, !par X FORM="FORMATTED",STATUS='UNKNOWN') WRITE(24,*) NZ,NELC WRITE(24,*) NLEV,MXE,NDRMET WRITE(24,*) (' 0',JLEV(I)+1,I=1,NLEV) WRITE(24,270) (EINL(I),I=1,NLEV) C INITIALIZE OMEGDR DO I=1,MXE DO J=1,NDRMET OMEGDR(J,I)=TZERO ENDDO ENDDO 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 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 IEK=IEK0 IF(IMODE.LT.0)IEK=IEK-1 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) PARAMETER (CONAS=6.476565E-8) !ALPHA**3/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,IQDT 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),ARDEC(NMLT),NEWAR,NTYP1 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 ARDEC(ILF)=TZERO WGHTZ=JLEV(ILF)+1 WGHTZ=WGHTZ*AZSQ 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 ARDEC(ILF)=ARDEC(ILF) !ARAD/Z**2 IN A.U. X +CONAS*(ELEV(ILF)-ELEV(ILI))**3*SS(N)/WGHTZ !ELEV IN RYD 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) PARAMETER (CONAS=6.476565E-8) !ALPHA**3/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,IQDT 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),ARDEC(NMLT),NEWAR,NTYP1 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 ARDEC(ILF)=TZERO WGHTZ=JLEV(ILF)+1 WGHTZ=WGHTZ*AZSQ 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 ARDEC(ILF)=ARDEC(ILF) !ARAD/Z**2 IN A.U. X +CONAS*(ELEV(ILF)-ELEV(ILI))**3*SS(N)/WGHTZ 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,IQDT 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,IREADK 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 ZK/SMTJK.DAT AND WRITE HEADER DATA TO FILE JBINJK C IF(IQDT.EQ.1) X OPEN(UNIT=16,FILE='smtjk.dat',STATUS='UNKNOWN', X FORM="UNFORMATTED") IF(IQDT.EQ.2) X OPEN(UNIT=16,FILE='zkmtjk.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 ZK/SMTIC.DAT AND WRITE HEADER DATA TO FILE JBINIC C IF(IQDT.EQ.1) X OPEN(UNIT=16,FILE='smtic.dat',STATUS='UNKNOWN', X FORM="UNFORMATTED") IF(IQDT.EQ.2) X OPEN(UNIT=16,FILE='zkmtic.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,IQDT 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,IREADK 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,IQDT 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,IREADK 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-Y) IMPLICIT COMPLEX*16 (Z) C C SUBROUTINE TO CALCULATE THE UNPHYSICAL K/S-MATRICES ON A NEW ENERGY C MESH FROM LINEAR INTERPOLATION OF THE UNPHYSICAL K/S-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,IQDT COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/ZRLS(MXRLS),ZRJ(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/ZRJ1(NOMCJ),ZRJ2(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) (ZRJ1(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) (ZRJ2(N),N=1,ICJMX) ELSE READ(18) (ZRJ1(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/S-MATRICES C IF(MOD(IEK-IEK0+1,2).EQ.0)THEN DO N=1,ICJMX ZSLOPE=(ZRJ2(N)-ZRJ1(N))/(EMESH(IEK)-EMESH(IEK-1)) ZRJ(N)=ZRJ1(N)+ZSLOPE*(EMSH(IEE)-EMESH(IEK-1)) END DO ELSE DO N=1,ICJMX ZSLOPE=(ZRJ1(N)-ZRJ2(N))/(EMESH(IEK)-EMESH(IEK-1)) ZRJ(N)=ZRJ2(N)+ZSLOPE*(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-Y) IMPLICIT COMPLEX*16 (Z) C C SUBROUTINE TO CALCULATE THE UNPHYSICAL K/S-MATRICES ON A NEW ENERGY C MESH FROM LINEAR INTERPOLATION OF THE UNPHYSICAL K/S-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,IQDT COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/ZRLS(MXRLS),ZRJ(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/ZRJ1(NOMCJ),ZRJ2(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) (ZRJ1(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) (ZRJ2(N),N=1,ICJMX) ELSE READ(18) (ZRJ1(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/S-MATRICES C IF(MOD(IEK-IEK0+1,2).EQ.0)THEN DO N=1,ICJMX ZSLOPE=(ZRJ2(N)-ZRJ1(N))/(EMESH(IEK)-EMESH(IEK-1)) ZRJ(N)=ZRJ1(N)+ZSLOPE*(EMSH(IEE)-EMESH(IEK-1)) END DO ELSE DO N=1,ICJMX ZSLOPE=(ZRJ1(N)-ZRJ2(N))/(EMESH(IEK)-EMESH(IEK-1)) ZRJ(N)=ZRJ2(N)+ZSLOPE*(EMSH(IEE)-EMESH(IEK-1)) END DO ENDIF C RETURN END C C ****************************************************************** C SUBROUTINE KMTRIC(IEK,IWJ0) IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C C DETERMINE THE K/S-MATRICES IN IC COUPLING FROM THE K/S-MATRICES C IN JK COUPLING C C INCLUDE 'PARAM' C PARAMETER (ZERO=(0.0,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,IQDT COMMON/BLOK5/ZRLS(MXRLS),ZRJ(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/ZTM(NOMCJ),ZA(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/ZRJ1(NOMCJ),ZRJ2(NOMCJ),IEK0 COMMON/HIGHPW/INOEXCH,JMNTWO,IREADK COMMON/NRBWRK/ZWORK(NMCJ),WORK(NMCJ) C IPOS(I,J,N)=J+N*(I-1)-(I*(I-1))/2 C C GENERATE K/S-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 ZWORK(ICCP)=ZERO 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 ZWORK(ICCP)=ZWORK(ICCP)+CVP*ZA(NN) ENDDO C NB ICCP CAN BE .GT. ICC NN=ICJM*(ICC-1)-NN+2*ICC DO ICCP=ICC+1,ICJM NN=NN+1 ZWORK(ICCP)=ZWORK(ICCP)+CVP*ZRJ(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) ZTM(N)=ZERO DO IVKP=1,NVEC(ILP) ITP=IVEC(IVKP,ILP) JVP=JVSV(ITP,INJP) IF(JVP.GT.0) THEN ICCP=ICHN(ITP,JVP,INLP,KVP) ZTM(N)=ZTM(N)+ZWORK(ICCP)*CVEC(IVKP,ILP) ENDIF ENDDO ENDDO ENDDO C C NOW STORE THE INTERMEDIATE COUPLING UNPHYSICAL K/S-MATRIX IN C RJ FOR THIS PARTIAL WAVE. C NMAX=N DO N=1,NMAX ZRJ(N)=ZTM(N) END DO C C WRITE THE UNPHYSICAL K/S-MATRIX IN IC COUPLING TO KMTIC C IF(IPRTIC.NE.0) THEN WRITE(12,300) FJTT(IWJ),JPTT(IWJ),ICJMX,NCHOJ(IWJ) 300 FORMAT(//46X,'K/S-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 ,(ZRJ(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/S-MATRIX TO FILE ZK/SMTIC.DAT C CASE IMODE.EQ.0.AND.JMNTWO.LT.0 C IF(IJBIN.NE.0)WRITE(16)(ZRJ(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 ZRJ1(N)=ZRJ(N) ENDDO ELSE N=0 DO N=1,NMAX ZRJ2(N)=ZRJ(N) ENDDO ENDIF ENDIF C RETURN END C C ****************************************************************** C SUBROUTINE KMTRJK(IEK,IWJ0) IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C C DETERMINE THE UNPHYSICAL K/S-MATRICES IN JK COUPLING FROM THE C THE UNPHYSICAL K/S-MATRICES IN LS COUPLING C C INCLUDE 'PARAM' C PARAMETER (ZERO=(0.0,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,IQDT COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/ZRLS(MXRLS),ZRJ(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/ZTM(NOMCJ),ZA(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,IREADK COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/INTER/ZRJ1(NOMCJ),ZRJ2(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/S-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 ZRJ(N)=ZERO ENDDO IF(ITCC.NE.0)THEN DO N=1,NMAX ZA(N)=ZERO 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) ZRJ(N)=ZRJ(N)+CLS(ICJ,IL)*CLS(ICJP,IL)*ZRLS(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) ZA(NN)=ZRJ(N) ENDDO ENDDO ENDIF C C WRITE THE UNPHYSICAL K/S-MATRIX IN JK COUPLING TO ZK/SMTJK C IF(IPRTJK.NE.0) THEN WRITE(11,300) FJTT(IWJ),JPTT(IWJ),ICJM,NCHOJ(IWJ) 300 FORMAT(//46X,'K/S-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 (ZRJ(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/S-MATRIX TO FILE ZK/SMTJK.DAT C CASE IMODE.EQ.0.AND.JMNTWO.LT.0 C IF(ITCC.EQ.0.AND.IJBIN.NE.0)WRITE(16)(ZRJ(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 ZRJ1(N)=ZRJ(N) ENDDO ELSE N=0 DO N=1,NMAX ZRJ2(N)=ZRJ(N) ENDDO ENDIF ENDIF C RETURN END C C ****************************************************************** C SUBROUTINE MQDTS(IEE,IWJ,NLVOP,NCHOP) IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C C DETERMINE THE PHYSICAL S-MATRIX FROM THE UNPHYSICAL S-MATRIX C C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (ZERO=(0.0,0.0)) PARAMETER (CONAS=6.476565E-8) !ALPHA**3/6 PARAMETER (BIG=150.) C COMMON/AUGER/AAUGER(NMLT),IAUGER 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,IQDT COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/ZRLS(MXRLS),ZRJ(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/ZTM(NOMCJ),ZA(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/DIPOLE/SSLS(NOMTM),SS(NOMLT),ARDEC(NMLT),NEWAR,NTYP1 COMMON/NRBWRK/ZWORK(NMCJ),WORK(NMCJ) COMMON/NRBQDT/QNMAX,FNUMAX C TPI=TWO*ACOS(-ONE) C IF(ITCC.EQ.0)NCC=NCHPJ(IWJ)-NCHOP IF(ITCC.NE.0)NCC=NCHPJX(IWJ)-NCHOP C IF(NCC.LT.1) GO TO 100 C C DETERMINE RADIATIVE WIDTH C IF(NTYP1*NEWAR.GT.0)THEN EC=AZSQ*EMSH(IEE) DO J=NLVOP+1,NLEV K=((J-1)*(J-2*IONE))/2 ARDEC(J)=TZERO WGHTZ=JLEV(J)+1 WGHTZ=WGHTZ*AZSQ DO I=1,J-1 K=K+1 DE=ELEV(J)-ELEV(I) IF(EC-DE.LE.ELEV(1))THEN ARDEC(J)=ARDEC(J)+CONAS*DE**3*SS(K)/WGHTZ ENDIF ENDDO ENDDO ENDIF C C STORE SCC FIRST C NOMAX=(NCHOP*(NCHOP+1))/2 NCMAX=(NCC*(NCC+1))/2 NN=NOMAX+NCC*NCHOP DO N=1,NCMAX NN=NN+1 ZA(N)=ZRJ(NN) !PACK LOWER HALF END DO C C ADD EXP(-2*PI*NU*I) DOWN THE DIAGONAL C N=1 ICT=NCHOP DO ICC=NCC,1,-1 ICT=ICT+1 IF(NTYP1.EQ.0.AND.IAUGER.LE.0.OR.FNU(ICT).GT.FNUMAX)THEN Z=DCMPLX(TZERO,-TPI*FNU(ICT)) ZEXP=EXP(Z) ZA(N)=ZA(N)-ZEXP ELSE IF(ITCC.EQ.0)ILV=ILVJ(ICT) IF(ITCC.NE.0)ILV=ILVJX(ICT) TR0=ELEV(ILV)/AZSQ-EMSH(IEE) TI0=TZERO IF(NTYP1.GT.0)TI0=-ARDEC(ILV) IF(IAUGER.GT.0)TI0=TI0-AAUGER(ILV) ZFNU=DCMPLX(ONE,TZERO)/SQRT(DCMPLX(TR0,TI0)) Z=DCMPLX(TZERO,-TPI)*ZFNU TR0=DBLE(Z) TR0=MIN(TR0,BIG) TI0=DIMAG(Z) Z=DCMPLX(TR0,TI0) ZEXP=EXP(Z) ZA(N)=ZA(N)-ZEXP ENDIF N=N+ICC END DO C C INITIALIZE SCO C N=0 NN=0 DO IC=NCHOP,1,-1 NN=NN+IC DO ICC=1,NCC N=N+1 NN=NN+1 CK ICT=ICC+NCHOP CK IF(FNU(ICT).LE.FNUMAX) THEN ZTM(N)=ZRJ(NN) CK ELSE CK ZTM(N)=ZERO CK ENDIF ENDDO ENDDO C CSTRTNBL CALL ZLUS(ZA,-NMCJ,NCC,WORK,IERR) IF (IERR.NE.0) THEN WRITE(6,600) STOP 'ERROR IN ZLUS' END IF CALL ZLUBS(ZA,ZTM,NCC,NCHOP,IERR) IF (IERR.NE.0) THEN WRITE(6,601) STOP 'ERROR IN ZLUBS' END IF CENDNBL C CSTRTBL CBL CALL ZSPTRF('L',NCC,ZA,IPIVOT,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,602) INFO CBL STOP 'ERROR IN BLAS ROUTINE ZSPTRF' CBL ENDIF CBL CALL ZSPTRS('L',NCC,NCHOP,ZA,IPIVOT,ZTM,NCC,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,603) INFO CBL STOP 'ERROR IN BLAS ROUTINE ZSPTRS' CBL ENDIF CENDBL C C C INITIALIZE SOC^T=SCO C N=0 NN=0 DO IC=NCHOP,1,-1 NN=NN+IC DO ICC=1,NCC N=N+1 NN=NN+1 CK ICT=ICC+NCHOP CK IF(FNU(ICT).LE.FNUMAX) THEN ZA(N)=ZRJ(NN) CK ELSE CK ZA(N)=ZERO CK ENDIF ENDDO ENDDO C C PACK SOO C N=NCHOP NN=NCHOP DO ICP=2,NCHOP NN=NN+NCC DO IC=ICP,NCHOP N=N+1 NN=NN+1 ZRJ(N)=ZRJ(NN) ENDDO ENDDO C C NOW FINISH UP THE CALCULATION OF THE PHYSICAL S-MATRIX (IN ZRJ) C 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 ZRJ(N)=ZRJ(N)-ZA(NI)*ZTM(NJ) ENDDO ENDDO NN=NN+NCC ENDDO C 100 RETURN C 600 FORMAT(' SR.MQDTS: ZLUS RETURNED WITH INFO =',I6) 601 FORMAT(' SR.MQDTS: ZLUBS RETURNED WITH INFO =',I6) CB602 FORMAT(//10X,10('*'),' SR. MQDTS: ZSPTRF RETURNED WITH INFO =',I6) CB603 FORMAT(//10X,10('*'),' SR. MQDTS: ZSPTRS RETURNED WITH INFO =',I6) 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-Y) IMPLICIT COMPLEX*16 (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 (ZERO=(0.0,0.0)) PARAMETER (TINY=1.0E-6) 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,IQDT COMMON/BLOK2/LT(NMTM),LST(NMTM),LPT(NMTM),NJV(NMTM), 1 FLT(NMTM),FST(NMTM),ENT(NMTM),FJT(NMTM,NMLV) COMMON/BLOK5/ZRLS(MXRLS),ZRJ(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/ZTM(NOMCJ),ZA(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),ARDEC(NMLT),NEWAR,NTYP1 COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/MEMORY/OMEM(MXOM+1),MPOS(0:NEMXF),ITMAX,JTMAX COMMON/NRBDR/PDR(NMCJ),OMEGDR(NMET,NEMXF),NDRMET COMMON/NRBQDT/QNMAX,FNUMAX 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 DIMENSION TM(NOMCJ) EQUIVALENCE (ZTM(1),TM(1)) 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 NCHOP=NCHOJ(IWJ) NCHOP1=NCHOP+1 DO I=NCHOP1,NCHPJX(IWJ) ILV=ILVJX(I) EE=ELEV(ILV)/AZSQ-EMSH(IEE) FNU(I)=ONE/SQRT(EE) ENDDO C C DETERMINE WEAKLY CLOSED CHANNELS FOR GAILITIS C IF(QNMAX.GT.TZERO)THEN ILVP=ILVJX(NCHOP1) DO I=NCHOP1,NCHPJX(IWJ) ILV=ILVJX(I) IF(EINL(ILV).GT.EINL(ILVP)+TINY)GO TO 12 IF(FNU(I).GT.QNMAX)THEN NCHOP=I ELSE GO TO 12 ENDIF ENDDO ENDIF C C INITIALIZE PARTIAL OMEGAS C 12 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) C C REWIND ANY SCRATCH FILES IF APPROPRIATE C 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 IF(NCHOJ(IWJ).LE.IONE) GO TO 270 C IF(IQDT.EQ.2)THEN C C DETERMINE THE PHYSICAL K MATRIX C CALL ZMQDTK(IEE,IWJ,NLVOP,NCHOP) C C CALCULATE THE S-MATRIX C CALL ZSMTRX(NCHOP) C ELSEIF(IQDT.EQ.1)THEN C C DETERMINE THE PHYSICAL S MATRIX C CALL MQDTS(IEE,IWJ,NLVOP,NCHOP) C ELSE WRITE(6,*)'***OMGIC: IQDT=',IQDT,' NOT SET CORRECTLY' STOP'***OMGIC: IQDT NOT SET CORRECTLY' ENDIF C C APPLY GAILITIS AVERAGE (INC DR) C IF(NCHOP.GT.NCHOJ(IWJ))THEN C CALL SQDTS(IEE,IWJ,NLVOP,NCHOP) C ELSE C C DETERMINE DR PROBABILITY AND/OR COLLISION STRENGTH (NON-GAILITIS CASE) C IF(NDRMET.GT.0)THEN DO I=1,NCHOJ(IWJ) PDR(I)=ONE ENDDO N=0 DO J=1,NCHOJ(IWJ) DO I=J,NCHOJ(IWJ) N=N+1 PP=ZRJ(N)*CONJG(ZRJ(N)) PDR(I)=PDR(I)-PP PDR(J)=PDR(J)-PP TM(N)=PP IF(I.EQ.J)THEN PDR(J)=PDR(J)+PP TM(N)=TM(N)+ONE-TWO*DBLE(ZRJ(N)) !NOT -ONE ENDIF ENDDO ENDDO ELSE N=1 DO J=NCHOJ(IWJ),1,-1 ZRJ(N)=ZRJ(N)-ONE !-T N=N+J END DO NOMAX=(NCHOJ(IWJ)*(NCHOJ(IWJ)+1))/2 DO N=1,NOMAX TM(N)=ZRJ(N)*CONJG(ZRJ(N)) !T**2 ENDDO N=1 DO J=NCHOJ(IWJ),1,-1 ZRJ(N)=ZRJ(N)+ONE !S N=N+J END DO ENDIF ENDIF 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 DETERMINE DR COLLISION STRENGTHS C IF(NDRMET.GT.0)THEN IP=MIN(NDRMET,NLVOP) DO ILI=1,IP JIMX=NCHJX(ILI) DO JI=1,JIMX ICI=ICHJX(ILI,JI) OMEGDR(ILI,IEE)=OMEGDR(ILI,IEE)+COEFF*PDR(ICI) ENDDO ENDDO ENDIF 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 DR COLLISION STRENGTHS TO FILE OMEGDR C 275 IF(NDRMET.GT.0)THEN WRITE(24,380) ENGMSH,(OMEGDR(N,IEE),N=1,NDRMET) ENDIF C C WRITE TOTAL COLLISION STRENGTHS TO FILE OMEGA, AND OMEGANOTOP (ITOP.NE.0) C 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-Y) IMPLICIT COMPLEX*16 (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 (ZERO=(0.0,0.0)) PARAMETER (TINY=1.0E-6) 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,IQDT 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/BLOK5/ZRLS(MXRLS),ZRJ(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/ZTM(NOMCJ),ZA(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) COMMON/BLOK10/JLEV(NMLT),JPLEV(NMLT),ELEV(NMLT),EINL(NMLT),NLEV COMMON/DIPOLE/SSLS(NOMTM),SS(NOMLT),ARDEC(NMLT),NEWAR,NTYP1 COMMON/ENG/EMESH(0:NEMXC),EMSH(0:NEMXF) COMMON/OMG/OMGPW(NOMLT),OMEGA(NOMLT) COMMON/NRBDR/PDR(NMCJ),OMEGDR(NMET,NEMXF),NDRMET COMMON/NRBQDT/QNMAX,FNUMAX 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 DIMENSION TM(NOMCJ) EQUIVALENCE (ZTM(1),TM(1)) 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 NCHOP=NCHOJ(IWJ) NCHOP1=NCHOP+1 DO I=NCHOP1,NCHPJ(IWJ) ILV=ILVJ(I) EE=ELEV(ILV)/AZSQ-EMSH(IEE) FNU(I)=ONE/SQRT(EE) ENDDO C C DETERMINE WEAKLY CLOSED CHANNELS FOR GAILITIS C IF(QNMAX.GT.TZERO)THEN ILVP=ILVJ(NCHOP1) DO I=NCHOP1,NCHPJ(IWJ) ILV=ILVJ(I) IF(EINL(ILV).GT.EINL(ILVP)+TINY)GO TO 12 IF(FNU(I).GT.QNMAX)THEN NCHOP=I ELSE GO TO 12 ENDIF ENDDO ENDIF C C INITIALIZE PARTIAL OMEGAS C 12 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) C C REWIND ANY SCRATCH FILES IF APPROPRIATE C 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 IF(NCHOJ(IWJ).LE.IONE) GO TO 270 C IF(IQDT.EQ.2)THEN C C DETERMINE THE PHYSICAL K MATRIX C CALL ZMQDTK(IEE,IWJ,NLVOP,NCHOP) C C CALCULATE THE S-MATRIX C CALL ZSMTRX(NCHOP) C ELSEIF(IQDT.EQ.1)THEN C C DETERMINE THE PHYSICAL S MATRIX C CALL MQDTS(IEE,IWJ,NLVOP,NCHOP) C ELSE WRITE(6,*)'***OMGJK: IQDT=',IQDT,' NOT SET CORRECTLY' STOP'***OMGJK: IQDT NOT SET CORRECTLY' ENDIF C C APPLY GAILITIS AVERAGE (INC DR) C IF(NCHOP.GT.NCHOJ(IWJ))THEN C CALL SQDTS(IEE,IWJ,NLVOP,NCHOP) C ELSE C C DETERMINE DR PROBABILITY AND/OR COLLISION STRENGTH (NON-GAILITIS CASE) C IF(NDRMET.GT.0)THEN DO I=1,NCHOJ(IWJ) PDR(I)=ONE ENDDO N=0 DO J=1,NCHOJ(IWJ) DO I=J,NCHOJ(IWJ) N=N+1 PP=ZRJ(N)*CONJG(ZRJ(N)) PDR(I)=PDR(I)-PP PDR(J)=PDR(J)-PP TM(N)=PP IF(I.EQ.J)THEN PDR(J)=PDR(J)+PP TM(N)=TM(N)+ONE-TWO*DBLE(ZRJ(N)) !NOT -ONE ENDIF ENDDO ENDDO ELSE N=1 DO J=NCHOJ(IWJ),1,-1 ZRJ(N)=ZRJ(N)-ONE !-T N=N+J END DO NOMAX=(NCHOJ(IWJ)*(NCHOJ(IWJ)+1))/2 DO N=1,NOMAX TM(N)=ZRJ(N)*CONJG(ZRJ(N)) !T**2 ENDDO N=1 DO J=NCHOJ(IWJ),1,-1 ZRJ(N)=ZRJ(N)+ONE !S N=N+J END DO ENDIF ENDIF 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 DETERMINE DR COLLISION STRENGTHS C 100 IF(NDRMET.GT.0)THEN ILI=0 DO I=1,NTLS IF(EC.LT.ENT(I))GO TO 101 NLVI=NJV(I) DO II=1,NLVI ILI=ILI+1 IF(ILI.GT.NDRMET)GO TO 101 JIMX=NCHJ(I,II) DO JI=1,JIMX ICI=ICHJ(I,II,JI) OMEGDR(ILI,IEE)=OMEGDR(ILI,IEE)+COEFF*PDR(ICI) ENDDO ENDDO ENDDO ENDIF C C WRITE PARTIAL OMEGAS TO SCRATCH FOR NON-DIPOLE TOP C 101 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 DR COLLISION STRENGTHS TO FILE OMEGDR C 275 IF(NDRMET.GT.0)THEN WRITE(24,380) ENGMSH,(OMEGDR(N,IEE),N=1,NDRMET) ENDIF C C WRITE TOTAL COLLISION STRENGTHS TO FILE OMEGA, AND OMEGANOTOP (ITOP.NE.0) C 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,IQDT 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,IQDT 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 ZK/SMTIC.DAT C IF(IQDT.EQ.1) XOPEN(UNIT=18,FILE='smtic.dat',STATUS='OLD',FORM="UNFORMATTED") IF(IQDT.EQ.2) XOPEN(UNIT=18,FILE='zkmtic.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,IQDT 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,IQDT 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 ZK/SMTJK.DAT C IF(IQDT.EQ.1) XOPEN(UNIT=18,FILE='smtjk.dat',STATUS='OLD',FORM="UNFORMATTED") IF(IQDT.EQ.2) XOPEN(UNIT=18,FILE='zkmtjk.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-Y) IMPLICIT COMPLEX*16 (Z) C C SUBROUTINE TO READ THE K-MATRIX FILE FROM z/k/smtls.001 etc FROM C THE 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 ZRLS(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) PARAMETER (ZERO=(0.0,0.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,IQDT 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/ZRLS(MXRLS),ZRJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK7/ZTM(NOMCJ),ZA(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,IREADK C c **** parallel **** character*1 filec,filed,fileu common/parainfo/iam,nproc include 'mpif.h' common/blockstop/istop c **** parallel **** c DIMENSION TEMP(NOMCL) EQUIVALENCE (TEMP(1),ZTM(1)) C SAVE IWJ,IRD DATA IWJ/0/,IRD/0/ C NPOS(I,J)=(J*(J-1))/2+I C C READ THE UNPHYSICAL K/S-MATRICES FOR THE LSPI PARTIAL WAVES C REQUIRED FOR THE INPUT JP SYMMETRY. C DO N=1,MXRLS ZRLS(N)=ZERO 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 IF(INOEXCH.GT.0)THEN !COMPLEX ZK OR S-MATRIX READ(IRD)(ZA(IC),IC=1,NRD) !N.B. UPPER HALF LS FROM STGF... ELSE !IF(INOEXCH.LT.0)THEN REAL K-MATRIX READ(IRD)(TEMP(N),N=1,NRD) DO N=1,NRD ZA(N)=DCMPLX(TEMP(N),TZERO) ENDDO ENDIF 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) ZRLS(N)=ZA(NN) END IF END DO ENDIF END DO END IF ENDDO ENDDO C C PRINT THE K/S-MATRIX TO FILE K/SMTLS.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(//42X,'K/S-MATRIX FOR L =',I3,' 2S+1 =',I3, 1 ' PARITY =',I3/44X,43('*')) DO I=1,NCHT(IW) WRITE(10,310) I,ITM(IW,I),LC(IW,I), 1 (ZRLS(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-Y) IMPLICIT COMPLEX*16 (Z) C C READ THE K-MATRICES IN LS COUPLING FROM z/k/smtls.001,002... FOR THE C IEKTH ENERGY NECESSARY TO GENERATE JP SYMMETRY IWJ0 AND STORE THEM IN C C ZRLS(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,IQDT 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/ZRLS(MXRLS),ZRJ(NOMCJ),LSPOS(0:NMLS) COMMON/BLOK7/ZTM(NOMCJ),ZA(NOMCJ),FNU(NMCJ),IPIVOT(NMCJ) 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/HIGHPW/INOEXCH,JMNTWO,IREADK COMMON/HOLDLS/ZRLSHLD(NMLSE*MXRLS+1) DIMENSION TEMP(NOMCL) EQUIVALENCE (TEMP(1),ZTM(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 ZRLS(N)=ZRLSHLD(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) IF(IREADK.GE.0)THEN READ(IRD)(ZRLS(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=N2+NRD READ(IRD)(ZRLSHLD(N),N=N1,N2) ENDDO ENDIF ELSE !READ REAL K-MX FROM kmtls.dat READ(IRD)(TEMP(N),N=1,NRD) DO N=1,NRD ZRLS(N+LSPOS(IL-1))=DCMPLX(TEMP(N),TZERO) ENDDO IF(IRD0.LT.50)THEN !READ EXTRA ENERGIES TO SAVE ON I/O N2=NMLSE*LSPOS(IL-1) DO M=1,MM READ(IRD)(TEMP(N),N=1,NRD) DO N=1,NRD ZRLSHLD(N2+N)=DCMPLX(TEMP(N),TZERO) ENDDO N2=N2+NRD ENDDO ENDIF ENDIF ENDDO C C PRINT THE K/S-MATRIX TO FILE K/SMTLS.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(//42X,'K/S-MATRIX FOR L =',I3,' 2S+1 =',I3, 1 ' PARITY =',I3,/44X,43('*')) DO I=1,NCHT(IW) WRITE(10,310) I,ITM(IW,I),LC(IW,I), 1 (ZRLS(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,IQDT 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,IREADK 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,IQDT 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,IREADK 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,IQDT 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,IQDT 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 SQDTS(IEE,IWJ,NLVOP,NCHOP) IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C C DETERMINE THE PHYSICAL S-MATRIX FROM THE UNPHYSICAL S-MATRIX. C THE CLOSED CHANNELS ARE ONLY THOSE ATTACHED TO THE LOWEST C CLOSED TARGET LEVEL I.E. NCC<