C N. R. BADNELL UoS v1.9 05/08/05 C C STGD C **** C C C N.B. THIS VERSION REQUIRES STGF BE RUN FIRST TO GENERATE THE ELASTIC C S-MATRICES. IT IS ONLY CONCERNED WITH THE CALCULATION OF DAMPING C CONSTANTS, UNLIKE MJS ORIGINAL WHICH INCORPORATED A STGF RUN. C C PROGRAM MAIN C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (MXTST=(MZTAR*MZTAR+MZTAR)/2) PARAMETER (NTMP=20) C COMMON/CTARG/ENAT(MZTAR),NAST,NELC,NZED,ISAT(MZTAR),LAT(MZTAR) COMMON/CINPUT/ICASE,ISLP(MZSLP) COMMON/CDAMP/SLINE(MXTST),NSLP(MZTAR) COMMON/CMESH/EMESH(MZMSH),MXE,IE,IET(MZTAR),NASTM COMMON/CNTRL/IPRINT,ITOP,LREC COMMON/CTEMP/TEMP(NTMP),FTEMP(NTMP),NTEMP COMMON/NRBSLP/NSPNQ(MZSLP),LRGLQ(MZSLP),NPTYQ(MZSLP) X,ICONAT(MZTAR,MZSLP),LLCHQ(MZCHF,MZSLP),KJQ(MZCHF,MZSLP) X,NCASE(MZSLP,MZTAR) C DIMENSION MSLP(MZSLP) DIMENSION KTA(MXTST),KTB(MXTST),SKINE(MXTST) C NAMELIST/STGD/IPRINT,ITOP,IOPT1,NTEMP C C C I/O UNITS C ********* C C 1 FOR THERMAL AVERAGED OMEGA (ON FILE ZUPSLN) C 2 FOR DIRECTORY FILE S00 C 3 FOR S-MATRIX & SYMMETRY DATA FILE SNM C 5 FOR PROGRAM CONTROL C 6 FOR PRINTED OUTPUT C 7 FOR VALUES OF OMEGA (ON FILE ZOMEGA) C C OPEN(5,FILE='dstgd') OPEN(6,FILE='routd',STATUS='UNKNOWN') OPEN(7,FILE='ZOMEGA',STATUS='UNKNOWN') REWIND(7) C C WRITE HEADING C WRITE(6,6000) C C READ DATA FROM UNIT 5 (NAMELIST/STGD/IPRINT,ITOP,IOPT1,NTEMP) C ********************* C C IPRINT = -1 FOR MINIMUM PRINT (DEFAULT) C = +1 FOR MAXIMUM PRINT C ITOP = 1 FOR TOP-UP C = 0 FOR NO TOP-UP (DEFAULT) C IOPT1 = 1 FOR ALL SLPI CASES (DEFAULT) C = 2 FOR SELECTED SLPI CASES C***FOLLOWED BY C IS, IL, IP FOR CASES SELECTED TERMINATING WITH C -1,-1,-1 C NTEMP = NUMBER OF TEMPERATURES C***FOLLOWED BY (AFTER ANY SYMMETRY INFO) C TEMP(IT) FOR IT=1,NTEMP (Z-SCALED TEMPERATURES) C C IPRINT=-1 ITOP=0 IOPT1=1 NTEMP=0 C READ(5,STGD) WRITE(6,601)IPRINT,ITOP,IOPT1 C C CASE OF IOPT1 = 2 C IF(IOPT1.EQ.2) THEN WRITE(6,602) KSLP=0 C 10 READ(5,*)IS,IL,IP IF(IL.GE.0)THEN KSLP=KSLP+1 IF(KSLP.GT.MZSLP)THEN WRITE(6,605)MZSLP STOP '***DIMENSION EXCEEDED, TOO MANY SYMMETRIES' ENDIF MSLP(KSLP)=10000*IABS(IS)+100*IL+IP IF(IS.LT.0)MSLP(KSLP)=-MSLP(KSLP) WRITE(6,603)MSLP(KSLP) GO TO 10 ENDIF C IF(KSLP.EQ.0)THEN WRITE(6,600) STOP '***ERROR: NO SYMMETRIES SPECIFIED BY IOPT1=2' ENDIF C ENDIF C WRITE(6,604) C C READ TEMPERATURES C IF(NTEMP.GT.NTMP)THEN WRITE(6,666)NTEMP,NTMP STOP ENDIF DO ITEMP=1,NTEMP READ(5,*)TEMP(ITEMP) FTEMP(ITEMP)=157894/TEMP(ITEMP) ENDDO C IF(NTEMP.GT.0)THEN OPEN(1,FILE='ZUPSLN',STATUS='unknown') REWIND(1) ENDIF C C READ S00 DIRECTORY FILE C *********************** C CALL RS00 C WRITE(7,620)NZED,NELC WRITE(7,620)NAST WRITE(7,710)(ENAT(N),N=1,NAST) C C START LOOP ON SLPI CASES C ************************ C IF(IPRINT.LT.0)WRITE(6,670) KASE=0 C DO III=1,ICASE C C FOR IOPT1=2, FIRST CHECK WHETHER THIS SLPI IS REQUIRED C IF(IOPT1.EQ.1)THEN KASE=KASE+1 IF(KASE.GT.MZSLP)THEN WRITE(6,605)MZSLP STOP '***DIMENSION EXCEEDED, TOO MANY SYMMETRIES' ENDIF MSLP(KASE)=ISLP(III) ELSE DO K=KASE+1,KSLP IF(ISLP(III).EQ.MSLP(K))THEN KASE=KASE+1 M=MSLP(KASE) MSLP(KASE)=ISLP(III) MSLP(K)=M GO TO 40 ENDIF ENDDO GO TO 2000 ENDIF C C READ S-MATRIX DATA FOR THIS SLPI CASE C ************************************* C 40 CALL RSFILE(III) C C SET INDICES AND WRITE ZS TO SCRATCH FILES C NCH2=0 DO IT=1,NASTM IF(ICONAT(IT,III).GT.0)THEN NCH1=NCH2+1 NCH2=NCH2+ICONAT(IT,III) NSLP(IT)=NSLP(IT)+1 N=NSLP(IT) NCASE(N,IT)=III ENDIF ENDDO C C **************** C 2000 ENDDO C C END OF SLPI LOOP C **************** C IF(IOPT1.EQ.2)THEN IF(KSLP.NE.KASE)THEN ! SYMMETRY CASES NOT FOUND WRITE(6,660) DO K=KASE+1,KSLP WRITE(6,661)MSLP(K) ENDDO ENDIF ENDIF KSLP=KASE C C C PRINT SUMMARY OF RESULTS THUS FAR C ********************************* C IF(IPRINT.GT.0)THEN WRITE(6,695) WRITE(6,699) DO IT=1,NASTM DO N=1,NSLP(IT) III=NCASE(N,IT) WRITE(6,698)IT,NSPNQ(III),LRGLQ(III),NPTYQ(III) ENDDO ENDDO ENDIF C C CALCULATE ZOMEGA C **************** C K=0 KMM=0 DO IT=2,NASTM DO JT=1,IT-1 K=K+1 IF(SLINE(K).GT.0.)THEN KMM=KMM+1 KTA(KMM)=IT KTB(KMM)=JT SKINE(KMM)=SLINE(K) ENDIF ENDDO ENDDO WRITE(7,620)KMM C C INITIAL WRITES TO ZPSLN C WRITE(1,1001)NZED,NELC,NAST,NTEMP,KMM DO I=1,NAST WRITE(1,1002)I,ISAT(I),LAT(I),ENAT(I) ENDDO C WRITE(6,695) WRITE(6,610) DO K=1,KMM IT=KTA(K) JT=KTB(K) WRITE(6,611)JT,IT,SKINE(K) WRITE(1,1003)JT,IT,SKINE(K) ENDDO WRITE(6,695) C DO NT=1,NTEMP WRITE(1,690)TEMP(NT) ENDDO C C *** COMPUTE DAMPING CONSTANTS C DO K=1,KMM IT=KTA(K) JT=KTB(K) C CALL OMDAMP(JT,IT,SKINE(K)) C ENDDO C C**** C CLOSE(7) CLOSE(1,STATUS='KEEP') WRITE(6,612) STOP 'STGD: NORMAL END' C C FORMATS C ******* C 600 FORMAT( //10X,'READS IOPT1 = 2 FOLLOWED BY TERMINATOR'//) 601 FORMAT(//5X,'DATA READ FROM UNIT 5'// 1 10X,'IPRINT = ',I2/ 3 10X,'ITOP = ',I2/ 6 10X,'IOPT1 = ',I2//) 602 FORMAT(/5X,'VALUES OF 1000*IS+100*IL+IP READ FOR IOPT1=2') 603 FORMAT(52X,I7) 604 FORMAT(/) 605 FORMAT(///20X,'TOO MANY SYMMETRIES'// 2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS MZSLP = ',I3//) 610 FORMAT(25X,'ALLOWED TRANSITIONS'// + 18X,'JT',8X,'IT',5X,'LINE STRENGTH'/) 611 FORMAT( 10X,2I10,1PE15.3) 612 FORMAT(//10X,23('*')/10X,'* FILE ZOMEGA WRITTEN *', + /10X,23('*')//) 620 FORMAT(2I5) 650 FORMAT(///1X,70('*')/1X,70('*')//20X,'ENERGIES AND TOTAL OMEGAS'/ 1 20X,25('*')/) 660 FORMAT(//10X,30('*')/ 1 10X,'NO DATA ON S-MATRIX FILE FOR'/ 2 10X,'10000*IS+100*IL+IP = '/) 661 FORMAT(10X,I10) 666 FORMAT(//1X,'***** NTEMP = ',I2,' IS LARGER THAN', + ' MAXIMUM VALUE OF ',I2,' ALLOWED BY DIMENSIONS *****'//) 670 FORMAT(//10X,'RUN WITH MINIMUM PRINT'//) 690 FORMAT(1PE12.2) 695 FORMAT(///1X,70('*')/1X,70('*')//) 698 FORMAT(13X,I2,9X,I2,5X,I2,5X,I2,11X,I2,5X,I2) 699 FORMAT(25X,' S-MATRICES READ FOR'// 1 12X,'NTARG',8X,'IS',3X,'IL/2J',4X,'IP'/) 6000 FORMAT(//30X,'PROGRAM STGD'/30X,12('*')// + 5 X,'CALCULATES COLLISION STRENGTHS OMEGA FOR', + ' ELECTRON-IMPACT LINE-BROADENING'//) C 710 FORMAT(5E16.6) 1001 FORMAT(5I5) 1002 FORMAT(3I5,1PE15.6) 1003 FORMAT(2I5,1PE15.3) C END C C********************************************************** C SUBROUTINE MERGE(ITA,ITB,KMAX,KPA,KPB) C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER(NMSH=2*MZMSH) C COMMON/CTARG/ENAT(MZTAR),NAST,NELC,NZED,ISAT(MZTAR),LAT(MZTAR) COMMON/CMESH/EMESH(MZMSH),MXE,IE,IET(MZTAR),NASTM COMMON/CZOM/EP(NMSH),ZOMEGA(NMSH),ZTOP1(NMSH),ZTOP2(NMSH) C DIMENSION EPA(NMSH),EPB(NMSH),KPA(NMSH),KPB(NMSH) C C INITIALISATION EA=ENAT(ITA) EB=ENAT(ITB) IETA=IET(ITA) IETB=IET(ITB) DO IE=IETA,MXE EPA(IE)=EMESH(IE)-EA ENDDO DO IE=IETB,MXE EPB(IE)=EMESH(IE)-EB ENDDO C C C SORT MERGED LEVELS IN ENERGY-ORDER IA=IETA IB=IETB K=0 1010 K=K+1 IF(EPA(IA).LT.EPB(IB))THEN KPA(K)=IA KPB(K)=0 EP(K)=EPA(IA) IF(IA.EQ.MXE)GOTO 1011 IA=IA+1 ELSE KPA(K)=0 KPB(K)=IB EP(K)=EPB(IB) IF(IB.EQ.MXE)GOTO 1011 IB=IB+1 ENDIF GOTO 1010 1011 KMAX=K C C FIND NEAREST NEIGHBOURS IA=IETA IB=IETB DO 1020 K=1,KMAX E=EP(K) IF(KPA(K).NE.0)THEN DO 1012 I=IB,MXE IF(EPB(I).GT.E)THEN IBB=I GOTO 1013 ENDIF 1012 CONTINUE 1013 IF(IBB.EQ.IETB)THEN KPB(K)=IB GOTO 1020 ENDIF IF(ABS(E-EPB(IBB-1)).LT.ABS(E-EPB(IBB)))THEN KPB(K)=IBB-1 IB=IBB-1 ELSE KPB(K)=IBB IB=IBB ENDIF ELSE DO 1014 I=IA,MXE IF(EPA(I).GT.E)THEN IAA=I GOTO 1015 ENDIF 1014 CONTINUE 1015 IF(IAA.EQ.IETA)THEN KPA(K)=IA GOTO 1020 ENDIF IF(ABS(E-EPA(IAA-1)).LT.ABS(E-EPA(IAA)))THEN KPA(K)=IAA-1 IA=IAA-1 ELSE KPA(K)=IAA IA=IAA ENDIF ENDIF 1020 CONTINUE C C ELIMINATE DUPLICATE PAIRS KM=KMAX DO 2010 K=2,KMAX IF(K.GT.KM)GOTO 2020 IF(KPA(K-1).EQ.KPA(K).AND.KPB(K-1).EQ.KPB(K))THEN EP(K-1)=.5*(EP(K-1)+EP(K)) KM=KM-1 DO 2000 J=K,KM EP(J)=EP(J+1) KPA(J)=KPA(J+1) KPB(J)=KPB(J+1) 2000 CONTINUE ENDIF 2010 CONTINUE 2020 KMAX=KM C C RETURN END C C********************************************************** SUBROUTINE OMDAMP(ITA,ITB,SKINE) C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C LOGICAL TOP C INCLUDE 'PARAM' C PARAMETER (MXDEG=8) PARAMETER (MXCST=(MXDEG*(MXDEG+1))/2) PARAMETER (MXTST=(MZTAR*MZTAR+MZTAR)/2) PARAMETER (NMSH=2*MZMSH) C PARAMETER (TZERO=0.0) PARAMETER (ZERO=DCMPLX(TZERO,TZERO)) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CHARACTER NAME*3,NUM(0:9)*1 C COMMON/CTARG/ENAT(MZTAR),NAST,NELC,NZED,ISAT(MZTAR),LAT(MZTAR) COMMON/CDAMP/SLINE(MXTST),NSLP(MZTAR) COMMON/CMESH/EMESH(MZMSH),MXE,IE,IET(MZTAR),NASTM COMMON/CNTRL/IPRINT,ITOP,LREC COMMON/CZOM/EP(NMSH),ZOMEGA(NMSH),ZTOP1(NMSH),ZTOP2(NMSH) COMMON/NRBSLP/NSPNQ(MZSLP),LRGLQ(MZSLP),NPTYQ(MZSLP) X,ICONAT(MZTAR,MZSLP),LLCHQ(MZCHF,MZSLP),KJQ(MZCHF,MZSLP) X,NCASE(MZSLP,MZTAR) C DIMENSION ZSA(MZMSH,MXCST),ZSB(MZMSH,MXCST),ZINC(NMSH) X,KPA(NMSH),KPB(NMSH) C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C IF(IPRINT.GE.0)WRITE(6,652)ITA,ITB,SKINE C C CHECK OK FOR TOP-UP IN L C CALL TOPOK(ITOP,ITA,ITB,TOP,LRGLX,LRGLN) IF(ITOP.NE.1)TOP=.FALSE. C C SET UP MERGED MESH C CALL MERGE(ITA,ITB,KMAX,KPA,KPB) C C DO K=1,KMAX ZTOP1(K)=ZERO ZTOP2(K)=ZERO ZOMEGA(K)=ZERO ENDDO C C ADJUST WEIGHT FACTORS FOR LS- OR JK-COUPLING C LTA=LAT(ITA) LTB=LAT(ITB) IF(NSPNQ(1).NE.0)THEN SA=ONE/DBLE(ISAT(ITA)) IDIP=1 IF(IPRINT.GE.0)WRITE(6,610) ELSE SA=ONE IDIP=2 IF(IPRINT.GE.0)WRITE(6,611) ENDIF C C LOOP OVER SYMMETRIES C c write(*,*)ita,itb DO NA=1,NSLP(ITA) C IIA=NCASE(NA,ITA) NSPNA=NSPNQ(IIA) LRGLA=LRGLQ(IIA) NPTYA=NPTYQ(IIA) NCHA=0 DO I=1,ITA-1 NCHA=NCHA+ICONAT(I,IIA) ENDDO KM=(ICONAT(ITA,IIA)*(ICONAT(ITA,IIA)+1))/2 C NAME='S'//NUM(IIA/10)//NUM(IIA-10*(IIA/10)) OPEN(3,FILE=NAME,STATUS='OLD',ACCESS='DIRECT',RECL=LREC) NREC=2 DO I=1,ITA-1 IF(ICONAT(I,IIA).GT.0)NREC=NREC+MXE-IET(I)+1 ENDDO DO IE=IET(ITA),MXE NREC=NREC+1 READ(3,REC=NREC)(ZSA(IE,K),K=1,KM) ENDDO C DO 100 NB=1,NSLP(ITB) C IIB=NCASE(NB,ITB) NSPNB=NSPNQ(IIB) LRGLB=LRGLQ(IIB) NPTYB=NPTYQ(IIB) C IF(NSPNB.NE.NSPNA.OR. 1 NPTYB.EQ.NPTYA.OR. 2 IABS(LRGLB-LRGLA).GT.IDIP.OR. 3 LRGLA.EQ.0.AND.LRGLB.EQ.0) GOTO 100 C IF(top.and.LRGLA.EQ.LRGLN.AND.LRGLB.EQ.LRGLN)GOTO 100 c write(*,*)'iia, iib: ',iia,iib C NCHB=0 DO I=1,ITB-1 NCHB=NCHB+ICONAT(I,IIB) ENDDO KM=(ICONAT(ITB,IIB)*(ICONAT(ITB,IIB)+1))/2 C NAME='S'//NUM(IIB/10)//NUM(IIB-10*(IIB/10)) OPEN(4,FILE=NAME,STATUS='OLD',ACCESS='DIRECT',RECL=LREC) NREC=2 DO I=1,ITB-1 IF(ICONAT(I,IIB).GT.0)NREC=NREC+MXE-IET(I)+1 ENDDO DO IE=IET(ITB),MXE NREC=NREC+1 READ(4,REC=NREC)(ZSB(IE,K),K=1,KM) ENDDO C C CONTINUE WITH CALCULATION OF OMEGA C IF(NSPNA.GT.0)THEN Q0=SA*NSPNA*(2*LRGLA+1)*(2*LRGLB+1) ELSEIF(NSPNA.LT.0)THEN Q0=2*(2*LRGLA+1)*(2*LRGLB+1) ELSE Q0=SA*(LRGLA+1)*(LRGLB+1) ENDIF C DO K=1,KMAX ZINC(K)=ZERO ENDDO C C LOOP OVER ALL POSSIBLE CHANNEL INTERACTIONS C KSA=0 DO IA=1,ICONAT(ITA,IIA) LA=LLCHQ(NCHA+IA,IIA) KA=KJQ(NCHA+IA,IIA) DO IAP=IA,ICONAT(ITA,IIA) KSA=KSA+1 LAP=LLCHQ(NCHA+IAP,IIA) KAP=KJQ(NCHA+IAP,IIA) C KSB=0 DO IB=1,ICONAT(ITB,IIB) LB=LLCHQ(NCHB+IB,IIB) IF(LA.NE.LB)THEN KSB=KSB+ICONAT(ITB,IIB)-IB+1 GO TO 110 ENDIF KB=KJQ(NCHB+IB,IIB) L=LA DO IBP=IB,ICONAT(ITB,IIB) KSB=KSB+1 LBP=LLCHQ(NCHB+IBP,IIB) IF(LAP.NE.LBP)GO TO 120 KBP=KJQ(NCHB+IBP,IIB) LP=LAP c write(6,*)'channels:',ia,iap,ib,ibp C IF(IA.EQ.IAP.AND.IB.EQ.IBP)THEN !DIAGONAL IF(NSPNA.NE.0)THEN !LS-COUPLING Q=Q0*WABS1(LTA,LTB,LRGLA,LRGLB,L)**2 !L.EQ.LP ELSE IF(KA+KB.EQ.0)GO TO 120 Q=Q0*(KA+1)*(KB+1) Q=Q *WABS2(LTA,LTB,KA,KB,2*L)**2 !L.NE.LP Q=Q *WABS2(KA,KB,LRGLA,LRGLB,1)**2 !K.NE.KP ENDIF DO K=1,KMAX ZINC(K)=ZINC(K)+Q*(ONE-ZSA(KPA(K),KSA) X *CONJG(ZSB(KPB(K),KSB))) c write(6,671)k,zinc(k) ENDDO ELSE !OFF-DIAG IF(NSPNA.NE.0)THEN !LS-COUPLING Q=Q0*WABS1(LTA,LTB,LRGLA,LRGLB,L) !L.NE.LP X *WABS1(LTA,LTB,LRGLA,LRGLB,LP) ELSE IF(KA+KB.EQ.0.OR.KAP+KBP.EQ.0)GO TO 120 Q=Q0*SQRT(DBLE((KA+1)*(KB+1)*(KAP+1)*(KBP+1))) Q=Q *WABS2(LTA,LTB,KA,KB,2*L) !L.NE.LP X *WABS2(LTA,LTB,KAP,KBP,2*LP) Q=Q *WABS2(KA,KB,LRGLA,LRGLB,1) !K.NE.KP X *WABS2(KAP,KBP,LRGLA,LRGLB,1) ENDIF DO K=1,KMAX ZINC(K)=ZINC(K)-TWO*Q*ZSA(KPA(K),KSA) X *CONJG(ZSB(KPB(K),KSB)) c write(6,671)k,zinc(k) c 671 format(i5,1pe12.3,e12.3) ENDDO ENDIF C IF(IPRINT.GE.0)THEN IF(NSPNA.NE.0)THEN WRITE(6,620)NSPNA,NPTYA,LRGLA,LRGLB,L,LP,Q ELSE WRITE(6,621)NPTYA,LRGLA,LRGLB,KA,KAP,KB,KBP,L,LP,Q ENDIF ENDIF C 120 ENDDO 110 ENDDO ENDDO ENDDO C DO K=1,KMAX IF(DBLE(ZINC(K)).LT.0)WRITE(6,699)DBLE(ZINC(K)),ITA,ITB, + LRGLA,LRGLB,K ZOMEGA(K)=ZOMEGA(K)+ZINC(K) ENDDO C IF(TOP)THEN IF(LRGLA.EQ.LRGLX.OR.LRGLB.EQ.LRGLX)THEN DO K=1,KMAX ZTOP2(K)=ZTOP2(K)+ZINC(K) ENDDO ELSEIF(LRGLA.EQ.LRGLX-1*IDIP.OR.LRGLB.EQ.LRGLX-1*IDIP)THEN DO K=1,KMAX ZTOP1(K)=ZTOP1(K)+ZINC(K) ENDDO ENDIF ENDIF C CLOSE(4) 100 ENDDO C CLOSE(3) C C END LOOP OVER SYMMETRIES C ENDDO C C THERMAL AVERAGE C CALL THEAV(ITA,ITB,KMAX,TOP) C C TOP-UP FOR OUTPUT ZOMEGA FILE C IF(TOP)THEN ETHR=ENAT(NAST)-MIN(ENAT(ITA),ENAT(ITB)) WRITE(6,660)ETHR C DO K=1,KMAX IF(EP(K).GE.ETHR)THEN ZOMTOP=ZOMEGA(K)+DCMPLX((DBLE(ZTOP2(K))**2)/ + DBLE(ZTOP1(K)-ZTOP2(K)),(DIMAG(ZTOP2(K))**2)/ + DIMAG(ZTOP1(K)-ZTOP2(K))) WRITE(6,661)K,EP(K),DBLE(ZOMEGA(K)),DBLE(ZTOP1(K)) + ,DBLE(ZTOP2(K)),DBLE(ZOMTOP) ENDIF ENDDO ENDIF C C PRINT RESULTS C IF(IPRINT.GE.0)WRITE(6,630) WRITE(7,700)ITA,ITB,KMAX,SLINE(ITA+((ITB-1)*(ITB-2))/2) DO K=1,KMAX WRITE(7,640)EP(K),ZOMEGA(K) IF(IPRINT.GE.0)WRITE(6,640)EP(K),ZOMEGA(K) ENDDO C RETURN C 610 FORMAT(25X,'ZOM CALCULATED FOR'// 1 5X,'NSPNA',5X,'NPTYA',5X,'LRGLA',5X,'LRGLB', 2 7X,'L',5X,'LP',9X,'Q'/) 611 FORMAT(25X,'ZOM CALCULATED FOR'// 1 5X,'NPTYA',5X,'LRGLA',5X,'LRGLB', 1 5X,'KA',4X,'KAP',5X,'KB',4X,'KBP', 2 7X,'L',5X,'LP',9X,'Q'/) 620 FORMAT(7X,I2,8X,I2,8X,I2,8X,I2,7X,I2,5X,I2,4X,F12.6) 621 FORMAT(7X,I2,8X,I2,8X,I2,6X,I2,5X,I2,5X,I2,5X,I2,6X,I2,5X,I2, X 4X,F12.6) 630 FORMAT(//6X,'EPS',7X,'ZOMEGA'//) 652 FORMAT(1X,69(1H*)// +10X,'TARGET STATES : -'/10X,'ITA = ',I5/10X,'ITB = ',I5/ + 10X,'LINE STRENGTH = ',1PE12.3//) 660 FORMAT(//10X,'ALL OPEN, E(THRESHOLD) = ',1PE12.3// + 9X, 'K EP(K) REAL(ZOM) REAL(ZTOP1) REAL(ZTOP2)', + ' REAL(ZOMTOP)'//) 661 FORMAT(I10,1PE12.4,4E12.3) 700 FORMAT(2I5,I10,1PE12.4) 640 FORMAT(0p,F12.6,1P,2E12.4) 699 FORMAT(//' ::::: REAL(ZINC).LT.0 ::::: REAL(ZINC) = ',1PE9.1/ + ' ::::: ITA = ',I3,', ITB = ',I3,', LRGLA = ',I3,', LRGLB = ',I3 + ,', K = ',I5) C END C C*************************************************************** C SUBROUTINE RS00 C C READS DATA INDEPENDENT OF SLPI FROM DIRECTORY FILE S00 C C THE FOLLOWING DATA ARE READ C NZED = NUCLEAR CHARGE C NELC = NUMBER OF ELECTRONS IN TARGET C NAST = NUMBER OF TARGET STATES C FOR I = 1,NAST - C ENAT(I) = TARGET ENERGIES C LAT(I) = TARGET ORBITAL ANGULAR MOMENTA C ISAT(I) = VALUES OF (2*S+1) FOR TARGET STATES C SLINE = (DIPOLE) LINE STRENGTHS BETWEEN TARGET STATES C THEN ALL SLP CASES: C LRGL2 = TOTAL ORBITAL ANGULAR MOMENTUM C NSPN2 = TOTAL (2*S+1) C NPTY2 = TOTAL PARITY C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (MXTST=(MZTAR*MZTAR+MZTAR)/2) C COMMON/CTARG/ENAT(MZTAR),NAST,NELC,NZED,ISAT(MZTAR),LAT(MZTAR) COMMON/CINPUT/ICASE,ISLP(MZSLP) COMMON/CDAMP/SLINE(MXTST),NSLP(MZTAR) COMMON/CMESH/EMESH(MZMSH),MXE,IE,IET(MZTAR),NASTM COMMON/CNTRL/IPRINT,ITOP,LREC COMMON/NRBSLP/NSPNQ(MZSLP),LRGLQ(MZSLP),NPTYQ(MZSLP) X,ICONAT(MZTAR,MZSLP),LLCHQ(MZCHF,MZSLP),KJQ(MZCHF,MZSLP) X,NCASE(MZSLP,MZTAR) C OPEN(2,FILE='S00',STATUS='OLD',FORM='UNFORMATTED') REWIND(2) C C READ TARGET INFO C READ(2)NZED,NELC,NAST,IONE IF(NAST.GT.MZTAR)THEN WRITE(6,610)NAST,MZTAR STOP '***INCREASE MZTAR' ENDIF C READ(2)(ENAT(I),I=1,NAST),(ISAT(I),LAT(I),I=1,NAST) X ,(SLINE(I),I=1,((NAST-IONE)*(NAST+1+IONE))/2) C C WRITE TARGET PROPERTIES C WRITE(6,650)NZED,NELC WRITE(6,655) DO J=1,NAST WRITE(6,660)J,ISAT(J),LAT(J),ENAT(J) ENDDO C C IF SLINE ALLOWS-FOR ELASTIC TRANSITIONS THEN REMAP. C IF(IONE.EQ.0)THEN K=0 DO I=2,NAST IM=I-1 DO J=1,IM K=K+1 SLINE(K)=SLINE(K+IM) ENDDO ENDDO ENDIF C C READ ENERGY MESH INFO C READ(2)MXE IF(MXE.GT.MZMSH)THEN WRITE(6,630)MXE,MZMSH STOP '*** INCREASE DIMENSION MZMSH' ENDIF READ(2)(EMESH(I),I=1,MXE) C C CALCULATE NASTM=NUMBER OF OPEN TARGET LEVELS C CALCULATE IET(IT)=STARTING MESH POINT FOR TARGET IT C INITIALISE NSLP(IT) C NASTM=NAST DO IT=1,NAST DO IE=1,MXE IF(EMESH(IE).GE.ENAT(IT))THEN IET(IT)=IE GO TO 913 ENDIF ENDDO NASTM=IT-1 GO TO 915 913 NSLP(IT)=0 ENDDO 915 CONTINUE C C READ AVAILABLE SYMMETRIES C ICASE=0 920 READ(2)NSPN2,LRGL2,NPTY2 IF(LRGL2.LT.0)THEN LREC=NPTY2 RETURN ENDIF C ICASE=ICASE+1 IF(ICASE.GT.MZSLP)THEN WRITE(6,605)MZSLP STOP '***TOO MANY SYMMETRIES' ENDIF ISLP(ICASE)=10000*IABS(NSPN2)+100*LRGL2+NPTY2 IF(NSPN2.LT.0)ISLP(ICASE)=-ISLP(ICASE) NSPNQ(ICASE)=NSPN2 LRGLQ(ICASE)=LRGL2 NPTYQ(ICASE)=NPTY2 GO TO 920 C 605 FORMAT(///20X,'TOO MANY SYMMETRIES'// 2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS MZSLP = ',I3//) 610 FORMAT(///20X,'TOO MANY TARGET STATES'// 1 10X,'VALUE READ FOR NAST IS',I3// 2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS MZTAR = ',I3//) 630 FORMAT(///10X,'MXE = ',I7,' IS GREATER THAN', + ' MZMSH = ',I7///) 650 FORMAT(/10X,'NUCLEAR CHARGE =',I3,', NUMBER OF TARGET ', 1 ' ELECTRONS =',I3/10X,53('*')//) 655 FORMAT(20X,'TARGET STATES -'/20X,15('*')// 1 10X,'INDEX',3X,'(2*S+1)/P',4X,'L/2J',7X,'SCALED ENERGY'/) 660 FORMAT(3X,3I10,7X,F12.6) C END C C*************************************************************** C SUBROUTINE RSFILE(KASE) C C READS S-MATRIX DATA AND SYMMETRY INFO FOR ONE SLPI CASE C C THE FOLLOWING DATA ARE READ - C NCHAN = NUMBER OF CHANNELS C FOR I = 1,NAST - C NCONAT(I) = NUMBER OF CHANNELS FOR TARGET STATE I C FOR I = 1,NCHAN - C L2P(I) = SMALL L FOR CHANNEL I C KJ(I) = K FOR CHANNEL I (JK-COUPLING ONLY) C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (MXDEG=8) PARAMETER (MXCST=(MXDEG*(MXDEG+1))/2) PARAMETER (MXTST=(MZTAR*MZTAR+MZTAR)/2) C CHARACTER NAME*3,NUM(0:9)*1 C COMMON/CTARG/ENAT(MZTAR),NAST,NELC,NZED,ISAT(MZTAR),LAT(MZTAR) COMMON/CMESH/EMESH(MZMSH),MXE,IE,IET(MZTAR),NASTM COMMON/CNTRL/IPRINT,ITOP,LREC COMMON/NRBSLP/NSPNQ(MZSLP),LRGLQ(MZSLP),NPTYQ(MZSLP) X,ICONAT(MZTAR,MZSLP),LLCHQ(MZCHF,MZSLP),KJQ(MZCHF,MZSLP) X,NCASE(MZSLP,MZTAR) C DIMENSION ZS(MXCST) DIMENSION NCONAT(MZTAR),L2P(MZCHF),KJ(MZCHF),ITARG(MZCHF) C DATA NUM /'0','1','2','3','4','5','6','7','8','9'/ C NSPN2=NSPNQ(KASE) LRGL2=LRGLQ(KASE) NPTY2=NPTYQ(KASE) C NAME='S'//NUM(KASE/10)//NUM(KASE-10*(KASE/10)) OPEN(3,FILE=NAME,STATUS='OLD',ACCESS='DIRECT',RECL=LREC) C READ(3,REC=1)(NCONAT(I),I=1,NAST),NCHAN IF(NCHAN.GT.MZCHF)THEN WRITE(6,500)NSPN2,LRGL2,NPTY2,NCHAN,MZCHF STOP '***INCREASE MZCHF' ENDIF READ(3,REC=2)(L2P(I),I=1,NCHAN),(KJ(I),I=1,NCHAN) C DO I=1,NAST IF(NCONAT(I).GT.MXDEG)THEN WRITE(6,640)NCONAT(I),MXDEG STOP '*** TOO MANY CHANNELS PER TARGET' ENDIF ICONAT(I,KASE)=NCONAT(I) ENDDO DO I=1,NCHAN LLCHQ(I,KASE)=L2P(I) KJQ(I,KASE)=KJ(I) ENDDO C C CHANNELS C C C WRITE SLPI AND CHANNEL DATA C IF(IPRINT.GE.0)THEN I=0 DO J=1,NAST K=NCONAT(J) IF(K.GT.0)THEN DO L=1,K I=I+1 ITARG(I)=J ENDDO ENDIF ENDDO IF(NSPN2.NE.0)THEN WRITE(6,600)NSPN2,LRGL2,NPTY2 WRITE(6,669) WRITE(6,670)(I,ITARG(I),L2P(I),I=1,NCHAN) ELSE WRITE(6,600)NSPN2,LRGL2,NPTY2 WRITE(6,668) WRITE(6,671)(I,ITARG(I),L2P(I),KJ(I),I=1,NCHAN) ENDIF ELSE WRITE(6,680)NSPN2,LRGL2,NPTY2 ENDIF C C READ S-MATRIX C IF(IPRINT.GT.0)THEN DO IE=1,MXE C WRITE(6,707)EMESH(IE) WRITE(6,659) NCH2=0 NREC0=2 DO IT=1,NASTM IF(IE.LT.IET(IT))GO TO 400 IF(NCONAT(IT).GT.0)THEN NREC=NREC0+IE-IET(IT)+1 M=(NCONAT(IT)*(NCONAT(IT)+1))/2 READ(3,REC=NREC)(ZS(K),K=1,M) NCH1=NCH2+1 NCH2=NCH2+NCONAT(IT) K=0 DO I=NCH1,NCH2 DO J=I,NCH2 K=K+1 WRITE(6,660)IT,I,J,ZS(K) ENDDO ENDDO NREC0=NREC0+MXE-IET(IT)+1 ENDIF ENDDO C 400 ENDDO ENDIF C CLOSE(3) C 500 FORMAT(///20X,'TOO MANY CHANNELS FOR (IS, IL, IP) = (', 1 3I3,')'//10X,'VALUE READ FOR NCHAN IS',I5// 2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS MZCHF =',I5//) 600 FORMAT(///1X,70('+')//15X,'S L/2J P =',3I3/15X,19('*')//) 640 FORMAT(///20X,'TOO MANY CHANNELS PER TARGET: ',I3,' EXCEEDS MXDEG' X,I3) 659 FORMAT(//' ........ ELEMENTS OF S-MATRIX'/) 660 FORMAT(10X,'IT = ',I3,', I = ',I3,', J = ',I3,', ZS = ', + 1P,2E12.4) 668 FORMAT(12X,'CHANNEL',2X,'TARGET',4X,'SMALL', 1 /12X,'INDEX ',3X,'INDEX ',5X,'L',5X,'K'/) 669 FORMAT(12X,'CHANNEL',2X,'TARGET',4X,'SMALL', 1 /12X,'INDEX ',3X,'INDEX ',5X,'L'/) 670 FORMAT(7X,I8,I9,I10) 671 FORMAT(7X,I8,I9,I10,I6) 680 FORMAT(15X,'S L/2J P = ',3I3) 707 FORMAT(//5X,'ETOT = ',F10.6/5X,17('*')) C RETURN END C C*********************************************************************** C SUBROUTINE THEAV(ITA,ITB,KMAX,TOP) C C C CALCULATES THERMALL-AVERAGED COLLISION STRENGTHS, ZUPSILON, C FOR ELECTRON-IMPACT LINE-BROADENING C C TOP-UP IN E C ZUPINC IS ESTIMATED CONTRIBUTION FROM EP.GT.EP(KMAX) C TOP-UP IN L C ZUP1 IS CONTRIBUTION FROM C LMAX-LMAX,LMAX-(LMAX-1),(LMAX-1)-LMAX C (LMAX-1),(LMAX-1);(LMAX-1),(LMAX-2);(LMAX-2),(LMAX-1) C ZUP2 IS CONTRIBUTION FROM C LMAX,LMAX; LMAX,(LMAX-1);(LMAX-1),LMAX C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (NMSH=2*MZMSH) PARAMETER (NTMP=20) C LOGICAL TOP C COMMON/CTEMP/TEMP(NTMP),FTEMP(NTMP),NT COMMON/CZOM/EP(NMSH),ZOMEGA(NMSH),ZTOP1(NMSH),ZTOP2(NMSH) DIMENSION EX1(NTMP),EX2(NTMP),EX3(NTMP) + , ZUPS(NTMP),ZUP1(NTMP),ZUP2(NTMP),ZUPINC(NTMP) C IF(NT.EQ.0)RETURN C C WRITE HEADINGS WRITE(6,601)ITA,ITB WRITE(1,1001)ITA,ITB C IF(KMAX.LT.3)THEN WRITE(6,640)KMAX RETURN ENDIF C IF(NT.GT.NTMP)THEN WRITE(6,660) NT=NTMP ENDIF C C INTEGRATION OF ZOMEGA C C FIRST TWO ENERGIES E1=EP(1) E2=EP(2) Z1=ZOMEGA(1) Z2=ZOMEGA(2) C C CONTINUE INTEGRATIONS DO IT=1,NT EX1(IT)=EXP(-E1*FTEMP(IT)) EX2(IT)=EXP(-E2*FTEMP(IT)) EXX=EX1(IT)-EX2(IT) ZUPS(IT)=Z1*EXX ENDDO DO K=3,KMAX E3=EP(K) Z3=ZOMEGA(K) DO IT=1,NT EX3(IT)=EXP(-E3*FTEMP(IT)) EXX=EX1(IT)-EX3(IT) ZUPS(IT)=ZUPS(IT)+Z2*EXX ENDDO E1=E2 E2=E3 Z1=Z2 Z2=Z3 DO IT=1,NT EX1(IT)=EX2(IT) EX2(IT)=EX3(IT) ENDDO ENDDO C C LAST POINT DO IT=1,NT ZUPINC(IT)=Z3*EX3(IT) EXX=EX1(IT)+EX3(IT) ZUPS(IT)=.5*(ZUPS(IT)+Z3*EXX) ENDDO C C TOP-UP IN L C IF(TOP)THEN C C FIRST TWO ENERGIES E1=EP(1) E2=EP(2) ZT11=ZTOP1(1) ZT21=ZTOP2(1) ZT12=ZTOP1(2) ZT22=ZTOP2(2) C C CONTINUE INTEGRATIONS DO IT=1,NT EX1(IT)=EXP(-E1*FTEMP(IT)) EX2(IT)=EXP(-E2*FTEMP(IT)) EXX=EX1(IT)-EX2(IT) ZUP1(IT)=ZT11*EXX ZUP2(IT)=ZT21*EXX ENDDO DO K=3,KMAX E3=EP(K) ZT13=ZTOP1(K) ZT23=ZTOP2(K) DO IT=1,NT EX3(IT)=EXP(-E3*FTEMP(IT)) EXX=EX1(IT)-EX3(IT) ZUP1(IT)=ZUP1(IT)+ZT12*EXX ZUP2(IT)=ZUP2(IT)+ZT22*EXX ENDDO E1=E2 E2=E3 ZT11=ZT12 ZT12=ZT13 ZT21=ZT22 ZT22=ZT23 DO IT=1,NT EX1(IT)=EX2(IT) EX2(IT)=EX3(IT) ENDDO ENDDO C C LAST POINT DO IT=1,NT EXX=EX1(IT)+EX3(IT) ZUP1(IT)=.5*(ZUP1(IT)+ZT13*EXX) ZUP2(IT)=.5*(ZUP2(IT)+ZT23*EXX) ENDDO C C ADD TOP-UP DO IT=1,NT A=(DBLE(ZUP2(IT))**2)/(DBLE(ZUP1(IT))-DBLE(ZUP2(IT))) B=(DIMAG(ZUP2(IT))**2)/(DIMAG(ZUP1(IT))-DIMAG(ZUP2(IT))) ZUPS(IT)=ZUPS(IT)+DCMPLX(A,B) ENDDO C ENDIF C C PRINT RESULTS DO IT=1,NT WRITE(6,670)TEMP(IT),ZUPS(IT),ZUPINC(IT) WRITE(1,671)ZUPS(IT) ENDDO C WRITE(6,690) DO IT=1,NT WRITE(6,695)TEMP(IT),ZUP1(IT),ZUP2(IT) ENDDO C C WRITE(6,650) C C FORMATS C 601 FORMAT(///1X,71('+')// + //10X,'THERMALLY-AVERAGED COLLISION STRENGTHS ZUPSILON'/ + 10X,'***********************************************'// + 10X,'TRANSITION I = ',I3,', J = ',I3// + 4X,'T/SZ**2',5X,'REAL AND IMAG. ZUPSILON',4X, + 'ZUPINC FOR TOP-UP IN E'/) 620 FORMAT(23X,I3,F12.6) 640 FORMAT(//5X,'***** KMAX = ',I2,', SHOULD BE AT LEAST 3', + ' FOR CALCULATION OF ZUPS *****'//) 650 FORMAT(/1X/10X,43('*')/ + 10X,'* THERMALLY-AVERAGED COLLISION STRENGTHS *'/ + 10X,'* WRITTEN TO FILE ZUPSILON',16X,'*'/ + 10X,43('*')//) 660 FORMAT(// '***** NUMBER OF TEMPERATURES REDUCED FROM NT = ',I3/ + 7X,'TO MAXIMUM MAXIMUM VALUE OF ',I3,' ALLOWED BY ', + 'PARAMETER NTMP') 670 FORMAT(1PE12.2,2(2X,2E12.3)) 671 FORMAT(1PE12.3,E12.3) 690 FORMAT(//17X,'ZUP1 FOR TOP-UP IN L',6X,'ZUP2 FOR TOP-UP IN L'/) 695 FORMAT(1PE12.2,2(2X,2E12.3)) 1001 FORMAT(2I5) C END C C********************************************************************* C SUBROUTINE TOPOK(ITOP,ITA,ITB,TOP,LRGLX,LRGLN) C C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (MXTST=(MZTAR*MZTAR+MZTAR)/2) C LOGICAL TOP C COMMON/CTARG/ENAT(MZTAR),NAST,NELC,NZED,ISAT(MZTAR),LAT(MZTAR) COMMON/CDAMP/SLINE(MXTST),NSLP(MZTAR) COMMON/NRBSLP/NSPNQ(MZSLP),LRGLQ(MZSLP),NPTYQ(MZSLP) X,ICONAT(MZTAR,MZSLP),LLCHQ(MZCHF,MZSLP),KJQ(MZCHF,MZSLP) X,NCASE(MZSLP,MZTAR) C C SEE IF LS- OR JK-COUPLING C IDIP=1 IF(NSPNQ(1).EQ.0)IDIP=2 C LRGLAX=0 LRGLN=10000 DO NA=1,NSLP(ITA) IIA=NCASE(NA,ITA) LRGL=LRGLQ(IIA) LRGLN=MIN(LRGLN,LRGL) LRGLAX=MAX(LRGLAX,LRGL) ENDDO LRGLBX=0 DO NB=1,NSLP(ITB) IIB=NCASE(NB,ITB) LRGL=LRGLQ(IIB) LRGLN=MIN(LRGLN,LRGL) LRGLBX=MAX(LRGLBX,LRGL) ENDDO C IF(LRGLAX.EQ.LRGLBX)THEN LRGLX=LRGLAX ELSE GOTO 1000 ENDIF IF(LRGLX.LE.2*IDIP)GOTO 1000 C IF(ITOP.EQ.0) RETURN C IIA=NCASE(1,ITA) IF(NSPNQ(IIA).GT.0)THEN IF(ISAT(ITA).EQ.1)THEN ISMIN=2 ELSE ISMIN=ISAT(ITA)-1 ENDIF ISMAX=ISAT(ITA)+1 ELSEIF(NSPNQ(IIA).LT.0)THEN ISMIN=NSPNQ(IIA) ISMAX=ISMIN ELSE ISMIN=0 ISMAX=0 ENDIF C DO 120 ILEV=1,2 IF(ILEV.EQ.1)THEN IT=ITA ELSE IT=ITB ENDIF DO 110 ISPIN=ISMIN,ISMAX,2 DO 100 LRGL=LRGLX,LRGLX-2*IDIP,-1*IDIP DO 40 N=NSLP(IT),1,-1 IIA=NCASE(N,IT) IF(NSPNQ(IIA).NE.ISPIN)GOTO 40 IF(LRGLQ(IIA).NE.LRGL)GOTO 40 IF(LAT(IT).EQ.0)GOTO 100 DO 30 M=N-1,1,-1 IIB=NCASE(M,IT) IF(NSPNQ(IIB).NE.ISPIN)GOTO 30 IF(LRGLQ(IIB).NE.LRGL)GOTO 30 GOTO 100 30 CONTINUE GOTO 1000 40 CONTINUE GOTO 1000 100 CONTINUE 110 CONTINUE 120 CONTINUE C TOP=.TRUE. RETURN C 1000 TOP=.FALSE. WRITE(6,600) RETURN C 600 FORMAT(//10X,' ***** INSUFFICIENT SL/JP VALUES FOR ', + 'TOP-UP *****'//) C END C C*********************************************************************** C REAL*8 FUNCTION WABS1(A,B,C,D,F) C C CALCULATES ABSOLUTE VALUE OF RACAH COEFFICIENT C W(A,B,C,D,1,F) C IMPLICIT REAL*8(W-Y) IMPLICIT INTEGER(A-F) C PARAMETER (TZERO=0.0) C IF(F.LT.IABS(A-C).OR.F.GT.IABS(A+C).OR. X F.LT.IABS(B-D).OR.F.GT.IABS(B+D).OR. X (A+B).EQ.0.OR.(C+D).EQ.0)THEN WABS1=TZERO RETURN ENDIF C IF(B-A)10,20,30 C 10 IF(D-C)11,12,13 20 IF(D-C)21,22,23 30 IF(D-C)31,32,33 C 11 X=(B+D+F+2)*(B+D+F+3)*(B+D-F+1)*(B+D-F+2) Y=(B+1)*(2*B+1)*(2*B+3)*(D+1)*(2*D+1)*(2*D+3) GOTO 40 12 X=(B+D+F+2)*(-B+D+F)*(B-D+F+1)*(B+D-F+1) Y=(B+1)*(2*B+1)*(2*B+3)*D*(D+1)*(2*D+1) GOTO 40 13 X=(-B+D+F-1)*(-B+D+F)*(B-D+F+1)*(B-D+F+2) Y=(B+1)*(2*B+1)*(2*B+3)*D*(2*D-1)*(2*D+1) GOTO 40 C 21 X=(B+D+F+2)*(-B+D+F+1)*(B-D+F)*(B+D-F+1) Y=B*(B+1)*(2*B+1)*(D+1)*(2*D+1)*(2*D+3) GOTO 40 22 X=B*(B+1)+D*(D+1)-F*(F+1) Y=B*(B+1)*(2*B+1)*D*(D+1)*(2*D+1) WABS1=.5*ABS(X)/SQRT(Y) RETURN 23 X=(B+D+F+1)*(-B+D+F)*(B-D+F+1)*(B+D-F) Y=B*(B+1)*(2*B+1)*D*(2*D-1)*(2*D+1) GOTO 40 C 31 X=(-B+D+F+1)*(-B+D+F+2)*(B-D+F-1)*(B-D+F) Y=B*(2*B-1)*(2*B+1)*(D+1)*(2*D+1)*(2*D+3) GOTO 40 32 X=(B+D+F+1)*(-B+D+F+1)*(B-D+F)*(B+D-F) Y=B*(2*B-1)*(2*B+1)*D*(D+1)*(2*D+1) GOTO 40 33 X=(B+D+F)*(B+D+F+1)*(B+D-F-1)*(B+D-F) Y=B*(2*B-1)*(2*B+1)*D*(2*D-1)*(2*D+1) GOTO 40 C C 40 WABS1=.5*SQRT(X/Y) C RETURN END C C*********************************************************************** C REAL*8 FUNCTION WABS2(A,B,C,D,F) C C CALCULATES |W(A,B,C,D,1,F)| WHERE W IS RACAH COEFFICIENT. C NRB: C*** AND *** A,B,C,D,E,F ARE INPUT TWICE THEIR ACTUAL VALUE. C IMPLICIT INTEGER(A-N) IMPLICIT REAL*8 (W) C PARAMETER (TZERO=0.0) C IF(F.LT.IABS(A-C).OR.F.GT.IABS(A+C).OR. X F.LT.IABS(B-D).OR.F.GT.IABS(B+D).OR. X (A+B).EQ.0.OR.(C+D).EQ.0)THEN WABS2=TZERO RETURN ENDIF C IF(B-A)10,20,30 C 10 IF(D-C)11,12,13 20 IF(D-C)21,22,23 30 IF(D-C)31,32,33 C 11 WSQ2=.25*DBLE((B+D+F+4)*(B+D+F+6)*(B+D-F+2)*(B+D-F+4))/ + DBLE((B+2)*(B+1)*(B+3)*(D+2)*(D+1)*(D+3))/4. GO TO 40 12 WSQ2=.25*DBLE((B+D+F+4)*(-B+D+F)*(B-D+F+2)*(B+D-F+2))/ + DBLE((B+2)*(B+1)*(B+3)*D*(D+2)*(D+1))/2. GO TO 40 13 WSQ2=.25*DBLE((-B+D+F-2)*(-B+D+F)*(B-D+F+2)*(B-D+F+4))/ + DBLE((B+2)*(B+1)*(B+3)*D*(D-1)*(D+1))/4. GO TO 40 C 21 WSQ2=.25*DBLE((B+D+F+4)*(-B+D+F+2)*(B-D+F)*(B+D-F+2))/ + DBLE(B*(B+2)*(B+1)*(D+2)*(D+1)*(D+3))/2. GO TO 40 22 G=-(B*(B+2)+D*(D+2)-F*(F+2)) WSQ2=.25*DBLE(G*G)/ + DBLE(B*(B+2)*(B+1)*D*(D+2)*(D+1)) GO TO 40 23 WSQ2=.25*DBLE((B+D+F+2)*(-B+D+F)*(B-D+F+2)*(B+D-F))/ + DBLE(B*(B+2)*(B+1)*D*(D-1)*(D+1))/2. GO TO 40 C 31 WSQ2=.25*DBLE((-B+D+F+2)*(-B+D+F+4)*(B-D+F-2)*(B-D+F))/ + DBLE(B*(B-1)*(B+1)*(D+2)*(D+1)*(D+3))/4 GO TO 40 32 WSQ2=.25*DBLE((B+D+F+2)*(-B+D+F+2)*(B-D+F)*(B+D-F))/ + DBLE(B*(B-1)*(B+1)*D*(D+2)*(D+1))/2. GO TO 40 33 WSQ2=.25*DBLE((B+D+F)*(B+D+F+2)*(B+D-F-2)*(B+D-F))/ + DBLE(B*(B-1)*(B+1)*D*(D-1)*(D+1))/4. GOTO 40 C C 40 if(wsq2.lt.0.)stop 'w**2.lt.0' WABS2=SQRT(WSQ2) C END