C  calculate BRW value of average energy transfer
C  AUTHOR: GILBERT.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TITLE(20)
      DOUBLE PRECISION M,NU,K,MLIGHT,MAV,MBG,MM,KTWVN
C THE FOLLOWING DATA ARE LOCAL EPSILON AND SIGMA VALUES
C  FOR INDIVIDUAL ATOMS
      DATA SIGH/3./,SIGC/3.2/,SIGO/3.2/,SIGF/3.2/,SIGCL/3.4/,
     1 SIGS/3.4/,SIGBR/3.6/,SIGI/4.0/,SIGN/3.2/,
     2 EPSH/6.5/,EPSC/20.3/,EPSO/20.3/,EPSF/20.3/,EPSCL/120./,
     3 EPSS/120./,EPSBR/190./,EPSI/230./,EPSN/20.3/
C     OPEN(UNIT=5,STATUS='OLD',FILE='BRW.DAT')
C     OPEN(UNIT=6,STATUS='NEW',FILE='BRW.OUT')
      WRITE(6,1000)
 1000 FORMAT(10X,20(2H *),/,
     1 ' BIASED RANDOM WALK MODEL B. VERSION: June 29 1991')
      READ(5,789) TITLE
 789  FORMAT(20A4)
      WRITE(6,789) TITLE
       DELR = 0.0005
      DELR = DELR*1E-10
 300  CONTINUE
      READ(5,*) EPS,SIGMA,M,T,NATOMS,MLIGHT
      IF(EPS.LT.0.) GO TO 987
      READ(5,*) EPSBG,SIGBG,MBG
      WRITE(6,1) EPSBG,SIGBG,MBG
   1  FORMAT(///,' BATH GAS: EPSILON (K) =',F10.2,/,
     1       '           SIGMA (A)   =',F10.2,/,
     2       '           MASS (AMU)  =',F10.2)
C FIND KT IN CM-1
      KTWVN = T/1.439
      WRITE(6,2) NATOMS
   2  FORMAT(' NUMBER OF DISTINCT ATOMS:',I5)
      MAV=0.
      EPSRL = 0.
      SIGRL=0.
      NTOT = 0
      WRITE(6,700)
  700 FORMAT(/,'  INPUT ATOMS:',/,
     1 4X,' MASS    STOICH.CFT.    LOCAL EPS.   LOCAL SIGMA')
      DO 4 I=1,NATOMS
      READ(5,*) MM,NST,EPSL,SIGL
      IF(EPSL.GT.0.) GO TO 704
C ASSIGN DEFAULT VALUES TO LOCAL EPSILON, SIGMA
      IF(MM.GT.3.5) GO TO 500
C  HYDROGEN OR AN ISOTOPE
      EPSL = EPSH
      SIGL = SIGH
      GO TO 704
 500  IF(MM.GT.13.5) GO TO 501
C  CARBON
      EPSL = EPSC
      SIGL = SIGC
      GO TO 704
 501  IF(MM.GT.15.) GO TO 502
C  NITROGEN
      EPSL = EPSN
      SIGL = SIGN
      GO TO 704
 502  IF(MM.GT.18.5) GO TO 503
C  OXYGEN
      EPSL = EPSO
      SIGL = SIGO
      GO TO 704
 503  IF(MM.GT.20.) GO TO 504
C  FLUORINE
      EPSL = EPSF
      SIGL = SIGF
      GO TO 704
 504  IF(MM.GT.34.) GO TO 505
C SULFUR
      EPSL = EPSS
      SIGL = SIGS
      GO TO 704
 505  IF(MM.GT.38.) GO TO 506
C  CHLORINE
      EPSL = EPSCL
      SIGL = SIGCL
      GO TO 704
 506  IF(MM.GT.81.) GO TO 507
C  BROMINE
      EPSL = EPSBR
      SIGL = SIGBR
      GO TO 704
C ONLY IODINE IS LEFT
 507  EPSL = EPSI
      SIGL = SIGI
 704  CONTINUE
      WRITE(6,703) MM,NST,EPSL,SIGL
  703 FORMAT(F10.2,I10,5X,F10.2,5X,F10.2)
      FNST = NST
      MAV = MAV + FNST*MM
      NTOT = NTOT+NST
      EPSRL = EPSRL + FNST*EPSL
      SIGRL = SIGRL + FNST*SIGL
    4 CONTINUE
      TOT = NTOT
      MAV = MAV/TOT
C MODEL A      MAV = MAV*MBG/(MAV+MBG)
C  MODEL B:
      MAV = 1./((1./(MAV*TOT-MAV)) + (1./MAV))
      ELOC = EPSRL/TOT
      SLOC = SIGRL/TOT
      ELOC = DSQRT(ELOC*EPSBG)
      SLOC = 0.5*(SLOC+SIGBG)      
      READ(5,*) ECOLL,NVIB,NU
      TSTAR = T/EPS
      OM22 = 1./(0.636 + 0.567*DLOG10(TSTAR))
      WRITE(6,3) EPS,SIGMA,M,T,OM22,MAV,MLIGHT
      ERAT = ECOLL*1.439/(T*NVIB)
    3   FORMAT(/,' EPSILON (K) =',T23,F10.0,/,' SIGMA (A) =',T23,
     1 F10.2,/,' RED. MASS (A.M.U.) =',T23,F10.2,/,
     2 ' TEMPERATURE (K) =',T23,F10.0,/,
     3 ' OMEGA 22 =',T23,F10.3,
     4 /,' AVERAGE ATOM/ATOM MASS =',F10.2,' AMU',
     5 /,' MASS LIGHTEST ATOM =',F10.2,' AMU')
      WRITE(6,101) ECOLL,NVIB,NU
 101  FORMAT(' COLLISION ENERGY (CM-1):',F10.0,/,
     1 ' (NOT USED IN MODEL B)',//,
     2 ' NO. VIBRATIONS:',I5, '   (NOT USED IN MODEL B)',
     3  /,' HIGHEST FREQUENCY (CM-1):', F10.0)
      NU = NU*3E10
      K = 39.5*(MLIGHT*1.66E-27)*NU*NU
      WRITE(6,31) SLOC,ELOC
  31  FORMAT(' LOCAL SIGMA, EPSILON',2F8.2)
      SLOC=SLOC*1E-10
      ELOC=ELOC*1.38E-23
C MODEL A      EBAR = ECOLL*1.986E-23/NVIB
C THIS PROGRAM IS FOR MODEL B
C  MODEL B:
      EBAR = 2.*KTWVN*1.986E-23
      EPS = EPS*1.38E-23
      SIGMA = SIGMA*1E-10
      DELTAX = DSQRT(2.*EBAR/K)
      C = 2.*3.14159*NU
C     WRITE(6,78) C
C 78  FORMAT(' C =',1P,E10.2,' S-1')
      M = M*1.66E-27
      E = 2.*1.38E-23*T
      D= SIGMA*OM22
      D = DMAX1(D,SIGMA)
      BOND=0.666667
      B = BOND*D
C  B IS AVERAGE IMPACT PARAMETER 
C     WRITE(6,51) B*1E10
C 51  FORMAT(//,' IMPACT PARAMETER (A):',F5.3)
      N=0
      FAVSLW=0.
      TAUC = 0.
      R = D
   5  SONR = SIGMA/R
      SONR6 = SONR**6
      SONR12 = SONR6*SONR6
      VEFF = 4.*EPS*(SONR12-SONR6) + E*((B/R)**2)
      IF(VEFF.GE.E) GO TO 10
      TAUC = TAUC +(1./SQRT(E-VEFF))
      R = R-DELR
      GO TO 5
  10  TAUC = TAUC*DELR*SQRT(2.*M)
C     WRITE(6,52) R*1E10
C 52  FORMAT(' TURNING PT (A) =',F6.3)
      WRITE(6,9) TAUC
  9   FORMAT(/,' COLLISION TIME (S) =',1P,E9.2)
      D= SLOC*OM22
      D = DMAX1(D,SLOC)
      BOND=0.666667
      B = BOND*D
      N=0
      FAVFST=0.
      E=EBAR
      R = D
   6  SONR = SLOC/R
      SONR6 = SONR**6
      SONR12 = SONR6*SONR6
      VEFF = 4.*ELOC*(SONR12-SONR6) + E*((B/R)**2)
      IF(VEFF.GE.E) GO TO 11
      FORCE = 4.*ELOC*(-12.*SONR12+6.*SONR6) -2.*E*((B/R)**2)
      FORCE = FORCE/R
      FAVFST=FAVFST+DABS(FORCE)
      N=N+1
      R = R-DELR
      GO TO 6
 11    CONTINUE
C      WRITE(6,52) R*1E10
      FORCE = 4.*ELOC*(-12.*SONR12+6.*SONR6) -2.*E*((B/R)**2)
      FORCE = FORCE/R
C      WRITE(6,21) FORCE
C  21   FORMAT(/,' FAST FORCE AT TURNING POINT (NEWTON) =',1P,E9.2)
      FFAST=DABS(FORCE)
      FAVFST = FAVFST/N
      F = 1.20
C      WRITE(6,1002) F
C1002 FORMAT(' SIGMA/R FOR :',F10.3)
      FTERM = DABS((4.)*(-12.*(F**12) + 6.*(F**6)))
      X = SLOC/F
      DELTAV = FTERM*ELOC*DELTAX/X
      IF(DELTAV.GT.0.5*EBAR) DELTAV = 0.5*EBAR
      EIDOT = (EBAR-DELTAV)*NU
      ECM = EIDOT/1.986E-23
      ECM2=ECM*ECM
      DVCM = DELTAV/1.986E-23
C     WRITE(6,86) DVCM
C86   FORMAT(' DELTA V (CM-1):',F10.1)
C      WRITE(6,120) ECM2
C120  FORMAT('  WITH MODIFIED DELTA V',1P,E10.2,
C    1 ' CM-2 S-2')
      EIDOT2=EIDOT*EIDOT
C     WRITE(6,71) FAVFST
C 71  FORMAT(' AVERAGE FAST FORCE =',1P,E9.2)
      FAV = FFAST
      FAV = DABS(FAV)
C     WRITE(6,127) FAV
C127  FORMAT(' AVERAGE FORCE USED:',1P,E9.2,' NEWTON')
      A = FAV/SQRT(MAV*1.66E-27*EBAR*0.5)
C     WRITE(6,2233) A
C2233 FORMAT(' CALC. A =',1P,E9.2,' S-1')
      S = 2.*A*EIDOT2*TAUC/(A*A+C*C)
      S = DSQRT(S)/1.986E-23
C  THIS VALUE OF S IS THAT WITHOUT THE SEMI-EMPIRICAL CORRECTION
      WRITE(6,246) S
  246  FORMAT(' S = ',
     1 F10.1,' CM-1',/,10x,'**********')
      IF(S.GT.1000.) WRITE(6,387)
 387  FORMAT(' WARNING: A VERY HIGH VALUE, METHOD MAY BE INVALID'
     1 ,/,' FOR THESE PARAMETERS')
      GO TO 300
987    CONTINUE   
C      CLOSE(5,STATUS='KEEP')
C      CLOSE(6,STATUS='KEEP')
      STOP
      END