C PROGRAM TO COMPUTE S=0,1 CROSS SECTIONS OF GENERALIZED HESS METHOD C FROM S-MATRIX SAVE TAPE. CLOSE COUPLING, ITYPE=1, 2, 7 C (SHOULD ALSO WORK FOR ITYPE = 5,6) C VERSION 4 (SINGLE SAVE TAPE); 22 DEC 94, SG C THESE CROSS SECTIONS CORRESPOND TO THE THEORY GIVEN BY: C MONCHICK & HUNTER, JCP 85, 713 (1986). C SHAEFER & MONCHICK, ASTRON ASTROPHYS 265, 859 (1992) C MONCHICK, JCP 101, 5566 (1994). C CROSS SECTIONS HERE, HOWEVER, ARE MULTIPLIED BY PI/WVEC**2 C AND HENCE HAVE UNITS OF DISTANCE SQUARED. C THE S=0,1 CROSS SECTIONS (VARIABLES WAR,WAI AND WRR,WRI) C CORRESPOND WITH COLLISION INTEGRALS OMEGA(0,0,0), OMEGA(1,1,1). C THE FORMER CORRESPONDS WITH W_A OF GHM THEORY, BUT W_R OF GHM C THEORY IS A MASS-WEIGHTED COMBINATION OF THESE. C C REQUIRES MOLSCAT GCLOCK ROUTINE FOR CPU TIME. 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=50000, MXJT=200, MXLN=100) DIMENSION SR(MXS),SI(MXS),L1(MXS),L2(MXS) DIMENSION JT(MXJT,MXLN,2),IST(MXJT,MXLN,2),NOP(MXJT,MXLN,2), 1 WV(MXLN),EREL(MXLN), 2 JACCMX(MXLN),WAR(MXLN),WAI(MXLN),WRR(MXLN),WRI(MXLN) LOGICAL LWV(MXLN) C C STORAGE FOR &INPUT LINE INPUT; DIMENSIONS FOR OFF-DIAG (NOT USED) PARAMETER (MXLINE=10) DIMENSION LINE(4,MXLINE),LTYPE(MXLINE) C C STORAGE FOR THE "LINE TABLE" (BUT NO "PRIME" LEVELS YET). C ALONG THE LINES OF THE ONE IN MOLSCAT PRBR C 1 - L(IN), 2 - L(FN), 3 - LTYPE(I.E., IQ) C 4 - INRG(L(IN)), 5 - INRG(L(FN)), 6- INX(L(IN)), 7 - INX(L(FN)) DIMENSION LN(7,MXLN) C C STORAGE FOR SAVE TAPE VALUES PARAMETER (MXCH=200,MXJLEV=500) DIMENSION J(MXCH),L(MXCH),WVEC(MXCH),IXB(MXCH),JLEV(MXJLEV) 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),ELEVEL(MXEL) CHARACTER*4 LABEL(20) LOGICAL LFMT,EQUAL,PF,EXISTS INTEGER PRNTLV 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 NLPRBR - NUMBER OF 'LINES' IN LINE INPUT C LINE - SETS OF 4 "LEVEL" DESCRIPTORS (ONLY 2 USED FOR NOW) C LTYPE - "Q" VALUES FOR THE LINES (DEFAULTS TO DELTA-J) 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 /INPUT/IIN,LFMT,PRNTLV,NLPRBR,LINE,LTYPE C DATA I00/0/, I01/1/ C C FORMATS FOR SAVE TAPE (MOLSCAT 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 FUNCTIONS EQUAL(X,Y)=ABS(X-Y).LT.1.D-5 EXISTS(I)=I.GT.0.AND.I.LE.NLEVEL 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=10 NLPRBR=0 LFMT=.FALSE. ISTOP=0 ISTOPX=0 DO 1001 IE=1,MXLN 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 DO 1002 IE=1,MXLN DO 1002 I=1,MXJT IST(I,IE,IA)=0 JT(I,IE,IA)=0 1002 NOP(I,IE,IA)=0 DO 1003 IE=1,MXEL 1003 ELEVEL(IE)=0.D0 DO 1004 I=1,MXLINE LTYPE(I)=-1 DO 1004 IA=1,2 1004 LINE(IA,I)=0 C READ(5,INPUT,END=1000) C 1000 IF (NLPRBR.LE.0) THEN WRITE(6,*) ' *** MUST SPECIFY INPUT &INPUT NLPRBR .GT. 0' STOP ENDIF IF (NLPRBR.GT.MXLINE) THEN WRITE(6,*) ' *** NLPRBR REDUCED TO MXLINE',NLPRBR,MXLINE NLPRBR=MXLINE ENDIF WRITE(6,601) IIN,(I,(LINE(II,I),II=1,2),I=1,NLPRBR) 601 FORMAT(' FROM MOLSCAT S-MATRICES ON UNIT',I4/ 1 ' REQUESTED SPECTRAL LINES ARE'/ 2 (6X,'LINE',I3,' IS BETWEEN LEVELS',2I4)) 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 IT=IIN C GET AND PROCESS HEADER RECORD FROM SAVE TAPE IF (LFMT) THEN READ(IT,100,END=9900) LABEL,ITYPE,NLEV,NQN,URED,IP ELSE READ(IT,END=9900) LABEL,ITYPE,NLEV,NQN,URED,IP ENDIF WRITE(6,603) IT,LABEL,ITYPE,NLEV,NQN,URED,IP 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.LT.3) THEN WRITE(6,*) ' *** NO IMPLEMENTATION FOR IPROGM.LT.3' STOP ENDIF IF (LFMT.AND.IP.GE.11) THEN WRITE(6,*) ' *** IPROGM INCONSISTENT W/ REQUESTED TAPE FORMAT' STOP ENDIF C NB CODE SHOULD WORK FOR ITYPE=5,6 ALSO IF (.NOT.(ITYPE.EQ.1.OR.ITYPE.EQ.2.OR. 1 ITYPE.EQ.7)) THEN WRITE(6,*) ' *** CODE NOT IMPLEMENTED FOR ITYPE',ITYPE STOP ENDIF C C GET JLEV INFORMATION FROM TAPE NSQ=NLEV*NQN IF (NSQ.GT.MXJLEV) THEN WRITE(6,*) ' *** NLEV*NQN.GT.MXJLEV',MXJLEV STOP ENDIF IF (LFMT) THEN READ(IT,101) (JLEV(I),I=1,NSQ) ELSE READ(IT) (JLEV(I),I=1,NSQ) ENDIF C C GET NLEVEL/ELEVEL INFORMATION IF (LFMT) THEN IF (IP.GE.3) READ(IT,102) NLEVEL,(ELEVEL(I),I=1,NLEVEL) ELSE READ(IT) NLEVEL,(ELEVEL(I),I=1,NLEVEL) ENDIF WRITE(6,604) 604 FORMAT(' BASIS SET DATA:'/) DO 1005 I=1,NLEV 1005 WRITE(6,606) I,ELEVEL(ABS(JLEV(NLEV*(NQN-1)+I))), 1 (JLEV(NLEV*(K-1)+I),K=1,NQN) 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,(ENERGY(I),I=1,NNRG) ELSE READ(IT) NNRG,(ENERGY(I),I=1,NNRG) ENDIF IF (NNRG.GT.MXNRG) THEN WRITE(6,*) ' *** NNRG.GT.MXNRG',NNRG,MXNRG STOP ENDIF WRITE(6,607) NNRG,(ENERGY(I),I=1,NNRG) 607 FORMAT(/' NUMBER OF COLLISION ENERGIES ON SAVE TAPE IS',I4/ 1 (' ',1P,6D13.4)) C C FIND APPROPRIATE ENERGY MATCHES; ENERGY = EREL + ELEVEL WRITE(6,*) WRITE(6,608) 608 FORMAT(' ***') NLINE=0 DO 1200 I=1,NLPRBR N1=LINE(1,I) N2=LINE(2,I) PF=.FALSE. IF (.NOT.(EXISTS(N1).AND.EXISTS(N2))) THEN WRITE(6,108) N1,N2 108 FORMAT(' *** WARNING. REQUESTED LINE FOR LEVELS',2I4, 1 ' CANNOT BE PROCESSED - OUTSIDE NLEVEL RANGE.') GO TO 1200 ENDIF C FOR ITYPE=1,2,7 TAKE J-VALUE FROM JLEV(LEV,1) JA=JLEV(N1) JB=JLEV(N2) WRITE(6,109) N1,N2,JA,JB 109 FORMAT(' *** REQUESTED LINE: LEVELS',2I4,', J-VALUES',2I4) IF (LTYPE(I).GE.0) THEN IQ=LTYPE(I) WRITE(6,107) IQ 107 FORMAT(6X,'Q TAKEN FROM &INPUT LTYPE =',I2) ELSE C IF 'LTYPE' NOT INPUT, CALCULATE FROM J-VALUES IQ=IABS(JA-JB) WRITE(6,106) IQ 106 FORMAT(6X,'Q COMPUTED FROM JA,JB =',I2) ENDIF C OTHER THAN DIPOLE AND RAMAN (Q=0,2) TRANSITIONS NOT IMPLEMENTED. IF (.NOT.((IQ.EQ.0.OR.IQ.EQ.1.OR.IQ.EQ.2))) THEN WRITE(6,613) 613 FORMAT(' ***** ONLY Q = 0,1,2 IMPLEMENTED -- SKIPPING') GO TO 1200 ENDIF NEPL=0 DO 1400 II=1,NNRG EA=ENERGY(II)-ELEVEL(N1) IF (EA.LE.0.D0) GO TO 1400 DO 1500 JJ=1,NNRG EB=ENERGY(JJ)-ELEVEL(N2) ED=(EA-EB)/MAX(EA,EB,1.D0) IF (.NOT.EQUAL(ED,0.D0)) GO TO 1500 IF (NLINE.GE.MXLN) THEN WRITE(6,614) N1,N2,II,JJ 614 FORMAT(' ***** WARNING. NO MORE TABLE SPACE FOR LEVELS =',2I3, 1 ', AND ENERGIES=',2I4) GO TO 1200 ENDIF C NB STORAGE IN LN FOLLOWS MOLSCAT/PRBR ROUTINE NLINE=NLINE+1 LN(1,NLINE)=N1 LN(2,NLINE)=N2 LN(3,NLINE)=IQ LN(4,NLINE)=II LN(5,NLINE)=JJ LN(6,NLINE)=0 LN(7,NLINE)=0 EREL(NLINE)=(EA+EB)/2.D0 NEPL=NEPL+1 WRITE(6,612) EREL(NLINE),LN(4,NLINE),LN(5,NLINE) 612 FORMAT(10X,'CALCULATED FOR RELATIVE K.E. =',F12.4, 1 ' (1/CM) USING ENERGIES',2I4) 1500 CONTINUE 1400 CONTINUE IF (NEPL.LE.0) WRITE(6,611) 611 FORMAT(8X,'***** NOTE. NO RELEVANT ENERGIES FOR THIS LINE') 1200 CONTINUE C IF (NLINE.LE.O) THEN WRITE(6,*) ' *** SORRY. NO ENERGY MATCHES WERE FOUND' WRITE(6,*) ' CHECK INPUT FILES AND/OR &INPUT' STOP ENDIF WRITE(6,608) WRITE(6,*) C THIS ENDS PROCESSING OF 'HEADER' INFORMATION FOR BOTH TAPES. C BELOW FOLLOWS 22 DEC 94 SG ANALYSIS THAT JDFMX IS GOVERNED BY C IQ+LAMBDA: NEED HERE LAMBDA=1 LIMIT. JDFMAX=0 DO 1700 I=1,NLINE 1700 JDFMAX=MAX(JDFMAX,LN(3,I)+1) 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 C FIRST, GET JTOT,INRG,M RECORD ----- 3000 IF (LFMT) THEN READ(IT,103,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M ELSE IF (IP.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 3000 ENDIF C ECHO INPUT IF (PRNTLV.GE.10) THEN IF (IP.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))) 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.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 2000 IC=1,NLINE JIN=JLEV(LN(1,IC)) JFN=JLEV(LN(2,IC)) DO 2000 IU=1,2 IF (LN(IU+3,IC).NE.INRG) GO TO 2000 if (prntlv.ge.15) write(6,621) ic,iu 621 format(5x,'processing S-matrices for line',i4,', iu =',i2) C C IF WE REACH HERE, WE WANT TO PROCESS THE (JTOT,INRG,M) BLOCK C FOR IC-TH LINE; LN(IU+5) HAS NO. SETS SAVED SO FAR. LN(IU+5,IC)=LN(IU+5,IC)+1 IF (LN(IU+5,IC).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 LN(IU+5,IC)=LN(IU+5,IC)-1 CALL PUSH(NLINE,LN,IST,JT,NOP,MXJT,MXLN,SR,SI,L1,L2,MXS, 1 JTOT,JDFMAX,ISTOP,PRNTLV) LN(IU+5,IC)=LN(IU+5,IC)+1 ENDIF C N.B. INX IS INDEX NO. OF CURRENT SET -- ALWAYS .LT. MXJT INX=LN(IU+5,IC) IST(INX,IC,IU)=ISTOP JT(INX,IC,IU)=JTOT IF (JTOT.GT.JTOLD) THEN IF (JTOLD.NE.-1.AND.PRNTLV.GE.3) 1 WRITE(6,625) JTOLD,(LN(1,IE),LN(2,IE),EREL(IE), 2 WAR(IE),WAI(IE),WRR(IE),WRI(IE),IE=1,NLINE) 625 FORMAT(/' *** FINISHED JTOT =',I4/ 1 (6X,'LINE',I3,' -',I3,', E =',F10.2,' S=0 IS', 2 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.LN(IU,IC)) 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 IF (NB.LE.0) THEN C SINCE NO RELEVANT BASIS FNS FOUND, UNDO SETTINGS FOR THIS SET. JT(INX,IC,IU)=0 IST(INX,IC,IU)=0 INX=INX-1 LN(IU+5,IC)=INX GO TO 2000 ENDIF 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)+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(NLINE,LN,IST,JT,NOP,MXJT,MXLN,SR,SI,L1,L2,MXS, 1 JTOT,JDFMAX,ISTOP,PRNTLV) GO TO 4001 ENDIF C STORE APPROPRIATE S-MATRIX ELEMENTS NS=0 ISTART=IST(INX,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)=NS ISTOP=ISTOP+NS ISTOPX=MAX(ISTOPX,ISTOP) C C ---------------------------------------------------------------- C PROCESS THIS SET W/PREVIOUSLY SAVED (3-IU) SETS C BELOW FOLLOWS 22 DEC 94 SG ANALYSIS THAT JDFMX IS GOVERNED BY C IQ+LAMBDA. JDFMX0=LN(3,IC) JDFMX1=LN(3,IC)+1 I1=INX 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" (I/F) LEVEL ... IO=3-IU I2=LN(IO+5,IC) 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 2000 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 2000 CONTINUE C ---------------------------------------------------------------- C THIS COMPLETES CURRENT SET OF S-MATRICES; GO FOR ANOTHER GO TO 3000 C C END OF FILE CONDITIONS C NORMAL END OF FILE 9000 WRITE(6,620) IT 620 FORMAT(' *** NORMAL END OF FILE ON UNIT',I4) 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) (LN(1,I),LN(2,I),EREL(I),JACCMX(I),WAR(I),WAI(I), 1 WRR(I),WRI(I),I=1,NLINE) 631 FORMAT(/' GENERALIZED HESS METHOD CROSS SECTIONS, ANGSTROMS**2,'/ 1 /' I - F ENERGY(1/CM) MAX(J) REAL(S=0) IMAG(S=0)', 2 ' REAL(S=1) IMAG(S=1)'/(I3,' -',I2,F13.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(NLINE,LN,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 (SINGLE INPUT TAPE -- VERSION 4). C ASSUMES THAT JTOT'S ARE STORED TO INCREASING ORDER C DIMENSION LN(7,MXCL),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,NLINE DO 1001 I=1,2 1001 NUMNOW=MAX(NUMNOW,LN(I+5,IE)) 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,NLINE DO 1102 I=1,2 IF (LN(I+5,IE).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,NLINE LN(I+5,IE)=MAX(0,LN(I+5,IE)-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,NLINE 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,NLINE 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 S-MATRIX ELEMENTS') RETURN END