!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine fdm_explicit_init 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, dt, volume_ratio_rock, time, density, density_rock_hydrous
  use global_parameters, only : irmax, dr, temperature_initial, iau, temp
  use material_parameters, only : mn, material_select, ice_water_rock, mass_ratio_water &
    & , mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock, mass_ratio_dehydration_rock
  use energy, only : set_energy, temp_old
  implicit none
  integer :: i
  real(rkind) :: coef, heat_decay, dtr &
    & , kappa, thermometric_conductivity &
    & , k, Tthermal_conductivity, irhoc, Tspecific_heat, Tdensity
  dtr = dt/(dr*dr)
  kappa = thermometric_conductivity(0)
  irhoc = 1._rkind/(Tspecific_heat(0)*Tdensity(0))
  k = Tthermal_conductivity(0)
  temp(0) = temp_old(0) + &
    & 6._rkind*k*irhoc*(temp_old(1)-temp_old(0))*dtr & 
    & + heat_decay(time,0)*dt*irhoc 
  do i = 1, irmax-1
    irhoc = 1._rkind/(Tspecific_heat(i)*Tdensity(i))
    k = Tthermal_conductivity(i)
    coef = k*irhoc*dtr/real(i) 
    temp(i) = temp_old(i) + coef*(real(i+1)*temp_old(i+1) &
      & -2._rkind*real(i)*temp_old(i)+real(i-1)*temp_old(i-1)) &
      & + heat_decay(time,i)*dt*irhoc 
  enddo
  temp(irmax) = temperature_initial(iau)
  return
end subroutine fdm_explicit_init