! OpenSC: a package for computing 3D space charge fields using convolution-based methods with integrated Green functions (IGFs). ! Boundary conditions in this version: (1) free-space, (2) rectangular pipe ! R.D. Ryne, February 2018 module open_spacecharge_core_mod use, intrinsic :: iso_fortran_env !uncomment the following for the version that uses the 1D FFT from Alan Miller's web page !and compile with: gfortran test_mod.f90 fast_fourier_am.f90 open_spacecharge_mod.f90 test_opensc.f90 !use fast_fourier_am !uncomment the following for the version that uses FFTW: use fft_interface_mod !$ use omp_lib implicit none ! Fortran 2008 integer, parameter, private :: sp = REAL32 integer, parameter, private :: dp = REAL64 integer, parameter, private :: qp = REAL128 !the following arrays are Green function arrays used by the rectangular pipe solver complex(dp), allocatable, private, dimension(:,:,:) :: cgrn1,cgrn2,cgrn3,cgrn4 integer, private :: g1ilo,g1jlo,g1klo,g2ilo,g2jlo,g2klo,g3ilo,g3jlo,g3klo,g4ilo,g4jlo,g4klo !lower bounds of grn1234 arrays contains !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ subroutine osc_freespace_solver(rho,gam0,delta,phi,efield,bfield,nlo,nhi,nlo_gbl,nhi_gbl,npad,idirectfieldcalc,igfflag) !note well: this assumes that osc_alloc_freespace_array has been called so that cgrn1 is already allocated! !use mpi implicit none integer, intent(in) :: igfflag,idirectfieldcalc real(dp), intent(in), dimension(3) :: delta integer, intent(in), dimension(3) :: nlo,nhi,nlo_gbl,nhi_gbl,npad real(dp), intent(in) :: gam0 real(dp), intent(in), dimension(:,:,:) :: rho real(dp), intent(out), dimension(:,:,:) :: phi real(dp), intent(out), dimension(:,:,:,:) :: efield,bfield real(dp) :: dx,dy,dz integer :: ilo,ihi,jlo,jhi,klo,khi integer :: ilo2,ihi2,jlo2,jhi2,klo2,khi2,iperiod,jperiod,kperiod complex(dp), allocatable, dimension(:,:,:) :: crho integer :: icomp,i,j,k,im1,ip1,jm1,jp1,km1,kp1 real(dp) :: gb0,xfac,yfac,zfac integer :: mprocs,myrank,ierr real(dp), parameter :: clight=299792458.d0 ! !call MPI_COMM_SIZE(MPI_COMM_WORLD,mprocs,ierr) !call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) myrank=0 if (igfflag == 0 .and. idirectfieldcalc.eq.1) then print *, 'Error: Direct field calc must use integrated Green function' print *, 'Aborting...' return endif dx=delta(1); dy=delta(2); dz=delta(3) ilo=nlo(1);ihi=nhi(1);jlo=nlo(2);jhi=nhi(2);klo=nlo(3);khi=nhi(3) ilo2=ilo; jlo2=jlo; klo2=klo if(.not.allocated(cgrn1))call osc_alloc_freespace_array(nlo,nhi,npad) iperiod=size(cgrn1,1); jperiod=size(cgrn1,2); kperiod=size(cgrn1,3) ihi2=ilo2+iperiod-1; jhi2=jlo2+jperiod-1; khi2=klo2+kperiod-1 allocate(crho(ilo2:ihi2,jlo2:jhi2,klo2:khi2)) !double-size complex array for the charge density write(6,*)'iperiod,jperiod,kperiod=',iperiod,jperiod,kperiod if(idirectfieldcalc.eq.0)then if(myrank.eq.0)write(6,*)'Solving for Phi' icomp=0 !compute/store the FFT of the charge density: call getrhotilde(rho,crho,ilo,jlo,klo) !compute/store the FFT of the green function: call osc_getgrnfree(cgrn1,delta,gam0,g1ilo,g1jlo,g1klo,npad,icomp,igfflag) !compute the convolution: (note that crho,cgrn1 are complex padded, phi is real not padded) call conv3d(crho,cgrn1,phi,ilo,jlo,klo,g1ilo,g1jlo,g1klo,ilo,jlo,klo,iperiod,jperiod,kperiod) !original version did not handle the edge values correctly when differencing; here is a first order approx: do k=lbound(phi,3),ubound(phi,3) km1=k-1; if(km1.lt.lbound(phi,3))km1=k kp1=k+1; if(kp1.gt.ubound(phi,3))kp1=k zfac=merge(2.d0,1.d0,(km1.eq.k).or.(kp1.eq.k)) do j=lbound(phi,2),ubound(phi,2) jm1=j-1; if(jm1.lt.lbound(phi,2))jm1=j jp1=j+1; if(jp1.gt.ubound(phi,2))jp1=j yfac=merge(2.d0,1.d0,(jm1.eq.j).or.(jp1.eq.j)) do i=lbound(phi,1),ubound(phi,1) im1=i-1; if(im1.lt.lbound(phi,1))im1=i ip1=i+1; if(ip1.gt.ubound(phi,1))ip1=i xfac=merge(2.d0,1.d0,(im1.eq.i).or.(ip1.eq.i)) efield(i,j,k,1)=-(phi(ip1,j,k)-phi(im1,j,k))/(2.d0*dx)*gam0*xfac efield(i,j,k,2)=-(phi(i,jp1,k)-phi(i,jm1,k))/(2.d0*dy)*gam0*yfac efield(i,j,k,3)=-(phi(i,j,kp1)-phi(i,j,km1))/(2.d0*dz)/gam0*zfac enddo enddo enddo endif if(idirectfieldcalc.eq.1)then if(myrank.eq.0)write(6,*)'Solving for Ex,Ey,Ez' if(.not.allocated(crho))allocate(crho(ilo2:ihi2,jlo2:jhi2,klo2:khi2)) !double-size array call getrhotilde(rho,crho,ilo,jlo,klo) if(.not.allocated(cgrn1))then write(6,*)'something wrong. freespace green function should have been allocated during initialization' stop endif icomp=1 call osc_getgrnfree(cgrn1,delta,gam0,g1ilo,g1jlo,g1klo,npad,icomp,igfflag) call conv3d(crho,cgrn1,efield(:,:,:,1),ilo,jlo,klo,g1ilo,g1jlo,g1klo,ilo,jlo,klo,iperiod,jperiod,kperiod) icomp=2 call osc_getgrnfree(cgrn1,delta,gam0,g1ilo,g1jlo,g1klo,npad,icomp,igfflag) call conv3d(crho,cgrn1,efield(:,:,:,2),ilo,jlo,klo,g1ilo,g1jlo,g1klo,ilo,jlo,klo,iperiod,jperiod,kperiod) icomp=3 call osc_getgrnfree(cgrn1,delta,gam0,g1ilo,g1jlo,g1klo,npad,icomp,igfflag) call conv3d(crho,cgrn1,efield(:,:,:,3),ilo,jlo,klo,g1ilo,g1jlo,g1klo,ilo,jlo,klo,iperiod,jperiod,kperiod) endif ! Set the magnetic field: ! Zero field if at rest gb0=sqrt((gam0+1.d0)*(gam0-1.d0)) if (gb0 == 0) then bfield = 0 return endif do k=lbound(phi,3),ubound(phi,3) do j=lbound(phi,2),ubound(phi,2) do i=lbound(phi,1),ubound(phi,1) bfield(i,j,k,1)=-efield(i,j,k,2)/clight/gb0/gam0 bfield(i,j,k,2)= efield(i,j,k,1)/clight/gb0/gam0 bfield(i,j,k,3)=0.d0 enddo enddo enddo end subroutine osc_freespace_solver !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine osc_alloc_freespace_array(nlo,nhi,npad) implicit none integer, intent(in), dimension(3) :: nlo,nhi,npad integer :: rilo,rihi,rjlo,rjhi,rklo,rkhi !rho dimensions integer :: cilo,cihi,cjlo,cjhi,cklo,ckhi !phi dimensions integer :: ipad,jpad,kpad integer :: g1ihi,g1jhi,g1khi rilo=nlo(1); rihi=nhi(1); rjlo=nlo(2); rjhi=nhi(2); rklo=nlo(3); rkhi=nhi(3) cilo=nlo(1); cihi=nhi(1); cjlo=nlo(2); cjhi=nhi(2); cklo=nlo(3); ckhi=nhi(3) ipad=npad(1); jpad=npad(2); kpad=npad(3) g1ilo=cilo-rihi; g1ihi=cihi-rilo; g1jlo=cjlo-rjhi; g1jhi=cjhi-rjlo; g1klo=cklo-rkhi; g1khi=ckhi-rklo if(.not.allocated(cgrn1))allocate(cgrn1(g1ilo:g1ihi+ipad,g1jlo:g1jhi+jpad,g1klo:g1khi+kpad)) return end subroutine osc_alloc_freespace_array !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function coulombfun(u,v,w,gam) result(res) implicit none real(dp) :: res real(dp) :: u,v,w,gam if(u.eq.0.d0 .and. v.eq.0.d0 .and. w.eq.0.d0)then res=0.d0 return endif res=1.d0/sqrt(u**2+v**2+(gam*w)**2) !coulomb ! res=u/(u**2+v**2+(gam*w)**2)**1.5d0 !x-electric field return end function coulombfun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function igfcoulombfun(u,v,w,gam,dx,dy,dz) result(res) implicit none real(dp) :: res real(dp) :: u,v,w,gam,dx,dy,dz real(dp) :: x1,x2,y1,y2,z1,z2 real(dp) :: x,y,z !-! real(dp), external :: lafun x1=u-0.5d0*dx x2=u+0.5d0*dx y1=v-0.5d0*dy y2=v+0.5d0*dy z1=(w-0.5d0*dz)*gam z2=(w+0.5d0*dz)*gam ! res=1.d0/sqrt(u**2+v**2+w**2) !coulomb ! res=u/(u**2+v**2+w**2)**1.5d0 !x-electric field res=lafun(x2,y2,z2)-lafun(x1,y2,z2)-lafun(x2,y1,z2)-lafun(x2,y2,z1)-lafun(x1,y1,z1)+ & lafun(x1,y1,z2)+lafun(x1,y2,z1)+lafun(x2,y1,z1) res=res/(dx*dy*dz*gam) end function igfcoulombfun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function lafun(x,y,z) result(res) ! lafun is the function involving log and atan in the PRSTAB paper (I should find a better name for this function) implicit none real(dp) :: res real(dp) :: x,y,z,r r=sqrt(x**2+y**2+z**2) res=-0.5d0*z**2*atan(x*y/(z*r))-0.5d0*y**2*atan(x*z/(y*r))-0.5d0*x**2*atan(y*z/(x*r)) & +y*z*log(x+r)+x*z*log(y+r)+x*y*log(z+r) end function lafun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function igfexfun(u,v,w,gam,dx,dy,dz) result(res) implicit none real(dp) :: res real(dp) :: u,v,w,gam,dx,dy,dz real(dp) :: x1,x2,y1,y2,z1,z2 real(dp) :: x,y,z real(dp), parameter :: ep=1.d-10 !-13 causes NaN's real(dp), parameter :: em=-1.d-10 !-! real(dp), external :: xlafun x1=u-0.5d0*dx x2=u+0.5d0*dx y1=v-0.5d0*dy y2=v+0.5d0*dy z1=(w-0.5d0*dz)*gam z2=(w+0.5d0*dz)*gam if(x1<0.d0.and.x2>0.d0.and.y1<0.d0.and.y2>0.d0.and.z1<0.d0.and.z2>0.d0)then res=xlafun(em,em,em) & -xlafun(x1,em,em)-xlafun(em,y1,em)-xlafun(em,em,z1)-xlafun(x1,y1,z1) & +xlafun(x1,y1,em)+xlafun(x1,em,z1)+xlafun(em,y1,z1)+xlafun(x2,em,em) & -xlafun(ep,em,em)-xlafun(x2,y1,em)-xlafun(x2,em,z1)-xlafun(ep,y1,z1) & +xlafun(ep,y1,em)+xlafun(ep,em,z1)+xlafun(x2,y1,z1)+xlafun(ep,y2,ep) & -xlafun(x1,y2,ep)-xlafun(ep,ep,ep)-xlafun(ep,y2,z1)-xlafun(x1,ep,z1) & +xlafun(x1,ep,ep)+xlafun(x1,y2,z1)+xlafun(ep,ep,z1)+xlafun(ep,ep,z2) & -xlafun(x1,ep,z2)-xlafun(ep,y1,z2)-xlafun(ep,ep,ep)-xlafun(x1,y1,ep) & +xlafun(x1,y1,z2)+xlafun(x1,ep,ep)+xlafun(ep,y1,ep)+xlafun(x2,y2,ep) & -xlafun(ep,y2,ep)-xlafun(x2,ep,ep)-xlafun(x2,y2,z1)-xlafun(ep,ep,z1) & +xlafun(ep,ep,ep)+xlafun(ep,y2,z1)+xlafun(x2,ep,z1)+xlafun(x2,ep,z2) & -xlafun(ep,ep,z2)-xlafun(x2,y1,z2)-xlafun(x2,ep,ep)-xlafun(ep,y1,ep) & +xlafun(ep,y1,z2)+xlafun(ep,ep,ep)+xlafun(x2,y1,ep)+xlafun(ep,y2,z2) & -xlafun(x1,y2,z2)-xlafun(ep,ep,z2)-xlafun(ep,y2,ep)-xlafun(x1,ep,ep) & +xlafun(x1,ep,z2)+xlafun(x1,y2,ep)+xlafun(ep,ep,ep)+xlafun(x2,y2,z2) & -xlafun(ep,y2,z2)-xlafun(x2,ep,z2)-xlafun(x2,y2,ep)-xlafun(ep,ep,ep) & +xlafun(ep,ep,z2)+xlafun(ep,y2,ep)+xlafun(x2,ep,ep) else res=xlafun(x2,y2,z2)-xlafun(x1,y2,z2)-xlafun(x2,y1,z2)-xlafun(x2,y2,z1) & -xlafun(x1,y1,z1)+xlafun(x1,y1,z2)+xlafun(x1,y2,z1)+xlafun(x2,y1,z1) endif res=res/(dx*dy*dz) return end function igfexfun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function igfeyfun(u,v,w,gam,dx,dy,dz) result(res) implicit none real(dp) :: res real(dp) :: u,v,w,gam,dx,dy,dz real(dp) :: x1,x2,y1,y2,z1,z2 real(dp) :: x,y,z real(dp), parameter :: ep=1.d-10 !should include em also, but probably doesn't matter !-! real(dp), external :: ylafun x1=u-0.5d0*dx x2=u+0.5d0*dx y1=v-0.5d0*dy y2=v+0.5d0*dy z1=(w-0.5d0*dz)*gam z2=(w+0.5d0*dz)*gam if(x1<0.d0.and.x2>0.d0.and.y1<0.d0.and.y2>0.d0.and.z1<0.d0.and.z2>0.d0)then res=ylafun(ep,ep,ep) & -ylafun(x1,ep,ep)-ylafun(ep,y1,ep)-ylafun(ep,ep,z1)-ylafun(x1,y1,z1) & +ylafun(x1,y1,ep)+ylafun(x1,ep,z1)+ylafun(ep,y1,z1)+ylafun(x2,ep,ep) & -ylafun(ep,ep,ep)-ylafun(x2,y1,ep)-ylafun(x2,ep,z1)-ylafun(ep,y1,z1) & +ylafun(ep,y1,ep)+ylafun(ep,ep,z1)+ylafun(x2,y1,z1)+ylafun(ep,y2,ep) & -ylafun(x1,y2,ep)-ylafun(ep,ep,ep)-ylafun(ep,y2,z1)-ylafun(x1,ep,z1) & +ylafun(x1,ep,ep)+ylafun(x1,y2,z1)+ylafun(ep,ep,z1)+ylafun(ep,ep,z2) & -ylafun(x1,ep,z2)-ylafun(ep,y1,z2)-ylafun(ep,ep,ep)-ylafun(x1,y1,ep) & +ylafun(x1,y1,z2)+ylafun(x1,ep,ep)+ylafun(ep,y1,ep)+ylafun(x2,y2,ep) & -ylafun(ep,y2,ep)-ylafun(x2,ep,ep)-ylafun(x2,y2,z1)-ylafun(ep,ep,z1) & +ylafun(ep,ep,ep)+ylafun(ep,y2,z1)+ylafun(x2,ep,z1)+ylafun(x2,ep,z2) & -ylafun(ep,ep,z2)-ylafun(x2,y1,z2)-ylafun(x2,ep,ep)-ylafun(ep,y1,ep) & +ylafun(ep,y1,z2)+ylafun(ep,ep,ep)+ylafun(x2,y1,ep)+ylafun(ep,y2,z2) & -ylafun(x1,y2,z2)-ylafun(ep,ep,z2)-ylafun(ep,y2,ep)-ylafun(x1,ep,ep) & +ylafun(x1,ep,z2)+ylafun(x1,y2,ep)+ylafun(ep,ep,ep)+ylafun(x2,y2,z2) & -ylafun(ep,y2,z2)-ylafun(x2,ep,z2)-ylafun(x2,y2,ep)-ylafun(ep,ep,ep) & +ylafun(ep,ep,z2)+ylafun(ep,y2,ep)+ylafun(x2,ep,ep) else res=ylafun(x2,y2,z2)-ylafun(x1,y2,z2)-ylafun(x2,y1,z2)-ylafun(x2,y2,z1) & -ylafun(x1,y1,z1)+ylafun(x1,y1,z2)+ylafun(x1,y2,z1)+ylafun(x2,y1,z1) endif res=res/(dx*dy*dz) end function igfeyfun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function igfezfun(u,v,w,gam,dx,dy,dz) result(res) implicit none real(dp) :: res real(dp) :: u,v,w,gam,dx,dy,dz real(dp) :: x1,x2,y1,y2,z1,z2 real(dp) :: x,y,z real(dp), parameter :: ep=1.d-10 !should include em also, but probably doesn't matter !-! real(dp), external :: zlafun x1=u-0.5d0*dx x2=u+0.5d0*dx y1=v-0.5d0*dy y2=v+0.5d0*dy z1=(w-0.5d0*dz)*gam z2=(w+0.5d0*dz)*gam if(x1<0.d0.and.x2>0.d0.and.y1<0.d0.and.y2>0.d0.and.z1<0.d0.and.z2>0.d0)then res=zlafun(ep,ep,ep) & -zlafun(x1,ep,ep)-zlafun(ep,y1,ep)-zlafun(ep,ep,z1)-zlafun(x1,y1,z1) & +zlafun(x1,y1,ep)+zlafun(x1,ep,z1)+zlafun(ep,y1,z1)+zlafun(x2,ep,ep) & -zlafun(ep,ep,ep)-zlafun(x2,y1,ep)-zlafun(x2,ep,z1)-zlafun(ep,y1,z1) & +zlafun(ep,y1,ep)+zlafun(ep,ep,z1)+zlafun(x2,y1,z1)+zlafun(ep,y2,ep) & -zlafun(x1,y2,ep)-zlafun(ep,ep,ep)-zlafun(ep,y2,z1)-zlafun(x1,ep,z1) & +zlafun(x1,ep,ep)+zlafun(x1,y2,z1)+zlafun(ep,ep,z1)+zlafun(ep,ep,z2) & -zlafun(x1,ep,z2)-zlafun(ep,y1,z2)-zlafun(ep,ep,ep)-zlafun(x1,y1,ep) & +zlafun(x1,y1,z2)+zlafun(x1,ep,ep)+zlafun(ep,y1,ep)+zlafun(x2,y2,ep) & -zlafun(ep,y2,ep)-zlafun(x2,ep,ep)-zlafun(x2,y2,z1)-zlafun(ep,ep,z1) & +zlafun(ep,ep,ep)+zlafun(ep,y2,z1)+zlafun(x2,ep,z1)+zlafun(x2,ep,z2) & -zlafun(ep,ep,z2)-zlafun(x2,y1,z2)-zlafun(x2,ep,ep)-zlafun(ep,y1,ep) & +zlafun(ep,y1,z2)+zlafun(ep,ep,ep)+zlafun(x2,y1,ep)+zlafun(ep,y2,z2) & -zlafun(x1,y2,z2)-zlafun(ep,ep,z2)-zlafun(ep,y2,ep)-zlafun(x1,ep,ep) & +zlafun(x1,ep,z2)+zlafun(x1,y2,ep)+zlafun(ep,ep,ep)+zlafun(x2,y2,z2) & -zlafun(ep,y2,z2)-zlafun(x2,ep,z2)-zlafun(x2,y2,ep)-zlafun(ep,ep,ep) & +zlafun(ep,ep,z2)+zlafun(ep,y2,ep)+zlafun(x2,ep,ep) else res=zlafun(x2,y2,z2)-zlafun(x1,y2,z2)-zlafun(x2,y1,z2)-zlafun(x2,y2,z1) & -zlafun(x1,y1,z1)+zlafun(x1,y1,z2)+zlafun(x1,y2,z1)+zlafun(x2,y1,z1) endif res=res/(dx*dy*dz*gam) !note the factor of gam in the denominator here, as needed for Ez end function igfezfun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function xlafun(x,y,z) result(res) implicit none real(dp) :: res real(dp) :: x,y,z,r r=sqrt(x**2+y**2+z**2) ! if (z+r == 0 ) print *, 'ERROR:x, r, x+r', x, y, z,r, x+r ! if (y+r == 0 ) print *, 'ERROR:y, r, y+r', x, y, z,r, y+r ! if (z+r == 0 ) print *, 'ERROR:z, r, z+r', x, y, z,r, z+r ! res=x*atan(y*z/(r*x))-z*atanh(r/y)-y*atanh(r/z) !res=z-x*atan(z/x)+x*atan(y*z/(x*r))-z*log(y+r)-y*log(z+r) res=z-x*atan(z/x)+x*atan(y*z/(x*r)) if (y+r /= 0) res = res -z*log(y+r) if (z+r /= 0) res = res -y*log(z+r) ! write(2,'(5(1pe12.5,1x))')x,y,z,r,res end function xlafun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function ylafun(x,y,z) result(res) implicit none real(dp) :: res real(dp) :: x,y,z,r r=sqrt(x**2+y**2+z**2) !res=x-y*atan(x/y)+y*atan(z*x/(y*r))-x*log(z+r)-z*log(x+r) res=x-y*atan(x/y)+y*atan(z*x/(y*r)) if (z+r /= 0) res = res -x*log(z+r) if (x+r /= 0) res = res -z*log(x+r) end function ylafun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function zlafun(x,y,z) result(res) implicit none real(dp) :: res real(dp) :: x,y,z,r r=sqrt(x**2+y**2+z**2) !res=y-z*atan(y/z)+z*atan(x*y/(z*r))-y*log(x+r)-x*log(y+r) res=y-z*atan(y/z)+z*atan(x*y/(z*r)) if (x+r /= 0) res = res -y*log(x+r) if (y+r /= 0) res = res -x*log(y+r) ! TEST !res = z*atan(x*y,(r*z)) - y*atanh(r/x) - x*atanh(r/y) ! Bad !res = z * atan(x*y/z*r) + y*log(1-x/r)/2 - y*log(1+x/r)/2 ! Works !res = -z**2 * (x/(x**2 + z**2) - y/(y**2 + z**2)) + z*atan(x*y/(z*r)) -y*log(x+r) - x*log(y+r) ! Works end function zlafun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine getrhotilde(rho,crho,ilo,jlo,klo) implicit none integer :: ilo,jlo,klo real(dp), dimension(ilo:,jlo:,klo:) :: rho complex(dp), dimension(ilo:,jlo:,klo:) :: crho integer :: ihi,jhi,khi,ihi2,jhi2,khi2,iperiod,jperiod,kperiod ihi=ubound(rho,1); jhi=ubound(rho,2); khi=ubound(rho,3) ihi2=ubound(crho,1); jhi2=ubound(crho,2); khi2=ubound(crho,3) iperiod=size(crho,1); jperiod=size(crho,2); kperiod=size(crho,3) crho(ilo:ihi2,jlo:jhi2,klo:khi2)=0.d0 !zero everything before storing rho in one octant crho(ilo:ihi,jlo:jhi,klo:khi)=rho(ilo:ihi,jlo:jhi,klo:khi) !fft the charge density: call ccfft3d(crho,crho,(/1,1,1/),iperiod,jperiod,kperiod,0) !can I do the fft in place? problematic in parallel code return end subroutine !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine osc_getgrnfree(cgrn,delta,gam,ilo_grn,jlo_grn,klo_grn,npad,icomp,igfflag) implicit none real(dp), dimension(3) :: delta integer, dimension(3) :: npad real(dp) :: dx,dy,dz,gam integer :: ilo_grn,jlo_grn,klo_grn,igfflag,icomp integer :: ilo_grn_gbl,jlo_grn_gbl,klo_grn_gbl !in a parallel code these would be in the arg list integer :: ipad,jpad,kpad,ishift,jshift,kshift,iperiod,jperiod,kperiod complex(dp), dimension(ilo_grn:,jlo_grn:,klo_grn:) :: cgrn integer :: i,j,k,ip,jp,kp real(dp) :: u,v,w real(dp) :: gval dx=delta(1); dy=delta(2); dz=delta(3) ilo_grn_gbl=ilo_grn; jlo_grn_gbl=jlo_grn; klo_grn_gbl=klo_grn iperiod=size(cgrn,1); jperiod=size(cgrn,2); kperiod=size(cgrn,3) ipad=npad(1); jpad=npad(2); kpad=npad(3) !this puts the Green function where it's needed so the convolution ends up in the correct location in the array ishift=iperiod/2-(ipad+1)/2; jshift=jperiod/2-(jpad+1)/2; kshift=kperiod/2-(kpad+1)/2 !$ print *, 'OpenMP Green function calc' !$OMP PARALLEL DO & !$OMP DEFAULT(FIRSTPRIVATE), & !$OMP SHARED(cgrn) do k=klo_grn,ubound(cgrn,3) kp=klo_grn_gbl+mod(k-klo_grn_gbl+kshift ,kperiod) !in the serial code klo_grn_gbl = klo_grn, similarly for i,j w=kp*dz do j=jlo_grn,ubound(cgrn,2) jp=jlo_grn_gbl+mod(j-jlo_grn_gbl+jshift ,jperiod) v=jp*dy do i=ilo_grn,ubound(cgrn,1) ip=ilo_grn_gbl+mod(i-ilo_grn_gbl+ishift ,iperiod) u=ip*dx if(igfflag.eq.0.and.icomp.eq.0)gval=coulombfun(u,v,w,gam) if(igfflag.eq.1.and.icomp.eq.0)gval=igfcoulombfun(u,v,w,gam,dx,dy,dz) if(igfflag.eq.1.and.icomp.eq.1)gval=igfexfun(u,v,w,gam,dx,dy,dz) if(igfflag.eq.1.and.icomp.eq.2)gval=igfeyfun(u,v,w,gam,dx,dy,dz) if(igfflag.eq.1.and.icomp.eq.3)gval=igfezfun(u,v,w,gam,dx,dy,dz) cgrn(i,j,k)= cmplx(gval,0.d0, dp) enddo enddo enddo !$OMP END PARALLEL DO call ccfft3d(cgrn,cgrn,(/1,1,1/),iperiod,jperiod,kperiod,0) end subroutine osc_getgrnfree !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine conv3d(crho,qgrn1,con,rilo,rjlo,rklo,g1ilo,g1jlo,g1klo,cilo,cjlo,cklo,iperiod,jperiod,kperiod) implicit none integer :: rilo,rjlo,rklo,g1ilo,g1jlo,g1klo,cilo,cjlo,cklo,iperiod,jperiod,kperiod !input arrays: complex(dp), dimension(rilo:,rjlo:,rklo:) :: crho complex(dp), dimension(g1ilo:,g1jlo:,g1klo:) :: qgrn1 real(dp), dimension(cilo:,cjlo:,cklo:) :: con ! local: complex(dp), allocatable, dimension(:,:,:) :: ccon real(dp) :: fpei,qtot,factr integer :: cihi,cjhi,ckhi print *, '---------conv3d' fpei=299792458.d0**2*1.d-7 ! this is 1/(4 pi eps0) qtot=1.d0 !fix later: 1.d-9 ! 1 nC allocate(ccon(cilo:cilo+iperiod-1,cjlo:cjlo+jperiod-1,cklo:cklo+kperiod-1)) ccon(:,:,:)=crho(:,:,:)*qgrn1(:,:,:) call ccfft3d(ccon,ccon,(/-1,-1,-1/),iperiod,jperiod,kperiod,0) !normalize: ! factr=hx*hy*hz/( (1.d0*iperiod)*(1.d0*jperiod)*(1.d0*kperiod) )*fpei*qtot factr= 1.d0/( (1.d0*iperiod)*(1.d0*jperiod)*(1.d0*kperiod) )*fpei*qtot !store final result in original size (not double size) real array: cihi=ubound(con,1) cjhi=ubound(con,2) ckhi=ubound(con,3) con(cilo:cihi,cjlo:cjhi,cklo:ckhi)=factr*real(ccon(cilo:cihi,cjlo:cjhi,cklo:ckhi),dp) end subroutine conv3d !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine osc_rectpipe_solver(rho,a,b,gam0,delta,umin,phi,efield,bfield,nlo,nhi,nlo_gbl,nhi_gbl,idirectfieldcalc,igfflag) !This assumes that cgrn1,...,cgrn2 have been allocated and calculated prior to entering this routine. !This assumes that the rho and phi arrays have the same dimensions, called ilo,ihi,jlo,jhi,klo,khi below. !The "gbl" variables are from the parallel code, so they are ignored in this serial code. implicit none integer, intent(in) :: igfflag,idirectfieldcalc integer :: ilo,ihi,jlo,jhi,klo,khi,ilo_rho_gbl,ihi_rho_gbl,jlo_rho_gbl,jhi_rho_gbl,klo_rho_gbl,khi_rho_gbl real(dp), intent(in) :: a,b,gam0 real(dp), intent(in), dimension(3) :: delta,umin integer, intent(in), dimension(3) :: nlo,nhi,nlo_gbl,nhi_gbl real(dp), intent(in), dimension(:,:,:) :: rho real(dp), intent(out), dimension(:,:,:) :: phi real(dp), intent(out), dimension(:,:,:,:) :: efield,bfield complex(dp), allocatable, dimension(:,:,:) :: crho integer :: iperiod,jperiod,kperiod real(dp) :: dx,dy,dz,xmin,ymin,zmin integer :: ilo2,ihi2,jlo2,jhi2,klo2,khi2 integer :: icomp real(dp), parameter :: clight=299792458.d0 integer :: i,j,k,im1,ip1,jm1,jp1,km1,kp1 real(dp) :: gb0,xfac,yfac,zfac dx=delta(1); dy=delta(2); dz=delta(3) xmin=umin(1); ymin=umin(2); zmin=umin(3) ilo=nlo(1);ihi=nhi(1);jlo=nlo(2);jhi=nhi(2);klo=nlo(3);khi=nhi(3) ilo2=ilo; jlo2=jlo; klo2=klo iperiod=size(cgrn1,1); jperiod=size(cgrn1,2); kperiod=size(cgrn1,3) ihi2=ilo2+iperiod-1; jhi2=jlo2+jperiod-1; khi2=klo2+kperiod-1 write(6,*)'iperiod,jperiod,kperiod=',iperiod,jperiod,kperiod if(idirectfieldcalc.eq.0)then write(6,*)'Solving for Phi' icomp=0 !compute/store the FFT of the charge density: if(.not.allocated(crho))allocate(crho(ilo2:ihi2,jlo2:jhi2,klo2:khi2)) !double-size array call getrhotilde(rho,crho,ilo,jlo,klo) write(6,*)'done computing FFT of the charge density' !compute/store the FFT of the green function: this has been moved to the main program if(.not.allocated(cgrn1))then write(6,*)'something wrong. rectpipe green functions should have been computed during initialization' stop endif call fftconvcorr3d(crho,cgrn1,cgrn2,cgrn3,cgrn4,phi,ilo,jlo,klo, & &g1ilo,g1jlo,g1klo,g2ilo,g2jlo,g2klo,g3ilo,g3jlo,g3klo,g4ilo,g4jlo,g4klo, & &ilo,jlo,klo,iperiod,jperiod,kperiod,dx,dy,dz,xmin,ymin,zmin) write(6,*)'finished convolutions and correlations' !original version did not handle the edge values correctly when differencing; here is a first order approx: do k=lbound(phi,3),ubound(phi,3) km1=k-1; if(km1.lt.lbound(phi,3))km1=k kp1=k+1; if(kp1.gt.ubound(phi,3))kp1=k zfac=merge(2.d0,1.d0,(km1.eq.k).or.(kp1.eq.k)) do j=lbound(phi,2),ubound(phi,2) jm1=j-1; if(jm1.lt.lbound(phi,2))jm1=j jp1=j+1; if(jp1.gt.ubound(phi,2))jp1=j yfac=merge(2.d0,1.d0,(jm1.eq.j).or.(jp1.eq.j)) do i=lbound(phi,1),ubound(phi,1) im1=i-1; if(im1.lt.lbound(phi,1))im1=i ip1=i+1; if(ip1.gt.ubound(phi,1))ip1=i xfac=merge(2.d0,1.d0,(im1.eq.i).or.(ip1.eq.i)) efield(i,j,k,1)=-(phi(ip1,j,k)-phi(im1,j,k))/(2.d0*dx)*gam0*xfac efield(i,j,k,2)=-(phi(i,jp1,k)-phi(i,jm1,k))/(2.d0*dy)*gam0*yfac efield(i,j,k,3)=-(phi(i,j,kp1)-phi(i,j,km1))/(2.d0*dz)/gam0*zfac enddo enddo enddo endif if(idirectfieldcalc.eq.1)then write(6,*)'**********ERROR********* rectpipe implementation presently requires directfieldcalc=0' endif ! set the magnetic field: gb0=sqrt((gam0+1.d0)*(gam0-1.d0)) do k=lbound(phi,3),ubound(phi,3) do j=lbound(phi,2),ubound(phi,2) do i=lbound(phi,1),ubound(phi,1) bfield(i,j,k,1)=-efield(i,j,k,2)/clight/gb0/gam0 bfield(i,j,k,2)= efield(i,j,k,1)/clight/gb0/gam0 bfield(i,j,k,3)=0.d0 enddo enddo enddo end subroutine osc_rectpipe_solver !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine osc_getgrnpipe(gam,a,b,delta,umin,npad) !makes use of arrays cgrn1,cgrn2,cgrn3,cgrn4 declared at the beginning of this module implicit none real(dp), dimension(3) :: delta,umin integer, dimension(3) :: npad integer :: ipad,jpad,kpad,ishift,jshift,kshift,iperiod,jperiod,kperiod real(dp) :: gam,a,b,hx,hy,hz,xmin,ymin,zmin real(dp) :: u,v,w integer :: i,j,k,ip,jp,kp ! real(dp), external ::rfun hx=delta(1); hy=delta(2); hz=delta(3) xmin=umin(1); ymin=umin(2); zmin=umin(3) ! iperiod=size(cgrn1,1); jperiod=size(cgrn1,2); kperiod=size(cgrn1,3) ipad=npad(1); jpad=npad(2); kpad=npad(3) !this puts the Green function where it's needed so the convolution ends up in the correct location in the array ishift=iperiod/2-(ipad+1)/2; jshift=jperiod/2-(jpad+1)/2; kshift=kperiod/2-(kpad+1)/2 ! !1: do k=g1klo,ubound(cgrn1,3) do j=g1jlo,ubound(cgrn1,2) do i=g1ilo,ubound(cgrn1,1) ip=g1ilo+mod(i-g1ilo+ishift,iperiod) jp=g1jlo+mod(j-g1jlo+jshift,jperiod) kp=g1klo+mod(k-g1klo+kshift,kperiod) w=kp*hz v=jp*hy u=ip*hx cgrn1(i,j,k)=cmplx(0.d0,0.d0, dp) cgrn1(i,j,k)=rfun(u,v,w,gam,a,b,hz,1,1) enddo enddo enddo call ccfft3d(cgrn1,cgrn1,(/1,1,1/),iperiod,jperiod,kperiod,0) !2: do k=g2klo,ubound(cgrn2,3) do j=g2jlo,ubound(cgrn2,2) do i=g2ilo,ubound(cgrn2,1) ip=g2ilo+mod(i-g2ilo+ishift,iperiod) kp=g2klo+mod(k-g2klo+kshift,kperiod) w=kp*hz v=2.d0*ymin+(j-g2jlo)*hy ! v=2.d0*ymin+j*hy u=ip*hx cgrn2(i,j,k)=cmplx(0.d0,0.d0, dp) cgrn2(i,j,k)=rfun(u,v,w,gam,a,b,hz,1,-1) enddo enddo enddo call ccfft3d(cgrn2,cgrn2,(/1,-1,1/),iperiod,jperiod,kperiod,0) !3: do k=g3klo,ubound(cgrn3,3) do j=g3jlo,ubound(cgrn3,2) do i=g3ilo,ubound(cgrn3,1) jp=g3jlo+mod(j-g3jlo+jshift,jperiod) kp=g3klo+mod(k-g3klo+kshift,kperiod) w=kp*hz v=jp*hy u=2.d0*xmin+(i-g3ilo)*hx ! u=2.d0*xmin+i*hx cgrn3(i,j,k)=cmplx(0.d0,0.d0, dp) cgrn3(i,j,k)=rfun(u,v,w,gam,a,b,hz,-1,1) enddo enddo enddo call ccfft3d(cgrn3,cgrn3,(/-1,1,1/),iperiod,jperiod,kperiod,0) !4: do k=g4klo,ubound(cgrn4,3) do j=g4jlo,ubound(cgrn4,2) do i=g4ilo,ubound(cgrn4,1) kp=g4klo+mod(k-g4klo+kshift,kperiod) w=kp*hz v=2.d0*ymin+(j-g4jlo)*hy ! v=2.d0*ymin+j*hy u=2.d0*xmin+(i-g4ilo)*hx ! u=2.d0*xmin+i*hx cgrn4(i,j,k)=cmplx(0.d0,0.d0, dp) cgrn4(i,j,k)=rfun(u,v,w,gam,a,b,hz,-1,-1) enddo enddo enddo call ccfft3d(cgrn4,cgrn4,(/-1,-1,1/),iperiod,jperiod,kperiod,0) return end subroutine osc_getgrnpipe !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ function rfun(u,v,w,gam,a,b,hz,i,j) result(res) !this is the IGF version (non-IGF line is commented out) !this version contains (i**m)*(j**n) implicit none real(dp) :: res real(dp) :: u,v,w,gam,a,b,hz integer :: i,j real(dp), parameter :: pi=acos(-1.0d0) real(dp) :: ainv,binv,piainv,pibinv real(dp) :: zfun,kapmn,term integer :: m,n ainv=1.0d0/a binv=1.0d0/b piainv=pi*ainv pibinv=pi*binv ! rfun=0.d0 res=0.d0 do m=1,5 do n=1,5 kapmn=sqrt((m*pi*ainv)**2+(n*pi*binv)**2) !!!!! zfun=exp(-kapmn*abs(gam*w)) !!!!! zfun=exp(-kapmn*abs(w)) zfun=(exp(-kapmn*abs(gam*w-hz))-2.d0*exp(-kapmn*abs(gam*w))+exp(-kapmn*abs(gam*w+hz)))/(hz**2*kapmn**2) ! zfun=(exp(-kapmn*abs(w-hz))-2.d0*exp(-kapmn*abs(w))+exp(-kapmn*abs(w+hz)))/(hz**2*kapmn**2) if(w.eq.0)zfun=zfun+2.d0/(hz*kapmn) term=(i**m)*(j**n)*cos(m*u*piainv)*cos(n*v*pibinv)*zfun/kapmn ! rfun=rfun+term res=res+term enddo enddo res=res*2.d0*pi*ainv*binv !here is the correct normalization end function rfun !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine fftconvcorr3d(crho,qgrn1,qgrn2,qgrn3,qgrn4,con,rilo,rjlo,rklo, & g1ilo,g1jlo,g1klo,g2ilo,g2jlo,g2klo,g3ilo,g3jlo,g3klo,g4ilo,g4jlo,g4klo, & cilo,cjlo,cklo,iperiod,jperiod,kperiod,hx,hy,hz,xmin,ymin,zmin) implicit none integer :: rilo,rjlo,rklo,g1ilo,g1jlo,g1klo,g2ilo,g2jlo,g2klo,g3ilo,g3jlo,g3klo,g4ilo,g4jlo,g4klo real(dp) :: hx,hy,hz,xmin,ymin,zmin integer :: cilo,cjlo,cklo,iperiod,jperiod,kperiod integer :: cihi,cjhi,ckhi real(dp) :: fpei,qtot,factr !input arrays: complex(dp), dimension(rilo:,rjlo:,rklo:) :: crho complex(dp), dimension(g1ilo:,g1jlo:,g1klo:) :: qgrn1 complex(dp), dimension(g2ilo:,g2jlo:,g2klo:) :: qgrn2 complex(dp), dimension(g3ilo:,g3jlo:,g3klo:) :: qgrn3 complex(dp), dimension(g4ilo:,g4jlo:,g4klo:) :: qgrn4 real(dp), dimension(cilo:,cjlo:,cklo:) :: con !local arrays: complex(dp), allocatable, dimension(:,:,:) :: ccon,ctmp fpei=299792458.**2*1.e-7 ! this is 1/(4 pi eps0) qtot=1.d0 !fix later: 1.d-9 ! 1 nC ! allocate(ccon(cilo:cilo+iperiod-1,cjlo:cjlo+jperiod-1,cklo:cklo+kperiod-1)) allocate(ctmp(cilo:cilo+iperiod-1,cjlo:cjlo+jperiod-1,cklo:cklo+kperiod-1)) ! !1st term: ccon(:,:,:)=crho(:,:,:)*qgrn1(:,:,:) call ccfft3d(ccon,ccon,(/-1,-1,-1/),iperiod,jperiod,kperiod,0) !2nd term: ctmp(:,:,:)=crho(:,:,:)*qgrn2(:,:,:) call ccfft3d(ctmp,ctmp,(/-1,1,-1/),iperiod,jperiod,kperiod,0) ccon(:,:,:)=ccon(:,:,:)-ctmp(:,:,:) !3rd term: ctmp(:,:,:)=crho(:,:,:)*qgrn3(:,:,:) call ccfft3d(ctmp,ctmp,(/1,-1,-1/),iperiod,jperiod,kperiod,0) ccon(:,:,:)=ccon(:,:,:)-ctmp(:,:,:) !4th term: ctmp(:,:,:)=crho(:,:,:)*qgrn4(:,:,:) call ccfft3d(ctmp,ctmp,(/1,1,-1/),iperiod,jperiod,kperiod,0) ccon(:,:,:)=ccon(:,:,:)+ctmp(:,:,:) !normalize: ! factr=hx*hy*hz/( (1.d0*iperiod)*(1.d0*jperiod)*(1.d0*kperiod) )*fpei*qtot factr= 1.d0/( (1.d0*iperiod)*(1.d0*jperiod)*(1.d0*kperiod) )*fpei*qtot !store final result in original size (not double size) real array: cihi=ubound(con,1) cjhi=ubound(con,2) ckhi=ubound(con,3) con(cilo:cihi,cjlo:cjhi,cklo:ckhi)=factr*real(ccon(cilo:cihi,cjlo:cjhi,cklo:ckhi), dp) return end subroutine fftconvcorr3d !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine osc_read_rectpipe_grn implicit none real(dp) :: apipeval,bpipeval,dxval,dyval,dzval,xminval,yminval,zminval,xmaxval,ymaxval,zmaxval,gamval integer :: iloval,jloval,kloval,ihival,jhival,khival,iperval,jperval,kperval open(unit=2,file='rectpipegrn.dat',form='unformatted',status='old') read(2)apipeval,bpipeval,dxval,dyval,dzval,xminval,yminval,zminval,xmaxval,ymaxval,zmaxval, & iloval,jloval,kloval,ihival,jhival,khival,iperval,jperval,kperval,gamval read(2)cgrn1,cgrn2,cgrn3,cgrn4 close(2) write(6,*)'Done reading rectpipe transformed Green functions with these parameters:' write(6,*)'apipe,bpipe=',apipeval,bpipeval write(6,*)'dx,dy,dz=',dxval,dyval,dzval write(6,*)'xmin,ymin,zmin=',xminval,yminval,zminval write(6,*)'xmax,ymax,zmax=',xmaxval,ymaxval,zmaxval write(6,*)'ilo,jlo,klo=',iloval,jloval,kloval write(6,*)'ihi,jhi,khi=',ihival,jhival,khival write(6,*)'convolution iperiod,jperiod,kperiod=',iperval,jperval,kperval write(6,*)'gamma=',gamval end subroutine osc_read_rectpipe_grn !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine osc_write_rectpipe_grn(apipe,bpipe,delta,umin,umax,nlo,nhi,gamma) implicit none real(dp) :: apipe,bpipe,gamma real(dp), dimension(3) :: delta,umin,umax integer, dimension(3) :: nlo,nhi real(dp) :: dx,dy,dz,xmin,ymin,zmin,xmax,ymax,zmax integer :: nxlo,nylo,nzlo,nxhi,nyhi,nzhi dx=delta(1); dy=delta(2); dz=delta(3) xmin=umin(1); ymin=umin(2); zmin=umin(3) xmax=umax(1); ymax=umax(2); zmax=umax(3) nxlo=nlo(1); nylo=nlo(2); nzlo=nlo(3) nxhi=nhi(1); nyhi=nhi(2); nzhi=nhi(3) open(unit=2,file='rectpipegrn.dat',form='unformatted',status='unknown') write(2)apipe,bpipe,dx,dy,dz,xmin,ymin,zmin,xmax,ymax,zmax,nxlo,nylo,nzlo,nxhi,nyhi,nzhi, & size(cgrn1,1),size(cgrn1,2),size(cgrn1,3),gamma write(2)cgrn1,cgrn2,cgrn3,cgrn4 close(2) write(6,*)'done writing rectpipe transformed Green functions' return end subroutine osc_write_rectpipe_grn !------------------------------------------------------------------------ !------------------------------------------------------------------------ !------------------------------------------------------------------------ !+ subroutine osc_alloc_rectpipe_arrays(nlo,nhi,npad) implicit none integer, intent(in), dimension(3) :: nlo,nhi,npad integer :: rilo,rihi,rjlo,rjhi,rklo,rkhi !rho dimensions integer :: cilo,cihi,cjlo,cjhi,cklo,ckhi !phi dimensions integer :: ipad,jpad,kpad integer :: g1ihi,g1jhi,g1khi,g2ihi,g2jhi,g2khi,g3ihi,g3jhi,g3khi,g4ihi,g4jhi,g4khi rilo=nlo(1); rihi=nhi(1); rjlo=nlo(2); rjhi=nhi(2); rklo=nlo(3); rkhi=nhi(3) cilo=nlo(1); cihi=nhi(1); cjlo=nlo(2); cjhi=nhi(2); cklo=nlo(3); ckhi=nhi(3) ipad=npad(1); jpad=npad(2); kpad=npad(3) g1ilo=cilo-rihi; g1ihi=cihi-rilo; g1jlo=cjlo-rjhi; g1jhi=cjhi-rjlo; g1klo=cklo-rkhi; g1khi=ckhi-rklo g2ilo=cilo-rihi; g2ihi=cihi-rilo; g2jlo=cjlo+rjlo; g2jhi=cjhi+rjhi; g2klo=cklo-rkhi; g2khi=ckhi-rklo g3ilo=cilo+rilo; g3ihi=cihi+rihi; g3jlo=cjlo-rjhi; g3jhi=cjhi-rjlo; g3klo=cklo-rkhi; g3khi=ckhi-rklo g4ilo=cilo+rilo; g4ihi=cihi+rihi; g4jlo=cjlo+rjlo; g4jhi=cjhi+rjhi; g4klo=cklo-rkhi; g4khi=ckhi-rklo if(.not.allocated(cgrn1))allocate(cgrn1(g1ilo:g1ihi+ipad,g1jlo:g1jhi+jpad,g1klo:g1khi+kpad)) if(.not.allocated(cgrn2))allocate(cgrn2(g2ilo:g2ihi+ipad,g2jlo:g2jhi+jpad,g2klo:g2khi+kpad)) if(.not.allocated(cgrn3))allocate(cgrn3(g3ilo:g3ihi+ipad,g3jlo:g3jhi+jpad,g3klo:g3khi+kpad)) if(.not.allocated(cgrn4))allocate(cgrn4(g4ilo:g4ihi+ipad,g4jlo:g4jhi+jpad,g4klo:g4khi+kpad)) end subroutine osc_alloc_rectpipe_arrays end module open_spacecharge_core_mod