/home/uryu/Cocal_svn_checkinout/Cocal/code/EOS/Subroutine/teos_initialize.f90

Go to the documentation of this file.

subroutine teos_initialize
!
  use phys_constant        !g,c,solmas,nnpeos
  use def_teos_parameter   !abc,abi,rhoi,qi,hi,nphase,rhoini_cgs,emdini_gcm1
  implicit none
!
  real(8) :: rho_0, pre_0, facrho, facpre, fac2, gg, cc, ss
  real(8) :: ene, pre, hh, rho
  integer :: ii, iphase, iofile, idum, i0
  real(8), external      :: fn_lagint
  real(8)                :: x4(4), f4(4)

!
! --  cgs to g = c = msol = 1 unit.
! --  assume pre = pre_0 dyn/cm^2 at rho = rho_0 gr/cm^3.
! --  typically pre_0 = 1.0d+37 dyn/cm^2
! --  and       rho_0 = 1.0d+16  gr/cm^3.
!
  facrho = (g/c**2)**3*solmas**2
  facpre = g**3*solmas**2/c**8

  open(850,file='teos_parameter.dat',status='old')
  nphase=0
  do
   read(850,*,iostat=iofile)  rho_0
   if (iofile > 0)  then
     write(6,*) "Problem reading file teos_parameter.dat...exiting"
     stop
   else if (iofile < 0) THEN
!    nphase is the total number of phases. The total number of points is nphase+1
     nphase = nphase - 1
     write(6,*) "Total number of lines is  =", nphase+1
     write(6,*) "Total number of phases is =", nphase
     exit
   else
     nphase = nphase + 1
   end if
!   write(6,*) nphase, rho_0
  end do
  close(850)

  open(850,file='teos_parameter.dat',status='old')
  do ii=0,nphase
    read(850,*) ene, pre, hh, rho
    enei(ii) = ene*facrho
    rhoi(ii) = rho*facrho
    prei(ii) = pre*facpre
    hi(ii)   = hh
    qi(ii)   = prei(ii)/rhoi(ii) 
!    rhoi(ii) = rhocgs(ii)*facrho
!    prei(ii) = precgs(ii)*facpre
!    write(6,'(i4,4e20.12)') ii, rhoi(ii), prei(ii), hi(ii), qi(ii)
  end do
  close(850)
!
  open(850,file='peos_parameter.dat',status='old')
  read(850,'(8x,1i5,es13.5)') idum, rhoini_cgs
  close(850)
  write(6,'(a21,e20.12)') "******* rhoini_cgs = ", rhoini_cgs
  rhoini_gcm1 = facrho*rhoini_cgs

!  rhoini_gcm1 = 1.73d-03
!  rhoini_cgs  = rhoini_gcm1/facrho

  call teos_lookup(rhoini_gcm1,rhoi,iphase)
  i0 = min0(max0(iphase-2,0),nphase-3)
  x4(1:4) = rhoi(i0:i0+3)
  f4(1:4) = qi(i0:i0+3)
  emdini_gcm1 = fn_lagint(x4,f4,rhoini_gcm1)

!  emdini_gcm1 = 1.132621660533359d-01  !FPS
!  emdini_gcm1 = 1.600980540676673d-01  !gamma=2
  write(6,'(a26,2e23.15)') "(rhoini_gcm1,emdini_gcm1)=", rhoini_gcm1, emdini_gcm1

end subroutine teos_initialize

---