C PARAMS.FOR PARAMS.FOR C The various parameterizations of physical processes needed C by drifter simulation program. These were separated out of C funcdrft.for because they are subject to tuning whereas those C routines remaining in funcdrft.for are stable. C 01 Mar 2012 C 07 May 2012 Lots of improvements, corrections, etc. C 07 Jan 2013 Add minority positive charge carrier C 03 Aug 2013 Make this compatible with hilumdrift C 11 Jan 2014 Add IV(7)=9 to FUNCTION PROB C 19 Mar 2014 Add IV(5)=8 to FUNCTION RCOMB C 29 Mar 2014 Improve FUNCTION RCMXB C234567890123456789012345678901234567890123456789012345678901234567890 C 1 2 3 4 5 6 7 REAL*8 FUNCTION VP(I) C Velocity of positive ions at edge I PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL C DATA PMOB/0.2/ !Positive ion mobility (mm^2/Vs) (Freeman's No.) C DATA PMOB/3.33/ !Positive ion mobility (much larger than actual) C DATA PMOB/2.5E1/ !Ridiculously large to see subtle effects C VP=PMOB*EDZ(I) C RETURN END REAL*8 FUNCTION VB(I) C Velocity of minority positive charge carriers at edge I PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL C VB=BMOB*EDZ(I) C RETURN END REAL*8 FUNCTION VX(I) C Velocity of negative oxygen ions at edge I PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL C see J.Chem. Phys V37 p2470 (1962) REAL*8 S,E C VX=-XMOB*EDZ(I) C RETURN END REAL*8 FUNCTION VM(I) C Electron velocity at edge I PARAMETER (N=2048) COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL LOGICAL FIRST DATA FIRST/.TRUE./ REAL*8 E,S REAL*8 E00 DATA E00/1000./ !Nominal electric field where VM=VE COMMON /VPARS/EaL,EbL,aeL,beL,ceL,AAL,BAL REAL*8 EaL,EbL,aeL,beL,ceL,AAL,BAL REAL*8 Ea(2),Eb(2),ae(2),be(2),ce(2),AA(2),BA(2) DATA Ea/10.,20./ DATA Eb/1.E4,2000./ DATA ae/1.0,1.0/ DATA be/0.0,0.365/ DATA ce/0.0,0.1/ DATA AA/5.E6,1.2E6/ C IF(FIRST) THEN II=IV(9) IF(II.EQ.1) THEN EMOB=4.2E5 Ea(1)=VE/EMOB AA(1)=VE BA(1)=VE ELSEIF(II.EQ.2) THEN EMOB=5.E4 AA(2)=EMOB*Ea(2) be(2)=DLOG(VE/AA(2))/DLOG(E00/Ea(2)) BA(2)=VE/(E00/Eb(2))**be(2) ELSE PRINT 100 100 FORMAT(' IV(9) is not set properly. Must be 1 or 2') ENDIF EaL=Ea(II) EbL=Eb(II) aeL=ae(II) beL=be(II) ceL=ce(II) AAL=AA(II) BAL=BA(II) C 1 2 3 4 5 6 7 C PRINT 200, EaL,EbL,aeL,beL,ceL,AAL/1.E6,BAL/1.E6 200 FORMAT(5X,'EaL,',8X,'EbL,',8X,'aeL,',8X,'beL,',8X,'ceL,',8X, 1 'AAL,',10X,'BAL = ',/, 2 5(F10.3,2X),2(F10.3,'M',2X)) FIRST=.FALSE. ENDIF C E=ABS(EDZ(I)) S=0. IF(E.NE.0.) S=EDZ(I)/E IF(E.LT.EaL) THEN VM=AAL*(E/EaL)**aeL ELSEIF((E.GE.EaL).AND.(E.LT.EbL)) THEN VM=BAL*(E/EbL)**beL ELSE VM=BAL*(E/EbL)**ceL ENDIF VM=-S*ESCAL*VM C RETURN END SUBROUTINE INITDIFF(I) C Initialize the diffusion constants COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL COMMON /DIFFUS/DFP,DFB,DFM,DFX,DFG REAL*8 DFP,DFB,DFM,DFX,DFG REAL*8 KBT !Thermal energy at 90K C KBT=1.0/40.0*90.0/273.0 DFP=KBT*PMOB DFB=KBT*BMOB DFM=KBT*EMOB DFX=KBT*XMOB DFG=KBT*XMOB IF(I.EQ.0) THEN DFP=0.0D0 DFB=0.0D0 DFM=0.0D0 DFX=0.0D0 DFG=0.0D0 ENDIF C RETURN END REAL*8 FUNCTION VPD(I) C For positive argon ions C calculate the average diffusion velocity at the upper edge of cell I C Fictitious cell 0 has the left electrode surface as the upper edge COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /DIFFUS/DFP,DFB,DFM,DFX,DFG REAL*8 DFP,DFB,DFM,DFX,DFG C VPD=0.0 J=I IF(I.EQ.0) J=1 IF((DFP.GT.1.D-10).AND.(QP(J).GT.1.D-10)) THEN IF((I.GT.0).AND.(I.LT.NN)) THEN VPD=2.0*(QP(I+1)-QP(I))/(QP(I+1)+QP(I)) ELSEIF(I.EQ.0) THEN VPD=+2.0 ELSEIF(I.EQ.NN) THEN VPD=-2.0 ENDIF VPD=-VPD*DFP/DZ ENDIF C RETURN END REAL*8 FUNCTION VBD(I) C For minority positives C calculate the average diffusion velocity at the upper edge of cell I C Fictitious cell 0 has the left electrode surface as the upper edge COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /DIFFUS/DFP,DFB,DFM,DFX,DFG REAL*8 DFP,DFB,DFM,DFX,DFG C VBD=0.0 J=I IF(I.EQ.0) J=1 IF((DFB.GT.1.D-10).AND.(QB(J).GT.1.D-10)) THEN IF((I.GT.0).AND.(I.LT.NN)) THEN VBD=2.0*(QB(I+1)-QB(I))/(QB(I+1)+QB(I)) ELSEIF(I.EQ.0) THEN VBD=+2.0 ELSEIF(I.EQ.NN) THEN VBD=-2.0 ENDIF VBD=-VBD*DFB/DZ ENDIF C RETURN END REAL*8 FUNCTION VMD(I) C For electrons C calculate the average diffusion velocity at the upper edge of cell I C Fictitious cell 0 has the left electrode surface as the upper edge COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /DIFFUS/DFP,DFB,DFM,DFX,DFG REAL*8 DFP,DFB,DFM,DFX,DFG C VMD=0.0 J=I IF(I.EQ.0) J=1 IF((DFM.GT.1.D-10).AND.(QM(J).GT.1.D-10)) THEN IF((I.GT.0).AND.(I.LT.NN)) THEN VMD=2.0*(QM(I+1)-QM(I))/(QM(I+1)+QM(I)) ELSEIF(I.EQ.0) THEN VMD=+2.0 ELSEIF(I.EQ.NN) THEN VMD=-2.0 ENDIF VMD=-VMD*DFM/DZ*ESCAL ! DFM is not already scaled by ESCAL ENDIF C RETURN END REAL*8 FUNCTION VXD(I) C For negative O2 molecular ions C calculate the average diffusion velocity at the upper edge of cell I C Fictitious cell 0 has the left electrode surface as the upper edge COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /DIFFUS/DFP,DFB,DFM,DFX,DFG REAL*8 DFP,DFB,DFM,DFX,DFG C VXD=0.0 J=I IF(I.EQ.0) J=1 IF((DFX.GT.1.D-10).AND.(QX(J).GT.1.D-10)) THEN IF((I.GT.0).AND.(I.LT.NN)) THEN VXD=2.0*(QX(I+1)-QX(I))/(QX(I+1)+QX(I)) ELSEIF(I.EQ.0) THEN VXD=+2.0 ELSEIF(I.EQ.NN) THEN VXD=-2.0 ENDIF VXD=-VXD*DFX/DZ ENDIF C RETURN END REAL*8 FUNCTION VGD(I) C For neutral O2 molecules C calculate the average diffusion velocity at the upper edge of cell I C Fictitious cell 0 has the left electrode surface as the upper edge C Assume neutral O2 molecules don't stick to electrodes COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ COMMON /DIFFUS/DFP,DFB,DFM,DFX,DFG REAL*8 DFP,DFB,DFM,DFX,DFG C VGD=0.0 J=I IF(I.EQ.0) J=1 IF((DFG.GT.1.D-10).AND.(QG(J).GT.1.D-10)) THEN IF((I.GT.0).AND.(I.LT.NN)) THEN VGD=2.0*(QG(I+1)-QG(I))/(QG(I+1)+QG(I)) ELSEIF(I.EQ.0) THEN VGD=+0.0 ELSEIF(I.EQ.NN) THEN VGD=-0.0 ENDIF VGD=-VGD*DFG/DZ ENDIF C RETURN END REAL*8 FUNCTION RCMBB(I) C Bulk recombination rate of electrons and minority positives C as function of E-field in cell I C Temporarily the same as RCOMB COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 E C RCMBB=RCOMB(I)*1.D0 C RETURN END REAL*8 FUNCTION RCMXB(I) C Bulk recombination rate of O2- with minority positive COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL C If RCOMX is scaled below Debye value, then so is RCMXB C RCMXB=RCOMX*(BMOB+XMOB)/(PMOB+XMOB) C RETURN END REAL*8 FUNCTION RCOMB(I) C Bulk recombination rate as function of E-field COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 PI,E0,YMIN,PMAX DATA PI/3.14159265/ REAL*8 SHINK,SHINS,E DATA SHINK/0.1/ !Shinsaka bulk recomb in mm^3/s, equivalent to SCRPTR=104275. C RCOMB=RCOM E=ABS(EZ(I)) IF(IV(5).EQ.0) THEN RCOMB=RCOM ELSEIF(IV(5).EQ.1) THEN IF(E.LT.60.) THEN RCOMB=SHINK*ESCAL/DZ**2 ENDIF ELSEIF(IV(5).EQ.2) THEN IF(E.LT.60.) THEN RCOMB=SHINK*ESCAL/DZ**2 ELSE RCOMB=(SHINK*ESCAL/DZ**2-RCOM)*EXP(-(E-60.)/40.)+RCOM ENDIF ELSEIF(IV(5).EQ.3) THEN IF(E.LT.60.) THEN RCOMB=SHINK*ESCAL/DZ**2*(1.5-E/60.) ELSE RCOMB=(SHINK*ESCAL/2/DZ**2-RCOM)*EXP(-(E-60.)/52.)+RCOM ENDIF ELSEIF(IV(5).EQ.4) THEN SHINS=SHINK/2. IF(E.LT.60.) THEN RCOMB=SHINS*ESCAL/DZ**2*(1.5-E/60.) ELSE RCOMB=(SHINS*ESCAL/2/DZ**2-RCOM)*EXP(-(E-60.)/52.)+RCOM ENDIF ELSEIF(IV(5).EQ.5) THEN IF(E.LT.15.) THEN RCOMB=SHINK*ESCAL/DZ**2*(0.8+1.7*E/15.) ELSE RCOMB=(2.5*SHINK*ESCAL/DZ**2-RCOM)*EXP(-(E-15.)/42.)+RCOM ENDIF ELSEIF(IV(5).EQ.6) THEN SHINS=SHINK/4. IF(E.LT.60.) THEN RCOMB=SHINS*ESCAL/DZ**2*(1.5-E/60.) ELSE RCOMB=(SHINS*ESCAL/2/DZ**2-RCOM)*EXP(-(E-60.)/52.)+RCOM ENDIF ELSEIF(IV(5).EQ.7) THEN SHINS=SHINK/8. IF(E.LT.60.) THEN RCOMB=SHINS*ESCAL/DZ**2*(1.5-E/60.) ELSE RCOMB=(SHINS*ESCAL/2/DZ**2-RCOM)*EXP(-(E-60.)/52.)+RCOM ENDIF ELSEIF(IV(5).EQ.8) THEN SHINS=SHINK/4. IF(E.LT.120.) THEN RCOMB=SHINS*ESCAL/DZ**2*(1.5-E/120.) ELSE RCOMB=(SHINS/2.*ESCAL/DZ**2-RCOM)*EXP(-(E-120.)/104.)+RCOM ENDIF ELSEIF(IV(5).EQ.9) THEN E0=300 ! half-way point between min and max YMIN=0.0017 ! Lowest value PMAX=0.7 ! YMAX is YMIN*10**PMAX RCOMB=PMAX*(0.5-ATAN(7.*DLOG10(E/E0))/PI) RCOMB=YMIN*10**RCOMB*ESCAL/DZ**2 ENDIF C RETURN END REAL*8 FUNCTION RATCH(I) C Bulk attachment rate of electrons to O2 molecules C as function of E-field in cell I C Modifications as of 13 June 2012 COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 E C E=ABS(EZ(I)) C IF(E.LT.11.977) E=11.977 ! Old version before 13 Jun 2012 IF(E.LT.60.) E=60. ! New version on 13 Jun 2012 C First rate constant in units of ppm^(-1) us^(-1) a al Zeitnitz et al., NIMA 545 (2005) 613 RATCH=201.6/E+0.168 C IF(E.LT.11.977) RATCH=17.0 ! Old version before 13 Jun 2012 C Now convert to units of mm^3/s C Also see G.Bakale et al., J.Phys.Chem 80 (1976) 2556. And see W.Schmidt book p52, Eq 14 C to convert from /Mole/s to cm^3/s RATCH=RATCH*1.D6*39.95/6.022E23/1.396*1.D3*1.D6 C Account for number in cell rather than number density RATCH=RATCH*ESCAL/DZ**2 C RETURN END REAL*8 FUNCTION RCOMX(I) C Bulk recombination rate of O2- with Ar+ COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 PARAMETER (N=2048) COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 HOLROYD LOGICAL FIRST DATA FIRST/.TRUE./ C IF(FIRST) THEN HOLROYD=EOEPS*(PMOB+XMOB) FIRST=.FALSE. ENDIF RCOMX=HOLROYD/DZ**2 C RCOMX=RCOMX*ESCAL ! Don't scale by ESCAL because C ! QP and QX don't scale with ESCAL C RETURN END REAL*8 FUNCTION PROB(I) C HV plateau curve, normalized to unity at infinite E-field C Returns one minus probability of initial recombination (geminate?) C Or probability of avoiding or surviving initial recombination PARAMETER (N=2048) COMMON /FIELD/ EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) REAL*8 E,E0 REAL*8 ALPHA,E1,E2 C EZI=ABS(EZ(I)) IF(IV(7).EQ.-1) THEN !Basic PROB (simple HV plateau curve) E0=VAR(7) PROB=1. IF(EZI.LT.E0) PROB=EZI/E0 ELSEIF(IV(7).EQ.0) THEN !Thomas and Imel parameterization E=EZI/VAR(7) PROB=E*DLOG(1.+1./E) ELSEIF(IV(7).EQ.1) THEN !New Scheme ALPHA=0.86 E1=7.7/0.269 E2=900./0.269 PROB=ALPHA/(1.+E1/EZI)+(1.-ALPHA)/(1.+E2/EZI) ELSEIF(IV(7).EQ.2) THEN !exponential ALPHA=0.86 E1=7.7/0.269 E2=900./0.269 PROB=ALPHA/(1.+E1/EZI)+(1.-ALPHA)/(1.+E2/EZI) PROB=(1.+3.*EXP(-0.269*EZI/5.))*PROB ELSEIF(IV(7).EQ.3) THEN !exponential modified coefficient ALPHA=0.86 E1=7.7/0.269 E2=900./0.269 PROB=ALPHA/(1.+E1/EZI)+(1.-ALPHA)/(1.+E2/EZI) PROB=(1.+2.*EXP(-0.269*EZI/5.))*PROB ELSEIF(IV(7).EQ.4) THEN !exponential modified denominator ALPHA=0.86 E1=7.7/0.269 E2=900./0.269 PROB=ALPHA/(1.+E1/EZI)+(1.-ALPHA)/(1.+E2/EZI) PROB=(1.+2.*EXP(-0.269*EZI/2.5))*PROB ELSEIF(IV(7).EQ.5) THEN !New Scheme ALPHA=0.861 E1=9.20/0.269 E2=1171./0.269 PROB=ALPHA/(1.+E1/EZI)+(1.-ALPHA)/(1.+E2/EZI) ELSEIF(IV(7).EQ.6) THEN !New Scheme ALPHA=0.000 BETA=0.861 E1=4.0 E2=32.0 E3=4000. PROB=ALPHA/(1.+E1/EZI) PROB=PROB+(1.-ALPHA)*BETA/(1.+E2/EZI) PROB=PROB+(1.-ALPHA)*(1.-BETA)/(1.+E3/EZI) C Good guess for plateau curve below critical E-field for 2 mCi Obninsk source C for all three fills. ELSEIF(IV(7).EQ.7) THEN ALPHA=0.25 BETA=0.80 E1=2.0 E2=35.0 E3=1500. PROB= ALPHA /(1.+E1/EZI) PROB=PROB+(1.-ALPHA)*BETA /(1.+E2/EZI) PROB=PROB+(1.-ALPHA)*(1.-BETA)/(1.+E3/EZI) ELSEIF(IV(7).EQ.8) THEN ALPHA=0.22 BETA=0.80 E1=1.5 E2=35.0 E3=1500. PROB= ALPHA /(1.+E1/EZI) PROB=PROB+(1.-ALPHA)*BETA /(1.+E2/EZI) PROB=PROB+(1.-ALPHA)*(1.-BETA)/(1.+E3/EZI) ELSEIF(IV(7).EQ.9) THEN ALPHA=0.3 BETA=0.80 E1=1.5 E2=30.0 E3=1500. PROB= ALPHA /(1.+E1/EZI) PROB=PROB+(1.-ALPHA)*BETA /(1.+E2/EZI) PROB=PROB+(1.-ALPHA)*(1.-BETA)/(1.+E3/EZI) ELSE PRINT 100 100 FORMAT(' IV(7) is not properly set, must be an integer -1-8.') ENDIF RETURN END