coordinate_grav_r_1D.f90

Go to the documentation of this file.

! r_coordinate for the field
!______________________________________________
module coordinate_grav_r_1D
  use phys_constant, only : nnrg, long
  use grid_parameter_1D, only : nrg, nrgin, rgin, rgmid, rgout, nrf
  implicit none
  real(long)  :: rg(0:nnrg), rginv(0:nnrg) ! radial grid points.  inv means 1/r
  real(long)  :: hrg(nnrg), hrginv(nnrg)   ! mid points.
  real(long)  :: drg(nnrg), drginv(nnrg)   ! intervals.
contains
subroutine grid_r_1D
! Grid points of r-coordinate
! --- Set up of meshes  ---  ( r(1) = dr ( not the center ) )
!
  implicit none
  real(long)  ::  drdr, drdrinv
  real(long) :: ratrr  ! ratio between subsequent step outside rgmid (or rather nrgin). \delta_{j+1} = k \delta_{j}
  real(long) :: rvdom, ratrrb, alge, dalge, error, rdet
  integer    :: nrout, ir
!
  drdr = (rgmid - rgin)/dble(nrgin)
  drdrinv = 1.0e0/drdr
  rvdom = rgout - rgmid
  nrout = nrg - nrgin
!
    drg(1:nrgin)    = drdr
    drginv(1:nrgin) = drdrinv
! 
    ratrr = 2.0e0
! 
  DO
    ratrrb = ratrr
    alge = -(ratrr**(nrout+1)-ratrr-rvdom*drdrinv*(ratrr-1.0e0))
    dalge = dble(nrout+1)*ratrr**nrout - 1.0e0 - rvdom*drdrinv
    ratrr = ratrr + alge/dalge
    error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
! 
    IF (abs(error)<=1.d-14) EXIT 
  END DO
! r-coordinate. interval dr
  DO ir = nrgin+1, nrg
    drg(ir) = ratrr*drg(ir-1)
    drginv(ir) = 1.0e0/drg(ir)
  END DO
! r-coordinate. grid ponit and mid point.
  rg(0) = 0.0e0
  rginv(0) = 0.0e0
  DO ir = 1, nrg
    rg(ir) = rg(ir-1) + drg(ir)
    hrg(ir) = (rg(ir) + rg(ir-1))*0.5e0
    rginv(ir) = 1.0e0/rg(ir)
    hrginv(ir) = 1.0e0/hrg(ir)
  END DO
!���ʓ����ɂ��邽�ߏ���ɂ‚�����
  write(6,"(' GR ', 4es23.15)") error, ratrr
  write(6,"('    ', 4es23.15)") rgout, rg(nrg)
!�‚������I���
!      WRITE(6,'(1p,2e12.4)') error, ratrr
!      WRITE(6,'(1p,2e12.4)') rgout, rg(nrg)
  rdet = (rg(nrg) - rgout)/rgout
  IF(dabs(rdet)>=1.e-10)then
    WRITE(6,*) ' bad coordinate GR ', rdet
    STOP
  END IF
end subroutine grid_r_1D
end module coordinate_grav_r_1D