Computational Methods for Kinetic Processes in Plasma Physics

Computational Methods for Kinetic Processes in Plasma Physics Ken Nishikawa National Space Science & Technology Center/UAH October 18, 2011 1/39

Context Electromagnetic code new_bbeps 1. f, a magnetized electromagnetic code 1. understanding the code 2. calculates 3 d phase space distribution, velocity moments, and entropy (PSDIST 13) add section to use this subroutine in newbbeps 1. f 3. how to add contour plots 4. call subroutine to plot phase space contour 5. understanding push 1 d

Kernel code: new_bbeps 1. f use use 1. 2. 3. 4. 5. 6. 7. 8. init 1 d bpush 1 d dpush 1 d fft 1 d field 1 d diag 1 d push 1 d init 1 d. mod bpush 1 d. mod dpush 1 d. mod fft 1 d. mod field 1 d. mod diag 1 d. mod bpush 1 d. mod init 1 mod. f bpush 1 mod. f dpush 1 mod. f fft 1 mod. f field 1 mod. f diag 1 mod. f bpush 1 mod. f init 1 lib. f (fort 77) bpush 1 lib. f dpush 1 lib. f fft 1 lib. f field 1 lib. f diag 1 lib. f bpush 1 lib. f understand the structure of the code understand the essential subroutines understand init 1 lib. f and init 1 mod. f understand diag 1 lib. f and diag 1 mod. f to implement new diagnostics calculates 3 d phase space distribution add diagnostics of phase velocity plot (contour in color like k-ω plot) how to add phase plot using dots how to push particles (collect charge)

!-----------------------------------! * * * periodic 1 d electromagnetic particle simulation kernel code * * ! this is a simple 1 d skeleton particle-in-cell code designed for ! exploring new computer architectures. it contains the critical pieces ! needed for depositing charge and current, advancing particles, and ! solving the fields. the code moves only electrons, with periodic ! electromagnetic forces obtained by solving maxwell's equation with ! fast fourier transforms. ! written by viktor k. decyk, ucla ! Fortran 90 for Macintosh G 3 ! copyright 1999, regents of the university of california ! update: july 14, 2011 ! modules below are included (ken) ! we need to add nxb, nmvf program bbeps 1 use init 1 d initialize particle position and velocity use bpush 1 d use dpush 1 d use push 1 d attention use fft 1 d use field 1 d use diag 1 d produce plots

implicit none ! idimp = dimension of phase space = 4 ! nmv = number of segments in v for velocity distribution integer : : idimp = 4, nmv = 40, ipbc = 1 ! default unit numbers integer : : iuin = 8, iuot = 18, iudm = 19, iup = 11, iua = 15 integer : : iue = 26 integer : : np, npx 1, nxh, nxeh, nx 1 integer : : nloop, itime 0, ntime, ltime, isign, irc integer : : it, itw real : : zero = 0. 0, time = 0. 0, tloop = 0. 0 real : : tpush = 0. 0, tdpost = 0. 0, tdjpost = 0. 0, tsort = 0. 0 real : : tfft = 0. 0, totpush = 0. 0, ts = 0. 0 real : : qbme, affp, dth, qi 0, omt, q 2 m 0, wp 0 real : : we, wf, wm, wef, wke, ws real, dimension(: , : ), pointer : : part real, dimension(: ), pointer : : qe, fxe real, dimension(: , : ), pointer : : cu, amu, fxyze, byze complex, dimension(: , : ), pointer : : eyz, byz complex, dimension(: ), pointer : : ffc, ffe

integer, dimension(: ), pointer : : mixup complex, dimension(: ), pointer : : sct real, dimension(: ), pointer : : pt integer, dimension(: ), pointer : : ip, npic real, dimension(: ), pointer : : sfield complex, dimension(: ), pointer : : pott real, dimension(: , : ), pointer : : vfield complex, dimension(: , : ), pointer : : vpott, vpotr real, dimension(: , : ), pointer : : fv, fvm real, dimension(: , : ), pointer : : wt character(len=10) : : cdrun character(len=32) : : fname 991 format (' T = ', i 7) 992 format (' field, kinetic, total energies = ', 3 e 14. 7) 993 format (' electric(l, t), magnetic energies = ', 3 e 14. 7) ! read namelist iuin = get_funit(iuin) open(unit=iuin, file='input 1', form='formatted', status='old') read (iuin, input 1) read parameters in input 1 ! override input data

idcode = 2 psolve = 1 ndim = 2 ! create string from idrun write (cdrun, '(i 10)') idrun cdrun = adjustl(cdrun) ! text output file iuot = get_funit(iuot) fname = 'output 1. '//cdrun open(unit=iuot, file=trim(fname), form='formatted', status='replace') ! np = total number of electrons in simulation np = npx + npxb; npx 1 = npx + 1 npx: ambient, npxb: beam nx = 2**indx; nxh = nx/2 nxe = nx + 4 ! ax =. 866025 if (inorder==LINEAR) then ! ax =. 912871 nxe = nx + 2 endif nxeh = nxe/2

! dimension for index and sorting arrays nx 1 = nx + 1 ! nloop = number of time steps in simulation nloop = tend/dt +. 0001 ! part(1, n) = position x of particle n ! part(2, n) = velocity vx of particle n initialized by subroutine DISTR 1 H ! part(3, n) = velocity vy of particle n ! part(4, n) = velocity vz of particle n allocate(part(idimp, np)) set the dimension ! in real space, qe(j) = charge density at grid point j ! in real space, fxe(j) = longitudinal force/charge at grid point j, ! that is, fxe is the convolution of the longitudinal electric field ! over the particle shape allocate(qe(nxe), fxe(nxe)) ! in real space, fxyze(i, j) = i component of force/charge at grid (j) ! that is, fxyze are the convolutions of the electric field ! over the particle shape

allocate(fxyze(3, nxe)) ! cu(i, j) = i component of current at grid (j). ! byze(i, j) = i component of magnetic field at grid (j). ! byze is the convolution of the magnetic field over the particle shape allocate(cu(2, nxe), byze(2, nxe)) ! in fourier space, exyz = transverse electric field ! in fourier space, bxyz = magnetic field allocate(eyz(2, nxeh), byz(2, nxeh)) ! ffc = form factor array for poisson solver allocate(ffc(nxh)) ! mixup = array of bit reversed addresses for fft ! sct = sine/cosine table for fft allocate(mixup(nxh), sct(nxh)) ! ! open graphics device irc = open_graphs(nplot) enable to use GKS graphic package ! initialize timer call wtimer(time, ltime, -1)

! initialize constants itime 0 = 0 itime = itime 0 ntime = itime + itime 0 qbme = qme affp = real(nx)/real(np) omt = sqrt(omy*omy + omz*omz) q 2 m 0 = qbme*qme*real(np)/real(nx) dth = 0. 0 ! debug ! dth =. 5*dt ! end debug wp 0 = q 2 m 0*affp plasma frequency ! set initial time t 0 = dt*real(itime 0) ! set default diagnostic file names if (ntp > 0) fpname = 'potk 1. '//cdrun if (nta > 0) faname = 'vpotk 1. '//cdrun if (nte > 0) fename = 'vpotrk 1. '//cdrun

! energy diagnostics if (ntw > 0) then allocate(wt((nloop-1)/ntw-(itime 0/ntw)+1, 7)) itw = 0 endif ! initialize electromagnetic fields byze = 0. 0 if (omt > 0. 0) then call baddext(byze, omy, omz, nx, inorder) call cguard(byze, nx, inorder) endif byz = cmplx(0. 0, 0. 0) eyz = cmplx(0. 0, 0. 0) cu = 0. 0 ! prepare fft tables call fft_init(mixup, sct, indx) ! calculate form factors call pois_init(ffc, ax, affp, nx)

! initialize density profile and velocity distribution ! background electrons if (npx > 0) call distr(part, 1, npx, vty, vtz, vx 0, vy 0, vz 0, npx, nx, &ipbc) ! beam electrons if (npxb > 0) call distr(part, npx 1, npxb, vtdx, vtdy, vtdz, vdx, vdy, vdz& &, npxb, nx, ipbc) Remark! in init 1 mod. f distr is assined by `distr 1 h’ ! distr => idistrh 1 initializes x and vx, vy, vz co-ordinates for ! magnetized 1 -2/2 d codes. ! calls DISTR 1 H (see the next page) ! initialize charge density to background qi 0 = -qme/affp ! fix guiding centers for electrons if (omt > 0. 0) call distr(part, byze, np, qbme, nx, ipbc, inorder) ! retard electron positions to deposit current ! call retard(part, np, dth, nx) ! sorting arrays if (sortime > 0) allocate(pt(np), ip(np), npic(nx 1)) ! initialize diagnostics ! diagnostic metafile iudm = get_funit(iudm)

! velocity diagnostic if (ntv > 0) then allocate(fv(2*nmv+2, 3), fvm(3, 3)) fv(1, : ) = 8. *max(vtx, vty, vtz) endif ! potential diagnostic if (ntp > 0) then allocate(sfield(nxe)) if (modesxp > nxh) modesxp = nxh allocate(pott(modesxp)) ! open output file if (nprec==0) then nprec = -1; iup = get_funit(iup) call bfopen(pott, modesxp, iup, nprec, trim(fpname)) endif ! vector potential or electromagnetic diagnostic if ((nta > 0). or. (nte > 0)) then allocate(vfield(2, nxe)) vfield = 0. 0 endif

! vector potential diagnostic if (nta > 0) then if (modesxa > nxh) modesxa = nxh allocate(vpott(2, modesxa)) ! open output file if (narec==0) then narec = -1; iua = get_funit(iua) call bfopen(vpott, modesxa, iua, narec, trim(faname)) endif ! electromagnetic diagnostic if (nte > 0) then if (modesxe > nxh) modesxe = nxh allocate(vpotr(2, modesxe)) ! open output file if (nerec==0) then nerec = -1; iue = get_funit(iue) call bfopen(vpotr, modesxe, iue, nerec, trim(fename)) endif

! record time call wtimer(time, ltime) write (iuot, *) 'initialization wall clock time = ', time, 'sec' ! ! * * * start main iteration loop * * * attention ! 500 if (nloop <= ntime) go to 2000 write (iuot, 991) ntime ! ! prepare electromagnetic diagnostic if (nte > 0) then it = ntime/nte if (ntime==nte*it) vfield = cu endif ! initialize current density to background call sguard(cu, zero, nx, inorder) ! deposit electron current call djpost(part, cu, np, qme, dth, tdjpost, nx, ipbc, inorder, djopt)

! initialize charge density to background call sguard(qe, qi 0, nx, inorder) ! deposit electron charge call dpost(part, qe, np, qme, tdpost, inorder, dopt) ! add guard cells for current call aguard(cu, nx, inorder) ! add guard cells call aguard(qe, nx, inorder) ! velocity diagnostic if (ntv > 0) then for example ntv=5 in input 1. whistler it = ntime/ntv if (ntime==ntv*it) then ! calculate paticle distribution function and moments call vdist(part, fvm, np, nmv) ! display velocity distributions call displayfv(fv, fvm, ' ELECTRON', ntime, nmv, 2, irc) if (irc==1) go to 2000 endif

! phase velocity diagnostic by Victor ! nmv = number of segments in v for velocity distribution ! integer : : idimp = 4, nmv = 40, ipbc = 1 ! phase space diagnostic if (nts > 0) then for example ntv=5 in input 1. whistler it = ntime/nts if (ntime==nts*it) then ! calculate phase space distribution function and moments call psdist(part, fpsm, nx, np, nmv) ! display phase space distributions call displayfps(fps, fpsm, ' ELECTRON', ntime, nmv, 999, 2, irc) if (irc==1) go to 2000 endif !----- for comparison ------!for comparison ! call vdist(part, fvm, np, nmv) ! display velocity distributions call displayfv(fv, fvm, ' ELECTRON', ntime, nmv, 2, irc)

! transform charge to fourier space isign = -1 call fft(qe, isign, mixup, sct, tfft, indx, inorder) ! potential diagnostic if (ntp > 0) then it = ntime/ntp if (ntime==ntp*it) then ! calculate potential in fourier space isign = 1 call pois(qe, sfield, isign, ffc, ws, nx, inorder) ! store selected fourier modes call gtmodes(sfield, pott, nx, modesxp, inorder) ! write diagnostic output call writebf(pott, modesxp, iup, nprec, order=LINEAR) ! transform potential to real space call fft(sfield, isign, mixup, sct, tfft, indx, inorder) call cguard(sfield, nx, inorder)

! display potential call displays(sfield, ' POTENTIAL', ntime, 999, 0, nx, irc, inorder& &) if (irc==1) go to 2000 endif ! transform current to fourier space isign = -1 call fft(cu, isign, mixup, sct, tfft, indx, inorder) ! electromagnetic diagnostic if (nte > 0) then it = ntime/nte if (ntime==nte*it) then ! calculate averaged radiative vector potential vfield = 0. 5*(vfield + cu) call avrpot(vfield, byz, ffc, ci, nx, inorder) ! store selected fourier modes call gtmodes(vfield, vpotr, nx, modesxe, inorder) ! write diagnostic output call writebf(vpotr, modesxe, iue, nerec, order=LINEAR)

! transform radiative vector potential to real space isign = 1 call fft(vfield, isign, mixup, sct, tfft, indx, inorder) call cguard(vfield, nx, inorder) ! display radiative vector potential call displayv(vfield, ' RADIATIVE VPOTENTIAL', ntime, 999, 0, 1, & &nx, irc, inorder) if (irc==1) go to 2000 endif ! calculate electromagnetic fields in fourier space if (ntime==0) then ! calculate initial darwin magnetic field call ibpois(cu, byz, ffc, ci, wm, nx, inorder) wf = 0. ! calculate initial darwin electric field allocate(amu(2, nxe), ffe(nxeh))

! deposit momentum flux call sguard(amu, zero, nx, inorder) call dmjpost(part, amu, np, qme, ts, inorder, djopt) call aguard(amu, nx, inorder) ! solve for darwin electric field isign = -1 call fft(amu, isign, mixup, sct, tfft, indx, inorder) ! byze = cu call dcuperp(byze, amu, nx, inorder) call epois_init(ffe, ax, affp, wp 0, ci, nx) call iepois(byze, eyz, ffe, ci, wf, nx, inorder) deallocate(amu, ffe) dth =. 5*dt ! calculate electromagnetic fields else call maxwel(eyz, byz, cu, ffc, ci, dt, wf, wm, nx, inorder) endif

! vector potential diagnostic if (nta > 0) then it = ntime/nta if (ntime==nta*it) then ! calculate vector potential in fourier space call avpot(byz, vfield, nx, inorder) ! store selected fourier modes call gtmodes(vfield, vpott, nx, modesxa, inorder) ! write diagnostic output call writebf(vpott, modesxa, iua, narec, order=LINEAR) ! transform vector potential to real space isign = 1 call fft(vfield, isign, mixup, sct, tfft, indx, inorder) call cguard(vfield, nx, inorder) ! display vector potential call displayv(vfield, ' VECTOR POTENTIAL', ntime, 999, 0, 1, nx, & &irc, inorder) if (irc==1) go to 2000 endif

! calculate longitudinal electric field in fourier space isign = -1 call pois(qe, fxe, isign, ffc, we, nx, inorder) ! add longitudinal and transverse electric fields call emfield(fxyze, fxe, eyz, ffc, nx, inorder) ! copy magnetic field call emfield(byze, byz, ffc, nx, inorder) ! transform force/charge to real space isign = 1 call fft(fxyze, isign, mixup, sct, tfft, indx, inorder) call cguard(fxyze, nx, inorder) ! transform magnetic field to real space isign = 1 call fft(byze, isign, mixup, sct, tfft, indx, inorder) ! add external magnetic field if (omt > 0. 0) call baddext(byze, omy, omz, nx, inorder) call cguard(byze, nx, inorder) ! push particles attention wke = 0. call push 3(part, fxyze, byze, omx, np, qbme, dth, wke, tpush, nx, ipbc, & &inorder, popt)

! sort electrons if (sortime > 0) then if (mod(ntime, sortime)==0) then call sortp(part, pt, ip, npic, tsort, inorder) endif ! energy diagnostic if (ntw > 0) then it = itime/ntw if (itime==ntw*it) then wef = we + wf + wm ws = wef + wke write (iuot, 992) wef, wke, ws write (iuot, 993) we, wf, wm itw = itw + 1 wt(itw, : ) = (/wef, wke, zero, ws, we, wf, wm/) endif itime = itime + 1 ntime = itime + itime 0 attention

call wtimer(tloop, ltime) time = time + tloop go to 500 2000 continue ! ! * * * end main iteration loop * * * ! ! energy diagnostic if (ntw > 0) then ts = t 0 + dt*real(ntw)*(itime 0 -(itime 0/ntw)*ntw) call reset_graphs call displayw(wt, ts, dt*real(ntw), itw, irc) endif ! accumulate timings write (iuot, *) 'electromagnetic code bbeps 1' write (iuot, *) 'main wall clock time = ', time, 'sec' totpush = tpush + tdpost + tdjpost write (iuot, *) 'electron push time = ', tpush, 'sec' write (iuot, *) 'electron charge deposit time = ', tdpost, 'sec'

write (iuot, *) 'electron current deposit time = ', tdjpost, 'sec' write (iuot, *) 'total electron push time = ', totpush, 'sec' write (iuot, *) 'electron sort time = ', tsort totpush = totpush + tsort write (iuot, *) 'total electron time = ', totpush, 'sec' write (iuot, *) 'total fft time=', tfft, 'sec' time = time - (totpush + tfft) write (iuot, *) 'other time=', time, 'sec' ! write final diagnostic metafile fname = 'diag 1. '//cdrun open(unit=iudm, file=trim(fname), form='formatted', status='replace') ! potential diagnostics if (ntp > 0) then nprec = nprec - 1 ceng = affp write (iudm, pot 1 d) endif

! vector potential diagnostics if (nta > 0) then narec = narec - 1 ceng = affp write (iudm, vpot 1 d) endif ! electromagnetic diagnostics if (nte > 0) then nerec = nerec - 1 ceng = affp write (iudm, em 1 d) endif ! write out input file write (iudm, input 1) write (iuot, *) ' * * * q. e. d. * * *' close(unit=iudm) close(unit=iuot) ! close graphics device call close_graphs stop end program bbeps 1

PUSH 3 call push 3(part, fxyze, byze, omx, np, qbme, dth, wke, tpush, nx, ipbc, & &inorder, popt) part: part(3, n) particle position and three velocities fxyze: smoothed total electric field array byze: smoothed magnetic field array omx: magnetic field electron cyclotron frequency in x np: np=npx +npxb qbme: = qme: charge on electron, in unit of e, usually -1 dt: time interval between successive calcualtions dth: dth = 0. 5*dt wke: wke = 0: ws = wef + wke tpush: tpush = 0 for diagnostic of speed of calculation nx: number of grid = 2**indx ipbc: ipbc=1 inorder: type of interpolation; 1=LINEAR, 2=QUADRATIC popt: particle optimization scheme: 1=STANDARD, 2= LOOKAHEAD

In bpush 1 lib. f c-----------------------------------c 1 d PIC library for pushing particles with magnetic field c and depositing current c bpush 1 lib. f contains procedures to process particles with magnetic c fields: c GJPOST 1 deposits current density, quadratic interpolation, STANDARD c optimization. c GSJPOST 1 deposits current density, quadratic interpolation, LOOKAHEAD c optimization. c GSJPOST 1 X deposits current density, quadratic interpolation, VECTOR c optimization. c GJPOST 1 L deposits current density, linear interpolation, STANDARD c optimization. c GSJPOST 1 L deposits current density, linear interpolation, LOOKAHEAD c optimization. c GSJPOST 1 XL deposits current density, linear interpolation, VECTOR c optimization. c GBPUSH 13 push particles with magnetic field, quadratic interpolation, c STANDARD optimization.

c GSBPUSH 13 push particles with magnetic field, quadratic interpolation, c LOOKAHEAD optimization. c GPUSH 1 L push particles with magnetic field, linear interpolation, c STANDARD optimization. c GSBPUSH 13 L push particles with magnetic field, linear interpolation, c LOOKAHEAD optimization. c RETARD 1 retard particle position a half time-step. c written by viktor k. decyk, ucla c copyright 1994, regents of the university of california c update: july 3, 2010

c-----------------------------------subroutine GJPOST 1(part, cu, qm, dt, nop, idimp, nxv, ipbc) c for 1 -2/2 d code, this subroutine calculates particle current density c using second-order spline interpolation c in addition, particle positions are advanced a half time-step c scalar version using guard cells c 27 flops/particle, 10 loads, 7 stores c input: all, output: part, cu c current density is approximated by values at the nearest grid points c cu(i, n)=qci*(. 75 -dx**2) c cu(i, n+1)=. 5*qci*(. 5+dx)**2 c cu(i, n-1)=. 5*qci*(. 5 -dx)**2 c where n = nearest grid point and dx = x-n c and qci = qm*vi, where i = y, z c part(1, n) = position x of particle n c part(2, n) = x velocity of particle n c part(3, n) = y velocity of particle n c part(4, n) = z velocity of particle n c cu(i, j+1) = ith component of current density at grid point j c qm = charge on particle, in units of e c dt = time interval between successive calculations

c nop = number of particles c idimp = size of phase space = 4 c nxv = second dimension of current array, must be >= nx+3 c ipbc = particle boundary condition = (0, 1, 2) = c (none, 2 d periodic, 2 d reflecting) implicit none integer nop, idimp, nxv, ipbc real part, cu, qm, dt dimension part(idimp, nop), cu(2, nxv) c local data integer j, nn real qmh, edgelx, edgerx, dxp, amx, dxl, vy, vz, dx qmh =. 5*qm c set boundary values edgelx = 0. 0 edgerx = real(nx) if (ipbc. eq. 2) then edgelx = 1. 0 edgerx = real(nx-1) endif

c find interpolation weights do 10 j = 1, nop nn = part(1, j) +. 5 dxp = part(1, j) - real(nn) nn = nn + 1 amx = qm*(. 75 - dxp*dxp) dxl = qmh*(. 5 - dxp)**2 dxp = qmh*(. 5 + dxp)**2 c deposit current vy = part(3, j) vz = part(4, j) cu(1, nn) = cu(1, nn) + vy*dxl cu(2, nn) = cu(2, nn) + vz*dxl cu(1, nn+1) = cu(1, nn+1) + vy*amx cu(2, nn+1) = cu(2, nn+1) + vz*amx cu(1, nn+2) = cu(1, nn+2) + vy*dxp cu(2, nn+2) = cu(2, nn+2) + vz*dxp c advance position half a time-step dx = part(1, j) + part(2, j)*dt

c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, j) part(2, j) = -part(2, j) endif c set new position part(1, j) = dx 10 continue return end subroutine GSJPOST 1(part, cu, qm, dt, nop, idimp, nxv, ipbc) c for 1 -2/2 d code, this subroutine calculates particle current density c using second-order spline interpolation c in addition, particle positions are advanced a half time-step

c scalar version using guard cells, integer conversion precalculation c 27 flops/particle, 10 loads, 7 stores c input: all, output: part, cu c current density is approximated by values at the nearest grid points c cu(i, n)=qci*(. 75 -dx**2) c cu(i, n+1)=. 5*qci*(. 5+dx)**2 c cu(i, n-1)=. 5*qci*(. 5 -dx)**2 c where n = nearest grid point and dx = x-n c and qci = qm*vi, where i = y, z c part(1, n) = position x of particle n c part(2, n) = x velocity of particle n c part(3, n) = y velocity of particle n c part(4, n) = z velocity of particle n c cu(i, j+1) = ith component of current density at grid point j c qm = charge on particle, in units of e c dt = time interval between successive calculations c nop = number of particles c idimp = size of phase space = 4 c nxv = second dimension of current array, must be >= nx+3 c ipbc = particle boundary condition = (0, 1, 2) =

c (none, 2 d periodic, 2 d reflecting) implicit none integer nop, idimp, nxv, ipbc real part, cu, qm, dt dimension part(idimp, nop), cu(2, nxv) c local data integer nnn, j, nn real dxn, qmh, edgelx, edgerx, dxp, amx, dxl, vy, vz, dx real dx 1, dx 2, dx 3 if (nop. lt. 1) return c begin first particle nnn = part(1, 1) +. 5 dxn = part(1, 1) - real(nnn) qmh =. 5*qm c set boundary values edgelx = 0. 0 edgerx = real(nx) if (ipbc. eq. 2) then edgelx = 1. 0 edgerx = real(nx-1) endif

c find interpolation weights do 10 j = 2, nop nn = nnn + 1 nnn = part(1, j) +. 5 dxp = dxn = part(1, j) - real(nnn) amx = qm*(. 75 - dxp*dxp) dxl = qmh*(. 5 - dxp)**2 dxp = qmh*(. 5 + dxp)**2 c deposit current vy = part(3, j-1) vz = part(4, j-1) dx 1 = cu(1, nn) + vy*dxl = cu(2, nn) + vz*dxl dx 2 = cu(1, nn+1) + vy*amx = cu(2, nn+1) + vz*amx dx 3 = cu(1, nn+2) + vy*dxp = cu(2, nn+2) + vz*dxp cu(1, nn) = dx 1 cu(2, nn) = dxl

cu(1, nn+1) = dx 2 cu(2, nn+1) = amx cu(1, nn+2) = dx 3 cu(2, nn+2) = dxp c advance position half a time-step dx = part(1, j-1) + part(2, j-1)*dt c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, j-1) part(2, j-1) = -part(2, j-1) endif c set new position part(1, j-1) = dx 10 continue

c deposit current for last particle nn = nnn + 1 amx = qm*(. 75 - dxn*dxn) dxl = qmh*(. 5 - dxn)**2 dxp = qmh*(. 5 + dxn)**2 c deposit current vy = part(3, nop) vz = part(4, nop) dx 1 = cu(1, nn) + vy*dxl = cu(2, nn) + vz*dxl dx 2 = cu(1, nn+1) + vy*amx = cu(2, nn+1) + vz*amx dx 3 = cu(1, nn+2) + vy*dxp = cu(2, nn+2) + vz*dxp cu(1, nn) = dx 1 cu(2, nn) = dxl cu(1, nn+1) = dx 2 cu(2, nn+1) = amx cu(1, nn+2) = dx 3 cu(2, nn+2) = dxp

c advance position half a time-step dx = part(1, nop) + part(2, nop)*dt c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, nop) part(2, nop) = -part(2, nop) endif c set new position part(1, nop) = dx return end SKIPPED

subroutine GBPUSH 13(part, fxyz, byz, omx, qbm, dtc, ek, idimp, nop, nx, n 1 xv, ipbc) c for 1 -2/2 d code, this subroutine updates particle co-ordinate and c velocities using leap-frog scheme in time and second-order spline c interpolation in space, with magnetic field. Using the Boris Mover. c scalar version using guard cells c 96 flops/particle, 1 divide, 19 loads, 4 stores c input: all, output: part, ek c velocity equations used are: c vx(t+dt/2) = rot(1)*(vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt) + c rot(2)*(vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt) + c rot(3)*(vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt) +. 5*(q/m)*fx(x(t))*dt) c vy(t+dt/2) = rot(4)*(vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt) + c rot(5)*(vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt) + c rot(6)*(vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt) +. 5*(q/m)*fy(x(t))*dt) c vz(t+dt/2) = rot(7)*(vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt) + c rot(8)*(vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt) + c rot(9)*(vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt) c where q/m is charge/mass, and the rotation matrix is given by: c rot(1) = (1 - (om*dt/2)**2 + 2*(omx*dt/2)**2)/(1 + (om*dt/2)**2) c rot(2) = 2*(omz*dt/2 + (omx*dt/2)*(omy*dt/2))/(1 + (om*dt/2)**2)

c rot(3) = 2*(-omy*dt/2 + (omx*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(4) = 2*(-omz*dt/2 + (omx*dt/2)*(omy*dt/2))/(1 + (om*dt/2)**2) c rot(5) = (1 - (om*dt/2)**2 + 2*(omy*dt/2)**2)/(1 + (om*dt/2)**2) c rot(6) = 2*(omx*dt/2 + (omy*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(7) = 2*(omy*dt/2 + (omx*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(8) = 2*(-omx*dt/2 + (omy*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(9) = (1 - (om*dt/2)**2 + 2*(omz*dt/2)**2)/(1 + (om*dt/2)**2) c and om**2 = omx**2 + omy**2 + omz**2 c the rotation matrix is determined by: c omy = (q/m)*by(x(t)), and omz = (q/m)*bz(x(t)), c position equations used are: c x(t+dt) = x(t) + vx(t+dt/2)*dt c fx(x(t)) is approximated by interpolation from the nearest grid points c fx(x) = (. 75 -dx**2)*fx(n)+. 5*(fx(n+1)*(. 5+dx)**2+fx(n-1)*(. 5 -dx)**2) c where n = nearest grid point and dx = x-n c similarly for fy(x), fz(x), by(x), bz(x) c part(1, n) = position x of particle n c part(2, n) = velocity vx of particle n c part(3, n) = velocity vy of particle n c part(4, n) = velocity vz of particle n

c fxyz(1, j+1, k+1) = x component of force/charge at grid (j, k) c fxyz(2, j+1, k+1) = y component of force/charge at grid (j, k) c fxyz(3, j+1, k+1) = z component of force/charge at grid (j, k) c that is, convolution of electric field over particle shape c byz(1, j+1, k+1) = y component of magnetic field at grid (j, k) c byz(2, j+1, k+1) = z component of magnetic field at grid (j, k) c that is, the convolution of magnetic field over particle shape c omx = magnetic field electron cyclotron frequency in x c qbm = particle charge/mass ratio c dt = time interval between successive calculations c dtc = time interval between successive co-ordinate calculations c kinetic energy/mass at time t is also calculated, using c ek =. 5*sum((vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt)**2 + c (vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt)**2 + c (vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt)**2) c idimp = size of phase space = 4 c nop = number of particles c nx = system length in x direction c nxv = second dimension of field arrays, must be >= nx+3 c ipbc = particle boundary condition = (0, 1, 2) =

c (none, 2 d periodic, 2 d reflecting) implicit none integer idimp, nop, nxv, ipbc real part, fxyz, byz, omx, qbm, dtc, ek dimension part(idimp, nop) dimension fxyz(3, nxv), byz(2, nxv) c local data integer j, nn real qtmh, edgelx, edgerx, dxp, amx, dxl, dx, dy, dz real ox, oy, oz, acx, acy, acz, omxt, omyt, omzt, omt, anorm real rot 1, rot 2, rot 3, rot 4, rot 5, rot 6, rot 7, rot 8, rot 9 double precision sum 1 qtmh =. 5*qbm*dt sum 1 = 0. 0 d 0 c set boundary values edgelx = 0. 0 edgerx = real(nx) if (ipbc. eq. 2) then edgelx = 1. 0 edgerx = real(nx-1) endif

do 10 j = 1, nop c find interpolation weights nn = part(1, j) +. 5 dxp = part(1, j) - real(nn) nn = nn + 1 amx =. 75 - dxp*dxp dxl =. 5*(. 5 - dxp)**2 dxp =. 5*(. 5 + dxp)**2 c find electric field dx = amx*fxyz(1, nn+1) + fxyz(1, nn)*dxl + fxyz(1, nn+2)*dxp dy = amx*fxyz(2, nn+1) + fxyz(2, nn)*dxl + fxyz(2, nn+2)*dxp dz = amx*fxyz(3, nn+1) + fxyz(3, nn)*dxl + fxyz(3, nn+2)*dxp c find magnetic field ox = omx oy = amx*byz(1, nn+1) + byz(1, nn)*dxl + byz(1, nn+2)*dxp oz = amx*byz(2, nn+1) + byz(2, nn)*dxl + byz(2, nn+2)*dxp c calculate half impulse dx = qtmh*dx dy = qtmh*dy dz = qtmh*dz

c half acceleration acx = part(2, j) + dx acy = part(3, j) + dy acz = part(4, j) + dz c time-centered kinetic energy sum 1 = sum 1 + (acx*acx + acy*acy + acz*acz) c calculate cyclotron frequency omxt = qtmh*ox omyt = qtmh*oy omzt = qtmh*oz c calculate rotation matrix omt = omxt*omxt + omyt*omyt + omzt*omzt anorm = 2. /(1. + omt) omt =. 5*(1. - omt) rot 4 = omxt*omyt rot 7 = omxt*omzt rot 8 = omyt*omzt rot 1 = omt + omxt*omxt rot 5 = omt + omyt*omyt rot 9 = omt + omzt*omzt

rot 2 = omzt + rot 4 = -omzt + rot 4 rot 3 = -omyt + rot 7 = omyt + rot 7 rot 6 = omxt + rot 8 = -omxt + rot 8 c new velocity dx = (rot 1*acx + rot 2*acy + rot 3*acz)*anorm + dx dy = (rot 4*acx + rot 5*acy + rot 6*acz)*anorm + dy dz = (rot 7*acx + rot 8*acy + rot 9*acz)*anorm + dz part(2, j) = dx part(3, j) = dy part(4, j) = dz c new position dx = part(1, j) + dx*dtc c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx

c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, j) part(2, j) = -part(2, j) endif c set new position part(1, j) = dx 10 continue c normalize kinetic energy ek = ek +. 5*sum 1 return end subroutine GSBPUSH 13(part, fxyz, byz, omx, qbm, dtc, ek, idimp, nop, nx, 1 nxv, ipbc) c for 1 -2/2 d code, this subroutine updates particle co-ordinate and c velocities using leap-frog scheme in time and second-order spline c interpolation in space, with magnetic field. Using the Boris Mover. c scalar version using guard cells, integer conversion precalculation, c and 1 d addressing

c 96 flops/particle, 1 divide, 19 loads, 4 stores c input: all, output: part, ek c velocity equations used are: c vx(t+dt/2) = rot(1)*(vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt) + c rot(2)*(vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt) + c rot(3)*(vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt) +. 5*(q/m)*fx(x(t))*dt) c vy(t+dt/2) = rot(4)*(vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt) + c rot(5)*(vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt) + c rot(6)*(vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt) +. 5*(q/m)*fy(x(t))*dt) c vz(t+dt/2) = rot(7)*(vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt) + c rot(8)*(vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt) + c rot(9)*(vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt) c where q/m is charge/mass, and the rotation matrix is given by: c rot(1) = (1 - (om*dt/2)**2 + 2*(omx*dt/2)**2)/(1 + (om*dt/2)**2) c rot(2) = 2*(omz*dt/2 + (omx*dt/2)*(omy*dt/2))/(1 + (om*dt/2)**2) c rot(3) = 2*(-omy*dt/2 + (omx*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(4) = 2*(-omz*dt/2 + (omx*dt/2)*(omy*dt/2))/(1 + (om*dt/2)**2) c rot(5) = (1 - (om*dt/2)**2 + 2*(omy*dt/2)**2)/(1 + (om*dt/2)**2) c rot(6) = 2*(omx*dt/2 + (omy*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(7) = 2*(omy*dt/2 + (omx*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(8) = 2*(-omx*dt/2 + (omy*dt/2)*(omz*dt/2))/(1 + (om*dt/2)**2) c rot(9) = (1 - (om*dt/2)**2 + 2*(omz*dt/2)**2)/(1 + (om*dt/2)**2)

c and om**2 = omx**2 + omy**2 + omz**2 c the rotation matrix is determined by: c omy = (q/m)*by(x(t)), and omz = (q/m)*bz(x(t)), c position equations used are: c x(t+dt) = x(t) + vx(t+dt/2)*dt c fx(x(t)) is approximated by interpolation from the nearest grid points c fx(x) = (. 75 -dx**2)*fx(n)+. 5*(fx(n+1)*(. 5+dx)**2+fx(n-1)*(. 5 -dx)**2) c where n = nearest grid point and dx = x-n c similarly for fy(x), fz(x), by(x), bz(x) c part(1, n) = position x of particle n c part(2, n) = velocity vx of particle n c part(3, n) = velocity vy of particle n c part(4, n) = velocity vz of particle n c fxyz(1, j+1, k+1) = x component of force/charge at grid (j, k) c fxyz(2, j+1, k+1) = y component of force/charge at grid (j, k) c fxyz(3, j+1, k+1) = z component of force/charge at grid (j, k) c that is, convolution of electric field over particle shape c byz(1, j+1, k+1) = y component of magnetic field at grid (j, k) c byz(2, j+1, k+1) = z component of magnetic field at grid (j, k) c that is, the convolution of magnetic field over particle shape c omx = magnetic field electron cyclotron frequency in x

c qbm = particle charge/mass ratio c dt = time interval between successive calculations c dtc = time interval between successive co-ordinate calculations c kinetic energy/mass at time t is also calculated, using c ek =. 5*sum((vx(t-dt/2) +. 5*(q/m)*fx(x(t))*dt)**2 + c (vy(t-dt/2) +. 5*(q/m)*fy(x(t))*dt)**2 + c (vz(t-dt/2) +. 5*(q/m)*fz(x(t))*dt)**2) c idimp = size of phase space = 4 c nop = number of particles c nx = system length in x direction c nxv = second dimension of field arrays, must be >= nx+3 c ipbc = particle boundary condition = (0, 1, 2) = c (none, 2 d periodic, 2 d reflecting) implicit none integer idimp, nop, nxv, ipbc real part, fxyz, byz, omx, qbm, dtc, ek dimension part(idimp, nop) dimension fxyz(3, nxv), byz(2, nxv) c local data integer nnn, nop 1, j, nn real dxn, qtmh, edgelx, edgerx, amx, dxl, dxp, dx, dy, dz

real ox, oy, oz, acx, acy, acz, omxt, omyt, omzt, omt, anorm real rot 1, rot 2, rot 3, rot 4, rot 5, rot 6, rot 7, rot 8, rot 9 double precision sum 1 = 0. 0 d 0 if (nop. lt. 1) go to 20 c begin first particle nnn = part(1, 1) +. 5 dxn = part(1, 1) - real(nnn) nop 1 = nop - 1 qtmh =. 5*qbm*dt c set boundary values edgelx = 0. 0 edgerx = real(nx) if (ipbc. eq. 2) then edgelx = 1. 0 edgerx = real(nx-1) endif do 10 j = 1, nop 1 c find interpolation weights nn = nnn + 1 nnn = part(1, j+1) +. 5

dxp = dxn = part(1, j+1) - real(nnn) amx =. 75 - dxp*dxp dxl =. 5*(. 5 - dxp)**2 dxp =. 5*(. 5 + dxp)**2 c find electric field dx = amx*fxyz(1, nn+1) + fxyz(1, nn)*dxl + fxyz(1, nn+2)*dxp dy = amx*fxyz(2, nn+1) + fxyz(2, nn)*dxl + fxyz(2, nn+2)*dxp dz = amx*fxyz(3, nn+1) + fxyz(3, nn)*dxl + fxyz(3, nn+2)*dxp c find magnetic field ox = omx oy = amx*byz(1, nn+1) + byz(1, nn)*dxl + byz(1, nn+2)*dxp oz = amx*byz(2, nn+1) + byz(2, nn)*dxl + byz(2, nn+2)*dxp c calculate half impulse dx = qtmh*dx dy = qtmh*dy dz = qtmh*dz c half acceleration acx = part(2, j) + dx acy = part(3, j) + dy acz = part(4, j) + dz

c time-centered kinetic energy sum 1 = sum 1 + (acx*acx + acy*acy + acz*acz) c calculate cyclotron frequency omxt = qtmh*ox omyt = qtmh*oy omzt = qtmh*oz c calculate rotation matrix omt = omxt*omxt + omyt*omyt + omzt*omzt anorm = 2. /(1. + omt) omt =. 5*(1. - omt) rot 4 = omxt*omyt rot 7 = omxt*omzt rot 8 = omyt*omzt rot 1 = omt + omxt*omxt rot 5 = omt + omyt*omyt rot 9 = omt + omzt*omzt rot 2 = omzt + rot 4 = -omzt + rot 4 rot 3 = -omyt + rot 7 = omyt + rot 7 rot 6 = omxt + rot 8 = -omxt + rot 8

c new velocity dx = (rot 1*acx + rot 2*acy + rot 3*acz)*anorm + dx dy = (rot 4*acx + rot 5*acy + rot 6*acz)*anorm + dy dz = (rot 7*acx + rot 8*acy + rot 9*acz)*anorm + dz part(2, j) = dx part(3, j) = dy part(4, j) = dz c new position dx = part(1, j) + dx*dtc c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, j) part(2, j) = -part(2, j) endif

c set new position part(1, j) = dx 10 continue c push last particle nn = nnn + 1 amx =. 75 - dxn*dxn dxl =. 5*(. 5 - dxn)**2 dxp =. 5*(. 5 + dxn)**2 c find electric field dx = amx*fxyz(1, nn+1) + fxyz(1, nn)*dxl + fxyz(1, nn+2)*dxp dy = amx*fxyz(2, nn+1) + fxyz(2, nn)*dxl + fxyz(2, nn+2)*dxp dz = amx*fxyz(3, nn+1) + fxyz(3, nn)*dxl + fxyz(3, nn+2)*dxp c find magnetic field ox = omx oy = amx*byz(1, nn+1) + byz(1, nn)*dxl + byz(1, nn+2)*dxp oz = amx*byz(2, nn+1) + byz(2, nn)*dxl + byz(2, nn+2)*dxp c calculate half impulse dx = qtmh*dx dy = qtmh*dy dz = qtmh*dz

c half acceleration acx = part(2, nop) + dx acy = part(3, nop) + dy acz = part(4, nop) + dz c time-centered kinetic energy sum 1 = sum 1 + (acx*acx + acy*acy + acz*acz) c calculate cyclotron frequency omxt = qtmh*ox omyt = qtmh*oy omzt = qtmh*oz c calculate rotation matrix omt = omxt*omxt + omyt*omyt + omzt*omzt anorm = 2. /(1. + omt) omt =. 5*(1. - omt) rot 4 = omxt*omyt rot 7 = omxt*omzt rot 8 = omyt*omzt rot 1 = omt + omxt*omxt rot 5 = omt + omyt*omyt rot 9 = omt + omzt*omzt rot 2 = omzt + rot 4

rot 4 = -omzt + rot 4 rot 3 = -omyt + rot 7 = omyt + rot 7 rot 6 = omxt + rot 8 = -omxt + rot 8 c new velocity dx = (rot 1*acx + rot 2*acy + rot 3*acz)*anorm + dx dy = (rot 4*acx + rot 5*acy + rot 6*acz)*anorm + dy dz = (rot 7*acx + rot 8*acy + rot 9*acz)*anorm + dz part(2, nop) = dx part(3, nop) = dy part(4, nop) = dz c new position dx = part(1, nop) + dx*dtc c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx

c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, nop) part(2, nop) = -part(2, nop) endif c set new position part(1, nop) = dx c normalize kinetic energy 20 ek = ek +. 5*sum 1 return end SKIPPED

!-----------------------------------! module bpush 1 d ! ! Fortran 90 interface to 1 d PIC Fortran 77 library bpush 1 lib. f ! bpush 1 mod. f contains interface procedures to process particles with ! magnetic fields: ! defines module bpush 1 d ! djpost => igjpost 1 deposits current density, with various ! interpolations and optimizations. ! calls GJPOST 1, GSJPOST 1 X, GJPOST 1 L, GSJPOST 1 L, or ! GSJPOST 1 XL ! push 3 => igbpush 13 push particles with magnetic field, with various ! interpolations and optimizations. ! calls GBPUSH 13, GSBPUSH 13, GBPUSH 13 L, or GSBPUSH 13 L ! retard => iretard 1 retard particle position a half time-step. ! calls RETARD 1 ! written by viktor k. decyk, ucla ! copyright 1999, regents of the university of california ! update: july 12, 2011 !

use globals, only: LINEAR, QUADRATIC, STANDARD, LOOKAHEAD, VECTOR use diag 1 d, only: wtimer implicit none private public : : LINEAR, QUADRATIC, STANDARD, LOOKAHEAD, VECTOR public : : wtimer public : : djpost, push 3, retard ! ! define interface to original Fortran 77 procedures ! interface subroutine GJPOST 1(part, cu, qm, dt, nop, idimp, nxv, ipbc) implicit none integer : : nop, idimp, nxv, ipbc real : : qm, dt real, dimension(idimp, nop) : : part real, dimension(2, nxv) : : cu end subroutine end interface

interface subroutine GSJPOST 1(part, cu, qm, dt, nop, idimp, nxv, ipbc) implicit none integer : : nop, idimp, nxv, ipbc real : : qm, dt real, dimension(idimp, nop) : : part real, dimension(2, nxv) : : cu end subroutine end interface subroutine GSJPOST 1 X(part, cu, qm, dt, nop, idimp, nxv, ipbc) implicit none integer : : nop, idimp, nxv, ipbc real : : qm, dt real, dimension(idimp, nop) : : part real, dimension(2*nxv) : : cu end subroutine end interface

interface subroutine GJPOST 1 L(part, cu, qm, dt, nop, idimp, nxv, ipbc) implicit none integer : : nop, idimp, nxv, ipbc real : : qm, dt real, dimension(idimp, nop) : : part real, dimension(2, nxv) : : cu end subroutine end interface subroutine GSJPOST 1 L(part, cu, qm, dt, nop, idimp, nxv, ipbc) implicit none integer : : nop, idimp, nxv, ipbc real : : qm, dt real, dimension(idimp, nop) : : part real, dimension(2, nxv) : : cu end subroutine end interface

interface subroutine GSJPOST 1 XL(part, cu, qm, dt, nop, idimp, nxv, ipbc) implicit none integer : : nop, idimp, nxv, ipbc real : : qm, dt real, dimension(idimp, nop) : : part real, dimension(2*nxv) : : cu end subroutine end interface subroutine GBPUSH 13(part, fxyz, byz, omx, qbm, dtc, ek, idimp, nop, n& &x, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : omx, qbm, dtc, ek real, dimension(idimp, nop) : : part real, dimension(3, nxv) : : fxyz real, dimension(2, nxv) : : byz end subroutine end interface

interface subroutine GSBPUSH 13(part, fxyz, byz, omx, qbm, dtc, ek, idimp, nop, & &nx, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : omx, qbm, dtc, ek real, dimension(idimp, nop) : : part real, dimension(3, nxv) : : fxyz real, dimension(2, nxv) : : byz end subroutine end interface subroutine GBPUSH 13 L(part, fxyz, byz, omx, qbm, dtc, ek, idimp, nop, & &nx, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : omx, qbm, dtc, ek real, dimension(idimp, nop) : : part real, dimension(3, nxv) : : fxyz real, dimension(2, nxv) : : byz end subroutine end interface

interface subroutine GSBPUSH 13 L(part, fxyz, byz, omx, qbm, dtc, ek, idimp, nop& &, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : omx, qbm, dtc, ek real, dimension(idimp, nop) : : part real, dimension(3, nxv) : : fxyz real, dimension(2, nxv) : : byz end subroutine end interface subroutine RETARD 1(part, dtc, idimp, nop, nx, ipbc) implicit none integer : : idimp, nop, nx, ipbc real : : dtc real, dimension(idimp, nop) : : part end subroutine end interface

! ! define generic interface to Fortran 90 library ! interface djpost module procedure igjpost 1 end interface ! interface push 3 module procedure igbpush 13 end interface ! interface retard module procedure iretard 1 end interface ! ! define Fortran 90 interface functions to Fortran 77 library ! contains !

subroutine igjpost 1(part, cu, nop, qm, dt, tdjpost, nx, ipbc, inorder, d& &jopt) ! deposit current implicit none integer : : nop, nx, ipbc integer, optional : : inorder, djopt real : : qm, dt, tdjpost real, dimension(: , : ), pointer : : part real, dimension(: , : ), pointer : : cu ! local data integer : : idimp, nxv, order, opt, ltime real : : tj idimp = size(part, 1) nxv = size(cu, 2) order = QUADRATIC if (present(inorder)) order = inorder opt = STANDARD if (present(djopt)) opt = djopt

! initialize timer call wtimer(tj, ltime, -1) select case(size(cu, 1)) case (2) if (order==LINEAR) then if (opt==LOOKAHEAD) then call GSJPOST 1 L(part, cu, qm, dt, nop, idimp, nxv, ipbc) else if (opt==VECTOR) then call GSJPOST 1 XL(part, cu, qm, dt, nop, idimp, nxv, ipbc) else call GJPOST 1 L(part, cu, qm, dt, nop, idimp, nxv, ipbc) endif else if (opt==LOOKAHEAD) then call GSJPOST 1(part, cu, qm, dt, nop, idimp, nxv, ipbc) else if (opt==VECTOR) then call GSJPOST 1 X(part, cu, qm, dt, nop, idimp, nxv, ipbc) else call GJPOST 1(part, cu, qm, dt, nop, idimp, nxv, ipbc) endif end select

! record time call wtimer(tj, ltime) tdjpost = tdjpost + tj end subroutine igjpost 1 ! subroutine igbpush 13(part, fxyz, byz, omx, nop, qbm, dtc, ek, tpush, & &nx, ipbc, inorder, popt) ! push particles with 1 -2/2 d electromagnetic fields implicit none integer : : nop, nx, ipbc integer, optional : : inorder, popt real : : omx, qbm, dtc, ek, tpush real, dimension(: , : ), pointer : : part real, dimension(: , : ), pointer : : fxyz, byz ! local data integer : : idimp, nxv, order, opt, ltime real : : tp idimp = size(part, 1) nxv = size(fxyz, 2) order = QUADRATIC

if (present(inorder)) order = inorder opt = STANDARD if (present(popt)) opt = popt ! initialize timer call wtimer(tp, ltime, -1) select case(size(byz, 1)) case (2) if (order==LINEAR) then if (opt==LOOKAHEAD) then call GSBPUSH 13 L(part, fxyz, byz, omx, qbm, dtc, ek, idimp, & &nop, nxv, ipbc) else call GBPUSH 13 L(part, fxyz, byz, omx, qbm, dtc, ek, idimp, n& &op, nxv, ipbc) endif else if (opt==LOOKAHEAD) then call GSBPUSH 13(part, fxyz, byz, omx, qbm, dtc, ek, idimp, n& &op, nxv, ipbc)

else call GBPUSH 13(part, fxyz, byz, omx, qbm, dtc, ek, idimp, no& &p, nxv, ipbc) endif end select ! record time call wtimer(tp, ltime) tpush = tpush + tp end subroutine igbpush 13 ! subroutine iretard 1(part, nop, dtc, nx, ipbc) ! retards particle positions half time-step implicit none integer : : nop, nx integer, optional : : ipbc real : : dtc real, dimension(: , : ), pointer : : part

! local data integer : : idimp, lpbc idimp = size(part, 1) lpbc = 1 if (present(ipbc)) lpbc = ipbc call RETARD 1(part, dtc, idimp, nop, nx, ipbc) end subroutine iretard 1 ! end module bpush 1 d

These are for new_beps. f (electrostatic code) c-----------------------------------c 1 d PIC library for pushing particles and depositing charge c push 1 lib. f contains procedures to process particles: c GPOST 1 deposits charge density, quadratic interpolation, STANDARD c optimization. c GSPOST 1 deposits charge density, quadratic interpolation, LOOKAHEAD c optimization c GSPOST 1 X deposits charge density, quadratic interpolation, VECTOR c optimization. c GPOST 1 L deposits charge density, linear interpolation, STANDARD c optimization. c GSPOST 1 L deposits charge density, linear interpolation, LOOKAHEAD c optimization. c GSPOST 1 XL deposits charge density, linear interpolation, VECTOR c optimization. c GPUSH 1 push particles, quadratic interpolation, STANDARD optimization. c GSPUSH 1 push particles, quadratic interpolation, LOOKAHEAD c optimization.

c GPUSH 1 L push particles, linear interpolation, STANDARD optimization. c GSPUSH 1 L push particles, linear interpolation, LOOKAHEAD optimization. c SORTP 1 X sort particles by grid, quadratic interpolation, memory c conserving algorithm. c SORTP 1 XL sort particles by grid, linear interpolation, memory c conserving algorithm. c DSORTP 1 X sort particles by grid, quadratic interpolation, high c performance algorithm. c DSORTP 1 XL sort particles by grid, linear interpolation, high c performance algorithm. c RMOVE 1 remove particles instead of reflecting at boundary. c written by viktor k. decyk, ucla c copyright 1994, regents of the university of california c update: july 12, 2011

c-----------------------------------subroutine GPOST 1(part, q, qm, nop, idimp, nxv) attention c for 1 d code, this subroutine calculates particle charge density c using second-order spline interpolation, periodic boundaries c scalar version using guard cells c 16 flops/particle, 4 loads, 3 stores c input: all, output: q c charge density is approximated by values at the nearest grid points c q(n)=qm*(. 75 -dx**2), q(n+1)=. 5*qm*(. 5+dx)**2, q(n-1)=. 5*qm*(. 5 -dx)**2 c where n = nearest grid point and dx = x-n c part(1, n) = position x of particle n c q(j) = charge density at grid point j c qm = charge on particle, in units of e c nop = number of particles c idimp = size of phase space = 2 c nxv = first dimension of charge array, must be >= nx+3 implicit none integer nop, idimp, nxv real part, q, qm

dimension part(idimp, nop), q(nxv) c local data integer j, nn real qmh, dx qmh =. 5*qm c find interpolation weights do 10 j = 1, nop nn = part(1, j) +. 5 dx = part(1, j) - real(nn) nn = nn + 1 c deposit charge q(nn) = q(nn) + qmh*(. 5 - dx)**2 q(nn+1) = q(nn+1) + qm*(. 75 - dx*dx) q(nn+2) = q(nn+2) + qmh*(. 5 + dx)**2 10 continue return end

subroutine GSPOST 1(part, q, qm, nop, idimp, nxv) c for 1 d code, this subroutine calculates particle charge density c using second-order spline interpolation, periodic boundaries c scalar version using guard cells, integer conversion precalculation c 16 flops/particle, 4 loads, 3 stores c input: all, output: q c charge density is approximated by values at the nearest grid points c q(n)=qm*(. 75 -dx**2), q(n+1)=. 5*qm*(. 5+dx)**2, q(n-1)=. 5*qm*(. 5 -dx)**2 c where n = nearest grid point and dx = x-n c part(1, n) = position x of particle n c q(j) = charge density at grid point j c qm = charge on particle, in units of e c nop = number of particles c idimp = size of phase space = 2 c nxv = first dimension of charge array, must be >= nx+3 implicit none integer nop, idimp, nxv real part, q, qm dimension part(idimp, nop), q(nxv)

c local data integer nnn, j, nn real dxn, qmh, dx c begin first particle nnn = part(1, 1) +. 5 dxn = part(1, 1) - real(nnn) qmh =. 5*qm c find interpolation weights do 10 j = 2, nop nn = nnn + 1 nnn = part(1, j) +. 5 dx = dxn = part(1, j) - real(nnn) c deposit charge q(nn) = q(nn) + qmh*(. 5 - dx)**2 q(nn+1) = q(nn+1) + qm*(. 75 - dx*dx) q(nn+2) = q(nn+2) + qmh*(. 5 + dx)**2 10 continue

c deposit charge for last particle nn = nnn + 1 q(nn) = q(nn) + qmh*(. 5 - dxn)**2 q(nn+1) = q(nn+1) + qm*(. 75 - dxn*dxn) q(nn+2) = q(nn+2) + qmh*(. 5 + dxn)**2 return end subroutine GSPOST 1 X(part, q, qm, nop, idimp, nxv) c for 1 d code, this subroutine calculates particle charge density c using second-order spline interpolation, periodic boundaries, c with short vectors over independent weights, c as in j. schwartzmeier and t. hewitt, proc. 12 th conf. on numerical c simulation of plasmas, san francisco, ca, 1987. c vectorized version with guard cells c 16 flops/particle, 4 loads, 3 stores c input: all, output: q c charge density is approximated by values at the nearest grid points c q(n)=qm*(. 75 -dx**2), q(n+1)=. 5*qm*(. 5+dx)**2, q(n-1)=. 5*qm*(. 5 -dx)**2 c where n = nearest grid point and dx = x-n

c part(1, n) = position x of particle n c q(j) = charge density at grid point j c qm = charge on particle, in units of e c nop = number of particles c idimp = size of phase space = 2 c nxv = first dimension of charge array, must be >= nx+3 implicit none integer nop, idimp, nxv real part, q, qm dimension part(idimp, nop), q(nxv) c local data integer npp, nn, npb, ipp, k, j, je, jb, n, i real amx, qmh, dx parameter(npp=512) dimension nn(3, npp), amx(3, npp) npb = npp if (npb. gt. nop) npb = nop ipp = real(nop - 1)/real(npb) + 1. qmh =. 5*qm

c outer loop over blocks of particles do 40 k = 1, ipp je = k*npb jb = je - npb if (je. gt. nop) npb = nop - jb c find interpolation weights do 10 j = 1, npb n = part(1, j+jb) +. 5 dx = part(1, j+jb) - real(n) n=n+1 amx(1, j) = qmh*(. 5 - dx)**2 nn(1, j) = n amx(2, j) = qm*(. 75 - dx*dx) nn(2, j) = n + 1 amx(3, j) = qmh*(. 5 + dx)**2 nn(3, j) = n + 2 10 continue

c deposit charge do 30 j = 1, npb cdir$ ivdep do 20 i = 1, 3 q(nn(i, j)) = q(nn(i, j)) + amx(i, j) 20 continue 30 continue 40 continue return end subroutine GPOST 1 L(part, q, qm, nop, idimp, nxv) c for 1 d code, this subroutine calculates particle charge density c using first-order linear interpolation, periodic boundaries c scalar version using guard cells c 7 flops/particle, 3 loads, 3 stores c input: all, output: q c charge density is approximated by values at the nearest grid points c q(n)=qm*(1. -dx) and q(n+1)=qm*dx c where n = nearest grid point and dx = x-n

c part(1, n) = position x of particle n c q(j) = charge density at grid point j c qm = charge on particle, in units of e c nop = number of particles c idimp = size of phase space = 2 c nxv = first dimension of charge array, must be >= nx+1 implicit none integer nop, idimp, nxv real part, q, qm dimension part(idimp, nop), q(nxv) c local data integer j, nn real dx c find interpolation weights do 10 j = 1, nop nn = part(1, j) dx = qm*(part(1, j) - real(nn)) nn = nn + 1

c deposit charge q(nn) = q(nn) + (qm - dx) q(nn+1) = q(nn+1) + dx 10 continue return end subroutine GSPOST 1 L(part, q, qm, nop, idimp, nxv) c for 1 d code, this subroutine calculates particle charge density c using first-order linear interpolation, periodic boundaries c scalar version using guard cells, integer conversion precalculation c 7 flops/particle, 3 loads, 3 stores c input: all, output: q c charge density is approximated by values at the nearest grid points c q(n)=qm*(1. -dx) and q(n+1)=qm*dx c where n = nearest grid point and dx = x-n c part(1, n) = position x of particle n c q(j) = charge density at grid point j c qm = charge on particle, in units of e c nop = number of particles

c idimp = size of phase space = 2 c nxv = first dimension of charge array, must be >= nx+1 implicit none integer nop, idimp, nxv real part, q, qm dimension part(idimp, nop), q(nxv) c local data integer nnn, j, nn real dxn, dx c begin first particle nnn = part(1, 1) dxn = qm*(part(1, 1) - real(nnn)) c find interpolation weights do 10 j = 2, nop nn = nnn + 1 nnn = part(1, j) dx = dxn = qm*(part(1, j) - real(nnn))

c deposit charge q(nn) = q(nn) + (qm - dx) q(nn+1) = q(nn+1) + dx 10 continue c deposit charge for last particle nn = nnn + 1 q(nn) = q(nn) + (qm - dxn) q(nn+1) = q(nn+1) + dxn return end

subroutine GSPOST 1 XL(part, q, qm, nop, idimp, nxv) c for 1 d code, this subroutine calculates particle charge density c using first-order linear interpolation, periodic boundaries, c with short vectors over independent weights, c as in j. schwartzmeier and t. hewitt, proc. 12 th conf. on numerical c simulation of plasmas, san francisco, ca, 1987. c vectorized version with guard cells c 16 flops/particle, 4 loads, 3 stores c input: all, output: q c charge density is approximated by values at the nearest grid points c q(n)=qm*(1. -dx) and q(n+1)=qm*dx c where n = nearest grid point and dx = x-n c part(1, n) = position x of particle n c q(j) = charge density at grid point j c qm = charge on particle, in units of e c nop = number of particles c idimp = size of phase space = 2 c nxv = first dimension of charge array, must be >= nx+1 implicit none integer nop, idimp, nxv

real part, q, qm dimension part(idimp, nop), q(nxv) c local data integer npp, nn, npb, ipp, k, j, je, jb, n, i real amx, dx parameter(npp=512) dimension nn(2, npp), amx(2, npp) npb = npp if (npb. gt. nop) npb = nop ipp = real(nop - 1)/real(npb) + 1. c outer loop over blocks of particles do 40 k = 1, ipp je = k*npb jb = je - npb if (je. gt. nop) npb = nop - jb

c find interpolation weights do 10 j = 1, npb n = part(1, j+jb) dx = qm*(part(1, j+jb) - real(n)) n=n+1 amx(1, j) = qm - dx nn(1, j) = n amx(2, j) = dx nn(2, j) = n + 1 10 continue c deposit charge do 30 j = 1, npb cdir$ ivdep do 20 i = 1, 2 q(nn(i, j)) = q(nn(i, j)) + amx(i, j) 20 continue 30 continue 40 continue return end

subroutine GPUSH 1(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) c for 1 d code, this subroutine updates particle co-ordinate and velocity c using leap-frog scheme in time and second-order spline interpolation c in space, with periodic boundary conditions. c scalar version using guard cells c 25 flops/particle, 5 loads, 2 stores c input: all, output: part, ek c equations used are: c v(t+dt/2) = v(t-dt/2) + (q/m)*fx(x(t))*dt, where q/m is charge/mass, c and x(t+dt) = x(t) + v(t+dt/2)*dt c fx(x(t)) is approximated by interpolation from the nearest grid points c fx(x) = (. 75 -dx**2)*fx(n)+. 5*(fx(n+1)*(. 5+dx)**2+fx(n-1)*(. 5 -dx)**2) c where n = nearest grid point and dx = x-n c part(1, n) = position x of particle n c part(2, n) = velocity vx of particle n c fx(j) = force/charge at grid point j, that is convolution of electric c field over particle shape c qbm = particle charge/mass c dt = time interval between successive calculations

c kinetic energy/mass at time t is also calculated, using c ek =. 125*sum((v(t+dt/2)+v(t-dt/2))**2) c idimp = size of phase space = 2 c nop = number of particles c nx = system length in x direction c nxv = dimension of field array, must be >= nx+3 c ipbc = particle boundary condition = (0, 1, 2) = c (none, 2 d periodic, 2 d reflecting) implicit none integer idimp, nop, nxv, ipbc real part, fx, qbm, dt, ek dimension part(idimp, nop), fx(nxv) c local data integer j, nn real qtm, edgelx, edgerx, dx double precision sum 1 qtm = qbm*dt sum 1 = 0. 0 d 0

c set boundary values edgelx = 0. 0 edgerx = real(nx) if (ipbc. eq. 2) then edgelx = 1. 0 edgerx = real(nx-1) endif do 10 j = 1, nop c find interpolation weights nn = part(1, j) +. 5 dx = part(1, j) - real(nn) nn = nn + 1 c find acceleration dx = (. 75 - dx*dx)*fx(nn+1) +. 5*(fx(nn)*(. 5 - dx)**2 + fx(nn+2)*( 1. 5 + dx)**2) c new velocity dx = part(2, j) + qtm*dx c average kinetic energy sum 1 = sum 1 + (part(2, j) + dx)**2 part(2, j) = dx

c new position dx = part(1, j) + dx*dt c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, j) part(2, j) = -part(2, j) endif c set new position part(1, j) = dx 10 continue c normalize kinetic energy ek = ek +. 125*sum 1 return end c end of gpush 1

subroutine GSPUSH 1(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) c for 1 d code, this subroutine updates particle co-ordinate and velocity c using leap-frog scheme in time and second-order spline interpolation c in space, with periodic boundary conditions. c scalar version using guard cells, integer conversion precalculation c 25 flops/particle, 5 loads, 2 stores c input: all, output: part, ek c equations used are: c v(t+dt/2) = v(t-dt/2) + (q/m)*fx(x(t))*dt, where q/m is charge/mass, c and x(t+dt) = x(t) + v(t+dt/2)*dt c fx(x(t)) is approximated by interpolation from the nearest grid points c fx(x) = (. 75 -dx**2)*fx(n)+. 5*(fx(n+1)*(. 5+dx)**2+fx(n-1)*(. 5 -dx)**2) c where n = nearest grid point and dx = x-n c part(1, n) = position x of particle n c part(2, n) = velocity vx of particle n c fx(j) = force/charge at grid point j, that is convolution of electric c field over particle shape c qbm = particle charge/mass c dt = time interval between successive calculations

c kinetic energy/mass at time t is also calculated, using c ek =. 125*sum((v(t+dt/2)+v(t-dt/2))**2) c idimp = size of phase space = 2 c nop = number of particles c nx = system length in x direction c nxv = dimension of field array, must be >= nx+3 c ipbc = particle boundary condition = (0, 1, 2) = c (none, 2 d periodic, 2 d reflecting) implicit none integer idimp, nop, nxv, ipbc real part, fx, qbm, dt, ek dimension part(idimp, nop), fx(nxv) c local data integer nnn, nop 1, j, nn real dxn, qtm, edgelx, edgerx, dx double precision sum 1

c begin first particle nnn = part(1, 1) +. 5 dxn = part(1, 1) - real(nnn) nop 1 = nop - 1 qtm = qbm*dt sum 1 = 0. 0 d 0 c set boundary values edgelx = 0. 0 edgerx = real(nx) if (ipbc. eq. 2) then edgelx = 1. 0 edgerx = real(nx-1) endif do 10 j = 1, nop 1 c find interpolation weights nn = nnn + 1 nnn = part(1, j+1) +. 5 dx = dxn = part(1, j+1) - real(nnn)

c find acceleration dx = (. 75 - dx*dx)*fx(nn+1) +. 5*(fx(nn)*(. 5 - dx)**2 + fx(nn+2)*( 1. 5 + dx)**2) c new velocity dx = part(2, j) + qtm*dx c average kinetic energy sum 1 = sum 1 + (part(2, j) + dx)**2 part(2, j) = dx c new position dx = part(1, j) + dx*dt c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, j) part(2, j) = -part(2, j) endif

part(1, j) = dx 10 continue c push last particle nn = nnn + 1 c find acceleration dx = (. 75 - dxn*dxn)*fx(nn+1) +. 5*(fx(nn)*(. 5 - dxn)**2 + fx(nn+2 1)*(. 5 + dxn)**2) c new velocity dx = part(2, nop) + qtm*dx c average kinetic energy sum 1 = sum 1 + (part(2, nop) + dx)**2 part(2, nop) = dx c new position dx = part(1, nop) + dx*dt c periodic boundary conditions if (ipbc. eq. 1) then if (dx. lt. edgelx) dx = dx + edgerx if (dx. ge. edgerx) dx = dx - edgerx

c reflecting boundary conditions else if (ipbc. eq. 2) then if ((dx. lt. edgelx). or. (dx. ge. edgerx)) then dx = part(1, nop) part(2, nop) = -part(2, nop) endif part(1, nop) = dx c normalize kinetic energy ek = ek +. 125*sum 1 return end continued!

!-----------------------------------! module push 1 d ! ! Fortran 90 interface to 1 d PIC Fortran 77 library push 1 lib. f ! push 1 mod. f contains interface procedures to process particles: ! defines module push 1 d ! dpost => igpost 1 deposits charge density, with various interpolations ! and optimizations. ! calls GPOST 1, GSPOST 1 X, GPOST 1 L, GSPOST 1 L, or ! GSPOST 1 XL ! push => igpush 1 pushes particles, with various interpolations and ! optimizations. attention ! calls GPUSH 1, GSPUSH 1, GPUSH 1 L, or GSPUSH 1 L ! sortp => isortp 1 x sort particles by grid with various interpolations. ! calls SORTP 1 X, or SORTP 1 XL ! written by viktor k. decyk, ucla ! copyright 1999, regents of the university of california ! update: july 12, 2011 !

use globals, only: LINEAR, QUADRATIC, STANDARD, LOOKAHEAD, VECTO use diag 1 d, only: wtimer implicit none private public : : LINEAR, QUADRATIC, STANDARD, LOOKAHEAD, VECTOR public : : wtimer public : : dpost, push, sortp, dpostgl, pushgl ! ! define interface to original Fortran 77 procedures ! interface subroutine GPOST 1(part, q, qm, nop, idimp, nxv) implicit none integer : : nop, idimp, nxv real : : qm real, dimension(idimp, nop) : : part real, dimension(nxv) : : q end subroutine end interface

interface subroutine GSPOST 1(part, q, qm, nop, idimp, nxv) implicit none integer : : nop, idimp, nxv real : : qm real, dimension(idimp, nop) : : part real, dimension(nxv) : : q end subroutine end interface subroutine GSPOST 1 X(part, q, qm, nop, idimp, nxv) implicit none integer : : nop, idimp, nxv real : : qm real, dimension(idimp, nop) : : part real, dimension(nxv) : : q end subroutine end interface

interface subroutine GPOST 1 L(part, q, qm, nop, idimp, nxv) implicit none integer : : nop, idimp, nxv real : : qm real, dimension(idimp, nop) : : part real, dimension(nxv) : : q end subroutine end interface subroutine GSPOST 1 L(part, q, qm, nop, idimp, nxv) implicit none integer : : nop, idimp, nxv real : : qm real, dimension(idimp, nop) : : part real, dimension(nxv) : : q end subroutine end interface

interface subroutine GSPOST 1 XL(part, q, qm, nop, idimp, nxv) implicit none integer : : nop, idimp, nxv real : : qm real, dimension(idimp, nop) : : part real, dimension(nxv) : : q end subroutine end interface subroutine GPUSH 1(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : qbm, dt, ek real, dimension(idimp, nop) : : part real, dimension(nxv) : : fx end subroutine end interface

interface subroutine GSPUSH 1(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : qbm, dt, ek real, dimension(idimp, nop) : : part real, dimension(nxv) : : fx end subroutine end interface subroutine GPUSH 1 L(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : qbm, dt, ek real, dimension(idimp, nop) : : part real, dimension(nxv) : : fx end subroutine end interface

interface subroutine GSPUSH 1 L(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) implicit none integer : : idimp, nop, nxv, ipbc real : : qbm, dt, ek real, dimension(idimp, nop) : : part real, dimension(nxv) : : fx end subroutine end interface subroutine SORTP 1 X(part, pt, ip, npic, idimp, nop, nx 1) implicit none integer : : idimp, nop, nx 1 real, dimension(idimp, nop) : : part real, dimension(nop) : : pt integer, dimension(nop) : : ip integer, dimension(nx 1) : : npic end subroutine end interface

interface subroutine SORTP 1 XL(part, pt, ip, npic, idimp, nop, nx 1) implicit none integer : : idimp, nop, nx 1 real, dimension(idimp, nop) : : part real, dimension(nop) : : pt integer, dimension(nop) : : ip integer, dimension(nx 1) : : npic end subroutine end interface subroutine DPOST 1 GL(part, q, sctx, qm, nop, idimp, nxh, nxvh) implicit none integer : : nop, idimp, nxh, nxvh real : : qm real, dimension(idimp, nop) : : part real, dimension(2*nxvh) : : q complex, dimension(nxvh) : : sctx end subroutine end interface

interface subroutine PUSH 1 GL(part, fx, sctx, qbm, dt, ek, idimp, nop, nxh, nxvh& &) implicit none integer : : nop, idimp, nxh, nxvh real : : qbm, dt, ek real, dimension(idimp, nop) : : part real, dimension(2*nxvh) : : fx complex, dimension(nxvh) : : sctx end subroutine end interface ! ! define generic interface to Fortran 90 library ! interface dpost module procedure igpost 1 end interface ! interface push module procedure igpush 1 end interface

! interface sortp module procedure isortp 1 x end interface ! interface dpostgl module procedure idpost 1 gl end interface ! interface pushgl module procedure ipush 1 gl end interface

! ! define Fortran 90 interface functions to Fortran 77 library ! contains ! subroutine igpost 1(part, q, nop, qm, tdpost, inorder, dopt) ! deposit charge implicit none integer : : nop integer, optional : : inorder, dopt real : : qm, tdpost real, dimension(: , : ), pointer : : part real, dimension(: ), pointer : : q ! local data integer : : idimp, nxv, order, opt, ltime real : : td idimp = size(part, 1) nxv = size(q) order = QUADRATIC

if (present(inorder)) order = inorder opt = STANDARD if (present(dopt)) opt = dopt ! initialize timer call wtimer(td, ltime, -1) if (order==LINEAR) then if (opt==LOOKAHEAD) then call GSPOST 1 L(part, q, qm, nop, idimp, nxv) else if (opt==VECTOR) then call GSPOST 1 XL(part, q, qm, nop, idimp, nxv) else call GPOST 1 L(part, q, qm, nop, idimp, nxv) endif else if (opt==LOOKAHEAD) then call GSPOST 1(part, q, qm, nop, idimp, nxv) else if (opt==VECTOR) then call GSPOST 1 X(part, q, qm, nop, idimp, nxv) else call GPOST 1(part, q, qm, nop, idimp, nxv) endif

! record time call wtimer(td, ltime) tdpost = tdpost + td end subroutine igpost 1 ! subroutine igpush 1(part, fx, nop, qbm, dt, ek, tpush, nx, ipbc, inorder, & &popt) ! push particles with 1 d electrostatic fields implicit none integer : : nop, nx, ipbc integer, optional : : inorder, popt real : : qbm, dt, ek, tpush real, dimension(: , : ), pointer : : part real, dimension(: ), pointer : : fx ! local data integer : : idimp, nxv, order, opt, ltime real : : tp idimp = size(part, 1) nxv = size(fx) order = QUADRATIC

if (present(inorder)) order = inorder opt = STANDARD if (present(popt)) opt = popt ! initialize timer call wtimer(tp, ltime, -1) if (order==LINEAR) then if (opt==LOOKAHEAD) then call GSPUSH 1 L(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) else call GPUSH 1 L(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) endif else if (opt==LOOKAHEAD) then call GSPUSH 1(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) else call GPUSH 1(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) endif

! record time call wtimer(tp, ltime) tpush = tpush + tp end subroutine igpush 1 ! subroutine isortp 1 x(part, pt, ip, nop, npic, tsort, inorder) ! sort particles by x grid using memory-conserving bin sort implicit none integer : : nop integer, optional : : inorder real : : tsort real, dimension(: , : ), pointer : : part real, dimension(: ), pointer : : pt integer, dimension(: ), pointer : : ip integer, dimension(: ), pointer : : npic ! local data integer : : idimp, nx 1, order, ltime real : : ts order = QUADRATIC if (present(inorder)) order = inorder idimp = size(part, 1); nx 1 = size(npic)

! initialize timer call wtimer(ts, ltime, -1) if (order==LINEAR) then call SORTP 1 XL(part, pt, ip, npic, idimp, nop, nx 1) else call SORTP 1 X(part, pt, ip, npic, idimp, nop, nx 1) endif ! record time call wtimer(ts, ltime) tsort = tsort + ts end subroutine isortp 1 x

! subroutine idpost 1 gl(part, q, nop, qm, nxh, tdpost) ! deposit charge using gridless method implicit none integer : : nop, nxh real : : qm, tdpost real, dimension(: , : ), pointer : : part real, dimension(: ), pointer : : q ! local data integer : : idimp, nxvh, ltime real : : td complex, dimension(size(q, 1)/2) : : sctx idimp = size(part, 1) nxvh = size(q, 1)/2 ! initialize timer call wtimer(td, ltime, -1) call DPOST 1 GL(part, q, sctx, qm, nop, idimp, nxh, nxvh) ! record time call wtimer(td, ltime) tdpost = tdpost + td end subroutine idpost 1 gl

! subroutine ipush 1 gl(part, fx, nop, qbm, dt, ek, nxh, tpush) ! push particles with 1 d electrostatic fields using gridless method implicit none integer : : nop, nxh real : : qbm, dt, ek, tpush real, dimension(: , : ), pointer : : part real, dimension(: ), pointer : : fx ! local data integer : : idimp, nxvh, ltime real : : tp complex, dimension(size(fx, 1)/2) : : sctx idimp = size(part, 1) nxvh = size(fx, 1)/2 ! initialize timer call wtimer(tp, ltime, -1) call PUSH 1 GL(part, fx, sctx, qbm, dt, ek, idimp, nop, nxh, nxvh) ! record time call wtimer(tp, ltime) tpush = tpush + tp end subroutine ipush 1 gl ! end module push 1 d

Summary in new_bbeps 1. f ! push particles attention wke = 0. call push 3(part, fxyze, byze, omx, np, qbme, dth, wke, tpush, nx, ipbc, & &inorder, popt) in push 1 lib. f for example subroutine GPOST 1(part, q, qm, nop, idimp, nxv) subroutine GPUSH 1(part, fx, qbm, dt, ek, idimp, nop, nxv, ipbc) subroutine SORTP 1 X(part, pt, ip, npic, idimp, nop, nx 1) in push 1 mod. f ! subroutine igpush 1(part, fx, nop, qbm, dt, ek, tpush, nx, ipbc, inorder, & &popt)