Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Module/radial_green_fn_grav_bhex

00001 ! --- Computation of the radial green's function hgfn. ---
00002 ! ________________________________________________________
00003 module radial_green_fn_grav_bhex_all
00004   use phys_constant, only : long
00005   use grid_parameter, only : nrg, nlg, rgin, rgout
00006   use coordinate_grav_r, only : rg, rginv, hrg, hrginv
00007   use make_array_3d
00008   implicit none
00009   real(long), pointer  ::  hgfn_nb(:,:,:), gfnsf_nb(:,:,:)
00010   real(long), pointer  ::  hgfn_di(:,:,:), gfnsf_di(:,:,:)
00011   real(long), pointer  ::  hgfn_dd(:,:,:), gfnsf_dd(:,:,:)
00012   real(long), pointer  ::  hgfn_nd(:,:,:), gfnsf_nd(:,:,:)
00013 contains
00014 subroutine allocate_hgfn_bhex_all
00015   implicit none
00016 !
00017   call alloc_array3d(hgfn_nb,1,nrg,0,nlg,0,nrg)
00018   call alloc_array3d(hgfn_di,1,nrg,0,nlg,0,nrg)
00019   call alloc_array3d(hgfn_dd,1,nrg,0,nlg,0,nrg)
00020   call alloc_array3d(hgfn_nd,1,nrg,0,nlg,0,nrg)
00021   call alloc_array3d(gfnsf_nb,0,nlg,0,nrg,1,4)
00022   call alloc_array3d(gfnsf_di,0,nlg,0,nrg,1,4)
00023   call alloc_array3d(gfnsf_dd,0,nlg,0,nrg,1,4)
00024   call alloc_array3d(gfnsf_nd,0,nlg,0,nrg,1,4)
00025 !
00026 end subroutine allocate_hgfn_bhex_all
00027 ! ----------------------------------------------
00028 ! Subroutine
00029 ! hgfn(irr,il,ir) --> for r' < r
00030 !                     for r  < r'
00031 ! _nb --> no boundary
00032 ! _di --> BH dirichlet, Asymptotics no boundary
00033 ! _dd --> BH dirichlet, Asymptotics Dirichlet
00034 ! _nd --> BH Neumann,   Asymptotics Dirichlet
00035 ! N.B.  rg(0) = rgin should not equal to zero
00036 ! ----------------------------------------------
00037 subroutine calc_hgfn_bhex_all
00038   implicit none
00039   integer     ::  ir, irr, il
00040   real(long)  ::  bhrad,  bhradinv,  bhsq
00041   real(long)  ::  radext, radextinv, radsq
00042   real(long)  ::  fac1, fac2, fac3, dil1, dil2
00043 !
00044   do ir  = 0, nrg
00045     do il = 0, nlg
00046       do irr = 1, nrg
00047         hgfn_nb(irr,il,ir) = 0.d0
00048         hgfn_di(irr,il,ir) = 0.d0
00049         hgfn_dd(irr,il,ir) = 0.d0
00050         hgfn_nd(irr,il,ir) = 0.d0
00051       end do
00052     end do
00053   end do
00054 ! 
00055 ! --  1. No boundary Green's fn.
00056 ! --  2. Dirichlet at rgin
00057 ! --  3. Drichlet at rgin and rgout
00058 ! --  4. Neumann at rgin and Dirichlet at rgout
00059 !
00060   bhrad     = rgin
00061   bhradinv  = 1.0d0/rgin
00062   radext    = rgout
00063   radextinv = 1.0d0/rgout
00064 !
00065   do ir  = 0, nrg
00066     do irr = 1, nrg
00067       if (hrg(irr).lt.rg(ir)) then
00068         do il  = 0, nlg
00069 ! -- 1.
00070           hgfn_nb(irr,il,ir) = hrg(irr)**il/rg(ir)**(il+1)
00071 ! -- 2.
00072           hgfn_di(irr,il,ir) =(hrg(irr)**il/rg(ir)**(il+1) &
00073           &                  - bhrad**(2*il+1)/(hrg(irr)*rg(ir))**(il+1))
00074 ! -- 3.
00075           fac1 = 1.0d0/(1.0d0 - (bhrad*radextinv)**(2*il+1))
00076           fac2 = (bhrad*radextinv)**il*radextinv
00077           hgfn_dd(irr,il,ir) = fac1*fac2 &
00078           &   *((hrg(irr)*bhradinv)**il - (bhrad*hrginv(irr))**(il+1)) &
00079           &   *((radext*rginv(ir))**(il+1) - (rg(ir)*radextinv)**il)
00080 ! -- 4.
00081           dil1 = dble(il)/dble(il+1)
00082           fac1 = 1.0d0/(1.0d0 + dil1*(bhrad*radextinv)**(2*il+1))
00083           fac2 = (bhrad*radextinv)**il*radextinv
00084           hgfn_nd(irr,il,ir) = fac1*fac2 &
00085           &   *((hrg(irr)*bhradinv)**il + dil1*(bhrad*hrginv(irr))**(il+1)) &
00086           &   *((radext*rginv(ir))**(il+1) - (rg(ir)*radextinv)**il)
00087         end do
00088       end if
00089       if (hrg(irr).ge.rg(ir)) then
00090         do il  = 0, nlg
00091 ! -- 1.
00092           hgfn_nb(irr,il,ir) = rg(ir)**il/hrg(irr)**(il+1)
00093 ! -- 2.
00094           hgfn_di(irr,il,ir) =(rg(ir)**il/hrg(irr)**(il+1) &
00095           &                  - bhrad**(2*il+1)/(hrg(irr)*rg(ir))**(il+1))
00096 ! -- 3.
00097           fac1 = 1.0d0/(1.0d0 - (bhrad*radextinv)**(2*il+1))
00098           fac2 = (bhrad*radextinv)**il*radextinv
00099           hgfn_dd(irr,il,ir) = fac1*fac2 &
00100           &   *((rg(ir)*bhradinv)**il - (bhrad*rginv(ir))**(il+1)) &
00101           &   *((radext*hrginv(irr))**(il+1) - (hrg(irr)*radextinv)**il)
00102 ! -- 4.
00103           dil1 = dble(il)/dble(il+1)
00104           fac1 = 1.0d0/(1.0d0 + dil1*(bhrad*radextinv)**(2*il+1))
00105           fac2 = (bhrad*radextinv)**il*radextinv
00106           hgfn_nd(irr,il,ir) = fac1*fac2 &
00107           &   *((rg(ir)*bhradinv)**il + dil1*(bhrad*rginv(ir))**(il+1)) &
00108           &   *((radext*hrginv(irr))**(il+1) - (hrg(irr)*radextinv)**il)
00109         end do
00110       end if
00111     end do
00112   end do
00113 !
00114 ! --- Green's functions for the surface term (TIMES R^2)
00115 !
00116 !     ia --> 1 for d phi/dr condition at rgin (inward normal). 
00117 !            2 for   phi(r) condition at rgin (inward normal).
00118 !            3 for d phi/dr at rgout (outward normal). 
00119 !            4 for   phi(r) at rgout (outward normal).
00120 !     gfnsf_(l,r,ia)
00121 !
00122   bhrad     = rgin
00123   bhradinv  = 1.0d0/rgin
00124   radext    = rgout
00125   radextinv = 1.0d0/rgout
00126   bhsq      = bhrad**2
00127   radsq     = radext**2
00128 !
00129   do ir = 0, nrg
00130     il = 0
00131 ! -- 1.
00132     gfnsf_nb(il,ir,1) = - bhrad*(bhrad*rginv(ir))
00133     gfnsf_nb(il,ir,2) = 0.0d0
00134     gfnsf_nb(il,ir,3) = radext
00135     gfnsf_nb(il,ir,4) = - dble(il+1)
00136 ! -- 2.
00137     gfnsf_di(il,ir,1) = - 0.0d0
00138     gfnsf_di(il,ir,2) = - bhrad*rginv(ir)
00139     gfnsf_di(il,ir,3) = radext*(1.0d0-bhrad*rginv(ir))
00140     gfnsf_di(il,ir,4) = - (1.0d0-bhrad*rginv(ir))
00141 ! -- 3.
00142     fac1 = 1.0d0/(1.0d0 - (bhrad*radextinv)**(2*il+1))
00143     fac2 = dble(2*il+1)*bhrad**(il-1)*radextinv**(il+1)
00144     fac3 = dble(2*il+1)*bhrad**il*radextinv**(il+2)
00145     gfnsf_dd(il,ir,1) = - 0.0d0
00146     gfnsf_dd(il,ir,2) = - fac1*fac2*(radext*rginv(ir) - 1.0d0)*bhsq
00147     gfnsf_dd(il,ir,3) = 0.0d0
00148     gfnsf_dd(il,ir,4) = - fac1*fac3*(1.0d0 - bhrad*rginv(ir))*radsq
00149 ! -- 4.
00150     dil1 = dble(il)/dble(il+1)
00151     dil2 = dble(2*il+1)/dble(il+1)
00152     fac1 = 1.0d0/(1.0d0 + dil1*(bhrad*radextinv)**(2*il+1))
00153     fac2 = dil2*(bhrad*radextinv)**il*radextinv
00154     fac3 = dble(2*il+1)*bhrad**il*radextinv**(il+2)
00155     gfnsf_nd(il,ir,1) = - fac1*fac2*(radext*rginv(ir) - 1.0d0)*bhsq
00156     gfnsf_nd(il,ir,2) = - 0.0d0
00157     gfnsf_nd(il,ir,3) = 0.0d0
00158     gfnsf_nd(il,ir,4) = - fac1*fac3 &
00159     &                 *(1.0d0 + dil1*bhrad*rginv(ir))*radsq
00160 !
00161     do il = 1, nlg
00162 ! -- 1.
00163       gfnsf_nb(il,ir,1) = - bhrad*(bhrad*rginv(ir))**(il+1)
00164       gfnsf_nb(il,ir,2) = - dble(il)*(bhrad*rginv(ir))**(il+1)
00165       gfnsf_nb(il,ir,3) = radext*(rg(ir)*radextinv)**il
00166       gfnsf_nb(il,ir,4) = - dble(il+1)*(rg(ir)*radextinv)**il
00167 ! -- 2.
00168       gfnsf_di(il,ir,1) = - 0.0d0
00169       gfnsf_di(il,ir,2) = &
00170       &    - (2.0d0*dble(il)+1.0d0)*(bhrad*rginv(ir))**(il+1)
00171       gfnsf_di(il,ir,3) = radext*(bhrad/radext)**il &
00172       &    *((rg(ir)*bhradinv)**il-(bhrad*rginv(ir))**(il+1))
00173       gfnsf_di(il,ir,4) = - dble(il+1)*(bhrad/radext)**il &
00174       &    *((rg(ir)*bhradinv)**il-(bhrad*rginv(ir))**(il+1))
00175 ! -- 3.
00176       fac1 = 1.0d0/(1.0d0 - (bhrad*radextinv)**(2*il+1))
00177       fac2 = dble(2*il+1)*bhrad**(il-1)*radextinv**(il+1)
00178       fac3 = dble(2*il+1)*bhrad**il*radextinv**(il+2)
00179       gfnsf_dd(il,ir,1) = - 0.0d0
00180       gfnsf_dd(il,ir,2) = - fac1*fac2 &
00181       &    *((radext*rginv(ir))**(il+1)-(rg(ir)*radextinv)**il)*bhsq
00182       gfnsf_dd(il,ir,3) = 0.0d0
00183       gfnsf_dd(il,ir,4) = - fac1*fac3 &
00184       &    *((rg(ir)*bhradinv)**il-(bhrad*rginv(ir))**(il+1))*radsq
00185 ! -- 4.
00186       dil1 = dble(il)/dble(il+1)
00187       dil2 = dble(2*il+1)/dble(il+1)
00188       fac1 = 1.0d0/(1.0d0 + dil1*(bhrad*radextinv)**(2*il+1))
00189       fac2 = dil2*(bhrad*radextinv)**il*radextinv
00190       fac3 = dble(2*il+1)*bhrad**il*radextinv**(il+2)
00191       gfnsf_nd(il,ir,1) = - fac1*fac2 &
00192       &    *((radext*rginv(ir))**(il+1)-(rg(ir)*radextinv)**il)*bhsq
00193       gfnsf_nd(il,ir,2) = - 0.0d0
00194       gfnsf_nd(il,ir,3) = 0.0d0
00195       gfnsf_nd(il,ir,4) = - fac1*fac3 &
00196       &    *((rg(ir)*bhradinv)**il+dil1*(bhrad*rginv(ir))**(il+1))*radsq
00197 !
00198     end do
00199   end do
00200 !
00201 end subroutine calc_hgfn_bhex_all
00202 end module radial_green_fn_grav_bhex_all
/home/uryu/Cocal_svn_checkinout/Cocal/code/Module/radial_green_fn_grav_bhex_all.f90
