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, nrg_1, r_surf
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 !
00303 subroutine grid_r_bns
00304 
00305   implicit none
00306   real(long)  ::  drdr, drdrinv
00307   real(long) :: ratrr  ! ratio between subsequent step outside rgmid (or rather nrgin). \delta_{j+1} = k \delta_{j}
00308   real(long) :: rvdom, ratrrb, alge, dalge, error, rdet, rdetf, rdet1
00309   real(long) :: rfdom, reso_min, drmin, drmininv, ratrf, ratrfb
00310   integer    :: nrout, ir, itry, nrf_bh
00311   real(long) :: dr1, dr1inv, tolerant_val = 5.0d-14
00312   real(long) :: dr2, dr2inv
00313   integer    :: nr_count, count, nr3
00314   character(len=3) :: char_bh
00315 
00316 !
00317 ! Region I r<=rgmid  for mpt3.
00318   drdr =  1.0d0/dble(nrf)
00319   drdrinv = 1.0d0/drdr
00320 !
00321   drg(1:nrgin+1)    = drdr
00322   drginv(1:nrgin+1) = drdrinv
00323 !
00324 ! the following if block only for mpt1,2:   Sets up the grid for r<3rgmid
00325 ! for mpt3 we have only the above lines i.e.  r<rgmid
00326 !
00327   if (rgin==0.0d0) then
00328     write(6,*) "REGION S, I, II, III:"
00329 !
00330     drdr = r_surf/dble(nrf)
00331     drdrinv = 1.0d0/drdr
00332 
00333 !   Region S
00334     drg(1:nrf)    = drdr
00335     drginv(1:nrf) = drdrinv
00336 
00337 !   Region II
00338     dr2 =  1.0e0/dble(nrg_1)            ! nrg_1 plays the role of nrf in BBH
00339     drg(nrg_1+1:nrgin)    = dr2
00340     drginv(nrg_1+1:nrgin) = 1.0e0/dr2
00341 
00342     if(r_surf<1.0d0) then
00343       if(nrf==nrg_1) then
00344         write(6,*) "nrf==nrg_1 but r_surf<1.0...exiting"
00345         stop
00346       end if
00347       do itry = 0, 1
00348         if (itry.eq.0) ratrr = 2.0e-01 ! good for finer grid
00349         if (itry.eq.1) ratrr = 2.0e0   ! good for finer grid
00350 !
00351         dr1 = drg(nrf)                 ! drg(nrf+1)= k drg(nrf)
00352         dr1inv = 1.0d0/dr1
00353         rvdom = 1.0d0 - r_surf
00354         nrout = nrg_1 - nrf
00355 
00356         DO
00357           ratrrb = ratrr
00358           alge = -(ratrr**(nrout+1)-ratrr-rvdom*dr1inv*(ratrr-1.0e0))
00359           dalge = dble(nrout+1)*ratrr**nrout - 1.0e0 - rvdom*dr1inv
00360           ratrr = ratrr + alge/dalge
00361           error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00362 !
00363           IF (abs(error)<=1.d-14) EXIT
00364         END DO
00365         write(6,*) 'region I: ratrr=', ratrr
00366 !
00367 !       r-coordinate. interval dr
00368         do ir = nrf+1, nrg_1
00369           drg(ir) = ratrr*drg(ir-1)
00370           drginv(ir) = 1.0e0/drg(ir)
00371         end do
00372 
00373         rg(0) = rgin
00374         DO ir = 1, nrgin
00375           rg(ir) = rg(ir-1) + drg(ir)
00376 !          write(6,*) ir, rg(ir), drg(ir)
00377         end do
00378 !
00379         rdet  = (rg(nrgin) - rgmid)/rgmid
00380         rdetf = 0.0d0
00381         if (rgin.lt.1.0d0) rdet1 =  rg(nrg_1) - 1.0d0
00382         if (rgin.lt.1.0d0) rdetf =  rg(nrf) - r_surf
00383         IF(dabs(rdet)< 1.e-10) then
00384           WRITE(6,*) 'good coordinate itry =', itry
00385           WRITE(6,*) 'good coordinate GR : rdet1, rdet, ratrr ', rdet1, rdet, ratrr
00386           WRITE(6,*) 'good coordinate FLUID : rdetf', rdetf
00387           exit
00388         end if
00389         IF(dabs(rdet)>=1.e-10)then
00390           WRITE(6,*) ' bad coordinate itry =', itry
00391           WRITE(6,*) ' bad coordinate GR : rdet1, rdet, ratrr ', rdet1, rdet, ratrr
00392           WRITE(6,*) ' bad coordinate FLUID : rdetf ', rdetf
00393           if (itry.eq.1) STOP
00394         END IF
00395       end do
00396     else
00397       write(6,*) "Region S = Region I"
00398       if(nrf.ne.nrg_1)  then
00399         write(6,*) "nrf is not equal to nrg_1, but r_surf=1...exiting"
00400         stop
00401       end if
00402     end if                   !  r_surf < 1.0
00403 !
00404 !   region III.
00405     drdr =  drg(nrgin)       !  In bhex first interval in region III: drg(nrgin+1)= k drg(nrgin)
00406                              !  If we want to change drdr we must also modify
00407                              !  LINE (*)
00408     drdrinv = 1.0d0/drdr
00409     rvdom = 2.0d0*rgmid
00410     nr3 = nrg_1               !  plays the role of nrf in BBH
00411     ratrr = 2.0e0
00412 !
00413     DO
00414       ratrrb = ratrr
00415       alge = -(ratrr**(nr3+1)-ratrr-rvdom*drdrinv*(ratrr-1.0e0))
00416       dalge = dble(nr3+1)*ratrr**nr3 - 1.0e0 - rvdom*drdrinv
00417       ratrr = ratrr + alge/dalge
00418       error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00419 !
00420       IF (abs(error)<=1.d-14) EXIT
00421     END DO
00422 
00423     write(6,*) 'region III: ratrr=', ratrr, "    dr1=", drdr, "    dr(nrf)=", drg(nrf)
00424 
00425 !   r-coordinate. interval dr:   drg(nrgin+1)= k drdr
00426     drg(nrgin+1)    = ratrr*drdr
00427     drginv(nrgin+1) = 1.0e0/drg(nrgin+1)
00428 
00429     DO ir = nrgin+2, nrgin+nr3
00430       drg(ir) = ratrr*drg(ir-1)               !   LINE (*)
00431       drginv(ir) = 1.0e0/drg(ir)
00432     END DO
00433 !
00434     rg(0) = rgin
00435     DO ir = 1, nrgin+nr3
00436       rg(ir) = rg(ir-1) + drg(ir)
00437     end do
00438     write(6,*)  "rg(nrgin+nr3)=", rg(nrgin+nr3), "      4rgmid=", 4.0d0*rgmid
00439   end if                         !  rgin==0
00440 !
00441 !  Now for r > 3rgmid if mpt1,2
00442 !  or      r > rgmid  if mpt3
00443 !
00444 do itry = 0, 3
00445   if (itry.eq.0) ratrr = 2.0e-01 ! good for finer grid
00446   if (itry.eq.1) ratrr = 2.0e0   ! good for finer grid
00447   if (itry.eq.2) ratrr = 10.0e0  ! good for coarser grid
00448   if (itry.eq.3) ratrr = 50.0e0  ! good for coarsest grid
00449 !
00450   drdr =  1.0d0/dble(nrf)
00451   drdrinv = 1.0d0/drdr
00452   rvdom = rgout - rgmid
00453   nrout = nrg - nrgin
00454   if (rgin==0.0d0) then
00455     drdr =  drg(nrgin+nr3)
00456     drdrinv = 1.0d0/drdr
00457     rvdom = rgout - 3.0d0*rgmid
00458     nrout = nrg - (nrgin + nr3)
00459   end if
00460 !
00461   DO
00462     ratrrb = ratrr
00463     alge = -(ratrr**(nrout+1)-ratrr-rvdom*drdrinv*(ratrr-1.0e0))
00464     dalge = dble(nrout+1)*ratrr**nrout - 1.0e0 - rvdom*drdrinv
00465     ratrr = ratrr + alge/dalge
00466     error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00467 !
00468     IF (abs(error)<=1.d-14) EXIT
00469   END DO
00470 
00471   write(6,*) 'region IV: ratrr=', ratrr
00472 
00473 ! r-coordinate. interval dr
00474   nr_count = nrgin+1
00475   if (rgin==0.0d0) nr_count = nrgin+nr3+1
00476   DO ir = nr_count, nrg
00477     drg(ir) = ratrr*drg(ir-1)
00478     drginv(ir) = 1.0e0/drg(ir)
00479   END DO
00480 ! r-coordinate. grid ponit and mid point.
00481   rg(0) = rgin
00482   if (rg(0).le.1.0d-14) then
00483     rginv(0) = 0.0e0
00484   else
00485     rginv(0) = 1.0d0/rg(0)
00486   end if
00487   DO ir = 1, nrg
00488     rg(ir) = rg(ir-1) + drg(ir)
00489     hrg(ir) = (rg(ir) + rg(ir-1))*0.5e0
00490     rginv(ir) = 1.0e0/rg(ir)
00491     hrginv(ir) = 1.0e0/hrg(ir)
00492 !      write(6,*)ir,rg(ir)
00493   END DO
00494 !      WRITE(6,'(1p,2e12.4)') error, ratrr
00495 !      WRITE(6,'(1p,2e12.4)') rgout, rg(nrg)
00496   rdet  = (rg(nrg) - rgout)/rgout
00497   rdetf = 0.0d0
00498   if (rgin.lt.1.0d0) rdet1 =  rg(nrg_1) - 1.0d0
00499   if (rgin.lt.1.0d0) rdetf =  rg(nrf) - r_surf
00500   IF(dabs(rdet)< 1.e-10) then
00501     WRITE(6,*) 'good coordinate itry =', itry
00502     WRITE(6,*) 'good coordinate GR : rdet1, rdet, ratrr ', rdet1, rdet, ratrr
00503     WRITE(6,*) 'good coordinate FLUID : rdetf ', rdetf
00504     exit
00505   end if
00506   IF(dabs(rdet)>=1.e-10)then
00507     WRITE(6,*) ' bad coordinate itry =', itry
00508     WRITE(6,*) ' bad coordinate GR : rdet1, rdet, ratrr', rdet1, rdet, ratrr
00509     WRITE(6,*) ' bad coordinate FLUID : rdetf ', rdetf
00510     if (itry.eq.3) STOP
00511   END IF
00512 end do
00513 !
00514 end subroutine grid_r_bns
00515 !
00516 !
00517 subroutine grid_r_bns_const
00518 
00519   implicit none
00520   real(long)  ::  drdr, drdrinv
00521   real(long) :: ratrr  ! ratio between subsequent step outside rgmid (or rather nrgin). \delta_{j+1} = k \delta_{j}
00522   real(long) :: rvdom, ratrrb, alge, dalge, error, rdet, rdetf, rdet1
00523   real(long) :: rfdom, reso_min, drmin, drmininv, ratrf, ratrfb
00524   integer    :: nrout, ir, itry, nrf_bh, nrbns
00525   real(long) :: dr1, dr1inv, tolerant_val = 5.0d-14
00526   real(long) :: dr2, dr2inv, rbns
00527   integer    :: nr_count, count, nr3
00528   character(len=3) :: char_bh
00529 
00530 !
00531 ! Region I r<=rgmid  for mpt3.
00532   drdr =  1.0d0/dble(nrf)
00533   drdrinv = 1.0d0/drdr
00534 !
00535   drg(1:nrgin+1)    = drdr
00536   drginv(1:nrgin+1) = drdrinv
00537 !
00538 ! the following if block only for mpt1,2:   Sets up the grid for r<3rgmid
00539 ! for mpt3 we have only the above lines i.e.  r<rgmid
00540 !
00541   if (rgin==0.0d0) then
00542     nrbns = 7*nrf
00543     rbns  = 7.0d0
00544     rvdom = rgout - rbns
00545     nrout = nrg - nrbns
00546 !
00547     if (nrout<0)  then
00548       write(6,*) "It must be nrg > 7*nrf....exiting"
00549       stop
00550     end if
00551 !   r-coordinate. interval dr
00552     DO ir = 1, nrbns
00553       drg(ir) = drdr
00554       drginv(ir) = 1.0e0/drg(ir)
00555     END DO
00556 !
00557     rg(0) = rgin
00558     DO ir = 1, nrbns
00559       rg(ir) = rg(ir-1) + drg(ir)
00560 !      write(6,*) ir, rg(ir), drg(ir)
00561     end do
00562 !
00563     rdet  = (rg(nrbns) - rbns)/rbns
00564     IF(dabs(rdet)< 1.e-10) then
00565       WRITE(6,*) 'good coordinate GR : rdet, ratrr ', rdet
00566     end if
00567     IF(dabs(rdet)>=1.e-10)then
00568       WRITE(6,*) ' bad coordinate GR : rdet, ratrr ', rdet
00569       stop
00570     END IF
00571   end if
00572 !
00573 !  Now for r > 3rgmid if mpt1,2
00574 !  or      r > rgmid  if mpt3
00575 !
00576 do itry = 0, 3
00577   if (itry.eq.0) ratrr = 2.0e-01 ! good for finer grid
00578   if (itry.eq.1) ratrr = 2.0e0   ! good for finer grid
00579   if (itry.eq.2) ratrr = 10.0e0  ! good for coarser grid
00580   if (itry.eq.3) ratrr = 50.0e0  ! good for coarsest grid
00581 !
00582   drdr =  1.0d0/dble(nrf)
00583   drdrinv = 1.0d0/drdr
00584   rvdom = rgout - rgmid
00585   nrout = nrg - nrgin
00586   if (rgin==0.0d0) then
00587     drdr =  drg(nrbns)
00588     drdrinv = 1.0d0/drdr
00589     rvdom = rgout - rbns
00590     nrout = nrg - nrbns
00591   end if
00592 !
00593   DO
00594     ratrrb = ratrr
00595     alge = -(ratrr**(nrout+1)-ratrr-rvdom*drdrinv*(ratrr-1.0e0))
00596     dalge = dble(nrout+1)*ratrr**nrout - 1.0e0 - rvdom*drdrinv
00597     ratrr = ratrr + alge/dalge
00598     error = 2.d0*(ratrr - ratrrb)/(ratrr + ratrrb)
00599 !
00600     IF (abs(error)<=1.d-14) EXIT
00601   END DO
00602 
00603   write(6,*) 'region II: ratrr=', ratrr
00604 
00605 ! r-coordinate. interval dr
00606   nr_count = nrgin+1
00607   if (rgin==0.0d0) nr_count = nrbns+1
00608   DO ir = nr_count, nrg
00609     drg(ir) = ratrr*drg(ir-1)
00610     drginv(ir) = 1.0e0/drg(ir)
00611   END DO
00612 ! r-coordinate. grid ponit and mid point.
00613   rg(0) = rgin
00614   if (rg(0).le.1.0d-14) then
00615     rginv(0) = 0.0e0
00616   else
00617     rginv(0) = 1.0d0/rg(0)
00618   end if
00619   DO ir = 1, nrg
00620     rg(ir) = rg(ir-1) + drg(ir)
00621     hrg(ir) = (rg(ir) + rg(ir-1))*0.5e0
00622     rginv(ir) = 1.0e0/rg(ir)
00623     hrginv(ir) = 1.0e0/hrg(ir)
00624 !      write(6,*)ir,rg(ir)
00625   END DO
00626 !
00627 !
00628   rdet  = (rg(nrg) - rgout)/rgout
00629   rdetf = 0.0d0
00630   if (rgin.lt.1.0d0) rdetf =  rg(nrf) - 1.0d0
00631   IF(dabs(rdet)< 1.e-10) then
00632     WRITE(6,*) 'good coordinate itry =', itry
00633     WRITE(6,*) 'good coordinate GR : rdet, ratrr ', rdet, ratrr
00634     WRITE(6,*) 'good coordinate FLUID : rdetf, ratrf', rdetf, ratrf
00635     exit
00636   end if
00637   IF(dabs(rdet)>=1.e-10)then
00638     WRITE(6,*) ' bad coordinate itry =', itry
00639     WRITE(6,*) ' bad coordinate GR : rdet, ratrr', rdet, ratrr
00640     WRITE(6,*) ' bad coordinate FLUID : rdetf, ratrf', rdetf, ratrf
00641     if (itry.eq.3) STOP
00642   END IF
00643 end do
00644 !
00645 end subroutine grid_r_bns_const
00646 !
00647 
00648 
00649 end module coordinate_grav_r
/home/uryu/Cocal_svn_checkinout/Cocal/code/Module/coordinate_grav_r.f90
