C N. R. BADNELL UoS v2.7 - QUB v1.1 07/10/20 C PROGRAM MAIN IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C REAL*4 TARRY(2) C C C CALCULATES BOUND-BOUND LINE STRENGTHS C C C C COMMON BLOCKS C ************* C C C C CD, DATA FROM D DATASET C == COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 C C C CB, DATA INDEPENDENT OF ENERGY FROM B DATASET C == COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C C C CASE, DATA FOR ENERGY-CASE C ==== LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM C C C CPERT, DATA FOR PERTURBATION POTENTIALS C ===== COMMON/CPERT/ 1 CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2),DMATV(MZCHF,MZCHF,2), 2 WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) C C C CSLPI, SLPI CASES C ===== COMMON/CSLPI/ 1 KOUNT,MSLP(100),MSLP1(100),MSLP2(100), 2 IS1,IL1,IP1,IS2,IL2,IP2,ISLP1,ISLP2,K9999 C C C CSCALE, SCALE FACTORS C ====== COMMON/CSCALE/ 1 AZ,AZP,AZ2 C C C CNTRL, CONTROL C ===== COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C COMMON/NRBZED/TZED C C C DIMENSION C ********* C DIMENSION EST1(MZEST),EST2(MZEST) DIMENSION GFL(MZEST,MZEST),GFV(MZEST,MZEST) C C C NAMELIST C ******** C NAMELIST/STGBB/IPRINT,IBUT,TINY,DELTA,AC C C C I/O UNITS C ********* C C UNIT FILE C C 1 DKK (D01 ETC) C 5 CONTROL OF PROGRAM C 6 PRINTED OUTPUT C 7 FVALUE C 20 B00 C 21 BSLPI1 (B01 ETC) C 22 BSLPI2 (B01 ETC) C 99 D00 C C C OPEN(5,FILE='dstgbb',STATUS='OLD') OPEN(6,FILE='routbb',STATUS='UNKNOWN') C C C WRITE HEADING C ************* C WRITE(6,598) WRITE(6,6002)MZCHF,MZEST,MZMNP,MZPTS,MZSLP,MZTET C C READ IPRINT AND IBUT FROM UNIT 5 C ******************************** C C IBUT=-1 FOR OLD D DATASETS (NO BUTTLE DATA). C IBUT= 0 FOR NEW D DATASETS, NOT TO BE PROCESSED. (DEFAULT) C IBUT=+1 FOR NEW D DATASETS, TO BE PROCESSED. C C IPRINT=-1 JUST WRITE F-VALUES TO FILE FVALUE. (DEFAULT) C IPRINT=1 AND WRITE LINE STRENGTHS TO UNIT6. C IPRINT=0 AND WRITE F-VALUES TO UNIT6 AS WELL. C IBUT=0 IPRINT=-1 TINY=1.E-6 DELTA=0.002 AC=1.E-5 C READ(5,STGBB) C WRITE(6,600)IPRINT,IBUT,TINY,DELTA,AC IF(IBUT.LE.0)WRITE(6,631) C C C READ DATA INDEPENDENT OF SLPI COMBINATION AND OF ENERGY C ******************************************************* C C READ DATA INDEPENDENT OF SLPI FROM B HEADER FILE, UNIT 20 CALL RB00 IF(IPERT.EQ.0)WRITE(6,632) OPEN(7,FILE='FVALUE',STATUS='UNKNOWN') REWIND(7) WRITE(7,700)NZED,NELC C TZED=1.0 IF(NZED.EQ.NELC)TZED=0.0 C C READ DATA FROM D HEADER FILE, UNIT 99 CALL RD00 C C C START PROCESSING SLPI CASES C *************************** C C READ FROM 5 FOR SLPI COMBINATIONS C TERMINATES WITH END-OF-FILE C OR -1 -1 -1 -1 -1 -1 C 9999 READ(5,*,end=9997)IS1,IL1,IP1, IS2,IL2,IP2 C IF(IL1.GE.0)GO TO 9998 9997 WRITE(7,710) 0,0,0, 0,0,0, 0 CLOSE(7,STATUS='KEEP') WRITE(6,640) CLOSE(20,STATUS='KEEP') CLOSE(99,STATUS='KEEP') C C SUN DUM=DTIME(TARRY) TIME=TARRY(1) TIME=TIME/60. WRITE(6,999)TIME 999 FORMAT(//' CPU TIME=',F9.3,' MIN') C STOP C C*********************************** C C C 9998 IF(IPRINT.EQ.1)WRITE(6,608)IS1,IL1,IP1,IS2,IL2,IP2 C C C OPEN FILES CALL FUNIT IF(K9999.EQ.9999)GOTO 9999 C C READ D FILE ON UNIT 1 CALL RDKKP C C READ B DATA INDEPENDENT OF ENERGY C.. FOR FIRST SLPI FROM UNIT 21 CALL RB1IND C.. FOR SECOND SLPI FROM UNIT 22 CALL RB2IND C C ARRAYS FOR IPERT = 1 IF(IPERT.EQ.1)CALL PERTAR C C DEFINE ARRAY LENGTHS ALLOWING FOR BUTTLE CORRECTION C IF(IBUT.LE.0)THEN KM1=MNP2D1 KM2=MNP2D2 ELSE KM1=MNP2D1+NCHND1 KM2=MNP2D2+NCHND2 ENDIF C C C START PROCESSING ENERGY CASES C ***************************** C C WRITE HEADING C WRITE(6,609)IS1,IL1,IP1,IS2,IL2,IP2 IF(IPRINT.EQ.1)WRITE(6,610) C C..START LOOP FOR ENERGY-CASE M1 C DO 100 M1=1,MXE1 C C READ ENERGY-DEPENDENT DATA FROM UNIT 21 CALL RB1DEP(KPM1) EST1(M1)=E1 C C POSITION UNIT 22 REWIND 22 CALL RB2IND C C..START LOOP FOR ENERGY-CASE M2 C DO 100 M2=1,MXE2 C C READ ENERGY-DEPENDENT DATA FROM 22 CALL RB2DEP(KPM2) EST2(M2)=E2 C C INITIALISE FOR ENERGY-CASE KPM=MIN(KPM1,KPM2) DE=E1-E2 DER=1./DE DLI=0. DVI=0. DLA=0. DLB=0. DVA=0. DVB=0. AMV=0. CALL FLAG C C..CALCULATE INNER REGIONS CONTRIBUTION, DVI AND DLI C DO K2=1,KM2 DVI=DVI+AVECT2(K2)*DVECTV(K2,M1) DLI=DLI+AVECT2(K2)*DVECTL(K2,M1) ENDDO DVI=-2.*DVI*DER C C..CALCULATE OUTER-REGION CONTRIBUTIONS C DO N2=1,NCHND2 DO N1=1,NCHND1 C C CALCULATE AMAT CONTRIBUTIONS AML=AMATL(N1,N2) IF(ABS(AML).GT.TINY)THEN IF(ABS(FNU1(N1)-FNU2(N2)).GT.DELTA.AND. ABS(DER).LE.1.E+4)THEN CC = CC1(N1)-CC2(N2) CALL S1S2(CC,P1(N1),PP1(N1),P2(N2),PP2(N2),PPP2(N2),S1,S2) CALL FINT(-2,DA) DAC=(-S2+TZED*4.*DA)*DER*DER DOR=-S1*DER+DAC IF(NSPND.EQ.0)AMV = AMATV(N1,N2) C = 0 IN LS COUPLING (ACCOUNTING FOR 'REVERSE' IN BP-RECUD) IF(AMV.EQ.0.) THEN S3 = DAC*AML ELSE CALL S1S2(-CC,P2(N2),PP2(N2),P1(N1),PP1(N1),PPP1(N1),S1,S3) S3=(S3-DA*4.*TZED)*DER*DER*AMV ENDIF DVA = S3 + DVA ELSE FLAGV=.FALSE. CALL FINT(1,DOR) ENDIF DLA=DLA+AML*DOR ENDIF C C CALCULATE BMAT CONTRIBUTIONS BML=BMATL(N1,N2) BMV=BMATV(N1,N2) IF(ABS(BML).GT.TINY)THEN IF(ABS(FNU1(N1)-FNU2(N2)).GT.DELTA)THEN S3=P1(N1)*PP2(N2)-PP1(N1)*P2(N2) DOV=-S3/(EPS1(N1)-EPS2(N2)) ELSE FLAGP=.FALSE. CALL FINT(0,DOV) ENDIF DLB=DLB+BML*DOV DVB=DVB+BMV*DOV ENDIF C ENDDO ENDDO C IF(.NOT.FLAGV)FLAGP=.FALSE. C C..CALCULATE PERTURBATION CONTRIBUTIONS C C INITIALISATION IF(FLAGP)THEN CALL DMAT DLAP=0. DLBP=0. DVAP=0. DVBP=0. AMV=0. C C START LOOP ON CHANNELS C DO N2=1,NCHND2 DO N1=1,NCHND1 C C AMAT CONTRIBUTIONS AML=AMATL(N1,N2) IF(ABS(AML).GT.TINY)THEN CC = CC1(N1)-CC2(N2) IF(NSPND.EQ.0)AMV=AMATV(N1,N2) IF(AMV.EQ.0.) THEN CALL S1S2(CC,P1(N1),PP1(N1),DP2(N2),DPP2(N2),DPPP2(N2),S12,S22) CALL S1S2(CC,DP1(N1),DPP1(N1),P2(N2),PP2(N2),PPP2(N2),S11,S21) DAC=-(S21+S22)*DER*DER DOR=-(S11+S12)*DER+DAC DVAP=DVAP-AML*DAC !+/- ELSE CALL S1S2(-CC,P2(N2),PP2(N2),DP1(N1),DPP1(N1),DPPP1(N1),S12,S22) CALL S1S2(-CC,DP2(N2),DPP2(N2),P1(N1),PP1(N1),PPP1(N1),S11,S21) DAC=-(S21+S22)*DER*DER DOR= (S11+S12)*DER+DAC DVAP=DVAP+AMV*DAC !+/- ENDIF DLAP=DLAP-AML*DOR !+/- ENDIF C C BMAT CONTRIBUTIONS BML=BMATL(N1,N2) BMV=BMATV(N1,N2) IF(ABS(BML).GT.TINY)THEN S3=P1(N1)*DPP2(N2)-PP1(N1)*DP2(N2) 1 +DP1(N1)*PP2(N2)-DPP1(N1)*P2(N2) DOV=-S3/(EPS1(N1)-EPS2(N2)) DLBP=DLBP+BML*DOV DVBP=DVBP+BMV*DOV ENDIF C C CMAT CONTRIBUTIONS CM=CMAT(N1,N2) IF(ABS(CM).GT.TINY)THEN CALL FINT(-3,SP) SP=SP*DER*DER*CM C IF(ABS(SP).GE.1.E-3)PRINT*,' N1,N2=',N1,N2, C + ' CM CONT. TO DLAP=',SP DLAP=DLAP+SP DVAP=DVAP+SP ENDIF C C DMAT CONTRIBUTIONS DML1=DMATL(N1,N2,1) DML2=DMATL(N1,N2,2) DMV1=DMATV(N1,N2,1) DMV2=DMATV(N1,N2,2) IF(ABS(DML1).GT.TINY)THEN CALL FINT(-2,SP) DLBP=DLBP+SP*DML1 DVBP=DVBP+SP*DMV1 ENDIF IF(ABS(DML2).GT.TINY)THEN CALL FINT(-3,SP) DLBP=DLBP+SP*DML2 DVBP=DVBP+SP*DMV2 ENDIF C ENDDO ENDDO C C ADD PERTURBATION CONTRIBUTIONS C DLP=DLAP+DLBP DVP=DVAP-2.*DVBP*DER C ELSE DLP=0. DVP=0. C ENDIF C C..ADD CONTRIBUTIONS AND CALCULATE LINE STRENGTH C DL=DLI+DLA+DLB+DLP SL=DL*DL if(nspnd.ne.0) sl=sl*nspnd GFL(M1,M2)=.333333333*DE*SL SL=SL*AZ2 IF(IPRINT.EQ.1)WRITE(6,620)M1,M2,DLI,DLA,DLB,DLP,DL,SL IF(FLAGV)THEN DVB=-2.*DVB*DER DV=DVI+DVA+DVB+DVP SV=DV*DV if(nspnd.ne.0) sv=sv*nspnd GFV(M1,M2)=.333333333*DE*SV SV=SV*AZ2 IF(IPRINT.EQ.1)WRITE(6,621)DVI,DVA,DVB,DVP,DV,SV ELSE GFV(M1,M2)=0. IF(IPRINT.EQ.1)WRITE(6,*) ENDIF 100 CONTINUE C C WRITE F-VALUES C IF(IPRINT.EQ.0)THEN WRITE(6,629) WRITE(6,630)((M1,M2,EST1(M1),EST2(M2),GFL(M1,M2),GFV(M1,M2), 1 M2=1,MXE2),M1=1,MXE1) ENDIF WRITE(7,710)IS1,IL1,IP1,IS2,IL2,IP2,mxe1*mxe2 WRITE(7,730)((M1,M2,GFL(M1,M2),GFV(M1,M2), 1 M2=1,MXE2),M1=1,MXE1) C C CLOSE FILES C CLOSE (1,STATUS='KEEP') CLOSE(21,STATUS='KEEP') CLOSE(22,STATUS='KEEP') C C GO TO NEXT SLPI CASE C GOTO 9999 C C FORMATS C ******* C 598 FORMAT(////1X,70(1H+)/1X,70(1H+)///5X, 1 'PROGRAM STGBB, BOUND-BOUND TRANSITION DATA'/5X,42(1H+)///) 600 FORMAT(5X,'DATA READ FROM UNIT 5:'// 1 10X,'IPRINT = ',I2/ 2 10X,'IBUT = ',I2/ 3 10X,'TINY = ',1PE9.2/ 4 10X,'DELTA = ',E9.2/ 5 10X,'AC = ',E9.2) 6002 FORMAT(//10X,'COMPILED FOR DIMENSIONS -'// + 15X,'CHANNELS MZCHF = ',I6/ + 15X,'BOUND ENERGIES/SYMMETRY MZEST = ',I6/ + 15X,'R-MATRIX POLES MZMNP = ',I6/ + 15X,'OUTER-REGION RADIAL POINTS MZPTS = ',I6/ + 15X,'S, L, PI CASES MZSLP = ',I6/ + 15X,'COEFFICIENTS FOR THETA MZTET = ',I6//) 608 FORMAT(//70(1H*)//10X,'READ FROM UNIT 5 -'/ 1 10X,'(IS1,IP1,IL1) = (',3I3,' )'/ 2 10X,'(IS2,IL2,IP2) = (',3I3,' )') 609 FORMAT(//70(1H*)//20X,'BOUND-BOUND TRANSITION DATA FOR'// 1 10X,'(IS1,IL1,IP1) = (',3I3,' )',5X, 2 '(IS2,IL2,IP2) = (',3I3,' )'//) 610 FORMAT(//' I J TYPE',3X,'DI',9X,'DA',9X, 1 'DB',9X,'DP',8X,'D',10X,'S'/) 613 FORMAT(//' READS MNP2 = ',I4,'WHICH IS LARGER THAN' 1,'MAXIMUM VALUE OF MZMNP ALLOWED BY DIMENSIONS'//) 620 FORMAT(2I3,2X,'L',1P3E11.3,E10.2,E11.3,E11.3) 621 FORMAT(8X,'V',1P3E11.3,E10.2,E11.3,E11.3/) 629 FORMAT(3X,'M1',3X,'M2',8X,'E1',13X,'E2', 1 13X,'GFL',10X,'GFV'/) 630 FORMAT(2I5,1P2E15.5,2E13.3) 631 FORMAT(//5X,'*******************************************'/ + 5X,'* CALCULATIONS WITHOUT BUTTLE CORRECTIONS *',/ + 5X,'*******************************************'//) 632 FORMAT(//5X,'*******************************'/ + 5X,'* CALCULATIONS WITH IPERT = 0 *'/ + 5X,'*******************************'//) 640 FORMAT(//10X,23(1H*)/10X,'* FILE FVALUE WRITTEN *'/ + 10X,23(1H*)//) 700 FORMAT(2I5,' F') 710 FORMAT(3I3,2x,3i3,i10) 730 FORMAT(2I5,1P2E12.3) C END C C********************************************************************* C SUBROUTINE DMAT IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 C LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C COMMON/CPERT/ 1 CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2),DMATV(MZCHF,MZCHF,2), 2 WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) C C CALCULATE ARRAYS DMATL, DMATV 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)-EPS2(I2)) DV=BMATV(J1,I2)*D DO I1=1,NCHND1 X=WMAT1(I1,J1,L)*BMATL(J1,I2) IF(ABS(X).GT.TINY.AND.ABS(FNU1(J1)-FNU2(I2)).LT.DELTA)THEN FLAGP=.FALSE. RETURN ENDIF DMATL(I1,I2,L)=DMATL(I1,I2,L)+X*D 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)-EPS2(J2)) X=BMATL(I1,J2)*WMAT2(J2,I2,L) IF(ABS(X).GT.TINY.AND.ABS(FNU1(I1)-FNU2(J2)).LT.DELTA)THEN FLAGP=.FALSE. RETURN ENDIF DMATL(I1,I2,L)=DMATL(I1,I2,L)-X*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********************************************************************* C SUBROUTINE FINT(MP,FI) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C C CALCULATES INTERGAL FROM RZERO TO INFINITY OF C F1*(R**MP)*F2 C LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C C CALCULATE INTERGAL FROM RZERO TO RTWO FI=0. DO 10 K=1,KPM 10 FI=FI+WBODE(K,MP)*FSD1(N1,K)*FSD2(N2,K) IF(.NOT.FLAGM)RETURN C C IF NECESSARY, ADD CONTRIBUTION FROM RTWO TO INFINITY IF(FLAG1(N1).AND.FLAG2(N2))THEN FI0=FI CALL TLAG3(MP,FI) CW IF(IPRINT.EQ.1)WRITE(6,6000)N1,N2,MP,FNU1(N1),FNU2(N2),FI0,FI ENDIF C RETURN C 6000 FORMAT(' TLAG3 CALLED: N1=',I3,' ,N2=',I3,' ,MP=',I3,' ,FNU1=', 1 F7.4, 2 ' ,FNU2=',F7.4/' BEFORE CALL, FI=',E12.4,' :AFTER CALL, FI=', 3 E12.4/) C END C C********************************************************************* C SUBROUTINE FLAG IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C C SETS UP ARRAYS EPS1, EPS2, FNU1, FNU2 C SETS FLAG1 AND FLAG2 C COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C COMMON/CPERT/ 1 CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2),DMATV(MZCHF,MZCHF,2), 2 WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) C COMMON/NRBZED/TZED C C SET FLAGM IF(KPM.LT.KP2)THEN FLAGM=.FALSE. DO 1 N1=1,NCHND1 EPS1(N1)=E1-ECH1(N1) 1 FNU1(N1)=SQRT(-1./EPS1(N1)) DO 2 N2=1,NCHND2 EPS2(N2)=E2-ECH2(N2) 2 FNU2(N2)=SQRT(-1./EPS2(N2)) GOTO 30 ELSE FLAGM=.TRUE. ENDIF C C SET FLAG1 DO 10 N1=1,NCHND1 EPS=E1-ECH1(N1) EPS1(N1)=EPS FNU1(N1)=SQRT(-1./EPS) EC=EPS*CC1(N1) FLAG1(N1)=.FALSE. IF(TZED.EQ.0.)THEN RI=RZERO ELSE IF(EC.LE.-1.)GO TO 10 RI=-(1.+SQRT(1.+EC))/EPS ENDIF IF(RI.GT.RTWO)THEN FLAG1(N1)=.TRUE. ELSE IF(RI.GT.RZERO)THEN KI=1+(RI-RZERO)/H ELSE KI=1 ENDIF IF(FSD1(N1,KP2).LE.AC*FSD1(N1,KI))THEN FLAG1(N1)=.FALSE. ELSE FLAG1(N1)=.TRUE. ENDIF ENDIF 10 CONTINUE C C SET FLAG2 DO 20 N2=1,NCHND2 EPS=E2-ECH2(N2) EPS2(N2)=EPS FNU2(N2)=SQRT(-1./EPS) EC=EPS*CC2(N2) FLAG2(N2)=.FALSE. IF(TZED.EQ.0.)THEN RI=RZERO ELSE IF(EC.LE.-1.)GO TO 20 RI=-(1.+SQRT(2.+EC))/EPS ENDIF IF(RI.GT.RTWO)THEN FLAG2(N2)=.TRUE. ELSE IF(RI.GT.RZERO)THEN KI=1+(RI-RZERO)/H ELSE KI=1 ENDIF IF(FSD2(N2,KP2).LE.AC*FSD2(N2,KI))THEN FLAG2(N2)=.FALSE. ELSE FLAG2(N2)=.TRUE. ENDIF ENDIF 20 CONTINUE C C C INITIALISE FLAGV AND FLAGP 30 FLAGV=.TRUE. IF(IPERT.EQ.1)THEN FLAGP=.TRUE. ELSE FLAGP=.FALSE. ENDIF C RETURN END C C********************************************************************* C SUBROUTINE FUNIT IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C C FINDS UNITS DKK, BSLPI1, BSLPI2 C COMMON/CSLPI/ 1 KOUNT,MSLP(100),MSLP1(100),MSLP2(100), 2 IS1,IL1,IP1,IS2,IL2,IP2,ISLP1,ISLP2,K9999 C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C CHARACTER BSLPI1*3,BSLPI2*3,DKK*3,NUM(0:9) DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C C INITIALISATION K9999=0 ISLP1=10000*IS1+100*IL1+IP1 ISLP2=10000*IS2+100*IL2+IP2 C C FIND UNIT DKK FOR SLPI COMBINATION DO 7 K=1,KOUNT IF(ISLP1.EQ.MSLP1(K).AND.ISLP2.EQ.MSLP2(K))THEN K1=K GOTO 9 ENDIF 7 CONTINUE DO 8 K=1,KOUNT IF(ISLP1.EQ.MSLP2(K).AND.ISLP2.EQ.MSLP1(K))THEN I=IS1 IS1=IS2 IS2=I I=IL1 IL1=IL2 IL2=I I=IP1 IP1=IP2 IP2=I I=ISLP1 ISLP1=ISLP2 ISLP2=I K1=K GOTO 9 ENDIF 8 CONTINUE C REQUIRED COMBINATION NOT FOUND ON UNIT DKK WRITE(6,604)IS1,IL1,IP1,IS2,IL2,IP2 K9999=9999 RETURN C C FIND UNIT MFB1 FOR FIRST SLPI 9 CONTINUE DO 50 K=1,KSLP IF(ISLP1.EQ.MSLP(K))THEN K21=K GOTO 51 ENDIF 50 CONTINUE WRITE(6,605)IS1,IL1,IP1 K9999=9999 RETURN C C FIND UNIT MFB2 FOR SECOND SLPI 51 CONTINUE DO 52 K=1,KSLP IF(ISLP2.EQ.MSLP(K))THEN K22=K GOTO 53 ENDIF 52 CONTINUE WRITE(6,605)IS2,IL2,IP2 K9999=9999 RETURN C C FILES NAMES 53 DKK='D'//NUM(K1/10)//NUM(K1-10*(K1/10)) BSLPI1='B'//NUM(K21/10)//NUM(K21-10*(K21/10)) BSLPI2='B'//NUM(K22/10)//NUM(K22-10*(K22/10)) IF(IPRINT.EQ.1)WRITE(6,607)BSLPI1,IS1,IL1,IP1,BSLPI2,IS2,IL2, + IP2,DKK C C OPEN FILES C..UNIT 1 FOR DKKP DATA OPEN(1,FILE=DKK,STATUS='OLD',FORM='UNFORMATTED',ERR=100) REWIND(1) C..UNIT 21 FOR FIRST BSLPI FILE 54 OPEN(21,FILE=BSLPI1,STATUS='OLD',FORM='UNFORMATTED',ERR=200) REWIND(21) C..UNIT 22 FOR SECOND BSLPI FILE 55 OPEN(22,FILE=BSLPI2,STATUS='OLD',FORM='UNFORMATTED',ERR=300) REWIND(22) C RETURN C 100 WRITE(6,6100)DKK,IS1,IL1,IP1,IS2,IL2,IP2 K9999=9999 GOTO 54 200 WRITE(6,6200)BSLPI1,IS1,IL1,IP1 K9999=9999 GOTO 55 300 WRITE(6,6200)BSLPI2,IS2,IL2,IP2 IF(K9999.EQ.9999)THEN CLOSE(1,STATUS='KEEP') CLOSE(21,STATUS='KEEP') ENDIF K9999=9999 RETURN C C C FORMATS 604 FORMAT(//3X,5(1H*),2X,'DATA ON DKK FILES NOT FOUND FOR'/ 1 10X,'IS1, IL1, IP1 = ',3I3,5X,'IS2, IL2, IP2 = ',3I3//) 605 FORMAT(//3X,5(1H*),' DATA FOR IS, IL,IP =',3I3, 1 ' NOT ON B DATASET') 607 FORMAT(//10X,'REQUIRED FILES ARE -'/ 1 10X,A3,' FOR IS, IL, IP = ',3I3/ 2 10X,A3,' FOR IS, IL, IP = ',3I3/ 3 10X,A3,' FOR COMBINATION'//) 6100 FORMAT(/3X,5(1H*),2X,'FILE ',A3,' NOT FOUND'/ 1 10X,'REQUIRED FOR COMBINATION IS1, IL1, IP1 = ',3I3/ 2 35X,'IS2, IL2, IP2 = ',3I3/) 6200 FORMAT(3X,5(1H*),2X,'FILE ',A3,' NOT FOUND, REQUIRED FOR' 1 ,' IS, IL, IP = ',3I3//) C END C C********************************************************************* C SUBROUTINE PERTAR IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 C COMMON/CPERT/ 1 CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2),DMATV(MZCHF,MZCHF,2), 2 WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C COMMON/NRBZED/TZED C C ARRAY CMAT 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********************************************************************* C SUBROUTINE RB00 IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C C READS DATA INDEPENDENT OF SLPI FROM FIRST B FILE, UNIT 20 C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C COMMON/CSCALE/ 1 AZ,AZP,AZ2 C COMMON/CSLPI/ 1 KOUNT,MSLP(100),MSLP1(100),MSLP2(100), 2 IS1,IL1,IP1,IS2,IL2,IP2,ISLP1,ISLP2,K9999 C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C C FIRST, DATA INDEPENDENT OF SLPI C OPEN(20,FILE='B00',STATUS='OLD',FORM='UNFORMATTED',ERR=100) C REWIND(20) READ(20)NZED,NELC READ(20)KP2 IF(KP2.GT.MZPTS)THEN WRITE(6,*)'INCREASE MZPTS TO ',KP2 STOP 'INCREASE MZPTS' ENDIF READ(20)RZERO,H READ(20)(WBODE(IR,0),IR=1,KP2) READ(20)IPERT WRITE(6,612)NZED,NELC IF(IPRINT.EQ.1)WRITE(6,630)RZERO,H,KP2 IF(IPRINT.EQ.1)WRITE(6,602) KSLP=0 C C THEN LIST OF SLPI CASES C 5 READ(20,end=509)IS,IL,IP IF(IL.GE.0)THEN KSLP=KSLP+1 IF=KSLP MSLP(KSLP)=10000*IS+100*IL+IP IF(IF.GE.10)THEN IF(IPRINT.EQ.1)WRITE(6,603)KSLP,IS,IL,IP,IF ELSE IF(IPRINT.EQ.1)WRITE(6,604)KSLP,IS,IL,IP,IF ENDIF GOTO 5 ENDIF 509 CLOSE(20) C C DATA FOR SCALING C AZ=MAX(NZED-NELC,1) AZP=1./AZ AZ2=AZP*AZP 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 RETURN C 100 WRITE(6,600) STOP C C FORMATS 600 FORMAT(//10X,40(1H*)//10X,'B HEADER FILE B00 NOT FOUND'// 1 10X,40(1H*)//) 602 FORMAT(//10X,'FROM B DATASET'/10X,'KSLP',3X,'IS',3X,'IL', 1 3X,'IP',10X,'FILE'/) 603 FORMAT(10X,I3,4X,I2,3X,I2,3X,I2,10X,'B',I2) 604 FORMAT(10X,I3,4X,I2,3X,I2,3X,I2,10X,'B0',I1) 612 FORMAT(//5X,'NZED = ',I3,' NELC = ',I3/5X,22(1H*)//) 630 FORMAT(/12X,' RZERO = ',F8.4/11X,' H =',F8.4/11X,' KP2 = ', 1 I4/) C END C C********************************************************************* C SUBROUTINE RB1DEP(KPM1) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) 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 C READ ENERGY-DEPENDENT DATA FROM UNIT 21 C LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/NRBZED/TZED DIMENSION XVECT(MZCHF) C V(X)=EQ+X*(Q2*TZED-X*CQ) C TOL0=-1.D-150 !INITIALIZE DEEPLY-CLOSED AS OFF (ZERO) C READ(21,ERR=100,END=100) E1,(XVECT(N),N=1,NCHND1) TOL0=ABS(TOL0) !SWITCH-ON AS CAN RENORM (NEW STGB/BMN) 100 READ(21) !WAS AVECT1 READ(21) (P1(N),N=1,NCHND1) READ(21) (PP1(N),N=1,NCHND1) READ(21) (PPP1(N),N=1,NCHND1) READ(21)KPM1 C IF(KPM1.GT.1)THEN READ(21)(FSD1(I,KPM1-1),FSD1(I,KPM1),I=1,NCHND1) Q=H*H*.083333333333 Q2=Q*2. C DO I=1,NCHND1 EQ=Q*(E1-ECH1(I)) CQ=Q*CC1(I) KP2T=KPM1 RTWOT=RZERO+H*(KPM1-1) IFLAG=0 C IF(ABS(FSD1(I,KP2T))+ABS(FSD1(I,KP2T-1)).LT.TOL0)THEN IFLAG=1 C C MOVE-IN AND REGENERATE UNTIL NON-VANISHING. NRB 16/06/97 C LL=(SQRT(ONE+FOUR*CC1(I))-ONE+PT02)/TWO FKNU=SQRT(-ONE/(E1-ECH1(I))) 1 CALL THETBB(RTWOT,FKNU,LL,T,TP,TD,TDP) IF(ABS(T)+ABS(TP).LT.TOL0)THEN KP2T=KP2T-4 IF(KP2T.LT.1)KP2T=1 RTWOT=(KP2T-1)*H+RZERO DO K=KP2T,KP2T+4 FSD1(I,K)=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(I,KP2T)=C1 FSD1(I,KP2T-1)=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(I,K)=((TWO-TEN*U2)*FSD1(I,K+1)-(ONE+U1)*FSD1(I,K+2)) + /(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(I,K)=(T2*FSD1(I,K+1)-T1*FSD1(I,K+2))/T3 C5=(T2*C4-T1*C3-Q*(FSD1(I,K)+TEN*FSD1(I,K+1) X +FSD1(I,K+2)))/T3 C3=C4 C4=C5 U1=U2 U2=U3 ENDDO T=FSD1(I,1) TP=(FSD1(I,2)-FSD1(I,1))/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(I,K)=FSD1(I,K)*W1 ENDDO ENDIF 210 CONTINUE ENDDO if(iprint.gt.1)then !check regeneration with stgb write(6,701)e1 write(6,705)rtwo,kp2,h do i=1,nchnd1 write(6,710)i,fsd1(i,1)/xvect(i),fsd1(i,kpm2)/xvect(i) enddo endif 701 format(//10x,' e = ',f10.6/11x,14('-')) 705 format(//' rtwo = ',f10.6,', kp2 = ',i4,', h = ',f10.6) 710 format(i5,4e15.6) ELSE READ(21)(FSD1(I,1),I=1,NCHND1) ENDIF C IF(IPERT.EQ.1) THEN READ(21)(DP1(N),N=1,NCHND1) READ(21)(DPP1(N),N=1,NCHND1) READ(21)(DPPP1(N),N=1,NCHND1) ENDIF C RETURN END C C********************************************************************* C SUBROUTINE RB1IND IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C C READS DATA INDEPENDENT OF ENERGY FROM FILE BSLPI1 ON UNIT 21 C COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 C COMMON/CPERT/ 1 CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2),DMATV(MZCHF,MZCHF,2), 2 WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C READ(21)IS,IL,IP READ(21)MNP2,NCHF C C CHECK ACCORD BETWEEN DKK AND BSLPI FILES IF(IS.NE.NSPND.OR. 1 IL.NE.LRGLD1.OR. 2 IP.NE.NPTYD1.OR. 3 MNP2.NE.MNP2D1)THEN WRITE(6,614)21,IS,IL,IP,MNP2 STOP ENDIF C READ(21)(ECH1(I),I=1,NCHF) READ(21)(CC1(I),I=1,NCHF) IF(IPERT.EQ.1)READ(21)(((WMAT1(I,J,L),I=1,NCHF),J=1,NCHF),L=1,2) READ(21)MXE1 C 614 FORMAT(//' *** MISMATCH ***'//15X,'UNIT MFA = ',I3, 1 ' HAS IS, IL, IP, MNP2 = ',4I5//) C RETURN END C C********************************************************************* C SUBROUTINE RB2DEP(KPM2) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) 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 C READ ENERGY-DEPENDENT DATA FROM UNIT 22 C LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 COMMON/NRBZED/TZED DIMENSION XVECT(MZCHF) C V(X)=EQ+X*(Q2*TZED-X*CQ) C TOL0=-1.D-150 !INITIALIZE DEEPLY-CLOSED AS OFF (ZERO) C READ(22,ERR=100,END=100) E2,(XVECT(N),N=1,NCHND2) TOL0=ABS(TOL0) !SWITCH-ON AS CAN RENORM (NEW STGB/BMN) 100 READ(22) (AVECT2(K2),K2=1,MNP2D2+NCHND2) READ(22) (P2(N),N=1,NCHND2) READ(22) (PP2(N),N=1,NCHND2) READ(22) (PPP2(N),N=1,NCHND2) READ(22)KPM2 C IF(KPM2.GT.1)THEN READ(22)(FSD2(I,KPM2-1),FSD2(I,KPM2),I=1,NCHND2) Q=H*H*.083333333333 Q2=Q*2. C DO I=1,NCHND2 EQ=Q*(E2-ECH2(I)) CQ=Q*CC2(I) KP2T=KPM2 RTWOT=RZERO+H*(KPM2-1) IFLAG=0 C IF(ABS(FSD2(I,KP2T))+ABS(FSD2(I,KP2T-1)).LT.TOL0)THEN IFLAG=1 C C MOVE-IN AND REGENERATE UNTIL NON-VANISHING. NRB 16/06/97 C LL=(SQRT(ONE+FOUR*CC2(I))-ONE+PT02)/TWO FKNU=SQRT(-ONE/(E2-ECH2(I))) 1 CALL THETBB(RTWOT,FKNU,LL,T,TP,TD,TDP) IF(ABS(T)+ABS(TP).LT.TOL0)THEN KP2T=KP2T-4 IF(KP2T.LT.1)KP2T=1 RTWOT=(KP2T-1)*H+RZERO DO K=KP2T,KP2T+4 FSD2(I,K)=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 FSD2(I,KP2T)=C1 FSD2(I,KP2T-1)=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) FSD2(I,K)=((TWO-TEN*U2)*FSD2(I,K+1)-(ONE+U1)*FSD2(I,K+2)) + /(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) FSD2(I,K)=(T2*FSD2(I,K+1)-T1*FSD2(I,K+2))/T3 C5=(T2*C4-T1*C3-Q*(FSD2(I,K)+TEN*FSD2(I,K+1) X +FSD2(I,K+2)))/T3 C3=C4 C4=C5 U1=U2 U2=U3 ENDDO T=FSD2(I,1) TP=(FSD2(I,2)-FSD2(I,1))/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 FSD2(I,K)=FSD2(I,K)*W1 ENDDO ENDIF 210 CONTINUE ENDDO if(iprint.gt.1)then !check regeneration with stgb write(6,701)e2 write(6,705)rtwo,kp2,h do i=1,nchnd2 write(6,710)i,fsd2(i,1)/xvect(i),fsd2(i,kpm2)/xvect(i) enddo endif 701 format(//10x,' e = ',f10.6/11x,14('-')) 705 format(//' rtwo = ',f10.6,', kp2 = ',i4,', h = ',f10.6) 710 format(i5,4e15.6) ELSE READ(22)(FSD2(I,1),I=1,NCHND2) ENDIF C IF(IPERT.EQ.1) THEN READ(22)(DP2(N),N=1,NCHND2) READ(22)(DPP2(N),N=1,NCHND2) READ(22)(DPPP2(N),N=1,NCHND2) ENDIF C RETURN END C C********************************************************************* C SUBROUTINE RB2IND IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C C READS DATA INDEPENDENT OF ENERGY FROM FILE BSLPI2 ON UNIT 22 C COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C COMMON/CPERT/ 1 CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2),DMATV(MZCHF,MZCHF,2), 2 WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) C READ(22)IS,IL,IP READ(22)MNP2,NCHF C CHECK ACCORD BETWEEN DKK AND BSLPI FILES IF(IS.NE.NSPND.OR. 1 IL.NE.LRGLD2.OR. 2 IP.NE.NPTYD2.OR. 3 MNP2.NE.MNP2D2)THEN WRITE(6,614)22,IS,IL,IP,MNP2 STOP ENDIF C READ(22)(ECH2(I),I=1,NCHF) READ(22)(CC2(I),I=1,NCHF) IF(IPERT.EQ.1)READ(22)(((WMAT2(I,J,L),I=1,NCHF),J=1,NCHF),L=1,2) READ(22)MXE2 C 614 FORMAT(//' *** MISMATCH ***'//15X,'UNIT MFA = ',I3, 1 ' HAS IS, IL, IP, MNP2 = ',4I5//) C RETURN END C C********************************************************************* C SUBROUTINE RD00 IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C C READ DATA FROM UNIT NFD C COMMON/CSLPI/ 1 KOUNT,MSLP(100),MSLP1(100),MSLP2(100), 2 IS1,IL1,IP1,IS2,IL2,IP2,ISLP1,ISLP2,K9999 C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C C OPEN DKKP HEADER FILE D00 ON UNIT 99 OPEN(99,FILE='D00',STATUS='OLD',FORM='UNFORMATTED',ERR=100) REWIND(99) C C READ(99)KOUNT IF(IPRINT.EQ.1)WRITE(6,600) DO 6 K=1,KOUNT READ(99)IS1,IL1,IP1,IS2,IL2,IP2 IF(K.LT.10)THEN IF(IPRINT.EQ.1)WRITE(6,601)K,IS1,IL1,IP1,IS2,IL2,IP2,K ELSE IF(IPRINT.EQ.1)WRITE(6,602)K,IS1,IL1,IP1,IS2,IL2,IP2,K ENDIF MSLP1(K)=10000*IS1+100*IL1+IP1 MSLP2(K)=10000*IS2+100*IL2+IP2 6 CONTINUE C RETURN C 100 WRITE(6,6000) STOP C C FORMATS 600 FORMAT(//10X,'FROM DKK DATASET'// 1 10X,'K',8X,'IS1',2X,'IL1',2X,'IP1',6X,'IS2',2X, 2 'IL2',2X,'IP2',10X,'FILE'/) 601 FORMAT(8X,I3,8X,I3,2X,I3,2X,I3,6X,I3,2X,I3,2X,I3,10X,'D0',I1) 602 FORMAT(8X,I3,8X,I3,2X,I3,2X,I3,6X,I3,2X,I3,2X,I3,10X,'D',I2) 6000 FORMAT(10X,40(1H*)//10X,'DKKP HEADER FILE D00 NOT FOUND'// 1 10X,40(1H*)//) C END C C*************************************************************** C SUBROUTINE RDKKP C C ADAPTED FROM C PROGRAM OF KTT FOR READING DIPOLE MATRIX ELEMENT DATA C AND KB/WE FOR REDUCED MEMORY - ERRORS CORRECTED BY NRB 10/03/99. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) COMMON/CD/ 3 DVECTL(MNPEXT,MZEST),DVECTV(MNPEXT,MZEST), 4 AMATL(MZCHF,MZCHF),AMATV(MZCHF,MZCHF), 4 BMATL(MZCHF,MZCHF),BMATV(MZCHF,MZCHF), 1 MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 2 MNP2D2,NCHND2,LRGLD2,NPTYD2 C COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC COMMON/CSCALE/ 1 AZ,AZP,AZ2 COMMON/CPERT/ 1 CMAT(MZCHF,MZCHF),DMATL(MZCHF,MZCHF,2),DMATV(MZCHF,MZCHF,2), 2 WMAT1(MZCHF,MZCHF,2),WMAT2(MZCHF,MZCHF,2) C DIMENSION DKL(MZCHF,MNPEXT),DKV(MZCHF,MNPEXT) !USE MZCHF NOT MZNR2 - NRB X ,CGC(MZCHF) C C READ(1) NOTERM,MNP2D2,NCHND2,LRGLD2,NPTYD2,NSPND, & MNP2D1,NCHND1,LRGLD1 C C CHECK DIMENSIONS FOR MNP2,MZCHF N=MAX(NCHND1,NCHND2) M=MAX(MNP2D1,MNP2D2) IF(M.GT.MZMNP.OR.N.GT.MZCHF) THEN WRITE(6,600)MNP2D1,NCHND1,MNP2D2,NCHND2 STOP ENDIF C C CHECK DIMENSION FOR "NR2". C IN SMALL CASES THE NO OF CHANNELS WILL BE LESS THAN THE NO OF BASIS C ORBITALS BUT IN LARGE CASES THE NO OF CHANNELS WILL BE LARGER. C THE BUTTLE CORRECTION NEEDS MZCHF. WHICHEVER IS USED (MZCHF OR MZNR2) C DIMENSIONS NEED TO BE TESTED-FOR. NRB IF(NOTERM.GT.MZCHF)THEN WRITE(6,601)NOTERM STOP ENDIF C C CALCULATE NPTYD1 C NPTYD1=1-NPTYD2+2*(NPTYD2/2) C C GET THE ENERGY DEPENDENT BOUND STATE DATA C READ(21)IS,IL,IP READ(21)MNP2,NCHF C C CHECK ACCORD BETWEEN DKK AND BSLPI FILES C IF(IS.NE.NSPND.OR. & IL.NE.LRGLD1.OR. & IP.NE.NPTYD1.OR. MNP2.NE.MNP2D1)THEN WRITE(6,614)21,IS,IL,IP,MNP2 STOP ENDIF READ(21)(ECH1(I),I=1,NCHF) READ(21)(CC1(I),I=1,NCHF) IF(IPERT.EQ.1)READ(21)(((WMAT1(I,J,L),I=1,NCHF),J=1,NCHF),L=1,2) READ(21)MXE1 C DO M1 = 1, MXE1 READ(21) E1 READ(21) (AVECT1(K1,M1),K1=1,MNP2D1+NCHND1) READ(21) READ(21) READ(21) READ(21) READ(21) IF(IPERT.EQ.1) THEN READ(21) READ(21) READ(21) ENDIF ENDDO REWIND 21 DO M1 = 1, MXE1 DO K2=1,MZMNP DVECTV(K2,M1)=0.D0 DVECTL(K2,M1)=0.D0 ENDDO ENDDO C C------THIS TEST PROGRAM TAKES IN THE PROCESSED DIPOLE MATRIX C------ELEMENTS C C IAIN1 = (MNP2D1 - 1) / NOTERM IBIN1 = ( MNP2D2 - 1) / NOTERM MCI = 0 DO IK = 1 , IAIN1 MCH = MCI + 1 MCI = MCI + NOTERM NCI = 0 DO JK = 1 , IBIN1 NCH = NCI + 1 NCI = NCI + NOTERM READ(1) ((DKL(I,J),J=1,NOTERM),I=1,NOTERM), & ((DKV(I,J),J=1,NOTERM),I=1,NOTERM) DO M1 = 1, MXE1 K1P = 1 DO K1=MCH,MCI AVL=AVECT1(K1,M1)*AZ AVV=AVECT1(K1,M1)*AZP K2P=1 DO K2=NCH,NCI DVECTV(K2,M1)=DVECTV(K2,M1)+AVV*DKV(K1P,K2P) DVECTL(K2,M1)=DVECTL(K2,M1)+AVL*DKL(K1P,K2P) K2P = K2P + 1 ENDDO K1P = K1P+1 ENDDO ENDDO ENDDO NCH = NCI + 1 NCI = MNP2D2 NCP = NCI - NCH + 1 READ(1) ((DKL(I,J),J=1,NCP),I=1,NOTERM), & ((DKV(I,J),J=1,NCP),I=1,NOTERM) DO M1 = 1, MXE1 K1P = 1 DO K1=MCH,MCI AVL=AVECT1(K1,M1)*AZ AVV=AVECT1(K1,M1)*AZP K2P=1 DO K2=NCH,NCI DVECTV(K2,M1)=DVECTV(K2,M1)+AVV*DKV(K1P,K2P) DVECTL(K2,M1)=DVECTL(K2,M1)+AVL*DKL(K1P,K2P) K2P = K2P + 1 ENDDO K1P = K1P+1 ENDDO ENDDO ENDDO MCH = MCI + 1 MCI = MNP2D1 MCP = MCI - MCH + 1 NCI = 0 DO JK = 1 , IBIN1 NCH = NCI + 1 NCI = NCI + NOTERM READ(1) ((DKL(I,J),J=1,NOTERM),I=1,MCP), & ((DKV(I,J),J=1,NOTERM),I=1,MCP) DO M1 = 1, MXE1 K1P = 1 DO K1=MCH,MCI AVL=AVECT1(K1,M1)*AZ AVV=AVECT1(K1,M1)*AZP K2P=1 DO K2=NCH,NCI DVECTV(K2,M1)=DVECTV(K2,M1)+AVV*DKV(K1P,K2P) DVECTL(K2,M1)=DVECTL(K2,M1)+AVL*DKL(K1P,K2P) K2P = K2P + 1 ENDDO K1P = K1P+1 ENDDO ENDDO ENDDO NCH = NCI + 1 NCI = MNP2D2 NCP = NCI - NCH + 1 READ(1) ((DKL(I,J),J=1,NCP),I=1,MCP), & ((DKV(I,J),J=1,NCP),I=1,MCP) DO M1 = 1, MXE1 K1P = 1 DO K1=MCH,MCI AVL=AVECT1(K1,M1)*AZ AVV=AVECT1(K1,M1)*AZP K2P=1 DO K2=NCH,NCI DVECTV(K2,M1)=DVECTV(K2,M1)+AVV*DKV(K1P,K2P) DVECTL(K2,M1)=DVECTL(K2,M1)+AVL*DKL(K1P,K2P) K2P = K2P + 1 ENDDO K1P = K1P+1 ENDDO ENDDO IF(IBUT.LT.0)GOTO 1000 C C READ BUTTLE PART C NCH=1 NCI=MNP2D2 MCH=MNP2D1+1 MCI=MNP2D1+NCHND1 READ(1)((DKL(I,J),I=1,NCHND1),J=NCH,NCI), & ((DKV(I,J),I=1,NCHND1),J=NCH,NCI) IF(IBUT.GT.0) THEN DO M1 = 1, MXE1 K1P = 1 DO K1=MCH,MCI AVL=AVECT1(K1,M1)*AZ AVV=AVECT1(K1,M1)*AZP K2P=1 DO K2=NCH,NCI DVECTV(K2,M1)=DVECTV(K2,M1)+AVV*DKV(K1P,K2P) DVECTL(K2,M1)=DVECTL(K2,M1)+AVL*DKL(K1P,K2P) K2P = K2P + 1 ENDDO K1P = K1P+1 ENDDO ENDDO ENDIF NCH=MNP2D2+1 NCI=MNP2D2+NCHND2 MCH=1 MCI=MNP2D1 READ(1)((DKL(J,I),J=1,NCHND2),I=MCH,MCI), & ((DKV(J,I),J=1,NCHND2),I=MCH,MCI) IF(IBUT.GT.0) THEN DO M1 = 1, MXE1 K1P = 1 DO K1=MCH,MCI AVL=AVECT1(K1,M1)*AZ AVV=AVECT1(K1,M1)*AZP K2P=1 DO K2=NCH,NCI DVECTV(K2,M1)=DVECTV(K2,M1)+AVV*DKV(K2P,K1P) DVECTL(K2,M1)=DVECTL(K2,M1)+AVL*DKL(K2P,K1P) K2P = K2P + 1 ENDDO K1P = K1P+1 ENDDO ENDDO ENDIF MCH=MNP2D1+1 MCI=MNP2D1+NCHND1 READ(1)((DKL(I,J),J=1,NCHND2),I=1,NCHND1), !NOT J=NCH,NCI - NRB & ((DKV(I,J),J=1,NCHND2),I=1,NCHND1) !NOT J=NCH,NCI - NRB IF(IBUT.GT.0) THEN DO M1 = 1, MXE1 K1P = 1 DO K1=MCH,MCI AVL=AVECT1(K1,M1)*AZ AVV=AVECT1(K1,M1)*AZP K2P=1 DO K2=NCH,NCI DVECTV(K2,M1)=DVECTV(K2,M1)+AVV*DKV(K1P,K2P) DVECTL(K2,M1)=DVECTL(K2,M1)+AVL*DKL(K1P,K2P) K2P = K2P + 1 ENDDO K1P = K1P+1 ENDDO ENDDO ENDIF C 1000 CONTINUE C-----READ ANGULAR COEFFICIENTS READ(1) MAXM1,(CGC(M),M=1,MAXM1) READ(1) ((AMATL(J,I),J=1,NCHND1),I=1,NCHND2), & ((BMATL(J,I),J=1,NCHND1),I=1,NCHND2), & ((BMATV(J,I),J=1,NCHND1),I=1,NCHND2) C DO I2=1,NCHND2 DO I1=1,NCHND1 BMATL(I1,I2)=AZ*BMATL(I1,I2) BMATV(I1,I2)=AZP*BMATV(I1,I2) ENDDO ENDDO C IF(NSPND.NE.0)RETURN C READ(1,END=2000) ((AMATV(J,I),J=1,NCHND1),I=1,NCHND2) RETURN C 2000 DO J=1,NCHND2 DO I=1,NCHND1 AMATV(I,J)=0 ENDDO ENDDO WRITE(6,620) C RETURN C 600 FORMAT(//1X,30(1H*)//' DIMENSION FOR MNP2 OR NCHF TOO' &,' SMALL'// ' MNP2D1, NCHND1 = ',I5,', ',I5/ & ' MNP2D2, NCHND2 = ',I5,', ',I5//1X,30(1H*)//) 601 FORMAT(//1X,30(1H*)//' MZCHF TOO SMALL FOR BUFFER'// & ' NEED AT LEAST ',I4//1X,30(1H*)//) 614 FORMAT(//' *** MISMATCH ***'//15X,'UNIT MFA =',I4, 1 ' HAS IS, IL, IP, MNP2 = ',4I5//) 620 FORMAT(' NO 2ND SET AMAT AVAILABLE: VELOCITY RESULTS UNCERTAIN') C END C C********************************************************************* C SUBROUTINE S1S2(CC,F1,FP1,F2,FP2,FPP2,S1,S2) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C FLIJ=(F1*FP2-FP1*F2)*RZERO FKIJ=CC/(RZERO*RZERO) FIFJ=F1*F2 S1=FLIJ+FIFJ S2=FKIJ*(FLIJ-FIFJ)-2*(F1*FPP2-FP1*FP2) C RETURN END C C********************************************************************* C SUBROUTINE THETBB(R,FNU,LL,T,TP,TD,TDP) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC COMMON/NRBZED/TZED LOGICAL FLAG C M=MZTET 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 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) 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)FNU,LL,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 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) C 610 FORMAT(//5X,55('*')/5X,'SUBROUTINE THETBB'/5X,'LL = ',I3, + ', FNU = ',F10.4/5X,'SUMMATIONS NOT CONVERGED WITH', + I4,' TERMS. ',' VALUE OF FNU TOO LARGE?'/5X,55('*')//) RETURN END C C********************************************************************* C SUBROUTINE TLAG3(MP,FI) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C C CALCULATES INTEGRAL FROM RTWO TO INFINITY OF C F1*(R**MP)*F2 C AND ADDS TO FI C LOGICAL FLAG1(MZCHF),FLAG2(MZCHF),FLAGP,FLAGV,FLAGM LOGICAL FLAG COMMON/CASE/ 1 N1,N2,E1,E2, 2 AVECT1(MNPEXT,MZEST),AVECT2(MNPEXT), 3 FSD1(MZCHF,MZPTS),FSD2(MZCHF,MZPTS), 4 P1(MZCHF),PP1(MZCHF),PPP1(MZCHF), 5 P2(MZCHF),PP2(MZCHF),PPP2(MZCHF), 6 DP1(MZCHF),DPP1(MZCHF),DPPP1(MZCHF), 7 DP2(MZCHF),DPP2(MZCHF),DPPP2(MZCHF), 8 EPS1(MZCHF),EPS2(MZCHF),FNU1(MZCHF),FNU2(MZCHF), 9 FLAG1,FLAG2,FLAGP,FLAGV,FLAGM C COMMON/CB/ 1 IPERT,KP2,KPM,KSLP,MXE1,MXE2,NELC,NZED, 2 H,RTWO,RZERO,WBODE(MZPTS,-3:1), 3 CC1(MZCHF),CC2(MZCHF),ECH1(MZCHF),ECH2(MZCHF) C COMMON/CNTRL/ 1 IPRINT,IBUT,TINY,DELTA,AC C COMMON/NRBZED/TZED C DIMENSION BB(2,MZTET),MSUM(2),FNUS(2),CS(2),PT(2) DIMENSION XLAG(30),WLAG(30) 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/ C C CS(1)=CC1(N1) CS(2)=CC2(N2) FNUS(1)=FNU1(N1) FNUS(2)=FNU2(N2) C C C USES LAGUERRE QUADRATURE WITH NLAG=6 C C CALCULATION OF COEFFICIENTS IN POWER SERIES EXPANSION C ADAPTED FROM SUBROUTINE THETA C DO 40 I=1,2 FNU=FNUS(I) F1=1./FNU F2=F1*F1 X=2.*RTWO*F1 Y=1./X FL=.5*(SQRT(1.+4.*CS(I))-1.) R1=1./RTWO A=TZED*FNU*R1-F1 B=TZED*LOG(RTWO)+RTWO*F2 C=TZED*R1+F2 D=-R1 YN=1. A1=FL-FNU*TZED A2=FL+FNU*TZED+1. BN=-2.*FNU-1. BET=1. GAM=0. W0=C S=1. U=0. CX=0. CY=0. N=0 FLAG=.FALSE. DO 20 MN=1,MZTET 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 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 C RETURN IF FUNCTION ZERO AT RTWO AND BEYOND IF(RTWO*F1.GT.150.)RETURN C C NOT CONVERGED WRITE(6,610)N STOP C C SUMMATIONS CONVERGED 30 PT(I)=S*EXP(-RTWO*F1) IF(TZED.GT.0)PT(I)=PT(I)*RTWO**FNU IF(ABS(PT(I)).LT.1.D-100)RETURN !ZERO FUNCTION AGAIN MSUM(I)=N 40 CONTINUE C C C CALCULATE INTEGRAL C ADAPTED FROM TLAG C FK=1./FNUS(1)+1./FNUS(2) C INITIALISE FOR QUADRATURE NS=NLAG/2 M=NS*(NS-1) M1=M+1 M2=M+NLAG T1=0. C START QUADRATURE DO 140 N=M1,M2 U=XLAG(N) R=RTWO+U/FK C CALCULATE THETA FUNCTIONS C CALCULATION ADAPTED FROM TANDTD DO 150 I=1,2 S=1. FNU=FNUS(I) X=2.*R/FNU Y=1./X AS=1. DO 170 L=1,MSUM(I) AS=AS*Y 170 S=S+BB(I,L)*AS DLR=LOG(R)*TZED IF(I.EQ.1)THEN TI=S TPI=-.5*X+FNU*DLR ELSE TJ=S TPJ=-.5*X+FNU*DLR ENDIF 150 CONTINUE C ADD TO SUM A1=(R**MP)*WLAG(N)*EXP(U+TPI+TPJ) 140 T1=(T1+TI*A1*TJ) C C ADD INCREMENT TO FI FI=FI+T1*FSD1(N1,KP2)*FSD2(N2,KP2)/(PT(1)*PT(2)*FK) C RETURN C 610 FORMAT(//30(1H*)//10X,'SUBROUTINE TLAG3'// 1 'SUMMATIONS NOT CONVERGED WITH ',I3,' TERMS'/ 2 /30(1H*)) END