C FUNCDRFT.FOR FUNCDRFT.FOR C Flow of charges in a liquid argon gap C 06 September 1999 C 04 Feb 2011 for Obninsk Plateau pub C 25 May 2011 Add O2 transport C 01 Mar 2012 Move parameterizations to a separate file params.for C 07 May 2012 Lots of improvements, corrections, etc. C 12 Jan 2013 Add minority positive charge carrier C 03 Aug 2013 Make compatible with hilumdrift C 25 Sep 2013 Get rid of hilmdrft.for. Use only funcdrft.for. C 16 Mar 2014 Neaten up the code a bit and normalize PROB to 2200 V/mm C23456789012345678901234567890123456789012345678901234567890123456789012 C 1 2 3 4 5 6 7 FUNCTION FUNCT(X) C COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) DATA IV/0,-1,0,1024,0,0,5,0,1,6*0/ DATA NV/10/,VAR/0.,0.,0.25,256.,0.,1.E5,0.001,0.2,0.001,0., 1 5*0./,DVR/15*0./ COMMON /FNAM/FNAME(2) DATA FNAME/4HFUNC,4HDRFT/ C PARAMETER (N=2048) !Divide gap up into N pieces PARAMETER (N1=N+1,M=2*N1) COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) C ND is not used. It pads the beginning of the common block so C that addresses of the REAL*8 variables are on byte addresses C which are multiples of 4. REAL*8 QP,QB,QM,QX,QG,QQ COMMON /CHRGDV/QPZ(0:1,0:N),QBZ(0:1,0:N),QMZ(0:1,0:N), 1 QXZ(0:1,0:N),QGZ(0:1,0:N) REAL*8 QPZ,QBZ,QMZ,QXZ,QGZ COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /CRNTS/JP(0:N),JB(0:N),JM(0:N),JX(0:N),JG(0:N),JQ(0:N) REAL*8 JP,JB,JM,JX,JG,JQ COMMON /JPHY/JDP(0:N),JDB(0:N),JDM(0:N),JDX(0:N), 1 JDG(0:N),JDT(0:N) REAL*8 JDP,JDB,JDM,JDX,JDG,JDT DATA NN/0/ ! Initial value DATA QP,QB,QM,QX,QG,QQ 1 /N*0.D0,N*0.D0,N*0.D0,N*0.D0,N*0.D0,N*0.D0/ DATA QPZ,QBZ,QMZ,QXZ,QGZ/M*0.D0,M*0.D0,M*0.D0,M*0.D0,M*0.D0/ DATA JP,JB,JM,JX,JG,JQ 1 /N1*0.D0,N1*0.D0,N1*0.D0,N1*0.D0,N1*0.D0,N1*0.D0/ COMMON /SLOPES/QPS(0:1,0:N),QBS(0:1,0:N),QMS(0:1,0:N), 1 QXS(0:1,0:N),QGS(0:1,0:N) REAL*8 QPS,QBS,QMS,QXS,QGS DATA QPS,QBS,QMS,QXS,QGS/M*0.D0,M*0.D0,M*0.D0,M*0.D0,M*0.D0/ COMMON /DPHY/DQP(N),DQB(N),DQM(N),DQX(N),DQG(N),DQQ(N) REAL*8 DQP,DQB,DQM,DQX,DQG,DQQ DATA DQP,DQB,DQM,DQX,DQG,DQQ 1 /N*0.D0,N*0.D0,N*0.D0,N*0.D0,N*0.D0,N*0.D0/ COMMON /LPARS/V0,A,D0 REAL*8 V0,A,D0 DATA V0,A,D0 /250.,0.269,8.586E10/ ! V0 is in Volts A is in mm COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL DATA VE/4.2D6/ !Electron drift speed in mm/sec (2mm per 476 nsec) DATA XMOB/0.078D0/ !Negative oxygen ion mobility (mm^2/Vs) (Holroyd EMail) COMMON /CRIT/DC,RC,JC,SCRPTR,RAT REAL*8 DC,RC,JC,SCRPTR,RAT COMMON /SIG/JT,NT,SG(10000) REAL*8 SG DATA NT,SG/10000,10000*0.D0/ DATA JT/1/ C NT Number of time samples C SG Induced transient current vs. time COMMON /CURNT/MZN,NZN,IZN,IDM,EZNN,SZNN,EZN(2000) C IDM is not used. It pads the beginning of the common block so C that the following REAL*8 variables have addresses which are C multiples of 4. REAL*8 EZNN,SZNN,EZN REAL*8 EZNNI DATA MZN/2000/,NZN/0/,IZN/0/ COMMON /CURNTS/CURP,CURM,CURDP,CURDM REAL*8 CURP,CURM,CURDP,CURDM DATA CURP,CURM,CURDP,CURDM/0.D0,0.D0,0.D0,0.D0/ COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 COMMON /TRIGGER/ITRIG,IFTRIG LOGICAL ITRIG,IFTRIG DATA ITRIG,IFTRIG/.FALSE.,.FALSE./ REAL*8 VDT,JSUM REAL*8 DIF DATA RAT,VDT,SCRPTR/0.D0,0.D0,0.D0/,DIF/0.00001D0/ DATA NN/0/ LOGICAL LOG COMMON /FILENO/IFNO DATA IFNO/0/ REAL*8 DFAC,VSUM,Z1,Z2,DN,AN,S,DQT,S1,S2,VEMAX,DRIFTVEL C C VAR(1) IV(1) Meaning C C VAR(2) IV(2) Meaning C V0 Potential across gap A C .LT.0 Reinitialize and start again C .EQ.0 No change C .GT.0 Number of Iterations C VAR(3) Multiplier for DT (to improve convergence) C IV(3) Total number of iterations (Read-only, do not set this) C VAR(4) Recombination rate C IV(4) Divide gap into NN pieces C VAR(5) Meaning C 1. Set first order conditions C 2. Read in data from file C 3. Refill charges from common /chargn/ with new no. of bins C -2. Write out data to file C IV(5) Flag to turn on higher recombination rate at low E-field C VAR(6) Instantaneous ionization (10000. is min. ionizing) C IV(6) Flag to inject instantaneous ionization C VAR(7) E0 for calculation of DI C IV(7) Choice for PROB C 0 Thomas and Imel C 1 original C 2 exponential C 3 exponential modified coefficient C 4 exponential modified denominator C 5 original with new parameters C 6 three-component Birks' Law C 7 three-component Birks' Law for all fills C VAR(8) PMOB for calculation of VP C VAR(9) ESCAL for scaling the electron velocity and SCRPTR C IV(9) Choice of electron velocity parameterization C .EQ.1 Simple parameterization as in analytic calcs C .EQ.2 Mozumder inspired parameterization C VAR(10) D0 ionization rate at infinite potential C IV(10) Flag for z-dependent ionization rate DI C .EQ.0 Flat ionization rate C .EQ.1 Ionization rate rises with z C .EQ.-1 Ionization rate falls with z C VAR(11) PPM - Parts per million of O2 molecules in LAr C IV(11) Flag to turn on diffusion C .EQ.0 No diffusion C .EQ.1 Diffusion C VAR(12) Fraction of positive charge carriers which are minority (B) C VAR(13) Fraction of PMOB which equals BMOB. VAR(13)=BMOB/PMOB C C Variables C N Number of cells in the LAr gap C DZ Size of cell (mm) C DT Iteration time interval (seconds) C A Distance across gap (mm) C V0 Voltage across gap (Volts) C QP(I) Number of positive ions in cell I per unit area (1 mm^2) C QB(I) Number of minority positive charge carriers in cell I per unit area C QM(I) Number of electrons in cell I per unit area C QX(I) Number of negative O2 ions in cell I per unit area C QG(I) Number of neutral O2 molecules in cell I per unit area C QQ(I) Net number of charges in cell I per unit area C EZ(I) Electric field at cell I (V/mm) C VV(I) Potential (Volts) at cell I C JP(I) Current density due to positive ions (No. per unit area per second) C JB(I) Current density due to minority positive charge carriers (No. per unit area per second) C JM(I) Current density due to electrons (No. per unit area per second) C JX(I) Current density due to O2 ions (No. per unit area per second) C JG(I) Current density due to O2 molecules (No. per unit area per second) C JQ(I) Total current density (No. per unit area per second) C Functions C DI(ITER,I) Number of ion pairs created per cubic mm per sec at cell I C VP(I) Drift velocity of positive ions at cell I (mm/s) C VB(I) Drift velocity of minority positive charge carriers (mm/s) C VM(I) Drift velocity of electrons at cell I (mm/s) C VX(I) Drift velocity of negative Oxygen ions C LOG=.FALSE. IF(NN.NE.IV(4)) THEN NN=IV(4) IF(NN.GT.N) NN=N LOG=.TRUE. ENDIF ESCAL=VAR(9) ITER=IV(2) IF((ITER.LT.0).OR.LOG) THEN EOEPS=1.603E-19/8.854E-15/1.51 !Proton charge divided by epsilon in V-mm DZ=A/DBLE(NN) RAT=0.D0 PPM=VAR(11) C Get HV plateau value at EZ=2200 V/mm EZ(1)=2200.D0 PROB2200=PROB(1) PRINT 50, PROB2200 50 FORMAT('PROB2200: ',F10.4) C The following will be overwritten if INSAVE is called DO I=1,NN QG(I)=6.022D23/39.95D0 ! Number of LAr molecules in a gram of LAr QG(I)=QG(I)*1.396D0 ! Number of LAr molecules in a cc of LAr QG(I)=QG(I)/1.D3 ! Number of LAr molecules in a mm^3 of LAr QG(I)=QG(I)*DZ ! Number of LAr molecules in bin of DZ by 1 mm^2 in X and Y QG(I)=QG(I)*PPM/1.0D6 ! Number of O2 molecules in bin of DZ by 1 mm^2 in X and Y ENDDO JT=1 NZN=0 IZN=0 IV(3)=0 !Counts total iterations ITER=0 LOG=.TRUE. ENDIF LOG= LOG.OR.(ABS(VAR(2)-V0 ).GT.DIF) 1 .OR.(ABS(VAR(3)-VDT ).GT.DIF) 2 .OR.(ABS(VAR(4)*ESCAL-SCRPTR).GT.DIF) 3 .OR.(ABS(VAR(8)-PMOB ).GT.DIF) 4 .OR.(ABS(VAR(13)*PMOB-BMOB).GT.DIF) IF(LOG) THEN V0=VAR(2) !Potential across the gap VDT=VAR(3) !Fraction of cell in step element PMOB=VAR(8) !Positive ion mobility BMOB=VAR(13)*PMOB !Minority positive charge mobility CALL INITDIFF(IV(11)) !Initialize diffusion coefficients D0=VAR(10) !Ionization rate-density at 2200 V/mm RB=VAR(12) !Fraction of min pos RP=1.D0-RB !Fraction of Ar+ AMOB=RP*PMOB+RB*BMOB !Weighted average mobility DC=4.D0*AMOB*V0**2/EOEPS/A**4 !Critical value of DI RC=2.D0*V0/EOEPS/A**2 !Critical number density SCRPTR=VAR(4)*ESCAL RCOM=SCRPTR*DC/RC**2/DZ**2 !Recombination rate JC=DC*A !Critical current density EZ(1)=V0/A RAT=D0*PROB(1)/PROB2200/DC IF(RAT.LE.1.0D0) THEN EDZ(0)=V0/A*SQRT(1.0D0+3.0D0*RAT) ! Linear approximation in sqrt ELSE EDZ(0)=V0/A*2.0D0*RAT**0.25 ENDIF VEMAX=ABS(VM(0)) C DT=DZ/VEMAX*VDT !Time for an electron to drift across DT=DZ/VE*VDT/ESCAL ! *** Substitute for hilumdrift C !VDT cells at nominal EZ PRINT 100, RAT,DT 100 FORMAT('RAT, DT =',1P,2D15.3) DTMDZ=DT*DZ DTODZ=DT/DZ DO L=0,N DO MN=0,1 QPS(MN,L)=0.D0 QBS(MN,L)=0.D0 QMS(MN,L)=0.D0 QXS(MN,L)=0.D0 QGS(MN,L)=0.D0 ENDDO ENDDO ENDIF C IF(VAR(5).GT.0.5) THEN !Set initial conditions IF(VAR(5).LT.1.5) THEN CALL PRECON ELSEIF(VAR(5).LT.2.5) THEN CALL INSAVE CALL REFILL ELSE CALL REFILL ENDIF VAR(5)=0. ENDIF C IF(VAR(5).LT.-1.5) THEN !Record current conditions CALL SAVEOUT CALL SAVEOTHER !Save info for debugging VAR(5)=0. ENDIF C IF(ABS(IV(6)).NE.0) THEN !Inject instantaneous ionization IV(6)=0 DO I=1,NN QP(I)=QP(I)+VAR(6)*DZ*DI(1,I)/D0 !One min. ionizing particle C This adds VAR(6) ion pairs per mm to each cell minus initial recombination. QM(I)=QM(I)+VAR(6)*DZ*DI(1,I)/D0 !with geminate recombination END DO JT=2 !Reset clock and start EZNN=0. SZNN=0. IF((NZN.GT.1).AND.(MZN.GT.1)) THEN EZNNI=EZN(1) !Ploy to avoid computer errors DO I=1,NZN !Average over last NZN<=MZN values of EDZ(NN) EZNN=EZNN+(EZN(I)-EZNNI) ENDDO IF(NZN.LT.MZN) THEN SZNN=(EZN(NZN)-EZN(1))/DBLE(NZN-1) ELSEIF(IZN.EQ.MZN) THEN SZNN=(EZN(IZN)-EZN(1))/DBLE(MZN-1) ELSE SZNN=(EZN(IZN)-EZN(IZN+1))/DBLE(MZN-1) ENDIF EZNN=EZNNI+EZNN/FLOAT(NZN)+DBLE(NZN-1)/2.*SZNN ENDIF ENDIF C DFAC=1.0D0 !Factor to soften slope IF(ITER.GT.0) THEN CALL EFIELD !First time to set E-field DO K=1,ITER !Iterate the configuration CALL EFIELD !This order is ************ CALL ITERT(K,DFAC) !backwards for a test ************ IF(ITRIG) CALL REPORTI(K) ENDDO IV(3)=IV(3)+ITER VSUM=V0+EDZ(0)*DZ/2. DO I=1,NN VSUM=VSUM-EDZ(I-1)*DZ VV(I)=VSUM ENDDO CALL EDGES(DFAC) DO I=0,NN DRIFTVEL=VP(I)+VPD(I) JLR=1 IF(DRIFTVEL.GE.0.0) JLR=0 JDP(I)=+DRIFTVEL*QPZ(JLR,I)/DZ DRIFTVEL=VB(I)+VBD(I) JLR=1 IF(DRIFTVEL.GE.0.0) JLR=0 JDB(I)=+DRIFTVEL*QBZ(JLR,I)/DZ DRIFTVEL=VM(I)+VMD(I) JLR=1 IF(DRIFTVEL.GE.0.0) JLR=0 JDM(I)=-DRIFTVEL*QMZ(JLR,I)/DZ DRIFTVEL=VX(I)+VXD(I) JLR=1 IF(DRIFTVEL.GE.0.0) JLR=0 JDX(I)=-DRIFTVEL*QXZ(JLR,I)/DZ DRIFTVEL=VGD(I) JLR=1 IF(DRIFTVEL.GE.0.0) JLR=0 JDG(I)=+DRIFTVEL*QGZ(JLR,I)/DZ JDT(I)=JDP(I)+JDM(I)+JDX(I) ENDDO ENDIF IV(2)=0 C RETURN END SUBROUTINE READINPARMS C Read in gap potential, gap width, and ionization rate COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /LPARS/V0,A,D0 REAL*8 V0,A,D0 CHARACTER*80 ASC REAL MUPLUS ! Ar+ mobility in mm^2/Vs - Favorite value is 0.1 C OPEN(UNIT=1,FILE='GapPotential.dat',STATUS='OLD',ERR=200) READ(1,100) ASC READ(1,100) ASC READ(1,100) ASC 100 FORMAT(A80) READ(1,110) V0,A,D0,MUPLUS,SCRPTR 110 FORMAT(E16.0,4E15.0) CLOSE(UNIT=1) VAR(2)=V0 VAR(4)=SCRPTR VAR(8)=MUPLUS VAR(10)=D0 C PRINT 120, VAR(2),A,VAR(10),VAR(8),VAR(4) C120 FORMAT(' V0,A,D0,MUPLUS,SCRPTR:',5E15.5) RETURN 200 PRINT 210 210 FORMAT(' Cannot read in GapPotential.dat') C RETURN END SUBROUTINE EFIELD C Recalculate E field. First do the constant term. C EZ is field at cell center, EDZ is field at cell boundaries. 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 /DPHY/DQP(N),DQB(N),DQM(N),DQX(N),DQG(N),DQQ(N) REAL*8 DQP,DQB,DQM,DQX,DQG,DQQ REAL*8 DSUM,SUM(N),DSMM,SMM(N) REAL*8 DSUN,SUN(N),DSNN,SNN(N) COMMON /LPARS/V0,A,D0 REAL*8 V0,A,D0 COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 C1,C2,S1,S2 REAL*8 DQD(N) C SUM(1)=QQ(1)/2.D0 SMM(1)=QQ(1) DSUM=SUM(1) DSMM=SMM(1) DQD(1)=DQP(1)-DQM(1) SUN(1)=DQD(1)/2.D0 SNN(1)=DQD(1) DSUN=SUN(1) DSNN=SNN(1) DO I=2,NN SUM(I)=SUM(I-1)+(QQ(I-1)+QQ(I))/2.D0 DSUM=DSUM+SUM(I) SMM(I)=SMM(I-1)+QQ(I) DSMM=DSMM+SMM(I) DQD(I)=DQP(I)-DQM(I) SUN(I)=SUN(I-1)+(DQD(I-1)+DQD(I))/2.D0 DSUN=DSUN+SUN(I) SNN(I)=SNN(I-1)+DQD(I) DSNN=DSNN+SNN(I) ENDDO DSMM=DSMM-SMM(NN)/2.D0 DSNN=DSNN-SNN(NN)/2.D0 C1=V0/A-DSUM*EOEPS*DZ/A C2=V0/A-DSMM*EOEPS*DZ/A S1= -DSUN*EOEPS*DZ/A S2= -DSNN*EOEPS*DZ/A C Now do the E field EDZ(0)=C2 DEZ(0)=S2 DO I=1,NN EZ(I) =C1+SUM(I)*EOEPS EDZ(I)=C2+SMM(I)*EOEPS DEZ(I)=S2+SNN(I)*EOEPS ENDDO C RETURN END SUBROUTINE ITERT(ITER,FRAC) C One iteration. Each iteration is a time step of DT. REAL*8 FRAC COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) 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 /DPHY/DQP(N),DQB(N),DQM(N),DQX(N),DQG(N),DQQ(N) REAL*8 DQP,DQB,DQM,DQX,DQG,DQQ COMMON /CHRGDV/QPZ(0:1,0:N),QBZ(0:1,0:N),QMZ(0:1,0:N), 1 QXZ(0:1,0:N),QGZ(0:1,0:N) REAL*8 QPZ,QBZ,QMZ,QXZ,QGZ COMMON /DCSC/DCP(N),DCB(N),DCM(N),DCX(N),DCG(N) REAL*8 DCP,DCB,DCM,DCX,DCG COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /CRNTS/JP(0:N),JB(0:N),JM(0:N),JX(0:N),JG(0:N),JQ(0:N) REAL*8 JP,JB,JM,JX,JG,JQ COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 COMMON /CURNTS/CURP,CURM,CURDP,CURDM REAL*8 CURP,CURM,CURDP,CURDM COMMON /CRIT/DC,RC,JC,SCRPTR,RAT REAL*8 DC,RC,JC,SCRPTR,RAT REAL*8 PART DATA PART/0.5D0/ REAL*8 DII,RP,RB,NPM,NPX,NBM,NBX,XMG C CURP=0. CURM=0. CURDP=0. CURDM=0. RB=VAR(12) ! Fraction of positive charge carriers which are minority RP=1.D0-RB DO I=1,NN NPM=RCOMB(I)*QP(I)*QM(I) ! Make neutral from QP and QM NPX=RCOMX(I)*QP(I)*QX(I) ! Make Neut O2 from QP and QX NBM=RCMBB(I)*QB(I)*QM(I) ! Make neutral from QB and QM NBX=RCMXB(I)*QB(I)*QX(I) ! Make Neut O2 from QB and QX XMG=RATCH(I)*QM(I)*QG(I) ! Make O2 ion from QM and QG DII=DI(ITER,I) DCP(I)=DTMDZ*(RP*DII-NPM-NPX) DCB(I)=DTMDZ*(RB*DII-NBM-NBX) DCM(I)=DTMDZ*( DII-XMG-NPM-NBM) DCX(I)=DTMDZ*(XMG-NPX-NBX) DCG(I)=DTMDZ*(NPX+NBX-XMG) CURDP=CURDP+DCP(I)+DCB(I) !Positive charge created in the gap CURDM=CURDM+DCM(I)+DCX(I) !Negative charge created in the gap ENDDO DO I=1,NN QP(I)=QP(I)+PART*DCP(I) QB(I)=QB(I)+PART*DCB(I) QM(I)=QM(I)+PART*DCM(I) QX(I)=QX(I)+PART*DCX(I) QG(I)=QG(I)+PART*DCG(I) QQ(I)=QP(I)+QB(I)-QM(I)-QX(I) ENDDO C IF((ITER/100)*100.EQ.ITER) CALL SLOPE(FRAC) CALL DRIFT(0.0D0) DO I=1,NN QP(I)=QP(I)+(1.D0-PART)*DCP(I)+DQP(I) IF(QP(I).LT.1.D-6) QP(I)=0.D0 QB(I)=QB(I)+(1.D0-PART)*DCB(I)+DQB(I) IF(QB(I).LT.1.D-6) QB(I)=0.D0 QM(I)=QM(I)+(1.D0-PART)*DCM(I)+DQM(I) IF(QM(I).LT.1.D-6) QM(I)=0.D0 QX(I)=QX(I)+(1.D0-PART)*DCX(I)+DQX(I) IF(QX(I).LT.1.D-6) QX(I)=0.D0 QG(I)=QG(I)+(1.D0-PART)*DCG(I)+DQG(I) IF(QG(I).LT.1.D-6) QG(I)=0.D0 QQ(I)=QP(I)+QB(I)-QM(I)-QX(I) ENDDO CURP=CURP/DT CURM=CURM/DT CURDP=CURDP/DT+JP(0)+JB(0) CURDM=CURDM/DT+JM(NN)+JX(NN) NREP=100000 IF((ITER/NREP)*NREP.EQ.ITER) 1 PRINT 100, ITER,CURP,CURM,CURDP,CURDM 100 FORMAT(' For ',I8,'th iteration,', 1 ' current at left/right/created in gap =', 2 1P,4D12.5) CALL EDGES(FRAC) DO I=0,NN JP(I)=+(VP(I)+VPD(I))*QPZ(0,I)/DZ JB(I)=+(VB(I)+VBD(I))*QBZ(0,I)/DZ JM(I)=-(VM(I)+VMD(I))*QMZ(1,I)/DZ JX(I)=-(VX(I)+VXD(I))*QXZ(1,I)/DZ JG(I)=+( +VGD(I))*QGZ(0,I)/DZ JQ(I)=JP(I)+JB(I)+JM(I)+JX(I) ENDDO C RETURN END SUBROUTINE DRIFT(DFAC) C Drift charges within each cell. REAL*8 DFAC 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 /SLOPES/QPS(0:1,0:N),QBS(0:1,0:N),QMS(0:1,0:N), 1 QXS(0:1,0:N),QGS(0:1,0:N) REAL*8 QPS,QBS,QMS,QXS,QGS COMMON /CRNTS/JP(0:N),JB(0:N),JM(0:N),JX(0:N),JG(0:N),JQ(0:N) REAL*8 JP,JB,JM,JX,JG,JQ COMMON /DPHY/DQP(N),DQB(N),DQM(N),DQX(N),DQG(N),DQQ(N) REAL*8 DQP,DQB,DQM,DQX,DQG,DQQ COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /SIG/JT,NT,SG(10000) REAL*8 SG COMMON /CURNT/MZN,NZN,IZN,IDM,EZNN,SZNN,EZN(2000) REAL*8 EZNN,SZNN,EZN COMMON /CURNTS/CURP,CURM,CURDP,CURDM REAL*8 CURP,CURM,CURDP,CURDM INTEGER UD,ZF,ZT REAL*8 DFZ,ADFZ,DQT,VZ C DO I=1,NN DQP(I)=0.D0 DQB(I)=0.D0 DQM(I)=0.D0 DQX(I)=0.D0 DQG(I)=0.D0 DQQ(I)=0.D0 ENDDO C Now drift for time DT across boundaries (edges) between cells DO I=1,NN-1 C First do positive Argon ions UD=1 DFZ=(VP(I)+VPD(I))*DTODZ IF(DFZ.GT.0.) UD=0 ZF=I+UD !Z-from ZT=I+1-UD !Z-to ADFZ=ABS(DFZ) DQT=ADFZ*QP(ZF)*(1.D0+QPS(UD,I)*(1.D0-DFAC*ADFZ)) DQP(ZF)=DQP(ZF)-DQT DQP(ZT)=DQP(ZT)+DQT C Now do minority positives UD=1 DFZ=(VB(I)+VBD(I))*DTODZ IF(DFZ.GT.0.) UD=0 ZF=I+UD !Z-from ZT=I+1-UD !Z-to ADFZ=ABS(DFZ) DQT=ADFZ*QB(ZF)*(1.D0+QBS(UD,I)*(1.D0-DFAC*ADFZ)) DQB(ZF)=DQB(ZF)-DQT DQB(ZT)=DQB(ZT)+DQT C Now do ionization electrons UD=1 DFZ=(VM(I)+VMD(I))*DTODZ IF(DFZ.GT.0.) UD=0 ZF=I+UD ZT=I+1-UD ADFZ=ABS(DFZ) DQT=ADFZ*QM(ZF)*(1.D0+QMS(UD,I)*(1.D0-DFAC*ADFZ)) DQM(ZF)=DQM(ZF)-DQT DQM(ZT)=DQM(ZT)+DQT C Next do negative O2 ions UD=1 DFZ=(VX(I)+VXD(I))*DTODZ IF(DFZ.GT.0.) UD=0 ZF=I+UD ZT=I+1-UD ADFZ=ABS(DFZ) DQT=ADFZ*QX(ZF)*(1.D0+QXS(UD,I)*(1.D0-DFAC*ADFZ)) DQX(ZF)=DQX(ZF)-DQT DQX(ZT)=DQX(ZT)+DQT C Finally do neutral O2 molecules UD=1 DFZ=VGD(I)*DTODZ IF(DFZ.GT.0.) UD=0 ZF=I+UD ZT=I+1-UD ADFZ=ABS(DFZ) DQT=ADFZ*QG(ZF)*(1.D0+QGS(UD,I)*(1.D0-DFAC*ADFZ)) DQG(ZF)=DQG(ZF)-DQT DQG(ZT)=DQG(ZT)+DQT ENDDO C Now do the special cases at z=0 and z=a (I=0, I=NN). C First do positive Argon ions VZ=VP(0)+VPD(0) JP(0)=0.D0 IF(VZ.LT.0.) THEN DFZ=-VZ*DTODZ DQT=DFZ*QP(1)*(1.D0-QPS(1,0)*(1.D0-DFAC*DFZ)) DQP(1)=DQP(1)-DQT JP(0)=-DQT CURP=CURP-DQT ENDIF VZ=VP(NN)+VPD(NN) JP(NN)=0.D0 IF(VZ.GE.0.) THEN DFZ=VZ*DTODZ DQT=DFZ*QP(NN)*(1.D0-QPS(0,NN)*(1.D0-DFAC*DFZ)) DQP(NN)=DQP(NN)-DQT JP(NN)=+DQT CURM=CURM+DQT ENDIF c Now do minority positives VZ=VB(0)+VBD(0) JB(0)=0.D0 IF(VZ.LT.0.) THEN DFZ=-VZ*DTODZ DQT=DFZ*QB(1)*(1.D0-QBS(1,0)*(1.D0-DFAC*DFZ)) DQB(1)=DQB(1)-DQT JB(0)=-DQT CURP=CURP-DQT ENDIF VZ=VB(NN)+VBD(NN) JB(NN)=0.D0 IF(VZ.GE.0.) THEN DFZ=VZ*DTODZ DQT=DFZ*QB(NN)*(1.D0-QBS(0,NN)*(1.D0-DFAC*DFZ)) DQB(NN)=DQB(NN)-DQT JB(NN)=+DQT CURM=CURM+DQT ENDIF C Now do ionization electrons VZ=VM(0)+VMD(0) JM(0)=0.D0 IF(VZ.LT.0.) THEN DFZ=-VZ*DTODZ DQT=DFZ*QM(1)*(1.D0-QMS(1,0)*(1.D0-DFAC*DFZ)) DQM(1)=DQM(1)-DQT JM(0)=+DQT CURP=CURP+DQT ENDIF VZ=VM(NN)+VMD(NN) JM(NN)=0.D0 IF(VZ.GE.0.) THEN DFZ=VZ*DTODZ DQT=DFZ*QM(NN)*(1.D0-QMS(0,NN)*(1.D0-DFAC*DFZ)) DQM(NN)=DQM(NN)-DQT JM(NN)=-DQT CURM=CURM-DQT ENDIF C Next do negative O2 ions and neutral O2 molecules VZ=VX(0)+VXD(0) JX(0)=0.D0 JG(0)=0.D0 IF(VZ.LT.0.) THEN DFZ=-VZ*DTODZ DQT=DFZ*QX(1)*(1.D0-QXS(1,0)*(1.D0-DFAC*DFZ)) DQX(1)=DQX(1)-DQT DQG(1)=DQG(1)+DQT ! Assume neutral O2 enters from cathode JX(0)=+DQT JG(0)=+DQT CURP=CURP+DQT ENDIF VZ=VX(NN)+VXD(NN) JX(NN)=0.D0 JG(NN)=0.D0 IF(VZ.GE.0.) THEN DFZ=VZ*DTODZ DQT=DFZ*QX(NN)*(1.D0-QXS(0,NN)*(1.D0-DFAC*DFZ)) DQX(NN)=DQX(NN)-DQT DQG(NN)=DQG(NN)+DQT JX(NN)=-DQT JG(NN)=-DQT CURM=CURM-DQT ENDIF C Tally the sum of all changes in the charge in each cell due C to charges drifting into and out of each cell. Ionization C and recombination always conserve the charge within a cell C so drifting is the only mechanism which changes the charge C accummulating in a cell. DO I=1,NN DQQ(I)=DQP(I)+DQB(I)-DQM(I)-DQX(I) ENDDO C Record induced current IF(JT.EQ.1) THEN NZN=NZN+1 IF(NZN.GT.MZN) NZN=MZN IZN=IZN+1 IF(IZN.GT.MZN) IZN=1 EZN(IZN)=EDZ(NN) !Save the last NZN<=MZN values of EDZ(NN) SG(1)=0. ELSEIF((JT.GT.1).AND.(JT.LE.NT)) THEN SG(JT)=EDZ(NN)-EZNN-FLOAT(JT-1)*SZNN !Induced current at z=a. JT=JT+1 ENDIF C RETURN END SUBROUTINE SLOPE(F) C Cells are numbered 1 to NN C Edges are numbered 0 to NN C Edge I is between cell I and cell I+1 C QPS(0,I) is relative change from center of cell I to edge I C QPS(1,I) is relative change from center of cell I+1 to edge I C F is a parameter to soften the slope to avoid instabilities. C QPS(UD,ZF/T) ZF=I+UD ZT=I+1-UD REAL*8 F 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 /SLOPES/QPS(0:1,0:N),QBS(0:1,0:N),QMS(0:1,0:N), 1 QXS(0:1,0:N),QGS(0:1,0:N) REAL*8 QPS,QBS,QMS,QXS,QGS INTEGER UD REAL*8 QT0,QT1,QTMP0,QTMP1 C C First work at boundaries between cells DO I=1,NN-1 C Ar+ ions QPS(0,I)=(QP(I+1)/QP(I)-1.D0)/2.D0 QPS(1,I)=(QP(I)/QP(I+1)-1.D0)/2.D0 C Minority positives QBS(0,I)=(QB(I+1)/QB(I)-1.D0)/2.D0 QBS(1,I)=(QB(I)/QB(I+1)-1.D0)/2.D0 C Ionization electrons QMS(0,I)=(QM(I+1)/QM(I)-1.D0)/2.D0 QMS(1,I)=(QM(I)/QM(I+1)-1.D0)/2.D0 C Negative O2 ions QXS(0,I)=(QX(I+1)/QX(I)-1.D0)/2.D0 QXS(1,I)=(QX(I)/QX(I+1)-1.D0)/2.D0 C Neutral O2 molecules QGS(0,I)=(QG(I+1)/QG(I)-1.D0)/2.D0 QGS(1,I)=(QG(I)/QG(I+1)-1.D0)/2.D0 ENDDO C Now do the special cases at z=0 and z=a C Ar+ ions QPS(0,0)=0.D0 !Not used QT1=-1.D0 QTMP0=QP(1) QTMP1=QP(2) IF(QTMP1.LT.3.D0*QTMP0) QT1=(1.D0-QTMP1/QTMP0)/2.D0 QPS(1,0)=QT1 QPS(1,NN)=0.D0 !Not used QT0=-1.D0 QTMP0=QP(NN-1) QTMP1=QP(NN) IF(QTMP0.LT.3.D0*QTMP1) QT0=(1.D0-QTMP0/QTMP1)/2.D0 QPS(0,NN)=QT0 C Minority positives QBS(0,0)=0.D0 !Not used QT1=-1.D0 QTMP0=QB(1) QTMP1=QB(2) IF(QTMP1.LT.3.D0*QTMP0) QT1=(1.D0-QTMP1/QTMP0)/2.D0 QBS(1,0)=QT1 QBS(1,NN)=0.D0 !Not used QT0=-1.D0 QTMP0=QB(NN-1) QTMP1=QB(NN) IF(QTMP0.LT.3.D0*QTMP1) QT0=(1.D0-QTMP0/QTMP1)/2.D0 QBS(0,NN)=QT0 C Ionization electrons QMS(0,0)=0.D0 !Not used QT1=-1.D0 QTMP0=QM(1) QTMP1=QM(2) IF(QTMP1.LT.3.D0*QTMP0) QT1=(1.D0-QTMP1/QTMP0)/2.D0 QMS(1,0)=QT1 QMS(1,NN)=0.D0 !Not used QT0=-1.D0 QTMP0=QM(NN-1) QTMP1=QM(NN) IF(QTMP0.LT.3.D0*QTMP1) QT0=(1.D0-QTMP0/QTMP1)/2.D0 QMS(0,NN)=QT0 C Negative O2 ions QXS(0,0)=0.D0 !Not used QT1=-1.D0 QTMP0=QX(1) QTMP1=QX(1) IF(QTMP1.LT.3.D0*QTMP0) QT1=(1.D0-QTMP1/QTMP0)/2.D0 QXS(1,0)=QT1 QXS(1,NN)=0.D0 !Not used QT0=-1.D0 QTMP0=QX(NN-1) QTMP1=QX(NN) IF(QTMP0.LT.3.D0*QTMP1) QT0=(1.D0-QTMP0/QTMP1)/2.D0 QXS(0,NN)=QT0 C Neutral O2 molecules QGS(0,0)=0.D0 !Not used QT1=-1.D0 QTMP0=QG(1) QTMP1=QG(2) IF(QTMP1.LT.3.D0*QTMP0) QT1=(1.D0-QTMP1/QTMP0)/2.D0 QGS(1,0)=QT1 QGS(1,NN)=0.D0 !Not used QT0=-1.D0 QTMP0=QG(NN-1) QTMP1=QG(NN) IF(QTMP0.LT.3.D0*QTMP1) QT0=(1.D0-QTMP0/QTMP1)/2.D0 QGS(0,NN)=QT0 C Now scale slopes to ease convergence DO I=0,NN DO UD=0,1 QPS(UD,I)=F*QPS(UD,I) !!! First and only appearance of 'F' in the program. QBS(UD,I)=F*QBS(UD,I) QMS(UD,I)=F*QMS(UD,I) QXS(UD,I)=F*QXS(UD,I) QGS(UD,I)=F*QGS(UD,I) ENDDO ENDDO C RETURN END SUBROUTINE EDGES(DF) C Calculate charge in cell extrapolated to the cell edge C Cells are numbered 1 to NN C Edges are numbered 0 to NN C Edge I is between cell I and cell I+1 C QPZ(0,I) is charge just to the left of edge I C QPZ(1,I) is charge just to the right of edge I C DF softens the calculation on top of some softening in C subroutine slope REAL*8 DF 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 /CHRGDV/QPZ(0:1,0:N),QBZ(0:1,0:N),QMZ(0:1,0:N), 1 QXZ(0:1,0:N),QGZ(0:1,0:N) REAL*8 QPZ,QBZ,QMZ,QXZ,QGZ COMMON /SLOPES/QPS(0:1,0:N),QBS(0:1,0:N),QMS(0:1,0:N), 1 QXS(0:1,0:N),QGS(0:1,0:N) REAL*8 QPS,QBS,QMS,QXS,QGS C DO I=1,NN-1 QPZ(0,I)=QP(I)*(1.D0+DF*QPS(0,I)) QPZ(1,I)=QP(I+1)*(1.D0+DF*QPS(1,I)) QBZ(0,I)=QB(I)*(1.D0+DF*QBS(0,I)) QBZ(1,I)=QB(I+1)*(1.D0+DF*QBS(1,I)) QMZ(0,I)=QM(I)*(1.D0+DF*QMS(0,I)) QMZ(1,I)=QM(I+1)*(1.D0+DF*QMS(1,I)) QXZ(0,I)=QX(I)*(1.D0+DF*QXS(0,I)) QXZ(1,I)=QX(I+1)*(1.D0+DF*QXS(1,I)) QGZ(0,I)=QG(I)*(1.D0+DF*QGS(0,I)) QGZ(1,I)=QG(I+1)*(1.D0+DF*QGS(1,I)) ENDDO QPZ(0,0)=0. QPZ(1,0)=QP(1)*(1.D0+DF*QPS(1,0)) QBZ(0,0)=0. QBZ(1,0)=QB(1)*(1.D0+DF*QBS(1,0)) QMZ(0,0)=0. QMZ(1,0)=QM(1)*(1.D0+DF*QMS(1,0)) QXZ(0,0)=0. QXZ(1,0)=QX(1)*(1.D0+DF*QXS(1,0)) QGZ(0,0)=0 QGZ(1,0)=QG(1)*(1.D0+DF*QGS(1,0)) QPZ(0,NN)=QP(NN)*(1.D0+DF*QPS(0,NN)) QPZ(1,NN)=0. QBZ(0,NN)=QB(NN)*(1.D0+DF*QBS(0,NN)) QBZ(1,NN)=0. QMZ(0,NN)=QM(NN)*(1.D0+DF*QMS(0,NN)) QMZ(1,NN)=0. QXZ(0,NN)=QX(NN)*(1.D0+DF*QXS(0,NN)) QXZ(1,NN)=0. QGZ(0,NN)=QG(NN)*(1.D0+DF*QGS(0,NN)) QGZ(1,NN)=0. C RETURN END SUBROUTINE PRECON C Institute 1st iterative solution PARAMETER (N=2048) !Divide gap up into N pieces 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 /CHRGDV/QPZ(0:1,0:N),QBZ(0:1,0:N),QMZ(0:1,0:N), 1 QXZ(0:1,0:N),QGZ(0:1,0:N) REAL*8 QPZ,QBZ,QMZ,QXZ,QGZ COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /CRNTS/JP(0:N),JB(0:N),JM(0:N),JX(0:N),JG(0:N),JQ(0:N) REAL*8 JP,JB,JM,JX,JG,JQ COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /CRIT/DC,RC,JC,SCRPTR,RAT REAL*8 DC,RC,JC,SCRPTR,RAT COMMON /LPARS/V0,A,D0 REAL*8 V0,A,D0 COMMON /VEL/VE,EMOB,XMOB,PMOB,BMOB,ESCAL REAL*8 VE,EMOB,XMOB,PMOB,BMOB,ESCAL C 1 2 3 4 5 6 7 REAL*8 KR,RK,SKR,Z0,Z1,SR,FR,ZC,DZ,RMOB,ZA0,ZA1,SRZA1,SQR0,SQR1 REAL*8 ZR0,ZR1,RP,RB REAL*8 SUMI1 C RB=VAR(12) RP=1.D0-RB AMOB=RP*PMOB+RB*BMOB RPJ=RP*PMOB/AMOB RBJ=RB*BMOB/AMOB SR=SQRT(RAT) FR=SQRT(SR) ZC=0.D0 IF(RAT.GT.1.D0) ZC=(1.D0-1.D0/FR)*A DZ=A/FLOAT(NN) RMOB=PMOB/EMOB/ESCAL IF(RAT.LT.1.D0) THEN KR=RK(RAT) SKR=SQRT(KR) ENDIF DO I=0,NN Z0=FLOAT(I)*A/FLOAT(NN) Z1=(FLOAT(I)-0.5D0)*A/FLOAT(NN) ZA0=Z0/A ZA1=Z1/A ZR0=(Z0-ZC)/(A-ZC) ZR1=(Z1-ZC)/(A-ZC) IF(RAT.LT.1.D0) THEN SRZA1=SR*ZA1 SQR0=SQRT(RAT*ZA0**2+KR) SQR1=SQRT(RAT*ZA1**2+KR) IF(I.GT.0) THEN IF(RAT.GT.0.D0) THEN VV(I)=V0*(1.D0+KR/SR*DLOG(SKR/(SRZA1+SQR1))-ZA1*SQR1) ELSE VV(I)=V0*(1.D0-ZA1) ENDIF EZ(I)=2.D0*(V0/A)*SQR1 QP(I)=RP*RAT*RC*ZA1/SQR1*DZ QB(I)=RB*RAT*RC*ZA1/SQR1*DZ QM(I)=RAT*JC/VE/ESCAL*(1.D0-ZA1)*DZ QQ(I)=QP(I)+QB(I)-QM(I) ENDIF EDZ(I)=2.D0*(V0/A)*SQR0 JP(I)=RPJ*RAT*JC*ZA0 JB(I)=RBJ*RAT*JC*ZA0 JM(I)=RAT*JC*(1.-ZA0) JQ(I)=RAT*JC ELSEIF(RAT.EQ.1.D0) THEN IF(I.GT.0) THEN VV(I)=V0*(1.D0-ZA1**2) EZ(I)=2.D0*(V0/A)*ZA1 QP(I)=RP*RC*DZ QB(I)=RB*RC*DZ QM(I)=JC/VE/ESCAL*(1.D0-ZA1)*DZ QQ(I)=QP(I)+QB(I)-QM(I) ENDIF EDZ(I)=2.D0*(V0/A)*ZA0 JP(I)=RPJ*JC*ZA0 JB(I)=RBJ*JC*ZA0 JM(I)=JC*(1.D0-ZA0) JQ(I)=JC ELSE IF(Z0.LT.ZC) THEN IF(I.GT.0) THEN VV(I)=V0 EZ(I)=2.D0*FR*SQRT(SCRPTR)*RMOB*(V0/A) !Replaced missing factor of 2 QP(I)=RP*SQRT(RAT/SCRPTR)*RC*DZ QB(I)=RB*SQRT(RAT/SCRPTR)*RC*DZ QM(I)=SQRT(RAT/SCRPTR)*RC*DZ QQ(I)=0.D0 ENDIF EDZ(I)=2.D0*FR*SQRT(SCRPTR)*RMOB*(V0/A) JP(I)=RPJ*SR*FR*JC*RMOB JB(I)=RBJ*SR*FR*JC*RMOB JM(I)=SR*FR*JC*(1.D0-RMOB) JQ(I)=SR*FR*JC ELSE IF(I.GT.0) THEN VV(I)=V0*(1.-ZR1**2) EZ(I)=2.D0*FR*(V0/A)*ZR1 QP(I)=RP*SR*RC*DZ QB(I)=RB*SR*RC*DZ QM(I)=SR*FR*JC/VE/ESCAL*(1.D0-ZR1)*DZ QQ(I)=QP(I)+QB(I)-QM(I) ENDIF EDZ(I)=2.D0*FR*(V0/A)*ZR0 JP(I)=RPJ*SR*FR*JC*ZR0 JB(I)=RBJ*SR*FR*JC*ZR0 JM(I)=SR*FR*JC*(1.D0-ZR0) JQ(I)=SR*FR*JC ENDIF ENDIF END DO C Check results SUMI1=0.D0 DO I=0,NN SUMI1=SUMI1+JQ(I) ENDDO SUMI1=SUMI1/FLOAT(NN+1) PRINT 100, SUMI1 100 FORMAT(' In PRECON ave current is ',1PD12.5) C RETURN END REAL*8 FUNCTION RK(R) C Find rk given r REAL*8 R REAL*8 SK,FK DATA ITER/10/ C IF(R.EQ.0.D0) THEN RK=0.25D0 RETURN ENDIF RK=(1.D0-R)/4.D0 !First guess DO I=1,ITER SK=RK/R RK=FK(SK) ENDDO C RETURN END REAL*8 FUNCTION FK(SK) C Find rk given k REAL*8 SK REAL*8 V,SV,SV1,T C V=1./SK SV=SQRT(V) SV1=SQRT(1.+V) T=SV1+DLOG(SV+SV1)/SV FK=1.D0/T**2 C RETURN END SUBROUTINE REFILL C NN is new number of bins C NNN is old number of bins PARAMETER (N=2048) COMMON /CHARGN/NNN,NDD,QPI(N),QBI(N),QMI(N),QXI(N),QGI(N),QQI(N) REAL*8 QPI,QBI,QMI,QXI,QGI,QQI COMMON /CHARGE/NN,ND,QP(N),QB(N),QM(N),QX(N),QG(N),QQ(N) REAL*8 QP,QB,QM,QX,QG,QQ LOGICAL LOG C MM=NN/NNN MR=MOD(NN,NNN) !Remainder IF((MM.LT.1).OR.(MR.NE.0)) THEN PRINT 20 20 FORMAT(' NN not divisible evenly by NNN in INSAVE') LOG=.TRUE. ENDIF IF(LOG) STOP II=0 DO I=1,NNN DO JJ=1,MM II=II+1 QP(II)=QPI(I)/DBLE(MM) QB(II)=QBI(I)/DBLE(MM) QM(II)=QMI(I)/DBLE(MM) QX(II)=QXI(I)/DBLE(MM) QG(II)=QGI(I)/DBLE(MM) QQ(II)=QQI(I)/DBLE(MM) ENDDO ENDDO C RETURN END SUBROUTINE SAVEOUT C Write out data to disk file 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 /CHARGN/NNN,NDD,QPI(N),QBI(N),QMI(N),QXI(N),QGI(N),QQI(N) REAL*8 QPI,QBI,QMI,QXI,QGI,QQI COMMON /CHRGDV/QPZ(0:1,0:N),QBZ(0:1,0:N),QMZ(0:1,0:N), 1 QXZ(0:1,0:N),QGZ(0:1,0:N) REAL*8 QPZ,QBZ,QMZ,QXZ,QGZ COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /CRNTS/JP(0:N),JB(0:N),JM(0:N),JX(0:N),JG(0:N),JQ(0:N) REAL*8 JP,JB,JM,JX,JG,JQ COMMON /SLOPES/QPS(0:1,0:N),QBS(0:1,0:N),QMS(0:1,0:N), 1 QXS(0:1,0:N),QGS(0:1,0:N) REAL*8 QPS,QBS,QMS,QXS,QGS COMMON /DPHY/DQP(N),DQB(N),DQM(N),DQX(N),DQG(N),DQQ(N) REAL*8 DQP,DQB,DQM,DQX,DQG,DQQ COMMON /LPARS/V0,A,D0 REAL*8 V0,A,D0 COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 COMMON /CRIT/DC,RC,JC,SCRPTR,RAT REAL*8 DC,RC,JC,SCRPTR,RAT COMMON /SIG/JT,NT,SG(10000) REAL*8 SG REAL*8 SGI(250) COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /CURNTS/CURP,CURM,CURDP,CURDM REAL*8 CURP,CURM,CURDP,CURDM COMMON /CURNT/MZN,NZN,IZN,IDM,EZNN,SZNN,EZN(2000) REAL*8 EZNN,SZNN,EZN CHARACTER*30 FILEN CHARACTER*5 EXT, FXT COMMON /FILENO/IFNO REAL*8 DSUM,DSMM CHARACTER*8 DATE CHARACTER*9 TIME CHARACTER*5 ZONE INTEGER DTVALUES(8) LOGICAL CHECK C C Save no. of bins and charges in temporary storage NNN=NN DO I=1,NN QPI(I)=QP(I) QBI(I)=QB(I) QMI(I)=QM(I) QXI(I)=QX(I) QGI(I)=QG(I) QQI(I)=QQ(I) ENDDO C C Get date and time CALL DATE_AND_TIME(DATE,TIME,ZONE,DTVALUES) C IFNO=IFNO+1 WRITE(EXT,50) 10000+IFNO WRITE(FXT,50) 10000+INT(V0) 50 FORMAT(I5) FILEN='driftdata'//FXT(2:5)//'.'//EXT(3:5) INQUIRE(FILE=FILEN,EXIST=CHECK) IF(CHECK) THEN PRINT 70, FILEN 70 FORMAT('Found existing file ',A30,' - Delete it.') OPEN(UNIT=1,FILE=FILEN,STATUS='OLD') CLOSE(UNIT=1,STATUS='DELETE') ENDIF OPEN(UNIT=1,FILE=FILEN,STATUS='NEW') WRITE(1,100) V0,A,D0,RAT,CURP,CURM,CURDP,CURDM 100 FORMAT(' Ion drift simulation data',5X,'V0=',F6.1,'V, A=', 1 F5.3,'mm,',' D0=',1P,D10.3,', rat=',0P,F12.5, 2 ', ','CURs=',1P,D12.5,', ',D12.5,', ',D12.5,', ',D12.5) WRITE(1,110) VAR(2),VAR(3),VAR(4),IV(3),NN, 1 IV(5),IV(7),IV(9),IV(10),IV(11) 110 FORMAT(' VAR2=',1P,E12.5,' VAR3=',E12.5,' VAR4=',E12.5, 1 ' ITER=',I10,' NN=',I4,' IV5=',I2,' IV7=',I2,' IV9=',I2, 2 ' IV10=',I2,' IV11=',I2) WRITE(1,120) VAR(7),VAR(8),VAR(9), 1 VAR(10),VAR(11),VAR(12),VAR(13), 2 DATE(1:4),DATE(5:6),DATE(7:8), 3 TIME(1:2),TIME(3:4),TIME(5:6) 120 FORMAT(' VAR7=',1P,E12.5,' VAR8=',E12.5,' VAR9=',E12.5, 1 ' VAR10=',E12.5,' VAR11=',E12.5, 2 ' VAR12=',E12.5,' VAR13=',E12.5,2X, 3 A4,'/',A2,'/',A2,1X,A2,':',A2,':',A2) WRITE(1,130) 130 FORMAT(' Comments:') WRITE(1,140) 140 FORMAT(1X,4X,'I',6X,'QP',9X,'QB',9X,'QM',9X,'QX',9X,'QG',9X, 1 'QQ',9X,'EZ',9X,'VV',9X, 2 'JP',9X,'JB',9X,'JM',9X,'JX',9X,'JQ') WRITE(1,190) JP(0),JB(0),JM(0),JX(0),JQ(0) 190 FORMAT(1X,' 0',88X,1P,5E11.3) DO I=1,NN WRITE(1,200) I,QP(I),QB(I),QM(I),QX(I),QG(I),QQ(I),EZ(I), 1 VV(I),JP(I),JB(I),JM(I),JX(I),JQ(I) 200 FORMAT(1X,I5,1P,13D11.3) ENDDO K=1 !Don't use SG(1) DO I=1,250 SGI(I)=0. DO J=1,40 K=K+1 IF(K.GT.10000) K=10000 SGI(I)=SGI(I)+SG(K) ENDDO ENDDO WRITE(1,300) VAR(6),VAR(3),SZNN 300 FORMAT(/,' Signal normalized to VAR6=',D13.5, 1 ' with time steps of VAR3=',D13.5,' SZNN=',D10.3) WRITE(1,320) SGI 320 FORMAT(25(2X,10F10.4,/)) WRITE(1,340) 340 FORMAT(/,' Signal as above but in 1/40th the time bins') WRITE(1,320) (1000.*SG(I),I=2,121) !Don't use SG(1) C Test C WRITE(1,360) MZN,NZN,IZN,EZNN,SZNN,EZN C 360 FORMAT(//,' MZN,NZN,IZN,EZNN,SZNN=',3I7,F10.4,D10.3,/, C 1 200(2X,10F10.4,/)) CLOSE(UNIT=1) C PRINT 400, FILEN 400 FORMAT(' Wrote data file ',A30) C RETURN END SUBROUTINE INSAVE C Read in data from disk file 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 /CHARGN/NNN,NDD,QPI(N),QBI(N),QMI(N),QXI(N),QGI(N),QQI(N) REAL*8 QPI,QBI,QMI,QXI,QGI,QQI COMMON /CHRGDV/QPZ(0:1,0:N),QBZ(0:1,0:N),QMZ(0:1,0:N), 1 QXZ(0:1,0:N),QGZ(0:1,0:N) REAL*8 QPZ,QBZ,QMZ,QXZ,QGZ COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /CRNTS/JP(0:N),JB(0:N),JM(0:N),JX(0:N),JG(0:N),JQ(0:N) REAL*8 JP,JB,JM,JX,JG,JQ COMMON /LPARS/V0,A,D0 REAL*8 V0,A,D0 COMMON /CRIT/DC,RC,JC,SCRPTR,RAT REAL*8 DC,RC,JC,SCRPTR,RAT COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) REAL*8 DIF DATA DIF/0.00001/ REAL*8 EZI,VVI,JPI,JBI,JMI,JXI,JGI,JQI CHARACTER*80 ASC LOGICAL LOG CHARACTER*30 FILEN CHARACTER*5 EXT,FXT COMMON /FILENO/IFNO C WRITE(EXT,50) 10000+IFNO WRITE(FXT,50) 10000+INT(V0) 50 FORMAT(I5) FILEN='driftdata'//FXT(2:5)//'.'//EXT(3:5) C Read in file previously written out OPEN(UNIT=1,FILE=FILEN,STATUS='OLD',ERR=400) READ(1,100) ASC 100 FORMAT(A80) C 1 2 3 4 5 6 7 READ(1,110) VAR2,VAR3,VAR4,NITER,NNN,IV5,IV7,IV9,IV10,IV11 110 FORMAT(6X,E12.5,6X,E12.5,6X,E12.5,6X,I10,4X,I4,5X,I2, 1 5X,I2,5X,I2,6X,I2,6X,I2) READ(1,120) VAR7,VAR8,VAR9,VAR10,VAR11,VAR12,VAR13 120 FORMAT(6X,E12.5,6X,E12.5,6X,E12.5,7X,E12.5,7X,E12.5, 1 7X,E12.5,7X,E12.5) READ(1,100) ASC READ(1,100) ASC IF(ABS(VAR(2)-VAR2).GT.DIF) PRINT 200 200 FORMAT(' VAR2 in INSAVE has changed') IF(ABS(VAR(3)-VAR3).GT.DIF) PRINT 210 210 FORMAT(' VAR3 in INSAVE has changed') IF(ABS(VAR(4)-VAR4).GT.DIF) PRINT 220 220 FORMAT(' VAR4 in INSAVE has changed') IF(ABS(VAR(7)-VAR7).GT.DIF) PRINT 230 230 FORMAT(' VAR7 in INSAVE has changed') IF(ABS(VAR(8)-VAR8).GT.DIF) PRINT 240 240 FORMAT(' VAR8 in INSAVE has changed') IF(ABS(VAR(9)-VAR9).GT.DIF) PRINT 250 250 FORMAT(' VAR9 in INSAVE has changed') IF(ABS(VAR(10)-VAR10).GT.DIF) PRINT 260 260 FORMAT(' VAR10 in INSAVE has changed') IF(ABS(VAR(11)-VAR11).GT.DIF) PRINT 261 261 FORMAT(' VAR11 in INSAVE has changed') IF(ABS(VAR(12)-VAR12).GT.DIF) PRINT 262 262 FORMAT(' VAR12 in INSAVE has changed') IF(ABS(VAR(13)-VAR13).GT.DIF) PRINT 263 263 FORMAT(' VAR13 in INSAVE has changed') C print 264, VAR(10),VAR10,VAR(10)-VAR10 264 FORMAT(' VAR(10) = ',E15.5,' VAR10 = ',E15.5, 1 ' DIFF = ',E15.5) LOG= (ABS(VAR(2)-VAR2).GT.DIF) 1 .OR.(ABS(VAR(3)-VAR3).GT.DIF) 2 .OR.(ABS(VAR(4)-VAR4).GT.DIF) 3 .OR.(ABS(VAR(7)-VAR7).GT.DIF) 4 .OR.(ABS(VAR(8)-VAR8).GT.DIF) 5 .OR.(ABS(VAR(9)-VAR9).GT.DIF) 6 .OR.(ABS(VAR(10)-VAR10).GT.DIF) 7 .OR.(ABS(VAR(11)-VAR11).GT.DIF) 8 .OR.(ABS(VAR(12)-VAR12).GT.DIF) 9 .OR.(ABS(VAR(13)-VAR13).GT.DIF) IF(LOG) THEN PRINT 280 280 FORMAT(' Stopping due to LOG = .TRUE.') CLOSE(UNIT=1) STOP ENDIF READ(1,290) J,JPI,JBI,JMI,JXI,JQI 290 FORMAT(1X,I5,88X,5E11.3) IV(3)=NITER DO I=1,NNN READ(1,300) J,QPI(I),QBI(I),QMI(I),QXI(I),QGI(I),QQI(I), 1 EZI,VVI,JPI,JBI,JMI,JXI,JQI 300 FORMAT(1X,I5,13E11.3) ENDDO PRINT 370, FILEN 370 FORMAT(' Successfully read in data file ',A30) CLOSE(UNIT=1) RETURN 400 PRINT 500, FILEN 500 FORMAT(' Error reading in data file ',A30) C RETURN END SUBROUTINE SAVEOTHER C Write out data to disk file 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 /CHRGDV/QPZ(0:1,0:N),QBZ(0:1,0:N),QMZ(0:1,0:N), 1 QXZ(0:1,0:N),QGZ(0:1,0:N) REAL*8 QPZ,QBZ,QMZ,QXZ,QGZ COMMON /JPHY/JDP(0:N),JDB(0:N),JDM(0:N),JDX(0:N), 1 JDG(0:N),JDT(0:N) REAL*8 JDP,JDB,JDM,JDX,JDG,JDT COMMON /FIELD/EZ(N),EDZ(0:N),DEZ(0:N),VV(N) REAL*8 EZ,EDZ,DEZ,VV COMMON /CRNTS/JP(0:N),JB(0:N),JM(0:N),JX(0:N),JG(0:N),JQ(0:N) REAL*8 JP,JB,JM,JX,JG,JQ COMMON /SLOPES/QPS(0:1,0:N),QBS(0:1,0:N),QMS(0:1,0:N), 1 QXS(0:1,0:N),QGS(0:1,0:N) REAL*8 QPS,QBS,QMS,QXS,QGS COMMON /DPHY/DQP(N),DQB(N),DQM(N),DQX(N),DQG(N),DQQ(N) REAL*8 DQP,DQB,DQM,DQX,DQG,DQQ COMMON /DCSC/DCP(N),DCB(N),DCM(N),DCX(N),DCG(N) REAL*8 DCP,DCB,DCM,DCX,DCG COMMON /LPARS/V0,A,D0 REAL*8 V0,A,D0 COMMON /OTHPMS/EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 REAL*8 EOEPS,DZ,DT,DTMDZ,DTODZ,RCOM,PROB2200 COMMON /CRIT/DC,RC,JC,SCRPTR,RAT REAL*8 DC,RC,JC,SCRPTR,RAT COMMON /SIG/JT,NT,SG(10000) REAL*8 SG REAL*8 SGI(250) COMMON /VARIBL/VAR(15),DVR(15),IV(15),NV,IL,LN(15) COMMON /CURNTS/CURP,CURM,CURDP,CURDM REAL*8 CURP,CURM,CURDP,CURDM COMMON /CURNT/MZN,NZN,IZN,IDM,EZNN,SZNN,EZN(2000) REAL*8 EZNN,SZNN,EZN COMMON /DIFFUS/DFP,DFB,DFM,DFX,DFG REAL*8 DFP,DFB,DFM,DFX,DFG CHARACTER*30 FILEN CHARACTER*5 EXT, FXT COMMON /FILENO/IFNO REAL*8 DSUM,DSMM REAL*8 VV1,VV2 BYTE FF DATA FF/"14/ REAL*8 DCPT,DCBT,DCMT,DCXT,DCGT,DQPT,DQBT,DQMT,DQXT,DQGT CHARACTER*8 DATE CHARACTER*9 TIME CHARACTER*5 ZONE INTEGER DTVALUES(8) LOGICAL CHECK C WRITE(EXT,50) 10000+IFNO WRITE(FXT,50) 10000+INT(V0) 50 FORMAT(I5) FILEN='driftother'//FXT(2:5)//'.'//EXT(3:5) INQUIRE(FILE=FILEN,EXIST=CHECK) IF(CHECK) THEN PRINT 70, FILEN 70 FORMAT('Found existing file ',A30,' - Delete it.') OPEN(UNIT=1,FILE=FILEN,STATUS='OLD') CLOSE(UNIT=1,STATUS='DELETE') ENDIF OPEN(UNIT=1,FILE=FILEN,STATUS='NEW') DCPT=0.D0 DCBT=0.D0 DCMT=0.D0 DCXT=0.D0 DCGT=0.D0 DQPT=0.D0 DQBT=0.D0 DQMT=0.D0 DQXT=0.D0 DQGT=0.D0 CALL DATE_AND_TIME(DATE,TIME,ZONE,DTVALUES) WRITE(1,120) DZ,DT,DTMDZ,DTODZ,RAT,IV(3), 1 DATE(1:4),DATE(5:6),DATE(7:8), 2 TIME(1:2),TIME(3:4),TIME(5:6) 120 FORMAT(1X,'DZ,DT,DTMDZ,DTODZ,RAT=',1P,5D13.5,1X,'ITER=',I10, 1 1X,A4,'/',A2,'/',A2,1X,A2,':',A2,':',A2) WRITE(1,140) 140 FORMAT(1X,4X,'I',6X,'DCP',8X,'DCB',8X,'DCM',8X,'DCX', 1 8X,'DCG',9X,'EZ',9X,'VP',8X,'VPD', 2 8X,'DQP',8X,'DQB',8X,'DQM',8X,'DQX',8X,'DQQ') DO I=1,NN DCPT=DCPT+DCP(I) DCBT=DCBT+DCB(I) DCMT=DCMT+DCM(I) DCXT=DCXT+DCX(I) DCGT=DCGT+DCG(I) DQPT=DQPT+DQP(I) DQBT=DQBT+DQB(I) DQMT=DQMT+DQM(I) DQXT=DQXT+DQX(I) DQGT=DQGT+DQG(I) WRITE(1,200) I,DCP(I),DCB(I),DCM(I),DCX(I),DCG(I),EZ(I), 1 VP(I),VPD(I),DQP(I),DQB(I),DQM(I),DQX(I),DQQ(I) 200 FORMAT(1X,I5,1P,13D11.3) ENDDO WRITE(1,205) 205 FORMAT(/) WRITE(1,210) DCPT,DCBT,DCMT,DCXT,DCGT,DQPT,DQBT,DQMT,DQXT 210 FORMAT(6X,1P,5D11.3,33X,4D11.3) DCPT=DCPT/DT DCBT=DCBT/DT DCMT=DCMT/DT DCXT=DCXT/DT DCGT=DCGT/DT DQPT=DQPT/DT DQBT=DQBT/DT DQMT=DQMT/DT DQXT=DQXT/DT DQGT=DQGT/DT WRITE(1,210) DCPT,DCBT,DCMT,DCXT,DCGT,DQPT,DQBT,DQMT,DQXT WRITE(1,212) DFP,DFB,DFM,DFX,DFG 212 FORMAT(/,' Diffusion constants:',1P,5E11.3) WRITE(1,220) FF 220 FORMAT(A1) WRITE(1,240) 240 FORMAT(1X,4X,'I',4X,'QPS(0)',5X,'QPS(1)',5X,'QPZ(0)', 1 5X,'QPZ(1)',6X,'EDZ',9X,'VP',8X,'VPD', 2 8X,'JDP',8X,'JDB',8X,'JDM',8X,'JDX',8X,'JDT') DO I=0,NN WRITE(1,300) I,QPS(0,I),QPS(1,I),QPZ(0,I),QPZ(1,I),EDZ(I), 1 VP(I),VPD(I),JDP(I),JDB(I),JDM(I),JDX(I),JDT(I) 300 FORMAT(1X,I5,1P,12D11.3) ENDDO C Calculate potential from E-field VV1=0. VV2=EDZ(0)/2. DO I=1,NN VV1=VV1+EZ(I) VV2=VV2+EDZ(I) ENDDO VV2=VV2-EDZ(NN)/2. VV1=VV1*DZ VV2=VV2*DZ WRITE(1,330) VV1,VV2 330 FORMAT(//,' Potentials ', 1P,2D20.10) CLOSE(UNIT=1) C PRINT 400, FILEN 400 FORMAT(' Wrote data file ',A30) C RETURN END