00001 subroutine gauge_Coulomb
00002   use grid_parameter, only : nrg, ntg, npg
00003   use def_emfield, only : vaxd, vayd, vazd
00004   use interface_poisson_solver
00005   use make_array_3d
00006   use interface_grgrad_4th_gridpoint
00007   use interface_grgrad_midpoint_type0
00008   implicit none
00009 
00010   real(8), pointer :: sou(:,:,:), gauge(:,:,:)
00011   real(8) :: dvaxdx, dvaxdy, dvaxdz, dvaydx, dvaydy, dvaydz, 
00012             dvazdx, dvazdy, dvazdz
00013   real(8) :: dgadx, dgady, dgadz
00014   integer :: irg, itg, ipg
00015 
00016   call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00017   call alloc_array3d(gauge,0,nrg,0,ntg,0,npg)
00018 
00019 
00020 
00021   do ipg = 1, npg
00022     do itg = 1, ntg
00023       do irg = 1, nrg
00024 
00025         call grgrad_midpoint_type0(vaxd,dvaxdx,dvaxdy,dvaxdz,irg,itg,ipg)
00026         call grgrad_midpoint_type0(vayd,dvaydx,dvaydy,dvaydz,irg,itg,ipg)
00027         call grgrad_midpoint_type0(vazd,dvazdx,dvazdy,dvazdz,irg,itg,ipg)
00028         sou(irg,itg,ipg) = -(dvaxdx + dvaydy + dvazdz)
00029 
00030       end do
00031     end do
00032   end do
00033 
00034   call poisson_solver(sou,gauge)
00035 
00036   do ipg = 0, npg
00037     do itg = 0, ntg
00038       do irg = 0, nrg
00039 
00040         call grgrad_4th_gridpoint(gauge,dgadx,dgady,dgadz,irg,itg,ipg)
00041 
00042         vaxd(irg,itg,ipg) = vaxd(irg,itg,ipg) + dgadx
00043         vayd(irg,itg,ipg) = vayd(irg,itg,ipg) + dgady
00044         vazd(irg,itg,ipg) = vazd(irg,itg,ipg) + dgadz
00045 
00046       end do
00047     end do
00048   end do
00049 
00050   deallocate(sou)
00051   deallocate(gauge)
00052 end subroutine gauge_Coulomb