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