C THIS ROUTINE IS THE MAIN PROGRAM FOR MOLSCAT VERSION 12, 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 12 - JUN 93 ----- 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 11 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 WAVEFUNCTIONS 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 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,RTURNM,LINE DIMENSIONS LIMITED BY VALUES ... PARAMETER(MXNRG=100,MXLN=200,MXTEMP=5,MXRTM=100) 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 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 ARRAY TO SAVE TURNING POINTS FROM DIFFERENT PARITY CASES C FOR IRMSET > 0 OPTION. DIMENSION RTURNM(MXRTM) 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 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(87),LOCN(87),INDX(87) 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'/ C DATA INDX/87*0/ C C DATA LABEL/10*' '/ DATA CWD/' ','(8-BYTE)'/ DATA CTIME/' '/,CDATE/' '/ DATA IPROGM/12/, PDATE/'(NOV 93)'/ 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 DATA NLABV/1,3,3,-1,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 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,' (1993)'/ 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 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 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 C CALL NAMLIS('&INPUT',INAMES,LOCN,INDX,87,IEOF) C IF(IEOF.EQ.1) GOTO 1040 C-------------------------------------------------------------- READ(5,INPUT,END=1040) C WRITE(6,120) 120 FORMAT(/'0 /INPUT/ DATA ARE --') WRITE(LABL,'(A80)') LABEL WRITE(6,130) LABL 130 FORMAT('0 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('0 SCRATCH CORE STORAGE ALLOCATION IS',I10,A8, 1 ' WORDS (',F8.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('0 INTEGRATOR REQUESTED BY INPUT VALUE INTFLG =',I3) 220 FORMAT('0***** ERROR - NO IMPLEMENTATION FOR THIS INTFLG' 1 ,' - RUN HALTED.') 240 FORMAT('0 COUPLED EQUATIONS SOLVED BY METHOD OF DEVOGELAERE.') 250 FORMAT('0 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('0 COUPLED EQUATIONS SOLVED BY WALKER-LIGHT R-MATRIX', 1 ' PROPAGATOR ALGORITHM'/'0 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('0',' RVFAC =',F7.2,' OVERRIDES INPUT RMID') 300 FORMAT('0 COUPLED EQUATIONS SOLVED BY LOG DERIVATIVE METHOD ', 1 'OF JOHNSON') 310 FORMAT('0 INTEGRATION PARAMETERS ARE RMIN =',F7.2,8X, 1 'STEPS = ',F7.1/33X,'RMAX =',F7.2) 320 FORMAT('0 CHANGING TO VARIABLE INTERVAL / VARIABLE STEP METHOD', 1 ' AT LONG RANGE'/'0 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('0 COUPLED EQUATIONS SOLVED BY DIABATIC ', 1 'MODIFIED LOG DERIVATIVE METHOD OF MANOLOPOULOS') 350 FORMAT('0 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('0 AIRY PARAMETERS ','RMID =',F10.4/ 2 33X,'DRAIRY=',F10.4/33X,'TOLHI=',F13.6/ 3 33X,'POWRX =',F8.2) 355 FORMAT('0 DRAIRY.LT.0 TAKES INITIAL AIRY STEP SIZE FROM' 1 ,' MODIFIED LOG-DERIVATIVE VALUE.') 356 FORMAT('0 TOLHI.GE.1 -- AIRY STEP SIZE INCREASED BY' 1 ,' FACTOR OF TOLHI AT EACH STEP') 357 FORMAT('0 TOLHI.LT.1 -- AIRY STEPS ADJUSTED TO MAINTAIN' 1 ,' APPROX. ACCURACY VIA PERTURBATION THEORY AND POWRX.') 370 FORMAT('0 EQUATIONS SOLVED BY WKB APPROXIMATION WITH GAUSS-' 1 ,'MEHLER INTEGRATION. SEE R. T PACK, JCP 60, 633 (1974).'/ 2 '0 NOTE THAT THIS IS IMPLEMENTED ONLY FOR ONE CHANNEL', 3 ' CASES, E.G., IOS CALCULATIONS.'/ 4 '0 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('0 IRXSET =',I3,' OPTION. RMAX ADJUSTED AUTOMATICALLY ', 1 'FOR EACH NEW JTOT,MVAL') IF(IRMSET.LE.0) GOTO 420 WRITE(6,390) IRMSET 390 FORMAT('0 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('0***** 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('0 ENERGY-INDEPENDENT MATRICES SAVED FROM A ', 1 'PREVIOUS RUN WILL BE READ FROM UNIT',I3) C*V12* OPEN(ISCRU,FILE='ISCRU',FORM='UNFORMATTED',STATUS='OLD') C***** GISS VERSION FOLLOWS OPEN(ISCRU, FORM='UNFORMATTED',STATUS='OLD') C OPEN(ISCRU,FORM='UNFORMATTED',STATUS='OLD',SHARED,READONLY) ELSE WRITE(6,450) ISCRU 450 FORMAT('0 ENERGY-INDEPENDENT MATRICES WILL BE SAVED ', 1 'TEMPORARILY ON UNIT',I3) C*V12* OPEN(ISCRU,FILE='ISCRU',FORM='UNFORMATTED',STATUS='UNKNOWN') C***** GISS VERSION FOLLOWS OPEN(ISCRU, FORM='UNFORMATTED',STATUS='UNKNOWN') ENDIF ENDIF C WRITE(6,470) URED 470 FORMAT('0 REDUCED MASS FOR COLLISION =',F14.9,' A.M.U.') IF(JTOTL.LT.0) JTOTL=0 IF(JTOTU.LT.JTOTL) JTOTU=999 WRITE(6,480) JTOTL,JTOTU,JSTEP 480 FORMAT('0 CONTROL DATA FOR TOTAL ANGULAR MOMENTUM IS'/ 1 7X,'JTOT FROM',I4,' TO',I6,' IN STEPS OF',I4) IF(JTOTU.GE.999) WRITE(6,490) NCAC,DTOL,OTOL 490 FORMAT('0 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.999.AND.NNRGPG.GT.1) WRITE(6,491) NNRGPG 491 FORMAT('0 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('0 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('0 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('0***** 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('0 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('0 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('0 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('0 NUMDER=.TRUE. POTENTIAL DERIVATIVE WILL BE COMPUTED', & ' NUMERICALLY FROM POTENTIAL.') WRITE(6,650) PRINT,ISIGPR,ITHROW 650 FORMAT('0 PRINT LEVEL (PRNTLV) =',I3,' OTHER PRINT CONTROLS', 1 ' ISIGPR =',I2,' ITHROW =',I2) WRITE(6,660) 660 FORMAT('0',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 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 *** 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('0****** 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 N.B 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) CALL CHKSTR(NUSED) IC1=IXNEXT WRITE(6,660) C EFIRST=ENERGY(1)*CINT MXP=0 CALL GCLOCK(TITIME) TTIME=TITIME-TFIRST WRITE(6,700) TTIME,NUSED 700 FORMAT('0 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('0 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 IF(M.LE.MXRTM) GO TO 735 WRITE(6,732) MXRTM 732 FORMAT(/' *** ERROR. EXCEEDED LIMIT ON RTURNM. MXRTM =',I5) STOP 735 PHI=PHILW+PHIST*DBLE(M-1) IF(MSET.GT.0.AND.(M.LT.MSET.OR.M.GT.MHI)) GO TO 980 IF(PRINT.LT.4) GOTO 760 IF(ITHROW.NE.0) WRITE(6,710) TIT IF(ITHROW.EQ.0) WRITE(6,740) TIT 740 FORMAT('0',120A1) WRITE(6,750) JTOT,M 750 FORMAT('0 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,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 INITIALISE RTURN FOR IRMSET > 0 OPTION C IF(JTOT.EQ.JTOTL) RTURNM(M)=RMIN IK=1 RTURN=RTURNM(M) 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 769 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 N.B. 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.LE.0) GOTO 770 IS10=IXNEXT IS11=IS10+NSQ IXNEXT=IS11+NSQ 770 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,THETA,PHI, 2 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 769 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)) 771,774,773 771 IF(IECONV(J).LT.-2*MXP) GOTO 775 WRITE(6,772) JTOT,J 772 FORMAT('0 * * * WARNING. JTOT =',2I5,'-TH ENERGY PREVIOUSLY ', 1 'FAILED TO CONVERGE.') LCALC=.TRUE. GOTO 775 773 IF(JTOTU.LT.999) GOTO 774 IF(IECONV(J).LT.NCAC*MXP) GOTO 774 GOTO 775 774 LCALC=.TRUE. 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 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 778 IF(PRINT.LT.4) GOTO 790 IF(ITHROW.NE.0) WRITE(6,710) TIT2 IF(ITHROW.EQ.0) WRITE(6,740) TIT2 WRITE(6,780) JTOT,M,J,ETOT 780 FORMAT('0 JTOT =',I5,' SYMMETRY BLOCK =',I4,' ENERGY(', 1 I3,') =',F18.9,' (1/CM)') C C FOR SURFACE SCATTERING AT SUBSEQUENT ENERGY, C GET CORRESPONDING THETA FOR PRINTING C 790 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('0 NOTE: K VECTORS PARALLEL TO SURFACE WERE CALCULATED ', 1 'FOR ENERGY(1)'/' SUBSEQUENT ENERGY(',I3,') =',F10.4, 2 ' CORRESPONDS TO THETA =',F10.4,' DEGREES') ENDIF C C TEMPORARY STORAGE FOR HEADER, FINDRX ... IT1=IXNEXT IT2=IT1+MXLAM IXNEXT=IT2+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) IF(ICODE.NE.1 .OR. IRMSET.LE.0) GOTO 810 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),MXLAM,NPOTL,XEPS,ITYPE,PRINT) IF(RVFAC.EQ.0.D0) GOTO 810 RMID=RVFAC*RTURN IF(PRINT.GE.3.AND.RSTOP.GT.RMAX) WRITE(6,799) RSTOP,RMAX 799 FORMAT(' DRIVER(11/01/89) RMID IGNORES RSTOP.GT.RMAX',2F8.2) IF(PRINT.GE.3) WRITE(6,800) RMID,RVFAC 800 FORMAT('0 RMID =',F7.2,' OBTAINED FROM RVFAC =',F6.3) C C RESET IXNEXT TO DELETE TEMPORARY STORAGE 810 IXNEXT=IT1 C C AND SOLVE COUPLED EQUATIONS. C PROPAGATORS ARE CALLED FROM SUBROUTINE STORAG 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) GOTO 910 IF(PRINT.GE.2) WRITE(6,900) JTOT,M,J,ETOT,TTIME 900 FORMAT('0 ****** NO OPEN CHANNELS FOR JTOT =',I5, 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 910 CONTINUE C CALL OUTPUT(JTOT,X(IS9),X(ISJ),X(IS8),X(IS7),X(IS0),X(IS1), 1 X(IS2),CONV,NOPEN,M,MXPAR,WGHT,J,RM,PRINT,TTIME, 2 ENERGY,X(IOUT),X(IXJLEV),ISST,IECONV,MINJT,MAXJT, 3 NLEV,NQN,OTOL,DTOL,KSAVE,ISIGU,IPARTU,ISAVEU,ISIGPR) 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('0 FINISHED JTOT =',I5,' M =',I4,' ENERGY(',I3, 1 ') =',F18.9,10X,'STEP TIME =',F8.2,' SECS') C 960 ICODE=2 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('0'/'0 CALCULATION TERMINATED BY CONVERGENCE OF TOTAL ', 1 'CROSS SECTIONS.'/'0 DIAGONAL AND OFF-DIAGONAL TOLERANCES WERE', 2 2F9.5,' NCAC =',I3) GOTO 1000 ENDIF C IF(PRINT.GE.2 .AND. PRINT.LT.5) WRITE(6,970) 970 FORMAT('0') C SAVE RTURN FOR USE IN SUBSEQUENT JTOT FOR SAME M RTURNM(M)=RTURN 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, 1 ISIGU,JSTEP) 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(1X/' 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 N.B. MX MAY HAVE BEEN REDUCED, SO USE MXSAVE FOR ALLOCATED STORAGE 1020 CALL CHKSTR(NUSED) WRITE(6,1030) IPROGM,PDATE,TOTIME,NUSED,MXSAVE 1030 FORMAT('0'/'0 ',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,'|'/' |',15X,'THIS RUN USED',F9.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 CHCK6I(N,JL,A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION JL(4,N),A(1) DATA EPS/7.D-6/ WRITE(6,600) 600 FORMAT('0 CHCK6I. INPUT FUNCTIONS WILL BE CHECKED FOR ', & 'ORTHOGONALITY.') NERR=0 DO 1000 I1=2,N DO 1000 I2=1,I1-1 C SEE IF SAME J-VALUE IF (JL(1,I2).NE.JL(1,I1)) GO TO 1000 C CHECK THAT NK AGREE NK1=2*JL(1,I1)+1 NK2=2*JL(1,I2)+1 3000 IF (NK1.EQ.NK2) GO TO 1001 WRITE(6,699) I1,I2,NK1,NK2 699 FORMAT('0 ***** CHCK6I ERROR. FOR LEVELS',2I4,', NK NOT EQUAL.', & 2I5) NERR=NERR+1 GO TO 1000 1001 SUM=0.D0 DO 1100 II=1,NK1 1100 SUM=SUM+A(JL(4,I1)+II)*A(JL(4,I2)+II) IF (ABS(SUM).LE.EPS) GO TO 1000 WRITE(6,698) I1,I2,SUM 698 FORMAT('0 ***** CHCK6I ERROR. LEVEL',2I4,' ARE NOT ORTHOGONAL.', & ' OVERLAP =',D12.4) NERR=NERR+1 1000 CONTINUE IF (NERR.LE.0) RETURN WRITE(6,697) NERR 697 FORMAT('0 *****'/' ***** CHCK6I. NUMBER OF ERRORS =',I4/ 1 ' ***** EXECUTION TERMINATING UNLESS CHCK6I MODIFIED'/ 2 ' *****') STOP END SUBROUTINE COLIM(A,NLA,NUA,TOL,N) C COMPUTES LIMITS FOR BAND OF SIGNIFICANT ELEMENTS IN COLUMNS OF A IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(1),NLA(1),NUA(1) C FIND THE TOLERANCE LIMITS FOR THE TOPS(BEGINNINGS) OF THE C COLUMNS OF A NP1 = N + 1 NM1 = N - 1 LIMLO = 1 C THIS LOOP IS OVER THE COLUMNS OF A DO 40 K=1,N LIMHI = LIMLO + NM1 C THIS LOOP STARTS AT THE TOP OF THE K-TH COLUMN DO 10 J=LIMLO,LIMHI IF(ABS(A(J)).LE.TOL) GO TO 10 NLA(K) = J-LIMLO+1 GO TO 20 10 CONTINUE C THIS IS REACHED ONLY IF ALL ELEMENTS IN THE K-TH COLUMN ARE TINY NLA(K) = N NUA(K) = 1 GO TO 40 20 CONTINUE C FIND LIMITS FOR BOTTOMS OF COLUMNS C THIS LOOP STARTS AT THE BOTTOM END OF THE K-TH COLUMN DO 30 J=1,N IF(ABS(A(LIMHI)).LE.TOL) GO TO 30 NUA(K) = NP1 - J GO TO 40 30 LIMHI = LIMHI - 1 40 LIMLO = LIMLO + N RETURN END SUBROUTINE CONVRG(J,SR,SI,SROLD,SIOLD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE C C SUBROUTINE TO ASSIST IN THE ESTIMATION OF CONVERGENCE ERRORS C DIMENSION SR(1),SI(1),SROLD(1),SIOLD(1) COMMON/DRIVE/STEST,STEPS,STABIL,XCONV,RMIN,RMAX,XEPS, 1 DRNOW,DRMAX,RMID,TOLHI,RTURN,VTOL,ESHIFT,ERED,RMLMDA, 2 NOPEN,JKEEP,ISCRU,MAXSTP C C CHARACTER*6 CNAMES C DIMENSION CNAMES(3),LOCN(3),INDX(3) C DATA CNAMES/'ICONVU','IVAR','DR'/ C DATA INDX/3*0/ NAMELIST/CONV/ICONVU,IVAR,DR C NOPSQ=NOPEN*NOPEN IF(J.GT.1) GOTO 200 ISCRU=0 IVAR=0 DR=0.1D0 C------------------------------------------------------------------ C LOCN(1)=LOC(ICONVU) C LOCN(2)=LOC(IVAR) C LOCN(3)=LOC(DR) C CALL NAMLIS('&CONV ',CNAMES,LOCN,INDX,3,IEOF) C------------------------------------------------------------------ READ(5,CONV) IF(ICONVU.LT.0) GOTO 200 NSQOLD=NOPSQ DO 100 I=1,NOPSQ SROLD(I)=SR(I) 100 SIOLD(I)=SI(I) IF(ICONVU.EQ.0) GOTO 300 REWIND ICONVU WRITE(ICONVU) NOPSQ WRITE(ICONVU) (SROLD(I),I=1,NOPSQ) WRITE(ICONVU) (SIOLD(I),I=1,NOPSQ) GOTO 300 200 ICONVU=IABS(ICONVU) IF(ICONVU.EQ.0) GOTO 300 REWIND(ICONVU) READ(ICONVU) NSQOLD READ(ICONVU) (SROLD(I),I=1,NSQOLD) READ(ICONVU) (SIOLD(I),I=1,NSQOLD) WRITE(6,601) ICONVU 601 FORMAT('0 CONVERGENCE TESTS: REFERENCE S-MATRIX READ IN FROM ', 1 'CHANNEL',I3) 300 IF(NOPSQ.NE.NSQOLD) GOTO 600 ERRSM=0.D0 ERRTP=0.D0 DO 400 I=1,NOPSQ DIF = (SR(I)-SROLD(I))**2 + (SI(I)-SIOLD(I))**2 ERRSM=ERRSM+DIF DIF = (SR(I)**2 - SROLD(I)**2 + SI(I)**2 - SIOLD(I)**2)**2 ERRTP=ERRTP+DIF 400 CONTINUE C ERRSM=SQRT(ERRSM/DBLE(NOPSQ)) ERRTP=SQRT(ERRTP/DBLE(NOPSQ)) XSM=LOG10(MAX(ERRSM,1.D-30)) XTP=LOG10(MAX(ERRTP,1.D-30)) WRITE(6,602) RMIN,RMAX,RMID,STEPS,DRNOW,ERRSM,XSM,ERRTP,XTP 602 FORMAT('0 FOR RMIN =',F7.2,' RMAX =',F7.2,' RMID =',F7.2, 1 ' STEPS =',F8.1,' DR =',F7.4/ 2 ' RMS CHANGE IN S-MATRIX ELEMENTS IS ',7X, 3 E12.5,5X,'LOG IS',F8.3/ 4 ' RMS CHANGE IN TRANSITION PROBABILITIES IS ', 5 E12.5,5X,'LOG IS',F8.3) C IF(IVAR.EQ.0) DRNOW=DRNOW+DRNOW IF(IVAR.EQ.0) STEPS=STEPS/2.D0 IF(IVAR.EQ.1) RMIN=RMIN+DR IF(IVAR.EQ.2) RMID=RMID-DR IF(IVAR.EQ.3) RMAX=RMAX-DR RETURN C 600 WRITE(6,605) NOPSQ,NSQOLD 605 FORMAT('0*** ERROR IN CONVRG - NUMBER OF OPEN CHANNELS HAS ', 1 ' CHANGED'/5X,'NOPSQ =',I5,6X,'NSQOLD =',I5) RETURN END SUBROUTINE CONVRT(PTPOT,NP,V,NLEG,ABSC,WT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C CONVRT TAKES A SET OF FUNCTION VALUES EVALUATED AT GAUSS-LEGENDRE C QUADRATURE POINTS, AND CONVERTS THEM INTO A LEGENDRE SERIES, C USING THE RECURSION RELATIONSHIP FOR LEGENDRE POLYNOMIALS. C C PTPOT IS THE INPUT ARRAY, EVALUATED AT NP POINTS C V IS THE OUTPUT ARRAY: V(N) CONTAINS THE COEFFICIENTS OF C THE (N-1)TH ORDER LEGENDRE POLYNOMIAL C C IF ABSC(1) NE 0.0 ON INPUT, THE SUBROUTINE ASSUMES THAT THE C ARRAYS ABSC AND WT ALREADY CONTAIN THE APPROPRIATE C GAUSS-LEGENDRE ABSCISSAE AND WEIGHTS; OTHERWISE, GAUSSP C (OR THE NAG ROUTINE D01BBF) IS CALLED TO OBTAIN THEM. C C EXTERNAL D01BAZ DIMENSION PTPOT(NP),V(NLEG),ABSC(NP),WT(NP) IF(NLEG.LE.NP) GOTO 5 WRITE(6,601)NLEG,NP 601 FORMAT(' *** ERROR IN CONVRT, NLEG =',I4,' > NP =',I4) STOP 5 IF(ABSC(1).EQ.0.0D0) THEN CALL GAUSSP(-1.D0,1.D0,NP,ABSC,WT) C IFAIL=0 C CALL D01BBF(D01BAZ,-1.0D0,1.0D0,0,NP,WT,ABSC,IFAIL) ENDIF DO 10 N=1,NLEG 10 V(N)=0.0D0 C DO 100 M=1,NP PTWT=PTPOT(M)*WT(M) X=ABSC(M) P0=1.0D0 P1=X DO 100 K=1,NLEG GOTO(30,40),K TEMP=(DBLE(2*K-3)*X*P1 - DBLE(K-2)*P0) / DBLE(K-1) P0=P1 P1=TEMP V(K) = V(K) + (DBLE(K)-0.5D0) * P1 * PTWT GOTO 100 30 V(1) = V(1) + 0.5D0 * PTWT GOTO 100 40 V(2) = V(2) + 1.5D0 * PTWT * X 100 CONTINUE RETURN END SUBROUTINE DASIZE(ILSU,MXREC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE MXUSED,IX PARAMETER (NREC=20000) DIMENSION IX(6,NREC) DIMENSION IR2(2),IS2(2) EQUIVALENCE (R,IR1,IR2(1)),(S,IS1,IS2(1)) COMMON /ASSVAR/IDA C C DYNAMIC STORAGE COMMON BLOCK ... C NEEDED FOR NIPR; PREVIOUSLY PASSED IN COMMON /INTPAC/ COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1) C DATA MAX/NREC/ C MXREC=MAX ILSU=999 WRITE(6,601) MAX 601 FORMAT( ' *** *** NUMBER OF SIMULATED RECORDS =',I7) RETURN C ENTRY DAOPEN MXUSED=0 WRITE(6,600) 600 FORMAT(/' *** *** IN-CORE DA SIMULATION ROUTINE HAS CONTROL.', 1 /' *** *** DA FILE WILL NOT BE USED.') C IF(NIPR.EQ.1 .OR. NIPR.EQ.2) GOTO 1000 WRITE(6,602) NIPR 602 FORMAT(' *** ERROR IN DASIZE/DAOPEN: NIPR =',I3,' INVALID') STOP 1000 RETURN C ENTRY DARD1(I1,I2,I3,I4,I5,I6) I1=IX(1,IDA) I2=IX(2,IDA) I3=IX(3,IDA) I4=IX(4,IDA) I5=IX(5,IDA) I6=IX(6,IDA) RETURN C ENTRY DAWR1(I1,I2,I3,I4,I5,I6) MXUSED=MAX0(MXUSED,IDA) IX(1,IDA)=I1 IX(2,IDA)=I2 IX(3,IDA)=I3 IX(4,IDA)=I4 IX(5,IDA)=I5 IX(6,IDA)=I6 RETURN C ENTRY DARD2(I1,I2,X1,X2) I1=IX(1,IDA) I2=IX(2,IDA) IF(NIPR.EQ.1) THEN IR1=IX(3,IDA) IS1=IX(4,IDA) ELSE IR2(1)=IX(3,IDA) IR2(2)=IX(4,IDA) IS2(1)=IX(5,IDA) IS2(2)=IX(6,IDA) ENDIF X1=R X2=S RETURN C ENTRY DAWR2(I1,I2,X1,X2) MXUSED=MAX0(MXUSED,IDA) IX(1,IDA)=I1 IX(2,IDA)=I2 R=X1 S=X2 IF(NIPR.EQ.1) THEN IX(3,IDA)=IR1 IX(4,IDA)=IS1 ELSE IX(3,IDA)=IR2(1) IX(4,IDA)=IR2(2) IX(5,IDA)=IS2(1) IX(6,IDA)=IS2(2) ENDIF RETURN C ENTRY DACLOS WRITE(6,610) MXUSED,MAX 610 FORMAT(/' *** IN-CORE DA SIMULATOR USED',I10,' OF THE',I10, 1 ' ALLOCATED RECORDS') RETURN END SUBROUTINE DASCAT(N, NSQ, MXLAM, NPOTL, 1 SR, SI, U, VL, IV, EINT, CENT, WVEC, L, NB, 2 P, Y1, Y2, Y3, Y4, 3 ICODE, IPRINT, IC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C *** --------------------------------------------------------------- C *** ROUTINE TO PERFORM A SCATTERING CALCULATION USING DAPROP. C *** ON EXIT SR AND SI CONTAIN THE S MATRIX. C *** SR IS USED INTERNALLY TO HOLD THE LOG DERIVATIVE MATRIX C *** IN ORDER TO ECONOMISE ON WORKSPACE. C *** --------------------------------------------------------------- C *** ICODE.EQ.2 FOR SUBSEQUENT ENERGIES. C *** C DIMENSION STATEMENTS FOR ARGUMENT LIST 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) 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 C DYNAMIC STORAGE COMMON BLOCK ... COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1) C C THE FOLLOWING VARIABLES FROM COMMON/DRIVE/ ARE USED WITH THIS C PROPAGATOR: STEPS,RMIN,RMAX,ERED,RMLMDA,NOPEN,ISCRU C LOGICAL IREAD,IWRITE 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 --------------------------------------------------------------- C C CALCULATE WAVEVECTORS AND STEP SIZE C 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) NSTEPS=WMAX*STEPS*(RMAX-RMIN)/PI RBEGIN=RMIN REND=RMAX IF (IWRITE) WRITE (ISCRU) RBEGIN,REND,NSTEPS GO TO 60 40 READ (ISCRU) RBEGIN,REND,NSTEPS 60 CONTINUE ISTART=0 C C PROPAGATE LOG DERIVATIVE MATRIX THROUGH THE SCATTERING REGION C --------------------------------------------------------------- IF(N.EQ.1) GOTO 90 CALL DAPROP(U, SR, N, & RBEGIN, REND, NSTEPS, IREAD, IWRITE, ISCRU, & Y1, Y2, Y3, Y4, & P, VL, IV, ERED, EINT, CENT, RMLMDA, & MXLAM, NPOTL, ISTART, NODES) C --------------------------------------------------------------- IF (IPRINT.GE.3) WRITE (6,1000) RBEGIN,REND,NSTEPS 1000 FORMAT('0 DAPROP. LOG DERIVATIVE MATRIX INTEGRATED FROM ', & F12.4,' TO ',F12.4,' IN ',I6,' STEPS.') C C SORT CHANNELS BY ASYMPTOTIC ENERGY C 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 GOTO 100 C C SPECIAL CASE FOR EFFICIENT SINGLE CHANNEL CALCULATIONS C 1/21/93 CHANGES TO DYNAMIC STORAGE: N.B IT5 FROM STORAG IS C PASSED AS ARGUMENT IC TO FOLLOW CODING IN EARLIER VERSIONS. C 90 NPT=NSTEPS+1 ISVMEM=IXNEXT IT1=IC IT2=IT1+NPT IT3=IT2+NPT IT4=IT3+NPT IT5=IT4+NPT IC1=IT5+NPT ITP=IT3 IC2=ITP+NPT*MXLAM IXNEXT=MAX0(IC1,IC2) NUSED=0 CALL CHKSTR(NUSED) CALL ODPROP(SR, X(IT1), X(IT2), X(IT3), X(IT4), X(IT5), & RBEGIN, REND, NPT, IREAD, IWRITE, ISCRU, & X(ITP), VL, IV, ERED, EINT, CENT, RMLMDA, & MXLAM, NPOTL, ISTART, NODES) C RESTORE STORAGE POINTER TO RECOVER TEMPORARY STORAGE. IXNEXT=ISVMEM C --------------------------------------------------------------- IF (IPRINT.GE.3) WRITE (6,1010) RBEGIN,REND,NSTEPS 1010 FORMAT('0 ODPROP. LOG DERIVATIVE INTEGRATED FROM ', & F12.4,' TO ',F12.4,' IN ',I6,' STEPS.') 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 DELRD(DR,CDIAG,COFF,TOL,DRMAX,E1,EN,RNOW,RMAX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C------------------------------------------------------------------- C THIS ROUTINE IS MODIFIED VERSION OF ROY GORDON'S QCPE PROGRAM. C THIS VERSION FOR SCATTERING CALCULATION C------------------------------------------------------------------- C ADJUST THE STEP SIZE TO TRY TO KEEP MAX(CDIAG,COFF) = TOL RTOL = 0.80D0*TOL C------------------------------------------------------------------- C FIND CORRECTION FACTOR FROM DIAGONAL PERTURBATIONS C------------------------------------------------------------------- IF (CDIAG .NE. 0.D0) GO TO 100 C------------------------------------------------------------------- C CASE IN WHICH DIAGONAL PERTURBATIONS VANISH C------------------------------------------------------------------- DFACT = 2.5D0 GO TO 110 C------------------------------------------------------------------- C DIAGONAL PERTURBATIONS VARY ROUGHLY AS THE FIFTH POWER OF STEP SIZE C------------------------------------------------------------------- 100 DFACT = (RTOL/CDIAG)**0.333D0 C------------------------------------------------------------------- C FIND CORRECTION FACTOR FROM OFF-DIAGONAL PERTURBATIONS C------------------------------------------------------------------- 110 IF (COFF .NE. 0.D0) GO TO 120 C------------------------------------------------------------------- C CASE IN WHICH OFF-DIAGONAL PERTURBATIONS VANISH C------------------------------------------------------------------- OFACT = 2.5D0 GO TO 130 C------------------------------------------------------------------- C OFF-DIAGONAL PERTURBATIONS VARY ROUGHLY AS CUBE OF STEP SIZE C------------------------------------------------------------------- 120 OFACT = (RTOL/COFF)**0.333D0 C------------------------------------------------------------------- C FIND MINIMUM FACTOR C------------------------------------------------------------------- 130 FACTOR = MIN(DFACT,OFACT) IF (EN .GT. 0.D0) GO TO 150 IF (E1 .GT. 0.D0) GO TO 140 C------------------------------------------------------------------- C THIS IS REACHED ONLY WHEN ALL CHANNELS ARE IN THEIR CLASSICALLY C FORBIDDEN REGIONS. THEN ACCURACY IS QUITE SENSITIVE TO CHANGES C IN STEP SIZE. C HENCE IN THIS REGION MAKE ONLY CAUTIOUS CHANGES IN STEP SIZE C------------------------------------------------------------------- IF (FACTOR .GT. 1.15D0) FACTOR = 1.15D0 GO TO 170 C------------------------------------------------------------------- C THIS IS REACHED WHEN SOME CHANNELS ARE CLASSICAL AND OTHERS NOT C------------------------------------------------------------------- 140 IF (FACTOR .GT. 1.20D0) FACTOR = 1.20D0 GO TO 170 C------------------------------------------------------------------- C THIS IS REACHED WHEN ALL CHANNELS ARE CLASSICAL. C THEN THE STEP SIZE IS OFTEN INCREASING RAPIDLY, AND ALSO THE C ACCURACY VARIES MORE SLOWLY WITH STEP SIZE. C THUS WE MAKE BOLDER INCREASES IN THE STEP SIZE, TO KEEP THE C CORRECTIONS OF THE SAME ORDER OF MAGNITUDE AS BEFORE C TEST TO SEE HOW FAR WE HAVE INTEGRATED C------------------------------------------------------------------- 150 IF (RNOW .GT. (0.10D0*RMAX)) GO TO 160 IF (FACTOR .GT. 1.6D0) FACTOR = 1.6D0 GO TO 170 C------------------------------------------------------------------- C CHOOSE FACTOR IN FAR AYSMPTOTIC REGION C------------------------------------------------------------------- 160 IF (FACTOR .GT. 2.5D0) FACTOR = 2.5D0 C------------------------------------------------------------------- C SET NEW STEP SIZE C------------------------------------------------------------------- 170 DR = DR*FACTOR C------------------------------------------------------------------- C CHECK AGAINST DRMAX AND EXCESSIVE GROWTH OF CLOSED CHANNELS C------------------------------------------------------------------- IF (EN .GE. 0.D0) GO TO 175 DREXP = 4.D0/SQRT(-EN) IF (DR .GT. DREXP) DR = DREXP 175 IF (DR .GT. DRMAX) DR = DRMAX IF (DR .LT. 1.0D-06*DRMAX) GO TO 180 RETURN 180 WRITE (6,1000) DR,CDIAG,COFF,TOL,DRMAX,E1,EN,RNOW,RMAX STOP C------------------------------------------------------------------- C FORMATS C------------------------------------------------------------------- 1000 FORMAT('0 * * * ERROR IN DELRD. STEP SIZE =',E20.6, 1 ' IS TOO SMALL'/'0',24X,'PARAMETERS PASSED ARE', 2 ' CDIAG, COFF, TOL, DRMAX, E1, EN, RNOW, RMAX'/ 4 25X,9(E10.3,1X)) C----------------***END-DELRD***------------------------------------ END SUBROUTINE DERMAT(IDER,W,N,R,P,VL,IV,CENT,RMLMDA, 1 MXLAM,NPOTL,NUMDER) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE RSAVE LOGICAL NUMDER C C EVALUATES THE IDER'TH DERIVATIVE OF THE POTENTIAL MATRIX AT RADIUS C W = VCOUPL + VCENT C ORDER OF THE REAL SYMMETRIC MATRIX W IS N C THE FULL MATRIX IS COMPUTED C VL IS THE PREVIOUSLY COMPUTED MATRIX OF THE COUPLING POTENTIAL C IV IS AN INDEX ARRAY MAPPING P ONTO VL, SUCH THAT VL(I) IS C A COEFFICIENT TO MULTIPLY P(IV(I)) C CENT(I) IS L*(L+1) FOR THE I-TH CHANNEL C RMLMDA IS THE SQUARE OF THE RATIO OF RM TO THE DEBROGLIE C WAVELENGTH AT RELATIVE ENERGY EPSILON C RMLMDA = 2.*URED*RM**2*EPSIL/HBAR**2 C RMLMDA MULTIPLIES THE POTENTIAL IN UNITS OF EPSIL C DIMENSION W(N,N),VL(1),IV(1),CENT(N),P(MXLAM) DATA DEL/1.D-3/, RSAVE/-999.D0/ C IF(NUMDER) GOTO 5 C C COMPUTE THE RADIAL PARTS OF THE POTENTIAL ANALYTICALLY CALL POTENL(IDER,MXLAM,NPOTL,IDUM1,R,P,IDUM2) C IDUM1 AND IDUM2 ARE DUMMY ARGUMENTS HERE. GOTO 14 C C NUMERICAL DERIVATIVE OPTION. C NOTE THAT IF IDER = 2 THIS ASSUMES THAT C THE POTENTIAL ITSELF IS ALREADY IS ALREADY IN THE FIRST C MXLAM ELEMENTS OF P. THIS IS NOT TRUE IF DERMAT HAS BEEN C CALLED MORE RECENTLY THAN WAVMAT, SO THE IDER = 2 CALL C MUST PRECEDE THE IDER = 1 CALL. C C FIRST SEE WHETHER DERMAT HAS BEEN CALLED BEFORE FOR THIS C VALUE OF R, AND IF SO SKIP POTENTIAL EVALUATIONS C 5 IF(R.EQ.RSAVE) GOTO 8 RSAVE=R RR=R-DEL CALL POTENL(0,MXLAM,NPOTL,IDUM1,RR,P(MXLAM+1),IDUM2) RR=R+DEL CALL POTENL(0,MXLAM,NPOTL,IDUM1,RR,P(2*MXLAM+1),IDUM2) C 8 DO 10 I=1,MXLAM P1=P(MXLAM+I) P2=P(2*MXLAM+I) IF(IDER.EQ.1) P(I) = (P2-P1)/(2.D0*DEL) 10 IF(IDER.EQ.2) P(I) = (P2+P1-2.D0*P(I)/RMLMDA)/(DEL*DEL) C 14 DO 15 I=1,MXLAM 15 P(I)=RMLMDA*P(I) C CALL WAVVEC(VL,P,IV,W,N,NPOTL) C C NOW COMPUTE THE DIAGONAL CONTRIBUTIONS W(I,I). C IF(IDER.EQ.1) RSQ=-2.D0/R**3 IF(IDER.EQ.2) RSQ= 6.D0/R**4 DO 20 I=1,N W(I,I) = W(I,I) + RSQ*CENT(I) 20 CONTINUE RETURN END SUBROUTINE DVFREE(UJ,UJP,UN,UNP,WRONS,L,N,WV,R,NB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DOUBLE PRECISION ASYMPTOTIC FUNCTIONS FOR MATCHING TO S-MATRIX. DIMENSION UJ(N),UJP(N),UN(N),UNP(N),WRONS(N),WV(N) DIMENSION L(N),NB(N) DO 3000 I=1,N NX=NB(I) DW=WV(NX) DARG=DW*R CALL RBES(L(NX),DARG,UJ(NX),UJP(NX),UN(NX),UNP(NX)) UJP(NX)=UJP(NX)*DW UNP(NX)=UNP(NX)*DW 3000 WRONS(NX)=(UJ(NX)*UNP(NX)-UJP(NX)*UN(NX))/SQRT(DW) RETURN END SUBROUTINE DVSCAT(N,NSQ,MXLAM,NPOTL, 1 SR,SI,A,VL,IV,EINT,CENT,WV,L,NB, 2 P,Y,YP,F,XM,YM,DIAG,ESHIFT,ICODE,PRINT) C C DEVOGELAERE INTEGRATION (DOUBLE PRECISION) C INCLUDING START ROUTINE, SUPPRESSION OF CLOSED-CHANNEL GROWTH C AND S-MATRIX DETERMINATION IN ASYMPTOTIC LIMIT. C FOLLOWS OUTLINE OF PAUL MCGUIRE'S PROGRAM. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL IREAD,IWRITE,END INTEGER L(N),NB(N),IV(1) INTEGER PRINT DIMENSION Y(NSQ,4),YP(NSQ,2),F(NSQ,4),A(NSQ),XM(NSQ),YM(NSQ), 1 DIAG(N),P(MXLAM),SR(NSQ),SI(NSQ),VL(2),EINT(N),CENT(N),WV(N) DIMENSION R(4) C C INDICES ON Y, YP, F ARE (ITH SOLN. COMP, NTH SOLN, KTH R-VALUE) C C COMMON FROM DRIVER COMMON/DRIVE/STEST,STEPS,STAB,CONV,RMIN,RMAX,DUMMY(8), 1 ERED,RMLMDA,NOPEN,JKEEP,ISCRU,MAXSTP C C STEPS IS NO. OF STEPS PER (SHORTEST) WAVELENGTH. C STAB IS NUMBER OF STEPS TAKEN BEFORE STABILIZATION. C C MAX. NUMBER OF TRIALS TO CONVERGE S-MATRIX IN ASYMPTOTIC REGION. DATA MXSTRY/20/ C IREAD = ICODE.EQ.2 .AND. ISCRU.GT.0 IWRITE = ICODE.EQ.1 .AND. ISCRU.GT.0 IF(IREAD .AND. PRINT.GE.5) WRITE(6,668) 668 FORMAT('0 DEVOGELAERE PROPAGATION WILL USE STORED INITIAL R,', 1 ' STEP SIZE AND POTENTIAL MATRICES') C C ZERO STORAGE . . . C NP1=N+1 DO 800 IJ=1,NSQ SR(IJ)=0.D0 800 SI(IJ)=0.D0 DO 900 I=1,2 DO 900 IJ=1,NSQ 900 YP(IJ,I)=0.D0 DO 1000 I=1,4 DO 1000 IJ=1,NSQ Y(IJ,I)=0.D0 1000 F(IJ,I)=0.D0 C NSTRY=0 RMSAVE=RMAX C C ********** START INTEGRATION ********** CALL DVSTRT(RMIN,STEPS,ERED,RMLMDA,N,MXLAM,NPOTL,NOPEN,PRINT, & A,DIAG,P,VL,IV,EINT,CENT,NB,WV,Y(1,2),YP(1,1),HH, & ISCRU,IREAD,IWRITE) IF (NOPEN.LE.0) GOTO 9000 NSTAB=STAB NSTAB=MAX0(NSTAB,1) H2=HH/2.D0 R(2)=RMIN NSTEP=1 C GET F(,,2) FROM Y(,,2) R4=R(2) IF(.NOT.IREAD) GOTO 1200 READ(ISCRU) A DO 1100 IJ=1,NSQ,NP1 1100 A(IJ)=A(IJ)-ESHIFT GOTO 1300 1200 CALL WAVMAT(A,N,R4,P,VL,IV,ERED,EINT,CENT,RMLMDA,DIAG, 1 MXLAM,NPOTL) IF(IWRITE) WRITE(ISCRU) A 1300 CALL DSYMM('L','L',N,N,1.D0,A,N,Y(1,2),N,0.D0,F(1,2),N) CALL DAXPY(NSQ,1.D0,F(1,2),1,Y(1,2),1) CALL DAXPY(NSQ,-H2,YP(1,1),1,Y(1,1),1) CALL DAXPY(NSQ,0.5D0*H2*H2,F(1,2),1,Y(1,1),1) C GET F(,,1) FROM THIS Y(,,1). NEEDS POTENTIAL AT R(1) R(1)=R(2)-H2 R4=R(1) IF(.NOT.IREAD) GOTO 1800 READ(ISCRU) A DO 1700 IJ=1,NSQ,NP1 1700 A(IJ)=A(IJ)-ESHIFT GOTO 1900 1800 CALL WAVMAT(A,N,R4,P,VL,IV,ERED,EINT,CENT,RMLMDA,DIAG, 1 MXLAM,NPOTL) IF(IWRITE) WRITE(ISCRU) A 1900 CALL DSYMM('L','L',N,N,1.D0,A,N,Y(1,1),N,0.D0,F(1,1),N) C C ********** MAIN BODY OF ITERATION ********** C PROPAGATE FROM (-1/2) AND (0) TO (1/2) AND (1). 2000 CONTINUE CALL DAXPY(NSQ,1.D0,Y(1,2),1,Y(1,3),1) CALL DAXPY(NSQ,H2,YP(1,1),1,Y(1,3),1) CALL DAXPY(NSQ,H2*H2*4.D0/6.D0,F(1,2),1,Y(1,3),1) CALL DAXPY(NSQ,-H2*H2/6.D0,F(1,1),1,Y(1,3),1) R(3)=R(2)+H2 R4=R(3) IF(.NOT.IREAD) GOTO 2200 READ(ISCRU) A DO 2100 IJ=1,NSQ,NP1 2100 A(IJ)=A(IJ)-ESHIFT GOTO 2300 2200 CALL WAVMAT(A,N,R4,P,VL,IV,ERED,EINT,CENT,RMLMDA,DIAG, 1 MXLAM,NPOTL) 2300 IF(IWRITE) WRITE(ISCRU) A CALL DSYMM('L','L',N,N,1.D0,A,N,Y(1,3),N,0.D0,F(1,3),N) CALL DAXPY(NSQ,1.D0,Y(1,2),1,Y(1,4),1) CALL DAXPY(NSQ,HH,YP(1,1),1,Y(1,4),1) CALL DAXPY(NSQ,HH*HH/6.D0,F(1,2),1,Y(1,4),1) CALL DAXPY(NSQ,HH*HH/3.D0,F(1,3),1,Y(1,4),1) R(4)=R(3)+H2 R4=R(4) IF(.NOT.IREAD) GOTO 2700 READ(ISCRU) A DO 2600 IJ=1,NSQ,NP1 2600 A(IJ)=A(IJ)-ESHIFT GOTO 2800 2700 CALL WAVMAT(A,N,R4,P,VL,IV,ERED,EINT,CENT,RMLMDA,DIAG, 1 MXLAM,NPOTL) IF(IWRITE) WRITE(ISCRU) A 2800 CALL DSYMM('L','L',N,N,1.D0,A,N,Y(1,4),N,0.D0,F(1,4),N) CALL DAXPY(NSQ,1.D0,YP(1,1),1,YP(1,2),1) CALL DAXPY(NSQ,HH/6.D0,F(1,2),1,YP(1,2),1) CALL DAXPY(NSQ,HH/6.D0,F(1,4),1,YP(1,2),1) CALL DAXPY(NSQ,HH*4.D0/6.D0,F(1,3),1,YP(1,2),1) NSTEP=NSTEP+1 C C ********** THIS ENDS DEVOGELAERE CYCLE ********** C NOPLOC=0 DO 2900 I=1,NSQ,NP1 2900 IF(A(I).LT.0.D0) NOPLOC=NOPLOC+1 C END=R4.GT.RMAX .AND. NOPLOC.GE.NOPEN IF(IREAD) READ(ISCRU) END IF(END) GOTO 3000 C C ********** STABILIZATION EVERY NSTAB STEPS ********** 4000 IF(NSTEP-NSTAB*(NSTEP/NSTAB).NE.0) GOTO 5000 IF (PRINT.GT.12) WRITE(6,673) R(4) 673 FORMAT(' STABILIZATION DONE AT R =',E12.4) C FIRST 2 COLS OF Y AND F AND ALSO A USED AS SCRATCH IN STABIL. CALL STABIL(N,NB,Y(1,4),YP(1,2),F(1,3),F(1,4), & A,Y(1,1),Y(1,2),F(1,1),F(1,2)) C C ********** RE-INITIALIZE FOR NEXT CYCLE OF INTEGRATION ********* 5000 R(1)=R(3) R(2)=R(4) IF(IWRITE) WRITE(ISCRU) IREAD CALL DCOPY(NSQ,YP(1,2),1,YP(1,1),1) CALL DCOPY(NSQ,Y(1,3),1,Y(1,1),1) CALL DCOPY(NSQ,Y(1,4),1,Y(1,2),1) CALL DCOPY(NSQ,F(1,3),1,F(1,1),1) CALL DCOPY(NSQ,F(1,4),1,F(1,2),1) DO 5200 IJ=1,NSQ YP(IJ,2)=0.D0 Y(IJ,3)=0.D0 5200 Y(IJ,4)=0.D0 GOTO 2000 C 3000 CONTINUE IF ((PRINT.GE.2.AND.NSTRY.LE.0) .OR. PRINT.GE.12) & WRITE(6,601) NSTEP,R(4) 601 FORMAT(' INTEGRATION REACHED ASYMPTOTIC LIMIT IN', & I5,' STEPS. R =',D12.4) C C ********** ASYMPTOTIC REGION - CALCULATE S-MATRIX ********** NOPSQ=NOPEN*NOPEN C USE FIRST 2 COLS OF Y AND F FOR REGULAR AND IRREGULAR BESSEL FNS. C AND DERIVATIVES - UJ, UJP, UN, AND UNP. C USE FIRST COL. OF A FOR WRONSKIAN/SQRT(WV). CALL DVFREE(Y(1,1),Y(1,2),F(1,1),F(1,2),A,L,NOPEN,WV,R(4),NB) C FORM TRANSPOSE OF X- AND Y- MATRICES DO 3200 J=1,NOPEN IJ=J DO 3100 I=1,NOPEN NX=NB(I) NY=NX+N*(J-1) XM(IJ)=(F(NX,2)*Y(NY,4)-F(NX,1)*YP(NY,2)) / A(NX) YM(IJ)=(Y(NX,2)*Y(NY,4)-Y(NX,1)*YP(NY,2)) / A(NX) 3100 IJ=IJ+NOPEN 3200 CONTINUE C GET K-MATRIX FROM SOLN TO LINEAR EQNS,REPLACES RHS DO 3300 I=1,NOPSQ 3300 A(I)=YM(I) CALL DGESV(NOPEN,NOPEN,XM,NOPEN,Y,A,NOPEN,IER) IF (IER.EQ.0) GOTO 3400 WRITE(6,688) 688 FORMAT('0 * * * WARNING. LOSS OF ACCURACY IN SOLVING FOR K.') CALL OUTERR(11) C C FORCE SYMMETRY ON K-MATRIX AND CALCULATE S MATRIX C 3400 CALL RSYM(NOPEN, A, STEST, PRINT) CALL KTOS(A,XM,YM,NOPEN) C C TEST FOR CONVERGENCE OF SR, SI C TEST=0.D0 DO 3500 I=1,NOPSQ TEST=MAX(TEST,ABS(SR(I)-XM(I)),ABS(SI(I)-YM(I))) SR(I)=XM(I) 3500 SI(I)=YM(I) C IF(IREAD) GOTO 9000 IF(TEST.GT.STEST) GOTO 3600 IF(PRINT.GE.2) WRITE(6,686) NSTRY,R(4),TEST 686 FORMAT(' S-MATRIX CONVERGED AFTER',I3,' TRIES IN ', & 'ASYMPTOTIC REGION. R =',D12.4,'. TEST =',E12.4) GOTO 9000 3600 IF(NSTRY.GT.0 .AND. PRINT.GE.2) WRITE(6,687) NSTRY,R(4),TEST 687 FORMAT(' S-MATRIX NOT CONVERGED AFTER',I3, & ' TRIES. R =',D12.4,'. LARGEST CHANGE =',D12.4) IF (NSTRY.LT.MXSTRY) GOTO 3700 C SET 'CONV' FLAG FOR OUTPUT ROUTINE. . . CONV=-1.D0 GOTO 9000 3700 NSTRY=NSTRY+1 RMAX=RMAX+STEPS*HH GOTO 5000 C C COMMON RETURN POINT C RESTORE RMAX 9000 RMAX=RMSAVE IF(IWRITE) WRITE(ISCRU) IWRITE RETURN END SUBROUTINE DVSTRT(R,STEPS,ERED,RMLMDA,N,MXLAM,NPOTL,NOPEN,PRINT, & W,DIAG,P,VL,IV,EINT,CENT,NB,WV,U,UP,HH,ISCRU,IREAD,IWRITE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C INTEGER PRINT LOGICAL SURF,IREAD,IWRITE DIMENSION W(N,N),DIAG(N),P(2),VL(2),IV(1),EINT(N),CENT(N),NB(N), & WV(N),U(N,N),UP(N,N) C C PROVIDE STARTING SOLUTION AND DERIVATIVE FOR DEVOGELAERE. C ALSO PICK STEP SIZE, HH. C STEPS IS NO. OF STEPS PER (SHORTEST) WAVELENGTH. C THIS IS SIMPLEST VERSION, SIMILAR TO THAT OF MCGUIRE. C RSAVE=R SURF=R.LT.0.D0 DRMIN=ABS(R/STEPS) C C * * * * * ORDER BASIS FUNCTIONS ON INCREASING INTERNAL ENERGY. DO 3000 I=1,N 3000 NB(I)=I IF (N.LE.1) GO TO 1000 NM1=N-1 DO 3100 I=1,NM1 IP1=I+1 DO 3100 J=IP1,N IF (EINT(NB(I)).LE.EINT(NB(J))) GO TO 3100 IT=NB(I) NB(I)=NB(J) NB(J)=IT 3100 CONTINUE C C * * * * * SEE THAT ALL CHANNELS (IN FREE BASIS) ARE CLOSED. IF(IREAD) GOTO 2000 1000 CALL WAVMAT(W,N,R,P,VL,IV,ERED,EINT,CENT,RMLMDA,DIAG, 1 MXLAM,NPOTL) IF(PRINT.GT.12) WRITE(6,698) R,(P(I),I=1,MXLAM) 698 FORMAT('0 POTENTIAL ARRAY AT R =',F10.4/(10(1X,D12.5))) DO 1100 I=1,N IF (PRINT.GT.12) WRITE(6,699) I,W(I,I) 699 FORMAT(' FOR CHANNEL',I4, ' V(RMIN) - E =',D13.4) IF (W(I,I).GT.0.D0) GO TO 1100 R=R-DRMIN IF(SURF .OR. R.GT.0.D0) GO TO 1000 WRITE(6,600) 600 FORMAT('0 * * * ERROR. RMIN LESS THAN ZERO IN DVSTRT.', 1 ' POTENTIAL MAY BE UNPHYSICAL') STOP 1100 CONTINUE C IF(R.NE.RSAVE) WRITE(6,602)RSAVE,R 602 FORMAT('0 * * * WARNING. DVSTRT HAS CHANGED RMIN FROM ',F6.2, & ' TO ',F6.2,' TO ENSURE THAT ALL CHANNELS ARE LOCALLY CLOSED') C C * * * * * INITIALIZE U, UP. 2000 DO 4000 I=1,N DO 4000 J=1,N U(I,J)=0.D0 UP(I,J)=0.D0 IF (I.EQ.NB(J)) UP(I,J)=1.D-8 4000 CONTINUE C * * * * * INITIALIZE NOPEN, WV. PICK STEP SIZE. NOPEN=0 BIG=0.D0 DO 5000 I=1,N DIF=ERED-EINT(I) IF (DIF.LE.0.D0) GO TO 5100 NOPEN=NOPEN+1 5100 WV(I)=SIGN(SQRT(ABS(DIF)),DIF) 5000 BIG=MAX(BIG,WV(I)) IF (NOPEN.LE.0) RETURN C CALCULATE STEP SIZE FROM LARGEST WVEC. HH=3.1416D0/(BIG*STEPS) IF(IWRITE) WRITE(ISCRU) R,HH IF(IREAD) READ(ISCRU) R,HH IF (PRINT.GE.2) WRITE(6,601) R,HH 601 FORMAT('0 INTEGRATION STARTED AT RMIN =',D12.4, & '. STEP SIZE =',D12.4) RETURN END SUBROUTINE EAVG(NT,T,NGP,E,NNRG,MXNRG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION E(1),T(1) DIMENSION A(20),W(20) C THIS ROUTINE SETS UP ENERGIES FOR NGP-POINT GAUSS-LAGUERRE INTEG. C AT SPECIFIED TEMPERATURES (DEG. KELVIN). DATA XK/.6950305D0/ DATA A/.585786437627D0, 3.414213562373D0, 2 0.415774556783D0, 2.294280360279D0, 6.289945082937D0, 3 0.322547689619D0, 1.745761101158D0, 4.536620296921D0, 4 9.395070912301D0, 0.263560319718D0, 1.413403059107D0, 5 3.596425771041D0, 7.085810005859D0, 12.640800844276D0, 6 6*0.D0/ DATA W/ 0.853553390593D0, 0.146446609407D0, 0.711093009929D0, 8 0.278517733569D0, 0.103892565016D-1, 0.603154104342D0, 9 0.357418692438D0, 0.388879085150D-1, 0.539294705561D-3, A 0.521755610583D0, 0.398666811083D0, 0.759424496817D-1, B 0.361175867992D-2, 0.233699723858D-4, 6*0.D0/ NGP=MAX0(2,MIN0(6,IABS(NGP))) IST=NGP*(NGP-1)/2-1 WRITE(6,600) NGP 600 FORMAT('0 ENERGY VALUES WILL BE GENERATED TO FACILITATE',I4, 1 '-POINT GAUSS-LAGUERRE INTEGRATION OVER BOLTZMANN DISTRIBUTION') NN=0 DO 1000 I=1,NT IF (NN+NGP.LE.MXNRG) GO TO 1010 WRITE(6,601) I,T(I) 601 FORMAT('0 * * * WARNING. NOT ENOUGH SPACE IN ENERGY() TO PROCESS 1TEMP(',I3,' ) =',F8.2) GO TO 1000 1010 XT=XK*T(I) WRITE(6,602) T(I),XT 602 FORMAT('0 FOR TEMP =',F8.2,' DEG. K =',F8.2,' (1/CM), THE 1AVERAGE IS APPROXIMATELY THE SUM OF') DO 1100 J=1,NGP EN=XT*A(IST+J) WT=A(IST+J)*W(IST+J) NN=NN+1 E(NN)=EN 1100 WRITE(6,603) WT,EN 603 FORMAT(15X,F13.8, ' * SIG( E =',F12.4,' ) ') 1000 CONTINUE NNRG=MIN0(MXNRG,MAX0(NNRG,NN)) RETURN END FUNCTION EPSUM(R,N,E,EVEC,WKS) C C FUNCTION TO EVALUATE THE EIGENPHASE SUM FROM THE R-MATRIX. C F02ABF DIAGONALISES THE N X N REAL SYMMETRIC R-MATRIX, C RETURNING THE EIGENVALUES IN E. C THE EIGENPHASE SUM IS THEN OBTAINED BY SUMMING ARCTANGENTS C OF THE EIGENVALUES. C THE RESULT IS RETURNED IN UNITS OF PI, AND IS SHIFTED TO BE C AS CLOSE AS POSSIBLE TO THE PREVIOUS EIGENPHASE SUM CALCULATED C (STORED IN SMLAST) C SEE ASHTON, CHILD AND HUTSON, J. CHEM. PHYS. 78, 4025 (1983). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE SMLAST DIMENSION R(N,N), E(N), EVEC(N,N), WKS(N) DATA PI/3.141592653589793238462643D0/ DATA SMLAST/10.D0/ C IF(N.EQ.1) GOTO 200 IFAIL=0 CALL F02ABF(R,N,N,E,EVEC,N,WKS,IFAIL) EPSUM=0.D0 DO 100 I=1,N X=ATAN(E(I)) EPSUM=EPSUM+X 100 CONTINUE GOTO 300 200 EPSUM=ATAN(R(1,1)) 300 EPSUM=EPSUM/PI DELTA=SMLAST-EPSUM+0.5D0 IF(DELTA.LE.0.D0) DELTA=DELTA-1.D0 IDEL=INT(DELTA) EPSUM=EPSUM+DBLE(IDEL) SMLAST=EPSUM RETURN END SUBROUTINE FINDRX(ENERGY,EINT,CENT,NNRG,N,CINT,RMAX,RSTOP, 1 NOPMAX,IRXSET,IPRINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C SUBROUTINE TO SCAN INPUT ENERGIES AND THRESHOLDS TO DETERMINE C A SAFE RMAX WHICH IS OUTSIDE THE CENTRIFUGAL BARRIER FOR C ALL COMBINATIONS. ALSO FIND THE LARGEST VALUE OF NOPEN C TO SAFEGUARD AGAINST SHRINKING THE BASIS SET TOO FAR. C DIMENSION ENERGY(NNRG),EINT(N),CENT(N) C NOPMAX=0 RSTOP=RMAX DO 200 J=1,NNRG NOPEN=0 ERED=ENERGY(J)*CINT DO 100 I=1,N DIF=ERED-EINT(I) IF(DIF.LT.0.D0) GOTO 100 NOPEN=NOPEN+1 IF(IRXSET.LE.0) GOTO 100 RCENT=SQRT(CENT(I)/DIF) RSTOP=MAX(RSTOP,RCENT) 100 CONTINUE 200 NOPMAX=MAX0(NOPMAX,NOPEN) IF(RSTOP.GT.RMAX .AND. IPRINT.GE.3) WRITE(6,601) RSTOP 601 FORMAT('0 RMAX INCREASED TO',F7.2,' FOR THIS PARITY CASE', 1 ' TO ENSURE THAT OPEN CHANNEL MATCHING'/' OCCURS BEYOND', 2 ' THE CENTRIFUGAL BARRIER FOR ALL ENERGIES') RETURN END SUBROUTINE GET102(MXLVL,NLEVEL,JLEVEL,ELEVEL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION JLEVEL(2,MXLVL),ELEVEL(MXLVL) WRITE(6,699) 699 FORMAT('0 GET102. DUMMY ROUTINE CALLED. TERMINAL ERROR.') STOP END SUBROUTINE IOSBIN(NVC,ITYPX,ATAU,MX,IASYMU,IPHIFX,IOSNG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE C C *** MODIFIED MAY 92 TO HANDLE ITYPE=103 (SG) *** C *** MODIFIED JAN/FEB 92 TO GENERALIZE ITYPE=2 HANDLING C *** MODIFIED FEB 88 TO CORRECT 'IHOMO' HANDLING OF ITYPE=5,6 CASES C WHERE POTENTIAL IS SYMMETRIC ABOUT THETA=PI/2 C *** AUG 86 ADD LM,LMAX ARGUMENTS TO IXQLF C *** UPDATED APR 86 TO MERGE BASIN-IOSBIN PROCESSING C C THIS IS IOS 'BASIS' ROUTINE FOR COMBINED MOLSCAT/IOS APR 86 C MODIFICATIONS MAY 1978 FOR ITYPE=5. C MODIFICATIONS SEPT 1985 FOR ITYPE=6, C INCLUDING ADDITION OF MMAX TO &INPUT. C C-----ENTRY IOSBIN READS &BASIS AND SETS UP BASIS DATA. C PARAMETERS ARE NVC (NO. VIB. CHANNELS), ITYPX RETURNS ROTOR TYPE C AND ATAU(MX) WHICH WILL HOLD ROTOR COEFF. FOR ITYPE=6 C DIMENSION ATAU(MX) C C-----ENTRY IOSBGP GETS LAM(MXLAM) INFORMATION FROM &POTL. C CAN THEN CHOOSE NGPT, LMAX, MMAX, AND SET UP GAUSS PTS/WTS C SPECIFICATIONS FOR ENTRY IOSBGP . . . DIMENSION LAM(MXLAM) LOGICAL ODD C C-----ENTRY IOSB1, CALLED AFTER STORAGE IS ALLOCATED, SETS UP PWGHT, VLI C ALSO IXQL AND LM. C SPECIFICATIONS . . . DIMENSION PWGHT(NGPT,LMAX), VLI(NGPT,MXXXXL) DIMENSION IXQL(NIXQL,NQL),LM(3,LMAX) C BELOW (TEMPORARY) TO CONTROL FLOW OF ALTERNATE ITYPE=3 CODE LOGICAL LNEW C C-----ENTRY IOSB2 IS CALLED JUST BEFORE INTEGRATOR. SETS UP VL, ETC. C SPECIFICATIONS FOR ENTRY IOSB2 . . . DIMENSION CENT(NVC),EINT(NVC),WVEC(NVC),VL(2),IVIX(2) DIMENSION LORB(NVC),JJJ(NVC),NB(NVC) C COMMON TO PASS ANGLES TO VRTP FOR "UNEXPANDED" (MXLAM=0) POTL CASE C N.B. 3RD ANGLE FOR ITYPE=3. IH0,IC0 TO SET IHOMO,ICNSYM IN VRTP C FACTOR IS 1./(VALUE OF LOWEST ANGULAR TERM) - DEPENDS ON ITYPE COMMON/ANGLES/COSANG(3),FACTOR,IH0,IC0 LOGICAL LVRTP C C-----ENTRY IXQLF RETURNS INDEX IN IXQL OF AN INPUT L,M1,M2,ICDE SYM. C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C TO CONTROL DEBUGGING OUTPUT OF COUPLING MATRIX C LDEBUG=.TRUE. CAN GIVE QUITE A BIT OF OUTPUT ! LOGICAL LDEBUG C C SPECIFICATIONS FOR MOLSCAT(&BASIS) COMPATIBILITY. . . DIMENSION ROTI(10),ALPHAE(2),BE(2),DE(2),A(2),B(2),C(2),WE(2), 1 WEXE(2),WT(2),ELEVEL(200) INTEGER JMIN,JMAX,NLEVEL,JLEVEL(400),J1MIN,J1MAX,J2MIN,J2MAX, 1 IDENT,JSTEP,J1STEP,J2STEP EQUIVALENCE (ROTI(1),BE(1),A(1)), (ROTI(3),ALPHAE(1),B(1)), 1 (ROTI(5),DE(1),C(1)),(JMIN,J1MIN),(JMAX,J1MAX), 2 (JSTEP,J1STEP), (ROTI(7),WE(1)),(ROTI(9),WEXE(1)) C COMMON /CMBASE/ ROTI,ELEVEL,EMAX,WT,SPNUC,JMIN,JMAX, 1 J2MIN,J2MAX,JSTEP,J2STEP,NLEVEL,JLEVEL,IDENT C C INTERNAL VERSION OF JLEVEL,ELEVEL IS ALSO USED . . . DIMENSION LEVV(200),EV(200) C C COMMON BLOCK TO COMMUNICATE WITH IOSOUT . . . COMMON /IOUTCM/ MAX,LEVV C C DYNAMIC STORAGE COMMON BLOCK ... COMMON /MEMORY/ MXXX,IXNEXT,NIPR,IVLFL,X(1) C C IOS GAUSS POINT CONTROL DIMENSION IOSNGP(3),IOSNG(3) C C ****************************************************************** C ** PROGRAM LIMITATION ** C ** DIMENSIONS FOR GAUSS POINTS ** C ** ------- ---------- ** C ** COULD USE /MEMORY/ FOR DYNAMIC STORAGE, BUT BELOW SHOULD ** C ** SUFFICE FOR MOST FEASIBLE CALCULATIONS. ** C ****************************************************************** DIMENSION COSA(400),GWT(400) C FOR ITYPE=3 TO HOLD PLM(LI,M,) AND COS(M*PHI) -- LNEW=.TRUE. CODE DIMENSION PL1(400),PL2(400),COSM(400) DATA MXGPT/400/ C C LIMIT FROM /CMBASE/ ELEVEL(MXLVL),JLEVEL(2*MXLVL) -- DATA MXLVL/200/ C AND EQUIVALENT INTERNAL ARRAYS DATA NVCMX/200/ DATA IZ/0/ DATA LNEW/.TRUE./,LDEBUG/.FALSE./ C C STATEMENT FUNCTION USED IN DETERMINING ITYPE=5,6 IHOMO SYMMETRY ODD(I,J)=(I-J)-2*((I-J)/2) .NE. 0 C C C SPECIFICATIONS FOR LEGENDRE STATEMENT FUNCTION . . . XLEG(I,TH)=SQRT(2.D0/DBLE(2*I+1))*PLM(I,0,TH) C N.B. PLM(L,M,COSTH) RETURNS A **NORMALIZED** ASSOC. LEG. POLY. C C NB. THE FOLLOWING VARIABLES ARE USED AS LIMITS (SOME ITYPE=5 ONLY) C MAX=HIGHEST J IN BASIS / MXK=HIGHEST K (SYM. TOP) IN BASIS. C LMBDMX=HIGHEST LAMBDA IN POTL / MUMX=HIGHEST MU IN POTL C LMMAX=HIGHEST 'L' IN SLLR,SLLI,QLT / MUMAX=HIGHEST 'M' C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C PI=ACOS(-1.D0) WRITE(6,666) 666 FORMAT('0 PROCESSED BY IOSBIN ROUTINE (FEB/MAY 92).', 1 ' MODIFIED/NEW ITYPE 102/103.') C C SET LOCAL (AND HENCE KEPT) VALUES FROM ARGUMENTS DO 1107 I=1,3 1107 IOSNGP(I)=IOSNG(I) IPHIFL=IPHIFX C INITIALIZE SIZE OF ATAU TO ZERO. MXA=MX MX=0 C SET DEFAULT IH0,IC0 WHICH MAY BE CHANGED IN VRTP IH0=0 IC0=0 C WRITE(6,620) ITYPX 620 FORMAT('0 INPUT ITYPE =',I4) C ITYP=ITYPX-10*(ITYPX/10) ITYPX=100+ITYP IF (ITYP.EQ.1) GO TO 1000 IF (ITYP.EQ.2) GO TO 2000 IF (ITYP.EQ.3) GO TO 3000 IF (ITYP.EQ.5) GO TO 5000 IF (ITYP.EQ.6) GO TO 6000 WRITE(6,699) ITYP 699 FORMAT('0 * * * ERROR. MOD(ITYPE,10) =',I3,' NOT SUPPORTED.') STOP C 1000 NVC=1 ASSIGN 6100 TO IGOTP EV(1)=0.D0 LEVV(1)=0 ILOFF=1 IF (NLEVEL.GT.0) GO TO 1200 WRITE(6,601) JMIN,JMAX,JSTEP 601 FORMAT('0 JLEVEL, NLEVEL CREATED FROM JMIN, JMAX, JSTEP =', & 3I5) JMIN=MAX0(JMIN,0) JMAX=MAX0(JMIN,JMAX) JSTEP=MAX0(JSTEP,1) MAX=0 NLEVEL=0 DO 1100 I=JMIN,JMAX,JSTEP IF (NLEVEL.GE.MXLVL) GO TO 1109 NLEVEL=NLEVEL+1 JLEVEL(NLEVEL)=I 1100 MAX=MAX0(MAX,I) GO TO 9000 1109 WRITE(6,698) 698 FORMAT('0 * * * WARNING. OUT OF SPACE IN JLEVEL. ', 1 'BASIS TRUNCATED.') NLEVEL=MXLVL GO TO 9000 1200 WRITE(6,602) NLEVEL,(JLEVEL(I),I=1,NLEVEL) 602 FORMAT('0 BASIS TAKEN FROM NLEVEL, JLEVEL INPUT. NO. OF LEVELS ', & '(NLEVEL) =',I4/(' ',20I5) ) MAX=0 DO 1300 I=1,NLEVEL IF (I.GT.MXLVL) GO TO 1109 1300 MAX=MAX0(MAX,JLEVEL(I)) GO TO 9000 C C>>SG --------- CODE REWORKED FEB 92 (ANTICIPATING H2-H CALCULATIONS) 2000 ILOFF=3 ASSIGN 6200 TO IGOTP C IF (NLEVEL.GT.0) GO TO 2100 C -------- GET 'VIBRATIONAL LEVELS' FROM SPECIAL SUBROUTINE ------- WRITE(6,697) 697 FORMAT('0 IOSBIN (FEB 92). NLEVEL.LT.0. AN APPROPRIATE', 1 ' SUBROUTINE MUST BE PROVIDED:'/ 2 35X,'GET102(MXLVL,NLEVEL,JLEVEL,ELEVEL)'/) C N.B. NLEVEL,JLEVEL,ELEVEL ARE AVAILABLE IN /CMBASE/ BUT MXLVL C IS A LOCAL DATUM GEARED TO DIMENSIONS IN COMMON CALL GET102(MXLVL,NLEVEL,JLEVEL,ELEVEL) NVC=NLEVEL MAX=0 DO 2190 I=1,NVC EV(I)=ELEVEL(I) LEVV(I)=JLEVEL(2*I) 2190 MAX=MAX0(MAX,JLEVEL(2*I-1)) C SKIP '2009' PRINT OUT OF LEVEL INFO/ DO IT IN GET102 IF DESIRED GO TO 9000 C C ----- GET 'VIBRATIONAL LEVELS' FROM JLEVEL ----- C CURRENT CODE DOES NOT ALLOW DUPLICATE VIB LEVELS. 2100 NVC=1 ITOP=2*NLEVEL WRITE(6,602) NLEVEL,(JLEVEL(I),I=1,ITOP) LEVV(1)=JLEVEL(2) MAX=JLEVEL(1) I=2 2102 IF (I.GT.NLEVEL) GO TO 2110 DO 2103 II=1,NVC IF (LEVV(II).NE.JLEVEL(2*I)) GO TO 2103 WRITE(6,693) I,JLEVEL(2*I-1),JLEVEL(2*I) 693 FORMAT(' IOSBIN (FEB 92). LEVEL',I4,' V,J =',2I4,' DUPLICATES' 1 ,' AN EARLIER VIB LEVEL.'/ 2 20X,'VIBRATIONAL VALUE IGNORED, HIGHER J-VALUE KEPT.') JXX=MAX0(JLEVEL(2*II-1),JLEVEL(2*I-1)) JLEVEL(2*II-1)=JXX MAX=MAX0(MAX,JXX) IF (I.LT.NLEVEL) GO TO 2120 C I.EQ.NLEVEL ==> REDUCE NLEVEL AND GET OUT NLEVEL=NLEVEL-1 GO TO 2110 C PULL DOWN LIST/ DECREASE NLEVEL/ GO BACK FOR NEW I-TH LEVEL 2120 DO 2121 J=I+1,NLEVEL ELEVEL(J-1)=ELEVEL(J) JLEVEL(2*J-3)=JLEVEL(2*J-1) 2121 JLEVEL(2*J-2)=JLEVEL(2*J) NLEVEL=NLEVEL-1 GO TO 2102 2103 CONTINUE C DUPLICATE VIB LEVEL NOT FOUND/ ADD THIS VIBRATIONAL LEVEL IF (NVC.LE.NVCMX) GO TO 2104 WRITE(6,694) NVCMX 694 FORMAT('0 ISOBIN -- ERROR. VIBRATIONAL LEVELS IN NLEVEL/JLEVEL' 1 ,' EXCEED NVCMX =',I4) STOP 2104 NVC=NVC+1 LEVV(NVC)=JLEVEL(2*I) MAX=MAX0(MAX,JLEVEL(2*I-1)) I=I+1 GO TO 2102 C 2110 WRITE(6,692) NVC 692 FORMAT('0 IOSBIN (FEB 92). NUMBER OF VIB. CHANNELS (NVC) =',I4) C C ----- GET ENERGY LEVELS ----- DO 2111 I=1,NVC IF (ELEVEL(I).EQ.0.) GO TO 2111 C IF ELEVEL() VALUES ARE SET (NON-ZERO), USE THEM GO TO 2290 2111 CONTINUE C IF WE REACH HERE, ALL ELEVEL ARE ZERO. C IF THERE IS ONLY ONE LEVEL, AND ENERGY()=0, WE ARE STILL OKEY IF (NVC.GT.1) GO TO 2280 C SET EV() FROM ELEVEL() AND WE ARE DONE. 2290 WRITE(6,691) 691 FORMAT('0 IOSBIN (FEB 92). VIBRATIONAL ENERGIES ', 1 'TAKEN FROM ELEVEL INPUT.') DO 2291 I=1,NVC 2291 EV(I)=ELEVEL(I) GO TO 2009 C OTHERWISE, SEE IF WE CAN CALCULATE ENERGIES FROM WE, WEXE 2280 IF (WE(1).GT.0.D0) GO TO 2200 WRITE(6,696) NVC,WE(1) 696 FORMAT('0 IOSBIN (FEB 92) CANNOT GET ENERGIES FROM ELEVEL ', & 'OR WE. NVC, WE =',I6,D14.4) STOP 2200 WRITE(6,603) WE(1) 603 FORMAT('0 TARGET ENERGY LEVELS (TAKING V = 0 AS ZERO ENERGY)', 1 ' COMPUTED FROM WE =',F10.4) IF (WEXE(1).NE.0.D0) WRITE(6,604) WEXE(1) 604 FORMAT(67X,'CORRECTED FOR WEXE =',F10.6) DO 2201 I=1,NVC FV=LEVV(I) EV(I)=WE(1)*FV-WEXE(1)*FV*(FV+1.D0) C STORE BACK IN JLEVEL,ELEVEL FOR ISAVEU OUTPUT PURPOSES. C>>SG JLEVEL(I)=LEVV(I) -->> THIS WAS USED FOR ISAVEU CAPABILITY IN C>>SG -->> IOSDRV AND WILL NO LONGER WORK THERE C ELEVEL() IS NEEDED IN IOSB2 TO GET EINT, ETC. 2201 ELEVEL(I)=EV(I) NLEVEL=NVC C ------ OUTPUT LEVV, EV ------ 2009 DO 2019 I=1,NVC 2019 WRITE(6,613) I,LEVV(I),EV(I) 613 FORMAT(' LEVEL',I4,' LEVV =',I4,' EV =',F12.4) GO TO 9000 C<>SG ----- ITYPE=3 CODE ADDED 5/6/92 (SG) 3000 ILOFF=3 ASSIGN 6300 TO IGOTP NVC=1 LEVV(1)=0 EV(1)=0. IF (NLEVEL.GT.0) GO TO 3901 WRITE(6,632) JMIN,JMAX,JSTEP,J2MIN,J2MAX,J2STEP 632 FORMAT('0 ROTATIONAL LEVELS FROM JMIN JMAX JSTEP'/ 1 ' ROTOR 1 -- ',3I5/' ROTOR 2 -- ',3I5) GO TO 3903 3901 ITOP=2*NLEVEL WRITE(6,633) NLEVEL,(JLEVEL(I),I=1,ITOP) 633 FORMAT('0 ROTATIONAL LEVELS FROM NLEVEL =',I4,' -- JLEVEL ='/ 1 (25I4)) DO 3902 I=1,NLEVEL JMIN=MIN0(JMIN,JLEVEL(2*I-1)) JMAX=MAX0(JMAX,JLEVEL(2*I-1)) J2MIN=MIN0(J2MIN,JLEVEL(2*I)) 3902 J2MAX=MAX0(J2MAX,JLEVEL(2*I)) 3903 MAX=MAX0(JMAX,J2MAX) IF (IDENT.GT.0) WRITE(6,634) IDENT 634 FORMAT('0 IDENTICAL PARTICLES SPECIFIED BY IDENT =',I3) GO TO 9000 C<>SG IN CASE ONLY LAMBDA=0 TERMS APPEAR IN THE POTENTIAL (E.G., IN C A BREATHING SPHERE TYPE VIBRATIONAL CALC) THE CODE BELOW WILL C (ERRONEOUSLY) SET LVRTP=.TRUE. HOWEVER, SINCE THE POTENTIAL C IS THEN SPHERICALLY SYMMETRIC, THIS OUGHT TO STILL WORK. ONE C MIGHT, WORRY, HOWEVER, ABOUT IHOMO SETTING, AND DOUBLE CHECK C BEFORE RUNNING SUCH A CASE C<>SG 5/11/92 C>>SG 5/11/92 N.B. THIS SHOULD BE REDONE TO TAKE ADVANTAGE OF NEW GAUSSP C>>SG 5/11/92 3500 IF (.NOT.LVRTP) GO TO 3540 IHOMO=1 IF (IH0.EQ.0) GO TO 3543 IHOMO=IH0 GO TO 3543 C ABOVE ALLOWS SETTING IN VRTP ROUTINE/ BELOW CHECKS INPUT L,M SYMS. 3540 IHOMO=2 DO 3542 L=1,MXLAM IF (ODD(LAM(2*L),LAM(2*L-1))) IHOMO=1 3542 CONTINUE 3543 THLO=-1.D0 THHI=1.D0 SFACT=1.D0 IF (IHOMO.EQ.1) GO TO 3541 WRITE(6,618) 618 FORMAT('0 * * * NOTE.'/' * * * NOTE. USE WILL BE MADE OF FACT ', & 'THAT POTENTIAL IS SYMMETRIC ABOUT THETA=PI/2'/' * * * NOTE.') THLO=0.D0 SFACT=2.D0 C NEXT GET ICNSYM (OLD CODE SHOULD ALWAYS STILL WORK OKEY) 3541 LMBDMX=MXXXXL MUMX=0 DO 3501 L=1,MXLAM 3501 MUMX=MAX0(MUMX,IABS(LAM(2*L))) C FIND ICNSYM WHICH IS PHI EQUIVALENT OF IHOMO IF (MUMX.GT.1) GO TO 3502 C ALLOW SETTING OF ICNSYM FOR 'MXSYM=0' (UNEXPANDED POTL) CASE ICNSYM=1 IF (.NOT.LVRTP .OR. IC0.EQ.0) GO TO 3503 ICNSYM=IC0 WRITE(6,654) ICNSYM 654 FORMAT('0 * * * NOTE. ICNSYM TAKEN FROM VRTP ROUTINE =',I4) GO TO 3503 3502 ICNSYM=MUMX 3506 DO 3504 L=1,MXLAM M=IABS(LAM(2*L)) IF (M-(M/ICNSYM)*ICNSYM .NE. 0) GO TO 3505 3504 CONTINUE GO TO 3503 3505 ICNSYM=ICNSYM-1 IF (ICNSYM.GT.1) GO TO 3506 3503 PHILO=0.D0 PHIHI=PI/DBLE(ICNSYM) SFACT=SFACT*DBLE(ICNSYM)*2.D0 C N.B. WE USE HERE FACT THAT POTENTIAL IS EVEN IN PHI SO INTEGRAL C IS TWICE THAT FROM 0 TO PI. THIS IS REFLECTED IN HAVING ONLY C COS (M*PHI) AND NOT SIN (M*PHI) IN PWGHT, ETC. IF (ICNSYM.GT.1) WRITE(6,658) ICNSYM 658 FORMAT('0 * * * NOTE.'/' * * * NOTE. USE WILL BE MADE OF',I4 & ,'-FOLD SYMMETRY ABOUT Z-AXIS.'/' * * * NOTE.') C DETERMINE NO. OF LAMBDA, MU SYMMETRIES (MXXXXL) MXXXXL=0 DO 3507 L=IZ,LMBDMX MTOP=MIN0(MUMX,L) DO 3507 M=IZ,MTOP,ICNSYM IF (IHOMO.EQ.2 .AND. ODD(L,M)) GO TO 3507 MXXXXL=MXXXXL+1 3507 CONTINUE C DETERMINE NO. OF GAUSSPOINTS FOR THETH(NGL) AND PHI (NGM) C ** N.B. SOME MANEUVERING IS NECESSARY SINCE GAUSSP MAY REDUCE NPT. WRITE(6,645) IPHIFL 645 FORMAT('0 * * * NOTE. IPHIFL (PHI INTEGRATION FLAG) =',I4) IPASS=0 IF (IOSNGP(1).GT.0) GO TO 3510 NGL=2*(MAX+LMBDMX) WRITE(6,656) NGL,MAX,LMBDMX 656 FORMAT('0 * * * NOTE. INPUT IOSNGP(1) .LE. 0 (DEFAULT). ', & 'NGL =',I4,' COMPUTED FROM MAX, LMBDMX =',2I4) GO TO 3531 3510 NGL=IOSNGP(1) WRITE(6,655) NGL,IOSNGP(1) 655 FORMAT('0 * * * NOTE. NGL =',I4,' TAKEN FROM', & ' &BASIS IOSNGP(1) =',I4) 3531 IF (IOSNGP(2).GT.0) GO TO 3532 NGM=MAX0(1,2*(MUMX+MXK)) WRITE(6,647) NGM,MXK,MUMX 647 FORMAT('0 * * * NOTE. INPUT IOSNGP(2) .LE. 0 (DEFAULT). NGM =', & I4, ' COMPUTED FROM MXK, MUMX =',2I4) GO TO 3511 3532 WRITE(6,646) IOSNGP(2) 646 FORMAT('0 * * * NOTE. NGM SET FROM &BASIS IOSNGP(2) =',I4) NGM=MAX0(1,IOSNGP(2)) 3511 CALL GAUSSP(THLO,THHI,NGL,COSA,GWT) IF (MUMX.GT.0 .OR. IPASS.GT.0 .OR.LVRTP) GO TO 3512 NGM=1 WRITE(6,650) 650 FORMAT('0 * * * NOTE.'/' * * * NOTE. POTENTIAL HAS NO PHI ', & 'DEPENDENCE. INTEGRAL DONE ANALYTICALLY.') 3512 IF (IPHIFL.NE.0) &CALL GAUSSP(PHILO,PHIHI,NGM,COSA,GWT) NGPT=NGM+NGL IF (NGPT.LE.MXGPT) GO TO 3513 WRITE(6,607) NGPT,MXGPT NGL=(DBLE(MXGPT)/DBLE(NGPT))*NGL NGM=(DBLE(MXGPT)/DBLE(NGPT))*NGM NGM=MAX0(NGM,1) IPASS=1 GO TO 3511 3513 CALL GAUSSP(THLO,THHI,NGL,COSA,GWT) WRITE(6,609) NGL,(COSA(I),GWT(I),I=1,NGL) IF (IPHIFL.EQ.0) GO TO 3515 CALL GAUSSP(PHILO,PHIHI,NGM,COSA(NGL+1),GWT(NGL+1)) WRITE(6,644) NGM,(COSA(NGL+I),GWT(NGL+I),I=1,NGM) 644 FORMAT('0 PHI INTEGRATION DONE BY ',I3,'-POINT GAUSSIAN ', & 'QUADRATURE. POINTS/WEIGHTS ARE'/(4(10X,2F10.6))) GO TO 3516 3515 IX=NGL FACTL=(PHIHI-PHILO)/DBLE(NGM) TH=-FACTL/2.D0 DO 3514 I=1,NGM IX=IX+1 TH=TH+FACTL GWT(IX)=FACTL 3514 COSA(IX)=TH WRITE(6,651) NGM,(COSA(NGL+I),GWT(NGL+I),I=1,NGM) 651 FORMAT('0 PHI INTEGRATION DONE BY ',I3,'-POINT GAUSS-MEHLER ', & 'CHEBYSCHEV) QUADRATURE. POINTS/WEIGHTS ARE'/(4(10X,2F10.6))) 3516 WRITE(6,653) SFACT 653 FORMAT('0 ABOVE WEIGHTS MULTIPLIED BY SYMMETRY FACTOR =',D16.8) NGPT=NGL*NGM C NEXT CHOOSE LMAX (LMMAX,MUMAX) IF (LMAX.LE.0) GO TO 3520 WRITE(6,610) LMAX LMMAX=LMAX GO TO 3523 3520 LMMAX=MIN0(NGL*IHOMO,2*MAX) WRITE(6,612) LMMAX C INPUT CAPABILITY ON MMAX ADDED IN VERSION 6. 3523 IF (MMAX.LE.0) GO TO 3525 MUMAX=MMAX WRITE(6,630) MMAX 630 FORMAT('0 MMAX TAKEN FROM &INPUT MMAX =',I5) GO TO 3524 3525 MUMAX=MIN0(NGM*ICNSYM,2*MXK,LMMAX) WRITE(6,631) MUMAX 631 FORMAT('0 * * * WARNING. MMAX=0 (DEFAULT). WILL USE HIGHEST', & ' VALUE CONSISTENT WITH IOSNGP(2), MMAX =',I4) C RESET LMAX TO REFLECT *NUMBER* OF LAMBDA,MU VALUES. C AND COUNT NQL (NUMBER OF QLT VALUES) 3524 LMAX=0 NQL=0 DO 3521 L=IZ,LMMAX MTOP=MIN0(MUMAX,L) DO 3521 M=IZ,MTOP,ICNSYM IF (IHOMO.EQ.2 .AND. ODD(L,M)) GO TO 3521 LMAX=LMAX+1 DO 3522 I=IZ,M,ICNSYM IF (IHOMO.EQ.2 .AND. ODD(L,I)) GO TO 3522 IX=2 IF (I.EQ.M) IX=1 NQL=NQL+IX 3522 CONTINUE 3521 CONTINUE NIXQL=3 RETURN C C ITYPE=3 CODE ADDED 5/6/92 (SG) C GET POTENL SYMS. (ITYPE=3 USES IHOMO,ICNSYM FOR IHOMO1,IHOMO2) 3330 L1MAX=0 L2MAX=0 LLMAX=0 IL=0 DO 3331 I=1,MXLAM L1MAX=MAX0(L1MAX,LAM(IL+1)) L2MAX=MAX0(L2MAX,LAM(IL+2)) LLMAX=MAX0(LLMAX,LAM(IL+3)) 3331 IL=IL+ILOFF MXXXXL=MAX0(MXXXXL,L1MAX,L2MAX,LLMAX) LVRTP=MXXXXL.LE.0 IF (.NOT.LVRTP) GO TO 3332 C ?? MXXXXL=1 THIS WILL BE TAKEN CARE OF IN 3336 LOOP IHOMO=1 ICNSYM=1 IF (IH0.EQ.0) GO TO 3333 IHOMO=IH0 WRITE(6,637) IHOMO 637 FORMAT('0 * * * NOTE. IHOMO (MOL 1) TAKEN FROM VRTP ROUTINE =',I4) 3333 IF (IC0.EQ.0) GO TO 3334 ICNSYM=IC0 WRITE(6,638) ICNSYM 638 FORMAT('0 * * * NOTE. IHOMO (MOL 2) TAKEN FROM VRTP ROUTINE =',I4) GO TO 3334 C FOR EXPANDED POTENTIAL (.NOT.LVRTP) GET IHOMO1,IHOMO2 FROM LAM 3332 IHOMO=2 ICNSYM=2 DO 3335 I=1,MXLAM IF (ODD(LAM(3*I-2),0)) IHOMO=1 3335 IF (ODD(LAM(3*I-1),0)) ICNSYM=1 3334 IM=1 IF (IHOMO.EQ.2) WRITE(6,639) IM 639 FORMAT('0 * * * NOTE.'/' * * * NOTE. USE WILL BE MADE OF FACT ', & 'THAT POTENTIAL IS SYMMETRIC ABOUT PI/2 FOR MOLECULE',I3) IM=2 IF (ICNSYM.EQ.2) WRITE(6,639) IM C>>SG (5/18/92) STORE IHOMO,ICNSYM BACK IN IH0,IC0 FOR USE IN IOSOUT IH0=IHOMO IC0=ICNSYM C COUNT L1,L2,LL SYMMETRIES (MXXXXL) C FOR IDENT PARTICLES, L1,L2<->L2,L1 MUST BOTH BE IN POTL SYMS MXXXXL=0 DO 3336 L1=IZ,L1MAX,IHOMO L2TOP=L2MAX IF (IDENT.GT.0) L2TOP=L1MAX DO 3336 L2=IZ,L2TOP,ICNSYM LLO=ABS(L1-L2) LHI=L1+L2 DO 3336 LL=LLO,LHI IF (ODD(L1+L2,LL)) GO TO 3336 MXXXXL=MXXXXL+1 3336 CONTINUE C SET INTEGRATION LIMITS AND GET GAUSS POINTS C CURRENT GAUSSP (5/6/92) DOES *ARBITRARY* NO PTS; C IF REQUEST EXCEEDS DIMENSIONS (MXGPT) TERMINATE. SFACT=1.D0 IF (IOSNGP(1)*IOSNGP(2)*IOSNGP(3).GT.0 .AND. 1 IOSNGP(1)+IOSNGP(2)+IOSNGP(3).LE.MXGPT) GO TO 3337 WRITE(6,636) IOSNGP,MXGPT 636 FORMAT('0 IOSBGP. ERROR. IOSNGP INPUT, ',3I5,' ILLEGAL OR ', 1 'EXCEEDS STORAGE (MXGPT) =',I5) STOP 3337 NGP1=IOSNGP(1) THLO=-1.D0 THHI=1.D0 IF (IHOMO.EQ.2) THEN THLO=0.D0 SFACT=SFACT*2.D0 ENDIF CALL GAUSSP(THLO,THHI,NGP1,COSA(1),GWT(1)) WRITE(6,609) NGP1,(COSA(I),GWT(I),I=1,NGP1) NGP2=IOSNGP(2) IST2=NGP1 THLO=-1.D0 IF (ICNSYM.EQ.2) THEN THLO=0.D0 SFACT=SFACT*2.D0 ENDIF CALL GAUSSP(THLO,THHI,NGP2,COSA(IST2+1),GWT(IST2+1)) WRITE(6,609) NGP2,(COSA(IST2+I),GWT(IST2+I),I=1,NGP2) WRITE(6,645) IPHIFL PHILO=0.D0 C CAN ALWAYS USE SYMMETRY V(-PHI)=V(PHI) TO REDUCE INTEGRAL C FROM (0,2*PI) TO (0,PI) -- CORRECT SFACT ACCORDINGLY PHIHI=PI SFACT=SFACT*2.D0 IST3=IST2+NGP2 NGM=IOSNGP(3) IF (IPHIFL.EQ.0) GO TO 3338 CALL GAUSSP(PHILO,PHIHI,NGM,COSA(IST3+1),GWT(IST3+1)) WRITE(6,644) NGM,(COSA(IST3+I),GWT(IST3+I),I=1,NGM) GO TO 3339 3338 IX=IST3 FACTL=(PHIHI-PHILO)/DBLE(NGM) TH=-FACTL/2.D0 DO 3342 I=1,NGM IX=IX+1 TH=TH+FACTL GWT(IX)=FACTL 3342 COSA(IX)=TH WRITE(6,651) NGM,(COSA(IST3+I),GWT(IST3+I),I=1,NGM) C SET NGPT AS PRODUCT OF THETA-1, THETA-2, PHI GRIDS 3339 NGPT=NGP1*NGP2*NGM WRITE(6,653) SFACT C RESET LMAX AND NQL=LMAX TO REFLECT NUMBER OF L1,L2,LL VALUES C USE ITYP=5,6 VARIABLES: LMMAX FOR L1, MUMAX FOR L2 IF (LMAX.GT.0) GO TO 3340 LMMAX=(NGP1-1)*IHOMO WRITE(6,640) LMAX,LMMAX,NGP1,IHOMO 640 FORMAT('0 &INPUT LMAX =',I4,' -- L1MAX =',I4,' CALCULATED FROM ', 1 ' NGP1 AND (SYMMETRY) IHOMO =',2I4) GO TO 3344 3340 LMMAX=LMAX WRITE(6,641) LMAX 641 FORMAT(' L1MAX TAKEN FROM &INPUT LMAX =',I4) 3344 IF (MMAX.GT.0) GO TO 3343 MUMAX=(NGP2-1)*ICNSYM WRITE(6,642) MMAX,MUMAX,NGP2,ICNSYM 642 FORMAT('0 &INPUT MMAX =',I4,' -- L2MAX =',I4,' CALCULATED FROM ', 1 ' NGP2 AND (SYMMETRY) ICNSYM=',2I4) GO TO 3345 3343 MUMAX=MMAX WRITE(6,643) MMAX 643 FORMAT(' L2MAX TAKEN FROM &INPUT MMAX =',I4) 3345 LMAX=0 DO 3341 L1=IZ,LMMAX,IHOMO L2TOP=MUMAX C IDENTICAL PARTICLES KEEP ONLY L1.GE.L2 IN LM(,) IF (IDENT.GT.0) L2TOP=L1 DO 3341 L2=IZ,L2TOP,ICNSYM LLO=ABS(L1-L2) LHI=L1+L2 DO 3341 LL=LLO,LHI IF (ODD(L1+L2,LL)) GO TO 3341 LMAX=LMAX+1 3341 CONTINUE NQL=LMAX NIXQL=2 RETURN C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ENTRY IOSB1(PWGHT,VLI,IXQL,LM,NGPT,LMAX,MXXXXL,NIXQL,NQL) C IF (ITYP.EQ.3) GO TO 4300 IF (ITYP.EQ.5 .OR. ITYP.EQ.6) GO TO 4500 C C N.B. FOR ITYPE=1,2 LMAX=NQL C PWGHT MULTIPLY SL(COS(THETA)) TO GET LEGENDRE COEFFICIENTS FACTL=-.5D0 DO 4100 L=1,LMAX IXQL(1,L)=L IXQL(2,L)=0 LM1=L-1 LM(1,L)=LM1 FACTL=FACTL+1.D0 DO 4101 NX=1,NGPT C N.B. WE KEEP EVEN AND ODD L FOR HOMONUCLEARS, BUT SET TO 0. IF NEC PWGHT(NX,L)=0.D0 IF (IHOMO.EQ.2 .AND. L-2*(L/2).EQ.0) GO TO 4101 PWGHT(NX,L)=FACTL*GWT(NX)*XLEG(L-1,COSA(NX)) 4101 CONTINUE 4100 CONTINUE C NEXT COMPUTE VLI DO 4200 NX=1,NGPT L=0 DO 4201 IL=1,MXXXXL VLI(NX,IL)=XLEG(L,COSA(NX)) 4201 L=L+IHOMO 4200 CONTINUE GO TO 4999 C C ITYPE=3 -- SETUP VLI 4300 I=0 IF (LNEW) GO TO 4993 C>>SG 5/21/92 BELOW IS OLD CODE - BYPASSED FOR LNEW=.TRUE. DO 4301 IX1=1,NGP1 DO 4301 IX2=1,NGP2 DO 4301 IX3=1,NGM C I COUNTS GAUSS POINTS TO NGPT. I=I+1 IL=0 DO 4301 L1=IZ,L1MAX,IHOMO L2TOP=L2MAX IF (IDENT.GT.0) L2TOP=L1MAX DO 4301 L2=IZ,L2TOP,ICNSYM LLO=ABS(L1-L2) LHI=L1+L2 DO 4301 LL=LLO,LHI IF (ODD(L1+L2,LL)) GO TO 4301 C IL COUNTS SYMMETRIES IN POTENTIAL TO L1MAX,L2MAX IL=IL+1 VLI(I,IL)=YRR(L1,L2,LL,COSA(IX1),COSA(IST2+IX2),COSA(IST3+IX3)) 4301 CONTINUE PIFACT=2.D0*PI*SFACT IL=0 DO 4302 L1=IZ,LMMAX,IHOMO L2TOP=MUMAX IF (IDENT.GT.0) L2TOP=L1 DO 4302 L2=IZ,L2TOP,ICNSYM LLO=ABS(L1-L2) LHI=L1+L2 DO 4302 LL=LLO,LHI IF (ODD(L1+L2,LL)) GO TO 4302 XLFACT=1.D0/(2.D0*LL+1.D0) IL=IL+1 IXQL(1,IL)=IL IXQL(2,IL)=0 LM(1,IL)=L1 LM(2,IL)=L2 LM(3,IL)=LL I=0 DO 4303 IX1=1,NGP1 DO 4303 IX2=1,NGP2 DO 4303 IX3=1,NGM I=I+1 4303 PWGHT(I,IL)=GWT(IX1)*GWT(IST2+IX2)*GWT(IST3+IX3)* 1 YRR(L1,L2,LL,COSA(IX1),COSA(IST2+IX2),COSA(IST3+IX3)) 2 *PIFACT*XLFACT 4302 CONTINUE GO TO 4998 C>>SG 5/21/92 ------- END OF OLD CODE C C NEW CODE 5/21/92 MUCH MORE EFFICIENT. YRR() ASSEMBLED AS NEEDED C AVOIDING RECALCULATION OF THRJ, PLM, ETC. 4993 DEN=SQRT(4.D0*PI)*2.D0*PI DO 4310 IL=1,MXXXXL DO 4310 IX=1,NGPT 4310 VLI(IX,IL)=0.D0 MTOP=MIN0(L1MAX,L2MAX) DO 4311 M=IZ,MTOP PTM=PARITY(M) XM=M DO 4312 IX=1,NGM COSM(IX)=COS(XM*COSA(IST3+IX))/DEN IF (M.EQ.0) GO TO 4312 COSM(IX)=COSM(IX)*(2.D0*PTM) 4312 CONTINUE IL=0 DO 4313 L1=IZ,L1MAX,IHOMO IF (L1.LT.M) GO TO 4317 XL1=L1 PTL1=PARITY(L1) DO 4314 IX=1,NGP1 4314 PL1(IX)=PLM(L1,M,COSA(IX))*PTL1 4317 L2TOP=L2MAX IF (IDENT.NE.0) L2TOP=L1MAX DO 4313 L2=IZ,L2TOP,ICNSYM IF (L2.LT.M) GO TO 4318 XL2=L2 PTL2=PARITY(L2) DO 4315 IX=1,NGP2 4315 PL2(IX)=PLM(L2,M,COSA(IST2+IX))*PTL2 4318 LLO=ABS(L1-L2) LHI=L1+L2 DO 4313 LL=LLO,LHI IF (ODD(L1+L2,LL)) GO TO 4313 IL=IL+1 IF (L1.LT.M .OR. L2.LT.M) GO TO 4313 XL=LL TJ=THRJ(XL1,XL2,XL,XM,-XM,0.D0)*(2.D0*XL+1.D0) I=0 DO 4316 IX1=1,NGP1 DO 4316 IX2=1,NGP2 DO 4316 IX3=1,NGM I=I+1 4316 VLI(I,IL)=VLI(I,IL)+PL1(IX1)*PL2(IX2)*COSM(IX3)*TJ 4313 CONTINUE 4311 CONTINUE C C NOW SET UP IXQL, LM, AND PWGHT C N.B. PIFACT IS CONSISTENT WITH AGG & CLARY, EQS. (19)-(20) C AND W/ GOLDFLAM & KOURI, EQS. (68), (69), (89), (121). C I.E. T(ANGLES)=(4*PI)*SUM(L1,L2,L) T(L1,L2,L)*YRR(L1,L2,L/ANGLES) C NEW CODE 5/21/92 MUCH MORE EFFICIENT; CALC YRR() LOCALLY PIFACT=2.D0*PI*SFACT DO 4320 IL=1,LMAX DO 4320 IX=1,NGPT 4320 PWGHT(IX,IL)=0.D0 MTOP=MIN0(LMMAX,MUMAX) IF (IDENT.GT.0) MTOP=LMMAX DO 4321 M=0,MTOP PTM=PARITY(M) XM=M DO 4322 IX=1,NGM COSM(IX)=COS(XM*COSA(IST3+IX))/DEN IF (M.EQ.0) GO TO 4322 COSM(IX)=COSM(IX)*(2.D0*PTM) 4322 CONTINUE IL=0 DO 4323 L1=IZ,LMMAX,IHOMO IF (L1.LT.M) GO TO 4324 XL1=L1 PTL1=PARITY(L1) DO 4325 IX=1,NGP1 4325 PL1(IX)=PLM(L1,M,COSA(IX))*PTL1 4324 L2TOP=MUMAX IF (IDENT.GT.0) L2TOP=L1 DO 4323 L2=IZ,L2TOP,ICNSYM IF (L2.LT.M) GO TO 4326 XL2=L2 PTL2=PARITY(L2) DO 4327 IX=1,NGP2 4327 PL2(IX)=PLM(L2,M,COSA(IST2+IX))*PTL2 4326 LLO=ABS(L1-L2) LHI=L1+L2 DO 4323 LL=LLO,LHI IF (ODD(L1+L2,LL)) GO TO 4323 IL=IL+1 C STORE IXQL, LM ONLY FOR M=0 PASS ONLY. IF (M.GT.0) GO TO 4328 IXQL(1,IL)=IL IXQL(2,IL)=0 LM(1,IL)=L1 LM(2,IL)=L2 LM(3,IL)=LL 4328 IF (L1.LT.M .OR. L2.LT.M) GO TO 4323 XL=LL C TJ=THRJ(XL1,XL2,XL,XM,-XM,0.D0)*PIFACT/(2.D0*XL+1.D0) C 2*L+1 FACTOR CANCELS THAT IN DEF OF YRR ??? TJ=THRJ(XL1,XL2,XL,XM,-XM,0.D0)*PIFACT I=0 DO 4329 IX1=1,NGP1 DO 4329 IX2=1,NGP2 DO 4329 IX3=1,NGM I=I+1 4329 PWGHT(I,IL)=PWGHT(I,IL)+PL1(IX1)*PL2(IX2)*COSM(IX3)*TJ 4323 CONTINUE 4321 CONTINUE C END OF M-LOOP - PWGHT NOW CONTAINS YRR; NEED TO MULT BY GAUSS WTS I=0 DO 4330 IX1=1,NGP1 DO 4330 IX2=1,NGP2 DO 4330 IX3=1,NGM I=I+1 WTFACT=GWT(IX1)*GWT(IST2+IX2)*GWT(IST3+IX3) DO 4330 IL=1,LMAX 4330 PWGHT(I,IL)=PWGHT(I,IL)*WTFACT C 4998 WRITE(6,659) (I,LM(1,I),LM(2,I),LM(3,I),I=1,LMAX) 659 FORMAT('0 BI-SPHERICAL HARMONICS FOR EXPANDING S-MATRIX ARE ', 1'AS FOLLOWS'/'0 INDX L1 L2 LL'/(' ',4I4)) GO TO 4999 C C ITYPE=5,6 -- COMPUTE VLI 4500 I=0 DO 4501 NX=1,NGL DO 4501 IX=1,NGM I=I+1 IL=0 DO 4501 L=IZ,LMBDMX MTOP=MIN0(MUMX,L) DO 4501 M=IZ,MTOP,ICNSYM IF (IHOMO.EQ.2 .AND. ODD(L,M)) GO TO 4501 IL=IL+1 VLI(I,IL)=PLM(L,M,COSA(NX))/SQRT(2.D0*PI) IF (M.NE.0) VLI(I,IL)=VLI(I,IL)*2.D0*COS(DBLE(M)*COSA(NGL+IX)) 4501 CONTINUE C SETUP PWGHT FACTL=1.D0/SQRT(2.D0*PI) IL=0 DO 4502 L=IZ,LMMAX MTOP=MIN0(L,MUMAX) DO 4502 M=IZ,MTOP,ICNSYM IF (IHOMO.EQ.2 .AND.ODD(L,M)) GO TO 4502 IL=IL+1 IV=0 DO 4503 IX=1,NGL DO 4503 NX=1,NGM IV=IV+1 4503 PWGHT(IV,IL)=GWT(IX)*GWT(NGL+NX)*PLM(L,M,COSA(IX))* & COS(DBLE(M)*COSA(NGL+NX))* 2 (SFACT*FACTL) 4502 CONTINUE I=0 IX=0 DO 4505 L=IZ,LMMAX MTOP=MIN0(MUMAX,L) DO 4505 M=IZ,MTOP,ICNSYM IF (IHOMO.EQ.2 .AND. ODD(L,M)) GO TO 4505 I=I+1 LM(1,I)=L LM(2,I)=M DO 4504 IL=1,I IF (LM(1,IL).NE.L) GO TO 4504 IX=IX+1 IXQL(1,IX)=I IXQL(2,IX)=IL IF (I.NE.IL) GO TO 4506 IXQL(3,IX)=0 GO TO 4504 4506 IXQL(3,IX)=1 IX=IX+1 IXQL(1,IX)=I IXQL(2,IX)=IL IXQL(3,IX)=2 4504 CONTINUE 4505 CONTINUE WRITE(6,657) (I,LM(1,I),LM(2,I),I=1,LMAX) 657 FORMAT('0 SPHERICAL HARMONIC SYMMETRIES FOR EXPANDING S-MATRIX ', 1 'ARE AS FOLLOWS'/'0 INDX L M'/(' ',2I4,I3)) WRITE(6,649) 649 FORMAT('0 BELOW ARE INDICES TO SYMMETRIES IN QLT'/ &'0 IN QLT LM1 L M LM2 L M CODE') DO 4507 I=1,NQL IL=IXQL(1,I) IX=IXQL(2,I) 4507 WRITE(6,648) I,IL,LM(1,IL),LM(2,IL),IX,LM(1,IX),LM(2,IX),IXQL(3,I) 648 FORMAT(' ',I7,I6,2I3,I6,2I3,I6) GO TO 4999 C 4999 RETURN C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ENTRY IOSB2(JTOT,LORB,JJJ,NB,CENT,EINT,CINT,WVEC,VL,IVIX,IP, & NVC,ERED,NPOTL,MXLAM,LAM,VLI,NGPT,MXXXXL) C XJ=JTOT XJ=XJ*(XJ+1.D0) DO 6666 I=1,NVC LORB(I)=JTOT JJJ(I)=I NB(I)=I CENT(I)=XJ EINT(I)=CINT*EV(I) DIF=ERED-EINT(I) WVEC(I)=SQRT(ABS(DIF)) IF (DIF.GE.0.D0) GO TO 6666 WVEC(I)=-WVEC(I) 6666 CONTINUE C GO TO IGOTP,(6100,6200,6300,6500) C C CHECK FOR CONSISTENT IVLFL 6100 IF (IVLFL.NE.0) GO TO 9999 DO 6101 I=1,MXLAM IL=LAM(I)/IHOMO + 1 C IVIX(I)=I 6101 VL(I)=VLI(IP,IL) C SET COSANG, FACTOR FOR MXLAM.LE.0 CASE COSANG(1)=COSA(IP) FACTOR=1.D0 GO TO 6900 C C>>SG -------------- NEW CODE --------------------->> C>>SG CHECK FOR NPOTL=1 (LVRTP) OR NPOTL=MXLAM (EXPANDED) C>>SG LATTER CASE UNCHANGED FROM VERSION 10 CODE 6200 IF (NPOTL.EQ.MXLAM .AND.(.NOT.LVRTP)) GO TO 6250 IF (NPOTL.EQ.1.AND.LVRTP) GO TO 6201 WRITE(6,670) NPOTL,MXLAM 670 FORMAT('0 IOSB2 (FEB 92) -- ERROR. LVRTP INCONSISTENT WITH', 1 ' NPOTL, MXLAM',2I6) STOP C C CHECK FOR CONSISTENT IVLFL 6201 IF (IVLFL.LE.0) GO TO 9999 ITOP=NVC*(NVC+1)*NPOTL/2 DO 6202 IX=1,ITOP IVIX(IX)=0 6202 VL(IX)=0.D0 DO 6203 L=1,MXLAM LLL=LAM(3*L-2) IL=LLL/IHOMO+1 C N.B. WE SHOULD HAVE LLL=0 AND IL=1 *** DEBUGGING ONLY *** IF (LLL.NE.0 .OR. IL.NE.1) WRITE(6,672) LLL,IL 672 FORMAT('0 IOSB2 (FEB 92) -- ERROR. LLL.NE.0 .OR IL.NE.1',2I6) IV=LAM(3*L-1) IVP=LAM(3*L) IVVP=0 DO 6204 IROW=1,NVC NV=LEVV(IROW) DO 6204 ICOL=1,IROW NVP=LEVV(ICOL) IVVP=IVVP+1 IF (.NOT.((NV.EQ.IV.AND.NVP.EQ.IVP).OR.(NV.EQ.IVP.AND.NVP.EQ.IV))) 1 GO TO 6204 C IF WE REACH BELOW, THIS ROW/COL CORRESPONDS TO CURRENT 'SYMMETRY' IX=(IVVP-1)*NPOTL+LLL+1 C SINCE NPOTL=1 AND LLL=0, SHOULD HAVE IX=IVVP *** DEBUGGING *** IF (IX.NE.IVVP) WRITE(6,673) IX,IVVP 673 FORMAT('0 IOSB2 (FEB 92) -- ERROR. IX.NE.IVVP FOR VL,IVIX',2I6) IVIX(IX)=L VL(IX)=VLI(IP,IL) 6204 CONTINUE 6203 CONTINUE C SET COSANG, FACTOR FOR VRTP CASE, AND RETURN COSANG(1)=COSA(IP) FACTOR=1.D0 GO TO 6900 C C CODE BELOW IS NPOTL=MXLAM (POTENTIAL EXPANDED IN LEGENDRE POLY'S) C THIS IS ESSENTIALLY CODE FROM MOLSCAT VERSION 9. C C CHECK FOR CONSISTENT IVLFL 6250 IF (IVLFL.NE.0) GO TO 9999 DO 6251 L=1,MXLAM LIX=L IL=LAM(3*L-2)/IHOMO + 1 LV1=LAM(3*L-1) LV2=LAM(3*L) DO 6252 IV=1,NVC DO 6252 IVP=1,IV C IVIX(LIX)=L VL(LIX)=0.D0 IF (LEVV(IV).EQ.LV1 .AND. LEVV(IVP).EQ.LV2) GO TO 6253 IF (LEVV(IV).EQ.LV2 .AND. LEVV(IVP).EQ.LV1) GO TO 6253 GO TO 6252 6253 VL(LIX)=VLI(IP,IL) 6252 LIX=LIX+MXLAM 6251 CONTINUE GO TO 6900 C C CHECK FOR CONSISTENT IVLFL 6300 IF (IVLFL.NE.0) GO TO 9999 IL=0 DO 6301 L1=IZ,L1MAX,IHOMO L2TOP=L2MAX IF (IDENT.GT.0) L2TOP=L1MAX DO 6301 L2=IZ,L2TOP,ICNSYM LLO=ABS(L1-L2) LHI=L1+L2 DO 6301 LL=LLO,LHI IF (ODD(L1+L2,LL)) GO TO 6301 IL=IL+1 DO 6302 I=1,MXLAM IF (L1.NE.LAM(3*I-2)) GO TO 6302 IF (L2.NE.LAM(3*I-1)) GO TO 6302 IF (LL.NE.LAM(3*I )) GO TO 6302 C IVIX(I)=I VL(I)=VLI(IP,IL) IF (LDEBUG) WRITE(6,635) I,L1,L2,LL,IL 635 FORMAT(' IOSB2. DEBUG. I,L1,L2,LL,IL',5I5) 6302 CONTINUE 6301 CONTINUE C FOR 'VRTP' CASE NEED TO SET COSANG(), FACTOR =(4*PI)**(3/2) FACTOR=(4.D0*PI)*SQRT(4.D0*PI) C CALCULATE IX1,IX2,IX3 FROM IP (# OF GAUSS POINT) IX3=IP IX1=(IX3-1)/(NGP2*NGM)+1 IX3=IX3-(IX1-1)*(NGP2*NGM) IX2=(IX3-1)/NGM+1 IX3=IX3-(IX2-1)*NGM COSANG(1)=COSA(IX1) COSANG(2)=COSA(IST2+IX2) COSANG(3)=COSA(IST3+IX3) GO TO 6900 C<>SG N.B. ITYPE=3 VALUES DIFFER FROM GOLFLAM-KOURI AND AGG-CLARY C ITYPE V AVGFCT C 1,2 1. 1. C 3 1/4*PI 1./SQRT(4*PI) C 5,6 1/4*PI 1./SQRT(4*PI) V=1.D0 IF (ITYPE.EQ.5.OR.ITYPE.EQ.6.OR.ITYPE.EQ.3) V=1.D0/(4.D0*PI) AVGFCT=SQRT(V) CINT=RMLMDA/EPSIL C INITIALIZE RSTART, IN CASE IRMSET.LE.0 AND FINDRM NOT CALLED RMINSV=RMIN RSTART=RMIN CALL GCLOCK (TITIME) C C PRINT LEVEL FOR SCATTERING CAN BE LESS THAN FOR IOS1 C IOSPR=MAX0(0,PRINT-10) C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C LOOP OVER ENERGIES. C DO 2000 IE=1,NNRG ICODE=1 IF(IE.GT.1 .AND. ISCRU.GT.0) ICODE=2 IF(ISCRU.GT.0) REWIND ISCRU WRITE(6,622) IE,ENERGY(IE) 622 FORMAT('1 IOSCLC (MAY 92). ENERGY(',I3,') =',F12.4,' (1/CM).') ERED=ENERGY(IE)*CINT IF (IE.EQ.1) EFIRST=ERED ESHIFT=ERED-EFIRST C A MORE SOPHISTICATED WAY OF SAVING RTURN IS PROBABLY WANTED, C BUT BELOW SHOULD WORK AS A TEMPORARY MEASURE RTURN=RMINSV C ZERO STORAGE DO 2100 I=1,NQL 2100 IEC(I)=0 DO 2109 IV=1,NVC DO 2109 IVP=1,NVC QLS(IV,IVP)=0.D0 SIGAV(IV,IVP)=0.D0 DO 2101 I=1,NQL 2101 QLT(IV,IVP,I)=0.D0 DO 2102 I=1,NGPT 2102 SIGTH(IV,IVP,I)=0.D0 2109 CONTINUE C C LOOP OVER PARTIAL WAVES C DO 3000 JTOT=JTOTL,JTOTU,JSTEP IF (PRINT.GT.1) WRITE(6,626) JTOT,IE,ENERGY(IE) 626 FORMAT('0 ***** PARTIAL WAVE =',I5,' FOR ENERGY(',I3,') = ', & F12.4,' *****') C IF (JTOTU.LT.999) GO TO 3001 C CHECK FOR CONVERGENCE C ONLY CHECK FOR QLT WHERE IXQL(NIXQL,IL).EQ.0 DO 3002 IL=1,LMAX IF (IXQL(NIXQL,IL).NE.0) GO TO 3002 IF (IEC(IL).LT.NCAC) GO TO 3001 3002 CONTINUE CALL GCLOCK(TJTIME) TIME=TJTIME-TITIME TITIME=TJTIME JTO=JTOT-JSTEP WRITE(6,620) IE,ENERGY(IE),LMAX,NCAC,TEST,JTOTL,JSTEP,JTO 620 FORMAT('1 ***** ***** ***** CALCULATION AT ENERGY(',I3,') =', 1 F10.2,' (1/CM) ', & ' TERMINATED DUE TO CONVERGENCE FOR',I4,' Q(L).'/ & 22X,'NCAC, TEST =',I4,2E12.4/ 2 22X,'PARTIAL WAVES',I4,' (',I4,' ) ',I5) WRITE(6,641) TIME 641 FORMAT('0 ***** ***** ***** TIME WAS',F9.2,' SEC.'/' ') GO TO 3009 3001 FACTL=(2*JTOT+1)*PI C C GET ANGLE-DEPENDENT SCATTERING / LOOP OVER GAUSS POINTS. DO 3100 IP=1,NGPT C C INITIALIZE SCAT VARIABLES VL, LORB, EINT, ETC. CALL IOSB2(JTOT,X(IXLORB),X(IXJJJ),NB,X(IXCENT),X(IXEINT),CINT, 1 WVEC,X(IXVL),X(IXIV),IP,NVC,ERED,NPOTL,MXLAM,LAMBDA,VLI, 2 NGPT,MXXXXL) C CONV=0.D0 RTURN=RSTART IF(ICODE.NE.1) GOTO 3005 CALL FINDRX(ENERGY(IE),X(IXEINT),X(IXCENT),1,NVC,CINT,RMAX,RSTOP, 1 NOPMAX,IRXSET,IOSPR) C IF(IRMSET.LE.0) GOTO 3005 C GET TEMPORARY STORAGE FOR FINDRM IT1=ICX IXNEXT=IT1+MXLAM CALL CHKSTR(NUSED) CALL FINDRM(X(IXSR),NVC,RSTART,RTURN,IK,X(IT1),X(IXVL),X(IXIV), 1 ERED,X(IXEINT),X(IXCENT),RMLMDA,X(IXSI),MXLAM,NPOTL, 2 XEPS,ITYPE,IOSPR) C RELEASE TEMPORARY STORAGE IXNEXT=IT1 IF(RVFAC.EQ.0.D0) GOTO 3005 RMID=RVFAC*RTURN IF(IOSPR.GE.3) WRITE(6,3003) RMID,RVFAC 3003 FORMAT('0 RMID =',F7.2,' OBTAINED FROM RVFAC =',F6.3) C C NOW READY TO SOLVE 'COUPLED' EQUATIONS; DONE AS CALL TO STORAG. C 3005 NV=NPOTL*NVC*(NVC+1)/2 CALL STORAG( INTFLG,NVC,MXLAM,NV,NPOTL, 1 IXJJJ,IXSR,IXSI,IXKMAT,IXVL,IXIV,IXEINT,IXCENT, 2 IXWV,IXLORB,IXNB, 3 ESHIFT,NOPMAX,DEEP,IK,ICODE,IOSPR, NUMDER) C C INITIALIZE TO UNIT S-MATRIX TO CLEAR 'NON-CLASSICAL' CHANNELS. C 4000 DO 4005 IV=1,NVC DO 4005 IVC=1,NVC DELVVP=0.D0 IF (IV.EQ.IVC) DELVVP=1.D0 SLR(IV,IVC,IP)=DELVVP 4005 SLI(IV,IVC,IP)=0.D0 IF (NOPEN.GT.0) GO TO 4009 WRITE(6,699) IP,NOPEN 699 FORMAT(' * * * NOTE. FOR ORIENTATION',I6,' NOPEN =',I3) GO TO 3100 4009 IF (NOPEN.LE.NVC) GO TO 4008 WRITE(6,698) IP,NOPEN,NVC 698 FORMAT(' * * * ERROR. FOR ORIENTATION',I6,' NOPEN.GT.NVC',2I6) GO TO 3100 4008 IF (CONV.GE.0.D0) GO TO 4007 WRITE(6,696) JTOT,IP 696 FORMAT('0 * * * WARNING. SLR,SLI,SIGTH NOT SET DUE TO LACK OF CON &VERGENCE FOR PART. WAVE',I4,' ORIENTATION',I5) GO TO 3100 C 4007 IF (PRINT.GE.15) WRITE(6,601) 601 FORMAT(' ') NNP=0 DO 4200 N=1,NOPEN IV=NB(N) WV=RM/WVEC(IV) C SET WVEC(IV) TO WAVENUMBER IN 1/ANGSTROMS FOR ISAVEU OUTPUT WVEC(IV)=1.D0/WV WV=WV*WV*FACTL DO 4200 NP=1,NOPEN IVP=NB(NP) NNP=NNP+1 DELVVP=0.D0 IF (IV.EQ.IVP) DELVVP=1.D0 C BELOW CHANGED APR 86 SINCE ONLY INDICES FOR SREAL,SIMAG ARE HERE SLR(IV,IVP,IP)=X(IXSR-1+NNP) SLI(IV,IVP,IP)=X(IXSI-1+NNP) C ACCUMULATE ANGLE-DEPENDENT TOTAL CROSS SECTION. ADD=DELVVP-SLR(IV,IVP,IP) ADD=(ADD*ADD+SLI(IV,IVP,IP)*SLI(IV,IVP,IP) )*WV SIGTH(IV,IVP,IP)=SIGTH(IV,IVP,IP)+ADD IF (PRINT.LT.15) GO TO 4200 WRITE(6,627) IP,IV,IVP,SLR(IV,IVP,IP),SLI(IV,IVP,IP), & ADD,SIGTH(IV,IVP,IP) 627 FORMAT(' FOR ORIENTATION',I6,' VIB LEVEL =',I2,' TO',I2, & ', SREAL, SIMAG =',2D14.6,' SIGTH ADD',D12.4,' = ',D12.4) 4200 CONTINUE 3100 CONTINUE C END OF LOOP OVER ORIENTATIONS C C INTEGRATE OVER ORIENTATIONS TO GET SLLR/SLLI C ** N.B. THESE ARE T-MATRIX COMPONENTS ** IF (PRINT.GE.20) WRITE(6,601) DO 3218 IV=1,NVC DO 3218 IVP=1,NVC DELVVP=0.D0 IF (IV.EQ.IVP) DELVVP=1.D0 DO 3218 L=1,LMAX SLLI(IV,IVP,L)=0.D0 SLLR(IV,IVP,L)=0.D0 DO 3208 NX=1,NGPT SLLR(IV,IVP,L)=SLLR(IV,IVP,L)+(DELVVP-SLR(IV,IVP,NX))*PWGHT(NX,L) 3208 SLLI(IV,IVP,L)=SLLI(IV,IVP,L)-SLI(IV,IVP,NX)*PWGHT(NX,L) IF (PRINT.GE.20) & WRITE(6,648) IV,IVP,L,SLLR(IV,IVP,L),SLLI(IV,IVP,L) 648 FORMAT(5X,3I5, 2D16.8) 3218 CONTINUE C C *** C>>SG MAY 92. CODE BELOW REPLACED BY CALL ISUTP AT STATEMENT NO. 3000 C SAVE SLLR/SLLI HRE / N.B. SLR/SLI MIGHT BE USEFUL LATER. C IF (ISU.LE.0) GO TO 3230 C WRITE(ISU,3231) JTOT,IE,ENERGY(IE) C3231 FORMAT(2I4,E16.8) C WRITE(ISU,3232) NOPEN,(NB(I),JTOT,WVEC(NB(I)),I=1,NOPEN) C3232 FORMAT(I4/(2I4,E16.8)) C WRITE(ISU,3233) (((SLLR(NB(IV),NB(IVP),L),IV=1,NOPEN),IVP=1,NOPEN) C & ,L=1,LMAX) C WRITE(ISU,3233) (((SLLI(NB(IV),NB(IVP),L),IV=1,NOPEN),IVP=1,NOPEN) C & ,L=1,LMAX) C3233 FORMAT(5E16.8) C *** C C COMPUTE QLS (QLOLD PREVIOUSLY) FOR 1ST (TOTALLY SYMMETRIC) CASE 3230 IF (PRINT.GE.10) WRITE(6,601) DO 3220 IV=1,NVC C SET WVEC(IV) TO (2*L+1)*PI/K**2 FOR USE IN GETTING QL'S C>>SG TRAP CLOSED CHANNELS (NEGATIVE WVEC) TO PREVENT ROUND-OFF PROBLEMS IF (WVEC(IV).LE.0.) GO TO 3220 WVEC(IV)=FACTL/(WVEC(IV)*WVEC(IV)) DO 3219 IVP=1,NVC DELVVP=0.D0 IF (IV.EQ.IVP) DELVVP=1.D0 SUMR=0.D0 SUMI=0.D0 DO 3209 NX=1,NGPT SUMR=SUMR+PWGHT(NX,1)*SLR(IV,IVP,NX) SUMI=SUMI+PWGHT(NX,1)*SLI(IV,IVP,NX) 3209 CONTINUE C>>SG BELOW SUFFERS FROM ROUND-OFF ERROR FOR IV=IVP CLOSED C>>SG TEST CASES GIVE V*(SUMR**2+SUMI**2)-DELVVP ABOUT 2.D-13 C>>SG BEST WAY TO FIX THIS IS PROBABLY TO TRAP *CLOSED* CHANNELS SUM2=(V*(SUMR*SUMR+SUMI*SUMI)-DELVVP)*WVEC(IV) QLS (IV,IVP)=QLS (IV,IVP)+SUM2 IF (PRINT.GE.10) WRITE(6,638) IV,IVP,SUM2, QLS(IV,IVP) 638 FORMAT(' FOR QLS( 0) VIB LEV =',I3,' TO',I3,15X, & 'ADD',D12.4,' =',D12.4) 3219 CONTINUE 3220 CONTINUE C C *** IN ACCUMULATING QL DIVERGENT CODE FOR ITYPE=1,2 AND ITYPE=5,6 C IF (ITYPE.EQ.1 .OR. ITYPE.EQ.2) GO TO 8881 IF (ITYPE.EQ.3) GO TO 8883 IF (ITYPE.EQ.5 .OR. ITYPE.EQ.6) GO TO 8885 STOP C C ACCUMULATE QL'S / TEST FOR CONVERGENCE 8881 BIGL=-1.D0 DO 3200 L=1,NQL LMP=L-1 BIGL=BIGL+2.D0 ITEST=0 IF (PRINT.GE.10 .AND. NVC.GT.1) WRITE(6,601) DO 3210 IV=1,NVC DO 3210 IVP=1,NVC TLLR=SLLR(IV,IVP,L) TLLI=SLLI(IV,IVP,L) TLLSQ=(TLLR*TLLR+TLLI*TLLI)*V*WVEC(IV)/BIGL QLT(IV,IVP,L)=QLT(IV,IVP,L)+TLLSQ XTEST=TEST(1) IF (L.GT.1 .OR. IV.NE.IVP) XTEST=TEST(2) IF (TLLSQ.GT.XTEST) ITEST=1 IF (PRINT.LT.10) GO TO 3210 WRITE(6,628) LMP,IV,IVP,TLLSQ,QLT(IV,IVP,L) 628 FORMAT(' FOR QLT(',I3,') VIB LEV =',I3,' TO',I3, & ' IOS T-MATRIX ADD',D12.4,' =',D12.4) 3210 CONTINUE C>>SG 5/12/92 STATEMENT BELOW SHOULD BE UNNECESSARY IF (JTOTU.LT.999) GO TO 3200 C SUPPRESS CONVERGENCE CHECK FOR LOW PARTIAL WAVES. IF (JTOT.LE.3*JSTEP*NCAC) GO TO 3200 IF (IXQL(NIXQL,L).NE.0) GO TO 3200 IEC(L)=IEC(L)+1 IF (ITEST.GT.0) IEC(L)=0 3200 CONTINUE GO TO 3000 C 8883 DO 8873 IL=1,NQL BIGL=(2*LM(3,IL)+1) C N.B. NVC=1 FOR ITYPE=3 TLLR=SLLR(1,1,IL) TLLI=SLLI(1,1,IL) TLLSQ=(TLLR*TLLR+TLLI*TLLI)*V*WVEC(1) * BIGL QLT(1,1,IL)=QLT(1,1,IL)+TLLSQ IF (PRINT.GE.10) WRITE(6,652) IL,LM(1,IL),LM(2,IL),LM(3,IL), 2 TLLSQ,QLT(1,1,IL) 652 FORMAT(' FOR QLT(',I3,'), L1,L2,L =',3I3,' ADD', 1 D12.4,' =',D12.4) XTEST=TEST(MIN0(2,IL)) IF (JTOT.LE.3*JSTEP*NCAC) GO TO 8873 C IF (IXQL(NIXQL,IL).NE.0) GO TO 8875 --- SHOULD ALL = 0 IEC(IL)=IEC(IL)+1 IF (TLLSQ.GT.XTEST) IEC(IL)=0 8873 CONTINUE GO TO 3000 C 8885 DO 8875 IL=1,NQL C N.B. NVC=1 FOR ITYPE=5,6 TLLR=SLLR(1,1,IXQL(1,IL)) TLLI=SLLI(1,1,IXQL(1,IL)) TLLR1=SLLR(1,1,IXQL(2,IL)) TLLI1=SLLI(1,1,IXQL(2,IL)) IF (IXQL(3,IL).EQ.2) GO TO 8865 C BELOW FOR REAL PART / ALSO FOR DIAGONAL CASES TLLSQ=(TLLR*TLLR1+TLLI*TLLI1)*V*WVEC(1) GO TO 8855 C BELOW FOR IMAGINARY PART 8865 TLLSQ=(TLLI*TLLR1-TLLR*TLLI1)*V*WVEC(1) 8855 QLT(1,1,IL)=QLT(1,1,IL)+TLLSQ IF (PRINT.GE.10) WRITE(6,651) IL,LM(1,IXQL(1,IL)),LM(2,IXQL(1,IL)) 1 ,LM(2,IXQL(2,IL)),IXQL(3,IL), 2 TLLSQ,QLT(1,1,IL) 651 FORMAT(' FOR QLT(',I3,'), L,M,M1 =',3I4,', CODE =',I2,' ADD', 1 D12.4,' =',D12.4) XTEST=TEST(MIN0(2,IL)) IF (JTOT.LE.3*JSTEP*NCAC) GO TO 8875 IF (IXQL(3,IL).NE.0) GO TO 8875 IEC(IL)=IEC(IL)+1 IF (TLLSQ.GT.XTEST) IEC(IL)=0 8875 CONTINUE GO TO 3000 C 3000 CALL ISUTP(ISU,ENERGY(IE),JTOTL,JSTEP,JTOT,NVC,NQL,QLS,QLT) C END OF LOOP OVER PARTIAL WAVES C CALL GCLOCK(TJTIME) TIME=TJTIME-TITIME TITIME=TJTIME WRITE(6,631) ENERGY(IE),JTOTL,JSTEP,JTOTU 631 FORMAT('1 ***** ***** ***** END OF CALCULATION FOR ENERGY =', 1 F12.4,' (1/CM) ***** ***** *****'/ & 22X,'PARTIAL WAVES',I4,' (',I4,' ) ',I5) WRITE(6,641) TIME C C END OF CALCULATION FOR THIS ENERGY / OUTPUT CROSS SECTIONS C MAKE SURE WE HAVE NUSED BY CALLING CHKSTR 3009 CALL CHKSTR(NUSED) WRITE(6,684) NUSED,MX 684 FORMAT('0',2(' *****'),' STORAGE SO FAR USED',I10,' OF THE', 1 I10,' AVAILABLE WORDS.') C>>SG C>>SG N.B. NVC SHOULD BE LOWERED TO NOUT=NOPEN (AS IN IOSOUT) C>>SG DO 3305 NX=1,NGPT IV=1 WRITE(6,632) NX,(LEFT,IV,IVP,SIGTH(IV,IVP,NX),IVP=1,NVC) 632 FORMAT('0 FOR ORIENTATION',I6,3(5X,A4 ,I2,',',I2,') =',1PE12.4) & /(23X,3(5X,A4,I2,',',I2,') =',1PE12.4))) IF (NVC.LE.1) GO TO 3008 DO 3007 IV=2,NVC 3007 WRITE(6,642) (LEFT,IV,IVP,SIGTH(IV,IVP,NX),IVP=1,NVC) C>>SG FORMAT CHANGED 2/6/92 TO ELIMINATE APPARENT COMPILER BUG 642 FORMAT(23X,3(5X, A4, I2,',',I2,') =',1PE12.4)/ 1 (23X,3(5X, A4, I2,',',I2,') =',1PE12.4))) 3008 DO 3305 IV=1,NVC DO 3305 IVP=1,NVC 3305 SIGAV(IV,IVP)=SIGAV(IV,IVP)+PWGHT(NX,1)*SIGTH(IV,IVP,NX)*AVGFCT WRITE(6,643) 643 FORMAT('0 AVERAGE OVER ORIENTATIONS') DO 3004 IV=1,NVC 3004 WRITE(6,642) (LEFT,IV,IVP,SIGAV(IV,IVP),IVP=1,NVC) C C CALL IOSOUT/IOSPB TO GET STATE TO STATE AND PR. BR. CROSS SECTIONS C N.B. ATAU, NEEDED ONLY FOR SIG6, IS STORED IN X(1), I.E., JLEV. C THIS IS PRETTY BAD CODING; BETTER PASSING OF ATAU TO SIG6 NEEDED C IATAU=1 CALL IOSOUT(ENERGY(IE),QLT,QLS,NVC,ITYPE,X(IATAU),LM,IXQL, 1 LMAX,NIXQL,NQL,JSTEP) IF(IFLS.GT.0) 1CALL IOSPB(ENERGY(IE),QLT,QLS,IFLS,LINE,LTYPE,ITYPE,NVC,LM,IXQL, 1 LMAX,NIXQL,NQL) C 2000 CONTINUE C C END OF LOOP OVER ENERGIES. C RETURN END SUBROUTINE IOSDRV(NNRG,NPR,ENERGY,JTOTL,JTOTU,JSTEP,TEST,NCAC, 1 IFLS,LINE,LTYPE,MXLN,INTFLG,ITYPE,LMAX,MMAX, 2 IPROGM,URED,LABEL,NUMDER, 3 LAMBDA,MXLAM,NPOTL,CINT,IRMSET,IRXSET,RVFAC, 4 DEEP,PRINT,NVC, ISAVEU,TITIME,RM,EPSIL,RMIN,RMAX) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C INTEGRATED MOLSCAT/IOS IMPLEMENTED APR 86 CAMBRIDGE, ENGLAND. C -- A GUTTED VERSION OF IOS1, INTERFACED TO CCP6 MOLSCAT. C THIS IS A DRIVER FOR THE IOS CODE; MOLSCAT/DRIVER CALLS C BASIN TO READ &BASIS, WHICH THEN CALLS IOSBIN. C DRIVER THEN CALLS POTENL TO GET &POTL DATA, C AND FINALLY CALLS IOSDRV TO SET UP AND PERFORM IOS CALCULATION. C INTEGER PRINT DIMENSION ENERGY(NNRG),TEST(2),LINE(2,MXLN),LTYPE(MXLN), 1 LAMBDA(MXLAM) LOGICAL NUMDER CHARACTER*80 LABEL C C LAST CHANGED 1/19/93. NEW DYNAMIC MEMORY HANDLING C ** VERSION 6 / OCT 85/ ADDS ITYPE=6 CAPABILITY C / ALSO ALLOWS "UNEXPANDED" POTL, V(R,ANGLES) C ** VERSION 5 / MAR 81/ ADDS INTFLG=4 (MOLSCAT V.8) C / JUNE 82/ REPLACES PLM WITH R. T PACK VERSION. C ** VERSION 4 / MAY. 78/ ADDS ITYPE=5 CODE. C / SEP. 78/ **TEMPORARY** ISAVEU CPABILITY C / APR. 79/ CHANGED FOR ISCRU (MOLSCAT V.7) COMPATABIL C ** VERSION 3 / DEC. 77/ IS TOTALLY NEW ORGANIZATION TO ACCOMMODATE C ITYPE=2 (VIBROTOR - ATOM) C ** VERSION 2 / OCT. 77/ ADDS WKB (R.T PACK) CAPABILITY ** C ** VERSION 1 / SEP. 77/ INTERFACE HOUSTON PROGRAM W/MOLSCAT. 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 POINTED W/ VL ARRAY. C NIPR IS NUMBER OF INTEGERS PER REAL; SHOULD BE 1 OR 2. C C CMBASE MODIFIED TO MATCH CURRENT SPECS IN MOLSCAT/BASIS C THIS IS USED ONLY TO GET LV,EV VALUES IN IOSBIN SO THEY CAN BE C WRITTEN TO ISAVEU HERE. DIMENSION EV(200),LV(400) COMMON /CMBASE/ DUM(214),IDUM(408) EQUIVALENCE (EV(1),DUM(11)),(LV(1),IDUM(8)) C C MUST INITIALIZE NUSED NON-NEGATIVE BEFORE CALL CHKSTR NUSED=0 WRITE(6,68) 68 FORMAT('0 IOSDRV ENTERED. SET-UP FOR INFINITE ORDER SUDDEN', 1 ' CALCULATION.') C C CONTINUE WITH SET-UP FOR IOS. PROCESS &POTL LAM(MXLAM) DATA C SET NGPT, LMAX AND GAUSS PTS/WTS. C N.B. LMAX/MMAX INITIALLY CONTAIN HIGHEST L,M VALUES C DESIRED FOR QLM. LMAX IS RESET TO EQUAL THE *NUMBER* OF L,M C VALUES IN LM,SLLR,SLLI,ETC. CALL IOSBGP(MXLAM,LAMBDA,MXXXXL,NGPT,LMAX,MMAX,NQL,NIXQL) C C RESERVE STORAGE FOR VARIABLES. IC HAS NEXT AVAIL LOC IN X() C STORAGE FOR SCATTERING VARIABLES . . . C SREAL(NVC,NVC),SIMAG(NVC,NVC),WVEC(NVC),EINT(NVC),CENT(NVC), C VL(NVC*(NVC+1)/2,MXLAM),JJJ(NVC),LORB(NVC),NB(NVC) C --ADDED APR 86-- KMAT(NVC,NVC),IV(NVC*(NVC+1)/2,MXLAM) C C V11 CODE EXPECTED IC TO BE STORAGE USED SO FAR ISVMEM=IXNEXT IC=IXNEXT-1 IXSR=IC+1 IXSR=IXNEXT IXSI=IXSR+NVC*NVC IXKMAT=IXSI+NVC*NVC IXWV=IXKMAT+NVC*NVC IXEINT=IXWV+NVC IXCENT=IXEINT+NVC IXVL=IXCENT+NVC NV=NVC*(NVC+1)*NPOTL/2 IXJJJ=IXVL+NV IXLORB=IXJJJ+NVC IXNB=IXLORB+NVC IXIV=IXNB+NVC IC=IXIV IF (IVLFL.GT.0) IC=IXIV+(NV+NIPR-1