!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine fdm_explicit_ice_rock 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, dt, volume_ratio_rock, time, year
  use global_parameters, only : irmax, dr, temperature_initial,iau, temp, radius, halfdr
  use heat_parameters, only : thermal_conductivity, thermometric_diffusivity 
  use material_parameters, only : mn, iw, ir, ii, radius_water, drock, material_select &
    & , radius_rock
  use ice_parameters, only : temperature_melting_I
  use energy, only : set_energy, temp_old
  implicit none
  integer :: i
  real(rkind) :: coef, heat_decay, coef_rad_rock, dtr &
    & , dri, drr, ir1, ir2, ki_drr, kr_dri, temp_B
  dtr = dt/(dr*dr)
  temp(0) = temp_old(0) + &
    & 6._rkind*thermometric_diffusivity(mn(0))*(temp_old(1)-temp_old(0))*dtr &
      & + heat_decay(time,0)*dt*thermometric_diffusivity(mn(0))/thermal_conductivity(mn(0))
  do i = 1, ir
    call material_select(i)
    coef = thermometric_diffusivity(mn(i))*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*thermometric_diffusivity(mn(i))/thermal_conductivity(mn(i))
  enddo

  dri = radius(ir+2) - radius_rock
  drr = radius_rock - radius(ir+1)
  ir1 = 1._rkind/radius(ir+1)
  ir2 = 1._rkind/radius(ir+2)
  ki_drr = thermal_conductivity(2)*drr
  kr_dri = thermal_conductivity(mn(ir))*dri
  temp_B = (ki_drr*temp_old(ir+2) + kr_dri*temp_old(ir+1))/(ki_drr+kr_dri)

  coef = thermometric_diffusivity(mn(ir))*dt*2._rkind/(dr*drr*(dr+drr))
  temp(ir+1) = temp_old(ir+1) + coef*((1._rkind+dr*ir1)*dr*temp_B &
    & +(dr+drr)*(drr-dr-radius(ir+1))*ir1*temp_old(ir+1)+(1._rkind-drr*ir1)*drr*temp_old(ir)) &
    & + heat_decay(time,ir+1)*dt*thermometric_diffusivity(mn(ir)) &
    & /thermal_conductivity(mn(ir))

  coef = thermometric_diffusivity(2)*dt*2._rkind/(dri*dr*(dri+dr))
  temp(ir+2) = temp_old(ir+2) + coef*((1._rkind+dri*ir2)*dri*temp_old(ir+3) &
    & +(dri+dr)*(dr-dri-radius(ir+2))*ir2*temp_old(ir+2)+(1._rkind-dr*ir2)*dr*temp_B)

  do i = ir+3, irmax-1
    call material_select(i)
    coef = thermometric_diffusivity(mn(i))*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*thermometric_diffusivity(mn(i)) &
      & /thermal_conductivity(mn(i))
  enddo
  temp(irmax) = temperature_initial(iau)
  return
end subroutine fdm_explicit_ice_rock