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

00001 ! r_coordinate for the field
00002 !______________________________________________
00003 module coordinate_grav_r
00004   use phys_constant, only : nnrg, long
00005   use grid_parameter, only : nrg, nrgin, rgin, rgmid, rgout, nrf
00006   implicit none
00007   real(long)  :: rg(0:nnrg), rginv(0:nnrg) ! radial grid points.  inv means 1/r
00008   real(long)  :: hrg(nnrg), hrginv(nnrg)   ! mid points.
00009   real(long)  :: drg(nnrg), drginv(nnrg)   ! intervals.
00010 contains
00011 subroutine grid_r
00012 ! Grid points of r-coordinate
00013 ! --- Set up of meshes  ---  ( r(1) = dr ( not the center ) )
00014 !
00015   implicit none
00016   real(long)  ::  drdr, drdrinv
00017   real(long) :: ratrr  ! ratio between subsequent step outside rgmid (or rather nrgin). \delta_{j+1} = k \delta_{j}
00018   real(long) :: rvdom, ratrrb, alge, dalge, error, rdet, rdetf
00019   real(long) :: rfdom, reso_min, drmin, drmininv, ratrf, ratrfb
00020   integer    :: nrout, ir, itry
00021 !
00022   drdr =  1.0d0/dble(nrf)
00023   drdrinv = 1.0d0/drdr
00024   rvdom = rgout - rgmid
00025   nrout = nrg - nrgin
00026 !
00027   drg(1:nrgin+1)    = drdr
00028   drginv(1:nrgin+1) = drdrinv
00029 ! 
00030   if (rgin.ge.1.0d-13.and.rgin.lt.1.0d0) then 
00031     rfdom = 1.0d0 - rgin
00032     ratrf = 0.1d0
00033     do 
00034       ratrfb = ratrf
00035       alge = -(ratrf**nrf - 1.0d0 - rfdom*drdrinv*(ratrf-1.0e0))
00036       dalge = dble(nrf)*ratrf**(nrf-1) - rfdom*drdrinv
00037       ratrf = ratrf + alge/dalge
00038       error = 2.d0*(ratrf - ratrfb)/(ratrf + ratrfb)
00039       if (abs(error) <= 1.0d-14) exit
00040     end do
00041     if (ratrf.gt.1.0d0) stop 'ratrf not good'
00042     drg(nrf) = drdr
00043     drginv(nrf) = 1.0d0/drg(nrf)
00044     do ir = nrf-1, 1, -1
00045       drg(ir) = drg(ir+1)*ratrf
00046       drginv(ir) = 1.0d0/drg(ir)
00047     end do
00048     do ir = nrf+1, nrgin
00049       drg(ir) = drg(nrf)
00050       drginv(ir) = 1.0d0/drg(ir)
00051     end do
00052   end if
00053 !
00054 do itry = 0, 3 
00055   if (itry.eq.0) ratrr = 2.0e-01   ! good for finer grid
00056   if (itry.eq.1) ratrr = 2.0e0   ! good for finer grid
00057   if (itry.eq.2) ratrr = 10.0e0  ! good for coarser grid
00058   if (itry.eq.3) ratrr = 50.0e0  ! good for coarsest grid
00059 ! 
00060   DO
00061     ratrrb = ratrr
00062     alge = -(ratrr**(nrout+1)-ratrr-rvdom*drdrinv*(ratrr-1.0e0))
00063     dalge = dble(nrout+1)*ratrr**nrout - 1.0e0 - rvdom*drdrinv
00064     ratrr = ratrr + alge/dalge
00065     error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00066 ! 
00067     IF (abs(error)<=1.d-14) EXIT 
00068   END DO
00069 ! r-coordinate. interval dr
00070   DO ir = nrgin+1, nrg
00071     drg(ir) = ratrr*drg(ir-1)
00072     drginv(ir) = 1.0e0/drg(ir)
00073   END DO
00074 ! r-coordinate. grid ponit and mid point.
00075   rg(0) = rgin
00076   if (rg(0).le.1.0d-14) then 
00077     rginv(0) = 0.0e0
00078   else 
00079     rginv(0) = 1.0d0/rg(0)
00080   end if
00081   DO ir = 1, nrg
00082     rg(ir) = rg(ir-1) + drg(ir)
00083     hrg(ir) = (rg(ir) + rg(ir-1))*0.5e0
00084     rginv(ir) = 1.0e0/rg(ir)
00085     hrginv(ir) = 1.0e0/hrg(ir)
00086 !      write(6,*)ir,rg(ir)
00087   END DO
00088 !      WRITE(6,'(1p,2e12.4)') error, ratrr
00089 !      WRITE(6,'(1p,2e12.4)') rgout, rg(nrg)
00090   rdet  = (rg(nrg) - rgout)/rgout
00091   rdetf = 0.0d0
00092   if (rgin.lt.1.0d0) rdetf =  rg(nrf) - 1.0d0
00093   IF(dabs(rdet)< 1.e-10) exit
00094   IF(dabs(rdet)>=1.e-10)then
00095     WRITE(6,*) ' bad coordinate itry =', itry
00096     WRITE(6,*) ' bad coordinate GR ', rdet, ratrr
00097     WRITE(6,*) ' bad coordinate FLUID ', rdetf, ratrf
00098     if (itry.eq.3) STOP
00099   END IF
00100 end do
00101 end subroutine grid_r
00102 !
00103 subroutine grid_r_bhex(char_bh)
00104 ! Grid points of r-coordinate
00105 ! --- Set up of meshes  ---  ( r(1) = dr ( not the center ) )
00106 !
00107   use def_bh_parameter, only : mass_pBH
00108   implicit none
00109   real(long)  ::  drdr, drdrinv
00110   real(long) :: ratrr  ! ratio between subsequent step outside rgmid (or rather nrgin). \delta_{j+1} = k \delta_{j}
00111   real(long) :: rvdom, ratrrb, alge, dalge, error, rdet, rdetf
00112   real(long) :: rfdom, reso_min, drmin, drmininv, ratrf, ratrfb
00113   integer    :: nrout, ir, itry, nrf_bh
00114   real(long) :: dr1, dr1inv, tolerant_val = 5.0d-14
00115   integer    :: nr_count, count
00116   character(len=3) :: char_bh
00117 !
00118   drdr =  1.0d0/dble(nrf)
00119   drdrinv = 1.0d0/drdr
00120 !
00121   drg(1:nrgin+1)    = drdr
00122   drginv(1:nrgin+1) = drdrinv
00123 ! 
00124 !  if (rgin.ge.1.0d-13.and.rgin.lt.1.0d0) then 
00125   if (rgin.lt.1.0d0) then 
00126 !
00127 do itry = 0, 3 
00128   if (itry.eq.0) ratrr = 1.01d0 ! good for finer grid
00129   if (itry.eq.1) ratrr = 1.1d0   ! good for finer grid
00130   if (itry.eq.2) ratrr = 2.0d0  ! good for coarser grid
00131   if (itry.eq.3) ratrr = 10.0d0  ! good for coarsest grid
00132 ! 
00133 !    dr1 = rgin/(0.5d0*dble(nrf))
00134     dr1 = rgin/(0.75d0*dble(nrf))
00135     if (char_bh.eq.'pBH') dr1 = 0.5d0*mass_pBH/(0.75d0*dble(nrf))
00136     dr1inv = 1.0d0/dr1
00137     rvdom = 1.0d0 - rgin
00138     nrout = nrf
00139 !
00140   count = 0
00141   DO
00142     count = count + 1
00143     ratrrb = ratrr
00144     alge = -(ratrr**(nrout)-1.0d0-rvdom*dr1inv*(ratrr-1.0d0))
00145     dalge = dble(nrout)*ratrr**(nrout-1) - rvdom*dr1inv
00146 !test    alge = -((ratrr**nrout-1.0d0)/(ratrr-1.0d0)-rvdom*dr1inv)
00147 !test    dalge = (dble(nrout-1)*ratrr**nrout-dble(nrout)*ratrr**(nrout-1)+1.0d0) &
00148 !test    &     / (ratrr-1.0d0)**2
00149     ratrr = ratrr + alge/dalge
00150     error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00151 ! 
00152     IF (abs(error)<=1.d-15) EXIT 
00153     IF (count.ge.1000) exit
00154   END DO
00155 !test  DO
00156 !test    count = count + 1
00157 !test    ratrrb = ratrr
00158 !test    alge = -(dlog(ratrr**nrout-1.0d0)-dlog(ratrr-1.0d0)-dlog(rvdom*dr1inv))
00159 !test    dalge = dble(nrout)*ratrr**(nrout-1)/(ratrr**nrout-1.0d0) &
00160 !test    &     - 1.0d0/(ratrr-1.0d0)
00161 !test    ratrr = ratrr + alge/dalge
00162 !test    error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00163 !test! 
00164 !test    IF (abs(error)<=1.d-15) EXIT 
00165 !test    IF (count.ge.1000) exit
00166 !test  END DO
00167 !
00168   IF (count.ge.1000) cycle
00169 ! r-coordinate. interval dr
00170     drg(1) = dr1
00171     drginv(1) = 1.0d0/dr1
00172     do ir = 2, nrf
00173       drg(ir) = ratrr*drg(ir-1)
00174       drginv(ir) = 1.0d0/drg(ir)
00175     end do
00176 !
00177     rg(0) = rgin
00178     DO ir = 1, nrf
00179       rg(ir) = rg(ir-1) + drg(ir)
00180     end do
00181 !
00182   rdet  = (rg(nrf) - 1.0d0)/1.0d0
00183   rdetf = 0.0d0
00184   if (rgin.lt.1.0d0) rdetf =  rg(nrf) - 1.0d0
00185   IF(dabs(rdet)< 2.e-14) then
00186     WRITE(6,*) 'good coordinate itry =', itry
00187     WRITE(6,*) 'good coordinate GR ', rdet, ratrr
00188     WRITE(6,*) 'good coordinate FLUID ', rdetf, ratrf
00189     exit
00190   end if
00191   IF(dabs(rdet)>=2.e-14)then
00192     WRITE(6,*) ' bad coordinate itry =', itry
00193     WRITE(6,*) ' bad coordinate GR ', rdet, ratrr
00194     WRITE(6,*) ' bad coordinate FLUID ', rdetf, ratrf
00195     if (itry.eq.3) STOP
00196   END IF
00197 !
00198 end do
00199 !
00200 ! region III
00201     drdr =  drg(nrgin)
00202     drdrinv = 1.0d0/drdr
00203     rvdom = 2.0d0*rgmid
00204     nrout = nrf
00205     ratrr = 2.0d0
00206 !
00207     DO
00208       ratrrb = ratrr
00209       alge = -(ratrr**(nrout+1)-ratrr-rvdom*drdrinv*(ratrr-1.0d0))
00210       dalge = dble(nrout+1)*ratrr**nrout - 1.0d0 - rvdom*drdrinv
00211       ratrr = ratrr + alge/dalge
00212       error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00213 ! 
00214       IF (abs(error)<=1.d-15) EXIT 
00215     END DO
00216 ! r-coordinate. interval dr
00217     DO ir = nrgin+1, nrgin+nrf
00218       drg(ir) = ratrr*drg(ir-1)
00219       drginv(ir) = 1.0d0/drg(ir)
00220     END DO
00221   end if
00222 !
00223 do itry = 0, 3 
00224   if (itry.eq.0) ratrr = 1.05d-00 ! good for finer grid
00225   if (itry.eq.1) ratrr = 2.0d0   ! good for finer grid
00226   if (itry.eq.2) ratrr = 10.0d0  ! good for coarser grid
00227   if (itry.eq.3) ratrr = 50.0d0  ! good for coarsest grid
00228 ! 
00229   drdr =  1.0d0/dble(nrf)
00230   drdrinv = 1.0d0/drdr
00231   rvdom = rgout - rgmid
00232   nrout = nrg - nrgin
00233 !  if (rgin.ge.1.0d-13.and.rgin.lt.1.0d0) then 
00234   if (rgin.lt.1.0d0) then 
00235     drdr =  drg(nrgin+nrf)
00236     drdrinv = 1.0d0/drdr
00237     rvdom = rgout - 3.0d0*rgmid
00238     nrout = nrg - (nrgin + nrf)
00239   end if
00240 !
00241   DO
00242     ratrrb = ratrr
00243 !test    alge = -(ratrr**(nrout+1)-ratrr-rvdom*drdrinv*(ratrr-1.0d0))
00244 !test    dalge = dble(nrout+1)*ratrr**nrout - 1.0d0 - rvdom*drdrinv
00245     alge = -(dble(nrout+1)*log(ratrr)-log(ratrr+rvdom*drdrinv*(ratrr-1.0d0)))
00246     dalge = dble(nrout+1)/ratrr - (1.0d0+rvdom*drdrinv)/ &
00247     &                             (ratrr+rvdom*drdrinv*(ratrr-1.0d0))
00248     ratrr = ratrr + alge/dalge
00249     error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00250 ! 
00251     IF (abs(error)<=1.d-16) EXIT 
00252   END DO
00253 !
00254 ! r-coordinate. interval dr
00255   nr_count = nrgin+1
00256 !  if (rgin.ge.1.0d-13.and.rgin.lt.1.0d0) nr_count = nrgin+nrf+1
00257   if (rgin.lt.1.0d0) nr_count = nrgin+nrf+1
00258   DO ir = nr_count, nrg
00259     drg(ir) = ratrr*drg(ir-1)
00260     drginv(ir) = 1.0d0/drg(ir)
00261   END DO
00262 ! r-coordinate. grid ponit and mid point.
00263   rg(0) = rgin
00264   if (rg(0).le.1.0d-14) then 
00265     rginv(0) = 0.0d0
00266   else 
00267     rginv(0) = 1.0d0/rg(0)
00268   end if
00269   DO ir = 1, nrg
00270     rg(ir) = rg(ir-1) + drg(ir)
00271     hrg(ir) = (rg(ir) + rg(ir-1))*0.5d0
00272     rginv(ir) = 1.0d0/rg(ir)
00273     hrginv(ir) = 1.0d0/hrg(ir)
00274 !      write(6,*)ir,rg(ir)
00275   END DO
00276 !      WRITE(6,'(1p,2e12.4)') error, ratrr
00277 !      WRITE(6,'(1p,2e12.4)') rgout, rg(nrg)
00278   rdet  = (rg(nrg) - rgout)/rgout
00279   rdetf = 0.0d0
00280   if (rgin.lt.1.0d0) rdetf =  rg(nrf) - 1.0d0
00281   IF(dabs(rdet)< tolerant_val) then
00282     WRITE(6,*) 'good coordinate itry =', itry
00283     WRITE(6,*) 'good coordinate GR ', rdet, ratrr
00284     WRITE(6,*) 'good coordinate FLUID ', rdetf, ratrf
00285     exit
00286   end if
00287   IF(dabs(rdet)>=tolerant_val)then
00288     WRITE(6,*) ' bad coordinate itry =', itry
00289     WRITE(6,*) ' bad coordinate GR ', rdet, ratrr
00290     WRITE(6,*) ' bad coordinate FLUID ', rdetf, ratrf
00291     if (itry.eq.3) STOP
00292   END IF
00293 end do
00294 ! testtest
00295 !do ir = 0, nrg
00296 !write(6,*) rg(ir)
00297 !end do
00298 !stop
00299 !
00300 end subroutine grid_r_bhex
00301 !
00302 end module coordinate_grav_r
/home/uryu/Cocal_svn_checkinout/Cocal/code/Module/coordinate_grav_r.f90
