subroutine sub_calc_temperature_tillotson
  use mod_variable_ec
  use mod_parameter_ec
  use mod_parameter_tillotson_eos
  use mod_variable
  use mod_parameter_t_e
  use mod_variable_t_e
  implicit none

  integer(4) :: i
  real(8) :: cp,rho0,ec
  real(8) :: ee
  real(8) :: ta_tbl,tb_tbl
  real(8) :: ea_tbl,eb_tbl

  do i=1,np
     if(imat(i) == 1) then
        rho0=any_rho0
     elseif(imat(i) == 2) then
        rho0=murch_rho0
     elseif(imat(i) == 3) then
        rho0=hryugu_rho0
     elseif(imat(i) == 4) then
        rho0=sryugu_rho0
     endif
     
     if(d(i) <= rho0) then
        ec=0.0e0_8
        iec(i) = 0
     else
        do
!           write(*,*) i,iec(i),d(i)
           if(rho_grid(imat(i),iec(i)) <= d(i) .and. &
                &  d(i) < rho_grid(imat(i),iec(i)+1)) then
              exit
           else
              if(d(i) > rho_grid(imat(i),iec(i))) then
                 iec(i)=iec(i)+1
                 if(iec(i) > N_RHO_GRID) then
                    write(*,*) 'ERROR',i,iec(i)
                    stop
                 endif
              else
                 iec(i)=iec(i)-1
                 if(iec(i) < 0) then
                    write(*,*) 'ERROR',i,iec(i)
                 endif
              endif
              cycle
           endif
        enddo
        ec=ec_grid(imat(i),iec(i)) &
             & +(ec_grid(imat(i),iec(i)+1)-ec_grid(imat(i),iec(i))) &
             & /(rho_grid(imat(i),iec(i)+1)-rho_grid(imat(i),iec(i))) &
             & *(d(i)-rho_grid(imat(i),iec(i)))
     endif

     !ec_tmp(i)=ec

     !---temptemp---
     !tt(i)=e(i)/cp
     !if(tt(i) > tt_pk(i)) tt_pk(i)=tt(i)
     !--------------

     ee=e(i)-ec

     do
        if(e_grid(imat(i),it(i)) <= ee .and. ee < e_grid(imat(i),it(i)+1)) then
           exit
        else
           if(ee > e_grid(imat(i),it(i))) then
              it(i)=it(i)+1
              if(it(i) > N_T_GRID_01+N_T_GRID_10) then
                 write(*,*) 'ERROR T_E_Table',i,it(i)
                 stop
              endif
           else
              it(i)=it(i)-1
              if(it(i) < 0) then
                 write(*,*) 'ERROR T_E_Table',i,it(i)
                 stop
              endif
           endif
           cycle
        endif
     enddo

     ta_tbl=t_grid(imat(i),it(i))
     tb_tbl=t_grid(imat(i),it(i)+1)
     ea_tbl=e_grid(imat(i),it(i))
     eb_tbl=e_grid(imat(i),it(i)+1)

     t(i)=ta_tbl+(tb_tbl-ta_tbl)/(eb_tbl-ea_tbl)*(ee-ea_tbl)
     
     !t(i)=(e(i)-ec)/cp

     if(t(i) > t_pk(i)) t_pk(i)=t(i)
     if(p(i) > p_pk(i)) p_pk(i)=p(i)

  enddo
  
  
  return
end subroutine sub_calc_temperature_tillotson