dw3d.F



[divxy] [getese] [getese3d] [getese3dfromrhophi] [geteserzfromrhophi] [gtlchg] [gtlchg3d] [gtlchg3dfromrho] [gtlchgrzfromrho] [multpole] [multpwrk] [pltfld3d] [prntpa3d] [rhodia] [rhodia3d] [sezax] [sezax3d] [sphiax] [sphiax3d] [srhoax] [srhoax3d]

#include top.h
 @(#) File DW3D.F, version $Revision: 3.173 $, $Date: 2012/01/04 17:20:20 $
 # Copyright (c) 1990-1998, The Regents of the University of California.
 # All rights reserved.  See LEGAL.LLNL for full text and disclaimer.
   This is part of the package W3D of the code WARP
   3d electrostatic PIC code, Cartesian geometry, for beam problems
   This file contains the diagnostic and plotting routines.
   Alex Friedman,  LLNL, (510) 422-0827
   David P. Grote, LLNL, (510) 423-7194

[w3dgen] [wxygen]
      subroutine prntpa3d(lprntpara)
      use Beam_acc
      use InMesh3d
      use Fields3d
      use InGen
      use Lattice
      use Picglb3d
      use Io
      use Ch_var
      logical(ISZ):: lprntpara

   Prints out various parameters to a plot frame and an output file or tty

      character(80):: textline
      real(kind=8):: xrbend,xbbend,xbendlen,xstralen,xz0bend

  Exit now if parameter are not to be printed
      if (.not. lprntpara) return

      --- Call script version of this routine.
      call callpythonfunc("printparameters3d","printparameters3d")

   20 format(1x,a,1pe11.4,a)
   30 format(1x,a,i8,a)

   Gather up various lattice qtys for printout

      if (nbend >= 1) then
         xrbend = bendrc(1)
         if (ndipo >= 1) then
           xbbend = dipoby(1)
         else
           xbbend = 0.
         endif
         xbendlen = bendze(1) - bendzs(1)
         xstralen = bendzs(1) - bendze(0)
         xz0bend = bendzs(1)
      else
         xrbend = 0.
         xbbend = 0.
         xbendlen = 0.
         xstralen = 0.
         xz0bend = 0.
      endif

   Write to plot frame
      write (textline,30) "Number of grid points in x = ",nx," "
      call remark(trim(textline))
      write (textline,30) "Number of grid points in y = ",ny," "
      call remark(trim(textline))
      write (textline,30) "Number of grid points in z = ",nz," "
      call remark(trim(textline))
      write (textline,20) "Grid spacing in x = ",dx," m"
      call remark(trim(textline))
      write (textline,20) "Grid spacing in y = ",dy," m"
      call remark(trim(textline))
      write (textline,20) "Grid spacing in z = ",dz," m"
      call remark(trim(textline))
      write (textline,20) "Bend radius = ",xrbend," m"
      call remark(trim(textline))
      write (textline,20) "Bending field = ",xbbend," T"
      call remark(trim(textline))
      write (textline,20) "Bend length = ",xbendlen," m"
      call remark(trim(textline))
      write (textline,20) "Straight section length = ",xstralen," m"
      call remark(trim(textline))
      write (textline,20) "Z at start of first bend = ",xz0bend," m"
      call remark(trim(textline))

   Write to text output file
      if (warpout > -1) then
        call edit (warpout, "nx")
        call edit (warpout, "ny")
        call edit (warpout, "nz")
        call edit (warpout, "dx")
        call edit (warpout, "dy")
        call edit (warpout, "dz")
      endif

      return
      end

[step3d]
      subroutine rhodia
      use Subtimersw3d
      use InGen,Only: fstype
      use InGen3d
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (fstype == 12) then
        call callpythonfunc("rhodiaregistered","fieldsolver")
      else if(solvergeom==XYZgeom) then
        call rhodia3d
      end if
!$OMP MASTER
      if (lw3dtimesubs) timerhodia = timerhodia + wtime() - substarttime
!$OMP END MASTER
      return
      end

[step3d]
      subroutine gtlchg
      use Subtimersw3d
      use InGen,Only: fstype
      use InGen3d
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (fstype == 12) then
        call callpythonfunc("gtlchgregistered","fieldsolver")
      else if(solvergeom==XYZgeom) then
        call gtlchg3d
      else if(solvergeom==RZgeom .or. solvergeom==XZgeom) then
        call gtlchgrz
      end if
!$OMP MASTER
      if (lw3dtimesubs) timegtlchg = timegtlchg + wtime() - substarttime
!$OMP END MASTER
      return
      end

[step3d]
      subroutine srhoax
      use Subtimersw3d
      use InGen,Only: fstype
      use InGen3d
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (fstype == 12) then
        call callpythonfunc("srhoaxregistered","fieldsolver")
      else if(solvergeom==XYZgeom) then
        call srhoax3d
      end if
!$OMP MASTER
      if (lw3dtimesubs) timesrhoax = timesrhoax + wtime() - substarttime
!$OMP END MASTER
      return
      end

[step3d]
      subroutine getese
      use Subtimersw3d
      use InGen,Only: fstype
      use InGen3d
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (fstype == 12) then
        call callpythonfunc("geteseregistered","fieldsolver")
      else if(solvergeom==XYZgeom) then
        call getese3d
      end if
!$OMP MASTER
      if (lw3dtimesubs) timegetese = timegetese + wtime() - substarttime
!$OMP END MASTER
      return
      end

[step3d]
      subroutine sphiax
      use Subtimersw3d
      use InGen,Only: fstype
      use InGen3d
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (fstype == 12) then
        call callpythonfunc("sphiaxregistered","fieldsolver")
      else if(solvergeom==XYZgeom) then
        call sphiax3d
      end if
!$OMP MASTER
      if (lw3dtimesubs) timesphiax = timesphiax + wtime() - substarttime
!$OMP END MASTER
      return
      end

[step3d]
      subroutine sezax
      use Subtimersw3d
      use InGen,Only: fstype
      use InGen3d
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (fstype == 12) then
        call callpythonfunc("sezaxregistered","fieldsolver")
      else if(solvergeom==XYZgeom) then
        call sezax3d
      end if
!$OMP MASTER
      if (lw3dtimesubs) timesezax = timesezax + wtime() - substarttime
!$OMP END MASTER
      return
      end

[gtlchg] [stepxy]
      subroutine gtlchg3d
      use Subtimersw3d
      use Picglb
      use InGen3d
      use Picglb3d
      use InMesh3d
      use Fields3d
      use Z_arrays
      use InDiag3d
      use Parallel

   Calculates the line charge density from rho.

      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (.not. lgtlchg3d .or. .not. associated(rho)) return

      call gtlchg3dfromrho(nxlocal,nylocal,nzlocal,
     &                     nxguardrho,nyguardrho,nzguardrho,
     &                     rho,dx,dy,dz,
     &                     zgrid,zmminlocal,l2symtry,l4symtry,
     &                     izproc==nzprocs-1)

!$OMP MASTER
      if (lw3dtimesubs) timegtlchg3d = timegtlchg3d + wtime() - substarttime
!$OMP END MASTER
      return
      end

[gtlchg3d]
      subroutine gtlchg3dfromrho(nxlocal,nylocal,nzlocal,
     &                           nxguardrho,nyguardrho,nzguardrho,
     &                           rho,dx,dy,dz,
     &                           zgrid,zmminlocal,l2symtry,l4symtry,islastproc)
      use Subtimersw3d
      use Picglb,Only: zbeam
      use Z_arrays,Only: linechg,zzmin,dzz,nzzarr
      use InDiag3d,Only: lgtlchg3d

      integer(ISZ):: nxlocal,nylocal,nzlocal
      integer(ISZ):: nxguardrho,nyguardrho,nzguardrho
      real(kind=8):: rho(-nxguardrho:nxlocal+nxguardrho,
     &                   -nyguardrho:nylocal+nyguardrho,
     &                   -nzguardrho:nzlocal+nzguardrho)
      real(kind=8):: dx,dy,dz,zgrid,zmminlocal
      logical(ISZ):: l2symtry,l4symtry,islastproc

   Calculates the line charge density from the rho argument.

      real(kind=8),pointer:: scrtch(:)
      real(kind=8):: zz,zg,wzg,dzi
      integer(ISZ):: iz,izg
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (.not. lgtlchg3d) return

      allocate(scrtch(0:nzlocal))

      if (nzlocal == 0) then
        dzi = 0.
      else
        dzi = 1./dz
      endif

      conversion factor to go from grid frame to beam frame
      zz = zgrid + zmminlocal - zzmin - zbeam

  First, sum up the charge density in each axial plane.

!$OMP PARALLEL DO
      do iz=0,nzlocal

        if (l2symtry) then
          --- 2-fold: account for only one side of rho being stored
          scrtch(iz) = sum(rho(0:nxlocal,0,iz)) +
     &                 2.*sum(rho(0:nxlocal,1:nylocal,iz))

        elseif (l4symtry) then
          --- 4-fold: account for only one quadrant of rho being stored
          scrtch(iz) = rho(0,0,iz) + 4.*sum(rho(1:nxlocal,1:nylocal,iz)) +
     &                               2.*sum(rho(1:nxlocal,0,iz)) +
     &                               2.*sum(rho(0,1:nylocal,iz))

        else
          --- no symmetries...
          scrtch(iz) = sum(rho(0:nxlocal,0:nylocal,iz))

        endif

      enddo
!$OMP END PARALLEL DO

  Then, linearly interpolate the data into the beam frame array, linechg.

      do iz=0,nzzarr

        zg = (iz*dzz - zz)*dzi
        izg = int(zg + nzlocal) - nzlocal
        wzg = zg - izg

        --- This is a bit tricky for the parallel version. Each processor
        --- must set linechg only in the range of 0 through nzlocal-1, except
        --- the last processor, which needs to set nzlocal as well. This is
        --- needed because of the overalap between nzlocal on one processor
        --- and 0 on the next processor. The last processor has no overlap
        --- at nzlocal.
        if (0 <= izg .and.
     &      (izg < nzlocal .or.
     &       (izg == nzlocal .and. islastproc))) then
          linechg(iz) = scrtch(izg)*dx*dy*(1. - wzg)
        else
          linechg(iz) = 0.
        endif
        if (0 < izg+1 .and. izg+1 <= nzlocal) then
          linechg(iz) = linechg(iz) + scrtch(izg+1)*dx*dy*wzg
        endif

      enddo

#ifdef MPIPARALLEL
      call parallelnonzerorealarray(linechg,nzzarr+1)
#endif

      deallocate(scrtch)

!$OMP MASTER
      if (lw3dtimesubs) timegtlchg3dfromrho = timegtlchg3dfromrho + wtime() - substarttime
!$OMP END MASTER
      return
      end

      subroutine gtlchgrzfromrho(nxlocal,nzlocal,
     &                           nxguardrho,nzguardrho,
     &                           rho,dx,dz,
     &                           zgrid,zmminlocal,islastproc)
      use Constant, Only: pi
      use Subtimersw3d
      use Picglb,Only: zbeam
      use Z_arrays,Only: linechg,zzmin,dzz,nzzarr
      use InDiag3d,Only: lgtlchg3d

      integer(ISZ):: nxlocal,nzlocal
      integer(ISZ):: nxguardrho,nzguardrho
      real(kind=8):: rho(-nxguardrho:nxlocal+nxguardrho,
     &                   -nzguardrho:nzlocal+nzguardrho)
      real(kind=8):: dx,dy,dz,zgrid,zmminlocal
      logical(ISZ):: l2symtry,l4symtry,islastproc

   Calculates the line charge density from the rho argument.

      real(kind=8),pointer:: scrtch(:)
      real(kind=8),pointer:: cellarea(:)
      real(kind=8):: zz,zg,wzg,dzi
      integer(ISZ):: ix,iz,izg
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()
      if (.not. lgtlchg3d) return

      allocate(scrtch(0:nzlocal))
      allocate(cellarea(0:nxlocal))

      if (nzlocal == 0) then
        dzi = 0.
      else
        dzi = 1./dz
      endif

      cellarea(0) = (pi*(0.5*0.5*dx*dx))/0.75
      do ix = 1, nxlocal
        cellarea(ix) = 2.*pi*ix*dx*dx
      enddo

      conversion factor to go from grid frame to beam frame
      zz = zgrid + zmminlocal - zzmin - zbeam

  First, sum up the charge density in each axial plane.

!$OMP PARALLEL DO
      do iz=0,nzlocal
        scrtch(iz) = sum(rho(0:nxlocal,iz)*cellarea)
      enddo
!$OMP END PARALLEL DO

  Then, linearly interpolate the data into the beam frame array, linechg.

      do iz=0,nzzarr

        zg = (iz*dzz - zz)*dzi
        izg = int(zg + nzlocal) - nzlocal
        wzg = zg - izg

        --- This is a bit tricky for the parallel version. Each processor
        --- must set linechg only in the range of 0 through nzlocal-1, except
        --- the last processor, which needs to set nzlocal as well. This is
        --- needed because of the overalap between nzlocal on one processor
        --- and 0 on the next processor. The last processor has no overlap
        --- at nzlocal.
        if (0 <= izg .and.
     &      (izg < nzlocal .or.
     &       (izg == nzlocal .and. islastproc))) then
          linechg(iz) = scrtch(izg)*(1. - wzg)
        else
          linechg(iz) = 0.
        endif
        if (0 < izg+1 .and. izg+1 <= nzlocal) then
          linechg(iz) = linechg(iz) + scrtch(izg+1)*wzg
        endif

      enddo

#ifdef MPIPARALLEL
      call parallelnonzerorealarray(linechg,nzzarr+1)
#endif

      deallocate(scrtch)
      deallocate(cellarea)

!$OMP MASTER
      if (lw3dtimesubs) timegtlchgrzfromrho = timegtlchgrzfromrho + wtime() - substarttime
!$OMP END MASTER
      return
      end

[padvncxy] [rhodia]
      subroutine rhodia3d
      use InGen3d
      use InDiag
      use InDiag3d
      use InMesh3d
      use Picglb
      use Picglb3d
      use Fields3d
      use Win_Moments
      use Z_Moments

      integer(ISZ):: iwin,ix,iy,iz,ii
      real(kind=8):: zwin,z0,z1,ww,dzi

   Sets rho window diagnostics (they need data from the 3d mesh)

      if (.not. lrhodia3d .or. .not. associated(rho)) return

      if (nz == 0) then
        dzi = 0.
      else
        dzi = 1./dz
      endif

      rhomid = 0.
      rhomax = 0.
      if (0 <= ix_axis .and. ix_axis <= nxlocal .and.
     &    0 <= iy_axis .and. iy_axis <= nylocal) then

!$OMP PARALLEL DO PRIVATE(zwin,iz,z1,z0)
        do iwin = 0,nzwind
          zwin = 0.5 * (zwindows(1,iwin) + zwindows(2,iwin))
          iz = (zwin + zbeam - zgrid - zmminlocal)*dzi
          if (zwindows(1,iwin).ne.zwindows(2,iwin) .and.
     &        0 <= iz .and. iz < nzlocal) then
            z1 = (zwin + zbeam - zgrid - zmminlocal)*dzi - iz
            z0 = 1. - z1
            rhomid(iwin) = z0*rho(ix_axis,iy_axis,iz) +
     &                     z1*rho(ix_axis,iy_axis,iz+1)
            do iy = 0,nylocal
              do ix = 0,nxlocal
                 rhomax(iwin) = max(rhomax(iwin),
     &                                z0*rho(ix,iy,iz) + z1*rho(ix,iy,iz+1))
              enddo
            enddo
          endif
        enddo
!$OMP END PARALLEL DO

      endif

#ifdef MPIPARALLEL
      call parallelnonzerorealarray(rhomid,nzwind+1)
      call parallelmaxrealarray(rhomax,nzwind+1)
#endif

   Sets rho grid diagnostics (they need data from the 3d mesh)
      if (ifzmmnt > 0) then

        rhomidz = 0.
        if (0 <= ix_axis .and. ix_axis <= nxlocal .and.
     &      0 <= iy_axis .and. iy_axis <= nylocal) then

          --- rho on axis
          do iz = 0,nzmmnt
             ii = (zmntmesh(iz) + zbeam - zgrid - zmminlocal)*dzi
             ww = (zmntmesh(iz) + zbeam - zgrid - zmminlocal)*dzi - ii
             if (ii >= 0 .and. ii < nzlocal)
     &         rhomidz(iz) = rho(ix_axis,iy_axis,ii)*(1. - ww) +
     &                       rho(ix_axis,iy_axis,ii+1)*ww
          enddo

        endif

        --- max rho at each z (vectorized)
        do iz = 0,nzmmnt
           rhomaxz(iz) = -LARGEPOS
        enddo
!$OMP PARALLEL DO PRIVATE(ii,ww)
        do iz = 0,nzmmnt
          ii = (zmntmesh(iz) + zbeam - zgrid - zmminlocal)*dzi
          ww = (zmntmesh(iz) + zbeam - zgrid - zmminlocal)*dzi - ii
          if (ii >= 0 .and. ii < nzlocal) then
            do iy = 0,nylocal
              do ix = 0,nxlocal
                rhomaxz(iz) = max(rhomaxz(iz),rho(ix,iy,ii)*(1. - ww) +
     &                                        rho(ix,iy,ii+1)*ww)
              enddo
            enddo
          elseif (ii == nzlocal) then
            do iy = 0,nylocal
              do ix = 0,nxlocal
                rhomaxz(iz) = max(rhomaxz(iz),rho(ix,iy,nzlocal))
              enddo
            enddo
          endif
        enddo
!$OMP END PARALLEL DO
      endif

#ifdef MPIPARALLEL
      call parallelnonzerorealarray(rhomidz,nzmmnt+1)
      call parallelmaxrealarray(rhomaxz,nzmmnt+1)
#endif

      return
      end

[srhoax] [stepxy]
      subroutine srhoax3d
   Sets 1d array for the charge density on the axis
      use InGen3d
      use InMesh3d
      use InDiag3d
      use Picglb
      use Picglb3d
      use Fields3d
      use Z_arrays
#ifdef MPIPARALLEL
      use Parallel
#endif

      integer(ISZ):: ix,iy,iz,iz0,izg
      real(kind=8):: zz,wzg,dzi

      if(.not. lsrhoax3d .or. .not. associated(rho)) return

      rhoax = 0.

      if (0 <= ix_axis .and. ix_axis <= nxlocal .and.
     &    0 <= iy_axis .and. iy_axis <= nylocal) then

        if (nz == 0) then
          dzi = 0.
        else
          dzi = 1./dz
        endif

        ix = ix_axis
        iy = iy_axis

#ifdef MPIPARALLEL
        iz0 = fsdecomp%iz(fsdecomp%izproc)
#else
        iz0 = 0
#endif

        zz = zzmin + zbeam - zgrid - zmmin

        if (nzlocal == nzzarr .and. abs(zz) <= dz*1.e-5) then
          do iz=0,nzlocal
            rhoax(iz+iz0) = rho(ix,iy,iz)
          enddo
        else
          do iz=0,nzzarr
            izg = int((iz*dzz + zz)*dzi) - iz0
            if (0 <= izg .and. izg < nzlocal) then
              wzg = (iz*dzz + zz)*dzi - iz0 - izg
              rhoax(iz) = rho(ix,iy,izg)*(1. - wzg) + rho(ix,iy,izg+1)*wzg
            else if (izg == nzlocal) then
              rhoax(iz) = rho(ix,iy,izg)
            endif
          enddo
        endif

      endif

#ifdef MPIPARALLEL
      call parallelnonzerorealarray(rhoax,nzzarr+1)
#endif

      return
      end

[getese] [stepxy]
      subroutine getese3d
      use Subtimersw3d
      use Beam_acc
      use Picglb
      use Picglb3d
      use InMesh3d
      use Fields3d
      use Moments
      use InPart
      use InGen
      use InGen3d
      use Z_arrays
      use InDiag3d

  Calculate the electrostatic energy, sum of rho*phi*dx*dy*dz/2
  Includes external fields: Uniform focusing and eears for cigars
  This does it in a partially vectorized manner.

      integer(ISZ):: ix,iy,iz,izb
      real(kind=8):: zm,zz,phiextun,wzb
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()

      if (.not. lgetese3d .or.
     &    .not. associated(rho) .or. .not. associated(phi)) then
!$OMP MASTER
        if (lw3dtimesubs) timegetese3d = timegetese3d + wtime() - substarttime
!$OMP END MASTER
        return
      endif

      ese = 0.
      zm = (zimax - zimin)*(1. - straight)*0.5

      conversion factor to go from grid frame to beam frame
      zz = zmminlocal + zgrid - zbeam

      do iy=0,nylocal
        do ix=0,nxlocal
          scrtch(ix,iy) = 0.0
        enddo
      enddo

      --- Loop order switched for the OpenMP parallel version - this
      --- avoids requiring the accumulation into scrtch to be put
      --- into a critical section.
!$OMP PARALLEL
!$OMP DO PRIVATE(phiextun,izb,wzb,ix,iy,iz)
#ifdef _OPENMP
      do iy=0,nylocal
#endif
      do iz=0,nzlocal-1
        phiextun = 0.
        if (ifeears == 1) then
          if (iz*dz+zz < -zm) phiextun=eears*.5*(iz*dz+zz+zm)**2
          if (iz*dz+zz >  zm) phiextun=eears*.5*(iz*dz+zz-zm)**2
        elseif (ifeears == 2) then
          --- This is not quite right for the parallel version since
          --- the longitudinal extent of eearsofz (the same as the extent
          --- of the particles) will not necessarily overlap the extent of
          --- the field solver.
          izb = (iz*dz+zz-zzmin)*dzzi
          wzb = (iz*dz+zz-zzmin)*dzzi - izb
          if (0 <= iz .and. iz < nzzarr)
     &      phiextun = eearsofz(izb)*(1. - wzb) + eearsofz(izb+1)*wzb
        endif
#ifndef _OPENMP
        do iy=0,nylocal
#endif
          do ix=0,nxlocal
            scrtch(ix,iy) = scrtch(ix,iy) + rho(ix,iy,iz)*(phi(ix,iy,iz)
     &                       - dedr*.5*((ix*dx+xmmin)**2 + (iy*dy+ymmin)**2)
     &                       - dexdx*.5*(ix*dx+xmmin)**2 
     &                       - deydy*.5*(iy*dy+ymmin)**2
     &                       + phiextun)
          enddo
        enddo
      enddo
!$OMP END DO
!$OMP END PARALLEL

      --- Apply appropriate multiplication factors for symmetries
      if (l4symtry) then
        --- Center point is unchanged
        --- Edges are doubled
        --- Bulk is quadrupoled
        scrtch(0,1:nylocal) = 2.*scrtch(0,1:nylocal)
        scrtch(1:nxlocal,0) = 2.*scrtch(1:nxlocal,0)
        scrtch(1:nxlocal,1:nylocal) = 4.*scrtch(1:nxlocal,1:nylocal)
      else if (l2symtry) then
        --- Edge is unchanged
        --- Bulk is doubled
        scrtch(1:nxlocal,0:nylocal) = 2.*scrtch(1:nxlocal,0:nylocal)
      endif

      do iy=1,nylocal
        do ix=0,nxlocal
          scrtch(ix,0) = scrtch(ix,0) + scrtch(ix,iy)
        enddo
      enddo

      do ix=0,nxlocal
        ese = ese + scrtch(ix,0)
      enddo

      ese = dx*dy*dz*0.5*ese

!$OMP MASTER
      if (lw3dtimesubs) timegetese3d = timegetese3d + wtime() - substarttime
!$OMP END MASTER
      return
      end

      subroutine getese3dfromrhophi(nxlocal,nylocal,nzlocal,
     &                              nxguardphi,nyguardphi,nzguardphi,
     &                              nxguardrho,nyguardrho,nzguardrho,
     &                              rho,phi,dx,dy,dz,l4symtry,l2symtry,ese)
      use Subtimersw3d
      integer(ISZ):: nxlocal,nylocal,nzlocal
      integer(ISZ):: nxguardphi,nyguardphi,nzguardphi
      integer(ISZ):: nxguardrho,nyguardrho,nzguardrho
      real(kind=8):: rho(-nxguardrho:nxlocal+nxguardrho,
     &                   -nyguardrho:nylocal+nyguardrho,
     &                   -nzguardrho:nzlocal+nzguardrho)
      real(kind=8):: phi(-nxguardphi:nxlocal+nxguardphi,
     &                   -nyguardphi:nylocal+nyguardphi,
     &                   -nzguardphi:nzlocal+nzguardphi)
      real(kind=8):: dx,dy,dz,ese
      logical(ISZ):: l4symtry,l2symtry

  Calculate the electrostatic energy, sum of rho*phi*dx*dy*dz/2
  This does it in a partially vectorized manner.

      integer(ISZ):: ix,iy,iz
      real(kind=8),pointer:: scrtch(:,:)
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()

      allocate(scrtch(0:nxlocal,0:nylocal))
      scrtch = 0.0

      --- Note that the z planes iz=0 and nzlocal may need skipping depending
      --- on the boundary conditions.

      do iz=0,nzlocal
        do iy=0,nylocal
          do ix=0,nxlocal
            scrtch(ix,iy) = scrtch(ix,iy) + rho(ix,iy,iz)*phi(ix,iy,iz)
          enddo
        enddo
      enddo

      --- Apply appropriate multiplication factors for symmetries
      if (l4symtry) then
        --- Center point is unchanged
        --- Edges are doubled
        --- Bulk is quadrupoled
        scrtch(0,1:nylocal) = 2.*scrtch(0,1:nylocal)
        scrtch(1:nxlocal,0) = 2.*scrtch(1:nxlocal,0)
        scrtch(1:nxlocal,1:nylocal) = 4.*scrtch(1:nxlocal,1:nylocal)
      else if (l2symtry) then
        --- Edge is unchanged
        --- Bulk is doubled
        scrtch(1:nxlocal,0:nylocal) = 2.*scrtch(1:nxlocal,0:nylocal)
      endif

      do iy=1,nylocal
        do ix=0,nxlocal
          scrtch(ix,0) = scrtch(ix,0) + scrtch(ix,iy)
        enddo
      enddo

      do ix=0,nxlocal
        ese = ese + scrtch(ix,0)
      enddo

      ese = dx*dy*dz*0.5*ese

      deallocate(scrtch)

!$OMP MASTER
      if (lw3dtimesubs) timegetese3dfromrhophi = timegetese3dfromrhophi + wtime() - substarttime
!$OMP END MASTER
      return
      end

      subroutine geteserzfromrhophi(nxlocal,nzlocal,
     &                              nxguardphi,nzguardphi,
     &                              nxguardrho,nzguardrho,
     &                              rho,phi,dx,dz,xmminlocal,ese)
      use Subtimersw3d
      use Constant
      integer(ISZ):: nxlocal,nzlocal
      integer(ISZ):: nxguardphi,nzguardphi
      integer(ISZ):: nxguardrho,nzguardrho
      real(kind=8):: rho(-nxguardrho:nxlocal+nxguardrho,
     &                   -nzguardrho:nzlocal+nzguardrho)
      real(kind=8):: phi(-nxguardphi:nxlocal+nxguardphi,
     &                   -nzguardphi:nzlocal+nzguardphi)
      real(kind=8):: dx,dz,xmminlocal,ese

  Calculate the electrostatic energy, sum of r*rho*phi*dx*dz/2
  This does it in a partially vectorized manner.

      integer(ISZ):: ix,iz
      real(kind=8):: r
      real(kind=8),pointer:: scrtch(:)
      real(kind=8):: substarttime,wtime
      if (lw3dtimesubs) substarttime = wtime()

      allocate(scrtch(0:nxlocal))
      scrtch = 0.0

      --- Note that the z planes iz=0 and nzlocal may need skipping depending
      --- on the boundary conditions.

      do iz=0,nzlocal
        do ix=0,nxlocal
          scrtch(ix) = scrtch(ix) + rho(ix,iz)*phi(ix,iz)
        enddo
      enddo

      do ix=0,nxlocal
        r = xmminlocal + ix*dx
        ese = ese + r*scrtch(ix)
      enddo

      ese = 2.*pi*dx*dz*0.5*ese

      deallocate(scrtch)

!$OMP MASTER
      if (lw3dtimesubs) timegeteserzfromrhophi = timegeteserzfromrhophi + wtime() - substarttime
!$OMP END MASTER
      return
      end

[sphiax] [stepxy]
      subroutine sphiax3d
   Sets 1d array for the potential on the axis
      use InGen3d
      use InMesh3d
      use InDiag3d
      use Z_arrays
      use Picglb
      use Picglb3d
      use Fields3d
#ifdef MPIPARALLEL
      use Parallel
#endif

      integer(ISZ):: ix,iy,iz,iz0,izg
      real(kind=8):: zz,wzg,dzi

      if(.not. lsphiax3d .or. .not. associated(phi)) return

      phiax = 0.

      if (0 <= ix_axis .and. ix_axis <= nxlocal .and.
     &    0 <= iy_axis .and. iy_axis <= nylocal) then

        if (nz == 0) then
          dzi = 0.
        else
          dzi = 1./dz
        endif

        ix = ix_axis
        iy = iy_axis
      
#ifdef MPIPARALLEL
        iz0 = fsdecomp%iz(fsdecomp%izproc)
#else
        iz0 = 0
#endif

        zz = zzmin + zbeam - zgrid - zmmin

        if (nz == nzzarr .and. abs(zz) <= dz*1.e-5) then
          do iz=0,nzlocal
            phiax(iz+iz0) = phi(ix,iy,iz)
          enddo
        else
          do iz=0,nzzarr
            izg = int((iz*dzz + zz)*dzi) - iz0
            if (0 <= izg .and. izg <= nzlocal) then
              wzg = (iz*dzz + zz)*dzi - iz0  - izg
              phiax(iz) = phi(ix,iy,izg)*(1. - wzg) + phi(ix,iy,izg+1)*wzg
            endif
          enddo
        endif

      endif

#ifdef MPIPARALLEL
      call parallelnonzerorealarray(phiax,nzzarr+1)
#endif

      return
      end

[sezax] [stepxy]
      subroutine sezax3d
   Sets 1d array for the space charge E field on the axis
      use InMesh3d
      use InGen3d
      use InDiag3d
      use Picglb
      use Picglb3d
      use Z_arrays
      use Fields3d
#ifdef MPIPARALLEL
      use Parallel
#endif

      integer(ISZ):: ix,iy,iz,iz0,izg
      real(kind=8):: wx,wy,zz,wzg,wx1,wy1,dzi

      if(.not. lsezax3d .or. .not. associated(phi)) return

      ezax = 0.

      if (0 <= ix_axis .and. ix_axis <= nxlocal .and.
     &    0 <= iy_axis .and. iy_axis <= nylocal) then

        if (nz == 0) then
          dzi = 0.
        else
          dzi = 1./dz
        endif

        ix = ix_axis
        wx = ix_axis - ix
        wx1 = 1. - wx
        iy = iy_axis
        wy = iy_axis - iy
        wy1 = 1. - wy

#ifdef MPIPARALLEL
        iz0 = fsdecomp%iz(fsdecomp%izproc)
#else
        iz0 = 0
#endif

        zz = zzmin + zbeam - zgrid - zmmin

        if ((xmmax == -xmmin .and.  ymmax == -ymmin) .or.
     &      (wx == 0. .and. wy == 0.)) then
          do iz=0,nzzarr
            izg = int((iz*dzz + zz)*dzi) - iz0
            if (0 <= izg .and. izg <= nzlocal-1) then
              wzg = (iz*dzz + zz)*dzi - iz0 - izg
              ezax(iz) = ((phi(ix,iy,izg-1)-phi(ix,iy,izg+1))*(1.-wzg) +
     &                    (phi(ix,iy,izg  )-phi(ix,iy,izg+2))*    wzg )*0.5*dzi
            else if (izg == nzlocal) then
              ezax(iz) = (phi(ix,iy,izg-1) - phi(ix,iy,izg+1))*0.5*dzi
            endif
          enddo
        else
          do iz=0,nzzarr
            izg = int((iz*dzz + zz)*dzi) - iz0
            if (0 <= izg .and. izg <= nzlocal) then
              wzg = (iz*dzz + zz)*dzi - iz0 - izg
          ezax(iz)=
     &          ((phi(ix  ,iy  ,izg-1)-phi(ix  ,iy  ,izg+1))*wx1*wy1*(1.-wzg)+
     &           (phi(ix+1,iy  ,izg-1)-phi(ix+1,iy  ,izg+1))*wx *wy1*(1.-wzg)+
     &           (phi(ix  ,iy+1,izg-1)-phi(ix  ,iy+1,izg+1))*wx1*wy *(1.-wzg)+
     &           (phi(ix+1,iy+1,izg-1)-phi(ix+1,iy+1,izg+1))*wx *wy *(1.-wzg)+
     &           (phi(ix  ,iy  ,izg  )-phi(ix  ,iy  ,izg+2))*wx1*wy1*    wzg +
     &           (phi(ix+1,iy  ,izg  )-phi(ix+1,iy  ,izg+2))*wx *wy1*    wzg +
     &           (phi(ix  ,iy+1,izg  )-phi(ix  ,iy+1,izg+2))*wx1*wy *    wzg +
     &           (phi(ix+1,iy+1,izg  )-phi(ix+1,iy+1,izg+2))*wx *wy *    wzg )
     &           *0.5*dzi
            else if (izg == nzlocal) then
              ezax(iz)=((phi(ix  ,iy  ,izg-1)-phi(ix  ,iy  ,izg+1))*wx1*wy1 +
     &                  (phi(ix+1,iy  ,izg-1)-phi(ix+1,iy  ,izg+1))*wx *wy1 +
     &                  (phi(ix  ,iy+1,izg-1)-phi(ix  ,iy+1,izg+1))*wx1*wy  +
     &                  (phi(ix+1,iy+1,izg-1)-phi(ix+1,iy+1,izg+1))*wx *wy)
     &                  *0.5*dzi
            endif
          enddo
        endif

      endif

#ifdef MPIPARALLEL
      call parallelnonzerorealarray(ezax,nzzarr+1)
#endif

      return
      end

[step3d]
      subroutine pltfld3d(fld,freqflag)
      use InDiag3d
      use InPltCtl3d
      use Timers
      character(3):: fld
      integer(ISZ):: freqflag

   Master control for doing field plots


      real(kind=8):: timetemp,wtime
      integer(ISZ):: i

      timetemp = wtime()

      if (lpltfld3d) then
        if (freqflag == ALWAYS) call parsestr("pltfld3d(1,always)")
        if (freqflag == SELDOM) call parsestr("pltfld3d(2,seldom)")
        return
      endif

      plottime = plottime + (wtime() - timetemp)

      return
      end

      subroutine divxy(pgroup,iz,ndiv,divx,divy,divvx,divvx2,divvy,divvy2,
     &                 wnpx,wnpy,itask)
      use ParticleGroupmodule
      use InMesh3d
      use InDiag
      use InPart
      use Picglb
      use Picglb3d
      use Fields3d
      type(ParticleGroup):: pgroup
      integer(ISZ):: iz,ndiv,itask
      real(kind=8):: divx(0:ndiv), divy(0:ndiv)
      real(kind=8):: divvx(0:ndiv), divvx2(0:ndiv)
      real(kind=8):: divvy(0:ndiv), divvy2(0:ndiv)
      real(kind=8):: wnpx(0:ndiv), wnpy(0:ndiv)

  Calculates rms Vx versus X and rms Vy versus Y and returns the data
  in the arrays divx and divy.  The arrays divvx, divvx2, divvy, and divvy2
  are scratch arrays that hold average vx, vx^2, vy and vy^2 versus X and Y
  respectively.
  Wnpx and wnpy are also scratch arrays that hold the numbers of points versus
  X and Y.
  The arrays divvx, divvx2, divvy, and divvy2 set up so that this routine can
  be called several times to build up data.  When itask is non-zero, the
  arrays are not cleared out.

      integer(ISZ):: i,ip,ix,iy
      real(kind=8):: ddxi,ddyi,wx,wy

  Clear everything if requested
      if (itask == 0) then
        do i=0,ndiv
          divvx(i) = 0.
          divvy(i) = 0.
          divvx2(i) = 0.
          divvx2(i) = 0.
          wnpx(i) = 0.
          wnpy(i) = 0.
        enddo
      endif

  Set inverse grid cell lengths
      ddxi = ndiv/(xplmax - xplmin)
      ddyi = ndiv/(yplmax - yplmin)

  loop through particles
      do ip=pgroup%ins(1),pgroup%ins(1)+pgroup%nps(1)-1

  if particle within z grid cell and not lost, accumulate
        if (zmesh(iz)+zbeam < pgroup%zp(ip) .and.
     &                           pgroup%zp(ip) < zmesh(iz+1)+zbeam
     &      .and. pgroup%uzp(ip) /= 0.) then
          ix = (pgroup%xp(ip) - xplmin)*ddxi
          iy = (pgroup%yp(ip) - yplmin)*ddyi
          if within x grid, accumulate
          if (ix >= 0 .and. ix < ndiv) then
            wx = (pgroup%xp(ip) - xplmin)*ddxi - ix
            wnpx(ix) = wnpx(ix) + 1. - wx
            wnpx(ix+1) = wnpx(ix+1) + wx
            divvx(ix+1)  = divvx(ix+1) +   pgroup%uxp(ip)/pgroup%uzp(ip)*wx
            divvx2(ix+1) = divvx2(ix+1) + (pgroup%uxp(ip)/pgroup%uzp(ip))**2*wx
            divvx(ix)  = divvx(ix) +   pgroup%uxp(ip)/pgroup%uzp(ip)*(1. - wx)
            divvx2(ix) = divvx2(ix) + (pgroup%uxp(ip)/pgroup%uzp(ip))**2*(1.-wx)
          endif
          if within y grid, accumulate
          if (iy >= 0 .and. iy < ndiv) then
            wy = (pgroup%yp(ip) - yplmin)*ddyi - iy
            wnpy(iy) = wnpy(iy) + 1. - wy
            wnpy(iy+1) = wnpy(iy+1) + wy
            divvy(iy)  = divvy(iy) +   pgroup%uyp(ip)/pgroup%uzp(ip)*(1.-wy)
            divvy2(iy) = divvy2(iy) + (pgroup%uyp(ip)/pgroup%uzp(ip))**2*(1.-wy)
            divvy(iy+1)  = divvy(iy+1) +   pgroup%uyp(ip)/pgroup%uzp(ip)*wy
            divvy2(iy+1) = divvy2(iy+1) + (pgroup%uyp(ip)/pgroup%uzp(ip))**2*wy
          endif
        endif
      enddo

  finally, calculate RMS quantities.
      do i=0,ndiv
        divx(i) = sqrt(max(1.e-50,divvx2(i)/dvnz(wnpx(i)) -
     &                            divvx(i)**2/dvnz(wnpx(i))))
        divy(i) = sqrt(max(1.e-50,divvy2(i)/dvnz(wnpy(i)) -
     &                            divvy(i)**2/dvnz(wnpy(i))))
      enddo

      return
      end

      subroutine multpole(lmod,nlmod,irpowmx,lcosex,lsinex,aper,xcen,ycen,
     &                    nmult,nres,tol)
      use InMesh3d
      use Multipole
      integer(ISZ):: nlmod, irpowmx
      real(kind=8):: aper, xcen, ycen 
      integer(ISZ):: lmod(nlmod) 
      logical(ISZ):: lcosex, lsinex 
      integer(ISZ):: nmult,nres
      real(kind=8):: tol

  Calculates the multipole moments as an axial function of z from the 
  electrostatic potential using a linear least squares fit to extract the 
  radial dependence.  Both the multipole expansion and the numerical 
  method employed is described in detail in Appendix A and B of HIF note 
  number 94-3 by Steve Lund.    
 
  Inputs:
 
     nlmod           number of azimuthal mode numbers to calculate moments 
     lmod(1:nlmod)   specific  azimuthal mode numbers to calculate moments 
     irpowmx         maximum radial power cutoff of moments 
     lcosex          log. flag -- calc. multipoles, cos(theta) expansion terms
     lsinex          log. flag -- calc. multipoles, sin(theta) expansion terms
     aper            aperture of moments calculation (for no apeture set to 1.)
     xcen            transverse x-coordinate moments are calculated about 
     ycen            transverse y-coordinate moments are calculated about 
     nmult           multipling factor on the number of radial points to fit
     nres            scaling for minimum radius to fit
     tol             tolerance for fitting
 
  Outputs (passed via common block) interpreter variables:
  [ All but nmom are actually calculated in the work routine multpwrk 
    called in multpole. ]
   
     nmom                  number of terms in multipole expansion 
     lazi(1:nmom)          vector of azimuthal mode numbers of moments 
     irpow(1:nmom)         vector of radial powers of moments 
     rmomcex(1:nmom,0:nzlocal)  multipoles, cos(theta) expansion terms on axial grid
     rmomsex(1:nmom,0:nzlocal)  multipoles, sin(theta) expansion terms on axial grid
  
 
  The multipole expansion is of the following schematic form 
 
                                                
 
                         
                     /--  /--
  phi(r,theta,iz) =  \--  \--   (  rmomcex(iz)_l,k cos(l theta) + 
                        l   k       
                                   rmomsex(iz)_l,k sin(l theta)   )(r/aper)^k
 
 
    where /-- denote constrained sums.  The azimuthal mode number sum l 
          \--
    is over the nlmod mode numbers specified by the vector lmod(1:nlmod), 
    and the sum over radial powers k ranges from l to irpowmx in steps of 
    2 (as a consequence of Maxwell's equations).  
          


      integer(ISZ):: l, k

      --- Set number of z grid points.
      nzmom = nzlocal

      --- find the total number of moments nmom to be calculated and 
      --- resize dynamic arrays in the common block Multipole 
      nmom = 0 
      do l = 1,nlmod
        do k = lmod(l),irpowmx,2
          nmom = nmom + 1
        enddo 
      enddo
      call gchange("Multipole",0)

      --- call work routine to do the actual multipole calculation.  This 
      --- structure is necessary to avoid addressing problems associated 
      --- with resizing multidimensional arrays
      call  multpwrk(lmod,nlmod,irpowmx,lcosex,lsinex,aper,xcen,ycen,
     &               nmult,nres,tol)

      return 
      end 

[multpole]
      subroutine multpwrk(lmod,nlmod,irpowmx,lcosex,lsinex,aper,xcen,ycen,
     &                    nmult,nres,tol)
      use Constant
      use InGen3d
      use InMesh3d
      use Picglb3d
      use Fields3d
      use Multipole
      integer(ISZ):: nlmod, irpowmx
      real(kind=8):: aper, xcen, ycen 
      integer(ISZ):: lmod(nlmod) 
      logical(ISZ):: lcosex, lsinex 
      integer(ISZ):: nmult,nres
      real(kind=8):: tol

  Work routine for the calculation of multipole moments as an axial 
  function of z.  See comments in routine multpole for variable and 
  general method descriptions.   
  Standard values for nmult=2, nres=10, tol=1.e-5.

     
      integer(ISZ):: nptsmx, nfitmx 
      parameter(nfitmx = 20, nptsmx = 500)
     
      integer(ISZ):: k, l, iz, imom, lfund, nfit, npts, ngmin 
      integer(ISZ):: np, ith, ipow, ifit, nth, ix, iy    
      real(kind=8):: ds, xsymin, ysymin, xsymax, ysymax
      real(kind=8):: dth, xc, yc, phiinter, chisq   
      real(kind=8):: rdat(nptsmx), basis(nptsmx,nfitmx)
      real(kind=8):: usvd(nptsmx,nfitmx), vsvd(nfitmx,nfitmx), wsvd(nfitmx)
      real(kind=8):: tmp(nptsmx) 
      real(kind=8):: azicosint(nptsmx), azisinint(nptsmx) 
      real(kind=8):: rmomcst(nfitmx), rmomsst(nfitmx)

      --- define constants 
      ds = sqrt((dx*dx+dy*dy)/2.)

      --- for each moment index i = 1 to nmom, specify azimuthal mode numbers 
      --- lmod(i) and radial powers irpow(i) consistent with the mode numbers 
      --- given in lmod(nlmod) and the maximum radial cutoff power irpowmx 
      imom = 0 
      do l = 1,nlmod
        do k = lmod(l),irpowmx,2
          imom = imom + 1
          lazi(imom)  = lmod(l)
          irpow(imom) = k 
        enddo 
      enddo 

      --- start moments calculation 
      imom = 0
      --- loop over each moment index with same azimuthal mode number l 
      do lfund = 1,nlmod 
        --- azimuthal mode number l   
        l = lmod(lfund)
        --- number of radial power terms nfit used in least-squares fit
        nfit = irpowmx/2 + 1 - (l+mod(irpowmx+1,2))/2
        if (nfit > nfitmx) then
          call kaboom("multpwrk: too many radial multipoles, inc. nfitmx")
        endif  
        --- calculate number radial data points npts, min grid index ngmin 
        npts  = nmult*nfit
        if (npts > nptsmx) then 
          call kaboom("multpwrk: too many radial data points, inc. nptsmx")
        endif 
        ngmin = (nres-1)*l/(2.*pi)
        if (ngmin < 2) ngmin = 2
        --- check for bound errors
        xsymin = xmmin 
        ysymin = ymmin 
        xsymax = xmmin + dx*nx
        ysymax = ymmin + dy*ny
        if (l2symtry) ysymin = -ymmax 
        if (l4symtry) then 
          xsymin = -xsymax 
          ysymin = -ysymax 
        endif  
        if ( (xcen + (npts-1+ngmin)*ds >= xsymax ) .or. 
     &       (xcen - (npts-1+ngmin)*ds <= xsymin)       ) then 
          call kaboom("multpwrk: data out of x-grid bound")
        endif   
        if ( (ycen + (npts-1+ngmin)*ds >= ysymax ) .or. 
     &       (ycen - (npts-1+ngmin)*ds <= ysymin)       ) then 
          call kaboom("multpwrk: data out of y-grid bound")
        endif   
        --- calculate npts radial data points and a npts x (kmax-l)/2 basis 
        --- array of radial powers basis for use in radial least squares fit
        do np = 1,npts
          --- assign radial data point 
          rdat(np) = (np-1+ngmin)*ds  
          --- fill (kmax-l)/2 vector column of basis matrix 
          ifit = 0
          do ipow = l,irpowmx,2
            ifit = ifit + 1
            basis(np,ifit) = rdat(np)**ipow 
          enddo 
        enddo 
        --- calculate singular value decomposition of basis matrix for 
        --- use in the linear least squares fit, first copy basis to usvd  
        do np = 1,npts 
          do ifit = 1,nfit
            usvd(np,ifit) = basis(np,ifit)  
          enddo 
        enddo 
        call dvdcmp(usvd,npts,nfit,nptsmx,nfitmx,wsvd,vsvd,tmp)
        --- loop over axial grid index iz, and calculate multipoles 
        --- with azimuthal mode number l at each index 
        do iz = 0,nzlocal
          --- azimuthal integrals of the potential at each radial data point 
          do np = 1,npts
            --- angular increments in integral 
            dth = ds/rdat(np)
            nth = int(2.*pi/dth)
            dth = 2.*pi/nth
            --- azimuthal integral 
            azicosint(np) = 0.
            azisinint(np) = 0.
            do ith = 0,nth-1
              --- transverse x and y coordinates to calculate potential  
              xc = xcen + rdat(np)*cos(dth*ith)
              yc = ycen + rdat(np)*sin(dth*ith)
              --- symmetry operations for 2-fold (l2symtry) and 
              --- 4-fold (l4symtry) symmetry about the y and x-y axes 
              if (l2symtry) then 
                if ((yc-ymmin) < 0.) yc = 2.*ymmin - yc 
              endif 
              if (l4symtry) then
                if ((xc-xmmin) < 0.) xc = 2.*xmmin - xc 
                if ((yc-ymmin) < 0.) yc = 2.*ymmin - yc 
              endif 
              --- bottom left grid indices of transverse grid box 
              --- bounding (xc,yc) 
              ix = int( ( xc - xmmin )/dx ) 
              iy = int( ( yc - ymmin )/dy ) 
              --- area interpolated potential at (xc,yc) 
              phiinter = 
     &          ( phi(ix  ,iy  ,iz)*
     &              (xmesh(ix+1) - xc       )*(ymesh(iy+1) - yc       ) +
     &            phi(ix+1,iy  ,iz)*
     &              (xc          - xmesh(ix))*(ymesh(iy+1) - yc       ) +
     &            phi(ix+1,iy+1,iz)*
     &              (xc          - xmesh(ix))*(yc          - ymesh(iy)) +
     &            phi(ix  ,iy+1,iz)*
     &              (xmesh(ix+1) - xc       )*(yc          - ymesh(iy))    )/ 
     &          (dx*dy) 
              --- add contributions for cosine and sine azimuthal integrals 
              azicosint(np) = azicosint(np) + phiinter*cos(l*dth*ith)
              azisinint(np) = azisinint(np) + phiinter*sin(l*dth*ith)
            enddo
            --- set azimuthal integrals to zero if not calculated 
            if (.not. lcosex) azicosint(np) = 0.
            if (.not. lsinex) azisinint(np) = 0.   
            --- rescale azimuthal integrals  
            azicosint(np) = (2./real(nth))*azicosint(np)
            azisinint(np) = (2./real(nth))*azisinint(np)
            if (l == 0) then 
              azicosint(np) = 0.5*azicosint(np)
              azisinint(np) = 0.
            endif  
          enddo  
          --- calculate vectors of moments, rmomcst and rmomsst for cosine 
          --- and sine expansion terms using linear least squares algorithm 
          call dvdfit(rdat,azicosint,npts,nptsmx,rmomcst,basis,nfit,nfitmx,
     &                usvd,wsvd,vsvd,tmp,tol,chisq)
          call dvdfit(rdat,azisinint,npts,nptsmx,rmomsst,basis,nfit,nfitmx,
     &                usvd,wsvd,vsvd,tmp,tol,chisq)
          --- load moments calculated in axial grid indexed vectors of 
          --- multipole moments 
          do ifit = 1,nfit
            rmomcex(imom+ifit,iz) = rmomcst(ifit)*aper**irpow(imom+ifit)
            rmomsex(imom+ifit,iz) = rmomsst(ifit)*aper**irpow(imom+ifit)  
          enddo
        enddo 
        imom = imom + nfit 
      enddo

      return 
      end