C PROGRAM SBE IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C PROGRAM FOR GENERALISED S.B.E. CROSS SECTIONS. C PROGRAM CALCULATES EQS (1) AND (2) OF THE PAPER BY FITZ, KOURI, C LIU, MCCOURT, EVANS AND HOFFMAN, J. PHYS. CHEM., 86 1087 (1982). C SEE ALSO, EQ (2.37) OF FITZ, ET AL., J. CHEM. PHYS. 74, 5022 C (1981). HOWEVER, IT THEN CONVERTS TO THE DEFINITIONS OF LIU, C ET AL., J. CHEM. PHYS. 71, 415 (1979) BY USING THE FUNCTION CFACT. C IT WOULD BE STRAIGHT-FORWARD TO MODIFY THIS LAST CONVERSION. C ADDITIONAL DISCUSSION OF THESE CROSS SECTIONS AND THEIR RELATION C TO EXPERIMENTAL OBSERVABLES CAN BE FOUND IN HUTSON AND MCCOURT, J. C CHEM. PHYS. 80, 1135 (1984.) C C PROGRAM BY SHELDON GREEN, NASA (WELL, PERHAPS ...) 1982 C MODIFIED BY JEREMY M. HUTSON TO ALLOW PLOTTING, CORRECT BUGS, C AND IMPROVE EFFICIENCY. C MODIFIED BY V. VESOVIC (1988) C PRODUCES PARTIAL WAVE OPACITIES TO CHANNEL 7 WITH 'INTERNAL' RESUM C MODIFIED BY C. BISSONNETTE (APRIL 1989) C FOR SUN FORTRAN 1.1 (RUNNING UNDER THE SUNOS 4.0) C OPTIONAL OUTPUT OF VALUES TO CHANNEL MXN. C MODIFIED MY MARK THACHUK IN DECEMBER 1989 C REMOVED ALL REMNANTS OF PLOTTING OPTION C ADD ICNTRL NAMELIST VARIABLE TO ALLOW FOR EXPLICIT CALC OF C JT1,JT2 AND JT2,JT1 CONTRIBUTIONS. C MODIFIED BY S GREEN (NOV 94) C MADE SAVE TAPE READS COMPATIBLE THROUGH VERSION 14 C &INPUT LFMT CONTROLS FORMATTED VS. UNFORMATTED (DEFAULT) C PARTIAL OPACITIES NOW CONTROLLED BY &INPUT IPARTU (DEFAULT 7) C MADE CLOCKS,DATE,TIME COMPATIBLE WITH MOLSCAT GCLOCK,GDATE,GTIME C C ----------------------------------------------------------------- C NAMELIST PARAMETERS ARE AS FOLLOWS: C NCALC - NUMBER OF ENERGIES FOR WHICH X-SECTIONS ARE TO BE C CALCULATED (DEFAULT = MXCALC, I.E., ALL ENERGIES) C ICALC - ARRAY OF POINTERS TO REQUIRED ENERGIES IN MOLSCAT C ENERGY ARRAY (DEFAULT=1,2,3,...) C ISU - CHANNEL FOR MOLSCAT S-MATRIX SAVE FILE (DEFAULT = 12) C MXSIG - IF .GT.0 (DEFAULT=0) MAX LEVELS FOR WHICH TO CALCULATE C CROSS SECTIONS C NKSET - NUMBER OF DIFFERENT CROSS-SECTIONS TO BE CALCULATED C DEFAULT IS NKSET=6 C KSET - ARRAY OF INTEGERS DESCRIBING TENSOR ORDERS INVOLVED C IN REQUIRED CROSS-SECTIONS. THE DEFAULT SET PROVIDED C IS AS FOLLOWS: C 1 - SOUND ABSORBTION C 2 - SIGMA(ETA,KAPPA=1) C 3 - SIGMA(ETA,KAPPA=2) C 4 - NMR SPIN-ROTATION C 5 - NMR DIPOLAR, DPR AND SBE SIGMA(T) C 6 - SBE SIGMA(ETA,T) C 7 - SBE SIGMA(T,ETA) C 8-11 THERMAL CONDUCTIVITY SBE C MXN - CHANNEL NUMBER FOR MORE PRECISELY REPORTED CROSS-SECTIONS C IPARTU- IF .GT. 0 WRITE PARTIAL OPACITIES TO THIS CHANNEL C NCCH - NUMBER OF CLOSED CHANNELS IN THE S-MATRIX C ONLY USE THIS WHEN THE NUMBER OF CLOSED CHANNELS IS C CONSTANT FOR ALL ENERGIES USED IN A GIVEN SBE+RESUM RUN. C THIS PARAMETER IS USEFUL FOR OMITTING THE CLOSED CHANNELS C FROM OUTPUT TO CHANNEL MXN. C IPRINT- A PRINT CONTROL (DEFAULT = 5) C C LFMT - TRUE/FALSE FOR FORMATTED (THROUGH VERSION 13)/ C /UNFORMATTED (VERSION 12 AND AFTER); DEFAULT IS .FALSE. C ICNTRL -DEFAULT = 2 ONLY DOES CALCS FOR JT2.GT.JT1 AND DOUBLES C OFF-DIAGONAL CONTRIBUTIONS C (WAS CALLED ICONTROL PREVIOUSLY) C ----------------------------------------------------------------- C PARAMETER (MXKSET=11) DIMENSION KSET(5,MXKSET),TTIME(4) C PARAMETER (MXCALC=30) LOGICAL LFMT COMMON /CNTRL/ NCALC,ICALC(MXCALC),ISU,MXN,NCCH,ICNTRL,IPARTU,LFMT C C FOR NAMLIS (NAMELIST SIMULATOR) -- C CHARACTER*6 INAMES C DIMENSION INAMES(12),LOCN(12),INDX(12) C CHARACTER*8 PDATE C ARRAY TO HOLD TIME AND DATE (SYSTEM DEPENDENT) C BELOW IS COMPATIBLE WITH MOLSCAT GDATE/GTIME ROUTINES C TO RUN ON THE VAX, USE: C CHARACTER CDATE*8,CTIME*8 CHARACTER CDATE*11,CTIME*9 C THE FOLLOWING IS PUT IN FOR FILENAME CONTROL - UNIV OF WATERLOO CHARACTER*1 FM(0:99) C DATA IPROGM/3/, PDATE/'(NOV 94)'/ C C DEFAULT CROSS SECTION LIST DATA KSET/ 0,0,0,0,0, 1,0,1,0,1, 2,0,2,0,2, 0,1,0,1,1, 0,2,0,2,2, & 2,0,0,2,2, 0,2,2,0,2, 1,0,1,2,1, & 1,2,1,2,1, 1,2,1,2,2, 1,2,1,2,3/ DATA TTIME/4*0.D0/ C NAMELIST /INPUT/ NCALC,ICALC,ISU,MXSIG,NKSET,KSET,IPRINT, 1 MXN,NCCH,ICNTRL,LFMT,IPARTU C C FOLLOWING ARE USED IN NAMLIS (NAMELIST SIMULATOR) C DATA INAMES/'NCALC','ICALC','ISU','MXSIG','NKSET','KSET', C 1 'IPRINT','MXN','NCCH','ICNTRL','LFMT','IPARTU'/ C DATA INDX/11*0/ C DATA (FM(I), I=0,99)/100*'F'/ C C MASK UNDERFLOW (MASK IS A MOLSCAT INTEFACE TO MACHINE DEP CODE) C CALL MASK C CALL GCLOCK(XKTIME) CALL GDATE(CDATE) CALL GTIME(CTIME) WRITE(6,100) IPROGM,PDATE,CDATE,CTIME 100 FORMAT('0'/'0 ',8('----- SBE -----')/' |',120X,'|'/' |',7X, 1 'TRANSPORT AND RELAXATION CROSS SECTION PROGRAM OF S. GREEN ', 2 '(AUG 82) - J.M. HUTSON (CCP6) VERSION',I2,1X,A8,6X,'|'/ 3 ' |',120X,'|'/' |',45X,'RUN ON ',A11, 4 ' AT ',A9,44X,'|'/' |',120X,'|'/' ',8('----- SBE -----') ) C C SET DEFAULTS C NCALC=MXCALC NKSET=6 DO 1000 I=1,MXCALC 1000 ICALC(I)=I ISU=12 MXSIG=0 IPRINT=5 MXN=0 IPARTU=7 NCCH=0 ICNTRL = 2 C C -------------------- C LOCN(1)=LOC(NCALC) C LOCN(2)=LOC(ICALC) C LOCN(3)=LOC(ISU) C LOCN(4)=LOC(MXSIG) C LOCN(5)=LOC(NKSET) C LOCN(6)=LOC(KSET) C LOCN(7)=LOC(IPRINT) C LOCN(8)=LOC(MXN) C LOCN(9)=LOC(NCCH) C LOCN(10)=LOC(ICNTRL) C LOCN(11)=LOC(LFMT) C LOCN(12)=LOC(IPARTU) C CALL NAMLIS('&INPUT',INAMES,LOCN,INDX,9,IEOF) C OPEN(5,STATUS='OLD',SHARED,READONLY) C WATSCI NEEDS TO HAVE THE NML IN NAMELIST READ C READ(5,NML=INPUT,END=2000) C ----------------------- C C OPEN(5,STATUS='OLD') --- this caused problems on rs/6000 READ(5,INPUT,END=2000) 2000 CONTINUE C CHECK FOR LEGAL ISU IF (ISU.LE.0 .OR. ISU.GE.100 .OR. ISU.EQ.5 .OR. 1 ISU.EQ.6 .OR.ISU.EQ.7) THEN WRITE(6,*) ' *** ILLEGAL SAVE TAPE UNIT ISU',ISU STOP ENDIF C C BELOW FOR UNIV OF WATERLOO MACHINE C ON THE ASSUMPTION THAT 'U' IS UNFORMATTED AND 'F' FORMATTED ... FM(ISU)='U' IF (LFMT) FM(ISU)='F' CALL IOINIT(.TRUE.,.FALSE.,.FALSE.,'FORT',.FALSE.,FM) C END OF WATERLOO EXTENSION C C CHECK VALIDITY OF SOME OTHER INPUTS IF (NCALC.GT.MXCALC) THEN WRITE(6,*) ' *** NCALC LARGER THAN MXCALC BEING REDUCED', 1 NCALC,MXCALC NCALC=MXCALC ENDIF IF (MXN.LT.0 .OR. MXN.GT.99 .OR. MXN.EQ.5 .OR. MXN.EQ.6 1 .OR. MXN.EQ.7 .OR. MXN.EQ.ISU) THEN WRITE(6,*) ' *** ILLEGAL MXN',MXN WRITE(6,*) ' WILL BE RESET TO ZERO' MXN=0 ENDIF IF (NKSET.GT.MXKSET) THEN WRITE(6,*) ' *** NKSET.GT.MXKSET BEING REDUCED',NKSET,MXKSET NKSET=MXKSET ENDIF IF (IPARTU.GT.0) WRITE(6,*) 1 ' PARTIAL OPACITIES WILL BE WRITTEN TO UNIT',IPARTU C C FINALLY, CALL SBEDRV TO DO PROCESSING ... CALL SBEDRV(MXSIG,NKSET,KSET,IPRINT,TTIME) C CALL GCLOCK(XMTIME) TTIME(4)=(XMTIME-XKTIME) TTIME(1)=TTIME(4) - TTIME(2) TTIME(2)=TTIME(2) - TTIME(3) WRITE(6,601) (TTIME(I),I=1,4) 601 FORMAT('0',13('----')/'0 TIMING SUMMARY (SECS):'/'0 MAIN PROGRAM', 1 T30,F12.3/' SBEADD',T30,F12.3/' ANGULAR MOMENTUM',T30,F12.3/ 2 ' TOTAL',T30,F12.3) STOP END FUNCTION CFACT(KSET,JA,JB) C C FUNCTION TO CONVERT CROSS-SECTION FROM FITZ, ET AL., C J. PHYS. CHEM., 86 1087 (1982) C DEFINITION TO LIU AND MCCOURT DEFINITION. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION KSET(5) FAC(I)=SQRT(DBLE(I+I+1)) X=-1.0D0 IF(KSET(1).EQ.KSET(3)) X=-X IF(KSET(2).EQ.KSET(4)) X=-X Y=FAC(KSET(5)) DO 10 K=1,4 10 IF(KSET(K).EQ.0) X=X*Y IF(KSET(2).EQ.0) X=X*FAC(JA) IF(KSET(4).EQ.0) X=X*FAC(JB) CFACT=X RETURN END SUBROUTINE HUBY(N,L,SR,SI) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C THIS ROUTINE APPLIES HUBY'S PHASE CORRECTION TO THE C ARTHURS-DALGARNO S-MATRICES. C DIMENSION L(N),SR(N,N),SI(N,N) C DO 1000 I1=1,N DO 1000 I =1,N IP=L(I1)-L(I) 2000 IF(IP.GE.0) GOTO 1001 IP=IP+4000 GOTO 2000 1001 MOD=IP-4*(IP/4) IF(MOD.EQ.0) GOTO 1000 IF(MOD-2) 2001,2002,2003 C BELOW FOR S=I * S 2001 T=SR(I1,I) SR(I1,I)=-SI(I1,I) SI(I1,I)=T GOTO 1000 C BELOW FOR S = - S 2002 SR(I1,I)=-SR(I1,I) SI(I1,I)=-SI(I1,I) GOTO 1000 C BELOW IS FOR S = - I * S 2003 T=-SR(I1,I) SR(I1,I)=SI(I1,I) SI(I1,I)=T 1000 CONTINUE C RETURN END SUBROUTINE OUTP(ENERGY,NLEV,JLEV,MXLV,NKSET,KSET,IPARTU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/CVEL/ X(11,20,20,150),JTMAX COMMON/CVEL1 / Y(11,20,20) DIMENSION JLEV(NLEV),KSET(5,NKSET),SIG(5000) C IF (IPARTU.NE.0) 1 WRITE(IPARTU,607) ENERGY 607 FORMAT('1',19('====')/'0 (2J+1) WEIGHTED PARTIAL OPACITIES AT ', + 'ENERGY',F12.6) NJT=JTMAX-1 DO 500 IJT=1,NJT JT=IJT-1 IF (IPARTU.NE.0) 1 WRITE(IPARTU,602)JT 602 FORMAT('0',19('====')/'0 OPACITIES FOR JTI =',I3/) IC=0 DO 10 IB=1,MXLV DO 10 IA=1,MXLV DO 10 I=1,NKSET IC=IC+1 10 SIG(IC)=X(I,IA,IB,IJT) C IF (IPARTU.NE.O) 1 CALL OUTSBE(NLEV,JLEV,MXLV,SIG,NKSET,KSET,JT,IPARTU) 500 CONTINUE C IC=0 DO 11 IB=1,MXLV DO 11 IA=1,MXLV DO 11 I=1,NKSET IC=IC+1 11 SIG(IC)=Y(I,IA,IB) IF (IPARTU.NE.0) 1 WRITE(IPARTU,603) 603 FORMAT('0',19('****')/'0 SUM OF THESE VALUES'/) IF (IPARTU.NE.O) 1 CALL OUTSBE(NLEV,JLEV,MXLV,SIG,NKSET,KSET,JT,IPARTU) RETURN END SUBROUTINE OUTSBE(NLEV,JLEV,MXLV,SIG,NKSET,KSET,IQUIET,ISU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION JLEV(NLEV),SIG(NKSET,MXLV,MXLV),KSET(5,NKSET) C SIMPLE OUTPUT ROUTINE FOR S.B.E. CROSS SECTIONS. CHARACTER*8 HOL DATA HOL/' K-SET ='/ C C MXOUT IS MAXIMUM NO. OF K-COLUMNS OUTPUT AT A TIME . . . DATA MXOUT/ 9/ C C STATEMENT FUNCTION FOR CONVERSION FN(I)=DBLE(I+I+1) C IF(IQUIET.EQ.0) WRITE(ISU,651) (K,(KSET(I,K),I=1,5),K=1,NKSET) 651 FORMAT('0 CROSS SECTIONS FOR FOLLOWING TENSOR ORDERS'/ 1 '0 SET KL1 KJ1 KL KJ BIG K'/ (6I5) ) C KLOW=1 KTOP=MIN0(MXOUT,NKSET) C 2000 WRITE(ISU,600) (HOL,K,K=KLOW,KTOP) 600 FORMAT('0'/'0 JI JF',10(2X,A8,I2) ) DO 1000 LF=1,MXLV JF=JLEV(LF) DO 1000 LI=1,MXLV JI=JLEV(LI) CONV=SQRT(FN(JF)/FN(JI)) 1000 WRITE(ISU,601) JI,JF,(SIG(K,LF,LI),K=KLOW,KTOP), CONV 601 FORMAT(' ',2I4,10F12.5) KLOW=KLOW+MXOUT KTOP=MIN0(KTOP+MXOUT,NKSET) C IF MORE K-COLUMNS ARE LEFT, GO BACK AND OUTPUT THEM. IF (KTOP.GE.KLOW) GO TO 2000 C RETURN END SUBROUTINE RESUB(K,LEVA,LEVB,JT1,JT2,LI,LI1,LF,LF1,CONTRB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/CVEL/ X(11,20,20,150),JTMAX COMMON/CVEL1 / Y(11,20,20) C INTRODUCE LOGICAL FUNCTION OK TO CHECK INDICES VS DIMENSIONS IN X,Y LOGICAL OK OK(LA,LB,KK,I1,I2)=LA.LE.20.AND.LB.LE.20.AND.KK.LE.11 1 .AND.I1.LE.150.AND.I2.LE.150 C JTMAX=MAX0(JTMAX,JT1) JTMAX=MAX0(JTMAX,JT2) CONTRB=0.5D0*CONTRB IJT1=1+JT1 IJT2=1+JT2 IF (OK(LEVA,LEVB,K,IJT1,IJT2)) GO TO 150 WRITE(6,601) LEVA,LEVB,K,JT1,JT2 601 FORMAT(/' *** RESUB. EXCEEDED DIMENSIONS OF X, Y'/ 1 ' LEVA,LEVB,K,JT1,JT2 =',5I6) STOP C 150 X(K,LEVA,LEVB,IJT1)=X(K,LEVA,LEVB,IJT1)+CONTRB X(K,LEVA,LEVB,IJT2)=X(K,LEVA,LEVB,IJT2)+CONTRB Y(K,LEVA,LEVB)=Y(K,LEVA,LEVB)+2.D0*CONTRB RETURN END SUBROUTINE RESUM(K,LEVA,LEVB,JT1,JT2,LI,LI1,LF,LF1 1 ,CONTRB,NKSET,MXLV,IFLAG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/CVEL/ X(11,20,20,150),JTMAX COMMON/CVEL1 / Y(11,20,20) C IF(IFLAG.EQ.1) GOTO 20 C JTMAX=0 C IF (MXLV.GT.20.OR.NKSET.GT.11) THEN WRITE(6,*) ' *** RESUM. INITIALIZATION OF X,Y EXCEEDS DIMS' WRITE(6,*) ' NKSET,MXLV =',NKSET,MXLV STOP ENDIF DO 10 I=1,NKSET DO 10 IA=1,MXLV DO 10 IB=1,MXLV Y(I,IA,IB)=0.D0 DO 10 IC=1,150 10 X(I,IA,IB,IC)=0.D0 C 20 CALL RESUB(K,LEVA,LEVB,JT1,JT2,LI,LI1,LF,LF1,CONTRB) IFLAG=1 RETURN END SUBROUTINE SBEADD(JT1,NB1,J1,L1,WV1,TR1,TI1,JT2,NB2,J2,L2,WV2, 2 TR2,TI2,ICODE,JLEV,NLEV,SIG,MXLV,KSET,NKSET, 3 NMATCH,IMATCH,TIME,IPRINT,IVV,KDIFMX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE IPR,PI,IFLAG DOUBLE PRECISION NMATCH, IMATCH C C THIS ROUTINE COMPUTES CONTRIBUTION TO SIG() FROM C TWO T-MATRICES (SUFFIX -1 AND -2) WHICH ARE INPUT IN SBEADD. C C * * * * * * * * * * * * * * * * * * * * * * * * * C C FROM WHAT I CAN SEE, THE SYMMETRY OF THE S MATRICES IS NOT USED HERE. C SO, YOU CAN USE NONSYMMETRIC S MATRICES (LIKE CS ONES) WITH THIS C CODE (INCLUDING THE DOUBLING) AND STILL GET THE RIGHT ANSWER. (MET) C C * * * * * * * * * * * * * * * * * * * * * * * * * C C DIMENSIONS FOR INPUT ARGUMENTS. . . DIMENSION JLEV(NLEV),SIG(NKSET,MXLV,MXLV),KSET(5,NKSET) C DIMENSION J1(NB1),L1(NB1),TR1(NB1,NB1),TI1(NB1,NB1),WV1(NB1), 1 J2(NB2),L2(NB2),TR2(NB2,NB2),TI2(NB2,NB2),WV2(NB2), 2 NMATCH(NB1),IMATCH(NB1,NB2) LOGICAL NEXTL,NEXTL1 C C ICODE = 0 FOR T1 = T2 C ICODE = 1 FOR PLOTTING OPTION: DO NOT DOUBLE CONTRIBUTION C = 2 FOR GENERAL CASE OF DIFFERENT S-MATRICES C CONTRIBUTION FROM THIS PAIR IS DOUBLED. C C TOLERANCE FOR SETTING CONTRIBUTION TO ZERO . . . . DATA EPS/1.D-15/ C C IF(IVV.EQ.1) GOTO 99 IPR=IPRINT PI=ACOS(-1.D0) IFLAG=0 RETURN C CHECK TO SEE IF NO CONTRIBUTION DUE TO TRIANGLE ON JT1,JT2,K 99 JDIFF=IABS(JT1-JT2) IF(JDIFF.GT.KDIFMX) RETURN C XJFCT=PI*DBLE((JT1+JT1+1)*(JT2+JT2+1)) C *** TO COMPENSATE FOR SUMMING OVER JT1 .GE. JT2 ONLY . . . IF(ICODE.EQ.2) XJFCT=2.D0*XJFCT C *** PROVIDED KJ AND KJ1 ARE OF THE SAME PARITY, XFCT=2.*XFCT WILL WORK C *** OTHERWISE CROSS SECTION MAY NOT BE REAL AND THIS CODE IS INVALID C IF(IPR.GE.20) WRITE(6,601) JT1,JT2 601 FORMAT('0 CONTRIBUTIONS FROM S-MATRICES WITH JTOT =',I3,' AND',I3/ ' '0 NK JA JB LILI1 LFLF1 S*S 3J 3J 9J', 2 ' 9J SIGMA'/) DO 3000 K=1,NKSET KK=KSET(5,K) IF(JDIFF.GT.KK) GOTO 3000 KL1=KSET(1,K) KJ1=KSET(2,K) KL=KSET(3,K) KJ=KSET(4,K) C C FIND OUT WHICH PARITY CASE THIS CROSS SECTION IS C QPAR = PARITY(KJ1+KJ) C C LOOP OVER ROWS/COLS OF 1ST S-MATRIX. C TO FIND THOSE ELEMENTS THAT CONTRIBUTE TO THE CROSS SECTION C EACH ONE DEFINES A JA(=J-PRIME) AND JB(=J). DO 1003 IR1=1,NB1 IM=0 LEVA=J1(IR1) C SKIP IF LEVEL IS HIGHER THAN WE ARE COLLECTING . . . IF(LEVA.GT.MXLV) GOTO 1003 JA=JLEV(LEVA) LF=L1(IR1) C ABOVE ASSUMES THAT 'J' IS 1ST INDEX IN JLEV(NLEV,NQN) DO 1002 IR2=1,NB2 IF(J2(IR2).NE.LEVA) GOTO 1002 LI=L2(IR2) LD=IABS(LF-LI) C C SATISFY TRIANGLE RELATIONS BETWEEN L'S AND K'S C IF(KL1.LT.LD .AND. KL.LT.LD) GOTO 1002 LS1=LD+KL1 LS2=LD+KL IF(2*(LS1/2).NE.LS1 .AND. 2*(LS2/2).NE.LS2) GOTO 1002 IM=IM+1 IMATCH(IR1,IM)=FLOAT(IR2) 1002 CONTINUE 1003 NMATCH(IR1)=FLOAT(IM) C DO 2005 IR1=1,NB1 IF(INT(NMATCH(IR1)).EQ.0) GOTO 2005 LEVA=J1(IR1) JA=JLEV(LEVA) C C XJ? ARE USED TO SIMPLIFY THE 9-J SYMBOLS WHEN ONE OF THE K'S IS 0 C XJA=1.D0/SQRT(DBLE((JA+JA+1)*(KK+KK+1))) LF1=L1(IR1) XLF1=1.D0/SQRT(DBLE((LF1+LF1+1)*(KK+KK+1))) NEXTL1=.FALSE. DO 2004 IC1=1,NB1 IF(INT(NMATCH(IC1)).EQ.0) GOTO 2004 LEVB=J1(IC1) JB=JLEV(LEVB) XJB=1.D0/SQRT(DBLE((JB+JB+1)*(KK+KK+1))) LF=L1(IC1) XLF=1.D0/SQRT(DBLE((LF+LF+1)*(KK+KK+1))) NEXTL=.FALSE. DENOM=SQRT(DBLE((LF+LF+1)*(LF1+LF1+1)))/WV1(IC1)**2 C CONVERT TO LIU-MCCOURT DEFINITION OF CROSS-SECTIONS TRANSL=XJFCT*DENOM*CFACT(KSET(1,K),JA,JB) CALL GCLOCK(XITIME) C C LOOP OVER ROWS/COLS OF 2ND S-MATRIX TO FIND MATCH WITH JA,JB C DO 2003 IM1=1,INT(NMATCH(IR1)) IR2=INT(IMATCH(IR1,IM1)) LI1=L2(IR2) X3J2=THREEJ(LF1,KL1,LI1) IF(ABS(X3J2).LT.EPS) GOTO 2003 C C USE SIMPLIFICATION OF 9-J SYMBOL WHEN ONE OF THE K'S IS 0 C IF(KJ1.EQ.0) THEN X9J2=PARITY(JA+KK-LI1-JT1)*SIXJN(LI1,JT2,LF1,JT1,JA,KK,NEXTL1,2) 1 *XJA NEXTL1=.TRUE. ELSE IF(KL1.EQ.0) THEN X9J2=PARITY(LF1+KK-JT2-JA)*SIXJ(JT1,JA,JT2,JA,LF1,KK)*XLF1 ELSE X9J2=XNINEJ(JA,LF1,JT1,JA,LI1,JT2,KJ1,KL1,KK) ENDIF IF(ABS(X9J2).LT.EPS) GOTO 2003 X3J2=X3J2*SQRT(DBLE(LI1+LI1+1)) C DO 2002 IM2=1,INT(NMATCH(IC1)) IC2=INT(IMATCH(IC1,IM2)) LI=L2(IC2) X3J1=THREEJ(LF,KL,LI) IF(ABS(X3J1).LT.EPS) GOTO 2002 IF(KJ.EQ.0) THEN X9J1=PARITY(JB+KK-LI-JT1)*SIXJN(LI,JT2,LF,JT1,JB,KK,NEXTL,1) 1 *XJB NEXTL=.TRUE. ELSE IF(KL.EQ.0) THEN X9J1=PARITY(LF+KK-JT2-JB)*SIXJ(JT1,JB,JT2,JB,LF,KK)*XLF ELSE X9J1=XNINEJ(JB,LF,JT1,JB,LI,JT2,KJ,KL,KK) ENDIF IF(ABS(X9J1).LT.EPS) GOTO 2002 X3J1=X3J1*SQRT(DBLE(LI+LI+1)) AFCT=X3J1*X3J2*X9J1*X9J2 C C IF QPAR IS EVEN, CALCULATE REAL PART OF SS*, IF QPAR IS ODD, CALCULATE C IMAGINARY PART OF SS* ONLY C IF( QPAR .GT. 0.0D0 ) THEN C SS IS REAL PART OF S(JF)*COMPLEX-CONJUGATE(S(JI)) SS=TR1(IR1,IC1)*TR2(IR2,IC2)+TI1(IR1,IC1)*TI2(IR2,IC2) IF(JA.EQ.JB .AND. LF.EQ.LF1 .AND. LI.EQ.LI1) SS=SS-1.D0 ELSE C SS IS IMAGINARY PART OF S(JF)*COMPLEX-CONJUGATE(S(JI)) SS = TI1(IR1,IC1)*TR2(IR2,IC2) - TR1(IR1,IC1)*TI2(IR2,IC2) ENDIF IF(ABS(SS).LT.EPS) GOTO 2002 CONTRB=PARITY(LI-LI1)*SS*AFCT*TRANSL SIG(K,LEVA,LEVB)=SIG(K,LEVA,LEVB)+CONTRB CALL RESUM(K,LEVA,LEVB,JT1,JT2,LI,LI1,LF,LF1,CONTRB 1 ,NKSET,MXLV,IFLAG) IF(IPR.LT.20) GOTO 2002 WRITE(6,602) K,JA,JB,LI,LI1,LF,LF1,SS,X3J1,X3J2,X9J1,X9J2,CONTRB IF(IPR.LT.30) GOTO 2002 WRITE(6,602) K,JA,JB,LI,LI1,LF,LF1,AFCT,TRANSL 602 FORMAT(7I3,5F9.5,F11.7) 703 FORMAT(10I3,F20.16) C 2002 CONTINUE 2003 CONTINUE C END OF LOOPS OVER 2ND S-MATRIX C CALL GCLOCK(XJTIME) TIME=TIME+(XJTIME-XITIME) C 2004 CONTINUE 2005 CONTINUE C END OF LOOPS OVER 1ST S-MATRIX C 3000 CONTINUE C END OF LOOP OVER K VALUES C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C RETURN END SUBROUTINE SBEDRV(MAXSIG,NKSET,KSET,IPRINT,TTIME) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C PROGRAM FOR SENFTLEBEN-BEENAKKER EFFECT CROSS SECTIONS, C USING THE FORMALISM OF D.E. FITZ (J.CHEM.PHYS. 74, 5022 (1981)) C C MODIFIED BY MARC TER HORST MARCH 9,1993 C CHANGED READ STATEMENTS TO READ V11'S UNFORMATTED S-MATRIX FILE C C MODIFIED BY MARK THACHUK ON APRIL 3, 1990 C TRIED TO SAVE SOME MEMORY. IN THE OLD VERSION OF SBE, ALL OF THE C S MATRICES WERE SAVED. LATER, FOR A GIVEN J, IT WOULD CHECK EACH C PREVIOUSLY STORED J ARRAY AND CHECK TO SEE IF IT COUPLED WITH THE C CALC BEING MADE WITH THE CURRENT J. FOR LARGE S MATRIX FILES C THE ARRAY SPACE NEEDED TO STORE ALL OF THESE ELEMENTS WAS HUGE. C HOWEVER, THE J'S ONLY COUPLE ACCORDING TO THE VALUE OF K C AND K'S ARE USUALLY 1, 2, OR 3. SO, YOU ONLY NEED TO STORE C THE LAST FEW J VALUES, NOT EVERY ONE. C NOTE: THIS CODE NOW ASSUMES J'S ARE IN INCREASING ORDER IN THE C S MATRIX FILE. THIS WAS NOT NECESSARY FOR THE OLD CODE. C THIS RESTRICTION IS WHAT ALLOWS THROWING AWAY THE OLDER J VALUES. C SO, THIS CODE NOW KEEPS ONLY 8 J VALUES AT ANY ONE TIME. WHEN C THE NUMBER OF J VALUES PASSES THIS, THE ARRAYS ARE MOVED C DOWN SO THAT THE TOTAL NUMBER STAYS AT 8. C NOTE: BECAUSE EACH VALUE OF J HAS 2 PARITIES ASSOCIATED WITH IT, C SAVING THE LAST 8 MATRICES ACTUALLY SAVES THE LAST 4 VALUES C OF J SO THAT COUPLING CAN OCCUR UP AND DOWN BY 3. I.E., CANNOT C HAVE K > 3 WITHOUT MAKING MODIFICATIONS TO THIS CODE. C C MODIFIED BY S GREEN NOV 94 C PROGRAM OBTAINS S MATRICES FROM (ISAVEU) TAPE OUTPUT BY VARIOUS C VERSIONS OF MOLSCAT. C SET NAMELIST &INPUT VARIABLE LFMT=.TRUE. FOR VERSION 11 AND BEFORE C DIMENSION KSET(5,NKSET),TTIME(4) C PARAMETER (MXCALC=30) LOGICAL LFMT COMMON /CNTRL/ NCALC,ICALC(MXCALC),ISU,MXN,NCCH,ICNTRL,IPARTU,LFMT C COMMON /CVEL1/ Y(11,20,20) C C STORAGE FOR SAVE TAPE INPUT. UP TO 100 (MXCH) CHANNELS / S-MATRIX C MXNRG, MXEL COMPATIBLE WITH MOLSCAT VERSION 14 BUILT-IN LIMITS PARAMETER (MXNRG=100, MXEL=1000) DIMENSION ENERGY(MXNRG),ELVL(MXEL) C PARAMETER (MXJT=8, MXBS=2000, MXCH=300, MXTM=400000, MXSG=5000, 1 MXJLEV=300) DIMENSION TR(MXTM),TI(MXTM),WV(MXBS),J(MXBS),L(MXBS) DIMENSION JLEV(MXJLEV),LEV(MXCH),LORB(MXCH),WVEC(MXCH) DIMENSION SIG(MXSG) DIMENSION JT(MXJT),ISTBS(MXJT),ISTTM(MXJT),NBS(MXJT) C CHARACTER*4 LABEL(20) C DATA ELVL/MXEL*0.D0/ C C FORMATS FOR SAVE TAPE . . . 100 FORMAT(20A4/3I4,F8.4,I4) 101 FORMAT(20I4) 102 FORMAT(I4/(5E16.8)) 103 FORMAT(2I4,E16.8,I4,E16.8,I4) 104 FORMAT(I4/(2I4,E16.8)) 114 FORMAT(I4) 105 FORMAT(5E16.8) C C SET KDIFMX AS MAX KSET(5,K) AND CHECK THAT JTMX IS ADEQUATE C KDIFMX=0 DO 300 K=1,NKSET IF(PARITY(KSET(2,K)+KSET(4,K)).GT.0.D0) GOTO 300 WRITE(6,690) K,KSET(2,K),KSET(4,K) 690 FORMAT('0 * * *'/'0 * * * WARNING, FOR',I4,'-TH KSET, KJ AND ' & ,'KJ1 NOT OF SAME PARITY --',2I4,' -- SIGMA MAY NOT BE REAL.'/ & '0 * * *') 300 KDIFMX=MAX0(KDIFMX,KSET(5,K)) IF (2*KDIFMX.GT.MXJT) THEN WRITE(6,*) ' *** SBEDRV. 2*KDIFMX.GT.MXJT',2*KDIFMX,MXJT WRITE(6,*) ' INADEQUATE S-MATRICES WILL BE STORED' STOP ENDIF C C NPASS IS NUMBER OF PASS THROUGH SAVE TAPE. ONE PASS / NCALC. NPASS=1 C C OPEN STATEMENTS FOR THE SAVE TAPE ON UNIT ISU; NO ERROR CHECKING IF (LFMT) THEN WRITE(6,*) ' *** LFMT SPECIFIES A FORMATTED TAPE' OPEN(ISU,STATUS='OLD',FORM='FORMATTED') ELSE WRITE(6,*) ' *** LFMT SPECIFIES AN UN FORMATTED TAPE' OPEN(ISU,STATUS='OLD',FORM='UNFORMATTED') ENDIF WRITE(6,650) ISU 650 FORMAT('0 COMPUTE CROSS SECTIONS FROM S-MATRICES' 1, ' STORED ON TAPE UNIT(',I3,' ).' ) WRITE(6,651) (K,(KSET(I,K),I=1,5),K=1,NKSET) 651 FORMAT('0 PROCESSING WILL BE DONE FOR FOLLOWING TENSOR ORDERS'/ 1 '0 SET KL1 KJ1 KL KJ BIG K'/ (6I5) ) C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C 1000 CALL GCLOCK(XITIME) XITYME=XITIME NMAX=0 C START READING TAPE REWIND ISU IF (LFMT) THEN READ(ISU,100) LABEL,ITYPE,NLEV,NQN,URED,IPROGM ELSE READ(ISU) LABEL,ITYPE,NLEV,NQN,URED,IPROGM ENDIF IF (LFMT.AND.IPROGM.GE.11) THEN WRITE(6,*) ' *** APPARENT ERROR IN SAVE TAPE FORMAT' STOP ENDIF IF(ITYPE.EQ.1 .OR. ITYPE.EQ.2 .OR. ITYPE.EQ.7) 1 GOTO 1991 WRITE(6,679) ITYPE 679 FORMAT(/' *** SORRY! NO IMPLEMENTATION FOR ITYPE =',I5) STOP 1991 IF(NPASS.EQ.1) THEN WRITE(6,600) LABEL,ITYPE,NLEV,NQN,URED,IPROGM IF (IPARTU.NE.0) 1 WRITE(IPARTU,600) LABEL,ITYPE,NLEV,NQN,URED,IPROGM 600 FORMAT('0 LABEL FROM SAVE TAPE =',20A4/'0 ITYPE =',I5,' NLEV =', 1 I5, ' NQN =',I5,' URED =',F8.4,' MOLSCAT VERSION',I4) ENDIF NSQ=NLEV*NQN IF (NSQ.GT.MXJLEV) THEN WRITE(6,*) ' *** NLEV*NQN EXCEES MXJLEV',NLEV,NQN,MXJLEV STOP ENDIF IF (LFMT) THEN READ(ISU,101) (JLEV(I),I=1,NSQ) ELSE READ(ISU) (JLEV(I),I=1,NSQ) ENDIF IF(NPASS.EQ.1) THEN WRITE(6,601)(JLEV(I),I=1,NSQ) IF (IPARTU.NE.0) 1 WRITE(IPARTU,601) (JLEV(I),I=1,NSQ) 601 FORMAT('0 JLEV ='/(' ',30I4)) ENDIF IF (LFMT) THEN IF(IPROGM.GE.3) READ(ISU,102) NLVL,(ELVL(I),I=1,NLVL) ELSE READ(ISU) NLVL,(ELVL(I),I=1,NLVL) ENDIF IF(NPASS.EQ.1.AND.IPROGM.GE.3) THEN WRITE(6,620) NLVL,(ELVL(I),I=1,NLVL) IF (IPARTU.NE.0) 1 WRITE(IPARTU,620) NLVL,(ELVL(I),I=1,NLVL) 620 FORMAT('0 THE',I3,' BASIS FUNCTION LEVEL ENERGIES ARE'/ & (1X,8F15.4)) ENDIF IF (LFMT) THEN READ(ISU,102) NNRG,(ENERGY(I),I=1,NNRG) ELSE READ(ISU) NNRG,(ENERGY(I),I=1,NNRG) ENDIF IF(NPASS.EQ.1) WRITE(6,602) NNRG,(ENERGY(I),I=1,NNRG) 602 FORMAT('0 TAPE CONTAINS',I4,' ENERGIES ='/(' ',8F15.4)) IF(NPASS.EQ.1) & WRITE(6,640) NCALC,(ICALC(I),ENERGY(ICALC(I)),I=1,NCALC) 640 FORMAT('0 CROSS SECTIONS WILL BE CALCULATED FOR',I4,' ENERGIES'/ & (3X,I3,F15.4)) C IF (ICALC(NPASS).LE.0.OR.ICALC(NPASS).GT.NNRG) THEN WRITE(6,641) NPASS,ICALC(NPASS) 641 FORMAT(' *** ILLEGAL ICALC(',I3,') =',I4,' WILL BE SKIPPED.') GO TO 1000 ENDIF C C SET MXLV - INPUT MXSIG CAN BE USED FOR LESS THAN NLEV LEVELS. MXLV=NLEV IF(MAXSIG.LE.0) GOTO 1002 IF(NPASS.EQ.1) WRITE(6,603) MAXSIG 603 FORMAT('0 * * * NOTE. MXLV SET FROM &INPUT MXSIG =',I4) MXLV=MIN0(MXLV,MAXSIG) 1002 ITOP=NKSET*MXLV*MXLV IF(ITOP.LE.MXSG) GOTO 1001 WRITE(6,699) MXSG,ITOP 699 FORMAT('0 * * * ERROR. TOO MANY LEVELS FOR STORAGE. MXSG=',2I4) STOP C ZERO STORAGE IN SIG AND INITIALIZE SBEADD (LAST PARM = 0) 1001 DO 1111 I=1,ITOP 1111 SIG(I)=0.D0 C SET NB1=1 SO THAT ARRAYS DO NOT HAVE ZERO DIMENSION UPON FIRST C CALL TO SBEADD. NB1=1 CALL SBEADD(JT1,NB1,J,L,WV,TR,TI,JT1,NB1,J,L,WV,TR,TI,IC,JLEV, & NLEV,SIG,MXLV,KSET,NKSET,TR,TI,TTIME,IPRINT,0,KDIFMX) C C BEGIN PROCESSING S-MATRIX INPUT. NSMAT COUNTS THESE. IF(ICALC(NPASS).LE.NNRG) GOTO 3000 WRITE(6,689) ICALC(NPASS),NNRG 689 FORMAT('0 REQUESTED ENERGY .GT. NNRG',2I8) GOTO 1900 C 3000 NSMAT=0 ICOUNT = 1 IXBS=0 IXTM=0 WRITE(6,611) 611 FORMAT('1 PROGRAM TO COMPUTE S.B.E. CROSS SECTIONS. ') C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C GET 'HEADER' RECORD FOR THIS SET OF S-MATRICES 2000 IF (LFMT) THEN READ(ISU,103,END=9000) JTOT,INRG,ECHK,ITMPEX,TMPWT,MTMP ELSE IF (IPROGM.LT.14) THEN READ(ISU,END=9000) JTOT,INRG,ECHK,ITMPEX,TMPWT,MTMP ELSE READ(ISU,END=9000) JTOT,INRG,ECHK,ITMPEX,TMPWT,MTMP,NOPEN ENDIF ENDIF C JTOT.LT.0 IS AN (OBSOLETE) "ROLOUT" MARKER IF(JTOT.LT.0) GOTO 2000 C C CHECK WHETHER THIS ENERGY IS PROCESSED. I.E. INRG=ICALC(NPASS). C NOTE: I DID NOT TEST THE MULTIPLE ENERGY FEATURE OF THIS CODE WHEN C I PORTED IT. (MET) C USUAL CHECK THAT ABS(ECHK-ENERGY(INRG)).LT.EPS NOT DONE HERE C IF (IPRINT.GT.5) THEN IF (IPROGM.LT.14) THEN WRITE(6,681) JTOT,MTMP,INRG,ECHK 681 FORMAT(' READ JTOT',I5,'.',I2,' ENERGY(',I2,') =',F12.3) ELSE WRITE(6,682) JTOT,MTMP,INRG,ECHK,NOPEN 682 FORMAT(' READ JTOT',I5,'.',I2,' ENERGY(',I2,') =',F12.3, 1 ' NOPEN =',I4) ENDIF ENDIF IF(INRG.EQ.ICALC(NPASS)) GOTO 1700 C C IF THIS ENERGY IS NOT USED, WE MUST SKIP OVER THIS SET OF DATA. C METHOD USED DEPENDS ON LFMT AND IPROGM IF (LFMT) THEN C SAVE TIME BY NOT CONVERTING VALUES, JUST SKIP CARDS READ(ISU,114,END=9002) NOPEN NSQ=NOPEN*NOPEN NCARD=2*((NSQ+4)/5)+NOPEN IF (IPRINT.GT.5) THEN WRITE(6,683) NOPEN,NCARD 683 FORMAT(' ---- SKIPPING THIS ENERGY. NOPEN/NCARD =',2I5) ENDIF DO 205 I=1,NCARD 205 READ(ISU,105,END=9002) ELSE IF (IPROGM.LT.14) THEN READ(ISU,END=9002) NOPEN,(LEV(1),LORB(1),WVEC(1),I=1,NOPEN) ELSE READ(ISU,END=9002) (LEV(1),LORB(1),WVEC(1),I=1,NOPEN) ENDIF NSQ=NOPEN*NOPEN IF (IXTM+NSQ.GT.MXTM) THEN WRITE(6,*) ' *** INADEQUATE STORAGE TO READ UNUSED S-MATRICES' WRITE(6,*) ' WILL TREAT AS AN EOF CONDITION' GO TO 9002 ENDIF CALL SREAD(ISU,NOPEN,TR(IXTM+1),IEND) IF (IEND.GT.0) GO TO 9002 CALL SREAD(ISU,NOPEN,TI(IXTM+1),IEND) IF (IEND.GT.0) GO TO 9002 ENDIF GOTO 2000 C C REACH HERE IF WE ARE CALCULATING CURRENT INRG 1700 NSMAT=NSMAT+1 IF( ICOUNT .LE. MXJT ) GOTO 1600 C IF(NSMAT.LE.MXJT) GOTO 1600 WRITE(6,694) MXJT 694 FORMAT('0 * * * ERROR. TOO MANY S-MATRICES / ENERGY FOR DIMENSION & MXJT =',I8) GOTO 9001 C 1600 IF (IPROGM.GE.14) THEN IF (NOPEN.GT.MXCH) THEN WRITE(6,*) ' *** NOPEN EXCEEDS MXCH',NOPEN,MXCH WRITE(6,*) ' JTOT,INRG,M',JTOT,INRG,MTMP GO TO 9001 ENDIF READ(ISU,END=9002) (LEV(I),LORB(I),WVEC(I),I=1,NOPEN) ELSE IF (LFMT) THEN READ(ISU,104,END=9002) NOPEN, 1 (LEV(I),LORB(I),WVEC(I),I=1,NOPEN) ELSE READ(ISU,END=9002) NOPEN,(LEV(I),LORB(I),WVEC(I),I=1,NOPEN) ENDIF IF(NOPEN.GT.MXCH) THEN WRITE(6,*) ' *** NOPEN EXCEEDS MXCH',NOPEN,MXCH C TERMINATE SINCE STORAGE IS PROBABLY CORRUPTED BY READ STOP ENDIF ENDIF C C SET POINTERS TO THIS SET OF INPUT S-MATRICES. . . C N.B. ISTBS AND ISTTM ARE 1 LESS THAN 1ST INDEX IN VECTOR. JT(ICOUNT) = JTOT NBS(ICOUNT) = NOPEN ISTBS(ICOUNT) = IXBS ISTTM(ICOUNT) = IXTM C DO 1400 I=1,NOPEN IXBS=IXBS+1 IF(IXBS.LE.MXBS) GOTO 1401 WRITE(6,695) MXBS 695 FORMAT('0 * * * ERROR. MXBS EXCEEDED.',I8) GOTO 9001 1401 J(IXBS)=LEV(I) L(IXBS)=LORB(I) 1400 WV(IXBS)=WVEC(I) C NSQ=NOPEN*NOPEN NMAX=MAX0(NMAX,NOPEN) C SG (NOV 94): I DON'T UNDERSTAND USE OF NMAX OR CALCULTION BELOW IF(IXTM+NSQ+NOPEN*NMAX .LE. MXTM) GOTO 1501 WRITE(6,696) MXTM 696 FORMAT('0 * * * ERROR. MXTM EXCEEDED.',I8) GOTO 9001 1501 IF (LFMT) THEN READ(ISU,105,END=9002) (TR(IXTM+I),I=1,NSQ) READ(ISU,105,END=9002) (TI(IXTM+I),I=1,NSQ) ELSE CALL SREAD(ISU,NOPEN,TR(IXTM+1),IEND) IF (IEND.GT.0) GO TO 9002 CALL SREAD(ISU,NOPEN,TI(IXTM+1),IEND) IF (IEND.GT.0) GO TO 9002 ENDIF IXTM=IXTM+NSQ IMCH=IXTM+1 C C CORRECT (ARTHURS-DALGARNO) S-MATRICES AS PER HUBY. C C CALL HUBY(NBS(NSMAT),L(ISTBS(NSMAT)+1), C X TR(ISTTM(NSMAT)+1),TI(ISTTM(NSMAT)+1)) CALL HUBY(NBS(ICOUNT), L(ISTBS(ICOUNT)+1), & TR(ISTTM(ICOUNT)+1), TI(ISTTM(ICOUNT)+1)) C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C NOW READY TO PROCESS TR, TI TO SIG. C C FIRST, PROCESS WITH SELF, I.E. 'DIAG' T1 = T2. C IXB1=ISTBS(NSMAT)+1 C IXT1=ISTTM(NSMAT)+1 C JT1=JT(NSMAT) C NB1=NBS(NSMAT) C IXB1 = ISTBS(ICOUNT) + 1 IXT1 = ISTTM(ICOUNT) + 1 JT1 = JT(ICOUNT) NB1 = NBS(ICOUNT) C IC=0 CALL GCLOCK(XLTIME) CALL SBEADD(JT1,NB1,J(IXB1),L(IXB1),WV(IXB1),TR(IXT1),TI(IXT1), 1 JT1,NB1,J(IXB1),L(IXB1),WV(IXB1),TR(IXT1),TI(IXT1), 2 IC,JLEV,NLEV,SIG,MXLV,KSET,NKSET,TR(IMCH),TI(IMCH), 3 TTIME(3),IPRINT,1,KDIFMX) CALL GCLOCK(XMTIME) TTIME(2)=TTIME(2)+(XMTIME-XLTIME) C C NEXT, PROCESS WITH PREVIOUSLY STORED T-MATRICES. C NS2=NSMAT NS2 = ICOUNT 9100 NS2=NS2-1 IF(NS2.LE.0) GOTO 9200 IXB2=ISTBS(NS2)+1 IXT2=ISTTM(NS2)+1 JT2=JT(NS2) NB2=NBS(NS2) IC = ICNTRL CALL GCLOCK(XLTIME) CALL SBEADD(JT1,NB1,J(IXB1),L(IXB1),WV(IXB1),TR(IXT1),TI(IXT1), 1 JT2,NB2,J(IXB2),L(IXB2),WV(IXB2),TR(IXT2),TI(IXT2), 2 IC,JLEV,NLEV,SIG,MXLV,KSET,NKSET,TR(IMCH),TI(IMCH), 3 TTIME(3),IPRINT,1,KDIFMX) C C IF IC=1 THEN DOUBLING IS NOT PERFORMED AND YOU MUST CALL SBEADD AGAI C TO HANDLE TO CASE FOR JT1 AND JT2 SWITCHED C C FROM WHAT I CAN SEE, DOUBLING WILL ALWAYS WORK FOR THE CROSS SECTION C THAT THIS CODE CALCULATES. THE OPTION FOR "NOT DOUBLING" IS LEFT IN C SO THAT THIS CAN BE TESTED. (MET) C IF( IC .EQ. 1 ) THEN C C GET CONTRIBUTION FROM JT1/JT2 SWITCHED C *** C *** THIS VERSION JUST DOUBLES THE CONTRIBUTION FOR FIRST CALL C *** TO SBEADD: CORRECT IF KJ AND KJ1 HAVE SAME PARITY C *** KJ+KJ1 ODD GIVES PURE IMAGINARY SIGMA, AND PROGRAM NEEDS C *** MODIFICATION. THESE WILL BE NEEDED FOR SOME TRANSPORT CASES. C *** C MODIFICATION TO TAKE INTO ACCOUNT (KJ+KJ1) ODD HAS BEEN MADE C (IT ACTUALLY OCCURS IN SUBROUTINE SBEADD) AND PROGRAM AUTO C SWITCHES BETWEEN REAL AND IMAGINARY CALCULATION, AS NEEDED. C IN OUTPUT, THE IMAGINARY CROSS SECTIONS ARE GIVEN AS REAL NUM'S C AND THESE NUMBERS SHOULD BE MULTIPLIED BY I TO GIVE THE C CORRECT PURE IMAGINARY RESULT. C CALL SBEADD(JT2,NB2,J(IXB2),L(IXB2),WV(IXB2),TR(IXT2),TI(IXT2), 1 JT1,NB1,J(IXB1),L(IXB1),WV(IXB1),TR(IXT1),TI(IXT1), 2 IC,JLEV,NLEV,SIG,MXLV,KSET,NKSET,TR(IMCH),TI(IMCH), 3 TTIME(3),IPRINT,1,KDFIMX) ENDIF CALL GCLOCK(XMTIME) TTIME(2)=TTIME(2)+(XMTIME-XLTIME) C C ALWAYS GO BACK FOR MORE TOTAL J VALUES. NOTE: IF THE J VALUES C ARE IN ASCENDING ORDER THEN THERE WILL BE A POINT WHERE THE C CONTRIBUTIONS FROM LOWER J VALUES ENDS AND GOING BACK IS SIMPLY C SPINNING YOUR WHEELS. SBEADD POPS BACK RIGHT AWAY IF THIS IS THE CASE. C THIS IS A BIT INEFFICIENT BUT SAFER, I SUPPOSE. C IT DOES TAKE A LOT OF MEMORY THOUGH. C THIS IS AN OLD COMMENT AND THE CURRENT VERSION SAVES MEMORY BY C ONLY KEEPING THE LAST 10 ARRAYS. C GOTO 9100 C C 9200 CONTINUE CALL GCLOCK(XJTIME) TIME=XJTIME-XITYME XITYME=XJTIME IF(IPRINT.GT.10) WRITE(6,614) JTOT,TIME 614 FORMAT(' TIME FOR JTOT =',I4,' WAS',F8.2,' SECONDS.') C C CHECK TO SEE IF ARRAY VALUES NEED TO BE SHIFTED C IF( ICOUNT .LT. MXJT ) THEN ICOUNT = ICOUNT + 1 ELSE ITMP1 = ISTTM(2) ITMP2 = ISTBS(2) DO 800 I = 1, MXJT-1 JT(I) = JT(I+1) NBS(I) = NBS(I+1) ISTTM(I) = ISTTM(I+1) - ITMP1 ISTBS(I) = ISTBS(I+1) - ITMP2 800 CONTINUE IXTM = IXTM - ITMP1 IXBS = IXBS - ITMP2 DO 810 I = 1, IXBS J(I) = J(ITMP2+I) L(I) = L(ITMP2+I) WV(I) = WV(ITMP2+I) 810 CONTINUE DO 820 I = 1, IXTM TR(I) = TR(ITMP1+I) TI(I) = TI(ITMP1+I) 820 CONTINUE ENDIF C GO BACK TO READ MORE INPUT S-MATRICES. GOTO 2000 C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C 9002 WRITE(6,688) 688 FORMAT('0'/'0 PREMATURE END OF FILE ENCOUNTERED ON TAPE.') GOTO 9000 C C OUT OF STORAGE BEFORE END OF FILE 9001 WRITE(6,698) 698 FORMAT('0'/'0 * * * WARNING.'/'0 * * * OUT OF STORAGE BEFORE ' 1 ,'END OF FILE. PROCESSING TERMINATED.'/ 2 '0 * * * WARNING.') C C EOF ON INPUT TAPE. CALC. COMPLETED FOR ENERGY(INRG=NPASS). 9000 WRITE(6,610) LABEL,NSMAT,ICALC(NPASS),ENERGY(ICALC(NPASS)) 610 FORMAT('0 S-MATRICES TAKEN FROM SAVE TAPE WITH LABEL =',20A4/ A '0 E.O.F. ON TAPE AFTER',I5, 2 ' S-MATRICES PROCESSED FOR ENERGY(',I3,') =',F11.4,' (1/CM).') WRITE(6,612) IXBS,MXBS,IXTM,MXTM 612 FORMAT('0 CHANNEL BASIS FUNCTIONS USED',I6,' OF',I6,' WORDS'/ & '0 S-MATRIX ELEMENTS USED',I7,' OF',I7,' WORDS') C CALL GCLOCK(XJTIME) TIME=XJTIME-XITIME WRITE(6,613) TIME 613 FORMAT('0'/'0 TIME FOR THIS ENERGY SET WAS',F9.2,' SECONDS.') C IF(MAXSIG.GT.0) WRITE(6,603) MAXSIG C C OUTPUT SIGMA MATRICES // CALL OUTSBE(NLEV,JLEV,MXLV,SIG,NKSET,KSET,0,6) ENE=ENERGY(ICALC(NPASS)) IF (MXN .GT. 0) THEN IF (NPASS .EQ. 1) THEN WRITE(MXN,16) NKSET 16 FORMAT(1X,I3) DO 12 I=1,NKSET 12 WRITE(MXN,17) (KSET(IA,I),IA=1,5) 17 FORMAT(1X,5I3) ENDIF WRITE(MXN,15) ENE NSQ1=NLEV*NQN JMAX1=JLEV(JLEV(NSQ1-NCCH)) JMIN1=JLEV(1) JSTEP1=JLEV(2)-JLEV(1) IF( JSTEP1 .EQ. 0 ) JSTEP1 = 1 WRITE(MXN,19) JMIN1,JMAX1,JSTEP1 15 FORMAT(1X,D18.10) 19 FORMAT(1X,3I3) DO 13 LF=1,MXLV-NCCH JF=JLEV(LF) DO 13 LI=1,MXLV-NCCH JI=JLEV(LI) 13 WRITE(MXN,14) JI,JF,(Y(I,LF,LI),I=1,NKSET) 14 FORMAT(1X,2I3,11D18.10) ENDIF CALL OUTP(ENE,NLEV,JLEV,MXLV,NKSET,KSET,IPARTU) C C IF ALL ENERGIES HAVE NOT BEEN DONE, GO BACK FOR NEXT. 1900 NPASS=NPASS+1 IF(NPASS.LE.NCALC) GOTO 1000 C RETURN END FUNCTION SIXJN(J1,J2,J5,J4,J3,J6,NEXTJ1,L) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE XJ6J,XJ1MIN,IVAL PARAMETER (MXDIM=200) LOGICAL NEXTJ1 DIMENSION XJ6J(MXDIM,2),IVAL(2),XJ1MIN(2) IF(NEXTJ1) GOTO 100 IVAL(L)=MXDIM CALL J6J( DBLE(J2),DBLE(J3), 1 DBLE(J4),DBLE(J5),DBLE(J6), 3 IVAL(L),XJ1MIN(L),XJ6J(1,L)) 100 IND=1+J1-INT(XJ1MIN(L)+0.1D0) SIXJN=0.D0 IF(IND.GE.1 .AND. IND.LE.IVAL(L)) SIXJN=XJ6J(IND,L) RETURN END FUNCTION SIXJ(J1,J2,J5,J4,J3,J6) C C CALCULATES 6-J SYMBOL: _(J1 J2 J3 )_ C (J4 J5 J6 ) C INTERFACE TO J6J ROUTINE. C MODIFIED BY S. GREEN 20 AUG 93; PASS DIMENSION OF XJ6J FOR CHECKING C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(MXDIM=200) DIMENSION XJ6J(MXDIM) IVAL=MXDIM CALL J6J( DBLE(J2),DBLE(J3), 1 DBLE(J4),DBLE(J5),DBLE(J6), 3 IVAL,XJ1MIN,XJ6J) IND=1+J1-INT(XJ1MIN+0.1D0) SIXJ=0.D0 IF(IND.GE.1 .AND. IND.LE.IVAL) SIXJ=XJ6J(IND) RETURN END SUBROUTINE J6J(J2,J3,L1,L2,L3,IVAL,J1MIN,D6J) IMPLICIT DOUBLE PRECISION (A-H,J-Z) DIMENSION D6J(2) DATA ZERO/0.D0/,TENTH/0.1D0/,HALF/0.5D0/,ONE/1.D0/,TWO/2.D0/, $ CONST/1.0D-12/ E(J1S)=SQRT((J1S-MJ23S)*(J23P1S-J1S)*(J1S-ML23S)*(L23P1S-J1S)) F(J1,JJP1)=(TWO*J1+ONE)*(JJP1*(FACT-JJP1-TWO*LLP1)+FACT2) C C THIS ROUTINE CALCULATES THE 6-J COEFFICIENTS FOR ALL PERMISSIBLE C VALUES OF J1 FOR FIXED VALUES OF J2, J3, L1, L2, AND L3 USING THE C RECURSIVE ALGORITHM OF K. SCHULTEN AND R. G. GORDON, J. MATH. PHYS. C VOL. 16, P. 1961, (1975). C PROGRAMMED BY D. E. FITZ, 10/22/79 C MODIFIED BY S. GREEN 20 AUG 93 TO TEST DIMENSION ON D6J C MXDIM=IVAL JJP2=J2*(J2+ONE) JJP3=J3*(J3+ONE) LLP1=L1*(L1+ONE) LLP2=L2*(L2+ONE) LLP3=L3*(L3+ONE) MJ23S=(J2-J3)**2 ML23S=(L2-L3)**2 J23P1S=(J2+J3+ONE)**2 L23P1S=(L2+L3+ONE)**2 FACT2=(LLP2-LLP3)*(JJP2-JJP3) FACT=JJP2+JJP3+LLP2+LLP3 J1MIN=MAX(ABS(J2-J3),ABS(L2-L3)) J1MAX=MIN(J2+J3,L2+L3) IVAL=INT(J1MAX-J1MIN+ONE+TENTH) IF (IVAL.GT.MXDIM) THEN WRITE(6,*) 'J6J: ARRAY D6J TOO SMALL. NEEDS ',IVAL,' BUT ONLY ', 1 MXDIM,' SUPPLIED' STOP ENDIF C C TEST FOR OTHER TRIANGULAR INEQUALITES. C IL1=INT(TWO*L1+TENTH) IL2=INT(TWO*L2+TENTH) IL3=INT(TWO*L3+TENTH) IJ2=INT(TWO*J2+TENTH) IJ3=INT(TWO*J3+TENTH) IF((IJ2.LE.IL1+IL3.AND.IJ2.GE.IABS(IL1-IL3)).AND. $ (IJ3.LE.IL1+IL2.AND.IJ3.GE.IABS(IL1-IL2))) GO TO 11 DO 12 I=1,IVAL 12 D6J(I)=ZERO RETURN C 11 INMID=(IVAL+3)/2 SGNV=J2+J3+L2+L3 SGN=ONE ISIGN=INT(SGNV+TENTH) IF(MOD(ISIGN,2).NE.0) SGN=-ONE D6J(1)=HALF C C UPWARD RECURSION. C IF(IVAL.EQ.1) GO TO 40 JJP1=J1MIN*(J1MIN+ONE) F1=F(J1MIN,JJP1) J1=J1MIN+ONE J1S=J1*J1 E2=E(J1S) IF(J1MIN.LT.TENTH) GO TO 15 D6J(2)=-F1*D6J(1)/(E2*J1MIN) GO TO 16 15 D6J(2)=-HALF*(LLP2+JJP2-LLP1)*D6J(1)/SQRT(JJP2*LLP2) 16 SCALE=D6J(2) IF(IVAL.EQ.2) GO TO 40 DO 21 IJ2=3,INMID JJP1=J1*(J1+ONE) F1=F(J1,JJP1) J1=J1+ONE E1=E2 J1S=J1*J1 E2=E(J1S) 21 D6J(IJ2)=-(F1*D6J(IJ2-1)+J1*E1*D6J(IJ2-2))/(E2*(J1-ONE)) SCALE=D6J(INMID) IEXC=5 IF(ABS(SCALE).GT.CONST) GO TO 18 INMID=INMID-1 SCALE=D6J(INMID) IEXC=3 GO TO 30 18 IF(IVAL.EQ.3) GO TO 40 C C DOWNWARD RECURSION. C 30 D6J(IVAL)=HALF J1=J1MAX J1S=J1*J1 JJP1=J1*(J1+ONE) F1=F(J1,JJP1) E1=E(J1S) D6J(IVAL-1)=-F1*D6J(IVAL)/(E1*(J1+ONE)) IEND=IVAL-INMID IF(IVAL.LE.IEXC) GO TO 31 DO 32 IJ2=2,IEND J1=J1-ONE E2=E1 J1S=J1*J1 JJP1=J1*(J1+ONE) E1=E(J1S) F1=F(J1,JJP1) 32 D6J(IVAL-IJ2)=-(J1*E2*D6J(IVAL-IJ2+2)+F1*D6J(IVAL-IJ2+1))/ $ (E1*(J1+ONE)) C C MATCH UPWARD AND DOWNWARD RECURSIVE RESULTS BY SCALING. C 31 SCALE=SCALE/D6J(INMID) DO 33 IJ2=INMID,IVAL 33 D6J(IJ2)=SCALE*D6J(IJ2) C C NORMALIZE RESULTS AND SET PHASE. C 40 SUM=ZERO DO 41 IJ2=1,IVAL J1=J1MIN+DBLE(IJ2-1) 41 SUM=SUM+(TWO*J1+ONE)*D6J(IJ2)**2 RNORM=ONE/SQRT(SUM*(TWO*L1+ONE)) IF((SGN*D6J(IVAL)).LT.ZERO) RNORM=-RNORM DO 42 IJ2=1,IVAL 42 D6J(IJ2)=D6J(IJ2)*RNORM RETURN END FUNCTION THREEJ (J1,J2,J3) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C COMPUTATION OF SPECIAL WIGNER 3J COEFFICIENT WITH C VANISHING PROJECTIONS. SEE EDMONDS, P. 50. C C THIS VERSION EVALUATES BINOM AND PARITY IN-LINE C SHOULD IMPROVE EFFICIENCY, ESPECIALLY ON CRAY; C ALSO GIVES IMPROVEMENT ON AMDAHL (SG: 20 DEC 92) C C STATEMENT FUNCTION FOR DELTA ASSOCIATED W/ RACAH AND SIXJ SYMBOLS C DELTA(I,J,K)= SQRT(1.D0/ ( BINOM(I+J+K+1,I+J-K) * C 1 BINOM(K+K+1,I-J+K) * DBLE(K+J-I+1) ) ) C I1=J1+J2+J3 IF (I1-2*(I1/2).NE.0) GO TO 8 1 I2=J1-J2+J3 IF (I2) 8,2,2 2 I3=J1+J2-J3 IF (I3) 8,3,3 3 I4=-J1+J2+J3 IF (I4) 8,4,4 4 I5=I1/2 I6=I2/2 SIGN=1.D0 IF (I5-2*(I5/2).NE.0) SIGN=-SIGN C 7 THREEJ=SIGN*DELTA(J1,J2,J3)*BINOM(I5,J1)*BINOM(J1,I6) C B1,B2 ARE BINOM ASSOCIATED W/ DELTA N=J1+J2+J3+1 M=J1+J2-J3 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B1=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 101 I = 1,MNM F = F+1.D0 C = (FN-F)*B 101 B = C/F B1 = B ENDIF N=J3+J3+1 M=J1-J2+J3 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B2=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 102 I = 1,MNM F = F+1.D0 C = (FN-F)*B 102 B = C/F B2 = B ENDIF DELTA=SQRT(1.D0/(B1*B2*(J3+J2-J1+1))) C B3=BINOM(I5,J1), B4=BINOM(J1,I6) N=I5 M=J1 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B3=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 103 I = 1,MNM F = F+1.D0 C = (FN-F)*B 103 B = C/F B3 = B ENDIF N=J1 M=I6 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B4=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 104 I = 1,MNM F = F+1.D0 C = (FN-F)*B 104 B = C/F B4 = B ENDIF THREEJ=SIGN*DELTA*B3*B4 RETURN 8 THREEJ=0.D0 RETURN END FUNCTION PARITY(I) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARITY=1.D0 IF((I/2)*2-I.NE.0) PARITY=-1.D0 RETURN END FUNCTION XNINEJ(IX1,IY1,IZ1,IX2,IY2,IZ2,IX3,IY3,IZ3) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XJ9J(200) DATA MXDIM/200/ C IVAL=MXDIM CALL J9J(DBLE(IX1),DBLE(IY1), 1 DBLE(IX2),DBLE(IY2),DBLE(IZ2), 2 DBLE(IX3),DBLE(IY3),DBLE(IZ3), 3 IVAL,Z1MIN,XJ9J) IND=1+IZ1-INT(Z1MIN+0.1D0) XNINEJ=0.D0 IF(IND.GE.1 .AND. IND.LE.IVAL) XNINEJ=XJ9J(IND) RETURN END SUBROUTINE J9J(J1,J2,J4,J5,J6,J7,J8,J9,IVAL,J3MIN,D9J) IMPLICIT DOUBLE PRECISION (A-H,J-Z) DIMENSION D9J(1),D6J3(200),D6J5(200),D6J7(200) DATA MXDIM6/200/ DATA ZERO/0.D0/,TENTH/0.1D0/,HALF/0.5D0/,ONE/1.D0/,TWO/2.D0/ C C THIS ROUTINE CALCULATES THE 9-J SYMBOLS BY SUMMATION OVER 6-J C SYMBOLS WHICH IN TURN ARE CALCULATED BY THE RECURSIVE METHOD C OF SCHULTEN AND GORDON. C PROGRAMMED BY D. E. FITZ, 22 OCT 79. C C MODIFIED BY M. L. DUBERNET, 15 SEP 93 AND J. M. HUTSON, 3 OCT 93 C TO ALLOW HALF-INTEGER ANGULAR MOMENTA C MODIFIED BY J. M. HUTSON, 3 OCT 93 TO CHECK D9J DIMENSION C MXDIM9=IVAL J3MIN=MAX(ABS(J1-J2),ABS(J6-J9)) J3MAX=MIN(J1+J2,J6+J9) IJ3N=INT(TWO*J3MIN+TENTH) IJ3X=INT(TWO*J3MAX+TENTH) IVAL=1+(IJ3X-IJ3N)/2 IF(IVAL.GT.MXDIM9) THEN WRITE(6,*) 'J9J: ARRAY D9J TOO SMALL. NEEDS ',IVAL,' BUT ONLY ', 1 MXDIM9,' SUPPLIED' STOP ENDIF C C TEST FOR TRIANGULAR INEQUALITIES. C D9J(1)=ZERO IF (IVAL.LE.0) RETURN IJ1=INT(TWO*J1+TENTH) IJ2=INT(TWO*J2+TENTH) IJ4=INT(TWO*J4+TENTH) IJ5=INT(TWO*J5+TENTH) IJ6=INT(TWO*J6+TENTH) IJ7=INT(TWO*J7+TENTH) IJ8=INT(TWO*J8+TENTH) IJ9=INT(TWO*J9+TENTH) DO 15 IJL=1,IVAL 15 D9J(IJL)=ZERO IF((IJ4-IABS(IJ7-IJ1))*(IJ7+IJ1-IJ4).LT.0) RETURN IF((IJ5-IABS(IJ8-IJ2))*(IJ8+IJ2-IJ5).LT.0) RETURN IF((IJ5-IABS(IJ6-IJ4))*(IJ6+IJ4-IJ5).LT.0) RETURN IF((IJ8-IABS(IJ9-IJ7))*(IJ9+IJ7-IJ8).LT.0) RETURN C IVAL7=MXDIM6 CALL J6J(J1,J9,J7,J8,J4,IVAL7,JMIN7,D6J7) C IVAL5=MXDIM6 CALL J6J(J6,J2,J5,J8,J4,IVAL5,JMIN5,D6J5) C JMIN=MAX(JMIN5,JMIN7) JMAX=MIN(J1+J9,J2+J6,J4+J8) IEND=INT(JMAX-JMIN+TENTH+ONE) I5=INT(JMIN-JMIN5+TENTH) I7=INT(JMIN-JMIN7+TENTH) C C LOOP RUNS OVER TWICE J3 TO ALLOW HALF-INTEGER VALUES C ITAB=0 DO 20 IJ3=IJ3N,IJ3X,2 ITAB=ITAB+1 J3=HALF*DBLE(IJ3) C IVAL3=MXDIM6 CALL J6J(J1,J9,J3,J6,J2,IVAL3,JMIN3,D6J3) I3=INT(JMIN-JMIN3+TENTH) SUM=ZERO C DO 10 I=1,IEND J=DBLE(I-1)+JMIN SGN=ONE ISIGN=INT(TWO*J+TENTH) IF(MOD(ISIGN,2).NE.0) SGN=-ONE SUM=SUM+(TWO*J+ONE)*SGN*D6J3(I+I3)*D6J5(I+I5)*D6J7(I+I7) 10 CONTINUE D9J(ITAB)=SUM 20 CONTINUE RETURN END SUBROUTINE SWRITE(IU,N,S) C THIS IS THE MOLSCAT ROUTINE FOR READING/WRITING THE S-MATRIX IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION S(N,N) WRITE(IU) ((S(I,J),J=1,I),I=1,N) RETURN C ENTRY SREAD(IU,N,S,IEND) IEND=0 READ(IU,END=9999) ((S(I,J),J=1,I),I=1,N) DO 1000 I=1,N DO 1000 J=1,I-1 1000 S(J,I)=S(I,J) RETURN 9999 IEND=1 RETURN END subroutine ioinit c dummy this for rs/6000 (sg, nov 94) return end