frz_mgrid_be.F



[Lphiberz] [clampphitophimaxberz] [cond_potmgberz] [cond_potmgbezerorz] [condbndyLphiberz] [condbndymgberz] [condbndymgintberz] [getlocaldens] [getlocalepsilon] [getlocalepszfac] [getlocaltemp] [multigridberzf] [multigridberzsolve] [relaxberz] [setupregionidsbe2d]

#include top.h
 @(#) File FRZ_MGRID_BE.F,
  version $Revision: 1.38 $, $Date: 2010/06/08 21:08:18 $
 # Copyright (c) 1990-1998, The Regents of the University of California.
 # All rights reserved.  See LEGAL.LLNL for full text and disclaimer.
   This is the RZ multigrid field sovler which is part of the FRZ
   package of WARP - it includes Boltzmann electrons.
   David P. Grote, LLNL, (510)423-7194

      subroutine multigridberzf(grid,accuracy)
      use GRIDtypemodule
      use BoltzmannElectrons
      use Conductor3d
      use InGen3d, only:solvergeom,RZgeom
      use Constant, only:eps0
      use Parallel
      use DKInterp, only:interpdk,icalceps,iusing2deps
      use FRZmgrid, only:lverbose
      use multigrid_common_base, only: nb_iters,maxerr
      use ifcore
      type(GRIDtype):: grid
      real(kind=8):: accuracy,rmax,zmax
       real(kind=8),allocatable:: iondensityarray(:,:)
       real(kind=8),allocatable:: epsilon(:,:),epszfac(:,:)
      real(kind=8):: aviondensity
       real(kind=8),allocatable:: electrontemparray(:,:)
      integer(4):: fff

      integer(ISZ):: bounds(0:5)
      integer(ISZ):: ixbemin,ixbemax,izbemin,izbemax,id
      integer(ISZ):: i,j, ip1
      logical(ISZ):: lrz,lanyuseparticledensity
      integer(ISZ):: nlevels
      real(kind=8):: wr,wz
      logical(ISZ):: check,lzavparticledensity
      call for_set_fpe(FPE_M_TRAP_INV)
      fff = for_set_fpe(15_4)

      --- copy boundary positions from bound0, boundnz, and boundxy
      bounds(0) = grid%ixlbnd
      bounds(1) = grid%ixrbnd
      bounds(4) = grid%izlbnd
      bounds(5) = grid%izrbnd

      --- Check for luseparticledensity
      lanyuseparticledensity = .false.
      lzavparticledensity = .false.
      do id=1,nberegions
        if (luseparticledensity(id) > 0) then
          lanyuseparticledensity = .true.
          if (luseparticledensity(id) ==  2) lzavparticledensity = .true.
          liondensitygrid2d(id) = .true.
          --- Shift min/max to grid frame
          iondensitygrid2d%xmin = iondensitygrid2d%xmin + grid%rmin
          iondensitygrid2d%xmax = iondensitygrid2d%xmax + grid%rmin
          iondensitygrid2d%ymin = iondensitygrid2d%ymin + grid%zmin
          iondensitygrid2d%ymax = iondensitygrid2d%ymax + grid%zmin
          --- Find the extrema of all of the regions
          iondensitygrid2d%xmin = min(xbemin(id),iondensitygrid2d%xmin)
          iondensitygrid2d%xmax = max(xbemax(id),iondensitygrid2d%xmax)
          iondensitygrid2d%ymin = min(zbemin(id),iondensitygrid2d%ymin)
          iondensitygrid2d%ymax = max(zbemax(id),iondensitygrid2d%ymax)
          --- Make sure that are within the grid bounds
          iondensitygrid2d%xmin = max(grid%rmin,iondensitygrid2d%xmin)
          iondensitygrid2d%xmax = min(grid%rmax,iondensitygrid2d%xmax)
          iondensitygrid2d%ymin = max(grid%zmin,iondensitygrid2d%ymin)
          iondensitygrid2d%ymax = min(grid%zmax,iondensitygrid2d%ymax)
          --- Shift min/max back to be relative to grid mins
          iondensitygrid2d%xmin = iondensitygrid2d%xmin - grid%rmin
          iondensitygrid2d%xmax = iondensitygrid2d%xmax - grid%rmin
          iondensitygrid2d%ymin = iondensitygrid2d%ymin - grid%zmin
          iondensitygrid2d%ymax = iondensitygrid2d%ymax - grid%zmin
        endif
        logelecdenmaxscale(id) = log(electrondensitymaxscale(id))
      enddo
 
      if (lanyuseparticledensity) then
        --- Force the min's and max's to coincide with field grid values
        --- Note that the min's and max's are relative to the field grid,
        --- so for example, the location if xmin in the lab frame is
        --- iondensitygrid2d%xmin + grid%rmin. This is done so that
        --- grid%rmin etc are not needed when iondensitygrid2d is used.
        --- The mins and maxs are rounded to ensure that they are at or
        --- outside of the BE region.
        wr = floor(iondensitygrid2d%xmin/grid%dr)
        iondensitygrid2d%xmin = grid%dr*nint(wr)
        wz = floor(iondensitygrid2d%ymin/grid%dz)
        iondensitygrid2d%ymin = grid%dz*nint(wz)
        wr = ceiling(iondensitygrid2d%xmax/grid%dr)
        iondensitygrid2d%xmax = grid%dr*nint(wr)
        wz = ceiling(iondensitygrid2d%ymax/grid%dz)
        iondensitygrid2d%ymax = grid%dz*nint(wz)
        --- Now dx and dz can be set.
        iondensitygrid2d%dx = grid%dr
        iondensitygrid2d%dy = grid%dz
        --- and now the number of cells
        iondensitygrid2d%nx =
     &  nint((iondensitygrid2d%xmax-iondensitygrid2d%xmin)/iondensitygrid2d%dx)
        iondensitygrid2d%ny =
     &  nint((iondensitygrid2d%ymax-iondensitygrid2d%ymin)/iondensitygrid2d%dy)
        --- This will calculate nx and nzlocal and make sure the input is OK.
        call setupgrid2dtype(iondensitygrid2d,check)
        if (.not. check) then
           call kaboom("ERROR: luseparticledensity is being used but the input data is inconsistent")
           return
        endif
      NOTE there is no analagous coding for electrontemperaturegrid2d as the user should
      set this up in a python script -- need to set all attributes of that object.

        --- Now, copy the rho data into the iondensitygrid2d grid
        ixbemin = nint(iondensitygrid2d%xmin/grid%dr)
         ixbemax = nint(iondensitygrid2d%xmax/grid%dr)
        ixbemax = ixbemin + iondensitygrid2d%nx
        izbemin = nint(iondensitygrid2d%ymin/grid%dz)
         izbemax = nint(iondensitygrid2d%ymax/grid%dz)
        izbemax = izbemin + iondensitygrid2d%ny

        Copy rho's into iondensitygrid2d.  If luseparticledensity = 2, then
         use the z-averaged value of rho.
        if (lzavparticledensity) then
           do i = ixbemin,ixbemax
              ip1 = i+1
            NOTE grid%rho is 1-based, not 0-based, hence 1 offsets
              aviondensity = 0
              do j=izbemin+2,izbemax
                 aviondensity=aviondensity+grid%rho(i+1,j)
              enddo
              aviondensity = aviondensity+0.5*(grid%rho(ip1,izbemin+1)
     &           + grid%rho(ip1,izbemax+1))
              aviondensity = aviondensity/max(1,izbemax-izbemin)
              iondensitygrid2d%grid(i,:)= aviondensity
           enddo
        else
           iondensitygrid2d%grid(:,:) = grid%rho(ixbemin+1:ixbemax+1,izbemin+1:izbemax+1)
        endif
        iondensitygrid2d%grid(:,:) = iondensitygrid2d%grid(:,:)/eps0
      endif
 

      --- Set flag for whether to do cylindrical or Cartesian.
      lrz = (solvergeom == RZgeom)
      nlevels = grid%nlevels

      --- Create arrays for polarization dielectric function.  These
      will be set to 1 as default; if we are interpolating any species,
      we calculate the dielectric in set_polarization, in w3d_interp.F
      Note we need these whether or not there are electrons at a 
      particular location, so this is not bounded by ixbemin,ixmbemax
      or izbemin,izbemax
 
      iusing2deps = 0
      if (maxval(interpdk) > 0 .and. icalceps > 0) then
        call set_polarization(grid%rho,
     &    grid%nr,grid%nz,grid%dr,grid%dz,grid%rmin,grid%zmin)
        iusing2deps = 1
      endif
       print*,"epsilon,epszfac", epsilon(0,0),epszfac(0,0)
      The epsilon defined as above is normalized to eps0, and epszfac
      is the ratio of the epsilon for z derivs to that for x derivs.  

      call multigridberzsolve(grid%nr,grid%nz,grid%nr,grid%nz,grid%dr,grid%dz,
     &                        grid%phi,grid%rho,
     &                        bounds,
     &                        grid%rmin,grid%zmin,
     &                        grid%mgparam,nb_iters,grid%ncmax,
     &                        nlevels,maxerr,accuracy,lverbose,
     &                        grid%npre,grid%npost,
     &                        lrz,lcndbndy,laddconductor,icndbndy,
     &                        grid%bndfirst,fsdecomp)

      return
      end

[multigridberzf]
      subroutine multigridberzsolve(nx,nz,nxlocal,nzlocal,dx,dz,phi,rho,
     &                              bounds,xmminlocal,zmminlocal,
     &                              mgparam,mgiters,mgmaxiters,
     &                              mgmaxlevels,mgerror,mgtol,mgverbose,
     &                              downpasses,uppasses,
     &                              lrz,lcndbndy,laddconductor,icndbndy,
     &                              bnd,fsdecomp)
      use Subtimersfrz
      use BNDtypemodule
      use Constant
      use BoltzmannElectrons
      use DKInterp, only: iusing2deps,epsilon2d
      use Decompositionmodule
      integer(ISZ):: nx,nz
      integer(ISZ):: nxlocal,nzlocal
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1)
      real(kind=8):: rho(0:nxlocal,0:nzlocal)
      integer(ISZ),allocatable:: regionid(:,:)
      real(kind=8):: dx,dz
      integer(ISZ):: bounds(0:5)
      real(kind=8):: xmminlocal,zmminlocal
      real(kind=8):: mgparam
      integer(ISZ):: mgiters,mgmaxiters,mgmaxlevels,mgverbose
      real(kind=8):: mgerror,mgtol
      integer(ISZ):: downpasses,uppasses
      logical(ISZ):: lrz,lcndbndy,laddconductor
      integer(ISZ):: icndbndy
      type(BNDtype):: bnd
      type(Decomposition):: fsdecomp

  Use the multigrid method for solving Poisson's equation on a 3-D Cartesian
  mesh. The fieldsolver allows internal conductors with subgrid scale
  resolution and includes the Boltzmann electron term.
 
  When the grid cells are rectangular, semi-coarsening is done until the
  grid cell dimensions are roughly equal. Roughly equal means that
    2/3 dx < dz < 4/3 dx
  This keeps (max(dz,dx) - min(dz,dx))/dx < 1/3.

      integer(ISZ):: i,ii,k,ix,iz,j
      real(kind=8),allocatable:: phisave(:,:)
      real(kind=8):: xminodx
      character(72):: errline
      integer(ISZ):: allocerror
      real(kind=8):: substarttime,wtime
      if (lfrztimesubs) substarttime = wtime()

      --- The parallel version does not yet work, it is just too complicated
      --- and not needed yet.
      if (fsdecomp%nxprocs*fsdecomp%nyprocs*fsdecomp%nzprocs > 1) then
        print*,"multigridberzsolve: does not yet work in parallel"
        call kaboom("multigridberzsolve: does not yet work in parallel")
        return
      endif

      --- Make sure that the conductor data is consistent.
      --- This is presumably already taken care of elsewhere.

!$OMP PARALLEL
!$OMP&PRIVATE(ii,i,k,ix,iz)

      --- Prepare rho by dividing it by -eps0
      if (iusing2deps > 0) then
         rho = -rho/(eps0*epsilon2d%grid)
      else
         rho = -rho/eps0
      endif
      iondensity = iondensity/eps0

      --- Calculate grid min in grid cells.
      --- Note that it needs to be a real number in case xmminlocal is not an
      --- integer multiple of dx.
      xminodx = xmminlocal/dx

      allocate(phisave(0:nxlocal,0:nzlocal),stat=allocerror)
      if (allocerror /= 0) then
        print*,"multigridrzsolve: allocation error ",allocerror,
     &         ": could not allocate phisave to shape ",nxlocal,nzlocal
        call kaboom("multigridrzsolve: allocation error")
        return
      endif

      --- Create and setup array which holds the BE region ids
      allocate(regionid(0:nxlocal,0:nzlocal))
      call setupregionidsbe2d(regionid,nxlocal,nzlocal,dx,dz,xmminlocal,zmminlocal)

      --- Main multigrid v-cycle loop. Calculate error each iteration since
      --- very few iterations are done.
      mgiters = 0
      mgerror = 2.*mgtol + 1.
      do while (mgerror > mgtol .and. mgiters < mgmaxiters)
        mgiters = mgiters + 1

        --- Save current value of phi
        phisave = phi(0:nxlocal,0:nzlocal)

        --- Do one vcycle.
        call vcycleberz(0,nx,nz,nxlocal,nzlocal,dx,dz,
     &                  xminodx,xmminlocal,zmminlocal,
     &                  phi,rho,regionid,bounds,mgparam,mgmaxlevels,
     &                  downpasses,uppasses,lrz,lcndbndy,icndbndy,bnd,fsdecomp)

        --- Calculate the change in phi (the error).

        --- This line seems to create a large temporary which can
        --- cause problems when memory is close to full. So it was replaced
        --- with the explicit loop below.
        mgerror = maxval(abs(phisave - phi))

        mgerror = 0.
!$OMP DO REDUCTION(MAX:mgerror)
        do iz=0,nzlocal
          do ix=0,nxlocal
            mgerror = max(mgerror,abs(phisave(ix,iz) - phi(ix,iz)))
          enddo
        enddo
!$OMP END DO

#ifdef MPIPARALLEL
        if (fsdecomp%nxprocs*fsdecomp%nyprocs*fsdecomp%nzprocs > 1) then
          --- calculate global sorerror
          call parallelmaxrealarraycomm(mgerror,1,fsdecomp%mpi_comm)
        endif
#endif

      enddo

      --- Set boundary conditions. This is only really needed for the
      --- Dirichlet boundaries, but this is convenient to call.
      call applyboundaryconditions3d(nxlocal,0,nzlocal,1,0,1,phi,1,
     &                               bounds,.true.,.false.)

      --- Make a print out.
      if (mgverbose >= 1 .or. mgerror > mgtol) then
        if (mgerror > mgtol) then
          call remark("MultigridBE-RZ: Maximum number of iterations reached")
        endif
        write(errline,20) mgerror,mgiters
  20    format("MultigridBE-RZ: Error converged to ",1pe11.3," in ",i5," v-cycles")
        call remark(errline)
      endif

      deallocate(phisave)
      deallocate(regionid)

      --- Undo the change of rho
      if (iusing2deps > 0) then
         rho = -rho*eps0*epsilon2d%grid
      else
         rho = -rho*eps0
      endif
      iondensity = iondensity*eps0

!$OMP END PARALLEL

      if (lfrztimesubs) timemultigridberzsolve = timemultigridberzsolve +
     &                                         wtime() - substarttime

      return
      end
      RECURSIVE subroutine vcycleberz(mglevel,nx,nz,nxlocal,nzlocal,dx,dz,
     &                                xminodx,xmminlocal,zmminlocal,
     &                                phi,rho,regionid,globalbounds,
     &                                mgparam,mgmaxlevels,downpasses,uppasses,
     &                                lrz,lcndbndy,icndbndy,bnd,fsdecomp)
      use BNDtypemodule
      use BoltzmannElectrons
      use Multigrid3d_diagnostic
      use Decompositionmodule
      use RRRRRR
      integer(ISZ):: mglevel
      integer(ISZ):: nx,nz
      integer(ISZ):: nxlocal,nzlocal
      real(kind=8):: dx,dz,xminodx,xmminlocal,zmminlocal
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1),rho(0:nxlocal,0:nzlocal)
      integer(ISZ):: regionid(0:nxlocal,0:nzlocal)
      integer(ISZ):: globalbounds(0:5)
      real(kind=8):: mgparam
      integer(ISZ):: mgmaxlevels,downpasses,uppasses
      type(BNDtype):: bnd
      logical(ISZ):: lrz,lcndbndy
      integer(ISZ):: icndbndy
      type(Decomposition):: fsdecomp

  Routine that does the v-cycle for multigrid. Note that it is recursive.

      integer(ISZ):: delx,delz
      real(kind=8):: dxsqi,dzsqi
      real(kind=8),allocatable:: phicoarse(:,:),rhocoarse(:,:)
      real(kind=8),allocatable:: Lphi(:,:),Lphicoarse(:,:)
      real(kind=8),allocatable:: phicoarsesave(:,:)
      integer(ISZ),allocatable:: regionidcoarse(:,:)
      integer(ISZ):: i,iszone=1
      integer(ISZ):: nxcoarse,nzcoarse
      integer(ISZ):: nxlocalcoarse,nzlocalcoarse
      real(kind=8):: dxcoarse,dzcoarse
      real(kind=8):: dxcoarsesqi,dzcoarsesqi
      real(kind=8):: xminodxcoarse
      real(kind=8):: mgscalecoarse
      integer(ISZ):: ixproc,izproc
      integer(ISZ):: localbounds(0:5),localboundsc(0:5)
      integer(ISZ):: lxoffsetall(0:fsdecomp%nxprocs-1)
      integer(ISZ):: rxoffsetall(0:fsdecomp%nxprocs-1)
      integer(ISZ):: lzoffsetall(0:fsdecomp%nzprocs-1)
      integer(ISZ):: rzoffsetall(0:fsdecomp%nzprocs-1)
      integer(ISZ):: lxoffset,rxoffset
      integer(ISZ):: lzoffset,rzoffset
      type(Decomposition):: whosendingdown,whosendingup
      type(Decomposition):: whosendingdownc,whosendingupc
      type(Decomposition):: coarsedecomp
      integer(ISZ):: allocerror
      integer(ISZ):: i1,i2,ic,ix,iz

      delx = 1
      delz = 1
      dxsqi = 1./dx**2
      dzsqi = 1./dz**2

      localbounds = globalbounds
      localboundsc = globalbounds

#ifdef MPIPARALLEL
      lxoffsetall = 0
      rxoffsetall = 0
      lzoffsetall = 0
      rzoffsetall = 0
      call mggetexchangepes(fsdecomp,globalbounds,nx,0,nz,
     &                      lxoffsetall,rxoffsetall,
     &                      0,0,
     &                      lzoffsetall,rzoffsetall,
     &                      whosendingdown,whosendingup)
      ixproc = fsdecomp%ixproc
      izproc = fsdecomp%izproc
      if (fsdecomp%ix(ixproc) > 0)          localbounds(0) = -1
      if (fsdecomp%ix(ixproc)+nxlocal < nx) localbounds(1) = -1
      if (fsdecomp%iz(izproc) > 0)          localbounds(4) = -1
      if (fsdecomp%iz(izproc)+nzlocal < nz) localbounds(5) = -1
#endif

      --- Do initial relaxation iterations
      do i=1,downpasses
        call relaxberz(mglevel,nxlocal,nzlocal,phi,rho,regionid,
     &                 dxsqi,dzsqi,dx,dz,xminodx,
     &                 localbounds,mgparam,lrz,lcndbndy,icndbndy,bnd,
     &                 fsdecomp,whosendingdown,whosendingup)
      enddo

      --- Recurse until the bottom layer is reached.
      --- This assumes that the parameters at the coarse levels are
      --- already calculated.
      if (associated(bnd%next) .and. mglevel < mgmaxlevels) then

        nxcoarse = bnd%next%nr
        nzcoarse = bnd%next%nz
        nxlocalcoarse = bnd%next%nr
        nzlocalcoarse = bnd%next%nz
        dxcoarse = bnd%next%dr
        dzcoarse = bnd%next%dz

        dxcoarsesqi = 1./dxcoarse**2
        dzcoarsesqi = 1./dzcoarse**2
        xminodxcoarse = xminodx*nxcoarse/nx

        lxoffset = 0
        rxoffset = 0
        lzoffset = 0
        rzoffset = 0

        --- Alloate new work space
        allocate(phicoarse(-1:nxlocalcoarse+1,-1:nzlocalcoarse+1),stat=allocerror)
        if (allocerror /= 0) then
          print*,"vcycleberz: allocation error ",allocerror,
     &           ": could not allocate phicoarse to shape ",
     &           nxlocalcoarse,nzlocalcoarse
          call kaboom("vcycleberz: allocation error")
          return
        endif
        allocate(rhocoarse(0:nxlocalcoarse,0:nzlocalcoarse),stat=allocerror)
        if (allocerror /= 0) then
          print*,"vcycleberz: allocation error ",allocerror,
     &           ": could not allocate rhocoarse to shape ",
     &           nxlocalcoarse,nzlocalcoarse
          call kaboom("vcycleberz: allocation error")
          return
        endif

        allocate(Lphi(-1:nxlocal+1,-1:nzlocal+1),stat=allocerror)
        if (allocerror /= 0) then
          print*,"vcycleberz: allocation error ",allocerror,
     &           ": could not allocate Lphi to shape ",nxlocal,nzlocal
          call kaboom("vcycleberz: allocation error")
          return
        endif
        allocate(Lphicoarse(-1:nxlocalcoarse+1,-1:nzlocalcoarse+1),stat=allocerror)
        if (allocerror /= 0) then
          print*,"vcycleberz: allocation error ",allocerror,
     &           ": could not allocate Lphicoarse to shape ",
     &           nxlocalcoarse,nzlocalcoarse
          call kaboom("vcycleberz: allocation error")
          return
        endif
 
      --- Create and setup array which holds the BE region ids
        allocate(regionidcoarse(0:nxlocalcoarse,0:nzlocalcoarse),
     &           stat=allocerror)
        if (allocerror /= 0) then
          print*,"vcycleberz: allocation error ",allocerror,
     &           ": could not allocate regionidcoarse to shape ",
     &           nxlocalcoarse,nzlocalcoarse
          call kaboom("vcycleberz: allocation error")
          return
        endif
        call setupregionidsbe2d(regionidcoarse,nxlocalcoarse,nzlocalcoarse,
     &                          dxcoarse,dzcoarse,xmminlocal,zmminlocal)

        --- Calculate the coarsened phi
        phicoarse = 0.
        call restrict2d(nx,nz,nxlocal,nzlocal,phi,1,1,
     &                  nxcoarse,nzcoarse,nxlocalcoarse,nzlocalcoarse,
     &                  phicoarse,1,1,
     &                  localbounds,localboundsc,lxoffset,lzoffset,.false.)
        call cond_potmgberz(bnd%next,nxlocalcoarse,nzlocalcoarse,phicoarse,
     &                      mglevel+1)
        call applyboundaryconditions3d(nxlocalcoarse,0,nzlocalcoarse,1,0,1,
     &                                 phicoarse,
     &                                 1,localbounds,.true.,.false.)
        if (mglevel == 1) then
        if (.not. associated(ppp)) then
        allocate(ppp(-1:nxlocal+1,-1:nzlocal+1),stat=allocerror)
        endif
        ppp = phi
        endif

        --- Calculate the coarsened Lphi, putting it into rhocoarse
        call Lphiberz(nxlocal,nzlocal,dxsqi,dzsqi,dx,dz,xminodx,phi,Lphi,regionid,
     &                mglevel,localbounds,lrz,lcndbndy,icndbndy,bnd)
        Lphi(0:nxlocal,0:nzlocal) = rho - Lphi(0:nxlocal,0:nzlocal)
        call applyboundaryconditions3d(nxlocal,0,nzlocal,1,0,1,Lphi,
     &                                 1,localbounds,.true.,.false.)
        if (mglevel == 1) then
        if (.not. associated(rrr)) then
        allocate(rrr(-1:nxlocal+1,-1:nzlocal+1),stat=allocerror)
        endif
        rrr = Lphi
        endif
        call restrict2d(nx,nz,nxlocal,nzlocal,Lphi,delx,delz,
     &                  nxcoarse,nzcoarse,nxlocalcoarse,nzlocalcoarse,
     &                  rhocoarse,0,0,
     &                  localbounds,localboundsc,lxoffset,lzoffset,.false.)

        deallocate(Lphi)

        --- Calculate L(R phi), adding it into rhocoarse
        call Lphiberz(nxlocalcoarse,nzlocalcoarse,
     &                dxcoarsesqi,dzcoarsesqi,dxcoarse,dzcoarse,xminodxcoarse,
     &                phicoarse,Lphicoarse,regionidcoarse,
     &                mglevel+1,localbounds,
     &                lrz,lcndbndy,icndbndy,bnd%next)
        rhocoarse = rhocoarse + Lphicoarse(0:nxlocalcoarse,0:nzlocalcoarse)
        call cond_potmgbezerorz(bnd%next,nxlocalcoarse,nzlocalcoarse,rhocoarse,
     &                          mglevel+1,0,0)

        --- Save the current coarsened phi since it is needed to apply
        --- the corrections after relaxations.
        allocate(phicoarsesave(-1:nxlocalcoarse+1,-1:nzlocalcoarse+1),stat=allocerror)
        if (allocerror /= 0) then
          print*,"vcycleberz: allocation error ",allocerror,
     &           ": could not allocate phicoarsesave to shape ",
     &           nxlocalcoarse,nzlocalcoarse
          call kaboom("vcycleberz: allocation error")
          return
        endif
        phicoarsesave = phicoarse

        --- Continue at the next coarsest level.
        call vcycleberz(mglevel+iszone,nxcoarse,nzcoarse,
     &                  nxlocalcoarse,nzlocalcoarse,
     &                  dxcoarse,dzcoarse,xminodxcoarse,xmminlocal,zmminlocal,
     &                  phicoarse,rhocoarse,
     &                  regionidcoarse,
     &                  globalbounds,mgparam,mgmaxlevels,downpasses,uppasses,
     &                  lrz,lcndbndy,icndbndy,bnd%next,fsdecomp)

        --- Add in the correction term.
        phicoarse = phicoarse - phicoarsesave
        if (lprintmgphimaxchange) then
          print*,"Max change in phi = ",maxval(abs(phicoarse))," at MG level ",mglevel
        endif
        call expand2d(nx,nz,nxlocal,nzlocal,phi,
     &                nxcoarse,nzcoarse,nxlocalcoarse,nzlocalcoarse,phicoarse,
     &                localbounds,lxoffset,lzoffset)
        call clampphitophimaxberz(nxlocal,nzlocal,dx,dz,phi,rho,regionid)

        deallocate(Lphicoarse)
        deallocate(phicoarsesave)
        deallocate(phicoarse)
        deallocate(rhocoarse)
        deallocate(regionidcoarse)

      endif

      --- Do final relaxation passes.
      do i=1,uppasses
        call relaxberz(mglevel,nxlocal,nzlocal,phi,rho,regionid,
     &                 dxsqi,dzsqi,dx,dz,xminodx,
     &                 localbounds,mgparam,lrz,lcndbndy,icndbndy,bnd,
     &                 fsdecomp,whosendingdown,whosendingup)
        call clampphitophimaxberz(nxlocal,nzlocal,dx,dz,phi,rho,regionid)
      enddo

      return
      end

[multigridberzsolve]
      subroutine relaxberz(mglevel,nxlocal,nzlocal,phi,rho,regionid,
     &                     dxsqi,dzsqi,dx,dz,xminodx,localbounds,
     &                     mgparam,lrz,lcndbndy,icndbndy,bnd,
     &                     fsdecomp,whosendingdown,whosendingup)
      use BNDtypemodule
      use CONDtypemodule
      use BoltzmannElectrons
      use Decompositionmodule
      integer(ISZ):: mglevel,nxlocal,nzlocal
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1),rho(0:nxlocal,0:nzlocal)
      integer(ISZ):: regionid(0:nxlocal,0:nzlocal)
      real(kind=8):: dxsqi,dzsqi,dx,dz,xminodx
      integer(ISZ):: localbounds(0:5)
      real(kind=8):: mgparam
      logical(ISZ):: lrz,lcndbndy
      integer(ISZ):: icndbndy
      type(BNDtype):: bnd
      type(Decomposition):: fsdecomp
      type(Decomposition):: whosendingdown
      type(Decomposition):: whosendingup

  This routine does one pass of point SOR with even-odd (red-black)
  ordering.  It makes calls to the routines which specify internal
  conductors. The routine also allows for a bent beam-pipe.
 
  The tranverse boundaries can either be held constant, have zero normal
  derivative, or be periodic.  When BOUNDXY is zero, the boundaries are held
  constant, when 1, they have zero normal derivative, and when 2, the
  boundaries are periodic.
 
  The longitudinal boundaries can either be held constant, have zero normal
  derivative, or be periodic.  When BOUND0 or BOUNDNZ is zero, the boundaries
  are held constant, when 1, they have zero normal derivative, and when 2, the
  boundaries are periodic.

      integer(ISZ):: parity,s_parity,e_parity
      integer(ISZ):: ix,iz,ixx,izz
      integer(ISZ):: ixmin,ixmax,izmin,izmax
      integer(ISZ):: ix1
      integer(ISZ):: ic,ii,id
      real(kind=8):: rr,epsilontmp,normfac,phiav
      real(kind=8):: rhoe(0:nxlocal,0:nzlocal),rhoeote(0:nxlocal,0:nzlocal)
      real(kind=8):: denom,Lphi,expo
      real(kind=8):: boltzfac(0:nxlocal,0:nzlocal)
      real(kind=8):: iondensitytmp,electemptmp,boltzfac1,epszfactmp
      type(CONDtype):: cnd
      integer,save:: printed=0
      real(kind=8):: phiold

      --- Put desired potential onto conductors in phi array.
      call cond_potmgberz(bnd,nxlocal,nzlocal,phi,mglevel)
      call condbndymgintberz(bnd,nxlocal,nzlocal,dx,dz,phi,mglevel)
      print*, "IN RELAXBE, iondensity =", iondensity(0,0), iondensity(4,4)

      --- Set starting and ending parity.
#ifdef MPIPARALLEL
      parity = + fsdecomp%ix(fsdecomp%ixproc)
     &         + fsdecomp%iz(fsdecomp%izproc)
      s_parity = mod(parity,2)
      e_parity = mod(s_parity+1,2)
#else
      s_parity = 0
      e_parity = 1
#endif

      --- Set the loop limits, including edges when appropriate.
      ixmin = 0
      ixmax = nxlocal
      izmin = 0
      izmax = nzlocal
      if (localbounds(0) == 0) ixmin = 1
      if (localbounds(1) == 0) ixmax = nxlocal-1
      if (localbounds(4) == 0) izmin = 1
      if (localbounds(5) == 0) izmax = nzlocal-1

      --- do loop to cover even and odd points
      do parity=s_parity,e_parity,e_parity-s_parity

        --- Save values just outside conductor surfaces. Only save phi at the
        --- subgrid points which are to be used at the current level of
        --- grid refinement.
        if (lcndbndy) then
          do ic = 1,bnd%nb_conductors
            if (ic == 1) cnd = bnd%cndfirst
            if (ic > 1) cnd = cnd%next
            do ii = 1,cnd%nbbnd
              ix = cnd%jj(ii) - 1
              iz = cnd%kk(ii) - 1
              cnd%phi0(ii) = phi(ix,iz)
            enddo
          enddo
        endif

        precompute Boltzmann factor and some other arrays
        do ix = ixmin,ixmax
           do iz = izmin,izmax
              rhoe(ix,iz) = 0.
              rhoeote(ix,iz) = 0.
              if (regionid(ix,iz) > 0) then
                id = regionid(ix,iz)
                call getlocaltemp(ix,iz,dx,dz,id,electemptmp)
                call getlocaldens(ix,iz,dx,dz,id,iondensitytmp)
                call getlocalepsilon(ix,iz,dx,dz,epsilontmp)
                if (iondensitytmp .ne. 0. .and.
     &              electemptmp .ne. 0) then
                  expo=(phi(ix,iz)-plasmapotential(id))/electemptmp
                  if (.not. luselinboltz) then
                     boltzfac(ix,iz) = exp(min(expo,logelecdenmaxscale(id)))
                  else
                     boltzfac(ix,iz) = 1.+ min(expo,logelecdenmaxscale(id))
                  endif
                   print*,ix,iz, epsilontmp
                  rhoe(ix,iz) = iondensitytmp*boltzfac(ix,iz)/epsilontmp
                  rhoeote(ix,iz) = rhoe(ix,iz)/(electemptmp)
                  if (lincludeiondensityinrho(id)) then
                    rhoe(ix,iz) = rhoe(ix,iz) - iondensitytmp/epsilontmp
                  endif
                  ix1 = ixmin + mod(ixmin + iz + parity,2)
                   if (ix == ix1+2 .and. iz==izmin+1 .and. printed .lt. 10) then
                    endif
                endif
              endif
           enddo
           Normalize to z-averaged Boltzmann factor if requestted
           if (lnormtoavboltzfac) then
              normfac=nzlocal/sum(boltzfac(ix,1:nzlocal))
              boltzfac(ix,izmin:izmax) = boltzfac(ix,izmin:izmax)*normfac
              rhoe(ix,izmin:izmax)=rhoe(ix,izmin:izmax)*normfac*
     &         (1.-boltzfac(ix,izmin:izmax)/nzlocal)
              rhoeote(ix,izmin:izmax)=rhoeote(ix,izmin:izmax)*normfac
            Note this doesn't quite do the conventional thing for linearized
            Boltzmann, but it saves doing interpolation of Te a second time
            and should be approximately right whenver linearization is
            valid
           endif
        enddo

        if (lrz) then
!$OMP DO
          do iz=izmin,izmax

            ix1 = ixmin + mod(ixmin + iz + parity,2)
            do ix=ix1,ixmax,2
              call getlocalepszfac(ix,iz,dx,dz,epszfactmp)
              rr = ix + xminodx
              if (rr == 0.) then
                denom = -4.*dxsqi - 2.*epszfactmp - rhoeote(ix,iz)
                Lphi =  4.*(phi(1,iz) - phi(0,iz))*dxsqi
     &                   + (phi(0,iz-1) - 2.*phi(0,iz) + phi(0,iz+1))*epszfactmp
     &                   - rhoe(ix,iz)
              else
                denom = -2.*dxsqi - 2.*epszfactmp - rhoeote(ix,iz)
                Lphi =  ((rr-0.5)*phi(ix-1,iz) - 2*rr*phi(ix,iz) +
     &                   (rr+0.5)*phi(ix+1,iz))*dxsqi/rr
     &                + (phi(ix,iz-1) - 2.*phi(ix,iz) + phi(ix,iz+1))*epszfactmp
     &                - rhoe(ix,iz)
              endif
              phi(ix,iz) = phi(ix,iz) - mgparam*(Lphi - rho(ix,iz))/denom
            enddo
          enddo
!$OMP END DO
        else
!$OMP DO
          do iz=izmin,izmax

            ix1 = ixmin + mod(ixmin + iz + parity,2)
            do ix=ix1,ixmax,2

              call getlocalepszfac(ix,iz,dx,dz,epszfactmp)
              denom = -2.*(dxsqi+epszfactmp)-rhoeote(ix,iz)
              Lphi =  (phi(ix-1,iz  )+phi(ix+1,iz  ))*dxsqi
     &             +  (phi(ix  ,iz-1)+phi(ix  ,iz+1))*epszfactmp
     &              -  phi(ix,iz)*2.*(dxsqi+epszfactmp) - rhoe(ix,iz)
              phiold = phi(ix,iz)
              phi(ix,iz) = phi(ix,iz) - mgparam*(Lphi - rho(ix,iz))/denom
               if (ix == ix1+2 .and. iz==izmin+1 .and. printed .lt. 1000) then
                print*, "printed,denom,Lphi,phiold,phi ", printed,denom,Lphi,phiold,phi(ix,iz)
                print*, "dxsqi,epszfacdz,rhoeote",dxsqi,epszfactmp,rhoeote(ix,iz)
               printed = printed + 1
               endif
            enddo
          enddo
           print*, "summed phi =", sum(phi)
!$OMP END DO
        endif
           

        --- Apply altered difference equation to the points near the
        --- surface of the conductor boundaries.
        if (lcndbndy) then
         call condbndymgberz(bnd,parity,nxlocal,nzlocal,phi,rho,regionid,
     &                       dxsqi,dzsqi,dx,dz,xminodx,mgparam,
     &                       mglevel,lrz,icndbndy)
        endif

        --- Put desired potential onto conductors in phi array.
        call cond_potmgberz(bnd,nxlocal,nzlocal,phi,mglevel)
        call condbndymgintberz(bnd,nxlocal,nzlocal,dx,dz,phi,mglevel)

        --- set phi in the guard planes
        --- This must be done inside the loop over parities so that the
        --- guard planes are updated with the most recent values.
        call applyboundaryconditions3d(nxlocal,0,nzlocal,1,0,1,phi,
     &                                 1,localbounds,.true.,.false.)

      --- end of loop over even and odd points
      enddo

      If we are normalizing Boltzman factor to its average, then 
      we need to fix the absolute value of phi.  For now do this by
      making the average of phi = 0.  THIS ONLY MAKES SENSE FOR 
      DOUBLY PERIODIC SYSTEM.
      if (lnormtoavboltzfac) then
         phiav = sum(phi(1:nxlocal,1:nzlocal))/(nxlocal*nzlocal)
         phi = phi - phiav
      endif
       print*, "end of relaxberez, phi = ", phi(ixmax,izmax),phi(ixmax,izmax-1)
      return
      end

[multigridberzsolve] [relaxberz]
      subroutine cond_potmgberz(bnd,nxlocal,nzlocal,phi,mglevel)
      use BNDtypemodule
      use CONDtypemodule
      type(BNDtype):: bnd
      integer(ISZ):: nxlocal,nzlocal,mglevel
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1)

  Set conductor points to the desired potential. The potential is used since
  at all levels, phi is being operated on directly.

      type(CONDtype):: cnd
      integer(ISZ):: ic,ii,ix,iz

!$OMP DO
      do ic = 1,bnd%nb_conductors
        if (ic == 1) cnd = bnd%cndfirst
        if (ic > 1) cnd = cnd%next
        do ii = 1,cnd%ncond
          ix = cnd%jcond(ii) - 1
          iz = cnd%kcond(ii) - 1
          phi(ix,iz) = cnd%voltage(ii)
        enddo
!$OMP END DO
      enddo

      return
      end

[Lphiberz] [multigridberzsolve]
      subroutine cond_potmgbezerorz(bnd,nxlocal,nzlocal,u,mglevel,delt,delz)
      use BNDtypemodule
      use CONDtypemodule
      type(BNDtype):: bnd
      integer(ISZ):: nxlocal,nzlocal,mglevel,delt,delz
      real(kind=8):: u(-delt:nxlocal+delt,-delz:nzlocal+delz)

  Set data at conductor points to zero.

      type(CONDtype):: cnd
      integer(ISZ):: ic,ii,ix,iz

!$OMP DO
      do ic = 1,bnd%nb_conductors
        if (ic == 1) cnd = bnd%cndfirst
        if (ic > 1) cnd = cnd%next
        do ii = 1,cnd%ncond
          ix = cnd%jcond(ii) - 1
          iz = cnd%kcond(ii) - 1
          u(ix,iz) = 0.
        enddo
!$OMP END DO
      enddo

      return
      end

[relaxberz]
      subroutine condbndymgberz(bnd,parity,nxlocal,nzlocal,phi,rho,regionid,
     &                          dxsqi,dzsqi,dx,dz,xminodx,
     &                          mgparam,mglevel,lrz,icndbndy)
      use BNDtypemodule
      use CONDtypemodule
      use BoltzmannElectrons
      type(BNDtype):: bnd
      integer(ISZ):: parity
      integer(ISZ):: nxlocal,nzlocal,mglevel
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1), rho(0:nxlocal,0:nzlocal)
      integer(ISZ):: regionid(0:nxlocal,0:nzlocal)
      real(kind=8):: dxsqi,dzsqi,dx,dz,xminodx,mgparam
      integer(ISZ):: icndbndy
      logical(ISZ):: lrz

  Uses adjusted difference equation to enforce sub-grid level placement of 
  conductor boundaries for points near conductor surface.
 
  Temporary variables pxm, pzm, pxp, and pzp hold
  phi(i-+1)-phi(i) at minus and plus one in each direction.
  These are changed when the finite difference in the appropriate direction
  includes the boundary condition.
 
  The Cx and Cz hold the numerator of the coefficients of phi(i,j,k).
  The delx and delz hold the denominator of the coefficients of the
  full finite difference of phi.
  For icndbndy==1, these coefficients are just 1. For icndbndy==2, they
  include the dels.
 
  The ppp factor (the minimum of the dels at each grid point) is included
  here to be consistent with its inclusion in the Lphi calculation.
  Because of the structure of the solver, terms included in Lphi
  must nearly cancel the same terms here, and so must be scaled
  the same. Note that is it still unknown why the ppp factor is needed
  at all for convergence (the solver diverges with ppp=1). Also note that
  here, ppp=1 for the base level since its rho is the raw rho and does
  not include Lphi.

      type(CONDtype):: cnd
      real(kind=8):: pik,pxm,pzm,pxp,pzp,denom,rhoe,Lphi,ppp,expo
      real(kind=8):: delx,delz,Cx,Cz,rr
      real(kind=8):: delxsqi,delzsqi
      integer(ISZ):: ic,ii,i1,i2
      integer(ISZ):: ix,iz,id
      real(kind=8):: dxi,dzi,dxm,dxp,dzm,dzp
      real(kind=8):: iondensitytmp,electemptmp,epsilontmp

      --- These should really be passed in
      dxi = sqrt(dxsqi)
      dzi = sqrt(dzsqi)

      --- loop over points near surface of conductors
      do ic = 1,bnd%nb_conductors
        if (ic == 1) cnd = bnd%cndfirst
        if (ic > 1) cnd = cnd%next
        if (parity == 0) then
          --- red is even
          i1 = 1
          i2 = cnd%nbbndred
        else
          --- black is odd
          i1 = cnd%nbbndred + 1
          i2 = cnd%nbbnd
        endif

!$OMP DO
        do ii = i1,i2

          ix = cnd%jj(ii) - 1
          iz = cnd%kk(ii) - 1
          rr = ix + xminodx

          --- Set temporaries with initial values.
          pik = cnd%phi0(ii)
          if (lrz) then
            if (rr == 0.) then
              pxm = 0.
              pxp = 4.*(phi(1,iz) - pik)
              Cx = 4.
            else
              pxm = (1.-0.5/rr)*(phi(ix-1,iz) - pik)
              pxp = (1.+0.5/rr)*(phi(ix+1,iz) - pik)
              Cx = 2. !(1.-0.5/rr) + (1.+0.5/rr)
            endif
          else
            pxm = phi(ix-1,iz) - pik
            pxp = phi(ix+1,iz) - pik
            Cx = 2.
          endif
          pzm = phi(ix,iz-1) - pik
          pzp = phi(ix,iz+1) - pik
          Cz = 2.
          delx = 1.
          delz = 1.
          ppp = 1.

          dxm = cnd%dxm(ii)*dxi
          dxp = cnd%dxp(ii)*dxi
          dzm = cnd%dzm(ii)*dzi
          dzp = cnd%dzp(ii)*dzi

          --- the point lower in x is inside the conductor
          if (0. < dxm .and. dxm < 1.) then
            if (lrz) then
              if (rr > 0.) then
                pxm = (1.-0.5/rr)*(cnd%volt0xm(ii) - pik)/dxm
                Cx = Cx + (1.-0.5/rr)*(1./dxm - 1.)
                if (icndbndy == 2) delx = delx - 0.5 + 0.5*dxm
                ppp = min(ppp,dxm)
              endif
            else
              pxm = (cnd%volt0xm(ii) - pik)/dxm
              Cx = Cx - 1. + 1./dxm
              if (icndbndy == 2) delx = delx - 0.5 + 0.5*dxm
              ppp = min(ppp,dxm)
            endif
          else if (-1. < dxm .and. dxm <= 0.) then
            pxm = 0.
            Cx = Cx - 1.
            delx = delx - 0.5 + (-dxm)
            if (-dxm > 0.) then
              ppp = min(ppp,-dxm)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- the point higher in x is inside the conductor
          if (0. < dxp .and. dxp < 1.) then
            if (lrz) then
              if (rr == 0.) then
                pxp = 4.*(cnd%volt0xp(ii) - pik)/dxp
                Cx = 4./dxp
              else
                pxp = (1.+0.5/rr)*(cnd%volt0xp(ii) - pik)/dxp
                Cx = Cx + (1.+0.5/rr)*(1./dxp - 1.)
              endif
            else
              pxp = (cnd%volt0xp(ii) - pik)/dxp
              Cx = Cx - 1. + 1./dxp
            endif
            if (icndbndy == 2) delx = delx - 0.5 + 0.5*dxp
            ppp = min(ppp,dxp)
          else if (-1. < dxp .and. dxp <= 0.) then
            pxp = 0.
            Cx = Cx - 1.
            delx = delx - 0.5 + (-dxp)
            if (-dxp > 0.) then
              ppp = min(ppp,-dxp)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- the point lower in z is inside the conductor
          if (0. < dzm .and. dzm < 1.) then
            pzm = (cnd%volt0zm(ii) - pik)/dzm
            Cz = Cz - 1. + 1./dzm
            if (icndbndy == 2) delz = delz - 0.5 + 0.5*dzm
            ppp = min(ppp,dzm)
          else if (-1. < dzm .and. dzm <= 0.) then
            pzm = 0.
            Cz = Cz - 1.
            delz = delz - 0.5 + (-dzm)
            if (-dzm > 0.) then
              ppp = min(ppp,-dzm)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- the point higher in z is inside the conductor
          if (0. < dzp .and. dzp < 1.) then
            pzp = (cnd%volt0zp(ii) - pik)/dzp
            Cz = Cz - 1. + 1./dzp
            if (icndbndy == 2) delz = delz - 0.5 + 0.5*dzp
            ppp = min(ppp,dzp)
          else if (-1. < dzp .and. dzp <= 0.) then
            pzp = 0.
            Cz = Cz - 1.
            delz = delz - 0.5 + (-dzp)
            if (-dzp > 0.) then
              ppp = min(ppp,-dzp)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- calculate the new phi based on the boundary conditions
          if (ppp < 1.) then
            if (mglevel == 0) ppp = 1.
            delxsqi = dxsqi/dvnz(delx)
            delzsqi = dzsqi/dvnz(delz)
            rhoe = 0.
            denom = - Cx*delxsqi - Cz*delzsqi
            if (regionid(ix,iz) > 0) then
              id = regionid(ix,iz)
              call getlocaldens(ix,iz,dx,dz,id,iondensitytmp)
              call getlocaltemp(ix,iz,dx,dz,id,electemptmp)
              call getlocalepsilon(ix,iz,dx,dz,epsilontmp)
              if (iondensitytmp .ne. 0. .and.
     &            electemptmp .ne. 0) then
                expo = (pik - plasmapotential(id))/electemptmp
                expo = min(expo,log(electrondensitymaxscale(id)))
                rhoe = iondensitytmp*exp(expo)/epsilontmp
                denom = denom - rhoe/electemptmp
                if (lincludeiondensityinrho(id)) then
                  rhoe = rhoe - iondensitytmp/epsilontmp
                endif
              endif
            endif
            Lphi = (pxm+pxp)*delxsqi + (pzm+pzp)*delzsqi - rhoe
            phi(ix,iz) = pik - mgparam*(Lphi*ppp - rho(ix,iz))/denom
          endif
        enddo
!$OMP END DO
      enddo

      return
      end

[relaxberz]
      subroutine condbndymgintberz(bnd,nxlocal,nzlocal,dx,dz,phi,mglevel)
      use BNDtypemodule
      use CONDtypemodule
      type(BNDtype):: bnd
      integer(ISZ):: nxlocal,nzlocal,mglevel
      real(kind=8):: dx,dz
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1)

  Set the potential at points just beyond an internal Neumann boundary.

      type(CONDtype):: cnd
      integer(ISZ):: ic,ii,i1,i2,ix,iz
      real(kind=8):: dxi,dzi,dxm,dxp,dzm,dzp

      --- These should really be passed in
      dxi = 1./dx
      dzi = 1./dz

      --- loop over points near surface of conductors
      do ic = 1,bnd%nb_conductors
        if (ic == 1) cnd = bnd%cndfirst
        if (ic > 1) cnd = cnd%next

!$OMP DO
        do ii = 1,cnd%nbbnd

          ix = cnd%jj(ii) - 1
          iz = cnd%kk(ii) - 1

          dxm = cnd%dxm(ii)*dxi
          dxp = cnd%dxp(ii)*dxi
          dzm = cnd%dzm(ii)*dzi
          dzp = cnd%dzp(ii)*dzi

          if (0. > dxm .and. dxm > -1.) phi(ix-1,iz) = phi(ix,iz)
          if (0. > dxp .and. dxp > -1.) phi(ix+1,iz) = phi(ix,iz)
          if (0. > dzm .and. dzm > -1.) phi(ix,iz-1) = phi(ix,iz)
          if (0. > dzp .and. dzp > -1.) phi(ix,iz+1) = phi(ix,iz)

          if (0. == dxm) phi(ix+1,iz) = phi(ix,iz)
          if (0. == dxp) phi(ix-1,iz) = phi(ix,iz)
          if (0. == dzm) phi(ix,iz+1) = phi(ix,iz)
          if (0. == dzp) phi(ix,iz-1) = phi(ix,iz)

        enddo
      enddo
!$OMP END DO

      return
      end

[multigridberzsolve]
      subroutine Lphiberz(nxlocal,nzlocal,dxsqi,dzsqi,dx,dz,xminodx,
     &                    phi,Lphi,regionid,
     &                    mglevel,localbounds,
     &                    lrz,lcndbndy,icndbndy,bnd)
      use BNDtypemodule
      use BoltzmannElectrons
      integer(ISZ):: nxlocal,nzlocal
      real(kind=8):: dxsqi,dzsqi,dx,dz,xminodx
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1)
      real(kind=8):: Lphi(-1:nxlocal+1,-1:nzlocal+1)
      real(kind=8):: boltzfac(0:nxlocal,0:nzlocal)
      real(kind=8):: iondenstmp,normfac,epszfactmp
      integer(ISZ):: regionid(0:nxlocal,0:nzlocal)
      integer(ISZ):: mglevel,localbounds(0:5)
      logical(ISZ):: lrz,lcndbndy
      integer(ISZ):: icndbndy
      type(BNDtype):: bnd

  Calculate the Lphi on the grid including the BE term.
  Note that this is note quite the residual since it does not have the
  source term included, and so will no go to zero at convergence.

      real(kind=8):: rhoe(-1:nxlocal+1,-1:nzlocal+1),expo,epsilontmp
      real(kind=8):: rr
      integer(ISZ):: ix,iz,ic,id
      integer(ISZ):: ixmin,ixmax,izmin,izmax
      real(kind=8):: iondensitytmp,electemptmp

      --- Set the loop limits, including edges when appropriate.
      ixmin = 0
      ixmax = nxlocal
      izmin = 0
      izmax = nzlocal
      if (localbounds(0) == 0) ixmin = 1
      if (localbounds(1) == 0) ixmax = nxlocal-1
      if (localbounds(4) == 0) izmin = 1
      if (localbounds(5) == 0) izmax = nzlocal-1

      --- These are set here to avoid complaints from valgrind about
      --- uninitialized variables.
      iondensitytmp = 0.
      electemptmp = 0.

      precompute Boltzmann factor and some other arrays
      do ix = ixmin,ixmax
         do iz = izmin,izmax
            rhoe(ix,iz) = 0.
            if (regionid(ix,iz) > 0) then
               id = regionid(ix,iz)
               call getlocaltemp(ix,iz,dx,dz,id,electemptmp)
               call getlocaldens(ix,iz,dx,dz,id,iondenstmp)
               call getlocalepsilon(ix,iz,dx,dz,epsilontmp)
               if (iondensitytmp .ne. 0. .and.
     &              electemptmp .ne. 0) then
                  expo=(phi(ix,iz)-plasmapotential(id))/electemptmp
                  if (.not. luselinboltz) then
                     boltzfac(ix,iz) = exp(min(expo,logelecdenmaxscale(id)))
                  else
                     boltzfac(ix,iz) = 1.+ min(expo,logelecdenmaxscale(id))
                  endif
                  rhoe(ix,iz) = iondensitytmp*boltzfac(ix,iz)/epsilontmp
                  if (lincludeiondensityinrho(id)) then
                    rhoe(ix,iz) = rhoe(ix,iz) - iondensitytmp/epsilontmp
                  endif
               endif
            endif
         enddo
           Normalize to z-averaged Boltzmann factor if requestted
         if (lnormtoavboltzfac) then
            normfac=nzlocal/sum(boltzfac(ix,1:nzlocal))
            boltzfac(ix,izmin:izmax) = boltzfac(ix,izmin:izmax)*normfac
            rhoe(ix,izmin:izmax)=rhoe(ix,izmin:izmax)*normfac*
     &         (1.-boltzfac(ix,izmin:izmax)/nzlocal)
            Note this doesn't quite do the conventional thing for linearized
            Boltzmann, but it saves doing interpolation of Te a second time
            and should be approximately right whenver linearization is
            valid
         endif
      enddo

      Lphi = 0.
      --- Calculate the Lphi.
      if (lrz) then
        --- RZ
!$OMP DO
        do iz=izmin,izmax

          do ix=ixmin,ixmax
            call getlocalepszfac(ix,iz,dx,dz,epszfactmp)
            rr = ix + xminodx
            if (rr == 0.) then
              Lphi(ix,iz) = 4*(phi(1,iz) - phi(0,iz))*dxsqi
     &                   + (phi(0,iz-1) - 2.*phi(0,iz) + phi(0,iz+1))*epszfactmp
     &                   - rhoe(ix,iz)
            else
              Lphi(ix,iz) = ((rr-0.5)*phi(ix-1,iz) - 2*rr*phi(ix,iz) +
     &                      (rr+0.5)*phi(ix+1,iz))*dxsqi/rr
     &                   + (phi(ix,iz-1) - 2.*phi(ix,iz) + phi(ix,iz+1))*epszfactmp
     &                   - rhoe(ix,iz)
            endif

          enddo
        enddo
!$OMP END DO
      else
        --- XZ
!$OMP DO
        do iz=izmin,izmax

          do ix=ixmin,ixmax

            call getlocalepszfac(ix,iz,dx,dz,epszfactmp)

            Lphi(ix,iz) = (phi(ix-1,iz) - 2.*phi(ix,iz) + phi(ix+1,iz))*dxsqi
     &                  + (phi(ix,iz-1) - 2.*phi(ix,iz) + phi(ix,iz+1))*epszfactmp
     &                  - rhoe(ix,iz)

          enddo
        enddo
!$OMP END DO
      endif

      --- Zero the Lphi inside conductors.
      call cond_potmgbezerorz(bnd,nxlocal,nzlocal,Lphi,mglevel,1,1)

      if (lcndbndy) then
        --- Calculate the Lphi near the conductor.
        call condbndyLphiberz(bnd,nxlocal,nzlocal,phi,Lphi,regionid,dxsqi,dzsqi,
     &                        dx,dz,xminodx,mglevel,lrz,icndbndy)
      endif

      call applyboundaryconditions3d(nxlocal,0,nzlocal,1,0,1,Lphi,1,
     &                               localbounds,.false.,.true.)

      return
      end

[Lphiberz]
      subroutine condbndyLphiberz(bnd,nx,nzlocal,phi,Lphi,regionid,dxsqi,dzsqi,
     &                            dx,dz,xminodx,mglevel,lrz,icndbndy)
      use BNDtypemodule
      use CONDtypemodule
      use BoltzmannElectrons
      type(BNDtype):: bnd
      integer(ISZ):: nx,nzlocal,mglevel
      real(kind=8):: phi(-1:nx+1,-1:nzlocal+1), Lphi(-1:nx+1,-1:nzlocal+1)
      integer(ISZ):: regionid(0:nx,0:nzlocal)
      real(kind=8):: dxsqi,dzsqi,dx,dz,xminodx
      integer(ISZ):: icndbndy
      logical(ISZ):: lrz

  Uses adjusted difference equation to enforce sub-grid level placement of 
  conductor boundaries for points near conductor surface.
 
  Temporary variables pxm, pzm, pxp, and pzp hold
  phi(i-+1)-phi(i) at minus and plus one in each direction.
  These are changed when the finite difference in the appropriate direction
  includes the boundary condition.
 
  The delx and delz hold the denominator of the coefficients of the
  full finite difference of phi.
  For icndbndy==1, these coefficients are just 1. For icndbndy==2, they
  include the dels.

      type(CONDtype):: cnd
      real(kind=8):: pik,pxm,pzm,pxp,pzp,denom,rhoe,ppp,expo
      real(kind=8):: delx,delz
      real(kind=8):: delxsqi,delzsqi
      real(kind=8):: rr
      integer(ISZ):: ic,ii
      integer(ISZ):: ix,iz,id
      real(kind=8):: dxi,dzi,dxm,dxp,dzm,dzp
      real(kind=8):: iondensitytmp,electemptmp,epsilontmp

      --- These should really be passed in
      dxi = sqrt(dxsqi)
      dzi = sqrt(dzsqi)

      --- loop over points near surface of conductors
      do ic = 1,bnd%nb_conductors
        if (ic == 1) cnd = bnd%cndfirst
        if (ic > 1) cnd = cnd%next

!$OMP DO
        do ii = 1,cnd%nbbnd

          ix = cnd%jj(ii) - 1
          iz = cnd%kk(ii) - 1
          rr = ix + xminodx

          --- Set temporaries with initial values.
          pik = phi(ix,iz)
          if (lrz) then
            if (rr == 0.) then
              pxm = 0.
              pxp = 4.*(phi(1,iz) - pik)
            else
              pxm = (1.-0.5/rr)*(phi(ix-1,iz) - pik)
              pxp = (1.+0.5/rr)*(phi(ix+1,iz) - pik)
            endif
          else
            pxm = phi(ix-1,iz) - pik
            pxp = phi(ix+1,iz) - pik
          endif
          pzm = phi(ix  ,iz-1) - pik
          pzp = phi(ix  ,iz+1) - pik
          delx = 1.
          delz = 1.
          ppp = 1.

          dxm = cnd%dxm(ii)*dxi
          dxp = cnd%dxp(ii)*dxi
          dzm = cnd%dzm(ii)*dzi
          dzp = cnd%dzp(ii)*dzi

          --- the point lower in x is inside the conductor
          if (0. < dxm .and. dxm < 1.) then
            if (lrz) then
              if (rr > 0.) then
                pxm = (1.-0.5/rr)*(cnd%volt0xm(ii) - pik)/dxm
                if (icndbndy == 2) delx = delx - 0.5 + 0.5*dxm
                ppp = min(ppp,dxm)
              endif
            else
              pxm = (cnd%volt0xm(ii) - pik)/dxm
              if (icndbndy == 2) delx = delx - 0.5 + 0.5*dxm
              ppp = min(ppp,dxm)
            endif
          else if (-1. < dxm .and. dxm <= 0.) then
            pxm = 0.
            delx = delx - 0.5 + (-dxm)
            if (-dxm > 0.) then
              ppp = min(ppp,-dxm)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- the point higher in x is inside the conductor
          if (0. < dxp .and. dxp < 1.) then
            if (lrz) then
              if (rr == 0.) then
                pxp = 4.*(cnd%volt0xp(ii) - pik)/dxp
              else
                pxp = (1.+0.5/rr)*(cnd%volt0xp(ii) - pik)/dxp
              endif
            else
              pxp = (cnd%volt0xp(ii) - pik)/dxp
            endif
            if (icndbndy == 2) delx = delx - 0.5 + 0.5*dxp
            ppp = min(ppp,dxp)
          else if (-1. < dxp .and. dxp <= 0.) then
            pxp = 0.
            delx = delx - 0.5 + (-dxp)
            if (-dxp > 0.) then
              ppp = min(ppp,-dxp)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- the point lower in z is inside the conductor
          if (0. < dzm .and. dzm < 1.) then
            pzm = (cnd%volt0zm(ii) - pik)/dzm
            if (icndbndy == 2) delz = delz - 0.5 + 0.5*dzm
            ppp = min(ppp,dzm)
          else if (-1. < dzm .and. dzm <= 0.) then
            pzm = 0.
            delz = delz - 0.5 + (-dzm)
            if (-dzm > 0.) then
              ppp = min(ppp,-dzm)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- the point higher in z is inside the conductor
          if (0. < dzp .and. dzp < 1.) then
            pzp = (cnd%volt0zp(ii) - pik)/dzp
            if (icndbndy == 2) delz = delz - 0.5 + 0.5*dzp
            ppp = min(ppp,dzp)
          else if (-1. < dzp .and. dzp <= 0.) then
            pzp = 0.
            delz = delz - 0.5 + (-dzp)
            if (-dzp > 0.) then
              ppp = min(ppp,-dzp)
            else
              ppp = min(ppp,1.-1.e-9)
            endif
          endif
          --- calculate the new phi based on the boundary conditions
          if (ppp < 1.) then
            delxsqi = dxsqi/dvnz(delx)
            delzsqi = dzsqi/dvnz(delz)
            rhoe = 0.
            if (regionid(ix,iz) > 0) then
              id = regionid(ix,iz)
              call getlocaldens(ix,iz,dx,dz,id,iondensitytmp)
              call getlocaltemp(ix,iz,dx,dz,id,electemptmp)
              call getlocalepsilon(ix,iz,dx,dz,epsilontmp)
              if (iondensitytmp .ne. 0. .and.
     &            electemptmp .ne. 0) then
                expo = (pik - plasmapotential(id))/electemptmp
                expo = min(expo,log(electrondensitymaxscale(id)))
                rhoe = iondensitytmp*exp(expo)/epsilontmp
                if (lincludeiondensityinrho(id)) then
                  rhoe = rhoe - iondensitytmp/epsilontmp
                endif
              endif
            endif
            Lphi(ix,iz) = (+ (pxm+pxp)*delxsqi
     &                     + (pzm+pzp)*delzsqi
     &                     - rhoe)*ppp

            if (0. > dxm .and. dxm >= -1.) Lphi(ix-1,iz  ) = 0.
            if (0. > dxp .and. dxp >= -1.) Lphi(ix+1,iz  ) = 0.
            if (0. > dzm .and. dzm >= -1.) Lphi(ix  ,iz-1) = 0.
            if (0. > dzp .and. dzp >= -1.) Lphi(ix  ,iz+1) = 0.

            if (0. == dxm) Lphi(ix+1,iz  ) = 0.
            if (0. == dxp) Lphi(ix-1,iz  ) = 0.
            if (0. == dzm) Lphi(ix  ,iz+1) = 0.
            if (0. == dzp) Lphi(ix  ,iz-1) = 0.

          endif
        enddo
!$OMP END DO
      enddo

      return
      end

[multigridberzsolve]
      subroutine clampphitophimaxberz(nxlocal,nzlocal,dx,dz,phi,rho,regionid)
      use BoltzmannElectrons
      integer(ISZ):: nxlocal,nzlocal
      real(kind=8):: dx,dz
      real(kind=8):: phi(-1:nxlocal+1,-1:nzlocal+1),rho(0:nxlocal,0:nzlocal)
      integer(ISZ):: regionid(0:nxlocal,0:nzlocal)

      integer(ISZ):: ix,iz,id
      real(kind=8):: expo,phimax,rhomax
      real(kind=8):: iondensitytmp,electemptmp

      --- Clamp phi to phimax. The phimax is calculated assuming that
      --- the max electron density will be the max of the iondensity and 
      --- and the ion particle density (rho).
      do iz=0,nzlocal
        do ix=0,nxlocal
          if (regionid(ix,iz) > 0) then
            id = regionid(ix,iz)
            call getlocaldens(ix,iz,dx,dz,id,iondensitytmp)
            call getlocaltemp(ix,iz,dx,dz,id,electemptmp)
            if (iondensitytmp == 0.) cycle
            if (luseparticledensity(id) > 0) then
              rhomax = max(electrondensitymaxscale(id),-rho(ix,iz)/iondensitytmp)
            else
              rhomax = electrondensitymaxscale(id)
            endif
            phimax = plasmapotential(id)+electemptmp*log(rhomax)
            if (phi(ix,iz) > phimax) phi(ix,iz) = phimax
          endif
        enddo
      enddo

      return
      end

[multigridberzsolve]
      subroutine setupregionidsbe2d(regionid,nxlocal,nzlocal,dx,dz,
     &                              xmminlocal,zmminlocal)
      use BoltzmannElectrons
      integer(ISZ):: nxlocal,nzlocal
      integer(ISZ):: regionid(0:nxlocal,0:nzlocal)
      real(kind=8):: dx,dy,dz,xmminlocal,zmminlocal

  Sets id of regions where the Boltzmann-electron distribution is used.

      integer(ISZ):: id,ixmin,ixmax,izmin,izmax,ix,iz

      regionid = 0
      do id=1,nberegions
        ixmin = (max(0.,xbemin(id) - xmminlocal))/dx
        ixmax = min(1.0*nxlocal,(xbemax(id) - xmminlocal)/dx)
        izmin = (max(0.,zbemin(id) - zmminlocal))/dz
        izmax = min(1.0*nzlocal,(zbemax(id) - zmminlocal)/dz)
        --- Constrain the mins and maxes to be within the BE region.
        if (xmminlocal+ixmin*dx < xbemin(id)) ixmin = ixmin + 1
        if (xmminlocal+ixmax*dx > xbemax(id)) ixmax = ixmax - 1
        if (zmminlocal+izmin*dz < zbemin(id)) izmin = izmin + 1
        if (zmminlocal+izmax*dz > zbemax(id)) izmax = izmax - 1
        regionid(ixmin:ixmax,izmin:izmax) = id
      enddo

      return
      end

[Lphiberz] [clampphitophimaxberz] [condbndyLphiberz] [condbndymgberz] [relaxberz]
      subroutine getlocaldens(ix,iz,dx,dz,id,iondensitytmp)
      generate interpolated local iondensity and electron temperature
      use BoltzmannElectrons
      integer(ISZ):: ix,iz,id
      real(kind=8):: dx,dz,reixx,reizz
      real(kind=8):: iondensitytmp

      real(kind=8):: xx,zz,wxx,wzz
      integer(ISZ):: ixx,izz

       print*,"ix,iz,id,liondensitygrid2d ", ix,iz,id,liondensitygrid2d
      iondensitytmp = 0.
      if (liondensitygrid2d(id)) then
        xx = ix*dx
        zz = iz*dz
        reixx = (xx - iondensitygrid2d%xmin)/iondensitygrid2d%dx
        reizz = (zz - iondensitygrid2d%ymin)/iondensitygrid2d%dy
        ixx = reixx
        izz = reizz
        wxx = reixx-ixx
        wzz = reizz-izz
        if (ixx+1 > iondensitygrid2d%nx) then
          ixx = iondensitygrid2d%nx-1
          wxx = 1.
        endif
        if (izz+1 > iondensitygrid2d%ny) then
          izz = iondensitygrid2d%ny-1
          wzz = 1.
        endif
        iondensitytmp =
     &       iondensitygrid2d%grid(ixx  ,izz  )*(1.-wxx)*(1.-wzz) +
     &       iondensitygrid2d%grid(ixx+1,izz  )*    wxx *(1.-wzz) +
     &       iondensitygrid2d%grid(ixx  ,izz+1)*(1.-wxx)*    wzz  +
     &       iondensitygrid2d%grid(ixx+1,izz+1)*    wxx *    wzz
      else if (iondensity(id) .ne. 0.) then
        iondensitytmp = iondensity(id)
      endif

      return
      end

[Lphiberz] [clampphitophimaxberz] [condbndyLphiberz] [condbndymgberz] [relaxberz]
      subroutine getlocaltemp(ix,iz,dx,dz,id,electemptmp)
      generate interpolated local electemp and electron temperature
      use BoltzmannElectrons
      integer(ISZ):: ix,iz,id
      real(kind=8):: dx,dz,reixx,reizz
      real(kind=8):: electemptmp

      real(kind=8):: xx,zz,wxx,wzz
      integer(ISZ):: ixx,izz

       print*,"ix,iz,id,electemperature ", ix,iz,id,electrontemperature(id)
      electemptmp = 0.
      if (lelectrontemperaturegrid2d(id)) then
        xx = ix*dx
        zz = iz*dz
        reixx = (xx - electrontemperaturegrid2d%xmin)/electrontemperaturegrid2d%dx
        reizz = (zz - electrontemperaturegrid2d%ymin)/electrontemperaturegrid2d%dy
        ixx = reixx
        izz = reizz
        wxx = reixx-ixx
        wzz = reizz-izz
        if (ixx+1 > electrontemperaturegrid2d%nx) then
           ixx = electrontemperaturegrid2d%nx-1
           wxx = 1.
        endif
        if (izz+1 > electrontemperaturegrid2d%ny) then
           izz = electrontemperaturegrid2d%ny-1
           wzz = 1.
        endif
        electemptmp =
     &       electrontemperaturegrid2d%grid(ixx  ,izz  )*(1.-wxx)*(1.-wzz) +
     &       electrontemperaturegrid2d%grid(ixx+1,izz  )*    wxx *(1.-wzz) +
     &       electrontemperaturegrid2d%grid(ixx  ,izz+1)*(1.-wxx)*    wzz  +
     &       electrontemperaturegrid2d%grid(ixx+1,izz+1)*    wxx *    wzz

      else if (electrontemperature(id) .ne. 0.) then
        electemptmp = electrontemperature(id)
      endif

      return
      end

[Lphiberz] [condbndyLphiberz] [condbndymgberz] [relaxberz]
      subroutine getlocalepsilon(ix,iz,dx,dz,epsilontmp)
      generate interpolated local epsilon2d and electron temperature
      use DKInterp, only:interpdk,icalceps,epsilon2d,iusing2deps
      use BoltzmannElectrons
      integer(ISZ):: ix,iz,id
      real(kind=8):: dx,dz,reixx,reizz
      real(kind=8):: epsilontmp

      real(kind=8):: xx,zz,wxx,wzz
      integer(ISZ):: ixx,izz

       print*,"ix,iz,id,lepsilongrid2d ", ix,iz,id,lepsilon2d2d
      epsilontmp = 0.
      if (iusing2depsɬ) then
        xx = ix*dx
        zz = iz*dz
        reixx = (xx - epsilon2d%xmin)/epsilon2d%dx
        reizz = (zz - epsilon2d%ymin)/epsilon2d%dy
        ixx = reixx
        izz = reizz
        wxx = reixx-ixx
        wzz = reizz-izz
        if (ixx+1 > epsilon2d%nx) then
           ixx = epsilon2d%nx-1
           wxx = 1.
        endif
        if (izz+1 > epsilon2d%ny) then
           izz = epsilon2d%ny-1
           wzz = 1.
        endif
        epsilontmp =
     &       epsilon2d%grid(ixx  ,izz  )*(1.-wxx)*(1.-wzz) +
     &       epsilon2d%grid(ixx+1,izz  )*    wxx *(1.-wzz) +
     &       epsilon2d%grid(ixx  ,izz+1)*(1.-wxx)*    wzz  +
     &       epsilon2d%grid(ixx+1,izz+1)*    wxx *    wzz
      else
        epsilontmp = 1.
      endif

      return
      end


[Lphiberz] [relaxberz]
      subroutine getlocalepszfac(ix,iz,dx,dz,epszfactmp)
      generate interpolated local epsilon2d and electron temperature
      use DKInterp, only:interpdk,icalceps,epszfac2d,iusing2deps
      use BoltzmannElectrons
      integer(ISZ):: ix,iz,id
      real(kind=8):: dx,dz,reixx,reizz
      real(kind=8):: epszfactmp

      real(kind=8):: xx,zz,wxx,wzz
      integer(ISZ):: ixx,izz

       print*,"ix,iz,id,iusing2deps ", ix,iz,id,lepsilon2d2d
      epszfactmp = 0.
      if (iusing2depsɬ) then
        xx = ix*dx
        zz = iz*dz
        reixx = (xx - epszfac2d%xmin)/epszfac2d%dx
        reizz = (zz - epszfac2d%ymin)/epszfac2d%dy
        ixx = reixx
        izz = reizz
        wxx = reixx-ixx
        wzz = reizz-izz
        if (ixx+1 > epszfac2d%nx) then 
           ixx = epszfac2d%nx-1
           wxx = 1.
        endif
        if (izz+1 > epszfac2d%ny) then 
           izz = epszfac2d%ny-1
           wzz = 1.
        endif
        epszfactmp =
     &       epszfac2d%grid(ixx  ,izz  )*(1.-wxx)*(1.-wzz) +
     &       epszfac2d%grid(ixx+1,izz  )*    wxx *(1.-wzz) +
     &       epszfac2d%grid(ixx  ,izz+1)*(1.-wxx)*    wzz  +
     &       epszfac2d%grid(ixx+1,izz+1)*    wxx *    wzz
      else
        epszfactmp = 1.
      endif
      epszfactmp = epszfactmp/(dz*dz)

      return
      end