00001 subroutine poisson_solver_no_dipole(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, hsinmpg, hcosmpg
00008 use weight_midpoint_grav, only : hwrtpg
00009 use make_array_1d
00010 use make_array_3d
00011 implicit none
00012
00013 integer :: im, il, nlgtmp
00014 integer :: ir, irr, it, itt, ip, ipp
00015
00016 real(long) :: pef, wok, pi4inv
00017 real(long), pointer :: sou(:,:,:),pot(:,:,:)
00018 real(long), pointer :: work1c(:,:,:)
00019 real(long), pointer :: work2c(:,:,:)
00020 real(long), pointer :: work3c(:,:,:)
00021 real(long), pointer :: work4c(:,:,:)
00022 real(long), pointer :: work1s(:,:,:)
00023 real(long), pointer :: work2s(:,:,:)
00024 real(long), pointer :: work3s(:,:,:)
00025 real(long), pointer :: work4s(:,:,:)
00026 real(long), pointer :: fac_l1(:)
00027
00028
00029 call alloc_array3d(work1c, 1, ntg, 1, nrg, 0, nlg)
00030 call alloc_array3d(work1s, 1, ntg, 1, nrg, 0, nlg)
00031 call alloc_array3d(work2c, 1, nrg, 0, nlg, 0, nlg)
00032 call alloc_array3d(work2s, 1, nrg, 0, nlg, 0, nlg)
00033 call alloc_array3d(work3c, 0, nlg, 0, nlg, 0, nrg)
00034 call alloc_array3d(work3s, 0, nlg, 0, nlg, 0, nrg)
00035 call alloc_array3d(work4c, 0, nlg, 0, nrg, 0, ntg)
00036 call alloc_array3d(work4s, 0, nlg, 0, nrg, 0, ntg)
00037 call alloc_array1d(fac_l1, 0, nlg)
00038 fac_l1(0:nlg) = 1.0d0
00039 fac_l1(1) = 0.0d0
00040
00041 work1c = 0.0
00042 work1s = 0.0
00043 do im = 0, nlg
00044 do irr = 1, nrg
00045 do itt = 1, ntg
00046 do ipp = 1, npg
00047 work1c(itt,irr,im) = work1c(itt,irr,im) &
00048 + hwrtpg(irr,itt,ipp)*sou(irr,itt,ipp)*hcosmpg(im,ipp)
00049 work1s(itt,irr,im) = work1s(itt,irr,im) &
00050 + hwrtpg(irr,itt,ipp)*sou(irr,itt,ipp)*hsinmpg(im,ipp)
00051 end do
00052 end do
00053 end do
00054 end do
00055
00056 work2c = 0.0d0
00057 work2s = 0.0d0
00058 do il = 0, nlg
00059 do im = 0, nlg
00060 pef = epsig(im)*facnmg(il,im)*fac_l1(il)
00061 do irr = 1, nrg
00062 do itt = 1, ntg
00063 wok = pef*hplmg(il,im,itt)
00064 work2c(irr,il,im) = work2c(irr,il,im) + wok*work1c(itt,irr,im)
00065 work2s(irr,il,im) = work2s(irr,il,im) + wok*work1s(itt,irr,im)
00066 end do
00067 end do
00068 end do
00069 end do
00070
00071 work3c = 0.0d0
00072 work3s = 0.0d0
00073 do ir = 0, nrg
00074 do il = 0, nlg
00075 do im = 0, nlg
00076 do irr = 1, nrg
00077 work3c(il,im,ir) = work3c(il,im,ir) &
00078 + work2c(irr,il,im)*hgfn(irr,il,ir)
00079 work3s(il,im,ir) = work3s(il,im,ir) &
00080 + work2s(irr,il,im)*hgfn(irr,il,ir)
00081 end do
00082 end do
00083 end do
00084 end do
00085
00086 work4c = 0.0d0
00087 work4s = 0.0d0
00088 do it = 0, ntg
00089 do ir = 0, nrg
00090 do im = 0, nlg
00091 do il = 0, nlg
00092 work4c(im,ir,it) = work4c(im,ir,it) &
00093 + work3c(il,im,ir)*plmg(il,im,it)
00094 work4s(im,ir,it) = work4s(im,ir,it) &
00095 + work3s(il,im,ir)*plmg(il,im,it)
00096 end do
00097 end do
00098 end do
00099 end do
00100
00101 pot = 0.0d0
00102 do ip = 0, npg
00103 do it = 0, ntg
00104 do ir = 0, nrg
00105 do im = 0, nlg
00106 pot(ir,it,ip) = pot(ir,it,ip) &
00107 + work4c(im,ir,it)*cosmpg(im,ip) &
00108 + work4s(im,ir,it)*sinmpg(im,ip)
00109 end do
00110 end do
00111 end do
00112 end do
00113
00114 pi4inv = 1.0d0/4.0d0/pi
00115 pot = - pi4inv*pot
00116
00117 deallocate(work1c)
00118 deallocate(work1s)
00119 deallocate(work2c)
00120 deallocate(work2s)
00121 deallocate(work3c)
00122 deallocate(work3s)
00123 deallocate(work4c)
00124 deallocate(work4s)
00125 deallocate(fac_l1)
00126
00127 end subroutine poisson_solver_no_dipole
00128
