C N. R. BADNELL UoS v2.16 - QUB v1.2 07/10/20 C C PROGRAM STGBF C C CALCULATE PHOTOIONIZATION CROSS SECTIONS, C AND PHOTORECOMBINATION (IPRINT=-2) COLLISION STRENGTHS (NRB). C MQDT MODE (WITH STGF) PLUS OPTIONAL PRECONVOLUTION (NRB). C BETA PARAMETERS NOW IN BP MODE AS WELL (TWG). C STGBF MUST BE PRECEDED BY PREBF WHICH PROCESSES THE DIPOLE C DATA TO SAVE STORAGE IN STGBF. C FURTHER STORAGE CAN BE SAVED BY USE OF THE PARAMETERS MZSA1,2,3 C SET MZSA1,2,3=1 FOR MINIMUM STORAGE. C SET MZSA1=MZEST FOR MAXIMUM EFFICIENCY. C SET MZSA2=MZMSH FOR MAXIMUM EFFICIENCY (IPRINT>-2 ONLY). C SET MZSA3=MZMSH FOR MAXIMUM EFFICIENCY (IPRINT=3 ONLY). C C I/O CHANNELS: C 1 INITIAL BOUND STATE DATA (FILES B); C 2 FINAL FREE STATE DATA (FILES F); C 3 DIPOLE DATA FOR THE INITIAL AND FINAL SLPI (FILES D); C 5 STANDARD INPUT. THE FOLLOWING DATA ARE READ: C IPRINT, IBUT C IS1, IL1, IP1, M11, M12 = S L PI OR 0 2J PI ... C IS2, IL2, IP2 - SET 0 0 2 FOR J=0 EVEN C IS2, IL2, IP2 - AS PI CONGR PI MODULO 2! C . . . C -1, -1, -1 :TERMINATING 1 INITIAL SLP C IS1, IL1, IP1, M11, M12 :AND STARTING ANOTHER C . . . C . . . C -1, -1, -1, -1, -1 :OR SIMPLY /EOF C 6 STANDARD OUTPUT. OUTPUT LEVEL IS CONTROLED BY IPRINT: C IPRINT=-3 AS -2 BUT RECOMBINATION PRODUCTION. C NO PERT, ONLY OMEGPR CALCULATED (AND C STORED INTERNALLY TO REDUCE I/O). C -2 CROSS SECTION FOR LEAVING THE RESIDUAL C ION IN EACH OF THE TARGET STATES (DR). C -1 NO PRINTING OF CROSS SECTIONS, LEVELS C ABOVE IONIZATION LIMIT SKIPPED. USED C FOR PRODUCTION RUNS. NO PERT. C 0 TOTAL CROSS SECTIONS TO A FINAL SLPI, C 1 FILES OPENED, QDT ONSET, C 2 PARTIAL CROSS SECTION TO EACH CHANNEL; C 3 BETA AND X-SECTN TO EACH TARGET STATE. C 9 BETA FILE. PHOTON ENERGIES, X-SECTN, BETA PARAMETERS. C 7 ARCHIVE XSCTN FILE. PHOTON ENERGIES, TOTAL CROSS SECTIONS. C 18 OMEGPR FILE, STORES RECOMBINATION OMEGAS (IPRINT=-2). C C PROGRAM MAIN C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMENT-OUT ABOVE LINE (THROUGHOUT CODE) AND ADD BACK INTO PARAM C WITH MZMSH0.LT.MZMSH IF THIS PROVES TOO LARGE. MZMSH0 IS KEPT C DISTINCT THROUGHOUT THE CODE - NRB. C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (IFAC=50) PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C SUN REAL*4 TARRY(2) C LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,FLAGV,SKIPPY, 1 QDT,NOQDT,CMPLT,FIRST,EXISTV,SAVE,QJUMP,FLAGL,EJUMP C CHARACTER NAME*6,NUM(0:9)*1 CHARACTER*3 FILE1,FILE2, FILE3*4 C COMMON/CACC/AC0,ACNUM,ACJWBK,ACZP,LACC COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/SCALE/AZ,AZP,AZ2,AZP2,K10,K20 COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/SLPI2/IS2V(3),IL2V(3),IP2V(3),KUT COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT COMMON/FLAGQ/QDT COMMON /TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) COMMON /VUSTOR/ VSTORE(-3:1,MZCHF,MZCHF),USTORE(-3:1,MZCHF,MZCHF), * X1(MZCHF,MZEST),X2(MZCHF),EXISTV(-3:1,MZCHF,MZCHF) COMMON/CSTOR/FAC(2,IFAC),ITOR,IS1,IL1,IP1,M11,MFIN,M1,IL2,M2,FLAGV X ,FLAGL C COMMON/NRBDR/OMEGPR(MZMET,MZMSH),NDRMET,IGAUGE COMMON/NRBFLG/TOL0V,IFLAGV COMMON/NRBPRT/IPERT1,IPERT2 COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBTMP/ETOT COMMON/NRBZED/TZED C COMMON/NRBDDD/DLVECT(MZCHF,MZEST),DVVECT(MZCHF,MZEST) X ,DLT(MZEST),DL2(MZCHF),DV2(MZCHF) X ,TOT(MZSA2,MZEST),DEE(MZEST),DEER(MZEST) X ,DLI(MZCHF,MZEST),DVI(MZCHF,MZEST) C DIMENSION EMESH(MZMSH),EMSH0(MZMSH0) C NAMELIST/STGBF/IPRINT,IBUT,IFLAGV,NDRMET,IGAUGE,TOL0V,IRDEC,QETEST X ,E0,EINCR,MXE,EWIDTH,AC,DELTA C C CX DATA AC,TINY,DELTA/2*1.E-4,1.E-2/ C C FPAAS = 4 * PI * ALPHA * A0**2 / MB DATA FPAAS/2.56789/, FILE3/'DVEC'/ C ALF3=ALPHA**3 DATA ALF3/3.88593E-7/ DATA NUM /'0','1','2','3','4','5','6','7','8','9'/ C C PICK-UP FROM BLOCK DATA C AC0=AC DELTA0=DELTA C C C OPEN(UNIT=5,FILE='dstgbf',FORM='FORMATTED') OPEN(UNIT=6,FILE='routbf',FORM='FORMATTED',STATUS='UNKNOWN') C WRITE(6,598)MZCHF,MZTAR,MZPTS,MZEST,MZMSH,MZMNP,MZTET, XMZDEG,MZSA1,MZSA2,MZSA3 C C C N.B. THE "MODE" OF OPERATION IN STGBF IS SET BACK IN STGF, SEE COMMENTS C THERE, BUT BASICALLY IQDT=0 IS THE ORIGINAL MODE OF OPERATON (NON-MQDT) C WITH THE ONLY QDT BEING THE GAILITIS AVERAGE. IQDT=1 IS AN MQDT MODE C WHICH TREATS **ALL** CLOSED CHANNELS AS OPEN. HENCE STGF SHOULD BE RUN C ON A COARSE ENERGY MESH THEN STGBF RUN ON A (MUCH) FINER ENERGY MESH AS C SPECIFIED BELOW BY E0, EINCR, MXE AS PER IMESH=1 OF STGF. C C IF THE NEXT 3 VARIABLES ARE NOT SET THEN THE COARSE MESH OF STGF IS USED. C E0 -- IS THE STARTING ENERGY IN Z-SCALED RYDBERGS. C EINCR -- IS THE INCREMENTAL ENERGY. C MXE -- IS THE NUMBER OF ENERGY MESH POINTS. C C IGAUGE .EQ. 0 USE LENGTH GAUGE FOR PR OMEGPR OUTPUT (DEFAULT). C .NE. 0 USE VELOCITY GAUGE FOR PHOTORECOMBINATION INSTEAD. C NDRMET .EQ. THE NUMBER OF INITIAL ION STATES FOR WHICH RECOMBINATION C (TO ATOM) IS CALCULATED. DEFAULT=1. C IRDEC .EQ. 0 NO RADIATIVE DECAYS IN I/QDT MODE (DEFAULT). C .NE. 0 INCLUDE RADIATIVE DECYAS IN I/QDT MODE. C .EQ. 1 BELL & SEATON C .EQ. 2 HICKMAN-ROBICHEAUX (P/MQDT ONLY) C EWIDTH .GT. 0 PRECONVOLUTE **TOTAL** PHOTORECOMB/IONIZATION WITH C LORENTZIAN OF THIS WIDTH (Z-SCALED RYD). SHOULD BE C AT LEAST 6 TIMES EINCR. C IF IRDEC.NE.1 USES HICKMAN-ROBICHEAUX (ELSE BELL&SEATON). C .LT. 0 DO NOT PRECONVOLUTE (DEFAULT). C .EQ. 0 RUNS THRU PRECONVOLUTION WITH ZERO WIDTH (CHECK). C IBUT .EQ. 0 NO BUTTLE CORRECTION TO ORBITALS (DEFAULT). C .NE. 0 APPLY BUTTLE CORRECTION TO ORBITALS. C IFLAGV .EQ. 0 SKIP VMAT INTEGRAL FROM RTWO TO INFINITY. C .NE. 0 KEEP VMAT INTEGRAL (DEFAULT). C TOL0V .EQ. SMALLEST NON-VANISHING THETA FUNCTION RETAINED BY OP/VMAT C SHOULD BE COMPARABLE WITH SMALLEST NUMBER REPRESENTABLE. C HOWEVER, NOT CRUCIAL SO IF UNSURE USE 1.D-70. C IPRINT -- SEE ABOVE I/O CHANNELS C DELTA .GT. 0 SWITCH-OFF USE OF ACCELERATION INTEGRAL WHEN NU1-NU2 C IS LESS THAN DELTAN (DEFAULT = 0.01) C AC .GT. 0 ACCURACY PARAMETER C.F. STGF (DEFAULT = 1.E-4). C C IPRINT=0 IBUT=0 IFLAGV=1 NDRMET=1 IGAUGE=0 TOL0V=1.D-150 IRDEC=0 QETEST=1.D-10 IMESH=1 !CURRENTLY, THE ONLY OPTION E0=-ONE EINCR=-ONE MXE=-1 EWIDTH=-ONE DELTA=-ONE AC=-ONE C READ(5,STGBF) C IF(DELTA.LE.TZERO)DELTA=DELTA0 IF(AC.LE.TZERO)AC=AC0 ACJWBK=(6.*AC)**.2 C IF(NDRMET.GT.MZMET)THEN WRITE(6,*)' ***DIMENSION ERROR, NDRMET=',NDRMET,' >MZMET',MZMET STOP ENDIF C WRITE(6,504) XIPRINT,NDRMET,IGAUGE,IRDEC,IMESH,IFLAGV,IBUT,AC,DELTA C IF(IBUT.NE.0)WRITE(6,500)IBUT ISAVE=MZSA2 SAVE = IPRINT .GT. -2 .AND. ISAVE .EQ. 1 IF(IPRINT.EQ.3.AND.MZSA3.EQ.1)THEN WRITE(6,*)' RECOMPILE WITH MZSA3=MZMSH FOR IPRINT=3' STOP ENDIF C FLAGL=.TRUE. IF(IPRINT.LT.-2)FLAGL=IGAUGE.EQ.0 C CALL RBF00(EMESH,MXE2) CALL RD00 C IF(NZED.EQ.NELC)THEN IF(IQDT.NE.0)THEN WRITE(6,*)' IQDT OPTION NOT CODED FOR NEUTRALS YET' STOP ENDIF TZED=0. ELSE TZED=1. ENDIF C IF(IQDT.NE.0.AND.IPRINT.EQ.3)THEN WRITE(6,*)' IQDT OPTION NOT CODED FOR IPRINT=3 YET' STOP ENDIF IF(IQDT.LT.0.AND.IPRINT.EQ.-1)THEN WRITE(6,*)' THIS COMBINATION OF IQDT AND IPRINT IS NOT' X ,' MEANINGFUL' STOP ENDIF C IF(IQDT.GT.0)THEN CALL MESH(E0,EINCR,MXE,EMSH0,MXE0,EMESH,MXE2,EPS) OPEN(70,STATUS='SCRATCH',FORM='UNFORMATTED') ENDIF C IFLAGBP=0 SIG0=FPAAS*AZP2*0.3333333333333333 OMEG0=ALF3*AZ2*0.3333333333333333 IPERT0=IPERT C IRECL=8*MZPTS*MZCHF OPEN(4,STATUS='SCRATCH',ACCESS='DIRECT',RECL=IRECL) OPEN(7,FILE='XSECTN',STATUS='UNKNOWN') REWIND(7) C IF(IPRINT.EQ.3)THEN OPEN(9,FILE='BETA',STATUS='UNKNOWN') M=30 N=40 DO K=1,6 M=M+1 NAME='CROSS'//NUM(K) OPEN(M,FILE=NAME,FORM='FORMATTED',STATUS='UNKNOWN') N=N+1 NAME='BETA0'//NUM(K) OPEN(N,FILE=NAME,FORM='FORMATTED',STATUS='UNKNOWN') ENDDO C OPEN(81,FILE='CROSSL',FORM='FORMATTED',STATUS='UNKNOWN') OPEN(82,FILE='CROSSV',FORM='FORMATTED',STATUS='UNKNOWN') OPEN(83,FILE='CROSSLV',FORM='FORMATTED',STATUS='UNKNOWN') OPEN(91,FILE='BETAL',FORM='FORMATTED',STATUS='UNKNOWN') OPEN(92,FILE='BETAV',FORM='FORMATTED',STATUS='UNKNOWN') OPEN(93,FILE='BETALV',FORM='FORMATTED',STATUS='UNKNOWN') ENDIF C IF(IPRINT.LE.-2) THEN OPEN(18,FILE='OMEGPR',STATUS='UNKNOWN') DO M=1,MZMSH DO I=1,MZMET OMEGPR(I,M)=TZERO ENDDO ENDDO IRECL=8*(2*NAST+1) ELSE IRECL=8*MZEST ENDIF OPEN(8,STATUS='SCRATCH',ACCESS='DIRECT',RECL=IRECL) C CX CALL BLOCK C C SET OR CLEAR FLAGS FOR PROCESSING ALL INITIAL STATES C EOF=.FALSE. CALL RD5L1(IS1,IL1,IP1,M11,M12) ALL1=EOF ALL2=.TRUE. IF(EOF) THEN WRITE(6,501) ELSE CALL RD5L2 C RETURNING NEW VALUES ALL2,EOF. BACKSPACE 5 BACKSPACE 5 ENDIF IF(ALL2) WRITE(6,502) C C START PROCESSING SLPI CASES, SET FIRST ENTRY FLAGS C N.B. PERT REQUIRES FLAGV=.TRUE. C OPEN(3,FILE=FILE3,STATUS='UNKNOWN',FORM='UNFORMATTED') REWIND(3) C FLAGV = IPRINT .NE. -1 IF(IPRINT.LT.-2.AND.FLAGL)FLAGV=.FALSE. IF(IQDT.GT.0.AND.IPERT.NE.0)THEN IFLAGV=0 WRITE(6,*)' INTEGRALS RTWO-INFINITY DROPPED, NOT YET CODED' X ,' FOR MQDT MODE (IQDT=1 & IPERT=1)' ENDIF C WRITE(6,503) RED=.TRUE. MORE1=.TRUE. IF(IPRINT.EQ.3)THEN C CALL FACTORIALS: FAC(1,1)=1 FAC(2,1)=1 DO I=2,IFAC FAC(1,I)=-FAC(1,I-1) FAC(2,I)=(I-1)*FAC(2,I-1)/16 ENDDO ITOR=IFAC*2 ENDIF C C START LOOP FOR INITIAL SLPI - AND CHECK ON CURRENT BETA (WE'93MAR9) C 1 CALL GET1(FILE1,IS1,IL1,IP1,M11,M12) IF(.NOT.OK) GOTO 5 MORE2=.TRUE. C C READ E1-INDEPENDENT DATA C CALL R1EIND(FILE1,M11,M12,M13,MNP2,NCHF1) C IF(M13 .EQ. 0) GO TO 5 IF(IS1.NE.0)THEN SIG1=SIG0/DBLE(2*IL1+1) OMEG1=IS1*OMEG0 ENDIF IF(IS1.EQ.0)THEN SIG1=SIG0/DBLE(IL1+1) OMEG1=OMEG0 IFLAGBP=1 ENDIF C STORE FUNCTIONS ETC FOR ALL INITIAL BOUND LEVELS MLAST=M13 KPMAX=0 DO 3 M1=M11,M13 C C READ E1-DEPENDENT DATA C CALL R1EDEP(M1,MNP2,NCHF1) C C WARNING IF LEVEL IS ABOVE IONIZATION LIMIT IF(EE1(M1).LE.0.0) GO TO 3 WRITE(6,515) C C BUT FOR PRODUCTION WORK WE SKIP IT IF(IPRINT.NE.-1) GO TO 3 MLAST=M1-1 IF (MLAST .LT. M11) GO TO 5 GO TO 300 3 KPMAX=MAX(KPMAX,KPMM(M1)) C 300 IF (MLAST .GT. MZEST) THEN WRITE(6,605) MLAST CALL ABORTT(10) END IF NCHND1=NCHF1 !SINCE HAVEN'T CALLED RDV YET CALL FLAG1(M11,MLAST) RED=.TRUE. MORE2=.TRUE. FIRST=.TRUE. DO M1=M11,MLAST DO M2=1,MZSA2 TOT(M2,M1)=TZERO ENDDO ENDDO C C START LOOP FOR FINAL SLPI C 2 CALL GET2(FILE2,IS1,IL1,IP1,IS2,IL2,IP2) IF(.NOT.OK) GOTO 4 CALL GET3DV(IS1,IL1,IP1,IS2,IL2,IP2) IF(.NOT.OK) GOTO 4 CALL RDV(3) C C CHECK ACCORD BETWEEN FILES B AND D C IF(IS1.NE.NSPND.OR.IL1.NE.LRGLD1.OR.IP1.NE.NPTYD1.OR.MNP2.NE. 1 MNP2D1.OR.NCHF1.NE.NCHND1) THEN WRITE(6,510) FILE1,IS1,IL1,IP1,MNP2 CALL ABORTT(8) ENDIF C C READ E2-INDEPENDENT DATA C 15 CALL R2EIND(FILE2) C IF(IPRINT.LE.-2) CALL TARG2(NCHND2) C C C CALCULATE CROSS SECTIONS C ************************ C WRITE(6,609)IS1,IL1,IP1,M11,MLAST,IS2,IL2,IP2 IF(IPRINT.EQ.0.OR.IPRINT.EQ.1) THEN WRITE(6,610) ELSE IF(IPRINT.EQ.2) THEN WRITE(6,611) ENDIF C IF(IPERT.NE.0) CALL PERTAR C C IF IQDT>0 MODE THEN FIRST GENERATE UNPHYSICAL MATRICES AT C THE ORIGINAL FINAL EMSH0(M,M=1,MXE0) ENERGIES C IF(IQDT.GT.0)THEN REWIND(70) NCHOLD=-1 NOQDT=.TRUE. QJUMP=.FALSE. EMIN0=EMESH(1)-EPS C DO M2=1,MXE0 SKIPPY=.FALSE. IQMSH(M2)=0 EJUMP=EMSH0(M2).LT.EMIN0 ETOT=EMSH0(M2) C CALL NIGEL(SAVE,IPERT0,QJUMP,NCHOLD,NOQDT,MLAST,FIRST X ,SIG1,OMEG1,EMSH0,MXE0,EJUMP) C IF(SKIPPY)IQMSH(M2)=M2 ENDDO C C REPOSITION SCRATCH FILE FOR NEXT CALL NIGEL. REWIND(70) ENDIF C C NOW START LOOP FOR FINAL ENERGY MESH FOR PHYSICAL DATA C NCHOLD=-1 NOQDT=.TRUE. QJUMP=IQDT.GT.0 EQSAVE=-99999. ETOT=-ONE EJUMP=.FALSE. C IF(QJUMP)CALL READXY(ETOT,QETEST,NCHND2,FLAGL,FLAGV,MFIN,IPERT0) C DO M2=1,MXE2 C ETOT=EMESH(M2) IF(QJUMP.AND.ETOT-EQSAVE.GT.QETEST)THEN CALL READXY(ETOT,QETEST,NCHND2,FLAGL,FLAGV,MFIN,IPERT0) EQSAVE=ETOT ENDIF C CALL NIGEL(SAVE,IPERT0,QJUMP,NCHOLD,NOQDT,MLAST,FIRST X ,SIG1,OMEG1,EMESH,MXE2,EJUMP) ENDDO C C END LOOP FOR FINAL ENERGY MESH C WRITE(6,503) CLOSE(2) FIRST=.FALSE. 4 IF(MORE2) GO TO 2 C C END LOOP FOR FINAL SLPI C C ARCHIVE IF(IPRINT.GE.-2.AND..NOT.FIRST)THEN DO M1=M11,MLAST CALL WARCHV(IPRINT,MLAST,IS1,IL1,IP1,M1,EE1(M1),TOT(1,M1) X ,EMESH,MXE2,IFLAGBP) ENDDO ENDIF IF(.NOT.CMPLT) WRITE(6,705)IS1,IL1,IP1 CLOSE(1) C C BETA PARAMETERS, KAB MARCH 1991 C IF(IPRINT.EQ.3)CALL BETA(NZED,NELC,ENAT,NAST,EMESH,MXE2,EE1,MLAST) 5 IF(MORE1) GOTO 1 C C END LOOP FOR INITIAL SLPI C IF(IPRINT.LE.-2) THEN WRITE(18,*)NZED,NELC WRITE(18,*)NAST,MXE2,NDRMET IF(IFLAGBP.EQ.0)WRITE(18,*)(ISAT(I),LAT(I),I=1,NAST) IF(IFLAGBP.NE.0)WRITE(18,*)(0,(LAT(I)+1)*(-1)**ISAT(I),I=1,NAST) WRITE(18,710)(AZP2*ENAT(I),I=1,NAST) DO IE=1,MXE2 WRITE(18,700)EMESH(IE),(OMEGPR(J,IE),J=1,NDRMET) ENDDO WRITE(6,702) ELSE WRITE(6,703) ENDIF C WRITE(6,697) CLOSE(7) C C SUN DUM=DTIME(TARRY) TIME=TARRY(1) C CRAY CRAY CALL SECOND(TIME) C TIME=TIME/60.0 WRITE(6,999) TIME 999 FORMAT(//1X,'CPU TIME=',F9.3,' MIN',5X) C C STOP C C FORMATS C ******* C 500 FORMAT(//' RUN WITH IBUT = ',I1/5X,17('*')//) 501 FORMAT(//20X,'FROM ALL INITIAL STATES AVAILABLE') 502 FORMAT(//20X,'TO ALL ALLOWED FINAL STATES') 503 FORMAT(//1X,79('+')/) 504 FORMAT(//5X,'DATA READ FROM UNIT 5:'// 1 10X,'IPRINT = ',I2/ 2 10X,'NDRMET = ',I2/ 3 10X,'IGAUGE = ',I2/ 4 10X,'IRDEC = ',I2/ 5 10X,'IMESH = ',I2/ 6 10X,'IFLAGV = ',I2/ 7 10X,'IBUT = ',I2/ 8 10X,'AC = ',1PE7.1/ 9 10X,'DELTA = ',E7.1//) 510 FORMAT(/'0*** MISMATCH'/4X,'FILE ',A3,' HAS IS,IL,IP,MNP2 = ', * 3I3,I5//) 515 FORMAT(/'0*** HIGHER LEVELS OF THE PRESENT SERIES ARE ', 1 'ABOVE THE IONIZATION LIMIT ***'//) 598 FORMAT(///2(/1X,79('+'))///28X,'PROGRAM STGBF' * /28X,13('+')// * 18X,'PHOTOIONIZATION CROSS SECTION DATA'// * ///10X,'COMPILED FOR DIMENSIONS -'// + 15X,'CHANNELS MZCHF = ',I6/ + 15X,'TARGET STATES MZTAR = ',I6/ + 15X,'OUTER-REGION RADIAL POINTS MZPTS = ',I6/ + 15X,'BOUND LEVELS MZEST = ',I6/ + 15X,'ENERGY-MESH POINTS MZMSH =',I7/ + 15X,'R-MATRIX POLES MZMNP = ',I6/ + 15X,'COEFFICIENTS FOR THETA MZTET = ',I6/ + 15X,'DEGENERATE CHANNELS MZDEG = ',I6/ + 15X,'SAVE SPACE 1 MZSA1 = ',I6/ + 15X,'SAVE SPACE 2 MZSA2 = ',I6/ + 15X,'SAVE SPACE 3 MZSA3 = ',I6/ +/) 605 FORMAT(/' FINDS MLAST =',I5,' WHICH IS LARGER THAN ', 1 'MAXIMUM VALUE OF MZEST ALLOWED BY DIMENSIONS'//) 609 FORMAT(//20X,'BOUND-FREE TRANSITION DATA FOR'// 1 10X,'(IS1,IL1,IP1) = (',3I3,' )',5X,'LEVELS',I3,' TO',I3/ 2 10X,'(IS2,IL2,IP2) = (',3I3,' )'//) 610 FORMAT(/' M1',8X,'DE',10X,'E2',7X,'LENGTH',8X,'VELOCITY'/) 611 FORMAT(/' M1',8X,'DE',10X,'E2',7X,'LENGTH',8X,'VELOCITY', * 3X,'PARTIAL CROSS SECTIONS'/) 697 FORMAT(///20X,'END OF STGBF'/20X,12('-')///) 700 FORMAT(1PE14.8,6(1PE11.3)/(14X,6(E11.3))) 702 FORMAT(///10X,'PARTIAL CROSS SECTION FILE XSECTN WRITTEN') 703 FORMAT(///10X,'ARCHIVE FILE XSECTN WRITTEN') 705 FORMAT(/' ***WARNING*** TOTAL CROSS SECTION SUM FOR STATE',3I3, * ' MAY BE INCOMPLETE'/1X,79('*')/) 710 FORMAT(5E16.6) C END C C*************************************************************** SUBROUTINE ABORTT(N) C C HANDEL(!) CONDITION CODE AND ABORT PROGRAM C WRITE(6,97)N GOTO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16) N C 1 PRINT 99,'CHECK DATASETS ACCESSED FOR BOUND AND FREE DATA' GOTO 50 2 PRINT 99,'REQUESTED TRANSITION NOT ALLOWED: ERROR IN TOTAL SPIN' GOTO 50 3 PRINT 99,'REQUESTED TRANSITION NOT ALLOWED: TOTAL L OUT OF RANGE' GOTO 50 4 PRINT 99,'REQUESTED TRANSITION NOT ALLOWED: ERROR IN PARITY' GOTO 50 5 PRINT 99,'INITIAL LEVEL SPECIFICATION ERROR' WRITE(6,98) GOTO 50 6 CONTINUE C C CASE DELETED GOTO 50 7 PRINT 99,'BOUND LEVEL DATA NOT ON B FILE: LEVEL TOO HIGH' GOTO 50 8 PRINT 99,'CHECK DATASETS ACCESSED FOR BOUND AND DIPOLE DATA' GOTO 50 9 PRINT 99,'CHECK DATASETS ACCESSED FOR FREE AND DIPOLE DATA' GOTO 50 10 PRINT 99,'RECOMPILE PROGRAM WITH LARGER DIMENSIONS' GOTO 50 11 PRINT 99,'RE-RUN STGB AND STGF WITH THE SAME IPERT' GOTO 50 12 PRINT 99,'DIRECTORY FILE B00 NOT FOUND' WRITE(6,96) GOTO 50 13 PRINT 99,'FILE NOT FOUND' WRITE(6,96) GOTO 50 14 PRINT 99,'RE-RUN STGF WITH LARGER QNMAX' GOTO 50 15 PRINT 99,'DIRECTORY FILE F00 NOT FOUND' WRITE(6,96) GOTO 50 16 PRINT 99,'DIRECTORY FILE D00 NOT FOUND' WRITE(6,96) GOTO 50 C 50 CALL EXIT (0) C 96 FORMAT(//5X,40('*')/10X,'THE STANDARD INPUT FILES ARE'/ 1 12X,'B00 F00 D00'/ 2 12X,'B01 F01 D01'/ 3 12X,'B02 F02 D02'/ 4 12X,' . . . '/ 5 12X,' . . . '/ 6 12X,' . . . '/ 7 12X,'AS MANY AS NECESSARY'/5X,40('*')/) 97 FORMAT(//' *** ABORT: CN',I3) 98 FORMAT(//10X,60('*')/10X,'THE FOLLOWING DATA ARE READ IN FREE ', 1 'FORMAT'/30X,'IS1, IL1, IP1, M11, M12'/10X,'FOR BOUND STATE ', 2 'WITH IS1,IL1,IP1 AND LEVELS RANGING FROM M11 AND M12'/ 3 10X,60('*')/) 99 FORMAT(' ',A/) END C*************************************************************** C REAL*8 FUNCTION ARGC(E,L,AC) C C CALCULATES ARG(GAMMA(L+1-I/K)) -1/K -(1/K)*LN(K) - L*PI/2 C NRB: NOT ACTUALLY CALLED BY INJWBK CASE TZED=0. C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NRBZED/TZED C IF(TZED.EQ.0.)THEN ARGC=-DBLE(L)*1.570796327 RETURN ENDIF C IF(E.GT.0)GOTO 10 ARGC=-(DBLE(L)+.25)*3.141592654 RETURN C 10 FK=SQRT(E) ET=1./FK IP=L+1 P=IP PP=IP*IP C IF(AC.LT.1.E-4)GOTO 100 A1=10.*SQRT(ET)-ET*ET IF(A1.GT.PP)GOTO 20 X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) GOTO 200 20 L1=IP IP=1.+SQRT(A1) P=IP PP=IP*IP X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) L2=IP-1 DO 30 I=L1,L2 30 A=A+ATAN(ET/DBLE(I)) GOTO 200 C 100 A1=35.*ET**.25-ET*ET IF(A1.GT.PP)GOTO 120 X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) GOTO 140 120 L1=IP IP=1.+SQRT(A1) P=IP PP=IP*IP X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) L2=IP-1 DO 130 I=L1,L2 130 A=A+ATAN(ET/DBLE(I)) 140 A=A+.000396825540*FK*E*(7.*(1.-3.*X)*XP1*XP1+ C 2.*E*(1.-10.*X+5.*X*X))*XP1**(-5) C 200 A1=FK*X*X*.1667*PP IF(A1.GT.AC)GOTO 210 A=A-FK*(2.-X)*.25*PP GOTO 220 210 A=A-.5*ET*LOG(XP1) 220 A2=(P-.5)*XH A1=A2*X*X IF(A1.GT.AC)GOTO 230 A=A+A2*(1.-X*.33333333) GOTO 240 230 A=A+(P-.5)*ATAN(XH) 240 ARGC=A+FK/(12.*(1.+X)) C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE BETA(NZED,NELC,ENAT,NAST,EMESH,MXE2,EE1,MLAST) C C --- TO CALCULATE BETA ASYMMETRY PARAMETER FOR PHOTOELECTRON IN LS. C C NZED,NELC = NUCLEAR CHARGE, ELECTRONS IN FINAL ION; C ENAT(1:NAST) = STATE ENERGIES OF FINAL(?) ION (UNSCALED RYDBERGS); C EMESH(1:MXE2)= FINAL ELECTRON ENERGIES (Z**2 SCALED RYDBERGS); C EE1(1:MLAST) = INITIAL BOUND STATE ENERGIES (Z**2 SCALED RYDBERGS). C C THE FOLLOWING PARAMETERS IN /CBETA/ MUST BE DEFINED ON INPUT: C ITARG,LTARG,L2P,KJ = TARGET STATE,TARGET L/2*J,ELECTRON L,TARGET 2*K C IN EACH CHANNEL, FOR UP TO 3 FINAL SYMMETRIES (IL1-1,IL1,IL1+1). C C THE FOLLOWING PARAMETERS IN /CSTOR/ MUST BE DEFINED ON INPUT: C IS1,IL1,IP1 = INITIAL STATE SYMMETRY; M11 = FIRST BOUND STATE. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C LOGICAL FLAGV,FLAGL PARAMETER(MXPSC=MZPTS*MZCHF*MZSA1+2*MZPTS*MZCHF,IFAC=50) PARAMETER(MXMSH=MXPSC/4, IIIFAC=2*IFAC) PARAMETER(MXDUM=1+MXPSC-MXMSH*4) DIMENSION ENAT(NAST),EMESH(MXE2),EE1(MLAST) DIMENSION DLR(MZCHF,-1:1),DLI(MZCHF,-1:1), * DVR(MZCHF,-1:1),DVI(MZCHF,-1:1), * NCHOP(-1:1),MSTART(MZTAR),MPOINT(MZTAR) COMMON/CSTOR/FAC(IIIFAC),ITOR,IS1,IL1,IP1,M11,MFIN,M1,IL2,M2,FLAGV X ,FLAGL * /CBETA/ITARG(MZCHF,-1:1),LTARG(MZCHF,-1:1),L2P(MZCHF,-1:1) C ,KJ(MZCHF,-1:1) C USE SPACE IN /PSC/ FOR BUFFERING OUTPUT TO BETA FILE ON UNIT 9: C MXMSH= LARGEST NUMBER OF ENERGIES WITHOUT EXCEEDING LENGTH OF /PSC/ COMMON/PSC/BETAR(2,MXMSH),XSECN(2,MXMSH),NDUMMY(MXDUM) C COMMON/NRBZED/TZED C C FPAAS = 4 * PI * ALPHA * A0**2 / MB (!) DATA FPAAS/2.56789/, CGCOEF/0.81649658/ C CG2(J1,J2) = SQRT(FAC(J2-J1+6)*FAC(J1-J2+6)*FAC(J1+J2-2)*5. * /FAC(J1+J2+8)) * (MOD((J1+J2)/2,4)-1) * FAC((J1+J2)/2+4) * /(FAC((J2-J1)/2+4)*FAC((J1-J2)/2+4)*FAC((J1+J2)/2)*4.) C =VCC(J1,J2,4,0,0,0,FAC,ITOR) -- CLEBSCH-GORDAN! WE'93MAR5TH. C C C --- LOOP OVER INITIAL BOUND STATE ENERGY (E1, UNSCALED RYDBERGS) C WRITE HEADER TO BETA FILE C AZ2=MAX(NZED-NELC,1) AZ2=AZ2*AZ2 DO 132 M1=M11,MLAST E1=AZ2*EE1(M1) DO 1 K=1,NAST 1 MSTART(K)=0 LAST=0 C DO 5 K=6,9,3 WRITE(K,700)NZED,NELC,NAST WRITE(K,701)(ENAT(I),I=1,NAST) 5 WRITE(K,700)IS1,IL1,IP1,M1 WRITE(6,703)E1,MXE2 C C --- LOOP OVER FINAL ELECTRON ENERGY (E2, UNSCALED RYDBERGS); C PHOTON ENERGY = W1*(UNSCALED RYDBERGS) C DO 100 M2=1,MXE2 E2=AZ2*EMESH(M2) W1 = E2-E1 IF(IS1.EQ.0)THEN SIG = (EMESH(M2)-EE1(M1))*FPAAS/((IL1+1)*3*AZ2) ELSE SIG = (EMESH(M2)-EE1(M1))*FPAAS/((2*IL1+1)*3*AZ2) ENDIF C C --- FOR EACH FINAL SYMMETRY (IL2), C RETRIEVE DIPOLE VECTORS (DLR...) AND OPEN CHANNELS (NCHOP); C DO 4 IL2=IL1-1,IL1+1 4 CALL DSTOR(1,DLR(1,IL2-IL1),DLI(1,IL2-IL1), * DVR(1,IL2-IL1),DVI(1,IL2-IL1),NCHOP(IL2-IL1)) C C --- LOOP OVER FINAL ION STATES (NA) C DO 21 NA=1,NAST IF(E2.LE.ENAT(NA)) GO TO 21 C E2P=E2-ENAT(NA)*(UNSCALED RYDBERGS) = CHANNEL ENERGY! XSECNL=0. XSECNV=0. BETALR=0. BETALI=0. BETAVR=0. BETAVI=0. C C PICK OUT COUPLED CHANNELS (I) C DO 17 IL2=IL1-1,IL1+1 IF(NCHOP(IL2-IL1).EQ.0) GO TO 17 DO 16 I=1,NCHOP(IL2-IL1) IF(ITARG(I,IL2-IL1).NE.NA) GO TO 16 C C EVALUATE CROSS SECTION (XSECNL,XSECNV) FOR IONIC STATE C XSECNL = (DLR(I,IL2-IL1)**2+DLI(I,IL2-IL1)**2)*SIG+XSECNL IF(FLAGV) *XSECNV = (DVR(I,IL2-IL1)**2+DVI(I,IL2-IL1)**2)*SIG+XSECNV C C EVALUATE BETA (BETALR,BETALI;BETAVR,BETAVI) FOR IONIC STATE C LA=L2P(I,IL2-IL1) LC=2*LA LR=2*IL2 LT=2*LTARG(I,IL2-IL1) IF(IS1.EQ.0)THEN LR=IL1+2*(IL2-IL1) LT=LTARG(I,IL2-IL1) KR=KJ(I,IL2-IL1) ENDIF ZETO = (NELC-NZED)/SQRT(E2-ENAT(NA)) H=LA+1 CALL COULGA(H,ZETO,CO,EPHA) C C PICK OUT COUPLED CHANNELS (J) C DO 15 JL2=IL1-1,IL1+1 IF(NCHOP(JL2-IL1).EQ.0) GO TO 15 DO 14 J=1,NCHOP(JL2-IL1) IF(ITARG(J,JL2-IL1).NE.NA) GO TO 14 LB=L2P(J,JL2-IL1) LD=2*LB IF(ABS(LC-LD).GT.4) GO TO 14 IF(LC+LD.EQ.0) GO TO 14 IF(MOD(LC+LD,4).NE.0) GO TO 14 IF(LC+LD+8.GT.ITOR) THEN M=(LC+LD)/2+4 PRINT *,' BETA REQUIRES FACTORIAL ARRAY OF SIZE IFAC = ',M STOP ENDIF H=LB+1 CALL COULGA(H,ZETO,CO,EPHB) LS=2*JL2 C C COMBINE CHANNELS C LN=2*IL1 IF(IS1.EQ.0)THEN LS=IL1+2*(JL2-IL1) LN=IL1 KS=KJ(J,JL2-IL1) ENDIF H = (LC+1)*(LD+1)*(LR+1)*(LS+1) C EXTRA=SQRT(H)*CGCOEF * VCC(LC,LD,4,0,0,0,FAC,ITOR) -- REPLACED: IF(IS1.EQ.0)THEN H=H*(KR+1)*(KS+1) MSIGN=1-MOD(2+KR+KR+KS+LT,4) EXTRA=MSIGN*SQRT(H)*CGCOEF * CG2(LC,LD) * * WRACAH(4,LR,2,IL1,LS,2,FAC,ITOR) * * WRACAH(4,KR,LS,1,KS,LR,FAC,ITOR) * * WRACAH(KR,LC,KS,LD,LT,4,FAC,ITOR) ELSE EXTRA=SQRT(H)*CGCOEF * CG2(LC,LD) * * WRACAH(LR,LC,LS,LD,LT,4,FAC,ITOR) * * WRACAH(2,LR,2,LS,LN,4,FAC,ITOR) IF(MOD(IL2,2).NE.0)EXTRA=-EXTRA !NEW IF(MOD(JL2,2).NE.0)EXTRA=-EXTRA !NEW ENDIF H = EPHB-EPHA COSR=COS(H) SINI=SIN(H) RR=COSR*DLR(I,IL2-IL1)-SINI*DLI(I,IL2-IL1) RT=COSR*DLI(I,IL2-IL1)+SINI*DLR(I,IL2-IL1) BETALR = (RR*DLR(J,JL2-IL1)+RT*DLI(J,JL2-IL1))*EXTRA + BETALR BETALI = (RT*DLR(J,JL2-IL1)-RR*DLI(J,JL2-IL1))*EXTRA + BETALI IF(.NOT.FLAGV) GO TO 14 RS=COSR*DVR(I,IL2-IL1)-SINI*DVI(I,IL2-IL1) RU=COSR*DVI(I,IL2-IL1)+SINI*DVR(I,IL2-IL1) BETAVR = (RS*DVR(J,JL2-IL1)+RU*DVI(J,JL2-IL1))*EXTRA + BETAVR BETAVI = (RU*DVR(J,JL2-IL1)-RS*DVI(J,JL2-IL1))*EXTRA + BETAVI 14 CONTINUE 15 CONTINUE 16 CONTINUE 17 CONTINUE C C DIVIDE BY CROSS SECTION AND MULTIPLY BY CONSTANT FACTOR C MSIGN = 1-MOD(LT+IL1+IL1,4) IF(IS1.EQ.0)MSIGN=1 IF(XSECNL.GT.0.0)THEN !new H=MSIGN*3*SIG/XSECNL ELSE !new H=0.0 !new ENDIF !new BETALR=BETALR*H BETALI=BETALI*H IF(FLAGV) THEN IF(XSECNV.GT.0.0)THEN !new H=MSIGN*3*SIG/XSECNV ELSE !new H=0.0 !new ENDIF !new BETAVR=BETAVR*H BETAVI=BETAVI*H ENDIF C C MSTART(NA)= FIRST ENERGY IN EMESH WITH OPEN CHANNELS FOR NA C MPOINT(NA) = POINTER FOR ENERGY RANGE IN BETA AND XSEC FOR NA C IF(MSTART(NA).EQ.0) THEN MPOINT(NA)=LAST LAST=LAST+(MXE2-M2+1) IF(LAST.GT.MXMSH) GO TO 20 MSTART(NA)=M2 ENDIF MP=MPOINT(NA) + M2-MSTART(NA)+1 BETAR(1,MP)=BETALR BETAR(2,MP)=BETAVR XSECN(1,MP)=XSECNL XSECN(2,MP)=XSECNV C 20 WRITE(6,704) W1,XSECNL,XSECNV,BETALR,BETAVR,BETALI,BETAVI,NA 21 CONTINUE 100 CONTINUE C C --- END OF LOOPS. WRITE BETA ON UNIT 9 C DO 221 K=1,NAST IF(MSTART(K).EQ.0) GO TO 221 WRITE(9,702) E1,MXE2-MSTART(K)+1,ENAT(K),K DO 200 I=MSTART(K),MXE2 M = MPOINT(K) + I-MSTART(K)+1 WRITE(9,704)AZ2*EMESH(I), * XSECN(1,M),XSECN(2,M),BETAR(1,M),BETAR(2,M) IF(K.EQ.1)THEN ERYD=AZ2*(EMESH(I)-EE1(1)) CT EEV=13.6058*ERYD WRITE(91,135)ERYD,BETAR(1,M) WRITE(92,135)ERYD,BETAR(2,M) WRITE(81,135)ERYD,XSECN(1,M) WRITE(82,135)ERYD,XSECN(2,M) WRITE(83,135)ERYD,XSECN(1,M),XSECN(2,M) WRITE(93,135)ERYD,BETAR(1,M),BETAR(2,M) ENDIF ERYD=AZ2*(EMESH(I)-ENAT(K)) CT EEV=13.6058*ERYD IF(K.LE.6)THEN MFILE=30+K NFILE=40+K WRITE(MFILE,135)ERYD,XSECN(1,M),XSECN(2,M) WRITE(NFILE,135)ERYD,BETAR(1,M),BETAR(2,M) ENDIF 200 CONTINUE 135 FORMAT(3F13.6) 221 CONTINUE 132 CONTINUE C 700 FORMAT(4I5) 701 FORMAT(5E14.6) 702 FORMAT(2(E14.6,I5)) 703 FORMAT(E14.6,I5/' E(PHOTON)/RYD XSECL/MB XSECV/MB', *' RE(BETAL) RE(BETAV) IM(BETAL) IM(BETAV) STATE') 704 FORMAT(1P,E14.6,6E10.2,I6) RETURN END C****************************************************************** C CX SUBROUTINE BLOCK C C PROVIDES BLOCK DATA C QUB90'FEB2/91MAR23, WE: BECAUSE OF IBM COMPILERS NOT CALLED AS CX SUBROUTINE (BLOCK TO AVOID LINKAGE PROBLEMS WITH LIBRARIES) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMON /CTRLBF/IPRINT,IBUT,TINY,DELTA,AC DATA AC,TINY,DELTA/2*1.E-4,1.E-2/ C 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 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 CX RETURN END C+++++++++++++++++++++++++++++++++++++++++++++ C******************************** SUBROUTINE COULGA(A1,B,C,D) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMON/NRBZED/TZED C C COULOMB PHASE C B=-1.0/K C C=MODULUS OF GAMMA FUNCTION C D=ARGUEMENT OF GAMMA FUNCTION C IF(TZED.EQ.0.)THEN C=0. D=0. RETURN ENDIF C C1=0.0 D1=0.0 K1= INT(A1) IF(K1.GT.9) GO TO 2 A2=A1-K1 DO 1 J=K1,9 A=A2+J D = ATAN2(B,A) A = LOG(A*A+B*B)*.5 C1=A+C1 1 D1=D+D1 A=10.0 +A2 GO TO 3 2 A=A1 3 C = LOG(A*A+B*B)*.5 D= ATAN2(B,A) F=C C = (A-0.5)*C-D*B D = (A-0.5)*D+F*B C C = C-A + 0.5*LOG(6.28318531) C = C-A + 0.918938534 D=D-B E=A*A+B*B A=A/E E=-B/E F=A*A-E*E G=2.0 *E*A H = 0.000841750842*F-0.000595238095 H1= 0.000841750842*G H2=H H=H*F-H1*G H1=H1*F+H2*G H=H+0.000793650794 H2=H H=H*F-H1*G H1=H1*F+H2*G H=H-0.00277777778 H2=H H=H*F-H1*G H1=H1*F+H2*G H=H+0.0833333333 C = EXP(H*A-H1*E+C-C1) D = H1*A+H*E+D-D1 RETURN END C************************************************ SUBROUTINE DA2(KEY,IREC,JDISC,LENGTH,ARRAY) C C TO STORE A LARGE ARRAY IN A DA FILE. CRAY VERSION. C C KEY = 1 FOR READ, = 2 FOR WRITE, C = 0 FOR FINDING NUMBER OF DA RECORDS A GIVEN ARRAY TAKES. C IREC= (ON CALL) POINTER TO FIRST DA RECORD FOR ARRAY, C = 0 FOR OPENING DA FILE (BY NAME), C =-1 FOR OPENING DA FILE (SCRATCH), C = (ON RETURN) POINTER TO NEXT AVAILABLE DA RECORD. C JDISC = DA FILE UNIT NUMBER. C ARRAY(LENGTH) = ARRAY TO READ OR WRITE. C IMPLICIT REAL*8 (A) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (LREC=512) DIMENSION ARRAY(LENGTH) C IF(IREC.GT.0) GO TO 10 IRECL=8*LREC IF(IREC.LT.0) THEN OPEN(JDISC,STATUS='SCRATCH',ACCESS='DIRECT',RECL=IRECL) GO TO 9 ENDIF IF(KEY.EQ.2) THEN OPEN(JDISC,STATUS='NEW',FILE='RK',ACCESS='DIRECT',RECL=IRECL) ENDIF IF(KEY.EQ.1) THEN OPEN(JDISC,STATUS='OLD',FILE='RK',ACCESS='DIRECT',RECL=IRECL) ENDIF 9 IREC=1 C 10 IF(LENGTH.EQ.0) RETURN I2=0 20 I1=I2+1 I2=MIN(I2+LREC,LENGTH) IF(KEY.EQ.0) GO TO 30 IF(KEY.EQ.2) THEN WRITE(JDISC,REC=IREC) (ARRAY(I),I=I1,I2) ELSE READ(JDISC,REC=IREC) (ARRAY(I),I=I1,I2) ENDIF 30 IREC=IREC+1 IF(I2.LT.LENGTH) GO TO 20 RETURN END C****************************************************************** C SUBROUTINE DIN(E2,NCHOP,M11,MFIN,DLI,DVI) C C CALCULATE INNER REGION CONTRIBUTION TO THE DIPOLE LENGTH AND C VELOCITY MATRIX ELEMENTS C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (TZERO=0.0) C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF) X,RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/RBASIS/VALUE(MZMNP),WMAT(MZCHF,MZMNP) COMMON/SCALE/AZ,AZP,AZ2,AZP2,K10,K20 COMMON/DVECT/DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST) C DIMENSION GP(MZCHF),AVECT2(MNPEXT) DIMENSION DLI(MZCHF,MZEST),DVI(MZCHF,MZEST) C DO JCHOP=1,NCHOP C C CALCULATE AVECT2 C DO N=1,NCHND2 GP(N)=FP2(N,JCHOP)-BSTO*F2(N,JCHOP) AVECT2(MNP2D2+N)=AZ*GP(N) ENDDO DO K=1,MNP2D2 AVECT2(K)=TZERO ENDDO C CCRAY CRAY DO N=1,NCHND2 CRAY DO K=1,MNP2D2 CRAY AVECT2(K)=AVECT2(K)+WMAT(N,K)*GP(N) CRAY ENDDO CRAY ENDDO CCRAY C DO K=1,MNP2D2 DO N=1,NCHND2 AVECT2(K)=AVECT2(K)+WMAT(N,K)*GP(N) ENDDO ENDDO C DO K=1,MNP2D2 AVECT2(K)=AVECT2(K)/(VALUE(K)-E2) ENDDO C C CHANGE F2 AND DERIVATIVES FROM 1ST ORDER TO 0TH ORDER FUNCTIONS C IF(IPERT.GT.0) THEN DO N=1,NCHND2 F2(N,JCHOP)=F2(N,JCHOP)-DF2(N,JCHOP) FP2(N,JCHOP)=FP2(N,JCHOP)-DFP2(N,JCHOP) FPP2(N,JCHOP)=FPP2(N,JCHOP)-DFPP2(N,JCHOP) ENDDO ENDIF C C CALCULATE INNER REGION CONTRIBUTION, DLI AND DVI C DO M1=M11,MFIN DLI(JCHOP,M1)=TZERO DVI(JCHOP,M1)=TZERO ENDDO DO M1=M11,MFIN DO K2=1,K20 DVI(JCHOP,M1)=DVI(JCHOP,M1)+AVECT2(K2)*DVECTV(K2,M1) DLI(JCHOP,M1)=DLI(JCHOP,M1)+AVECT2(K2)*DVECTL(K2,M1) ENDDO ENDDO C ENDDO C RETURN END C********************************************************************* C SUBROUTINE DMAT(M1) C C SETS UP ARRAYS FOR PERTURBATION CALCULATIONS C C NRB: C VECTOR OR SCALAR VERSIONS. C N.B. SEPARATING LENGTH AND VELOCITY ONLY GIVES A SMALL TIME SAVING C (WHEN EVALUATING ONE OR THE OTHER) NOT A FACTOR OF 2. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C PARAMETER (MNPEXT=MZMNP+MZCHF) COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CPERT/CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2) X,DMATV(MZCHF,MZCHF,2) 1 ,WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) CCRAY CRAY DO L=1,2 CRAY DO J=1,NCHND2 CRAY DO I=1,NCHND1 CRAY DMATL(I,J,L)=TZERO CRAY DMATV(I,J,L)=TZERO CRAY ENDDO CRAY ENDDO CRAY DO J1=1,NCHND1 CRAY DO I2=1,NCHND2 CRAY D=ONE/(EPS1(J1,M1)-EPS2(I2)) CRAY DO I1=1,NCHND1 CRAY DMATL(I1,I2,L)=DMATL(I1,I2,L)+WMAT1(I1,J1,L)*BMATL(J1,I2)*D CRAY DMATV(I1,I2,L)=DMATV(I1,I2,L)+WMAT1(I1,J1,L)*BMATV(J1,I2)*D CRAY ENDDO CRAY ENDDO CRAY ENDDO CRAY DO J2=1,NCHND2 CRAY DO I1=1,NCHND1 CRAY D=ONE/(EPS1(I1,M1)-EPS2(J2)) CRAY DO I2=1,NCHND2 CRAY DMATL(I1,I2,L)=DMATL(I1,I2,L)-BMATL(I1,J2)*WMAT2(J2,I2,L)*D CRAY DMATV(I1,I2,L)=DMATV(I1,I2,L)-BMATV(I1,J2)*WMAT2(J2,I2,L)*D CRAY ENDDO CRAY ENDDO CRAY ENDDO CRAY ENDDO CCRAY C DO L=1,2 DO I2=1,NCHND2 DO I1=1,NCHND1 DMATL(I1,I2,L)=TZERO DMATV(I1,I2,L)=TZERO ENDDO DO J1=1,NCHND1 D=ONE/(EPS1(J1,M1)-EPS2(I2)) DL=BMATL(J1,I2)*D DV=BMATV(J1,I2)*D DO I1=1,NCHND1 DMATL(I1,I2,L)=DMATL(I1,I2,L)+WMAT1(I1,J1,L)*DL DMATV(I1,I2,L)=DMATV(I1,I2,L)+WMAT1(I1,J1,L)*DV ENDDO ENDDO DO J2=1,NCHND2 DO I1=1,NCHND1 D=ONE/(EPS1(I1,M1)-EPS2(J2)) DMATL(I1,I2,L)=DMATL(I1,I2,L)-BMATL(I1,J2)*WMAT2(J2,I2,L)*D DMATV(I1,I2,L)=DMATV(I1,I2,L)-BMATV(I1,J2)*WMAT2(J2,I2,L)*D ENDDO ENDDO ENDDO ENDDO C RETURN END C********************************************************************* C SUBROUTINE DMQDT(DLVECT,DVVECT,FLAGL,FLAGV,DL2,DV2,M1) C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C CBL PARAMETER (MWORK=100*MZDEG) PARAMETER (TZERO=0.0) PARAMETER (ZERO=(0.0,0.0)) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (BIG=150.0) PARAMETER (TINE=1.0E-6) C LOGICAL FLAGV,QDT,FLAGL C C NRB: C EVALUATE DETAILED MULTI-STATE QDT. C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF), X RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAGQ/QDT COMMON/TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) C COMMON/NRBKHI/ZKHICC(MZDEG,MZDEG),ZKHIOC(MZCHF,MZDEG),ZVAL(MZDEG) X ,ZDL(MZCHF),ZDV(MZCHF),ZDL0(MZCHF),ZDV0(MZCHF) CBL X,ZWL(MZDEG,MZDEG),ZWR(MZDEG,MZDEG),ZWORK(MWORK),RWORK(2*MZDEG) COMMON/NRBZDD/ZDLM(MZCHF,MZEST),ZDVM(MZCHF,MZEST) C DIMENSION IPIV(MZDEG) DIMENSION DLVECT(MZCHF),DVVECT(MZCHF),DL2(MZCHF),DV2(MZCHF) DIMENSION XMOLD(MZCHF,MZCHF),YMOLD(MZCHF,MZCHF) C EQUIVALENCE (XMOLD,DFP2),(YMOLD,DFPP2) C CBL DATA LWMAX/MWORK/ C C C C ALL CHANNELS C NCHF=NCHND2 C C FIRST DETERMINE NCC (NUMBER OF CLOSED CHANNELS) C DO J=NCHF,1,-1 IF(EPS2(J).GE.TZERO) THEN NCC=NCHOP-J NCHOP=J GO TO 1 ENDIF ENDDO C C CASE OF NO REAL OPEN CHANNELS C NCHOP=0 RETURN C 1 IF(QDT)THEN C C DIVIDE INTO WEAK (NWC) AND STRONGLY CLOSED (NCC) C NCHOPP=NCHOP+1 DO I=NCHOPP,NCHF IF(ABS(ECH2(I)-ECH2(NCHOPP)).GT.TINE)GO TO 2 NQ=I ENDDO 2 NWC=NQ-NCHOP NCHOP=NCHOP+NWC NCC=NCC-NWC ENDIF C C CASE ALL CHANNELS OPEN (OR AT LEAST NO STRONGLY CLOSED) C IF(NCC.EQ.0)THEN IF(FLAGL)THEN DO J=1,NCHF ZDL0(J)=ZDLM(J,M1) DL2(J)=ZDL0(J)*CONJG(ZDL0(J)) ENDDO ENDIF IF(FLAGV)THEN DO J=1,NCHF ZDV0(J)=ZDVM(J,M1) DV2(J)=ZDV0(J)*CONJG(ZDV0(J)) ENDDO ENDIF RETURN ENDIF C PI=ACOS(-ONE) TPI=TWO*PI C C INITIALIZE KHICC; USE (KHI) = CONJG(2I(Y)-1) = (2(XM)-1,2(YM)) C IF(NCC.GT.MZDEG) THEN WRITE(6,601)NCC CALL ABORTT(10) ENDIF DO N2=1,NCC DO N1=1,NCC ZKHICC(N1,N2)=DCMPLX(TWO*XM(NCHOP+N1,NCHOP+N2), X TWO*YM(NCHOP+N1,NCHOP+N2)) ENDDO ZKHICC(N2,N2)=ZKHICC(N2,N2)-ONE ENDDO C C ZKHICC-ZFDEC C DO N=1,NCC NNN=NCHOP+N C C RADIATIVE DECAYS C N.B. WE CAN SUM OVER ALL FINAL STATES, NOT JUST BOUND, SINCE WE C WE ARE NOT RELYING ON UNITARITY. C IF(IRDEC.EQ.0)THEN TPINU=FNU2(NNN)*TPI ZFDEC=EXP(DCMPLX(TZERO,-TPINU)) ELSEIF(IRDEC.EQ.1)THEN !BELL&SEATON T=ARDEC(NTARG(NNN))*(FNU2(NNN)**3)/TWO T=MIN(T,BIG) FDEC=EXP(T) TPINU=FNU2(NNN)*TPI ZFDEC=FDEC*EXP(DCMPLX(TZERO,-TPINU)) ELSE !HICKMAN-ROBICHEAUX FR=ONE/FNU2(NNN)**2 FI=-ARDEC(NTARG(NNN))/TPI ZFNU=DCMPLX(ONE,TZERO)/SQRT(DCMPLX(FR,FI)) Z=DCMPLX(TZERO,-TPI)*ZFNU FR=DBLE(Z) FI=DIMAG(Z) FR=MIN(FR,BIG) ZFDEC=EXP(DCMPLX(FR,FI)) ENDIF C ZKHICC(N,N)=ZKHICC(N,N)-ZFDEC ENDDO C C INVERT C CSTRTNBL CALL ZVERT(ZKHICC,MZDEG,NCC,IPIV) CENDNBL C CSTRTBL CBL CALL ZGETRF(NCC,NCC,ZKHICC,MZDEG,IPIV,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,602) INFO CBL STOP CBL ENDIF CBL CALL ZGETRI(NCC,ZKHICC,MZDEG,IPIV,ZWORK,MWORK,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,603) INFO CBL STOP CBL ENDIF CBL LWOPT=INT(ZWORK(1)) CBL IF (LWOPT.GT.MWORK) WRITE(6,604) LWOPT CENDBL C C ZKHIOC*ZKHICC C DO N2=1,NCC DO N1=1,NCHOP ZKHIOC(N1,N2)=ZERO ENDDO DO K=1,NCC DO N1=1,NCHOP ZKHIOC(N1,N2)=ZKHIOC(N1,N2)+DCMPLX(TWO*XM(N1,NCHOP+K), X TWO*YM(N1,NCHOP+K))*ZKHICC(K,N2) ENDDO ENDDO ENDDO C C DIPOLE MATRIX: WE WILL ALREADY HAVE IT, RECOVERED AND C INTERPOLATED IN READXY. C IF(FLAGL)THEN DO J=1,NCHF ZDL(J)=ZDLM(J,M1) ENDDO ENDIF IF(FLAGV)THEN DO J=1,NCHF ZDV(J)=ZDVM(J,M1) ENDDO ENDIF C IF(.NOT.QDT)THEN C C FORM PHYSICAL DIPOLE MATRIX C IF(FLAGL)THEN DO J=1,NCHOP Z=ZDL(J) DO K=1,NCC Z=Z-ZKHIOC(J,K)*ZDL(NCHOP+K) ENDDO DL2(J)=Z*CONJG(Z) ENDDO ENDIF IF(FLAGV)THEN DO J=1,NCHOP Z=ZDV(J) DO K=1,NCC Z=Z-ZKHIOC(J,K)*ZDV(NCHOP+K) ENDDO DV2(J)=Z*CONJG(Z) ENDDO ENDIF C ELSE C C STORAGE ORIGINAL XM,YM C DO J=1,NCHF DO I=1,NCHF XMOLD(I,J)=XM(I,J) YMOLD(I,J)=YM(I,J) ENDDO ENDDO C C NOW CONTRACT C DO J=1,NCHOP DO I=1,J ZKHI=DCMPLX(XMOLD(I,J),YMOLD(I,J)) DO K=1,NCC ZKHI=ZKHI-ZKHIOC(I,K)*DCMPLX(XMOLD(J,NCHOP+K) X ,YMOLD(J,NCHOP+K)) ENDDO XM(I,J)=DBLE(ZKHI) XM(J,I)=XM(I,J) YM(I,J)=DIMAG(ZKHI) YM(J,I)=YM(I,J) ENDDO ENDDO C C FORM CONTRACTED DIPOLE MATRIX C IF(FLAGL)THEN DO I=1,NCHOP Z=ZDL(I) DO K=1,NCC Z=Z-ZKHIOC(I,K)*ZDL(NCHOP+K) ENDDO ZDL0(I)=Z DL2(I)=Z*CONJG(Z) ENDDO ENDIF IF(FLAGV)THEN DO I=1,NCHOP Z=ZDV(I) DO K=1,NCC Z=Z-ZKHIOC(I,K)*ZDV(NCHOP+K) ENDDO ZDV0(I)=Z DV2(I)=Z*CONJG(Z) ENDDO ENDIF C ENDIF C RETURN C 601 FORMAT(//' ***NUMBER OF CLOSED CHANNELS, NCC = ',I2/4X, 1 'LARGER THAN DIMENSION VALUE OF DEG = MZDEG'//) CB602 FORMAT(//10X,10('*'),' SR. DMQDT: ZGETRF RETURNED WITH INFO =',I2) CB603 FORMAT(//10X,10('*'),' SR. DMQDT: ZGETRI RETURNED WITH INFO =',I2) CB604 FORMAT(//10X,10('*'),' SR. DMQDT: MWORK SHOULD BE LARGER THAN ', CB X 'WORK(1) = ',F5.0,' FOR OPTIMAL PERFORMANCE') END C********************************************************************* C SUBROUTINE DOUT(DLA,DLB,DVA,DVB,DLP,DVP,M1,JCHOP,FLAGP,FLAGV,DER) C C CALCULATES THE OUTER REGION CONTRIBUTIONS TO THE DIPOLE LENGTH AND C VELOCITY MATRIX ELEMENTS INCLUDING, IF REQUIRED, PERTURBATION C CONTRIBUTIONS. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) C PARAMETER (TZERO=0.0) PARAMETER (HALF=0.5) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF) X,RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/CPERT/CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2) X,DMATV(MZCHF,MZCHF,2) 1 ,WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBTMP/E2 COMMON/NRBZED/TZED C LOGICAL FLAGP,FLAGV,ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT X ,SKIPPY C C NRB: C WHEN N1 APPROACHES N2, CANCELLATION ERROR CAN OCCUR (ZERO/ZERO) IN C THE "ACCELERATION" FORM OF THE INTEGRALS. IN ADDITION, WHEN N2 IS A C S+KC STATE THE WRONSKIAN AT INFINITY MAY NOT BE NEGLIGIBLE. IN THE C LATTER CASE A LARGER DELTA APPEARS TO BE NECESSARY WHEN PRECONVOLUTING. C THE REASON FOR THIS IS NOT CLEAR. THERE IS NO "PROBLEM" WHEN EWIDTH=ZERO. C DELTAM=DELTA DELTAP=DELTA C DLA=TZERO DLB=TZERO DVA=TZERO DVB=TZERO DLP=TZERO DVP=TZERO AMV=TZERO C C DO N2=1,NCHND2 DO N1=1,NCHND1 C C IF PRECONVOLUTING USE CONSTANT DELTA-E NOT DELTA-N. C IF(EWIDTH.GE.TZERO)THEN DELTAM=DELTA*(FNU1(N1,M1)/THREE)**3 IF(FNU2(N2).LT.TZERO)THEN DELTAP=-1.0E10 ELSE DELTAP=1.0E10 ENDIF ENDIF C AML=AMATL(N1,N2) BML=BMATL(N1,N2) BMV=BMATV(N1,N2) DELTAN=FNU2(N2)-FNU1(N1,M1) C C CALCULATE AMAT CONTRIBUTIONS C IF(ABS(AML).GT.TINY) THEN IF(DELTAN.LT.-DELTAM.OR.DELTAN.GT.DELTAP) THEN IF(KKP) THEN CALL S1S2(P1(N1,M1),PP1(N1,M1),CC1(N1),F2(N2,JCHOP), X FP2(N2,JCHOP),FPP2(N2,JCHOP),CC2(N2),S1,S2) ELSE CALL S1S2(F2(N2,JCHOP),FP2(N2,JCHOP),CC2(N2),P1(N1,M1), X PP1(N1,M1),PPP1(N1,M1),CC1(N1),S1,S2) ENDIF DA=VMAT(-2,N1,N2,M1)*RK(N2,JCHOP) IF(N2.EQ.JCHOP) DA=DA+UMAT(-2,N1,N2,M1) DAC=(-S2+TZED*FOUR*DA)*DER*DER DOR=-S1*DER+DAC IF(NSPND.EQ.0)AMV = AMATV(N1,N2) C = 0 IN LS COUPLING (ACCOUSNTING FOR 'REVERSE' IN BP-RECUD) IF(AMV.EQ.0.) THEN S3 = DAC*AML ELSE IF(KKP) THEN CALL S1S2(F2(N2,JCHOP),FP2(N2,JCHOP),CC2(N2),P1(N1,M1), X PP1(N1,M1),PPP1(N1,M1),CC1(N1),S1,S3) ELSE CALL S1S2(P1(N1,M1),PP1(N1,M1),CC1(N1),F2(N2,JCHOP), X FP2(N2,JCHOP),FPP2(N2,JCHOP),CC2(N2),S1,S3) ENDIF S3=(S3-DA*FOUR*TZED)*DER*DER*AMV ENDIF DVA = S3 + DVA ELSE FLAGV=.FALSE. DOR=VMAT(1,N1,N2,M1)*RK(N2,JCHOP) IF(N2.EQ.JCHOP) DOR=DOR+UMAT(1,N1,N2,M1) ENDIF DLA=DLA+AML*DOR ENDIF C C CALCULATE BMAT CONTRIBUTIONS C IF(ABS(BML).GT.TINY)THEN IF(DELTAN.LT.-DELTAM.OR.DELTAN.GT.DELTAP) THEN S3=P1(N1,M1)*FP2(N2,JCHOP)-PP1(N1,M1)*F2(N2,JCHOP) DOV=-S3/(EPS1(N1,M1)-EPS2(N2)) ELSE FLAGP=.FALSE. DOV=VMAT(0,N1,N2,M1)*RK(N2,JCHOP) IF(N2.EQ.JCHOP) DOV=DOV+UMAT(0,N1,N2,M1) ENDIF DLB=DLB+BML*DOV DVB=DVB+BMV*DOV ENDIF C ENDDO ENDDO C C C PERTURBATION CONTRIBUTIONS C IF(.NOT.FLAGV) FLAGP=.FALSE. IF(.NOT.FLAGP) THEN DLP=TZERO DVP=TZERO RETURN ENDIF C CALL DMAT(M1) C DLAP=TZERO DLBP=TZERO DVAP=TZERO DVBP=TZERO AMV=TZERO C DO N2=1,NCHND2 DO N1=1,NCHND1 C AML=AMATL(N1,N2) BML=BMATL(N1,N2) BMV=BMATV(N1,N2) CM=CMAT(N1,N2) DML1=DMATL(N1,N2,1) DML2=DMATL(N1,N2,2) DMV1=DMATV(N1,N2,1) DMV2=DMATV(N1,N2,2) C C AMAT CONTRIBUTION C IF(ABS(AML).GT.TINY) THEN IF(NSPND.EQ.0)AMV=AMATV(N1,N2) IF(AMV.EQ.0.) THEN IF(KKP) THEN CALL S1S2(P1(N1,M1),PP1(N1,M1),CC1(N1),DF2(N2,JCHOP), X DFP2(N2,JCHOP),DFPP2(N2,JCHOP),CC2(N2),S12,S22) CALL S1S2(F2(N2,JCHOP),FP2(N2,JCHOP),CC2(N2),DP1(N1,M1), X DPP1(N1,M1),DPPP1(N1,M1),CC1(N1),S11,S21) ELSE CALL S1S2(DF2(N2,JCHOP),DFP2(N2,JCHOP),CC2(N2),P1(N1,M1), X PP1(N1,M1),PPP1(N1,M1),CC1(N1),S12,S22) CALL S1S2(DP1(N1,M1),DPP1(N1,M1),CC1(N1),F2(N2,JCHOP), X FP2(N2,JCHOP),FPP2(N2,JCHOP),CC2(N2),S11,S21) ENDIF DAC=-(S21+S22)*DER*DER DOR=-(S11+S12)*DER+DAC DVAP=DVAP-AML*DAC !+/- ELSE IF(KKP) THEN CALL S1S2(DF2(N2,JCHOP),DFP2(N2,JCHOP),CC2(N2),P1(N1,M1), X PP1(N1,M1),PPP1(N1,M1),CC1(N1),S12,S22) CALL S1S2(DP1(N1,M1),DPP1(N1,M1),CC1(N1),F2(N2,JCHOP), X FP2(N2,JCHOP),FPP2(N2,JCHOP),CC2(N2),S11,S21) ELSE CALL S1S2(P1(N1,M1),PP1(N1,M1),CC1(N1),DF2(N2,JCHOP), X DFP2(N2,JCHOP),DFPP2(N2,JCHOP),CC2(N2),S12,S22) CALL S1S2(F2(N2,JCHOP),FP2(N2,JCHOP),CC2(N2),DP1(N1,M1), X DPP1(N1,M1),DPPP1(N1,M1),CC1(N1),S11,S21) ENDIF DAC=-(S21+S22)*DER*DER DOR= (S11+S12)*DER+DAC DVAP=DVAP+AMV*DAC !+/- ENDIF DLAP=DLAP-AML*DOR !+/- ENDIF C C BMAT CONTRIBUTION C IF(ABS(BML).GT.TINY) THEN S3=P1(N1,M1)*DFP2(N2,JCHOP)-PP1(N1,M1)*DF2(N2,JCHOP)+ X DP1(N1,M1)*FP2(N2,JCHOP)-DPP1(N1,M1)*F2(N2,JCHOP) DOV=-S3/(EPS1(N1,M1)-EPS2(N2)) DLBP=DLBP+BML*DOV DVBP=DVBP+BMV*DOV ENDIF C C CMAT CONTRIBUTION C IF(ABS(CM).GT.TINY) THEN SP=VMAT(-3,N1,N2,M1)*RK(N2,JCHOP) IF(N2.EQ.JCHOP) SP=SP+UMAT(-3,N1,N2,M1) SP=SP*DER*DER*CM DLAP=DLAP+SP DVAP=DVAP+SP ENDIF C C DMAT CONTRIBUTION C IF(ABS(DML1).GT.TINY) THEN SP=VMAT(-2,N1,N2,M1)*RK(N2,JCHOP) IF(N2.EQ.JCHOP) SP=SP+UMAT(-2,N1,N2,M1) DLBP=DLBP+SP*DML1 DVBP=DVBP+SP*DMV1 ENDIF IF(ABS(DML2).GT.TINY) THEN SP=VMAT(-3,N1,N2,M1)*RK(N2,JCHOP) IF(N2.EQ.JCHOP) SP=SP+UMAT(-3,N1,N2,M1) DLBP=DLBP+SP*DML2 DVBP=DVBP+SP*DMV2 ENDIF C ENDDO ENDDO C DLP=DLAP+DLBP DVP=DVAP-TWO*DVBP*DER C RETURN END C********************************************************************* C SUBROUTINE DSQDT(DLVECT,DVVECT,FLAGL,FLAGV,DL2,DV2,QJUMP) C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C CBL PARAMETER (MWORK=100*MZDEG) PARAMETER (TZERO=0.0) PARAMETER (ZERO=(0.0,0.0)) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (BIG=150.0) C LOGICAL FLAGV,QDT,BMQDT,FLAGL,QJUMP,SKIPPY C C NRB: C EVALUATE SINGLE-STATE QDT, DETAILED AND AVERAGED. C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF), X RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAGQ/QDT COMMON/TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) C COMMON/NRBKHI/ZKHICC(MZDEG,MZDEG),ZKHIOC(MZCHF,MZDEG),ZVAL(MZDEG) X ,ZDL(MZCHF),ZDV(MZCHF),ZDL0(MZCHF),ZDV0(MZCHF) CBL X,ZWL(MZDEG,MZDEG),ZWR(MZDEG,MZDEG),ZWORK(MWORK),RWORK(2*MZDEG) COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) C DIMENSION DLVECT(MZCHF),DVVECT(MZCHF),DL2(MZCHF),DV2(MZCHF) DIMENSION XMOLD(MZCHF,MZCHF),YMOLD(MZCHF,MZCHF) C EQUIVALENCE (XMOLD,DFP2),(YMOLD,DFPP2) C CBL DATA LWMAX/MWORK/ C C C ALL CHANNELS C NCHF=NCHND2 C C SEE IF THERE WAS A PRIOR CONTRACTION VIA DMQDT C BMQDT=IQDT.GT.0.AND.(NCHOP.LT.NCHF.OR.QJUMP) C C FIRST DETERMINE NWC (NUMBER OF WEAKLY CLOSED CHANNELS) C DO J=NCHOP,1,-1 IF(EPS2(J).GE.TZERO) THEN NWC=NCHOP-J NCHOP=J GOTO 2 ENDIF ENDDO C C CASE OF NO REAL OPEN CHANNELS C NCHOP=0 RETURN C C CASE ALL CHANNELS OPEN C 2 IF(NWC.EQ.0)RETURN C PI=ACOS(-ONE) TPI=TWO*PI C C DIAGONALIZE KHICC; USE (KHI) = CONJG(2I(Y)-1) = (2(XM)-1,2(YM)) C IF(NWC.GT.MZDEG) THEN WRITE(6,601)NWC CALL ABORTT(10) ENDIF DO N2=1,NWC DO N1=1,NWC ZKHICC(N1,N2)=DCMPLX(TWO*XM(NCHOP+N1,NCHOP+N2), X TWO*YM(NCHOP+N1,NCHOP+N2)) ENDDO ZKHICC(N2,N2)=ZKHICC(N2,N2)-ONE ENDDO C CSTRTNBL CALL ZEIGEN(ZKHICC,ZVAL,NWC,AC) CENDNBL C CSTRTBL CBL CALL ZGEEV('N','V',NWC,ZCHICC,MZDEG,ZVAL,ZWL,MZDEG,ZWR,MZDEG, CBL X ZWORK,MWORK,RWORK,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,'("ZGEEV CALLED, INFO=",I5)') INFO CBL STOP CBL ENDIF CBL LWOPT=INT(ZWORK(1)) CBL IF (LWOPT.GT.LWMAX) THEN CBL LWMAX=LWOPT CBL WRITE(6,'("ZGEEV: OPTIMAL WORK SPACE LENGTH =",I5)') LWOPT CBL ENDIF CBL DO N1=1,NWC CBL DO N2=1,NWC CBL ZKHICC(N2,N1)=ZWR(N2,N1) CBL ENDDO CBL ENDDO CENDBL C C AFTER THE LAST CALL, THE DIAGONALIZATION MATRIX (C)=(ZKHICC). C NOW TRANSFORM ZKHIOC C DO N2=1,NWC DO N1=1,NCHOP ZKHIOC(N1,N2)=ZERO ENDDO DO K=1,NWC DO N1=1,NCHOP ZKHIOC(N1,N2)=ZKHIOC(N1,N2)+DCMPLX(TWO*XM(N1,NCHOP+K), X TWO*YM(N1,NCHOP+K))*ZKHICC(K,N2) ENDDO ENDDO ENDDO C C RADIATIVE DECAYS C IF(IRDEC.EQ.1)THEN T=ARDEC(NTARG(NCHOP+1))*(FNU2(NCHOP+1)**3) T=MIN(T,BIG) FDEC=EXP(T) ELSE FDEC=ONE ENDIF C C TRANSFORM DIPOLE MATRIX ELEMENT VIA D(KHI)=Y**(T*)D(R). C IF THERE HAS BEEN A PRIOR CONTRACTION THEN WE ALREADY HAVE D(KHI). C IF(FLAGL)THEN IF(.NOT.BMQDT)THEN DO I=1,NCHOP+NWC ZDL0(I)=ZERO DO K=1,NCHOP+NWC ZDL0(I)=ZDL0(I)+DLVECT(K)*DCMPLX(-YM(K,I),XM(K,I)) ENDDO ENDDO ENDIF DO N=1,NWC ZDL(N)=ZERO DO K=1,NWC ZDL(N)=ZDL(N)+ZDL0(NCHOP+K)*ZKHICC(K,N) ENDDO ENDDO ENDIF IF(FLAGV) THEN IF(.NOT.BMQDT)THEN DO I=1,NCHOP+NWC ZDV0(I)=ZERO DO K=1,NCHOP+NWC ZDV0(I)=ZDV0(I)+DVVECT(K)*DCMPLX(-YM(K,I),XM(K,I)) ENDDO ENDDO ENDIF DO N=1,NWC ZDV(N)=ZERO DO K=1,NWC ZDV(N)=ZDV(N)+ZDV0(NCHOP+K)*ZKHICC(K,N) ENDDO ENDDO ENDIF C C IF(QDT)THEN C C FINALLY COMPLETE THE GAILITIS SUM C IF(FLAGL)THEN DO J=1,NCHOP TZ=TZERO DO K=1,NWC VV=ABS(ZVAL(K)) IF((ONE-VV).LT.AC)THEN VV=FDEC-ONE DO M=1,NCHOP VV=VV+ABS(ZKHIOC(M,K))**2 ENDDO ELSE VV=FDEC-VV**2 ENDIF TZ=TZ+(ABS(ZKHIOC(J,K))*ABS(ZDL(K)))**2/VV ENDDO DL2(J)=DL2(J)+TZ Z=ZERO DO K1=1,NWC-1 DO K2=K1+1,NWC Z=Z+ZKHIOC(J,K1)*ZDL(K1)*CONJG(ZKHIOC(J,K2)* X ZDL(K2))/(FDEC-ZVAL(K1)*CONJG(ZVAL(K2))) ENDDO ENDDO DL2(J)=DL2(J)+TWO*DBLE(Z) ENDDO ENDIF IF(FLAGV) THEN DO J=1,NCHOP TZ=TZERO DO K=1,NWC VV=ABS(ZVAL(K)) IF((ONE-VV).LT.AC)THEN VV=FDEC-ONE DO M=1,NCHOP VV=VV+ABS(ZKHIOC(M,K))**2 ENDDO ELSE VV=FDEC-VV**2 ENDIF TZ=TZ+(ABS(ZKHIOC(J,K))*ABS(ZDV(K)))**2/VV ENDDO DV2(J)=DV2(J)+TZ Z=ZERO DO K1=1,NWC-1 DO K2=K1+1,NWC Z=Z+ZKHIOC(J,K1)*ZDV(K1)*CONJG(ZKHIOC(J,K2)* X ZDV(K2))/(FDEC-ZVAL(K1)*CONJG(ZVAL(K2))) ENDDO ENDDO DV2(J)=DV2(J)+TWO*DBLE(Z) ENDDO ENDIF C ELSE C C DETAILED C TPINU=FNU2(NCHOP+1)*TPI ZFDEC=SQRT(FDEC)*EXP(DCMPLX(TZERO,-TPINU)) IF(FLAGL)THEN DO J=1,NCHOP Z=ZDL0(J) DO K=1,NWC Z=Z-ZKHIOC(J,K)*ZDL(K)/(ZVAL(K)-ZFDEC) ENDDO DL2(J)=Z*CONJG(Z) ENDDO ENDIF IF(FLAGV)THEN DO J=1,NCHOP Z=ZDV0(J) DO K=1,NWC Z=Z-ZKHIOC(J,K)*ZDV(K)/(ZVAL(K)-ZFDEC) ENDDO DV2(J)=Z*CONJG(Z) ENDDO ENDIF C ENDIF C C IF REALLY MQDT RESTORE OLD XM,YM C IF(BMQDT)THEN DO J=1,NCHF DO I=1,NCHF XM(I,J)=XMOLD(I,J) YM(I,J)=YMOLD(I,J) ENDDO ENDDO ENDIF C RETURN C 601 FORMAT(//' ***NUMBER OF DEGENERATE WC CHANNELS, NWC = ',I2/4X, 1 'LARGER THAN DIMENSION VALUE OF DEG = MZDEG'//) END C************************************************* SUBROUTINE DSTOR(KEY,DLR,DLI,DVR,DVI,NCHOP) C C STORES/RETRIEVES DIPOLE VECTORS USING BUFFER AND DIRECT ACCESS FILE. C BUFFER STORES VECTORS FOR ALL INITIAL ENERGIES (M1); C FILE STORES BUFFER FOR ALL FINAL ENERGIES (M2) AND SYMMETRIES (IL2) C BUT IS OVERWRITTEN WHEN INITIAL SYMMETRY (IS1,IL1,IP1) CHANGES. C C KEY = 1 FOR RETRIEVE, = 2 FOR STORE, C DLR,DLI,DVR,DVI,NCHOP ARE THE DIPOLE VECTORS OF SIZE NCHOP C (LENGTH REAL,LENGTH IMAG,VELOCITY REAL,VELOCITY IMAG) C AND ARE INPUT IF KEY=2, OR RETURNED IF KEY=1. C ALL PARAMETERS IN /CSTOR/ MUST BE DEFINED ON INPUT IF KEY=2 C (MFIN,FLAGV ARE NOT REQUIRED IF KEY=1). C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C LOGICAL FLAGV,FLAGL PARAMETER (MXBUF=4*MZCHF*MZEST,IFAC=50) DIMENSION DLR(MZCHF),DLI(MZCHF),DVR(MZCHF),DVI(MZCHF) COMMON/CSTOR/FAC(2,IFAC),ITOR,IS1,IL1,IP1,M11,MFIN,M1,IL2,M2,FLAGV X ,FLAGL DIMENSION BUFFER(MXBUF) DIMENSION NLEN(MZSA3,-1:1),NREC(MZSA3,-1:1),NCH(MZSA3,-1:1) SAVE BUFFER,NLEN,NREC,NCH,IREC, JS1,JL1,JP1,IR C DATA IREC/-1/, JS1,JL1,JP1/3*-1/, JDISC/99/ C IF(KEY.NE.1.AND.KEY.NE.2) RETURN IDEL=IL2-IL1 IF(IS1.EQ.0)IDEL=IDEL/2 C C --- INITIALIZE IF NEW INITIAL SYMMETRY (IS1,IL1,IP1) C IF(IS1.NE.JS1.OR.IL1.NE.JL1.OR.IP1.NE.JP1) THEN IF(KEY.EQ.1) THEN PRINT *,' *** ERROR IN DSTOR *** READING BEFORE WRITING ', * IS1,IL1,IP1,JS1,JL1,JP1 STOP ENDIF JS1=IS1 JL1=IL1 JP1=IP1 DO 2 L=-1,1 DO 2 I=1,MZSA3 2 NCH(I,L)=0 IF(IREC.NE.-1) IREC=1 IR=2 IF(FLAGV) IR=4 ENDIF C C --- DEFINE STORE/RETRIEVE POINTERS C IT IS ASSUMED THAT NCHOP IS CONSTANT FOR ALL INITIAL ENERGIES (M1) C IF(KEY.EQ.1) NCHOP=NCH(M2,IDEL) IF(NCHOP.EQ.0) RETURN C WE KL=IR*(M11-M1)*NCHOP -- CORR'93MAR7 MARSEILLE: KL = (M1-M11)*NCHOP*IR KV=KL+2*NCHOP IF(KEY.EQ.2) THEN C C --- STORE C IF(NCH(M2,IDEL).EQ.0) NCH(M2,IDEL)=NCHOP IF(NCH(M2,IDEL).NE.NCHOP) THEN PRINT *,' *** ERROR IN DSTOR *** CHANNELS NOT CONSISTENT ', * IS1,IL1,IP1,IL2,M2,NCH(M2,IDEL),NCHOP STOP ENDIF DO 12 I=1,NCHOP BUFFER(KL+I)=DLR(I) 12 BUFFER(KL+I+NCHOP)=DLI(I) IF(IR.EQ.4) THEN DO 14 I=1,NCHOP BUFFER(KV+I)=DVR(I) 14 BUFFER(KV+I+NCHOP)=DVI(I) ENDIF IF(M1.EQ.MFIN) THEN LENGTH=IR*(MFIN-M11+1)*NCHOP NLEN(M2,IDEL)=LENGTH NREC(M2,IDEL)=MAX(IREC,1) CALL DA2(KEY,IREC,JDISC,LENGTH,BUFFER) ENDIF ELSE C C --- RETRIEVE C LENGTH=NLEN(M2,IDEL) IF(KL+IR*NCHOP.GT.LENGTH) THEN NCHOP=0 RETURN ENDIF IREC=NREC(M2,IDEL) CALL DA2(KEY,IREC,JDISC,LENGTH,BUFFER) DO 22 I=1,NCHOP DLR(I)=BUFFER(KL+I) 22 DLI(I)=BUFFER(KL+I+NCHOP) IF(IR.EQ.4) THEN DO 24 I=1,NCHOP DVR(I)=BUFFER(KV+I) 24 DVI(I)=BUFFER(KV+I+NCHOP) ENDIF ENDIF RETURN END C C*********************************************************************** C SUBROUTINE FLAG1(M11,MLAST) C C INITIATE FLAGS FOR INITIAL CHANNELS TO INDICATE DROP FROM RTWO TO C INFINITY C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/CACC/AC0,ACNUM,ACJWBK,ACZP,LACC COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAGIN/FLG1(MZCHF,MZEST),FLG2(MZCHF),EXIST1(MZCHF,MZEST), 1 EXIST2(MZCHF),DROP COMMON/NRBZED/TZED C LOGICAL FLG1,FLG2,DROP,EXIST1,EXIST2 C C...MJS MOD 27.3.87 DROP=.TRUE. DO 20 M1=M11,MLAST IF(KPMM(M1).LT.KP2)THEN DO 10 N=1,NCHND1 EP1=EE1(M1)-ECH1(N) EPS1(N,M1)=EP1 FNU1(N,M1)=SQRT(-1./EP1) EXIST1(N,M1)=.FALSE. 10 FLG1(N,M1)=.FALSE. ELSE C...ENDMOD C DO 1 N=1,NCHND1 EP1=EE1(M1)-ECH1(N) EPS1(N,M1)=EP1 FN1=SQRT(-ONE/EP1) FNU1(N,M1)=FN1 EXIST1(N,M1)=.FALSE. EC=EP1*CC1(N) FLG1(N,M1)=.FALSE. IF(TZED.EQ.TZERO)THEN RI=RZERO ELSE IF(EC.LE.-ONE) GO TO 1 RI=-(SQRT(ONE+EC)+ONE)/EP1 ENDIF IF(RI.GT.RTWO) THEN FLG1(N,M1)=.TRUE. DROP=.FALSE. IF(IPRINT.GT.0) WRITE(6,600)N,RI,RTWO ELSE KI=1 IF(RI.GT.RZERO) KI=1+(RI-RZERO)/H MA=1 IF(MZSA1 .NE. 1) MA=M1 FLG1(N,M1) = ABS(FSD1(KP2,N,MA)/FSD1(KI,N,MA)).GE.AC IF(FLG1(N,M1)) DROP=.FALSE. ENDIF 1 CONTINUE END IF 20 CONTINUE C RETURN 600 FORMAT(/' *** TRANSITION DATA FOR THE NEXT INITIAL STATE MAY', * ' BE INACCURATE ***'/10X,'CHANNEL',I3, * ' (RINF =',F9.2,').GT.(RTWO =',F9.2,')'/) END C****************************************************************** C SUBROUTINE FLAG2(E2) C C INITIATE FLAGS FOR FINAL CHANNELS TO INDICATE DROP FROM RTWO TO C INFINITY C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C LOGICAL FLG1,FLG2,DROP,EXIST1,EXIST2,EXISTV,QDT,SKIPPY C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAGIN/FLG1(MZCHF,MZEST),FLG2(MZCHF),EXIST1(MZCHF,MZEST), 1 EXIST2(MZCHF),DROP COMMON /VUSTOR/ VSTORE(-3:1,MZCHF,MZCHF),USTORE(-3:1,MZCHF,MZCHF), * X1(MZCHF,MZEST),X2(MZCHF),EXISTV(-3:1,MZCHF,MZCHF) COMMON/FLAGQ/QDT C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBZED/TZED C DATA TNY/-1.E-6/ C C DO N=1,NCHND2 EP2=E2-ECH2(N) IF(EP2.GT.TNY.AND.EP2.LT.0.)EP2=TZERO C NRB: SURELY ALL WE NEED IS THIS C IF(EP2.LT.TZERO)THEN C RATHER THAN IF(N.GT.NCHOP.OR.((QDT.OR.IQDT.NE.0).AND.EP2.LT.TZERO)) THEN C SINCE WE ALWAYS WANT EFFECTIVE N WHEN EP2 IS NEGATIVE...... C FN=SQRT(-ONE/EP2) ELSE FN=-ONE ENDIF EPS2(N)=EP2 FNU2(N)=FN EXIST2(N)=.FALSE. ENDDO C IF(DROP) THEN DO N=1,NCHND2 FLG2(N)=.FALSE. ENDDO RETURN ENDIF C DO N=1,NCHND2 EP2=EPS2(N) IF(N.LE.NCHOP) THEN C C WC CHANNELS DROPPED BECAUSE NO ASYMPTOTIC EXPANSION AVAILABLE C NRB: SURELY ALL WE NEED IS THIS C FLG2(N)=EP2.GE.TZERO C RATHER THAN FLG2(N)=EP2.GE.TZERO.OR.(.NOT.QDT.AND.IQDT.EQ.0) C SINCE IN NON-QDT MODE WE CAN'T HAVE NEGATIVE EP2 HERE (N.LE.NCHOP)...... ELSE EC=EP2*CC2(N) FLG2(N)=.FALSE. IF(TZED.EQ.TZERO)THEN RI=RZERO ELSE IF(EC.LE.-ONE) GO TO 1 RI=-(ONE+SQRT(ONE+EC))/EP2 ENDIF IF(RI.GT.RTWO) THEN FLG2(N)=.TRUE. IF(IPRINT.GT.0) WRITE(6,600)N,RI,RTWO ELSE KI=1 IF(RI.GT.RZERO) KI=(RI-RZERO)/H+KI FLG2(N) = ABS(C(KP2,N)/C(KI,N)).GE.AC ENDIF ENDIF 1 ENDDO C RETURN C 600 FORMAT(/' *** TRANSITION DATA FOR THE NEXT ENERGY MAY BE ', * 'INACCURATE ***'/10X,'CHANNEL',I3,' (RINF =',F9.2, * ').GT.(RTWO =',F9.2,')'/) END C C*********************************************************************** C SUBROUTINE GET1(FILE,IS1,IL1,IP1,M11,M12) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C C GET AN INITIAL STATE SLPI FILE AND RANGE OF BOUND LEVELS C COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT COMMON/LIST1/KSLP1,MSLP1(100) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/SLPI2/IS2V(3),IL2V(3),IP2V(3),KUT CHARACTER FILE*3,NUM(0:9) LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ SAVE KASE C IF(ALL1) THEN IF(RED) THEN KASE=1 RED=.FALSE. ELSE KASE=KASE+1 ENDIF IF(KASE.EQ.KSLP1) MORE1=.FALSE. IS1=MSLP1(KASE)/10000 IL1=(MSLP1(KASE)-IS1*10000)/100 IP1=MSLP1(KASE)-(IS1*100+IL1)*100 MF=KASE ELSE CALL RD5L1(IS1,IL1,IP1,M11,M12) IP1 = MOD(IP1,2) CALL RD5L2 IF(EOF) THEN MORE1=.FALSE. ELSE CALL RD5L1(IQ,JUNK1,JUNK2,JUNK3,JUNK4) IF(EOF.OR.JUNK1.EQ.-1) THEN MORE1=.FALSE. ELSE BACKSPACE 5 ENDIF ENDIF ISLP1=(IS1*100+IL1)*100+IP1 DO 1 K=1,KSLP1 IF(ISLP1.NE.MSLP1(K)) GO TO 1 MF=K GO TO 2 1 CONTINUE WRITE(6,500)IS1,IL1,IP1 OK=.FALSE. RETURN ENDIF C 2 FILE='B'//NUM(MF/10)//NUM(MF-10*(MF/10)) IF(IPRINT.NE.0) WRITE(6,501)FILE C C CMPLT IS CLEARED IF TOTAL X-SECTION SUM IS UNSAFE, WILL WARN AT END CMPLT = ALL2 .OR. IL1.EQ.0 .OR. KUT.EQ.3 OK=.TRUE. RETURN C 500 FORMAT(//' *** DATA FOR ',3I3,' NOT ON B DATASET'//) 501 FORMAT(' OPENING FILE ',A3) END C C*************************************************************** C SUBROUTINE GET2(FILE,IS1,IL1,IP1,IS2,IL2,IP2) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C C GET A FINAL STATE SLPI FILE C COMMON/SLPI2/IS2V(3),IL2V(3),IP2V(3),KUT COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/LIST2/KSLP2,MSLP2(100) CHARACTER FILE*3, NUM(0:9) LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ SAVE KRRNT C IF(ALL2) THEN IF(RED) THEN IL2=IL1+1 IF(IS1.EQ.0) IL2=IL2+1 IS2=IS1 IP2=1-IP1+2*(IP1/2) RED=.FALSE. ELSE IL2=IL2-1 IF(IS1.EQ.0) IL2=IL2-1 ENDIF IF(IL2.EQ.0.OR.IL2.EQ.IL1-1) MORE2=.FALSE. IF(IL1.EQ.0) MORE2=.FALSE. IF(IS1.EQ.0.AND.(IL2.EQ.IL1-2.OR.IL2.EQ.1)) MORE2=.FALSE. ELSE IF(RED) THEN KRRNT=1 RED=.FALSE. ELSE KRRNT=KRRNT+1 ENDIF IS2=IS2V(KRRNT) IL2=IL2V(KRRNT) IP2=IP2V(KRRNT) IF(IS2.NE.IS1) THEN WRITE(6,505)IS1,IS2 CALL ABORTT(2) ENDIF IF(IL2.GT.IL1+2 .OR.IL2.LT.IL1-2 .OR.IL2.LT.0) THEN WRITE(6,506)IL1,IL2 CALL ABORTT(3) ENDIF IF(IP2+2*(IP2/2).NE.1-IP1+2*(IP1/2)) THEN WRITE(6,507)IP1,IP2 CALL ABORTT(4) ENDIF IF(KRRNT.EQ.KUT) MORE2=.FALSE. ENDIF C ISLP2=(IS2*100+IL2)*100+IP2 DO 3 K=1,KSLP2 IF(ISLP2.NE.MSLP2(K)) GO TO 3 FILE='F'//NUM(K/10)//NUM(K-10*(K/10)) IF(IPRINT.NE.0) WRITE(6,501)FILE OK=.TRUE. GO TO 4 3 CONTINUE C C DATA NOT FOUND. CLEAR CMPLT TO GIVE WARNING AT THE END WRITE(6,500)IS2,IL2,IP2 OK=.FALSE. CMPLT=.FALSE. 4 RETURN C C 500 FORMAT(//' *** DATA FOR ',3I3,' NOT ON F DATASET'//) 501 FORMAT(' OPENING FILE ',A3) 505 FORMAT(//10X,'IS1 =',I4,10X,'IS2 =',I4) 506 FORMAT(//10X,'IL1 =',I4,10X,'IS2 =',I4) 507 FORMAT(//10X,'IP1 =',I4,10X,'IP2 =',I4) END C C********************************************************************* C SUBROUTINE GET3DV(IS1,IL1,IP1,IS2,IL2,IP2) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C C GET SLPI COMBINATION FILE C C MODIFICATION TO FIND KKP PARAMETER (MXTRAN=99) COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/LIST0/KOUNT,MSLPL(MXTRAN),MSLPR(MXTRAN) LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT C ISLP1=(IS1*100+IL1)*100+IP1 ISLP2=(IS2*100+IL2)*100+IP2 DO 1 K=1,KOUNT IF(ISLP1.NE.MSLPL(K) .OR.ISLP2.NE.MSLPR(K)) GO TO 1 MF=K KKP=.TRUE. GOTO 3 1 CONTINUE DO 2 K=1,KOUNT IF(ISLP1.NE.MSLPR(K) .OR.ISLP2.NE.MSLPL(K)) GO TO 2 MF=K KKP=.FALSE. GOTO 3 2 CONTINUE C WRITE(6,500)IS1,IL1,IP1,IS2,IL2,IP2 OK=.FALSE. WRITE(6,502) C 3 RETURN C 500 FORMAT(/'0*** DATA FOR ',3I3,3X,3I3, * ' COMBINATION NOT FOUND ON D DATASET'//) 501 FORMAT(' OPENING FILE ',A4) 502 FORMAT(' ***WARNING*** TOTAL CROSS SECTION SUMMATION MAY BE'/ 1 16X,'INCOMPLETE FOR THIS INITIAL STATE.'//) END C*************************************************************** C SUBROUTINE INJWBK(E,L) C C COMPUTES ARRAY D WHICH IS USED FOR CALCULATION OF JWBK FUNCTIONS. CNRB C NEUTRAL CASE ADDED C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMON/CJWBK/D(15) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/NRBZED/TZED C D(1)=E C IF(TZED.EQ.0)THEN FK=SQRT(E) D(2)=FK IF(L.GT.0)THEN D(3)=1./FK C=DBLE(L*(L+1)) D(4)=C SC=SQRT(C) D(5)=SC D(6)=(C+.125)/SC D(7)=E*C D(12)=6.*E*C D(14)=-C*C D(15)=-L*1.5707963 ELSE D(4)=0.0 D(15)=0.0 ENDIF RETURN ENDIF C IF(L.GT.0)GOTO 10 C C CASE OF L.EQ.0 D(4)=0. IF(E.EQ.0)GOTO 30 C CASE OF L.EQ.0 AND E.GT.0 FK=SQRT(E) D(2)=FK D(3)=1./FK GOTO 30 C 10 IF(E.GT.0)GOTO 20 C C CASE OF L.GT.0 AND E.EQ.0 C=DBLE(L*(L+1)) D(4)=C SC=SQRT(C) D(5)=SC D(6)=(C+.125)/SC D(13)=6.*C D(14)=-C*C GOTO 30 C C CASE OF L.GT.0 AND E.GT.0 20 FK=SQRT(E) D(2)=FK D(3)=1./FK C=DBLE(L*(L+1)) D(4)=C SC=SQRT(C) D(5)=SC D(6)=(C+.125)/SC A=1.+E*C D(7)=A A=3.*A D(8)=A-1. D(9)=A+1. D(10)=FK*C D(11)=-4.*E D(12)=-9.+2.*A D(13)=6.*C D(14)=-C*C C C TERM IN ARG GAMMA ETC 30 D(15)=ARGC(E,L,AC) C RETURN END C C*************************************************************** C SUBROUTINE MATMUL(A,B,C,N) C C CALCULATES C=A*B C C NRB: VECTOR OR SCALAR VERSIONS C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (ZERO=0.0) C DIMENSION A(MZCHF,MZCHF),B(MZCHF,MZCHF),C(MZCHF,MZCHF) 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 MATSQ(A,B,N) C C CALCULATES B=A*A C C NRB: VECTOR OR SCALAR VERSIONS C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (ZERO=0.0) C DIMENSION A(MZCHF,MZCHF),B(MZCHF,MZCHF) CCRAY CRAY DO J=1,N CRAY DO I=1,N CRAY B(I,J)=ZERO CRAY ENDDO CRAY ENDDO CRAY DO K=1,N CRAY DO J=1,N CRAY DO I=1,N CRAY B(I,J)=B(I,J)+A(I,K)*A(K,J) CRAY ENDDO CRAY ENDDO CRAY ENDDO CCRAY C---- DO J=1,N DO I=1,N B(I,J)=ZERO ENDDO DO K=1,N DO I=1,N B(I,J)=B(I,J)+A(I,K)*A(K,J) ENDDO ENDDO ENDDO C---- RETURN END C********************************************************************** C SUBROUTINE MESH(E0,EINCR,MXE,EMSH0,MXE0,EMESH,MXE2,EPS) C C NRB: C GENERATE A NEW ENERGY MESH FOR THE FINAL STATE I.E. ELECTRON CONTINUUM, C GIVEN THE MESH AT WHICH THE UNPHYSICAL S-MATRIX HAS BEEN TABULATED BY C STGF. IT INCLUDES THE ORIGINAL ENERGIES ("IMESH=1"), BUT NEED NOT. C THIS IS RELATIVELY SIMPLE AT THE MOMENT, BUT CAN BE BUILT UPON. C FOR NOW, THE NEW MESH IS GENERATED USING E0, EINCR, MXE AS IMESH=1 OF C STGF AND INPUT VIA THE STGBF NAMELIST. C ONE COULD ENVISAGE A DQN, QNMAX ETC AS WELL AS PER IMESH=2. C NO EXTRAPOLATION OUTSIDE OF THE STGF MESH IS ALLOWED, THE REQUESTED C RANGE IS CURTAILED IF NECESSARY. C C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER(TZERO=0.0) C DIMENSION EMESH(MZMSH),EMSH0(MZMSH0) C C FIRST TRANSFER INPUT MESH FROM STGF TO EMSH0. THIS SHOULD BE COARSE, C TYPICALLY A FEW HUNDRED ENERGIES, BUT LET'S CHECK SPACE FIRST: C MXE0=MXE2 IF(MXE0.GT.MZMSH0)THEN WRITE(6,*)' INPUT ENERGY MESH TRUNCATED, MAX ENERGY IS ' X ,EMESH(MZMSH0) IF(MXE0.GT.5000)THEN WRITE(6,*)' DO YOU REALLY NEED ALL THESE ENERGIES? WE ARE' X ,' USING MQDT AFTER ALL!' ELSE WRITE(6,*)' INCREASE MZMSH0 TO AT LEAST',MXE0 ENDIF MXE0=MZMSH0 ENDIF C C NOW TRANSFER C DO I=1,MXE0 EMSH0(I)=EMESH(I) ENDDO C C DO A QUICK CHECK TO SEE IF A NEW MESH HAS BEEN SET, IF NOT THEN C JUST USE THE STGF MESH. C IF(E0.LT.TZERO.OR.MXE.LE.0.OR.EINCR.LE.TZERO)RETURN C C NOW OVERWRITE EMESH WITH NEW MESH AND REDEFINE MXE2 C EPS=EINCR*0.1 EMAX=E0+(MXE-1)*EINCR C C FIRST FIND CORRESPONDING STARTING POINT IN EMSH0 C IF(E0.LT.EMSH0(1)-EPS)THEN WRITE(6,*)' ***WARNING: STARTING ENERGY INCREASED TO',EMSH0(1) X ,' RE-RUN STGF WITH LOWER ENERGY IF INSUFFICIENT' E0=EMSH0(1) IF(E0.GT.EMAX+EPS)STOP 'INVALID ENERGY RANGE REQUESTED' IT=2 ELSE DO I=MXE0,1,-1 IF(E0.GT.EMSH0(I)-EPS)THEN E0=EMSH0(I) EMESH(1)=E0 IT=I+1 GO TO 5 ENDIF ENDDO ENDIF C C NOW LOOK AT FINISHING POINT C 5 IF(EMAX.GT.EMSH0(MXE0)+EPS)THEN WRITE(6,*)' ***WARNING: FINAL ENERGY DECREASED TO',EMSH0(MXE0) X ,' RE-RUN STGF WITH HIGHER ENERGY IF INSUFFICIENT' EMAX=EMSH0(MXE0) ELSE MX00=MXE0 DO I=1,MX00 IF(EMAX.LT.EMSH0(I)+EPS)THEN EMAX=EMSH0(I) MXE0=I GO TO 7 ENDIF ENDDO ENDIF C 7 I2=1 10 IF(I2.LT.MZMSH)THEN T=EMESH(1)+I2*EINCR I2=I2+1 IF(T+EPS.GE.EMSH0(IT))THEN T=EMSH0(IT) IT=IT+1 ENDIF EMESH(I2)=T IF(T.LT.EMAX-EPS)GO TO 10 GO TO 20 ELSE WRITE(6,*)' INPUT ENERGY MESH TRUNCATED, MAX ENERGY IS ' X ,EMESH(I) WRITE(6,*)' INCREASE MZMSH, OR RE-RUN STGF OVER A NARROWER' X ,' ENERGY RANGE' ENDIF C 20 MXE2=I2 C RETURN END C************************************************************************ C SUBROUTINE NIGEL(SAVE,IPERT0,QJUMP,NCHOLD,NOQDT,MLAST,FIRST X ,SIG1,OMEG1,EMESH,MXE2,EJUMP) C C INNER LOOP OVER INITIAL AND FINAL ENERGIES FOR FIXED SYMMETRIES. C CALLED TWICE IN IQDT>0 MODE. FIRSTLY ON THE COARSE ENERGY MESH C (EMSH0) TO GENERATE UNPHYSICAL S- AND D-MATRICES, AND WRITE THEM TO C SCRATCH FILE. SECONDLY ON THE FINE ENERGY MESH (EMESH2), TO READ AND C INTERPOLATE THE UNPHYSICAL MATRICES, TO GENERATE PHYSICAL DATA. C N.B. MXE0 SHOULD BE SMALL, ONLY ORDER 100 ENERGIES SO LIMITED I/O. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (IFAC=50) PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,FLAGV,FLAGP, 1 QDT,NOQDT,CMPLT,FIRST,EXISTV,SAVE,QJUMP,FLAGL,EJUMP,SKIPPY C COMMON/CACC/AC0,ACNUM,ACJWBK,ACZP,LACC COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/SCALE/AZ,AZP,AZ2,AZP2,K10,K20 COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/SLPI2/IS2V(3),IL2V(3),IP2V(3),KUT COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT COMMON/FLAGQ/QDT COMMON /TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) COMMON /VUSTOR/ VSTORE(-3:1,MZCHF,MZCHF),USTORE(-3:1,MZCHF,MZCHF), * X1(MZCHF,MZEST),X2(MZCHF),EXISTV(-3:1,MZCHF,MZCHF) COMMON/CSTOR/FAC(2,IFAC),ITOR,IS1,IL1,IP1,M11,MFIN,M1,IL2,M2,FLAGV X ,FLAGL C COMMON/NRBDR/OMEGPR(MZMET,MZMSH),NDRMET,IGAUGE COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBTMP/E2 C COMMON/NRBDDD/DLVECT(MZCHF,MZEST),DVVECT(MZCHF,MZEST) X ,DLT(MZEST),DL2(MZCHF),DV2(MZCHF) X ,TOT(MZSA2,MZEST),DEE(MZEST),DEER(MZEST) X ,DLI(MZCHF,MZEST),DVI(MZCHF,MZEST) C DIMENSION EMESH(MZMSH) C SAVE NCHOPS C C MS=M2 IF (SAVE) MS=1 IPERT=IPERT0 E2=EMESH(M2) C C READ E2-DEPENDENT DATA; SKIP CASE OF NO OPEN CHANNEL C IF(.NOT.QJUMP)CALL R2EDEP(M2,E2,EJUMP) IF(EJUMP)RETURN IF(.NOT.OK) GOTO 99 C C PRESERVE THE VALUE OF NCHOP AS I/QDT MODE CHANGES NCHOP TO C THE NUMBER OF TRUE OPEN CHANNELS. C IF(.NOT.QJUMP)NCHOPS=NCHOP NCHOP=NCHOPS C C C DETECT CHANGE IN OPEN/WC CHANNEL STRUCTURE C IF(IPRINT.GT.0)THEN FLAGV=.TRUE. IF(NCHOP.NE.NCHOLD)THEN IF(.NOT.QDT) THEN NOQDT=.TRUE. IF(IQDT.LE.0)THEN WRITE(6,621) NCHOP NCHOLD=NCHOP ENDIF ELSE IF(NOQDT) THEN WRITE(6,622) NOQDT=.FALSE. ENDIF ENDIF ENDIF C CALL FLAG2(E2) FLAGP = IPERT.NE.0 C C SKIP FINAL ENERGY BELOW INITIAL BOUND STATE C MFIN=MLAST DO M1=M11,MLAST E1=EE1(M1) IF(E1.GT.E2)THEN MFIN=M1-1 GO TO 131 END IF DE=E2-E1 DER=ONE/DE DEE(M1)=DE DEER(M1)=DER ENDDO C 131 IF(IQDT.GT.0.AND..NOT.QJUMP)WRITE(70)M11,MFIN IF(MFIN.LT.M11)GO TO 99 C C CALCULATE INNER REGION CONTRIBUTION, DLI,DVI C IF(.NOT.QJUMP)CALL DIN(E2,NCHOP,M11,MFIN,DLI,DVI) C C CALCULATE OUTER REGION CONTRIBUTION C C IF (NON-RESOLVED) PRODUCTION RUN USE OPTIMISED LOOP C IF (IPRINT.EQ.-1.AND.IPERT.EQ.0) THEN C IF(.NOT.QJUMP)CALL OPDOUT(DEER,M11,MFIN,DLI,DLVECT) CALL OPSTRA(DEE,M11,MFIN,DLVECT,DL2,DLT,QJUMP,FLAGL) C DO M1=M11,MFIN TOT(MS,M1)=TOT(MS,M1)+SIG1*DLT(M1) ENDDO GO TO 136 END IF C C OTHERWISE C DO 132 M1=M11,MFIN C NCHOP=NCHOPS IF(QJUMP)GO TO 155 IF(MZSA1 .EQ. 1) * READ(4,REC=M1)((FSD1(K,N1,1),K=1,KPMM(M1)),N1=1,NCHND1) DER=DEER(M1) IF(KKP)DER=-DER C C CLEAR FLAGS FOR EXISTANCE OF V INTEGRALS C DO I=-3,1 IF(I.NE.-1) THEN DO N1=1,NCHND1 DO N2=1,NCHND2 EXISTV(I,N1,N2)=.FALSE. ENDDO ENDDO ENDIF ENDDO C DO JCHOP=1,NCHOP IF(IOMIT(JCHOP).EQ.0)THEN C C CALCULATE OUTER REGION CONTRIBUTION, DLA,DLB,DLP; DVA,DVB,DVP C CALL DOUT(DLA,DLB,DVA,DVB,DLP,DVP,M1,JCHOP,FLAGP,FLAGV,DER) C C ADD CONTRIBUTIONS FROM ALL TERMS C IF(FLAGL)DLVECT(JCHOP,M1)=DLI(JCHOP,M1)+DLA+DLB+DLP IF(FLAGV) THEN DVIB=-TWO*DER*(DVI(JCHOP,M1)+DVB) DVVECT(JCHOP,M1)=DVIB+DVA+DVP ENDIF C ENDIF ENDDO C C CHANGE TO S-MATRIX NORMALIZATION AND COMPUTE DIPOLE MATRIX SQUARED C 155 CALL STRANS(DLVECT(1,M1),DVVECT(1,M1),FLAGL,FLAGV,DLTOT,DVTOT X ,DL2,DV2,QJUMP,M1) C C SKIP IF JUST WRITING UNPHYSICAL DATA C IF(IQDT.GT.0.AND..NOT.QJUMP)GO TO 132 C C ELSE CALCULATE PHYSICAL DATA C DE=DEE(M1) C SIG=DE*SIG1 OMEG=DE**3*OMEG1 IF(IPRINT.GT.-2) THEN TOT(MS,M1)=TOT(MS,M1)+SIG*DLTOT C C PRINTING RESULTS C IF(IPRINT.GT.-1) THEN IF(IQDT.GT.0.AND.IPRINT.GT.0.AND.NCHOP.NE.NCHOLD)THEN WRITE(6,621) NCHOP NCHOLD=NCHOP ENDIF IF(IPRINT.GE.0)WRITE(6,615)M1,AZ2*DE,AZ2*E2,SIG*DLTOT X ,SIG*DVTOT C C EXTRA PRINTING FOR ... IF(IPRINT.GT.1.AND.NCHOP.GT.1.AND..NOT.QDT) THEN IF(FLAGV.AND.FLAGL) THEN WRITE(6,616) (J,SIG*DL2(J),SIG*DV2(J),J=1,NCHOP) ENDIF IF(FLAGL.AND..NOT.FLAGV)THEN WRITE(6,617) (J,SIG*DL2(J),J=1,NCHOP) ENDIF IF(FLAGV.AND..NOT.FLAGL)THEN WRITE(6,617) (J,SIG*DV2(J),J=1,NCHOP) ENDIF ENDIF ENDIF C C ELSE C C CREATE OR UPDATE A SCRATCH RECORD C CALL WSCRA(IPRINT,M1,M2,FIRST,SIG,OMEG,DLTOT,DL2,DV2 X ,FLAGL,FLAGV,MXE2) C ENDIF C 132 CONTINUE C C C RETURN IF JUST WRITING UNPHYSICAL DATA C 136 IF(IQDT.GT.0.AND..NOT.QJUMP)RETURN C C ELSE CALCULATE PHYSICAL DATA C IF (SAVE) THEN C CALL WSCRA(IPRINT,MLAST,M2,FIRST,SIG,OMEG,TOT,DL2,DV2 X ,FLAGL,FLAGV,MXE2) C DO M1=M11,MLAST TOT(1,M1)=TZERO ENDDO RETURN END IF C C WHEN THE POINT IS SKIPPED FOR ANY BOUND LEVEL WRITE DUMMY LINE C TO SUM FILE WHEN IPRINT=-2 C 99 IF(IPRINT .EQ. -2)THEN IF(FIRST) THEN DO M1=MFIN+1,MLAST CALL WSCRA(-1986,M1,M2,FIRST,SIG,OMEG,TOT,DL2,DV2 X ,FLAGL,FLAGV,MXE2) ENDDO ENDIF ENDIF C RETURN C 615 FORMAT(1X,I3,2(1X,F12.6),1P,2E14.5) 616 FORMAT((49X,I3,1P,2E14.5)) 617 FORMAT((49X,I3,1P,E20.7)) 621 FORMAT(58X,'NCHOP=',I4) 622 FORMAT(' ENTER QDT REGION') END C******************************************************************** C SUBROUTINE NUMB(M1,MA,NCHF11,XVECT,TOL0B) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TZERO=0.0) PARAMETER (PT02=0.02) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (FOUR=4.0) PARAMETER (TEN=10.0) PARAMETER (TWELVE=12.0) C PARAMETER(NCHF1=MZCHF+1) C COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/CHAN/ECH(NCHF1),EPS(NCHF1),FKNU(NCHF1),CC(NCHF1), 1 RINF(NCHF1),LLCH(NCHF1),ITARG(NCHF1),NCHF,NCHOP0,NCHOP1 COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/NRBZED/TZED DIMENSION XVECT(NCHF11) C V(X)=EQ+X*(Q2*TZED-X*CQ) C Q=H*H/TWELVE Q2=Q*TWO KPM=KPMM(M1) E1=EE1(M1) C DO I=1,NCHF11 C EQ=Q*(E1-ECH1(I)) CQ=Q*CC1(I) KP2T=KPM RTWOT=RZERO+H*(KPM-1) IFLAG=0 C IF(ABS(FSD1(KP2T,I,MA))+ABS(FSD1(KP2T-1,I,MA)).LT.TOL0B)THEN IFLAG=1 C C MOVE-IN AND REGENERATE UNTIL NON-VANISHING. NRB 16/06/97 C LLCH(NCHF1)=(SQRT(ONE+FOUR*CC1(I))-ONE+PT02)/TWO FKNU(NCHF1)=SQRT(-ONE/(E1-ECH1(I))) 1 CALL THETA(RTWOT,NCHF1,T,TP,TD,TDP) IF(ABS(T)+ABS(TP).LT.TOL0B)THEN KP2T=KP2T-4 IF(KP2T.LT.1)KP2T=1 RTWOT=(KP2T-1)*H+RZERO DO K=KP2T,KP2T+4 FSD1(K,I,MA)=TZERO ENDDO IF(KP2T.EQ.1)GO TO 210 GO TO 1 ENDIF C1=T C2=C1-H*TP C3=TD C4=C3-H*TDP FSD1(KP2T,I,MA)=C1 FSD1(KP2T-1,I,MA)=C2 ENDIF C R=RTWOT U1=V(ONE/R) R=R-H U2=V(ONE/R) IF(IFLAG.EQ.0)THEN DO K=KP2T-2,1,-1 R=R-H U3=V(ONE/R) FSD1(K,I,MA)=((TWO-TEN*U2)*FSD1(K+1,I,MA)-(ONE+U1)* + FSD1(K+2,I,MA))/(ONE+U3) U1=U2 U2=U3 ENDDO ELSE C C GENERATE THETA AND THETADOT: NRB 16/06/97 C DO K=KP2T-2,1,-1 R=R-H U3=V(ONE/R) T1=(ONE+U1) T2=(TWO-TEN*U2) T3=(ONE+U3) FSD1(K,I,MA)=(T2*FSD1(K+1,I,MA)-T1*FSD1(K+2,I,MA))/T3 C5=(T2*C4-T1*C3-Q*(FSD1(K,I,MA)+TEN*FSD1(K+1,I,MA) X +FSD1(K+2,I,MA)))/T3 C3=C4 C4=C5 U1=U2 U2=U3 ENDDO T=FSD1(1,I,MA) TP=(FSD1(2,I,MA)-FSD1(1,I,MA))/H TD=C4 TDP=(C3-C4)/H AMAX=ONE/MAX(ABS(T),ABS(TP),ABS(TD),ABS(TDP)) T=T*AMAX TP=TP*AMAX TD=TD*AMAX TDP=TDP*AMAX W1=AMAX/SQRT(T*TDP-TP*TD) W1=W1*XVECT(I) DO K=1,KP2T FSD1(K,I,MA)=FSD1(K,I,MA)*W1 ENDDO ENDIF 210 CONTINUE ENDDO C RETURN END C********************************************************** C SUBROUTINE NUMFC(E2,I,INOUT,C1,C2) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TZERO=0.0) PARAMETER (PT02=0.02) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (FOUR=4.0) PARAMETER (TEN=10.0) PARAMETER (TWELVE=12.0) C PARAMETER(NCHF1=MZCHF+1) C LOGICAL SKIPPY C COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CHAN/ECH(NCHF1),EPS(NCHF1),FKNU(NCHF1),CC(NCHF1), 1 RINF(NCHF1),LLCH(NCHF1),ITARG(NCHF1),NCHF,NCHOP0,NCHOP1 COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBZED/TZED C V(X)=EQ+X*(Q2*TZED-X*CQ) C TOL0F=1.D-150 C IF(IOMIT(I).NE.0)THEN DO K=1,KPMAX S(K,I)=TZERO C(K,I)=TZERO ENDDO ENDIF C Q=H*H/TWELVE EQ=Q*(E2-ECH2(I)) CQ=Q*CC2(I) Q2=Q*TWO C IF(INOUT.EQ.0)THEN KP2T=KP2 RTWOT=RTWO IFLAG=0 IF(ABS(C1)+ABS(C2).LT.TOL0F)THEN IFLAG=1 C C MOVE-IN AND REGENERATE UNTIL NON-VANISHING. C NOTE: INPUT C1,C2 ARE ACTUALLY THETA-FUNCTIONS. NRB 14/06/97 C LLCH(NCHF1)=(SQRT(ONE+FOUR*CC2(I))-ONE+PT02)/TWO FKNU(NCHF1)=SQRT(-ONE/(E2-ECH2(I))) C 1 CALL THETA(RTWOT,NCHF1,T,TP,TD,TDP) C IF(ABS(T)+ABS(TP).LT.TOL0F)THEN KP2T=KP2T-4 IF(KP2T.LT.1)KP2T=1 RTWOT=(KP2T-1)*H+RZERO DO K=KP2T,KP2T+4 C(K,I)=TZERO ENDDO IF(KP2T.EQ.1)RETURN GO TO 1 ENDIF C1=T C2=C1-H*TP C3=TD C4=C3-H*TDP ENDIF C C(KP2T,I)=C1 C(KP2T-1,I)=C2 C R=RTWOT U1=V(ONE/R) R=R-H U2=V(ONE/R) IF(IFLAG.EQ.0)THEN DO K=KP2T-2,1,-1 R=R-H U3=V(ONE/R) C(K,I)=((TWO-TEN*U2)*C(K+1,I)-(ONE+U1)*C(K+2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ELSE C GENERATE S AND C: NRB 14/06/97 DO K=KP2T-2,1,-1 R=R-H U3=V(ONE/R) T1=(ONE+U1) T2=(TWO-TEN*U2) T3=(ONE+U3) C(K,I)=(T2*C(K+1,I)-T1*C(K+2,I))/T3 C5=(T2*C4-T1*C3-Q*(C(K,I)+TEN*C(K+1,I)+C(K+2,I)))/T3 C3=C4 C4=C5 U1=U2 U2=U3 ENDDO T=C(1,I) TP=(C(2,I)-C(1,I))/H TD=C4 TDP=(C3-C4)/H AMAX=ONE/MAX(ABS(T),ABS(TP),ABS(TD),ABS(TDP)) T=T*AMAX TP=TP*AMAX TD=TD*AMAX TDP=TDP*AMAX W1=AMAX/SQRT(T*TDP-TP*TD) DO K=1,KP2T C(K,I)=C(K,I)*W1 ENDDO ENDIF ELSE C(1,I)=C1 C(2,I)=C2 R=RZERO U1=V(ONE/R) R=R+H U2=V(ONE/R) DO K=3,KPMAX R=R+H U3=V(ONE/R) C(K,I)=((TWO-TEN*U2)*C(K-1,I)-(ONE+U1)*C(K-2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ENDIF C RETURN END C********************************************************** C SUBROUTINE NUMFO(E2,I,INOUT,S1,S2,C1,C2) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (TEN=10.0) PARAMETER (TWELVE=12.0) C LOGICAL SKIPPY C COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBZED/TZED C V(X)=EQ+X*(Q*TZED-X*CQ) C IF(IOMIT(I).NE.0)THEN DO K=1,KPMAX S(K,I)=TZERO C(K,I)=TZERO ENDDO ENDIF C Q=H*H/TWELVE EQ=Q*(E2-ECH2(I)) CQ=Q*CC2(I) Q=Q*TWO KPM=KPMAX S(1,I)=S1 S(2,I)=S2 R=RZERO U1=V(ONE/R) R=R+H U2=V(ONE/R) C IF(INOUT.EQ.0)THEN C(1,I)=C1 C(2,I)=C2 DO K=3,KPM R=R+H U3=V(ONE/R) D3=ONE/(ONE+U3) D2=(TWO-TEN*U2)*D3 D1=(ONE+U1)*D3 S(K,I)=D2*S(K-1,I)-D1*S(K-2,I) C(K,I)=D2*C(K-1,I)-D1*C(K-2,I) U1=U2 U2=U3 ENDDO ELSE DO K=3,KPM R=R+H U3=V(ONE/R) S(K,I)=((TWO-TEN*U2)*S(K-1,I)-(ONE+U1)*S(K-2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO C(KP2,I)=C1 C(KP2-1,I)=C2 R=RTWO U1=V(ONE/R) R=R-H U2=V(ONE/R) DO K=KP2-2,1,-1 R=R-H U3=V(ONE/R) C(K,I)=((TWO-TEN*U2)*C(K+1,I)-(ONE+U1)*C(K+2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ENDIF C RETURN END C********************************************************************* C SUBROUTINE OPDOUT(DEER,M11,MFIN,DLI,DLVECT) C C CALCULATES THE OUTER REGION CONTRIBUTIONS TO THE DIPOLE LENGTH C MATRIX ELEMENTS C ONLY ENTERED FOR PRODUCTION RUNS C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (ICH2=MZCHF*MZCHF) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF) X,RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/CPERT/CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2) X,DMATV(MZCHF,MZCHF,2) 1 ,WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT COMMON /VUSTOR/ VSTORE(-3:1,MZCHF,MZCHF),USTORE(-3:1,MZCHF,MZCHF), * X1(MZCHF,MZEST),X2(MZCHF),EXISTV(-3:1,MZCHF,MZCHF) COMMON/N1N2/N1N2A,N1N2B,N1AM(ICH2),N2AM(ICH2),N1BM(ICH2), 1 N2BM(ICH2) C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) C LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT LOGICAL EXISTV,SKIPPY C DIMENSION DEER(MZEST),DLI(MZCHF,MZEST),DLVECT(MZCHF,MZEST) DIMENSION FLIJ(MZEST),FIFJ(MZEST),S1(MZEST),S2(MZEST),S3(MZEST) DIMENSION MACONT(ICH2),MBCONT(ICH2),MAC(MZEST,ICH2),MBC(MZEST,ICH2 X) DIMENSION GP1(MZEST),GPP1(MZEST),GPPP1(MZEST),GDEER(MZEST),GEP(MZE XST) DIMENSION GDL(MZEST,MZCHF) C C NRB: C SEE NOTES IN DOUT ABOUT VARIOUS DELTAS AND EWIDTH. C DELTAM=DELTA DELTAP=DELTA C C FIRST ADD IN THE INTERNAL CONTRIBUTION C DO JCHOP=1,NCHOP DO M1=M11,MFIN DLVECT(JCHOP,M1)=DLI(JCHOP,M1) ENDDO ENDDO C C ADD IN THE PART WHICH DOES NOT VECTORISE OVER BOUND LEVELS C THE LOOPS OVER NCHAND1 AND NCHAND2 HAVE BEEN COMBINED INTO C A SINGLE LOOP OVER N1N2A FOR AMAT CONTRIBUTIONS AND OVER C N1N2B FOR BMATL CONTRIBUTIONS. SEE SUBROUTINE RDV. C C SIMILARLY IN THE NEXT SECTIONS THE LOOP OVER M1 IS INDEXED C SO THAT THE CONTRIBUTION WHEN FNU1-FNU2>DELTA CAN BE C VECTORIZED. C DO N12A=1,N1N2A MACONT(N12A)=0 ENDDO DO N12B=1,N1N2B MBCONT(N12B)=0 ENDDO C DO M1=M11,MFIN IF(MZSA1.EQ.1) X READ(4,REC=M1)((FSD1(K,N1,1),K=1,KPMM(M1)),N1=1,NCHND1) D4=FOUR*DEER(M1)*DEER(M1) C DO N12A=1,N1N2A N1=N1AM(N12A) N2=N2AM(N12A) C IF(EWIDTH.GE.TZERO)THEN DELTAM=DELTA*(FNU1(N1,M1)/THREE)**3 IF(FNU2(N2).LT.TZERO)THEN DELTAP=-1.0E10 ELSE DELTAP=1.0E10 ENDIF ENDIF DELTAN=FNU2(N2)-FNU1(N1,M1) C AML=AMATL(N1,N2) IF(DELTAN.LT.-DELTAM.OR.DELTAN.GT.DELTAP) THEN AMD4=AML*D4 VMAD4=OPVMAT(-2,N1,N2,M1)*AMD4 IF(N2.LE.NCHOP)THEN DLVECT(N2,M1)=DLVECT(N2,M1)+(USTORE(-2,N1,N2)*AMD4) ENDIF DO JCHOP=1,NCHOP DLVECT(JCHOP,M1)=DLVECT(JCHOP,M1)+VMAD4*RK(N2,JCHOP) ENDDO MACONT(N12A)=MACONT(N12A)+1 MAC(MACONT(N12A),N12A)=M1 ELSE VMA=OPVMAT(1,N1,N2,M1)*AML IF(N2.LE.NCHOP)THEN DLVECT(N2,M1)=DLVECT(N2,M1)+USTORE(1,N1,N2)*AML ENDIF DO JCHOP=1,NCHOP DLVECT(JCHOP,M1)=DLVECT(JCHOP,M1)+VMA*RK(N2,JCHOP) ENDDO ENDIF ENDDO C DO N12B=1,N1N2B N1=N1BM(N12B) N2=N2BM(N12B) C IF(EWIDTH.GE.TZERO)THEN DELTAM=DELTA*(FNU1(N1,M1)/THREE)**3 IF(FNU2(N2).LT.TZERO)THEN DELTAP=-1.0E10 ELSE DELTAP=1.0E10 ENDIF ENDIF DELTAN=FNU2(N2)-FNU1(N1,M1) C BML=BMATL(N1,N2) IF(DELTAN.LT.-DELTAM.OR.DELTAN.GT.DELTAP) THEN MB12=MBCONT(N12B) MB12=MB12+1 MBCONT(N12B)=MB12 MBC(MB12,N12B)=M1 ELSE VMB=OPVMAT(0,N1,N2,M1)*BML IF(N2.LE.NCHOP)THEN DLVECT(N2,M1)=DLVECT(N2,M1)+USTORE(0,N1,N2)*BML END IF DO JCHOP=1,NCHOP DLVECT(JCHOP,M1)=DLVECT(JCHOP,M1)+VMB*RK(N2,JCHOP) ENDDO ENDIF ENDDO C ENDDO C RZERO2=ONE/(RZERO*RZERO) C C NOW ADD IN THE CONTRIBUTIONS WHICH VECTORIZE OVER BOUND LEVELS C FNU1-FNU2>DELTA C C IF (KKP) THEN C C DO N12A=1,N1N2A N1=N1AM(N12A) N2=N2AM(N12A) AML=AMATL(N1,N2) FKIJ=(CC1(N1)-CC2(N2))*RZERO2 MA12=MACONT(N12A) C C GATHER C DO MA=1,MA12 M1=MAC(MA,N12A) GP1(MA)=P1(N1,M1) GPP1(MA)=PP1(N1,M1) GDEER(MA)=DEER(M1) ENDDO DO JCHOP=1,NCHOP DO MA=1,MA12 S3(MA)=(GP1(MA)*FP2(N2,JCHOP))-(GPP1(MA)*F2(N2,JCHOP)) FLIJ(MA)=S3(MA)*RZERO FIFJ(MA)=GP1(MA)*F2(N2,JCHOP) S1(MA)=FLIJ(MA)+FIFJ(MA) S2(MA)=(FKIJ*(FLIJ(MA)-FIFJ(MA)))-(TWO*((GP1(MA)*FPP2(N2 X ,JCHOP))-(GPP1(MA)*FP2(N2,JCHOP)))) GDL(MA,JCHOP)= X AML*GDEER(MA)*(S1(MA)-(S2(MA)*GDEER(MA))) ENDDO ENDDO DO JCHOP=1,NCHOP DO MA=1,MA12 M1=MAC(MA,N12A) DLVECT(JCHOP,M1)=DLVECT(JCHOP,M1)+GDL(MA,JCHOP) ENDDO ENDDO ENDDO C DO N12B=1,N1N2B N1=N1BM(N12B) N2=N2BM(N12B) BML=BMATL(N1,N2) MB12=MBCONT(N12B) DO MB=1,MB12 M1=MBC(MB,N12B) GP1(MB)=P1(N1,M1) GPP1(MB)=PP1(N1,M1) GEP(MB)=EPS1(N1,M1)-EPS2(N2) ENDDO DO JCHOP=1,NCHOP DO MB=1,MB12 S3(MB)=(GP1(MB)*FP2(N2,JCHOP))-(GPP1(MB)*F2(N2,JCHOP)) GDL(MB,JCHOP)=BML*S3(MB)/GEP(MB) ENDDO ENDDO DO JCHOP=1,NCHOP DO MB=1,MB12 M1=MBC(MB,N12B) DLVECT(JCHOP,M1)=DLVECT(JCHOP,M1)-GDL(MB,JCHOP) ENDDO ENDDO ENDDO C C ELSE C C DO N12A=1,N1N2A N1=N1AM(N12A) N2=N2AM(N12A) AML=AMATL(N1,N2) FKIJ=(CC1(N1)-CC2(N2))*RZERO2 MA12=MACONT(N12A) C C GATHER C DO MA=1,MA12 M1=MAC(MA,N12A) GP1(MA)=P1(N1,M1) GPP1(MA)=PP1(N1,M1) GPPP1(MA)=PPP1(N1,M1) GDEER(MA)=DEER(M1) ENDDO DO JCHOP=1,NCHOP DO MA=1,MA12 S3(MA)=(GP1(MA)*FP2(N2,JCHOP))-(GPP1(MA)*F2(N2,JCHOP)) FLIJ(MA)=S3(MA)*RZERO FIFJ(MA)=GP1(MA)*F2(N2,JCHOP) S1(MA)=FLIJ(MA)-FIFJ(MA) S2(MA)=FKIJ*(FLIJ(MA)+FIFJ(MA))-(TWO*((F2(N2,JCHOP)* X GPPP1(MA))-(FP2(N2,JCHOP)*GPP1(MA)))) GDL(MA,JCHOP)= X AML*GDEER(MA)*(S1(MA)-(S2(MA)*GDEER(MA))) ENDDO ENDDO DO JCHOP=1,NCHOP DO MA=1,MA12 M1=MAC(MA,N12A) DLVECT(JCHOP,M1)=DLVECT(JCHOP,M1)+GDL(MA,JCHOP) ENDDO ENDDO ENDDO C DO N12B=1,N1N2B N1=N1BM(N12B) N2=N2BM(N12B) BML=BMATL(N1,N2) MB12=MBCONT(N12B) C C GATHER C DO MB=1,MB12 M1=MBC(MB,N12B) GP1(MB)=P1(N1,M1) GPP1(MB)=PP1(N1,M1) GEP(MB)=EPS1(N1,M1)-EPS2(N2) ENDDO DO JCHOP=1,NCHOP DO MB=1,MB12 S3(MB)=(GP1(MB)*FP2(N2,JCHOP))-(GPP1(MB)*F2(N2,JCHOP)) GDL(MB,JCHOP)=BML*S3(MB)/GEP(MB) ENDDO ENDDO DO JCHOP=1,NCHOP DO MB=1,MB12 M1=MBC(MB,N12B) DLVECT(JCHOP,M1)=DLVECT(JCHOP,M1)-GDL(MB,JCHOP) ENDDO ENDDO ENDDO C C END IF C C RETURN END C********************************************************************** C SUBROUTINE OPSTRA(DEE,M11,MFIN,DLVECT,DL2,DLTOT,QJUMP,FLAGL) C C CHANGE TO S-MATRIX NORMALIZATION; RETURN THE SQUARES OF DIPOLE C MATRIX ELEMENTS AND THE SUMMED TOTAL C IMPLICIT REAL*8 (A-H,O-Y) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TZERO=0.0) C LOGICAL QJUMP,SKIPPY,FLAGL C COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF) X,RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) C DIMENSION DEE(MZEST),DLVECT(MZCHF,MZEST),DLTOT(MZEST) DIMENSION DLSR(MZCHF,MZEST),DLSI(MZCHF,MZEST),DL2(MZCHF) C SAVE NCHOPS C C IF(.NOT.QJUMP)THEN C C FIRST TIME THRU, TRANSFORM DIPOLE MATRIX ELEMENT VIA D(KHI)=Y**(T*)D(R) C AND WRITE IT TO SCRATCH FILE, THEN RETURN. HERE (Y)=-(YM)-I(XM). C PUT THE REAL AND IMAGINARY PARTS IN DLSR AND DLSI VECTORS C IF(IQDT.GT.0)WRITE(70)FLAGL C DO M1=M11,MFIN DLTOT(M1)=TZERO ENDDO C CCRAY CRAY DO J=1,NCHOP CRAY DO M1=M11,MFIN CRAY DLSR(J,M1)=TZERO CRAY DLSI(J,M1)=TZERO CRAY ENDDO CRAY ENDDO CRAY DO I=1,NCHOP CRAY DO J=1,NCHOP CRAY DO M1=M11,MFIN CRAY DLSR(J,M1)=DLSR(J,M1)-DLVECT(I,M1)*YM(I,J) CRAY DLSI(J,M1)=DLSI(J,M1)+DLVECT(I,M1)*XM(I,J) CRAY ENDDO CRAY ENDDO CRAY ENDDO CRAY DO J=1,NCHOP CRAY DO M1=M11,MFIN CRAY DLTOT(M1)=DLTOT(M1)+DLSR(J,M1)**2+DLSI(J,M1)**2 CRAY ENDDO CRAY ENDDO CCRAY C C--- DO M1=M11,MFIN DO J=1,NCHOP DLSR(J,M1)=TZERO DLSI(J,M1)=TZERO ENDDO DO I=1,NCHOP DO J=1,NCHOP DLSR(J,M1)=DLSR(J,M1)-DLVECT(I,M1)*YM(J,I) DLSI(J,M1)=DLSI(J,M1)+DLVECT(I,M1)*XM(J,I) ENDDO ENDDO ENDDO DO M1=M11,MFIN DO J=1,NCHOP DLTOT(M1)=DLTOT(M1)+DLSR(J,M1)**2+DLSI(J,M1)**2 ENDDO ENDDO C--- DO M1=M11,MFIN DLTOT(M1)=DLTOT(M1)*DEE(M1) ENDDO C IF(IQDT.GT.0)THEN DO M1=M11,MFIN WRITE(70)(DLSR(J,M1),J=1,NCHOP) WRITE(70)(DLSI(J,M1),J=1,NCHOP) ENDDO ENDIF C ENDIF C C IF(IQDT.GT.0)THEN C NRB: C NOT VECTORIZED (BUT FAST) C IF(.NOT.QJUMP)NCHOPS=NCHOP C DO M1=M11,MFIN NCHOP=NCHOPS C IF(EWIDTH.GE.TZERO)THEN CALL PMQDT(DLVECT(1,M1),DLVECT(1,M1),.TRUE.,.FALSE.,DL2,DL2,M1) ELSE CALL DMQDT(DLVECT(1,M1),DLVECT(1,M1),.TRUE.,.FALSE.,DL2,DL2,M1) ENDIF C DLTOT(M1)=TZERO DO J=1,NCHOP DLTOT(M1)=DLTOT(M1)+DL2(J) ENDDO DLTOT(M1)=DLTOT(M1)*DEE(M1) ENDDO C ENDIF C RETURN END C********************************************************************* REAL*8 FUNCTION OPVMAT(I,N1,N2,M1) C C CALCULATES V INTEGRAL, INTEGRAL OF C(N2)*R**I*P(N1) OR C THETA(N2)*R**I*P(N1) C IMPLICIT REAL*8 (A-H,O-Y) C IMPLICIT COMPLEX*16 (Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TZERO=0.0) PARAMETER (HALF=0.5) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (FOUR=4.0) C PARAMETER (NCHF1=MZCHF+1) PARAMETER (MNPEXT=MZMNP+MZCHF) C LOGICAL FLG1,FLG2,EXIST1,EXIST2,EXISTV,DROP,QDT,BMQDT,SKIPPY C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CHAN/ECH(NCHF1),EPS(NCHF1),FKNU(NCHF1),CC(NCHF1), 1 RINF(NCHF1),LLCH(NCHF1),ITARG(NCHF1),NCHF,NCHOP0,NCHOP1 COMMON/FLAGIN/FLG1(MZCHF,MZEST),FLG2(MZCHF),EXIST1(MZCHF,MZEST), 1 EXIST2(MZCHF),DROP COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON /VUSTOR/ VSTORE(-3:1,MZCHF,MZCHF),USTORE(-3:1,MZCHF,MZCHF), * X1(MZCHF,MZEST),X2(MZCHF),EXISTV(-3:1,MZCHF,MZCHF) COMMON/CACC/AC0,ACNUM,ACJWBK,ACZP,LACC COMMON/CTHET/BB(NCHF1,MZTET),BG(NCHF1,MZTET),MSUM(NCHF1) COMMON/CJWBK/D(15) COMMON/FLAGQ/QDT COMMON/CVMAT/TI(MZCHF,MZEST),TIP(MZCHF,MZEST),TID(MZCHF,MZEST), + TIDP(MZCHF,MZEST),PIN(MZCHF,MZEST), + TF(MZCHF),TFP(MZCHF),TFD(MZCHF),TFDP(MZCHF),PFN(MZCHF), * BB2(MZCHF,MZTET),BB1(MZCHF,MZTET,MZEST),D2(15,MZCHF), * MSUM1(MZCHF,MZEST),MSUM2(MZCHF) C COMMON/NRBFLG/TOL0V,IFLAGV COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) C DIMENSION TEMP(MZPTS) DATA NLAG/6/ C C C C IF N2 IS WC CHANNEL CHECK CONVERGENCE OF INTEGRAL C DELTAN=DELTA BMQDT=QDT.OR.IQDT.NE.0 IF(BMQDT)DELTAN=HALF IF(BMQDT.AND.N2.LE.NCHOP.AND.EPS2(N2).LT.TZERO) THEN IF(FNU2(N2).LE.FNU1(N1,M1)+DELTAN) THEN IF(IQDT.LE.0)THEN WRITE(6,600) N1,N2,FNU1(N1,M1),FNU2(N2) CALL ABORTT(14) ENDIF SKIPPY=.TRUE. OPVMAT=TZERO RETURN ENDIF ENDIF C C RZERO TO RTWO C KPM=KPMM(M1) OPVMAT=TZERO UM=TZERO MA=M1 IF(MZSA1.EQ.1)MA=1 DO K=1,KPM TEMP(K)=FSD1(K,N1,MA)*WBODE(K,I) ENDDO DO K=1,KPM OPVMAT=OPVMAT+TEMP(K)*C(K,N2) ENDDO IF(N2.LE.NCHOP) THEN DO K=1,KPM UM=UM+TEMP(K)*S(K,N2) ENDDO USTORE(I,N1,N2)=UM ENDIF C IF(.NOT.(FLG1(N1,M1).AND.FLG2(N2))) GO TO 2 C C SEE IF FLAG HAS BEEN SET TO SKIP RTWO TO INFINITY C IF(IFLAGV.EQ.0)GO TO 2 C C RTWO TO INFINITY C NCHOP0=NCHOP C C SET UP INITIAL CHANNEL IF NECESSARY C IF(EXIST1(N1,M1)) THEN MS=MSUM1(N1,M1) MSUM(N1)=MS DO K=1,MS BB(N1,K)=BB1(N1,K,M1) ENDDO ELSE EPS(N1)=EPS1(N1,M1) FKNU(N1)=FNU1(N1,M1) LLCH(N1)=(SQRT(ONE+FOUR*CC1(N1))-ONE+DELTA)*HALF CALL THETA(RTWO,N1,T,TP,TD,TDP) IF(ABS(T).LT.TOL0V)THEN X1(N1,M1)=TZERO !NRB ELSE MS=MSUM(N1) MSUM1(N1,M1)=MS DO K=1,MS BB1(N1,K,M1)=BB(N1,K) ENDDO X1(N1,M1)=FSD1(KP2,N1,MA)/T TI(N1,M1)=T TIP(N1,M1)=TP TID(N1,M1)=TD TIDP(N1,M1)=TDP PIN(N1,M1)=(RTWO**(.5*FKNU(N1)))*EXP(-.5*RTWO 1 /FKNU(N1)) ENDIF EXIST1(N1,M1)=.TRUE. ENDIF C C SET UP FINAL CHANNEL C N21=NCHND1+1 IF(EXIST2(N2)) THEN IF(N2.LE.NCHOP) THEN DO K=1,15 D(K)=D2(K,N2) ENDDO ELSE MS=MSUM2(N2) MSUM(N21)=MS DO K=1,MS BB(N21,K)=BB2(N2,K) ENDDO ENDIF ELSE EP2=EPS2(N2) EPS(N21)=EP2 L=(SQRT(ONE+FOUR*CC2(N2))-ONE+DELTA)*HALF LLCH(N21)=L IF(N2.LE.NCHOP) THEN FKNU(N21)=SQRT(EP2) CALL INJWBK(EP2,L) DO K=1,15 D2(K,N2)=D(K) ENDDO ELSE FKNU(N21)=FNU2(N2) CALL THETA(RTWO,N21,T,TP,TD,TDP) IF(ABS(T).LT.TOL0V)THEN X2(N2)=TZERO !NRB ELSE X2(N2)=C(KP2,N2)/T MS=MSUM(N21) MSUM2(N2)=MS DO K=1,MS BB2(N2,K)=BB(N21,K) ENDDO TF(N2)=T TFP(N2)=TP TFD(N2)=TD TFDP(N2)=TDP PFN(N2)=(RTWO**(HALF*FNU2(N2)))*EXP(-HALF*RTWO/FNU2(N2)) ENDIF ENDIF EXIST2(N2)=.TRUE. ENDIF C XN1=X1(N1,M1) XN2=X2(N2) IF(ABS(XN1).LT.TOL0V.OR.ABS(XN2).LT.TOL0V)GO TO 2 !NRB IF(N2.LE.NCHOP) THEN FKNU(N1)=FNU1(N1,M1) FKNU(N21)=SQRT(EPS2(N2)) CALL ZSLAG(N21,N1,NLAG,-I,ZS3,ZSD3) OPVMAT=OPVMAT+DBLE(ZS3)*XN1 C STORE UMAT FOR RTWO TO INFINITY, SPIN OFF FROM ZSLAG USTORE(I,N1,N2)=UM+DIMAG(ZS3)*XN1 ELSE IF(I.NE.0)THEN CALL TLAG(N1,N21,NLAG,-I,T1,T2,T3) ELSE T1=HALF*((PIN(N1,M1)/PFN(N2))* X (TI(N1,M1)*TIDP(N1,M1)-TID(N1,M1)*TIP(N1,M1))+ X (PFN(N2)/PIN(N1,M1))* X (TF(N2)*TFDP(N2)-TFD(N2)*TFP(N2))) ENDIF OPVMAT=OPVMAT+T1*XN1*XN2 ENDIF C 2 RETURN 600 FORMAT(//10X,'*** OUTER REGION INTEGRAL WILL NOT CONVERGE ***'// 1 26X,'BOUND',10X,'FREE'/26X,5('-'),10X,4('-')/ 2 ' CHANNEL',18X,I3,12X,I3/ 3 ' EFFECTIVE QUANTUM NUMBER',F10.4,5X,F10.4/) END C C********************************************************************* C SUBROUTINE PERTAR C C SETS UP CMAT ARRAY FOR PERTURBATION CALCULATIONS C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CPERT/CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2) X,DMATV(MZCHF,MZCHF,2) 1 ,WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/NRBZED/TZED C DO I2=1,NCHND2 DO I1=1,NCHND1 CMAT(I1,I2)=0.0 ENDDO ENDDO C DO I2=1,NCHND2 DO J1=1,NCHND1 T=(CC1(J1)-CC2(I2))*AMATL(J1,I2) DO I1=1,NCHND1 CMAT(I1,I2)=CMAT(I1,I2)+WMAT1(I1,J1,1)*T ENDDO ENDDO ENDDO C DO I2=1,NCHND2 DO J2=1,NCHND2 DO I1=1,NCHND1 CMAT(I1,I2)=CMAT(I1,I2) X -(4.*TZED+CC1(I1)-CC2(J2))*AMATL(I1,J2)*WMAT2(J2,I2,1) ENDDO ENDDO ENDDO C RETURN END C********************************************************************* C SUBROUTINE PMQDT(DLVECT,DVVECT,FLAGL,FLAGV,DL2,DV2,M1) C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C CBL PARAMETER (MWORK=100*MZDEG) PARAMETER (TZERO=0.0) PARAMETER (ZERO=(0.0,0.0)) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (FOUR=4.0) PARAMETER (BIG=150.0) PARAMETER (TINE=1.0E-6) C LOGICAL FLAGV,QDT,FLAGL,SKIPPY C C NRB: C EVALUATE TOTAL PRE-CONVOLUTED PHOTOIONIZATION USING MULTI-STATE QDT C C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF), X RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAGQ/QDT COMMON/TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) C COMMON/NRBKHI/ZKHICC(MZDEG,MZDEG),ZKHIOC(MZCHF,MZDEG),ZVAL(MZDEG) X ,ZDL(MZCHF),ZDV(MZCHF),ZDL0(MZCHF),ZDV0(MZCHF) CBL X,ZWL(MZDEG,MZDEG),ZWR(MZDEG,MZDEG),ZWORK(MWORK),RWORK(2*MZDEG) COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBZDD/ZDLM(MZCHF,MZEST),ZDVM(MZCHF,MZEST) C DIMENSION IPIV(MZDEG) DIMENSION DLVECT(MZCHF),DVVECT(MZCHF),DL2(MZCHF),DV2(MZCHF) DIMENSION XMOLD(MZCHF,MZCHF),YMOLD(MZCHF,MZCHF) C EQUIVALENCE (XMOLD,DFP2),(YMOLD,DFPP2) C CBL DATA LWMAX/MWORK/ C C C C ALL CHANNELS C NCHF=NCHND2 C C FIRST DETERMINE NCHOP (NUMBER OF OPEN CHANNELS) C DO J=NCHF,1,-1 IF(EPS2(J).GE.TZERO) THEN NCHOP=J GO TO 1 ENDIF ENDDO C C CASE OF NO REAL OPEN CHANNELS (THIS WOULD BE PHOTOABSORPTION) C NCHOP=0 RETURN C 1 IF(QDT)THEN WRITE(6,*)' ***SR.PMQDT: WARNING, NO GAILITIS AVERAGE & NO ' X,'DAMPING ALLOWED/USED FOR PRECONVOLUTED TOTAL PHOTOIONIZATION' ENDIF C PI=ACOS(-ONE) TPI=TWO*PI C C INITIALIZE KHICC; USE (KHI) = CONJG(2I(Y)-1) = (2(XM)-1,2(YM)) C IF(NCHF.GT.MZDEG) THEN WRITE(6,601)NCHF CALL ABORTT(10) ENDIF C DO N2=1,NCHF DO N1=1,NCHF ZKHICC(N1,N2)=DCMPLX(TWO*XM(N1,N2), X TWO*YM(N1,N2)) ZKHIOC(N1,N2)=ZKHICC(N1,N2) ENDDO ZKHICC(N2,N2)=ZKHICC(N2,N2)-ONE ZKHIOC(N2,N2)=ZKHICC(N2,N2) ENDDO C C ZKHICC-ZFDEC; ZKHIOC+ZFDEC C (NOTE THE SIGN DIFFERENCE ON THE WIDTH COMPARED TO DMQDT) C DO N=NCHOP+1,NCHF TPINU=FNU2(N)*TPI Z=DCMPLX(TZERO,-TPINU) IF(EWIDTH.EQ.TZERO)THEN ZFDEC=EXP(Z) ELSEIF(IRDEC.EQ.1)THEN !BELL&SEATON T=MIN(EWIDTH*FNU2(N)**2/TWO,BIG) ZB=DCMPLX(ONE,-T) ZFDEC=EXP(Z*ZB) ELSE !HICKMAN-ROBICHEAUX (& CASE IRDEC=0) ZFNU=DCMPLX(ONE,TZERO)/SQRT(DCMPLX(ONE/FNU2(N)**2,-EWIDTH)) Z=DCMPLX(TZERO,-TPI)*CONJG(ZFNU) FR=DBLE(Z) FI=DIMAG(Z) FR=MIN(FR,BIG) ZFDEC=EXP(DCMPLX(FR,FI)) C C ALPHA CORRECTION FACTOR (NEGLIGIBLE). C REQUIRES COMPLEX DIGAMMA FUNCTION ZPSI. C CALP Z=-ONE/ZFNU**2 CALP ZA=DCMPLX(ONE,TZERO) CALP L=(SQRT(ONE+FOUR*CC2(N))-ONE+TINE)/TWO CALP DO I=1,L CALP ZA=ZA*(ONE+I*I*Z) CALP ENDDO C CALP Z=ZFNU+L+1 CALP ZY1=ZPSI(Z) CALP Z=ZFNU-L CALP ZY2=ZPSI(Z) CALP Z=TWO*LOG(ZFNU) C CALP ZG=ZA*(ZY1+ZY2-Z) CALP ZG=ZG/TPI C CALP ZALF=IMAG(ZG)/ZA CALP ZALF=CONJG(ZALF) C CALP ZFDEC=ZFDEC*(ONE+ZALF)+ZALF C C END ALPHA CORRECTION FACTOR C ENDIF C ZKHICC(N,N)=ZKHICC(N,N)-ZFDEC ZKHIOC(N,N)=ZKHIOC(N,N)+ZFDEC ENDDO C C INVERT C CSTRTNBL CALL ZVERT(ZKHICC,MZDEG,NCHF,IPIV) CENDNBL C CSTRTBL CBL CALL ZGETRF(NCHF,NCHF,ZKHICC,MZDEG,IPIV,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,602) INFO CBL STOP CBL ENDIF CBL CALL ZGETRI(NCHF,ZKHICC,MZDEG,IPIV,ZWORK,MWORK,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,603) INFO CBL STOP CBL ENDIF CBL LWOPT=INT(ZWORK(1)) CBL IF (LWOPT.GT.MWORK) WRITE(6,604) LWOPT CENDBL C C RE-ZERO-OUT DIAGONAL OF OMITTED CHANNELS (SHOULD BE UNNECESSARY HERE, C BUT JUST TO BE SAFE INCASE SOMEONE GOES AND USES THEM LATER). C DO N=1,NCHF IF(IOMIT(N).NE.0)THEN ZKHICC(N,N)=ZERO ZKHIOC(N,N)=ZERO ENDIF ENDDO C C DIPOLE MATRIX: WE WILL ALREADY HAVE IT, RECOVERED AND C INTERPOLATED IN READXY. C IF(FLAGL)THEN DO J=1,NCHF ZDL(J)=ZDLM(J,M1) ENDDO ENDIF IF(FLAGV)THEN DO J=1,NCHF ZDV(J)=ZDVM(J,M1) ENDDO ENDIF C C FORM PHYSICAL DIPOLE MATRIX C IF(FLAGL)THEN DO J=1,NCHF ZDL0(J)=ZERO ENDDO DO K=1,NCHF DO J=1,NCHF ZDL0(J)=ZDL0(J)+ZKHICC(J,K)*ZDL(K) ENDDO ENDDO DO J=1,NCHF ZVAL(J)=ZERO ENDDO DO K=1,NCHF DO J=1,NCHF ZVAL(J)=ZVAL(J)+ZKHIOC(J,K)*ZDL0(K) ENDDO ENDDO DLTOT=TZERO DO J=1,NCHF DLTOT=DLTOT+DBLE(CONJG(ZDL(J))*ZVAL(J)) DL2(J)=TZERO ENDDO DL2(1)=DLTOT ENDIF C IF(FLAGV)THEN DO J=1,NCHF ZDV0(J)=ZERO ENDDO DO K=1,NCHF DO J=1,NCHF ZDV0(J)=ZDV0(J)+ZKHICC(J,K)*ZDV(K) ENDDO ENDDO DO J=1,NCHF ZVAL(J)=ZERO ENDDO DO K=1,NCHF DO J=1,NCHF ZVAL(J)=ZVAL(J)+ZKHIOC(J,K)*ZDV0(K) ENDDO ENDDO DVTOT=TZERO DO J=1,NCHF DVTOT=DVTOT+DBLE(CONJG(ZDV(J))*ZVAL(J)) DV2(J)=TZERO ENDDO DV2(1)=DVTOT ENDIF C NCHOP=1 C RETURN C 601 FORMAT(//' ***NUMBER OF CHANNELS, NCHF = ',I2/4X, 1 'LARGER THAN DIMENSION VALUE OF DEG = MZDEG'//) CB602 FORMAT(//10X,10('*'),' SR. PMQDT: ZGETRF RETURNED WITH INFO =',I2) CB603 FORMAT(//10X,10('*'),' SR. PMQDT: ZGETRI RETURNED WITH INFO =',I2) CB604 FORMAT(//10X,10('*'),' SR. PMQDT: MWORK SHOULD BE LARGER THAN ', CB X 'WORK(1) = ',F5.0,' FOR OPTIMAL PERFORMANCE') END C C*********************************************************************** C SUBROUTINE R1EDEP(M1,MNP2,NCHF) C C READ INITIAL ENERGY DEPENDENT DATA C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) C COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF) X,RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/NRBPRT/IPERT1,IPERT2 DIMENSION XVECT(MZCHF) C TOL0B=-1.D-150 !INITIALIZE DEEPLY-CLOSED AS OFF (ZERO) C READ(1,ERR=100,END=100) EE1(M1),(XVECT(N),N=1,NCHF) TOL0B=ABS(TOL0B) !SWITCH-ON AS CAN RENORM (NEW STGB/BMN) 100 READ(1) READ(1) (P1(N,M1),N=1,NCHF) READ(1) (PP1(N,M1),N=1,NCHF) READ(1) (PPP1(N,M1),N=1,NCHF) READ(1)KPM KPMM(M1)=KPM MA=M1 IF(MZSA1 .EQ. 1) MA=1 IF(KPM.GT.1)THEN READ(1)(FSD1(KPM-1,I,MA),FSD1(KPM,I,MA),I=1,NCHF) CALL NUMB(M1,MA,NCHF,XVECT,TOL0B) ELSE READ(1)(FSD1(1,I,MA),I=1,NCHF) ENDIF IF (MZSA1.EQ.1) WRITE(4,REC=M1)((FSD1(K,N1,1),K=1,KPM),N1=1,NCHF) IF(IPERT1.EQ.1) THEN READ(1)(DP1(N,M1),N=1,NCHF) READ(1)(DPP1(N,M1),N=1,NCHF) READ(1)(DPPP1(N,M1),N=1,NCHF) ENDIF C RETURN END C C******************************************************************* C SUBROUTINE R1EIND(FILE1,M11,M12,M13,MNP2,NCHF) C C READS INITIAL ENERGY INDEPENDENT DATA C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C CHARACTER*3 FILE1 C COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/CPERT/CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2) X,DMATV(MZCHF,MZCHF,2) 1 ,WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/NRBPRT/IPERT1,IPERT2 C OPEN(1,FILE=FILE1,STATUS='OLD',FORM='UNFORMATTED',ERR=90) REWIND(1) GOTO 6 90 WRITE(6,696) FILE1 CALL ABORTT(13) C C READ E1-INDEPENDENT DATA 6 READ(1)IS,IL,IP READ(1)MNP2,NCHF IF(MNP2.GT.MZMNP) THEN WRITE(6,613)MNP2 CALL ABORTT(10) ENDIF IF(NCHF.GT.MZCHF) THEN WRITE(6,614)NCHF CALL ABORTT(10) ENDIF READ(1)(ECH1(I),I=1,NCHF) READ(1)(CC1(I),I=1,NCHF) IF(IPERT1.EQ.1)READ(1)(((WMAT1(I,J,L),I=1,NCHF),J=1,NCHF),L=1,2) READ(1)MXE1 C IF(M11.GT.MXE1) THEN WRITE(6,505)M11,MXE1 C CALL ABORTT(7) ENDIF C C SKIP PRECEDING LEVELS; NREC HOLDS THE NUMBER OF RECORDS FOR A LEVEL M13=MIN(M12,MXE1) NREC=7 !7 NOT 6 - NRB IF(IPERT1.EQ.1) NREC=NREC+3 M1SKIP=(M11-1)*NREC DO 16 I=1,M1SKIP 16 READ(1) DUMMY C RETURN 505 FORMAT(//' BOUND LEVEL REQUESTED',I4,10X, * 'NUMBER OF LEVELS FOR THIS SLPI',I4) 613 FORMAT(//' READS MNP2 =',I5,' WHICH IS LARGER THAN ', 1 'MAXIMUM VALUE OF MZMNP ALLOWED BY DIMENSIONS'//) 614 FORMAT(//' READS NCHF =',I5,' WHICH IS LARGER THAN MAXIMUM ', *'VALUE OF MZCHF (LESS 1 FOR FINAL STATE) ALLOWED BY DIMENSIONS'//) 696 FORMAT(/' *** CANNOT OPEN FILE ',A/) END C********************************************************************** C SUBROUTINE R2EDEP(M2,E2,EJUMP) C C READS FINAL ENERGY DEPENDENT DATA C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) CBL PARAMETER (LWORK=MZCHF*MZCHF) PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT,QDT,EJUMP X ,SKIPPY C COMMON/FLAGQ/QDT COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF) X,RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT C COMMON/NRBPRT/IPERT1,IPERT2 COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) C DIMENSION IPIV(MZCHF) CBL DIMENSION WORK(LWORK) CBL EQUIVALENCE (DF2,WORK) C C READ(2) QDT IF(QDT.AND.IQDT.LE.0) IPERT=0 CT READ(2) NCHOP !TEMP FOR OLD STGF CT DO I=1,NCHOP CT IOMIT(I)=0 CT ENDDO READ(2) NCHOP,(IOMIT(I),I=1,NCHOP) OK = NCHOP.GT.0 IF(.NOT.OK) GO TO 90 CW READ(2) ((RK(I,J),J=1,NCHOP),I=1,NCHND2) READ(2) ((F2(I,J),J=1,NCHOP),I=1,NCHND2) READ(2) ((FP2(I,J),J=1,NCHOP),I=1,NCHND2) READ(2) ((FPP2(I,J),J=1,NCHOP),I=1,NCHND2) C IF(EJUMP)GO TO 10 C C THE TRANSFORMATION MATRIX IS (Y)=-(K+I)/(1+K**2) C OR (Y)=-1/2I(1+CONJG(S)) C C XM=(1+RK**2)**(-1) CSTRTNBL CALL MATSQ(RK,XM,NCHOP) CENDNBL CSTRTBL CBL CALL DGEMM('N','N',NCHOP,NCHOP,NCHOP,ONE,RK,MZCHF,RK,MZCHF, CBL X TZERO,XM,MZCHF) CENDBL C DO I=1,NCHOP XM(I,I)=XM(I,I)+ONE ENDDO C CSTRTNBL CALL VERT(XM,MZCHF,NCHOP,IPIV) CENDNBL CSTRTBL CBL CALL DGETRF(NCHOP,NCHOP,XM,MZCHF,IPIV,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,602) INFO CBL STOP CBL ENDIF CBL CALL DGETRI(NCHOP,XM,MZCHF,IPIV,WORK,LWORK,INFO) CBL IF (INFO.NE.0) THEN CBL WRITE(6,603) INFO CBL STOP CBL ENDIF CBL IF (WORK(1).GT.LWORK) WRITE(6,604) WORK(1) CENDBL C DO I=1,NCHOP IF(IOMIT(I).NE.0)XM(I,I)=TZERO ENDDO C C YM=RK*XM CSTRTNBL CALL MATMUL(RK,XM,YM,NCHOP) CENDNBL CSTRTBL CBL CALL DGEMM('N','N',NCHOP,NCHOP,NCHOP,ONE,RK,MZCHF,XM,MZCHF, CBL X TZERO,YM,MZCHF) CENDBL C C WRITE UNPHYSICAL S-MATRIX TO SCRATCH (NCHOP=NCHND2 HERE) C IF(IQDT.GT.0)THEN WRITE(70)NCHOP,E2,M2 WRITE(70)((XM(I,J),I=J,NCHOP),J=1,NCHOP) WRITE(70)((YM(I,J),I=J,NCHOP),J=1,NCHOP) ENDIF C 10 DO I=1,NCHOP READ(2)INOUT,S1,S2,C1,C2 IF(.NOT.EJUMP)CALL NUMFO(E2,I,INOUT,S1,S2,C1,C2) ENDDO DO I=NCHOP+1,NCHND2 READ(2)INOUT,C1,C2 IF(.NOT.EJUMP)CALL NUMFC(E2,I,INOUT,C1,C2) ENDDO IF(IPERT2.NE.0) THEN READ(2,ERR=92)((DF2(I,J),J=1,NCHOP),I=1,NCHND2) READ(2)((DFP2(I,J),J=1,NCHOP),I=1,NCHND2) READ(2)((DFPP2(I,J),J=1,NCHOP),I=1,NCHND2) ENDIF RETURN C C UNABLE TO USE PERTURBATION IN STGF, SO TEMPORARILY SET IPERT=0 C 92 BACKSPACE 2 IF(EJUMP)RETURN IPERT=0 IF(IPRINT.GT.0) WRITE(6,601) C 90 RETURN C 601 FORMAT(' *** SET IPERT=0 TEMPORARILY FOR THE NEXT ENERGY ***') CB602 FORMAT(//10X,10('*'),' SR.R2EDEP: DGETRF RETURNED WITH INFO =',I2) CB603 FORMAT(//10X,10('*'),' SR.R2EDEP: DGETRI RETURNED WITH INFO =',I2) CB604 FORMAT(//10X,10('*'),' SR.R2EDEP: LWORK SHOULD BE LARGER THAN ', CB X 'WORK(1) = ',F5.0,' FOR OPTIMAL PERFORMANCE') END C C********************************************************************* C SUBROUTINE R2EIND(FILE2) C C READS FINAL ENERGY INDEPENDENT DATA C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) C CHARACTER*3 FILE2 C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/RBASIS/VALUE(MZMNP),WMAT(MZCHF,MZMNP) COMMON/CPERT/CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2) X,DMATV(MZCHF,MZCHF,2) 1 ,WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CBETA/ITARG(MZCHF,-1:1),LTARG(MZCHF,-1:1),L2P(MZCHF,-1:1) C ,KJ(MZCHF,-1:1) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/NRBPRT/IPERT1,IPERT2 C OPEN(2,FILE=FILE2,STATUS='OLD',FORM='UNFORMATTED') REWIND(2) READ(2)IS,IL,IP READ(2)MNP2,NCHF IF(MNP2.GT.MZMNP) THEN WRITE(6,613)MNP2 CALL ABORTT(10) ENDIF IF(NCHF.GT.MZCHF) THEN WRITE(6,614)NCHF CALL ABORTT(10) ENDIF READ(2)(ECH2(I),I=1,NCHF) C C --- EXTRA INPUT FROM STGF FOR BETA (KAB, MARCH 1991) K=LRGLD2-LRGLD1 IF(IS.EQ.0)K=K/2 IF(IPRINT.EQ.3.AND.ABS(K).LE.1) THEN READ(2)(CC2(I),I=1,NCHF),(L2P(I,K),I=1,NCHF),(LTARG(I,K),I=1,NCHF) * ,(ITARG(I,K),I=1,NCHF),(KJ(I,K),I=1,NCHF) ELSE READ(2)(CC2(I),I=1,NCHF) ENDIF READ(2)(VALUE(I),I=1,MNP2) READ(2)((WMAT(I,J),J=1,MNP2),I=1,NCHF) IF(IPERT2.NE.0)READ(2)(((WMAT2(I,J,L),I=1,NCHF),J=1,NCHF),L=1,2) C C CHECK ACCORD BETWEEN FILES F AND D C IF(IS.NE.NSPND.OR. 1 IL.NE.LRGLD2.OR. 2 IP.NE.NPTYD2.OR. 3 MNP2.NE.MNP2D2.OR.NCHF.NE.NCHND2)THEN WRITE(6,510)FILE2,IS,IL,IP,MNP2 CALL ABORTT(9) ENDIF C RETURN 510 FORMAT(//' *** MISMATCH'/4X,'FILE ',A3,' HAS IS,IL,IP,MNP2 = ', 1 4I4//) 613 FORMAT(//' READS MNP2 = ',I4,' WHICH IS LARGER THAN ', 1 'MAXIMUM VALUE OF MZMNP ALLOWED BY DIMENSIONS'//) 614 FORMAT(//' READS NCHF = ',I4,' WHICH IS LARGER THAN MAXIMUM ', *'VALUE OF MZCHF (LESS 1 FOR FINAL STATE) ALLOWED BY DIMENSIONS'//) END C********************************************************************** C SUBROUTINE RBF00(EMESH,MXE2) C C READ DIRECTORY FILES ON B AND F DATASETS C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMON/LIST1/KSLP1,MSLP1(100) COMMON/LIST2/KSLP2,MSLP2(100) COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/SCALE/AZ,AZP,AZ2,AZP2,K10,K20 COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON /TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) C COMMON/NRBPRT/IPERT1,IPERT2 COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) C LOGICAL SKIPPY C DIMENSION EMESH(MXE2) C C INITIAL BOUND STATE DATA C OPEN(1,FILE='B00',STATUS='OLD',FORM='UNFORMATTED',ERR=9) C REWIND(1) READ(1)NZED,NELC READ(1)KP2 READ(1)RZERO,H READ(1)(WBODE(IR,0),IR=1,KP2) READ(1)IPERT1 IF(IPERT1.EQ.0)WRITE(6,500) WRITE(6,506)NZED,NELC IF(IPRINT.GT.0)WRITE(6,503)RZERO,H,KP2 IF(IPRINT.GT.0)WRITE(6,501) KSLP1=0 C C AND LIST OF SLPI CASES c 5 READ(1,END=998)IS,IL,IP IF(IL.GE.0)THEN KSLP1=KSLP1+1 MSLP1(KSLP1)=(IS*100+IL)*100+IP IF(IPRINT.GT.0)WRITE(6,502)KSLP1,IS,IL,IP GOTO 5 ENDIF 998 CLOSE(1) C C DATA FOR SCALING c AZ=MAX(NZED-NELC,1) AZP=1./AZ AZ2=AZ*AZ AZP2=AZP*AZP C C CHECK DIMENSION FOR KP2 IF(KP2.GT.MZPTS)THEN WRITE(6,600)KP2,MZPTS STOP ENDIF C C COMPLETE CALCULATION OF WBODE ARRAY C RTWO=RZERO+H*(KP2-1) R=RZERO-H DO K=1,KP2 R=R+H AR=1./R AR2=AR*AR AR3=AR2*AR WBODE(K,-3)=WBODE(K,0)*AR3 WBODE(K,-2)=WBODE(K,0)*AR2 WBODE(K,1)=R*WBODE(K,0) ENDDO C C FINAL FREE STATE DATA C OPEN(2,FILE='F00',STATUS='OLD',FORM='UNFORMATTED',ERR=10) C REWIND(2) READ(2)NZED2,NELC2 IF(IRDEC.EQ.0) !TEMP XREAD(2)NAST,(ENAT(I),I=1,NAST),(ISAT(I),LAT(I),I=1,NAST) IF(IRDEC.NE.0) !TEMP XREAD(2)NAST,(ENAT(I),I=1,NAST),(ISAT(I),LAT(I),I=1,NAST) X,(ARDEC(I),I=1,NAST) READ(2)KP22 READ(2)RZERO2,H2 READ(2)WBODE2 READ(2)IPERT2,IQDT IF(NZED2.NE.NZED.OR. 1 NELC2.NE.NELC.OR. 2 RZERO2.NE.RZERO.OR. 3 KP22.NE.KP2.OR. 4 H2.NE.H) THEN WRITE(6,504)NZED,NZED2,NELC,NELC2,RZERO,RZERO2,KP2,KP22,H,H2 CALL ABORTT(1) ENDIF IF(IPERT2.NE.IPERT1) THEN WRITE(6,507)IPERT1,IPERT2 IPERT=0 ELSE IPERT=IPERT1 ENDIF READ(2)MXE2 C C CHECK DIMENSION FOR MXE C IF(MXE2.GT.MZMSH)THEN WRITE(6,610)MXE2,MZMSH STOP ENDIF IF(IPRINT.EQ.3.AND.MXE2.GT.MZSA3)THEN WRITE(6,617)MXE2,MZSA3 STOP ENDIF C READ(2)(EMESH(I),I=1,MXE2) READ(2)BSTO C C UN-SCALE TARGET ENERGIES TO BE WRITTEN TO ARCHIVE FILE C DO I=1,NAST ENAT(I)=AZ2*ENAT(I) ENDDO IF(IPRINT.GT.0)WRITE(6,505) KSLP2=0 C C AND LIST OF SLPI CASES 6 READ(2,END=997)IS,IL,IP KSLP2=KSLP2+1 MSLP2(KSLP2)=(IS*100+IL)*100+IP IF(IPRINT.GT.0)WRITE(6,502)KSLP2,IS,IL,IP GOTO 6 997 CLOSE(2) C RETURN C 9 CALL ABORTT(12) 10 CALL ABORTT(15) C 500 FORMAT(17X,'(CALCULATION WITHOUT MUTIPOLE POTENTIALS)') 501 FORMAT(//10X,'FROM STGB DATASET'/10X,'KSLP IS IL IP') 502 FORMAT(10X,I3,4X,I2,3X,I2,3X,I2) 503 FORMAT(/7X,'RZERO =',F9.4/7X,'H =',F8.4/7X,'KP2 =',I5/) 504 FORMAT(/30X,'B',20X,'F'/10X,'NZED',2I20/10X,'NELC',2I20/ 1 10X,'RZERO',2F20.2/10X,'KP2',2I20/10X,'H',2F20.2/) 505 FORMAT(//10X,'FROM STGF DATASET'/10X,'KSLP IS IL IP') 506 FORMAT(//14X,'NUCLEAR CHARGE =',I3, * ', NUMBER OF TARGET ELECTRONS =',I3/14X,52('-')//) 507 FORMAT(//'*WARNING*','STGB',5X,'STGF'/5X,'IPERT',I4,5X,I4/ X10X,'PERTURBATIONS SWITCHED-OFF'//) 600 FORMAT(///10X,'KP2 =',I6, * ' IS GREATER THAN MZPTS = ',I5///) 610 FORMAT(///10X,'MXE2 =',I7, * ' IS GREATER THAN MZMSH = ',I7///) 617 FORMAT(///10X,'MXE2 =',I7, * ' IS GREATER THAN MZSA3 = ',I7///) C END C C********************************************************************* C SUBROUTINE RD00 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C C READ DIRECTORY FILE ON D DATASET (UNLESS A FILE D00 EXISTS) C PARAMETER (MXTRAN=99) COMMON/LIST0/KOUNT,MSLPL(MXTRAN),MSLPR(MXTRAN) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC C OPEN(3,FILE='D00',STATUS='OLD',FORM='UNFORMATTED',ERR=4) GO TO 5 4 OPEN(3,FILE='D',STATUS='OLD',FORM='UNFORMATTED',ERR=9) 5 READ(3)KOUNT IF(IPRINT.GT.0)WRITE(6,500) DO 1 K=1,MIN(KOUNT,MXTRAN) READ(3)IS1,IL1,IP1,IS2,IL2,IP2 IF(IPRINT.GT.0)WRITE(6,501)K,IS1,IL1,IP1,IS2,IL2,IP2 MSLPL(K)=(IS1*100+IL1)*100+IP1 1 MSLPR(K)=(IS2*100+IL2)*100+IP2 CLOSE(3) RETURN 9 CALL ABORTT(16) C 500 FORMAT(//10X,'FROM D DATASET'//10X,'K',8X,'IS1 IL1 IP1', * 6X,'IS2 IL2 IP2'/) 501 FORMAT(2(8X,I3),2I5,6X,I3,2I5) C END C C********************************************************************** C SUBROUTINE RD5L1(IS1,IL1,IP1,M11,M12) C C READ INITIAL STATE SPECIFICATION FROM UNIT 5 C VALID INPUT IS EITHER IS1,IL1,IP1,M11,M12 OR EOF C COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT CHARACTER CARD*80 C READ(5,*,END=1,ERR=2)IS1,IL1,IP1,M11,M12 IF(M11.EQ.0) GO TO 3 IF(M12.LT.M11) M12=M11 GO TO 4 C 1 EOF=.TRUE. 3 M11=1 M12=9999 4 RETURN C 2 BACKSPACE 5 READ(5,*)CARD WRITE(6,*)' *** THE FOLLOWING INPUT IS ENCOUNTERED WHEN ' WRITE(6,*)' AN INITIAL STATE SPECIFICATION IS EXPECTED:' WRITE(6,'(A)')CARD CALL ABORTT(5) END C C******************************************************************* C SUBROUTINE RD5L2 C C READ UNIT 5 FOR SECOND SLPI LIST, OR SET ALL2 FLAG C COMMON/SLPI2/IS2V(3),IL2V(3),IP2V(3),KUT COMMON/FLAG0/ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT LOGICAL ALL1,ALL2,RED,MORE1,MORE2,OK,EOF,KKP,CMPLT C KUT=0 DO 2 I=1,4 READ(5,*,END=4) IS,IL,IP IF(I.GT.3) GO TO 6 IP2V(I)=MOD(IP,2) IL2V(I)=IL IS2V(I)=IS IF(IL.EQ.-1) GO TO 5 IF(ALL2) GO TO 8 CW .TRUE. IN FLAG-SETTING FIRST CALL FROM MAIN. 2 KUT=I GO TO 8 4 EOF=.TRUE. 5 IF(KUT.EQ.0) GO TO 9 GO TO 8 6 IF(IL.EQ.-1) GO TO 8 WRITE(6,500)IS2V(3),IL2V(3),IP2V(3) BACKSPACE 5 8 ALL2=.FALSE. 9 RETURN 500 FORMAT(/' UNIT 5 INPUT: TERMINATOR EXPECTED AFTER IS2,IL2,IP2=', 1 3I2,' NEXT LINE SAVED'/) END C C*************************************************************** C SUBROUTINE RDV(ITAPE7) C C ADAPTED FROM C PROGRAM OF KTT FOR READING DIPOLE MATRIX ELEMENT DATA C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (ICH2=MZCHF*MZCHF) COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/DVECT/DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST) COMMON/SCALE/AZ,AZP,AZ2,AZP2,K10,K20 COMMON/N1N2/N1N2A,N1N2B,N1AM(ICH2),N2AM(ICH2),N1BM(ICH2), 1 N2BM(ICH2) C C READ(ITAPE7) MNP2D2,NCHND2,LRGLD2,NPTYD2,NSPND, 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,K20,K10,M11,MLAST C C READ(ITAPE7)((DVECTL(K2,M1),K2=1,K20),M1=M11,MLAST) READ(ITAPE7)((DVECTV(K2,M1),K2=1,K20),M1=M11,MLAST) C C-----READ ANGULAR COEFFICIENTS READ(ITAPE7) ((AMATL(J,I),J=1,NCHND1),I=1,NCHND2), 1 ((BMATL(J,I),J=1,NCHND1),I=1,NCHND2), 2 ((BMATV(J,I),J=1,NCHND1),I=1,NCHND2) C C FORM INDEXING ARRAYS TO CONDENSE LOOPS OVER NCHND1 AND C NCHND2 SATLISFYING CONDITIONS THATL AMATL(N1,N2) IS C NOT NEGLIGIBLE AND BMATL(N1,N2) IS NOT NEGLIGIBLE. C THESE WILL BE USED IN SUBROUTINE OPDOUT TO ENABLE C VECTORIZATION. C N1N2A=0 N1N2B=0 DO 1 N1=1,NCHND1 DO 1 N2=1,NCHND2 IF (ABS(AMATL(N1,N2)) .GT. TINY) THEN N1N2A=N1N2A+1 N1AM(N1N2A)=N1 N2AM(N1N2A)=N2 END IF IF (ABS(BMATL(N1,N2)) .GT. TINY) THEN N1N2B=N1N2B+1 N1BM(N1N2B)=N1 N2BM(N1N2B)=N2 END IF 1 CONTINUE C IF(NSPND.NE.0)RETURN C READ(ITAPE7,ERR=2,END=3) ((AMATV(J,I),J=1,NCHND1),I=1,NCHND2) RETURN C 2 BACKSPACE(ITAPE7) !ATTEMPT SALVAGE 3 DO J=1,NCHND2 DO I=1,NCHND1 AMATV(I,J)=0 ENDDO ENDDO WRITE(6,620) C RETURN C 620 FORMAT(' NO 2ND SET AMAT AVAILABLE: VELOCITY RESULTS UNCERTAIN') C END C************************************************************************* C SUBROUTINE READXY(ETOT,QETEST,NCHOP,FLAGL,FLAGV,MFIN,IPERT) C C NRB: C READ AND INTERPOLATE UNPHYSICAL XM,YM,ZDL,ZDV MATRICES C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C LOGICAL FLAGL,FLAGV,FLAGV1,FLAGV2,SKIPPY C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (MXBUF=MZCHF*MZEST) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF), X RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/DVECT/BUF(2*MNPEXT*MZEST) C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBZDD/ZDLM(MZCHF,MZEST),ZDVM(MZCHF,MZEST) C DIMENSION XM1(MZCHF,MZCHF),YM1(MZCHF,MZCHF) DIMENSION XM2(MZCHF,MZCHF),YM2(MZCHF,MZCHF) DIMENSION DXL1(MZCHF,MZEST),DYL1(MZCHF,MZEST) DIMENSION DXL2(MZCHF,MZEST),DYL2(MZCHF,MZEST) DIMENSION DXV1(MZCHF,MZEST),DYV1(MZCHF,MZEST) DIMENSION DXV2(MZCHF,MZEST),DYV2(MZCHF,MZEST) C EQUIVALENCE (XM1,F2),(YM1,FP2),(XM2,FPP2),(YM2,DF2) C EQUIVALENCE (XMOLD,DFP2),(YMOLD,DFPP2) EQUIVALENCE (DXL1(1,1),BUF(1)),(DYL1(1,1),BUF(MXBUF+1)) EQUIVALENCE (DXL2(1,1),BUF(2*MXBUF+1)),(DYL2(1,1),BUF(3*MXBUF+1)) EQUIVALENCE (DXV1(1,1),BUF(4*MXBUF+1)),(DYV1(1,1),BUF(5*MXBUF+1)) EQUIVALENCE (DXV2(1,1),BUF(6*MXBUF+1)),(DYV2(1,1),BUF(7*MXBUF+1)) C SAVE E1,E2,MFIN1,MFIN2,FLAGV1,FLAGV2,M11 C C BELOW ARE ALREADY SAVED BY BEING EQUIVALENCED TO VARIABLES IN COMMON, C USE OF SAVE THEN GENERATES ERRORS ON SOME COMPILERS C C SAVE XM1,YM1,XM2,YM2,DXL1,DYL1,DXL2,DYL2,DXV1,DYV1,DXV2,DYV2 C C C N.B. WE ALLOW FOR THE POSSIBILTY OF A CHANGE IN FLAGV AND MFIN, C BUT FLAGL AND M11 DO NOT CHANGE IN THIS LOOP. C NCHOP=NCHND2 HERE. C 5 IF(ETOT.GE.TZERO)THEN C C SEE IF WE NEED TO READ NEW DATA, OR JUST INTERPOLATE EXISTING DATA C IF(ETOT.GT.E1-QETEST.AND.ETOT.LT.E2+QETEST) GO TO 10 C C MOVE OLD DATA C E1=E2 FLAGV1=FLAGV2 MFIN1=MFIN2 DO J=1,NCHOP DO I=J,NCHOP XM1(I,J)=XM2(I,J) YM1(I,J)=YM2(I,J) ENDDO ENDDO IF(FLAGL)THEN DO M1=M11,MFIN1 DO J=1,NCHOP DXL1(J,M1)=DXL2(J,M1) DYL1(J,M1)=DYL2(J,M1) ENDDO ENDDO ENDIF IF(FLAGV1)THEN DO M1=M11,MFIN1 DO J=1,NCHOP DXV1(J,M1)=DXV2(J,M1) DYV1(J,M1)=DYV2(J,M1) ENDDO ENDDO ENDIF ENDIF C C READ NEW DATA C XM,YM (NCHOP SHOULD NEVER CHANGE FOR FIXED SYMMETRIES) C SKIP TO NEXT ENERGY IF UNPHYSICAL DATA COULD NOT BE COMPUTED C 7 SKIPPY=.TRUE. READ(70,END=10)NCHXX,E2,M20 IF(NCHXX.NE.NCHOP)STOP 'READXY ERROR' SKIPPY=IQMSH(M20).GT.0 READ(70)((XM2(I,J),I=J,NCHOP),J=1,NCHOP) READ(70)((YM2(I,J),I=J,NCHOP),J=1,NCHOP) C C ZDL,ZDV C READ(70,END=10)M11,MFIN2 DO M1=M11,MFIN2 READ(70)FLAGL IF(FLAGL)THEN READ(70)(DXL2(J,M1),J=1,NCHOP) READ(70)(DYL2(J,M1),J=1,NCHOP) ENDIF IF(IPRINT.EQ.-1.AND.IPERT.EQ.0)THEN FLAGV2=.FALSE. ELSE READ(70)FLAGV2 IF(FLAGV2)THEN READ(70)(DXV2(J,M1),J=1,NCHOP) READ(70)(DYV2(J,M1),J=1,NCHOP) ENDIF ENDIF ENDDO C IF(SKIPPY)GO TO 7 C C IF(ETOT.LT.TZERO)THEN IF(MZMNP.LT.3*MZCHF)THEN WRITE(6,*)' DIMENSION ERROR: INCREASE MZMNP TO AT LEAST' X ,' 3*MZCHF=',3*MZCHF STOP ENDIF E1=99999. RETURN ENDIF C C SEE IF THIS ENERGY IS ENOUGH C IF(ETOT.GT.E2+QETEST)GO TO 5 C C FORM NEW MATRICES, JUST LINEAR INTERPOLATION C 10 FLAGV=FLAGV1.AND.FLAGV2 MFIN=MIN(MFIN1,MFIN2) IF(SKIPPY)THEN C LAST ENERGY BAD SO T1=ONE T2=TZERO ELSE T1=(E2-ETOT)/(E2-E1) T2=(E1-ETOT)/(E2-E1) ENDIF DO J=1,NCHOP DO I=J,NCHOP XM(I,J)=T1*XM1(I,J)-T2*XM2(I,J) YM(I,J)=T1*YM1(I,J)-T2*YM2(I,J) XM(J,I)=XM(I,J) YM(J,I)=YM(I,J) ENDDO ENDDO IF(FLAGL)THEN DO M1=M11,MFIN DO J=1,NCHOP TR=T1*DXL1(J,M1)-T2*DXL2(J,M1) TI=T1*DYL1(J,M1)-T2*DYL2(J,M1) ZDLM(J,M1)=DCMPLX(TR,TI) ENDDO ENDDO ENDIF IF(FLAGV)THEN DO M1=M11,MFIN DO J=1,NCHOP TR=T1*DXV1(J,M1)-T2*DXV2(J,M1) TI=T1*DYV1(J,M1)-T2*DYV2(J,M1) ZDVM(J,M1)=DCMPLX(TR,TI) ENDDO ENDDO ENDIF C RETURN END C********************************************************************** C SUBROUTINE S1S2(PI,PPI,CCI,PJ,PPJ,PPPJ,CCJ,S1,S2) C C CALCULATES SURFACE TERMS S1 AND S2 IN DIPOLE LENGTH AND ACCELERATION C CONVERSION FORMULAE C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H C FLIJ=(PI*PPJ-PPI*PJ)*RZERO FKIJ=(CCI-CCJ)/(RZERO*RZERO) FIFJ=PI*PJ S1=FLIJ+FIFJ S2=FKIJ*(FLIJ-FIFJ)-2.*(PI*PPPJ-PPI*PPJ) C RETURN END C********************************************************************* C SUBROUTINE STRANS(DLVECT,DVVECT,FLAGL,FLAGV,DLTOT,DVTOT,DL2,DV2 X ,QJUMP,M1) C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C CBL PARAMETER (MWORK=100*MZDEG) PARAMETER (TZERO=0.0) PARAMETER (ZERO=(0.0,0.0)) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C LOGICAL FLAGV,QDT,QJUMP,FLAGL,SKIPPY C C CHANGE TO S-MATRIX NORMALIZATION; RETURN THE SQUARES OF DIPOLE MATRIX C ELEMENTS AND THE SUMMED TOTAL. C NRB: C EVALUATE SINGLE AND MULTI-STATE QDT, DETAILED AND AVERAGED, C MQDT IS INTERPOLATED. ALSO, OPTIONALLY PRECONVOLUTE TOTAL MQDT. C COMMON/FMATCH/P1(MZCHF,MZEST),PP1(MZCHF,MZEST),PPP1(MZCHF,MZEST), 1 F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FPP2(MZCHF,MZCHF), 2 DP1(MZCHF,MZEST),DPP1(MZCHF,MZEST),DPPP1(MZCHF,MZEST), 3 DF2(MZCHF,MZCHF),DFP2(MZCHF,MZCHF),DFPP2(MZCHF,MZCHF) X,RK(MZCHF,MZCHF), 4 XM(MZCHF,MZCHF),YM(MZCHF,MZCHF) COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/FLAGQ/QDT COMMON/TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) C COMMON/NRBKHI/ZKHICC(MZDEG,MZDEG),ZKHIOC(MZCHF,MZDEG),ZVAL(MZDEG) X ,ZDL(MZCHF),ZDV(MZCHF),ZDL0(MZCHF),ZDV0(MZCHF) CBL X,ZWL(MZDEG,MZDEG),ZWR(MZDEG,MZDEG),ZWORK(MWORK),RWORK(2*MZDEG) COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) C DIMENSION DLVECT(MZCHF),DVVECT(MZCHF),DL2(MZCHF),DV2(MZCHF) DIMENSION DLSR(MZCHF),DLSI(MZCHF),DVSR(MZCHF),DVSI(MZCHF) C IF(QJUMP)GO TO 100 C C FIRST TIME THRU, TRANSFORM DIPOLE MATRIX ELEMENT VIA D(KHI)=Y**(T*)D(R) C AND WRITE IT TO SCRATCH FILE, THEN RETURN. HERE (Y)=-(YM)-I(XM). C PUT THE REAL AND IMAGINARY PARTS IN DLSR,DVSR AND DLSI,DVSI VECTORS C IF(IQDT.GT.0)WRITE(70)FLAGL IF(FLAGL)THEN DO J=1,NCHOP DLSR(J)=TZERO DLSI(J)=TZERO ENDDO DO J=1,NCHOP DO I=1,NCHOP DLSR(J)=DLSR(J)-DLVECT(I)*YM(I,J) DLSI(J)=DLSI(J)+DLVECT(I)*XM(I,J) ENDDO ENDDO IF(IQDT.GT.0)THEN WRITE(70)(DLSR(J),J=1,NCHOP) WRITE(70)(DLSI(J),J=1,NCHOP) ENDIF ENDIF IF(IQDT.GT.0)WRITE(70)FLAGV IF(FLAGV) THEN DO J=1,NCHOP DVSR(J)=TZERO DVSI(J)=TZERO ENDDO DO J=1,NCHOP DO I=1,NCHOP DVSR(J)=DVSR(J)-DVVECT(I)*YM(I,J) DVSI(J)=DVSI(J)+DVVECT(I)*XM(I,J) ENDDO ENDDO IF(IQDT.GT.0)THEN WRITE(70)(DVSR(J),J=1,NCHOP) WRITE(70)(DVSI(J),J=1,NCHOP) ENDIF ENDIF C IF(IQDT.GT.0)RETURN C C STORE DIPOLE MATRICES FOR BETA (KAB, MARCH 1991) C NRB: ARE THE SIGNS CORRECT SINCE ORIGINALS ARE TRANSFORMED BY Y NOT Y**(T*)? C THE ORIGINAL SIGNS ARE KEPT FOR DSTOR HENCE THE FIRST SIGN SWAP, C THEN RETURNED FOR MQDT MODE. C IF(IPRINT.EQ.3) THEN DO J=1,NCHOP !TO GET ORGINAL CODE (Y) DLSR(J)=-DLSR(J) DVSR(J)=-DVSR(J) ENDDO CALL DSTOR(2,DLSR,DLSI,DVSR,DVSI,NCHOP) DO J=1,NCHOP !TO GET BACK TO Y**(T*) DLSR(J)=-DLSR(J) DVSR(J)=-DVSR(J) ENDDO ENDIF C C COMPUTE MODULUS SQUARES OF DIPOLE MATRIX ELEMENTS AND SUMMED TOTAL C DLTOT=TZERO IF(FLAGL)THEN DO J=1,NCHOP DL2(J)=DLSR(J)**2+DLSI(J)**2 DLTOT=DLTOT+DL2(J) ENDDO ENDIF DVTOT=TZERO IF(FLAGV) THEN DO J=1,NCHOP DV2(J)=DVSR(J)**2+DVSI(J)**2 DVTOT=DVTOT+DV2(J) ENDDO ENDIF C 100 IF((ABS(IPRINT).NE.2.OR..NOT.QDT).AND.IQDT.EQ.0) RETURN C IF(EWIDTH.GE.TZERO.AND.IQDT.GT.0)THEN C C TOTAL, PRECONVOLUTED C CALL PMQDT(DLVECT,DVVECT,FLAGL,FLAGV,DL2,DV2,M1) C DLTOT=DL2(1) DVTOT=DV2(1) RETURN ENDIF C IF(IQDT.GT.0) THEN C C TREAT ALL CLOSED-CHANNELS AS OPEN C CALL DMQDT(DLVECT,DVVECT,FLAGL,FLAGV,DL2,DV2,M1) C IF(.NOT.QJUMP)RETURN C C COMPUTE NEW TOTALS IF A CONTRACTION HAS TAKEN PLACE C IF(QDT)THEN DLTOT=TZERO IF(FLAGL)THEN DO J=1,NCHOP DLTOT=DLTOT+DL2(J) ENDDO ENDIF DVTOT=TZERO IF(FLAGV) THEN DO J=1,NCHOP DVTOT=DVTOT+DV2(J) ENDDO ENDIF ENDIF C ENDIF C IF(IQDT.LE.0.OR.QDT)THEN C C TREAT ONLY THOSE ATTACHED TO LOWEST TARGET STATE AS CLOSED. C CALL DSQDT(DLVECT,DVVECT,FLAGL,FLAGV,DL2,DV2,QJUMP) C ENDIF C C C CHECK THAT THE TOTAL AGREES WITH WHAT WAS OBTAINED BEFORE C CASE QDT=.TRUE. AND NO RADIATIVE DECAYS, OR C GENERATE NEW TOTALS CASE RADIATIVE DECAYS AND/OR QDT=.FALSE. C TC=AC IF(IQDT.GT.0)TC=10.*AC IF(FLAGL)THEN CHKTOT=TZERO DO J=1,NCHOP CHKTOT=CHKTOT+DL2(J) ENDDO IF(IRDEC.NE.0.OR..NOT.QDT)THEN DLTOT=CHKTOT ELSE IF(ABS(DLTOT-CHKTOT).GT.TC)WRITE(6,602)'LENGTH',DLTOT,CHKTOT ENDIF ENDIF IF(FLAGV) THEN CHKTOT=TZERO DO J=1,NCHOP CHKTOT=CHKTOT+DV2(J) ENDDO IF(IRDEC.NE.0.OR..NOT.QDT)THEN DVTOT=CHKTOT ELSE IF(ABS(DVTOT-CHKTOT).GT.TC)WRITE(6,602)'VELOCITY',DVTOT,CHKTOT ENDIF ENDIF C RETURN C 602 FORMAT(' ***WARNING: QDT VALUES OF DIPOLE ',A,' MATRIX ELEMENT', 1 ' SQUARED DO NOT AGREE'/' SUMMATION OVER OPEN PLUS WC CHANNELS', 2 ' GIVES',E17.8/' SUMMATION OVER S-MATRIX OPEN CHANNELS GIVES', 3 E17.8) END 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 ADAPTED FROM STGF. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (NCHF1=MZCHF+1) C COMMON/CTHET/BB(NCHF1,MZTET),BG(NCHF1,MZTET),MSUM(NCHF1) 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 = ((DLR+R*F2)*S+CX)*E C RETURN END C C********************************************************** C SUBROUTINE TARG2(NCHF2) C C THIS SUBROUTINE SORTS OUT THE NUMBER OF COUPLED CHANNELS ATTACHED C TO EACH TARGET STATE AND STORES THE RESULTS IN ARRAY NCONAT C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (ROUND=1.E-5) COMMON/SCALE/AZ,AZP,AZ2,AZP2,K10,K20 COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON /TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) KUT=1 DO 1 N=1,NAST NCONAT(N)=0 ETARG=ENAT(N)*AZP2 DO 2 I=KUT,NCHF2 IF(ABS(ETARG-ECH2(I)).GT.ROUND) THEN KUT=I GOTO 1 ELSE IF(ABS(ETARG-ECH2(I)).GT.ROUND*1.E-3) X WRITE(6,*)'****SR.TARG2: WARNING, POSSIBLE MIS-MATCH IN' X,' DECODING CHANNEL TARGETS -',I,N ENDIF NCONAT(N)=NCONAT(N)+1 NTARG(I)=N 2 CONTINUE KUT=NCHF2+1 1 CONTINUE 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 PARAMETER (MZMSH0=MZMSH) C PARAMETER (NCHF1=MZCHF+1) C COMMON/CHAN/ECH(NCHF1),EPS(NCHF1),FKNU(NCHF1),CC(NCHF1), 1 RINF(NCHF1),LLCH(NCHF1),ITARG(NCHF1),NCHF,NCHOP,NCHOP1 COMMON/CEN/ETOT,MXE,NWT,NZ COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/CTHET/BB(NCHF1,MZTET),BG(NCHF1,MZTET),MSUM(NCHF1) COMMON/NRBZED/TZED LOGICAL FLAG 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 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 ADAPTED FROM STGF. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (NCHF1=MZCHF+1) C COMMON/CTHET/BB(NCHF1,MZTET),BG(NCHF1,MZTET),MSUM(NCHF1) COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H 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 REAL*8 FUNCTION UMAT(I,N1,N2,M1) C C CALCULATES U INTEGRAL, INTEGRAL OF S(N2)*R**I*P(N1) C REMEMBER THAT UMAT CAN NEVER EXIST DUE TO THE DELTA-FUNCTION, BUT THE C CONTRIBUTION FROM RTWO TO INFINITY IS CALCULATED BY A PREVIOUS CALL C TO VMAT C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/FLAGIN/FLG1(MZCHF,MZEST),FLG2(MZCHF),EXIST1(MZCHF,MZEST), 1 EXIST2(MZCHF),DROP COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON /VUSTOR/ VSTORE(-3:1,MZCHF,MZCHF),USTORE(-3:1,MZCHF,MZCHF), * X1(MZCHF,MZEST),X2(MZCHF),EXISTV(-3:1,MZCHF,MZCHF) LOGICAL FLG1,FLG2,EXIST1,EXIST2,EXISTV,DROP C UMAT=0. MA=M1 IF(MZSA1.EQ.1)MA=1 DO 1 K=1,KPMM(M1) 1 UMAT=UMAT+FSD1(K,N1,MA)*WBODE(K,I)*S(K,N2) IF(FLG1(N1,M1).AND.FLG2(N2)) UMAT=UMAT+USTORE(I,N1,N2) RETURN END C**************************************************************** C SUBROUTINE VERT(V,LV,N,W) C C ________________________________________________________ C | | C | INVERT A GENERAL MATRIX | C | | C | INPUT: | C | | C | V --ARRAY CONTAINING MATRIX | C | | C | LV --LEADING (ROW) DIMENSION OF ARRAY V | C | | C | N --DIMENSION OF MATRIX STORED IN ARRAY V | C | | C | W --INTEGER WORK ARRAY WITH AT LEAST N-1 | C | ELEMENTS | C | | C | OUTPUT: | C | | C | V --INVERSE | C | | C | BUILTIN FUNCTIONS: ABS | C |________________________________________________________| C REAL*8 V(LV,1),S,T INTEGER W(1),I,J,K,L,M,N,P IF ( N .EQ. 1 ) GOTO 110 L = 0 M = 1 10 IF ( L .EQ. N ) GOTO 90 K = L L = M M = M + 1 C --------------------------------------- C |*** FIND PIVOT AND START ROW SWAP ***| C --------------------------------------- P = L IF ( M .GT. N ) GOTO 30 S = ABS(V(L,L)) DO 20 I = M,N T = ABS(V(I,L)) IF ( T .LE. S ) GOTO 20 P = I S = T 20 CONTINUE W(L) = P 30 S = V(P,L) V(P,L) = V(L,L) IF ( S .EQ. 0. ) GOTO 120 C ----------------------------- C |*** COMPUTE MULTIPLIERS ***| C ----------------------------- V(L,L) = -1. S = 1./S DO 40 I = 1,N 40 V(I,L) = -S*V(I,L) J = L 50 J = J + 1 IF ( J .GT. N ) J = 1 IF ( J .EQ. L ) GOTO 10 T = V(P,J) V(P,J) = V(L,J) V(L,J) = T IF ( T .EQ. 0. ) GOTO 50 C ------------------------------ C |*** ELIMINATE BY COLUMNS ***| C ------------------------------ IF ( K .EQ. 0 ) GOTO 70 DO 60 I = 1,K 60 V(I,J) = V(I,J) + T*V(I,L) 70 V(L,J) = S*T IF ( M .GT. N ) GOTO 50 DO 80 I = M,N 80 V(I,J) = V(I,J) + T*V(I,L) GOTO 50 C ----------------------- C |*** PIVOT COLUMNS ***| C ----------------------- 90 L = W(K) DO 100 I = 1,N T = V(I,L) V(I,L) = V(I,K) 100 V(I,K) = T K = K - 1 IF ( K .GT. 0 ) GOTO 90 RETURN 110 IF ( V(1,1) .EQ. 0. ) GOTO 120 V(1,1) = 1./V(1,1) RETURN 120 WRITE(6,*) 'ERROR: MATRIX HAS NO INVERSE' STOP END C********************************************************************* C REAL*8 FUNCTION VMAT(I,N1,N2,M1) C C CALCULATES V INTEGRAL, INTEGRAL OF C(N2)*R**I*P(N1) OR C THETA(N2)*R**I*P(N1) C IMPLICIT REAL*8 (A-H,O-Y) C IMPLICIT COMPLEX*16 (Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TZERO=0.0) PARAMETER (HALF=0.5) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (FOUR=4.0) C PARAMETER (NCHF1=MZCHF+1) PARAMETER (MNPEXT=MZMNP+MZCHF) C LOGICAL FLG1,FLG2,EXIST1,EXIST2,EXISTV,DROP,QDT,BQDT,SKIPPY C COMMON /CDKKP/ AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF) X,BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), * MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/CTRLBF/IPRINT,IBUT,TINY,DELTA,AC COMMON/CHAN1/ECH1(MZCHF),EPS1(MZCHF,MZEST),FNU1(MZCHF,MZEST) X,CC1(MZCHF) COMMON/CHAN2/ECH2(MZCHF),EPS2(MZCHF),FNU2(MZCHF),CC2(MZCHF) COMMON/CHAN/ECH(NCHF1),EPS(NCHF1),FKNU(NCHF1),CC(NCHF1), 1 RINF(NCHF1),LLCH(NCHF1),ITARG(NCHF1),NCHF,NCHOP0,NCHOP1 COMMON/FLAGIN/FLG1(MZCHF,MZEST),FLG2(MZCHF),EXIST1(MZCHF,MZEST), 1 EXIST2(MZCHF),DROP COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON /VUSTOR/ VSTORE(-3:1,MZCHF,MZCHF),USTORE(-3:1,MZCHF,MZCHF), * X1(MZCHF,MZEST),X2(MZCHF),EXISTV(-3:1,MZCHF,MZCHF) COMMON/CACC/AC0,ACNUM,ACJWBK,ACZP,LACC COMMON/CTHET/BB(NCHF1,MZTET),BG(NCHF1,MZTET),MSUM(NCHF1) COMMON/CJWBK/D(15) COMMON/FLAGQ/QDT COMMON/CVMAT/TI(MZCHF,MZEST),TIP(MZCHF,MZEST),TID(MZCHF,MZEST), + TIDP(MZCHF,MZEST),PIN(MZCHF,MZEST), + TF(MZCHF),TFP(MZCHF),TFD(MZCHF),TFDP(MZCHF),PFN(MZCHF), * BB2(MZCHF,MZTET),BB1(MZCHF,MZTET,MZEST),D2(15,MZCHF), * MSUM1(MZCHF,MZEST),MSUM2(MZCHF) C COMMON/NRBQDT/EWIDTH,IQDT,IOMIT(MZCHF),SKIPPY,IQMSH(MZMSH0) COMMON/NRBFLG/TOL0V,IFLAGV COMMON/NRBTMP/E2 C DATA NLAG/6/ C C C C C C CHECK IF INTEGRAL HAS BEEN CALCULATED C IF(EXISTV(I,N1,N2)) THEN VMAT=VSTORE(I,N1,N2) RETURN ENDIF C C RZERO TO RTWO C VMAT=TZERO MA=M1 IF(MZSA1.EQ.1)MA=1 DO K=1,KPMM(M1) VMAT=VMAT+FSD1(K,N1,MA)*WBODE(K,I)*C(K,N2) ENDDO C C SEE IF FLAG HAS BEEN SET TO SKIP RTWO TO INFINITY C IF(.NOT.(FLG1(N1,M1).AND.FLG2(N2))) GOTO 2 IF(IFLAGV.EQ.0)GO TO 2 C C IF N2 IS WC CHANNEL CHECK CONVERGENCE OF INTEGRAL. C DELTAN=HALF BQDT=QDT.OR.IQDT.NE.0 IF(BQDT.AND.N2.LE.NCHOP.AND.EPS2(N2).LT.TZERO) THEN IF(FNU2(N2).LT.FNU1(N1,M1)+DELTAN) THEN IF(IQDT.LE.0)THEN WRITE(6,600) E2,N1,N2,FNU1(N1,M1),FNU2(N2) CALL ABORTT(14) ENDIF c call abortt(14) !temp SKIPPY=.TRUE. VMAT=TZERO RETURN ENDIF ENDIF C C RTWO TO INFINITY C NCHOP0=NCHOP C C SET UP INITIAL CHANNEL IF NECESSARY C IF(EXIST1(N1,M1)) THEN MS=MSUM1(N1,M1) MSUM(N1)=MS DO K=1,MS BB(N1,K)=BB1(N1,K,M1) ENDDO ELSE EPS(N1)=EPS1(N1,M1) FKNU(N1)=FNU1(N1,M1) LLCH(N1)=(SQRT(ONE+FOUR*CC1(N1))-ONE+DELTA)*HALF CALL THETA(RTWO,N1,T,TP,TD,TDP) IF(ABS(T).LT.TOL0V)THEN X1(N1,M1)=TZERO !NRB: FUNCTION & INTEGRAL ZERO ELSE MS=MSUM(N1) MSUM1(N1,M1)=MS DO K=1,MS BB1(N1,K,M1)=BB(N1,K) ENDDO X1(N1,M1)=FSD1(KP2,N1,MA)/T TI(N1,M1)=T TIP(N1,M1)=TP TID(N1,M1)=TD TIDP(N1,M1)=TDP PIN(N1,M1)=(RTWO**(HALF*FKNU(N1)))*EXP(-HALF**RTWO X /FKNU(N1)) ENDIF EXIST1(N1,M1)=.TRUE. ENDIF C C SET UP FINAL CHANNEL C N21=NCHND1+1 IF(EXIST2(N2)) THEN IF(N2.LE.NCHOP) THEN DO K=1,15 D(K)=D2(K,N2) ENDDO ELSE MS=MSUM2(N2) MSUM(N21)=MS DO K=1,MS BB(N21,K)=BB2(N2,K) ENDDO ENDIF ELSE EP2=EPS2(N2) EPS(N21)=EP2 L=(SQRT(ONE+FOUR*CC2(N2))-ONE+DELTA)*HALF LLCH(N21)=L IF(N2.LE.NCHOP) THEN FKNU(N21)=SQRT(EP2) CALL INJWBK(EP2,L) DO K=1,15 D2(K,N2)=D(K) ENDDO ELSE FKNU(N21)=FNU2(N2) CALL THETA(RTWO,N21,T,TP,TD,TDP) IF(ABS(T).LT.TOL0V)THEN X2(N2)=TZERO !NRB: FUNCTION & INTEGRAL ZERO ELSE X2(N2)=C(KP2,N2)/T MS=MSUM(N21) MSUM2(N2)=MS DO K=1,MS BB2(N2,K)=BB(N21,K) ENDDO TF(N2)=T TFP(N2)=TP TFD(N2)=TD TFDP(N2)=TDP PFN(N2)=(RTWO**(HALF*FNU2(N2)))*EXP(-HALF*RTWO/FNU2(N2)) ENDIF ENDIF EXIST2(N2)=.TRUE. ENDIF C XN1=X1(N1,M1) XN2=X2(N2) IF(ABS(XN1).LT.TOL0V.OR.ABS(XN2).LT.TOL0V)GO TO 2 !NRB: FUNCTION & INTEGRAL ZERO IF(N2.LE.NCHOP) THEN FKNU(N1)=FNU1(N1,M1) FKNU(N21)=SQRT(EPS2(N2)) CALL ZSLAG(N21,N1,NLAG,-I,ZS3,ZSD3) VMAT=VMAT+DBLE(ZS3)*XN1 C STORE UMAT FOR RTWO TO INFINITY, SPIN OFF FROM ZSLAG USTORE(I,N1,N2)=DIMAG(ZS3)*XN1 ELSE IF(I.NE.0)THEN CALL TLAG(N1,N21,NLAG,-I,T1,T2,T3) ELSE T1=HALF*((PIN(N1,M1)/PFN(N2))* X (TI(N1,M1)*TIDP(N1,M1)-TID(N1,M1)*TIP(N1,M1))+ X (PFN(N2)/PIN(N1,M1))* X (TF(N2)*TFDP(N2)-TFD(N2)*TFP(N2))) ENDIF VMAT=VMAT+T1*XN1*XN2 ENDIF C 2 VSTORE(I,N1,N2)=VMAT EXISTV(I,N1,N2)=.TRUE. C RETURN 600 FORMAT(//10X,'*** OUTER REGION INTEGRAL WILL NOT CONVERGE FOR' X ,' E2=',E14.8,' ***'// X 26X,'BOUND',10X,'FREE'/26X,5('-'),10X,4('-')/ X ' CHANNEL',18X,I3,12X,I3/ X ' EFFECTIVE QUANTUM NUMBER',F10.4,5X,F10.4/) END C********************************************************************** C SUBROUTINE WARCHV(IPRINT,MLAST,IS1,IL1,IP1,M1,E1,TOTAL,EMESH,MXE2 X ,IFLAGBP) C C WRITES A BLOCK IN THE ARCHIVE FILE XSECTN ON UNIT 7 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C PARAMETER (TINNY=1.0E-10) C COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON/SCALE/AZ,AZP,AZ2,AZP2,K10,K20 COMMON /TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) C DIMENSION SIGMAL(MZTAR),SIGMAV(MZTAR), TOTAL(MZSA2),TOTM1(MZEST) DIMENSION EMESH(MXE2) C LOGICAL QDT C WRITE(7,700)NZED,NELC,NAST IF(IFLAGBP.EQ.0)WRITE(7,*)(ISAT(I),LAT(I),I=1,NAST) IF(IFLAGBP.NE.0)WRITE(7,*)(0,(LAT(I)+1)*(-1)**ISAT(I),I=1,NAST) WRITE(7,701)(ENAT(I),I=1,NAST) C MXENOZ=MXE2 ISAVE=MZSA2 C C SKIP ZEROS IF ARCHIVE IS REQUIRED C IF(IPRINT.NE.-2) THEN IF(ISAVE .NE. 1)THEN DO 10 M2=1,MXE2 TOT=TOTAL(M2) IF(TOT .GT. TINNY) GO TO 4 MXENOZ=MXENOZ-1 10 CONTINUE ELSE DO 20 M2=1,MXE2 READ(8,REC=M2)(TOTM1(I),I=1,MLAST) IF(TOTM1(M1).GT.TINNY) GOTO 4 MXENOZ=MXENOZ-1 20 CONTINUE ENDIF END IF 4 WRITE(7,700)IS1,IL1,IP1,M1 WRITE(7,702)AZ2*E1,MXENOZ WRITE(7,703)0. C C PARTIAL CROSS SECTION FOR EACH TARGET STATE C IF(IPRINT.EQ.-2) THEN IREC0=MXE2*(M1-1) C DO 1 M2=1,MXE2 IREC=IREC0+M2 E2=EMESH(M2) C READ(8,REC=IREC)QDT,(SIGMAL(I),SIGMAV(I),I=1,NAST) C C PRINT CROSS SECTIONS ASSOCIATED WITH THE FIRST NTARO STATES C DO 5 I=NAST,1,-1 IF(SIGMAL(I)+SIGMAV(I).LE.TINNY) GO TO 5 NTARO=I GOTO 6 5 CONTINUE WRITE(6,605)AZ2*E2,TINNY NTARO=0 6 WRITE(7,800)AZ2*E2,NTARO,QDT DO I=1,NTARO WRITE(7,801)SIGMAL(I),SIGMAV(I) ENDDO 1 CONTINUE WRITE(6,603)IS1,IL1,IP1,M1 C C OR TOTAL CROSS SECTION C ELSE IF (ISAVE .NE. 1) THEN DO M2=MXE2-MXENOZ+1,MXE2 E2=EMESH(M2) DE=AZ2*(E2-E1) TOT=TOTAL(M2) WRITE(7,704)DE,TOT ENDDO WRITE(6,604)IS1,IL1,IP1,M1 ELSE DO M2=MXE2-MXENOZ+1,MXE2 E2=EMESH(M2) DE=AZ2*(E2-E1) READ(8,REC=M2)(TOTM1(I),I=1,MLAST) WRITE(7,704)DE,TOTM1(M1) ENDDO WRITE(6,604)IS1,IL1,IP1,M1 ENDIF RETURN C 700 FORMAT(4I5) 701 FORMAT(5E14.6) 702 FORMAT(E14.6,I5) 703 FORMAT(F12.6) 704 FORMAT(1P,E14.6,E10.3) 603 FORMAT(/' CROSS SECTIONS TO TARGET STATES WRITTEN FOR STATE', * 3I2,' LEVEL',I2,' IN ARCHIVE FILE'/1X,79('+')/) 604 FORMAT(/' TOTAL CROSS SECTION DATA FOR STATE ',3I3, * ' LEVEL',I3,' WRITTEN TO ARCHIVE FILE'/1X,79('+')/) 605 FORMAT(' FOR FINAL ENERGY =',F13.6, * ' NO CROSS SECTIONS GREATER THAN',1P,E9.1) 800 FORMAT(F12.6,I3,L2) 801 FORMAT(1P,E12.5,E13.5) END C**************************************************** REAL*8 FUNCTION WRACAH(J1,J2,L2,L1,J3,L3,FCT,MFD) C OUT REAL FUNCTION SJS(J1,J2,J3,L1,L2,L3,FCT,MFD) C C NORDITA1966(F.P.L.4) WITH MODIFICATIONS UCL'72MAR20,DL'82JUL28W.E. C THE SIX QUANTUM NUMBER ARGUMENTS HAVE TWICE THEIR PHYSICAL VALUE; C FACTORIALS MUST BE SUPPLIED AS FCT(I)=(I/2-1)!/16**(I/2-1),I=4,M,2 C (FCT(2)=0!=1), AND PHASE FACTORS AS FCT(I)=MOD(I+1,4)-1,I=1,MFD,2. C DOUBLE-PRECISION FCT, OMEGA: IMPLICIT REAL*8 (F,O,W) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C DIMENSION FCT(MFD) OMEGA = 0.0 IF (J1+J2.LT.J3) GO TO 700 IF (J1+L2.LT.L3) GO TO 700 IF (J2+L3.LT.L1) GO TO 700 IF (J3+L1.LT.L2) GO TO 700 IF (ABS(J1-J2).GT.J3) GO TO 700 IF (ABS(J1-L2).GT.L3) GO TO 700 IF (ABS(J2-L3).GT.L1) GO TO 700 IF (ABS(J3-L1).GT.L2) GO TO 700 IJ0 = J1 + J2 + J3 + 2 IJ1 = J1 + L2 + L3 + 2 IJ2 = L1 + J2 + L3 + 2 IJ3 = L1 + L2 + J3 + 2 C IF (MOD(IJ0,2)+MOD(IJ1,2)+MOD(IJ2,2)+MOD(IJ3,2).NE.0) GO TO 700 IWMIN = MAX(IJ0,IJ1,IJ2,IJ3)+2 ID1 = IJ0 + IJ1 - J1-J1 + 2 ID2 = IJ0 + IJ2 - J2-J2 + 2 ID3 = IJ0 + IJ3 - J3-J3 + 2 IWMAX = MIN(ID1,ID2,ID3)-2 IF (IWMAX.LT.IWMIN) GO TO 700 IF (IWMAX.GT.MFD) GO TO 702 DO 701 IW=IWMIN,IWMAX,2 701 OMEGA = -FCT(IW-1)*FCT(IW)/(FCT(ID1-IW)*FCT(ID2-IW)*FCT(ID3-IW)* * FCT(IW-IJ0)*FCT(IW-IJ1)*FCT(IW-IJ2)*FCT(IW-IJ3)) + OMEGA IJ0 = IJ0+2 IJ1 = IJ1+2 IJ2 = IJ2+2 IJ3 = IJ3+2 OMEGA = OMEGA * SQRT( C (FCT(ID1-IJ0)*FCT(ID2-IJ0)*FCT(ID3-IJ0)/FCT(IJ0))* C (FCT(ID1-IJ1)*FCT(ID2-IJ1)*FCT(ID3-IJ1)/FCT(IJ1))* C (FCT(ID1-IJ2)*FCT(ID2-IJ2)*FCT(ID3-IJ2)/FCT(IJ2))* C (FCT(ID1-IJ3)*FCT(ID2-IJ3)*FCT(ID3-IJ3)/FCT(IJ3)))/16 C PHASE FACTOR BETWEEN SJS AND W: IF(MOD(J1+J2+L1+L2,4).NE.0) OMEGA=-OMEGA 700 WRACAH = OMEGA RETURN 702 WRITE(6,'("0FCT.SJS/WRACAH: FACTORIAL ARRAY TOO SHORT")') STOP END C********************************************************************** C SUBROUTINE WSCRA(IPRINT,M,M2,FIRST,SIG,OMEG,DLTOT,DL2,DV2 X ,FLAGL,FLAGV,MXE2) C C WRITES OR UPDATES THE RECORD M2 IN THE SCRATCH FILE ON UNIT 8 C NOT USED WHEN IPRINT=-3, RECOMBINATION PRODUCTION. OMEGPR STORED C INTERNALLY, NO SIGMAL OR SIGMAV. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C COMMON/MISC/WBODE(MZPTS,-3:1),EE1(MZEST),BSTO,IPERT COMMON /TARG/ARDEC(MZTAR),ENAT(MZTAR),NCONAT(MZTAR),NZED,NELC,NAST X,IRDEC,ISAT(MZTAR),LAT(MZTAR),NTARG(MZCHF) COMMON/PSC/FSD1(MZPTS,MZCHF,MZSA1),S(MZPTS,MZCHF),C(MZPTS,MZCHF) X,NCHOP COMMON/FLAGQ/QDT COMMON /NRBDR/OMEGPR(MZMET,MZMSH),NDRMET,IGAUGE C DIMENSION DL2(MZCHF),DV2(MZCHF), TOT(MZEST),DLTOT(1,MZEST) DIMENSION SIGMAL(MZTAR),SIGMAV(MZTAR),PARTL(MZTAR),PARTV(MZTAR) C LOGICAL FIRST,FLAGV,QDT,FLAGL C C EITHER CREATE/UPDATE SUM FILE FOR ARCHIVE C IF(IPRINT.GT.-2) THEN C IF(FIRST) THEN WRITE(8,REC=M2)(DLTOT(1,M1),M1=1,M) ELSE READ(8,REC=M2)(TOT(M1),M1=1,M) DO M1=1,M TOT(M1)=TOT(M1)+DLTOT(1,M1) ENDDO WRITE(8,REC=M2)(TOT(M1),M1=1,M) ENDIF C C OR CROSS SECTIONS FOR LEAVING THE RESIDUAL ION IN EACH STATE C TAKING INTO ACCOUNT THE CASE OF IPRINT=-1986: QDT + DUMMY LINE C ELSE IF(IPRINT.EQ.-1986) THEN C IREC=MXE2*(M-1)+M2 WRITE(8,REC=IREC)QDT,(0.0,I=1,2*NAST) C ELSE C IREC=MXE2*(M-1)+M2 K=1 NASTY=NAST IF(IPRINT.EQ.-3)NASTY=NDRMET DO 28 I=1,NASTY SIGMAL(I)=0. SIGMAV(I)=0. IF(K.GT.NCHOP) GOTO 28 IF(NCONAT(I).EQ.0) GOTO 28 DO J=1,NCONAT(I) IF(FLAGL) SIGMAL(I)=SIGMAL(I)+DL2(K) IF(FLAGV) SIGMAV(I)=SIGMAV(I)+DV2(K) K=K+1 ENDDO IF(IGAUGE.LE.0.AND.I.LE.NDRMET)OMEGPR(I,M2)=OMEGPR(I,M2)+ X OMEG*SIGMAL(I) IF(IGAUGE.GT.0.AND.I.LE.NDRMET)OMEGPR(I,M2)=OMEGPR(I,M2)+ X OMEG*SIGMAV(I) SIGMAL(I)=SIG*SIGMAL(I) SIGMAV(I)=SIG*SIGMAV(I) 28 CONTINUE C IF(IPRINT.EQ.-3)RETURN C C AND CREATE/UPDATE RECORD FOR IT C IF(FIRST) THEN WRITE(8,REC=IREC)QDT,(SIGMAL(I),SIGMAV(I),I=1,NAST) ELSE READ(8,REC=IREC)QDT,(PARTL(I),PARTV(I),I=1,NAST) DO I=1,NAST PARTL(I)=PARTL(I)+SIGMAL(I) PARTV(I)=PARTV(I)+SIGMAV(I) ENDDO WRITE(8,REC=IREC)QDT,(PARTL(I),PARTV(I),I=1,NAST) ENDIF C ENDIF RETURN END C****************************************************************** C SUBROUTINE ZEIGEN(A,VAL,N,DELTA) C C DIAGONALISATION OF COMPLEX SYMMETRIC MATRIX C USING METHOD DESCRIBED BY M.J.SEATON, COMPUTER JOURNAL, VOL.12, C PAGE 156, 1969. C IMPLICIT REAL*8 (A-H,O-Z) C COMPLEX*16 A,X,VAL,C,S,H,P,Q,CIM,ZERO,ZONE c logical bflag C INCLUDE 'PARAM' C PARAMETER (ZERO=(0.0,0.0)) PARAMETER (ZONE=(1.0,0.0)) PARAMETER (CIM=(0.0,1.0)) PARAMETER (TZERO=0.0) PARAMETER (QUART=0.25) PARAMETER (HALF=0.5) PARAMETER (TWO=2.0) C DIMENSION A(MZDEG,MZDEG),X(MZDEG,MZDEG),VAL(MZDEG) C DATA IROTM/20/ C SQ(H)=H* CONJG(H) C C TEST FOR QUICK RETURN C IF(N.LE.0)RETURN C IF(N.EQ.1)THEN VAL(1)=A(1,1) A(1,1)=ZONE RETURN ENDIF C IF(N.EQ.2)THEN Q=A(1,2) P = (A(1,1)-A(2,2))*HALF FL= ABS(P*P+Q*Q) V=-QUART*LOG(SQ(P-CIM*Q)/FL) T= CONJG(P)*Q+CONJG(Q)*P U=-QUART * ATAN2(T/FL,(SQ(P)-SQ(Q))/FL) S= SIN(U)*COSH(V) +CIM*COS(U)*SINH(V) C= COS(U)*COSH(V) -CIM*SIN(U)*SINH(V) H=(S*P+C*Q)*S*TWO VAL(1)=A(1,1)-H VAL(2)=A(2,2)+H A(1,1)=C A(1,2)=S A(2,1)=-S A(2,2)=C RETURN ENDIF C C GENERAL CASE C DELTA2=DELTA*DELTA NM1=N-1 C C INITIAL D1,D2 AND X C D1=TZERO H=ZERO DO I=1,N X(I,I)=ZONE D1=D1+SQ(A(I,I)) H=H+A(I,I) ENDDO D1=D1-SQ(H)/DBLE(N) C D2=TZERO DO I=1,NM1 IP1=I+1 DO J=IP1,N X(I,J)=ZERO X(J,I)=ZERO D2=D2+SQ(A(I,J)) ENDDO ENDDO c c identity test c bflag=.true. if(abs(d1).lt.delta2.and.d2.lt.delta2)then c write(6,*)d1,d2 if(d1.eq.tzero)return !must return bflag=.false. endif C D1 = (N-1)*(D1*HALF+D2)*DELTA2/DBLE((N+1)) C C BEGIN ROTATIONS C DO 1010 IROT=1,IROTM DO 1000 IP=1,NM1 IPP1=IP+1 DO 1000 IQ=IPP1,N C C ROTATION CONSTANTS C Q=A(IP,IQ) P = (A(IP,IP)-A(IQ,IQ))*HALF FL = ABS(P*P+Q*Q) BETA = LOG(SQ(P-CIM*Q)/FL)*HALF T=( CONJG(P)*Q+P* CONJG(Q))/FL D=(SQ(P)-SQ(Q))/FL U=-QUART*ATAN2(T,D) C T=TZERO D=TZERO DO I=1,N IF(I.NE.IP.AND.I.NE.IQ)THEN D=D+SQ(A(IP,I))+SQ(A(IQ,I)) T=T+SQ(A(IP,I)+CIM*A(IQ,I)) ENDIF ENDDO T=D-T FN= SQRT(D*D-T*T) GAMMA=LOG((D+T)/FN) C C ITERATION FOR V C V0=-HALF*(BETA+GAMMA) DO ITERV=1,100 V=V0-(FL*SINH((V0+BETA)*TWO)+FN*SINH(V0+GAMMA))/ * (TWO*FL*COSH((V0+BETA)*TWO)+FN*COSH(V0+GAMMA)) IF(ABS(V-V0).LT.DELTA*0.1) GO TO 580 V0=V ENDDO WRITE(6,5010) DELTA C C NEW A,X AND D2 C 580 V=HALF*V S= SIN(U)*COSH(V) +CIM*COS(U)*SINH(V) C= COS(U)*COSH(V) -CIM*SIN(U)*SINH(V) H = (S*P+C*Q)*S*TWO A(IP,IP)=A(IP,IP)-H A(IQ,IQ)=A(IQ,IQ)+H A(IP,IQ)=(C*C-S*S)*Q+TWO*C*S*P D2=D2-SQ(A(IQ,IP))+SQ(A(IP,IQ)) A(IQ,IP)=A(IP,IQ) C DO I=1,N H=X(I,IP) X(I,IP)=H*C-X(I,IQ)*S X(I,IQ)=H*S+X(I,IQ)*C IF(I.NE.IP.AND.I.NE.IQ)THEN H=A(IP,I) A(IP,I)=C*H-A(IQ,I)*S A(IQ,I)=S*H+A(IQ,I)*C D2=D2-SQ(A(I,IP))-SQ(A(I,IQ))+SQ(A(IP,I))+SQ(A(IQ,I)) A(I,IP)=A(IP,I) A(I,IQ)=A(IQ,I) ENDIF ENDDO C C TEST CONVERGENCE IF(D2.LT.D1) GO TO 610 C C END ROTATIONS C 1000 CONTINUE C C RECALCULATE D2 C D2=TZERO DO I=1,NM1 IP1=I+1 DO J=IP1,N D2=D2+SQ(A(I,J)) ENDDO ENDDO C 1010 CONTINUE c if(bflag) xWRITE(6,5020)IROTM,DELTA C C EIGENVALUES AND EIGENVECTORS C 610 DO I=1,N VAL(I)=A(I,I) DO J=1,N A(I,J)=X(I,J) ENDDO ENDDO C RETURN C 5010 FORMAT(//10X,'*** SUBROUTINE ZEIGEN ***'// + ' WARNING - NO CONVERGENCE IN ITERATIONS FOR V'/ + ' ACCURACY PARAMETER DELTA = ',1PE9.2/ + ' NEXT VALUE OF AVERAGED OMEGA MAY BE INCORRECT'//) 5020 FORMAT(//11X,'*** SUBROUTINE ZEIGEN ***' + /' NO CONVERGENCE FOR IROTM =',I3/' DELTA =',1PE11.2/ + ' NEXT VALUE OF AVERAGED OMEGA MAY BE INCORRECT'//) C END C*************************************************************** C CN SUBROUTINE ZPHIN(I,ZR,N1,N2,ZAI,ZPI) SUBROUTINE ZPHI(ZR,ZAI,ZPI) C C COMPUTES AMPLITUDE ZAI AND PHASE ZPI OF COULOMB FUNCTION ZPHI C FOR COMPLEX-RADIAL CO-ORDINATE ZR. C USES DATA IN ARRAY D WHICH IS HELD IN COMMON/CJWBK/ C AND SHOULD HAVE BEEN COMPUTED IN SUBROUTINE INJWBK. C THE STRUCTURE OF ZPHI IS SIMILAR TO THAT OF SUBROUTINE JWBK. CNRB C NEUTRAL CASE ADDED C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16 (Z) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C CN PARAMETER (MX15N=15*MZCHF) PARAMETER (MX15N=15) C COMMON/CJWBK/D(MX15N) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC COMMON/NRBZED/TZED C CN DIMENSION ZAI(30),ZPI(30),ZR(30) DIMENSION ZAI(1),ZPI(1),ZR(1) C C FOR ZPHI, COMMENT-OUT FOR ZPHIN - NRB C N1=1 N2=1 I=1 C C END ZPHI C NP = 0 C JKB=(I-1)*15 E=D(JKB+1) C=D(JKB+4) C IF(TZED.EQ.0)THEN FK=D(JKB+2) IF(C.EQ.0)THEN DO N=N1,N2 NP=NP+1 ZAI(NP)=1.0/SQRT(FK) ZPI(NP)=FK*ZR(NP) ENDDO ELSE DO N=N1,N2 NP=NP+1 ZX=1./ZR(NP) ZXSQ=ZX*ZX ZW=E-C*ZXSQ ZA1=0.0625*(ZX/ZW)**3 ZCC=ZA1*(D(JKB+14)*ZXSQ+D(JKB+12))*ZX ZWH=SQRT(ZW) Z=ZR(NP)*ZWH ZSQ=Z*Z ZP=(ZSQ-0.125-.2083333*C/ZSQ)/Z+D(JKB+15) ZA1=1./(ZR(NP)*FK) ZS=D(JKB+5)*ZA1 ZG=Z*ZA1 ZP=ZP+D(JKB+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) ZPI(NP)=ZP ZAI(NP)=(1.-ZCC)/SQRT(ZWH) ENDDO ENDIF RETURN ENDIF C IF(E.GT.0) THEN IF(C.GT.0) THEN DO N=N1,N2 NP = NP+1 ZX=1./ZR(NP) ZW=E+ZX*(2.-C*ZX) ZWH=SQRT(ZW) Z=ZR(NP)*ZWH FK=D(JKB+2) ZRK=ZR(NP)*FK ZRMC=ZR(NP)-C ZALP=Z+ZRK CK=D(JKB+10) C COMPUTE PHASE ZP=Z+D(JKB+15) C LOG TERM ZB=FK*ZALP IF(ABS(ZB).GT.ACJWBK) THEN ZP=ZP+D(JKB+3)*LOG(1.+ZB) ELSE ZB=-ZB ZP=ZP+ZALP*((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.) END IF C ARCTAN TERM ZA1=1./(ZR(NP)*D(JKB+7)) ZS=D(JKB+5)*(Z-FK*ZRMC)*ZA1 ZG=(CK*Z+ZRMC)*ZA1 ZP=ZP+D(JKB+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) C CAP PHI TERM ZP=ZP+((5.*ZRMC/(Z*Z))-(Z*D(JKB+9)+ZRK*D(JKB+8)+CK)/ X (ZALP*D(JKB+7)))/(24.*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=ZA1*(((D(JKB+14)*ZX+D(JKB+13))*ZX+D(JKB+12))*ZX+ X D(JKB+11)) ZPI(NP)=ZP ZAI(NP)=(1.-ZCC)/SQRT(ZWH) ENDDO ELSE C CASE OF C.EQ.0 AND E.GT.0 C DO N=N1,N2 NP = NP+1 ZX=1./ZR(NP) ZW=2.*ZX+E ZWH=SQRT(ZW) Z=ZR(NP)*ZWH FK=D(JKB+2) ZRK=ZR(NP)*FK ZALP=Z+ZRK C COMPUTE PHASE ZP=Z+D(JKB+15) ZB=FK*ZALP IF(ABS(ZB).GT.ACJWBK) THEN ZP=ZP+D(JKB+3)*LOG(1.+ZB) ELSE ZB=-ZB ZP=ZP+ZALP*((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.) END IF ZP=ZP+1/(4.*ZALP)+(5.*ZR(NP)/(Z*Z)-2.*(Z+ZALP)/ZALP)/(24.*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=ZA1*(-4.*E-3.*ZX) ZPI(NP)=ZP ZAI(NP)=(1.-ZCC)/SQRT(ZWH) ENDDO ENDIF ELSE IF(C.EQ.0) THEN DO N=N1,N2 NP = NP+1 ZX=1./ZR(NP) ZW=2.*ZX ZWH=SQRT(ZW) Z=ZR(NP)*ZWH ZP=2.*Z*(1.+.046875*ZX)+D(JKB+15) ZWMQ=1./SQRT(ZWH) ZET=(1.+.0234375*ZX)*ZWMQ ZAI(NP)=ZET ZPI(NP)=ZP ENDDO ELSE C C CASE OF E.EQ.0 AND C.GT.0 C DO N=N1,N2 NP = NP+1 ZX=1./ZR(NP) ZW=ZX*(2.-C*ZX) ZWH=SQRT(ZW) Z=ZR(NP)*ZWH ZRMC=ZR(NP)-C C COMPUTE PHASE ZP=2.*Z+D(JKB+15) ZA1=1./ZR(NP) ZS=D(JKB+5)*Z*ZA1 ZG=ZRMC*ZA1 ZP=ZP+D(JKB+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) ZP=ZP-(3.*ZR(NP)+C)/(24.*(ZRMC+ZR(NP))*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=((D(JKB+14)*ZX+D(JKB+13))*ZX-3.)*ZX*ZA1 ZAI(NP)=(1.-ZCC)/SQRT(ZWH) ZPI(NP)=ZP ENDDO ENDIF ENDIF RETURN END C*************************************************************** C SUBROUTINE ZSLAG(I,J,NLAG,LIJ,ZS3,ZSD3) C C CALCULATES S INTEGRALS USING LAGUERRE QUADRATURE. ADAPTED FROM STGF. C IMPLICIT REAL*8 (A-H,O-Y) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C IMPLICIT COMPLEX*16 (Z) C PARAMETER (NCHF1=MZCHF+1) C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/KP0,KP1,KP2,KPMM(MZEST),KPMAX,RZERO,RONE,RTWO,H COMMON/CHAN/ECH(NCHF1),EPS(NCHF1),FKNU(NCHF1),CC(NCHF1), 1 RINF(NCHF1),LLCH(NCHF1),ITARG(NCHF1),NCHF,NCHOP,NCHOP1 C ZFK=DCMPLX(FKNU(I),1./FKNU(J)) X=RTWO B=SQRT(2.*X) ZB=(0.,1.)/B ZG=.5*ZFK*B C ZA1=1./ZG ZA2=1.+ZG ZA2=ZA2*ZA2 ZA3=ZB*ZG C NS=NLAG/2 M=NS*(NS-1) N1=M+1 N2=M+NLAG C ZS3=0. ZSD3=0. DO 10 N=N1,N2 U=XLAG(N) ZET=ZA1*(SQRT(ZA2+ZA3*U)-1.) ZMUM=-.5*ZET/(ZB*(1.+ZET*ZG)) ZET=ZET*ZET ZR=ZET*X CALL ZPHI(ZR,ZAI,ZPI) CALL ZTHETA(J,ZR,ZTAJ,ZTDAJ,ZTPJ) ZB1=(ZR**(-LIJ))*WLAG(N)*ZMUM* 1 EXP((0.,1.)*ZPI+ZTPJ+U)*ZAI ZS3=ZS3+ZB1*ZTAJ 10 ZSD3=ZSD3+ZB1*ZTDAJ C RETURN END C*************************************************************** C SUBROUTINE ZTHETA(I,ZR,ZTA,ZTDA,ZTP) C C CALCULATES THETA AND THETAD FOR CHANNEL I AND COMPLEX ZR C THETA = ZTA*CEXP(ZTP) C THETAD = ZTDA**CEXP(ZTP) C ZTP = FNUI*CLOG(ZR) - ZR/FNUI C ADAPTED FROM STGF. C IMPLICIT REAL*8 (A-H,O-Y) C INCLUDE 'PARAM' C PARAMETER (MZMSH0=MZMSH) C IMPLICIT COMPLEX*16 (Z) C PARAMETER (NCHF1=MZCHF+1) C COMMON/CTHET/BB(NCHF1,MZTET),BG(NCHF1,MZTET),MSUM(NCHF1) COMMON/NRBZED/TZED C MI=MSUM(I) FNUI=BB(I,1) E=BG(I,1) F2=BG(I,2) Z=2.*ZR/FNUI ZY=1./Z ZAS=1. ZS=BB(I,2) ZX=0. DO 10 M=3,MI ZAS=ZAS*ZY ZX=BG(I,M)*ZAS+ZX 10 ZS=BB(I,M)*ZAS+ZS C ZDLR=LOG(ZR)*TZED ZTP=FNUI*ZDLR-.5*Z ZTA=ZS ZTDA=((ZR*F2+ZDLR)*ZS+ZX)*E C RETURN END C***************************************************************** C SUBROUTINE ZVERT(V,LV,N,W) C C ________________________________________________________ C | | C | INVERT A GENERAL COMPLEX MATRIX | C | | C | INPUT: | C | | C | V --COMPLEX ARRAY CONTAINING MATRIX | C | | C | LV --LEADING (ROW) DIMENSION OF ARRAY V | C | | C | N --DIMENSION OF MATRIX STORED IN ARRAY V | C | | C | W --INTEGER WORK ARRAY WITH AT LEAST N-1 | C | ELEMENTS | C | | C | OUTPUT: | C | | C | V --INVERSE | C | | C | BUILTIN FUNCTIONS: CABS | C |________________________________________________________| C COMPLEX*16 V(LV,1),Y,Z REAL*8 S,T INTEGER W(1),I,J,K,L,M,N,P IF ( N .EQ. 1 ) GOTO 110 L = 0 M = 1 10 IF ( L .EQ. N ) GOTO 90 K = L L = M M = M + 1 C --------------------------------------- C |*** FIND PIVOT AND START ROW SWAP ***| C --------------------------------------- P = L S = ABS(V(L,L)) IF ( M .GT. N ) GOTO 30 DO 20 I = M,N T = ABS(V(I,L)) IF ( T .LE. S ) GOTO 20 P = I S = T 20 CONTINUE W(L) = P 30 IF ( S .EQ. 0. ) GOTO 120 Y = V(P,L) V(P,L) = V(L,L) C ----------------------------- C |*** COMPUTE MULTIPLIERS ***| C ----------------------------- V(L,L) = -1. Y = 1./Y DO 40 I = 1,N 40 V(I,L) = -Y*V(I,L) J = L 50 J = J + 1 IF ( J .GT. N ) J = 1 IF ( J .EQ. L ) GOTO 10 Z = V(P,J) V(P,J) = V(L,J) V(L,J) = Z IF ( ABS(Z) .EQ. 0. ) GOTO 50 C ------------------------------ C |*** ELIMINATE BY COLUMNS ***| C ------------------------------ IF ( K .EQ. 0 ) GOTO 70 DO 60 I = 1,K 60 V(I,J) = V(I,J) + Z*V(I,L) 70 V(L,J) = Y*Z IF ( M .GT. N ) GOTO 50 DO 80 I = M,N 80 V(I,J) = V(I,J) + Z*V(I,L) GOTO 50 C ----------------------- C |*** PIVOT COLUMNS ***| C ----------------------- 90 L = W(K) DO 100 I = 1,N Y = V(I,L) V(I,L) = V(I,K) 100 V(I,K) = Y K = K - 1 IF ( K .GT. 0 ) GOTO 90 RETURN 110 IF ( ABS(V(1,1)) .EQ. 0. ) GOTO 120 V(1,1) = 1./V(1,1) RETURN 120 WRITE(6,*) 'ERROR: MATRIX HAS NO INVERSE' STOP END