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