!--------1---------2---------3---------4---------5---------6---------7---------8
function Tdensity(i)
  use global_constants, only : rkind, density_rock, density_ice, density_water, density_rock_hydrous, density_ice0
  use material_parameters, only : mn, mass_ratio_rock &
    & , mass_ratio_ice, mass_ratio_water &
    & , mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock, mass_ratio_dehydration_rock
  implicit none
  integer,intent(in) :: i
  real(rkind) :: Tdensity
  if(mn(i)==4) then
    Tdensity = density_rock
  elseif(mn(i)==6) then
    Tdensity = 1._rkind/(mass_ratio_rock(i)/density_rock &
      & +mass_ratio_ice(i)/density_ice)
  elseif(mn(i)==7) then
    Tdensity = 1._rkind/(mass_ratio_rock(i)/density_rock &
      & +mass_ratio_ice(i)/density_ice0+mass_ratio_water(i)/density_water)
  elseif(mn(i)==15) then
    Tdensity = 1._rkind/(mass_ratio_anhydrous_rock(i)/density_rock &
      & + mass_ratio_hydrous_rock(i)/density_rock_hydrous &
      & + mass_ratio_dehydration_rock(i)/density_rock &
      & + mass_ratio_water(i)/density_water)
  else
    write(6,*) "Error at thermometric_diffusivity_rock"
  endif
  return
end function Tdensity
!--------1---------2---------3---------4---------5---------6---------7---------8
function thermometric_conductivity(i)
  use global_constants, only : rkind, density, density_rock &
    & , density_ice, density_water, density_rock_hydrous, density_ice0, thermometric_conductivity_ice0
  use material_parameters, only : mn, mass_ratio_rock &
    & , mass_ratio_ice, mass_ratio_water &
    & , mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock, mass_ratio_dehydration_rock
  implicit none
  integer,intent(in) :: i
  real(rkind) :: thermometric_conductivity, thermometric_conductivity_rock &
    & , volume_ratio_rock, volume_ratio_ice, volume_ratio_water &
    & , thermometric_conductivity_ice, thermometric_conductivity_water, Tdensity
  if(mn(i)==4) then
    thermometric_conductivity = thermometric_conductivity_rock(i)
  elseif(mn(i)==6) then
    volume_ratio_rock = Tdensity(i)*mass_ratio_rock(i)/density_rock
    thermometric_conductivity = thermometric_conductivity_rock(i)*volume_ratio_rock &
      & + thermometric_conductivity_ice(i)*(1._rkind-volume_ratio_rock) 
  elseif(mn(i)==7) then
    volume_ratio_rock = Tdensity(i)*mass_ratio_rock(i)/density_rock
    volume_ratio_ice = Tdensity(i)*mass_ratio_ice(i)/density_ice0
    volume_ratio_water = Tdensity(i)*mass_ratio_water(i)/density_water
    thermometric_conductivity = thermometric_conductivity_rock(i)*volume_ratio_rock &
      & + thermometric_conductivity_ice0*volume_ratio_ice &
      & + thermometric_conductivity_water(i)*volume_ratio_water 
  elseif(mn(i)==15) then
    volume_ratio_rock = Tdensity(i)*mass_ratio_rock(i)/density_rock
    thermometric_conductivity = thermometric_conductivity_rock(i)*volume_ratio_rock &
      & + thermometric_conductivity_water(i)*(1._rkind-volume_ratio_rock)
  else
    write(6,*) "Error at thermometric_diffusivity_rock"
  endif
  return
end function thermometric_conductivity
!--------1---------2---------3---------4---------5---------6---------7---------8
function thermometric_conductivity_rock(i)
  use global_constants, only : rkind &
    & , thermometric_conductivity_rockA, thermometric_conductivity_rockB 
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: thermometric_conductivity_rock, temp, itemp
  temp = temp_old(i)
  itemp = 1._rkind/temp
  thermometric_conductivity_rock = thermometric_conductivity_rockA &
   & + thermometric_conductivity_rockB*itemp
  return
end function thermometric_conductivity_rock
!--------1---------2---------3---------4---------5---------6---------7---------8
function thermometric_conductivity_ice(i)
  use global_constants, only : rkind &
    & , thermometric_conductivity_iceA
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: thermometric_conductivity_ice, temp, itemp
  temp = temp_old(i)
  itemp = 1._rkind/temp
  thermometric_conductivity_ice = thermometric_conductivity_iceA 
  return
end function thermometric_conductivity_ice
!--------1---------2---------3---------4---------5---------6---------7---------8
function thermometric_conductivity_water(i)
  use global_constants, only : rkind &
    & , thermometric_conductivity_waterA
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: thermometric_conductivity_water, temp, itemp
  temp = temp_old(i)
  itemp = 1._rkind/temp
  thermometric_conductivity_water = thermometric_conductivity_waterA 
  return
end function thermometric_conductivity_water
!--------1---------2---------3---------4---------5---------6---------7---------8
function Tspecific_heat(i)
  use global_constants, only : rkind, specific_heat_ice0
  use material_parameters, only : mn, mass_ratio_rock &
    & , mass_ratio_ice, mass_ratio_water
  implicit none
  integer,intent(in) :: i
  real(rkind) :: Tspecific_heat, specific_heat_rock &
    & , specific_heat_ice, specific_heat_water
  if(mn(i)==4) then
    Tspecific_heat = specific_heat_rock(i)*mass_ratio_rock(i)
  elseif(mn(i)==6) then
    Tspecific_heat = specific_heat_rock(i)*mass_ratio_rock(i) &
     & + specific_heat_ice(i)*(1._rkind-mass_ratio_rock(i))
  elseif(mn(i)==7) then
    Tspecific_heat = specific_heat_rock(i)*mass_ratio_rock(i) &
     & + specific_heat_ice0*mass_ratio_ice(i) &
     & + specific_heat_water(i)*mass_ratio_water(i)
  elseif(mn(i)==15) then
    Tspecific_heat = specific_heat_rock(i)*mass_ratio_rock(i) &
     & + specific_heat_water(i)*(1-mass_ratio_rock(i))
  else
    write(6,*) "Error at thermal_conductivity_rock"
  endif
  return
end function Tspecific_heat
!--------1---------2---------3---------4---------5---------6---------7---------8
function specific_heat_rock(i)
  use global_constants, only : rkind &
    & , specific_heat_rockA, specific_heat_rockB, specific_heat_rockC
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: specific_heat_rock, temp, itemp
  temp = temp_old(i)
  itemp = 1._rkind/temp_old(i)
  specific_heat_rock = specific_heat_rockA + specific_heat_rockB*temp - specific_heat_rockC*itemp*itemp
  return
end function specific_heat_rock
!--------1---------2---------3---------4---------5---------6---------7---------8
function specific_heat_ice(i)
  use global_constants, only : rkind &
    & , specific_heat_iceA
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: specific_heat_ice, temp, itemp
  temp = temp_old(i)
  itemp = 1._rkind/temp_old(i)
  specific_heat_ice = specific_heat_iceA
  return
end function specific_heat_ice
!--------1---------2---------3---------4---------5---------6---------7---------8
function specific_heat_water(i)
  use global_constants, only : rkind &
    & , specific_heat_waterA
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: specific_heat_water, temp, itemp
  temp = temp_old(i)
  itemp = 1._rkind/temp_old(i)
  specific_heat_water = specific_heat_waterA
  return
end function specific_heat_water
!--------1---------2---------3---------4---------5---------6---------7---------8
function Tthermal_conductivity(i)
  use global_constants, only : rkind, density, density_rock_hydrous&
    & , density_rock, density_ice, density_water, density_rock_hydrous, density_ice0, thermometric_conductivity_ice0
  use material_parameters, only : mn, mass_ratio_rock &
    & , mass_ratio_ice, mass_ratio_water &
    & , mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock, mass_ratio_dehydration_rock
  implicit none
  integer,intent(in) :: i
  real(rkind) :: Tthermal_conductivity, thermal_conductivity_rock &
    & , volume_ratio_rock, volume_ratio_ice, volume_ratio_water &
    & , thermal_conductivity_ice, thermal_conductivity_water, Tdensity
  if(mn(i)==4) then
    Tthermal_conductivity = thermal_conductivity_rock(i)*volume_ratio_rock
  elseif(mn(i)==6) then
    volume_ratio_rock = Tdensity(i)*mass_ratio_rock(i)/density_rock
    Tthermal_conductivity = thermal_conductivity_rock(i)*volume_ratio_rock &
     & + thermal_conductivity_ice(i)*(1._rkind-volume_ratio_rock)
  elseif(mn(i)==7) then
    volume_ratio_rock = Tdensity(i)*mass_ratio_rock(i)/density_rock
    volume_ratio_ice = Tdensity(i)*mass_ratio_ice(i)/density_ice0
    volume_ratio_water = Tdensity(i)*mass_ratio_water(i)/density_water
    Tthermal_conductivity = thermal_conductivity_rock(i)*volume_ratio_rock &
     & + thermometric_conductivity_ice0*volume_ratio_ice &
     & + thermal_conductivity_water(i)*volume_ratio_water
  elseif(mn(i)==15) then
    volume_ratio_rock = Tdensity(i)*mass_ratio_rock(i)/density_rock
    Tthermal_conductivity = thermal_conductivity_rock(i)*volume_ratio_rock &
     & + thermal_conductivity_water(i)*(1._rkind-volume_ratio_rock)
  else
    write(6,*) "Error at thermal_conductivity_rock"
  endif
  return
end function Tthermal_conductivity
!--------1---------2---------3---------4---------5---------6---------7---------8
function thermal_conductivity_rock(i)
  use global_constants, only : rkind, density  
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: thermal_conductivity_rock, specific_heat_rock &
    & , thermometric_conductivity_rock
  thermal_conductivity_rock = thermometric_conductivity_rock(i) & 
   & *specific_heat_rock(i)*density(4)
  return
end function thermal_conductivity_rock
!--------1---------2---------3---------4---------5---------6---------7---------8
function thermal_conductivity_ice(i)
  use global_constants, only : rkind, density  
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: thermal_conductivity_ice, specific_heat_ice &
    & , thermometric_conductivity_ice
  thermal_conductivity_ice = thermometric_conductivity_ice(i) & 
   & *specific_heat_ice(i)*density(2)
  return
end function thermal_conductivity_ice
!--------1---------2---------3---------4---------5---------6---------7---------8
function thermal_conductivity_water(i)
  use global_constants, only : rkind, density 
  use energy, only : temp_old
  implicit none
  integer,intent(in) :: i
  real(rkind) :: thermal_conductivity_water, specific_heat_water &
    & , thermometric_conductivity_water
  thermal_conductivity_water = thermometric_conductivity_water(i) & 
   & *specific_heat_water(i)*density(1)
  return
end function thermal_conductivity_water