!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine fdm_explicit_water_rock(temp_old,dtr)
!--------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, temp, radius
  use heat_parameters, only : thermal_conductivity, thermometric_diffusivity 
  use material_parameters, only : ir, dr_rock, radius_rock, material_check, iw, mn
  use ice_parameters, only : temperature_melting_I
  implicit none
  integer :: i
  real(rkind),intent(in) :: temp_old(0:irmax), dtr
  real(rkind) :: coef, heat_decay
  call material_check(0)
  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))
  if(ir /= 0) then
    do i = 1, ir-1
      call material_check(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
    call boundary
  else
   write(999,*) 'Stop Checking ir @fdm_expliciti_water_rock'
  endif
  return
end subroutine fdm_explicit_water_rock
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine boundary 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, dt, pi, time, year
  use global_parameters, only : irmax, dr, temp, radius, halfdr
  use heat_parameters, only : thermal_conductivity, thermometric_diffusivity 
  use material_parameters, only : mn, ir, dr_rock, radius_rock, material_check
  use ice_parameters, only : temperature_melting_I
  use energy, only : energy_rock, flux_rock, rad_energy, time_variation
  implicit none
  real(rkind) :: dtdr, dtdr2, idr, heat_decay, fpik, fpirhoc, dtdrb &
    & , deltat, deltar, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, ia532 &
    & , radius_irh, radius3, radius4, radius5 &
    & , radius_irhp, radius3p, radius4p, radius5p, a01, a02
  fpik = 4._rkind*pi*thermal_conductivity(15)
  fpirhoc = fpik/thermometric_diffusivity(15)
  idr = 1._rkind/dr
  radius_irh = radius(ir) - halfdr
  radius3 = (radius_rock**3._rkind - radius_irh**3._rkind)/3._rkind
  radius4 = (radius_rock**4._rkind - radius_irh**4._rkind)*0.25_rkind
  radius5 = (radius_rock**5._rkind - radius_irh**5._rkind)*0.2_rkind
  radius_irhp = radius(ir) + halfdr
  radius3p = (radius_rock**3._rkind - radius_irhp**3._rkind)/3._rkind
  radius4p = (radius_rock**4._rkind - radius_irhp**4._rkind)*0.25_rkind
  radius5p = (radius_rock**5._rkind - radius_irhp**5._rkind)*0.2_rkind
  a1 = fpirhoc*(radius4-radius_rock*radius3)
  a2 = fpirhoc*0.5_rkind*((radius5+(2._rkind*radius_irh-radius_rock) &
    & *radius_rock*radius3) -2._rkind*radius_irh*radius4)
  a3 = radius_rock - radius_irh
  a4 = -dt*fpik*radius_irh**2._rkind
  a5 = dt*fpik*radius_rock**2._rkind
  a6 = energy_rock(1) + heat_decay(time+dt,0)*4._rkind*pi*radius3*dt
  a7 = idr
  a8 = -temp(ir-1)*idr
  a9 = radius(ir) - radius_rock
  a10 = (radius(ir)*radius(ir)-2._rkind*radius_irh*(radius(ir)-radius_rock) &
    & -radius_rock*radius_rock)*0.5_rkind
  a11 = temperature_melting_I
  ia532 = 1._rkind/(a5*a3-a2)
  deltat = a8/a7 + a11 - a10*a6*ia532
  deltar = 1._rkind/a7 - a9 - a10*(a1-a4-a5)*ia532
  dtdr = deltat/deltar
  dtdr2 = ((a1-a4-a5)*dtdr-a6)*ia532
  temp(ir) = a11 + a9*dtdr + a10*dtdr2
  dtdrb = dtdr + a3*dtdr2
  rad_energy(ir) = (a6 - energy_rock(1))/dt
  rad_energy(ir+1) = 0._rkind
  time_variation(ir) = (a4*dtdr + a5*dtdrb)/dt + rad_energy(ir)
  time_variation(ir+1) = 0._rkind
  energy_rock(1) = a4*dtdr + a5*dtdrb + a6
  a01 = fpirhoc*(radius4p-radius_rock*radius3p)
  a02 = fpirhoc*0.5_rkind*(radius5p-2._rkind*radius_irh*radius4p &
    & +(2._rkind*radius_irh-radius_rock)*radius_rock*radius3p)
  energy_rock(0) = a01*dtdr + a02*dtdr2
  flux_rock = -thermal_conductivity(15)*dtdrb
  if(radius_rock > radius(ir+1)) then
    temp(ir+1) = a11 + (radius(ir+1)-radius_rock)*dtdr &
      & + (radius(ir+1)**2._rkind-2._rkind*radius_irh*(radius(ir+1)-radius_rock) &
      & -radius_rock*radius_rock)*0.5_rkind*dtdr2
  elseif(radius_rock <= radius(ir+1)) then
    temp(ir+1) = temperature_melting_I
  endif
end subroutine boundary