C THIS ROUTINE IS THE MAIN PROGRAM FOR MOLSCAT VERSION 14 C WITH DYNAMIC SPACE ALLOCATION CAPABILITY. C C INCREASE MXDIM AS NECESSARY TO PROVIDE SUFFICIENT WORKSPACE. C IXNEXT,NIPR ARE INITIALIZED IN DRIVER C IVLFL IS ALSO SET IN DRIVER; COULD BE CHANGED BY BASIN ROUTINES C PARAMETER (MXDIM=250000) DOUBLE PRECISION X DIMENSION X(MXDIM) COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X C MX=MXDIM C C CALL PRINCIPAL MOLSCAT/BOUND SUBROUTINE C CALL DRIVER STOP END SUBROUTINE DRIVER C*********************************************************************** C C ------ MOLSCAT - J.M. HUTSON AND S.GREEN - VERSION 14 - JUL 94 ----- C C MAIN DRIVER FOR QUANTUM MOLECULAR SCATTERING PROGRAM C C REVISION HISTORY SINCE VERSION 7 OF SHELDON GREEN'S QCPE PROGRAM C (MAY 79): C C VARIOUS NEW PROPAGATORS HAVE BEEN ADDED SINCE EARLY VERSIONS. C THE COMPLETE LIST IN VERSION 12 IS: C C INTFLG =-1 : WKB METHOD FOR SINGLE CHANNEL, SINGLE TURNING POINT C INTFLG = 2 : DEVOGELAERE'S PROPAGATOR C INTFLG = 3 : WALKER-LIGHT R-MATRIX PROPAGATOR C INTFLG = 4 : HYBRID LOG-DERIVATIVE / VIVS (VIVAS) PROPAGATOR C INTFLG = 5 : JOHNSON'S LOG-DERIVATIVE PROPAGATOR C INTFLG = 6 : MANOLOPOULOS'S DIABATIC MODIFIED C LOG-DERIVATIVE PROPAGATOR C INTFLG = 7 : MANOLOPOULOS'S QUASIADIABATIC MODIFIED C LOG-DERIVATIVE PROPAGATOR C INTFLG = 8 : ALEXANDER-MANOLOPOLOUS MODIFIED LOG-DERIVATIVE C AIRY PROPAGATOR (HIBRIDON) C VERSION 8: CHANGES MADE BY CHRIS ASHTON (1982) AND JEREMY HUTSON C (1982-4) AT WATERLOO AND CAMBRIDGE UNIVERSITIES. C C (1) ENTIRE PROGRAM CONVERTED TO DOUBLE PRECISION C C (2) GORDON ALGORITHM (INTFLG=1) REMOVED. C C (3) LOOP OVER "PARITY CASES" IN DRIVER HAS BEEN MADE EXPLICIT C FOR CLARITY. C C (4) EIGENPHASE SUM CALCULATION AND RESONANCE SEARCH OPTION C INCORPORATED. NEW OUTPUT CHANNEL (KSAVE) WITH OPTIONAL C UNFORMATTED OUTPUT ON CHANNEL ISAVEU. C C (5) COLLISION TYPE ITYPE=10*N+7 HAS BEEN ADDED, C FOR AN ATOM HITTING A DIATOMIC VIB-ROTOR, WHERE THE C POTENTIAL MATRIX IS CONSTRUCTED BY DOING PROPERLY THE C AVERAGING OF POTENTIAL TERMS OVER (V,J) AND (V',J') DIATOM C INTERNAL STATES. C C (6) COLLISION TYPE ITYPE=8 ADDED, FOR ELASTIC SCATTERING OF ATOMS C FROM CORRUGATED SURFACES. USES SUBROUTINE SURBAS TO SET UP C THE BASIS SET. THE LOOPS IN DRIVER OVER JTOT AND M ARE USED C TO LOOP OVER ANGLES THETA AND PHI RESPECTIVELY. C C (7) THE STORAGE OF THE COUPLING ARRAY VL HAS BEEN REARRANGED. THE C METHOD OF CONSTRUCTING POTENTIAL MATRICES FROM IT HAS BEEN C CHANGED, AND IN PARTICULAR A NEW INDEXING ARRAY IV HAS BEEN C INTRODUCED. C C*********************************************************************** C C VERSION 9 (APR 86): JMH AND SG CODES UNIFIED C C (9) IOS CODE RE-INCORPORATED FROM SG'S PROGRAM. C IT IS ACCESSED BY SETTING ITYPE = 100 + 'ITYPE' C C (10) MANOLOPOULOS'S DIABATIC AND ADIABATIC MODIFIED LOG-DERIVATIVE C PROPAGATORS ADDED (INTFLG=6 AND 7 RESPECTIVELY). C C*********************************************************************** C C SG VERSION 10 (AUG 91): C C (10) NEW PRBR/IOSPB FOR OFF-DIAGONAL LINESHAPE CROSS SECTIONS, C WITH HAS IN-CORE SIMULATION OF DIRECT ACCESS FILES. C OUTPUT CROSS-SECTIONS NOW MULTIPLIED BY JSTEP (FOR JTOT). C C (11) ALEXANDER/MANOLOPOULOS MODIFIED LOG-DERIVATIVE/AIRY PROPAGATOR C ADDED AS INTFLG=8. INTERFACED BY TIM PHILLIPS (NASA/GISS) C C VERSION 11 (JUN 92): JMH AND SG CODES INTEGRATED AGAIN. C C (12) LOOP OVER ENERGY IN DRIVER MODIFIED TO SIMPLIFY PARALLELISATION C C (13) ISAVEU OUTPUT MODIFIED TO USE UNFORMATTED WRITES C C (14) USAGE OF LINEAR ALGEBRA AND BLAS ROUTINES UNIFIED C C AND THE FOLLOWING ENHANCEMENTS ADDED FROM JMH'S CODE: C C (15) BASE9 INTERFACE ADDED C C (16) POTENL ENHANCED TO EVALUATE RADIAL STRENGTH FUNCTIONS BY C QUADRATURE FOR ITYPE=1, 2, 5 AND 6. C C (17) CODE ADDED TO CALCULATE ASYMMETRIC TOP ENERGIES AND WAVEFUNCTION C FROM ROTATIONAL CONSTANTS. MECHANISM FOR SELECTING ASYMMETRIC C TOP STATES TO BE INCLUDED GENERALISED C C (18) CODE FOR ATOM-SPHERICAL TOP SCATTERING ADDED C C*********************************************************************** C C VERSION 12 (MAY 93) C C (19) DYNAMIC STORAGE HANDLING COMPLETELY REORGANIZED. C C (20) VECTOR/MATRIX ROUTINES RATIONALIZED TO USE LAPACK AND BLAS. C C (21) IV() ARRAY USED ONLY FOR 'NON-TRIVIAL' CASES. C C (22) OPTION TO WRITE VL ARRAY TO DISC TO AVOID EXCESSIVE MEMORY USE. C C (23) SOME CODE FOR COUPLING VL MATRIX ELEMENTS MODIFIED TO AVOID C UNNECESSARY RECALCULATION OF NJ COEFFICIENTS C C*********************************************************************** C C VERSTION 13 (SG EXPERIMENTAL VERSION) APR 94, BUT CONTAINED IN V14 C C (24) IV ARRAY INTRODUCED FOR ITYPE=2 CASES C C (25) EXPANDED POTENL CAPABILITIES C C (26) BIGGER DIMENSIONS: /CMBASE/ ...ELEVEL(1000),...,JLEVEL(4000),... C ALSO ADD ISYM(10),ISYM2(10); REORGANIZED ORDERING C C (27) CHANGES TO CALLING SEQUENCE FOR OUTINT/OUTPCH; IEXCH NOW CORRECT C ON ISAVEU TAPE; BASE/OUTPCH RECOGNIZE CS SIGMA WHICH ARE C NOT COMPLETE. C C*********************************************************************** C C VERSION 14 (JUL 94) C C (28) ISAVEU TAPE FORMAT CHANGE: NOPEN WITH JTOT,INRG,...,M,NOPEN REC C C (29) FILE='FILENAME' REMOVED FROM OPEN STATEMENTS C C (30) FLAG FOR NCAC,DTOL,OTOL INCREASED TO JTOTU=999999 C C (31) RESTART ABILITIES (IRSTRT) FROM ISAVEU C C (32) ITYPE=4 CODE (ASYMMETRIC TOP - LINEAR ROTOR) ADDED C C (33) COMMON /CMBASE/ ALTERED TO ALLOW MORE SPACE FOR LEVELS AND C INTRODUCE EXTRA INPUT VARIABLES FOR HANDLING FUTURE EXTENSIONS. C THIS CHANGE REQUIRES SIMILAR CHANGES IN BASE9 ROUTINES. C C (34) HANDLING OF TOTAL ENERGIES CHANGED IN PRESSURE BROADENING WITH C IFEGEN=2 OPTION: AVOID CALCULATING S MATRICES THAT ARE NOT USED. C C*********************************************************************** C C EXTERNAL UNITS FOR MASSES ARE ATOMIC MASS UNITS (CARBON MASS/12) C EXTERNAL UNITS FOR ENERGIES ARE WAVENUMBERS C EXTERNAL UNITS FOR LENGTH RM ARE ANGSTROMS C ALL OTHER LENGTHS ARE IN UNITS OF RM C C INTFLG CONTROLS METHOD OF SOLVING EQUATIONS. NPOTL AND MXLAM C FOR SUM OVER ANGULAR DEPENDENCE OF POTENTIAL, NQN IS NO. OF C QUANTUM NUMBERS NECESSARY TO DESCRIBE COLLISION PARTNERS. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C ***** PROGRAM DIMENSION LIMITATIONS ***** C ENERGY,TEMP,LINE DIMENSIONS LIMITED BY VALUES ... PARAMETER(MXNRG=100,MXLN=200,MXTEMP=5) C INTEGER EUNITS,PRNTLV,PRINT,SHRINK C C ARRAY TO HOLD TIME AND DATE C INTEGER CTIME(2),CDATE(4) CHARACTER CTIME*9,CDATE*11 C C TYPES FOR COMMON/LDVVCM/ LOGICAL IALFP,IV,IVP,IVPP,NUMDER,ISHIFT,IDIAG,IPERT,ISYM LOGICAL LCALC,ALDONE LOGICAL IREAD,IWRITE LOGICAL LWARN C C DOUBLE PRECISION LABEL(10) CHARACTER*80 LABEL CHARACTER*80 LABL CHARACTER*1 TITLE(80),TIT(120),TIT2(120),BL CHARACTER*8 PDATE CHARACTER*8 CWD(2) EQUIVALENCE (LABL,TITLE(1)) C C FOLLOWING ARRAYS ALL HAVE DIMENSION MXNRG. MXNRG IS THE MAXIMUM C ALLOWED NUMBER OF TOTAL ENERGIES PER RUN. DIMENSION ENERGY(MXNRG) DIMENSION IECONV(MXNRG),ISST(MXNRG),MINJT(MXNRG),MAXJT(MXNRG) C C VARIABLES DIMENSIONED FOR NO. OF LINES IN PRES. BROAD. CALC. C N.B. PRBRIN STILL MAX NO. LINES = 2*MXLN DESPITE OFF-DIAG CHANGES DIMENSION LINE(2*MXLN),LTYPE(MXLN) EQUIVALENCE (ILSU,IPRBRU), (NLPRBR,IFLS) C DIMENSION TEMP(MXTEMP) C C VARIABLES TO TEST PARTIAL WAVE CONVERGENCE DIMENSION TEST(2) EQUIVALENCE (TEST(1),DTOL),(TEST(2),OTOL) C DIMENSION NLABV(9) C C DYNAMIC STORAGE COMMON BLOCK ... COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1) C MX,IXNEXT ARE MAX AND NEXT AVAILABLE LOCATION IN X() ARRAY C IVLFL FLAGS WHETHER IV() ARRAY IS USED AS POINPER W/ VL ARRAY. C NIPR IS NUMBER OF INTEGERS PER REAL; SHOULD BE 1 OR 2. C E.G. FOR IBM R*8/I*4, NIPR=2. AN INTEGER ARRAY OF DIM. N C CAN BE STORED IN A REAL ARRAY OF DIMENSION (N+NIPR-1)/NIPR. C C COMMON BLOCK FOR COMMUNICATING WITH COUPLED EQUATION SOLVERS COMMON/DRIVE/STEST,STEPS,STABIL,CONV,RSTART,RSTOP,XEPS, 1 DR,DRMAX,RMID,TOLHI,RTURN,VTOL,ESHIFT,ERED,RMLMDA, 2 NOPEN,JKEEP,ISCRU,MAXSTP C C EXTRA COMMON BLOCK FOR LDVIVS COMMON/LDVVCM/XSQMAX,ALPHA1,ALPHA2,IALPHA,IALFP,IV,IVP,IVPP, 1 NUMDER,ISHIFT,IDIAG,IPERT,ISYM,IREAD,IWRITE C C COMMON BLOCK FOR WKB INTEGRATOR COMMON/WKBCOM/NGMP(3) C C COMMON BLOCK TO SUBROUTINE OUTPUT FOR USE IN RESONANCE SEARCHES COMMON/EIGSUM/EPSM(5) C C COMMON BLOCK FOR AIRPRP ARGUMENTS IN MANOLOPOLOUS/ALEXANDER C PROPAGATOR COMMON/HIBRIN/POWRX,DRAIRY,IABSDR C COMMON/VLSAVE/IVLU C NAMELIST /INPUT/ LABEL,RMIN,RMAX,IRMSET,IRXSET,URED,ISCRU,ISIGPR 1 ,ITHROW,STEST,NNRG,ENERGY,DNRG,JTOTL,JTOTU,JSTEP,MSET,MHI,NCAC 2 ,PRNTLV,INTFLG,MXSIG,STEPS,STABIL,NTEMP,NGAUSS,TEMP,EUNITS 3 ,ISIGU,IPARTU,ILSU,IPRBRU,IFLS,NLPRBR,LINE,IFEGEN,LTYPE,MAXSTP 4 ,TOLHI,RVIVAS,RVFAC,XSQMAX,ALPHA1,ALPHA2,IALPHA 5 ,IALFP,IV,IVP,IVPP,NUMDER,ISHIFT,IDIAG,IPERT,ISYM 6 ,ISAVEU,DTOL,OTOL,KSAVE,DR,DRNOW,DRMAX,RMID,VTOL,ICONV 7 ,THETLW,THETST,PHILW,PHIST,MXPHI,SHRINK,LASTIN 8 ,MMAX,LMAX,NGMP 9 ,VMAX,TMAX,TOLLO,CTOL,UTEST,TOLER,TOL,MXXX,MNNN A ,POWRX,DRAIRY,IABSDR,NNRGPG,IRSTRT C EQUIVALENCE (MXPAR,MXPHI), (RMID,RVIVAS), (DR,DRNOW), 1 (TOL,TOLER,TOLHI) C C NGPT,LMAX, MMAX, AND NGMP(3) ARE VARIABLES ADDED FOR C COMPATIBILITY WITH THE IOS PROGRAMS C VARIABLES VMAX,...,MNNN ADDED FOR COMPATIBILITY WITH S.GREEN CODE C (MOSTLY GORDON INTEGRATOR). ALSO TOL, TOLER, DRNOW C C RMIN IS THE RADIUS AT WHICH THE INTEGRATION IS BEGUN C RMAX IS THE OUTER RADIUS TO WHICH THE INTEGRATION MUST EXTEND C MAXSTP IS MAX NO. OF STEPS IN RADIAL INTEGRATION (INTFLG=3 ONLY) C C ARRAYS FOR NAMELIST SIMULATOR C CHARACTER*6 INAMES C DIMENSION INAMES(88),LOCN(88),INDX(88) C C DATA INAMES/'LABEL','RMIN','RMAX','IRMSET','IRXSET', C 1 'URED','ISCRU','ISIGPR', C 1 'ITHROW','STEST','NNRG','ENERGY','DNRG', C 2 'JTOTL','JTOTU','JSTEP','MSET','MHI','NCAC', C 2 'PRNTLV','INTFLG','MXSIG','STEPS','STABIL', C 3 'NTEMP','NGAUSS','TEMP','EUNITS','ISIGU','IPARTU','ILSU', C 4 'IPRBRU','IFLS','NLPRBR','LINE','IFEGEN','LTYPE','MAXSTP', C 4 'TOLHI','RVIVAS','RVFAC','XSQMAX','ALPHA1','ALPHA2','IALPHA', C 5 'IALFP','IV','IVP','IVPP','NUMDER','ISHIFT','IDIAG','IPERT', C 6 'ISYM','ISAVEU','DTOL','OTOL','KSAVE','DR','DRNOW','DRMAX', C 7 'RMID','VTOL','ICONV','THETLW','THETST','PHILW','PHIST', C 8 'MXPHI','SHRINK','LASTIN','MMAX','LMAX','NGMP','VMAX', C 9 'TMAX','TOLLO','CTOL','UTEST','TOLER','TOL','MXXX','MNNN' C A 'PWRX','DRAIRY','IABSDR','NNRGPG','IRSTRT'/ C DATA INDX/88*0/ C C DATA LABEL/10*' '/ DATA CWD/' ','(8-BYTE)'/ DATA CTIME/' '/,CDATE/' '/ DATA IPROGM/14/, PDATE/'(AUG 94)'/ DATA TITLE/80*' '/, BL/' '/ DATA TIT/120*'='/, TIT2/120*'-'/ C DATA LTYPE/MXLN*-1/ C C NLABV ARRAY CONTAINS NUMBER OF LABELS PER SYMMETRY TERM FOR EACH C VALUE OF ITYPE (ITYPE=4 ADDED JUL 94 TRP/SG) DATA NLABV/1,3,3,4,2,2,5,2,1/ C C THE PHYSICAL CONSTANTS USED ARE COMBINED IN THE SINGLE NUMBER BFCT. C BFCT IS 0.5*(HBAR**2) IN UNITS OF (ATOMIC MASS UNITS)*(WAVENUMBERS) C *(ANGSTROMS**2). C THE FOLLOWING VALUE IS FROM THE 1973 PHYSICAL CONSTANTS. DATA BFCT/16.857630D0/ C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CALL ROUTINE TO MASK FLOATING-POINT UNDERFLOW. CALL MASK C C STORE VALUE OF MX IN CASE IT NEEDS TO BE RESET; C NEEDED IN FUTURE CODE WHICH USES MAXMAX/MX TO ALLOCATE C 'PERMANENT' STORAGE FOR A RUN W/ MULTIPLE (LASTIN=0) INPUT DECKS MXSAVE=MX C 100 MX=MXSAVE CALL GCLOCK(TFIRST) CALL GDATE(CDATE) CALL GTIME(CTIME) WRITE(6,110) IPROGM,PDATE,CDATE,CTIME,IPROGM,PDATE 110 FORMAT(2X,8('----MOLSCAT----')/' |',120X,'|'/' |',24X, 1 'COUPLED CHANNEL MOLECULAR SCATTERING PROGRAM OF J. M. HUTSON ', 2 'AND S. GREEN',23X,'|'/' |',29X,'VERSION 1 BY S. GREEN ', 3 '(NOV 1973); THIS IS VERSION',I3,1X,A8,29X,'|'/ 4 ' |',120X,'|'/' |',44X,'RUN ON ',A11,2X, 5 'AT ',A9,44X,'|'/' |',120X,'|'/2X,8('----MOLSCAT----')// 6 2X,'PUBLICATIONS RESULTING FROM THE USE OF THIS PROGRAM SHOULD ', 7 'REFER TO'/2X,'J. M. HUTSON AND S. GREEN, MOLSCAT COMPUTER ', 8 'CODE, VERSION',I3,1X,A8 / 9 2X,'DISTRIBUTED BY COLLABORATIVE COMPUTATIONAL PROJECT NO. 6 ', A 'OF THE SCIENCE AND ENGINEERING RESEARCH COUNCIL (UK)') C C INITIALIZE STORAGE PARAMETERS IN /MEMORY/ NIPR=2 IXNEXT=1 C SET IVLFL TO 1 TO ENSURE STORAGE COMPATIBILITY W/ VERSION 11 IVLFL=1 C SET NUSED.LT.0 AND CALL CHKSTR TO RESET COUNTER FOR EACH &INPUT. NUSED=-1 CALL CHKSTR(NUSED) C C SET INITIAL VALUES BEFORE READ(5,INPUT) . . . C LWARN=.FALSE. IOSFLG=0 NGMP(1)=8 NGMP(2)=1 NGMP(3)=16 NNRG=0 NNRGPG=1 DNRG=0.D0 NTEMP=0 NGAUSS=3 JSTEP=1 JTOTL=-1 JTOTU=-1 MSET=0 MHI=0 MXSIG=0 ISIGPR=0 ITHROW=0 DTOL=0.3D0 OTOL=.005D0 NCAC=4 ISIGU = 0 IPARTU=0 ISAVEU=0 KSAVE=0 ILSU=11 IFLS=0 IFEGEN=0 ICONV=0 INTFLG=4 RMIN=0.8D0 RMAX=10.D0 STEST=1.D-4 STEPS=10.D0 STABIL=5.D0 ISCRU=0 IRMSET=9 IRXSET=1 DR=2.D-2 RMID=9999.D0 RVFAC=0.D0 DRMAX=5.D0 VTOL=1.D-06 MAXSTP=10000 TOLHI=0.001D0 XSQMAX=1.D04 ALPHA1=1.D0 ALPHA2=1.5D0 IALPHA=6 IALFP=.FALSE. IV=.TRUE. IVP=.FALSE. IVPP=.FALSE. NUMDER=.FALSE. ISHIFT=.FALSE. IDIAG=.FALSE. IPERT=.TRUE. ISYM=.TRUE. EUNITS=0 PRNTLV=0 MXPHI=1 THETLW=0.D0 THETST=0.D0 PHILW=0.D0 PHIST=0.D0 SHRINK=1 LASTIN=1 PI=ACOS(-1.D0) POWRX=3.D0 DRAIRY=-1.D0 IABSDR=0 IRSTRT=0 C C READ &INPUT DATA. C OPEN(5,STATUS='OLD',SHARED,READONLY) C---------------------------------------------------------------- C ARRAYS FOR NAMELIST SIMULATOR C LOCN(1)=LOC(LABEL) C LOCN(2)=LOC(RMIN) C LOCN(3)=LOC(RMAX) C LOCN(4)=LOC(IRMSET) C LOCN(5)=LOC(IRXSET) C LOCN(6)=LOC(URED) C LOCN(7)=LOC(ISCRU) C LOCN(8)=LOC(ISIGPR) C LOCN(9)=LOC(ITHROW) C LOCN(10)=LOC(STEST) C LOCN(11)=LOC(NNRG) C LOCN(12)=LOC(ENERGY) C LOCN(13)=LOC(DNRG) C LOCN(14)=LOC(JTOTL) C LOCN(15)=LOC(JTOTU) C LOCN(16)=LOC(JSTEP) C LOCN(17)=LOC(MSET) C LOCN(18)=LOC(MHI) C LOCN(19)=LOC(NCAC) C LOCN(20)=LOC(PRNTLV) C INDX(20)=4 C LOCN(21)=LOC(INTFLG) C LOCN(22)=LOC(MXSIG) C LOCN(23)=LOC(STEPS) C LOCN(24)=LOC(STABIL) C LOCN(25)=LOC(NTEMP) C LOCN(26)=LOC(NGAUSS) C LOCN(27)=LOC(TEMP) C LOCN(28)=LOC(EUNITS) C INDX(28)=4 C LOCN(29)=LOC(ISIGU) C LOCN(30)=LOC(IPARTU) C LOCN(31)=LOC(ILSU) C LOCN(32)=LOC(IPRBRU) C LOCN(33)=LOC(IFLS) C LOCN(34)=LOC(NLPRBR) C LOCN(35)=LOC(LINE) C LOCN(36)=LOC(IFEGEN) C LOCN(37)=LOC(LTYPE) C LOCN(38)=LOC(MAXSTP) C LOCN(39)=LOC(TOLHI) C LOCN(40)=LOC(RVIVAS) C LOCN(41)=LOC(RVFAC) C LOCN(42)=LOC(XSQMAX) C LOCN(43)=LOC(ALPHA1) C LOCN(44)=LOC(ALPHA2) C LOCN(45)=LOC(IALPHA) C LOCN(46)=LOC(IALFP) C LOCN(47)=LOC(IV) C LOCN(48)=LOC(IVP) C LOCN(49)=LOC(IVPP) C LOCN(50)=LOC(NUMDER) C LOCN(51)=LOC(ISHIFT) C LOCN(52)=LOC(IDIAG) C LOCN(53)=LOC(IPERT) C LOCN(54)=LOC(ISYM) C DO 115 I=46,54 C 115 INDX(I)=3 C LOCN(55)=LOC(ISAVEU) C LOCN(56)=LOC(DTOL) C LOCN(57)=LOC(OTOL) C LOCN(58)=LOC(KSAVE) C LOCN(59)=LOC(DR) C LOCN(60)=LOC(DRNOW) C LOCN(61)=LOC(DRMAX) C LOCN(62)=LOC(RMID) C LOCN(63)=LOC(VTOL) C LOCN(64)=LOC(ICONV) C LOCN(65)=LOC(THETLW) C LOCN(66)=LOC(THETST) C LOCN(67)=LOC(PHILW) C LOCN(68)=LOC(PHIST) C LOCN(69)=LOC(MXPHI) C LOCN(70)=LOC(SHRINK) C INDX(70)=4 C LOCN(71)=LOC(LASTIN) C LOCN(72)=LOC(MMAX) C LOCN(73)=LOC(LMAX) C LOCN(74)=LOC(NGMP) C LOCN(75)=LOC(VMAX) C LOCN(76)=LOC(TMAX) C LOCN(77)=LOC(TOLLO) C LOCN(78)=LOC(CTOL) C LOCN(79)=LOC(UTEST) C LOCN(80)=LOC(TOLER) C LOCN(81)=LOC(TOL) C LOCN(82)=LOC(MXXX) C LOCN(83)=LOC(MNNN) C LOCN(84)=LOC(POWRX) C LOCN(85)=LOC(DRAIRY) C LOCN(86)=LOC(IABSDR) C LOCN(87)=LOC(NNRGPG) C LOCN(88)=LOC(IRSTRT) C C CALL NAMLIS('&INPUT',INAMES,LOCN,INDX,88,IEOF) C IF(IEOF.EQ.1) GOTO 1040 C-------------------------------------------------------------- READ(5,INPUT,END=1040) C WRITE(6,120) 120 FORMAT(//' /INPUT/ DATA ARE --') WRITE(LABL,'(A80)') LABEL WRITE(6,130) LABL 130 FORMAT(/' RUN LABEL = ',A80) DO 140 IST=1,80 IF(TITLE(IST).NE.BL) GOTO 150 140 CONTINUE GOTO 190 150 DO 160 IND=1,80 IF(TITLE(81-IND).NE.BL) GOTO 170 160 CONTINUE GOTO 190 170 IND=81-IND NST=(119-IND+IST)/2 TIT(NST)=BL TIT2(NST)=BL DO 180 I=IST,IND NST=NST+1 TIT(NST)=TITLE(I) TIT2(NST)=TITLE(I) 180 CONTINUE TIT(NST+1)=BL TIT2(NST+1)=BL C 190 AMXKB=MX/128.D0 IF (NIPR.EQ.1.OR.NIPR.EQ.2) THEN WRITE(6,200) MX,CWD(NIPR),AMXKB 200 FORMAT(/' SCRATCH CORE STORAGE ALLOCATION IS',I10,A8, 1 ' WORDS (',F10.2,' KBYTES)') WRITE(6,202) NIPR 202 FORMAT(2X,I1,' INTEGER(S) CAN BE STORED IN EACH WORD.') ELSE WRITE(6,204) NIPR 204 FORMAT(/' *** ILLEGAL NIPR =',I10) ENDIF C PRINT=PRNTLV C C PROCESS INTFLG -- REQUESTED PROPAGATOR -- AND ITS INPUT DATA. C WRITE(6,210) INTFLG 210 FORMAT(/' INTEGRATOR REQUESTED BY INPUT VALUE INTFLG =',I3) 220 FORMAT(/' ***** ERROR - NO IMPLEMENTATION FOR THIS INTFLG' 1 ,' - RUN HALTED.') 240 FORMAT(/' COUPLED EQUATIONS SOLVED BY METHOD OF DEVOGELAERE.') 250 FORMAT(/' INTEGRATION PARAMETERS ARE RMIN =',F7.2/ 1 30X,'RMAX =',F7.2/30X,'STEST =',D11.2/30X,'STEPS =', 2 F6.1,' (PER WAVELENGTH)'/30X,'STABIL =',F6.1,' (STEPS PER', 3 ' STABILIZATION)') 270 FORMAT(/' COUPLED EQUATIONS SOLVED BY WALKER-LIGHT R-MATRIX', 1 ' PROPAGATOR ALGORITHM'//' PARAMETERS ARE',5X,'RMIN =', 2 F7.2,8X,'DR = ',G8.2/21X,'RMAX =',F7.2,8X, 3 'VTOL =',D9.2/21X,'RMID =',F7.2,8X,'MAXSTP =',I9) 271 FORMAT(/' RVFAC =',F7.2,' OVERRIDES INPUT RMID') 300 FORMAT(/' COUPLED EQUATIONS SOLVED BY LOG DERIVATIVE METHOD ', 1 'OF JOHNSON') 310 FORMAT(/' INTEGRATION PARAMETERS ARE RMIN =',F7.2,8X, 1 'STEPS = ',F7.1/33X,'RMAX =',F7.2) 320 FORMAT(/' CHANGING TO VARIABLE INTERVAL / VARIABLE STEP METHOD', 1 ' AT LONG RANGE'//' INTEGRATION PARAMETERS ARE RVIVAS =', 2 F7.2,8X,'DR =',G8.2/ 3 33X,'RMAX =',F7.2,8X,'DRMAX =',F8.2/ 4 56X,'ALPHA1 = ',F7.2/33X,'XSQMAX =',G7.1,8X,'ALPHA2 = ',F7.2/ 5 33X,'TOLHI =',G7.1,8X,'IALPHA =',I8/33X,'ISHIFT =',L7,8X, 6 'IV =',L8/33X,'IPERT =',L7,8X,'IVP =',L8/33X, 7 'IALFP =',L7,8X,'IVPP =',L8/33X,'ISYM =',L7,8X, 8 'NUMDER =',L8) 340 FORMAT(/' COUPLED EQUATIONS SOLVED BY DIABATIC ', 1 'MODIFIED LOG DERIVATIVE METHOD OF MANOLOPOULOS') 350 FORMAT(/' COUPLED EQUATIONS SOLVED BY QUASIADIABATIC ', 1 'MODIFIED LOG DERIVATIVE METHOD OF MANOLOPOULOS') 352 FORMAT(33X,'IABSDR =',I4) 353 FORMAT(33X,'OVERRIDES STEPS PARAMETER WITH DR =',F9.3) 354 FORMAT(/' AIRY PARAMETERS ','RMID =',F10.4/ 2 33X,'DRAIRY=',F10.4/33X,'TOLHI=',F13.6/ 3 33X,'POWRX =',F8.2) 355 FORMAT(/' DRAIRY.LT.0 TAKES INITIAL AIRY STEP SIZE FROM' 1 ,' MODIFIED LOG-DERIVATIVE VALUE.') 356 FORMAT(/' TOLHI.GE.1 -- AIRY STEP SIZE INCREASED BY' 1 ,' FACTOR OF TOLHI AT EACH STEP') 357 FORMAT(/' TOLHI.LT.1 -- AIRY STEPS ADJUSTED TO MAINTAIN' 1 ,' APPROX. ACCURACY VIA PERTURBATION THEORY AND POWRX.') 370 FORMAT(/' EQUATIONS SOLVED BY WKB APPROXIMATION WITH GAUSS-' 1 ,'MEHLER INTEGRATION. SEE R. T PACK, JCP 60, 633 (1974).'/ 2 /' NOTE THAT THIS IS IMPLEMENTED ONLY FOR ONE CHANNEL', 3 ' CASES, E.G., IOS CALCULATIONS.'/ 4 /' INTEGRATION PARAMETERS ARE RMIN =',D15.4/ 5 30X,'STEST =',D14.4/30X,'NGMP =',I6,' (',I2,')',I3) C IF(INTFLG.EQ.2) THEN WRITE(6,240) C STABIL=MIN(STABIL,STEPS/2.D0) WRITE(6,250) RMIN,RMAX,STEST,STEPS,STABIL GO TO 380 ENDIF C IF(INTFLG.EQ.3) THEN WRITE(6,270) RMIN,DR,RMAX,VTOL,RMID,MAXSTP IF(RVFAC.GT.0.D0 .AND. IRMSET.GT.0) WRITE(6,271) RVFAC GO TO 380 ENDIF C IF(INTFLG.EQ.4 .OR. INTFLG.EQ.5) THEN IF(IDIAG) THEN IV=.TRUE. IVP=.TRUE. IVPP=.TRUE. ISHIFT=.TRUE. IPERT=.TRUE. ENDIF IF(INTFLG.EQ.5) RVIVAS=RMAX WRITE(6,300) WRITE(6,310) RMIN,STEPS,RVIVAS IF(INTFLG.EQ.4) WRITE(6,320) RVIVAS,DR,RMAX,DRMAX,ALPHA1,XSQMAX, 1 ALPHA2,TOLHI,IALPHA,ISHIFT,IV,IPERT,IVP,IALFP,IVPP,ISYM,NUMDER GO TO 380 ENDIF C IF(INTFLG.EQ.6) THEN WRITE(6,340) WRITE(6,310) RMIN,STEPS,RMAX GO TO 380 ENDIF C IF(INTFLG.EQ.7) THEN WRITE(6,350) WRITE(6,310) RMIN,STEPS,RMAX GO TO 380 ENDIF C IF(INTFLG.EQ.8) THEN CALL MHAACK(6) WRITE(6,310) RMIN,STEPS,RMAX WRITE(6,352) IABSDR IF(IABSDR.EQ.1) WRITE(6,353) DR WRITE(6,354) RMID,DRAIRY,TOLHI,POWRX IF(RVFAC.GT.0.D0.AND.IRMSET.GT.0) WRITE(6,271) RVFAC IF(DRAIRY.LT.0.D0) WRITE(6,355) IF(TOLHI.GE.1.D0) THEN WRITE(6,356) ELSE WRITE(6,357) ENDIF GO TO 380 ENDIF C IF(INTFLG.EQ.-1) THEN WRITE(6,370) RMIN,STEST,NGMP GO TO 380 ENDIF C WRITE(6,220) STOP C 380 JKEEP=-1 XEPS=-1.D0 DEEP=1.D30 IF(IRXSET.GT.0) WRITE(6,381) IRXSET 381 FORMAT(/' IRXSET =',I3,' OPTION. RMAX ADJUSTED AUTOMATICALLY ', 1 'FOR EACH NEW JTOT,MVAL') IF(IRMSET.LE.0) GOTO 420 WRITE(6,390) IRMSET 390 FORMAT(/' IRMSET =',I3,' OPTION. RMIN CHOSEN AUTOMATICALLY ', 1 'FOR EACH NEW JTOT') C C XEPS IS SUCH THAT AIRY(XEPS) APPROX. EQUALS 10**(-IRMSET) C XEPS=(-1.5D0*LOG(4.D0*SQRT(PI)* 1 10.D0**(-IRMSET)))**(2.D0/3.D0) C>>SG 1/18/93 BELOW REMOVED AT SUGGESTION OF JMH C IF(ISCRU.EQ.0 .AND. NNRG.NE.1) SHRINK=0 IF(INTFLG.NE.3 .OR. SHRINK.NE.1) GOTO 420 DEEP=2.D0+XEPS**1.5D0/1.5D0 WRITE(6,400) 400 FORMAT(22X,'AND DEEPLY CLOSED CHANNELS ', 1 'DROPPED IN LONG-RANGE REGION') IF(NNRG.NE.1 .AND. ISCRU.NE.0) WRITE(6,410) 410 FORMAT(22X,'NOTE THAT BASIS SET CONTRACTION IS PERFORMED FOR ', 1 'ENERGY(1),'/22X,'SO THAT SUBSEQUENT ENERGIES MUST NOT BE ', 2 'SIGNIFICANTLY HIGHER.') C 420 ISAV=0 IF(JTOTL.EQ.JTOTU .AND. MSET.GT.0) ISAV=1 IF(ISCRU.LT.0) ISAV=-ISAV ISCRU=IABS(ISCRU) C IF(ISCRU.EQ.0) THEN IF(NNRG.GT.1.OR.NTEMP.GT.0) WRITE(6,430) 430 FORMAT(/' ***** WARNING - NO SCRATCH FILE SPECIFIED BY ISCRU ', 1 'PARAMETER - FULL CALCULATION WILL BE DONE AT EACH ENERGY') ELSE IF(ISAV.EQ.-1) THEN WRITE(6,440) ISCRU 440 FORMAT(/' ENERGY-INDEPENDENT MATRICES SAVED FROM A ', 1 'PREVIOUS RUN WILL BE READ FROM UNIT',I3) OPEN(ISCRU,FORM='UNFORMATTED',STATUS='OLD') ELSE WRITE(6,450) ISCRU 450 FORMAT(/' ENERGY-INDEPENDENT MATRICES WILL BE SAVED ', 1 'TEMPORARILY ON UNIT',I3) OPEN(ISCRU,FORM='UNFORMATTED',STATUS='UNKNOWN') ENDIF ENDIF C WRITE(6,470) URED 470 FORMAT(/' REDUCED MASS FOR COLLISION =',F14.9,' A.M.U.') IF(JTOTL.LT.0) JTOTL=0 IF(JTOTU.LT.JTOTL) JTOTU=999999 WRITE(6,480) JTOTL,JTOTU,JSTEP 480 FORMAT(/' CONTROL DATA FOR TOTAL ANGULAR MOMENTUM IS'/ 1 7X,'JTOT FROM',I4,' TO',I6,' IN STEPS OF',I4) IF(JTOTU.GE.999999) WRITE(6,490) NCAC,DTOL,OTOL 490 FORMAT(/' JTOT SERIES WILL BE TERMINATED WHEN MAX CHANGE IN ', 1 'CROSS SECTIONS IS LESS THAN TOLERANCE FOR NCAC =',I3, 2 ' CONSECUTIVE JTOT'/25X, 3 'DIAGONAL (DTOL) AND OFF-DIAGONAL (OTOL) TOLERANCES ARE',2F9.5) IF(JTOTU.GE.999999.AND.NNRGPG.GT.1) WRITE(6,491) NNRGPG 491 FORMAT(/' N.B. CONVERGENCE CHECKING IS DONE FOR ENERGY GROUPS', 1 ' OF NNRGPG =',I4) IF(MSET.GT.0 .AND. MHI.LE.0) MHI=MSET IF(MSET.GT.0) WRITE(6,500) MSET,MHI 500 FORMAT(/' CALCULATIONS WILL BE FOR SYMMETRY BLOCK ("PARITY ', 1 'CASES")',I4,' TO',I4) C C PROCESS TOTAL ENERGIES C CALL ECNV(EUNITS,EFACT) IF(NNRG.GT.0 .AND. DNRG.EQ.0.D0 .AND. ABS(EFACT-1.D0).GT.1.D-3 1 .AND. ICONV.EQ.0) WRITE(6,510) (ENERGY(I),I=1,NNRG) 510 FORMAT(/' INPUT ENERGY LIST IS'/(16X,7D16.6)) IF(NTEMP.LE.0) GOTO 520 C OVERRIDE ENERGY INPUT WITH TEMP INPUT NTEMP=MIN0(NTEMP,MXTEMP) CALL EAVG(NTEMP,TEMP,NGAUSS,ENERGY,NNRG,MXNRG) NPR=NNRG GOTO 590 520 ISRCH=0 NPR=NNRG C C PROCESS A NEGATIVE INPUT NNRG FOR RESONANCE SEARCH OPTION C IF(NNRG.GE.0 .OR. DNRG.EQ.0.D0 .OR. JTOTL.NE.JTOTU .OR. 1 MSET.LE.0 .OR. KSAVE.LE.0) GOTO 530 ISRCH=1 NNRG=5*(IABS(NNRG)/5) MXN=5*(MXNRG/5) NNRG=MIN0(NNRG,MXN) NNRGPG=5 NPR=5 C 530 NNRG=MIN0(MXNRG,NNRG) NPR=MIN0(MXNRG,NPR) IF(NNRG.GT.0) GOTO 550 WRITE(6,540) 540 FORMAT(/' ***** ERROR - NO INPUT ENERGIES SPECIFIED - RUN HALTED') STOP 550 IF(NNRG.LE.1 .OR. (DNRG.EQ.0.D0 .AND. ICONV.EQ.0)) GOTO 570 DO 560 I=2,NPR 560 ENERGY(I)=ENERGY(1)+(I-1)*DNRG 570 DO 580 I=1,NPR 580 ENERGY(I)=ENERGY(I)*EFACT 590 WRITE(6,600) NNRG 600 FORMAT(/' CONTROL DATA FOR TOTAL ENERGIES. CALCULATIONS WILL ', 1 'BE PERFORMED FOR',I4,' VALUES') DO 610 I=1,NPR ENEV=ENERGY(I)/8065.5410D0 610 WRITE(6,620) I,ENERGY(I),ENEV 620 FORMAT(7X,'ENERGY NO.',I4,' =',F17.9,' (1/CM) =',F17.12,' E.V.') C IF(ISRCH.EQ.1) WRITE(6,630) 630 FORMAT(/' RESONANCE SEARCH OPTION. ONLY FIRST 5 ENERGIES ', 1 'GIVEN. OTHERS WILL BE DETERMINED INTERACTIVELY.') C IF(IFLS.GT.0 .AND. IFEGEN.GT.0) WRITE(6,640) 640 FORMAT(/' THESE ENERGY VALUES WILL BE USED AS RELATIVE (CENTER', 1 ' OF MASS) VALUES AND LIST MAY BE MODIFIED ACCORDINGLY.') C IF(NUMDER) WRITE(6,641) 641 FORMAT(/' NUMDER=.TRUE. POTENTIAL DERIVATIVE WILL BE COMPUTED', & ' NUMERICALLY FROM POTENTIAL.') WRITE(6,650) PRINT,ISIGPR,ITHROW 650 FORMAT(/' PRINT LEVEL (PRNTLV) =',I3,' OTHER PRINT CONTROLS', 1 ' ISIGPR =',I2,' ITHROW =',I2) WRITE(6,660) 660 FORMAT(/' ',30('====')) C C INITIALIZE BASIS (BASIN/IOSBIN) C COMBINED MOLSCAT (BASIN) AND IOS (IOSBIN) -- APR 86 C IOSBIN GRABS STORAGE IN ATAU=JLEV=X (ITYPE=6 ONLY). MAX AVAILABLE C PASSED INITIALLY IN NLEV; SET6I/IOSBIN MUST UPDATE C IC ACCORDINGLY. N.B. IOS CASE ALSO USES NLEV TO PASS 'NVC' C FROM BASIN/IOSBIN TO IOSDRV. C BASIN TAKES STORAGE FOR JLEV=X, AND ALSO RESETS IC ACCORDINGLY; C FOR THIS CASE, NLEV INITIALIZED TO MAXIMUM AVAILABLE IN X(). IXJLEV=IXNEXT NLEV=MX C IXNEXT REMOVED FROM ARGUMENT LIST: JMH, 10 NOV 93 CALL BASIN(NLEV,X(IXJLEV),URED,NQN,NLABV(9),MXPAR,ITYPE,IOSFLG) C BASE ROUTINE INCREMENTS IXNEXT BY AMOUNT OF STORAGE IN JLEV. CALL CHKSTR(NUSED) WRITE(6,660) C C INITIALIZE POTENTIAL. C ILAM=IXNEXT MXLAM=NIPR*(MX-ILAM+1) CALL POTENL(-1,MXLAM,NPOTL,X(ILAM),RM,EPSIL,ITYPE) C THIS READS (5, POTL). RM AND EPSIL ARE SET HERE. C RM IS A LENGTH PARAMETER (IN ANGSTROMS) C EPSIL IS AN ENERGY PARAMETER IN WAVENUMBERS. ITYP=MOD(ITYPE,10) C INCREMENT IXNEXT FOR STORAGE TAKEN FOR LAM(NLABV,MXLAM) IXNEXT=IXNEXT+(MXLAM*NLABV(ITYP)+NIPR-1)/NIPR CALL IVCHK(IVLFL,PRNTLV,ITYPE,NLABV,MXLAM,NPOTL,X(ILAM)) WRITE(6,660) C C COMPUTE SOME DIMENSIONLESS PARAMETERS C C RMLMDA IS THE SQUARE OF THE RATIO OF RM TO DEBROGLIE WAVELENGTH RMLMDA=URED*RM*RM*EPSIL/BFCT C CINT IS THE FACTOR TO REDUCE THE ROTATIONAL CONSTANTS CINT = RMLMDA/EPSIL C C *** IF(IOSFLG.LE.0) GOTO 670 C *** C *** THIS IS WHERE IOS CODE DIVERGES - CALL IOS CODE AND SKIP TO EXIT C *** IF (IRSTRT.NE.0) THEN WRITE(6,*) ' *** RESTART REQUESTED WITH IOS RUN - NOT ALLOWED' WRITE(6,*) ' *** MODIFY INPUT DECK AND RESUBMIT' STOP ENDIF CALL IOSDRV(NNRG,NPR,ENERGY,JTOTL,JTOTU,JSTEP,TEST,NCAC, 1 IFLS,LINE,LTYPE,MXLN,INTFLG,ITYPE,LMAX,MMAX, 2 IPROGM,URED,LABL,NUMDER, 3 X(ILAM),MXLAM,NPOTL,CINT,IRMSET,IRXSET,RVFAC, 4 DEEP,PRINT,NLEV,ISAVEU,TFIRST,RM,EPSIL,RMIN,RMAX) CALL GCLOCK(TLAST) TOTIME=TLAST-TFIRST GOTO 1020 C C PROCESS PRESSURE-BROADENING LINE-SHAPE INPUT PARAMETERS. C 670 IF(IFLS.GT.0) THEN CALL PRBRIN(IFLS,LINE,LTYPE,MXLN,ILSU,NNRG,ENERGY,MXNRG,IFEGEN, 1 X(IXJLEV),PRINT) IF(IFEGEN.GT.0) NPR=NNRG WRITE(6,660) IF(KSAVE.EQ.0) GOTO 690 WRITE(6,680) IFLS,KSAVE 680 FORMAT(/' ****** WARNING. IFLS =',I3,' AND KSAVE =',I3,' ARE ', 1 'INCOMPATIBLE. KSAVE IS RESET TO ZERO') KSAVE=0 ENDIF C C INITIALIZE OUTPUT ROUTINE. C OUTPUT TAKES AN ADDITIONAL AMOUNT OF STORAGE C FOR SIG AT X(IXNEXT) AND INCREASES IXNEXT ACCORDINGLY. C 690 IOUT=IXNEXT C NOTE THAT IXNEXT WILL BE CHANGED BY OUTINT CALL OUTINT(LABL,ENERGY,NNRG,NLEV,NQN,X(IXJLEV),X(IOUT),IXNEXT, 1 IECONV,URED,ITYPE,KSAVE,ISST,MINJT,MAXJT,ISIGU,IPARTU,ISAVEU, 2 IPROGM,MXSIG,ISIGPR,JSTEP,IRSTRT) CALL CHKSTR(NUSED) IC1=IXNEXT C PROCESS RESTART REQUEST ... MXP=0 CALL RESTRT(IRSTRT,ISAVEU,JTOTL,JSTEP,MXPAR,MSET,MHI, 1 LABEL,ITYPE,NLEV,NQN,URED,IPROGM, 2 X(IXJLEV),NNRG,ENERGY,MXNRG, 3 X(IOUT),ISST,IECONV,MINJT,MAXJT,ISIGU,IPARTU,KSAVE, 4 OTOL,DTOL,X(IC1),X(IC1),MRSTRT,IERST,MXP,PRINT) WRITE(6,660) C EFIRST=ENERGY(1)*CINT CALL GCLOCK(TITIME) TTIME=TITIME-TFIRST WRITE(6,700) TTIME,NUSED 700 FORMAT(/' INITIALIZATION DONE. TIME WAS',F7.2,' CPU SECS.',I10, 1 ' WORDS OF STORAGE USED.') IF(PRINT.LT.4) WRITE(6,710) TIT 710 FORMAT('1',120A1) IF(PRINT.GE.4.AND.ITHROW.EQ.0) WRITE(6,720) 720 FORMAT('1') C C ************** LOOP OVER JTOT VALUES BEGINS HERE. ****************** C DO 990 JTOT=JTOTL,JTOTU,JSTEP IF(PRINT.GE.1 .AND. PRINT.LE.4) WRITE(6,730) JTOT 730 FORMAT(/' ANGULAR MOMENTUM JTOT =',I4/2X,7('****')) THETA=THETLW+THETST*DBLE(JTOT) C C *************** LOOP OVER SYMMETRY BLOCKS BEGINS HERE ************** C DO 980 M=1,MXPAR PHI=PHILW+PHIST*DBLE(M-1) IF(MSET.GT.0.AND.(M.LT.MSET.OR.M.GT.MHI)) GO TO 980 IF (IRSTRT.GE.2.AND.JTOT.EQ.JTOTL.AND.M.LT.MRSTRT) THEN WRITE(6,736) M,IRSTRT 736 FORMAT(' *** SKIPPING MVALUE =',I3,' DUE TO IRSTRT =',I3) GO TO 980 ENDIF IF(PRINT.LT.4) GOTO 760 IF(ITHROW.NE.0) WRITE(6,710) TIT IF(ITHROW.EQ.0) WRITE(6,740) TIT 740 FORMAT(/' ',120A1) WRITE(6,750) JTOT,M 750 FORMAT(/' TOTAL ANGULAR MOMENTUM, JTOT =',I5,' SYMMETRY', 1 ' BLOCK =',I4) 760 CONTINUE C C CHOOSE BASIS FUNCTIONS C CALL BASE (JTOT,X(IXJLEV),N,X,X,CINT,X,X,X,X,MXLAM,NPOTL,X(ILAM), 1 X,WGHT,IEXCH,THETA,PHI,M,.TRUE.,EFIRST,NLEV,PRINT) C C MOLD IS A REMNANT OF THE PREVIOUS "PARITY CASE" PROCESSING. C MXP IS USED IN CONVERGENCE CHECKING, MOLD IS PASSED TO PRBR C MOLD=-M IF(M.EQ.MXPAR.AND.N.LE.0) MOLD=0 MXP=MAX0(MXP,IABS(MOLD)) IF(M.EQ.MXPAR) MOLD=0 C C N IS THE NUMBER OF BASIS FUNCTIONS C SKIP THIS JTOT,M IF NO CHANNELS C C IF(N.LE.0) GOTO 980 <<- SG: FIXES ISIGU BUG IF(N.LE.0) GOTO 770 NSQ = N*N C C ALLOCATE STORAGE FOR COUPLED EQUATION SOLVER. C C ALLOCATE STORAGE COMMON TO ALL SCATTERING. . . C IS0-IS9 ARE SREAL,SIMAG,K-MATRIX,VL,IV,EINT,CENT,WVEC,L,NBASIS C NOTE THAT INTEGER ARRAYS OF LENGTH N ARE NOT REDUCED BY NIPR C IC1 IS IXNEXT AFTER ALLOCATIONS OF BASIN, POTENL, OUTINT ... ISJ=IC1 IS0=ISJ+N IS1=IS0+NSQ IS2=IS1+NSQ IS3=IS2+NSQ NV=N*(N+1)/2 IF(IVLU.EQ.0) NV=NV*NPOTL IS4=IS3+NV IS5=IS4 IF(IVLFL.GT.0) IS5=IS4+(NV+NIPR-1)/NIPR IS6=IS5+N IS7=IS6+N IS8=IS7+N IS9=IS8+N IXNEXT=IS9+N C C SET UP SOME STORAGE POINTERS FOR LATER USE IN CONVRG C IF(ICONV.GT.0) THEN IS10=IXNEXT IS11=IS10+NSQ IXNEXT=IS11+NSQ ENDIF IC2=IXNEXT CALL CHKSTR(NUSED) C IXNEXT/IC2 REFLECT STORAGE ALWAYS NEEDED FOR THIS JTOT,PARITY. C C SET UP BASIS FUNCTIONS IN ALLOCATED STORAGE C CALL BASE(JTOT,X(IXJLEV),N,X(ISJ),X(IS8),CINT,X(IS5),X(IS6), 1 X(IS3),X(IS4),MXLAM,NPOTL,X(ILAM),X(IS7),WGHT,IEXCH,THETA, 2 PHI,M,.FALSE.,EFIRST,NLEV,PRINT) C C CHECK THAT RMAX IS BEYOND CENTRIFUGAL BARRIER C CALL FINDRX(ENERGY,X(IS5),X(IS6),NPR,N,CINT,RMAX,RSTOP, 1 NOPMAX,IRXSET,PRINT) IF(INTFLG.EQ.5) RVIVAS=RSTOP RSTART=RMIN C C ****************** LOOP OVER ENERGIES BEGINS HERE ****************** C 770 NELOOP=(NNRG+NNRGPG-1)/NNRGPG JHI=0 ICODE=0 ALDONE=.TRUE. DO 966 IEL=1,NELOOP JLO=JHI+1 JHI=MIN(JHI+NNRGPG,NNRG) C C SEE WHETHER THIS BLOCK OF ENERGIES CAN BE SKIPPED C LCALC=.FALSE. DO 775 J=JLO,JHI IF(IECONV(J).LT.0 .AND. IECONV(J).GT.-2*MXP) THEN WRITE(6,772) JTOT,J 772 FORMAT(/' * * * WARNING. JTOT =',2I5,'-TH ENERGY PREVIOUSLY ', 1 'FAILED TO CONVERGE.') LCALC=.TRUE. ELSEIF(IECONV(J).EQ.0) THEN LCALC=.TRUE. ELSEIF(IECONV(J).GT.0) THEN IF(JTOTU.LT.999999 .OR. IECONV(J).LT.NCAC*MXP) LCALC=.TRUE. ENDIF 775 CONTINUE C IF(.NOT.LCALC) GOTO 966 ALDONE=.FALSE. DO 960 J=JLO,JHI IF(N.LE.0) THEN CALL OUTSIG(ISIGU,M,MXPAR,J,ENERGY,MINJT,MAXJT,X(IOUT)) GOTO 960 ENDIF C C IF THIS IS A PRESSURE BROADENING CALC AND THIS S-MATRIX C WILL NOT BE USED, SKIP IT C IF(IFLS.GT.0 .AND. IFEGEN.GE.2) THEN CALL PRBCNT(J,X(ISJ),N,IUSE) IF(IUSE.EQ.0) THEN LWARN=.TRUE. IF(PRINT.GE.2) WRITE(6,777) JTOT,M,J,ENERGY(J) 777 FORMAT(/' ****** S MATRIX FOR JTOT =',I5,' M =',I4,3X, 1 'ENERGY(',I3,') =',F18.9,/9X,'WILL NOT BE USED ', 2 'IN PRESSURE BROADENING CALCULATION: SKIPPING') IF(IECONV(J).GE.0) IECONV(J)=IECONV(J)+1 CALL OUTSIG(ISIGU,M,MXPAR,J,ENERGY,MINJT,MAXJT,X(IOUT)) GOTO 960 ENDIF ENDIF C IF (IRSTRT.EQ.3.AND.JTOT.EQ.JTOTL.AND.M.EQ.MRSTRT.AND.J.LT.IERST) 1 GO TO 960 ETOT=ENERGY(J) ERED=ETOT*CINT IF(ICODE.EQ.0) THEN EFIRST=ERED ICODE=1 ENDIF ESHIFT=ERED-EFIRST C C ICODE CONTROLS WHETHER POTENTIAL INFORMATION IS READ FROM CHANNEL C ICODE=1 CALCULATES INFORMATION AND STORES IT C ICODE=2 (SET AFTER 1ST ENERGY) READS STORED INFORMATION C IF(PRINT.GE.4) THEN IF(ITHROW.EQ.0) THEN WRITE(6,740) TIT2 ELSE WRITE(6,710) TIT2 ENDIF WRITE(6,780) JTOT,M,J,ETOT 780 FORMAT(/' JTOT =',I5,' SYMMETRY BLOCK =',I4,' ENERGY(', 1 I3,') =',F18.9,' (1/CM)') ENDIF C C FOR SURFACE SCATTERING AT SUBSEQUENT ENERGY, C GET CORRESPONDING THETA FOR PRINTING C IF(ITYPE.EQ.8 .AND. J.NE.1) THEN SINTH=SIN(THETA*PI/180.D0) SINTH=SINTH**2*ENERGY(1)/ETOT IF(SINTH.GT.1.D0) GOTO 960 THETJ=ASIN(SQRT(SINTH))*180.D0/PI WRITE(6,795) J,ETOT,THETJ 795 FORMAT(/' NOTE: K VECTORS PARALLEL TO SURFACE WERE ', 1 'CALCULATED FOR ENERGY(1)'/' SUBSEQUENT ENERGY(',I3,') =', 2 F10.4,' CORRESPONDS TO THETA =',F10.4,' DEGREES') ENDIF C C TEMPORARY STORAGE FOR HEADER, FINDRM C IT1=IXNEXT IT2=IT1+MXLAM IT3=IT2+N IT4=IT3+N IT5=IT4+N IXNEXT=IT5+N CALL CHKSTR(NUSED) C CALL HEADER(X(IS1),X(IS2),N,NSQ,X(IT1),X(IS3),X(IS4),X(IS5), 1 X(IS6),X(IT2),MXLAM,NPOTL,ICODE,ISAV,EFIRST) C IF(ICODE.EQ.1 .AND. IRMSET.GT.0) THEN C FOR IRMSET > 0 OPTION, CHOOSE APPROPRIATE RMIN RSTART=RMIN CALL FINDRM(X(IS1),N,RSTART,RTURN,IK,X(IT1),X(IS3),X(IS4),ERED, 1 X(IS5),X(IS6),RMLMDA,X(IT2),X(IT3),X(IT4),X(IT5),MXLAM,NPOTL, 2 IRMSET,ITYPE,PRINT) IF(RVFAC.NE.0.D0) THEN RMID=RVFAC*RTURN IF(PRINT.GE.3.AND.RSTOP.GT.RMAX) WRITE(6,799) RSTOP,RMAX 799 FORMAT(' RMID OBTAINED FROM RTURN EVEN THOUGH RSTOP.GT.RMAX', 1 2F8.2) IF(PRINT.GE.3) WRITE(6,800) RMID,RVFAC 800 FORMAT(/' RMID =',F7.2,' OBTAINED FROM RVFAC =',F6.3) ENDIF ELSE RTURN=RMIN IK=1 ENDIF C C RESET IXNEXT TO RECOVER TEMPORARY STORAGE FROM HEADER AND FINDRM C IXNEXT=IT1 C C SOLVE COUPLED EQUATIONS. C PROPAGATORS ARE CALLED FROM SUBROUTINE STORAG C CALL STORAG(INTFLG,N,MXLAM,NV,NPOTL, 1 ISJ,IS0,IS1,IS2,IS3,IS4,IS5,IS6,IS7,IS8,IS9, 2 ESHIFT,NOPMAX,DEEP,IK,ICODE,PRINT,NUMDER) C CALL GCLOCK(TJTIME) TTIME=(TJTIME-TITIME) TITIME=TJTIME C IF(NOPEN.GT.0) THEN C RESET ICODE TO ALLOW "SUBSEQUENT ENERGY" CALCULATIONS ICODE=2 ELSE IF(PRINT.GE.4) WRITE(6,900) JTOT,M,J,ETOT,TTIME 900 FORMAT(/' ****** NO OPEN CHANNELS FOR JTOT =',I5,3X, 1 'M =',I4,' ENERGY(',I3,') =',F18.9,10X,'STEP TIME =', 2 F6.2,' SECS') IF(IECONV(J).GE.0) IECONV(J)=IECONV(J)+1 GOTO 960 ENDIF C C FORCE IRSTRT=0 SO THAT ISAVEU WILL BE UPDATED. IRSTRX=0 CALL OUTPUT(JTOT,X(IS9),X(ISJ),X(IS8),X(IS7),X(IS0),X(IS1), 1 X(IS2),CONV,NOPEN,M,MXPAR,WGHT,IEXCH,J,RM,PRINT,TTIME, 2 ENERGY,X(IOUT),X(IXJLEV),ISST,IECONV,MINJT,MAXJT, 3 NLEV,NQN,OTOL,DTOL,KSAVE,ISIGU,IPARTU,ISAVEU,ISIGPR,IRSTRX) C IF(ICONV.GT.0) CALL CONVRG(J,X(IS0),X(IS1),X(IS10),X(IS11)) IF(IECONV(J).LT.0 .OR. IFLS.LE.0) GOTO 940 C C TEMPORARY STORAGE FOR PRBR -- THESE ARE INTEGERS, COULD USE NIPR IT1=IXNEXT IT2=IT1+N IT3=IT2+N IT4=IT3+N IXNEXT=IT4+N CALL CHKSTR(NUSED) CALL PRBR(JTOT,MOLD,NOPEN,J,RM, 1 X(IS9),X(ISJ),X(IS8),X(IS7), 2 X(IS0),X(IS1),X(IT1),X(IT2),X(IT3),X(IT4), 3 X(IXJLEV),MXPAR,WGHT,PRINT,ILSU) C RECOVER TEMPORARY STORAGE ... IXNEXT=IT1 C 940 IF(PRINT.GE.5) WRITE(6,950) JTOT,M,J,ETOT,TTIME 950 FORMAT(/' FINISHED JTOT =',I5,' M =',I4,' ENERGY(',I3, 1 ') =',F18.9,10X,'STEP TIME =',F8.2,' SECS') C 960 CONTINUE C C RESONANCE SEARCH OPTION - GENERATE NEXT 5 ENERGIES C IF(ISRCH.EQ.0) GOTO 964 CALL NEXTE(ENERGY(JLO),EPSM,ENEW,DNRG,KSAVE) IF(JHI.EQ.NNRG) GOTO 964 IF(ENEW.LE.0.D0) GOTO 1000 JST=JHI+1 JND=JHI+5 WRITE(6,600) NNRG DO 962 JJ=JST,JND ENERGY(JJ)=ENEW+(JJ-JST)*DNRG ENEV=ENERGY(JJ)/8065.541D0 WRITE(6,620) JJ,ENERGY(JJ),ENEV 962 CONTINUE 964 CONTINUE C 966 CONTINUE C C ******************** END OF LOOP OVER ENERGIES ********************* C IF(ALDONE) THEN WRITE(6,968) DTOL,OTOL,NCAC 968 FORMAT(///' CALCULATION TERMINATED BY CONVERGENCE OF TOTAL ', 1 'CROSS SECTIONS.'//' DIAGONAL AND OFF-DIAGONAL TOLERANCES ', 2 'WERE ',2F9.5,' NCAC =',I3) GOTO 1000 ENDIF C IF(PRINT.GE.2 .AND. PRINT.LT.5) WRITE(6,970) 970 FORMAT(/) C RESTORE ERED TO FIRST ENERGY VALUE. ERED = EFIRST 980 CONTINUE C C ****************** END OF LOOP OVER SYMMETRY BLOCKS **************** C IF(IFLS.GT.0) CALL PRBOUT(JSTEP) 990 CONTINUE C C ******************** END OF LOOP OVER JTOT VALUES ****************** C C END OF RUN BOOKKEEPING C 1000 CALL OUTPCH(X(IOUT),ENERGY,NNRG,MINJT,MAXJT,ISIGPR,LABL,ISIGU, 1 LWARN) IF(IFLS.GT.0) WRITE(6,710) TIT IF(IFLS.GT.0) CALL PRBOUT(JSTEP) IF(IFLS.GT.0) CALL DACLOS CALL GCLOCK(TLAST) TOTIME=TLAST-TFIRST C MAKE SURE WE HAVE NUSED FOR KSAVE BY CALLING CHKSTR CALL CHKSTR(NUSED) IF(KSAVE.GT.0) WRITE(KSAVE,1010) TOTIME,TTIME,NUSED 1010 FORMAT(/' TOTAL CPU =',F9.2,' SECS LAST CYCLE =', 1 F8.2,' SECS NUSED =',I8) C C *** IOS CALCULATION (IOSFLG.GT.0) REJOINS CODE BELOW C ASCERTAIN 'HIGH-WATER' MARK IN STORAGE FROM CHKSTR. C MX MAY HAVE BEEN REDUCED, SO USE MXSAVE FOR ALLOCATED STORAGE C 1020 CALL CHKSTR(NUSED) WRITE(6,1030) IPROGM,PDATE,TOTIME,NUSED,MXSAVE 1030 FORMAT(///' ',8('----MOLSCAT----')/' |',120X,'|'/' |',13X, 1 'COUPLED CHANNEL MOLECULAR SCATTERING PROGRAM OF J. M. HUTSON ', 2 'AND S. GREEN, VERSION',I3,1X,A8,13X,'|'/ 3 ' |',120X,'|'/' |',13X,'THIS RUN USED',F11.2,' CPU SECS ', 4 'AND',I10,' OF THE ALLOCATED',I10,' WORDS OF STORAGE',14X, 5 '|'/' |',120X,'|'/' ',8('----MOLSCAT----') ) IF(LASTIN.EQ.0) GOTO 100 1040 RETURN END SUBROUTINE AIRPRP (Z, W, TMAT, VECNOW, VECNEW, + EIGOLD, EIGNOW, HP, Y1, Y2, CC, Y4, XF, REND, DRNOW, EN, + TOLAI, POWR, ESHIFT, NCH, ITWO, IREAD, IWRITE, IPRINT, $ ISCRU, P, MXLAM, VL, IV, RMLMDA, EINT, CENT, NPOTL) C * AIRY ZEROTH-ORDER PROPAGATOR FROM R=XF TO R=REND * FOR REFERENCE SEE M. ALEXANDER, "HYBRID QUANTUM SCATTERING ALGORITHMS * J. CHEM. PHYS. 81, 4510 (1984) * AND M. ALEXANDER AND D. MANOLOPOULOS, "A STABLE LINEAR * REFERENCE POTENTIAL ALGORITHM FOR SOLUTION ..." * J. CHEM. PHYS. 86, 2044 (1987) * AUTHOR: MILLARD ALEXANDER * CURRENT REVISION DATE: 4-FEB-1991 * ---------------------------------------------------------------------- * ADAPTED TO MOLSCAT 4/91 BY TRP@NASAGISS * ADAPTED TO MOLSCAT VERSION 11 BY JMH, JUN 92 *----------------------------------------------------------------------- * DEFINITION OF VARIABLES IN CALL LIST: * Z: MATRIX OF MAXIMUM DIMENSION NCH*NCH * ON ENTRY Z CONTAINS THE INITIAL Z-MATRIX AT R=XF * ON RETURN Z CONTAINS THE Z-MATRIX AT R=REND * W, TMAT, VECNOW * , VECNEW: SCRATCH MATRICES OF DIMENSION AT LEAST NCH*NCH * EIGOLD, EIGNOW * , HP, Y1, Y2 * , CC, Y4: SCRATCH VECTORS OF DIMENSION AT LEAST NCH * XF: ON ENTRY: CONTAINS INITIAL VALUE OF INTERPARTICLE D * ON EXIT: CONTAINS FINAL VALUE OF INTERPARTICLE DIS * THIS IS EQUAL TO REND IF NORMAL TERMINATI * OTHERWISE AN ERROR MESSAGE IS PRINTED * DRNOW: ON ENTRY: CONTAINS INITIAL INTERVAL SIZE * ON EXIT: CONTAINS FINAL INTERVAL SIZE * EN: COLLISION ENERGY IN ATOMIC UNITS * TOLAI: PARAMETER TO DETERMINE STEP SIZES * IF TOLAI .LT. 1, THEN ESTIMATED ERRORS ARE USED TO * DETERMINE NEXT STEP SIZES FOLLOWING THE PROCEDURE O * IN M.H. ALEXANDER, "HYBRID QUANTUM SCATTERING ALGOR * IF TOLAI .GE. 1, THEN STEP SIZES ARE CONTROLLED BY * ALGORITHM: DRNEXT = TOLAI * DRNOW * POWR: POWER AT WHICH STEP SIZES INCREASE C C * LOGICAL VARIABLES: * ISYM: IF .TRUE., PROPAGATION ASSUMES SYMMETRY OF Y MATRIX * ---------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL ISYM INTEGER I, IEND, IERR, ITWO, IZERO, KSTEP, MAXSTP, : NCH, NPT, NSKIP * REAL CDIAG, CMAX, COFF, DRFIR, DRMID, DRNOW, EN, ESHIFT, FACT, * : ONE, POWR, REND, RLAST, RMIN, RNEW, RNEXT, RNOW, ROLD, * : SPCMN, SPCMX, TOLAI, XF, XLARGE, ZERO * REAL CC, EIGNOW, EIGOLD, HP,Y1, Y2, Y4 * REAL TMAT, VECNEW, VECNOW, W, Z EXTERNAL CORR, TRNSFM, OUTMAT, POTENT, DAXPY, DCOPY, : SYMINV, SPROPN, DSCAL, TRNSP, WAVEIG * MATRIX DIMENSIONS (ROW DIMENSION = NCH, MATRICES STORED COLUMN BY CO DIMENSION Z(1), W(1), TMAT(1), VECNOW(1), VECNEW(1) * VECTORS DIMENSIONED NCH DIMENSION EIGOLD(1), EIGNOW(1), HP(1), Y1(1), Y2(1), CC(1), Y4(1) DIMENSION P(1),VL(1),IV(1),EINT(1),CENT(1) C DATA IZERO, IONE, ZERO, ONE /0, 1, 0.D0, 1.D0/ DATA ISYM /.TRUE./ C * ---------------------------------------------------------------------- C ERED = EN RMIN = XF SPCMX = 0.D0 SPCMN = 0.D0 IF (ITWO .GT. 0) GO TO 60 SPCMN = REND - RMIN * DETERMINE LOCAL WAVEVECTORS AT RMIN TO USE IN ESTIMATING SECOND DERIV * HP AND Y1 ARE USED AS SCRATCH VECTORS HERE CALL WAVEIG (W, EIGOLD, HP, Y1, RMIN, NCH, P, MXLAM, VL, IV, 1 RMLMDA, ERED, EINT, CENT, NPOTL) * LOCAL WAVEVECTORS AT RMIN ARE RETURNED IN EIGOLD DRFIR = DRNOW DRMID = DRNOW * 0.5D0 RLAST = XF ROLD = XF RNOW = RLAST + DRMID RNEXT = RLAST + DRNOW * DEFINE LOCAL BASIS AT RNOW AND CARRY OUT TRANSFORMATIONS * VECNEW IS USED AS SCRATCH MATRIX AND Y1 IS USED AS SCRATCH VECTOR HER CALL POTENT (W, VECNOW, VECNEW, EIGNOW, HP, Y1, + RNOW, DRNOW, EN, XLARGE, NCH, $ P, MXLAM, VL, IV, RMLMDA, ERED, EINT, CENT, NPOTL) * VECNOW IS TRANSFORMATION FROM FREE BASIS INTO LOCAL BASIS * IN FIRST INTERVAL * E.G. P1=VECNOW ; SEE EQ.(23) OF * M.H. ALEXANDER, "HYBRID QUANTUM SCATTERING ALGORITHMS ..." * STORE VECNOW IN TMAT CALL DCOPY (NCH*NCH, VECNOW, 1, TMAT, 1) * DETERMINE APPROXIMATE VALUES FOR DIAGONAL AND OFF-DIAGONAL * CORRECTION TERMS CALL CORR (EIGNOW, EIGOLD, HP, DRNOW, DRMID, XLARGE, CDIAG, : COFF, NCH) MAXSTP = ( (REND-XF) / DRNOW ) * 5 XF = REND IF (IPRINT.GT.40) THEN WRITE (6, 40) 40 FORMAT(/' ** AIRY PROPAGATION (NO DERIVATIVES):') WRITE (6, 50) 50 FORMAT(' STEP RNOW', 5X, 5HDRNOW, 5X, 5HCDIAG, 6X, 4HCOFF) END IF 60 IEND = 0 IF (ITWO .LT. 0) GO TO 70 IF (ITWO .EQ. 0) WRITE(ISCRU) MAXSTP IF (ITWO. GT. 0) READ(ISCRU) MAXSTP * WRITE OR READ RELEVANT INFORMATION CALL OUTMAT (TMAT, EIGOLD, HP, ESHIFT, DRNOW, RNOW, : NCH, NCH, ITWO, ISCRU) C C * START AIRY PROPAGATION C * ---------------------------------------------------------------------- 70 DO 200 KSTEP = 1, MAXSTP NSTEP=KSTEP C * TRANSFORM LOG-DERIV MATRIX FROM LOCAL BASIS IN LAST INTERVAL TO * LOCAL BASIS IN PRESENT INTERVAL. SEE EQ.(23) OF * M.H. ALEXANDER, "HYBRID QUANTUM SCATTERING ALGORITHMS ..." * W IS USED AS SCRATCH MATRIX HERE, AND Y1 IS SCRATCH ARRAY C CALL TRNSP ( TMAT, NCH) CALL TRNSFM ( TMAT, Z, W, NCH, .FALSE., ISYM ) C * TMAT IS NO LONGER NEEDED * SOLVE FOR LOG-DERIVATIVE MATRIX AT RIGHT-HAND SIDE OF * PRESENT INTERVAL. THIS USES NEW ALGORITHM OF MANALOPOULOS AND ALEXAN * NAMELY * (N) (N) -1 (N) (N) * Z = - Y [ Y + Z ] Y + Y * N+1 2 1 N 2 4 * WHERE Y , Y , AND Y ARE THE (DIAGONAL) ELEMENTS OF THE "IMBEDDING * 1 2 4 * PROPAGATOR DEFINED IN ALEXANDER AND MANOLOPOULOS * DETERMINE THESE DIAGONAL MATRICES FOR PROPAGATION OF LOG-DERIV MATRIX * EQS. (38)-(44) OF M. ALEXANDER AND D. MANOLOPOULOS, "A STABLE LINEAR * REFERENCE POTENTIAL ALGORITHM FOR SOLUTION ..." C CALL SPROPN ( DRNOW, EIGOLD, HP, Y1, Y4, Y2, NCH) C * SET UP MATRIX TO BE INVERTED * NSKIP IS SPACING BETWEEN DIAGONAL ELEMENTS OF MATRIX STORED COLUMN BY C NSKIP = NCH + 1 CALL DAXPY (NCH, ONE, Y1, 1, Z, NSKIP) C * INVERT (Y + Z ) * 1 N C CALL SYMINV (Z, NCH, NCH, IERR) CALL DSYFIL ('U', NCH, Z, NCH) IF (IERR .GT. NCH) THEN WRITE (6, 80) 80 FORMAT (' *** INSTABILITY IN SYMINV IN AIRPRP.') STOP END IF C * -1 * EVALUATE - Y ( Y + Z ) Y * 2 1 N 2 * IN THE NEXT LOOPS EVALUATE THE FULL, RATHER THAN LOWER TRIANGLE C NPT = 1 DO 85 I = 1, NCH FACT = Y2(I) CALL DSCAL (NCH, FACT, Z(NPT), 1) NPT = NPT + NCH 85 CONTINUE * -1 * Z NOW CONTAINS ( Y + Z ) Y , THIS IS G(N-1,N) IN THE LOCAL BASI * 1 N 2 C DO 110 I = 1, NCH FACT = - Y2(I) CALL DSCAL (NCH, FACT, Z(I), NCH) 110 CONTINUE C * ADD ON Y * 4 CALL DAXPY (NCH, ONE, Y4, 1, Z, NSKIP) C IF (ITWO .GT. 0) GO TO 160 C C * OBLIGATORY WRITE OF STEP INFORMATION IF DEVIATIONS FROM LINEAR * POTENTIAL ARE UNUSUALLY LARGE * THIS IS ONLY DONE IF TOLAI .LT. 1, IN WHICH CASE THE LARGEST CORRECTI * IS USED TO ESTIMATE THE NEXT STEP C IF (TOLAI .LT. 1.) THEN CMAX = MAX (CDIAG, COFF) IF (CMAX .GT. (5. * TOLAI)) THEN WRITE (6,125) 125 FORMAT : (' ** ESTIMATED CORRECTIONS LARGER THAN 5*TOLAI IN AIRPRP') IF (KSTEP .EQ. 1) THEN WRITE (6, 130) 130 FORMAT (' THE INITIAL VALUE OF DRNOW (SPAC*FSTFAC) IS', : ' PROBABLY TOO LARGE') ELSE WRITE (6, 140) 140 FORMAT : (' CHECK FOR DISCONTINUITIES OR UNPHYSICAL OSCILLATIONS', : /,' IN YOUR POTENTIAL') END IF IF (IPRINT.LT.41) THEN WRITE (6, 50) WRITE (6,150) KSTEP, RNOW, DRNOW, CDIAG, COFF END IF END IF END IF C C * WRITE OUT INFORMATION ABOUT STEP JUST COMPLETED C IF (IPRINT.GT.40) THEN WRITE (6,150) KSTEP, RNOW, DRNOW, CDIAG, COFF 150 FORMAT (I6, 4E10.3) END IF C C * GET SET FOR NEXT STEP C 160 IF (IEND .EQ. 1) GO TO 250 IF (ITWO .GT. 0) GO TO 180 C C * IF TOLAI .LT. 1, PREDICT NEXT STEP SIZE FROM LARGEST CORRECTION C IF (TOLAI .LT. 1.) THEN C * NOTE THAT THE FOLLOWING STATEMENT IS SLIGHTLY DIFFERENT FROM EQ. (30 * OF M.H. ALEXANDER, "HYBRID QUANTUM SCATTERING ALGORITHMS ... AND THA * THE STEP-SIZE ALGORITHM IS ONLY APPROXIMATELY RELATED TO ANY REAL * ESTIMATE OF THE ERROR COFF AND CDIAG SHOULD BE APPROXIMATELY TOLAI, S * FROM EQ. (27): * DRNOW(AT N+1) = (12 TOLAI/KBAR(N+1)W(N+1)-TILDA')**(1/3) * WHICH IS APPROXIMATELY = (12 TOLAI/KBAR(N)W(N)-TILDA')**(1/3) * = ((12 COFF/KBAR W-TILDA') (TOLAI/COFF))**(1/3 * = DRNOW(AT N) (TOLAI/COFF)**(1/3) * OR FROM EQ. (29): * DRNOW = DRNOW (TOLAI/CDIAG)**(1/3) * THEN, USING THE LARGER ERROR AND ALLOWING POW TO VARY: C CSG>> FACTOR=(TOLAI/CMAX) ** (1./POWR) CSG LIMIT INCREMENT/DECREMENT FOR STABILITY ... IF (FACTOR.GT.2.D0) FACTOR=2.D0 IF (FACTOR.LE.0.1D0) FACTOR=1.D-1 DRNOW = DRNOW * FACTOR CSG<< DRNOW = DRNOW * (TOLAI/CMAX) ** (1. / POWR) ELSE C * IF TOLAI .GE. 1, THEN * MINIMUM STEP SIZE IS FIRST INTERVAL WIDTH C IF (KSTEP .EQ. 1) SPCMN = DRNOW C * AND NEXT STEP SIZE IS TOLAI * PRESENT STEP SIZE C DRNOW = TOLAI * DRNOW END IF C * DRNOW IS STEP SIZE IN NEXT INTERVAL C RLAST = RNEXT RNEXT = RNEXT + DRNOW IF (RNEXT .LT. REND) GO TO 170 IEND = 1 RNEXT = REND DRNOW = RNEXT - RLAST 170 RNEW = RLAST + 0.5D0 * DRNOW IF (KSTEP .GT. 1 .AND. IEND .NE. 1) THEN IF (TOLAI .LT. 1) THEN IF (DRNOW .LT. SPCMN) SPCMN = DRNOW END IF IF (DRNOW .GT. SPCMX) SPCMX = DRNOW END IF DRMID = RNEW - RNOW C C * RESTORE EIGENVALUES C CALL DCOPY (NCH, EIGNOW, 1, EIGOLD, 1) C * DEFINE LOCAL BASIS AT RNEW AND CARRY OUT TRANSFORMATIONS * TMAT IS USED AS SCRATCH MATRIX AND Y1 IS USED AS SCRATCH VECTOR HERE C CALL POTENT (W, VECNEW, TMAT, EIGNOW, HP, Y1, + RNEW, DRNOW, EN, XLARGE, NCH, $ P, MXLAM, VL, IV, RMLMDA, ERED, EINT, CENT, NPOTL) C C * DETERMINE MATRIX TO TRANSFORM LOG-DERIV MATRIX INTO NEW INTERVAL * SEE EQ. (22) OF M.H. ALEXANDER, "HYBRID QUANTUM SCATTERING ALGORITHMS C CALL DGEMUL(VECNEW, NCH, 'N', VECNOW, NCH, 'T', TMAT, NCH, 1 NCH, NCH, NCH) CALL DCOPY (NCH*NCH, VECNEW, 1, VECNOW, 1) C C * RESTORE RADIUS VALUES C RNOW = RNEW C C * DETERMINE APPROXIMATE VALUES FOR DIAGONAL AND OFF-DIAGONAL * CORRECTION TERMS C CALL CORR (EIGNOW, EIGOLD, HP, DRNOW, DRMID, XLARGE, CDIAG, : COFF, NCH) IF (ITWO .LT. 0) GO TO 200 IF (IEND .EQ. 1) RNOW = - RNOW C C * WRITE OR READ RELEVANT INFORMATION C 180 CALL OUTMAT (TMAT, EIGOLD, HP, ESHIFT, DRNOW, RNOW, : NCH, NCH, ITWO, ISCRU) IF (ITWO .EQ. 0) GO TO 200 C C * NEGATIVE RNOW IS CUE FOR LAST STEP IN SECOND ENERGY CALCULATION C IF (RNOW .GT. 0.D0) GO TO 200 RNOW = - RNOW IEND = 1 C C * GO BACK TO START NEW STEP C 200 CONTINUE C C * THE FOLLOWING STATEMENT IS REACHED ONLY IF THE INTEGRATION HAS * NOT REACHED THE ASYMPTOTIC REGION IN MAXSTP STEPS C WRITE (6,210) MAXSTP, RNEXT C IF(IPRINT.GT.40) WRITE(6,210) MAXSTP,RNEXT C 210 FORMAT (' *** AIRY PROPAGATION NOT FINISHED IN', I4, : ' STEPS: R-FIN SET TO', F8.4,' ***',/) XF = RNEXT 250 CONTINUE IF (ITWO .LT. 0) GO TO 260 CALL OUTMAT (VECNOW, EIGOLD, HP, ESHIFT, DRNOW, XF, NCH, : NCH, ITWO, ISCRU) C C * TRANSFORM LOG-DERIV MATRIX INTO FREE BASIS. TRANSFORMATION MATRIX IS * JUST VECNOW-TRANSPOSE; SEE EQ.(24) OF M.H. ALEXANDER, "HYBRID QUANTUM * SCATTERING ALGORITHMS ..." 260 CALL TRNSFM (VECNOW, Z, W, NCH, .FALSE., ISYM ) C IF (IPRINT.LT.41) GO TO 318 IF (ITWO .LT. 0) WRITE (6,280) IF (ITWO .EQ. 0) WRITE (6,290) IF (ITWO .GT. 0) WRITE (6,300) 280 FORMAT (' ** AIRY PROPAGATION - FIRST ENERGY;', : ' TRANSFORMATION MATRICES NOT WRITTEN') 290 FORMAT (' ** AIRY PROPAGATION - FIRST ENERGY;', : ' TRANSFORMATION MATRICES WRITTEN') 300 FORMAT (' ** AIRY PROPAGATION - SECOND ENERGY;', : ' TRANSFORMATION MATRICES READ') WRITE (6,305) RMIN, REND, TOLAI, NSTEP WRITE (6,310) SPCMN, SPCMX, POWR 305 FORMAT (' RBEGIN =', F7.3, ' REND =', F7.3, : ' TOLAI =', 1PE8.1, ' NINTERVAL =', I3) 310 FORMAT (' DR-MIN =', F7.3, ' DR-MAX =', F8.3, : ' POWER =', F4.1) C 318 CONTINUE IF(IPRINT.LT.35) GO TO 319 IF (ITWO .LT. 0) WRITE (6,280) IF (ITWO .EQ. 0) WRITE (6,290) IF (ITWO .GT. 0) WRITE (6,300) 319 CONTINUE C IF(IPRINT.LT.3) GO TO 320 WRITE (6, 315) RMIN, REND, SPCMN, SPCMX, NSTEP 315 FORMAT (' ** AIRY: RSTART =' ,F7.3,' REND =',F7.3, : ' DRMIN =',F7.3, ' DRMAX =',F7.3,' NSTEP =', I4) 320 CONTINUE RETURN END SUBROUTINE AIRYMP (X, FTHETA, FPHI, XMMOD, XNMOD) * SUBROUTINE TO RETURN THE MODULI AND PHASES OF THE AIRY FUNCTIONS AND * DERIVATIVES * AUTHOR: MILLARD ALEXANDER * CURRENT REVISION DATE: 23-SEPT-87 * ---------------------------------------------------------------------- * VARIABLES IN CALL LIST: * X ARGUMENT OF AIRY FUNCTIONS * FTHETA, XMMOD ON RETURN: CONTAIN THE (DOUBLE PRECISION) * PHASE AND MODULUS OF AI(X) AND BI(X) ( * BELOW). * FPHI, XNMOD ON RETURN: CONTAIN THE (DOUBLE PRECISION) * PHASE AND MODULUS OF AI'(X) AND BI'(X) * BELOW). * ---------------------------------------------------------------------- * FOR NEGATIVE X * ---------------------------------------------------------------------- * THE MODULI AND PHASES ARE DEFINED BY * AI(-X) = M(X) COS[THETA(X)] * BI(-X) = M(X) SIN[THETA(X)] * AI'(-X) = N(X) COS[PHI(X)] * BI'(-X) = N(X) SIN[PHI(X)] * IN OTHER WORDS * 2 2 2 * M(X) = SQRT[ AI(X) + BI(X) ] * 2 2 2 * N(X) = SQRT[ AI'(X) + BI'(X) ] * THETA(X) = ATAN [ BI(X) / AI(X) ] * PHI(X) = ATAN [ BI'(X) / AI'(X) ] * TO DETERMINE THESE MODULI AND PHASES WE USE THE SUBROUTINE * SCAIRY, WRITTEN BY D. MANOLOPOULOS (SEPT. 1986) * THIS SUBROUTINE RETURNS THE FOLLOWING QUANTITIES: * SCAI, SCBI, SCAIP, SCPIB, AND ZETA, WHERE C FOR X .LT. -5.0 C AI(X) = SCAI * COS(ZETA) + SCBI * SIN(ZETA) C BI(X) = SCBI * COS(ZETA) - SCAI * SIN(ZETA) C AI'(X) = SCAIP * COS(ZETA) + SCBIP * SIN(ZETA) C BI'(X) = SCBIP * COS(ZETA) - SCAIP * SIN(ZETA) C WHERE ZETA = (2/3) * (-X) ** (3/2) + PI/4 C C FOR -5.0 .LE. X .LE. 0.0 C C AI(X) = SCAI C BI(X) = SCBI C AI'(X) = SCAIP C BI'(X) = SCBIP C AND ZETA = 0 * ---------------------------------------------------------------------- * FOR POSITIVE X * ---------------------------------------------------------------------- * THE MODULI AND PHASES ARE DEFINED BY * AI(X) = M(X) SINH[THETA(X)] * BI(X) = M(X) COSH[THETA(X)] * AI'(X) = N(X) SINH[PHI(X)] * BI'(X) = N(X) COSH[PHI(X)] * IN OTHER WORDS * 2 2 2 * M(X) = SQRT[ BI(X) - AI(X) ] * 2 2 2 * N(X) = SQRT[ BI'(X) - AI'(X) ] * THETA(X) = ATANH [ AI(X) / BI(X) ] * PHI(X) = ATANH [ AI'(X) / BI'(X) ] * HERE THE THE EXPONENTIALLY SCALED AIRY FUNCTIONS * AI(X), AI'(X), BI(X), BI'(X) ARE: * AI(X) = AI(X) * EXP[ZETA] * AI'(X) = AI'(X) * EXP[ZETA] * BI(X) = BI(X) * EXP[-ZETA] * BI'(X) = BI'(X) * EXP[-ZETA] * TO DETERMINE THESE MODULI AND PHASES WE USE THE SUBROUTINE * SCAIRY, WRITTEN BY D. MANOLOPOULOS (SEPT. 1986) * THIS SUBROUTINE RETURNS THE FOLLOWING QUANTITIES: * SCAI, SCBI, SCAIP, SCPIB, AND ZETA * IN TERMS OF WHICH THE EXPONENTIALLY SCALED AIRY FUNCTIONS ARE DEFINED * AI(X) = SCAI * EXP(-ZETA) * BI(X) = SCBI * EXP(+ZETA) * AI'(X) = SCAIP * EXP(-ZETA) * BI'(X) = SCBIP * EXP(+ZETA) * WHERE ZETA = (2/3) * X ** (3/2) * * ---------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION X, FTHETA, FPHI, XMMOD, XNMOD, SCAI, : SCBI, SCAIP, SCBIP, ZETA, RATIO CALL SCAIRY (X, SCAI, SCBI, SCAIP, SCBIP, ZETA) IF ( X .LE. 0.D0) THEN XMMOD = SQRT( SCAI ** 2 + SCBI ** 2) XNMOD = SQRT( SCAIP ** 2 + SCBIP ** 2) FTHETA = ATAN2 (SCBI, SCAI) FPHI = ATAN2 (SCBIP, SCAIP) IF (X .LT. (-5.0D0) ) THEN FTHETA = FTHETA - ZETA FPHI = FPHI - ZETA END IF ELSE XMMOD = SQRT( - SCAI ** 2 + SCBI ** 2) XNMOD = SQRT( - SCAIP ** 2 + SCBIP ** 2) RATIO = SCAI / SCBI FTHETA = 0.5 * LOG ( (1.D0 + RATIO) / (1.D0 - RATIO) ) RATIO = SCAIP / SCBIP FPHI = 0.5 * LOG ( (1.D0 + RATIO) / (1.D0 - RATIO) ) END IF RETURN END SUBROUTINE ASROT(J,EVEC,H,EVAL,WKS,NH) IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL TD,EVLIST DIMENSION EVEC(NH,NH),H(NH,NH),EVAL(NH),WKS(NH) DIMENSION WT(2),ELEVEL(1000),JLEVEL(4000),ISYM(10),ISYM2(10), 1 ROTI(2) COMMON /CMBASE/ A(2),B(2),C(2),DJ,DJK,DK,DT,ROTI, 1 ELEVEL,EMAX,WT,SPNUC,NLEVEL,JLEVEL,JMIN,JMAX,JSTEP,ISYM,J2MIN, 1 J2MAX,J2STEP,ISYM2,JHALF,IDENT,MXJL,MXEL DATA EVLIST/.FALSE./ C C DO THE ACTUAL CALCULATION FOR A GIVEN J C ALPHA=0.5D0*(A(1)+B(1)) BETA=C(1)-ALPHA GAMMA=0.25D0*(A(1)-B(1)) TD = A(1).EQ.B(1) .AND. B(1).EQ.C(1) C JJ=J*(J+1) NK=J+J+1 DO 100 IR=1,NK KR=IR-J-1 DO 100 IC=1,IR KC=IC-J-1 TERM=0.D0 IF(KR.EQ.KC) THEN TERM=ALPHA*DBLE(JJ)+BETA*DBLE(KC*KC) 1 -DJ*DBLE(JJ*JJ)-DJK*DBLE(JJ*KC*KC)-DK*KC**4 IF(TD) TERM=TERM+0.5D0*DT*DBLE(-3*JJ*(JJ-2)+30*(JJ-2)*KC*KC 1 -35*KC*KC*(KC*KC-1)) ELSEIF(KR-KC.EQ.2) THEN KMID=(KR+KC)/2 TERM=GAMMA*SQRT(DBLE((JJ-KR*KMID)*(JJ-KC*KMID))) ELSEIF(KR-KC.EQ.4 .AND. TD) THEN TERM=1.25D0*DT*SQRT(DBLE((JJ-KC*(KC+1))*(JJ-(KC+1)*(KC+2))) 1 *DBLE((JJ-(KC+2)*(KC+3))*(JJ-(KC+3)*(KC+4)))) ENDIF H(IR,IC)=TERM 100 CONTINUE IFAIL=0 CALL F02ABF(H,NH,NK,EVAL,EVEC,NH,WKS,IFAIL) C WRITE(6,603) J,(EVAL(IC),IC=1,NK) 603 FORMAT('0 CALCULATED ROTATIONAL LEVELS FOR J =',I3,' ARE'/ 1 (8X,9F12.5)) C C IF THE RAW EIGENVECTORS ARE DEGENERATE, THEY MAY NOT HAVE C PROPER SYMMETRY. SEEK DEGENERATE SETS AND FORCE SYMMETRY ON THEM. C ALSO PRINT SPHERICAL TOP SYMMETRY LABELS IF ANY DEGENERATE SETS C ARE PRESENT. C CALL DMSYM(J,NK,EVAL,EVEC,H,WKS) C IF(EVLIST) THEN WRITE(6,604) 604 FORMAT('0 EIGENVECTOR COEFFICIENTS:') DO 200 IR=1,NK KR=IR-J-1 WRITE(6,605) J,KR,(EVEC(IR,IC),IC=1,NK) 605 FORMAT(2I4,9F12.8/(8X,9F12.8)) 200 CONTINUE ENDIF RETURN END SUBROUTINE AXSCAT(N, NSQ, MXLAM, NPOTL, 1 SR, SI, U, VL, IV, EINT, CENT, WVEC, L, NB, 2 P, Y1, Y2, Y3, Y4, VECNOW, VECNEW, EIGOLD, EIGNOW, HP, 3 ICODE, IPRINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C --------------------------------------------------------------- C THIS IS TIM PHILLIP'S INTERFACE OF ALEXANDER HIBRIDON C CODES TO MOLSCAT: SCATTERING CALC USING DAPROP AND THEN C AIRPRP. ON EXIT SR AND SI CONTAIN THE S-MATRIX. C SR IS USED INTERNALLY TO HOLD THE LOG DERIVATIVE MATRIX C ICODE.EQ.2 FOR SUBSEQUENT ENERGIES. C --------------------------------------------------------------- C REORGANIZED BY SG (2/2/93): CORRECTS NSTEPS=0 PROBLEM, BUT C ALSO CALCULATES SOMEWHAT DIFFERENT STEP SIZES FROM EARLIER CODE. C --------------------------------------------------------------- C C DIMENSION STATEMENTS FOR ARGUMENT LIST C DIMENSION U(NSQ),Y1(N),Y2(N),Y3(N),Y4(N) DIMENSION P(MXLAM),VL(2),IV(2),SR(NSQ),SI(NSQ), & EINT(N),CENT(N),WVEC(N),L(N),NB(N) DIMENSION VECNOW(NSQ),VECNEW(NSQ),EIGOLD(N),EIGNOW(N),HP(N) C LOGICAL IREAD,IWRITE, LLD,LAIRY C C COMMON BLOCKS TO COMMUNICATE WITH PROPAGATORS C THE FOLLOWING VARIABLES FROM COMMON/DRIVE/ ARE USED WITH THIS C PROPAGATOR: STEPS,RMIN,RMAX,ERED,RMLMDA,NOPEN,ISCRU,TOLHI,RMID, C AND SOMETIMES DR C COMMON/DRIVE/STEST,STEPS,STABIL,CONV,RMIN,RMAX,XEPS,DR, 1 DRMAX,RMID,TOLHI,RTURN,VTOL,ESHIFT,ERED,RMLMDA, 2 NOPEN,JKEEP,ISCRU,MAXSTP C COMMON/HIBRIN/POWRX,DRAIRY,IABSDR C C SET UP TO USE UNIT (ISCRU) IREAD = ICODE.EQ.2 .AND. ISCRU.GT.0 IWRITE = ICODE.EQ.1 .AND. ISCRU.GT.0 C --------------------------------------------------------------- IF((.NOT.IREAD) .AND. IWRITE) THEN ITWO = 0 ELSE IF(IREAD .AND. (.NOT.IWRITE)) THEN ITWO = 1 ELSE IF((.NOT.IREAD) .AND. (.NOT.IWRITE)) THEN ITWO = -1 ELSE WRITE(6,*) ' ILLEGAL IREAD/IWRITE COMBINATION ' WRITE(6,*) ' BOTH SIMULTANEOUSLY TRUE ' STOP END IF C-------------------------------------------------------------------- C C DECIDE WHICH CALCULATIONS TO DO. C ON THE ASSUMPTION THAT RMIN.LT.RMAX THERE ARE THREE CASES C 1) RMID.LE.RMIN.LT.RMAX C 2) RMIN.LT.RMID.LT.RMAX C 3) RMIN.LT.RMAX.LE.RMID C CODE BELOW SETS FOLLOWING LLD LAIRY RSWTCH C CASE 1 F T RMIN C CASE 2 T T RMID C CASE 3 T F RMAX C INTEGRATION RANGES ARE THEN DAPROP: RMIN -> RSWTCH C AIRY: RSWTCH -> RMAX C-------------------------------------------------------------------- RBEGIN=RMIN REND=RMAX LLD=RBEGIN.LT.RMID LAIRY=RMID.LT.REND RSWTCH=MIN(REND,RMID) RSWTCH=MAX(RSWTCH,RBEGIN) C C CALCULATE WAVEVECTORS AND STEP SIZE WMAX=0.D0 NOPEN=0 DO 20 I=1,N DIF=ERED-EINT(I) WVEC(I)=SIGN(SQRT(ABS(DIF)),DIF) WMAX=MAX(WMAX,WVEC(I)) NB(I)=I IF (DIF.GT.0.D0) NOPEN=NOPEN+1 20 CONTINUE IF (NOPEN.EQ.0) RETURN C IF (IREAD) GO TO 40 PI=ACOS(-1.D0) DRLD=PI/(WMAX*STEPS) IF (IABSDR .EQ. 1 .AND. DR .GT. 0.D0) DRLD=DR NSTEPS=(RSWTCH-RBEGIN)/DRLD IF (IWRITE) WRITE (ISCRU) RBEGIN,RSWTCH,REND,DRLD,NSTEPS GO TO 60 40 READ (ISCRU) RBEGIN,RSWTCH,REND,DRLD,NSTEPS 60 CONTINUE C C SET REND FOR YTOK, AND RESET RSWTCH IN CASE WE DON'T CALL AIRPRP RYON=REND REND=RSWTCH RSWTCH=RBEGIN C LLD=LLD .AND. NSTEPS.GT.0 IF (LLD) THEN RSWTCH=REND C PROPAGATE LOG DERIVATIVE THROUGH FIRST SEGMENT. C ISTART=0 REQUESTS INITIALIZATION OF LOG-DERIVATIVE MATRIX ISTART=0 CALL DAPROP(U, SR, N, 1 RBEGIN, REND, NSTEPS, IREAD, IWRITE, ISCRU, 2 Y1, Y2, Y3, Y4, 3 P, VL, IV, ERED, EINT, CENT, RMLMDA, 4 MXLAM, NPOTL, ISTART, NODES) C ------------------------------------------------------------- IF (IPRINT.GE.3) WRITE (6,1000) RBEGIN,RSWTCH,NSTEPS 1000 FORMAT(' AXSCAT. LOG DERIVATIVE MATRIX INTEGRATED FROM ', & F12.4,' TO ',F12.4,' IN ',I6,' STEPS.') ELSE C INITIALIZE LOG-DERIVATIVE MATRIX IF DAPROP NO CALLED DO 42 I=1,NSQ 42 SR(I)=0.D0 DO 43 I=1,NSQ,N+1 43 SR(I)=1.D30 IF (IPRINT.GE.3) WRITE (6,1010) 1010 FORMAT(' AXSCAT. DAPROP NOT CALLED: LOG DERIVATIVE MATRIX ', & 'INITIALIZED.') ENDIF C C USE AIRY PROPAGATOR FOR THE REMAINDER OF THE SCATTERING REGION C IF (.NOT.LAIRY) GO TO 41 REND=RYON DRA=DRLD IF (DRAIRY .GT. 0.D0) DRA = DRAIRY CALL AIRPRP(SR,U,SI,VECNOW,VECNEW,EIGOLD,EIGNOW,HP, 1 Y1,Y2,Y3,Y4,RSWTCH, 2 REND,DRA,ERED,TOLHI,POWRX,ESHIFT,N, 3 ITWO,IREAD,IWRITE,IPRINT,ISCRU,P,MXLAM,VL,IV,RMLMDA, 4 EINT,CENT,NPOTL) C C SORT CHANNELS BY ASYMPTOTIC ENERGY C 41 CONTINUE IF (N.EQ.1) GO TO 100 NM1=N-1 DO 80 I=1,NM1 IP1=I+1 DO 80 J=IP1,N IF (EINT(NB(I)).LE.EINT(NB(J))) GO TO 80 IT=NB(I) NB(I)=NB(J) NB(J)=IT 80 CONTINUE C C CALCULATE K AND S MATRICES C 100 CALL YTOK(NB,WVEC,L,N,NOPEN,Y1,Y2,Y3,Y4,SR,SI,U,REND) CALL KTOS(U,SR,SI,NOPEN) RETURN END SUBROUTINE BAS9IN(PRTP,IBOUND) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE CHARACTER*8 PRTP(4),QNAME(10) LOGICAL LEVIN,EIN,LCNT DIMENSION ROTI(12),ELEVEL(1000),JLEVEL(4000), 1 ISYM(10),ISYM2(10),WT(2) DIMENSION JLEV(1),VL(1),IV(1),CENT(1),J(1),L(1),LAM(1) COMMON/CMBASE/ROTI,ELEVEL,EMAX,WT,SPNUC,NLEVEL,JLEVEL, 1 JMIN,JMAX,JSTEP,ISYM,J2MIN,J2MAX,J2STEP,ISYM2,JHALF,IDENT 2 ,MXJL,MXEL C C BAS9IN IS CALLED ONCE FOR EACH SCATTERING SYSTEM (USUALLY ONCE C PER RUN) AND CAN READ IN ANY BASIS SET INFORMATION NOT CONTAINED C IN NAMELIST BLOCK &BASIS. IT MUST ALSO HANDLE THE FOLLOWING C VARIABLES AND ARRAYS: C C PRTP SHOULD BE RETURNED AS A CHARACTER STRING DESCRIBING THE C COLLISION TYPE C IDENT CAN BE SET>0 IF A COLLISION OF IDENTICAL PARTICLES IS C BEING CONSIDERED AND SYMMETRISATION IS REQUIRED. C HOWEVER, THIS WOULD REQUIRE EXTRA CODING IN IDPART. C IBOUND CAN BE SET>0 IF THE CENTRIFUGAL POTENTIAL IS NOT OF THE C FORM L(L+1)/R**2; IF IBOUND>0, THE CENT ARRAY MUST BE C RETURNED FROM ENTRY CPL9 C IBOUND=1 PRTP(1)=' BODY-F' PRTP(2)='IXED ATO' PRTP(3)='M-DIATOM' PRTP(4)=' ' RETURN C ENTRY SET9(LEVIN,EIN,NLEV,JLEV,NQN,QNAME,MXPAR,NLABV) C C SET9 IS CALLED ONCE FOR EACH SCATTERING SYSTEM. IT SETS UP: C MXPAR, THE NUMBER OF DIFFERENT SYMMETRY TYPES ("PARITY CASES") C NLEVEL AND JLEVEL, UNLESS LEVIN IS .TRUE.; C JLEV AND NLEV; C ELEVEL, UNLESS EIN IS .TRUE. C IF THE LOGICAL VARIABLES ARE .TRUE. ON ENTRY, THE CORRESPONDING C QUANTITIES WERE INPUT EXPLICITLY IN NAMELIST BLOCK &BASIS. C IF EIN IS .FALSE., THE MOLECULAR CONSTANTS MUST HAVE BEEN SUPPLIED C IN THE &BASIS ARRAY ROTI: THE PROGRAMMER MAY USE THESE IN ANY WAY C HE LIKES, BUT SHOULD OUTPUT THEM HERE FOR CHECKING. C NOTE THAT JLEVEL CONTAINS JUST THE QUANTUM NUMBERS NECESSARY TO C SPECIFY THE THRESHOLD ENERGY (AND ELEVEL CONTAINS THE CORRESPONDING C ENERGIES WHEREAS JLEV CONTAINS ALL THE CHANNEL QUANTUM NUMBERS EXCEPT C THE ORBITAL L, WHICH MAY BE A SUPERSET. THE LAST COLUMN OF THE JLEV C ARRAY CONTAINS A POINTER TO THE ENERGY IN THE ELEVEL ARRAY. C MXPAR=2 NLABV=1 IF(LEVIN) GOTO 220 NLEVEL=0 NLEV=0 DO 210 I=JMIN,JMAX,JSTEP NLEVEL=NLEVEL+1 JLEVEL(NLEVEL)=I C NL IS NUMBER OF SETS OF INTERNAL QUANTUM NUMBERS FOR THIS LEVEL NL=1+MIN(J2MAX,I) NLEV=NLEV+NL 210 CONTINUE GOTO 230 220 WRITE(6,602) 602 FORMAT('0 BASIS FUNCTIONS TAKEN FROM &BASIS (JLEVEL) INPUT') C C IF NLEV AND NLEVEL ARE DIFFERENT, IT MAY BE NECESSARY TO BUILD UP JLEV C IN A DIFFERENT ORDER AND REARRANGE IT LATER - SEE SET3 CODING IN SETBAS 230 NQN=3 QNAME(1)=' J ' QNAME(2)=' |K| ' C LOOP OVER LEVELS AGAIN, THIS TIME SETTING UP JLEV II=0 IJ=0 DO 250 I=JMIN,JMAX,JSTEP II=II+1 DO 250 K=0,MIN(J2MAX,I) IJ=IJ+1 JLEV(IJ)=I JLEV(NLEV+IJ)=K JLEV(NLEV*(NQN-1)+IJ)=II 250 CONTINUE C IF(EIN) GOTO 280 WRITE(6,604) ROTI(1) 604 FORMAT('0 ENERGY LEVELS CALCULATED FROM B =',F10.5) C DO 270 I=1,NLEVEL JI=JLEVEL(I) ELEVEL(I)=ROTI(1)*DBLE(JI*(JI+1)) 270 CONTINUE RETURN C 280 WRITE(6,605) 605 FORMAT('0 ENERGY LEVELS TAKEN FROM &BASIS (ELEVEL) INPUT') RETURN C ENTRY BASE9(LCNT,N,JTOT,ICODE,JLEV,NLEV,NQN,J,L) C C BASE9 IS CALLED EITHER TO COUNT THE ACTUAL NUMBER OF CHANNEL BASIS C FUNCTIONS OR ACTUALLY TO SET THEM UP (IN THE J AND L ARRAYS). C IT IS CALLED FOR EACH TOTAL J (JTOT) AND PARITY CASE (ICODE). C IF LCNT IS .TRUE. ON ENTRY, JUST COUNT THE BASIS FUNCTIONS. OTHERWISE, SET C UP J (POINTER TO JLEV) AND L (ORBITAL ANGULAR MOMENTUM) FOR EACH CHANNEL. C THIS MUST TAKE INTO ACCOUNT JTOT AND ICODE. C ONE IMPORTANT OVERALL EFFECT IS THAT THE THRESHOLD ENERGY IS IN C ELEVEL(JLEV(NLEV*(NQN-1)+J(I)). CHECK THIS! C N=0 DO 320 I=1,NLEV K=JLEV(NLEV+I) IF(K.GT.JTOT) GOTO 320 IF(K.EQ.0 .AND. ICODE.EQ.1) GOTO 320 N=N+1 IF(LCNT) GOTO 310 J(N)=I L(N)=JTOT 310 CONTINUE 320 CONTINUE RETURN C ENTRY CPL9(N,ICODE,NPOTL,LAM,MXLAM,NLEV,JLEV,J,L,JTOT, 1 VL,IV,CENT,IBOUND,IEXCH,IPRINT) C C CPL9 IS CALLED AFTER BASE9 FOR EACH JTOT AND ICODE, TO SET UP THE C POTENTIAL COUPLING COEFFICIENTS VL. C IF IBOUND>0, IT ALSO SETS UP THE CENTRIFUGAL COEFFICIENTS CENT. C INDICES SPECIFYING THE MXLAM DIFFERENT POTENTIAL SYMMETRIES ARE IN C THE FIRST XX*MXLAM ELEMENTS OF LAM; THE STRUCTURE OF THE LAM ARRAY C (AND THE VALUE OF XX) IS CHOSEN BY THE PROGRAMMER, AND MUST C CORRESPOND WITH THAT USED IN SUBROUTINE POTENL. C NPOTL IS THE NUMBER OF DIFFERENT POTENTIAL TERMS WHICH CONTRIBUTE TO C EACH MATRIX ELEMENT (SEE SUBROUTINE WAVVEC). IT SOMETIMES SAVES C A SIGNIFICANT AMOUNT OF SPACE IF IT CAN BE LESS THAN MXLAM. C ROOT2=SQRT(2.D0) NPOTL=MXLAM DO 550 LL=1,MXLAM LM=LAM(LL) NNZ=0 I=LL DO 540 ICOL=1,N JC=JLEV(J(ICOL)) KC=JLEV(J(ICOL)+NLEV) DO 540 IROW=1,ICOL JR=JLEV(J(IROW)) KR=JLEV(J(IROW)+NLEV) VL(I)=0.D0 IV(I)=LL IF(KR.NE.KC) GOTO 510 ISUM=LM+JR+JC IF(ISUM-2*(ISUM/2).NE.0 .OR. LM.GT.JC+JR) GOTO 540 VL(I)=PARITY(KC)*SQRT(DBLE(JC*(JC+1)*JR*(JR+1)))*THREEJ(JC,LM,JR) 1 *THRJ(DBLE(JC),DBLE(LM),DBLE(JR),DBLE(KC),0.D0,DBLE(-KC)) IF(VL(I).NE.0.D0) NNZ=NNZ+1 GOTO 540 510 IF(LM.GE.0 .OR. JR.NE.JC .OR. IABS(KR-KC).GT.1) GOTO 540 VL(I)=SQRT(DBLE((JC*(JC+1)-KR*KC)*(JTOT*(JTOT+1)-KR*KC))) IF(JC.EQ.0 .NEQV. JR.EQ.0) VL(I)=ROOT2*VL(I) IF(VL(I).NE.0.D0) NNZ=NNZ+1 540 I=I+NPOTL IF(NNZ.EQ.0) WRITE(6,612) JTOT,LL 612 FORMAT(' * * * NOTE. FOR JTOT =',I4,', ALL COUPLING', 1 ' COEFFICIENTS ARE 0.0 FOR POTENTIAL SYMMETRY',I4) 550 CONTINUE C C NOW THE CENTRIFUGAL POTENTIAL DO 570 I=1,N JC=JLEV(J(I)) KC=JLEV(J(I)+NLEV) CENT(I)=DBLE(JTOT*(JTOT+1)+JC*(JC+1)-2*KC*KC) 570 CONTINUE RETURN C ENTRY DEGEN9(J1,J2,RESULT) C C DEGEN9 IS CALLED TO OBTAIN THE DEGENERACY FACTOR FOR THE DENOMINATOR C OF A CROSS-SECTION CALCULATION; IT DOES NOT MATTER FOR BOUND STATES. C C JI=JLEVEL(J1) C RESULT=DBLE(2*JI+1) RETURN END SUBROUTINE BASE (JTOT, JLEV, N, J, L, CINT, EINT, CENT, VL, IV, & MXLAM, NPOTL, LAM, WVEC, WGHT, IEXCH, THETA, PHI, & ICODE, LCNT, ERED, NLEVV, PRINT) C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C . CHANGED APR 1986 TO COMBINE BASIN & IOSBIN HANDLING. C . CHANGED JAN 1988 TO ALLOW USER-DEFINED BASIS (ITYPE=9) C . CHANGED MAR 1993 TO SET MPLMIN=TRUE FOR ITYPE=23,IDENT=1 C . CHANGED NOV 1993 TO USE IXNEXT DIRECTLY IN PLACE OF IC. C . IC REMOVED FROM ARGUMENT LIST OF BASIN AND SET6. C . CHANGED JAN 1994 TO USE IV() INDEXING FOR ITYP=10*N + 2 C . APR 1994 TO INCLUDE IEXCH IN PARAMETER LIST C NEW ORDERING OF ITYPE=23 'PARITY CASES' C . CHANGED JUL 1994 FOR VERSION 14 C . AUG 1994 TO INTEGRATE ITYPE = 4 CODE TO VERSION 14 C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE C BELOW IS BEGINNING OF A LIMITED SAVE LIST C SAVE ITP,ITYPE,IOFF,IBOUND,JZCSFL,EMAXK,WTM,IGO,IGODG,MJMX C IDENT,WT,JMAX,JMIN -- IN CMBASE C NLEV,NQN,MPLMIN -- IN PRBASE C C BASE HANDLES QUANTUM BASIS SETS FOR SCATTERING CALCULATION. C C ITYPE DESCRIBES TARGET-PROJECTILE TYPES. C CURRENT IMPLEMENTATION FOR C ITYPE=1 LINEAR RIGID ROTOR HIT BY AN ATOM C ITYPE=2 DIATOMIC VIB-ROTOR HIT BY AN ATOM C ITYPE=3 LINEAR RIGID ROTOR - LINEAR RIGID ROTOR C ITYPE=4 RIGID ASYMMETRIC TOP HIT BY LINEAR RIGID ROTOR C ITYPE=5 NEAR-SYMMETRIC TOP RIGID ROTOR HIT BY AN ATOM C ITYPE=6 ASYMMETRIC TOP RIGID ROTOR HIT BY AN ATOM. C ITYPE=7 DIATOMIC VIB-ROTOR HIT BY AN ATOM, WHERE A FULL C SET OF EXPECTATION VALUES OF THE INTERMOLECULAR C POTENTIAL BETWEEN (V,J) AND (V',J') DIATOM INTERNAL C STATES IS SUPPLIED C ITYPE=8 ATOM-SURFACE SCATTERING C C ITYPE=ITYPE+10 FOR EFFECTIVE POTENTIAL METHOD OF RABITZ. C ITYPE=ITYPE+20 FOR COUPLED STATES OF MCGUIRE-KOURI. C ITYPE=ITYPE+30 FOR DLD METHOD OF DEPRISTO AND ALEXANDER. C C C ENTRY BASIN C READS AND PROCESSES &BASIS DATA TO DESCRIBE ASYMPTOTIC LEVELS. C QUANTUM NOS. AND INDEXING ARE IN JLEVEL(NLEVEL) AND C JLEV(NLEV,NQN), IN DIFFERENT FORMATS. C ASYMPTOTIC ENERGIES ARE IN ELEVEL(NLEVEL). C C MAIN ENTRY BASE C SETS UP BASIS FOR EACH PARTIAL CALCULATION FROM ASYMPTOTIC C LEVELS (STORED IN JLEV) COUPLED WITH COLLISION ORBITAL ANGULAR C MOMENTUM. C LCNT=.TRUE. MEANS ONLY COUNT NUMBER OF CHANNELS C LCNT=.FALSE. MEANS SET UP BASIS FUNCTIONS IN ALLOCATED STORAGE. C ICODE (=1...MXPAR) IS AN INDEX FOR THE CURRENT SYMMETRY BLOCK. C IPAR AND IEXCH SUBDIVIDE ICODE=1,MXPAR INTO C 1) PARITY, IPAR=0 (EVEN), 1 (ODD) C 2) EXCHANGE SYM., IEXCH=0 (NO EXCHANGE), 1 (ODD), 2 (EVEN). C IT IS NECESSARY TO SET FOLLOWING -- C ASYMPTOTIC LEVEL IN J, ORBITAL ANGULAR MOMENTUM IN L, C ASYMPTOTIC ENERGY IN EINT, CENTRIFUGAL ENERGY IN CENT, C AND COUPLING MATRIX ELEMENTS IN VL. C C EXTRA FEATURE ADDED AT UNIV. OF WATERLOO MAY 82. ARRAY IV C IS AN INDEX ARRAY SUCH THAT VL(I) IS A COEFFICIENT FOR C TERM NUMBER IV(I) IN THE POTENTIAL ARRAY RETURNED BY C SUBROUTINE POTENL. THE INTRODUCTION OF THIS ARRAY ONLY C MAKES ANY REAL DIFFERENCE FOR ITYPE=10*N+7, FOR WHICH IT C ENABLES LARGE ECONOMIES IN STORAGE FOR THE VL ARRAY. C ** 1/27/93 IV ARRAY USED IF AND ONLY IF IVLFL.GT.0 ** C NPOTL IS THE NUMBER OF "CHUNKS" OF SIZE N*(N+1)/2 WHICH C COMPRISE VL AND IV. NPOTL=MXLAM EXCEPT FOR ITYPE=10*N+7 AND 8. C FOR ITYPE=10*N+7, NPOTL IS EQUAL TO K+1, WHERE K IS THE C INDEX OF THE HIGHEST ORDER LEGENDRE POLYNOMIAL ACTUALLY C PRESENT IN THE POTENTIAL. C FOR ITYPE=8, NPOTL=1. C C ADDITIONAL CHANGE MADE MAY 82. THE VL ARRAY IS NOW STORED C SO THAT THE POTENTIAL SYMMETRY TERM IS MOST RAPIDLY VARYING, C RATHER THAN THE CHANNEL INDICES AS BEFORE. THIS IS TO KEEP C PAGING TO A MINIMUM IN SUBROUTINE WAVMAT. C C ENTRY DEGEN PROVIDES DEGENERACY INFORMATION FOR USE IN OUTPUT. C DIMENSION CENT(NLEVV),EINT(NLEVV),WVEC(NLEVV),VL(NLEVV),IV(NLEVV) LOGICAL LCNT,EIN,LEVIN,MPLMIN,LSIG INTEGER J(NLEVV),L(NLEVV),LAM(MXLAM),JLEV(NLEVV) INTEGER PRINT,EUNITS CHARACTER*8 QNAME(10),QTYPE(10),PRTP(4,9),PTP(2) DIMENSION WTM(2),IOSNGP(3) C C DYNAMIC STORAGE COMMON BLOCK ... COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1) C C COMMON BLOCK FOR BASIS DATA C 2 AUG 94 V14 VERSION OF CMBASE; ISYM NO LONGER EQUIV J2MAX C DIMENSIONS OF JLEVEL,ELEVEL SET HERE AND HELD IN /CMBASE/MXJL,MXEL PARAMETER (MXJLVL=4000,MXELVL=1000) DIMENSION ROTI(12),ALPHAE(2),BE(2),DE(2),A(2),B(2),C(2),WE(2), 1 WEXE(2),WT(2),ELEVEL(MXELVL),EEE(1016) DIMENSION JLEVEL(MXJLVL),ISYM(10),ISYM2(10) EQUIVALENCE (ROTI(1),BE(1),A(1)), (ROTI(3),ALPHAE(1),B(1)), 1 (ROTI(5),DE(1),C(1)),(JMIN,J1MIN),(JMAX,J1MAX),(JSTEP,J1STEP), 2 (ROTI(1),EEE(1)), (ROTI(7),WE(1),DJ ), (ROTI(8),DJK), 3 (ROTI(9),WEXE(1),DK), (J2MAX,KSET) COMMON /CMBASE/ ROTI,ELEVEL,EMAX,WT,SPNUC, NLEVEL,JLEVEL,JMIN, 1 JMAX,JSTEP,ISYM,J2MIN,J2MAX,J2STEP,ISYM2,JHALF,IDENT,MXJL,MXEL C COMMON /PRBASE/ ITYPX,NQN,NLEV,MVALUE,IPTY,MPLMIN C COMMON/VLSAVE/IVLU C C ARRAYS FOR NAMELIST &BASIS C CHARACTER*6 BNAMES(39) C DIMENSION LOCN(39),INDX(39) C C V14: ISYM2, JHALF ADDED TO NAMELIST NAMELIST /BASIS/ ROTI,JMIN,JMAX,JSTEP,ITYPE 1 ,NLEVEL,JLEVEL,ELEVEL,EMAX,EMAXK,BE,ALPHAE,DE,A,B,C,WE,WEXE 2 ,J1MAX,J1MIN,J2MAX,J2MIN,J1STEP,J2STEP 3 ,WT,IDENT,SPNUC,EUNITS,IASYMU,JZCSMX,IBOUND,JZCSFL 4 ,IOSNGP,IPHIFL,ISYM,ISYM2,KSET,IVLU,JHALF C DATA QTYPE/ ' J ', ' K ', 1 ' PRTY ', ' J1 ', ' J2 ', ' J12 ', 2 ' V ', ' TAU ',' ','SIG INDX'/ DATA PRTP/' LINEAR',' RIGID R','OTOR - ',' ATOM. ', 1 ' DIATOM','IC VIB-R','OTOR - ',' ATOM. ', 2 ' LINEAR',' ROTOR -',' LINEAR ','ROTOR. ', 3 ' ASYMME','TRIC TOP',' - LINEA','R ROTOR.', 4 'ATOM - N','EAR SYM.',' TOP RIG','ID ROTOR', 5 ' ASYMME','TRIC TOP',' - ATOM ',' ', 6 4*' ', 7 ' ATOM -',' CORRUGA','TED SURF','ACE ', 8 4*' '/ DATA PTP/', ODD ',', EVEN '/ C C C DATA BNAMES/'ROTI','JMIN','JMAX','JSTEP','ITYPE', C 1 'NLEVEL','JLEVEL','ELEVEL','EMAX','EMAXK', C 2 'BE','ALPHAE','DE','A','B','C','WE','WEXE', C 2 'J1MAX','J1MIN','J2MAX','J2MIN','J1STEP','J2STEP', C 3 'WT','IDENT','SPNUC','EUNITS','IASYMU','JZCSMX','IBOUND', C 4 'JZCSFL','IOSNGP','IPHIFL','ISYM','KSET','IVLU','ISYM2', C 5 'JHALF'/ C DATA INDX/39*0/ C C SET UP BASIS FUNCTIONS C IEXCH=0 IPAR=0 C IF (ITYPE.EQ.8) GO TO 5208 IF (ITP.EQ.9) GO TO 5209 IF (ITYPE.LE.10) GO TO 5201 IF (ITYPE.LE.20) GO TO 5202 IF (ITYPE.LE.30) GO TO 5203 C C CODE FOR DECOUPLED L-DOMINANT APPROX OF ALEXANDER C 5204 N=0 LAMBDA=ICODE-1 DO 4001 I=1,NLEV JI=JLEV(I) LI=JTOT+LAMBDA-JI IF (LI.LT.IABS(JTOT-JI) .OR. LI.GT.(JTOT+JI) ) GO TO 4001 N=N+1 IF (LCNT) GO TO 4001 J(N)=I L(N)=LI 4001 CONTINUE IF (LCNT) GO TO 5000 GO TO 8000 C C CODE BELOW FOR MCGUIRE-KOURI J-Z CONSERVING COUPLED STATES APPROX. C 5203 N=0 C>SG CODE BELOW IS REVISED APR 94 -- NEW ORDERING OF PARITY CASES IF (MPLMIN) THEN IF (IDENT.EQ.0) THEN MVALUE=ICODE-1 IF (ITYPE.EQ.25.AND.ISYM(1).NE.-1) THEN IBLOCK=1+MVALUE/(MJMX+1) MVALUE=MOD(MVALUE,MJMX+1) KREQ=IBLOCK-1 ENDIF ELSE IEXCH=2-MOD(ICODE,2) MVALUE=(ICODE+1)/2-1 IF (WT(IEXCH).EQ.0.D0) THEN IF (PRINT.GE.3) WRITE(6,690) JTOT,ICODE,PTP(IEXCH) GO TO 5000 ENDIF ENDIF ELSE C CODE BELOW IS FOR MPLMIN=.FALSE. (NOT USED, BUT COULD BE REVIVED) IF (IDENT.EQ.0) THEN WRITE(6,*) ' *** BASE (APR 94). MPLMIN=.FALSE. .AND. IDENT.' 1 ,'EQ.0 ARE NOT ALLOWED' STOP ELSE ICD=(ICODE+1)/2 MVALUE=ICD/2 IF (ICD-2*(ICD/2).EQ.0) MVALUE=-MVALUE ENDIF ENDIF C SET IPAR (=1 FOR MVALUE=0, =2 OTHERWISE) IPAR=1 IF (MVALUE.NE.0) IPAR=2 C