========== beginning of macroem.f ================================= C*********************************************************************** C* * C* ########### HIDENEK1.FORT ########### * C* * C* LARGE TIME-AND-SPACE SCALES, IMPLICIT PARTICLE CODE. * C* * C*********************************************************************** C* * C* >> THIS VERSION << * C* HOMOGENEOUS PLASMA, 1-D FOURIER SPACE SOLVER. * C* TWO SPECIES ASSUMED: (1) IONS, (2) ELECTRONS. * C* * C* * C* AUTHOR: MOTOHIKO TANAKA * C* NATIONAL INSTITUTE FOR FUSION SCIENCE * C* NAGOYA 464-01, JAPAN * C* * C* REFERENCES: J. COMPUT. PHYSICS, VOL 79, 209 (1988). * C J. COMPUT. PHYSICS, IN PRESS (1992). * C* * C******************************************** VERSION: 11.23.91 ******** C * C TIME-SPACE RESTRICTIONS TO THIS MACROSCALE PARTICLE * C CODE. EXPERIENCE RECOMMENDS THE FOLLWING VALUES: * C (1) WCI*DT < 0.2 (FULMOV, ESCORR) * C (2) (1.-GIMPL)*KZ*V*DT < 0.1 (ESCORR) * C (*) WPE*DT >> 1, KZ*DEBYE >> 1. * C * C----------------------------------------------------------------------- C* * C* * << E/M FIELD SOLVER >> * C* FIELD-PARTICLE COUPLED EQUATIONS ARE SOLVED. * C* (ET & EL PARTS ARE NOT SEPARATED) * C* * C* * << PARTICLE MOVER >> * C* ELECTRONS ... PARALLEL MOTION & E*B DRIFT. * C* IONS ........ FULLY KINETIC. * C* * C* * CHOOSE IMPLICIT PARAMETERS : 0.5 < AIMPL, GIMPL < 1. * C* TYPICALLY, AIMPL= GIMPL= 0.6. * C* * C*********************************************************************** C* * C* -0.0----------------------ZMAX-- * C* I I * C* I X X X X X I X FIELD AND PARTICLE GRIDS. * C* I I * C* ---------------------------------> Z * C* 1 2 ... ... MZ 1 * C* * C*********************************************************************** C* * C* FILE 10..... LOAD * C* FILE 12..... DUMP * C* FILE 13..... PDATA * C* FILE 15..... FDATA * C* * C*********************************************************************** C PARAMETER (MZ=128,NPRT=12800) C*** CHARACTER*8 IDATE,TTL*64 COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) COMMON/KFIELD/ EKX(MZ),EKY(MZ),EKZ(MZ),BKX(MZ),BKY(MZ),BKZ(MZ) C NAMELIST/NAMEL1/ KSTART,NTIMES,CPTOT,EPSLN1, * ZMAX,AZ,DT,AIMPL,GIMPL,QSPEC,WSPEC,NSPEC,NSP, * VDR,VBEAM,VETH,TETI,TANIS,WCEWPE,THB NAMELIST/NAMEL2/ VLIMA,VLIMB,BMIN,EMIN,NPWRT,NFWRT,NPLOT,NHA, * NHIST,KDUMP NAMELIST/NAMEL3/ NTITLE,TTL C C C*********************************************************************** C* READ-IN CONTROL PARAMETERS. * C*********************************************************************** C CALL CLOCKS (CPUI) C WRITE(6,*) '### MACRO1H.FORT (FULL-IMPLICIT) ###' READ(5,NAMEL1) C** IF(KSTART.NE.0) GO TO 21 C** 25 READ(5,NAMEL2) READ(5,NAMEL3) C WRITE(6,*) '*INPUT DATA BY &NAMEL1, &NAMEL2 ...' WRITE(6,NAMEL1) WRITE(6,NAMEL2) WRITE(6,NAMEL3) WRITE(6,'(1H ,//,3X,A64)') TTL C C CALL TIME(ITIME) CALL DATE(IDATE) TS= ITIME/1000.0 IH= IFIX(TS/3600.0) TM= (TS-IH*3600.0)/60.0 IM= IFIX(TM) IS= (TM-FLOAT(IM))*60.0 WRITE(6,603) IH,IM,IS,IDATE 603 FORMAT(1H0,' **TIME .....',I2,':',I2,':',I2,' DATE .....',A8,/) C C C*********************************************************************** C* INITIAL LOADING OF PARTICLES AND FIELDS. * C*********************************************************************** C CALL INIT (KSTART) WRITE(6,*) '**INITIALIZATION COMPLETED**' C C C*********************************************************************** C* MAIN PART. * C*********************************************************************** C 100 IT= IT+1 T = T +WCEWPE*DT/WSPEC(1) CALL TRANS C CALL CLOCKS (CPU1) CPU1= CPU1/60.0 IF((CPU1.LT.CPTOT).AND.(IT.LE.NTIMES)) GO TO 100 C C C*********************************************************************** C* FINAL PROCEDURE TO CLOSE THE JOB. * C*********************************************************************** C ISTOP= 1 C CALL CLOCKS (CPUF) TCP= CPUF-CPUI WRITE(6,17) TCP,IT 17 FORMAT(1H0,'**TIME USED (TOTAL)',F9.2,' SECS, IT=',I5,'*****',/) C C----------------------------------------------------------------------- C* RESTART FILE. C----------------------------------------------------------------------- C NB1= 10 NB2= 62 +256*23 NB3= 18*MZ CALL RESTRT (IT,NB1,ZMAX,NB2,EX,NB3,2) CALL EXIT C C C----------------------------------------------------------------------- C* RESTART PROCEDURE. C----------------------------------------------------------------------- C 21 NB1= 10 NB2= 62 +256*23 NB3= 18*MZ CALL RESTRT (IT,NB1,ZMAX,NB2,EX,NB3,1) C GO TO 25 END C C C C----------------------------------------------------------------------- SUBROUTINE TRANS C----------------------------------------------------------------------- PARAMETER (MZ=128,NPRT=12800) C***** COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) COMMON/KFIELD/ EKX(MZ),EKY(MZ),EKZ(MZ),BKX(MZ),BKY(MZ),BKZ(MZ) COMMON/SRIMPL/ QIX(MZ),QIY(MZ),QIZ(MZ),QJX(MZ),QJY(MZ),QJZ(MZ), * QEX(MZ),QEY(MZ),QEZ(MZ), QI(MZ), QB(MZ), QE(MZ), * QJH(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) COMMON/EMITER/ ASE,ASB,ASEL,WE,WB,WL,ITERM,ITERS COMMON/TIMING/ CTIME C* REAL*4 MUE COMMON/PARTI/ ZI(NPRT),VXI(NPRT),VYI(NPRT),VZI(NPRT) COMMON/PARTE/ ZE(NPRT),MUE(NPRT),VPE(NPRT),VHE(NPRT) C C*** CALL CLOCKS (CT1) C*** IF(MOD(IT,NPLOT).EQ.0) THEN WRITE(6,776) IT 776 FORMAT(1H0,'**** IT=',I4,' *****') END IF C C IF(IT.EQ.0) THEN DTSAV= DT GDTSAV=GDT HDTSAV=HDT DT=0. GDT=0. HDT=0. C IPC= -1 CALL FULMOV (ZI,VXI,VYI,VZI,QIX,QIY,QIZ,QJX,QJY,QJZ,QI, * QSPEC(1),WSPEC(1),NSPEC(1),IPC) CALL DRMOVE (ZE,MUE,VPE,VHE,QEX,QEY,QEZ,QJH,QE, * QSPEC(2),WSPEC(2),NSPEC(2),IPC) C CALL EMFLD0 CALL INITVE (ZE,MUE,VPE,VHE,NSPEC(2)) C DT= DTSAV HDT=HDTSAV GDT=GDTSAV GO TO 700 END IF C C C*********************************************************************** C* 1. ACCUMULATE SOURCE TERMS AT T= T(N) (IPC=1). * C* NEW TIME CYCLE STARTS HERE. * C*********************************************************************** C C IPC= 1 CALL FULMOV (ZI,VXI,VYI,VZI,QIX,QIY,QIZ,QJX,QJY,QJZ,QI, * QSPEC(1),WSPEC(1),NSPEC(1),IPC) CALL DRMOVE (ZE,MUE,VPE,VHE,QEX,QEY,QEZ,QJH,QE, * QSPEC(2),WSPEC(2),NSPEC(2),IPC) C C C*********************************************************************** C* 2. FIELD SOLVE BY FULL IMPLICIT METHOD. * C*********************************************************************** C* EPSLN1 : CRITERIA OF THE ITERATION CONVERGENCE (EMFILD, ESCORR). C CALL EMFILD C C C----------------------------------------------------------------------- C* CORRECTION TO THE LONGITUDINAL PART IS MADE SEPARATELY OUTSIDE C* THE FIELD ITERATION LOOP (A GOOD APPROXIMATION). C----------------------------------------------------------------------- C CALL FULMV2 (ZI,VXI,VYI,VZI,QI,QSPEC(1),WSPEC(1),NSPEC(1)) CALL DRMOV2 (ZE,MUE,VPE,VHE,QE,QSPEC(2),WSPEC(2),NSPEC(2)) C CALL ESCORR C C C*********************************************************************** C* 3. UPDATE PARTICLES : X(N) TO X(N+1) (IPC=0). * C*********************************************************************** C IPC= 0 CALL FULMOV (ZI,VXI,VYI,VZI,QIX,QIY,QIZ,QJX,QJY,QJZ,QI, * QSPEC(1),WSPEC(1),NSPEC(1),IPC) CALL DRMOVE (ZE,MUE,VPE,VHE,QEX,QEY,QEZ,QJH,QE, * QSPEC(2),WSPEC(2),NSPEC(2),IPC) C C C----------------------------------------------------------------------- C REWRITE : EX0 <--- EX. C----------------------------------------------------------------------- C DO 300 K= 1,MZ EX0(K)= EX(K) EY0(K)= EY(K) EZ0(K)= EZ(K) BX0(K)= BX(K) BY0(K)= BY(K) BZ0(K)= BZ(K) 300 CONTINUE C C 700 CALL DIAG1 (ZI,VXI,VYI,VZI,NSPEC(1),1) CALL DIAG1 (ZE,MUE,VPE,VHE,NSPEC(2),2) C CALL DIAG2 C*** CALL CLOCKS (CT2) CTIME= CT2 -CT1 C*** C RETURN END C C C----------------------------------------------------------------------- BLOCK DATA C----------------------------------------------------------------------- COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C DATA DT/20./,PI/3.1415927/,ZMAX/120./,AZ/5./, * NSP/2/,NSPEC/12800,12800,2*0/, * WSPEC/50.0, 1.0, 2*0./,QSPEC/1.0, -1.0, 2*0./,VBEAM/4*0./, * VETH/0.2/,TETI/1.0/,TANIS/10./,WCEWPE/0.7/,THB/0.0/, * AIMPL/0.60/,GIMPL/0.60/,EPSLN1/1.0E-3/, * VLIMA/0.6/,VLIMB/0.2/,BMIN/1.0E-7/,EMIN/1.0E-7/, * NHA/50/,NPLOT/1000/,NHIST/1000/,NPWRT/9999/,NFWRT/9999/ END C C C C----------------------------------------------------------------------- SUBROUTINE FULMOV (Z,VX,VY,VZ,QIX,QIY,QIZ,QJX,QJY,QJZ,QI, * QMULT,WMULT,NPR,IPC) C----------------------------------------------------------------------- C*********************************************************************** C* (ANALYTICAL) PARTICLE MOVER WITH V*B ROTATION. * C*********************************************************************** C PARAMETER (MZ=128,NPRT=12800) DIMENSION Z(NPRT),VX(NPRT),VY(NPRT),VZ(NPRT) DIMENSION QIX(MZ),QIY(MZ),QIZ(MZ),QJX(MZ),QJY(MZ),QJZ(MZ),QI(MZ) C INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) C COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C DIMENSION EXA(MZ),EYA(MZ),EZA(MZ),BXA(MZ),BYA(MZ),BZA(MZ) DIMENSION RZG(NPRT),VXJ(NPRT),VYJ(NPRT),VZJ(NPRT) C C C----------------------------------------------------------------------- H = DT*QMULT/WMULT HE= 0.5*H C IF(IPC.NE.0) GO TO 200 C----------------------------------------------------------------------- C*********************************************************************** C* 1. MOVE PARTICLES: * C*********************************************************************** C C* DEFINE E,B OF TIME LEVEL (N+AIMPL). C DO 1 K= 1,MZ EXA(K)= AIMPL*EX(K) +(1.-AIMPL)*EX0(K) EYA(K)= AIMPL*EY(K) +(1.-AIMPL)*EY0(K) EZA(K)= AIMPL*EZ(K) +(1.-AIMPL)*EZ0(K) BXA(K)= AIMPL*BX(K) +(1.-AIMPL)*BX0(K) BYA(K)= AIMPL*BY(K) +(1.-AIMPL)*BY0(K) BZA(K)= AIMPL*BZ(K) +(1.-AIMPL)*BZ0(K) 1 CONTINUE C C DO 100 L= 1,NPR RZ= HZI*(Z(L) +GDT*VZ(L)) +1.0001 KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C EXI= EXA(KR)*FZR +EXA(K)*FZC +EXA(KL)*FZL EYI= EYA(KR)*FZR +EYA(K)*FZC +EYA(KL)*FZL EZI= EZA(KR)*FZR +EZA(K)*FZC +EZA(KL)*FZL C C RZ= HZI*(Z(L) +HDT*VZ(L)) +1.0001 KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C BXI= BXA(KR)*FZR +BXA(K)*FZC +BXA(KL)*FZL BYI= BYA(KR)*FZR +BYA(K)*FZC +BYA(KL)*FZL BZI= BZA(KR)*FZR +BZA(K)*FZC +BZA(KL)*FZL BSQI= BXI*BXI +BYI*BYI +BZI*BZI C VLX= ( VX(L) +HE*(EXI +VY(L)*BZI -VZ(L)*BYI * +HE*(EYI*BZI -EZI*BYI)) )/(1.+HE**2*BSQI) VLY= ( VY(L) +HE*(EYI +VZ(L)*BXI -VX(L)*BZI * +HE*(EZI*BXI -EXI*BZI)) )/(1.+HE**2*BSQI) VLZ= ( VZ(L) +HE*(EZI +VX(L)*BYI -VY(L)*BXI * +HE*(EXI*BYI -EYI*BXI)) )/(1.+HE**2*BSQI) C C* NEW POSITION AND VELOCITY. C Z(L)= Z(L) +DT*(VZ(L) +0.5*H*(EZI +VLX*BYI -VLY*BXI)) C VX(L)= VX(L) +H*(EXI +VLY*BZI -VLZ*BYI) VY(L)= VY(L) +H*(EYI +VLZ*BXI -VLX*BZI) VZ(L)= VZ(L) +H*(EZI +VLX*BYI -VLY*BXI) 100 CONTINUE C CALL BUNDRY (Z,VX,VY,VZ,NPR,KSP) RETURN C C C*********************************************************************** C* 2. ACCUMULATE MOMENTS FOR THE E.M. FIELD (SRIMP1 - 2). * C*********************************************************************** C 200 DO 300 L= 1,NPR C RZ= HZI*(Z(L) +GDT*VZ(L)) +1.0001 KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C EXS= EX0(KR)*FZR +EX0(K)*FZC +EX0(KL)*FZL EYS= EY0(KR)*FZR +EY0(K)*FZC +EY0(KL)*FZL EZS= EZ0(KR)*FZR +EZ0(K)*FZC +EZ0(KL)*FZL C C RZ= HZI*(Z(L) +HDT*VZ(L)) +1.0001 KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C BXS= BX0(KR)*FZR +BX0(K)*FZC +BX0(KL)*FZL BYS= BY0(KR)*FZR +BY0(K)*FZC +BY0(KL)*FZL BZS= BZ0(KR)*FZR +BZ0(K)*FZC +BZ0(KL)*FZL C VXJ(L)= VX(L) +GIMPL*H*(EXS +VY(L)*BZS -VZ(L)*BYS) VYJ(L)= VY(L) +GIMPL*H*(EYS +VZ(L)*BXS -VX(L)*BZS) VZJ(L)= VZ(L) +GIMPL*H*(EZS +VX(L)*BYS -VY(L)*BXS) C RZG(L)= (Z(L) +GDT*(VZ(L) +0.5*H*(EZS +VX(L)*BYS -VY(L)*BXS)))/HZ * +1.0001 300 CONTINUE C C CALL BUNDR2 (RZG,VXJ,VYJ,VZJ,NPR,KSP) C CALL SRIMP2 (QI,RZG,QMULT,NPR) CALL SRIMP1 (QIX,QIY,QIZ,RZG,VXJ,VYJ,VZJ,QMULT,NPR) CALL SRIMP1 (QJX,QJY,QJZ,RZG, VX, VY, VZ,QMULT,NPR) C C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE DRMOVE (Z,MUE,VY,VH,QEX,QEY,QEZ,QJH,QE, * QMULT,WMULT,NPR,IPC) C----------------------------------------------------------------------- C********************************************************* 11/30/87 **** C* E*B DRIFT MOTION(PERPENDICULAR), -MUE*GRAD.B FORCE (PARALLEL) * C* MUE(L) HOLDS THE MAGNETIC MOMENT, VH(L) THE PARALLEL VELOCITY. * C* (NORMAL VX-VZ ARE CONVERTED INTO (MUE,VPERP,VPARA) BY SUB.INIT) * C*********************************************************************** C PARAMETER (MZ=128,NPRT=12800) REAL*4 MUE DIMENSION Z(NPRT),MUE(NPRT),VY(NPRT),VH(NPRT) DIMENSION QEX(MZ),QEY(MZ),QEZ(MZ),QJH(MZ),QE(MZ) C INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C DIMENSION EXA(MZ),EYA(MZ),EZA(MZ),BXA(MZ),BYA(MZ),BZA(MZ), * BXH(MZ),BYH(MZ),BZH(MZ) DIMENSION RZG(NPRT),VXJ(NPRT),VYJ(NPRT),VZJ(NPRT) C C C----------------------------------------------------------------------- H= DT*QMULT/WMULT IF(IPC.NE.0) GO TO 200 C----------------------------------------------------------------------- C*********************************************************************** C* 1. MOVE PARTICLES: * C*********************************************************************** C C* DEFINE E,B OF TIME LEVELS (N+AIMPL), AND (N+GIMPL). C DO 1 K= 1,MZ EXA(K)= AIMPL*EX(K) +(1.-AIMPL)*EX0(K) EYA(K)= AIMPL*EY(K) +(1.-AIMPL)*EY0(K) EZA(K)= AIMPL*EZ(K) +(1.-AIMPL)*EZ0(K) BXA(K)= AIMPL*BX(K) +(1.-AIMPL)*BX0(K) BYA(K)= AIMPL*BY(K) +(1.-AIMPL)*BY0(K) BZA(K)= AIMPL*BZ(K) +(1.-AIMPL)*BZ0(K) C BXH(K)= 0.5*(BX(K) +BX0(K)) BYH(K)= 0.5*(BY(K) +BY0(K)) BZH(K)= 0.5*(BZ(K) +BZ0(K)) 1 CONTINUE C C DO 100 L= 1,NPR K= PZC(INT(HZI*Z(L) +1.0001)) BSQRT= SQRT(BX0(K)**2 +BY0(K)**2 +BZ0(K)**2) C RZ= HZI*(Z(L) +GDT*VH(L)*BZ0(K)/BSQRT) +1.0001 KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C C AEX= EXA(KR)*FZR +EXA(K)*FZC +EXA(KL)*FZL AEY= EYA(KR)*FZR +EYA(K)*FZC +EYA(KL)*FZL AEZ= EZA(KR)*FZR +EZA(K)*FZC +EZA(KL)*FZL ABX= BXA(KR)*FZR +BXA(K)*FZC +BXA(KL)*FZL ABY= BYA(KR)*FZR +BYA(K)*FZC +BYA(KL)*FZL ABZ= BZA(KR)*FZR +BZA(K)*FZC +BZA(KL)*FZL ABSQ= ABX**2 +ABY**2 +ABZ**2 C HBX= BXH(KR)*FZR +BXH(K)*FZC +BXH(KL)*FZL HBY= BYH(KR)*FZR +BYH(K)*FZC +BYH(KL)*FZL HBZ= BZH(KR)*FZR +BZH(K)*FZC +BZH(KL)*FZL HBSQ= HBX**2 +HBY**2 +HBZ**2 C EPARA= (AEX*ABX +AEY*ABY +AEZ*ABZ)/SQRT(ABX**2 +ABY**2 +ABZ**2) VHH= VH(L) +0.5*H*EPARA VH(L)= VH(L) +H*EPARA C Z(L)= Z(L) +DT*(VHH*HBZ/SQRT(HBSQ) +(AEX*ABY -AEY*ABX)/ABSQ) 100 CONTINUE C CALL BUNDRY (Z,MUE,VY,VH,NPR,KSP) RETURN C C C*********************************************************************** C* 2. ACCUMULATE THE MOMENTS. * C*********************************************************************** C 200 DO 300 L= 1,NPR C K= PZC(INT(HZI*Z(L) +1.0001)) BSQRT= SQRT(BX0(K)**2 +BY0(K)**2 +BZ0(K)**2) RZ= HZI*(Z(L) +GDT*VH(L)*BZ0(K)/BSQRT) +1.0001 C KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C C EXS= EX0(KR)*FZR +EX0(K)*FZC +EX0(KL)*FZL EYS= EY0(KR)*FZR +EY0(K)*FZC +EY0(KL)*FZL EZS= EZ0(KR)*FZR +EZ0(K)*FZC +EZ0(KL)*FZL BXS= BX0(KR)*FZR +BX0(K)*FZC +BX0(KL)*FZL BYS= BY0(KR)*FZR +BY0(K)*FZC +BY0(KL)*FZL BZS= BZ0(KR)*FZR +BZ0(K)*FZC +BZ0(KL)*FZL BSSQ= BXS**2 +BYS**2 +BZS**2 BS0= SQRT(BSSQ) C EPARS= (EXS*BXS +EYS*BYS +EZS*BZS)/BS0 VHH= VH(L) + 0.5*H*EPARS VHG= VH(L) +GIMPL*H*EPARS C VXS= (EYS*BZS-EZS*BYS)/BSSQ VYS= (EZS*BXS-EXS*BZS)/BSSQ VZS= (EXS*BYS-EYS*BXS)/BSSQ C VXJ(L)= VHG*BXS/BS0 +VXS VYJ(L)= VHG*BYS/BS0 +VYS VZJ(L)= VHG*BZS/BS0 +VZS C RZG(L)= HZI*(Z(L) +GDT*(VHH*BZS/BS0 +VZS)) +1.0001 300 CONTINUE C C----------------------------------------------------------------------- C* VPARA MUST BE DEFINED SEPARATELY AT IT= 0. C* (NORMAL VX-VZ ARE REQUIRED.) C----------------------------------------------------------------------- C IF((IT.EQ.0).AND.(IPC.EQ.-1)) THEN DO 400 L= 1,NPR VXJ(L)= MUE(L) VYJ(L)= VY(L) VZJ(L)= VH(L) 400 CONTINUE END IF C C CALL BUNDR2 (RZG,VXJ,VYJ,VZJ,NPR,KSP) C CALL SRIMP2 ( QE,RZG, QMULT,NPR) CALL SRIMP3 (QJH,RZG,VH,QMULT,NPR) CALL SRIMP1 (QEX,QEY,QEZ,RZG,VXJ,VYJ,VZJ,QMULT,NPR) C C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE FULMV2 (Z,VX,VY,VZ,QI,QMULT,WMULT,NPR) C----------------------------------------------------------------------- C PARAMETER (MZ=128,NPRT=12800) DIMENSION Z(NPRT),VX(NPRT),VY(NPRT),VZ(NPRT) DIMENSION QI(MZ) C INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) C COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C DIMENSION EXA(MZ),EYA(MZ),EZA(MZ),BXA(MZ),BYA(MZ),BZA(MZ) DIMENSION RZL(NPRT),VXL(NPRT),VYL(NPRT),VZL(NPRT) C C C*********************************************************************** C* 1. MOVE PARTICLES: * C*********************************************************************** C C* DEFINE E,B OF TIME LEVEL (N+AIMPL). C DO 1 K= 1,MZ EXA(K)= AIMPL*EX(K) +(1.-AIMPL)*EX0(K) EYA(K)= AIMPL*EY(K) +(1.-AIMPL)*EY0(K) EZA(K)= AIMPL*EZ(K) +(1.-AIMPL)*EZ0(K) BXA(K)= AIMPL*BX(K) +(1.-AIMPL)*BX0(K) BYA(K)= AIMPL*BY(K) +(1.-AIMPL)*BY0(K) BZA(K)= AIMPL*BZ(K) +(1.-AIMPL)*BZ0(K) 1 CONTINUE C H = DT*QMULT/WMULT HE= 0.5*H C C DO 100 L= 1,NPR RZ= HZI*(Z(L) +GDT*VZ(L)) +1.0001 KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C EXI= EXA(KR)*FZR +EXA(K)*FZC +EXA(KL)*FZL EYI= EYA(KR)*FZR +EYA(K)*FZC +EYA(KL)*FZL EZI= EZA(KR)*FZR +EZA(K)*FZC +EZA(KL)*FZL C C RZ= HZI*(Z(L) +HDT*VZ(L)) +1.0001 KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C BXI= BXA(KR)*FZR +BXA(K)*FZC +BXA(KL)*FZL BYI= BYA(KR)*FZR +BYA(K)*FZC +BYA(KL)*FZL BZI= BZA(KR)*FZR +BZA(K)*FZC +BZA(KL)*FZL BSQI= BXI**2 +BYI**2 +BZI**2 C VLX= ( VX(L) +HE*(EXI +VY(L)*BZI -VZ(L)*BYI * +HE*(EYI*BZI -EZI*BYI)) )/(1.+HE**2*BSQI) VLY= ( VY(L) +HE*(EYI +VZ(L)*BXI -VX(L)*BZI * +HE*(EZI*BXI -EXI*BZI)) )/(1.+HE**2*BSQI) VLZ= ( VZ(L) +HE*(EZI +VX(L)*BYI -VY(L)*BXI * +HE*(EXI*BYI -EYI*BXI)) )/(1.+HE**2*BSQI) C RZL(L)= (Z(L) +DT*(VZ(L) +0.5*H*(EZI +VLX*BYI -VLY*BXI)))/HZ * +1.0001 100 CONTINUE C C CALL BUNDR2 (RZL,VXL,VYL,VZL,NPR,KSP) CALL SRIMP2 ( QI,RZL,QMULT,NPR) C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE DRMOV2 (Z,MUE,VY,VH,QE,QMULT,WMULT,NPR) C----------------------------------------------------------------------- C PARAMETER (MZ=128,NPRT=12800) REAL*4 MUE DIMENSION Z(NPRT),MUE(NPRT),VY(NPRT),VH(NPRT) DIMENSION QE(MZ) C INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) C COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C DIMENSION EXA(MZ),EYA(MZ),EZA(MZ),BXA(MZ),BYA(MZ),BZA(MZ), * BXH(MZ),BYH(MZ),BZH(MZ) DIMENSION RZL(NPRT),VXL(NPRT),VYL(NPRT),VZL(NPRT) C C*********************************************************************** C* 1. MOVE PARTICLES: * C*********************************************************************** C C* DEFINE E,B OF TIME LEVELS (N+AIMPL), AND (N+GIMPL). C DO 1 K= 1,MZ EXA(K)= AIMPL*EX(K) +(1.-AIMPL)*EX0(K) EYA(K)= AIMPL*EY(K) +(1.-AIMPL)*EY0(K) EZA(K)= AIMPL*EZ(K) +(1.-AIMPL)*EZ0(K) BXA(K)= AIMPL*BX(K) +(1.-AIMPL)*BX0(K) BYA(K)= AIMPL*BY(K) +(1.-AIMPL)*BY0(K) BZA(K)= AIMPL*BZ(K) +(1.-AIMPL)*BZ0(K) C BXH(K)= 0.5*(BX(K) +BX0(K)) BYH(K)= 0.5*(BY(K) +BY0(K)) BZH(K)= 0.5*(BZ(K) +BZ0(K)) 1 CONTINUE C C H= DT*QMULT/WMULT C DO 100 L= 1,NPR K= PZC(INT(HZI*Z(L) +1.0001)) BSQRT= SQRT(BX0(K)**2 +BY0(K)**2 +BZ0(K)**2) RZ= HZI*(Z(L) +GDT*VH(L)*BZ0(K)/BSQRT) +1.0001 C KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C AEX= EXA(KR)*FZR +EXA(K)*FZC +EXA(KL)*FZL AEY= EYA(KR)*FZR +EYA(K)*FZC +EYA(KL)*FZL AEZ= EZA(KR)*FZR +EZA(K)*FZC +EZA(KL)*FZL ABX= BXA(KR)*FZR +BXA(K)*FZC +BXA(KL)*FZL ABY= BYA(KR)*FZR +BYA(K)*FZC +BYA(KL)*FZL ABZ= BZA(KR)*FZR +BZA(K)*FZC +BZA(KL)*FZL C HBX= BXH(KR)*FZR +BXH(K)*FZC +BXH(KL)*FZL HBY= BYH(KR)*FZR +BYH(K)*FZC +BYH(KL)*FZL HBZ= BZH(KR)*FZR +BZH(K)*FZC +BZH(KL)*FZL ABSQ= ABX**2 +ABY**2 +ABZ**2 HBSQ= HBX**2 +HBY**2 +HBZ**2 C EPARA= (AEX*ABX +AEY*ABY +AEZ*ABZ)/SQRT(ABSQ) VHH= VH(L) +0.5*H*EPARA C RZL(L)= (Z(L) +DT*(VHH*HBZ/SQRT(HBSQ) +(AEX*ABY -AEY*ABX)/ABSQ)) * /HZ +1.0001 100 CONTINUE C C CALL BUNDR2 (RZL,VXL,VYL,VZL,NPR,KSP) CALL SRIMP2 ( QE,RZL,QMULT,NPR) C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE INITVE (Z,VX,VY,VZ,NPR) C----------------------------------------------------------------------- C*********************************************************************** C* PUT MAGNETIC MOMENT IN VX(L), PARALLEL VELOCITY IN VZ (WITH SIGN) * C*********************************************************************** C PARAMETER (MZ=128,NPRT=12800) DIMENSION Z(NPRT),VX(NPRT),VY(NPRT),VZ(NPRT) C INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C C DO 100 L= 1,NPR RZ= HZI*Z(L) +1.0001 C KK= RZ K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZ -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C BXA= BX(KR)*FZR +BX(K)*FZC +BX(KL)*FZL BYA= BY(KR)*FZR +BY(K)*FZC +BY(KL)*FZL BZA= BZ(KR)*FZR +BZ(K)*FZC +BZ(KL)*FZL BSQRT= SQRT(BXA*BXA+BYA*BYA+BZA*BZA) C C* # COMPUTE PARALLEL VELOCITY ### C VX(L)= 0. VY(L)= 0. VZ(L)= (VX(L)*BXA +VY(L)*BYA +VZ(L)*BZA)/BSQRT 100 CONTINUE C C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE BUNDRY (Z,VX,VY,VZ,NPR,KSP) C----------------------------------------------------------------------- PARAMETER (MZ=128,NPRT=12800) C C* ## BOUNDARY CONDITIONS ### C DIMENSION Z(NPRT),VX(NPRT),VY(NPRT),VZ(NPRT) COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C C DO 100 L= 1,NPR IF(Z(L).GE.ZMAX) THEN Z(L)= Z(L)-ZMAXE C ELSE IF(Z(L).LE.0.) THEN Z(L)= Z(L)+ZMAXE END IF 100 CONTINUE C C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE BUNDR2 (RZM,VX,VY,VZ,NPR,KSP) C----------------------------------------------------------------------- C* ## BOUNDARY CONDITIONS ### C PARAMETER (MZ=128,NPRT=12800) DIMENSION RZM(NPRT),VX(NPRT),VY(NPRT),VZ(NPRT) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C C AZMAX= MZ +1.0001 AZMAXE= 0.99999*MZ C DO 100 L= 1,NPR C *VOCL STMT, IF(10) C**** IF(RZM(L).GE.AZMAX) THEN RZM(L)= RZM(L) -AZMAXE ELSE IF(RZM(L).LE.1.0001) THEN RZM(L)= RZM(L) +AZMAXE END IF 100 CONTINUE C C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE SRIMP1 (QJX,QJY,QJZ,RZM,VXJ,VYJ,VZJ,QMULT,NPR) C----------------------------------------------------------------------- PARAMETER (MZ=128,NPRT=12800) C*** INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) DIMENSION QJX(MZ),QJY(MZ),QJZ(MZ) DIMENSION RZM(NPRT),VXJ(NPRT),VYJ(NPRT),VZJ(NPRT) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C C*** C DO 100 K= 1,MZ QJX(K)= 0. QJY(K)= 0. QJZ(K)= 0. 100 CONTINUE C C DO 200 L= 1,NPR KK= RZM(L) K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZM(L) -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C QGAMX = QMULT*VXJ(L) QGAMY = QMULT*VYJ(L) QGAMZ = QMULT*VZJ(L) C QJX(M,KL)= QJX(M,KL) + QGAMX*FZL QJX(M,K )= QJX(M,K ) + QGAMX*FZC QJX(M,KR)= QJX(M,KR) + QGAMX*FZR C QJY(M,KL)= QJY(M,KL) + QGAMY*FZL QJY(M,K )= QJY(M,K ) + QGAMY*FZC QJY(M,KR)= QJY(M,KR) + QGAMY*FZR C QJZ(M,KL)= QJZ(M,KL) + QGAMZ*FZL QJZ(M,K )= QJZ(M,K ) + QGAMZ*FZC QJZ(M,KR)= QJZ(M,KR) + QGAMZ*FZR 200 CONTINUE C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE SRIMP2 (QQ,RZM,QMULT,NPR) C----------------------------------------------------------------------- PARAMETER (MZ=128,NPRT=12800) C*** DIMENSION QQ(MZ),RZM(NPRT) INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C C*** C DO 100 K= 1,MZ 100 QQ(K)= 0. C C DO 200 L= 1,NPR KK= RZM(L) K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZM(L) -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C QQ(M,KL)= QQ(M,KL) + QMULT*FZL QQ(M,K )= QQ(M,K ) + QMULT*FZC QQ(M,KR)= QQ(M,KR) + QMULT*FZR 200 CONTINUE C RETURN END C C C----------------------------------------------------------------------- SUBROUTINE SRIMP3 (QJH,RZM,VH,QMULT,NPR) C----------------------------------------------------------------------- PARAMETER (MZ=128,NPRT=12800) C*** DIMENSION QJH(MZ),RZM(NPRT),VH(NPRT) C INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C C*** C DO 100 K= 1,MZ 100 QJH(K)= 0. C C DO 200 L= 1,NPR KK= RZM(L) K = PZC(KK) KL= PZL(KK) KR= PZR(KK) C ZZ= RZM(L) -KK -0.5001 FZL= 0.5*(0.5-ZZ)*(0.5-ZZ) FZC= 0.75-ZZ*ZZ FZR= 0.5*(0.5+ZZ)*(0.5+ZZ) C QMULVH= QMULT*VH(L) QJH(M,KL)= QJH(M,KL) + QMULVH*FZL QJH(M,K )= QJH(M,K ) + QMULVH*FZC QJH(M,KR)= QJH(M,KR) + QMULVH*FZR 200 CONTINUE C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE EMFLD0 C----------------------------------------------------------------------- PARAMETER (MZ=128) C REAL*4 KZ INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) COMMON/FACTR/ KZ(MZ),FZC0(MZ),SGN(MZ),ISB(MZ) C COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) COMMON/KFIELD/ EKX(MZ),EKY(MZ),EKZ(MZ),BKX(MZ),BKY(MZ),BKZ(MZ) COMMON/SRIMPL/ QIX(MZ),QIY(MZ),QIZ(MZ),QJX(MZ),QJY(MZ),QJZ(MZ), * QEX(MZ),QEY(MZ),QEZ(MZ), QI(MZ), QB(MZ), QE(MZ), * QJH(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C DIMENSION QQ(MZ),QQX(MZ),QQY(MZ),QQZ(MZ) C C C*********************************************************************** C* 1. ELECTRIC AND MAGNETIC FIELD LOADING (T=0). * C*********************************************************************** C C DO 100 K= 1,MZ QQX(K)= QIX(K) +QEX(K) QQY(K)= QIY(K) +QEY(K) QQZ(K)= QIZ(K) +QEZ(K) QQ(K)= QI(K) +QE(K) 100 CONTINUE C IFZ= -1 CALL FFT84 (QQX,IFZ) CALL FFT84 (QQY,IFZ) CALL FFT84 (QQZ,IFZ) CALL FFT84 ( QQ,IFZ) C C*** C DO 200 N= 1,MZ RQFILT= FZC0(N)/Q0 C QQX(N)= RQFILT*QQX(N) QQY(N)= RQFILT*QQY(N) QQZ(N)= RQFILT*QQZ(N) QQ(N)= RQFILT* QQ(N) 200 CONTINUE C C *VOCL LOOP, NOVREC (BKX,BKY,EKZ) C* DO 300 N= 2,MZ/2 NI= 2*N NR= NI-1 RKZ= 1./KZ(NI) C BKX(NR)= QQY(NI)*RKZ BKX(NI)= -QQY(NR)*RKZ C BKY(NR)= -QQX(NI)*RKZ BKY(NI)= QQX(NR)*RKZ C EKZ(NR)= QQ(NI)*RKZ EKZ(NI)= -QQ(NR)*RKZ 300 CONTINUE C C C*********************************************************************** C* 2. FILTER THE FIELD TO GET ONE FELT BY THE PARTICLES * C*********************************************************************** C DO 400 N= 1,MZ EZ(N)= FZC0(N)*EKZ(N) BX(N)= FZC0(N)*BKX(N) BY(N)= FZC0(N)*BKY(N) 400 CONTINUE C IFZ= 1 CALL FFT84 (EZ,IFZ) CALL FFT84 (BX,IFZ) CALL FFT84 (BY,IFZ) C C DO 500 K= 1,MZ BX(K)= BXC +BX(K) BY(K)= BYC +BY(K) BZ(K)= BZC C EX0(K)= EX(K) EY0(K)= EY(K) EZ0(K)= EZ(K) BX0(K)= BX(K) BY0(K)= BY(K) BZ0(K)= BZ(K) 500 CONTINUE C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE EMFILD C----------------------------------------------------------------------- C*********************************************************************** C* E*B CURRENT (OF ELECTRONS) IS INCLUDED IN THE LHS OF THE * C* FIELD-PARTICLE COUPLED EQUATION. THE LONGITUDINAL AND * C* TRANSVERSE PARTS ARE NOT SEPARATED IN THIS VERSION. * C*********************************************************************** C PARAMETER (MZ=128) C REAL*4 KZ INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) COMMON/FACTR/ KZ(MZ),FZC0(MZ),SGN(MZ),ISB(MZ) C COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) COMMON/KFIELD/ EKX(MZ),EKY(MZ),EKZ(MZ),BKX(MZ),BKY(MZ),BKZ(MZ) COMMON/SRIMPL/ QIX(MZ),QIY(MZ),QIZ(MZ),QJX(MZ),QJY(MZ),QJZ(MZ), * QEX(MZ),QEY(MZ),QEZ(MZ), QI(MZ), QB(MZ), QE(MZ), * QJH(MZ) C COMMON/EKWORK/ EKX1(MZ),EKY1(MZ),EKZ1(MZ),BKX1(MZ),BKY1(MZ), * BKZ1(MZ) COMMON/EMITER/ ASE,ASB,ASEL,WE,WB,WL,ITERM,ITERS C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C DIMENSION SXI(MZ),SYI(MZ),SZI(MZ),EX1(MZ),EY1(MZ),EZ1(MZ) C C* EKX, EX0 .... KEEP; FIELD OF THE TIME STEP (N). C* EKX1, EX .... WORK; USED FOR TIME STEP (N+1). C C*********************************************************************** C* 1. CALCULATE THE SOURCE TERMS. * C*********************************************************************** C----------------------------------------------------------------------- C GUESS TO THE ITERATION (USE PREVIOUS-STEP VALUE). C----------------------------------------------------------------------- C ITERM= 0 WEM= 0. ASE0= 0. ASB0= 0. C DO 10 N= 1,MZ EKX1(N)= EKX(N) EKY1(N)= EKY(N) EKZ1(N)= EKZ(N) 10 CONTINUE C 1000 ITERM = ITERM +1 C C----------------------------------------------------------------------- C* CURRENT DENSITY JX(N+AIMPL) -JX0, ETC. ARE CALCULATED. C----------------------------------------------------------------------- C HDTQW= HDT*QWI BCSQR= SQRT(BXC**2 +BYC**2 +BZC**2) C DO 100 K= 1,MZ EXA = AIMPL*EX(K) +(1.-AIMPL)*EX0(K) EYA = AIMPL*EY(K) +(1.-AIMPL)*EY0(K) EZA = AIMPL*EZ(K) +(1.-AIMPL)*EZ0(K) BXA = AIMPL*BX(K) +(1.-AIMPL)*BX0(K) BYA = AIMPL*BY(K) +(1.-AIMPL)*BY0(K) BZA = AIMPL*BZ(K) +(1.-AIMPL)*BZ0(K) C BXH = 0.5*(BX(K) +BX0(K)) BYH = 0.5*(BY(K) +BY0(K)) BZH = 0.5*(BZ(K) +BZ0(K)) C BSQA= BXA**2 +BYA**2 +BZA**2 BSQH= BXH**2 +BYH**2 +BZH**2 BSQ0= BX0(K)**2 +BY0(K)**2 +BZ0(K)**2 BSA= SQRT(BSQA) BSH= SQRT(BSQH) BS0= SQRT(BSQ0) C VIX= QJX(K) +HDTQW*( QI(K)*EXA +QJY(K)*BZA -QJZ(K)*BYA * +HDTQW*QI(K)*(EYA*BZA -EZA*BYA) ) VIY= QJY(K) +HDTQW*( QI(K)*EYA +QJZ(K)*BXA -QJX(K)*BZA * +HDTQW*QI(K)*(EZA*BXA -EXA*BZA) ) VIZ= QJZ(K) +HDTQW*( QI(K)*EZA +QJX(K)*BYA -QJY(K)*BXA * +HDTQW*QI(K)*(EXA*BYA -EYA*BXA) ) C EHA= (EXA*BXA +EYA*BYA +EZA*BZA )/BSA EH0= (EX0(K)*BX0(K) +EY0(K)*BY0(K) +EZ0(K)*BZ0(K))/BS0 C C SXI(K)= QIX(K) +QEX(K) * +GDT*QWI*( QI(K)*(EXA -EX0(K)) -AQI0*EX(K) * +(VIY*BZA -VIZ*BYA)/(1.+HDTQW**2*BSQA) * -(QJY(K)*BZ0(K) -QJZ(K)*BY0(K)) ) * * +QJH(K)*(BXH/BSH -BX0(K)/BS0) * +GDT*QWE*QE(K)*( EHA*BXH/BSH -EH0*BX0(K)/BS0 ) * +QE(K)*((EYA*BZA -EZA*BYA )/BSQA * -(EY0(K)*BZ0(K)-EZ0(K)*BY0(K))/BSQ0) * -AIMPL*QE0*EY(K)/BCSQR C SYI(K)= QIY(K) +QEY(K) * +GDT*QWI*( QI(K)*(EYA -EY0(K)) -AQI0*EY(K) * +(VIZ*BXA -VIX*BZA)/(1.+HDTQW**2*BSQA) * -(QJZ(K)*BX0(K) -QJX(K)*BZ0(K)) ) * * +QJH(K)*(BYH/BSH -BY0(K)/BS0) * +GDT*QWE*QE(K)*( EHA*BYH/BSH -EH0*BY0(K)/BS0 ) * +QE(K)*((EZA*BXA -EXA*BZA )/BSQA * -(EZ0(K)*BX0(K)-EX0(K)*BZ0(K))/BSQ0) * +AIMPL*QE0*EX(K)/BCSQR C SZI(K)= QIZ(K) +QEZ(K) * +GDT*QWI*( QI(K)*(EZA -EZ0(K)) -AQI0*EZ(K) * +(VIX*BYA -VIY*BXA)/(1.+HDTQW**2*BSQA) * -(QJX(K)*BY0(K) -QJY(K)*BX0(K)) ) * * +QJH(K)*(BZH/BSH -BZ0(K)/BS0) * +GDT*QWE*( QE(K)*(EHA*BZH/BSH -EH0*BZ0(K)/BS0) * -AQE0*EZ(K) ) * +QE(K)*((EXA*BYA -EYA*BXA )/BSQA * -(EX0(K)*BY0(K)-EY0(K)*BX0(K))/BSQ0) 100 CONTINUE C IFZ= -1 CALL FFT84 (SXI,IFZ) CALL FFT84 (SYI,IFZ) CALL FFT84 (SZI,IFZ) C C C*********************************************************************** C* 2. THE ELECTRIC FIELD IS SOLVED. * C* <> IN 1-D CASE, (EX,EY)... TRANSVERSE. * C* EZ....... LONGITUDINAL. * C*********************************************************************** C C DO 200 N= 1,MZ AF= 1. -AIMPL*(1.-AIMPL)*(KZ(N)*DT)**2 DTRQ= FZC0(N)*DT/Q0 C NS= ISB(N) SXI(N)= -DTRQ*SXI(N) +AF*EKX(N) +SGN(N)*DT*KZ(N)*BKY(NS) SYI(N)= -DTRQ*SYI(N) +AF*EKY(N) -SGN(N)*DT*KZ(N)*BKX(NS) C SZI(N)= -DTRQ*SZI(N) +EKZ(N) 200 CONTINUE C C DO 230 N= 1,MZ FILT2= FZC0(N)**2 AA= 1. +(ADT*KZ(N))**2 +FILT2*AGDT2*QQWI CC= FILT2*QSPEC(2)*ADT/BCSQR C EX1(N)= (AA*SXI(N) -CC*SYI(N))/(AA**2 +CC**2) EY1(N)= (CC*SXI(N) +AA*SYI(N))/(AA**2 +CC**2) C EZ1(N)= SZI(N)/(1. +FILT2*AGDT2*(QQWI+QQWE)) 230 CONTINUE C EX1(2)= 0. EY1(2)= 0. EZ1(2)= 0. C C C----------------------------------------------------------------------- C THIS AVERAGING PROCEDURE (OF TWO ITERATION VALUES) IS NECESSARY C TO HAVE CONVERGENCE. C----------------------------------------------------------------------- C CMI= 0.5 IF(WEM.LT.10.*EPSLN1) CMI= 1. -0.05*WEM/EPSLN1 C DO 260 N= 1,MZ EKX1(N)= CMI*EX1(N) +(1.-CMI)*EKX1(N) EKY1(N)= CMI*EY1(N) +(1.-CMI)*EKY1(N) EKZ1(N)= CMI*EZ1(N) +(1.-CMI)*EKZ1(N) 260 CONTINUE C C C*********************************************************************** C* 3. THE MAGNETIC FIELD IS CALCULATED (K-SPACE). * C*********************************************************************** C BKZ= 0. IN 1-D CASE. C EP1= AIMPL*DT EP0= (1.-AIMPL)*DT C DO 300 N= 1,MZ NS= ISB(N) BKX1(N)= BKX(N) -SGN(N)*KZ(N)*(EKY1(NS)*EP1 +EKY(NS)*EP0) BKY1(N)= BKY(N) +SGN(N)*KZ(N)*(EKX1(NS)*EP1 +EKX(NS)*EP0) 300 CONTINUE C BKX1(2)= 0. BKY1(2)= 0. C C C*********************************************************************** C* 4. FILTER THE FIELD TO GET ONE FELT BY THE PARTICLES * C*********************************************************************** C C DO 400 N= 1,MZ EX(N)= FZC0(N)*EKX1(N) EY(N)= FZC0(N)*EKY1(N) EZ(N)= FZC0(N)*EKZ1(N) BX(N)= FZC0(N)*BKX1(N) BY(N)= FZC0(N)*BKY1(N) 400 CONTINUE C C IFZ= 1 CALL FFT84 (EX,IFZ) CALL FFT84 (EY,IFZ) CALL FFT84 (EZ,IFZ) CALL FFT84 (BX,IFZ) CALL FFT84 (BY,IFZ) C C DO 430 K= 1,MZ BX(K)= BXC +BX(K) BY(K)= BYC +BY(K) BZ(K)= BZC 430 CONTINUE C C C*********************************************************************** C* 5. CONVERGENCE CHECK OF THE INNER LOOP. * C*********************************************************************** C ASE= 0. ASB= 0. DO 500 K= 1,MZ ASE= ASE +EX(K)**2 +EY(K)**2 +EZ(K)**2 ASB= ASB +(BX(K)-BXC)**2 +(BY(K)-BYC)**2 +(BZ(K)-BZC)**2 500 CONTINUE C ASE= ASE/FLOAT(MZ) ASB= ASB/FLOAT(MZ) C WE= 0. WB= 0. IF(ASE.GT.1.E-13) WE= ABS((ASE-ASE0)/ASE) IF(ASB.GT.1.E-13) WB= ABS((ASB-ASB0)/ASB) C WEM= AMAX1(WE,WB) ASE0= ASE ASB0= ASB C IF(MOD(IT,NPLOT).EQ.0) THEN WRITE(6,777) ITERM,WE,WB,ASE,ASB 777 FORMAT(1H ,' **ITERM=',I3,', WE,WB=',1P2E12.3, * ', ASE,ASB=',2E12.3) END IF C IF((WEM.GT.EPSLN1).AND.(ITERM.LT.20)) GO TO 1000 C ####### LOOP 1000 ENDS HERE ######## C RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE ESCORR C----------------------------------------------------------------------- C******** IMPLICIT ES SOLVER ******************************************* C* CORRECTION TO THE SCALAR POTENTIAL FIELD. * C*********************************************************************** C PARAMETER (MZ=128) C*** REAL*4 KZ INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) COMMON/FACTR/ KZ(MZ),FZC0(MZ),SGN(MZ),ISB(MZ) C COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) COMMON/KFIELD/ EKX(MZ),EKY(MZ),EKZ(MZ),BKX(MZ),BKY(MZ),BKZ(MZ) COMMON/SRIMPL/ QIX(MZ),QIY(MZ),QIZ(MZ),QJX(MZ),QJY(MZ),QJZ(MZ), * QEX(MZ),QEY(MZ),QEZ(MZ), QI(MZ), QB(MZ), QE(MZ), * QJH(MZ) COMMON/EKWORK/ EKX1(MZ),EKY1(MZ),EKZ1(MZ),BKX1(MZ),BKY1(MZ), * BKZ1(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) COMMON/EMITER/ ASE,ASB,ASEL,WE,WB,WL,ITERM,ITERS DIMENSION POT(MZ),POTK(MZ),PTK1(MZ),QW(MZ),QRZ(MZ) C C*********************************************************************** C* 1. FIND THE CORRECTION SOURCE TERM. * C*********************************************************************** C C DO 10 K= 1,MZ POT(K)= 0. POTK(K)= 0. 10 CONTINUE C HDTQW= HDT*QWI ITERS= 0 WL= 1. ASEL0= 0. C 1 ITERS = ITERS + 1 C C DO 100 K= 1,MZ EZL = -(POT(PZR(K))-POT(PZL(K)))/HZ2 C BXA = AIMPL*BX(K) +(1.-AIMPL)*BX0(K) BYA = AIMPL*BY(K) +(1.-AIMPL)*BY0(K) BZA = AIMPL*BZ(K) +(1.-AIMPL)*BZ0(K) BSQA= BXA**2 +BYA**2 +BZA**2 C BXH = 0.5*(BX(K) +BX0(K)) BYH = 0.5*(BY(K) +BY0(K)) BZH = 0.5*(BZ(K) +BZ0(K)) BSQH= BXH**2 +BYH**2 +BZH**2 C C* FOR EL. C EHA= EZL*BZA/SQRT(BSQA) ETX= -EZL*BYA/BSQA ETY= EZL*BXA/BSQA C C* FOR IONS. C EX1= HDTQW*BSQA*ETX/(1.+HDTQW**2*BSQA) EY1= HDTQW*BSQA*ETY/(1.+HDTQW**2*BSQA) C EZI= EZL +HDTQW*(EX1*BYA -EY1*BXA) C C C----------------------------------------------------------------------- C* CORRECTION DUE TO DIFFERENCE OF TWO PIVOTS ARE NOT MADE. C----------------------------------------------------------------------- C QRZ(K)= AHDT2*( QWI*QI(K)*EZI +QWE*QE(K)*EHA*BZH/SQRT(BSQH)) 100 CONTINUE C C DO 130 K= 1,MZ QW(K)= QI(K) +QE(K) -(QRZ(PZR(K))-QRZ(PZL(K)))/HZ2 130 CONTINUE C C IFZ=-1 CALL FFT84 (QW,IFZ) C C C*********************************************************************** C* 2. CALCULATE THE CORRECTION TO THE SCALAR POTENTIAL FIELD. * C*********************************************************************** C DO 200 N= 3,MZ AKSQ= KZ(N)**2 AQQW= AHDT2*(QQWI +QQWE)*AKSQ C PTK1(N)= ( FZC0(N)*(QW(N)/Q0 +FZC0(N)*AQQW*POTK(N)) * +SGN(N)*KZ(N)*EKZ1(ISB(N)) )/(AKSQ +FZC0(N)**2*AQQW) 200 CONTINUE C PTK1(1)= 0. PTK1(2)= 0. C C C----------------------------------------------------------------------- C* AVOID TOO MUCH CORRECTION AT A TIME. C----------------------------------------------------------------------- C CPI= 0.6 IF((WL.LT.10.*EPSLN1).OR.(ITERS.EQ.1)) CPI= 1.0 C DO 230 N= 1,MZ POTK(N)= CPI*PTK1(N) +(1.-CPI)*POTK(N) 230 CONTINUE C C C*********************************************************************** C* 3. CONVERGENCE CHECK OF POT (ELECTRIC FIELD ENERGY). * C*********************************************************************** C ASEL= 0. DO 300 N= 1,MZ ASEL= ASEL +(EKZ1(N) +SGN(N)*KZ(N)*POTK(ISB(N)))**2 300 CONTINUE C ASEL= ASEL/FLOAT(MZ) C WL= 0. IF(ASEL.GT.1.E-13) THEN WL= ABS((ASEL-ASEL0)/ASEL) ASEL0= ASEL END IF C C C*********************************************************************** C* 4. THE SCALAR POTENTIAL FIELD IN X-SPACE (POT). * C*********************************************************************** C DO 400 N= 1,MZ POT(N)= FZC0(N)*POTK(N) 400 CONTINUE C IFZ= 1 CALL FFT84 (POT,IFZ) C IF(MOD(IT,NPLOT).EQ.0) THEN WRITE(6,444) ITERS,WL,ASEL 444 FORMAT(1H ,' **ITERS=',I3,', WL,ASEL= ',1P2E12.3) END IF C IF((WL.GT.EPSLN1).AND.(ITERS.LT.20)) GO TO 1 C ############ LOOP 1 ENDS HERE ########### C C C C*********************************************************************** C* 5. FINAL STEP TO ADVANCE TO THE NEW TIME STEP. * C* CORRECTION TO THE ELECTRIC FIELD IS MADE FOR BOTH EX & EKX. * C*********************************************************************** C C DO 500 K= 1,MZ EZ(K)= EZ(K) -(POT(PZR(K)) -POT(PZL(K)))/HZ2 500 CONTINUE C C DO 530 N= 1,MZ EKX(N)= EKX1(N) EKY(N)= EKY1(N) EKZ(N)= EKZ1(N) +SGN(N)*KZ(N)*POTK(ISB(N)) C BKX(N)= BKX1(N) BKY(N)= BKY1(N) 530 CONTINUE C C IF(MOD(IT,NPLOT).EQ.0) THEN WRITE(6,777) IT WRITE(6,778) ITERM,ITERS,ASE,ASB,ASEL,WE,WB,WL C 777 FORMAT(1H ,'----- IT=',I5,' -------------------------------------- *-------------------------------------------------------------') 778 FORMAT(1H ,'*ITERM, ITERS=',2I3, * /,' ,,=',1P3E12.3, * ' >>> WE,WB,WL=',3E12.3,' <<<',/) END IF C C RETURN END C C C C---------------------------------------------------------------------- SUBROUTINE INIT (KSTART) C----------------------------------------------------------------------- C* PARAMETER (MZ=128,NPRT=12800) C* INTEGER*4 PZR,PZC,PZL COMMON/TABLE/ PZR(-MZ:2*MZ),PZC(-MZ:2*MZ),PZL(-MZ:2*MZ) C REAL*4 KZ COMMON/FACTR/ KZ(MZ),FZC0(MZ),SGN(MZ),ISB(MZ) C COMMON/FIELDS/ EX(MZ),EY(MZ),EZ(MZ),BX(MZ),BY(MZ),BZ(MZ) COMMON/FIELD0/ EX0(MZ),EY0(MZ),EZ0(MZ),BX0(MZ),BY0(MZ),BZ0(MZ) COMMON/KFIELD/ EKX(MZ),EKY(MZ),EKZ(MZ),BKX(MZ),BKY(MZ),BKZ(MZ) C COMMON/PARM1/ IT,LDEC,NSP,NSPEC(2),NFWRT,NPWRT,NHA,NPLOT,NHIST COMMON/PARM2/ ZMAX,HZ,HZI,HZ2,ZMAXE,AZ,QSPEC(2),WSPEC(2),VETH, * TETI,TANIS,WCEWPE,THB,PI,AIMPL,GIMPL,EPSLN1, * Q0,QI0,QE0,AQI0,AQE0,QWI,QWE,AQWI,AQWE,QQWI,QQWE, * ADT,GDT,AHDT2,AGDT2,VTHX(2),VTHZ(2),VDR(2), * VBEAM(2),EFE,EFB,ETOT0,T,DT,HDT,BXC,BYC,BZC, * VLIMA,VLIMB,BMIN,EMIN,EDEC(256,23) C REAL*4 MUE COMMON/PARTI/ ZI(NPRT),VXI(NPRT),VYI(NPRT),VZI(NPRT) COMMON/PARTE/ ZE(NPRT),MUE(NPRT),VPE(NPRT),VHE(NPRT) C C****** C HZ= ZMAX/MZ HZI= 0.99999/HZ HZ2= 2.0*HZ ZMAXE= 0.99999*ZMAX C PI2= 2.0*PI C C C----------------------------------------------------------------------- C* INDEX TABLES USED FOR PARTICELS AND FIELDS. C----------------------------------------------------------------------- C DO 20 K=1,MZ PZL(K)= K-1 PZC(K)= K PZR(K)= K+1 20 CONTINUE C PZL(1)= MZ PZR(MZ)= 1 C DO 23 K= 1,MZ PZL(K-MZ)= PZL(K) PZC(K-MZ)= PZC(K) PZR(K-MZ)= PZR(K) 23 CONTINUE C DO 24 K= 1,MZ PZL(K+MZ)= PZL(K) PZC(K+MZ)= PZC(K) PZR(K+MZ)= PZR(K) 24 CONTINUE C C* ION TEMPERATURE ANISOTROPY..... TANIS. C VITH= VETH/SQRT(TETI*WSPEC(1)) VTHX(1)= SQRT(TANIS)*VITH VTHZ(1)= VITH VTHX(2)= VETH VTHZ(2)= VETH C THB1= PI*THB/180.0 SN= SIN(THB1) CS= COS(THB1) BXC= 0. BYC= WCEWPE*SN BZC= WCEWPE*CS C C C----------------------------------------------------------------------- C* WAVENUMBER TABLES (IGNORE N=2). C----------------------------------------------------------------------- C TWO STORAGES OF KZ(N) CORRESPOND TO ONE FOURIER MODE: C N: 1 2 3 4 5 6 7... C MODE : 0 MZ 1 1 2 2 3... C C DO 730 N= 2,MZ/2 NI= 2*N NR= NI-1 KZ(NR)= PI2*(N-1)/ZMAX KZ(NI)= KZ(NR) 730 CONTINUE C KZ(1)= 0. KZ(2)= PI C DO 830 N= 3,MZ HKZ= 0.5*KZ(N)*HZ C0= (SIN(HKZ)/HKZ)**3/(1.0-SIN(HKZ)**2+(2./15.)*SIN(HKZ)**4) FZC0(N)= C0*EXP(-(KZ(N)*AZ)**8) 830 CONTINUE C FZC0(1)= 1. FZC0(2)= 0. C WRITE(6,801) WRITE(6,803) FZC0 801 FORMAT(1H0,'** FZC0 **',/) 803 FORMAT(1H ,1P10E12.4) C C DO 860 N= 1,MZ/2 NI= 2*N NR= NI-1 SGN(NR)= 1. SGN(NI)= -1. ISB(NR)= NI ISB(NI)= NR 860 CONTINUE C C WRITE(6,*) '** SET FFT2 **' CALL SETFFT C C C ### PARTICLE LOADING ######################################## C IF(KSTART.NE.0) RETURN C IT= 0 T = 0. LDEC= 0 C C CALL PARTLD (ZI,VXI,VYI,VZI,1,NSPEC(1)) CALL PARTLD (ZE,MUE,VPE,VHE,2,NSPEC(2)) C Q0= ABS(Q0)/FLOAT(MZ) C C CALL CLOCKS (CPUIL) WRITE(6,*) ' PARTICLE LOAD..... COMPLETED.' WRITE(6,*) ' ** CPU (PARTICLE LOAD) ',CPUIL C C C----------------------------------------------------------------------- C* CONSTANTS USED IN FIELD SOLUTION. C----------------------------------------------------------------------- C C* <> Q0 IS DEFINED JUST ABOVE. C QI0= QSPEC(1)*Q0 QE0= QSPEC(2)*Q0 QWI= QSPEC(1)/WSPEC(1) QWE= QSPEC(2)/WSPEC(2) QQWI= QSPEC(1)*QWI QQWE= QSPEC(2)*QWE C ADT= AIMPL*DT GDT= GIMPL*DT HDT= 0.5*DT AHDT2= ADT*HDT AGDT2= ADT*GDT C AQI0= AIMPL*QI0 AQE0= AIMPL*QE0 AQWI= AIMPL*QWI AQWE= AIMPL*QWE C C C----------------------------------------------------------------------- C INITIAL FIELDS (RESET). C----------------------------------------------------------------------- C DO 70 K= 1,MZ EKX(K)= 0. EKY(K)= 0. EKZ(K)= 0. BKX(K)= 0. BKY(K)= 0. BKZ(K)= 0. C EX(K)= 0. EY(K)= 0. EZ(K)= 0. BX(K)= BXC BY(K)= BYC BZ(K)= BZC C EX0(K)= 0. EY0(K)= 0. EZ0(K)= 0. BX0(K)= BXC BY0(K)= BYC BZ0(K)= BZC 70 CONTINUE C C CALL CLOCKS (CPU) WRITE(6,781) CPU 781 FORMAT(1H0,'**PARTICLE LOAD COMPLETED** CPU=',F9.2,' SECS', * //,'**INIT CALL TO TRANS (MOVER, FIELD, DIAG)***') C CALL TRANS C RETURN END ========== end of macroem.f =================================