00001 subroutine poisson_solver(sou, pot)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, npg, ntg, nlg
00004 use radial_green_fn_grav, only : hgfn
00005 use legendre_fn_grav, only : plmg, hplmg, facnmg, epsig
00006 use trigonometry_grav_theta
00007 use trigonometry_grav_phi, only : sinmpg, cosmpg
00008 use weight_midpoint_grav, only : hwrtpg
00009 use make_array_3d
00010 implicit none
00011 integer :: im, il
00012 integer :: ir, irr, it, itt, ip, ipp
00013 real(long), pointer :: sou(:,:,:),pot(:,:,:)
00014 real(long) :: pef, wok, pi4inv
00015 real(long), pointer :: work1c(:,:,:)
00016 real(long), pointer :: work2c(:,:,:)
00017 real(long), pointer :: work3c(:,:,:)
00018 real(long), pointer :: work4c(:,:,:)
00019 real(long), pointer :: work1s(:,:,:)
00020 real(long), pointer :: work2s(:,:,:)
00021 real(long), pointer :: work3s(:,:,:)
00022 real(long), pointer :: work4s(:,:,:)
00023
00024 call alloc_array3d(work1c, 1, nrg, 1, ntg, 0, nlg)
00025 call alloc_array3d(work1s, 1, nrg, 1, ntg, 0, nlg)
00026 call alloc_array3d(work2c, 1, nrg, 0, nlg, 0, nlg)
00027 call alloc_array3d(work2s, 1, nrg, 0, nlg, 0, nlg)
00028 call alloc_array3d(work3c, 0, nrg, 0, nlg, 0, nlg)
00029 call alloc_array3d(work3s, 0, nrg, 0, nlg, 0, nlg)
00030 call alloc_array3d(work4c, 0, nrg, 0, ntg, 0, nlg)
00031 call alloc_array3d(work4s, 0, nrg, 0, ntg, 0, nlg)
00032
00033 work1c = 0.0
00034 work1s = 0.0
00035 do im = 0, nlg
00036 do ipp = 1, npg
00037 work1c(1:nrg,1:ntg,im) = work1c(1:nrg,1:ntg,im) &
00038 & + 0.5*hwrtpg(1:nrg,1:ntg,ipp)*sou(1:nrg,1:ntg,ipp) &
00039 & *(cosmpg(im,ipp)+cosmpg(im,ipp-1))
00040 work1s(1:nrg,1:ntg,im) = work1s(1:nrg,1:ntg,im) &
00041 & + 0.5*hwrtpg(1:nrg,1:ntg,ipp)*sou(1:nrg,1:ntg,ipp) &
00042 & *(sinmpg(im,ipp)+sinmpg(im,ipp-1))
00043 end do
00044 end do
00045
00046 work2c = 0.0d0
00047 work2s = 0.0d0
00048 do il = 0, nlg
00049 do im = 0, nlg
00050 pef = epsig(im)*facnmg(il,im)
00051 do itt = 1, ntg
00052 wok = pef*hplmg(il,im,itt)
00053 work2c(1:nrg,il,im) = work2c(1:nrg,il,im) + wok*work1c(1:nrg,itt,im)
00054 work2s(1:nrg,il,im) = work2s(1:nrg,il,im) + wok*work1s(1:nrg,itt,im)
00055 end do
00056 end do
00057 end do
00058
00059 work3c = 0.0d0
00060 work3s = 0.0d0
00061 do ir = 0, nrg
00062 do il = 0, nlg
00063 do irr = 1, nrg
00064 work3c(ir,il,0:nlg) = work3c(ir,il,0:nlg) &
00065 & + work2c(irr,il,0:nlg)*hgfn(irr,il,ir)
00066 work3s(ir,il,0:nlg) = work3s(ir,il,0:nlg) &
00067 & + work2s(irr,il,0:nlg)*hgfn(irr,il,ir)
00068 end do
00069 end do
00070 end do
00071
00072 work4c = 0.0d0
00073 work4s = 0.0d0
00074 do it = 0, ntg
00075 do im = 0, nlg
00076 do il = 0, nlg
00077 work4c(0:nrg,it,im) = work4c(0:nrg,it,im) &
00078 & + work3c(0:nrg,il,im)*plmg(il,im,it)
00079 work4s(0:nrg,it,im) = work4s(0:nrg,it,im) &
00080 & + work3s(0:nrg,il,im)*plmg(il,im,it)
00081 end do
00082 end do
00083 end do
00084
00085 pot = 0.0d0
00086 do ip = 0, npg
00087 do im = 0, nlg
00088 pot(0:nrg,0:ntg,ip) = pot(0:nrg,0:ntg,ip) &
00089 & + work4c(0:nrg,0:ntg,im)*cosmpg(im,ip) &
00090 & + work4s(0:nrg,0:ntg,im)*sinmpg(im,ip)
00091 end do
00092 end do
00093
00094 pi4inv = 1.0d0/4.0d0/pi
00095 pot = - pi4inv*pot
00096
00097 deallocate(work1c)
00098 deallocate(work1s)
00099 deallocate(work2c)
00100 deallocate(work2s)
00101 deallocate(work3c)
00102 deallocate(work3s)
00103 deallocate(work4c)
00104 deallocate(work4s)
00105
00106 end subroutine poisson_solver
00107
