00001
00002
00003 module radial_green_fn_grav
00004 use phys_constant, only : long
00005 use grid_parameter, only : nrg, nlg, rgmid, rgout
00006 use coordinate_grav_r, only : rg, rginv, hrg, hrginv
00007 use make_array_3d
00008 implicit none
00009 real(long), pointer :: hgfn(:,:,:), gfnsf(:,:,:)
00010 contains
00011 subroutine allocate_hgfn
00012 implicit none
00013
00014 call alloc_array3d(hgfn,1,nrg,0,nlg,0,nrg)
00015
00016 end subroutine allocate_hgfn
00017 subroutine allocate_hgfn_bhex
00018 implicit none
00019
00020 call alloc_array3d(hgfn,1,nrg,0,nlg,0,nrg)
00021 call alloc_array3d(gfnsf,0,nlg,0,nrg,1,4)
00022
00023 end subroutine allocate_hgfn_bhex
00024
00025
00026
00027
00028 subroutine calc_hgfn
00029 implicit none
00030 integer :: ir, irr, nn, status
00031
00032 do ir = 0, nrg
00033 do nn = 0, nlg
00034 do irr = 1, nrg
00035 hgfn(irr,nn,ir) = 0.e0
00036 end do
00037 end do
00038 end do
00039
00040 do ir = 1, nrg
00041 do irr = 1, nrg
00042 IF (hrg(irr)<rg(ir)) THEN
00043 do nn = 0, nlg
00044
00045 hgfn(irr,nn,ir) = (hrg(irr)/rg(ir))**nn/rg(ir)
00046 end do
00047 end IF
00048 IF (hrg(irr)>=rg(ir)) THEN
00049 do nn = 0, nlg
00050
00051 hgfn(irr,nn,ir) = (rg(ir)/hrg(irr))**nn/hrg(irr)
00052 end do
00053 end IF
00054 end do
00055 end do
00056
00057 ir = 0
00058 nn = 0
00059 do irr = 1, nrg
00060 hgfn(irr,nn,ir) = 1.0e0/hrg(irr)
00061 end do
00062 end subroutine calc_hgfn
00063 end module radial_green_fn_grav