! ! program pstgb.f ! ! Parallel code derived from the serial code stgb.f ! ! References : ! ! 1. Use of the R matrix method for bound-state calculations. ! I. General theory ! M J Seaton J. Phys. B: At. Mol. Phys. 18 No 11 (14 June 1985) 2111-2131 ! ! ! 2. The vicinity of an R-matrix pole ! P G Burke and M J Seaton ! J. Phys. B: At. Mol. Phys. 17 No 20 (28 October 1984) L683-L687 ! ! ! Notes : the code is standalone (with PARAM), though using the BLAS ! can give over a factor of 2 speed up if accelerated ! Lapack/ATLAS/ACML libraries are implemented on your machine. ! ! use : mpirun -np XX pstgb.x (where XX is the number of processors) ! !**************************************************************************** module comm_interface implicit none include "mpif.h" public comm_init ! Initialize MPI public comm_barrier ! MPI barrier public comm_finalize ! Terminate MPI integer, public :: iam integer, public :: nproc SAVE private integer :: mpicom CONTAINS !----------------------------------------------------------------------------- subroutine comm_init() implicit none integer :: ier mpicom = MPI_COMM_WORLD call mpi_init(ier) call mpi_comm_rank(mpicom, iam, ier) call mpi_comm_size(mpicom, nproc, ier) return end subroutine comm_init !----------------------------------------------------------------------------- subroutine comm_barrier() implicit none integer :: ier call mpi_barrier(mpicom, ier) return end subroutine comm_barrier !----------------------------------------------------------------------------- subroutine comm_finalize() implicit none integer :: ier call mpi_finalize(ier) return end subroutine comm_finalize !----------------------------------------------------------------------------- end module comm_interface c PSTGB.f :: Parallel version of stgb. 07 OCT 20 c C.P. Ballance, D.C. Griffin, N.R. Badnell c v1.4/2.20 c C*********************************************************************** C N. R. BADNELL UoS V2.21 29/06/21 C C PROGRAM STGB C C BOUND-STATE CALCULATIONS C C RECENT MODIFICATIONS C NAMELIST I/O C DEEPLY CLOSED CHANNELS ZEROED-OUT C EXTENDED TO HANDLE NEUTRAL ATOMS C DARC DSTGH.DAT (AND H.DAT VIA DARC INTERFACE) C C C PROGRAM MAIN c use comm_interface, only : iam,nproc,comm_init,comm_barrier, A comm_finalize C IMPLICIT REAL*8 (A-H,O-Z) integer,allocatable :: ist(:),ilt(:),ipt(:) C INCLUDE 'PARAM' C C SUN REAL*4 TARRY(2) C C PARAMETER (N1TAR=10*MZTAR) C LOGICAL EX C CHARACTER*9 ID CHARACTER*4 ZZNN CHARACTER*2 IDNUM C COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CVECT/AAA(MZCHF,MZCHF),BBB(MZCHF,MZCHF),CCC(MZCHF,MZCHF), 1 P(MZCHF,MZCHF),Q(MZCHF,MZCHF),XVECT(MZCHF) COMMON/CBODE/WBODE(MZPTS) COMMON/COULSC/FS(MZCHF,MZPTS),FSP(MZCHF),FC(MZCHF,MZPTS) 1 ,FCP(MZCHF) COMMON/CPOT/LAMP(MZCHF,MZCHF),BW(MZCHF,MZCHF) COMMON/CALP/ASS(MZCHF,MZCHF),ASC(MZCHF,MZCHF),ACS(MZCHF,MZCHF) 1 ,ACC(MZCHF,MZCHF) COMMON/CNTRL/IPRINT,IRAD,IPERT COMMON/CITER/ECH1,XSTORE(MZEST),K0,KSTORE(MZEST) COMMON/CDEGEN/NASTD,NASTR,NLEV(MZTAR),NCNATR(MZTAR),ENATR(MZTAR) COMMON/CID/WT(MZCHF,MZEST),FNST(MZCHF,MZEST),ISLP,IDNST,IDNST1, 1 IDSLP(MZTAR),IDSLP1(N1TAR),NCORE(MZLP1) COMMON/NRBWRD/IWORD,KFLAG COMMON/NRBZED/TZED C DIMENSION XITER(5) DIMENSION AVECT(MZMNP+MZCHF) DIMENSION PB(MZCHF),PBP(MZCHF),PBPP(MZCHF) DIMENSION BVECT(MZCHF) DIMENSION MSLP(MZSLP), 1 XMINST(MZEST),XMAXST(MZEST),BDXST(MZEST), 2 MXEST(MZEST),ERDST(MZSLP,MZEST),EST(MZEST) DIMENSION IDSAT(MZTAR),IDLAT(MZTAR),IDPAT(MZTAR), 1 IDSAT1(N1TAR),IDLAT1(N1TAR),IDPAT1(N1TAR) DIMENSION ID(MZEST),IDDD(MZEST) DIMENSION IDNUM(0:99) DIMENSION ITARG1(MZSLP) C NAMELIST/STGB/IPRINT,IRAD,IPERT,AC,RONE,IOPT2,RTWOIN,IDNST,IDNST1 X,IORDER C DATA NZERO,NONE,NTWO/0,1,2/ CHARACTER NAME*3,NUM(0:9)*1 CHARACTER*7 filex DATA NUM /'0','1','2','3','4','5','6','7','8','9'/ DATA IDNUM/'00','01','02','03','04','05','06','07','08','09', 1 '10','11','12','13','14','15','16','17','18','19', 2 '20','21','22','23','24','25','26','27','28','29', 3 '30','31','32','33','34','35','36','37','38','39', 4 '40','41','42','43','44','45','46','47','48','49', 5 '50','51','52','53','54','55','56','57','58','59', 6 '60','61','62','63','64','65','66','67','68','69', 7 '70','71','72','73','74','75','76','77','78','79', 8 '80','81','82','83','84','85','86','87','88','89', 9 '90','91','92','93','94','95','96','97','98','99'/ C LOGICAL WARN COMMON/CWARN/WARN WARN=.TRUE. XNUMAX=0. KFLAG=0 C C CALL BLOCK DATA AS SUBROUTINE TO AVOID LINKAGE PROBLEMS C C call comm_init() C call cpu_time(timei) C CALL BLOCK C C I/O UNITS C ********* C C 5 FOR PROGRAM CONTROL C 6 FOR PRINTED OUTPUT C 7 FOR ENERGY LEVELS. FILE ELEV. C 8 FOR DIPOLE DATA INDEPENDENT OF SLPI. FILE B00. C 9 FOR DIPOLE DATA DEPENDENT ON SLPI. FILES B01 ETC. C 10 FOR R-MATRIX DATA. FILE H C 11 FOR LEVEL IDENTIFICATIONS. FILE ID. C 12 FOR LEVEL IDENTIFICATIONS. FILE E. C 13 FOR LEVEL IDENTIFICATIONS. FILE T. C 14 RESERVED C OPEN(5,FILE='dstgb') i1=(iam+1)/10 i2=((iam+1)-10*((iam+1)/10)) filex='routb'//NUM(i1)//NUM(i2) OPEN(6,FILE=filex,STATUS='UNKNOWN') C INQUIRE(FILE='H.DAT',EXIST=EX) IF(EX)THEN IWORD=1 OPEN(10,FILE='H.DAT',FORM='UNFORMATTED',STATUS='OLD') ELSE INQUIRE(FILE='DSTGH.DAT',EXIST=EX) IF(EX)THEN IWORD=2 OPEN(10,FILE='DSTGH.DAT',FORM='UNFORMATTED',STATUS='OLD') ELSE WRITE(6,*)'***ERROR: NO SUITABLE H.DAT FILE FOUND!' STOP '***ERROR: NO SUITABLE H.DAT FILE FOUND!' ENDIF ENDIF C C C WRITE NEWS C ********** C WRITE(6,6000) WRITE(6,6002)MZCHF,MZTAR,MZLMX,MZLP1,MZPTS,MZEST,MZMNP,MZSLP X ,MZTET,MZNRG C C C READ DATA FROM UNIT 5 C ********************* C C IPRINT = 0 FOR MINIMUM PRINT C IPRINT = 3 FOR MAXIMUM PRINT C IRAD = 0 FOR NO RADIATIVE DATA ON UNITS 8 AND 9 C IRAD = 1 FOR RADIATIVE DATA ON UNITS 8 AND 9 C IRAD = 2 FOR DATA ON UNITS 8 AND 9 AND NO WRITES TO 7 C RTWOIN = OPTIONAL INPUT OF RTWO, CASE IRAD.GT.0 C AS OPPOSED TO DEFAULT USE OF ENTIRE MESH. C IPERT = 0 FOR NO PERTURBATION C IPERT = 1 FOR FIRST-ORDER PERTURBATION C AC FOR ACCURACY REQUIRED C RONE = DEBUG PARAMETER (NORMALLY TAKE RONE = 1.) C IOPT2 = 1 TO SCAN FOR BOUND STATES C IOPT2 = 2 TO READ ESTIMATES OF BOUND-STATE ENERGIES C IOPT2 = -1 OR -2 AS ABOVE FOLLOWED BY C NASTD C (NLEV(N),N=1,NASTD) C IS, IL, IP FOLLOWED BY C XMIN, XMAX, BDX FOR IOPT2 = 1 OR C MXE AND (E, M=1,MXE) FOR IOPT2 = 2 C -1, -1, -1 FOR TERMINATION C IORDER = 1 PROCESS SYMMETRIES IN ORDER LISTED BY IOPT2 (DEFAULT) C = 0 PROCESS SYMMETRIES IN THE ORDER THEY APPEAR IN H.DAT C (THIS IS THE ORIGINAL AND IS SLIGHTLY FASTER BUT IT IS C MORE CONVENIENT FOR THE DAMPING CODES TO BE ABLE TO C LIST THE STGB SYMMETRIES IN ENERGY ORDER.) C C INITIALIZE NAMELIST C IPRINT=0 IRAD=0 IPERT=0 AC=1.0E-5 RONE=1.0 IOPT2=1 RTWOIN=-1.0 IDNST=0 IDNST1=-1 IORDER=1 C READ(5,STGB) C WRITE(6,600)IPRINT,IRAD,IPERT,AC,RONE,RTWOIN NASTD=0 IF(IOPT2.LT.0)THEN IOPT2=ABS(IOPT2) READ(5,*)NASTD IF(NASTD.GT.MZTAR)THEN WRITE(6,605)NASTD,MZTAR STOP END IF READ(5,*)(NLEV(N),N=1,NASTD) WRITE(6,606)NASTD,(NLEV(N),N=1,NASTD) END IF IF(IOPT2.EQ.1)WRITE(6,601) IF(IOPT2.EQ.2)WRITE(6,602) C 500 READ(5,*,END=502)IS,IL,IP C IF(IL.EQ.-1)GOTO 502 C KSLP=KSLP+1 KKCOUNT=0 IFLAGAA=0 c**** allocate(IST(0:nproc-1),ILT(0:nproc-1),IPT(0:nproc-1)) c**** DO JJJ=0,nproc-1 READ(5,*,END=502)IS,IL,IP READ(5,*)XMINST(1),XMAXST(1),BDXST(1) if((iam).eq.JJJ)then KSLP=1 MSLP(KSLP)=10000*IS+100*IL+IP C write(0,'(I5,3F9.5,I5)')MSLP(KSLP),XMINST(1),XMAXST(1),BDXST(1) C X ,iam endif if(iam.eq.0)then IST(JJJ)=IS ILT(JJJ)=IL IPT(JJJ)=IP endif KKCOUNT=KKCOUNT+1 ENDDO if(KKCOUNT.ne.nproc)then write(6,*)'number of symmetries must equal procs!' call comm_barrier() call comm_finalize() stop'number of symmetries must equal procs!' endif if(IOPT2.ne.1)then write(6,*)' IOPT2 eq 1 , only choice at the moment' call comm_barrier() call comm_finalize() stop' IOPT2 eq 1 , only choice at the moment' endif C IF(IOPT2.EQ.1)THEN C READ(5,*)XMINST(KSLP),XMAXST(KSLP),BDXST(KSLP) C ELSE C READ(5,*)MXEST(KSLP) C DO 501 M=1,MXEST(KSLP) C 501 READ(5,*)ERDST(KSLP,M) C ENDIF C GOTO 500 502 CONTINUE WRITE(6,603)(MSLP(K),K=1,KSLP) WRITE(6,604) CALL ACSUB C C READ R-MATRIX DATA INDEPENDENT OF SLPI C ************************************** C CALL READ1 C IF(NZED.EQ.NELC)TZED=0. IF(NZED.NE.NELC)TZED=1. C CALL SCALE1 C INITIAL WRITES TO UNITS 7 AND 8 IF(IRAD.LT.2)THEN i1=(iam+1)/10 i2=((iam+1)-10*((iam+1)/10)) filex='ELEV'//NUM(i1)//NUM(i2) OPEN(7,FILE=filex,FORM='FORMATTED',STATUS='UNKNOWN') REWIND(7) WRITE(7,700)IPRINT,IRAD,IPERT,NZED,NELC WRITE(7,701)AC WRITE(7,702)RONE WRITE(7,703) C C **** CODING TO READ DATA FOR LEVEL IDENTIFICATION C IF(IDNST1.GE.0)THEN filex='ID'//NUM(i1)//NUM(i2) OPEN (11,FILE=filex,STATUS='NEW') ZZNN=IDNUM(NZED)//IDNUM(NELC) WRITE(11,1110)ZZNN,KSLP C C IF IDNST = 0 TAKE THE STATES IN THE ORDER READ FROM FILE H C IF (IDNST .NE. 0) THEN IF (IDNST .NE. NASTR) THEN WRITE(6,'('' ERROR IN DATA. NUMBER OF ID TARGETS '', 1 ''DIFFERENT FROM NAST READ FROM H FILE.''/ 2 '' IDNST = '',I3,'' NAST = '',I3)')IDNST,NASTR CALL EXIT(0) END IF C READ(5,*)(IDSAT(I),IDLAT(I),IDPAT(I),I=1,IDNST) ! READ5******** C END IF C IF (IDNST1 .GT. 10*MZTAR) THEN WRITE(6,'('' NUMBER OF N+1 ELECTRON ID STATES EXCEEDS '', 1 ''10 TIMES MZTAR. INCREASE MZTAR TO'')')IDNST1/10+1 CALL EXIT(0) END IF C IF (IDNST1 .GT. 0) THEN C READ(5,*)(IDSAT1(I),IDLAT1(I),IDPAT1(I),I=1,IDNST1) !READ5*** C END IF C DO 1100 I=1,IDNST IDSLP(I)=10000*IDSAT(I)+100*IDLAT(I)+IDPAT(I) 1100 CONTINUE DO 1101 I=1,IDNST1 IDSLP1(I)=10000*IDSAT1(I)+100*IDLAT1(I)+IDPAT1(I) 1101 CONTINUE C READ(5,*)LRANG1,(NCORE(L),L=1,LRANG1) !READ5***** C DO 1102 LP1=LRANG1+1,LRANG2 NCORE(LP1)=LP1-1 1102 CONTINUE C WRITE(6,'(/'' DATA FOR USE IN IDENTIFYING BOUND LEVELS'')') WRITE(6,'(/'' ID TARGET STATES'')') WRITE(6,607)(I,IDSAT(I),IDLAT(I),IDPAT(I),I=1,IDNST) WRITE(6,'(/'' ID N+1 ELECTRON STATES'')') WRITE(6,607)(I,IDSAT1(I),IDLAT1(I),IDPAT1(I),I=1,IDNST1) WRITE(6,'(/'' HIGHEST PRINCIPAL QUANTUM NUMBER ASSIGNED'', 1 '' IN THE ID STATES FOR EACH L'')') WRITE(6,'(10I4)')(NCORE(L),L=1,LRANG2) C ENDIF C C ***** END OF EXTRA CODING ***** C ENDIF C IF(IRAD.GT.0)THEN H=ACNUM IF(RTWOIN.LE.RZERO)THEN KP2=4*((MZPTS-1)/4)+1 RTWO=(KP2-1)*H+RZERO ELSE RTWO=RTWOIN N=(RTWO-RZERO)/H N=4*((N-1)/4)+5 KP2=N IF(KP2.GT.MZPTS)THEN KP2=4*((MZPTS-1)/4)+1 RTWO=(KP2-1)*H+RZERO ENDIF H=KP2-1 H=(RTWO-RZERO)/H ENDIF CALL BODE if(iam.eq.0)then OPEN(8,FILE='B00',STATUS='UNKNOWN',FORM='UNFORMATTED') REWIND(8) WRITE(8)NZED,NELC WRITE(8)KP2 WRITE(8)RZERO,H WRITE(8)(WBODE(IR),IR=1,KP2) WRITE(8)IPERT endif ENDIF C C READ R-MATRIX FILE FOR SLPI CASE C ******************************** C KASE=1 1000 CALL READ2 C C TERMINATE IF EOF ON H.DAT C IF(MORE2.EQ.-777)GO TO 3000 C ISLP=10000*NSPN2+100*LRGL2+NPTY2 IF(IORDER.EQ.0)THEN C CHECK WHETHER THIS CASE REQUIRED DO 1001 K=KASE,KSLP IF(ISLP.NE.MSLP(K))GOTO 1001 KK=K GOTO 1002 1001 CONTINUE C CASE NOT FOUND GOTO 2000 C CASE FOUND ELSE KK=KASE IF(ISLP.EQ.MSLP(KASE))GO TO 1002 IF(MORE2.EQ.0)THEN WRITE(6,650) WRITE(6,651)MSLP(KASE) WRITE(6,652) KASE=KASE+1 IF(KASE.GT.KSLP)GO TO 3000 REWIND(10) DO I=1,5 READ(10) ENDDO ENDIF GO TO 1000 ENDIF C 1002 IS=NSPN2 IL=LRGL2 IP=NPTY2 IF(IOPT2.EQ.1)THEN XMIN=XMINST(KK) XMAX=XMAXST(KK) BDX=BDXST(KK) ELSE MXE=MXEST(KK) DO 1003 M=1,MXE 1003 EST(M)=ERDST(KK,M) ENDIF C RE-ORDER STORED UNIT-5 DATA IF(KK.NE.KASE)THEN MSLP(KK)=MSLP(KASE) IF(IOPT2.EQ.1)THEN XMINST(KK)=XMINST(KASE) XMAXST(KK)=XMAXST(KASE) BDXST(KK)=BDXST(KASE) ELSE MXEST(KK)=MXEST(KASE) DO 1004 M=1,MXEST(KK) 1004 ERDST(KK,M)=ERDST(KASE,M) ENDIF ENDIF C C PROCESS SLPI CASE C ***************** C CALL SCALE2 ECH1=ECH(1) C C SCAN FOR ENERGIES IF IOPT2.EQ.1 C ******************************* C IF(IOPT2.EQ.1)THEN MXE=0 IF(XMAX.LT.0)XMAX=1./SQRT(ECH1+1./(XMAX**2)) C RANGE OF EPSILON EPSMIN=-1./(XMIN*XMIN) EPSMAX=-1./(XMAX*XMAX) C RANGE OF ETOT EMIN=ECH1+EPSMIN EMAX=ECH1+EPSMAX C SCAN FOR FIRST POLE. ASSUME POLES IN DESCENDING ORDER. X0=XMIN E0=EMIN C FIND FIRST POLE ABOVE EMIN DO K=MNP2,1,-1 K0=K EP0=VALUE(K) IF(EP0.GT.E0)GOTO 20 ENDDO C CHECK NEAREST POLE BELOW EMIN 20 IF(K0.LT.MNP2)THEN KM=K0+1 EPM=VALUE(KM) IF(ABS(EPM-E0).LT.ABS(EP0-E0))THEN EP0=EPM K0=KM ENDIF ENDIF C NEXT LOWER POLE IF(K0.LT.MNP2)THEN XPM=1./SQRT(ECH1-VALUE(K0+1)) ELSE XPM=0. ENDIF C START LOOP FOR SCAN C CASE OF NEAREST POLE IN THE CONTINUUM 30 IF(EP0.GT.ECH1)THEN X1=XMAX GOTO 40 ENDIF XP0=1./SQRT(ECH1-EP0) C FIND RANGE OF X FOR POLE AT EP0 C NEXT POLE AT EP1 K1=K0-1 EP1=VALUE(K1) C CASE OF EP1 IN THE CONTINUUM IF(EP1.GT.ECH1)THEN X1=XMAX GOTO 40 ENDIF C EFFECTIVE QUANTUM NUMBER FOR NEAREST POLE IS XP0 C EFFECTIVE QUANTUM NUMBER FOR NEXT POLE XP1=1./SQRT(ECH1-EP1) C CALCULATE RANGE X1=0.5*(XP0+XP1) IF(X1.GT.XMAX)X1=XMAX C INTERVAL FOR SCAN 40 IF(IPRINT.EQ.1)THEN WRITE(6,609)K0,EP0 IF(ECH1.GT.EP0)THEN WRITE(6,611)XP0 ELSE WRITE(6,612) ENDIF ENDIF C CALL SCAN CALL SCAN(X0,X1,BDX,XPM) IF(MXE.GT.MZEST)THEN MXE=MZEST WRITE(6,608)MXE GO TO 50 ENDIF C CHECK WHETHER X1=XMAX IF(X1.NE.XMAX)THEN X0=X1 EP0=EP1 XPM=XP0 XP0=XP1 K0=K1 GOTO 30 ENDIF C SCAN COMPLETED FOR CASE OF IOPT2=1 WRITE(6,610)MXE ENDIF C C INITIALISE FOR CASE OF IOPT2.EQ.2 C ********************************* C IF(IOPT2.EQ.2)THEN I0=MNP2 DO M=1,MXE E=EST(M) XSTORE(M)=1./SQRT(ECH1-E) DO K=I0,1,-1 K1=K IF(VALUE(K).GT.E)GOTO 110 ENDDO 110 IF(K1.LT.MNP2)THEN IF(ABS(VALUE(K1+1)-E).LT.ABS(VALUE(K1)-E))K1=K1+1 ENDIF KSTORE(M)=K1 I0=K1 ENDDO ENDIF C C FINAL CALCULATION OF ENERGIES C ***************************** C 50 IF(IPERT.EQ.0)WRITE(6,641) C FOR IRAD.EQ.1 WRITE TO UNITS 8 AND 9 IF(IRAD.GT.0)THEN C WRITE(8)IS,IL,IP !CPB 2005 need to gather IS,IL,IP from all procs if(iam.eq.0)then do JJJ=0,nproc-1 write(8)IST(JJJ),ILT(JJJ),IPT(JJJ) C write(0,*)IST(JJJ),ILT(JJJ),IPT(JJJ) enddo endif deallocate(IST,ILT,IPT) NAME='B'//NUM((iam+1)/10)//NUM((iam+1)-10*((iam+1)/10)) OPEN(9,FILE=NAME,STATUS='UNKNOWN',FORM='UNFORMATTED') REWIND(9) WRITE(9)IS,IL,IP WRITE(9)MNP2,NCHF WRITE(9)(ECH(I),I=1,NCHF) WRITE(9)(CC(I),I=1,NCHF) IF(IPERT.EQ.1)WRITE(9)(((CF(I,J,L),I=1,NCHF),J=1,NCHF),L=1,2) WRITE(9)MXE ENDIF IF(IPRINT.EQ.1)WRITE(6,620) C C START ENERGY LOOP C DO 60 M=1,MXE C X=XSTORE(M) XR=X KR=KSTORE(M) CALL UTRANS(KR) ICONV=0 X0=X DO J=1,5 CALL BOUND(KR,X) XITER(J)=X IF(ABS(X-X0).LT.AC)THEN JR=J GOTO 56 ENDIF X0=X ENDDO C ITERATIONS FOR ENERGY ICONV=1 C WRITE FOR IPRINT.EQ.1 56 IF(IPRINT.EQ.1)THEN IF(IOPT2.EQ.1)THEN WRITE(6,630)M,XR ELSE WRITE(6,634)M,XR ENDIF ENDIF C WRITE DEPENDING ON CONVERGENCE IF(ICONV.EQ.0)THEN IF(IPRINT.EQ.1)WRITE(6,631)(XITER(J),J=1,JR) ELSE IF(IPRINT.EQ.1)THEN WRITE(6,632)(XITER(J),J=1,5) ELSE WRITE(6,635)M,(XITER(J),J=1,5) ENDIF ENDIF C OBTAIN VECTOR CALL VECT(KR,X) C C STORE EFFECTIVE QUANTUM NUMBERS AND WEIGHTS DO I=1,NCHF FNST(I,M)=FKNU(I) WT(I,M)=100.*(XVECT(I)**2) ENDDO C IF(IPRINT.EQ.1)WRITE(6,633)(XVECT(J),J=1,NCHF) C SAVE FOR WRITING TO 7 XSTORE(M)=XITER(JR) C C WRITE DATA FOR RADIATIVE CALCULATIONS TO UNIT 9 C *********************************************** C IF(IRAD.GT.0)THEN C C ETOT X=XSTORE(M) ETOT=ECH1-1./(X*X) WRITE(9)ETOT,(XVECT(N),N=1,NCHF) XNUMAX=MAX(XNUMAX,X) C C AVECT C BVECT=Q*XVECT CALL MTV(Q,XVECT,BVECT,MZCHF,NCHF) C C AVECT=((VALUE-ETOT)**-1)*WMAT(TRANSPOSE)*BVECT DO K=1,MNP2 S=0. DO J=1,NCHF S=S+WMAT(J,K)*BVECT(J) ENDDO AVECT(K)=S/(VALUE(K)-ETOT) ENDDO C C NEXT TWO STATEMENTS ADDED TO HANDLE BUTTLE CORRECTIONS DO J=1,NCHF AVECT(MNP2+J)=BVECT(J) ENDDO C C IN NEXT STATEMENT, '+NCHF' ADDED TO HANDLE BUTTLE CORRECTIONS WRITE(9)(AVECT(K),K=1,MNP2+NCHF) C C NORMALISED ZERO-ORDER FUNCTIONS AT RZERO AND THEIR DERIVATIVES AR=1./RZERO DO N=1,NCHF PB(N)=FS(N,1)*XVECT(N) PBP(N)=FSP(N)*XVECT(N) PBPP(N)=(AR*(AR*CC(N)-2.)-EPS(N))*PB(N) ENDDO WRITE(9)(PB(N),N=1,NCHF) WRITE(9)(PBP(N),N=1,NCHF) WRITE(9)(PBPP(N),N=1,NCHF) C C FIRST-ORDER CORRECTIONS TO FUNCTIONS AT RZERO IF(IPERT.EQ.1)THEN AR2=AR*AR AR3=AR2*AR DO I=1,NCHF PBPP(I)=0. ENDDO DO J=1,NCHF DO I=1,NCHF IF(LAMP(I,J).EQ.2)THEN PBPP(I)=PBPP(I)-BW(I,J)*AR2*PB(J) ELSEIF(LAMP(I,J).EQ.3)THEN PBPP(I)=PBPP(I)-BW(I,J)*AR3*PB(J) ENDIF ENDDO ENDDO DO I=1,NCHF PB(I)=0. PBP(I)=0. ENDDO DO J=1,NCHF DO I=1,NCHF PB(I)=PB(I)+(FC(I,1)*ASS(I,J)-FS(I,1)*ACS(I,J))*XVECT(J) PBP(I)=PBP(I)+(FCP(I)*ASS(I,J)-FSP(I)*ACS(I,J))*XVECT(J) ENDDO ENDDO DO I=1,NCHF PBPP(I)=PBPP(I)+(AR*(AR*CC(I)-2.)-EPS(I))*PB(I) ENDDO ENDIF C C NORMALISED OUTER-REGION ZERO-ORDER FUNCTIONS DO KP=1,KP2 DO II=1,NCHF FS(II,KP)=FS(II,KP)*XVECT(II) ENDDO ENDDO ACSQ=AC*AC R=RTWO H2=H*2. DO KP=KP2-2,1,-2 R=R-H2 SUM=0. DO II=1,NCHF SUM=SUM+FS(II,KP)**2 ENDDO IF(SUM*R.GT.ACSQ)THEN KPM=KP+2 GOTO 67 ENDIF ENDDO KPM=1 67 IF(IPRINT.GT.0)WRITE(6,671)M,KPM C WRITE(9)KPM IF(KPM.GT.1)THEN WRITE(9)(FS(II,KPM-1),FS(II,KPM),II=1,NCHF) ELSE WRITE(9)(FS(II,1),II=1,NCHF) ENDIF C C WRITE FIRST-ORDER CORRECTIONS IF(IPERT.EQ.1)THEN WRITE(9)(PB(N),N=1,NCHF) WRITE(9)(PBP(N),N=1,NCHF) WRITE(9)(PBPP(N),N=1,NCHF) ENDIF C C END OF IRAD LOOP ENDIF C C C END ENERGY LOOP C 60 CONTINUE C C C IDENTIFY LEVELS C *************** C IF (IRAD .LT. 2.AND.IDNST1.GE.0) THEN CALL IDENTFY(ID,IDDD) ITARG1(KASE)=ITARG(1) END IF C C C FINAL WRITES TO UNITS 6, 7, 8 AND 9 C ************************************ C IF(IRAD.LT.2) THEN WRITE(7,704)IS,IL,IP,ITARG(1) WRITE(7,705)MXE,ECH1 IF(IDNST1.GE.0)THEN WRITE(11,704)IS,IL,IP WRITE(11,705)MXE ENDIF ENDIF WRITE(6,636) DO 70 M=1,MXE X=XSTORE(M) E=ECH1-1./(X*X) WRITE(6,637)M,X,E IF(IRAD.LT.2.AND.IDNST1.GE.0)WRITE(11,1111)M,ID(M),E,IDDD(M) IF(IRAD.LT.2)WRITE(7,706)E,X 70 CONTINUE C C WRITE INFORMATION FOR IDENTIFICATIONS IF(IDNST1.LT.0)THEN WRITE(6,660)(M,M=1,NCHF) WRITE(6,662) DO 71 M=1,MXE WRITE(6,661)M,(WT(J,M),J=1,NCHF) 71 CONTINUE WRITE(6,663)(M,M=1,NCHF) WRITE(6,662) DO 72 M=1,MXE WRITE(6,664)M,(FNST(J,M),J=1,NCHF) 72 CONTINUE ENDIF C C FINAL ORGANISATION IF(IRAD.GT.0)WRITE(6,653) IF(IRAD.GT.0)CLOSE(9) IF(KASE.EQ.KSLP)GOTO 3000 C WRITE FOR CASE OF MORE2.EQ.0 IF(MORE2.EQ.0.AND.IORDER.EQ.0)THEN IF(IRAD.LT.2)WRITE(7,707) if(iam.eq.0)CLOSE(8) CLOSE(7,STATUS='KEEP') WRITE(6,680) IF(KASE.NE.KSLP)THEN WRITE(6,650) DO I=KASE+1,KSLP WRITE(6,651)MSLP(I) ENDDO WRITE(6,652) ENDIF GO TO 1997 ENDIF C SLPI CASE COMPLETED KASE=KASE+1 IF(KASE.GT.KSLP)GO TO 3000 C IF(IORDER.NE.0)THEN REWIND(10) DO I=1,5 READ(10) ENDDO GO TO 1000 ENDIF C C CHECK WHETHER MORE DATA ON R-MATRIX FILE C **************************************** C 2000 IF(MORE2.EQ.0)THEN WRITE(6,650) DO I=KASE,KSLP WRITE(6,651)MSLP(I) ENDDO WRITE(6,652) GOTO 3000 ENDIF C C READ NEXT CASE ON R-MATRIX FILE GOTO 1000 C C TERMINATION C *********** 3000 IF(IRAD.LT.2)THEN WRITE(7,707) CLOSE(7,STATUS='KEEP') IF(IDNST1.GE.0)WRITE(11,707) WRITE(6,680) ENDIF IF((IRAD.GT.0).and.(iam.eq.0))THEN WRITE(8)-1,-1,-1 WRITE(8)NINT(XNUMAX) CLOSE(8) ENDIF 1997 IF(IRAD.LT.2.AND.IDNST1.GE.0) THEN CALL IDFORM(ITARG1,KSLP,ENAT,NAST,ISAT,LAT,IDNST1) CLOSE(11,STATUS='KEEP') ENDIF C C SUN C C call comm_barrier() C call cpu_time(timef) time=timef-timei write(6,999) time/60.,nproc 999 format(//1x,'CPU TIME=',f9.3,' MIN -- processors=:',i4) C C call comm_finalize() C stop C C C C STOP C C FORMATS C ******* 600 FORMAT(5X,'DATA READ FROM UNIT 5'/ 1 10X,'IPRINT = ',I2/ 2 10X,'IRAD = ',I2/ 3 10X,'IPERT = ',I2/ 4 10X,'AC = ',E9.2/ 5 10X,'RONE = ',F6.4/ 6 10X,'RTWOIN = ',F7.2) 601 FORMAT(10X,'IOPT2 = 1 FOR SCAN') 602 FORMAT(10X,'IOPT2 = 2 FOR INPUT ENERGY ESTIMATES') 603 FORMAT(10X,'VALUES OF 10000*IS+100*IL+IP',//(35X,I10)) 604 FORMAT(/) 605 FORMAT(10X,30('*')/10X,'NASTD = ',I5, + ' IS LARGER THAN MZTAR = ',I4/10X,30('*')//) 606 FORMAT(10X,'IOPT2 INITIALLY NEGATIVE',/ + 10X,'NASTD =',I3/10X,'NLEV =',20I3) 607 FORMAT(1X,I2,3X,3I3) 608 FORMAT(//' ***WARNING: LOWEST',I5,' ROOTS DETERMINED, ' X,'SCAN INCOMPLETE, RECOMPILE WITH LARGER MZEST') 609 FORMAT(//' NEAREST POLE, K0 =',I4,' EP0 = ',F10.6) 610 FORMAT(//10X,'SCAN FOR ROOTS COMPLETED, ', 1 I3,' ROOTS FOUND'/) 611 FORMAT(' XP0 = ',F10.6/) 612 FORMAT(' POLE ABOVE LIMIT') 620 FORMAT(//10X,'EFFECTIVE QUANTUM NUMBERS AND VECTORS'/ 1 10X,37('*')/) 630 FORMAT(/1X,'LEVEL',I3/10X, 1 'NEGLECTING MULTIPOLE POTENTIAL, EFFECTIVE QUANTUM NUMBER IS'/ 2 10X,F10.5/10X, 3 'INCLUDING MULTIPOLE POTENTIAL, ITERATIONS FOR EFFECTIVE' 4 ,' QUANTUM NUMBER ARE') 631 FORMAT(10X,5F10.5) 632 FORMAT(10X,5F10.5,' ** NOT CONVERGED **') 633 FORMAT(10X,'X VECTOR'/(10X,10F10.5)) 634 FORMAT(/1X,'LEVEL',I3/10X, 1 'INPUT ESTIMATE FOR EFFECTIVE QUANTUM NUMBER IS'/ 2 10X,F10.5/10X, 3 'INCLUDING MULTIPOLE POTENTIAL, ITERATIONS FOR EFFECTIVE' 4 ,' QUANTUM NUMBER ARE') 635 FORMAT(//10X,30('*')/ 1 10X,'ITERATIONS NOT CONVERGED FOR LEVEL ',I3/ 2 10X,'XITER = ',5F10.5/ 3 10X,30('*')//) 636 FORMAT(//10X,'EFFECTIVE QUANTUM NUMBERS FNU AND ENERGIES E'/ 1 10X,44('-')// 2 10X,'(FNU IS RELATIVE TO LOWEST STATE IN EXPANSION.'/ 3 11X,'E IS Z-SCALED ENERGY RELATIVE TO GROUND STATE.)'// 4 17X,'M',4X,'FNU',7X,'E'//) 637 FORMAT(15X,I3,F10.5,E16.6) 640 FORMAT(//3X,10('*'),3X,'TERMINATION AT END OF DATA ON' 1 ,' R-MATRIX FILE',3X,10('*')//) 641 FORMAT(//10X,40('*')/ 1 10X,10('*'),' RUN WITH IPERT = 0 ',10('*')/ 2 10X,40('*')//) 650 FORMAT(//10X,30('*')/ 1 10X,'NO DATA ON R-MATRIX FILE FOR'/ 2 10X,'10000*IS+100*IL+IP = '/) 651 FORMAT(10X,I10) 652 FORMAT(/10X,30('*')//) 653 FORMAT(//10X,'IRAD DATA WRITTEN TO FILES'/) 660 FORMAT(//6X,'CHANNEL PERCENTAGE WEIGHTS FOR CHANNEL NUMBER = '/ 1 (6X,8I8)) 661 FORMAT(I3,3X,8(2PE8.1)/6X,8E8.1) 662 FORMAT(/) 663 FORMAT(//6X,'EFFECTIVE QUANTUM NUMBERS, CHANNEL NUMBER = '/ 1 (6X,12I7)) 664 FORMAT(I3,3X,12F7.3/(6X,12F7.3)) 671 FORMAT(5X,'M = ',I4,', NUMBER OF POINTS FOR FUNCTION', + ' OUTPUT, KPM = ',I5) 680 FORMAT(10X,38('*')/10X,'* ENERGY LEVELS WRITTEN TO FILE ELEV *'/ + 10X,38('*')//) 6000 FORMAT(1X,70('+')//10X,'PROGRAM STGB'/10X,12('*')// 1 3X,'HISTORY OF MODIFICATIONS'// 2 5X,'- CONTROL OF PRINT LEVEL'/ 3 5X,'- FILES FOR CALCULATION OF RADIATIVE DATA'/ 4 5X,'- CONTROL OF INCLUSION OF PERTURBATION POTENTIALS'/ 5 5X,'- PROGRAM OPENS FILES B00, B01, ....'/ 6 5X,'- THE DATA ON THESE FILES ARE FOR NUMEROV INTEGRATIONS' 7 ,' IN STGBB,BF'/ 8 5X,'- ONLY ONE WARNING FOR REQUIRED RTWO LARGER THAN', 9 ' MAXIMUM ALLOWED'/ 1 5X,'- NEARLY DEGENERATE TARGET LEVELS CAN BE COMBINED'/ 2 5X,'- NAMELIST INPUT'/ 3 5X,'- DIMENSIONS FROM INCLUDE, NO PRE-PROCESSING'/ 4 5X,'- NUMERICAL PROBLEMS ASSOCIATED WITH SIMULTANEOUS WEAK'/ 5 5X,' AND STRONGLY CLOSED CHANNELS SOLVED' 6 //1X,70('+')//) 6002 FORMAT(//10X,'COMPILED FOR DIMENSIONS -'// + 15X,'CHANNELS MZCHF =',I6/ + 15X,'TARGET STATES MZTAR =',I6/ + 15X,'MULTIPOLES MZLMX =',I6/ + 15X,'SMALL L VALUES MZLP1 =',I6/ + 15X,'OUTER-REGION RADIAL POINTS MZPTS =',I6/ + 15X,'NUMBER BOUND-STATE ENERGIES MZEST =',I6/ + 15X,'R-MATRIX POLES MZMNP =',I6/ + 15X,'S, L, PI CASES MZSLP =',I6/ + 15X,'COEFFICIENTS FOR THETA MZTET =',I6/ + 15X,'TERMS IN BUTTLE FIT MZNRG =',I6//) 700 FORMAT(3I3,I9,I4) 701 FORMAT(1PE11.3) 702 FORMAT(F10.6) 703 FORMAT(' 2') 704 FORMAT(I3,2I4,5X,I2) 705 FORMAT(I3,8X,F14.6) 706 FORMAT(1PE16.8,0PF12.7) C 707 FORMAT(2X,'0',2(3X,'0')) 707 FORMAT(2X,'0',2X,'-1',3X,'0') 1110 FORMAT(1X,A4,I4) 1111 FORMAT(1X,I3,2X,A9,1PE14.6,I5) C C END C C********************************************************************* C SUBROUTINE ACSUB IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC C ACNUM=(24.*AC)**.166666667 ACJWBK=(6.*AC)**.2 ACZP=16.*AC LACC=0 IF(AC.LT.1.E-3)LACC=2 IF(AC.LT.1.E-4)LACC=4 C RETURN END C C********************************************************************* C SUBROUTINE ALPHAB C C CALCULATES ALPHA INTEGRALS C IMPLICIT REAL*8 (A-H,O-Y) C INCLUDE 'PARAM' C COMMON/CALP/ASS(MZCHF,MZCHF),ASC(MZCHF,MZCHF),ACS(MZCHF,MZCHF) 1 ,ACC(MZCHF,MZCHF) COMMON/CNTRL/IPRINT,IRAD,IPERT COMMON/COULSC/FS(MZCHF,MZPTS),FSP(MZCHF),FC(MZCHF,MZPTS) 1 ,FCP(MZCHF) COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CPOT/LAMP(MZCHF,MZCHF),BW(MZCHF,MZCHF) COMMON/CBODE/WBODE(MZPTS) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC C C C INITIALISE ALPHA TO ZERO DO I=1,NCHF DO J=I,NCHF ASS(J,I)=0. ASC(J,I)=0. ACS(J,I)=0. ACC(J,I)=0. ENDDO ENDDO C C C CONTRIBUTION FROM RZERO TO RTWO R=RZERO-H C START LOOP ON RADIAL POINTS DO K=1,KP2 W=WBODE(K) R=R+H X=1./R C CLOSED-CLOSED PART DO I=NCHOP1,NCHF DO J=I,NCHF LIJ=LAMP(J,I) IF(LIJ.NE.1)THEN WIJ=W*BW(J,I)*X**LIJ T1=WIJ*FS(J,K) ASS(J,I)=ASS(J,I)+FS(I,K)*T1 ASC(J,I)=ASC(J,I)+FC(I,K)*T1 ACS(J,I)=ACS(J,I)+FS(I,K)*WIJ*FC(J,K) ENDIF ENDDO ENDDO ENDDO C C C C ASYMPTOTIC INTEGRALS FOR R=RTWO TO INFINITY IF(IPRINT.EQ.2)WRITE(6,750) 750 FORMAT(/' J,I AND ASS(J,I), ASC(J,I), ACS(J,I), ACC(I,J)'/) C C CLOSED-CLOSED PART 180 IF(NCHOP.EQ.NCHF)GOTO 300 DO I=NCHOP1,NCHF DO J=I,NCHF LIJ=LAMP(J,I) IF(LIJ.NE.1)THEN NLAG=2*LIJ+LACC BIJ=BW(J,I) CALL TLAG(I,J,NLAG,LIJ,T1,T2,T3) ASS(J,I)=ASS(J,I)+T1*BIJ ASC(J,I)=ASC(J,I)+T2*BIJ ACS(J,I)=ACS(J,I)+T3*BIJ IF(IPRINT.EQ.2) WRITE(6,760)J,I,ASS(J,I),ASC(J,I),ACS(J,I) 760 FORMAT(2I5,4E14.6) ENDIF ENDDO ENDDO C C C SYMMETRISE ALPHA 300 IF(NCHF.EQ.1)RETURN DO I=2,NCHF K=I-1 DO J=1,K ASS(J,I)=ASS(I,J) ACC(J,I)=ACC(I,J) ACS(J,I)=ASC(I,J) ASC(J,I)=ACS(I,J) ENDDO ENDDO C C RETURN END C C********************************************************************* C SUBROUTINE BLOCK C C PROVIDES C BLOCK DATA C CALLED AS SUBROUTINE TO AVOID LINKAGE PROBLEMS WITH LIBRARIES C C DATA FOR QUADRATURES - C LAGUERRE AND LEGENDRE QUADRATURES WITH NUMBERS OF POINTS C N = 2, 4, 6, 8 AND 10 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) C DATA XLAG/ 1 .58578644,3.4142136, 2 .32254769,1.7457611,4.5366203,9.3950709, 3 .22284660,1.1889321,2.9927363,5.7751436,9.8374674, 4 15.982874, 5 .17027963,.90370178,2.2510866,4.26670017,7.0459054, 6 10.758516,15.7406786,22.8631317, 7 .13779347,.72945455,1.8083429,3.4014337,5.5524961, 8 8.3301527,11.8437858,16.279258,21.996586,29.920697/ DATA WLAG/ 1 .85355339,.14644661, 2 .60315410,.35741869,.38887909E-1,.53929471E-3, 3 .45896467,.41700083,.11337338,.10399197E-1, 4 .26101720E-3,.89854791E-6, 5 .36918859,.41878678,.17579499,3.3343492E-2,2.7945362E-3, 6 9.0765088E-5,8.4857467E-7,1.0480012E-9, 7 .30844112,.40111993,.21806829,6.2087456E-2,9.5015170E-3, 8 7.5300839E-4,2.8259233E-5,4.2493140E-7,1.8395648E-9, 9 9.9118272E-13/ DATA XLEG,WLEG/.577350269, 1 .339981044,.861136312, 2 .238619186,.661209386,.932469514, 3 .183434642,.525532410,.796666477,.960289856, 4 .148874339,.433395394,.679409568,.865063367,.973906529, 5 1., 6 .652145159,.347854845, 7 .467913935,.360761573,.171324492, 8 .362683783,.313706646,.222381034,.101228536, 9 .295524225,.269266719,.219086363,.149451349,.066671344/ C RETURN END C C********************************************************************* C SUBROUTINE BODE C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CBODE/WBODE(MZPTS) C W=.31111111*H WBODE(1)=W WBODE(KP2)=W W=1.4222222*H M=KP2-1 DO 10 K=2,M,2 10 WBODE(K)=W W=.53333333*H M=KP2-2 DO 20 K=3,M,4 20 WBODE(K)=W IF(KP2.EQ.5)RETURN W=.62222222*H M=KP2-4 DO 30 K=5,M,4 30 WBODE(K)=W C RETURN END C C********************************************************************* C SUBROUTINE BOUND(KP,XR) C C ITERATION FOR ENERGY C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (ZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/COULSC/FS(MZCHF,MZPTS),FSP(MZCHF),FC(MZCHF,MZPTS) 1 ,FCP(MZCHF) COMMON/CALP/ASS(MZCHF,MZCHF),ASC(MZCHF,MZCHF),ACS(MZCHF,MZCHF) 1 ,ACC(MZCHF,MZCHF) COMMON/CRMAT/RMAT(MZCHF,MZCHF),RMATD(MZCHF,MZCHF) 1,UMAT(MZCHF,MZCHF),GAMMA COMMON/CVECT/A(MZCHF,MZCHF),B(MZCHF,MZCHF),C(MZCHF,MZCHF), 1 P(MZCHF,MZCHF),Q(MZCHF,MZCHF),XVECT(MZCHF) COMMON/CNTRL/IPRINT,IRAD,IPERT C DIMENSION SP(MZCHF),SQ(MZCHF),SPD(MZCHF),SQD(MZCHF) C C ETOT=ECH(1)-1./(XR*XR) C DE=VALUE(KP)-ETOT DELE=2.*AC/(XR**3) C CALL POINTB CALL COULB CALL RMATB2(KP) C DO I=1,NCHF SP(I)=FS(I,1) SQ(I)=FSP(I)-BSTO*FS(I,1) SPD(I)=FC(I,1) SQD(I)=FCP(I)-BSTO*FC(I,1) ENDDO C DO J=1,NCHF DO I=1,NCHF P(I,J)=0. Q(I,J)=0. ENDDO ENDDO DO I=1,NCHF P(I,I)=SP(I) Q(I,I)=SQ(I) ENDDO IF(IPERT.EQ.1)THEN CALL ALPHAB DO J=1,NCHF DO I=1,NCHF P(I,J)=P(I,J)+SPD(I)*ASS(I,J)-SP(I)*ACS(I,J) Q(I,J)=Q(I,J)+SQD(I)*ASS(I,J)-SQ(I)*ACS(I,J) ENDDO ENDDO ENDIF C C COMPUTE C=P-T*Q CSTRTNBL CALL MATMUL(RMAT,Q,C,MZCHF,NCHF) CENDNBL CSTRTBL C CALL DGEMM('N','N',NCHF,NCHF,NCHF,ONE,RMAT,MZCHF,Q,MZCHF, C X ZERO,C,MZCHF) CENDBL C DO J=1,NCHF DO I=1,NCHF C(I,J)=-C(I,J) ENDDO C(J,J)=SP(J)+C(J,J) ENDDO C C COMPUTE A=U(TRANSPOSE)*C CSTRTNBL CALL MTMLT(UMAT,C,A,MZCHF,NCHF) CENDNBL CSTRTBL C CALL DGEMM('T','N',NCHF,NCHF,NCHF,ONE,UMAT,MZCHF,C,MZCHF, C X ZERO,A,MZCHF) CENDBL C C COMPUTE C=TD*Q+T*QD-PD DO J=1,NCHF DO I=1,NCHF C(I,J)=RMATD(I,J)*SQ(J)+RMAT(I,J)*SQD(J) ENDDO ENDDO DO I=1,NCHF C(I,I)=C(I,I)-SPD(I) ENDDO C C COMPUTE B=U(TRANSPOSE)*C CSTRTNBL CALL MTMLT(UMAT,C,B,MZCHF,NCHF) CENDNBL CSTRTBL C CALL DGEMM('T','N',NCHF,NCHF,NCHF,ONE,UMAT,MZCHF,C,MZCHF, C X ZERO,B,MZCHF) CENDBL C C CHANGE LAST ROW OF B DO J=1,NCHF B(NCHF,J)=DE*B(NCHF,J)+A(NCHF,J)+GAMMA*UMAT(J,NCHF)*SQD(J) ENDDO C C CHANGE LAST ROW OF A DO J=1,NCHF A1=0. DO K=1,NCHF A1=A1+UMAT(K,NCHF)*Q(K,J) ENDDO A(NCHF,J)=DE*A(NCHF,J)-GAMMA*A1 ENDDO C C CALL SOLV(A,B,NCHF,DELE,FB,XVECT,ICONV) C C C WRITE MESSAGE IF NOT CONVERGED IN SOLV IF(ICONV.EQ.1)WRITE(6,600)ETOT,XR,FB C ETOT=ETOT+FB C XR=1./SQRT(ECH(1)-ETOT) C RETURN C 600 FORMAT(//5X,10('*'), 'WARNING, ITERATIONS IN SOLVE NOT CONVERGED'/ + 17X,'INPUT TO SOLV : ETOT = ',E14.6,' XR = ',F9.5/ + 17X,'OUTPUT FROM SOLV : ETOT INCREMENT = ',E12.4/ + 17X,'SHOULD BE O.K. IF FINAL ITERATIONS FOR ENERGY CONVERGE'/ + 5X,10('*')//) END C C********************************************************************* C SUBROUTINE BSCT(XA,XB,DA,DB,IW) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/CITER/ECH1,XSTORE(MZEST),K0,KSTORE(MZEST) C C INITIALISATIONS X1=XA X2=XB D1=DA D2=DB IFL1=1 IFL2=1 IF(D1.LT.0.)IFL1=-1 IF(D2.LT.0.)IFL2=-1 C C START OUTER LOOP FOR INTERVAL HALVING DO I=1,5 XH=.5*(X1+X2) DH=DETERM(XH) IFLH=1 IF(DH.LT.0.)IFLH=-1 IF(IFL1*IFLH.LT.0.)THEN X2=XH D2=DH IFL2=IFLH ELSE X1=XH D1=DH IFL1=IFLH ENDIF C C START INNER LOOP FOR BISECTIONS T=MAX(ABS(D1),ABS(D2)) DO J=1,5 XR=(X1*(D2/T)-X2*(D1/T))/((D2-D1)/T) IF(ABS(XR-X2).LT.AC)THEN IF(IW.EQ.1)THEN DO M=1,MXE IF(ABS(XR-XSTORE(M)).LT.AC)go to 1000 ENDDO ENDIF MXE=MXE+1 IF(MXE.GT.MZEST)go to 1000 XSTORE(MXE)=XR KSTORE(MXE)=K0 go to 1000 ENDIF DR=DETERM(XR) IFLR=1 IF(DR.LT.0.)IFLR=-1 IF(IFL2*IFLR.LT.0.)THEN X1=X2 D1=D2 IFL1=IFL2 ENDIF X2=XR D2=DR IFL2=IFLR ENDDO ENDDO C C LOOPS COMPLETED, NO ROOTS FOUND C WRITE(6,600)XA,XB IF(IW.EQ.0)THEN WRITE(6,610) ELSE WRITE(6,620) ENDIF C 1000 RETURN C 600 FORMAT(' ***** BISECTIONS NOT CONVERGED FOR ROOT BETWEEN' 1,' XA =',F10.4,' AND XB =',F10.4,' *****') 610 FORMAT(' ***** BSCT CALLED FROM SCAN') 620 FORMAT(' ***** BSCT CALLED FROM QUAD') C END C C********************************************************************* C REAL*8 FUNCTION BUT0(NBUT,U) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C LOGICAL POLE COMMON/CBUT/FKN(0:MZNRG),UKN(0:MZNRG) C BUT0=0. C C CASE OF U.GT.0.04 IF(U.GT..04)THEN FK=SQRT(U) POLE=.FALSE. DO 10 N=0,NBUT IF(ABS(FK-FKN(N)).GT..3)THEN BUT0=BUT0+1./(U-UKN(N)) ELSE POLE=.TRUE. D1=FK-FKN(N) ENDIF 10 CONTINUE IF(POLE)THEN D2=D1**2 D=.33333333*D1*(1.+.066666667*D2*(1.+.0952381*D2)) BUT0=2.*BUT0+(D+1./(2.*FK-D1))/FK ELSE BUT0=2.*BUT0+TAN(FK)/FK ENDIF C C SUM FOR U.LT..04 ELSE DO 20 N=0,NBUT 20 BUT0=BUT0+1./(U-UKN(N)) C C CASE OF U.LT..04 AND U.GT.-.04 IF(U.GT.-.04)THEN BUT0=2.*BUT0+1.+.33333333*U*(1.+.4*U) C C CASE OF U.LT.-.04 ELSE FK=SQRT(-U) BUT0=2.*BUT0+TANH(FK)/FK ENDIF C ENDIF C RETURN END C C********************************************************************* C SUBROUTINE CBUT0(NBUT,U,B,C) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CBUT/FKN(0:MZNRG),UKN(0:MZNRG) C LOGICAL POLE C C B=0. C=0. C C CASE OF U.GT.0.04 IF(U.GT..04)THEN FK=SQRT(U) POLE=.FALSE. DO 10 N=0,NBUT IF(ABS(FK-FKN(N)).GT..3)THEN A=1./(U-UKN(N)) B=B+A C=C-A**2 ELSE POLE=.TRUE. D1=FK-FKN(N) ENDIF 10 CONTINUE IF(POLE)THEN D2=D1**2 D=.33333333*D1*(1.+.066666667*D2*(1.+.0952381*D2)) A=1./(2.*FK-D1) BB=(D+A)/FK D=.33333333*(1.+D2*(.2+.031746032*D2)) C=2.*C+.5*(D-A**2-BB)/U B=2.*B+BB ELSE T=TAN(FK) TK=T/FK B=2.*B+TK C=2.*C+.5*(1.+T**2-TK)/U ENDIF C C SUM FOR U.LT..04 ELSE DO 20 N=0,NBUT A=1./(U-UKN(N)) B=B+A 20 C=C-A**2 C C CASE OF U.LT..04 AND U.GT.-.04 IF(U.GT.-.04)THEN B=2.*B+1.+.33333333*U*(1.+.4*U) C=2.*C+.33333333*(1.+U*(.8+.48571429*U)) C C CASE OF U.LT.-.04 ELSE FK=SQRT(-U) T=TANH(FK) TK=T/FK B=2.*B+TK C=2.*C+.5*(1.-T**2-TK)/U ENDIF C ENDIF C RETURN END C C********************************************************************* C SUBROUTINE COULB C C CALCULATION OF COULOMB FUNCTIONS C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/COULSC/FS(MZCHF,MZPTS),FSP(MZCHF),FC(MZCHF,MZPTS) 1 ,FCP(MZCHF) COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CNTRL/IPRINT,IRAD,IPERT C C TOL0=1.0D-150 IFLAG=0 C DO 210 I=NCHOP1,NCHF C KP2T=KP2 RTWOT=RTWO C 205 CALL THETA(RTWOT,I,T,TP,TD,TDP) C FS(I,KP2T)=T FC(I,KP2T)=TD C C CHECK THAT WHITTAKER FUNCTION IS NOT VANISHINGLY SMALL - NRB 24/05/96 C IF(ABS(T).LT.TOL0.OR.ABS(TD).LT.TOL0)THEN IFLAG=1 KP2T=KP2T-4 IF(KP2T.LT.1)KP2T=1 RTWOT=(KP2T-1)*H+RZERO C IF(KP2T.LT.KP0)STOP 'SR.COUL: KP2T .LT. KP0' DO 206 K=KP2T,KP2T+4 FS(I,K)=0.0 FC(I,K)=0.0 206 CONTINUE IF(KP2T.EQ.1)GO TO 210 GO TO 205 ENDIF C IF(IPRINT.GT.1.AND.KP2T.LT.KP2)WRITE(6,725)I,RTWOT,RTWO C FSP(I)=TP FCP(I)=TDP C CALL NUMT(EPS(I),CC(I),RTWOT,H,KP2T,KP0,I) C 210 CONTINUE C IF(IFLAG.GT.0.AND.IPRINT.GT.1)WRITE(6,720) C IF(IPRINT.GT.1)THEN WRITE(6,701)EPS(1) WRITE(6,705)RTWO,KP2,H WRITE(6,700) DO 230 I=1,NCHF 230 WRITE(6,710)I,FS(I,1),FSP(I),FC(I,1),FCP(I) ENDIF C 700 FORMAT(//10X,'COULOMB FUNCTIONS S,SP,C AND CP'/) 701 FORMAT(//10X,' E = ',F10.6/11X,14('-')) 705 FORMAT(//' RTWO = ',F10.6,', KP2 = ',I4,', H = ',F10.6) 710 FORMAT(I5,4E15.6) 720 FORMAT(3X,'I',3X,'RTWOT',8X,'RTWO') 725 FORMAT(I5,2F10.4) C RETURN END C C********************************************************************* C REAL*8 FUNCTION DET(A,MZ,N) C C CALCULATES DETERMINANT WITH FULL PIVOTING C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (ZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (DETX=1.d200) C DIMENSION A(MZ,*) C DET=ONE C DO M=1,N-1 C P=ZERO DO J=M,N DO I=M,N Q=ABS(A(I,J)) IF(Q.GT.P)THEN P=Q IP=I JP=J ENDIF ENDDO ENDDO C IF(IP.NE.M)THEN DET=-DET DO J=M,N Q=A(IP,J) A(IP,J)=A(M,J) A(M,J)=Q ENDDO ENDIF IF(JP.NE.M)THEN DET=-DET DO I=M,N Q=A(I,JP) A(I,JP)=A(I,M) A(I,M)=Q ENDDO ENDIF C P=A(M,M) DET=DET*P if(abs(det).gt.detx)det=det/detx P=ONE/P DO J=M+1,N Q=A(M,J)*P DO I=M+1,N A(I,J)=A(I,J)-A(I,M)*Q ENDDO ENDDO C ENDDO C DET=DET*A(N,N) C RETURN END C C********************************************************************* C REAL*8 FUNCTION DETERM(X) C C CALCULATES DETERMINANT OF C A = DE * (P - R * Q) C WHERE C DE = VALUE(KP) - ETOT C USES ZERO-ORDER FUNCTIONS P AND Q C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (ZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/COULSC/FS(MZCHF,MZPTS),FSP(MZCHF),FC(MZCHF,MZPTS) 1 ,FCP(MZCHF) COMMON/CRMAT/RMAT(MZCHF,MZCHF),RMATD(MZCHF,MZCHF) 1,UMAT(MZCHF,MZCHF),GAMMA COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CITER/ECH1,XSTORE(MZEST),K0,KSTORE(MZEST) DIMENSION A(MZCHF,MZCHF),B(MZCHF,MZCHF),Q(MZCHF) C ETOT=ECH1-1./(X*X) DE=VALUE(K0)-ETOT C CALL POINTB CALL COULB CALL RMATB1(K0) C C FUNCTIONS Q DO J=1,NCHF Q(J)=FSP(J)-BSTO*FS(J,1) ENDDO C C A=P-T*Q DO J=1,NCHF QJ=Q(J) DO I=1,NCHF A(I,J)=-RMAT(I,J)*QJ ENDDO ENDDO DO I=1,NCHF A(I,I)=A(I,I)+FS(I,1) ENDDO C C C B=UT*A C CSTRTNBL CALL MTMLT(UMAT,A,B,MZCHF,NCHF) CENDNBL CSTRTBL C CALL DGEMM('T','N',NCHF,NCHF,NCHF,ONE,UMAT,MZCHF,A,MZCHF, C X ZERO,B,MZCHF) CENDBL C C CHANGE LAST ROW OF B DO J=1,NCHF B(NCHF,J)=B(NCHF,J)*DE-GAMMA*UMAT(J,NCHF)*Q(J) ENDDO C DETERM=DET(B,MZCHF,NCHF) C RETURN END C C*********************************************************************** C SUBROUTINE IDENTFY(ID,IDDD) C C TO IDENTIFY THE LEVELS AND WRITE FILE 11 IN THE OPACITY C ID FORMAT C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (N1TAR=10*MZTAR) CHARACTER*1 IDNUM1,IPNT CHARACTER*2 IDNUM,IDSEQ,IDSEQ1,IDTARG CHARACTER*9 ID,EXCLAM COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CDEGEN/NASTD,NASTR,NLEV(MZTAR),NCNATR(MZTAR),ENATR(MZTAR) COMMON/CID/WT(MZCHF,MZEST),FNST(MZCHF,MZEST),ISLP,IDNST,IDNST1, 1 IDSLP(MZTAR),IDSLP1(N1TAR),NCORE(MZLP1) COMMON/CIDFOR/IPAT(MZTAR) DIMENSION IDNUM(0:99),IDNUM1(0:9),IDTARG(MZCHF),IDTEMP(MZEST), 1 ITASLP(MZTAR),WTEMP(MZCHF,MZEST),FN(MZEST), 2 MTEMP(MZEST),MDUP(MZEST),NDUP(MZEST),JID(MZEST), 3 IDD(MZEST),NJMAX(MZCHF),NDJM(MZCHF,MZEST),DEFECT(MZCHF) C DIMENSION IDOUBT(MZEST) !NOT USED IN THIS VERSION DIMENSION ITAR(MZCHF) DIMENSION IFN(MZCHF),MHIGH(MZCHF),MLOW(MZCHF) DIMENSION NOT(MZEST),ID(MZEST),IDDD(MZEST) C DATA IDNUM/'00','01','02','03','04','05','06','07','08','09', 1 '10','11','12','13','14','15','16','17','18','19', 2 '20','21','22','23','24','25','26','27','28','29', 3 '30','31','32','33','34','35','36','37','38','39', 4 '40','41','42','43','44','45','46','47','48','49', 5 '50','51','52','53','54','55','56','57','58','59', 6 '60','61','62','63','64','65','66','67','68','69', 7 '70','71','72','73','74','75','76','77','78','79', 8 '80','81','82','83','84','85','86','87','88','89', 9 '90','91','92','93','94','95','96','97','98','99'/ DATA IDNUM1/'0','1','2','3','4','5','6','7','8','9'/ DATA IPNT/'.'/ DATA EXCLAM/' ! '/ ccc DATA SMALL/1.0E-01/ reduced to avoid '!' by KAB, Apr 1997 DATA SMALL/1.0E-02/ DATA WTEST/5.0E+00/ C CARBON SEQUENCE USED WTEST=10., FEV TO FEII USED WTEST=5. C IDSEQ=IDNUM(NELC) IDSEQ1=IDNUM(NELC+1) C C THE CHANNELS MUST FIRST BE IDENTIFIED IN TERMS OF THE C ID TARGET STATES. THESE CAN BE IN A SLIGHTLY DIFFERENT ORDER C FROM THE TARGET STATES READ FROM FILE H WHICH ARE IN ENERGY C ORDER. C I=0 DO 1 NST=1,NASTR K=NCNATR(NST) IF (K .EQ. 0) GO TO 1 DO 2 L=1,K I=I+1 ITAR(I)=NST IPAT(NST)=MOD(LLCH(I)+ISLP,2) 2 CONTINUE ITASLP(NST)=10000*ISAT(NST)+100*LAT(NST)+IPAT(NST) 1 CONTINUE CW WRITE(6,*)NASTR,(ITASLP(I),I=1,NASTR) C IF (IDNST .NE. 0) THEN DO 3 IDN=1,IDNST IDTA=IDSLP(IDN) DO 4 NST=1,NASTR IF (IDTA .EQ. ITASLP(NST)) THEN IDTEMP(NST)=IDN ITASLP(NST)=0 GO TO 3 END IF 4 CONTINUE 3 CONTINUE CW WRITE(6,*)(IDTEMP(I),I=1,NASTR) C C STORE THE TARGET IDENTIFICATION FOR EACH CHANNEL C DO 5 J=1,NCHF IDT=IDTEMP(ITAR(J)) IDTARG(J)=IDNUM(IDT) 5 CONTINUE CW WRITE(6,*)(IDTARG(J),J=1,NCHF) ELSE DO 52 J=1,NCHF IDT=ITAR(J) IDTARG(J)=IDNUM(IDT) 52 CONTINUE END IF C DO 6 M=1,MXE JID(M)=0 IDDD(M)=0 6 CONTINUE C C IDENTIFY LOWEST LEVELS BELONGING TO NEXT HIGHER SEQUENCE C M1=0 DO 7 I1=1,IDNST1 IF (ISLP .EQ. IDSLP1(I1)) THEN M1=M1+1 ID(M1)=IDSEQ1//IDNUM(I1)//IPNT//IDNUM(0)//IDNUM(0) CW WRITE(6,*)M1,ID(M1) END IF 7 CONTINUE MSTART=M1+1 C C FORM A WORKING COPY OF THE WEIGHT TABLE C DO 8 M=MSTART,MXE DO 9 J=1,NCHF WTEMP(J,M)=WT(J,M) 9 CONTINUE 8 CONTINUE C C START MAIN IDENTIFICATION LOOPS. START WITH FIRST CHANNEL. C DO 10 J=1,NCHF C C C START AT HIGHEST LEVEL AND RECORD ALL WEIGHTS GT. WTEST C K=0 DO 11 M=MXE,MSTART,-1 IF (WTEMP(J,M) .GE.WTEST) THEN K=K+1 FN(K)=FNST(J,M) MTEMP(K)=M END IF 11 CONTINUE KLAST=K C C COLLECT MORE FNU WITH WT LT. WTEST BY CHOOSING LEVELS WHERE THIS C CHANNEL HAS WEIGHT GT. TWICE THAT OF HIGHER CHANNELS STARTING C AFTER THE LAST RECORDED LEVEL. C IGNORE A LEVEL WHERE THE WEIGHT IS VERY SMALL AS IT IS EITHER C ALREADY ASSIGNED OR SHOULD HAVE BEEN ASSIGNED TO AN N+1 ELECTRON C BOUND STATE. C IF (KLAST .GT. 0) THEN MLAST=MTEMP(KLAST)-1 ELSE MLAST=MXE END IF DO 12 M=MLAST,MSTART,-1 WTEM1=WTEMP(J,M)*2. IF (WTEM1 .LE. SMALL) GO TO 12 DO 13 J1=J+1,NCHF IF (WTEM1 .LT. WTEMP(J1,M)) GO TO 12 13 CONTINUE K=K+1 FN(K)=FNST(J,M) MTEMP(K)=M 12 CONTINUE KLAST=K CW WRITE(6,*)KLAST,(FN(I),I=1,KLAST) C C SEARCH FOR DUPLICATED FNU AND CHECK FOR POSSIBLY MISSING FNU C FORM ARRAY NDUP(N) THE NUMBER OF DUPLICATED FNU WHERE N COUNTS C THE NUMBER OF DIFFERENT FNU C IF (KLAST .EQ. 0) GO TO 10 F1=FN(1) N=1 NDUP(N)=1 DO 15 K=2,KLAST F2=FN(K) IF (F1-F2 .GT. .8) THEN IF (F1-F2 .GT. 1.5) THEN WRITE(6,'('' POSSIBLE MISSING LEVEL FOR CHANNEL'' 1 ,I3,''BETWEEN FNU='',F5.3,''AND'',F5.3)')J,F2,F1 END IF N=N+1 NDUP(N)=1 F1=F2 ELSE NDUP(N)=NDUP(N)+1 END IF 15 CONTINUE C C CONSTRUCT ID FILE FOR CHANNEL J. C SELECT MOST LIKELY LEVEL WHERE DUPLICATION OCCURS BY A PROCESS C OF ELIMINATION SO THAT THE MOST LIKELY LEVEL IS NAMED IN MDUP(1) C WHERE INITIALLY THE LENGTH OF ARRAY MDUP IS NDUP(N) FOR NTH FNU C NTOTAL=N CW WRITE(6,*)NTOTAL,(NDUP(N),N=1,NTOTAL) KSTART=1 DO 16 N=1,NTOTAL IDT=0 MDUP(1)=MTEMP(KSTART) IF (NDUP(N) .GT. 1) THEN ND1=NDUP(N) C C REDUCE DUPLICATION BY ELIMINATING LEVELS WHERE TWICE WT FOR THIS C IS NOT MAX. CF HIGHER CHANNELS. ND2 WILL BE NEW NUMBER OF LEVELS. C K=KSTART ND2=0 DO 17 ND=1,ND1 M=MTEMP(K) WTEM1=WTEMP(J,M)*2. DO 18 J1=J+1,NCHF IF (WTEM1 .LT. WTEMP(J1,M)) GO TO 19 18 CONTINUE ND2=ND2+1 MDUP(ND2)=M 19 K=K+1 17 CONTINUE C C IF STILL DUPLICATION CHECK WHETHER ONE OF THE LEVELS C BELONGS TO A LOWER CHANNEL BY LOOKING AT IDOUBT. C IF NO LEVEL LEFT KEEP ND2 LEVELS C **** THIS FACILITY NOT IN USE IN THIS VERSION C IF (ND2 .EQ. 0) THEN K=KSTART DO 91 ND=1,ND1 M=MTEMP(K) ND2=ND2+1 MDUP(ND2)=M K=K+1 91 CONTINUE END IF IF (ND2 .GT. 1) THEN ND3=0 DO 20 ND=1,ND2 M=MDUP(ND) ND3=ND3+1 MDUP(ND3)=M 20 CONTINUE C C IF STILL DUPLICATION ELIMINATE ANY LEVELS WHERE THIS CHANNEL C WEIGHT IS LESS THAN HALF OF THE MAX AMONGST THE DUPLICATING C LEVELS C IF (ND3 .EQ. 0) MDUP(1)=0 IF (ND3 .GT. 1) THEN WMAX=0. DO 23 ND=1,ND3 M=MDUP(ND) IF (WMAX .LT. WTEMP(J,M)) THEN WMAX=WTEMP(J,M) MMAX=M END IF 23 CONTINUE HALFWM=.5*WMAX ND4=0 DO 24 ND=1,ND3 M=MDUP(ND) IF (WTEMP(J,M) .GT. HALFWM) THEN ND4=ND4+1 MDUP(ND4)=M END IF 24 CONTINUE CW WRITE(6,*)' ND4',ND4,(MDUP(NN),NN=1,ND4) C C IF STILL DUPLICATION FIND THE CHANNELS WITH WEIGHT GREATER C THAN HALF THE MAX. WEIGHT OF THE CHANNELS ABOVE J. C IS THERE A COMMON PERTURBER FOR 2 OR MORE LEVELS? C TRY TO SELECT ONE OF THESE LEVELS ACCORDING TO CHANNEL L C OR CHANNEL TARGET OTHERWISE JUST CHOOSE LEVEL WITH HIGHEST C WEIGHT C IF (ND4 .GT. 1) THEN DO 25 J1=J+1,NCHF NJMAX(J1)=0 25 CONTINUE DO 26 ND=1,ND4 M=MDUP(ND) WTEM1=WTEMP(J+1,M) JM=J+1 DO 27 J1=J+2,NCHF IF (WTEM1 .LT. WTEMP(J1,M)) THEN WTEM1=WTEMP(J1,M) JM=J1 END IF 27 CONTINUE WTEMH=.5*WTEM1 DO 61 J1=J+1,NCHF IF (WTEMH .LT. WTEMP(J1,M)) THEN NJMAX(J1)=NJMAX(J1)+1 NDJM(J1,NJMAX(J1))=M END IF 61 CONTINUE 26 CONTINUE NDJMAX=1 JM=0 DO 28 J1=J+1,NCHF IF (NJMAX(J1) .GT. NDJMAX) THEN NDJMAX =NJMAX(J1) JM=J1 END IF 28 CONTINUE C IF (NDJMAX .GT. 1) THEN DO 90 ND=1,NDJMAX M=NDJM(JM,ND) MDUP(ND)=M 90 CONTINUE M=NDJM(JM,1) IF ((FNST(J,M)-FNST(JM,M)) .LT. .8) THEN IF (L2P(J) .EQ. L2P(JM) .OR. 2 IDTARG(J) .EQ. IDTARG(JM)) THEN DO 50 ND=2,NDJMAX MDUP(ND-1)=MDUP(ND) 50 CONTINUE NDJMAX=NDJMAX-1 END IF END IF IF (NDJMAX .GT. 1) THEN M=MDUP(1) WTEM1=WTEMP(J,M) DO 51 ND=2,NDJMAX M=MDUP(ND) WTEM2=WTEMP(J,M) IF (WTEM1 .LT. WTEM2) THEN WTEM1=WTEM2 MDUP(1)=M END IF 51 CONTINUE END IF IF (MDUP(1) .NE. MMAX) IDT=5 ELSE MDUP(1)=MMAX IDT=5 END IF END IF END IF END IF END IF C C STORE THE ID C C IF MDUP(1)=0 NO LEVEL IDENTIFIED IN THIS LOOP C C M=MDUP(1) IF (M .GT. 0) THEN JID(M)=J IDD(M)=IDT DO 35 J1=1,NCHF WTEMP(J1,M)=0. 35 CONTINUE END IF 30 KSTART=KSTART+NDUP(N) 16 CONTINUE 10 CONTINUE C C NOW SHOULD HAVE JID(M) FOR ALL M C FILL ID FILE C C ZEROISE THE MLOW(J) ARRAY WHICH STORES LOWEST LEVEL FOR C CHANNEL J WHICH HAS BEEN IDENTIFIED C AND THE MHIGH(J) ARRAY WHICH STORES HIGHEST LEVEL FOR C CHANNEL J C NEWTRY=0 42 DO 31 J=1,NCHF MLOW(J)=0 MHIGH(J)=0 31 CONTINUE C NOTID=0 DO 32 M=MSTART,MXE IF (JID(M) .NE. 0) THEN J=JID(M) IDDD(M)=J IF (MLOW(J) .EQ. 0) THEN MLOW(J)=M MHIGH(J)=M L=L2P(J) IFNJ=NCORE(L+1)+1 IFN(J)=IFNJ DEFECT(J)=IFNJ-FNST(J,M) ID(M)=IDSEQ//IDTARG(J)//IPNT//IDNUM(IFNJ)//IDNUM1(L) 1 //IDNUM1(IDD(M)) ELSE L=L2P(J) IFNJ=IFN(J)+1 NOTE=0 33 DEF=IFNJ-FNST(J,M) IF (DEFECT(J)-DEF .GT. .8) THEN NOTE=1 IFNJ=IFNJ+1 WRITE(6,'('' LEVEL MISSING FOR CHANNEL'',I3)')J GO TO 33 END IF IF (ABS(DEFECT(J)-DEF) .GT. .2)IDD(M)=5 ID(M)=IDSEQ//IDTARG(J)//IPNT//IDNUM(IFNJ)//IDNUM1(L) 1 //IDNUM1(IDD(M)) IFN(J)=IFNJ IF (NOTE .EQ. 0) THEN DO 36 MA=MHIGH(J),M WTEMP(J,MA)=0. 36 CONTINUE END IF MHIGH(J)=M END IF ELSE NOTID=NOTID+1 NOT(NOTID)=M ID(M)=EXCLAM END IF 32 CONTINUE C C C IF ANY LEVELS WERE NOT IDENTIFIED HAVE ANOTHER TRY C IF (NEWTRY .EQ. 1 .OR. NOTID .EQ. 0) GO TO 43 DO 37 N=1,NOTID M=NOT(NOTID) 39 JMAX=0 WMAX=SMALL DO 38 J=1,NCHF IF (WMAX .LT. WTEMP(J,M)) THEN WMAX=WTEMP(J,M) JMAX=J END IF 38 CONTINUE IF (JMAX .GT. 0) THEN M1=MLOW(J) M2=MHIGH(J) IF (M .GT. M2) THEN IF (M2 .GT. 0) THEN DEL=FNST(JMAX,M)-FNST(JMAX,M2) IF (DEL .LT. .9 .OR. DEL .GT. 1.2) THEN WTEMP(JMAX,M)=0. GO TO 39 END IF MMB=M MMA=MHIGH(JMAX) MHIGH(JMAX)=M ELSE MMA=0 MMB=0 MLOW(JMAX)=M MHIGH(JMAX)=M END IF ELSE IF (M .LT. M1) THEN DEL=FNST(JMAX,M1)-FNST(JMAX,M) IF (DEL .LT. .9 .OR. DEL .GT. 1.2) THEN WTEMP(JMAX,M)=0. GO TO 39 END IF MMA=M MMB=MLOW(JMAX) MLOW(JMAX)=M END IF JID(M)=JMAX IDD(M)=4 DO 40 JB=1,NCHF WTEMP(JB,M)=0. 40 CONTINUE DO 41 MB=MMA,MMB WTEMP(JMAX,MB)=0. 41 CONTINUE END IF 37 CONTINUE NEWTRY=1 GO TO 42 C WRITE INFORMATION FOR IDENTIFICATIONS 43 WRITE(6,660)(M,M=1,NCHF) WRITE(6,662) DO 71 M=1,MXE 71 WRITE(6,661)M,(WT(J,M),J=1,NCHF) WRITE(6,663)(IDTARG(J),L2P(J),J=1,NCHF) WRITE(6,662) DO 72 M=1,MXE 72 WRITE(6,664)M,ID(M),(FNST(J,M),J=1,NCHF) 660 FORMAT(//6X,'CHANNEL PERCENTAGE WEIGHTS FOR CHANNEL NUMBER = '/ 1 (6X,8I8)) 661 FORMAT(1X,I3,3X,8(2PE8.1)/7X,8E8.1/7X,8E8.1) 662 FORMAT(/) 663 FORMAT(//6X,'EFFECTIVE QUANTUM NUMBERS, CHANNEL ID = '/ 1 (16X,8(2X,A2,'L=',I2))) 664 FORMAT(1X,I3,2X,A9,2X,8F7.3/(17X,8F7.3)/ 1 (17X,8F7.3)) RETURN END C C************************************************************************ C SUBROUTINE IDFORM(ITARG1,KSLP,ENAT,NAST,ISAT,LAT,IDNST1) C C CONVERT ID FILE (11) TO E FILE (12) BY KAB, APRIL 1997 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CIDFOR/IPAT(MZTAR) DIMENSION ITARG1(KSLP),ENAT(NAST),SAMPLE(10),ISAT(NAST),LAT(NAST) CHARACTER*1 DOT,COLON DATA KOUNT/0/ C OPEN(12,FILE='E',FORM='FORMATTED',STATUS='UNKNOWN') C REWIND 11 C READ(11,1000,END=9) NZ,NE,NS C NZNE=NZ*100+NE ZSQ=(MAX(NZ-NE,1))**2 WRITE(12,1001) NZ,NE,'E' C OPEN(13,FILE='T',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(13,1001) NZ,NE,'T' WRITE(13,1001) NAST DO M=1,NAST WRITE(13,1013) ISAT(M),LAT(M),IPAT(M),ZSQ*ENAT(M),M ENDDO WRITE(13,1001) IDNST1 C WRITE(6,1001) NZ,NE,'E' KOUNT=2 KASE=0 C 1 READ(11,*,END=10) IS,IL,IP IF(IL.EQ.-1) GO TO 10 READ(11,*) NLEV WRITE(12,1010) IS,IL,IP,NLEV KASE=KASE+1 ECH1=ENAT(ITARG1(KASE)) WRITE(12,1005) ITARG1(KASE),ZSQ*ECH1,NAST KOUNT=KOUNT+NLEV+2 C DO 5 M=1,NLEV READ(11,1002,END=9,ERR=9) I,NELC,IT,DOT,N,L,IC,ESCALE,IDDD RYD=ESCALE*ZSQ IF(ECH1.LE.ESCALE) GO TO 9 EFFN=1.0/SQRT(ECH1-ESCALE) IF(DOT.EQ.'!') IC=9 COLON=' ' IF(IC.NE.0) COLON=':' IF(DOT.EQ.'!') THEN ISAMPL=0 DO 3 N=ITARG1(KASE),NAST IF(ENAT(N).LE.ESCALE) GO TO 4 ISAMPL=ISAMPL+1 SAMPLE(ISAMPL)=1.0/SQRT(ENAT(N)-ESCALE) IF(ISAMPL.EQ.10) GO TO 4 3 CONTINUE 4 WRITE(12,1003) I,'!',0,0,0,COLON,RYD,(SAMPLE(N),N=1,ISAMPL) ELSE IF(NELC.EQ.NE) THEN IF(ENAT(IT).LE.ESCALE) GO TO 9 EFFM=1.0/SQRT(ENAT(IT)-ESCALE) WRITE(12,1113) I,'T',IT,N,L,COLON,RYD,EFFN,EFFM,IDDD ELSE IF(NELC.EQ.NE+1) THEN WRITE(12,1003) I,'C',IT,0,0,COLON,RYD,EFFN ELSE GO TO 9 ENDIF 5 CONTINUE GO TO 1 C 9 WRITE(6,*)'ERROR - DATA PROBLEM:',I,NELC,IT,DOT,N,L,IC,KOUNT STOP 10 WRITE(12,1010) 0,-1,0,0 WRITE(6,*)'E FILE COMPLETE: ',KOUNT,' RECORDS.' RETURN C 1000 FORMAT(1X,2I2,I4) 1001 FORMAT(2I5,4X,A1) 1002 FORMAT(I4,2X,2I2,A1,I2,2I1,E14.6,I5) 1003 FORMAT(I5,2X,A1,2I3,I2,A1,1PE13.5,0P10F10.5) 1005 FORMAT(I5,E14.6,I6,'T') 1010 FORMAT(4I5) 1013 FORMAT(3I5,1PE13.5,3X,'T',I3) 1113 FORMAT(I5,2X,A1,2I3,I2,A1,1PE13.5,0P2F10.5,I5) END C*************************************************************** C SUBROUTINE INTBUT(KK,XX,BUTTLE) C C DETRMINE BUTTLE CORRECTION FROM DARC C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C C FX(X)=((X-X2)/(X1-X2))*((X-X3)/(X1-X3))*Y1+((X-X1)/(X2-X1))* X ((X-X3)/(X2-X3))*Y2+((X-X1)/(X3-X1))*((X-X2)/(X3-X2))*Y3 C MM=NBUTD(KK) C IF (MM.EQ.2) THEN X1=EBUTD(1,KK) X2=EBUTD(2,KK) Y1=CBUTD(1,KK) Y2=CBUTD(2,KK) BUTTLE=(XX-X2)/(X1-X2)*Y1+(XX-X1)/(X2-X1)*Y2 RETURN ENDIF C DO M=1,MM IF(XX.LT.EBUTD(M,KK))GO TO 10 ENDDO II=MM GO TO 20 C 10 CONTINUE II=M 20 CONTINUE C IF (II.EQ.1) THEN X1=EBUTD(1,KK) X2=EBUTD(2,KK) X3=EBUTD(3,KK) Y1=CBUTD(1,KK) Y2=CBUTD(2,KK) Y3=CBUTD(3,KK) BUTTLE=FX(XX) RETURN ENDIF C IF (II.EQ.MM) THEN X1=EBUTD(MM-2,KK) X2=EBUTD(MM-1,KK) X3=EBUTD(MM ,KK) Y1=CBUTD(MM-2,KK) Y2=CBUTD(MM-1,KK) Y3=CBUTD(MM ,KK) BUTTLE=FX(XX) RETURN ENDIF C X1=EBUTD(II-1,KK) X2=EBUTD(II ,KK) X3=EBUTD(II+1,KK) Y1=CBUTD(II-1,KK) Y2=CBUTD(II ,KK) Y3=CBUTD(II+1,KK) BUTTLE=FX(XX) RETURN C END C*************************************************************** C SUBROUTINE INTBUTD(KK,XX,BUTTLED) C C DETRMINE ENERGY DERIVATIVE OF BUTTLE CORRECTION FROM DARC C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C C FXD(X)=((X-X2+X-X3)/((X1-X2)*(X1-X3)))*Y1+((X-X1+X-X3)/((X2-X1)* X (X2-X3)))*Y2+((X-X1+X-X2)/((X3-X1)*(X3-X2)))*Y3 C MM=NBUTD(KK) C IF (MM.EQ.2) THEN X1=EBUTD(1,KK) X2=EBUTD(2,KK) Y1=CBUTD(1,KK) Y2=CBUTD(2,KK) BUTTLED=(Y2-Y1)/(X2-X1) RETURN ENDIF C DO M=1,MM IF(XX.LT.EBUTD(M,KK))GO TO 10 ENDDO II=MM GO TO 20 C 10 CONTINUE II=M 20 CONTINUE C IF (II.EQ.1) THEN X1=EBUTD(1,KK) X2=EBUTD(2,KK) X3=EBUTD(3,KK) Y1=CBUTD(1,KK) Y2=CBUTD(2,KK) Y3=CBUTD(3,KK) BUTTLED=FXD(XX) RETURN ENDIF C IF (II.EQ.MM) THEN X1=EBUTD(MM-2,KK) X2=EBUTD(MM-1,KK) X3=EBUTD(MM ,KK) Y1=CBUTD(MM-2,KK) Y2=CBUTD(MM-1,KK) Y3=CBUTD(MM ,KK) BUTTLED=FXD(XX) RETURN ENDIF C X1=EBUTD(II-1,KK) X2=EBUTD(II ,KK) X3=EBUTD(II+1,KK) Y1=CBUTD(II-1,KK) Y2=CBUTD(II ,KK) Y3=CBUTD(II+1,KK) BUTTLED=FXD(XX) RETURN C END C C********************************************************************* C SUBROUTINE LU(A,N,IR,JC) C C LOWER-UPPER TRIANGULAR DECOMPOSITION OF MATRIX A, C WITH FULL PIVOTING. C C A = L * U C WHERE L IS UNIT LOWER TRIANGULAR AND U UPPER TRIANGULAR. C C RETURNS L BELOW DIAGONAL IN A, U ON AND ABOVE DIAGONAL IN A. C IR AND JC SPECIFY ROW AND COLUMN INTERCHANGES. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C DIMENSION A(MZCHF,MZCHF),IR(MZCHF),JC(MZCHF) C C C INITIALISE IR,JC C DO K=1,N IR(K)=K JC(K)=K ENDDO C C FIRST PIVOT AND FIRST INTERCHANGES. C M=1 CALL PIV(A,N,M,IR,JC) C C FIRST COLUMN OF L C T=1./A(1,1) DO I=2,N A(I,1)=A(I,1)*T ENDDO C C MAIN LOOP C DO M=2,N C C NEW PIVOT AND INTERCHANGES C CALL PIV(A,N,M,IR,JC) C C CONTINUE CALCULATION OF L AND U C C - ROW OF U C DO J=M,N DO K=1,M-1 A(M,J)=A(M,J)-A(M,K)*A(K,J) ENDDO ENDDO C C - COLUMN OF L C DO K=1,M-1 DO I=M+1,N A(I,M)=A(I,M)-A(I,K)*A(K,M) ENDDO ENDDO T=1./A(M,M) DO I=M+1,N A(I,M)=A(I,M)*T ENDDO ENDDO C RETURN END C C*************************************************************** C SUBROUTINE MATMUL(A,B,C,MZ,N) C C CALCULATES C=A*B C C NRB: VECTOR OR SCALAR VERSIONS C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (ZERO=0.0) C DIMENSION A(MZ,*),B(MZ,*),C(MZ,*) C CCRAY CRAY DO J=1,N CRAY DO I=1,N CRAY C(I,J)=ZERO CRAY ENDDO CRAY ENDDO CRAY DO K=1,N CRAY DO J=1,N CRAY DO I=1,N CRAY C(I,J)=C(I,J)+A(I,K)*B(K,J) CRAY ENDDO CRAY ENDDO CRAY ENDDO CCRAY C---- DO J=1,N DO I=1,N C(I,J)=ZERO ENDDO DO K=1,N DO I=1,N C(I,J)=C(I,J)+A(I,K)*B(K,J) ENDDO ENDDO ENDDO C---- RETURN END C C********************************************************************* C SUBROUTINE MTMLT(A,B,C,MZ,N) C C CALCULATES C=A(TRANSPOSE)*B C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION A(MZ,*),B(MZ,*),C(MZ,*) C DO J=1,N DO I=1,N S=0. DO K=1,N S=S+A(K,I)*B(K,J) ENDDO C(I,J)=S ENDDO ENDDO C RETURN END C C********************************************************************* C SUBROUTINE MTV(A,B,C,MZ,N) C C CALCULATES C=A*B WHERE A IS A SQUARE MATRIX AND C B AND C ARE COLUMN VECTORS. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C DIMENSION A(MZ,*),B(*),C(*) C DO I=1,N C(I)=0. ENDDO DO J=1,N DO I=1,N C(I)=C(I)+A(I,J)*B(J) ENDDO ENDDO C RETURN END C C********************************************************************* C SUBROUTINE NUMT(E,C,R1,HP,N1,N2,I) C C NUMEROV INTEGRATION OF COULOMB FUNCTIONS. C THE INTERVAL HP IS POSITIVE. C INTEGRATION FROM TABULAR POINT N1 TO TABULAR POINT N2. C FUNCTIONS THETA,THETP STORED IN FS,FSP C FUNCTIONS THETAD,THETADP STORED IN FC,FCP C STARTS WITH FUNCTIONS AND DERIVATIVES AT N1 STORED IN FS, FSP C CALCULATES FUNCTIONS AT ALL POINT TO N2 AND DERIVATIVE AT C THE POINT N2. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/COULSC/FS(MZCHF,MZPTS),FSP(MZCHF),FC(MZCHF,MZPTS) 1 ,FCP(MZCHF) COMMON/CTHET/BB(MZCHF,MZTET),BG(MZCHF,MZTET),MSUM(MZCHF) COMMON/NRBZED/TZED C V(X)=EQ+X*(Q2*TZED-X*CQ) C N21=N2-N1 C C RENORMALISE FOR CASE OF N2.EQ.N1 IF(N21.NE.0)GOTO 5 W=FS(I,N1)*FCP(I)-FSP(I)*FC(I,N1) IF(W.EQ.0.0)RETURN W1=1./SQRT(W) FS(I,N1)=FS(I,N1)*W1 FC(I,N1)=FC(I,N1)*W1 FSP(I)=FSP(I)*W1 FCP(I)=FCP(I)*W1 BB(I,2)=BB(I,2)*W1 NM=MSUM(I) DO 2 M=3,NM BB(I,M)=BB(I,M)*W1 2 BG(I,M)=BG(I,M)*W1 RETURN C C INTEGRATIONS FOR N2.NE.N1 5 IP=IABS(N21) IS=N21/IP H=HP*IS K=N1 IP=IP-1 C F1=FS(I,N1) F1P=FSP(I) G1=FC(I,N1) G1P=FCP(I) C C FUNCTIONS AT K=(N1+IS) Q=H*H EQ=E*Q Q2=2.*Q CQ=C*Q C X1=1./R1 CX=C*X1 HX=H*X1 A1=-2.*HX*HX*H U1P=A1*(1.*TZED-CX) A1=-A1*HX U1PP=A1*(2.*TZED-3.*CX) A1=-A1*HX*6. U1PPP=A1*(1.*TZED-2.*CX) U1=V(X1) C R2=R1+H X2=1./R2 U2=V(X2) C A2=1.+U2*.033333333 B1=1.+U1*(.025*U1-.46666667)-.13333333*U1P-.025*U1PP 1 +.0027777778*(4.*U1*U1P-U1PPP) C1=H*(1.+U1*(.0027777778*U1-.13333333) 1 -.05*U1P-.0083333333*U1PP) P2=.033333333*Q D1=Q*(-.46666667+.025*U1+.025*U1+.0027777778*U1P) E1=H*Q*(-.13333333+.0027777778*U1) C F2=(B1*F1+C1*F1P)/A2 G2=(B1*G1+C1*G1P+D1*F1+E1*F1P-P2*F2)/A2 K=K+IS FS(I,K)=F2 FC(I,K)=G2 C C U3=U2 U2=U1 F3=F2 F2=F1 G3=G2 G2=G1 IF(IP.EQ.0)GOTO 20 C C CONTINUE INTEGRATION EQ=EQ*.083333333 Q2=Q2*.083333333 CQ=CQ*.083333333 U2=U2*.083333333 U3=U3*.083333333 Q=Q*.083333333333 R3=R2 C DO 10 M=1,IP U1=U2 U2=U3 F1=F2 F2=F3 G1=G2 G2=G3 R3=R3+H X3=1./R3 U3=V(X3) D3=1./(1.+U3) D2=(2.-10.*U2)*D3 D1=(1.+U1)*D3 F3=D2*F2-D1*F1 G3=D2*G2-D1*G1-Q*D3*(F3+10.*F2+F1) K=K+IS FC(I,K)=G3 10 FS(I,K)=F3 C U2=12.*U2 U3=12.*U3 EQ=12.*EQ Q2=12.*Q2 CQ=12.*CQ Q=12.*Q C C CALCULATE FINAL DERIVATIVE 20 H=-H CX=C*X3 HX=H*X3 A1=-2.*HX*HX*H U3P=A1*(1.*TZED-CX) A1=-A1*HX U3PP=A1*(2.*TZED-3.*CX) A1=-A1*HX*6. U3PPP=A1*(1.*TZED-2.*CX) A2=1.+U2*.033333333 B3=1.+U3*(.025*U3-.46666667)-.13333333*U3P-.025*U3PP 1 +.0027777778*(4.*U3*U3P-U3PPP) C3=H*(1.+U3*(.0027777778*U3-.13333333)-.05*U3P 1 -.0083333333*U3PP) P2=.033333333*Q D3=Q*(-.466666666+.025*U3+.0027777778*U3P) E3=H*Q*(-.13333333+.0027777777*U3) F3P=(A2*F2-B3*F3)/C3 G3P=(A2*G2-B3*G3+P2*F2-D3*F3-E3*F3P)/C3 C C RE-NORMALISE CLOSED-CHANNEL FUNCTIONS C == MODIFICATION FOR VAX CALCULATION OF C == W1=1./SQRT(F3*G3P-F3P*G3) AMAX=1./MAX(ABS(F3),ABS(G3),ABS(F3P),ABS(G3P)) AF3=F3*AMAX AG3=G3*AMAX AF3P=F3P*AMAX AG3P=G3P*AMAX W1=AMAX/SQRT(AF3*AG3P-AF3P*AG3) C == END MODIFICATION FSP(I)=F3P*W1 FCP(I)=G3P*W1 IP=IP+2 DO 30 J=1,IP FS(I,J)=FS(I,J)*W1 30 FC(I,J)=FC(I,J)*W1 C RE-NORMALISE COEFFICIENTS BB(I,2)=BB(I,2)*W1 NM=MSUM(I) DO 40 M=3,NM BB(I,M)=BB(I,M)*W1 40 BG(I,M)=BG(I,M)*W1 C RETURN END C C********************************************************************* C SUBROUTINE PIV(A,N,M,IR,JC) C C PIVOTING OF A FOR I=M,N AND J=M,N C C RETURN A WITH INTERCHANGES AND INTERCHANGE INFORMATION C IN IR AND JC C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C DIMENSION A(MZCHF,MZCHF),IR(MZCHF),JC(MZCHF) C C C FIND PIVOT ELEMENT C P=0. DO J=M,N DO I=M,N Q=ABS(A(I,J)) IF(Q.GT.P) THEN P=Q IP=I JP=J ENDIF ENDDO ENDDO C C NEW IR,JC C I=IR(M) IR(M)=IR(IP) IR(IP)=I J=JC(M) JC(M)=JC(JP) JC(JP)=J C C PERFORM INTERCHANGES C IF(IP.NE.M) THEN DO J=1,N S=A(M,J) A(M,J)=A(IP,J) A(IP,J)=S ENDDO ENDIF C IF(JP.NE.M) THEN DO I=1,N S=A(I,M) A(I,M)=A(I,JP) A(I,JP)=S ENDDO ENDIF C RETURN END C C********************************************************************* C SUBROUTINE POINTB C C CALCULATES CHANNELS ENERGIES, NUMBER OF OPEN CHANNELS C AND TABULAR POINTS C NRB: C NEUTRAL CASE ADDED C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C LOGICAL WARN COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CNTRL/IPRINT,IRAD,IPERT COMMON/CWARN/WARN COMMON/NRBZED/TZED C KP0=1 C C FOR IRAD.EQ.1, HOLD OLD KP2, RTWO AND H IF(IRAD.GT.0)THEN KP2OLD=KP2 RTWOLD=RTWO HOLD=H ENDIF C C CHANNEL ENERGIES EPS AND NUMBER OF OPEN CHANNELS NCHOP IF(IPRINT.GT.2)WRITE(6,630) NCHOP=0 DO 100 I=1,NCHF E=ETOT-ECH(I) IF(E.GE.0.)THEN NCHOP=NCHOP+1 ELSE FKNU(I)=1./SQRT(-E) ENDIF EPS(I)=E IF(IPRINT.GT.2)WRITE(6,640)I,EPS(I),FKNU(I) 100 CONTINUE C C STOP IF NCHOP.GT.0 IF(NCHOP.GT.0)THEN WRITE(6,620)ETOT,NCHOP STOP ENDIF NCHOP1=NCHOP+1 C C CALCULATION OF RTWO C INCLUDES CALCULATION OF INNER POINTS OF INFLECTION RINF C FOR CLOSED CHANNELS RTWO IS EQUAL TO THE OUTER POINT OF C INFLECTION, EXCEPT FOR THE CASE OF FNU.LT.(LL+1) RTWO=RZERO DO 150 I=NCHOP1,NCHF IF(TZED.EQ.0)THEN RINF(I)=0. C R2=RZERO ELSE FNU=FKNU(I) FLP1=LLCH(I)+1 IF(FNU.GE.FLP1)GOTO 145 R2=FLP1*(FLP1+SQRT(FLP1)) RINF(I)=0. IF(R2.LT.10.)R2=10. IF(RTWO.LT.R2)RTWO=R2 GOTO 155 145 A1=SQRT(FNU*FNU-CC(I)) R2=FNU*(FNU+A1) RINF(I)=FNU*(FNU-A1) IF(RTWO.LT.R2)RTWO=R2 155 IF(IPRINT.GT.2)WRITE(6,640)I,R2 ENDIF 150 CONTINUE IF(RONE.GT.RTWO)RTWO=RONE C C CASE OF IRAD.EQ.1 IF(IRAD.GT.0)THEN IF(IPRINT.GT.1)WRITE(6,650)ETOT,RTWO C CHECK VALUE OF RTWO IF(RTWO.GT.RTWOLD)THEN IF(WARN)WRITE(6,610)ETOT,RTWOLD WARN=.FALSE. ENDIF C RE-INSTATE OLD KP2,RTWO AND H KP2=KP2OLD RTWO=RTWOLD H=HOLD RETURN ENDIF C C FIND INTERVAL AND TABULAR POINTS BETWEEN RZERO AND RTWO IF(RTWO.LE.RZERO)THEN KP2=1 H=0. RETURN ENDIF WM=0. DO 170 I=1,NCHF C=CC(I) E=EPS(I) X=1./RZERO W=ABS(E+X*(2.*TZED-C*X)) IF(W.GT.WM)WM=W X=1./RTWO W=ABS(E+X*(2.*TZED-C*X)) IF(W.GT.WM)WM=W IF(C.LT.RZERO.OR.C.GT.RTWO)GOTO 170 W=ABS(E+TZED/C) IF(W.GT.WM)WM=W 170 CONTINUE C H=ACNUM/SQRT(WM) C..... MODIFICATION 2.2.87 IF(ABS(RTWO-RZERO).LT.(2.*H))THEN RTWO=RZERO KP2=1 H=0. RETURN ENDIF C..... END MODIFICATION N=(RTWO-RZERO)/H N=4*((N-1)/4)+5 KP2=N C CHECK THAT DIMENSIONS NOT EXCEEDED IF(KP2.GT.MZPTS)THEN KP2=4*((MZPTS-1)/4)+1 RTWO=(KP2-1)*H+RZERO IF(WARN)WRITE(6,600)ETOT,KP2 WARN=.FALSE. ENDIF H=KP2-1 H=(RTWO-RZERO)/H C C WEIGHTS FOR BODE RULE INTEGRATION CALL BODE C 600 FORMAT(//3X,' ** WARNING ** , FOR ETOT = ', 1 1PE11.4,' KP2 REDUCED TO MAXIMUM VALUE OF ',I4, 2 3X,' ** NO MORE WARNINGS GIVEN **'//) 610 FORMAT(//1X,' ** WARNING ** , FOR ETOT=',E13.5, 1 ' RTWO REDUCED TO MAXIMUM VALUE OF ',F6.0, 2 1X,' ** NO MORE WARNINGS GIVEN **'//) 620 FORMAT(//5X,20('*')/5X,'ETOT = ',F10.4,'GIVES ',I3,' OPEN', + ' CHANNELS'/5X,20('*')/) 630 FORMAT(/3X,'I',4X,'EPS(I)',3X,'FKNU(I)'/3X,'I',6X,'R2') 640 FORMAT(I4,F10.4,F10.4) 650 FORMAT('ETOT=',F10.4,3X,'RTWO=',F6.1) C RETURN END C C********************************************************************* C SUBROUTINE QUAD(XM,XA,XB,DM,DA,DB) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C C SEARCH FOR ROOTS BY FITTING TO QUADRATIC C C C INITIALISATIONS X1=XM X2=XA X3=XB D1=DM D2=DA D3=DB IFL1=1 IFL2=1 C IFL3=1 IF(D1.LT.0)IFL1=-1 IF(D2.LT.0)IFL2=-1 C IF(D3.LT.0)IFL3=-1 C C C START LOOP DO 10 I=1,5 C C FIT TO QUADRATIC AND FIND EXTREMUM U1=X1-X2 U3=X3-X2 Y1=D1-D2 Y3=D3-D2 U=.5*(Y3*U1*U1-Y1*U3*U3)/(Y3*U1-Y1*U3) C C SELECT END POINTS OF RANGE IF(U.LT.0.)THEN X3=X2 D3=D2 C IFL3=IFL2 ELSE X1=X2 D1=D2 IFL1=IFL2 ENDIF C C NEW X2 AND D2 X2=X2+U D2=DETERM(X2) IFL2=1 IF(D2.LT.0.)IFL2=-1 C C LOOK FOR ROOT IF(IFL1*IFL2.LE.0.)THEN CALL BSCT(X1,X2,D1,D2,1) CALL BSCT(X2,X3,D2,D3,1) WRITE(6,605) RETURN ENDIF C 10 CONTINUE C C LOOP COMPLETED C C NO ROOTS FOUND WRITE(6,610)XM,XA,XB,DM,DA,DB C RETURN C 605 FORMAT(/5X,'*** TWO CALLS OF BISECT FROM QUAD'/) 610 FORMAT(5X,23('*'),' NO ROOTS FOUND BY QUAD'/15X,' XM = ', + F10.4,' XA = ',F10.4,' XB = ', F10.4/ + 15X, ' DM = ',E12.4,' DA = ',E12.4,' DB = ',E12.4/) C END C C********************************************************************* C SUBROUTINE READ1 C C READS DATA INDEPENDENT OF SLPI FROM R-MATRIX FILE, FILA C C THE FOLLOWING DATA ARE READ _ C NZ = NUCLEAR CHARGE C NELC = NUMBER OF ELECTRONS IN TARGET C NAST = NUMBER OF TARGET STATES C LRANG2 = TOTAL NUMBER OF SMALL L/BIG K VALUES C LAMAX = MAXIMUM LAMBDA FOR MULTIPOLE POTENTIALS C RA = R-MATRIX RADIUS C BSTO = LOGARITHMIC DERIVATIVE C FOR I = 1,NAST - C ENAT(I) = TARGET ENERGIES C LAT(I) = TARGET ORBITAL ANGULAR MOMENTA C ISAT(I) = VALUES OF (2*S+1) FOR TARGET STATES C FOR I = 1,3 AND L = 1,LRANG2 - C COEFF(I,L) = BUTTLE CORRECTION FITS C FOR I=1,NBUTD(L) FOR L=1,LRANG2 DARC BUTTLE CORRECTION C EBUTD(I,L),CBUTD(I,L) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C LOGICAL EX C COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZ,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/NRBWRD/IWORD,KFLAG COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C C READ AND DIMENSION CHECKS C IF(IWORD.EQ.2)GO TO 50 READ(10)NELC,NZ,LRANG2,LAMAX,NAST,RA,BSTO C IF(LRANG2.LT.0)THEN !FROM DARC DSTG2 VIA DTO3. LRANG2=-LRANG2 KFLAG=-1 ENDIF C C IF(NAST.GT.MZTAR)THEN WRITE(6,610)NAST STOP '*** DIMENSION EXCEEDED: INCREASE MZTAR' ENDIF C READ(10)(ENAT(I),I=1,NAST) READ(10)(LAT(I),I=1,NAST) READ(10)(ISAT(I),I=1,NAST) C IF(LRANG2.GT.MZLP1)THEN WRITE(6,620)LRANG2,MZLP1 STOP '*** DIMENSION EXCEEDED: INCREASE MZLP1' ENDIF C READ(10)((COEFF(I,L),I=1,3),L=1,LRANG2) C INQUIRE(FILE='DBUT.DAT',EXIST=EX) IF(EX)THEN IWORD=-1 IB=14 OPEN(14,FILE='DBUT.DAT',FORM='UNFORMATTED') GO TO 51 ENDIF C RETURN C C DARC 50 READ(10) READ(10)NELC,NZ,NRANG2,LRANG2,NAST,RA,BSTO C IF(NAST.GT.MZTAR)THEN WRITE(6,610)NAST STOP '*** DIMENSION EXCEEDED: INCREASE MZTAR' ENDIF C READ(10)(ENAT(I),I=1,NAST) READ(10)(LAT(I),I=1,NAST) C DO I=1,NAST LAT(I)=ABS(LAT(I))-1 ISAT(I)=0 ENDDO C IF(LRANG2.GT.MZLP1)THEN WRITE(6,620)LRANG2,MZLP1 STOP '*** DIMENSION EXCEEDED: INCREASE MZLP1' ENDIF C IB=10 C 51 DO L=1,LRANG2 READ(IB)NBUTD(L) IF(NBUTD(L).GT.MZNRG)THEN WRITE (6,699) NBUTD(L),MZNRG STOP '*** DIMENSION EXCEEDED: INCREASE MZNRG' ENDIF READ(IB)(EBUTD(I,L),I=1,NBUTD(L)) READ(IB)(CBUTD(I,L),I=1,NBUTD(L)) ENDDO C IF(IB.EQ.14)CLOSE(14) C RETURN C 610 FORMAT(///20X,'TOO MANY TARGET STATES'// 1 10X,'VALUE READ FOR NAST IS ',I3// 2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS MZTAR'//) 620 FORMAT(///20X,'TOO MANNY BUTTLE COEFFICIENTS'// 1 10X,'VALUE READ FOR LRANG2 IS ',I3// 2 10X,'MAXIMUM VALUE ALLOWED BY DIMENSIONS IS MZLP1 =',I3//) 699 FORMAT(//' ******* NBUT = ',I4,' LARGER THAN ', + 'MZNRG = ',I4//) C END C C********************************************************************* C SUBROUTINE READ2 C C READS R-MATRIX DATA FOR ONE SLPI CASE, FROM FILA C C THE FOLLOWING DATA ARE READ - C LRGL2 = TOTAL ORBITAL ANGULAR MOMENTUM C NSPN2 = TOTAL (2*S+1) C NPTY2 = TOTAL PARITY C NCHAN = NUMBER OF CHANNELS C MNP2 = NUMBER OF R-MATRIX POLES C MORE2 = ZERO TO TERMINATE SLPI CASES C FOR I = 1,NAST - C NCONAT(I) = NUMBER OF CHANNELS FOR TARGET STATE I C FOR I = 1,NCHAN - C L2P(I), KJ(I) = SMALL L AND BIG jK FOR CHANNEL I C = K-1, 0 FROM DSTG2/DARC C FOR I = 1,NCHAN AND N = 1,NCHAN AND M = 1,LAMAX - C CF(I,N,M) = COEFFICIENTS IN MULTIPOLE POTENTIALS C FOR I = 1,MNP2 - C VALUE(I) = R-MATRIX POLE ENERGIES C FOR K = 1,NCHAN AND I = 1,MNP2 - C WMAT(K,I) = R-MATRIX AMPLITUDES C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CDEGEN/NASTD,NASTR,NLEV(MZTAR),NCNATR(MZTAR),ENATR(MZTAR) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZ,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/NRBWRD/IWORD,KFLAG COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C IF(IWORD.NE.2)THEN READ(10,END=999)LRGL2,NSPN2,NPTY2,NCHAN,MNP2,MORE2 ELSE !DARC READ(10,END=999)LRGL2,NPTY2,NCHAN,MNP2,NCFGP,LAMAX MORE2=1 LRGL2=ABS(LRGL2)-1 NSPN2=0 IF(NPTY2.EQ.-1)THEN NPTY2=1 ELSE NPTY2=0 ENDIF KFLAG=-2 ENDIF C C DIMENSION CHECKS C IF(NCHAN.GT.MZCHF)THEN WRITE(6,600)NSPN2,LRGL2,NPTY2,NCHAN STOP '*** DIMENSION EXCEEDED: INCREASE MZCHF' ENDIF IF(MNP2.GT.MZMNP)THEN WRITE(6,610)NSPN2,LRGL2,NPTY2,MNP2 STOP '*** DIMENSION EXCEEDED: INCREASE MZMNP' ENDIF IF(LAMAX.GT.MZLMX)THEN WRITE(6,620) STOP '*** DIMENSION EXCEEDED: INCREASE MZLMX' ENDIF C READ(10)(NCONAT(I),I=1,NASTR) ! NAST? READ(10)(L2P(I),I=1,NCHAN) IF(LAMAX.GT.0) XREAD(10)(((CF(I,N,M),I=1,NCHAN),N=1,NCHAN),M=1,LAMAX) C C C IF(MNP2.GE.0)THEN READ(10)(VALUE(I),I=1,MNP2) READ(10)((WMAT(K,I),K=1,NCHAN),I=1,MNP2) ELSE MNP2=-MNP2 READ(10)(VALUE(I),I=1,MNP2) KK=0 13 READ(10)ILOW,IUPPER,NDIV READ(10)((WMAT(I-ILOW+1+KK,K),K=1,NCHAN),I=ILOW,IUPPER) KK=KK+IUPPER-ILOW+1 IF(ILOW.NE.1) GO TO 13 ENDIF C IF(IWORD.EQ.2)THEN !DARC DO I=1,NCHAN L2P(I)=2*L2P(I) !2*kappa IF(L2P(I).LT.0)L2P(I)=-L2P(I)-1 !K, analogue of L L2P(I)=L2P(I)-1 !K-1 ENDDO ENDIF IF(KFLAG.LT.0)THEN !L2P(I)=K-1 DO I=1,NCHAN LLCH(I)=(L2P(I)+1)/2 !CHANNEL ORB ANG K2P(I)=LLCH(I) IF(LLCH(I)*2.NE.L2P(I)+1)K2P(I)=-K2P(I)-1 !KAPPA ENDDO ELSE DO I=1,NCHAN LLCH(I)=L2P(I) !CHANNEL ORB ANG ENDDO ENDIF C DO I=1,NASTR NCNATR(I)=NCONAT(I) ENDDO C C GROUP TOGETHER CHANNELS BELONGING TO DEGENERATE LEVELS C AND PRESERVE THE READ VALUES OF NCONAT IN /CDEGEN/ IF(NASTD.GT.0)THEN N1=1 DO I=1,NASTD NCON=0 N2=NLEV(I)+N1-1 DO IN=N1,N2 NCON=NCON+NCNATR(IN) ENDDO NCONAT(I)=NCON N1=N2+1 ENDDO END IF RETURN C 999 MORE2=-777 !ALLOW EXECUTION TO CONTINUE RETURN C 600 FORMAT(///20X,'TOO MANY CHANNELS FOR (IS, IL, IP) = (', 1 3I3,')'//10X,'VALUE READ FOR NCHAN IS',I3// 2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS CHF = MZCHF'//) 610 FORMAT(///20X,'TOO MANY R-MATRIX STATES FOR (IS, IL, IP) = (', * 3I3,')'//10X,'VALUE READ FOR MNP2 IS',I5// 3 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS MNP = MZMNP'//) 620 FORMAT(///20X,'TOO MANY MULTIPOLES'// 1 10X,'VALUE READ FOR LAMAX IS',I3// 2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IF LMX = MZLMX'//) C END C C********************************************************************* C SUBROUTINE RMATB1(KP) C C THE R-MATRIX IS C R = T + S/DE C WHERE C DE = VALUE(KP) - ETOT C THIS SUBROUTINE CALCULATES T IN RMAT. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C LOGICAL NEWBUT C COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CRMAT/RMAT(MZCHF,MZCHF),RMATD(MZCHF,MZCHF) 1,UMAT(MZCHF,MZCHF),GAMMA COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CLOGB/NEWBUT COMMON/NRBWRD/IWORD,KFLAG COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C C DO I=1,NCHF DO J=I+1,NCHF RMAT(J,I)=0. ENDDO ENDDO C C INITIALISE DIAGONAL ELEMENTS TO BUTTLE CORRECTION C IF(IWORD.EQ.1)THEN IF(NEWBUT)THEN RA2=RZERO*RZERO DO I=1,NCHF L=L2P(I)+1 NBUT=COEFF(3,L) RMAT(I,I)=COEFF(1,L)*BUT0(NBUT,COEFF(2,L)+RA2*EPS(I)) ENDDO ELSE DO I=1,NCHF L=L2P(I)+1 E=EPS(I) RMAT(I,I)=COEFF(1,L)+E*(COEFF(2,L)+E*COEFF(3,L)) ENDDO ENDIF ELSE !DARC DO I=1,NCHF L=L2P(I)+1 E=EPS(I) CALL INTBUT(L,E,BUTT) RMAT(I,I)=BUTT ENDDO ENDIF C DO K=1,KP-1 V=1./(VALUE(K)-ETOT) DO I=1,NCHF VI=V*WMAT(I,K) DO J=I,NCHF RMAT(J,I)=RMAT(J,I)+VI*WMAT(J,K) ENDDO ENDDO ENDDO DO K=KP+1,MNP2 V=1./(VALUE(K)-ETOT) DO I=1,NCHF VI=V*WMAT(I,K) DO J=I,NCHF RMAT(J,I)=RMAT(J,I)+VI*WMAT(J,K) ENDDO ENDDO ENDDO C C SYMMETRIZE C DO I=2,NCHF DO J=1,I-1 RMAT(J,I)=RMAT(I,J) ENDDO ENDDO C RETURN END C C********************************************************************* C SUBROUTINE RMATB2(KP) C C THE R-MATRIX IS C R = T + S/DE C WHERE C DE = VALUE(KP) - ETOT C THIS SUBROUTINE CALCULATES T IN RMAT AND ENERGY DERIVATIVE C OF T IN RMATD. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C LOGICAL NEWBUT C COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CRMAT/RMAT(MZCHF,MZCHF),RMATD(MZCHF,MZCHF) 1,UMAT(MZCHF,MZCHF),GAMMA COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CLOGB/NEWBUT COMMON/NRBWRD/IWORD,KFLAG COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C DE=VALUE(KP)-ETOT C DO I=1,NCHF DO J=I+1,NCHF RMAT(J,I)=0. RMATD(J,I)=0. ENDDO ENDDO C C INITIALISE DIAGONAL ELEMENTS TO BUTTLE CORRECTION C IF(IWORD.EQ.1)THEN IF(NEWBUT)THEN RA2=RZERO*RZERO DO I=1,NCHF L=L2P(I)+1 NBUT=COEFF(3,L) U=COEFF(2,L)+RA2*EPS(I) CALL CBUT0(NBUT,U,B,C) RMAT(I,I)=B*COEFF(1,L) RMATD(I,I)=C*COEFF(1,L)*RA2 ENDDO ELSE DO I=1,NCHF L=L2P(I)+1 E=EPS(I) RMAT(I,I)=COEFF(1,L)+E*(COEFF(2,L)+E*COEFF(3,L)) RMATD(I,I)=COEFF(2,L)+2.*E*COEFF(3,L) ENDDO ENDIF ELSE !DARC DO I=1,NCHF L=L2P(I)+1 E=EPS(I) CALL INTBUT(L,E,BUTT) CALL INTBUTD(L,E,BUTTD) RMAT(I,I)=BUTT RMATD(I,I)=BUTTD ENDDO ENDIF C DO K=1,KP-1 V=1./(VALUE(K)-ETOT) DO I=1,NCHF VI=V*WMAT(I,K) DO J=I,NCHF VIJ=VI*WMAT(J,K) RMAT(J,I)=RMAT(J,I)+VIJ RMATD(J,I)=RMATD(J,I)+VIJ*V ENDDO ENDDO ENDDO DO K=KP+1,MNP2 V=1./(VALUE(K)-ETOT) DO I=1,NCHF VI=V*WMAT(I,K) DO J=I,NCHF VIJ=VI*WMAT(J,K) RMAT(J,I)=RMAT(J,I)+VIJ RMATD(J,I)=RMATD(J,I)+VIJ*V ENDDO ENDDO ENDDO C DO I=2,NCHF DO J=1,I-1 RMAT(J,I)=RMAT(I,J) RMATD(J,I)=RMATD(I,J) ENDDO ENDDO C RETURN END C C********************************************************************* C SUBROUTINE SCALE1 C C CONVERTS R-MATRIX TARGET DATA TO Z-SCALED FORM C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C LOGICAL NEWBUT C C COMMON BLOCKS FROM ASYMPTOTIC ROUTINE COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CPOT/LAMP(MZCHF,MZCHF),BW(MZCHF,MZCHF) COMMON/CENAT1/ENAT1 COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/CLOGB/NEWBUT COMMON/CBUT/FKN(0:MZNRG),UKN(0:MZNRG) C C COMMON BLOCK FROM SUBROUTINE READ C NOTE USE OF NZED IN PLACE OF NZ COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/NRBWRD/IWORD,KFLAG COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C COMMON BLOCK FOR DEGENERATE TARGET LEVELS COMMON/CDEGEN/NASTD,NASTR,NLEV(MZTAR),NCNATR(MZTAR),ENATR(MZTAR) C DIMENSION INDEX(MZTAR),AVENGY(MZTAR) c IWT=0 C C Z-SCALING FACTORS C AZ=MAX(NZED-NELC,1) AZ1=ONE/AZ AZ2=AZ1*AZ1 TAZ2=2*AZ2 AZAZ=AZ*AZ C C C Z-SCALED TARGET ENERGIES, RELATIVE TO ZERO FOR TARGET GROUND STATE C ENAT1=ENAT(1) COLD IF(NAST.EQ.1)GOTO 20 DO I=1,NAST ENAT(I)=TAZ2*(ENAT(I)-ENAT1) ENDDO C C TEST FOR DEGENERACY AND FORM AVERAGED TARGET ENERGY LEVELS C IF (NASTD.GT.0) THEN NASTR=0 DO N=1,NASTD NASTR=NASTR+NLEV(N) ENDDO IF(NASTR.NE.NAST)THEN WRITE(6,640)NASTR,NAST,(NLEV(N),N=1,NASTD) STOP 'LEVEL INCONSISTANCY ON GROUPING...' END IF AVENGY(1)=0.0 N1=NLEV(1)+1 INDEX(1)=NLEV(1) DO 14 J=2,NASTD N=NLEV(J) N2=N+N1-1 INDEX(J)=N2 ENSUM=0.0 DO 15 IN=N1,N2 ENSUM=ENAT(IN)+ENSUM 15 CONTINUE AVENGY(J)=ENSUM/N N1=N2+1 14 CONTINUE 20 ENAT(1)=0. NASTR=NAST C COLD AVENGY(1)=0.0 COLD N1=NLEV(1)+1 COLD INDEX(1)=NLEV(1) N1=1 DO J=1,NASTD N=NLEV(J) N2=N+N1-1 INDEX(J)=N2 ENSUM=TZERO WTSUM=TZERO DO IN=N1,N2 IF(IWT.LT.0)THEN WT=LAT(IN)+1 ELSEIF(IWT.EQ.0)THEN WT=1 ELSE WT=ISAT(IN)*(2*LAT(IN)+1) ENDIF ENSUM=ENAT(IN)*WT+ENSUM WTSUM=WTSUM+WT ENDDO AVENGY(J)=ENSUM/WTSUM N1=N2+1 ENDDO ENDIF C COL20 ENAT(1)=0. NASTR=NAST C C RA AND BSTO RZERO=RA*AZ BSTO=BSTO/RZERO C C BUTTLE CORRECTION C IF(IWORD.EQ.1)THEN IF(COEFF(3,1).GT.-10000)THEN NEWBUT=.FALSE. AA=RZERO*AZ2 DO M=1,3 AA=AA*AZAZ DO L=1,LRANG2 COEFF(M,L)=COEFF(M,L)*AA ENDDO ENDDO ELSE NEWBUT=.TRUE. DO L=1,LRANG2 NBUT=-INT(COEFF(3,L))/10000 IF(NBUT.GT.MZNRG)THEN WRITE(6,699)NBUT,MZNRG STOP '*** DIMENSION EXCEEDED: INCREASE MZNRG' ENDIF COEFF(3,L)=NBUT COEFF(1,L)=RZERO*COEFF(1,L) ENDDO C C INITIALISE FKN AND UKN PI=ACOS(-ONE) G=-PI/TWO DO I=0,MZNRG G=G+PI FKN(I)=G UKN(I)=G*G ENDDO ENDIF ELSE !DARC DO L=1,LRANG2 DO N=1,NBUTD(L) EBUTD(N,L)=EBUTD(N,L)*TAZ2 CBUTD(N,L)=CBUTD(N,L)*RZERO ENDDO ENDDO ENDIF C C WRITE TARGET PROPERTIES C WRITE(6,650)NZED,NELC WRITE(6,655) DO J=1,NAST WRITE(6,660)J,ISAT(J),LAT(J),ENAT(J) ENDDO C C IF THERE ARE DEGENERATE LEVELS PRINT THE AVERAGED ENERGIES C PUT THEM IN /CINPUT/ AND PRESERVE THE ENERGIES READ FROM C THE H FILE IN /CDEGEN/ C IF(NASTD.NE.0)THEN DO I=1,NAST ENATR(I)=ENAT(I) ENDDO INR1=1 DO IND=1,NASTD N=NLEV(IND) INR2=INR1+N-1 IF(N.GT.1)WRITE(6,664)INR1,INR2 INR1=INR2+1 ENDDO WRITE(6,665)IWT AV1=AVENGY(1) DO J=1,NASTD ENAT(J)=AVENGY(J)-AV1 WRITE(6,666)INDEX(J),J,ENAT(J) ENDDO C NAST=NASTD END IF C WRITE(6,670)RZERO,BSTO WRITE(6,680)AC C 640 FORMAT(///5X,'*****INCORRECT DATA ',/ 1 5X,'EXPECT ',I3, 'LEVELS FROM NLEV DATA BUT NAST IS', 2 I3,/5X,'NLEV=',20I3) 650 FORMAT(' ',10X,'NUCLEAR CHARGE =',I3,', NUMBER OF TARGET ', 1 ' ELECTRONS =',I3/11X,53('*')//) 655 FORMAT(20X,'TARGET STATES -'/20X,15('*')// 1 10X,'INDEX',5X,'2*S+1',5X,'TOTAL L',5X,'SCALED ENERGY'/ 2 20X,'OR P OR 2*J'/) 660 FORMAT(3X,3I10,7X,F12.6) 664 FORMAT(/' LEVELS ',I3,' TO ',I3,' ARE COMBINED') 665 FORMAT(//' IWT=',I2,8X,'EQUIVALENT TARGET STATES -'/15X,26('*')// 1 12X,'OLD INDEX',5X,'NEW INDEX',7X,'SCALED ENERGY'//) 666 FORMAT(3X,2I14,9X,F12.6) 670 FORMAT(//5X,'RZERO = ',F10.4,3X,'BSTO = ',F10.4/) 680 FORMAT(5X,'AC = ',E12.4) 699 FORMAT(//' ******* NBUT = ',I4,' LARGER THAN ', + 'MZNRG = ',I4//) RETURN END C C********************************************************************* C SUBROUTINE SCALE2 C C CONVERTS R-MATRIX DATA TO Z-SCALED FORM C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C C COMMON BLOCKS FROM ASYMPTOTIC ROUTINE COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 COMMON/CPOT/LAMP(MZCHF,MZCHF),BW(MZCHF,MZCHF) C C COMMON BLOCK FROM SUBROUTINE READ C NOTE USE OF NZED IN PLACE OF NZ COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CENAT1/ENAT1 COMMON/CNTRL/IPRINT,IRAD,IPERT COMMON/NRBWRD/IWORD,KFLAG COMMON/DBUT/EBUTD(MZNRG,MZLP1),CBUTD(MZNRG,MZLP1),NBUTD(MZNRG) X ,K2P(MZCHF) C C STATISTICAL WEIGHT * 2 (OR * 4 FOR NON-EXCHANGE CASE) C LRGL2=2*J FOR B.P. C C C NWT=(2*LRGL2+1) IF(NSPN2.EQ.0)THEN NWT=NWT+1 ELSE IF(NSPN2.LT.0)NWT=-2*NWT NWT=NWT*NSPN2*2 ENDIF C C Z-SCALING FACTORS C AZ=MAX(NZED-NELC,1) AZ1=ONE/AZ AZ2=AZ1*AZ1 TAZ2=2*AZ2 AZHR=ONE/SQRT(AZ) AZAZ=AZ*AZ C C CHANNELS C NCHF=NCHAN DO I=1,NCHF LL=LLCH(I) CC(I)=DBLE(LL*(LL+1)) ENDDO I=0 DO 50 J=1,NAST K=NCONAT(J) IF(K.EQ.0)GOTO 50 DO L=1,K I=I+1 ITARG(I)=J ECH(I)=ENAT(J) ENDDO 50 CONTINUE C C R-MATRIX C VALUE AND WMAT C DO N=1,MNP2 DO I=1,NCHF WMAT(I,N)=WMAT(I,N)*AZHR ENDDO ENDDO IF(IPRINT.GT.2)WRITE(6,635) DO N=1,MNP2 VALUE(N)=TAZ2*(VALUE(N)-ENAT1) IF(IPRINT.GT.2)WRITE(6,640)N,VALUE(N),(WMAT(I,N),I=1,NCHF) ENDDO C C COEFFICIENTS IN POTENTIAL C LMX=LAMAX IF(IPERT.NE.1)GOTO 160 IF(IPRINT.GT.2)WRITE(6,655) C C LAMP(I,J) AND BW(I,J) C DO J=1,NCHF DO I=1,NCHF LAMP(I,J)=1 ENDDO ENDDO IF(LMX.EQ.0)GOTO 160 C C DIPOLE C IABW=IABS(IWORD) A1=ONE/RZERO A2=A1/2 DO J=2,NCHF J1=J-1 DO 130 I=1,J1 CF(I,J,1)=CF(I,J,1)*IABW BU1=-CF(I,J,1) IF(ABS(BU1).LT.1.D-6)GOTO 130 LAMP(I,J)=2 BW(I,J)=BU1 P=A2*BU1 IF(IPRINT.GT.2)WRITE(6,660)I,J,P 130 CONTINUE ENDDO IF(LMX.EQ.1)GOTO 160 C C QUADRUPOLE C IF(IPRINT.GT.2)WRITE(6,665) A2=A2*A1 DO J=1,NCHF DO 150 I=1,J CF(I,J,2)=CF(I,J,2)*IABW BU2=-CF(I,J,2) IF(ABS(BU2).LT.1.D-6)GOTO 150 LAMP(I,J)=3 BU2=BU2*AZ BW(I,J)=BU2 CF(I,J,2)=AZ*CF(I,J,2) P=BU2*A2 IF(IPRINT.GT.2)WRITE(6,660)I,J,P 150 CONTINUE ENDDO C 160 CONTINUE C C SYMMETRIZE C DO I=1,NCHF DO J=I,NCHF LAMP(J,I)=LAMP(I,J) BW(J,I)=BW(I,J) ENDDO ENDDO C C WRITE SLPI AND CHANNEL DATA C WRITE(6,600)NSPN2,LRGL2,NPTY2 IF(KFLAG.GE.0)THEN WRITE(6,669) WRITE(6,670)(I,ITARG(I),LLCH(I),I=1,NCHF) ELSE WRITE(6,673) WRITE(6,672)(I,ITARG(I),LLCH(I),K2P(I),I=1,NCHF) ENDIF C 600 FORMAT(///80('+')//12X,'(2S+1) L/2J P =',3I3/12X,24('*')//) 635 FORMAT(//' N, VALUE(N) AND WMAT(I,N)'/) 640 FORMAT(I5,E12.3,5X,9E12.3,(/5X,9E12.3)) 655 FORMAT(/' PERTURBATION P FOR MULTIPOLE POTENTIAL'// 1 ' - DIPOLE PART'/) 660 FORMAT(2I5,E12.4) 665 FORMAT(/' - QUADRUPOLE PART'/) 669 FORMAT(12X,'CHANNEL',2X,'TARGET',4X,'SMALL', 1 /12X,'INDEX ',3X,'INDEX ',5X,'L'/) 670 FORMAT(7X,I8,I9,I10) 672 FORMAT(7X,I8,I9,I10,I7) 673 FORMAT(12X,'CHANNEL',2X,'TARGET',4X,'SMALL', 1 /12X,'INDEX ',3X,'INDEX ',5X,'L',4X,'KAPPA'/) C RETURN END C C********************************************************************* C SUBROUTINE SCAN(X0,X1,BDX,XPM) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (NDUST=1000) C COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CNTRL/IPRINT,IRAD,IPERT COMMON/CITER/ECH1,XSTORE(MZEST),K0,KSTORE(MZEST) DIMENSION FNUST(NDUST),DETST(NDUST) C C INITIALISE FOR NUMBER OF POINTS AND INTERVAL C N=1+(X1-X0)/BDX N=(X1-X0)/BDX IF(IPRINT.EQ.1.AND.N.GT.NDUST)THEN WRITE(6,*)'TOO MANY SCAN POINTS, IPRINT HAS BEEN RESET=0' IPRINT=0 ENDIF IF(N.LT.5)N=5 SDX=N SDX=(X1-X0)/SDX C C INITIALISE FOR SCAN CALL UTRANS(K0) XM=MAX(X0-SDX,0.5*(X0+XPM)) DM=DETERM(XM) IFLM=1 IF(DM.LT.0.)IFLM=-1 XA=X0 DA=DETERM(XA) IFLA=1 IF(DA.LT.0.)IFLA=-1 IF(IFLA*IFLM.LT.0.)THEN KROOT=1 ELSE KROOT=2 SLA=(DA-DM)/(DA+DM) ENDIF C C SCAN FROM X=X0 TO X1 XB=X0 DO I=1,N IF(IPRINT.EQ.1)THEN FNUST(I)=XA DETST(I)=DA ENDIF XB=XB+SDX DB=DETERM(XB) IFLB=1 IF(DB.LT.0.)IFLB=-1 C LOOK FOR ROOT IF(IFLA*IFLB.LE.0.)THEN CALL BSCT(XA,XB,DA,DB,0) KROOT=1 ELSE KROOT=KROOT+1 SLB=(DB-DA)/(DB+DA) IF(KROOT.GT.2)THEN IF(SLA.LE.0..AND.SLB.GT.0.)CALL QUAD(XM,XA,XB,DM,DA,DB) ENDIF ENDIF IF(MXE.GT.MZEST)go to 1000 XM=XA DM=DA IFLM=IFLA XA=XB DA=DB IFLA=IFLB IF(KROOT.GT.1)SLA=SLB ENDDO C IF(IPRINT.EQ.1)WRITE(6,600)(FNUST(I),DETST(I),I=1,N),XA,DA C 1000 RETURN C 600 FORMAT(5X,'VALUES OF EFFECTIVE QUANTUM NUMBER AND DETERMINANT'/ 1 /5(0PF10.4,1PE9.1)) C END C C********************************************************************* C SUBROUTINE SOLV(A,B,N,D,F,X,ICONV) C C A AND B ARE SQUARE MATRICES WITH N ROWS AND N COLUMNS. C X IS A COLUMN VECTOR WITH N COLUMNS. C THE EIGENVALUE PROBLEM C (A - F*B)*X = 0 C IS SOLVED FOR THE SMALLEST EIGENVALUE F, C TO AN ACCURACY OF D IN F. C C THE ORIGINAL MATRIX A IS DESTROYED. C THE MATRIX B IS RETURNED UNCHANGED. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C DIMENSION A(MZCHF,MZCHF),B(MZCHF,MZCHF),X(MZCHF) DIMENSION IR(MZCHF),JC(MZCHF) DIMENSION C(MZCHF,MZCHF),Y(MZCHF) C C CASE OF N.EQ.1 C IF(N.EQ.1)THEN F=A(1,1)/B(1,1) X(1)=1. RETURN ENDIF C C CALL ROUTINE FOR LOWER-UPPER DECOMPOSITION OF A, C WITH FULL PIVOTING AND INTERCHANGES. C CALL LU(A,N,IR,JC) C C CALCULATE C = (A**-1)*B, WITH INTERCHANGES C - INTERCHANGES C DO J=1,N DO I=1,N C(I,J)=B(IR(I),JC(J)) ENDDO ENDDO C C - DIVIDE BY A C DO J=1,N DO I=2,N DO K=1,I-1 C(I,J)=C(I,J)-A(I,K)*C(K,J) ENDDO ENDDO ENDDO DO I=1,N Y(I)=1./A(I,I) ENDDO DO J=1,N DO I=N,1,-1 DO K=I+1,N C(I,J)=C(I,J)-A(I,K)*C(K,J) ENDDO C(I,J)=C(I,J)*Y(I) ENDDO ENDDO C C INITIAL ESTIMATE OF EIGENVALUE AND EIGENVECTOR C F0=0. X(N)=1. DO I=N-1,1,-1 S=0. DO J=I+1,N S=S-A(I,J)*X(J) ENDDO X(I)=S/A(I,I) ENDDO C C START ITERATION LOOP C DO M=1,10 C CALCULATE Y = C*X C DO I=1,N Y(I)=0. ENDDO DO J=1,N DO I=1,N Y(I)=Y(I)+C(I,J)*X(J) ENDDO ENDDO C C NEW ESTIMATES OF EIGENVALUE AND EIGENVECTOR C C - EIGENVALUE C S=0. T=0. DO I=1,N YI=Y(I) S=S+X(I)*YI T=T+YI*YI ENDDO F=S/T C C - EIGENVECTOR C DO I=1,N X(I)=F*Y(I) ENDDO C C TEST CONVERGENCE C IF(ABS(F-F0).LT.D)THEN ICONV=0 GOTO 130 ENDIF F0=F ENDDO C C END OF ITERATION LOOP C C NOT CONVERGED AFTER 10 ITERATIONS C ICONV=1 C C INTERCHANGES FOR X. C 130 DO J=1,N Y(JC(J))=X(J) ENDDO DO J=1,N X(J)=Y(J) ENDDO C RETURN END C C********************************************************************* C SUBROUTINE TANDTD(R,I,TA,TDA,TP) C C CALCULATES THETA AND THETAD =THETA DOT FOR R REAL C THETA = TA*EXP(TP) C THETAD = TDA*EXP(TP) C TP = FNU*LOG(R) - R/FNU C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CTHET/BB(MZCHF,MZTET),BG(MZCHF,MZTET),MSUM(MZCHF) COMMON/NRBZED/TZED C MI=MSUM(I) S=BB(I,2) CX=0. E=BG(I,1) F2=BG(I,2) FNU=BB(I,1) C X=2.*R/FNU Y=1./X AS=1. DO 10 L=3,MI AS=AS*Y S=S+BB(I,L)*AS 10 CX=CX+BG(I,L)*AS C DLR=LOG(R)*TZED TP=-.5*X+FNU*DLR TA=S TDA=E*((DLR+R*F2)*S+CX) C RETURN END C C********************************************************************* C SUBROUTINE THETA(R,I,T,TP,TD,TDP) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/CTHET/BB(MZCHF,MZTET),BG(MZCHF,MZTET),MSUM(MZCHF) COMMON/NRBZED/TZED LOGICAL FLAG C C FNU=FKNU(I) LL=LLCH(I) M=FNU+LL+12 IF(TZED.EQ.0)M=MZTET C IF(M.GT.MZTET)THEN WRITE(6,600)ETOT,I,FNU,M STOP ENDIF C F1=1./FNU F2=F1*F1 X=2.*R*F1 Y=1./X FL=LL C R1=1./R A=TZED*FNU*R1-F1 B=TZED*LOG(R)+R*F2 C=TZED*R1+F2 D=-R1 E=.5*FNU**3 C BB(I,1)=FNU BB(I,2)=1. BG(I,1)=E BG(I,2)=F2 C YN=1. A1=FL-FNU*TZED A2=FL+FNU*TZED+1. BN=-2.*FNU-1. C BET=1. GAM=0. W0=C C S=1. U=0. CX=0. CY=0. C C N=0 FLAG=.FALSE. DO 20 MN=3,M N=N+1 A1=A1+1. A2=A2-1. BN=BN+2. CN=A1*A2 DN=BN*TZED+F1*CN GAM=CN*GAM+DN*BET BET=CN*BET YN=YN*Y U=U+BET*YN CY=CY+GAM*YN AN=N AN=1./AN GAM=GAM*AN BET=BET*AN S=S+BET*YN CX=CX+GAM*YN W1=C*S*S+D*(S*CY-U*CX) BB(I,MN)=BET BG(I,MN)=GAM IF(W1.EQ.0.0)GO TO 30 IF(ABS((W1-W0)/W1).LT.AC)THEN IF(FLAG)THEN GO TO 30 ELSE FLAG=.TRUE. END IF ELSE FLAG=.FALSE. END IF W0=W1 20 CONTINUE C IF(R*F1.GT.150.)THEN T=0.0 TP=0.0 TD=0.0 TDP=0.0 RETURN ENDIF C C NOT CONVERGED WRITE(6,610)ETOT,I,FNU,N STOP C C SUMMATIONS CONVERGED C 30 P=EXP(-.5*R*F1) IF(TZED.GT.0)P=P*R**(.5*FNU) C-NRB IF(ABS(P).GT.1.D-100)THEN CFACT=1./P ELSE CFACT=0.0 ENDIF C-NRB BB(I,2)=BB(I,2)*CFACT DO 40 J=3,N+2 BB(I,J)=BB(I,J)*CFACT 40 BG(I,J)=BG(I,J)*CFACT C T=P*S TP=P*(A*S+D*U) TD=P*E*(B*S+CX) TDP=P*E*((A*B+C)*S+B*D*U+A*CX+D*CY) N2=N+2 MSUM(I)=N2 C 600 FORMAT(//5X,55('*')/5X,'SUBROUTINE THETA'/5X,'ETOT = ',E12.5, + ', FNU(',I2,') = ',F10.4/5X,'M = ',I4,' WHICH IS LARGER THAN', + ' MAXIMUM VALUE OF TET=MZTET'/5X,'ALLOWED BY DIMENSIONS. ' + ,'VALUE OF FNU TOO LARGE.'/5X,55('*')//) 610 FORMAT(//5X,55('*')/5X,'SUBROUTINE THETA'/5X,'ETOT = ',E12.5, + ', FNU(',I2,') = ',F10.4/5X,'SUMMATIONS NOT CONVERGED WITH', + I4,' TERMS. ',' VALUE OF FNU TOO LARGE.'/5X,55('*')//) RETURN END C C********************************************************************* C SUBROUTINE TLAG(I,J,NLAG,LIJ,T1,T2,T3) C C CALCULATES T INTEGRALS FOR CHANNELS I,J USING C LAGUERRE QUADRATURE WITH NLAG POINTS C THE INTEGRALS ARE - C T1 FOR (THETAI,THETAJ) C T2 FOR (THETADI,THETAJ) C T3 FOR (THETAI,THETADJ) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CTHET/BB(MZCHF,MZTET),BG(MZCHF,MZTET),MSUM(MZCHF) COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2 C C CAP K FNUI=BB(I,1) FNUJ=BB(J,1) FK=1./FNUI+1./FNUJ C INITIALISE FOR QUADRATURE NS=NLAG/2 M=NS*(NS-1) N1=M+1 N2=M+NLAG T1=0. T2=0. T3=0. C START QUADRATURE DO 40 N=N1,N2 U=XLAG(N) R=RTWO+U/FK C CALCULATE THETA FUNCTIONS CALL TANDTD(R,I,TI,TDI,TPI) CALL TANDTD(R,J,TJ,TDJ,TPJ) C ADD TO SUM C++ VAX MOD C A1=(R**(-LIJ))*WLAG(N)*EXP(U+TPI+TPJ) A1=(R**(-LIJ))*WLAG(N) U2=.5*U AI=EXP(U2+TPI) AJ=EXP(U2+TPJ) TI=TI*AI TDI=TDI*AI TJ=TJ*AJ TDJ=TDJ*AJ C++ END MOD T1=T1+TI*A1*TJ T2=T2+TDI*A1*TJ 40 T3=T3+TI*A1*TDJ F1=1./FK T1=T1*F1 T2=T2*F1 T3=T3*F1 C RETURN END C C********************************************************************* C SUBROUTINE UTRANS(KP) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CRMAT/RMAT(MZCHF,MZCHF),RMATD(MZCHF,MZCHF) 1,UMAT(MZCHF,MZCHF),GAMMA C C A1=WMAT(1,KP) S2=A1*A1 DS2=ABS(A1) C DO 60 J=1,NCHAN-1 S1=S2 DS1=DS2 A1=WMAT(J+1,KP) S2=S1+A1*A1 DS2=SQRT(S2) A1=A1/(DS1*DS2) DO 40 I=1,J 40 UMAT(I,J)=WMAT(I,KP)*A1 UMAT(J+1,J)=-DS1/DS2 DO 50 I=J+2,NCHAN 50 UMAT(I,J)=0. 60 CONTINUE C A1=1./DS2 DO 70 I=1,NCHAN 70 UMAT(I,NCHAN)=WMAT(I,KP)*A1 C GAMMA=S2 C RETURN END C C********************************************************************* C SUBROUTINE VECT(KR,X) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (ZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CVECT/A(MZCHF,MZCHF),B(MZCHF,MZCHF),C(MZCHF,MZCHF), 1 P(MZCHF,MZCHF),Q(MZCHF,MZCHF),XVECT(MZCHF) COMMON/CRMAT/RMAT(MZCHF,MZCHF),RMATD(MZCHF,MZCHF) 1,UMAT(MZCHF,MZCHF),GAMMA COMMON/CINPUT/ 3 BSTO,RA, 4 CF(MZCHF,MZCHF,MZLMX),COEFF(3,MZLP1),ENAT(MZTAR), 5 WMAT(MZCHF,MZMNP),VALUE(MZMNP), 1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2, 2 ISAT(MZTAR),LAT(MZTAR),L2P(MZCHF),NCONAT(MZTAR) COMMON/CHAN/ECH(MZCHF),LLCH(MZCHF),EPS(MZCHF),FKNU(MZCHF) 1 ,CC(MZCHF),RINF(MZCHF),ITARG(MZCHF),NCHF,NCHOP,NCHOP1 C ETOT=ECH(1)-1./(X*X) DE=VALUE(KR)-ETOT DER=1./DE DER2=DER*DER DO J=1,NCHF DO I=1,NCHF T=WMAT(I,KR)*WMAT(J,KR) RMAT(I,J)=RMAT(I,J)+T*DER RMATD(I,J)=RMATD(I,J)+T*DER2 ENDDO ENDDO C CSTRTNBL CALL MATMUL(RMATD,Q,A,MZCHF,NCHF) CALL MTMLT(Q,A,B,MZCHF,NCHF) CENDNBL CSTRTBL C CALL DGEMM('N','N',NCHF,NCHF,NCHF,ONE,RMATD,MZCHF,Q,MZCHF, C X ZERO,A,MZCHF) C CALL DGEMM('T','N',NCHF,NCHF,NCHF,ONE,Q,MZCHF,A,MZCHF, C X ZERO,B,MZCHF) CENDBL C A1=0. DO I=1,NCHF A1=A1+XVECT(I)*XVECT(I) ENDDO C DO J=1,NCHF DO I=1,NCHF A1=A1+XVECT(I)*B(I,J)*XVECT(J) ENDDO ENDDO C A1=1./SQRT(A1) DO I=1,NCHF XVECT(I)=XVECT(I)*A1 ENDDO C RETURN END