!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine fdm_explicit_water_ice
!--------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 &
    & , water_reach_ice, mass_ratio_ice, mass_ratio_water
  use ice_parameters, only : temperature_melting_I
  use energy, only : set_energy, temp_old
  implicit none
  integer :: i, iiw
  real(rkind) :: coef, heat_decay, coef_rad_rock, dtr &
    & , dr1, dr2, ir1
  dtr = dt/(dr*dr)
  if(.NOT.water_reach_ice) then
    radius_water = radius(iw+1)
    mn(iw+1) = 11
    call material_select(iw+1)
    water_reach_ice = .TRUE.
    temp(iw+1) = temperature_melting_I
  endif
  iiw = iw+2
  dr1 = radius(iiw)-radius_water
  dr2 = dr
  ir1 = 1._rkind/radius(iiw)
  call material_select(iiw)
  coef = thermometric_diffusivity(mn(iiw))*dt*2._rkind/(dr1*dr2*(dr1+dr2))
  temp(iiw) = temp_old(iiw) + coef*((1._rkind+dr1*ir1)*dr1*temp_old(iiw+1) &
    & +(dr1+dr2)*(dr2-dr1-radius(iiw))*ir1*temp_old(iiw)+(1._rkind-dr2*ir1)*dr2*temperature_melting_I)&
    & + heat_decay(time,iiw)*dt*thermometric_diffusivity(mn(iiw)) &
    & /thermal_conductivity(mn(iiw))
  do i = iiw+1, 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)
  if(ir/=iw)  temp(ir+2:iw) = temperature_melting_I
  call fdm_explicit_water_rock(temp_old,dtr)
  return
end subroutine fdm_explicit_water_ice