C PROGRAM TO COMPUTE S=0,1 CROSS SECTIONS OF GENERALIZED HESS METHOD C VERSION FOR TWO S-MATRIX SAVE TAPES, ONE FOR EACH VIB LEVEL C VERSION 3 (MULTIPLE ENERGY CALCULATIONS) 2 DEC 94, SG C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C STORAGE FOR S-MATRICES, L-VALUES, POINTERS; LAST DIMENSION C IS FOR INITIAL-FINAL SPECTRAL LEVEL PARAMETER (MXS=70000, MXJT=300, MXCL=25) DIMENSION SR(MXS),SI(MXS),L1(MXS),L2(MXS) DIMENSION JT(MXJT,MXCL,2),IST(MXJT,MXCL,2),NOP(MXJT,MXCL,2), 1 ICALC(MXCL,2),INX(MXCL,2), WV(MXCL),ECALC(MXCL), 2 JACCMX(MXCL),WAR(MXCL),WAI(MXCL),WRR(MXCL),WRI(MXCL) LOGICAL LWV(MXCL) C C STORAGE FOR SAVE TAPE VALUES, SECOND DIMENSION IS FOR TWO TAPES PARAMETER (MXCH=200,MXJLEV=500) DIMENSION J(MXCH),L(MXCH),WVEC(MXCH),IXB(MXCH),JLEV(MXJLEV,2) DIMENSION SREAL(MXCH*MXCH),SIMAG(MXCH*MXCH) C DIMS OF ELEVEL (IN COMMON/CMBASE/) AND ENERGY FIXED BY MOLSCAT PARAMETER (MXEL=1000,MXNRG=100) DIMENSION ENERGY(MXNRG,2),ELEVEL(MXEL,2) DIMENSION ITYPE(2),URED(2),NLEV(2),NQN(2),NNRG(2),IP(2),NLEVEL(2), 1 IIN(2) CHARACTER*4 LABEL(20,2) LOGICAL LFMT,EQUAL,LEOF(2) INTEGER PRNTLV C C LSPEC(2) HAS THE 'SPECTRAL' LEVEL NOS. CORRESPONDING TO JLEV(,NQN) DIMENSION LSPEC(2) C C C NAMELIST PARAMETERS: C IIN - UNIT OF S-MATRIX SAVE TAPES (DEFAULT TO 10, 11) C LFMT - FORMATTED (.TRUE.) / UNFORMATTED (.FALSE.) SAVE TAPE C NCALC - NUMBER OF ENERGIES AT WHICH TO CALCULATE C ECALC - ENERGIES AT WHICH TO CALCULATE C LSPEC - LEVEL NUMBERS OF INITIAL/FINAL SPECTRAL LEVEL C IQ - TENSOR ORDER (DEFAULTS TO ABS(JIN-JFN)) C PRNTLV - PRINT LEVEL CONTROL C .GE.3 ACCUMULATED CROSS SECTIONS AFTER EACH JTOT C .GE.5 ECHO INPUT SETS USED IN CALC; PARTIAL JTOT CONTRIBS C .GE.10 ECHO ALL INPUT SETS C NAMELIST /GHM/IIN,LFMT,NCALC,ECALC,PRNTLV,LSPEC,IQ C DATA I00/0/, I01/1/ C C FORMATS FOR SAVE TAPE (FOR VERSION 11 AND EARLIER) 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) ) 105 FORMAT(5E16.8) C C STATEMENT FUNCTION EQUAL(X,Y)=ABS(X-Y).LT.1.D-5 C C ----------------- FIRST EXECUTABLE STATEMENT -------------- C WRITE(6,600) 600 FORMAT(' PROGRAM TO CALCULATE GENERALIZED HESS METHOD', 1 ' CROSS SECTIONS') C C INITIALIZE AND SET DEFAULT VALUES FOR NAMELIST CONSTANTS CALL GCLOCK(TSTART) PI=ACOS(-1.D0) PRNTLV=3 IIN(1)=10 IIN(2)=11 IQ=-1 LFMT=.FALSE. ISTOP=0 ISTOPX=0 NCALC=-1 DO 1001 IE=1,MXCL ECALC(IE)=0.D0 JACCMX(IE)=-1 LWV(IE)=.FALSE. WAR(IE)=0.D0 WAI(IE)=0.D0 WRR(IE)=0.D0 1001 WRI(IE)=0.D0 DO 1002 IA=1,2 LEOF(IA)=.FALSE. DO 1003 IE=1,MXEL 1003 ELEVEL(IE,IA)=0.D0 DO 1002 IE=1,MXCL ICALC(IE,IA)=0 INX(IE,IA)=0 DO 1002 I=1,MXJT IST(I,IE,IA)=0 JT(I,IE,IA)=0 1002 NOP(I,IE,IA)=0 C READ(5,GHM,END=1000) C 1000 IF (NCALC.LE.0) THEN WRITE(6,*) ' *** NCALC.LE.0 *** MUST INPUT NCALC,ECALC' STOP ENDIF WRITE(6,601) IIN,LSPEC,(IE,ECALC(IE),IE=1,NCALC) 601 FORMAT(' FROM MOLSCAT S-MATRICES ON UNITS',2I4/ 1 ' REQUESTED SPECTRAL LINE BETWEEN LEVELS',2I4/ 2 ' REQUESTED COLLISION KINETIC ENERGIES (1/CM) ARE'/ 3 (8X,'ENERGY(',I3,') =',F12.4)) WRITE(6,*) IF (LFMT) THEN WRITE(6,*) ' LFMT SPECIFIES A FORMATTED SAVE TAPE' ELSE WRITE(6,*) ' LFMT SPECIFIES AN UNFORMATTED SAVE TAPE' ENDIF WRITE(6,*) ' PRINT LEVEL CONTROL (PRNTLV) IS',PRNTLV C DO 8000 IU=1,2 IT=IIN(IU) C GET AND PROCESS HEADER RECORD FROM SAVE TAPE IF (LFMT) THEN READ(IT,100,END=9900) (LABEL(I,IU),I=1,20),ITYPE(IU),NLEV(IU), 1 NQN(IU),URED(IU),IP(IU) ELSE READ(IT,END=9900) (LABEL(I,IU),I=1,20),ITYPE(IU),NLEV(IU), 1 NQN(IU),URED(IU),IP(IU) ENDIF WRITE(6,603) IT,(LABEL(I,IU),I=1,20),ITYPE(IU),NLEV(IU),NQN(IU), 1 URED(IU),IP(IU) 603 FORMAT(/' INPUT FROM UNIT (',I3,'):'/' LABEL =',20A4/ 2 ' ITYPE =',I3,5X,'NLEV, NQN =',2I4,5X,'URED =',F8.4, 3 5X,'IPROGM =',I3) IF (IP(IU).LT.3) THEN WRITE(6,*) ' *** NO IMPLEMENTATION FOR IPROGM.LT.3' STOP ENDIF IF (LFMT.AND.IP(IU).GE.11) THEN WRITE(6,*) ' *** IPROGM INCONSISTENT W/ REQUESTED TAPE FORMAT' STOP ENDIF IF (.NOT.(ITYPE(IU).EQ.1.OR.ITYPE(IU).EQ.2.OR. 1 ITYPE(IU).EQ.7)) THEN WRITE(6,*) ' *** CODE NOT IMPLEMENTED FOR ITYPE',ITYPE STOP ENDIF C C GET JLEV INFORMATION FROM TAPE NSQ=NLEV(IU)*NQN(IU) IF (NSQ.GT.MXJLEV) THEN WRITE(6,*) ' *** NLEV*NQN.GT.MXJLEV',MXJLEV STOP ENDIF IF (LFMT) THEN READ(IT,101) (JLEV(I,IU),I=1,NSQ) ELSE READ(IT) (JLEV(I,IU),I=1,NSQ) ENDIF C C GET NLEVEL/ELEVEL INFORMATION IF (LFMT) THEN IF (IP(IU).GE.3) READ(IT,102) NLEVEL(IU),(ELEVEL(I,IU), 1 I=1,NLEVEL(IU)) ELSE READ(IT) NLEVEL(IU),(ELEVEL(I,IU),I=1,NLEVEL(IU)) ENDIF WRITE(6,604) 604 FORMAT(' BASIS SET DATA:'/) DO 1005 I=1,NLEV(IU) 1005 WRITE(6,606) I,ELEVEL(ABS(JLEV(NLEV(IU)*(NQN(IU)-1)+I,IU)),IU), 1 (JLEV(NLEV(IU)*(K-1)+I,IU),K=1,NQN(IU)) 606 FORMAT(' LEVEL',I4,' ENERGY',F10.3,' QUANTUM NOS. ',10I4) C C GET COLLISION ENERGIES FROM SAVE TAPE IF (LFMT) THEN READ(IT,102) NNRG(IU),(ENERGY(I,IU),I=1,NNRG(IU)) ELSE READ(IT) NNRG(IU),(ENERGY(I,IU),I=1,NNRG(IU)) ENDIF IF (NNRG(IU).GT.MXNRG) THEN WRITE(6,*) ' *** NNRG.GT.MXNRG',NNRG(IU),MXNRG STOP ENDIF WRITE(6,607) NNRG(IU),(ENERGY(I,IU),I=1,NNRG(IU)) 607 FORMAT(/' NUMBER OF COLLISION ENERGIES ON SAVE TAPE IS',I4/ 1 (' ',1P,6D13.4)) WRITE(6,*) WRITE(6,608) 608 FORMAT(' ***') C SEE IF THERE IS AN APPROPRIATE ENERGY = ECALC + ELEVEL IE=1 2101 DO 2100 I=1,NNRG(IU) IF (EQUAL(ELEVEL(LSPEC(IU),IU)+ECALC(IE),ENERGY(I,IU))) THEN ICALC(IE,IU)=I WRITE(6,609) ECALC(IE),ICALC(IE,IU),ENERGY(ICALC(IE,IU),IU) 609 FORMAT(' *** CALC FOR KINETIC ENERGY =', 1 F12.4,' WILL USE E(',I3,') =',F12.4) GO TO 2102 ENDIF 2100 CONTINUE C NO MATCH. REMOVE FROM LIST, DECREASING NCALC WRITE(6,610) ECALC(IE) 610 FORMAT(' *** *** NO MATCHING TOTAL E FOR K.E. =',F12.4/ 1 ' *** *** REMOVED FROM LIST') NCALC=NCALC-1 IF (IE.LE.NCALC) THEN DO 2103 I1=IE,NCALC ECALC(I1)=ECALC(I1+1) DO 2103 IO=1,2 2103 ICALC(I1,IO)=ICALC(I1+1,IO) ENDIF IE=IE-1 2102 IE=IE+1 IF (IE.LE.NCALC) GO TO 2101 IF (NCALC.LE.O) THEN WRITE(6,*) ' *** SORRY. NO ENERGY MATCHES WERE FOUND' WRITE(6,*) ' CHECK INPUT FILES AND/OR &GHM ECALC' STOP ENDIF C CHECK THAT WE DON'T USE SAME ENERGY FOR 2 CALCS. IF (NCALC.LE.1) GO TO 8000 DO 2202 I1=1,NCALC-1 DO 2202 I2=I1+1,NCALC IF (ICALC(I1,IU).NE.ICALC(I2,IU)) GO TO 2202 WRITE(6,*) ' *** CANNOT USE USE SAME S-MATRICES FOR TWO ENERGIES',I1,I2 STOP 2202 CONTINUE 8000 WRITE(6,608) C THIS ENDS PROCESSING OF 'HEADER' INFORMATION FOR BOTH TAPES. C CHECK J-VALUES AND COMPARE WITH IQ JIN=JLEV(LSPEC(1),1) JFN=JLEV(LSPEC(2),2) WRITE(6,611) LSPEC(1),JIN,LSPEC(2),JFN 611 FORMAT(/' INITIAL LEVEL ',I3,' IS J(IN) =',I3/ 1 ' FINAL LEVEL ',I3,' IS J(FN) =',I3) IF (IQ.LT.0) THEN IQ=ABS(JIN-JFN) WRITE(6,612) IQ 612 FORMAT(' TENSOR ORDER (TAKEN FROM SPECTRAL LEVELS) =',I2) ELSE WRITE(6,613) 613 FORMAT(' TENSOR ORDER (TAKEN FROM &GHM IQ) =',I2) IF (IQ.NE.ABS(JIN-JFN)) WRITE(6,*) ' *** WARNING .NE. DELTA(J)' ENDIF C BELOW FOLLOWS 22 DEC 94 SG ANALYSIS THAT JDFMX IS GOVERNED BY C IQ+LAMBDA: NEED HERE LAMBDA=1 LIMIT. JDFMX0=ABS(IQ) JDFMX1=ABS(IQ)+1 JDFMAX=MAX(JDFMX0,JDFMX1) JTOLD=-1 C C START LOOP OVER JTOT / ENERGY VALUES ON SAVE TAPE C ALTERNATE TAPES EVERY TIME WE FIND A 'USEABLE' SET OF S-MATRICES C N.B. S-MATRICES FOR LEVEL LSPEC(IU) COME FROM UNIT IIN(IU) C 3000 DO 7000 IU=1,2 IF (LEOF(IU)) GO TO 7000 IT=IIN(IU) C FIRST, GET JTOT,INRG,M RECORD ----- 2000 IF (LFMT) THEN READ(IT,103,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M ELSE IF (IP(IU).LT.14) THEN READ(IT,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M ELSE READ(IT,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M,NOPEN ENDIF ENDIF C N.B. JTOT .LT. 0 INDICATES A OBSOLETE 'CHECKPOINT' RECORD. IF (JTOT.LT.0) THEN WRITE(6,*) ' *** ROLLOUT MARKER (OBSOLETE) ENCOUNTERED',IT,JTOT GO TO 2000 ENDIF C ECHO INPUT IF (PRNTLV.GE.10) THEN IF (IP(IU).LT.14) THEN WRITE(6,615) IT,JTOT,M,INRG,ECHK,IEXCH,WT 615 FORMAT(' UNIT(',I2,'): JTOT ',I4,'.',I3,' E(',I3,')=', 1 F10.3,' IEX,WT=',I2,F6.2) ELSE WRITE(6,616) IT,JTOT,M,INRG,ECHK,IEXCH,WT,NOPEN 616 FORMAT(' UNIT(',I2,'): JTOT ',I4,'.',I3,' E(',I3,')=', 1 F10.3,' IEX,WT,NOP=',I2,F7.2,I5) ENDIF ENDIF C C CHECK THAT ENERGY(INRG) CORRESPONDS WITH HEADER RECORD. . . IF (.NOT.EQUAL(ECHK,ENERGY(INRG,IU))) THEN WRITE(6,699) INRG,IT,ECHK 699 FORMAT(' *** WARNING.',I4,'-TH ENERGY ON UNIT(',I3,')', 1 ' ECHECK =',D16.8) ENDIF C C GET J,L,WVEC VALUES ----- IF (LFMT) THEN READ(IT,104,END=9100) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN) IF (NOPEN.GT.MXCH) THEN WRITE(6,698) MXCH,IT,JTOT,INRG,M 698 FORMAT(/' *** ERROR. NO. OF OPEN CHANNELS EXCEEDS',I6, 1 ' IT,JTOT,INRG,M =',I4,I7,3I4) STOP ENDIF ELSE IF (IP(IU).GE.14) THEN IF (NOPEN.GT.MXCH) THEN WRITE(6,698) MXCH,IT,JTOT,INRG,M STOP ENDIF READ(IT,END=9100) (J(I),L(I),WVEC(I),I=1,NOPEN) ELSE READ(IT,END=9100) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN) IF (NOPEN.GT.MXCH) THEN WRITE(6,698) MXCH,IT,JTOT,INRG,M STOP ENDIF ENDIF ENDIF C C FINALLY GET S-MATRICES --- NSQ=NOPEN*NOPEN IF (LFMT) THEN C WE SHOULD ARRANGE TO "SKIP" IF INRG.NE.ICALC C CF CODE IN /MOLV14/UTIL/SBE_SAVE.F READ(IT,105,END=9100) (SREAL(I),I=1,NSQ) READ(IT,105,END=9100) (SIMAG(I),I=1,NSQ) ELSE CALL SREAD(IT,NOPEN,SREAL,IEND) IF (IEND.GT.0) GO TO 9100 CALL SREAD(IT,NOPEN,SIMAG,IEND) IF (IEND.GT.0) GO TO 9100 ENDIF C C SEE IF WE NEED THIS ENERGY FOR AN ECALC(IE). DO 2001 IE=1,NCALC IF (ICALC(IE,IU).EQ.INRG) THEN IC=IE GO TO 2002 ENDIF 2001 CONTINUE C WE DON'T NEED THIS SET OF S-MATRICES; GO BACK FOR ANOTHER GO TO 2000 C C IF WE REACH HERE, WE WANT TO PROCESS THE (JTOT,INRG,M) BLOCK C CHOOSE THE ELEMENTS REQUIRED FOR SPECTRAL LEVELS LSPEC(IU) C N.B. INX IS INDEX NO. OF CURRENT SET -- ALWAYS .LT. MXJT 2002 INX(IC,IU)=INX(IC,IU)+1 IF (INX(IC,IU).GT.MXJT) THEN WRITE(6,*) ' *** NO. JTOT-S EXCEEDS MXJT',MXJT IF (PRNTLV.LT.3) WRITE(6,*) ' *** CURRENT JTOT =',JTOT C TRY TO PUSH DOWN STORAGE INX(IC,IU)=INX(IC,IU)-1 CALL PUSH(NCALC,INX,IST,JT,NOP,MXJT,MXCL,SR,SI,L1,L2,MXS, 1 JTOT,JDFMAX,ISTOP,PRNTLV) INX(IC,IU)=INX(IC,IU)+1 ENDIF IST(INX(IC,IU),IC,IU)=ISTOP JT(INX(IC,IU),IC,IU)=JTOT IF (JTOT.GT.JTOLD) THEN IF (JTOLD.NE.-1.AND.PRNTLV.GE.3) 1 WRITE(6,625) JTOLD,(IE,WAR(IE),WAI(IE),WRR(IE),WRI(IE), 2 IE=1,NCALC) 625 FORMAT(/' *** FINISHED JTOT =',I4,/ 1 (6X,'FOR E(',I3,'), S=0 IS',2F9.2,' S=1 IS',2F9.2)) JTOLD=JTOT ENDIF C C SET UP INDICES OF BASIS FNS WHICH MATCH LINE, LSPEC(IU) NB=0 DO 4000 I1=1,NOPEN IF (J(I1).NE.LSPEC(IU)) GO TO 4000 NB=NB+1 IXB(NB)=I1 IF (LWV(IC)) THEN IF (.NOT.EQUAL(WVEC(I1),WV(IC))) THEN WRITE(6,692) WVEC(I1),WV(IC),JTOT,INRG,M,I1 692 FORMAT(' WVEC MISMATCH',1P,2D12.4,4I4) STOP ENDIF ELSE WV(IC)=WVEC(I1) LWV(IC)=.TRUE. ENDIF 4000 CONTINUE WVFACT=PI/(WV(IC)*WV(IC)) ITRY=0 C CHECK THAT THERE IS ROOM FOR NB*NB ITEMS IN SR,SI,... 4001 IF (IST(INX(IC,IU),IC,IU)+NB*NB.GT.MXS) THEN C THIS IS AN ERROR CONDITION C TRY TO REMOVE LOWER JTOTS AND PUSH EVERYONE DOWN WRITE(6,*) ' *** S-MATRICES WILL EXCEED STORAGE, MXS',MXS IF (PRNTLV.LT.3) WRITE(6,*) ' *** CURRENT JTOT =',JTOT ITRY=ITRY+1 IF (ITRY.GT.1) THEN WRITE(6,*) ' *** UNSUCCESSFUL. NUMBER OF TRIES =',ITRY STOP ENDIF CALL PUSH(NCALC,INX,IST,JT,NOP,MXJT,MXCL,SR,SI,L1,L2,MXS, 1 JTOT,JDFMAX,ISTOP,PRNTLV) GO TO 4001 ENDIF NS=0 ISTART=IST(INX(IC,IU),IC,IU) C I2 COUNTS ACROSS COLS/ I1 COUNTS DOWN ROWS DO 4010 I2=1,NB IXS=(IXB(I2)-1)*NOPEN DO 4010 I1=1,NB NS=NS+1 SR(ISTART+NS)=SREAL(IXS+IXB(I1)) SI(ISTART+NS)=SIMAG(IXS+IXB(I1)) L1(ISTART+NS)=L(IXB(I1)) L2(ISTART+NS)=L(IXB(I2)) 4010 CONTINUE NOP(INX(IC,IU),IC,IU)=NS ISTOP=ISTOP+NS ISTOPX=MAX(ISTOPX,ISTOP) C C ---------------------------------------------------------------- C PROCESS THIS SET W/PREVIOUSLY SAVED (3-IU) SETS I1=INX(IC,IU) IND1=IST(I1,IC,IU)+1 JT1=JT(I1,IC,IU) JACCMX(IC)=MAX(JACCMX(IC),JT1) IF (PRNTLV.LT.10.AND.PRNTLV.GE.5) 1 WRITE(6,616) IT,JTOT,M,INRG,ECHK,IEXCH,WT,NOPEN IF (PRNTLV.GE.5) WRITE(6,617) NOP(I1,IC,IU),IND1 617 FORMAT(12X,'PROCESSING',I3, 1 ' S-MATRIX ELEMENTS; START IN STORAGE =',I5) C INDICES FOR "OTHER" TAPE IO=3-IU I2=INX(IC,IO) 7100 IF (I2.LE.0) THEN C MAY WANT TO TRY TO CHECK THAT WE'VE GOTTEN TO EITHER JT2=0 C OR BEYOND MAX(JDFMX0,JDFMX1) GO TO 7000 ENDIF IND2=IST(I2,IC,IO)+1 JT2=JT(I2,IC,IO) JDIF=ABS(JT1-JT2) IF (JDIF.LE.JDFMAX.AND.PRNTLV.GE.5) 1 WRITE(6,618) I2,JT2,NOP(I2,IC,IO),IND2 618 FORMAT(' OTHER SET',I4,' JT, NS, IST',I5,I4,I6) C IN CODE BELOW NB WE MUST HAVE 1ST SPECTRAL LEVEL ASSOC W/ JT1 C AND 2ND SPECTRAL LEVEL ASSOC W/ JT2 IF (JDIF.LE.JDFMX0) THEN IF (IU.EQ.1) THEN CALL ADD0(JT1,NOP(I1,IC,IU),SR(IND1),SI(IND1),L1(IND1), 1 L2(IND1),JT2,NOP(I2,IC,IO),SR(IND2),SI(IND2), 2 L1(IND2),L2(IND2),JIN,JFN,IQ,ADD0R,ADD0I) ELSE CALL ADD0(JT2,NOP(I2,IC,IO),SR(IND2),SI(IND2),L1(IND2), 1 L2(IND2),JT1,NOP(I1,IC,IU),SR(IND1),SI(IND1), 2 L1(IND1),L2(IND1),JIN,JFN,IQ,ADD0R,ADD0I) ENDIF ADD0R=ADD0R*WVFACT ADD0I=ADD0I*WVFACT WAR(IC)=WAR(IC)+ADD0R WAI(IC)=WAI(IC)+ADD0I IF (PRNTLV.GE.5) WRITE(6,619) I00,ADD0R,ADD0I,WAR(IC),WAI(IC) 619 FORMAT(' >>> OMEGA-',I1,' REAL/IMAG ADD',1P,2D13.4, 1 ' SUM',2D13.4) ENDIF IF (JDIF.LE.JDFMX1) THEN IF (IU.EQ.1) THEN CALL ADD1(JT1,NOP(I1,IC,IU),SR(IND1),SI(IND1),L1(IND1), 1 L2(IND1),JT2,NOP(I2,IC,IO),SR(IND2),SI(IND2), 2 L1(IND2),L2(IND2),JIN,JFN,IQ,ADD1R,ADD1I) ELSE CALL ADD1(JT2,NOP(I2,IC,IO),SR(IND2),SI(IND2),L1(IND2), 1 L2(IND2),JT1,NOP(I1,IC,IU),SR(IND1),SI(IND1), 2 L1(IND1),L2(IND1),JIN,JFN,IQ,ADD1R,ADD1I) ENDIF ADD1R=ADD1R*WVFACT ADD1I=ADD1I*WVFACT WRR(IC)=WRR(IC)+ADD1R WRI(IC)=WRI(IC)+ADD1I IF (PRNTLV.GE.5) WRITE(6,619) I01,ADD1R,ADD1I,WRR(IC),WRI(IC) ENDIF I2=I2-1 GO TO 7100 C ---------------------------------------------------------------- C 9000 WRITE(6,620) IT 620 FORMAT(' *** NORMAL END OF FILE ON UNIT',I4) LEOF(IU)=.TRUE. IF (LEOF(3-IU)) GO TO 9999 C 7000 CONTINUE C THIS COMPLETES A SET OF S-MATRICES FOR EACH TAPE; GO FOR ANOTHER GO TO 3000 C C END OF FILE CONDITIONS C NORMAL END OF FILE 9999 CALL GCLOCK(TEND) TIME=TEND-TSTART WRITE(6,630) TIME,ISTOPX,MXS 630 FORMAT(/' *** GHM *** NORMAL TERMINATION.'/' *** CPU TIME WAS', 1 F9.2,' SEC'/ 2 ' *** STORAGE FOR SAVED MATRIX ELEMENTS WAS',I7, 3 '; MAX ALLOWED =',I7) WRITE(6,631) (ECALC(I),JACCMX(I),WAR(I),WAI(I),WRR(I),WRI(I), 1 I=1,NCALC) 631 FORMAT(/' GENERALIZED HESS METHOD CROSS SECTIONS, ANGSTROMS**2,'/ 1 /' ENERGY(1/CM) MAX(J) REAL(S=0) IMAG(S=0)', 2 ' REAL(S=1) IMAG(S=1)'/(F15.4,I8,4F13.4)) STOP C 9100 WRITE(6,*) ' *** PREMATURE END OF FILE READING UNIT',IT WRITE(6,*) ' *** ACCUMULATED GHM CROSS SECTIONS TO THIS POINT ***' GO TO 9999 9900 WRITE(6,*) ' *** UNABLE TO READ FROM UNIT',IT STOP END SUBROUTINE PUSH(NCALC,INX,IST,JT,NOP,MXJT,MXCL,SR,SI,L1,L2,MXS, 1 JTOT,JDFMAX,ISTOP,IPRINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C ROUTINE TO PUSH DOWN STORAGE STACKS FOR GHM CODE C THIS IS "VIB" CODE -- TWO INPUT TAPES -- VERSION 3. C ASSUMES THAT JTOT'S ARE STORED TO INCREASING ORDER C DIMENSION INX(MXCL,2),IST(MXJT,MXCL,2),JT(MXJT,MXCL,2), 1 NOP(MXJT,MXCL,2),SR(MXS),SI(MXS),L1(MXS),L2(MXS) C WRITE(6,*) ' *** PUSH: ATTEMPT TO RECYCLE STORAGE' C NUMNOW IS CURRENT (MAX) NUMBER OF SAVED SETS NUMNOW=0 DO 1001 IE=1,NCALC DO 1001 I=1,2 1001 NUMNOW=MAX(NUMNOW,INX(IE,I)) IF (NUMNOW.LE.0) THEN WRITE(6,*) ' *** PUSH: NO CURRENT SAVED SETS ?' STOP ENDIF C DECIDE HOW MANY OF THESE CAN BE ELIMINATED, NELIM NELIM=0 DO 1101 IC=1,NUMNOW DO 1102 IE=1,NCALC DO 1102 I=1,2 IF (INX(IE,I).LT.IC) GO TO 1102 C BE CONSERVATIVE; COMPARE WITH JDFMAX+1 IF (ABS(JTOT-JT(IC,IE,I)).LE.JDFMAX+1) GO TO 1103 1102 CONTINUE C IF WE REACH HERE, THIS SET NUMBER CAN BE ELIMINATED NELIM=IC 1101 CONTINUE 1103 IF (NELIM.LE.0) THEN WRITE(6,*) ' *** PUSH: NO SAVED SETS CAN BE ELIMINATED.' STOP ELSE IF (IPRINT.GE.3) WRITE(6,601) NELIM,NUMNOW,JDFMAX 601 FORMAT(' *** PUSH: WILL ELIMINATE',I4,' OF',I4,' SAVED SETS,', 1 ' BASED ON JDFMAX =',I3) NUMNEW=NUMNOW-NELIM C KEEP AT LEAST ONE SET (THIS SHOULD NOT HAPPEN) IF (NUMNEW.LE.0) THEN IF (IPRINT.GE.3) WRITE(6,*) ' *** PUSH: KEEP AT LEAST ONE SET' NELIM=NELIM-1 NUMNEW=NUMNEW+1 ENDIF ENDIF C DO 1201 I=1,2 DO 1201 IE=1,NCALC INX(IE,I)=MAX(0,INX(IE,I)-NELIM) DO 1201 IC=1,NUMNEW IST(IC,IE,I)=IST(IC+NELIM,IE,I) JT(IC,IE,I)=JT(IC+NELIM,IE,I) NOP(IC,IE,I)=NOP(IC+NELIM,IE,I) IF (IPRINT.GE.5) WRITE(6,699) IC,IE,I,IST(IC,IE,I),JT(IC,IE,I) 699 FORMAT(5X,'NEW IST(',I4,I3,I2,') =',I6,' JTOT =',I4) 1201 CONTINUE C FIND THE SMALLEST SR,SI,... WHICH IS STILL NEEDED. C NB IST() IS 1 LESS THAN THE INDEX IN STORAGE ISXX=MXS+1 DO 1301 I=1,2 DO 1301 IE=1,NCALC DO 1301 IC=1,NUMNEW IF (IST(IC,IE,I).GT.0) ISXX=MIN(ISXX,IST(IC,IE,I)) 1301 CONTINUE IF (ISXX.LE.0) THEN WRITE(6,*) ' *** PUSH: CANNOT REDUCE S-MATRIX STORAGE' STOP ENDIF C MODIFY IST() POINTERS DO 1302 I=1,2 DO 1302 IE=1,NCALC DO 1302 IC=1,NUMNEW 1302 IST(IC,IE,I)=MAX(IST(IC,IE,I)-ISXX,0) C SHIFT IN STORAGE ISTOP=ISTOP-ISXX DO 1303 IC=1,ISTOP SR(IC)=SR(IC+ISXX) SI(IC)=SI(IC+ISXX) L1(IC)=L1(IC+ISXX) 1303 L2(IC)=L2(IC+ISXX) IF (IPRINT.GE.3) WRITE(6,602) ISXX 602 FORMAT(' *** PUSH: REMOVED',I6,' SAVED ELEMENTS') RETURN END